煤层采动地表变形特征及数值模拟研究

Numerical Simulation Research on Surface Deformation Characteristics and Numerical Simulation of Coal Seams

摘要: 为探究缓倾斜煤层开采地表变形规律及采深、倾角的影响机理，以邯郸峰峰矿区大社矿2^#、4^#煤层为研究对象，采用FLAC^3D构建三维模型，动态模拟双煤层分层开采全过程，分析不同开采阶段地表沉降及竖直位移变化规律，研究不同采深、倾角作用下的地表变形特征，并结合现场实测数据进行验证。结果表明：双煤层叠加开采地表沉降呈“椭圆形拱状分布”，最大沉降位于采空区正上方，向四周呈空间衰减趋势；随着煤层采深的增加，最大沉降值从1436 mm降至905 mm，降幅达37%，沉降影响范围减小42%；煤层倾角对沉降范围、程度、中心位置及形态影响显著，倾角每增加8˚，上山边界沉降量减少9%，下山边界增加11%。现场验证显示，FLAC^3D模型预测沉降值与实测值相对误差 ≤ 0.39%，精度可靠。研究成果为华北型缓倾斜煤层群的开采设计与生态保护提供量化依据。

Abstract: To investigate the law of surface deformation during gently inclined coal seam mining and the mechanism of coupling effects between mining depth and dip angle, this study focuses on the 2^# and 4^# coal seams of Dashe Mine in Fengfeng Mining Area, Handan. A three-dimensional model was established via FLAC^3D to dynamically simulate the full process of layered mining of the double coal seams. The evolution laws of surface subsidence and vertical displacement at different mining stages were analyzed, and the characteristics of surface deformation under the coupling effects of varying mining depths and dip angles were examined, with validation against field-measured data. The results indicate that surface subsidence induced by the superimposed mining of double coal seams exhibits an “elliptical arch-shaped distribution”, with the maximum subsidence occurring directly above the goaf and showing spatial attenuation toward the surroundings. As the mining depth of the coal seams increases, the maximum subsidence value decreases from 1436 mm to 905 mm, representing a 37% reduction, while the affected range of subsidence shrinks by 42%. The coal seam dip angle exerts a significant influence on the scope, magnitude, central position, and morphology of subsidence: for each 8˚ increase in dip angle, subsidence at the uphill boundary decreases by 9%, whereas subsidence at the downhill boundary increases by 11%. Field validation demonstrates that the relative error between the subsidence values predicted by the FLAC^3D model and the measured values is ≤ 0.39%, confirming high reliability. The research findings provide a quantitative reference for mining design and ecological protection of North China-type gently inclined coal seam groups.

1. 引言

煤炭矿产资源的开采活动不可避免地会对上部覆岩的稳定性造成影响。煤炭资源开采后，覆岩失稳引发变形破坏，促使岩体应力重分布[1]，进而引发地表沉降与移动变形，对煤矿安全生产、周边建筑及基础设施安全构成严重威胁[2]。

近年来，国内外学者针对采动地表变形开展了广泛的研究。在理论层面，基于岩土力学的经典理论，王仲仁等系统剖析了材料从弹性变形到塑性变形、断裂的力学机制，为采动岩体应力应变分析提供理论支撑[3]。然而，Akshay在研究开采过程中的覆盖层土体的塑性和弹性变形时，发现该理论在应对复杂地质条件和开采工艺时，常常难以全面、准确地描述地表变形过程[4]。随着计算技术发展，FLAC^3D等数值模拟软件被广泛应用于煤矿开采研究[5]。戴华阳利用FLAC^3D模拟煤矿地表移动变形规律并求解参数[6]；张海波则运用该软件分析铁矿采空区稳定性[7]。

