基于强度折减的边坡稳定安全系数有限元迭代解法

周桂云

(1.金陵科技学院建筑工程学院,江苏南京 211169;2.中国人民解放军理工大学工程兵工程学院,江苏南京 210007)

基于强度折减的边坡稳定安全系数有限元迭代解法

周桂云1,2

(1.金陵科技学院建筑工程学院,江苏南京 211169;2.中国人民解放军理工大学工程兵工程学院,江苏南京 210007)

讨论有限元强度折减法的失稳判据,结果表明:失稳判据对安全系数影响较大,对临界滑动面形状和位置的确定影响不大。在此基础上,提出基于强度折减的边坡稳定安全系数有限元迭代解法,先用强度折减法搜索滑动面,再按滑动面已知边坡用迭代解法求解稳定安全系数。该方法充分结合了强度折减法和迭代解法的优点,求解边坡稳定安全系数时效率快、精度高,可为实际工程参考应用。

强度折减法;迭代解法;边坡稳定安全系数;有限元法

20世纪60年代以来,有限元法被广泛应用于边坡稳定分析,已成为近年来边坡稳定分析的新趋势。有限元法克服了极限平衡法的不足,并有其独特的优点:①能够充分考虑土体的非线性本构关系、土体的非均质性和各向异性等复杂特性;②能模拟施工过程,可考虑变形体中泥化夹层的渗流效应、孔隙水压力等;③能给出边坡内的应力场、位移场,可以了解边坡的逐步破坏机理和塑性区的开展情况;④能适应复杂的边界条件,比如复杂的几何形状、结构物的作用、有支护的边坡稳定;⑤可以真实考虑水库水位骤降、地震荷载等对边坡稳定性的影响。Duncan[1]对有限元法在边坡和堤坝稳定分析中的应用研究进行了详细的归纳总结,按照所采用的本构模型将它们分为4类,即:①线弹性模型;②折线弹性模型;③双曲线弹性模型;④弹塑性或弹黏塑性模型。从这些总结可以看出,大部分研究工作是对边坡进行有限元分析,对计算范围内各单元上积分点的应力进行强度判别,凡应力状态达到剪切破坏或拉破坏判别标准的部位称之为破坏区,根据破坏区的分布位置和范围的大小可以对边坡的稳定性进行评价,为边坡的治理、施工方法提供依据。随着计算机的发展,有限元法作为一种强有力的数值计算方法,在边坡稳定分析中已经取得了大量的成果。目前主要的方法有2类:①有限元迭代解法;②有限元强度折减法。迭代解法和强度折减法在边坡稳定安全系数的定义上是相同的,不同的是迭代解法只适用于滑动面预先知道的情况下求解边坡的稳定安全系数,而抗剪强度折减法主要用于解决滑动面未知的情况(当然也可以用于滑动面已知的情况)。对滑动面未知情况可以利用破坏时刻的塑性区分布图来确定滑动面的位置。2种有限元法各有优缺点,本文在对有限元强度折减法的失稳判据进行讨论的基础上,结合有限元强度折减法和迭代解法,提出基于强度折减的边坡稳定安全系数有限元迭代解法,使得求解边坡稳定安全系数时速度快且精度高。

1 有限元强度折减法的失稳判据讨论

有限元强度折减法将强度折减技术与弹塑性有限元方法相结合,首先进行强度折减系数下边坡的弹塑性有限元分析,得到边坡内的应力场、应变场、位移场,然后以位移、广义剪应变等描述变形程度的某种物理量作为评判指标,定量地描绘边坡的潜在塑性破坏区域及其破坏程度与发展趋势,据此基于一定的经验评判准则确定边坡的极限平衡状态,并将由此所确定的相应强度折减系数作为边坡的稳定安全系数[2-4]。

式中:K为安全系数;c和φ为土体的实际抗剪强度参数;ce和φe为土体达到极限状态时的抗剪强度参数。

