基于光滑粒子流方法的双向聚能张拉爆破数值模拟

2022-08-01朱淳胡南翔张晓虎何满潮陶志刚尹乾孟庆祥

中南大学学报（自然科学版） 2022年6期

朱淳，胡南翔，张晓虎，何满潮，陶志刚，尹乾，孟庆祥

(1.河海大学地球科学与工程学院，江苏南京，210098；2.中国矿业大学(北京)深部岩土力学与地下工程国家重点实验室，北京，100083；3.中国科学院武汉岩土力学研究所，岩土力学与工程国家重点实验室，湖北武汉，430071；4.贵州工程应用技术学院土木建筑工程学院，贵州毕节，551700)

在采矿和岩土工程中，普遍存在爆破开挖岩体和保护围岩完整性这一相互矛盾但又必须解决的问题。钻孔内的炸药不仅会对预开挖的岩体进行爆破，而且不可避免地会对围岩造成损伤和破坏。为了最大程度地降低爆破过程对围岩的损伤和破坏，便产生了定向爆破技术。其中，聚能爆破产生的冲击射流具有良好方向性，被广泛应用于各种地下工程中[1-3]。目前聚能控制爆破可归纳为2类[4]：

1)对孔壁施加不均布爆炸荷载，使得孔壁在预定形成断面的部位优先形成裂纹[5]；

2)在孔壁断面优先形成某种形状的缺口，使得裂纹在缺口处优先形成并发展[6-7]。

国内外学者们大量研究了岩石材料在爆炸荷载作用下的动态裂纹扩展行为，XU 等[8]运用有限元和离散元相耦合的方法(CCDM)研究了爆炸粒子在固体中的运动行为，发现爆炸过程存在2个粒子运动区域，即停滞区和运动区，这些区域会随着粒子速度改变而发生变化；ZHU 等[9-10]利用AUTO2D 程序模拟了圆柱体岩石中心爆炸过程中的裂纹扩张行为，发现了爆破孔周围会形成由剪应力导致的破碎区，而拉应力控制着径向裂纹的扩展，在距离边界较近的地方会因为应力波的反射形成环形裂纹；MA 等[11]利用LS-DYNA 软件模拟了爆破诱发岩石裂纹扩展的过程，发现加载速率对裂纹的形成具有显著影响，聚能装置可以较好地控制裂纹的扩展行为；CHO 等[12]通过实验和数值模拟研究了聚能爆破过程中裂纹的动力扩展行为，并分析了引导孔对裂纹扩展行为的影响；WEI 等[13]运用有限单元方法研究了爆破作用下岩石的损伤发展过程；NING等[14]运用非连续变形分析方法(DDA)模拟了爆炸过程块体的运动行为；DEHGHAN 等[15]利用ANSYS AUTODYN软件模拟了应力波导致的裂纹扩展过程；ONEDERRA等[16]利用Hybrid Stress Blasting Model 成功地预测了岩石在爆破过程中损伤区域的边界和形状。ZHU等[17]利用多物理场模型模拟了煤柱爆破损伤和瓦斯气体流动之间的相互作用过程；YILMAZ 等[18]运用三维有限差分方法模拟了岩体在爆破荷载下的力学行为，研究了不同应变率以及加载速率下损伤区域的变化规律；郭德勇等[19-21]在对双孔聚能爆破过程中的爆炸应力波叠加效应基础上，采用LS-DYNA 软件构建了双孔聚能爆破数值分析模型，研究了爆炸应力波的传播特征以及裂纹扩展机制；梁洪达等[22]采用数值软件建立了双向张拉聚能爆破模型，研究了不同聚能爆破模式下应力波的传播和裂纹开裂规律，发现联孔爆破的裂纹扩展长度更长；邓永兴等[23]运用LS-DYNA软件模拟了螺旋管聚能药包的破岩机理，证明了螺旋管药包具有良好的工程应用效果；吴波等[24]运用LSDYNA 软件对椭圆双极性聚能爆破结构进行了计算优化，发现药罩厚度为0.15 cm时，破岩效果最好；WANG 等[25]运用有限元和离散元相耦合的方法研究了含节理岩体在爆炸荷载作用下的裂纹扩展行为。上述模拟岩石爆破过程的数值方法可以归结为3种，分别是：

