第 14 章：源项离散、松弛及其它细节（Discretization of the Source Term, Relaxation, and Other Details）

14.1 源项离散（Source Term Discretization）

源项（汇与源）出现在流动和输运现象诸多问题的控制方程中，例子包括湍流模型方程、化学反应、辐射换热、传质以及多相流方程等。源项不仅影响问题的物理特性，还会影响数值计算的稳定性。然而若处理得当，源项能提升算法的鲁棒性。一条通用建议是：源项中的负值（汇）隐式处理，正值（源）显式评估。

源项的处理思路可以通过考察守恒方程在质心为 \(C\)、体积为 \(V_C\) 的单元上的离散形式来阐明（显式写出源项，图 14.1），方程为

\[ a_C \phi_C + \sum_{F \in nb(C)} a_F \phi_F = Q_C V_C \tag{14.1} \]

其中 \(Q_C V_C\) 是源项在单元 \(C\) 上的积分。一般地源项是因变量 \(\phi\) 的函数，关系式记作

\[ Q_C = Q(\phi_C) \tag{14.2} \]

在这种形式下，源项可以基于已有的 \(\phi\) 值显式算出——在迭代过程中，这些 \(\phi\) 值来自上一轮迭代。当 \(Q_C\) 是常数或相对其它项较小时，这种显式做法是可以接受的；但当 \(Q_C\) 随 \(\phi_C\) 的变化幅度较大时，迭代的收敛速度会变差。这种情况下，可以通过对 \(\phi_C\) 做 Taylor 展开来线性化 \(Q_C\)。用上标 \(\circ\) 标记上一轮迭代的值，当前迭代源项 \(Q_C\) 可写为

\[ Q(\phi_C) = Q^\circ_C + \left.\frac{\partial Q}{\partial \phi_C}\right|^\circ (\phi_C - \phi^\circ_C) = \left.\frac{\partial Q^\circ_C}{\partial \phi_C}\right|^\circ \phi_C + Q^\circ_C - \left.\frac{\partial Q^\circ_C}{\partial \phi_C}\right|^\circ \phi^\circ_C \tag{14.3} \]

右端第一项是隐式部分（依赖当前 \(\phi_C\)），第二项是显式部分（基于上一轮迭代值计算）。

在有限容积的语境下，方程 (14.1) 右端的源项积分可写为

\[ Q_C V_C = \int_{V_C} Q\,dV = \int_{V_C} \left.\frac{\partial Q^\circ_C}{\partial \phi_C}\right|^\circ \phi\,dV + \int_{V_C}\left(Q^\circ_C - \left.\frac{\partial Q^\circ_C}{\partial \phi_C}\right|^\circ \phi^\circ_C\right)dV = Flux_C \phi_C V_C + Flux_V V_C \tag{14.4} \]

代回 (14.1) 得到代数方程

\[ \left[a_C - Flux_C\right] \phi_C + \sum_{F\in nb(C)} a_F \phi_F = Flux_V \tag{14.5} \]

按此写法，源项的隐式部分 \(Flux_C\)（定义见 (14.3)）必须为负值以保证对角占优，否则可能不满足 Scarborough 准则而发散。此外，对于正定的 \(\phi\)，显式部分 \(Flux_V\)（见 (14.3)）必须为正以保证 \(\phi\) 的预测值非负。

例 1 在涉及辐射换热的问题中，能量方程的源项形式为

\[ Q_T = A\left(T_1^4 - T^4\right) \]

其中 \(A\) 为常数，\(T\) 为网格点温度，\(T_1\) 为不变的背景温度。将其在一个质心为 \(C\)、体积为 \(V_C\) 的单元上积分，并比较不同线性化方案对收敛性的影响。

解：

\[ \int_{V_C} Q_T dV = Q_C^T V_C = Flux_V + Flux_C\,T_C \]

源项 \(Q_T\) 可有多种任意的线性化选择。

方案 1 \(Flux_C = 0\)，\(Flux_V = A(T_1^4 - T_C^4)V_C\)。这一做法可能导致求解发散：当 \(T_C > T_1\) 时 \(Flux_V\) 为负，迭代过程中可能出现非物理的负绝对温度。

方案 2 对上一轮迭代的温度 \(T_C^\circ\) 做 Taylor 展开：

\[ Q_C^T = Q^\circ_C{}^T + \left.\frac{dQ_C^T}{dT_C}\right|^\circ (T_C - T_C^\circ) = A(T_1^4 - T_C^{\circ 4}) - 4A T_C^{\circ 3}(T_C - T_C^\circ) \]

与 (14.3) 对比可得

\[ Flux_C = -4A T_C^{\circ 3} V_C,\qquad Flux_V = A\left(T_1^4 + 3 T_C^{\circ 4}\right) V_C \]

