非线性有限元及弹塑性力学讲解课件

第一章第一章 非线性代数方程组的数非线性代数方程组的数 值解法值解法 1.1 1.1 直接迭代法直接迭代法 1.2 1.2 牛顿法和修正牛顿法牛顿法和修正牛顿法 1.3 1.3 拟牛顿法拟牛顿法 1.4 1.4 增量方法增量方法 1.5 1.5 增量弧长法增量弧长法 2000.31哈尔滨工业大学 王焕定教授制作 非线性问题可分为三类：非线性问题可分为三类：材料非线性材料非线性 不管那类非线性问题，最终都归结为一组非不管那类非线性问题，最终都归结为一组非 线性方程线性方程( (a a)=0)=0，a a为待求的未知量。为待求的未知量。 对许多问题，用某些方法可将对许多问题，用某些方法可将( (a a)=0)=0改造成改造成 ( (a a) =) =P P( (a a) )- -R R= =K K( (a a) ) a a - -R R=0=0 的形式。的形式。 对非线性问题的方程对非线性问题的方程( (a a)=0)=0，一般只能用数一般只能用数 值方法求近似解答。值方法求近似解答。 、几何、几何 非线性非线性 和边界非线性。和边界非线性。 我们只讨论前两类问题。我们只讨论前两类问题。 其其实质是，用一系列线性实质是，用一系列线性 方程组的解去逼近所讨论非线性方程组的解。方程组的解去逼近所讨论非线性方程组的解。 本章将简单介绍有限元分析中常见的各种求解本章将简单介绍有限元分析中常见的各种求解 非线性方程组的数值方法。非线性方程组的数值方法。 2000.32哈尔滨工业大学 王焕定教授制作 1.1 1.1 直接迭代法直接迭代法 当用某些方法将当用某些方法将( (a a)=0)=0改造成迭代格式改造成迭代格式 ( (a a) =) =P P( (a a) )- -R R= =K K( (a a) ) a a - -R R=0=0 后后 a a1 1 = = K K( (a a 0 0 ) ) - -1 1R R 如果问题是收敛的，如果问题是收敛的， a a 1 1 将比将比a a 0 0 有所改善。有所改善。 a a n n+1+1= = K K( (a a n n ) ) - -1 1 R aR a n n =a=an n+1 +1 - - a a n n 当设范数为当设范数为 或设范数为或设范数为 收敛条件则为收敛条件则为 ，设一初始未知量，设一初始未知量a a 0 0 ，则由它可得则由它可得 如此反如此反 复迭代可得复迭代可得 2000.33哈尔滨工业大学 王焕定教授制作 如果考虑到每步迭代如果考虑到每步迭代 ( (a a n n ) =) =P P( (a a n n ) )- -R R= =K K( (a a n n ) ) a a n n - -R R00 将将( (a a n n ) )视为不平衡力（或失衡力）并作为衡量视为不平衡力（或失衡力）并作为衡量 收敛的标准收敛的标准 应指出的是，对单变量情况，如讲义图示，应指出的是，对单变量情况，如讲义图示， 直接迭代实质是直接迭代实质是“割线割线”法法 1.1 1.1 直接迭代法直接迭代法 返首页返首页 ，则收敛条件也可改为，则收敛条件也可改为 ，一定条件下这种一定条件下这种 迭代过程是收敛的迭代过程是收敛的 ，但对多自由度情况，由于，但对多自由度情况，由于 未知量通过矩阵未知量通过矩阵K K( (a a n n ) )的元素互相耦合，在迭代的元素互相耦合，在迭代 过程中往往出现不稳定现象。过程中往往出现不稳定现象。 2000.34哈尔滨工业大学 王焕定教授制作 1.2 1.2 牛顿法和修正牛顿法牛顿法和修正牛顿法 如果将非线性方程如果将非线性方程( (a a) =0) =0在在a a n n 附近展开，则附近展开，则 又如果又如果 ( (a a) ) n n 的逆存在，则的逆存在，则aan n 近似等于 近似等于 记记 K K T T( (a a n n )=)= ( (a a) ) n n ，P P n n = =( (a a n n ) ) aa n n - ( (a a) ) n n -1-1 ( (a a n n ) ) 则则 aa n n - -K K T T( (a a n n ) ) -1 -1 P Pn n ， a an n+1 +1 = =a a n n + +aa n n 切线矩阵切线矩阵 不平衡力不平衡力 如此逐步计算，即可得到非线性方程的解答如此逐步计算，即可得到非线性方程的解答 ，这就是牛顿，这就是牛顿- -拉夫森法。拉夫森法。 ( (a a) =) =( (a a n n )+ )+ ( (a a) )n n a a n n+ + 。 =0=0 或用求和约定可写为或用求和约定可写为 2000.35哈尔滨工业大学 王焕定教授制作 1.2 1.2 牛顿法和修正牛顿法牛顿法和修正牛顿法 牛顿法要每步都计算切线矩阵牛顿法要每步都计算切线矩阵K K T T （也称刚度也称刚度 ）并解线性方程组，虽精度高，但工作量也大）并解线性方程组，虽精度高，但工作量也大 。 其中其中 n n 的作用是改变切线矩阵的作用是改变切线矩阵K K T T 的主对角元素的主对角元素 ，使奇异性或病态得到改善。，使奇异性或病态得到改善。更多的改进方法更多的改进方法 可参看沈聚敏钢筋混凝土有限元与板壳极限可参看沈聚敏钢筋混凝土有限元与板壳极限 分析等。分析等。 此外，在某些非线性问题此外，在某些非线性问题( (如理想塑性和软化如理想塑性和软化 塑性问题塑性问题) )中用牛顿法，迭代过程中切线矩阵可中用牛顿法，迭代过程中切线矩阵可 能是奇异的或病态的能是奇异的或病态的 ，为了克服这一现象，可，为了克服这一现象，可 有多种处理方法，其一是按下式来求有多种处理方法，其一是按下式来求 2000.36哈尔滨工业大学 王焕定教授制作 1.2 1.2 牛顿法和修正牛顿法牛顿法和修正牛顿法 如果在计算的每一步内，矩阵如果在计算的每一步内，矩阵K K T T 都用初始近都用初始近 似解似解K K T T 0 0 计算，在这种情况下，仅第一步迭代需计算，在这种情况下，仅第一步迭代需 要完全求解一个线性方程组，如果将要完全求解一个线性方程组，如果将K K T T 0 0 三角分三角分 解并存储起来，而以后各步迭代中采用公式解并存储起来，而以后各步迭代中采用公式 则只需对上式右端项中的则只需对上式右端项中的 进行回代就行进行回代就行 了。这种方法称为修正的牛顿法。了。这种方法称为修正的牛顿法。 为了提高修正牛顿法的收敛速度可采用某些为了提高修正牛顿法的收敛速度可采用某些 过量修正技术。讲义上作了简要介绍，请大家过量修正技术。讲义上作了简要介绍，请大家 自己看。自己看。 返首页返首页 2000.37哈尔滨工业大学 王焕定教授制作 1.3 1.3 拟牛顿法拟牛顿法 拟牛顿法的主要思想是：拟牛顿法的主要思想是： 首先设首先设( (K K T T ) ) n n+1+1可写成如下修正形式 可写成如下修正形式 接着设接着设( (K K T T ) ) n n+1+1必须满足如下所谓拟牛顿方程 必须满足如下所谓拟牛顿方程 ( (K KT T ) ) n n+1+1=( =(K K T T) ) n n +(+( K KT T) )n n 式中式中( ( K KT T) )n n 称为修正矩阵。称为修正矩阵。 由此可建立拟牛顿法迭代格式由此可建立拟牛顿法迭代格式( (略去了下标略去了下标T)T) 2000.38哈尔滨工业大学 王焕定教授制作 要用拟牛顿法，还需给出修正矩阵的计算。要用拟牛顿法，还需给出修正矩阵的计算。 推导修正矩阵算式的思路是：推导修正矩阵算式的思路是： ( (u u n n ) )和和( (v v n n ) )是秩是秩1 1（或秩（或秩2 2，讲义为秩，讲义为秩2 2）的列向）的列向 量，将修正矩阵量，将修正矩阵代入代入拟牛顿方程可得拟牛顿方程可得 设设 ( ( K KT T) )n n = =( (u u n n )()(v v n n) )T T 如果取如果取( (v v n n )=()=( a a) ) n n ，则当则当( ( a a) ) n n (0)(0)时时 (K K T T) ) n n + +( (u u n n )()(v v n n) )T T ( ( a a) ) n n =(=( ) ) n n ( (u u n n ) )=( (v v n n) )T T( ( a a) ) n n -1-1 ( ( ) ) n n - -( (K K T T) ) n n( ( a a) ) n n 假设假设( (v v n n) )T T( ( a a) ) n n 00，则有则有 ( ( K KT T) ) n n = =( a a) ) n n) )T T( ( a a) ) n n -1-1 ( ( ) ) n n- - ( (K KT T ) )n n( ( a a) ) n n ( ( a a n n) )T T 当当( ( a a) ) n n =(0)=(0)时，迭代已收敛，时，迭代已收敛，( ( K KT T) ) n n =(0) =(0)。 2000.39哈尔滨工业大学 王焕定教授制作 讲义上的内容比这里说明的多，但基本思路讲义上的内容比这里说明的多，但基本思路 是一样的。关于秩是一样的。关于秩2 2的算法，请大家自己看。的算法，请大家自己看。 1. 1.设设( (a a) ) 0 0 求求( (K K T T) ) 0 0 ； 对秩对秩1 1算法来说，实际使用的步骤为：算法来说，实际使用的步骤为： 2. 2.求求 ( ( a a) ) 0 0 = =- -(K K T T) ) 0 0 -1-1 ( ( ) ) 0 0 ；a a 1 1 = =a a 0 0 + +aa 0 0 5. 5.计算计算( (K K T T) ) 1 1 = =( (K K T T) ) 0 0 + +( ( K KT T ) ) 0 0 ; ;( (K K T T) ) 0 0 = = ( (K K T T) ) 1 1 . . 3. 3. 计算计算( ( ) ) 0 0 ； a a 0 0 = =a a1 1 。 6.6.重复第重复第2 2步，直到达到精度要求为止。步，直到达到精度要求为止。 4. 4. 计算计算 ( ( K KT T ) ) 0 0 = (= ( ) ) 0 0 - -( (K K T T ) ) 0 0( ( a a) ) 0 0 ( ( a a0 0) )T T / /( ( a a0 0) )T T( ( a a) ) 0 0 ； 讲义上的塞尔曼公式可用逆矩阵定义验证。讲义上的塞尔曼公式可用逆矩阵定义验证。 2000.310哈尔滨工业大学 王焕定教授制作 返首页返首页 讲义上给出了讲义上给出了 三种方法的对比，三种方法的对比， 指出了选用什么算法指出了选用什么算法 应考虑的因素。应考虑的因素。 2000.311哈尔滨工业大学 王焕定教授制作 求解非线性方程组的另一类方法是增量方法求解非线性方程组的另一类方法是增量方法 。 当问题的性质与加载的历史有关时，如弹塑当问题的性质与加载的历史有关时，如弹塑 性问题，则必须采用增量方法。性问题，则必须采用增量方法。 1.4 1.4 增量方法增量方法 使用这种方法需要知道使用这种方法需要知道“荷载荷载”项项( (R R) )为零时问为零时问 题的解题的解( (a a) ) 0 0 。 在实际问题中，在实际问题中，( (R R) )经常代表真实荷经常代表真实荷 载，载，( (a a) )0 0 代表结构位移。在问题的初始状态， 代表结构位移。在问题的初始状态， 它们均为零。它们均为零。 这种从问题的初值开始，随着荷这种从问题的初值开始，随着荷 载列阵载列阵( (R R) )按增量形式逐渐增大，研究按增量形式逐渐增大，研究( (a a) ) i i 的变的变 化规律的方法，称为增量方法。化规律的方法，称为增量方法。 2000.312哈尔滨工业大学 王焕定教授制作 设设“荷载荷载”( (R R) )在任一增量步的值为在任一增量步的值为( (R R) )， 为荷载增量因子，为荷载增量因子，( (R R) )为标准荷载列阵为标准荷载列阵, ,则非则非 线性方程线性方程 ( ( (a a)= (0)= (0)可写为可写为 引入切线矩阵且略去高阶小量后可改写为引入切线矩阵且略去高阶小量后可改写为 ( ( (a, a,)=()=(P P( (a a) )- -( (R R)=(0)=(0) 若若 + + 时的解答为时的解答为( (a a)+)+( (a a) )，象牛顿法一象牛顿法一 样，将样，将( ( ( (a a) )+ +( (a a), , + + ) ) 按按TaylorTaylor级数级数 展开展开，则可得，则可得 2000.313哈尔滨工业大学 王焕定教授制作 设荷载增量因子设荷载增量因子分别取分别取如下值如下值 ( (a a) )m m+1+1=( =(a a) )m m + + ( (a a) )m m 则荷载则荷载( (R R) )可分成可分成MM级，第级，第mm级荷载为级荷载为 m m( (R R) )， ，其其 增量为增量为( ( m m+1+1 - - mm)( )(R R)= )=m m( (R R) )。 。 由此可得由此可得 ( (a a) )m m= =K K T T (a a) )m m , , m m) )-1 -1 m m( (R R) ) 但是，这样做的每一步都将产生误差，结果使但是，这样做的每一步都将产生误差，结果使 解答漂移。讲义上简单介绍了四种解决漂移的解答漂移。讲义上简单介绍了四种解决漂移的 方法，下面仅对混合法作简单说明，其他方法方法，下面仅对混合法作简单说明，其他方法 请大家自行阅读。请大家自行阅读。 2000.314哈尔滨工业大学 王焕定教授制作 所谓混合法是指，在增量法每一增量步进行所谓混合法是指，在增量法每一增量步进行 自修正的迭代计算。其自修正的迭代计算。其mm增量增量步步n n次迭代的计算次迭代的计算 公式为公式为 在实际计算中，对于在实际计算中，对于 mmMM-1-1的各增量步的的各增量步的 计算，可以只进行少许几次计算，可以只进行少许几次( (例如例如3 3次次) )迭代，而迭代，而 对于对于m=Mm=M - - 1 1，即最后的一个荷载增量，需耍使即最后的一个荷载增量，需耍使 用较多次迭代，以使近似解更接近于真解。用较多次迭代，以使近似解更接近于真解。 返首页返首页 自修正自修正 不平衡力不平衡力 用混合法求解时，所选取的荷载增量的步长用混合法求解时，所选取的荷载增量的步长 可以比普通增量算法的步长大一些。可以比普通增量算法的步长大一些。 2000.315哈尔滨工业大学 王焕定教授制作 用迭代法或增量法进行极限分析时，在极值用迭代法或增量法进行极限分析时，在极值 点附近往往可能不收敛。这时可用增量弧长法点附近往往可能不收敛。这时可用增量弧长法 来解决来解决。 为便于理解，以杆单向拉伸为例加以说明。为便于理解，以杆单向拉伸为例加以说明。 1.5 1.5 增量弧长法增量弧长法 增量弧长法的基本思想是：将增量弧长法的基本思想是：将 作为独立变量作为独立变量 ，在每个增量步进行自修正法平衡迭代，在迭，在每个增量步进行自修正法平衡迭代，在迭 代过程中自动控制荷载因子代过程中自动控制荷载因子 的取值。也即的取值。也即 前步结果前步结果 本步本步n n次增量次增量 2000.316哈尔滨工业大学 王焕定教授制作 如图所示，矢径可表达为如图所示，矢径可表达为 u u u 因为因为 由于弧长法引入了如下约束方程由于弧长法引入了如下约束方程 因此因此 由此可得由此可得 2000.317哈尔滨工业大学 王焕定教授制作 由矢量代数和约束方程可得由矢量代数和约束方程可得 因此因此 若记若记 也即也即 则则 2000.318哈尔滨工业大学 王焕定教授制作 将其代入约束方程，可得将其代入约束方程，可得 式中系数为式中系数为 上述式子是从简单情况推出的，如果除上述式子是从简单情况推出的，如果除 外外 均理解为矩阵，即为一般情况的弧长法方程。均理解为矩阵，即为一般情况的弧长法方程。 一元二次方程有两个根，应取一元二次方程有两个根，应取 和和 间成间成 锐角的根。据此可建立判别条件，具体推导这锐角的根。据此可建立判别条件，具体推导这 里从略了。里从略了。 2000.31