刘晓明 ，罗周全，徐纪成

(1. 中南大学 资源与安全工程学院，湖南 长沙，410083；2. 中国矿业大学 煤炭资源与安全开采国家重点实验室，江苏 徐州，221116；3. 中南大学 现代分析测试中心，湖南 长沙，410083)

隐患空区动力失稳过程是一个极其复杂的力学系统，但从局部而言，都是岩石的变形与破裂问题。随着岩石力学的发展，各种数值计算方法如有限元法、有限差分法、边界元法、离散元法、拉格朗日元法和流形元法等广泛应用于各种岩石力学工程问题中[1−4]。而岩层移动与破裂研究的数值计算方法主要包括 2类：一类是建立在连续介质力学基础上的数值计算方法，如有限元法、边界单元法等；另一类则是非连续介质力学分析方法，如离散元法。近年来，细观力学分析方法为岩石破裂过程的数值模拟方法的发展提供了很好的思路。岩石破裂过程分析系统 RFPA(rock failure process analysis)是基于弹性损伤统计模型的数值模拟工具，能较好地模拟岩石破裂失稳过程，目前已广泛应用于岩石破裂机制以及岩石力学与工程问题的研究，取得了许多成果[5−6]。隐患空区失稳根据失稳类型可以分为静力失稳和动力失稳。静力失稳即在原岩应力场作用下进行空区开挖，由于空区周围应力集中致使空区顶板或两帮出现拉伸破坏或剪切破坏；动力失稳指在静力作用的基础上，受到周围爆破动荷载的作用，致使空区失稳。与静力失稳相比，空区动力失稳模式虽然更复杂，仿真难度更大，但是，更能反映实际空区破坏过程。此外，隐患空区形态边界复杂，如采用矩形、梯形或圆等规则形状进行简化，再用简化的模型进行空区动力失稳仿真，空区的破坏过程与规律往往与实际不符，不能如实反映实际复杂隐患空区动力失稳过程，因此，有必要继承隐患空区三维信息获取的结果，在此基础上进行动力失稳仿真模型构建。华锡集团铜坑矿92号矿体的开采技术条件复杂，井下隐患空区空间分布复杂，而爆破地震波作用下空区发生动力失稳的趋势日益明显。为保障92号矿体的安全高效开采，揭示空区动力失稳的内在机制，有效预防空区塌陷事故的发生，在铜坑矿隐患空区 CMS(cavity monitoring system)探测并获取准确的三维信息的基础上，借助RFPA提供的岩石破裂过程动力仿真分析功能，同时考虑开挖过程和爆破震动作用对空区稳定性的影响，对空区动力失稳过程进行仿真，以便为开展空区灾害预防与控制、空区处理方案的制定、应力监测点布置、资源开采布局等工作提供基础依据。

1 隐患空区三维信息获取

1.1 工程概况

柳州华锡集团铜坑矿主要开采对象为细脉带、91号和92号矿体，这三大矿体在空间位置上重叠产出。目前细脉带、91号等矿体开采即将结束，正在大量开采92号矿体及细脉带及91号矿体边缘薄贫矿体。在长期开采中形成了大量的采空区，局部地段空区的冒落或坍塌已构成威胁矿山生产的重大安全隐患。为能采取有效措施控制因空区引发的地质灾害，急需对现有空区进行彻底探查和分析，以掌握空区空间分布、实际边界、体积等基本情况，并在此基础上对空区群危险性进行分析，确定不稳定或随着周围矿体的开采将不稳定的空区，以便为空区综合治理与控制、资源有效回收提供技术支撑[7−8]。

1.2 隐患空区三维信息获取

空区激光探测系统(CMS)因其能准确测定井下危险采空区及其他硐室空间的三维形态而被广泛应用于国内外金属地下矿山[9−12]。因此，研究采用CMS对铜坑矿92号矿体隐患空区进行探测，输出的点云数据经过滤、拼接及坐标、格式转换，采用三维矿业软件Surpac生成空区三维实体模型，并以此为基础开展模型计算、三维可视化及剖面输出等应用。隐患空区三维信息获取及可视化过程如图1所示[13]。

图1 隐患空区三维信息获取及可视化过程Fig.1 Process of 3D information acquisition of disaster cavity and its visualization

隐患空区三维空间信息主要包括空区位置形态、体积、顶板面积及边界等，这些信息将随时间发生变化，因此，准确、及时掌握隐患空区信息可为矿山进行采矿贫损指标控制、开采设计、空区管理、安全评价及灾害预测等提供可靠的基础性依据。

1.3 计算区域的选取及其特征分析

