基于PFC2D的某矿区反倾岩质边坡稳定性分析

何 倩，薛喜成

(1.徐州中国矿大岩土工程新技术发展有限公司，江苏 徐州 221000；2.西安科技大学地质与环境学院，陕西 西安 710054)

0 引 言

目前，对于边坡稳定性评价的方法主要有定性分析和定量评价[1]，定量评价多采用极限平衡法和数值分析法，极限平衡法通常不能真实地反映边坡失稳时的应力场和位移场[2]，存在一定局限性，而数值分析法越来越受到国内外学者的重视，常见的数值模拟方法有有限差分法、有限单元法和离散单元法。颗粒流法(Particle Flow Code，PFC)是基于离散单元模型框架的一种定量评价方法，于1971年Cundall[3]首次提出，适用于研究非连续介质力学问题，可体现多相介质的不同物理关系，颗粒之间可彼此分离，能够反映岩块或土体之间的滑动、平移和分离等大位移。通过定义颗粒为有质量和表面的刚性体和颗粒之间的接触模型，模拟颗粒的运动以及颗粒之间的相互作用，反映运动机理、过程和结果[4]。贺续文等[5]采用PFC建立含密集节理的岩体边坡模型，研究了节理连通率对边坡破坏形式的影响。Camones等[6]利用PFC研究裂隙岩体的机理以及裂纹扩展在岩质边坡阶梯状滑移中的作用。张小勇等[7]采用PFC建立顺层页岩边坡模型，分析其稳定状态并对其失稳过程进行研究。杨天鸿等[8]研究了露天矿边坡的动态稳定性评价分析预测方法。王培涛等[9]将强度折减法引入离散元法中研究边坡稳定性，将颗粒间粘结强度及摩擦系数进行折减，对边坡定性预测并监测其破坏失稳过程。Zheng等[10]利用UDEC模拟反倾层状岩体边坡的倾倒破坏，研究影响边坡稳定的因素。张志飞等[11]基于PFC模拟反倾层状岩质边坡变形破坏过程，分析变形破坏机制的影响。本文采用离散元二维颗粒流程序(PFC2D)，研究某反倾岩质边坡的稳定状态，并预测其失稳过程，分析破坏过程，为后期边坡加固处理以及矿区开采安全提供参考。

1 工程地质概况

矿区位于燕山山前平原与燕山山脉接壤的地带，采矿方法为台阶式露天开采。由于石灰矿开采的影响，使地层长期缓慢移动，坡体周围分布不同程度的地表裂缝，一些甚至出现了沿倾向的裂缝，使沿走向的裂缝连通起来，说明裂缝对地层稳定的影响在逐渐加剧，对矿区安全开采极为不利。矿区内出露地层主要为奥陶系、石炭系、二叠系岩层及第四系冲洪积、残坡积物，矿区内发育一定规模的断层，矿体赋存于奥陶系上马家沟组中段地层中，边坡上部岩性主要为浅灰～土黄色泥质灰岩，中下部为灰～浅灰褐色厚层状白云质灰岩和豹皮状灰岩。边坡示意见图1。由现场调查分析，该边坡坡顶存在拉张裂缝，边坡产状SSE167°∠60°，坡面存在少量破碎岩体，易发生岩体滚落。

图1 边坡示意

2 边坡模型的建立

2.1 单轴压缩试验标定

PFC中岩体的参数是细观参数，这些参数与实际材料的物理力学参数没有确切的对应关系，不能将试验得到的宏观力学参数直接应用于PFC的模拟之中。目前，在PFC中一般采用试错法进行细观参数的确定，其确立的过程需要根据岩石的弹性模量、抗压强度、泊松比等宏观参数进行不断地调整，使细观参数模拟的结果与岩石的宏观参数相对应，此时得到的颗粒细观参数就是边坡的稳定性分析时所需要的细观参数。

采用PFC内置的单轴压缩模拟程序确定颗粒的细观参数[4，12]，颗粒的半径越小，越接近真实的情况。但受限于计算机的运行速度及内存容量等因素，可适当调整颗粒半径[13]，单轴压缩颗粒半径与后续边坡模型颗粒半径一致。刚度模型采用线性刚度模型，滑动模型采用摩尔-库伦模型，粘结模型采用平行粘结模型。通过岩体的物理力学参数标定，得到PFC模型细观参数(见表1)。

表1 PFC模型细观参数

通过单轴压缩试验得到边坡岩体的应力应变关系见图2。计算得到数值模拟单轴压缩试验下岩体的弹性模量和抗压强度的模拟值，与宏观力学试验实测值对比见表3。从表3可知，模拟值与实测值的误差百分比在2%左右[14]，符合数值模拟要求，可用于后续的边坡稳定性分析之中。

图2 应力-应变关系

表2 数值模拟值与宏观力学实测值对比

2.2 边坡数值模型建立

