基于修正Mathews 稳定图法与FLAC3D 的阿舍勒铜矿深部回采方案优化研究

黄 聪，魏超城，丘永富

（1.长沙矿山研究院有限责任公司，湖南 长沙 410012；2.金属矿山安全技术国家重点实验室，湖南 长沙 410012；3.新疆哈巴河阿舍勒铜业股份有限公司，新疆 哈巴河 836700；4.新疆金宝矿业有限公司，新疆 富蕴 836199）

0 引言

经过十多年的开采，阿舍勒铜矿浅部资源逐渐枯竭，矿山逐步转入深部开采阶段。受到深部“三高一扰动”的复杂力学环境影响，岩体的力学特征及其工程也相应发生了明显变化。岩爆、冒顶、采空区失稳等地质灾害程度加剧，对矿山深部资源的开采产生了巨大威胁。因此，选择合理的采场结构参数和回采顺序对于控制围岩变形、预防采空区失稳有着积极作用。

近年来，众多专家学者对采场结构参数优化与回采顺序选择进行了深入研究。周科平[1]将遗传算法应用于采场结构参数的优化，并阐述了其基本原理。彭康等[2]通过工业试验验证了使用数值模拟来确定采场参数的可行性。赵永等[3]利用Mathews 稳定图法对赤峰红岭铅锌矿的采空区进行了稳定性及稳定概率评价。徐荃等[4]采用数值模拟与Mathews稳定图法相结合的方式对某金矿深部采场布置形式和结构参数进行了优化研究。孙健等[5]针对地下矿山采场结构参数应用Mathews 图解法对急倾斜采场暴露面冒落问题进行了研究。尹土兵等[6]认为采场暴露面积的确定直接影响采场结构参数的合理布置，采用阶段空场嗣后充填采矿法时，需要同时讨论采场顶板和充填体侧帮的暴露面积。特鲁曼等[7]收集了大量工程实例，与加拿大数据库合并分析后，提出了一种修正的Mathews 图解法用于预测采场帮壁稳定性。

以上研究对于矿山回采方案优化具有一定的指导意义，但由于阿舍勒矿区构造复杂，并不能直接套用。因此，使用Mathews 稳定图法与数值模拟相结合的方法，以阿舍勒铜矿深部采场参数、回采顺序等为优化对象，分析不同回采顺序与采场结构条件下，开采范围内的矿岩介质应力、位移、塑性区的变化过程，以期为矿山生产提供依据。

1 工程概况

阿舍勒铜矿矿体倾角为55o～75o，矿体产状总体走向近南北，倾向东，分布范围为水平方向3#勘探线～11#勘探线，空间范围+300～+50 m，其矿体赋存条件如图1 和图2 所示。目前阿舍勒铜矿的深部生产中段主要为+400～0 m 共9 个中段，中段高度为50 m，采用大直径深孔空场嗣后充填采矿法开采，采场垂直矿体走向布置，长度为矿体厚度，采场宽度20 m，开采顺序为由北向南隔一采一。

图1 +150 m 中段地质平面图Fig.1 Geological plane of +150 m middle section

图2 5#勘探线剖面图Fig.2 Profile of exploration line 5

阿舍勒铜矿深部开采目前存在贫化损失率居高不下、采场稳定性较差以及地压控制困难等诸多难题，选取具有代表性的+150 m 中段作为研究对象，旨在通过优化采场结构、选择合理开采顺序，达到控制地压、降低贫化损失率的目的，实现矿山安全高效生产。

2 基于修正Mathews 稳定图法采场尺寸确定

Mathews 图解法是MATHEWS 等基于大量工程实例在1980 年提出的一种矿山开采设计方法，最初应用于对采场稳定性进行预测。经过大量学者的研究与补充，目前Mathews 图解法已成为一种相对简单的岩体分级系统，被广泛应用于国外矿山采场设计。根据阿舍勒铜矿深部采场所使用的采矿方法，采用修正Mathews 稳定性图解方法分别分析采场顶板和采场上盘围岩区域的岩体稳定性[8]。

2.1 修正Mathews 稳定图法

