2.4 完整的 Navier–Stokes 方程组
2.4 完整的 Navier–Stokes 方程组
前面几节分别推导了质量、动量和能量守恒定律。现在可将它们组合成一个方程组，以便更清楚地考察其中的各项。为此，回到式（2.2）所表达的矢量物理量的一般守恒定律。出于稍后将说明的原因，引入两个通量矢量，即 $ \vec{F}_{c} $ 和 $ \vec{F}_{v} $。第一个矢量 $ \vec{F}_{c} $ 与流体中各物理量的对流输运有关，通常称为对流通量矢量；不过，对于动量方程和能量方程，它还分别包含压力项 $ \vec{p}\vec{n} $（式（2.5））和 $ p(\vec{v}\cdot\vec{n}) $（式（2.11））。第二个矢量称为黏性通量矢量 $ \vec{F}_{v} $，包含
黏性应力和热扩散。此外，定义源项 $ \vec{Q} $，其中包括由体力和体积加热引起的全部体积源项。在此基础上，将相应的张量或矢量与单位法向量 $ \vec{n} $ 作内积，便可结合式（2.3）、（2.5）和（2.13），将式（2.2）整理为
$$ \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 v\\ \rho w\\ \rho E\end{bmatrix}.\tag{2.20} $$
对流通量矢量为
$$ \vec{F}_{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\end{bmatrix}\tag{2.21} $$
其中，逆变速度 V，即垂直于面元 dS 的速度，定义为速度矢量与单位法向量的内积：
$$ V\equiv\vec{v}\cdot\vec{n}=n_{x}u+n_{y} v+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_{y}\tau_{xy}+n_{z}\tau_{xz}\\n_{x}\tau_{y x}+n_{y}\tau_{yy}+n_{z}\tau_{y z}\\n_{x}\tau_{zx}+n_{y}\tau_{zy}+n_{z}\tau_{zz}\\n_{x}\Theta_{x}+n_{y}\Theta_{y}+n_{z}\Theta_{z}\end{bmatrix},\tag{2.23} $$
其中，
$$ \begin{aligned}\Theta_{x}&=u\tau_{xx}+ v\tau_{xy}+w\tau_{xz}+k\frac{\partial T}{\partial x},\\\Theta_{y}&=u\tau_{y x}+ v\tau_{yy}+w\tau_{y z}+k\frac{\partial T}{\partial y},\\\Theta_{z}&=u\tau_{zx}+ v\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,y}\\ \rho f_{e,z}\\ \rho\vec{f}_{e}\cdot\vec{v}+\dot{q}_{h}\end{bmatrix}.\tag{2.25} $$
对于 Newton 流体，即黏性应力满足关系式（2.15）的流体，上述方程组 $ ((2.19)-(2.25)) $ 称为 Navier–Stokes 方程。它们描述通过边界 $ \partial\Omega $ 发生的质量、动量和能量交换，即通量；该边界围成一个固定在空间中的控制体 $ \Omega $（见图 2.1）。我们根据守恒定律，推导出了 Navier–Stokes 方程的积分形式。应用 Gauss 定理，可将式（2.19）改写为微分形式 [7]。由于文献中经常使用微分形式，为完整起见，第 A.1 节也列出了这一形式。
在某些情况下，例如透平机械应用或地球物理问题中，控制体会绕某一轴旋转，通常为匀速旋转。此时，应将 Navier–Stokes 方程变换到旋转参考系中。因此，必须在源项 $ \vec{Q} $ 中加入 Coriolis 力和离心力的影响 [8]。相应的 Navier–Stokes 方程形式见第 A.4 节。另一些情况下，控制体可能发生平移或变形，例如研究流固相互作用时。此时，Navier–Stokes 方程（2.19）中必须增加一个描述面元 dS 相对于固定坐标系运动的项 [9]。此外，还必须满足所谓的几何守恒律 [10–12]。相应表述见第 A.5 节。
2.4.1 完全气体的表述
在纯空气动力学问题中，通常可以合理地假设工作流体为量热完全气体，其状态方程具有以下形式 $ [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}+ v^{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 的单位为 Kelvin（K）。因此，在 $ T = 288 \, K $ 时可得 $ \mu = 1.78 \cdot 10^{-5} \, kg/ms $。对于气体，导热系数 k 的温度依赖性与 $ \mu $ 类似；对于液体，k 则基本保持不变。正因如此，空气通常采用以下关系：
$$ k=c_{p}\frac{\mu}{Pr}\tag{2.31} $$
此外，一般假设 Prandtl 数 Pr 在整个流场中为常数。对于空气，Prandtl 数取 Pr = 0.72。
2.4.2 真实气体的表述
处理真实气体时，问题会更加复杂，因为除了流体动力学，还需要建立热力学过程和化学反应的模型。真实气体流动的实例包括燃烧模拟、再入飞行器周围的高超声速流动，以及汽轮机内的流动。
原则上，可采用两种不同的方法求解这一问题。第一种方法适用于气体同时处于化学平衡和热力学平衡的情况，此时存在唯一的状态方程。控制方程（2.19）保持不变，只需采用曲线拟合，从查询表中插值得到压力、温度、黏度等数值 $ [15-18] $。但在实际问题中，气体更常处于化学非平衡和／或热力学非平衡状态，必须据此建模。
为便于说明，考虑一种由 N 个不同组分组成的气体混合物。Damköhler 数定义为流动驻留时间与化学反应时间之比。当该数为有限值时，模型中必须考虑有限速率化学反应，以描述化学反应引起的组分生成与消耗。下文还将假设，流体动力学和化学反应的时间尺度与空间尺度，远大于热力学过程的相应尺度。因此，气体在热力学上处于平衡状态，在化学上则处于非平衡状态。为模拟这种气体混合物的行为，需要在 Navier–Stokes 方程中增加描述 N 个组分的 $ (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 v\\\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}+ v\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_{y}=u\tau_{y x}+ v\tau_{yy}+w\tau_{y z}+k\frac{\partial T}{\partial y}+\rho\sum_{m=1}^{N}h_{m} D_{m}\frac{\partial Y_{m}}{\partial y},\\&\Theta_{z}=u\tau_{zx}+ v\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_{y,m}=\rho D_{m}\frac{\partial Y_{m}}{\partial y},\\&\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,y}\\ \rho f_{e,z}\\ \rho\vec{f}_{e}\cdot\vec{v}+\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 l=1,2,\cdots,N_{R}.\tag{2.41} $$
其中，$ v_{lm}^{\prime} $ 和 $ v_{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),\\&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}_{l m}=(\nu_{l m}^{\prime\prime}-\nu_{l m}^{\prime})\left(K_{\mathrm{f}l}\prod_{n=1}^{N}C_{n}^{\nu_{l n}^{\prime}}-K_{\mathrm{b}l}\prod_{n=1}^{N}C_{n}^{\nu_{l n}^{\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} $$
更多细节可参见前述参考文献。文献 $ ^{[28]} $ 也对描述化学反应流的方程进行了详细综述，同时给出了通量的 Jacobian 矩阵及其特征值。
真实气体的另一个实际应用，是透平机械中的蒸汽流动模拟，或要求更高的湿蒸汽流动模拟 $ [29-39] $。对于后者，蒸汽中混有水滴，因此属于多相流。模拟时，既可以求解一组附加输运方程，也可以沿若干条流线跟踪水滴运动。这些模拟在现代汽轮机叶栅设计中具有十分重要的应用。例如，分析叶片周围的流动，有助于理解凝结引起的超临界激波和流动不稳定现象；这些现象会给叶片带来附加动载荷，并导致效率下降。
2.4.3 Navier–Stokes 方程的简化
薄剪切层近似
对于高 Reynolds 数的绕体流动，也就是边界层相对于特征尺寸很薄的情况，可对 Navier–Stokes 方程（2.19）进行简化。一个必要条件是不存在大范围的边界层分离。在这种情况下，可以认为，只有流动物理量沿物体表面法向的梯度（图 2.4 中的 $ \eta $ 方向）对黏性应力有贡献 [40, 41]。计算剪应力张量，即式（2.14）和（2.15）时，则忽略其他坐标方向的梯度（图 2.4 中的 $ \xi $ 方向）。这就是所谓的 Navier–Stokes 方程薄剪切层（thin shear layer，TSL）近似。采用 TSL 简化的目的，是降低黏性项数值计算的开销，同时在这些假设成立的条件下，使解仍保持足够的精度。从实际计算的角度，也可以解释 TSL 近似的合理性。对于高 Reynolds 数流动，
为充分解析边界层，壁面法向的网格必须非常细密。受计算机内存与运算速度限制，其他方向只能采用粗得多的网格，因此，其梯度计算的数值精度明显低于法向的精度。为完整起见，第 A.6 节列出了 TSL 方程。由于这一简化不能充分解析二次流，例如叶栅中的二次流，因此 TSL 近似通常仅用于外部空气动力学。
抛物化 Navier–Stokes 方程
当以下三个条件同时成立时：
-
流动为定常流动，即 $ \partial\vec{W}/\partial t=0 $；
-
流体主要沿一个主方向运动，例如不能存在边界层分离；
-
横向流动分量可以忽略；
控制方程（2.19）便可简化为所谓的抛物化 Navier–Stokes（parabolized Navier–Stokes，PNS）方程 $ [8, 42–44] $。上述条件允许在黏性应力项，即式（2.15）中，将 u、v、w 对流向的导数设为零。此外，还可从式（2.23）的黏性通量矢量中，删去沿流向的黏性应力张量 $ \overline{\overline{\tau}} $ 分量、黏性应力做功 $ (\overline{\overline{\tau}} \cdot \vec{v}) $ 分量，以及导热 $ k \nabla T $ 分量。连续性方程和对流通量，即式（2.21），保持不变。详情见第 A.7 节。考虑图 2.5 所示的情形，主流方向与 x 坐标方向一致。可以证明，PNS 近似得到的是一组抛物型与椭圆型混合方程。其中，流向动量方程和能量方程为抛物型，因此可沿 x 方向推进求解；y 方向和 z 方向的动量方程为椭圆型，必须在各个 x 截面内迭代求解。因此，PNS 方法的主要优点，是大幅降低流动求解的复杂程度：将完整的三维流场求解，转化为
Euler 方程
$$ \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 方程用于流动模拟的频率已相对较低。