第 4 章：离散化过程（The Discretization Process）
4.1 The Discretization Process
本章以「数值求解偏微分方程」为切入，开宗明义地给出离散化过程的定义：要解一个偏微分方程，本质上就是在若干特定位置上求出因变量 ϕ 的值，再用这些离散值在求解域上重建 ϕ 的分布。这些特定位置被称为 grid elements 或 grid nodes，它们来自原始几何体被切分成「互不重叠的离散单元」的过程——这一过程称为 meshing（网格划分）。得到的节点或变量通常布置在单元质心（cell centroid）处或顶点上，具体取决于所采用的离散格式。所有数值方法的共同点，都是用离散值替代偏微分方程的连续精确解。由于 ϕ 的分布被离散化，因此把「将控制方程转换为 ϕ 离散值的代数方程组」这一过程称为离散化过程，而把完成这一转换的具体手段称为离散化方法。ϕ 的离散值一般通过一组代数方程求得，这些方程把相邻网格单元上的值联系在一起；这些离散后的（也称代数化的）方程，由 ϕ 所满足的守恒方程推导而来。ϕ 值算出后，再对数据做后处理以提取所需信息。
整章内容在图 4.1 中以「流程图」形式被概括为四个分支：从 Physical Domain 出发做 Domain Modeling（几何建模），从 Physical Phenomena 出发做 Physical Modeling（物理建模），二者汇合得到定义在 Computational Domain 上的 Set of Governing Equations；这条方程组再分两条路径——Domain Discretization（域离散化，产出网格，按结构网格 Cartesian/Non-Orthogonal、Block Structured、Unstructured、Chimera 等分类）以及 Equation Discretization（方程离散化，按方法分 Finite Difference、Finite Volume、Finite Element、Boundary Element）；两者汇合得到 System of Algebraic Equations；最后由 Solution Method（直接法、迭代法、Multigrid、Coupled/Uncoupled 等）求解出 Numerical Solutions。
为了让这些抽象概念落地，作者用一个具体算例贯穿全章：图 4.2 给出的「带铜底散热器的微处理器」传热问题。铜底作 heat spreader（热扩散片），上方接 microprocessor，下方接 heat sink。本书将沿用这一算例来引出离散化过程中的各种概念。考虑到本书的目的是介绍一种求解物理过程的数值方法，且方法要落到计算机程序中，因此整章的叙述都按「从数值到实现」的思路展开，例如会谈到内部单元、边界单元、变量等如何在代码中存储。书中将使用一套基于 Matlab® 的程序 uFVM 作为开发工具，演示这些方法在实现层面的细节；同时也会从用户和开发者两个角度介绍基于有限体积法的开源工业级代码 OpenFOAM®，主线仍是把数值方法过渡到实现细节。
4.1.1 Step I: Geometric and Physical Modeling
对物理现象的建模在某种意义上是科学事业的核心。一般而言，一个物理现象只有当它能被数学化、且该数学化能够被检验和验证时，才算真正被理解。就本书目的而言，要在两个层面上做建模：一是物理域的几何，二是感兴趣的物理现象本身。两个层面上，与研究目标无关或不必关心的细节都会被忽略或简化。例如，三维域可以简化为二维描述；可以利用对称性缩小研究域的尺寸；某些情况下，还可以把物理构件直接拿掉，代之以合适的数学描述。
在图 4.2 那个微处理器—散热器的例子里，第一步建模就要同时简化其物理与几何。Heat sink 和 processor 被替换为边界条件：前者指定 heat sink 的估算温度，后者指定 processor 的预期工作温度。物理域被建模为二维计算域，因为 heat sink 厚度方向的温度变化可以忽略。对铜底中稳态的热流和温度分布，只考虑热传导。建模过程的产物是一组线性的（或当 k 随 T 变化时为非线性的）偏微分方程，对应于能量方程的简化形式
其中 k 是 heat spreader base 的导热系数，\(\dot{q}\) 是单位体积的源/汇。
4.1.2 Step II: Domain Discretization
物理域的几何离散化的产物是网格（mesh），守恒方程最终在网格上求解。域离散化的具体操作，是把整个域切分成互不重叠的、彼此贴合的离散单元（cell 或 element），构成一套网格系统。这项工作可由多种技术完成，产出的网格类型也多种多样。这些网格按若干特征分类：结构性（structure）、正交性（orthogonality）、分块（blocks）、单元形状（cell shape）、变量排布（variable arrangement）等。无论如何，网格都由顶点（vertices）和面（faces）共同界定的离散单元组成。要让网格能作为方程离散化的有用平台，还需要网格单元的拓扑信息及一些派生的几何信息——单元—单元关系、面—单元关系、面/单元质心与面积/体积、面的法向等。这些信息通常由基础网格数据推导而来。对某些网格拓扑（如结构网格），可由单元下标直接推出拓扑细节；对另一些（如非结构网格），则必须构建并存储为列表供后续检索。
图 4.3a 展示了一个简单二维域：包含一个体积（二维下是面积）和三类边界，分别对应 microprocessor 的加热/冷却、heat sink 和 heat spreader base。域被图 4.3b 中的一组简单网格所离散；网格边界被分成三个 patch（Patch#1、Patch#2、Patch#3），各 patch 用于施加不同的物理边界条件。该网格由 25 个互不重叠的单元构成，单元几何由 40 个点（cell 顶点）定义，单元之间共有 66 个面（二维下为线段），其中 34 个是内部面。Step III 中由守恒方程离散化所得的代数方程，对计算域内每个单元都写一个，解用单元场表示，值定义在每个单元的质心处。本例中单元取正方形，其他形状（如三角形，见图 4.3c）也可以。
网格可从不同角度描述。最基本的层面，网格就是一个由一维、二维或三维空间中代表位置的点（vertices 或 points）组成的列表。网格也代表被切分成互不重叠单元的离散域，单元可以是任意凸多面体。单元由面（faces）完整地围成，这些面通常与相邻单元共享（边界处除外）。单元既可用界定它的点定义，也可用围成它的面定义。存储在列表中的面分两类：（i）内部面（interior faces）——由两个单元共享（或连接）；（ii）边界面（boundary faces）——与域边界重合，只属于一个单元。内部面可由单元拓扑信息推导；边界面则必须由用户显式提供，因为它们定义域的物理边界。二维中面由其两个端点描述；三维中端点描述 edge，由 edge 围出面。内部面法向方向一般由两侧相邻单元的拓扑关系决定；边界面法向总是指向域外。图 4.4 展示了网格的三类组成（图 a 是 vertices，b 是 faces，c 是 elements）。进一步，边界面的组织方式是按其所属 patch 归到不同的列表中。
4.1.3 Mesh Topology
离散化过程中，偏微分方程在每个单元上积分，得到一组代数方程；每条方程把该单元的变量值与其相邻单元上的值联系起来。这些代数方程再被组装到全局矩阵和向量中，每个方程的系数都按对应单元下标所在的行、列位置存放。逐单元对方程的积分称为「局部组装」（local assembly），把所有单元的贡献合并成整套方程组的过程称为「全局组装」（global assembly）。方程的离散化是「按相邻单元」推导出来的，但方程的全局组装需要把单元映射到实际下标。具体的实施细节将在后续章节展开；这里只在最基础的层面、以 connectivity list（连通性表）的形式，简单介绍元素、面、顶点的拓扑信息——它们是实施上述过程的「原料」。
元素连通性（element connectivity）把局部组装矩阵与全局矩阵联系起来，使一个单元写出的方程与其它单元的方程相容。一般会建立单元—单元、单元—面、单元—顶点三类连通性，分别把单元与相邻单元、围成它的面、界定它的顶点联系起来。以图 4.4 为例，单元 9 的连通性如图 4.5 所示。
对任意形状的单元，更高效的做法是「按面循环」来组装通量项。这种情况下，面两侧的单元信息必须随时可用——这由面连通性（face connectivity）提供。对每个面，存储共享它的两个单元以备计算。面的方向按「法向从单元 1（owner）指向单元 2（neighbour）」来定义。边界面只约束一个单元（定义为单元 1），故其法向总是指向域外。图 4.6 给出面 12 的连通性。
顶点连通性（vertex connectivity）常用于后处理和梯度计算。如图 4.7 所示，顶点连通性一般包含共享该顶点的单元列表和面列表。
局部下标与全局下标的对应关系用图 4.8 简要说明：一个五单元网格，先在单元 3 上组装出「局部」代数方程，再通过元素连通性把局部下标转换为全局下标，最终把局部矩阵的贡献放入全局矩阵。
Example 1 对一个简单网格（图 4.9，含 4 个单元）演示了元素连通性，并把它表示成全局矩阵：1 → 2, 3；2 → 1, 3, 4；3 → 1, 2, 4；4 → 2, 3。再写成 4×4 的稀疏矩阵形式
（稀疏模板示意，列出的是非零模式而非具体值——原书图 4.9 的 Example 1 解展示的就是这种「邻接表→稀疏矩阵」的过程。）
4.1.4 Step III: Equation Discretization
Step III 把控制偏微分方程转换为一套代数方程，计算域内每个单元一条。这些代数方程再被组装为全局矩阵和向量，写成
其中未知变量 T 在每个内部单元和计算域边界上定义。T 的边界值一般由指定的边界条件给出。为此，必须为 T（一般也为每个控制方程）定义一个单元场（element field）。如图 4.10 所示，单元场是一组定义在每个单元质心处的值，由「interior element field」表示——它是一个大小等于内部与边界单元总数之和的数组。
方程离散化这一步在计算域的每个单元上执行，得到一条代数关系，把该单元上的变量值与相邻单元上的值联系起来。这条代数方程由微分方程离散化得到；对于本章算例，方程就是以温度 T 为未知量的能量方程。在有限体积法中，离散化的做法是：先把微分方程在控制体（控制容积 / 单元）上积分，得到「半离散」的方程形式；再通过「在网格单元之间施加 profile（剖面假设）」近似因变量的变化，得到最终的离散形式。只有少量网格单元参与到某条离散方程中——这正是所选剖面的「分段性」的直接推论。某个网格点上的 T 只会影响其紧邻邻域内的 T 分布。当网格单元数增多时，离散方程的解应当趋近于原微分方程的精确解：网格单元越密，相邻单元之间 T 的变化越小，因而剖面假设的具体形式也变得不重要。
对同一个微分方程，可写出的离散方程形式不止一种；不过只要网格单元数足够多，所有方法给出的解应当一致。形式上的差异来自剖面假设和推导方法的不同。
下面以有限体积法对方程离散化步骤作示例：对图 4.11 中控制体 C 上的能量方程进行离散化。首先在 C 上对式 (4.1) 积分，恢复出第三章给出的积分守恒形式
然后用散度定理把体积分转换为面积分
此式即为单元 C 上的热平衡，是原偏微分方程的积分形式，不含任何近似。把面积分用「对面求和」代替，得到
其中 f 是面质心处的求积点。这一变换是「第一次近似」——式 (4.4) 的积分被数值近似为各面质心处的通量之和，这是一个二阶近似（后续章节会证明）。
把求和展开得到四项相加的形式
对图 4.12 中面 f1 而言，其面积向量和温度梯度为
其中 \(x_C\) 是单元 C 质心的 x 坐标，\(x_{F_1}\) 是单元 F1 质心的 x 坐标，\(\Delta y_{f_1}\) 是面 f1 的面积，\(\mathbf{S}_{f_1}\) 是由 C 向外的面积向量，\(\nabla T_{f_1}\) 是面 f1 质心处的温度梯度。代入后，式 (4.6) 中的第一项变为
要继续推进，需要 C 与 F1 之间 T 的剖面假设。若假设 T 在两者之间呈线性变化，则面 f1 处的 x 方向梯度可写为
于是式 (4.8) 可近似为
更一般地
其中
对其余三个面作同样推导，得到
代回式 (4.6) 得
或更紧凑地写成
其中
对域内所有单元都可类似地写出方程，得到一组代数方程，可用多种直接或迭代方法求解。以图 4.13 中的单元 C 为例，式 (4.15) 给出 \(T_C\) 与其四个邻居 \(T_{F_1}, T_{F_2}, T_{F_3}, T_{F_4}\)（全局下标为 \(T_9\)、\(T_{10}\)、\(T_4\)、\(T_8\)、\(T_{15}\)）之间的关系。对边界单元同样写出方程，所得的方程组如图 4.14 所示；其矩阵形式即式 (4.2)，其中 A 是系数矩阵，\([T]\) 是解向量，b 是无法归入 A 的项。式 (4.2) 的求解方法将在下一节给出。
最后需要指出：有限体积法在精度、鲁棒性等方面的性质会在后续章节回顾，包括更细致地考察扩散项的有限体积离散（本章只对矩形 Cartesian 网格作了推导）。
4.1.5 Step IV: Solution of the Discretized Equations
微分方程的离散化产出一组离散代数方程，必须通过求解得到 T 的离散值。这些方程的系数可能与 T 无关（即线性），也可能依赖 T（即非线性）。求解这一代数方程组的技术与离散化方法无关，代表「通向解的不同路径」。本书涉及的线性代数方程组，解的唯一性是有保证的；因此只要所用方法能给出解，就一定是所求的解。对于同一组离散方程，所有「能走到解」的方法都会给出相同的解。
代数方程组的求解方法大致分为直接法和迭代法两类。
4.1.5.1 Direct Methods
直接法对系数给定的方程组（如式 (4.2)）施加一套相对复杂的算法（比迭代法复杂），只需一次即可求得解。直接法的一个例子是矩阵求逆
只要 \(A^{-1}\) 存在，就能保证 \([T]\) 有解。然而，求逆一个 N×N 矩阵的运算量是 \(O(N^3)\)，计算代价高昂。因此实际中几乎不用求逆。线性系统有更高效的方法。对于本书关心的离散方法，A 是稀疏矩阵；对结构网格，A 还是带状矩阵。对于某些方程（如纯扩散），A 还是对称的。矩阵运算可以利用 A 的这些特殊结构来设计高效算法，这些方法将在第 10 章中讨论。
总体上，直接法在 CFD 中很少使用，因为其计算和存储开销太大。当今大多数工业 CFD 问题包含数十万单元，即便简单问题每个单元也有 5–10 个未知数，因此 A 非常大，多数直接法对此类大规模问题不再实用。进一步，A 常常是非线性的——直接法必须嵌在一个外层迭代循环中以更新 A 中的非线性项；这意味着直接法会被反复调用，计算时间更为可观。
4.1.5.2 Iterative Methods
迭代法用「猜测—修正」的过程逐步精化估计解，反复求解离散系统。下面以最朴素的 Gauss-Seidel 迭代法为例。整体求解循环可写为：
(a) 猜出 T 在域内所有网格单元上的离散值。
(b) 逐个访问网格单元，用下式更新 T
更新 \(T_C\) 时需要邻居值，按当前已知的值代入即可。因此，已经访问过的网格单元会用上更新后的 T；还没访问的，则使用旧值。
(c) 扫过整个域直到所有网格单元都被覆盖，这完成一次迭代。
(d) 检查是否满足收敛判据。判据可以是「T 在所有网格点上最大变化小于 1%」之类的阈值。满足则停止；否则回到 (b) 重复。
上述迭代过程并不保证对任意 \(a_C\) 和 \(a_{NB}\) 组合都收敛。对线性问题，只要满足 Scarborough 准则，过程就保证收敛。Scarborough 准则要求对所有网格点
且至少在一点上严格小于 1。满足该准则的矩阵具有对角占优。
Gauss-Seidel 方案的存储需求极小：只需存放 T 在网格单元上的离散值。系数可以「按需」现场计算，因为更新某点的 T 时不需要保存整套系数矩阵。这种迭代结构对非线性问题尤其合适：若系数依赖 T，可以在迭代推进时用 T 的当前值更新系数。不过，Gauss-Seidel 在实际 CFD 中也很少使用——系统规模一大，收敛速率会降到不可接受的低。第 10 章将介绍代数多重网格（algebraic multigrid）方法，用于加速迭代法的收敛速率、改进其性能。
4.1.6 Other Types of Fields
除上面介绍的单元场外，还有其他场用于不同目的。两种常见的场是 face field 和 vertex field，下面简单介绍。
face field 由定义在面中心处的值的数组组成。如图 4.15 所示，face field 包含 interior faces 的数组和各个 patch faces 的数组。face field 常用于定义面质量通量，供求解对流和流动问题时使用。
vertex field（图 4.16）的变量存放在顶点上，按 interior vertices 和 patch vertices 归组。vertex field 常用于后处理，在某些情况下也用于梯度计算。
4.2 Closure
本章概览了离散化过程，并沿途强调了开发一套 CFD 代码所需的基本要素。后续章节会逐一拆解这些要素，期间开发 uFVM 并学习工业开源 CFD 库 OpenFOAM®。接下来两章将进一步展开有限体积网格和有限体积离散化两方面的内容。
本章个人批注
本章是 Moukalled 整本书的方法论总览，把「求解一个 PDE」这件事拆成四个高度抽象的步骤：几何/物理建模 → 域离散化（出网格）→ 方程离散化（出代数方程组）→ 求解代数方程组。从「写代码」的视角看，这套四步法的好处是：步骤之间是数据流关系——前一步的产物（建模后的问题、网格、离散方程）就是后一步的输入，每一步都对应一个独立的软件模块。这也是为什么作者引入了 uFVM 和 OpenFOAM® 作为「开发载体」：他想让读者把每一步落到具体的数据结构和代码上，而不是停留在公式。Step III 末尾那 10 个公式（4.5–4.16）实际上是后续「推导一般网格上有限体积格式」的最简单原型——对矩形 Cartesian 网格、线性剖面、稳态能量方程做了最直白的演示，但读者应当意识到这套推导模式会推广到任意多面体、非正交网格、瞬态问题。
阅读时我想留意的几个点：(1) Section 4.1.3 里讲「element / face / vertex connectivity」是这一节的核心——这是把网格从「几何对象」变成「程序可操作的数据结构」的关键步骤，后续的 face loop、owner/neighbour 约定、对角占优的判定都依赖这些连通性表。OpenFOAM 里 polyMesh 那一套 owner/neighbour 数组、face-cells、cell-faces 链接表，源头就在这里。(2) Scarborough 准则的引入很早（4.1.5.2）——作者想要传达的是：离散方程的系数矩阵性质决定了迭代法能不能收敛；这一步先给出一个粗糙但通用的判据（第 10 章会用多重网格来加速收敛）。（3）章节末尾的「其他类型的场」是一个看似不起眼但很重要的伏笔：单元场用于「变量在单元质心处的值」（collocated 或 cell-centered 思路），face field 用于「通量」（通量必然在面上），vertex field 用于后处理——这套「三种场」是 OpenFOAM 字段体系的基础。综合来说，第 4 章真正想传递给读者的不是公式，而是「一个 PDE 求解器在概念上和软件上应该长什么样」。
与上下章的衔接（一段话）
第 3 章「数学描述物理现象」在守恒律的微分形式和积分形式上做了铺垫——第 4 章开头就把「积分形式」拿出来作为「有限体积法在控制体上做积分」的起点（即式 4.3、4.4），让第 3 章的抽象定义直接成为本章的「算子」。第 4 章末尾预告「接下来两章将展开有限体积网格和有限体积离散化」——这正好对应第 5 章「Mesh Generation and Manipulation」和第 6 章「Diffusion»，前者接续 4.1.2–4.1.3 讲网格的生成、变换与质量评价，后者接续 4.1.4 对方程离散化中「扩散项」的最一般推导（包括非正交、多维、变系数情形）。从全书结构看，第 4 章处于「从原理过渡到实现」的转折点：前三章都在建立数学和物理基础，从第 4 章起全书的语调转向「写代码」；到第 10 章会再次回到「求解」环节（直接法、迭代法、多重网格），把 4.1.5 中点到为止的收敛性话题展开；到全书末尾再回到「对流项离散」「压力—速度耦合」「湍流模型」等具体子主题，因此第 4 章在全书里既是入口也是目录。