基于三维数值分析的大桥磷矿采空区稳定性研究

王飞飞，邹 平，孟中华，李爱兵，刘正宇，虎万杰，马 增

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

矿山采空区的稳定性一直是矿山安全生产的关键性因素，关于矿山采空区稳定性问题，有众多学者做了相关研究[1-8]，并提出了处理采空区的一些方法与技术，例如：崩落法处理采空区、充填法处理采空区及封闭处理采空区。由于地下采空区具有隐伏性强、空间分布特征规律性差、采空区顶板冒落以及采空区塌陷情况难以预测等特点。因此，准确地分析采空区稳定性是关乎矿山能否安全生产的重要因素。

关于采空区的稳定性，大多学者主要研究的是影响采空区稳定的因素，对采空区矿柱在上部岩体作用下状态的研究尚少。本文采用MIDAS-GTS/NX有限元数值模拟软件，建立了采空区—矿柱—矿山三维力学分析模型，通过非线性静力分析得到了矿柱剪应力，矿柱轴向应力与矿山塑性区分布。研究结果可以为相似矿山提供参考。

1 矿山概况

1.1 矿山地质

大桥矿区磷矿矿体呈灰黑色薄层状，产于震旦系上统灯影组麦地坪段含磷白云岩中部，含矿层位岩性特征为：上、下部灰白色薄至中厚层状含磷砂屑白云岩-粉晶白云岩，中部为深灰—灰黑色胶磷矿磷块岩。矿层呈层状产出，层位稳定，矿体沿大渡河东岸陡坡地段分布，呈南北向出露，矿体长 4 250 m，厚度 0.42～1.02 m，平均厚 0.57 m，倾向 111°～194°，倾角 5°～10°。矿层底板为灰白色薄层至中厚层状含磷白云岩夹燧石条带，顶板为灰白色薄层至中厚层状含磷白云岩。矿层为深灰、灰黑色磷块岩，地表及深部矿层厚度变化规律不明显，矿层层位稳定。矿层围岩主要为含磷白云岩，与矿体界线清楚。

1.2 矿山采空区

大桥磷矿已经过三十多年的有组织无规划、采富弃贫等不同程度的开采作业，形成了数量巨大的、不规则的采空区。形成采空区面积约为45万m2，留下了许多不规则矿柱，部分已经变形，空区顶板出现少量冒落现象。经现场勘查，近些年新增采空区48万m3左右，但通过人工电耙配合，使用前期废石、废碴充填接顶，后阶段在荒石中搅拌入5%水泥通过人工电耙配合充填接顶，已治理约20万m3，目前尚存在大面积的采空区未治理。随着矿山开采工作的继续，很难保证原有的老采空区稳定性。因此，矿山老采空区亟待开展矿柱稳定性研究工作。

2 建立模型

2.1 非线性分析理论

非线性有限元解法是一种通过累积增量解收敛于正确解的方法，计算过程如图1所示。

图1 增量法Fig.1 Incremental method

图中tfext和t+Δtfext为t时刻和t+Δt时刻的外力，两个时刻的解和增量解关系如式(1)：

t+Δtμ=tμ+Δμ

(1)

式中：Δμ为时间增量Δt产生的增量解，累积增量解表达式如式(2)：

Δμi+1=Δμi+δμi+1

(2)

式中：Δμi为直到第i次迭代计算的累计增量解；δμi+1为第i+1次迭代计算时的增量解。δμi+1通过式(3)计算得到，Ki+1为切线刚度矩阵。

(3)

gi为残余力、非平衡力。非平衡力gi表达式如下：

gi=t+Δtfext-fint,j

(4)

式中t+Δtfext为外力，fint,j为内力。对方程(2～4)进行重复迭代计算，直到满足收敛准则为止，收敛准则包括力准则、位移准则和能量准则等。

为了改善上述的迭代求解问题，MIDAS -NX采用线性搜索功能来解决该问题，线性搜索算法(如图2所示)是在把增量δμi+1累计到增量解时采用一个比例值概念η来改善精度，此时累计增量解计算如式(5)：

Δμi+1=Δμi+ηδμi+1

(5)

假定Δμi+1满足平衡状态，利用总势能不变原理，可将线性搜索问题归结为解总势能对η导数为零时的η值问题，如式(6)所示：

(6)

图2 线性搜索算法Fig.2 Linear search algorithm

假定能量对η的导数是线性变化，则满足方程(6)的η计算见式(7)：

(7)

η=0或1时的斜率表达式如式(8)、(9)：

(8)

(9)

虽然线性搜索法不能十分准确地满足实际情况，按照式(8)、(9)计算s(η)一般不为零，但在MIDAS-NX中上面所述的过程会重复进行直到s(ηi)/s(η=0)的值低于设定值为止，即达到预定的计算精度。

2.2 采空区矿柱模型建立

