基于FLAC3D模型动态更新的自然崩落法开采过程模拟技术*

方传峰，余传玉，史存丁，冯兴隆

(1.中南大学 资源与安全工程学院，湖南 长沙 410083；2.云南迪庆有色金属有限责任公司，云南 香格里拉，674400)

0 引言

自然崩落法作为一种低成本，高效率，以及安全性好的采矿方法，适用于大型低品位矿山的开采，自1895年在Pewabic铁矿首次应用起，技术和装备不断发展，在国内外得到了大量应用[1-2]。该方法是通过简单的拉底和切割工程，在拉底空间上方利用矿体自身的软弱结构面，在重力及次生构造应力作用下，诱导矿体失稳崩落，从而省却大量凿岩爆破工程及费用。

自然崩落法开采过程主要受矿床开采技术条件，以及拉底和出矿过程控制，通常采用经验推导法、理论计算法、相似物理试验法以及数值模拟法等分析和确定上述因素对崩落矿岩时空演化特征的影响[3-5]。随着计算机技术的发展，数值模拟法已成为国内外研究学者用于预测自然崩落力学响应特征的有效途径。王涛等[6]采用PFC2D研究矿岩在崩落过程中的裂纹扩展及接触力分布情况；王连庆等[7]以某镍铜矿的地质条件及矿岩物理力学性质为背景,采用PFC2D数值模拟的方法分析了自然崩落法的崩落规律；孙闯等[8]采用FLAC3D软件，在获取泥岩的粘结力软化和内摩擦角硬化参数的基础上预测大跨度的泥质顶板冒落失稳； Vyazmensky等[9]采用实例统计、经验图分析与有限—离散元软件数值模拟研究地表沉降的影响因素；Woo等[10]采用FLAC3D与干涉合成孔径雷达(InSAR)分析、评价与预测矿块崩落法引起的地表变形，并修正数值模拟参数。

PFC离散元方法不能进行大规模开采模型的计算和分析[11]，FLAC3D虽可以分析大规模问题，但在以往的应用研究中，未引入网格动态生成机制，不能真实反映矿岩崩落和出矿过程的动态性，也不能及时生成底部结构承重矿堆，从而影响了对自然崩落过程力学响应研究的准确性。本文基于FLAC3D软件设计工作流程，并利用软件内嵌FISH逻辑语言进行编程与模拟，及时消除已崩落单元且同时在底部结构生成相应质量矿堆，以实现开挖过程的动态模拟。以西南某铜矿为例，通过建立二维矿山模型，对其进行拉底模拟，以分析显式崩落流程设计的有效性。

1 崩落过程动态模拟技术流程

以实际矿山崩落为动态破坏过程。传统连续介质算法模拟自然崩落，虽可识别各步拉底对应崩落区，但未做破坏矿岩单元的动态更新，既不能体现崩落过程的动态性，更重要的是导致模拟结果偏于保守，与实际工程相比，无论在空间范围，还是在崩落量上，均出现较大偏差，制约了生产计划和监测方案的准确制定。在对自然崩落法开采过程特征进行分析的基础上，提出和设计了基于FLAC3D的过程动态模拟算法，编制了相应的代码。

1.1 模型动态更新实现思想

首先，按照传统FLAC3D力学分析步骤和方法建立模拟模型。

其次，在设定每一步骤拉底范围后，采用摩尔-库伦模型模拟计算其力学性态，并按照考虑位移、最大主应力或塑性状态的单指标判断准则，或同时考虑多个因素的综合指标判断准则，确定潜在破坏和需要动态处理的单元。

最后，利用FLAC3D编写FISH函数修改单元的本构模型，实现破坏单元的消除和崩落面和底部结构之间单元的显现，即将发生破坏的单元类型修改为空模型(null)，原属于空模型的单元修改为其他本构类型，并修改其力学参数以表征松散矿堆。

1.2 模型动态更新基本原则

基于FLAC3D单元的消除与再现[12]，程序主要遵从3个基本原则，实现模型动态更新和崩落过程模拟。

1)质量守恒原则

矿岩破坏后发生冒落，若忽略矿岩的吸水性且考虑出矿作用，则质量恒定，满足式(1)：

m崩=m出+m堆

(1)

