6.2 隐式时间推进格式

6.2 隐式时间推进格式

在式（6.2）中取 $ \beta\neq0 $，可得到各种隐式时间积分格式。对定常流动问题，取 $ \omega=0 $ 的隐式格式最适合；非定常流动见 6.3 节。于是，式（6.2）简化为：

$$ \frac{(\Omega\bar{M})_{I}}{\Delta t_{I}}\Delta\vec{W}_{I}^{n}=-\beta\vec{R}_{I}^{n+1}-\left(1-\beta\right)\vec{R}_{I}^{n}.\tag{6.26} $$

可见，隐式形式产生了一组关于时刻 $ (t+\Delta t) $ 未知流动变量的非线性方程。求解式（6.26）需要新时间层的残差 $ \vec{R}^{n+1} $；由于 $ \vec{W}^{n+1} $ 未知，无法直接计算。不过，可在当前时间层附近，将式（6.26）的残差 $ \vec{R}^{n+1} $ 线性化，即：

$$ \vec{R}_{I}^{n+1}\approx\vec{R}_{I}^{n}+\left(\frac{\partial\vec{R}}{\partial\vec{W}}\right)_{I}\Delta\vec{W}^{n},\tag{6.27} $$

其中，$ \partial\vec{R}/\partial\vec{W} $ 称为通量 Jacobian。应指出，它往往根据空间离散 $ \vec{R}^{n} $ 的粗略近似导出。例如，高阶迎风离散通常只用一阶迎风格式构造通量 Jacobian。不过，为获得最佳效率和稳健性，它仍须体现空间离散的主要特征。

将式（6.27）的线性化代替式（6.26）中的 $ \vec{R}^{n+1} $，得到以下隐式格式：

$$ \left[\frac{(\Omega\bar{M})_{I}}{\Delta t_{I}}+\beta\left(\frac{\partial\vec{R}}{\partial\vec{W}}\right)_{I}\right]\Delta\vec{W}^{n}=-\vec{R}_{I}^{n}.\tag{6.28} $$

式（6.28）左端方括号内的项称为隐式算子或系统矩阵，右端相应称为显式算子。解的空间精度只由显式算子决定。

隐式算子构成大型稀疏非对称分块矩阵，其分块维数等于总单元数（单元中心格式）或总网格点数（单元顶点格式）。以下分别讨论结构与非结构网格上的形式。如 3.2 节所述，质量矩阵 $ \bar{M} $ 可换为单位矩阵，而不影响定常解。式（6.28）的参数 $ \beta $ 通常取 1，得到一阶时间精度；取 $ \beta = 1/2 $ 可得二阶时间精度，但不推荐这样做，因为 $ \beta = 1 $ 的格式稳健得多，且时间精度对定常问题并不重要。

控制方程具有刚性时，源项必须包含在隐式算子中。式（6.27）的残差线性化自然会纳入源项，产生与式（6.10）相同的形式。文献 $ ^{[12]} $ 证明，只要 $ \partial\vec{Q}/\partial\vec{W} $ 的全部特征值均非正，上述隐式格式（6.28）对任意时间步长都稳定。

求解线性方程组 $ (6.28) $，需要对隐式算子求逆，即对一个极大型矩阵求逆。原则上有两种方式。第一种直接求逆，采用 Gauss 消元或某种直接稀疏矩阵方法 $ [21, 22] $。不过，由于内存和计算量过大，这种方式不适用于实际问题 $ [23] $。

第二种方式采用迭代方法，3.2.2 节已提到最常用的方法。大致可分两类。第一类将隐式算子分成若干部分，称为因子分解；各因子的构造使其比原算子更易求逆。第二类用 Krylov 子空间方法对隐式算子求逆。此时，通常令 $ \Delta t \rightarrow \infty $，将隐式时间推进格式（6.28）变成 Newton 方法，称为 Newton–Krylov 方法。

下文先讨论隐式算子的矩阵结构，再考察式（6.28）中通量 Jacobian $ \partial\vec{R}/\partial\vec{W} $ 的计算，随后详细介绍三种最常用的迭代方法，最后讨论隐式 Runge–Kutta 格式。

6.2.1 隐式算子的矩阵形式

根据式（6.27），残差 $ \vec{R}^{n+1} $ 的线性化可写为：

$$ \vec{R}^{n+1}\approx\vec{R}^{n}+\sum_{m=1}^{N_{F}}\left\{\frac{\partial}{\partial\vec{W}}\left[(\vec{F}_{\mathrm{c}}-\vec{F}_{\mathrm{v}})_{m}\Delta S_{m}\right]\Delta\vec{W}^{n}\right\}-\frac{\partial(\Omega\vec{Q})}{\partial\vec{W}}\Delta\vec{W}^{n}\tag{6.29} $$

其中，$ N_{F} $ 为控制体 $ \Omega $ 的面数，参见式（4.2）或式（5.2）。于是，通量 Jacobian 为：

$$ \frac{\partial\vec{R}}{\partial\vec{W}}=\sum_{m=1}^{N_{F}}\frac{\partial(\vec{F}_{c})_{m}}{\partial\vec{W}}\Delta S_{m}-\sum_{m=1}^{N_{F}}\frac{\partial(\vec{F}_{v})_{m}}{\partial\vec{W}}\Delta S_{m}-\frac{\partial(\Omega\vec{Q})}{\partial\vec{W}}.\tag{6.30} $$

必须强调，通量 Jacobian 应视为作用于增量 $ \Delta\vec{W} $ 的算子。如上所述，式（6.30）的对流通量与黏性通量未必必须与显式算子中的通量相同。

由于结构与非结构网格上的隐式算子明显不同，下面分别讨论。

结构网格上的隐式算子

为导出式（6.28）的系统矩阵，考察图 6.1 的一维网格。假设空间离散采用具有对偶控制体的单元顶点格式（4.2.3 节）及简单通量平均，参见图 4.8。不含黏性通量时，残差为：

$$ \vec{R}_{i}^{n}=(\vec{F}_{c})_{i+1/2}\Delta S_{i+1/2}+(\vec{F}_{c})_{i-1/2}\Delta S_{i-1/2}-\Omega_{i}\vec{Q}_{i}.\tag{6.31} $$

此外，式（6.30）中面 $ m = i + \frac{1}{2} $ 的对流通量导数可写为：

$$ \begin{aligned}\frac{\partial(\vec{F}_{c})_{i+1/2}}{\partial\vec{W}}&=\frac{\partial}{\partial\vec{W}}\left\{\frac{1}{2}\left[\vec{F}_{c}(\vec{W}_{i+1}^{n})+\vec{F}_{c}(\vec{W}_{i}^{n})\right]\right\}\\&=\frac{1}{2}\left[(\bar{A}_{c})_{i+1}+(\bar{A}_{c})_{i}\right],\end{aligned}\tag{6.32} $$

其中，$ \bar{A}_{c} $ 为对流通量 Jacobian（A.9 节）。由式（6.29）、（6.31）、（6.32），时刻 $ (t+\Delta t) $ 的残差近似为：

$$ \begin{aligned}\vec{R}_{i}^{n+1}\approx\vec{R}_{i}^{n}+&\frac{1}{2}\left[(\bar{A}_{c})_{i+1}\Delta\vec{W}_{i+1}^{n}+(\bar{A}_{c})_{i}\Delta\vec{W}_{i}^{n}\right]\Delta S_{i+1/2}\\&+\frac{1}{2}\big[(\bar{A}_{c})_{i-1}\Delta\vec{W}_{i-1}^{n}+(\bar{A}_{c})_{i}\Delta\vec{W}_{i}^{n}\big]\Delta S_{i-1/2}-\frac{\partial(\Omega_{i}\vec{Q}_{i})}{\partial\vec{W}}.\end{aligned}\tag{6.33} $$

最后，取 $ \beta = 1 $ 并采用集中质量矩阵，可由式（6.26）导出隐式格式：

