采空区三维建模及稳定性数值模拟研究
2018-12-05容宇
容 宇
(云南驰宏锌锗股份有限公司 会泽矿业分公司, 云南 曲靖市 654211)
0 引 言
采矿区顶板及围岩的稳定性是金矿矿床高效开采面临的主要问题。在深部开采过程中事故大多有多发性和突发性特征,原因是矿山深部岩体始终处于高地应力、高渗透压、高温、爆破扰动,即“三高一扰动”的不良地质条件。目前,某铜矿进行了部分矿房的回采,而充填工作刚开始,因此存留了大量的采空区,加重了地应力对回采作业的干扰。采矿区存在的危害主要变现为与采空区毗邻的穿脉巷道发生非均匀变形破坏、片帮或冒顶;地应力分布不均,采场部分区域应力集中不易识别和控制;围岩蓄能较大,岩爆风险加大等等。正是因为这些采空区的存在,还可能造成更多灾难性事故的发生。采空区稳定分析对于井下应力控制、指导生产具有重要作用。数值模拟方法是最为常用的采空区稳定分析方法,其可以提供较为直观的定量描述,识别弹塑性区域,对于指定合理的生产方案十分有益。
本次数值模拟研究采用ANSYS三维有限元数值模拟软件建立采空区、地层、断层和巷道单元模型,划分网格后再导入FLAC3D有限差分数值模拟软件中进行计算。通过定量计算,分析空区开挖导致-120 m中段及各纵剖面的沉降(Z方向的位移)和弹塑性区域分布。结合采空区形成后在-120 m阶段所引起的沉降变化,以及提取相应方向的塑性区结果,为后续采空区周边的巷道开挖提供参数指导。
1 模型前处理
1.1 模型建立
该矿模型以井下实际情况为基础,选取既能满足计算要求且范围最小的计算模型,同时考虑对应力状态影响较大的地层、断层的影响,忽略地下水及爆破等外在影响因素,揭示该矿采空区对上部地层的影响规律。
模型范围为:X坐标范围为49600~50323.2684;Y坐标范围为8650~9050。确定研究范围后即对模型进行建立,模型分外部围岩、内部断层影响带、内部采空区、内部-120 m水平开拓巷道。
(1)地层模型。结合矿山所提供的地形地貌图等CAD图件,通过地表范围划设地表点网格,然后提取各网格点坐标,建立地表曲面。地层建成后,依据采空区范围及计算要求选取底部水平标高为-240 m,形成含地层的围岩体。
(2)断层模型。通过地质图和采空区调查结果得到断层附存参数,包括位置坐标和倾向倾角等(见表1),对断层进行建模,其中,F25走向NNW~NEE(近EW),F26走向NEE~EW。
表1 F25、F26断层参数
(3)采空区模型。采空区包括测量采空区和推测采空区,对于已经测量过的采空区可通过三维建模方法直接构建,如图1所示。对于受条件限制无法测量的采空区和已经充填的采空区,主要是依据上下穿脉巷道的位置和该矿纵剖面图确定其宽度和高度,从而大致对这些采空区进行建模,其结果如图2所示,其中有L1、L8、L14、L16、L18采场。地层、断层、采空区模型构建后将其组合,如图3所示。
图1 已测采空区建立模型
图2 未测采空区模型
图3 某铜矿所有采空区模型
(4)巷道模型的建立。巷道模型依据矿方提供的-120 m实测图,建立上下盘主运输巷道及该矿区域内穿脉巷道和通风巷道模型,建立的巷道和采空区空间三维图如图4所示。
图4 巷道及采空区三维模型
1.2 网格划分
对上述建立的模型进行有限元计算,对该矿模型通过限制单元尺寸的方法对模型进行均匀划分,使模型网格划分均匀精细。最终划分的网格模型如图5所示。模型共划分226748个节点和1303313个单元,满足计算要求。
图5 该矿网格模型
1.3 参数选择
计算过程中将地层、矿体及断层均视为单一连续的介质,为保证模拟结果的准确性,根据岩层的岩性对个部分逐一设定物理力学参数,试验参数来源于原有试验报告,具体参数如表2所示,充填体的体积模量为0.25 GPa,剪切模量为0.08 GPa。
表2 计算参数
2 结果分析
2.1 初始应力结果分析
取粘聚力和抗拉强度,采用更改参数的弹塑性求解初始应力的方法,对该矿采空区的初始位移及应力情况进行计算,其Z方向位移计算结果见图6。结果表明,在自重应力的作用下,Z方向的位移呈分层分布,最大沉降位移出现在地表,为-1.3 m。
图6 初始应力Z方向位移
2.2 开挖结果分析
初始应力计算完成后,对各方向的位移、速度及塑性区清零,改为实际参数对该矿进行采空区开挖模拟,计算其开挖后的位移状态,开挖后的Z方向位移表面整体图如图7所示,最大位移出现在被断层切割的岩体中,为7 cm。
图7 开挖后Z方向位移三维图
2.2.1 地层位移模拟结果
图8为-120 m平面的Z方向位移云图,其中显示的范围主要为采空区正上方区域。采场开挖后,-120 m平面的最大沉降位移为9.732 mm。
图8 采空区开挖后-120 m平面Z方向位移图
为了解模型内部位移及塑性的变化规律,对模型进行切片分析,分别截取-120 m平面、垂直巷道和经过采空区的纵剖面,纵剖面选择结果见图9所示。同时,提取3个纵剖面的Z方向位移图,结果如图10所示。
图9 纵剖面选取位置示意
图10 Z方向位移图
由图10可知,采空区范围内对应的Z方向最大位移分别为2.95,4.248,4.499 cm,最大位移发生在暴露面积较大的采空区顶板处。
2.2.2 塑性区域模拟结果
为了解模型的塑性破坏情况,提取了地表面、-120 m平面和3个纵剖面的塑性区图,模拟结果见图11。由图11可知,-120 m平面采空区范围内未进入塑性破坏,只有断层经过的区域有塑性破坏;3个纵剖面在采空区附近均有部分区域进入塑性破坏状态。虽然采空区上部并未形成联通的塑性区,但由于采空区L1、L2、L3体量相对更大,塑性区域更为明显,上部岩层需要开挖巷道或回采时,要远离该区域。
图11 地层塑性区图
3 结 论
采空区稳定分析对于井下应力控制、指导生产具有重要作用。数值模拟方法是最为常用的稳定分析方法,本研究采用ANSYS软件建立采空区、地层、断层和巷道单元模型,划分网格再由FLAC3D软件进行计算,分析空区开挖导致的-120 m中段及各纵剖面的沉降(Z方向的位移)和弹塑性区域分布,研究结果表明:
(1)采空区形成后,-120 m平面,最大Z方向位移(沉降位移)为0.973 cm;1#、2#、3#纵剖面,对应的最大Z方向位移分别为2.95,4.248,4.499 cm。
(2)计算得到了各截面的塑性区结果,采空区模拟开挖工况下的塑性区均未联通,对-120 m水平的塑性区范围影响有限。
(3)结合采空区形成后在-120 m阶段所引起的沉降变化,以及提取相应方向的塑性区结果,为后续采空区巷道开挖提供指导。