6.1 显式时间推进格式

第 6 章 时间离散

采用线方法，即对控制方程（2.19）分别进行空间离散和时间离散，并为每个控制体写出方程，得到以下关于时间的耦合常微分方程组：

$$ \frac{\mathrm{d}\left(\Omega\bar{M}\vec{W}\right)_{I}}{\mathrm{d}t}=-\vec{R}_{I}.\tag{6.1} $$

式（6.1）中，$ \Omega $ 为体积，$ \vec{R} $ 为残差，$ \bar{M} $ 为质量矩阵，下标 I 表示相应控制体。必须对方程组（6.1）进行时间积分，以获得定常解（$ \vec{R}_{I}=0 $），或再现非定常流动随时间的变化。

3.2 节已简要讨论方程组（6.1）的求解。各种显式和隐式方法都可从一个基本的非线性格式导出。对静止网格，该格式为：

$$ \frac{(\Omega\bar{M})_{I}}{\Delta t_{I}}\Delta\vec{W}_{I}^{n}=-\frac{\beta}{1+\omega}\vec{R}_{I}^{n+1}-\frac{1-\beta}{1+\omega}\vec{R}_{I}^{n}+\frac{\omega}{1+\omega}\frac{(\Omega\bar{M})_{I}}{\Delta t_{I}}\Delta\vec{W}_{I}^{n-1},\tag{6.2} $$

其中，

$$ \Delta\vec{W}_{I}^{n}=\vec{W}_{I}^{n+1}-\vec{W}_{I}^{n}\tag{6.3} $$

表示解的增量（修正量）。上标 n 和 $ (n+1) $ 表示时间层。因此，$ \vec{W}^{n} $ 为当前时刻 t 的流动解，$ \vec{W}^{n+1} $ 为时刻 $ (t+\Delta t) $ 的解。参数 $ \beta $ 和 $ \omega $ 决定离散类型（显式或隐式）及时间精度。例如，要达到二阶时间精度，必须满足式（3.6）的条件。

以下较详细地介绍最常用的显式与隐式时间推进方法，并说明如何计算特定格式允许的最大时间步长。此外，还将讨论在结构网格及非结构网格上的适当实现。最后一节讨论非定常流动问题的时间精确求解。

6.1 显式时间推进格式

显式格式从已知解 $ \vec{W}^{n} $ 出发，利用相应残差 $ \vec{R}^{n} $，获得时刻 $ (t + \Delta t) $ 的新解。换言之，新解 $ \vec{W}^{n+1} $ 只依赖已经知道的数值。因此，显式格式十分简单，易于实现。

如 3.2.1 节所述，在式（6.2）中取 $ \beta = 0 $、$ \omega = 0 $，即可导出基本显式格式：

$$ \bar{M}_{I}\Delta\vec{W}_{I}^{n}=-\frac{\Delta t_{I}}{\Omega_{I}}\vec{R}_{I}^{n},\tag{6.4} $$

这称为前向 Euler 近似。对定常问题或单元中心离散，可将质量矩阵 $ \bar{M} $ 作集中处理，即用单位矩阵代替。

最常用、应用最广泛的显式方法，是多级（Runge–Kutta）时间推进格式及其变体——混合多级格式。以下介绍这两种方法。

6.1.1 多级格式（Runge–Kutta）

Jameson 等首先提出显式多级格式的概念 $ [1] $。该格式通过若干步，即所谓“级”，推进解，可视为按式（6.4）依次进行修正。对离散控制方程（6.1），在质量矩阵集中处理后，m 级格式为：

$$ \begin{aligned}\vec{W}_{I}^{(0)}&=\vec{W}_{I}^{n},\\\vec{W}_{I}^{(1)}&=\vec{W}_{I}^{(0)}-\alpha_{1}\frac{\Delta t_{I}}{\Omega_{I}}\vec{R}_{I}^{(0)},\\\vec{W}_{I}^{(2)}&=\vec{W}_{I}^{(0)}-\alpha_{2}\frac{\Delta t_{I}}{\Omega_{I}}\vec{R}_{I}^{(1)},\\&\vdots\\\vec{W}_{I}^{n+1}&=\vec{W}_{I}^{(m)}=\vec{W}_{I}^{(0)}-\alpha_{m}\frac{\Delta t_{I}}{\Omega_{I}}\vec{R}_{I}^{(m-1)}.\end{aligned}\tag{6.5} $$

