基于SPH算法的脉冲射流破岩应力波效应数值分析

司　鹄， 薛永志(1.重庆大学 资源及环境科学学院，重庆　400044；2.重庆大学 煤矿灾害动力学与控制国家重点实验室，重庆　400044)



基于SPH算法的脉冲射流破岩应力波效应数值分析

司鹄1,2， 薛永志1,2(1.重庆大学 资源及环境科学学院，重庆400044；2.重庆大学 煤矿灾害动力学与控制国家重点实验室，重庆400044)

摘要：采用光滑粒子流体动力学(SPH)方法分析脉冲射流破岩的动力学过程，可以很好地避免传统有限元方法(FEM)在处理大变形问题时的网格畸变问题。引入J-H-C模型利用SPH方法建立了脉冲射流破岩的数值计算模型。利用该方法模拟了脉冲射流在破岩过程中应力波的形成、传播及衰减的过程。通过计算分析应力波在时间上和空间上的变化规律，绘出了应力波在传播过程中的波形图，并分析了脉冲射流冲击煤岩的过程中应力波传播的时间及空间特性；另外还发现应力波在不同岩石介质中的传播特性是不同的，对不同岩石的破坏效果也有较大差异，结合岩石材料特性及波的传播特性对计算结果进行了分析解释。数值计算的结果与实际冲击情况相吻合，对脉冲射流破岩的工程应用具有一定的指导意义。

关键词：脉冲；水射流；应力波；SPH

与普通的连续水射流不同，脉冲水射流能通过水锤效应对材料的高穿透力和冲击作用从一开始就使材料碎裂并使裂纹迅速扩散，从而为人们开辟了一条破碎硬脆性材料的新思路。尤其值得一提的是，脉冲射流在靶体表面产生的冲击压力大约为一般连续射流的滞止压力的1.5倍～2.5倍，从而非常显著地减小了切割比能，脉冲水射流的这些优越之处受到了越来越多的研究人员的重视[1-3]。

目前关于脉冲射流破岩作用主要形成了以下几种学说：冲击作用、水锤作用、水楔作用、空化作用、疲劳破坏作用等[4-7]。对于脉冲射流应力波效应的研究仍具有较大空间。随着计算机硬件和软件的发展，数值模拟技术越来越受到人们关注。加之实验环境特殊，以及试验设备和经费的局限，使得用数值模拟的方法来研究脉冲射流破岩的应力波效应表现出了更好的经济性与可行性。目前数值模拟的方法已经对诸多复杂的物理过程进行了有效的模拟，在水射流破岩领域也已发展为一种成熟有效的研究手段。倪红坚等[8]根据连续介质力学和连续损伤力学和细观损伤力学理论，给出了高压水射流破岩系统中流体和岩石的控制方程及适用于水射流破岩全过程分析的岩石损伤模型。刘佳亮等[9]运用ALE算法对高压水射流冲击高围压岩石的损伤演化过程进行了数值模拟，认为围压对岩石轴向损伤影响较大径向损伤影响较小。Junkara等[10]将磨料射流简化为单颗磨料粒子以不同的速度、角度射向物体进行了试验分析与数值模拟，通过试验现象结果与数值分析结果的比较，得出了合理的最优入射角和最优入射速度。刘佳亮等[11]基于流-固耦合罚函数算法，从连续损伤力学和细观损伤力学角度出发，模拟了磨料水射流破碎岩石的过程，得出了岩石的破坏主要以冲击破坏为主，破坏具有明显的局部效应等结论。

本文采用光滑粒子流体动力学方法(SPH)建立了脉冲射流冲击岩石的数值模型，摒弃了有限元网格， 直接利用离散点来构造近似函数，能方便处理大变形和应力应变局部化等难题，更加接近脉冲射流破坏岩石的真实过程[12]。基于模型模拟了脉冲射流冲击岩石过程中应力波的产生、传播及消退的完整过程；突破性的得出了应力波在岩石介质中传播过程的波形图；通过多种分析方法结合的方式揭示了应力波在岩石中的传播特性。

1SPH算法控制方程

