有限元分析：非线性分析求解算法——牛顿迭代法
结构有限元分析中，涉及到弹塑性材料、结构大变形大转动以及接触问题都是非线性问题，需要使用非线性工况进行求解。
Abaqus、OptiStruct及Ansys Mechanical等有限元求解器多使用牛顿迭代法，即Newton-Raphson迭代法，来求解非线性问题。
很多结构非线性培训资料都会用下面这张图来讲解牛顿迭代法，但说实话，笔者第一次看到这张图时，只知道是逐步逼近准确解，知其然但不知其所以然，不知道怎么来的。
牛顿迭代法有以下三种变形：常规牛顿迭代法，修正的牛顿迭代法以及牛顿下山法。有限元求解器一般会根据实际情况自动选择最合适的那种，帮助文档中一般不会讲得很细。
下面请听我一一道来。
1、常规牛顿迭代法
若不存在加速度和阻尼，有限元算法可以简化为求解平衡方程 KX=F ，其中刚度矩阵K和载荷F为已知条件，通常需要计算刚度矩阵K的逆矩阵从而求解位移矩阵X。
常规牛顿迭代法如下图所示，注意曲线L(u)并不是已知的，而是逐步求解出来的，下面详细说明求解过程。
假设分析步总载荷为100N，非线性算法中一般会分多次加载，比如分5次加载每次增加20N，5次加载过程称为5个增量步；每次加载时需要多次迭代才会达到平衡状态，称为迭代步，迭代时使用的算法就是牛顿迭代法。
假设在第n个迭代步时，需要加载大小为 f 的力，下面开始迭代，迭代顺序依次为：A>B>C>D>E。
- 初始位移为上图中的A点， 对应横坐标为u_{0} ，根据材料参数和结构形状可以计算出A点对应的刚度矩阵 K_{0} ，根据方程 K_{0}X=f 可得到位移 u_{1} ，从而确定B点坐标( u_{1} ，f )；
- 根据B点位移 u_{1} 、材料参数以及结构形状可计算出此状态下结构的刚度矩阵 K_{1} ，根据方程 f_{1}=K_{1}u_{1} 求得C点坐标( u_{1} ，f_{1} )以及残差 R_{1}=f-f_{1} ;
- 以C点为起点，对应刚度矩阵为 K_{1} ，加载大小为 R_{1} 的力，可得到D点坐标( u_{2} ，f )；
- 根据D点位移 u_{2} 、材料参数以及结构形状可计算此状态下结构的刚度矩阵 K_{2} ，根据方程 f_{2}=K_{2}u_{2} 求得E点坐标( u_{2} ，f_{2} )以及残差 R_{2}=f-f_{2} ，若残差 R_{2} 小于设定误差值，则认为已经达到收敛，若残差 R_{2} 太大，则继续迭代。
觉得不错的话，点个赞吧，更多有限元相关硬核分享，请点击头像，查看我的主页 @熊库辛
从数学上来说，牛顿迭代法是数值计算方法，具体可翻阅数值分析教材，下面简单说明其数学原理。
设已知方程 f(x)=0 有近似根 x_{k} (假定 f'(x_{k})\ne0 )，根据泰勒级数将函数 f(x) 在 x_{k} 点展开，有
f(x)\approx f(x_{k})+f'(x_{k})(x-x_{k})
于是方程 f(x)=0 可近似地表示为
f(x_{k})+f'(x_{k})(x-x_{k})=0
根据该方程可求出迭代公式 x_{k+1}=x=x_{k}-\frac{f(x_{k})}{f'(x_{k})} ，其中k=0,1,...，取接近于方程根的一个初始值，然后开始迭代，一般经过几次迭代就能得到一个精度较高的解。
常规牛顿迭代法的优点是收敛速度快。
2、修正的牛顿迭代法
常规牛顿迭代法每次迭代都是求解新的刚度矩阵K以及K的逆矩阵，计算量较大，如果刚度矩阵总是使用初始的 K_{1} ，并且保持不变，则可以不求解求逆矩阵，大大减少了计算量，这就是修正的牛顿迭代法，示意图如下。
3、牛顿下山法
常规牛顿迭代法要求初始值必须在准确值附近才会收敛，初始值不合适可能会导致结果不收敛。
为克服收敛问题，引入判断条件 \left| f(x_{k+1}) \right|<\left| f(x_{k}) \right| ，满足此要求的算法称为下山法。
常规牛顿迭代法迭代公式为
\bar{x}_{k+1}=x_{k}-\frac{f(x_{k})}{f'(x_{k})} (1)
若对计算结果进行加权平均作为新的改进值，
x_{k+1}=\lambda \bar{x}_{k+1}+(1-\lambda)x_{k} (2)
其中 0<\lambda\leq 1 ，称为下山因子。
将(1)式代入(2)式有
x_{k+1}=x_{k}-\frac{\lambda f(x_{k})}{f'(x_{k})} ， k=0,1,2...
选择下山因子时从 \lambda=1 开始，逐次将 \lambda 减半进行试算，直到能使下降条件 \left| f(x_{k+1}) \right|<\left| f(x_{k}) \right| 满足为止。
通过引入下山因子，保证了牛顿迭代法的收敛性。
因此牛顿下山法收敛速度快，一定收敛，而且对初始值要求不高，相比常规牛顿迭代法更稳健。
更多有限元相关硬核分享，请点击头像，查看我的主页 @熊库辛
码字不易，点个赞吧。
参考文献：
[1]李庆扬，王能超，易大义. 数值分析[M]. 第5版. 北京: 清华大学出版社, 2008:222-226.
[1]曾攀. 有限元分析及应用[M]. 北京: 清华大学出版社, 2004:325-327.
