非线性有限元及弹塑性力学讲解,哈工大版---1

2000.3 哈尔滨工业大学 王焕定教授制作 9
讲义上的内容比这里说明的多，但基本思路 是一样的。关于秩2的算法，请大家自己看。 讲义上的塞尔曼公式可用逆矩阵定义验证。 对秩1算法来说，实际使用的步骤为： 1.设(a)0求(KT)0 ； 2.求 (Δa)0=-[(KT)0]-1(Ψ)0；a1=a0+Δa0 3. 计算(ΔΨ)0； a0=a1 。 4. 计算 (ΔKT)0 = [(ΔΨ)0-(KT)0(Δa)0] (Δa0)T/(Δa0)T(Δa)0； 5.计算(KT)1=(KT)0+(ΔKT)0 ;(KT)0= (KT)1. 6.重复第2步，直到达到精度要求为止。 2000.3 哈尔滨工业大学 王焕定教授制作
j
1.2 牛顿法和修正牛顿法
牛顿法要每步都计算切线矩阵KT（也称刚度） 并解线性方程组，虽精度高，但工作量也大。 此外，在某些非线性问题(如理想塑性和软化 塑性问题)中用牛顿法，迭代过程中切线矩阵可 能是奇异的或病态的 ，为了克服这一现象，可 有多种处理方法，其一是按下式来求 n n 1 n ai (Ψi , j (ak ) ij ) Ψj (ak ) 其中μn的作用是改变切线矩阵KT的主对角元素， 使奇异性或病态得到改善。更多的改进方法可 参看沈聚敏《钢筋混凝土有限元与板壳极限分 析》等。 2000.3 6 哈尔滨工业大学 王焕定教授制作
1
设荷载增量因子λ分别取如下值
0 0 1 2 M 1
则荷载(R)可分成M级，第m级荷载为λm(R)，其 增量为(λm+1-λm)(R)=Δλm(R)。
由此可得 Δ(a)m=[KT((a)m,λm)]-1Δλm(R) (a)m+1=(a)m+Δ(a)m 但是，这样做的每一步都将产生误差，结果使 解答漂移。讲义上简单介绍了四种解决漂移的 方法，下面仅对混合法作简单说明，其他方法 请大家自行阅读。 2000.3 14 哈尔滨工业大学 王焕定教授制作n Pmຫໍສະໝຸດ n Ψ(am ,2000.3
哈尔滨工业大学 王焕定教授制作
15
1.5 增量弧长法
用迭代法或增量法进行极限分析时，在极值 点附近往往可能不收敛。这时可用增量弧长法 来解决。 增量弧长法的基本思想是：将λ作为独立变量， 在每个增量步进行自修正法平衡迭代，在迭代 过程中自动控制荷载因子λ的取值。也即 n n 1 n n am ( K Tm ) ( R m Pm ) 前步结果
返首页
1.2 牛顿法和修正牛顿法
如果在计算的每一步内，矩阵KT都用初始近 似解KT0计算，在这种情况下，仅第一步迭代需 要完全求解一个线性方程组，如果将KT0三角分 解并存储起来，而以后各步迭代中采用公式
a ( K T (a )) Ψ(a ) n 则只需对上式右端项中的 Ψ(a ) 进行回代就行
n 0 1 n
了。这种方法称为修正的牛顿法。 为了提高修正牛顿法的收敛速度可采用某些 过量修正技术。讲义上作了简要介绍，请大家 自己看。
2000.3 哈尔滨工业大学 王焕定教授制作 7
1.3 拟牛顿法
拟牛顿法的主要思想是： 首先设(KT)n+1可写成如下修正形式 (KT)n+1=(KT)n+(ΔKT)n 式中(ΔKT)n称为修正矩阵。 接着设(KT)n+1必须满足如下所谓拟牛顿方程 K n1 (a n1 a n ) Ψ(a n1 ) Ψ(a n ) 由此可建立拟牛顿法迭代格式(略去了下标T)
2000.3
因此 i 1 i2 um am am u i i i u um1 um am1 u 由此可得 i i 1 i i r r r a m1 i a im 1 j R 由于弧长法引入了如下约束方程 i i 1 rm rm l 2
2000.3 哈尔滨工业大学 王焕定教授制作 18
将其代入约束方程，可得 i 1 2 i 1 a(m ) 2bm c 0 式中系数为 a 1 ( 1i 1 ) T ( 1i 1 ) i i 1 T i 1 i b m ( 1 ) [( 2 ) um ]
若λ+Δλ时的解答为(a)+Δ(a)，象牛顿法一 样，将(Ψ ((a)+Δ(a)),λ+Δλ) 按Taylor级数 展开，则可得
引入切线矩阵且略去高阶小量后可改写为
Ψ(a , ) Ψ,a Δ a Ψ, Δ 0
Δ a K T (a , ) R Δ
2000.3 哈尔滨工业大学 王焕定教授制作 13
哈尔滨工业大学 王焕定教授制作 17
u a
i m
i 1 m
am
由矢量代数和约束方程可得 i i i i r r r ( um ) 2 ( im ) 2 l 2
i 1 i i i i 2 r ( r r ) ( r r ) l
第一章 非线性代数方程组的数 值解法
1.1 直接迭代法 1.2 牛顿法和修正牛顿法 1.3 拟牛顿法 1.4 增量方法 1.5 增量弧长法
非线性问题可分为三类：材料非线性 、几何 非线性和边界非线性。我们只讨论前两类问题。 不管那类非线性问题，最终都归结为一组非 线性方程Ψ(a)=0，a为待求的未知量。 对许多问题，用某些方法可将Ψ(a)=0改造成 Ψ(a) =P(a)-R=K(a) a -R=0 的形式。 对非线性问题的方程Ψ(a)=0，一般只能用数 值方法求近似解答。其实质是，用一系列线性 方程组的解去逼近所讨论非线性方程组的解。 本章将简单介绍有限元分析中常见的各种求解 非线性方程组的数值方法。 2000.3 2 哈尔滨工业大学 王焕定教授制作
因此
i i i r (2r r ) 0 也即 i 1 i 1 i 1 i (am ia m j R ) [(am 2um )ia (2im im 1 ) j R )] 0 若记 i i i i 1 ( K Tm ) 1 R i2 ( K Tm ) 1 Pm 则 i 1 n in 1 n i 1 n i 1 1 aam ( Tm1) m Pm ) m K (R m 2
从上述公式可见，求 的工作量是很大的， i i 为此，可令 和 r 相互垂直，也即 r i i r r 0 这样做迭代的轨迹很接近圆弧，而计算工作量 可减少很多。
n n n am1 am am am u n 本步n次增量 为便于理解，以杆单向拉伸为例加以说明。
n Pm
n Ψ(am ,
n m)
2000.3
哈尔滨工业大学 王焕定教授制作
16
i i i 如图所示，矢径可表达为 r um i a m j R i 1 i 1 i 1 r um i a m j R 因为
i i i c ( 21 )T [( 21 ) 2um ] i 上述式子是从简单情况推出的，如果除 m 外 均理解为矩阵，即为一般情况的弧长法方程。 i i 一元二次方程有两个根，应取 r 和 r 间成 锐角的根。据此可建立判别条件，具体推导这 里从略了。 2000.3 19 哈尔滨工业大学 王焕定教授制作
返首页
所谓混合法是指，在增量法每一增量步进行 自修正的迭代计算。其m增量步n次迭代的计算 公式为 自修正 n n 1 n 不平衡力 am ( K Tm ) ( R m Pm )
n n n am1 am am m ) 1 在实际计算中，对于 m＜M-1的各增量步的 计算，可以只进行少许几次(例如3次)迭代， 而对于m=M-1，即最后的一个荷载增量，需耍 使用较多次迭代，以使近似解更接近于真解。 用混合法求解时，所选取的荷载增量的步长 可以比普通增量算法的步长大一些。
1.1 直接迭代法
当用某些方法将Ψ(a)=0改造成迭代格式 Ψ(a) =P(a)-R=K(a) a -R=0 后，设一初始未知量a0 ，则由它可得 a1= K(a0)-1R 如果问题是收敛的， a1将比a0有所改善。 如此反 复迭代可得 an+1= K(an)-1R Δan=an+1- an n a max a i 当设范数为
2000.3 哈尔滨工业大学 王焕定教授制作 4
1.2 牛顿法和修正牛顿法
如果将非线性方程Ψ(a) =0在an 附近展开，则 Ψ(a) =Ψ(an)+ [Ψ’(a)]n Δan+。。。 =0 或用求和约定可写为 n n Ψi (a j ) Ψi (a ) Ψi , j (a )a j 0 (i, j 1,2, ) 又如果[Ψ’(a)]n的逆存在，则Δan 近似等于 切线矩阵 Δan≈-[Ψ’(a)]n-1Ψ(an) 不平衡力 记 KT(an)=[Ψ’(a)]n，Pn =Ψ(an) 则 Δan≈-KT(an)-1 Pn ， an+1=an+Δan 如此逐步计算，即可得到非线性方程的解答， 这就是牛顿-拉夫森法。王焕定教授制作 2000.3 5 哈尔滨工业大学
n n ain (Kij )1Ψj (ak )
ain1 ain ain
n n
8
K ijn1 (an1 an ) Ψi (an1 ) Ψi (an ) j j k k
2000.3
K
n 1
哈尔滨工业大学 王焕定教授制作
K K
要用拟牛顿法，还需给出修正矩阵的计算。 推导修正矩阵算式的思路是： 设 (ΔKT)n=(un)(vn)T (un)和(vn)是秩1（或秩2，讲义为秩2）的列向 量，将修正矩阵代入拟牛顿方程可得 [(KT)n +(un)(vn)T](Δa)n=(ΔΨ)n 假设(vn)T(Δa)n≠0，则有 (un)=[(vn)T(Δa)n]-1[(ΔΨ)n-(KT)n(Δa)n] 如果取(vn)=(Δa)n，则当(Δa)n≠(0)时 (ΔKT)n =[((Δa)n)T(Δa)n]-1[(ΔΨ)n(KT)n(Δa)n] (Δa n)T 当(Δa)n=(0)时，迭代已收敛，(ΔKT)n =(0)。