由于SPH算法不使用网格，所以能在拉格朗日式下处理大变形问题。SPH的基础是差值理论，在其中任一宏观变量均可由一组无序点上的只表示为积分差值得到。质点近似函数定义为[13-16]：

(1)

式中W是核函数，核函数W用辅助函数θ进行定义：

(2)

式中：d是空间维数，h是光滑长度，光滑长度随时间和空间变化。

SPH中常用的光滑核是三次B样条，定义为：

(3)

取初始时刻粒子α的坐标为xi在t时刻质点移动到另一位置，粒子α的坐标变为xj它是初始坐标的函数， 即物体运动Lagrange描述为:

xj=xj(xi，t)

(4)

在流体动力学中， SPH方法求解流体力学问题时主要采用如下方程来描述流体的运动和状态。

质点运动位置方程为：

(5)

质量守恒方程为：

(6)

动量守恒方程为：

(7)

能量守恒方程为：

(8)

式中，ρ(xi)为粒子i的密度，单位kg/m3；mj为粒子j的质量，单位kg；σαβ为应力张量，单位Pa；v(xi)为粒子i的速度，单位m/s；E(xi)为粒子i的单位质量内能，单位J；∏为人为黏性力，单位N；H为人为热流，单位J/s；W为核函数。

2几何模型及状态方程

建立脉冲射流破岩的三维几何模型。考虑到模型的对称性，实际建模中只建立1/4模型，对称面采用SPH_SYMMETRY_PLANE对称约束，如图所示。射流为φ4 mm×50 mm的圆柱体模型，划分为3 200个光滑粒子。岩石为50 mm×50 mm×30 mm的长方体模型，划分为192 000个光滑粒子。为切合实际工程中岩石无限大的情况，岩石模型周围采用非反射边界条件，底部采用全断面约束，流体采用自由边界[17-18]。

图1　脉冲射流破岩几何模型Fig.1 The model of pulsed water jet breaking rock

为了满足岩石大变形、高应变率、高拉压效应的工况，采用Johnson-Holmquist-Concrete非线性本构关系，强度以规范化等效应力描述如下[19-20]：

(9)

通过等效塑性应变和塑性体积应变累积得到的损伤因子表示如下：

(10)

式中，ΔεP是等效塑性应变增量；ΔμP为等效体积应变增量。

本文是用的岩石材料有煤岩、砂岩、花岗岩，相应主要力学参数如表1所示。

表1　岩石的材料参数

模型中水视为完全塑性材料，选择NULL材料模型，并赋予Gruneisen状态方程：

(γ0+αμ)E

(11)

式中：C是冲击波速度与质点速度变化曲线的截距；S1、S2、S3是冲击波速度与质点速度变化曲线的斜率系数；γ0为Gruneisen常数；a是γ0和μ=ρ/ρ0-1的一阶体积修正量，采用半经验半理论的公式表示如下：

(12)

射流的相关物理力学参数如表2所示[21]。

表2　流体参数

3数值结果分析

3.1脉冲射流冲击煤岩应力波的演化

图2所示为速度500 m/s的脉冲射流冲击煤岩时应力波的形成、传播、衰减过程以及岩石破碎坑形成的时序演化情况。图(a)为t=0 μs即冲击前的煤岩模型，在t=10 μs时，岩石表面出现洼坑，在岩石表面形成初始应力如图(b)所示，可以看出此时的应力值很大且影响范围比较集中，这是由于在岩石破坏前接触点周围的岩石变形处于持续累积中，能量也在不断积聚，所以应力值比破坏后时要大很多。之后应力波以射流与岩石的接触点为中心，呈球面波的形式向周围传播，影响范围逐步扩大，如(c)-(f)所示。同时不难看出：射流与岩石接触点的竖直方向相比水平方向的应力值更大，因为射流下方的岩石受到射流的直接挤压作用(岩石冲击后密度变化见图3)，势能积聚高于周围岩石，所以出现了这种现象。也正是由于该原因，岩石竖直方向上的破坏速度远大于水平方向，与实际工程中钻孔深度大于直径的情况也是相符的。从图2中(d)-(f)还可以发现，应力波出现明显的衰减现象，图2(f)表现得尤为明显，该现象是应力波在岩石中传播时，变形在粒子中的传递造成了能量损耗导致的。

