0 引言

近年来，随着矿产资源的不断开发，开采的技术越来越先进，露天采场也越来越深，边坡的高度和扰动面积逐渐增大，故而引发的滑坡事故增多，而且滑坡规模增大，破坏力强，严重威胁人类社会的安全，因此，边坡的稳定性越来越成为中国矿山防灾的重点之一。含弱层顺倾岩质边坡受地质环境、地形地貌、岩层结构等方面存在较大影响，导致滑动模式种类繁多，变形破坏机理越发复杂。基于长期的现场工程实践发现，在云南省鹤庆县等地多发育含多弱层顺倾岩质边坡，在凝灰岩或花岗岩等硬岩中夹有一些性质软弱、有一定厚度的软弱带或软弱结构面，如泥岩或泥灰岩，软弱夹层强度较低，遇水易崩解，在地质构造活动或人为因素扰动下，极易发生岩间错动，从而演化为沿软弱夹层的斜坡崩滑失稳。因此深入研究充水作用下含软弱夹层顺倾岩质边坡的变形、滑动破坏演化全过程，为以后预测和防治非煤露天矿山滑坡提供理论基础和数据参考。

谭龙金对缓倾顺层岩质边坡中的地下水作用进行分析，推导出了边坡稳定系数和临界充水高度的表达式，建立了边坡滑移破坏判据（谭龙金等， 2014）。霍逸康等（2023）基于UDEC软件，设计了不同边坡角、坡高、软硬层厚度比及结构面特征等 32 组正交实验，将边坡安全系数及相对位移作为敏感性分析指标，发现坡角和岩层倾角对稳定系数影响较大，结构面力学参数和坡高对边坡位移影响较大。江巍等（2023）基于连续-非连续法，模拟坡脚岩体劣化条件下的边坡变形破坏过程，得到了坡脚岩体劣化条件下软硬互层反倾岩质边坡两种破坏机制。宋洋等（2023）开展多种力学模型试验，分析了不同降雨强度、不同岩层倾角及有无支护条件下边坡岩体内部的物理量发展模式及变化特征。许明等（2023）研究了多维多参数地震动工况下层状岩质边坡的地震动力响应和失稳破坏模式，发现地质结构和地震动特性对边坡破坏模式起到决定作用。王来贵等（2023）以抚顺西露天矿为研究对象，采用 GDEM 软件模拟了重力作用下含软弱夹层边坡的变形破坏规律，发现了边坡整体破坏模式为拉裂-滑移-剪断变形。王志颖等（2023）开展大量的含软弱夹层顺层岩质边坡模型试验，并结合数值模拟试验结果，从多种角度分析了地震作用下含软弱夹层顺层岩质斜坡动力响应规律（宋健等，2023；信春雷等，2023）。

由以上研究可知，目前国内外众多学者通过模型试验、数值模拟、理论推导及现场监测等方面对含软弱夹层岩质边坡开展了大量的科研工作，并取得了丰富的研究成果，但在实际工程当中，边坡的稳定性状态只是其滑动、变形、破坏整体演化过程中的一个阶段，需从整体进行全过程分析。因此，本文基于连续-非连续数值模拟软件 GDEM 对某露天矿含软弱夹层顺倾岩质边坡在不同充水状态下边坡的变形演化过程进行模拟研究，研究结果为采后矿区闭坑充水治理及矿坑水蓄能循环利用提供了依据。

1 某露天矿山基本概况

某露天矿区位于滇西北鹤庆县城南部 172°方向，平距 47.5 km，隶属大理白族自治州鹤庆县西邑镇北衙村公所（图1）。矿区地处马鞍山山脉东侧，锅厂河南西部，北衙向斜轴部，地势较低，属区域水文地质单元的径流—排泄区。西部为近南北走向的马鞍山脉，构成澜沧江与金沙江的区域地表分水岭；中部为北衙山间盆地，构成矿区主体；东侧为近南北走向的桅杆坡—笔架山—锅盖条带状山地，构成北衙山间盆地与锅厂河的次级地表分水岭。区内四面环山，中部低凹平坦，自南向北缓倾，西部山体高陡。总体地形西高东低，北高南低，矿区最高点为焦石洞北营村北东侧山坡，标高约 2490 m，最低点为北东部边缘的锅厂河河谷，标高约 1640 m，最大相对高差 850 m，属构造剥蚀—溶蚀高中山地貌。南东部锅厂河 1640 m 标高可视为矿区最低侵蚀基准面标高。现露采矿坑坑底标高已达 1541 m，已采至北衙组（T₂b）岩溶含水层水位以下，保有的资源量、储量主要分布于矿区最低侵蚀基准面之下，地形有利于地表水、地下水汇集，矿坑水不能自流排放。