1)离散单元法(DEM)[8]；

2)基于连续介质假说的有限差分方法(FDM)和有限单元法(FEM)[9-24]；

3)离散元和有限元相结合的FEM-DEM方法[25]。

基于连续介质假说的方法更偏向于解决小变形问题，对于断裂、破坏和爆炸等大变形过程的模拟不够准确。离散单元方法适合解决非连续问题，但是模型的物理参数，如节理分布以及滑移准则等在实际情况中较难获得。近年来，一种无网格的光滑粒子流体动力学方法(SPH)可以处理固体和流体的大变形破坏问题，被迅速应用到各种大变形破坏过程的模拟中。光滑粒子流体动力学方法的基本思想是将计算区域内物质离散为一系列具有质量、速度和能量的粒子，然后通过核函数进行估值，从而求得不同时刻和位置时的各种物理力学量。由于这种方法不存在网格划分的过程，避免了在大变形过程中网格所产生的扭曲缠绕问题，因此，该方法存在较大优势。

随着第三次矿业科学技术变革在中国的发生，何满潮等[26-30]提出了切顶短臂梁理论，该理论的三大关键技术之一便是聚能爆破切顶技术。目前双向聚能爆破方法因为成本较高，难以将其大规模推广。为此，本文作者提出采用光滑粒子流方法研究双向聚能张拉爆破技术，通过模拟聚能爆破中裂纹萌生到扩展的过程，探究在60°，90°和150°聚能角度下聚能爆破裂纹的扩展特征，并与现场工程的爆破效果进行对比，证明了此方法的实用性和可行性，为巷道和隧道工程中岩体的定向爆破开挖提供了可靠的参考依据。

1 光滑粒子流方法(SPH method)

SPH 方法的核心是插值理论。在插值理论的基础上，首先，采用核函数近似方法将函数及其梯度用核函数的积分表示，把微分形式的连续体方程转化为积分形式的方程；随后，采用离子近似方法将连续形式的积分方程转化为离散形式的方程[31]。它通过一系列均匀分布、携带各种物理量的光滑粒子来求解各种边界下的偏微分方程，各质点的相互作用借助插值函数来描述，这些粒子之间不需要网格连接，而是某一粒子和周围2h范围内的粒子依据核函数估计值来发生作用，空间中任意域内质点r处的核估计值f(x)可表示为

式中：h为SPH粒子光滑长度，是核宽度的一种度量，决定核函数的区域影响半径；dV(r′)为r′位置粒子的体积；W(|r-r′|,h)为核函数，通常使用辅助函数θ(x)来进行定义，即

式中：d为空间维数；x=|r-r′|，为粒子间的距离，r粒子的矢量位置。SPH中最常用的光滑核函数是三次B样条曲线，即

式中：θ(u)为光滑核函数；D为归一化常量，由空间维数确定；

取初始时刻粒子α的坐标为x，它是初始坐标的函数，即物体运动的Lagrange描述为

在SPH 方法中，质量守恒方程、动量方程以及能量方程可表示如下：

式中：ρ为密度；m为质量；v为速度；E为能量；σ为应力张量；xi为第i个粒子的位置向量；α1和β1为空间指标；rij为空间2 个粒子i和j之间的距离；W为核函数；Aij为关于核函数的一个算子。

从上述方程可知，SPH 方法是一种纯拉格朗日方法，适用于求解高速碰撞等动态大变形问题。在SPH 方法中，求解精度虽然依赖于质点的排列规则度[32]，但相对于其他有网格的计算方法，对质点的排列要求远低于对网格的要求。因此，SPH方法可以避免网格大变形导致的精度损失问题，同时它可以处理不同介质边界的情况。

2 模型计算参数的确定

2.1 花岗岩模型参数

