第 6 章：有限体积网格（The Finite Volume Mesh）
6.1 Domain Discretization
物理域的离散化（mesh generation）在 FVM 实施中产生一个计算网格（图 6.1），后续的守恒方程就在其上求解。近几十年域离散化的方法与技术发生了巨大变化，如今大多已自动化。FVM 上下文对网格系统有若干特殊属性要求，作者以"在结构与非结构三角网格上计算变量 ϕ 的梯度"作为贯穿全章的最小工作样例来阐述这些属性。在展开 FVM 网格属性之前，作者先介绍网格系统通常应满足的特征（这些特征也用于支撑梯度计算）；这些特征既适用于结构网格也适用于非结构三角网格。
一般而言，几何域既可以用结构网格也可以用非结构网格离散。在结构网格中，三维单元由其局部下标 (i, j, k) 唯一定义；这种索引体系使得拓扑信息隐式地嵌入网格结构中，从而带来编码、缓存、向量化的效率优势，但代价是几何灵活性受限。提升结构网格灵活性的两条途径：（1）多块拼接（多 block 结构网格，各 block 独立或联合生成）；（2）放弃结构网格的隐式拓扑，改用非结构网格 + 显式 connectivity 表 + 几何实体编号的方案。结构网格曾长期是数值仿真的主流，过去 20 年非结构网格才逐渐流行。自动网格生成的研究兴趣自 1970 年代初起急剧上升，因为问题规模增大后人工网格生成过于耗时；早期方法是半自动（操作员先在域内手工布置点，再由计算机在第二步生成网格），如今已实现点和单元的完全自动生成。现代 CFD 代码大多能同时处理非结构网格与多种混合多块网格。OpenFOAM® 使用非结构网格，但也能使用 conforming 与 non-conforming 多块网格。本书将在非结构网格的语境下展开 FVM，但在刻画非结构有限体积网格的属性时，会与结构网格做对比——这种"用结构网格作参照"的做法贯穿全章，使每引入一个非结构网格上的概念都能立即看到结构网格上的对应版本。
6.2 The Finite Volume Mesh
为使后续讨论具体，作者把对 FVM 网格属性的讨论锚定在一个简单问题——单元场（element field）的梯度计算。先在结构网格上计算梯度，再到非结构网格上做同样的计算，两种路径的差异能澄清若干关键问题。
6.2.1 Mesh Support for Gradient Computation
作者采用 Green-Gauss 定理作为梯度计算的方法（这是 Chap. 9 主题之一），理由是它对各种拓扑与网格（结构/非结构、正交/非正交）都适用且形式直接。出发点是把单元 C（质心为 C，体积为 \(V_C\)）上 ϕ 的平均梯度定义为
（式 (6.1)）。然后用散度定理把体积分转化为面积分，得到
（式 (6.2)），其中 \(d\mathbf{S}\) 是面向外的表面向量元。在离散面存在的情况下，式 (6.2) 可改写为对所有面的求和
（式 (6.3)）。接下来用中点积分法则把面 f 上的面积分近似为"面质心处 ϕ 的插值 × 面积"，得到
（式 (6.4)，求和遍历 C 的所有邻接面 f）。
由式 (6.4) 与图 6.2 可以看出，要在单元 C 上计算梯度的平均值，需要两类信息：（1）面的面积与方向 \(\mathbf{S}_f\)；（2）邻接单元信息以及 ϕ 在单元质心处的取值 \(\phi_C\)、\(\phi_F\)。第二类信息用于在面质心处插值 ϕf。沿两侧节点 ϕ 值的分布假设本质上为梯度计算引入了近似；无论如何 ϕf 都必须在每个面质心处被计算——它不能"绕过"被离散面，直接从单元质心值得到。假设 ϕ 沿跨面的两个单元 C、F 之间的变化为线性，则面质心处 ϕ 的近似值 \(\hat{\phi}_f\) 可写为
\(\hat{\phi}_f = g_F \phi_F + g_C \phi_C\)
（式 (6.5)）。一种计算权因子 \(g_F\) 与 \(g_C\) 的方式是按两侧单元的体积加权
（式 (6.6)）。章末注脚提到，书中后续会介绍其他插值方案——本章 6.5.2.3 给出了三维情形下的几种替代方案。
Example 1 给出一个二维例子：对图 6.3 所示单元（中心单元 C 体积 \(V_C = 37.8\)，6 个面的 \(\mathbf{S}_f\) 在表 6.1 给出），分别在 Case 1（ϕ 恒为 1）与 Case 2（ϕ 在各面上取 10, 9, 5, 3, 4, 8）下计算梯度。Case 1 期望得到 \(\nabla \phi_C = \mathbf{0}\)：把 \(\phi_f = 1\) 代入式 (6.4) 后梯度化简为
；按表 6.1 数值计算后向量和恰为零向量。作者指出这其实是闭单元上一个普遍性质——只要所有面 \(\mathbf{S}_f\) 都朝外（或都朝内），其向量和必为零（这是 Gauss 定理在 ϕ = 常数下的特例）。Case 2 代入 \(\phi_f = (10, 9, 5, 3, 4, 8)\) 后得到
。
6.3 Structured Grids
对规则结构网格，域内每个内部单元连接到相同数目的邻接单元——这一性质使"邻接关系"成为网格结构内禀的一部分。邻接单元（图 6.4）可用 x、y、(z) 三个方向上的下标 (i, j, k) 来识别，并通过对相应下标做 ±1 直接访问。这种索引体系把拓扑信息嵌入到网格结构中，因此可以降低内存使用（无需为每单元显式存储邻居列表），同时也带来更高的编码效率、缓存利用率与向量化能力。结构网格在 FVM 与 FDM 的发展早期被广泛使用，并塑造了 CFD 编码对"按下标访问相邻数据"的习惯。
在结构网格中，可为每个计算单元关联一个有序下标集 (i, j, k)，每个下标在固定范围内变化、且独立于其他下标；相邻单元的下标差恰为 ±1——这是结构网格的"局部性"判据。若 i、j、k 方向上分别有 \(N_i\)、\(N_j\)、\(N_k\) 个单元，则域内单元总数为 \(N_i N_j N_k\)。三维情形下单元为六面体（hexahedron），6 面 8 顶点，每个内部单元有 6 个邻居；二维情形下单元为四边形（quadrilateral），4 面 4 顶点，每个内部单元有 4 个邻居。
6.3.1 Topological Information
全局下标（global index）通常用于在整个计算域上建立完整的方程组（特别是当生成的代数方程组要交给通用稀疏矩阵求解器时），而局部下标（local index）用于定义单元的局部 stencil（对离散化过程很有用）。在结构网格中局部下标与全局下标可以互换使用，因为两者之间可以相互直接转换（图 6.5），无需维护额外的查找表——这正是结构网格的便利之处。
二维情形下局部下标 (i, j) 与全局下标 n 的关系为
\(n = i + (j - 1) N_i, \qquad 1 \le i \le N_i, \; 1 \le j \le N_j\)
（式 (6.7)）。而 (i, j) 的邻居单元对应的全局下标为
（式 (6.8)）。
三维情形下关系为
（式 (6.9)），邻居 (i, j, k) 的全局下标相应地为 n、n+1、n-1、n+\(N_i\)、n-\(N_i\)、n+\(N_i N_j\)、n-\(N_i N_j\)（式 (6.10)）。这种关系大大简化了系数访问与方程组求解：按单元局部 stencil 构造的系数可以直接用于全局方程组，无需在局部与全局下标之间做显式翻译。几何场与各守恒场的访问同理。
Example 2 在 5×7 结构网格中求 (3, 4) 处单元的全局下标：\(N_i = 5, N_j = 7\) → 由式 (6.7) n = 3 + (4-1)·5 = 18；由式 (6.8) 邻居的局部/全局下标分别为 (2,4)→17, (4,4)→19, (3,3)→13, (3,5)→23。
6.3.2 Geometric Information
结构网格上访问单元周围的局部几何信息十分简单（图 6.6）。对单元 (i, j)，其周围被存储的面为 \(S_1(i, j)\)、\(S_2(i, j)\)、\(S_1(i+1, j)\) 与 \(S_2(i, j+1)\)——即只存储 i 与 j 较大方向上的两个面，下标较小方向的面则通过取负得到。由于任一单元的面都必须朝外（图 6.6），下标较小方向的面与存储的面反向：
（式 (6.11)）；而对其他面则用正向：
（式 (6.12)）。
结构网格上，单元场（二维 / 三维）也定义为大小 \([N_x] \times [N_y]\) 或 \([N_x] \times [N_y] \times [N_z]\) 的数组；因此访问单元值及其邻居同样简单。也可以把多维单元场定义为一维数组——二维时大小 \([N_i N_j]\)，三维时大小 \([N_i N_j N_k]\)。在多重网格（multi-grid）系统中，用全局下标把多维场按一维数组存储可显著节省计算机内存。
6.3.3 Accessing the Element Field
结构网格上访问单元场与使用单元下标一样简单。在二维中，\(\phi(i, j)\) 或 \(\phi_{i, j}\) 是单元 (i, j) 处的 ϕ 值。如图 6.7a 所示，(i, j) 处单元的 ϕ 在邻接单元处的值分别为 \(\phi_{i+1, j}\)、\(\phi_{i-1, j}\)、\(\phi_{i, j+1}\)、\(\phi_{i, j-1}\)。如上所述，用式 (6.4) 计算梯度需要在有限体积的每个面上计算 ϕ 值（对伪单元而言前后面具有相同值，因此不参与计算）。所以除 (i, j) 处的 ϕ 外，还需要 (i+1, j)、(i-1, j)、(i, j+1)、(i, j-1) 处的 ϕ 值。在结构网格上这些信息都现成可用，面上的 ϕ 由共享该面的两个单元的质心处 ϕ 简单插值得到。用式 (6.5)，面 (i+1/2, j) 上用局部下标表示的插值可写为
\(\phi_{i+1/2, j} = g_{i+1/2, j} \phi_{i+1, j} + (1 - g_{i+1/2, j}) \phi_{i, j}\)
（式 (6.13)）。线性插值的细节将在本章稍后给出。如图 6.7a，单元 (i, j) 上的梯度用局部下标可写为
（式 (6.14)）。若用全局下标（图 6.7b），则
（式 (6.15)）。
需要说明的是，\(\mathbf{S}\) 是控制体面上朝外的法向量。除域边界外，控制体的面被两个单元共享；因此一个单元的"朝外"恰是另一个单元的"朝内"。为避免在接口处重复存储面向量，约定只存储一个向量，其方向取为 i 或 j 增大的方向。结构网格的索引特性使"选取正确方向"无需额外存储信息：单元 (i, j) 的面下标大于 i 或 j 的取正号，小于 i 或 j 的取负号——这正是式 (6.11) 与 (6.14) 中出现负号的原因。
6.3.3.1 Discretization Indexing
除局部下标与全局下标外，还有一种称为 discretization indexing 的下标体系，场的值与几何量都用其位置或邻居值来定义——在 FVM 实际算法中这种命名最贴近"人写代码"时的命名习惯（e, w, n, s, NW, NE, SW, SE）。如图 6.7c，单元 (i, j) 上的梯度用 discretization 下标可写为
（式 (6.16)）。
相应的算法可写为：循环所有单元 (i, j)，先把单元梯度初始化为零；再循环该单元各面，根据 \(\mathbf{S}_f\) 朝外/朝内方向把 flux_f = phi_f · S_f 加或减到单元梯度上；最后把所有 flux 之和除以单元体积得到单元梯度。这种"按单元逐个循环"的写法是 6.4.1 节"按整个域循环"算法的前身——结构网格上两者都可行。
6.4 Unstructured Grids
非结构网格在"单元类型选择"与"单元加密位置选择"两方面都比结构网格更灵活——前者能使用三角形、四边形、五边形、多面体等多种单元形状，后者在 CFD 实际应用中意味着局部加密、boundary layer 处理、近壁面各向异性等手段都能直接表达。但灵活性是有代价的：在非结构网格中，单元、节点、面等几何量都按顺序编号——这意味着无法仅靠下标直接关联各种实体；必须显式定义局部 connectivity——从某个单元的几何量开始，逐步建立邻接关系。如图 6.8，单元 9 的邻居无法从其下标直接推出；其邻接面（或面上的节点）同样无法像在结构网格中那样从下标推出。因此需要为邻接单元、面、节点等提供详尽的拓扑信息（connectivity 表），以补充全局下标无法表达的关系。这是 FVM 实施在非结构网格上要解决的核心数据结构问题。
6.4.1 Topological Information (Connectivities)
如图 6.9，拓扑信息通过显式构造的局部下标（图 6.9a）和全局下标（图 6.9b）来定义几何组件的连接关系（单元-单元、单元-面、面-单元、单元-节点等）。为此，单元、面、节点的数据结构现在都包含以局部与全局下标表示的邻接信息——这些连接关系在结构网格上由下标关系隐式给出，在非结构网格上必须显式存储。梯度计算算法可基于图 6.9a 的 discretization 下标写成
（式 (6.17)），或用图 6.9a 中各面局部编号表示为
（式 (6.18)）。基于图 6.9b 中的面全局下标，梯度关系也可写为
（式 (6.19)）——注意面 23、15、11、10 项前的负号。
局部面向量 \(\mathbf{S}\) 总是被假设为朝外方向；然而实际存储的面向量并不一定朝外——因为每个面上只存储一个法向量。观察图 6.9 可知这些特定存储的面向量实际上指向单元 9 内部，所以要加负号。结构网格中"正确方向"可由下标关系直接获得（6.3 节），非结构网格中面的法向方向必须以某种方式存储。这点将在 6.4.1 的 face owner / neighbour 约定中给出具体方案。为把方向信息纳入，引入一个 sign 函数，梯度方程可改写为
（式 (6.20)）。第一项求和遍历"邻居在 k 外"的面（即 k 是 owner 之一），第二项求和遍历"邻居与 k 同侧"的面（即 k 是 neighbour）。从单元 9 的视角看，存储的 \(\mathbf{S}_{16}\)、\(\mathbf{S}_{22}\) 指向外部，所以正号；\(\mathbf{S}_{23}\)、\(\mathbf{S}_{15}\)、\(\mathbf{S}_{11}\)、\(\mathbf{S}_{10}\) 实际指向内部，所以负号。
对每个面来说，跨该面的两个单元就决定了该面的拓扑。面的方向可用"按特定顺序对两个单元下标"的约定来标准化：界面法向取为"从单元 1 指向单元 2"——OpenFOAM® 中分别称它们为 owner 与 neighbour 单元（图 6.10）。这一约定在 OpenFOAM® 与 uFVM 的实现中都被沿用。若以单元 2 的视角看这个界面，则应乘以负号。因此对单元 9 的邻接面，其 connectivity 信息定义如图 6.11 所示。
梯度本可对每个单元逐一计算（这与 6.3.3.1 在结构网格上按单元循环的算法相对应），但每个跨面单元对上的 flux \(\phi_f \mathbf{S}_f\) 对两侧单元除符号外完全相同，所以更高效的做法是按"对整个域"（如图 6.12 所示）计算梯度：循环所有面，把算出的 flux 加到 owner 单元的梯度上、减到 neighbour 单元的梯度上。因此非结构网格上计算梯度场的算法可写为：声明梯度数组并初始化为 0；循环所有内部面，对每个面计算 flux_f = phi_f · S_f，把 flux_f 加到 owner 单元的梯度上、-flux_f 加到 neighbour 单元的梯度上；循环所有边界面，把 flux_f 加到 owner 单元的梯度上；循环所有单元，把梯度除以单元体积。该算法实质上对计算域内每个单元按式 (6.20) 得到梯度。同样的算法也可用于结构网格以降低计算成本（让每个内部面只被访问一次，而不是被两侧单元各访问一次）。
Example 3 要求对图 6.13/6.14 中的非结构网格写出单元 1、5 与面 1、7、11、23 的 connectivity 数组。约定：对单元，邻居按共享面下标递增顺序存储，内部面先于边界面，二者都按下标递增；对一个面，owner 取下标较小者，neighbour 取下标较大者；边界面只有 owner，没有 neighbour；面 \(\mathbf{S}\) 的方向由 owner 指向 neighbour。结果：单元 5 的邻居为 [6, 3]、面为 [12, 13, 15, 16]；单元 1 的邻居为 [2, 6, 3, 8, 4]、面为 [1, 3, 5, 7, 14]；面 1 owner=1, neighbour=2；面 7 owner=1, neighbour=8；面 11 owner=7, neighbour=9；面 23（边界）owner=8, neighbour=-1。
6.5 Geometric Quantities
除拓扑数据外，有限体积网格还包含其几何实体的信息：单元的体积、面的面积、单元和面的质心、面与连接 owner-neighbour 单元质心向量的对齐度（图 6.15）等。本节后续将介绍这些几何量中若干的计算方法；先描述网格生成中可用的单元类型，再说明几何信息的计算技术。
6.5.1 Element Types
在 FVM 网格中，单元在三维网格中是一个多面体（图 6.16），在二维网格中是一个多边形（图 6.17）。三维网格中最常用的形状有四面体、六面体、棱柱，特殊情况下还有一般多面体。这些三维单元所对应的面类型（即二维单元类型）也多种多样，最常用的是四边形、三角形和五边形；某些应用中也使用一般多边形。值得指出的是，在二维网格中"单元的体积"被视为"二维单元面积 × 平面外方向单位长度"；因此二维网格上"单元体积"的计算方法与三维网格上"面面积"的计算方法完全相同。离散化过程中产生的其他几何量将在需要时给出。
6.5.2 Computing Surface Area and Centroid of Faces
三维 FVM 网格中单元面的一般形状是多边形，但最常用的是三角形和四边形。所有多边形类型上，\(\mathbf{S}\) 与质心的计算流程都相同：先在多边形内构造一个点，该点是定义多边形的全部点的平均——这就是多边形的几何中心 \(\mathbf{x}_G = (x_G, y_G, z_G)\)，它仅在某些特殊形状（包括三角形）下才与多边形质心 \(\mathbf{x}_{CE} = (x_{CE}, y_{CE}, z_{CE})\) 重合。因此，k 个点构成的多边形的几何中心为
以几何中心为顶点、与多边形每条边形成一个三角形（图 6.18）。对每个三角形（其几何中心与质心重合）都能直接求出质心与面积；这些三角形的面积求和即为多边形的总面积。多边形质心的求法是：对每个子三角形取"面积加权的几何中心"再按多边形面积归一化：
6.5.2.1 Surface of a Triangle
三角形面积用向量积计算：两个向量的向量积的大小代表由这两个向量张成的平行四边形的面积；故三角形的面积是这两个向量向量积大小的一半。把图 6.19 中三角形三个顶点 1, 2, 3 的位置向量分别记为 \(\mathbf{r}_1\)、\(\mathbf{r}_2\)、\(\mathbf{r}_3\)，则该三角形的面向量为
（式 (6.23)）；面积大小为
\(S = \sqrt{S_x^2 + S_y^2 + S_z^2}\)
（式 (6.24)）。
要判断面向量是否朝外，可计算其与"从单元质心 \(\mathbf{x}_{CE}\) 指向面质心 \(\mathbf{x}_{ce}\)"的位置向量的点积：点积为正则面向量朝外，为负则朝内。同样的方法也可用来判断二维情形下面向量的方向。对二维网格，单元的"面积"代表"控制体在平面外方向上具有单位深度时的体积"——即 \(V_{2D} = S_{2D} \cdot 1\)。因此二维三角单元的体积为
（式 (6.25)）。注意：若三角形顶点 1, 2, 3 按逆时针定向，则带符号的体积（或面积）为正；否则为负。取式 (6.25) 右端的绝对值总能得到正确的体积值。
Example 4 用表 6.2 的 5 顶点坐标（\((1, 6.4), (2.4, 4.0), (2, 0.2), (0.4, 0), (0, 4.0)\)）求该多边形质心与面积。k=5 时的几何中心为 \(x_G = (1 + 2.4 + 2 + 0.4 + 0)/5 = 1.16\)、\(y_G = (6.4 + 4 + 0.2 + 0 + 4)/5 = 2.92\)。多边形被分解为 5 个以 G 为顶点的三角形，表 6.3 列出 5 个子三角形的质心坐标（x-centroid 1.52, 1.85333, 1.18666, 0.52, 0.72；y-centroid 4.44, 2.37333, 1.04, 2.30666, 4.44），表 6.4 给出 5 个子三角形的面积（2.244, 2.14, 2.62, 2.104, 1.932），其和为多边形面积 \(S_t = 2.244 + 2.14 + 2.62 + 2.104 + 1.932 = 11.04\)。按面积加权求质心：
，
。几何中心 \((1.16, 2.92)\) 与质心 \((1.175, 2.826)\) 之间的差别是显然的——这是非三角形多边形上两者不重合的实例。
6.5.2.2 Volume and Centroid of Elements
一般多面体的体积与质心计算思路在概念上很简单，与多边形的"几何中心 + 子棱锥"思路完全对应——把多边形面变成多面体面、把子三角形变成子棱锥。先计算多面体单元的几何中心，把多面体分解为多个多边形棱锥。如图 6.21，每个多边形棱锥由几何中心作顶点、单元的一个多边形面作底面、侧面都是三角形。对一个多边形棱锥，体积与质心都能直接计算：体积 \(V = (1/3) \cdot \text{底面积} \cdot \text{高}\)；底面即单元的一个面；棱锥的质心（自底面质心起算）位于"连接底面质心与棱锥顶点的连线"的 1/4 处——具体推导是把棱锥的体积分拆为多个平行底面切片积分的极限，每个切片的质心在该切片中点处，由此得出 1/4 系数。多面体单元的体积等于这些棱锥体积之和；质心则按各棱锥质心做体积加权平均。数学上为
（与式 (6.21) 同一形式，只是对单元顶点），
（棱锥质心 = 0.75 × 底面质心 + 0.25 × 棱锥顶点，即"1/4 距底面 / 3/4 距底面"关系），
（式 (6.26)）。注意 \(V_{\text{pyramid}} = d_{Gf} \cdot \mathbf{S}_f / 3\) 的几何含义：\(d_{Gf} \cdot \mathbf{S}_f\) 给出"以 \(\mathbf{S}_f\) 为底、\(d_{Gf}\) 为法向高"的平行六面体体积（对朝外法向取正值），除以 3 即得棱锥体积。
由于多边形棱锥的质心与体积容易计算，上述方法能精确计算一般多面体单元的体积与质心；这也是 6.6 节 uFVM 与 OpenFOAM® 实际代码中"按面遍历、累加棱锥体积与体积加权质心"算法的理论根据。
6.5.2.3 Face Weighting Factor
考虑图 6.22 所示的一维 FVM 网格系统。ϕ 在控制体质心 C 与 F 上的值已知，要用来计算界面 f 上的 ϕf。简单线性插值公式为
\(\hat{\phi}_f = g_f \phi_F + (1 - g_f) \phi_C\)
（式 (6.27)），其中
（式 (6.28)）。这个简单公式在多维情形下并不直接可用——在二、三维时，几何权因子的定义不止一种选择。一种选择是按两侧体积定义权因子
\(g_f = \frac{V_C}{V_C + V_F}\)
（式 (6.29)）。但这种方法在某些情形下会给出错误结果，例如图 6.23 所示轴对称网格系统。另一个困难是当 C、f、F 三点不共线时（图 6.24a）。对这种情形更好的方案如图 6.24b，是按"到面的法向距离" \(d_{Cf'}\)、\(d_{fF'}\) 做插值
（式 (6.30)），其中 \(\mathbf{e}_f\) 为面的单位向量
（式 (6.31)）。
Example 5 对图 6.25 与表 6.5 的两个三角单元比较式 (6.29)（按体积加权）与式 (6.30)（按法向距离加权）给出的 \(g_f\)。两单元顶点坐标分别为 C 单元 \((0, 0), (1.2, 0.4), (1, 1)\) 与 F 单元 \((1.2, 0.4), (1, 1), (2, 0.1)\)。按式 (6.29) 得 \(V_C = 0.4\)、\(V_F = 0.42\) → \(g_f = 0.4/0.82 = 0.4878\)。按式 (6.30)：\(C = (0.7333, 0.4666)\)、\(F = (1.4, 0.5)\)、面质心 \(f = (1.1, 0.7)\)、\(\mathbf{S}_f = 0.6 \mathbf{i} + 0.2 \mathbf{j}\)（按 \(\mathbf{S}_f = (y_3 - y_2) \mathbf{i} - (x_3 - x_2) \mathbf{j}\) 算得，\((y_3 - y_2) = 0.6, (x_3 - x_2) = -0.2\) 故取负号后 \(\mathbf{S}_f = 0.6 \mathbf{i} + 0.2 \mathbf{j}\)）、\(S_f = \sqrt{0.6^2 + 0.2^2} = \sqrt{0.4} = 0.6325\)、
；
、
；
、
→ \(g_f = 0.4217 / (0.4217 + 0.2215) = 0.4217 / 0.6432 = 0.6556\) 。两种方法的结果差异是明显的（0.4878 vs 0.6556）——这正说明在多维非正交情形下选错插值方案会得到相当不同的 \(g_f\)，进而影响梯度计算精度。
6.6 Computational Pointers
6.6.1 uFVM
uFVM 中，所有几何与拓扑数据的处理都集中在一个名为 cfdProcessOpenFoamMesh 的例程中；该例程在用 cfdReadOpenFoamMesh 读入 OpenFOAM® 网格后立即执行。读入原生 OpenFOAM® 网格需要按以下顺序读多个文件：points、faces、owners、neighbours、boundaries（这种"先点后面再 owner/neighbour"的顺序是 OpenFOAM® polyMesh 格式的标准）。Listing 6.1 给出了 cfdProcessOpenFoamMesh 中对"基础面几何"（centroid、area、\(\mathbf{S}_f\)）的处理。
Listing 6.1 的算法分三层循环：外层 for iFace = 1 : numberOfFaces 遍历所有面；中层 for iNode ∈ iNodes 累加各节点质心得到该面的"粗略中心" centre = centre / numberOfiNodes（这是 6.5.2 公式 (6.21) 对面顶点的直接应用）；内层 for iTriangle = 1 : numberOfiNodes 把面拆为 numberOfiNodes 个虚拟三角形（每条边加 centre 构成一个），对每个子三角形算 centroid = (point1+point2+point3)/3（式 (6.21) k=3）与 \(\mathbf{S}\) = 0.5 × cross(point2-point1, point3-point1)（式 (6.23)）与面积 cfdMagnitude(local_Sf)。三层循环结束后做归一化：centroid = centroid / area（即除以总面积 S_f，按式 (6.22) 归一），把 centroid、Sf、area 写回 theMesh.faces(iFace)。
Listing 6.2 给出"基础单元几何"的处理。其核心是按"棱锥法"算单元体积与质心：先用单元各面的 centroid 平均出"粗略中心" centre = centre / length(iFaces)（对应 6.5.2.2 公式 (6.26) 第一式）；再遍历单元各面，按 localFaceSign = theMesh.elements(iElement).faceSign(iFace) 给 Sf 加正负号（即面的方向，与 6.4.1 的 sign 函数对应），用 \(V = \mathbf{S}_f \cdot \mathbf{C}_f / 3\) 算该面对应棱锥的体积（\(\mathbf{C}_f\) 是从 centre 到面 centroid 的向量，等价于 6.5.2.2 中的 \(d_{Gf}\)，但符号取 faceSign 后自动对齐外法向）；用
\(0.75 \cdot \text{localFace.centroid} + 0.25 \cdot \text{centre}\)
算该棱锥的质心（对应 6.5.2.2 公式 (6.26) 第二式）；把棱锥的"体积 × 质心"累加到 localVolumeCentroidSum，棱锥体积累加到 localVolumeSum；最后用 centroid = localVolumeCentroidSum / localVolumeSum、volume = localVolumeSum 得到单元质心与体积。
6.6.2 OpenFOAM®
OpenFOAM® 采用非结构网格平台，因此其所有几何实体——单元、体积、面积、质心以及面权因子——都需要被显式求值并存储。本节概述在 OpenFOAM® 中评估这些几何量所涉及的几何关系。
6.6.2.1 Area and Centroid of Faces
为评估一般多边形面的面积与中心，把面分解为一系列三角面，把各三角形部分的性质累加得到整个多边形面的度量——这与 6.5.2 公式 (6.22) 描述的策略相同，只是 OpenFOAM® 用代码而非数学形式实现。因此对每个计算单元都要应用式 (6.21)–(6.23)。OpenFOAM® 在文件 $FOAM_SRC/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C 中构造面中心并计算其面积，对应的函数如 Listing 6.3 所示 makeFaceCentresAndAreas(const pointField& p, vectorField& fCtrs, vectorField& fAreas)，它有三个参数：第一个 const pointField& p 为 const 参数代表从文件读入的数据；第二、三个参数 vectorField& fCtrs、vectorField& fAreas 分别代表返回的、与域内面数等长的"面中心列表"与"面面积列表"。
pointField 数据是所有网格顶点的列表（OpenFOAM® 中顶点定义为 point 类型，本质是三维向量），每个顶点用三个空间坐标定义。所需的第二份数据是面定义——在 OpenFOAM® 中由 faceList（Listing 6.4）给出，它是一个 face 对象的列表，每个 face 内部存的是"定义该面的顶点 id"列表。Listing 6.4 中用 forAll(fs, facei) 遍历所有面，nPoints = f.size() 是该面的顶点数。
然后对每个面执行循环（Listing 6.4）：在循环内先读出描述该面的点数与各点的 id 以便访问对应顶点的坐标。当面只由 3 个顶点定义时，直接计算其质心位置与面积（Listing 6.5）以提高效率并避免舍入误差相关问题：此时面中心 fCtrs[facei] = (1.0/3.0) * (p[f[0]] + p[f[1]] + p[f[2]]) 用式 (6.21) 取 k=3 直接求；面积 fAreas[facei] = 0.5 * ((p[f[1]] - p[f[0]]) ^ (p[f[2]] - p[f[0]])) 用式 (6.23) 计算（其中符号 ^ 代表向量积，对应 6.5.2.1 中的 \(\times\)）。
对一般多边形，先把面分解为三角形；为此 OpenFOAM® 借"所有顶点平均"先估计一个面中心 fCentre（Listing 6.6）——这一步对应 6.5.2 公式 (6.21)，对当前面所有顶点取平均得到几何中心 \(\mathbf{x}_G\)。
Listing 6.7 表明对所有面做循环，把面分解为子三角形并计算其几何中心与面积：内层 for (label pi = 0; pi < nPoints; pi++) 遍历当前面的所有顶点；每步用 nextPoint = p[f[(pi + 1) % nPoints]] 取下一个顶点（用 % nPoints 实现首尾相连）；子三角形的质心 c = p[f[pi]] + nextPoint + fCentre（系数 1/3 在后面乘，即对应 6.5.2 公式 (6.21) 形式但暂不归一）；面法向量 n = (nextPoint - p[f[pi]]) ^ (fCentre - p[f[pi]])、面积大小 a = mag(n)（系数 1/2 在后面乘，对应 6.5.2 公式 (6.23) 但暂不归一）；把 n、a、a*c 累加到 sumN、sumA、sumAc。
如 Listing 6.8 所示，最后做归一化：在退化面情形（即 sumA < ROOTVSMALL）下，为面设置一个 rescue 值（fCtrs[facei] = fCentre; fAreas[facei] = vector::zero;）；否则按式 (6.22) 归一——fCtrs[facei] = (1.0/3.0) * sumAc / sumA（即 1/3 系数补上 6.5.2 公式 (6.21) 的归一化）、fAreas[facei] = 0.5 * sumN（即 1/2 系数补上 6.5.2 公式 (6.23) 的归一化）。这一"先累加、最后归一"的写法与 6.5.2 公式 (6.22) 形式完全一致，是把 1/3 与 1/2 系数"延迟到最后"应用的实现技巧。
6.6.2.2 Volume and Centroid of Elements
在完成面法向量、面积、中心的计算后，就可以计算单元的度量。与面的方法类似，多面体单元的基本思路是把它分解为多个四面体之和（具体在 OpenFOAM® 中是"以面 centroid + 单元几何中心为顶点的棱锥"——见 6.5.2.2 公式 (6.26) 的几何推导）。计算单元体积与质心的过程定义在文件 `
x_G$）——cEst[celli] /= nCellFaces[celli]。
得到 \(x_G\) 后，按 6.5.2.2 公式 (6.26) 与 Listing 6.11 计算各棱锥的体积与单元质心：对每个面算 pyr3Vol = max(fAreas[facei] & (fCtrs[facei] - cEst[own[facei]]), VSMALL)，其物理含义为 \(3 \cdot V_{\text{pyramid}}\)（与 6.5.2.2 公式 (6.26) 第三式 \(V_{\text{pyramid}} = d_{Gf} \cdot \mathbf{S}_f / 3\) 一致，
\(fAreas[facei] \cdot (fCtrs[facei] - cEst[own[facei]])\)
即 \(\mathbf{S}_f \cdot d_{Gf}\)，\(VSMALL\) 是一极小正数用于"防止 0 体积"）；算面棱锥的质心 pc = (3.0/4.0) * fCtrs[facei] + (1.0/4.0) * cEst[own[facei]]（与 6.5.2.2 公式 (6.26) 第二式
一致）；把 pyr3Vol * pc 累加到 cellCtrs[own[facei]]、把 pyr3Vol 累加到 cellVols[own[facei]]；对 neighbour 端做对称处理——唯一差别是 dot product 中 \(fCtrs[facei] - cEst\) 顺序反转（fAreas[facei] & (cEst[nei[facei]] - fCtrs[facei])），这是因为对 neighbour 单元而言 \(\mathbf{d}_{Gf}\) 从 cell centre 指向面 centroid 的反向。
如 Listing 6.12 所示，最后用 cellCtrs /= cellVols 与 cellVols *= (1.0/3.0) 得到单元质心与体积的最终值——也就是说，Listing 6.11 中实际累加的是 \(3 \cdot V_{\text{pyramid}}\)，因此最后要再除 3；同时 cellCtrs 此时累加的是
，除以 \(\sum 3 V_{\text{pyramid}}\) 才得到体积加权质心
。fCtrs 对应 \(\mathbf{x}_{CE}\)；"cEst - fCtrs" 对应图 6.21 中的距离向量 \(\mathbf{d}_{Gf}\)——这一对应关系在 Listing 6.11 中被直接使用。
uFVM 与 OpenFOAM® 中网格数据结构的细节将在第 7 章给出。
6.7 Closure
本章给出了定义有限体积网格的几何数据。本章强调，FVM 网格不仅仅是"非重叠单元与节点的集合"；它还包括所有几何量及其拓扑信息。所有这些信息的集合构成了本书记述的方程离散化方法（即 FVM）所需的基础设施。
6.8 Exercises
本章个人批注
本章是 Moukalled 团队对"FVM 网格"的系统开篇——第 5 章给出 FVM 离散化的两步走骨架（守恒方程在单元上积分 + 通量在面上求和），本章则把隐式假设的"网格支撑"显式化。从"5.9 The Mesh Support"那段一两句的"未来章节主题"到本章六大节，体量翻了几倍，但叙事依然是"先讲为什么需要 mesh support、再讲结构网格上的简单情形、再讲非结构网格上的复杂情形、再讲几何量算法、最后给出 uFVM 与 OpenFOAM® 的实现参照"。我读下来特别有印象的几点：
- gradient 作为最小工作样例：作者不直接展开守恒方程离散，而是用"ϕ 的梯度"作为最小例子来揭示 FVM 网格需要哪些信息——面质心、面向量、邻接单元列表、跨面插值权因子。这是一种"由最小可工作反例倒推必要支撑"的写法，对工程实现很有借鉴价值。
- 结构网格 vs 非结构网格的 contrast：结构网格靠 (i, j, k) 隐式表达拓扑，邻居/面都靠下标算出来；非结构网格必须显式存 connectivity 表（element→neighbors/faces/nodes, face→owner/neighbor）。前者省内存、易向量化但几何灵活性差；后者灵活但存储/计算都更贵——这正是结构网格"长期主流 + 近 20 年被非结构网格反超"的根本原因。
- 面的方向性（face sign）：在 6.3.2 作者用下标差来"白送"出符号；6.4.1 则必须显式引入 sign 函数（Eq. 6.20）。这是非结构网格实现最容易踩坑的地方——存储面向量只存一份（owner→neighbour 方向），使用时若以 neighbour 视角就要加负号。OpenFOAM 的
faceSign字段本质就是把这个 sign 表存起来。 - 几何中心 vs 质心：6.5.2 明确指出二者仅在特殊形状（含三角形）下重合。这是个细节但很关键——"用顶点平均的几何中心"作为棱锥顶点时，可保证多边形能严格分解为非重叠子三角形；如果用质心则一般做不到。
- face weighting factor：6.5.2.3 给的三个备选方案（按距离、按体积、按法向距离）有微妙差异，Example 5 算出 0.4878 vs 0.6556 的差距——这是个很值得记住的对比，提醒在多维非正交情形下不要想当然用"按体积加权"。
- uFVM / OpenFOAM® 实现对照：两个代码都用"棱锥分解"算体积与质心（Eq. 6.26），但 uFVM 在 Listing 6.2 用的是"按面 centroid 累加的 volume-weighted face-pyramid centre"，而 OpenFOAM 在 Listing 6.11 用的是"用 owner / neighbour 分别累加、最后归一"。两段代码思想相同、形式不同——这是一个很好的"同一算法的两种实现风格"对照实例。
- 个人疑问：（1）为什么 Example 1 的 Case 1 中"6 个 \(\mathbf{S}_f\) 之和恰为 0"？这其实是闭单元上的一个恒等式（Gauss 定理在 ϕ=常数下的特例），但作者只用一句"this is actually a property of the surfaces of closed elements"带过——值得在 Chap. 9 重新审视。（2）6.5.2.2 的 Eq. 6.26 把"棱锥的质心"写成 0.75 (x_CE)_f + 0.25 (x_G)_pyramid，这是基于"质心位于 1/4 处（自底面起算）"的结论——这个 1/4 系数实际只对"底面为平面的均匀棱锥"成立；如果底面是非平面多边形这个推导要重看。（3）6.4.1 Eq. 6.20 引入的 sign 函数在 OpenFOAM 实际代码里是
faceSign(iFace)字段，但 Listing 6.2 那里也是localFaceSign——两个名字指同一物，没在文中点明。 - 与上下章的衔接：第 5 章"5.9 The Mesh Support"小节明确预告"网格描述将是下两章的主题"；第 6 章就是下两章的第一章。下一章（第 7 章）"The Finite Volume Mesh in OpenFOAM and uFVM"应当展开 OpenFOAM 的 LDU addressing 细节与 uFVM / OpenFOAM 的具体 mesh 数据结构——这正是 6.6 节预告的内容（"The mesh data structure for uFVM and OpenFOAM® will be described in detail in Chap. 7"）。
与上下章的衔接（一段话）
第 5 章末尾用"5.9 The Mesh Support"预告 FVM 网格将是下两章的主题，本章即按此承诺展开——先讨论域离散化的概念（结构 vs 非结构、为什么需要非结构、自动化历程），再以"梯度计算"为最小例子引出 mesh support 所需的具体信息（面面积、面向量、邻接单元、跨面插值权因子）。结构网格因下标体系隐式编码了拓扑，在 6.3 节被作为"易处理情形"先行讨论，给出下标转换公式 (6.7)–(6.10) 与 Example 2；非结构网格（6.4 节）则显式引入了 connectivity 表（element→neighbors/faces/nodes, face→owner/neighbour）与 face sign 概念，并通过 Eq. (6.20) 给出带符号的梯度公式，使结构网格上的"白送"符号在非结构网格上变为显式存储。6.5 节进一步把几何量（面面积、面质心、单元体积、单元质心、面权因子）按"棱锥分解 + 面积加权质心"统一处理，并讨论了一维 / 二维 / 三维的不同复杂度。6.6 节给 uFVM 与 OpenFOAM® 的实现指路（cfdProcessOpenFoamMesh、makeFaceCentresAndAreas、makeCellCentresAndVols），明确预告 6.6.2.2 中"LDU addressing"将在第 7 章介绍。下一章（第 7 章）将按 6.6 节预告的 LDU addressing 主题展开 mesh data structure 的细节，把 owner / neighbour / upper / lower 数组具体化，回应 6.4.1 与 6.6 节埋下的伏笔——因此第 6 章本质上是为第 7 章"具体 mesh data structure"做几何与拓扑的概念准备。Exercises 部分（6.8）不涉及具体细节，跳过不涵盖。