有限元强度折减法

有限元强度折减法1 背景1974年，Smith & Hobbs[1]使用有限元方法分析了φu=0条件下的边坡稳定性并与Taylar[2]的结果进行对比，得到了很好的一致性；1975年，Zienkiewicz等[3]考虑c’、φ’进行有限元边坡稳定性分析，其结果与圆弧滑面解有较好吻合；1980年Griffiths[4]验证了一系列具有不同材料特性和形状的边坡稳定性并通过与Bishop& Morgenstern[5]的结果进行了对比确定了数据的可靠性；此后也有研究证实了利用有限元方法进行边坡稳定性分析的可靠性[6,7,8,9]；在文献[9]中，引入一些案例证明了有限元强度折减法的准确性，并证明了有限元强度折减法在分析非均质边坡时相对于传统方法的优越性。

2001年，郑颖人等[10]把有限元强度折减法引入国内，并对此进行了后续研究[11,12,13,14]。

相较于一些传统的边坡稳定型分析方法，有限元强度折减法有以下几个优点[9]：(1)不必假设滑面的位置和形状，当土体自身强度不足以抵抗剪应力时土体失稳会自然发生。

(2)由于有限元强度折减法中没有条分的概念，因此也不必假设条间力，在整体失稳之前土体都处于整体稳定状态。

(3)使用有限元方法能够查看破坏过程。

2 有限元强度系数折减法1.模型参数边坡模型主要包括六个参数，分别是：膨胀角ψ、内摩擦角φ’、黏聚力c’、弹性模量E’、泊松比υ’、重度γ。

膨胀角影响土体屈服后的体积变形，若ψ<0，则土体屈服后体积减小，若ψ>0则体积增大，ψ=0则体积不变。

ψ=φ的情况被称之为关联流动法则，但是此时ψ值通常高于实验观测值，特别是在侧限条件下会提高土的承载力预测值。

边坡稳定型问题通常是处于无侧限条件下，此时膨胀角的选取不再重要[9]，因此文献[9]选取ψ=0条件下的非关联流动法则，并且通过案例分析可以得出此膨胀角的选取可以得出准确的安全系数以及滑动面。

c’和φ’指Mohr-Coulomb准则中边坡土体的有效黏聚力和内摩擦角；E’和υ’是土体材料的弹性参数，这两个参数对土体稳定性分析的影响较小；γ是土体的重度。

应用有限元方法进行边坡稳定性分析中最重要的三个参数是c’、φ’、和γ。

2.屈服条件(1)Mohr-Coulomb准则Mohr-Coulomb准则用大小主应力表示如式(1)所示：σ1′−σ3′2=σ1′+σ3′2sinφ′−c cosφ′(1)其中， σ1′、σ3′分别指土中一点的大小主应力。

在主应力空间中，如果不考虑σ1、σ2、σ3之间的大小关系，屈服面是一个不等角六棱锥，在π平面上是一个等边不等角六边形。

(2)D-P准则(4) (6) D-P 准则可以写成式(2)形式：−βI 1+√J 2=k f (2)其中I 1为第一应力不变量、J 2为第二偏应力不变量，β和k f 为试验常数。

在主应力空间中其屈服面为一个圆锥，在π平面上是一个圆形。

3)D-P 准则转换为Mohr-Coulomb 准则首先引入参数b ，如式(3)所示： b =σ2−σ3σ1−σ3(3) 则，I 1和√J 2分别可转化为式(4)：I 1=3(σ1+σ3)2+(b −12)(σ1−σ3)√J 2=λ√31−σ3) 其中λ=√1−b +b 2将其带入(2)，得式(5)：σ1−σ3= 1.5β2−bβ+λ√3(σ1+σ3)+k f β2−bβ+λ√3 (5)与式(1)对比可知两个准则之间的转换关系如式(6)所示：sinφb = 1.5β2−bβ+λ√3 2c b cosφb =k f β2−bβ+λ√3 因此，当b=0时，即外角点外接DP 圆的两个试验常数分别如式(7)所示，当b=1时，即内角点外接DP 圆的两个试验常数分别如式(8)所示。

β=√3(3−sinφ)k f =√3(3−sinφ) (7) β=√3(3+sinφ)k f =√3(3+sinφ) (8)21σ2≥≥σ2σ1≥b P 圆内角点外接D3.安全系数的定义(1)Mohr-Coulomb 准则中的安全系数1955年，Bishop [15]首先在边坡稳定性分析中提出了抗剪强度折减的概念，在有限元强度折减法中通过将坡体的强度参数：黏聚力c 和内摩擦角φ同时除一个折减系数F t ，得到一组新的c’和φ’值，作为一个新的强度参数输入进行试算，(10) 当计算不收敛时，对应的F t 即为所求的安全系数，此时坡体达到极限状态，发生剪切破坏。

c’=c/F tφ’=arctan(tanφ/F t )(2)D-P(Drucker-Prager)准则中的安全系数取F t 为D-P 准则中的强度折减系数，则D-P 准则可以表示为式(9)，−βF t I 1+√J 2=k f F t (9) (3)不同屈服条件下安全系数转换[13]首先引入Mohr-Coulomb 等面积圆屈服准则，在π平面上，其屈服面是一个圆，并且面积与Mohr-Coulomb 准则的不等角六边形相等，Mohr-Coulomb 等面积圆屈服准则中的试验参数如式(10)所示： −β=k f 式中θδ=arcsin −23Asinφ+[49A 2sin 2φ−4(sin 2φ3+1)(A 23−1)]122(sin 2φ3−1)，A =√26√3简称外接圆屈服准则为DP1准则，其试验常数分别为β1，k f1；Mohr-Coulomb 等面积圆屈服准则为DP2准则，其试验常数分别为β2，k f2。

