周期冲击载荷下巷道顶板开裂机理数值模拟

王学滨1,刘桐辛2,田 锋2,钱帅帅2

(1.辽宁工程技术大学 计算力学研究所,辽宁 阜新 123000; 2.辽宁工程技术大学 力学与工程学院,辽宁 阜新 123000)

摘 要:周期冲击载荷在采矿、土木等工程中较为常见。在周期冲击荷载作用下,巷道围岩会在极短时间内产生大量裂缝,块体会以猛烈方式脱离巷道围岩,造成岩爆、垮塌等灾害。自主开发的拉格朗日元和离散元耦合的连续-非连续方法主要包括4个计算模块:应力应变模块、节点分离模块、接触力求解模块和运动方程求解模块。鉴于单元畸变和局部自适应阻尼可能导致单元弹射出模型时速度失真,为了准确模拟弹射现象,在该连续-非连续方法的基础上,对弹射单元进行了刚化处理。该方法的实质是确保弹射单元在碰撞前后作刚体运动,这是通过对节点速度取平均、消除单元应力实现的。对周期冲击载荷作用下巷道围岩的变形—开裂—运动过程进行了模拟。模型尺寸为40 m×40 m,被剖分成160×160个正方形单元,巷道尺寸为6 m×6 m。模型的左、右侧面均为透射边界,即应力波经过时不会发生反射。通过在巷道顶板布置多个测点,获取了最大主应力在整个模拟过程中的演化规律。该过程包括巷道开挖前模型的平衡过程、巷道开挖后模型的平衡过程及应力波在模型中传播过程。研究发现,当应力波在模型中传播时,弹射单元任一节点的水平速度呈现上升—稳定—衰减的变化过程,这与不进行刚化处理的结果相比更符合实际;巷道顶板左、右对称线上大部分测点的最大主应力呈现近似正弦波动上升—衰减—稳定的变化过程。在某一节点分离前,最大主应力首先呈现有规律的波动;然后,波动幅度突然增大,并伴随着剧烈震荡。阐明了顶板的拉裂机理:当应力波传入模型之后,顶板某一位置的最大主应力的波峰和波谷随着时间的增加呈增加的趋势,反射的多个拉应力波的累积作用不断提升该位置最大主应力的波峰,致使拉裂。另外,还初步分析了周期冲击载荷频率对单元弹射初速度和剪、拉裂缝区段数目的影响。分析发现,随着频率的增加,脱离巷道围岩的单元的平均弹射初速度增加,剪、拉裂缝区段数目减小。

关键词:周期冲击载荷;巷道围岩;拉裂;弹射初速度;刚化处理;频率

中图分类号:TD32;O347

文献标志码:A

文章编号:0253-9993(2021)10-3106-10

移动阅读

收稿日期:2020-10-25

修回日期:2020-12-03

责任编辑:郭晓炜

DOI:10.13225/j.cnki.jccs.2020.1674

基金项目:国家自然科学基金面上资助项目(51874162)

作者简介:王学滨(1975—),男,黑龙江双鸭山人,教授,博士生导师。E-mail: wxbbb@263.net

引用格式:王学滨,刘桐辛,田锋,等.周期冲击载荷下巷道顶板开裂机理数值模拟[J].煤炭学报,2021,46(10):3106-3115.

WANG Xuebin,LIU Tongxin,TIAN Feng,et al.Numerical simulation of tunnel roof cracking mechanism under periodic impact loads[J].Journal of China Coal Society,2021,46(10):3106-3115.

Numerical simulation of tunnel roof cracking mechanism under periodic impact loads

WANG Xuebin1,LIU Tongxin2,TIAN Feng2,QIAN Shuaishuai2

(1.Institute of Computational MechanicsLiaoning Technical UniversityFuxin 123000,China; 2.College of Mechanics and EngineeringLiaoning Technical UniversityFuxin 123000,China)

