第 2 章：流体力学控制方程

流动模型

控制体 空间固定 随流体运动
有限 积分守恒 积分非守恒
微元 微分守恒 微分非守恒

这里所谓“守恒型”和“非守恒型”是指同一个物理守恒定律采用了不同的数学表达方式。

二者在流场足够光滑时可以相互转换。

自变量与因变量

在控制方程中，自变量是空间坐标

(x,y,z)(x,y,z)

与时间

tt

流场中的物理量则是这些自变量的函数，例如

ρ=ρ(x,y,z,t)\rho=\rho(x,y,z,t) u=u(x,y,z,t)u=u(x,y,z,t) v=v(x,y,z,t)v=v(x,y,z,t) w=w(x,y,z,t)w=w(x,y,z,t) p=p(x,y,z,t)p=p(x,y,z,t) T=T(x,y,z,t)T=T(x,y,z,t)

因此压强 pp 同样是因变量，而不是预先给定的常数。

这些因变量并不全部相互独立。例如对于理想气体，

p=ρRTp=\rho RT

因此 ρ,p,T\rho,p,T 之间存在状态方程约束。

对于三维可压缩流，通常只需要五个独立的流场未知量，对应：

  • 一个质量方程；
  • 三个动量方程；
  • 一个能量方程。

如果采用守恒变量，常写作

U=[ρρuρvρwρE]\boxed{ \mathbf U= \begin{bmatrix} \rho\\ \rho u\\ \rho v\\ \rho w\\ \rho E \end{bmatrix} }

这也是CFD所偏爱的解向量

需要注意：

“数值计算中直接求解的变量”和“流场中的因变量”并不是完全相同的概念。

例如压强 pp 是因变量，但在可压缩流守恒格式中，通常不直接推进 pp，而是由 ρ,ρE,ρu,ρv,ρw\rho,\rho E,\rho u,\rho v,\rho w 和状态方程恢复得到。

物质导数

对于某个流场变量

ϕ=ϕ(x,y,z,t)\phi=\phi(x,y,z,t)

如果观察的是一个随流体运动的流体微团，那么它的位置本身满足

x=x(t),y=y(t),z=z(t)x=x(t),\quad y=y(t),\quad z=z(t)

因此沿该流体质点观察到的 ϕ\phi 的变化率为

DϕDt=∂ϕ∂t+u∂ϕ∂x+v∂ϕ∂y+w∂ϕ∂z\frac{D\phi}{Dt} = \frac{\partial\phi}{\partial t} + u\frac{\partial\phi}{\partial x} + v\frac{\partial\phi}{\partial y} + w\frac{\partial\phi}{\partial z}

即

DDt=∂∂t+V⋅∇\boxed{ \frac{D}{Dt} = \frac{\partial}{\partial t} + \mathbf V\cdot\nabla }

其中

∂ϕ∂t\frac{\partial\phi}{\partial t}

表示固定空间位置处由于流场本身随时间变化造成的变化，称为局部变化；

而

V⋅∇ϕ\mathbf V\cdot\nabla\phi

表示流体质点运动到空间中具有不同 ϕ\phi 值的位置而产生的变化，称为迁移变化（对流变化）。

因此物质导数并不是一种新的微积分运算，本质上来自多元函数的链式法则。

但更准确地说：

物质导数是沿流体质点在时空中的运动轨迹所取的方向导数。

守恒型方程的数学特征

判断一个方程是不是守恒型，最直接的办法不是看它是否“来自守恒定律”，而是看它能否写成

∂U∂t+∇⋅F=S\boxed{ \frac{\partial U}{\partial t} + \nabla\cdot\mathbf F = S }

其中：

  • UU 是单位体积内的某个守恒量；
  • F\mathbf F 是该守恒量的通量；
  • SS 是源项。

记方程中守恒量的单位为 xx

U : x⋅m−3x \cdot m^{-3}

F : x⋅m−2⋅s−1x \cdot m^{-2} \cdot s^{-1}

S : x⋅m−3⋅s−1x \cdot m^{-3} \cdot s^{-1} (守恒量在 单位体积/固定空间微元 的变化速率)

