有限元强度折减法

1、有限元强度折减法1背景1974年，Smith & Hobbs1 使用有限元方法分析了血=0条件下的边坡稳定性并与Taylar的结果进行对比，得到了很好的一致性；1975年，Zienkiewicz等3 考虑c' ©进行有限元边坡稳定性分析，其结果与圆弧滑面解有较好吻合；1980年Griffiths验证了一系列具有不同材料特性和形状的边坡稳定性并通过与 Bishop & Morge nsterN5的结果进行了对比确定了数据的可靠性；此后也有研究证 实了利用有限元方法进行边坡稳定性分析的可靠性Q"9】；在文献9中,引入一些案例证明了有限元强度折减法的准确性

2、，并证明了有限元强度折减法在分析非均 质边坡时相对于传统方法的优越性。2001年，郑颖人等10把有限元强度折减法 引入国内，并对此进行了后续研究11,12,13,14。相较于一些传统的边坡稳定型分析方法，有限元强度折减法有以下几个优点9:(1) 不必假设滑面的位置和形状，当土体自身强度不足以抵抗剪应力时土体 失稳会自然发生。(2) 由于有限元强度折减法中没有条分的概念，因此也不必假设条间力，在 整体失稳之前土体都处于整体稳定状态。(3) 使用有限元方法能够查看破坏过程。2有限元强度系数折减法1. 模型参数边坡模型主要包括六个参数，分别是：膨胀角职内摩擦角黏聚力c'弹性模量E'、

3、泊松比u'重度丫。膨胀角影响土体屈服后的体积变形，若书0则土体屈服后体积减小，若书0 则体积增大，书=0则体积不变。书=©的情况被称之为关联流动法则，但是此时书值通常高于实验观测值，特别是在侧限条件下会提高土的承载力预测值。边坡稳定型问题通常是处于无侧限条件下， 此时膨胀角的选取不再重要9，因此文献9 选取书=0条件下的非关联流动法则，并且通过案例分析可以得出此膨胀角的选 取可以得出准确的安全系数以及滑动面。c和©指Mohr-Coulomb准则中边坡土体的有效黏聚力和内摩擦角；E'和u、是土体材料的弹性参数，这两个参数对土体稳定性分析的影响较小；丫是土体的重

