结构求解器设计：从单元积分到非线性、动力与接触 / Structural Solver Design from Element Integration to Contact and Dynamics

📅 创建时间：2026-07-23
🏷️ 标签：#FEA #结构求解器 #矩阵组装 #非线性 #接触
📚 前置知识：[[01-production-solver-architecture]] [[../marine-cae/07-fea-from-physics-to-matrix]]
📚 相关知识：[[../marine-cae/09-analysis-types]] [[../marine-cae/11-discrete-systems-and-matrices]]

---

1. 结构求解器究竟在做什么

有限元把连续结构分成单元，用有限个自由度近似位移场。生产程序需要重复完成：

网格与属性
→ 自由度编号
→ 单元积分
→ 局部矩阵和向量
→ 全局组装
→ 约束与载荷
→ 线性/非线性/特征值求解
→ 位移、应力和工程量恢复

不同结构分析共用前半段，但外层方程、状态历史和瓶颈并不相同。

2. 自由度编号

典型自由度：

• 实体单元：节点平移；
• 梁、壳单元：平移与转动；
• 热-结构耦合：位移加温度；
• 混合或不可压缩单元：可能增加压力等内部变量；
• 接触：可能引入乘子或约束变量。

DofManager 需要输出每个单元的方程号数组：

element 42 nodes = [8, 11, 19, 20]
element dofs = [
  ux8, uy8, uz8,
  ux11, uy11, uz11,
  ux19, uy19, uz19,
  ux20, uy20, uz20
]

编号还影响：

• 稀疏矩阵非零分布；
• 直接法 Fill-in；
• SpMV 的向量访问局部性；
• MPI 子域边界大小。

因此常使用 RCM、Nested Dissection 或图分区后的局部重编号。

3. 单元 Kernel

以静力实体单元为例：

Ke = ∫ Bᵀ D B dΩ
fe = ∫ Nᵀ b dΩ + ∫ Nᵀ t dΓ

程序对每个积分点执行：

1. 读取节点坐标和当前状态；
2. 计算形函数及自然坐标导数；
3. 计算 Jacobian、行列式和逆；
4. 得到物理坐标中的梯度与 B 矩阵；
5. 调用材料模型得到应力和切线；
6. 累加局部残差与切线矩阵；
7. 输出需要保存的积分点状态。

单元 Kernel 通常计算密集，但全局组装会产生不规则写入。优化时要分别计时，不能把两者混为“组装”。

4. 稀疏模式预生成

结构网格固定时，矩阵非零结构通常在求解前确定：

遍历单元
→ 收集单元自由度之间的邻接
→ 去重并排序
→ 生成 CSR rowPtr / colIdx
→ 建立局部条目到全局 value 位置的映射

预先建立 Scatter Map，可以避免每次组装都在 CSR 行中搜索列号。代价是额外索引内存，需要在内存占用与组装速度间权衡。

接触激活、单元删除或自适应网格会改变拓扑，此时需要动态模式、预留空间或阶段性重建。

5. 全局组装策略

多个单元共享节点，并行组装存在写冲突。

方法	优点	局限
原子加	实现直接，GPU 可用	冲突高时吞吐下降
图着色	同色单元无写冲突	颜色间同步，预处理和重排成本
线程局部矩阵	无共享写冲突	内存开销大，最终归并昂贵
COO 生成后排序归并	高度并行，适合 GPU	临时数据量大
子域私有组装	适合 MPI 域分解	界面自由度需要通信与累加
Matrix-Free	避免存储和组装全局矩阵	预条件、边界和复杂材料实现更难

选择依据是单元类型、平均节点度、矩阵块结构、硬件和后续求解算法。

6. 约束怎样进入系统

6.1 消元

固定自由度从系统中移除或将对应行列修改为单位形式。优点是系统小且稳定；缺点是多点约束实现更复杂。

6.2 罚函数

向约束方向加入很大的刚度。实现简单，但过大的罚参数会恶化条件数，过小则产生明显约束误差。

6.3 拉格朗日乘子

增加未知量：

[ K  Cᵀ ] [u] = [f]
[ C   0 ] [λ]   [g]

约束可精确满足，但矩阵通常变为对称不定，不能继续假设适合普通 CG。

6.4 多点约束与静力凝聚

从属自由度可以写成主自由度线性组合；单元内部自由度可以通过 Schur Complement 消去。凝聚减少全局未知量，但增加局部计算。

7. 线性静力

K u = f

在线弹性、约束充分时，K 通常对称正定：

• 中小规模、多载荷且内存足够：稀疏 Cholesky；
• 大规模：CG + AMG/IC/Domain Decomposition；
• 多右端项：优先复用符号分解、数值因子或预条件器。

