边坡稳定性分析中的有限元极限平衡法＊

郭子仪 范振华 朱云升 张谢东 曾格华

（内蒙古交通设计研究院有限责任公司1） 呼和浩特 010010） （武汉理工大学交通学院2） 武汉 430063）

0 引 言

当前，边坡稳定性分析方法主要有2类：一类为建立在刚体极限平衡理论上的极限平衡法；另一类为以有限元法为代表的数值分析方法．极限平衡法能给出物理意义明确的边坡稳定安全系数以及可能的破坏面．以有限元为代表的数值分析方法不但能考虑坡体变形对其稳定性的影响，而且能通过边坡内应力应变分布分析边坡失稳过程等，但很难给出明确的稳定安全系数以及可能的破坏面．边坡应力应变一般是按未破坏时的边界条件计算出来，而实际上边坡内任何超过抗剪强度或抗拉强度的应力状态都是不能稳定的，一旦发生局域破坏，应力将会重新调整，边坡稳定性安全系数也随之而有所变动．因此，简单采用传统的极限平衡法或者有限元法都不能真实模拟边坡失稳破坏过程中稳定安全系数的变化［1－2］．

本文将传统极限平衡方法和有限元方法相结合，提出有限元极限平衡法用于边坡的稳定性分析，采用有限元方法计算边坡岩土体内的应力应变分布，求出边坡潜在滑动面上的应力，通过极限平衡法的概念，求出该潜在滑动面对应的稳定安全系数，通过边坡潜在滑动面上的下滑力和抗滑力在边坡失稳破坏过程中的变化规律探寻浅变质岩边坡变形破坏机理．

1 边坡失稳破坏有限元模拟原理

边坡发生失稳破坏并不是瞬间发生整体滑动或剧滑失稳破坏，而是一个由局部破坏逐渐扩展以至贯通形成滑面的渐进蠕滑过程．当边坡呈现局域性的集中应力超过材料的强度值时会发生局部破坏，而后发生应力释放、应力转移和应力重新调整．在破坏区域的邻近区域所受到的影响最大，该邻域可能由原先没有超过强度值转变为超过强度值而发生破坏，并进行应力释放，又把多余的载荷转加到其他区域．在不断地发生应力释放、转移和调整过程中，边坡破坏面也不断地延伸，最终将有2种可能性：（1）是破坏面完全贯穿，滑体将在滑床上开始作加速运动；（2）是破坏面没有贯穿，在伸展到某一区域后就停止扩展，其前方区域的应力应变均未超过强度值．前一种情况，显然安全系数小于1，而在后一种情况安全系数大于1．由于应力调整的结果，按最终应力状态计算出的安全系数必然与按初始边界条件进行应力分析计算所得的安全系数不一致，按最终应力状态计算出的安全系数更能反映边坡稳定性的实际状况［3－5］．

模拟边坡失稳破坏过程在数学处理上不易实现，采用边坡滑动面上的法向和切向弹簧模拟滑坡体和滑床之间的接触摩擦问题，通过弹簧刚度变化模拟滑动面上拉应力开裂和岩土体剪应力超过其抗剪强度导致滑动面上应力重分配．将边坡滑动面上切向弹簧力超过岩土体抗剪强度的区域定义为剪切塑性区，法向弹簧为拉力的区域定义为拉裂区，因边坡潜在滑动面上的应力重分配会导致拉裂区和塑性区在边坡失稳破坏过程中逐步扩展，通过边坡潜在滑动面上塑性区的变化发展趋势模拟边坡失稳破坏过程．

2 有限元极限平衡法分析边坡稳定性的实现

采用传统极限平衡方法和有限元分析相结合的方法分析边坡稳定性最关键的问题包括3个方面：（1）浅变质岩边坡潜在滑动面的确定；（2）模拟潜浅变质岩边坡蠕滑过程中潜在滑动面上应力重分配；（3）如何根据潜在滑动面上的应力计算边坡稳定安全系数［6－7］．

2．1 计算边坡的选择

某高速公路JK104＋820～SK104＋840段最大路堑边坡高度为37m，该路堑左侧边坡完成边坡削方施工后，左侧中线105m附近出现山体沿路线前进方向开裂，裂缝左侧方向沿山间凹地近似垂直路线方向延伸，右侧沿与路线近似37°方向向下发展，滑坡界限清楚，钢筋混凝土框架局部断裂，框架整体发生下滑，观测显示裂缝不断向两侧及山顶发展．

根据深层侧向位移监测成果，该滑坡后缘拉裂缝与中线水平距离最长192m，滑体平均厚度21．4m，最厚处达到39．5m，平均宽度为85m，滑体面积约为3万m2，滑体体积约为51万m3，属于大型滑坡．滑动面发育在强风化与中风化的接触面附近，主要物质成分为褐黄色－黄色强风化硅化绢云母板岩，裂隙极发育，其间充填粘土，钻探岩芯随含泥量的不同而呈泥柱状和砂状，泥柱状天然状态呈可塑－软塑状．

