跳转至

第 7 章：OpenFOAM® 与 uFVM 中的有限体积网格（The Finite Volume Mesh in OpenFOAM® and uFVM）

7.1 uFVM

uFVM 是一个非结构化有限体积教学代码，作者写它的目的有两个：一是演示构成一个 CFD 程序的各种数值技术与算法；二是它的数值方法在许多方面与 OpenFOAM® [1] 相似，因此可以作为载体来理解和展示 OpenFOAM® 的内部实现。uFVM 使用的主要数据结构总体上镜像了 OpenFOAM® 的设计，特别是在网格场（mesh fields）和边界条件方面；同时，作者刻意保留 uFVM 与 OpenFOAM® 之间的若干差异，用以突出 CFD 编码者的多种可选实现方案，从而更好地呈现某些实现细节。本节将围绕 cavity 测试算例中的 polyMesh 文件夹展开。

7.1.1 一个 OpenFOAM® 测试算例

uFVM 能够读取任何一个 OpenFOAM® 测试算例中所附带的 OpenFOAM® 网格。一个 OpenFOAM® 测试算例本质上是一个目录，其内部一般至少包含三个子文件夹。图 7.1 展示了 cavity 算例的内部结构。第一个是“时间”目录族，每个时间步对应一个子目录，用于存放该时刻场变量的初始化值、边界条件信息以及该时间步的计算结果；子目录的名称即为对应的时间，例如 cavity 算例中速度场 U 与压力场 p 的初始化分别来自 0/U 与 0/p 文件。第二个是 system 目录，其中至少包含三个文件：controlDict 用于设置仿真的起止时间、时间步长与数据输出参数等总体控制参数；fvSchemes 用于定义离散格式；fvSolution 用于设定求解算法与松弛因子。第三个是 constant 目录，其中包含描述流体物理属性的文件（如 transportProperties）和描述网格系统的 polyMesh 子目录。为了理解 uFVM 中的有限体积网格结构，后续讨论将聚焦在 polyMesh 文件夹，因为构造有限体积网格所需的全部信息都定义在它之中。

7.1.2 polyMesh 文件夹

polyMesh 子目录包含 points、faces、owners、neighbours、boundary 五个核心文件。

points 文件是一个向量列表，记录每个网格顶点的坐标，列表中的第一个向量即为顶点 0，第二个向量为顶点 1，依此类推，其格式如 Listing 7.1 所示（先给顶点数，再给出每个顶点的 (x y z) 三元组）。一个实际的 points 文件示例如 Listing 7.2，头部声明顶点总数为 1074，随后以每行 (x y z) 的方式列出顶点，例如 (32 16 0.9377383239)、(33.9429245 16.11834526 0.9377383239) 等。

faces 文件是面的列表，每一面用其在 points 列表中的顶点下标序列来描述，列表第一项对应面 0，第二项对应面 1，等等，其格式如 Listing 7.3 所示（先给面总数，再逐面给出构成该面的顶点数与顶点下标）。Listing 7.4 给出一个示例：头部声明面总数为 3290，然后逐行列出每面，例如 4(36 573 589 52)、4(41 578 634 97)，其中行首数字表示构成该面的顶点数。

owners 文件存放每一面的 owner 单元下标，其在列表中的位置即该面的索引：列表中的第一个值就是面 0 的 owner，第二个值是面 1 的 owner，依此类推；owner 总数等于面总数（即内部面加边界面的总和）。owners 的元素数量等于最大 owner 下标值，即单元数。owners 文件的格式如 Listing 7.5 与示例 Listing 7.6 所示。域中的单元总数 nCells 可在 owners 文件的表头注释中找到，Listing 7.7 显示了一个典型表头 nPoints:1074 nCells:918 nFaces:3290 nInternalFaces:1300。

neighbours 文件是邻居单元下标列表，其元素个数基本等于内部面的个数，格式如 Listing 7.8 所示，示例如 Listing 7.9，头部声明邻居总数为 1300（例如 22、68、29、96、31、34 等）。

