基于DIMINE-MIDAS/GTS的采空区稳定性分析

任青云

(湖南有色金属研究院，长沙 410100)

地下矿山开采的重点安全问题是开采形成的采空区对后续开采产生的岩层移动问题，开采过程中伴随的地层沉降及地压活动[1-2]。矿山开采面临的重要安全问题主要表现在，开采区域的不断扩大，地层随之进行持续的运动，由于开采现场会出现不同的地质条件，例如溶洞、断层、空区等可能会产生移动变形的大范围空间，此时会造成地下矿山的地压显现，导致顶板出现大变形，继而影响围岩发生冒落和片帮等现象。因此地下矿山长期开采导致的采空区稳定性是影响矿山生命财产安全及后续开采的重要因素[3-4]。

影响采空区的稳定性及矿山地表位移沉降的主要原因有开采矿体的埋深、采场的高度，矿体的厚度、倾角、岩石力学特性及矿山采用的采矿工艺。采用大型有限元数值模拟[5-6]前处理模式可以直观地表示出地下采空区的分布情况及整改矿山的轮廓，其后处理模式能够计算出矿山的地表位移沉降、采空区是否发生拉应力破坏、矿柱是否产生主应力集中及塑性区变形破坏情况。本文根据南川河石灰岩矿矿区地质资料，结合矿体情况及开采设计，对南川河石灰岩矿进行取样试验，采用数值模拟方法分析矿山开采后岩层移动变形情况，并结合传统的理论计算方法对矿柱、顶板稳定性进行分析，以此确定无支护状态下采空区跨度，从而评判采空区的稳定性及地表位移沉降的影响范围和程度。

1 工程概况

南川河石灰岩矿位于浏阳市澄潭江镇南川河畔，距浏阳市城区直距 21 km，公司始建于1972年，原采用露天开采，1996年开始地下开采，主要采用斜坡道开拓，对角式通风。矿山长期使用的采矿方法是房柱采矿法。矿山总的回采顺序是沿矿体倾斜由上而下循阶段进行，按边探边采的原则进行开采。采场回采顺序是后退式开采，矿房内采用分层开采。矿山2014年之前开采130 m中段，2014年至今开采100 m中段，目前矿山正在+100 m中段北部进行生产。

由于石灰岩的岩性较软，并且矿山开采历史悠久，地下空区存在极大的不稳定性及不确定性，为了掌握空区对下部矿体回采的影响，必须针对采空区进行稳定性研究和分析。根据矿山的实际生产现状，针对空区进行了现场调查，并通过统计分析获得矿山不同开采水平的采空区规模、尺寸。此次现场调查针对+130 m中段和+100 m中段共统计出130个采空区，+100 m中段采空区部分调查数据见表1。现场调查发现，矿山部分采空区已形成采空区群；大部分空区宽度约12～14 m，采空区未发生冒顶片帮事故，侧壁处理较好，有极少量的浮石[7]。溶洞构造溶蚀形成，并有软塑性、硬塑溶洞泥质充填，与地表水源连通的构造少。

表1 +100 m中段采空区特征统计表(部分)

2 矿岩力学参数室内试验

2.1 矿岩力学参数试验

为了合理评价南川河石灰岩矿工程地质条件，为后续数值模拟提供依据，开展了矿岩力学参数室内试验，测试的内容包括矿岩的容重、单轴抗压强度、抗拉强度、弹性模量、泊松比等基本物理力学参数。针对矿山实际情况，选取了石灰岩岩样进行加工，并对各岩样进行了相关参数测定。岩样总计22个，其中5个用于开展单轴压缩试验、5个用于巴西劈裂试验、12个用于变角剪试验。

根据数值模拟计算分析需要的基本数据，主要进行了如下四项内容的实验室试验：1)矿、岩的天然容重；2)矿、岩的单轴抗压力学特性；3)岩体的劈裂抗拉强度试验；4)岩体的泊松比与弹性模量试验。

试验开始前，先对岩样加工。首先，按照标准岩石力学试验所要求的规格，制成试样尺寸为Φ50 mm×100 mm的标准岩石试样；然后擦干试样，在自然状态下风干；最后，试验在250 t全数字型液压伺服刚性岩石力学试验系统(RMT150，图1)及30 t万能材料试验机上完成。

图1 RMT150型试验设备Fig.1 RMT150 test equipment

2.1.1 单轴抗压试验

试样尺寸为直径Ф50 mm，高径比约为2∶1，岩石单轴抗压试验结果如表2所示。

表2 岩石单轴抗压强度及静力受压弹性模量和泊松比测试

2.1.2 抗拉强度试验

岩石的抗拉强度试验如图2所示，各种试验结果见表3和表4。

图2 巴西劈裂试验Fig.2 Brazilian splitting experiment

表3 各试样劈裂抗拉试验参数测定值

表4 变角剪试验结果

根据单轴抗压试验和单轴抗拉试验，将室内试验结果汇总如表5所示。

表5 矿岩力学参数试验结果汇总表

2.2 岩体强度参数分析