其中，$ \alpha_{1}\ldots\alpha_{m} $ 为各级系数，记号 $ \vec{R}_{I}^{(k)} $ 表示利用第 k 级的解 $ \vec{W}_{I}^{(k)} $ 计算残差。

与经典 Runge–Kutta 格式不同，这里只存储第零级解和最新残差，以减少内存需求。可调节各级系数，增大特定空间离散允许的最大时间步长并改善稳定性 $ [2-4] $。相容性只要求 $ \alpha_{m}=1 $。这种 Runge–Kutta 格式修改的一个后果，是只有满足 $ \alpha_{m-1}=1/2 $ 时才能实现二阶时间精度；否则，多级格式只有一阶时间精度。

上述多级方法 $ (6.5) $ 特别适合结构网格及非结构网格上的迎风空间离散。中心离散则采用下述混合多级方法更为高效。表 6.1 给出三至五级格式用于一阶及二阶迎风空间离散的优化系数组 $ [2] $。实际经验表明，只要流场中存在强激波，不论空间离散阶数如何，都宜采用一阶格式的系数。原因是，为防止解振荡，任何高阶格式都会在激波处切换为一阶；而强激波处的残差对定常收敛的影响最大。

所有显式格式的主要缺点，是时间步长 $ (\Delta t) $ 受到控制方程的特征性质及网格几何的严格限制。6.1.4 节介绍最大允许时间步长的计算；确定时间步长和 CFL 数的理论问题将在 10.3 节的稳定性分析中讨论。

6.1.2 混合多级格式

将显式多级格式 $ (6.5) $ 用于方程组 $ (6.1) $ 时，若不在每一级重新计算黏性通量与耗散，就能明显降低计算量。此外，可混合不同级的耗散项，提高格式稳定性。Martinelli [5] 以及 Mavriplis 和 Jameson [6] 提出这类方法，称为混合多级格式。只要各级系数经过仔细优化，其稳健性可与基本多级格式相当。

以一种常用的五级混合格式为例，其耗散项在奇数级计算，通常记为 $ (5,3) $ 格式。首先，将空间离散分为两部分：

$$ \vec{R}_{I}=(\vec{R}_{c})_{I}-(\vec{R}_{\mathrm{d}})_{I}.\tag{6.6} $$

第一部分 $ \vec{R}_{c} $ 包含对流通量的中心离散，可以采用变量平均或通量平均，并包含源项。第二部分 $ \vec{R}_{d} $ 由黏性通量和数值耗散组成。例如，对采用人工耗散的中心格式（4.3.1 节或 5.3.1 节），可取：

$$ (\vec{R}_{c})_{I}=\sum_{k=1}^{N_{F}}\left[\vec{F}_{c}(\vec{W}_{\mathrm{av}})\Delta S\right]_{k}-(\vec{Q}\Omega)_{I}, $$

$$ (\vec{R}_{\mathrm{d}})_{I}=\sum_{k=1}^{N_{F}}\left[\vec{F}_{\mathrm{v}}\Delta S+\vec{D}\right]_{k}, $$

其中，$ \vec{W}_{av} $ 为面 k 左右两侧流动变量的算术平均。

按式（6.6）分裂残差后，$ (5,3) $ 格式可写为：

