非线性有限元及弹塑性力学chap1-典尚设计

1、第一章 非线性代数方程组的数值解法1.1 直接迭代法1.2 牛顿法和修正牛顿法1.3 拟牛顿法1.4 增量方法1.5 增量弧长法2000.3哈尔滨工业大学 王焕定教授制作2 非线性问题可分为三类：材料非线性 不管那类非线性问题，最终都归结为一组非线性方程(a)=0，a为待求的未知量。 对许多问题，用某些方法可将(a)=0改造成(a) =P(a)-R=K(a) a -R=0的形式。 对非线性问题的方程(a)=0，一般只能用数值方法求近似解答。 、几何非线性 和边界非线性。 我们只讨论前两类问题。 其实质是，用一系列线性方程组的解去逼近所讨论非线性方程组的解。本章将简单介绍有限元分析中常见的各种求

2、解非线性方程组的数值方法。2000.3哈尔滨工业大学 王焕定教授制作31.1 直接迭代法 当用某些方法将(a)=0改造成迭代格式(a) =P(a)-R=K(a) a -R=0后a1= K(a0)-1R如果问题是收敛的， a1将比a0有所改善。an+1= K(an)-1R an=an+1- an当设范数为或设范数为收敛条件则为 ，设一初始未知量a0 ，则由它可得 如此反复迭代可得2000.3哈尔滨工业大学 王焕定教授制作4 如果考虑到每步迭代(an) =P(an)-R=K(an) an -R0将(an)视为不平衡力（或失衡力）并作为衡量收敛的标准 应指出的是，对单变量情况，如讲义图示，直接迭代实

3、质是“割线”法1.1 直接迭代法返首页 ，则收敛条件也可改为 ，一定条件下这种迭代过程是收敛的 ，但对多自由度情况，由于未知量通过矩阵K(an)的元素互相耦合，在迭代过程中往往出现不稳定现象。2000.3哈尔滨工业大学 王焕定教授制作51.2 牛顿法和修正牛顿法 如果将非线性方程(a) =0在an 附近展开，则又如果(a)n的逆存在，则an 近似等于 记 KT(an)=(a)n，Pn =(an)an-(a)n-1(an)则 an-KT(an)-1 Pn ， an+1=an+an切线矩阵不平衡力 如此逐步计算，即可得到非线性方程的解答，这就是牛顿-拉夫森法。(a) =(an)+ (a)n an+

4、。 =0或用求和约定可写为2000.3哈尔滨工业大学 王焕定教授制作61.2 牛顿法和修正牛顿法 牛顿法要每步都计算切线矩阵KT（也称刚度）并解线性方程组，虽精度高，但工作量也大。其中n的作用是改变切线矩阵KT的主对角元素，使奇异性或病态得到改善。更多的改进方法可参看沈聚敏钢筋混凝土有限元与板壳极限分析等。 此外，在某些非线性问题(如理想塑性和软化塑性问题)中用牛顿法，迭代过程中切线矩阵可能是奇异的或病态的 ，为了克服这一现象，可有多种处理方法，其一是按下式来求2000.3哈尔滨工业大学 王焕定教授制作71.2 牛顿法和修正牛顿法 如果在计算的每一步内，矩阵KT都用初始近似解KT0计算，在这种

5、情况下，仅第一步迭代需要完全求解一个线性方程组，如果将KT0三角分解并存储起来，而以后各步迭代中采用公式则只需对上式右端项中的 进行回代就行了。这种方法称为修正的牛顿法。 为了提高修正牛顿法的收敛速度可采用某些过量修正技术。讲义上作了简要介绍，请大家自己看。返首页2000.3哈尔滨工业大学 王焕定教授制作81.3 拟牛顿法拟牛顿法的主要思想是： 首先设(KT)n+1可写成如下修正形式 接着设(KT)n+1必须满足如下所谓拟牛顿方程(KT)n+1=(KT)n+(KT)n式中(KT)n称为修正矩阵。由此可建立拟牛顿法迭代格式(略去了下标T)2000.3哈尔滨工业大学 王焕定教授制作9 要用拟牛顿法

6、，还需给出修正矩阵的计算。推导修正矩阵算式的思路是：(un)和(vn)是秩1（或秩2，讲义为秩2）的列向量，将修正矩阵代入拟牛顿方程可得 设 (KT)n=(un)(vn)T如果取(vn)=(a)n，则当(a)n(0)时(KT)n +(un)(vn)T(a)n=()n(un)=(vn)T(a)n-1()n-(KT)n(a)n假设(vn)T(a)n0，则有(KT)n =(a)n)T(a)n-1()n-(KT)n(a)n (a n)T当(a)n=(0)时，迭代已收敛，(KT)n =(0)。2000.3哈尔滨工业大学 王焕定教授制作10 讲义上的内容比这里说明的多，但基本思路是一样的。关于秩2的算法，

