1. 引言
煤层底板突水是我国煤矿开采过程中常见且严重的地质灾害之一，对矿工生命安全和矿井正常生产构成重大威胁[1] [2]。随着煤矿开采深度的增加，底板突水[3]-[5]的风险不断上升，因此，准确预测煤层底板的破坏深度和突水潜力成为确保矿井安全的重要课题。传统的预测方法包括经验公式法和现场监测法，这些方法虽然简单实用，但在应对复杂地质条件和大规模采掘活动时，常常显得不足。
近年来，数值模拟方法因其能够考虑复杂的地质条件和应力变化，成为研究煤层底板破坏行为的有效工具。FLAC3D (快速拉格朗日差分法三维模拟)作为一种先进的数值模拟软件，广泛应用于岩土工程和矿井安全领域。通过FLAC3D [6]，可以模拟煤层底板在开采过程中应力场和位移场的变化，从而预测底板的破坏深度[7] [8]，为制定有效的突水防治措施提供科学依据。
本研究以朱仙庄煤矿为例，应用FLAC3D模拟煤层底板的采动破坏过程，结合经验公式法，对底板破坏深度进行预测和分析[9] [10]，以期为煤矿的安全生产和防突措施提供参考。
2. 经验公式法预测底板导水破坏深度
2.1. 经验公式法
参照《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》中相关公式计算朱仙庄矿煤层底板采动后破坏深度值。仅使用工作面斜长做参数，底板采动破坏深度h为：
H = 0.7007 + 0.1079L
式中：L为壁式工作面斜长，m。
考虑采场深度、倾角以及工作面斜长时，则有：
H = 0.0086H + 0.1666α + 0.1079L − 4.3679
式中：H为煤层深度，m；α为煤层倾角，(˚)。
经验公式计算以III1031工作面参数为依据，H取918 m，α为15˚，L为190 m。计算得1031工作面理论破坏深度为26.4 m。
2.2. 首采面情况
III1031工作面位于矿井南部III3采区10煤层一区段。东南至朱仙庄矿井田边界，西北至III3采区回风上山98 m左右，上邻II1059工作面(设计面)，下邻III1033综采工作面(设计面)。工作面回采上限为−760 m，下限为−918 m，走向长1210 m，倾斜宽为190 m，煤层倾角为8˚~26˚，煤层厚度0.75~3.53 m，平均厚度为2.6 m。
3. FLAC3D数值模拟法预测导水破坏深度
FLAC3D是快速拉格朗日差分分析(Fast Lagrangian Analysis of Continua)的简称，是一种分析流体运动的方法，它是由美国ITASCA咨询集团公司开发的三维快速拉格朗日分析程序。三维快速拉格朗日法基于三维显示有限差分法，广泛应用于模拟岩土或其他材料的三维力学特性。三维快速拉格朗日分析采用了显示有限差分格式来求解场的控制微分方程，并且应用了混合单元离散模型。在计算时将计算区域划分为若干六面体单元，单元网格可以随着材料的变形而变形，而且每个单元在给定的边界条件下遵循指定的线性或非线性本构关系。FLAC采用离散元的动态松驰法求解，应用了节点位移连续的条件，利用拖带坐标系分析大变形问题，并可用差分格式按时步积分求解，不但可以对连续介质进行大变形分析，而且还能模拟岩土体沿某一软弱面产生的滑动变形。FLAC能较好的模拟地质材料在达到强度极限或屈服极限时发生的破坏或塑性流动的力学行为，分析渐进破坏和失稳，特别适用于模拟大变形。
3.1. 模型的建立和本构关系的选取
根据地面钻孔柱状图资料，朱仙庄矿10煤底板到奥陶纪灰岩有多种岩性，其中主要以砂岩、泥岩和灰岩为主。当研究区域的空间尺寸过大，地层结构横向与纵向跨度大且岩性复杂多样，且研究区域的地层包含着许多岩性不同的岩石层，如果将研究区域的所有地层都建立模型，这将带来巨大的工作量和计算量，这在实际工程应用中不太实用，所以在实际工程应用中，有必要将地层结构中物理力学性质相近的岩层进行复合，简化研究地层结构。根据朱仙庄矿钻孔综合柱状图，研究地层结构化为16层，即从顶部第1层到底部16层。
1. 边界条件
底部边界采用全约束边界条件，即X、Y和Z向位移均为零；左右边界、前后边界为约束边界条件，即X向、Y向均固定水平位移为零，Z向为自由边界；模型上边界为自由边界，不予约束。
2. 本构关系的选取
煤系底板岩体本构关系被选为理想弹塑性本构模型，采用Mohr-Coulomb塑性本构模型和Mohr-Coulomb屈服准则。
3. 计算模型的建立
(1) 模拟模型地质情况
数值模拟以朱仙庄矿采区1031工作面为工程背景，10煤平均采厚2.4 m，倾角约15˚。模拟工作面标高−918 m，工作面走向长1210 m，倾斜长190 m，顶板采用自由垮落方式管理。
(2) 模型参数
以1031工作面为参考，模型力学参数(见表1)，考虑边界作用效应的影响和网格的剖分，模型沿走向方向左侧预留100 m，走向方向右侧预留90 m，即模型沿走向长1400 m。沿倾向方向两侧各预留55 m，即模型沿倾向方向长300 m。其中x轴为走向方向，y轴为倾向方向，z轴为垂直方向。最终模型的尺寸为：1400 m × 300 m × 280 m。
Table 1. Model mechanical parameters
表1. 模型力学参数
岩性 |
体积模量(GPa) |
切变模量(GPa) |
容重(kg/m3) |
内聚力(MPa) |
内摩擦角 |
抗拉强度(MPa) |
厚度(m) |
砂岩1 |
17.23 |
16.76 |
2.66 |
6.36 |
40.6 |
3.3 |
10 |
泥岩1 |
0.86 |
0.69 |
2.66 |
0.24 |
29 |
2.3 |
16 |
砂岩2 |
17.23 |
16.76 |
2.66 |
6.36 |
40.6 |
2.6 |
10.8 |
10煤 |
0.6 |
0.18 |
1.33 |
0.9 |
29 |
1.6 |
2.4 |
泥岩2 |
0.86 |
0.69 |
2.66 |
0.24 |
29 |
2.6 |
7.2 |
砂岩3 |
17.23 |
16.76 |
2.66 |
6.36 |
40.6 |
12.2 |
22.2 |
泥岩3 |
0.86 |
0.69 |
2.66 |
0.24 |
29 |
2.6 |
41.6 |
灰岩1 |
18.32 |
11.22 |
2.766 |
6.2 |
40 |
6.26 |
47.2 |
泥岩4 |
0.86 |
0.69 |
2.66 |
0.24 |
29 |
2.62 |
16.7 |
灰岩2 |
18.33 |
11.22 |
2.766 |
6.2 |
40 |
6.97 |
10.7 |
泥岩6 |
0.86 |
0.69 |
2.66 |
0.24 |
29 |
2.61 |
8.1 |
灰岩3 |
18.33 |
11.22 |
2.766 |
6.2 |
40 |
6.98 |
27 |
泥岩6 |
0.86 |
0.69 |
2.66 |
0.24 |
29 |
1.74 |
26.2 |
灰岩4 |
18.33 |
11.22 |
2.766 |
6.2 |
40 |
6.98 |
10.2 |
泥岩7 |
18.33 |
11.22 |
2.766 |
6.2 |
40 |
6.06 |
8.4 |
砂岩4 |
17.23 |
16.76 |
2.66 |
6.36 |
40.6 |
13.6 |
6 |
(3) 地质模型(见图1)
Figure 1. Numerical model diagram of 1031 working face
图1. 1031工作面数值模型图
3.2. 数值模拟结果分析
在初始应力场建立完以后，模型顶部施加基岩竖向荷载，对煤层开挖进行模拟，首先对自重应力下产生的速度和位移进行清零操作。设计开挖方式为分段开挖，开挖步距20 m，总共开挖10次，共开挖200 m。然后对开挖结果进行分析，主要是分析由于开挖引起的竖向位移变化，垂直应力场变化及塑性区破坏特征。同时根据开挖产生的煤层顶板塑性区范围和应力分布两个方面综合分析确定底板破坏带发育高度，并与底板破坏带的预计公式进行对比，以此检验本次数值模拟的有效性。
1. 位移场变化规律
随着采煤工作面模型的不断推进，煤层竖直方向基岩不断发生变形破坏，从而使煤层底板围岩向采空区方向位移。为了研究煤层下方各岩层的竖向位移特征及其发生变形破坏时的位移发展演化过程，所以需要对开采过程进行数值模拟。由于工作面走向采长较大，因此本次模拟选取开挖步距为20 m，主要分析推进距离分别为20 m、40 m、60 m、80 m、100 m、120 m、140 m、160 m、180 m、200 m时的覆岩竖向位移变化规律，通过数值模拟计算得出不同推进距离下的竖向位移如图2所示。
(a) 推进20 m (b) 推进40 m
(c) 推进60 m (d) 推进80 m
(e) 推进100 m (f) 推进120 m
(g) 推进140 m (h) 推进160 m
(i) 推进180 m (j) 推进200 m
Figure 2. Vertical displacement diagram of different advancing distance of working face
图2. 工作面不同推进距离竖向位移图
由图2可知，当煤层开挖后，本身存在的煤层区域被挖掉，致使该区域原有初始地应力平衡被破坏，导致该10煤层的底板出现了垂直应力，该区域在底板的垂直应力作用下必将引起位移的产生，所以在10煤层的底板出现了垂直向上的竖直位移，煤层底板的表现特征为直接底板向上隆起。工作面推进至80 m时，由于开挖卸荷，底板发生底鼓产生向上的位移量，最大值为44.6 mm。随着工作面不断推进，最后工作面推进至180 m结束后，煤层底板向上底鼓位移达到最大值71.6 mm。综上各开挖阶段位移值变化规律，可知底板竖向位移值随着工作面推进距离的增大而增大，底板竖向位移变化呈现出拱形的特征。
2. 应力场变化规律
天然状态下，基岩–煤层各岩体处于平衡状态，当煤层进行开采后形成采空区，这种平衡状态将会被破坏，从而使围岩应力重新分布。采动后当岩石强度低于其所能承受的极限压力时，就会发生岩层变形破坏，其应力分布规律是随着工作面不断推进而发生变化的，因此要分析煤层开采条件下的应力场演化规律，同样要分析工作面推进20 m、40 m、60 m、80 m、100 m、120 m、140 m、160 m、180 m、200 m时竖向垂直应力场的分布规律，其应力场分布数值模拟计算结果如图3所示。
(a) 推进20 m (b) 推进40 m
(c) 推进60 m (d) 推进80 m
(e) 推进100 m (f) 推进120 m
(g) 推进140 m (h) 推进160 m
(i) 推进180 m (j) 推进200 m
Figure 3. The distribution law of vertical stress field at different advancing distances of working face
图3. 工作面不同推进距离垂直应力场分布规律图
由图3可知，随着1031工作面沿走向方向不断的推进，10煤层原始存在的平衡应力被破坏，导致煤层底板出现了最大拉应力的状态。从图中可以得到应力值为开挖进尺推进到80 m时的最大拉应力为2.29 MPa，开挖进尺推进到200 m时的最大拉应力为2.27 MPa。
4. 塑性区变化规律
为了研究各破坏区的发育演化规律，主要分析工作面推进20 m、40 m、60 m、80 m、100 m、120 m、140 m、160 m、180 m、200 m时的塑性区变化规律，其塑性区分布规律数值模拟计算结果见图4。
煤层开采后原始地应力被破坏，岩层为了达到一种新的应力平衡，在平衡的过程中，随着开挖进尺的不断推进，研究煤层底板底鼓塑性区分布特征，可直观分析得出底板破坏带发育高度。开采煤层后，由于底板所受破坏程度不同，所以煤层底板下塑性破坏区破坏程度不同。如图4，当工作面推进至40 m时开始出现塑性破坏，推进180 m时塑性破坏达到最大值，最大值为26.2 m，这时岩层破坏都为拉伸破坏，继续推进至180 m时，破坏深度并未向深部发展。至此可推测，继续沿走向面开挖的情况下，底板破坏深度发育的最大值为26.2 m。
运用经验公式法预测底板破坏深度值为26.4 m，运用数值模拟的手段获得研究区内10煤在采动作用下，底板最大破坏深度约为26.2 m，两种预测结果相近，有待实践验证。本次研究取两者较大值，即研究区内10煤层开采条件下底板破坏深度为26.4 m。
(a) 推进20 m (b) 推进40 m
(c) 推进60 m (d) 推进80 m
(e) 推进100 m (f) 推进120 m
(g) 推进140 m (h) 推进160 m
(i) 推进180 m (j) 推进200 m
Figure 4. The distribution law of plastic zone in different advancing distances of working face
图4. 工作面不同推进距离塑性区分布规律图
5. 总结
(1) 随着开采活动进行，应力主要集中在工作面两端，煤层工作面顶部以及底部均有应力拱产生。
(2) 综合分析位移云图和塑性图，预测得出三水平首采面底板破坏深度值为26.2 m，这与经验公式法所预测的破坏值26.4 m相近。本次研究取两者较大值，即三水平1031工作面内10煤层开采条件下底板破坏深度为26.4 m。