郭子钰 刘向峰 王来贵 王 勇2

（1.辽宁工程技术大学力学与工程学院，辽宁 阜新 123000；2.中煤科工集团沈阳设计研究院有限公司，辽宁 沈阳 110015）

近年来，露天采矿规模迅速扩张，开采技术越来越先进，露天采场越来越深，边坡的高度和扰动面积逐渐增大，引发的滑坡事故增多[1-2]，许多露天煤矿边坡[3-4]都经历过变形，当变形量累积到一定程度时，萌生裂纹，裂纹之间相互扩展连接，再当坡体被裂纹贯穿后，就在重力或其他外界因素扰动下发生失稳滑动。由于滑坡成因和诱因十分复杂，边坡的变形、破坏、滑动过程也具有阶段性，处于不同变形阶段的边坡，距离边坡整体破坏的时间间隔也不同，所以滑坡时间预报和空间预报至今还是一个世界性的难题[5]。

目前，曹兰柱等[6]应用数值模拟的方法，模拟了宝日希勒露天煤矿的应力应变和位移变化规律，结果表明南帮边坡失稳时，边坡失稳以剪切作用为主，拉伸破坏相伴发生。王玉凯等[7]采用模型试验对比天然和饱和工况下的排土场变形破坏过程，应用数字散斑和测点追踪等软件，监测位移分布特征及演化规律，结果表明排土场变形破坏形态呈现滑塌—牵引—推移态势。张志飞等[8]采用颗粒流程序（PFC）对反倾层状岩质边坡变形破坏过程进行模拟，得到了边坡岩层的主要变形破坏方式为弯曲变形、折断破坏，变形首先在坡顶发生，破坏从坡脚开始，边坡在变形破坏过程中表现为上部受拉下部受压的特征。张硕等[9]对坡体裂缝发育与演化及变形破坏模式进行模型试验，结果表明坡体裂缝体系演化机制为临空面附近产生张拉裂缝且向后扩展—坡体侧翼产生剪张裂缝—后缘产生贯通性的张拉裂缝。杨天鸿等[10]通过位移监测手段在弱层流变的基础之上研究阜新海州露天矿边坡变形破坏机理，划分出了滑移变形、溃屈变形、隆起变形和剪切变形。

上述学者通过数值模拟、物理模型、位移监测等手段对边坡变形、破坏、滑动进行研究并取得了丰硕的研究成果，但是在实际工程当中，边坡变形、破坏、滑动过程是一个连续不间断的一个灾害演化，因此在前人研究的基础上本项目应用连续—非连续方法（CDEM），以抚顺西露天矿为研究对象，研究露天矿南帮边坡变形、破坏、滑动全过程，从南帮边坡的变形孕育到裂纹破坏，再到最后失稳滑动进行系统的分析和讨论，从而对南帮边坡以及相似工程背景的露天矿边坡提供数据参考。

1 CDEM方法

中科院力学研究所提出了基于连续介质力学的离散元方法，简称CDEM[11-13]。CDEM的数值模型可以分成块体部分、界面部分，如图1所示[14]，真实界面是块体真实存在的不连续面，如断层、节理等，其物理力学参数与块体的参数一致。虚拟界面在块体中并不存在，主要有两个作用：①把相邻块体连接，传递力学信息，②为裂纹的发育扩展提供潜在的备选通道。本研究应用虚拟界面和块体单元来描述边坡裂纹的扩展以及滑动，真实界面转换为块体单元进行计算，在坡面区域将单元细化，最大限度提供裂纹的扩展路径。

CDEM方法的核心计算公式为

式中，ü(t)、u̇(t)和u(t)分别是单元内所有节点的加速度列阵、速度列阵和位移列阵；M、C、K和F（t）分别为单元质量矩阵、阻尼矩阵、刚度矩阵和节点外部荷载列阵。

CDEM采用显式欧拉前向差分法进行问题的求解，整个计算过程中通过不平衡率表征系统受力的平衡程度，计算流程如图2所示。

2 变形破坏滑动机理

根据文献[15]可知，边坡失稳滑动的机理分为变形、破坏以及滑动，在变形阶段满足弹性准则，变形分为法向变形和切向变形；在破坏阶段满足摩尔库伦准则，分为拉张破坏和剪切破坏；在滑动阶段满足动力学方程，分为稳定滑动和非稳定滑动，即负反馈和正反馈系统。

变形阶段：满足τ＜σ·tanϕ+c，其中τ为弱面剪应力，σ为弱面法向应力，c为弱面的黏聚力，ϕ为内摩擦角。

