用有限元强度折减法评价节理岩质边坡稳定性

徐国良，张振飞，裴伦培，管宏梓，付厚起，董方鹏

(山东省地矿局第一地质大队(山东省第一地质矿产勘查院)，山东 济南 250000)

0 引言

边坡稳定性评价是地质灾害防治工作中的重要研究内容，其评价结果合理与否影响防治工作部署，事关人民群众生命财产安全。

常用的边坡稳定性计算方法主要为基于刚体极限平衡理论的极限平衡法[1]，该方法经过长时间发展，又根据不同的条分方式、条间力假定衍生出多种分析方法，像简单(瑞典)条分法、毕肖普法、简布法、斯宾塞法、Morgenstern-Price法、陈祖煜的通用条分法等，核心原理均是将岩土体视为刚体，根据静力平衡原理、滑体抗滑力与下滑力的比值来评价边坡稳定性。极限平衡法能给出边坡稳定性系数这一直观的评价指标[2]，但是需要事先假定滑动面的形状和位置，并且只考虑在滑动面上的极限平衡状况，而忽略岩体自身变形对边坡所造成的影响，此外分析时做了部分条件假定和简化，在边界条件和计算参数的选取方面也存在争议，因此计算精度受到一定影响。

近年来，有限元数值计算逐渐应用于边坡稳定性评价。与极限平衡法相比，有限元法理论体系更为严格，无需预先假定滑动面的形状及位置[3]，通过前处理建模，可模拟各种岩体类型以及其中不同产状、特性的软弱结构面，可模拟各种复杂的应力—应变本构关系，迅速用图形表示计算结果等。不少学者对有限元法分析边坡问题作了有益探讨[4-11]，但是已有的研究倾向于将有限元法作为传统极限平衡法的辅助论证手段使用，对节理岩质边坡稳定性分析关注不多，且对有限元建模过程论述不细，本文结合一个具体算例，对有限元模拟在节理岩质边坡稳定性评价中的应用作系统论述，讨论有限元模拟在边坡稳定性评价中的积极意义。

1 基本原理

边坡有限元稳定性评价有两个重要分支：有限元极限平衡法和有限元强度折减法，前者从传统极限平衡法演变而来，在某些情况下，滑动面上剪应力方向可能与滑动方向趋势不一致，使得安全系数计算值偏小[12]，而后者则可规避这一问题，并且保持了有限元法在模拟复杂问题上的优点，因此本文选取后者即有限元强度折减法作为边坡稳定评价的核心算法。

有限元强度折减法定义了一个强度折减系数，即边坡内岩土体在外部不变的荷载下所提供的抗剪强度最大值与在边坡内由外荷载产生的实际剪应力的比值[13]。外部荷载在极限状况下产生的剪应力与岩土体抵御荷载所施加的最低抗剪强度应当相等，当假设边坡内部所有岩土体单元抗剪强度的施展程度一致时，抗剪强度折减系数F'就等同于传统定义的整体边坡稳定安全系数Fs。在使用计算机进行计算时，通过不断降低岩体或结构面单元的黏聚力及内摩擦角，使模拟对象达到失稳状态，此时折减系数与极限平衡法给出的边坡稳定性安全系数在概念上是一样的，并且同时可得到潜在滑动面。

折减后的强度参数用下式表示：

式中：c'和φ'分别为岩体初始黏聚力和内摩擦角；cm和φm分别为岩体折减后的黏聚力和内摩擦角，F'为强度折减系数。

本文在进行有限元计算时采用知名工程模拟有限元软件ABAQUS。根据强度折减法原理，在ABAQUS中实现强度折减过程为，将强度折减系数F'定义为一个场变量；定义材料模型参数；指定场变量大小；对模型加载自重影响，建立天然应力平衡状态；线性增加F'，当数值不收敛时结束计算，按照边坡失稳标准确定安全系数。

2 工程算例

2.1 节理岩质边坡概况

某露天开采建筑用花岗岩矿在常年开采后于采坑西北侧形成了一坡面走向约65°的高陡岩质边坡，坡面倾向SEE，边坡水平方向长度约430m，垂向高度78m，总体坡角约75°。闭坑后坡面裸露岩石受风化作用影响出现多次局部滑塌。选取走向155°的典型剖面进行分析，根据工程地质资料，边坡岩土从上到下分别为残积层砂质黏土、花岗岩(图1)。

残积层砂质黏土为黄褐色，为基岩风化残积而成，以粉黏粒为主，含石英砂砾，稍湿—湿，硬—可塑，遇水易软化。厚度0.05～0.1m，平均0.08m。

