近水平煤层开采地表移动边界的预测研究

2023-10-07孙学阳袁官思龙苗霖田杜荣军

煤矿安全 2023年9期

孙学阳，袁官思龙，苗霖田，杜荣军，程正

（1.西安科技大学地质与环境学院，陕西西安 710054；2.陕西省煤炭绿色开发地质保障重点实验室，陕西西安 710054；3.煤炭资源勘查与综合利用重点实验室，陕西西安 710054）

获取采煤沉陷特征值是研究煤矿开采地表移动变形规律和揭示其形成机理的基础，现行规范和规程[1-2]根据采煤沉陷特征值对开采沉陷盆地边界和损毁程度做了界定。理论分析方、相似材料模拟实验、数值模拟是目前地表移动和变形预测的3 种方法。理论分析受到相关变量的限制，相似材料模拟实验受到实验条件的限制，20 世纪60年代R W Clough 首次将有限元引入土石坝的稳定性分析中，从此数值模拟技术成为了岩土工程领域主流的研究和设计方法[3]。因此对于采煤沉陷特征值的获取和地表移动边界的预测研究中，数值模拟具有优势。

地表损毁边界的主要影响因子为下沉、倾斜、水平变形。胡振琪等[4]以理论分析结合现场调查，提出了以耕地损毁临界下沉值、临界水平变形值和临界倾斜值中的最小值圈定耕地损毁范围；邓伟男[5]利用FISH 语言代码，提取模型主断面上的下沉、水平移动、倾斜、水平变形及曲率值来绘制曲线图，为数值模拟中数据的处理提供了参考；文献[6-8]采用FLAC3D软件进行模拟开采时，通过在走向和倾向主断面上布设监测点的方式获得其采煤沉陷特征值；LI 等[9]从应力场和位移场的角度出发，结合地表影响范围的放大比，分析了地表沉降影响区随采空区不断扩大而发生的动态变化；WANG 等[10]研究了在FLAC3D建立含复杂地层边界的三维数值计算模型的方法，并通过MATLAB 软件导出地表沉降和水平移动数据，进而由计算公式转化出倾斜、水平变形和曲率的等值线图来直观反映地表移动变形的影响范围；蔡音飞等[11]利用Mathematica 软件对FLAC3D模拟结果进行再次分析得到地表移动变形量，并通过列表绘图函数ListPlot 绘制各次试验的变形图对地表受损情况进行评价。

基于此，从沉陷盆地全面积的角度，运用FISH 语言提取计算采煤沉陷特征值并进行运用，圈定了霍宝干河煤矿多工作面开采后地表移动盆地最外边界、地表危险移动边界、地表移动裂缝边界；研究结果不仅为开采沉陷数值模拟数据后处理提供了新思路，还为生态恢复治理规划提供参考。

1 地表移动变形数据的提取

地下岩土体的移动变形是采煤沉陷问题的核心,是地面沉陷的原因[12]。沉陷就是移动变形的总体描述，地表沉陷状态按移动变形性质、方向、方式的不同进行描述，可分为垂直方向的移动变形和水平方向的移动变形。工程中通常采用下沉、倾斜、曲率，水平移动、水平变形等5 个采煤沉陷特征值进行定量评价[13]。无论是实测数据分析还是形变预计，采煤沉陷特征值都是基于下沉和水平位移2 个移动分量[14]。因此先将水平移动分量和垂直分量提取出来，再进行相应计算就能得到倾斜、水平变形和曲率。

在FLAC3D软件中，力作用在单元中心坐标上，位移产生在网格节点上。针对这一特征，利用FISH 语言对六面体网格主导的计算模型中各个网格单元的特定节点进行遍历，进行数据的提取并保存为Surfer 能够识别的数据格式。

除此之外，FLAC3D中提供了丰富的内置函数，用户可以运用它们进行简单的程序设计，并得到相应结果。针对采煤沉陷问题，有必要先对FLAC3D中模型网格的节点编号和单元编号进行一下梳理，同时也是提取数据的代码的编写依据。

为了便于阐述，只建立了3× 3× 1 的模型，利用内置函数gp.pos(gp)、data.label.create(v)、data.label.text(lp)可以将单元编号和每个单元的节点编号在FLAC3D中显示出来。FLAC3D中单元和节点编号示意图如图1，图中：红色数字为每个单元的编号；蓝色数字为1 号单元的节点编号。规律如下：单元编号在走向上（x方向）依次排列，在倾向上（y方向）按等差数列排列；节点编号在单个单元上遵循右手螺旋法则。

图1 FLAC3D 中单元和节点编号示意图Fig.1 Schematic diagrams of unit and node number in FLAC3D

在FLAC3D中，对特定范围的指定有2 种方式，自定义分组和自动分面；遍历的对象也有2 种，一种是对单元进行遍历，另外一种是对节点进行遍历。