通过分析铜坑矿 92号矿体探测空区群的三维空间信息，发现92号矿体西南采区内的空区之间互相重叠、贯通，稳定性较差，邻近频繁爆破生产作业诱发采空区失稳破坏的趋势日益明显。为此，结合目前铜坑矿的开采现状，选择92号矿体西南采区内的采空区及其周边矿岩作为本次动力失稳仿真技术研究的计算区域。

因空区三维形态复杂，构建复杂空区三维数值计算模型难度大，为此，选用二维 RFPA，选择典型空区剖面进行分析，相比三维数值分析，二维计算结果安全性更高。因空区主要沿勘探线分布，为选取典型空区进行动力失稳过程模拟仿真，分别沿 205-2，203-2，201-2和202-2这4条勘探线截取剖面，分析剖面上空区动力失稳的过程，其中剖面位置如图2所示。不同剖面上隐患空区的分布特征见表1。

图2 隐患空区及计算剖面位置Fig.2 Disaster cavity and section position

2 隐患空区RFPA数值计算模型

2.1 计算模型

目前，RFPA-dynamic尚未开发与其他三维建模软件或数值软件数据转换接口，只能在软件内建立计算模型，软件提供了完整的用于建立各种直线、曲线、不规则多边形、矩形和圆等工具，方便建立各种边界不规则的隐患空区。因此，充分结合205-2，203-2，201-2和202-2这4条勘探线剖面上空区形态边界，在RFPA内利用不规则多边形完成复杂边界隐患空区数值计算模型的建立，综合考虑模型尺寸及软件计算能力，针对不同模型选取不同的尺寸及单元尺寸，具体参数见表2，最后建立的RFPA计算模型如图3所示。

表1 不同剖面上隐患空区的分布特征Table 1 Distribute characteristic in different sections

表2 RFPA数值计算模型参数Table 2 RFPA numerical calculation model parameters

2.2 基元赋值

RFPA通过 Monte-Carlo方法和统计描述相结合对基元进行初始化赋值[14−15]。设网格中所有基元的弹性模量最大阈值为E0，φ(E)为具有某弹性模量E基元的分布值，则弹性模量Weibull分布函数的积分为

其中：φ(E)为具有弹性模量E的基元的统计量。由上式统计分布构成的基元组成一个样本空间，在阈值E0不变的情况下，由于m的差别，积分空间分布不完全一样。这些基元构成的岩石类介质在宏观上性质可能大体一致(E0相同)，但是，由于细观结构的无序性，使得基元的空间排列方式有显著不同。这种细观上的无序性正好体现了岩石类介质独特的离散性特征。

图3 铜坑矿隐患空区RFPA数值计算模型Fig.3 RFPA numerical calculation models of Tongkeng Mine disaster cavity

一般物理空间随机分布的无序性可以通过Monte-Carlo方法来实现，其产生方法是：产生一组在(0, 1)区间上均匀分布的随机数序列{iγ≤1|i=1, 2, …,n}；对于任何iγ，存在一个与{iγ≤1|i=1, 2, …,n}相对应的随机数序列{Ei|i=1, 2, …,n}。由随机数序列{iγ}映射一组弹性模量参数序列{iγ}。这一组基元弹性模量参数随机序列逐一赋予网络中的每一个基元，其他力学参数(强度、泊松比、容重)同样用此种方法赋值。这种既具有统计性、又具有随机性(无序性)的基元力学参数赋值方法满足了网格中基元力学参数非均匀性和随机性的要求。

根据铜坑矿92号矿体的矿岩的基本条件，采用的基元参数见表3。

表3 岩体物理力学参数Table 3 Physical and mechanical parameters of rock mass

2.3 初始应力与边界条件

静力计算边界条件为：根据铜坑矿的原岩应力场分布，采用上部受压、侧向受原岩应力场作用的平面应变模型，原岩应力场的水平方向应力为25.4 MPa，垂直方向应力根据模型边界以上岩石自重求取。

动力计算边界条件为：在静力计算边界条件的基础上，在空区底板设置爆破地震波。

2.4 爆破震动荷载输入

将现场实测爆破地震波时间−速度数据进行简化，并转换为RFPA能够接受的应力波应力，如图4所示。

图4 爆破震动应力时程曲线Fig.4 Schedule curve of blasting vibration stress

在计算软件中选择静力和动力计算模块，先设置静力加载条件，再设置动力加载条件。其中动力加载可以通过应力、位移、波速和波加速度这4种方式实现。本文选择应力加载方式，在空区底部施加动力荷载。

3 隐患空区动力失稳过程模拟仿真

