有限元强度折减法

有限元强度折减法
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′
2
sinφ′−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=λ√3
1−σ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 圆内角点外接D
3.安全系数的定义
(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 。

岩体重度为25kN/m 3，弹性模量1×1010Pa ，泊松比0.2，黏聚力1MPa ，内摩擦角38°，两组节理参数相同，重度为17kN/m 3，弹性模量1×107Pa ，泊松比0.3，黏聚力0.12MPa ，内摩擦角24°。

按照二维平面应变问题建立有限元模型，按照连续介质处理。

通过有限元强度折减，求得坡体破坏时的运动矢量如图8所示，滑动面如图9(a)所示，它是最先贯通的塑性区，塑性区贯通并不等于破坏，当塑性区贯通后继续发展到一定程度，岩体发生整体破坏，同时出现第二条贯通的塑性面，如图9(b)所示。

求得的稳定安全系数如表4所列，其中，极限平衡方法计算结果是根据最先贯通的那一条滑动面求得的。

图 7 岩质边坡节理
图 9坡体破坏时的运动矢量图 图 8 极限状态时的塑性区
例五，存在接触问题的边坡稳定性分析[12]
当边坡中存在如图10所示的硬性结构面时，不能按照例四中软弱结构面的方法进行处理，可以采用接触单元来模拟硬性接触面的不连续性。

按照Mohr-Coulomb定律来定义接触面上的摩擦行为，如式11所示，则其接触面上的安全系数定义如式11所示。

τ=c+σtanφ，σ≥0 F t=c/c’=tanφ/tanφ’ (11)
图10 无充填的硬性结构面
图11所示为两个直线滑面组成的折线型滑体ABMCD。

岩体重度γ=20kN/m3，弹性模量E=109Pa。

滑块ABCD面积433m2，滑面AB=20m，倾角为15°，AD=25m，DC=19.32m，BC=19.82m；滑块BCM面积196.5m2，滑面BM=28.03m，倾角为45°，CM=19.82m。

CM面上施加有线性变化的面荷载，P M=400kPa，P C=0。

图11 折线型平面滑动岩质边坡
在滑动面AB，BM上布置接触单元，坡体达到极限状态后的破坏滑动如图12所示，并把有限元计算结果，与传统极限平衡方法Spencer法进行对比，接触单元的相关力学参数以及两种计算结果对比如表5所列。

图12 坡体达到极限状态后的破坏滑动
另外，在单元划分的过程中，在两个滑动面的交汇处形成了尖角，在尖角处形成较大的应力集中，求解时会产生病态方程。

为了避免这些建模问题，需要在实体模型上，使用线的倒角来使尖角光滑化，或者在曲率突然变化的区域使用更细的网格。

例六，泥岩层上粉质粘土边坡计算分析
坡体材料力学参数为弹性模量40MPa、泊松比0.28、重度19kN/m3、黏聚力20kPa、内摩擦角25°，地基材料力学参数为弹性模量400MPa、泊松比0.23、重度24kN/m3、黏聚力400kPa、内摩擦角32°。

模型几何尺寸及边界条件如图13所示：
图13 模型几何尺寸及边界条件
通过ABAQUS有限元软件计算，当边坡的安全系数为1.3时，计算不收敛，通过Slide软件利用瑞典条分法计算得到的边坡安全稳定系数为1.295，瑞典条分法计算的滑面和有限元计算的塑性区如图14所示，可见两种方法计算出的安全系数和滑面吻合性较好。

5 参考文献
[1] Smith, I. M. & Hobbs, R. (1974). Finite element analysis of centrifuged and built-up slopes. Ge Âotechnique 24, No. 4, 531-559:
[2]Taylor, D. W. (1937). Stability of earth slopes. J. Boston Soc. Civ. Eng. 24, 197-246:
[3]Zienkiewicz, O. C., Humpheson, C. & Lewis, R. W.(1975). Associated and non-associated viscoplasticityand plasticity in soil mechanics. Geotechnique 25,671-689:
[4]Griffiths, D. V. (1980). Finite element analyses of walls, footings and slopes. PhD thesis, University of Manchester.
[5]Bishop, A. W. & Morgenstern, N. R. (1960). Stability coefficients for earth slopes. Geotechnique 10, 129-150:
[6] Griffiths, D. V. (1989). Computation of collapse loads in geomechanics by finite elements. Ing Arch 59, 237-244:
[7] Potts, D. M., Dounias, G. T. & Vaughan, P. R. (1990). Finite element analysis of progressive failure of Carsington embankment. Ge Âotechnique 40, No. 1, 79-102:
[8] Matsui, T. & San, K.-C. (1992). Finite element slope stability analysis by shear strength reduction technique. Soils Found. 32, No. 1, 59-70:
[9]Griffiths D V,Lane P A Slope stability analysis by finite elements [J] .Geotechnique, 1999, 49( 3): 387-403.
[10] 郑颖人等, 边坡稳定分析的一些进展地下空间, 2001(04)
[11] 郑颖人, 赵尚毅与张鲁渝, 用有限元强度折减法进行边坡稳定分析中国工程科学, 2002(10)
[12] 赵尚毅, 郑颖人与邓卫东, 用有限元强度折减法进行节理岩质边坡稳定性分析岩石力学与工程学报, 2003(02)
[13] 赵尚毅等, 用有限元强度折减法求边坡稳定安全系数岩土工程学报, 2002(03)
[14] 郑颖人与赵尚毅, 有限元强度折减法在土坡与岩坡中的应用. 岩石力学与工程学报, 2004(19).
[15] Bishop, A. W. (1955). The use of the slip circle in the stability analysis of slopes. Geotechnique 5, No. 1, 7-17:



















。