磁窑沟煤矿黄土沟壑地貌开采沉陷规律研究

2022-07-05杨玉龙宋选民

山西煤炭 2022年2期

杨玉龙，宋选民

(1.山西河曲晋神磁窑沟煤业有限公司, 山西忻州 036500；2.太原理工大学原位改性采矿教育部重点实验室，太原 030024)

目前，我国能源结构中煤炭仍占据着主导地位，绝大部分煤炭资源来自于井工煤矿，其工作面的回采在不同程度上引起了地表沉陷变形[1]，使得煤炭资源开发与矿区环境保护之间的矛盾日益凸显[2-3]。河东煤田是山西省重要产煤基地之一，煤炭资源赋存较浅、煤层厚度较大、上覆黄土沟壑错综复杂，工作面回采后对矿区地表产生较大影响，制约了矿区经济社会的可持续发展[4-5]。

针对黄土沟壑地貌的开采沉陷问题，国内外众多专家学者取得了大量有益研究成果。ALAEE et al[6]基于模糊理论，建立了有关地表开采沉陷的模糊隶属度函数，并于伊朗Tabas煤矿进行了现场验证。KOPEC et al[7]对波兰上西西里亚煤田开采引发的地表沉陷进行了实测研究，分析了其对地表道路、房屋等建筑物的影响。魏国武[8]以大佛寺煤矿30201工作面为背景，总结得出矿山开采沉陷影响因子，构建了矿山开采沉陷危险性系数计算模型。吴文敏等[9]研究了黄土高原丘陵地貌对地表移动变形的影响，得出受采动影响沟谷两侧向下坡方向滑移，沟谷底部的地表会发生抬升，滑移分量比值与坡度存在线性递减的关系。赵兵朝等[10]研究了生态脆弱区厚松散层薄基岩、黄土沟壑及山区丘陵下煤层开采的地表损害特征，并以基采比和基载比为关键参数，揭示了生态脆弱区煤层开采损害特征形成的机理。陈超等[11]根据采动地裂缝的力学成因进行分类，系统梳理了不同类型采动地裂缝形成机理研究进展，分析了采动地裂缝的影响因素与研究方法，讨论了相关的关键科学问题。张磊等[12]运用数值模拟方法，研究了黄土覆盖区域煤层开采引起的地表沉陷特征，认为受采动及上覆岩层自重的影响，地表沉陷变形随采掘范围的增加而加大。贺国伟等[13]基于支持向量机，构建了适用于黄土矿区的最大下沉预计模型，提高了黄土地表最大下沉值预计精度和可靠性。李春意等[14]采用现场实测、理论分析和离散元数值模拟相结合的方法，分析了地表沉陷的基本特征，构建了基于坐标变换的倾向地表下沉和水平移动拟合函数模型，确定了地表沉陷的角量参数和预计参数。温建生[15]以柳湾矿61115工作面地表沉陷实测资料为基础，分析了地表移动变形动静态规律，得出地表下沉集中在活跃期并且下沉速度较大，基岩移动角偏大而综合边界角偏小。

综上所述，对于煤层开采影响下的黄土沟壑地表沉陷预测的研究需要根据矿井实际地质生产条件综合分析判断。本文以磁窑沟煤矿13101工作面生产地质条件为工程背景，运用理论分析及数值模拟相结合的方法，分析其开采引起的地表移动变形规律，对该矿厚煤层综放开采条件下的地表沉陷预判具有一定的指导意义。

1 工程背景

1.1 13101工作面工程地质条件

13101工作面位于磁窑沟矿二水平一盘区，主采13#煤层，位于井田南部，上部为10#煤层采空区，平均间距为20 m；东部为13102工作面(未形成)，以东为13#煤层未开采的实体煤；南部为山西华鹿阳坡泉煤矿，留设20 m保护煤柱，该煤矿目前未开采13#煤层;西部为13#煤层未开采的实体煤，与留设的寺也村村庄保护煤柱相接；北部为131盘区辅助运输大巷/较大大巷/回风大巷，以北为13#煤一盘区，为未开采的实体煤。工作面宽度135.6 m，长度为1 579.4 m，煤层总厚9.40～11.10 m，平均10.25 m，煤层倾角2.8°～4.6°，平均3.7°，采用综合机械化放顶煤工艺进行回采。

1.2 13101工作面地表典型特征