这是理想方案：\(Flux_V\) 为正、\(Flux_C\) 为负，引入的隐式程度最合适，收敛速度也最佳。

14.2 代数方程的欠松弛（Under-Relaxation of the Algebraic Equations）

如前几章所述，离散过程的最终产物是一组形如 (14.1) 的代数方程，其中 \(a_F\) 是邻点系数（如图 14.2 所示），代表邻点变量 \(\phi_F\) 对中心变量 \(\phi_C\) 的影响；\(b_C\) 是方程右端，通常包含源项及来自其它变量的贡献；\(a_C\) 是代数方程的主对角系数，包含各种影响（空间离散效应、瞬态效应等）。方程组 (14.1) 通常是对角占优的。

在用迭代法求解该代数方程组时，常常希望减缓因变量在迭代之间的变化幅度。这对非线性问题的收敛是有益的，同时也能在初始猜测场远离解时避免发散。非线性可来自网格的非正交性、源项的存在以及方程本身的非线性特性等。一种常用的"放慢"（"松弛"）迭代过程中变量变化的方法即松弛法。许多 CFD 代码中标准使用的松弛法是 Patankar [1] 的隐式欠松弛法（第 8 章已简述）。文献中还提出其它欠松弛方法，如 E-Factor 法 [2] 和假瞬态法 [3]。Van Doormaal 与 Raithby [2] 指出这些松弛方法之间存在一定联系：任一方法中的欠松弛都可与其它方法中的欠松弛建立对应——它们都同样地减缓邻点和源对被松弛单元值的影响。换言之，欠松弛对相关单元中的源项与空间系数的作用是等价的。下面介绍其中一些松弛方法。

14.2.1 欠松弛方法（Under-Relaxation Methods）

解的松弛可以在每次迭代得到新解后显式地执行，也可以通过在求解之前把松弛效应纳入方程来隐式地实现。下面分别概述这两种做法。

14.2.2 显式欠松弛（Explicit Under-Relaxation）

在显式欠松弛方法中，每次迭代结束后遍历计算域内所有单元，按下式修正各单元 \(C\) 的预测值 \(\phi_C^{new,predicted}\)：

\[ \phi_C^{new,used} = \phi_C^{old} + \alpha_\phi \left(\phi_C^{new,predicted} - \phi_C^{old}\right) \tag{14.6} \]

其中 \(\alpha_\phi\) 为松弛因子。无论显式还是隐式松弛，\(\alpha_\phi\) 的取值含义如下：

1. \(\alpha_\phi < 1\)：欠松弛，可能减慢收敛速度但增强稳定性，即降低解的发散或振荡概率。
2. \(\alpha_\phi = 1\)：无松弛，迭代中预测的值直接用于下一轮。
3. \(\alpha_\phi > 1\)：过松弛，有时可加速收敛，但通常会降低稳定性。

显式欠松弛用于 SIMPLE 算法中的压力松弛（下一章介绍）。此外，在流体物性依赖于解并需迭代更新的问题中，显式欠松弛对促进收敛常常是必要的，例子包括湍流中的湍流黏度、可压缩流的密度、HR 格式下计算得到的界面值等。它也可用于松弛守恒方程中的单独项，例如源项，甚至某些情形下变量梯度本身。

14.2.3 隐式欠松弛方法（Implicit Under-Relaxation Methods）

本类方法中已发展出多种方案。第 8 章介绍过的 Patankar 方法 [1] 是标准方法，本节为完整起见再做概述。其它的方法还包括 E-factor 法与假瞬态法。

14.2.3.1 Patankar 欠松弛（Patankar's Under-Relaxation）

如上所述，迭代解的欠松弛可由松弛因子 \(\alpha_\phi\) 通过 (14.6) 表达。为简化隐式松弛的记号，将 (14.6) 改写为

\[ \phi_C = \phi^\circ_C + \alpha_\phi \phi_C^{new, iteration} - \phi^\circ_C \tag{14.7} \]

其中 \(\phi^\circ_C\) 是上一轮迭代的 \(\phi_C\)。在 Patankar 的松弛方法中，将 (14.7) 中的 \(\phi_C^{new, iteration}\) 用 (14.1) 等价表达式替换，得到

\[ \phi_C = \phi^\circ_C + \alpha_\phi \left[\frac{\sum_{F\in nb(C)} a_F \phi_F + b_C}{a_C} - \phi^\circ_C\right] \tag{14.8} \]

整理后方程变为

\[ \frac{a_C}{\alpha_\phi} \phi_C + \sum_{F\in nb(C)} a_F \phi_F = b_C + \frac{1 - \alpha_\phi}{\alpha_\phi} a_C \phi^\circ_C \tag{14.9} \]