boundary 文件列出域的所有边界，每一类边界上的面集合称为一个 patch，并赋以名字；每个 patch 在文件中声明其 type、nFaces（面上的面数）和 startFace（该 patch 的第一个面在面列表中的起始下标），其通用格式如 Listing 7.10 所示。一个 wall 类型 patch 的实例如 Listing 7.11：type wall、nFaces 100、startFace 1300。

7.1.3 uFVM 网格

在 uFVM 中，读取一个 OpenFOAM® 网格的脚本是 cfdReadOpenFoamMesh。其读取顺序为：先读取 points 文件，把 (x, y, z) 坐标存入 struct nodes 数组；再读取 faces 文件，把面所对应的顶点下标存入 struct faces 数组；然后从 boundary 文件读取面 patch 的信息；最后读取 owners 和 neighbours 文件，并组合生成 struct elements。读取得到的数据随后在脚本 cfdProcessOpenFoamMesh 中被进一步处理，计算出额外的几何与拓扑信息。

以 elbow 网格为例，Listing 7.12 给出了 uFVM 在读取后能够显示的信息摘要：numberOfNodes = 1074，numberOfFaces = 3290，numberOfElements = 918，numberOfInteriorFaces = 1300，numberOfBoundaries = numberOfPatches = 6，numberOfBElements = numberOfBFaces = 1990。该网格可以用 cfdPlotMesh 命令可视化（图 7.2）。

struct nodes 中存放的信息以节点 1 为例如 Listing 7.13 所示：centroid 是该节点的坐标向量，index 是节点编号，iFaces 是与该节点相连的面下标列表（如 [172 328 1355 1386 1677 1891 1893]），iElements 是与该节点相连的单元下标列表。

struct faces 中存放的信息以面 3 为例如 Listing 7.14 所示：iNodes 是定义该面的顶点下标；index 为该面的编号；iOwner 与 iNeighbour 分别是该面的 owner 与 neighbor 单元下标；centroid 为该面的形心；Sf 为该面的面向量；area 为该面的面积；T 是连接 owner 与 neighbor 单元形心的距离向量；gf 为几何因子（geometric factor）；CN 是从 owner 单元形心指向面形心的距离向量；walldist 是 owner 单元形心到最近壁面的法向距离；iOwnerNeighbourCoef 与 iNeighbourOwnerCoef 是后续插值所需的系数；对于边界面，其 neighbor 索引被设为 -1。

struct elements 中存放的信息以单元 20 为例如 Listing 7.15 所示：包含 iNeighbours（邻居单元下标列表）、iFaces（相邻面下标列表）、iNodes（节点下标列表）三个索引列表；以及 volume（单元体积）、faceSign（一个长度为该单元面数的列表，正负号标记该单元在对应面上是 owner (+1) 还是 neighbor (-1)）、numberOfNeighbours（邻居个数）、centroid（单元形心）。Listing 7.16 给出在网格上高亮显示所选单元（如单元 20 与 300）的 cfdPlotElements 命令。Listing 7.17 展示另一个示例单元 300 的属性。元素与面的下标顺序保持同步：单元按其在网格中的顺序列出，面按其在网格中的顺序列出，且边界面定义在面列表的尾部。元素与面的下标顺序对应关系是：单元索引与该单元对应的面在面列表中的顺序一致。

struct boundaries 存储各个边界 patch 的信息。每个 boundary 数组包含该 patch 第一个边界面的起始下标、属于该 patch 的边界面数量、物理类型与名字。以 patch 1 为例（Listing 7.18），其 userName 为 wall-4，type 为 wall，numberOfBFaces = 100，startFace = 1301。该 patch 的第一个边界面信息可通过 Listing 7.19 获得：注意该边界面的下标是 1300+1，因为 MATLAB® 数组从 1 开始，而 C 语言数组从 0 开始。要遍历某个 patch 的所有面，需要用到其 startFace 与 numberOfBFaces，Listing 7.20 给出遍历 patch 2 的面的示例脚本，其中关键辅助函数 cfdGetFaceIndicesForBoundaryIndex 的实现见 Listing 7.21：它读取 boundary 的 numberOfBFaces 与 startFace，然后返回下标范围 [theStartFace : theStartFace+theNumberOfBFaces-1]。