目前关于煤层采动地表变形的研究多集中于单一煤层或急倾斜煤层；针对煤层倾角、开采深度等多因素耦合作用对地表沉降变化规律的影响，现有研究仍较为薄弱。陈凯等针对浅埋煤层单一采深条件下的地表沉降计算方法进行了探讨，但未涉及缓倾斜煤层中倾角与采深的耦合作用[8]；戴华阳等聚焦急倾斜煤层的岩层移动机理，而对缓倾斜双煤层叠加开采的沉降特征研究不足[9]。因此，本研究以峰峰矿区大社矿为研究对象，运用FLAC^3D软件构建三维数值模型，模拟2^#、4^#双煤层的动态开采过程，并结合现场实测数据进行验证。重点分析地表沉降的空间分布、范围及程度变化规律，揭示不同开采深度和煤层倾角对地表沉降的影响机制[10]。研究成果将为缓倾斜煤层的开采设计、矿区生态保护与安全生产提供科学依据，并为相关领域的研究提供参考[11]。

2. 研究区工作概况

研究区域位于邯郸市峰峰矿区大社矿东南部。其中，2^#煤层赋存于二叠系下统山西组，倾角8˚，埋深359.05~433.38 m，采厚5.4 m；该煤层分布稳定，结构复杂，顶板为砂质页岩，底板为泥岩，含1~2层夹矸。4^#煤层位于2^#煤层下方30~50 m，层间岩性为粉砂岩，倾角8˚，埋深385.85~470.18 m，采厚1.3 m，顶板为灰岩，底板为中细粒砂岩。

3. 研究区数值模拟

3.1. 模型建立及参数选取

为探究煤层开采诱发的地面沉降变形机制，本研究基于FLAC^3D数值模拟平台，结合研究区地质资料与地层特征，采用Mohr-Coulomb本构模型开展沉降模拟[12]，公式为：

$f_a={\sigma }_1-{\sigma }_2\frac{\left(1+\sin \phi \right)}{\left(1-\sin \phi \right)}-2c\sqrt{\frac{1+\sin \phi }{1-\sin \phi }}$ (1)

式(1)中：${\sigma }_1$ 为最大主应力；${\sigma }_2$ 为最小主应力；$c$ 为黏聚力；$\phi$ 为内摩擦角；$f_a$ > 0为判断土体是否引起了剪切损伤的临界条件。

Table 1. Rock physical and mechanical parameters

表1. 岩石物理力学参数

地层岩性	体积模量/GPa	黏聚力/MPa	剪切模量/GPa	抗拉强度/M	内角/(˚)	密度/(kg/m^3)
粉质黏土	0.28	0.02	0.01	0.02	28	1800
中粒砂岩	5.45	4.35	4.65	4.60	35	2550
砂岩	5.25	4.15	3.62	4.10	33	2580
粉砂岩1	4.60	4.82	3.60	3.80	32	2510
泥岩	4.40	2.40	2.80	2.40	30	2400
中粒砂岩	6.20	4.86	3.96	5.20	35	2650
细砂岩1	5.30	4.00	3.20	3.80	32	2580
2^#煤	1.50	1.60	1.75	0.50	29	1450
细砂岩2	5.30	4.00	3.20	3.80	32	2580
粉砂岩2	4.90	4.80	3.10	4.50	33	2610
4^#煤	1.50	1.60	1.75	0.50	29	1450
粉砂岩3	4.90	4.80	3.10	4.50	33	2610

根据地层结构，划分为12个岩组，设置岩层沿倾向、走向、竖直方向分别为x、y、z轴，构建一个800 m × 1000 m × 500 m的三维数值模型[13]，生成20,480个单元和23,595个节点，数值模型如图1所示。为消除边界效应，在煤层开采初始阶段，沿倾向和走向分别预留100 m、200 m宽的保护煤柱[14]。模型边界条件设置为：x轴与y轴方向施加位移及应力约束，顶部为自由边界。煤层采用长壁式综采全高开采工艺以及全部垮落法管理顶板，设定100 m开采步距，按“先2^#煤后4^#煤”顺序实施分层开挖，每层开挖待节点应力收敛平衡后推进，单煤层累计开挖6次，合计开挖12次。模拟过程中所涉及的岩石物理力学参数见表1。

Figure 1. Diagram of the numerical model

图1. 数值模型图

3.2. 模拟结果分析

3.2.1. 竖直位移变化分析

