, 陈佳耀, 周健, 邱贤阳中南大学 资源与安全工程学院，湖南 长沙 410083
收稿日期：2015-12-04
基金项目：“十二五”国家科技支撑计划项目 (2013BAB02B05)。
作者简介：史秀志 (1966-)，男，河北邢台人，中南大学教授，博士生导师。
摘要：为揭示露天坑干堆排尾堆高 (H) 对采空区跨度 (D) 与境界顶柱高 (h) 的影响，采用FLAC3D对不同跨度的双空区在尾砂回填过程中进行开挖模拟，并结合曲线拟合、传统破坏判据及厚度折减法，分别对地下双空区的受力、塑性区及关键点位移等特征进行分析及数学预测.研究表明：当D不变时，拉、压应力及塑性区面积均随H的增大而增加，且在D＞20 m后位移变化率明显增大；当H不变时，监测位移和塑性区随D增大而增加，且定量预测出D=26 m时，上下空区临界破坏对应的H分别为182.44，189.85 m.当D与H均变化时，双空区的塑性区逐渐呈现“蝴蝶状”，上部空区的塑性区域明显大于下部，且在D＞20 m，H＞150 m时，上部空区出现与露天坑底贯穿现象.
关键词：露天坑排尾双空区境界顶柱数值模拟数学预测
Mathematical Prediction on Effect of Tailing Discharge of Open-Pit on Stability of Underground Mined Out Area
SHI Xiu-zhi
, CHEN Jia-yao, ZHOU Jian, QIU Xian-yangSchool of Resources and Safety Engineering, Central South University, Changsha 410083, China
Corresponding author: SHI Xiu-zhi, E-mail: csublasting@163.com
Abstract: To reveal the impact of the height of the tailing dry stacked gangue (H) on the span of double gob area (D) and the roof thickness (h), the simulation software FLAC3D was used to simulate the excavation process of the double gob area of different spans during the tailing-backfilling process. The stress, plastic zones and the displacement of monitoring points of double gob area were analyzed by combining with the hyperbolic fitting equations, traditional failure criterion and the thickness reduction methods. When D is kept constant, tensile and compressive stress and the plastic zone area increase with the increase of H, and the displacement of all monitoring points increase significantly when the span is larger than 20 m. When H is kept constant, the displacement and the plastic zone area increase with the increase of D, and it can be quantitatively predicted that when the span is up to 26 m, the critical damage height of tailing dry stacked gangue is 182.44 and 189.85 m, respectively. When D and H are changed, the plastic zones of the double gob area take on a butterfly-shape, and the plastic zone area of the upper gob is greater than that of the lower part. A cutting-through failure will occur when the span exceeds 20 m and the tailing height exceeds 150 m.
Key Words: tailing discharge into open-pitdouble gob areaboundary pillarnumerical simulationmathematical prediction
露天转地下开采矿山的地下采空区稳定性一直是矿山生产和安全工程的一项重要研究内容.近年来，许多露天转地下开采的矿山利用露天坑自身提供的空间，采用尾砂压滤干堆技术进行采坑的回填，该技术也是金属矿山应对尾矿库容积不足、尾矿库建设成本大等问题提出的切实可行、安全环保的新方案.这虽然对矿山绿化和土地再利用提供了有效途径，但同时也给地下采空区和预留的境界顶柱带来未知的负担，复杂双空区的稳定性问题更是难以全面预测，随着回采的深部延伸，采场顶板的稳定性成为威胁生产连续性的巨大障碍.故对干堆排尾下复杂采空区的深部开采矿体的稳定性进行研究，以及对采空区顶板的预测和趋势评价具有重要的工程实际意义.
国内外专家针对露天转地下回采的稳定性问题做了很多相关研究.其中，林杭等[1]借鉴强度折减法对采空区临界安全顶板进行预测；吴启红等[2]采用专家评定法和三相线性隶属函数对各种类型变量评价并建立采空区稳定性二维模糊评判模型；张耀平等[3]采用FLAC3D并结合现代仿真技术研究空区稳定性；叶振华[4]利用有限元程序对露天坑排放尾砂后的变形和渗流进行模拟研究；王新民等[5]利用层次分析和模糊数学法研究露天转地下矿山的最佳开采模式.然而，露天坑排尾条件下的采空区稳定性是一个受多因素影响、伴随时空变异的复杂动态系统，故通过简单的理论分析很难预估合理的采场参数和安全经济的境界顶柱厚度，同时大多数分析方法对复杂情况下采空区沉降量所做的监测工作较少，据统计，矿山由于选择不合理参数或错误预测所带来的经济损失较大.因此，实现多空区稳定性数据的合理评估对于回采的安全可持续影响重大.
本文就某金属矿利用自身露天坑干堆尾砂的情况，结合其地下采场双层空区的复杂情况，借助FLAC3D建立三维地质力学模型进行空区稳定性分析，对不同跨度的双空区在尾砂回填过程中的情况进行开挖模拟.通过监测关键点位移、塑性区特征、应力变形等，结合相关曲线特征拟合方程、传统破坏判据及厚度折减法，预测并分析尾砂堆高、采空区跨度和空区稳定性的关系，并在实际工程中加以应用，为复杂空区的稳定性问题提供一种新思路.
1 工程概况某金属矿目前已由露天转地下开采，预计在未来5年内对采坑进行尾砂回填，回填时间跨度较大，预计堆高最高可达200 m.矿体周边分布2种低强度的岩体，中间含有较倾斜的大理岩夹层，一直延伸至露天坑口，平均厚度8 m，见图 1.该矿利用上向水平分层嗣后充填法进行回采，先采2层并胶结充填1层，层高3 m，回采过程中最大空区高达6 m，夹石上下均含有开采的空区；定义双空区跨度为D，尾砂堆高为H，境界顶柱高度为h.露天坑口已形成大面积采坑空间，并集聚大量积水，回填的尾砂与水混合形成低强度泥沙，造成复杂的地质构造.
图 1(Fig. 1)
 | 图 1 矿体剖面示意图及监测点布置 (单位：m)Fig.1 Ore body profile sketch map and layout of monitoring points |
