地表移动影响范围与地质采矿条件关系数值模拟

2022-11-04张琰君贺福帅阎跃观朱元昊

科学技术与工程 2022年27期

张琰君，贺福帅，阎跃观*，朱元昊

(1.中国矿业大学(北京)地球科学与测绘工程学院，北京 100083； 2.中国电建集团中南勘测设计研究院有限公司，长沙 410000)

煤炭开采造成了地表不同程度的破坏，严重影响了地面建构筑物的安全。FLAC3D是基于拉格朗日算法的有限元差分软件，在模拟稳定性、振动、较大的非线性形变问题时具有独特优势，目前FLAC3D已广泛应用于注浆加固及充填开采[1-2]、边坡稳定[3]、采空区稳定性分析[4-5]、开采沉陷[6-7]等多个方面并且日趋成熟。

在开采沉陷方面，许多学者已经通过数值模拟的手段探究了各种地质开采条件下的覆岩及地表移动变形规律。刘剑等[8]运用FLAC3D探究了余吾煤业在大采宽开采条件下的沉陷规律，与实测值对比得到了较高的模拟准确性。余学义等[9]以FLAC3D数值模拟为技术手段探究了厚松散层大采高的地表移动变形规律，得出了地表下沉量与覆岩岩性成正比，下沉速度与覆岩岩性成反比的结论。张天军等[10]利用FLAC3D模拟了增加开采工作面及增加开采深度对对变电站及其周围的影响并得到了一系列成果。庞会等[11]以谢桥矿为研究背景，通过构建FLAC3D煤层开采模型得到了地表移动影响范围，总结了不同推进距离下覆岩运移规律。但是学者们仅针对单一特殊地质采矿因素研究地表移动范围，缺乏多种因素综合考虑研究地表移动影响范围。

研究认为松散层厚度、基岩厚度、覆岩岩性、工作面开采尺寸等是影响地表移动变形规律的重要因素[9-12]。鉴于此，现通过建立FLAC3D三维数值模型，运用控制变量法探究地表移动影响范围随岩层厚度、松散层厚度、工作面采长及岩层物理力学性质改变的不同变化趋势，得到地表移动影响范围L(10 mm下沉边界)与地质采矿条件的相关关系方程，为矿区的三下采煤及建构筑物保护提供参考。

1 FLAC3D实例模拟

1.1 观测站概述

A工作面开采4 #煤层，走向长度1 000 m，倾向宽200 m，平均开采深度500 m，黄土覆盖厚度20 m，煤层平均厚度2.3 m，煤层倾角2°，属于近水平煤层。采煤方式为走向长壁综采放顶煤开采，顶板管理方法为全部垮落法及强制放顶。观测站布设简图见图1。

图1 观测站布设图Fig.1 Layout diagram of observation station

1.2 FLAC3D模型建立

以A工作面实际地层构造为背景建立三维数值模型，模型尺寸为3 000 m×1 000 m×515.3 m，开挖尺寸为1 000 m×200 m×2.3 m。模型采用Mohr-Coulomb计算准则，边界采用位移约束方式，x方向边界采用x方向约束，y方向边界采用y方向约束，底部约束3个方向的位移，顶面为自由面，计算初始地应力场按岩体自重应力场考虑，表1为模型各层的物理力学参数，FLAC3D模型及布点情况见图2。

1.3 模拟精度分析

从FLAC3D中提取地表移动的下沉值，汇总后在表2列出了部分下沉值较大点位的模拟精度情况(地表下沉值是数值模型平衡后所提取的结果，因此可将其等同于实际煤层开采稳定后的沉降值)。根据地表岩移资料，与稳定后的实测下沉值对比其模拟最大绝对误差不超过50 mm，最大相对误差不超过5%，绝对误差与相对误差均较小，说明FLAC3D数值模拟结果较为准确。

图2 A工作面数值模型Fig.2 Numerical model of working face A

表1 模型岩层物理力学参数Table 1 Physical and mechanical parameters of working face A

表2 模拟与实测下沉值对比Table 2 Comparison of simulated and measured subsidence

2 L与地质采矿条件关系模拟

在上述A工作面模拟精度较好的情况下，为避免结果出现偶然性，因此同时选取了不同地质条件下的另一矿区开采工作面共同探究地表移动影响范围L与地质采矿条件关系，图3为地表移动影响范围示意图。

同实验选取的另一对象为B工作面，B工作面开采3 #煤层，走向长度1 000 m，倾向宽200 m，平均开采深度313.5 m，黄土覆盖厚度15.5 m，煤层平均厚度4.4 m，水平煤层。采煤方式为走向长壁综采放顶煤开采，顶板管理方法为全部垮落法。根据B工作面所处煤岩层地质条件及布点情况构建数值模拟计算模型，见图4。模型尺寸为3 000 m×1 000 m×327.9 m，开挖尺寸为1 000 m×200 m×4.4 m，模型仍采用Mohr-Coulomb计算准则，与实测地表移动变形量对比结果相近的最佳岩石力学参数如表3所示。

