Abaqus 隐式求解到底在算什么？从 Newton-Raphson 法说起

2月前浏览189

在Abaqus 中的应力值到底是怎么计算出来的？一文中我们介绍到Abaqus 非线性隐式问题（Standard）的核心算法是Newton-Raphson迭代法（以下均称“N-R 法”）。如果你翻开传统的有限元教材（如图1），关于该算法的介绍往往充斥着密密麻麻的多维矩阵推导，对许多初学者来说，这可能有点跳跃，难以理解。今天这篇文章我希望用通俗易懂的形式为大家讲清楚N-R迭代法。

图1 有限元教材中N-R法介绍

N-R法是用来算什么的？

一句话回答：它是用来解复杂方程的（求方程的根）。

我们在初中学的方程，比如
，一眼看出。这叫线性方程。

但工程中，遇到的多是非线性方程，比如：

这种方程你无法用公式直接把提出来算，只能靠“猜”。而在Abaqus非线性分析中，方程变成了（刚度矩阵 × 位移 = 外力，结构刚度也是变化的），非线性方程更加复杂。

N-R法的核心目的，就是教计算机如何“有策略地猜”，并且用最快的速度猜中准确答案。

计算机是怎么“猜”的？

牛顿法的核心几何思想只有七个字：“以直代曲，顺藤摸瓜”。

想象你在一座地形复杂的山上，你的任务是走到海拔为0的海岸线（即求解
）。那怎样走最快呢？

定起点：随便找个位置站好，测出当前高度。
探竹竿：顺着脚下的坡度（导数
），用一根笔直的长竹竿向山下斜探过去。
定新点：竹竿扎中海平面的位置，就是你下一步要去的坐标。
重复走：走到，顺着新坡度再次探出竹竿……重复几次，你就能迅速且精准地踩中海岸线！

这就是神奇的N-R迭代公式：

：你现在所在的横坐标（当前的猜测值）。

：你现在的海拔高度（当前的误差值，目标是 0）。

：你脚下的坡度（函数的导数）。

：竹竿扎中海平面的新位置（下一次的猜测值）。

算例：求解

如何用N-R法求解方程：？

第一步：准备函数和它的导数（坡度）

函数：

导数：

第二步：代入牛顿法公式

每次往下走一步的公式就是：

第三步：开始迭代（猜）

第0次（瞎猜）：我猜答案是。
此时误差为（距离 0 还差很远）。
第 1 次迭代：根据公式算出下一步该去哪：

此时误差为（离 0 近多了！）。
第 2 次迭代：

以 1.65 为起点继续走：

此时误差约为 0.0479 （误差不到十分之一了！）。
第 3 次迭代：

继续走：

此时误差约为 0.00028（接近0！）。

图2 算例N-R法求解过程

结论：你看，一开始瞎猜是2.5，仅仅用了3次套公式，得出的结果就和计算器按出来的几乎一模一样！这就是N-R法闻名天下的原因——收敛速度极快。

映射回有限元：

Abaqus中N-R法在算什么？

理解了上面这个简单的例子，你就能明白 Abaqus中N-R法在计算什么了。

Abaqus中N-R法相关变量

：物体的变形（位移）。
：不平衡力（外部拉力减去物体内部的抵抗力）。Abaqus 的目标是让不平衡力等于 0。
：物体的刚度（你拉它一下，它反抗的强烈程度）。

Abaqus中N-R法迭代过程

给定一个外力，Abaqus 猜一个变形量。
算一下现在内外力平不平衡？如果不平衡（误差值就是残余力）。
Abaqus 计算一下结构现在的“刚度”（坡度）。
残余力除以刚度，得出下一步该增加多少变形量。
一次次逼近，直到内外力平衡（误差小于规定公差），这一步就算收敛了，然后继续算下一增量步的受力！

(注：回头看有限元书上的多维矩阵，其实就是把这里 1 维的变成了包含几十万个节点的数组，思想是一样的！)

图3 从N-R法映射到 Abaqus 隐式求解

疑问：N-R法有没有可能“翻车”？

看到这里你可能会问：既然牛顿法这么快，那是不是每次都能找到答案？迭代后的误差有没有可能不仅没变小，反而越算越大？

绝对有可能！这就是你的 Abaqus 经常报错“不收敛”的罪魁祸首。

在数值计算中，误差不仅可能变大，甚至可能变成天文数字。这种现象在数学上称为发散。N-R法极易在以下三种情况下“翻车”：

踩到了“平地”（刚度趋于 0）

数学现象：如果迭代时刚好走到一个平缓的山谷底部，切线坡度几乎等于 0。公式中除以一个极小数，会导致下一步被“踢”到无穷远。
工程对应：当材料发生开裂、屈服或彻底断裂的瞬间，丧失了抵抗变形的能力，刚度瞬间突变为 0 甚至负数。此时软件会算出巨大的位移增量，导致误差失控（常伴随“数值奇异 Numerical Singularity”报错）。

2. 陷入了“无限死循环”（震荡现象）

数学现象：由于函数形状特殊，顺着切线走到另一边后，新的切线又把你送回起点附近，在 A 点和 B 点之间来回横跳。
工程对应：在接触（Contact）分析中极常见。软件发现穿透，施加巨大排斥力排斥力过大把物体弹开判定无接触并撤销排斥力物体再次靠拢穿透……陷入死循环。

3. 初始位置“猜得太离谱”

工程对应：在仿真中，这通常是因为时间步长（Time Increment）设置过大。步子迈得太大导致状态发生突变，初始猜测完全偏离真实解，切线顺着错误的坡度越跑越远。

N-R法失效了怎么办？

Abaqus中面对Newton-Raphson法发散失效采用了两种措施：

限制最大迭代次数：比如每次增量步默认最多允许迭代 15 次。如果发现误差不断变大（出现 1U, 2U... 等警告）或陷入横跳震荡，Abaqus 会立刻强制停止该步的尝试。
减小步长：当检测到发散时，Abaqus 会自动退回上一个稳定状态，将时间步长切碎（比如把 0.1 秒缩小一半变成 0.05 秒）。这就相当于让每次前进的步伐变小，使初始猜测点离真实答案更近一点，试图避开危险的“突变区”重新开始迭代。如果还是不行，就继续切碎（也就是大家常见的报错直至放弃）。