7.3 大涡模拟
7.3 大涡模拟
Smagorinsky 早在 1963 年就在气象学（大气环流）中采用了所谓的 LES 方法 $ [98] $。Deardorff 于 1970 年给出了 LES 的首个工程应用，即湍流槽道流 $ [99] $。随后 Schumann 对该方法作了扩展和改进 $ [100] $。20 世纪 80 年代，湍流模拟的研究重点从 LES 转向 DNS，但仍开展了一些重要的 LES 研究，例如 Bardina 等人 $ [101] $ 以及 Moin 和 Kim $ [102] $ 的工作。20 世纪 90 年代初，对 LES 的关注再度上升 $ [103–109] $。如今，LES 越来越多地用于燃烧室等在物理过程和几何形状上都较复杂的工程流动，实例见文献 $ [110, 110–128] $。这一趋势无疑受到低成本、高性能计算服务器普及的推动。此外，今天的工程师常常需要处理标准湍流模型无法解决的流动问题。在某些情况下，平均流动的频率与湍流脉动频率处于同一数量级，此时时间平均不再有效，必须采用 LES 或 DNS。
LES 基于以下观察：小尺度湍流结构比大涡具有更强的普适性。因此，其思路是计算携带能量的大尺度结构对动量和能量输运的贡献，同时对数值格式无法解析的小尺度结构的影响进行建模。由于小尺度结构更加均匀且更具普适性，所谓的亚格子尺度模型可以比 RANS 方程所用的湍流模型简单得多。
LES 对控制方程进行三维非定常求解。与基于 RANS 方程的湍流建模相比，LES 在流向 $ (50 \leq x^{+} \leq 150) $ 和横向 $ (15 \leq z^{+} \leq 40) $ 也要求较高的网格分辨率。不过，LES 的计算成本远低于 DNS。解析外层所需的网格点（单元）数与 $ Re^{0.4} $ 成正比 [129]；在黏性底层，分辨率则必须按 $ Re^{1.8} $ 增加。因此，与 DNS 所需的 $ Re^{9/4} $ 相比，LES 可用于至少高一个数量级的雷诺数。为了进一步降低网格分辨率要求，可以将 LES 与近似壁面模型结合（7.3.4 节），或与 RANS 模型耦合（7.3.5 节）。两种方法均可在合理的计算成本下实现工程问题的 LES。
为了准确解析高波数的湍流脉动，空间离散格式必须在波数空间中具有相应的性能（例如参见 [130] 或 [131]）。因此，人们常采用谱方法。不过，谱方法仅适用于具有（准）周期边界且几何形状简单的区域，因而有限差分和有限体积空间离散方法日益普及。中心差分格式已被证明比迎风格式更适合。这是因为迎风格式无论精度阶数如何，都会因其固有的数值阻尼，在湍流频谱的相当大范围内耗散过多能量 [132,133]。不同离散方法的数值误差讨论见文献 [134–136]。
在非结构网格上实现 LES 方法 $ [137-147] $ 尤具挑战，但这使得处理极其复杂的几何形状、移动边界或动态网格自适应成为可能。研究重点是为混合单元网格发展数值效率高的高阶空间离散方法。
LES 的入门介绍见文献 [25，第 269—336 页]、[148–152] 或 [113]，研究现状综述见文献 [153–155]。
7.3.1 空间滤波
LES 基于空间滤波操作，将任意流动变量 U 分解为滤波后的部分 $ \overline{U} $（大尺度、已解析）与滤波尺度以下的部分 $ U' $（未解析），即
$$ U=\overline{U}+U^{\prime}\tag{7.76} $$
空间位置 $ \vec{r}_{0} $ 处的滤波变量定义为
$$ \overline{{U}}(\vec{r}_{0},t)=\int_{D}U(\vec{r},t)\;G(\vec{r}_{0},\vec{r},\Delta)\mathrm{d}\vec{r},\tag{7.77} $$
其中，$ \Omega $ 表示整个流动区域，G 为滤波函数，$ \vec{r} $ 为位置矢量。滤波函数决定小尺度结构的形式和大小，其取决于差值 $ \vec{r}_{0}-\vec{r} $ 及滤波宽度 $ \Delta=(\Delta_{1}\Delta_{2}\Delta_{3})^{1/3} $，其中 $ \Delta_{i} $ 为第 i 个空间坐标方向上的滤波宽度。最常用的滤波函数如下（见图 7.3）：
#### 1. 顶帽滤波器
$$ G=\left\{\begin{aligned}&1/\Delta^{3}&\text{若 }|(x_{0})_{i}-x_{i}|\leq\Delta_{i}/2\\ &0&\text{其他情况。}\end{aligned}\right.\tag{7.78} $$
2. 锐 Fourier 截断滤波器
$$ G=\prod_{i=1}^{3}\frac{\sin\left(\frac{\pi}{\Delta_{i}}[(x_{0})_{i}-x_{i}]\right)}{\pi\left[(x_{0})_{i}-x_{i}\right]}.\tag{7.79} $$
3. Gaussian 滤波器
$$ G=\left(\frac{6}{\pi\Delta^{2}}\right)^{3/2}\exp\left(\frac{-6\left\|\vec{r}_{0}-\vec{r}\right\|_{2}^{2}}{\Delta^{2}}\right).\tag{7.80} $$
顶帽滤波器和 Gaussian 滤波器不仅平滑滤波宽度以下的小尺度，也会平滑大尺度脉动。截断滤波器只影响截断波数以下所对应的小尺度。在实际应用中，Gaussian 滤波器总是与锐 Fourier 截断结合使用。适用于单元尺寸变化网格的滤波器见文献 $ ^{[156, 157]} $。
7.3.2 滤波后的控制方程
必须将式（7.76）、（7.77）定义的空间滤波作用于 Navier–Stokes 方程，以去除小尺度湍流。滤波宽度 $ \Delta $ 和滤波函数均视为自由参数。实际上，通常不对控制方程进行显式滤波，而是假定网格及离散误差共同定义滤波器 G。显式滤波的讨论见文献 $ [158, 159] $。
不可压缩 Navier–Stokes 方程
对于 Newton 流体的不可压缩流动，滤波后的控制方程（7.6）为
$$ \begin{aligned}\frac{\partial\overline{ v}_{i}}{\partial x_{i}}&=0,\\\frac{\partial\overline{ v}_{i}}{\partial t}+\frac{\partial}{\partial x_{j}}(\overline{ v}_{i}\overline{ v}_{j})&=-\frac{1}{\rho}\frac{\partial\overline{p}}{\partial x_{i}}+\nu\nabla^{2}\overline{ v}_{i}-\frac{\partial\tau_{ij}^{S}}{\partial x_{j}},\end{aligned}\tag{7.81} $$
其中，$ \nu $ 表示运动黏度。式（7.81）描述携带能量的大尺度运动在时间和空间中的演化。对流项的非线性使方程中出现所谓的亚格子尺度（SGS）应力张量：
$$ \tau_{i j}^{S}=\overline{{ v_{i} v_{j}}}-\overline{{ v}}_{i}\overline{{ v}}_{j},\tag{7.82} $$
它描述未解析尺度的影响。为了封闭方程，必须对 SGS 张量建模（见 7.3.3 节）。
SGS 张量可以分解为三个部分 $ [160] $，即
$$ \tau_{ij}^{S}=L_{ij}+C_{ij}+\tau_{ij}^{SR}.\tag{7.83} $$
各部分具有以下物理意义：
$$ L_{i j}=\overline{{\bar{ v}_{i}\bar{ v}_{j}}}-\bar{ v}_{i}\bar{ v}_{j}\tag{7.84} $$
称为 Leonard 应力项，表示大尺度涡之间产生小尺度湍流的相互作用。只有这一项能够由滤波后的速度场 $ v_{i} $ 显式计算。此外，交叉应力项
$$ C_{i j}=\overline{{\bar{ v}_{i} v_{j}^{\prime}}}+\overline{{ v_{i}^{\prime}\overline{{ v}}_{j}}}\tag{7.85} $$
描述大尺度涡与小尺度涡之间的相互作用。最后，
$$ \tau_{ij}^{SR}=\overline{ v_{i}^{\prime} v_{j}^{\prime}}\tag{7.86} $$
称为 SGS Reynolds 应力张量，反映小尺度结构之间的相互作用。上述分解（7.83）已不再使用，主要原因是 $ L_{ij} $ 和 $ C_{ij} $ 在 Galilean 变换下不具有不变性。1
1 Galilean 不变性：彼此作匀速平移的所有参考系均等价。返回正文
可压缩 Navier–Stokes 方程
将 LES 用于可压缩流动时，必须对式（7.1）同时采用 Favre 平均（7.1.2 节）和空间滤波，否则滤波后的 Navier–Stokes 方程会包含密度与速度、温度等其他变量的乘积。因此，式（7.1）中的速度分量、能量和温度分解为
$$ U=\tilde{U}+U^{\prime\prime}.\tag{7.87} $$
空间位置 $ \vec{r}_{0} $ 处的滤波变量为
$$ \tilde{U}(\vec{r}_{0},t)=\frac{\overline{\rho U}}{\overline{\rho}}=\frac{1}{\overline{\rho}}\int_{D}\rho(\vec{r},t)U(\vec{r},t)G(\vec{r}_{0},\vec{r},\Delta)\mathrm{d}\vec{r},\tag{7.88} $$
其中，上横线表示按式（7.77）进行的滤波。Favre 滤波 Navier–Stokes 方程（7.1）为 $ [149, 153] $
$$ \begin{gathered}\frac{\partial\overline{{\rho}}}{\partial t}+\frac{\partial}{\partial x_{j}}(\overline{{\rho}}\tilde{ v}_{j})=0,\\\frac{\partial\overline{\rho}\tilde{ v}_{i}}{\partial t}+\frac{\partial(\overline{\rho}\tilde{ v}_{j}\tilde{ v}_{i})}{\partial x_{j}}+\frac{\partial\overline{p}}{\partial x_{i}}-\frac{\partial\hat{\sigma}_{i j}}{\partial x_{j}}=-\frac{\partial\tau_{i j}^{S F}}{\partial x_{j}}+\frac{\partial}{\partial x_{j}}\left(\overline{\sigma}_{i j}-\hat{\sigma}_{i j}\right),\\\frac{\partial\overline{\rho\tilde{e}}}{\partial t}+\frac{\partial(\overline{\rho}\tilde{ v}_{j}\tilde{e})}{\partial x_{j}}+\frac{\partial\hat{q}}{\partial x_{j}}+\overline{p}\tilde{S}_{k k}-\hat{\sigma}_{i j}\tilde{S}_{i j}=-\mathcal{A}-\mathcal{B}-\mathcal{C}+\mathcal{D}\end{gathered}\tag{7.89} $$
其中，各项为
$$ \mathcal{A}=\frac{\partial}{\partial x_{j}}\left[\overline{\rho}(\widetilde{ v_{j}e}-\widetilde{ v_{j}}\widetilde{e})\right]\quad-\text{亚格子尺度热通量的散度} $$
$$ \mathcal{B}=\frac{\partial}{\partial x_{j}}\left[\overline{{q}}_{j}-\hat{q}_{j}\right]\quad-\text{SGS 热扩散的散度} $$
$$ \mathcal{C}\;=\;\left[\overline{{p S_{k k}}}-\overline{{p}}\tilde{S}_{k k}\right]\quad-\text{SGS 压力—膨胀} $$
$$ \mathcal{D}=\left[\overline{\sigma_{ij}S_{ij}}-\hat{\sigma}_{ij}\tilde{S}_{ij}\right]\quad-\text{SGS 黏性耗散} $$
以及
$$ \begin{gathered}\overline{\sigma}_{ij}=\overline{2\mu S_{ij}}+\overline{\left(\mu_{\mathrm{B}}-\frac{2\mu}{3}\right)\delta_{ij}S_{kk}},\\\hat{\sigma}_{i j}=2\tilde{\mu}\tilde{S}_{i j}+\left(\tilde{\mu}_{\mathrm{B}}-\frac{2\tilde{\mu}}{3}\right)\delta_{i j}\tilde{S}_{k k},\\\tilde{S}_{i j}=\frac{1}{2}\left(\frac{\partial\tilde{ v}_{i}}{\partial x_{j}}+\frac{\partial\tilde{ v}_{j}}{\partial x_{i}}\right),\\\overline{{q}}_{j}=-k\frac{\overline{{\partial T}}}{\partial x_{j}},\quad\tilde{q}_{j}=-\tilde{k}\frac{\partial\tilde{T}}{\partial x_{j}}.\end{gathered}\tag{7.90} $$
式（7.89）、（7.90）中，e 为单位质量内能，$ \tilde{S}_{ij} $ 是 Favre 滤波应变率张量，$\tau_{ij}^{SF}=\overline{\rho}(\widetilde{v_i v_j}-\tilde v_i\tilde v_j)$ 表示 Favre 平均亚格子尺度应力。此外，$ \mu $、$ \mu_{B} $ 和 k 分别表示分子黏度、体积黏度和热导率。最后，$ \tilde{\mu} $、$ \tilde{\mu}_{B} $、$ \tilde{k} $ 分别是滤波温度 $ \tilde{T} $ 下的对应值。
式（7.89）右端包含必须建模的项。在动量方程中，对 SGS 应力 $ \tau_{ij}^{SF} $ 作近似处理，而第二项，即 $ (\overline{\sigma}_{ij} - \hat{\sigma}_{ij}) $，通常被忽略。在能量方程中，A 项可以通过 SGS 应力表示 [161]，B 项可以忽略，C 项和 D 项可采用文献 [162] 提出的方法建模。
7.3.3 亚格子尺度建模
亚格子尺度模型的主要任务是模拟大尺度与亚格子尺度之间的能量输运。平均而言，能量从大尺度向小尺度传递，即湍流级串过程。因此，亚格子尺度模型必须提供适当的能量耗散机制。不过，在某些情况下，能量也从小尺度流向大尺度，这一过程称为反向散射。模型显然也应计及这一效应。反向散射模型的讨论例如见文献 $ ^{[163]} $。
过去已经提出了多种亚格子尺度模型，相关研究仍在继续。这些模型可分为两大类。第一类显式模拟 SGS 张量 $ \tau_{ij}^{S} $，其必要条件是空间离散格式产生的数值耗散必须远小于亚格子尺度耗散。多数显式 SGS 模型基于下面介绍的涡黏性概念。随后介绍作为所有亚格子尺度模型基础的 Smagorinsky 模型，并简要讨论所谓的动态亚格子尺度模型的基本原理。文献 [164] 最近比较了六种不同的显式亚格子尺度模型。
第二类通过适当离散对流通量，隐式模拟 SGS 应力，因此省去 $ \tau_{ij}^{S} $。这类模型称为单调积分大涡模拟（MILES）。Boris 等人 [165] 最先提出这一方法，Grinstein 和 Fureby [166,167] 则推动了其应用。这种情况下，要求数值耗散能够正确模拟亚格子尺度耗散。文献 [168,169] 的研究表明，这并不容易做到。尽管如此，MILES 已较为成功地应用于多种流动问题 [170–177]。
涡黏性模型
这些显式模型能够表示小尺度的整体耗散效应，却不能重现能量交换的局部细节。对于不可压缩流动，涡黏性模型将 SGS 应力与大尺度应变率 $ \overline{S}_{ij} $ 联系如下：
$$ \tau_{i j}^{S}-\frac{\delta_{i j}}{3}\tau_{k k}^{S}=-2\nu_{\mathrm{T}}\overline{{S}}_{i j}.\tag{7.91} $$
应变率 $ \overline{S}_{ij} $ 采用滤波后的速度分量，由式（7.3）计算。为降低计算成本，涡黏度 $ \nu_{T} $ 通常由代数关系计算。SGS 应力的各向同性部分 $ (\tau_{kk}^{S}) $ 可以并入滤波压力 [178]、另行建模 [179]，或直接忽略。
$$ \tau_{i j}^{S F}-\frac{\delta_{i j}}{3}\tau_{k k}^{S F}=-2\overline{\rho}\nu_{\mathrm{T}}\tilde{S}_{i j}+\left(\frac{2\overline{\rho}\nu_{\mathrm{T}}}{3}\right)\frac{\partial\tilde{ v}_{k}}{\partial x_{k}}\delta_{i j}.\tag{7.92} $$
应变率张量 $ \tilde{S}_{ij} $ 的各分量见式（7.90）。
Smagorinsky SGS 模型
Smagorinsky 模型 $ [98] $ 基于平衡假设，即小尺度结构完全且瞬时地耗散从大尺度接收的全部能量。该代数模型为
$$ \nu_{\mathrm{T}}=(C_{s}\Delta)^{2}|\overline{{S}}|,\tag{7.93} $$
其中，$ |\overline{S}|=(2\overline{S}_{ij}\overline{S}_{ij})^{1/2} $ 为应变率张量的大小，$ C_{s} $ 表示 Smagorinsky 常数。Lilly [180] 得到的理论值为 $ C_{s}=0.18 $，但 Smagorinsky 常数取决于流动类型。例如，在剪切流动中，$ C_{s} $ 必须减小到约 0.1。式（7.93）中的滤波宽度 $ \Delta $ 通常取平均网格尺寸的两倍，即 $ \Delta=2(\Delta x_{1}\Delta x_{2}\Delta x_{3})^{1/3} $。
为计及近壁区小尺度结构增长减弱的效应，必须减小涡黏度 $ \nu_{T} $。因此，采用 Van Driest 阻尼将 Smagorinsky 模型（式（7.93））修正为
$$ \nu_{\mathrm{T}}=\left[C_{s}\Delta\left(1-e^{- y^{+}/25}\right)\right]^{2}|\overline{S}|,\tag{7.94} $$
其中，$ y^{+} $ 表示无量纲壁面距离（壁面距离的计算例如见文献 [48]）。
Smagorinsky 模型计算成本低且易于实现，但存在几个严重缺点：
-
在具有平均剪切的层流区中过度耗散；
-
在近壁区及层流—湍流转捩处需要特殊处理；
-
参数 $ C_{s} $ 没有唯一的定义；
-
不模拟能量反向散射过程。
针对这些不足，人们提出了其他多种方法（例如见文献 $ ^{[149]} $），其中下面介绍的动态模型十分流行。
动态 SGS 模型
动态 SGS 模型采用与 Smagorinsky 模型（式（7.93））相同的关系，计算式（7.91）或（7.92）中的涡黏度 $ \nu_{T} $。区别在于，将事先调整的 Smagorinsky 常数替换为随空间和时间动态变化的参数。因此，
$$ \nu_{\mathrm{T}}=C_{d}(\vec{r},t)\;\Delta^{2}|\overline{S}|.\tag{7.95} $$
参数 $ C_{d} $ 根据最小湍流尺度的能量含量计算。为此，Germano 等人 [181] 提出采用第二个滤波器，即所谓的测试滤波器 $ \hat{\Delta} $。测试滤波器的宽度必须大于作用于控制方程的滤波器 $ \Delta $ 的宽度（通常取 $ \hat{\Delta}=2\Delta $）。将测试滤波器作用于已滤波的方程，得到所谓的亚测试尺度应力 $ \tau_{ij}^{ST} $。
$$ \tau_{i j}^{S T}=\widehat{\overline{{ v_{i} v_{j}}}}-\hat{\overline{{\nu}}}_{i}\hat{\overline{{\nu}}}_{j}.\tag{7.96} $$
亚测试尺度应力通过 Germano 恒等式 [182] 与 SGS 应力 $ \tau_{ij}^{S} $（式（7.83））关联：
$$ \hat{L}_{i j}=\tau_{i j}^{S T}-\hat{\tau}_{i j}^{S}=\widehat{\overline{{ v}}_{i}\overline{{ v}}_{j}}-\hat{\overline{{\nu}}}_{i}\hat{\overline{{\nu}}}_{j},\tag{7.97} $$
其中，$ \hat{L}_{ij} $ 表示与测试滤波器对应的 Leonard 应力，反映长度介于滤波宽度 $ \Delta $ 与测试滤波宽度 $ \hat{\Delta} $ 之间的尺度对 Reynolds 应力的贡献。
采用涡黏性方法（式（7.91））与式（7.95），表示式（7.97）中的亚测试尺度应力和 SGS 应力，得到
$$ \hat{L}_{i j}-\frac{\delta_{i j}}{3}L_{k k}=-2C_{d}M_{i j}\tag{7.98} $$
其中，
$$ M_{i j}=\hat{\Delta}^{2}|\hat{\overline{{S}}}|\hat{\overline{{S}}}_{i j}-\left[\Delta^{2}|\overline{{S}}|\overline{{S}}_{i j}\right]^{\wedge}.\tag{7.99} $$
记号 $[\,\cdot\,]^{\wedge}$ 表示对方括号内的整个项进行测试滤波。参数 $ C_{d} $ 可通过 Lilly [183] 的最小二乘最小化方法，由式（7.98）求得，结果为
$$ C_{d}(\vec{r},t)=-\frac{1}{2}\frac{L_{i j}M_{i j}}{M_{m n}M_{m n}}.\tag{7.100} $$
上述表达式（7.100）在数学上并不一致，因为式（7.98）中将参数 $ C_{d} $ 移到了测试滤波器之外。因此，在实际应用中，对式（7.100）的分子和分母沿均匀方向作集合平均，即
$$ C_{d}(\vec{r},t)=-\frac{1}{2}\frac{\langle L_{i j}M_{i j}\rangle}{\langle M_{m n}M_{m n}\rangle}.\tag{7.101} $$
改进的动态 SGS 模型例如由 Ghosal 等人 $ [184] $、Carati 等人 $ [185] $、Piomelli 和 Liu $ [186] $，以及 Held $ [113] $ 提出。
7.3.4 壁面模型
对于高雷诺数 $ (Re > 10^{6}) $ 的壁面受限流动，LES 的计算成本仍然高于工程应用所能承受的水平。这是因为，为了适当地解析壁面层，需要过多的网格点（单元）。为了降低成本，可以通过指定外部流动速度与壁面应力之间的关联，对壁面层建模。这与 RANS 模拟中采用壁面函数十分相似。其基本假设是近壁区与外部区域之间的相互作用较弱，文献 [187,188] 的研究支持这一假设。
早期壁面模型的实现基于以下假设：壁面层的动力学具有普适性，因此可以用广义壁面律近似。基本上，这些模型采用对数律（见文献 $ [100, 189–191] $）。Balaras 等人 $ [192] $ 最近提出了一种新的分区方法。在该双层模型中，滤波后的 Navier–Stokes 方程（7.81）求解到壁面上方的第一个网格点；从该点到壁面的区域，则在加密的嵌入网格上求解二维边界层方程。随后利用嵌入网格上的解，为 LES 指定作为边界条件的壁面剪切应力。Balaras 等人 $ [192] $ 的分区方法允许将第一个点置于 $ 20 < y^{+} < 100 $ 区域，显著减少网格规模及计算时间。该方法已成功用于平面槽道、方管和旋转槽道中的湍流，后来也用于分离流动的 LES，并取得令人鼓舞的结果 $ [193–196] $。各类壁面模型的综述与评估见文献 $ [197] $。
7.3.5 分离涡模拟
虽然上述壁面模型有助于大幅减少网格点（单元）数，但对于复杂工程构型，LES 的成本仍然过高。因此，Spalart 最近提出了另一种方法，即分离涡模拟（DES），旨在模拟高雷诺数的大规模分离流动 $ [198, 199] $。该方法结合 RANS 与 LES：采用高度拉伸的网格及 RANS 湍流模型（主要为 7.2.1 节的 Spalart–Allmaras 模型或 7.2.3 节的 Menter SST 模型）解析附着边界层；在壁面区域以外，则采用各向同性网格进行 LES，以捕捉脱离壁面的三维涡。因此，DES 试图在同一框架中结合两种方法的优势。
在 Spalart–Allmaras 湍流模型中实现 DES 时，将式（7.36）、（7.38）—（7.40）中的壁面距离 d 替换为 DES 长度尺度：
$$ l=\min(d,\;C_{DES}\;\Delta).\tag{7.102} $$
该长度尺度取决于控制体积的最大尺寸，即 $ \Delta = \max(\Delta x, \Delta y, \Delta z) $。常数 $ C_{DES} $ 在一定程度上取决于流动类型。对于均匀湍流，$ C_{DES} = 0.65 $ 被认为是最佳值 [200]。
另一方面，对于跨声速和超声速射流，建议取 $ C_{DES} = 0.1 $ [201]。式（7.102）的长度尺度 l 的定义保证：在边界层内部，$ d < C_{DES} \Delta $，因此 l = d，从而恢复原有的 RANS 模型；在边界层外部，则有 $ l = C_{DES} \Delta $，Spalart–Allmaras 模型用作 LES 的单方程 SGS 模型（参见式（7.91）和（7.24））。对控制方程进行时间积分时，必须调整全局时间步长，以解析分离涡的时间尺度。这通常意味着所用时间步长会远远超过边界层区域中显式格式的稳定性限制。因此，采用保证时间精度的隐式格式更加高效，例如 6.3 节介绍的双时间推进方法。DES 方法的进一步说明及模拟示例见上述文献或 [202–209]。