Abstract:Periodic impact loads are common in mining and civil engineering. The roadway surrounding rock will produce a large number of cracks in a very short time under periodic impact loads, and rock blocks will be ejected from the roadway surrounding rock in a violent way, causing rockbursts and collapses, etc. The continu-um-discontinuum method has been developed, which the Lagrangian element method and the discrete element method are combined. The method includes four modules: the stress-strain module, the nodal separation module, the solving contact force module and the solving motion equation module. The element distortion and local adaptive damping may lead to inaccurate velocities for ejected elements from a model. To accurately model the phenomenon, the ejected elements are treated as rigid bodies based on the continuum-discontinuum method. The present treatment can ensure the rigid body movement of the ejected element without contact, which is achieved by averaging the velocities of nodes and eliminating the stresses of elements. The deformation-cracking-movement process of the roadway surrounding rock under periodic impact loads is modeled. The size of the model is 40 m×40 m. The model is divided into 160×160 square elements, and the size of the roadway is 6 m×6 m. The left and right sides of the model are transmitting boundaries so that no reflection occurs when the stress wave passes through. A few monitored points are arranged at the roadway roof, and the evolution of their maximum principal stresses in the whole calculation process is presented. The process includes the balanced process of the model before excavation, the balanced process of the model after excavation and the propagation process of stress wave in the model. The following results are obtained. When the stress wave propagates in the model, the horizontal velocities of nodes of ejected elements exhibit an increase, followed by a constant and then a decay. This is more practical than the result that the ejected elements are not treated as rigid bodies, and the maximum principal stresses of most monitored points at the symmetric line between the left and right parts of the roadway roof exhibit an approximately fluctuating increase, followed by a decay until a constant is reached. Before a node separate, the maximum principal stress regularly fluctuates, followed by a sudden and violent fluctuation with a large amplitude. The mechanism of roof cracking is expounded. After the stress wave is introduced into the model, the peaks and troughs of maximum principal stresses at the positions of roof increase with time, and the accumulation of several reflected tensile stress waves continuously elevates the peaks of maximum principal stresses at positions, resulting in cracking. In addition, the effects of frequency on the initial ejected velocity and the number of shear and tensile crack segments are preliminarily analyzed. With an increase of frequency, the averaged initial ejected velocity of each element ejected from the roadway surrounding rock increases, and the number of shear and tensile crack segments decreases.

Key words:periodic impact load;tunnel surrounding rock;tensile cracking;initial ejected velocity;rigid treatment;frequency

采矿、土木等工程中岩层周期破断、断层周期黏滑和反复爆破都能产生周期冲击载荷。周期冲击载荷是诱发岩爆、垮塌和片帮等灾害的重要因素。在周期冲击载荷作用下,巷道顶板和两帮会在极短时间内产生大量相互联通的裂缝,进而将巷道围岩切割成众多块体,它们会以猛烈方式脱离巷道围岩,造成灾害[1-4]

轻微岩爆发生时,可以听见清脆的噼啪、撕裂声,偶有爆裂声,岩块弹射初速度小于1 m/s;而强烈岩爆发生时,可以听见似炸药爆破的爆裂声,声响强烈,岩块弹射初速度处于5~10 m/s,极易造成重大危害。目前,关于岩块弹射初速度已有不少来自实验室和现场的数据积累[5-6]。例如,王之东等[5]在单轴压缩条件下对带方孔的3种岩石试样的能量释放过程进行了观测,其中岩块弹射初速度处于1.6~21.0 m/s;在加拿大和南非的岩爆现场[6],岩块弹射初速度处于7.65~12.60 m/s。对岩块弹射初速度估算的方法主要包括理论方法和数值模拟研究方法。例如,在前者方面,秦剑锋等[7]基于岩板的屈曲失稳理论估算的岩块弹射初速度在9 m/s;在后者方面,陈滔等[6]基于能量守恒原理和破坏前后模型的应变能差估算的岩块弹射初速度处于8.16~13.60 m/s。客观地讲,岩爆过程较为复杂,受多种因素影响。在理论上想准确估算岩块弹射初速度极为困难。而且,在基于连续方法(例如,有限元法)的数值模拟中,仅能通过应变能来估算岩块弹射初速度,而不能模拟出岩块从围岩中脱离的过程,进而难以准确估算岩块弹射初速度。

对于巷道的拉裂机理,目前已取得了一些重要进展。例如,李夕兵等[8]采用PFC2D对动载作用下深部巷道围岩的动力响应进行了数值模拟,发现巷道顶板裂纹数量增加是由静态的拉应力与应力波到达顶板后反射的拉应力波叠加所引起;ZHU等[9-10]采用Autodyn2D对爆破诱发的巷道围岩破坏过程进行了数值模拟,发现从自由边界反射的拉应力波引起了距离自由边界一定距离的环向裂纹。

周期冲击载荷3要素是幅值、频率和持续时间。周期冲击载荷3要素的影响研究一直受到重视[11]。前人已对周期冲击载荷频率的影响开展了一定的研究。载荷频率的取值范围在不同文献中不尽相同[12-13]。例如,闫长斌等[12]采用FLAC在动载作用下对地下巷道群的频率影响进行了数值模拟,频率处于5~100 Hz;陈国祥等[13]采用FLAC对半正弦波作用下巷道围岩的破坏过程进行了数值模拟,频率处于5~50 Hz。现有的研究表明,载荷频率对岩样及岩石结构的动力响应都有一定的影响[13-18]。例如,陈国祥等[13]发现随着频率的减小,巷道两帮围岩的最大垂直和水平应力增加;宫凤强等[16]在常规静载和“预静载+扰动”条件下对岩石断裂特性的频率影响进行了实验,发现随着频率的增加,断裂韧度呈线性减小的趋势;SU等[17-18]在真三轴加载条件下对含孔洞岩样径向应力梯度的频率影响进行了实验,发现随着频率的增加,岩爆更易发生。目前,该方面研究主要集中在实验方面和基于连续方法的数值模拟方面,而基于连续-非连续方法的研究还十分少见。

