残矿矿柱回采方案数值模拟与优化研究

陈祉颖，刘志祥，李夕兵，黄麟淇，王春元，李炎峰，张建涛

(中南大学 资源与安全工程学院，湖南 长沙，410083)

随着地下开采技术的不断发展，矿床开采会导致高地压聚集在围岩及矿柱上而诱发大规模冒落及大变形坍塌等矿山灾害，对矿柱与充填体稳定性亟需深入研究[1-2]。在赤峰柴矿上部矿床留矿法开采过程中，存在大量间柱和顶底残柱，随时间推移，地压转移已造成局部顶底残柱破碎，矿床开采往深部延伸，上部残矿资源变形与沉降问题将越来越严重，若不及时回收，顶底残矿资源开采难度将升高，开采安全性变差。为最大化地开发利用矿产资源，提高矿石回收率，在保证采场安全、稳定的前提下有效回收暂留残矿柱，需要制定科学合理的采场结构参数和开采方案。

针对矿柱回采的稳定性问题，国内外学者主要从经验公式、理论分析、模糊可靠度分析和数值模拟等方面进行研究[3-6]，陈顺满等[7]运用普氏理论和BIENIAWSKI 矿柱强度公式推导出方形矿柱的安全系数计算公式，用以优化矿柱宽度，提高矿石回收率；ZHOU等[8]基于压力拱理论方法和结构动力学，建立了评估作用于与动态拆除支柱相邻的支柱上的载荷的力学模型，结果表明通过增加矿柱回收时间可减少动态矿柱扰动；CHEN等[9]基于条形矿柱应力应变本构关系和尖点突变理论，分析了多条形矿柱-顶板结构系统的非线性动力失稳机理，推导出不同类型结构系统维持稳定的矿柱宽度表达式。史秀志等[10]通过突变理论推导与数值模拟相结合的方法，得出条形矿柱最优宽度为5 m，具有较好的稳定性。于清军等[11]应用理论计算和ABAQUS 数值模拟软件对6 种桃形矿柱稳定性进行分析，得出桃形矿柱坡角为60°、总高为15 m 的方案为最佳方案。莫卿等[12]采用模糊评价理论，得出不同开采技术条件下的采场最优结构参数。LIU 等[13]通过FLAC3D数值模拟，选出合理的矿柱框架构建方案，并利用模糊综合评判法进行验证；LI 等[14]利用DIC 数字成像系统，进行了双中段的深部房柱型矿柱回采实验，分析得到回采过程中模型的沉降过程分为稳定位移及不稳定位移两个阶段，第一个回采矿柱上方的沉降最大；寇永渊等[15]利用FLAC3D数值模拟软件，揭示二矿区井下充填体及地表在整个回采过程中的动态地压分布规律，优化选择了矿柱回采方案。李俊平等[16]采用FLAC3D数值模拟了矿柱回采及深部卸压开采的应力变化过程，提出间隔矿柱抽采法回采间柱及“V”型松动爆破并上盘底板松动爆破隔断开采两种方式。

综上所述，为了保障矿柱安全高效回采同时使残矿资源得到最大化利用，本文作者根据柴胡栏子金矿I-4号井残矿资源的赋存状况，以确保生产安全前提下提高出矿量，提出4种不同残矿柱资源回收方案，采用ABAQUS 数值模拟方法，揭示间柱及顶底柱回收的应力位移时空演化特征，分析不同矿柱留设方案的回采率、应力演化规律、位移特征和塑性应变变化，利用模糊数学法优化整体残矿回收模式，并选出残矿资源最佳回收方案。

1 工程概况

1.1 矿区开采背景

柴胡栏子金矿区位于内蒙古自治区赤峰市，区内分为Ⅰ号和Ⅱ号两个矿化带，其中Ⅰ号矿化带长约1 100 m，宽50～200 m，Ⅰ号矿化带内共控制或探明矿体13 条，I-4 号矿体呈脉状、不规则透镜状、扁豆状及囊状产出，矿体控制长为460 m，控制斜深为270 m，矿体走向90°，倾向南，倾角为59°～70°。矿体厚度变化范围为0.37～14.44 m，平均厚度为2.89 m，厚度变化系数98.85%。矿体Au品位为1.17～23.10 g/t，三中段(标高810 m)：Au品位平均为4.94 g/t，四中段(标高770 m)：Au 品位平均为4.97 g/t，五中段(标高730 m)：Au品位平均为5.12 g/t，六中段(标高700 m)：Au品位平均为4.02 g/t。

