第 15 章：不可压缩流体流动的计算（Fluid Flow Computation: Incompressible Flows）
15.1 主要困难（The Main Difficulty）
本章开头先重述前述章节处理一般标量输运方程时所使用的方法对于求解 Navier-Stokes 方程并不充分：需要一种算法能够处理压力-速度耦合。书上把连续性方程 (15.1) 与动量方程 (15.2) 用矩阵形式 (15.3) 重新写出，强调其左下方块为零、典型的鞍点（saddle point）特征——这种结构使压力和速度场不能用一般迭代法直接求得，因此必须另外推导一个压力方程。
作者给出的解决思路是把分块矩阵 A 分解为 LU 三角阵，其中 \(-BF^{-1}B^T\) 即 Schur 补；这正是 Patankar 与 Spalding 提出的 SIMPLE（Semi-Implicit Method for Pressure Linked Equations）算法所采用的方法。该算法把 Navier-Stokes 拆成动量方程与压力方程两个方程，按顺序分别求解：先以迭代前的压力场求解动量方程，所得速度场满足动量守恒但不一定满足质量守恒；再用该速度场构造压力方程，求解压力修正量，最后同时修正压力和速度，使新场满足质量守恒。新一轮迭代再以更新后的压力场重新开始，直到两个守恒条件都被满足。
具体地，由 (15.5)–(15.6) 形式的分裂可得四条操作步骤：解 \(F v^* = b\)、解 \(-B D^{-1} B^T p^* = -B v^*\)、更新 \(v = v^* - D^{-1} B^T p^*\)、更新 \(p = p^*\)（其中 \(F^{-1}\) 用对角元 \(D^{-1}\) 近似，星号表示当前迭代的中间量）。这一分裂即 SIMPLE 系列算法的雏形。
15.2 初步推导（A Preliminary Derivation）
作者把上节思想放到一维稳态、均匀网格上具体推导。先给出简化连续性方程 (15.7) 与动量方程 (15.8)，并以元素 C 为控制体（如图 15.1）。接下来几节依此离散动量方程、连续性方程，说明棋盘问题及其解决方法，最后给出压力修正方程与一维 SIMPLE 算法。
15.2.1 动量方程的离散（Discretization of the Momentum Equation）
对 (15.8) 在元素 C 上做体积分，用散度定理把对流与扩散项的体积分转化为面积分 (15.10)，再用单高斯点表示面通量，得到半离散形式 (15.11)–(15.12)，其中 \(F_m\) 记对流通量，\(D_m\) 记扩散通量。
按前述章节的高分辨率与中心差分方法离散对流和扩散项，可得形如 (15.13) 的代数方程
。压力项的处理则留到连续性方程离散之后。
15.2.2 连续性方程的离散（Discretization of the Continuity Equation）
将 (15.7) 在元素 C 上积分，再用散度定理转化为面积分 (15.15)；按各面通量求和得
(15.16)。
15.2.3 棋盘问题（The Checkerboard Problem）
压力梯度项的离散有两种做法。其一，单高斯体积分 + 中心差分得 (15.17)–(15.18)。其二，用散度定理转化为面积分 (15.19)–(15.20)。若对压力采用线性插值剖面，可化为只含交替网格点压力的形式 (15.21)，即压力梯度依赖 pE 和 pW 的差而非 pC 与邻居。类似地，连续性方程在常密度情形下化为 \(u_E - u_W = 0\) (15.22)，即依赖同样两个交替点。这意味着任何非物理的"之字形（zigzag）"压力-速度分布（如 p = (10, -100, 10, -100, 10)、u = (1, 10, 1, 10, 1) 模式）都会被数值格式识别为均匀场（图 15.2）：对 W、C、E 三个控制体分别积分 (15.18) 与 (15.22) 均给出零，从而无法被检出。多维情形也类似，棋盘型压力-速度场始终是方程组的解——这为引入交错网格提供了动机。
15.2.4 交错网格（The Staggered Grid）
前节问题的根源是压力与速度场的解耦。如果把各变量存储在不同位置，使压力梯度（在动量方程中）和速度（在连续性方程中）都无需插值即可计算，即可强制二者耦合。图 15.3 给出 1D 交错网格：速度 u 存储在单元面 e、w 上，压力与其他标量存储在单元形心 C 上。
在该布置下，元素 C 的离散连续性方程直接是 (15.23) \(\dot m_e + \dot m_w = 0\) 或 \(u_e - u_w = 0\)，无须插值；动量方程在面 e 上的控制体上积分得 (15.24)，其中压力梯度直接取自相邻两点 \(p_E - p_C\)，无插值。因此之字形压力-速度解将不再被数值格式"接受"——它们会被检出并被算法排除。
15.2.5 压力修正方程（The Pressure Correction Equation）
本节的推导依据 Patankar 与 Spalding 的原始 SIMPLE 论文。起始方程为 (15.25) 连续性、(15.26) 动量。给定上一步压力 \(p^{(n)}\) 与速度 \(u^{(n)}\)，先用上一步压力解动量方程得到满足动量守恒的中间速度 \(u^*\)（星号），即 (15.27)。由于压力不准确，\(u^*\) 通常不满足连续性方程。
引入修正场 \((u', p')\)，令 \(u = u^* + u'\)、\(p = p^* + p'\)（(15.28)），面通量修正 \(\dot m_f = \dot m_f^* + \rho u'_f S_{x,f}\)（(15.29)）。把它代入连续性方程得 (15.31) \(\dot m'_e + \dot m'_w = -(\dot m^*_e + \dot m^*_w)\)，即用当前质量守恒误差驱动修正场。质量流量与修正量的具体形式由 (15.32)–(15.33) 给出。
为把压力引入，把 (15.26) 改写为紧凑形式
（(15.34)，其中 (15.35) 定义 \(H_e, B_e, D_e\)）。中间速度 \(u^*\) 满足 (15.36)；把 (15.36) 从 (15.34) 中减去即得修正场方程 (15.37)
，类似地对面 w 得 (15.38)。
把 (15.33) 代入 (15.31) 并把 \(u'_e, u'_w\) 用 (15.37)–(15.38) 替换，可得扩散型压力修正方程 (15.40)；离散化得 (15.42)，最终整理为标准五点形式 (15.43) \(a'_{p_C} p'_C + a'_{p_E} p'_E + a'_{p_W} p'_W = b'_{p_C}\) ，系数为 (15.44)。其中带下划线的 \(H_e(u')\)、\(H_w(u')\) 项在收敛时为零，对最终解无影响；不同近似这些项的方式就区分出不同的算法，原始 SIMPLE 直接将其忽略。对一维常截面情形 \(\Delta y\) 可取 1 并从方程中略去。
15.2.6 交错网格上的 SIMPLE 算法（The SIMPLE Algorithm on Staggered Grid）
SIMPLE 算法利用动量方程与压力修正方程，按以下顺序迭代：(1) 给定压力与速度初猜 \(p^{(n)}, u^{(n)}\)；(2) 求解动量方程 (15.27) 得到 \(u^*_f\)；(3) 用动量满足的速度场更新质量流量 \(\dot m^*_f\)；(4) 用新质量流量求压力修正方程得到 \(p'\)；(5) 按 (15.45) 修正压力、速度与质量流量，使新场满足连续性；(6) 把 \(u^{**}, p^*\) 作为新一轮初始值；(7) 回到步骤 2 直至收敛。这种顺序求解而非联立求解的方式称为 segregated approach。
Example 1 — 管网流动：图 15.4 给出一段水管网，各管段质量流量与压降的关系为 \(\dot m = D \Delta p\)，其中 \(D_A=0.5, D_B=D_F=0.4, D_C=D_E=0.3, D_D=0.19, D_G=0.1875, D_H=0.35\) ，外加固定压力 \(p_1=400, p_2=350, p_4=50, p_5=300, p_7=80, p_9=200\) 与固定流量 \(m_I=50\)。作者假设初值 \(p^{(n)}_3=300, p^{(n)}_6=200, p^{(n)}_8=120\)，按动量方程算出各质量流量，再检验节点 3、6、8 的连续性残差分别为 24、14、7（非零）。写出压力修正方程 (系数矩阵非零，分别对应三个内点) 并求解，得 \(p'_3 = p'_6 = p'_8 = -20\)。更新质量流量与压力后，第二轮节点 3、6、8 的连续性残差全部为 0，压力修正方程右端项也都为 0，迭代一轮即收敛。
15.2.7 二维交错笛卡尔网格上的压力修正方程（Pressure Correction Equation in Two Dimensional Staggered Cartesian Grids）
在二维情形需三套网格：u-分量、v-分量、压力（及其他标量）各占一套。把 1D 推导推广到二维即可。对图 15.6 的元素 C，压力修正方程为 (15.46)，系数见 (15.47)。注意 N、S 方向的系数含 \(\Delta x_C\)（与 E、W 方向含 \(\Delta y_C\) 对偶）。
Example 2：图 15.7 中给 \(u_w=50, v_s=20, p_N=0, p_E=10\)，密度 1，动量方程 \(u_e = -d_e(p_E-p_C), v_n=-d_n(p_N-p_C)\)，其中 \(d_e=1, d_n=0.25\)，网格 \(\Delta x=\Delta y=1\)。设 \(p^{(n)}_C=100\)，由动量方程得 \(u^*_e=90, v^*_n=25\)，即 \(\dot m^*_e=90, \dot m^*_n=25\)。连续性残差为 \(90 - 50 + 25 - 20 = 45\)。在二维例题中由于 \(\dot m^*_e + \dot m^*_n - \dot m^*_w - \dot m^*_s = 45 \ne 0\) ，压力修正方程为 \(1.25 p'_C = -45\)，故 \(p'_C = -36\)。修正后 \(u^{**}_e = 90 - 36 = 54\), \(v^{**}_n = 25 - 0.25 \cdot 36 = 16\), \(p^{**}_C = 64\)。第二轮连续性残差为 0，压力修正方程右端为 0，一步收敛。
15.2.8 三维交错笛卡尔网格上的压力修正方程（Pressure Correction Equation in Three Dimensional Staggered Cartesian Grid）
图 15.8 给出三维元素：u、v、w 分别存储在 (e, w)、(n, s)、(t, b) 面上。压力修正方程扩展为 7 点形式 (15.48)，系数 (15.49)。
15.3 交错网格的缺点（Disadvantages of the Staggered Grid）
虽然交错网格支撑了 SIMPLE 的诞生，但它有显著缺陷：
- 在二维、三维中分别需要 3、4 套交错网格（每套对应一个速度分量），既占内存又增加几何信息管理开销。
- 对非笛卡尔特别是非结构网格，staggering 本身变得不自然；图 15.9 表明在曲线网格中某一面可能与交错速度方向对齐，使速度梯度计算失效。
- 替代方案有协变（covariant）或逆变（contravariant）曲线分量（图 15.10a, b），但曲线坐标系下动量方程的离散会引入额外的非守恒项，扩散项处理也更复杂。
- 另一种思路是把每个方向上的笛卡尔分量交错到所有面（图 15.12），使每面都有三个速度分量，但 2D 下动量方程数翻倍、3D 下三倍，对非结构网格更是没有明显的"交错方向"，只能改变压力与速度所用网格大小，或把所有速度分量交错到所有面，变量数爆炸。
- 因此转向 cell-centered collocated（同心型）网格（图 15.13）：所有变量存储于同一形心，而质量通量（标量）存于面；它可视为一种特殊的 contravariant 分量，只是其值通过 Rhie-Chow 插值自定义计算——下一节即讨论。
15.4 Rhie-Chow 插值（The Rhie-Chow Interpolation）
1983 年 Rhie 和 Chow 提出一种插值方法使 SIMPLE 可在 collocated 网格上运行。其思想是在线性插值得到的单元面速度之上加一项"耗散项"，该耗散项等于两种压力梯度估计之差（图 15.14）。这两种估计基于不同的网格 stencil。该方法等价于在面 f 上构造一个 pseudo-momentum 方程，其系数由相邻形心 C、F 的动量方程系数线性插值得到，压力梯度则用小 stencil 计算，从而模拟了交错网格的紧邻压力-速度耦合。
从单元 C、F 的 x-动量离散方程 (15.50) 出发，要写面 f 上的类似方程 (15.51)，但其系数不能直接计算，故用线性插值得到 (15.52)：\(H_f[u] = (H_C[u] + H_F[u])/2\)，\(B_f, D_f\) 同理。代入即得 (15.53)，即重建后的 pseudo-momentum 方程。
书中用 (15.54) 记所有带横杠的量为线性插值 \(h_f = g_C h_C + g_F h_F\)（几何因子 \(g_C, g_F\) 与 f 相对 C、F 的位置相关）。把 (15.52) 代入 \(H_f\) 得 (15.55)；二者之差（(15.56)）被证明是二阶精度。把 (15.55) 代回 (15.53) 即得面速度 (15.57)：
15.5 一般推导（General Derivation）
在推导多维 collocated 压力修正方程之前，先给出多维动量方程的离散形式。
15.5.1 动量方程的离散（The Discretized Momentum Equation）
从动量方程 (15.64) 出发，三处下划线分别是非稳态项、对流项、扩散项；按前述章节方法离散。应力张量的第二部分（速度梯度转置）通过散度定理 (15.65) 转化为面积分 (15.66)；压力梯度按 (15.67) 或 (15.68) 处理；体力按 (15.69) 积分。
用一阶 Euler 离散非稳态项、HR 格式（延迟修正）离散对流项、扩散通量拆为隐式对齐项 + 显式交叉扩散项，最终得向量形式的离散动量方程 (15.70)，其中系数 (15.71) 由面通量 (15.72) 和单元通量 (15.73) 构成。
虽然代数形式 (15.70) 是线性的，但其系数依赖于 v 与 p 自身——这种非线性用迭代处理，每步开始用上一轮值计算系数。这种系数大幅变化会使 v 产生大变动，影响收敛甚至发散。为此可采用 under-relaxation，引入松弛因子 \(\lambda_v\)（Patankar 隐式松弛）后得 (15.74)，重新定义系数为 (15.75)，可得简洁形式 (15.76)。为推导 collocated 压力修正方程，把压力梯度从源项中显式提出，得 \(b_v^C = -V_C (\nabla p)_C + \hat b_v^C\)（(15.77)），代入 (15.76) 得 (15.78)。定义向量算子 (15.79) 后，最终得到动量方程的简洁形式 (15.80)：\(v_C + H_C[v] = -D_v^C (\nabla p)_C + B_v^C\)，后文推导会反复使用。
15.5.2 Collocated 压力修正方程（The Collocated Pressure Correction Equation）
与交错网格类似，从迭代初值 \((v^{(n)}, \dot m^{(n)}, p^{(n)})\) 出发，先解动量方程 (15.80) 得满足动量守恒的中间速度场 \(v^*\)（(15.81)）。该速度场不满足连续性，因为压力、速度用了上一轮值。引入修正场 \((v', p', \dot m')\)，(15.82) 给出精确场与中间场的关系。
把 (15.82) 代入 (15.25) 连续性方程得 (15.83)
，其中面速度用 Rhie-Chow 插值 (15.84) 计算。若
（质量守恒），则右端为 0、修正场为 0；否则需要非零修正。
把速度修正写成与压力修正的关系：用 (15.80) 减去 (15.81) 得 (15.85) \(v'_C + H_C[v'] = -D_v^C (\nabla p')_C\)，对面 F 类似 (15.86)；面质量流量修正 (15.87) 中面速度修正为 (15.88)。代入 (15.83) 后得 (15.89)，下划线部分代表邻居速度修正对当前速度修正的影响。把 (15.85)–(15.86) 内插到面得 (15.90)，代回 (15.89) 即得压力修正方程的扩散型形式 (15.91)，显式展开为 (15.92)。
下划线项 \(H_f[v']\) 的处理决定了算法的可解性：SIMPLE 直接把它忽略（不影响终解只影响收敛路径）。其余项按第 8 章各向异性扩散方法处理。把 \(D_v^f (\nabla p')_f \cdot S_f\) 改写为
（(15.93)–(15.94)），其中 \(S'^f_f = D_f \cdot S_f\)。再把 \(S'^f_f\) 分解为正交分量 \(E_f\) 与切向分量 \(T_f\)：(15.95)–(15.96)。\(T_f\) 项（来自非正交网格）可被忽略（不影响终解）或保留（按 OpenFOAM 的非正交内部循环处理）；保留时该项需要在每轮迭代中更新。
忽略非正交贡献，线性化后 (15.97) \(D_f (p'_F - p'_C)\)。代回 (15.91) 得代数形式 (15.98)，系数见 (15.99)。原始 SIMPLE 即把 (15.99) 中下划线项直接忽略。质量流量 \(\dot m^*_f\) 在 (15.99) 中用 Rhie-Chow 插值 + 最新速度场算得 (15.100)。
解出 \(p'\) 后，按 (15.101) 同时更新压力、速度与质量流量（注意只对压力做 under-relaxation，速度和质量流量的更新因 \(p'\) 已保证质量守恒而不再松弛）。
15.5.3 \(D_f\) 项的计算（Calculation of the \(D_f\) Term）
(15.96) 中的 \(E_f\) 可按第 8 章任一正交修正方法构造，从而给出不同的 \(D_f\) 表达式。
15.5.3.1 Minimum Correction Approach：把 \(S'^f_f\) 代入 (8.68) 得 \(E_f = e_{CF} \cdot S'^f_f \cdot e_{CF}\)（(15.102)），结合 (15.94)、(8.64) 得 (15.103)，最终 (15.104)：
15.5.3.2 Orthogonal Correction Approach：由 (8.69) 得 \(E_f = |S'^f_f| e_{CF}\)（(15.105)），结合 (15.94)、(8.64) 得 (15.106)。
15.5.3.3 Over-Relaxed Approach：由 (8.64)、(8.70) 得 \(E_f = (S_f \cdot S'^f_f / (d_{CF} \cdot S'^f_f)) d_{CF}\) （(15.107)），结合 (15.94) 得 (15.108)。三种方法任选其一即可。
15.5.4 Collocated SIMPLE 算法（The Collocated SIMPLE Algorithm）
按图 15.16 流程：
- 在时间步 \(t+\Delta t\)，以 \(t\) 时刻收敛值 \((m^{(n)}, v^{(n)}, p^{(n)})\) 为初猜。
- 解动量方程 (15.70) 得 \(v^*\)。
- 用 Rhie-Chow 插值 (15.100) 更新面质量流量 \(\dot m^*_f\)。
- 装配压力修正方程 (15.98) 并解得 \(p'\)。
- 用 (15.101) 修正 \(p, v, \dot m\) 到满足连续性的 \(p^*, v^{**}, \dot m^{**}\)。
- 把新值作为下一轮初猜，回到步骤 2 直至收敛。
- 把 \(t+\Delta t\) 时刻的解设为收敛值；推进到下一时间步；回到步骤 1。
Example 3：图 15.17 中 \(p_W=100, p_N=20, p_E=50\)，s 面入口 \(v_s=20\)、零压力梯度，密度 1。动量方程 \(u_C = -d_x(p_e-p_w), v_C = -d_y(p_n-p_s)\)，\(d_x=1, d_y=0.25\)，\(\Delta x=\Delta y=1\)，初猜 \(p^{(n)}_C=70\)。由零梯度得 \(p^{(n)}_s = p^{(n)}_C = 70\)。计算插值压力
。代入动量方程得 \(u^*_C = -1(60-85)=25, v^*_C = -0.25(45-70)=6.25\)。用 Rhie-Chow 计算面通量 \(\dot m^*_e = 20, \dot m^*_n = 12.5, \dot m^*_w = -30, \dot m^*_s = -20\) 。连续性残差 \(20 + 12.5 - 30 - 20 = -17.5\)。构造压力修正方程 \(2.25 p'_C - p'_E - 0.25 p'_N - p'_W = 17.5\)。若 E、N、W 修正均为零，则 \(p'_C = 7.78\)，\(p^{**}_C = 77.78\)；\(v^{**}_C = v^*_C + 0.9725 = 7.2225\)。
15.6 边界条件（Boundary Conditions）
至少一个面在边界 patch 上的元素称 boundary element，其边界面的处理对精度与鲁棒性至关重要。压力基算法能否成功取决于动量方程与压力修正方程边界条件是否正确实施。
Rhie-Chow 插值在边界面上需修改：边界面上不存在像内面那样的两侧平均，平均化为 \(\bar h_b = h_C\)（(15.109)，b 表边界）。代入得 (15.110)
，相应的 Rhie-Chow 插值为 (15.111)，质量流 (15.112)。下面先讲动量边界条件，再讲压力修正边界条件。
15.6.1 动量方程的边界条件（Boundary Conditions for the Momentum Equation）
半离散动量方程 (15.113) 在边界元素上须按边界类型修改面通量项 (15.114)–(15.116)。压力项无论按元素离散或面离散，展开后都需边界值，故统一采用面离散形式以反映边界压力对解的影响。
15.6.1.1 Wall Boundary Conditions
No-Slip Wall：速度满足 \(v_b = v_\text{wall}\)（静止壁为 0）。这并非简单 Dirichlet 条件，因为除了速度等于壁面速度，还要保证法向通量为零并考虑剪应力；图 15.19 表明需保证剪应力沿切向并使壁面速度等于边界速度。
壁面作用于流体的力 (15.117) 可分解为法向 \(F_\perp\)（应为 0）与切向 \(F_\parallel\)；因此 \(F_b = s_\text{wall} S_b\)（(15.118)），其中
（(15.119)）。\(v_k\) 是 v 在切向的分量，\(d_\perp\) 是 C 形心到壁面法向距离 (15.120)，\(v_k\) 用 (15.121) 由 v 减去其法向分量得到。把 (15.121) 代入 (15.119) 得 (15.122)，最终层流边界力为 (15.123)。u、v、w 三分量的动量系数修正分别见 (15.124)–(15.126)。
边界压力 \(p_b\) 由内场外推：(15.127) \(p_b = p_C + \nabla p^{(n)}_C \cdot d_{Cb}\)（Taylor 截断）；或由 Rhie-Chow 插值（质量流和速度在壁面为零）得 (15.128)–(15.129)，化简为 (15.130)；写成 \(S'^b_b = E_b + T_b\) 后 (15.131)，利用 \(E_b\) 在 Cb 方向得 (15.132)–(15.133)。
Slip Wall：壁面剪应力为零（图 15.20），故 \(F_b = 0\)。\(p_b\) 仍按 (15.127) 或 (15.133) 计算。动量系数修改为 (15.134)：\(a^0_{vC}\) 不变，\(b_v^C\) 减去 \(p_b S_b\)。
15.6.1.2 Inlet Boundary Conditions
三种类型：(i) 规定速度，(ii) 规定静压 + 速度方向，(iii) 规定总压 + 速度方向。
Specified Velocity：\(v_b, \dot m_b\) 已知，\(p_b\) 用 (15.135) 外推。边界速度相关项作为显式源加到 \(b_v^C\)，并把 \(a^0_{vF=b}\) 加到 \(a_v^C\) 上 (15.136)。
Specified Pressure and Velocity Direction：\(p_b\) 已知，直接用于计算 \(\nabla p_C\)；\(\dot m_b\) 由连续性给出；给定方向单位向量 \(e_v\)，\(\dot m_b = \rho_b |v_b| e_v \cdot S_b\) 得 (15.137)，然后按 (15.136) 修改系数。
Specified Total Pressure and Velocity Direction：总压
（(15.138)）。\(\dot m_b\) 由连续性给出；用 (15.137) 求 \(|v_b|\) 后按速度已知情形处理（系数按 (15.136) 修改）。
15.6.1.3 Outlet Boundary Conditions
Specified Static Pressure：假定充分发展即出口速度沿面向量方向梯度为零，相当于零通量条件，系数修改 (15.139)。为保证只在出流方向梯度为零，用
（(15.140)），然后用 Taylor 展开 \(v_b = v_C + \nabla v_b \cdot d_{Cb}\)（(15.141)），系数修正 (15.142)。
Specified Mass Flow Rate：不可压缩情形下等价于指定法向速度；设方向与 C 相同 \(e_b = e_C\)，\(|v_b| = \dot m_b / (\rho_b e_C \cdot S_b)\)（(15.144)–(15.143)），按指定速度处理（系数按 (15.136)）。
Fully Developed Outlet Flow：假定法向梯度为零，\(v_b\) 用 (15.140)–(15.141) 算得作为已知；\(p_b\) 由 (15.145) 外推；系数按 (15.142) 修改。
15.6.1.4 Symmetry Boundary Condition
对称面上标量反射，法向梯度为零；对速度向量则反射后法向分量为 0、切向分量保持，即
（(15.146)）。法向分量 (15.147)；力 \(F_b\) 分解为 \(F_\parallel = 0, F_\perp = \sigma_\perp S_b\)（(15.148)），其中法向应力 (15.149)，最终 \(F_b\) 为 (15.150)。压力梯度沿法向为零：(15.151)，故
（(15.152)），\(p_b = p_C + \nabla p_b \cdot d_{Cb}\)（(15.153)）。u、v、w 三分量系数修改分别见 (15.154)–(15.156)。该处理同样适用于粘性流中的滑移壁条件。
本节未涵盖全部动量边界条件，但已覆盖最常用类型。
15.6.2 压力修正方程的边界条件（Boundary Conditions for the Pressure Correction Equation）
边界元素上连续性方程 (15.157) 写为 (15.158)，其中 \(\dot m^*_b\) 和 \(\dot m'_b\) 需按边界类型计算。对内面 \(\dot m_f, \dot m'_f\) 由 (15.100)、(15.101) 定义；对边界 b 仅边界单元贡献平均值，即 (15.109) 代入 (15.100)、(15.101) 得 (15.159)。
实现边界条件时需计算 \(\dot m^*_b, \dot m'_b, p_b, p'_b\)。三种类型：(i) 规定质量流量（如壁面或速度入口），\(\dot m'_b = 0\)，无需修改压力修正方程，但 \(p_b\) 需由内场计算；(ii) 规定压力，\(p'_b = 0\)，需对压力修正方程施加强制 Dirichlet；(iii) 隐式耦合关系（如总压入口），需把隐式关系显式化后代入压力修正方程。
15.6.2.1 Wall Boundary Condition
无论无滑移（图 15.19）还是滑移（图 15.20），壁面 \(\dot m_b = 0\)，故 \(\dot m'_b = 0\)——等价于规定零通量，压力修正方程无需修改。壁面 \(p_b\) 用 (15.127) 或 (15.133) 或低阶外推 (15.160) 计算。
15.6.2.2 Inlet Boundary Conditions
Specified Velocity：\(\dot m_b\) 已知故 \(\dot m'_b = 0\)，与壁面相同直接从压力修正方程中删除该项；\(p_b\) 用 (15.127) 或 (15.133) 或 (15.160) 外推。
Specified Pressure and Velocity Direction：\(p_b\) 已知故 \(p'_b = 0\)，但 \(\dot m'_b \ne 0\)。入口作为 Dirichlet 条件，\(p'\) 方程系数按 (15.161) 加上 \(q_b D^v_C\) 项。
Specified Total Pressure and Velocity Direction：总压 (15.138) 改写为关于 \(\dot m_b\) 的形式 (15.162)；由 Taylor 展开 (15.163) 与对 (15.162) 求偏导得 \(p'_b\) 表达式 (15.164)。代入 (15.159) 得 \(\dot m'_b\)（(15.165)），把它代入连续性方程 (15.158) 得 \(a^0_{pC}\) 的修正 (15.166)。
15.6.2.3 Outlet Boundary Conditions
Specified Pressure：\(p'_b = 0\)，但 \(\dot m'_b\) 按 (15.167) 计算；速度方向按惯例取为上风速度 \(v_C\) 方向，\(a^0_{pC}\) 修正为 (15.168)。
Specified Mass Flow Rate：\(\dot m'_b = 0\)，直接从压力修正方程中删除，不需修改系数；由 (15.159) 知 \(\dot m'_b = 0\) 意味着 \(p'_b = p'_C\)。
Fully Developed Outlet Flow：\(v_b\) 由零法向梯度视为已知，\(\dot m_b\) 已知故 \(\dot m'_b = 0\)；\(p_b\) 由内场外推。由于速度是迭代更新的，这种处理不能保证除收敛外的整体守恒。惯用做法是修正出口 \(\dot m_b\) 使整体守恒 (15.169)：先把入流总流量 \(\sum \dot m_\text{in}\) 与出流总流量 \(\sum \dot m_\text{out}\) 算出，按 (15.169) 把各出流 \(\dot m_\text{out}\) 重新调整，但要求出口远离回流区。
15.6.2.4 Symmetry Boundary Condition
对称面上质量流量为零，故 \(\dot m'_b = 0\)，与壁面相同直接删除该项；\(p_b\) 由 (15.127) 或 (15.133) 或 (15.160) 外推。
15.6.2.5 The Relative Nature of Pressure
对不可压缩流、所有边界上规定法向速度的情形，存在一个困难：动量方程中只含压力梯度，没有办法确定压力的绝对水平，只有压力差有物理意义；这使系数矩阵 A 奇异，直接求解失败。解决办法是在域内某点把压力设为规定值，其余压力即相对于该值计算。
15.7 SIMPLE 系列算法（The SIMPLE Family of Algorithms）
SIMPLE [13] 把压力、速度用 segregated 方式处理：通过构造压力修正方程、把连续性方程中的速度场替换为压力项实现耦合。在推导中 \(H_f[v']\) 项被忽略（保留则方程难以处理）；这一忽略不影响终解（收敛时为 0）但影响收敛路径——该项较大时会导致发散或缓慢收敛。为抵消，SIMPLE 对压力做 under-relaxation：\(p = p^* + \lambda_p p'\)，其中 \(\lambda_p\) 通常取 \(1 - \lambda_v\)（见 15.8）。
尽管用松弛，SIMPLE 的收敛率仍依赖问题。研究者提出改进方案，最终形成 SIMPLE-like 系列：SIMPLEC [17]、SIMPLER [3]、PISO [18]、SIMPLEX [5]、PRIME [19]、SIMPLEM [20]、SIMPLEST [21]。Moukalled 与 Darwish [22] 统一了不可压与可压缩情形下的这些算法的公式；Darwish 等 [23]、Jang 等 [24] 评估了它们的性能。本章重点讨论 SIMPLEC（Van Doormal 与 Raithby）和 PISO（Issa）——前者用加权平均估计速度修正，把 \(H_f[v']\) 项改写为幅度更小的 \(\tilde H_f[v']\) 后忽略；后者则通过 split operator 把 \(H_f[v']\) 部分回收。其他算法（包括 SIMPLEM、SIMPLER、SIMPLEX、SIMPLEST）也都将 \(H_f[v']\) 忽略，但通过对动量方程或 \(D_v\) 算子做不同修改实现改进。其中 PISO 等价于一步 SIMPLE + 一步或多步 PRIME，因此后者也一并讨论。
PRIME（Pressure Implicit Momentum Explicit）[19] 把动量方程显式求解。这一显式处理的理由是动量方程的迭代对整体流场收敛贡献小，而压力场的正确求解对整体收敛影响最大。SIMPLEST [21] 把动量方程系数拆为对流与扩散两部分：对流显式、扩散隐式，从而改变 \(D_v\) 与 \(H\)。SIMPLEM（SIMPLE-Modified）[20] 把压力修正方程放在动量方程之前求解——用旧速度场求压力，使压力修正优于速度修正，从而交换了 SIMPLE 的优劣。SIMPLER（SIMPLE-Revised）[3] 额外推导一个方程直接计算压力，再用 SIMPLE-like 压力修正方程更新速度——单独求压力是因为用预测压力修正更新速度后该速度不再满足动量方程，需另求压力匹配速度以重新满足动量。SIMPLEX [5] 旨在使收敛率不随网格加密而下降，与 SIMPLE 的差异在 \(D_v\) 的计算方式：它用额外一组方程（基于压力差的空间分布在网格加密下变化不大的假设）解得 \(D_v\)。所有这些算法最初都在 segregated 网格上推导，但也适用 collocated 网格。
15.7.1 SIMPLEC 算法（The SIMPLEC Algorithm）
SIMPLEC（SIMPLE-Consistent）[17] 是 SIMPLE 的修正版本，假设 C 点速度修正是相邻点速度修正的加权平均。数学上即 (15.170)，用 H 算子写成 (15.171)：\(H_C[v'] \approx v'_C H_C[1]\)。
因此不是忽略 \(H_C[v']\)，而是用 (15.171) 中的近似值代入，把速度修正 (15.85) 写成 \((1 + H_C[1]) v'_C = -D_v^C (\nabla p')_C\)，得 (15.172)：\(v'_C = -\tilde D_v^C (\nabla p')_C\)。这一修正使压力无需松弛，且速度修正比 SIMPLE 更好地满足动量方程——因而收敛率更高。除把忽略项改为 \(\tilde H_C[v' - v'_C]\)、把 \(D_v\) 改为 \(\tilde D_v\) 外，SIMPLEC 的步骤与 SIMPLE 相同。
15.7.2 PRIME 算法（The PRIME Algorithm）
PRIME（PRessure Implicit Momentum Explicit）[19] 把动量方程显式求解：
（(15.176)）。显式处理的理由是：动量方程的迭代对整体流场收敛贡献小，而压力场的正确求解对整体收敛影响最大。用 (15.177) 定义的修正场，修正后满足 (15.178)。减一下可得 (15.179)，代入连续性方程得 (15.180)。
PRIME 中忽略的项 \(H_C[v^* - v^{(n)}] + H_C[v']\) 在 \(H_C[v']\) 与 \(H_C[v^* - v^{(n)}]\) 符号相反时会比 SIMPLE 忽略的 \(H_C[v']\) 小：\(H_C[v'] = H_C[v^{**} - v^*]\) 是为满足连续性而作的修正，\(H_C[v^* - v^{(n)}]\) 是为满足动量而作的修正——一般二者反向，因而 PRIME 忽略的总项幅度更小。又因为动量方程是显式求解，无需 under-relaxation，从而提高算法稳定性。
15.7.3 PISO 算法（The PISO Algorithm）
PISO [18, 25] 把 \(H_C[v']\) 项通过两步或多步修正过程回收。第一步与 SIMPLE 相同——从 (15.83) 求解 \(v'\) 并忽略 \(H_C[v']\)，得到满足连续性的 \(v^{**}, p^*\)；用它们重新计算动量方程系数并显式求解得新速度 \(v^{***}\)；用 Rhie-Chow 算 \(\dot m^{***}_f\)。第二步修正：速度修正写成 (15.181)，可整理为 (15.182)。用 (15.181) + Rhie-Chow 在 C、F 间插值，得到新的压力修正方程 (15.183)，其中下划线项再次忽略——该修正步可重复多次以逐步回收更多 \(H_C[v']\)。
由此 PISO 可看作一步 SIMPLE + 一步或多步 PRIME 的组合，结合了 SIMPLE 的隐式性与 PRIME 的稳定性。Collocated PISO 算法流程（图 15.27）：
- 在 \(t+\Delta t\) 时刻以 \(t\) 时刻收敛值 \((p^{(n)}, u^{(n)}, \dot m^{(n)})\) 为初猜。 SIMPLE Step：2. 隐式求解动量方程 (15.70) 得 \(v^*\)；3. 用 Rhie-Chow (15.100) 更新面通量 \(\dot m^*_f\)；4. 装配压力修正方程 (15.98) 并解得 \(p'\)；5. 用 (15.101) 修正得满足连续性的 \(m^{**}, v^{**}, p^*\)。 PRIME Step(s)：6. 用最新速度、压力场计算动量方程系数并显式求解；7. 用 Rhie-Chow 更新面通量；8. 装配压力修正方程 (15.183) 并求解；9. 类似 (15.101) 修正 \(p, v, \dot m\)；10. 回到 6 重复所需次数；11. 把 \(v^{**}, \dot m^{**}, p^*\) 作为下一轮初猜；12. 回到 2 直到收敛。
- 设置 \(t+\Delta t\) 时刻解为收敛值；推进时间步；回到 1 重复直至最后时间步。
15.8 v 与 p' 的最佳欠松弛因子值（Optimum Under-Relaxation Factor Values for v and p0）
为促进 SIMPLE 收敛，对动量与连续性方程分别用 \(\lambda_v, \lambda_p\) 做欠松弛。速度修正直接为 (15.184) \(v'_C = -D_v^C (\nabla p')_C\)（无松弛）；压力修正则松弛以使 (15.184) 满足精确速度修正方程
（(15.185)）。
联立两式得 (15.186) \(\lambda_p = 1 + H_C[v'] / v'_C\)。SIMPLEC 消除了对压力松弛的需要并给出最佳加速率；用 SIMPLEC 的近似把 C 点速度修正写为相邻点速度修正的加权平均 (15.187)。从 (15.70)–(15.73) 知 \(a_v^C\) 由 (15.188) 给出，稳态极限下收敛于 (15.189)
。代入 (15.187) 得 (15.190)，再代入 (15.186) 即得 (15.191) \(\lambda_p \approx 1 - \lambda_v\)。经验表明 SIMPLE 算法在欠松弛因子满足 (15.191) 时性能与 SIMPLEC 相当。
15.9 Rhie-Chow 插值对各项的处理（Treatment of Various Terms with the Rhie-Chow Interpolation）
15.9.1 欠松弛项的处理（Treatment of the Under-Relaxation Term）
在 collocated 网格 + Rhie-Chow 插值下，解会依赖于动量方程欠松弛因子 \(\lambda_v\)。为消除这一依赖，需修正 Rhie-Chow 插值。把欠松弛动量方程写为 (15.192)（其中 \(b_v^C\) 是把压力与欠松弛源项剥离后的源项，\(v^{(n)}_C\) 是上一迭代值）；交错网格下的对应形式为 (15.193)。
Rhie-Chow 插值通过构造面伪动量方程模拟交错网格形式——这正是它能成功的原因。因此指导原则是：任何对 Rhie-Chow 的修正都应使其形式与交错网格形式相似 (15.194)。把 (15.194) 右端首项求平均 (15.195)，代入 (15.194) 得修正后的面速度 (15.196)：
。不修正欠松弛对界面速度的影响会使解依赖于欠松弛因子。
15.9.2 非稳态项的处理（Treatment of the Transient Term）
用 backward Euler 离散非稳态项时，离散动量方程为 (15.197)；交错网格对应 (15.198)。用 Rhie-Chow 构造面伪动量方程 (15.199)，求平均 (15.200)。代入 (15.84) 得 (15.201)：
。不修正非稳态项对界面速度的影响会使解依赖于时间步长且在小步长时出现振荡。该修正仅对一阶 Euler 离散有效；对更高阶时间格式可按相同原则类推。
15.9.3 体力项的处理（Treatment of the Body Force Term）
交错网格下体力项的 stencil 与压力梯度项完全一致（图 15.28）。collocated 网格下体力、速度、动量都在同一位置计算——为使体力离散保留与压力相同的 stencil，需对体力项重分配。离散动量方程为 (15.202)（双横杠表示两步平均）：先在面 f 上由 (15.203) \(B^f_v = g_C B^C_v + (1-g_C) B^F_v\) 算面体力；再把这些面值在形心上求平均（图 15.29）。
为推导单元中心上的体力平均值，按 [26] 考察图 15.30 的一维静止流体：压力梯度与体力平衡 (15.204) \(0 = -\nabla p_f + B^f_v\)，故 C、F 间压力关系为 (15.205)–(15.206)，体力大小 \(B_f = \rho_f g\)（(15.207)）。对单元 C 同样有 (15.208) \(\nabla p_C = B^C_v\)。但压力梯度用面积分 (15.209) 计算时把 (15.206) 代入并化简（(15.210)）即得 (15.211)：
第二个要求是面速度满足交错网格方程形式 (15.212)（其中 \(b_v^C\) 是从 (15.71) 中剥离压力与体力后的源项）。把系数求平均 (15.213) 代入 (15.179)，得修正后的面速度 (15.214)：
，其中 \(B^f_v\) 按 (15.215) 计算。这一额外处理提高了对体力变化重要（如自由面流）情形的鲁棒性。
15.9.4 欠松弛、非稳态与体力项的合并处理（Combined Treatment of Under-Relaxation, Transient, and Body Force Terms）
三者应一并处理。利用叠加原理可得到界面速度的完整修正 (15.216)：
15.10 计算实现（Computational Pointers）
15.10.1 uFVM
uFVM 中压力修正方程在脚本 cfdAssembleMdotTerm 中实现。Listing 15.1 给出核心算法：压力修正方程的系数由对各内面的通量线性化组装；mdot 场（面质量流率）按 (15.216) 计算，分解为 9 项（I–IX）逐步组装得到面速度：
- I: 插值速度场 \(\bar v_f\)；
- II, III: 压力梯度的修正项
； - IV, V: 体力的平均与重分配 \(\bar D_v^f (B^f_v - \bar B^f_v)\)； - VI, VII: 非稳态通量 \((\bar a_f \bar D_v^f / V_f)(\bar v_f - v^{(n)}_f)\)； - VIII, IX: 松弛修正项 \((1 - \lambda_V)(\bar v^{(n)}_f - v^{(n)}_f)\)。
清单最后把所有项组装为全局通量并输出 \(\dot m_f\)、矩阵系数等。
15.10.2 OpenFOAM®
下面用前述数值技术在 OpenFOAM® [27] 上开发不可压缩 Navier-Stokes 求解器。
15.10.2.1 Pressure Correction SIMPLE Solvers
基于 SIMPLE 算法构造多个求解器：
simpleFoam（非 OpenFOAM 自带）：基础版，仅含 SIMPLE 最基本形式；simpleFoamImproved：扩展基础代码以改进松弛处理；simpleFoamTransient：在稳态 simpleFoam 上加瞬态能力；simpleFoamBuoyancy：加体力处理；simpleFoamCompressible：可压缩版（第 16 章）；simpleFoamTurbulent：湍流处理（第 17 章）。
simpleFoam 首先给出。Listing 15.2 显示 #include 宏定义——fvCFD.H 包含构建任何 OpenFOAM 应用所需的类定义，createMesh.H 实例化网格类（Listing 15.3）。然后在 createFields.H 中定义压力场 p（Listing 15.4），同时定义压力修正场 pp（Listing 15.5，构造时把 p 的边界类型复制到 pp），并把所有 fixedValue 边界置零（Listing 15.6）。速度场 U 由输入文件读入，质量通量 mDot 由插值速度构造（Listing 15.7）。热物理属性仅指定运动粘度 nu（Listing 15.8）。
SIMPLE 主体循环：解动量方程 (15.217)–(15.219)（Listing 15.9；用 fvm::div(mDot, U) 表示对流，fvm::laplacian(nu, U) 表示扩散）；新速度一般不满足连续性，需装配压力修正方程 (15.220)（(15.221) 定义 D），Listing 15.10 用 fvm::laplacian(DUf, pp) 实现；DUf 由 Listing 15.11 通过线性插值得到；Listing 15.12 用 Rhie-Chow 计算面质量流 mDot；Listing 15.13 解压力修正方程；Listing 15.14 用 ppEqn.flux() 更新 mDot 使其满足连续性，flux 函数简化版见 Listing 15.15；Listing 15.16 修正压力（用 URF 即 \(\lambda_p\)）与速度。
simpleFoamImproved 把 Rhie-Chow 扩展以考虑速度松弛：按 (15.196) 在 mdot 中加松弛项（Listing 15.17，需定义 mdotf.prevIter() 与 U.prevIter()）。
simpleFoamTransient 进一步扩展以考虑非稳态项：按 (15.201) 在 mdot 中加非稳态项（Listing 15.18）；主循环用 pimpleControl 加瞬态外循环（Listing 15.19）。
simpleFoamBuoyancy 在 simpleFoamTransient 基础上加：(i) 能量方程（Listing 15.20）、(ii) 动量方程中的体力源项（Listing 15.21）、(iii) Rhie-Chow 插值中考虑体力重分配（Listing 15.22）。
15.11 小结（Closure）
本章介绍了 collocated 网格上求解不可压缩流问题的 segregated 压力基方法。Rhie-Chow 插值之所以能在 collocated 网格上成功，是因为它在面 f 上构造一个 pseudo-momentum 方程，其压力梯度 stencil 与交错网格的紧邻 stencil 相同。本章还详细讨论了动量方程与压力修正方程中最常见边界条件的实现细节。下一章将把压力基方法推广到全速域可压缩流体流动。
15.12 习题（Exercises）
（本章习题涵盖 1D 管网、1D 多孔介质流、1D 不可压缩、2D 交错网格 SIMPLE、1D Forchheimer 多孔介质、Rhie-Chow 面速度手算以及 OpenFOAM® 的 SIMPLEC/PISO 实现练习。本节不展开具体内容。）
本章个人批注
本章是 Moukalled 等教材的核心之一，把压力-速度耦合问题从"为什么需要压力修正方程"讲到了"怎么在 OpenFOAM 里写出来"，对计算流体力学的工程实现有提纲挈领的作用。读完一遍后我个人最深的感受是：作者用同一条主线——交错网格 → 压力修正方程推导 → Rhie-Chow → 算法族（SIMPLEC/PRIME/PISO）——把"为什么这样设计"和"如何落地"串在了一起，避免了只讲操作而不见来龙去脉。
关于 15.2.3 棋盘问题：作者清楚地把"两种压力梯度离散方式（体积分 vs 面积分）"最终都化为了"只含 pE 和 pW 的差"这一本质结果，并用 p = (10, -100, 10, -100, 10) 的具体数值示例说明这种解在多维上也可能被错认为均匀场。这是全章最直观的物理论证之一。
15.2.5 推导压力修正方程时，作者把 \(H_e(u')\) 等"含速度修正的修正项"明确标注为下划线——它们在收敛时为零、不影响最终解。不同算法（SIMPLE、SIMPLEC、PISO、PRIME）的差异本质上就在于对这一项的处理策略：忽略（SIMPLE）、加权后忽略（SIMPLEC）、通过 split operator 多次回收（PISO）、作为动量-连续性两修正之差而部分抵消（PRIME）。把握住这条主线，再看 15.7 一节就不会被各种算法变体淹没。
15.6 节 197 行的压力修正边界条件虽然占了大半篇幅，但其实是把 (15.161)、(15.166)、(15.168) 三种 \(a^0_{pC}\) 修正公式反复应用，思路是一致的：规定的量（\(\dot m_b\) 或 \(p_b\)）直接置零其修正，未规定的量则由连续性或外推给出。15.6.2.5 关于"压力相对性"的提醒虽然简短，但点出了一个不可压缩算法普遍需要处理的奇异性问题。
15.10 的 OpenFOAM 实现其实是对前 9 节理论的一次"反推"——看到 fvm::laplacian(DUf, pp) 与 (15.98) 的对应、ppEqn.flux() 与 (15.101) 修正通量的对应，理论就不再抽象。
未覆盖的源文习题部分（15.12）按 skill 规则跳过，不在内概述中展开。
与上下章的衔接（一段话）
第 14 章刚结束对时间项（包括非稳态格式与 Source Term Relaxation）的讨论；本章把焦点转到了不可压缩 Navier-Stokes 的压力-速度耦合这一更困难的问题上。从结构上看，本章建立在前面所有章节之上：第 8 章的扩散项离散为 \(D_f\) 的最小/正交/超松弛三种算法提供了 15.5.3 的直接素材；第 9 章的非正交梯度公式被引为 \(\nabla p_f\) 计算依据；第 13、14 章的非稳态与松弛格式则在本章 15.9 的 Rhie-Chow 扩展里被显式纳入。这一衔接既回顾又预告：作者用一整章把"压力基算法在 collocated 网格上能不能工作、怎么工作"讲透，作为后续第 16 章可压缩流（所有速度均可压）与第 17 章湍流模型（同一 SIMPLE/SIMPLEC/PISO 框架上叠加输运方程）的过渡。