其物理结构是

控制体内的积累+通过边界的净流出=内部产生\text{控制体内的积累} + \text{通过边界的净流出} = \text{内部产生}

守恒型最重要的数学特征是：空间输运项以一个整体通量的散度形式出现。

例如

∇⋅(ρuV)\nabla\cdot(\rho u\mathbf V)

是守恒型的动量输运项，而

ρV⋅∇u\rho\mathbf V\cdot\nabla u

则不是直接的守恒形式。

非守恒型到守恒型到底转换了什么

以 xx 方向动量方程中的惯性项为例：

ρDuDt\rho\frac{Du}{Dt}

它描述的是：

跟随某个流体微团时，该微团的 xx 方向动量变化。

将物质导数展开：

ρDuDt=ρ∂u∂t+ρV⋅∇u\rho\frac{Du}{Dt} = \rho\frac{\partial u}{\partial t} + \rho\mathbf V\cdot\nabla u

但 CFD 中更希望得到

∂(ρu)∂t+∇⋅(ρuV)\frac{\partial(\rho u)}{\partial t} + \nabla\cdot(\rho u\mathbf V)

因为：

ρu\rho u

才是单位体积的 xx 方向动量。

因此整个转换的目标其实从一开始就已经确定：

把“速度的物质变化”改写为“动量在固定控制体中的积累 + 动量穿过控制面的通量”。

第一步使用乘积法则，把

ρ∂u∂t\rho\frac{\partial u}{\partial t}

改写成含有

∂(ρu)∂t\frac{\partial(\rho u)}{\partial t}

的形式。

目的不是单纯的代数变形，而是为了构造出真正的动量密度的时间变化率。

第二步将

ρV⋅∇u\rho\mathbf V\cdot\nabla u

改写成含有

∇⋅(ρuV)\nabla\cdot(\rho u\mathbf V)

的形式。

目的则是构造动量通量的散度。

这两个变换都会额外产生一些项，而这些额外项最后恰好组合成

u[∂ρ∂t+∇⋅(ρV)]u \left[ \frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\mathbf V) \right]

括号中的部分正是连续方程，因此等于零。

于是得到

ρDuDt=∂(ρu)∂t+∇⋅(ρuV)\boxed{ \rho\frac{Du}{Dt} = \frac{\partial(\rho u)}{\partial t} + \nabla\cdot(\rho u\mathbf V) }

所以：

非守恒型到守恒型的转换，本质上是利用质量守恒，把“跟随质点观察的变化”重新解释成“固定控制体中的积累与流出”。

为什么 CFD 特别重视守恒形式

对于一个网格单元，如果方程写成

∂U∂t+∇⋅F=S\frac{\partial U}{\partial t} + \nabla\cdot\mathbf F = S

那么对网格单元积分并使用 Gauss 定理后，空间项就会变成控制面上的通量。

相邻两个网格单元共享同一个面，因此一个单元流出的量恰好是另一个单元流入的量。

于是内部通量天然相互抵消。

这意味着数值离散之后仍然可以非常自然地保持

总流出=总流入\text{总流出}=\text{总流入}

这也是有限体积法能够天然保持离散守恒的根本原因。

对于光滑流动，守恒型和非守恒型通常数学等价。

但出现激波等间断后，经典导数可能不存在，此时非守恒形式中的乘积和导数已经不能随意使用普通微积分规则。

而守恒形式仍然可以通过积分意义定义弱解。

因此对于可压缩流和激波问题：

守恒形式并不仅仅是“另一种写法”，而是正确描述间断解的重要基础。

CFD中5个控制方程更常见写法就是：

∂U∂t+∂F∂x+∂G∂y+∂H∂z=0.\frac{\partial \mathbf U}{\partial t} + \frac{\partial \mathbf F}{\partial x} + \frac{\partial \mathbf G}{\partial y} + \frac{\partial \mathbf H}{\partial z} =0.

这些 U\mathbf U、F\mathbf F、G\mathbf G、H\mathbf H 均以列向量的形式给出。

