数值分析基础 — 有限元、有限体积、有限差分理论与求解器选择

1. 偏微分方程分类与数值方法综述

1.1 二阶PDE的分类

一般形式的二阶线性PDE（二维，$u = u(x,y)$）：$A u_{xx} + B u_{xy} + C u_{yy} + \ldots = 0$，判别式 $\Delta = B^2 - 4AC$ 决定类型：

PDE类型直接影响数值方法选择：椭圆型适合直接求解器；抛物型适合时间推进（隐式稳定）；双曲型适合显式格式（但需稳定性条件）。

1.2 三大数值方法对比

2. 有限差分法（FDM）

2.1 常用差分格式

泰勒展开推导差分公式：

一阶迎风格式（对流主导问题稳定）：

2.2 稳定性——冯诺依曼分析

显式格式的稳定性条件（以热扩散方程为例）：

违反此条件会导致数值振荡并发散。隐式格式（如Crank-Nicolson）无条件稳定，但每时间步需求解线性方程组。

2.3 FDM的局限性

• 复杂几何难以处理（需要贴体坐标变换，增加实现难度）
• 非结构网格支持差（规则/结构化网格为主）
• 边界条件施加不如FEM/FVM灵活
• 但在均匀矩形网格+简单边界情况下，FDM最高效、代码最简洁

3. 有限元法（FEM）数学基础

3.1 从强形式到弱形式

以泊松方程为例（强形式）：$-\nabla^2 u = f$ 在 $\Omega$ 上，$u = 0$ 在 $\partial\Omega$ 上。

乘以任意测试函数 $v$（满足 $v = 0$ 在 $\partial\Omega$）并在 $\Omega$ 上积分：

通过分部积分（Green定理）降低对 $u$ 的可微性要求：

弱形式只要求 $u$ 有一阶偏导，而强形式要求二阶偏导存在。

3.2 Galerkin离散——形函数

将解近似为有限维子空间中的函数：$u_h = \sum_{j=1}^n u_j N_j(x)$，用同样的函数作为测试函数（Galerkin法），代入弱形式，得到线性方程组 $\mathbf{K}\mathbf{u} = \mathbf{f}$：

3.3 高阶单元与收敛率

$p$ 阶完全多项式单元对光滑解的收敛率：

$h$ 是网格尺寸，$p$ 是多项式阶次，$C$ 是与网格无关的常数。对于足够光滑的解，提高 $p$（h-p 加密）比单纯加密网格（h 加密）更高效。这就是谱元法、p型FEM的出发点。

3.4 常用单元及其特性

4. 有限体积法（FVM）

4.1 FVM的守恒形式

对守恒律 $\frac{\partial \phi}{\partial t} + \nabla \cdot \mathbf{F}(\phi) = S$ 在控制体 $V_P$ 上积分，并用高斯散度定理：

对所有面上的通量求和（以一维为例，$P$ 是当前格，$E/W$ 是东/西邻格）：

$F_e$、$F_w$ 是界面通量，需要插值计算（此处是各种格式——迎风、中心、QUICK的区别所在）。

4.2 通量插值格式对比

4.3 OpenFOAM的FVM实现

OpenFOAM是开源FVM求解器，其代码结构高度模块化：

// OpenFOAM的扩散+对流方程（C++对象语法）
solve(
    fvm::ddt(rho, phi)           // 时间导数
  + fvm::div(phi, U)             // 对流（使用插值）
  - fvm::laplacian(nu, U)        // 扩散
  == fvOptions(rho, phi)         // 源项
);

这种高级语法直接对应PDE，极大降低了实现新求解器的门槛。

5. 线性方程组求解器

5.1 直接求解器

基于LU分解：$\mathbf{A} = \mathbf{L}\mathbf{U}$（$\mathbf{L}$ 下三角，$\mathbf{U}$ 上三角），求解 $\mathbf{A}\mathbf{x} = \mathbf{b}$ 分为：

1. 分解阶段：$\mathbf{A} = \mathbf{L}\mathbf{U}$，计算量 $O(n^3)$（一般矩阵）或 $O(n \cdot b^2)$（带宽 $b$ 的带状矩阵）
2. 前代/回代：$\mathbf{L}\mathbf{y} = \mathbf{b}$，$\mathbf{U}\mathbf{x} = \mathbf{y}$，计算量 $O(n^2)$

稀疏直接求解器（MUMPS, PARDISO, SuperLU）通过填充最小化重排序（如Metis, AMD算法）大幅减少非零元素填入，使实际计算量远低于密集矩阵的 $O(n^3)$。3D结构问题约为 $O(n^{2})$。

5.2 迭代求解器

求解 $\mathbf{A}\mathbf{x} = \mathbf{b}$ 的迭代法：

• CG（共轭梯度）：仅适用于对称正定矩阵（SPD），结构静力学的刚度矩阵满足此条件；收敛性好，每步计算一次矩阵-向量乘法
• GMRES（广义最小残差）：适用于非对称矩阵（CFD压力方程、非对称刚度矩阵）；需要存储Krylov子空间向量，内存随迭代步增加（重启策略：GMRES(m)）
• BiCGSTAB：适用于非对称矩阵，内存需求固定（不需要重启），CFD中广泛使用

5.3 预处理器——提速的关键

迭代求解的收敛速度取决于矩阵条件数 $\kappa(\mathbf{A})$。预处理器 $\mathbf{M} \approx \mathbf{A}^{-1}$ 将问题变换为 $\mathbf{M}^{-1}\mathbf{A}\mathbf{x} = \mathbf{M}^{-1}\mathbf{b}$，使条件数减小：

6. 特征值求解器

6.1 常用特征值求解算法

7. 数值精度与误差管理

7.1 Richardson外推

设关键量在三套网格上的值为 $f_1$（细）、$f_2$（中）、$f_3$（粗），细化比 $r = h_2/h_1 = h_3/h_2$，数值精度阶次 $p$：

外推"精确"值：

7.2 网格收敛指数（GCI）

细网格上的GCI（相对误差上界）：

$e_{12} = (f_2 - f_1)/f_1$ 是相对误差，$F_s = 1.25$ 是安全系数（标准建议值）。一般要求 $GCI < 3\%$ 可认为网格收敛。

7.3 错误的来源与分类

8. 求解器选择指南