有限元强度折减法不需要事先确定滑动面的形状和位置,并可以按一定的规则搜索出最危险临界滑动面,因而能解决复杂的工程问题。然而,该方法的关键问题是如何根据有限元计算结果来判断边坡是否处于破坏状态,即如何定义边坡的失稳判据,至今工程界还没有统一的标准。目前,主要有3种失稳判据:

a.以有限元迭代求解过程的不收敛作为边坡失稳的判据[5-9]。在非线性有限元迭代求解时,给出迭代收敛标准,若在给定的迭代次数和收敛准则内有限元计算仍未收敛,则认为发生破坏。迭代收敛准则有位移准则、失衡力准则和能量准则3种。

b.以塑性区(或等效塑性应变)从坡脚到坡顶贯通作为边坡失稳的判据[10-12]。根据强度折减弹塑性有限元分析结果,绘制坡体内的塑性应变(广义剪应变)分布图,当某一幅值的塑性应变(广义剪应变)开展区贯通时则认为边坡失稳破坏。

c.以坡体内特征部位位移发生突变作为边坡失稳的判据[13-14]。边坡的变形破坏总会具有一定的位移特性,位移结果是边坡失稳最直观的表达形式。根据有限元计算结果建立特征点的位移和折减系数关系曲线,若位移发生突变且无限发展则认为边坡失稳破坏。

通过研究分析,上述3种失稳判据各有优缺点:由于岩土体的复杂性,影响迭代计算不收敛的因素很多,如预先定义的收敛标准、本构方程形式、荷载步长、初始地应力、土体材料特性等,且收敛准则的定义具有一定的人为任意性,所以以迭代过程的不收敛作为失稳判据是不合理的,适用性差;以塑性区贯通作为失稳判据具有较明确的物理意义,然而塑性应变的大小受到材料类型的影响,而且以某一定量值作为标准缺乏广泛的适用性,边坡是在塑性应变贯通前失稳还是贯通后失稳,缺乏依据;以坡体内特征部位位移发生突变作为失稳判据的物理意义明确且直观,然而,其缺陷是究竟用哪个点以及哪个方向的位移没有统一的认识,如何从曲线上给出安全系数也没有明确的方法。考虑到实用性与简便性,文献[15][16]建议联合采用特征部位位移突变和塑性区贯通作为失稳判据。

为了分析失稳判据对稳定安全系数和滑动面的影响,假定一均质土坡如图1所示,土的计算参数为:弹性模量E=2.0×105Pa,泊松比 ν=0.25,密度ρ=2000kg/m3,黏聚力c=1.0×104Pa,摩擦角 φ=20°。

图1 均质土坡计算简图

有限元网格剖分如图2所示,由有限元计算可知,坡体是由下往上逐渐贯通破坏的,所以选择处于塑性区范围内的A点为特征点。以有限元迭代求解过程的不收敛为失稳判据时求得稳定安全系数为1.1236,以塑性区贯通作为失稳判据时求得稳定安全系数为1.1494,以特征部位位移发生突变为失稳判据求得稳定安全系数为1.1628。图3为3种失稳判据下边坡失稳时的塑性区分布。

图2 均质土坡有限元模型

经对比分析可知,有限元计算迭代是否收敛与边坡失稳没有必然的关系,迭代不收敛时边坡不一定发生失稳;塑性区贯通是破坏的必要条件,但不是充分条件,塑性区贯通并不一定意味着破坏;位移突变是破坏的充分条件,边坡破坏总有一定的位移特征。所以建议结合塑性区贯通和位移突变来判别边坡是否破坏。计算结果表明,不同失稳判据求得的稳定安全系数不同,且失稳判据对稳定安全系数的影响较大。但是塑性区分布相差不大,若以塑性区分布图上塑性应变较大点的连线作为临界滑动面,则3种失稳判据下的临界滑动面很接近,说明失稳判据对临界滑动面形状和位置的确定影响不大。

图3 3种失稳判据下的塑性区分布

2 基于强度折减的边坡稳定安全系数有限元迭代解法