$$ \begin{aligned}\Bigl\{\frac{\Omega_{i}}{\Delta t_{i}}\bar{I}+&\frac{1}{2}\left[(\bar{A}_{c})_{i}\Delta S_{i+1/2}+(\bar{A}_{c})_{i}\Delta S_{i-1/2}\right]-\frac{\partial(\Omega_{i}\vec{Q}_{i})}{\partial\vec{W}}\\&+\frac{1}{2}\left[(\bar{A}_{c})_{i+1}\Delta S_{i+1/2}\right]+\frac{1}{2}\left[(\bar{A}_{c})_{i-1}\Delta S_{i-1/2}\right]\Bigr\}\Delta\vec{W}^{n}=-\vec{R}_{i}^{n},\end{aligned}\tag{6.34} $$

其中，$ \bar{I} $ 为单位矩阵。应注意，每个矩阵 $ \bar{A}_{c} $ 的单位法向量，与对应面积 $ \Delta S $，均在控制体的同一侧计算。可见，隐式算子与空间离散采用相同的三点模板，即节点 i−1、i、$ i + 1 $。为直观展示系统矩阵，将隐式算子中与中心节点 i 相关的全部项记为 D，即：

$$ \mathbf{D}\equiv\frac{\Omega_{i}}{\Delta t_{i}}\bar{I}+\frac{1}{2}\left[(\bar{A}_{c})_{i}\Delta S_{i+1/2}+(\bar{A}_{c})_{i}\Delta S_{i-1/2}\right]-\frac{\partial(\Omega_{i}\vec{Q}_{i})}{\partial\vec{W}}, $$

将与下风节点 $ (i+1) $ 相关的项记为 U，即：

$$ \mathbf{U}\equiv\frac{1}{2}\left(\bar{A}_{\mathrm{c}}\right)_{i+1}\Delta S_{i+1/2}, $$

将与上风节点 $ (i-1) $ 相关的项记为 L，即：

$$ \mathbf{L}\equiv\frac{1}{2}\left(\bar{A}_{\mathrm{c}}\right)_{i-1}\Delta S_{i-1/2}, $$

对图 6.1 网格的全部八个节点写出式（6.34），便得到图右侧的 $ 8 \times 8 $ 分块三对角矩阵。由于一维有三个守恒方程，每个 L、D、U 块均为 $ 3 \times 3 $ 矩阵。

同样的思路适用于多维。例如，若二维空间离散采用图 6.2 的五点模板，则得到

分块五对角矩阵，见图 6.2 右侧。节点排序使 i 下标比 j 下标变化更快，对应于 FORTRAN 中的 $ \text{mat}(i,j) $。注意，第二条副对角线与主对角线相隔八个元素，即 i 方向的节点总数。三维若采用图 4.9(b) 的七点空间模板，系统矩阵则为分块七对角矩阵。

因此，结构网格的系统矩阵始终具有规则、稀疏的带状形式。还应指出，式（6.28）中的以下项可能造成困难：

$$ \frac{(\Omega\bar{M})_{I}}{\Delta t_{I}}\tag{6.35} $$

这一项始终位于主对角线上。时间步长较大时，隐式算子的对角占优性减弱，某些迭代矩阵求逆格式，例如 Gauss–Seidel 方法，可能失效。

非结构网格上的隐式算子

从结构网格转到非结构网格，系统矩阵的外观完全改变。用图 6.3 的小型非结构网格说明。假设空间离散采用 5.2.1 节的单元中心格式，只使用最近邻数值，例如一阶迎风格式，见 5.3.2、5.3.3 节。于是，单元 2 的模板包括单元 18、10、13。所得系统矩阵见图 6.3 右侧。显然，

由于网格单元的编号通常是任意的，中位对偶格式的节点编号也是如此，矩阵不会具有规则模式。只有主对角线始终存在，至少包含 $ (6.35) $，还可能包含源项的导数。

系统矩阵非零元素近似随机的分布并不理想。它会减慢 Gauss–Seidel 等迭代求逆方法的收敛，也使 Krylov 子空间方法的不完全 LU 分解（ILU）等预处理技术难以有效发挥作用。因此，已有方法通过重编号单元或节点，大幅缩小系统矩阵的带宽，使非零元素聚集在主对角线附近。最著名的是反向 Cuthill–McKee（RCM）算法 $ [24, 25] $。将其用于图 6.3 的网格，所得单元编号与系统矩阵见图 6.4。可见，矩阵带宽明显缩小，结构也更规则。

另有一些重编号策略，以减少缓存未命中或便于数值格式向量化为目标。综述例如见文献 $ ^{[26]} $。

6.2.2 通量 Jacobian 的计算

根据所用空间离散格式的类型，式（6.28）中通量 Jacobian $ \partial\vec{R}/\partial\vec{W} $ 的解析计算可能非常复杂，甚至

无法完成。为清楚说明概念，先推导无黏流动的通量 Jacobian，再讨论向 Navier–Stokes 方程的扩展。

中心格式

中心空间离散的通量 Jacobian 最容易构造。如图 6.1 的实例所示，$ \partial\vec{R}/\partial\vec{W} $ 由对流通量 Jacobian 组成（参见式（6.34）），后者可解析导出，见 A.9 节。人工黏性通常以简化形式纳入，不含非线性压力传感器（式（4.55））。6.2.3 节将再讨论这一点。

通量矢量分裂格式

以某种通量矢量分裂格式（4.3.2 节）为基础时，通量 Jacobian 的计算更复杂。以 Steger–Warming 格式 $ [27] $ 为例。已有研究 $ [28] $ 表明，对于隐式算子的各种迎风离散，Steger–Warming 分裂优于 van Leer 通量矢量分裂（式（4.60））等方法。

Steger–Warming 通量矢量分裂的基本思想，是将对流通量分成正、负两部分，即：

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

通量定义为：

$$ \vec{F}_{c}^{\pm}=\bar{A}_{S W}^{\pm}\vec{W}=\left(\bar{T}\bar{\Lambda}^{\pm}\bar{T}^{-1}\right)\vec{W}.\tag{6.37} $$

其中，$ \bar{A}_{SW}^{\pm} $ 为正／负 Steger–Warming 通量分裂 Jacobian，$ \bar{T} $ 为右特征矢量矩阵，$ \bar{T}^{-1} $ 为左特征矢量矩阵，$ \bar{\Lambda}^{\pm} $ 为正／负特征值组成的对角矩阵，参见 A.11 节。特征值矩阵定义为 [27]：

$$ \bar{\Lambda}^{\pm}=\frac{1}{2}(\bar{\Lambda}_{\mathrm{c}}\pm|\bar{\Lambda}_{\mathrm{c}}|),\tag{6.38} $$

其中，$ \bar{\Lambda}_{c} $ 由式（A.84）给出。

采用式（6.36）的分裂，式（6.29）中通量 Jacobian 与增量 $ \Delta\vec{W}^{n} $ 的乘积为：

$$ \frac{\partial\vec{R}_{I}}{\partial\vec{W}}\Delta\vec{W}^{n}=\sum_{m=1}^{N_{F}}\left[\frac{\partial(\vec{F}_{c}^{+}\Delta S)_{m}}{\partial\vec{W}_{\mathrm{L},m}}\Delta\vec{W}_{\mathrm{L},m}^{n}+\frac{\partial(\vec{F}_{c}^{-}\Delta S)_{m}}{\partial\vec{W}_{\mathrm{R},m}}\Delta\vec{W}_{\mathrm{R},m}^{n}\right].\tag{6.39} $$

其中，$ \Delta\vec{W}_{L,m}^{n} $ 与 $ \Delta\vec{W}_{R,m}^{n} $ 分别为面 m 左右状态的增量。结构网格上，可用 MUSCL 方法（式（4.46））计算左右状态；非结构网格则可采用

5.3.3 节介绍的重构方法。不过，精度提高会使模板变宽，增大系统矩阵的带宽。因此，为降低数值复杂性，式（6.39）通常只采用一阶精度近似。一种折中，是以更高精度重构左右状态，但仍采用一阶格式的模板计算导数 [29]。

继续讨论式（6.39）中导数 $ \partial\vec{F}_{c}^{\pm}/\partial\vec{W} $ 的计算。以面 m 上的正通量为例：

