5.3 对流通量的离散

5.3 对流通量的离散

前几节讨论了空间离散方法的一般问题及所需数据结构。下文进一步介绍对流通量计算的具体实现。

如 3.1.5 节所述，在有限体积方法框架内，主要可选择以下格式：

• 中心格式；

• 通量矢量分裂格式；

• 通量差分分裂格式；

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

• 波动分裂格式。

首先较详细地介绍非结构网格上的中心离散，因为它与结构网格上的中心离散有明显不同。迎风格式的基本原理在结构网格与非结构网格上则相同，具体见 4.3.2–4.3.4 节。不过，非结构网格新增的内容是解的重构，它用于获得控制体面上的流动变量。因此，5.3.3 节将较详细地讨论常用重构方法。限于篇幅，本章不讨论波动分裂方法，其相关文献见 3.1.5 节。

5.3.1 采用人工耗散的中心格式

中心格式的基本思想，是按式（5.20）对控制体面两侧的守恒变量取算术平均，再计算面上的对流通量。这样会造成解的奇偶解耦，即离散方程产生两组独立的解，并在激波处产生振荡，因此必须增加人工耗散以保持稳定。人工耗散由二阶差分与四阶差分混合组成。Jameson 等首先在结构网格上实现了用于 Euler 方程的这一格式 $ [44] $。根据作者姓名的首字母，它也简称 JST 格式。

JST 格式在非结构网格上的实现，以 Laplace 算子表示二阶差分，以 Laplace 算子的再次作用表示四阶差分 $ [27, 45] $。为降低计算量，用伪 Laplace 算子代替真正的 Laplace 算子。基于这一简化，文献 $ [46] $ 首先提出二维形式，文献 $ [47] $ 随后作了改进，后来又由文献 $ [2] $ 扩展至三维。该格式采用距离加权，使任意网格上线性变化函数的伪 Laplace 算子均为零。对单元 I 中一般标量 U，伪 Laplace 算子为：

$$ L(U_{I})=\sum_{J=1}^{N_{A}}\theta_{I J}\left(U_{J}-U_{I}\right),\tag{5.24} $$

其中，$ N_{A} $ 为相邻控制体的个数。对中位对偶格式，需将单元下标换成节点下标 $ (i,j) $。式（5.24）的求和最适合像通量计算那样，通过遍历面（单元中心格式）或遍历边（中位对偶格式）的循环实现。几何权重 $ \theta $ 定义为：

$$ \theta_{I J}=1+\Delta\theta_{I J}\tag{5.25} $$

这些权重由求解优化问题得到 $ [2] $。优化问题采用 Lagrange 乘子法求解，因此几何权重可由下式获得：

$$ \Delta\theta_{I J}=\lambda_{x,I}(x_{J}-x_{I})+\lambda_{ y,I}( y_{J}- y_{I})+\lambda_{z,I}(z_{J}-z_{I}),\tag{5.26} $$

其中，x、y、z 为单元形心的 笛卡尔 坐标；对中位对偶格式，则为节点坐标。每个单元（节点）的 Lagrange 乘子 $ \lambda $ 按下式计算 [2]：

$$ \begin{aligned}\lambda_{x}&=\frac{R_{x}a_{11}+R_{ y}a_{12}+R_{z}a_{13}}{d},\\\lambda_{ y}&=\frac{R_{x}a_{21}+R_{ y}a_{22}+R_{z}a_{23}}{d},\\\lambda_{z}&=\frac{R_{x}a_{31}+R_{ y}a_{32}+R_{z}a_{33}}{d}\end{aligned}\tag{5.27} $$

其中的系数为：

$$ \begin{aligned}a_{11}&=I_{yy}I_{zz}-I_{yz}^{2},\\a_{12}&=I_{xz}I_{yz}-I_{xy}I_{zz},\\a_{13}&=I_{xy}I_{yz}-I_{xz}I_{yy},\\a_{21}&=I_{xz}I_{yz}-I_{xy}I_{zz},\\a_{22}&=I_{xx}I_{zz}-I_{xz}^{2},\\a_{23}&=I_{xy}I_{xz}-I_{xx}I_{yz},\\a_{31}&=I_{xy}I_{yz}-I_{xz}I_{yy},\\a_{32}&=I_{xz}I_{xy}-I_{xx}I_{yz},\\a_{33}&=I_{xx}I_{yy}-I_{xy}^{2},\\d&=I_{xx}I_{yy}I_{zz}-I_{xx}I_{yz}^{2}-I_{yy}I_{xz}^{2}-I_{zz}I_{xy}^{2}+2I_{xy}I_{xz}I_{yz}.\end{aligned}\tag{5.28} $$

对单元 I，一阶矩为：

$$ \begin{aligned}R_{x,I}&=\sum_{J=1}^{N_{A}}(x_{J}-x_{I}),\\R_{ y,I}&=\sum_{J=1}^{N_{A}}( y_{J}- y_{I}),\\R_{z,I}&=\sum_{J=1}^{N_{A}}(z_{J}-z_{I}).\end{aligned}\tag{5.29} $$

二阶矩为：

