6.3 非定常流动的求解方法

6.3 非定常流动的求解方法

非定常流动模拟如今在许多工程领域十分重要，例如叶轮机械静止与旋转部件的相互作用、活塞发动机、流固耦合、直升机空气动力学、气动声学、湍流的 DNS 或 LES、爆轰等。显然，模拟既要高效，精度也要满足问题需求。

对某些非定常应用，当时间尺度与空间尺度除以特征值相当，即物理过程要求的 CFL 数为 1 的量级时，显式格式是最佳选择，例如气动声学、DNS、LES。显式 Runge–Kutta 格式在这类应用中十分常用。由于整体物理现象演化远慢于局部解的变化，需要在很长的物理时间区间积分。为准确实现这一点，显式格式必须具有三阶及更高阶的时间分辨精度，因而需采用不同于 6.1 节所介绍的 Runge–Kutta 方法，例如见文献 [122, 123]。

另一些情况下，物理时间尺度远大于空间尺度除以特征值，例如颤振及转静子相互作用，CFL 数可取数百甚至数千，而不损害模拟精度。显然，此时隐式格式更适合。下面讨论常用于非定常流动的双时间推进方法，其与 6.2.6 节隐式 Runge–Kutta 格式的比较见文献 $ [124-126] $。

双时间推进以式（6.2）基本非线性格式的二阶时间精度版本为基础。为此，在式（6.2）取 $ \beta = 1 $、$ \omega = 1/2 $，得到：

$$ \frac{3(\Omega\bar{M})_{I}^{n+1}\vec{W}_{I}^{n+1}-4(\Omega\bar{M})_{I}^{n}\vec{W}_{I}^{n}+(\Omega\bar{M})_{I}^{n-1}\vec{W}_{I}^{n-1}}{2\Delta t}=-\vec{R}_{I}^{n+1},\tag{6.79} $$

其中，$ \Delta t $ 为全局物理时间步长，$ \bar{M} $ 为质量矩阵。式（6.79）对式（6.1）的时间导数作三点后向差分（后向 Euler）近似。求解这一非线性方程组，可采用 Newton 方法或时间推进方法。后者为：

$$ \frac{\partial}{\partial t^{*}}\left(\Omega_{I}^{n+1}\vec{W}_{I}^{*}\right)=-\vec{R}_{I}^{*}(\vec{W}^{*}),\tag{6.80} $$

其中，$ \vec{W}^{*} $ 为 $ \vec{W}^{n+1} $ 的近似，$ t^{*} $ 为伪时间变量。注意，时间导数中没有质量矩阵。非定常残差定义为：

$$ \vec{R}_{I}^{*}(\vec{W}^{*})=\vec{R}_{I}(\vec{W}^{*})+\frac{3}{2\Delta t}(\Omega\bar{M})_{I}^{n+1}\vec{W}_{I}^{*}-\vec{Q}_{I}^{*}.\tag{6.81} $$

将式（6.80）时间推进过程中保持不变的全部项归入源项，即：

$$ \vec{Q}_{I}^{*}=\frac{2}{\Delta t}(\Omega\bar{M})_{I}^{n}\vec{W}_{I}^{n}-\frac{1}{2\Delta t}(\Omega\bar{M})_{I}^{n-1}\vec{W}_{I}^{n-1}.\tag{6.82} $$

网格运动和／或变形时，式（6.80）中新控制体的体积 $ \Omega^{n+1} $ 必须满足几何守恒律，详细说明及文献见 A.5 节。

式（6.80）的定常解对应于新时间层的流动变量，即 $ \vec{W}^{*} = \vec{W}^{n+1} $。伪时间中的定常状态满足 $ \vec{R}_{I}^{*} = 0 $，因此式（6.79）成立。前述任一显式或隐式时间推进格式，都可在伪时间中求解式（6.80）。下面分别讨论显式多级与隐式格式中的双时间推进实现。

6.3.1 显式多级格式的双时间推进