本文采用删除法建立基于PFC2D软件的边坡模型[4，12]，具体步骤如下：①在AutoCAD中使用多线段命令绘制边坡的几何模型，边坡模型宽330 m，高170 m，保存为DXF格式文件；②在PFC2D中通过Geometry import命令导入DXF文件并通过wall命令生成墙体，使用ball distribute命令在矩形范围内生成一定半径范围的颗粒集合，删除边坡墙体外的多余颗粒；③对填充好的边坡模型中的颗粒进行分组并进行初始化应力场，设置循环为最大不平衡力达最大接触力的1.0×10-5时颗粒集合达到平衡状态，消除悬浮颗粒；④设置边坡模型的顶面为自由面，将底面及左右两侧采用fix命令进行刚性墙约束，生成边坡模型。颗粒之间采用平行粘结接触，粒径范围0.4～0.5 m之间，共生成颗粒49 432个。

在泥质灰岩层中布置4个测量圆，从坡脚至坡顶坐标分别为测量圆1(207，70)、测量圆2(235，100)、测量圆3(265，127.5)和测量圆4(294.2，149)，监测记录这些位置在坡体滑动过程中的应力和孔隙率的变化状态，并对图中A～H颗粒X方向位移进行监测。监测点A、B在边坡前缘，监测点C、D在边坡中前部，监测点E、F在边坡中后部，监测点G、H在边坡后缘。建立的边坡模型见图3。深色颗粒为白云质灰岩层，浅色颗粒为泥质灰岩层。

图3 边坡离散元模型(单位：m)

3 岩质边坡稳定性分析

3.1 边坡运动特征

在模拟时，对边坡施加重力作用，同时对模型中颗粒的位移和速度清零，以便于后续对边坡位移和速度的记录，循环至边坡达到稳定状态并记录边坡破坏过程中的应力、孔隙率以及裂纹数量的变化情况。边坡随时步进行的破坏过程见图4。从图4可知：

图4 边坡破坏过程模拟

(1)计算运行至3×103时步时，边坡在重力作用下整体发生向下的位移变形，边坡滑体中上部岩体的抗剪强度较弱，岩体较破碎，颗粒发生局部剪切破坏，位移明显，坡脚处颗粒累积一定压力后发生位移，并且边坡出现明显的滑动带，颗粒最大位移0.1 m。

(2)运行至3×104时步时，坡脚处受到滑体中部以及坡后缘处不断累积的重力挤压作用产生应力集中，坡体颗粒冲破坡脚的锁固作用，颗粒向前滑出，坡体中后部岩体受到拉张作用产生1条明显的拉张裂纹，出现变形破坏的趋势，颗粒最大位移达2.18 m。

(3)运行至1.5×105时步时，边坡中后部及后缘处拉张裂纹逐渐增加并逐步扩展，岩体间产生大的缝隙并失去支撑，变形加剧，坡脚发生挤压破坏，颗粒滑出边坡坡面，并逐渐在坡脚处堆积，颗粒最大位移达12.63 m。

(4)运行至5×105时步时，边坡滑体部分不断下滑，坡脚处颗粒滑出后逐渐堆积，坡体中部剪切破坏与后缘拉张裂缝相互贯通，变形严重，原边坡后缘处一处平台小边坡的颗粒全部向前滑动，导致裂纹扩展，岩体破碎，颗粒堆积，颗粒最大位移达55.12 m。

(5)运行至1×106时步时，颗粒最大滑动距离达71.79 m，边坡滑体相互贯通，向下滑移，最终坡体整体破坏。

图5 边坡位移时程

A～H监测点颗粒位移曲线见图5。从图5可知，边坡颗粒在重力作用下同时发生位移，上层的颗粒产生的位移最大，而下层颗粒的位移相对较小，是由于坡体前缘岩体滑动，坡中部发生剪切破坏，后缘岩体失去支撑向下滑动，边坡上部颗粒所受应力大于下部颗粒，导致边坡上层颗粒的位移变化显著，而下层颗粒位移较小，即颗粒B位移最大，颗粒A位移最小。边坡前缘的颗粒较后缘的颗粒产生的位移大，是由于边坡下部颗粒粘结破坏而向下滑动，导致上部颗粒失去支撑产生向下的加速度和动能，推动边坡前部颗粒继续向前滑动，因此前缘颗粒的位移较大，后缘颗粒位移相对较小，即颗粒B位移最大，颗粒G位移最小。边坡位移曲线总体变化趋势大致相同，随着时间的推进，位移不断增加，位移曲线的斜率先迅速增加，然后逐渐变缓，表明颗粒的位移速度先增加，后在摩擦力的作用下位移速度逐渐减小，与边坡破坏的规律相符。