图2　射流冲击煤岩过程中应力波的演化Fig.2 Timing sequences of stresswave for coal impacted by jet

图3　煤岩受冲击后密度变化Fig.3 Sequences of density for coal under impacted

该模型模拟的冲击孔的形状及尺寸与文献[6](见图4)及文献[14]的实验结果相符较好，岩石侧面的破碎裂纹与球面波型的破坏效果基本一致，很好地证明了该模型的科学性与可行性。

图4　岩石冲蚀破碎图Fig.4 Rock broken under pulsed jet

3.2应力波对脉冲射流速度的响应特征

为了研究脉冲射流的速度对应力波效应的影响特性，在控制脉冲射流长度、直径、耙距不变的情况下，分别用速度v=100 m/s，300 m/s，500 m/s，700 m/s的脉冲射流冲击相同材料参数的煤岩模型，观察应力波的特征。图5所示为t=30 μs时，各速度射流冲击下煤岩模型中应力波的分布状况。不难看出，射流冲击速度越大应力波的影响范围就越大，传播速度也越快。同时还可以发现，v=100 m/s时，岩石模型中应力波影响范围极小且岩石几乎没有发生破坏。这说明岩石的破坏存在一个门限应力值，该应力值对应一个门限射流速度，当射流速度小于一定值时，未达到煤岩的破坏强度，煤岩不会发生破坏。用相应速度的脉冲射流对煤岩相似材料进行冲击试验，结果如图6所示。v=100 m/s时，无明显破坏，随着射流速度的提高，材料破坏范围和破坏深度都随之增大，与模拟结果及推论相符较好。

图5　t=30 μs，应力波随射流速度的变化Fig.5 T=30 μs，effective stress evolution with velocity

图6　不同速度射流冲击煤岩相似材料破坏结果Fig.6 Damage evolution with velocity

统计不同速度脉冲射流冲击下煤岩模型所受应力的峰值，应力峰值与冲击速度之间的关系如图7中曲线a所示。由曲线可知，射流速度为v=100 m/s时，煤岩所受应力峰值为0.000 293 Mbar，射流速度为v=300 m/s，v=500 m/s及v=700 m/s时，应力峰值均为0.000 336 Mbar。根据岩石力学知识可知，煤岩在破坏前一直处于能量累计状态，应力值不断增加，当应力值达到破坏强度时发生破坏，此时积聚的能量得到释放，应力随之减小。由图6(a)可知，射流速度为v=100 m/s时岩石未能遭到破坏，应力峰值未达到煤岩的破坏强度。同时再由曲线中射流速度为v=300 m/s，v=500 m/s，v=700 m/s时的应力峰值均为0.000 336 Mbar可知，该三种情况下的煤岩应力峰值都已达到破坏强度，岩石发生了破坏，所以不难得出该煤岩模型的破坏强度约为σ=3.36×10-4Mbar。

图7　应力峰值随脉冲射流速度的变化情况Fig.7 Peak pressure evolution with jet velocity

综上所述，在脉冲射流冲击下，当射流速度不足以破坏煤岩时，应力峰值随冲击速度增大而增大；射流速度可以破坏煤岩后，应力峰值不再随射流速度的变化而变化，且该此时的应力值可以反映岩石的破坏强度。图7中的曲线b，c，d表示出了应力峰值与射流速度的几种可能关系。

3.3应力波对传播距离的响应特征

为准确的描述传播距离对应力波效应的影响，在所建模型中由近及远依次选择1-20个参考单元，分析它们在射流冲击过程中应力的发展规律。参考点1距射流轴线的距离为D=7.5 cm，后面每个单元距离逐个递增ΔD=0.625 cm，单元的具体分布见图8。

图8　参考单元分布图Fig.8 The position of selected elements

通过计算得出每一个选定位置上对应的参考单元在整个冲击过程中的应力变化，而后结合各参考点距射流轴线的距离能够绘得很好地反映应力波传播过程的位置效应的“位置-时间-应力”云图，见图9。

图9　煤岩模型位置-时间-应力云图Fig.9 Cloud chart of position-time-stress for coal model