图3 地表移动影响范围Fig.3 Influence range of surface movement

图4 B工作面数值模型Fig.4 Numerical model of working face B

表3 B工作面模型物理力学参数Table 3 Physical and mechanical parameters of working face B

2.1 L与岩层厚度关系

在上述模拟效果较好的条件下，采用控制变量法以采长1 000 m、松散层情况及覆岩岩性保持不变为前提，改变岩层层厚探究地表移动影响范围与其关系，实验分别改变A工作面泥岩厚度和B工作面粉砂岩1厚度，厚度变化范围0～120 m，从FLAC3D中提取各次计算结果并绘制地表移动影响范围变化量ΔL与岩层厚度变化量Δm关系见图5。

从图5中实验结果可获得，ΔL与Δm近似呈线性关系，地表移动影响范围随着岩层厚度的增加而增加，属正相关关系，抽象后可拟合得A工作面ΔL与Δm的线性回归方程为：ΔL=0.68Δm，B工作面ΔL与Δm的线性回归方程为：ΔL=1.23Δm。

图5 地表移动影响范围变化量与岩层厚度变化量关系Fig.5 Relationship between variation of influence range of the surface movement and variation of bedrock thickness

2.2 松散层厚度对L的影响

在采长为1 000 m、岩层厚度和覆岩岩性保持不变的前提下，改变松散层厚度探究其对地表移动影响范围的影响，A工作面松散层厚度从0 m增加至120 m，B工作面松散层厚度从15.5 m增加到124.5 m，绘制地表移动影响范围变化量ΔL与松散层厚度变化量Δh关系见图6。

图6 地表移动影响范围变化量与松散层厚度变化量关系Fig.6 Relationship between variation of influence range of surface movement and variation of loose layer thickness

从图6中可以看出，ΔL与Δh整体近似呈线性关系，地表移动影响范围随着松散层厚度的增加而增加，属正相关关系，抽象后可拟合得A工作面ΔL与Δh的线性回归方程为：ΔL=0.71Δh，但地表移动影响范围在松散层为120 m时发生了突增，B工作面ΔL与Δh的线性回归方程为：ΔL=0.42Δh。

2.3 L与采长关系

分别改变A工作面和B工作面开挖长度(依据启动距为0.25～0.5倍采深，确定A、B工作面起始长度分别为300 m和200 m)研究地表移动影响范围与工作面采长的关系，其余岩层、松散层厚度及覆岩岩性均不变，A工作面每次推进距离100 m，B工作面每次推进距离200 m，统计各次开挖结果并绘制地表移动影响范围变化量ΔL与采长变化量Δs关系见图7。

从图7可以看出，两个工作面表现出一致的规律性，地表移动影响范围与采长呈负相关关系，ΔL与Δs整体呈线性关系，抽象后可拟合得A工作面ΔL与Δs的线性回归方程为：ΔL=-0.11Δs，B工作面ΔL与Δs的线性回归方程为：ΔL=-0.05Δs。

2.4 L与覆岩岩性关系

在采长1 000 m的基础上保持岩层厚度与松散层厚度不变，改变岩层的物理力学参数探究地表移动影响范围与覆岩岩性的关系，实验为了具有代表性选择改变关键层参数，将原来的较坚硬的关键层分别替换为砂质泥岩和泥岩，实验中分别改变了A工作面的中粗粒砂岩层和B工作面的粗砂岩1岩层，改变及替换的岩层物理参数见表4。

获得A和B工作面对应岩层下的地表移动影响范围见表5，通过表5可知将A工作面中粗粒砂岩替换为砂质泥岩和泥岩后，地表移动影响范围分别增加了4%和6.6%，而B工作面替换岩层后地表移动影响范围分别增加了2.5%和3.4%，地表移动影响范围受岩石力学性质影响轻微。

图7 地表移动影响范围变化量与采长变化量关系Fig.7 Relationship between variation of influence range of the surface movement and variation of mining length

表4 改变岩层物理力学性质Table 4 Changing physical and mechanical properties of rock

表5 改变岩层后地表移动影响范围变化Table 5 Changes of influence range of surface movement

2.5 L与地质采矿条件关系

由上述实验可知，地表移动影响范围L与岩层厚度、松散层厚度、采长呈线性关系，同时又会受到岩层力学性质的影响，因此可以得到L与地质采矿条件关系式为

L=am+bh+cs+d

(1)

式(1)中：a、b、c为岩层厚度、松散层厚度、采长斜率参数，与地质采矿条件有关；d为与岩层力学性质相关的改正数。

