4.3 对流通量的离散

4.3 对流通量的离散

前面各节讨论了空间离散的一般方法。本节将详细介绍如何近似对流通量。

如第 3.1.5 节所述，在有限体积方法中，基本可以选择：

• 中心格式；

• 通量矢量分裂格式；

• 通量差分分裂格式；

• 总变差减小（TVD）格式；

• 波动分裂格式。

为控制篇幅，本章仅介绍最重要且应用最广的方法，不详细描述基本格式的全部可能修改，而是给出相应参考文献。

在详细介绍各种离散格式之前，先说明左状态、右状态以及模板（也称计算分子）的含义。

某些单元中心与对偶控制体格式要求将流动变量插值到控制体面。图 4.8 以 i 方向网格为例说明。中心格式（见下一节）采用的一种方法，是用面左右两侧相同数量的值进行线性插值，即插值以该面为中心。基于 Euler 方程特征性质的离散——迎风格式——则用非对称公式，分别从面左侧和右侧插值流动变量。所得两个值称为左状态和右状态，再用它们计算穿过该面的对流通量（见式（4.18）、（4.23）、（4.38）或（4.43））。除 TVD 格式外，插值公式几乎都基于 Van Leer 的 MUSCL（守恒律的单调上游中心格式）方法 [30]。对于一般流动变量 U，公式为：

$$ \begin{aligned}&U_{\mathrm{R}}=U_{I+1}-\frac{\epsilon}{4}\left[(1+\hat{\kappa})\Delta_{-}+(1-\hat{\kappa})\Delta_{+}\right]U_{I+1},\\&U_{\mathrm{L}}=U_{I}\quad+\frac{\epsilon}{4}\left[(1+\hat{\kappa})\Delta_{+}+(1-\hat{\kappa})\Delta_{-}\right]U_{I}.\\ \end{aligned}\tag{4.46} $$

前向差分算子 $ (\Delta_{+}) $ 和后向差分算子 $ (\Delta_{-}) $ 定义为：

$$ \begin{aligned}&\Delta_{+}U_{I}=U_{I+1}-U_{I},\\&\Delta_{-}U_{I}=U_{I}-U_{I-1}.\\ \end{aligned}\tag{4.47} $$

下标按需要作相应平移。对于采用对偶控制体的单元顶点格式，只要将 I 替换为节点下标 i，上述关系仍成立。将参数 $ \epsilon $ 设为零，即可得到一阶迎风离散。参数 $ \hat{\kappa} $ 决定插值的空间精度。当 $ \epsilon = 1 $ 且 $ \hat{\kappa} = -1 $ 时，式（4.46）给出完全单边的流动变量插值，在等间距网格上得到二阶迎风近似。$ \hat{\kappa} = 0 $ 对应二阶偏迎风线性插值。取 $ \hat{\kappa} = 1/3 $，则得到三点插值公式，在有限体积框架内（参见 [31, 49]）构成二阶偏迎风格式，其截断误差小于 $ \hat{\kappa} = -1 $ 和 $ \hat{\kappa} = 0 $ 格式。最后，取 $ \hat{\kappa} = 1 $ 时，MUSCL 方法退化为纯中心格式，即变量平均。实践中最常采用 $ \hat{\kappa} = 0 $ 和 $ \hat{\kappa} = 1/3 $ 的格式。

如果流动区域内存在强梯度，必须为 MUSCL 插值（4.46）加入限制器函数（简称限制器），以抑制解的非物理振荡。第 4.3.5 节将进一步讨论限制器。

模板或计算分子，指参与残差、梯度等计算的单元形心或网格点的集合。例如，若分别按式（4.15）、（4.20）或（4.35）、（4.40）平均控制体面通量，则二维中得到由以下单元／点构成的五点模板：

$$ \begin{array}{llll}(I,J)&(I+1,J)&(I,J+1)&(I-1,J)&(I,J-1).\end{array} $$

三维中则得到包含以下单元／点的七点模板：

$$ \begin{aligned}&(I,J,K)\quad(I+1,J,K)\quad(I,J+1,K)\quad(I-1,J,K)\\&(I,J-1,K)\quad(I,J,K-1)\quad(I,J,K+1).\end{aligned} $$

图 4.9 展示了两种模板。注意，在笛卡尔网格上，它们对应于以下方向一阶导数的二阶中心差分近似：

$ i-, j- $ 以及 k 方向。因此，对控制方程微分形式采用中心差分的有限差分格式，会给出相同结果。

4.3.1 带人工耗散的中心格式

与其他离散方法相比，带人工耗散的中心格式非常简单。它容易在单元中心格式或两种单元顶点格式中实现，因此得到广泛应用。

中心格式的基本思想，是对控制体面两侧的守恒变量作算术平均，计算该面上的对流通量。由于这会允许解的奇偶解耦（离散方程产生两个相互独立的解）以及激波处的过冲，必须加入类似黏性通量的人工耗散以保持稳定。Jameson 等人 $ [7] $ 首次将该格式用于 Euler 方程；根据作者姓氏的首字母，也称 JST 格式。

中心格式分辨间断和边界层的精度通常不如迎风格式，但计算成本低得多。因此，研究者尝试在维持较低计算量的同时提高精度，例如降低边界层内的人工耗散 $ [32, 33] $，或提高激波分辨能力 $ [34–36] $。文献 $ [36] $ 的另一思想，是利用对流通量的 Jacobian 矩阵，分别缩放各守恒方程的耗散。这种成功的方法称为矩阵耗散格式，可视为原始标量格式与迎风格式之间的折中。基本格式还使用单一的压力传感器，在间断处由二阶精度切换到一阶精度，以防止流动变量的非物理振荡。Jameson 在文献 $ [37] $ 中提出了 SLIP（对称受限保正）格式的概念，其中采用限制器。

对每个守恒方程分别施加限制器。文献 [38] 简要介绍了上述方法，并比较了二维无黏与黏性流的结果。

标量耗散格式

穿过控制体面的对流通量（式（2.21）），分别按式（4.16）、（4.21）、（4.36）或（4.41）用变量平均近似。随后，为保持稳定，在中心通量上加入人工耗散 [7, 39]。因此，面 $ (I + 1/2, J, K) $ 上的总对流通量为：

$$ \begin{array}{r}{\big(\vec{F}_{c}\Delta S\big)_{I+1/2,J,K}\approx\vec{F}_{c}(\vec{W}_{I+1/2,J,K})\Delta S_{I+1/2,J,K}-\vec{D}_{I+1/2,J,K},}\end{array}\tag{4.48} $$

流动变量按下式平均（另见图 4.8）：

$$ \vec{W}_{I+1/2,J,K}=\frac{1}{2}\left(\vec{W}_{I,J,K}+\vec{W}_{I+1,J,K}\right).\tag{4.49} $$

对于采用对偶控制体的单元顶点格式，改用节点下标 $ (i,j,k) $。为简单起见，下文将 $ (I+1/2,J,K) $ 简记为 $ (I+1/2) $。人工耗散通量由自适应二阶和四阶差分混合组成，它们分别通过对一阶和三阶差分算子求和得到。

$$ \begin{aligned}\vec{D}_{I+1/2}=&\ \hat{\Lambda}_{I+1/2}^{S}\left[\epsilon_{I+1/2}^{(2)}(\vec{W}_{I+1}-\vec{W}_{I})\right.\\&\left.-\epsilon_{I+1/2}^{(4)}(\vec{W}_{I+2}-3\vec{W}_{I+1}+3\vec{W}_{I}-\vec{W}_{I-1})\right].\end{aligned}\tag{4.50} $$

由式（4.50）可见，该格式在二维中具有紧致的九点模板，在三维中具有十三点模板。

耗散用各坐标方向对流通量 Jacobian 谱半径之和进行尺度调整：

$$ \hat{\Lambda}_{I+1/2}^{S}=(\hat{\Lambda}_{\mathrm{c}}^{I})_{I+1/2}+(\hat{\Lambda}_{\mathrm{c}}^{J})_{I+1/2}+(\hat{\Lambda}_{\mathrm{c}}^{K})_{I+1/2}.\tag{4.51} $$

