基于重度增加法在土坡稳定性计算上的改进★

王少华 许小明 蔡 瑾

(湖北三峡职业技术学院，湖北宜昌 443000)

边坡稳定性分析理论经过多年的研究发展已经逐步成熟，极限平衡法主要表现在其力学模型简单，无法反应模型内部应力应变情况以及破坏机理，且在安全系数计算方面操作繁琐;有限元法虽然克服了极限平衡法以上各方面的不足，也能计算各种复杂边界的边坡，但不能和安全系数建立直接联系而给出一个量化的评价标准，在一定程度上限制了其应用范围。本文主要以重力加载法确定的潜在软弱面为指定滑弧，计算土坡极限平衡状态下的安全系数。

1 相关软件介绍

ANSYS程序是一个功能强大、灵活的设计分析优化分析软件包，可以浮动运行于从PC，NT工作站直至巨型机的各类计算机及操作系统中，与众多先进CAD软件共享数据接口，数据文件在其所有的产品系列和工作平台上均兼容。

Slide是一款评价岩质或土质边坡安全系数或者失效概率的二维极限平衡程序，滑面兼具圆弧和非圆弧形式，可指定已知滑面或者驱动程序使之自动搜索滑面。

本文在研究过程中均采用ANSYS，Slide进行分析。

2 流动法则及相应的本构模型选取

本文假定被研究模型为理想弹塑性材料，其遵循关联流动法则，更符合岩土体的实际情况。

3 研究步骤

大变形有限元重度增加法与极限平衡综合的基本思路分为三步［1-3］:第一步，通过大变形有限元法重度增加法确定边坡的最危险滑动面及坡体内应力应变情况;第二步，按极限平衡条分法原理将滑坡体分条，分析条块的受力情况，并推导安全系数公式以及编制计算程序;第三步，根据最危险滑动面位置和形状拟合成滑动面方程，将滑动面方程以及岩土体材料输入计算程序，精确求解安全系数。研究方法用结构示意图可以简单的表示，如图1所示。

4 有限元法潜在滑动面的确定［4，5］

边坡滑动面以及潜在滑动面［6］作为边坡稳定分析的一个重要环节，如从传统极限平衡法的角度来看，确定边坡临界滑动面时，搜索范围的指定具有任意性，并且是基于一定假设的基础之上，不一定能够搜索到真正的危险滑移面。

根据本文思想，如通过数值模拟来确定滑动面位置所在，则失稳准则的选择在整个分析过程中显得尤为重要。判定边坡失稳的理论方法目前概括起来只有三种:1)有限元数值迭代计算不收敛;2)特征部位特征点位移突变;3)等效塑性应变从坡脚至坡顶的贯通。由于塑性区贯通是破坏的必要非充分条件，那么，危险滑移面必定发生在塑性区贯通面上，本文尝试通过有限元重力加载法找到该塑性区贯通面。

图1 研究方法结构示意图

具体问题具体分析，以下用模型实例来探讨在重力增加的过程中，塑性区的扩展(渐变)过程。算例基本资料:淹锅沙坝某均质土坡，边坡模型以及物理力学参数如表1所示。

表1 物理力学参数表

通过大变形有限元重度增加法在有限元软件ANSYS中模拟边坡破坏过程:随着重力加速度的增大，坡内塑性区不断重分布并不断增大，最大塑性应变逐渐形成塑性应变等值脊线密集区，并不断向上发展直至贯通，形成明显的滑动带，如图2中MX标志处。在重度增加到1．82g时，塑性区已基本贯通，但未形成滑移面，在重度增加至1．84g时，塑性区域内形成了贯通的滑移面，但不能判断是否已经发生了大变形，而在重度增加到1．86g时，滑动面上部开始闭合，说明大变形已经发生，故在1．82～1．84之间的任意一个数值都可能作为超载储备安全系数，如果人为观察确定塑性区刚刚贯通时的数值这个操作是难以实现的，所以单以塑性区贯通为标志来确定安全系数值是不合理的，从而也间接证明了重度增加法的塑性区贯通判据只能作为边坡失稳的必要非充分条件，与前人研究结果一致。形成贯通的滑移面即为潜在滑动面，可以作为极限平衡计算中的指定滑带。

5 极限平衡法原理

对于圆弧滑动面，简化Bishop法是具有足够精度的工程实用计算方法;对于任意形状滑动面，只有严格按条分法才能得到合理的安全系数，而且条间力函数的变化对安全系数值影响也很小。鉴于此，本文参考极限平衡法中简化Bishop法来推导安全系数计算公式不失一般性。

