第 8 章：空间离散——扩散项（Spatial Discretization: The Diffusion Term）

8.1 二维矩形域中的扩散（Two-Dimensional Diffusion in a Rectangular Domain）

本章详细描述扩散项（由空间 Laplacian 算子表达）的离散。扩散与对流是两种不同的物理现象，从数值角度看必须分别处理、采用不同的插值剖面与权衡。本节先把注意力集中在一个具有规则 Cartesian 网格的简单矩形域上,讨论稳态扩散方程在源项存在时的离散。

被离散的稳态扩散方程写作 $\nabla\cdot(\Gamma_\phi\nabla\phi) = Q_\phi$（式 8.1），其中 $\phi$ 是标量（如温度、化学组分质量分数、湍动能等），$Q_\phi$ 是域内每体积 $\phi$ 的生成率，$\Gamma_\phi$ 是扩散系数。更一般地可写成 $\nabla\cdot\mathbf{J}_{\phi,D} = Q_\phi$（式 8.2），其中 $\mathbf{J}_{\phi,D} = -\Gamma_\phi\nabla\phi$（式 8.3）。

沿用第 5 章的一阶段离散,式 8.1 写成

\[ \sum_{f\sim nb(C)}(-\Gamma_\phi\nabla\phi)_f\cdot\mathbf{S}_f = Q_{\phi,C}V_C \]

（式 8.4）。把 $e/w/n/s$ 四个面展开就得到式 8.5 的四面通量形式。对一个均匀 Cartesian 网格,四个面的外法向面积矢量为:$\mathbf{S}_e = +(\Delta y)_e\mathbf{i}$、$\mathbf{S}_w = -(\Delta y)_w\mathbf{i}$、$\mathbf{S}_n = +(\Delta x)_n\mathbf{j}$、$\mathbf{S}_s = -(\Delta x)_s\mathbf{j}$（式 8.6）。

以东面 $e$ 为例,扩散通量

\[ \mathbf{J}_{\phi,D}|_e\cdot\mathbf{S}_e = -\Gamma_{\phi,e}(\partial\phi/\partial x)|_e(\Delta y)_e \]

（式 8.7）。按照之前对单个积分点的离散守恒方程,可把面通量写成三段 $\text{Flux}_e = \text{FluxC}_e\phi_C + \text{FluxF}_e\phi_F + \text{FluxV}_e$ （式 8.8）——要得到这三个系数,需要一个描述两侧单元形心之间 $\phi$ 变化规律的剖面（profile）。

假设 $\phi$ 在两单元形心之间线性变化（图 8.2），则东面 $e$ 处沿 $\mathbf{i}$ 方向的梯度可以写成

\[ (\partial\phi/\partial x)|_e = (\phi_E - \phi_C)/(dx)_e \]

（式 8.9）。把它代回式 8.8,得到沿面 $e$ 的扩散通量离散形式

\[ \text{Flux}_T^e = \Gamma_{\phi,e}(\Delta y)_e/(dx)_e\,(\phi_C - \phi_E) = \text{FluxC}_e\phi_C + \text{FluxF}_e\phi_E + \text{FluxV}_e \]

（式 8.10）。

定义几何扩散系数 $g_{Diffe} = |\mathbf{S}_e|/|\mathbf{d}_{CE}| = \mathbf{S}_e/\mathbf{d}_{CE}$ （式 8.11）后，式 8.10 的三个系数变为 $\text{FluxC}_e = \Gamma_{\phi,e}g_{Diffe}$、$\text{FluxF}_e = -\Gamma_{\phi,e}g_{Diffe}$、$\text{FluxV}_e = 0$（式 8.12）。对西面 $w$ 同样有 $\text{FluxC}_w = \Gamma_{\phi,w}g_{Diffw}$、$\text{FluxF}_w = -\Gamma_{\phi,w}g_{Diffw}$、$\text{FluxV}_w = 0$（式 8.14）以及 $g_{Diffw} = |\mathbf{S}_w|/|\mathbf{d}_{CW}|$。北面 $n$ 与南面 $s$ 的对应关系完全类比（式 8.15–8.17）。

把四个面的离散通量代回式 8.5,扩散方程的代数形式为

\[ a_C\phi_C + \sum_{F=a_N,a_S,a_E,a_W} a_F\phi_F = b_C \]

（式 8.18），其中 $a_F = \text{FluxF}_f = -\Gamma_{\phi,f}g_{Difff}$、

\[ a_C = \sum_{f\sim nb(C)}\text{FluxC}_f = -(a_E+a_W+a_N+a_S) \]

、

\[ b_C = Q_{\phi,C}V_C - \sum_{f\sim nb(C)}\text{FluxV}_f \]

（式 8.19、8.21）。紧凑形式即

\[ a_C\phi_C + \sum_{F=NB(C)} a_F\phi_F = b_C \]

,下标 $F$ 标邻居单元（$E,W,N,S$），下标 $f$ 标邻居面（$e,w,n,s$）。

8.2 离散方程的若干评注（Comments on the Discretized Equation）

一个合格的离散方法应当让最终的代数方程反映出原守恒方程的特性。前一章已经介绍了离散技巧的性质,这里补充两条离散化系数必须满足的额外规则。

8.2.1 零和规则（The Zero Sum Rule）

回头看离散过程,所做的第一个主要近似就是假定 $\phi$ 在共享某面的两单元形心之间线性变化。读者可能会问:为什么用一阶剖面而不用更高阶?为了回答这个问题,考虑一个无源项的一维构型,此时离散方程化简为 $a_C\phi_C + a_E\phi_E + a_W\phi_W = 0$（式 8.22），其中 $a_E = -\Gamma_{\phi,e}g_{Diffe}$、$a_W = -\Gamma_{\phi,w}g_{Diffw}$、$a_C = -(a_E+a_W)$（式 8.23）。在没有源/汇的情况下,$\phi$ 沿域仅由扩散（按 Fourier 定律,沿 $\phi$ 减小方向）传递,因此 $\phi_e$ 或 $\phi_w$ 应当落在 $\phi_C$ 与 $\phi_E$（或 $\phi_C$ 与 $\phi_W$）之间——线性剖面天然保证这一点。二阶剖面（如抛物线型）可能让面值高于或低于两侧形心值,这是非物理的;更高阶剖面也有同样的问题。所以,若离散格式要保证物理结果,应使用线性剖面。需要补充的是,所选剖面在单元尺寸减小时其影响会变小,因为当尺寸趋近零时,所有近似都应回到相同的解析解。

进一步,无源项时多维热传导方程化简为 $\nabla\cdot(\Gamma_\phi\nabla\phi) = 0$（式 8.24），这意味着 $\phi$ 与 $\phi + \text{const}$ 都是守恒方程的解。一个一致的离散方法应当通过其代数方程反映出这一性质,即代数方程应同时满足

\[ a_C\phi_C + \sum a_F\phi_F = 0 \]

与

\[ a_C(\phi_C + \text{const}) + \sum a_F(\phi_F + \text{const}) = 0 \]

（式 8.25），由此得到

\[ a_C + \sum a_F = 0 \]

（式 8.26）或等价地

\[ \sum a_F/a_C = -1 \]

（式 8.27）。式 8.25/8.26 实际由式 8.19 直接满足,有/无源项时都成立。

由此 $\phi_C$ 可视为其邻居的加权和;在无源项时 $\phi_C$ 始终被邻居值 $\phi_F$ 包夹（bounded）。当存在源项 $S/C \neq 0$ 时,$\phi_C$ 不必再被包夹,允许 over/undershoot,这是完全物理的——其越界幅度由 $S/C$ 与邻居节点系数 $a_F$ 的相对大小决定。

8.2.2 异号规则（The Opposite Signs Rule）

