缓倾斜厚大矿体崩落法开采对临近构筑物的影响

徐 曲，王玉平，曹华兵，杜建华，冯 杰

(1.湖北省应急救援中心，湖北 武汉 430070；2.武汉盛世金安安全科技有限公司，湖北 武汉430080；3.武钢资源集团程潮矿业有限公司，湖北 鄂州 436050；4.中南财经政法大学信息与安全工程学院，湖北 武汉 430073)

地下矿产资源采出诱发的岩层沉陷、地表开裂是采掘区域岩体初始应力平衡被打破后，应力的转移、释放、再平衡的过程，随开采范围、深度加大围岩发生渐进式破坏、变形，由采掘区域向地表逐渐发展，进而导致地表、井巷工程、地表构筑物的破坏，其发展规律受矿体空间形态、采掘工艺、工程与水文地质条件等多种因素的影响。岩层移动、破坏从广义上来讲属于开采沉陷学的范畴，按照矿业行业划分，煤矿开采领域相关的研究成果较丰富，理论性强，这主要是由于煤矿多属于沉积矿层，矿体、岩层多呈层状分布，岩层变形、破坏规律性较为一致。而金属矿床多在地壳活动过程中岩浆沿地壳断裂、薄弱地带的侵入、喷出形成，其具有空间形态多变、围岩种类繁多、强度多变、分布无序、岩体各级结构面错综复杂等特点。因此，概率积分法、砌体梁经典理论、关键层理论等方法均不能很好地应用于金属矿山开采引起的地表沉陷、岩层破坏的预测。

金属矿山开采引起的岩层移动、破坏规律是矿岩体空间分布、强度、结构面特征以及采掘工艺等一系列因素综合影响的结果，对其研究是巨大的系统工程问题。由于金属矿山开采引起的地表沉陷问题的复杂性，导致我国早先矿山设计在岩层移动范围的测算上出现误差，进而导致了众多矿山关键工程破坏，例如金川集团二矿区14行风井破坏事故,程潮铁矿地表运矿巷道、措施井开裂等。

针对金属矿山开采过程中地表沉陷问题的研究方法主要包括经典理论分析、数值模拟、物理模型试验等方法，其中三维数值模拟计算方法对于适合于复杂空间形态多因素协同作用诱发的地表沉陷、岩层破坏问题的研究最为有效。因此，本文以程潮铁矿无底柱分段崩落法开采为背景，通过构建涵盖矿体、围岩体分布、采掘工艺、初始应力条件等信息的高精度三维空间模型，并通过准确的岩体力学实验，选取适合的本构模型，且辅以地表沉降监测措施，可较为完善地解决复杂条件金属矿山开采诱发的地表沉陷的测算问题。

1 数值模型构建

1.1 模型结构

本研究采用三维有限差分软件(FLAC3D)，根据矿岩体空间分布特征、强度特征、采掘顺序建立数值计算模型，按照时间顺序真实地模拟矿体开采，从而得出研究矿区内的岩层移动和破坏规律。

程潮铁矿由多个分矿体组成，分矿体呈叠瓦状分布，多呈透镜状、扁豆状、脉状、囊状分布，矿体厚度变异大、夹岩复杂，矿体整体走向为北西西，倾向为南西，倾角为20°～30°，最大分矿体长度超过1 000 m，厚度超过100 m。

程潮铁矿实际开采过程中分为东西两个作业区域，西区矿体埋深深、厚度较大，东区矿体埋深浅、厚度相对较小。本文研究的地表构筑物位于西区，受东区矿体开采的影响较小，因此为了提升计算效率，本次建立的三维数值模型仅考虑了西区矿体和围岩，主要针对西区开采对地表临近构筑物的影响进行研究。根据地质勘探数据生成的程潮铁矿区矿体空间分布模型，见图1。

图1 程潮铁矿区矿体实体空间分布模型图

根据西区矿体赋存位置以及地表投影范围，确定矿体覆盖区域尺寸如下：走向方向为1 200 m，倾向方向为800 m，矿体最大埋深接近1 200 m。本文中仅对-500 m水平以上矿体开采进行研究，故根据-500 m水平以上矿体的三维几何外轮廓，按照移动角40°(保守取值)计算模型边界。数值模型取矿体走向为

X

轴，倾向为

Y

轴，垂直方向为

Z

轴，最终确定数值模型三维尺寸为：

X

方向3 000 m，

Y

方向2 800 m，

Z

方向1 200 m，见图2。

图2 程潮铁矿区数值计算模型示意图

数值模型单元尺寸的划分以及与地质模型、开采工艺的对应关系分析如下。

