均质土坡渐进性破坏分析及变分法验证①

刘 芸，李 林

(1.同济大学土木工程学院地下建筑与工程系，上海200092;2.岩土及地下工程教育部重点实验室(同济大学)，上海200092)

0 引 言

随着我国基础建设投资力度的增大，边坡工程的稳定问题也日渐突出.每年我国发生的边坡失稳灾害不计其数，严重威胁人民的生命财产安全［1］.因此，为加强对边坡及其作用形式的认识，正确理解边坡渐进性失稳的机理具有重大的现实指导意义.

就目前而言，分析边坡稳定性的主要方法有两种，分别为传统的基于极限平衡法和有限元数值模拟法［2～4］.但前者的分析方法中，需要做出许多近似假设［5～6］，因此结果中往往得不到滑移体内土体的变形特性以及应力迁移规律，也得不到支护结构对土坡变形及稳定性的影响.近年来，计算机技术不断发展，基于有限元的强度折减系数法也开始受到重视.它可以灵活地结合强度折减法和三维弹塑性有限元分析法，在已知评判指标的前提下，通过不断地调整稳定性折减系数，求出土坡发生破坏前提下的最小稳定安全系数，进而分析判定边坡稳定性.而边坡稳定性分析的变分法以变分原理和极限平衡法为基础研究边坡的稳定性问题［7］.上述两种方法在模拟边坡渐进破坏过程以及确定边坡滑裂面位置方面具有不可比拟的优势.

本文以上海地区某均质土坡工程为背景，构建二维有限元数值模型，模拟土坡的破裂区域的渐进性破坏发展过程，从而根据塑性区范围搜索出潜在危险的破裂滑动面.以变分法为基础，研究平面应变状态下边坡稳定的平衡方程，求解出最危险滑动面的函数表达式，并与数值模拟结果进行了对比验证.

1 强度折减法

1.1 强度折减法基本原理

从Zienkiewicz［8］(1975)提出强度折减概念至今，绝大多数研究者都是沿着他的最初的思路进行探讨:

土的实际强度指标为c 与φ，折减系数为Fr，则强度折减为:

令

则

参数φr与cr即为折减的强度指标，将它们代入有限元计算，若Fr＞1，则有限元计算的位移就比实际大，应力水平也大，破坏单元要多.令Fr从1.0 开始逐步增大，则算得的位移逐步增大，破坏单元逐步增多，乃至最后达到整体失稳.

边坡失稳的判据直接影响强度折减法计算的准确性，现有的失稳判据大致可分为3 种［9～11］:①折减后的强度参数使得计算不能收敛;②坡体内塑性区从坡脚到坡顶贯通;③坡体中的特征点位移或应变发生突变且无限发展.

1.2 边坡渐进性破坏计算思路

根据有限元强度折减法和土坡渐进性破坏规律，动态模拟边坡渐进破坏过程.具体建模流程为:

①确定边坡的地质概况、土体物理力学参数、初始应力场，根据上述信息构建相对应的数值计算模型.

②在折减系数Fr=1 时，进行边坡的弹塑性力学计算，利用屈服接近度指标判断坡体单元是否发生破损，如果没有破损区域的发生，则相应酌情增加折减系数Fr的取值，再重新进行弹塑性力学分析，直至土坡发生局部破损为止停止.

③确定土坡渐进性破坏过程中的塑性破损区，由折减系数Fr折减局部破损区内的强度参数，然后将新的折减系数代入前次的折减有限元计算程序中，重新进行弹塑性力学模拟计算.

④参照步骤①～③，随着折减次数的不断增多以及折减系数Fr的逐渐增大，土坡塑性破损区不断发生扩展，当土坡塑性破损区完全贯通时，土质边坡即达到极限失稳状态，至此土坡的渐进性破坏结束.

2 算例分析

2.1 有限元模型构建

上海地区某土坡高H=10.0m，土体容重为18 kN·m－3，粘聚力c=28kPa，内摩擦角5°.采用四节点平面应变单元，边坡左右两侧设置为水平约束，底面设置为水平和竖向约束，建立二维有限元模型.模型如图1 所示.

图1 计算模型示意图

2.2 结果分析

按照边坡渐进破坏计算步骤，首先取折减系数Fr=1，此时模型所有单元都未出现破损，如图2 所示.

图2 未折减时土坡变形分布

调用第1 次折减计算完成的程序进行第2 次折减计算，此时折减系数Fr=1.02(图3(a))，针对根据该折减模型，再次进行弹塑性力学计算.

