2.4 完整Navier-Stokes方程组

前几节中，我们分别推导了质量、动量与能量守恒定律。现在我们将它们汇总成一个方程组，以更清晰地总览其中涉及的各项。为此，回到式(2.2)所示的一般向量守恒定律。出于稍后要解释的原因，我们引入两个通量向量$\vec{F}_{c}$与$\vec{F}_{v}$。第一个$\vec{F}_{c}$与流体中各量的对流输运相关，通常称为对流通量向量，尽管对动量与能量方程而言，它还包含压力项$p\vec{n}$（式(2.5)）与$p(\vec{\nu}\cdot\vec{n})$（式(2.11)）。第二个通量向量$\vec{F}_{v}$称为粘性通量向量，包含粘性应力与热扩散。此外，定义源项$\vec{Q}$，它包含由体力与体积加热产生的所有体源。考虑与单位法向量$\vec{n}$的标量积后，可将式(2.2)与式(2.3)、(2.5)、(2.13)写为  

$$\frac{\partial}{\partial t}\int_{\Omega}\vec{W}\mathrm{d}\Omega+\oint_{\partial\Omega}(\vec{F}_{\mathrm{c}}-\vec{F}_{\mathrm{v}})\mathrm{d}S=\int_{\Omega}\vec{Q}\mathrm{d}\Omega.\tag{2.19}$$  

所谓守恒变量向量$\vec{W}$在三维中由五个分量组成  

$$\vec{W}=\begin{bmatrix}\rho\\ \rho u\\ \rho\nu\\ \rho w\\ \rho E\end{bmatrix}.\tag{2.20}$$  

对流通量向量为  

$$\vec{F}_{c}=\begin{bmatrix}\rho V\\ \rho u V+n_{x}p\\ \rho\nu V+n_{\gamma}p\\ \rho w V+n_{z}p\\ \rho H V\end{bmatrix}\tag{2.21}$$  

其中逆变速度$V$——即垂直于面元$\mathrm{d}S$的速度——定义为速度向量与单位法向量的标量积  

$$V\equiv\vec{\nu}\cdot\vec{n}=n_{x}u+n_{\gamma}\nu+n_{z}w.\tag{2.22}$$  

式(2.21)中的总焓$H$由式(2.12)给出。粘性通量向量由式(2.14)得  

$$\vec{F}_{\mathrm{v}}=\begin{bmatrix}0\\n_{x}\tau_{xx}+n_{\gamma}\tau_{x\gamma}+n_{z}\tau_{xz}\\n_{x}\tau_{\gamma x}+n_{\gamma}\tau_{\gamma\gamma}+n_{z}\tau_{\gamma z}\\n_{x}\tau_{zx}+n_{\gamma}\tau_{zy}+n_{z}\tau_{zz}\\n_{x}\Theta_{x}+n_{\gamma}\Theta_{\gamma}+n_{z}\Theta_{z}\end{bmatrix},\tag{2.23}$$  

其中  

$$\begin{aligned}\Theta_{x}&=u\tau_{xx}+\nu\tau_{xy}+w\tau_{xz}+k\frac{\partial T}{\partial x},\ \Theta_{\gamma}=u\tau_{\gamma x}+\nu\tau_{\gamma\gamma}+w\tau_{\gamma z}+k\frac{\partial T}{\partial\gamma},\ \Theta_{z}=u\tau_{zx}+\nu\tau_{zy}+w\tau_{zz}+k\frac{\partial T}{\partial z}\end{aligned}\tag{2.24}$$  

分别描述粘性应力与热传导在流体中所做的功项。最后，源项为  

$$\vec{Q}=\begin{bmatrix}0\\ \rho f_{e,x}\\ \rho f_{e,\gamma}\\ \rho f_{e,z}\\ \rho\vec{f}_{e}\cdot\vec{\nu}+\dot{q}_{h}\end{bmatrix}.\tag{2.25}$$  