13101工作面地面标高+1 026～+1 135 m，工作面标高+889～+921 m，盖山厚度122～215 m;松散层厚度40.4～108.0 m，平均74.2 m;基岩厚度65.5～119.0 m，平均95.0 m。工作面地表地形为山西西北部黄土高原中-低山区地形，地面植被稀少，冲刷剧烈，沟谷成“V”字形羽毛状分布，地表大部分为新生界黄土及红土覆盖，且陡峭悬崖有较小规模的黄土崩塌现象。地表地形总体北高南低，山坡和沟谷分布广泛，大部分为旱地和荒地。距切眼1 500～1 580 m地表为沟谷，大气降水汇集后形成的积水会在短时间内流入北部的磁窑沟-洞沟一带沟壑，其地表形貌如图1所示。

图1 13101工作面地表

2 概率积分法预计参数及计算分析

2.1 相关参数选取

13101工作面地表变形理论分析中，首先需要根据工作面实际生产、地质条件，依据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》的相关规定，对下沉系数、水平移动系数、主要影响角正切等参数进行确定。

1)下沉系数。由于13101工作面上覆10#煤层均已采空，且两层煤平均间距为20 m，故根据以下公式计算：

qf=(1+a)qc.

(1)

式中:qf为复采下沉系数；a为下沉活化系数，取0.2；qc为初采下沉系数，取0.55。代入磁窑沟矿相关数据可知，13101工作面上覆岩层下沉系数为0.66。

2)水平移动系数b。重复采动条件下，水平移动系数与初次采动相同，即：bf=bc，故取bf为0.25。

3)主要影响范围角正切tanβf。重复采动时tanβf较初次采动增加0.3～0.8。对于中硬岩层可按式(2)计算：

tanβf=tanβc+0.06236lnh-0.017.

(2)

式中：tanβf为重采时主要影响范围角正切；tanβc为初采时主要影响范围角正切，取2.16；h为第二层煤的采深，取平均埋深164.4 m。带入磁窑沟矿相关数据可知，13101工作面上覆岩层主要影响范围角正切tanβf为2.48。

2.2 地表变形计算结果

根据以上分析所得相关系数，对13101工作面回采影响下地表变形进行理论预测，包括地表最大下沉值、最大水平移动值、最大倾斜变形值、最大曲率变形值及最大水平变形值，各参数的具体计算过程如下。

1)地表最大下沉值计算。

充分采动

Wcm=dqcosα.

(3)

非充分采动

Wbm=dqncosα.

(4)

式中:Wcm为充分采动条件下地表最大下沉值，mm；Wbm为非充分采动条件下地表最大下沉值，mm；q为充分采动条件下的下沉系数，取0.66；d为煤层法向开采厚度，取值10.25 m；α为煤层倾角，(°)；n为地表充分采动系数,其中

(5)

式中：n1和n3大于1时取1；k1、k3为与覆岩岩性有关的系数，坚硬岩层：k1、k3=0.7；中硬岩层：k1、k3=0.8；软弱岩层k1、k3=0.9；l1、l3为倾向及走向工作面长度，m；h0为工作面平均开采深度，m。13101工作面宽度135.6 m，长度为1 579.4 m，平均埋深164.4 m，工作面上覆岩层以砂岩为主，为中硬岩层，代入上式计算可知13101工作面地表最大下沉值为5.48 m。

2)地表最大水平移动值计算。沿煤层走向方向上的最大水平移动

Um=bWmax.

(6)

式中，Um为最大水平移动值，mm；b为水平移动系数，取0.25；Wmax为最大下沉值，mm。代入13101工作面地表最大下沉值5.48 m计算可知，沿煤层走向的最大水平移动为1.37 m。

沿煤层倾斜方向的最大水平移动

Um=(b+0.7P)Wmax，

(7)

式中:P为覆岩综合评价系数；h为表土层厚度，m；α为煤层倾角，(°)。磁窑沟13101工作面，上覆松散层平均厚度为69.4 m，平均埋深164.4 m，煤层平均倾角3.7°，地表最大下沉值为5.48 m。代入上式计算可知13101工作面地表沿煤层倾向方向上的最大水平移动为1.37 m。

3)最大倾斜变形值计算

(8)

4)最大曲率变形值计算

(9)

式中:kmax为最大曲率变形，mm-1。磁窑沟13101工作面地表最大下沉值为5.48 m，带入上式计算可知，最大曲率变形1.89 mm-1。

5)最大水平变形值计算

(10)

