大采深综放开采覆岩移动规律离散元数值模拟研究

余学义，王昭舜，杨 云，毛旭魏

( 1.西安科技大学 能源学院，陕西 西安 710054；2.西安科技大学 教育部西部矿井开采与灾害防治重点实验室，陕西 西安 710054 )

实现科学采矿已然成为新时期煤矿开采的必然要求[1-2]。近年来，随着煤矿机械化程度的不断提高，不少煤矿已采用大采深综放开采的生产方式，高强度的资源开发必然带来开采沉陷等破坏。故而大采深综放开采地表移动变形规律的研究就显得尤为重要，而此项研究的首要条件就是了解上覆岩层的移动变形规律。

早在20世纪初，便有专家通过简单的力学原理解释了一些矿山压力现象，并初步形成了矿山压力假设。随着压力拱假说，悬臂假说和铰接岩石假说的发展，矿山压力理论已经形成了一个系统。在19世纪80年代，钱鸣高[3-4]等提出控制全部或部分上覆地层的地层称为关键地层；茅献彪[5]等研究了采场覆盖层中关键层的破坏规律和关键层负荷的变化规律；许家林[6-12]等提出了基于关键层控制的局部充填开采技术以及影响深部开采关键层对地表沉陷影响的一些因素；黄庆享[13-15]等研究了厚沙层地质条件下关键层上部的载荷分布，传递规律以及破断规律；余学义[16-18]等对大采高覆岩灾变规律进行了研究，分析了影响下沉的采矿因素，研究得出大采深综放开采地表移动变形规律；李化敏[19]等研究了特厚煤层综放开采中的放煤理论。

众多学者经过长期研究，还针对近距离厚煤层组工作面覆岩破坏规律、特厚煤层不同开采阶段“弱-弱”结构覆岩破坏高度、厚煤层综放开采支架围岩关系、顶板运动、窄煤柱围岩控制及含夹矸特厚煤层综放工作面顶煤破碎机理等开展了大量研究[20-21]。

到目前为止，众多学者已经进行了大量研究工作，但是大采深综放开采条件下，随着工作面开采数量的增加，地表从极不充分-非充分-充分采动过程中覆岩的演化规律对地表移动变形的影响鲜有深入研究。本文基于陈家沟煤矿特殊的地质条件，研究厚埋深大采高综放多工作面开采过程中的覆岩演化规律以及对地表移动变形的影响。

1 试验矿井地质条件

通过多个钻孔探测，采区地层自上而下有表土层的第四系全新统，总厚度为268.75 m的上第三系甘肃群，涵盖直罗组与延安组的中侏罗统，下侏罗统富县组与上三叠统延安群。详情如图1，表1所示。

陈家沟煤矿主采5煤层，工作面平均宽120 m，煤层倾角6°，分层开采高度12 m。平均开采深度为506 m，最大开采深度达到580 m，即大采深( 煤层埋深500～1 000 m )煤矿[22]。

图1 岩层柱状图Fig.1 Histogram of rock formation

表1 陈家沟煤矿井田地层Table 1 Field layer profiles of Chenjiagou Mine

根据钻孔取芯分析，按岩性分类，上覆岩为中硬岩性，且泥岩厚度占较大比例。由于泥岩本身具有很好的防水性，其已成为覆岩中良好的隔水层。

2 地表移动观测结果分析

工作面测点布置如图2所示。8512工作面累计回采长度2 235 m，走向观测线共计35个有效工作测点，点号分别为Z1～Z35，倾向观测线共计35个工作测点，点号分别为Q1～Q35；8513工作面累计推进1 710 m，走向观测线共计118个工作测点，点号分别为zx1～zx118，倾向观测线共计40个工作测点，点号分别为qx1～qx40。

8512工作面自2013-11—2015-11观测结束共计观测15次。8513工作面从2016-05-16开始，截止到2019-01-13，累计观测23次。8512工作面走向最大下沉达到1.124 m，8513工作面走向最大下沉达到2.226 m，由于2个工作面走向均属超充分采动，故将倾向地表下沉作为研究重点，现对2个工作面倾向的观测数据进行分析。

图2 工作面测点布置Fig.2 Layout of measuring points in working face

8512工作面倾向开采宽度120 m，开采宽深比小于1/4，属于极不充分采动，根据多次观测数据，绘制倾向观测点下沉曲线如图3所示，未出现下沉盆地。8513工作面根据观测数据，绘制观测点下沉曲线如图4所示。倾向累计工作面宽度为246 m，属于非充分采动。