$$ \begin{aligned}I_{x x,I}&=\sum_{J=1}^{N_{A}}(x_{J}-x_{I})^{2},\\I_{ y y,I}&=\sum_{J=1}^{N_{A}}( y_{J}- y_{I})^{2},\\I_{z z,I}&=\sum_{J=1}^{N_{A}}(z_{J}-z_{I})^{2},\\&I_{x y,I}=\sum_{J=1}^{N_{A}}(x_{J}-x_{I})( y_{J}- y_{I}),\\I_{xz,I}&=\sum_{J=1}^{N_{A}}(x_{J}-x_{I})(z_{J}-z_{I}),\\I_{ y z,I}&=\sum_{J=1}^{N_{A}}( y_{J}- y_{I})(z_{J}-z_{I}).\\ \end{aligned}\tag{5.30} $$

在严重畸变的网格上，几何权重 $ (5.25) $ 可能使 Laplace 算子的近似不满足正性，进而失去稳定性。因此，建议将权重截断至区间 $ [0,2] $，见文献 $ [46] $。不过，这会损害离散精度。进一步讨论亦见文献 $ [12] $。

四阶差分用 Laplace 算子的再次作用计算，即将式（5.24）中的 U 换为 $ L(U) $。因此，单元 I 的人工耗散项最终为：

$$ \begin{aligned}\vec{D}_{I}=&\sum_{J=1}^{N_{A}}(\hat{\Lambda}_{c})_{IJ} \epsilon_{I J}^{(2)}\theta_{I J}\left(\vec{W}_{J}-\vec{W}_{I}\right)\\&-\sum_{J=1}^{N_{A}}(\hat{\Lambda}_{c})_{IJ} \epsilon_{I J}^{(4)}\theta_{I J}\left[L(\vec{W}_{J})-L(\vec{W}_{I})\right].\end{aligned}\tag{5.31} $$

加上人工耗散项后，式（5.2）的方程组成为：

$$ \Omega_{I}\frac{\mathrm{d}\vec{W}_{I}}{\mathrm{d}t}=-\left[\sum_{m=1}^{N_{F}}(\vec{F}_{c}-\vec{F}_{ v})_{m}\Delta S_{m}\right]+\vec{D}_{I}+\vec{Q}_{I}\Omega_{I},\tag{5.32} $$

其中，$ N_{F} $ 表示控制体的面数；这一数量可能与相邻控制体数不同，例如将一个四边形面分成两个三角形时。

式（5.31）的二阶项和四阶项由对流通量 Jacobian 的谱半径缩放。根据文献 [28]，单元 I 的谱半径可按下式计算：

$$ (\hat{\Lambda}_{c})_{I}=\sum_{m=1}^{N_{F}}\left(|V_{m}|+c_{m}\right)\Delta S_{m},\tag{5.33} $$

其中，$ V_{m} $ 为逆变速度（2.22），$ c_{m} $ 为声速。二者均由面上的平均流动变量计算。控制体面上的谱半径为：

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

用基于压力的传感器，在激波处关闭四阶差分，在流场光滑区域关闭二阶差分。因此，式（5.31）中的系数 $ \epsilon_{IJ}^{(2)} $ 和 $ \epsilon_{IJ}^{(4)} $ 定义为：

$$ \begin{aligned}\epsilon_{IJ}^{(2)}&=k^{(2)}\max(\Upsilon_{I},\Upsilon_{J}),\\\epsilon_{IJ}^{(4)}&=\max\left[0,(k^{(4)}-\epsilon_{IJ}^{(2)})\right]\end{aligned}\tag{5.35} $$

压力传感器为：

$$ \Upsilon_{I}=\frac{\left|\sum_{J=1}^{N_{A}}\theta_{I J}\left(p_{J}-p_{I}\right)\right|}{\sum_{J=1}^{N_{A}}(p_{J}+p_{I})}.\tag{5.36} $$

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

如 4.3.1 节所述，若将式（5.31）的谱半径 $ (\hat{\Lambda}_{c})_{IJ} $ 换为矩阵，可提高上述格式的精度。这种所谓矩阵耗散格式 [48] 在非结构网格上的实现与结构网格相同，缩放矩阵按式（4.59）定义。矩阵耗散格式在三维混合单元网格上的应用，例如见文献 [49]。

应注意，对三角形／四面体以外的单元，常用的显式 Runge–Kutta 时间离散与中心格式结合时，会出现严重的稳定性问题 $ [50] $。原因在于用 Laplace 算子的再次作用表示四阶差分。改进方法是用左右状态之差（参见 4.3 节）近似四阶差分 $ [50] $，即：

$$ \vec{D}_{I}=\sum_{J=1}^{N_{A}}(\hat{\Lambda}_{c})_{I J}\epsilon_{I J}^{(2)}\theta_{I J}\left(\vec{W}_{J}-\vec{W}_{I}\right)-\sum_{J=1}^{N_{A}}4(\hat{\Lambda}_{c})_{I J}\epsilon_{I J}^{(4)}\left(\vec{W}_{\mathrm{L}}-\vec{W}_{\mathrm{R}}\right).\tag{5.37} $$

在四边形／六面体网格上，该方法得到的模板与相应结构网格格式相同。左右状态可采用 5.3.3 节所述线性重构等方法计算。

5.3.2 迎风格式

至少就目前而言，非结构网格上的迎风格式似乎比上述中心格式应用广泛得多。事实上，Roe 通量差分分裂格式 $ [51] $ 是非结构网格上使用最广泛的方法。与中心格式相比，它能更准确地分辨边界层，且对网格畸变较不敏感，这解释了 Roe 格式的吸引力。不过，性能改善的代价是计算量增加；若还必须使用限制器抑制解的振荡（5.3.5 节），增加的计算量会相当显著。

