基于FLAC3D 对开采引起地表沉降规律研究

王志凯 ,宋文龙

(1.五矿有色金属股份有限公司,北京 100043;2.鑫联环保科技股份有限公司,北京 100027)

0 引言

地下开采造成了岩土体内部原有应力平衡破坏,往往导致井下采场围岩冒落、变形、移动甚至断裂[1-2]。围岩的破坏和移动有时会波及地表,致使地表发生不同程度的变形和破坏,进而严重影响井下人员的作业安全[3]。因此,如何利用一定的理论和数值方法对矿山开采引起地表沉降规律进行深入分析研究是保障矿区安全和人员作业安全具有重要意义。

目前,在我国主要采用经验公式法、工程类比法、现场监测法等方法对浅部矿产开采引起地表沉降影响范围和影响程度的研究[4]。但是随着开采深度的不断加深,深部开采涉及“三高一扰动”问题,同时井下环境恶劣,上述方法常带来各种弊端。随着新技术、新理论、新方法的发展,尤其计算技术的应用,对井下开采造成地表移动与变形影响研究取得了许多进展。刘宝琛等[5]采用随机介质理论分析了开采中岩石移动问题。李增琪等[6]通过层间递推公式将多层体矿压和岩层移动的计算得到了岩层移动的三维层体模型。唐春安等[7]研究开发RFPA2D能实现开采过程中覆岩变形、离层、断裂及地表沉陷规律模拟。

为此,本文基于FLAC3D运用有限差分法从地表垂直位移、水平位移、斜率与曲率三个维度对某金矿开采引起地表沉降问题进行数值分析,探究其引起地表沉降的规律,以期为该矿安全开采提供一定保障。

1 工程背景

某金矿已开采至-730 m 中段,开拓深度已至-1 030 m 水平,深部矿体倾斜深至-1 500 m 以下。Ⅰ、Ⅴ号矿体是该金矿主矿体,两者占矿山资源总量的96.42%。Ⅰ号矿体上部走向长约240 m,下部300 m,平均厚度30 m,倾角平均30°。先后开拓了-10～-580 m 的17 个中段,中段高度20～50 m,-10 m 以上标高作为永久保安矿柱。

Ⅴ号矿体为中间厚上下两边薄的透镜状矿体,平均厚度80 m,上部走向长约240 m,下部走向300 m,倾角30°。先后开拓了-480～-730 m 六个阶段,阶段高度50 m。

根据该矿的实际开采情况,确定研究对象为Ⅰ号和Ⅴ号整个矿体,模型简化Ⅰ号矿体赋存高度-10～-580 m,走向长300 m,平均厚度30 m,倾角30°。Ⅴ号矿体赋存高度-480～-730 m,走向长240 m,厚度80 m,倾角30°。

2 矿体的模型建立与参数的选取

2.1 矿体模型的建立

根据矿体赋存条件计算方案中计算区域取1 800 m×1 500 m×800 m。模型采用岩土工程数值模拟软件建立模型及划分网格,共划分802 872 单元,831 480 个节点,计算模型示意图如图1 所示。

图1 计算模型示意图

2.2 模拟参数的选取

2.2.1 参数与本构模型的选取

根据工程实际情况,选定模拟计算采用的岩体力学参数如表1 所示,数值计算中主要采用摩尔-库伦本构模型。

表1 岩体力学参数

2.2.2 边界条件

边界应力与约束条件是计算模型的重要内容,直接影响计算结果的可靠性及精度。应用应力解除法测量某金矿地应力得出的最大水平主应力、最小水平主应力、垂直主应力关系式:

式中σh.max—最大水平主应力;

σh.min—最小水平主应力;

σv—垂直主应力,MPa;

H—测点埋深,m。

采用基于应力边界的快速地应力场拟合方法,应用此法形成的初始应力场如图2 所示,与实测地应力场有较高的吻合度。

图2 初始应力场云图

3 开采过程地表沉降规律分析

3.1 地表竖直位移特征及分析

(1)地表垂直位移特征分析

