梁 艳 李同春
(河海大学 水利水电学院,南京 210098)
有限元强度折减法是Zienkiewicz于1975年提出的一种边坡稳定分析方法,与传统极限平衡法相比而言,有限元强度折减法通过对某一强度折减系数下的边坡进行弹塑性有限元分析,得到其应力场、应变场及位移场;然后,根据一定的失稳判据对边坡的稳定情况进行判断.不需要事先确定滑裂面的位置和形状,最危险临界滑裂面也可以按一定的规则搜索出,因而能解决复杂的工程问题.随着计算机技术的迅速发展以及各种大型成熟软件的推出,有限元强度折减法逐渐成为边坡稳定分析研究的热点,得到学术界及工程界的广泛关注.现行的边坡失稳判据主要有以下几种:①以有限元迭代求解过程的不收敛作为边坡失稳的判据[1];②以塑性区(或等效塑性应变)从坡脚到坡顶贯通作为边坡失稳的判据[2];③以坡体内特征部位位移发生突变作为边坡失稳的判据[3-7].
失稳判据的选择是有限元强度折减法的关键问题,很多学者对此进行研究,但到目前为止没有形成统一的认识.本文利用有限元强度折减法对典型边坡算例进行分析,研究了非线性迭代方法对计算结果的影响,并对位移曲线拐点状态及塑性区判据进行了探讨分析.
强度折减法中边坡稳定的基本原理基于强度储备概念,其安全系数定义为:使边坡刚好达到临界破坏状态时,对岩土体的抗剪强度进行折减的程度,即定义安全系数为岩土体的实际抗剪强度与临界破坏时折减后剪切强度的比值.强度折减法的关键是利用公式(1)和(2)来调整岩土体的强度指标c和φ,然后对边坡稳定性进行数值分析,不断的增加折减倍数,反复计算,直至其达到临界破坏,此时得到的折减倍数即为安全系数Fs,其中折减倍数的倒数即为折减系数K.
式中,c为土的粘聚力,也称内聚力;φ为土的内摩擦角,即抗剪强度线的倾角;cF为折减后的粘结力;φF为折减后的摩擦角;Ftrial为折减倍数.
以经典边坡为例[1,8],算例边坡几何尺寸及有限元模型,如图1所示.
图1 边坡算例几何尺寸及有限元网格模型
其具体材料参数见表1.表1中γ表示土体重度,E为土体弹性模量,υ表示土体泊松比.本算例属于平面应变问题,边坡土体采用理想弹塑性 Mohr-Coulomb屈服准则及非关联流动法则.边界条件设左右两边为水平约束,底边为固定约束,坡面为自由边界.采用四节点网格进行分析,在坡面坡脚处网格划分较密,有限元模型节点数4 961,单元数4 800.
表1 算例土体力学参数表
对于该算例,根据不同方法给出的安全系数如表2所示.
表2 不同分析方法得到的安全系数值
对此经典边坡算例进行有限元强度折减法计算分析,并对不同的判据进行讨论.
有限元计算的收敛性与非线性方程的解法及收敛容差紧密相关.非线性问题的求解是通过对荷载步内不平衡量进行线性迭代,最终获得其收敛解的一个过程.在实际的有限元分析中,主要有两种求解方法[10]:常刚度迭代法与变刚度法.将最大迭代步数设为1 000步,一次性施加重力荷载,采用两种求解方法对经典边坡进行计算,当收敛容差为0.01时,变刚度法在安全系数等于1.08时不收敛,常刚度法安全系数等于1.36时不收敛.以安全系数为1.08为例,常刚度与变刚度迭代法收敛曲线走势如图2所示.
图2 K=1.08时常刚度与变刚度迭代收敛曲线走势
由图可得当安全系数等于1.08时,常刚度迭代法收敛量级随迭代次数增加而减小直至达到收敛标准;而变刚度法收敛量级随迭代次数增加先减小后又增大,逐渐远离收敛标准,已经不收敛,但此时不收敛并不代表完全屈服,对应塑性区分布如图3所示.
图3 折减系数为1.08时的塑性区分布图
在此基础上,通过对其他10组不同参数边坡进行计算,对其收敛走势进行研究得:当土体大部分屈服时,变刚度法基本不收敛,无法计算出土体破坏时的安全系数.采用常刚度法,在边坡接近破坏时,迭代次数常常很大,收敛较慢,但基本能保证收敛.
当采用常刚度法进行计算时,不同收敛容差下得到的安全系数变化较大,见表3.
表3 常刚度迭代法在不同收敛容差下安全系数比较
从上面的结果可以看出以有限元计算不收敛作为判据时,计算的结果与所用的迭代计算方法和收敛容差有很大关系.安全系数的确定具有较大的人为任意性.除此之外有限元计算收敛性受多重因素综合影响,如有限元计算模型,计算单元类型,地应力影响,计算边界等,目前也缺乏行之有效的方法来消除这些影响.在这种情况下,有限元计算收敛性和边坡极限平衡状态之间似乎没有直接的一一对应关系,计算不收敛并不意味着边坡已达到极限平衡状态,而计算收敛也不能表明边坡是安全的,以此为判据不够合理.
将边坡整体作为一个结构来讨论,则其位移曲线的拐点应该有两种情况:第1种是结构从弹性状态变为非线性状态;第2种是结构从非线性状态变为失稳状态.边坡的变形破坏总具有一定的位移特性,因此有限元计算的位移结果是边坡失稳最直观的表达.目前以位移作为失稳判据的方法是建立每次有限元计算的某个部位的位移或者最大位移与折减系数的关系曲线,以曲线上的拐点作为边坡处于临界破坏状态的临界点.以边坡顶部临空面端点A及坡脚B为特征点,研究A点竖直位移及B点水平位移随安全系数变化规律.若按变刚度法,在边坡临近破坏时,计算不能保证收敛,且边坡接近破坏时刚度趋向于零,无法求出位移.常刚度法计算一方面能保证收敛,另一方面能算出位移曲线的拐点,因此采用常刚度法计算,位移-安全系数曲线如图4所示,可知在1.113和1.176两处均有拐点.
图4 A点竖向位移与B点水平位移随安全系数变化曲线
第1种情况不一定存在,与边坡初始状态相关,如图5所示.
图5 不同材料边坡位移随安全系数变化曲线
当材料强度高时,边坡起始状态处于弹性状态,两个明显的拐点,当材料强度低时,边坡起始状态出于非线性状态,只有一个明显的拐点.但第1个拐点与边坡稳定定义不符;第2个拐点出现后用常规的基于小变形假定的有限元解法无法求解,即在第2个拐点之前计算得到的数据是准确的,而拐点之后边坡已经滑动,此时得到的位移是强行按照小变形假定来计算,并不是其真实值,但仍能体现边坡位移的发展趋势,拐点的物理意义明确.因此判断此边坡的安全系数为1.176.
由于土体是弹塑性的,当应力达到一定程度时,土体便会发生塑性破坏,岩土体的塑性破坏与塑性区出现、扩展及其分布紧密相关.以塑性区贯通为失稳判据,得出安全系数为1.150,如图6~7所示.边坡破坏时,其塑性应变必然是贯通的,但是在用有限元法计算时,力与位移之间的关系并不是精确推导出来的,而是利用每一单元中近似的位移函数得到节点位移,然后计算高斯点应变与应力,输出时将积分点结果线性外推至单元节点上,因此由高斯点外推得到节点塑性应变值,在滑动面未知情况下塑性区分布在一条带上,当塑性应变区贯通时,塑性破坏点并不一定贯通,还可以继续发展,因此以塑性区贯通作为失稳判据难以准确把握临界点.
基于以上分析,塑性区贯通判据在实际应用中可作为位移拐点判据(尤其当位移-折减系数曲线只有一个拐点时)的补充判据,结合位移场及塑性区变化进行对比分析,在位移曲线出现拐点时,对应塑性区全部贯通,此时对应的折减系数即为安全系数.
鉴于目前有限元强度折减法的失稳判据不统一,进行了讨论,得出以下结论:
1)采用有限元计算不收敛作为判据时,计算的结果与所采用的迭代计算方法和收敛容差有很大关系.安全系数的确定具有较大的人为任意性,以此为判据不够合理.
2)完全接近失稳时变刚度法很难收敛或不收敛,此时不收敛不代表完全屈服,常刚度法基本上能收敛,但收敛速度很慢或很难达到完全收敛,但采用此法能出现位移拐点.
3)以塑性区贯通作为失稳判据难以准确把握临界点,在实际应用中可作为位移拐点判据的补充判据.
4)以位移曲线拐点作为判据时,曲线可能会出现两个拐点,拐点的情况与边坡初始状态相关.建议以位移出现拐点与塑性区全部贯通相结合作为边坡失稳判据.
[1] 赵尚毅,郑颖人,时卫民,等.用有限元强度折减法求边坡稳定安全系数[J].岩土工程学报,2002(3):343-346.
[2] 栾茂田,武亚军,年廷凯.强度折减有限元法中边坡失稳的塑性区判据及其应用[J].防灾减灾工程学报,2003(3):1-8.
[3] 周桂云,李同春.饱和-非饱和非稳定渗流作用下岩质边坡稳定性分析[J].水电能源科学,2006(5):79-82+101-102.
[4] 曹泽伟.重力坝深层抗滑稳定分析[D].南京:河海大学,2011.
[5] Chen Lihong,Shu Yu,Zhang Hongtao.Discussion of Criteria of Shear Strength Reduction FEM[A]//1st International Conference on Civil Engineering,Architecture and Building Materials[C].Haikou,China:Trans Tech Publications,2011.
[6] 段庆伟,陈祖煜,王玉杰,等.重力坝抗滑稳定的强度折减法探讨及应用[J].岩石力学与工程学报,2007(S2):4510-4517.
[7] 宋二祥.土工结构安全系数的有限元计算[J].岩土工程学报,1997,19(2):1-7.
[8] 周翠英,刘祚秋,董立国,等.边坡变形破坏过程的大变形有限元分析[J].岩土力学,2003(4):644-647,652.
[9] 顾芳芳.边坡稳定分析的若干问题研究[D].南京:河海大学,2010.
[10]Ian M Smith,D V Griffiths.Programming the Finite Element Method,Third Edition[M].US:John Wiley &Sons,Inc,2004.