$$ \frac{\partial(\vec{F}_{c}^{+}\Delta S)_{m}}{\partial\vec{W}_{\mathrm{L},m}}=\frac{\partial}{\partial\vec{W}_{\mathrm{L},m}}\left[(\bar{A}_{\mathrm{S W}}^{+}\vec{W})_{\mathrm{L},m}\Delta S_{m}\right] $$

由式（6.37）、（6.38），可得：

$$ \begin{aligned}\frac{\partial(\vec{F}_{\mathrm{c}}^{+}\Delta S)_{m}}{\partial\vec{W}_{\mathrm{L},m}}&=\frac{\Delta S_{m}}{2}\left[(\bar{A}_{c})_{\mathrm{L},m}+|(\bar{A}_{c})_{\mathrm{L},m}|\right]\\&\quad+\frac{\Delta S_{m}}{2}\left[\frac{\partial(\bar{A}_{c})_{\mathrm{L},m}}{\partial\vec{W}_{\mathrm{L},m}}+\frac{\partial|(\bar{A}_{c})_{\mathrm{L},m}|}{\partial\vec{W}_{\mathrm{L},m}}\right]\vec{W}_{\mathrm{L},m}^{n}.\end{aligned}\tag{6.40} $$

负通量也有类似式（6.40）的表达式。可见，第一项由对流通量 Jacobian 组成（见 A.9 节），不难计算；第二项却涉及矩阵元素的导数。虽可手工或用符号代数软件获得解析导数，但生成的代码很长，计算效率低 [30]。也可假设矩阵 $ \tilde{A}_{c} $ 局部为常量，从而忽略式（6.40）的第二项。不过，根据隐式格式类型的不同，这可能严重限制 CFL 数 [31]。

计算式（6.39）的导数 $ \partial\vec{F}_{c}^{\pm}/\partial\vec{W} $，还可采用源代码自动微分，例如 ADIFOR [32]、ADIC [33]、OpenAD [34]、Adept [35]，或有限差分方法，例如见文献 [23, 30]。矢量 $ \vec{F} $ 的第 i 分量相对于依赖变量 $ \vec{X} $ 第 j 分量的导数，可近似为：

$$ \frac{\partial f_{i}}{\partial x_{j}}\approx\frac{f_{i}(\vec{X}+h_{j}\vec{e}^{j})-f_{i}(\vec{X})}{h_{j}},\tag{6.41} $$

其中，$ \vec{e}^{j} $ 为第 j 个标准基矢量。Dennis 和 Schnabel [36] 建议步长 $ h_{j} $ 取：

$$ h_{j}=\sqrt{\epsilon}\max\left\{\left|x_{j}\right|,\operatorname{typ}x_{j}\right\}\operatorname{sign}(x_{j})\tag{6.42} $$

其中，$ \epsilon $ 为机器精度，typ $ x_{j} $ 表示 $ x_{j} $ 的典型量级。

高效数值计算 Jacobian 矩阵的建议亦见文献 [37, 38]。

通量差分分裂格式

对 Roe 通量差分分裂格式（4.3.3 节，式（4.91）），式（6.29）中通量 Jacobian 与增量的乘积可写为：

$$ \begin{aligned}\frac{\partial\vec{R}_{I}}{\partial\vec{W}}\Delta\vec{W}^{n}=&\sum_{m=1}^{N_{F}}\frac{\Delta S_{m}}{2}\left\{(\bar{A}_{\mathrm{c}})_{\mathrm{L},m}\Delta\vec{W}_{\mathrm{L},m}^{n}+(\bar{A}_{\mathrm{c}})_{\mathrm{R},m}\Delta\vec{W}_{\mathrm{R},m}^{n}\right.\\&\left.-\frac{\partial}{\partial\vec{W}_{\mathrm{L},m}}\left[|\bar{A}_{\mathrm{Roe}}|_{m}\left(\vec{W}_{\mathrm{R},m}^{n}-\vec{W}_{\mathrm{L},m}^{n}\right)\right]\Delta\vec{W}_{\mathrm{L},m}^{n}\right.\\&\left.-\frac{\partial}{\partial\vec{W}_{\mathrm{R},m}}\left[|\bar{A}_{\mathrm{Roe}}|_{m}\left(\vec{W}_{\mathrm{R},m}^{n}-\vec{W}_{\mathrm{L},m}^{n}\right)\right]\Delta\vec{W}_{\mathrm{R},m}^{n}\right\}.\end{aligned}\tag{6.43} $$

与通量矢量分裂相似，式（6.43）包含对流通量 Jacobian，以及 Roe 矩阵 $ \bar{A}_{Roe} $ 的导数。文献 [31] 给出了这些导数，但相应公式很复杂，亦见文献 [30]，因此按上述方式数值计算 $ \partial\vec{F}_{c}/\partial\vec{W} $ 更合适。也可假设 Roe 矩阵局部为常量，以简化式（6.43）[39]。

$$ \begin{aligned}\frac{\partial\vec{R}_{I}}{\partial\vec{W}}\Delta\vec{W}^{n}&\approx\sum_{m=1}^{N_{F}}\frac{\Delta S_{m}}{2}\left\{(\bar{A}_{\mathrm{c}})_{\mathrm{L},m}\Delta\vec{W}_{\mathrm{L},m}^{n}+(\bar{A}_{\mathrm{c}})_{\mathrm{R},m}\Delta\vec{W}_{\mathrm{R},m}^{n}\right.\\&\left.-|\bar{A}_{\mathrm{Roe}}|_{m}\left(\Delta\vec{W}_{\mathrm{R},m}^{n}-\Delta\vec{W}_{\mathrm{L},m}^{n}\right)\right\}.\end{aligned}\tag{6.44} $$

与 Steger–Warming 通量矢量分裂不同，上述近似线性化 $ (6.44) $ 只会轻微降低隐式格式的性能 $ [31] $。

黏性流动

对 Navier–Stokes 方程，隐式算子还须考虑黏性通量。式（6.30）的导数 $ \partial\vec{F}_{v}/\partial\vec{W} $，即黏性通量 Jacobian，通常不容易得到；黏性通量矢量包含流动变量的导数，使问题进一步复杂。因此，要么用有限差分（式（6.41））计算，要么采用简化形式。

对 Navier–Stokes 方程的薄剪切层（TSL）近似（参见 2.4.3 节和 A.6 节），若假设动力黏度及热导率局部为常量，即可解析求得黏性通量 Jacobian。根据 A.10 节，式（6.29）的黏性通量项成为：

$$ \frac{\partial(\vec{F}_{\mathrm{v}}\Delta S)_{m}}{\partial\vec{W}}\Delta\vec{W}^{n}\approx\left[(\bar{A}_{\mathrm{v}}^{*})_{\mathrm{R},m}\Delta\vec{W}_{\mathrm{R},m}^{n}-(\bar{A}_{\mathrm{v}}^{*})_{\mathrm{L},m}\Delta\vec{W}_{\mathrm{L},m}^{n}\right]\Delta S_{m}.\tag{6.45} $$

其中，$ \bar{A}_{v}^{*} $ 为式（A.71）或式（A.75）的黏性通量 Jacobian，但不含空间算子 $ \partial_{\psi}() $，参见式（A.74）、（A.79）。除动力黏度采用算术平均外，其他变量均采用左状态或右状态计算 Jacobian。若左右状态按一阶精度计算，式（6.45）会在隐式算子中产生二阶中心差分近似。

6.2.3 交替方向隐式格式

交替方向隐式（ADI）格式是最早的迭代隐式方法之一 $ [40] $，只能用于结构网格。它将式（6.28）的隐式算子近似分裂，即因子分解，为二维的两个因子或三维的三个因子。每个因子包含计算空间某一方向上对流与黏性通量的线性化。三维形式为 $ [40, 41] $：