式中：m崩为崩落矿石质量；m出为出矿量；m堆为出矿后矿堆质量。

2)体积平衡原则

由于破碎后矿岩自身的碎胀性，体积会膨胀。为满足底部结构存在足够体积的空区顺利生成矿堆，则体积应满足式(2)：

V空+V崩≥V堆

(2)

式中：V空为已存在空区体积；V崩为待崩落矿岩体积；V堆为待沉积矿堆体积。

3)拱形发育原则

底部结构沉积的矿堆呈散体状态，矿堆按照一定的自然安息角呈拱形发育。程序根据现有空区平面形态和位置，按设定的安息角在一定拱形范围内重现单元，其形态见图1。

1.崩落区；2.已开挖拉底层；3.已开挖聚矿槽；4.已开挖进路；5.矿堆。图1 矿堆示意Fig.1 Mine heap schematic

1.3 动态过程模拟流程

基于上述3项原则，引入人工控制变量，编制Fish程序，具体动态过程模拟流程见图2。

图2 显式崩落实现概图Fig.2 The overview of realization of explicit caving

图2为显式崩落实现概图，首先确定每步拉底范围，其次对第N步模拟拉底，运行平衡后，根据破坏判断准则选定破坏单元，基于三项原则实现破坏单元的消除，并在底部结构生成相应质量矿堆，实现该拉底步显式崩落操作后，再次平衡模型。随后进入第N+1步拉底模拟，以此周而复始，最终实现各拉底步的动态崩落。

2 应用研究

2.1 模型构建及地应力反演

依据西南某铜矿地形等高线构建X方向长1 000 m，厚度2 m，地形复杂的二维数值模型，并根据井巷工程布置构建拉底层、聚矿槽、出矿进路等地下工程，最终模型见图3。

由于研究内容为地下开挖问题，模拟选用的本构模型为摩尔-库伦模型。显式崩落模拟需岩体参数、矿堆参数、接触面参数，其中岩体参数采用胡克-布朗(Hoek-Brown)强度折减准则并利用RocLab软件，通过输入完整岩体参数，选择符合矿山实际的岩体条件GSI(地质强度指标)、mi(完整岩石参数)、D(扰动因数)最终确定；矿堆为散体，故强度参数设置为0，弹性参数参照文献[13-14]所提供曲线，根据本次模拟特征近似确定；接触面按照常规取值[15]，具体参数见表1～表3。

图3 矿山模拟模型Fig.3 The model of mine

表1 矿山岩体参数Table 1 The mechanical parameters of mine rock mass

表2 矿堆岩体参数Table 2 The mechanical parameters of ore heap rock mass

表3 接触面参数Table 3 The mechanical parameters of interface

该铜矿地应力的测量主要集中在进路水平，将模型X与Y方向滚支固定，底端完全固定，地表作为自由面，施加初始应力，进行地应力反演。最终进路水平地应力反演值与实际监测值对比见表4。

表4 初始地应力反演结果Table 4 The result of back analysis of origin geostatic stress

模拟过程需人为设定若干参数，且模拟采用的出矿模式为相对出矿量，所需参数见表5。

表5 人工设定参数Table 5 Manual setting parameter

2.2 崩落区及矿堆演化特征

平衡完毕后的模型按图3模拟开挖：开挖始自拉底层中段，跨度15 m，随后每步开挖均向拉底层两侧各推进15 m，聚矿槽与出矿进路滞后拉底推进线30 m开挖，各步骤开挖后运行前均按图2所示流程处理前一步模拟结果，实现显式崩落。共模拟15步，各步模拟结果汇总见图4～图6。

图4 崩落区演化特征Fig.4 Evolution characteristics of caving area

图4中，随拉底向两侧推进，崩落区整体呈穹隆状向上、向四周延伸，该现象与工程实际经验相符；拉底第13步时崩透地表，随后的拉底过程中地表出现大面积破坏，基于地表散体间的摩擦效果，将与水平呈45°范围内的破坏区作为最终崩落区。

图5 矿堆形态演化特征Fig.5 Evolution characteristics of ore heap