位置-时间-应力云图可以直观的看出应力波在空间及时间上的变化关系。对云图分别在D=8.5 cm，D=10.5 cm，D=12.5 cm，D=14.5 cm处取截面，得到不同位置处参考单元的应力变化规律，如图10所示。观察曲线1，D=8.5 cm的参考单元在t=0 μs到t=150 μs的过程中，应力先快速增大，到达峰值后再迅速降为0。这与之前分析的能量的积累与释放过程相对应，最终降为0是由于该单元距射流轴线较近，岩石被破坏时该单元被冲垮导致的。再观察曲线2，3，4对应的参考单元，开始阶段由于应力波的传播与衰退，应力值也都出现先增大后减小的过程。由于这些单元距射流轴线较远，不在射流的直接破坏范围内，单元只受到了不同程度的损伤并未被冲垮，所以在应力波衰退后应力值并没有降低为零，而是基于射流冲击到不同深度时传播到该点的多个应力波相互交错等因素的干扰显现出了无规律的波动现象。

图10　脉冲射流冲击煤岩应力波位置效应图Fig.10 Presentation for Stress of elements with different position

对位置-时间-应力云图分别在t=20 μs，t=30 μs，t=40 μs及t=50 μs处取截面，得到不同时刻应力波在煤岩材料中传播的波形图，如图11所示。曲线1表示t=10 μs时的应力波波形，应力波波峰在D=7.5 cm处；随着冲击作用的进行，应力波也在不断向远处传播，t=20 μs到t=40 μs过程中应力波波峰不断转移，对应图中的曲线2，3，4。同时从图中还可以发现应力波振幅A1>A2>A3>A4，进一步印证了应力波的传播中的衰退现象。统计20个参考单元的应力波峰值，结合距射流中心的距离得出单元应力峰值随距离的关系如图11所示。可以发现单元峰值在开始阶段随着距离的增大快速下降；当与射流轴线的距离超过一定值时，由于应力受到的影响因素变得复杂，持续处于低值波动状态，所以其峰值的下降幅度也变得平缓。

图11　不同时刻的应力波波形图Fig.11 Stress wave pattern evolution with time

图12　单元应力峰值随距离变化关系Fig.12 Peak Stress of element evolution with distance

3.4应力波对传播介质的响应特征

为了研究介质对应力波传播的影响特性，用速度v=500 m/s的脉冲射流分别冲击煤岩、石灰岩和花岗岩三种材料，冲击结果如图13所示。图中列出了t=20 μs和t=30 μs时各岩石模型中应力波的分布情况。可以看出石灰岩和花岗岩中应力波的影响范围明显大于煤岩。 鉴于ρgranite>ρlimestone>ρcoal，花岗岩密度最大，粒子排列最紧密，变形在粒子之间传播时耗散的能量最低，所以应力波能传播到更远的距离。此外计算还得出花岗岩的峰值应力为σgranite=8.30×10-3Mbar，石灰岩的峰值应力为σlimestone=6.58×10-3Mbar，煤岩的峰值应力为σcoal=3.36×10-4Mbar，反映出了三者的破坏强度之间的关系为：σgranite>σlimestone>σcoal，该结果和实际工程中的材料参数关系也是相吻合的。

图13　应力波在不同岩石中的传播Fig.13 Stress wave in different rocks

4结论

(1) 采用光滑粒子流体动力学方法(SPH)建立了脉冲射流冲击煤岩的数值模型，摒弃了有限单元的方法，避免了出现网格畸变的问题，很好地对脉冲射流冲击煤岩过程中应力波的传播问题进行了研究，模拟效果和实际试验结果相符较好，说明模型具有较好的科学性与合理性。

(2) 模拟研究了脉冲射流的速度对应力波的影响特性，结合试验发现随着射流速度的增加，应力波的影响范围随之增加，应力峰值出现先增加后不变的发展趋势；同时还发现当脉冲射流足以破碎岩石时，应力峰值和射流参数无关，只与岩石材料有关，且该应力峰值能够反应岩石的破坏强度。

(3) 模拟得出了应力波与传播距离之间的关系，同时结合选取单元的位置信息得到了不同时刻应力波的波形图；图形显示了应力波传播过程中峰值变化趋势，直观的反映出了应力波在传播过程中的衰退规律。