图3 8512工作面倾向下沉曲线Fig.3 Sinking curves of 8512 working face

图4 8513工作面倾向下沉曲线Fig.4 Sinking curves of 8513 working face

通过分析观测数据可知，在8512首采工作面开采后，下沉变化过程呈规律性增大且下沉量均不大。根据8513工作面观测结果发现，在qx24点、qx27点处下沉发生突变，而其正处于8512工作面的采空区上覆岩层，说明8513工作面的开采对8512工作面采空区产生扰动，导致8512工作面上覆岩层发生活化，造成下沉量增加。在地表发育过程中由于采深较大，当回采结束后，采空区内围岩应力处于新的平衡，处于稳定阶段；当再次开采，扰动范围内的工作面会打破这种应力平衡，已达到平衡的采空区会被活化，从而导致二次沉降，这与观测数据吻合，故而在大采深开采过程中，地表达到稳态的周期会较长。

由此得出：8512工作面开采后由于是单一工作面开采，此时地表处于极不充分采动状态，受关键层支承上部压力的影响，传递到地表表现出沉降量较小，下沉不明显，沉降影响范围也较小。在8513工作面开采后，地表下沉逐渐增大，且部分区域产生扰动，此时地表处于非充分采动状态。故地表由极不充分采动向非充分采动发育时，下沉会有突变，且下沉速度突变，即是关键层断裂失稳所致。

此外，利用地表移动盆地的主断面上实测下沉曲线求得地表移动角参数。具体结果详见表2。

表2 实测角参数Table 2 Measured angular parameters

3 覆岩结构控制条件分析

3.1 关键层分析

该矿区煤层覆盖层的综合岩性为中硬型。较硬的地层主要分布在侏罗纪地层。上第三系甘肃群地层一般较软，结构疏松，抗压强度低，遇水易开裂和风化；第四系地层薄而软；特厚煤层中煤5层的顶部为砂岩，属于较厚的硬岩地层；在中侏罗统第二和第三层中，有几层较厚且坚硬的岩性，这些地层在覆盖层的运动和变形中起着控制作用。使用固定支撑梁的力学模型来分析关键层的应力状态。

设其挠曲方程为

则满足边界条件：

利用伽辽金法通过解算，可得固支梁的最大挠度为

式中，I为惯性矩；E为岩体的弹性模量，GPa；q为岩体上的载荷，kN/m2。

3.2 关键层的确定

上覆岩层的荷载q系由岩梁自重及其上覆岩层共同作用，因此计算各岩层承载的荷载可按式( 4 )计算，详细参数见表3。

表3 钻孔部分岩层力学参数Table 3 Mechanical parameters of some rock strata in borehole

式中，E1，E2，…，En为各层岩层的弹性模量，MPa；h1，h2，…，hn为各层岩层的厚度，m；γ1，γ2，…，γn为各层岩层的容重，kN/m2；( qn)1为第1层岩层所受上方第2层至第n层的载荷，kPa。

计算结果见表4。

表4 上覆岩层承载载荷Table 4 Bearing load of overlying strata

计算得出坚硬岩层为第3层砾岩、第6层砂质泥岩、第11层细粒砂岩和第19层粉砂岩，并计算破断距综合得出第19层为关键层1，第11层为关键层2，关键层3为第3层砾岩。

计算结果表明，8512工作面开采后关键层2以下的岩层均断裂破坏，关键层2起着控制覆岩沉陷的作用，但是关键层2以上岩层在移动过程中还会出现离层结构，到达地表的沉陷量小于关键层2的挠曲量。

4 覆岩移动演化规律模拟

4.1 模拟试验数据分析

数值模拟采用3DEC块体离散元数值模拟软件，模拟基于陈家沟煤矿八采区地质采矿条件的多工作面开采覆岩移动变形规律。

结合地质资料和岩石力学试验结果，得出各岩层物理力学参数并采用莫尔-库仑塑性模型来进行计算,得到

式中，K为体积模量；E为杨氏模量；G为剪切模量；v为泊松比。

本模型岩层节理面材料模型采用库仑滑动模型( 完全弹塑性 )，各节理面物理力学参数见表5。

表5 各节理面力学参数Table 5 Mechanical parameters of each joint surface

续 表

4.2 试验模型建立

