落凼矿分段空场嗣后充填法采场结构参数研究

贾学元，乔登攀，何名声，王 俊，李广涛，杨天雨

(1.昆明理工大学 国土资源工程学院，昆明 650093；2.凉山矿业股份有限公司，四川 凉山 615141)

采矿方法是矿山安全生产的核心，采矿方法的选择直接决定了矿山生产运营过程中的效益及安全。现阶段随着国家政策的支持及矿山自我认识的提高，我国地下矿山采用充填法的比例逐渐提升[1]。分段空场嗣后充填采矿法采用两步骤进行开采，先采一步骤矿柱，隔一采一，对一步骤矿柱进行充填，充填体达到一定强度，进行二步骤矿房的开采，对二步骤矿房进行充填，二步骤矿房是在一步骤充填后留下的充填体矿柱支撑下进行开采，一步骤矿柱的充填质量是二步骤矿房安全开采的保障，该方法经济高效，同时又能保障井下安全作业。

采场结构参数是影响采场稳定性的重要因素之一[2]。在采矿方案设计时，若设计的矿柱尺寸过大、矿房跨度过小，将导致矿石损失大、回收率低；若设计的矿柱尺寸过小，矿房稳定性降低、生产成本提高，且易发生大规模地压活动[3]。国内学者对采场结构参数优化进行了大量研究，初步范围的确定有Mathew法[4]，单一数值模拟手段有FLAC3D[5]、ANSYS[6]、秩和比法[7]，采用两种数值模拟结合的有ANSYS-FLAC3D[8]，SUBPAC-FLAC3D[9]，考虑多种因素优化采场结构参数的有响应面法[10]、CRITIC法[11]、AHP法[12]和TOPSIS法[12]。

针对拉拉铜矿落凼矿区深部矿段地下开采遇到的高低压导致采场结构稳定性变差问题，研究基于FLAC3D数值模拟对采场结构参数进行优化。FLAC3D在数值计算过程中功能十分强大，但对于复杂赋存形态、不规则地质体的前处理建模上有不足之处[8]，研究用RHINO软件建模，导入FLAC3D数值模拟软件进行计算，克服了单一采用FLAC3D法存在的建模和网格划分困难等问题。

1 工程概况

落凼矿区具有多年露天开采历史，已形成了大、小两个东、西露天采坑。矿山近年来开展了落凼矿区深部的探矿工作，扩大了矿权范围，为延长矿山开采服务年限提供了资源条件。地下开采矿体形态较为复杂，矿体倾角多在15°～40°，矿体厚度多在2～32 m。

现阶段，落凼矿处于露天转地下开采的关键时期，矿体厚度在7～8 m以上水平至倾斜矿体和矿体厚度在30 m以上的倾斜及急倾斜矿体采用分段空场嗣后充填法进行开采，研究选取矿体厚度在30 m以上的倾斜及急倾斜矿体用FLAC3D数值模拟软件进行分析，优化采场结构参数，为进一步开采提供科学的指导。

2 建模方案

考虑到倾角厚度的变化及建模工作量，本次模拟选用平均倾角45°、厚度30 m作为矿体产状，为了克服FLAC3D建模、网格划分困难，选用RHINO建模后导入FLAC3D进行计算，如图1所示。为优化分段空场嗣后充填法采场结构参数，本次数值模拟设计两个方案进行对比，如表1所示。

图1 RHINO建模区域示意图Fig.1 Schematic diagram of RHINO modeling area

表1 数值模拟方案对比

方案一Ⅰ步骤矿柱沿走向、倾向均有布置，矿房宽度11 m，Ⅱ步骤矿房沿走向布置，矿房宽度25 m，先回采Ⅰ步骤矿柱，胶结充填后回采Ⅱ步骤矿房并充填。本次建模区域中包含二次暴露的Ⅰ步骤矿柱12个，其中沿走向6个，沿倾向6个。建模区域如图2所示，总单元数(zone)为901 050，节点数(grid point)为898 571。方案二与方案一的区别在于Ⅱ步骤矿房宽度为20 m。

图2 由RHINO导入FLAC3D建模区域示意图Fig.2 Schematic diagram of modeling area imported from RHINO to FLAC3D

2.1 计算模型的选择

本次数值模拟选用的本构模型为mohr-coulomb，各种力学性质岩体视为各向同性弹塑性连续介质，并在数值模拟计算过程中，应用摩尔-库伦塑性本构模型表示这种各向同性的塑性连续介质，其中摩尔-库伦屈服准则为：