上面的推导已经表明,系数 $a_C$ 与 $a_F$ 是异号的。这在物理上意味着:当 $\phi_F$ 增大/减小时,$\phi_C$ 应当相应地增大/减小,本质上是有界性(boundedness)的要求。保证该性质的一个充分条件就是主系数与邻居系数异号;若不满足,有界性可能无法被强制。

8.3 边界条件（Boundary Conditions）

众所周知,任何常/偏微分方程的解析解都需要由相应的边界条件确定若干常数,不同的边界条件会得到不同的解——尽管一般方程不变。数值解也遵循同一约束,必须正确且精确地实现边界条件,否则由数值近似引入的细微变化就会让所求问题的解出错。

对于扩散问题,本章按与各项的相关性,讨论四类边界条件:Dirichlet、Neumann、mixed 以及 symmetry。边界条件施加在边界单元上——边界单元有一个或多个面位于域边界上。$\phi$ 的离散值既存储在边界单元的形心处,也存储在边界面形心处。

设 $C$ 为图 8.3 所示边界单元的形心,该单元有一边界面的形心为 $b$、外指面积为 $\mathbf{S}_b$。按前几章的离散方法,作用于单元 $C$ 的离散过程给出

\[ \sum_{f\sim nb(C)} \mathbf{J}_{\phi,D}|_f\cdot\mathbf{S}_f = Q_{\phi,C}V_C \]

（式 8.28）。内部面的通量按之前离散,边界通量则按"对 $\phi_C$ 线性化"的原则离散,即

\[ \mathbf{J}_b^{\phi,D}\cdot\mathbf{S}_b = \text{Flux}_T^b = -\Gamma_{\phi,b}(\nabla\phi)_b\cdot\mathbf{S}_b = \text{FluxC}_b\phi_C + \text{FluxV}_b \]

（式 8.29）。

边界条件的指定有两种方式:一是直接指定未知的边界值 $\phi_b$,二是指定边界通量 $\mathbf{J}_b^{\phi,D}$。下面利用式 8.18 给出不同边界条件下扩散问题边界单元的离散方程。

8.3.1 Dirichlet 边界条件（Dirichlet Boundary Condition）

Dirichlet 条件指定边界上的 $\phi$ 值,即 $\phi_b = \phi_\text{specified}$（式 8.30, 图 8.3）。对这种情况,

\[ \text{Flux}_T^b = -\Gamma_{\phi,b}(\nabla\phi)_b\cdot\mathbf{S}_b = -\Gamma_{\phi,b}|\mathbf{S}_b|/|\mathbf{d}_{Cb}|\,(\phi_b - \phi_C) = \text{FluxC}_b\phi_C + \text{FluxV}_b \]

,给出 $\text{FluxC}_b = \Gamma_{\phi,b}g_{Diffb} = a_b$、 $\text{FluxV}_b = -\Gamma_{\phi,b}g_{Diffb}\phi_b = -a_b\phi_b$ （式 8.32），其中 $g_{Diffb} = \mathbf{S}_b/\mathbf{d}_{Cb}$（式 8.33）。

对图 8.3 所示的单元 $C$,$a_E$ 系数为零,离散方程简化为 $a_C\phi_C + a_W\phi_W + a_N\phi_N + a_S\phi_S = b_C$ （式 8.34），其中 $a_E = 0$、$a_W = \text{FluxF}_w = -\Gamma_{\phi,w}g_{Diffw}$、$a_N = \text{FluxF}_n = -\Gamma_{\phi,n}g_{Diffn}$、$a_S = \text{FluxF}_s = -\Gamma_{\phi,s}g_{Diffs}$、

\[ a_C = \text{FluxC}_b + \sum_{f\sim nb(C)}\text{FluxC}_f \]

（式 8.35）;

\[ b_C = Q_{\phi,C}V_C - [\text{FluxV}_b + \sum_{f\sim nb(C)}\text{FluxV}_f] \]

（式 8.36）。

关于该边界离散方程可以做出如下重要观察:

系数 $a_b$ 大于其他邻居系数,因为 $b$ 离 $C$ 更近,对 $\phi_C$ 影响更大。
系数 $a_C$ 仍然是所有邻居系数的代数和（包含 $a_b$）。这意味着对边界单元有 $\sum_{F=NB(C)}|a_F|/|a_C| < 1$,这是满足 Scarborough 准则的第二个必要条件,从而保证线性方程组在任一次迭代下都能通过迭代法收敛。
乘积 $a_b\phi_b$（即 $\text{FluxV}_b$）被移到了方程右端、成为 $b_C$ 的一部分,因为它不包含未知量。

8.3.2 Neumann 边界条件（Von Neumann Boundary Condition）

当边界上指定的是通量（或法向梯度,图 8.4）时,称 Neumann 条件。此时指定的通量为

\[ -(\Gamma_\phi\nabla\phi)_b\cdot\mathbf{i} = q_b \]

（式 8.37）,它实质上就是

\[ \mathbf{J}_b^{\phi,D}\cdot\mathbf{S}_b = -(\Gamma_\phi\nabla\phi)_b\cdot|\mathbf{S}_b|\mathbf{i} = q_b|\mathbf{S}_b| = \text{FluxC}_b\phi_C + \text{FluxV}_b \]

（式 8.38），其中 $\text{FluxC}_b = 0$、$\text{FluxV}_b = q_b\mathbf{S}_b = q_b(\Delta y)_C$（式 8.39）。这里假定通量分量在与其所作用的坐标系同向时为正。

把 $q_b$ 直接代入式 8.18,得到边界单元 $C$ 的离散方程 $a_C\phi_C + a_W\phi_W + a_N\phi_N + a_S\phi_S = b_C$ （式 8.40），其中 $a_E = 0$、$a_W = -\Gamma_{\phi,w}g_{Diffw}$、$a_N = -\Gamma_{\phi,n}g_{Diffn}$、$a_S = -\Gamma_{\phi,s}g_{Diffs}$、

\[ a_C = \sum_{f\sim nb(C)}\text{FluxC}_f = -(a_W+a_N+a_S) \]

、

\[ b_C = Q_{\phi,C}V_C - [\text{FluxV}_b + \sum_{f\sim nb(C)}\text{FluxV}_f] \]

（式 8.41）。

关于该离散方程可作以下重要说明:

Neumann 条件不会让 $a_C$ 成为支配性系数。
若 $q_b$ 与 $S/C$ 均为零,则 $\phi_C$ 仍由邻居值包夹。否则 $\phi_C$ 可以超过（或低于）邻居 $\phi$ 值,这是允许的——例如当 $\phi$ 表示温度时,$q_b$ 代表施加在边界上的热通量,若在边界处加热,临近边界的温度理应高于内部。
一旦求出 $\phi_C$,可用 $\phi_b = (\Gamma_{\phi,b}g_{Diffb}\phi_C - q_b)/(\Gamma_{\phi,b}g_{Diffb})$ （式 8.42）回算边界值 $\phi_b$。
最后,Neumann 条件可视为有限体积法的自然边界条件:当通量指定为零时,对该面无需任何离散处理;而 Dirichlet 条件（指定零值）则仍需进行离散。

8.3.3 混合边界条件（Mixed Boundary Condition）

混合边界条件（图 8.5）指边界信息通过对流换热系数 $h$ 和外部 $\phi$ 值 $\phi_\infty$ 给出,即

\[ \mathbf{J}_b^{\phi,D}\cdot\mathbf{S}_b = -(\Gamma_\phi\nabla\phi)_b\cdot\mathbf{S}_b = -h(\phi_\infty - \phi_b)(\Delta y)_C \]

（式 8.43），可改写为

\[ -(\Gamma_{\phi,b}\mathbf{S}_b/dx_b)(\phi_b - \phi_C) = -h(\phi_\infty - \phi_b)\mathbf{S}_b \]

（式 8.44）。由该式解出 $\phi_b$:

\[ \phi_b = \frac{h\phi_\infty + (\Gamma_{\phi,b}/dx_b)\phi_C}{h + \Gamma_{\phi,b}/dx_b} \]

（式 8.45）。代回式 8.43 后,通量方程变为

\[ \mathbf{J}_b^{\phi,D}\cdot\mathbf{S}_b = -\underbrace{\frac{h(\Gamma_{\phi,b}/dx_b)}{h + \Gamma_{\phi,b}/dx_b}\mathbf{S}_b}_{R_{eq}}(\phi_\infty - \phi_C) = \text{FluxC}_b\phi_C + \text{FluxV}_b \]

（式 8.46），其中 $\text{FluxC}_b = R_{eq}$、$\text{FluxV}_b = -R_{eq}\phi_\infty$（式 8.47）。

将式 8.47 代入边界单元 $C$ 的离散方程,得到修改后的形式 $a_C\phi_C + a_W\phi_W + a_N\phi_N + a_S\phi_S = b_C$ （式 8.48），其中 $a_E = 0$、$a_W = -\Gamma_{\phi,w}g_{Diffw}$、$a_N = -\Gamma_{\phi,n}g_{Diffn}$、$a_S = -\Gamma_{\phi,s}g_{Diffs}$、

\[ a_C = \text{FluxC}_b + \sum_{f\sim nb(C)}\text{FluxC}_f \]

、

\[ b_C = Q_{\phi,C}V_C - [\text{FluxV}_b + \sum_{f\sim nb(C)}\text{FluxV}_f] \]

（式 8.49）。

8.3.4 对称边界条件（Symmetry Boundary Condition）

沿对称边界,标量 $\phi$ 的法向通量为零。因此对称边界条件等价于通量取零的 Neumann 条件,即 $\text{FluxC}_b = \text{FluxV}_b = 0$。因此对称边界条件下的修改方程可以直接由式 8.40/8.41 令 $q_b = 0$ 得到,即 $a_C\phi_C + a_W\phi_W + a_N\phi_N + a_S\phi_S = b_C$ （式 8.50），其中 $a_E = 0$、$a_W = -\Gamma_{\phi,w}g_{Diffw}$、$a_N = -\Gamma_{\phi,n}g_{Diffn}$、$a_S = -\Gamma_{\phi,s}g_{Diffs}$、

\[ a_C = \sum_{f\sim nb(C)}\text{FluxC}_f = -(a_W+a_N+a_S) \]

、

\[ b_C = Q_{\phi,C}V_C - \sum_{f\sim nb(C)}\text{FluxV}_f \]

（式 8.51）。

8.4 界面扩散系数（The Interface Diffusivity）

在前面的离散方程中,$\Gamma_{\phi,e}$、$\Gamma_{\phi,w}$、$\Gamma_{\phi,n}$、$\Gamma_{\phi,s}$ 表示 $\Gamma_\phi$ 在单元 $e/w/n/s$ 面上的值。当 $\Gamma_\phi$ 随位置变化时,其值在单元形心 $E,W,\ldots$ 上已知,需要某种规则由这些点上的值来求取界面值。下述讨论在 $\Gamma_\phi$ 处处相等时不适用。

把能量方程作为具体例子会更清楚:此时扩散系数代表材料的导热率,$\phi$ 代表温度。非均匀导热率出现在非均质材料和/或导热率随温度变化的情形。对一般 $\phi$ 的控制方程,$\Gamma_\phi$ 按相同方式处理。在湍流中,$\Gamma_\phi$ 可能代表湍流黏度或湍流导热率,$\Gamma_\phi$ 可能有显著变化——因此一个正确的非均匀 $\Gamma_\phi$ 公式非常重要。

最简单的方法是假设 $\Gamma_\phi$ 在 $C$ 与任一邻居之间线性变化。对东面有 $\Gamma_{\phi,e} = (1-g_e)\Gamma_{\phi,C} + g_e\Gamma_{\phi,E}$ （式 8.52），其中插值因子 $g_e = d_{Ce}/(d_{Ce}+d_{eE})$（式 8.53）。对一个 Cartesian 网格,若界面在两格点正中,$g_e = 0.5$,$\Gamma_{\phi,e}$ 即为 $\Gamma_{\phi,C}$ 与 $\Gamma_{\phi,E}$ 的算术平均。其他面有类似定义（式 8.54）。每个单元每个面的这些系数只需计算一次;对 Cartesian 网格,用距离和用体积得到的插值因子相同。

但这种最基本的方法在某些情况下不正确——例如它不能正确处理复合材料中可能出现的导热率突变。幸运的是,有一个同等简单但更好的替代方案。开发这一替代方案的核心理念是:界面上导热率的局部值并非主要关注点,主要目标是获得界面处扩散通量 $\mathbf{J}_{\phi,D}$ 的良好表示[1]。

对图 8.6 所示的一维问题,假设单元 $C$ 由导热率 $\Gamma_{\phi,C}$ 的材料组成,单元 $E$ 由导热率 $\Gamma_{\phi,E}$ 的材料组成。对 $C$ 与 $E$ 之间的非均质平板,稳态一维（无源）分析给出（界面 $e$ 两侧的通量相同）:

\[ \mathbf{J}_{\phi,D}|_e\cdot\mathbf{S}_e = \frac{\phi_C - \phi_e}{(dx)_{Ce}/\Gamma_{\phi,C}} = \frac{\phi_e - \phi_E}{(dx)_{eE}/\Gamma_{\phi,E}} = \frac{\phi_C - \phi_E}{(dx)_{Ce}/\Gamma_{\phi,C} + (dx)_{eE}/\Gamma_{\phi,E}} \]

（式 8.55）。由此得到平板的等效导热率为:

\[ (dx)_{CE}/\Gamma_{\phi,e} = (dx)_{Ce}/\Gamma_{\phi,C} + (dx)_{eE}/\Gamma_{\phi,E} \quad\Rightarrow\quad 1/\Gamma_{\phi,e} = (1-g_e)/\Gamma_{\phi,E} + g_e/\Gamma_{\phi,C} \]

（式 8.56）。当界面位于 $C$ 与 $E$ 正中（$g_e = 0.5$）时,式 8.56 化为

\[ \Gamma_{\phi,e} = 2\Gamma_{\phi,C}\Gamma_{\phi,E}/(\Gamma_{\phi,C}+\Gamma_{\phi,E}) \]

（式 8.57）——即 $\Gamma_{\phi,C}$ 与 $\Gamma_{\phi,E}$ 的调和平均而非算术平均。

需要强调:对间断扩散系数,调和平均插值在一维扩散中是精确的;但用于多维时仍有一项重要优势——共轭界面的处理无需特殊处理。固体单元与流体单元可以被视为同一域的一部分,各自在形心处存放不同的扩散系数;通过按调和平均计算面的扩散系数,共轭界面处的扩散通量就是正确的。

Example 1（例 8.4.1）展示了一个完整的二维矩形双材料稳态导热算例,边界条件为:左边界 T = 320 K(温度指定,Dirichlet),下边界热通量 $q_b = 100\,\text{W/m}^2$(Neumann),右边界零通量,上边界对流换热 $h = 20\,\text{W/m}^2\!\cdot\!\text{K}$、$T_\infty = 300\,\text{K}$(mixed)。材料 1 区域 $k_1 = 10^{-3}\,\text{W/m}\!\cdot\!\text{K}$,材料 2 区域 $k_2 = 10^2\,\text{W/m}\!\cdot\!\text{K}$,网格非均匀（$3\times 3$ 单元）。该例详细展示了如何用调和平均处理两侧不同材料界面、如何为每个单元写出代数方程,并用 Gauss-Seidel 迭代法求解。算例以能量守恒校验收尾:顶/左/底三边界的总热流量之和

