9.3 残差光顺

9.3 残差光顺

通过优化级系数 $ [5-7] $，可在一定程度上改变显式多级时间推进格式（6.1.1、6.1.2 节）的最大 CFL 数和收敛性质。Jameson 和 Baker $ [8] $ 引入残差光顺技术，旨在使显式格式具有隐式特征，从而显著提高最大允许 CFL 数。残差光顺的另一个目的是增强对残差高频误差分量的阻尼，这对成功应用多重网格方法尤为重要。

残差光顺可以用显式、隐式或混合方式实现 $ [9, 10] $，通常在显式时间推进格式（式（6.5）、（6.7））的每一级应用。先将计算得到的残差 $ \vec{R}^{(k)} $ 替换为光顺残差 $ \vec{R}^{*} $，再更新解 $ \vec{W}^{(k)} $。下面讨论常用的隐式残差光顺（IRS）在结构网格和非结构网格上的实现。

9.3.1 结构网格上的中心 IRS

IRS 的标准三维形式为

$$ \begin{aligned}-\epsilon^{I}\vec{R}_{I-1,J,K}^{*}+(1+2\epsilon^{I})\vec{R}_{I,J,K}^{*}-\epsilon^{I}\vec{R}_{I+1,J,K}^{*}&=\vec{R}_{I,J,K},\\-\epsilon^{J}\vec{R}_{I,J-1,K}^{**}+(1+2\epsilon^{J})\vec{R}_{I,J,K}^{**}-\epsilon^{J}\vec{R}_{I,J+1,K}^{**}&=\vec{R}_{I,J,K}^{*},\\-\epsilon^{K}\vec{R}_{I,J,K-1}^{***}+(1+2\epsilon^{K})\vec{R}_{I,J,K}^{***}-\epsilon^{K}\vec{R}_{I,J,K+1}^{***}&=\vec{R}_{I,J,K}^{**},\end{aligned}\tag{9.5} $$

其中，$ \vec{R}_{I,J,K}^{*}, \vec{R}_{I,J,K}^{**} $ 和 $ \vec{R}_{I,J,K}^{***} $ 分别表示 I、J、K 方向的光顺残差；参数 $ \epsilon^{I}, \epsilon^{J} $ 和 $ \epsilon^{K} $ 表示三个计算坐标方向的光顺系数。式（9.5）的隐式算子类似于二阶中心差分，因此称为中心隐式残差光顺（CIRS）。式（9.5）的隐式方程组采用 Thomas 算法对三对角矩阵求逆得到解。

光顺系数通常定义为对流通量 Jacobian 谱半径的函数 $ [11] $，旨在使各坐标方向仅采用满足稳定性和良好误差阻尼所需的光顺量。文献 $ [12] $ 提出了一种适宜的二维公式：

$$ \begin{aligned}&\epsilon^{I}=\max\left\{\frac{1}{4}\left[\left(\frac{\sigma^{*}}{\sigma}\frac{1}{1+\Psi r}\right)^{2}-1\right],0\right\},\\&\epsilon^{J}=\max\left\{\frac{1}{4}\left[\left(\frac{\sigma^{*}}{\sigma}\frac{1}{1+\Psi/r}\right)^{2}-1\right],0\right\},\\ \end{aligned}\tag{9.6} $$

其中，$ \sigma^{*}/\sigma $ 表示采用和不采用光顺的格式之 CFL 数的比值。变量 r 为对流谱半径之比（式（4.53）），即 $ r = \hat{\Lambda}_{c}^{J}/\hat{\Lambda}_{c}^{I} $。参数 $ \Psi \approx 0.125 $ 保证光顺操作的线性稳定性。

三维光顺系数可用与式（9.6）相似的表达式计算 [12]：

$$ \begin{aligned}\epsilon^{I}&=\max\left\{\frac{1}{4}\left[\left(\frac{\sigma^{*}}{\sigma}\frac{1}{1+\Psi(r^{JI}+r^{KI})}\right)^{2}-1\right],0\right\}\\\epsilon^{J}&=\max\left\{\frac{1}{4}\left[\left(\frac{\sigma^{*}}{\sigma}\frac{1}{1+\Psi(r^{KJ}+r^{IJ})}\right)^{2}-1\right],0\right\}\\\epsilon^{K}&=\max\left\{\frac{1}{4}\left[\left(\frac{\sigma^{*}}{\sigma}\frac{1}{1+\Psi(r^{IK}+r^{JK})}\right)^{2}-1\right],0\right\},\end{aligned}\tag{9.7} $$

其中，

