非线性方程组解法

非线性方程组可以表示为:
\pmb \psi(\pmb a)=\pmb P(\pmb a)-\pmb Q=0    (1)
在位移为基本未知量的有限元分析中,\pmb a是结点位移向量,\pmb Q是结点载荷向量。对于线性方程组\pmb K \pmb a=\pmb Q,由于\pmb K是常数矩阵,可以没有困难直接求解。对于非线性方程组需要迭代求解。

1 直接迭代法

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

3.PNG

1.1直接迭代法的求解步骤

(1)方程可以改写为\pmb K(\pmb a) \pmb a=\pmb Q
选择一个初始的试探解\pmb a=\pmb a^{(0)}代入上式可以得到初始的\pmb K^{(0)},接着可以得到第一步迭代的位移解\pmb a^{(1)}=\pmb (K^{(0)})^{-1} \pmb Q,由此类推得到第n次迭代后的近似解\pmb a^{(n)}=\pmb (K^{(n-1)})^{-1} \pmb Q,一直到误差的某种范数小于某个规定的容许小量er,即||e||=||\pmb a^{(n)}-\pmb a^{(n-1)}||≤er,上述迭代终止。

1.2常系数矩阵的直接迭代法

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

4.PNG

利用下式求a^{(1)}的修正量\Delta a^{(1)}
\Delta a^{(1)}=\pmb (K^{(0)})^{-1}(\pmb Q-(\pmb K^{(1)}a^{(1)}) )    (2)
由此可以得到a^{(2)} = a^{(1)}+\Delta a^{(1)}
同理可以得到\Delta a^{(n-1)} = (\pmb K^{(0)})^{-1}(\pmb Q-(\pmb K^{(n-1)} a^{(n-1)}) )
a^{(n)} = a^{(n-1)}+\Delta a^{(n-1)}
非线性材料重新形成\pmb K^{(n-1)}的工作量远小于对其进行分解求逆的工作量(用于迭代的刚度矩阵和用于计算内力的刚度矩阵注意分开)

2 Newton-Raphson方法

一般情况下,(1)不能被精确地满足,即\pmb \psi(a^{(n)}) \neq 0,为了得到进一步的近似解a^{(n+1)}(假设a^{(n+1)}为某单自由度方向上的位移) ,将\pmb \psi(a^{(n+1)})表示成在a^{(n)}附近的仅保留线性项的Taylor展开式
\pmb \psi(a^{(n+1)})=\pmb \psi(a^{(n)})+(\frac{d\pmb\psi}{da})_{n}\Delta a^{(n)}=0  (3)
且有a^{(n+1)} = a^{(n)}+\Delta a^{(n)}
式中d\psi/da是切线矩阵(tangent matrix),即
\frac{d\pmb\psi}{da}=\frac{d\pmb P}{da}=\pmb K_T(a)
\Delta a^{(n)} = (\pmb K_T^{(n)})^{-1}(\pmb Q-\pmb P^{(n)})
\pmb K_T^{(n)}=\pmb K_T(a^{(n)}) \pmb P^{(n)}=\pmb P(a^{(n)})
其迭代过程如图所示

5.PNG

为了克服N-R方法对于每次迭代需要形成并求逆一个新的切线矩阵所带来的麻烦,可以仿照直接迭代法采用一种修正的N-R方法。即\pmb K_T^{(n)}=\pmb K_T^{(0)})如图所示
6.PNG

3 增量法

