有限元强度折减法

1、有限元强度折减法1背景1974年，Smith & Hobbs1 使用有限元方法分析了 血=0条件下的边坡稳定性 并与Taylar2的结果进行比照，得到了很好的一致性；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'

3、弹性模量E'、泊松比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?

5、(? + ?3?) + ?2-?+巧_-?+3与式(1)比照可知两个准那么之间的转换关系如式 所示：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平面上是一个

6、圆形3)D-P准那么转换为Mohr-Coulomb准那么 首先引入参数b,如式(3)所示：b =贝?1和二?分别可转化为式(4):_ 3(? + ?=0?02> 001外角点外接D芦 T> 0203=0 b = 1=3b=0(,)2> 310?2 > 10?0 T 03 内角点外接D P圆02*3. 平安系数的定义(1)Mohr-Coulomb准那么中的平安系数1955年, Bishop15首先在边坡稳定性分析中提出了抗剪强度折减的概念， 在 有限元强度折减法中通过将坡体的强度参数：黏聚力 c和内摩擦角©同时除一 个折减系数Ft,得到一组新的c'和&#

7、169;值，作为一个新的强度参数输入进行试算，当计算不收敛时，对应的Ft即为所求的平安系数，此时坡体到达极限状态，发生 剪切破坏。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 积

8、圆屈服准那么中的试验参数如式(10)所示：sin © (9)等面3 =V3(乜？?？?？?？)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= ?+ ?，

9、所以 W =选?+?2 = ?= ?(?)由此可知，n是©的函数，当©取不同值时可以得到不同的 n值如表1所 列：表1不同内摩擦角时的 n值<p/C)】0203()40.ia)L 165L233.301L367<p/ f)5060708090L4281*48()L52Ih 5461,5554. 失稳判据目前两个比拟主流的失稳判据分别是有限元计算中力不平衡和位移的不收敛以及广义塑性应变或者等效塑性应变从坡脚到坡顶贯穿。Griffiths9和郑颖人11121314都使用计算不收敛作为失稳判据。Griffiths9提出，当在用户定义的最大迭代数目下计算仍不收敛时， 那

10、么没有任 何一种应力分布方式可以同时满足 Mohr-Coulomb准那么以及整体稳定，这种情况 可看做边坡失稳判据。边坡失稳与数值计算不收敛同时发生，并伴随着极大的节 点位移，并以1000作为最大的迭代步数。郑颖人14提出，有限元的计算迭代过程就是寻找外力和内力到达平衡状态 的过程，整个迭代过程直到一个适宜的收敛标准得到满足才停止。可见，如果边坡失稳破坏，滑面上将产生没有限制的塑性变形，有限元程序无法从有限元方程 组中找到一个既能满足静力平衡又能满足应力-应变关系和强度准那么的解，此时 不管是从力的收敛标准，还是从位移的收敛标准来判断有限元计算都不收敛。3案例分析例一，不含地基的均质边坡9该边

11、坡如图1所示，有限元程序采用Mohr-Coulomb失效准那么，建立平面应 变条件下八节点四边形单元减缩积分计算模型，其强度参数为=2,边坡坡度为26.57° (2:，坡底水平，其边界条件为坡底约束竖直方向位移与水平方向位 移,左侧约束水平方向位移，其余面为自由面。施加重力荷载后使平安系数从到 逐步变化直至计算不收敛1 仲-一一一+冷H rr图1不含地基的均质边坡每一个平安系数对应的迭代次数如表 2所列，当真正的平安系数接近时需要 更多的迭代次数。表2例一计算结果FOSE rHUK/y h'Iterationso-8o n037921 000 3KIJO1-20042220b

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

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

14、u2/Cu1时，潜在滑面为结 构软弱面。图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 = 02500图5不同软弱层强

15、度时的平安系数例三，不同坡度边坡平安系数计算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.992LO38伦采川蚓按風風

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

17、角75°平均间 距10m。岩体重度为25kN/m3,弹性模量1xi00Pa,泊松比，黏聚力1MPa,内 摩擦角38°两组节理参数相同，重度为17kN/m3,弹性模量1X10Pa,泊松比， 黏聚力，内摩擦角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(11)图11所示为两个直线滑面组成的折线型滑体 ABM

19、CD。岩体重度丫 =20kN/m, 弹性模量E=109Pa。滑块ABCD面积433m2，滑面AB=20m，倾角为15° ,AD=25m ,，; 滑块BCM面积2，滑面，倾角为45°,。CM面上施加有线性变化的面荷载，PM=400kPa, Pc=0。11折线型平面滑动岩质边坡在滑动面AB，BM上布置接触单元，坡体到达极限状态后的破坏滑动如图12所示，并把有限元计算结果，与传统极限平衡方法 Spencer法进行比照，接触单元的相关力学参数以及两种计算结果比照方表5所列图12坡体到达极限状态后的破坏滑动表5例五计算结果有限元强度折减法Spencer 法1.00- l&O

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

21、24kN/m3、 黏聚力400kPa内摩擦角32 °。模型几何尺寸及边界条件如图13所示：p地基图13模型几何尺寸及边界条件通过ABAQUS有限元软件计算，当边坡的平安系数为时，计算不收敛，通 过Slide软件利用瑞典条分法计算得到的边坡平安稳定系数为，瑞典条分法计算 的滑面和有限元计算的塑性区如图14所示，可见两种方法计算出的平安系数和 滑面吻合性较好。瑞典条分法计算的滑动面5参考文献1 Smith, I. M. & Hobbs, R. (1974). Finite element analysis of centrifuged and built-up slopes. G

22、e ?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 plasticity in soil mechanics. Geotechnique 25,671-689:4 Griffiths,

23、 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 59, 237