2．2 边坡潜在滑动面的确定

边坡潜在滑动面的确定可根据强度折减法结合现场地质勘探资料和深部位移测斜资料拟合分析得出［8－14］．边坡岩土体材料强度系数折减后，经过非线性有限元分析，此时边坡内将出现一塑性区，见图1，塑性应变等值线示意图见图2．

图1 主滑动面塑性区分布图

图2 最危险滑动面示意图

滑坡塑性区主要出现在强风化层和弱风化层交界面上，属于深层滑坡．以计算坡面贯通的塑性区作为潜在滑动面，结合边坡钻探地质资料，以强风化层、软弱夹层为控制点采用拟合法对计算潜在滑动面进行修正，得出合理的潜在滑动面．

2．3 滑动面上应力重分配及稳定安全系数

本文在有限元分析过程中对浅变质岩边坡潜在滑动面上施加切向和法向弹簧模拟滑坡体和滑床之间的接触摩擦问题．假定滑动面上法向弹簧单元出现拉应力即认为该点出现拉裂缝，此时切向弹簧的应力即为该点抗剪强度，计算过程中通过改变弹簧刚度将超过土体抗剪强度的那部分剩余应力平均分配到其它弹簧上，由此模拟滑坡体在滑动过程中潜在滑动面上应力变化．在此过程中，大于最大切向力的弹簧力得到逐步释放，剩余力由其他弹簧承担，通过反复迭代计算即可实现浅变质岩边坡在渐进破坏过程中潜在滑动面上的应力重分配．将边坡滑动面上的塑性区基本贯通时各节点的法向弹簧力和切向弹簧力代入传统极限平衡法边坡稳定安全系数计算公式中计算边坡稳定安全系数，见式（1）．这样，一方面可以考虑边坡渐进破坏过程中内部应力应变状态对边坡稳定性的影响，另一方面可以利用极限平衡理论边坡稳定安全系数明确定义计算边坡稳定安全系数．

式中：Fni和Fτi分别为有限元计算得到的破坏面上第i单元的法向力和切向力；Li为第i单元的长度．

2．4 边坡主滑动面有限元模型建立

以边坡主滑动方向上的剖面作为主滑动剖面建立二维有限元模型，在滑动面每个节点上施加2个弹簧单元，其中1个弹簧单元与该节点处滑动面的切线方向相同，另1个弹簧单元与该节点处滑动面的法线方向相同，可以通过潜在滑动面上相邻两节点的坐标计算出这2节点中任意1个节点的切线方向和法线方向．在边坡潜在滑动面上某节点施加的切向弹簧单元应力模拟该点滑动面上的剪应力大小，法向弹簧单元应力模拟该点潜在滑动面上正应力大小．滑动面每个节点上的2个弹簧单元一端与节点相连，另一端的位移自由度和转动自由度全部施加约束．通过这两个弹簧单元的位移大小可以较好地模拟潜在滑体和滑床之间的相对滑动趋势．

2．5 材料参数和荷载的确定

本文计算材料抗剪强度参数主要根据室内试验结果和参数反演法结果综合取值确定，根据潜在滑动面以上坡体岩土层处于天然状态时不同厚度进行加权平均得出潜在滑动面以上岩土体综合容重、弹性模量和泊松比，边坡荷载仅考虑潜在滑动面以上滑体的自重作用，根据滑体综合容重来计算自重，计算参数见表1．

表1 计算参数表

3 计算结果分析

3．1 边坡应力和变形分析

通过有限元计算可以得出主滑动剖面第一主应力云分布图，合成位移图以及剪应力云图，分别见图3～5．同时，可以得到主滑动面内拉应力区分布和关键点拉应力值，见表2，表中所列拉应力区是指从起点到终点连续出现拉应力的区域．

表2 主滑动坡面拉应力值分布及其高程

图3 主滑动剖面第一主应力云

图4 主滑动剖面合成位移云

图5 主滑动剖面剪应力云

从上述表中边坡在自重作用下的第一主应力计算结果来看，主滑动剖面的表面和潜在潜在滑动面上不同高程处均会出现连续拉应力区，拉应力基本上分布在边坡中后缘部位，由于边坡岩土体抗拉强度往往较低，在这些拉应力作用下，边坡后缘往往会率先产生拉裂缝．潜在滑动面拉裂区平均拉应力均小于坡面拉裂区平均拉应力，表明边坡拉裂区均是由坡面产生，随着滑体滑动逐步向坡体内部扩展．