该矿山现阶段采用浅孔留矿法开采，空区用废石进行回填。在上部矿床开采过程中，存在大量间柱和顶底残柱，顶底柱高度为3～5 m。根据残矿资源调查结果，早期开采剩余残矿资源主要分布在I-4 号矿脉，残矿资源量为77 998.27 t，金金属质量为261.57 kg，该部分残矿资源形成时间长，局部地压显现明显，急待回收。图1 所示为I-4 号矿体可回收残矿资源纵剖面示意图。

图1 I-4号矿体可回收残矿资源纵剖面示意图Fig.1 Schematic longitudinal section of recoverable residual ore resources from orebody I-4

1.2 岩石稳定性能与岩体力学参数

研究岩石的力学性能与质量分级对掌握残矿中岩石稳定性至关重要，并为采场结构参数及采矿工艺优化提供了理论依据[17]。按照国际岩石力学学会(ISRM)的力学试验标准对取自柴胡栏子金矿的试样进行加工，并开展室内岩石力学实验。试样分为上盘围岩、下盘围岩和矿体三类，主要进行单轴抗压试验、抗拉试验和抗剪强度试验，如图2所示。其中单轴抗压强度较高，柴胡栏子金矿的矿岩均属于坚硬的脆性岩石且上盘围岩的强度高于下盘围岩和矿体的强度，下盘围岩的强度最低。但岩石稳定性不仅与强度有关，而且节理面对其影响巨大。

图2 现场取样试件Fig.2 Sampling specimens on site

结合多种现有的岩体质量分级方法，以RMR分级法为基础，参照节理岩体的岩石力学分类指标表，对取自于赤峰柴矿的钻孔岩芯的完整岩石强度、岩石质量指标(RQD)、节理间距、节理状态、地下水条件和节理方向对RMR 的修正值等6个分级参数依次进行打分并累加，得到赤峰柴矿地下残矿回收地点的岩体质量[18]。矿区岩体分级结果显示，柴胡栏子金矿I 号井一—五中段10～20号勘探线采场岩体质量等级基本为III 级，少部分区段有F1 断层破碎带穿过，岩体质量等级为IV级。

根据岩体质量和室内岩石力学试验，以及工程勘探修正后得到岩体力学参数如表1所示。

表1 岩体力学参数Table 1 Mechanics parameters of main rocks mass

2 残矿柱回采数值模拟

2.1 计算模型建立

由于矿块开挖后采空区上覆岩层往往存在较大的变形与位移，岩体为非线性弹塑性介质，因此，选用非线性求解较强的有限元软件ABAQUS进行模拟计算[19]。鉴于井下开采技术条件复杂多变，为便于建模和分析计算，作出如下假设：

1) 矿体和充填体为局部均质、各向同性材料；

2) 将矿体和充填体视为非线性弹塑性介质，并且服从摩尔-库仑屈服准则；

3) 分析中涉及所有物理量均与时间无关；

4) 不考虑应变硬化(或软化)。

模型选取柴矿I-4号矿脉三中段至六中段，高为160 m，范围为沿矿体走向方向(X方向)长420 m，竖直方向(Y方向)上地表左侧取342 m，右侧取274 m，围岩宽为50 m。共包括3个部分，按照材料特征分为围岩、矿体和废石充填体，其中开采矿房部分按碎石充填处理，充填材料参数如表1所示。

2.2 地应力参数确定

岩体地应力是进行岩体工程问题数值计算的最基本也是最重要的初始条件。采用岩芯新型声发射法测试赤峰柴矿地应力。根据过往工程经验及理论分析，垂直主应力随深度h呈线性增大，即

式中：σv为垂直主应力，MPa；h为埋深，m。

水平主应力随深度呈现线性增长关系，为方便开采设计，利用线性回归法将各埋深点所测地应力进行拟合，得出赤峰柴矿岩体中地应力随深度变化的计算公式：

式中：σmax为最大水平主应力，MPa；σmin为最小水平主应力，MPa。

因此，该模型采用Mohr-Coulomb 准则[20]，在模型两侧和底部施加约束，矿体顶部施加均匀荷载，整体模型施加重力荷载。本次模型高度约为420 m，距离地表(+915 m)最小埋深为65 m，最大埋深为292 m，将地应力测试结果转化为X和Y方向应力，代入模型中进行计算。