1.1.1 单元尺寸划分

(1) 精细网格划分区域：该区域主要为矿石赋存区域，其形状为平行六面体，依据矿体

X

方向覆盖范围确定左、右边界，根据矿体

Z

方向覆盖范围确定上、下边界，再根据矿体倾角与矿体在

Y

方向的投影边界确定前、后边界。该区域

Z

方向最小单元尺寸根据矿山实际生产中的开采分段高度(17.5 m)确定，矿体区域采用六面体单元划分网格，单元网格

X

、

Y

、

Z

三向尺寸比例控制在1∶2范围内，即六面体单元最长边长度不大于最短边长度的2倍，单元网格

X

、

Y

、

Z

三向尺寸依次为33 m、28 m、17.5 m。(2) 外部围岩区域：整个模型扣除上述精细网格划分区域之外均为外部围岩区域。该区域按照六面体单元划分网格，为了提高计算速度，单元网格尺寸适当放大，单元网格

X

、

Y

、

Z

三向尺寸依次为50 m、50 m、17.5 m，模型共划分357 840个单元、378 288个结点。

1.1.2 数值模型与地质模型的对应关系

数值模型建立后，首先导出每个单元的节点三维坐标，进而求取单元格质心坐标；然后将每个单元格的质心坐标导入地质模型，从而获得该区域矿岩类型，并以此为依据进行单元分组及属性赋值。

1.1.3 开采工艺模拟

依据开采顺序，由上至下模拟矿体开采，每次开采整个分段高度(17.5 m)的矿体，矿石采出后顶部围岩冒落区域按照椭球体计算，冒落区域最大宽度

l

与开采范围的关系依据下面公式确定：

L

<50 m时

l

=(1

.

1～1

.

5)

L

(1)

L

>50 m时

l

=(1

.

5～2

.

0)

L

(2)

式中：

L

为矿体长度(m);

l

为冒落区域最大宽度(m)。

1.2 计算参数

数值模型建模过程中已将模型边界扩展至开采影响区域之外，所以计算采取位移边界条件，即模型四周限制水平位移，底部限制垂直位移；模型顶面为地表，无需施加边界条件。

本构模型选取摩尔-库伦(Mohr-Coulomb)强度准则，岩体力学参数在岩体力学试验的基础上采用Hoek-Brown准则进行折减，模拟计算采用的岩体力学和变形参数见表1。

表1 程潮铁矿西区岩体的力学和变形参数

冒落区域为散体，其强度属性随时间发生变化，其强度计算公式如下：

ρ

=1 600+800(1-e-1.25)

(3)

E

=15+200(1-e-1.25)

(4)

υ

=0.05+(1-e-1.25)

(5)

式中：

ρ

为材料的密度(kg/m)；

E

为材料的弹性模量(MPa)；

υ

为材料的泊松比；

t

为时间(a)。

程潮北断层是区域稳定的重要影响因素，其力学参数见表2。

表2 程潮北断层的力学参数

模型应力条件根据矿区实测施加，第一主应力

σ

位于水平面内，方向与矿体走向一致；最小主应力

σ

位于水平面内，方向垂直矿体走向；第二主应力

σ

位于水平面垂直方向。埋深

H

位置的主应力分量计算公式如下：

(6)

式中：

γ

为岩层的平均容重(kN/m)，取值为28 kN/m；

H

为埋深(m)。

2 模拟计算结果与分析

2.1 模拟计算流程

针对程潮铁矿矿体开采诱发的岩层移动，已有研究表明，矿区地表变形和地表构筑物的破坏规律与岩体力学性质、矿体空间形态、采空区范围、开采顺序密切相关，因此在计算过程中，依据实际开采范围、开采顺序、无底柱分段崩落法采矿作业流程对矿体开采诱发的地表变形和地表构筑物的破坏规律进行了数值模拟计算。

根据矿区岩体强度和初始应力条件拟合开采前的原岩应力状态，在此基础上按照时间顺序模拟不同阶段的矿石开采。首采水平为-290 m，本研究的最大开采深度至-500 m水平，在此区间内按照17.5 m标高设置开采分段，共包括13个开采阶段。

2.2 地表错动区域范围判别准则

依据前冶金部《冶金矿山测量规范》规定的标准，将竖向倾斜率

S

=2.5 mm/m、水平拉伸率

T

=1.5 mm/m作为地表临界变形值用来判定地表错动区域范围。

2.3 模拟计算结果分析

本文采用开源GIS软件Cesium将模拟计算得到的

X

方向地表水平位移分量和

Y

方向地表水平位移分量合成为地表水平位移，进而根据下面公式(7)、(8)计算地表变形诱发的竖向倾斜率