对于牛顿流体（即粘性应力满足式(2.15)），上述方程组$((2.19)-(2.25))$称为Navier-Stokes方程。它们描述质量、动量与能量通过固定在空间中的控制体$\Omega$边界$\partial\Omega$的交换（通量）（见图2.1）。我们依据守恒定律推导了Navier-Stokes方程的积分形式。利用Gauss定理，式(2.19)可改写为微分形式$[7]$。由于文献中常见微分形式，为完整起见，我们在A.1节给出。

 

在某些情况下，例如涡轮机械应用或地球物理中，控制体绕某轴旋转（通常为匀速旋转）。此时Navier-Stokes方程需变换到旋转参考系中，因此源项$\vec{Q}$必须加入Coriolis力与离心力的影响$[8]$。相应形式见A.4节。在其他情况下，控制体可能发生平移或变形，例如研究流固耦合时。此时，Navier-Stokes方程(2.19)需要增加一项，用以描述面元$\mathrm{d}S$相对于固定坐标系的相对运动$[9]$。此外，还必须满足所谓的几何守恒定律$[10-12]$。我们在A.5节给出相应表述。

 

Navier-Stokes方程在三维中是由五个方程组成的系统，对应五个守恒变量$\rho$、$\rho u$、$\rho v$、$\rho w$与$\rho E$。但它们包含七个未知流场变量：$\rho$、u、v、w、E、p与T。因此需要提供两条附加方程，用以给出状态变量之间的热力学关系。例如，将压力表示为密度与温度的函数，并将内能或焓表示为压力与温度的函数。除此之外，还必须给出动力粘度系数$\mu$与导热系数$k$作为流体状态函数的关系，以完成整个方程组。显然，这些关系取决于所考虑的流体类型。下面我们将给出两种常见情形下的闭合方法。

在纯气动问题中，通常可合理假设工质为热量学意义上的理想气体（calorically perfect gas），其状态方程为$[13,14]$  

$$p=\rho R T,\tag{2.26}$$  

其中$R$为比气体常数。焓由下式给出  

$$h=c_{p}T.\tag{2.27}$$  

将压力用守恒变量表示会更方便。为此，将式(2.12)（总焓与总能量关系）与状态方程(2.26)结合。代入焓表达式(2.27)，并利用  

$$R=c_{p}-c_{\mathrm{v}},\quad\gamma=\frac{c_{p}}{c_{\mathrm{v}}},\tag{2.28}$$  

最终得到压力  

$$p=(\gamma-1)\rho\biggl[E-\frac{u^{2}+\nu^{2}+w^{2}}{2}\biggr].\tag{2.29}$$  

温度随后借助关系式(2.26)计算。对于理想气体，动力粘度系数$\mu$对温度依赖很强，而对压力依赖很弱。常用Sutherland公式。对空气（SI单位）有  

$$\mu=\frac{1.45T^{3/2}}{T+110}\cdot10^{-6},\tag{2.30}$$  

其中温度$T$以开尔文（K）计。因此在$T=288\,K$时，得到$\mu=1.78\cdot10^{-5}\,kg/ms$。导热系数$k$在气体情形下对温度的依赖与$\mu$相似；而在液体情形下$k$几乎为常数。因此，对空气一般使用  

$$k=c_{p}\frac{\mu}{Pr}\tag{2.31}$$  

并常假设Prandtl数$Pr$在全流场中为常数。对空气，$Pr=0.72$。

当需要处理真实气体时，问题会更复杂。原因在于除了流体力学之外，还必须对热力学过程与化学反应进行建模。真实气体流动的例子包括：燃烧模拟、再入飞行器的高超声速绕流，或蒸汽涡轮内流动。

 

原则上可采用两种不同方法。第一种方法适用于气体在化学与热力学上处于平衡的情况。这意味着存在唯一的状态方程，此时控制方程(2.19)保持不变，只是压力、温度、粘度等的数值通过曲线拟合从查表数据插值得到$[15-18]$。但在实践中，气体更常处于化学与/或热力学非平衡状态，因此必须作相应建模。

 