鉴于单元畸变和局部自适应阻尼可能导致单元弹射出模型时速度失真,为了准确模拟弹射现象,在自主开发的拉格朗日元与离散元耦合的连续-非连续方法[19]的基础上,对弹射单元进行了刚化处理,研究了周期冲击载荷作用下巷道围岩的变形-开裂-运动过程,还初步分析了周期载荷频率对单元弹射初速度和剪、拉裂缝区段数目的影响。

1 连续-非连续方法简介

1.1 原始方法简介

拉格朗日元和离散元耦合的连续-非连续方法主要包括4个计算模块[19]:应力应变模块、节点分离模块、接触力求解模块和运动方程求解模块。

在应力应变模块中,采用了混合离散方法,通过节点的速度利用高斯定理求解单元的应力和应变,可以避免沙漏问题。

在节点分离模块中,通过引入虚拟裂缝模型,处理应变软化问题。分别选取最大拉应力准则和莫尔-库仑准则作为拉裂和剪裂判据。当节点分离后,虚拟裂缝产生,虚拟裂缝面之间存在法向及切向黏聚力。通过同时引入Ⅰ型和Ⅱ型断裂能计算法向及切向黏聚力。当法向或切向黏聚力降至0时,虚拟裂缝成为真实裂缝。

在接触力求解模块中,采用了基于空间划分的接触检测方法和基于势的接触力求解方法,具有接触检测效率较高、无需对“角-角”接触问题进行特殊处理的优势。

在运动方程求解模块中,采用中心差分方法求解节点的速度和位移,具有计算效率高、计算精度较高的优势。

1.2 脱离模型的单元作刚体平动的必要性及处理方法

网格法(有限元法、有限差分法等)容易出现网格或单元畸变问题。本文的连续-非连续方法属于网格法。当某单元弹射出模型后,该单元的4个节点的速度一般并不相同,而且,各节点的运动是相互独立的,所以,该单元的变形、运动规律将极为复杂,从而可能导致该单元发生畸变,这会使应力、应变和由应力引起的弹性力等计算错误,从而会进一步加剧该单元畸变。当某单元弹射出模型后,在与其他单元碰撞之前和之后,该单元应仅受重力作用,而在水平方向不受力,这样,该单元的各节点水平速度vx应是常量。一旦该单元发生畸变,将不能保证这一点。另外,本文方法使用局部自适应阻尼(与常见的黏性阻尼不同),由于其方向取决于节点速度的方向,而大小取决于不平衡力(由应力引起的弹性力和重力等构成)的大小,这也可能导致vx不是常量。

为此,笔者提出一种处理方法以确保弹射单元各节点速度的计算结果正常,即在与其他单元碰撞之前和之后,保持恒定。

对于弹射单元,首先,需找到应力状态较为接近于自由状态的临界时步数目。考虑到该单元刚被弹射出时常处于压缩状态,因而其刚进入完全拉伸状态时的应力状态较为接近于自由状态。为此,将该单元最小主应力σ1>0(在本文中,σ3σ1σ3为最大主应力)对应的时步数目作为临界时步数目。然后,将该单元4个节点的速度矢量取平均,获得平均速度矢量,并赋给这些节点,并对该单元的应力清零,即清除弹性力。各节点的新速度v0取为

(1)

式中,vi(i=1~4)为该单元各节点的原速度。

上述处理迫使弹射单元在碰撞之前和之后作刚体平动。由于弹射单元的应力已被清除,单元畸变不会进一步发展。而且,由于各节点只受重力(或为不平衡力的唯一成分)和局部自适应阻尼力(仅位于垂直方向上,其大小取决于重力的大小,其方向取决于节点垂直方向速度vy),各节点在水平方向上将作匀速直线运动。这样,即可避免节点速度的计算过于失真,从而确保岩块弹射现象的模拟结果尽量真实。

2 计算模型及方案

本文的计算模型(以下简称模型)建立依据某大巷[20]所处地质条件,岩性中硬。

模型尺寸为40 m×40 m,被剖分成160×160个正方形单元。压应力波施加后的力学模型如图1(a)所示。应当指出,模型的左、右侧面均为透射边界,即应力波经过时不会发生反射。

参数取值:面密度ρ=2 700 kg/m2,弹性模量E=20 GPa,泊松比μ=0.3,抗拉强度σt=5 MPa,法向接触刚度Kn=10 GPa,莫尔-库仑准则中的黏聚力c=20 MPa、内摩擦角φ=40°,摩擦因数μ0=0.1(用于摩擦力的计算[19]),Ⅰ型断裂能型断裂能 000 N/m,局部自适应阻尼系数α=0.2,时步长度Δt=1.979 24×10-5 s,Δt小于临界时步长度(7.919 6×10-5 s),以确保数值稳定性。

在巷道顶板,布置了5个测点,在巷道左帮,布置了6个监测单元,具体如图1(b)所示。计算在平面应变、大变形条件下进行。

图1 力学模型、测点和监测单元
Fig.1 Mechanical model,monitored points and monitored elements