针对有限元强度折减法的上述问题,为了避免由于失稳判据对安全系数的影响,笔者提出改进方法,即基于强度折减的边坡稳定安全系数有限元迭代解法。强度折减法之所以方便实用是因为此法不需要事先假定滑动面的形式和位置,然而该方法的关键是失稳判据的定义,且失稳判据对稳定安全系数影响较大,对临界滑动面位置影响不大,而迭代法适用于已知滑动面的情况求解边坡的稳定安全系数。因此,基于强度折减的边坡稳定安全系数有限元迭代解法的思路为:对未知滑动面的边坡先用有限元强度折减法分析,得到临界滑动面,再用有限元迭代解法按已知滑动面情况求解稳定安全系数。在已知滑动面的有限元迭代求解时,关于滑动面的单元形式在文献[17]中有详细的讨论,提出了常规的矩形(平面)或立方体(空间)单元的统一模式。对厚度为零的滑动面采用较小厚度取代零厚度使得问题有解。计算分析结果表明,一般可近似地取结构面的厚度t=1/100l,其中l为与滑动面相连的实体单元的最小边长。滑动面的材料参数取与滑动面相连的2种材料中较小尺寸相对应的材料参数。滑动面上的应力强度准则一般采用摩尔-库仑准则,同时满足相应的抗拉条件。

3 算例分析

3.1 算例1

图4为一垂直坡覆盖了一宽度为B的刚性条形基础,基础中部施加垂直向下的集中荷载F,不考虑自重,右侧和底部施加约束,其理论极限荷载值为F/(Bcu)=2.0,其中cu为土体的黏聚力。由强度折减法的抗滑稳定安全系数定义有F/[B(cu/K*)]=2.0,则该坡体理论安全系数为K*=2.0Bcu/F。本次计算中取F=1500N,B=5m,cu=200Pa,摩擦角 φ=0°,弹性模量E=1.0×105Pa,泊松比ν=0.495。理论破坏面与水平面成45°角方向,理论安全系数K*=1.3333。

图4 垂直剪问题模型尺寸

先用有限元强度折减法搜索滑动面,未知滑动面网格剖分如图5所示,计算所得塑性区分布如图6所示。根据折减法计算所得的塑性区分布确定滑动面,作为已知滑动面问题用有限元迭代解法求解稳定安全系数,已知滑动面网格剖分如图7所示,取迭代初始稳定安全系数K1=1.0,经过4次迭代计算,结果为K3=K4=1.3309,与理论解仅相差1.8%,破坏时的位移矢量如图8所示,破坏滑移变形如图9所示。

图5 未知滑动面网格剖分

图6 塑性区分布

图7 已知滑动面网格剖分

图8 位移矢量

图9 破坏变形

3.2 算例2

采用图1所示的均质土坡,土的计算参数如下:弹性模量E=2.0×105Pa,泊松比 ν=0.25,密度 ρ=2000kg/m3,黏聚力c=1.0×104Pa,摩擦角 φ=20°。

未知滑动面计算时网格剖分如图10所示,这里网格剖分较密,可以比较准确地确定滑动面,塑性区分布如图11所示,根据塑性区分布图确定滑动面后剖分网格如图12所示。经过3次迭代满足误差要求,稳定安全系数K=1.1732,与文献[18]中不同网格计算所得稳定安全系数基本一致。破坏时的位移矢量如图13所示,破坏变形如图14所示。

图10 未知滑动面网格剖分

图11 塑性区分布

图12 已知滑动面网格剖分

图13 位移矢量

图14 破坏变形

4 结 语

对有限元强度折减法的失稳判据进行讨论,建议结合使用塑性区贯通和位移突变作为边坡失稳判据,并认为失稳判据对稳定安全系数影响较大,但对搜索临界滑动面的形状和位置影响不大。在此基础上,提出基于强度折减的边坡稳定安全系数有限元迭代解法。算例结果表明,该方法方便有效,充分利用了强度折减法能够方便准确地搜索滑动面的优点和迭代法计算稳定安全系数效率高、精度高的优点,是一种求解边坡稳定安全系数的快速、简便方法,可为实际工程应用提供参考。

[1]DUNCAN J M.State of the art:limit equilibrium and finite element analysis of slopes[J].Journal of Geotechnical Engineering,ASCE,1996,122(7):477-596.