S

和水平拉伸率

T

。假设点1、点2两点距离为

D

，点1的地表竖向沉降量为

D

，点2的地表竖向沉降量为

D

，则地表变形诱发的竖向倾斜率

S

的计算公式如下：

(7)

点1的地表水平位移为

D

，点2的地表水平位移为

D

，则地表变形诱发的水平拉伸率

T

的计算公式如下：

(8)

限于篇幅，本文仅对关键开采阶段的数据进行分析论述，-360 m水平矿体开采完毕后，该矿区地表变形的模拟计算结果如下：

X

方向为水平面矿体水平走向方向，

Y

方向为水平面矿体垂直走向方向，

Z

方向为竖直方向；矿体开采至-360 m水平时地表

X

、

Y

方向位移云图分别见图3(a)和图3(b)，将

X

、

Y

方向的地表水平位移分量进行矢量合成，得到地表水平位移合成见图3(c)，图3(d)为地表竖向沉降等值线云图。

图3 程潮铁矿区地表变形云图

通过提取不同开采阶段措施井所在位置的水平位移、垂直位移，可根据公式(7)、(8)计算得到对应位置的地表变形诱发的竖向倾斜率

S

和水平拉伸率

T

。

程潮铁矿区措施井所在位置各开采阶段地表变形数据汇总于表3，各开采阶段地表变形曲线见图4。

表3 程潮铁矿区措施井对应位置各开采阶段地表变形统计表

图4 程潮铁矿区措施井对应位置各开采阶段地表变形曲线图

由表3和图4可以看出：

(1) 按照开采中段计算，-360 m水平矿体开采完成后，措施井对应位置地表变形诱发的竖向倾斜率为0.80 mm/m、水平拉伸率为0.78 mm/m，上述两值均小于规范要求的限值，即说明开采至-360 m水平，措施井处于安全状态。

(2) -430 m水平矿体开采后，措施井对应位置地表变形诱发的竖向倾斜率为1.18 mm/m、水平拉伸率为1.50 mm/m，水平拉伸率达到临界值，即说明矿体开采至-430 m水平时，地表错动区域边界达到措施井附近。

(3) 但由于措施井为贯穿至地下采矿作业区域的圆筒结构，而岩层变形是由下部逐渐向上部围岩扩展，即由开采区域至地表的变形呈递减规律，在地表变形达到临近破坏限值前，措施井井筒地下部分变形已经超过限值。通过对措施井井筒地下部位不同开采阶段的地表变形分析可知，-360 m至-375 m水平矿体开采后，措施井井筒在地表200 m以下区域岩层水平拉伸率超过了1.5 mm/m限值，措施井井筒破坏区域首先在地表200 m以下范围产生，且随着开采深度的增加，其破坏程度加剧，破坏范围向措施井井筒上部扩展。

(4) -500 m水平矿体开采后,措施井对应位置地表变形诱发的竖向倾斜率为1.38 mm/m、水平拉伸率为1.60 mm/m，水平拉伸率超过了临界值，即说明开采至-500 m水平，措施井进入地表错动区域，其破坏区域严重。

3 结 论

本文采用三维仿真计算方法对程潮铁矿区内由于矿体地下开采引起的地表变形规律进行了研究，并对不同开采水平地表范围和典型剖面地表移动角进行了预测，得到的结论如下：

(1) -342 m水平以上矿体开采对于措施井所在位置地表变形的影响甚微，对措施井稳定性无影响。

(2) 由于-360 m至-412 m水平矿体向下盘方向扩展，-360 m水平以下矿体开采对措施井所在位置地表变形的影响加剧。开采至-375 m水平时，措施井所在位置地表变形诱发的水平拉伸率为0.9 mm/m，尚未超过允许极限值，但由于措施井为贯通至地下采矿作业区域的井筒结构，根据模拟计算结果得出：-360 m至-375 m水平矿体开采后措施井井筒在地表200 m以下区域岩层水平拉伸率超过1.5 mm/m限值，措施井井筒破坏区域首先在地表200 m以下范围产生，且随着开采深度的增加，其破坏程度加剧，破坏范围向措施井井筒上部扩展。

(3) 开采至-430 m水平时，地表变形诱发的水平拉伸率达到1.5 mm/m，措施井进入地表错动区域范围，措施井井筒破坏范围扩展至地面，由于-430 m至-500 m水平矿体宽度变窄，此范围矿体开采对于措施井稳定性的影响逐渐减弱。

(4) 措施井所在位置的地表错动角在开采至-500 m水平时为50°左右，远低于设计值(约60°)。