\[ Q_\text{Top}+Q_\text{Left}+Q_\text{Bottom} = 50.058 - 0.053398 - 50 = 0.004602 \approx 0 \]

,与零的微小偏差归因于计算中保留的有效位数有限。

8.5 非 Cartesian 正交网格（Non-Cartesian Orthogonal Grids）

接下来考虑正交但不沿 $x,y$ 轴取向的网格。如图 8.8 所示,这种网格可以由图 8.1 的 Cartesian 网格旋转一定角度得到。对这种网格,所得到的离散方程应当与 Cartesian 网格得到的完全相同——对相似的边界条件,也应得到相同的解。

仍考虑稳态导热方程 $\nabla\cdot\mathbf{J}_{\phi,D} = Q_\phi$（式 8.58）。其离散形式仍为

\[ \sum_{f\sim nb(C)} (-\Gamma_\phi\nabla\phi)_f\cdot\mathbf{S}_f = Q_{\phi,C}V_C \]

（式 8.59）。

对面 $e$ 离散,有

\[ \mathbf{J}_{\phi,D}|_e\cdot\mathbf{S}_e = -\Gamma_{\phi,e}(\nabla\phi\cdot\mathbf{n})_e\mathbf{S}_e = -\Gamma_{\phi,e}(\partial\phi/\partial n)|_e\mathbf{S}_e \]

（式 8.60），其中

\[ (\nabla\phi\cdot\mathbf{n})_e = (\partial\phi/\partial n)|_e \]

（式 8.61）代表 $\phi$ 在面 $e$ 处沿 $\mathbf{n}$ 方向的梯度。再次假设 $\phi$ 沿 $n$ 坐标轴线性变化,该梯度可以写成

\[ (\partial\phi/\partial n)|_e = (\phi_E - \phi_C)/d_{CE} \]

（式 8.62）。其他项的离散过程与 Cartesian 网格完全相同,最终得到的离散方程也完全一致。

8.6.1 非正交性（Non-orthogonality）

在前述构型中,通量方向都垂直于面。而一般的结构化曲线网格与非结构化网格都是非正交的[2–5],因此面积矢量 $\mathbf{S}_f$ 与连接共享该面的两单元形心的矢量 $\mathbf{d}_{CF}$ 并不共线（图 8.9）。此时法向梯度不能写成 $\phi_F$ 与 $\phi_C$ 的函数,因为它在垂直于 $\mathbf{d}_{CF}$ 的方向上还有分量。

正交网格上,法向梯度为

\[ (\nabla\phi\cdot\mathbf{n})_f = (\partial\phi/\partial n)_f = (\phi_F-\phi_C)/|\mathbf{r}_F-\mathbf{r}_C| = (\phi_F-\phi_C)/d_{CF} \]

（式 8.63）,因为 $\mathbf{d}_{CF}$ 与 $\mathbf{n}$（面法向单位矢量）对齐。在非正交网格上[6,7],要使梯度只涉及 $\phi_F$ 与 $\phi_C$,方向必须沿连接 $C$ 与 $F$ 的直线。若 $\mathbf{e}$ 表示沿 $C,F$ 连线的单位矢量,则

\[ \mathbf{e} = (\mathbf{r}_F-\mathbf{r}_C)/|\mathbf{r}_F-\mathbf{r}_C| = \mathbf{d}_{CF}/d_{CF} \]

（式 8.64）,$\mathbf{e}$ 方向的梯度为

\[ (\nabla\phi\cdot\mathbf{e})_f = (\partial\phi/\partial e)_f = (\phi_F-\phi_C)/|\mathbf{r}_F-\mathbf{r}_C| = (\phi_F-\phi_C)/d_{CF} \]

（式 8.65）。

为了在非正交网格上实现通量的线性化,需要把面积矢量 $\mathbf{S}_f$ 写成两个矢量之和 $\mathbf{S}_f = \mathbf{E}_f + \mathbf{T}_f$（式 8.66），其中 $\mathbf{E}_f$ 沿 $\mathbf{d}_{CF}$ 方向,使得一部分扩散通量可以写成 $\phi_F$ 与 $\phi_C$ 的函数:

\[ (\nabla\phi)_f\cdot\mathbf{S}_f = \underbrace{(\nabla\phi)_f\cdot\mathbf{E}_f}_{\text{正交贡献}} + \underbrace{(\nabla\phi)_f\cdot\mathbf{T}_f}_{\text{非正交贡献}} = \mathbf{E}_f(\partial\phi/\partial e)_f + (\nabla\phi)_f\cdot\mathbf{T}_f = \mathbf{E}_f(\phi_F-\phi_C)/d_{CF} + (\nabla\phi)_f\cdot\mathbf{T}_f \]

（式 8.67）。式 8.67 右端第一项与正交网格上的贡献类似,涉及 $\phi_F$ 与 $\phi_C$;第二项称为交叉扩散(cross-diffusion)或非正交扩散[8],由网格的非正交性引起。下面讨论 $\mathbf{S}_f$ 分解的几种可选方案。

8.6.2 最小修正方法（Minimum Correction Approach）

如图 8.10 所示,该方法对 $\mathbf{S}_f$ 的分解使式 8.67 的非正交修正尽可能小——通过让 $\mathbf{E}_f$ 与 $\mathbf{T}_f$ 互相正交实现。非正交性越大,扩散通量中 $\phi_F$ 与 $\phi_C$ 贡献的比重越小。此时 $\mathbf{E}_f$ 由下式计算:

\[ \mathbf{E}_f = (\mathbf{e}\cdot\mathbf{S}_f)\mathbf{e} = |\mathbf{S}_f|\cos\theta\,\mathbf{e} \]

（式 8.68）。

8.6.3 正交修正方法（Orthogonal Correction Approach）

如图 8.11 所示,该方法无论网格非正交程度如何,都让涉及 $\phi_F$ 与 $\phi_C$ 的项与正交网格上的贡献完全相同。为此,$\mathbf{E}_f$ 被定义为 $\mathbf{E}_f = \mathbf{S}_f\mathbf{e}$（式 8.69）。

8.6.4 过松弛方法（Over-Relaxed Approach）

该方法强制让涉及 $\phi_F$ 与 $\phi_C$ 的项的贡献随非正交程度的增加而增加。如图 8.12 所示,这通过令 $\mathbf{T}_f$ 垂直于 $\mathbf{S}_f$ 实现。数学上 $\mathbf{E}_f$ 由下式计算:

\[ \mathbf{E}_f = \frac{\mathbf{S}_f}{\cos\theta}\mathbf{e} = \frac{S_f^2}{\mathbf{S}_f\cos\theta}\mathbf{e} = \frac{\mathbf{S}_f\cdot\mathbf{S}_f}{\mathbf{e}\cdot\mathbf{S}_f}\mathbf{e} \]

（式 8.70）。

总结:非正交网格上,单元面 $f$ 处的扩散通量不能仅写成跨面两节点值的函数,必须加上一项反映非正交性的项——文献中称为"交叉扩散",由下式给出:

\[ (\nabla\phi)_f\cdot\mathbf{T}_f = (\nabla\phi)_f\cdot\mathbf{S}_f - \mathbf{E}_f \cdot \nabla\phi_f = \begin{cases} (\nabla\phi)_f\cdot(\mathbf{n} - \cos\theta\,\mathbf{e})|\mathbf{S}_f| & \text{最小修正} \\ (\nabla\phi)_f\cdot(\mathbf{n} - \mathbf{e})|\mathbf{S}_f| & \text{正交修正} \\ (\nabla\phi)_f\cdot(\mathbf{n} - \mathbf{e}/\cos\theta)|\mathbf{S}_f| & \text{过松弛} \end{cases} \]