在2.1～2.3节已使用控制变量法得出a、b、c系数，将2.4节中模拟结果代入式(1)可得，A工作面L与地质采矿条件关系式为：L=0.68m+0.71h-0.11s+d，结合地质采矿参数，其中当关键层为中粗粒砂岩时，d=239.4 m；当替换为砂质泥岩后，d=258.4 m；当替换为泥岩后，d=270.4 m。B工作面L与地质采矿条件关系式为：L=1.23m+0.42h-0.05s+d，其中当关键层为粗砂岩时，d=31 m；当替换为砂质泥岩后，d=40 m；当替换为泥岩后，d=43 m，可以看出d随着岩层变软弱而逐渐增大。

3 L与地质采矿条件关系讨论

3.1 L与岩层厚度关系分析

由上文可知不同地区的两工作面模拟结果相似，证实利用FLAC对地表移动变形的研究是可信且正确的，因此下文各小节中只选择其中一处工作面展开讨论。在采长1 000 m、覆岩岩性及松散层不变的情况下，以B工作面为例，改变粉砂岩厚度将得到的下沉值绘制为下沉曲线见图8。

从图8中可以看出，随着粉砂岩厚度的增加，地表最大下沉值在逐渐减小但地表影响范围略有增大，曲线基本形态相似且逐渐变平缓。基岩厚度变大提高了岩层的稳定性使得地表下沉值变小，同时岩层变厚使得采深更大，根据地表影响范围与采深呈正相关关系，采深越大地表影响范围也越大，这与岩层变厚地表移动影响范围变大的实验结果是一致的。

图8 粉砂岩厚度变化时的下沉曲线Fig.8 Subsidence curve of siltstone thickness variation

3.2 松散层厚度对L的影响分析

在采长1 000 m、岩层厚度和覆岩岩性不变的前提下，以A工作面为例改变松散层厚度得到各次下沉结果绘制曲线如图9所示。

图9 松散层厚度变化时的下沉曲线Fig.9 Subsidence curve of loose layer thickness variation

从图9中可以看出松散层0～100 m下沉曲线形态变化较小，在120 m时下沉曲线产生明显变化，表现为下沉值增大，地表影响范围显著增大，这可能是由于巨厚松散层导致荷载过大破坏了基岩的稳定性，使得地表再次发生沉陷。地表移动影响范围随松散层厚度增大而增大符合厚松散层地质条件下影响范围广的特性，实验进一步证实了已知的地表移动变形规律。

3.3 L与采长关系分析

以B工作面为例，将得到的各次地表移动变形值提取统计后，绘制下沉曲线如图10所示。

从图10中可以看出，随着工作面不断推进最大下沉点也在向前移动，下沉值不断增大，最终地表形成一个碗状盆地。当开挖分别为200、400 m和600 m时下沉曲线拐点附近切线与水平线夹角对应为31°、47°和50°，推进过程中拐点附近切线与水平线夹角不断增大且角度增长量逐渐减小，下沉曲线越来越陡峭，这与地表移动影响范围不断减小且逐渐稳定是吻合的。

图10 不同采长下的下沉曲线Fig.10 Subsidence curves under different mining lengths

分别将开挖200、400、600 m的结果绘制应力云图见图11。通过应力云图可以发现开采完毕后岩体稳定形成了压力拱，随着开挖长度的不断增加，拱顶部分不断延伸，而竖向均布载荷作用下拱的合理轴线是二次抛物线，因此实际压力拱的稳定性在降低，这就导致开采长度越大地表移动盆地活动越剧烈，下沉曲线越陡峭，地表移动影响范围反而越小，合理解释了地表移动影响范围随采长变化的原理。

图11 不同开挖长度下的应力云图Fig.11 Stress nephogram under different excavation lengths

3.4 L与覆岩岩性关系分析

在控制采长为1 000 m和岩层、松散层厚度不变条件下，替换关键层后将A工作面和B工作面各次实验结果绘制下沉曲线如图12所示。

图12 改变岩层力学参数下沉曲线对比Fig.12 Comparison of subsidence curves of changing rock mechanical parameters

从图12(a)下沉曲线中可以看出，将A工作面较坚硬的中粗粒砂岩关键层分别替换为较软弱的砂质泥岩和泥岩后下沉值明显增大，最大下沉值接近煤层开挖厚度(2 300 mm)，表明没有关键层支撑将会导致岩层全部垮落并压实，将会对地面造成极大的影响，相应的参数d也较大(250 m左右)。从图12(b)下沉曲线中可以看出将B工作面较坚硬的关键层替换为砂质泥岩和泥岩后，下沉值明显增大但下沉最大值未达到采厚(4 400 mm)，原因是除粗砂岩外还有粉砂岩关键层存在，这种非单一关键层的支撑效应对地表移动破坏起到了很好的保护作用，参数d较小(40 m左右)。