[2]ZIENKIEWICZ O C,HUMPHESON C,LEWIS R W.Associated and non-associated visco-plasticity and plasticity in soil mechanics[J].Geotechnique,1975,25(4):671-689.

[3]宋二祥.土工结构安全系数的有限元计算[J].岩土工程学报,1997,19(2):1-7.

[4]张鲁渝,郑颖人,赵尚毅,等.有限元强度折减系数法计算土坡稳定安全系数的精度研究[J].水利学报,2003,34(1):21-27.

[5]GRIFFITH D V,LANE P A.Slope stability analysis by finite element[J].Geotechnique,1999,49(3):387-403.

[6]DAWSON E M,ROTH W H,DRESCHER A.Slope stability analysis by strength reduction[J].Geotechnique,1999,49(6):835-840.

[7]赵尚毅,郑颖人,时卫明,等.用有限元强度折减法求边坡稳定安全系数[J].岩土工程学报,2002,24(3):343-346.

[8]赵尚毅.用有限元强度折减法及其在土坡与岩坡中的应用[D].重庆:后勤工程学院,2004.

[9]郑颖人,赵尚毅.用有限元强度折减法求滑(边)坡支挡结构的内力[J].岩石力学与工程学报,2004,23(20):3552-3558.

[10]郑宏,李春光,李焯芬,等.求解安全系数的有限元法[J].岩土工程学报,2002,24(5):626-628.

[11]连镇营,韩国城,孔宪京.强度折减有限元法研究开挖边坡的稳定性[J].岩土工程学报,2001,23(4):406-411.

[12]栾茂田,武亚军,年延凯.强度折减有限元法中边坡失稳的塑性区判据及其应用[J].防灾减灾工程学报,2003,23(3):1-8.

[13]宋二祥.土工结构安全系数的有限元计算[J].岩土工程学报,1997,19(2):1-7.

[14]周桂云,李同春.饱和-非饱和非稳定渗流作用下岩质边坡稳定性分析[J].水电能源科学,2006,24(5):79-82.

[15]姚纬明,李同春,任旭华,等.带软弱结构面的弹塑性有限元分析[J].河海大学学报:自然科学版,1999,27(3):34-38.

[16]刘金龙,栾茂田,赵少飞,等.关于强度折减有限元方法中边坡失稳判据的讨论[J].岩土力学,2005,26(8):1345-1348.

[17]李同春,卢智灵,姚纬明,等.边坡抗滑稳定安全系数的有限元迭代解法[J].岩石力学与工程学报,2003,22(3):446-450.

[18]张丹.边坡稳定及抗滑桩内力分析方法研究[D].南京:河海大学,2004.

Finite element iteration method for safety factor of slope stability based on strength reduction method

ZHOU Gui-yun1,2(1.Institute of Civil Engineering and Architecture of Jinling Institute of Technology,Nanjing211169,China;2.Engineering Institute of PLAUniversity of Science and Technology,Nanjing210007,China)

The instability criterion of strength reductionmethodwas discussed.The instability criterionhad a great influence onthe safety factor,but has a little influence on the shape and position of the critical sliding surface.On this basis,the finite element iteration method for the safety factor of slopesbased on the strength reductionmethod was proposed.The critical sliding surface was searched by use of the strength reductionmethod,and then the safety factor of given slopes of the sliding surface was solved by the finite element iteration method.A conclusion is drawn that the proposed method fully combines the advantages of the strength reduction method and the iteration method.It is fast,accurate and practical in solving the safety factor for the slope stability.

strength reduction method;iteration method;safety factor of slopes;finite element method

TU311.3

A

1006-7647(2010)03-0058-04

10.3880/j.issn.1006-7647.2010.03.015

中国博士后科学基金(20080431343);金陵科技学院科研基金(jit-n-2009-014)

周桂云(1980—),女,江苏泰州人,讲师,博士,从事土木工程教学及科研工作。E-mail:guiyunzhouhhu@yahoo.com.cn

book=94,ebook=275

2009-07-17 编辑:高建群)