由于自然岩体内部有很多缺陷，比如节理、裂隙和孔洞等，因此自然岩体与小试样岩石的力学性质是有差异的，导致实验室测得小试样岩石力学参数不能直接作为数值模拟的力学参数。因此，岩体力学工程界，大多数依据室内岩石力学实验，采用理论力学的方法，估算出岩体的力学参数，用于有限元数值模拟计算。

根据HOEK，CARRANZA-TORRES提出的岩体破坏经验准则：

(1)

岩体的mb、S和a由GSI和上式计算的mi确定。

(2)

(3)

(4)

根据确定的mb、S和a值，可以计算出岩体的单轴抗压强度和抗拉强度：

σc=σci·sa

(5)

(6)

岩体的内摩擦角φ和黏聚力c可以由以下公式进行计算：

(7)

(8)

针对南川河石灰岩矿工程地质条件，取石灰岩GSI=70，根据室内岩石力学参数试验结果，利用HOEK和CARRANZA-TORRES提出的岩体破坏经验准则及岩体分类RMR值，采用Roclab软件进行计算。选取岩体的力学参数如表6所示。

表6 岩体物理力学参数汇总表

3 采空区稳定性数值分析

3.1 三维矿体模型的建立

制约矿山地表位移沉降及采空区稳定性的重要因素，包括矿区地表的地形地质条件及井下采空区的分布状态，因此在三维有限元数值模拟计算过程中，缺少完整的地形模块，不能真实地反映矿山的实际情况。本研究通过基于DIMINE强度的前处理功能，利用二维地质剖面图或者MAPGIS地形数据库生成三维地质模型，对勘探线剖面进行初步的处理，仅保留勘探线与矿体解译线[8]。将矿山剖面图经过CAD处理之后导入到DIMINE中，进行线串运算使勘探线“立起来”，对立起来的三维线文件进行坐标运算，使各点坐标为其真实坐标值。通过上述基本工作的完成，可以联接各矿体解译线剖面，对于比较复杂的矿体解译线，可以构造中线，通过中线与剖面线来生成实体，矿体实体模型见图3。在DIMINE中生成的矿山实体模型，无法直接用于有限元数值模拟的运算，需要将此模型导入到MIDAS/GTS中，利用MIDAS/GTS的网格划分功能，结合采矿工艺，进一步建立采区模型，划分网格单元，进行数值模拟分析。

图3 南川河石灰岩矿三维实体模型Fig.3 3D solid model of Nanchuanhe limestone mine

3.2 理论基础及相关参数设定

根据矿山的地层参数，确定有限元数值模拟的岩体力学参数见表6。计算过程矿岩及废石均采用摩尔-库伦((Mohr-Coulomb)屈服准则，该屈服准则的控制方程为：

(9)

最大拉应力屈服准则函数为：

ft=σ3-σt

(10)

式中：c—黏聚力，MPa；φ—内摩擦角，(°)；σ1—最大主应力，MPa；σ3—最小主应力，MPa；σt—抗拉强度，MPa。MIDAS/GTS有限元计算程序，默认当fs大于0时，表明岩体已发生了剪切破坏；当ft大于0时，表明岩体已发生了拉伸破坏。

为了模型建立与分析结果的准确性，本数值模拟进行了以下假设：模型的矿岩体为均质、各向同性材料，模拟的三维网格模型底部为固定约束，四周为限制水平约束，只在竖直方向有位移变形。初始阶段考虑地应力与自重的影响，位移清零，表明矿山在原始状态已处于稳定的固结沉降。计算模拟不考虑地震、爆破等动荷载的影响，不分析地下水的影响；矿山地质条件较复杂，矿岩稳固性在空间分布上具有较大的随机性，矿体为断层分割，且矿岩和围岩中均有穿插，这些地质现象的影响在岩体参数折减时已做考虑，在模拟过程中不予另行考虑。

根据地下采空区对矿山围岩的影响范围，设置有限元数模网格模型的边界固定条件如下：1)模型X方向的边界为X方向固定边界，取u=0，v≠0，w≠0(u为X轴方向的位移，v为Y轴方向的位移，w为Z轴方向的位移)；2)模型Y方向边界为Y方向固定边界，取u≠0，v=0，w≠0；3)模型底部为全固定边界，取u=0，v=0，w=0；4)模型地表部分为自由面，不限制位移约束。矿体边界约束情况及三维网络质量情况见图4和图5。

图4 三维网格模型+边界条件Fig.4 3D grid model + boundary conditions

图5 网格质量检查Fig.5 Grid quality inspection

3.3 数值模拟计算结果分析

3.3.1 位移场分析

采空区的稳定性分析，位移是评判矿山稳定性的一个重要指标，因为大量的地表变形，会使矿山发生地表塌陷，且位移变形量超过岩体的变形限值，会发生采空区顶板冒顶片帮事故。从图6～9的位移云图得出，该矿山产生的位移变化情况，在矿体的顶板位移值小于0，说明由于产生了采空区，发生了顶板下沉变形情况，但随着矿体的开挖过程进行，在采空区的底部及巷道的底板，竖直方向位移大于0，是因为矿体开挖深度加大，采空区及巷道周围发生了主应力集中现象，水平方向的应力大于垂直方向的应力，使巷道底板产生上拱的位移变形。