在 (14.9) 中，松弛因子 \(\alpha_\phi\) 修改对角系数和右端而不改变方程的数学形式。由于 \(\alpha_\phi < 1\)，欠松弛增大了代数系统的对角占优性，从而增强迭代线性求解器的稳定性。这是相比显式方法的一个显著优点。

不过要注意的是，隐式松弛施加的关系正比于对角系数。因此对角系数越大松弛越强，亦即较小控制体积上的松弛更显著。下一节将演示这一性质。

14.2.3.2 E-Factor 松弛（E-Factor Relaxation）

E-Factor 方法 [2] 是 Patankar 方法的重新表述。先将 (14.1) 改写为

\[ a_C \phi_C = b_C - \sum_{F\in nb(C)} a_F \phi_F \tag{14.10} \]

对右端做欠松弛，得到

\[ a_C \phi_C = \alpha_\phi \left[b_C - \sum_{F\in nb(C)} a_F \phi_F\right] + (1 - \alpha_\phi) a_C \phi^\circ_C \tag{14.11} \]

将松弛因子替换为 \(\frac{E_\phi}{1 + E_\phi}\)，则 (14.11) 变为

\[ \frac{E_\phi}{1+E_\phi} \left[b_C - \sum_{F\in nb(C)} a_F \phi_F\right] + \frac{1}{1+E_\phi} a_C \phi^\circ_C \tag{14.12} \]

\[ \Rightarrow a_C \left(1 + \frac{1}{E_\phi}\right) \phi_C + \sum_{F\in nb(C)} a_F \phi_F = b_C + \frac{1}{E_\phi} a_C \phi^\circ_C \tag{14.13} \]

按此写法，欠松弛效应可以直观地解释为某个人工瞬态时间尺度在每次求解器迭代中推进 \(\phi_C\)。时间步长 \(\Delta t\) 与特征时间步 \(\Delta t^*\) 满足

\[ \Delta t = E_\phi \Delta t^* \tag{14.14} \]

其中

\[ \Delta t^* = \frac{\rho_C V_C}{a_C} \tag{14.15} \]

式 (14.15) 中 \(\rho_C\) 是单元 \(C\) 中流体的密度。特征时间步 \(\Delta t^*\) 对应于把 \(\phi_C\) 的变化扩散/对流过该单元所需的时间。因此 E-factor 相当于单元的 CFL 数。

由 (14.15) 可见，E-Factor 松弛的时间步推进依赖于单元体积：小单元中解的推进慢于大单元。这对稳态求解的收敛速度可能不利——实际中常在边界附近使用高拉伸比的细长网格，使关键区域的推进时间步远小于其它区域。Patankar 松弛法也有同样的特点。

\(E_\phi\) 与 \(\alpha_\phi\) 的关系为

\[ E_\phi = \frac{1}{1 - \alpha_\phi} \tag{14.16} \]

一般 \(E_\phi\) 在 4–10 之间取值，对应 \(\alpha_\phi\) 在 0.75–0.9 之间。

例 2 下图给出贴壁边界附近的示意网格。单元 \(A\) 与 \(D\) 代表贴壁细长网格，体积比约为 \(V_C^A/V_C^D \approx 0.1\)。此类网格下对角系数通常由扩散项主导，约为 \(a_C^A/a_C^D \approx 2\)，因为 \(A\) 含一个边界面。计算欠松弛因子 0.8 时两单元的相对伪瞬态时间步（图 14.3）。

解：先由所给松弛因子求等效 E-factor：

\[ E = \frac{1}{1 - \alpha} = \frac{1}{1 - 0.8} = 5 \]

各单元的伪时间步为

\[ \Delta t = E \Delta t^* = 5 \frac{\rho_C V_C}{a_C} \]

因此 \(A\) 与 \(D\) 的相对伪时间步为

\[ \frac{\Delta t_A}{\Delta t_D} = \frac{V_C^A}{a_C^A} / \frac{V_C^D}{a_C^D} = \frac{V_C^A}{V_C^D} \cdot \frac{a_C^D}{a_C^A} \approx 0.1 \times \frac{1}{2} = 0.05 \]

即 \(D\) 单元的伪时间步约为 \(A\) 单元的 20 倍。

14.2.3.3 假瞬态松弛（False Transient Relaxation）

假瞬态松弛法 [3] 是对欧拉一阶隐式瞬态法的修改，其中使用上一轮迭代值代替旧时间步值。与欧拉法相同，通过加入伪瞬态项 \(a^o_C \phi_C\) 到对角系数，并将伪旧时间步项 \(a^o_C \phi^\circ_C\) 加到右端，使代数方程的对角占优性增大。修改后的方程为