$$ \begin{aligned}\left\{\frac{\Omega}{\Delta t}\bar{I}+\frac{\partial\left[(\vec{F}_{\mathrm{c}}^{I}-\vec{F}_{\mathrm{v}}^{I})\Delta S^{I}\right]_{i+1/2}}{\partial\vec{W}}+\frac{\partial\left[(\vec{F}_{\mathrm{c}}^{I}-\vec{F}_{\mathrm{v}}^{I})\Delta S^{I}\right]_{i-1/2}}{\partial\vec{W}}\right\}\\\left\{\frac{\Omega}{\Delta t}\bar{I}+\frac{\partial\left[(\vec{F}_{\mathrm{c}}^{J}-\vec{F}_{\mathrm{v}}^{J})\Delta S^{J}\right]_{j+1/2}}{\partial\vec{W}}+\frac{\partial\left[(\vec{F}_{\mathrm{c}}^{J}-\vec{F}_{\mathrm{v}}^{J})\Delta S^{J}\right]_{j-1/2}}{\partial\vec{W}}\right\}\\\left\{\frac{\Omega}{\Delta t}\bar{I}+\frac{\partial\left[(\vec{F}_{\mathrm{c}}^{K}-\vec{F}_{\mathrm{v}}^{K})\Delta S^{K}\right]_{k+1/2}}{\partial\vec{W}}+\frac{\partial\left[(\vec{F}_{\mathrm{c}}^{K}-\vec{F}_{\mathrm{v}}^{K})\Delta S^{K}\right]_{k-1/2}}{\partial\vec{W}}\right.\\ \left.\quad-\frac{\partial(\Omega\vec{Q})}{\partial\vec{W}}\right\}\Delta\vec{W}^{n}=-\vec{R}^{n}.\end{aligned}\tag{6.46} $$

校注：原书各方向因子均含 $\Omega/\Delta t$，但未显示因子分解所需的归一化比例。三个因子直接相乘与式（6.28）的算子并不具有相同尺度；使用时应核对归一化约定。译文保留原式。

为清晰起见，式（6.46）在不需要处略去节点下标 i、j、k。上标 I、J、K 表示通量矢量或面面积对应于计算空间中的哪一坐标。例如，$ \Delta S_{i+1/2,j,k}^{I} $ 对应于图 4.1(b) 的 $ \Delta S_{2} $。用单元下标替换节点下标，即可将式（6.46）用于单元中心离散。

式（6.46）的对流与黏性通量导数，可按前一小节的方法计算。ADI 方法传统上与采用人工耗散的中心格式结合。此时，每个因子都具有类似式（6.34）的形式；源项的线性化只包含在一个因子中。为使格式稳健高效，必须在隐式算子中纳入人工耗散项的线性化 $ [42-44] $。通常，将谱半径和耗散系数视为与 $ \vec{W} $ 无关，以简化线性化。于是，根据式（4.48），I 方向的因子例如为：

$$ \left\{\begin{aligned}\frac{\Omega}{\Delta t}\bar{I}+\frac{\partial\left[(\vec{F}_{\mathrm{c}}^{I}-\vec{F}_{\mathrm{v}}^{I})\Delta S^{I}\right]_{i+1/2}}{\partial\vec{W}}+\frac{\partial\left[(\vec{F}_{\mathrm{c}}^{I}-\vec{F}_{\mathrm{v}}^{I})\Delta S^{I}\right]_{i-1/2}}{\partial\vec{W}}\\ -\frac{\partial\left(\vec{D}_{IM}^{I}\Delta S^{I}\right)_{i+1/2}}{\partial\vec{W}}-\frac{\partial\left(\vec{D}_{IM}^{I}\Delta S^{I}\right)_{i-1/2}}{\partial\vec{W}}\end{aligned}\right\} $$

其中，JST 格式的隐式人工耗散项为，参见式（4.50）：

$$ \begin{aligned}&\frac{\partial(\vec{D}_{IM}^{I})_{i+1/2}}{\partial\vec{W}}\Delta\vec{W}^{n}\approx\hat{\Lambda}_{i+1/2}^{S}\left[(\epsilon_{IM}^{(2)})_{i+1/2}\left(\Delta\vec{W}_{i+1}^{n}-\Delta\vec{W}_{i}^{n}\right)\right.\\&\left.-(\epsilon_{IM}^{(4)})_{i+1/2}\left(\Delta\vec{W}_{i+2}^{n}-3\Delta\vec{W}_{i+1}^{n}+3\Delta\vec{W}_{i}^{n}-\Delta\vec{W}_{i-1}^{n}\right)\right].\\ \end{aligned}\tag{6.47} $$

其他方向的隐式耗散项定义类似。隐式算子也可仅保留二阶差分 $ [43, 44] $，使各因子的矩阵从分块五对角变为分块三对角，见图 6.1。不过，文献 $ [44] $ 指出，这会限制格式稳定性。

式（6.46）隐式算子的求逆分三步进行，二维为两步，即：

$$ \begin{array}{l}{\left(\mathbf{D}^{I}+\mathbf{L}^{I}+\mathbf{U}^{I}\right)\Delta\vec{W}^{(1)}=-\vec{R}^{n},}\\ {\left(\mathbf{D}^{J}+\mathbf{L}^{J}+\mathbf{U}^{J}\right)\Delta\vec{W}^{(2)}=\Delta\vec{W}^{(1)},}\\ {\left(\mathbf{D}^{K}+\mathbf{L}^{K}+\mathbf{U}^{K}\right)\Delta\vec{W}^{n}=\Delta\vec{W}^{(2)},}\end{array}\tag{6.48} $$

其中，D 为对角项，L 为下对角项，U 为上对角项。每一步都需直接求解分块三对角矩阵；若式（6.47）中 $ \epsilon_{IM}^{(4)} > 0 $，则为分块五对角矩阵。为降低计算量，Pulliam 和 Chaussee [45] 提出对角化 ADI 格式，将由对流与黏性通量 Jacobian 组成的矩阵块变为对角矩阵。因此，只需对非分块三对角或五对角矩阵求逆，可大幅节省计算量和内存 [43]。对角化虽严格而言只对 Euler 方程成立，但也可用于黏性流动 [43]。此时，从隐式算子中舍去黏性通量的线性化，例如式（6.45），或以黏性特征值（式（6.19））近似。

隐式算子的分裂引入所谓因子分解误差，即式（6.28）基本格式的隐式算子与分解后算子的差。在 ADI 格式中，该误差项由 $ (\Delta t)^{N} $ 缩放，N 为空间维数。这一项使三维 ADI 格式失去无条件稳定性 [46]；不过，若在隐式算子中包含人工黏性的四阶差分，稳定性会改善 [44]。事实上，ADI 格式已成功用于求解各种三维问题 $ [43, 47] $。Rosenfeld 和 Yassour $ [48] $ 提出一种多块网格上的有趣实现，以隐式方式处理块边界。

时间步长 $ \Delta t $ 可按 6.1.4 节的方法，用式（6.14）计算。最佳 CFL 数在 20 至 50 之间，取决于具体流动的物理性质。

6.2.4 LU 对称 Gauss–Seidel 格式

隐式 LU 对称 Gauss–Seidel（LU-SGS）格式，也称 LU 对称逐次超松弛（LU-SSOR）格式，因数值复杂性低、内存需求小而广泛应用，两者均可与显式多级格式相当。此外，它易于在向量及并行计算机上实现，也同时适用于结构和非结构网格。

LU-SGS 的起源可追溯至 Jameson 和 Turkel 的工作 $ [49] $，他们将隐式算子分解成对角占优的下、上三角因子。Yoon 和 Jameson $ [50–52] $ 首次提出 LU-SGS 方法本身，将其作为求解式（6.28）未分解隐式格式的松弛方法。Rieger 和 Jameson $ [53] $ 进一步发展并用于三维黏性流场。此后，研究者将其用于结构网格 $ [31, 54, 56–63] $ 及非结构网格 $ [64–71] $ 上的黏性流动，也常用于化学反应流模拟 $ [72–78] $。

LU-SGS 格式采用 Steger–Warming 一阶精度通量矢量分裂的简化形式（见 6.2.2 节），线性化式（6.29）的对流通量。无论显式算子如何离散，该线性化形式都相同。此外，它将式（6.28）的隐式算子分解为以下三部分：

$$ \left(\mathbf{D}+\mathbf{L}\right)\mathbf{D}^{-1}\left(\mathbf{D}+\mathbf{U}\right)\Delta\vec{W}^{n}=-\vec{R}_{I}^{n}.\tag{6.49} $$

