王伟

摘要：为研究复杂采空区的稳定性，开展工程地质调查，利用RMR分级法和Q系统分级法对各类型岩体进行评价分级，采用Flac3D软件对采空区形成过程及采空区稳定性进行模拟计算和预测分析。研究结果表明：上盘围岩与矿体岩体质量均良好，下盘围岩岩体质量中等。通过应力等值线云图、位移等值线云图和塑性区共同分析得出：方案1中2号和3号采空区的围岩由于剪切应力过于集中，剪切应力大于岩体抗剪能力，使围岩进入了塑性状态;最大位移出现在方案2的2号采空区顶板。水平位移最大值出现在采空区的顶面及底面中心处，且有贯通的趋势。方案2中塑性区虽未贯通，但采空区已有塑性区将对采空区的稳定造成一定安全隐患。

关键词：采空区;数值模拟;稳定性;Q系统分级法;RMR分级法

中图分类号：TD32          文章编号：1001-1277（2021）06-0041-06

文献标志码：Adoi：10.11792/hj20210608

引 言

近年来由于房柱采矿法、全面采矿法及留矿采矿法等空场类方法广泛应用于采矿工程中[1]，形成了大量的地下采空区，影响矿山的安全生产，特别是一些大型复杂采空区群条件下的岩层。因此，开展地下采空区稳定性研究，对保证矿山回采安全具有重要的意义。

专家学者对采空区的稳定性分析进行了大量的研究，周宗红等[2]基于平衡拱理论，对跑马坪铅锌矿采空区顶板临界冒落面积进行计算分析，并采用三维有限元数值模拟方法，分析采空区顶板岩层变形破坏机制，提出相应控制方法。张耀平等[1]利用Flac3D软件对龙桥铁矿采空区形成过程及采空区稳定性进行模拟计算和预测分析，并对后续开采提出了建议。李想等[3]利用Flac3D软件模拟某铁矿中深部矿房回采后的采空区应力分布情况及下阶段回采过程对上阶段采空区的影响，得出采空区的稳定性状况及下阶段回采对上阶段采空区的影响较小的结论，并预测了采空区的发展。对地下采空区进行精确探测技术的研究和采空区稳定性安全评价技术的研究显得尤为重要，国内外很多学者在该方面取得了很多富有成效的研究成果[4-15]。

某铁矿区的工程地质条件为基岩块状坚硬—半坚硬岩区和第四系碎粒散体松散岩区，矿区内大面积分布基岩块状坚硬—半坚硬岩区，主要岩性为安山质凝灰岩，出露厚度1 000 m左右。安山质凝灰岩属于侏罗系，节理、裂隙发育，每米一般可见2～5条，密集段可见5～10条，裂隙以闭合为主，裂隙宽0.1～8.0 mm，矿体大部分围岩及顶底板均为安山质凝灰岩。第四系碎粒散体松散岩区远离矿体，对矿床开采影响不大。

开采过程中，地表形成9个露天采坑，井下形成30个采空区。在30个井下采空区中，Fe1号矿体形成11个采空区，Fe2号矿体形成13个采空区，Fe3号矿体形成6个采空区。采空区最长442 m，最宽21.3 m，采高最高18.8 m，地下采空区中无积水。已有采空区的稳定直接关系到矿山深部开采安全，一旦发生大范围采空区冒落垮塌事故，势必会造成严重的经济损失，以及人员伤亡和设备损坏。为此，亟待开展矿山地下采空区稳定性研究工作，掌握采空区的应力场、位移场和塑性区分布情况，为矿山安全高效回采奠定基础。

1 矿（岩）体质量分级

1.1 采场围岩结构面调查及分析

1.1.1 结构面调查

按照测线法和体积密度法，使用罗盘和皮尺对影响采空区稳定性的矿（岩）体中的节理、裂隙产状、规模、密度、形态、地下水状况等内容进行现场调查，共调查了8条测线。结构面调查测线分布情况如表1所示。

1.1.2 结果及分析

1）结构面极点等密图。利用罗盘和皮尺，调查统计测线周围的结构面倾向、倾角、迹长、类型、粗糙度、淋水情况、风化程度等影响岩体稳定性的因素，采用Dips软件，分别对上盘围岩、矿体、下盘围岩3种岩体进行统计分析，绘制3种岩体结构面的极点等密图、优势结构面赤平投影图，如图1～3所示。

根据调查结果及Dip软件分析得到不同优势结构面产状，如表2所示。从表2可以看出，下盘围岩优势结构面较多，岩体完整性较差。

2）RQD估算值。根据调查结果，利用Farmer密度估计法、体积节理统计估算法、精测线估算法3种方法对不同类型矿（岩）体进行RQD值估算。其中，Farmer密度估计法得到的RQD值为范围值，体积节理统计估算法及精测线估算法得到的RQD值为具体值。首先根据后2种方法得到平均值，再和第1种方法进行比较。3种方法得到的RQD值如表3所示。

