采空区稳定性数值模拟研究

王永奇,戴 兵

(1.贵州开磷集团矿业总公司， 贵州 贵阳 550302；2.中南大学资源与安全工程学院, 湖南 长沙 410083)

0 引 言

矿产资源安全、高效开采对国民经济的发展至关重要。在上世纪中后期，房柱法、留矿法、阶段矿房法等空场法[1]在地下矿山中广泛运用，取得良好开采效果。受经济因素影响和空区安全隐患意识淡薄，部分采用空场法开采后残留的空区未进行处理。由于部分采空区高度大、跨度大，且空区在受到长时间的风化、爆破振动的作用下，空区周围岩体强度大幅降低，经常出现采空区顶板坍塌、冒落现象，严重影响矿山后续开采中人员、设备的安全。

采空区稳定性制约矿山的发展，因此对其进行研究具有重要意义。目前，对采空区稳定性研究主要通过力学理论计算[2]、数学模型[3,4]和数值计算。由于岩体中存在大量节理、裂隙，通过理论计算所得结果往往存在较大误差，且理论计算相对复杂，难以得到准确的计算结果。数学模型通过数学方法对采空区进行分类评价，虽能获得分类结果，但不能直观体现空区顶板受力状态，且分类结果与实际存在一定误差。而数值模拟计算能根据空区围岩实际情况对力学参数进行调整，模拟精度相对较高。因此本文以开磷矿区内空区为研究对象，拟通过数值模拟技术对采空区稳定性进行研究。

1 矿山概况

开阳磷矿自1958年建矿以来，已连续生产了50多年。以往矿山使用的采矿方法为空场法，采后采空区未进行处理。矿区形成的采空区，导致山体塌陷和滑坡。矿体延伸15 km，从地表可见绵延15 km的山体塌陷区。矿体走向南北，倾向为280°～290°，倾角27°～30°。受F314断层影响，断层下盘矿体产状由北向南逐渐变陡，在W11～W11+1线矿体出现倒转。矿体厚度5.8～7.5 m，南薄北厚，厚度变化不大，平均6.0 m左右。

2 ANSYS数值软件

ANSYS数值软件采用有限元计算程序，集结构、热、流体、电磁场、声场和耦合场分析于一体，能处理线性和非线性复杂力学计算问题，在开发之初便得到广泛应用[5]。

有限单元法基本思路是将连续结构离散成单元，单元上设定节点，使用变分原理建立有限元方程，求解离散域中的自由度问题。有限元分析求解过程主要可分为以下几个步骤[6]。

(1)对结构进行离散化。通过将求解的连续结构化解与分割，转换成有限单元和节点，保证相邻单元力学性能参数连续。

(2)位移插值函数选择。通过对单元位移函数的假设，利用节点位移来表示单元体位移、形变和应力，以对连续体问题进行分析。位移矩阵如下：

{f}=[N]{δ}e

(1)

式中，{f}—单元中任意一点的位移，

[N]—行函数，

{δ}—单元节点位移。

(3)单元力学特性。首先推导出单元应变(利用节点位移)表达式：

{ε}=[B]{δ}e

(2)

式中，{ε}—单元应变，

[B]—单元应变矩阵。

根据本构方程，得出单元应力：

{σ}=[D][B]{δ}e

(3)

式中，[D]—弹性矩阵(与单元材料相关)。

最后根据变分原理，推导出单元节点力关于节点位移的关系式：

[F]=[k]e{δ}e

(4)

式中，[k]e—单元刚度矩阵：

[k]e=∭[B]T[D][B]dxdydz

(5)

(4)根据(1)、(2)、(3)、(4)、(5)平衡方程，建立整体连续结构的平衡方程：

[K]{σ}=[F]

(6)

式中：[K]为总刚度矩阵。

(5) 计算求解节点位移和单元应力。有限单元法求解程序的内部过程如图1所示。

图1 有限单元法求解程序的内部过程

3 采空区稳定性的数值模拟

本文主要对开磷矿区3#采空区和6#采空区稳定性进行数值模拟，3#空区尺寸为120 m×50 m×17.5 m；6#采空区主要由A、B、C三个小采空区构成，三个空区尺寸分别为85 m×72 m×30 m,50 m×30 m×15 m、30 m×30 m×30 m，空区埋深均为230 m。

3.1 基本假定

假定矿体为理想弹塑性体，矿体和围岩为局部均质、各向同性的材料。考虑到采空区跨度极大，计算按平面应变问题进行。由于岩石具有脆性，分析中涉及到的所有物理量均与时间无关。

3.2 岩石物理力学参数

