跳转至

第15章 流体计算：不可压缩流（Fluid Flow Computation: Incompressible Flows）

作者

本章由 F. Moukalled、L. Mangani 和 M. Darwish 合著。本章是全书 Algorithms 模块（Ch 15-16）的开篇，专门处理不可压缩 N-S 方程的数值求解。N-S 方程是流体动力学的"圣杯"，不可压缩情形下"无压力方程" + "压力-速度强耦合"两大难题，使 SIMPLE 系列算法成为工业 CFD 的事实标准。

内容概述

本章 145+ 页篇幅是全书 20 章中最长的章节，覆盖了不可压缩 N-S 方程求解的完整算法体系：SIMPLE、SIMPLEC、PRIME、PISO，以及 Rhie-Chow 插值、动量插值、压力边界条件等关键技术。

§15.1 主要困难 —— 不可压缩 N-S 方程的两大核心难题：

  1. 压力-速度强耦合 —— 动量方程中的 \(\nabla p\) 项使速度依赖压力；但压力不显式出现在连续性方程 \(\nabla \cdot \mathbf{v} = 0\) 中。
  2. 无压力方程 —— 没有直接可解的"压力 PDE"。

在数学上表现为"鞍点问题"（saddle point problem）：线性化后

\[\begin{pmatrix} A & B^T \\ B & 0 \end{pmatrix} \begin{pmatrix} v \\ p \end{pmatrix} = \begin{pmatrix} f_b \\ 0 \end{pmatrix} \tag{15.3}\]

其中 0 在右下角，使得该系统无法用标准的迭代方法（Jacobi、Gauss-Seidel）直接求解。Schur 补 \(-B A^{-1} B^T\) 的存在使得必须"分块求解"，这就是 SIMPLE 算法的数学基础。

§15.2 预备推导 —— 在 1D 均匀网格上演示 SIMPLE 算法的推导。

  • 1D 连续性 + 动量方程 (Eq. 15.7-15.8)

$\(\frac{\partial (\rho u)}{\partial x} = 0 \tag{15.7}\)$

$\(\frac{\partial (\rho u u)}{\partial x} = \frac{\partial}{\partial x} \left( \mu \frac{\partial u}{\partial x} \right) - \frac{\partial p}{\partial x} \tag{15.8}\)$

  • 1D 离散动量方程 —— 形成 \(a_C u_C + \sum_F a_F u_F = b_C - (\partial p / \partial x)_C V_C\)
  • 1D 压力修正方程推导 —— 通过"先动量、后连续性"两步迭代