(4) 应力波的传播特性与介质的材料属性相关；脉冲射流分别冲击煤岩、石灰岩和花岗岩的结果显示：随着岩石密度的增加，应力波传播范围随之增加，应力峰值也将变大；同时应力峰值反映出岩石的破坏强度关系，即σgranite>σlimestone>σcoal，与实际相符。

参 考 文 献

[ 1 ] 胡寿根,丁胜．脉冲高压水射流工作原理及研究现状[J]．华东工业大学学报, 1997, 19(2):1-9．

HU Shou-gen, DING Sheng．The principle and status of high pressure pulsed water jet[J]．J East China University of Technology, 1997, 19(2):1-9．

[ 2 ] 李晓红, 卢义玉, 向文英．水射流技术及在矿业工程中的应用[M]．重庆:重庆大学出版社, 2007．

[ 3 ] 胡东，唐川林，张凤华．脉冲气液射流冲蚀特性实验分析[J]．振动与冲击，2013，32(11)：141-144.

HU Dong，TANG Chuan-lin，ZHANG Feng-hua. Erosion characteristic of a pulsed air-water jet[J].Journal of Vibration and Shock，2013，32(11)：141-144.

[ 4 ] 田方宝, 林缅．水射流辅助破岩机理研究(1):气泡空蚀[J]．力学与实践, 2007, 29(1):29-33．

TIAN Fang-bao, LIN Mian．Studies on the mechanism of water jet-assisted drilling technology(1):cavitation and erosion[J]．Mechanics in Engineering, 2007, 29(1):29-33．

[ 5 ] 卢义玉, 李晓红, 向文英．空化水射流破碎岩石的机理研究[J]．岩土力学, 2005, 26(8):1233-1237．

LU Yi-yu, LI Xiao-hong, XIANG Wen-ying．Rock erosion mechanism of cavitating water jets[J]．Rock and Soil Mechanics, 2005, 26(8):1233-1237．

[ 6 ] 田方宝, 林缅．水射流辅助破岩机理研究(2):水滴撞击[J]．力学与实践, 2007, 29(2):34-39．

TIAN Fang-bao, LIN Mian．Studies on the mechanism of water jet-assisted drilling technology(2):cavitation and erosion[J]．Mechanics in Engineering, 2007, 29(2):34-39．

[ 7 ] 卢义玉, 张赛, 刘勇, 等．脉冲水射流破岩过程中的应力波效应分析[J]．重庆大学学报, 2012．35(1):117-124．

LU Yi-yu, ZHANG Sai, LIU Yong, et al．Analysis on stress effect during the process of rock breaking by pulsed jet[J]．Journal of Chongqing University, 2012．35(1):117-124．

[ 8 ] 倪红坚, 王瑞和, 张延庆．高压水射流作用下岩石破碎机理及过程的数值模拟研究[J]．应用数学和力学, 2005, 26(12):102-105．

NI Hong-jian, WANG Rui-he, ZHANG Yan-qing．Numerical simulation study on rock breaking mechanism and process under high pressure water jet[J]．Applied Mathematics and Mechanics, 2005, 26(12):102-105．

[ 9 ] 刘佳亮, 司鹄．高压水射流破碎高围压岩石损伤场的数值模拟[J]．重庆大学学报:自然科学版, 2011, 34(4):40-46．

LIU Jia-liang, SI Hu．Numerical simulation on damage field of high pressure water jet breaking rock under high ambient pressure[J]．Journal of Chongqing University:Natural Science Edition, 2011, 34(4):40-46．

[10] Junkara M, Jurisevica B, Fajdigab M, et al. Finite element analysis of single-particle impact in abrasive water jet machining[J]．International Journal of Impact Engineering, 2006, 32 (2): 1095-1112．

[11] 司鹄, 谢延明, 杨春和．磨料水射流作用下岩石损伤场的数值模拟[J]．岩土力学, 2011, 32(3):935-940．

SI Hu, XIE Yan-ming, YANG Chun-he．Numerical simulation of rock damage field under abrasive water jet[J]．Rock and Soil Mechanics, 2011, 32(3):935-940．