隐患空区动力失稳仿真的整个计算过程分为静力计算和动力计算。动力分析以静力计算为前提，即在完成静力计算的基础上，再施加动力时程进行计算。静力分析时，首先进行初始应力平衡计算，然后，根据空区形成次序分步模拟空区形成过程，每步形成一个空区并计算至平衡状态。动力计算在静力计算结束后进行，通过在空区底板施加爆破震动波，再进行动力计算，直至空区破坏，最终获得空区从形成到失稳破坏的过程模拟仿真。

3.1 勘探线205-2上的空区

图5所示为勘探线205-2空区动力扰动下失稳破坏过程的剪应力分布。从图5可见：随着T506空区的开挖，空区围岩应力发生重分布，T505空区形成后应力进一步增大，尤其表现在间柱上；施加爆破动载后，空区围岩受力状态有所改变，首先在 T506空区右帮出现裂纹；随着应力波的逐步向上传播，裂纹不断往 T505空区方向扩展。由于空区间柱最小宽度仅为5.5 m，裂纹最先在该位置贯通，致使空区间柱完全破坏。

图5 勘探线205-2剖面空区动力失稳过程仿真Fig.5 Dynamic destabilization numerical simulations of cavity in section 205-2

3.2 勘探线201-2上的空区

图6所示为勘探线201-2空区动力扰动下失稳破坏过程的剪应力分布。从图6可见：

(1) 由于S1空区与其他空区间距高达82 m，因此，其他空区的开挖对其影响较小。动荷载作用后，空区左右两帮周围围岩出现局部裂纹，但并未扩展。

(2) T111与T112空间上呈现上下斜交，间柱最小宽度仅为6 m，静力作用下随着空区的逐步开挖，空区围岩应力发生重分布，间柱出现应力集中。

(3) 动荷载作用后，T112空区右帮首先出现局部片帮跨落现象，随着爆破应力波的进一步作用，片帮范围逐步加大，并有往下扩展的趋势。

(4) 动荷载进一步作用后，T110-111顶板靠T112处出现裂纹，扩展与 T112片帮区域连通，最后，整个间柱出现断裂破坏。

图6 勘探线201-2剖面空区动力失稳过程仿真Fig.6 Dynamic destabilization numerical simulation of cavity in section 201-2

3.3 勘探线202-2上的空区

图7所示为勘探线202-2空区动力扰动下失稳破坏过程的剪应力分布。从图7可见：

(1) 随着空区的逐步开挖，空区围岩应力发生重分布，尤其在空区开挖第3步后，由于空区分布密集，空区间柱跨度小，使得空区间柱应力集中。

(2) 空区开挖第 4步后，空区间柱应力集中进一步加强，并出现局部破坏，如T214与T213，T215与T216的间柱均不同程度出现裂纹，此时，裂纹并未扩散。

(3) 动荷载作用后，应力集中程度得到加强，间柱裂纹出现扩展，并沿着间柱薄弱地带发展，致使T214与T213间柱首先发生破坏。

(4) 动荷载进一步作用后，T215与T216以及T211与 T212的间柱依次出现破坏，致使整个空区向下垮塌。

图7 勘探线202-2剖面空区动力失稳过程仿真Fig.7 Dynamic destabilization numerical simulation of cavity in section 202-2

3.4 勘探线203-2上的空区

图8所示为勘探线203-2空区动力扰动下失稳破坏过程的剪应力分布。从图8可见：

图8 勘探线203-2剖面空区动力失稳过程仿真Fig.8 Dynamic destabilization numerical simulation of cavity in section 203-2

(1) 模型中T214-316与T312呈上下叠加式，即T314，T315和T316空区的顶板为T312空区的底板，相交处为最危险处，容易出现断裂。

(2) 空区开挖第 3步后，空区间柱出现明显的应力集中现象，T314，T315和T316空区顶板出现局部破坏。

(3) 动荷载作用后，空区围岩应力出现二次分布。T312空区靠T314底板出现裂纹，随着动荷载的作用，T312空区右帮出现局部破坏。

(4) 动荷载进一步作用后，原有裂纹继续扩展，最后空区间柱全部垮塌破坏。

4 结论

(1) 在研究 RFPA基本思想和岩石破裂过程仿真原理的基础上，提出了运用RFPA开展隐患空区动力失稳模拟仿真的方法。

(2) 在对铜坑矿隐患空区现场实测并准确获取相关信息的基础上，结合隐患空区危险性可视化分级结果，对处于205-2，203-2，201-2和202-2这4条勘探线上的空区空间分布特征进行了分析，为空区动力失稳过程模拟仿真提供便利。

(3) 模拟了空区从形成到动力扰动作用下的破坏过程，获得了空区破坏过程的裂纹扩展规律，形成了金属矿隐患空区动力失稳过程模拟方法。