在露天采场西帮钻探所达深度范围内，场地岩土体主要分为 3 个主层，其中边坡基底岩层为花岗片麻岩，基岩之上是玄武岩风化作用后的太古风化壳（软弱结构面），厚度在 15～20 m 之间，软弱结构面之上为玄武岩层，其夹杂着中风化凝灰岩及未风化花岗岩，厚度在80～200 m范围内，地表岩层主要为强、中风化凝灰岩，广泛分布于露天采场东西边帮，图2 为露天采场西帮所揭露出软弱夹层带范围内的工程地质剖面图。

2 连续-非连续方法（CDEM）

2.1 基本原理

连续-非连续方法是基于连续介质力学的离散元方法，基于Lagrangian能量系统建立控制方程，通过动态松弛法迭代求解的动态显式解算法。该方法将离散元与有限元方法进行深度融合，用于描述材料的渐进破坏全过程。在材料内部进行有限元计算，当材料单元达到相应的破坏准则（最大拉应力准则和 Mohr-Coulomb准则）时，材料单元发生断裂，生成两个有限单元，在材料边界进行离散元计算，分析各单元间的切向弹簧和法向弹簧的断裂情况，从而实现连续性非连续的统一描述。因此，连续-非连续方法可以模拟材料变形、断裂、破坏的全过程，也可以描述材料内部裂纹萌生、扩展、贯通的全过程。

连续-非连续方法的数值模型主要由块体和界面两个部分组成，块体用于描述材料的连续变形特征，由一个以上的单元组成。界面包含真实和虚拟界面，真实界面分为块体边界、结构面及断层等；虚拟界面是指相邻有限单元间的接触面，通过接触面上的弹簧拉伸断裂来描述材料破裂过程的不连续性。具体如图3所示。

a—真实界面；b—虚拟界面

2.2 CDEM的基本方程

基于 Lagrangian 能量系统建立控制方程，拉格朗日能量系统的方程表达式为：

式（1）中，L 为拉格朗日函数，Q_i 为系统非保守力（N），写成 L = Π_m + Π_e + Π_f，其中 Π_m、Π_e、Π_f分别为系统的动能、弹性能和势能（J）。

对拉格朗日系统下的单元展开讨论，单元的能量泛函为

式（2）中，ρ 为单元密度（kg/m³），u_i 为单元节点位移（m），${\dot{u}}_i$为节点速度（m/s），σ_ij 为单元的平均应力张量（N/mm²），V为单元体积（m³）。非保守力包括边界外力及阻尼力（N）。

采用分部积分及一系列转化，最终得到速度和节点位移的函数：

式（3）是CDEM计算模型的基本微分方程，以积分区域作为计算单元，在积分区域中基于指定的速度和位移，定义由节点位移和速度表达的任意一点的值。通过求解方程得到单元节点速度和位移，即可获得整体的速度和位移。该方程与微分方程等价，可从中求得离散元传统的表达式。

最常用的拉格朗日控制方程是质点弹簧系统，而连续-非连续方法（CDEM）是基于拉格朗日方程系统，该方法的控制方程为质点运动方程，节点受力情况应用式（4）所示：

式（4）中：mⁱ 指第 i 个节点的质量，${\overrightarrow{u}}^i、{\dot{\overrightarrow{u}}}^i、{\ddot{\overrightarrow{u}}}^i$ 分别表示节点 i的位移、速度、加速度矢量，c_m表示节点 i 的阻尼，$m^i{\ddot{\overrightarrow{u}}}^i$表示节点 i的惯性力，$c_m{\dot{\overrightarrow{u}}}^i$ 表示节点 i的阻尼力，$\overrightarrow{F}$表示节点外力，${\overrightarrow{F}}_b{}^d$ 表示块体的变形力，${\overrightarrow{F}}_b{}^c$ 表示块体的刚度阻尼力，${\overrightarrow{F}}_j{}^d$表示界面的变形力，${\overrightarrow{F}}_j{}^c$表示界面的阻尼力。

