第 5 章：有限体积法（The Finite Volume Method）
5.1 Introduction
本章以抽象方式介绍有限体积法（FVM），点明其在计算流体力学（CFD）中的流行源自高度的离散化灵活性。FVM 虽比有限差分法和有限元法出现得晚，却在 1970 年代初经 Imperial College 的 Spalding 教授团队（Patankar、Gosman、Runchal 等）系统发展后，承担起流体流动与相关输运现象的数值模拟任务。其灵活性与流行度来自两个事实：一是离散化直接在物理空间进行，物理坐标系与计算坐标系之间无需变换；二是采用 collocated arrangement（变量同位布置），令该方法适合在复杂几何中求解流动问题。这些发展使 FVM 在保留数学表述简洁性的同时，把应用范围扩展到各类工程问题。FVM 数值方法的一个重要特征是其"镜照物理"的能力——数值格式本身反映守恒律的积分性质以及被离散项的特性。本节末尾预告：下文先推导一般标量方程的半离散化形式，再讨论离散化方法应具备的性质与若干指导原则，章末讨论若干与 FVM 相关的议题；把半离散化方程转化为代数方程的工作将是后续若干章的主题。
5.2 The Semi-Discretized Equation
本节进入 FVM 离散化的第 1 步：先把守恒方程在网格所分出的单元（有限体积）上积分，再用 Gauss 定理把对流项与扩散项的体积分转化为面积分；随后用积分点（integration point, ip）把面积分与体积分近似为离散量。为便于理解后文高级离散格式，作者以二维输运问题为例演示这一流程。第 3 章给出一般标量变量 ϕ 的守恒方程为
即 transient 项 + 对流项 = 扩散项 + 源项。稳态形式把 transient 项略去，得
。把此式在图 5.1 所示的单元 C 上积分，得到体积分形式；再用散度定理把对流项与扩散项的体积分换成面积分，得
。该式中黑体表示向量，\((\cdot)\) 为点积，\(Q_\phi\) 表示源/汇项，\(\mathbf{S}\) 为面向量，\(\mathbf{v}\) 为速度向量，\(\phi\) 为守恒量，\(\partial V_C\) 为单元 C 上的面积分。
5.2.1 Flux Integration Over Element Faces
把对流通量与扩散通量分别记为 \(\mathbf{J}_{\phi,C} = \rho \mathbf{v} \phi\) 与 \(\mathbf{J}_{\phi,D} = \Gamma_\phi \nabla \phi\)，二者和为总通量 \(\mathbf{J}_\phi = \mathbf{J}_{\phi,C} + \mathbf{J}_{\phi,D}\) 。把单元 C 上的面积分替换为 C 各面上通量之和，可把对流、扩散、总通量分别写成
等三式。这一转换将"面通量在面上积分"化为"在单元各面 f 上逐面求值"，是该方法具有守恒性的根本原因之一。
要进一步离散，需把每个面上的面积分与源项体积分求值。采用 Gauss 积分，可把面 f 上的积分写成
，其中 ip 指积分点，ip(f) 为面 f 上积分点个数。源项体积分的高斯形式为
，其精度同样取决于积分点个数 ip 与权重 \(\omega_{\text{ip}}\)。图 5.3 展示了体积分的若干选项：一积分点 Gauss 积分（\(\text{ip} = \omega_{\text{ip}} = 1\)，积分点位于单元质心，二维三维皆可用，二阶精度）；二维四积分点 Gauss 积分（\(\text{ip} = 4\)，权重为单元一维权重的乘积，\(\omega_{\text{ip}} = 1/4\)，由一维节点构造二维节点）；九积分点 Gauss 积分（\(\text{ip} = 9\)）。精度随积分点增加而提高，计算成本也相应上升。
精度取决于积分点个数 ip 与权重 \(\omega_{\text{ip}}\)。图 5.2 展示面积分的若干选项：（a）简单均值积分（梯形法则），面质心处一积分点，\(\text{ip} = \omega_{\text{ip}} = 1\)，二阶精度，二三维皆可用；（b）二维下两积分点（\(\text{ip} = 2\)），位置为 \(\xi_1 = (3 - \sqrt{3})/6\) 与 \(\xi_2 = (3 + \sqrt{3})/6\)（\(\xi\) 为沿面归一化距离），权重 \(\omega_1 = \omega_2 = 1/2\)，三阶精度；（c）三积分点（\(\text{ip} = 3\)），位置为 \(\xi_1 = (5 - \sqrt{15})/10\)、\(\xi_2 = 1/2\)、\(\xi_3 = (5 + \sqrt{15})/10\)，权重 \(\omega_1 = 5/18\)、\(\omega_2 = 4/9\)、\(\omega_3 = 5/18\)。计算成本随积分点个数上升。以 ip(f) 表示面 f 上积分点个数，对流项与扩散项的通用离散形式分别为
与
。
5.2.2 Source Term Volume Integration
源项采用体积分的高斯形式
，其精度同样取决于积分点个数 ip 与权重 \(\omega_{\text{ip}}\)。图 5.3 展示体积分的若干选项：一积分点 Gauss 积分（\(\text{ip} = \omega_{\text{ip}} = 1\)，积分点位于单元质心，二阶精度，二三维皆可用）；二维四积分点 Gauss 积分（\(\text{ip} = 4\)，权重为单元一维权重的乘积，\(\omega_{\text{ip}} = 1/4\)，由一维节点构造二维节点）；九积分点 Gauss 积分（\(\text{ip} = 9\)）。精度随积分点增加而提高，计算成本也相应上升。
5.2.3 The Discrete Conservation Equation for One Integration Point
虽然原则上上述各项可以用任意多个积分点离散，FVM 的惯例是只用一个积分点以取得二阶精度——这是精度、灵活性、方法简洁性与计算成本之间的良好折中。采用中点积分近似，图 5.4 所示单元 C 的稳态半离散有限体积方程最终简化为
。离散化第 2 步的目标是把此方程转化为代数方程，即把面通量与体通量用相邻单元质心处的变量值表示。这一对通量的线性化是第 2 步离散化的核心。
5.2.4 Flux Linearization
如图 5.2a 所示，面通量可拆分为线性部分与非线性部分：线性部分取决于跨面的两个节点处 \(\phi\) 值（即 \(\phi_C\) 与 \(\phi_F\)），非线性部分为无法用 \(\phi_C\)、\(\phi_F\) 表达的部分。所得方程为
，三项分别对应单元 C 与 F 的线性化系数与非线性项。\(\text{Flux}^C_f\)、\(\text{Flux}^F_f\)、\(\text{Flux}^V_f\) 的取值取决于被离散项及其离散格式。
将上式代入式 (5.15) 的左端，并对所有单元面重复，得
。体通量的线性化（图 5.3a）写成 \(\phi\) 在单元节点值 \(\phi_C\) 的线性函数： \(Q_{\phi,C} V_C = \text{Flux}^T_C = \text{Flux}^C \phi_C + \text{Flux}^V\) 。在源项为常数的情况下，体通量（即式 (5.15) 右端）简化为 \(\text{Flux}^C = 0\)、\(\text{Flux}^V = Q_{\phi,C} V_C\)。
把上述线性化形式代回式 (5.15)，得到所求代数关系：
，其中方程系数与通量线性化系数的关系为
、\(a_F = \text{Flux}^F_f\)、
。
Example 1 演示了对流项在速度场为正向、用迎风格式 \(\phi_f = \phi_C\) 下的线性化系数：\(J_{\phi,f} = (\rho \mathbf{v} \phi)_f\)，故
。带入 \(\phi_f = \phi_C\) 即得 \(\text{Flux}^C_f = \dot{m}_f\)、\(\text{Flux}^F_f = 0\)、\(\text{Flux}^V_f = 0\)，总通量为 \(\text{Flux}^T_f = \dot{m}_f \phi_C + 0 \cdot \phi_F + 0\) 。
5.3 Boundary Conditions
域边界处面通量的求值通常不需要假设分布，而直接代入。边界条件类型繁多，但一般标量最常用的两类是 Dirichlet 与 Neumann 边界条件，在数学上分别对应"指定值"（第一类）与"指定通量"（第二类）。
5.3.1 Value Specified (Dirichlet Boundary Condition)
考虑标量 \(\phi\) 通过 inlet 被对流输送的情形。若扩散可忽略，边界条件为 \(\phi_b = \phi_{b,\text{specified}}\)。对图 5.5 所示的边界面，边界通量用已知的 \(\phi_b\) 值直接求值：
，故 \(\text{Flux}^C_b = 0\)、 \(\text{Flux}^V_b = \dot{m}_f \phi_{b,\text{specified}}\) 。
5.3.2 Flux Specified (Neumann Boundary Condition)
考虑图 5.6 所示情形：单元 C 的边界面 b 物理上是一面墙，\(\phi\) 的通量被指定。数学上等价于写
。其中 \(q_{b,\text{specified}}\) 是用户指定的已知量，表示单位面积通量。故 \(\text{Flux}^C_b = 0\)、\(\text{Flux}^V_b = q_{b,\text{specified}} S_b\)。以上两类边界条件的处理方法及其它边界条件将在后续章节详述。
5.4 Order of Accuracy
前文已讨论，面通量与体源项按中值方法——即面质心（中点法则）与单元质心处的值——求值。这一处理加上对 \(\phi\) 在空间上分布的假设（即 \(\phi = \phi(x)\)）决定了离散化过程的精度。本章方法假设 \(\phi\) 在空间上线性变化：
，其中 \(\phi_C = \phi(\mathbf{x}_C)\)。
5.4.1 Spatial Variation Approximation
变量 \(\phi = \phi(x)\) 在图 5.7 所示单元内的空间变化可用围绕 \(\mathbf{x}_C\) 的 Taylor 展开描述为
。式中 \((\mathbf{x} - \mathbf{x}_C)^n\) 表示向量 \((\mathbf{x} - \mathbf{x}_C)\) 与自身的 \(n\) 次张量积，所得为 \(n\) 阶张量；运算符 \((:)\) 为两个 2 阶张量的内积，\((\vdots\vdots)\) 为两个 3 阶张量的内积，一般地 \(\overbrace{\vdots \dots \vdots}^{(n-1)\text{ times}}\) 是两个 \(n\) 阶张量的内积，结果皆为标量。把式 (5.27) 给出的假设分布与 Taylor 展开式 (5.28) 相比较，可知误差正比于 \(|(\mathbf{x} - \mathbf{x}_C)^2|\)，对应二阶空间精度。
5.4.2 Mean Value Approximation
中值近似的精度可通过把变量 \(\phi(x)\) 在图 5.8 所示质心为 C 的单元上积分得到：
把积分拆开后，第二项在 C 为质心的情况下为零，略去末项则引入二阶误差，故中值近似是二阶精度的。
单元 C 在面 f（图 5.9）上的对流通量计算为
把 \(\phi\) 在面质心处 Taylor 展开后逐项积分，
因 \(\mathbf{x}_f\) 为面质心而等于零，余项为 \(O(|(\mathbf{x} - \mathbf{x}_f|^2)\)，故对流通量二阶精度。扩散通量同理：
综上，二阶方法要求面积分与单元内 \(\phi\) 的假设分布均为二阶精度。更高的精度可通过提升面积分精度或 \(\phi\) 的假设分布精度来实现。Lilek 和 Peric 在二维下提出了一种用 Simpson 法则实现的四阶精度面积分：
。对图 5.10 所示面 f，\(\phi\) 需要在三个积分点（面质心 \(\text{ip}_2\) 与面两端顶点 \(\text{ip}_1\)、\(\text{ip}_3\)）处取值。要得到形式上的四阶精度离散，单元内 \(\phi\) 的假设分布也应四阶精度，即
。此时梯度计算至少要二阶精度，Hermite 至少要一阶精度。
5.5 Transient Semi-Discretized Equation
对非稳态问题，需保留式 (5.1) 的时间项，除体积分外还要做时间积分。积分后的方程为
对固定网格（单元体积与面积不随时间变化），第一项可化为
，其中
。代回原式得简化方程
再用中点法则即得
往后推进需要对变量随时间变化的方式做某种假设。
5.6 Properties of the Discretized Equations
随单元尺寸趋于零，数值解应趋近一般守恒方程（如式 (5.2)）的精确解，与求单元 \(\phi\) 值所用的插值分布无关。但因为使用了有限体积，离散方程必须具备若干性质以确保解场有意义。这些性质下文逐一讨论。
5.6.1 Conservation
从物理角度看，被输运变量（一般是守恒量，如质量、能量等）必须在离散后的求解域内仍然守恒，否则结果可能失真。这一性质为 FVM 所固有，因为单元面上的通量积分基于共享该面的两个单元的值。任何两个单元共有的面上，一个单元离开该面的通量恰等于另一单元通过同一面进入的通量——即大小相等方向相反（图 5.11）。具备此性质的方法称为守恒。
5.6.2 Accuracy
精度指数值解与精确解的接近程度。但在多数情形下待求问题的精确解未知，因此无法直接比较。替代方案是用截断误差来度量精度。上述各项在第 1 步离散化中的误差为 \(O(|(\mathbf{x} - \mathbf{x}_f)|^2)\)，即二阶精度。这意味着网格点数加倍时离散误差将降为原来的 1/4。离散格式的截断误差是其各项中最大的那一项。离散误差并不给出某一网格下误差的具体值，但能告诉我们误差随网格加密减小的速度——误差阶数越高，随网格加密衰减越快。
5.6.3 Convergence
这里处理的守恒方程是非线性的，因此需要迭代求解。以一个初始猜测出发，反复施加求解算法，把上一轮迭代末的解作为下一轮的初始猜测。理想情况下，若解不再随迭代推进而变化即称收敛；实际中则以连续两次迭代之间的变化低于一个微小量 \(\epsilon\) 作为收敛判据。一般而言，收敛被用来指任何方法求得解的过程。有时"收敛"也被用来指获得网格无关解，即解不再随进一步网格加密而变化。
5.6.4 Consistency
若某一代数方程在某一点处，当时间步与网格间距趋于零（即离散误差趋于零）时，数值解趋于原偏微分方程的精确解，则称该代数方程与该偏微分方程相容（consistent）。为此，离散误差应为 \(\Delta t\) 和/或 \(\Delta x\) 的幂函数。若离散误差用 \(\Delta x / \Delta t\) 表示，则相容性要求 \(\Delta x\) 趋于零的速度快于 \(\Delta t\)。
5.6.5 Stability
稳定性指被迭代求解器解析的离散方程组的特性——即所得代数方程组在各种初、边值条件下能否被求解。从这一意义上，稳定性更属于方程组的性质而非离散化过程本身的性质。前章已述，线性方程组稳定且收敛的充分条件是满足 Scarborough 准则，即系数矩阵对角占优。对瞬态问题，稳定的数值格式保证解中的误差随时间推进保持有界。显式或隐式格式的选用直接影响数值方法的稳定性：显式方法的稳定性靠限制时间步长来保证；隐式方法的稳定性则可通过亚松弛（欠松弛因子或伪瞬态方法）增强。
5.6.6 Economy
经济性是 CFD 代码开发与应用中的重要考量。它指计算真实流动问题所需的时间不应过于昂贵。
5.6.7 Transportiveness
流体输运所具有的方向性是众所周知的，它体现在一般标量输运方程（式 (5.2)）类型的变化上——在某些条件下方程会变为双曲型。其对有限体积方程的影响（图 5.12）可这样说明：在均匀速度与扩散率的流场中，若单元 C 内有恒定源，则常 \(\phi\) 标量等值线的形状将受对流与扩散强度之比，即 Péclet 数的影响，其定义为
。\(\text{Pe} = 0\) 时 \(\phi\) 的输运由扩散主导，扩散行为为椭圆型，此时 \(\phi\) 的等值线为圆，单元 C 处的 \(\phi\) 受周围节点 W 与 E 的共同影响。随对流效应增强（\(\text{Pe}\) 增大），圆形等值线变为椭圆，影响 C 处 \(\phi\) 的区域向流动方向偏移。因此对高 Pe 流动，C 节点处的变化对上游节点影响微弱或不产生影响，而对下游节点影响显著。所选离散格式若不满足这一要求，可能引发不稳定（即非物理振荡）。
5.6.8 Boundedness of the Interpolation Profile
仅满足守恒性并不能保证原偏微分方程的其它重要性质在离散方程中得到保持。例如从物理考虑出发，在无源情形下，域内守恒变量 \(\phi\) 的值应被定义域边界处的值所界定（即 bounded）。离散方程
在节点 \(\phi_C\) 值满足以下约束时满足此要求： \(\min_i (\phi_{F_i}) \le \phi_C \le \max_i (\phi_{F_i})\) ，其中 \(F_i\) 为 C 的第 \(i\) 个邻居，\(N\) 为邻居总数。这可以通过审慎地控制离散格式与其线性化来实现，将在后续章节详述。
5.7 Variable Arrangement
虽然 cell-centered 变量排布是 OpenFOAM® 所采用且在 FVM 中通常更受偏好（图 5.13a），vertex-centered [15]（图 5.13b）变量排布方法也曾被使用。vertex-centered 排布又分为两类，分别产生重叠单元与对偶网格；因为对偶网格法更受欢迎，下面把它作为 vertex-centered 方法的代表加以介绍。
5.7.1 Vertex-Centered FVM
在 vertex-centered 排布中，流动变量存储于顶点（或网格点），单元通过引入对偶网格与对偶单元 [16]（图 5.14）的概念，围绕变量位置构造。引入对偶网格后，可用多种方式在网格点周围构造单元：二维下的一种方式（图 5.14a）是把共享该网格点的各单元质心连接起来；图 5.14b 中的另一种方式是把周围单元的质心与它们各面的质心相连。三维情形下同样的构造也可用——围绕网格点的单元由周围单元的质心、面质心、棱质心适当连接而成，所得单元互不重叠。
vertex-centered 变量排布允许以顶点变量显式地在单元上定义分布。在此情形下变量代表点值，单元内的变化可用形函数或插值分布计算。这一方法对各种网格拓扑都能精确求取面通量，但因顶点不一定位于单元质心，单元积分的精度较低。此外，对偶网格的生成增加了矩阵规模，存储需求上升。进一步，在 cell-vertex 格式中边界条件的处理（图 5.14b）需要额外处理，以确保在多个网格块共享的边界点处获得一致的解。其主要缺点是需要网格所基于的单元类型集合能定义形函数。
对偶单元的 vertex-centered 格式在固壁边界处还有另一缺点：单元只有部分形成，存储变量值的节点位于壁面（图 5.14b）。在常规单元中，面通量积分得到的残差存于单元内、存放变量的主点；但此情形下残差将存于壁面边界处，与 cell-centered 格式相比，在该处引入了离散误差的增大。在尖锐棱角与分支切割处也可能出现额外复杂性。
5.7.2 Cell-Centered FVM
cell-centered 变量排布是目前 FVM 中最常用的排布方式。它把变量及其相关量存于网格单元（或单元）的质心，单元与离散单元一致。一般而言，由于所有量都在单元与面质心处求值，变量值与其均值之差为 \(O(\Delta x^2)\)，故该方法是二阶精度的。单元内的变化可用 Taylor 展开重构。cell-centered 形式的另一优势是允许使用一般多边形单元而无需预定义形函数，从而便于实施完整的多重网格策略。但该方法有两个重要缺点：对非接续单元的处理方式，以及扩散项在非正交单元上的离散方式。前者影响精度，后者影响鲁棒性，二者都受网格质量影响。
考虑图 5.15 所示的两单元排布：在 C 和 F 处定义值的任何均值都会定义在 \(f'\)（连接两质心的线段中点）而非 f（面质心）。因此任何使用此插值的离散格式都不具有 \(O(\Delta x^2)\) 的精度。
对 cell-centered 格式，离散误差强依赖于网格的光滑度。图 5.16 展示了这种情况下的结果：物理情形（图 5.16a）为两个水平横截面为菱形的圆柱之间的环形空腔，内圆柱维持均匀高温 \(T_h\)，外圆柱维持低温 \(T_c\)；温差在封闭腔内造成密度变化，形成浮升力，建立起流场。用图 5.16a 所示网格系统，用 FVM 数值求解速度与温度分布，所得等温线如图 5.16b 所示。仔细观察图 5.16b 可在水平中线附近看到一些小的扭结（kinks）。扭结处局部放大后网格与等温线如图 5.16c、d 所示。扭结的出现是因为网格存在斜率不连续（slope discontinuity），由此造成的离散误差无论如何加密网格都不能减小。这一零阶误差在 cell-vertex 格式中不会出现。尽管如此，对充分光滑的网格，cell-centered 排布仍可达到二阶或更高精度。
另一个影响 cell-centered FVM 的问题是非正交性在扩散项离散中的处理。这将在第 8 章详述，此处只做简要说明。扩散项离散（图 15.17）中，单元 C 与 F 质心连线方向的单位向量 \(\mathbf{e}\) 与 C、F 共享面法向单位向量 \(\mathbf{n}\) 之间的夹角 \(\theta\) 影响隐式程度，进而影响扩散项离散方法的鲁棒性。对离散而言，扩散项通常写为
。\(\theta\) 越大，显式项越大，离散方法的鲁棒性越差。
综合而言：对稳态计算，vertex-centered 对偶单元格式与 cell-centered 格式在域内部数值上非常相似。只有在扭曲网格上 vertex-centered 格式才表现优于 cell-centered 格式；其它情形下用 cell-centered 排布更利于在计算机代码中直接实现。
5.8 Implicit Versus Explicit Numerical Methods
如 5.1 节所述，一旦积分点个数与线性化类型确定，cell-centered 有限体积离散方法所得方程组以待求因变量在单元质心处的值为未知数。这些未知数的组织与求解方式决定了所用计算方法的类别。数值求解格式常被称为显式或隐式。
显式数值方法中，因变量通过已知值直接算出——任何离散算子都能直接基于实际变量值计算。隐式数值方法中，因变量被视为未知数，组装为耦合方程组后用专门的直接法或迭代法求解算法求解。CFD 中的守恒方程是非线性的，隐式方法比显式方法更受偏好。完成离散化与线性化后，最后一步是求解所得代数方程组。用直接方程求解器（参见前章）不可行；迭代法更经济，是最常用甚至唯一可用的方法。与直接法不同，迭代法从初始猜测出发，用本轮迭代末的结果作为下一轮迭代的初始猜测，不断逼近解。中间解不断接近最终解时称序列收敛，否则称发散。这一过程是通用的，无论是对瞬态问题的一个时间步求解还是对稳态问题的最终解都适用。事实上，用迭代法获得稳态解比用时间推进直到稳态更经济。更多关于矩阵迭代求解器的内容将在第 10 章讨论。
5.9 The Mesh Support
至此 FVM 数值方法基础及其性质已建立，本节概述所需的网格支撑信息。前述推导中隐式地假设了若干几何与拓扑信息可用。虽然 FVM 网格的描述将是下两章的主题，但根据目前已述内容，有必要列出其主要特征的简要清单。
FVM 作为一种离散化数值技术的实施，要求网格支撑提供从几何到拓扑的一系列信息。对单元，需要的信息包括：单元的索引、质心、邻接面列表、相邻单元列表。对面，需要：面的索引、质心、面向量、相邻单元列表（内部面 2 个，边界 1 个）以及定义面的顶点列表。此外还需要关于计算域边界的信息，即定义每个边界 patch 的边界面。
另一个要解决的问题是单元面法向的方向性。前述推导中假设单元所有面的法向都指向外侧，但计算域内所有面的法向不可能都指向外侧。一般而言任一单元总有部分面法向朝外，部分朝内。这意味着对面符号（face sign）应做恰当处理。
以上各项的定义与使用方式将是下章的主题。
5.10 Computational Pointers
本书将用基于 Matlab® 的有限体积 CFD 教学代码 uFVM 与可求解工业级问题的开源有限体积代码 OpenFOAM® [17] 来举例说明实现细节并落实数值流程。除另作说明外，关于这两套代码所采用实现技术、方法、数据结构的评注，将出现在后续每章的"Computational Pointers"小节中。此外，书中还为读者提供了若干 uFVM 测试用例（testDiffusion、testAdvection、testFlow、testSource 等），可从前文提到的网站地址下载。
5.10.1 uFVM
uFVM 是为学术目的开发的非结构三维有限体积代码，用 Matlab® 编写，因此便于读者按需剖解。Matlab® 环境允许用户在运行算例时随时展开代码所用的各种数据结构。此外，因其面向学术与教学用途，编写实现细节时把代码清晰度列为高优先级。尽管如此，其数值与算法在许多方面与工业 CFD 代码所用相似，故对有志开发 CFD 代码的人来说 uFVM 仍是有用的指引。
uFVM 中，离散方程组以矩阵形式存储，每行代表一个单元的离散方程，非对角系数代表相邻单元的相互影响。该矩阵的存储形式一般利用"每行的非零元素个数很少"这一事实，即每行的非零元素个数等于该行对应单元的邻居数。以图 5.18 所示算例为例，网格含 7 个单元，单元 3 有 4 个邻居，故单元 3 的离散方程除对角系数外还含 4 个系数。对大型网格，很多非对角元素为零的原因由此可知。在 uFVM、OpenFOAM® 乃至所有 CFD 代码中，都采用稀疏矩阵存储所得方程组的系数。在 uFVM 中，矩阵 A 存为数组的数组，每行的首元素为对角元素，其余元素存放非对角系数。右端项存于单独的数组 B 中。图 5.18 展示了这一结构。
5.10.2 OpenFOAM®
OpenFOAM®（Open source Field Operation And Manipulation）是一个面向对象的 C++ 框架，可用于构建以有限体积离散为中心的多种连续介质力学求解器。OpenFOAM® 还自带若干即用求解器、工具与应用。其核心是一组允许程序员以高编码层次操作网格、几何与离散技术的对象类。表 5.1、表 5.2、表 5.3 与表 5.4 列出主要 OpenFOAM® 类及其功能——这些类构成基于 OpenFOAM® 的应用与工具开发的基础，使程序员在广泛复用代码的同时构建各种算法。
OpenFOAM® 的另一特征是采用运算符重载（operator overloading），让算法能以自然方式表达。例如一般标量 \(\phi\) 的输运方程
即 unsteady 项 + convection 项 = diffusion 项 + source/sink 项。在 OpenFOAM® 中基本写为 Listing 5.1 所示脚本。
其中 fvm::div 算子以控制单元面定义的通量系数场为输入、以单元质心定义的变量场 phi 为输入，返回包含 LHS 矩阵与 RHS 源的方程组，二者代表 convection 算子的离散。这些 LHS 矩阵与 RHS 向量为各算子分别生成，再按需相加减，得到域内每个单元上的最终方程组（图 5.19）。
命名空间 fvm:: 与 fvc:: 分别允许隐式与显式地求值各种算子。显式算子 fvc（finite volume calculus）基于实际场值返回等价场；例如 fvc::div(\(\phi\)) 返回一个等价 geometricField，其中每个单元含 \(\phi\) 散度值。隐式算子 fvm 用系数矩阵定义隐式有限体积离散；例如 fvm::laplacian(\(\phi\)) 返回一个 fvMatrix，其中所有系数都基于 \(\phi\) 的拉普拉斯算子的有限体积离散。fvm:: 与 fvc:: 算子的作用是构造网格每个单元上代表式 (5.43) 离散形式的方程组 LHS 与 RHS。离散过程产出的方程组可用图 5.19 所示矩阵形式表示。
在 OpenFOAM® 中，系数矩阵存储于类 lduMatrix 与专用类 fvMatrix。系数采用 ldu 稀疏存储格式：所有系数存于三个主数组 diagonal、upper、lower 中。对角系数存于数组内，可通过 diag() 函数访问，维度等于计算域内单元数；上三角与下三角系数各存于一数组（图 5.20）。源或右端项存于名为 source 的特定向量，维度也等于计算域内单元数。
值得指出的是，虽然 fvMatrix 类被定义为"template"，其向量化特化并不暗示矩阵本身是块耦合形式——它只是一个标量方程向量，按标准分离方式求解。这由 fvMatrix 类内部定义的代码（Listing 5.2）证实：其返回矩阵对角向量系数的成员函数无论模板类型如何，返回对象均为标量场。
uFVM 与 OpenFOAM® 中系数存储技术的更多细节将在第 7 章给出。
5.11 Closure
本章对 FVM 做了整体概述。离散化过程分为两步：第 1 步产出一组半离散化方程。本章还讨论了离散化过程的一些指导原则，以保证所得离散方程具备若干期望属性。在对第 2 步离散化过程中将用到的数值技术做详细描述之前，下章将聚焦于计算域的离散化与若干相关议题。
本章个人批注
本章是整书的"方法论开篇"——前面四章把读者安置在"守恒方程 + 网格"的工作台前，第五章则正式给出 FVM 的两步离散化骨架，并讨论离散方程应满足的若干属性。作为第一次系统读 FVM 离散化的人，我留意到 Moukalled 团队的几个强调点：
- 离散化的两步走：（1）把守恒方程在单元上积分并用 Gauss 定理转化为面积分，再用积分点（Gauss 积分）把积分转为离散求和——这一步得到的是"半离散化方程"；（2）把面/体通量用相邻单元质心变量值表达（线性化），得到代数方程。本章覆盖的是第 1 步与第 2 步的"骨架"（线性化形式），第 2 步的具体细节（差分格式、插值方式）留给后续章节。
- 二阶精度的双重来源：（1）面/体积分用中点法则（中值近似），（2）单元内 \(\phi\) 分布用线性假设。两者都贡献 \(O(\Delta x^2)\) 截断误差——任何一项精度更差都会拖累整体精度。Lilek-Peric 的 Simpson 法则四阶面积分 + 四阶 \(\phi\) 分布则是同等量级提升的两项配合。
- 守恒性是 FVM 的招牌：因为面通量基于共享该面的两个单元，因此内部面的通量在两侧天然相等相反；这就是 FVM 在传热、燃烧等守恒量问题中流行的根本原因。FVM 的 collocated arrangement 与"无需物理-计算坐标变换"是其实用工程地位的另两个支柱。
- 离散方程的属性清单：守恒性、精度、收敛、相容性、稳定性、经济性、输运性（transportiveness）、有界性。这些不只是理论口号——后续章节在设计具体差分格式时，往往要回到这份清单上来做权衡（例如：满足输运性与有界性要付出怎样的精度代价？）。
- vertex-centered vs cell-centered：Moukalled 的偏好是 cell-centered（也是 OpenFOAM 的选择），理由是实现更直接、能用任意多边形单元。但要注意 cell-centered 在非接续单元上的二阶精度丢失以及非正交网格上扩散项离散显式化导致的鲁棒性下降——这些短板在第 8 章会再次出现。
- uFVM 与 OpenFOAM®：本章 Computational Pointers 节实际上等于"本书配套代码导读"。uFVM 提供 Matlab 透明性，OpenFOAM 提供工业级 C++ 框架（ldu 稀疏存储、fvm::/fvc:: 运算符重载）。从"教学代码 → 工业代码"的视角读这两节会很有收获——很多 CFD 代码架构的基本选择（每行稀疏矩阵、面通量面 owner/neighbour）在这里都已经埋下伏笔。
- 与"模型验证"的关联：我在想，本章给出的精度分析（截断误差正比于 \((\mathbf{x} - \mathbf{x}_C)^2\)）是对均匀光滑分布而言的；当网格本身存在斜率不连续时（如图 5.16 的扭结），截断误差中出现零阶项，再加密网格也无济于事——这正是计算域离散化质量直接决定数值精度的反例。换言之，FVM 的精度不仅取决于格式，也取决于网格——这是实践中网格独立性研究（grid independence study）的理论根据。
与上下章的衔接（一段话）
第 4 章从"什么是离散化"切入，把离散化流程拆解为四大步骤（建模、网格、方程、求解），用微处理器散热器的具体例子串联。第 5 章进入方程离散化的细节，正式介绍 FVM 的两步离散化骨架与离散方程应具备的属性。读到这里，前后衔接就清晰了：第 4 章是"为什么离散化"和"离散化的宏观流程"，第 5 章是"如何做 FVM 离散化"的入门。下一章（第 6 章）将按本章预告的那样转入"计算域的离散化与若干相关议题"，即网格数据结构与几何信息——这正是 5.9 节"Mesh Support"小节预告的内容。整本书的叙事因此呈现"先理论骨架（FVM 两步走）→ 再实施细节（网格数据结构）→ 再各项细节（对流项、扩散项、源项、时间项的离散化）"的清晰层次。本章为后续章节定下了"用单元质心、面通量、稀疏矩阵"这套词汇与符号体系，是全书的术语奠基章。