基于FLAC3D金川二矿区双中段充填法开采回采顺序研究

赵 杨， 李向东， 朱远乐， 王 涛*

(1.武汉大学水资源与水电工程国家重点实验室， 武汉 430072；2.长沙矿山研究院有限责任公司， 长沙 410012)

金川矿区采用下向进路胶结充填采矿法.充填采矿法于20世纪30年代起源于加拿大，经过近一个世纪的发展，已经成为一种高强度、高效率、高回采率的采矿方法.矿山深部开采，特别是高地应力、大范围开采下，采场失稳风险高、潜在风险高[1].金川矿区深部开采采用蜂窝式结构驱动技术实现了大面积连续开采，提高了矿区生产效率，金川矿区充填法采矿技术处于世界领先水平[2].在深部开采中，充填采矿法有利于地压平稳释放及减少岩爆的发生，与其他采矿技术相比具有显著优势，在矿山开采工程中得到了广泛认同[3].充填采矿法在20世纪60、70年代被引入我国，广泛应用于有色金属矿山开采，获得了长足发展.陶干强等[4]结合国内某锡矿山，以有限元法为手段进行充填法采矿中采场结构参数优化，对开采方案进行综合评价.于世波等[5]以某铁矿工程为研究对象，将现场监测数据与数值模拟结果对比，研究充填体对围岩变形控制作用.叶义成等[6]以上横山矿为例，构建物理相似模型，研究回采顺序对巷道变形的影响.尚振华等[7]运用有限差分方法对我国某大型铜矿多中段回采顺序优化展开研究.赵千里等[8]采用有限差分方法对金川矿区竖井开挖展开研究，研究发现深部衬管应力变形大，证明了释放位移法施工技术有助于减小衬管变形.江权等[9]提出了运用激光扫描技术的变形监测方法，并运用于金川二矿区巷道变形监测中，根据监测资料研究开采支护合理时机.回采顺序优化是充填采矿法研究的重点问题之一，数值模拟是研究回采顺序优化的重要手段.

为了进一步提高采矿效率，诸多矿区开辟多个工作面，采用双中段以及多中段的开采方案.多中段开采过程中会产生大面积水平矿柱，水平矿柱的安全及回采问题成为矿区生产安全问题的关键影响因素.黄昌兴等[10]针对彝良某铅锌矿多中段水平矿柱的回采问题展开研究，提出六边形-矩形联合进路布置方式来提高水平矿柱稳定性.程龙[11]以江西省宜丰新庄铜铅锌矿区为例，展开渗流场和应力场共同作用下防水矿柱留设宽度研究.吴顺川等[12]以甘肃某矿为例，用有限差分法研究水平矿柱回采过程中矿柱力学响应及沉降规律.水平矿柱稳定性是多中段开采稳定性研究的重点.

金川二矿区采用大面积不留矿柱充填法同时开采两中段时，遇到了诸多问题.该矿区矿体埋藏较深，地应力较高，对开采稳定造成严重不良影响[13].双中段开采中，两中段之间会产生水平矿柱，且随着开采的进行，水平矿柱逐渐变薄，刚度减小，承载力降低，水平矿柱的“抵顶作用”会造成开采巷道出现严重变形[14].金川矿区分期开采，致使前期开采留下大体积充填体，矿山进入深部开采后，上覆大体积充填体的稳定性是影响矿山生产安全的重要因素[15].

针对金川二矿区开采形成的水平矿柱、大体积充填体失稳风险及高地应力下双中段同时开采围岩稳定性问题，本文以有限差分数值模拟的手段对回采顺序对开采稳定性的影响进行研究，并比较选择了开采顺序，为改善金川矿区开采作业环境，保证矿山生产安全性提供了理论依据和数据上的支持.数值模拟方法目前已成为工程问题研究的重要手段[16-17].以往关于金川矿区的研究由于模型粗糙和计算方法的限制等问题，存在数值模拟结果与实际有所偏差的可能[18-19].本文借助于三维地质模型软件GOCAD中的离散光滑插值法建立了精细的矿体FLAC3D模型，针对长沙矿山研究院采矿中心提出的中段开采对照方案，中段的开采采取从左至右、从中间到两边、从两边到中间依次开采各个盘区的方案，通过数值计算综合比较不同开采方案下充填体、水平矿柱稳定性，选取最优开采方案.