花岗岩未风化，新鲜面灰白色，中粒结构，块状构造。造岩矿物主要为石英(30%)、钾长石(50%)、斜长石(15%)、黑云母(5%)，矿物粒径在3～5mm。结构构造、矿物成分及色泽未发生变化；坡面岩石微风化；新鲜面岩石未风化，锤击声清脆、有回弹，属坚硬岩。

1—花岗岩；2—节理组；3—产状图1 节理岩质边坡地质剖面图

矿区内无较大断裂及破碎带发育，岩体中的结构面主要为节理，由于陡崖常年临空，卸荷作用明显，使得节理对岩体破坏的主控作用显得较为明显。根据实测工程地质资料，主要节理发育情况如下：节理组Ⅰ产状172°∠27°；延伸较长，贯通边坡，从上到下集中发育3处，每处节理组真厚度为0.5～0.6m，节理密度7条/m，单条节理宽3～8mm，裂面大多平直，旁侧发育羽裂；内有少量充填物，遇水不变质且有一定黏结力；两壁无明显位移。节理组Ⅱ产状168°∠74°，性质与节理组Ⅰ一致，为同一应力作用产生的节理系。

在边坡安全处避开节理组间隔采取5个新鲜完整的较大尺寸岩块试件进行二次取样、加工、磨平。节理组按含节理岩体进行取样，在节理组地表出露处清理场地露出新鲜节理组，厘清节理组上下盘，旋转好方位后搭设结实支架，用小型取心钻机切穿节理组上下盘进行取样，沿节理走向间隔取样5件。为防止试件开裂、震碎，加热PVC管后靠其冷缩箍紧试件，带管进行加工。通过三轴应力试验，分别得到岩块和节理组在不同应力状态下的强度和变形特性，绘制强度包络线和应变关系曲线，取得抗剪断峰值强度、泊松比、弹性模量等参数(表1)。

表1 岩块与节理组力学参数

2.2 传统极限平衡法计算

首先采用极限平衡法中的条分法对算例的岩质边坡进行定量计算。根据《滑坡防治工程勘查规范》(GBT32864—2016)中滑动面为折线的边坡公式进行计算：

式中：

Rn=(Wn((1-rU)cosαn-Asinαn)-RDn)tgφn+CnLn

Tn=Wn(Asinαn+Acosαn)+TDn

式中：F—稳定性安全系数；ψi—第i块段的剩余下滑力传递至第i+1块段时的传递系数(j=i)。Wi—第i条块的重量(kN/m)；Li—第i条块的滑面长度(m)；αi—第i条块的滑面倾角；ci—第i条块的黏聚力(kPa)；φi—第i条块的内摩擦角(°)；rU—孔隙压力比；A—地震加速度；TDi—孔隙水压力产生的平行滑面压力；RDi—孔隙水压力产生的垂直滑面压力。

本文不考虑岩体中水的渗流作用影响。使用Geostudio中的极限平衡法模块Slope/W进行计算得F=1.26。

2.3 ABAQUS有限元模拟分析

2.3.1 算例节理组的模拟

结构面强度参数较岩块低很多，因此对岩质边坡的稳定性来说，起控制作用的是结构面强度，而不是岩块[14]。岩体结构非常复杂，在建立模型时完全细致的反映全部结构面特征是没有必要的，工程实践证明，稳定性分析时只需考虑对岩体起控制作用的2～3组结构面便可，算例中的节理组是与坡面同向的贯通性结构面，属可能导致边坡失稳、起重要控制作用的结构面。有学者采用无厚度硬接触单元模拟结构面[10]，如果结构面两侧实体接触紧密且未发生明显位移，比如剪性节理，使用该单元则效果较好；如果结构面在一定范围内发育密集，且结构面两壁裂开一定距离，内部有低强度填隙物，使用该单元模拟将失真。本算例中由于节理组与岩块强度差异较大，且节理组在厚度方向具有一定规模，因此在ABAQUS中采用具有一定厚度的、低强度实体单元模拟，按连续介质处理，节理组与岩块之间创建相互接触作用，通过不断折减结构面强度参数，以模拟达到失稳状态。

2.3.2 弹塑性本构模型

边坡岩体应力-应变本构关系宜使用弹塑性分析模型[15]。对边坡岩体，当进入塑性状态时，其非线性和塑性变形具有不可恢复的特性，应力增量和应变增量之间的关系用下式描述：

式中：{dσ}和{dε}分别为应力增量和应变增量，[De]为弹性矩阵，[Dp]为塑性矩阵，[Dep]为弹塑性矩阵，g为塑性势函数，f为屈服函数。

2.3.3 屈服准则及流动法则