3 充水条件下含软弱夹层顺倾岩质边坡变形破坏机理研究

3.1 参数取值及工况介绍

采用 AutoCAD 软件绘制广东省开平市某露天矿山西帮二维地质模型，形成面域后导入犀牛Rhino 软件进行分组和网格划分，并采用 Rhino TO GDEM 插件将含地层信息的二维地质模型导入 GDEM，该数值模型采用二维平面应力模型，如图3 所示，二维数值模型边坡高度为 420 mm，底部横向宽度为 2050 mm。各层界面均为完全弹性接触，在模型底部约束竖向位移，且不考虑自重的影响。该模型全部采用三角形单元，共划分2827个三角形网格单元，3617 个单元节点。单元模型采用摩尔-库伦模型，接触面采用脆性断裂模型，为了提高数值模拟精度和计算速度，将模型坡面单元细化，提高坡面区域的收敛性和精度。为记录不同条件下的相关运动参数，以便分析边坡失稳破坏演化机理，分别在凝灰岩层坡顶、坡中和坡底各设置一个监测点，各监测点坐标值为：A1（x：862，y：868）；A2（x： 1176，y：677）；A3（x：1417，y：592）（图4）。

数值模型主要岩层为花岗片麻岩、玄武岩及凝灰岩，并含两层软弱夹层，通过前期收集的该露天矿山工程地质资料，并结合其他类似矿山经验对各岩层参数进行综合取值，如表1所示。

针对现场实际情况而言，矿坑充水的作用分为两大类：一是在水的渗透作用下，弱化岩层和软弱夹层的力学参数；二是为坡面骨架提供压力，起到压坡脚的效果。本次采用连续-非连续方法 （CDEM）数值模拟不同充水条件下含软弱夹层顺倾岩质边坡滑动、变形、破坏演化全过程，主要分为 3 种工况：全水-水位线、边坡高度的1/2水位线、边坡高度的1/3水位线、无水（图5）。

3.2 数值模拟结果分析

通过GDEM软件模拟4种不同充水工况下的含软弱夹层顺倾岩质边坡破坏、滑动、至最终失稳过程，最终各工况下 x方向、y方向位移结果如图6、图7所示（其中，x方向向右为正，y方向向上为正）。由图6 可知，在无水工况下，受边坡在自重作用，凝灰岩岩层坡脚和坡顶、玄武岩坡顶存在少量张拉裂缝，裂缝弥散分布，但未持续贯通形成宏观破裂面，凝灰岩层坡脚岩体在持续受压情况下，x 方向位移值最大，x_max为 74.3 mm；在 1/3 水位线工况下，可以明显发现，边坡整体横向位移均有所减小，在充水作用下，凝灰岩岩层坡脚未见明显裂缝，在该层及玄武岩顶部裂缝相对聚集，且凝灰岩岩层坡脚 x 方向位移值低于 2.5 mm，x 方向最大位移值出现在凝灰岩岩层顶部，x_max为 5.39 mm；在 1/2 水位线工况下，凝灰岩岩层 1/2高度以下未产生新生裂纹，但玄武岩层的裂缝不断增加，同样 x 方向最大位移值出现在凝灰岩岩层坡顶，x_max为4.23 mm；在满水-水位线工况下，与前几种工况相对，横向裂缝数量明显减少，主要以竖向裂缝发育为主，凝灰岩岩层坡顶和坡底位移值相差不大，坡中的横向位移最大，x_max 为2.75 mm。由上述分析可知，不同充水工况下（无水-1/3 水位线-1/2 水位线-满水水位线），x 方向最大位移值依次出现在坡底-坡顶-坡顶-坡中，横向裂纹数量及最大横向位移值也依次减少，说明矿区充水后，在水压力作用下，抑制边坡坡面横向裂缝发育，从而横向位移也有所减小。

a—无充水工况；b—1/3水位线；c—1/2水位线；d—全水-水位线