若遇降水，滑坡土体处于一定程度的饱水状态，除土体自重增加以外，土体的抗剪强度也要下降，滑坡体的稳定性降低，上述拉应力区的任何部位都有出现拉裂缝的可能性．边坡坡体内一旦产生拉裂缝，裂缝会伴随着滑坡体滑动过程逐步扩展，加速滑坡体下滑过程，最终在坡体内形成连续剪切破坏面，导致边坡失稳破坏．边坡上部向下蠕滑过程中对坡脚部分岩土体起到推挤作用，使得坡脚位置承受较大的压应力，因此，边坡下部主要承受压应力作用．

计算剖面最大合成位移出现在人工开挖临空面顶部和坡体后缘之间，人工开挖临空面上出现挤出变形．因此，人工削方临空面直接影响滑体位移走向，临空面过高过陡会导致过大位移产生，坡体除了沿潜在滑动面产生主滑动破坏外，还有可能由于临空面位移过大产生浅表层次生滑动．

从计算剖面剪应力云图分布来看，坡体内最大剪应力出现在潜在滑动面上，沿着潜在滑动面方向边坡的剪切应力逐渐增大，尤其是坡脚部位剪应力较大，当潜在滑动面上剪切应力超过了该点的抗剪强度，则在该点产生剪切破坏，潜在滑动面上出现剪切破坏点后随着滑坡体滑动过程的逐步发展，潜在滑动面上发生应力重分配，破坏点逐步增多形成连续破坏面．

3．2 滑动面受力特征及稳定安全系数分析

考虑边坡潜在滑动面上部分区域剪应力超过岩土体抗剪强度而产生应力重分配的受力特点，进行主滑动面应力重分配的迭代计算，可以得出失稳边坡潜在滑动面上下滑力和抗滑力的变化规律以及边坡稳定安全系数变化规律．

计算剖面潜在滑动面经过每次迭代计算过程均发生了应力重分配，导致潜在滑动面上塑性区随着迭代次数增加而逐步扩展，直至基本贯通，迭代计算过程中，边坡各计算剖面上抗滑力和下滑力也随之发生变化，边坡稳定安全系数逐步减小，滑坡体产生完全失稳破坏．潜在滑动面上抗滑力、下滑力以及稳定安全系数变化见表3．

表3 计算剖面潜在滑动面受力表

通过潜在滑动面上的法向和切向弹簧模拟滑坡体和滑床间的接触摩擦问题，采用本文有限元极限平衡法对边坡进行计算分析可知：

1）计算剖面潜在滑动面上出现明显分段，潜在滑动面上边坡后缘顶部产生拉应力集中，出现“卸荷拉裂段”，边坡前缘坡脚处或者人工开挖临空面上产生剪应力集中或者边坡坡顶后缘由于坡顶拉裂产生蠕滑而出现剪切塑性区，定义为“蠕变剪滑段”，边坡中部应力接近而尚未超过该处岩体抗剪强度，仅有变形而未剪滑，为边坡稳定区，该段对保持斜坡整体稳定起到锁固作用，定义为“锁固段”．随着各计算剖面迭代计算次数增加，潜在滑动面上“卸荷拉裂段”、“蠕变剪滑段”和“锁固段”的分布长度在逐步改变，“剪变蠕滑段”分布长度随着潜在滑动面上应力重分配过程逐步增大，超过岩土体抗剪强度的那部分剪应力逐步向“锁固段”转移，当锁固段岩体内的剪应力积聚到超过了该处岩体的抗剪强度时，“锁固段”则逐步转换成“剪变蠕滑段”，“锁固段”分布长度则逐步减小，当“锁固段”全部变成“剪变蠕滑段”后坡体剪切滑动面贯通，边坡产生整体失稳．

2）边坡初始稳定安全系数均等于1左右，处于临界极限平衡状态，潜在滑动面上均分布有长度不等的剪切塑性区．路堑边坡开挖后，形成了较高的临空面，改变了坡体受力状态，由于边坡处于断层破碎影响带内，强风化层岩石结构松散，风化裂隙发育，岩质较软，坡体上水田内蓄积大量积水，地下水埋藏浅，致使岩土体抗剪强度降低．滑坡体裂缝已发展到宽2～20cm，长3～84m，滑坡周界裂缝仍在不断发展，已基本连通．深层侧向位移监测数据表明，滑坡体上不同部位的13个孔均有位移，滑坡体一直存在滑动速率．