2.3 矿柱回收方案设计

方案设计以最大矿柱回收率为准则，回收率大于等于70%为佳。考虑到地应力影响，需首先进行地应力平衡，采矿过程可认为是准静力过程，不考虑冲击，采用静力分析步进行模拟仿真，第一步模拟矿房开采充填，之后从上到下依次模拟各中段残柱的回采充填。开采时矿柱受上覆岩层与充填体压力作用，垂直荷载是破坏的主要原因，将水平矿柱当作固定梁，跨度较大，承载充填体与上盘围岩重力作用，垂直矿柱视为桩，承载围岩与充填体侧向挤压作用，因此，为选取最佳的顶底柱回采方案，按照顶底矿柱不同留设位置共建立4 种矿柱回采方案，如图3 所示。方案1 为阶梯支撑矿柱回采，基本回收全部间柱，留底柱成阶梯状，其矿柱整体回收率为74.00%，见图3(a)；方案2 为T 型支撑矿柱回采，如图3(b)所示，T 型矿柱支撑上下盘围岩，其余矿柱全部回采，其矿柱回收率为79.06%；方案3所留部分矿柱为Z型矿柱和2 个竖条状矿柱，其余矿柱全部回采，如图3(c)所示，矿柱回收率为76.96%；方案4所留矿柱为Z 型矿柱和3 个竖条状矿柱，见图3(d)，其余矿柱全部回采，矿柱回收率为76.08%。

图3 矿柱回采方案示意图Fig.3 Schematic diagrams of pillar recovery program

3 数值模拟结果分析

3.1 不同矿柱回采方案应力分析

图4所示为各回采方案开采步骤中所留矿柱的应力变化。开采初期，矿柱所受应力较小；随着矿床三中段至六中段的开采，矿柱拉应力均不高于抗拉强度，矿柱未受到拉伸破坏，而部分矿柱主要呈压剪破坏，开采初期充填体保持着弹性工作状态，矿柱作为主要承载体系，拉应力与剪切应力在四中段开采时会剧烈释放，随着充填体被不断压紧密实，所留矿柱上的应力将进一步加剧。

图4 各采矿方案的留设矿柱各分析步的最大应力Fig.4 The maximum stress of each analysis step of left pillars for each mining plan

3.1.1 垂直应力演化规律

各方案的矿柱垂直应力云图如图5所示。由图5 可知：4 种方案均有小区域拉应力产生，但所留矿柱所受拉应力均不超过抗拉强度5.81 MPa，未受到拉伸破坏。方案1所留矿柱呈阶梯状，在开采过程中，矿柱上应力逐渐向充填体转移，充填区域应力变为拉应力，开采完成后应力主要集中于顶底板间的间柱上，充填体受拉有上下移动的趋势，致使顶底板间的间柱部分区域处于受拉状态，其拐角处外侧受拉严重，最大拉应力发生在五中段的间柱上，并且所留矿柱整体柔性较好，在回采过程中压应力均匀增加，矿柱未受到压缩破坏；相反，方案2 形成了T 型的支撑结构，结构刚度大，在采矿回填的过程中，上覆岩层的压力全都作用于四中段底板处，底板弯曲下沉产生拉应力，间柱主要承受压缩应力，且支撑矿柱间距较大，易产生应力集中，最大压应力高达72.79 MPa；方案3 所留矿柱有单独竖条矿柱与Z 型矿柱共同组成支撑结构，对充填体变形阻碍小，拉应力主要集中在顶板与间柱连接的拐角处，最大为5.77 MPa，随着回采充填过程释放压应力，最大压应力主要出现在采场的角隅部分，少量分布于五中段的间柱上，最大压应力为56.93 MPa。方案4 所留矿柱应力主要集中于长竖条柱上，最大拉应力出现在五中段顶板与间柱拐角处外侧帮，最大压应力在内侧帮处，上覆回采充填体压力集中作用于五中段底柱，并传递给相邻的2个短间柱。综上可知，方案4的拉应力与方案1的压应力最低。

图5 回采后各方案矿柱模型垂直应力云图Fig.5 Vertical stress cloud maps of pillar model of each scheme after mining

3.1.2 剪切应力演化规律

