矿震扰动下立体采空区群动力响应研究

2018-06-14姜立春罗恩民

金属矿山 2018年5期

姜立春罗恩民

（1.华南理工大学土木与交通学院，广东广州510640；2.华南理工大学安全科学与工程研究所，广东广州510640）

地下金属矿山在长期开采过程中遗留大量采空区，随着开采进行将形成采空区群[1-6]。随着开采的深入进行，频繁的采矿开拓作业导致采动应力叠加严重，并导致岩体的承载应力逐渐增加和屈服损伤日益严重，易引发顶板冒落、围岩片帮、冲击地压等地质灾害，当受地质灾害引发的矿震应力波扰动时，极易改变采空区群原始的应力状态，空区群的整体稳定性会受到弱化，严重影响采空区群的安全。因此，开展矿震扰动下采空区群的动态响应规律研究具有非常重要的意义。

目前国内外学者对动荷载作用下地下金属矿山采空区群的动力响应规律开展了研究。唐礼忠等[7]采用应力波方法分析了采空区围岩受到侧崩时的动力响应，采用有限差分软件FLAC3D模拟了某金属矿典型采空区围岩受到动力扰动时的动力特征；贾楠等[8]采用三维动力有限差分法分析了动力扰动下地下金属矿山采空区围岩非线性演化的规律；付建新等[9]利用相似材料模型试验，验证了在外力扰动作用下采空区顶板和矿柱应变的演变规律；W C Zhu等[10]数值模拟了动态扰动波形对地下采空区围岩破坏区发育的动态影响；SK REDDY等[11]采用有限差分的方法在爆破作用下研究了障碍柱的应力分布规律。总体来说，目前金属矿山采空区动力响应问题研究尚处于起步阶段，主要研究方法为数值模拟，动态扰动主要指的是爆破荷载的扰动。然而，针对矿震应力波扰动作用下采空群动力响应的研究尚不多见，亟待开展该领域研究方法，减少巨大工程灾害发生的概率。

以往关于采空区受到矿震作用的研究，主要从单元采空区和矿震机理的角度进行研究[12-13]，关于采空区群在矿震作用下的动态响应规律研究目前较少，而采空区群又在矿山地质中频频出现。为此，笔者在前期工作中借助结构动力学的离散化手段，提出水平和纵向采空区群动力响应模型分析方法，研究单行和单列采空区群的动力响应问题，取得良好成效[13-15]。在此基础上，以某金属矿山8号矿体1 100～1 250 m中段36～42勘探线为工程背景，构建动力响应模型，运用动力学理论，引入剪切力作用系数矩阵B，采用Newmark-β分析方法，利用Matlab工具包和自编程，研究矿震扰动下多自由度系统的动力响应规律，并与数值分析方法结果进行对比，分析该方法的合理性，进而为采空区群矿震防治提供理论指导。

1　工程概况

某金矿位于秦岭褶皱系南秦岭印支褶皱带凤县-镇安褶皱束的北缘，自上而下分为泥盆系中统的王家楞组（D2W）和古道岭组（D2g），矿床赋存于古道岭组（D2g1）下段的含金角砾岩地带，上覆岩层和围岩为王家楞组（D2W）和古道岭组（D2g）的钠长石化板岩、绢云板岩和结晶灰岩。矿体与围岩为过渡接触，矿体边界与角砾岩边界一致时，矿体与角砾岩的边界线容易辨认[14]。矿山由KT5、KT8、KT9等6个矿体组成。

选取8号矿体1 100～1 250 m中段36～42勘探线间的3行、3列的复杂采空区群作为研究对象，如图1所示。单元采空区跨度40 m；间柱宽10 m、高35 m；顶板厚5 m、底板厚10 m；围岩厚度20 m。岩体物理力学参数见表1。

随着采矿开拓作业的深入或残采的进行，矿区内在区域应力场原始应力发生改变处于失调不稳的异常状态，若在空区顶板处积累了一定能量后其极限应力区发生脆性破坏，导致顶板覆岩大面积垮落至空区22顶板上产生了岩层震动，造成了采空区群的矿震灾害，图2为立体采空区矿震示意图。