各开挖步骤矿堆堆积模拟效果见图5，矿堆于崩透地表前所呈形态为顶端带有一定自然安息角的弧形，一定高度以下崩落区完全填充，同时崩落下的矿堆未侵占未崩落单元。当第13步崩透地表后，引起地表大范围破坏，图4中第13步开挖对应的崩落区将按照与地表轮廓平行的状态进行填充，最终实现地表下沉效果。因此，通过图2所示的显式崩落流程可近似实现矿体的崩落下降效果和矿堆在底部结构的堆积效果。

图6 体积对比统计Fig.6 The plot of volume contrast statistics

图6为各项体积统计图，理论矿堆体积按初始出矿率换算所得，底部空隙体积为累计可填充矿堆体积，实际崩落体积为最终在底部结构实际生成的矿堆体积。工程模拟过程中，由于碎胀系数与出矿量的相对关系，理论矿堆体积始终大于崩落体积；当底部空隙体积大于理论矿堆体积时，崩落下的散体矿石可按计划出矿，此时理论矿堆体积与实际矿堆体积相当；当模拟至15步时，底部空隙体积已小于理论矿堆体积，此时认为矿堆已将底部填满，故将第15步的崩落矿量在原有出矿率基础上多出矿10%，所以实际矿堆体积在15步中较理论矿堆小。由图6可知所设计流程能够自动实现对矿堆是否完全填充空隙的判断，进而动态调整出矿率。

3 2种方式崩落高度与体积对比分析

为更好理解崩落区不做沉降处理的传统崩落模拟方式与显式崩落模拟对崩落过程的不同影响，选用同一模型通过传统崩落模拟方式进行拉底模拟，对比2种方式在崩落高度与崩落体积2方面的差别。

以拉底层上端为基准水平，将每一步崩落区最高值作为崩落高度，图7为2种模拟方式的崩落高度与实际监测顶板高度对比。

由于工程与地质原因，顶板监测数据有若干缺失，但并不影响2种模拟方式对比。显式崩落模拟相较传统模拟，与监测数据吻合度更高，更能反映顶板崩落高度的实际变化趋势。随拉底推进，显式崩落速度明显快于传统崩落，其崩透地表时预测拉底步仅比实际崩透地表晚一步。传统模拟方式顶板崩落高度变化明显滞后于显式崩落和实际监测值。原因是显式崩落移除了每步满足破坏准则的单元，进而促进了崩落区进一步向上发展，所以与传统崩落模拟方式相比，崩落速度更快且更符合实际情况。

图7 崩落高度曲线Fig.7 The plot of caving height

由图8所示，与传统崩落模拟方式相比，显式崩落模拟方法诱发的崩落量更大，原因是在每一步运行时，及时移除了上一步的崩落区单元，促进了与上一步崩落区紧贴的围岩发生更大范围的破坏。因此，从显式崩落的高度与体积均大于传统方式可知：显式崩落模拟能够体现崩落区消除后所产生的空区对下一步崩落的促进效果，该效果有助于提高对崩落预测的准确性。

图8 2种模拟方式崩落体积Fig.8 The plot of caving volume of two simulation methods

4 结论

1)设计并优化了显式崩落算法流程并利用FLAC3D内嵌的编程语言，实现了自然崩落法中崩落体与围岩的脱离，并完成了矿堆在矿柱上方的负载施加，有利于分析底部结构的稳定性。

2)实现了崩落区形态的自动识别，并引入矿堆碎胀系数、出矿量、自然安息角等变量，便于对矿山工程实际问题进行特定分析。

3)设计流程能自动监测矿堆体积与底部空隙体积的相对关系，进而实现出矿量的动态调整。

4)传统与显式2种模拟结果与监测数据对比，突出了崩落区消除后的空区对下一步崩落的促进效果和显式崩落模拟的准确性。

5)虽然该优化算法用于解决自然崩落法模拟过程中存在的问题，但其思路也适用于煤矿开挖等与开挖崩落相关的工程过程的模拟。

[1] WOO K S, EBERHARDT E, ELMOl D, et al. Empirical investigation and characterization of surface subsidence related to block cave mining[J]. International Journal of Rock Mechanics & Mining Sciences, 2013, 61(61):31-42.

[2] ELMO D, STEAD D. An Integrated Numerical Modelling-Discrete Fracture Network Approach Applied to the Characterisation of Rock Mass Strength of Naturally Fractured Pillars[J]. Rock Mechanics & Rock Engineering, 2010, 43(1):3-19.