$$ \begin{aligned}\vec{W}_{I}^{(0)}&=\vec{W}_{I}^{n},\\\vec{W}_{I}^{(1)}&=\vec{W}_{I}^{(0)}-\alpha_{1}\frac{\Delta t_{I}}{\Omega_{I}}\left[\vec{R}_{\mathrm{c}}^{(0)}-\vec{R}_{\mathrm{d}}^{(0)}\right]_{I},\\\vec{W}_{I}^{(2)}&=\vec{W}_{I}^{(0)}-\alpha_{2}\frac{\Delta t_{I}}{\Omega_{I}}\left[\vec{R}_{\mathrm{c}}^{(1)}-\vec{R}_{\mathrm{d}}^{(0)}\right]_{I},\\\vec{W}_{I}^{(3)}&=\vec{W}_{I}^{(0)}-\alpha_{3}\frac{\Delta t_{I}}{\Omega_{I}}\left[\vec{R}_{\mathrm{c}}^{(2)}-\vec{R}_{\mathrm{d}}^{(2,0)}\right]_{I},\\\vec{W}_{I}^{(4)}&=\vec{W}_{I}^{(0)}-\alpha_{4}\frac{\Delta t_{I}}{\Omega_{I}}\left[\vec{R}_{\mathrm{c}}^{(3)}-\vec{R}_{\mathrm{d}}^{(2,0)}\right]_{I},\\\vec{W}_{I}^{n+1}&=\vec{W}_{I}^{(0)}-\alpha_{5}\frac{\Delta t_{I}}{\Omega_{I}}\left[\vec{R}_{\mathrm{c}}^{(4)}-\vec{R}_{\mathrm{d}}^{(4,2)}\right]_{I},\end{aligned}\tag{6.7} $$

其中，

$$ \begin{aligned}\vec{R}_{\mathrm{d}}^{(2,0)}&=\beta_{3}\vec{R}_{\mathrm{d}}^{(2)}+(1-\beta_{3})\vec{R}_{\mathrm{d}}^{(0)},\\\vec{R}_{\mathrm{d}}^{(4,2)}&=\beta_{5}\vec{R}_{\mathrm{d}}^{(4)}+(1-\beta_{5})\vec{R}_{\mathrm{d}}^{(2,0)}.\end{aligned}\tag{6.8} $$

式（6.7）和式（6.8）的级系数 $ \alpha_{m} $ 与混合系数 $ \beta_{m} $，在表 6.2 中分别给出中心格式及迎风格式的取值。这两组系数专门针对多重网格方法优化（9.4 节）。上述混合多级格式的性质将在 10.3 节进一步讨论。

另一种常用做法，是只在前两级计算耗散项 $ \vec{R}_{d} $，不作混合。著名的 (5,2) 格式

常与中心空间离散结合，采用表 6.2 的级系数。不过，$ (5,2) $ 格式对黏性流动和多重网格的适用性不如 $ (5,3) $ 格式。

6.1.3 源项的处理

在某些情形下，式（4.2）或式（5.2）的源项 Q 占主导地位，采用化学反应或湍流模型时常会遇到这种情况。问题在于，大源项会使流动变量在空间和时间上迅速变化。强源项引起的变化，其时间尺度远小于流动方程的时间尺度，因而显著增加控制方程的刚性。刚性定义为 Jacobian 矩阵 $ \partial\vec{R}/\partial\vec{W} $ 的最大与最小特征值之比，也可视为最大与最小时间尺度之比。

将上述显式多级格式或其他纯显式格式用于刚性方程组时，必须大幅缩小时间步长以保持时间积分稳定，因而定常收敛很慢。更严重的是，显式格式甚至可能无法得到正确的解 $ [7] $。Curtiss 和 Hirschfelder $ [8] $ 建议采用隐式方法处理源项。为说明这一做法，将式（6.4）的基本显式格式改写为以下形式，参见式（4.2）或式（5.2）：

$$ \begin{array}{l}{\displaystyle\frac{\Omega_{I}}{\Delta t_{I}}\Delta\vec{W}_{I}^{n}=-\left[\sum_{k=1}^{N_{F}}(\vec{F}_{\mathrm{c}}^{n}-\vec{F}_{\mathrm{v}}^{n})_{k}\Delta S_{k}-\Omega_{I}\vec{Q}_{I}^{n+1}\right]}\\ {\quad=-(\vec{R}_{Q})_{I}^{n}+\Omega_{I}\vec{Q}_{I}^{n+1},}\end{array}\tag{6.9} $$

这里改在新时间层 $ (n+1) $ 计算源项。为简单起见，式（6.9）略去质量矩阵。由于时刻 $ (n+1) $ 的源项未知，必须作近似。为此，在当前时间层 n 附近将源项线性化，得到：