(1)

ft=σt-σ3

(2)

式中，σ1，σ3，σt分别为最大主应力、最小主应力和抗拉强度；c，φ分别为黏聚力和内摩擦角，当fs>0时，材料将发生剪切破坏；ft<0时，材料将发生拉伸破坏。

2.2 力学参数的取值

本次数值模拟选取的力学参数如表2所示。

表2 计算岩石力学参数选取结果表

2.3 边界条件

落凼矿区矿体埋深较浅(300 m)，故竖直应力取自重应力，水平应力随深度按梯度变化。

σv=γH

(3)

(4)

式中，γ为矿岩容重，γ=ρg，g取9.81 m/s2，N/m3；H为埋深，m；μ为泊松比。

模型的四个侧面和底面均采用固定速度约束。

2.4 方案步骤模拟

矿房矿柱布置如图3所示。方案模拟开采过程共分为48步，见表3。

图3 模拟区域矿房矿柱布置示意图Fig.3 Layout of room and pillar in simulated area

步骤模拟计算初始地应力计算模拟三维模型原岩应力状态第一步(开挖1)开挖沿走向Ⅰ步骤胶结矿柱A1第二步(充填1)充填沿走向Ⅰ步骤胶结矿柱A1第三步(开挖2)开挖沿倾向Ⅰ步骤胶结矿柱B1第四步(充填2)充填沿倾向Ⅰ步骤胶结矿柱B1第五步(开挖3)开挖沿走向Ⅰ步骤胶结矿柱A7第六步(充填3)充填沿走向Ⅰ步骤胶结矿柱A7第七步(开挖4)开挖沿走向Ⅱ步骤矿房A4第八步(充填4)充填沿走向Ⅱ步骤矿房A4︙︙第四十六步(充填23)充填沿走向Ⅱ步骤矿房A17第四十七步(开挖24)开挖沿走向Ⅱ步骤矿房A18第四十八步(充填24)充填沿走向Ⅱ步骤矿房A18备注胶结矿柱胶结充填、矿房非胶结尾砂充填

2.5 数值模拟计算结果分析

分段空场嗣后充填法按表3所示的回采顺序对矿体进行模拟开挖、充填，Ⅰ步骤胶结矿柱和Ⅱ步骤矿房开挖结束后按表2中的强度设计、尾砂充填体参数进行赋值。分析胶结矿柱在区段回采过程中最危险状态即一侧临空另一侧受尾砂主动压力作用时的塑性区分布，并统计塑性区体积，从而判断胶结充填体的稳定性。

以下为两个方案Ⅰ步骤矿柱的塑性区分布情况，由于篇幅限制，无法列出所有开挖方案对应的所有矿柱塑性区分布图，每个方案只列举沿走向和倾向各两幅图展示。

所有开挖完成后，对沿走向Ⅰ步骤矿柱A8进行塑性区统计，开挖17沿走向Ⅰ步骤矿柱A8塑性区分布如图4所示，方案一塑性区总体积6 848.8 m3，矿柱总体积23 499.7 m3，塑性区占比29.14%；方案二塑性区总体积6 580.6 m3，矿柱总体积21 308.7 m3，塑性区占比30.88%，方案二比方案一塑性区占比大，稳定性较差。

图4 开挖17沿走向Ⅰ步骤矿柱A8塑性区分布Fig.4 Distribution map of plastic zone of Ⅰ step pillar A8 along strike in excavation 17

所有开挖完成后，对沿走向Ⅰ步骤矿柱A14进行塑性区统计，方案一塑性区总体积5 380.8 m3，矿柱总体积23 499.7 m3，塑性区占比22.90%；方案二塑性区总体积5 334.1 m3，矿柱总体积21 308.7 m3，塑性区占比25.03%，方案二比方案一塑性区占比大，稳定性较差。开挖23沿走向Ⅰ步骤矿柱A14塑性区分布如图5所示。

图5 开挖23沿走向Ⅰ步骤矿柱A14塑性区分布Fig.5 Distribution map of plastic zone of Ⅰ step pillar A14 along strike in excavation 23

所有开挖完成后，对沿倾向Ⅰ步骤矿柱B1进行塑性区统计，方案一塑性区总体积2 300.1 m3，矿柱总体积15 549.6 m3，塑性区占比14.79%；方案二塑性区总体积2 386.3 m3，矿柱总体积12 931.0 m3，塑性区占比18.45%，方案二比方案一塑性区占比大，稳定性较差。开挖8沿倾向Ⅰ步骤矿柱B1塑性区分布如图6所示。