把DP1准则表示为f 1=√J 2=β1I 1+k f1，DP2准则可表示为f 2=√J 2=β2I 1+k f2。

令η=β1\β2=k f1\k f2=f(φ)，f 1=β1I 1+k f1=ηβ2I 1+ηk f2，所以f 1f 2=ηβ2I 1+ηk f2β2I 1+k f2=η=f(φ)。

由此可知，η是φ的函数，当φ取不同值时可以得到不同的η值如表1所列：表 1 不同内摩擦角时的η值4.失稳判据目前两个比较主流的失稳判据分别是有限元计算中力不平衡和位移的不收敛以及广义塑性应变或者等效塑性应变从坡脚到坡顶贯通。

Griffiths [9]和郑颖人[11,12,13,14]都使用计算不收敛作为失稳判据。

Griffiths [9]提出，当在用户定义的最大迭代数目下计算仍不收敛时，则没有任何一种应力分布方式可以同时满足Mohr-Coulomb 准则以及整体稳定，这种情况可看做边坡失稳判据。

边坡失稳与数值计算不收敛同时发生，并伴随着极大的节点位移，并以1000作为最大的迭代步数。

郑颖人[14]提出，有限元的计算迭代过程就是寻找外力和内力达到平衡状态的过程，整个迭代过程直到一个合适的收敛标准得到满足才停止。

可见，如果边坡失稳破坏，滑面上将产生没有限制的塑性变形，有限元程序无法从有限元方程组中找到一个既能满足静力平衡又能满足应力-应变关系和强度准则的解，此时不管是从力的收敛标准，还是从位移的收敛标准来判断有限元计算都不收敛。

3 案例分析例一，不含地基的均质边坡[9]该边坡如图1所示，有限元程序采用Mohr-Coulomb失效准则，建立平面应变条件下八节点四边形单元减缩积分计算模型，其强度参数为φ’=20°，c’/γh=0.05。

边坡坡度为26.57°(2:1)，坡底水平，其边界条件为坡底约束竖直方向位移与水平方向位移，左侧约束水平方向位移，其余面为自由面。

施加重力荷载后使安全系数从0.8到1.4逐步变化直至计算不收敛图 1 不含地基的均质边坡每一个安全系数对应的迭代次数如表2所列，当真正的安全系数接近时需要更多的迭代次数。

表 2 例一计算结果当安全系数为1.4时，无量纲位移E’δmax/γH2突变，并且此时计算无法收敛，在此情况下有限元计算结果与Bishop & Morgenstern[5]给出的结果吻合良好，如图2所示。

图 2 安全系数与无量纲位移边坡失稳时(FOS=1.4)节点位移矢量和网格变形如图3(a)和图3(b)所示，由此可得到边坡的潜在滑动面。

图3安全系数为1.4计算不收敛时边坡变形(a)节点位移矢量(b)网格变形例二，有软弱层的不排水黏性土边坡在本案例中，使用Tresca准则(φu=0)进行总应力分析。

边坡几何形状如图4所示，地基厚度与边坡高度相同，该边坡有一个软弱层，在有限元计算中，令其抗剪强度(C u2)在一定范围内变化但其周围土体抗剪强度保持C u1/γH=0.25不变。

利用有限元方法计算该边坡的安全系数结果如图5所示，对于均质边坡情况，C u2/C u1=1，有限元计算结果与Taylor[2]的结论很接近，随着软弱层的强度逐渐减小，在C u2/C u1≈0.6时，结果发生了明显的变化。

分别假定圆弧滑面和穿过软弱面的三段线滑面并利用Janbu法计算安全系数，可见在C u2/C u1≈0.6处也发生了滑动机制的转换，当C u2/C u1>0.6时，潜在滑面形状为圆弧，当C u2/C u1<0.6时，潜在滑面为结构软弱面。

图6更加清晰的展示了这一现象，图6(a)为均质边坡(C u2/C u1=1)时的潜在滑面，可见此时的滑面形状为圆弧滑面，与Taylor[2]的预测相同；图6(c)为软弱层强度只有其周围土体20%( C u2/C u1=0.2)时的潜在滑面，此时潜在滑面沿软弱层发展；图6(b)为软弱层强度只有其周围土体60%( C u2/C u1=0.6)时的潜在滑面，此时圆弧滑面和沿软弱层的三段线式滑面都有可能发展，至少存在两种明显的滑动机制。

图4有软弱层的不排水黏性边坡图 5 不同软弱层强度时的安全系数图 6 不同软弱层强度下的网格变形(a) C u2/C u1=1.0 (b) C u2/C u1=0.6 (c) C u2/C u1=0.2例三，不同坡度边坡安全系数计算[13]，验证Mohr-Coulomb等面积圆屈服准则均质边坡,坡高H=20m,土容重γ=25kN/m3,黏聚力c=42kPa,内摩擦角φ=17°,求坡角β分别为30°,35°,40°,45°,50°时边坡的安全系数。

计算结果如表3所列。

表 3 安全系数计算结果从表中计算结果可以看出,采用外接圆屈服准则计算的安全系数比传统的方法大许多,采用莫尔-库仑等面积圆屈服准则计算的结果与传统极限平衡方法(Spencer 法)计算的结果十分接近,说明采用莫尔-库仑等面积圆屈服准则来代替莫尔-库仑不等角六边形屈服准则是可行的,这样使计算大为方便。

而采用外接圆屈服准则计算的安全系数要比莫尔-库仑等面积圆屈服准则计算的结果大η(1.21)倍。

例四，存在两组节理面的岩质边坡稳定性分析[12]如图7所示，岩体中存在两组方向不同的软弱结构面，贯通率100%，第一组软弱结构面倾角为30°，平均间距10m ；第二组软弱结构面倾角75°，平均间距10m 。