$$ \vec{\boldsymbol{Q}}^{n+1}\approx\vec{\boldsymbol{Q}}^{n}+\frac{\partial\vec{\boldsymbol{Q}}}{\partial\vec{\boldsymbol{W}}}\Delta\vec{\boldsymbol{W}}^{n}.\tag{6.10} $$

将式（6.10）代入式（6.9）并整理，得 $ [9, 10] $：

$$ \left[\frac{1}{\Delta t_{I}}\bar{I}-\left(\frac{\partial\vec{Q}}{\partial\vec{W}}\right)_{I}\right]\Delta\vec{W}_{I}^{n}=-\frac{1}{\Omega_{I}}\left[(\vec{R}_{Q})_{I}^{n}-\Omega_{I}\vec{Q}_{I}^{n}\right],\tag{6.11} $$

其中，$ \bar{I} $ 为单位矩阵。式（6.11）称为点隐式形式，因为左端方括号内的隐式算子只依赖控制体 $ \Omega_{I} $ 自身的数值。与式（6.9）比较可见，原来的标量时间步长 $ \Delta t $ 现在变成一个矩阵。因此，每个流动方程按与其特征值对应的独立参数缩放，从而补偿不同时间尺度间的差异，缓解源项引起的时间步长限制。

将上述点隐式方法用于式（6.5）的多级格式，第 k 级为：

$$ \vec{W}_{I}^{(k)}=\vec{W}_{I}^{(0)}-\left[(\vec{R}_{Q})_{I}^{(k-1)}-\Omega_{I}\vec{Q}_{I}^{(0)}\right]\left[\frac{\Omega_{I}}{\alpha_{k}\Delta t_{I}}\bar{I}-\left(\frac{\partial\vec{Q}}{\partial\vec{W}}\right)_{I}\right]^{-1}\tag{6.12} $$

校注：原书把残差列矢量写在逆矩阵左侧，且该矩阵中的源项导数未乘控制体积。按式（6.11），应核对逆矩阵左乘残差，以及源项导数的体积因子；译文保留原书写法。

其中，$ \vec{R}_{Q} $ 定义见式（6.9）。式（6.7）的混合多级格式也有类似表达式。源项对稳定性影响的详细研究见文献 [11–13]。

更复杂的方法是采用隐式 Runge–Kutta 格式处理刚性源项，最早由 Rosenbrock $ [14] $ 和 Butcher $ [15] $ 描述。这类格式比式（6.12）的方法具有更好的数值稳定性，但每一级都要对大型线性方程组求逆。隐式 Runge–Kutta 格式及其应用见 6.2.6 节。

6.1.4 最大时间步长的确定

每种显式时间推进格式，都只有在时间步长 $ \Delta t $ 不超过某一值时才稳定。要保持稳定，格式必须满足所谓 CFL 条件 $ [16] $，即数值方法的依赖域必须包含偏微分方程的依赖域。对基本显式格式 $ (6.4) $，这意味着时间步长不得大于信息穿过空间离散模板所需的时间。因此，对一维线性对流方程，条件为：

$$ \Delta t=\sigma\frac{\Delta x}{\left|\Lambda_{c}\right|},\tag{6.13} $$

其中，$ \Delta x/|\Lambda_{c}| $ 为信息穿过尺寸为 $ \Delta x $ 的单元所需的时间，传播速度为 $ \Lambda_{c} $。速度 $ \Lambda_{c} $ 对应于对流通量 Jacobian 的最大特征值。正系数 $ \sigma $ 为 CFL 数，其大小取决于时间推进格式的类型和参数，以及空间离散的形式。10.3 节将用两个模型问题考察 $ \sigma $ 的依赖关系。表 6.1 和表 6.2 列出各种多级格式及离散方法的 CFL 数。

对线性模型方程，可借助 von Neumann 稳定性分析确定最大时间步长（10.3 节）。不过，对多维和非线性控制方程，只能近似计算。以下介绍结构及非结构网格上，无黏与黏性流动的时间步长估算公式。

结构网格上的时间步长

Euler 方程

在结构网格上，时间步长 $ \Delta t $ 对控制体 $ \Omega_{I} $ 可由以下近似关系确定 [17–19]：