§15.3 SIMPLE 算法 (Semi-Implicit Method for Pressure Linked Equations) —— Patankar-Spalding 1972 提出的经典算法。

  • 算法步骤 (Algorithm 15.1)：
  • 设定初始压力场 \(p^* = p^n\)
  • 求解动量方程得到 \(u^*, v^*, w^*\)
  • 求解压力修正方程 (Pressure Correction Equation) 得到 \(p'\)
  • 修正压力 \(p^{n+1} = p^* + p'\)
  • 修正速度 \(u^{n+1} = u^* + u'\)，其中 \(u' = -\frac{V_C}{a_C} \nabla p'\)
  • 求解其他标量（湍流量、温度、组分）

  • 返回步骤 2 直到收敛

  • 压力修正方程推导 —— 在 collocated 网格上，\(u' = -\frac{1}{a_C} \frac{\partial p'}{\partial x}\)（沿面需要 \(\mathbf{d}\) 与 \(\mathbf{S}_f\) 的几何因子修正）。

  • 完整压力修正方程 (Eq. 15.20)

$\(\sum_f \frac{\rho V_P}{a_P} (\nabla p')_f \cdot \mathbf{S}_f = -\sum_f \rho \mathbf{v}_f^* \cdot \mathbf{S}_f\)$

右侧是"质量流量残差"，左侧是"压力修正扩散项"。

  • SIMPLE 的不足 —— 压力修正采用"完全隐式"，导致压力更新过于保守，需要较小的亚松弛因子（\(\alpha_p = 0.3\)）。这带来收敛慢的问题。

§15.4 SIMPLEC 算法 (SIMPLE Consistent) —— van Doormaal-Raithby 1984 改进。

  • 核心改进 —— 在推导 \(u'\) 时"精确"计算动量系数 \(a_C\)，避免 SIMPLE 中的"完全隐式"近似。
  • 结果 —— 压力更新更准确，可以用更大的亚松弛因子（\(\alpha_p = 0.7\)），收敛速度提升 2-3 倍。
  • OpenFOAM 默认算法 —— simpleFoam 默认使用 SIMPLEC。

§15.5 PRIME 算法 (Pressure Implicit with Momentum Explicit) —— Raithby-Schneider 1979 提出。

  • 核心思想 —— 显式处理动量项、隐式处理压力项
  • 优势 —— 收敛速度比 SIMPLE 更快
  • 局限 —— 稳定性差（动量显式处理），需要更小的 \(\Delta t\)
  • 应用 —— 较少用于现代工业 CFD，OpenFOAM 也无内置支持

§15.6 PISO 算法 (Pressure-Implicit with Splitting of Operators) —— Issa 1986 提出。

  • 核心思想 —— 把 SIMPLE 的"压力修正 - 速度修正"循环执行 1-2 次额外迭代
  • 优势 —— 对瞬态 N-S 方程特别有效，是 OpenFOAM pimpleFoam 的核心
  • 算法步骤 (Algorithm 15.4)：
  • 动量预测 \(u^*\)
  • 压力修正 1（SIMPLE 步骤）
  • 速度修正 1
  • 压力修正 2（额外一次）
  • 速度修正 2

  • （可选）压力修正 3

  • PISO 的优势 —— 对每个物理时间步 \(\Delta t\)，通过 1-2 次额外的压力修正使解更准确，因此整体时间步长可以更大。

§15.7 Rhie-Chow 插值 —— Ch 7 提到的 collocated 网格"棋盘压力场"问题解。

  • 问题 —— 在 collocated 网格上，\(p, u, v\) 都存于 cell 形心，压力梯度 \(\nabla p\) 的中心差分会产生"棋盘状"解（checkerboard pressure field）。
  • 解决 —— Rhie-Chow 1983 提出的"动量加权插值"（Momentum Weighted Interpolation）：

$\((\mathbf{v}_f)^* = \overline{(\mathbf{v})}_f - \overline{\left( \frac{V}{a} \right)}_f (\nabla p)_f\)$

其中 \(\overline{(\cdot)}_f\) 是面 \(f\) 上的简单线性插值，\((\nabla p)_f\) 是用 Ch 9 的方法计算的面梯度。关键是 \(\overline{(V/a)}_f\) 是基于动量系数的"加权"插值，而不是简单的几何插值。

  • 效果 —— 通过动量加权项，抑制了"棋盘压力场"模式，使 collocated 网格 + SIMPLE 算法稳定收敛。

§15.8 动量插值 (Momentum Interpolation) —— 进一步推广 Rhie-Chow 思想。

  • 普通线性插值 (Eq. 15.60)

$\((\nabla p)_f = g_P (\nabla p)_P + (1 - g_P) (\nabla p)_{NB}\)$

  • 动量插值 (Eq. 15.61)

$\((\nabla p)_f = \left( \frac{V}{a} \right)_f^{-1} \left[ \overline{ \left( \frac{V}{a} \right) } (\nabla p)_P + \overline{ \left( \frac{V}{a} \right) } (\nabla p)_{NB} \right] \cdot (\mathbf{S}_f / S_f)\)$

在边界处需要特别处理（边界 \(\partial p / \partial n\) 已知 vs 未知）。

§15.9 边界条件 —— 不可压缩 N-S 方程的常见边界条件： - 入口 (Inlet) —— 给定 \(u, v\) (Dirichlet)，压力外推 (zeroGradient) - 出口 (Outlet) —— 给定 \(p\) (Dirichlet)，速度外推 (zeroGradient) - 壁面 (Wall) —— 速度无滑移 \(u = 0\)（Dirichlet），压力 zeroGradient - 对称面 (Symmetry) —— 法向速度 0、切向速度外推

§15.10 Rhie-Chow 扩展 —— 把 Rhie-Chow 插值推广到瞬态、亚松弛、体积力项。

§15.11 计算指针 —— uFVM 与 OpenFOAM 的实现： - OpenFOAM simpleFoam（稳态不可压）：SIMPLEC + Rhie-Chow - OpenFOAM pimpleFoam（瞬态不可压）：PISO + Rhie-Chow - uFVM cfdSolveSimple(...) / cfdSolvePiso(...) 调用

§15.12 闭包 —— 总结 SIMPLE 算法家族的核心思想。

核心方程与概念

  • 不可压缩 N-S 方程 (Eq. 15.1-15.2)

$\(\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0 \tag{15.1}\)$

$\(\frac{\partial (\rho \mathbf{v})}{\partial t} + \nabla \cdot (\rho \mathbf{v} \mathbf{v}) = -\nabla p + \nabla \cdot (\mu \nabla \mathbf{v}) + \rho \mathbf{g} \tag{15.2}\)$

对不可压牛顿流体，\(\rho = \text{const}\)，故 \(\nabla \cdot \mathbf{v} = 0\)。

  • 鞍点问题 (Eq. 15.3)

$\(\begin{pmatrix} A & B^T \\ B & 0 \end{pmatrix} \begin{pmatrix} v \\ p \end{pmatrix} = \begin{pmatrix} f_b \\ 0 \end{pmatrix} \tag{15.3}\)$

\(A\) 是动量方程离散矩阵，\(B\) 是连续性方程离散矩阵，\(B^T\) 是压力梯度算子。

  • Schur 补 (Eq. 15.4)

$\(A = \begin{pmatrix} F & B^T \\ B & 0 \end{pmatrix} = \begin{pmatrix} F & 0 \\ B & -BF^{-1}B^T \end{pmatrix} \begin{pmatrix} I & F^{-1}B^T \\ 0 & I \end{pmatrix} = LU \tag{15.4}\)$

其中 \(-BF^{-1}B^T\) 是 Schur 补矩阵。

  • SIMPLE 算法的矩阵分裂 (Eq. 15.5-15.6)

$\(\begin{pmatrix} I & D^{-1}B^T \\ 0 & I \end{pmatrix} \begin{pmatrix} v^* \\ p^* \end{pmatrix} = \begin{pmatrix} v^n \\ p^n \end{pmatrix} \tag{15.5}\)$

$\(\begin{pmatrix} F & 0 \\ B & -BD^{-1}B^T \end{pmatrix} \begin{pmatrix} v^* \\ p^* \end{pmatrix} = \begin{pmatrix} f_b \\ 0 \end{pmatrix} \tag{15.6}\)$

其中 \(D^{-1}\) 是 \(F^{-1}\) 的对角近似。

  • 压力修正方程 (Eq. 15.20)

$\(\sum_f \frac{\rho V_P}{a_P} (\nabla p')_f \cdot \mathbf{S}_f = -\sum_f \rho \mathbf{v}_f^* \cdot \mathbf{S}_f \tag{15.20}\)$

左侧是压力修正的二阶导数，右侧是"质量流量残差"。

  • 速度修正 (Eq. 15.22)

$\(u'_P = -\frac{V_P}{a_P} (\nabla p')_P \tag{15.22}\)$

即"用压力梯度修正速度"。

  • Rhie-Chow 插值 (Eq. 15.58)

$\(\mathbf{v}_f^* = \overline{\mathbf{v}}_f - \overline{\left( \frac{V}{a} \right)}_f (\nabla p^*)_f \tag{15.58}\)$

其中 \(\overline{\mathbf{v}}_f\) 是简单线性插值，\(\overline{(V/a)}_f\) 是动量系数加权的插值。

  • 动量插值 (Momentum Interpolation) (Eq. 15.61)

$\((\nabla p)_f = \overline{(\nabla p)}_f - \left( \overline{ \left( \frac{V}{a} \right) }_f - \frac{V_P}{a_P} g_P - \frac{V_{NB}}{a_{NB}} g_{NB} \right) (\nabla p)_\text{face linear}\)$

该式在 Rhie-Chow 基础上进一步修正面梯度，确保压力-速度解耦的彻底抑制。

关键结论

  1. SIMPLE 算法家族是工业不可压缩 N-S 求解的"事实标准"：SIMPLE、SIMPLEC、PRIME、PISO 各有优势，工业 CFD 代码（OpenFOAM、ANSYS Fluent、STAR-CD）都基于这些算法。
  2. Rhie-Chow 插值是 collocated 网格 + SIMPLE 算法的"必备配套"：没有 Rhie-Chow，collocated 网格会产生"棋盘压力场"使解物理无效。
  3. 压力修正方程是 SIMPLE 算法的"核心"：把"无压力方程"问题转化为"通过连续性残差构造压力修正"问题。
  4. SIMPLEC 比 SIMPLE 收敛快 2-3 倍：因为压力更新更准确，可以用更大的亚松弛因子。OpenFOAM simpleFoam 默认采用 SIMPLEC。
  5. PISO 是瞬态不可压缩 N-S 的"工业标配"：通过每个时间步的 1-2 次额外压力修正，使解更准确、整体时间步长可以更大。OpenFOAM pimpleFoam 默认采用 PISO。
  6. 压力亚松弛 \(\alpha_p = 0.3\) 是工业典型：压力更新要慢，否则容易发散。SIMPLEC 可以用 \(\alpha_p = 0.5-0.7\)。

挑战和开放性问题

  1. 全耦合求解 (Fully Coupled Solver) —— 相比 SIMPLE 的"分块迭代"，全耦合方法（如 block-coupled、Newton-Krylov）在大规模问题上收敛更快，但实现复杂。OpenFOAM 6+ 引入 preconditioner + multiField 求解器。
  2. Schur 补预条件子的工程实现 —— SIMPLE 的数学基础是 Schur 补，但工业实现中很少显式使用 Schur 补预条件子。Block preconditioner 是更工程化的选择。
  3. 高 \(Re\) 流动的 SIMPLE 收敛 —— 当 \(Re > 10^5\) 时，SIMPLE 收敛极慢，需要特殊处理（低 \(Re\) 初始化、多重网格、Vanka smoother）。
  4. Vanka smoother for SIMPLE —— Vanka smoother 是基于 SIMPLE 矩阵结构的 block smoother，在多重网格加速 SIMPLE 中常用。本书未展开。
  5. 多物理场耦合的 SIMPLE 扩展 —— 流-热-化-力多场耦合时，SIMPLE 算法需要扩展（如 SIMPLE-T for 热、SIMPLE-Chem for 化学）。
  6. GPU 加速的 SIMPLE —— 工业级 CFD 的 GPU 加速是 2016 年后的研究热点，AMG-based pressure solver 的 GPU 实现是关键。
  7. 机器学习辅助的 SIMPLE 加速 —— 现代研究用 NN 加速压力修正方程的求解，本书未涉及。

个人反思与批判性分析

本章是全书最长（145+ 页）、最难、最重要的章节，是工业 CFD 求解器设计的"圣经"。

从写作特点看：

  • 从 1D 推导到 3D 非结构 —— 严格遵守"由简到繁"的工程递进：1D staggered 网格（§15.2）→ 1D collocated 网格（§15.3）→ 多维 collocated 非结构网格（§15.7）。读者可以在每一步都看到"前一步如何推广到下一步"。
  • SIMPLE 算法的"分块迭代"本质 —— 作者明确指出 SIMPLE = "动量预测 + 压力修正 + 速度修正 + 标量更新"4 步循环。这种"分块"结构是工业 CFD 求解器的核心架构。
  • 4 大算法 (SIMPLE / SIMPLEC / PRIME / PISO) 的并列展开 —— 每种算法的"动量 - 压力 - 速度"循环步骤都详细给出，使读者能够直接对照实现。
  • Rhie-Chow 插值的"动量加权"思想 —— 作者清晰解释了为什么普通线性插值会产生"棋盘压力场"，而动量加权插值能抑制它。这一节是 Ch 7 "collocated 网格"承诺的兑现。
  • OpenFOAM simpleFoam / pimpleFoam 的对应 —— 把抽象的算法映射到具体的 OpenFOAM 求解器，是工程实用性的体现。

但本章也存在以下不足：

  • 多重网格加速 SIMPLE 在工业大规模算例中是核心，本书未深入。OpenFOAM 的 GAMG 是 SIMPLE 压力求解的工业实现，但 Ch 10 仅给出 AMG 算法，本章未详细说明 AMG 与 SIMPLE 的结合。
  • Block-coupled / 全耦合求解 是 2010 年代后的研究热点，本书未涉及。
  • 动量插值（Momentum Interpolation）的推导较为简略，读者需要参考 Rhie-Chow 1983 原始文献。
  • 多物理场耦合的 SIMPLE 扩展 （如 SIMPLE-T for 共轭传热、SIMPLE-Chem for 化学反应）未展开。
  • 非牛顿流体的 SIMPLE —— 黏度 \(\mu(\nabla \mathbf{v})\) 依赖速度梯度，导致动量方程高度非线性。工业上需要专门的"非牛顿 SIMPLE"（如 SIMPLER-NN）。
  • 机器学习加速 是 2020 年代后的研究热点，本书出版年代限制未涉及。

总体而言，本章是 FVM 不可压缩 N-S 求解的"权威教材式"参考章节，与 Patankar 1980、Ferziger-Peric、Issa 1986 的对应章节形成对标。Moukalled 的优势在于"理论（鞍点问题 + Schur 补） + 工程（4 大 SIMPLE 算法 + Rhie-Chow） + 实现（OpenFOAM simpleFoam/pimpleFoam）"三位一体的连贯叙述。本章是工业 CFD 求解器开发者的"必读圣经"。

重要参考文献

  • [X1] Patankar S.V., Spalding D.B. (1972) A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. International Journal of Heat and Mass Transfer, 15(10): 1787-1806. （SIMPLE 算法原始文献）
  • [X2] van Doormaal J.P., Raithby G.D. (1984) Enhancements of the SIMPLE method for predicting incompressible fluid flows. Numerical Heat Transfer, 7(2): 147-163. （SIMPLEC 原始文献）
  • [X3] Issa R.I. (1986) Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of Computational Physics, 62(1): 40-65. （PISO 算法原始文献）
  • [X4] Rhie C.M., Chow W.L. (1983) Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA Journal, 21(11): 1525-1532. （Rhie-Chow 插值原始文献）
  • [X5] Raithby G.D., Schneider G.E. (1979) Elliptic systems: finite-difference method II. Handbook of Numerical Heat Transfer (Minkowycz et al. eds.). Wiley. （PRIME 算法）
  • [X6] Ferziger J.H., Perić M. (2002) Computational Methods for Fluid Dynamics, 3rd Edition. Springer. （SIMPLE 算法深入参考）
  • [X7] Versteeg H.K., Malalasekera W. (2007) An Introduction to Computational Fluid Dynamics: The Finite Volume Method, 2nd Edition. Pearson.
  • [X8] Moukalled F., Darwish M. (2001) A unified formulation of the segregated class of algorithms for the incompressible Navier-Stokes equations. Computers & Fluids, 30(5): 561-588. （SIMPLE 系列算法的统一矩阵推导）
  • [X9] Jasak H. (1996) Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows. PhD Thesis, Imperial College. （OpenFOAM 早期开发者的博士论文，Rhie-Chow 实现的原始设计）
  • [X10] Weller H.G., Tabor G., Jasak H., Fureby C. (1998) A tensorial approach to computational continuum mechanics using object-oriented techniques. Computers in Physics, 12(6): 620-631. （OpenFOAM 早期论文）
  • [X11] OpenFOAM Foundation (2014) The Open Source CFD Toolbox — User Guide. （OpenFOAM simpleFoam / pimpleFoam 实现）
  • [X12] Vanka S.P. (1986) Block-implicit multigrid solution of Navier-Stokes equations in primitive variables. Journal of Computational Physics, 65(1): 138-158. （Vanka smoother 原始文献）
  • [X13] Hutchinson B.R., Raithby G.D. (1986) A multigrid method based on the additive correction strategy. Numerical Heat Transfer, Part A, 9(5): 511-537. （多重网格加速 SIMPLE）