Arc-length Continuation Algorithm for Nonlinear Finite Element Equations

YIN Youquan1, DI Yuan1,†, YAO Zaixing2

1. College of Engineering, Peking University, Beijing 100871; 2. Beijing Lisuan Technology Co., Ltd, Beijing 100013; † Corresponding author, E-mail: diyuan@mech.pku.edu.cn

Abstract Stability analysis of engineering structures requires tracing equilibrium path of the structure when member’s buckling or material softening occurs. In nonlinear finite element analysis, the traditional Newton method fails at limit point and bifurcation point. The arc-length continuation method can overcome these numerical difficulties. To develop a nonlinear finite element code for stability analysis, the standard iteration formulation of Newton method is presented for the arc-length continuation method. Two practical formulations of the arc-length continuation method and their relationships with the standard form are also discussed. The applicability of the practical formulation is examined by the finite element analysis of stability for a slope. Key words arc-length continuation algorithm; equilibrium path curve; nonlinear analysis; Newton iteration method; finite element method

在结构非线性分析中, 结构平衡路径曲线上存在两种类型的临界点: 极值点和分叉点, 二者往往分别对应于结构的极值点型失稳和分叉点型失稳[1]。当运用有限元分析方法研究结构和连续体的非线性性质时, 应该追踪平衡路径, 进行增量计算,并识别极值点和分叉点, 这就需要使用延拓算法去求解非线性方程组[2]。

非线性有限元分析中的荷载增量法, 采用的是荷载控制的 Newton 方法, 其中包含延拓方法的思想。在平衡路径的极限点(力的转向点)附近, 切线

刚度矩阵趋于奇异, 计算分析的稳定性和收敛性相当差, 以至无法进行, 于是引入位移控制技术。在位移转向点附近, 位移控制方法亦失效, 为此 Sabir等[3]采用荷载控制和位移控制相互转换的控制技术。然而, 数值稳定、计算效率更高的方法还是

Riks[4]和 Wempner[5]建议的弧长法。该方法能够同时控制结构的变形与荷载, 使分析沿着结构的平衡路径顺利进行。经 Crisfield[6]和 Ramm[7]的改进, 弧长延拓方法已成为结构非线性分析中平衡路径追踪最主要的算法。

弧长延拓方法初期主要用于弹性结构的大变形和稳定性分析[8–12], 近年来开始应用于软化材料三维连续体稳定性的有限元分析[1,13–16]。弧长延拓算法的编程实现往往要根据原有非线性有限元程序的特点, 采用不同的 Newton 迭代格式。本文研究弧长延拓方法 Newton 迭代的标准格式和两种编程实现的实用格式, 并分析它们之间的关系。

1 有限元系统的非线性方程组

由虚功原理导出的有限元系统的平衡方程[17]如下:

ψ( σ )  P( σ )  R  0, (1) P( σ )  ct e B T σdv , (2) e

R  cT  NpT d V  ce T  N T qds 。 (3) e e Se其中, P( σ ) 和 R分别是物体应力矢量σ 和载荷矢量p , q 的等效节点力; N 和 ce分别是有限单元的形函数和选择矩阵; B是单元应变和节点位移之间的转换矩阵; ψ( σ )  0表示物体应力和外载处于平衡状态; ψ( σ )  0表示物体处于不平衡(失衡)状态,将 ψ ( σ )  R P ( σ )称为失衡力矢量。 设有限元系统的自由度总数为N , 则系统节点位移矢量a是N 维矢量: a [ a a , ..., an ]T。 (4) 1, 2等效节点力P( σ ) 和 R也是 N 维矢量, 如果将式(1)和(2)中应力 σ 用应变 ε 表示, 进而用位移表示, 可得到以位移a为变量的平衡方程: ψ( a )  P( a )  R  0。 (5)等效节点力 P()与应力  之间是线性关系。在线弹性材料情况下, P(a)与 a 之间也是线性关系。对于弹塑性材料, 应力是应变的非线性函数,也是位移的非线性函数。于是, P( a ) 是 a 的非线性函数, 也就是说, 式(5)是非线性方程组。

2 早期的延拓算法——载荷增量法

方程式(5)是以全量形式给出的, 通常需要用Newton 迭代算法进行数值求解。Newton 迭代具有

a 0局部收敛性的特征, 要求迭代给定的位移初值 离真解较近, 才能保证迭代解序列收敛到真解。当选择接近真解的初值难以做到时, 往往采用下面的延 拓方法[2]来处理。引入载荷参数 , 并嵌入式(5)中, 得到含参数 的平衡方程组: ψ( a ,  )  P( a )  R  0。(6)通常参数 是事先指定的, 并且逐渐增大, 将其离散为 0 0    ... m   ...   1, 1 2 m 1 M (7)     m。m m 1其中, m为施加的载荷步数, 为载荷参数的增量。随着  的增加, 逐步对每个增量  求解在m1下的节点位移am1, 称为求解式(6)的载荷增量法。在每一增量步内, 采用 Newton 迭代格式, 其初值取为前一增量步的解am。只要增量步足够小,就会满足局部收敛性条件。式(6)的 Jacobi 矩阵(即系统的切线刚度矩阵)和 Newton 迭代公式分别为ψ J   ce T  BDT B dvce , (8) a ep e

 a0  0,  (9) a n 1  an  J 1 ψ( a  an )。m其中, Dep是弹塑性本构矩阵。