2　动力响应模型

2.1　基本假设

（1）立体采空区群岩体构造完整，不考虑断层、节理和水的渗流的影响；

（2）立体采空区群顶（底）板、间柱和围岩未发生严重变形。

2.2　模型构建

将采空区群及围岩作为一个系统进行考虑，采空区群系统离散化处理后，考虑到上覆岩层和围岩水平地应力的影响，上覆岩层等效为自重应力。围岩和间柱共3行，每行有4个元素；顶板共2行，每行有3个元素；系统共含18个元素。在平面剖面方向上，1个质点为1个自由度，系统共有18个自由度单元（图3）。

根据结构动力学的D’Alembert原理，采空区系统动力响应的运动微分方程为：

式中，M、C、K、Q分别为质量矩阵、阻尼矩阵、刚度矩阵、剪切力矩阵，分别为加速度矩阵、速度矩阵和位移矩阵，λ为剪切力作用因子矩阵，B为剪切力作用系数矩阵。

3　参数确定

（1）质量矩阵M。根据文献和现场工况，3行、3列立体采空区群各质点的质量矩阵为

（2）刚度矩阵K。根据文献[14-16]可得刚度系数k为：

式中，E为岩体弹性模量，；A为单元采空区围岩的水平截面积，A=1 050 m2；L为单元采空区的高度。

利用求得的刚度系数借助于Matlab自编程可求得刚度矩阵K。

（3）频率系数ω。在Rayleigh阻尼的计算过程中需要考虑结构的振动快慢，即自振频率ω，ω通过求解特征方程式求得：

式中，是振幅矩阵。

（4）阻尼矩阵C。计算式（6）求得模型的18阶频率系数，选取第3阶频率系数78.96 Hz和第7阶频率系数147.95 Hz，计算岩体的Rayleigh阻尼，最终获取系统的阻尼矩阵C。Rayleigh阻尼计算式如下：

式中，和分别是质量阻尼比例系数和刚度阻尼比例系数。

工程中通常取岩体阻尼比，则和为

（5）剪切力矩阵Q。外部荷载扰动下，采空区群内部的顶（底）板、间柱与周边围岩接触面的剪切应力，可通过广义胡克定律进行求解，由文献[14-16]可得剪切力矩阵Q。剪切应力的大小受多种因素影响，受既有工程技术水平发展的限制，精确计算剪切力存在一定的困难。本研究引入剪切力作用因子矩阵进行修正，的取值参考文献[14]。修正后的剪切力矩阵为：

（6）矿震等效冲击荷载矩阵P。若在空区顶板处积累了一定能量后其极限应力区发生脆性破坏，导致顶板覆岩大面积垮落至空区22顶板上产生了冲击荷载，冲击荷载的研究参考文献[14]，本研究中冲击荷载即为矿震等效冲击荷载矩阵P，P的表达式为：

4　类动力响应方程分析方法

对于3行、3列采空区群动力的响应问题，剪切力矩阵和矿震等效冲击荷载矩阵构成复杂，类动力响应运动响应方程计算困难，而且它是一个对时间函数非常敏感的刚性方程，采用显示方法计算时，必须取非常小的时间步长[17]。相比之下，隐式方法对时间步长要求较低，且计算结果更具稳定性和精确性[17]。因此，本文拟采用隐式方法Newmark-β法进行[18]。计算程序如下。

（1）Newmark-β法的稳定性分析表明[17-18]，当β≥0.5，γ≥0.25（0.5+β）2时，它是稳定的（β和γ为算法参数）。设爆破荷载步长Δt=0.001，参数γ=0.5，β=0.25，Newmark-β法中各积分常数为

（4）求解时刻的位移、加速度和速度

