0 引 言
近年来，我国东部煤炭资源储量逐渐匮乏，为了满足对煤炭开采产量的需求，煤炭开采的重心开始从东部向西部转移。构造简单、开采条件优越的厚煤层在西部煤田占重要比例，是我国能源发展中重要的战略资源[1-3]。但是在西部地区煤层开采过程中，经常遇到与东部地区不一样的地质条件，其中最具代表性的是具有胶结性差、遇水泥化、强度低且易崩解等特点的弱胶结地层。弱胶结地层多以砂岩为主，其胶结性质、水理性质等力学性质与常规砂岩或软岩还具有显著的差别，大部分处在含水层中[4-6]。覆岩内部裂隙受工作面开采扰动的影响，逐渐发育形成导水、导砂通道，其次，由于弱胶结地层多处于含水层中，其结构遇水易溶解且强度较低，容易形成水砂混合物涌入工作面，造成严重的经济损失[7]，同时水砂混合物的涌入使得含水层水压降低，造成砂层厚度减少，破坏了上覆岩层的整体平衡，对矿井工作面的安全生产造成很大威胁。
煤矿发生溃水溃砂灾害与上覆含水层的规模与性质、开采方式、开采厚度、覆岩破坏形式等因素相关，很多专家学者从不同角度分别研究了溃水溃砂事故的发生机理以及防治手段。任胜文[8]利用水砂源、动力源、通道、流动空间4个主控因素的相互作用得到了溃水溃砂的内在机理。隋旺华等[9-10]通过建立采掘溃砂安全评判模型，利用溃砂通道宽度及含水层水压从深层次进一步得到溃水溃砂机理。文献[11-12]分别利用自主研制的溃水溃砂气液联动相似模拟试验装置以及采动覆岩溃水溃砂灾害模拟试验系统，从不同的角度分析了溃水溃砂通道的发育过程及溃水溃砂的区域特征，得到形成溃水溃砂通道的整个过程。贾后省等[13]利用相似模拟、理论分析、现场实测等综合研究方法，从纵向贯通裂隙张开闭合影响机制的角度进行解释，得到纵向贯通裂隙张开闭合规律的4个阶段，对工作面控制溃水溃砂的隐患具有理论意义。
综上所述，大量学者通过相似模拟试验及理论分析研究了溃水溃砂灾害的过程，并得到了相关控制因素及水砂混合流动的规律，然而煤层上覆含水层的水压、煤层开采厚度等多种因素都会造成煤矿溃水溃砂灾害的发生。由于西部某矿开采的21-1煤上覆岩层内部有古近系砂砾岩含水层，该砂砾岩含水层主要以砂砾岩为主，为裂隙-孔隙承压含水层，补给条件差、富水性弱，是开采21-1煤时的充水水源。因此，笔者利用离散元数值模拟软件PFC3D，模拟不同采厚情况下工作面开采时上覆岩层的垮落情况，从溃水溃砂通道发育过程及含水层水压变化2个角度，研究分析弱胶结地层在不同采厚情况下发生溃水溃砂灾害的可能性，为该矿保持砂砾岩含水层的完整性及防治溃水溃砂事故的发生提供理论依据。
1 颗粒流数值模拟原理
20世纪80年代，Cundall等[14]提出了颗粒流离散元方法(Discrete Element Method)，以此作为基础结合牛顿第二定律，Itasca公司开发了可模拟非连续性颗粒物质结构及运动定律的DEM软件PFC3D。离散元PFC模拟软件的计算原理是运用时步迭代的方式，使得颗粒与颗粒接触重复使用力-位移定律、牛顿第二定律，具体计算关系如图1所示。PFC颗粒流离散元软件具有多种接触模型，如刚度模型、滑移模型、接触黏结模型和平行黏结模型[15]。由于平行黏结模型能够满足对岩体破裂模拟的要求，因此更适用于模拟密实材料，如岩体类材料，如图2所示。
图1 PFC计算循环示意
Fig.1 Schematic diagram of PFC calculation cycle
kn、ks—颗粒的法向、切向接触刚度；Rmax、Rmin—最大、最小颗粒的半径；Fn、Fs—颗粒之间的法向、剪切接触力；
平行黏结颗粒的法向、剪切刚度
图2 平行黏结模型原理[15]
Fig.2 Schematic diagram of parallel bonding model
在模拟水砂混合流动时，砂颗粒与水之间不仅存在相对速度，还受到水黏滞性的作用，水砂两相之间将产生作用力。当水相速度大于砂粒速度时，砂颗粒受到水的拖曳力，水也受到砂颗粒所带的阻力作用。因此，假设砂颗粒全部为球形颗粒、不可压缩性流体，且不考虑温度等外在因素，则水对于砂颗粒的拖曳力的大小Fd[16-17]为
(1)
式中：Cd0为与颗粒雷诺数相关的拖曳力系数；rp为颗粒半径；ρf为流体密度；vf为流体速度；vp为颗粒速度。
此外，水压梯度力和浮力的作用将影响颗粒流动，设压力梯度为
p，则水压梯度力Fp和浮力Fa分别为[16]
p
(2)
(3)
其中：g为重力加速度。颗粒在水流体中除了受到拖曳力、水压梯度力和浮力外，还受到重力、做变速直线运动产生的Basset力和旋转产生的Magnus力等，因此，忽略其他作用力，只考虑砂颗粒在水的黏滞性作用下产生的力及流体对颗粒的拖曳力与浮力，其中在PFC软件中，假定每一个颗粒的中心都被平均施加同一个流体单元内的拖曳力、浮力等。
2 溃水溃砂数值模拟
2.1 宏-细观参数的选取
根据矿井提供的宏观参数，依据宏-细观参数转换公式[18-19]，经过一系列复杂的计算转换为细观参数。
1)弹性模量经验公式：
E/Ec=a+bln(kn/ks)
(4)
式中：E为弹性模量，GPa；Ec为杨氏模量，GPa；kn/ks为颗粒之间法向与切向的接触刚度比，即Krat；a、b 为常数，取值分别为1.652、-0.395。
2)泊松比经验公式：
ν=cln(kn/ks)+d
(5)
式中：ν为泊松比；c、d为常数，取值分别为0.209、0.111。
3)单轴抗压强度回归性分析：
(6)
式中：σc为抗压强度，
为平行连接法向连接强度，
为平行连接切向连接强度，MPa；a=-0.965，b=2.292，c=1.327。
4)抗拉强度回归性分析：
(7)
式中：σt为抗拉强度，MPa；d=-0.174，e=0.463，f=0.289。
最终计算得到的岩层宏-细观参数见表1。
表1 岩层宏-细观物理力学性质参数
Table 1 Macro-meso physical and mechanical parameters of rock strata
岩层岩性厚度/m宏观参数泊松比弹性模量/MPa抗拉强度/MPa黏聚力/MPa内摩擦角/(°)细观参数KratEmod/GPakn/GPaks/GPak′n/GPak′s/GPaJ10U砂质泥岩3.00.157.41.427.6036.41.204.69.27.65.74.7J10M泥岩3.20.157.21.366.3239.41.204.59.07.45.64.6J10D细粒砂岩3.00.1417.43.065.7738.81.1410.821.618.813.511.8J9泥岩2.00.157.21.366.3239.41.204.59.07.45.64.6J8中粒砂岩6.00.1815.42.584.3534.31.3910.120.214.512.69.0J7砂质泥岩3.50.157.41.427.6036.41.204.69.27.65.74.7J6中粒砂岩6.50.1815.42.584.3535.31.3910.120.214.512.69.0J5砂质泥岩4.00.157.41.607.6035.41.204.69.27.65.74.7J4中粒砂岩7.30.1815.42.584.3535.31.3910.120.214.512.69.0J3细粒砂岩6.70.1417.43.065.7738.81.1410.1821.618.813.511.8J2煤层8.80.276.91.332.4242.62.135.110.24.76.32.9J1中粒砂岩12.00.1815.42.584.3535.31.3910.120.214.512.69.1
注：Emod为接触模量；k′n为平行黏结法向刚度；k′s为平行黏结切向刚度。
2.2 溃水溃砂模型的建立
根据现场所给地质柱状图，构建三维模型，如图3所示。模型长120 m，宽120 m，高65 m，共计10层，由于煤层平均角度为6°，则视为近水平煤层，将不再设置煤层倾角。模型上部边界增加水压1.25 MPa，由于模型只允许颗粒沿着垂直方向进行移动，因此，底部边界设置为固定边界，限制颗粒的垂直方向上的位移，前、后、左、右4个边界也是固定不可移动的。为了更好地展示水砂混合物流动特征，模型第10层设置为弱胶结砂层，选择线性黏结模型进行设置，其中在第10层设置下(10D)、中(10M)、上(10U)3层不同粒径的砂层，10D与10U两层颗粒最大粒径为0.8 m，最小粒径为0.5 m，10M层颗粒最大粒径为0.5 m，最小粒径为0.3 m，模型其余部分颗粒最大粒径为1.66 m，最小粒径为1 m，各部分颗粒粒径比例尺寸较为合理[20]。
图3 溃水溃砂模拟模型
Fig.3 Water and sand inrush simulation model
根据PFC3D模拟软件特点，在模型上方含水层布置测量圆[21-23]，监测上方边界水压在不同煤层采厚情况下的变化情况，测量圆布置如图4所示。模型测量圆粒径均为8 m，测量圆布置长为120 m、宽120 m、高8 m。测量圆只布置在含水层内部。测量圆将始终处于稳定状态，前、后、顶、底、左、右共6侧边界固定不可移动。当颗粒沿垂直方向进行移动时，颗粒将会穿过每层测量圆，因此，测量圆将会根据颗粒的移动来监测含水层水压的变化规律。选取测量圆布置区域的中间一行测量圆监测数据进行分析，便于观察水压变化。
图4 水压测量圆布置
Fig.4 Hydraulic pressure measurement circle layout
2.3 溃水溃砂模拟方案
根据不同煤层采厚情况下的含水层水压变化及溃水溃砂可能性，设置模拟方案见表2。
表2 溃水溃砂模拟设置方案
Table 2 Setup scheme table of water and sand inrush simulation
方案煤层采厚/m开采步距/m方案一2方案二4方案三6方案四815、16、14、16、18、10、10、9、12
3 溃水溃砂模拟结果分析
3.1 溃水溃砂灾害发生可能性分析
煤层开采厚度为2、4 m时，上覆岩层垮落水砂混合物流动如图5所示。上覆岩层随着开采步距的增加逐渐发生垮落，但由于煤层开采厚度为2、4 m，分别约占煤层厚度的1/4与1/2，因此，垮落高度不会延伸至砂层，其内部裂隙向上延伸至基本顶。由于上部边界水压的作用，砂层底板内部裂隙逐渐发育，使得水砂混合物沿裂隙通道下渗，砂层出现下沉现象。模拟开采结束时，工作面开采至120 m，整体模型趋于稳定，此时砂层底板产生的裂隙尺寸较小导致水砂混合物出现堆积现象，砂层不再继续下沉。因此，煤层开采厚度为2、4 m时，工作面正常开采不受任何影响，出现溃水溃砂现象的概率微小。
图5 采厚2、4 m时水砂混合物流动示意
Fig.5 Flow diagram of water-sand mixture when mining height is 2 m and 4 m
煤层开采厚度为6、8 m时，上覆岩层垮落、水砂混合物流动如图6所示。由图6a—图6d可得，此时开采煤层厚度约为煤层厚度的2/3，随着工作面不断向前推进，煤层上覆岩层内部裂隙逐渐向上发育，砂层底板裂隙发育至砂层基本底附近，此时少许裂隙向下发育与上覆岩层内部裂隙出现贯通现象，但因为形成的贯通裂隙尺寸较小，较大砂颗粒容易堵塞贯通裂隙通道，使得水砂混合物出现堆积现象，较小砂颗粒则在水流的作用下通过还未堵塞的裂隙通道流入至工作面内。因此，当煤层开采厚度为6 m时，工作面有可能发生溃水溃砂的现象。
图6 采厚6、8 m时水砂混合物流动示意
Fig.6 Flow diagram of water-sand mixture when mining height is 6 m and 8 m
如图6e—图6h所示，煤层开采厚度为8 m时，此时开采厚度约等于整个煤层厚度，上覆岩层逐渐垮落，竖向裂隙逐渐增多，砂层底板在砂颗粒的作用下，内部裂隙发育，贯通裂隙逐渐增多，但由于贯通裂隙的尺寸较小，容易被较大砂颗粒堵塞，水砂混合物则通过尺寸较大的贯通裂隙流入工作面，对工作面的安全生产造成巨大威胁。因此，当煤层开采厚度为8 m时，工作面极易出现溃水溃砂的现象。
3.2 含水层水压分析
通过对砂层位置布置测量圆监测含水层水压变化，进而分析工作面出现溃水溃砂灾害的可能性。此次测量圆只监测模型z轴方向上的压力，其中正负号只表示方向，向下为负，选择工作面推进距离为15、45、99、120 m时的监测结果(图7)。
由图7a、图7b可得，当煤层采厚为2、4 m时，含水层水压随着工作面的不断推进逐渐降低，主要是由砂层底板裂隙逐渐发育使得水砂混合物下渗造成，但此时砂层底板裂隙尺寸较小，水砂混合物极易堵塞裂隙通道，使得其不再继续下沉，含水层水压不再继续降低，最终含水层水压整体趋势为随着工作面推进距离的不断增加而不断减小，此时煤层开采厚度较薄，不会减小为0，因此，工作面发生溃水溃砂的概率微小。
当煤层开采厚度为6 m，煤层开采厚度约为煤层厚度的2/3时，含水层水压变化如图7c所示。模拟结束时，工作面推进至120 m，覆岩内部与砂层底板贯通裂隙的形成，使得水砂混合物流入工作面内，此时含水层水压持续降低至0。由图7d可得，当煤层开采厚度为8 m，煤层开采厚度约为煤层厚度，贯通裂隙逐渐增多，尺寸逐渐增大，含水层水压进一步降低，但模型长度方向28 m至92 m处保持为0。因此，煤层开采厚度为6、8 m时，工作面将有较大的概率发生溃水溃砂等灾害，其中开采厚度为8 m时，概率最大。
图7 含水层水压变化示意
Fig.7 Schematic diagram of aquifer water pressure change
4 结 论
1)随着煤层开采厚度的逐渐增加，覆岩垮落的高度也不断增加，内部裂隙向上逐渐发育，当煤层开采厚度为2、4 m时，工作面发生溃水溃砂灾害概率微小。
2)煤层厚度为6、8 m时，工作面有较大概率发生溃水溃砂等灾害，因此，当煤层开采厚度为4 m时，既可保证工作面的高效开采，也可保证工作面的安全开采。
3)从含水层水压变化的角度进行分析，随着工作面的不断推进，含水层水压整体变化趋势是降低的，而煤层开采厚度为6、8 m时，模拟结束时，含水层水压均降低为0，表明此时含水层水砂混合物下渗，工作面出现溃水溃砂灾害的概率较大。
[1] 张洪彬, 田成林, 孙 赑, 等.浅埋煤层弱胶结顶板破断规律数值模拟研究[J].山东科技大学学报：自然科学版, 2015, 34(2):36-40.
ZHANG Hongbin, TIAN Chenglin, SUN Bi,et al.Numerical simulation of roof break laws in weakly cemented shallow seam[J].Journal of Shandong University of Science and Technology：Natural Science, 2015, 34(2):36-40.
[2] HUANG Xuchao, FAN Cheng, WANG Rui, et al.Study on gas emission prediction of steeply inclined coal seam applied horizontal slice mining[J].Journal of Engineering Science and Technology Review, 2019, 12(2):17-26.
[3] 孙利辉, 纪洪广, 杨本生.西部典型矿区弱胶结地层岩石的物理力学性能特征[J].煤炭学报, 2019, 44(3):866-874.
SUN Lihui, JI Hongguang, YANG Bensheng.Physical and mechanical characteristic of rocks with weakly cemented strata in Western representative mining area[J].Journal of China Coal Society, 2019, 44(3):866-874.
[4] 孙利辉, 纪洪广, 蒋 华, 等.弱胶结地层条件下垮落带岩层破碎冒落特征与压实变形规律试验研究[J].煤炭学报, 2017, 42(10):2565-2572.
SUN Lihui, JI Hongguang, JIANG Hua, et al.Experimental study on characteristics of broken caving and regularity of compaction deformation of rocks in caving zone in the weakly cemented strata[J].Journal of China Coal Society, 2017, 42(10):2565-2572.
[5] 任建喜, 景 帅, 刘朝科.弱胶结砂岩破坏规律及其声发射特性试验研究[J].煤炭技术, 2019, 38(6):1-4.
REN Jianxi, JING Shuai, LIU Chaoke.Experimental study on deformation and failure laws and acoustic emission characteristics of weak cemented sandstone[J].Coal Technology, 2019, 38(6):1-4.
[6] 纪洪广, 蒋 华, 宋朝阳, 等.弱胶结砂岩遇水软化过程细观结构演化及断口形貌分析[J].煤炭学报, 2018, 43(4):993-999.
JI Hongguang, JIANG Hua, SONG Chaoyang,et al.Analysis on the microstructure evolution and fracture morphology during the softening process of weakly cemented sandstone[J].Journal of China Coal Society, 2018, 43(4):993-999.
[7] 隋旺华，董青红，蔡光桃，等.采掘溃砂机理与预防[M].北京:地质工程出版社，2008.
[8] 任胜文.大采深煤层弱胶结厚层砾岩溃水溃砂灾害研究[J].煤炭科学技术, 2019, 47(9):249-255.
REN Shengwen.Study on disaster of water and sand inrush of weakly cemented thick conglomerate on deep mining coal seam[J].Coal Science and Technology, 2019, 47(9):249-255.
[9] 隋旺华, 梁艳坤, 张改玲, 等.采掘中溃水溃砂机理研究现状及展望[J].煤炭科学技术, 2011, 39(11):5-9.
SUI Wanghua, LIANG Yankun, ZHANG Gailing, et al.Study status and outlook of risk evaluation on water inrush and sand inrush mechanism of excavation and mining[J].Coal Science and Technology, 2011, 39(11):5-9.
[10] 隋旺华, 刘佳维, 高炳伦, 等.采掘诱发高势能溃砂灾变机理与防控研究与展望[J].煤炭学报, 2019, 44(8):2419-2426.
SUI Wanghua, LIU Jiawei, GAO Binglun,et al.A review on disaster mechanism of quicksand with a high potential energy due to mining and its prevention and control[J].Journal of China Coal Society, 2019, 44(8):2419-2426.
[11] 赵启峰, 张 农, 韩昌良, 等.浅埋薄基岩含水层下煤层开采溃水溃砂相似模拟实验研究[J].采矿与安全工程学报, 2017, 34(3):444-451.
ZHAO Qifeng, ZHANG Nong, HAN Changliang, et al.Simulation experiment of water-sand inrush during the mining of the shallow coal seam under roof aquifer with thin bedrock[J].Journal of Mining and Safety Engineering, 2017, 34(3):444-451.
[12] 郭惟嘉, 王海龙, 陈绍杰, 等.采动覆岩涌水溃砂灾害模拟试验系统研制与应用[J].岩石力学与工程学报, 2016, 35(7):1415-1422.
GUO Weijia, WANG Hailong, CHEN Shaojie, et al.Development and application of simulation test system for water and sand inrush across overburden fissures due to coal mining[J].Chinese Journal of Rock Mechanics and Engineering, 2016, 35(7):1415-1422.
[13] 贾后省, 马念杰, 赵希栋.浅埋薄基岩采煤工作面上覆岩层纵向贯通裂隙“张开—闭合”规律[J].煤炭学报, 2015, 40(12):2787-2793.
JIA Housheng, MA Nianjie, ZHAO Xidong.“Open-close” law of longitudinal transfixion cracks in shallow buried coal face with thin bedrock[J].Journal of China Coal Society, 2015, 40(12):2787-2793.
[14] CUNDALL P A, STRACK O D L.Discussion:A discrete numerical model for granular assemblies[J].Geotechnique, 1980, 30(3):331-336.
[15] Itasca Consulting Group Inc PFC3D.Particle flow code in 3 Dimensions, version 5.0[M].MN,USA; Itasca Consulting Group Inc,2017.
[16] 肖 成.基于颗粒流的弱胶结裂隙岩体水沙两相流特性研究[D].徐州：中国矿业大学, 2016.
[17] 钟江城, 周宏伟, 赵宇峰, 等.浅埋煤层开采溃水溃砂两相流的耦合数值研究[J].工程力学, 2017, 34(12):229-238.
ZHONG Jiangcheng, ZHOU Hongwei, ZHAO Yufeng,et al.The two-phase flow of water-sand inrush under shallow coal seam mining:a coupled numerical study[J].Engineering Mechanics, 2017, 34(12):229-238.
[18] 王 刚, 王 锐, 武猛猛, 等.火区下近距离煤层开采有害气体入侵灾害防控技术[J].煤炭学报, 2017, 42(7):1765-1775.
WANG Gang, WANGRui, WU Mengmeng, et al.Prevention and control technology of harmful gas intrusion in close-up coal seam under fire area[J].Journal of China Coal Society, 2017, 42(7):1765-1775.
[19] WANG Tao, ZHOU Weibo, CHEN Jinhua, et al.Simulation of hydraulic fracturing using particle flow method and application in a coal mine[J].International Journal of Coal Geology，2014,121：1-13.
[20] 石 崇, 徐卫亚.颗粒流数值模拟技巧与实践[M].北京:中国建筑工业出版社，2015:24.
[21] WANG Gang, FAN Cheng, XU Hao, et al.Determination of long horizontal borehole height in roofs and its application to gas drainage[J].Energies,2018, 11：2647.
[22] 武猛猛, 王 刚, 王 锐, 等.浅埋采场上覆岩层孔隙率的时空分布特征[J].煤炭学报, 2017, 42(S1):112-121.
WU Mengmeng, WANG Gang, WANG Rui, et al.Space-time porosity distribution on overlying strata above a shallow seam[J].Journal of China Coal Society, 2017, 42(S1):112-121.
[23] 王 锐, 修 毓, 王 刚, 等.基于颗粒流理论的煤与瓦斯突出数值模拟研究[J].山东科技大学学报：自然科学版, 2016,35(4):52-61
WANG Rui, XIU Yu, WANG Gang,et al.Numerical simulation of coal and gas outburst based on particle flow theory[J].Journal of Shandong University of Science and Technology：Natural Science, 2016,35(4):52-61.