对于 Jacobi 矩阵是非奇异(detj ≠ 0)的情况, 在增量步内迭代可以收敛到真解 am+1。当m= M, M =1 时, 就得到原方程组式(5)的解。这种引入载荷因子, 用增量方法求解非线性方程组的思想,就是延拓方法的思想, 载荷增量法是早期的(或经典的)延拓方法。

3 弧长延拓算法

对于某些强非线性问题, 例如金属结构大变形引起曲屈和突跳以及材料软化引起岩土结构的不稳定 [1,8,18], 载荷参数 不再是单调的, 式(7)不再成立。同时, 在载荷参数 转向点附近, 式(6)的Jacobi 矩阵 J 可能是奇异或病态的。这种情况下,采用载荷增量法求解式(6)将会失效。弧长延拓算法[4–5]不再将载荷参数  视为事前已知和单调增加的参数, 而是将它视为与位移分量同样待定的未知变量。于是, 式(6)成为含 N  1个未知变量( a 和 )的N个方程。为此需要引进一个辅助参数 s (通常是 N  1维空间解曲线的弧长), 并

增加一个与位移矢量 a、荷载参数 和弧长 s 关联的约束方程

B( a , λ , s )=0, (10)才能对问题求解。于是, 问题变为求解(N + 1)个变量的(N + 1) 阶方程组:

 ψ( a , ) 