7.1.4 单元场（The Element Fields）

除 struct mesh 中存储的所有数据外，求解模型所需的额外信息以及感兴趣物理场的值也需要存储，以便随时访问。三类 locale（节点、单元、面）以及三种数据类型（标量、向量、张量）的组合对应不同种类的场；uFVM 中的场首先按 locale 分类定义。

构造一个单元场的脚本形式为 Listing 7.22：cfdSetupMeshField(theUserName, theLocale, theType, theTimeStep)，其中四个参数分别是：场的名字（theUserName）、几何位置 locale（Elements / Faces / Nodes）、数组元素类型（Scalar / Vector）以及时间步标识（如 Step0、Step1 等）。Listing 7.23 给出一个示例：cfdSetupMeshField('U:water', 'Elements', 'Vector', 'Step0')，它建立一个定义于单元、当前时间步的向量场。

如图 7.4 所示，该数组的长度等于 numberOfElements + numberOfBoundaryFaces：前半段存放内部单元的值，后半段按 patch 顺序依次存放各边界面的值，这些边界值即该场的边界条件。

Listing 7.24 给出初始化 patch 1 上 UField 边界值的脚本：先通过 theMesh 拿到第 iPatch 个 boundary 的 startFace 与 numberOfBFaces，再计算对应的边界单元下标区间 iElementStart = numberOfElements + iFaceStart - numberOfInteriorFaces 到 iElementStart + numberOfBFaces - 1，最后通过 UField.phi(iBElements, :) = cfdComputeFormulaAtLocale('[1;0;0]', 'BPatch1', 'Vector') 将其设为 [1;0;0]。其中 cfdComputeFormulaAtLocale（Listing 7.25）在指定 locale 上对给定表达式求值，返回长度适当、类型匹配的数组。

7.1.5 面场（The Face Fields）

面场的构造脚本与单元场类似，只是 theLocale 设为 'Faces'（Listing 7.26）。如图 7.5 所示，面场数组的长度等于 numberOfFaces，即 numberOfInteriorFaces 加上所有边界面的总数；其内部面的下标从 1 到 numberOfInteriorFaces，边界面随后按 patch 顺序排列。访问某一 patch 的边界面的脚本如 Listing 7.27，其中关键的下标换算 iElementStart = numberOfElements + iFaceStart - numberOfInteriorFaces 与单元场一致；在此基础上，标量场的边界值可通过 phi_b = phi(iBElements)（Listing 7.28）取得。

7.1.6 节点场（The Node Field）

节点场是一个由面所引用的顶点下标列表来确定的列表。每个面由其顶点在 points 列表中的下标指代（图 7.6），面在列表中的位置即该面的下标：列表中的第一项是面 0，第二项是面 1，依此类推。

7.1.7 在 uFVM 网格上工作（遍历单元）

在离散化与求解循环中，遍历单元、内部面、边界面、边界单元乃至边界 patch 是常见的操作。本子节先讨论遍历单元：由于单元数已知且单元下标从 1 到 numberOfElements（OpenFOAM® 中为 0 到 numberOfElements-1），且单元场的下标方式相同，遍历脚本非常简单，如 Listing 7.29 所示。访问边界单元的脚本如 Listing 7.30 所示。

7.1.8 在 uFVM 网格上工作（遍历面）

面数组的构造使得所有内部面位于 1 到 numberOfInteriorFaces 的下标范围，边界面紧随其后并按所属 patch 顺序排列。因此遍历内部面的脚本如 Listing 7.31 所示；遍历所有边界面的脚本如 Listing 7.32 所示；若只关心某一 patch 的边界面，则使用该 patch 在 boundary 中定义的 startFace 与 nFaces，循环如 Listing 7.33 所示。子例程 cfdGetFaceIndicesForBoundaryIndex 可直接返回向量 startFace : startFace+nFaces-1。

7.1.9 计算 Gauss 梯度

在 uFVM 中计算单元场的 Gauss 梯度需要综合使用前述若干例程。函数 cfdComputeGradientGauss0（Listing 7.34）展示了 Gauss 梯度的实现细节。