顶板岩性为白云岩，有轻微破碎，节理不太发育。顶板岩石密度为2.7 t/m3，抗拉强度为5.5 MPa; 单轴抗压强度为159 MPa，弹性模量为30，泊松比为0.27， 内聚力(抗剪强度)为37.49 MPa，内磨擦角为32.27°，岩石较稳定；矿石岩性为褐色磷矿石，较为致密坚硬。密度为3.3 t/m3， 抗拉强度为4.5 MPa，单轴抗压强度为148 MPa，弹性模量为29，泊松比为0.25， 内聚力(抗剪强度)为36.67 MPa，内磨擦角为41.94°，岩石较稳定；底板砂岩呈青灰色，有轻微破碎，节理不发育，其密度为2.7 t/m3，抗拉强度为3 MPa，单轴抗压强度为110 MPa，弹性模量为18，泊松比为0.23，内聚力(抗剪强度)为29.78 MPa，内磨擦角为42.56°，岩石较稳定。底板页岩呈酱红色，层理明显、节理发育，岩石易风化，遇水易泥化，其密度为2.7 t/m3，抗拉强度为2.7 MPa，单轴抗压强度为1 MPa，弹性模量为9.5，泊松比为0.39，内聚力(抗剪强度)为14.09 MPa，内磨擦角为42.56°，岩石不稳定。

3.3 建立计算模型

采用ANSYS前处理分别对3#、6#采空区建模。为尽可能反映采空区对围岩应力场分布的影响，所建模型的尺寸大小为5倍采空区尺寸。模型建立完成后进行网格划分，3#空区和6#空区模型划分网格后分别如图2、图3所示。

图2 3#采空区模型单元

图3 6#采空区模型单元

3.4 边界条件与荷载

施加模型边界条件主要考虑构造应力的作用。首先对3#、6#采空区模型X轴方向的前后两个面进行法向约束，然后再对Y轴方向的前后两个面进行约束，最后对两个采空区模型的底面进行约束。

本模型主要分析空区围岩的稳定性，因此主要考虑在重力场的作用下围岩的应力应变特征。

原岩自重应力可分别按下式进行计算：

(7)

(8)

式中：γi为上覆第层岩体重度，kN/m3；为上覆岩体泊松比；Hi为上覆岩体分层厚度，m。

3.5 模拟结果分析

3#空区和6#空区X、Y、Z向的应力分布，分别如图4、图5所示。

图4 3#采空区等效应力图

模拟结果表明，3#采空区的顶板主要受到35.5 MPa的压应力，并在空区中央出现拉应力集中现象，最大拉应力为13.6 MPa，远大于顶板岩体抗拉强度，空区中央岩体发生拉伸破坏。从空区中心向空区两端延伸，拉应力逐渐减小，直至变成压应力，并在空区顶板与侧帮相交位置出现压应力集中，最大压应力为84.7 MPa。空区底板与侧帮相交位置拉应力作用较为明显，最大拉应力为16.6 MPa，远大于底板岩体的抗拉强度，发生拉伸破坏。对于采空区的凸出部位主要表现为拉应力，最大拉应力为11.6 MPa，大于岩体抗拉强度，也发生拉伸破坏。而采空区的凹进部位压应力作用较为明显，最大压应力为84.7 MPa。

图5 6#采空区等效应力图

6#采空区由A、B、C三个采空区构成，在三个采空区上方岩体主要受到拉应力作用，最大拉应力为20.136 MPa，均大于岩体抗拉强度，说明6#采空区上部及四周已发生拉伸破坏，这与实际观测结果相符。从图5中还可看出，A采空区顶板围岩受到拉应力大于B、C空区，这是由于A采空区体积大、跨度大所致。

4 结 论

采用ANSYS有限元分析软件，分别对用沙坝矿3#、6#采空区周围岩体的稳定性进行数值分析，得到如下结论：

(1) 3#采空区和6#采空区上方岩体受拉应力作用明显，所受最大拉应力分别为13.581 MPa和20.136 MPa，大于岩石抗拉强度，说明3#、6#空区均不稳定，易发生拉伸破坏；

(2) 采空区附近岩体主要受拉伸破坏，6#采空区围岩所受拉应力大于3#采空区，说明采空区体积增大，采空区上部围岩所受拉应力相对增大。

参考文献：

[1]冯长根, 李俊平, 于文远, 等. 东桐峪金矿空场处理机制研究[J]. 黄金, 2002, 23(10): 11-15.

[2]于学馥, 郑颖人. 地下工程围岩稳定分析[M]. 北京: 煤炭工业出版社, 1983.

[3]唐胜利,唐 皓,郭 辉. 基于BP神经网络的空洞型采空区稳定性评价研究[J]. 西安科技大学学报, 2012, 32(2): 234-238.

[4]程爱宝, 王新民, 刘洪强. 灰色层次分析法在地下采空区稳定性评价中的应用[J]. 金属矿山, 2011(2): 17-21.

[5]魏海波, 吴 敏. 边坡的有限元分析及ANSYS软件对边坡开挖的模拟[J]. 云南水力发电, 2004, 20(4): 42-44.

[6]孙永刚. 有限元基础教学的探索[J]. 科技信息, 2010(24): 517-518.