求解器基础 - CAE 软件的"引擎" / Solvers as the Computational Engine of CAE Software
日期：2026-04-29 标签：#求解器 #线性求解 #非线性 #牛顿法 #并行计算 前置知识：CAE 计算机辅助工程
本章目标
- 理解求解器的工作原理
- 掌握线性方程组求解方法
- 理解非线性问题的求解
- 了解特征值分析方法
- 认识并行计算在求解中的应用
第1部分：求解器概述
1.1 什么是求解器？
求解器 = Solver
CAE 软件中负责核心数值计算的模块，是软件的"引擎"。
┌─────────────────────────────────────────────────────────────────────────────┐
│ CAE 软件架构 │
├─────────────────────────────────────────────────────────────────────────────┤
│ │
│ 前处理 → ┌─────────────────┐ → 后处理 │
│ │ 求解器 │ │
│ │ │ │
│ │ 方程组求解 │ │
│ │ 迭代收敛控制 │ │
│ │ 并行计算管理 │ │
│ └─────────────────┘ │
│ │
│ 求解器 = 数值计算核心，决定分析能力和性能 │
│ │
└─────────────────────────────────────────────────────────────────────────────┘1.2 求解器分类
┌─────────────────────────────────────────────────────────────────────────────┐
│ 求解器分类体系 │
├─────────────────────────────────────────────────────────────────────────────┤
│ │
│ 按问题类型分： │
│ ├── 结构求解器 │
│ │ ├── 静力学求解器 │
│ │ ├── 动力学求解器 │
│ │ └── 非线性求解器 │
│ ├── CFD 求解器 │
│ │ ├── 压力修正法 (SIMPLE) │
│ │ ├── 密度基求解器 │
│ │ └── 涡流求解器 │
│ ├── 热求解器 │
│ └── 电磁求解器 │
│ │
│ 按数学方法分： │
│ ├── 直接求解器 │
│ │ └── 高斯消元、LU分解 │
│ └── 迭代求解器 │
│ ├── 共轭梯度法 (CG) │
│ ├── GMRES │
│ └── 多重网格法 │
│ │
└─────────────────────────────────────────────────────────────────────────────┘第2部分：线性方程组求解
2.1 有限元基本方程
回顾：有限元分析最终归结为求解线性方程组。
[K]{u} = {F}
展开形式：
┌ K₁₁ K₁₂ K₁₃ ┐ ┌ u₁ ┐ ┌ F₁ ┐
│ K₂₁ K₂₂ K₂₃ │ · │ u₂ │ = │ F₂ │
└ K₃₁ K₃₂ K₃₃ ┘ └ u₃ ┘ └ F₃ ┘
刚度矩阵 节点位移 外载荷物理含义：
- K = 刚度矩阵（描述材料对变形的抵抗能力）
- u = 未知位移
- F = 外力
2.2 直接求解法
原理：通过消元直接得到精确解。
高斯消元法：
步骤：把系数矩阵化为上三角，然后回代
示例：
┌ 2 1 │ u₁ ┐ ┌ 8 ┐ ← 方程组
│ 1 3 │ u₂ │ = │ 18 │
└ │ │ └ ┘
步骤1：消元（第二行 - 0.5×第一行）
┌ 2 1 │ u₁ ┐ ┌ 8 ┐
│ 0 2.5 │ u₂ │ = │ 14 │
└ │ │ └ ┘
步骤2：回代
u₂ = 14 / 2.5 = 5.6
u₁ = (8 - 1×5.6) / 2 = 1.2
答案：u₁ = 1.2, u₂ = 5.6LU 分解法：
原理：把 [K] = [L][U]，其中 L 是下三角，U 是上三角
[K] = ┌ 2 1 ┐
└ 1 3 ┘
分解为：
[L] = ┌ 1 ┐ [U] = ┌ 2 1 ┐
└ 0.5 1 ┘ └ 0 2.5 ┘
然后：
1. 求 y: [L]{y} = {F}
2. 求 u: [U]{u} = {y}直接法的特点：
| 优点 | 缺点 |
|---|---|
| 求解精确（数值误差仅来自舍入） | 内存占用大 O(n²) |
| 对病态矩阵稳定 | 计算量 O(n³) |
| 不需要迭代 | 大规模问题效率低 |
2.3 迭代求解法
原理：从一个初始估计开始，逐步逼近真解。
共轭梯度法 (CG)：
核心思想：在每次迭代中，选择一个与之前所有方向共轭的方向
算法流程：
1. 初始化
x₀ = 初始估计（通常为0）
r₀ = b - Ax₀ （残差）
p₀ = r₀ （搜索方向）
2. 迭代（直到收敛）
for k = 0, 1, 2, ...
αₖ = (rₖ · rₖ) / (pₖ · Apₖ) ← 计算步长
xₖ₊₁ = xₖ + αₖ · pₖ ← 更新解
rₖ₊₁ = rₖ - αₖ · Apₖ ← 更新残差
if ||rₖ₊₁|| < ε: 收敛，退出
βₖ = (rₖ₊₁ · rₖ₊₁) / (rₖ · rₖ) ← 计算方向更新因子
pₖ₊₁ = rₖ₊₁ + βₖ · pₖ ← 更新搜索方向收敛判定：
残差范数 ||r|| / ||b|| < 容差（通常 10⁻⁶）
或者：
位移增量 ||Δu|| / ||u|| < 容差迭代法的特点：
| 优点 | 缺点 |
|---|---|
| 内存效率高 | 需要迭代，可能不收敛 |
| 计算量约 O(n·iter) | 收敛速度依赖问题特性 |
| 适合稀疏矩阵 | 需要预条件子加速 |
2.4 预条件子
为什么需要预条件子？
病态矩阵的例子：
┌ 1.00 0.99 ┐
└ 0.99 0.98 ┘
行列式接近0，条件数很大
直接求解也可能不准，迭代法收敛很慢
预条件子 = 变换矩阵，让问题"好解"一些常见预条件子：
1. 对角预条件子（最简单）
M = diag(K)
2. 不完全LU分解 (ILU)
近似 LU 分解，丢弃小元素
3. 代数多重网格 (AMG)
层级粗化，加速收敛第3部分：非线性求解
3.1 什么是非线性？
线性 vs 非线性：
线性问题：
F = K · u
K 是常数，与 u 无关
非线性问题：
F = K(u) · u
K 是 u 的函数，随位移变化非线性来源：
┌─────────────────────────────────────────────────────────────────────────────┐
│ 非线性类型 │
├─────────────────────────────────────────────────────────────────────────────┤
│ │
│ 1. 材料非线性 │
│ 材料应力-应变关系不是直线 │
│ │
│ σ │
│ ↑ │
│ │ ╱ │
│ │ ╱ │
│ │ ╱ │
│ │ ╱ │
│ │ ╱ ← 塑性变形 │
│ │ ╱ │
│ └──────────────────→ ε │
│ │
│ 2. 几何非线性 │
│ 大变形，位移影响几何形状 │
│ │
│ 小变形： │
│ ═════════ │
│ │
│ 大变形： │
│ ═════════════╗ │
│ ╚════ 变形后几何形状完全不同 │
│ │
│ 3. 边界非线性 │
│ 接触、摩擦、螺栓预紧 │
│ │
│ ┌─────┐ │
│ │ │ ← 接触对 │
│ └─────┘ │
│ ↓ │
│ 接触/分离状态变化 │
│ │
└─────────────────────────────────────────────────────────────────────────────┘3.2 牛顿-拉夫森法
原理：通过切线刚度迭代逼近非线性解。
迭代公式：
{K_T(uₖ)} · Δuₖ = F - K(uₖ) · uₖ
↑ ↑
切线刚度 残差力
uₖ₊₁ = uₖ + Δuₖ算法流程：
┌─────────────────────────────────────────────────────────────────────────────┐
│ 牛顿-拉夫森迭代流程 │
├─────────────────────────────────────────────────────────────────────────────┤
│ │
│ 1. 初始化 │
│ u₀ = 0 │
│ 施加总载荷 F │
│ │
│ 2. 增量加载 │
│ for 载荷步 i = 1 to N: │
│ │
│ ΔF = F · i/N - F · (i-1)/N │
│ │
│ for 迭代 j = 1 to max_iter: │
│ │ │
│ ├── 计算 K_T(uⱼ) ← 切线刚度 │
│ ├── 计算 R = F - K(uⱼ)·uⱼ ← 残差力 │
│ ├── if ||R|| < ε: 收敛，退出内循环 │
│ ├── 求解 Δu = K_T⁻¹ · R │
│ └── uⱼ₊₁ = uⱼ + Δu │
│ │
│ 3. 收敛后保存 u │
│ │
│ 4. 输出结果 │
│ │
└─────────────────────────────────────────────────────────────────────────────┘收敛判定：
| 准则 | 条件 | 说明 |
|---|---|---|
| 力残差 | ||R|| / ||F|| < 0.001 | 残余力很小 |
| 位移增量 | ||Δu|| / ||u|| < 0.001 | 位移变化很小 |
| 应变能 | ΔE / E < 0.001 | 能量变化很小 |
图示：
F
↑
│ ╱│
│ ╱ │
│ ╱ │
│ ╱ │
│╱ │ ← 实际曲线
│ ╲
│ ╲
│ ╲
│ ╲ 切线
│ ╲╱
│ ╲
│ ╲___
│ ‾‾‾‾‾‾‾‾ 近似曲线
└────────────────────────→ u
每次迭代：计算切线 → 求解增量 → 更新位移
迭代直到近似曲线与实际曲线足够接近3.3 收敛困难的原因
常见不收敛情况：
1. 载荷步太大
跳过了关键点（如载荷峰值）
2. 材料问题
塑性流动、软化导致不收敛
3. 几何不稳定
屈曲、 snap-through
4. 接触问题
剧烈接触/分离
5. 模型错误
约束不足、载荷施加错误调试策略：
┌─────────────────────────────────────────────────────────────────────────────┐
│ 收敛调试策略 │
├─────────────────────────────────────────────────────────────────────────────┤
│ │
│ 策略1：减小载荷步长 │
│ 把 10 步改为 50 步，缓慢加载 │
│ │
│ 策略2：开启线搜索 │
│ 自动调整迭代步长 │
│ │
│ 策略3：使用弧长法 │
│ 追踪复杂路径（如屈曲） │
│ │
│ 策略4：检查模型 │
│ 约束是否充分、载荷是否合理 │
│ │
│ 策略5：简化问题 │
│ 先求解线性问题，确认基本正确 │
│ │
└─────────────────────────────────────────────────────────────────────────────┘第4部分：特征值分析
4.1 特征值问题
方程：
[K]{φ} = λ[M]{φ}物理含义：
- K = 刚度矩阵
- M = 质量矩阵
- λ = 特征值（ω² = 固有频率的平方）
- {φ} = 振型向量（固有振型）
应用：
- 模态分析（固有频率、振型）
- 屈曲分析（临界载荷）
- 动力响应分析
4.2 求解方法
Lanczos 方法（最常用）：
特点：
- 专门用于稀疏矩阵
- 可以只求前几阶模态
- 效率高
算法：
1. 生成 Lanczos 向量
2. 在降维空间中求解
3. 反变换得到原空间振型子空间迭代法：
特点：
- 同时迭代多个向量
- 适合求多阶模态
算法：
1. 初始化向量组
2. 交替迭代：
- 用 K⁻¹ 乘（逆迭代）
- 用 M 正交化（质量归一）
3. 收敛后得到特征对4.3 模态分析示例
问题：求解悬臂梁的前三阶固有频率
已知：
长度 L = 1m
截面 0.1m × 0.2m
材料 钢 (E = 200 GPa, ρ = 7850 kg/m³)
理论解：
f₁ = 1.875² · EI / (2π · L² · ρA) ≈ 12 Hz
f₂ = 4.694² · EI / (2π · L² · ρA) ≈ 75 Hz
f₃ = 7.855² · EI / (2π · L² · ρA) ≈ 210 Hz
数值解（FEA）：
f₁ = 11.8 Hz
f₂ = 74.2 Hz
f₃ = 207.5 Hz
误差 < 2%，验证正确第5部分：并行计算
5.1 为什么需要并行？
问题复杂度增长：
100万自由度 → 10秒
1000万自由度 → 1000秒 ≈ 17分钟（单机困难）
1亿自由度 → 100000秒 ≈ 27小时（必须并行）并行加速比：
┌─────────────────────────────────────────────────────────────────────────────┐
│ 并行加速比 │
├─────────────────────────────────────────────────────────────────────────────┤
│ │
│ 理想情况： │
│ S(p) = p │
│ │
│ 实际情况： │
│ S(p) = p / (1 + σ(p-1)) ← Amdahl定律 │
│ │
│ 其中 σ 是不可并行部分的比例 │
│ │
│ σ = 5%: 8核加速 ≈ 5.7倍 │
│ σ = 10%: 8核加速 ≈ 4.7倍 │
│ σ = 20%: 8核加速 ≈ 3.5倍 │
│ │
└─────────────────────────────────────────────────────────────────────────────┘5.2 并行计算框架
MPI（Message Passing Interface）：
多进程，每个进程独立内存
进程0 ──── 通信 ──── 进程1
│ │
│ ──── 通信 ──── │
│ │
进程2 ──── 通信 ──── 进程3
优点：扩展性好，适合集群
缺点：编程复杂OpenMP：
多线程，共享内存
线程0 ──┐
线程1 ──┼── 共享内存
线程2 ──┤
线程3 ──┘
优点：编程简单，共享数据方便
缺点：内存共享限制扩展性CUDA / GPU加速：
CPU：少量强大核心
GPU：大量简单核心（数千个）
适合：大规模矩阵运算、蒙特卡洛、神经网络
求解器中的 GPU 加速：
- 矩阵-向量乘
- 稀疏矩阵运算
- 预条件子应用5.3 分布式求解
域分解方法：
┌─────────────────────────────────────────────────────────────────────────────┐
│ 域分解并行 │
├─────────────────────────────────────────────────────────────────────────────┤
│ │
│ 把网格分成多个区域，每个区域分配给一个处理器 │
│ │
│ ┌────────┬────────┐ │
│ │ 区域1 │ 区域2 │ ← 处理器1 │
│ │ CPU 1 │ CPU 2 │ │
│ ├────────┼────────┤ │
│ │ 区域3 │ 区域4 │ ← 处理器3 │
│ │ CPU 3 │ CPU 4 │ │
│ └────────┴────────┘ │
│ ↑ │
│ │ │
│ 边界信息交换 │
│ │
└─────────────────────────────────────────────────────────────────────────────┘第6部分：主流求解器介绍
6.1 结构求解器
| 求解器 | 公司 | 特点 |
|---|---|---|
| Ansys Mechanical | Ansys | 通用、结构完整 |
| Abaqus Standard | SIMULIA | 隐式、非线性强 |
| Abaqus Explicit | SIMULIA | 显式、接触、高速 |
| NX Nastran | Siemens | 航空航天、背书 |
6.2 CFD 求解器
| 求解器 | 公司 | 特点 |
|---|---|---|
| ANSYS Fluent | Ansys | 工业主流 |
| ANSYS CFX | Ansys | 耦合强 |
| STAR-CCM+ | Siemens | 自动网格 |
| OpenFOAM | 开源 | 免费、定制 |
6.3 求解器技术对比
┌─────────────────────────────────────────────────────────────────────────────┐
│ 隐式 vs 显式求解 │
├─────────────────────────────────────────────────────────────────────────────┤
│ │
│ 隐式求解 (Implicit)： │
│ ├─ 特点：K(u)·Δu = R │
│ ├─ 优点：时间步长大、静态问题高效 │
│ ├─ 缺点：每次迭代要解方程组 │
│ └─ 代表：Abaqus Standard、Ansys Mechanical │
│ │
│ 显式求解 (Explicit)： │
│ ├─ 特点：M·ü = F - K·u │
│ ├─ 优点：无须求解方程组、接触容易处理 │
│ ├─ 缺点：时间步长受限制、动态问题专用 │
│ └─ 代表：Abaqus Explicit、LS-DYNA │
│ │
│ 选择建议： │
│ - 静态/准静态 → 隐式 │
│ - 高速碰撞/成形 → 显式 │
│ - 复杂接触 → 看软件擅长 │
│ │
└─────────────────────────────────────────────────────────────────────────────┘核心总结
总结1：求解器类型
线性求解：
直接法：高斯消元、LU分解（精确、内存大）
迭代法：CG、GMRES（省内存、需收敛）
非线性求解：
牛顿-拉夫森：切线刚度迭代
弧长法：处理屈曲等复杂路径
特征值求解：
Lanczos：稀疏矩阵、求低阶模态
子空间迭代：求多阶模态总结2：收敛控制
收敛准则：
- 力残差 < 容差
- 位移增量 < 容差
- 能量变化 < 容差
加速收敛：
- 预条件子
- 自动时间步长
- 弧长法总结3：并行策略
单机多核：OpenMP
多机集群：MPI
GPU加速：CUDA
域分解：每个区域独立求解章节测试
测试1：基本方程
有限元分析的三大基本方程之一是什么？各符号含义？
测试2：直接法 vs 迭代法
直接求解法和迭代求解法各有什么优缺点？
测试3：牛顿法
牛顿-拉夫森法中，"切线刚度"是什么？为什么要用切线刚度？
测试4：收敛判定
非线性分析中，如何判断迭代是否收敛？
测试5：显式 vs 隐式
什么情况下应该选择显式求解器而不是隐式求解器？
参考答案
测试1答案
基本方程 [K]{u} = {F}，其中 [K] 是刚度矩阵，{u} 是节点位移向量，{F} 是外载荷向量。
测试2答案
- 直接法：精确、稳定，但内存占用大、计算量 O(n³)；适合中小规模问题
- 迭代法：内存效率高、计算量约 O(n·iter)，但需要迭代、可能不收敛；适合大规模问题
测试3答案
切线刚度是当前位移状态下的瞬时刚度矩阵。非线性问题中，刚度随位移变化，用切线刚度能更准确地描述当前状态的力学响应。
测试4答案
通过力残差、位移增量或能量变化与容差比较。当这些值小于设定的容差（如 0.001）时，认为迭代收敛。
测试5答案
应该选择显式求解器的情况：
- 高速碰撞、爆炸、冲击问题
- 大规模接触问题
- 成形加工（如冲压）
- 瞬态动力学问题
原因：显式求解器不需要在每个时间步求解线性方程组，对接触问题处理更稳定，虽然时间步长受限但这类问题天然需要小步长。
相关笔记
- [[03 - CAE 计算机辅助工程]] - 仿真基础
- [[04 - 前处理与后处理]] - 前后处理
- [[08 - 技术栈与开发语言]] - 求解器开发技术
下一步学习
- [ ] 阅读 06 - 主流工业软件介绍
学习状态：完成