李鸿飞，汪为平，肖益盖，刘海林，刘涛影

(1.中钢集团马鞍山矿山研究院有限公司，安徽 马鞍山 243000；2.华唯金属矿产资源高效循环利用国家工程研究中心有限公司，安徽 马鞍山243000；3.中南大学 资源与安全工程学院，长沙 410083)

近年来，随着采空区治理[1]的需要，采空区实测技术已经发展成熟，同时国内外许多学者对采空区稳定性数值模拟技术进行了广泛的研究，罗周全[2]和刘海林[3]分别采用Surpac和3DMine建立了实测采空区块体模型进行稳定性分析，韦文蓬[4]采用3DMine、Midas和FLAC3D对不同尺寸下的理想采场进行方案优化模拟，罗广强[5]根据单一溜井测量点云进行了溜井稳定性分析等。采用理想模型虽然可以对采空区稳定性规律进行研究，但实际矿山采空区数量较多、层次结构复杂，构建实测数值模拟模型困难，进而导致数值模拟时模型的精度较差和计算效率低。

本文介绍了一套高效、精准的实测采空区稳定性分析技术，以北辛庄矿区为例，通过对采空区群三维实测，经3DMine对模拟区域多组剖面线切割及修正，导入Midas建模和网格划分，利用Midas-FLAC3D耦合建模，最后采用FLAC3D采空区群数值模拟分析。

1 采空区群三维实测及建模

北辛庄矿区位于山东省兰陵县，属于倾斜急倾斜，薄至中厚矿体，倾角45°～67°，矿体厚度1.08～12.47 m，平均厚度5.48 m，矿山采空区分布在-40 m、-80 m、-108 m和-165 m四个中段内。

采空区群稳定性模拟的准确性取决于实际测量模型的准确性及矿山地质准确建模。本次采空区测量采用FARO三维激光扫描仪，测量的数据经FARO专用软件Scene和studio处理成DXF点云文件，将DXF点云文件导入3DMine中进行封装及三维建模，全矿区采空区群三维模型见图1。

图1 矿山实测采空区群三维模型Fig.1 Three-dimensional model of mine measured goaf group

由于矿山采用留矿法开采，采空区内留有大量残矿，扫描点云封装后的实体模型与实际采空区轮廓存在误差，对采空区群模型进行多组剖面线切割，保留剖面线并结合采场实际开采现状对其进行修正，模拟区域实测采空区群实体模型及修正多组剖面线见图2。

Midas的几何建模和网格划分功能强大，FLAC3D为常用的岩土工程数值模拟工具，二者结合可以高效的进行岩土工程数值模拟[6-12]。将在3DMine中修正后的采空区多组剖面线.DXF文件导入至Midas中，进行剖面线平面建模，再将各剖面放样实体模型；采用地形数据生成器生成模拟区域地表模型及地层模型。对实体模型进行网格划分，划分时，采空区周边线条进行的尺寸控制参数为3，全矿模型边界尺寸控制参数为<23，采用三角网格划分方式，模型长485 m、宽286 m、高约310 m，共划分57 970节点和614 640单元，矿区实体模型及其网格划分见图3。

图3 模拟矿区实体模型和网格Fig.3 Solid and mesh model of simulated mining area

将划分网格的模型导出GTS中型.fpn文件，转换为.fLAC3D文件，将转换后的网格文件导入到FLAC3D中，见图4，并进一步进行实测采空区群稳定性分析，从而实现Midas到FLAC3D的耦合。

图4 FLAC3D数值模拟模型Fig.4 FLAC3D numerical simulation model

2 采空区群稳定性分析

2.1 本构模型

FLAC3D是岩土数值模拟常用工具，在进行数值模拟时，将岩体材料假定为理想弹塑性材料。模拟采用了Mohr-Coulomb强度准则，强度准则和塑性之间的解析关系式如下。Mohr-Coulomb准则对fs=0定义：

(1)

式中，Nφ、NΦ分别表示岩体材料的内摩擦角和剪胀角；c表示岩体材料的黏聚力，并且有：

(2)

(3)

(4)

(5)

(6)

式(3)～(6)中：

(7)

式(7)中，K为岩体体积模量，μ为岩体泊松比，E为岩体弹性模量，G为岩体剪切模量。

2.2 模拟参数确定

通过现场取样，进行室内岩石力学试验，测试项目包括：1)岩石物理性质(密度)试验；2)岩石单轴抗压强度试验；3)岩石单轴变形试验；4)岩石抗拉强度试验；5)抗剪断试验。获取了矿岩的密度、单轴抗压强度、抗拉强度、弹性模量、泊松比、黏聚力和内摩擦角等参数，并采用广义Hoek-Brown强度折减准则进行折减[13-16]，折减后的矿岩物理力学参数见表1。

表1 矿岩物理力学参数表Table 1 Mechanical parameters of ore body and surrounding rock

2.3 模拟结果分析

模型导入FLAC3D中后，首先施加自重应力进行求解，收敛后，进行采空区开挖求解。分别绘制出采空区群周边最大主应力云图，采空区群围岩垂直变形云图，采空区群塑性破坏区分布云图。见图5。

模拟结果显示此种方法的数值分析效率较高，模型在5 432步求解后达到平衡状态；拉应力区主要分布于采空区上下盘和顶板，大小与暴露面积正相关，顶板及上盘可能会造成采空区局部围岩掉落；压应力主要集中于采空区群间柱及靠近上盘的采空区底部区域和靠近下盘的采空区顶部区域，最大压应力5.54 MPa，未超模拟矿岩最小抗压强度；采空区上下盘由于受拉出现最大2.18 cm垂直位移；采空区群整体未出现大规模塑性区，仅部分中段由于采空区形状不规则导致局部出现shear-p型塑性破坏区。

2.4 模型验证

上述数值模拟虽然可以快速对采空区稳定性进行分析，但模型的准确性需要验证。本模型验证通过对比数值模拟地表位移量与实际地表监测位移量二者之间的误差来判定。

FLAC3D在模拟开挖和采场充填时，仅能监测到节点位移。因此，根据采空区群地表网格划分现状，设置网度为5 m×6 m共30个监测点组成的地表位移数值模拟位移监测网，并将监测网附近最近节点设置为位移监测点，监测点坐标及数值模拟最大水平位移和垂直位移见表2。

将模拟地表监测点最大垂直位移进行三维拟合，见图6。与矿山实际地表位移量数据进行对比，结果显示地表模拟位移量与实际监测值误差不超过10%，地表未出现大的变形，符合矿山实际情况。因此本次模拟所构建的矿山采空区群数值模型较为准确。

表2 数值模拟位移监测点详细信息表Table 2 Details of the displacement of monitoring points

图6 地表监测垂直位移量图Fig.6 Monitored vertical displacement

3 结论

1)基于Faro-3Dmine及Midas-FLAC3D耦合的稳定性分析技术通过三维扫描仪精确扫描空区轮廓并根据实际采场情况进行模型修正，建立了精准的三维采空区模型，再采用有限差分法进行稳定性分析并通过地表监测数据对模型进行了有效性检验。结果显示此方法可以快速构建精准的复杂采空区群模型，高效的对复杂采空群稳定性进行分析，为矿山复杂采空区群安全治理提供了可靠的理论依据。

2)稳定性分析显示实测复杂采空区群采空区形状不规则会导致局部应力集中；采空区顶板及上盘采空区局部变形较大；上下中段贯通的采空区暴露面积大，应力较为集中且变形较大，应对其重点监控并优先治理。