从图6竖向位移云图可以看出，+130 m中段矿体开采时，采空区顶板竖向位移最大值为-1.74 mm，矿体开采底板巷道发生向上最大位移为+2.39 mm。

图6 +130 m中段顶板位移云图Fig.6 Displacement cloud image of the +130 m middle section roof

图7为矿柱的位移云图，+130 m中段采空区矿柱位移变形为受压沉降，最大位移值为-1.8 mm，+100 m中段矿柱由于+130 m中段矿体已全部回采，+100 m矿柱所受侧向压力大于竖向压力，变形为部分矿柱产生上拱位移，最大位移值为+2.5 mm。

图7 矿柱位移云图Fig.7 Cloud diagram of pillar displacement

图8为+100 m中段矿体开采后采空区顶板产生的位移沉降变形云图，目前在+100 m中段回采时采空区顶板最大位移沉降为-2.4 mm，采空区底板巷道上拱位移值为+2.94 mm。

图8 +100 m中段空区顶板位移云图Fig.8 Displacement cloud image of the +100 m middle section roof

采空区形成的地表沉降有微弱的影响，地表沉降值不足1 mm，该位移沉降整体上可以忽略，采空区地表沉降等值线范围见图9，位移变形值数据见表7。

图9 地表沉降位移云图与沉降等值线图Fig.9 Surface displacement cloud map and subsidence contour map

表7 采空区位移变化值计算结果

3.3.2 应力场分析

矿体开采后对围岩造成的影响是一个非线性的不可逆的加载过程，它使得处于初始地应力状态下的围岩进行应力重新分布，最后达到新的平衡，矿体开挖引起径向应力释放、切向应力增加，导致矿体开采顶板和底板产生拉应力；从应力云图可以看出矿体开采区域与矿柱接触处发生较大的应力集中，周边最大。采空区影响范围一定深度后，逐步变化为恢复到原岩应力状态。采空区形成后的主应力值都小于0，表明主应力(压应力)为受压状态。

图10 最小主应力云图Fig.10 The minimum principal stress cloud diagram

图11 最大主应力云图Fig.11 The maximum principal stress cloud diagram

根据计算结果图10、图11得出结论：在采空区的顶板和底板产生拉应力破坏，其中+130 m中段采空区顶板产生的最大拉应力约为2.27 MPa；+100 m中段采空区顶板产生的最大拉应力为1.85 MPa，已接近岩体的抗拉强度，表明在采空区的顶板部分已发生拉应力破坏，容易产生冒顶片帮事故。从最大应力云图可以得出：在矿柱与空区的边界围岩出现压应力集中现象，最大压应力集中发生在矿山+100 m中段东北方向的采空区与13号矿柱交界处，产生的最大压应力值为12.66 MPa，+130 m中段的(老采空区)矿柱发生最大压应力集中值为8.88 MPa，表明矿柱发生受压破坏。应力结果见表8。

表8 最大、最小主应力值计算结果

3.3.3 塑性区分析

图12为分析得出的采空区塑性变形情况，从计算结果可以得出，在采空区的附近发生了塑性区变形破坏，塑性区面积占总面积约44.5%，最大塑性变形为2.06×10-3，矿体开采采空区周边矿柱仅有少量单元发生塑性破坏，塑性区变化主要出现在采空区顶板围岩。

图12 塑性区云图Fig.12 Cloud image of plastic zone

4 结论

利用DIMINE强大的前处理功能，建立了南川河石灰岩矿的三维地质模型，导入到MIDAS/GTS，形成了包含采空区、矿体、地表地形的三维有限元数值模型。利用基于DIMINE-MIDAS/GTS的耦合技术，分析了矿体开采后矿山的位移变形、应力集中现象及塑性区破坏情况，分析论证该矿山采空区的稳定性，结果表明：

1)+100 m中段采空区顶板产生的最大位移沉降为-2.4 mm，采空区的底板与巷道底板发生最大正向位移+2.94 mm。采空区形成的地表沉降有微弱的影响，沉降值不足1 mm，表明该矿山开采对地表影响较小。

2)采空区顶板发生了拉应力破坏。+130 m中段老采空区最大拉应力为2.27 MPa，接近围岩的抗拉强度值，采空区容易发生冒顶片帮事故。+100 m中段采空区与13号矿柱接触处产生了最大压应力值集中12.66 MPa，矿柱处于压应力破坏。

3)采空区开采后产生的塑性区变形主要分布在采空区的四周边界，产生的塑性区变形范围占总面积的44.5%，其中最大塑性区变形量为2.06×10-3，采空区周边矿柱仅有少量单元发生塑性破坏，塑性区主要存在于采空区顶板。

4)综合考虑岩体位移、应力、塑性形变等计算结果，矿体开采采空区顶板产生少量的竖直位移，地表发生微量的沉降，部分采空区顶板发生塑性变形，只在+130 m中段老采空区发生拉应力破坏，表明矿山采空区基本处于稳定状态。