3）从潜在滑动面上应力重分配的重复迭代计算结果可看出，各计算剖面的潜在滑动面上剪切塑性区均随着迭代次数增加而逐步扩展，锁固段长度逐渐减小，当锁固段的剪切应变能赋存到与锁固段的抗剪强度相等时，此时整个坡体处于稳定不平衡状态，滑体沿着潜在滑动面作时而停止时而运动的缓慢移动．在这交替运动过程中，潜在滑动面岩体受到强烈的扰动，使得岩体的抗剪强度骤然下降，当锁固段失稳断裂后与坡体的后缘拉裂段和前缘的滑移段形成贯通面，滑坡体就形成了完整的滑动面，坡体缓慢地向边坡前缘软弱结构面寻找出滑口，前缘临空面的岩体开始与母体脱离，滑向坡下，滑坡体产生完全失稳破坏．

4）随着迭代计算次数的增加，潜在滑动面上下滑力基本保持稳定，抗滑力随着潜在滑动面上剪切塑性区的扩展而逐步减小，稳定安全系数也随之减小，滑坡体逐步产生蠕变滑动．

3．3 计算结果对比分析

按照上述传递系数法的基本原理对王家寨和南约沟滑坡各计算剖面进行稳定安全系数进行计算，结果见表4．

表4 边坡稳定安全系数对比表

从表中可以看出，传递系数法由于没有考虑滑坡体内部和应力状态，计算的主滑动剖面稳定安全系数大于有限元极限平衡法计算结果。该边坡本来属于潜在不稳定体，处于极限平衡状态，如果采用传递系数法计算的稳定安全系数进行稳定性评价，则可以认为是稳定体．采用有限元极限平衡法计算各剖面的稳定安全系数，考虑了边坡潜在滑动面上的应力在边坡蠕变滑动过程中的重分配问题，这种方法计算的稳定安全系数更符合边坡的实际情况．

4 结 论

1）采用弹塑性有限元强度折减和浅变质岩边坡现场地质勘探资料相结合的方法确定浅变质岩风化层边坡潜在滑动面是合理可行的．

2）采用浅变质岩风化层边坡潜在滑动面上的法向和切向弹簧单元模拟滑坡体和滑床之间的接触摩擦问题能够真实反映浅变质岩风化层边坡渐进破坏的基本过程和特点．

3）浅变质岩风化层边坡潜在滑动面在边坡渐进破坏过程中发生应力重分配，坡破坏从边坡后缘拉裂开始，随着潜在滑动面应力重分配过程的发生，潜在滑动面上的剪切塑性区逐步扩展，直至贯通，浅变质岩风化层边坡主要破坏模式为渐进蠕滑拉裂破坏．

4）通过具体工程实例分析可知，本文提出的有限元极限平衡法分析边坡稳定性是一种有效可行的方法．

［1］黄昌乾，丁恩保．边坡工程常用稳定性分析方法［J］．水电站设计，1999，15（1）：45－48．

［2］朱云升，李洪亮．道路滑坡稳定性分析有限元方法研究［J］．公路交通科技，2007，24（4）：55－59．

［3］戴自航，卢才金．边坡失稳机理的力学解释［J］．岩土工程学报，2006，28（10）：104－107．

［4］江学平，于远忠．边坡渐进失稳破坏的概率分析［J］．西南科技大学学报，2007，22（2）：78－82．

［5］LIU Y Q，LI H B，DAI H C，et al．Preliminary study on the progressive failure of a layered rock slope under explosions［J］．Key Engineering Materials，2006，306－308：1461－1466．

［6］朱典文，唐小兵，谢支钢．边坡稳定性分析的数值解研究［J］．武汉理工大学学报：交通科学与工程版，2003，22（4）：856－859．

［7］EBERHARDT E，STEAD D，COGGAN J S．Numerical analysis of initiation and progressive failure in natural rock slopes－the 1991Randa rocksfide［J］．International Journal of RockMechanics and Mining Sciences，2004，41：69－87．

［8］DAWSON E M，ROTH W H，DRESCHER A．Slope stability anaysis by strength reduction ［J］．Geotechnique，1999，49（6）：835－840．

［9］郑颖人，赵尚毅．有限元强度折减法在土坡和岩坡中的应用［J］．岩石力学与工程学报，2004，23（19）：3381－3388．

［10］MATSUI T，SAN K C．Finite element slope stability analysis by shear strength reduction technique［J］．Soils and Foundation，1992，32（1）：59－70．

［11］靳付成．岩质边坡失稳判据及临界滑动面确定方法研究［D］．南京：河海大学，2007．

［12］张金保．推移式滑坡对施工中桥梁的影响以及处理对策［J］．青海交通科技，2013（4）：55－58．

［13］岳尚全．岩体稳定和区域稳定数值模拟边界条件研究方法［J］．岩石力学与工程学季，1999（2）：262－266．

［14］张鲁渝，时卫民，郑颖人．平面应变条件下土坡稳定有限元分析［J］．岩土工程学报，2002（4）：487－490．