非线性方程组可以表示为:
在位移为基本未知量的有限元分析中,是结点位移向量,是结点载荷向量。对于线性方程组,由于是常数矩阵,可以没有困难直接求解。对于非线性方程组需要迭代求解。
1 直接迭代法
只适用于与变形历史无关的非线性问题,例如非线性弹性问题,利用形变理论分析的弹塑性问题,稳态蠕变问题等。对于依赖于变形历史的非线性问题,应力需要由应变所经历的路径决定,直接迭代法不适用。例如加载路径不断变化或涉及卸载和反复加载等弹塑性问题。此时要利用增量理论。
可以指出,当是凸曲线(如图所示),其中是标量,即系统为单自由度情形,通常解收敛。
1.1直接迭代法的求解步骤
(1)方程可以改写为
选择一个初始的试探解代入上式可以得到初始的,接着可以得到第一步迭代的位移解,由此类推得到第次迭代后的近似解,一直到误差的某种范数小于某个规定的容许小量,即,上述迭代终止。
1.2常系数矩阵的直接迭代法
为避免每次迭代对刚度矩阵进行求逆计算,一般可以将刚度矩阵设定为常数进行迭代过程,单自由度系统下的常刚度直接迭代法如图所示
利用下式求的修正量
由此可以得到
同理可以得到
非线性材料重新形成的工作量远小于对其进行分解求逆的工作量(用于迭代的刚度矩阵和用于计算内力的刚度矩阵注意分开)
2 Newton-Raphson方法
一般情况下,(1)不能被精确地满足,即,为了得到进一步的近似解(假设为某单自由度方向上的位移) ,将表示成在附近的仅保留线性项的Taylor展开式
且有
式中是切线矩阵(tangent matrix),即
其迭代过程如图所示
为了克服N-R方法对于每次迭代需要形成并求逆一个新的切线矩阵所带来的麻烦,可以仿照直接迭代法采用一种修正的N-R方法。即如图所示
3 增量法
核心思想是将载荷分解成若干增量步,即,相应的位移也分为同样的步数,即,每两步之间的增长量称之为增量。增量解法的一般做法是假设第步的载荷和相应的位移为已知,然后将载荷增加为,如果每一步的载荷增量足够小,则解的收敛性可以保证。
使用增量法的(1)式改写为如下形式(位移均以单自由度为例)
,其中是用来表示载荷变化的参数,将上式对求导,可以得到
其中为定义的切线刚度矩阵。
(4)式是一典型的常微分方程组问题,下面介绍的是有限元中对这类方程组的求解方法
1.欧拉方法(单自由度为例)
显然为了满足精度要求,必须是足够小的量。使用Runge-Kutte方法可以改善精度
此方法会导致解的漂移(与真实曲线上的解产生较大的误差)
为了克服每一增量步解的误差可能导致的解的漂移,可以将(5)式改写为
这里解释一下
此方法称为考虑平衡矫正的欧拉增量方法,即将前增量步的载荷和内力的不平衡量合并到当前增量步求解,一定程度上避免了解的漂移。
2.N-R方法
为了改进欧拉法的精度,现在更多采用N-R方法,如果采用N-R方法,是在每一增量步内进行迭代
则对于的次增量步的第次迭代可以表示为
表示成前述的N-R形式为
采用mN-R迭代
4 加速收敛的方法(Aitken加速法)
利用mN-R方法求解非线性方程组时,可以避免每次迭代重新形成和求逆切线矩阵,但是降低了收敛速度,特别是曲线突然趋于平坦的情况
有Aitken加速的迭代和无Aitken加速的迭代如图所示
核心原理是将通过初始切线刚度得到的第m+1个增量步的第n个迭代解利用局部割线刚度转换出新的迭代解,以说明此转换过程
可以改写为
被称作加速因子
Aitken加速收敛方法是每隔一次迭代进行一次加速(为了保证精度的同时提高收敛速度)。
推广到个自由度系统,Aitken方法可以表示为
其中是对角矩阵,其元素
5载荷增量步长的自动选择
在研究载荷增量步长自动选择的方法时,首先是假设载荷的分布模式是给定的,变化的只是他的幅值,在此情况下,外载荷可以表示为:
等效结点载荷向量可以表示为:
是载荷幅值,载荷分布实际是的分布
求解上述载荷分布方程依然是一个N-R迭代过程,设
是经过n+1次迭代修正后得到的数值
下面是几种常用的自动选择载荷步长的方法
5.1规定”本步刚度参数“的变化量以控制载荷增量
此方法对于计算结构的极限载荷很有效,利用本步刚度参数可以使步长调整的比较合理,并可以减少总的增量步数,适合于计算由理想弹塑性材料组成结构的极限载荷情况。
5.1.1 ”本步刚度参数“的概念
第i增量步结构的总体刚度可用下式度量
初始结构总体刚度的度量是
其中和是载荷向量和按弹性分析得到的位移向量。
5.1.2 增量步长的自动选择
(1) 利用”本步刚度参数“的规定变化量自动选择增量步长。
称为第i步的”本步刚度参数“,它代表结构本身的刚度性质,与载荷增量大小无关,当结构处于完全弹性时,,随着塑性区扩大,结构变软,逐渐减小,到达极限载荷时,
对于比例加载的情况:
是极限载荷参数,是时的结点载荷向量
利用本步刚度参数可以使步长调整得比较合理,并可以减少总的增量步
(2) 增量步长的自动选择
利用(1)得到的”本步刚度参数“的规定变化量自动选择增量步长,具体的算法步骤如下
1、求弹性极限载荷参数
先施加任意载荷,假定结构为完全弹性求解,求出结构的最大等效应力,则有
是材料的初始屈服应力
2、给定第一步载荷参数增量 的值可以事先给定,用求解第一增量步
3、给定第二及以后各增量步的刚度参数变化的预测值,其大小决定步长的大小,例如可在0.05~0.2之间选择,并给定刚度的最小允许值,则每步载荷参数增量为
然后用求解第增量步
然后可以通过(1)中的公式计算出”本步刚度参数“
从图中可见是在给定预测值情况下,利用线性外推得到的,为了得到更好的适应复杂的曲线,可以使用二次外推。
5.2规定某个结点的位移增量以确定载荷增量
以迭代为例,它可以表示为
在载荷增量步长自动控制的求解方法中,可以表示成
其中
以上两式中:
因此(7)可以改写为
其中
是在第次迭代中由某个规定的约束条件来确定的载荷因子增量的第次修正量。在现在的方法中,这个条件就是某个结点的位移增量的大小,例如规定中的最大分量是给定的,此条件可以表示为
其中是的规定值,是除第g个元素为1,其余元素为零的向量,具体迭代算法步骤如下:
(1)计算对于节点载荷模式的位移模式,即
(2)计算对于不平衡结点力向量的位移增量修正值和次迭代后位移增量修正值的全量,即
其中是待定载荷因子增量的修正值
(3)利用条件(8)确定,由于
从上式可以得到
这样就确定了第次迭代后的载荷因子增量
载荷因子增量求出之后,可以知道每一步的外载荷,从而使用法迭代得到位移,进而求出应变,应力和当前增量步的内力等物理量。并检验收敛准则是否满足,如未满足,回到步骤二进行新的一次迭代,直至收敛准则满足为止。关于每一个增量步某个指定结点位移增量本身的选择,通常的方法是第1个增量步可以由某个给定的载荷因子增量(例如令,其中p_e是弹性极限载荷因子,可取5~10),通过求解得到,以后各增量步的可由下式确定,即
其中
6 小结
以上是王勖成的有限单元法给出的非线性方程组解法以及一些提高运算效率的策略,如有补充或者理解偏差,请联系指正。