不更新切向量, 减少了大量高阶矩阵行列式的计算,能够减小一个弧长步内的计算量, 同时也能保证实用格式的收敛性。参考文献

式中, ym是弧长s=sm 的解, s  sm  sm是本弧长1步的弧长增量。2) 迭代修正。  y 1 E y , n (44)   ψ( ym  ∆y ) n 1 n   y  y ( J *( y m )) 1 ,   0 式中, ψ ( y  y n ) 是第 n 次迭代步的失衡力矢m量。失约弧长项取值为 0, 表明

τT  ye τ T  yn τ T  y n 1 τ T  y N  s。 (45) m m m m式 (45)相当于设定约束条件为 B(y, s)= τ T y m  s  0。在每一迭代步, 将近似解增量y n投影到τm 方向, 可得到弧长增量。本文算法采用初应力迭代格式, 将两步合并完成:   y 0  0,   ψ ( y )  Pn  (46) y n 1 ( J)* 1 m m  , n 0, 1, 2,,   s 式中, P n  P(  σ n  σ n )  c T BT( σ n  σn )dA ; E e E e

n ( εn σ 是按照材料本构关系  Dep d ), 由计算得到的应力增量; σ n 是按照准弹性关系E ( σ  D ε ( D )m ε), 由 ε计算出的应力增量; E ep Pn为初应力等效节点力矢量。初值P 取为零, 得0 1到的第一次近似解就是 Euler 预报解, y = y E。当失约弧长项取为s 时, 有τ T y 1  τ T yn τ T y n 1 m m m τ T  y N  s , 相当于约束条件为m τT  y  s  0。 (47) m上述两种实用格式分别是简化 Newton 法的失衡力迭代格式和初应力迭代格式, 且将约束方程简化为τT  y  s  0。在迭代过程中, 这两种实用格m τ τ( y  y n ),式用 代替 虽然略去高阶小量, 但在m m弧长步内不再更新切向量, 从而简化了算法, 利于编程实现, 也能够减小计算量。这两种实用格式采用更简单的约束方程式(47)来替代式(20), 因而迭代的收敛性是有保障的。

6 弧长法在边坡稳定分析中的应用

[16]方建瑞等 将弧长算法应用于边坡的稳定性研究, 边坡材料使用 Duncan-chang 模型, 这是一 种稳定材料模型, 因而其边坡失稳研究采用的是失衡机制。

基于失稳机制研究边坡的稳定性, 需要建立平衡路径曲线。对于边坡这种多自由度有限元系统,需要引入广义力和广义位移来描述平衡路径。

设边坡的等效节点荷载为 R, 引入荷载参数,荷载做功的变分为

δw  R δ a   δ( R  a)。 (48)如果定义广义力为F   , 则相应广义位移为U  R a。平衡路径曲线F F ( U )能直观地分析边坡有限元系统的稳定性。

基于初应力迭代格式的弧长法, 编制有限元分析程序, 计算某边坡失稳机制下的安全储备系数(安全系数)Kcr。该边坡有限元网格见图 2, 边坡高度为 9.5 m, 坡角为 64 , 边坡仅受重力荷载作用,岩土重度为 24 kn/m3, 弹性模量为 28.7 GPA, 泊松比为 0.27。单元选用平面应变的三角形网格, 两侧边界上节点的水平位移被约束, 底面节点固定。

采用强度折减法分析该边坡的稳定性, 强度折减系数设为。边坡材料服从 Drucker-prager 屈服准则: f  J  ( ) I k (  )  0, (49) 2 1式中, I1 和 J1分别是应力张量的第一不变量和偏应力张量的第二不变量; 材料强度参数k ( )和 ( )是塑性内变量 的函数, 其峰值分别为 k = 61.50 kpa和 = 0.066。选用等效塑性应变为塑性内变量, 边坡材料的应变软化可通过下式的负幂次形式来考虑: k (  )  (  ) 1  (3e 2 1) 。 (50) k  4选取强度折减系数 1.00 , 弧长增量s  0.1 ,采用初应力格式的弧长法有限元程序进行增量分析。当荷载参数  1.44 时, 边坡达到临界状态,

此时边坡部分的变形如图 3 所示(边坡坡顶 A 点的水平位移为0.01 mm, 竖向位移为 0.20 mm), Mises应力和静水应力分布如图 4 和 5 所示。分别取折减

系数 = 1.25, 1.45, 1.50, 依次使用弧长延拓算法,计算得到相应的平衡路径曲线, 如图 6 所示。曲线上极值点(即临界点)的纵坐标即是稳定性临界点的载荷参数cr。当荷载参数 =1 时, 对应于边坡上作用有全部重力载荷的情况。图 6 中,  = 1.45 对应平衡路径曲线极值点的纵坐标cr  1, 即全部重力荷载作用下、材料强度折减系数 = 1.45 时, 该边坡处于稳定状态的临界点, 即该边坡的安全系数Kcr=1.45。

7 结语

材料软化和构件屈曲往往引起工程结构的不稳定, 对这种稳定性问题的研究需要追踪结构的平衡路径。当采用有限元方法进行数值计算时, 问题最终都归结为一非线性代数方程的求解。如果直接采用 Newton 迭代法求解非线性有限元方程组, 在平衡路径曲线上的临界点附近, 由于切线刚度矩阵趋于奇异而无法得到收敛解。弧长延拓算法能够很好地解决这一问题。在非线性有限元程序中编程实现弧长延拓算法时, 需要较为详细的迭代格式。本文给出弧长延拓算法 Newton 迭代和简化 Newton 迭代的标准格式,并讨论了失衡力迭代和初应力迭代两种实用格式及其与标准格式之间的关系。如果弧长增量步是一阶小量, 则这两种实用格式与标准格式在约束方程上的差是二阶小量。由于这两种实用格式在弧长步内

图 2 边坡算例的有限元网格Fig. 2 Finite element mesh of the slope

图 6 不同强度折减系数对应的平衡路径曲线Fig. 6 Equilibrium path curves for different strength reduction factors 

图 5边坡静水应力的分布Fig. 5 Hydrostatic stress distribution of the slope

图 4 边坡 Mises 应力的分布Fig. 4 Von Mises stress distribution of the slope

图 3 边坡的变形Fig. 3 Deformed mesh of the slope