图3、4 分别为I #矿体和V #矿体-730 阶段开采结束地表沉降位等直线图。从图中可以看出,I #矿体开采结束后,地表沉陷位置比较集中,地表最大垂直位移达到了148 mm,影响范围沿走向方向约900 m,垂直于走向方向约1 100 m。最大下沉量出现在矿山的上盘位置,约为矿体走向方向的中心位置。随着开采向深部进行,沉降中心逐渐向V #矿体上部移动,V #矿体-730 阶段开采结束后地表下沉幅度和影响范围进一步增大;地表沉陷中心略向V#矿体上方偏移,沉陷盆地形状基本没有什么变化,最大下沉量达到166.5 mm,增大幅度达到18.5 mm,影响范围沿走向方向约900 m、垂直走向方向约1 300 m。表明地表沉降主要由于浅部资源开采引起,深部开采对地表沉降的影响不大。通过等高线图还可以很清楚地看到地表沉降形成凹陷区域的平面图,随V #矿体开采深度的增加,地表沉降逐渐扩大,并且沉降中心逐渐向V #矿体上方移动。等值线在垂直于走向方向由浅部开采到深部开采过程中逐渐变疏,A 较B 位置较密,可知在A 位置地表倾斜变化率较大,此位置地标建筑物更容易产生破坏。

图3 I#矿体开采结束后地表沉降等值线图

图4 V#矿体-730 阶段开采结束地表沉降等直线图

3.2 地表水平位移特征及分析

通过图5、图6 可以看出,随着各矿体相续开采,地表水平移动范围不断扩大,该金矿最大水平主应力垂直于走向方向,形成了以采空区为中心,并且地表水平位移向空区两侧不断减小的水平移动带,水平移动范围基本呈现对称分布。I #矿体开采结束水平移动最大值约为170 mm,V #矿体-730 阶段开采结束后正向移动最大值约为180 mm,增加10 mm,幅度不大。A～C 区间水平移动值逐渐减小,以拉伸变形为主,A 区域水平移动值最大,A～B 区域逐渐减小,表明在拉伸变形到达一定平衡时转而以压缩变形为主。B 区域等值线较密表明此处地表建筑受压缩变形破坏影响较大。

图5 I#矿体开采结束后地表x 水平位移等值线图

图6 V#矿体-730 阶段开采结束地表x 水平位移等值线图

3.3 地表斜率与曲率的数值微分计算

通过在地表坐标设置历史记录点,将数值计算后地表的垂直位移和水平位移自动记录并保存,才能进一步计算地表变形曲率K、倾斜i、水平变形值ε。记录点的间距跟矿体的平开采深度呈正相关。开采深度越大,记录点的间距也越大。在开采深度大于500 m 时,选取的点间距应为50 m,所以在地表每50 m 设置一个记录点,记录水平和垂直位移值。通过仿真软件计算结束后生成的x、z两个方向位移值,得出的地表所有点的结果如表2 所示。

表2 地表监测点的水平和竖直位移值

数值微分问题指已知函数y=f(x)的一组离散值(xk,f(xk))(k=1,2,…,n),求f(x)在节点处的导数或微分值。对离散点求导数,使用数值微分方法比较方便。

根据向前差商或向后差商,分别得到微分方程:

对式(4)(5)取算术平均值,并用二阶中心差商近似f″(x0),便得到f″(x0)的数值微分公式(6)。

通过求式(6)误差渐进展开式。根据Taylor 公式可得:

令G(h)为上面两式等号左端的第二项,对于一阶和二阶导数可建立外推公式:

计算时,步长分别取100 m 和50 m,即可按照所得的Gm+1(h)值计算出相应的一阶导数和二阶导数,即地表变形的倾斜和曲率。

计算后所得地表最大倾斜0.6 mm/m,最大曲率0.005 2 mm/m2,最大水平变形0.404 mm/m。对比《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程》中建筑物保护等级及允许变形值分类。整个地表出现的变形值都没有超过一般砖木结构允许的临界变形值。因此,充填开采后地表的变形不会对地表建(构)筑物的安全造成较大的危害。

4 结论

(1)通过FLAC3D对地表竖直位移分析得到I#矿体开采结束后,地表沉陷位置比较集中,沉陷盆地形状较规整,地表最大下沉量达到了144.9 mm,影响范围沿走向方向约900 m,垂直于走向方向约1 100 m。V#矿体-730 阶段开采结束后,最大下沉量达到166.5 mm,增大幅度达到21.6 mm,影响范围沿走向方向约900 m、垂直走向方向约1 300 m。

(2)从地表水平位移得到,I #矿体开采结束正向移动最大值约162 mm,V #矿体-730 阶段开采结束后正向移动最大值约175 mm。

(3)通过地表斜率与曲率的数值微分计算得到,地表最大倾斜0.6 mm/m,最大曲率0.005 2 mm/m2,最大水平变形0.404 mm/m,均没有超过一般砖木结构允许的临界变形值。因此,开采后对采空区进行及时充填,地表的变形不会对地表建(构)筑物的安全造成较大的危害。