矩形巷道矿震应力波扩散规律谱元法分析*

唐杰兵,鞠文君,焦建康,陈法兵

(1.煤炭科学研究总院 开采研究分院，北京 100013；2.中煤科工开采研究院有限公司，北京 100013；3.天地科技股份有限公司 开采设计事业部，北京 100013)

0 引言

矿震应力波触发巷道冲击是一种常见的现象[1]，研究矿震应力波传播特性的方法较多，由于应力波在介质中传播是一个入射、折射、反射、散射及它们之间相互干扰等的复杂过程，当传播过程中遇到异质性岩体，应力波会发生叠加、产生角波、P波转化为S波以及振幅增大等复杂现象[2]，在理论分析上，采用数学、力学等学科较难建立相应的冲击地压复杂现象的数学模型，较难捕捉到应力波传播过程特征[3-4]。实验室试验也是如此，霍普金森杆试验研究动载应力波的比较多[5]，也有相似材料模拟的动载应力波研究[6-7](震源激发采用爆炸物或摆锤)，其无法真实模拟原岩应力条件和介质的力学环境，或离心试验模型[8](人造重力场)框架周边橡胶垫层只能减轻边界效应，所以这些现场或实验室试验无法解决应力波的边界影响问题，且不便进行参数敏感性分析，现场或实验室试验中只能设置有限的监测点，震源模拟的可靠性缺乏可依靠的研究。在定量化研究中，主要的研究手段还是数值模拟方法，可用于应力波传播的数值模拟软件：FALC2D[9]，UDEC[10]，PFC[11]，ABAQUS[12]等，可以合理地设置应力波的吸收边界条件及数据监测等，但是这些数值模拟及理论分析中常用的是简谐二维波，SPECFEM2D程序包谱元法数值模拟在模拟应力波传播上的应用，弥补了上述几种数值模拟、现场或实验室试验以及理论分析的不足，谱元法的破裂震源更接近于现场实际情况，且SPECFEM2D可用于模拟煤矿巷道等小尺度范围内的矿震应力波扩散[13]。冲击地压巷道支护设计需要同时考虑动静作用，以往的研究侧重于巷道围岩的静载分布，对巷道周围应力波扩散规律研究还不多见。本文采用SPECFEM2D程序包谱元法数值模拟方法，分析应力波在巷道周围扩散过程中的波场快照、质点最大振动速度PPV在巷道围岩表面以及巷道围岩分布规律以及应力波记录等。研究结果对于冲击巷道危险性分析及其合理支护设计具有借鉴意义。

1 谱元法及其计算流程

谱元法(Spectral-element method,SEM)最初是Patera 在1984年流体力学计算中提出[14]，Gharti等[15]在1998年把拉格朗日谱元法用于地震波模拟，推动了应力波的数值模拟研究。谱元法的二维或三维弹性波波动方程强形式可以转化为弱形式，使用基于高斯-勒让德-洛巴托多项式(Gauss-Lobatto-Legendre,GLL)插值点[16]的拉格朗日高阶正交多项式基函数获得对角质量矩阵[17]，大大简化了算法实现和存储量，并在时间域内采用2阶精度的显式Newmark中心差分递推格式求解，有利于计算机并行计算(Parallel computing,MPI)的实现[18]。谱元法能用较少的单元获得与其他方法相同的精度，且具有最优张量乘积技术的逐元算法，可处理复杂地质结构模型。与有限元法的区别在于基函数的选取和质量矩阵中积分计算方法的不同，它兼具了有限元的处理边界和结构的灵活性和谱方法[19]的快速收敛特性。SPECFEM2D谱元法程序包可采用4节点或9节点网格单元对各向同性或各向异性介质进行弹性波场模拟、伴随波场模拟以及全波形反演等，该程序包能模拟声波和弹性波在局部、区域范围内各种介质如流体、弹性体、粘弹性体、各向异性及多孔介质中的扩散过程。另外，SPECFEM2D/3D程序包具有与其他程序进行耦合分析的扩展性。

SPECFEM2D的非结构化网格可由外部程序CUBIT，Gmsh，GID等软件生成，由程序包自带的xmeshfem2D命令也可以进行简单的网格划分，并由SCOTH 进行网格节点的负载分配，并通过转换程序或程序包自带的xmeshfem2D命令生成SPECFEM2D的网格文件、节点文件、材料属性文件、应力波吸收单元文件等程序运行所需的参数输入，然后程序初始化，得到离散方程组中的质量矩阵、刚度矩阵、吸收边界条件矩阵及震源的等效力项，然后用xspcfem2D命令进行迭代求解，计算流程如图1所示。

图1 SPECFEM2D数值模拟流程Fig.1 Procedure of SPECFEM2Dnumerical simulation

2 谱元法数值模型建立

图2 网格建立Fig.2 Grid establishment

图3 SEM网格分块Fig.3 SEM grid database partition

表1 材料属性Table 1 Material properties

3 谱元法数值模拟结果分析