对于下沉值和水平移动值仅对节点进行遍历即可；倾斜值和水平变形值先对单元进行遍历，然后在具体编号的单元内通过节点指针找到相应编号节点的坐标和位移，曲率需要同时在2 个单元内通过节点指针找到相应编号节点的坐标和位移。最后，通过FISH 语言编写出倾斜、水平变形和曲率的函数来提取数据并进行文本文件的保存，为绘制等值线提供数据基础。

1.1 下沉值数据提取及保存

采用自动分面命令，模型表面将自动被命名为“Top i”，这就是下沉值的遍历范围。垂直移动分量为负，由下沉值正负号规定可知：正值表示测点下沉，负值表示测点上升，所以要在变量前面添加负号。另外，FLAC3D的默认单位为m，下沉值的单位需要转换成mm，具体代码如下：

1.2 水平移动值数据提取及保存

水平移动值遍历范围同上，它的正负号规定：指向矿层走向断面走向方向的移动为正，反之为负；指向倾斜断面上山方向的移动为正，下山方向的为负。水平移动值涉及倾向和走向两个方向的分量，在FISH 变量中分别用dx 和dy 表示。单位为mm，具体代码如下：

1.3 倾斜值数据提取及保存

倾斜值的计算涉及同一单元的相邻两节点，根据每个单元的节点编号规则，对提取倾斜值的代码设计如下：走向上倾斜值为每个单元的8 号节点和6 号节点（或7 号节点和4 号节点）的下沉值之差与它们之间沿x轴水平距离的比值；倾向上倾斜值为8 号节点和7 号节点（或6 号节点和4号节点）的下沉值之差与它们之间沿y轴水平距离的比值[15]。

倾斜值反映了地表移动盆地沿某一方向的坡度，通常以i 表示，下列代码中ix 代表走向上的倾斜值，iy 代表倾向上的倾斜值。它的正负号规定：指向矿层走向断面走向方向的倾斜为正，反之为负；指向倾斜断面上山方向的倾斜为正，下山方向的为负。单位需转换为mm/m，具体代码如下（以走向为例）：

1.4 水平变形值数据提取及保存

水平变形简言之就是两点间单位长度上的水平移动差，用相邻两点的水平移动差与两点间水平距离的比值来表示。如图1（b），走向上水平变形同样用8 号和6 号节点（或7 号和4 号节点）的x轴位移分量之差与它们之间沿x轴水平距离的比值表示；倾向上水平变形用6 号和4 号节点（或8号和7 号节点）的y轴位移分量之差与它们之间沿y轴方向水平距离的比值表示。竖直方向的变化发生倾斜，水平方向的变化发生水平变形。所以参考提取倾斜值的代码，将内置函数gp.disp.z()在走向上改成gp.disp.x()，在倾向上改成gp.disp.y()，节点编号和坐标函数对应修改即可。它的正负号规定：单位长度线段的拉伸变形为正，压缩变形为负[16]，单位为mm/m。

1.5 曲率值数据提取及保存

曲率反映了地表垂直断面上观测线的弯曲程度，它可以用相邻线段的倾斜差与两线段中点间的水平距离之比来表示。

如图1，走向上的曲率为相邻单元8 号节点和6 号节点连成线段的倾斜值之差与两线段x轴中点间水平距离的比值；倾向上的曲率为相邻单元的6 号节点和4 号节点连成线段的倾斜值之差与两线段y轴中点间水平距离的比值。它的正负号规定：上凸表示地表下沉曲线为正，下凹则为负。单位需转换成mm/m2，具体代码如下（以倾向为例）：

其中a、b、c、d为网格单元编号，e为沿走向方向的网格单元个数，g=a-e。

2 近水平煤层开采沉陷模拟

霍宝干河煤矿为近水平煤层，第四系黄土平均厚100 m，属于厚松散层，煤层平均厚度为4.5m，开采深度为460～500 m。

模型尺寸为2 000 m ×2 050 m ×600 m，共分为8 层，平均采深为480 m，采厚4.5 m，煤层倾角为0°，共计116 850 个单元，132 396 个节点；数值计算模型分为3 个工作面进行开采模拟，分别是1101 首采工作面，1102 工作面和1103 工作面。3 个工作面尺寸均为600 m ×210 m，工作面之间均留设10 m 宽的护巷煤柱。模型剖切图如图2，岩层物理力学参数见表1。

表1 岩层物理力学参数Table 1 Physical and mechanical parameters of rock strata

由于地表移动观测站只设置在1101 首采工作面，为了保证3 个工作面全采后数值模拟结果的正确性以及用FISH 语言提取特征值的合理性，有必要先将该工作面下沉盆地主断面的特征值进行提取，1101 工作面主断面特征值如图3～图5。

图4 1101 工作面主断面倾斜值和水平变形值Fig.4 Inclination values and horizontal deformation values of main section of 1101 working face

图5 1101 工作面主断面曲率值Fig.5 Curvature values of main section of 1101 working face

将图3～图5 中各特征值曲线的极值进行统计，1101 工作面地表移动变形量极值见表2。