Mathews 稳定性系数的计算公式见式（1）。Mathews稳定性系数与水力半径的相关关系如图3 所示。

图3 Mathews 稳定性系数与水力半径关系Fig.3 Relationship between Mathews stability coefficient and hydraulic radius

式中：N为Mathews 稳定性系数；Q'为修正的Q值；A为岩石应力系数；B为节理方位系数；C为重力调整系数。

应用Mathews 稳定性图表方法需要计算两个参数，即稳定性系数和形状因素，其中，稳定性系数反映了在一定的应力条件下岩体自立的能力；形状因素考虑了单独采场暴露表面的尺寸和形状，形状系数是暴露面的水力半径。

Mathews 稳定性图解方法采用了修正的NGI 隧道质量指标Q，与Q值不同的是，Q'值中的应力折减系数（SRF）和节理渗水折减系数（Jw）都为1.0，Q'值的计算公式见式（2）。

A值考虑高应力影响降低岩体稳定。A值为完整岩体的单轴抗压强度与平行开挖面的最大诱导应力的比值。A值与 σc/σi成线性关系，变化范围为0.1～1.0。岩石应力系数A的图解如图4 所示。

图4 岩石应力系数A 的图解Fig.4 Diagram of rock stress coefficient A

B值要考虑不连续面的方向影响，根据控制性节理与采场表面的相对方位确定，结构面与开挖面的夹角为90°时，B系数被赋值1，不连续结构面与开挖表面的夹角为20°时，B值为0.2，节理方位系数B的图解如图5 所示。

图5 节理方位系数B 的图解Fig.5 Diagram of joint azimuth coefficient B

受重力影响，采场顶板稳定性小于侧帮，重力调整系数C考虑了重力对采场暴露表面崩落、滑落等稳定性影响，重力调整系数C和采场表面倾角β的关系由式（3）确定。重力调整系数C的图解如图6所示。水力半径用于分析表面积除以暴露墙周长的比值，如图7 所示，计算公式见式（4）。

图6 重力调整系数C 的图解Fig.6 Diagram of gravity adjustment coefficient C

图7 水力半径确定方法的图解Fig.7 Diagram of hydraulic radius determination method

2.2 阿舍勒铜矿采场结构参数优化

2.2.1 稳定性系数和水力半径

修正的Q值的计算结果见表1。岩石应力系数A见表2。矿山深部以厚大、急倾斜矿体为主，矿体倾向东，倾角60°～80°。其上盘围岩主要为黄铁矿，矿体为铜硫矿，下盘围岩主要为凝灰岩。B值根据控制性节理与采场表面的相对方位确定，见表3。

表1 修正的Q 值计算结果Table 1 Calculation results of corrected Q-value

表2 岩石应力系数ATable 2 Rock stress coefficient A

表3 节理方位系数BTable 3 Joint azimuth coefficient B

受重力的影响，采场顶板的稳定性小于侧帮，重力调整系数C考虑了重力对采场暴露表面崩落、滑落等稳定性的影响。

1）采场顶板。采场顶板为水平状即β=0◦，重力调整系数C计算见式（5）。

2）采场上盘。采场上盘倾角一般为70°，即β=70◦，重力调整系数C计算见式（6）。

根据求得的Mathews 稳定性系数计算参数，按照式（1）计算得出采场顶板和采场上盘岩体的稳定性系数，运用图解法得出与稳定性系数对应的容许水力半径取值区间，结果见表4。

表4 稳定性系数计算结果Table 4 Calculation results of stability coefficient

2.2.2 采场顶板稳定性分析

矿山深部主要采用大直径深孔采矿法，采场垂直矿体走向布置，中段高度为50 m，矿体最大厚度为50 m，回采规划中采场宽度以12 m 为主，根据水力半径计算方法，计算得出不同采场宽度和长度条件下顶板的水力半径，计算结果如图8 所示。

图8 不同采场结构下顶板的水力半径Fig.8 Hydraulic radius of roof under different stope structures

在对深部岩体质量指标Q值计算的基础上，根据Mathews 图解法，在中段高度一定的条件下，对不同采场宽度和长度尺寸下采场顶板的容许水力半径进行计算分析，得出水力半径值≤4.8 m 对应的采场尺寸，其采场顶板处于无支护临界稳定区，整体稳定状态较好，满足采场安全生产要求。

