韦文蓬，梁东强，冯思成，李鸿飞

(1.中信大锰矿业有限责任公司 大新锰矿分公司，广西 大新 532315；2.中钢集团 马鞍山矿山研究总院股份有限公司，安徽 马鞍山 243000)

广西大新锰矿存在采空区体积庞大，存在时间较长的问题，采空区稳定性问题成为矿山亟待解决的难题[1-3]，常规的采空区模糊定性分析不能准确反映采空区稳定性，本文采用3DMine[4-5]、Surfer、Rhino、ANSYS和Flac3d[6-7]等工具，结合矿山采空区实际情况，对矿山采空区的稳定性进行了详细的数值模拟分析。矿区采空区平面分布见图1。

图1 大新锰矿采空区分布平面示意

1 数值建模

以3DMine为基础，Surfer软件和Rhino造型软件为依托，将所构建的地表、矿体、空区和矿柱经过这些若干软件不同功能的耦合以及各种不同格式的转换，最终组合出ANSYS-workbench中可识别的犀牛模型，利用ANSYS进行矿区模型构建和网格划分，经转换后导入Flac3d之中，完成建模。数值建模流程见图2。

图2 大新锰矿采空区稳定性数值模拟建模分析流程

1.1 3Dmine建模

采用3DMine软件的表面建模功能构建地表面模型，并生成面模型层面点集，见图3。

图3 地表模型和层面点集

1.2 Surfer点数据处理及Rhino建模

将在Excel预处理后的点数据导入Surfer，生成等高线文件，经surfer多次处理后，导出txt格式文件，将txt文件导入Rhino中，生成矿体表面、空区表面等面模型，生成组合面模型。将组合模型导出iges文件。Surfer点处理和Rhino建模结果见图4。

图4 surfer数据处理和Rhino建模

1.3 ANSYS建模与网格划分

将Rhino导出的iges格式文件里的组合模型导入ANSYS软件中，利用‘Generate’、‘Extrude’、‘Slice’、‘Boolean’等功能生成包含矿柱、矿体以及空区的实体模型，再采用Mesh模块进行模型网格最优化划分，ANSYS建模和网格划分结果见图5。

图5 实体建模与网格划分

将网格划分后的模型导出dat格式单元和节点文件，经数据转换工具生成Flac3d文件，在Flac3d中打开此文件，即可进行数值模拟计算。

2 Flac3d有限差分法数值模拟

2.1 破坏判据

2.1.1 围岩强度理论

岩体破坏主要分为拉伸破裂和剪切破坏两种类型。在低围压下呈拉伸破裂，围压增高后，出现剪切破裂。由于单用一种强度理论，不可能解释所有情况下岩石的破坏问题，针对岩体在不同应力条件下发生的破坏现象和破坏形式，采用的判据有抗拉破坏准则、Mohr-Coulomb准则和抗拉强度准则进行评判。

材料的抗剪强度用Mohr-Coulomb准则[8]来表示：

τ=σ·tanφ+C

(1)

式(1)中，τ为材料的抗剪强度；σ为斜截面上的正应力；φ为材料的内摩擦角；C为材料的凝聚力。

材料在拉应力状态下时，采用拉破坏强度准则。最大拉应力准则常用于脆性材料的受拉破坏，设σ1，材料抗拉强度为f1，则准则为：

f=σ1-f1=0

(2)

如果拉应力超过材料的抗拉强度(f≥0)，材料将发生拉破坏。

2.1.2 容许极限位移量判据

甯瑜琳[9]给出了岩体位移量与其对地下采场大跨度开挖不进行支护或临时支护条件下岩体稳定性关系，见表1。

表1 岩体位移量与岩体稳定性关系

2.1.3 塑性区贯通法

陈欣等[10]给出了塑性区贯法失稳的判据表达形式：

φ(Ai)<0

(3)

式(3)中，φ为屈服面；Ai为由局部区φj组成的机构。

屈服区域是由几个局部屈服区φj连通而形成的，可以表示为：Aj=φj∪φm∪…φn，即只要某几个局部屈服区φj能构成一个破坏区域使工程失去整体稳定性，则这几个局部屈服区即可构成一个Ai，当屈服区连通形成任意一个Ai时，结构整体失稳。

2.2 矿岩物理力学参数确定

通过现场取样，进行了室内岩石力学试验，获取了矿区顶板、底板、矿体和夹层的密度、抗压强度、抗拉强度、弹性模量、泊松比、体积模量、剪切模量、粘聚力、内摩擦角等参数，采用多种方法对岩石力学参数进行折减。数值模拟模型物理力学参数见表2。

表2 折减后矿岩体物理力学参数汇总

2.3 初始条件

本研究计算模型上表面为地表，边界条件的采用自由边界，底面采用固定约束，x和y方向边界施加铰支约束，即模型在计算过程中不产生边界法线方向的位移，而可允许其它方向的位移[11]。在模拟计算的过程中将自重应力场设置为初始应力场。

2.4 结果分析

2.4.1 中部采区顶板及矿柱数值模拟结果分析

中部采区14线应力分析结果和位移分析结果见图6～7，塑性区分布状况见图8。

图6 中部采区主应力云图

图7 中部采空区位移云图

图8 中部采空区塑性区分布

中部采区数值模拟结果显示：矿块回采后，上部覆岩拉应力最大1.1 MPa，下沉量最大为13.6 mm；矿柱出现压应力最大值10.2 MPa，空区底鼓变形最大值为15.4 mm。空区上下盘围岩水平位移最大值为12.1 mm。回采结束后，采场顶板和矿柱区域因受较大的剪应力进入塑性状态，且下部中段个别矿柱内塑性区已贯通，矿柱将屈服破坏。

2.4.2 西北采区顶板及矿柱数值模拟结果分析

西北采区应力分析结果和位移分析结果见图9～10，塑性区分布状况见图11。

图11 西北采区塑性区分布

图9 西北采区主应力云图

西北采区数值分析结果显示：西北采区矿块回采结束后，采空区顶底板拉应力集中，最大值为1.01 MPa，矿柱内出现了不同程度的应力集中分布区，最大压应力为19.11 MPa。最大下沉量出现在上盘区域，且下沉量为8.6 mm。地表下沉量最大值为7.8 mm。当矿床回采结束后，采场出现了塑性区，下部的点柱塑性区已贯通，点柱将屈服破坏。

图10 西北采空区位移云图

3 结 论

1)基于3Dmine-Rhino-Ansys-Flac3d分析方法可以对采空区稳定性进行较为合理的半定量分析，从而为矿区治理采空区的实践提供理论依据。

2)稳定性分析表明：中部采空区和西北采空区顶板及上部覆岩呈现拉应力，矿柱内呈现压应力集中现象，采空区间柱塑性区贯通，可能会造成失稳，建议矿区应减少采空区暴露时间，尽快对采空区进行充填或加固。