主要成本可能是单元组装、矩阵分解、SpMV、预条件器或应力恢复，不能只看求解器名称判断。

8. 非线性静力

目标是恢复内外力平衡：

r(u) = f_ext - f_int(u) = 0
K_t(u_k) Δu = r(u_k)
u_{k+1} = u_k + α Δu

一次增量通常包含：

1. 从已提交状态开始；
2. 计算内力和一致切线；
3. 施加约束；
4. 求解位移增量；
5. 线搜索或阻尼；
6. 更新试算状态；
7. 检查力、位移和能量准则；
8. 收敛后提交材料与接触历史。

性能必须拆成：

总时间 =
增量数
× 每增量平均 Newton 次数
× 每次 Newton 的组装与线性求解成本

如果 Newton 次数异常高，应先检查切线一致性、载荷步、接触和模型病态，而不是先优化 SpMV。

9. 模态与屈曲

模态

K φ = λ M φ
λ = ω²

通常只求低阶特征对，使用 Lanczos、Arnoldi 或 LOBPCG。Shift-and-Invert 会反复求解偏移线性系统，因子复用和内存是关键。

线性屈曲

(K + λ K_g) φ = 0

需要先完成预应力静力分析，再建立几何刚度。屈曲系数是理想化指标，不能替代含缺陷的非线性失稳分析。

10. 隐式动力

M ü + C u̇ + K u = f(t)

Newmark 或 Generalized-α 把每个时间步转化为有效刚度系统：

K_eff = K_t + a0 M + a1 C

每步可能需要 Newton 迭代。时间步由精度和高频控制决定，不只是稳定性。优化应检查：

• K、M、C 是否可复用；
• 有效矩阵结构是否固定；
• 预条件器多久重建一次；
• 输出频率是否远高于工程需求。

11. 显式动力

使用对角集中质量时：

a_n = M⁻¹(f_ext - f_int)

无需解大型线性方程，但稳定时间步受最小单元控制。主要热点是：

• 单元内力；
• 接触搜索与接触力；
• 状态更新；
• 高频结果输出。

畸变小单元可能让时间步极小。质量缩放能提高速度，但必须评估惯性与能量比例，不能只报告加速比。

12. 接触

接触通常包含两类问题：

1. 几何搜索：可能接触的面或点在哪里；
2. 约束执行：怎样阻止穿透并处理摩擦。

搜索使用包围盒、BVH、空间哈希等；约束可用罚函数、乘子或增广拉格朗日。接触集合变化会引起刚度突变、负载不均和非线性困难。

诊断应记录：

• 候选对与活动接触对数量；
• 搜索、局部接触计算和线性求解时间；
• 最大穿透、接触力和平衡；
• 接触区域在 Rank 间的分布。

13. 应力恢复和后处理

位移求解完成后仍需：

• 积分点应变与应力；
• 节点外推与平均；
• 主应力、von Mises 和安全系数；
• 多工况包络；
• 结果压缩和写出。

这些操作数据并行度高，适合 SIMD、OpenMP 和 GPU，但常受内存带宽和输出限制。[[/06-other/01-projects/profile/C++/2-StructKernelBench/2-StructKernelBench]] 提供了相关实验路线。

14. 观察、解释与下一步

观察到什么	说明什么	下一步检查什么
矩阵奇异	可能存在刚体模态、失效单元或不定约束系统	约束充分性、零主元位置和乘子方程
直接法内存爆炸	Fill-in 而非原始非零元主导内存	排序、消元树、峰值预测和块结构
CG 不收敛	系统可能非 SPD，或条件数与缩放很差	矩阵对称性、特征性质和预条件器
Newton 迭代多	切线、增量或接触活动集可能不稳定	切线一致性、线搜索和载荷步历史
组装扩展差	写冲突、负载不均或分配开销占主导	原子冲突、线程时间和临时内存
显式计算极慢	最小单元、接触或输出控制总成本	稳定时间步分布、接触占比和写出频率
GPU 加速不明显	模型太小、数据往返或带宽算子占主导	时间线、传输字节和 Roofline

15. 正确性检查

• 静力：外力与反力平衡；
• 能量：外功、内能、动能和接触能关系合理；
• 网格：关键位移和应力随细化趋于稳定；
• 动力：时间步收敛，频率和相位正确；
• 模态：质量归一化、模态正交和刚体模态数量正确；
• 非线性：路径、极限点和材料历史合理；
• 接触：穿透在允许范围，接触力与载荷路径一致；
• 并行：不同线程或 Rank 数下工程量在容差内一致。

下一篇：[[03-cfd-solver-design]]