例如，在 I 方向（以上标 I 表示），单元面 $ (I + 1/2) $ 处的谱半径由平均获得：

$$ (\hat{\Lambda}_{c}^{I})_{I+1/2}=\frac{1}{2}\left[(\hat{\Lambda}_{c}^{I})_{I}+(\hat{\Lambda}_{c}^{I})_{I+1}\right].\tag{4.52} $$

其计算公式为：

$$ \hat{\Lambda}_{\mathrm{c}}=\left(\left|V\right|+c\right)\Delta S,\tag{4.53} $$

其中，V 为逆变速度（2.22），c 为声速。

采用基于压力的传感器，在激波处关闭会引起强烈解振荡的四阶差分。该传感器还会在流场光滑区域关闭二阶差分，将耗散降低到尽可能小的水平。因此，式（4.50）的系数 $ \epsilon^{(2)} $ 与 $ \epsilon^{(4)} $ 定义为：

$$ \begin{array}{l}\epsilon_{I+1/2}^{(2)}=k^{(2)}\max(\Upsilon_{I},\Upsilon_{I+1}),\\\epsilon_{I+1/2}^{(4)}=\max\left[0,(k^{(4)}-\epsilon_{I+1/2}^{(2)})\right]\end{array}\tag{4.54} $$

压力传感器为：

$$ \Upsilon_{I}=\frac{\left|p_{I+1}-2p_{I}+p_{I-1}\right|}{p_{I+1}+2p_{I}+p_{I-1}}.\tag{4.55} $$

参数的典型取值是 $ k^{(2)} = 1/2 $ 和 $ 1/128 \leq k^{(4)} \leq 1/64 $。

为减少跨黏性剪切层的人工耗散，可以重新定义式（4.50）的尺度因子 $ [32, 33] $：

$$ \begin{aligned}\hat{\Lambda}_{I+1/2}^{S}&=\frac{1}{2}\left[(\phi^{I}\hat{\Lambda}_{c}^{I})_{I}+(\phi^{I}\hat{\Lambda}_{c}^{I})_{I+1}\right],\\\hat{\Lambda}_{J+1/2}^{S}&=\frac{1}{2}\left[(\phi^{J}\hat{\Lambda}_{c}^{J})_{J}+(\phi^{J}\hat{\Lambda}_{c}^{J})_{J+1}\right],\\\hat{\Lambda}_{K+1/2}^{S}&=\frac{1}{2}\left[(\phi^{K}\hat{\Lambda}_{c}^{K})_{K}+(\phi^{K}\hat{\Lambda}_{c}^{K})_{K+1}\right].\end{aligned}\tag{4.56} $$

用它们代替式（4.51）。与方向有关的系数 $ \phi $ 由下列关系给出：

$$ \begin{array}{l}{\phi^{I}=1+\max\left[\left(\frac{\hat{\Lambda}_{c}^{J}}{\hat{\Lambda}_{c}^{I}}\right)^{\sigma},\left(\frac{\hat{\Lambda}_{c}^{K}}{\hat{\Lambda}_{c}^{I}}\right)^{\sigma}\right],}\\ {\phi^{J}=1+\max\left[\left(\frac{\hat{\Lambda}_{c}^{I}}{\hat{\Lambda}_{c}^{J}}\right)^{\sigma},\left(\frac{\hat{\Lambda}_{c}^{K}}{\hat{\Lambda}_{c}^{J}}\right)^{\sigma}\right],}\\ {\phi^{K}=1+\max\left[\left(\frac{\hat{\Lambda}_{c}^{I}}{\hat{\Lambda}_{c}^{K}}\right)^{\sigma},\left(\frac{\hat{\Lambda}_{c}^{J}}{\hat{\Lambda}_{c}^{K}}\right)^{\sigma}\right].}\end{array}\tag{4.57} $$

参数 $ \sigma $ 通常取 1/2 或 2/3。对于长边沿流动方向的控制体，这种表述减小了沿短边方向的耗散项尺度。

矩阵耗散格式

为通过降低数值耗散提高精度，可修改上述 JST 格式，使其更接近迎风格式。其思想是，以矩阵——对流通量 Jacobian——取代标量 $ \hat{\Lambda}^{S} $，用于缩放耗散项 [36]。这样，每个方程均由相应特征值恰当地缩放。因此，式（4.50）变为：

$$ \begin{array}{l}{\vec{D}_{I+1/2}=\left|\bar{A}_{\mathrm{c}}\right|_{I+1/2}\left[\epsilon_{I+1/2}^{(2)}\left(\vec{W}_{I+1}-\vec{W}_{I}\right)\right.}\\ {\quad\left.-\epsilon_{I+1/2}^{(4)}\left(\vec{W}_{I+2}-3\vec{W}_{I+1}+3\vec{W}_{I}-\vec{W}_{I-1}\right)\right].}\end{array}\tag{4.58} $$

尺度矩阵对应于对流通量 Jacobian $ (\bar{A}_{c} = \partial \vec{F}_{c} / \partial \vec{W}) $，其对角化时采用特征值的绝对值：

$$ |\bar{A}_{c}|=\bar{T}|\bar{\Lambda}_{c}\Delta S|\bar{T}^{-1}.\tag{4.59} $$

右特征向量矩阵 $ (\bar{T}) $、左特征向量矩阵 $ (\bar{T}^{-1}) $ 及特征值对角矩阵 $ \bar{\Lambda}_{c} $ 见第 A.11 节。必须在驻点和声速线处对特征值作限制，以防耗散变为零。文献 [36] 给出了高效计算 $ |\bar{A}_{c}| $ 与 $ \vec{W} $ 乘积的方法。注意，取 $ \epsilon^{(2)} = 1/2 $、$ \epsilon^{(4)} = 0 $，便得到一阶精度的完全迎风格式。

如前所述，这里的目标是构造一种精度接近迎风格式、但计算成本仅比标量耗散方法略高（约 15%–20%）的格式。文献 [38] 报道了它与通量矢量分裂格式（CUSP 和 AUSM）的比较结果。

4.3.2 通量矢量分裂格式

通量矢量分裂方法可视为第一层次的迎风格式，因为它只考虑波的传播方向。该类格式将对流通量矢量分成两部分：或者按某些特征变量的符号分裂，或者分为对流部分与压力部分。著名的 Van Leer 通量矢量分裂格式 $ [40] $ 属于第一类，基于特征分解。较新的方法采用第二种思路，例如 Liou 等人的对流迎风分裂法（AUSM）$ [41, 42] $，以及 Jameson 的对流迎风分裂压力（CUSP）格式 $ [43, 44] $。类似方法还包括 Edwards 的低扩散通量分裂格式（LDFSS）$ [45] $ 和 Rossow 的基于 Mach 数的对流—压力分裂（MAPS）格式 $ [46, 47] $。

通量矢量分裂格式只能在单元中心格式（第 4.2.1 节）或对偶控制体单元顶点格式（第 4.2.3 节）中实现。其优势是，相较于采用标量人工耗散的中心格式，计算量仅适度增加，而激波和边界层的分辨能力明显改善。不过，矩阵耗散格式（式（4.58））也能得到精度相当的结果 $ [38] $。由于某些数值困难，许多研究者对基本格式，特别是 AUSM，提出了多种修改，相关发展仍在继续。下面介绍 Van Leer、AUSM 和 CUSP 格式的基础，指出最重要的修改，并提供相应文献。

Van Leer 格式

Van Leer 通量矢量分裂格式 $ [40] $ 基于对流通量的特征分解。该方法在贴体网格上的扩展见文献 $ [48, 49] $。

对流通量分为正、负两部分，即：

$$ \vec{F}_{\mathrm{c}}=\vec{F}_{\mathrm{c}}^{+}+\vec{F}_{\mathrm{c}}^{-},\tag{4.60} $$

分裂依据是控制体面法向 Mach 数（例如 $ (I + 1/2) $ 处，见图 4.8）：