由于地应力不均衡和充填体、矿柱或围岩有相对位移趋势，造成残留矿柱中弹性应变能量集中，这4 种方案间柱都受到明显剪应力集中作用。图6所示为矿柱回采方案的剪切应力云图。从图6 可知：主要在矿房顶柱与间柱的交接处及矿房和围岩交接处发生剪切应力破坏。方案1中，随着三中至六中段的开采，充填体和矿柱或围岩交接处剪切破坏加剧，最大剪切应力逐渐变大，Z型支撑结构对充填体变形阻碍小，六中段间柱作为桩，造成局部区域切应力集中，最大剪切应力为25.61 MPa。方案2的顶板贯穿于整个矿房，上覆岩层的压力由顶板传递至竖条间柱承载，在间柱与顶板连接处受挤压作用集中显著，集聚能量大，产生剪切破坏，稳定性较差，最大剪切应力为30.4 MPa。方案3和方案4 的矿柱留设结构合理，应力释放更加均匀，残矿开采充填后，避免应力集中造成较大破坏，最大剪切应力分别为23.18 MPa 和23.52 MPa。因此，方案中，矿柱的剪应力由大到小顺序依次为方案2、方案1、方案4和方案3。

图6 回采后各方案矿柱模型剪切应力云图Fig.6 Shear stress cloud maps of pillar model of each scheme after mining

3.2 不同矿柱回采方案位移场分析

由于矿体结构部分区域的应力已释放，从应力上无法直观地展示其受扰动前后应力状态对矿岩的影响，而位移是研究开挖过程中岩体变化的重要参数，矿体结构的应力释放区域通常是位移较大区域，因此，分析位移变化以研究开采过程对矿体的影响。矿柱回采后，图7所示为矿床各开采步骤下的水平与垂直位移变化。由图7可知：开采前位移小于0.01 m，方案1的矿柱结构为阶梯型矿柱结构，水平与垂直位移变化最大，方案2～4的水平位移在5～7 cm 间变化，垂直位移在20～35 cm间变化，垂直位移比水平位移变化更为显著。由于垂直位移过大导致顶板断裂、冒落甚至地表沉陷，顶板的垂直位移对采场稳定性的影响要比水平位移重要。

图7 各方案矿体模型的各分析步的最大位移变化Fig.7 The maximum displacement change of each analysis step of ore body model of each scheme

矿床整体模型的垂直位移变化云图如图8 所示。由图8可知：当矿柱被全部开采后，构造应力和重力作用下会出现上盘沉降、下盘凸起现象，方案1、方案2 与方案4 的起始矿柱起到较好的支撑作用，残柱开采后，受压沉降，沉降范围从开采区域中部逐渐向周围减小，最大位移出现在中部区域，方案1的阶梯型结构稳定性较差，位移变形显著区域最大，水平与垂直位移最大；方案3的结构整体性较好，位移较均匀。方案1至方案4的最大垂直位移分别59.08、34.83、24.65 和35.69 cm，相比之下，方案3的位移最小。

在某种程度上，人是决定工程能否取得成功的关键因素。这是因为几乎所有项目都是需要通过人来进行组织与管理的。拥有一支高素质、高水平的队伍，是打造高品质项目的前提。因此，在现场管理中必须做好人员管理，具体措施包括。

图8 矿柱模型垂直位移云图Fig.8 Vertical displacement cloud maps of mine pillar model

3.3 不同矿柱回采方案塑性应变分析

由以上受力分析可知，随着开采的进行，塑性应变也发生相应变化，当残留矿柱开采后，二中段以下间柱和顶底板产生较大面积的塑性应变，由于所留矿柱或回采矿柱后附近充填体经历反复的采动影响，更易造成塑形破坏，充填体产生部分拉伸破坏区；本文重点关注的所留保护矿柱，易在间柱和顶底柱转折处发生剪切应力集中现象，造成塑形破坏，而充填体刚度小、变形大，减小了矿柱的应力集中程度，起到矿柱保护作用。矿柱塑性应变变化云图如图9所示。从图9可知：方案1的塑性应变最大，在顶底柱与间柱转弯处产生塑性破坏，分布于各中段矿柱上，矿柱极有可能发生断裂造成失稳，危险性最高，方案1与方案4的最大塑性应变均发生在顶底柱与围岩连接处，方案2与方案3 的最大塑性应变均发生在顶柱与间柱交界处，方案3的塑性应变整体比另3种方案的小，稳定性最强，更有利于开采的稳定性。