为说明起见，考虑由$N$种不同组分构成的气体混合物。对有限Damköhler数（定义为流动驻留时间与化学反应时间之比）而言，需要在模型中加入有限速率化学反应，以描述由化学反应导致的组分生成/消耗。下面进一步假设：流体力学与化学反应的时间尺度和空间尺度相对热力学过程要大得多。因此我们认为气体在热力学上处于平衡，但在化学上处于非平衡。为模拟这种混合气体行为，必须在Navier-Stokes方程上增加$(N-1)$条组分输运方程$[19-27]$。于是形式上得到与式(2.19)相同的系统，但守恒变量向量$\vec{W}$、通量向量$\vec{F}_{c}$与$\vec{F}_{v}$以及源项$\vec{Q}$都要扩展，加入$(N-1)$条组分方程。回顾式(2.20)至(2.25)，此时守恒变量向量为  

$$\vec{W}=\left[\begin{array}{c}\rho\\\rho u\\\rho\nu\\\rho w\\\rho E\\\rho Y_{1}\\\vdots\\\rho Y_{N-1}\end{array}\right].\tag{2.32}$$  

对流与粘性通量向量变为  

$$\vec{F}_{\mathrm{c}}=\begin{bmatrix}\rho V\\ \rho u V+n_{x}p\\ \rho v V+n_{y}p\\ \rho w V+n_{z}p\\ \rho H V\\ \rho Y_{1}V\\ \vdots\\ \rho Y_{N-1}V\end{bmatrix},\quad\vec{F}_{\mathrm{v}}=\begin{bmatrix}0\\ n_{x}\tau_{x x}+n_{y}\tau_{x y}+n_{z}\tau_{x z}\\ n_{x}\tau_{y x}+n_{y}\tau_{y y}+n_{z}\tau_{y z}\\ n_{x}\tau_{z x}+n_{y}\tau_{z y}+n_{z}\tau_{z z}\\ n_{x}\Theta_{x}+n_{y}\Theta_{y}+n_{z}\Theta_{z}\\ n_{x}\Phi_{x,1}+n_{y}\Phi_{y,1}+n_{z}\Phi_{z,1}\\ \vdots\\ n_{x}\Phi_{x,N-1}+n_{y}\Phi_{y,N-1}+n_{z}\Phi_{z,N-1}\end{bmatrix},\tag{2.33}$$  

其中  

$$\begin{aligned}&\Theta_{x}=u\tau_{xx}+\nu\tau_{xy}+w\tau_{xz}+k\frac{\partial T}{\partial x}+\rho\sum_{m=1}^{N}h_{m}D_{m}\frac{\partial Y_{m}}{\partial x},\ \Theta_{\gamma}=u\tau_{\gamma x}+\nu\tau_{\gamma\gamma}+w\tau_{\gamma z}+k\frac{\partial T}{\partial\gamma}+\rho\sum_{m=1}^{N}h_{m}D_{m}\frac{\partial Y_{m}}{\partial\gamma},\ \Theta_{z}=u\tau_{zx}+\nu\tau_{zy}+w\tau_{zz}+k\frac{\partial T}{\partial z}+\rho\sum_{m=1}^{N}h_{m}D_{m}\frac{\partial Y_{m}}{\partial z},\ \Phi_{x,m}=\rho D_{m}\frac{\partial Y_{m}}{\partial x},\ \Phi_{\gamma,m}=\rho D_{m}\frac{\partial Y_{m}}{\partial\gamma},\ \Phi_{z,m}=\rho D_{m}\frac{\partial Y_{m}}{\partial z}.\end{aligned}\tag{2.34}$$  

源项变为  

$$\vec{Q}=\begin{bmatrix}0\\ \rho f_{e,x}\\ \rho f_{e,\gamma}\\ \rho f_{e,z}\\ \rho\vec{f}_{e}\cdot\vec{\nu}+\dot{q}_{h}\\ \dot{s}_{1}\\ \vdots\\ \dot{s}_{N-1}\end{bmatrix}.\tag{2.35}$$  

在上述表达式(2.32)–(2.35)中，$Y_{m}$为质量分数，$h_{m}$为焓，$D_{m}$为第$m$种组分的有效二元扩散系数。此外，$\dot{s}_{m}$为由化学反应导致的第$m$种组分变化率。注意混合物总密度$\rho$等于各组分密度$\rho Y_{m}$之和。因此，由于总密度被视为独立量，仅剩$(N-1)$个独立密度$\rho Y_{m}$；剩余的质量分数$Y_{N}$由  

