数值分析的基础 — FEMFVMFDM的理論精度

PDE与数值解法的全景

CAE的所有分析工作，归根结底是用数值方法求解偏微分方程（PDE）。为了在计算机上处理PDE，需要将连续体离散化为有限个点、单元或网格。其中FDM、FEM、FVM三种方法是主要的求解方式。

1.1 PDE的分类

二阶线性PDE根据系数的性质可分为3种：

1.2 各方法的哲学差异

有限差分法（FDM）

2.1 泰勒展开与差分近似

泰勒展开式：

由此推导一阶导数的近似：

二阶导数的二阶中心差分：

2.2 显式法的稳定性条件 — CFL条件与冯·诺依曼稳定性分析

一维扩散方程 $\partial u/\partial t = \alpha \partial^2 u/\partial x^2$ 用FTCS（前向时间·中心空间）格式离散化：

冯·诺依曼稳定性分析：将 $u_i^n = G^n e^{ikx_i}$ 代入得放大因子 $G$：

稳定性条件 $|G| \leq 1$ 给出：

平流方程 $\partial u/\partial t + c\partial u/\partial x = 0$ 的CFL（Courant-Friedrichs-Lewy）条件：

物理含义：一个时间步内平流距离不能超过网格宽度。

2.3 高阶差分方案对比

有限元法（FEM）的数学基础

3.1 加权残差法与伽辽金法

椭圆型PDE的例子：$-\nabla \cdot (k\nabla u) = f$ in $\Omega$，边界条件: $u = g$ on $\Gamma_D$。

强形式: $r(u) = -\nabla \cdot (k\nabla u) - f = 0$ 在所有点都成立。

加权残差法： 对于权函数 $v$：

弱形式（用格林定理分部积分）：

伽辽金法： 从同一空间选择近似解和权函数 — $u \approx \sum_j u_j N_j(\mathbf{x})$, $v = N_i(\mathbf{x})$：

这以矩阵形式表示为 $\mathbf{K}\mathbf{u} = \mathbf{f}$。

3.2 形状函数的种类

代表性的形状函数：

• 拉格朗日插值基函数： $N_i$ 在节点 $j \neq i$ 处为零，在节点 $i$ 处为1。$C^0$ 连续性。结构·热分析的标准选择。
• 埃尔米特（Hermite）基函数： 节点值包含变位和梯度。$C^1$ 连续性。梁·板分析中使用。
• 塞伦迪皮蒂单元： 具有边中点节点的高阶拉格朗日单元（8节点四边形等）。

一维拉格朗日形状函数（2节点线性单元，$\xi \in [-1,1]$）：

3节点二阶单元（$\xi \in [-1,1]$）：

3.3 等参数单元与数值积分

物理坐标 $\mathbf{x}$ 与参考坐标 $\boldsymbol{\xi}$ 的变换（等参数变换）：

积分的变换：

高斯-勒让德求积（$n_g$ 点）：

可精确积分至 $2n_g - 1$ 次多项式（$n_g$ 点高斯积分的精度）。

有限体积法（FVM）

4.1 守恒律的积分形式

一般性保存方程（平流-扩散）的积分形式（应用于每个控制体积CV）：

$\phi$：通用变量（密度·动量·能量等），$\Gamma$：扩散系数，$S_\phi$：源项。

物理意义：「CV中积累的量 = 流入平流通量 - 流出平流通量 + 扩散通量 + 体积源」。

4.2 单元中心 vs 单元顶点

4.3 界面通量的插值方案

界面 $f$ 处标量值 $\phi_f$ 的计算方案比较：

$F_f = \mathbf{v}_f \cdot \mathbf{S}_f$ 是界面通量。

高分辨率方案（TVD: Total Variation Diminishing）用通量限制器 $\psi(r)$ 兼顾精度与有界性：

线性方程组求解器

FEM/FVM最终需要求解大规模稀疏矩阵方程 $\mathbf{A}\mathbf{x} = \mathbf{b}$。根据问题规模和性质选择求解器。

5.1 直接法

LU分解： $\mathbf{A} = \mathbf{L}\mathbf{U}$（$\mathbf{L}$：下三角，$\mathbf{U}$：上三角）。通过前代和回代求解。

Cholesky分解（对称正定矩阵用）： $\mathbf{A} = \mathbf{L}\mathbf{L}^T$。计算量约为LU分解的一半。结构分析静力学中常用。

直接法的复杂度（带状矩阵，半带宽 $b$）：

• 内存：$O(n \cdot b)$
• 计算量：$O(n \cdot b^2)$
• 三维问题中半带宽 $b \sim \sqrt{n}$，故计算量 $O(n^2)$ — 大规模问题需迭代法