为了研究陈家沟煤矿大采深综放开采覆岩移动变形规律，采用3DEC三维计算模型，模型具有6个自由面，水平方向和底部边界采取薄层的方法进行速度为零的固定。模型尺寸为1 000 m×30 m×570 m，本文模拟单一煤层开采，煤层埋深525 m，煤层为水平煤层，煤层厚度26 m，设计4个开采工作面，分别为8512，8513，8511，8514工作面，工作面宽度均为120 m，留设区段煤柱6 m，各单元在运动时一般不发生弹跳，这是由于运动时动能转化成热能耗散掉的缘故，因此煤岩块的运动是不可逆的过程，为避免岩块在平衡位置振动，应加入阻尼来吸收系统的动能，本文采用了自适应阻尼。三维数值计算模型如图5所示。

图5 3DEC三维计算模型Fig.5 3D computing model of 3DEC

4.3 模拟结果分析

4.3.1 851 2工作面模拟结果分析

在模拟过程中，在4个工作面的覆岩和地表运动稳定后，可以获得覆岩移动变形图、地表沉降云图和测量线数据。通过覆岩移动变形图和云图，可以获得覆岩和地表沉降的分布，地表的最大沉降点，地表沉降的范围和裂隙破裂带的范围；地表沉降曲线可以通过对地表测量线数据的分析得出。

关键层1初次破断以及周期破断时块体的位置变化，如图6所示，随着8512工作面的推进，垮落高度呈上升趋势，开采至10 m时，垮落高度6 m；开采至20，30 m 时，垮落高度分别增加到14.50 m 和23.40 m。工作面持续推进，直接顶继续垮落，由于岩层间的互相挤压，当开采至70 m时，关键层1初次破断，此时监测到地表下沉0.54 m；继续开采至27.10 m，失稳现象严重，第1次周期来压，由图6可知，下部岩层产生了约5 m的超前断裂，此时监测到地表下沉0.77 m。

图6 关键层1首次破断Fig.6 First break of the key stratum 1

模拟8512工作面开采结束，如图7所示，开采尺寸为120 m，平均采高9 m，工作面宽深比为0.24，此时地表处于极不充分采动状态。覆岩中关键层1发生破断，破断距分别为70.0 m和27.1 m。从煤层顶板上方向上33.4 m为覆岩垮落带，再向上76.9 m的范围为裂缝带，即裂缝带发育高度为110.3 m，裂采比为12.2，裂缝带向上直达地表范围为覆岩的弯曲下沉带。由于开采范围小，开采对覆岩扰动破坏影响到关键层2下部，关键层2以下的岩层均断裂破坏，关键层2表现为微量的挠曲变形，但是关键层2以上岩层在移动过程中还会出现离层结构，其上覆岩土层表现为弯曲下沉。根据模拟观测结果表明在开采一个工作面后，覆岩移动变形规律与关键层2挠曲形态吻合，再次证明，地表处于极不充分采动状态是由于关键层起控制作用。

图7 8512工作面开采结束Fig.7 End of 8512 working face mining

4.3.2 851 3工作面模拟结果分析

8513工作面在开采过程中，关键层2初次破断以及周期破断时块体的位置变化如图8所示，工作面在持续推进过程中，上覆岩层向上发育迅速。当开采至146.6 m时，关键层2底部初次破断，整体表现为局部破断的挠曲变形，此时地表下沉1.18 m；继续推进70.2 m，累计推进213.8 m时，关键层2周期破断，此时监测地表下沉1.22 m。工作面初次来压是70.0 m，周期来压是21.7 m，这与现场实测及理论计算结果均基本吻合。

图8 关键层2首次破断Fig.8 First break of the key stratum 2

模拟8513工作面开采结束，如图9所示，开采尺寸为240 m，平均采高12 m，工作面宽深比为0.48，地表处于非充分采动状态。覆岩移动破坏向上传递致使关键层2中下部发生破断，关键层2初次破断距与周期破断距分别为146.6 m和70.2 m。覆岩垮落带高度增至54.8 m，再向上100 m的范围为裂缝带，即裂缝带发育高度为154.8 m，其裂采比为12.9，关键层整体表现为局部破断的挠曲变形。裂缝带向上379.1 m直达地表为弯曲下沉带范围。由于8512，8513工作面之间区段煤柱的支承作用，在覆岩一定高度范围内，破断的岩层呈现“波浪形”分布，由于关键层2在自重和上覆岩层载荷作用下产生变形，与下部岩层运动相反，因此在地表处于非充分采动状态下，关键层2与其下部岩层之间产生了离层，关键层2略有失稳。此处亦证实了在地表移动观察过程中下沉值发生突变的原因。