屈服准则是描述岩体中应力质点由弹性开始进入塑性状态时各分量之间关系的数学表达式，其合理选用是岩坡稳定性分析过程的重要环节，直接影响分析的可靠性和计算精度。长期工程实践表明，Mohr-coulomb屈服准则(莫尔-库伦准则)能够较好的描述岩土材料的破坏行为，因此算例中的岩块及节理均采用该准则：

式中：I1—应力张量的第一不变量；J2—应力偏张量的第二不变量；θσ—剪裂角；c—黏聚力；φ—内摩擦角。

由于算例为无约束的高陡岩质边坡，在进行稳定性计算时，坡体的体积变形所造成的影响并不明显，因此算例将忽略剪胀角，即认为有限元数值计算满足非关联流动法则。

2.3.4 失稳判别准则

边坡失稳判据是有限元强度折减法的关键，目前有限元分析中边坡失稳判据主要有以下2种：

(1)折减后岩体强度参数在规定的迭代次数内不再收敛[16]。

(2)存在从坡面贯通至坡顶的塑性区[17]。郑颖人等研究表明有限元数值计算不收敛时必然表示塑性区发生贯通或者位移发生突变[18]，而塑性区贯通则只是边坡失稳的必要条件，塑性区贯通后，需进一步观察边坡质点的位移情况。本文将采用迭代不收敛结合塑性区贯通并出现明显位移两种判据来分析边坡失稳情况。

2.4 ABAQUS有限元建模过程

在ABAQUS软件中对实体建模的步骤为：创建部件→创建截面材料与截面特征→装配部件→定义分析步→建立相互作用→定义荷载及边界条件→网格划分→提交任务→可视化操作。

2.4.1 创建部件

在部件模块中，进入图形编辑界面，按实测剖面进行边坡几何轮廓的绘制，也可直接导入由CAD创建的*DXF文件。

边坡轮廓绘制完成后，按实测情况建立节理组与岩块的相对位置(图2)。本文算例存在两组与坡面同向的节理组，在剖面中交汇后将边坡分为八个岩块，分别对节理组和岩块建立集合，目的是方便后续对部件属性进行统一赋值。

图2 岩块与节理部件图

2.4.2 创建截面材料与截面特征

在属性模块，分别建立节理与岩块的材料属性。岩体的破坏是一个逐渐发生的过程，由最初的弹性状态过渡到塑性流动，最终达到极限破坏状态，因此岩体兼具弹性与塑性特征。在弹性设置中分别输入节理组与岩块的弹性模量和泊松比，在莫尔-库伦塑性设置中需要指定两种材料的黏聚力和内摩擦角随场变量的变化情况。进入莫尔-库伦塑性命令框，将场变量数量设为1，再在塑性部分依据强度折减公式设置随场变量变化的内摩擦角和剪胀角，场变量此时的意义即为强度折减系数，初始取值一般小于1，目的是为放大强度，便于搜寻极限破坏状态。本算例将初始值设为0.5，采取线性增加的方式分段直线模拟，增量为0.25，终值为3[19]。前文已述，不考虑坡体的体积变形，因此剪胀角取值为零。同理，依据强度折减公式设置随场变量变化的黏聚力。特别的，ABAQUS中不指定单位，用户只需保证所选单位协调一致即可，如采用国际单位制。

2.4.3 装配部件

在该模块中，将节理组与各岩块进行装配，至此，先前分别操作的各部件统一到共同界面。

2.4.4 定义分析步

创建通用静力分析步，初始增量步由1改为0.1，目的是通过减小初始增量步，使迭代次数增加，以达到收敛目的。在方程求解器中将矩阵存储设为非对称分析，选择非对称分析的原因是ABAQUS中的莫尔-库伦本构模型使用的是非关联流动法则，屈服面和塑性势面是不同的，因此刚度矩阵不对称，要用非对称器进行求解。

2.4.5 建立相互作用

首先定义各个接触表面。在管理器中定义8处主面(图3标红部分)，之后定义从面。主从面的定义一般遵循：网格精细的作为从面、刚度较小的作为从面，算例中将岩块表面作为主面，节理面作为从面。

图3 所定义的主面示意图

在该模块中，执行Instance&Create(实例&创建)，将节理组与各岩块进行装配，至此，先前分别操作的各部件统一到共同操作界面。

其次定义带摩擦的接触属性。在接触属性选项里，将力学中的切向和法向行为进行加载，法向行为设为硬接触；对切向行为，在静摩擦-动摩擦系数选项中进行设置。

最后定义接触对。将分析步设为初始状态，表示接触自初始状态就已存在，再选择主面，后选择从面，分别建立岩块与节理组之间的接触对，对应8组主从面，算例中共建立8组接触对。

2.4.6 定义荷载及边界条件