\[ \left(a_C + a^o_C\right) \phi_C + \sum_{F\in nb(C)} a_F \phi_F = b_C + a^o_C \phi^\circ_C \tag{14.17} \]

其中

\[ a^o_C = \frac{\rho_C V_C}{\Delta t} \tag{14.18} \]

\(a^o_C\) 与欧拉一阶隐式离散瞬态项所得的瞬态系数等价，\(\rho_C\) 是密度、\(V_C\) 是单元体积、\(\Delta t\) 是用户定义的假时间步长。当 \(\Delta t\) 取大值时，所加项可忽略，欠松弛效应也可忽略，方程的解与未松弛时相同。当 \(\Delta t\) 取极小值时，\(a^o_C\) 变得很大并主导其它项，解被强烈欠松弛，\(\phi_C\) 的变化极小（\(\phi_C \approx \phi^\circ_C\)）。

除了能在整个计算域中一致地推进解外，假瞬态法还能保证即使在对角系数为零的极端情况下也向对角系数提供非零贡献。

并不存在通用的最优欠松弛因子选择规则——某种情形下适用的值在另一种情形下可能不奏效。此外，不同方程可使用不同的欠松弛因子；计算域内同一时刻也不必使用同一欠松弛值；甚至同一变量的欠松弛值在迭代之间也可以变化。对于第 15、16 章将介绍的 SIMPLE 算法，Raithby 与 Schneider [4] 推导了速度与压力场欠松弛因子之间的最优关系，将在下一章给出。

14.3 方程的残差形式（Residual Form of the Equation）

迄今为止，离散代数方程一直写成"直接"或"标准"形式

\[ a_C \phi_C + \sum_{F\in nb(C)} a_F \phi_F = b_C \tag{14.19} \]