选取计算总时间T=0.8 s，循环式（12）至式（16）计算步骤，借助于Matlab自编程对式（3）进行迭代求解，获取爆破荷载作用下各顶（底）板、间柱和围岩的动力响应参数位移矩阵和速度矩阵。

5　数值分析模型

5.1　基本假设

（1）立体采空区群岩体构造完整，不考虑断层、节理和水的渗流的影响

（2）立体采空区群岩体具有均质各向同性的半无限体。

5.2　数值模型构建

根据矿山现场实际情况结合圣维南定理合理简化边界，利用Midas/GTS软件构建采空区3行3列的区域数值模型（如图4）。该模型长为380 m，宽为30 m，左侧高为400 m，右侧高为300 m。

在数值模拟过程中，岩体破坏准则服从摩尔-库仑（Mohr-Coulomb）准则，仅考虑岩体的初始地应力，岩体力学参数参考表1。数值模拟过程中，为避免爆炸振动波在模型边界反射回来影响结果，在模型左右、底部设置粘性边界，前后设置位移条件约束。为了提高数值模拟的精度，对模型进行了网格细分（图4）。

矿震等效荷载施加于如图2在空区22顶板的上向临空面上。

6　动力响应特征分析

利用动力响应模型法计算时，由于动力响应模型是将采空区的顶（底）板和间柱经过离散化处理得到的18个自由度单元质体，每一条质体对应一个动力响应时程曲线，得出动力响应模型有18条动力响应时程曲线。而数值模型法计算时，要经过特征值分析和时程分析获得了动力响应特征参数，但是，数值模型的顶（底）板和间柱是连续分布的，不同节点对应不同的动力响应时程曲线。为了与动力响应模型对比分析，需要对数值模型的顶（底）板和间柱的数值模拟结果进行加权平均处理。

6.1　位移响应特征分析

图5所示动态响应模型法和数值模型法得到的位移响应时程曲线对比分析可以看出它们的变化趋势趋于一致，峰值都处于同一数量级。矿震冲击荷载激励瞬间，采空区的顶（底）板和间柱的位移瞬时到达峰值（见图5（a）～（f）），直接作用岩体的位移是在塑性区域内发生有规律的振动（见图5（c））。相邻中段的顶（底）板和间柱的位移极值随着与矿震冲击荷载的距离增大而减小。由于岩体内部阻尼的能量耗散和剪切力的共同作用，位移的幅值随时间增大而逐渐衰减，最后实现平衡。

图5（c）中，在t=0.0～0.1 s时，岩体振动强烈，振幅较大；t=0.1～0.2 s时岩体位移的振幅忽高忽低，位移时程曲线出现跳跃式的衰减；t≥0.2 s时，位移振幅逐渐衰减，采空区群最终恢复平衡。

6.2　速度响应特征分析

从图6中动力响应模型法与数值模拟法的速度时程曲线对比分析可以看出它们的变化趋势出现相似的规律。首先达到采空区的顶（底）板和间柱的位移迅速到达峰值，由于岩体内部的阻尼和剪切力的共同作用，位移的幅值随时间增大而逐渐衰减，最后实现平衡。直接作用岩体的位移是在塑性区域内发生有规律的振动，相邻中段的顶（底）板和间柱的位移极值随着与矿震冲击荷载的距离增大而减小。

图6中，在t=0.0～0.1 s时岩体动力响应强烈；t=0.1～0.2 s时岩体速度的振幅忽高忽低，速度时程曲线出现跳跃式的衰减；t≥0.3 s速度响应振幅逐渐衰减，采空区群最终恢复平衡。

综合图5和图6分析可知，2种方法存在一定的时间差，主要原因归结于动力响应模型法的能量传递瞬间完成，而数值模拟法的能量是逐步完成的。也可看出直接作用岩体的位移是在塑性区域内发生有规律的振动。同时，动力响应模型法和数值模型法计算得出的顶（底）板和间柱的位移和速度时程曲线变化趋势趋于一致，具有相类似的变化规律。从而在一定程度上揭示了矿震扰动下立体采空区的动力响应规律。