图2 不同重度幅值下边坡内塑性区发展及分布图

6 算例分析

下面针对边坡来进行计算分析，验证本文方法的可靠性。

工程背景:淹锅沙坝某土坡，坡高35 m，坡长90 m，坡角45°，假三维有限元模型共有6 813个单元，9 026个节点，模型边界采用固定约束，坡体材料选用M-C模型，材料参数见表1。

由于土体所发生的位移远远小于边坡几何尺寸，这是经典有限元常用的小变形假设［9］。在有限元软件计算中打开大变形分析，选择大变形关联流动法则来模拟分析，加载子步100，最大循环50次，收敛容差0．05。通过编程计算得出各自结果见图3～图6。

图3 强度折减FS=1.243塑性区贯通图

图4 极限平衡FS=1.116失稳滑动面图

图7中1号滑弧为有限元重度增加法确定的滑面，2号滑弧为有限元强度折减法确定的滑面，3号滑弧为极限平衡法确定的滑面。由于边坡的最危险滑动面应穿过塑性区上最大塑性应变的峰值点，由此借助AutoCAD寻找滑弧圆心(－4．2，63．85)半径63．987，以及滑弧剪出口坐标(0，0)，(81．341，50)拟合最危险滑动面方程，输入已编制的程序，计算安全系数。从以上三条滑动面位置对比发现最危险滑面位置基本一致，只是由于重度增加法与强度折减法得到的塑性区大小不同，主要是施加荷载的方式不同引起的坡体内部网格单元达到屈服状态的先后顺序不一而导致的。

图5 原重度增加FS=1.432塑性区贯通图

图6 原重度增加法塑性应变等值线图

图7 CAD描图三条滑动面位置对比

表2 不同方法计算安全系数的对比

从表2安全系数上对比分析，综合法与传统极限平衡简化Bishop法的结果最为接近，与强度折减法的结果相比次之，而与原重度增加法的结果相差最大。宏观评价以上各种方法对安全系数精度的影响，极限平衡法是建立在对最危险滑面形状的一系列假定，并视岩土体为刚体的前提之下，通过复杂的内置计算程序搜索最小安全系数对应的滑面，其本身与实际情况就有一定的出入。而有限元法自身的精度影响和收敛容差的设置，会出现安全系数计算误差，而综合法兼顾两者之长，取长避短。表3验证了运用综合法分析的结果较为理想，并得出原重度增加法的分析结果误差较大。

表3 最危险滑面确定的难易度比较

7 结语

1)本文方法不但使滑动面的确定简单明了，且避免了以往传统极限平衡法中对滑动面位置的假设、搜索和安全系数试算等一系列复杂的过程，同时避免了有限元软件中命令流、计算程序等复杂的调试过程以及人为设置带来的偏差。2)有限元重度增加法和极限平衡法的综合法在分析求解安全系数方面思路简洁清晰，且能够建立安全系数与稳定分析的直接联系，具有一定的参考价值。3)有限元重度增加法在计算安全系数方面存在较大的误差，本文分析重度增加法安全系数误差主要来源于重度在滑面方向的投影分力，并因此对重度增加法安全系数计算公式加以改进，减小了误差。4)依然从重度增加法安全系数计算式建立抗滑力增量与下滑力之间的不等式，从而确定该方法的适用范围。

［1］谢永利．大变形固结理论及其有限元法［M］．北京:人民交通出版社，1998．

［2］任青文．非线性有限单元法［M］．南京:河海大学出版社，2003．

［3］唐春安，唐烈先．岩土破裂过程分析RFPA离心加载法［J］．岩土工程学报，2007，29(1):31-32．

［4］邓建辉．CAD技术在边坡稳定性分析中的应用［D］．武汉:中国科学院武汉岩土力学研究所，1992．

［5］林 鹏，卓家寿．边坡稳定性广角度分析的等值面法［J］．水利学报，2000(8):75-79．

［6］王成华，夏绪勇，李广信．基于应力场的土坡临界滑动面的蚂蚁算法搜索技术［J］．岩石力学与工程学报，2003，22(5):813-819．

［7］傅荟璇．MATLAB应用设计［M］．北京:机械工业出版社，2010．

［8］Hibbitt H D，Marcal P V，Rice J R．A finite element formulation for problems of large strain and large displacement［J］．International Journal of Solids and Structure，1970(6):1069-1086．

［9］孙广忠．工程地质理论与实践［M］．北京:科学出版社，1997．