陈建君等[21]将冲击载荷简化为半正弦波,采用的半正弦压应力波P(N)表达式为

(2)

式中,Pmax为压应力波幅值,取7.3 MPa;频率f=ω/π;ω为圆频率;N为时步数目。

计算可分为3个过程:

(1)对开挖前的模型进行计算,直到模型处于静力平衡状态(所用N=12 000)。

(2)开挖尺寸为6 m×6 m的巷道(所用N=4 000);此后,继续计算,直到模型处于静力平衡状态(所用N=4 000)。

(3)当N=20 000时,在模型的上端面施加竖直向下的周期冲击载荷(半正弦压应力波),此后,继续计算;当N=64 000时最后1个压应力波波尾传入模型,此后一段时间之内计算仍在继续。在最后1个压应力波在围岩中逐渐消失的过程中,由于围岩应力场的略微调整,一些裂缝仍有可能进一步发展,而且,脱离围岩的单元仍在运动。因此,多计算一段时间是必要的。

共采用4个计算方案。方案1~4的f分别为15,25,35及45 Hz。文献[22]通过现场观测发现,震动波的频率主要集中于50 Hz以内。本文选取的f涵盖了上述范围。

3 计算结果及分析

3.1 周期冲击载荷作用下巷道围岩变形-开裂-运动过程

3.1.1 多个压应力波冲击下剪、拉裂缝的时空分布

以方案1为例简单分析多个压应力波冲击下剪、拉裂缝的时空分布。

图2为方案1的剪裂缝与最大主应力σ3的时空分布规律,黑色线段代表剪裂缝区段。图3为方案1的拉裂缝与最大主应力σ3的时空分布规律,黑色线段代表拉裂缝区段。由图2,3可以发现,首先,第1个压应力波传至巷道两帮后,巷道两帮产生剪裂缝,并逐渐发展形成V形坑,其内产生少量拉裂缝(图2(a),(b)和图3(a),(b))。然后,后续压应力波陆续传入模型,巷道两帮既有V形坑外若干新的V形坑形成,其内仍产生少量拉裂缝,与此同时,巷道顶板产生拉裂缝,并发展形成层裂结构,巷道两帮脱离围岩的一些单元弹入巷道(图2(b),(c)和图3(b),(c))。文献[1]采用真三轴试验机对含预制矩形巷道的立方体岩样进行压缩实验,发现巷道两帮出现岩块弹射现象,岩块明显堆积于巷道底板;文献[23]利用水泥类膨胀胶凝材料与水反应体积骤增的特点对巷道围岩进行单次冲击实验,发现巷道顶板和两帮均发生破坏,且两帮存在V形坑。上述现象的条件尽管与本文有所不同,但上述现象与本文的结果(图2(b),(c)和图3(b),(c))基本一致。最终,最后1个压应力波传出模型后,即在围岩中消失后,剪、拉裂缝停止发展。本文方法在处理开裂、接触和摩擦等非线性问题时不可避免存在微小误差,随着计算的进行,这种误差可能被放大,这会导致剪、拉裂缝的分布不具有严格的对称性(例如,图2(c),(d)和图3(c),(d))。这种现象是非线性数值模拟方法的共性。应当指出,2个单元之间的裂缝称之为1个裂缝区段,裂缝区段的形状为四边形。若干裂缝区段连在一起构成裂缝。考虑到单元脱离围岩后裂缝将变得很大,图2~3仅显示了各边长度均不大于1个单元边长的裂缝区段。

图2 方案1的剪裂缝和σ3的时空分布规律
Fig.2 Spatiotemporal distributions of shear cracks and σ3 of Scheme 1

图3 方案1的拉裂缝和σ3的时空分布规律
Fig.3 Spatiotemporal distributions of tensile cracks and σ3 of Scheme 1

3.1.2 巷道左帮监测单元的右下角节点vx演化

图4为方案1的1~6号监测单元的右下角节点vx-N曲线。由图4可以发现,在压应力波传入模型之后,随着N的增加,2~6号监测单元的右下角节点vx(对于任一监测单元,弹入巷道后其4个节点vx均相同)呈现上升—稳定—衰减的变化过程,而1号监测单元的右下角节点vx呈现波动—稳定的变化过程。应当指出,当N=20 000时第1个压应力波开始传入模型,当N=64 000时最后1个压应力波波尾传入模型。

图4 方案1的1~6号监测单元的右下角节点vx-N曲线
Fig.4 Evolution of horizontal velocities of lower right corner nodes of monitored elements 1-6 with time steps of Scheme 1

vx有上升趋势的阶段,2~6号监测单元的右下角节点vx分别上升至最大值,此阶段位于压应力波传入模型之后,vx上升的过程是监测单元逐渐脱离围岩的过程。例如,当N=28 120~31 990时,2号监测单元的右下角节点vx由0.396 m/s增至峰值6.8 m/s,此时,2号监测单元脱离巷道围岩(图2(b))。期间,vx的最大值就是弹射初速度,分别为6.8,5.8,9.7,9.5和10.3 m/s,上述弹射初速度的平均值,即平均弹射初速度,为8.43 m/s。这位于强烈岩爆[24]发生时平均岩块弹射初速度范围(5.0~10.0 m/s)之内。

