Question

非线性方程组解法

Answer 1

非线性方程组可以表示为：

在位移为基本未知量的有限元分析中， 是结点位移向量， 是结点载荷向量。对于线性方程组 ，由于 是常数矩阵，可以没有困难直接求解。对于非线性方程组需要迭代求解。

只适用于与变形历史无关的非线性问题，例如非线性弹性问题，利用形变理论分析的弹塑性问题，稳态蠕变问题等。对于依赖于变形历史的非线性问题，应力需要由应变所经历的路径决定，直接迭代法不适用。例如加载路径不断变化或涉及卸载和反复加载等弹塑性问题。此时要利用增量理论。
可以指出，当 是凸曲线（如图所示），其中 是标量，即系统为单自由度情形，通常解收敛。

(1)方程可以改写为
选择一个初始的试探解 代入上式可以得到初始的 ，接着可以得到第一步迭代的位移解 ，由此类推得到第 次迭代后的近似解 ，一直到误差的某种范数小于某个规定的容许小量 ，即 ，上述迭代终止。

为避免每次迭代对刚度矩阵 进行求逆计算，一般可以将刚度矩阵设定为常数进行迭代过程，单自由度系统下的常刚度直接迭代法如图所示

一般情况下，（1）不能被精确地满足，即 ，为了得到进一步的近似解 （假设 为某单自由度方向上的位移） ，将 表示成在 附近的仅保留线性项的Taylor展开式

且有
式中 是切线矩阵(tangent matrix)，即



其迭代过程如图所示

核心思想是将载荷分解成若干增量步，即 ，相应的位移也分为同样的步数，即 ，每两步之间的增长量称之为增量。增量解法的一般做法是假设第 步的载荷 和相应的位移 为已知，然后将载荷增加为 ，如果每一步的载荷增量 足够小，则解的收敛性可以保证。
使用增量法的（1）式改写为如下形式（位移均以单自由度为例）
,其中 是用来表示载荷变化的参数，将上式对 求导，可以得到


其中 为定义的切线刚度矩阵。
（4）式是一典型的常微分方程组问题，下面介绍的是有限元中对这类方程组的求解方法
1.欧拉方法(单自由度为例)

显然为了满足精度要求， 必须是足够小的量。使用Runge-Kutte方法可以改善精度
此方法会导致解的漂移（与真实曲线上的解产生较大的误差）
为了克服每一增量步解的误差可能导致的解的漂移，可以将（5）式改写为

这里解释一下
此方法称为考虑平衡矫正的欧拉增量方法，即将前增量步的载荷和内力的不平衡量合并到当前增量步求解，一定程度上避免了解的漂移。

2.N-R方法
为了改进欧拉法的精度，现在更多采用N-R方法，如果采用N-R方法，是在每一增量步内进行迭代
则对于 的 次增量步的第 次迭代可以表示为

表示成前述的N-R形式为
采用mN-R迭代

利用mN-R方法求解非线性方程组时，可以避免每次迭代重新形成和求逆切线矩阵，但是降低了收敛速度，特别是 曲线突然趋于平坦的情况
有Aitken加速的迭代和无Aitken加速的迭代如图所示

在研究载荷增量步长自动选择的方法时，首先是假设载荷的分布模式是给定的，变化的只是他的幅值，在此情况下，外载荷可以表示为：

等效结点载荷向量可以表示为：

是载荷幅值，载荷分布实际是 的分布

求解上述载荷分布方程依然是一个N-R迭代过程，设

是经过n+1次迭代修正后得到的 数值
下面是几种常用的自动选择载荷步长的方法

此方法对于计算结构的极限载荷很有效，利用本步刚度参数可以使步长调整的比较合理，并可以减少总的增量步数，适合于计算由理想弹塑性材料组成结构的极限载荷情况。

第i增量步结构的总体刚度可用下式度量

初始结构总体刚度的度量是

其中 和 是载荷向量和按弹性分析得到的位移向量。

(1) 利用”本步刚度参数“的规定变化量自动选择增量步长。

称为第i步的”本步刚度参数“，它代表结构本身的刚度性质，与载荷增量大小无关，当结构处于完全弹性时， ，随着塑性区扩大，结构变软， 逐渐减小，到达极限载荷时，
对于比例加载的情况：
　　　　　
是极限载荷参数， 是 时的结点载荷向量
利用本步刚度参数可以使步长调整得比较合理，并可以减少总的增量步
(2) 增量步长的自动选择
利用(1)得到的”本步刚度参数“的规定变化量自动选择增量步长，具体的算法步骤如下
1、求弹性极限载荷参数
先施加任意载荷 ，假定结构为完全弹性求解，求出结构的最大等效应力 ，则有

是材料的初始屈服应力
2、给定第一步载荷参数增量 的值可以事先给定，用 求解第一增量步
3、给定第二及以后各增量步的刚度参数变化的预测值 ，其大小决定步长的大小，例如可在0.05~0.2之间选择，并给定刚度的最小允许值 ，则每步载荷参数增量为

然后用 求解第 增量步
然后可以通过(1)中的公式计算出”本步刚度参数“

以 迭代为例，它可以表示为

在载荷增量步长自动控制的求解方法中， 可以表示成

其中

以上两式中:
　　　　　　　
因此（7）可以改写为
其中
是在第 次迭代中由某个规定的约束条件来确定的载荷因子增量 的第 次修正量。在现在的方法中，这个条件就是某个结点的位移增量的大小，例如规定 中的最大分量 是给定的，此条件可以表示为

其中 是 的规定值， 是除第g个元素为1，其余元素为零的向量，具体迭代算法步骤如下：
(1)计算对于节点载荷模式 的位移模式 ，即

(2)计算对于不平衡结点力向量 的位移增量修正值 和 次迭代后位移增量修正值的全量 ，即


其中 是待定载荷因子增量的修正值
(3)利用条件(8)确定 ，由于

从上式可以得到

这样就确定了第 次迭代后的载荷因子增量
载荷因子增量求出之后，可以知道每一步的外载荷，从而使用 法迭代得到位移，进而求出应变，应力和当前增量步的内力等物理量。并检验收敛准则是否满足，如未满足，回到步骤二进行新的一次迭代，直至收敛准则满足为止。 关于每一个增量步某个指定结点位移增量 本身的选择，通常的方法是第1个增量步可以由某个给定的载荷因子增量 (例如令 ，其中p_e是弹性极限载荷因子， 可取5~10)，通过求解得到 ，以后各增量步的 可由下式确定，即

其中

以上是王勖成的有限单元法给出的非线性方程组解法以及一些提高运算效率的策略，如有补充或者理解偏差，请联系指正。