[1] 何满潮, 李学元, 刘斌, 等. 深部开采岩体力学研究[J]. 岩石力学与工程学报, 2005, 24(16): 2803−2813.HE Manchao, LI Xueyuan, LIU Bin, et al. Study on rock mechanics in deep mining engineering[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(16): 2803−2813.

[2] 林杭, 曹平, 李江腾, 等. 基于Surpac的Flac3D三维模型自动构建[J]. 中国矿业大学学报, 2008, 37(3): 339−342.LIN Hang, CAO Ping, LI Jiangteng, et al. Automatic generation of Flac3Dmodel based on Surpac[J]. Journal of China University of Mining & Technology, 2008, 37(3): 339−342.

[3] 王涛, 陈晓玲, 杨建. 基于3DGIS和3DEC的地下洞室围岩稳定性研究[J]. 岩石力学与工程学报, 2005, 24(19): 3476−3481.WANG Tao, CHEN Xiaoling, YANG Jian. Study on stability of underground cavern based on 3DGIS and 3DEC[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(19):3476−3481.

[4] Tang C A, Tang S B. Applications of rock failure process analysis (RFPA) method[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2011, 3(4): 352−372.

[5] 朱万成, 唐春安, 黄志平, 等. 静态和动态载荷作用下岩石劈裂破坏模式的数值模拟[J]. 岩石力学与工程学报, 2005, 24(1):1−7.ZHU Wancheng, TANG Chunan, HUANG Zhiping, et al.Numerical simulation on splitting failure mode of rock under static and dynamic loadings[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(1): 1−7.

[6] ZHU Wancheng, TANG Chunan. Numerical simulation on shear fracture process of concrete using mesoscopic mechanical model[J]. Constructions and Building Materials, 2002, 16(8):453−463.

[7] 高峰. 顶板诱导崩落机理及次生灾变链式效应控制研究[D].长沙: 中南大学资源与安全工程学院, 2009: 20−24.GAO Feng. Study on mechanism of induced caving roof and control of its secondary catastrophe chain effect[D]. Changsha:Central South University. School of Resources and Safety Engineering, 2009: 20−24.

[8] 苏家红. 大范围隐患区资源开采安全控制技术研究[D]. 长沙:中南大学资源与安全工程学院, 2008: 1−2.SU Jiahong. Study on safety control to mining under large scale hidden trouble zone[D]. Changsha: Central South University.School of Resources and Safety Engineering, 2008: 1−2.

[9] LUO Zhouquan, LIU Xiaoming, ZHANG Bao, et al. Cavity 3D modeling and correlative techniques based on cavity monitoring[J]. Journal of Central South University of Technology, 2008, 15(5): 639−644.

[10] LIU Xiling, LI Xibing, LI Faben, et al. 3D cavity detection technique and its application based on cavity auto scanning laser system[J]. Journal of Central South University of Technology,2008, 15(4): 285−288.

[11] LIU Xiaoming, LUO Zhouquan, YUAN Wenni, et al. CMS applications in metal mines of China[C]// Proceedings of Pacific-Asia Workshop on Computational Intelligence and Industrial Application. Wuhan, China, 2009: 353−356.

[12] 刘晓明, 罗周全, 孟稳权, 等. 深井采场大规模垮塌三维探测及可视化计算[J]. 中南大学学报: 自然科学版, 2011, 42(1):158−163.LIU Xiaoming, LUO Zhouquan, MENG Wenquan, et al. 3D monitoring and visible calculation of large-scale collapse area in deep mine[J]. Journal of Central South University: Science and Technology, 2011, 42(1): 158−163.

[13] 刘晓明, 罗周全, 袁雯妮, 等. 基于 CMS的隐患空区三维特征信息获取[J]. 科技导报, 2011, 29(5): 32−36.LIU Xiaoming, LUO Zhouquan, YUAN Wenni, et al. 3D information acquisition of disaster cavity based on CMS[J].Science & Technology Review, 2011, 29(5): 32−36.

[14] 谢林茂, 朱万成, 王述红, 等. 含孔洞岩石试样三维破裂过程的并行计算分析[J]. 岩土工程学报, 2011, 33(9): 1447−1455.XIE Linmao, ZHU Wancheng, WANG Shuhong, et al.Three-dimensional parallel computing on failure process of rock specimen with a pre-existing circular opening[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(9): 1447−1455.

[15] 黄志平, 唐春安, 马天辉, 等. 卸载岩爆过程数值试验研究[J].岩石力学与工程学报, 2011, 30 (S1): 3120−3128.HUANG Zhiping, TANG Chunan, MA Tianhui, et al. Numerical test investigation on unloading rockburst processes[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(S1):3120−3128.