ψ*( a ,  , s)  s)  0。 (11) B ( a ,  , 辅助参数s可选用N +1维空间中解曲线 a(s)和(s)的弧长。弧长微分的平方为(ds )2 d a Td a  (d )2 , (12)弧长参数s 是单调上升的, 将其离散化为 0= s s s   sm  sm   sm , 0 1 1 1 (13)  s  s s s m。m m 1然后, 按弧长增量 s 逐步地对方程组增量求解。弧长延拓算法的核心问题是构建具体的约束方程式(10)。约束方程可以采用各种不同的形式。Riks[19]和 Belytschko 等[20]在每个弧长步内采用的方程为

B( a ,  , s )  aT  a  (  )2 (  s)2  0, (14)其中, s是在位移荷载(N + 1)维空间中给定的弧长增量。式(14)是式(12)的有限增量表述形式。这时,扩展后的方程组式(11)的 Jacobi 矩阵为 J R  J*   。 (15)  T 2 a 2 J *显然, 矩阵 不总是正则的。在使用 Newton迭代算法时, 迭代初值总是取  a0  0, 0 0 ( a0 a , 0   ), 相应的 J *是奇异的。m 1 m m 1 m Crisfield[6] 给出对位移的控制要求ata= (s)2, 从而写出如式(10)的约束方程: ( a ,  , )  aT  a  (   0。B s s)2 (16)式(16)是一个球形路径控制的约束方程。为避免方程组式(11)的 Jacobi 矩阵失去对称性, 需要先导出求解  的二次代数方程, 求出后再求解a。在建立二次代数方程时, 方程式的系数含J 1。武际可等[8]通过位移–荷载参数(N+1)维空间内解曲线的单位切向矢量来建立约束方程。首先,构造一个(N+1)维矢量: v(a,  )  [ v 1, v 2,  , vn , vn ]T , (17) 1

 ψ ψˆ  ψ ψ , vi  ( 1) i det a ψ, , , , , ,    a a a  1 2 i N (18)

式中, “^”表示划去该列。矢量v( a ,  ) 是 N+1 维空间中解曲线 a(s)和(s)的切向矢量, 而解曲线的单位切向矢量为v τ( a ,  )  , (19) || v ||式中, ||v||代表矢量 v 的模。对第m个弧长步s, 约束方程为 a ( s) B( a,  , s )  τT( s )   s  0 , (20)  (s) 

式中, a  a ( s )  a ,  ( s )  , s  sm  sm。m m 1这时, 方程组式(11)的 Jacobi 矩阵可写为ψ ψ a   J  R J* B    B   τt  a   

 R1    R2    J  。 (21)    RN  1  N     N 1在临界点a a cr,   cr 附近, 矩阵 J 是奇异或*病态的, 但是矩阵J 是非奇异的(可证明 det J *= ||v||>0)。因此, 平衡路径上当遇到转向点时, 可以毫不困难地求解。本文根据式(20)给出的约束方程, 讨论弧长延拓算法 Newton 迭代的标准格式和有限元程序编程实现时可采用的实用格式。

4 迭代算法的标准格式

为了表述简洁紧凑, 本文采用新的变量: a y     at T   。 (22)

式(11)和 Jacobi 矩阵 J*可分别写为ψ ( y)  ψ*( y , s)   0, (23) B ( y , s)  * ψ J*( y , s)  。 (24) y

按弧长步增量求解。设s  sm时的解 ym已经得到, 在本增量步内需用 Newton 法求s  sm1时的解ym+1。设已经得一个近似解 ymn , n 表示迭代次数, 1 m1表示弧长增量步数。为了求得一个更好的解,将式(11)左端在 ymn 附近 Taylor 展开, 并截至线性1项, 可得* 0  ψ ( y, ) ψ * ( y n J * ( y n y  ymn s 1) 1)( 1) 。 (25) m  m 其中, y 的改进解为 ym , 迭代初值y 取自前一弧n+1 0 1 m 1长步的解 ym, 可得到 Newton 法标准的迭代公式:  y  ym , 0 m 1 (26)  n 1  1   y  yn ( J *( yn )) ψ *( yn )。m 1 m 1 m 1 m 1如果在式(26)中, 用矩阵J *( ym ) 代替 J *( yn ) , m 1就得到简化 Newton 法的迭代公式:  y0  y , m 1 m (27)  n  1  1  y  yn ( J *( y )) ψ *( yn )。m 1 m 1 m m 1如果仅迭代一次, 得到解y : 1 m 1 y1  y ( J *( y )) 1 ψ *( y ), (28) m 1 m m m即是自修正 Euler 方法的解[18]。在使用 Newton 法迭代(式(26))和简单 Newton法的迭代(式(27))时, 要根据式(20)和(11)正确地给J y n ) J ym ) ψ y n )出 *( , *( 和*( 的表达式。矩阵J* m 1 m 1的左上角子矩阵 J 是有限元系统的切线刚度矩阵,将其加边扩展为 J*( J *称为广义切线刚度矩阵)。ψ ymn )类似地, *(称为广义失衡力矢量, 是N + 1 1维矢量, 前 N 个分量构成的矢量是通常意义下的ψ( y n )失衡力 , 第 N +1个分量则为失约弧长m 1 B ( yn )。随着迭代步数 n 的增加, 广义失衡力矢m 1 ( y n 1  yn )  0, ym n+1 ym量趋于零,则收敛到真解 。m 1 m 1 1 1 Newton 法和简化 Newton 法的迭代公式也可以用增量形式给出。考虑到n 0 n n  y 1= y  y  y  y , m  m 1 m (29)  n 1 n 1 n 1 ,   y  y0  y  y  y m 1 m 1 m Newton 法增量形式的迭代公式为y 0=0, (30)  +1=   *) 1 *(  n n n  y y ( J ψ y y )。y yn m m简化 Newton 法增量形式的迭代公式为y 0= 0,  (31) n +1 n *  1 * n y = y ( J ) ψ ( y  y )。ym m 根据式(6), 失衡力矢量ψ 通常可写为ψ ( ym  yn )= P ( a  an )  (   n )R m m =P ( σm  σn )  ( m  n )R =P( σ ) P (  σ n )  R  n R m m =ψ ( ym )  ψ ( y n ), (32a)广义失衡力ψ *(y  y n)可写为m n ψ ( y  y ) * n m ψ( y  y )  ) m n  B(y  y  m  ψ( y ) ψ ( yn )  m (32b)  。n n τ T( y  y ) y  s m ψ(  y n ) n P (  σ )  nr   ct  BT σ dv  R。 (33) n n e e

式中, σ n是按弹塑性材料本构关系, 由 ε n 计算的应力增量矢量。采用简化Newton法时, 迭代的每一步都是按广义切线刚度 J*(ym)重新分配广义失衡力 ψ *(y m y n ) , 以便恢复平衡和满足约束条件。这种方法称为失衡力迭代方法, 更确切地说是广义失衡力迭代方法。参照式(27)和(29), 简化 Newton 法增量形式的迭代公式(30)可以改写为J*( y )( y n  1  yn )  ψ *( y  yn )  0。(34) m m其中, J *(ym)是 y  ym时的广义切线刚度矩阵, 它是( N +1) × ( N +1)的常数矩阵: J (a )  R  J *( ym )  m 。 (35)  τ T( ym) 其中, 左上角的子矩阵J ( am )是通常意义下的切线刚度矩阵, 为 N × N 的常数矩阵 , J ( am ) c T  B T DB d Vc , 其中积分项中的 DE 是在 y= e E e

e ym 情况下的弹塑性本构矩阵( Dep )m , 它在当前弧长步迭代过程中保持不变, 具有弹性本构矩阵的特点( σ  DE ε ), 但它又不是材料的弹性矩阵, 故称为准弹性矩阵。这时有σ n  D ε n  DB  a  Dbca, (36) E E E e E e J (a )  an  cT  B T  σ n d V  P ( σ n ), m e E E

en (37) P (  σ )   n R。  J * ( y ) y n  E  m T(

  τ y )yn  m