跳转至

第10章 线性代数方程组求解（Solving the System of Algebraic Equations）

作者

本章由 F. Moukalled、L. Mangani 和 M. Darwish 合著。本章是 Ch 5-9 离散化结果的"求解阶段" —— 把 FVM 离散化得到的稀疏线性系统 \(\mathbf{A} \mathbf{x} = \mathbf{b}\) 通过直接法或迭代法求解。本章是 FVM 离散化与 CFD 求解器的"桥梁"，详细展开 Gauss 消元、LU 分解、Thomas 算法（三对角）、Jacobi、Gauss-Seidel、SOR、共轭梯度（CG）、GMRES、BiCGSTAB、AMG（代数多重网格）等主流算法。

内容概述

本章 100+ 页篇幅是全书 20 章中数学密度最高的章节之一。FVM 离散化得到的线性系统是大型稀疏对角占优矩阵，矩阵规模 = 控制体数（典型工业算例 10⁶-10⁹）、非零元素数 = 控制体数 × 平均邻接数（典型 5-10 倍控制体数）。这种"大规模 + 稀疏"特性决定了求解算法的选择范围。

§10.1 引言 —— \(\mathbf{A} \mathbf{x} = \mathbf{b}\) 的矩阵形式 (Eq. 10.2)。CFD 求解器面临的两大挑战：(a) 矩阵规模大（10⁶+ 行）；(b) 非线性（\(\mathbf{A}\) 本身依赖 \(\mathbf{x}\)，每步迭代需更新）。

§10.2 直接法 (Direct Methods) —— Gauss 消元、LU 分解、Thomas 算法（三对角）、Pentadiagonal 算法。 - Gauss 消元 (Eq. 10.3-10.7) —— 通过前向消元 + 回代求解。 - LU 分解 \(\mathbf{A} = \mathbf{L} \mathbf{U}\) —— 把矩阵分解为下三角 + 上三角，再两次前向 / 回代求解。 - Thomas 算法 (TDMA) —— 三对角矩阵的 \(\mathcal{O}(N)\) 求解，是 1D 问题的标准算法。 - Pentadiagonal —— 5 对角矩阵的 \(\mathcal{O}(N)\) 求解，用于 2D / 3D 结构网格。

直接法在 CFD 中极少使用（除非极小规模），原因：(a) 大规模稀疏矩阵的 LU 分解会填入（fill-in）大量非零元素，破坏稀疏性；(b) 非线性问题的 \(\mathbf{A}\) 随 \(\mathbf{x}\) 变化，每步迭代都需重新 LU 分解。

§10.3 基本迭代法 (Basic Iterative Methods) —— Jacobi、Gauss-Seidel、SOR（Successive Over-Relaxation）。 - Jacobi —— \(\mathbf{x}^{k+1} = \mathbf{D}^{-1} (\mathbf{b} - (\mathbf{L} + \mathbf{U}) \mathbf{x}^k)\)。每步更新使用上一次迭代的 \(\mathbf{x}^k\)，并行性好。 - Gauss-Seidel —— \(\mathbf{x}^{k+1} = (\mathbf{D} + \mathbf{L})^{-1} (\mathbf{b} - \mathbf{U} \mathbf{x}^k)\)。每步立即使用新值 \(\mathbf{x}^{k+1}\)，收敛更快但顺序依赖。 - SOR —— 引入松弛因子 \(\omega \in (0, 2)\) 加速收敛。\(\omega > 1\) 为 over-relaxation，\(\omega < 1\) 为 under-relaxation。

§10.4 Krylov 子空间方法 —— 工业 CFD 的"标准求解器"。

  • 共轭梯度法 (CG) —— 适用于对称正定（SPD）矩阵。热传导方程、压力 Poisson 方程（不可压流 SIMPLE 算法的压力修正方程）属于此类。
  • 预处理 CG (PCG) —— 通过预处理矩阵 \(\mathbf{M}\) 加速：求解 \(\mathbf{M}^{-1} \mathbf{A} \mathbf{x} = \mathbf{M}^{-1} \mathbf{b}\)。常见预处理：Jacobi 预处理、SSOR 预处理、ILU 预处理。
  • GMRES (Generalized Minimal Residual) —— 适用于非对称矩阵。N-S 方程动量方程属于此类。
  • BiCGSTAB —— 双共轭梯度稳定法，适用于非对称矩阵，比 GMRES 内存少但稳定性略差。
  • 不完全 LU 分解 (ILU) —— 保留稀疏性 + 部分填入的 LU 分解作为预处理。

