3.2 时间离散

3.2 时间离散

如本章开头所述，求解 Euler 方程和 Navier–Stokes 方程的绝大多数数值格式采用线方法，即分别进行空间与时间离散。这种方法具有最大的灵活性，因为可以根据待求问题的需要，方便地选择不同程度的近似，

分别用于对流通量、黏性通量和时间积分。因此，本书采用这一方法。关于将时间和空间离散相结合的其他方法，可参阅文献 $ ^{[63]} $。

将线方法应用于控制方程（2.19），并对每个控制体写出方程后，可得到一组关于时间的耦合常微分方程：

$$ \frac{\mathrm{d}(\Omega\bar{M}\vec{W})}{\mathrm{d}t}=-\vec{R}.\tag{3.3} $$

为使表达清晰，这里省略了单元下标。式（3.3）中，$ \Omega $ 表示控制体的体积；$ \vec{R} $ 表示包含源项的完整空间（有限体积）离散，即残差。残差是守恒变量 $ \vec{W} $ 的非线性函数。最后，$ \bar{M} $ 是所谓的质量矩阵。对于单元顶点格式，它将控制体内 $ \vec{W} $ 的平均值与相应内部节点及相邻节点上的点值联系起来 [212, 213]。对于单元中心格式，可将质量矩阵替换为单位矩阵，而不损害格式的时间精度。在均匀网格上应用单元顶点格式时也是如此，因为此时节点与控制体的形心重合。质量矩阵只依赖网格，并使微分方程组（3.3）相互耦合。对于不要求时间精度的定常问题，可将质量矩阵“集中化”，即替换为单位矩阵，从而避免代价较高的 $ \bar{M} $ 求逆，并使方程组（3.3）解耦。这里必须认识到：在定常状态下，解的精度完全由残差的近似阶次决定。因此，只有在将单元顶点（中点对偶）格式用于非定常流时，质量矩阵才变得重要。

假定网格静止，可将体积 $ \Omega $ 和质量矩阵移到时间导数之外。随后，可以用下列非线性格式近似时间导数 [63]：

$$ \frac{\Omega\bar{M}}{\Delta t}\Delta\vec{W}^{n}=-\frac{\beta}{1+\omega}\vec{R}^{n+1}-\frac{1-\beta}{1+\omega}\vec{R}^{n}+\frac{\omega\Omega\bar{M}}{(1+\omega)\Delta t}\Delta\vec{W}^{n-1}\tag{3.4} $$

其中，

$$ \Delta\vec{W}^{n}=\vec{W}^{n+1}-\vec{W}^{n}\tag{3.5} $$

为解的修正量。上标 n 和 $ (n+1) $ 表示时间层（n 为当前时间层），$ \Delta t $ 表示时间步长。当满足条件

$$ \beta=\omega+\frac{1}{2}\tag{3.6} $$

时，式（3.4）具有二阶时间精度；否则，时间精度降为一阶。通过设置参数 $ \beta $ 和 $ \omega $，可以得到显式（$ \beta = 0 $）或隐式时间推进格式。下面将简要讨论这两大类格式，第 6 章再作详细介绍。

3.2.1 显式格式

在式（3.4）中令 $ \beta=0 $、$ \omega=0 $，可得到一种基本的显式时间积分格式。此时，时间导数采用前向差分近似，而残差只在当前时间层计算（基于已知流动量），即：

$$ \Delta\vec{W}^{n}=-\frac{\Delta t}{\Omega}\vec{R}^{n},\tag{3.7} $$

这里已对质量矩阵作集中化处理。这是单级格式，因为仅需计算一次残差，就能得到新解 $ \vec{W}^{n+1} $。式（3.7）的格式没有实际应用价值，因为它只有与一阶迎风空间离散相结合时才稳定。