在本次模拟中，煤矿顶板考虑采用花岗岩替代常见的砂岩进行计算。虽然在煤矿中顶板多为砂岩，但是对于赋存条件好、较完整的砂岩顶板，其弹性模量、抗拉强度、剪切强度、抗压强度等与花岗岩的相关参数较接近。此外，砂岩的JH2模型对应的参数获得较困难，而花岗岩的JH2模型参数有较多的文献可以参考。表1所示为花岗岩基本力学参数。

表1 花岗岩基本力学参数Table 1 Basic mechanical parameters of granite

为了研究花岗岩在动态荷载下的应力以及裂纹扩展规律，利用Johnson-Holmquist 模型研究爆炸荷载作用下花岗岩聚能爆破过程的裂纹扩展过程。Johnson-Holmquist 模型适用于短脉冲强动载的条件，可模拟大应变、高应变率和高压条件下的损伤破坏问题。炸药的爆破过程是一种短脉冲、高强度加载问题，同时爆炸过程中也存在大应变、高应变率和高压条件的损伤破坏，因此，本文采用JH2模型结合光滑粒子流方法进行双向聚能张拉爆破过程的研究。另外JH2模型本身自带损伤失效模型，可方便直观地研究裂纹的扩展过程，JH2模型的相关参数详见文献[33]。

2.2 炸药参数

采用动力分析软件LS-DYNA进行聚能预裂爆破二维数值模拟，炸药采用2号岩石乳化炸药，炸药材料模型采用MAT_HIGH_EXPLO_SIVE_BURN 材料模型和JWL(Jones-Wilkins-Lee)状态方程[34]。

含能量形式的JWL-EOS方程如下：

式中：PJWL和V分别为爆轰产物的压力和相对比容；常数A和B分别为326 GPa 和5.8 GPa；常数R1，R2和w分别为5.81，1.56和0.57；能量密度E0为2.67×106kJ/m3。

2.3 聚能管材料参数

聚能材料选用PVC 管材，本构方程选用塑性随动模型[34]，其力学参数如表2所示。塑性随动模型是各向同性、随动硬化或各向同性和随动硬化的混合模型，其与应变率相关，可考虑失效。通过在0(仅随动硬化)和1(仅各向同性硬化)间调整硬化参数来选择各向同性或随动硬化。应变率采用Cowper-Symonds 模型，选用与应变率有关的因数σY表示屈服应力，如式(10)所示。

表2 聚能管材力学参数表Table 2 Mechanical parameter of cumulative pipe

式中：p和C为与应变率有关的常数；βh为硬化参数；σ0为初始屈服强度；EP为塑性硬化模量；εeffP为有效塑性应变。

2.4 聚能爆破物理模型

2.4.1 物理模型

双向聚能张拉爆破技术利用聚能装置在设定方向上产生张拉作用，致使岩石产生定向破裂。当炸药引爆后，双向聚能装置对爆轰产物产生瞬时抑制作用，爆轰过程产生的高温、高压和高速气体优先从聚能孔释放，促使径向初始裂纹萌生，气体进入到裂纹中，形成强有力的“气楔”。一方面，“气楔”会导致裂纹尖端不断延伸，在垂直聚能方向上产生反射拉应力集中，促使岩体沿设定方向张拉开裂；另一方面，聚能管对爆轰气体具有一定的抑制作用，减少了爆轰产物对孔壁周边岩体的破坏，抑制了非聚能方向裂纹的扩展。当使用双向聚能张拉装置进行联孔爆破或间隔爆破时，爆破应力波将在爆孔间产生叠加应力场，提高爆破孔间的拉张应力。若爆孔间距适当，爆孔间将形成光滑的定向控制破裂面，实现精确控制爆破的目的。图1给出了双向聚能张拉爆破原理。

2.4.2 双向聚能张拉爆破裂纹开裂准则

根据断裂力学理论，当裂缝端部应力强度因子KI大于岩石的断裂韧度KIC时开始起裂。图1(b)所示为双向聚能张拉爆破情况下的断裂力学模型，裂纹扩展过程中其尖端处的应力强度因子KI[34]为

图1 双向聚能爆破原理Fig.1 Principle of bidirectional cumulative blasting