破坏阶段：满足σ3＞σt，则块体发生拉伸破坏，裂纹的方向与最小主应力σ3垂直，其中σ3为最小主应力，σt为块体的抗拉强度；满足σ1≥σ3tan2(π /4+ϕ/2)+2ctan(π /4+ϕ/2)，则考虑剪切破坏，裂纹方向与σ3夹角为π/4+ϕ/2，其中σ1为最大主应力。

2.1 单元变形演化规律

在变形阶段时，边坡单元为线弹性模型，根据变形机理发生法向和切向变形，对应到数值计算中就是节点位置坐标发生更改，本节讨论单元变形的演化过程。结合CDEM方法，将边坡划分单元网格进行迭代计算，单元发生法向和切向变形，如图3所示，从图中可以明显地看出块体单元图3（a）到块体单元图3（b）再到块体单元图3（c）的变化。块体单元图3（a）中共有5个小三角形单元，节点数为6个，块体单元图3（c）三角形单元也有5个，节点数也为6个，但是单元图3（a）的三角形单元开始阶段都是正三角形单元，受到重力作用，导致各个节点受力不同，最终使得单元形状发生改变，变得扁长，最后形成单元图3（c）的形态，块体单元图3（a）的重心向右下方向移动一小段距离。在变形阶段时，边坡内各个单元的平均应力均满足τ＜σ·tanϕ+c，经过法向和切向变形，逐渐累积变形能量。总之，在变形阶段地质体受到重力作用影响，会发生切向和法向变形，边坡地质体会产生微小变形，映射到数值计算中块体单元在变形阶段，单元间是连续的，节点与节点之间不会发生分离，在单元内部都是共享节点，单元的节点由于受到外力作用，导致位置坐标发生更改，单元的形状就会发生微小的改变。

2.2 单元破坏演化规律

如图4所示，在变形阶段计算平衡后即进入到破坏阶段，在单元界面设置接触弹簧（接触面），对应的接触弹簧满足拉破坏准则或剪破坏准则。

满足破坏准则时界面发生断裂，也就是当变形量逐渐增加，达到一定阈值时，单元间界面发生破坏，边界面上的接触弹簧断裂，形成裂纹，由变形到破坏，即单元由连续状态变为非连续状态。由于单元之间所受到的应力不同，导致接触面破坏分成拉张破坏和剪切破坏，对应的裂纹也分成拉张裂纹和剪切裂纹。即单元间发生破坏，裂纹就会沿着单元的边界面进行扩展，即虚拟界面扩展。

在单元变形之后继续发展，当满足拉张破坏条件时，单元的节点会发生分离，单元与单元之间出现缝隙如图4（b），即单元破坏。单元破坏是一个时间点，在这个特定时间点内单元破坏，即接触面满足拉张破坏准则，单元由连续状态变成非连续状态。随着外力的扰动，连续区域内的单元满足变形准则，继续发生法向和切向变形，导致非连续区域单元间的缝隙越来越大，随着单元应力变化，节点的应力也在改变，当相邻单元间再次满足破坏准则时，接触界面又再次发生断裂，就再次萌生裂纹，当裂纹逐渐连接成一条大裂纹时，就标着破坏区域完全贯穿破坏（图4（c））。如果继续发展就有可能发生滑动，就进入到滑动阶段。

图5为沿着单元界面的剪切破坏演化，和拉张破坏一样，当满足剪切破坏条件时，单元发生剪切破坏如图5（b）所示。在图5（b）的左侧开始剪开一个小裂纹，在小裂纹的周围单元是连续状态，满足变形准则，发生法向和切向变形，变形量继续增加，当下一个界面满足破坏条件时，则界面弹簧断裂，形成裂纹，以此类推实现剪切裂纹的扩展，形成图5（c）中的形态。

2.3 单元滑动演化规律

对于滑动阶段则是在破坏的基础上延续，边坡破坏产生的大裂纹上下贯穿坡体，形成滑体，受到重力或其他扰动作用，滑体沿着弱层面向下发生滑动。滑动阶段的单元演化主要是以滑体在破坏路径上整体滑动（即裂纹扩展方向），局部单元产生微小变形，如果再次达到界面破坏准则，则继续开裂形成新的裂纹，新的裂纹继续扩展连接，形成新的离散区域，继续滑动。由于技术原因，在单元中无法监测阻尼力变化，故监测速度变化以等效阻尼力变化。

南帮滑移方向为x正方向，以x正方向为例，式可以看出，下滑力P(t)基本不再改变，如果a增大，阻尼力减小，就可以等效阻尼力随着速度增大而减小，系统为正反馈系统。