$$ r^{J I}=\hat{\Lambda}_{c}^{J}/\hat{\Lambda}_{c}^{I},\quad r^{K I}=\hat{\Lambda}_{c}^{K}/\hat{\Lambda}_{c}^{I}, $$

其余类似。参数 $ \Psi $ 的典型值为 0.0625。

CFL 数比值 $ \sigma^{*}/\sigma $ 的最大值取决于光顺系数及空间离散格式的类型。对于中心格式（4.3.1 节），该比值为

$$ \frac{\sigma^{*}}{\sigma}\leq\sqrt{1+4\epsilon}.\tag{9.8} $$

实际可达到 $ \sigma^{*}/\sigma\approx2 $（$ \epsilon=0.8 $）。更高比值会削弱时间推进格式的阻尼。对于迎风空间离散，没有像式（9.8）这样简单的条件，因为 $ \sigma^{*}/\sigma $ 的最大值还取决于级系数。不过，通常无黏流动的 CFL 数（见表 6.1、6.2）可提高约两倍，黏性流动采用表 6.2 的混合格式时可提高至五倍。图 9.2、9.3 展示了结构网格上中心 IRS 对无黏及黏性流动收敛的影响。为了计及 IRS 所增加的计算量，比较采用 CPU 时间。可以看出，求解时间能够显著缩短，尤其是

对于黏性流动。将 IRS 与多重网格（9.4 节）结合，还可以节省更多计算成本。

借助 IRS，可以抵消式（6.18）中黏性谱半径 $ \hat{\Lambda}_{ v} $ 对时间步长的限制。最大时间步长按式（6.14）计算，不计 $ \hat{\Lambda}_{ v} $。在黏性谱半径占主导的区域，采用更高的光顺系数 $ \epsilon^{I}, \epsilon^{J}, \epsilon^{K} $。基于黏性谱半径的光顺系数可由以下公式计算 [12,13]：

$$ \begin{aligned}\epsilon_{ v}^{I}&=C\left(\frac{\sigma^{*}}{\sigma}\right)\frac{\hat{\Lambda}_{ v}^{I}}{\hat{\Lambda}_{c}^{I}+\hat{\Lambda}_{c}^{J}+\hat{\Lambda}_{c}^{K}},\\\epsilon_{ v}^{J}&=C\left(\frac{\sigma^{*}}{\sigma}\right)\frac{\hat{\Lambda}_{ v}^{J}}{\hat{\Lambda}_{c}^{I}+\hat{\Lambda}_{c}^{J}+\hat{\Lambda}_{c}^{K}},\\\epsilon_{ v}^{K}&=C\left(\frac{\sigma^{*}}{\sigma}\right)\frac{\hat{\Lambda}_{ v}^{K}}{\hat{\Lambda}_{c}^{I}+\hat{\Lambda}_{c}^{J}+\hat{\Lambda}_{c}^{K}},\end{aligned}\tag{9.9} $$

其中，常数 $ C \approx 5/4 $。最终光顺系数取式（9.7）至（9.9）各系数的最大值。

9.3.2 非结构网格上的中心 IRS

在非结构网格上，将 Laplace 算子（见式（5.24））作用于残差即可实现 CIRS。因此，控制体 I 中的光顺残差 $ \vec{R}_{I}^{*} $ 由以下隐式关系得到 [14]：

$$ \vec{R}_{I}^{*}+\sum_{J=1}^{N_{A}}\epsilon\left(\vec{R}_{I}^{*}-\vec{R}_{J}^{*}\right)=\vec{R}_{I}.\tag{9.10} $$

求和包括全部 $ N_{A} $ 个相邻控制体。采用 Jacobi 迭代求解式（9.10），得到 $ \vec{R}_{I}^{*} $。适宜的光顺系数为 $ 0.3 \leq \epsilon \leq 0.8 $；采用这些 $ \epsilon $ 值，CFL 数及时间步长可提高两至五倍。由于矩阵具有对角占优性，Jacobi 迭代约两步即可收敛。

9.3.3 结构网格上的迎风 IRS

前述 CIRS 方法对亚声速和跨声速流动效果较好，在黏性占主导的区域也有帮助。不过，将 CIRS 与迎风空间离散格式结合时，其误差阻尼特性较差；出现强激波时，经 CIRS 加速的多级格式的稳健性也会下降。因此，人们发展了所谓的迎风隐式残差光顺（UIRS）方法 $ [15, 16] $，特别适用于高马赫数流动。