从表3可以看出：上盘围岩RQD综合取值为93.1 %，岩体质量很好;矿体RQD综合取值为99.5 %，岩体质量很好;下盘围岩RQD综合取值為85.1 %，岩体质量好。从调查结果看，下盘围岩较破碎，优势节理面较多，且部分交叉，切割破碎情况相对严重，调查部分裂隙较多，且多张开，充填物较多，软硬不一，表面潮湿，因此应对下盘围岩部分地段岩体质量进行折减。矿体RQD综合取值最高，3条穿脉中节理裂隙不发育，且分布不均匀，除一条较为明显的断层露头外，其余没有特别明显的节理面，岩体完整，岩体质量极好。

1.2 矿（岩）体稳定性分级方法

采用Bieniawski的地质力学RMR分级法和Barton的Q系统分级法对该铁矿地下矿（岩）体的稳定性进行分级评价。根据各岩层参数，RMR分类结果表明：上盘围岩属于Ⅱ类岩石，质量良;矿体属于Ⅱ类岩石，质量良;下盘围岩属于Ⅲ类岩石，质量中等。

Q系统分级法考虑的因素与RMR分级法考虑因素比较接近，但是其得分计算方法却是乘积法，即对六因素进行如下计算[4]：

Q=RQDJn×JrJa×JwSRF（1）

式中：RQD为岩石质量指标;Jn为节理组数系数;Jr为节理粗糙度系数（最不利的不连续面或节理组）;Ja为节理蚀变度（变异）系数（最不利的不连续面或节理组）;Jw为节理渗水折减系数;SRF为应力折减系数。其中RQD与Jn之比可粗略表示岩石的块度;Jr与Ja之比表示嵌合岩块的抗剪强度;而Jw与SRF之比反映了岩石的主动应力。

各类型岩体参数，Q系统分级的分类结果表明：上盘围岩属于Ⅱ类岩石，质量良;矿体属于Ⅱ类岩石，质量良;下盘围岩属于Ⅲ类岩石，质量中等。

2 计算方案确定

根据弹性力学，采空区之间的相互作用一般集中在相邻的采空区之间，若采空区相距较远，则相互之间的影响较小，因此，根据该铁矿的勘探线确定以下2种计算方案：

1）方案1。Ⅴ号勘探线所在剖面上存在着4个采空区，即1号采空区、2号采空区、3号采空区和4号采空区。

2）方案2。Ⅱ号勘探线所在剖面上存在着3个采空区，即2号采空区、5号采空区和7号采空区。

模拟采空区布置方案如图4所示。从图4可以看出，这2个方案中采空区数量多，且位置相对较近，若这2个方案中的采空区是稳定的，则其他采空区亦可认定为稳定的。

对以上2个方案进行Flac3D数值模拟，分析采空区之间的相互关系，通过应力场、位移场和塑性区范围判断其断面上采空区之间的稳定性。

3 计算模型与参数确定

3.1 模型建立

由于矿（岩）体形态与地层结构复杂，直接采用Flac3D软件建模不方便，因此采用Abaqus软件提取各剖面图的关键点依次生成线、面、体的方法生成包含矿体和围岩的数值模型，并自动进行网格剖分。通过接口程序，将模型数据文件转化成Flac3D软件可以使用的数据格式，导入Flac3D软件，从而实现Flac3D网格剖分的直观化、形象化和自动化。

根据设计采场围岩受采动影响范围的大小，模型沿东西方向为z方向，沿南北方向为x方向，垂直水平面方向为y方向。分别建立了上盘围岩、下盘围岩和矿体模型，方案1模型网格图如图5-a）所示，整体模型尺寸为370 m×370 m×280 m，从上至下分别是1号采空区、2号采空区、3号采空区和4号采空区。方案2模型网格图如图5-b）所示，整体模型尺寸为370 m×370 m×60 m，从上至下分别是2号采空区、5号采空区和7号采空区。初始模型采用四面体单元剖分，重点剖分矿体采空区部分，精度以确保无畸变单元为原则并在局部适当加密。

本次三维数值计算将岩性简化为围岩和矿体2种介质类型，材料为理想弹塑性，本构模型为Mohr-Coulomb 准则。

3.2 计算参数确定

根据采场围岩岩体质量分级结果，并参考周边矿山矿（岩）体参数，确定计算模型所采用的矿（岩）体物理力学参数如表4所示。

4 模拟结果及分析

4.1 主应力分布规律