4.3 节介绍的结构网格迎风格式都可用于非结构网格，无需改变基本方法。只有左右状态的计算（式（5.22）），即解的重构，以及限制函数的计算，需要新的表达形式。因此，这里仅讨论重构和限制器。各种迎风方法的具体内容见 4.3.2–4.3.4 节。采用中位对偶方法在非结构网格上实现 Roe 格式的实例，例如见文献 $ ^{[33]} $。

5.3.3 解的重构

如 4.3.2–4.3.4 节所示，迎风格式需要指定控制体面左右两侧的流动状态。式（5.37）的修正人工耗散格式也有同样要求。

最简单的处理是假设解在各控制体内部为常量。此时，左右状态就是左右控制体中计算所得的流动变量。例如，对中位对偶格式（图 5.9），可令：

$$ \begin{array}{l}U_{\mathrm{L}}=U_{i},\\U_{\mathrm{R}}=U_{j}\end{array}\tag{5.38} $$

其中，U 为某个标量流动变量。这种空间离散仅有一阶精度。对黏性流动，一阶精度解的耗散过大，会导致剪切层过度增长，因此计算黏性流动必须采用更高空间精度的方法。

若允许解在控制体内部变化，就能获得二阶及更高阶的空间精度。二阶方法是最常用的高阶方法，其假设解在控制体内线性变化。为计算左右状态，必须重构所假设的解分布。下文讨论最常用的线性与二次重构方法。各种线性重构技术的比较见文献 $ [39] $。

基于 MUSCL 方法的重构

获得二阶精度的一种办法，是将 MUSCL 方法 $ [52] $ 扩展至非结构网格。用于中位对偶格式时，该方法为每条边 ij 生成两个“虚拟”节点 $ i^{\prime} $ 和 $ j^{\prime} $ [53–57]。如图 5.12 所示，将边 ij 向两端各延伸一个边长，两个虚拟节点位于延长线的端点。将周围单元（图 5.12 中的灰色单元）的解插值至虚拟节点后，即可按 MUSCL 公式（4.46）计算左右状态。因此有：

$$ \begin{aligned}U_{\mathrm{R}}&=U_{j}-\frac{1}{4}\left[(1+\hat{\kappa})\Delta_{-}+(1-\hat{\kappa})\Delta_{+}\right]U_{j},\\U_{\mathrm{L}}&=U_{i}+\frac{1}{4}\left[(1+\hat{\kappa})\Delta_{+}+(1-\hat{\kappa})\Delta_{-}\right]U_{i}\end{aligned}\tag{5.39} $$

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

$$ \begin{aligned}&\Delta_{+}U_{i}=U_{j}-U_{i},\quad\Delta_{-}U_{i}=U_{i}-U_{i^{\prime}},\\&\Delta_{+}U_{j}=U_{j^{\prime}}-U_{j},\quad\Delta_{-}U_{j}=U_{j}-U_{i}.\\ \end{aligned}\tag{5.40} $$

存在强间断时，必须按 4.3.5 节的方法，为 MUSCL 插值 $ (5.39) $ 增加限制函数。

该方法的一个缺点，是必须为每条边存储包含虚拟节点的单元。另一个概念上的缺点，是每个控制体未重构出唯一的梯度。此外，在边界处，其中一个虚拟点可能位于物理域外，从而造成困难。

分片线性重构

Barth 和 Jespersen 在文献 $ [30] $ 中提出一种与有限元格式密切相关的重构方法。这里假设解在控制体内分片线性分布。于是，单元中心格式的左右状态可由下式求得：

$$ \begin{aligned}U_{\mathrm{L}}&=U_{I}+\Psi_{I}\left(\nabla U_{I}\cdot\vec{r}_{\mathrm{L}}\right),\\U_{\mathrm{R}}&=U_{J}+\Psi_{J}\left(\nabla U_{J}\cdot\vec{r}_{\mathrm{R}}\right),\end{aligned}\tag{5.41} $$

其中，$ \nabla U_{I} $ 为单元中心 I 处 $ U\left(=[\partial U/\partial x,\partial U/\partial y,\partial U/\partial z]^{T}\right) $ 的梯度，$ \Psi $ 为限制函数（参见 5.3.5 节）。矢量 $ \vec{r}_{L} $ 与 $ \vec{r}_{R} $ 从单元形心指向面中点，见图 5.13(a)。

同一方法也适用于中位对偶格式 $ [30] $，即：

$$ \begin{aligned}U_{\mathrm{L}}&=U_{i}+\frac{1}{2}\Psi_{i}\left(\nabla U_{i}\cdot\vec{r}_{ij}\right),\\U_{\mathrm{R}}&=U_{j}-\frac{1}{2}\Psi_{j}\left(\nabla U_{j}\cdot\vec{r}_{ij}\right).\end{aligned}\tag{5.42} $$

根据图 5.9 或图 5.13(b)，

$$ \vec{r}_{ij}=\vec{r}_{j}-\vec{r}_{i}\tag{5.43} $$

表示从节点 i 指向节点 j 的矢量。

容易看出，Barth–Jespersen 方法相当于在面两侧的相邻单元中心／节点处进行 Taylor 展开，并只保留线性项。线性重构在规则网格上具有形式二阶精度 $ [39] $。只要梯度 $ \nabla U $ 的计算没有误差，该格式便能在任意网格上精确重构线性函数。线性重构可能是重构方法中应用最广泛的一种。