§10.5 代数多重网格 (AMG) —— 现代工业 CFD 求解器的"标配"加速器。多重网格的核心思想：在粗网格上消除低频误差，在细网格上消除高频误差，从而全局收敛。AMG 的"代数"特性使其无需几何信息，适合非结构网格。

  • 几何多重网格 (GMG) vs 代数多重网格 (AMG) —— GMG 需要规则网格层（结构网格），AMG 仅需矩阵（适合非结构）。
  • V-cycle / W-cycle / F-cycle —— 多重网格的不同循环方式。V-cycle 最常用，W-cycle 收敛更慢但更稳健。
  • GAMG (Geometric-Algebraic Multigrid) —— OpenFOAM 的多重网格实现，结合了几何与代数方法。

§10.6 矩阵预条件 (Preconditioning) —— 详细介绍各类预条件子：Jacobi、ILU、SSOR、AMG、Block preconditioner（多物理场耦合问题）。

§10.7 计算指针 —— uFVM 与 OpenFOAM 的具体实现： - OpenFOAM 通过 system/fvSolution 配置：solver PCG / GAMG / smoothSolver / PBiCG，preconditioner DIC / DILU / FDIC / GAMG，tolerance 1e-6，relTol 0.01。 - uFVM 通过 cfdSolveEquation(A, b, x, solverType, ...) 调用。

§10.8 收敛性诊断 —— 残差下降、收敛速度、停滞判断、亚松弛因子调节。

§10.9 闭包 —— 总结线性求解器选择规则。

核心方程与概念

本章是"线性代数 + 数值分析"的密集章节，按"直接法 → 迭代法 → 预条件 → 多重网格"四层列举。

一、矩阵与稀疏性（§10.1）

  • 线性系统 (Eq. 10.1)

$\(\mathbf{A} \mathbf{x} = \mathbf{b} \tag{10.1}\)$

\(\mathbf{A}\)：稀疏对角占优矩阵 \(\mathbf{x}\)：待求向量（控制体形心变量） \(\mathbf{b}\)：源项向量

  • 稀疏度

$\(\text{sparsity} = \frac{\text{非零元素数}}{N^2} \sim \mathcal{O}(1/N)\)$

对 \(N = 10^6\) 算例，sparsity ~ \(10^{-6}\)，即 \(\mathbf{A}\) 中 99.9999% 元素为零。

二、直接法（§10.2）

  • Gauss 消元 (Eq. 10.7) —— 对 \(N \times N\) 矩阵，\(\mathcal{O}(N^3)\) 操作，\(\mathcal{O}(N^2)\) 内存。
  • LU 分解 \(\mathbf{A} = \mathbf{L} \mathbf{U}\) —— 一次分解 + 多次回代，对 \(m\) 个右端向量总成本 \(\mathcal{O}(N^3/3 + m N^2)\)。
  • Thomas 算法 (TDMA) —— 对三对角矩阵 \(\mathcal{O}(N)\) 操作，\(\mathcal{O}(N)\) 内存。1D 问题的标准算法。

三、基本迭代法（§10.3）

  • Jacobi 迭代 (Eq. 10.31)

$\(x_i^{k+1} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij} x_j^k \right) \tag{10.31}\)$

收敛条件：\(\mathbf{A}\) 严格对角占优 \(\sum_{j \neq i} |a_{ij}| < |a_{ii}|\)。

  • Gauss-Seidel 迭代 (Eq. 10.35)

$\(x_i^{k+1} = \frac{1}{a_{ii}} \left( b_i - \sum_{j < i} a_{ij} x_j^{k+1} - \sum_{j > i} a_{ij} x_j^k \right) \tag{10.35}\)$

收敛速度通常比 Jacobi 快 2 倍，但顺序依赖难以并行。

  • SOR 迭代

$\(x_i^{k+1} = (1 - \omega) x_i^k + \frac{\omega}{a_{ii}} \left( b_i - \sum_{j < i} a_{ij} x_j^{k+1} - \sum_{j > i} a_{ij} x_j^k \right) \tag{10.42}\)$

最优松弛因子 \(\omega_{opt} = \frac{2}{1 + \sqrt{1 - \rho^2}}\)，其中 \(\rho\) 是 Jacobi 迭代矩阵的谱半径。

四、Krylov 子空间方法（§10.4）

  • CG 算法（对称正定情形） (Eq. 10.51-10.58)