vx稳定的阶段,2~6号监测单元的右下角节点vx保持不变。期间,这些监测单元未与其他单元发生碰撞,即在水平方向上不受力,因而没有产生能量损耗。由于巷道的空间有限,这些监测单元一定会与其他单元发生碰撞,因此这些监测单元的右下角节点vx无法总保持不变。这些监测单元的位置和弹射初速度的不同导致vx发生变化的时刻不同,有的位于压应力波传出模型之前,而有的位于压应力波传出模型之后。

vx有衰减趋势的阶段,2~6号监测单元的右下角节点vx呈现总体衰减的趋势。此阶段存在vx反向现象。首先,当N=74 661~74 760时,2号监测单元的右下角节点vx由6.8 m/s降至-5.16 m/s。这是因为该监测单元与其他单元发生碰撞导致反向并产生能量损耗(图2(c));然后,该监测单元不断与其他单元发生碰撞,导致其右下角节点vx发生多次反向现象,并呈衰减趋势(图2(d))。

图5为刚化处理前后2号监测单元的右下角节点vx-N曲线。为了呈现对弹射单元进行刚化处理前后的差异,对方案1重新进行了计算。由此可以发现,刚化处理前2号监测单元的右下角节点vx杂乱无章,这并不是该单元与其他单元碰撞的结果,而刚化处理后2号监测单元在未与其他单元碰撞时,可以保持恒定的vx,这较为符合该单元在水平方向上不受力的事实,这在一定程度上说明了刚化处理的正确性。

图5 刚化处理前后2号监测单元的 右下角节点vx-N曲线
Fig.5 Evolution of horizontal velocities of lower right corner nodes of monitored elements 2 with time steps before and after rigid treatments

3.1.3 巷道顶板测点σ3的演化规律及开裂机理

以方案1为例,简单分析巷道顶板各测点的σ3N的演化规律。图6为方案1的1~5号测点的σ3-N曲线。由图6可以发现,在压应力波传入模型(N=20 000)之前,首先,各测点的σ3稳定在-27 MPa,这对应于巷道开挖之前模型的静力平衡状态;随后,受巷道开挖的影响,各测点的σ3先经历震荡上升、后逐渐衰减至稳定的变化过程。在压应力波传入模型之后,大部分测点的σ3呈现近似正弦波动上升—衰减—稳定的变化过程。

图6 方案1的1~5号测点的σ3-N曲线
Fig.6 Evolution of σ3 of monitored points 1-5 with time steps of Scheme 1

下面以4号测点为例,详细分析近似正弦波动上升阶段σ3的演化规律。在此阶段中,在σ3-N曲线上,可隐约观察到13次波动,这是因为方案1中施加了13个压应力波,各压应力波均会对该测点的σ3产生影响。图7为方案1的4号测点的σ3-N曲线。由图7可以发现:

图7 方案1的4号测点的σ3-N曲线
Fig.7 Evolution of σ3 of the monitored point 4 with time steps of Scheme 1

(1)σ3首先呈现有规律的波动,然后波动幅度突然增大,并伴随着剧烈震荡。σ3-N曲线的波峰和波谷均有随着N的增加而增加的趋势。

N=20 460~23 760时,σ3出现第1次波动。这表明,第1个压应力波传至并逐渐传过4号测点。

N=20 460~22 130时,σ3由-15.8 MPa减小至-30.94 MPa。这表明第1个压应力波的峰前部分逐渐传至该测点。应当指出,当N=20 800~20 970时,σ3有小幅度震荡。随后,随着N的增加,σ3继续下降,σ3-N曲线斜率的绝对值明显小于N=20 800之前的。也就是说,σ3小幅度震荡前后σ3-N曲线的斜率有所不同,前面曲线更陡,而后面更平缓。σ3小幅震荡的原因是第1个压应力波的前缘经由巷道顶板表面反射的拉应力波传至该测点。此后,第1个压应力波的后续部分继续传过该测点,并与反射的拉应力波发生叠加,致使σ3-N曲线的斜率改变。

N=22 131~22 480时,σ3-N曲线出现第1个波谷,这是因为第1个压应力波的峰值部分刚传至该测点,与此同时,该测点仍位于第1个压应力波反射的拉应力波中。

N=22 481~23 760时,第1个压应力波的峰后部分逐渐传至该测点,σ3由-30.04 MPa增大至-11 MPa。当N=23 760时,第1个压应力波已完全传过该测点,与此同时,上述拉应力波尚未完全传出该测点,σ3-N曲线出现第1个波峰。第2个压应力波立即传至该测点,导致σ3下降。应当指出,与4号测点刚受第1个压应力波影响时(N=20 460)的σ3相比,第1个压应力波刚传出该测点时(N=23 760)的σ3更大,这是因为第1个压应力波反射的拉应力波的作用。

