第 13 章：时间离散：瞬态项（Temporal Discretization: The Transient Term）
13.1 引言（Introduction）
本节先界定瞬态问题的时间离散思路。稳态情形无须对瞬态项做处理；瞬态模拟则要在空间离散之外增设时间离散。空间离散沿空间域定义，离散方程在控制容积上组装；瞬态离散沿时间坐标进行，对瞬态项的导数（差分法）或积分（有限体积法）进行评估（图 13.1）。
变量 \(\phi\) 的瞬态行为由如下方程控制：
其中 \(L(\phi)\) 是空间算子，包含所有非瞬态项（对流、扩散、源等），\(\partial(\rho\phi)/\partial t\) 为瞬态算子，二者都示于图 13.1。
将 (13.1) 在单元 \(C\) 上（图 13.2）积分，得到
经过对单元质心的空间离散，可写成
其中 \(V_C\) 是单元体积，\(L(\phi)_C^t\) 是空间离散算子在某个参考时间 \(t\) 上的值，其代数形式为
方程 (13.3) 中，当 \(t \to \infty\) 时回到稳态离散方程；通过时间推进达到稳态时（\(\phi_C^{t+\Delta t} = \phi_C^t\)）亦然，从而保证时间推进得到的稳态解与直接做稳态求解所得结果一致。
瞬态项的离散传统上沿用有限差分法（参考文献 [1–3]）：用 Taylor 展开将 \(\partial(\rho\phi)/\partial t\) 用离散节点值表达。本章还介绍一种更贴近有限体积思路的离散方法——把 \(\partial(\rho\phi)/\partial t\) 在一个时间单元 [4] 上积分，转化为面通量，与对流格式的构造方式相同，只是离散方向沿时间轴。
13.2.1 前向 Euler 格式（Forward Euler Scheme）
为了评估瞬态项，需要对函数在时间方向做 Taylor 展开。对前向情形，把函数 \(T\) 在 \(t\) 处的值通过 \(t+\Delta t\) 处的 \(T\) 与其导数表达为 Taylor 级数：
从 \(\Delta t^2\) 阶开始截断，可将一阶导数写成
这是一阶离散（除以 \(\Delta t\) 才得到梯度近似）。用 \((\rho\phi)\) 替换 \(T\)，把结果代入 (13.3)，离散方程变为
(13.7) 的瞬态模板（图 13.4）表明，\(t+\Delta t\) 处 \((\rho\phi)_C\) 的计算不需要解方程组；可以基于前一时间步的值显式算出 \(t+\Delta t\) 处的 \(\phi_C\)，因为所有空间项都在旧时刻 \(t\) 上评估。该格式属于显式瞬态格式 [5–12]。显式瞬态格式的主要特点是按时间推进求解，每一时间层无须解线性方程组，计算效率高、并行化容易。然而商业代码采用较少，原因与 \(\Delta t\) 的大小限制有关，将在下节讨论。
将空间算子的离散代数关系代入 (13.7)，得到完整代数方程：
其中
为简化记号，全章约定：前一时间步的变量用上标 \(^\circ\) 表示，前两时间步的变量用上标 \(^{\circ\circ}\)；当前时间步的变量不加角标（除瞬态项乘以 \(\phi_C\) 的系数 \(a_C^\delta\) 外）。这样 (13.8)–(13.9) 简化为
其中
方程 (13.10) 可改写为
明显表明当前时间步的 \(\phi\) 是通过显式关系计算的，无需求解方程组。
13.2.2 前向 Euler 格式的稳定性（Stability of the Forward Euler Scheme）
差分方程解的收敛性与稳定性最早由 Courant、Friedrichs 与 Lewy [13] 给出。他们证明差分方程的解要收敛于偏微分方程的解，数值格式必须使用初始数据中所有影响解的信息。这一要求后来被称为 CFL 条件。
CFL 条件可简单地理解为系数应满足的"异号规则"扩展到瞬态系数上：与 \(\phi_F\) 视为 \(\phi_C\) 的"空间邻居"类似，\(\phi_C^\circ\) 是 \(\phi_C\) 的"时间邻居"，异号规则同样适用。注意到对角系数是 \(a_C^\delta\)，而"时间邻居"的系数是 \((a_C + a_C^\circ)\)，异号要求为
13.2.2.1 瞬态对流情形（Transient-Advection Case）
对图 13.5 所示一维纯对流问题（流向右），在元素面用 upwind 格式插值所有变量时，离散方程中元素 \(C\) 的 \(a_C\) 与 \(a_C^\circ\) 系数为
CFL 条件要求
即
对对流主导问题，定义对流 CFL 数：
数值稳定性要求
13.2.2.2 瞬态扩散情形（Transient-Diffusion Case）
对纯扩散问题，CFL 数不同。以图 13.6 所示一维纯扩散问题为例，采用线性插值剖面时元素 \(C\) 的 \(a_C\) 与 \(a_C^\circ\) 系数为
CFL 条件要求
即
对网格均匀且扩散系数为常数的情形，(13.21) 简化为
对扩散主导问题，定义扩散 CFL 数：
稳定性要求
13.2.2.3 瞬态对流-扩散情形（Transient-Convection-Diffusion Case）
对多维非稳态对流-扩散问题（图 13.7），基于第 12 章推导，方程 (13.14) 中的系数为
把 (13.25) 代入 (13.13)，CFL 条件变为
即
方程 (13.27) 是显式瞬态格式稳定性的通用要求。一维纯扩散（均匀网格、密度 \(\rho\) 与扩散系数 \(\Gamma\) 为常数）与一维纯对流（upwind、流向从左到右）情形都可作为 (13.27) 的特例恢复：
这一稳定性约束非常严格，使瞬态问题求解必须使用极小的时间步：尽管每步计算量小于求解方程组，但 CFL 条件要求的时间步数极大。而且 (13.27) 还表明，提高空间精度（减小网格）会进一步缩小最大允许时间步。
接下来将看到，这种约束不适用于隐式格式，因为其瞬态项总是具有正确的符号。
13.2.3 后向 Euler 格式（Backward Euler Scheme）
推导后向 Euler 格式时，把函数 \(T\) 在 \(t-\Delta t\) 处的值用 \(t\) 处的 \(T\) 与其导数通过 Taylor 级数表达：
整理后得到一阶导数：
用 \((\rho\phi)\) 替换 \(T\)，代入 (13.3)，离散方程变为
引入空间算子的代数关系和前述记号后，完整的瞬态代数方程为
系数为
(13.33) 的模板（图 13.8）显示空间算子与新的时间系数处于同一时间层，在新时间层上求 \(\phi\) 场需要解方程组。该格式属于隐式格式 [5–12]。
由 (13.34) 可知 \(a_C\) 与 \(a_C^\circ\) 符号相反，保证 \(\phi_C\) 被当前时间步 \(t\) 上的空间邻居值与前一个时间步 \(t-\Delta t\) 的"时间邻居"值所界定。格式对任意 \(\Delta t\) 都稳定，可以用大步长快速推进。但它并非理想格式，因为精度低（一阶）；若要获得高精度只能减小 \(\Delta t\)，这与大步长的高效形成矛盾。
13.2.4 Crank-Nicolson 格式（Crank-Nicolson Scheme）
Crank-Nicolson（CN）格式 [2, 14] 把函数 \(T\) 在 \(t-\Delta t\) 与 \(t+\Delta t\) 处的值用 \(t\) 处的 \(T\) 与其导数表达为 Taylor 级数 [式 (13.35)]。
把两式相减，得到一阶导数：
注意此时导数的精度为 \(O(\Delta t^2)\)，因为二阶导数被完全消去。
把 (13.36) 的时间导数代入 (13.3)：
引入空间算子的代数关系和前述记号后，瞬态代数方程为
其中
(13.38) 的模板（图 13.9）显示该格式是显式的，因为 \((\rho\phi)^{t+\Delta t}\) 可以仅用旧时刻值计算；不过现在需要两个旧时刻，空间算子则在其中一个旧时刻上求值。
对 CN 格式的稳定性，可对原始方程稍作改写后分析。利用以下近似：
方程 (13.38) 变为 [式 (13.41)]。稳定性条件变为
对图 13.5 所示一维瞬态对流问题（用 upwind 离散对流项），(13.42) 给出
写成对流 CFL 数形式：
更大的 CFL 上限令人欣慰，但更重要的是精度提升——二阶导数误差被消去，无需极小的时间步即可获得高精度解。精度分析将在后文给出。
13.2.5 实现细节（Implementation Details）
CN 格式也可由前向 Euler 与后向 Euler 瞬态格式相加得到 [4]：
前向 Euler：
后向 Euler：
两者相加：
由此得到一种简单的两步实现：第一步用后向 Euler 隐式地求 \((\rho\phi)^t\)：
第二步显式地求 \(t+\Delta t\) 的 CN 值：
推导中假设瞬态时间步 \(\Delta t\) 被分为两个相等的局部时间步 \(\Delta t^{local}\)，且 \(\Delta t^{local}\) 等于设定时间步的一半。
需要注意的是，CN 格式虽然二阶精度，但仍是显式格式，受 CFL 类条件约束。
13.2.6 Adams-Moulton 格式（Adams-Moulton Scheme）
二阶 Adams-Moulton 格式 [15, 16] 需要把 \(T\) 在 \(t-\Delta t\) 与 \(t-2\Delta t\) 处的值用 \(t\) 处的 Taylor 展开 [式 (13.50)]。
把两个方程线性组合消去二阶导数，得到一阶导数：
代入 (13.3)：
展开空间项后，得到代数方程：
其中
显然 \(a_C^{\circ\circ}\) 系数为正，意味着 \(\phi_C^{\circ\circ}\) 增大时 \(\phi_C\) 会减小。这被 \(a_C^\delta\) 的大系数抵消（符号正确）。因此格式稳定但不有界，在某些情形会出现非物理解振荡。
Example 1
固体球（体积 \(1\,m^3\)）导热系数极高，球内温度梯度可忽略；初始温度 \(T_h\)，周围温度 \(T_1\)，密度 \(\rho\)、比热 \(c\)、表面积 \(A_s\)、对流换热系数 \(h_1\)。无辐射时能量方程为
定义无量纲温度 \(\phi = (T - T_1)/(T_h - T_1)\)，方程变为
设 \(h_1 A_S / (\rho c V) = 1\)，时间步 \(0.1\)，比较解析解与一阶显式、一阶隐式、二阶 CN 格式（两步实现）在 \(0.1\)、\(0.2\)、\(0.3\) 的结果。
求解：方程退化为 \(d\phi/dt = -\phi\)，\(\phi(0)=1\)。分离变量并应用初值得解析解 \(\phi(t) = e^{-t}\)。\(0.1\)、\(0.2\)、\(0.3\) 处的解析值分别为 \(0.9048\)、\(0.8187\)、\(0.7408\)。
数值求解取 \(V=1\)、\(L(\phi^n) = -\phi^n\)、\(L(\phi^{n+1}) = -\phi^{n+1}\)。误差定义为 \(|\phi_{numerical} - \phi_{exact}|\)。
一阶显式格式 \(\phi^{t+\Delta t} = (1 - \Delta t)\phi^t\)：
- \(\phi^{0.1} = 0.9,\ \phi^{0.2} = 0.81,\ \phi^{0.3} = 0.729\) - 误差：\(4.8 \times 10^{-3}\)、\(8.7 \times 10^{-3}\)、\(1.18 \times 10^{-2}\)
一阶隐式格式 \(\phi^{t+\Delta t} = \phi^t/(1+\Delta t)\)：
- \(\phi^{0.1} = 0.9091,\ \phi^{0.2} = 0.8264,\ \phi^{0.3} = 0.7513\) - 误差：\(4.3 \times 10^{-3}\)、\(7.7 \times 10^{-3}\)、\(1.05 \times 10^{-2}\)
二阶 CN 格式用 (13.48)–(13.49)，把总时间步 \(\Delta t\) 分为两个 \(\Delta t/2\)：
- \(\phi(0.05) = \phi(0)/(1+0.05) = 0.95238\)
- \(\phi_{CN}(0.1) = 2 \cdot 0.95238 - 1 = 0.90476\)
-
\(\phi(0.15) = 0.90476/(1+0.05) = 0.861678\)
\(\phi_{CN}(0.2) = 2 \cdot 0.861678 - 0.90476 = 0.81859\) -
\(\phi(0.25) = 0.81859/(1+0.05) = 0.779615\)
\(\phi_{CN}(0.3) = 2 \cdot 0.779615 - 0.81859 = 0.7406\) - 误差：\(7.55 \times 10^{-5}\)、\(1.366 \times 10^{-4}\)、\(1.854 \times 10^{-4}\)
13.3.1 一阶瞬态格式（First Order Transient Schemes）
下文将通过分别采用"上风格式"和"下风格式"时间插值剖面 [14, 17] 与 [4, 18]，构造一阶隐式与显式 Euler 瞬态格式。
13.3.2 一阶隐式 Euler 格式（First Order Implicit Euler Scheme）
一阶隐式 Euler 格式采用一阶"upwind"插值剖面 [14, 17]。如图 13.11 所示，时间单元面上 \(\rho\phi\) 的值取上游单元质心值：
代入 (13.57)：
即一阶隐式 Euler 格式。线性化为
13.3.2.1 数值扩散（Numerical Diffusion）
作为一阶格式，与对流格式的经验一样，预计会产生数值扩散。其大小可通过在时刻 \(t\) 做 Taylor 展开恢复原控制方程来估算。\((\rho\phi)^{t-\Delta t}\) 在 \(t\) 处展开：
整理得
其中第二项为"数值扩散项"。代入离散方程：
方程右端被加上一个数值扩散项，其大小与时间步成正比，类似于对流项 upwind 格式的数值扩散行为。因此格式虽然无条件稳定，但在大时间步下得到的解本质上是稳态解。
13.3.3 一阶显式 Euler 格式（First Order Explicit Euler Scheme）
一阶显式 Euler 格式采用一阶"downwind"插值剖面 [4, 18]。如图 13.12 所示，时间单元面上 \(\rho\phi\) 的值取下游单元质心值：
代入 (13.57)：
即一阶显式 Euler 格式。线性化为
注意此时新时刻在 \(t+\Delta t\) 上，(13.67) 的空间算子在时刻 \(t\) 评估。因此可以完全计算右端，无需解线性代数方程组即可得到 \(\rho\phi\) 在 \(t+\Delta t\) 的值。这就是显式格式，对应于"\(\rho\phi\) 在整个时间步内保持不变"的假设。
13.3.3.1 数值反扩散（Numerical Anti-Diffusion）
同样在 \(t\) 时刻做 Taylor 展开：
整理：
代入 (13.67)：
右侧二阶导项为负号，表现为"数值反扩散"，对剖面有压缩效应，与对流 downwind 格式类似。反扩散项也随时间步变化：当与 upwind 对流格式联用且 \(CFL_{conv}=1\) 时，对流格式的数值扩散与显式 Euler 的数值反扩散大小相等、符号相反，几乎对消得到精确解。但这并不实际，因为实际多维网格很难保证 \(CFL_{conv}=1\)。
与反扩散行为相关的是数值不稳定性，随 \(\Delta t\) 增大而加剧，对时间步形成很强的限制。这可以通过"负邻居系数规则"加以评估。
13.3.4 二阶瞬态 Euler 格式（Second Order Transient Euler Schemes）
与对流格式类似，二阶瞬态格式可以用线性插值剖面构造：对称剖面（中心差分）得到 Crank-Nicolson 格式 [2]；上风格式（二阶上风格式）[4, 19, 20] 得到 Adams-Moulton 格式 [15, 16]，即隐式二阶上风格式（SOUE）。
13.3.5 Crank-Nicholson 中心差分剖面（Crank-Nicholson Central Difference Profile）
用"上风格点"和"下风格点"之间的线性插值得到 \(\rho\phi\)，得到图 13.13 所示的 Crank-Nicholson 格式。对均匀时间步，可写成：
代入 (13.57)：
CN 格式的线性化系数为
图 13.9 模板表明该格式是显式的：\(t+\Delta t\) 处的值由 \(t\) 与 \(t-\Delta t\) 处的值显式算出。因此稳定性受 CFL 限制。
同样，与有限差分公式类似，可以用 (13.48)–(13.49) 改写为两步实现：先做一步一阶隐式 Euler，再做一步修正显式 Euler（外推形式）。
13.3.5.1 数值精度（Numerical Accuracy）
在 \(t\) 附近对 \((\rho\phi)^{t+\Delta t}\) 与 \((\rho\phi)^{t-\Delta t}\) 做 Taylor 展开 [式 (13.75) 与 (13.76)]，相减得到：
代入 (13.73)：
表明格式为二阶精度。三阶导项为色散项，会导致不稳定性。
13.3.6 二阶上风格式 SOUE（Second Order Upwind Euler Scheme）
采用图 13.14 所示的二阶"上风格"插值剖面，界面 \(\rho\phi\) 的近似为
代入 (13.57)：
即隐式二阶上风格式（SOUE）。该格式需要储存两个旧时刻的 \(\rho\phi\) 值，线性化系数为
13.3.6.1 数值精度（Numerical Accuracy）
通过 Taylor 级数评估可知该格式为二阶精度。在 \(t\) 附近展开 \((\rho\phi)^{t-\Delta t}\) 与 \((\rho\phi)^{t-2\Delta t}\) [式 (13.82)–(13.83)]，将 (13.82) 乘 4 与 (13.83) 相减，得到 SOUE：
将 (13.84) 与 (13.80) 合并，恢复的 \(\rho\phi\) 方程为
只有三阶数值色散项，没有数值扩散。
13.3.7 有限体积方法的初始条件（Initial Condition for the FV Approach）
有限体积公式的实现是直接的，例外在于初始时间步。如图 13.15 所示，第一个时间单元是时间边界单元，没有上游邻居。下单元面的 \(\rho\phi\) 值直接用于面，导致得到的梯度只有正确数值的一半——这是因为梯度计算用的是 \(\phi_C^{t_{initial}+\Delta t/2}\) 与 \(\phi_C^{t_{initial}}\) 之差（相隔半步），却除以整个 \(\Delta t\)，造成不可忽略的初始误差。
这可通过考察一阶隐式 Euler 格式在第一个时间单元上的离散方程看出。(13.57) 给出 \(\rho\phi\) 离散方程：
第一个时间单元的上风格插值产生的梯度等于 \((t_{initial}+\Delta t/2)\) 与 \(t_{initial}\) 处 \(\rho\phi\) 之差除以 \(\Delta t\)；而对正常单元（图 13.16），梯度等于 \((t_{initial}+3\Delta t/2)\) 与 \((t_{initial}+\Delta t/2)\) 处 \(\rho\phi\) 之差除以 \(\Delta t\)。两者差异显著，使用 (13.86) 的梯度作为起步会在初始产生大误差，影响后续步骤。如果采用类似图 13.16 的网格，则有限差分与有限体积方法的解基本相同。
采用此方案后，第一个时间单元（覆盖 \([t_{initial}+\Delta t/2;\, t_{initial}+3\Delta t/2]\) ）的上风格面值为
代入 (13.57)：
与任意内部单元的方程形式相同。
Example 2
使用 CN 格式（(13.73)）与 SOUE 格式重新做 Example 1，时间步 \(0.05\)，求 \(0.1\)、\(0.2\)、\(0.3\) 处的值。
解析解与 Example 1 相同。因为需要两个旧值，第一个时间步用一阶后向 Euler。数值求解用 (13.61)、(13.73)、(13.80) 化简：
按所述实现，第一个时间单元覆盖 \([0.025;\, 0.075]\)，第二个覆盖 \([0.075;\, 0.125]\)，依此类推。
CN 数值解：
- \(\phi_{EU}(0.05) = \phi(0)/(1+0.05) = 0.95238\)
- \(\phi_{CN}(0.15) = 0.95238 - 2 \cdot 0.05 \cdot 0.90476 = 0.861904\) - \(\phi_{CN}(0.2) = 0.90476 - 2 \cdot 0.05 \cdot 0.861904 = 0.81857\) - \(\phi_{CN}(0.25) = 0.861904 - 2 \cdot 0.05 \cdot 0.81857 = 0.780047\) - \(\phi_{CN}(0.3) = 0.81857 - 2 \cdot 0.05 \cdot 0.780047 = 0.74056\) - 误差：\(4 \times 10^{-5}\)、\(1.3 \times 10^{-4}\)、\(2.4 \times 10^{-4}\)
误差显示二阶精度。本例与 Example 1 的数值差异来自计算保留的小数位数。
SOUE 数值解：
- \(\phi_{EU}(0.05) = \phi(0)/(1+0.05) = 0.9524\)
\(\phi_{SOUE}(0.1) = (4 \cdot 0.9524 - 1)/3.1 = 0.90632\)
\(\phi_{SOUE}(0.15) = (4 \cdot 0.90632 - 0.9524)/3.1 = 0.86219\)
\(\phi_{SOUE}(0.2) = (4 \cdot 0.86219 - 0.90632)/3.1 = 0.82014\)
\(\phi_{SOUE}(0.25) = (4 \cdot 0.82014 - 0.86219)/3.1 = 0.780119\)
\(\phi_{SOUE}(0.3) = (4 \cdot 0.780119 - 0.82014)/3.1 = 0.74204\) - 误差：\(1.5 \times 10^{-3}\)、\(1.44 \times 10^{-3}\)、\(1.24 \times 10^{-3}\)
解为二阶精度，但精度低于 CN。
13.4.1 有限差分方法的非均匀时间步（FD）（Non-Uniform Time Steps with FD）
13.4.1.1 Crank-Nicolson 格式（FD）
如图 13.17 所示，非均匀时间步下的 CN 格式通过 Taylor 展开把 \(t+\Delta t\) 与 \(t-\Delta t^\circ\) 处 \(\rho\phi\) 用 \(t\) 处的值与导数表达 [式 (13.89) 与 (13.90)]。把 (13.89) 乘 \((\Delta t^\circ)^2\)、(13.90) 乘 \((\Delta t)^2\) 并相减，得到一阶导数：
把 (13.91) 的梯度代入 (13.3)，得到非均匀时间步 CN 格式的离散方程：
展开空间项后，代数方程为
时间相关系数为
均匀时间步下恢复 (13.39) 的系数。
13.4.2 Adams-Moulton 或 SOUE 格式（FD）（Adams-Moulton or SOUE Scheme FD）
如图 13.18 所示，非均匀时间步下的 Adams-Moulton 格式（即 SOUE 格式）通过 Taylor 展开把 \(t-\Delta t\) 与 \(t-\Delta t-\Delta t^\circ\) 处 \(\rho\phi\) 用 \(t\) 处的值与导数表达 [式 (13.95)–(13.96)]。把 (13.95) 乘 \((\Delta t+\Delta t^\circ)^2/\Delta t^2\) 与 (13.96) 相减，得到二阶的一阶导数（即 SOUE）：
把 (13.97) 的梯度代入 (13.3)，离散方程为
展开空间项后，代数方程为
时间相关系数为
均匀时间步下恢复 (13.54) 的系数。
13.4.3 有限体积方法的非均匀时间步（Non-Uniform Time Steps with FV）
沿用有限体积方法中的术语：时间单元的大小记为 \(\Delta t\)，相邻时间单元质心之间的距离记为 \(dt\)。均匀时间步下两者相等——有限差分与有限体积方法下两个连续计算场之间的时间都是 \(\Delta t = dt\)。变时间步下，对有限差分方法 \(\Delta t\) 不变，但对有限体积方法 \(dt = (\Delta t + \Delta t^\circ)/2\)，从而产生不同的公式。
与有限差分方法一样，变时间步下当前与旧时间步的值都影响格式的插值剖面，因此影响有限体积离散。这与在结构化非均匀网格上写对流格式剖面类似。下面以 CN 与 SOUE 格式为例说明该步骤。其他剖面可类推。
13.4.4 Crank-Nicolson 格式（FV）（Crank-Nicolson Scheme FV）
CN 格式的界面 \(\rho\phi\) 值由跨越界面的主节点的 \(\rho\phi\) 平均得到（图 13.19）：
代入 (13.57) 得到 \(\rho\phi\) 离散方程：
非均匀时间步 CN 格式的线性化系数为
与恒定时间步情形一样，方法为显式，需要储存两个旧时刻的值。\(\Delta t = \Delta t^\circ = \Delta t^{\circ\circ}\) 时退化为均匀时间步形式。
13.4.5 Adams-Moulton 或 SOUE 格式（FV）（Adams-Moulton or SOUE Scheme FV）
采用 (11.84) 给出的二阶"上风格"插值剖面，图 13.20 所示 \(t+\Delta t/2\) 与 \(t-\Delta t/2\) 处的界面 \(\rho\phi\) 值为
用该剖面代入 (13.57)，单元 \(C\)（图 13.2）上的 (13.1) 离散形式为
非均匀时间步 SOUE 格式的线性化系数为
与恒定时间步情形一样，方法为隐式，每步需要求解方程组。\(\Delta t = \Delta t^\circ = \Delta t^{\circ\circ}\) 时退化为 (13.80)。
13.5.1 uFVM
uFVM 中瞬态项的离散遵循有限体积方法并在隐式框架下实现。Listing 13.1 展示了一阶后向 Euler 格式产生的瞬态通量组装：
% theDensityField = cfdGetMeshField(['Density' theFluidTag]);
% density = theDensityField.phi(iElements);
% density_old = theDensityField.phi_old(iElements);
%
% volumes = [theMesh.elements(iElements).volume]';
%
% theFluxes.FLUXCE(iElements) = volumes .* density / dt;
% theFluxes.FLUXCEOLD(iElements) = - volumes .* density_old / dt;
% theFluxes.FLUXTE(iElements) = theFluxes.FLUXCE .* phi;
% theFluxes.FLUXTEOLD(iElements) = theFluxes.FLUXCEOLD .* phi_old;
Listing 13.1 由一阶隐式 Euler 格式产生的瞬态通量组装。
13.5.2 OpenFOAM®
OpenFOAM® 中显式与隐式时间导数 [21] 由命名空间 fvm、fvc 与对应函数 fvc::ddt(rho, phi) 与 fvm::ddt(rho, phi) 定义。还提供一阶与二阶上风格 Euler 格式以及二阶 Crank-Nicholson 格式，后者按两步方式实现。
瞬态格式文件位于目录 $FOAM_SRC/finiteVolume/finiteVolume/ddtSchemes。基类 ddtScheme<Type> 由所有时间离散格式继承。
一阶 Euler 格式在类 EulerDdtScheme 中实现。该类声明在基类 ddtScheme<Type> 之上，如 Listing 13.2 所示：
template<class Type>
class EulerDdtScheme
:
public ddtScheme<Type>
对应的 fvc 与 fvm 命名空间实现在文件 EulerDdtScheme.C 中。Euler 格式的隐式求值通过 Listing 13.3 的函数定义：
template<class Type>
tmp<fvMatrix<Type> >
EulerDdtScheme<Type>::fvmDdt
(
const volScalarField& rho,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
如 Listing 13.4 所示，该函数第一步定义 fvMatrix，仅填对角系数向量：
{
tmp<fvMatrix<Type> > tfvm
(
new fvMatrix<Type>
(
vf,
vf.dimensions()*dimVol/dimTime
)
);
fvMatrix<Type>& fvm = tfvm();
分配好存储对角系数和源项所需空间后，定义并计算需要存储的值。如 Listing 13.5 所示，先定义时间步的倒数 rDeltaT，再计算 \(a^\circ = \rho_C V_C / \Delta t\) 并存入 fvm.diag()。源项贡献 \(a^{\circ\circ} = -\rho_C^\circ V_C / \Delta t\) 乘以 vf 的旧值并存入 fvm.source()，其中 vf 是应用时间格式时的通用变量：
scalar rDeltaT = 1.0/mesh().time().deltaTValue();
fvm.diag() = rDeltaT*rho*mesh().V();
fvm.source() = rDeltaT*rho.oldTime()*vf.oldTime().internalField()*mesh().V();
return tfvm;
Listing 13.6 表明 OpenFOAM® 也允许 Euler 格式下显式求非稳态项，返回一个 GeometricField 对象，其中存有 \((\rho\phi - \rho^\circ\phi^\circ)/\Delta t\) 的值（此处为单位体积）：
tmp<GeometricField<Type, fvPatchField, volMesh> >
EulerDdtScheme<Type>::fvcDdt
(
const GeometricField<Type, fvPatchField, volMesh>& vf
)
return tmp<GeometricField<Type, fvPatchField, volMesh> >
(
new GeometricField<Type, fvPatchField, volMesh>
(
ddtIOobject,
rDeltaT*(vf*rho - vf.oldTime()*rho.oldTime())
)
);
SOUE 格式在 OpenFOAM® 中由类 backwardDdtScheme 实现。类定义如 Listing 13.7，位于基类之上：
template<class Type>
class backwardDdtScheme
:
public fv::ddtScheme<Type>
{
// Private Member Functions
//- Return the current time-step
scalar deltaT_() const;
//- Return the previous time-step
scalar deltaT0_() const;
此时 OpenFOAM® 使用当前与前一时间步的信息。一般情形下两个时间步不同，需要用两个变量储存。隐式时间离散方式与一阶瞬态格式相似，由 Listing 13.8 至 Listing 13.10 的函数定义：
template<class Type>
tmp<fvMatrix<Type> >
backwardDdtScheme<Type>::fvmDdt
(
const volScalarField& rho,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
Listing 13.9 计算乘以因变量当前值、旧值和旧旧值的瞬态系数，并把对对角系数的贡献存入 fvm.diag()：
scalar rDeltaT = 1.0/deltaT_();
scalar deltaT = deltaT_();
scalar deltaT0 = deltaT0_(vf);
scalar coefft = 1 + deltaT/(deltaT + deltaT0);
scalar coefft00 = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
scalar coefft0 = coefft + coefft00;
fvm.diag() = (coefft*rDeltaT)*rho.internalField()*mesh().V();
在 Listing 13.10 所示函数最后部分，非稳态项的贡献被计算并存入 fvm.source()：
fvm.source() = rDeltaT*mesh().V()*
(
coefft0*rho.oldTime().internalField()
*vf.oldTime().internalField()
- coefft00*rho.oldTime().oldTime().internalField()
*vf.oldTime().oldTime().internalField()
);
对比 SOUE 格式的系数与 (13.100)，可以清楚地看到 OpenFOAM® 对非稳态项采用有限差分方法，时间导数由 Taylor 级数展开近似。
13.6 收尾（Closure）
本章介绍了非稳态守恒方程瞬态项的离散。讨论了两类通用方法：基于有限差分离散的方法；以及在时间单元上积分守恒方程的有限体积方法。给出了一阶完全隐式与完全显式瞬态格式，研究了 CN 与 SOUE 格式在均匀与变时间步下的高阶近似。下一章将讨论源项离散、代数方程组松弛以及其他相关细节。
本章个人批注
读本章最直观的感受是：与第 11 章对流、第 12 章 HR 格式紧密耦合。瞬态格式的"上风格式"剖面对应第 11 章 upwind，二者 Taylor 展开思路一脉相承。
特别值得注意的"反扩散"现象（第 13.3.3.1 节）：显式 Euler 与对流 upwind 联用、CFL=1 时数值扩散与反扩散几乎对消，得到几乎精确的解——这是教科书给出的"完美搭档"。但作者紧接着承认这一组合不实用（实际多维网格很难保证 CFL=1），又揭示了不稳定性的根源——反扩散项随 Δt 增大而加剧，对 Δt 形成极强限制。这一点对我很重要：它解释了为什么商业 CFD 代码普遍回避显式格式（结合第 13.2.2.3 的通用稳定性约束），也更深刻地说明了"无条件稳定"不等于"任意大 Δt 都准确"——后向 Euler 精度低、CN 受 CFL 约束、Adams-Moulton 稳定但不振荡。隐式格式的选择本质上是 CFL、稳定性和精度之间的权衡。
第二个让我印象深刻的点是第 13.3.7 节初始条件的处理。第一个时间单元是边界单元，上风格剖面在面 t_{initial}+Δt/2 处取 t_{initial}（而非 t_{initial}+Δt），导致梯度只有正确值的一半。这种"边界单元"在时间维度上的特殊性，与第 8 章扩散项中"边界处的内部单元"处理思路类似，但时间维度的对称性更强（不像空间边界有几何不对称）。解决方法是引入"虚拟单元" t_{initial}，让第一个时间单元覆盖 [t_{initial}+Δt/2; t_{initial}+3Δt/2]，这样方程形式与内部单元一致。FVM 与 FDM 在均匀网格下结果相同，但在变时间步下不等价——作者清楚地区分了 Δt（单元大小）和 dt（相邻单元质心距）。
第三个观察：第 13.2.5 给出 CN 的两步实现——后向 Euler + 修正显式 Euler——与第 13.4 节非均匀时间步 CN 的 FD 实现形成对比。后者直接修改插值系数，时间导数的精度仍为 O(Δt²) 但空间导数不在时间单元中心。这体现了 FVM 在变步长下的灵活性。
第 13.5 节 Computational Pointers 揭示了 uFVM 与 OpenFOAM® 实现的关键差异：OpenFOAM® 在 backwardDdtScheme 中实际上是用 Taylor 展开（FDM 思路）实现非均匀 SOUE（与方程 (13.100) 对比可直接看出），而 uFVM 严格遵循 FVM 思路。两者殊途同归，但 OpenFOAM® 的实现需要用 deltaT、deltaT0、deltaT00 三个时间步变量，而 uFVM 用 dt 即可——这反映了两种离散哲学的差别。
最后，Example 1 与 Example 2 的对比数据：Example 1 中 CN 误差 7.55·10⁻⁵、1.366·10⁻⁴、1.854·10⁻⁴（二阶精度）；Example 2 中 CN 误差 4·10⁻⁵、1.3·10⁻⁴、2.4·10⁻⁴，SOUE 误差 1.5·10⁻³、1.44·10⁻³、1.24·10⁻³（也是二阶精度但比 CN 差约 30 倍）。SOUE 的色散项 (Δt²/3) 比 CN 的色散项 (Δt²/6) 系数大一倍——这是 SOUE 精度差的核心原因。
与上下章的衔接（一段话）
本章承接第 12 章高分辨率（HR）对流格式，并指向第 14 章源项离散与代数方程组松弛。第 12 章讨论的是空间方向的有界对流格式（NVF/NVD/TVD 框架），本章则将相同的"上风格 vs 下风格 vs 中心差分"剖面思想沿时间方向类比：瞬态项的离散在结构化时间网格上自然适合有限差分（Taylor 展开），但作者特别把它"再实现"为有限体积法——通过把 ρϕ 在时间单元 [t - Δt/2; t + Δt/2] 上积分，得到与对流项离散完全对称的 Flux_T 表达式，从而把第 11 章建立的"通量线性化"工具（Flux_C、Flux_C^°、Flux_V）扩展到时间维度。这种对称性也体现在"显式 vs 隐式 vs 二阶"三档精度梯度上：与对流格式的 upwind / CD / SOU 完全类比。CN 格式由前向 + 后向 Euler 求和得到（第 13.2.5），是这一类比的极端例子。第 13.3.7 节对初始时间单元（"边界单元"）的处理与第 8 章扩散项中边界内部单元的处理思路相似——都是离散方法在边界处的特殊处理，体现了"内部模板 vs 边界模板"这一贯穿全书的主题。第 13.6 节 Closure 明确指出"下一章将讨论源项离散、代数方程组松弛及其他相关细节"，与第 14 章衔接。