式中：a0为初始聚能方向的裂纹深度；r为药卷半径；P为炸药粒子充满爆破孔时的压力；F为应力强度因子修正系数，F随裂纹扩展长度增大而增大，当a较小时，F变化较大，当a远远大于r时，F则变化较小且趋近于1。当聚能射流接触初始裂纹尖端，即a=a0时，裂纹尖端的应力强度因子仍然采用式(11)计算。

当以KIC表示岩石的断裂韧度时，将应力强度因子的公式稍作变形，即可得到裂纹起裂及扩展的准则为

当a=a0时，P的数值即为聚能方向张裂纹初始起裂对应的爆炸气体应满足的压力；当a≥a0时，裂纹起裂以后，后续爆炸气体将因为“气楔”作用进一步迫使裂纹扩展，而裂纹扩展导致爆炸气体压力下降，为保证裂纹持续扩展，爆炸气体瞬时压力需满足以上的起裂准则。同时，在非聚能方向上，a0=0，式(12)去掉分子中的第2 项即可得到在非聚能方向裂纹起裂的准则为

对比式(12)和(13)可知，非聚能方向上裂纹起裂的压力要大于聚能方向上的裂纹起裂压力。首先，爆生气体的射流侵彻作用有利于裂纹在聚能方向上优先发育，从而产生定向裂纹；其次，聚能方向裂纹的起裂压力相比非聚能方向的要小；最后，聚能装置在一定程度上抑制了非聚能方向上裂纹的萌生扩展。上述3点对于提高定向爆破的精度和控制围岩完整性具有良好的技术支撑作用。

3 不同聚能角度下爆破过程模拟

基于以上力学参数，结合现场聚能爆破的要求，利用LS-DYNA软件中的SPH求解器，模拟在3种聚能角度α(60°，90°，150°)，β为24°情况下单孔双向张拉聚能爆破的过程(α及β的取值示意如图1所示)，再现了聚能爆破过程中等效应力以及裂纹扩展的演化过程。

单孔情况下，针对聚能角度小于180°时定向开裂效果的相关研究较少，而这种角度下对于某些具有明显爆破开挖边界要求的工程具有非常重要的意义，同样煤矿顶板也会有高度起伏的情况，上述情况均需考虑聚能角度小于180°情况下的双向聚能张拉爆破的效果，因此本文α选用60°，90°和150°，考虑了大部分顶板的起伏极端情况，具有明显的工程意义。

本文作者在模拟β处于24°～160°之间单孔开裂过程的研究中发现，β越大，裂纹的线性越差，在聚能尖端产生分叉裂纹的概率越高，这种情况与实际情况下的爆破开裂效果不符。当β≤24°的情况下，可在聚能方向上生成定向单一裂纹，与实际较符合。β越小，定向开裂效果越好，但是对于光滑粒子流方法中粒子的配置要求更高，难以实现，因此本文β取24°，可方便、精准地模拟双向聚能张拉爆破的定向开裂过程。

以α=60°情况下双向拉伸聚能爆破过程为例，介绍其损伤裂纹扩展、粒子运动以及等效应力演化的过程。图2所示为α=60°情况下双向拉伸聚能爆破过程中岩体损伤裂纹演变过程云图。由图2可知，在爆破时间为0.02 ms 时，巨大的爆破压力均匀施加在聚能部位以外的岩石表面，致使此环状区域的岩石进入剪切破坏阶段，形成破碎区，同时在2个聚能方向形成2个较深的初始损伤裂纹。在0.03 ms 时，此2 个方向的裂纹进一步向深部发展，同时在338.1°(此角度为从爆孔中心水平向右出发的射线逆时针转动到裂纹位置所旋转的角度，后文中的裂纹角度均按此计算)方向也产生1 个裂纹。随着爆炸波不断由爆孔向外传播，在0.13 ms 时，这3 条裂纹几乎同时开裂延伸到距离爆孔中心400 mm处。在0.20 ms时，此3条裂纹基本上完成贯通。在0.50 ms 时，267.0°方向形成了1个较短的裂纹，同时爆孔局部区域的损伤发展几乎停滞，裂纹宽度也开始加大，同时等效应力也不断减小甚至消失。