图6（a）为局部剪切破坏最终状态，图6（b）和图6（c）为滑动演化的单元形态。在剪切破坏结束时，裂纹完全贯穿破坏区域，由图6（a）到图6（b）在主裂纹区域相邻的界面都已经完全分离，界面上的节点都为非共享节点，受到重力作用，再经过迭代计算，连续区域发生法向和切向变形，非连续区域每个单元的下滑力大于阻尼力，导致沿着裂纹面方向滑动，如图6（c）所示，在滑移的过程中伴随着变形和进一步的破坏。单元滑动演化过程又细分为稳定滑动和非稳定滑动，图6只是滑动演化，没有体现出稳定或非稳定滑动，在实际计算当中，设置监测点监测速度变化来描述稳定或非稳定滑动，判断依据根据速度曲线斜率来判断稳定或非稳定滑动。

3 开挖作用下南帮边坡数值模拟

3.1 南帮滑坡模式分析

国内许多学者对抚顺西露天矿南帮进行滑坡分析[16-18]，结合抚顺西露天矿的工程地质因素等，本研究提出南帮的滑坡模式主要是滑移—剪断模式，如图7所示。南帮现在形成一个大滑体，其潜在滑动面埋深较大且其倾角（30°左右）大于坡角（19°～27°），构成“隐伏型”顺倾坡[19]。

从图7中可知，岩层自上而下为凝灰岩区域、弱层区域、玄武岩区域以及花岗片麻岩区域，南帮滑坡大概分成五部分，分别为拉张变形区域、沉陷变形区域、滑移变形区域、隆起变形区域以及剪切变形区域。边坡顶部区域受到重力作用，使得顶部发生沉降变形差，形成地裂缝，在拉张变形区域，凝灰岩上部区域受到边坡下部区域的拉力，产生拉张变形，使其发生沉陷变形，凝灰岩中部区域沿着弱层向下滑移，在坑底处，受到北帮的影响，滑动产生了抑制作用故产生隆起变形，最后在坑底区域剪出。

3.2 计算工况及模型参数

结合抚顺西露天矿现状，建立E1600南帮计算模型，以左下角为原点，底边长为2 050 m，高为990 m（顶部标高为+126 m），单元为2 827个，节点为1 488个，如图8所示。

南帮地层自上而下分布依次为凝灰岩（分3次开挖）、弱层、玄武岩、弱层和花岗片麻岩。

开挖作用下边界条件：考虑重力场，左右两侧边界x方向固定约束，底部边界x、y方向固定约束，上部为自由边界。南帮从顶部向下开挖分3次，第一次开挖130 m，第二次开挖160 m，第三次开挖140 m。

南帮监测点设置如图8中所示，共计6个。

南帮岩土的物理力学参数根据矿区试验资料确定算例的计算参数，如表1所示。

?

3.3 南帮变形、破坏、滑动规律演化

图9为开挖作用下，随着时步增加，南帮边坡的变形破坏滑动过程，细分为变形压密、变形稳定、破坏、第二次破坏滑移、第三次破坏滑移5个阶段，当满足收敛平衡条件后，就进入到下一个阶段，后续以这5个阶段来描述南帮边坡的变形破坏滑动过程。

3.3.1 位移变化曲线

（1）坡顶区域。由图10可知，坡顶区域x方向和y方向位移的总体变化趋势，靠近坡面位置随着开挖深度增加而逐级增大。在变形压密阶段，由于N1的片麻岩强度高于N2处的弱层，故N2的x方向变形量大，而N1处基本不变；y方向N2处的沉降量也高于N1处的。在第一次开挖变形稳定阶段，x和y方向的数值基本不变，这时主要是积蓄能量，在破坏时释放。进入到破坏阶段，由于产生裂纹，释放能量，导致岩体会产生回弹效应，使得位移曲线产生波动，x方向和y方向的变形量回弹，后来在重力的作用下，继续压实。在第二次开挖时，由于上部岩体被挖空，下部会释放压力，又使得岩体回弹，当第三次开挖时，又反复一次。经过3次的开挖，导致坡顶区域N2处产生x方向2 m，y方向-2.6 m的位移值；N1处x方向基本不改变，y方向-1.6 m。

（2）坡中区域。由图11可知，坡中区域经过3次的开挖发生大规模滑动，坡中区域的N3和N4位移曲线差别很大，曲线间隔为滑动距离。在变形压密阶段，x方向基本不动，主要是y方向发生沉降变形，此时积蓄能量；进入第一次开挖稳定变形阶段，x和y方向在原来的基础上稳定不变，当到达破坏时释放能量；在破坏阶段，有回弹现象，由于在坡中位置，故第一次开挖，对N3和N4的影响不大，基本和坡顶处的位移一致；在第二次开挖时，坡中位置完全临空，N4更靠近坡面位置，这时N4的x和y方向位移迅速增加；到第三次开挖时，滑动得就更加明显。