2.2.3 采场上盘稳定性分析

由于深部开采中段高度固定为50 m，矿体最大厚度为50 m，上盘围岩的倾角平均为70°，根据水力半径计算公式，分别计算不同采场宽度尺寸下，上盘围岩的水力半径（以上盘斜长计算），计算结果见表5。

表5 不同采场上盘围岩的水力半径Table 5 Hydraulic radius of hanging wall surrounding rock in different stopes

结合表4 采场上盘围岩稳定性系数对应的容许水力半径值可知，当采场上盘的水力半径≤4.9 m，上盘围岩处于稳定区或无支护过渡区，则满足采场安全生产需要。即在采场高度一定的条件下，采场宽度≤12 m，可以保证矿体最厚处（矿房长度50 m）的采场上盘围岩整体稳定状态较好，满足采场安全生产要求。

3 深部开采回采顺序优化

根据分段空场嗣后充填法的特点，综合利用了CAD、Midas GTS、FLAC3D建立深部矿体三维数值模型，基于深部开采的实际情况、前期开展的工程地质调查和室内岩石力学实验结果，对不同的开采方案进行了数值模拟。在地质条件、采场结构尺寸一定的条件下，使用FLAC3D软件对各开采顺序方案进行模拟计算，比较采场附近的应力和位移，分析各回采方案的优劣，从而得到合理的回采顺序[9-14]。

3.1 数值模拟模型建立

深部+150 m 中段是研究对象，在网格划分时将150 m 矿体部分网格细化，并适当控制最大单元和最小单元尺寸间的比例以防止计算出现较大的偏差。

计算模型x方向563 m，y方向315 m，z方向250 m，单元大小矿体2m×5m×5m，围岩4m×5m×5m，模型单元总数31.2 万单元，计算模型如图9所示。

图9 数值模拟模型Fig.9 Numerical simulation model

3.2 模拟方案

阿舍勒铜矿深部采用从上至下的中段开采顺序，在中段开采顺序和采场结构参数确定的条件下进行采场回采顺序的模拟分析。根据矿床开采技术条件和大直径深孔空场嗣后充填法的特点，采场垂直矿体走向布置，矿房长度为矿体厚度，中段高度50 m，为分析不同的采场回采顺序对采场稳定性的影响，以+150 m 中段为模拟对象，共设计了四种不同的回采方案[15]（图10 和表6），即从中央向两端推进回采（方案1）、从两端向中央推进回采（方案2）、从北端向南端推进回采（方案3）、从南端向北端推进回采（方案4），每种回采方案矿房、矿柱的回采方式均为隔三采一。采场编号如图10 所示。

表6 模拟方案详情Table 6 Details of simulation scheme

图10 采场编号示意图Fig.10 Schematic diagram of stope number

3.3 计算结果及比较分析

对于数值模拟计算的结果，主要从应力分布状态、位移变化状态、塑性区分布状态等三个方面进行对比分析。

3.3.1 采场应力分布状态的比较分析

根据各方案数值模拟计算的结果，全部采场回采结束后的最大主应力等值线分布图如图11 所示。由图11 可知，在采场回采过程中，最大主应力集中区主要出现在两侧帮中部，最大主应力值-1.62～-0.95 MPa，侧帮岩体应力表现为受压状态，应力集中区域在采场上盘，采场易产生下盘岩体受压底鼓及上盘岩体受压冒落的现象。

图11 采场最大主应力等值线图Fig.11 Contour map of maximum principal stress in stope

根据数值模拟计算结果，对各方案中采场在回采过程中产生的最大主应力值进行统计，采场最大主应力值统计见表7。由表7 可知，方案4 中各采场的最大主应力值均小于其他几种方案，表明方案4中采场的最大主应力处于较低的应力状态。由各采场最大主应力变化趋势可以看出，方案1、方案2 中采场最大主应力值变化区间较大，部分采场集中应力值较大，方案3、方案4 中采场最大主应力值变化趋势相似，变化区间较小。由此表明方案1、方案2中局部采场应力集中突出现象明显，方案3、方案4中采场压应力变化较为平缓，应力在采场回采中分布较为均匀。