核心思想是将载荷分解成若干增量步,即\pmb Q_0,\pmb Q_1,\pmb Q_2,\pmb Q_3,\pmb Q_4......,相应的位移也分为同样的步数,即\pmb a_0,\pmb a_1,\pmb a_2,\pmb a_3,\pmb a_4......,每两步之间的增长量称之为增量。增量解法的一般做法是假设第m步的载荷\pmb Q_m和相应的位移a_m为已知,然后将载荷增加为\pmb Q_{m+1}(=\pmb Q_{m}+\Delta \pmb Q_{m}),如果每一步的载荷增量\Delta \pmb Q_{m}足够小,则解的收敛性可以保证。
使用增量法的(1)式改写为如下形式(位移均以单自由度为例)
\pmb \psi(a)=\pmb P(a)- \lambda \pmb Q_0=0,其中\lambda是用来表示载荷变化的参数,将上式对\lambda求导,可以得到
\frac{d\pmb P}{da}\frac{da}{d\lambda}-\pmb Q_0=\pmb K_T \frac{da}{d\lambda}=0
\frac{da}{d\lambda}=\pmb K_T^{-1}(a)\pmb Q_0    (4)
其中\pmb K_T为定义的切线刚度矩阵。
(4)式是一典型的常微分方程组问题,下面介绍的是有限元中对这类方程组的求解方法
1.欧拉方法(单自由度为例)
a_{m+1}-a_m=\Delta a_m=\pmb K_T^{-1}(a_m)\pmb Q_0 \Delta \lambda_m= (\pmb K_T)^{-1}_{m} \Delta \pmb Q_m    (5)
显然为了满足精度要求,\Delta \lambda_m必须是足够小的量。使用Runge-Kutte方法可以改善精度
此方法会导致解的漂移(与真实曲线上的解产生较大的误差)
为了克服每一增量步解的误差可能导致的解的漂移,可以将(5)式改写为
a_{m+1}-a_m=\Delta a_m=\pmb K_T^{-1}(a_m)(\pmb Q_{m+1}-\pmb P(a_m))   (6)
这里解释一下(\pmb Q_{m+1}-\pmb P(a_m))=(\pmb Q_{m+1}-\pmb Q_{m}+\pmb Q_{m}-\pmb P(a_m))
此方法称为考虑平衡矫正的欧拉增量方法,即将前增量步的载荷和内力的不平衡量合并到当前增量步求解,一定程度上避免了解的漂移。

7.PNG

2.N-R方法
为了改进欧拉法的精度,现在更多采用N-R方法,如果采用N-R方法,是在每一增量步内进行迭代
则对于\lambdam+1次增量步的第n+1次迭代可以表示为
\psi_{m+1}^{(n+1)}=\pmb P(a_{m+1}^{(n+1)})-\pmb Q_{m+1}=\pmb P(a_{m+1}^{(n)})-\pmb Q_{m+1}+(\pmb K_T^n)_{m+1}\Delta a_m^n=0   (6)
表示成前述的N-R形式为\Delta a_m^n=(\pmb K_T^n)_{m+1}^{-1}(\pmb Q_{m+1}-\pmb P(a_{m+1}^{(n+1)}))
采用mN-R迭代
(\pmb K_T^n)_{m+1}=(\pmb K_T^0)_{m+1}=\pmb K_T(a_m)

8.PNG
9.PNG

4 加速收敛的方法(Aitken加速法)

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

10.PNG

11.PNG