通过对采动过程中所有节点应力变化进行监测，得到煤层开采应力平衡后的位移变化情况。根据煤层的走向，选取了开挖200 m、400 m和600 m的模拟结果。通过研究每一步开挖所引发的塑性变形，得出竖直位移剖面图(见图2)。

从竖直位移剖面图中可以看出，地层对掘进的响应是复杂的，涉及应力释放、地层沉降及变形等多个地质力学过程。随着2^#和4^#煤层开采的深入，岩层的竖直位移持续发生着变化，开采活动对地表的影响呈现出明显的空间分布特征，地表沉降呈“椭圆形拱状分布”，位移值在掘进区域上方形成高值区，且随着开挖深度的增加，高值区的范围和强度均有所增大；位移的影响范围随着掘进深度的增加而扩大，在600 m深度时，地表对开采活动的响应较为敏感，位移梯度更为显著，掘进引起的地表位移影响范围显著大于在200 m深度时。

(a) 2^#煤开挖200 m竖直位移剖面图 (b) 2^#煤开挖400 m竖直位移剖面图

(c) 2^#煤开挖600 m竖直位移剖面图 (d) 4^#煤开挖200 m竖直位移剖面图

(e) 4^#煤开挖400 m竖直位移剖面图 (f) 4^#煤开挖600 m竖直位移剖面图

Figure 2. Vertical displacement profile of 2^# and 4^# coal seams mining

图2. 2^#、4^#煤层开采竖直位移剖面图

Figure 3. Coal seam mining settlement curve

图3. 煤层开采沉降曲线图

由图3可知，当2^#、4^#煤层开挖完成，地层应力释放和岩体结构调整达到相对稳定状态后，地面沉降的最大值为5708.2 mm，最大位移产生在采空区的顶板处，向地表依次减小。地面沉降表现为从沉降中心向两侧逐渐减小，且随着距离沉降中心的增加，沉降值逐步降低。

3.2.2. 垂直应力变化分析

为探究采煤过程中围岩垂直应力的分布与演化规律，通过数值模拟获取了2^#与4^#煤层分别推进至200 m、400 m和600 m时的垂直应力剖面云图。由图4可知，煤层采动使周围岩体应力发生不同程度的变化，应力状态总体呈“马鞍状”分布，沿走向表现出较好的对称性。工作面上方一定范围内覆岩垂直应力降低，且随着工作面开挖距离的增大，应力减少区范围也增大。当2^#煤层推进到200 m时，初步形成应力等值线拱形态，在煤层采空区正上方竖向应力较低，而在开切眼及煤壁附近出现应力集中现象，最大垂直应力可达24.86 MPa。随着2^#开挖至600 m时，采空区中部因前期垮落岩块逐渐压实，形成局部应力升高区，关键层作为主要承载结构，控制顶板持续弯曲下沉，其下方软弱岩层与已垮落岩体共同构成协同支撑体系[15]；当4^#煤层推进到600 m时，覆岩中的应力等值线拱形态趋于稳定，采动影响范围进入相对稳定的阶段；此时关键层承担主要上部荷载，采空区内垮落岩体充填裂隙、承担部分岩体压力，在一定程度上缓解了覆岩的移动变形。

(a) 2^#开挖200 m时垂直应力剖面图 (b) 2^#开挖400 m时垂直应力剖面图

(c) 2^#开挖600 m时垂直应力剖面图 (d) 4^#开挖200 m时垂直应力剖面图

(e) 4^#开挖400 m时垂直应力剖面图 (f) 4^#开挖600 m时垂直应力剖面图

Figure 4. Profile diagram of vertical stress variation during coal seam mining

图4. 煤层开采垂直应力变化剖面图

4. 现场实测数据验证

为评估FLAC^3D数值模型的预测精度，对模拟预测的地表沉降值进行了现场实测验证[16]。选取了研究区域内4个沉降观测点(1027、1028、1029、1030)，记录并计算出观测点的实际沉降值。通过将FLAC^3D模型预测的沉降值与实测值进行对比，并计算相对误差，量化了预测值与实测值之间的差异。FLAC^3D数值模拟沉降值以及观测点实际沉降值统计结果如表2所示。

Table 2. Actual settlement values and predicted settlement values