算法框架： 1. 初始 \(\mathbf{x}^0, \mathbf{r}^0 = \mathbf{b} - \mathbf{A} \mathbf{x}^0, \mathbf{p}^0 = \mathbf{r}^0\) 2. 迭代：\(\alpha_k = \frac{(\mathbf{r}^k)^T \mathbf{r}^k}{(\mathbf{p}^k)^T \mathbf{A} \mathbf{p}^k}\)，\(\mathbf{x}^{k+1} = \mathbf{x}^k + \alpha_k \mathbf{p}^k\) 3. \(\mathbf{r}^{k+1} = \mathbf{r}^k - \alpha_k \mathbf{A} \mathbf{p}^k\)，\(\beta_k = \frac{(\mathbf{r}^{k+1})^T \mathbf{r}^{k+1}}{(\mathbf{r}^k)^T \mathbf{r}^k}\)，\(\mathbf{p}^{k+1} = \mathbf{r}^{k+1} + \beta_k \mathbf{p}^k\) 4. 检查收敛

收敛速度：理论上 \(k\) 步内精确求解（对 \(k \leq N\)），但实际因舍入误差需更多步。条件数 \(\kappa(\mathbf{A})\) 决定收敛速度。

  • 预处理 CG (PCG)

求解等效系统

$\(\mathbf{M}^{-1} \mathbf{A} \mathbf{x} = \mathbf{M}^{-1} \mathbf{b} \tag{10.59}\)$

其中 \(\mathbf{M}\) 是 \(\mathbf{A}\) 的近似逆（\(\mathbf{M}^{-1} \mathbf{A}\) 的条件数远小于 \(\mathbf{A}\)）。

  • GMRES 算法（非对称情形）—— 通过 Arnoldi 过程构造 Krylov 子空间的正交基，在子空间内求解最小残差问题。内存需求随迭代步线性增长（重启 GMRES 可缓解）。

  • BiCGSTAB —— 通过双共轭梯度 + 残差最小化步骤，构造稳定的非 Krylov 子空间迭代。内存需求固定，比 GMRES 节省。

  • ILU(0) / ILU(k) 预处理 (Eq. 10.66) —— 对 \(\mathbf{A}\) 做 LU 分解，丢弃绝对值小于阈值的小元素，保持稀疏性。OpenFOAM 的 DILU、FDIC 即此类。

五、AMG 多重网格（§10.5）

  • 多重网格的核心思想

  • 细网格：消除高频误差

  • 粗网格：消除低频误差

  • 套娃：粗网格之上还有更粗网格

  • V-cycle 流程 (Fig. 10.10)

  • 前光滑（细网格上几次迭代）

  • 计算残差 \(\mathbf{r}\)
  • 限制到粗网格 \(\mathbf{R} \mathbf{r}\)
  • 粗网格求解
  • 延拓回细网格 \(\mathbf{P} \mathbf{e}\)

  • 后光滑

  • OpenFOAM GAMG 实现

solver GAMG;
tolerance 1e-7;
relTol 0;

关键结论

  1. FVM 离散化得到的线性系统是"大型稀疏对角占优"矩阵：求解算法的选择需考虑矩阵规模、稀疏性、对称性、病态程度。
  2. 直接法在工业 CFD 中极少使用：原因是大规模稀疏矩阵的 fill-in 破坏稀疏性，且非线性问题需频繁重分解。TDMA 例外（用于 1D 子问题）。
  3. Krylov 子空间方法是工业 CFD 的标准选择：CG 用于 SPD（压力、温度），GMRES / BiCGSTAB 用于非对称（动量、湍流量）。
  4. 预条件子决定 Krylov 方法的实用价值：没有好的预条件子，Krylov 方法在大规模问题上收敛极慢。ILU、AMG、Block preconditioner 是工业级预条件子。
  5. AMG 是非结构网格 CFD 求解器的"标配"加速器：相比 CG/ILU，AMG 把收敛速度从 \(\mathcal{O}(N)\) 量级降至 \(\mathcal{O}(N^{0.5})\)，对工业规模算例至关重要。
  6. OpenFOAM 通过 fvSolution 配置求解器：用户在 solver、preconditioner、tolerance、relTol 4 个参数中选择，决定整个算例的求解效率。