将初始状态分析步类别设为力学、类型为位移/转角，然后选定算例中除坡面外的左右边界后，在弹出的边界条件编辑中勾选U1，目的是限定本实例沿X轴方向的平移自由度。同理，再次创建边界条件，选择底部边界后，在弹出的边界条件编辑中勾选U1，U2，以限定本实例沿X，Y轴方向的平移自由度。在节理、岩块的分量2中各输入-19.36，-26.71(容重值)，以此来模拟重力荷载(图4)。

图4 边界控制及荷载图

2.4.7 网格划分

在网格模块中对部件进行网格划分，以将连续的岩块与节理系离散化为有限数量的单元体。算例对节理及岩块均采用三结点线性平面应变三角形单元(图5)，对全局种子设置时，节理系的种子密度约为岩块的2倍。岩块共划分计算网格906个，节理系共划分计算网格1978个。

图5 岩块与节理有限数量网格图

2.4.8 控制场变量变化

在左侧模型树中对模型的关键字进行编辑，以控制场变量变化。联合使用*Initial conditions与*Field命令[19]，在前者命令中设置type=field定义场变量初始值。定位第一个分析步语句“*Step, name=load, nlgeom=NO, unsymm=YES”，在其前面、点划线之后输入：

*Initial conditions,type=field,variable=1

slope-1.set-1,0.5

定位第二个分析步语句“*Step, name=reduce, nlgeom=NO, unsymm=YES”，在其后输入以下命令：

*Field,VARIABLE=1

slope-1.set-1,3

其中，“slope-1.set-1”为折减部分的点集合名称，算例中该集合为节理组的点集合，从模型树的装配模块中可手动定位点集合的名称。

2.4.9 提交任务

在任务模块中提交计算。

(1)安全系数

进行场变量输出。本文算例中，折减分析步在时间为0.362494时便无法收敛，计算终止，这表明边坡强度在折减到一定程度后到达失稳状态。

将场变量与X方向位移作为输出变量。切换到单元与节点选项，将坡面右上角的顶点作为分析位移的质点，调用Combine函数，绘制FV1与U1关系图(图6)。

图6 FV1随U1(场变量随X方向位移)变化关系图

从图6可以看出，以数值计算不收敛作为该边坡失稳的评价标准，对应的FV1为1.41，即边坡稳定性安全系数Fs=1.41；以坡顶质点出现明显位移作为评价标准，对应的FV1为1.32，即边坡稳定性安全系数Fs=1.32，两种评价标准均表示边坡目前处于稳定状态。与极限平衡法求得的安全系数1.26相比，表明有限元法得到的数据偏安全。

(2)塑性区变化情况

在场输出命令框上方可以看到不收敛前运算的帧数，本算例共形成88帧，逐帧演示可以看到塑性状态在节理系中的演变过程，将其形成PEMAG积分点的等效塑性应变图(图7)，为了突出显示塑性应变效果，隐去边坡网格。从图看以看出，T=0.1时，缓倾的第二组节理在临坡处出现屈服现象，然后一直向上延伸；T=0.2594时，陡倾的节理组在上部开始出现屈服；T=0.3245时，两组节理共同作用，在整个岩坡首次实现贯通；T=0.3625时，节理系区域塑性加强并造成岩坡失稳。图7d为计算终止时塑性应变图，从该图可以很清楚识别出潜在滑动面位置，其形态表现为折线。

图7 积分等效塑性应变区演化图

极限平衡法分析边坡稳定性的时候，需事先假定滑动面的位置和形状，然后再进行计算，而有限元模拟通过对岩块和节理组进行离散化，用数值分析的方法求解其应力应变，得到贯通的塑性应变区，所指示的潜在滑动面直观可视并且物理意义明确。

3 结论

(1)通过具体算例，利用强度折减法对节理岩质边坡在ABAQUS中实现稳定性分析进行了论述，并提供了详细的建模过程。

(2)通过分析坡顶典型质点场变量与X方向的位移关系，以数值计算不收敛作为边坡失稳评价标准的稳定性安全系数为1.41，以坡顶质点出现明显位移作为失稳标准的稳定性安全系数为1.32；极限平衡法所求取的安全系数为1.26；两种方法所求取的安全系数较为接近，其中极限平衡法安全系数略小的原因是其在理论上没有考虑节理与岩石内部的应力、应变关系，未将边坡破坏的发展过程纳入计算，所以计算结果偏小，有限元法在理论上更加符合实际情况。

(3)通过积分等效塑性应变区演化图，直观的得到了边坡内部塑性应变情况，并识别出了潜在滑动面。

(4)对结构面更为复杂的岩质边坡，为最大程度模拟边坡真实情况，在建立有限元模型时所采用的参数，包括计算范围、边界条件、网格划分及收敛标准等需进行进一步反演研究，其与极限平衡法计算结果的比较也需进一步研究。