1、有限元强度折减法1 背景1974 年，Smith & Hobbs1使用有限元方法分析了 u=0 条件下的边坡稳定性并与 Taylar2的结果进行对比，得到了很好的一致性；1975 年，Zienkiewicz等 3考虑 c、 进行有限元边坡稳定性分析，其结果与圆弧滑面解有较好吻合；1980 年 Griffiths4验证了一系列具有不同材料特性和形状的边坡稳定性并通过与 Bishop& Morgenstern5的结果进行了对比确定了数据的可靠性；此后也有研究证实了利用有限元方法进行边坡稳定性分析的可靠性 6,7,8,9；在文献 9中，引入一些案例证明了有限元强度折减法的准确性，并证明了有限元强度折
2、减法在分析非均质边坡时相对于传统方法的优越性。2001 年，郑颖人等 10把有限元强度折减法引入国内，并对此进行了后续研究 11,12,13,14。相较于一些传统的边坡稳定型分析方法，有限元强度折减法有以下几个优点 9：(1) 不必假设滑面的位置和形状，当土体自身强度不足以抵抗剪应力时土体失稳会自然发生。(2) 由于有限元强度折减法中没有条分的概念，因此也不必假设条间力，在整体失稳之前土体都处于整体稳定状态。(3) 使用有限元方法能够查看破坏过程。2 有限元强度系数折减法1.模型参数边坡模型主要包括六个参数，分别是：膨胀角 、内摩擦角 、黏聚力 c、弹性模量 E、泊松比 、重度 。膨胀角影响土
3、体屈服后的体积变形，若 0 则体积增大，=0 则体积不变。= 的情况被称之为关联流动法则，但是此时 值通常高于实验观测值，特别是在侧限条件下会提高土的承载力预测值。边坡稳定型问题通常是处于无侧限条件下，此时膨胀角的选取不再重要 9，因此文献9选取 =0 条件下的非关联流动法则，并且通过案例分析可以得出此膨胀角的选取可以得出准确的安全系数以及滑动面。c和 指 Mohr-Coulomb 准则中边坡土体的有效黏聚力和内摩擦角；E和是土体材料的弹性参数，这两个参数对土体稳定性分析的影响较小； 是土体的重度。应用有限元方法进行边坡稳定性分析中最重要的三个参数是 c、 、和 。2.屈服条件(1)Mohr-
4、Coulomb 准则Mohr-Coulomb 准则用大小主应力表示如式(1)所示：(1)132 =1+32 其中， 分别指土中一点的大小主应力。在主应力空间中，如果不考 1、 3虑 之间的大小关系，屈服面是一个不等角六棱锥，在 平面上是一1、 2、 3个等边不等角六边形。(4)(6)(2)D-P 准则D-P 准则可以写成式(2) 形式：(2)1+2=其中 I1 为第一应力不变量、J 2 为第二偏应力不变量， 和 kf为试验常数。在主应力空间中其屈服面为一个圆锥，在 平面上是一个圆形。3)D-P 准则转换为 Mohr-Coulomb 准则首先引入参数 b，如式(3)所示：(3)=2313则， 和
5、 分别可转化为式(4)：1 21=3(1+3)2 +(12)(13)2=3(13)其中 = 1+2将其带入(2)，得式 (5)：(5)13= 1.52+3(1+3)+ 2+3与式(1)对比可知两个准则之间的转换关系如式(6) 所示：= 1.52+32= 2+3因此，当 b=0 时，即外角点外接 DP 圆的两个试验常数分别如式(7)所示，当 b=1 时，即内角点外接 DP 圆的两个试验常数分别如式(8)所示。(7)=23(3)， = 63(3)(8)=23(3+)， = 63(3+)1 322121311 2b0外 角 点 外 接 DP圆内 角 点 外 接 圆3.安全系数的定义(1)Mohr-C
6、oulomb 准则中的安全系数(10)1955 年，Bishop 15首先在边坡稳定性分析中提出了抗剪强度折减的概念，在有限元强度折减法中通过将坡体的强度参数：黏聚力 c 和内摩擦角 同时除一个折减系数 Ft，得到一组新的 c和 值，作为一个新的强度参数输入进行试算，当计算不收敛时，对应的 Ft 即为所求的安全系数，此时坡体达到极限状态，发生剪切破坏。c=c/Ft=arctan(tan/Ft)(2)D-P(Drucker-Prager)准则中的安全系数取 Ft 为 D-P 准则中的强度折减系数，则 D-P 准则可以表示为式(9)，(9)1+2=(3)不同屈服条件下安全系数转换 13首先引入 M
7、ohr-Coulomb 等面积圆屈服准则，在 平面上，其屈服面是一个圆，并且面积与 Mohr-Coulomb 准则的不等角六边形相等， Mohr-Coulomb 等面积圆屈服准则中的试验参数如式(10)所示： = 3( 3)= 33式中=23+49224(23+1)(231)122(231) ， =(92)63简称外接圆屈服准则为 DP1 准则，其试验常数分别为 1，k f1；Mohr-Coulomb 等面积圆屈服准则为 DP2 准则，其试验常数分别为 2，k f2。 把 DP1 准则表示为 ，DP2 准则可表示为 。令1=2=11+1 2=2=21+2=12=kf1kf2=f()， ，所以1
8、=11+1=21+2。由此可知， 是 的函数，当 取不同值时可以得到12=21+221+2=()不同的 值如表 1 所列：表 1 不同内摩擦角时的 值4.失稳判据目前两个比较主流的失稳判据分别是有限元计算中力不平衡和位移的不收敛以及广义塑性应变或者等效塑性应变从坡脚到坡顶贯通。Griffiths 9和郑颖人11,12,13,14都使用计算不收敛作为失稳判据。Griffiths9提出，当在用户定义的最大迭代数目下计算仍不收敛时，则没有任何一种应力分布方式可以同时满足 Mohr-Coulomb 准则以及整体稳定，这种情况可看做边坡失稳判据。边坡失稳与数值计算不收敛同时发生，并伴随着极大的节点位移，
9、并以 1000 作为最大的迭代步数。郑颖人 14提出，有限元的计算迭代过程就是寻找外力和内力达到平衡状态的过程，整个迭代过程直到一个合适的收敛标准得到满足才停止。可见，如果边坡失稳破坏，滑面上将产生没有限制的塑性变形，有限元程序无法从有限元方程组中找到一个既能满足静力平衡又能满足应力-应变关系和强度准则的解，此时不管是从力的收敛标准，还是从位移的收敛标准来判断有限元计算都不收敛。3 案例分析例一，不含地基的均质边坡 9该边坡如图 1 所示，有限元程序采用 Mohr-Coulomb 失效准则，建立平面应变条件下八节点四边形单元减缩积分计算模型，其强度参数为 =20，c/h=0.05。边坡坡度为
10、26.57(2:1)，坡底水平，其边界条件为坡底约束竖直方向位移与水平方向位移，左侧约束水平方向位移，其余面为自由面。施加重力荷载后使安全系数从 0.8 到 1.4 逐步变化直至计算不收敛图 1 不含地基的均质边坡每一个安全系数对应的迭代次数如表 2 所列，当真正的安全系数接近时需要更多的迭代次数。表 2 例一计算结果当安全系数为 1.4 时，无量纲位移 Emax/H2 突变，并且此时计算无法收敛，在此情况下有限元计算结果与 Bishop & Morgenstern5给出的结果吻合良好，如图 2 所示。图 2 安全系数与无量纲位移边坡失稳时(FOS=1.4) 节点位移矢量和网格变形如图 3(a
11、)和图 3(b)所示，由此可得到边坡的潜在滑动面。图 3 安全系数为 1.4 计算不收敛时边坡变形(a)节点位移矢量 (b)网格变形例二，有软弱层的不排水黏性土边坡在本案例中，使用 Tresca 准则( u=0)进行总应力分析。边坡几何形状如图4 所示，地基厚度与边坡高度相同，该边坡有一个软弱层，在有限元计算中，令其抗剪强度(C u2)在一定范围内变化但其周围土体抗剪强度保持 Cu1/H=0.25 不变。利用有限元方法计算该边坡的安全系数结果如图 5 所示，对于均质边坡情况，C u2/Cu1=1，有限元计算结果与 Taylor2的结论很接近，随着软弱层的强度逐渐减小，在 Cu2/Cu10.6
12、时，结果发生了明显的变化。分别假定圆弧滑面和穿过软弱面的三段线滑面并利用 Janbu 法计算安全系数，可见在 Cu2/Cu10.6 处也发生了滑动机制的转换，当 Cu2/Cu10.6 时，潜在滑面形状为圆弧，当 Cu2/Cu10.6时，潜在滑面为结构软弱面。图 6 更加清晰的展示了这一现象，图 6(a)为均质边坡(C u2/Cu1=1)时的潜在滑面，可见此时的滑面形状为圆弧滑面，与 Taylor2的预测相同；图 6(c)为软弱层强度只有其周围土体 20%( Cu2/Cu1=0.2)时的潜在滑面，此时潜在滑面沿软弱层发展；图 6(b)为软弱层强度只有其周围土体 60%( Cu2/Cu1=0.6)
13、时的潜在滑面，此时圆弧滑面和沿软弱层的三段线式滑面都有可能发展，至少存在两种明显的滑动机制。图 4 有软弱层的不排水黏性边坡图 5 不同软弱层强度时的安全系数图 6 不同软弱层强度下的网格变形 (a) Cu2/Cu1=1.0 (b) Cu2/Cu1=0.6 (c) Cu2/Cu1=0.2例三，不同坡度边坡安全系数计算 13，验证 Mohr-Coulomb 等面积圆屈服准则均质边坡,坡高 H=20m,土容重 =25kN/m3,黏聚力 c=42kPa,内摩擦角 =17,求坡角 分别为 30,35,40,45,50时边坡的安全系数。计算结果如表 3 所列。表 3 安全系数计算结果从表中计算结果可以看
14、出,采用外接圆屈服准则计算的安全系数比传统的方法大许多,采用莫尔-库仑等面积圆屈服准则计算的结果与传统极限平衡方法(Spencer 法)计算的结果十分接近 ,说明采用莫尔- 库仑等面积圆屈服准则来代替莫尔- 库仑不等角六边形屈服准则是可行的,这样使计算大为方便。而采用外接圆屈服准则计算的安全系数要比莫尔-库仑等面积圆屈服准则计算的结果大 (1.21)倍。例四，存在两组节理面的岩质边坡稳定性分析 12如图 7 所示，岩体中存在两组方向不同的软弱结构面，贯通率 100%，第一组软弱结构面倾角为 30，平均间距 10m；第二组软弱结构面倾角 75，平均间距 10m。岩体重度为 25kN/m3，弹性模
15、量 11010Pa，泊松比 0.2，黏聚力图 Error! Main Document Only. 岩质边坡节理图 Error! Main Document Only.坡体破坏时的运动矢量图图 Error! Main Document Only. 极限状态时的塑性区1MPa，内摩擦角 38，两组节理参数相同，重度为 17kN/m3，弹性模量1107Pa，泊松比 0.3，黏聚力 0.12MPa，内摩擦角 24。按照二维平面应变问题建立有限元模型，按照连续介质处理。通过有限元强度折减，求得坡体破坏时的运动矢量如图 8 所示，滑动面如图 9(a)所示，它是最先贯通的塑性区，塑性区贯通并不等于破坏，当
16、塑性区贯通后继续发展到一定程度，岩体发生整体破坏，同时出现第二条贯通的塑性面，如图 9(b)所示。求得的稳定安全系数如表 4 所列，其中，极限平衡方法计算结果是根据最先贯通的那一条滑动面求得的。表 4 例 四计算结 果例 五，存在接触问题的边坡稳定性分析 12当 边坡中存在如图 10 所示的硬性结构面时，不能 按照例四中软弱结构面的方法进行处理，可以 采用接触单元来模拟硬性接触面的不连续性。 按照 Mohr-Coulomb 定律来定义接触面上的摩 擦行为，如式 11 所示，则其接触面上的安全系 数定义如式 11 所示。=c+tan，0 Ft=c/c=tan/tan (11)图 10 无充填的硬
17、性结构面图 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 面上施加有线性变化的面荷载，PM=400kPa，P C=0。图 11 折线型平面滑动岩质边坡在滑动面 AB，BM 上布置接触单元，坡体达到极限状态后的破坏滑动如图12 所示，并把有限元计算结果，与传统极限平衡方法 Spencer 法进行对比，接触单元的相关力学
18、参数以及两种计算结果对比如表 5 所列。图 12 坡体达到极限状态后的破坏滑动表 5 例五计算结果另外，在单元划分的过程中，在两个滑动面的交汇处形成了尖角，在尖角处形成较大的应力集中，求解时会产生病态方程。为了避免这些建模问题，需要在实体模型上，使用线的倒角来使尖角光滑化，或者在曲率突然变化的区域使用更细的网格。例六，泥岩层上粉质粘土边坡计算分析坡体材料力学参数为弹性模量 40MPa、泊松比 0.28、重度 19kN/m3、黏聚力 20kPa、内摩擦角 25，地基材料力学参数为弹性模量 400MPa、泊松比0.23、重度 24kN/m3、黏聚力 400kPa、内摩擦角 32。模型几何尺寸及边界
19、条件如图 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 sl
20、opes. Ge otechnique 24, No. 4, 531-559:2Taylor, D. W. (1937). Stability of earth slopes. J. Boston Soc. Civ. Eng. 24, 197-246:3Zienkiewicz, O. C., Humpheson, C. & Lewis, R. W.(1975). Associated and non-associated viscoplasticityand plasticity in soil mechanics. Geotechnique 25,671-689:4Griffiths, D.
21、 V. (1980). Finite element analyses of walls, footings and slopes. PhD thesis, University of Manchester.5Bishop, 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 ele
22、ments. 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
23、. Soils Found. 32, No. 1, 59-70:9Griffiths 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:
