基于FLAC3D某铁矿采空区稳定性数值模拟研究

龙 斌 徐 超 戴 军 卿自强

（1.长沙矿山研究院有限责任公司；2.金属矿山安全技术国家重点实验室）

随着我国经济水平的持续提升，我国的非煤矿山获得了较大程度的发展，矿石开采规模的扩大，导致了我国许多矿区开始逐渐出现大面积的采空区。随着矿山向深部开采，地压增大，采空区在高应力作用下，容易发生坍塌事故。残留大量的采空区没有进行及时处理，给矿山开采带来了严重的安全隐患，同时给矿山作业人员和设备带来严重的威胁，对自然环境造成一定程度的破坏［1-2］。因此，如何有效地治理采空区，避免因采空区失稳给矿山带来灾害，是目前矿山安全生产的亟需解决的技术难题。

目前，国内一些学者对地下矿山采空区稳定性展开了相关研究工作。张佳男等［3］通过对矿山现状采空区情况进行调查、统计和分析，并采用理论计算和数值模拟分析相结合的方式，开展矿山采空区稳定性研究工作，评估采空区安全性，为保证矿山深部安全生产提供技术支撑。周科平等［4］基于某矿山采空区群精细探测数据，采用Geomagic、Midas GTS 与FLAC3D的耦合建模以及数值模拟分析方法，对复杂采空区群进行稳定性分析；同时建立了复杂采空区群的安全分级评价体系，开展了复杂采空区群安全分级评价。谢伟等［5］以湖南某钨矿为研究背景，基于Mathews 稳定性图解法理论模型，对不同中段的采空区稳定性进行了理论计算，并依据矿山采空区现状，采用数值模拟法分析了当前采空区下各区域的主应力、位移和塑性区的分布，评估当前状态下主采空区的稳定状态。余正方等［6］采用三维激光扫描方法，对大红山铁矿Ⅱ-1 头部采空区进行了精细探测，并采用Surpac-Midas/GTS-FLAC3D耦合建模与数值模拟方式，对采空区进行稳定性数值分析。

本研究以某铁矿为例，利用3Dmine-Midas/GTS软件耦合技术，建立了相关地下采场分析数值模型，并采用FLAC3D软件进行了模拟计算，通过分析采空区围岩、顶板、矿柱的应力、应变、塑性区情况，研究矿山南区采空区的稳定性，计算结果可为矿山企业下一步的采空区治理方案提供依据。

1 工程概况

某铁矿采用地下开采方式，目前矿山南区470 m标高以上已开采结束，现有开拓系统为平硐+斜坡道开拓，采用汽车运输，后期开采410，350 m 等中段。矿山目前采用浅孔留矿法回采，平底结构铲运机出矿。采场留有顶柱，采场之间留设有矿柱，跨度较大的采场留设有点柱。南区现有采空区主要位于530，500，470 m中段，500 m中段采空区顶板与530 m中段底板之间留有6～15 m高顶柱，470 m中段采空区顶板与500 m 中段底板之间留有4～10 m 高顶柱。530 m中段现有采空区20 个，编号为S53C01～S53C20，采空区总体积为36.50×104m3；500 m 中段现有采空区20个，编号为S50C01～S50C20，采空区总体积为22.36×104m3；470 m 中 段 现 有 采 空 区15 个，编 号 为S47C01～S47C15，采空区总体积为18.20×104m3。

2 数值模拟过程

2.1 模型材料参数

根据数值模拟计算分析所需要的基本数据，开展了矿岩力学参数室内试验，主要进行了4 项试验：矿岩容重、矿岩单轴抗压力学试验、劈裂抗拉试验、岩体的泊松比与弹性模量试验。根据现场工程地质调查，基于Hoek-Brown 准则利用Roclab1.0 软件对岩石力学参数进行修正后确定岩体力学参数，见表1。

?

2.2 三维模型建立

本次研究的主要内容为矿山南区开采产生的采空区的稳定性，结合矿山实际开采情况，确定本次研究的范围为X=64 800～65 500 m，Y=65 100～65 800 m，Z=300 m～地表；南区模型尺寸为700 m×700 m×（300-地表高度）m（长×宽×高），数值模型建立如图1所示。

模型采用位移约束边界条件，模型4个侧面为水平位移约束，底面为竖直位移约束，即在模型底部施加Z方向位移约束，4 个侧面分别施加X，Y方向位移约束，地表为自由面。

在地下工程数值分析时，计算成功的关键不仅仅需要准确的计算模型，初始地应力场也十分重要。根据该矿区的地形特点和地质构造条件，原岩应力以自重应力场为主，可以不考虑构造应力的存在，其地应力场按自重应力场分析，应力场各主应力均随深度呈线性增加。因此，本次模拟中初始地应力场仅按自重应力场考虑，如图2所示。

2.3 模拟步骤

矿山采用浅孔留矿法回采，实际生产中段间自上而下开采，中段内后退式回采。为了较真实地反映出矿山开采对岩体的扰动影响，数值模拟过程开挖顺序和矿山实际一致。本次数值模拟的过程如下：

（1）在未进行地下开采的情况下，生成初始应力及初始位移。

（2）清除初始位移。

（3）模拟采用浅孔留矿法回采，对采场进行开挖，形成采空区。

（4）模拟开采顺序自上而下逐个中段开挖，南区由南向北开挖。

3 数值模拟计算结果

针对矿区南区复杂的采空区赋存情况，从位移、应力以及塑性区3个方面合理准确地分析模拟结果，从而综合评价南区采空区的稳定性。通过切取剖面分析南区采空区的位移、应力、塑性区状态，本文以勘探线为剖面进行分析。

3.1 南区位移分析

地下开采时，开采区周围岩体受到扰动，引起岩体变形甚至破坏。采空区顶板、上盘围岩大面积冒落是引发地表塌陷、冲击地压灾害的主要因素，给地表设施以及地下生产活动带来极大的安全危害。因此，有必要对各采空区顶板、上盘围岩的变形进行分析，根据采空区的变形程度判断采空区的稳定性，为空区治理提供参考依据。模拟南区530，500，470 m中段开采至目前现状后，部分勘探线剖面围岩变形如图3所示。

从图3 可以看出，矿体开采后，围岩产生向空区一侧的变形，变形较大区域为顶板。同时，可以看出围岩的变形扩展到了地表，使得地表产生了一定程度变形。影响采空区稳定性的主要是沉降变形，南区各采空区围岩沉降变形情况如图4所示。

从图4 可以看出，南区530，500，470 m 中段开采后围岩产生的沉降位移很小，基本在0.1～1.8 mm，由于矿山工程地质条件好，围岩和矿体均为坚硬的岩体，破坏时的变形很小。

3.2 南区压应力分析

矿体开采前，地下岩层处于应力平衡状态。矿体的开采破坏了原有的平衡状态，引起围岩应力重新分布，出现应力状态改变和应力集中。如果围岩新应力场中的应力没有超过围岩的承载能力，围岩将会自行达到新的平衡；否则，围岩将可能产生破坏，出现岩爆、片帮、冒顶等现象，影响采空区的稳定性和矿山正常安全生产。模拟南区530，500，470 m中段开采至目前现状后，围岩压应力如图5所示。

从图5 可看出，矿体开采后，在采空区底部角隅处出现了明显的压应力集中，其次在采空区间的矿柱所受到的压应力相较于开采前也出现了明显的增加。南区各采空区围岩最大压应力统计如图6所示。

从图6可以看出，各中段采空区围岩的最大压应力未出现较大的变化。总体来看，围岩压应力值往深部逐渐增大，最大压应力值未超过岩体的抗压强度，因此围岩不会出现压裂破坏。

3.3 南区拉应力分析

矿体回采后，围岩由于受到开采扰动，可能出现拉应力。通常岩石的抗拉强度为抗压强度的1/4～1/25（平均1/10），因此围岩出现拉应力时需引起重视。模拟南区530，500，470 m 中段开采至目前现状后，围岩拉应力如图7所示。

从图7 可以得出，较多采空区围岩产生了拉应力，其拉应力区域在顶板中间或上盘靠近顶板处。南区各采空区围岩最大拉应力统计如图8 所示。可见，大部分采空区围岩均产生了不同大小的拉应力，拉应力值在0～1.44 MPa，均小于岩体的抗拉强度。因此，围岩不会产生拉伸破坏。同时，对比分析各采空区拉应力，顶板跨度大的采空区，围岩应力相对要更大，比如S50C02、S50C14 等。上、下中段层叠采空区的顶板拉应力值除与采空区跨度有关外，还与顶板厚度相关，主要是因为顶板需承受自身的重力。

3.4 南区塑性区分析

塑性区可以直观分析矿体回采后围岩的破坏规律及采空区稳定性情况，模拟矿体开采后，围岩塑性区分布如图9 所示，可以看出，S53C05 采空区顶板和与S53C06 间的矿柱产生了连续贯通塑性破坏区，表明S53C05 采空区过去出现过顶板冒落、矿柱破坏，可能导致上覆岩层塌陷，以及引起周边S53C06 采空区失稳。现场可见南区S53C05、S53C06采空区上覆岩层已塌陷至地表。530 m中段其余采空区围岩未出现塑性区或仅有零星塑性区，对采空区稳定性影响很小。500，470 m中段围岩未出现塑性区，说明稳定性较好。

4 结 论

以某铁矿为例，利用3Dmine-Midas/GTS 软件耦合技术，建立了相关地下采场分析数值模型，并采用三维数值模拟软件FLAC3D软件进行了模拟计算，并对矿山南区采空区围岩、顶板、矿柱的应力、应变、塑性区情况进行了分析，得出如下结论：

（1）南区530，500，470 m 中段开采后围岩产生的沉降位移很小，基本在0.1～1.8 mm。

（2）南区各中段采空区围岩的最大压应力未出现较大的变化。总体来看，围岩压应力值往深部逐渐增大，最大压应力值未超过岩体的抗压强度，因此围岩不会出现压裂破坏。

（3）南区大部分采空区围岩均产生了不同大小的拉应力，拉应力值在0～1.44 MPa，均小于岩体的抗拉强度。因此，围岩不会产生拉伸破坏。

（4）530 中段S53C05 采空区顶板和与S53C06 间的矿柱产生了连续贯通塑性破坏区，南区其余中段采空区围岩未出现塑性区或仅有零星塑性区，对采空区稳定性影响很小。

（5）综合各采空区围岩位移、应力、塑性区情况，南区除S53C05、S53C06采空区稳定性较差外，其余采空区总体的稳定性较好。