Jameson [127] 首先用显式多级格式实现双时间方法，并采用局部时间推进与多重网格加速。其主要优点，是物理时间步长不受通常显式方法的限制，可仅根据流动物理选择。额外存储也只用于源项 $ \vec{Q}^{*} $，因此很有吸引力。求解伪时间问题（6.80）的 m 级显式格式为：

$$ \begin{aligned}\vec{W}_{I}^{(0)}&=(\vec{W}_{I}^{*})^{l},\\\vec{W}_{I}^{(1)}&=\vec{W}_{I}^{(0)}-\frac{\alpha_{1}\Delta t_{I}^{*}}{\Omega_{I}^{n+1}}\vec{R}_{I}^{*}(\vec{W}^{(0)}),\\\vec{W}_{I}^{(2)}&=\vec{W}_{I}^{(0)}-\frac{\alpha_{2}\Delta t_{I}^{*}}{\Omega_{I}^{n+1}}\vec{R}_{I}^{*}(\vec{W}^{(1)}),\\&\vdots\\(\vec W_I^*)^{l+1}&=\vec W_I^{(0)}-\frac{\alpha_m\Delta t_I^*}{\Omega_I^{n+1}}\vec R_I^*(\vec W^{(m-1)}).\end{aligned}\tag{6.83} $$

其中，l 为当前伪时间层，$ (l+1) $ 为新伪时间层。时间推进的初值可取 $ (\vec{W}^{*})^{l} = \vec{W}^{n} $，也可由之前的物理时间步外推，例如文献 [128] 的取法：

$$ (\vec{W}_{I}^{*})^{l}=\vec{W}^{n}+\frac{3\vec{W}^{n}-4\vec{W}^{n-1}+\vec{W}^{n-2}}{2}.\tag{6.84} $$

持续推进，直到 $ (\vec{W}_{I}^{*})^{l+1} $ 以足够精度逼近 $ \vec{W}_{I}^{n+1} $；通常以残差 $ \vec{R}_{I}^{*} $ 下降两个或三个数量级为准。之后再进行下一物理时间步。伪时间步长 $ \Delta t^{*} $ 按 6.1.4 节的方法计算。

Arnone 等 [128] 指出，物理时间步长 $ \Delta t $ 与伪时间步长 $ \Delta t^{*} $ 相当或更小时，多级格式（6.83）会失稳。Melson 等 [129] 证明，失稳由式（6.81）的以下项引起：

$$ \frac{3}{2\Delta t}(\Omega\bar{M})_{I}^{n+1}\vec{W}_{I}^{*} $$

这一项在 $ \Delta t $ 较小时变得显著。他们建议隐式处理这一项。因此，需修改式（6.83）的多级格式，使第 k 级成为 [129]：

$$ \begin{aligned}&\vec{W}_{I}^{(k)}=\vec{W}_{I}^{(0)}-\frac{\alpha_{k}\Delta t_{I}^{*}}{\Omega_{I}^{n+1}}\left[\bar{I}+\frac{3}{2\Delta t}\alpha_{k}\Delta t_{I}^{*}\bar{M}^{n+1}\right]^{-1}\\&\quad\cdot\left[\vec{R}_{I}(\vec{W}^{(k-1)})+\frac{3}{2\Delta t}(\Omega\bar{M})_{I}^{n+1}\vec{W}_{I}^{(0)}-\vec{Q}_{I}^{*}\right].\end{aligned}\tag{6.85} $$

同样的方法也可用于混合多级格式（见 6.1.2 节）。上述形式（6.85）对任意物理时间步长 $ \Delta t $ 都稳定 [129]。

对单元中心格式，可将式（6.85）的质量矩阵 $ \bar{M}^{n+1} $ 集中处理，即以单位矩阵替代，而不降低解精度。这样，以下项变成标量：

$$ \left[\bar{I}+\frac{3}{2\Delta t}\alpha_{k}\Delta t_{I}^{*}\bar{M}^{n+1}\right]^{-1} $$