1 金川矿区开采条件

位于甘肃金昌市的金川矿区是中国三大多金属共生矿之一，目前已有年产300万t的能力，是我国重要的镍钴生产基地.二矿区是金川镍矿床矿体重要矿区之一，其在规模、矿石品位、储量上在国内均占据显著优势.

1.1 金川矿区地质构造

金川矿区位于龙首山北东段北东侧，矿床主要以白云质大理岩、黑云母片麻岩及云母石英片岩为基底.如图1所示，二矿区位于金川矿区东侧，与一矿区以F16断层为界.二矿区岩体全长约3 000 m，总体走向为北西50°，倾向南西，倾角50°至80°.二矿区主矿体为1#和2#矿体，其中1#矿体位于该矿区西部，富矿体长约1 300余m，镍金属量占该矿区76.5%.本文以金川二矿区1#矿体为研究对象.

图1 金川矿区地质构造图[2]Fig.1 Geological structure of Jinchuan mining area[2]

1.2 金川矿区岩土力学参数

岩石力学参数的选取对数值模拟的准确性、现实可用性起着重要作用[20].金川矿区地质构造复杂，因变质作用影响，致使各类岩体力学性质变化剧烈，因此该地区岩体力学参数确定困难.本文参照各科研院所深入细致的调查结果，经过反复对比分析，选取表1中的符合FLAC3D软件计算特点数据作为计算参数.

表1 矿区岩体及充填体力学参数

1.3 金川矿区开采现状

金川二矿区采用大面积不留矿柱充填开采方式自上而下逐层开采，在垂直方向上矿体共划分为五个中段，每个中段又划分为若干分段，每个分段为20 m；在水平方向上矿体每100 m划分为一个盘区，盘区又划分为若干条带.该矿区开采工作分期进行，一期工程目标开采1 250 m以上的矿体，包括1 250 m中段和1 300 m中段，于1998年完成回采工作；二期工程已经完成1 150 m中段的回采工作，计划采用多中段同时开采的方式进行1 000 m中段和850 m中段回采，本文将对该问题展开研究.

2 数值模拟模型

2.1 计算几何模型

本文采用在岩土力学计算中广泛应用的FLAC3D程序进行数值模拟，其采用三维快速拉格朗日差分方法求解.为了模拟真实的回采环境、准确地再现矿体复杂边界，采用了三维地质模型软件GOCAD中的离散光滑插值法进行模拟.根据二维CAD数据，建立不同高程矿体边界控制曲线，在曲线上建立控制点，然后设置边界约束，由线生成面，进一步生成矿体模型，在矿体模型的建立过程中多次运用离散光滑插值方法.最终将其转化成相应计算网格(SGRID)模型，再通过FLAC3D中FISH语言编制成的接口程序提取模型网格节点坐标，导入计算软件中，建立最终计算模型(图2).该模型底部高程100 m，顶部高程1 500 m，共有196 000个单元，206 681个节点.模型中Z方向为竖直向上方向，X方向为矿体走向，Y方向垂直于矿体走向.

图2 计算模型Fig.2 Calculation model

2.2 本构模型

本文中采用以Mohr-Coulomb模型为屈服准则的理想弹塑性本构模型，用于模型围岩、矿体和充填体，FLAC3D以Mohr-Coulomb屈服准则(图3)为基础建立其本构模型，弹性计算准则为如下胡克定律增量表达式：

(1)

(2)

(3)

塑性屈服准则分为剪切屈服准则和拉伸屈服准则，公式表达分别为式(4)和(5).流动法则如式(6)所示.

Fs=σ1-σ3Nφ+2cNφ，

(4)

Ft=σt-σ3，

(5)

h=σ3-σt+ap(σ1-σp)，