其流程是：

  1. 初始化输出数组 phiGrad = zeros(numberOfElements+numberOfBFaces, 3)。
  2. 内部面贡献：取出内部面下标 iFaces = 1:numberOfInteriorFaces 与边界面的下标 iBFaces = numberOfInteriorFaces+1:numberOfFaces；取出所有内部面的 owner 与 neighbor 索引、表面向量 Sf、几何因子 gf；用线性插值计算内部面处的 φ 值 phi_f = gf.*phi(iNeighbours) + (1-gf).*phi(iOwners)；再遍历每张内部面，将 phi_f * Sf 同时累加到 owner 与 neighbor 单元上（注意相邻单元符号相反）。
  3. 边界面贡献：取出所有边界面的 owner iOwners_b、边界 φ 值 phi_b = phi(iBElements)、边界表面向量 Sb；遍历每张边界面，将 phi_b * Sb 累加到其 owner 单元上。
  4. 求平均：将每个单元上的累加值除以该单元的体积 volumes，得到该单元形心处的平均梯度。
  5. 边界梯度赋值：将每个边界单元（即每个边界面在 phi 数组中对应的槽位）的梯度直接设为对应内部 owner 单元的梯度，即 phiGrad(iBElements, :) = phiGrad(iOwners_b, :)。

Listing 7.35 给出另一种边界贡献的实现思路：不是简单遍历所有边界面，而是先遍历 patch，再在每个 patch 内部遍历其各自的边界面。两种写法在数学上等价，但 Listing 7.35 的版本与 patch 概念的组织方式更贴合，便于扩展。

7.2 OpenFOAM®

OpenFOAM® [1] 使用有限体积、单元中心离散格式，并基于所谓 face-addressing 存储方式处理非结构网格数据。这种数据结构的目的是为非结构网格的定义提供最大灵活性，从而允许使用任意多面体形状的网格单元。如图 7.7 所示，多面体是由若干平面多边形在棱处相互连接而成的三维实体；四面体（四个三角形面）和六面体（六个四边形面）就是常见的两种多面体。能够描述任意三维形状并将其作为有限体积单元用于方程离散，为网格生成带来了多方面的优势与灵活性。

描述多面体网格的有效方式是 face addressing（面寻址）。在 face addressing 中，单元的形状对离散过程没有影响，形成系数时需要使用全局寻址（关于从局部与全局两个角度形成系数的细节将在后续讨论）。

在 OpenFOAM® 中，points、faces、elements 的数据分别存储在若干列表（数组）中（图 7.8）。points 列表存放三维空间坐标向量，对应实际网格的顶点，单位为米；每个顶点有一个标签，由其在列表中的位置决定，由于使用 C++ 实现，标签计数从零开始。faces 列表存放顶点标签序列，相邻的两个顶点用一条边相连；面列表的组织方式是所有内部面排在前，随后依次是第一个边界对应的面、第二个边界对应的面，等等；边界上的面集合也称为 patch。需要记住的是：内部面属于两个单元，而边界面只属于一个单元。

elements（或 cells）列表由索引构成：列表中的位置即单元的索引，每个位置上第一个索引是该单元的面数，其后的索引是该单元的各个面。OpenFOAM® 可以读取任何能够写出所需强制文件的网格生成软件所生成的网格。前面已经说明，网格名为 polyMesh，必须用一组规定的文件来定义，这些文件放置在 constant/polyMesh 目录下，包括 points、faces、owner、neighbour、boundary 五个文件。

OpenFOAM® 中域的整个边界由 boundary 文件描述。boundary 文件是所有已定义 patch 的列表，每个 patch 用一个字典条目声明，其语法如 Listing 7.36：先给出 patch 名字，再给出 type、nFaces 与 startFace。Listing 7.37 展示一个含两个 patch（inlet 与 outlet）的 boundary 文件示例。