2 数值模拟与计算分析2.1 建模根据矿山赋存条件和工程地质调查，获取相关岩体物理力学参数，如表 1所示.结合图 1中工程示意，建立1:1比例模型，合计4 821个单元，9 750个节点.
表 1(Table 1)
 表 1 岩体物理力学参数Table 1 Mechanical parameters of rock mass
| 表 1 岩体物理力学参数 Table 1 Mechanical parameters of rock mass |
工程实践表明顶板中心和两端是最易破坏点，故在上下2个采空区分别布置3个监测点，并分别命名为1#，2#，3#，4#，5#，6#监测点，其中1#，4#监测点分别布置于空区顶板中心位置，其余各点分别布置于空区两端.利用自带elastic命令获取初始应力，null命令开挖实体模型形成2个高度为6 m、跨度为D的模型矿体，空区垂直间隔距离为8 m，对尾砂堆积以5 m为单位进行堆高模拟，同时监测各点的位移、塑性区变化，利用history write功能进行记录，并做定量分析.
模拟分析分为2步： ① 利用函数拟合法得到系列模拟数据，包括空区跨度、境界顶柱高度等. ② 将 ① 中不同组数据进行逐层堆积尾砂模拟，通过监测各点位移、塑性区面积、应力值等数据绘制曲线，预测破坏情况；同时，利用自带FISH语言编写的程序，定义破坏判据值，通过判据分析单元体的破坏状态，研究尾砂堆高过程对空区影响，观察空区破坏情况，综合模拟结果分析围岩变形和破坏情况[6-7].
2.2 分析与拟合为进一步探究安全境界顶柱高度与空区跨度间的函数关系，选取合理的模拟试验值，结合文献[1]中的强度折减参数拟定多组临界安全顶柱和空区跨度数据，选择具有代表性的数组作为函数样本，绘制曲线见图 2.其中，模拟过程保证较大安全系数和研究范围，并进行拟合处理.
图 2(Fig. 2)
 | 图 2 境界顶柱高度与空区跨度关系Fig.2 Relation between span and roof thickness of gob |