挑战和开放性问题

  1. 病态系统的收敛性 —— 当 \(\kappa(\mathbf{A})\) 极大（如高 \(Re\) 流动、湍流的各向异性网格），即使有预条件子，Krylov 方法也可能需要数千步才能收敛。Block preconditioner、Physics-based preconditioner 是研究热点。
  2. 多物理场耦合 —— 流-热-化-力多场耦合时，\(\mathbf{A}\) 是 block 结构（每个物理量一个 block），传统预条件子效果差。需用 Schur 补预条件、Physics-based 块预条件。
  3. 并行可扩展性 —— 在 10⁴+ 核上求解时，AMG 的粗网格求解成为瓶颈（仅 1 个核承担）。需要并行粗化、Hypre 的 BoomerAMG 等。
  4. GPU 加速求解器 —— 工业级 CFD 的 GPU 加速是 2016 年后的研究热点，AMGX、cuPDE 等 GPU-native 求解器陆续出现。本书出版年代限制未涉及。
  5. 自适应求解策略 —— 对不同物理场选择不同求解器（如压力用 AMG、动量用 GMRES、湍流量用 PCG），需要 trade-off。

个人反思与批判性分析

本章是全书"数学密度最高"的一章，100+ 页篇幅覆盖了线性求解器的全部主流算法。从写作特点看：

  • 直接法 → 迭代法 → Krylov → AMG 的递进展开 —— 这种结构既照顾了"线性代数基础"的复习，又直接进入"工业 CFD 求解器"的实现细节。读者可以在递进过程中逐步建立对求解器设计的整体认识。
  • CG / GMRES / BiCGSTAB 的算法框架详细给出 —— 这对工业 CFD 求解器开发者是非常有价值的参考。
  • OpenFOAM fvSolution 配置的对应 —— 把抽象的算法映射到具体的 OpenFOAM 配置选项，是工程实用性的体现。

但本章也存在以下不足：

  • Krylov 子空间方法的收敛性理论展开不够 —— CG / GMRES 的收敛速度取决于 \(\kappa(\mathbf{A})\) 和特征值分布，但本书未给出严格的"为什么在病态问题上需要 AMG"的分析。
  • AMG 的具体算法（smoothed aggregation、classical AMG）展开不充分 —— OpenFOAM 默认的 GAMG 是 "Geometric-Algebraic" 混合方法，但本书未详细给出 smoothed aggregation 或 classical AMG 的具体算法步骤。
  • 多物理场耦合的 Block preconditioner —— 这是工业 CFD 的核心难点，本书仅在 §10.6 简略提及，未给出 SIMPLE 算法的压力-速度耦合如何用 Schur 补预条件加速。
  • 并行可扩展性 —— 在 10⁴+ 核上，AMG 的粗网格求解与 Krylov 的全局通信是瓶颈，本书未深入。

总体而言，本章是"线性代数 + 数值分析 + 工业 CFD 求解器"三方面的"百科全书"式参考。对需要进入 OpenFOAM 二次开发（自定义求解器、湍流模型、化学反应模型）的读者是必读章节。

重要参考文献

  • [X1] Saad Y. (2003) Iterative Methods for Sparse Linear Systems, 2nd Edition. SIAM. （Krylov 子空间方法圣经）
  • [X2] Hestenes M.R., Stiefel E. (1952) Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6): 409-436. （CG 原始论文）
  • [X3] Saad Y., Schultz M.H. (1986) GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 7(3): 856-869. （GMRES 原始论文）
  • [X4] van der Vorst H.A. (1992) Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 13(2): 631-644. （BiCGSTAB 原始论文）
  • [X5] Wesseling P. (2004) An Introduction to Multigrid Methods. RT Edwards. （多重网格方法的经典教材）
  • [X6] Stüben K. (2001) A review of algebraic multigrid. Journal of Computational and Applied Mathematics, 128(1-2): 281-309. （AMG 综述）
  • [X7] Briggs W.L., Henson V.E., McCormick S.F. (2000) A Multigrid Tutorial, 2nd Edition. SIAM. （多重网格入门经典）
  • [X8] Trottenberg U., Oosterlee C.W., Schüller A. (2001) Multigrid. Academic Press. （多重网格参考书）
  • [X9] Barrett R., et al. (1994) Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. SIAM. （Krylov 与预条件子的模板实现）
  • [X10] OpenFOAM Foundation (2014) OpenFOAM: The Open Source CFD Toolbox — Programmer's Guide. （OpenFOAM fvSolution 实现）
  • [X11] Trefethen L.N., Bau D. (1997) Numerical Linear Algebra. SIAM. （数值线性代数标准教材）
  • [X12] Demmel J.W. (1997) Applied Numerical Linear Algebra. SIAM. （矩阵计算与稳定性）