表2 1101 工作面地表移动变形量极值Table 2 Extreme values of surface movement and deformation in 1101 working face

从表2 中可以看出用FISH 语言提取的各特征值与现场实测资料相吻合，验证了该数值模拟的正确性及用FISH 语言提取地表移动变形数据的合理性。

3 应用实例

3.1 Surfer 绘图流程

Surfer 软件的主要优势是利用其插值功能进行等值线图的绘制。用它来绘制等值线图的前提是数据文件格式要转换成“.grd”文件格式[17]。

由于Surfer 的插值功能，即使数据不是等间距的，它依然能够作图。如果不进行数据的插值，只是绘制原始数据的图，可以在“.text”转“.grd”的对话框中选择距离平方反比法或克里金插值法，本实例所有等值线图的绘制过程皆采用克里金插值法。

运行提取采煤沉陷特征值的各段FISH 代码会得到相应数据文件，按上述步骤操作之后新建等值线图即可，按照项目的要求进行加密、着色等设置。除此之外，Surfer 中还可以添加边界文件，可用于工作面轮廓线的绘制，为测量地表各边界与采空区边界线间的距离提供了便利。

3.2 地表移动边界成图

地表必然受到煤炭开采的影响，提前预判这一影响的1 个重要环节就是对变形后的地表进行等值线的绘制[18]。实测法获取地表移动盆地角值参数存在困难，建立地表移动观测站成本高，观测数据匮乏；类比法获取角值参数又受到地质采矿条件的制约；统计模型法通过统计经验公式获取角值参数的方法又受到实测数据和地质采矿条件的双重制约。

综上所述，随着1101 工作面、1102 工作面和1103 工作面的开采，将整个模型通过FISH 代码提取采煤沉陷数据并绘制等值线图，从数值模拟的角度预测近水平煤层多工作面开采后地表移动边界。地表下沉等值线图如图6，地表移动边界等值线图如图7，图中：蓝色线框代表的是工作面轮廓，红色和紫色等值线为临界值等值线。

图6 地表下沉等值线图Fig.6 Contour map of surface subsidence

图7 地表移动边界等值线图Fig.7 Contour diagrams of surface movement boundary

3.2.1 地表移动盆地最外边界

地表移动盆地最外边界以下沉10 mm 确定，该边界在走向方向上距离采空区边界518 m，在倾向方向上距离采空区边界502 m。该边界以内的区域都属于开采沉陷影响区，这个边界以外的区域一般认为不受开采沉陷的影响。

3.2.2 地表危险移动边界

地表危险移动边界按照现行规范[1]中对砖混结构建筑物临界变形值的规定：倾斜变形值i=±3 mm/m，水平变形值ε=±2 mm/m，曲率K=±0.2 mm/m2来确定。

图7（a）和图7（b）中的倾斜临界值等值线为红色椭圆线，图7（c）和图7（d）走向上曲率最大值为0.045 6，倾向上曲率最大值为0.036 4 均未达到曲率临界值，因此未能绘制出曲率的临界等值线，图7（e）和图7（f）中紫色等值线为水平变形临界值等值线；由此可以看出走向上地表危险移动边界与采空区边界左侧相距299 m，右侧相距291 m，倾向上的地表危险移动边界与采空区边界上侧相距308 m，下侧相距316 m。上下左右间距不相等的原因是由于工作面推进方向和工作面开采顺序不同造成的。地表破坏的程度和范围在等值线图上的直观反映就是倾斜和曲率大的区域[19]。

3.2.3 地表移动裂缝边界

开采沉陷引起的土体拉伸破坏临界值与土体变形值的定量关系如式（1）[20]：

式中： εxc为土体产生裂缝的水平变形临界值，mm/m；h为土体埋深，m； ρ为平均密度，t/m3；φ为摩擦角，（°）； μ为泊松比；E为弹性模量，MPa；c为黏聚力，MPa。

εxc随土体的埋深增加而增大，土体中任意一点的水平拉伸变形值εx≥εxc时，该处土体产生拉伸破坏，否则仅产生连续的沉陷变形。地表（h=0）处土体破坏极限状态的水平变形临界值εh=0为最小，可通过式（2）计算。

将黄土的物理力学参数E=10 MPa、μ=0.3、c=0.02 MPa、φ=15°代入式(2)得到地表产生拉伸变形的水平变形临界值为εh=0=2.15 mm/m。

将该临界值的等值线绘制于图7（e）和图7（f）中，以红色等值线标示出，可以看出，走向上左侧裂缝边界距离采空区边界498 m，右侧裂缝边界距离采空区边界484 m；倾向上上侧裂缝边界距离采空区边界476 m，下侧裂缝边界距离采空区边界493 m。同样，上下左右间距不相等的原因是由于工作面推进方向和工作面开采顺序不同造成的。红色等值线不仅确定了裂缝边界，而且圈定了产生裂缝区域的范围，如图7（e）和图7（f）。