表2. 实际沉降值与预测沉降值

沉降观测点	监测周期(月)	实际沉降值(mm)	FLAC^3D预测沉降值(mm)	相对误差(%)
1027	6	4970.000	4975.131	0.10
1028	6	5050.000	5066.315	0.32
1029	6	5120.000	5140.213	0.39
1030	6	4876.000	4890.512	0.29

由表2可知，FLAC^3D模型对1027、1028等监测点沉降值的预测结果与实测值的相对误差仅为0.10%~0.39%，体现了较高的预测精度。同时，长期监测显示，开采结束后6个月地面沉降趋于稳定，这与模型模拟的应力收敛时间吻合，进一步验证了模型的可靠性。FLAC^3D数值模拟对区域经验依赖性较小，但对地层岩体物理力学参数的准确性及合理性要求较高，使其能够为工程实践提供可靠的预测结果[17]。

5. 不同影响因素下地表沉降变化规律

矿产资源的开采正朝着更深采深和更大倾角的方向发展[18]。以峰峰矿区大社矿为例，当前煤层开采工作面主要为8˚~25˚的缓倾斜煤层，但部分煤矿煤层工作面已达到25˚~45˚的倾斜煤层[19]。在这一背景下，研究不同采深和倾角条件下地表沉降的变化规律，对于确保矿山开采的安全性、稳定性以及环境保护具有重要的理论价值和现实意义[20]。

5.1. 不同因素模拟方案设计

为了深入探讨不同采深和倾角条件下煤层开采对地表沉降的影响特征。研究选择煤层上山方向采煤边界中点(P)、开采区域沉降最大位置(O)以及下山方向采煤边界中点(Q)作为研究对象(见图5)，通过模拟不同采深和倾角条件下的地表沉降变化，研究各因素对地表沉降的影响。

Figure 5. Schematic diagram of the position of points P and Q and observation points

图5. P、Q点和观测点位置示意图

具体模拟方案设计如下：

方案一：在煤层开采倾角为8˚时，分别模拟不同采深(200 m、250 m、300 m、400 m)条件下的地表沉降情况，得到不同采深条件下的地表沉降云图(见图6)；分析不同采深条件下，上山方向采煤边界中点(P)、开采区域沉降最大位置(O)和下山方向采煤边界中点(Q)的沉降值随煤层采深变化的趋势(见图7)。

方案二：在煤层开采深度为300 m时，分别模拟不同倾角(0˚、16˚、24˚、32˚)条件下的地表沉降情况，得到不同煤层倾角下的开采区域沉降分布图(见图8)；分析不同倾角条件下，上山方向采煤边界中点(P)、开采区域沉降最大位置(O)和下山方向采煤边界中点(Q)的沉降值随煤层倾角变化的趋势(见图9)。

5.2. 不同采深煤层开采地表沉降分布及变化趋势

(a) 煤层采深为200 m时地表沉降云图 (b) 煤层采深为250 m时地表沉降云图

(c) 煤层采深为300 m时地表沉降云图 (d) 煤层采深为400 m时地表沉降云图

Figure 6. Cloud map of surface subsidence in mining area

图6. 开采区域地表沉降云图

结合图6和图7，采深对地表沉降的影响呈现“范围缩小、峰值衰减”的协同效应。具体表现为：随着煤层采深的增加，地表沉降程度呈现规律性递减，这是由于深部开采覆岩厚度增大，岩体自重压力增强，抑制了变形向上传递。地面沉降呈现从中心向两侧逐步衰减的特征，颜色梯度的变化直观地反映了沉降程度的减小。沉降的影响范围随着采深的增加而扩大，在200 m~400 m采深范围内，最大沉降值从1436 mm降至905 mm，降幅达37%；因深部覆岩约束更强，变形向周边扩散的范围受限，沉降影响范围减小42%；沉降区域从开采中心向外扩展，开采活动对地表稳定性的影响逐渐减弱。

Figure 7. Variation of subsidence values in the mining area

图7. 开采区域沉降值变化曲线图

由图7可知，煤层采深是影响地表沉降的关键因素，在不同煤层采深条件下的地表沉降变化规律如下：

