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

2021-11-10王学滨刘桐辛钱帅帅

煤炭学报 2021年10期

王学滨，刘桐辛，田锋，钱帅帅

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

采矿、土木等工程中岩层周期破断、断层周期黏滑和反复爆破都能产生周期冲击载荷。周期冲击载荷是诱发岩爆、垮塌和片帮等灾害的重要因素。在周期冲击载荷作用下，巷道顶板和两帮会在极短时间内产生大量相互联通的裂缝，进而将巷道围岩切割成众多块体，它们会以猛烈方式脱离巷道围岩，造成灾害[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)所示。应当指出，模型的左、右侧面均为透射边界，即应力波经过时不会发生反射。

在巷道顶板，布置了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个单元边长的裂缝区段。

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为例，简单分析巷道顶板各测点的σ3随N的演化规律。图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和拉裂缝区段数目Nt随N的演变规律，统计的Ns和Nt包括图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的Ns和Nt的最终稳定值均随着f的增加而减小。已有研究表明，应力波的频率越高，在介质中传播衰减程度越大[25]。在本文模型中，采用局部自适应阻尼，应力波在其中传播，能量自然也会衰减。上述Ns和Nt的最终稳定值依赖于f的数值结果可以在一定程度上得到解释。

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

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