上述格式需要计算单元中心或节点处的梯度，可采用下文 5.3.4 节介绍的 Green–Gauss 方法或最小二乘方法。非结构网格上限制函数的实现将在 5.3.5 节详细介绍。

基于节点加权过程的线性重构

Frink [26] 证明，对单元中心格式，在纯三角形或纯四面体网格上，线性重构 $ (5.41) $ 无需显式计算梯度。原因在于这类单元具有两个不变的几何性质。第一，从节点穿过单元形心的直线总会经过对面中点。第二，单元形心至面中点的距离，是面中点至对面节点距离的四分之一；三角形则为三分之一。因此，单元中心的梯度可用简单有限差分近似 [26]。例如，在图 5.7(a) 的面中点 $ F_{3} $ 重构解时，式 $ (5.41) $ 成为：

$$ U_{L/R}=U_{\mathrm{C}}+\frac{\Psi_{\mathrm{C}}}{4}\left[\frac{1}{3}(U_{1}+U_{2}+U_{4})-U_{3}\right]\tag{5.44} $$

其中，$ U_{C} $ 为单元形心处的数值，$ U_{1} $、$ U_{2} $ 等为节点值，$ \Psi $ 为限制器。

Frink 提出两种不同方法确定节点值。第一种基于距离倒数加权：周围单元对某个节点的贡献，与该节点至单元形心的距离 $ [26, 58] $ 成反比，即：

$$ U_{i}=\left(\sum_{J=1}^{N_{A}}\theta_{i J}U_{J}\right)\Bigg/\left(\sum_{J=1}^{N_{A}}\theta_{i J}\right),\tag{5.45} $$

其中，权重为 $ \theta_{iJ}=1/r_{iJ} $。距离按下式计算：

$$ r_{i J}=\sqrt{(x_{J}-x_{i})^{2}+( y_{J}- y_{i})^{2}+(z_{J}-z_{i})^{2}}.\tag{5.46} $$

下标 J 和 i 分别表示单元形心与节点。该方法所得重构的精度低于二阶。不过，Frink 指出，至少对无黏流动不必使用限制器 $ [26] $，因而可大幅降低计算量。

第二种方法基于 Holmes 和 Connell [46] 以及 Rausch 等 [47] 的二维工作，后来由 Frink [2] 扩展至三维。这里对式（5.45）的权重 $ \theta_{ij} $ 作如下定义：若解线性变化，计算出的节点值必须精确。这与伪 Laplace 算子（5.24）的计算要求相同，因此权重也相同，由式（5.25）至（5.30）给出。

只需将单元形心坐标 $ x_{I} $、$ y_{I} $、$ z_{I} $ 替换为节点坐标 $ x_{i} $、$ y_{i} $、$ z_{i} $。由于能精确计算线性函数的节点值，格式具有形式二阶精度。为在畸变网格上保证正性，仍须将权重限制在区间 [0, 2] [46]，但这会降低重构精度。Frink 和 Pirzadeh [4] 还报告了这种重构用于 Navier–Stokes 方程时出现的某些异常表现。

分片二次重构

若用多项式重构获得高于二阶的精度，必须在面两侧相邻单元中心／节点处的截断 Taylor 展开中保留更多项。基于 Barth 和 Frederickson 的工作 $ [59] $，Barth 提出 k 阶精确重构格式的概念 $ [60] $，即对 k 次多项式精确的重构。Barth 方法定义多项式时保证均值守恒，换言之，重构多项式的平均值等于控制体中的平均解。这样可确保重构过程中质量、动量和能量守恒。该方法以 k=3 在中位对偶格式中实现，多项式系数通过最小二乘方法求得。Mitchell 和 Walters $ [61] $ 以及 Mitchell $ [62] $ 也在单元中心格式中采用了类似思路。不过，这些方法的计算量过大，数据结构复杂，因而未能广泛应用。

Delanaye 和 Essers $ [63] $ 以及 Delanaye $ [64] $ 为单元中心格式发展了一种特殊的二次重构，比 Barth 方法的计算效率更高。左右状态采用截断至二次项的 Taylor 展开近似 $ [63, 64] $：

$$ \begin{aligned}U_{\mathrm{L}}&=U_{I}+\Psi_{I,1}\left(\nabla U_{I}\cdot\vec{r}_{\mathrm{L}}\right)+\frac{1}{2}\Psi_{I,2}\left(\vec{r}_{\mathrm{L}}^{T}\bar{H}_{I}\vec{r}_{\mathrm{L}}\right),\\U_{\mathrm{R}}&=U_{J}+\Psi_{J,1}\left(\nabla U_{J}\cdot\vec{r}_{\mathrm{R}}\right)+\frac{1}{2}\Psi_{J,2}\left(\vec{r}_{\mathrm{R}}^{T}\bar{H}_{J}\vec{r}_{\mathrm{R}}\right).\end{aligned}\tag{5.47} $$

其中，$ \bar{H}_{I} $ 为 Hessian（海森）矩阵，即：

$$ \bar H_I=\begin{bmatrix}\partial_{xx}^2U&\partial_{xy}^2U&\partial_{xz}^2U\\\partial_{xy}^2U&\partial_{yy}^2U&\partial_{yz}^2U\\\partial_{xz}^2U&\partial_{yz}^2U&\partial_{zz}^2U\end{bmatrix}_I,\tag{5.48} $$