按照上节所述的计算流程，此时可以对折减系数Fr进行反复自动搜索，折减计算过程中土坡破损区演化如图3 所示.在第7 次折减计算时(此时折减系数Fr=1.15)，有限元模拟计算结果正好显示出土坡完全出现贯通塑性破损区，此时土坡处于极限失稳状态.由图可见潜在破裂滑动面在起始时期变化发展比较缓慢，而在后期阶段滑动面破裂贯通区快速发展，这一演化过程也符合边坡滑坡的发展特点.

图3 折减计算过程中破损区演化

3 变分法验证

边坡稳定性分析的变分法计算，其关键在于利用最危险滑动面函数和最危险滑动面上的正应力分布函数建立土条的静力平衡方程，通过构造相应泛函并求解其极值，得到最危险滑动面函数和最危险滑动面上的正应力分布函数的表达式，进而研究边坡的稳定性问题.

3.1 土坡稳定性分析模型

为验证基于强度折减法的均质土坡最危险滑裂面位置，以简单边坡坡角A 为原点建立直角坐标系，如图4 所示.

图4 直角坐标系边坡计算模型

边坡坡比1:n，高度H，均质边坡土体黏聚力c，内摩擦角φ，容重γ，坡面函数ym(x)可写为两分段函数，假设滑动面函数为y(x).取微元体土条进行受力分析，有

则:

土体重力:

式中，σ，τ 分别为沿滑动面函数y(x)分布的正应力和剪应力，α 为τ 与水平面的夹角，l 为滑动面长度.

本文考虑均质边坡在土体自重、滑动面上的正应力和剪应力作用下达到平衡，则平衡方程为:

滑动面贯穿边坡，端点A(x1，y1)，C(x2，y2)，土体屈服准则选用Mohr－Coulomb 准则，即

代入式(7)，有

参考变分原理［12］和文献［13］，对长度和应力无量纲化，记

则有:

边坡分析模型转化为图5 所示.

图5 边坡计算模型

3.2 最危险滑动面变分法求解

由平衡方程知，最危险滑动面函数Y(X)及沿该滑动面分布的正应力函数σ(X)应满足式(10).记泛函

式中λ1，λ2为待定常数.边坡最危险滑动面问题即转化为泛函的极值问题:寻找Y(X)及σ(X)使G 取极值.由变分原理，上述问题的解Y(X)、σ(X)应满足欧拉方程.

由式(14)可以确定滑动面及其正应力分布的控制微分方程，选用Mohr－Coulomb 屈服准则时，边坡最危险滑动面与最危险滑动面上的正应力分布无关.求解式(14)得:

由式(15)即可确定滑动面函数Y(X)和沿该滑动面分布的正应力函数σ(X)，分2 种情况讨论极值函数.

(1)当λ2=0 时，由式(15)得

当μ=tanφ 为常数，即边坡为均质情况时，式(16)确定的最危险滑动面函数Y(X)为直线，表达式为

为便于计算，引入关于未知极点O(X0，Y0)的极坐标系，如图5 所示.令

式中，

代之入式(15)中

即

且由式(18)和式(19)得:

式中C1，C2为积分常数.

式(23)和式(24)即为极坐标系下最危险滑动面及与之相对应的正应力表达式，为方便编程计算，将式(10)按照式(18)、(19)进行坐标变换，得

且由边界条件，得

边坡坡面函数表达式为:

边坡B 点坐标为(nH，H)，(n，1)，(Rm(θ)，arctan((λ1－nλ2)/(1+λ2))).假设边坡最危险滑动面一端点坐标A 为(0，0)，即最危险滑动面经过坡角［14］，则

得

联立式(25)和式(26)，得

式中:

式(29)待求量共有5 个:待定常数λ1，λ2，积分常数C1，C2及未知坐标θ2，给定均质边坡坡比1:n、坡高H、土体黏聚力c、内摩擦角φ 和容重γ 参数后，按照式(29)5 个方程即可求解.

3.3 结果分析

输入均质土坡已知参数:坡比1:n，坡高H，土体黏 聚 力 c，内 摩 擦 角 φ 和 容 重 γ，利 用Mathematicas 软件编程求解式(29)，结果见表1.

表1 土坡参数及计算结果

利用坐标变换式(18)和(19)将边坡最危险滑动面转换到直角坐标系中.由此，即利用变分法得到了均质边坡的最危险滑动面，极坐标系下边坡最危险滑动面的表达式(23)说明其为对数螺旋面.并将理论结果与数值模拟结果进行对比，如图6 所示.