图6 开挖8沿倾向Ⅰ步骤矿柱B1塑性区分布Fig.6 Distribution map of plastic zone of Ⅰ step pillar B1 along dip in excavation 8

所有开挖完成后，对沿倾向Ⅰ步骤矿柱B6进行塑性区统计，方案一塑性区总体积878.1 m3，矿柱总体积15 549.6 m3，塑性区占比5.65%；方案二塑性区总体积1 620.0 m3，矿柱总体积12 931.0 m3，塑性区占比12.53%，方案二比方案一塑性区占比大，稳定性较差。开挖24沿倾向Ⅰ步骤矿柱B6塑性区分布如图7所示。

图7 开挖24沿倾向Ⅰ步骤矿柱B6塑性区分布Fig.7 Distribution map of plastic zone of Ⅰ step pillar B6 along dip in excavation 24

部分沿走向矿柱塑性区统计结果如表4所示，为了方案对比便利，绘制成柱状图，如图8所示。在进行开挖时，两个方案A15矿柱塑性区占比近乎相同，其它矿柱方案二比方案一塑性区占比要大，说明方案一比方案二更加稳定。

表4 沿走向矿柱塑性区统计结果

图8 沿走向矿柱塑性区占比Fig.8 Proportion of plastic zone along strike pillar

部分沿倾向矿柱塑性区统计结果如表5所示，为了方案对比便利，绘制成柱状图，如图9所示。在进行开挖时，方案二矿柱整体要比方案一矿柱塑性区占比大，说明方案一比方案二更加稳定。

图9 沿倾向矿柱塑性区占比Fig.9 Proportion of plastic zone along dip pillar

表5 沿倾向矿柱塑性区统计结果

2.6 方案综合分析

矿体厚度在30 m以上的倾斜及急倾斜矿体选用分段空场嗣后充填法数值模拟，该类型矿体采用方案一和方案二数值模拟计算。对比两个方案的塑性区占比，对两个方案的稳定性情况进行综合分析，两个方案的数值模拟如前所示。

由数值模拟分析结果来看，沿走向矿柱比沿倾向矿柱塑性区占比更大，破坏更严重，塑性区集中在采场顶部和临空面，矿柱整体性较好。随开采区域的增大，顶板暴露面积增大，作用于胶结矿柱的松脱地压增大，超过胶结矿柱顶部设计强度时，顶部发生破坏。在模拟区域回采过程中，胶结矿柱受四周矿房回采的影响多次暴露，围岩应力重新分布过程中，胶结矿柱应力状态随之变化，由矿岩包裹的三维力学状态转变为单侧暴露的三维力学状态，胶结矿柱变形向空区内发展，临空侧单元体应力状态沿加载路径向强度曲线靠近，部分单元体应力状态达到极限应力状态发生破坏。

由方案一和方案二对比可知，矿柱在50 m×11 m参数不变条件下，在整体上方案二(矿房50 m×20 m)沿走向和倾向上矿柱塑性区占比明显比方案一(矿房50 m×25 m)大，稳定性更差，故选择方案一。

3 结论

1)随着采深的加大，拉拉铜矿落凼矿深部开采按原矿柱、矿房的采场结构参数设计势必导致地压显现频繁，需要优化采场结构参数指导进一步的开采。

2)落凼矿深部采用分段空场嗣后充填法开采，塑性区并未贯穿胶结矿柱，胶结矿柱整体性较好，采用方案一(Ⅰ步骤矿柱50 m×11 m，Ⅱ步骤矿房50 m×25 m)和方案二(Ⅰ步骤矿柱50 m×11 m，Ⅱ步骤矿房50 m×20 m)的采场结构参数均不会出现失稳破坏。胶结矿柱塑性区主要分布于胶结矿柱顶部和临空面，而顶部出现塑性区的主要原因是充填体在分层充填过程中顶部选用的充填配比强度较低导致。

3)通过FLAC3D数值模拟软件对采区进行参数优化，综合安全稳定性方面考虑，最终选择了Ⅰ步骤矿柱50 m×11 m，Ⅱ步骤矿房50 m×25 m为最佳采场结构参数，为拉拉铜矿落凼深部矿区的进一步开采提供了科学的指导，也对采用分段空场嗣后充填法的相似矿山有一定的借鉴意义。