(6)

式中,

在FLAC3D中规定σ1≤σ2≤σ3，故单元应力只处于σ3-σ1=0左侧，其左侧区域又被剪切屈服准则线(Fs=0)、拉伸屈服准则线(F-t=0)和流动法则线(h=0)划分为三个区域，其中区域1为弹性区，区域2为剪切塑性区，区域3为拉伸塑性区.

图3 FLAC3D中的 Mohr-Coulomb屈服准则Fig.3 FLAC3D Mohr-Coulomb failure criterion

2.3 初始地应力条件

金川二矿区矿体埋深大，地应力对采矿稳定性影响较大.因此确定地应力参数对数值模拟真实性影响较大，本文参考北京科技大学蔡美峰等[21]的研究成果，该成果已经被广泛应用，其准确性得到了验证.金川二矿区地应力场以水平构造应力为主，最大主应力方向为水平，垂直主应力基本以上覆岩体重力作用为主，水平主应力约为垂直主应力的1.69～2.27倍.对地应力测量结果线性回归得出以下主应力随深度变化的规律，其中，应力数值以压应力为负，拉应力为正.

σH=-0.098-0.0507H,

(7)

σh=0.015-0.0200H,

(8)

σV=0.208-0.0254H,

(9)

式中,σH为最大水平主应力，单位为MPa；σh为最小水平主应力，单位为MPa；σV为垂直主应力，单位为MPa；H为测点埋深，单位为m.

3 数值模拟结果

为了研究二矿区双中段合理回采顺序，设置三种中段开采方案，分别是从方案1左到右依次开采(图4(a))、方案2从两边到中间依次开采(图4(b))和方案3从中间到两边依次开采(图4(c)).

图4 开采方案Fig.4 Mining schemes

开采过程中，两中段相应位置同时进行，岩体稳定后，进行回填作业.首先依次对1 300 m中段、1 250 m中段和1 150 m中段进行开采和回填，然后对1 000 m中段和850 m按照设计方案进行开采回填模拟直至水平矿柱厚度为20 m.通过综合对比1 000 m水平矿柱、1 250m以上大体积充填体及围岩弹塑性状态、应力、沉降状况等来综合评价各方案的安全性.

3.1 应力计算结果分析

图5、图6中的应力(单位：MPa)的计算结果表明：三种采矿方案下，方案1水平矿柱1 000 m处应力释放更均匀，另外，三种方案下矿体边界处均有小区域应力集中现象.从图7、图8中可以看出，运用方案1采矿下1 250 m以上充填体应力状态显著优于方案2、方案3，应力释放更均匀，且无应力集中现象.三种采矿方案下，水平矿柱1 000 m处、充填体1 250 m处均无拉应力产生.

图5 水平矿柱1 000 m处最大压应力云图Fig.5 The cloud picture of maximum compressive stress of horizontal pillar at 1 000 m

图6 水平矿柱1 000 m处最小压应力云图Fig.6 The cloud picture of minimum compressive stress of horizontal pillar at 1 000 m

图7 充填体1 250 m处最大压应力云图Fig.7 The cloud picture of maximum compressive stress of backfill at 1 250 m

图8 充填体1 250 m处最小压应力云图Fig.8 The cloud picture of maximum compressive stress of backfill at 1 250 m

图9、图10为典型垂直截面处应力云图，由图可知，方案1对周围岩体扰动相对来说更小，整体应力释放更加均匀.三种方案均有小区域充填体产生拉应力，最大值分别为0.32 MPa、0.25 MPa、0.39 MPa，均小于充填体抗拉强度0.6 MPa，其中方案2产生拉应力的区域更大.

图9 截面处最大压应力云图Fig.9 The cloud picture of maximum compressive stress at section

图10 截面处最小压应力云图Fig.10 The cloud picture of minimum compressive stress at section

3.2 沉降计算结果分析