从编程角度，有必要简要介绍处理网格并允许访问特定数据的 C++ 类。处理网格“底层”结构的基础类叫 primitiveMesh，它是一个包装几何信息的通用类，不假设任何特定的离散形式，是关于网格的低层信息的基础类。Listing 7.38 列出了 primitiveMesh 类的若干成员函数，如 cellCells()、pointCells()、cells()、cellCentres()、faceCentres()、cellVolumes()、faceAreas()。该类本身并不识别边界或域接口；这些信息在从 primitiveMesh 派生的 polyMesh 类中定义。除了拥有 primitiveMesh 的所有属性外，polyMesh 还引入了边界定义与边界信息的处理（Listing 7.39）：它提供 owner()、neighbour() 等函数获取内部面的 owner/neighbour，还提供 V0() 返回上一时间步的单元体积，V00() 返回上上时间步的单元体积，delta() 返回面间的 delta 向量作为 surfaceVectorField。

在 primitiveMesh 与 polyMesh 之上，fvMesh 类继承自 polyMesh，增加了有限体积离散所需的数据与函数，可访问寻址信息、边界信息以及特定的网格数据。OpenFOAM® 将边界网格分解为 patches 存储在 polyBoundaryMesh 类下的 polyPatchList 中；与内部网格类似，从 polyBoundaryMesh 派生出 fvBoundaryMesh，继承其功能并扩展以包含有限体积离散所需的数据与函数。对于内部离散，存在 primitivePatch、polyPatch、fvPatch 的类层次结构。fvPatch 用于在有限体积离散中实现边界条件。图 7.9 给出 OpenFOAM® 中基本网格描述的示意图。

读取网格需要 fvMesh 类，并配合特殊的构造函数（Listing 7.40）。include 语句用于在读取网格前完成初始化。一旦以 mesh 变量构造好 fvMesh，就可以操作网格并提取必要数据，如 Listing 7.41 所示：通过 mesh.C() 获得单元形心向量场 volVectorField C，通过 mesh.V() 获得单元体积 volScalarField V，通过 mesh.Cf() 获得面形心 surfaceVectorField Cf。

要遍历单元体积，需先调用专用函数 mesh.V()（Listing 7.42），然后用 forAll 宏进行遍历（Listing 7.43）。forAll 宏的定义如 Listing 7.44 所示，展开为 for (Foam::label i=0; i<(list).size(); i++)。其他网格信息的访问如 Listing 7.45 所示。对于边界 patch，访问过程类似：每个 patch 有自己的类，所有 patch 列表定义在 fvBoundaryMesh 类中，访问方式如 Listing 7.46：先取 mesh.boundary()，再用 forAll 遍历各 patch，对每个 patch 用 patch.Sf()、patch.magSf()、patch.nf()、patch.Cf()、patch.faceCells() 等函数访问该 patch 的具体面数据。

7.2.1 场与内存

在 OpenFOAM® 的通用框架内，可以定义不同类型与大小的列表、数组与一般容器。对于给定的网格与计算结构，定义一个能够将场、列表与向量直接与网格相关联的特定类将非常有用。满足这种需求的类是模板类 GeometricField。使用该类定义的每份数据都与网格维度严格绑定，包括边界与内部域（单元数、内部面数或内部顶点数）。

模板类 GeometricField 根据网格的三种主要特征存储以下数据结构：

  • volField：定义于单元中心的场；
  • surfaceField：定义于单元面的场；
  • pointField：定义于单元顶点的场。

该类还继承以下属性：

  • Dimensions（量纲）：OpenFOAM® 通过 GeometricField 为每个场关联一个量纲（米、千克、秒等），用以表征变量的物理含义。基于此，所有涉及 GeometricField 的运算只能在同量纲的场之间进行（例如速度只能与速度相加，不能与压力相加），否则运行时就会出错。此外，将已有场组合生成新场时，新场的量纲会由 OpenFOAM® 自动生成——具体做法是把生成新场所用的代数关系同样作用于参与运算的场的量纲上（例如质量场除以体积场，结果是密度量纲）。
  • InternalField（内部场）：一个容量等于内部网格属性（单元中心、顶点或面）大小的仓库，存储所定义场的内部信息。
  • BoundaryField（边界场）：包含所定义变量在边界上的所有相关信息。会建立一个 patch 列表，每个场为边界的 patches 定义，名字为 GeometricBoundaryField。可以在整个边界集合上操作，也可以在使用 fvPatchField 的特定 patch 上操作。
  • Mesh（网格）：作为与网格严格绑定的类，每个 GeometricField 都包含对相应网格的引用。
  • Time Values and Previous Values（时间步与前时间步值）：用于在仿真过程中处理该特定场。为达到二阶时间精度，会保存前两个时间步的信息。