图2 聚能角度为60°时聚能爆破岩体损伤裂纹演变过程云图Fig.2 Cloud diagram of damage cracks evolution in rock when cumulative angle is 60°

图3所示为α=60°时粒子局部运动的演变过程云图。由图3可以看出，在0.03 ms 之前，岩石基质裂纹损伤主要是因为爆炸波对聚能方向岩体的张拉以及对非聚能方向岩体的剪切作用；从0.20 ms 开始，炸药粒子已经开始进入到聚能管的聚能孔内；在0.50 ms 时，炸药粒子通过损伤裂纹进入到岩石基质中，此时岩石的损伤裂缝宽度也相应增加。这说明在聚能爆炸后期，炸药粒子或爆生气体会进入到损伤裂纹中，进一步促进裂纹在宽度和长度方向上扩展，这种加速致裂的作用称为爆生气体或炸药粒子的“气楔”作用。

图3 聚能角度为60°时聚能爆破粒子局部演变过程云图Fig.3 Cloud diagram of local particles evolution process when cumulative angle is 60°

图4所示为聚能角度为60°时聚能爆破岩体等效应力演变过程云图。由图4可以看出等效应力最大值始终发生在爆破孔周围附近，在0.02 ms 时，高压爆炸气体在爆破孔周围形成了高应力带，最高约为1.72 GPa。岩石介质在聚能方向形成了高应力集中点，这些点周围的岩石介质因为受到聚能的侵彻作用萌生了初始裂纹。随后爆炸波向周边扩散，在0.03 ms时，119.8°和181.5°这2个聚能方向上的2个初始裂纹已经开始明显扩展。同时，由于裂纹扩展使得能量持续释放，高应力带也随即减弱，在0.13 ms 时338.1°方向出现了部分低应力带。随着爆炸应力波继续传播，从0.20 ms 开始低应力带不断伸长，且宽度也随之增加。

图4 聚能角度为60°时聚能爆破岩体等效应力演变过程云图Fig.4 Cloud diagram of von mises stress evolution in rock when cumulative angle is 60°

图5所示为α=60°情况下聚能爆破过程中岩体产生裂纹的长度和角度关系曲线。由图5可以看出，3 条主裂纹的初始方向分别发生在2 个聚能方向(119.8°和181.5°)和非聚能方向(338.1°)的位置。在3 个主裂纹贯通时，119.8°方向裂纹向103.9°方向发生了偏移，偏移量为15.9°，同时181.5°裂纹向200.2°方向发生了偏移，偏移量为18.7°，而338.1°方向裂纹则偏移至339.1°，偏移量仅为1.0°。

图5 聚能角度为60°时聚能爆破岩体裂纹的扩展长度与角度关系曲线Fig.5 Relation ship curve between crack extension length and angle when cumulative angle is 60°

图6(a)和6(b)所示分别为裂纹扩展长度及角度随时间的变化曲线。由图6(a)可知，在0～0.20 ms之间，3 条主裂纹扩展速度趋于一致，达到了2.85 km/s。在0.20 ms 以后，3 条主裂纹的扩展速度大幅度降低，平均速度只有0.16 km/s，仅为初始扩展速度的5.6%。最终，3条主裂纹均实现了贯通，而267.0°方向的次生裂纹长度只发展至111 mm。由图6(b)可知，0.20 ms 后裂纹的开裂方向基本保持稳定，说明双向拉伸聚能爆破对于定向裂纹的形成具有较好的控制效果。

图6 聚能角度为60°时聚能爆破岩体裂纹的扩展特征曲线Fig.6 Propagation characteristic curve of crack in rock mass when cumulative angle is 60°

α=90°及α=150°情况下损伤裂纹的发展规律、炸药粒子的运动规律以及等效应力演化规律与α=60°情况下基本一致，3 种聚能角度下的裂纹扩展特征如表3所示。