（式 8.71）。对正交网格,$\mathbf{n}$ 与 $\mathbf{e}$ 共线,$\mathbf{n}$ 与 $\mathbf{e}$ 之间的夹角 $\theta$ 为零,交叉扩散项为零。当交叉扩散不为零时,因为它无法写成 $\phi_F$ 与 $\phi_C$ 的函数,被作为源项加到单元代数方程的右端。

上述三种方法都正确且满足式 8.59;它们之间的区别仅在于在非正交网格上的精度与稳定性。研究表明,即使在网格高度非正交时,过松弛方法也是最稳定的。但三种方法最终离散扩散项的一般形式是相同的。

8.6.5 交叉扩散项的处理（Treatment of the Cross-Diffusion Term）

交叉扩散项无法写成节点值的函数。因此,采用延后修正（deferred correction）方式处理:用当前梯度场计算其值,作为源项加到代数方程的右端。梯度在主网格点（单元形心）上计算,界面值通过插值得到。

8.6.6 梯度计算（Gradient Computation）

一维或正交多维计算域上,$\phi$ 的梯度可以显式写成形心 $\phi$ 值的函数。在非正交域中,扩散通量的计算则更复杂——梯度的非正交分量不能被线性化写为节点值的函数,而必须移到右端显式求值。这意味着必须对梯度进行求值,以便将非正交贡献纳入离散方程。一种广泛使用的单元梯度计算方法是 Green-Gauss 定理（梯度定理）:对任意闭合体积 $V$（其表面为 $\partial V$）有

\[ \int_V \nabla\phi\,dV = \oint_{\partial V} \phi\,d\mathbf{S} \]

（式 8.72），其中 $d\mathbf{S}$ 是外指微元面积矢量。

为得到该式的离散形式,应用中值定理:体平均梯度与左端体积分的关系为

\[ \overline{\nabla\phi}_V = (1/V)\int_V \nabla\phi\,dV \]

（式 8.73）。结合式 8.72/8.73,图 8.9 所示单元 $C$ 上的体平均梯度为

\[ \overline{\nabla\phi}_C = (1/V_C)\oint_{\partial V_C} \phi_f d\mathbf{S}_f \]

（式 8.74）。

接下来把单元格面上的积分近似为面值乘以面积。于是 $\overline{\nabla\phi}_C$（即 $\nabla\phi_C$）由

\[ \overline{\nabla\phi}_C = (1/V_C)\sum_{f\sim nb(C)} \phi_f \mathbf{S}_f \]

（式 8.75）给出。单元面 $f$ 上的梯度可由共享该面的两单元形心梯度的加权平均得到:

\[ \overline{\nabla\phi}_f = g_C\overline{\nabla\phi}_C + g_F\overline{\nabla\phi}_F \]

（式 8.76），其中 $g_C$、$g_F$ 是与面 $f$ 相对于 $C$、$F$ 形心位置相关的几何插值因子。

8.6.7 非正交网格的代数方程（Algebraic Equation for Non-orthogonal Meshes）

将面积矢量 $\mathbf{S}_f$ 分解为 $\mathbf{E}_f$ 与 $\mathbf{T}_f$,并把相应表达式代入扩散通量的半离散方程,有:

\[ \sum_{f\sim nb(C)} \mathbf{J}_f^{\phi,D}\cdot\mathbf{S}_f = \sum_{f\sim nb(C)} (-\Gamma_\phi\nabla\phi)_f\cdot(\mathbf{E}_f + \mathbf{T}_f) = \underbrace{\sum_{f\sim nb(C)} (-\Gamma_\phi\nabla\phi)_f\cdot\mathbf{E}_f}_{\text{正交可线性化}} + \underbrace{\sum_{f\sim nb(C)} (-\Gamma_\phi\nabla\phi)_f\cdot\mathbf{T}_f}_{\text{非正交不可线性化}} \]

\[ = \sum_{f\sim nb(C)} -\Gamma_{f}\mathbf{E}_f(\phi_F-\phi_C)/d_{CF} + \sum_{f\sim nb(C)} (-\Gamma_\phi\nabla\phi)_f\cdot\mathbf{T}_f \]

\[ = \sum_{f\sim nb(C)} \Gamma_f g_{Difff}(\phi_C - \phi_F) + \sum_{f\sim nb(C)} (-\Gamma_\phi\nabla\phi)_f\cdot\mathbf{T}_f \]

\[ = \left(\sum_{f\sim nb(C)} \text{FluxC}_f\right)\phi_C + \sum_{f\sim nb(C)} \text{FluxF}_f\phi_F + \sum_{f\sim nb(C)} \text{FluxV}_f \]

（式 8.77），其中几何扩散系数 $g_{Difff} = \mathbf{E}_f/d_{CF}$（式 8.78）。

把该形式的扩散通量代入并展开,得到非结构化/结构化非正交网格上的扩散方程最终形式为

\[ a_C\phi_C + \sum_{F=NB(C)} a_F\phi_F = b_C \]

（式 8.79），其中:

\[ a_F = \text{FluxF}_f = -\Gamma_{\phi,f}g_{Difff} \]

\[ a_C = \sum_{f\sim nb(C)} \text{FluxC}_f = -\sum_{f\sim nb(C)} \text{FluxF}_f = \sum_{f\sim nb(C)} \Gamma_f g_{Difff} \]

\[ b_C = Q_{\phi,C}V_C - \sum_{f\sim nb(C)} \text{FluxV}_f = Q_{\phi,C}V_C + \sum_{f\sim nb(C)} (-\Gamma_\phi\nabla\phi)_f\cdot\mathbf{T}_f \]

（式 8.80）。注意右端非正交项的符号变化。

Example 2（例 8.6.7）针对一个多边形单元 $C$ 与其邻居 $F$（图 8.13）,给定 $\phi = x^2+y^2+x^2y^2$,要求 (1) 在 $C$ 形心 $(1.75, 2)$ 处用数值与解析两种方法计算 $\nabla\phi_C$;(2) 计算 $\nabla\phi_F$ 的解析值;(3) 用 $\nabla\phi_C$ 数值与 $\nabla\phi_F$ 解析值插值得 $\nabla\phi_{f_1}$ 近似,与解析值比较;(4) 用三种方法（最小修正、正交修正、过松弛）将 $\nabla\phi_{f_1}\cdot\mathbf{S}_{f_1}$ 表示为 $\phi_C$、$\phi_F$ 的函数。该例完整演示了非正交网格上梯度计算的 Green-Gauss 实施——包括面矢量方向校正、单元形心梯度、界面梯度插值,以及三种分解方案下 $\mathbf{E}_f$、$\mathbf{T}_f$ 的具体数值。

8.6.8 非正交网格的边界条件（Boundary Conditions for Non-orthogonal Grids）

非正交网格上边界条件的处理与正交网格类似,差异主要与非正交扩散贡献有关,简述如下。

8.6.8.1 Dirichlet 边界条件

当 $\phi$ 由用户在边界上指定时（图 8.14）,边界离散过程与正交网格相同。但与内部面一样,需要考虑交叉扩散——只要面积矢量与连接单元形心与边界形心的矢量不共线,就会出现。沿边界面的扩散通量离散为:

\[ \mathbf{J}_b^{\phi,D}\cdot\mathbf{S}_b = -\Gamma_{\phi,b}(\nabla\phi)_b\cdot\mathbf{S}_b = -\Gamma_{\phi,b}(\nabla\phi)_b\cdot(\mathbf{E}_b + \mathbf{T}_b) = -\Gamma_{\phi,b}\frac{\phi_b - \phi_C}{d_{Cb}}\mathbf{E}_b - \Gamma_{\phi,b}(\nabla\phi)_b\cdot\mathbf{T}_b = \text{FluxC}_b\phi_C + \text{FluxV}_b \]

（式 8.81），其中 $\text{FluxC}_b = \Gamma_{\phi,b}g_{Diffb}$、