该矩阵在单元形心 I 处计算。变量 $ \Psi_{I,1} $ 和 $ \Psi_{I,2} $ 分别为线性项与二次项使用的两种不同限制函数 [63]。二次重构在规则网格上具有三阶精度；在任意网格上，由于误差项抵消，至少具有二阶精度 [64]。不过，要达到这些性质，式（5.47）中的梯度 $ \nabla U $ 必须至少以二阶精度计算，Hessian 矩阵必须至少以一阶精度计算。为此，将 Green–Gauss 梯度计算与基于最小二乘的二阶导数近似结合 [63, 64]，得到计算效率较高的格式。

但与线性重构相比，其内存和 CPU 时间的额外开销仍相当显著。该方法采用由面邻居和节点邻居组成的固定模板。图 5.14 同时显示了模板和 Green–Gauss 梯度计算的积分路径。要确定二次多项式的全部系数，模板至少须提供六个数值，三维则为十个。为保持二次重构的精度，必须允许解沿面线性变化，而不能视为常量。这意味着二维面上必须在两个所谓 Gauss 求积点重构解（参见图 5.14），三角形面则需三个点，并在控制体面上分段计算通量积分 [38]。

如今，二次及更高阶重构的研究重点转向 ENO（本质无振荡），尤其是 WENO（加权 ENO）格式 $ [65–85] $。其思想是采用方向各异的离散模板，为每个控制体构造多个不同的重构多项式。ENO 格式选择其中振荡最小的多项式；WENO 方法则对各多项式加权组合，权重与振荡程度成反比。这样，解在光滑区域具有高阶空间精度，同时又不会将真实间断抹平。WENO 格式的良好入门介绍见文献 $ [73] $。

5.3.4 梯度的计算

分片线性与二次重构的讨论还留下一个问题：如何确定梯度。计算黏性通量也需要速度分量和温度的梯度（5.4 节）。以下介绍两种方法：第一种基于 Green–Gauss 定理，第二种采用最小二乘法。

Green–Gauss 方法

该方法将标量函数 U 的梯度，近似为 U 与外指单位法向量乘积在某个控制体 $ \Omega' $ 表面上的积分，即：

$$ \nabla U\approx\frac{1}{\Omega^{\prime}}\int_{\partial\Omega^{\prime}}U\vec{n}\mathrm{d}S.\tag{5.49} $$

中位对偶格式

Barth 和 Jespersen $ [30] $ 从 Galerkin 有限元方法导出了 Green–Gauss 方法的一种离散形式，后来由 Barth $ [86] $ 扩展至三维。他们将式 $ (5.49) $ 应用于共享某节点的单元所组成的区域，并证明这种方法可写成与基于边的数据结构相容的形式。不过，这仅适用于三角形／四面体网格上的中位对偶格式。所得公式为：

$$ \nabla U_{i}\approx\frac{1}{\Omega}\sum_{j=1}^{N_{F}}\frac{1}{2}(U_{i}+U_{j})\vec{n}_{i j}\Delta S_{i j}.\tag{5.50} $$

这里，式（5.49）中的 $ \Omega^{\prime} $ 等于中位对偶控制体的体积 $ \Omega $。求和遍及与节点 i 相连的全部 $ N_{F} $ 条边。此外，$ \vec{n}_{ij} $ 为式（5.23）的平均单位法向量，$ \Delta S_{ij} $ 为总面面积。式（5.50）同时适用于二维和三维。应注意，在边界处必须修改求和方式，才能得到相容的近似 [87]，亦见 8.10 节。

单元中心格式

单元中心格式也可采用 Green–Gauss 方法。单元形心 I 处的梯度为：

$$ \nabla U_{I}\approx\frac{1}{\Omega}\sum_{J=1}^{N_{F}}\frac{1}{2}(U_{I}+U_{J})\vec{n}_{I J}\Delta S_{I J},\tag{5.51} $$

求和遍及体积为 $ \Omega $ 的单元的所有面。式（5.51）中，$ \vec{n}_{IJ} $ 为单位法向量，$ \Delta S_{IJ} $ 为面面积。

混合单元网格

采用式（5.50）或式（5.51）计算 Green–Gauss 梯度的主要吸引力，是其与通量计算（例如式（5.19））相似，因此无需增加数据结构即可重构梯度。主要缺点是，这两种近似在混合单元网格上会失效。文献 [50] 证明，梯度可能严重失准，不同单元类型交接处尤其如此。对中位对偶格式，可将式（5.49）的控制体 $ \Omega^{\prime} $ 仍取为与节点 i 相连的全部单元的并集，以解决这一问题。对图 5.8 的情形，二维梯度为：

$$ \nabla U_{i}\approx\frac{1}{\Omega^{\prime}}\sum_{j=1}^{N_{O}}\frac{1}{2}(U_{j}+U_{j+1})\vec{n}_{j}\Delta S_{j}\tag{5.52} $$

其中，i=0，外侧面数 $ N_{O}=6 $ 由点 $ P_{1}-P_{6} $ 确定，且 j=6 时 $ j+1=1 $。此外，$ \vec{n}_{j} $ 和 $ \Delta S_{j} $ 分别为单元外侧面的单位法向量与面积。在三维情况下，可采用：

$$ \nabla U_{i}\approx\frac{1}{\Omega^{\prime}}\sum_{j=1}^{N_{O}}\frac{1}{3}(U_{j,1}+U_{j,2}+U_{j,3})\vec{n}_{j}\Delta S_{j},\tag{5.53} $$