核心原理是将通过初始切线刚度\pmb K_T^0得到的第m+1个增量步的第n个迭代解\Delta a_{m}^{n}利用局部割线刚度转换出新的迭代解\Delta \widetilde a_{m}^{n},以n=0说明此转换过程
K_s \Delta a_{m}^{0}=(K_T)_m(\Delta a_{m}^{0}-\Delta a_{m}^{1})
可以改写为
\frac{(K_T)_m}{K_s}=\frac{\Delta a_{m}^{0}}{(\Delta a_{m}^{0}-\Delta a_{m}^{1})}=\alpha^{(1)}
K_s \Delta \widetilde a_{m}^{1}=(K_T)_m\Delta a_{m}^{1}
\Delta \widetilde a_{m}^{1}=\frac{(K_T)_m}{K_s}\Delta a_{m}^{1}=\alpha^{(1)}\Delta a_{m}^{1}
\alpha^{(1)}被称作加速因子
Aitken加速收敛方法是每隔一次迭代进行一次加速(为了保证精度的同时提高收敛速度)。
推广到N个自由度系统,Aitken方法可以表示为
\pmb a_{m+1}^{(n+1)}=\pmb a_{m+1}^{n}+\pmb \alpha^{(n)}\Delta \pmb a_{m}^{n}
其中\pmb \alpha^{(n)}是对角矩阵,其元素\alpha_i^{(n)}=\left\{\begin{matrix}1 \\ \frac{\Delta a_{i,m}^{n-1}}{\Delta a_{i,m}^{n-1}-\Delta a_{i,m}^{n}}\end{matrix}\right.

5载荷增量步长的自动选择

在研究载荷增量步长自动选择的方法时,首先是假设载荷的分布模式是给定的,变化的只是他的幅值,在此情况下,外载荷可以表示为:
^t \overline {\pmb F} ={^t}P \overline{\pmb F_0}
等效结点载荷向量可以表示为:
^t \pmb Q_l ={^t}P \pmb Q_0
P是载荷幅值,载荷分布实际是P(t)的分布
^{t+\Delta t}P = {^t}P+\Delta P
求解上述载荷分布方程依然是一个N-R迭代过程,设^{t+\Delta t}P^{(0)}={^t}P
^{t+\Delta t}P^{(n+1)} = ^{t+\Delta t}P^{(n)}+\Delta P^{(n)}
^{t+\Delta t}P^{(n+1)}是经过n+1次迭代修正后得到的^{t+\Delta t}P数值
下面是几种常用的自动选择载荷步长的方法

5.1规定”本步刚度参数“的变化量以控制载荷增量

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

5.1.1 ”本步刚度参数“的概念

第i增量步结构的总体刚度可用下式度量
S_p^{(i)*}=\frac{\Delta \pmb Q^{(i)T}\Delta \pmb Q^{(i)}}{\Delta \pmb a^{(i)T} \Delta \pmb Q^{(i)}}
初始结构总体刚度的度量是
S_p^{(0)*}=\frac{\pmb Q^{(e)T}\pmb Q^{(e)}}{\pmb a_e^{T} \pmb Q_e}
其中\pmb Q_e\pmb a_e是载荷向量和按弹性分析得到的位移向量。

5.1.2 增量步长的自动选择

(1) 利用”本步刚度参数“的规定变化量自动选择增量步长。
S_p^{(i)}= \frac{S_p^{(i)*}}{S_p^{(0)*}}
S_p^{(i)}称为第i步的”本步刚度参数“,它代表结构本身的刚度性质,与载荷增量大小无关,当结构处于完全弹性时,S_p^{(i)}=1,随着塑性区扩大,结构变软,S_p^{(i)}逐渐减小,到达极限载荷时,S_p^{(i)}=0
对于比例加载的情况:
\pmb Q_e=p_e\pmb Q_0     \Delta \pmb Q^{(i)}=\Delta p_i\pmb Q_0
p_e是极限载荷参数,\pmb Q_0p=1时的结点载荷向量
利用本步刚度参数可以使步长调整得比较合理,并可以减少总的增量步
(2) 增量步长的自动选择
利用(1)得到的”本步刚度参数“的规定变化量自动选择增量步长,具体的算法步骤如下
1、求弹性极限载荷参数p_e
先施加任意载荷p \pmb Q_0,假定结构为完全弹性求解,求出结构的最大等效应力\overline \sigma_{max},则有
p_e=\frac{\sigma_{s0}}{\overline \sigma_{max}}
\sigma_{s0}是材料的初始屈服应力
2、给定第一步载荷参数增量\Delta p_1=p_e/N, N的值可以事先给定,用p_1=p_e+\Delta p_1求解第一增量步
3、给定第二及以后各增量步的刚度参数变化的预测值\Delta \widetilde S_p,其大小决定步长的大小,例如可在0.05~0.2之间选择,并给定刚度的最小允许值S_p^{min}(S_p^{min}\neq 0),则每步载荷参数增量为
\Delta p_i=\Delta p_{i-1} \frac{min(\Delta \widetilde S_p,S_p^{(i-1)}-S_p^{min})}{\lvert S_p^{(i-2)}-S_p^{(i-1)}\rvert}    (i=2,3,......)
然后用p_i=p_{i-1}+\Delta p_i求解第i增量步
然后可以通过(1)中的公式计算出”本步刚度参数“

12.PNG

从图中可见\Delta p_i是在给定预测值\Delta \widetilde S_p情况下,利用\Delta p_{i-1}/{\lvert S_p^{(i-2)}-S_p^{(i-1)}\rvert}线性外推得到的(y=kx),为了得到更好的适应复杂的曲线,可以使用二次外推。

5.2规定某个结点的位移增量以确定载荷增量

mN-R迭代为例,它可以表示为
^\tau \pmb K_{eq}\Delta \pmb a^{(n)}=\Delta \pmb Q^{(n)}   (n=0,1,2,...)
在载荷增量步长自动控制的求解方法中,\Delta \pmb Q^{(n)}可以表示成
\Delta \pmb Q^{(n)}=^{t+\Delta t} \pmb Q_l^{(n+1)}-^{t+\Delta t} \pmb Q_i^{(n)}    (7)
其中^{t+\Delta t} \pmb Q_l^{(n+1)}=^{t+\Delta t} \pmb Q_l^{(n)}+\Delta \pmb Q_l^{(n)}=(^{t+\Delta t}p^{(n)}+\Delta p^{(n)})\pmb Q_0
^{t+\Delta t} \pmb Q_i^{(n)}=\sum_{e} \int_{V_e}\pmb B^T {^{t+\Delta t}\pmb \sigma^{(n)}}dV
以上两式中:
^{t+\Delta t}p^{(0)}={^t}p       ^{t+\Delta t}\pmb \sigma^{(0)}={^t}\pmb \sigma
因此(7)可以改写为\Delta \pmb Q^{(n)}=\Delta \pmb Q_u^{(n)}+\Delta \pmb Q_l^{(n)}=\Delta \pmb Q_u^{(n)}+\Delta p^{(n)}\pmb Q_0
其中\Delta \pmb Q_u^{(n)}={^{t+\Delta t}}p^{(n)}\pmb Q_0-{^{t+\Delta t}} \pmb Q_i^{(n)}
\Delta p^{(n)}是在第n次迭代中由某个规定的约束条件来确定的载荷因子增量\Delta p的第n次修正量。在现在的方法中,这个条件就是某个结点的位移增量的大小,例如规定\Delta \pmb a中的最大分量\Delta a_g是给定的,此条件可以表示为
\Delta a_g^{(n)}=\pmb b^T \Delta \pmb a^{(n)}=\left\{\begin{matrix}\Delta l(n=0)\\0 (n=1,2,...)\end{matrix}\right.(8)
其中\Delta l\Delta a_g的规定值,\pmb b是除第g个元素为1,其余元素为零的向量,具体迭代算法步骤如下:
(1)计算对于节点载荷模式\pmb Q_0的位移模式\pmb a_0,即
\pmb a_0={^\tau} \pmb K_{eq}^{-1}\pmb Q_0
(2)计算对于不平衡结点力向量\Delta \pmb Q_u^{(n)}的位移增量修正值\Delta \pmb a_u^{(n)}n次迭代后位移增量修正值的全量\Delta \pmb a^{(n)},即
\Delta \pmb a_u^{(n)}={^\tau} \pmb K_{eq}^{-1}\Delta \pmb Q_u^{(n)}
\Delta \pmb a^{(n)}=\Delta \pmb a_u^{(n)}+\Delta p^{(n)}\pmb a_0
其中\Delta p^{(n)}是待定载荷因子增量的修正值
(3)利用条件(8)确定\Delta p^{(n)},由于
\Delta a_g^{(n)}=\pmb b^T \Delta \pmb a^{(n)}=\pmb b^T \Delta \pmb a_u^{(n)}+\Delta p^{(n)}\pmb b^T\pmb a_0=\left\{\begin{matrix}\Delta l(n=0)\\0 (n=1,2,...)\end{matrix}\right.(8)
从上式可以得到
\Delta p^{(n)}=\frac{\Delta a_g^{(n)}-\pmb b^T \Delta \pmb a_u^{(n)}}{\pmb b^T\pmb a_0}
这样就确定了第n次迭代后的载荷因子增量
载荷因子增量求出之后,可以知道每一步的外载荷,从而使用mN-R法迭代得到位移,进而求出应变,应力和当前增量步的内力等物理量。并检验收敛准则是否满足,如未满足,回到步骤二进行新的一次迭代,直至收敛准则满足为止。关于每一个增量步某个指定结点位移增量\Delta l本身的选择,通常的方法是第1个增量步可以由某个给定的载荷因子增量\Delta p_1(例如令\Delta p_1=p_e/N,其中p_e是弹性极限载荷因子,N可取5~10),通过求解得到\Delta l_1,以后各增量步的\Delta l_i可由下式确定,即
\Delta l_i=\sqrt{\frac{r_0}{r_{i-1}}}\Delta l_{i-1}
其中r_0=4,5,6

6 小结

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,001评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,210评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 161,874评论 0 351
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,001评论 1 291
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,022评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,005评论 1 295
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,929评论 3 416
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,742评论 0 271
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,193评论 1 309
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,427评论 2 331
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,583评论 1 346
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,305评论 5 342
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,911评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,564评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,731评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,581评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,478评论 2 352

推荐阅读更多精彩内容