（3）坡底区域。由图12可知，在坡底区域经过3次开挖，坡底区域N5和N6的x方向位移曲线差别很大，曲线间隔为滑动距离。其中N6为坡底，在整个过程x方向基本不动，N5为坡趾处，到第三次开挖才对其产生影响，x方向位移迅速增长；y方向的变化情况：在变形压密阶段，发生沉降变形，经过开挖，底部产生回弹效应，由于上部和中部的滑体经过3次的开挖已经产生滑动，使坡底处有大量的剪切裂纹，坡底也隆起形成底鼓。

3.3.2 速度变化曲线

坡中区域位移量最大，故监测N4点的x方向和y方向速度变化，来描述整体边坡在开挖情况下的滑动情况。本研究通过速度曲线斜率来描述边坡的滑动特征。算例的一时步需要的时间为1×10-4s，共计5.6 s。

由图13（a）可知N4点x方向的速度变化情况。在压密阶段，寻求计算平衡状态，故速度波动较大。第一次开挖稳定阶段，速度的变化基本不大，这时在积蓄变形能量，进入到破坏阶段，由于裂纹萌生，释放能量，速度增加，阻尼力也增加，增加到一定值时曲线斜率＜0，速度又逐渐降低，属于负反馈系统。到第二次开挖时，岩体挖空，导致N4区域完全临空，速度又迅速增加，由于阻尼力的影响，最后速度不变，这时x方向处于临界状态。在进入到第三次开挖阶段，滑体开始并未整体滑移，N4的x方向速度逐渐降低，随着时间推移，裂纹越来越多，上下部完全贯穿，滑面形成，曲线斜率＞0，最后速度迅速增大，演化成滑坡灾害，这时滑坡属于正反馈系统，边坡失稳，最后x方向位移12.5 m。

由图13（b）可知N4点y方向的速度变化情况。在变形压密阶段，主要是发生沉陷变形，y方向的速度变化较大；压密阶段结束后，速度基本不再变化，一直到第一次开挖稳定，随后进入到破坏阶段。由于裂纹的影响，裂纹萌生释放变形能，岩体向上回弹，y方向的速度向上也迅速增加，而受到阻尼力的影响速度又降低，这时在y方向也是属于负反馈系统。第二次开挖时，由于裂纹继续延伸扩展，继续释放能量，岩体又一次回弹一部分，速度也跟着增加，随后在重力和阻尼力的作用下，速度又向下增加，此时还是属于负反馈系统，阻尼力一直阻碍速度变化。进入到第三次开挖时，由于开始产生岩体裂纹并完全贯穿，故速度逐渐减低，但是随着裂纹增多，破裂区域逐渐连接，滑面形成，出现滑坡灾害，速度就一直增大，边坡失稳，最后y方向的位移为-11.5 m。

3.3.3 位移变化云图

分析不同时段的云图变化来描述灾变过程，云图选取12 000时步、16 000时步、26 000时步、36 000时步和56 000时步分别对应变形压密阶段最后、第一次开挖稳定阶段最后、破坏阶段最后、第二次开挖破坏滑移阶段最后和第三次开挖破坏滑移阶段最后，如图14所示。

图14（a）为变形压密最后阶段，整体来看，整个地层处于压实状态，由于地层的岩层不同，故在凝灰岩区域x方向变形较大，变形值为0.971 m。图14（b）为第一次开挖稳定状态最后，此时经过开挖作用，变形已经基本稳定，把地层的凝灰岩区域挖空，使得玄武岩临空，故在坡面位置x方向变形较大，变形值为0.677 m。图14（c）为破坏阶段最后，可以明显地看出，裂纹布满坡面位置，尤其边坡顶部处有大量的拉张裂纹，并向下延伸，但并没贯穿边坡，x方向位移逐渐向坡中区域靠拢，这是坡中由凝灰岩层和弱层构成，力学参数较弱，故变形值较大，变形值为2.147 m。图14（d）为第二次开挖滑移结束阶段，明显地看出，裂纹在破坏阶段的基础上继续扩展连接，形成新的大裂纹，裂纹也逐渐的向坡底处延伸，从图中也可以看出坡底形成大量的剪切裂纹，整体的位移情况也在延续破坏阶段的位移，并逐渐形成滑体轮廓，在此阶段坡中的位移值仍是最大，位移值为5.257 m。图14（e）为第三次开挖滑移最后，可以看出整个滑体沿着弱层面向下滑移，坡顶、坡中和坡底已经布满裂纹，最后边坡失稳滑动，最大位移值为16.04 m还是在坡中，这是由于在第二次开挖结束时裂纹已经基本贯穿边坡，经第三次开挖，使得边坡中下部临空，没有岩体阻碍滑动，形成滑体后，迅速向中下滑动，但是由于北帮和底部岩层的影响，使其坡中区域向下滑移在底部推挤，底部也形成底鼓。