构造因子时，使 L 只包含严格下三角部分，U 只包含严格上三角部分，D 包含对角项。应注意，因子数始终相同，与空间维数无关。

式（6.49）的 LU-SGS 系统矩阵可分两步求逆，即前向和后向扫掠：

$$ \begin{array}{l}{(\mathbf{D}+\mathbf{L})\Delta\vec{W}^{(1)}=-\vec{R}_{I}^{n},}\\ {(\mathbf{D}+\mathbf{U})\Delta\vec{W}^{n}=\mathbf{D}\Delta\vec{W}_{I}^{(1)}}\end{array}\tag{6.50} $$

其中，$ \vec{W}^{n+1} = \vec{W}^{n} + \Delta\vec{W}^{n} $。结构与非结构网格上的 L、D、U 算子及求逆过程在一些方面不同，下面分别讨论。

结构网格上的 LU-SGS

结构网格的算子定义为，见文献 $ [50–52, 55, 73] $：

$$ \begin{aligned}\mathbf{L}=&\ \left(\bar{A}^{+}+\bar{A}_{\mathrm{v}}\right)_{i-1}\Delta S_{i-1/2}^{I}+\left(\bar{A}^{+}+\bar{A}_{\mathrm{v}}\right)_{j-1}\Delta S_{j-1/2}^{J}\\&+\left(\bar{A}^{+}+\bar{A}_{\mathrm{v}}\right)_{k-1}\Delta S_{k-1/2}^{K},\\\mathbf{U}=&\ \left(\bar{A}^{-}-\bar{A}_{\mathrm{v}}\right)_{i+1}\Delta S_{i+1/2}^{I}+\left(\bar{A}^{-}-\bar{A}_{\mathrm{v}}\right)_{j+1}\Delta S_{j+1/2}^{J}\\&+\left(\bar{A}^{-}-\bar{A}_{\mathrm{v}}\right)_{k+1}\Delta S_{k+1/2}^{K},\\\mathbf{D}=&\ \frac{\Omega}{\Delta t}\bar{I}+\left(\bar{A}^{-}-\bar{A}_{\mathrm{v}}\right)\Delta S_{i-1/2}^{I}+\left(\bar{A}^{-}-\bar{A}_{\mathrm{v}}\right)\Delta S_{j-1/2}^{J}\\&+\left(\bar{A}^{-}-\bar{A}_{\mathrm{v}}\right)\Delta S_{k-1/2}^{K}+\left(\bar{A}^{+}+\bar{A}_{\mathrm{v}}\right)\Delta S_{i+1/2}^{I}\\&+\left(\bar{A}^{+}+\bar{A}_{\mathrm{v}}\right)\Delta S_{j+1/2}^{J}+\left(\bar{A}^{+}+\bar{A}_{\mathrm{v}}\right)\Delta S_{k+1/2}^{K}-\frac{\partial(\Omega\vec{Q})}{\partial\vec{W}}.\end{aligned}\tag{6.51} $$

为便于阅读，式（6.51）只列出不同于 i、j、k 的节点下标，单元中心格式则为单元下标。$ \Delta S $ 的上标 i、j、k 表示计算空间方向。正／负通量 Jacobian $ \bar{A}^{\pm} $ 和黏性通量 Jacobian $ \bar{A}_{v} $ 的单位法向量，与相应面面积 $ \Delta S $ 在控制体的同一侧计算。注意，这里假设单位法向量指向控制体外部；而一些文献假设相对两面的单位法向量指向同一方向。

式（6.51）的黏性通量 Jacobian 可数值计算，也可按式（6.45）以 TSL 近似替代。无论实际边界层朝向如何，都可在所有计算坐标方向应用 TSL 近似。进一步简化，是按文献 [65] 的建议，用黏性谱半径（式（6.19））代替黏性通量 Jacobian，即 $ \bar{A}_{v}\Delta S\approx\hat{\Lambda}_{v} $。

分裂后的对流通量 Jacobian $ \bar{A}^{\pm} $ 构造为：正矩阵的特征值全部非负，负矩阵的特征值全部非正。一般定义为 [49]：

$$ \bar{A}^{\pm}\Delta S=\frac{1}{2}\left(\bar{A}_{\mathrm{c}}\Delta S\pm r_{A}\bar{I}\right),\quad r_{A}=\omega\hat{\Lambda}_{\mathrm{c}},\tag{6.52} $$

其中，$ \bar{A}_{c} $ 为对流通量 Jacobian（A.9 节），$ \hat{\Lambda}_{c} $ 为其谱半径，由式（4.53）或式（6.15）给出。注意，忽略 $ \bar{A}_{c} $ 的导数后，式（6.52）与式（6.40）相似。因子 $ \omega $ 是超松弛参数，也决定隐式耗散量，因而影响收敛性质。可在区间 $ 1 < \omega \leq 2 $ 内选取。较大的 $ \omega $ 提高 LU-SGS 的稳定性，但可能减慢定常收敛。式（6.52）对 Jacobian $ \bar{A}^{\pm} $ 的定义保证系统矩阵对角占优，这对迭代求逆过程（6.50）的效率与稳健性十分重要。

式（6.52）的分裂配合平均面矢量，可以简化对角算子 D 的计算：

$$ \begin{aligned}\mathbf{D}&=\left[\frac{\Omega}{\Delta t}+\omega\left(\hat{\boldsymbol{\Lambda}}_{c}^{I}+\hat{\boldsymbol{\Lambda}}_{c}^{J}+\hat{\boldsymbol{\Lambda}}_{c}^{K}\right)\right]\bar{\boldsymbol{I}}\\&\quad+2\left(\bar{\boldsymbol{A}}_{\mathrm{v}}^{I}\Delta\boldsymbol{S}^{I}+\bar{\boldsymbol{A}}_{\mathrm{v}}^{J}\Delta\boldsymbol{S}^{J}+\bar{\boldsymbol{A}}_{\mathrm{v}}^{K}\Delta\boldsymbol{S}^{K}\right)-\frac{\partial(\Omega\vec{\boldsymbol{Q}})}{\partial\vec{\boldsymbol{W}}}.\end{aligned}\tag{6.53} $$

对流通量 Jacobian 的谱半径 $ \hat{\Lambda}_{c} $ 见式（6.15）。面面积与法向量按式（6.16），在相应 I、J、K 方向平均。下文将看到，这一近似可明显减少运算次数和内存需求。

LU-SGS 方法的特色，是式（6.50）中前向及后向扫掠的实施方式。二维情况下，沿计算空间的对角线 $ (i+j)= $ 常数扫掠。图 6.5 显示前向扫掠，即式（6.50）的第一行。这样，L、U 算子的非对角项从之前的扫掠中已知，图中以叉号表示。三维情况下，在 $ i+j+k= $ 常数的平面上求逆隐式算子，见图 6.6。因此，LU-SGS 格式可写为：

$$ \begin{aligned}\mathbf{D}\;\Delta\vec{W}_{i,j,k}^{(1)}&=-\vec{R}_{i,j,k}^{n}-\mathbf{L}\;\Delta\vec{W}^{(1)},\\\mathbf{D}\;\Delta\vec{W}_{i,j,k}^{n}&=\mathbf{D}\;\Delta\vec{W}_{i,j,k}^{(1)}-\mathbf{U}\;\Delta\vec{W}^{n}.\end{aligned}\tag{6.54} $$

式（6.54）表明，唯一需要求逆的是对角项 D。因此，LU-SGS 将稀疏带状矩阵的求逆转化为分块对角矩阵的求逆。若进一步用黏性谱半径近似式（6.53）的黏性通量 Jacobian，除源项外，D 将成为对角矩阵。因此，与其他隐式格式，如前述 ADI 格式相比，其计算量很小。此外，对角平面上各节点／单元的对角算子可独立求逆，易于向量化。对角平面上节点／单元的下标，可由以下伪代码确定 $ [79] $：