采用MIDAS-GTS/NX有限元模拟软件建立采空区—矿柱—矿山三维力学分析模型，数值模型大小与原型一致。矿柱与岩体均采用实体单元模拟，采用摩尔—库伦屈服准则[9]，每个采空区的埋深均参照实际情况建立，矿柱埋深100～200 m，采空区高度1～2 m。采用四面体单元划分三维模型，矿柱与围岩采用实体单元模拟。计算模型划分节点总数有121 483个，单元总数有651 810个，模型为各项同性均质理想的线弹塑性岩土体。模型底部与四周采用固定边界，上部采用自由面边界。三维模型分为三层，上层为白云岩围岩，中部为含有采空区的矿体，下层为与上层同样的岩体。矿山的岩体力学参数如表1所示。鉴于矿山采空区较多，限于篇幅，本文仅对大桥磷矿的5-2、9-2与320三个采空区进行分析。所建立的三维矿山力学分析模型如图3所示。

表1 材料物理力学参数

图3 三维矿山数值模型Fig.3 Three-dimensional mine numerical model

3 数值模拟结果与分析

对建立的不同采空区三维数值模型进行非线性静力分析，得到了矿柱剪应力云图，矿柱轴应力云图与矿山塑性区分布图。

3.1 矿柱剪应力

矿柱剪应力云图如图4所示。

在图4(a)中，图例中从上到下剪应力由高到低，5-2采空区矿柱最大剪应力为34.9 MPa，最小剪应力为10 MPa。矿山5-2采空区埋深100～200 m。由图可知不同矿柱的剪应力大小差异较大。有的较大横截面的矿柱剪应力较大，而有的较小的横截面的矿柱剪应力较小，可以得到矿柱剪应力大小与矿柱横截面的大小没有直接的关系。剪应力较大或较小的矿柱集中在一个区域内，说明矿柱剪应力的大小与上覆岩层厚度有直接联系，较厚的上覆岩层的重力较大导致矿柱承受的剪切力较大。最大剪力产生在小矿柱中间位置处，这与现场矿柱的破坏形式非常吻合，验证了数值模拟计算得到的结果是可信的。

由图4(b)可知，9-2采空区矿柱最大剪应力为39.9 MPa，最小剪应力为14 MPa。矿山9-2采空区埋深260～360 m，与5-2采空区相比，其埋深增加了160 m左右，故9-2采空区矿柱的剪应力较大。9-2采空区中，有个别矿柱的剪应力较大，而周边的矿柱剪应力较小，可能是因为矿柱分布导致上覆压力对个别矿柱产生集中力，最终极个别的矿柱剪应力较大。9-2采空区矿柱剪应力在整体采空区分布较为均匀，主要是因为该采空区上部山体高低起伏度不大，对该采空区产生较为均匀的压力。

由图4(c)可知，320采空区矿柱最大剪应力为38.9 MPa，最小剪应力为14.8 MPa。矿山320采空区最大埋深在150～230 m范围内，其埋深比9-2采空区小，比5-2采空区大，故320采空区矿柱最大剪应力比9-2采空区矿柱剪应力小，比5-2采空区矿柱剪应力大。

通过比较图4中三个采空区矿柱剪应力大小与分布可知，采空区剪应力大小顺序为：9-2采空区>320采空区>5-2采空区，相比较5-2采空区矿柱较为稳定。采空区矿柱剪应力与矿柱的埋深有直接的联系。

3.2 矿柱轴向应力

采空区矿柱轴向应力云图如图5所示。

由图5(a)可知，5-2采空区矿柱的最大轴向应力为84.44 MPa，最小轴向应力为22.93 MPa。由图5(b)可知，9-2采空区矿柱的最大轴力为100.75 MPa，最小轴力为39.56 MPa。由图5(c)可知，320采空区矿柱的最大轴力为98.78 MPa，最小轴力为42.19 MPa。通过比较以上三个采空区矿柱的轴力，可知5-2采空区较为稳定，其它两个采空区具有不安全可能性，因此，要加强监测与防护措施。

3.3 矿山塑性区

矿山塑性区分布如图6所示。

图4 采空区矿柱剪应力Fig.4 Shear stress of pillar in goaf

图5 矿柱轴向应力云图Fig.5 Axial stress cloud map of the pillar

图6 矿山塑性区分布Fig.6 Distribution map of plastic zone in mine

大桥磷矿矿体赋存倾角5°～10°，矿体上盘有向山体内侧移动的趋势，不会向山体外侧移动，即不会发生山体垮塌堵塞大渡河的危险。由图6可知，矿山没有产生大面积的塑性区分布，仅在矿体处有局部塑性区。整体而言，山体稳定性良好。

4 结论

本文采用MIDAS-GTS/NX有限元数值模拟软件，建立了采空区—矿柱—矿山三维力学分析模型，采用非线性静力分析方法对大桥磷矿采空区开展了三维数值模拟研究，得到以下结论：

1)矿柱剪应力大小与矿柱横截面的大小没有直接的关系，矿柱剪应力的大小与上覆岩层厚度有直接联系。采空区剪应力大小顺序为：9-2采空区>320采空区>5-2采空区，相比较5-2采空区矿柱较为稳定。

2)采空区埋深对采空区矿柱的剪应力与矿柱轴向应力均有很大的影响，上部岩体的重力较小，矿柱产生的剪应力与轴向应力较小。

3)矿山没有产生大面积的塑性区分布，仅在矿体处有局部塑性区。整体而言，山体稳定性良好。