这里假设全部面均为三角形，或本来如此，或由分解得到。单元中心格式也可采用同样的处理：控制体 $ \Omega^{\prime} $ 的表面由第一层和第二层相邻单元的形心定义 [63, 64]，见图 5.14。梯度按式（5.52）或式（5.53）计算，只需将节点下标换成单元下标。

这种基于单元的方法有一个明显缺点：必须增加数据结构，以建立中心节点／形心与控制体 $ \Omega^{\prime} $ 外侧面的联系。因此，它不再具有网格透明性，即不再独立于单元信息，也无法采用高效的收集—散布循环。这使下述最小二乘技术在混合单元网格上更具吸引力。

在三角形或四面体网格上采用基于边／面的实现 $ (5.50) $ 或 $ (5.51) $，在混合单元网格上采用基于单元的方法 $ (5.52) $ 或 $ (5.53) $，Green–Gauss 方法至少具有一阶精度 $ [64] $，且满足相容性，即线性函数梯度的计算误差仅为舍入误差。一阶精度已足以用于线性重构。二次重构需要任意网格上的二阶精度，可从一阶梯度近似中减去估计的截断误差获得 $ [64] $。

最小二乘方法

Barth $ [37, 86] $ 首先提出用最小二乘方法计算梯度。以中位对偶格式说明：对与中心节点 i 相连的每条边，采用一阶 Taylor 展开近似。沿边 ij 的解变化可计算为：

$$ \left(\nabla U_{i}\right)\cdot\vec{r}_{i j}=U_{j}-U_{i},\tag{5.54} $$

其中，$ \vec{r}_{ij} $ 由式（5.43）给出，是从节点 i 指向节点 j 的矢量（见图 5.9 或图 5.13(b)）。将式（5.54）用于与节点 i 相连的全部边，得到以下超定线性方程组：

$$ \left[\begin{array}{c c c}\Delta x_{i1}&\Delta y_{i1}&\Delta z_{i1}\\ \Delta x_{i2}&\Delta y_{i2}&\Delta z_{i2}\\ \vdots&\vdots&\vdots\\ \Delta x_{i j}&\Delta y_{i j}&\Delta z_{i j}\\ \vdots&\vdots&\vdots\\ \Delta x_{i N_{A}}&\Delta y_{i N_{A}}&\Delta z_{i N_{A}}\end{array}\right]\left[\begin{array}{c}\partial_{x}U\\ \partial_{ y}U\\ \partial_{z}U\end{array}\right]_{i}=\left[\begin{array}{c}\theta_{1}\left(U_{1}-U_{i}\right)\\ \theta_{2}\left(U_{2}-U_{i}\right)\\ \vdots\\ \theta_{j}\left(U_{j}-U_{i}\right)\\ \vdots\\ \theta_{N_{A}}(U_{N_{A}}-U_{i})\end{array}\right]\tag{5.55} $$

其中，$ \Delta(\cdot)_{ij} = (\cdot)_j - (\cdot)_i $，$ \partial_m(\cdot) = \partial(\cdot)/\partial m $。此外，$ N_A $ 为通过边与 i 相连的相邻节点 j 的个数，$ \theta_j $ 为某个加权系数。权重可以依赖几何和／或解，例如见文献 [41]；但实际通常将 $ \theta_j $ 设为 1。为方便，将方程组（5.55）简写为：

$$ \bar{A}\vec{x}=\vec{b}.\tag{5.56} $$

由式（5.56）求梯度矢量 $ \vec{x} $，需要对矩阵 $ \vec{A} $ 求逆。为避免病态问题，尤其是拉伸网格上的病态问题，Anderson 和 Bonhaus 建议用 Gram–Schmidt 过程，将 $ \vec{A} $ 分解为正交矩阵 $ \bar{Q} $ 与上三角矩阵 $ \bar{R} $ 的乘积 [88]。文献 [50] 最近将这一方法扩展至三维。因此，式（5.56）的解直接由下式得到：

$$ \vec{x}=\bar{R}^{-1}\bar{Q}^{T}\vec{b}.\tag{5.57} $$

用带双下标的小写字母表示矩阵元素，矩阵 $ \vec{A} = [\vec{a}_{1}, \vec{a}_{2}, \vec{a}_{3}] $ 的 Gram–Schmidt 正交化可写为 $ \vec{Q} = [\vec{q}_{1}, \vec{q}_{2}, \vec{q}_{3}] $，其中：

$$ \begin{aligned}\vec{q}_{1}&=\frac{1}{r_{11}}\vec{a}_{1},\\\vec{q}_{2}&=\frac{1}{r_{22}}\left(\vec{a}_{2}-\frac{r_{12}}{r_{11}}\vec{a}_{1}\right),\\\vec{q}_{3}&=\frac{1}{r_{33}}\left[\vec{a}_{3}-\frac{r_{23}}{r_{22}}\vec{a}_{2}-\left(\frac{r_{13}}{r_{11}}-\frac{r_{12}}{r_{11}}\frac{r_{23}}{r_{22}}\right)\vec{a}_{1}\right].\end{aligned}\tag{5.58} $$

上三角矩阵 $ \bar{R} $ 的各元素为：

