非线性有限元及弹塑性力学讲解(哈工大 )

1、第一章非线性代数方程组的数值解法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，一般只能用数值方法求近似解答。其实质是，用一系列线性方程组的解去逼近所讨论非线性方程组的解。本章将简单介绍有限元分析中常见的各种求解非线性方程组的数值方法。1.1 直接迭代法当用某些方法将(a)=0改造成迭代格式(a) =P(a)-R=K(a) a-

2、R=00后，设一初始未知量a，则由它可得-101a= K(a)R01如果问题是收敛的，a将比a有所改善。如此反复迭代可得n+1n-1nn+1na= K(a)R a=a-ana=maxai当设范数为nnTn1/2a=(a)a或设范数为nn收敛条件则为aa011.1 直接迭代法如果考虑到每步迭代(an) =P(an)-R=K(an) an-R0n将(a)视为不平衡力（或失衡力）并作为衡量收敛的标准n (a)R01应指出的是，对单变量情况，如讲义图示，直接迭代实质是“割线”法1.2 牛顿法和修正牛顿法如果将非线性方程(a) =0在a附近展开，则(a) =(an)+ (a)n an+。=0或用求和约定

3、可写为nni(aj)=i(a)+i,j(a)aj+=0 (i,j=1,2,)n 又如果(a)n的逆存在，则a近似等于n-1na-(a)n(a)n(a)=(a)n，Pn=(则an-KT(an)-1 Pn，an+1=an+an如此逐步计算，即可得到非线性方程的解答，这就是牛顿-拉夫森法。jn如果将非线性方程(a) =0在a1.2 牛顿法和修正牛顿法牛顿法要每步都计算切线矩阵KT（也称刚度）并解线性方程组，虽精度高，但工作量也大此。外，在某些非线性问题(如理想塑性和软化塑性问题)中用牛顿法，迭代过程中切线矩阵可能是奇异的或病态的，为了克服这一现象，可有多种处理方法，其一是按下式来求nnn1ai=(i

4、,j(ak)+ij)j(ak)n其中的作用是改变切线矩阵KT的主对角元素，使奇异性或病态得到改善。更多的改进方法可参看沈聚敏钢筋混凝土有限元与板壳极限分析等。1.2 牛顿法和修正牛顿法如果在计算的每一步内，矩阵KT都用初始近0似解KT计算，在这种情况下，仅第一步迭代需要完全求解一个线性方程组，如果将KT0三角分解并存储起来，而以后各步迭代中采用公式a=(KT(a) (a)n则只需对上式右端项中的 (a)进行回代就行了。这种方法称为修正的牛顿法。为了提高修正牛顿法的收敛速度可采用某些过量修正技术。讲义上作了简要介绍，请大家自己看。7n01n1.3 拟牛顿法拟牛顿法的主要思想是：n+1首先设(KT

5、)可写成如下修正形式nn+1n(KT)=(KT)+(KT)式中(KT)n称为修正矩阵。n+1接着设(KT)必须满足如下所谓拟牛顿方程n+1n+1nn+1nK(aa)=(a)(a)由此可建立拟牛顿法迭代格式(略去了下标T)a=(K)j(a)nin1ijnkKij(ajn+1n+1aj)=i(ak)i(ak)nn+1in+1a=a+anniniKn+1=+Knn要用拟牛顿法，还需给出修正矩阵的计算。推导修正矩阵算式的思路是：设(KT)n=(un)(vn)Tnn(u)和(v)是秩1（或秩2，讲义为秩2）的列向量，将修正矩阵代入拟牛顿方程可得Tnnn nn(KT)+(u)(v)(a)=()Tnn假设(

6、v)(a)0，则有nnTn-1nnn(u)=(v)(a)()-(KT)(a)nnn如果取(v)=(a)，则当(a)(0)时(KT)n =(a)n)T(a)n-1()n-nnTn(KT)(a) (a )当(a)n=(0)时，迭代已收敛，(KT)n =(0)。讲义上的内容比这里说明的多，但基本思路是一样的。关于秩2的算法，请大家自己看。讲义上的塞尔曼公式可用逆矩阵定义验证。对秩1算法来说，实际使用的步骤为：001.设(a)求(KT)；-10010002.求(a)=-(KT)()；a=a+a3. 计算()0；a0=a1 。00004. 计算(KT)= ()-(KT)(a)0T0T0(a)/(a)(a

7、)；100015.计算(KT)=(KT)+(KT);(KT)= (KT).6.重复第2步，直到达到精度要求为止。1.4 增量方法求解非线性方程组的另一类方法是增量方法。使用这种方法需要知道“荷载”项(R)为零时问题0的解(a)。在实际问题中，(R)经常代表真实荷0 载，(a)代表结构位移。在问题的初始状态，它们均为零。这种从问题的初值开始，随着荷i载列阵(R)按增量形式逐渐增大，研究(a)的变化规律的方法，称为增量方法。当问题的性质与加载的历史有关时，如弹塑性问题，则必须采用增量方法。设“荷载”(R)在任一增量步的值为(，为荷载增量因子，(为标准荷载列阵,则非线性方程(a)= (0)可写为(a

8、,)=(P(a)-()=(0)若+时的解答为(a)+(a)，象牛顿法一样，将(a)+(a),+) 按Taylor级数展开，则可得引入切线矩阵且略去高阶小量后可改写为(a, )+,aa+, +=0a=KT(a, ) 1设荷载增量因子分别取如下值0=012M=1则荷载(R)可分成M级，第m级荷载为m()，其增量为(m+1-m)()=m()。由此可得-1(a)m=KT(a)m,m)m(a)m+1=(a)m+(a)m但是，这样做的每一步都将产生误差，结果使解答漂移。讲义上简单介绍了四种解决漂移的方法，下面仅对混合法作简单说明，其他方法请大家自行阅读。所谓混合法是指，在增量法每一增量步进行自修正的迭代计

9、算。其m增量步n次迭代的计算公式为在实际计算中，对于mM-1的各增量步的计算，可以只进行少许几次(例如3次)迭代，而对于m=M-1，即最后的一个荷载增量，需耍使用较多次迭代，以使近似解更接近于真解。1.5 增量弧长法用迭代法或增量法进行极限分析时，在极值点附近往往可能不收敛。这时可用增量弧长法来解决。增量弧长法的基本思想是：将作为独立变量，在每个增量步进行自修正法平衡迭代，在迭代过程中自动控制荷载因子的取值。也即nn1nnam=(KTm)(mPm)nPmnn=(am,m)an+1m=a+nmn为便于理解，以杆单向拉伸iii如图所示，矢径可表达为r=umia+mjRi+1i+1i+1r=umia

10、+mjR因为ii+1ii+1i+1r=rr=amia+mjR由于弧长法引入了如下约束方程ii+12rm=rm=lu=aam因此i+1i+2um=amami+1ii+1um=um+am由此可得imi+1m由矢量代数和约束方程可得iiii2i22r=rr=(um)+(m)=li+1iiii2r=(r+r)(r+r)=l因此iiir(2r+r)=0也即i+1i+1i+1i(amia+mjR)(am+2um)iaii+1+(2m+m)jR)=0若记 =(Ki1iTm)1 =(Ki+1n(mmi2iTm)P1im则aa=(=K)ni+1mmni+11Tm1+ in+1Pm2)将其代入约束方程，可得i+1

11、2i+1a(m)+2bm+c=0式中系数为i+1Ti+1a=1+(1)(1)ii+1Ti+1ib=m+(1)(2)+umc=()()+2ui上述式子是从简单情况推出的，如果除m外均理解为矩阵，即为一般情况的弧长法方程。ii一元二次方程有两个根，应取r和r间成锐角的根。据此可建立判别条件，具体推导这里从略了。i+1T2i+12im从上述公式可见，求的工作量是很大ii的，为此，可令r和r相互垂直，也即iirr=0这样做迭代的轨迹很接近圆弧，而计算工作量可减少很多。i+1m由相互垂直的条件可得iTi+1(um)2i+1m=iTi+1i(um)1+m综上所述，弧长法求解步骤为：11)选定荷载参考值，和本步荷载因子m，解12121T1得am，由l=(m)(1)(1)求弧长。2)修改切线刚度矩阵并三角化。检查对角元，正定则应行列式为i3)与上一步同时，求不平衡力Pm=(am,m).ii+1i+14)由和Pm求1和2。5)计算荷载