基于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年生