1) 在采深变化初期，上山方向采煤边界中点(P)的沉降值变化较为平缓，但随着采深超过300 m后，沉降值变化速率显著加快，整体呈现凸曲线形态。

2) 开采区域沉降最大位置(O)的沉降值与煤层采深呈负相关关系，整体变化为凸曲线形态。随着采深的增加，沉降值变化速率增大，但沉降值逐渐减小。

3) 下山方向采煤边界中点(Q)沉降值在初期变化平缓，随采深增大变化速率逐渐加快，整体呈凸曲线形态，且始终大于上山方向采煤边界中点(P)沉降值；随着采深的增大，P、O、Q三点的沉降值逐渐趋于一致。

5.3. 不同倾角煤层开采地表沉降分布及变化趋势

图8模拟结果表明，煤层倾角是影响地表沉降范围、程度、中心位置和形态的关键因素之一，通过改变覆岩重心偏移方向影响沉降分布。当倾角为0˚时(图8(a))，沉降中心对称分布于采空区正上方，而随着倾角增加到32˚ (图8(d))，地表沉降区域呈现出从集中到分散的变化趋势，沉降中心向下山方向偏移28 m，因下山方向覆岩在重力分力(F = Gsinθ)作用下更易发生剪切破坏，裂隙发育范围比上山方向大35%；地表沉降的程度和形状也随着倾角的变化而变化，倾角为0˚时沉降程度较大，且沉降区域形状呈现较为规则的圆形，而随着倾角的增大，沉降程度有所减小，且沉降区域形状变得不规则。

(a) 煤层倾角为0˚时地表沉降云图 (b) 煤层倾角为16˚时地表沉降云图

(c) 煤层倾角为24˚时地表沉降云图 (d) 煤层倾角为32˚时地表沉降云图

Figure 8. Cloud map of surface subsidence in mining area

图8. 开采区域地表沉降云图

Figure 9. Variation of subsidence values in the mining area

图9. 开采区域沉降值变化曲线图

由图9可知，地表沉降的变化趋势与煤层倾角之间存在显著的相关性。随着煤层倾角的增加，不同位置的沉降值呈现出差异化变化：

1) 倾角每增加8˚，上山边界(P点)沉降量减少9%，下山边界(Q点)沉降量增加11%。

2) 关键层承载重心下移，与煤柱协同支撑抑制变形，上山方向采煤边界中点沉降值逐渐减小，两者呈负相关关系；地表最大沉降值逐渐减小，两者呈负相关关系；下山方向采煤边界中点沉降值呈正相关关系，关键层受重力分力(F = Gsinθ)易剪切破断，难成稳定承载结构，塑性变形累积使沉降增加。

6. 结论

1) 缓倾斜双煤层叠加开采引发的地表沉降呈现显著的“椭圆形拱状分布”特征，沉降最大值出现在采空区正上方，且向四周呈空间衰减趋势。随着开采深度的增加，地表位移高值区的范围与强度均呈扩大态势，凸显了采动影响的空间延展性。

2) 采深对地表沉降存在“范围缩小、峰值衰减”的协同效应。在200~400 m采深范围内，最大沉降值从1436 mm降至905 mm，降幅达37%，沉降影响范围缩减42%；地表沉降最大位置的沉降值与采深呈负相关，且随采深增大，沉降值变化速率加快并显著减小；上、下山方向采煤边界中点沉降值变化呈现出先平缓后加速的趋势，整体表现为凸曲线。

3) 煤层倾角对地表沉降的范围、程度和形态有显著影响。倾角增大时，沉降区域从集中变为分散，沉降中心向下山方向偏移，形态趋于不规则；沉降值在上山方向采煤边界呈负相关，下山方向则呈正相关；倾角每增加8˚，上山边界沉降量减少9%，下山边界沉降量增加11%。

4) FLAC^3D模型在地表沉降预测中表现出高精度与可靠性，与现场实测值的相对误差 ≤ 0.39%，沉降稳定时间与模拟收敛时间吻合，为工程实践提供了可靠的理论支持。

NOTES

^*通讯作者。