N=20 460~46 889时,σ3规律性波动,期间,波峰和波谷随着N的增加均有增加趋势。波峰变化的原因同前所述。波谷变化的原因为:反射的多个拉应力波会对传至该测点的压应力波有一定的衰减作用,致使压应力波的作用不再强烈。

N=46 890~47 040时,与此前几次有规律的波动相比,σ3的波动幅度突然增大,震荡加剧,这与该测点下方节点的分离(介质的拉裂)有关。该测点下方节点的分离将使σ3向围岩的深部转移,从而将提升该节点的σ3。此外,该测点下方节点的分离将对应力波的已有传播产生影响,因而该测点σ3的震荡将加剧。

(2)节点发生分离后,σ3的震荡幅度有衰减趋势。

N=53 190时,4号测点的σ3达到σt(5 MPa),顶板在此处将发生拉裂。下面,阐明顶板的拉裂机理。若σ3的波峰不发生改变,则顶板不可能被拉裂。σ3的波峰随着N的增加而增加使顶板拉裂成为可能。当σ3极小时,波峰的σ3为负,代表压缩;当σ3极大时,波峰的σ3为正,代表拉伸。前一个压应力波刚完全传过该测点之时,恰是下一个压应力波刚传至该测点之时。下一个压应力波刚传入该测点时的σ3总比前一个压应力波的大,这说明该测点的最大σ3伴随着应力波的不断传入和传出而越来越高,在不断累积,这种累积只能源于反射的多个拉应力波的作用,即反射的拉应力波每通过一次该测点将其σ3提升一次,直至达到σt

此后,σ3的震荡规律变得复杂。与此同时,震荡幅度有衰减的趋势,直至稳定在零值附近。水平展布的裂缝会阻断应力波的传播,并引起应力波的反射,从而引起该测点σ3的剧烈震荡。经由巷道顶板表面反射的拉应力波产生的拉伸应力应在接近垂直方向上。所以,顶板将产生水平展布的拉伸裂缝。与此同时,若干传入模型的压应力波相当于增加了模型所受的垂直方向应力。由此,顶板将产生垂直方向展布的拉伸裂缝。

此外,由图6还可以发现,从整体上看,在近似正弦波动上升阶段,随着N的增加,大部分测点的σ3波动上升。越靠近巷道顶板表面,σ3的波动幅度越小,其中,1号测点的最小,5号测点的最大。下面对其原因进行解释。当第1个压应力波的前缘传至1号测点不久,在巷道顶板表面发生反射;随后,由于1号测点离巷道顶板表面很近,反射的拉应力波立即传至1号测点,由于拉应力波σ3的绝对值与此时传至1号测点的压应力波的σ3相差应不大,二者叠加将造成压应力波的能量极大衰减,即对1号测点的影响大幅降低。当压应力波的前缘经由巷道顶板表面反射的拉应力波传至5号监测节点时,由于该拉应力波在有阻尼的介质中已传播了一定距离,其能量将有所衰减,传至该测点的压应力波的σ3的绝对值将大于该拉应力波的,即该拉应力波对5号测点的影响将不如对1号测点的大。

3.1.4 频率的影响

(1)对剪、拉裂缝区段数目的影响。

图8,9分别为方案1~4的剪裂缝区段数目Ns和拉裂缝区段数目NtN的演变规律,统计的NsNt包括图2~3中显示与否的裂缝区段。

由图8可以发现,各方案的Ns均呈现恒为0—近似阶梯增长—恒定的变化趋势。只标注了方案1的Ns-N曲线和Nt-N曲线。由图9可见,各方案的Nt均呈现恒为0—近似阶梯增长—恒定的变化趋势。

图8 方案1~4的Ns-N曲线
Fig.8 Evolution of the number of shear crack segments with time steps of Schemes 1-4

图9 方案1~4的Nt-N曲线
Fig.9 Evolution of the number of tensile crack segments with time steps of Schemes 1-4

在恒定为0阶段,第1个压应力波的前缘尚未传至巷道顶板表面。

在近似阶梯增长阶段,Nt时增时稳。以方案1(f=15 Hz)为例。当N=23 200~27 520时,Nt-N曲线存在若干个增长阶梯,期间,巷道两帮V形坑内产生拉裂缝。当N=27 521~46 300时,Nt稳定在280附近,期间,V形坑内拉裂缝几乎停止发展,同时,巷道顶板尚未产生拉裂缝。在N=46 301之后,Nt快速增长,期间,反射的多个拉应力波的累积作用不断提升顶板某一位置最大主应力的波峰,巷道顶板产生拉裂缝。

在恒定阶段,方案1~4的Nt分别最终稳定在2 681,1 110,741和546。期间,模型中不再有压应力波传入,拉裂缝停止发展。由此可见,随着f的增加,Nt的最终稳定值减小。

