工业求解中的方程与矩阵：看懂 CFD 和 FEA 的计算核心 / Equations and Matrices at the Core of CFD and FEA Solvers
📅 创建时间：2026-07-20 🏷️ 标签：#稀疏矩阵 #线性方程 #特征值 #非线性系统 📚 前置知识：[[10-fluid-structure-interaction]]
1. 为什么要先看矩阵性质
同样写成 Ax=b，矩阵可能对称或非对称、正定或不定、稀疏或稠密。矩阵性质决定能否使用某种算法，也决定并行效率。
物理模型 → 离散方法 → 矩阵结构 → 求解算法 → 硬件映射不能先决定“用 GPU”，再倒推应该使用什么算法。
2. 线性方程组
Ax = b在 SAM 中，A 可能是刚度矩阵；在 MarineFlow 中，可能是压力、速度或湍流变量的离散矩阵。
典型计算操作：
- 稀疏矩阵向量乘
y=Ax - 向量点积
- 向量线性组合
- 三角方程求解
- 矩阵分解
- 预条件器应用
3. 稀疏矩阵存储
CSR
保存：
- 非零值数组
- 列索引数组
- 每行起始位置
适合通用稀疏矩阵、逐行遍历和 SpMV，是 CPU/GPU 稀疏计算常用格式。
CSC
按列组织，某些直接分解和列操作更方便。
BSR
按小块存储。结构有限元每个节点包含多个自由度时，矩阵天然具有块结构，BSR 可以减少索引开销并提高局部数据复用。
ELL/SELL
把每行非零项填充到接近统一长度，适合 GPU 规则访问，但行长度差异大时会浪费空间。
4. 对称性
线性弹性结构刚度矩阵通常对称：
K = Kᵀ对称矩阵只需存储一半，并可使用专门算法。
CFD 对流离散、摩擦接触、非保守载荷和某些多物理耦合通常产生非对称矩阵。
5. 正定、半正定与不定
对称正定 SPD
对任意非零向量：
xᵀAx > 0约束充分的线性弹性刚度矩阵通常 SPD，可使用 Cholesky 和 CG。
半正定
未约束结构存在刚体模态，刚度矩阵可能半正定并出现零特征值。
不定矩阵
拉格朗日乘子、鞍点问题、屈曲和某些耦合系统会产生正负特征值并存的矩阵，需要更一般算法。
6. 条件数
κ(A) = ||A|| · ||A⁻¹||条件数大意味着小误差可能被放大，迭代收敛通常较慢。
造成病态的原因：
- 网格尺寸差异过大
- 材料刚度相差巨大
- 过强惩罚参数
- 单元质量差
- 几乎重复的约束
- CFD 网格高度非正交
预条件器的目标可以理解为把原系统变成条件更好的等价系统。
7. CFD 中的矩阵
压力泊松型方程
通常接近椭圆问题，稀疏、邻域耦合，常适合多重网格。
速度和标量输运
对流项使矩阵可能非对称。上风格式更稳定，但会增加数值扩散。
湍流和相分数
同样是输运型系统，但系数非线性强，可能需要限制器和更强稳定化。
CFD 通常不会组装一个包含所有变量的巨大单体矩阵，而是分步求解压力、速度、湍流和相分数。
8. FEA 中的矩阵
静力
Ku=f模态
Kφ=λMφ屈曲
Kφ=-λKgφ隐式动力
时间离散后形成有效刚度：
Keff = K + a0M + a1C非线性
每轮 Newton 迭代形成切线矩阵：
Kt Δu = R9. 特征值问题
特征值问题不直接求全部矩阵元素对应的所有特征值。大型模型通常只求低阶若干模态。
常见计算核：
- 稀疏矩阵向量乘
- 解移位线性系统
- 向量正交化
- Rayleigh 商
当模态数量增加时，正交化可能从次要成本变成瓶颈。
10. 显式动力为什么不解大矩阵
若质量矩阵集中为对角矩阵：
a = M⁻¹(fext-fint)每个自由度可以直接更新，不需全局矩阵分解。这使显式动力具有很好的并行性。
代价是稳定时间步非常小：
Δtcrit ≈ 最小单元特征长度 / 材料波速11. 矩阵性质与算法速查
| 系统性质 | 常用方法 |
|---|---|
| 稠密小矩阵 | LU、QR、SVD |
| 稀疏 SPD | Cholesky、CG、AMG |
| 稀疏对称不定 | LDLᵀ、MINRES |
| 稀疏非对称 | LU、GMRES、BiCGSTAB |
| 多右端项 | 分解复用、Block Krylov |
| 广义特征值 | Lanczos、Arnoldi、LOBPCG |
| 局部显式更新 | SIMD、GPU Kernel |
下一篇：[[12-serial-and-algorithmic-acceleration]]