式中:εmax为最大水平变形，mm/m。磁窑沟13101工作面地表最大下沉值为5.48 m，带入上式计算可知，最大水平变形值0.47 mm/m。

3 黄土沟壑地表三维数值模型

3.1 数值模型及边界条件

借助Midas-GTS/ANSYS～FLAC3D耦合建模方法,可以实现复杂地质模型的构建，模型可将各岩层、地质构造按照等高线图真实地表现出来，大大提高了模拟预计分析的准确度和可靠性。选取13101工作面中部作为建模区域，以该区域内勘探钻孔为岩层赋存依据，由于煤层开采所引起的地表下沉范围大于工作面长度(136 m)，因此数值模型沿工作面方向分别向外扩展100 m，故数值模型在工作面方向的长度为430 m，工作面推进方向上模型长度为500 m，模型垂直方向按照区内实际钻孔勘测结果进行建模，共包含岩层30层，计算单元327 495个，计算节点342 380个，如图2所示。

图2 三维数值模型

基于上述模型，对13101工作面进行回采，根据实际情况设置边界条件，模型包含5个位移边界，即模型四周平面固定对应法向变形，模型底面固定3个主方向的变形。

3.2 数值模拟本构及参数

由于本次数值模拟中所涉及材料均为岩土类，因此选用FLAC3D软件中内嵌摩尔-库伦本构，计算过程中需用户赋值的材料参数包括：容重、抗拉强度、内聚力、内摩擦角、弹性模量及泊松比，由区内钻孔数据可知，13101工作面上覆岩层大致可分为9种岩性，各岩性及参数汇总见表1所示。

表1 上覆岩层岩性及参数

3.3 数值模拟方案

首先根据前述参数对模型材料进行赋参，计算初始平衡，因区内钻孔勘测结果中，13#煤层上部存在10-1#煤层采空区，且目前该采空区已压实稳定，因此将10-1#煤层所对应模型设置为采空区模拟参数，故在本章数值模拟过程中不需对10-1#煤层采空区进行开采模拟。初始平衡后对13#煤层进行分步开挖，初次开挖50 m，以模拟工作面的前期开挖工况;其后，以20 m为开挖步距对剩余工作面进行回采模拟，直至工作面回采至450 m，总计开挖模拟21次。从工作面回采过程中的地表下沉量、水平移动量、剪切应变及剪应变增量4个方面出发，分别提取工作面推进不同长度时的地表各项参数模拟结果，对地表变形机理展开分析。

4 数值模拟结果及分析

数值模拟共开挖21次，由于篇幅所致，以工作面推进50 m、150 m、250 m、350 m及450 m为例，提取各参数模拟结果，对13101综放工作面回采引起的黄土沟壑地表沉陷变形进行分析。

4.1 地表下沉量

工作面推进过程中，黄土沟壑地表下沉量分布云图如图3所示。由图3可知，随着工作面推进至50 m、150 m、250 m、350 m以及450 m时，工作面上方地表最大下沉量为1.50 m、5.08 m、6.12 m、6.41 m以及6.55 m，可见随着工作面的推进，地表下沉量不断增加，但下沉增速由0～50 m的30.02 mm/m，下降至350～450 m的1.38 mm/m，可见采空区覆岩运移趋于稳定。在任一推进距离，采区地表垂直方向沉陷出现“中部大两侧小”的现象，这主要是由于采空区两侧为实体煤，对两侧覆岩运移具有限制作用，导致采空区中部覆岩以竖向直接下沉为主，两侧覆岩向中部倾斜。空区地表西部黄土坡沿垂直方向下沉大于东侧，推测是由于西侧山坡较为陡峭，可能出现黄土陡坡失稳滑落，影响地表构筑物安全。黄土“山脊”位移较黄土“山谷”位移大，则表明黄土沟谷地表的防护重点应在边坡山脊处。

(a)推进50 m

4.2 地表水平移动量

工作面的回采会引起上覆地表产生水平移动，将其沿工作面长度和推进两个方向进行分解，分别提取分布云图，分析其水平移动变形特征。

1)工作面长度方向。工作面推进过程中，黄土沟壑地表沿工作面长度方向水平移动分布云图如图4所示。