综上所述,方案1~4的NsNt的最终稳定值均随着f的增加而减小。已有研究表明,应力波的频率越高,在介质中传播衰减程度越大[25]。在本文模型中,采用局部自适应阻尼,应力波在其中传播,能量自然也会衰减。上述NsNt的最终稳定值依赖于f的数值结果可以在一定程度上得到解释。

(2)对平均弹射初速度的影响。

由图4,10可以计算得到,方案1~4的平均弹射初速度分别为8.43,8.56,10.14和11.07 m/s。由此可见,随着f的增加,平均弹射初速度增加。应当指出,方案1~4传入模型的压应力波分别为13,22,30和33个。传入的压应力波数量越多,则巷道围岩对即将弹射单元的推动作用越频繁,这会导致平均弹射初速度越大。

图10 方案2~4的监测单元的右下角节点vx-N曲线
Fig.10 Evolution of horizontal velocities of lower right corner nodes of monitored elements with time steps of Schemes 2-4

4 结 论

(1)鉴于单元畸变和局部自适应阻尼可能导致单元弹射出模型时速度失真,为了准确模拟弹射现象,在原始方法的基础上对弹射单元进行了刚化处理,确保了脱离模型的单元在碰撞前后作刚体运动。结果表明,弹射单元任一节点的水平速度呈现上升—稳定—衰减的变化趋势,这与不进行刚化处理的结果相比更符合实际。

(2)在压应力波传入模型之后,巷道顶板左、右对称线上大部分测点的最大主应力呈现近似正弦波动上升—衰减—稳定的变化过程。在某一节点分离前,最大主应力首先呈现有规律的波动;然后,波动幅度突然增大,并伴随着剧烈震荡。

(3)当应力波传入模型之后,顶板某一位置的最大主应力的波峰和波谷随着时间的增加而有增加的趋势,反射的多个拉应力波的累积作用不断提升该位置最大主应力的波峰,致使拉裂,此即为巷道顶板拉裂机理。

(4)周期载荷频率对单元弹射初速度和剪、拉裂缝区段数目有一定影响。随着频率的增加,脱离巷道围岩的单元的平均弹射初速度增加,剪、拉裂缝区段数目减小。

参考文献(References):

[1] 宫凤强,伍武星,李天斌,等.深部硬岩矩形隧洞围岩板裂破坏的试验模拟研究[J].岩土力学,2019,40(6):2085-2098.

GONG Fengqiang,WU Wuxing,LI Tianbin,et al.Simulation experimental study of spalling failure of surrounding rock of rectangular tunnel of deep hard rock[J].Rock and Soil Mechanics,2019,40(6):2085-2098.

[2] 李夕兵,陶明,宫凤强,等.冲击载荷作用下硬岩层裂破坏的理论和试验研究[J].岩石力学与工程学报,2011,30(6):1081-1088.