下面给出使用 GeometricField 类访问场主要属性的示例。第一个示例要求定义两个变量 U（速度场）与 T（温度场），定义在网格的单元中心。这通过使用 GeometricField 的专用模板完成，模板支持标量、向量与张量数据类型。Listing 7.47 给出构造脚本（一个构造函数示例）：两个场都链接到网格，并需要指定四个参数——(i) 场名，(ii) 场量纲，(iii) 初始化值，(iv) 边界条件；这里对整套边界使用零阶外推（zeroGradient）。

类似构造函数也适用于定义在面上的变量。例如，要定义单元面处的质量通量（体积流量）场 mdot，脚本如 Listing 7.48 所示，它通过对速度在面上的插值结果与面面积向量取点积构造而成。该场表示控制体面处的体积流量；在不可压缩流且密度为 1 时，也表示质量通量场。场定义完成后，就可以通过脚本访问网格各部分的具体数据。

7.2.2 InternalField 数据

访问内部场数据的脚本如 Listing 7.49 所示：对 T.internalField() 与 U.internalField() 使用 forAll，依次取得每个 cell 的标量或向量。

7.2.3 BoundaryField 数据

访问边界场数据的脚本如 Listing 7.50 所示：先取 U.boundaryField()（类型为 GeometricBoundaryField），外层 forAll 遍历 patch，内层 forAll 遍历每个 patch 的面，取得对应的面场值。更紧凑的形式如 Listing 7.51 所示。

7.2.4 lduAddressing

OpenFOAM® 在其离散循环与系数存储中专门使用面寻址。OpenFOAM® 的网格也支持任意多面体单元，每个多面体单元可以有任意数量的面，每一张内部面对应一个邻居单元。OpenFOAM® 中系数的存储基于面寻址方案：系数按内部面顺序存储，通过与内部面关联的 owner/neighbor 索引来访问单元及其系数。owner 总是下标较小的单元，neighbor 是下标较大的单元；对边界面，owner 是该面所附的单元，neighbor 通过设其索引为 -1 表示不存在。owner 或 neighbor 的下标列表从而确定了各积分算子中单元到单元系数的装配顺序。

上述方案称为 lduAddressing，由 lduMatrix 类实现（图 7.10）。lduMatrix 包含 5 个数组，分别是对角、Upper、Lower 系数以及 owner 的 lower 索引、neighbor 的 upper 索引。在这种方案下，owner 对应矩阵的下三角部分（lower addressing），neighbor 对应上三角部分（upper addressing）。给定一面，对于 owner 单元而言，lower 与 upper addressing 分别提供该面通量系数在矩阵中存储的列与行；对于 neighbor 单元则正好相反。以图 7.10 左上域为例，内部面 4 的相关信息存储在 lower()、upper()、lowerAddr()、upperAddr() 数组的第 5 行（C++ 从 0 开始计数）。其 owner 是单元 2（lowerAddr），neighbor 是单元 4（upperAddr），在单元 2 的代数方程中乘以 φ4 的系数存储在 upper() 的第 5 行，而在单元 4 的代数方程中乘以 φ2 的系数存储在 lower() 的第 5 行。

由此 lduAddressing 提供了与面相关的非对角系数的地址信息。这意味着当矩阵运算主要基于遍历所有面时，计算效率很高；但要直接访问矩阵某个特定行列元素则较为困难且低效。一个例子是对每行求非对角系数之和：

\[ a_C = \sum_{n \neq b(C)} a_n \quad (7.1) \]