$$ \begin{aligned}&\text{遍历平面 }\mathrm{plane}=1,\ldots,\mathrm{nplanes}:\\&\quad\text{遍历 }k=1,\ldots,\mathrm{kmax}:\\&\qquad\text{遍历 }j=1,\ldots,\mathrm{jmax}:\\&\qquad\quad\text{遍历 }i=1,\ldots,\mathrm{imax}:\\&\qquad\qquad\text{若 }i+j+k=\mathrm{plane}+2,\text{ 则存储下标。}\\&\qquad\quad\text{结束 }i\text{ 循环。}\\&\qquad\text{结束 }j\text{ 循环。}\\&\quad\text{结束 }k\text{ 循环。}\\&\text{结束平面循环。}\end{aligned} $$

对角平面的数量为 $ n_{planes} = \max + j\max + k\max - 2 $。显然，上述代码还可优化，以提高计算效率。

为避免显式计算并存储 L、U 中的对流通量 Jacobian，可用通量的 Taylor 展开替代乘积 $ \bar{A}^{\pm}\Delta\vec{W}^{n} $ [53]。由式（6.52），可写为：

$$ (\bar{A}^{\pm}\Delta S)\Delta\vec{W}^{n}\approx\frac{1}{2}\left(\Delta\vec{F}_{\mathrm{c}}\Delta S\pm r_{A}\bar{I}\Delta\vec{W}^{n}\right)\tag{6.55} $$

其中，对流通量的增量为：

$$ \Delta\vec{F}_{\mathrm{c}}=\vec{F}_{\mathrm{c}}^{n+1}-\vec{F}_{\mathrm{c}}^{n}.\tag{6.56} $$

由于沿对角平面扫掠时 $ \vec{F}_{c}^{n+1} $ 已知，式（6.55）的简化才成为可能。这进一步大幅降低 LU-SGS 的计算量。

时间步长 $ \Delta t $ 可按 6.1.4 节的方法，用式（6.14）计算。不过，Rieger 和 Jameson [53] 指出，当 $ \Delta t \to \infty $ 时，隐式 LU-SGS 格式（6.49）相当于近似 Newton 迭代。因此，实际定常流动通常采用 $ 10^{4} $ 至 $ 10^{6} $ 量级的 CFL 数；此时由超松弛参数 $ \omega $ 控制收敛。非定常流动可采用下文 6.3 节的形式，也可采用文献 [80] 的修正 LU-SGS 格式。

非结构网格上的 LU-SGS

对中位对偶格式，算子为 $ [65-67] $：

$$ \begin{aligned}\mathbf{L}&=\sum_{j\in L(i)}\left[\bar{A}_{j}^{+}+(\bar{A}_{\mathrm{v}})_{j}\right]\Delta S_{ij},\\\mathbf{U}&=\sum_{j\in U(i)}\left[\bar{A}_{j}^{-}-(\bar{A}_{\mathrm{v}})_{j}\right]\Delta S_{ij},\\\mathbf{D}&=\frac{\Omega_{i}}{\Delta t_{i}}\bar{I}+\frac{\omega}{2}(\hat{\Lambda}_{\mathrm{c}})_{i}+\sum_{j=1}^{N_{F}}(\bar{A}_{\mathrm{v}})_{i}\Delta S_{ij}-\frac{\partial(\Omega_{i}\vec{Q}_{i})}{\partial\vec{W}}.\end{aligned}\tag{6.57} $$

其中，$ L(i) $ 和 $ U(i) $ 表示节点 i 属于下三角或上三角矩阵部分的最近邻，$ \Delta S_{ij} $ 为边 ij 对应的面面积（见图 5.9），$ N_{F} $ 为控制体 $ \Omega_{i} $ 的面数。对流通量谱半径 $ (\hat{\Lambda}_{c})_{i} $ 按式（6.21）计算，黏性通量 Jacobian $ \bar{A}_{v} $ 仍可用其谱半径近似 [65]。此时，对角算子为：

$$ \mathbf{D}=\frac{\Omega_{i}}{\Delta t_{i}}\bar{I}+\frac{\omega}{2}(\hat{\Lambda}_{\mathrm{c}})_{i}+(\hat{\Lambda}_{\mathrm{v}})_{i}-\frac{\partial(\Omega_{i}\vec{Q}_{i})}{\partial\vec{W}},\tag{6.58} $$

其中，$ (\hat{\Lambda}_{\mathrm{v}})_{i} $ 按式（6.21）计算。单元中心格式也可得到类似式（6.57）、（6.58）的公式，主要区别是 L、U 算子对单元面求和，而非对相连的边求和。

式（6.57）的集合 $ L(i) $ 和 $ U(i) $ 应发挥与结构网格对角平面相同的作用。因此，必须将节点或单元分层，使其满足 [65]：

• 当前层的节点 i 或单元 I，与流动变量已经更新的层相连；否则 LU-SGS 会退化为 Jacobi 迭代；

• 同一层内的节点或单元彼此不相连；否则格式无法向量化。

中位对偶格式可按文献 $ [65] $ 的过程生成层，单元中心格式的相应方法见文献 $ [81, 82] $。

适当定义集合 $ L(i) $ 与 $ U(i) $ 后，可采用以下两步求逆过程 $ [65–67] $：

$$ \begin{aligned}\mathbf{D}\Delta\vec{W}_{i}^{(1)}&=-\vec{R}_{i}^{n}-\sum_{j\in L(i)}\frac{1}{2}\Big[(\Delta F_{\mathrm{c}}^{(1)})_{j}\Delta S_{ij}+(r_{A}^{*})_{j}\bar{I}\Delta\vec{W}_{j}^{(1)}\Big],\\\mathbf{D}\Delta\vec{W}_{i}^{n}&=\mathbf{D}\Delta\vec{W}_{i}^{(1)}-\sum_{j\in U(i)}\frac{1}{2}\left[(\Delta F_{\mathrm{c}}^{n})_{j}\Delta S_{ij}-(r_{A}^{*})_{j}\bar{I}\Delta\vec{W}_{j}^{n}\right],\end{aligned}\tag{6.59} $$

其中，黏性 Jacobian 由谱半径近似，正／负 Jacobian 按式（6.55）线性化。因子 $ (r_{A}^{*})_{j} $ 定义为：

$$ \begin{aligned}(r_{A}^{*})_{j}=&\omega\left(|\vec{v}_{j}\cdot\vec{n}_{ij}|+c_{j}\right)\Delta S_{ij}\\&+\frac{\Delta S_{i,j}}{\|\vec{r}_{j}-\vec{r}_{i}\|_{2}}\left[\max\left(\frac{4}{3\rho_{j}},\frac{\gamma_{j}}{\rho_{j}}\right)\left(\frac{\mu_{L}}{\mathrm{Pr}_{L}}+\frac{\mu_{T}}{\mathrm{Pr}_{T}}\right)_{j}\right]\end{aligned}\tag{6.60} $$

其中，$ \|\vec{r}_{j}-\vec{r}_{i}\|_{2} $ 为边 ij 的长度。

时间步长可按显式格式的式（6.20）计算。不过，应略去黏性特征值，因为黏性项已包含在隐式算子中，不再像显式格式那样缩小时间步长。定常流动的 CFL 数可取 $ 10^{4} $ 至 $ 10^{6} $，再用超松弛参数 $ \omega $ 调节收敛速度与稳健性。

6.2.5 Newton–Krylov 方法

首先，将式（6.28）的隐式格式改写为：

$$ \bar{J}\Delta\vec{W}^{n}=-\vec{R}^{n},\tag{6.61} $$

其中，$ \bar{J} $ 为隐式算子（系统矩阵）。前文已知，$ \bar{J} $ 是大型稀疏、通常非对称的矩阵。前几节讨论了两种将 $ \bar{J} $ 分解为若干因子的方法，各因子比 $ \bar{J} $ 本身更易求逆。但由于因子分解误差，以及 $ \vec{R}^{n+1} $ 的近似线性化，定常收敛仅为线性。要获得 Newton 方法求解非线性方程时的二次收敛，必须满足四个条件：

• 残差的线性化必须精确；

• 必须准确地对 $ \bar{J} $ 求逆；

• 时间步长必须为 $ \Delta t\to\infty $；

• 初始解必须在某种意义上接近最终解。

显然，主要困难在于完整系统矩阵的线性化和求逆。