$$ \Delta t_{I}=\sigma\frac{\Omega_{I}}{(\hat{\Lambda}_{c}^{I}+\hat{\Lambda}_{c}^{J}+\hat{\Lambda}_{c}^{K})_{I}}.\tag{6.14} $$

多级格式的 CFL 数 $ \sigma $ 见表 6.1，混合格式见表 6.2。三个网格方向的对流通量 Jacobian（A.9）的谱半径为：

$$ \begin{aligned}\hat{\boldsymbol{\Lambda}}_{c}^{I}&=\left(|\vec{v}\cdot\vec{n}^{I}|+c\right)\Delta S^{I},\\\hat{\boldsymbol{\Lambda}}_{c}^{J}&=\left(|\vec{v}\cdot\vec{n}^{J}|+c\right)\Delta S^{J},\\\hat{\boldsymbol{\Lambda}}_{c}^{K}&=\left(|\vec{v}\cdot\vec{n}^{K}|+c\right)\Delta S^{K}.\end{aligned}\tag{6.15} $$

式（6.15）的法向量与面面积，由控制体在相应方向上两个相对面的对应数值平均得到。例如，若对偶控制体（4.2.3 节）的方向如图 4.1(b) 所示，则 I 方向取：

$$ \vec{n}_{i,j,k}^{I}=\frac{1}{2}(\vec{n}_{1}-\vec{n}_{2}),\quad\Delta S_{i,j,k}^{I}=\frac{1}{2}(\Delta S_{1}+\Delta S_{2}).\tag{6.16} $$

J、K 方向也有类似表达式。

式（6.14）给出的是只对某一控制体有效的局部时间步长。若只求定常解，可用局部时间步长加快收敛，参见 9.1 节；若时间精度重要，则所有控制体必须采用统一的全局时间步长，即：

$$ \Delta t=\min_{I}\left(\Delta t_{I}\right),\tag{6.17} $$

其中，最小值取遍所有控制体。

Navier–Stokes 方程

对黏性流动，计算 $ \Delta t $ 时还须包含黏性通量 Jacobian（A.10）的谱半径。它们会在边界层内严重限制最大时间步长。时间步长可按下式计算 [5, 18, 19]：

$$ \Delta t_{I}=\sigma\frac{\Omega_{I}}{(\hat{\Lambda}_{\mathrm{c}}^{I}+\hat{\Lambda}_{\mathrm{c}}^{J}+\hat{\Lambda}_{\mathrm{c}}^{K})_{I}+C(\hat{\Lambda}_{\mathrm{v}}^{I}+\hat{\Lambda}_{\mathrm{v}}^{J}+\hat{\Lambda}_{\mathrm{v}}^{K})_{I}}.\tag{6.18} $$

乘在黏性谱半径前的常数，中心空间离散通常取 C=4，一阶迎风离散取 C=2，二阶迎风离散取 C=1。若采用涡黏湍流模型，则黏性谱半径为 $ [18, 19] $：

$$ \hat{\Lambda}_{\mathrm{v}}^{I}=\max\left(\frac{4}{3\rho},\frac{\gamma}{\rho}\right)\left(\frac{\mu_{\mathrm{L}}}{\mathrm{Pr}_{\mathrm{L}}}+\frac{\mu_{\mathrm{T}}}{\mathrm{Pr}_{\mathrm{T}}}\right)\frac{(\Delta S^{I})^{2}}{\Omega}\tag{6.19} $$

其他方向类似。式（6.19）中，$ \mu_{L} $ 与 $ \mu_{T} $ 分别为层流及湍流动力黏度系数，$ \mathrm{Pr}_{L} $ 与 $ \mathrm{Pr}_{T} $ 分别为层流及湍流 Prandtl 数。表 6.1、6.2 的 CFL 数同样适用于黏性流动。式（6.7）的 (5,3) 混合格式对黏性流动尤其高效。

非结构网格上的时间步长

已有若干方法估算非结构网格上的最大时间步长，下面介绍两种不同方法。

方法 1

一种经过验证、与结构网格实现十分接近的方法为 $ [6] $：

$$ \Delta t_{I}=\sigma\frac{\Omega_{I}}{(\hat{\Lambda}_{c}+C\hat{\Lambda}_{\mathrm{v}})_{I}},\tag{6.20} $$