[3] 赵文. 岩石力学[M]. 长沙：中南大学出版社, 2014:187-188.

[4] BRADVB H G, BROWNE T. 地下采矿岩石力学[M]. 佘诗刚, 朱万成, 赵文，等译. 北京: 科学出版社, 2011.

[5] 薛东杰, 周宏伟, 任伟光, 等. 浅埋煤层超大采高开采柱式崩塌模型及失稳[J]. 煤炭学报, 2015, 40(4):760-765.

XUE Dongjie, ZHOU Hongwei, REN Weiguang, et al. Instability of pillar collapse model and generation of cracks in thick coal seam mining at shallow depth[J]. Journal of China Coal Society, 2015,40(4):760-765.

[6] 王涛, 盛谦, 熊将. 基于颗粒流方法自然崩落法数值模拟研究[J]. 岩石力学与工程学报, 2007, 26(S2):4202-4207.

WANG Tao，SHENG Qian，XIONG Jiang．Research on numerical simulation based on particle flow method[J]．Chinese Journal of Rock Mechanics and Engineering，2007,26(S2):4202-4207.

[7] 王连庆, 高谦, 王建国, 等. 自然崩落采矿法的颗粒流数值模拟[J]. 北京科技大学学报, 2007, 29(6):557-561.

WANG Lianqing, GAO Qian, WANG Jiaguo, et al. Particle flow code numerical simulation of natural caving mining method[J]. University of Science & Technology Beijing, 2007, 29(6):557-561.

[8] 孙闯, 张向东, 张涛,等. 深部大跨度泥质顶板剪切冒落失稳区预测研究[J]. 中国安全生产科学技术, 2016, 12(1):23-27.

SUN Chuang, ZHANG Xiangdong, ZHANG Tao, et al. Prediction on shear caving instability region of deep and large span muddy roof[J]. Journal of Safety Science and Technology, 2016,12(1):23-27.

[9] VYAZMENSKY A, ELMO D, STEA D. Role of Rock Mass Fabric and Faulting in the Development of Block Caving Induced Surface Subsidence[J]. Rock Mechanics & Rock Engineering, 2010, 43(5):533-556.

[10] WOO K S, EBERHARDT E, RABUS B, et al. Integration of field characterisation, mine production and InSAR monitoring data to constrain and calibrate 3-D numerical modelling of block caving-induced subsidence[J]. International Journal of Rock Mechanics & Mining Sciences, 2012, 53(9):166-178.

[11] 曾庆田,刘科伟,严体,等.基于多数值模拟方法联合的自然崩落法开采研究[J].黄金科学技术,2015,23(1):66-73.

ZENG Qingtian, LIU Kewei, YAN Ti, et al. Study on natural caving mining method based on multi-numerical simulation method[J]. Gold Science and Technology, 2015, 23(1): 66-73.

[12] 王涛, 韩煊, 赵先宇, 等. FLAC3D数值模拟方法及工程应用[M]. 中国建筑工业出版社, 2015:170-171.

[13] 陆永青, 计三有, 戴诗亮,等. 散体材料弹性参数的估计[J]. 武汉理工大学学报(交通科学与工程版), 1997(2):178-185.

LU Yongqing, JI Sanyou, DAI Shiliang, et al. Estimate of elastic parameters for granular material[J]. Journal of Wuhan Transportation University, 1997, 21(2): 178-185.

[14] 贾林刚,刘卓然.充填开采充填率与地表移动规律的数值模拟研究[J].辽宁工程技术大学学报(自然科学版),2015,34(9):1010-1015.

JIA Lingang，LIU Zhuoran．Study on numerical simulation of filling rate and surface subsidence in backfill mining[J]．Journal of Liaoning Technical University (Traffic Science and Engineering Edition)，2015，34(9)：1010-1015.

[15] 曹帅,杜翠凤,母昌平,等.崩落法转充填法采矿地表移动二维离散元程序数值模拟及其规律验证[J].岩土力学,2015,36(6):1737-1751.

CAO Shuai，DU Cuifeng，MU Changping，et al．UDEC based modelling of mining surface movement due to transforming from block caving to sublevel filling and its law verification[J]．Rock and Soil Mechanics, 2015, 35(6) :1737-1751.