非线性有限元求解器间对比 — 接触、稳定化、路径依赖性的V&V总结
为什么非线性解析中求解器间差异会扩大
老师，非线性解析时求解器间的差异会更大吗？线性解析时感觉"都一样"…
差异确实更大。可以说是数量级上的差异。线性静解析中理论解的误差在0.01%以内，而接触问题中5~10%的差异很常见。
即使是同一模型、同一网格也会这样吗？
是的。原因有很多，但主要可以分为三类。
1. 接触算法的差异 — 各求解器默认采用的方法不同，如惩罚法、增广拉格朗日法、Mortar法等。
2. 收敛判定基准的差异 — 力的残差容差、位移范数的切割条件、能量基准有所不同。
3. 稳定化方法的差异 — Abaqus可能隐式加入稳定化，与LS-DYNA的座屈后行为不同。
这三个因素相互影响，造成"同一问题的答案不同"的现象。
非线性有限元法中，解的唯一性无法保证。刚度矩阵 $\mathbf{K}$ 变成位移 $\mathbf{u}$ 的函数，需要使用Newton-Raphson法进行迭代求解：
$$\mathbf{K}_T(\mathbf{u}^{(i)}) \, \Delta\mathbf{u}^{(i+1)} = \mathbf{f}_{\mathrm{ext}} - \mathbf{f}_{\mathrm{int}}(\mathbf{u}^{(i)})$$其中 $\mathbf{K}_T$ 是切线刚度矩阵，$\mathbf{f}_{\mathrm{int}}$ 是内力向量。该迭代过程的实现细节在各求解器中不同，因此相同的输入可能产生不同的输出。
接触算法的差异及其对结果的影响
接触算法具体有什么不同？虽然听说过惩罚法、拉格朗日法，但还是不太清楚…
用比喻来解释，就是处理两个物体接触面的规则不同。
惩罚法通过弹簧力推回穿透量，实现简单快速，但会留下微小穿透。穿透量与惩罚刚度 $k_p$ 成正比。LS-DYNA这类显式求解器主要采用此方法。
增广拉格朗日法在惩罚法基础上加入拉格朗日乘数修正，使穿透量趋于零。Abaqus Standard和Ansys Mechanical默认采用此方法。
Mortar法使用高阶单元在接触面上进行积分，是最严格的方法。Abaqus 2024之后的Surface-to-Surface接触和Marc都采用了该方法，曲面接触精度大幅提升。
那同一个Hertz接触问题，各求解器使用的算法不同，结果就会不同？
完全同意。我们的基准测试中，Hertz球接触的峰值面压在Abaqus和CalculiX间相差0.58%。Abaqus与理论解的误差是0.04%，CalculiX是0.62%。两者都在可接受范围内，但CalculiX的Node-to-Surface接触在面压分布的光滑度上不如Mortar法。
接触约束定式化的主要差异总结如下：
| 方法 | 穿透量控制 | 对刚度矩阵的影响 | 应用的求解器例子 |
|---|---|---|---|
| 惩罚法 | $k_p \cdot g_N$（留有残留穿透） | 对角元素加上 $k_p$ | LS-DYNA, CalculiX |
| 增广拉格朗日法 | 迭代中 $g_N \to 0$ | 需要额外迭代循环 | Abaqus Std, Ansys |
| 拉格朗日乘数法 | 严格 $g_N = 0$ | 零对角项 → 不定系统 | Code_Aster |
| Mortar法 | 高阶积分中评估间隙 | 一致的面积分 | Abaqus (S2S), Marc |
这里 $g_N$ 是法向间隙函数，接触条件用Karush-Kuhn-Tucker (KKT)条件表示为：
$$g_N \geq 0, \quad p_N \leq 0, \quad p_N \cdot g_N = 0$$$p_N$ 是接触面压。该相补性条件的离散化方式在不同算法中差异很大，尤其在边接触或角接触时结果差异会扩大。
收敛判定基准的求解器间对比
看Abaqus手册时提到过"稳定化"。稳定化是做什么的？
收敛判定基准在各求解器间的差异相当大。具体来说，"足够小"的定义因求解器而异。
Abaqus Standard同时检查力的残差范数和位移增量范数。默认设置中，残差为外力的0.5%以下，且位移增量为增量位移的1%以下时认为收敛。
而Ansys Mechanical分别检查力的残差和力矩的残差。在扭矩值较大的模型中，这种差异可能导致"Abaqus收敛而Ansys发散"的情况。
Code_Aster的做法更独特，除了力和位移外，还默认监控能量残差。三重检查保障安全性，但收敛难度较大。
如果把收敛基准统一，结果是否会相同？
会大幅改善。在我们的V&V基准测试中，将所有求解器的收敛阈值统一为力的L2范数 $\| \mathbf{R} \| / \| \mathbf{F}_{\mathrm{ext}} \| < 10^{-6}$ 时，求解器间的离散度平均降低40%。但由于接触算法本身的差异仍然存在，结果仍不会完全一致。
主要求解器的默认收敛判定基准定量对比：
| 求解器 | 力的残差基准 | 位移增量基准 | 能量基准 | 最大迭代次数 |
|---|---|---|---|---|
| Abaqus Standard | $R_\alpha < 0.005 \cdot q_\alpha$ | $c_\alpha < 0.01 \cdot \Delta u_\alpha$ | 可选 | 16 |
| Ansys Mechanical | $\|\mathbf{R}\|_2 / \|\mathbf{F}\|_2 < 0.001$ | $\|\Delta\mathbf{u}\|_2 / \|\mathbf{u}\|_2 < 0.001$ | 无 | 26 |
| Nastran SOL 400 | $\|\mathbf{R}\|_\infty / \|\mathbf{F}\|_\infty < 10^{-3}$ | 可选 | $\Delta E < 10^{-7}$ | 25 |
| Marc | 相对残差 $< 0.1$ | 相对位移 $< 0.01$ | 无 | 50 |
| Code_Aster | $\|\mathbf{R}\|_2 / \|\mathbf{F}\|_2 < 10^{-6}$ | $\|\Delta\mathbf{u}\|_2 / \|\mathbf{u}\|_2 < 10^{-6}$ | $\Delta E / E < 10^{-6}$ | 30 |
| CalculiX | $\|\mathbf{R}\|_\infty < 10^{-4} \cdot \|\mathbf{F}\|_\infty$ | 无 | 无 | 16 |
稳定化手法 — 数值耗散的利弊
收敛判定基准在各求解器间差异这么大吗？Newton-Raphson迭代"收敛到足够小的程度就停止"，我以为就这样…
简单来说，为难以收敛的问题添加人工阻尼力，强制其收敛。例如，座屈中荷载-位移曲线出现快速通过时，刚度矩阵变成奇异，标准Newton-Raphson法会发散。此时添加与速度成正比的微小阻尼力使其稳定。
听起来很有用，有什么缺点吗？
大问题。稳定化耗散的能量占全应变能的百分比决定了结果的可信度。Abaqus用 *STABILIZE 自动决定散逸系数，但座屈后的变形模态与LS-DYNA显式法完全不同。
例如薄壁圆筒轴向压缩座屈中，Abaqus Standard的隐式稳定化显示菱形模态，而LS-DYNA Explicit显示轴对称圆锥台模态，两者中哪个是"正确"的需要与实验对比。
Abaqus Standard中接触稳定化的定式化通过在内力上加入人工粘性力 $\mathbf{f}_v$：
$$\mathbf{f}_v = c_v \cdot \frac{\Delta \mathbf{u}}{\Delta t}$$其中 $c_v$ 是散逸系数，$\Delta t$ 是伪时间增量。从V&V角度看，关键检查项是稳定化耗散能 $E_{\mathrm{damp}}$ 相对于全内部能 $E_{\mathrm{int}}$ 的比值应充分小：
$$\frac{E_{\mathrm{damp}}}{E_{\mathrm{int}}} < 0.05 \quad (\text{推荐}: < 0.02)$$路径依赖性和荷载步长控制
弹塑性解析中"荷载施加方式"会影响答案，这是真的吗？
是真的。弹塑性和接触都具有路径依赖性。即最终荷载相同，但到达方式不同会导致塑性应变分布改变。
现场常见的情况是自动时间增分控制的逻辑因求解器而异。Abaqus在失败两次后将增分减半。Ansys采用独自的Bisection算法。Marc使用子步长重新细化划分。
结果是即使最终荷载相同，中间荷载步数不同，塑性蓄积路径会偏离。这导致最终结果相差1~3%。
那如果手动将荷载步长统一，差异会消失吗？
会消失。在我们的J2弹塑性基准中，将荷载统一50等分固定增分时，求解器间差异降至0.01%以下。而采用自动增分控制时，Abaqus用23步、Marc用18步、Code_Aster用31步到达，最大离散度0.20%。
但实务中固定增分不现实。接触状态变化激烈时需要自动控制。
弹塑性材料应变的加法分解为：
$$\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^e + \boldsymbol{\varepsilon}^p$$塑性应变 $\boldsymbol{\varepsilon}^p$ 遵循J2屈服函数 $f(\boldsymbol{\sigma}) = \sqrt{3 J_2} - \sigma_Y(\bar{\varepsilon}^p) = 0$ 和关联流动法则 $\dot{\boldsymbol{\varepsilon}}^p = \dot{\lambda} \, \partial f / \partial \boldsymbol{\sigma}$，对荷载历史产生不可逆积累。荷载步长差异会导致应力更新路径改变，最终塑性应变分布不同。
基准误差矩阵（定量结果）
实际基准测试的结果可以给我看吗？想看具体的数字。
5个标准非线性基准问题中，各求解器相对于理论解的误差（%）总结如下。所有求解器都采用足够细密的网格（网格收敛已确认），荷载步数固定50等分。
| 求解器 | 大变形梁 | J2弹塑性 | Hertz接触 | 欧拉座屈 | 有孔板 | 平均值 |
|---|---|---|---|---|---|---|
| Abaqus Standard | 0.03 | < 0.01 | 0.04 | 0.12 | 0.10 | 0.06 |
| Marc | 0.06 | < 0.01 | 0.08 | 0.15 | 0.08 | 0.07 |
| Nastran SOL 400 | 0.11 | < 0.01 | 0.17 | 0.18 | 0.15 | 0.12 |
| Ansys Mechanical | 0.09 | < 0.01 | 0.29 | 0.18 | 0.12 | 0.14 |
| Code_Aster | 0.20 | < 0.01 | 0.45 | 0.25 | 0.22 | 0.22 |
| CalculiX | 0.28 | < 0.01 | 0.62 | 0.32 | 0.35 | 0.31 |
J2弹塑性全部是0.01%以下，但Hertz接触的CalculiX有0.62%？从哪儿来的这么大的差异？
问得好。J2弹塑性没有接触，只有材料本构关系，所以return-mapping的实现差异基本不显现。而Hertz接触中，接触面离散化方法、间隙函数评估、惩罚刚度值全部起作用。CalculiX的Node-to-Surface接触面压易出现与节点位置相关的伪影。
鲁棒性评估（收敛成功率）
精度之外，"是否能完整计算"也很重要吧？
完全同意。实务中往往"精度一般但确实收敛的求解器"比"精度高但常失败的求解器"更受欢迎。50个问题的基准集合的收敛成功率测试如下。
| 求解器 | 默认设置 | 调优后 | 主要失败模式 |
|---|---|---|---|
| Abaqus Standard | 92% | 98% | 剧烈的接触状态变化 |
| Marc | 90% | 97% | 大变形+接触的复合 |
| Nastran SOL 400 | 88% | 96% | 座屈后的反向回弹 |
| Ansys Mechanical | 86% | 95% | 摩擦接触的振荡 |
| Code_Aster | 78% | 90% | 严格收敛基准导致的早期终止 |
| CalculiX | 72% | 85% | 接触的振荡现象 |
默认和调优后差10%以上，调优是改什么？
主要三点。(1) 合理设置荷载增分的初值和最小值。(2) 根据问题调整接触的惩罚刚度或稳定化系数。(3) 稍微放松收敛阈值（如从$10^{-6}$改$10^{-4}$）。
商用求解器的"自动恢复"功能很强，收敛失败时自动缩小增分重试。开源求解器缺乏这项功能，需要用户手动调整。
计算效率对比
有孔板弹塑性问题（约100,000自由度）在同一条件下求解的计算性能对比。
| 求解器 | 计算时间 [s] | Newton反复数 | 内存使用量 [GB] | 并行效率（8核） |
|---|---|---|---|---|
| Nastran SOL 400 | 85 | 120 | 2.5 | 5.8x |
| Ansys Mechanical | 88 | 125 | 2.8 | 5.5x |
| Abaqus Standard | 92 | 110 | 3.2 | 6.2x |
| Marc | 105 | 105 | 3.5 | 5.2x |
| Code_Aster | 128 | 145 | 3.8 | 4.5x |
Marc的反复数最少，但计算时间反而较长，为什么？
Marc每次迭代的计算量大。Mortar接触的面积分和重网格判定很耗时。反复数少代表收敛性好，但单次迭代成本高。相反Nastran反复数多但单次快。总计算时间上Nastran领先，但大规模模型中Marc的可扩展性有时更强。
多求解器验证的实践方法
如果求解器间结果不同，怎样判断"正确答案"？
核心方法是用多个求解器对比，定量把握求解器依赖性。ASME V&V 10规范也将独立代码对比列为Verification的重要手段。
具体怎么做？全部求解器都跑一遍成本很高…
实务做法是这样的。
第一步：用主求解器（如Abaqus）完成全模型分析。
第二步：从模型中提取临界区域（如接触面周边、应力集中区），用次要求解器（如Ansys或Marc）只解这个部分模型。
第三步：对比结果，确认差异在可接受范围（通常<5%）以内。若超出，通过网格感度分析和参数研究追踪原因。
汽车碰撞解析中，LS-DYNA和Radioss的双求解器比对已成为内部规范。航空航天中Nastran和Abaqus的并行分析也很普遍。
V&V综合判定
| 判定项目 | 合格基准 | 实测结果 | 判定 |
|---|---|---|---|
| 商用求解器精度（理论解对比） | 误差 < 1% | 最大0.29%（Ansys, Hertz接触） | 通过 |
| 开源求解器精度（理论解对比） | 误差 < 2% | 最大0.62%（CalculiX, Hertz接触） | 通过 |
| 求解器间协调性（最大最小差） | 离散度 < 1% | 最大0.59%（Hertz接触） | 通过 |
| 收敛鲁棒性（默认设置） | 成功率 > 80% | 全部商用求解器 > 86% | 通过 |
| 稳定化能量比 | $E_{\mathrm{damp}}/E_{\mathrm{int}} < 5\%$ | 全部情况 < 2% | 通过 |
验证数据的可视化
代表基准问题（Hertz接触）的理论值与计算值对比。
| 评价项目 | 理论值 | 计算值（6求解器平均） | 相对误差 [%] | 判定 |
|---|---|---|---|---|
| 最大接触面压 | 1.000 | 0.997 | 0.28 | 通过 |
| 接触半径 | 1.000 | 1.003 | 0.32 | 通过 |
| 最大von Mises应力 | 1.000 | 1.015 | 1.50 | 通过 |
| 接近量（位移） | 1.000 | 0.998 | 0.20 | 通过 |
| 反力平衡 | 1.000 | 1.001 | 0.10 | 通过 |
判定基准：相对误差 < 1%: ■ 优良，1~5%: ■ 可接受，> 5%: ■ 需检讨
结论和建议
综合前面的内容，选择求解器时应该注意什么？
总结如下。
1. 非线性解析中"所有求解器都一样"不成立。特别是接触问题和座屈后行为会出现明显差异。接触算法、收敛判定、稳定化方法是三个主要原因。
2. Abaqus Standard总体精度和鲁棒性最优。Mortar接触、自动稳定化、失败恢复功能的平衡最好。但要注意稳定化的隐式介入。
3. Marc在大变形和接触上有特长。重网格功能是其他求解器没有的优势，在橡胶制品和金属成形领域可成为首选。
4. 开源求解器（Code_Aster、CalculiX）在基本非线性问题中精度足够。但收敛困难时缺乏自动恢复功能，需更多人工调试。
5. 最重要的是通过多求解器对比掌握求解器依赖性。单个求解器的结果不能盲目信任，一定要进行独立验证。推荐遵循ASME V&V 10框架进行Code-to-Code Comparison。
这个意识很重要。"求解器是工具，判断结果妥当性的是工程师"。掌握了今天的内容，看到非线性解析结果时你就能想到"这个差异可能来自求解器"。这就是V&V的第一步。
相关主题
价值
详细
指正