其中，$ \hat{\Lambda}_{c} $ 和 $ \hat{\Lambda}_{v} $ 分别为控制体全部面上的对流谱半径之和与黏性谱半径之和。与结构网格一样，通常取 $ 1 \leq C \leq 4 $。对单元中心格式，谱半径定义为 [6]：

$$ \begin{aligned}&(\hat{\Lambda}_{\mathrm{c}})_{I}=\sum_{J=1}^{N_{F}}\left(|\vec{v}_{I J}\cdot\vec{n}_{I J}|+c_{I J}\right)\Delta S_{I J},\\&(\hat{\Lambda}_{\mathrm{v}})_{I}=\frac{1}{\Omega_{I}}\sum_{J=1}^{N_{F}}\left[\max\left(\frac{4}{3\rho_{I J}},\frac{\gamma_{I J}}{\rho_{I J}}\right)\left(\frac{\mu_{\mathrm{L}}}{\mathrm{Pr}_{\mathrm{L}}}+\frac{\mu_{\mathrm{T}}}{\mathrm{Pr}_{\mathrm{T}}}\right)_{I J}\left(\Delta S_{I J}\right)^{2}\right].\\ \end{aligned}\tag{6.21} $$

控制体面上的流动变量由算术平均得到。

方法 2

式（6.21）预测的谱半径过大，混合单元网格尤其如此，因而时间步长小于稳定性实际需要的值。文献 [20] 的实现可更准确地估算时间步长，即：

$$ \Delta t_{I}=\sigma\frac{\Omega_{I}}{(\hat{\Lambda}_{c}^{x}+\hat{\Lambda}_{c}^{y}+\hat{\Lambda}_{c}^{z})_{I}+C(\hat{\Lambda}_{\mathrm{v}}^{x}+\hat{\Lambda}_{\mathrm{v}}^{y}+\hat{\Lambda}_{\mathrm{v}}^{z})_{I}}\tag{6.22} $$

其中，对流谱半径为：

$$ \begin{array}{l}\hat{\Lambda}_{c}^{x}=(|u|+c)\Delta\hat{S}^{x},\\\hat{\Lambda}_{c}^{y}=(| v|+c)\Delta\hat{S}^{y},\\\hat{\Lambda}_{c}^{z}=(|w|+c)\Delta\hat{S}^{z}\end{array}\tag{6.23} $$

黏性谱半径为，假设采用涡黏湍流模型：

$$ \hat{\Lambda}_{\mathrm{v}}^{x}=\max\left(\frac{4}{3\rho},\frac{\gamma}{\rho}\right)\left(\frac{\mu_{\mathrm{L}}}{\mathrm{Pr}_{\mathrm{L}}}+\frac{\mu_{\mathrm{T}}}{\mathrm{Pr}_{\mathrm{T}}}\right)\frac{(\Delta\hat{S}^{x})^{2}}{\Omega},\text{等}\tag{6.24} $$

变量 $ \Delta\hat{S}^{x} $、$ \Delta\hat{S}^{y} $、$ \Delta\hat{S}^{z} $ 分别为控制体在 y–z、x–z、x–y 平面上的投影，按下式计算：

$$ \begin{aligned}\Delta\hat{S}^{x}&=\frac{1}{2}\sum_{J=1}^{N_{F}}|S_{x}|_{J},\\\Delta\hat{S}^{y}&=\frac{1}{2}\sum_{J=1}^{N_{F}}|S_{y}|_{J},\\\Delta\hat{S}^{z}&=\frac{1}{2}\sum_{J=1}^{N_{F}}|S_{z}|_{J},\end{aligned}\tag{6.25} $$

其中，$ S_{x} $、$ S_{y} $、$ S_{z} $ 分别为面矢量的 x、$ \gamma- $、z 分量；该面矢量为 $ \vec{S} = \vec{n} \cdot \Delta S $。

CFL 数一般与结构网格相同，因此仍可采用表 6.1、6.2 的值。同样，可用局部时间推进加快定常收敛。模拟非定常流动需要的全局时间步长，仍可按式（6.17）确定。