图7为不同工况下y方向的位移云图，由图7可知，在无水工况下，最大竖向位移位于凝灰岩岩层坡顶，y_max为45.91 mm；在1/3水位线工况下，凝灰岩岩层坡顶竖向位移仍然最大，但在地表水入渗影响，坡脚处竖向位移较上一工况有所增加，但总体来看，竖向位移存在明显下降，y_max为 5.73 mm；在 1/ 2水位线工况下，可以明显发现，在凝灰岩岩层坡中以下区域云图颜色较深，与其余位置相比，位移值明显偏高，y_max为 6.16 mm；在满水-水位线工况下，与无水工况相比，变形最大区域基本一致，但位移值明显较小，y_max为 7.33 mm。综上所述，在地表水作用下会减缓边坡坡面裂缝发育，无水与充水条件相比，竖向位移存在明显下降趋势，但是充水作用易导致水进入裂缝及软弱夹层中，导致局部岩体吸水鼓胀，使y方向位移出现回弹效应，竖向位移值存在小幅度增加。

图8～图10 为 A1、A2、A3 监测点 x 方向及 y 方向的时步-位移曲线。由图可知，各监测点在不同充水条件下的时步-位移曲线基本保持一致趋势，且均表现出完整的蠕变三阶段，即：减速蠕变阶段、等速蠕变阶段及加速蠕变阶段。减速蠕变阶段为应变（位移）随着时间延续而增加，但增速有所减缓；等速蠕变阶段为应变（位移）随着时间而增加，增速基本保持趋于一个固定值；加速蠕变阶段为应变（位移）随着时间而增加，增速不断增加。由图8a 可知，各监测点在减速蠕变阶段前期在自重作用下，边坡内部空隙及裂隙压密，导致不同工况下横向位移存在明显负向变形，该阶段下的位移值均小于 0.5 mm，进入等速蠕变阶段后，蠕变速率及横向位移基本保持不变；随着时步不断增加，逐渐进入加速蠕变阶段，此时边坡受自重及充水作用下，裂缝持续扩展及贯通，横向位移明显增加，1/3工况下监测点A1率先进入加速蠕变阶段，其次为1/2工况、满水工况，最后为无水工况，图9a、图10a 中 A2、A3 监测点变形趋势基本一致。由图8b～图10b 可知，不同充水工况下，各监测点在减速蠕变初期的竖向位移量基本一致，集中于 1.5~2 mm 之间；进入等速蠕变阶段后，竖向位移及变形速率无明显增加，各工况下边坡基本处于动态平衡状态，边坡内部裂缝扩展缓慢，但未连接贯通；随着时步不断增加，凝灰岩岩层裂缝持续增多，竖向裂缝相连，形成宏观破裂面，最后进入加速蠕变阶段，竖向位移显著增加，最终边坡发生失稳破坏。其中，A1监测点进入加速蠕变阶段先后顺序为：满水工况-1/3 水位工况-1/2 水位工况-无水工况；A2及A3监测点进入加速蠕变阶段时间基本一致。

a—无充水工况； b—1/3水位线；c—1/2水位线；d—全水-水位线

a—无充水工况； b—1/3水位线；c—1/2水位线；d—全水-水位线

a—x方向位移值；b—y方向位移值

a—x方向位移值；b—y方向位移值

a—x方向位移值；b—y方向位移值

4 结论

本文应用连续-非连续方法（CDEM）动态模拟不同充水条件下含软弱夹层顺倾岩质边坡滑动破坏全过程，通过设置不同监测点，系统分析边坡横向及竖向位移变化规律，主要结论如下：

（1）不同充水工况下（无水-1/3水位线-1/2水位线-满水水位线），x 方向最大位移值依次出现在坡底-坡顶-坡顶-坡中，横向裂纹数量及最大横向位移值也依次减少，说明矿区充水后，在水压力作用下，抑制边坡坡面横向裂缝发育，从而横向位移及竖向位移也有所减小。

（2）无水及充水条件下，含软弱夹层顺倾岩质边坡同样表现出完整的蠕变三阶段理论，即减速蠕变阶段、等速蠕变阶段及加速蠕变阶段。

（3）1/3充水工况易导致边坡变形提前进入加速蠕变阶段，1/2充水工况及满水工况变形曲线基本一致，无水条件下等速蠕变持续时间最长。

（4）矿坑充水后起到压坡脚效果，但需合理设置水位高度，即能抑制边坡横向及竖向变形，同时也需考虑过高水位对边坡岩体的软化效果。

参考文献