$$Y_{N}=1-\sum_{m=1}^{N-1}Y_{m}.\tag{2.36}$$  

得到。

为得到压力$p$的表达式，先假设各组分都表现为理想气体，即  

$$p_{m}=\rho Y_{m}\frac{R_{u}}{W_{m}}T,\tag{2.37}$$  

其中$R_{u}$为通用气体常数，$W_{m}$为分子量。结合Dalton定律  

$$p=\sum_{m=1}^{N}p_{m},\tag{2.38}$$  

可写为  

$$p=\rho R_{u}T\sum_{m=1}^{N}\frac{Y_{m}}{W_{m}}.\tag{2.39}$$  

需要强调的是：由于气体在热力学上处于平衡，所有组分具有相同温度$T$。温度需由下式迭代求解$[22,28]$  

$$e=\sum_{m=1}^{N}\left[Y_{m}\left(h_{f,m}^{0}+\int_{T_{\mathrm{ref}}}^{T}c_{p,m}\mathrm{d}T\right)\right]-\frac{p}{\rho}.\tag{2.40}$$  

混合气体内能$e$由式(2.6)得到。式(2.40)中的$h_{f,m}^{0}$、$c_{p,m}$与$T_{ref}$分别表示第$m$种组分的生成焓、定压比热与参考温度。上述量以及组分的导热系数$k$与动力粘度$\mu$的取值由曲线拟合确定$[20,22,24]$。

最后尚需建模的是式(2.35)中的化学源项$\dot{s}_{m}$。对于包含$N$种组分、$N_{R}$个基元反应的反应集合，其速率方程可写为一般形式  

$$\sum_{m=1}^{N}\nu_{lm}^{\prime}C_{m}\underset{K_{bl}}{\overset{K_{fl}}{\rightleftharpoons}}\sum_{m=1}^{N}\nu_{lm}^{\prime\prime}C_{m}\qquad for\quad l=1,2,\cdots,N_{R}.\tag{2.41}$$  

其中$\nu_{lm}^{\prime}$与$\nu_{lm}^{\prime\prime}$分别为第$l$个正/逆反应中第$m$种组分的化学计量系数；$C_{m}$为第$m$种组分的摩尔浓度（$C_{m}=\rho Y_{m}/W_{m}$）；$K_{fl}$与$K_{bl}$分别为第$l$步反应的正/逆反应速率常数。它们由经验Arrhenius公式给出  

$$\begin{aligned}&K_{\mathrm{f}}=A_{\mathrm{f}}T^{B_{\mathrm{f}}}\exp(-E_{\mathrm{f}}/R_{u}T),\quad K_{\mathrm{b}}=A_{\mathrm{b}}T^{B_{\mathrm{b}}}\exp(-E_{\mathrm{b}}/R_{u}T),\end{aligned}$$  

其中$A_{f}$与$A_{b}$为Arrhenius系数，$E_{f}$与$E_{b}$为活化能，$B_{f}$与$B_{b}$为常数。第$l$个反应对第$m$种组分摩尔浓度的变化率为  

$$\dot{C}_{lm}=(\nu_{lm}^{\prime\prime}-\nu_{lm}^{\prime})\left(K_{\mathrm{f}l}\prod_{n=1}^{N}C_{n}^{\nu_{ln}^{\prime}}-K_{\mathrm{b}l}\prod_{n=1}^{N}C_{n}^{\nu_{ln}^{\prime\prime}}\right).\tag{2.42}$$  

因此结合式(2.42)，第$m$种组分的总变化率为  

$$\dot{s}_{m}=W_{m}\sum_{l=1}^{N_{R}}\dot{C}_{lm}.\tag{2.43}$$  

更多细节可参见上述引用文献。关于化学反应流动控制方程的详细概览，以及通量的Jacobian矩阵与特征值，也可参见文献$^{[28]}$。

 