所谓 Krylov 子空间方法，是特别适合大型线性方程组求解的一类迭代技术。已有多种方法用于 CFD 矩阵的求逆，例如共轭梯度平方（CGS）法 $ [83] $、稳定双共轭梯度（Bi-CGSTAB）法 $ [84] $、无转置准最小残差（TFQMR）法 $ [85] $。不过，最成功的是广义最小残差（GMRES）技术，最初由 Saad 和 Schulz 提出 $ [86, 87] $，此后又由一些研究者改进和扩展 $ [88–91] $。由于应用广泛，下面重点介绍 GMRES，但大部分讨论也适用于其他 Krylov 子空间方法。

GMRES 方法

如 3.2.2 节所述，亦见 A.12 节，GMRES 在由 m 个标准正交矢量（搜索方向）构成的空间内，最小化全局残差的范数 $ \|\vec{J}\Delta\vec{W}^{n}+\vec{R}^{n}\| $。这些矢量张成式（3.10）的 Krylov 子空间 $ K_{m} $。算法可概括为：

1. 猜测初始解 $ \vec{W}_{0}^{n} $，计算初始残差矢量 $ \vec{r}_{0}=\bar{J}\Delta\vec{W}_{0}^{n}+\vec{R}^{n} $；

2. 通过 Gram–Schmidt 正交化生成 m 个搜索方向；

3. 求解最小化问题；

4. 按下式构造式（6.61）的近似解 $ \Delta\vec{W}^{n}=\Delta\vec{W}_{0}^{n}+\vec{\gamma}_{m} $：

内存需求随搜索方向数线性增长，因此实际将 m 限制在 10 至 40。这可能不足以得到收敛解 $ \Delta\vec{W}^{n} $，于是必须重启 GMRES：令 $ \Delta\vec{W}_{0}^{n} = \Delta\vec{W}^{n} $，计算 $ \vec{r}_{0} $，再从第 2 步继续。文献 $ ^{[92]} $ 指出，搜索方向数不宜固定；当全局残差范数低于规定容限时，应减小 m。这样，可在 Newton 迭代后期节省大量运算。

通量 Jacobian 的计算

GMRES 及其他 Krylov 子空间方法可以避免显式计算和存储通量 Jacobian $ \partial\vec{R}/\partial\vec{W} $。原因是，它们只需要形式为 $ \bar{J}\Delta\vec{W}^{n} $ 的矩阵—矢量乘积，不需要显式给出矩阵 $ \bar{J} $。通量 Jacobian 与解增量的乘积可用简单有限差分近似：

$$ \frac{\partial\vec{R}}{\partial\vec{W}}\Delta\vec{W}^{n}\approx\frac{\vec{R}(\vec{W}^{n}+h\Delta\vec{W}^{n})-\vec{R}(\vec{W}^{n})}{h}\tag{6.62} $$

这只需计算两次残差。为尽量减小数值误差，必须谨慎选择步长 h [36]。一种特别合适的取法为 [93]：

$$ h=\frac{\sqrt{\epsilon}}{\left\|\Delta\vec{W}^{n}\right\|_{2}^{2}}\max\left\{\left|d\right|,\operatorname{typ}\left[\vec{W}^{n}\cdot\left|\Delta\vec{W}^{n}\right|\right]\right\}\operatorname{sign}(d),\tag{6.63} $$

其中，$ \epsilon $ 为机器精度，d 为内积 $ \vec{W}^{n} \cdot \Delta \vec{W}^{n} $，$ |\Delta \vec{W}^{n}| $ 为矢量 $ \Delta \vec{W}^{n} $ 的各分量取绝对值后所得矢量，typ U 为 U 的典型量级。有限差分近似除节省内存和运算外，还有更重要的优点：能方便地对高阶残差 $ \vec{R}^{n} $ 进行数值上准确的线性化，包括边界条件、限制器、源项等。因此，只需适度成本就能实现 Newton 格式的二次收敛。这类格式称为 Newton–Krylov 方法 [93–100]。

预处理

Krylov 子空间方法的效率很大程度上取决于良好的预处理器，其目的是使系统矩阵 $ \bar{J} $ 的特征值聚集在 1 附近。因此，求解式（3.11）的左或右预处理系统，而非直接求解式（6.61）。结合式（6.61）、（6.62）及条件 $ \Delta t \rightarrow \infty $，采用左预处理的 Newton–Krylov 方法为：

$$ \bar{P}_{\mathrm{L}}\frac{\vec{R}(\vec{W}^{n}+h\Delta\vec{W}^{n})-\vec{R}(\vec{W}^{n})}{h}=-\bar{P}_{\mathrm{L}}\vec{R}^{n}\tag{6.64} $$

即左预处理形式；而

$$ \begin{array}{l}\displaystyle\frac{\vec{R}(\vec{W}^{n}+h\bar{P}_{R}\Delta\vec{W}^{*})-\vec{R}(\vec{W}^{n})}{h}=-\vec{R}^{n},\\ \displaystyle\bar{P}_{R}^{-1}\Delta\vec{W}^{n}=\Delta\vec{W}^{*}\end{array}\tag{6.65} $$

为右预处理形式。两种方法的主要区别是：左预处理缩放残差 $ \vec{R}^{n} $，而右预处理不缩放。监测 Krylov 方法的收敛时必须注意这一点。

显然，预处理器应尽可能接近系统矩阵 $ (\bar{P}_{L,R} \approx \bar{J}^{-1}) $ 的逆，但另一方面，又应能以较小计算量求逆。因此，必须在 Krylov 方法的收敛速度和预处理矩阵求逆所需时间之间寻求最佳折中。最成功的预处理器之一，是允许不同填充层次的不完全 LU 分解方法 $ [101, 102] $；通常采用零填充，记为 $ \mathrm{ILU}(0) $。ILU 对黏性湍流，即刚性方程，尤其高效 $ [103] $。为在非结构网格上获得良好性能，必须重排系统矩阵元素以减小带宽。6.2.1 节已讨论 RCM 重编号策略 $ [24, 25] $。

ILU 预处理的严重缺点，是必须计算并存储矩阵 $ \bar{J} $ 的元素，见 6.2.2 节。因此，一些作者建议用 LU-SGS 格式作预处理器 $ [104, 105] $。结合式（6.62），即可完全避免构造和存储 $ \bar{J} $。不过，文献 [103] 的几个二维算例表明，以 CPU 时间衡量，GMRES 与 LU-SGS 的组合不如 GMRES 与 ILU(0) 的组合。尽管如此，LU-SGS 仍是有吸引力的替代方案，尤其在三维中，与多重网格结合时效果最好。

启动问题

隐式 Newton–Krylov 方法的时间步长为无穷大。不过，Newton 迭代开始时宜采用较小时间步长，因为初始流动解通常远离定常状态，即远离下述非线性方程的根：

$$ \vec{R}(\vec{W})=0, $$

这可能导致 Newton 迭代失败。一种改进是所谓切换演化松弛技术 $ [106] $，在式（6.61）中保留 $ \Omega/\Delta t $。时间步长按显式格式的方法计算，即式（6.14）或式（6.20），但略去黏性特征值。CFL 数 $ \sigma $ 从较小初始值出发，随残差二范数的下降而增大，即：

$$ \sigma^{n+1}=\sigma^{n}\frac{\|\vec{R}^{n-1}\|_{2}}{\|\vec{R}^{n}\|_{2}}.\tag{6.66} $$

于是，迭代过程（6.61）起初具有线性收敛；CFL 数很大时，趋于 Newton 方法的二次收敛。时间项还增强 $ \bar{J} $ 的对角占优性，其贡献与 $ \Omega/\Delta t $ 成反比，有助于稳定迭代。文献 [95] 建议对式（6.66）的 $ \sigma^{n+1} $ 作截断，使每次最多增大两倍，减小时不超过十倍。

另一种解决 Newton 启动问题的办法是网格序列法：先在一系列较粗网格上求初始解，再插值到细网格。也可用计算成本低且稳健的迭代格式获得初始猜测。例如，先采用 LU-SGS 驱动的多重网格格式，再切换到以 LU-SGS 为预处理器的 GMRES。对黏性湍流，这尤其有意义，因为多重网格在初始阶段之后通常收敛变慢，但此时整体流动解已接近定常状态。