由表3可知，不同聚能角度爆破下，在聚能方向上均有裂纹产生、扩展，并最终实现了贯通。聚能方向上裂纹的偏移量随着聚能角度增大而降低，α=60°时，聚能方向(119.8°和181.5°)的平均偏移量为17.30°，α=90°时，聚能方向(89.3°和178.8°)的平均偏移量为10.75°，α=150°时，聚能方向(91.5°和237.9°)的平均偏移量为6.55°。因此，随着聚能角度越大，双向聚能张拉爆破的精度相对越好。总体而言，在双向聚能张拉爆破后期，岩体裂纹均能沿着稳定方向扩展，而且对爆孔周围产生较小的损伤。

表3 不同聚能角度下裂纹扩展特征对比Table 3 Comparison of crack propagation characteristics under different cumulative angles

聚能方向裂纹的萌生是由于爆炸所产生高速射流粒子对接触围岩产生的侵蚀作用，而裂纹的进一步扩展，则是爆生气体在裂缝内的不断深入，以“气楔”作用的形式推动裂纹的发展。另外，在3种聚能角度下，均在聚能角平分线的反方向上出现了裂纹萌生和扩展，这是因为高压气体不断推进聚能方向裂纹扩展，也在相邻的围岩上施加了作用力，如图7所示，作用力分别为F1和F2，在聚能孔周边区域，必然会形成2 个相反方向的力偶，使得聚能孔周边区域以拉应力为主，而围岩的外侧边界则以压应力为主。随着聚能方向裂纹不断扩展，也会导致更多气体进入裂纹，同时F1和F2随之升高，最终聚能孔附近的拉应力不断变大，当超过围岩抗拉强度时，必然会形成开裂。

图7 非聚能方向开裂示意图Fig.7 Sketch map for cracking in non-cumulative direction

若F1和F2大小、方向和作用位置随着裂纹的开裂和偏移而变化较小，则在2个聚能方向角平分线的负方向与聚能孔的交点处附近必然拉应力最大。若此拉应力超过围岩的抗拉强度，则最终会在这个交点位置产生开裂，从而形成唯一1 条与2个聚能方向角平分线相反方向的伴生裂纹。

若F1和F2大小、方向和作用位置随着裂纹的开裂和偏移而变化较大，则裂纹不仅会在2个聚能方向角平分线的负方向与聚能孔的交点处附近生成，而且聚能孔附近的最大拉应力点的位置也会随着F1和F2的变化发生移动，裂纹也会在这些超出围岩抗拉强度的最大拉应力点生成，如在模拟90°聚能方向的爆破过程中，在另外1 个非聚能方向(285.1°方向)上也出现了裂纹扩展的现象。

4 现场应用

选取中煤集团姚桥煤矿某个工作面进行煤矿顶板定向预裂的现场试验，该工作面净高为2.8 m，采用走向长壁放顶煤、顶板全部垮落的采煤方法，姚桥煤矿试验现场地质参数如表4所示。

表4 姚桥煤矿试验现场地质参数Table 4 Geological parameters of experimental site in Yaoqiao coal mine

综合考虑各个岩层的岩性及厚度，根据理论研究和综合现场爆破材料情况，现场采用乳化炸药进行定向预裂爆破，以此探究60°，90°和150°聚能角度下单孔聚能张拉爆破的裂纹扩展特征。

由于聚能爆破后较难直接观测到聚能方向裂缝的完整形貌，因此现场实验中通常采用CXK6矿用钻孔窥视仪对聚能爆破裂缝进行监测，图8给出了不同角度单孔聚能爆破的裂缝成像效果图。由图8可知：在聚能爆破之后，在相应的聚能方向上生成了连续的裂缝，即生成了连续的破裂面，因为直接顶为极易垮落的泥岩，强度较低，所以聚能方向生成的裂缝呈曲线，且不光滑。60°和90°聚能爆破产生的裂缝走向较不规则，而150°聚能爆破产生的裂缝走向较规则，为光滑平直。总体而言，60°，90°和150°单孔聚能爆破均在聚能方向产生了定向裂缝，验证了双向聚能张拉爆破在不同聚能角度下良好的定向成缝效果。