边坡速度矢量见图6。从图6可知，运行至5×103时步时，边坡坡体产生明显的滑动带，滑体部分产生一定速度，其中靠近滑动面下部颗粒的速度相对较大，最大速度为1.01 m/s；运行至1×105时步时，由于坡脚处遭受上覆岩体的重力作用和挤压作用，最大速度出现在边坡坡脚处，为2.36 m/s，而滑体与滑床接触位置的颗粒速度相对较小，由于中部岩体较破碎，导致滑体中部速度明显增加，平均速度0.8～1.0 m/s；运行至2×105时步时，边坡滑体部分速度较均匀，平均速度0.5～1.0 m/s，而坡脚处颗粒受到压力挤出，颗粒速度最大，为5.58 m/s，在边坡中后部以及后缘处由于产生了拉张裂缝使岩体失去一定支撑，导致岩体在重力作用下发生倾倒，颗粒产生一定大小的速度，平均速度1.5～2.0 m/s，方向偏向于坡面内；运行4×105时步及至6×105时步时，滑体颗粒最大速度均出现在坡脚，分别为3.02 m/s和2.94 m/s，颗粒速度呈现出逐渐减小的趋势，是因为边坡上部颗粒受到下部的阻碍，将动能传递至坡脚颗粒，坡脚颗粒被上覆颗粒推出，但在摩擦力的作用下动能逐渐降低，速度逐渐减小；继续运行至1×106时步时，颗粒平均速度减小至0.2～0.6 m/s，后缘处有个别颗粒速度相对较大，为1.36 m/s，是因为边坡变形破坏过程中可能存在不稳定岩块崩落等现象，导致速度突变。

图6 边坡速度矢量

3.2 岩质边坡破坏机理

边坡的破坏过程中设立测量圆监测边坡滑动破坏过程中的应力变化，应力变化见图7。从图7可知，边坡底部颗粒由于受到重力作用和挤压作用，颗粒向下滑动并在坡脚堆积。测量圆1中X和Y方向的应力显著增加，其中Y方向的应力主要来源于上覆土压力；测量圆2处于边坡中前部，由于受到上覆岩土层重力作用向下挤压以及测量圆内的岩土层持续向下滑动的作用，使得圆内的应力不断波动，但总体仍呈现下降的趋势；测量圆3处于边坡的中后部，上覆岩土层逐渐滑出测量圆范围，所受的竖向压力和挤压力逐渐减小，因此测量圆内X和Y方向应力逐渐减小，最终当岩土体全部滑出测量圆范围后应力趋近于零；测量圆4处于边坡后缘顶部，边坡顶部颗粒产生拉张裂缝并持续向下滑动，因此测量圆内颗粒应力明显减小，最终趋于零。

图7 测量圆应力变化

不同位置岩体孔隙率随时间变化见图8。边坡前缘测量圆1和边坡中前部测量圆2中颗粒受到重力作用向下滑出，因此孔隙率在前期减小，但上覆岩体在失去下部支撑后迅速跟进，使得颗粒压密，孔隙率反弹；边坡中后部测量圆3处岩体剪切破坏，颗粒逐渐减少，孔隙率增大；坡顶处测量圆4中颗粒向下滑动直至全部滑出，孔隙率增大。

图8 不同位置岩体孔隙率变化

在PFC软件中，若颗粒间实际拉力或切向力大于相应的粘结强度，则颗粒间的粘结就会破坏，而发生破坏的粘结强度就是颗粒间产生的裂纹[4]。边坡破坏过程中产生的裂纹数量随时间变化见图9。从图9可知，裂纹数量随时步先急剧增加后逐步减缓，边坡坡脚首先产生大量细小裂纹，坡脚粘结破坏，边坡中后部及坡顶处颗粒失去下部支撑向下滑动，导致坡体产生几处贯通的拉张裂缝，随后裂纹在边坡滑体中逐步扩展，导致边坡粘结大量破坏并最终失稳。最大裂纹数达3 723个，并且裂纹数量还在缓慢发展。

图9 裂纹数量随时间变化

4 结 语

本文采用离散元二维颗粒流程序(PFC2D)，对某矿区中的反倾岩质边坡的稳定性进行研究，得出以下结论：

(1)对比边坡岩土体标定的细观力学参数与实测宏观力学参数值，PFC建模所标定的细观力学参数符合边坡后续的稳定性分析的要求。

(2)该方法能够在边坡失稳破坏过程中直观地反映边坡颗粒运动过程以及运动过程中所搜索出的滑动面。

(3)边坡中部破碎岩层首先发生剪切破坏，坡脚处颗粒因重力及上覆颗粒挤压的作用，颗粒向外滑出，导致坡顶失去底部岩体的支撑产生拉张裂缝，颗粒向下滑动，最终边坡整体失稳破坏为倾倒失稳破坏。

(4)矿区在开采时应保证一定坡度比例进行开挖，边坡适当设置防护网，拦截崩落碎石，并结合实际情况，设置监测设施，保证矿区安全开采。