图6 最危险滑裂面对比图

从图6 中可以看出由强度折减法确定的土坡最危险滑裂面与采用变分法计算的解析结果趋势大致相同，均为圆弧形滑裂面.但需要指出的是，运用变分法计算的最危险滑裂面位置要更靠近坡面，并且圆弧形滑裂面的曲率较强度折减法计算所得曲率要稍大.对比分析结果还表明，采用强度折减法模拟均质土坡渐进破坏过程及确定其最危险破裂面是准确、可行的.

4 结 论

本文以上海地区某均质土坡工程为背景，构建了二维有限元数值模型，动态模拟了土坡的塑性区发展和渐进破坏过程.以变分法为基础，研究平面应变状态下边坡稳定的平衡方程，求解出最危险滑动面的函数表达式，并与数值模拟结果进行了对比验证.得出以下主要结论:

(1)随着有限元模型折减次数的不断增加，稳定性折减系数Fr 也在逐渐增大，土质边坡塑性破损区不断动态发展破坏，边坡的破裂面不断发展形成，土体的稳定性逐渐下降，当土坡滑裂面贯通时，土坡达到极限失稳状态，渐进破坏完成.

(2)边坡稳定性分析的变分法计算，其关键在于利用最危险滑动面函数和最危险滑动面上的正应力分布函数建立土条的静力平衡方程，通过构造相应泛函并求解其极值，得到最危险滑动面函数和最危险滑动面上的正应力分布函数的表达式.

(3)由强度折减法确定的土坡最危险滑裂面与采用变分法计算的解析结果趋势大致相同，均为圆弧形滑裂面.运用变分法计算的最危险滑裂面位置要更靠近坡面，且圆弧形滑裂面的曲率较强度折减法计算所得曲率要稍大.

［1］ 郑颖人，赵尚毅，时为民.边坡稳定性分析的一些新进展［J］.地下空间，2001，21(4):263－271.

［2］ 丁桦，张均锋，郑哲敏.关于边坡稳定分析的通用条分法的探讨［J］.岩石力学与工程学报，2004，23(21):3684－3688.

［3］ 张士兵，王练柱，张建.边坡稳定性弹塑性大变形有限元强度折减分析［J］.岩石力学与工程学报，2004，23(1):4463－4467.

［4］ 于玉贞，林鸿州，李荣建，等.非稳定渗流条件下非饱和土边坡稳定分析［J］.岩土力学，2008，29(11):2892－2898.

［5］ 许建聪，尚岳全，陈侃福，等.顺层滑弹塑性接触有限元稳定性分析［J］.岩石力学与工程学报，2005，24(13):2231－2236.

［6］ DUNCAN J M.Fellow State of Art:Limit Equilibrium and Finite Element Analysis of Slopes［J］.Journal of Geotechnical Engineering，ASCE，1996，122(7):577－596.

［7］ YAMAGAMI T，UETA Y.Search for Noncircular Slip Surfaces by the Morgenstern Preice Method［J］.Proc.of 6th Int.Conf.on Numerical Methods in Geomechanics，1988，17(1):1335－1340.

［8］ ZIENKIEWICZ O，HUMPHESON C，LEWIS R.Associated and Non－associated Visco－plasticity in Soil Mechanice［J］.Geotechnique，1975，5(4):671 ～689.

［9］ GRIFFITHS D V，LANE P A.Slope Stability Analysis by Finite Elements［J］.Geotechnique，1999，49(3):387－403.

［10］ 陈力华，靳晓光.有限元强度折减法中边坡三种失效判据的适用性研究［J］.土木工程学报，2012，45(9):136－146.

［11］ 梁庆国，李德武.有限元强度折减法中边坡三种失效判据的适用性研究［J］.岩土力学，2008，29(11):3053－3 058.

［12］ Dov Leshchinsky，San K C.Pseudostatic Seismic Stability of Slopes:Design Charts［J］.Journal of Geotechnical Engineering，1994，120(9):1514－1532.

［13］ Shen W Y，Teh C I.Practical Solution for Group Stiffness Analysis of Piles［J］.Journal of Geotechnical and Geoenvironmental Engineering，2002，128(8):692－698.

［14］ 柳厚祥，廖雪，李宁.公路边坡稳定性分析的二维变分方法［J］.中国公路学报，2007，20(4):7－11.