此时使用面寻址无法直接遍历每行的非对角元素，要做这种求和必须遍历所有面（因为只能按面进行），如 Listing 7.52 所示：循环遍历所有面，对 ac[l[faceI]] 与 ac[u[faceI]] 分别减去对应的 Lower[faceI] 与 Upper[faceI]。这里的 ac 是非对角系数之和，l 与 u 分别是 upper 与 lower addressing，Lower 与 Upper 是对应的系数。每行的求和不是顺序的，仅取决于网格的 owner-neighbor 编号。

总体而言，lduAddressing 引入了更复杂的矩阵操作处理过程，这一复杂性也会反映在线性求解器的实现中，但换来的是更高的计算速度。

7.2.5 计算梯度

OpenFOAM® 中用于计算 Green-Gauss 梯度的脚本如 Listing 7.53 所示。该实现位于 Foam::fv::gaussGrad<Type>::gradf 函数中：

  1. 定义 GradType 为向量与 Type 的外积类型。
  2. 获取 mesh 引用，建立临时输出场 tgGrad，其量纲为 ssf.dimensions()/dimLength，边界条件类型为 zeroGradientFvPatchField<GradType>。
  3. 取 mesh.owner() 与 mesh.neighbour()，以及面向量 mesh.Sf()。
  4. 取输出场的内部引用 igGrad 与输入场 issf。
  5. 内部面循环：对每个内部面，计算 Sfssf = Sf[facei]*issf[facei]，将其加到 owner 单元上，从 neighbor 单元上减去（同量但反号）。
  6. 边界 patch 循环：对每个 patch，取其 pFaceCells = mesh.boundary()[patchi].faceCells()、面向量 pSf = mesh.Sf().boundaryField()[patchi]、边界场 pssf = ssf.boundaryField()[patchi]；对每个 patch 内的面，将 pSf[facei]*pssf[facei] 累加到对应的 owner 单元上。
  7. 将内部累加值除以 mesh.V()，得到体积平均梯度。
  8. 调用 gGrad.correctBoundaryConditions()，返回 tgGrad。

作者强调，上述简短介绍既不是要替代 OpenFOAM® 用户手册 [1]，也不是替代 C++ 手册；其目的是向读者介绍有助于快速理解 OpenFOAM® 框架的总体思路与一般概念。尽管介绍较为概括，其中描述的语法已经涉及 OpenFOAM® 中为编写完整求解器所必需的主要数据，相关内容将在后续章节展开。

7.3 网格转换工具

有多种工具能够将各种格式的网格文件转换为 OpenFOAM® 格式，其中部分工具列举如下 [1]：

  • ansysToFoam：从 I-DEAS 导出的 ANSYS 输入网格文件转换为 OpenFOAM® 格式；
  • cfx4ToFoam：将 CFX 4 网格转换为 OpenFOAM® 格式；
  • datToFoam：读入 datToFoam 网格文件并输出 points 文件，常与 blockMesh 配合使用；
  • fluent3DMeshToFoam：将 Fluent 网格转换为 OpenFOAM® 格式；
  • fluentMeshToFoam：将 Fluent 网格转换为 OpenFOAM® 格式，可处理多域与域边界；
  • foamMeshToFluent：将 OpenFOAM® 网格写成 Fluent 网格格式；
  • foamToStarMesh：读入 OpenFOAM® 网格，写出 PROSTAR（v4）bnd/cel/vrt 格式；
  • foamToSurface：读入 OpenFOAM® 网格，将边界写成 surface 格式；
  • gambitToFoam：将 GAMBIT 网格转换为 OpenFOAM® 格式；
  • gmshToFoam：读入 Gmsh 所写的 .msh 文件；
  • ideasUnvToFoam：I-Deas unv 格式的网格转换；
  • kivaToFoam：将 KIVA 网格转换为 OpenFOAM® 格式；
  • mshToFoam：转换 Adventure 系统生成的 .msh 文件；
  • netgenNeutralToFoam：转换 Netgen v4.4 所写的中性格式；
  • plot3dToFoam：Plot3d 网格（ascii/formatted 格式）转换；
  • sammToFoam：将 STAR-CD（v3）SAMM 网格转换为 OpenFOAM® 格式；
  • star3ToFoam：将 STAR-CD（v3）PROSTAR 网格转换为 OpenFOAM® 格式；
  • star4ToFoam：将 STAR-CD（v4）PROSTAR 网格转换为 OpenFOAM® 格式；
  • tetgenToFoam：转换 tetgen 所写的 .ele、.node、.face 文件；
  • writeMeshObj：用于网格调试，将网格写成三个独立的 OBJ 文件，可用 javaview 等查看。