拟合曲线采用二次多项式：
 | (1) |
史秀志等[8]结合数值模拟法，对尾砂填充下的安全顶板高度和空区跨度进行位移应力变形等分析，得到两者之间的一系列数组，本文根据文献[1, 8]提供的有效模拟数据，进行某安全系数下的厚度折减调整，作为境界顶柱高度的选择标准，通过设置最大不平衡率R=5×10-7判定系统平衡状态，通过拟定曲线进行破坏状态的判定，即拟合尾砂堆积过程中两者安全表达式为
 | (2) |
选取一系列不同跨度、顶板高度作为模拟数据，利用上述拟合二次型函数，其中空区跨度为8，11，14，17，20，23，26 m时，得到对应的安全境界顶柱高分别为7，9，12，15，18，21，24 m，以上7组数据将作为模拟第一步的参数.
Duncan[6]根据Konder建议采用曲线方程拟合土体应变-应力曲线，取得了较好的效果.类似地，本文利用空区位移和尾砂堆高之间的关系，通过拟合曲线方程[7-11]，根据曲线的突变特征，假设空区监测点位移与尾砂堆高满足泰勒级数曲线的形式：
 | (3) |
令
 | (4) |
 | (5) |
 | (6) |
 | (7) |
 | (8) |
 | (9) |
当Δ>0时，空区处于稳定状态；Δ=0时，空区处于临界失稳状态；当Δ＜0时，空区处于失稳破坏状态；故Δ=0所对应的H值为“破坏点”.
3 模拟结果与分析模拟过程中通过监测应力应变观察空区周边破坏情况，由等值云图了解尾砂堆积影响范围及对周边和地表建筑的作用.通过记录塑性区分布、关键点位移等来分析采空区变形情况[12-19].另外，结合统计的塑性区面积拟合得到的关键点位移和影响因子之间的相关关系，对尾砂逐渐堆积的过程进行系统预测研究.
3.1 应力分析不同D的双空区在同一尾砂堆高 (H=160 m) 下的最大主应力分布情况见图 3.由于空区位于矿体中部，最大主应力在各跨度空区周边呈近对称分布.当D较小时，上部空区顶板中心位置附近出现了较小的压应力区，空区两帮出现剪切应力，底板出现较小压应力；下部空区应力分布较为简单，空区顶板和底板位置出现拉应力，两帮出现较小的剪切应力.随着D的增加，应力分布区域发生较大变化，上部空区受到的压应力范围逐渐增大，顶板中心出现了局部拉应力集中现象，在D=26 m的空区顶板处拉应力达0.542 MPa，上部空区两帮受到剪切应力，D=26 m时剪切应力达0.684 MPa；下部空区顶底板拉应力值和压应力的范围均变大，且在下部空区底板中心位置出现了拉应力集中现象.在任意空区跨度条件下，下部空区受到的拉应力最大值均小于上部空区拉应力最大值，同时影响范围也不及下部空区.跨度为8，14，20，26 m对应的最大拉应力分别为0.545，0.941，1.278，1.452 MPa，有明显的增大趋势.
图 3(Fig. 3)
 | 图 3 不同跨度下最大主应力分布图Fig.3 Maximum stress distribution under different spans (a)—D=8 m；(b)—D=14 m; (c)—D=20 m；(d)—D=26 m. |
同一跨度的双空区 (D=23 m) 在不同尾砂堆高H条件下的最大主应力分布见图 4.H不同时，最大应力云图也呈近对称分布.当H=0 m时，上下部空区均有小范围的拉应力区，拉应力值为0.076 MPa，两帮有较小的压应力.随着H的增大，各应力分布区域发生较大变化，上部空区的拉应力增大，且最大拉应力出现在顶板中心位置附近，两帮受到的剪切应力逐渐变大；下部空区受到拉应力增大，且下部空区底板位置出现受力最大值.尾砂堆高H分别为0，50，100，150 m时，空区的最大应力值变化明显，分别为0.075，0.081，0.278，0.452 MPa.根据应力数据的统计分析知，下部空区受到的拉应力比上部空区的小.随着H的增大，双层空区受到的最大应力均明显增大.
图 4(Fig. 4)
 | 图 4 不同尾砂堆高下最大主应力分布图Fig.4 Maximum stress distribution under different height of tailings pile (a)—H=0 m；(b)—H=50 m; (c)—H=100 m；(d)—H=150 m. |