4、度。应用有限元方法进行边坡稳定性分析中最重要的三个参数是c' ©、和Y2. 屈服条件(1) Mohr-Coulomb 准则Mohr-Coulomb准则用大小主应力表示如式(1)所示:?-?32?+?'2? ?其中，？?、？?分别指土中一点的大小主应力。在主应力空间中，如果不考虑?、?、?之间的大小关系，屈服面是一个不等角六棱锥，在n平面上是一个 等边不等角六边形。(2) D-P准则?-?3?-?3(? - ?)2 + (?-?v?=(?-V3?)其中？?= V1- ? ?将其带入(2),得式：?- ?=1.5?(? + ?3?) + ?2-?+巧_-?+3与式(1)对

5、比可知两个准则之间的转换关系如式 所示：1.5sin ? = ?-? - ?二？ 2心 -?2c?qos (?= ?-2V3因此，当b=0时，即外角点外接DP圆的两个试验常数分别如式 所示，当 b=1时，即内角点外接DP圆的两个试验常数分别如式(8)所示。2?卩=v3(3-?) 2?卩=v3(3+?，?=?=6?v3(3-?)6?v3(3+?)(8)D-P准则可以写成式(2)形式：-B1 + Vv2 = ?(2)其中Ii为第一应力不变量、J2为第二偏应力不变量，B和kf为试验常数。在 主应力空间中其屈服面为一个圆锥，在n平面上是一个圆形3)D-P准则转换为Mohr-Coulomb准则 首先引入

6、参数b,如式(3)所示：b =贝?1和二?分别可转化为式(4):_ 3(? + ?=002> 001外角点外接D芦 T> 0203=0 b = 1=3b=0(,)2> 3102 > 100 T 03 内角点外接D P圆02*3. 安全系数的定义(1)Mohr-Coulomb准则中的安全系数1955年, Bishop15首先在边坡稳定性分析中提出了抗剪强度折减的概念， 在 有限元强度折减法中通过将坡体的强度参数：黏聚力 c和内摩擦角©同时除一 个折减系数Ft,得到一组新的c'和©值，作为一个新的强度参数输入进行试算，当计算不收敛时，对应的Ft即

7、为所求的安全系数，此时坡体达到极限状态，发生 剪切破坏。c' =c/F ©' =arctan(tant) © /F(2) D-P(Drucker-Prager)准则中的安全系数取Ft为D-P准则中的强度折减系数，则 D-P准则可以表示为式(9),?-?+ "?=鬲(3) 不同屈服条件下安全系数转换13首先引入Mohr-Coulomb等面积圆屈服准则，在n平面上，其屈服面是 圆，并且面积与 Mohr-Coulomb准则的不等角六边形相等，Mohr-Coulomb 积圆屈服准则中的试验参数如式(10)所示：sin © (9)等面3 =V3(乜

8、？?？?？?？)v3ccos ©(10)?=-v3?124 2 2?2.一？? +1)(1)式中?= ? 9? 3玄宁-1)简称外接圆屈服准则为DP1准则，其试验常数分别为 3,kf1 ； Mohr-Coulomb 等面积圆屈服准则为DP2准则，其试验常数分别为 血，kf2。把DP1准则表示为 f1 = V?= ?!?+ ?1， DP2 准则可表示为 f2 = V?= ?+ ?2 o 令 f1?+?2n =3 2=kf1kf2=f( ©,) f1 = ?+ ?1= ?+ ?，所以 W =选?+?2 = ?= ?(?)由此可知，n是©的函数，当©取不同值时

9、可以得到不同的 n值如表1所 列：表1不同内摩擦角时的 n值<p/C)】0203()40.ia)L 165L233.301L367<p/ f)5060708090L4281*48()L52Ih 5461,5554. 失稳判据目前两个比较主流的失稳判据分别是有限元计算中力不平衡和位移的不收 敛以及广义塑性应变或者等效塑性应变从坡脚到坡顶贯通。Griffiths9和郑颖人11121314都使用计算不收敛作为失稳判据。Griffiths9提出，当在用户定义的最大迭代数目下计算仍不收敛时， 则没有任 何一种应力分布方式可以同时满足 Mohr-Coulomb准则以及整体稳定，这种情况 可看做

10、边坡失稳判据。边坡失稳与数值计算不收敛同时发生，并伴随着极大的节 点位移，并以1000作为最大的迭代步数。郑颖人14提出，有限元的计算迭代过程就是寻找外力和内力达到平衡状态 的过程，整个迭代过程直到一个合适的收敛标准得到满足才停止。可见，如果边坡失稳破坏，滑面上将产生没有限制的塑性变形，有限元程序无法从有限元方程 组中找到一个既能满足静力平衡又能满足应力-应变关系和强度准则的解，此时 不管是从力的收敛标准，还是从位移的收敛标准来判断有限元计算都不收敛。3案例分析例一，不含地基的均质边坡9该边坡如图1所示，有限元程序采用Mohr-Coulomb失效准则，建立平面应变条件下八节点四边形单元减缩积分

11、计算模型，其强度参数为©' =2氏/ 丫 h=d05边坡坡度为26.57 ° (2:1)坡底水平，其边界条件为坡底约束竖直方向位移与水平 方向位移，左侧约束水平方向位移，其余面为自由面。施加重力荷载后使安全系 数从0.8到1.4逐步变化直至计算不收敛1 仲-一一一+冷巴詈常图1不含地基的均质边坡每一个安全系数对应的迭代次数如表 2所列，当真正的安全系数接近时需要 更多的迭代次数。表2例一计算结果FOSE rHUK/y h'Iterationso-8o n037921 000 3KIJO1-20042220b300453411-350 5447921-401