图9 851 3工作面开采结束Fig.9 End of 8513 working face mining

4.3.3 模拟与实测数据分析

( 1 ) 模拟与实测数据对比分析

8512 工作面开采后，实测地表最大下沉为1.124 m，数值模拟结果显示地表最大下沉为1.156 m，下沉系数为0.128，与实测实时下沉率0.124基本吻合。8513工作面开采结束后，由于关键层2的控制作用，下沉发育至地表后下沉量2.2 m，如图10所示，与现场实测的数据2.226 m基本吻合。模拟下沉量误差70 mm左右。

图10 模拟开采与实测数据对比Fig.10 Comparison of simulated mining and measured data

( 2 ) 模拟与实测数据回归分析

统计学中的假设检验和回归分析法可以研究2个变量之间是否存在相关关系及数量规律，并使用数学公式的方法来表达变量之间的关系[23]。本文应用该方法分析，模拟开采2个工作面后得到的下沉数据与该工作面在现场实测的数据，验证其是否存在相关性，并检验模拟开采数据对实测数据有无显著性影响。

对观测值与模拟值进行回归分析，回归模型见表6，并由方差分析法分析2组数据，判断系数R2值均接近最大值1，表明2变量间的线性关系密切，2组变量即为正线性相关，如图11所示；F值均小于F临界值，即在显著水平0.05下，数值模拟数据对实测下沉数据无显著性影响，说明模拟开采下沉数据与地表观测下沉数据基本吻合。也表明在数值模拟中参数的选取是合理的。

表6 模拟结果与实测结果回归模型汇总Table 6 Summary of simulation and measured regression models

图11 模拟数据与实测数据回归曲线Fig.11 Regression curves of simulated data and measured data

由前2个工作面模拟开采数据与实测值的对比分析可知，数值模拟在参数的选取以及模拟地质情况上是合理的、可取的，故而进行第3个工作面及第4个工作面开采，从而获取覆岩及地表移动变形的相应参数是必要的、科学的。

4.3.4 多工作面模拟结果分析

模拟8511工作面开采完成，如图12所示，开采规模达到360 m，开采高度为12 m，工作面宽深比为0.72，地表仍处于非充分采动状态。覆岩中关键层2完全断裂，断裂距离为70.2 m，关键层2上方有2层超厚的泥岩，会抑制导水裂缝带发育，破裂不会遍及整个岩体，而是裂缝仅发展到一定高度。覆岩垮落带的高度为54.8 m，由于厚砂质泥岩的抑制作用，且表现为同步弯曲，故砂质泥岩以上为弯曲下沉带，导水裂缝带的发育高度为196 m( 开采高度的16.3倍 )，弯曲下沉带为向上336.7 m。关键层2和上覆的泥岩相隔较大。关键层2破裂后，由于自重的原因逐渐压实，该断面煤柱的支护影响范围大于8513工作面开采结束时的支护影响范围，最大地面沉降为3.8 m。

图12 8511工作面开采结束Fig.12 End of 8511 working face mining

图13 8514工作面开采结束Fig.13 End of 8514 working face mining

模拟8514工作面开采完毕，如图13所示，开采规模达到480 m，开采高度为12 m，工作面宽深比为0.96，此时，采区已被充分开采，基本达到充分采动状态。“三带”发育高度与8511工作面开采完成后基本一致，只是下沉盆地范围增大。覆岩中关键层2与其下伏岩层之间的离层闭合，这是由于：一方面，开采尺寸的继续增大，关键层的铰接结构范围在横向扩展的同时也伴随着整体挠曲下沉；另一方面，关键层下伏岩层破断后持续压实压密到一定程度后，不再发生继续下沉，随着工作面的推进，覆岩发育转向横向发育，破坏范围扩大，下沉盆地也随着扩大。当关键层的铰接结构挠曲下沉量达到下伏岩层的最大允许下沉值时，离层不再存在。此时，在充分采动状态下，关键层2及其覆岩的压力通过岩层传递至采空区，所有这些都由该段煤柱承担，煤柱应力达到强度极限从而产生塑性变形，失去了煤柱的支承作用，完全被煤屑充填密实后，关键层3对地表起绝对控制作用，上覆岩层的运动基本达到了稳定状态，模拟结果表明最大地面沉降为5.7 m。