5.2 迭代法 — 共轭梯度法（CG法）

对称正定矩阵 $\mathbf{A}$ 的CG法算法：

初值: x₀（通常 x₀ = 0）
残差: r₀ = b - Ax₀
搜索方向: p₀ = r₀

for k = 0, 1, 2, ...
  α_k = (r_k^T r_k) / (p_k^T A p_k)   // 步长
  x_{k+1} = x_k + α_k p_k              // 解的更新
  r_{k+1} = r_k - α_k A p_k            // 残差的更新
  β_k = (r_{k+1}^T r_{k+1}) / (r_k^T r_k)
  p_{k+1} = r_{k+1} + β_k p_k          // 搜索方向更新

收敛条件：$\|r_k\| / \|b\| < \varepsilon$（通常 $\varepsilon = 10^{-6}$ 〜 $10^{-8}$）

5.3 预处理（Preconditioning）

条件数 $\kappa(\mathbf{A})$ 较大（病态）的矩阵收敛缓慢。用预处理矩阵 $\mathbf{M} \approx \mathbf{A}$ 求解等价系统 $\mathbf{M}^{-1}\mathbf{A}\mathbf{x} = \mathbf{M}^{-1}\mathbf{b}$。

5.4 稀疏矩阵的存储格式

高效存储大规模稀疏矩阵的数据结构：

• CSR（Compressed Sparse Row）： 按行存储非零元素的列指标和值。行运算（矩阵-向量乘积）高速。FEM·FVM中标准。
• CSC（Compressed Sparse Column）： 按列存储。LU分解合适。
• COO（Coordinate List）： (行, 列, 值)的列表。组装容易但运算慢。

$n$ 节点、$nnz$ 非零元素的CSR矩阵内存：$3 \times nnz + (n+1)$ 整数 + $nnz$ 浮点数。

特征值求解器

6.1 幂法与反迭代法

幂法（Power method）： 求最大特征值的最简单方法：

收敛速率：$|\lambda_2/\lambda_1|$（第二大特征值与最大特征值之比）。特征值密集时收敛非常慢。

反迭代法（Inverse iteration）： 用位移 $\mu$ 的反迭代求与指定值 $\lambda \approx \mu$ 最接近的模式：

6.2 Lanczos法

大规模对称特征值问题中应用最广的方法。对Krylov子空间 $\mathcal{K}_m(\mathbf{A}, \mathbf{v}_1) = \text{span}\{\mathbf{v}_1, \mathbf{A}\mathbf{v}_1, \ldots, \mathbf{A}^{m-1}\mathbf{v}_1\}$ 进行三对角化：

$\mathbf{T}_m$ 是 $m \times m$ 三对角矩阵，其特征值（Ritz值）是 $\mathbf{A}$ 特征值的良好近似。即使 $n = 10^6$ 自由度的系统，$m = 100 \sim 300$ 个Lanczos向量也足以达到足够精度。

6.3 子空间反迭代法与ARPACK

实用的大规模特征值求解器多数基于ARPACK（Arnoldi PACKage），实现了对非对称矩阵也适用的隐式重启Arnoldi法（IRAM）。

基于ARPACK的特征值求解器例子：

• SLEPc（OpenFOAM, FEniCS等）
• LAPACK/ScaLAPACK（中小规模）
• FEAST特征值求解器（特定频带特征值的高并行计算）

数值精度与误差管理

7.1 误差的种类

7.2 Richardson外推与GCI

网格收敛研究中用3个网格解估计真解的Richardson外推法（Richardson Extrapolation）：

从两个网格（$h_1 < h_2$，细化比 $r = h_2/h_1$）得到：

网格收敛指数（GCI: Grid Convergence Index）： Roache（1994）提出的标准误差评价指标：

$e_{12} = (f_2 - f_1)/f_1$（相对误差），$F_s = 1.25$（安全系数，Roache推荐）。

GCI < 5%时"收敛充分"被普遍认可。ASME V&V 20-2009标准也采纳了这一指标。

7.3 数值陷阱与对策

最优求解方法选择指南

8.1 问题类型别求解器选择表

8.2 并行计算可扩展性

阿姆达尔定律： 可并行化比例 $p$ 时，$N$ 处理器的理论最大加速比：

$p = 0.95$（95%可并行化）时，理论极限 $N \to \infty$ 仅 $1/(1-0.95) = 20$ 倍。

大规模FEM/FVM并行计算的注意事项：

• 域分割（Domain Decomposition）： MeTIS、SCOTCH等图分割。最小化halo交换通信成本。
• AMG的并行化： 粗化阶段需进程间通信。可扩展性有限。
• 实际并行效率： 实务上50〜70%。通信支配时更低。

相关交互式工具