6.2.6 隐式 Runge–Kutta 格式

求解关于时间的常微分方程组（6.1），一般 m 级格式可写为：

$$ \vec{W}_{I}^{n+1}=\vec{W}_{I}^{n}-\frac{\Delta t}{(\Omega\bar{M})_{I}}\sum_{k=1}^{m}b_{k}\vec{R}_{I}^{(k)},\tag{6.67} $$

其中，每一级残差 $ \vec{R}_{I}^{(k)} $ 都是时间层及其他中间残差加权和的函数，即：

$$ \begin{array}{l}{\vec{R}_{I}^{(k)}=\vec{R}_{I}\left(t^{n}+c_{k}\Delta t,\vec{W}_{I}^{n}-\frac{\Delta t}{(\Omega\bar{M})_{I}}\sum_{l=1}^{m}a_{k l}\vec{R}_{I}^{(l)}\right)}\\ {\quad=\vec{R}_{I}\left(t^{n}+c_{k}\Delta t,\vec{W}_{I}^{(k)}\right).}\end{array}\tag{6.68} $$

系数 $ a_{kl} $、$ b_{k} $、$ c_{k} $ 可方便地写成所谓 Butcher 表 [15]：

$$ \begin{array}{c|cccc}c_1&a_{11}&a_{12}&\cdots&a_{1m}\\c_2&a_{21}&a_{22}&\cdots&a_{2m}\\\vdots&\vdots&\vdots&\ddots&\vdots\\c_m&a_{m1}&a_{m2}&\cdots&a_{mm}\\\hline&b_1&b_2&\cdots&b_m\end{array}\tag{6.69} $$

系数 $ c_{k} $ 定义中间时间层，$ a_{kl} $ 为第 k 级的权重，$ b_{k} $ 为最终级的权重。这些系数决定格式的精度、稳定性和效率。不同取值可构成不同 Runge–Kutta 格式。例如，取 $ a_{kl} = 0 $，对全部 $ l \geq k $ 均如此，得到纯显式方法，如下表所示的经典四阶 Runge–Kutta 格式：

$$ \begin{array}{c|cccc}0&0&0&0&0\\\frac12&\frac12&0&0&0\\\frac12&0&\frac12&0&0\\1&0&0&1&0\\\hline&\frac16&\frac13&\frac13&\frac16\end{array}\tag{6.70} $$

另一方面，若 $ a_{kl} \neq 0 $ 至少对某些 $ l \geq k $ 成立，即表（6.69）的对角线或上三角部分存在非零系数，则得到隐式 Runge–Kutta 格式。此时，式（6.68）的第 k 级中间残差 $ \vec{R}_{I}^{(k)} $ 也隐式依赖本级解，还可能依赖后续级的解。CFD 中通常避免与后续级耦合。仅在下三角部分具有非零 $ a_{kl} $ 的格式，称为对角隐式 Runge–Kutta（DIRK）格式；若所有对角系数 $ a_{kk} $ 还都相同，则称为单一对角隐式 Runge–Kutta（SDIRK）格式。

系数必须满足一定条件。首先，物理上合理的格式要求 $ c_{k} = \sum_{l} a_{kl} $。此外，还必须满足以下阶数条件：

$$ \sum_{k=1}^{m}b_{k}=1\tag{6.71} $$

保证一阶精度；

$$ \sum_{k=1}^{m}b_{k}c_{k}=\frac{1}{2}\tag{6.72} $$

保证二阶精度；

$$ \begin{array}{l}\displaystyle\sum_{k=1}^{m}b_{k}c_{k}^{2}=\frac{1}{3}\\\displaystyle\sum_{k=1}^{m}\sum_{l=1}^{m}b_{k}a_{kl}c_{l}=\frac{1}{6}\end{array}\tag{6.73} $$

保证三阶精度；以及

$$ \begin{array}{l}\displaystyle\sum_{k=1}^{m}b_{k}c_{k}^{3}=\frac{1}{4}\\\displaystyle\sum_{k=1}^{m}\sum_{l=1}^{m}b_{k}c_{k}a_{kl}c_{l}=\frac{1}{8}\\\displaystyle\sum_{k=1}^{m}\sum_{l=1}^{m}b_{k}a_{kl}c_{l}^{2}=\frac{1}{12}\\\displaystyle\sum_{k=1}^{m}\sum_{l=1}^{m}\sum_{n=1}^{m}b_{k}a_{kl}a_{ln}c_{n}=\frac{1}{24}\end{array}\tag{6.74} $$

保证四阶精度。满足这些条件后，可进一步优化系数，以获得所需精度、稳定性及阻尼性质。

隐式 Runge–Kutta 格式的式（6.68）中，中间残差 $ \vec{R}_{I}^{(k)} $ 隐式依赖尚未计算的残差，即 $ \vec{R}_{I}^{(l)} $，此时 $ l \geq k $。因此，每一级都须求解一组非线性方程，通常采用 Newton 方法。以对角隐式格式为例，第 k 级中间解为：

$$ \vec{W}_{I}^{(k)}=\vec{W}_{I}^{n}-\frac{\Delta t}{(\Omega\bar{M})_{I}}\left[\sum_{l=1}^{k-1}a_{k l}\vec{R}_{I}^{(l)}+a_{k k}\vec{R}_{I}^{(k)}\right].\tag{6.75} $$

式（6.75）的隐式项 $ \vec{R}_{I}^{(k)} $ 是 $ \vec{W}_{I}^{(k)} $ 与 $ (t^{n} + c_{k} \Delta t) $ 的函数，可用式（6.27）线性化为：

$$ \vec{R}_{I}^{(k)}\approx\vec{R}_{I}^{n}+\bar{J}_{I}\Delta\vec{W}^{(k)}\tag{6.76} $$

其中，$ \Delta\vec{W}^{(k)}=\vec{W}^{(k)}-\vec{W}^{n} $，$ \vec{J}=\partial\vec{R}/\partial\vec{W} $ 为 Jacobian 矩阵。代入式（6.75），得到第 k 级解的表达式：

$$ \left[\frac{(\Omega\bar{M})_{I}}{\Delta t}+a_{k k}\bar{J}_{I}\right]\Delta\vec{W}^{(k)}=-\sum_{l=1}^{k-1}a_{k l}\vec{R}_{I}^{(l)}-a_{k k}\vec{R}_{I}^{n}.\tag{6.77} $$

它在形式上等价于式（6.28），可用前几小节的任一隐式方法求解。

当微分方程（6.1）具有刚性，且必须进行准确的时间推进时，适合采用隐式 Runge–Kutta 格式。刚性可由湍流建模、高度拉伸的网格单元，或真实气体效应造成。Cooper 和 Sayfy [107, 108] 提出的思路，是将方程分为刚性项和非刚性项：前者隐式处理，后者同时显式积分。这类多级方法称为加性 Runge–Kutta 格式。与全部项隐式处理相比，这种分裂可能大幅简化式（6.77）的 Jacobian $ \bar{J} $，从而明显降低计算量。Kennedy 和 Carpenter [109] 为此设计了优化格式，其耦合显式—隐式系数见表 6.3、6.4。两种方法均为六级，具有四阶时间精度。为便于耦合，它们采用相同的 $ b_{k} $、$ c_{k} $ 系数。第二组最终级权重 $ \hat{b}_{k} $ 在两种方法中都定义一个三阶格式，称为嵌入格式，用于估计积分误差并调整时间步长。时间误差计算为：

$$ \epsilon_{I}=\frac{\Delta t}{(\Omega\bar{M})_{I}}\sum_{k=1}^{m}(b_{k}-\hat{b}_{k})\vec{R}_{I}^{(k)}.\tag{6.78} $$

表 6.5–6.7 列出其他可能特别有意义的对角隐式 Runge–Kutta 格式系数。实现与计算结果实例见文献 $ [113–121] $。

表 6.5 单一对角隐式三级四阶 Runge–Kutta 格式的系数 $ [110] $。参数 $ a=(2/\sqrt{3})\cos(\pi/18) $。该格式即使采用大时间步长，也具有良好的阻尼性质。