多级时间推进格式（Runge–Kutta 格式）非常流行，它通过多个级推进解 $ [143] $，并在中间状态计算残差。各级残差用相应系数加权。可以优化这些系数，以扩大稳定域并改善格式的阻尼特性，进而提高收敛速度和稳健性 $ [143, 214, 215] $。此外，根据级系数与级数的选择，多级格式可以达到二阶或更高阶时间精度。研究者还设计了专门的 Runge–Kutta 格式，在尽可能增大允许时间步长的同时，保持 TVD 和 ENO 空间离散方法的性质 $ [216] $。

显式多级时间推进格式可以与任何空间离散格式结合，容易在串行、向量及并行计算机上实现。显式格式计算成本低，内存需求也较小。然而，由于稳定性限制，最大允许时间步长受到严格约束。尤其对于黏性流动和高度拉伸的网格单元，向定常状态的收敛会明显减慢。此外，对于刚性方程组（例如真实气体模拟、湍流模型）或刚性源项，达到定常状态可能需要极长时间。更糟的是，显式格式甚至可能失稳，或产生虚假的定常解 $ [217] $。

如果只关心定常解，可以选择或组合多种收敛加速技术。第一种也是非常常用的方法是局部时间步进：在每个控制体内，以当地最大允许时间步长推进解。这样可以显著加快向定常状态的收敛，但中间瞬态解不再具有时间精度。另一种方法是特征时间步进，它不仅使用随位置变化的时间步长，还用各自的时间步长积分连续性、动量和能量方程。文献 [218] 展示了这一思想在二维 Euler 方程中的潜力。另一项与特征时间步进类似的加速技术是 Jacobi 预处理 $ [219–221] $。它本质上是在 Runge–Kutta 格式的每一级进行点隐式 Jacobi 松弛。Jacobi 预处理可以看作一种将全部波分量（通量 Jacobian 的特征值）缩放到相同有效速度的时间步进，也为基本显式格式引入了隐式成分。

另一种非常常用的加速方法，是在显式格式中引入一定程度的隐式处理，从而增大可用的最大时间步长；它称为隐式残差光顺或残差平均 $ [222, 223] $。在结构网格上，该方法需要为每个守恒变量求解一个三对角矩阵系统；在非结构网格上，通常用 Jacobi 迭代实现矩阵求逆。标准隐式残差光顺可以使时间步长增大到原来的 2–3 倍。研究者还开发了其他隐式残差光顺技术，例如与迎风空间离散配合使用的迎风隐式残差光顺方法 $ [224] $。相较于标准方法，它允许使用显著更大的时间步长，并提高时间推进过程的稳健性 $ [225] $。另一个方法是隐显式残差光顺 $ [226, 227] $，旨在改善较大时间步长下时间离散的阻尼特性。

最后需要介绍的，也可能是最重要的收敛加速技术，是多重网格方法。该方法于 20 世纪 60 年代由俄罗斯的 Fedorenko $ [228] $ 和 Bakhvalov $ [229] $ 提出，用于求解椭圆型边值问题；随后 Brandt $ [230, 231] $ 对其进行了发展和推广。多重网格的思想基于这样的观察：迭代格式通常能非常有效地消除解中的高频误差（即控制体之间的振荡），但对低频（即全局）误差的抑制较差。因此，在某一网格上推进解后，将解传递到更粗的网格；原来的低频误差在那里部分转变为高频误差，因而又能被迭代求解器有效衰减。这一过程在一系列逐渐变粗的网格上递归进行，每一级多重网格有助于抑制某一频带的误差成分。到达最粗网格后，逐级汇集解的修正量，并将其插值回最初的细网格，再更新解。重复整个多重网格循环，直到解的变化小于给定阈值。为进一步加快收敛，可以从粗网格开始多重网格过程，进行若干循环后，将解传递到更细的网格，再次进行多重网格循环；依次重复，直至到达最细网格。这种方法称为全多重网格（FMG）$ [231] $。