3.2 位移分析不同跨度双层空区在不同尾砂堆高H下的监测点位移情况见图 5，模拟过程中对各空区的顶板及两帮分别进行记录.可得，当D小于空区高度，即D＜6m时，上下空区的位移变化小且平稳，尾砂堆积过程中，空区稳定性较好，最大位移为8.41 mm，出现在上部空区顶板位置；当D>6 m时，各空区的位移变化速率明显增加.随着D的增加，各监测点位移有不同程度的上升趋势，且D＜20 m时，各位移变化速率相对较小，D>20 m时，各位移变化率明显增大，可以定性地分析出在H=180 m，D=26 m的上部空区出现了垮塌现象，在H=190 m，D=26 m的下部空区也出现了垮塌现象，且下部空区的位移量小于上部空区.模拟过程中检测的数据常用于定性分析，并不能对数据进行正确的估计，故需要进一步的函数定量分析进行优化.
图 5(Fig. 5)
 | 图 5 不同尾砂堆高条件下监测点位移图Fig.5 Displacements of monitoring points under different height of tailings pile |
表 2是在图 5监测点位移的基础上，通过式 (3) 提供的拟合曲线函数方程D(H)=a0+a1(H)+a2(H2)+a3(H3)+a4(H4)，定量计算尾砂堆积过程中的破坏情况.结果显示，在D＜26 m时，拟合的平均相关系数R2=0.997 4，且判据Δ≠0，即在H＜200 m时，式 (8) 无解，空区可以基本判断为未发生破坏；在D≥26 m时，由于H>180 m出现了较大位移值，故拟合程度出现了轻微的变动，空区出现了破坏现象，根据破坏判据Δ=0，得出在H=182.44 m，D=26 m时的上部空区出现了破坏，位移值为3 254 mm；当H=189.85 m，D=26 m时的下部空区出现了破坏，位移值为2 564 mm.表 2的函数结果能较为准确地计算出动态系统中某值对应的稳定情况，能为综合的位移监测提供相关途径.
表 2(Table 2)
表 2 各顶板中心监测点在不同尾砂堆高下曲线拟合结果Table 2 Fitting results of monitoring points of upper gob under different height of tailings pile
| 表 2 各顶板中心监测点在不同尾砂堆高下曲线拟合结果 Table 2 Fitting results of monitoring points of upper gob under different height of tailings pile |
3.3 塑性区及空区破坏特征分析不同D的空区在同一尾砂堆高 (H=150 m) 条件下的塑性区分布情况见图 6.在空区高度和空区间距相同的条件下，D较小时，上部空区主要发生剪切破坏，中间夹层出现极小塑性区分布，下部空区两帮受到拉剪破坏；随着D的增大，空区周边分布的塑性区面积明显增大，渐形成“蝴蝶状”塑性区分布，且当D>20 m时，破坏形式由拉剪破坏变为剪切破坏，且上下空区塑性区开始在两帮位置贯通，并开始延伸至境界顶柱上部，上部空区有破坏的危险.当D较大时，见图 6d，空区间的中间夹层没有过大的塑性区，上部空区面临着垮塌危险，但间柱和下部空区则属于相对“安全”的位置，此时间柱对下部硐室的“保护”作用明显.
图 6(Fig. 6)
 | 图 6 不同跨度条件下塑性区分布图Fig.6 Plastic zone distributions under different spans (a)—D=8 m；(b)—D=14 m; (c)—D=20 m；(d)—D=26 m. |