动量方程为什么能够生成动能方程

动量方程本质上是

质量×加速度=力\text{质量}\times\text{加速度} = \text{力}

将三个方向的动量方程分别乘以对应的速度 u,v,wu,v,w 后，左侧出现

uDuDt,vDvDt,wDwDtu\frac{Du}{Dt}, \quad v\frac{Dv}{Dt}, \quad w\frac{Dw}{Dt}

利用链式法则：

uDuDt=DDt(u22)u\frac{Du}{Dt} = \frac{D}{Dt}\left(\frac{u^2}{2}\right)

三个方向相加后就得到

DDt(u2+v2+w22)\frac{D}{Dt} \left( \frac{u^2+v^2+w^2}{2} \right)

即单位质量动能的变化率。

因此：

动量方程乘速度，本质上就是把“力学中的牛顿第二定律”转化成“功—能关系”。

这也是为什么动量方程能够用来构造动能方程。

总能量、动能和内能之间的消去

单位质量总能量可以写成

E=e+V22E=e+\frac{V^2}{2}

其中：

  • ee 是内能；
  • V2/2V^2/2 是宏观平动动能。

因此

总能量方程−动能方程=内能方程\boxed{ \text{总能量方程} - \text{动能方程} = \text{内能方程} }

书中所谓利用动量方程进行“消去”，实际上不是把某些物理作用删除，而是：

先用动量方程识别出总能量变化中属于宏观动能的那部分，再从总能量账本中扣除。

因此最后剩下的就是改变内能的机制。

例如体积力做功

ρf⋅V\rho\mathbf f\cdot\mathbf V

同时存在于总能量方程和动能方程中，因此相减后消失。

这不意味着重力“不做功”，而意味着：

体积力做的这部分功已经归入宏观动能，不直接属于内能变化。

压力项则不会完全消失。

相减后留下的

−p∇⋅V-p\nabla\cdot\mathbf V

表示压缩和膨胀引起的内能变化。

流体被压缩时，

∇⋅V<0\nabla\cdot\mathbf V<0

因此

−p∇⋅V>0-p\nabla\cdot\mathbf V>0

内能增加。

黏性作用同样不会全部消失，最后会留下黏性耗散项。

它表示

有序的宏观机械能→无序的分子热运动\text{有序的宏观机械能} \rightarrow \text{无序的分子热运动}

因此属于内能增加的一部分。

像这样衍生出来的这些单独内能、单独静焓、单独总焓方程，也是同样的，并非独立的方程，而是利用其它方程在能量的总体上做消去。

黏性应力为什么与速度梯度有关

速度场

V=(u,v,w)\mathbf V=(u,v,w)

的速度梯度并不是一个单独的量，而是一个张量：

∇V\nabla\mathbf V

其中包含所有形式的

∂ui∂xj\frac{\partial u_i}{\partial x_j}

速度散度

∇⋅V\nabla\cdot\mathbf V

只是速度梯度张量的迹：

∇⋅V=∂u∂x+∂v∂y+∂w∂z\nabla\cdot\mathbf V = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z}

因此，正应力中出现速度散度并不意味着它“不依赖速度梯度”；恰恰相反，速度散度本身就是速度梯度的一部分。

对于牛顿流体，黏性应力真正与速度梯度的对称部分有关：

D=12[∇V+(∇V)T]\mathbf D = \frac12 \left[ \nabla\mathbf V + (\nabla\mathbf V)^T \right]

其中 D\mathbf D 描述流体微团的形变速率。

速度梯度中的反对称部分主要对应刚体旋转，它只改变流体微团的方向，而不改变形状，因此不会产生通常意义上的黏性形变应力。

所以可以把牛顿流体的基本逻辑理解为：

速度梯度→流体微团形变速率→黏性应力\boxed{ \text{速度梯度} \rightarrow \text{流体微团形变速率} \rightarrow \text{黏性应力} }

其中：

  • 速度散度主要描述体积膨胀和压缩；
  • 交叉速度梯度主要描述剪切形变。