如前所述，多重网格方法最初用于求解椭圆型边值问题（Poisson 方程），在此类问题中效率很高。Jameson 首先提出将多重网格用于 Euler 方程的求解 $ [222, 232] $。该方法以所谓的全近似存储（FAS）格式 $ [231] $ 为基础，将多重网格直接应用于非线性控制方程。如今，多重网格已成为求解 Navier–Stokes 方程的标准加速技术。结构网格上的实现实例见文献 $ [233–239] $，非结构网格上的实例见文献 $ [240–250] $。尽管不如求解椭圆型微分方程时那样快，许多结果仍表明，多重网格可以使 Euler 或 Navier–Stokes 方程的求解加快 5–10 倍。图 3.7 给出了跨声速流动的实例。近期研究还发现，将控制方程分解为双曲型与椭圆型部分，可以实现更快的收敛 $ [251] $。第 9.4 节将再次讨论多重网格方法。

3.2.2 隐式格式

在式（3.4）中令 $ \beta\neq0 $，即可得到一族隐式时间积分格式。对于非定常流模拟，具有 $ \beta=1 $ 和 $ \omega=1/2 $ 的三点隐式后向差分格式非常常用，它具有二阶时间精度。此时，该格式通常用于所谓的双时间步进方法 [212, 213, 252–254]，即在每一个物理时间步内，以伪时间求解一个定常问题。

对于定常流动问题，采用 $ \omega=0 $ 的格式更合适，因为内存需求较低。将式（3.4）中的残差 $ \vec{R}^{n+1} $ 在当前时间层线性化，可得到格式：

$$ \left(\bar{M}\frac{\Omega}{\Delta t}+\beta\frac{\partial\vec{R}}{\partial\vec{W}}\right)\Delta\vec{W}^{n}=-\vec{R}^{n}.\tag{3.8} $$

其中，$ \partial\vec{R}/\partial\vec{W} $ 称为通量 Jacobian，是一个大型稀疏矩阵。式（3.8）左侧括号内的表达式也称为隐式算子。如前所述，将质量矩阵 $ \bar{M} $ 替换为单位矩阵，不会影响定常解。通常将式（3.8）的参数 $ \beta $ 设为 1，得到一阶时间精度；取 $ \beta=1/2 $ 时，则得到二阶时间精度。然而，不建议这样做，因为 $ \beta=1 $ 的格式稳健得多，而且对定常问题而言，时间精度本来就不起作用。

隐式格式相对于显式格式的主要优势是，可以采用显著更大的时间步长，而不损害时间积分过程的稳定性。事实上，当 $ \Delta t \to \infty $ 时，格式（3.8）转化为具有二次收敛性的标准 Newton 方法。不过，二次收敛的前提是通量 Jacobian 包含残差的完整线性化。隐式格式的另一个重要优势，是在刚性方程组和／或刚性源项问题中具有更好的稳健性和收敛速度；这类问题常见于真实气体模拟、湍流建模及高度拉伸网格（高 Reynolds 数流动）。另一方面，隐式格式的收敛越快（按时间步数或迭代次数衡量）、稳健性越好，每一时间步或迭代的计算量通常也越大。因此，采用多重网格加速的显式格式可能同样高效，甚至更高效。此外，隐式格式的向量化和并行化显著比显式格式困难。

对每个控制体写出式（3.8）的隐式格式后，会得到一个大型线性方程组；为获得更新量 $ \Delta\vec{W}^{n} $，必须在每个时间步 $ \Delta t $ 求解它。这一任务可以采用直接法或迭代法完成。

直接法利用 Gauss 消元或某种稀疏矩阵直接法 $ [255, 256] $，对式（3.8）左侧进行精确求逆。尽管在结构网格 $ [257–260] $ 和非结构网格 $ [261] $ 上均展示了二次收敛，直接法仍不适用于三维问题，因为它需要极大的计算量和内存。

因此，对于较大的网格或三维问题，唯一实际可行的途径是迭代法。在每个时间步内，用某种迭代求逆方法求解线性系统，得到 $ \Delta\vec{W}^{n} $。为降低内存需求并增强对角占优性，通常通过对右侧一阶空间精度的离散作线性化来构造通量 Jacobian $ \partial\vec{R}/\partial\vec{W} $。这一近似主要带来两个后果：无法达到 Newton 格式的二次收敛，且最大时间步长受到限制。然而，每次迭代的计算量明显降低，因而可得到计算效率较高的格式。

