厚煤层一般指厚度超过3.5 m的煤层，特厚煤层一般指厚度超过8 m的煤层[1]。特厚煤层矿井数量较多，主要分布于我国西部矿区。由于一次采出煤层厚度较大，上覆岩层变形破坏更为剧烈，在地表形成的采动裂缝更为普遍[2]。随着我国煤炭资源开采逐渐向西部转移，对特厚煤层下开采引起的地表破坏研究显得尤为重要。相关学者采用相似材料模拟、数值模拟、理论分析等方法对特厚煤层下开采岩层运动规律及其控制问题进行了研究。比如戴华阳等[3]以窑街矿区为原型，采用相似材料模拟研究了特厚急倾斜煤层岩层移动机理和地表移动规律，揭示了特厚急倾斜煤层开采具有显著的重复开采和变方向传播的特征；许家林等[4-5]对特大采高工作面开采覆岩结构形态及其运动规律进行了研究，得出不同的覆岩结构形式下，工作面矿压显现规律的差异。另外，针对坡体下煤层开采时斜坡变形的破坏机理[6-9]、裂缝发育规律[10-12]、地表损害及防治[13]等，相关学者也进行了大量的研究工作，取得了丰硕的成果。
现有研究主要是针对特厚煤层开采或者坡体下开采所进行的单一研究，对特厚煤层开采引起的斜坡变形研究相对较少，且覆岩运动及地表裂缝发育是一个非连续大变形离散介质问题，以往的方法不能很好地研究其动态演化规律。20世纪末，Peter Cundall提出了颗粒流理论(又称粒子流理论)，即Particle Flow Code(PFC)[14]。最初，这种方法是研究颗粒介质特性的一种工具，它将物体划分成数百个颗粒单元，利用局部模拟结果来研究计算边值问题。随着计算机技术的不断发展，颗粒流方法逐渐被用于模拟固体力学，解决实际工程问题[15]，大量研究证明颗粒流方法更适用于研究煤层开采引起的沉陷问题。比如周健等[15]归纳总结了颗粒流模拟方法产生的背景，比较了与其他模拟方法的异同之处，并介绍了其在岩土工程领域的应用实例；LIAN et al[16]、WANG et al[17]、侯恩科等[18]采用颗粒流方法研究了煤层开采过程中地表裂缝发育规律，验证了PFC模拟地表非连续变形的可行性。以上研究成果和方法手段对特厚煤层采动诱发斜坡变形破坏研究具有重要的参考价值。
本文以西部某矿8092工作面为研究背景，采用现场调查、数值模拟与理论分析相结合的方法，分析了特厚煤层条件下开采斜坡变形破坏特征和机理，旨在为颗粒流方法在此方面的应用、采动斜坡变形评价及地面保护提供借鉴。
1 工程概况
西部某矿区地势起伏较大，大量冲沟发育，地表有厚黄土层覆盖，为典型的黄土沟壑区。所研究8092工作面走向长400 m，倾向长127 m，开采煤层为8#煤层，平均煤厚8.4 m，采深170～250 m，平均煤层倾角为6°.工作面顶板为粉砂岩，厚11.6 m，底板为粉砂岩，厚30 m，地表黄土层最大厚度为85 m.工作面覆岩组成见表1.工作面采用综采放顶煤，走向长壁后退式采煤法，全部垮落法管理顶板，采掘进度约为2.5 m/d.图1为工作面布置平面图。
表1 工作面覆岩组成
Table 1 Composition of coal seam overburden
序号岩性厚度/m深度/mS10黄土层5～8585S9泥岩26111S8砂质泥岩32143S7中粗粒砂岩28171S6砂质泥岩34205S5中粗粒砂岩20225S4泥岩13238S3粉砂岩11.6249.6S2煤层8.4258S1粉砂岩30288
图1 工作面布置平面图
Fig.1 Working face location
井田内东低西高，地表西侧坡体较陡，平均坡度为24°，坡顶坡底起伏落差最大高约75 m.随着井下开采的不断进行，位于采空区上方的覆岩发生移动变形，开采影响向上传播，导致地表出现大量裂缝。依据调查结果，坡体面发育有多条裂缝，裂缝宽度为0.1～0.2 m，坡顶和坡底出现明显的台阶状裂缝，其水平张开量和垂直错动量较大，严重影响生态环境(见图2).
图2 地表裂缝
Fig.2 Surface cracks
2 颗粒流数值模拟及结果分析
2.1 颗粒流程序原理
颗粒流程序PFC(particle flow code)是一种离散元程序，与基于连续介质力学的方法不同，颗粒流理论将颗粒间的微观力学特性与宏观力学特性联系在一起，通过相互作用反映颗粒的宏观状态，因此，PFC中没有像有限元法那样的网格和网格概念[12]。颗粒流理论通过两个或多个颗粒连接形成任意组合，利用力-位移定律与牛顿第二定律交替进行模拟块体运动问题，图3为其计算循环过程[14]。
图3 颗粒流程序(PFC)的计算循环过程
Fig.3 The calculation cycle process of particle flow program (PFC)
2.2 数值模型设计及参数标定
2.2.1数值模型设计
为了消除数值模型的边界效应，依据实际现场情况，以图1中的A线为地质剖面建立二维颗粒流采场模型，设定模型最大尺寸为864 m×288 m，见图4.模型共计229 193个颗粒，上覆岩层以表1为依据进行分层设置，通过删除煤层所对应的颗粒模拟开采。在模型计算过程中，约束模型左、右边界水平侧向位移，底部施加固支条件，上边界是自由边界，以重力加载到模型上(重力加速度为9.8 m/s2).开切眼距离模型左边界243 m，每次推进20 m用于模拟煤层沿走向开采，每次开挖后数值模型不平衡力平均比小于等于1×10-5则为平衡。图4为采用PFC2D建立的数值模型。
图4 颗粒流数值模型
Fig.4 Particle flow code numerical model
2.2.2参数标定
对于岩石断裂问题，绝大部分学者采用平行粘结模型来进行研究。但是研究发现该模型的颗粒集合体往往具有较小的压拉比，与实际岩石类材料不符，平直节理模型则克服了这一缺陷，更适用于研究岩石力学性质[17，19]。因此本文中的数值计算模型采用平直节理模型定义。颗粒流模型的微观参数标定采用“试错法”，具体步骤为：对岩土体进行数值模拟单轴压缩试验，通过不断调整模型的微观参数，使得室内实验结果与数值模拟结果相匹配，最终确定相应的模型参数。表2为不同岩层岩体标定的微观参数。
表2 岩层微观力学参数
Table 2 Micromechanical parameters of rock strata
符号描述黄土层S10砂质泥岩S6、S8中粗砂岩 S5、S7泥岩S4、S9煤层S2粉砂岩S1、S3γ/(kN·m-3)容重182425241624Rmin/cm颗粒最小半径404040404040Rmax/Rmin颗粒半径比1.61.61.61.61.61.6E*/GPa平节理有效模量0.414.22213.5428.2K*平节理刚度比222222σc/MPa平节理抗拉强度平均值与标准差0.1/0.0250.8/0.21/0.250.8/0.20.25/0.062 51.2/0.3c/MPa平节理粘聚力平均值与标准差4/120/520/520/510/2.520/5
“层面效应”在数值模拟中是不可忽视的，考虑层面影响的计算模型才更符合工程实际。本文考虑了岩层间界面对特厚煤层采动诱发斜坡变形破坏的影响。岩层间界面作为接触面处理，层面厚度为最小颗粒半径。首先通过group命令对不同岩层进行分层，其次在CAD中绘制层面几何模型，通过geometry命令将其导入PFC程序中的数值模型，将层面置于岩层与岩层之间，最后通过cmat命令对不同层面赋予了相应的力学参数，参数选取参考文献[20-21]进行取值，大致为两侧岩层力学参数的1/10.
2.3 斜坡变形破坏规律
工作面开采阶段，斜坡下覆岩受采动影响出现不同的破坏情况，如图5所示。
图5 工作面推进过程中上覆岩层破坏演化
Fig.5 Overlying strata collapse evolution during working face advance
工作面推进至80 m时，直接顶发生初次垮落，垮落高度为14.4 m，此时坡体前缘受沉陷区域拉应力破坏出现拉裂。随着工作面继续推进，直接顶随采随落，采空区面积逐渐增大，当推进至100 m时，直接顶发生周期垮落，周期垮落步距为20 m.
工作面推进至140 m时(坡底)，由于黄土层的载荷作用，上覆岩层岩块不断被压碎，充填采空区，覆岩破坏范围则呈拱形向上发育，此时斜坡面受到前方土体牵引作用出现拉张裂缝R1、R2.工作面推进至180 m时，直接顶完全垮落至煤层底板上，裂缝R1、R2继续发育，移动变形沿坡体向上发育形成裂缝R3.由于坡体上部所受拉应力增大，坡顶后缘受中下部土体牵引作用出现明显的拉张裂缝R4.
工作面推进至220 m时，坡顶受坡体牵引作用出现裂缝R5.垮落带中部覆岩被压实，坡底裂缝则与工作面联通，形成贯通性采动裂缝R6.工作面回采结束后，坡顶裂缝向下延伸，坡体失稳，并出现向下滑动的趋势，地表出现正台阶状破坏。
根据模拟结果可知，工作面开采阶段，斜坡底部率先受到采动影响，之后移动变形沿坡体由下往上发展。工作面回采结束后在地表共产生6条明显裂缝，裂缝主要存在于斜破底部、斜坡面和斜坡顶，裂缝为正向裂缝，坡顶裂缝(R5)最大宽度为1.17 m，坡底裂缝与工作面联通，坡体面产生3处明显裂缝。
2.4 采场支承压力变化规律分析
工作面推进过程中上覆岩层逐渐垮落至采空区引起采场应力分布，采场应力的变化是覆岩结构动态演化和引起地表采动破坏的直接体现[17]。图6为颗粒流模型采场支承压力与煤层不同开采距离的关系，图中测量圆布置于煤层底板上，自模型左边界向右布置，半径为4.2 m，共102个。
工作面回采阶段，拉应力升高区位于采空区正上方，支承压力的升高区位于采空区前方和后方，采空区前方和后方的最大支承压力值逐渐增大，且采空区前方最大支承压力位置不断向深部转移。工作面回采80 m时，直接顶发生初次垮落，所有监测点最大支承压力值为9.40 MPa，此时坡体前缘受拉应力作用出现拉张破坏。工作面回采100 m时，直接顶发生周期性垮落，所有监测点最大支承压力值为10.39 MPa，采空区前方支承压力峰值增幅较小，这一特征与直接顶传递覆岩载荷能力直接相关，此时坡体未出现明显变形破坏，采空区上方支承应力变化相对于煤体中的应力变化较为复杂。工作面回采180 m时(坡体正下方)，监测点支承压力峰值达到最大值11.22 MPa，这是因为直接顶完全垮落在煤层底板上，直接顶承载结构的失稳致使其失去了传递载荷的能力，且由于坡体面拉应力增幅较大，坡体裂缝继续发育，坡体后缘受牵引作用出现新的裂缝。工作面回采220 m时，直接顶上方覆岩出现整体切落，覆岩裂缝与坡底裂缝联通，此时支承压力峰值急剧下降。工作面回采结束后，坡体后缘出现台阶状破坏，煤体中的支承压力出现回升，所有监测点最大支承压力值为8.57 MPa.图6表明覆岩稳定性对采空区上方的支承压力变化有显著影响。
图6 工作面不同推进距离支承压力分布
Fig.6 Distribution of abutment pressure with different
propulsion distances of working face
2.5 斜坡位移分布规律
煤层开挖会引起上覆岩层的位移、变形和破坏。图7分别为直接顶初次垮落时和工作面开采结束后的覆岩垂直位移场与水平位移场分布图。
图7 开采引起的垂直位移场与水平位移场
Fig.7 Horizontal and vertical displacement caused by mining
工作面推进至80 m时，直接顶发生初次垮落，采空区周围岩层受到水平应力与自重作用向采空区方向移动，采空区上方出现离层现象，由于开挖空间较小，地表黄土层受采动影响较小，覆岩下沉量自下往上逐渐递减。随着工作面推进，覆岩不断填充采空区，坡体则受滑移影响偏向临空面，呈现出采空区正上方覆岩与坡体水平移动方向相反的趋势，裂缝带发育至地表形成采动裂缝。
根据地表不同位置的监测点位移数据，绘制位移曲线。图8为工作面回采过程中的地表下沉曲线。工作面开采初期，下沉曲线较为平缓，坡底先受到采动影响，下沉量较大，坡体受采动影响较小。随着工作面开采范围的扩大，地表的下沉量与影响范围逐渐增大，受坡体滑移影响，最大下沉点偏向上坡方向，下沉曲线呈现出不以采空区为中心的“V”字型，整体表现出盆地边缘偏小，盆地中央偏大的趋势。工作面回采结束后，地表最大下沉量为9 020 mm.
图8 地表下沉曲线
Fig.8 Surface subsidence curve
图9为工作面回采过程中的地表水平移动曲线。工作面开采初期，水平移动曲线较为平缓，出现两个方向相反的极值。随着工作面不断推进，地表水平移动量与地表影响范围逐渐增大，下坡方向的极值向盆地边缘移动。当工作面推进到坡底时，坡体滑移影响加剧，下坡方向的地表水平移动值和影响范围较大，上坡方向的地表水平移动值和影响范围较小，且由于坡体面出现裂缝，导致此时地表对应的监测点水平移动值出现突变。工作面回采结束后，地表最大水平移动量为4 170 mm.从图9可以明显看出，工作面回采阶段，地表水平移动曲线关于零值点呈现出非反对称状态，下坡方向下沉值为零的点，水平移动值并不为零，甚至有很大的水平移动值，表明坡体下开采，地表受采动和地形的影响。地表的移动向量是指向采空区的移动分量和沿坡面指向下坡方向的矢量合成，地表水平移动存在着指向下坡方向的水平滑移量，这与作者在文献[22]中采用相似模拟实验所得的结论是一致的，表明PFC数值模拟软件对地下采煤引起的覆岩破坏及地表移动变形研究是可信且正确的。
图9 地表水平移动曲线
Fig.9 Surface horizontal displacement curve
根据地表岩移资料，8092工作面开采结束后，实测地表最大下沉值为8 920 mm，最大水平移动量为3 800 mm.8#煤开采结束后数值模拟结果地表最大下沉值为9 020 mm，最大水平移动量为4 170 mm，通过对比可以得出数值模拟结果与实际相差不大。
3 斜坡变形破坏机理
通过对工作面回采阶段采空区上方覆岩的变形破坏规律、支承压力变化规律及位移分布规律进行分析，得出特厚煤层条件下开采斜坡变形具有阶段性，可以将斜坡的变形过程分为以下几个阶段：
1) 坡底拉裂破坏：工作面开采初期，采空区上覆岩层由于自重力下沉，地表黄土层随之发生弯曲变形，同时坡体向采空区内侧移动。斜坡底部受采空区岩土体的拉应力率先产生裂缝，如图10(a).
2) 坡体牵引破坏：工作面推进至坡体下方时，地表最大沉陷值位于采空区中心，坡底沉陷程度加深，受拉应力增大作用裂缝继续延伸。采动影响由坡底向坡体延伸，坡体受前方岩土体牵引作用产生拉张裂缝，破坏模式为牵引式破坏，如图10(b).
图10 斜坡变形破坏过程
Fig.10 Slope deformation failure process
3) 坡顶滑动破坏：随着工作面继续推进，采空区范围不断增大，上方坡体受下方坡体牵引作用向临空面移动，坡顶产生拉张裂缝，同时坡底土体受采空区作用向斜坡内移动，坡底受到双向的挤压，裂缝破坏加剧。工作面回采结束后，坡底裂缝与采空区联通形成贯通性裂缝，坡体面与坡顶裂缝不断延伸，坡体失稳，并出现向下滑动的趋势，如图10(c).
4 结论
本文以西部某矿区8092工作面为背景，采用颗粒流软件PFC建立了复杂斜坡条件下的数值计算模型，并分析了工作面推进过程中的斜坡变形破坏特征及机理。主要结论如下：
1) 煤层开采阶段，采空区上方覆岩破坏范围呈拱形向上发育。斜坡底部率先受到采动影响产生裂缝，之后移动变形沿坡体由下往上发展，在坡体面与坡顶产生裂缝。煤层开采结束后，坡底裂缝与采空区联通形成贯通性裂缝。
2) 直接顶周期垮落前，支承压力分布较为规律，随着工作面不断推进，采空区上方支承应力变化相对于煤体中的应力变化较为复杂，采空区前方最大支承压力位置不断向深部转移，且工作面推进至坡体正下方时，监测点支承压力峰值达到最大，之后支承压力峰值会急剧下降。采场应力的变化反映了覆岩结构的动态演化过程。
3) 斜坡下开采，地表下沉曲线呈现出不以采空区为中心的“V”字型，地表水平移动曲线关于零值的点呈现出非反对称状态。由于地表水平移动存在着指向下坡方向的水平滑移量，下坡方向下沉值为零的点，水平移动值并不为零，甚至有很大的水平移动值。
4) 特厚煤层条件下开采斜坡变形具有阶段性，斜坡区经历了坡底拉裂破坏—坡体牵引破坏—坡顶滑动破坏的过程。颗粒流理论最大优势在于可以直观看出覆岩结构与裂隙发育的动态演化过程，更适用于研究采动覆岩非连续变形，对于采动裂缝的预测研究具有重要意义。
[1] 杨永康.特厚煤层大采高综放采场覆岩移动规律及围岩控制研究[M].北京:煤炭工业出版社,2014.
[2] 高超，徐乃忠，刘贵.特厚煤层综放开采地表沉陷预计模型算法改进[J].煤炭学报,2018,43(4):939-944.
GAO C,XU N Z,LIU G.Improvement on prediction model algorithm for surface subsidence of extra thick seam using fully-mechanized top coal caving method[J].Journal of China Coal Society,2018,43(4):939-944.
[3] 戴华阳，郭俊廷，易四海，等.特厚急倾斜煤层水平分层开采岩层及地表移动机理[J].煤炭学报,2013,38(7):1109-1115.
DAI H Y,GUO J T,YI S H,et al.The mechanism of strata and surface movements induced by extra-thick steeply inclined coal seam applied horizontal slice mining[J].Journal of China Coal Society,2013,38(7):1109-1115.
[4] 许家林，鞠金峰.特大采高综采面关键层结构形态及其对矿压显现的影响[J].岩石力学与工程学报,2011,30(8):1547-1556.
XU J L,JU J F.Structural morphology of key stratum and its influence on strata behaviors in fully-mechanized face with super-large mining height[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(8):1547-1556.
[5] 许家林，朱卫兵，鞠金峰.浅埋煤层开采压架类型[J].煤炭学报,2014,39(8):1625-1634.
XU J L,ZHU W B,JU J F.Supports crushing types in the longwall mining of shallow seams[J].Journal of China Coal Society,2014,39(8):1625-1634.
[6] 汤伏全，田家琦，于根.地下开采诱发滑坡的机制分析[J].辽宁工程技术大学学报：自然科学版,1991,13(S1):16-17.
TANG F Q,TIAN J Q,YU G.Mechanism analysis of landslide induced by underground mining[J].Journal of Liaoning Technical University (Natural Science),1991,13(S1):16-17.
[7] 王创业，李俊鹏，张琪，等.地下采动条件下露天矿边坡覆岩移动规律研究[J].矿业研究与开发,2018,38(4):75-78.
WANG C Y,LI J P,ZHANG Q,et al.Study on the movement law of overlying strata in an open-pit slope during underground mining[J].Mining Research and Development,2018,38(4):75-78.
[8] 邓洋洋，陈从新，夏开宗，等.地下采矿引起的程潮铁矿东区地表变形规律研究[J].岩土力学,2018,39(9):3385-3394.
DENG Y Y,CHEN C X,XIA K Z,et al.Ground movement and deformation caused by underground mining in eastern area of Chengchao iron mine[J].Rock and Soil Mechanics,2018,39(9):3385-3394.
[9] 李腾飞，李晓，李守定，等.地下采掘诱发斜坡失稳破坏机制研究--以武隆鸡尾山崩滑为例[J].岩石力学与工程学报,2012,31(Z2):3803-3810.
LI T F,LI X,LI S D,et al.Mechanism of slope failure induced by undermining-acase study of jiweishan rockslide in wulong[J].Chinese Journal of Rock Mechanics and Engineering,2012,31(Z2):3803-3810.
[10] 史文兵，黄润秋，赵建军，等.山区平缓采动斜坡裂缝成因机制研究[J].工程地质学报,2016,24(5):768-774.
SHI W B,HUANG R Q,ZHAO J J,et al.Genetic mechanism of mining cracks of gently inclined slope in mountains[J].Journal of Engineering Geology,2016,24(5):768-774.
[11] 单晓云，姜耀东，王乐杰，等.地下采煤对巍山山体裂缝影响的有限元分析[J].煤炭学报,2008,33(1):23-27.
SHAN X Y,JIANG Y D,WANG L J,et al.Finite element analysis of effect of underground coal mining on landslide of the Wei Mountain[J].Journal of China Coal Society,2008,33(1):23-27.
[12] 崔杰，王兰生，王卫，等.采空区边坡变形破裂演化机制研究[J].采矿与安全工程学报,2008,25(4):409-414.
CUI J,WANG L S,WANG W,et al.Deformation and fracturing mechanism of goaf slope[J].Journal of Mining & Safety Engineering.2008,25(4):409-414.
[13] 赵兵朝，孙浩，郭亚欣，等.湿陷性黄土覆盖区煤层开采地表裂缝发育规律研究[J].矿业研究与开发,2021,41(6):106-110.
ZHAO B C,SUN H,GUO Y X.Study on the cracks development law of collapsible loess surface in coal mining[J].Mining Research and Development.2021,41(6):106-110.
[14] ZHANG Y J,YAN Y G,DAI H Y,et al.Stability and Force Chain Characteristics of “nclined Step Cutting Body” in Stope[J].Applied Sciences,2022,11:10276.
[15] 周健，池永，池毓蔚，等.颗粒流方法及PFC2D程序[J].岩土力学,2000,21(3):271-274.
ZHOU J,CHI Y,CHI Y W,et al.The method of particle flow and PFC2DCode[J].Rock and Soil Mechanics,2000,21(3):271-274.
[16] LIAN X G,ZHANG Y J,YUAN H Y,et al.Law of movement of discontinuous deformation of strata and ground with a thick loess layer and thin bedrock in long wall mining[J].Applied Sciences,2020,10(8):2874.
[17] WANG C L,ZHANG C S,ZHAO X D,et al.Dynamic structural evolution of overlying strata during shallow coal seam longwall mining[J].International Journal of Rock Mechanics and Mining Sciences,2018,103:20-32.
[18] 侯恩科，从通，谢晓深，等.基于颗粒流的浅埋双煤层斜交开采地表裂缝发育特征[J].采矿与岩层控制工程学报,2020,2(1):16-24.
HOU E K,CONG T,XIE X S,et al.Ground surface fracture development characteristics of shallow double coal seam staggered mining based on particle flow[J].Journal of Mining and Strata Control Engineering,2020,2(1):16-24.
[19] 陈鹏宇.PFC2D模拟裂隙岩石裂纹扩展特征的研究现状[J].工程地质学报,2018,124(2):253-264.
CHEN P Y.Research progress on PFC2Dsimulation of crack propagation characteristics of cracked rock[J].Journal of Engineering Geology,2018,124(2):253-264.
[20] 夏磊，曾亚武，张森.层理面细观参数对层状岩体强度特性影响的研究[J].长江科学院院报,2016,33(7):68-75.
XIA L,ZENG Y W,ZHANG S.Influence of meso-mechanical parameters of bedding plane on strength characteristics of layered rock mass[J].Journal of Yangtze River Scientific Research Institute,2016,33(7):68-75.
[21] 王培涛，杨天鸿，于庆磊，等.节理边坡岩体参数获取与PFC2D应用研究[J].采矿与安全工程学报,2013,30(4):560-565.
WANG P T,YANG T H,YU Q L,et al.On obtaining jointed rock slope geo-parameters and the application of PFC2D[J].Journal of Mining & Safety Engineering,2013,30(4):560-565.
[22] 戴华阳，翟厥成，胡友健.山区地表移动的相似模拟实验研究[J].岩石力学与工程学报,2000,19(4):501-504.
DAI H Y,ZHAI J C,HU Y J.Testing study on surface displacement of mountainous region with similar material[J].Chinese Journal of Rock Mechanics and Engineering,2000,19(4):501-504.