表7 采场最大主应力值Table 7 Maximum principal stress value of stope

3.3.2 采场位移变化状态的比较分析

根据数值模拟的计算结果，各方案采场回采结束后所产生的位移等值线分布图如图12 所示。由图12 可知，各采场的最大位移均位于采场的两侧帮和上盘岩体，侧帮位移值为19.9～44.9 mm，上盘岩体位移值为16.2～51.9 mm，顶板产生的位移则较小，位移值为4.6～37.8 mm（表8～表10）。

表8 采场侧帮最大位移值Table 8 Maximum displacement value of stope side slope

表9 采场上盘岩体最大位移值Table 9 Maximum displacement value of hanging wall rock mass in stope

表10 采场顶板最大位移值Table 10 Maximum displacement of stope roof

图12 采场位移等值线图Fig.12 Contour map of stope displacement

由表8～表10 可知，方案1 中回采采场在中间时侧帮、上盘和顶板位移值均较小，但随着采场向两端的推进，位移值明显增加，尤其端部采场各区域位移值较大，方案1 中各采场上盘、顶板的位移较其他方案均较大。方案2、方案3 各采场侧帮、上盘和顶板位移值较接近，在中部和端部的采场中位移变化较方案4 偏大。方案4 中采场位移较其他方案均较小，表明方案4 中采场位移变化较小，有利于采场的稳定，从位移变化的角度考虑，方案4 明显优于其他方案。

3.3.3 采场塑性区分布状态的比较分析

根据数值模拟的计算结果，各方案采场回采结束后所产生的塑性区分布图如图13～图16 所示。由图13～图16 可知，方案1 和方案2 中采场充填体的塑性区面积较接近，呈现大面积的贯通性拉张破坏，上盘和顶板岩体表现为剪切破坏，且方案2 中上盘岩体剪切破坏面积明显大于方案1，采场侧帮以拉张破坏为主，这与充填体强度较小有较大的关系，采场充填体较易受拉垮落。

图13 方案1 采场回采结束后塑性区分布图Fig.13 Distribution map of plastic zone after stope mining in scheme 1

图14 方案2 采场回采结束后塑性区分布图Fig.14 Distribution map of plastic zone after stope mining in scheme 2

图15 方案3 采场回采结束后塑性区分布图Fig.15 Distribution map of plastic zone after stope mining in scheme 3

图16 方案4 采场回采结束后塑性区分布图Fig.16 Distribution map of plastic zone after stope mining in scheme 4

方案3、方案4 中采场充填体塑性区面积相近，采场充填体的塑性区面积小于方案1、方案2，且塑性破坏区的也未完全贯通，这一结果有利于改善充填体内的受力状态，减少充填体的破坏程度，增强矿体二步骤回采时充填体的自稳能力。所以从塑性区分布的角度考虑，方案3、方案4 明显优于其他方案。

通过从应力分布状态、位移变化状态和塑性区分布状态三个方面对上述四种模拟方案进行对比分析可知，方案4 最优，其次为方案3。因此，建议深部采场中段回采顺序为由南端向北端推进回采。

4 结语

采用Mathews 稳定性图解法，根据深部采场的岩体特征，对采场顶板和上盘岩体的稳定性进行了计算分析。当深部采场的中段高度为50 m，采场长度为矿体最大厚度时，采场宽度为12 m 即可处于稳定状态。根据水力半径计算结果，只要采场长度小于50 m，采场顶板和上盘围岩处于无支护临界稳定区，整体稳定状态较好，满足采场安全生产要求。

通过三维数值模拟计算方法，对深部采场的回采顺序进行了模拟对比分析，共设计了四种回采方案，即从中央向两端推进回采、从两端向中央推进回采、从南端向北端推进回采、从北端向南端推进回采，对采场回采过程中的应力分布状态、位移变化特点和塑性区特征进行了对比分析，方案4（从南端向北端推进回采）最优，方案3 次之。因此，建议深部采场中段回采顺序为由矿体南端向北端推进依次回采。