\[ \text{FluxV}_b = -\Gamma_{\phi,b}g_{Diffb}\phi_b - \Gamma_{\phi,b}(\nabla\phi)_b\cdot\mathbf{T}_b \]

（式 8.82），$g_{Diffb} = \mathbf{E}_b/d_{Cb}$。代入式 8.79 得到修改后的系数: $a_F = \text{FluxF}_f$、

\[ a_C = \sum_{f\sim nb(C)}\text{FluxC}_f + \text{FluxC}_b \]

、

\[ b_C = Q_{\phi,C}V_C - \text{FluxV}_b - \sum_{f\sim nb(C)}\text{FluxV}_f \]

（式 8.83）。

8.6.8.2 Neumann 边界条件

非正交网格上的 Neumann 条件与正交网格相同:用户指定的边界通量直接作为源项加入,边界单元的代数方程由式 8.40 给出,系数由式 8.41 修正。

8.6.8.3 混合边界条件

对混合边界条件（图 8.5）,记对流换热系数为 $h_\infty$、外部 $\phi$ 值为 $\phi_\infty$,边界处扩散通量可写为:

\[ \mathbf{J}_b^{\phi,D}\cdot\mathbf{S}_b = -\Gamma_{\phi,b}\frac{\phi_b - \phi_C}{d_{Cb}}\mathbf{E}_b - \Gamma_{\phi,b}(\nabla\phi)_b\cdot\mathbf{T}_b = -h(\phi_\infty - \phi_b)\mathbf{S}_b \]

（式 8.84）。由此解出 $\phi_b$:

\[ \phi_b = \frac{h\mathbf{S}_b\phi_\infty + (\Gamma_{\phi,b}\mathbf{E}_b/d_{Cb})\phi_C - \Gamma_{\phi,b}(\nabla\phi)_b\cdot\mathbf{T}_b}{h\mathbf{S}_b + \Gamma_{\phi,b}\mathbf{E}_b/d_{Cb}} \]

（式 8.85）。代回式 8.84,通量方程变为:

\[ \mathbf{J}_b^{\phi,D}\cdot\mathbf{S}_b = -\underbrace{\frac{h\mathbf{S}_b(\Gamma_{\phi,b}\mathbf{E}_b/d_{Cb})}{h\mathbf{S}_b + \Gamma_{\phi,b}\mathbf{E}_b/d_{Cb}}}_{a_b}(\phi_\infty - \phi_C) - \underbrace{\frac{h\mathbf{S}_b\Gamma_{\phi,b}(\nabla\phi)_b\cdot\mathbf{T}_b}{h\mathbf{S}_b + \Gamma_{\phi,b}\mathbf{E}_b/d_{Cb}}}_{\mathbf{S}_b^{\phi,CD}} = \text{FluxC}_b\phi_C + \text{FluxV}_b \]

（式 8.86），其中

\[ \text{FluxC}_b = (h\mathbf{S}_b(\Gamma_{\phi,b}\mathbf{E}_b/d_{Cb}))/(h\mathbf{S}_b + \Gamma_{\phi,b}\mathbf{E}_b/d_{Cb}) \]

、

\[ \text{FluxV}_b = -\text{FluxC}_b\phi_\infty - (h\mathbf{S}_b\Gamma_{\phi,b}(\nabla\phi)_b\cdot\mathbf{T}_b)/(h\mathbf{S}_b + \Gamma_{\phi,b}\mathbf{E}_b/d_{Cb}) \]

（式 8.87）。最终边界单元的修改系数为 $a_F = \text{FluxF}_f - \Gamma_{\phi,f}\mathbf{E}_f/d_{Cf}$ 、

\[ a_C = \text{FluxC}_b + \sum_{f\sim nb(C)}\text{FluxC}_f \]

、

\[ b_C = Q_{\phi,C}V_C - \text{FluxV}_b - \sum_{f\sim nb(C)}\text{FluxV}_f \]

（式 8.88）。

8.7 偏斜（Skewness）

在评估一般 $\phi$ 离散方程的诸多项时,常常需要估计其在单元面上的值。面上的估计值应当是整个面的平均值。在离散过程的多个步骤中,都假设变量在节点之间线性变化。如果把这一假设延伸到"沿面方向也线性变化",则任何变量 $\phi$ 的平均值应在面形心处取值。常见做法是使用线性插值剖面,并把面值估计在面与跨面两节点连线的交点处。当网格偏斜时,这条连线不一定穿过面形心[9,10]。

图 8.15 给出示意:线段 $[CF]$ 与面的交点为 $f'$,与面形心 $f$ 不重合。为保持离散方法整体精度为二阶,所有面积分都需要在点 $f$ 处进行。因此需要一个偏斜修正——从插值得到的 $f'$ 值换算到 $f$ 处的值。修正通过对 $\phi$ 在 $f'$ 处做 Taylor 展开实现:

\[ \phi_f = \phi_{f'} + (\nabla\phi)_{f'}\cdot\mathbf{d}_{f'f} \]

（式 8.89），其中 $\mathbf{d}_{f'f}$ 是从交点 $f'$ 指向面形心 $f$ 的矢量。

8.8 各向异性扩散（Anisotropic Diffusion）

至此介绍的扩散方程都假设材料没有传质方向偏好,各方向扩散系数相同,即各向同性介质。当介质的扩散系数与方向有关时,扩散被称为各向异性的[11–16]。如第 3 章所述,某些固体是各向异性的,半离散扩散方程变为

\[ \sum_{f\sim nb(C)} (-\boldsymbol{\kappa}_\phi\cdot\nabla\phi)_f\cdot\mathbf{S}_f = S_{\phi,C}V_C \]

（式 8.90），其中 $\boldsymbol{\kappa}_\phi$ 是二阶对称张量。

假设一般三维情形,经数学处理[17]后,左端可改写为:

\[ (-\boldsymbol{\kappa}_\phi\cdot\nabla\phi)_f\cdot\mathbf{S}_f = -\begin{pmatrix}\kappa_{11} & \kappa_{12} & \kappa_{13} \\ \kappa_{21} & \kappa_{22} & \kappa_{23} \\ \kappa_{31} & \kappa_{32} & \kappa_{33}\end{pmatrix}_f \begin{pmatrix}\partial\phi/\partial x \\ \partial\phi/\partial y \\ \partial\phi/\partial z\end{pmatrix}_f \cdot\mathbf{S}_f \]

（式 8.91）。执行矩阵乘法,式 8.91 化为:

\[ (-\boldsymbol{\kappa}_\phi\cdot\nabla\phi)_f\cdot\mathbf{S}_f = -\begin{pmatrix}\kappa_{11}\partial\phi/\partial x + \kappa_{12}\partial\phi/\partial y + \kappa_{13}\partial\phi/\partial z \\ \kappa_{21}\partial\phi/\partial x + \kappa_{22}\partial\phi/\partial y + \kappa_{23}\partial\phi/\partial z \\ \kappa_{31}\partial\phi/\partial x + \kappa_{32}\partial\phi/\partial y + \kappa_{33}\partial\phi/\partial z\end{pmatrix}_f \cdot \begin{pmatrix}S_x \\ S_y \\ S_z\end{pmatrix}_f \]

（式 8.92）。继续整理得:

\[ (-\boldsymbol{\kappa}_\phi\cdot\nabla\phi)_f\cdot\mathbf{S}_f = -\begin{pmatrix}\partial\phi/\partial x & \partial\phi/\partial y & \partial\phi/\partial z\end{pmatrix}_f \begin{pmatrix}\kappa_{11} & \kappa_{21} & \kappa_{31} \\ \kappa_{12} & \kappa_{22} & \kappa_{32} \\ \kappa_{13} & \kappa_{23} & \kappa_{33}\end{pmatrix}_f \begin{pmatrix}S_x \\ S_y \\ S_z\end{pmatrix}_f = -(\nabla\phi)_f \cdot (\boldsymbol{\kappa}_\phi\cdot\mathbf{S})_f = -(\nabla\phi)_f \cdot \mathbf{S}'_f \]

（式 8.93）。

把式 8.93 代回式 8.90,扩散方程的新形式为

\[ \sum_{f\sim nb(C)} (-\nabla\phi)_f\cdot\mathbf{S}'_f = Q_{\phi,C}V_C \]

（式 8.94）。显然,这种形式下,只要把 $\Gamma_\phi$ 设为 1、$\mathbf{S}_f$ 替换为 $\mathbf{S}'_f$,前文给出的离散流程就直接适用。也就是说,同一套代码即可用于各向同性与各向异性扩散问题。

8.9 迭代求解过程的下松弛（Under-Relaxation of the Iterative Solution Process）

对一般扩散问题,$\Gamma_\phi$ 可能是未知量 $\phi$ 的函数,网格也可能高度非正交、有较大的交叉扩散项——而交叉扩散项采用延后修正处理。因此 $\phi$ 在迭代之间会出现较大变化,随之而来的是大源项、大系数变化,可能引起迭代过程发散。这种发散通常由系数与交叉扩散项所引入的非线性导致——源项强烈依赖于尚未收敛的当前解场。为促进收敛、稳定迭代过程,减缓 $\phi$ 在迭代之间的变化是可取的,这通过一种称为下松弛（under-relaxation）的技术强制实现。

引入下松弛的方法有多种,这里描述其中一种,其他方法在第 14 章给出。在一般离散方程

\[ a_C\phi_C + \sum_{F=NB(C)} a_F\phi_F = b_C \]

（式 8.95）上做推导,式 8.95 可改写为

\[ \phi_C = (-\sum_{F=NB(C)} a_F\phi_F + b_C)/a_C \]

（式 8.96）。

设 $\phi^*_C$ 表示来自上一次迭代的 $\phi_C$ 值。在右端加上再减去 $\phi^*_C$,式 8.96 变为

\[ \phi_C = \phi^*_C + [ (-\sum_{F=NB(C)} a_F\phi_F + b_C)/a_C - \phi^*_C ] \]

（式 8.97）。其中括号内的表达式代表当前迭代对 $\phi_C$ 产生的变化。通过引入松弛因子 $\lambda_\phi$ 修改该变化,有:

\[ \phi_C = \phi^*_C + \lambda_\phi\left[ (-\sum_{F=NB(C)} a_F\phi_F + b_C)/a_C - \phi^*_C \right] \]

（式 8.98），等价形式为:

\[ (a_C/\lambda_\phi)\phi_C + \sum_{F=NB(C)} a_F\phi_F = b_C + [(1-\lambda_\phi)/\lambda_\phi]a_C\phi^*_C \]

（式 8.99）。首先需要注意:在收敛时,$\phi_C$ 与 $\phi^*_C$ 相等,与所使用的松弛因子无关——这由式 8.98 给出,因为收敛时 $\phi_C$ 满足原方程 (式 8.95)。任何松弛方案都应满足这一性质。

根据 $\lambda_\phi$ 的取值,方程可以是下松弛（$0 < \lambda_\phi < 1$）或过松弛（$\lambda_\phi > 1$）。CFD 应用中通常使用下松弛:$\lambda_\phi$ 接近 1 几乎不进行下松弛,接近 0 则产生强下松弛,$\phi_C$ 逐次变化极小。

最佳下松弛因子与问题相关,无统一规则。影响 $\lambda_\phi$ 取值的因素包括:问题的类型、方程组规模（域内网格点数量）、网格间距及其扩展率、所采用的迭代方法等。通常根据经验或初步计算确定 $\lambda_\phi$。此外,无须在整个计算域使用相同的下松弛值,迭代之间也可以变化。

式 8.99 可改写为式 8.95 的形式,此时中心系数变为 $a_C \to a_C/\lambda_\phi$（式 8.100）,源项加上新的一项 $b_C \to b_C + [(1-\lambda_\phi)/\lambda_\phi]a_C\phi^*_C$ （式 8.101）。该松弛方法在稳定非线性问题求解中扮演重要角色。

8.10.1 uFVM

在 uFVM 中,内部面上的扩散项离散由函数 cfdAssembleDiffusionTermInterior 完成,其核心如 Listing 8.1 所示:先由 cfdInterpolateFromElementsToFaces('Average', gamma) 取得面扩散系数 gamma_f,再取面几何量 gDiff_f、面积矢量 Sf、切向分量 Tf、owner/neighbor 单元下标;然后把面扩散系数组装为 $\text{FluxC1}_f = \gamma_f g_{Difff}$、$\text{FluxC2}_f = -\gamma_f g_{Difff}$;非正交项则存入 $\text{FluxV}_f = \gamma_f \cdot \text{dot}(\text{grad}_f, \mathbf{T}_f)$ ;面总通量 $\text{FluxT}_f$ 按式 8.15 形式组装。FLUXC1f 和 FLUXC2f 即正文中的 $\text{FluxC}_f$、$\text{FluxF}_f$,分别是 owner 和 neighbor 单元的系数。

扩散项的设置在 cfdProcessOpenFoamMesh.m 中完成,Listing 8.2 给出 $g_\text{Diff}$ 系数的计算:由 dCF = element2.centroid - element1.centroid 求出形心连线,再计算单位矢量 $\mathbf{e}_{CF} = \mathbf{d}_{CF}/|\mathbf{d}_{CF}|$ 、$\mathbf{E} = \text{face.area}\cdot\mathbf{e}_{CF}$,最后 $g_\text{Diff} = |\mathbf{E}|/|\mathbf{d}_{CF}|$、$\mathbf{T} = \mathbf{S}_f - \mathbf{E}$。

边界条件对方程的影响也需要纳入,Listing 8.3 给出 Dirichlet 条件下边界面扩散项的组装:取出该 patch 的 iBFaces 与对应 iBElements,由 $\text{FLUXC1}_b = \gamma(iBE)g_\text{Diff}$、$\text{FLUXC2}_b = -\gamma(iBE)g_\text{Diff}$、 $\text{FLUXV}_b = -\gamma(iBE)\cdot\text{dot}(\text{grad}(iBE), \mathbf{T}_b)$ 完成。其它边界条件类型的实现可在 cfdAssembleDiffusionTerm.m 中查阅。

对每个内部面与边界面计算出线性化系数后,即可装入全局（稀疏）矩阵——uFVM 在 cfdAssembleIntoGlobalMatrixFaceFluxes 函数中完成。Listing 8.4 展示内部面系数装入稀疏 LHS 矩阵和 RHS 向量的过程:对每个内部面同时更新 owner 单元（$a_C$、邻居项 $a_{NB}$、$b_C$）和 neighbor 单元（$a_C$、$a_{NB}$、$b_C$）的系数;owner 和 neighbor 装配时符号相反,原因是面的面积矢量指向 neighbor 单元而从 owner 单元指出去。

8.10.2 OpenFOAM®

在 OpenFOAM®[18] 中,扩散项可显式或隐式地求值:用 fvc::laplacian(gamma, phi) 显式求值,返回每个单元上 $\phi$ 的 Laplacian 场,加到方程组右端;用 fvm::laplacian(gamma, phi) 隐式求值,返回一个按式 8.19 组装系数的 fvMatrix,加到左端,同时把非正交项加到右端。Laplacian 算子的定义位于 $FOAM_SRC/finiteVolume/finiteVolume/laplacianSchemes/gaussLaplacianScheme 目录下的 gaussLaplacianScheme.C/H 和 gaussLaplacianSchemes.C 文件中。

fvm::laplacian 的实现（Listing 8.5）先构造 $\gamma\cdot|\mathbf{S}_f|$ 场 gammaMagSf,然后调用 fvmLaplacianUncorrected 装配出式 8.18 的基本形式;若所选 snGrad 方案是 "corrected" 的,则额外装配面通量修正 faceFluxCorrectionPtr 并把它通过 fvc::div 加到源项,实现非正交贡献。

fvmLaplacianUncorrected（Listing 8.6）创建 fvMatrix 对象并按 Laplacian 算子对称性只填上三角 fvm.upper();边界条件装配时,对 coupled 边界使用 gradientInternalCoeffs(pDeltaCoeffs) 与 gradientBoundaryCoeffs(pDeltaCoeffs),否则使用 gradientInternalCoeffs()、gradientBoundaryCoeffs()。

OpenFOAM® 中装配是直接在全局系数中进行的——如第 5、6、7 章所述,存储在 fvm.upper()、fvm.lower()、fvm.diag() 三个数组中。主要离散部分在 fvmLaplacianUncorrected 中定义,先创建 fvMatrix 对象,再把额外对角系数（这里依赖 Laplacian 算子返回对称矩阵的性质）填到上三角。deltaCoeffs.internalField() 即 $g_\text{Diff}$ 场,gamma 是扩散系数场。对角系数在 fvm.negSumDiag() 中按式 8.18 由上、下三角系数的负和装配得到（Listing 8.7 给出 lduMatrix::negSumDiag() 的实现）。

边界条件在 OpenFOAM® 中严格按式 8.36、8.41 实现,边界系数存储在 internalCoeffs（$\text{FluxC}_b$）和 boundaryCoeffs（$\text{FluxV}_b$）中（如第 7.6 节所述）。fvmUncorrected 只包含正交离散;非正交贡献按 Listing 8.8 通过 faceFluxCorrectionPtr 与 fvc::div 加到源项,正好实现式 8.77 的最后一项。snGrad 类把非正交项包装在 correction 函数中（Listing 8.9 给出 correctedSnGrad::correction 与 fullGradCorrection）,后者通过 mesh.nonOrthCorrectionVectors()（即式 8.66 中的 $\mathbf{T}$ 矢量）实现非正交修正。

Laplacian 的离散类型在 fvSchemes 文件的 system 目录中指定(Listing 8.10):laplacian(gamma, phi) Gauss linear corrected;——其中 Gauss（唯一选项）定义标准 Gauss 离散,得到式 8.18;linear 指计算 $\gamma$ 在面上插值的方式;corrected 描述所用的非正交修正类型。OpenFOAM® 中边界条件实现的更多细节将在后续章节给出。

8.11 小结（Closure）

本章描述了扩散方程的离散。讨论了若干相关议题:正交与非正交网格系统的使用、边界条件的实现以及下/过松弛。下一章将集中讨论梯度场的计算。

8.12 习题（Exercises）

(Exercises skipped per Constraint #13.)

本章个人批注

第 8 章是 Moukalled 等人 FVM 著作的核心方法论章节,主题是扩散项的空间离散。从章节结构上看,作者遵循了一个"由简入繁"的清晰递进路径:8.1 节先在二维 Cartesian 网格上把线性剖面 + 二阶中心差分这一最基本的扩散离散范式呈现出来;8.2 节用零和规则与异号规则两条原则对离散方程的物理意义进行"事后审查",这其实是回到 Patankar《Numerical Heat Transfer and Fluid Flow》中Scarborough 准则的两个核心要求;8.3 节把 8.1 节的范式推广到四类边界条件;8.4 节通过调和平均替代算术平均,把单 $\Gamma$ 范式扩展到间断 $\Gamma$;8.5–8.6 节把整套离散从 Cartesian 推广到正交旋转网格、再推广到任意非正交结构化/非结构化网格;8.7 节针对偏斜网格的面形心与连心线交点不一致给出 Taylor 展开修正;8.8 节把标量 $\Gamma$ 升级为二阶张量 $\boldsymbol{\kappa}$;8.9 节则处理非线性 $\Gamma(\phi)$ 与大交叉扩散带来的迭代发散问题。

对我自己而言,这一章最值得反复咀嚼的细节是 8.4 节"调和平均"那一段。从偏微分方程的角度看,$\Gamma$ 在 $\phi$ 的 Laplacian 里是一个系数,在单元交界面上取算术平均还是调和平均似乎只是工程取舍。但作者把这一点和一维稳态、无源的情形严格连接起来:对一维间断系数问题,调和平均是精确的,算术平均反而是近似的,且对大比例系数差(如 $k_\text{材料2}/k_\text{材料1} = 10^5$)在 $C/E$ 之间的 $1\times 1$ 单元上,算术平均会给出严重的误差(这是他们 Example 1 几乎全用调和平均的根本原因)。这一论据的简洁性让人印象深刻:不是"调和平均在某种意义下更好",而是"在它能精确的一维问题里它是精确的"。再加上"共轭传热界面无需特殊处理"这一工程优势,调和平均在 FVM 教材中几乎成为标配——这在 Moukalled 这本书里被反复强化。

8.6 节非正交网格的处理则体现了 FVM 与 FDM/FEM 的一个根本区别:FVM 的所有推导都以面积矢量 $\mathbf{S}_f$ 为核心,而非以节点/自由度为核心。因此当网格非正交时,$\mathbf{S}_f$ 与 $\mathbf{d}_{CF}$ 不再共线,作者把 $\mathbf{S}_f$分解为沿 $\mathbf{d}_{CF}$ 的 $\mathbf{E}_f$ 与切向 $\mathbf{T}_f$,从而把面通量写成"可线性化(只涉及 $\phi_C, \phi_F$)"+ "延后修正(用旧梯度场显式求值)"两部分。这一分解是任意的——三种方法(最小修正、正交修正、过松弛)都是合法的,区别只在于精度与稳定性;过松弛方法在高度非正交网格上稳定性最好。这种"分解选择 + 延后修正"的思想在 OpenFOAM® 的 nonOrthogonalCorrection 框架里被完整保留(8.10.2 节)。

8.7 节的偏斜修正与 8.8 节的各向异性扩散则展示了 FVM 的"通用性"。前者解决"面值应在面形心 $f$ 处取值,但线性插值给出的是连心线交点 $f'$ 处的值"这一几何不一致;后者则通过把"$\Gamma_\phi$ 设为 1、$\mathbf{S}_f$ 替换为 $\boldsymbol{\kappa}_\phi\cdot\mathbf{S}_f$"实现同一套代码求解各向同性与各向异性问题——这是 Moukalled 等人在 8.8 节末尾强调的核心设计思想。

最后,8.9 节的下松弛技术是非线性迭代的标准武器,作者清楚指出"最佳 $\lambda_\phi$ 与问题相关、无统一规则",并提示影响因素——这在 CFD 实际应用中确实如此,没有放之四海皆准的松弛因子。

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

从全书的章节顺序看,第 7 章聚焦"网格数据结构与 OpenFOAM®/uFVM 实现",为后续的离散化方法提供几何与拓扑载体;第 8 章转向扩散项的空间离散——这是 FVM 三大物理项(扩散、对流、瞬态)中最基本也是最稳定的一项;第 9 章紧接着将讨论梯度场的计算(8.6.6 节中已经用到了 Green-Gauss 梯度定理,8.7 节的偏斜修正涉及 $\nabla\phi$ 在 $f'$ 处的 Taylor 展开,8.8 节的各向异性扩散也涉及梯度项的处理),可以说第 9 章是第 8 章的自然延续;第 10 章再转入代数方程组的求解(8.9 节的下松弛是其中的关键预处理);第 11 章则开始对流项的离散——这与本章的扩散项形成对照(作者在前言里也明确指出"对流与扩散代表两种不同的物理现象,必须分别处理")。所以第 8 章在全书中处于承上启下的枢纽位置:向上承接网格,向下开启各项离散。