$$ (M_{n})_{I+1/2}=\left(\frac{V}{c}\right)_{I+1/2},\tag{4.61} $$

其中，V 为逆变速度（2.22），c 为声速。对于对偶控制体单元顶点格式，需要将单元下标替换为节点下标。

首先，按式（4.19）或（4.39），将流动变量 $ \rho $、u、v、w 和 p 插值到控制体面。随后，用左状态计算正通量，用右状态计算负通量。对流 Mach 数 $ (M_{n})_{I+1/2} $ 由下式获得 [40]：

$$ (M_{n})_{I+1/2}=M_{\mathrm{L}}^{+}+M_{\mathrm{R}}^{-},\tag{4.62} $$

分裂 Mach 数定义为：

$$ M_{\mathrm{L}}^{+}=\left\{\begin{array}{l l}M_{\mathrm{L}}&\text{当}M_{\mathrm{L}}\geq+1,\\ \frac{1}{4}(M_{\mathrm{L}}+1)^{2}&\text{当}|M_{\mathrm{L}}|<1,\\ 0&\text{当}M_{\mathrm{L}}\leq-1,\end{array}\right.\tag{4.63} $$

以及：

$$ M_{\mathrm{R}}^{-}=\left\{\begin{array}{l l}0&\text{当}M_{\mathrm{R}}\geq+1,\\ \frac{1}{4}(M_{\mathrm{R}}-1)^{2}&\text{当}|M_{\mathrm{R}}|<1,\\ M_{\mathrm{R}}&\text{当}M_{\mathrm{R}}\leq-1.\end{array}\right.\tag{4.64} $$

校注：原书式（4.64）的亚声速分支写为正号，译文保留原式；它与负向 Mach 数分裂的符号约定存在疑点，亦请对照式（4.67）的负质量通量。

Mach 数 $ M_{L} $ 和 $ M_{R} $ 分别根据左、右状态计算，即：

$$ M_{\mathrm{L}}=\frac{V_{\mathrm{L}}}{c_{\mathrm{L}}},\quad M_{\mathrm{R}}=\frac{V_{\mathrm{R}}}{c_{\mathrm{R}}}.\tag{4.65} $$

当 $ |M_{n}| < 1 $（亚声速流）时，正、负通量部分为：

$$ \vec{F}_{\mathrm{c}}^{\pm}=\left[\begin{aligned}f_{\mathrm{mass}}^{\pm}\\ f_{\mathrm{mass}}^{\pm}\left[n_{x}(-V\pm2c)/\gamma+u\right]\\ f_{\mathrm{mass}}^{\pm}\left[n_{y}(-V\pm2c)/\gamma+ v\right]\\ f_{\mathrm{mass}}^{\pm}\left[n_{z}(-V\pm2c)/\gamma+w\right]\\ f_{\mathrm{energy}}^{\pm}\end{aligned}\right].\tag{4.66} $$

质量通量与能量通量分量定义为：

$$ \begin{aligned}&f_{\mathrm{mass}}^{+}=+\rho_{\mathrm{L}}c_{\mathrm{L}}\frac{(M_{\mathrm{L}}+1)^{2}}{4},\\&f_{\mathrm{mass}}^{-}=-\rho_{\mathrm{R}}c_{\mathrm{R}}\frac{(M_{\mathrm{R}}-1)^{2}}{4},\\&f_{\mathrm{energy}}^{\pm}=f_{\mathrm{mass}}^{\pm}\left\{\frac{[(\gamma-1)V\pm2c]^{2}}{2(\gamma^{2}-1)}+\frac{u^{2}+ v^{2}+w^{2}-V^{2}}{2}\right\}_{\mathrm{L/R}}.\\ \end{aligned}\tag{4.67} $$

对于超声速流，即 $ |M_{n}| \geq 1 $，通量由下式计算：

$$ \begin{aligned}\vec{F}_{c}^{+}&=\vec{F}_{c},\quad\vec{F}_{c}^{-}=0&&\text{当 }M_n\geq+1,\\\vec{F}_{c}^{+}&=0,\quad\vec{F}_{c}^{-}=\vec{F}_{c}&&\text{当 }M_n\leq-1.\end{aligned}\tag{4.68} $$

左、右状态通常按式（4.46）的 MUSCL 方法 $ [30] $ 计算。当流场包含激波等间断时，高阶格式 $ \hat{\kappa} = -1 $、$ \hat{\kappa} = 0 $ 和 $ \hat{\kappa} = 1/3 $ 需要采用限制器。详见第 4.3.5 节。

Van Leer 通量矢量分裂格式用于 Euler 方程时表现很好。然而，多项 Navier–Stokes 方程研究 $ [50, 51] $ 表明，动量和能量方程的分裂误差会抹平边界层，并导致驻点温度及壁面温度不准确。因此，文献 $ [52] $ 建议修改边界层法向的动量通量；文献 $ [53] $ 对能量通量提出了类似修正。两种修改结合，可以消除分裂误差，从而显著改善解的精度 $ [54] $。

AUSM

对流迎风分裂法（AUSM）由 Liou 和 Steffen $ [41] $ 以及 Liou $ [55] $ 提出。随后，Wada 和 Liou $ [56] $ 对其修改，并更名为 AUSMD/V。最后，Liou [42, 57] 提出了称为 AUSM$ ^{+} $ 的改进版本。

该方法基于这样的观察：对流通量矢量（2.21）由物理上不同的两部分组成，即对流部分和压力部分：

$$ \vec{F}_{c}=V\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho w\\ \rho H\end{bmatrix}+\begin{bmatrix}0\\ n_{x}p\\ n_{y}p\\ n_{z}p\\ 0\end{bmatrix}.\tag{4.69} $$

式（4.69）的第一项表示由逆变速度 V 对流输运的标量；压力项则受声波传播速度控制。由此，对流项采用纯迎风离散，按 V 的符号选择左状态或右状态，即使亚声速流也如此。另一方面，亚声速情况下压力项同时包含两个状态，只有在超声速流中才变为完全迎风。

按照 $ [41] $ 的基本 AUSM 表述，引入式（4.61）的对流 Mach 数 $ (M_{n})_{I+1/2} $。于是，可将控制体面 $ (I+1/2) $ 或 $ (i+1/2) $ 上的对流通量改写为：

$$ (\vec{F}_{c})_{I+1/2}=(M_{n})_{I+1/2}\begin{bmatrix}\rho c\\ \rho c u\\ \rho c v\\ \rho c w\\ \rho c H\end{bmatrix}_{\mathrm{L/R}}+\begin{bmatrix}0\\ n_{x}p\\ n_{y}p\\ n_{z}p\\ 0\end{bmatrix}_{I+1/2},\tag{4.70} $$

其中，

$$ (\bullet)_{\mathrm{L}/\mathrm{R}}=\left\{\begin{array}{l l}(\bullet)_{\mathrm{L}}&\text{当}M_{I+1/2}\geq0,\\ (\bullet)_{\mathrm{R}}&\text{否则}\end{array}\right.\tag{4.71} $$

与 Van Leer 通量矢量分裂格式类似，对流 Mach 数根据式（4.62）及（4.63）至（4.65），计算为左右分裂 Mach 数之和。左、右状态（流动物理量 $ \rho $、u、v、w、p、H）同样通过按式（4.19）或（4.39）分别向控制体面插值得到。插值采用 MUSCL 方法 [30]，公式见式（4.46）。当流场包含激波等强梯度时，全部高阶 MUSCL 格式（$ \hat{\kappa} = -1 $、$ \hat{\kappa} = 0 $ 和 $ \hat{\kappa} = 1/3 $）都需要限制器。详见第 4.3.5 节。

控制体面 $ (I+1/2) $ 上的压力由分裂获得 [41]：

$$ p_{I+1/2}=p_{\mathrm{L}}^{+}+p_{\mathrm{R}}^{-},\tag{4.72} $$

分裂压力在文献 $ [40] $ 中给出：

$$ p_{\mathrm{L}}^{+}=\left\{\begin{array}{l l}p_{\mathrm{L}}&\text{当}M_{\mathrm{L}}\geq+1,\\ \frac{p_{\mathrm{L}}}{4}(M_{\mathrm{L}}+1)^{2}(2-M_{\mathrm{L}})&\text{当}|M_{\mathrm{L}}|<1,\\ 0&\text{当}M_{\mathrm{L}}\leq-1,\end{array}\right.\tag{4.73} $$

以及：

$$ p_{\mathrm{R}}^{-}=\left\{\begin{array}{l l}0&\text{当}M_{\mathrm{R}}\geq+1,\\ \frac{p_{\mathrm{R}}}{4}(M_{\mathrm{R}}-1)^{2}(2+M_{\mathrm{R}})&\text{当}|M_{\mathrm{R}}|<1,\\ p_{\mathrm{R}}&\text{当}M_{\mathrm{R}}\leq-1.\end{array}\right.\tag{4.74} $$

对于 $ \left|M_{L/R}\right| < 1 $，也可以采用下列较低阶展开：

$$ p_{\mathrm{L/R}}^{\pm}=\frac{p_{\mathrm{L/R}}}{2}\left(1\pm M_{\mathrm{L/R}}\right).\tag{4.75} $$

应当注意，AUSM 还可以写成：

$$ \begin{aligned}(\vec{F}_{c})_{I+1/2}=&\frac{1}{2}(M_{n})_{I+1/2}\left\{\left[\begin{array}{c}\rho c\\\rho cu\\\rho cv\\\rho cw\\\rho c H\end{array}\right]_{\mathrm{L}}+\left[\begin{array}{c}\rho c\\\rho cu\\\rho cv\\\rho cw\\\rho c H\end{array}\right]_{\mathrm{R}}\right\}\\&-\frac{1}{2}\left|(M_{n})_{I+1/2}\right|\left\{\left[\begin{array}{c}\rho c\\\rho cu\\\rho cv\\\rho cw\\\rho c H\end{array}\right]_{\mathrm{R}}-\left[\begin{array}{c}\rho c\\\rho cu\\\rho cv\\\rho cw\\\rho c H\end{array}\right]_{\mathrm{L}}\right\}\\&+\left[\begin{array}{c}0\\n_{x}(p_{\mathrm{L}}^{+}+p_{\mathrm{R}}^{-})\\n_{y}(p_{\mathrm{L}}^{+}+p_{\mathrm{R}}^{-})\\n_{z}(p_{\mathrm{L}}^{+}+p_{\mathrm{R}}^{-})\\0\end{array}\right].\end{aligned}\tag{4.76} $$

式（4.76）右侧第一项为左、右状态的 Mach 数加权平均，类似于通量平均式（4.15）、（4.20）或（4.35）、（4.40）。第二项具有耗散性质，用标量 $ \left|(M_{n})_{I+1/2}\right| $ 缩放。

实践证明，AUSM 能清晰分辨强激波，并准确计算边界层。然而，原始 AUSM $ [41, 55] $ 在激波处及流动与网格对齐时会产生局部压力振荡 $ [58] $。

因此，文献 [58, 59] 建议在激波处切换到 Van Leer 格式。当对流 Mach 数 $ (M_{n})_{I+1/2} $ 趋于零时，式（4.76）的耗散项也趋于零，格式无法衰减扰动。为解决流动与网格对齐的问题，文献 [58] 建议将耗散项的尺度修改为：

$$ \left|(M_{n})_{I+1/2}\right|=\left\{\begin{array}{l l}{\left|(M_{n})_{I+1/2}\right|}&{\text{当}\left|(M_{n})_{I+1/2}\right|>\delta,}\\ {\displaystyle\frac{(M_{n})_{I+1/2}^{2}+\delta^{2}}{2\delta}}&{\text{当}\left|(M_{n})_{I+1/2}\right|\leq\delta,}\end{array}\right.\tag{4.77} $$

其中，$ \delta $ 是一个小量（$ 0 < \delta \leq 0.5 $），从而始终保留足够的数值耗散。为维持 AUSM 对边界层的精度，可沿壁面法向减小参数 $ \delta $，采用与中心格式式（4.57）相同的思想。

文献 [42, 57] 提出了 AUSM$ ^{+} $，进一步改善基本 AUSM 在激波附近的表现。修改包括新的 Mach 数分裂和压力分裂，分别替代式（4.63）、（4.64）及（4.73）、（4.74）。

CUSP 格式

对流迎风分裂压力（CUSP）格式的概念与 AUSM 很相似。不过，CUSP 的优势是能写成通量平均（无 AUSM 的加权）减去耗散项；这一特点对在显式混合多级格式中的实现至关重要。此外，由于尺度因子不同，CUSP 在流动与网格对齐时表现优于 AUSM。CUSP 由 Jameson $ [43, 60, 61] $ 提出，后由 Tatsumi 等人 $ [44, 62] $ 修改，可在单元中心或（单元顶点）对偶控制体空间离散中实现。

穿过控制体面的对流通量（式（2.21）），分别按式（4.15）、（4.20）、（4.35）或（4.40）采用通量算术平均近似。随后，从中心通量中减去耗散项以稳定计算。因此，面 $ (I + 1/2) $ 上的总对流通量为：

$$ (\vec{F}_{\mathrm{c}})_{I+1/2}=\frac{1}{2}\left[\vec{F}_{\mathrm{c}}(\vec{W}_{\mathrm{R}})+\vec{F}_{\mathrm{c}}(\vec{W}_{\mathrm{L}})\right]-\vec{D}_{I+1/2}.\tag{4.78} $$

对于对偶控制体离散，改用 $ (i+1/2) $。耗散项由状态矢量差和通量矢量差的线性组合构成，可表示为：

$$ \begin{aligned}\vec{D}_{I+1/2}=\frac{1}{2}(\alpha^{*}c)_{I+1/2}&\left\{\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho w\\ \rho\phi\end{bmatrix}_{\mathrm{R}}-\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho w\\ \rho\phi\end{bmatrix}_{\mathrm{L}}\right\}\\+\frac{1}{2}\beta_{I+1/2}&\left\{\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}_{\mathrm{R}}-\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}_{\mathrm{L}}\right\}.\end{aligned}\tag{4.79} $$

式（4.79）中的 $ \phi $ 有两种选择：将 $ \phi $ 设为总能量 E，或设为总焓 H。前者称为 E-CUSP 格式 [63]，后者相应称为 H-CUSP 格式。采用 $ \phi = H $ 的表述保持总焓 [44]，因而适用于无黏流。左（L）、右（R）状态的计算与 MUSCL 方法 [30] 类似，采用带限制器的插值式（4.118）至（4.121）。式（4.79）的两个因子 $ \alpha^{*}c $ 和 $ \beta $ 定义为：

$$ \alpha^{*}c=\left\{\begin{array}{l l}{|V|}&{\text{当}\beta=0,}\\ {-(1+\beta)\Lambda^{-}}&{\text{当}\beta>0\text{且}0<M_{n}<1,}\\ {+(1-\beta)\Lambda^{+}}&{\text{当}\beta<0\text{且}-1<M_{n}<0}\\ {0}&{\text{当}|M_{n}|\geq1,}\end{array}\right.\tag{4.80} $$

以及：

$$ \beta=\left\{\begin{array}{l l}{+\max\left(0,\frac{V+\Lambda^{-}}{V-\Lambda^{-}}\right)}&{\text{当}0\leq M_{n}<1,}\\ {-\max\left(0,\frac{V+\Lambda^{+}}{V-\Lambda^{+}}\right)}&{\text{当}-1<M_{n}<0,}\\ {\operatorname{sign}(M_{n})}&{\text{当}|M_{n}|\geq1}\end{array}\right.\tag{4.81} $$

其中，$ M_{n}=V/c $。式（4.80）与（4.81）的逆变速度 V，通过速度分量的算术平均计算，即：

$$ V_{I+1/2}=\frac{1}{2}\left[\left(u_{I}+u_{I+1}\right)n_{x}+\left(v_{I}+v_{I+1}\right)n_{y}+\left(w_{I}+w_{I+1}\right)n_{z}\right].\tag{4.82} $$

计算 $ M_{n} $ 时，声速 c 也由算术平均量求得。

式（4.80）、（4.81）的正、负特征值 $ \Lambda^{+} $ 与 $ \Lambda^{-} $，是 Roe 矩阵 [64] 的特征值；该矩阵将在第 4.3.3 节讨论。文献 [61] 给出的特征值为：

$$ \Lambda^{\pm}=\frac{\gamma+1}{2\gamma}\tilde{V}\pm\sqrt{\left(\frac{\gamma-1}{2\gamma}\tilde{V}\right)^{2}+\frac{\tilde{c}^{2}}{\gamma}},\tag{4.83} $$

其中，$ \tilde{V} $ 为逆变速度，$ \gamma $ 为比热比，$ \tilde{c} $ 为声速。式（4.83）中需要在控制体面上计算的流动变量，由 Roe 平均 [64] 得到：

$$ \begin{aligned}\tilde{u}_{I+1/2}&=\frac{u_{\mathrm{L}}\sqrt{\rho_{\mathrm{L}}}+u_{\mathrm{R}}\sqrt{\rho_{\mathrm{R}}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}},\\\tilde{ v}_{I+1/2}&=\frac{ v_{\mathrm{L}}\sqrt{\rho_{\mathrm{L}}}+v_{\mathrm{R}}\sqrt{\rho_{\mathrm{R}}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}},\\\tilde{w}_{I+1/2}&=\frac{w_{\mathrm{L}}\sqrt{\rho_{\mathrm{L}}}+w_{\mathrm{R}}\sqrt{\rho_{\mathrm{R}}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}},\\\tilde{H}_{I+1/2}&=\frac{H_{\mathrm{L}}\sqrt{\rho_{\mathrm{L}}}+H_{\mathrm{R}}\sqrt{\rho_{\mathrm{R}}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}},\\\tilde{c}_{I+1/2}&=\sqrt{(\gamma-1)\left(\tilde{H}-\frac{\tilde{u}^{2}+\tilde{v}^{2}+\tilde{w}^{2}}{2}\right)_{I+1/2}},\\\tilde{V}_{I+1/2}&=\tilde{u}_{I+1/2}n_{x}+\tilde{v}_{I+1/2}n_{y}+\tilde{w}_{I+1/2}n_{z}.\end{aligned}\tag{4.84} $$

因子 $ \alpha^{*}c $ 和 $ \beta $ 的定义使超声速流的对流通量完全迎风，即 $ \alpha^{*}c=0 $ 和 $ \beta=\operatorname{sign}(M_{n}) $。另一方面，亚声速流（$ \beta=0 $）的耗散用 $ |V| $ 缩放，这有利于黏性层计算。对于长宽比较大的单元，显式时间推进通常需要增大沿单元长边方向的数值耗散，以保持稳健性；可采用类似式（4.57）的谱半径比实现。详见文献 [38]，其中还比较了 CUSP 与标量、矩阵人工耗散格式（第 4.3.1 节）。

将式（4.84）的 Roe 平均替换为算术平均，并将式（4.83）的特征值替换为 $ \Lambda^{\pm}=V\pm c $，可以显著简化 CUSP 的实现。引入定义：

$$ \alpha^{*}c=\alpha c-\beta~V,\tag{4.85} $$

则修改后格式的因子为 $ [44, 61] $：

$$ \alpha=\left\{\begin{array}{l l}{|M_{n}|}&{\text{当}|M_{n}|\geq\delta,}\\ {\displaystyle\frac{M_{n}^{2}+\delta^{2}}{2\delta}}&{\text{当}|M_{n}|<\delta,}\end{array}\right.\tag{4.86} $$

以及：

$$ \beta=\left\{\begin{array}{l l}{\max\left(0,2M_{n}-1\right)}&{\text{当}0\leq M_{n}<1,}\\ {\min\left(0,2M_{n}+1\right)}&{\text{当}-1<M_{n}<0,}\\ {\operatorname{sign}(M_{n})}&{\text{当}|M_{n}|\geq1.}\end{array}\right.\tag{4.87} $$

$ M_{n}=V/c $ 中的逆变速度仍按式（4.82）计算。式（4.86）的参数 $ \delta $ 旨在防止驻点处耗散消失，但似乎并非总有必要。这些简化得到一种计算效率很高的格式，不过，与采用 Roe 平均的原始表述相比，激波分辨能力略有下降。

4.3.3 通量差分分裂格式

通量差分分裂格式通过求解 Riemann（激波管）问题，根据通常不连续的左、右状态，计算控制体面上的对流通量。这一思想由 Godunov $ [65] $ 首先提出。与通量矢量分裂不同，通量差分分裂不仅考虑波（信息）的传播方向，还考虑波本身。为降低 Godunov 格式精确求解 Riemann 问题的计算量，研究者提出了近似 Riemann 求解器，例如 Osher 等人 $ [66] $ 和 Roe $ [64] $ 的方法。特别是 Roe 方法，因边界层流动计算精度高、激波分辨能力好，得到广泛应用。因此，下一部分将详细介绍 Roe 求解器。

Roe 迎风格式

Roe 近似 Riemann 求解器可在单元中心格式或对偶控制体格式中实现。它将控制体面上的通量差分解为各波贡献之和，同时确保 Euler 方程的守恒性质。在面 $ (I+1/2) $ 或 $ (i+1/2) $ 上，通量差写为 [64]：

$$ (\vec{F}_{\mathrm{c}})_{\mathrm{R}}-(\vec{F}_{\mathrm{c}})_{\mathrm{L}}=(\bar{A}_{\mathrm{Roe}})_{I+1/2}(\vec{W}_{\mathrm{R}}-\vec{W}_{\mathrm{L}}).\tag{4.88} $$

其中，$ \bar{A}_{Roe} $ 为 Roe 矩阵，L 和 R 分别表示左、右状态（见图 4.8）。Roe 矩阵与对流通量 Jacobian $ \bar{A}_{c} $（见第 A.9 节）相同，只是将流动变量替换为 Roe 平均变量。当 Roe 平均按下列公式由左、右状态计算时 [64, 67]，式（4.88）的通量差关系精确成立：

$$ \begin{aligned}\tilde{\rho}&=\sqrt{\rho_{\mathrm{L}}\rho_{\mathrm{R}}},\\\tilde{u}&=\frac{u_{\mathrm{L}}\sqrt{\rho_{\mathrm{L}}}+u_{\mathrm{R}}\sqrt{\rho_{\mathrm{R}}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}},\\\tilde{ v}&=\frac{ v_{\mathrm{L}}\sqrt{\rho_{\mathrm{L}}}+ v_{\mathrm{R}}\sqrt{\rho_{\mathrm{R}}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}},\\\tilde{w}&=\frac{w_{\mathrm{L}}\sqrt{\rho_{\mathrm{L}}}+w_{\mathrm{R}}\sqrt{\rho_{\mathrm{R}}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}},\\\tilde{H}&=\frac{H_{\mathrm{L}}\sqrt{\rho_{\mathrm{L}}}+H_{\mathrm{R}}\sqrt{\rho_{\mathrm{R}}}}{\sqrt{\rho_{\mathrm{L}}}+\sqrt{\rho_{\mathrm{R}}}},\\\tilde{c}&=\sqrt{(\gamma-1)\left(\tilde{H}-\tilde{q}^{2}/2\right)},\\\tilde{V}&=\tilde{u}n_{x}+\tilde{ v}n_{y}+\tilde{w}n_{z},\\\tilde{q}^{2}&=\tilde{u}^{2}+\tilde{ v}^{2}+\tilde{w}^{2}.\end{aligned}\tag{4.89} $$

将 Roe 矩阵的对角化形式 $ \bar{A}_{Roe} = \bar{T}\bar{\Lambda}_{c}\bar{T}^{-1} $ 代入式（4.88），可更清楚地说明 Roe 格式中的波分解：

$$ (\vec{F}_{\mathrm{c}})_{\mathrm{R}}-(\vec{F}_{\mathrm{c}})_{\mathrm{L}}=\bar{T}\bar{\Lambda}_{\mathrm{c}}(\vec{C}_{\mathrm{R}}-\vec{C}_{\mathrm{L}}).\tag{4.90} $$

左特征向量矩阵 $ (\bar{T}^{-1}) $、右特征向量矩阵 $ (\bar{T}) $ 和特征值对角矩阵 $ (\bar{\Lambda}_{c}) $，均采用 Roe 平均（4.89）计算。式（4.90）中，特征变量 $ \vec{C} $ 表示波幅，特征值 $ \Lambda_{c} $ 是近似 Riemann 问题的相应波速，右特征向量则是波本身。

根据上述讨论，控制体面上的对流通量按下式计算 $ [64] $：

$$ (\vec{F}_{c})_{I+1/2}=\frac{1}{2}\left[\vec{F}_{c}(\vec{W}_{\mathrm{R}})+\vec{F}_{c}(\vec{W}_{\mathrm{L}})-|\bar{A}_{\mathrm{Roe}}|_{I+1/2}(\vec{W}_{\mathrm{R}}-\vec{W}_{\mathrm{L}})\right].\tag{4.91} $$

$ \left|\bar{A}_{Roe}\right| $ 与左、右状态之差的乘积，可按下列方式高效计算：

$$ |\bar{A}_{\mathrm{Roe}}|(\vec{W}_{\mathrm{R}}-\vec{W}_{\mathrm{L}})=|\Delta\vec{F}_{1}|+|\Delta\vec{F}_{2,3,4}|+|\Delta\vec{F}_{5}|,\tag{4.92} $$

其中，

$$ |\Delta\vec{F}_{1}|=|\tilde{V}-\tilde{c}|\left(\frac{\Delta p-\tilde{\rho}\tilde{c}\Delta V}{2\tilde{c}^{2}}\right)\begin{bmatrix}1\\ \tilde{u}-\tilde{c}n_{x}\\ \tilde{ v}-\tilde{c}n_{y}\\ \tilde{w}-\tilde{c}n_{z}\\ \tilde{H}-\tilde{c}\tilde{V}\end{bmatrix},\tag{4.93} $$

$$ \begin{aligned}|\Delta\vec{F}_{2,3,4}|&=|\tilde{V}|\left\{\left(\Delta\rho-\frac{\Delta p}{\tilde{c}^{2}}\right)\begin{bmatrix}1\\ \tilde{u}\\ \tilde{ v}\\ \tilde{w}\\ \tilde{q}^{2}/2\end{bmatrix},\right.\\&\left.+\tilde{\rho}\begin{bmatrix}0\\ \Delta u-\Delta V n_{x}\\ \Delta v-\Delta V n_{y}\\ \Delta w-\Delta V n_{z}\\ \tilde{u}\Delta u+\tilde{ v}\Delta v+\tilde{w}\Delta w-\tilde{V}\Delta V\end{bmatrix}\right\},\end{aligned}\tag{4.94} $$

$$ |\Delta\vec{F}_{5}|=|\tilde{V}+\tilde{c}|\left(\frac{\Delta p+\tilde{\rho}\tilde{c}\Delta V}{2\tilde{c}^{2}}\right)\begin{bmatrix}1\\ \tilde{u}+\tilde{c}n_{x}\\ \tilde{ v}+\tilde{c}n_{y}\\ \tilde{w}+\tilde{c}n_{z}\\ \tilde{H}+\tilde{c}\tilde{V}\end{bmatrix}.\tag{4.95} $$

跳跃量定义为 $ \Delta(\bullet)=(\bullet)_{\mathrm{R}}-(\bullet)_{\mathrm{L}} $，Roe 平均变量见式（4.89）。

左、右状态由式（4.46）的 MUSCL 格式 $ [30] $ 确定。当流场存在间断时，全部高阶格式 $ (\hat{\kappa} = -1, \hat{\kappa} = 0, \text{及} \hat{\kappa} = 1/3) $ 都必须加入限制器（第 4.3.5 节）。

由于式（4.88）的表述，Roe 近似 Riemann 求解器在定常膨胀情况下会产生非物理的膨胀激波；此时 $ (\vec{F}_{\mathrm{c}})_{\mathrm{L}} = (\vec{F}_{\mathrm{c}})_{\mathrm{R}} $，但 $ \vec{W}_{L} \neq \vec{W}_{R} $。此外，还可能出现“激波鼓包”现象，即强弓形激波前方沿驻点线增长的扰动 [68, 69]，另见文献 [70] 的讨论。根本困难在于原始格式无法识别声速点。为解决这一问题，用 Harten 熵修正 [71, 72] 修改特征值 $ |\Lambda_{c}| = |\tilde{V}\pm\tilde{c}| $ 的绝对值：

$$ \left|\Lambda_{\mathrm{c}}\right|=\left\{\begin{array}{ll}\left|\Lambda_{\mathrm{c}}\right|&\text{当}\left|\Lambda_{\mathrm{c}}\right|>\delta\\\frac{\Lambda_{\mathrm{c}}^{2}+\delta^{2}}{2\delta}&\text{当}\left|\Lambda_{\mathrm{c}}\right|\leq\delta,\end{array}\right.\tag{4.96} $$

其中，$ \delta $ 为一个小量，可方便地取为当地声速的某一比例，例如 1/10 或 1/20。为避免线性波 $ |\Delta\vec{F}_{2,3,4}| $ 在 $ \tilde{V}\to0 $ 时消失（例如驻点或流动与网格对齐处），上述修改也可用于 $ |\tilde{V}| $。

在真实气体模拟中，Roe 求解器相对于中心格式或通量矢量分裂格式的一项明显缺点，是必须相应修改 Roe 矩阵及平均公式，这可能相当复杂。平衡和非平衡真实气体流动的表述实例，见文献 [73–76] 及其所引文献。文献 [77] 最近介绍了适用于任意可压缩及不可压缩流体的 Roe 格式实现。另请注意，第 9.5 节给出了经过低 Mach 数预处理的 Roe 格式（见式（9.82））。

4.3.4 总变差减小格式

总变差减小（TVD）格式的思想最早由 Harten $ [78] $ 研究。TVD 格式旨在防止流动解产生新的极值。对于标量守恒方程，解的总变差定义为：

$$ \mathrm{TV}\equiv\sum_{I}\left|U_{I+1}-U_{I}\right|\tag{4.97} $$

校注：本节原文将“对称 TVD”与“迎风 TVD”对特征速度符号的依赖关系写反，与第 3.1.5 节及本节后续说明不一致。通常依赖传播方向者为迎风格式；此处保留原文表述并注明这一疑点。

TVD 格式的主要条件是，总变差随时间减小。这意味着解的最大值不得增加，最小值不得减小，因此时间演化中不得产生新的局部极值。由此，具有 TVD 性质的离散能够准确分辨强激波，而不出现例如采用标量或矩阵人工耗散的中心格式（第 4.3.1 节）所产生的虚假振荡。

TVD 格式实现为对流通量平均加上满足 TVD 条件的耗散项（通量受限耗散）$ [78, 79] $。原文将依赖特征速度符号的耗散项称为对称 TVD 格式 $ [80, 81] $，将不依赖者称为迎风 TVD 格式 $ [82–86] $。经验表明，迎风 TVD 格式的精度优于对称 TVD 格式 $ [87] $，尤其适用于超声速和高超声速流场 $ [88] $，也能准确分辨边界层 $ [54] $，特别是在采用文献 $ [89] $ 的修改后。

迎风 TVD 格式

在这一框架中，穿过控制体面 $ (I+1/2) $（见图 4.8）的对流通量可表示为：

$$ (\vec{F}_{c})_{I+1/2}=\frac{1}{2}\left[(\vec{F}_{c})_{I+1}+(\vec{F}_{c})_{I}\right]+\frac{1}{2}\;\bar{T}_{I+1/2}\vec{\Theta}_{I+1/2}.\tag{4.98} $$

对于采用对偶控制体的单元顶点格式（第 4.2.3 节），相应下标为 $ (i+1/2) $、$ (i+1) $ 等。矩阵 $ \bar{T} $ 包含 Jacobian $ \bar{A}_{c} = \partial\vec{F}_{c}/\partial\vec{W} $ 的右特征向量，其元素见第 A.11 节。式（4.98）中的 $ \Theta $ 考虑特征速度的方向，

控制差分算子的迎风方向。矢量 $ \vec{\Theta} $ 的第 l 个分量定义为（参见 [85]）：

$$ \Theta_{I+1/2}^{l}=\frac{1}{2}\psi\left(\Lambda_{I+1/2}^{l}\right)\left(\Psi_{I+1}^{l}+\Psi_{I}^{l}\right)-\psi\left(\Lambda_{I+1/2}^{l}+\chi_{I+1/2}^{l}\right)\Delta C_{I+1/2}^{l},\tag{4.99} $$

其中，$ \Lambda^{l} $ 为对角矩阵 $ \bar{\Lambda}_{c} $ 的各特征值（见第 A.11 节），$ \Psi $ 为限制器函数（式（4.122））。此外，

$$ \chi_{I+1/2}^{l}=\frac{1}{2}\psi\left(\Lambda_{I+1/2}^{l}\right)\cdot\left\{\begin{array}{l l}{\frac{\left(\Psi_{I+1}^{l}-\Psi_{I}^{l}\right)}{\Delta C_{I+1/2}^{l}}}&{\text{当}\Delta C_{I+1/2}^{l}\neq0,}\\ {0}&{\text{当}\Delta C_{I+1/2}^{l}=0,}\end{array}\right.\tag{4.100} $$

最后，$ \Delta C^{l} $ 为特征变量差的各元素，即：

$$ \Delta\vec{C}_{I+1/2}=\bar{T}_{I+1/2}^{-1}\left(\vec{W}_{I+1}-\vec{W}_{I}\right)\tag{4.101} $$

其中，$ \bar{T}^{-1} $ 为左特征向量矩阵。Harten 熵修正 [71, 72]，即：

$$ \psi(z)=\left\{\begin{array}{ll}|z|&\text{当}|z|>\delta_{1},\\\displaystyle\frac{z^{2}+\delta_{1}^{2}}{2\delta_{1}}&\text{当}|z|\leq\delta_{1}\end{array}\right.\tag{4.102} $$

防止 $ \psi(z) $ 在 $ |z|\rightarrow0 $ 时变为零。参数 $ \delta_{1} $ 最好表示为速度分量及声速的函数 [85]：

$$ (\delta_{1})_{I+1/2}=\delta\left(|u_{I+1/2}|+| v_{I+1/2}|+|w_{I+1/2}|+c_{I+1/2}\right),\tag{4.103} $$

其中，$ 0.05 \leq \delta \leq 0.5 $。面 $ (I + 1/2) $ 上的原始变量，可由 I 与 $ (I + 1) $ 处的状态采用 Roe 平均（4.89）或简单算术平均获得。限制器函数 $ \Psi $ 防止强梯度附近出现虚假解，将在下一节介绍。必须强调，上述迎风 TVD 格式并未采用 MUSCL 方法获得高阶精度。

可以证明，当式（4.99）、（4.100）的限制器函数 $ \Psi $ 设为零时，迎风 TVD 方法恰为一阶空间精度 [90]；这发生在间断处。其他情况下，上述迎风 TVD 格式在光滑流动区域具有二阶精度。

4.3.5 限制器函数

二阶和高阶迎风空间离散需要采用限制器（限制器函数），以防止大梯度区域（如激波处）出现振荡和虚假解。因此，至少需要一个保单调格式，也就是要求流场中的最大值不得增加，最小值不得减小，且时间演化中不得产生新的局部极值。换言之，如果初始数据单调，解必须保持单调。保单调格式的条件相当严格，TVD 条件则更加严格，因此常用局部极值减小（LED）条件 $ [61] $ 代替：只要求模板范围内的局部极值减小。

然而，根据 Godunov 定理，高阶线性格式（例如 MUSCL 方法）不可能保单调 $ [91] $。因此，必须采用非线性限制器函数，构造保单调或 TVD 离散。图 4.10 展示了这一点：用式（4.98）的迎风 TVD 格式，分别在采用和不采用限制器的情况下计算 NACA 0012 翼型的二维跨声速绕流。可以清楚看到，无限制器时，翼型上下表面激波附近的解存在大幅振荡；远离激波时，受限解与未受限解则几乎相同。

限制器的作用，是减小将流动变量插值到控制体面时使用的斜率（即 $ (U_{I+1}-U_{I})/\Delta x $），从而约束解的变化。

在强间断处，限制器必须将斜率降为零，以防产生新的极值。对于 MUSCL 和 TVD 格式，这意味着在大梯度紧邻区域恢复（单调的）一阶迎风格式（式（4.46）中 $ \epsilon = 0 $）。限制器的最后一项要求很明显：在光滑流动区域，应恢复原来的无限制离散，使数值耗散尽可能低。图 4.11 示意了限制器对左、右状态插值的影响：它减小 I 处局部极小值位置的斜率，并改变单元 $ (I + 1) $、$ (I + 2) $ 的斜率，以得到单调解。必须认识到，面上左右状态之间仍可能存在差异，而且通常确实如此。

下面介绍四种成熟且经实践检验的限制器函数，分别考虑用于二阶 MUSCL、CUSP 和迎风 TVD 格式的限制器。

MUSCL 插值的限制器函数

Van Leer 的 MUSCL 方法 $ [30] $，通过在必要时用限制器减小式（4.47）的差分 $ \Delta_{+}U_{I} $ 和 $ \Delta_{-}U_{I} $，转化为保单调格式。引入斜率限制器 $ \Phi^{\pm} $ 后，式（4.46）的 MUSCL 插值修改为（另见图 4.8）：

$$ \begin{aligned}&U_{\mathrm{R}}=U_{I+1}-\frac{1}{4}\Big[(1+\hat{\kappa})\Phi_{I+1/2}^{+}\Delta_{-}+(1-\hat{\kappa})\Phi_{I+3/2}^{-}\Delta_{+}\Big]U_{I+1},\\&U_{\mathrm{L}}=U_{I}+\frac{1}{4}\Big[(1+\hat{\kappa})\Phi_{I+1/2}^{-}\Delta_{+}+(1-\hat{\kappa})\Phi_{I-1/2}^{+}\Delta_{-}\Big]U_{I}.\\ \end{aligned}\tag{4.104} $$

式（4.46）的参数 $ \epsilon $ 已设为 1。斜率限制器是相邻解变化量之比的函数，即 $ \Phi_{I+1/2}^{\pm} = \Phi(r_{I+1/2}^{\pm}) $，其定义为 [1]：

$$ \begin{aligned}&r_{I+1/2}^{+}=\frac{U_{I+2}-U_{I+1}}{U_{I+1}-U_{I}},\\&r_{I+1/2}^{-}=\frac{U_{I}-U_{I-1}}{U_{I+1}-U_{I}},\text{等}\\ \end{aligned}\tag{4.105} $$

现在，用 $ r_{L} $ 代替 $ r_{I-1/2}^{+} $，用 $ r_{R} $ 代替 $ r_{I+3/2}^{-} $，即：

$$ \begin{aligned}r_{\mathrm{R}}&=\frac{U_{I+1}-U_{I}}{U_{I+2}-U_{I+1}}=\frac{\Delta_{-}}{\Delta_{+}}U_{I+1},\\r_{\mathrm{L}}&=\frac{U_{I+1}-U_{I}}{U_{I}-U_{I-1}}=\frac{\Delta_{+}}{\Delta_{-}}U_{I},\end{aligned}\tag{4.106} $$

则式（4.104）可写为：

$$ \begin{aligned}&U_{\mathrm{R}}=U_{I+1}-\frac{1}{4}\left[(1+\hat{\kappa})r_{\mathrm{R}}\Phi(1/r_{\mathrm{R}})+(1-\hat{\kappa})\Phi(r_{\mathrm{R}})\right](U_{I+2}-U_{I+1}),\\&U_{\mathrm{L}}=U_{I}+\frac{1}{4}\left[(1+\hat{\kappa})r_{\mathrm{L}}\Phi(1/r_{\mathrm{L}})+(1-\hat{\kappa})\Phi(r_{\mathrm{L}})\right](U_{I}-U_{I-1}).\\ \end{aligned}\tag{4.107} $$

若只考虑具有下列对称性的斜率限制器，式（4.107）还可简化：

$$ \Phi(r)=\Phi(1/r).\tag{4.108} $$

由此，受限 MUSCL 插值式 $ (4.104) $ 变为 [92]：

$$ \begin{array}{l}U_{\mathrm{R}}=U_{I+1}-\frac{1}{2}\Psi_{\mathrm{R}}\left(U_{I+2}-U_{I+1}\right),\\U_{\mathrm{L}}=U_{I}+\frac{1}{2}\Psi_{\mathrm{L}}\left(U_{I}-U_{I-1}\right)\end{array}\tag{4.109} $$

限制器函数定义为：

$$ \Psi_{\mathrm{L/R}}=\frac{1}{2}\left[\left(1+\hat{\kappa}\right)r_{\mathrm{L/R}}+\left(1-\hat{\kappa}\right)\right]\Phi_{\mathrm{L/R}}.\tag{4.110} $$

式（4.110）的斜率限制器 $ \Phi $ 可以有不同形式，并针对 $ \hat{\kappa} $ 的特定值设计，从而得到尽可能精确、同时稳定且保单调的 MUSCL 格式。

$ \hat{\kappa}=0 $ 的 MUSCL 格式

对于 $ \hat{\kappa}=0 $ 的二阶偏迎风格式，一种尤其合适的组合为 [93]：

$$ \Phi(r)=\frac{2r}{r^{2}+1}.\tag{4.111} $$

此时，函数 $ \Psi(r) $ 对应 Van Albada 限制器 [94]：

$$ \Psi(r)=\frac{r^{2}+r}{1+r^{2}},\tag{4.112} $$

结合式（4.109），得到下列左、右状态表达式：

$$ \begin{aligned}&U_{\mathrm{R}}=U_{I+1}-\frac{1}{2}\delta_{\mathrm{R}},\\&U_{\mathrm{L}}=U_{I}+\frac{1}{2}\delta_{\mathrm{L}}.\\ \end{aligned}\tag{4.113} $$

函数 $ \delta $ 对两种状态在形式上相同，写为：

$$ \delta=\frac{a(b^{2}+\epsilon)+b(a^{2}+\epsilon)}{a^{2}+b^{2}+2\epsilon}.\tag{4.114} $$

系数 a、b 对左、右状态分别定义为：

$$ \begin{aligned}&a_{\mathrm{R}}=\Delta_{+}U_{I+1},b_{\mathrm{R}}=\Delta_{-}U_{I+1},\\&a_{\mathrm{L}}=\Delta_{+}U_{I},\quad b_{\mathrm{L}}=\Delta_{-}U_{I}\\ \end{aligned}\tag{4.115} $$

差分算子 $ \Delta_{\pm} $ 由式（4.47）给出。式（4.114）的附加参数 $ \epsilon $ 防止光滑流动区域内的小尺度振荡激活限制器 [93]；有时，这是获得完全收敛定常解所必需的。参数 $ \epsilon $ 宜取为与局部网格尺度成正比，例如三维中与 $ \Omega^{1/3} $ 成正比 [93, 95]。若状态变量 U 使用有量纲物理单位，还需对 $ \epsilon $ 作额外尺度调整。可以证明，在流动光滑变化时，式（4.113）与 $ \hat{\kappa}=0 $ 的原始无限制 MUSCL 格式（4.46）相同，因此不影响解的精度。另一方面，局部极值处的函数 $ \delta $ 变为零，将精度按要求降为一阶。

$ \hat{\kappa}=1/3 $ 的 MUSCL 格式

对于 $ \hat{\kappa}=1/3 $ 的三点二阶偏迎风 MUSCL 格式，研究者构造了另一种限制器。其斜率限制器为：

$$ \Phi(r)=\frac{3r}{2r^{2}-r+2}.\tag{4.116} $$

此时，函数 $ \Psi(r) $ 对应 Hemker 和 Koren 的限制器 [96]。沿用前一种情况的推导，得到面上左右状态 $ (I+1/2) $ 的公式，与式（4.113）相同，但此时 [93]：

$$ \delta=\frac{(2a^{2}+\epsilon)b+(b^{2}+2\epsilon)a}{2a^{2}+2b^{2}-ab+3\epsilon}.\tag{4.117} $$

系数 a、b 及参数 $ \epsilon $ 的定义保持不变。

CUSP 格式的限制器

在 CUSP 格式（第 4.3.2 节）中，按 $ [44] $ 以二阶精度计算左（L）、右（R）状态：

$$ \begin{array}{l}{U_{\mathrm{R}}=U_{I+1}-\frac{1}{2}L\left(\Delta U_{I+3/2},\Delta U_{I-1/2}\right),}\\ {U_{\mathrm{L}}=U_{I}+\frac{1}{2}L\left(\Delta U_{I+3/2},\Delta U_{I-1/2}\right),}\end{array}\tag{4.118} $$

其中，

$$ \begin{aligned}\Delta U_{I-1/2}&=U_{I}-U_{I-1},\\\Delta U_{I+3/2}&=U_{I+2}-U_{I+1}.\end{aligned}\tag{4.119} $$

式（4.118）、（4.119）中，U 表示依赖变量，$ L() $ 表示受限平均：

$$ L(\Delta_{1},\Delta_{2})=\frac{1}{2}\Psi(\Delta_{1},\Delta_{2})\left(\Delta_{1}+\Delta_{2}\right),\tag{4.120} $$

限制器本身定义为：

$$ \Psi(\Delta_{1},\Delta_{2})=1-\left|\frac{\Delta_{1}-\Delta_{2}}{\left|\Delta_{1}\right|+\left|\Delta_{2}\right|+\epsilon}\right|^{\sigma},\tag{4.121} $$

其中，$ \sigma $ 为正系数，通常取 2；常数 $ \epsilon $ 用于防止除以零（例如取 $ \epsilon = 10^{-20} $）。若 $ \Delta_{1} $ 与 $ \Delta_{2} $ 异号但大小相等，则限制器为 $ \Psi = 0 $，意味着左、右状态近似只有一阶精度。

应当指出，也可以用 $ \hat{\kappa}=0 $ 的 MUSCL 格式和式（4.113）至（4.115）的 Van Albada 限制器，替代上述关系。

TVD 格式的限制器

与前面不同，这里的限制器作用于特征变量 $ \vec{C} $，而非守恒变量或原始变量。文献 [85] 给出了一种特别合适的限制器函数：

$$ \Psi_{I}^{l}=\frac{\Delta C_{I-1/2}^{l}\Delta C_{I+1/2}^{l}+|\Delta C_{I-1/2}^{l}\Delta C_{I+1/2}^{l}|}{\Delta C_{I-1/2}^{l}+\Delta C_{I+1/2}^{l}+\epsilon},\tag{4.122} $$

其中，$ \Delta C_{I+1/2}^{l} $ 为控制体面 $ (I+1/2) $ 上的特征变量差（式（4.101））。分母中的正常数 $ \epsilon\approx10^{-20} $ 防止除以零。在高梯度区域，限制器函数变为零，结合式（4.99）和（4.98），得到一阶迎风格式。在流动变量光滑变化的区域，$ \Psi_{I}^{l}=C_{I}^{l}-C_{I-1}^{l} $，式（4.98）的迎风 TVD 格式保持二阶精度。