与 CIRS 不同，UIRS 计及对流特征值 $ \bar{\Lambda}_{c} $ 的正负。其思路是仅沿 Euler 方程特征线方向光顺残差，即 dx/dt = const. = $ \Lambda_{c} $，避免对上游残差产生非物理影响。UIRS 需要将残差变换为特征变量（参见 A.11 节），这样便可根据相应特征值的正负，独立光顺残差矢量的各分量。残差第 l 个分量的一维隐式算子定义为 [15,16]

$$ \begin{array}{r l}{-\epsilon^{I}(R^{*})_{I-1}^{l}+(1+\epsilon^{I})(R^{*})_{I}^{l}=(R^{c})_{I}^{l}}&{\quad\text{当}\quad\Lambda_{c}^{l}\geq0,}\\ {(1+\epsilon^{I})(R^{*})_{I}^{l}-\epsilon^{I}(R^{*})_{I+1}^{l}=(R^{c})_{I}^{l}}&{\quad\text{当}\quad\Lambda_{c}^{l}<0,}\end{array}\tag{9.11} $$

其中，$ \vec{R}^{c} $ 表示变换为特征变量的残差，即

$$ \vec{R}^{c}=\bar{T}^{-1}\vec{R}.\tag{9.12} $$

采用 Thomas 算法求解三对角方程组，得到光顺残差 $ \vec{R}^{*} $；若 $ \Lambda_{c}^{I} $ 不变号，则为双对角方程组。随后，将残差 $ \vec{R}^{*} $ 变换回守恒变量，即可更新解。

UIRS 最理想的特性是使多级格式具有很好的阻尼性能，尤其与迎风空间离散结合时。即使光顺系数很大，解误差的阻尼能力仍然保持，甚至得到改善。对于一维 Euler 方程，已经证明 $ \epsilon = 500 $、$ \sigma^{*} = 1000 $ 等参数值能够产生稳定且非常快速的显式多级格式 [15,16]，亦见源程序的 analysis/mstage。然而，多维 UIRS 的实现存在困难。由于对流通量 Jacobian 不能沿全部坐标方向同时对角化（见 A.11 节），必须针对每个计算坐标分别执行式（9.12）的变换及式（9.11）的光顺。坐标分裂使最大光顺系数降低至 $ 2 \leq \epsilon \leq 6 $，但与 CIRS 相比，定常收敛仍明显加快 [16,17]。与多重网格结合时，收敛速度和稳健性的改善最大 [16,18]。

在多维情况下，光顺系数按特征值缩放。例如，二维时可以采用以下公式：

$$ \begin{aligned}\epsilon^{I}&=\epsilon\cdot\min\left[\frac{\Lambda_{c}^{I}}{\Lambda_{c}^{J}},1\right],\\\epsilon^{J}&=\epsilon\cdot\min\left[\frac{\Lambda_{c}^{J}}{\Lambda_{c}^{I}},1\right].\end{aligned}\tag{9.13} $$

采用光顺的格式之 CFL 数与系数 $ \epsilon $ 之间的关系为

$$ \frac{\sigma^{*}}{\sigma}\leq1+C\epsilon.\tag{9.14} $$

常数 C 取决于空间离散类型。中心格式取 C = 1，一阶或二阶迎风格式取 C = 2。

为了避免变换为特征变量的计算成本，人们提出了简化 UIRS 方法 $ [16–18] $。其 I 方向表达式为

$$ \begin{aligned}&-\epsilon^{I}\vec{R}_{I-1}^{*}+(1+\epsilon^{I})\vec{R}_{I}^{*}=\vec{R}_{I}\quad&\text{当}\quad M^{I}>1,\\&-\epsilon^{I}\vec{R}_{I-1}^{*}+(1+2\epsilon^{I})\vec{R}_{I}^{*}-\epsilon^{I}\vec{R}_{I+1}^{*}=\vec{R}_{I}\quad&\text{当}\quad|M^{I}|\leq1,\\&(1+\epsilon^{I})\vec{R}_{I}^{*}-\epsilon^{I}\vec{R}_{I+1}^{*}=\vec{R}_{I}\quad&\text{当}\quad M^{I}<-1.\end{aligned}\tag{9.15} $$

马赫数 M 根据投影到相应计算坐标方向（此处为 I 方向）的速度确定。由于运算次数少，简化 UIRS 尤其适合三维流动。文献 [16] 表明，与 CIRS 相比，达到定常状态所需的 CPU 时间可减少一半。图 9.4 给出另一示例，可以看出 UIRS 有助于稳定多重网格格式并加速收敛。采用 CIRS 的格式无法得到完全收敛的解，若不采用迎风延拓，甚至无法正常工作。