7.4 Closure

本章通过解释 uFVM 与 OpenFOAM® 这两个 CFD 代码的若干实现属性，概述了 FVM 在计算机代码中的实现方式：从数据结构、内存管理方案、到算例设置这两个维度对两种代码作了对比介绍。下一章将详细介绍应用于扩散通量的有限体积第二步离散。

7.5 Exercises

注：本节为章末练习题，按精读工作流约定（Constraint #13）不在内容概述范围内。

本章个人批注

本章是 ch06（理论层 mesh 数据结构）之后的工程实现篇——把 ch06 介绍过的 owners/neighbours/face-addressing 概念落到两个具体的代码框架里：教学用的 MATLAB uFVM 和工业级的 OpenFOAM® C++ 实现。两者的核心数据结构（points/faces/owners/neighbours + boundary patches）几乎一样，但表达层与抽象层差异很大，这种"理论同源、实现分叉"的对照对理解 CFD 代码非常有帮助。

对我个人而言最有价值的部分是 7.2.4 lduAddressing 那段——之前读 ch06 时只是知道 OpenFOAM® "用面寻址存系数"，但没具体想清楚它怎么存。这一章把 lduMatrix 的五个数组（diag/upper/lower + lowerAddr/upperAddr）讲清楚了，而且作者特别强调了"按面循环快、按行循环慢"这个 trade-off——owner 是下三角、neighbor 是上三角，符号约定也是在这套 ldu scheme 里统一的。这对后面读 SIMPLE/PISO 求解器的 C++ 源码会有直接的帮助。

另一个值得记下来的是 7.1.4 关于"场数组长度 = nElements + nBFaces"的设计。这种"把内部单元和边界值拼到一个连续数组里"的思路在 OpenFOAM® 的 GeometricField.InternalField/BoundaryField 里有完全相同的影子（只不过用对象层次封装得更彻底）。它本质上是用一个一维连续内存同时表达"内部 + 边界 patch 顺序"这两层语义，是 CFD 代码把"逻辑多块、物理一块"网格高效跑起来的一个核心技巧。

最后是 7.2.5 的 gaussGrad::gradf 实现——它和 7.1.6 的 cfdComputeGradientGauss0 在算法层面完全一致：累加面通量、按体积除、修正边界条件。差异主要在 C++ 的对象访问方式（mesh.Sf()、mesh.owner() 等）和模板（任意 scalar/vector/tensor 类型）。把这个对应关系记清楚，后面梯度章节再出现 Gauss 梯度公式时就不会被模板参数搞晕。

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

本章位于 Part I 的最后一章，紧接 ch06（理论层 mesh 数据结构）与 Part II（离散化）之间，起到了从"理论数据模型"到"工业代码实现"的桥梁作用。ch06 抽象地定义了 FVM 网格所需的拓扑与几何元素（点、面、单元、owner/neighbour、boundary patch），本章则把这些抽象实体分别落到 uFVM（教学）与 OpenFOAM®（工业）两套具体实现上，让读者通过对照看清"同样的 mesh 在不同抽象层级下的代码长什么样"。本章末尾的 7.4 Closure 明确预告下一章将进入 Part II 第一步离散——扩散项的空间离散（对应 ch08 Spatial Discretization: The Diffusion Term），这意味着从 ch08 起，章节重心从"几何/数据"完全转向"算子/代数"；本章的 7.2.5 gaussGrad 实际上已经为后续梯度章节做了一个小预演，lduAddressing 那段则提前为 ch10 求解代数系统（矩阵装配的 ldu 视图）埋好伏笔。