如图14可知，随着开采工作面的增加，除下沉值下沉曲线逐步显现下沉盆地，边界角、移动角也随着开采工作面的增加而减小，下沉盆地在扩大。下沉曲线由不对称分布逐渐变为对称分布，最大下沉点在采空区中央。

图14 模拟开采后各工作面下沉曲线Fig.14 Simulated subsidence curves of each working face after mining

综上所述：分析4个工作面开采模拟结果可知，随着开采工作面数量的增加，开采宽度增加，宽深比增大，上覆岩层也在逐步发生变化，8512工作面开采后，宽深比为0.24，属于极不充分采动状态，由于关键层的控制作用，地表沉陷小；8513工作面开采后，宽深比为0.48，地表从极不充分采动变为非充分采动，由于关键层失稳，地表急剧下沉，覆岩内部结构失稳；开采8511工作面后，宽深比为0.72，地表仍处于非充分采动状态，沉降值仍在增加，关键层完全破裂，由于厚泥岩的抑制作用，裂缝只发展到一定深度；8514工作面开采后，宽深比为0.96，关键层2的铰接结构范围已经完全密实不再沉降，覆岩的发育转向横向发展，破坏范围扩大，从而沉降盆地扩大，此时，地表基本达到充分采动状态，关键层3未产生破断，对地表起绝对控制作用，最大下沉值达到5.7 m。由8512与8513工作面对比来看，不同采高情况下，下沉量的变化较大，覆岩破坏程度较大，垮落带与导水裂缝带在采高12 m的条件下发育较大，比采高9 m条件下均有所增加。采高越大，覆岩垮落带高度增加，导水裂缝带向上发育尺寸基本不变。

4.4 对比分析

从开采至今已完成2个工作面的开采任务，地表状态也由最初的极不充分采动发育至如今的非充分采动，未来地表将达到充分采动状态。故而现将地表发育过程中模拟值与实测值进行对比说明。

由表2，7可知，地表在达到稳态后，地表下沉值达到最大，极不充分状态下由于覆岩关键层控制，地表下沉量较小；非充分采动状态下覆岩关键层破断失稳，地表下沉速度加速，地表产生较宽的台阶裂缝；充分采动状态下，覆岩由纵向发育变为横向发育，地表表现为开采沉陷区扩大，导水裂缝带范围增大。模拟得出达到充分采动时，下沉系数为0.475，最大下沉值为5.7 m；由地表移动角量参数可知，地表在发育过程中，随着下沉盆地的扩大，移动角及边界角均有约3°的减小，主要影响范围角也相应有较大幅度的减小，这表明下沉盆地范围扩大较明显，影响范围扩大；地表由极不充分采动向非充分采动发育过程中，角度变动较大，由非充分采动向充分采动发育过程中，角度变化较小，最大下沉角逐步接近90°。

表7 数值模拟角参数Table 7 Numerical simulation of angular parameters

5 结 论

( 1 ) 由地表移动观测数据分析，8512工作面开采后地表处于极不充分采动，8513工作面开采后地表处于非充分采动，地表由极不充分采动向非充分采动发育时，8512工作面采空区活化，产生二次沉降，造成下沉量突变，且下沉速度突变，即由关键层断裂失稳所致。

( 2 ) 应用关键层理论结合3DEC数值模拟方法，模拟多工作面开采关键层断裂失稳与地表移动的关系。结果表明：工作面宽深比为0.24时，地表极不充分采动状态下由于覆岩中关键层控制垮落，地表表现为沉降较小，裂缝宽度较窄，地表破坏范围较小；工作面宽深比为0.48时，地表处于非充分采动状态下覆岩关键层由离层—失稳—破断，地表也相应的表现为沉降突变、裂缝宽度增加并产生台阶裂缝且破坏范围增加；工作面宽深比为0.72时，裂采比为16.3倍采高，地表仍处于非充分采动状态，工作面宽深比为0.96时，基本达到充分采动状态下，覆岩由纵向发育变为横向发育，地表表现为开采沉陷区扩大，主要影响范围角减小，影响范围增大。

( 3 ) 模拟结果得出地表达到充分采动时，下沉系数为0.475，最大下沉值为5.7 m，大采深综放开采采高越大，垮落带高度越大，导水裂缝带向上发育尺寸基本不变，同时应用统计学回归分析与方差分析方法将地表移动观测数据与数值模拟数据做假设检验，确定了模拟参数的可靠性。为该矿实现科学开采提供了理论依据。