2种方案采空区最大主应力场分布如图6、图7所示。从图6和图7可以看出：矿体回采后，各个采空区周围岩体出现条带状应力等值线，且等值线倾向下坡方向。方案1中2号采空区右侧顶部和4号采空区以下出现了最大主应力，最大主应力值为4.4 MPa;方案2中2号采空区、5号采空区和7号采空区之间出现了最大主应力，特别是在5号采空区与7号采空区之间较为明显，最大主应力值为0.5 MPa。最大主应力为负值，表明采空区受到了拉应力，拉应力的出现对采空区稳定不利，但不会导致其产生破坏。

从应力场角度分析：2号和3号采空区的围岩由于剪切应力过于集中，大于岩体自身的剪切应力，使围岩进入了塑性状态，采空区的顶板塑性状态大多也是由于受到较大的剪切应力产生的。这就对采空区的稳定性造成一定影响，且不排除采空区局部地段发生失稳塌陷的可能。因此，随着开采深度的增加，将对采空区的稳定形成一定的威胁，故必须加强采空区现场监测工作，一旦发生异常及时采取有效的控制措施，确保矿山深部开采安全。

4.2 位移变形分布规律

2种方案的采空区位移场分布如图8、图9所示。从图8可以看出：矿体回采后，由于围岩应力转移并重新分布，采空区周围岩体位移呈现条带状等值线云图。1号采空区的顶板产生向下的竖向位移，为0.045 m，表明该位置出现了沉降现象;底板产生向上的竖向位移，为0.015 m，表明该位置出现了轻微的隆起。2号采空区的顶板产生向下的竖向位移，为0.045 m，表明该位置出现了沉降现象;底板产生向上的竖向位移，为0.02 m，表明该位置出现了轻微的隆起。3号采空区的顶板产生向下的竖向位移，为0.015 m，表明该位置出现了沉降现象;底板产生向上的竖向位移，为0.04 m，表明该位置出现了轻微的隆起。4号采空区的顶板产生向下的竖向位移，为0.015 m，表明该位置出现了沉降现象;底板产生向上的竖向位移，为0.03 m，表明该位置出现了轻微的隆起。

从图9可以看出：当矿体回采后，由于围岩应力转移并重新分布，采空區周围岩体位移呈现条带状等值线云图。2号采空区的顶板产生向下的竖向位移，为0.08 m，表明该位置出现了沉降的趋势;底板产生向上的竖向位移，为0.03 m，表明该位置出现轻微的隆起。5号采空区的顶板产生向下的竖向位移，为0.07 m，表明该位置出现了沉降的趋势;底板产生向上的竖向位移，为0.03 m，表明该位置出现轻微的隆起。7号采空区的顶板产生向下的竖向位移，为0.07 m，表明该位置出现了沉降的趋势;底板产生向上的竖向位移，为0.07 m，表明该位置出现轻微的隆起。

从位移场角度分析：由于原有应力平衡状态被矿体的开挖所打破，采空区顶板、底板受高垂直、水平应力所挤压，产生变形位移向采空区移动，在采空区的顶板、底板以外的围岩区域形成位移等值线拱，随着距采空区的距离增大，拱径逐渐变大，而位移逐渐减小。水平位移最大值出现在采空区的顶板及底板中心处。因此，建议矿山加强采空区的现场监测工作，及时采取科学有效的采空区治理措施，为矿山深部安全开采奠定基础。

4.3 塑性区分布规律

2种方案采空区塑性区分布如图10、图11所示。从图10可以看出，塑性区主要出现在2号采空区和3号采空区之间，且有贯通的趋势。该塑性区是由剪切应力引起的，其体积为82.28 m3。

从图11可以看出：塑性区主要出现在2号采空区、5号采空区和7号采空区的顶板和底板，并未出现贯通的趋势，对于这3个采空区来说，虽然未出现贯通，但是顶板和底板都出现了塑性区，对采空区稳定不利。该塑性区是由剪切应力引起的，其体积为323.54 m3，但采空区之间并未有任何影响。

从塑性区角度分析：采空区已有塑性区将对采空区的稳定性造成一定的安全隐患，这对矿山深部矿体的开采极为不利，需要引起矿山工作人员的高度重视，做好现场围岩的监测管理工作。

同时，在矿山应做好深部矿体回采的安全管理工作，如优化爆破参数和回采顺序等，确保矿山采场安全高效的开采。

5 结 论

1）采用RMR分级法和Q系统分级法对该铁矿地下矿岩的稳定性进行分级评价，结果表明：上盘围岩与矿体均属于Ⅱ级，岩体质量良好。下盘围岩属于Ⅲ级，岩体质量中等。参照周边矿体及围岩情况确定数值模拟计算参数。

2）方案1中2号和3号采空区的围岩由于剪切应力过于集中，大于岩体抗剪能力，使围岩进入了塑性状态;采空区顶板塑性状态大多也是由于受到较大的剪切应力产生的。