12、476lOOt)当安全系数为1.4时，无量纲位移E' nSax/ yH突变，并且此时计算无法收敛， 在此情况下有限元计算结果与Bishop & MorgensterH5给出的结果吻合良好，女口图2所示。Bishop & Morgertstern (I960)FOS = 1 380oa1V2?4V6FOS图2安全系数与无量纲位移边坡失稳时(F0S=1.4)节点位移矢量和网格变形如图3(a)和图3(b)所示，由此可得到边坡的潜在滑动面。(9)图3安全系数为1.4计算不收敛时边坡变形 (a)节点位移矢量(b)网格变形例二，有软弱层的不排水黏性土边坡在本案例中，使用Tresca

13、准则(U=0)进行总应力分析。边坡几何形状如图4所示，地基厚度与边坡高度相同，该边坡有一个软弱层，在有限元计算中，令其 抗剪强度(Cu2)在一定范围内变化但其周围土体抗剪强度保持Cui/ 丫 H=0.2不变。利用有限元方法计算该边坡的安全系数结果如图5所示，对于均质边坡情况，Cu2/Cu1=1，有限元计算结果与Taylor的结论很接近，随着软弱层的强度逐 渐减小，在Cu2/Cu1- 0.6寸，结果发生了明显的变化。分别假定圆弧滑面和穿过软 弱面的三段线滑面并利用Janbu法计算安全系数，可见在 Cu2/Cu1 - 0.6处也发生 了滑动机制的转换，当Cu2/Cu1>0.6时，潜在滑面形状

14、为圆弧，当Cu2/Cu1<0.6时， 潜在滑面为结构软弱面。图6更加清晰的展示了这一现象，图6(a)为均质边坡(Cu2/Cu1=1)时的潜在滑面，可见此时的滑面形状为圆弧滑面，与Taylor2的预测相同；图6(c)为软弱层强度只有其周围土体 20%( Cu2/Cu1=0.2)时的潜在滑面，此 时潜在滑面沿软弱层发展；图6(b)为软弱层强度只有其周围土体60%( Cu2/Cu1=0.6) 时的潜在滑面，此时圆弧滑面和沿软弱层的三段线式滑面都有可能发展，至少存在两种明显的滑动机制。1 p3 -1614121 -0 80硏04 -02Taylor (1937)FOS = 1 -47唧YH =

15、02500图5不同软弱层强度时的安全系数图例三，不同坡度边坡安全系数计算13，验证Mohr-Coulomb等面积圆屈服准均质边坡，坡高H=20m,土容重丫 =25kN/m,黏聚力c=42kPa，内摩擦角 © =17求坡角B分别为30° ,35 ° ,40 ° ,45时边坡的安全系数。计算结果如表 3所列。表3安全系数计算结果坡角P /(°)-安全系数iw Pr 法|H K UlMho卩仕301.7SL47l.妙1.463351.62L341. 318402L22L 1531.212451.36L 121.0621. 115501291.060.

16、992LO38伦采川蚓按風風肘准叭-2采川克尔一压仑等面积關W鴨准Wl从表中计算结果可以看出，采用外接圆屈服准则计算的安全系数比传统的方 法大许多，采用莫尔-库仑等面积圆屈服准则计算的结果与传统极限平衡方法(Spenee法)计算的结果十分接近，说明采用莫尔-库仑等面积圆屈服准则来代替莫 尔-库仑不等角六边形屈服准则是可行的，这样使计算大为方便。而采用外接圆屈 服准则计算的安全系数要比莫尔-库仑等面积圆屈服准则计算的结果大n (21)倍。例四，存在两组节理面的岩质边坡稳定性分析12如图7所示，岩体中存在两组方向不同的软弱结构面，贯通率 100%，第一 组软弱结构面倾角为30°平均间距1