在结构网格上，主要采用交替方向隐式（ADI）格式 $ [262–265] $、（线）Jacobi 或 Gauss–Seidel 松弛格式 $ [266–270] $，特别是下—上对称 Gauss–Seidel 格式（LU-SGS；也称 LU-SSOR，即下—上对称逐次超松弛）$ [271–275] $。这些方法均将隐式算子分解为若干部件之和或之积，各部件更容易求逆。由于存在分解误差（与原矩阵之间的差异），且通量 Jacobian 也作了简化，将线性系统求得非常精确并不划算。事实上，ADI 和 LU-SGS 方法在每个时间步只进行一次迭代。

非结构网格上的隐式迭代法大多基于 Gauss–Seidel 松弛格式 $ [276–279] $。可以用红黑 Gauss–Seidel 方法改善收敛，其在非结构网格上的扩展见文献 $ [280–282] $。将 LU-SGS 格式用于非结构网格也是一种尤其值得关注的选择 $ [27, 283, 284] $，因为其内存需求和计算量都很低。

鉴于线隐式方法在结构网格上的成功，一些研究尝试将其用于非结构网格 $ [285, 286] $。其思路是构造连续的线，使每个网格点或每个网格单元（对于单元中心格式）仅被访问一次，即 Hamilton 巡回路径 $ [287] $。这些线主要沿坐标方向布置，但在边界及其他必要位置折返，因此被戏称为“蛇形线”。随后，采用三对角求解器对式（3.8）左侧求逆。后来发现，线的折返会减慢收敛；为克服这一问题，将每条线分成多条短线段（linelet）$ [288] $。然而，这种方法在向量计算机上的性能

较差。短线段的思想还被用于改善显式格式在高度拉伸的黏性非结构网格上的收敛，通过在横跨边界层的方向使用隐式求解器实现 $ [248] $。

更复杂的迭代技术，是以更全局的方式处理线性方程组的 Krylov 子空间方法。其发展由 Hestenes 和 Stiefel 提出的一种高效大型稀疏线性系统迭代格式——共轭梯度法——推动 [289]。原始共轭梯度法只适用于 Hermitian 正定矩阵，但对于一个 $ n \times n $ 矩阵，它至多经过 n 次迭代便可收敛。此后，研究者提出了多种 Krylov 子空间方法，用于求解 CFD 应用中出现的任意非奇异矩阵系统，例如共轭梯度平方（CGS）[290]、稳定双共轭梯度（Bi-CGSTAB）[291] 和无转置准最小残差（TFQMR）方法 [292]。

不过，应用可能最广泛的是 Saad 和 Schultz 提出的广义最小残差（GMRES）方法 $ [293] $。将式（3.8）的隐式格式改写为：

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

则 $ \bar{J} $ 是一个大型稀疏非对称矩阵（对应左侧）。从初始猜测 $ \Delta\vec{W}_{0} $ 出发，GMRES(m) 方法寻找解 $ \Delta\vec{W}^{n} $，其形式为 $ \Delta\vec{W}^{n}=\Delta\vec{W}_{0}^{n}+\vec{y}_{m} $，其中 $ \vec{y}_{m} $ 属于 Krylov 子空间：

$$ \begin{aligned}\mathcal{K}_{m}&\equiv\operatorname{span}\{\vec{r}_{0},\bar{J}\vec{r}_{0},\bar{J}^{2}\vec{r}_{0},\ldots,\bar{J}^{m-1}\vec{r}_{0}\}\\\vec{r}_{0}&=\bar{J}\Delta\vec{W}_{0}^{n}+\vec{R}^{n}.\end{aligned}\tag{3.10} $$

使残差 $ \|\bar{J}\Delta\vec{W}^{n}+\vec{R}^{n}\| $ 最小。参数 m 指定 Krylov 子空间的维数，也就是搜索方向 $ (\bar{J}^{i}\vec{r}_{0}) $ 的数量。由于必须存储全部方向，m 通常取 10–40；对病态矩阵（例如湍流或真实气体模拟中出现的矩阵），需要较大的取值。如果 m 次子迭代后仍未收敛，必须重新启动 GMRES。GMRES 的内存需求显著高于 Bi-CGSTAB 或 TFQMR，但它更稳健、收敛更平滑，通常也更快。文献 $ ^{[294]} $ 对各种方法作了非常详细的比较。

然而，与其他共轭梯度类方法一样，对于 CFD 问题，预处理是必不可少的。此时求解：

$$ (\bar{P}_{\mathrm{L}}\bar{J})\Delta\vec{W}^{n}=-\bar{P}_{\mathrm{L}}\vec{R}^{n},\quad\text{或}\quad\bar{J}\bar{P}_{R}(\bar{P}_{R}^{-1}\Delta\vec{W}^{n})=-\vec{R}^{n}\tag{3.11} $$

以取代式（3.9）的系统。矩阵 $ \bar{P}_{L} $ 和 $ \bar{P}_{R} $ 分别表示左预处理矩阵和右预处理矩阵。预处理矩阵应尽可能接近 $ \bar{J}^{-1} $，

使特征值聚集在 1 附近；同时，它当然还应容易求逆。一种特别高效的预处理方法是零填充的不完全下—上分解 $ [295] $，记为 ILU(0)。关于与 GMRES 配合使用的不同预处理技术，可参阅 $ [198, 296–299] $。

由于 GMRES 需要较多内存来存储搜索方向，有时还要存储预处理矩阵，因此避免显式构造和存储通量 Jacobian $ \partial\vec{R}/\partial\vec{W} $ 是有益的。这可以通过所谓的无矩阵方法实现。其思想基于如下观察：GMRES 及某些其他 Krylov 子空间方法只使用下列形式的矩阵—向量乘积：

$$ \frac{\partial\vec{R}}{\partial\vec{W}}\Delta\vec{W}^{n}, $$

它可简单地用有限差分近似为：

$$ \frac{\partial\vec{R}}{\partial\vec{W}}\Delta\vec{W}^{n}=\frac{\vec{R}(\vec{W}+\epsilon\Delta\vec{W}^{n})-\vec{R}(\vec{W})}{\epsilon},\tag{3.12} $$

因此只需计算残差。参数 $ \epsilon $ 必须谨慎选择，以使数值误差最小（例如见 [300] 或 [301]）。无矩阵方法的另一项、也是更重要的优势是，可方便地在隐式格式中使用高阶残差 $ \vec{R}^{n} $ 的数值精确线性化。因此，能够以适中的成本获得 Newton 格式的二次收敛。这时称为 Newton–Krylov 方法 [212, 301–305]。实践表明，在所有 Krylov 子空间方法中，GMRES 最适合无矩阵实现 [306]。一个值得关注的选择是用 LU-SGS 格式为无矩阵 GMRES 作预处理。由于 LU-SGS 同样无须显式存储通量 Jacobian，内存需求可以进一步降低。文献 [307] 在非结构网格上的三维无黏流和层流计算中展示了这一方法的计算效率。

隐式格式的收敛也可用多重网格改善，基本有两种途径。第一，在隐式格式内部使用多重网格，作为每个时间步产生的线性方程组的求解器 $ (3.9) $，或作为某种共轭梯度类方法的预处理器 $ [308, 309] $。第二，将隐式格式本身作为 FAS 多重网格方法中的光顺器，并将 FAS 直接应用于控制方程 $ [282, 310–313] $。一些研究表明，至少对于纯空气动力学问题，较“简单”的隐式格式（如 Gauss–Seidel）与多重网格结合，在 CPU 时间方面可能比 GMRES 等方法更高效 $ [282, 302] $。