同一跨度 (D=20 m) 空区在不同尾砂堆高H情况下塑性区分布见图 7.尾砂堆高H=0 m时，上下空区均为很小的塑性区面积，稳定性好.随着尾砂H的增大，塑性区面积逐渐增大，堆高H在0~50 m范围内塑性区发展较平缓；H>50 m后，塑性区面积出现突变，即堆高H对空区的稳定性影响明显增大，下部空区的塑性区面积始终小于上部空区，且呈现明显的“蝴蝶状”.当H>150 m时，上下空区塑性区局部贯通，两空区均有破坏危险，下部空区主要破坏形式为两帮拉剪破坏，而上部空区则为两帮剪切破坏，上部破坏要先于下部空区，空区间夹层塑性区范围小，在上部空区破坏后，中间夹层对空区有相对“保护”作用，夹层和空区均处于相对“安全”状态.
图 7(Fig. 7)
 | 图 7 不同尾砂堆高条件下最大主应力分布图Fig.7 Distribution of maximum principal stress under different heights of tailings pile (a)—H=0 m；(b)—H=50 m; (c)—H=100 m；(d)—H=150 m. |
4 结论1) 尾砂堆高H相同，D较小时，上部空区顶板主要分布压应力，两帮出现剪切应力，下部空区出现较小拉应力.随着D的增大，上部空区由压应力变为拉应力集中区，分布范围较广，上下部空区两帮均受到不同大小的剪切应力，下部空区底板出现拉应力集中区.随着H的增大，双层空区应力均有不同程度增大的趋势，且速率不同，尾砂堆积越高，对空区稳定性影响越大，对周边的影响范围越广.
2) 跨度D小于空区高度时，双层空区位移值均较小，变化平缓，随着D的增加，空区各点位移增大；当D，H相同时，上部空区顶板中心监测点位移值始终大于下部空区位移值，且满足拟合曲线函数，通过破坏判据确定下部空区破坏滞后于上部空区，且定量预测出跨度D=26 m时，上下空区破坏对应的H分别为182.44，189.85 m.
3) D相同时，塑性区面积随H的增大而增大，且H>50 m后增大明显，H>150 m时会出现同境界顶柱局部贯通，对下部空区稳定性产生影响；H相同时，塑性区随D增大而增大，且在跨度D>20 m时出现空区贯通现象，上下空区塑性区呈现“蝴蝶状”，空区间隔处于相对“安全”.
参考文献
| [1] | 林杭, 曹平, 李江腾, 等. 采空区临界安全顶板预测的厚度折减法[J].煤炭学报, 2009, 34(1): 53–57. ( Lin Hang, Cao Ping, Li Jiang-teng, et al. The thickness reduction method in forecasting the critical safety roof thickness of gob area[J].Journal of China Coal Society, 2009, 34(1): 53–57.) |
| [2] | 吴启红, 彭振斌, 陈科平, 等. 矿山采空区稳定性二级模糊综合评判[J].中南大学学报 (自然科学版), 2010, 41(2): 661–667. ( Wu Qi-hong, Peng Zhen-bin, Chen Ke-ping, et al. Synthetic judgment on two-stage fuzzy of stability of mine gob area[J].Journal of Central South University (Science and Technology), 2010, 41(2): 661–667.) |
| [3] | 张耀平, 曹平, 袁海平, 等. 复杂采空区稳定性数值模拟分析[J].采矿与安全工程学报, 2010, 27(2): 233–238. ( Zhang Yao-ping, Cao Ping, Yuan Hai-ping, et al. Numerical simulation on stability of complicated goaf[J].Journal of Mining & Safety Engineering, 2010, 27(2): 233–238.) |
| [4] | 叶振华. 露天采场排放尾矿后下部矿体开采安全性研究[J].矿冶, 2006, 15(4): 9-12–26. ( Ye Zhen-hua. Study on safety of lower deposit mining under the condition of discharging tailings in the upper strip pit[J].Mining & Metallurgy, 2006, 15(4): 9-12–26.) |
| [5] | 王新民, 赵建文, 张钦礼, 等. 露天转地下最佳开采模式[J].中南大学学报 (自然科学版), 2012, 43(4): 1434–1439. ( Wang Xin-min, Zhao Jian-wen, Zhang Qin-li, et al. Mining model of transition from open-pit to underground mining[J].Journal of Central South University (Science and Technology), 2012, 43(4): 1434–1439.) |
| [6] | Duncan J M. State of the art limit equilibrium and finite element analysis of slopes[J].Journal of Geotechnical Engineering, ASCE, 1996, 122(7): 577–596.DOI:10.1061/(ASCE)0733-9410(1996)122:7(577) |
| [7] | Zhang H R. Researches and applications on geostatistical simulation and laboratory modeling of mine ventilation network and gas drainage zone[J].Process Safety and Environmental Protection, 2015, 94: 55–64.DOI:10.1016/j.psep.2014.10.003 |
| [8] | 史秀志, 黄刚海, 张舒, 等. 基于FLAC3D的复杂条件下露天转地下开采空区围岩变形及破坏特征[J].中南大学学报 (自然科学版), 2011, 42(6): 1710–1718. ( Shi Xiu-zhi, Huang Gang-hai, Zhang Shu, et al. Goaf surrounding rock deformation and failure features using FLAC3Din underground mining shifted from open-pit in complex situation[J].Journal of Central South University (Science and Technology), 2011, 42(6): 1710–1718.) |
| [9] | Zhou J, Li X B. Identification of large-scale goaf instability in underground mine using particle swarm optimization and support vector machine[J].International Journal of Mining Science and Technology, 2013, 23(5): 701–707.DOI:10.1016/j.ijmst.2013.08.014 |
| [10] | Nomikos P P, Sofianos A I, Tsoutrelis C E. Structural response of vertically multi-jointed roof rock beams[J].International Journal of Rock Mechanics and Mining Sciences, 2002, 39(1): 79–94.DOI:10.1016/S1365-1609(02)00019-9 |
| [11] | Yu Z Y, Yang S Q, Qin Y, et al. Experimental study on the goaf flow field of the "U+I"type ventilation system for a comprehensive mechanized mining face[J].International Journal of Mining Science and Technology, 2015, 25(6): 1003–1010.DOI:10.1016/j.ijmst.2015.09.019 |
| [12] | Deng Q W, Liu X H, Lu C, et al. Numerical simulation of spontaneous oxidation zone distribution in goaf under gas stereo drainage[J].Procedia Engineering, 2013, 52: 72–78.DOI:10.1016/j.proeng.2013.02.108 |
| [13] | Hao T X, Jin Z C, Li F. Optimization of goaf gas drainage parameters based on numerical simulation studying fracture in overlying strata[J].Procedia Engineering, 2012, 43: 269–275.DOI:10.1016/j.proeng.2012.08.046 |
| [14] | Qin Z Y, Yuan L, Guo H, et al. Investigation of longwall goaf gas flows and borehole drainage performance by CFD simulation[J].International Journal of Coal Geology, 2015, 150/151: 51–63.DOI:10.1016/j.coal.2015.08.007 |
| [15] | Xu P, Mao X B, Zhang M X, et al. Safety analysis of building foundations over old goaf under additional stress from building load and seismic actions[J].International Journal of Mining Science and Technology, 2014, 24(5): 713–718.DOI:10.1016/j.ijmst.2014.03.030 |
| [16] | Qin Y P, Liu W, Yang W W, et al. Numerical simulation study of spontaneous combustion in goaf based on non-Darcy seepage[J].Procedia Engineering, 2011, 26: 486–494.DOI:10.1016/j.proeng.2011.11.2196 |
| [17] | Li Z X, Lu Z L, Wu Q, et al. Numerical simulation study of goaf methane drainage and spontaneous combustion coupling[J].Journal of China University of Mining and Technology, 2007, 17(4): 503–507.DOI:10.1016/S1006-1266(07)60134-5 |
| [18] | Li Z X, Huang Z A, Zhang A R, et al. Numerical analysis of gas emission rule from a goaf of tailing roadway[J].Journal of China University of Mining and Technology, 2008, 18(2): 164–167.DOI:10.1016/S1006-1266(08)60035-8 |
| [19] | He X G, Zhang R W, Pei X D, et al. Numerical simulation for determining three zones in the goaf at a fully-mechanized coal face[J].Journal of China University of Mining and Technology, 2008, 18(2): 199–203.DOI:10.1016/S1006-1266(08)60042-5 |