3.3.4 塑性区变化云图

分析不同时段塑性区云图来描述灾变过程。故云图选取16 000时步、26 000时步、36 000时步和56 000时步分别对应第一次开挖稳定阶段最后、破坏阶段最后、第二次开挖破坏滑移阶段最后和第三次开挖破坏滑移阶段最后。如图15边坡整体塑性区云图所示。图中色带的不同值分别表示边坡的破坏状态：0表示坡体未破坏；1表示当前拉坏；2表示当前剪坏；4表示过去拉坏；8表示过去剪坏。

图15（a）为第一次开挖变形稳定最后，塑性区已经出现，并沿着弱层区域向下扩展，在坡顶区域出现拉张破坏塑性区，但无裂纹产生，当前剪切破坏塑性区沿着弱层向下。图15（b）为破坏阶段最后，可以看出裂纹已经在剪切和拉张塑性区处产生，裂纹也随着塑性区的扩展而向下扩展。图15（c）为第二次开挖破坏滑移最后，可以看出当前剪切破坏基本贯穿整个边坡，塑性区面积也增加明显，裂纹也布满坡面位置，坡底有大量的剪切裂纹，坡顶有大量的拉张裂纹，在第二次开挖时边坡就已经基本破坏，但是由于下部还有岩层的存在，故没有发生大规模滑移，如果再次开挖，裂纹完全贯穿边坡就会发生滑坡灾害。图15（d）为第三次开挖破坏滑移最后，看出边坡沿着弱层面向下滑移，滑体中部的滑移量最大，整个滑体最后推挤在坡底处，而坡底由于受到剪切裂纹和上部岩体的挤压也形成底鼓。为了更直观地分析开挖过程边坡地表的破裂程度，定义地表破裂度为

式中，α为地表破裂度，S1为地表的塑性区面积，S2为地表的整体面积。通过对塑性区云图的处理，将破裂度画成曲线表示，如图16所示。

由图16可知，边坡坡面的破裂度随着时步的增加而增加。在变形压密阶段，此时为弹性阶段，故没有塑性区产生。但是在第一次开挖变形稳定时，塑性区开始出现，这就表明岩体开始出现塑性区域，但无裂纹产生，破裂度达到0.36。进入到破坏阶段，破裂度再次上升，达到0.5，这时边坡已经破坏一半，但是裂纹并没有贯穿边坡。随后进入到第二次开挖，破裂度再次攀升，达到0.8，这时坡面已经几乎完全破坏，坡面也已经被裂纹基本贯穿，此时边坡非常危险。再次开挖，当裂纹完全贯穿，就有很大概率发生失稳滑动，进入到第三次开挖后，破裂度到达0.9以上，边坡发生失稳。

4 结论

（1）南帮边坡在变形压密阶段到第三次开挖阶段前，均属于负反馈系统，属于稳定滑动。第三次开挖阶段之后，顶部裂纹继续增多，但不发生大规模滑移，但中部弱层和底部区域发生失稳滑动，中部和底部属于正反馈系统，滑体沿着弱层面向下滑移，最后推挤在坡底，坡底受到上部的压力，再加上剪切裂纹的存在和北帮的阻碍，使得坡底形成底鼓，南帮的滑坡模式为滑移—剪切式。

（2）变形阶段，边坡整体向下压密变形，越靠近坡面，x方向和y方向的变形位移越大；破坏阶段，边坡顶部受拉，萌生拉张裂纹，沿着弱层面扩展，坡底受剪，随着开采深度增加，裂纹数目增多，裂纹与裂纹之间相互连接，坡顶和坡底裂纹逐渐连接，最终贯穿边坡模型；滑动阶段，滑体沿着裂纹面发生失稳滑动，最后在坡底剪出。

（3）随着开挖深度增加，南帮弱层出露越多，边坡的临空区域增多，引发滑坡几率越高，当开挖深度达到标高为-164 m后，南帮边坡出现大规模滑坡。