不过，对单元顶点空间离散必须考虑质量矩阵，否则小物理时间步长下多级格式会失稳。为避免代价高昂的 $ \bar{M}^{n+1} $ 求逆，Venkatakrishnan [130] 以及 Venkatakrishnan 和 Mavriplis [131] 建议将式（6.85）修改为：

$$ \begin{aligned}&\vec{W}_{I}^{(k)}=\vec{W}_{I}^{(0)}-\frac{\alpha_{k}\Delta t_{I}^{*}}{\Omega_{I}^{n+1}}\left[1+\frac{3}{2\Delta t}\alpha_{k}\Delta t_{I}^{*}\beta\right]^{-1}\\&\quad\cdot\left[\vec{R}_{I}^{*}(\vec{W}^{(k-1)})-\frac{3}{2\Delta t}\Omega_{I}^{n+1}\beta\vec{W}_{I}^{(k-1)}\right].\end{aligned}\tag{6.86} $$

校注：原书修正格式的括号中减去了含上一中间解的项，但零级解位于逆因子之外。按式（6.80）、（6.81）的固定点要求，定常伪时间解应满足非定常残差为零；此式应核对零级解相关的补偿项。译文保留原式。

参数 $ \beta $ 可用于稳定时间推进格式。实际经验表明，取 $ \beta = 2 $ 已足够 [130, 131]。

在伪时间中用显式多级格式求解的双时间方法应用广泛。以局部时间推进（针对 $ t^{*} $）和多重网格加速时，计算效率最高，因为多重网格只需少数循环就能迅速收敛到式（6.80）的定常解。结构网格应用实例见文献 [127–129]、[132–134]；非结构网格实现例如见文献 [130, 131, 135]。

6.3.2 隐式格式的双时间推进

在伪时间 $ t^{*} $ 中采用隐式格式求解式（6.80），步骤与 6.2 节相同。首先，将式（6.80）写成非线性隐式格式：

$$ \frac{\partial}{\partial t^{*}}\left(\Omega_{I}^{n+1}\vec{W}_{I}^{*}\right)=-(\vec{R}_{I}^{*})^{l+1}\tag{6.87} $$

其中，$ (l+1) $ 为新伪时间层。再次注意，时间导数中没有 $ \bar{M} $。式（6.81）定义的非定常残差，可在伪时间中线性化为：

$$ (\vec{R}^{*})^{l+1}\approx(\vec{R}^{*})^{l}+\frac{\partial\vec{R}^{*}}{\partial\vec{W}^{*}}\Delta\vec{W}^{*},\tag{6.88} $$

其中，$ \Delta\vec{W}^{*}=(\vec{W}^{*})^{l+1}-(\vec{W}^{*})^{l} $；通量 Jacobian 定义为：

$$ \frac{\partial\vec{R}^{*}}{\partial\vec{W}^{*}}=\frac{\partial\vec{R}}{\partial\vec{W}}+\frac{3}{2\Delta t}(\Omega\bar{M})^{n+1}.\tag{6.89} $$

将上述线性化代入式（6.87），得到未分解的隐式格式 $ [136] $：

$$ \left[\left(\frac{1}{\Delta t_{I}^{*}}+\frac{3}{2\Delta t}\right)(\Omega\bar{M})_{I}^{n+1}+\left(\frac{\partial\vec{R}}{\partial\vec{W}}\right)_{I}\right]\Delta\vec{W}^{*}=-(\vec{R}_{I}^{*})^{l}.\tag{6.90} $$

校注：式（6.80）、（6.87）的伪时间导数明确不含质量矩阵，但原书式（6.90）在伪时间项中又乘了质量矩阵。按前式直接线性化，伪时间贡献应为 $\Omega_I^{n+1}\bar I/\Delta t_I^*$，物理时间贡献才包含质量矩阵。译文保留原式并注明这一不一致。

6.2 节的任一方法均可求解方程组（6.90）。时间精确隐式方法的详细讨论见文献 $ [137] $，具体实现例如见文献 $ [136, 138, 139] $。