LI Xibing,TAO Ming,GONG Fengqiang,et al.Theoretical and experimental study of hard rock spalling fracture under impact dynamic loading[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(6):1081-1088.

[3] XIA M,GONG F Q.Effects of loading waveforms on rock damage using particle simulation method[J].Journal of Central South University,2018,25(7):1755-1765.

[4] 朱哲明,李元鑫,周志荣,等.爆炸荷载下缺陷岩体的动态响应[J].岩石力学与工程学报,2011,30(6):1157-1167.

ZHU Zheming,LI Yuanxin,ZHOU Zhirong,et al.Dynamic response of defected rock under blasting load[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(6):1157-1167.

[5] 王之东,黎立云,陈滔,等.矿柱岩爆模型试验中能量释放研究[J].岩土力学,2018,39(2):177-185,208.

WANG Zhidong,LI Liyun,CHEN Tao,et al.Study of energy release in model tests on pillar rockburst[J].Rock and Soil Mechanics,2018,39(2):177-185,208.

[6] 陈滔,黎立云,邓建辉,等.岩爆问题中岩块弹射速度的数值计算与实验研究[J].四川大学学报(工程科学版),2014,46(1):26-31.

CHEN Tao,LI Liyun,DENG Jianhui,et al.Numerical simulation and experimental research on rock ejection velocity in rockburst[J]. Journal of Sichuan University(Engineering Science Edition),2014,46(1):26-31.

[7] 秦剑峰,卓家寿.岩爆问题中块体速度探讨[J].岩土力学,2011,32(5):1365-1368.

QIN Jianfeng,ZHUO Jiashou.A discussion on rock velocity in rockburst[J].Rock and Soil Mechanics,2011,32(5):1365-1368.

[8] 李夕兵,廖九波,赵国彦,等.动力扰动下高应力巷道围岩动态响应规律[J].科技导报,2012,30(22):48-54.

LI Xibing,LIAO Jiubo,ZHAO Guoyan,et al.Dynamic response of surrounding rock in highly-stressed.tunnel under dynamic disturbance[J].Science & Technology Review,2012,30(22):48-54.

[9] ZHU Z M,MOHANTYB B,XIE H P.Numerical investigation of blasting-induced crack initiation and propagation in rocks[J]. International Journal of Rock Mechanics & Mining Sciences,2007,44:412-424.

[10] 朱哲明,李元鑫,周志荣,等.爆炸荷载下缺陷岩体的动态响应[J].岩石力学与工程学报,2011,30(6):1157-1167.

ZHU Zheming,LI Yuanxin,ZHOU Zhirong,et al.Dynamic response of defected rock under blasting load[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(6):1157-1167.

[11] 金解放,李夕兵,常军然,等.循环冲击作用下岩石应力应变曲线及应力波特性[J].爆炸与冲击,2013,33(6):613-618.

JIN Jiefang,LI Xibing,CHANG Junran,et al.Stress-strain curve and stress wave characteristics of rock subjected to cyclic impact loadings[J].Explosion and Shock Waves,2013,33(6):613-618.

[12] 闫长斌,徐国元.动荷载对竖向排列地下硐室群稳定性影响分析[J].中国铁道科学,2006(3):27-33.

YAN Changbin,XU Guoyuan.Analysis on the influence of the dynamic loads on the stability for vertical arranged underground caverns[J].China Railway Science,2006(3):27-33.

[13] 陈国祥,窦林名,高明仕,等.动力挠动对回采巷道冲击危险的数值模拟[J].采矿与安全工程学报,2009,26(2):153-157.

CHEN Guoxiang,DOU Linming,GAO Mingshi,et al.Numerical simulation of dynamic vibration affecting rock burst in mining gateway caused by tremor[J].Journal of Mining & Safety Engineering,2009,26(2):153-157.

[14] 朱珍德,孙林柱,王明洋.不同频率循环荷载作用下岩石阻尼比试验与变形破坏机制细观分析[J].岩土力学,2010,31(S1):8-12.

ZHU Zhende,SUN Linzhu,WANG Mingyang.Damping ratio experiment and mesomechanical analysis of deformation failure mechanism on rock under different frequency cyclic loadings[J].Rock and Soil Mechanics,2010,31(S1):8-12.

[15] 葛科,赵玲,鲍东杰.岩石在循环加载下的力学特性[J].科技信息,2013(18):8-9.

GE Ke,ZHAO Ling,BAO Dongjie.Mechanical properties of rock under cyclic loading[J].Science & Technology Information,2013(18):8-9.

[16] 宫凤强,张乐,李夕兵,等.不同预静载硬岩在动力扰动下断裂特性的试验研究[J].岩石力学与工程学报,2017,36(8):1841-1854.

GONG Fengqiang,ZHANG Le,LI Xibing,et al.Experimental study on fracture behaviors of hard rock under dynamic disturbance with different pre-static loads[J].Journal of Mining & Safety Engineering,2017,36(8):1841-1854.

[17] SU G S,JIANG J Q,ZHAI S B,et al.Influence of tunnel axis stress on strainburst:An experimental study[J].Rock Mechanics and Rock Engineering,2017,50:1551-1567.

[18] SU G S,ZHAI S B,JIANG J Q,et al.Influence of radial stress gradient on strainbursts:An experimental study[J].Rock Mechanics and Rock Engineering,2017,50(10):2659-2676.

[19] 郭翔.基于势接触力的二维连续-非连续方法研究[D].阜新:辽宁工程技术大学,2017.

GUO Xiang.Studies of the two-dimensional continuum-discontinuum method based on potential contact force[D].Fuxin:Liaoning Technical University,2017.

[20] 杜卫东.煤矿岩爆事故防治技术与管理[J].能源技术与管理,2015,40(4):126-128.

DU Weidong.Prevention and control technology and management of coal mine rock burst accidents[J].Energy Technology and Management,2015,40(4):126-128.

[21] 陈建君,马鹏,余伊河,等.动载方向对巷道冲击影响的数值模拟研究[J].煤,2013,22(11):3-5,49.

CHEN Jianjun,MA Peng,YU Yihe,et al.Numerical simulation study on the effect of roadway burst induced by dynamic disturbance at different directions[J].Coal,2013,22(11):3-5,49.

[22] 韩晓亮.爆破动载作用下出矿巷道稳定性研究[D].赣州:江西理工大学,2016.

HAN Xiaoliang.Study on the stability of ore-drawing roadway based on blasting dynamic load[D].Ganzhou:Jiangxi University of Science and Technology,2016.

[23] 董锁堂.动载作用下巷道围岩变形试验研究[D].邯郸:河北工程大学,2016.

DONG Suotang.Experimental study on deformation of surrounding rock of roadway under dynamic loading[D].Handan:Hebei University of Engineering,2016.

[24] 冯夏庭,肖亚勋,丰光亮,等.岩爆孕育过程研究[J].岩石力学与工程学报,2019,38(4):649-673.

FENG Xiating,XIAO Yaxun,FENG Guangliang,et al.Study on the development process of rockbursts[J].Journal of Mining & Safety Engineering,2019,38(4):649-673.

[25] 陈颙,黄庭芳,刘恩儒.岩石物理学[M].合肥:中国科学技术大学出版社,2009.