0 引 言
随着东部矿区煤炭资源的枯竭，中部受资源与环境约束的矛盾加剧，我国煤炭资源开发重点已快速向西部地区转移。西部地区煤层赋存的主要特点是厚煤层与特厚煤层分布广泛，煤炭资源储量丰富，但这些煤层多为倾斜煤层且煤层透气性较低[1]。煤层采掘后，下伏煤层经历了应力集中、开采卸压、采空区应力恢复等应力变化，并随之产生不同形式的塑性破坏。而倾斜特厚煤层分层开采带来的大尺度煤层活动空间，导致下伏煤层的应力、塑性破坏区分布规律显示出与常规薄、中厚煤层不同的规律。因此掌握采动后下伏煤层的应力分布规律及破坏特征能够为精准判别特厚煤层采空区瓦斯富集区提供一定理论指导。
目前，围绕煤层下伏煤岩应力分布及破坏特征，国内外专家学者做了大量的研究工作，取得了丰硕的理论成果。前苏联B.Д.СnесареВ[2]最早对工作面底板岩层采动变形破坏规律进行了研究，李树刚等[3]研究了上保护层开采后的下伏煤岩体应力、渗流等变化规律。文献[4-5]建立了大倾角煤层走向长壁开采岩体结构空间模型，深入分析了非对称煤岩的破坏失稳过程；蒋力帅等[6]通过工作面支承压力场演化特征的实例仿真分析实现了采动应力场与采空区压实承载作用的耦合分析；程志恒等[7]研究了采动条件下工作面顶底板应力的卸压程度以及覆岩断裂带的动态演化特征，文献[8-9]通过分析底板不同位置垂直应力，得到了采空区底板不同深度岩体裂隙演化规律。虽然从工作面采深[10-12]、承压水上底板扰动[13-14]、采空区压实效应[15-16]等方面对底板应力分布特征开展了大量研究，但煤层倾角作为对采空区底板应力分布及破坏特征的重要影响因素，尤其是在特厚煤层分层开采这一工程背景下，其作用规律有待于进一步研究。
以新疆某特厚煤层主采工作面为例，首先采用岩体质量评分法对煤层顶底板岩性进行分析，确定其岩体力学参数。进而采用FLAC3D建立不同倾角的物理模型，对工作面开采进行模拟计算，分析煤层倾角对下伏煤层应力分布及破坏特征的影响规律，为精准判定采空区瓦斯富集区提供一定理论依据。
1 数值模拟方案
1.1 模拟原型及岩体力学参数确定
新疆某特厚煤层矿井的主采工作面平均埋深380 m，设计走向1 072 m，倾向长度104.7 m，煤厚32.6～36.9 m，平均厚33.3 m，煤层倾角21°～25°，自然安息角35°。工作面采用走向长壁后退式综合机械化放顶煤采煤法，采高3 m，放顶6 m，采放比为1∶2，下伏煤层厚24 m。
为模拟不同煤层倾角下伏煤层应力分布差异，以开采层赋存条件为基础，煤层倾角设定为0°、15°、30°。煤岩层岩体力学参数首先由实验室测试获取，进而采用Hoek-Brown强度准则计算得到，见表1。
表1 模型各层岩体力学参数
Table 1 Rock mechanics parameters of each seam
岩性体积模量/GPa剪切模量/GPa黏聚力/MPa内摩擦角/(°)抗拉强度/MPa密度/(kg·m-3)煤3.111.602.2029.501.611 460泥岩2.711.552.9533.002.872 080细砂岩7.503.463.6534.005.102 620粉砂岩9.072.964.1935.502.212 840中砂岩3.391.754.2736.003.602 660粗砂岩4.861.276.0437.003.402 430
1.2 不同煤层倾角物理模型构建
根据试验矿井实际的生产技术条件，数值计算中物理模型平均埋深设定为380 m，工作面长度设定为100 m，推进总距离设定为400 m，每步推20 m。
煤层开采后，根据应力的影响范围及围岩破坏情况，可将煤层化分为塑性区、弹性区和原始岩层应力区。如果模型边界取在煤层弹、塑性区域内，边界受采动影响将产生水平方向位移，而物理模型边界采用的是位移约束条件，可能会导致岩体应力模拟结果产生较大误差。当模型侧边界取在原始岩层应力区时，受开采影响相对较小，则得到的模拟结果可以认为接近真实情况[13]。
因此对于物理模型尺寸，可通过计算弹、塑性区沿走向及倾向的长度，确定物理模型沿X(走向)、Y(倾向)方向的长度；通过计算裂隙高度以及底板破坏深度可确定物理模型沿Z(垂向)方向的长度。以倾角为0°为例，利用A.H.威尔逊理论公式计算塑性区长度为
s=0.004 92MH=16.83 m
(1)
其中：H为采深；M为采厚。为确保计算正确性，塑性区长度s取20 m，采场沿走向前后两端、左右两端各预留10倍塑性区长度[17]，即物理模型长800 m，宽500 m。
根据《煤矿防治水规定》选取“两带”最大高度的经验公式，计算垮落带Hc、断裂带高度Hf，见式(2)。
(2)
根据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程》中给出工作面底板破坏深度h计算公式。
h=0.008 5H+0.166 5α+0.207 9L-4.357 9
(3)
其中：α为煤层倾角；L为工作面长。计算得到工作面底板破坏深度为19.66 m。为确保计算正确性，断裂带顶部预留约50 m，底板破坏带底部预留50 m，因此Z方向长度为200 m。通过上述计算，确定物理模型的长、宽、高分别为800、500、200 m。同理，可计算其他煤层倾角的物理模型尺寸，见表2，构建模型时各煤岩层厚度保持一致，不同煤层倾角采场物理模型如图1所示。
表2 不同煤层倾角物理模型尺寸
Table 2 Physical model sizes of different coal seam dip angles
煤层倾角/(°)长/m宽/m高/m顶部荷载/MPa08005002006.25158005353234.75308005334233.55
L1—倾角0°、15°、30°时分别为200、323、423 m；L2—倾角0°、15°、30°时分别为500、535、533 m
图1 不同煤层倾角采场物理模型
Fig.1 Physical model of stope with different coal seam inclination
1.3 模型边界条件及本构模型的选取
不同倾角模型的前后左右4个边界均采用位移约束，模型上部边界为自由边界，模型下部边界为位移全约束边界，即Z=0 m平面上的节点，沿X、Y、Z三向固定；模型顶端施加等效荷载见表2，模拟上覆岩层至地表土层自重压力，重力加速度取9.8 m/s2；为更好地贴近真实采场矿压显现规律，模型在开采煤层顶板设置接触面以防止顶板压入底板。
模型初始化时，煤岩材料屈服破坏判据采用Mohr-Coulomb强度准则，见式(4)；该准则认为煤岩中某个平面上的法向应力影响煤岩抗剪强度，当该切面上的法向应力与剪应力符合某种组合时，煤岩将在切面上产生断裂带；当断裂带贯通时，即煤岩发生破裂[18]，相应的发生剪切屈服的条件为fs<0。
(4)
式中：c为煤岩体黏聚力；φ为内摩擦角；σ1、σ3为第1、第3主应力。
由于FLAC3D模拟的是连续介质，无法真实反映垮落岩体与母岩分离的过程，因此比较符合基本顶初次来压前的采空区情况；当回采距离大于初次来压步距后，随着采空区的不断扩大，垮落岩体会逐渐充满整个回采空间，在此之后，垮落岩体受上覆岩层压力作用而被逐渐压实，支承能力不断提高，最终形成支承上覆岩层的支承体。因此回采至初次来压前，可采用直接开挖法；而回采距离超过初次来压步距之后，直接开挖法会导致回采空间中缺少与采空区底板直接接触垮落岩体，导致上覆岩层的重力荷载无法有效传递至采空区底板，产生较大的计算误差。可在计算过程通过识别拉伸破坏区来确定垮落带范围，进而挖去整个垮落带，然后对垮落带及采空区进行充填处理，即可使采空区顶底板连接起来，进而可以针对采空区底板应力分布及破坏情况进行有效模拟计算。
由上述分析可知，煤层的推进模拟可通过开挖→识别垮落带范围→挖去垮落带→充填垮落带及采空区→下一步开挖的模拟计算步骤来实现，因此需首先确定垮落带破碎岩体的充填模型及岩体力学参数。FLAC3D中的双屈服模型，考虑了剪切及拉涨屈服函数，同时加入了模型体积屈服函数，比较适合模拟破碎岩体的特性，可作为破碎岩体的本构模型。
2 破碎岩体的力学参数
对于破碎岩体的力学参数，可首先基于采空区压实理论，对其应力-应变特性进行分析，进而构建简单块体模型，通过反演迭代并与理论解进行匹配来获得力学参数。
2.1 采空区压实理论
采空区垮落矸石的渐进压实过程及其承载力学特性对采空区与底板煤体应力状态有着极其重要的影响[6]。对松散破碎岩体压实过程中的岩体力学响应特性，SALAMON[19]提出的破碎岩体压缩过程的应力-应变关系，见式(5)，被国内外研究人员广泛认可并应用[20-22]。
(5)
其中： σ为垮落岩体的所受垂直应力；E0为岩体初始切线模量；ε为垂直应力作用下垮落岩体的应变；εm为岩体碎胀后最大应变，可由碎胀系数b，通过式(6)求得。
εm=(b-1)/b
(6)
E0可以由式(7)计算获得：
(7)
其中，σc为单轴抗压强度。由式(5)—式(7)得垮落岩体应力-应变表达式
(8)
2.2 采空区垮落岩体岩体力学参数
依据该煤矿的实际条件可知碎胀系数b=1.349，单轴抗压强度σc=43.91 MPa，代入式(6)至式(8)，计算得最大应变εm=0.26、参数E0=54.72 MPa以及压实过程中的应力-应变关系，见式(9)。
(9)
为得到合理的破碎岩体参数且使其应力-应变关系能够与式(9)匹配，首先建立长、宽为1 m、高为2 m的简单块体模型，赋予双屈服模型，底部固支，顶部自由移动，其他4个面约束横向位移。在其顶部表面施加竖向速率模拟加载过程，通过迭代反演获得拟合度比较合适的岩体力学参数。经反复模拟计算，得到破碎岩体岩体力学参数见表3，反演匹配结果如图2所示。
表3 破碎岩体岩体力学参数
Table 3 Mechanical parameters of backfilling rock mass
泊松比体积模量/GPa剪切模量/GPa黏聚力/MPa内摩擦角/(°)0.230.290.080.1017
图2 岩体力学参数反演结果
Fig.2 Inversion results of mechanical parameters
3 数值模拟结果分析
3.1 煤层倾角变化对下伏煤层应力分布的影响
工作面推进400 m时，煤层倾角为0°、15°、30°的采场底板应力分布如图3所示，为优化显示效果，依据对称性原理，选取X为200～400 m。分析可知，红色区域为采空区周围邻近巷道的卸压区，采空区内部的绿色区域为应力恢复区，而工作面及采空区两侧两巷道外部蓝色区域则为应力集中区。由此可见，工作面开采使围岩应力重新分布，沿走向在工作面前方产生超前支承压力区，沿倾向在采空区两侧产生应力集中区，采空区内部则在上覆岩层作用下压实，出现应力回升的现象。
随着煤层倾角增大，倾向上部卸压区范围增大，倾向下部卸压区范围减小，工作面前方支承应力峰值逐渐增大，采空区应力恢复“O”形圈不再以采空区中部走向线对称分布，而向倾向下部偏移，当煤层倾角为30°时，偏移距离达到最大为19 m。这是由于垮落后的矸石沿煤层底板滑移会造成回采空间非均匀充填，且煤层倾角越大，充填倾向中下部的岩块堆砌的高度越高，产生较强的支承作用并导致垮落空间的相对减少，垮落不充分，而受载采空区底板煤体受力较大；倾向上部区域顶板垮落空间则相对较充足且无法得到垮落岩体的有效支承，因而垮落较为充分，而受载采空区底板煤体受力较小，因此整体上呈现出“O”形圈中部偏离工作面中部，且煤层倾角越大，偏移程度越大(图3)。
图3 不同倾角采场底板应力分布
Fig.3 Stress distribution diagram with different inclination angles
为获得准确的采空区底板应力分布规律，使用FLAC3D内置hist函数，对采空区底板不同深度煤层应力进行监测。测线(监测点构成的线)布置如图4所示，其中，α为煤层倾角。
图4 监测线布置
Fig.4 Monitoring line layout
采场上部、中部及下部位置(即测线AA1、BB1、CC1)提取的垂直应力分布如图5所示，工作面推进方向为正方向。可见，工作面前方和采空区的垂直应力分布随着煤层倾角的变化存在差异。煤层倾角变化对工作面前方支承应力影响显著，随着煤层倾角增大，在倾向上、中部，工作面超前支承应力集中系数呈现减小趋势，在倾向下部则呈现增大趋势。
同时，图5可反映出采空区由于垮落带压实造成的应力恢复现象。采空区应力恢复在走向上有着显著差异，距离工作面较远的采空区中部应力恢复程度较高，距离采空区边界越近即越靠近工作面应力恢复程度越低。随着煤层倾角增大，在采场倾向上、中部，采空区内部应力恢复程度均逐渐减小，而在采场倾向下部，采空区内部应力恢复程度则逐渐增大。
图5 沿走向采场底板应力分布
Fig.5 Stope floor stress distribution along strike
相较于水平煤层工作面，煤层倾角较大的工作面中部首先开始裂断和垮落，接着是工作面的倾向上部，然后是工作面倾向下部。与此同时，由于垮落岩体对采空区下部顶板产生有效支承而对上部顶板的支承出现效果较弱甚至“负约束效应”，引发了采空区中上部底板受上覆岩层荷载较小，下部的荷载较大的非均衡现象。因此随着煤层倾角增大，采空区下部和下端侧向煤柱承担了较多的上覆岩层的压力，采空区则产生沿倾向中上部应力恢复程度减小，沿倾向下部应力恢复增大的现象。
3.2 煤层倾角对底板不同深度煤层应力分布的影响
提取图4中采空区底板下方1、12、24 m测线数据(即测线DD1、EE1、FF1)，得到下伏煤层不同深度煤体应力分布，如图6所示(倾向下端为正方向)。
图6 下伏煤层沿倾向不同深度应力分布
Fig.6 Stress distribution of different depths along tendency of underlying coal seam
采场底板下方24 m处相较于12、1 m处垂直应力分布整体更加趋于平缓，卸压程度更小，但煤层倾角的变化对应力分布影响依然存在。由图6可知，倾向上端侧向应力集中系数均小于倾向下端侧向应力集中系数；随着煤层倾角增大，倾向上、下端侧向应力集中系数呈减小趋势，采空区应力恢复程度则增大。
由此可知，随着埋深增加，倾向上下端侧向应力集中系数降低，采空区内部应力恢复程度增大，整体更加趋于原始应力。这是因为随着深度的增加，受煤层倾角影响，倾向上下端煤体不同埋深逐渐成为影响垂直应力分布的主要因素，下端煤体埋深大于上端煤体，因而呈现下侧压力大于上侧压力。采动过后的下伏煤层虽然仍处于卸压状态，但在围岩支承和采空区压实承载作用下回升至原始应力状态。
3.3 煤层倾角变化对下伏煤层塑性破坏特征的影响
图7为倾角0°、15°、30°的工作面推进400 m时，沿走向和倾向顶底板塑性破坏区分布。
图7 不同煤层倾角工作面推进400 m时顶底板塑性区分布
Fig.7 Plastic distribution of roof and floor when working face with different coal seam inclination is pushed 400 m
由图7可知：工作面回采后，在采空区底板煤体内产生1个大的塑性破坏区，在工作面下方及倾斜上、下端底板煤体内各自产生1个具有一定破坏深度和范围的塑性区，且沿倾向下端底板塑性破坏区的深度与采空区底板岩体内的塑性破坏区深度相接近。
随着煤层倾角增大，采空区底板煤体内最大塑性破坏深度位置逐渐远离采空区中部，向倾向下端偏移，当煤层倾角为0°、15°、30°时，偏离距离分别为0、9.8、15.6 m，呈现增大趋势。
如图8所示，沿倾向上下端、工作面及采空区底板煤体内塑性破坏深度均呈现随煤层倾角增大而增大趋势，这是由于底板煤体在30°左右时，所受到的剪切应力最大，更容易发生滑移破坏。沿倾向，在采空区两侧的侧向支承压力的作用下，形成1个沿倾向长度大于工作面长度塑性破坏区；由图8反映出倾向上端破坏深度始终小于下端破坏深度，但倾向上端和下端破坏深度均随着煤层倾角增大而增大。
图8 工作面及沿倾向底板煤体塑性区破坏深度
Fig.8 Failure depth of plastic zone of coal body in working face and along dip
4 结 论
1)基于采空区压实理论获得了较为精确的垮落岩体的应力-应变关系，通过迭代反演获得了破碎岩体的岩体力学参数。
2)垮落后的矸石沿煤层底板移动会造成回采空间非均匀充填，引发了采空区中上部底板受上覆岩层荷载较小，下部的荷载较大的非均衡现象。由于采空区下部和下端侧向煤柱承担了较多的上覆岩层的压力，因此随着煤层倾角增大，倾向上部的卸压区范围增大而倾向下部卸压区范围减小，采空区应力恢复“O”形圈不再以采空区中部走向线对称分布，而向倾向下部偏移。
3)对于下伏煤层不同深度煤体，倾向上下端侧向应力集中系数都呈现随煤层倾角增大而减小的趋势，且倾向上端侧向应力集中系数均小于倾向下端侧向应力集中系数；采空区应力恢复程度在下伏煤层不同深度处均随煤层倾角增大而增大。受煤层倾角影响，沿倾向，上下端煤体的埋深成为影响垂直应力分布的主要因素，呈现倾向下端压力大于倾向上端压力。
4)分析了下伏煤层塑性区分布规律。工作面回采后，采空区及采空区四周巷道底板均产生具有一定范围和破坏深度的塑性破坏区。各个区域内的塑性破坏区深度均随煤层倾角的增大而增大，其中上端破坏深度始终小于下端破坏深度，采空区底板最大塑性破坏深度所在位置随倾角增大逐渐远离采空区中部，向倾向下端偏移。
[1] 李树刚，徐培耘，林海飞，等.倾斜煤层卸压瓦斯导流抽采技术研究与工程实践[J].采矿与安全工程学报，2020,37(5):1001-1008.
LI Shugang,XU Peiyun,LIN Haifei,et al.Technology research and engineering practice of pressure-relief gas diversion extraction in inclined coal seams[J].Journal of Mining and Safety Engineering，2020,37(5):1001-1008.
[2] 周世宁.煤层瓦斯赋存与流动理论[M].北京:煤炭工业出版社,1999.
[3] 李树刚,索 亮,林海飞,等.不同间距上保护层开采卸压效应UDEC数值模拟[J].辽宁工程技术大学学报(自然科学版),2014,33(3):294-297.
LI Shugang,Suo Liang,LIN Haifei,et al.Numerical simulation of pressure relief effects of protective layer mining in different distances[J]. Journal of Liaoning Technical University(Natural Science Edition),2014,33(3):294-297.
[4] 罗生虎,伍永平,刘孔智,等.大倾角大采高综采工作面煤壁非对称受载失稳特征[J].煤炭学报,2018,43(7):1829-1836.
LUO Shenghu,WU Yongping,LIU Kongzhi,et al.Asymmetric load and instability characteristics of coal wall at large mining height fully-mechanized face in steeply dipping seam[J]. Journal of China Coal Society,2018,43(7):1829-1836.
[5] 伍永平,胡博胜,解盘石,等.大倾角工作面飞矸冲击损害及其控制[J].煤炭学报,2018,43(10):2694-2702.
WU Yongping,HU Bosheng,XIE Panshi,et al.Impact damage of flying gangue in steeply dipping seams and its control[J]. Journal of China Coal Society,2018,43(10):2694-2702.
[6] 蒋力帅,武泉森,李小裕,等.采动应力与采空区压实承载耦合分析方法研究[J].煤炭学报，2017,42(8):1951-1959.
JIANG Lishuai,WU Quansen,LI Xiaoyu.Numerical simulation on coupling method between mining-induced stress and gob compression [J]. Journal of China Coal Society,2017,42(8):1951-1959.
[7] 程志恒,齐庆新,李宏艳,等.近距离煤层群叠加开采采动应力-裂隙动态演化特征实验研究[J].煤炭学报,2016,41(2):367-375.
CHENG Zhiheng,QI Qingxin,LI Hongyan,et al.Evolution of the superimposed mining induced stress-fissure field under extracting of close distance coal seam group[J]. Journal of China Coal Society，2016,41(2):367-375.
[8] 张春雷.煤层群上行开采层间裂隙演化及卸压空间效应[D].北京：中国矿业大学(北京),2017.
[9] 郭 良,张春雷.采空区底板不同深度岩体裂隙演化规律[J].煤矿安全，2019,50(2):204-207.
GUO Liang,ZHANG Chunlei.Fracture evolution law of rock mass at different depths of floor in gob[J]. Safety in Coal Mines,2019,50(2):204-207.
[10] 翟晓荣.矿井深部煤层底板采动效应的岩体结构控制机理研究[D].淮南：安徽理工大学,2015.
[11] 李家卓,谢广祥,王 磊,等.深部煤层底板岩层卸荷动态响应的变形破裂特征研究[J].采矿与安全工程学报,2017,34(5):876-883.
LI Jiazhuo,XIE Guangxiang,WANG Lei,et al.Mechanical mechanism of dynamic fracture evolution of coal floor unloading in deep mining[J].Journal of Mining and Safety Engineering,2017,34(5):876-883.
[12] 苗 飞.大采深工作面煤层底板采动破坏及特征研究[D].石家庄：河北工程大学,2017.
[13] 李振峰.承压水上采空区底板岩层扰动规律研究[D].焦作：河南理工大学,2015.
[14] 宋文成,梁正召,赵春波.承压水上开采沿工作面倾向底板力学破坏特征[J].岩石力学与工程学报,2018,37(9):2131-2143.
SONG Wencheng,LIANG Zhengzhao,ZHAO Chunbo. Mechanical failure characteristics of mining floor along working face inclination above confined water[J].Chinese Journal of Rock Mechanics and Engineering,2018,37(9):2131-2143.
[15] 王 猛,李文峰.基于Salamon模型的采空区冒落矸石加载机理分析[J].煤炭科技,2015,41(2):50-54.
WANG Meng,LI Wenfeng.Analysis on the loading mechanism of caving gangues in gob based on Salamon model[J].Coal Science and Technology Magazine,2015,41(2):50-54.
[16] 王文学,王四巍,刘海宁,等.采后覆岩裂隙岩体应力恢复的时空特征[J].采矿与安全工程学报,2017,34(1):127-133.
WANG Wenxue,WANG Siwei,LIU Haining,et al.The space and time characteristics of the cover stress re-establishment of the fractured rock mass in the goaf after coal mining[J]. Journal of Mining and Safety Engineering,2017,34(1):127-133.
[17] 钱鸣高,石平五,许家林.矿山压力与岩层控制[M].徐州：中国矿业大学出版社,2010.
[18] 连民杰,李俊平.矿山岩石力学[M].北京：冶金工业出版社,2011.
[19] SALAMON M D G. Mechanism of caving in longwall coal mining[C]//Proc. of the 31st US. Symposium of Rock Mechanics. Golden West Virginia University Press,1990:161-168.
[20] SAEEDI G,SHAHRIAR K,REZAI B.Numerical modelling of out-of-seam dilution in longwall retreat mining[J].International Journal of Rock Mechanics and Mining Sciences，2010,47(4):533-543.
[21] XUE J,WANG H,ZHOU W,et al. Erratum to:Experimental research on overlying strata movement and fracture evolution in pillarless stress-relief mining[J]. International Journal of Coal Science & Technology，2016,3(3):350.
[22] 白庆升,屠世浩,袁 永,等.基于采空区压实理论的采动响应反演[J].中国矿业大学学报，2013,42(3):355-361.
BAI Qingsheng,TU Shihao,YUAN Yong,et al.Back analysis of mining induced responses on the basis of goaf compaction theory[J].Journal of China University of Mining and Technology，2013,42(3):355-361.
O