由图4可知，随着工作面推进至50 m、150 m、250 m、350 m以及450 m时，工作面上方地表东侧沿工作面长度方向最大水平位移值为0.67 m、2.37 m、2.92 m、3.00 m以及3.09 m。随着工作面的推进，地表水平位移值不断增加，但变形增速由0～50 m的13.28 mm/m，下降至350～450 m的0.78 mm/m，可见东侧覆岩运移趋于稳定。工作面上方地表西侧沿工作面长度方向最大水平位移值为0.86 m、3.12 m、3.84 m、3.92 m以及3.99 m，随着工作面的推进，地表水平位移值不断增加，但变形增速由0～50 m的17.26 mm/m，下降至350～450 m的0.71 mm/m，可见西侧覆岩运移趋于稳定。同时，在工作面长度方向上，采区地表沉陷出现“中部小两侧大”的现象，这主要是由于采空区两侧为实体煤，对两侧覆岩运移具有限制作用，导致采空区中部覆岩以竖向直接下沉为主，两侧覆岩向中部倾斜，出现沿着工作面长度方向的水平位移。采空区地表西部黄土坡沿工作面长度方向水平位移大于东侧，推测是由于西侧山坡较为陡峭，可能出现黄土陡坡失稳滑落，影响地表构筑物安全。黄土“山脊”位移较黄土“山谷”位移大，则表明黄土沟谷地表的防护重点应在边坡山脊处。

(a)推进50 m

2)工作面推进方向。工作面推进过程中，黄土沟壑地表沿工作面推进方向水平移动分布云图如图5所示。

(a)推进50 m

由图5可知，当工作面推进至50 m、150 m以及250 m时，工作面上方地表沿工作面推进方向最大水平位移值为0.71 m、2.54 m以及2.31 m，可见随着工作面的推进，地表水平位移值不断增加并趋于稳定；同时工作面上方覆岩向采空区方向运移，导致地表出现超前于工作面推进距离的变形。由于西侧黄土山较高，向采空区后方偏移量较大，已出现黄土陡坡滑塌危险。当工作面推进至350 m与450 m时，采空区后侧上方地表出现向采空区方向偏移趋势，偏移量分别为0.17 m及0.26 m，同时西侧陡坡位移量增加，东侧缓坡亦出现少量位移。

4.3 地表最大剪切应变

工作面推进过程中，黄土沟壑地表最大剪应变分布云图如图6所示。由图6可知，随着工作面推进至50 m、150 m、250 m、350 m以及450 m时，工作面上方地表剪切应变范围不断增加，说明地表失稳、破坏范围逐渐增加。如图6所示，采空区地表中部剪切应变较小，四周剪切应变量大，说明采空区四周围岩均出现向采空区中部倾斜，造成地表出现盆地状变形现象；同时采空区上覆地表形貌复杂，对剪切应变分布产生一定的影响，集中体现为地表山坡顶端剪切应变加大，坡中与坡底剪切应变较小，说明黄土坡整体稳定性较差，在工作面推进时以及推进后，应当对黄土坡采取相应治理措施，防止黄土坡整体翻覆影响地表结构。

(a)推进50 m

4.4 地表最大剪应变增量

工作面推进过程中，黄土沟壑地表最大剪应变增量分布云图如图7所示。由图7可知，工作面推进50 m、150 m、250 m、350 m以及450 m时，采空区上方地表最大剪应变增量范围不断增加，且最大剪切应变增量集中分布于采空区地表中部沟底位置，说明采空区中部沟底应变量最大，最容易发生破坏，当工作面推进时应当注重对工作面中部沟底区域地表的维护与及时处理。同时可见，两侧黄土坡剪切应变增量集中于坡顶，表明随着工作面推进，地表两侧黄土坡易出现整体倾覆，在工作面推进时以及推进后应当注重对黄土坡进行保护，防止出现黄土坡整体滑塌。

(a)推进50 m

综合概率积分计算结果及数值模拟结果，由于13101工作面地表为黄土沟壑地貌，因此在概率积分计算时选取了煤层平均埋深作为计算参数，所获得的值可以认为是地表移动变形的平均值。而在数值模拟部分，通过真实地表形貌特征建立了三维数值模型，获取了13101地表移动变形的具体数据。本文概率积分法计算结果中，13101工作面地表平均沉陷值为5.48 m，平均水平移动值为1.37 m，在回采450 m数值模型中，于最大下沉值位置取平行于工作面的剖面，计算得到该剖面地表各节点下沉的平均值为5.59 m，同样的方法可以计算得到水平移动平均值的数值模拟结果为1.64 m。因此，概率积分计算结果与数值模拟结果相互验证，保证了研究结果的可靠性。