图9 各方案矿柱回采后模型塑性应变云图Fig.9 Cloud map of model plastic zone after mining of each project pillar

4 矿柱回收方案综合分析优选

根据前面4种方案的计算结果发现，不同矿柱回收方案各有优缺点，凭主观难以判别，为此，基于层次分析法和模糊数学理论，对矿柱回收方案进行评价。首先建立各方案的综合评价指标体系，运用层次分析法客观确定各因素的权重，并通过模糊数学理论建立模糊综合评价，最终选取最佳残留矿柱回采方案。

4.1 综合评判定指标

针对4个方案的模拟计算结果，选取最大压应力、最大拉应力、最大位移、回采率、塑性应变、剪应力6个指标，各矿柱留设方案的主要指标比较如表2所示。

表2 各采矿方案的主要指标比较Table 2 Comparison of main indexes of mining schemes

4.2 量化指标的隶属度

选择的定量指标为各方案的最大压应力、最大拉应力、最大位移、回采率、塑性应变和最大剪应力指标。定量指标的隶属度由隶属函数法确定，非定量指标采用相对二元比较法确定。

假如每个方案有n个方案和m个指标，目标特征值矩阵可以写成：

式中：i=1，2，…，m；j=1，2，…，n；n=4；m=6。

目标相对优属计算如下。

对越大越优指标，用以下公式进行规范化处理：

对越小越优指标，用以下公式进行规格化处理：

得到目标相对优属度矩阵R6为

4.3 层次分析法

应用AHP 分析决策问题时，首先要把问题条理化、层次化，构造出一个有层次的结构模型，这些层次可以分为三类：最高层(目的层)、中间层(准则层)和最底层(方案层)[21]。递阶层次结构中的层次数与问题的复杂程度及需要分析的详尽程度有关，一般的层次数不受限制，每一层次中各元素所支配的元素一般不超过9个。

运用层次分析法可得到6 个指标的判断矩阵R：

求得最大特征值λmax=6.32，可计算一致性，判断矩阵一致性检验的公式为

其中：RC为一致性比率；IC为不一致程度的量化指标；IR平均随机一致性指数；n为判断矩阵的阶数。

通常认为，RC＜0.1 时通过一致性检验，此时，认为判断矩阵的一致性是可接受的，否则，应对判断矩阵作适当修正。根据表3得，RC=0.052，低于0.100，满足一致性检验。

表3 平均随机一致性指标Table 3 Average random consistency indicators

因此，利用方根法可得影响矿柱留设方案选择的指标权重向量为：ω=[0.036，0.064，0.146，0.495，0.175，0.084]。

4.4 矿柱留设方案的优化选择

赤峰柴矿残矿资源开采中，采用加权平均模型对上述4种方案从最大压应力、最大拉应力、最大位移、回采率、塑性应变、剪应力6个指标，进行综合评价。4 种方案的权重为：W=[0.730 4，0.829 5，0.949 5，0.826 6]。根据最大隶属度原则，W越大，方案越好，即各方案优先顺序为方案3、方案2、方案4和方案1。因此，方案3为回收残矿资源最佳方案。

5 结论

1) 针对赤峰柴矿残矿资源开采技术特征，提出4个残矿资源整体开采技术方案，即阶梯支撑矿柱回采、T 型支撑矿柱回采、Z 型加2 个竖条型支撑矿柱回采和Z 型加3 个竖条型支撑矿柱回采方案。

2) 矿柱回采后，间柱应力显著增加，且在间柱和顶底柱转折处易发生应力集中，造成塑形破坏，承载能力差。

3) 顶板垂直位移比水平位移对采场稳定性的影响显著，垂直位移变化范围为25～60 cm。开采前位移小于0.01 m，开采后顶底板出现沉降，沉降范围从中部逐渐向周围减小，最大位移出现在开采中部区域，“Z 型与2 个竖条型支撑”矿柱回采方案的位移最小。

4) 选取最大压应力、最大拉应力、最大垂直位移、回采率、塑性应变和剪应力作为指标，采用模糊数学优选法进行比较，得到方案优劣顺序为：方案3(Z型加2个竖条型支撑矿柱回采)，方案2(T 型支撑矿柱回采)，方案4(Z 型加3 个竖条型支撑矿柱回采)，方案1(阶梯支撑矿柱回采)。因此，赤峰柴矿采用“Z 型加2 个竖条型支撑矿柱回采”方案回收残矿资源最佳。