真实气体的另一个实际例子是蒸汽的模拟，或在涡轮机械应用中更具挑战性的湿蒸汽模拟$[29-39]$。在后者情形下，蒸汽与水滴混合，属于多相流：既可以求解一组额外的输运方程，也可以沿若干流线追踪水滴。这类模拟在现代蒸汽涡轮叶栅设计中非常重要。例如，对叶片绕流的分析可帮助理解由凝结引发的超临界激波以及流动不稳定性——它们会导致叶片额外的动态载荷并造成效率损失。

下面我们考虑Navier-Stokes方程(2.19)的三种常见简化形式。此处我们将注意力限制在每种近似背后的物理理由上。前两种简化形式的方程在附录中给出。

在模拟高雷诺数下的绕体流动时（即边界层相对于特征尺度很薄），Navier-Stokes方程(2.19)可以简化。一个必要条件是：不存在大范围的边界层分离。此时可预期，只有沿物体表面法向方向的流动量梯度（图2.4中的$\eta$方向）对粘性应力有贡献$[40,41]$。另一方面，在计算剪应力张量（式(2.14)与式(2.15)）时，可忽略其他坐标方向（图2.4中的$\xi$）上的梯度。这称为Navier-Stokes方程的薄剪切层（TSL）近似。TSL修正的动机在于：粘性项的数值计算将更省计算量，同时在假设成立的范围内仍能保持足够精度。从实践角度也可为TSL近似辩护：在高雷诺数流动中，为正确解析边界层，网格在壁面法向必须非常细；受限于计算机内存与速度，其他方向不得不采用更粗网格，这会导致相对于法向方向而言梯度计算的数值精度显著降低。为完整起见，TSL方程在A.6节给出。由于二次流（例如叶栅中的二次流）无法得到恰当分辨，TSL简化通常只用于外流气动问题。

当满足以下三个条件时：

 

控制方程(2.19)可简化为一种称为抛物化Navier-Stokes（PNS）方程的形式$[8,42-44]$。上述条件允许我们在粘性应力项（式(2.15)）中令$u$、$v$与$w$对主流向的导数为零。此外，还将主流向上的粘性应力张量$\overline{\tau}$分量、由其所做功$(\overline{\tau}\cdot\vec{v})$以及热传导$k\nabla T$从式(2.23)的粘性通量向量中略去。连续方程以及对流通量（式(2.21)）保持不变。细节见A.7节。对图2.5所示主流方向与$x$坐标一致的情形，可证明PNS近似得到一组抛物/椭圆混合型方程：主流向的动量方程与能量方程为抛物型，因此可沿$x$方向推进求解；而$y$与$z$方向的动量方程为椭圆型，必须在每个$x$截面内迭代求解。因此，PNS方法的主要收益在于显著降低了求解复杂性——从完整三维场问题变为一系列二维问题。PNS方程的典型应用包括管道与通道内流动的计算，以及利用空间推进方法模拟定常超声速流动$[45-48]$。

如前所述，Navier-Stokes方程描述粘性流体的行为。在许多情况下，完全忽略粘性效应是合理近似，例如高雷诺数流动中边界层相对于物体尺度非常薄。此时可直接从式(2.19)中去掉粘性通量向量$\vec{F}_{v}$，得到  

$$\frac{\partial}{\partial t}\int_{\Omega}\vec{W}\mathrm{d}\Omega+\oint_{\partial\Omega}\vec{F}_{\mathrm{c}}\mathrm{d}S=\int_{\Omega}\vec{Q}\mathrm{d}\Omega.\tag{2.44}$$  

其余各项仍由同样的关系式(2.20)-(2.22)与式(2.25)给出。这种简化形式称为Euler方程。它们描述无粘流体中流动量的纯对流。若以守恒形式表述（如上式），Euler方程能够准确表示激波、膨胀波以及三角翼（锐前缘）上的涡等重要现象。此外，Euler方程过去并且现在仍是发展离散方法与边界条件的重要基础。

 

不过需要指出的是：由于如今即使个人计算机也具备相当的计算能力，同时对模拟质量的要求提高，Euler方程用于流动模拟的情况相对已经不多。