3）采空区顶板、底板受高垂直、水平应力所挤压，产生变形位移向采空区移动，最大位移出现在方案2中的2號采空区顶板。水平位移最大值出现在采空区的顶板及底板中心处。

4）方案1中塑性区出现在2号采空区和3号采空区之间，且有贯通的趋势。方案2中塑性区虽未贯通，但其采空区已有塑性区将对采空区的稳定性造成一定的安全隐患，这对矿山深部矿体的开采极为不利，需要引起矿山工作人员的高度重视，做好现场围岩的监测管理工作。

[参考文献]

[1] 张耀平，曹平，袁海平，等.复杂采空区稳定性数值模拟分析[J].采矿与安全工程学报，2010，27（2）：233-238.

[2] 周宗红，侯克鹏，任凤玉.跑马坪铅锌矿采空区稳定性分析及控制方法[J].采矿与安全工程学报，2013，30（6）：863-867.

[3] 李想，朱为民，谢晓斌.FLAC3D在某铁矿采空区稳定性分析中的应用[J].中国矿山工程，2012，41（1）：25-29.

[4] 李夕兵，李地元，赵国彦，等.金属矿地下采空区探测、处理与安全评判[J].采矿与安全工程学报，2006，23（1）：24-29.

[5] 宫凤强，李夕兵，董陇军，等.基于未确知测度理论的采空区危险性评价研究[J].岩石力学与工程学报，2008，27（2）：323-330.

[6] 童立元，刘松玉，邱钰，等.高速公路下伏采空区问题国内外研究现状及进展[J].岩石力学与工程学报，2004，23（7）：1 198-1 202.

[7] 刘敦文，褚夫蛟，任高峰，等.硬岩矿深部开采过渡区厚度优选与数值分析[J].东北大学学报（自然科学版），2015，36（4）：571-575.

[8] 李俊平，王晓光，王红星，等.某铁矿采空区处理与卸压开采方案研究[J].安全与环境学报，2015，15（1）：137-141.

[9] 程久龙，胡克峰，王玉和，等.探地雷达探测地下采空区的研究[J].岩土力学，2004（增刊1）：79-82.

[10] 陶帅，王学滨，潘一山，等.基于摩尔-库仑模型的非线性本构模型的开发及其在应变局部化中的应用[J].岩土力学，2011，32（增刊1）：403-407.

[11] 罗周全，谭浪浪，邓俏，等.采空区失稳模式数值分析及实测验证[J].矿业研究与开发，2012，32（3）：89-92.

[12] 王金安，尚新春，刘红，等.采空区坚硬顶板破断机理与灾变塌陷研究[J].煤炭学报，2008，33（8）：850-855.

[13] 吴启红，万世明，彭文祥.一种多层采空区群稳定性的综合评价法[J].中南大学学报（自然科学版），2012，43（6）：2 324-2 330.

[14] 吴爱祥，王贻明，胡国斌.采空区顶板大而积冒落的空气冲击波[J].中国矿业大学学报，2007，36（4）：473-477.

[15] 张敏思，朱万成，侯召松，等.空区顶板安全厚度和临界跨度确定的数值模拟[J].采矿与安全工程学报，2012，29（4）：543-548.

Stability analysis of complex goaf based on Flac3D

Wang Wei

（China National Gold Group Co.，Ltd.）

Abstract：In order to study the stability of complex goaf and carry out engineering geological survey，RMR classification method and Q system classification method are used to evaluate and classify each part of rock mass，and Flac3D numerical simulation software is used to simulate and predict the formation process and stability of goaf.The classification results show that the quality of the surrounding rock mass of the upper wall and the ore body is good，and the qua-lity of the surrounding rock mass of the lower wall is medium.Through the analysis of stress isoline cloud map，displacement isoline cloud map and plastic zone distribution map，it is concluded that the surrounding rock of goaf No.2 and No.3 in Scheme 1 has a plastic state because the shear stress is too concentrated and higher than that of the rock mass itself，the maximum displacement appears in the roof of No.2 goaf，and the maximum horizontal displacement appears in the center of the top and bottom of the goaf and has a tendency of connection.Although the plastic zone in Scheme 2 is not connected，the existing plastic zone will impose certain potential safety hazard on the stability of the goaf.

Keywords：goaf;numerical simulation;stability;Q system classification method;RMR classification method

收稿日期：2020-12-27; 修回日期：2021-02-05

作者简介：王 伟（1987—），男，山西孝义人，高级工程师，硕士，从事科研管理及采矿技术研究等工作;北京市东城区安定门外大街9号，中国黄金集团有限公司，100011;E-mail：wangwei@chinagoldgroup.com