[12] 宋祖厂, 陈建民, 刘丰．基于SPH算法的高压水射流破岩机理数值模拟[J]．石油矿场机械, 2009, 38(12):39-43．

SONG Zu-chang, CHEN Jian-min, LIU Feng．Numerical simulation for high pressure water jet breaking rock mechanism based on SPH algorithm [J]．Oil Field Equipment, 2009, 38(12):39-43．

[13] 刘更, 刘天祥, 谢琴．无网格法及其应用[M]．西安:西北工业大学出版社, 2005．

[14] 李裕春, 时党勇, 赵远．ANSYS10．0/ LS-DYNA 理论基础与工程实践[M] ．北京:中国水利水电出版社, 2006 ．

[15] 林晓东, 卢义玉, 汤积仁,等．基于SPH-FEM 耦合算法的磨料水射流破岩数值模拟[J]．振动与冲击, 2014, 33(18):170-176．

LIN Xiao-dong, LU Yi-yu, TANG Ji-ren,et al．Numerical simulation of abrasive water jet breaking rock with SH-FEM coupling algorithm [J]．Journal of Vibration and Shock, 2014, 33(18):170-176．

[16] Libersky L D, Petschek A G．High Strain Lagrangian Hydro dynamics [J]． JCP, 1993, 109:67-71．

[17] Holmquist T J, Johnson G R．A computational constitutive model for concrete subjected to large strains, high strain rates, and high pressure [J]．Journal of Applied Mechanics, 2011, 78(5):1-9．

[18] Holmquist T J, Templeton D W, Bishnoi K D．Constitutive modeling of aluminum nitride for large strain high-strain rate and high-pressure applications[J]．Computer Methods in Applied Mechcanics and Engineering, 2006, 1995:110-132．

[19] 巫绪涛, 李耀, 李和平．混凝土H-J-C本构模型参数的研究[J]．应用力学学报, 2010, 27(2):340-344．

WU Xu-tao, LI Yao, LI He-ping．Research on the material constants of the H-J-C dynamics constitutive model for concrete [J]．Chinese Journal of Applied Mechanics, 2010, 27(2):340-344．

[20] 王建明, 宫文军, 高娜．基于ALE法的磨料水射流加工数值模拟[J]．山东大学学报:工学版, 2010, 40(1):48-52．

WANG Jian-ming, GONG Wen-jun, GAO Na． Numerical simulation for the abrasive water jet machining based on the ALE algorithm [J]．Journal of Shandong University:Engineering Science, 2010, 40(1):48-52．

Numerical analysis for stress wave effects of rock broken under pulse jets

SIHu1,2,XUEYong-zhi1,2(1. College of Resources and Environmental Science, Chongqing University, Chongqing 400044, China;2. State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing 400044, China)

Abstract:Using the method of smooth particle hydrodynamics (SPH) to analyze the dynamic process of pulse jets’ breaking rock is an efficient way to avoid mesh distortion problems occurring in the traditional finite element method (FEM) to deal with large deformation problems. Here, a numerical calculation model of pulse jets’ breaking rock was built by introducing the J-H-C model and using the SPH method. The formation，propagation and attenuation of stress waves during the process of rock broken were simulated with this method. Through analyzing stress wave changing laws in time and space, stress wave patterns figures in the propagation process were drawn. The time domain characteristics and spacial characteristics of stress waves propagation during pulse jets’ breaking rock were analyzed. It was shown that the propagation characteristics of stress waves and their damage offect in different rock mediums are different. The calculation results were explained according to rock properties and wave propagation characteristics. The results of numerical simulation agreed well with those of the actual cases and provided a guide for engineering applications.

Key words:pulse; water jet; stress wave; SPH

中图分类号:TD231.62

文献标志码：A

DOI:10.13465/j.cnki.jvs.2016.05.023

通信作者薛永志 男，硕士生，1992年生

收稿日期：2015-01-27修改稿收到日期：2015-03-10

基金项目：国家自然科学基金资助(51274259);重庆市研究生科技创新项目资助(CYS14015)

第一作者 司鹄 女，博士，教授，1964年生