如图4所示，浅色代表应力波振幅为正值，深色代表应力波振幅为负值，方框部分为应力波监测测点，为35 m×35 m的矩形范围，每个测点相隔1 m，巷道内部无测点，无巷道时布置1 332个测点，有巷道时布置1 310个测点，展示了应力波在0.028,0.06,0.076和0.143 s扩散时的波场快照。可以看出，CPML边界应力波吸收层达到了效果。在0.028 s时，震源开始扩散；0.06 s时，P波和S波开始分离；0.076 s时，P波经过巷道，由于巷道与周围围岩体的波阻抗不同，应力波在巷道周围扩散，值得注意的是当应力波不在震源与巷道的直线方向上时，巷道的存在并不会引起应力波扩散的变化，因此，巷道周围的应力波扩散样式与震源类型和它与巷道的相对位置有关；0.143 s时，S波经过巷道，可见在S波经过巷道后，应力波在巷道周围扩散的比在P波经过巷道时剧烈，S波携带了大量的能量，是矿震应力波导致巷道破坏的主要原因。

图4 应力波在有无巷道时扩散过程波场快照Fig.4 Wave field snapshot of stress wave diffusion process with or without roadway

图5 无巷道时巷道周围PPV分布规律Fig.5 Distribution laws of PPV around roadway without roadway

图6 有巷道时巷道周围PPV分布规律Fig.6 Distribution laws of PPV around roadway

如图7所示，在巷道围岩表面布置23个测点，在周围围岩及内部布置28个测点，箭头从1开始依次表示了监测点的序列号，每个测点相隔1 m。如图8所示，巷道围岩表面PPV最大值分布规律为：左帮>顶板>右帮>底板，分别为测点5、测点18、测点15和测点7。较无巷道时，PPV值增大倍数分别为2.4，1.7，1.7和1.4。如图9所示，巷道围岩内部PPV在迎波侧较无巷道时应力波扩散的PPV值大，背波侧除了巷道围岩表面的监测点PPV值大于无巷道时的PPV值，其他均小于无巷道时的PPV值，可见，在背波侧出现了1个应力波扩散的屏蔽区。巷道围岩界面上入射压缩波在自由表面反射成卸载波、再与入射卸载波相互作用，将形成拉应力所产生的卸载效应。在巷道围岩高静载接近煤岩体的极限强度时，极小的应力波扰动引起的静载增量都可能导致冲击地压的发生。

图7 巷道围岩测点布置示意Fig.7 Schematic diagram for layout of monitoring points in surrounding rock of roadway

图8 巷道围岩表面PPV分布规律Fig.8 Distribution law of PPV on surrounding rock surface of roadway

图9 巷道围岩内部PPV分布规律Fig.9 Distribution laws of PPV in surrounding rock of roadway

由图10所示，水平方向X和垂直方向Z的应力波的S波幅值均远大于P波幅值。巷道左帮监测点5表明顶底板X方向应力波振幅放大明显，尾波波动幅度比较大，应力波持时显著增大，Z方向应力波与无巷道存在时相差较小。巷道底板监测点11表明Z方向应力波振幅放大明显，尾波波动幅度比较大，应力波持时显著增大，X方向应力波与无巷道存在时相差较小。受篇幅限制，在同时分析了巷道周围其他监测点应力波幅值大小及衰减快慢后，发现应力波在巷道周围扩散的速度矢量方向有一定的规律，也即巷道顶底板主要受水平方向应力波作用，巷道两帮主要受垂直方向应力波作用，巷道两帮径向约束的最小主应力较小，幅值较大的垂直方向应力波与帮部垂直集中应力叠加，其表现形式为巷道两帮围岩体向开挖空间抛射，出现层裂现象以及巷道肩窝处剪切滑移等，这与某矿现场出现的矿震应力波破坏情况相吻合，如图11所示。

图10 巷道围岩表面测点5和测点11 X及Z方向应力波波形Fig.10 Stress wave waveforms in X and Z directions of No.5 and No.11 measuring points on surrounding rock surface of roadway

图11 矿震应力波致灾现场图示Fig.11 Field diagram of disaster caused by mine earthquake stress wave

4 结论

1)质点最大振动速度PPV>1 m/s在巷道迎波侧的分布范围较大，在背波侧形成了一定范围的应力波屏蔽区。

2)矿震应力波在遇到不同波阻抗的围岩界面时，产生应力波的叠加现象，巷道围岩表面PPV最大值发生左帮上半部分，是无巷道影响时的2.4倍，深入巷道左帮、底板围岩内部的测点有巷道时PPV值反而没有无巷道时大。

3)巷道两帮垂直方向矿震应力波幅值较水平方向大，顶底板水平方向应力波幅值较垂直方向大。且幅值较大的垂直或水平方向应力波尾波幅值较大，持续时间长。

4)PPV值在巷道围岩内形成了升高区和屏蔽区，在巷道围岩集中应力区，较小的应力波扰动极可能诱发冲击显现，数值模拟有助于加深理解巷道冲击地压的矿震应力波诱发机理，并为下一步合理的冲击巷道支护设计提供了参考。