这也是 OpenFOAM® 中使用的形式。方程 (14.19) 也可以改写为"修正"或"残差"形式——重新整理各项以便求解使方程成立所需的修正量。设 \(\phi^\circ_C\) 与 \(\phi'_C\) 分别为 \(\phi_C\) 的上一轮迭代值与使之满足 (14.19) 所需的修正量，则

\[ \phi_C = \phi^\circ_C + \phi'_C \tag{14.20} \]

将 (14.20) 代入 (14.19) 整理得

\[ a_C (\phi^\circ_C + \phi'_C) + \sum_{F\in nb(C)} a_F (\phi^\circ_F + \phi'_F) = b_C \tag{14.21} \]

或

\[ a_C \phi'_C + \sum_{F\in nb(C)} a_F \phi'_F = b_C - \left[a_C \phi^\circ_C + \sum_{F\in nb(C)} a_F \phi^\circ_F\right] \tag{14.22} \]

式 (14.22) 的右端代表了场 \(\phi^\circ_C\) 上方程的残差误差。记单元 \(C\) 上该残差为 \(Res^\circ_C\)，(14.22) 变为

\[ a_C \phi'_C + \sum_{F\in nb(C)} a_F \phi'_F = Res^\circ_C \tag{14.23} \]

注意：对精确场而言 \(Res^\circ_C\) 为零。

虽然 (14.23) 与 (14.19) 在数学上等价，但前者在数值上有一个优点：当 \(\phi\) 较大而其变化较小时，残差形式下的数值误差略小于标准形式。

14.3.1 Patankar 欠松弛的残差形式（Residual Form of Patankar's Under-Relaxation）

Patankar 隐式松弛方程的残差形式可由将 (14.19) 改写为如下形式得到

\[ a_C (\phi^\circ_C + \phi'_C) = \alpha_\phi \left[b_C - \sum_{F\in nb(C)} a_F (\phi^\circ_F + \phi'_F)\right] + (1 - \alpha_\phi) a_C \phi^\circ_C \tag{14.24} \]

化简为

\[ a_C \phi'_C + \alpha_\phi \sum_{F\in nb(C)} a_F \phi'_F = \alpha_\phi \left[b_C - \left(a_C \phi^\circ_C + \sum_{F\in nb(C)} a_F \phi^\circ_F\right)\right] \tag{14.25} \]

注意到上式右端代表原方程的残差，故 (14.25) 可写为

\[ \frac{a_C}{\alpha_\phi} \phi'_C + \sum_{F\in nb(C)} a_F \phi'_F = Res^\circ_C \tag{14.26} \]

即以残差形式对代数方程做欠松弛，只需要修改对角系数。

14.4 残差与求解收敛（Residuals and Solution Convergence）

在任何迭代求解过程中，重要的是能够判断何时可以认为解已经足够好、何时可以估计误差低于某一容差、甚至守恒方程已被满足的精度如何。具备回答上述任何问题的工具是任何 CFD 代码的重要组成部分。这也可表述为：如何在不知道最终解的前提下评估解场的收敛程度。为此，多年来已提出多种指标——从简单地监测某点处变量在多次迭代间的变化，到监测某积分量（如阻力系数、总质量流、壁面剪切应力等）的演化，再到更常见的监测某种方程残差。其挑战在于：方法必须能适用于宽范围的流动参数以及各种几何与边界条件。

14.4.1 残差（Residuals）

在求解 (14.1) 所代表的离散方程组时，平衡方程的误差通过定义单元残差来量化

\[ Res^\circ_C = b_C - \left[a_C \phi_C + \sum_{F\in nb(C)} a_F \phi_F\right] \tag{14.27} \]

显然当解被求出、方程被满足时，\(Res^\circ_C\) 为零。利用 \(Res^\circ_C\) 可以构造全计算域上的多种残差指标，下面逐一介绍。

14.4.2 绝对残差（Absolute Residual）

按 (14.27) 的定义，残差可正可负。由于符号并不重要，通常使用 \(Res^\circ_C\) 的绝对值 \(R^\circ_C\) 来判断解是否收敛。若 \(R^\circ_C\) 随迭代下降则解在收敛，反之则在发散。点 \(C\) 处的 \(R^\circ_C\) 定义为

\[ R^\circ_C = \left|b_C - \left[a_C \phi_C + \sum_{F\in nb(C)} a_F \phi_F\right]\right| \tag{14.28} \]

14.4.3 最大残差（Maximum Residual）

当全计算域上绝对残差的最大值

\[ R^{\circ, max}_C = \max_{\text{all cells}} \left|b_C - \left[a_C \phi_C + \sum_{F\in nb(C)} a_F \phi_F\right]\right| = \max_{\text{all cells}} R^\circ_C \tag{14.29} \]

下降到某一小量 \(\epsilon\) 以下时，即认为解已收敛：

\[ R^{\circ, max}_C \leq \epsilon \;\Rightarrow\; \text{解已收敛} \tag{14.30} \]

14.4.4 均方根残差（Root-Mean Square Residual）

另一种用作收敛指标的参数是各单元绝对残差平方和的均方根 \(R^{\circ, rms}_C\)，数学上为

\[ R^{\circ, rms}_C = \sqrt{\frac{\sum_{C\in\text{all cells}} \left|b_C - \left[a_C \phi_C + \sum_{F\in nb(C)} a_F \phi_F\right]\right|^2}{\text{number of elements}}} = \sqrt{\frac{\sum_{C\in\text{all cells}} \left(R^\circ_C\right)^2}{\text{number of elements}}} \tag{14.31} \]

收敛判据写作

\[ R^{\circ, rms}_C \leq \epsilon \;\Rightarrow\; \text{解已收敛} \tag{14.32} \]

14.4.5 残差的归一化（Normalization of the Residual）

绝对残差的水平与变量 \(\phi\) 强相关。不同变量下的 \(R^\circ_C\) 量级不同，使得难以直接判断解是否收敛。一种更合理的做法是将各残差按其对应的最大通量进行缩放。注意到 \(a_C\) 代表单元上各通量之和，把残差对局部 \(\phi\) 值做相对化（除以全计算域上 \(a_C \phi_C\) 的最大值）即可得到相对误差

\[ R^{\circ, scaled}_C = \frac{\left|a_C \phi_C + \sum_{F\in nb(C)} a_F \phi_F - b_C\right|}{\max_{\text{all cells}} |a_C \phi_C|} \tag{14.33} \]

当全计算域上标度后的绝对残差最大值下降至小量 \(\epsilon\) 以下时，即认为解已收敛：

\[ \max_{\text{all cells}} R^{\circ, scaled}_C \leq \epsilon \;\Rightarrow\; \text{解已收敛} \tag{14.34} \]

工程中通常要求 \(\epsilon\) 对标度残差而言在 \(10^{-3}\) 至 \(10^{-5}\) 量级或更小才视为收敛。

除了使用绝对或标度残差之外，监测某些积分量（如前所述）也是有益的——必须在确认解真正收敛之后再下结论，因为未收敛的解可能产生误导。

14.5 计算实现指引（Computational Pointers）

本节讨论 uFVM 与 OpenFOAM® 中源项线性化与松弛技术的实现细节。

14.5.1 uFVM

uFVM 中源项的线性化与组装在函数 cfdAssembleSourceTerm（Listing 14.1）中实现。该函数依赖用户给定的线性化方式——用户须提供源项的常数部分 \(S_b\) 与线性化部分 \(S_c\)。这两项随后按 (14.4) 累加到 FLUXV 与 FLUXC 中。

值得指出的是，uFVM 中方程以残差形式求解，因此源项的总量被放进 FLUXTE，而不是只放其常数部分。

theEquationField = cfdGetMeshField(theEquationName);
phi = theEquationField.phi(iElements);
%
Sb = cfdComputeFormulaAtLocale(theTerm.Sb,'Interior Elements')';
Sc = cfdComputeFormulaAtLocale(theTerm.Sc,'Interior Elements')';
%
volume = [theMesh.elements.volume];
%
% Assemble Source Term
%
pos = zeros(1,size(phi));
pos(Sc<0) = 1;
theFluxes.FLUXCE = -pos .* Sc .* volume;
theFluxes.FLUXTE = -(Sb +Sc .*phi) .* volume;

Listing 14.1：源项的线性化与实现

uFVM 中实现了 Patankar 的隐式欠松弛方法。由于方程以残差形式求解，对方程做欠松弛（Listing 14.2）只需修改其对角系数（参 (14.26)）。

function cfdApplyURF(theEquationName)
%===================================================
% written by the CFD Group @ AUB, Fall 2006
%===================================================
theEquation = cfdGetModel(theEquationName);
urf = theEquation.urf;
theCoefficients = cfdGetCoefficients;
theCoefficients.ac = theCoefficients.ac/urf;
cfdSetCoefficients(theCoefficients);

Listing 14.2：Patankar 隐式欠松弛方法的实现

14.5.2 OpenFOAM®

14.1 节讨论了通用变量 \(\phi\) 的输运方程中源项的处理。所述的线性化（或隐式处理）可视为对系数矩阵施加一个人工时间步，从而影响解推进的特征时间。此外它增大了对角占优性，当负源项发生较大变化时系统能自适应地调整时间步以捕捉被模拟现象的特征，从而增强代数方程组的求解鲁棒性。这与不做线性化（显式处理）形成对比：后者需对整个方程组施以更重的欠松弛，且松弛因子通常并不最优。

在源项的离散上 OpenFOAM® [5] 使用隐式 fvm:: 与显式 fvc:: 算子。fvc:: 算子的实现位于目录 $FOAM_SRC/finiteVolume/finiteVolume/fvc/ 的 fvcSup.H 与 fvcSup.C 文件中。但通常做法是直接以不调用 fvc:: 算子的方式定义显式源项。例如对一个与主变量无直接关系、不能线性化的标量输运方程

\[ \nabla\cdot(\rho \mathbf{U} \phi) - \nabla\cdot(k \nabla \phi) = a U^2 \tag{14.35} \]

在 OpenFOAM® 中可按 Listing 14.3 实现：

fvMatrix<scalar> phiEqn
(
    fvm::div(mDot,phi) - laplacian(k,phi) == a*magSqr(U)
);

Listing 14.3：不调用 fvc 算子直接定义显式源项

此处源项无需任何专门的包装函数或算子。

fvm:: 函数的实现位于目录 $FOAM_SRC/finiteVolume/finiteVolume/fvm/ 的 fvmSup.H 与 fvmSup.C 文件中。线性化源项的离散在 fvm::Sp 函数中设置。按照 (14.5)，源项的隐式部分被加到系数矩阵的主对角上。Sp 函数的定义见 Listing 14.4。

template<class Type>
Foam::tmp<Foam::fvMatrix<Type> >
Foam::fvm::Sp
(
    const DimensionedField<scalar, volMesh>& sp,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

       tmp<fvMatrix<Type> > tfvm
       (
           new fvMatrix<Type>
           (
               vf,
               dimVol*sp.dimensions()*vf.dimensions()
           )
       );
       fvMatrix<Type>& fvm = tfvm();

       fvm.diag() += mesh.V()*sp.field();

       return tfvm;
}

Listing 14.4：定义与实现 Sp 函数所用的脚本

值得指出的是，Sp 函数对源项的处理与线性化形式的符号无关。这意味着在线性化项斜率为正的情况下，该操作可能导致解的发散，因为它会破坏代数方程组的对角占优性。因此务必保证只在隐式处理产生负的线性化斜率时才使用它。当线性化源项的斜率在计算域不同区域可正可负时，应把负贡献作隐式处理、正贡献作显式处理。OpenFOAM® 为此提供了一种特殊源项函数 fvm::SuSp，其隐式/显式处理自动完成。该函数脚本见 Listing 14.5。

template<class Type>
Foam::tmp<Foam::fvMatrix<Type> >
Foam::fvm::SuSp
(
    const DimensionedField<scalar, volMesh>& susp,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

      tmp<fvMatrix<Type> > tfvm
      (
          new fvMatrix<Type>
          (
              vf,
              dimVol*susp.dimensions()*vf.dimensions()
          )
      );
      fvMatrix<Type>& fvm = tfvm();

      fvm.diag() += mesh.V()*max(susp.field(), scalar(0));

      fvm.source() -= mesh.V()*min(susp.field(), scalar(0))
          *vf.internalField();

      return tfvm;
}

Listing 14.5：定义与实现 SuSp 函数所用的脚本

在此函数中，对角向量与源项向量都根据线性化源项斜率的局部符号分别填入。max/min 函数的调用实现了"按符号选择性离散"。例如，对计算域中任一单元，若线性化源项的斜率为负，则对源向量的贡献为零（即 min(SuSp.field(), scalar(0)) = 0），对角向量的贡献则相反。

OpenFOAM® 中的欠松弛方法既包括 Patankar 隐式法，也包括显式变量松弛。具体而言，隐式欠松弛只对 fvMatrix 对象（即实际的有限体积离散矩阵）施加；显式松弛只对 GeometricField 对象定义。

式 (14.6) 所述的显式松弛可在 GeometricField.C（位于 $FOAM_SRC/OpenFOAM/fields/GeometricFields/GeometricField 目录）中找到。GeometricField 类内执行此任务的专门函数（Listing 14.6）为：

template<class Type, template<class> class PatchField, class GeoMesh>
void    Foam::GeometricField<Type,   PatchField,   GeoMesh>::relax(const
scalar alpha)
{
     if (debug)
     {
         InfoIn
         (
             "GeometricField<Type, PatchField, GeoMesh>::relax"
             "(const scalar alpha)"
         ) << "Relaxing" << endl << this->info() << " by " << alpha <<
endl;
    }
    operator==(prevIter() + alpha*(*this - prevIter()));
}

Listing 14.6：GeometricField 类中显式欠松弛函数脚本

其中 == 运算符按 (14.6) 由当前值与上一轮迭代值定义 GeometricField 本身的新值。

一般地，要对一个变量做显式松弛，首先将其值存入 prevIter() 数组，然后执行计算得到新的预测值，最后施加松弛。例如以压力 p 作为 GeometricField 变量，应写（Listing 14.7）：

volScalarField p
(
    IOobject
    (
         "p",
         runTime.timeName(),
         mesh,
         IOobject::MUST_READ,
         IOobject::AUTO_WRITE
    ),
    mesh
);
    any preliminary operation
p.storePrevIter();
p=      Perform the operation for the new predicted pressure
p.relax();

Listing 14.7：压力场的显式欠松弛

在此情形下松弛因子的值直接从 fvSolution 字典读取。若开发者需要某个固定值（如 0.5），应将 Listing 14.5 末行替换为 p.relax(0.5);。

Patankar 松弛直接作用于系数矩阵，在 OpenFOAM® 中于 fvMatrix 类内实现。$FOAM_SRC/finiteVolume/fvMatrices/fvMatrix 目录下的 fvMatrix.C 文件含有 Listing 14.8 所列的隐式松弛定义。

template<class Type>
void Foam::fvMatrix<Type>::relax(const scalar alpha)
{
    if (alpha <= 0)
    {
        return;
    }

Listing 14.8：隐式松弛定义脚本

由于系数矩阵对角占优性要求的限制，函数定义较长，但相关核心代码为（Listing 14.9）：

      Field<Type>& S = source();
      scalarField& D = diag();

     // Store the current unrelaxed diagonal for use in updating the
source
    scalarField D0(D);

      // ... then relax
      D /= alpha;
      // Finally add the relaxation contribution to the source.
      S += (D - D0)*psi_.internalField();

Listing 14.9：隐式欠松弛脚本片段

在第一部分中，先检查线性化源项的符号，建立对角与源向量的引用。原始对角存在 D0 标量场中，然后除以松弛因子 alpha。再向矩阵的源向量加上一项附加贡献，其中 psi_ 是与 fvMatrix 类相关联的变量。源项与 (14.4) 右端形式上略有不同，但定义完全一致的贡献，如下式所示：

\[ \frac{(1 - \alpha)}{\alpha} a_C \phi^\circ_C = \frac{a_C}{\alpha} \phi^\circ_C - a_C \phi^\circ_C \tag{14.36} \]

14.6 本章小结（Closure）

本章讨论了通用守恒方程中源项的处理以及多种用于欠松弛代数方程组的方法，还介绍了若干用于检查收敛性的残差指标。下一章将专门讨论不可压缩流动问题的求解。

本章个人批注

本章是 Moukalled FVM 教材中少有的"小细节"专门章节，处理三件对求解稳定性与收敛性影响很大的事：源项的隐式/显式分裂、迭代过程中的欠松弛以及残差监测。

源项离散（14.1）的核心思想是"负隐正显"——隐式部分 \(Flux_C\) 必须为负、显式部分 \(Flux_V\) 必须为正，前者保证对角占优性（Scarborough 准则），后者保证正定变量的非负性。作者用辐射换热 \(Q_T = A(T_1^4 - T^4)\) 做了正反对照（Option 1 全显式导致 \(T_C > T_1\) 时 \(Flux_V\) 变负、可能算出负温度；Option 2 围绕 \(T_C^\circ\) 展开则保证两个条件都满足）。这在 CFD 教材里是少见的把"为什么不能简单全显式"讲清楚的地方——很多教材只给出"线性化"一笔带过。

欠松弛（14.2）的三种方法——Patankar 隐式、E-Factor、假瞬态——其实有共同本质：都是给代数方程加一项"对角强化"项，从而改善对角占优性。E-Factor 用 \(\frac{1}{1-\alpha}\) 重新表示松弛因子，可以被解释为某种伪时间步（\(\Delta t = E_\phi \Delta t^*\)）；假瞬态直接加 \(a^o_C = \rho V/\Delta t\) 到对角，与一阶隐式 Euler 完全同构。书中例 2 用贴壁细网格 \(V_A/V_D \approx 0.1\)、\(a_A/a_D \approx 2\) 算出 \(\Delta t_A/\Delta t_D \approx 0.05\)，是少见的把"边界层附近小网格拖慢整体收敛速度"定量化的例子——也是 Patankar 法和 E-factor 法共同的弱点。

残差形式（14.3、14.3.1）的关键观察是：把方程 (14.19) 改写成修正量形式 (14.23) 后右端就是 \(Res^\circ_C\)；而 Patankar 隐式松弛作用于残差形式时只需改对角系数（14.26），比直接形式更简洁。这是 uFVM 选择在残差形式下求解方程的根本原因（14.5.1.2 明确写出了"对角系数除以 urf"）。

残差监测（14.4）提供了四种收敛判据：绝对残差（14.28）、最大残差（14.30）、均方根残差（14.32）、归一化残差（14.34）。归一化残差通过除以 \(\max|a_C \phi_C|\) 把不同变量、不同量级的残差放到同一尺度下比较，是实践中最常用的指标；\(10^{-3}\) 到 \(10^{-5}\) 的判据是经验值。但作者特别强调："还应同时监测积分量（如阻力系数、总质量流）"——残差小不代表解对，因为残差只是逐点方程的不平衡度，积分量才是真正关心的物理量。这是 CFD 老手才会强调的注意事项。

计算实现方面（14.5），OpenFOAM® 的 fvm::Sp 与 fvm::SuSp 是同源但行为不同的源项算子——Sp 不论符号一律加到对角；SuSp 根据局部符号自动分流。14.5.2.1 末段提到"斜率为正时用 Sp 可能破坏对角占优性"，这与 14.1 节的"负隐正显"原则完全一致，可以看作是同一原则在 OpenFOAM® API 层的体现。

与第 8 章的关联：8.4 节已讲过 Patankar 隐式松弛的简要版本，本章是其详尽论述（变量/算子记号、为什么增大对角占优、对小控制体积的偏好）。与后续章节的关联：14.2 末预告 SIMPLE 算法下速度/压力欠松弛因子间的最优关系（Raithby-Schneider [4]）将在第 15、16 章给出；14.4 节的残差定义会被后续章节的求解器示例反复引用。

与上下章的衔接（一段话）

紧接第 13 章（Temporal Discretization: The Transient Term）讲完瞬态项的离散，第 14 章回过来处理代数方程层面"求解更稳"的三类细节——源项离散中的隐式/显式分裂、迭代的欠松弛（Patankar/E-Factor/假瞬态三种实现）以及残差监测（绝对/最大/均方根/归一化）。从结构上看，本章是离散过程结束之后、压力-速度耦合算法之前的过渡：13 章把瞬态项的算子形式给齐了，本章则讨论一旦代数方程组列好之后如何"喂给"迭代求解器更稳。第 15 章（The SIMPLE Algorithm）会立刻用上这里的 Patankar 隐式松弛（压力修正方程的处理）以及残差监测（SIMPLE 算法的收敛判据），并在 Raithby-Schneider 关系中引用 14.2 末段的最优欠松弛因子。第 16 章（Extensions of the SIMPLE Algorithm）继续推广 SIMPLE 类算法，第 17 章（turbulence）则进入湍流模型方程——而湍流源项（如 \(k\)-\(\epsilon\) 模型的 \(\epsilon\) 方程源项）正是 14.1 节线性化方法的应用对象。习题部分（14.7）从符号代数（Exercise 1、2、4、5）到工程实现（Exercise 6 用 uFVM/OpenFOAM 求解 Fithigh-Nagumo 模式生成方程）覆盖了本章所有核心方法，按惯例不纳入精读范围。