1. 引言
矿山压力(简称“矿压”)显现是煤炭地下开采过程中的核心力学现象，其分布与演化规律直接关系到工作面稳定性、瓦斯涌出控制以及地表沉陷预测[1] [2]。随着采煤工作面的推进，上覆岩层原有的应力平衡状态被打破，应力重新分布，并在采空区周围形成支承压力带，具体可划分为卸压区、应力集中区和原始应力区。其中，卸压区的形成是由于煤壁边缘在集中应力作用下被压酥，产生大量裂隙，导致其承载能力显著下降；而应力集中区则因煤体处于双向或三向应力状态，其强度增大，承受着远高于原岩应力的峰值压力。这一过程的显著标志是顶板的周期性破断，即初次来压与周期来压。来压步距是表征顶板稳定性和矿压强度的关键参数。
更为重要的是，卸压区内煤体的大范围塑性破坏与裂隙发育，会使其透气性相较于原始煤层显著提升(据研究可提高数十至数千倍)，这不仅直接影响瓦斯抽采效率，也为瓦斯在采空区的运移和积聚创造了条件，构成了安全生产的重大隐患[3]。同时，采动引起的应力重分布会随着时间的推移向上传递，最终引发地表移动与变形，形成采煤沉陷区。因此，精确揭示采动过程中煤岩体的应力分布、塑性区范围及卸压区特征，是实现矿井灾害精准防治的理论基础[4]。
在岩土工程与采矿科学领域，数值模拟已成为揭示复杂地质力学过程不可或缺的工具。其中，基于显式有限差分法的FLAC3D (Fast Lagrangian Analysis of Continua in 3 Dimensions)软件，因其在处理大变形、非线性材料及复杂边界条件方面的独特优势，在国际学术界与工程界享有盛誉。
FLAC3D内置了丰富的本构模型(如Mohr-Coulomb、应变软化/硬化、Hoek-Brown等)，能够精确模拟岩土体从弹性到塑性破坏的全过程。其“命令驱动”模式为构建复杂三维模型和实现多工况参数化分析提供了极大灵活性[5]。在采矿工程应用中，FLAC3D被广泛用于模拟巷道围岩稳定性、工作面矿压规律、底板突水机制以及地表沉陷预测等。尽管FLAC3D在矿压模拟方面取得了丰硕成果，但多数研究侧重于应力场本身的分布规律，或将模型参数简化为经验值，未能充分结合研究矿井特定煤岩层的实测力学参数，从而限制了模型预测的精准性。
基于前述研究背景及FLAC3D数值模拟方法的适用性，本文的核心在于将现场实测的煤岩力学参数与数值模拟技术相结合，对山西省朔州市平鲁区后安煤炭有限公司后安煤矿(以下简称“后安煤矿”) 110203工作面的矿压分布规律进行定量化、可视化的深入研究。本文旨在达成四个具体目标：
(1) 构建基于实测力学参数的高精度三维数值模型
以实验室试验获取的11#煤层及其直接顶板(页岩)的真实力学参数为核心输入，建立能够真实反映井下地质条件的FLAC3D三维数值模型。该模型将充分考虑覆岩多层结构，确保模型在几何形态与材料属性上的合理性，为后续的精确模拟奠定坚实基础。
(2) 揭示不同来压步距下的矿压分布与塑性区演化规律
通过模拟110203工作面在顶板初次来压(推进20 m)与周期来压(推进30 m)两个关键阶段的采动过程，系统分析煤壁前方支承压力的分布特征。重点获取垂直应力云图与屈服塑性区分布图，从而精确确定卸压区宽度、应力集中系数及峰值应力位置，揭示不同来压步距对围岩应力场与破坏范围的影响机制。
(3) 定量分析采动应力对地表沉陷及煤层透气性的潜在影响
基于模拟得到的应力重新分布与塑性区范围，分析采动影响向上传递的规律，为地表沉陷边界的圈定提供力学依据。同时，结合卸压区煤体破碎、裂隙发育的特征，定量探讨其对煤层透气性的提升效应，评估其对工作面瓦斯涌出与抽采治理的潜在影响，为瓦斯灾害防治提供理论参考。
(4) 提出针对性的工程建议并验证模型的实用性
将数值模拟的主要结论与现场观测现象或已有经验数据进行对比分析，验证模型的可靠性。最终，将研究成果转化为对矿井安全生产具有直接指导意义的对策建议，具体包括：优化工作面支护参数、确定合理的瓦斯抽采巷位置、划定沉陷区域等，为实现矿区的安全高效开采与绿色可持续发展提供科学支撑。
2. 模型建立
数值模型的合理性是确保模拟结果准确可靠的前提。为精确模拟工作面的采动影响，同时兼顾计算效率，模型尺寸的确定需遵循以下原则：其范围应足以包含主要的应力影响区，并有效消除边界效应。根据该矿地质条件最终确定模型的几何尺寸为：长(X方向) × 宽(Y方向) × 高(Z方向) = 206 m × 191 m × 53.8 m。该模型在高度上完整涵盖了11#煤层、其直接顶板(页岩)、上覆的粉砂岩、砂岩、中粒砂岩等关键岩层，直至地表等效载荷作用面。
采用FLAC3D内置的六面体网格进行划分。为提高计算精度，对110203工作面采空区、开挖巷道及煤壁前方应力集中区域进行了局部网格加密处理。在采空区附近，网格尺寸控制在1.5 m × 1.5 m × (层厚)左右，以确保能够精确捕捉应力梯度和塑性区的演化；而在模型外围区域，则采用相对稀疏的网格，从而在保证结果准确性的前提下，显著提升计算效率。最终生成的模型共计约12.8万个网格单元，14.2万个节点。模型的边界条件设置如下：
模型底部(Z = 0)：施加固定铰支约束，即限制所有(X, Y, Z)方向的位移。
模型四周(X、Y方向)：施加法向位移约束，即限制模型边界在法线方向的移动，以模拟无限远岩体对模型区域的约束作用。
模型顶部(Z = 53.8 m)：为模拟上覆岩层的自重应力，施加一个均匀分布的垂直压应力σv。其大小根据上覆岩层的平均容重(约25 kN/m3)与埋深计算得出。
模型的初始应力场主要由岩体自重产生。在计算开始时，首先在重力加速度(g = 9.81 m/s2)作用下生成初始地应力场，并使模型在仅受重力作用下达到平衡，以此作为模拟开采前原岩应力状态的基准。侧压力系数λ根据矿区地应力场特征，取值范围设定为0.8~1.2。
岩土材料的本构关系是数值模拟的核心。考虑到煤矿岩层通常表现为弹塑性行为，且在剪切破坏中表现出明显的应力阈值，本研究选用在岩土工程中广泛应用且物理意义明确的Mohr-Coulomb弹塑性本构模型。该模型通过内聚力c和内摩擦角φ来定义材料的剪切强度，能够较好地模拟岩体在采动应力下的屈服与破坏过程。
模型各岩层的物理力学参数主要来源于实验室实测结果(对于11#煤层及其直接顶板)，部分上覆岩层的参数参考了矿区地质勘探报告中的典型值。所有参数统一整理于表1中，确保了模型输入数据的准确性与代表性。
Table 1. Mechanical parameters table of strata in numerical model
表1. 数值模型岩层力学参数表
岩层名称 |
弹性模量(GPa) |
泊松比 |
密度(kg/m3) |
内聚力(MPa) |
内摩擦角(˚) |
抗拉强度(MPa) |
中粒砂岩 |
10.0 |
0.28 |
2600 |
7.0 |
27 |
2.5 |
砂岩 |
8.6 |
0.30 |
2450 |
3.5 |
27 |
1.8 |
粉砂岩 |
8.6 |
0.28 |
2300 |
3.5 |
18 |
1.5 |
11#煤 |
0.76 |
0.24 |
1230 |
1.54 |
24 |
0.5 |
11#煤顶板 |
5.43 |
0.22 |
2630 |
6.39 |
36 |
2.0 |
粘土岩 |
5.0 |
0.30 |
2500 |
3.0 |
21 |
1.2 |
砂页岩 |
5.0 |
0.29 |
2500 |
3.5 |
17 |
1.0 |
为真实还原井下回采过程，模拟采用分步开挖(null模型)的方法。模拟过程遵循以下步骤：
(1) 初始平衡：在施加所有边界条件和重力后，运行模型直至达到初始力学平衡状态。
(2) 分步开挖：模拟工作面的连续推进。首先，根据初次来压步距20 m的工况，一次性开挖对应长度的煤体单元，并将其材料模型设置为null以模拟采空。计算至再次平衡后，记录应力、位移场数据。
(3) 周期来压模拟：在初次来压的基础上，继续向前开挖煤体，使总推进距离达到周期来压步距30 m，再次计算至平衡并记录结果。这一过程模拟了顶板周期性断裂的关键阶段。
(4) 数值稳定性控制：为确保计算的收敛性与稳定性，采用了FLAC3D默认的自适应全局阻尼与自动时间步长。同时，对关键监测点(如煤壁前方应力、顶板位移)的历史变量进行跟踪，以判断计算的稳定性。
3. 数值试验方案
3.1. 模拟工况设置
基于现场矿压观测确定的顶板来压步距，设置以下两种核心模拟工况：
工况A (顶板初次来压)：模拟工作面自开切眼开始回采，累计推进20 m时的采动状态。此工况旨在揭示基本顶初次断裂前后的应力分布与围岩破坏特征。
工况B (顶板周期来压)：在工况A的基础上，继续模拟工作面推进，直至累计推进30 m (即一个周期来压步距)时的采动状态。此工况用于分析基本顶周期性破断对工作面前方支承压力及塑性区演化的影响规律。
所有工况均遵循上述分步开挖与平衡计算流程，确保每个工况的模拟结果均代表一个稳定的力学状态。
3.2. 观测内容与数据提取
为定量分析两种工况下的矿压规律，设定以下关键观测与数据提取内容：
(1) 垂直应力场：沿工作面推进方向(Y方向)，在煤层中部提取一条垂直应力剖面。重点观测以下数值：
① 应力峰值(σmax)及其距煤壁的距离(Lpeak)。
② 卸压区宽度(Lrelief)：从煤壁至垂直应力恢复至原岩应力值(σ0)处的距离。
③ 应力集中系数(K = σmax/σ0)。
(2) 塑性区分布：输出模型的塑性区(Shear-n, Shear-p, Tension-n)云图。重点分析以下区域：
① 塑性区的范围、形态及主要破坏类型(剪切破坏或拉伸破坏)。
② 煤壁前方及顶底板中塑性区的贯通与发育情况，评估潜在冒落与片帮风险。
(3) 位移场：监测关键点(如煤壁前方顶底板、地表对应点)的垂直位移，为分析顶板下沉与地表沉陷提供数据支持。
3.3. 参数敏感性分析
为评估关键力学参数对模拟结果的影响程度，并增强研究结论的鲁棒性，设计单因素敏感性分析方案。在基准模型的基础上，依次对以下参数进行±20%的浮动，并重新运行工况B的模拟，观察卸压区宽度(Lrelief)和应力峰值(σmax)的变化率：
(1) 11#煤的弹性模量(E)：影响煤体的变形能力及应力分布的梯度。
(2) 11#煤的内聚力(c)与内摩擦角(φ)：共同控制煤体的剪切强度，对塑性区范围尤为敏感。
(3) 直接顶板的弹性模量(E)：影响顶板的弯曲刚度与破断步距，间接影响应力传递。
通过比较各参数变动后结果与基准结果的差异，识别出对矿压分布最敏感的参数，为现场岩性勘探与力学参数测试的侧重点提供参考。
3.4. 模型验证方法
由于获取完整的现场实测数据存在困难，模型验证通过以下多种途径进行：
(1) 定性对比：将模拟得到的塑性区形态(如顶板的垮落、煤壁的片帮)与矿压观测记录中描述的围岩破坏现象进行对比，检验模型在反映宏观破坏模式上的合理性。
(2) 半定量对比：将模拟计算得到的卸压区宽度(14~15 m)、应力峰值位置(16~17 m)与经验值进行对比，评估其是否处于合理的范围之内。
(3) 趋势一致性检验：检验模拟结果是否符合经典的矿压理论，例如，周期来压的应力集中程度是否普遍高于初次来压，以及卸压区随工作面推进是否呈现规律性变化等。
通过上述综合验证，可在较大程度上保证本数值模型及其模拟结果的可靠性与工程参考价值。
4. 结果与分析
4.1. 应力分布特征
(a) 初次来压(20 m) (b) 周期来压(30 m)
A——采煤机割煤；B——伪顶；C——直接顶垮落区；D——老顶。
Figure 1. Comparison of vertical stress distribution cloud diagrams in front of the working face coal wall
图1. 工作面煤壁前方垂直应力分布云图对比
图1分别展示了初次来压与周期来压时，模型中部剖面的垂直应力云图。对应力曲线进行量化提取与分析，得到关键参数对比如下：
(1) 卸压区宽度(Lrelief)：在初次来压工况下，卸压区宽度约为14 m；而在周期来压工况下，该宽度扩大至约15 m。这表明随着工作面推进和采空区范围增大，煤壁前方被压酥破坏的煤体范围有所扩展。
(2) 应力峰值位置(Lpeak))：初次来压时，垂直应力峰值点位于煤壁前方约16 m处；周期来压时，峰值点略微后移至煤壁前方约17 m处。峰值位置的迁移反映了高支承压力区随着采动影响加剧而向煤体深部转移的趋势。
(3) 应力集中系数(K)：初次来压时的应力集中系数约为1.8，而周期来压时增至约2.1。这一结果印证了周期来压期间顶板活动更为剧烈，作用于煤体上的支承压力显著增强的现场观测事实。
4.2. 塑性区演化规律
图2展示了两种工况下的塑性区分布。塑性区范围直观地揭示了围岩的破坏状况。在初次来压时，塑性区主要集中于煤壁前方的卸压带内以及直接顶的垮落区(图2(a)中C区)。煤体以剪切破坏(Shear-p)为主。在周期来压时，塑性区范围明显扩大。不仅煤壁前方的剪切破坏区向深部延伸，顶板中的塑性区也显著增高、增宽，老顶(图2(b)中D区)开始出现塑性破坏，表明其参与并影响了周期来压的矿压显现。这种塑性区的扩展是导致周期来压强度大于初次来压的直接原因。
塑性区的形态与范围为确定巷道合理护巷煤柱宽度和锚杆索支护的有效长度提供了关键依据，建议将主要支护结构锚固在塑性区以外的稳定岩体中。
(a) 初次来压(20 m) (b) 周期来压(30 m)
A——采煤机割煤；B——伪顶；C——直接顶垮落区；D——老顶。
Figure 2. Comparison of plastic zone distribution in surrounding rock of working face
图2. 工作面围岩塑性区分布对比
4.3. 对地表沉陷预测的启示
通过监测模型地表点的垂直位移发现，在周期来压工况下，地表的最大下沉值相较于初次来压工况增大了约25%。模拟显示的沉陷盆地中心与采空区中心基本对应。基于模拟得到的卸压区边界和顶板塑性破坏高度，可以推断出地表的主要影响边界。建议将煤壁对应地表位置向采空区外侧延伸15~20 m的范围划定为重点监测与预复垦区域，此范围内的地表将经历最显著的拉伸与压缩变形。
4.4. 对煤层透气性演化的影响分析
卸压区的塑性破坏与大量裂隙的生成，将极大改变煤体的渗透特性。本模拟显示的14~15 m的卸压区宽度，为布置本煤层瓦斯抽采钻孔提供了关键参数。在此区域内，由于煤体由三向应力状态转变为双向或单向应力状态，其体积膨胀(扩容效应)，裂隙率大幅增加。此状态下煤层的透气性系数可比原始煤层提高数十倍至数百倍。因此，将瓦斯抽采钻孔的终孔位置设计在卸压区范围内，并利用周期来压后透气性增大的时机强化抽采，可望显著提高瓦斯抽采效率，实现煤与瓦斯的协同安全开采。
4.5. 敏感性分析结论
以工况B为例，参数敏感性分析结果表明，煤体内摩擦角(φ)和内聚力(c)是对卸压区宽度(Lrelief)和塑性区范围最敏感的参数。当c和φ值降低20%时，Lrelief增大约18%，塑性区显著扩大。这提示在现场生产中，若遇煤质松软、强度较低的区段，应提前加强支护并注意瓦斯涌出变化。煤层弹性模量(E)主要影响应力分布的梯度，对其峰值大小和位置有次要影响。顶板弹性模量的变化对煤壁前方的应力集中系数有较明显的影响，弹性模量越高，应力集中越显著。
此分析结论强调，在矿井地质勘探与煤岩力学性质测试中，应优先确保煤体强度参数(c, φ)和顶板刚度测定的准确性，以提升矿压预测与控制的可靠性。
5. 结论与工程意义
5.1. 主要结论
本章通过构建基于实测力学参数的FLAC3D三维数值模型，对后安煤矿110203综采工作面在顶板初次来压与周期来压过程中的矿压分布规律进行了系统模拟与分析，得出以下主要结论：
(1) 模拟结果清晰地揭示了采动影响下煤壁前方依次形成的卸压区、应力集中区和原岩应力区。定量分析表明，初次来压时卸压区宽度约为14 m，应力峰值位于煤壁前方16 m处；而周期来压时，卸压区宽度扩大至约15 m，应力峰值位置略微后移至17 m处，且应力集中系数由1.8增至2.1，证实了周期来压的剧烈性。
(2) 塑性区分布显示，周期来压相较于初次来压，其破坏范围显著扩大，不仅煤壁前方的剪切破坏区向深部延伸，顶板中的塑性区也显著增高、增宽，老顶开始参与并影响矿压显现。这为顶板管理与围岩控制提供了直接的力学依据。
(3) 参数敏感性分析表明，煤体的内聚力c和内摩擦角φ是对卸压区宽度和塑性区范围最敏感的力学参数，其微小变化会引起围岩破坏状态的显著改变。其次是顶板的弹性模量，它主要影响应力集中的程度。
(4) 模拟结果表明，14~15 m的卸压区是煤层透气性显著增加的区域，为瓦斯高效抽采提供了最佳区位。同时，基于塑性区范围和应力传递路径，可将煤壁对应地表外侧15~20 m的范围划定为地表沉陷监测与复垦治理的重点区域。
5.2. 工程意义
本研究的结果对后安煤矿及类似条件矿井的安全生产与绿色开采具有直接的指导意义：
(1) 建议将本煤层瓦斯抽采钻孔的密集布置区和终孔位置确定在煤壁前方15~20 m的卸压带内，并充分利用周期来压后煤层透气性倍增的时机实施强化抽采，以有效治理瓦斯灾害。
(2) 应根据周期来压期间更大的应力集中和更宽的塑性区范围，来优化工作面的液压支架工作阻力和巷道的锚杆索支护参数，确保支护系统能够有效控制顶板并抑制煤壁片帮。
(3) 建议将沉陷区的发育与工作面回采进度相结合，优先对模拟预测的重点沉陷范围进行动态监测与整治。