图11说明在开采过程中，1 000 m水平矿柱中部出现较大沉降，主要因为矿柱中部跨度更大.方案1中水平矿柱沉降(单位：m)最小，方案1水平矿柱沉降最大值为0.62 m，方案2、方案3中最大位移值分别为1.27 m、1.32 m.图12说明1 250 m充填体在边缘出现较大沉降，主要因为充填体位于接触面下方，围岩对其支撑作用小.对充填体支撑作用小.方案1中充填体沉降最小，方案1充填体沉降最大值为0.62 m，方案2、方案3最大值分别为0.69 m、0.71 m.

图11 水平矿柱1 000 m处沉降云图Fig.11 The cloud picture of settlement of horizontal pillar at 1 000 m

图12 充填体1 250 m处沉降云图Fig.12 The cloud picture of settlement of backfill at 1 250 m

图13中，方案2、方案3的1 250 m以上充填体沉降显著大于方案1，以1 000 m水平矿柱沉降来看，方案1的优越性更加明显.方案1中该截面最大沉降为0.68 m，方案2、方案3最大沉降分别为1.21 m、1.26 m.在850 m中段处，充填体因高水平构造应力作用，会发生向上隆起现象.从数值来看，方案1最大隆起为0.33 m，方案2、方案3最大隆起分别为0.44 m、0.46 m，且从分布范围来看，方案1发生隆起区域更小，整体隆起幅度小.

图13 截面处沉降云图Fig.13 The cloud picture of settlement at section

计算中，在1 000 m水平矿柱沉降较大区域中轴线设置三个监测点，位置如图2 (d)所示，监测了其在采矿过程中沉降变化，将三点监测结果取均值后绘制曲线图(图14).方案1中水平矿柱沉降在开采过程中一直小于方案2、方案3，且方案1中沉降变化均匀，没有发生陡然沉降.方案2、方案3在前期开挖过程中，矿柱突发约1 m沉降，开采过程中要引起重视.

图14 水平矿柱1 000 m沉降变化曲线Fig.14 The settlement curve of horizontal pillar at 1 000 m

3.3 弹塑性状态计算结果分析

图15是矿体1 000 m高程处截面，展示了水平矿柱底部岩体弹塑性状态.图16为充填体1 250 m高程处截面，展示了大体积充填体底部岩体弹塑性状态.从结果来看，1 000 m水平矿柱塑性区主要集中在大跨度区域中部，其中方案1塑性区整体比方案2、方案3面积小且没有贯通.1 250 m充填体状态三者均在充填体边缘产生塑性区，方案2、方案3塑性区比方案1更加贯通，面积更大，危险性更高.综合比较下，方案1即从左至右依次开采1 000 m中段与850 m中段更有利于开采的稳定.

图15 水平矿柱1 000 m处弹塑性状态分布图Fig.15 The elastic-plastic state distribution of horizontal pillar at 1 000 m

图16 充填体1 250 m处弹塑性状态分布图Fig.16 The elastic-plastic state distribution of backfill at 1 250 m

4 结论

本文以金川二矿区开采为例，采用有限差分法进行矿区开采中双中段开采回采顺序研究及水平矿柱和大体积充填体稳定性计算研究，综合对比分析三种开采方案计算结果，得出以下结论.

1) 相比其他方案，方案1围岩应力释放更均匀，水平矿柱、大体积充填体沉降更小，开采过程中变形稳定，均匀增长.方案2和方案3水平矿柱、充填体沉降远大于方案1，变形过程也出现急速变化，特别是方案2水平矿柱出现大面积塑性区.经过综合对比不同方案下采场稳定性表明方案1即从左至右依次开采各个盘区沉降较小，应力释放更均匀，更有利于保证采矿作业的安全，最终推荐金川公司采用该方案进行1 000 m中段和850 m中段的同时开采.

2) 1 000 m中段和850 m中段同时开采过程中，应重点监测1 000 m水平矿柱大跨度区域沉降，建议采取细分盘区，分步开采的方案，或在850 m中段内设立垂直保安矿柱.1 250 m以上大体积充填体易在边界发生失稳，可采用锚杆支护加强围岩对充填体的支护作用.