7、请大家自己看。 1.设(a)0求(KT)0 ； 对秩1算法来说，实际使用的步骤为： 2.求 (a)0=-(KT)0-1()0；a1=a0+a0 5.计算(KT)1=(KT)0+(KT)0 ;(KT)0= (KT)1. 3. 计算()0； a0=a1 。 6.重复第2步，直到达到精度要求为止。 4. 计算 (KT)0 = ()0-(KT)0(a)0 (a0)T/(a0)T(a)0； 讲义上的塞尔曼公式可用逆矩阵定义验证。2000.3哈尔滨工业大学 王焕定教授制作11返首页讲义上给出了三种方法的对比，指出了选用什么算法应考虑的因素。2000.3哈尔滨工业大学 王焕定教授制作12 求解非线性方程组的

8、另一类方法是增量方法。 当问题的性质与加载的历史有关时，如弹塑性问题，则必须采用增量方法。1.4 增量方法使用这种方法需要知道“荷载”项(R)为零时问题的解(a)0。 在实际问题中，(R)经常代表真实荷载，(a)0 代表结构位移。在问题的初始状态，它们均为零。 这种从问题的初值开始，随着荷载列阵(R)按增量形式逐渐增大，研究(a)i的变化规律的方法，称为增量方法。2000.3哈尔滨工业大学 王焕定教授制作13 设“荷载”(R)在任一增量步的值为(R)， 为荷载增量因子，(R)为标准荷载列阵,则非线性方程 (a)= (0)可写为 引入切线矩阵且略去高阶小量后可改写为(a,)=(P(a)-(R)=

9、(0)若+时的解答为(a)+(a)，象牛顿法一样，将( (a)+(a),+) 按Taylor级数展开，则可得2000.3哈尔滨工业大学 王焕定教授制作14 设荷载增量因子分别取如下值(a)m+1=(a)m+(a)m则荷载(R)可分成M级，第m级荷载为m(R)，其增量为(m+1-m)(R)=m(R)。 由此可得(a)m=KT(a)m,m)-1m(R)但是，这样做的每一步都将产生误差，结果使解答漂移。讲义上简单介绍了四种解决漂移的方法，下面仅对混合法作简单说明，其他方法请大家自行阅读。2000.3哈尔滨工业大学 王焕定教授制作15 所谓混合法是指，在增量法每一增量步进行自修正的迭代计算。其m增量步

10、n次迭代的计算公式为 在实际计算中，对于 mM-1的各增量步的计算，可以只进行少许几次(例如3次)迭代，而对于m=M-1，即最后的一个荷载增量，需耍使用较多次迭代，以使近似解更接近于真解。返首页自修正不平衡力 用混合法求解时，所选取的荷载增量的步长可以比普通增量算法的步长大一些。2000.3哈尔滨工业大学 王焕定教授制作16 用迭代法或增量法进行极限分析时，在极值点附近往往可能不收敛。这时可用增量弧长法来解决。 为便于理解，以杆单向拉伸为例加以说明。1.5 增量弧长法 增量弧长法的基本思想是：将作为独立变量，在每个增量步进行自修正法平衡迭代，在迭代过程中自动控制荷载因子的取值。也即前步结果本步

11、n次增量2000.3哈尔滨工业大学 王焕定教授制作17 如图所示，矢径可表达为uuu因为 由于弧长法引入了如下约束方程因此由此可得2000.3哈尔滨工业大学 王焕定教授制作18 由矢量代数和约束方程可得因此若记也即则2000.3哈尔滨工业大学 王焕定教授制作19将其代入约束方程，可得式中系数为上述式子是从简单情况推出的，如果除 外均理解为矩阵，即为一般情况的弧长法方程。 一元二次方程有两个根，应取 和 间成锐角的根。据此可建立判别条件，具体推导这里从略了。2000.3哈尔滨工业大学 王焕定教授制作20这样做迭代的轨迹很接近圆弧，而计算工作量可减少很多。 从上述公式可见，求 的工作量是很大的，为此，可令 和 相互垂直，也即2000.3哈尔滨工业大学 王焕定教授制作21 由相互垂直的条件可得 综上所述，弧长法求解步骤为：1)选定荷载参考值 ，和本步荷载因子 ，解得 ，由 求弧长。2)修改切线刚度矩阵并三角化。检查对角元，正定则加载，负定则加负荷载。若矩阵对应行列式为零，达到极限荷载。4)由 和 求 和