17、0m;第二组软弱结构面倾角75°平均间 距10m。岩体重度为25kN/m3,弹性模量1xi00Pa,泊松比0.2,黏聚力1MPa, 内摩擦角38°两组节理参数相同，重度为17kN/m3,弹性模量1X 10Pa,泊松比 0.3,黏聚力0.12MPa,内摩擦角24°按照二维平面应变问题建立有限元模型，按照连续介质处理。通过有限元强 度折减，求得坡体破坏时的运动矢量如图8所示，滑动面如图9(a)所示，它是最先贯通的塑性区，塑性区贯通并不等于破坏，当塑性区贯通后继续发展到一定程 度，岩体发生整体破坏，同时出现第二条贯通的塑性面，如图9(b)所示。求得的稳定安全系数如表4所

18、列，其中，极限平衡方法计算结果是根据最先贯通的那一 条滑动面求得的。<0图7岩质边坡节理图9坡体破坏时的运动矢量图图8极限状态时的塑性区表4例四计算结果计算方袪安全爲数相限尤法（外接倜屈眼准则）1.62有限无法（莫尔库仑等面积阀尿眼准则）1.33极限乎衛庁法朽pgnwr法）136例五，存在接触问题的边坡稳定性分析12当边坡中存在如图10所示的硬性结构面时，不能按照例四中软弱结构面的 方法进行处理，可以采用接触单元来模拟硬性接触面的不连续性。按照Mohr-Coulomb 定律来定义接触面上的摩擦行为，如式11所示，则其接触面上的安全系数定义如式11所示。T=tan © /tan(

19、11)图11所示为两个直线滑面组成的折线型滑体 ABMCD。岩体重度丫 =20kN/m, 弹性模量E=109Pa。滑块ABCD面积433m2，滑面AB=20m，倾角为15° ,AD=25m， DC=19.32m, BC=19.82m;滑块 BCM 面积 196.5m2，滑面 BM=28.03m，倾角为Pc=0。45°, CM=19.82m。在滑动面AB，BM上布置接触单元，坡体达到极限状态后的破坏滑动如图12所示，并把有限元计算结果，与传统极限平衡方法 Spencer法进行对比，接触单元的相关力学参数以及两种计算结果对比如表5所列图12坡体达到极限状态后的破坏滑动表5例五计

20、算结果有限元强度折减法Spencer 法1.00- I&0 kPa* y>= 30'2.11v 320 kPa*30工332.33£ = 1 bO kPa* 於=斗523Er = 0 kPa,= 453.0K2.4另外，在单元划分的过程中，在两个滑动面的交汇处形成了尖角， 在尖角处 形成较大的应力集中，求解时会产生病态方程。为了避免这些建模问题，需要在 实体模型上，使用线的倒角来使尖角光滑化，或者在曲率突然变化的区域使用更 细的网格。例六，泥岩层上粉质粘土边坡计算分析坡体材料力学参数为弹性模量 40MPa、泊松比0.28、重度19kN/m3、黏聚力 20kPa、

21、内摩擦角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). Fini

22、te element analysis of centrifuged and built-up slopes. Ge ?otech ni que 24, No. 4, 531-559:2 Taylor, D. W. (1937). Stability of earth slopes. J. Boston Soc. Civ. Eng. 24, 197246:3 Zie nkiewicz, O. C., Humpheso n, C. & Lewis, R. W.(1975). Associated and nonassociated viscoplasticityand plasticit

23、y in soil mechanics. Geotechnique 25,671-689:4 Griffiths, D. V . (1980). Finite element analyses of walls, footings and slopes. PhD thesis, Uni versity of Man chester.5 Bishop, A. W. & Morgenstern, N. R. (1960). Stability coefficients for earth slopes. Geotech nique 10, 129-150:6 Griffiths, D. V. (1989). Computation of collapse loads in geomechanicsby finite elements. I ng Arch 5