$$ \begin{aligned}&r_{11}=\sqrt{\sum_{j=1}^{N_{A}}(\Delta x_{ij})^{2}},\\&r_{12}=\frac{1}{r_{11}}\sum_{j=1}^{N_{A}}\Delta x_{ij}\Delta y_{ij},\\&r_{22}=\sqrt{\sum_{j=1}^{N_{A}}(\Delta y_{ij})^{2}-r_{12}^{2}},\\&r_{13}=\frac{1}{r_{11}}\sum_{j=1}^{N_{A}}\Delta x_{ij}\Delta z_{ij},\\&r_{23}=\frac{1}{r_{22}}\left(\sum_{j=1}^{N_{A}}\Delta y_{ij}\Delta z_{ij}-\frac{r_{12}}{r_{11}}\sum_{j=1}^{N_{A}}\Delta x_{ij}\Delta z_{ij}\right),\\&r_{33}=\sqrt{\sum_{j=1}^{N_{A}}(\Delta z_{ij})^{2}-(r_{13}^{2}+r_{23}^{2})}.\\ \end{aligned}\tag{5.59} $$

由式（5.57）至（5.59），节点 i 处的梯度可写为边差的加权和：

$$ \nabla U_{i}\equiv\vec{x}=\sum_{j=1}^{N_{A}}\vec{w}_{ij}\theta_{j}(U_{j}-U_{i})\tag{5.60} $$

权重矢量 $ \vec{w}_{ij} $ 定义为：

$$ \vec{w}_{ij}=\left[\begin{array}{c}\alpha_{ij,1}-\frac{r_{12}}{r_{11}}\alpha_{ij,2}+\beta\alpha_{ij,3},\\\alpha_{ij,2}-\frac{r_{23}}{r_{22}}\alpha_{ij,3},\\\alpha_{ij,3}\end{array}\right].\tag{5.61} $$

式（5.61）中的各项为：

$$ \begin{aligned}\alpha_{ij,1}&=\frac{\Delta x_{ij}}{r_{11}^{2}},\\\alpha_{ij,2}&=\frac{1}{r_{22}^{2}}\left(\Delta y_{ij}-r_{12}r_{11}\Delta x_{ij}\right),\\\alpha_{ij,3}&=\frac{1}{r_{33}^{2}}\left(\Delta z_{ij}-r_{23}r_{22}\Delta y_{ij}+\beta\Delta x_{ij}\right),\end{aligned}\tag{5.62} $$

校注：原书式（5.62）写出乘积 $r_{12}r_{11}$ 和 $r_{23}r_{22}$。按式（5.58）、（5.59）的 QR 消元关系，这两处应核对为比值 $r_{12}/r_{11}$ 和 $r_{23}/r_{22}$。译文保留原式。

其中，

$$ \beta=\frac{r_{12}r_{23}-r_{13}r_{22}}{r_{11}r_{22}}.\tag{5.63} $$

最小二乘方法用于单元中心格式时，形式不变，只需将节点换成单元形心。实例见文献 $ ^{[16]} $。

最小二乘方法在一般网格上具有一阶精度 $ [64] $，也满足相容性：不论单元类型如何，线性函数梯度的计算误差都仅为舍入误差。因此，它特别适合混合单元网格。其计算量与 Green–Gauss 方法相当，因为在一次面／边循环中，只需进行矢量与标量的乘法（式 $ (5.60) $）。不过，必须在每个节点预计算并存储六个元素（式 $ (5.59) $），它们属于上三角矩阵 $ \bar{R} $。文献 $ [89] $ 的详细研究表明，加权系数 $ \theta_{j} $ 在式 $ (5.55) $ 和式 $ (5.60) $ 中必须取节点 i、j 间距离的倒数，类似权重 $ \theta_{iJ} $ 在式 $ (5.45) $ 中的取值，才能在高度拉伸且同时弯曲的网格上准确近似非线性函数的梯度。不过，对三角形（四面体）网格上的单元中心格式，这一做法并无帮助 $ [89] $。

经验还表明，对中位对偶格式，若黏性壁面上采用棱柱或六面体单元，最小二乘方法需要特别注意。考察图 5.15，假设要计算节点 i 处速度分量的梯度。可见，只有边 ij 的贡献有效，因为其他通过边与 i 相连的节点处均有 $ u = v = w = 0 $。为增加模板的支持，可插入虚拟边 [50]，即图 5.15 中的虚线。虚拟边能明显提高离散格式的精度和稳健性。应强调，它们仅用于梯度重构，不用于通量计算。

5.3.5 限制函数

二阶及更高阶迎风空间离散需要使用所谓限制器或限制函数，避免在大梯度区域（例如激波处）产生振荡和伪解。目标至少是使格式保持单调性，即时间演化过程中，流场最大值不得增大，最小值不得减小，且不能产生新的局部极值。4.3.5 节已针对结构网格迎风离散讨论这一点。

在非结构网格上，限制器用于减小控制体面左右状态重构所用的梯度。强间断处限制函数必须为零，使格式退化为能保证单调性的一阶迎风格式。限制器取零时得到式（5.38）的常量重构。当然，在光滑流动区域应保留原来的无限制重构，以尽量减少数值耗散。下文介绍两种广泛使用的限制函数：Barth–Jespersen 限制器 [30] 和 Venkatakrishnan 限制器 [90, 91]。最近提出的多维限制器概念见文献 [92] 和 [93]。

Barth–Jespersen 限制器

文献 $ ^{[30]} $ 首次介绍了非结构网格上限制函数的实现。对中位对偶格式，节点 i 处定义为：

$$ \Psi_{i}=\min_{j}\left\{\begin{aligned}\min\left(1,\frac{U_{\max}-U_{i}}{\Delta_{2}}\right)\quad&\text{当}\Delta_{2}>0,\\ \min\left(1,\frac{U_{\min}-U_{i}}{\Delta_{2}}\right)\quad&\text{当}\Delta_{2}<0,\\ 1\quad&\text{当}\Delta_{2}=0\end{aligned}\right.\tag{5.64} $$

其中采用以下简记：

$$ \begin{aligned}\Delta_{2}&=\frac{1}{2}\left(\nabla U_{i}\cdot\vec{r}_{ij}\right),\\U_{\max}&=\max\left(U_{i},\max_{j}U_{j}\right),\\U_{\min}&=\min\left(U_{i},\min_{j}U_{j}\right).\end{aligned}\tag{5.65} $$

式（5.64）和式（5.65）中的 $ \min_{j} $ 或 $ \max_{j} $，分别表示节点 i 的全部直接邻居 j，即通过边与 i 相连的节点中的最小值或最大值。此外，图 5.9 或图 5.13(b) 中的边矢量 $ \vec{r}_{ij} $ 按式（5.43）定义，$ U_{j} $ 为某个相邻节点 j 的标量。单元中心格式也有类似公式，只需将节点下标换成单元下标，并取：

$$ \Delta_{2}=\nabla U_{I}\cdot\vec{r}_{\mathrm{L}},\tag{5.66} $$

其中，$ \vec{n}_{L} $ 为从单元形心指向对应单元面中点的矢量。为避免式（5.64）中除以很小的 $ \Delta_{2} $，宜将 $ \Delta_{2} $ 修正为 $ \operatorname{Sign}(\Delta_{2})(|\Delta_{2}| + \omega) $，其中 $ \omega $ 约为机器精度 [90]。

Barth 限制器能保证解的单调性，但耗散较大，往往会抹平间断。另一个问题是，光滑流动区域的数值噪声也会激活限制器，通常导致无法完全收敛至定常状态 $ [39, 90] $。因此，Venkatakrishnan 限制函数得到了更广泛的应用。

Venkatakrishnan 限制器

Venkatakrishnan 限制器 $ [90, 91] $ 因其较好的收敛性质而得到广泛使用。它在顶点 i 处将重构梯度 $ \nabla U $ 乘以下列因子：

$$ \Psi_{i}=\min_{j}\left\{\begin{array}{cl}\displaystyle\frac{1}{\Delta_{2}}\left[\frac{(\Delta_{1,\max}^{2}+\epsilon^{2})\Delta_{2}+2\Delta_{2}^{2}\Delta_{1,\max}}{\Delta_{1,\max}^{2}+2\Delta_{2}^{2}+\Delta_{1,\max}\Delta_{2}+\epsilon^{2}}\right]&\text{当}\Delta_{2}>0,\\\displaystyle\frac{1}{\Delta_{2}}\left[\frac{(\Delta_{1,\min}^{2}+\epsilon^{2})\Delta_{2}+2\Delta_{2}^{2}\Delta_{1,\min}}{\Delta_{1,\min}^{2}+2\Delta_{2}^{2}+\Delta_{1,\min}\Delta_{2}+\epsilon^{2}}\right]&\text{当}\Delta_{2}<0,\\\qquad\qquad\qquad\qquad\qquad1&\text{当}\Delta_{2}=0,\end{array}\right.\tag{5.67} $$

其中，

$$ \begin{aligned}\Delta_{1,\max}&=U_{\max}-U_{i},\\\Delta_{1,\min}&=U_{\min}-U_{i}.\end{aligned}\tag{5.68} $$

这里，$ U_{max} $ 和 $ U_{min} $ 为所有周围节点 j 以及节点 i 自身的最小值／最大值。$ U_{max} $、$ U_{min} $、$ \Delta_{2} $ 的定义见式（5.65）。参数 $ \epsilon^{2} $ 用于控制限制程度。将 $ \epsilon^{2} $ 设为零会得到完全限制，但可能使收敛停滞；反之，若 $ \epsilon^{2} $ 很大，限制函数将接近 1，即几乎不作限制，解中可能出现振荡。实践表明，$ \epsilon^{2} $ 应与局部长度尺度成比例，即：

$$ \epsilon^{2}=(K\Delta h)^{3},\tag{5.69} $$

其中，K 为 $ \mathcal{O}(1) $ 的常数，$ \Delta h $ 例如可取控制体体积的立方根，二维则取面积的平方根。应注意，限制函数（5.67）必须用无量纲量定义。式（5.69）中系数 K 对激波分辨的影响见图 5.16。

可见，完全限制的解 $ (K=0) $ 与 K=5 的解相同。不过，显式时间推进格式在 K=0 时仅收敛约三个数量级，而 K=5 时收敛至机器零（图 5.17）。图 5.16 还显示，随着 K 增大，解逐渐趋于无限制状态，表现为激波处过冲不断增大。

计算上述任一限制函数的工作量都较大。为求 $ U_{max} $、$ U_{min} $ 及限制器 $ \Psi $ 本身，需要两次边循环（单元中心格式为面循环）和一次节点循环（单元中心格式为单元循环）。此外，对每个流动变量，均须按节点或单元分别存储 $ U_{max} $、$ U_{min} $ 和 $ \Psi $。