基于SPH-FEM 耦合算法的磨料水射流破巖數值模擬

林曉東，盧義玉，湯積仁，敖 翔，張 磊

（1．重慶大學 煤礦災害動力學與控制國家重點實驗室，重慶 400030；2．重慶大學 復雜煤氣層瓦斯抽采國家地方聯合工程實驗室，重慶 400030）

磨料水射流是磨料與高速流動的水，或者與高壓水互相混合而形成的液固兩相介質射流，因其切割破碎作業效率高、作業過程沒有熱反應區、不發生化學反應等優點，被廣泛運用在石油化工、機械加工、采礦、隧道開挖等行業中［1－4］。尤其在巖石破碎領域，文獻［5］提出磨料水射流聯合機械破碎硬巖方法，取得了較好的效果。在該方法中磨料水射流除了沖蝕去除巖石外，其對沖孔周邊巖石的損傷對后續的機械破巖具有重要意義。

磨料射流破巖是一個涉及諸多因素的非線性沖擊動力學問題，具有瞬時強值動載荷、大變形及高應變率等特點。受現有理論研究能力和試驗條件的限制，利用理論方法或實驗手段進一步探索磨料射流沖蝕損傷破巖機理的難度非常大。然而，隨著計算機技術和計算理論的發展，可以運用數值模擬計算手段對上述問題進行分析研究，這一方法已在大量復雜問題的求解中取得成功。孫清德等［6］運用動態非線性有限元法和Hoffman破碎準則對高壓水射流破巖規律進行了模擬，認為提高射流沖擊速度、噴射直徑、橫移速度、射流束數，在35°至40°之間合理選擇射流入射角可以提高射流破碎巖石的效率；宋祖廠等［7］運用SPH算法模擬了高壓水射流破巖過程，針對水射流破巖能量轉化、水射流沖擊力、水力破巖演化進行了分析；盧義玉等［8］運用SPH算法模擬了脈沖水射流破巖過程中的應力波效應，得到了不同巖性巖石在脈沖射流應力波作用下的破壞形式；劉佳亮等［9］運用ALE算法對高壓水射流沖擊高圍壓巖石的損傷演化過程進行了數值模擬，認為圍壓對巖石軸向損傷影響較大徑向損傷影響較小。以上模擬方法，主要是針對高壓水射流的模擬。Kumar等［10］應用有限元法分析了單顆磨料粒子以不同角度沖擊鈦合金的損傷機理；王明波等［11－12］應用動態非線性有限元法對單顆粒磨料沖擊巖石破巖效果進行了模擬并對其破巖過程和破巖機理進行了分析，該方法沒有考慮到沖蝕過程是連續的磨料粒子的作用，同時也忽略了水的貢獻。Anwar等［13］應用有限元法模擬了磨料粒子束磨損Ti6Al4V材料的過程，并研究了磨痕的形成機理，但該方法忽略的水的作用。司鵠等［4］采用ALE流－固耦合罰函數算法模擬了磨料水射流破煤巖的過程，但只定性地分析了巖石的損傷場，并未對損傷范圍進行量化。

本文采用SPH-FEM耦合算法，將模型考慮成混有磨料的水柱與巖石作用，磨料水射流用SPH模擬，巖石采用Lagrange有限元模擬，充分考慮磨料和水對巖石的共同作用，與真實情形更加接近。數值模擬模擬磨料水射流沖擊作用下巖石沖蝕坑的形成過程，并研究巖石破碎坑周圍的損傷范圍。

1 SPH-FEM耦合算法

磨料水射流沖蝕損傷破巖過程中，高速射流出現高壓和大變形問題，采用傳統的拉格朗日法模擬磨料水射流容易出現網格畸變而導致計算終止［14］。SPH算法不使用單元，而是使用固定質量的可動點，所以沒有網格畸變。但是SPH方法相比拉格朗日法其計算精度不高，因此巖石采用FEM模擬，能夠得到相對準確的損傷特性。SPH-FEM耦合算法既可以很好的模擬瞬間大沖擊、大變形、高應變率等問題，又克服了SPH計算精度低的問題。

1.1 SPH基本理論

SPH方法的基礎是插值理論［15－17］。在SPH中任一宏觀變量都能方便的借助于一組無序點上的值表示成積分插值計算得到。核近似函數為：

式（1）中，f（x）為三維坐標向量 x的函數；Ω為點 x的支持域；x－x′為粒子間距離；h為SPH粒子的光滑長度，光滑長度隨時間和空間變化；W（x－x′，h）為核函數通常使用輔助函數θ（x－x′）進行定義，即

式（2）中，d為空間維數。

SPH中最常用的光滑核是三次B樣條曲線，即

式（3）中，C為歸一化常量，由空間維數確定

聯系式（1）（2）用粒子近似方法將連續形式的積分方程轉換成離散型的方程形式為：

式（4）中，ρi為粒子 i的密度；mi為位于 i處的粒子質量。

1.2 SPH-FEM耦合處理

SPH粒子在LS-DYNA中被視為特殊的節點單元，控制參數為節點編號、質量以及空間位置［14］。SPH粒子與FEM耦合通過罰函數方式將質點的力作用于有限元單元表面，因此SPH-FEM耦合采用點面接觸的接觸方式。在本文SPH與FEM耦合處理中，將SPH粒子定義為從節點，將與SPH粒子接觸界面上的有限元單元表面定義為主面。接觸耦合算法如圖1所示［18］。

圖1 SPH-FEM耦合算法Fig．1 Coupling algorithm of SPH-FEM

2 計算模型描述

為了對問題進行一定的簡化，在建立模型之前作如下假設：①磨料水射流破巖過程僅僅涉及水、磨料、巖石三種物質；②水粒子與磨料粒子為等直徑的球體；③磨料粒子隨機分布在水粒子之間，磨料粒子與水粒子具有相同的速度；④巖石為連續介質體，忽略孔隙介質的影響。

2.1 磨料水射流模型

磨料水射流破巖過程中，水、磨料均采用SPH算法模擬，在LS-DYNA中采用NULL材料模型，并對水、磨料材料模型賦予 Mie-Grueisen狀態方程（式（5）），其中磨料選用陶粒。水與磨料本構模型參數見表1。

式（5）中，E是單位體積內能，C是vs－vp曲線的截距，s1、s2和 s3是 vs－vp曲線的斜率系數，γ0是 Mie-Grueis-en常數，a是一介體積修正量。

表1 水、磨料本構模型參數Tab．1 Parameters in M ie-Grueisen equation of state for water and abrasive

磨料水射流包含水與磨料兩種成分，如何實現磨料與水兩種成分的隨機混合是磨料水射流建模的關鍵。本文按照如下思路建立磨料水射流模型：

第一步，根據磨料水射流中磨料的質量濃度α確定磨料水射流中磨料與水的體積比β。

式（6～7）中，ma、mw分別為磨料水射流中磨料與水的質量；Va、Vw分別為磨料水射流中磨料與水的體積；ρa、ρ0分別為磨料水射流中磨料與水的密度。

第二步，根據假設建模所用的水粒子與磨料粒子為等直徑的球體，則由二者的體積比β可以確定出二者的個數比γ。

式（8～9）中，na、nw分別為磨料水射流模型中磨料粒子與水粒子的個數；d為磨料粒子與水粒子的直徑。

第三步，根據磨料水射流沖蝕靶距處的射流寬度（直徑）D［19］，確定單層模型直徑上的粒子的總數 Nd，LS-DYNA軟件由Nd便確定了單層SPH粒子總數Ns；根據模擬的需要確定磨料水射流模型長度L，確定模型層數Nc，確定SPH粒子總數N。

式（10）中，D0為磨料噴嘴出口直徑；x為沿軸向距噴嘴出口的距離，即靶距。

第四步，確定SPH粒子總數后，根據磨料粒子與na水粒子個數比γ確定磨料水射流模型中磨料粒子數量與水粒子的數量nw。

第五步，利用EXCEL程序隨機抽樣命令從SPH粒子中抽取數量為na的樣本，對其賦予磨料粒子的屬性，剩余的nw個粒子賦予水粒子的屬性。

第六步，將磨料粒子和水粒子的信息寫入 LS-DYNA關鍵字文件中，完成磨料水射流建模。

2.2 巖石模型

巖石采用H－J－C材料模型。由于H－J－C材料模型等效屈服強度是拉力、壓力、應變率及損傷的函數，而拉力和壓力是體積應變的函數，損傷累積是塑性體積應變、等效應變、拉力及壓力的函數，所以能夠滿足大變形、高應變率、高拉壓效應的巖石工況［20－21］。

H－J－C模型的強度模型以規范化等效應力描述：

損傷因子通過等效塑性應變和塑性體積應變累積得到：

式（14～15）中，σ*T／fc；σ為等效應力，σ*≤SMAX，SMAX為材料所能承受的最大強度；P為單元內的靜壓；T為材料的最大拉伸強度；ε·為應變率；ε0為參考應變率；fc為材料的抗壓強度；A、B、C、N、D1、D2為材料常數；D為損傷度（0＜D＜1），且 D1（P*＋T*）D2≥EFMIN，EFMIN為材料的最小斷裂應變；Δεp和Δμp分別代表在一個積分步長內單元的等效塑性應變和塑性體積應變。單元的變形分為抗壓和拉伸兩種情況。巖石的本構模型參數見表2，其中巖石材料為灰巖。

表2 灰巖的本構模型參數Tab．2 H－J－C parameters of limestone

2．3 幾何模型及邊界條件

根據以上建模方法建立如圖2所示幾何模型，由于磨料水射流破巖過程基本是軸對稱的，為了簡化計算，本次數值模擬建立1／4模型［4］。磨料水射流為Ф4 mm×70 mm的1／4圓柱體，并對對稱面采用 SPH＿SYMMETRY＿PLANE對稱約束。巖石模型的幾何尺寸為25 m×25 mm×25 mm，為了模擬巖石無限大的工況，對模型底面和外圍兩個面采用NON＿REFLECTING非反射約束，對稱面采用CONSTRAINED＿GLOBAL對稱約束。磨料水射流模型包含1 680個SPH粒子，其中磨料粒子（紅色）153個，水粒子（藍色）1 527個，巖石模型包含288 000個六面體有限元。

圖2 磨料水射流破巖模型Fig．2 Themodel of abrasive water jet breaking rock

3 計算結果分析

根據所述的建模方法對磨料水射流破巖過程進行模擬。圖3給出了速度為250 m／s、磨料濃度30%的磨料水射流破巖效果圖。

3.1 沖孔形態分析

圖3 磨料水射流破巖效果圖Fig．3 Effect picture of abrasive water jet breaking rock

圖4為射流速度250 m／s、磨料濃度30%不同時刻巖體沖蝕坑截面形狀演化圖。沖蝕坑形態的演化之所以會如圖4中所示，主要是因為：磨料水射流與巖石發生接觸，首先在沖蝕作用中心區域的巖石單元迅速破壞失效并刪除，形成初始孔徑，隨著沖蝕過程的延續，一部分粒子繼續加深沖蝕深度，與此同時與孔底撞擊后向周圍發散飛濺的另一部粒子對初始孔孔壁巖體進行沖蝕并擴大孔徑——擴孔，這便形成了“V”形剖面；但是用于擴孔的粒子能量低，孔徑很快便趨于穩定，隨著時間的進程，在沖蝕坑孔徑相對穩定的情況下深度卻在不斷增加，便形成了由成“V”形剖面和圓形截面組成的“子彈”體，如圖4（d）所示，與圖5實驗中沖孔形態ICT掃描結果以及文獻［9，22］中的結果均吻合較好，表明本文建立的磨料水射流沖蝕損傷巖石模型是正確可行的，可用于后續的磨料水射流破巖損傷特性的研究。

圖4 沖蝕坑截面形狀演化圖Fig．4 Evolution chart of erosion pit sectional shape

圖5 沖孔形態ICT掃描實驗結果圖Fig．5 Results of ICT scanning experiment

3.2 巖石損傷分析

磨料水射流破巖是一個“沖蝕－損傷－沖蝕－損傷……”交替進行的動態過程。磨料粒子對巖石表面的沖擊作用，在巖石中形成由徑向－中間裂紋和側向－擴展裂紋構成的裂紋系統，其中側向－擴展裂紋向巖石自由表面擴展，實現巖石的宏觀破碎，形成沖蝕坑，而徑向－中間裂紋則部分殘留于巖石之中，造成對沖蝕坑周圍巖石的損傷［23］。在磨料水射流聯合鉆頭鉆進的技術中，巖石沖孔后的損傷范圍尤其是沖蝕坑徑向損傷范圍對后續的機械巖石破碎影響重大。圖6給出了巖石在速度250 m／s磨料濃度30%磨料水射流沖蝕下不同時刻的宏觀損傷云圖。

圖6 巖石損傷變量D分布圖Fig．6 Distribution of damage variable D in rock

為了對磨料水射流沖蝕損傷范圍進行量化，如圖7所示，在X－Z剖面上從距巖石頂端自由面深4 mm距軸線4．5 mm處沿徑向方向由內及外依次選取相鄰的17個單元作為跟蹤測量點，給出各單元損傷值隨時間變化曲線，并繪制t＝200μs時刻各單元損傷值隨距離變化曲線。由圖8～9可知在磨料水射流沖蝕作用下處于沖蝕坑范圍內單元A損傷值迅速達到1并完全失效；由于沖蝕能量的衰減劇烈，處于沖蝕坑外的單元由內到外損傷值迅速下降，距離軸線足夠遠的單元K損傷值為0并未受到損傷，其余各單元均有不同程度的損傷。由此可以得出：① 單元A與單元B的交界面即為沖蝕坑的孔壁，其距沖孔軸心O的距離即為沖孔半徑，沖孔半徑約為0．45 cm；② 單元J與單元K的交界面即為沖蝕損傷區外邊界，損傷半徑約為0．90 cm。

圖7 X－Z剖面單元選取示意圖Fig．7 Schematic view of selected units in the X－Z cross-section

圖8 損傷因子D隨時間的變化關系Fig．8 D evolution with time

圖9 損傷因子D隨距離的變化關系Fig．9 D evolution with distance

按照上述方法模擬了濃度30%不同速度的磨料水射流沖蝕巖石，表3為不同速度磨料水射流作用下的沖蝕坑半徑和損傷半徑。由表3可知磨料水射流作用下巖石的損傷半徑與沖蝕坑半徑隨著射流速度減小而減小，兩者之比在1．8－2．2之間。

表3 不同速度磨料水射流作用下巖石損傷、沖蝕坑半徑Tab．3 Rock dam age，erosion pit radius under the action of abrasive water jet at different speeds

4 結 論

（1）針對磨料水射流沖蝕損傷破巖過程中因為高速射流高壓和大變形特征導致網格畸變而使計算終止的問題，提出SPH-FEM耦合算法，充分考慮了磨料粒子與水形成射流后共同對巖石作用的工況。磨料水射流采用SPH模擬，巖石采用FEM模擬，既解決了高速射流因為高壓和大變形導致模擬計算困難的問題，又充分利用了FEM計算精度高的特性。

（2）模擬分析了磨料水射流沖蝕作用下，巖體沖蝕坑截面形狀演化過程：首先在沖蝕作用中心區域形成初始孔徑，進而由于磨料水射流的擴孔作用形成“V”形剖面，隨著沖蝕過程的延續，沖蝕坑不斷的加深，最終形成了由成“V”形剖面和圓形截面組成的“子彈”體。

（3）采用H－J－C損傷材料模型，通過定義跟蹤測量點的方式研究了磨料濃度30%不同速度磨料水射流巖石的損傷范圍，得到磨料水射流作用下巖石的損傷半徑與沖蝕坑半徑隨著射流速度減小而減小，兩者之比在1．8－2．2范圍之間的結論。

（4）由于磨料水射流破巖作用機理復雜，實驗過程透明度差，本文提出的數值模擬方法旨在為研究磨料水射流破巖提供一種研究的方法，對于不同的磨料濃度和射流速度以及不同的巖性有待今后進一步的研究。

［1］李曉紅，盧義玉，向文英．水射流理論及在礦業工程中的應用［M］．重慶：重慶大學出版社，2007．

［2］牛繼磊，李根生，宋劍，等．磨料射流射孔增產技術研究與應用［J］．石油鉆探技術，2003，31（5）：55－57．NIU Ji-lei，LIGen-sheng，SONG Jian，etal．，Investigation and application of abrasive water jet perforation to enhance oil production［J］．Petroleum Drllng Technques，2003，31（5）：55－57．

［3］李根生，沈忠厚．高壓水射流理論及其在石油工程中應用研究進展［J］．石油勘探與開發，2005，32（1）：96－99．LIGen-sheng，SHEN Zhong-hou．Advances in seorches and appliecation ofwater jet theory in peeroleum enginearing［J］．Petrcleum Exploration and Development，2005，31（5）：96－99．

［4］司鵠，謝延明，楊春和．磨料水射流作用下巖石損傷場的數值模擬［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．

［5］Lu Y Y，Tang J R，Ge Z L，et al．Hard rock drilling technique with abrasive water jet assistance ［J］．International Journal of Rock Mechanics and Mining Sciences．2013，60：47－56．

［6］孫清德，汪志明，于軍泉，等．高壓水射流破巖規律的數值模擬研究［J］．巖土力學，2005，26（6）：978－982．SUN Qing-de，WANG Zhi-ming，YU Jun-quan，et al．A disquisition on breaking mechanismof high pressure jet impacting on rock［J］．Rock and Soil Mechanics，2005，26（6）：978－982．

［7］宋祖廠，陳建民，劉豐．基于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．

［8］盧義玉，張賽，劉勇，等．脈沖水射流破巖過程中的應力波效應分析［J］．重慶大學學報，2012，35（1）：117－124．LU Yi-yu，ZHANG Sai，LIU Yong，et al．Analysis on stress wave effect during the process of rock breaking by pulsed jet［J］．Journal of Chongqing University，2012，35（1）：117－124．

［9］劉佳亮，司鵠．高壓水射流破碎高圍壓巖石損傷場的數值模擬［J］．重慶大學學報，2011，34（4）：40－46．LIU Jia-liang，SIHu．Numerical simulation on damage field of high pressure water jet breaking rock under high ambient pressure［J］．Journal of Chongqing University，2011，34（4）：40－46．

［10］Kumar N，Shukla M．Finite elementanalysis ofmulti-particle impact on erosion in abrasive water jetmachining of titanium alloy［J］．Journal of Computer and Applied Mathematics，2012，263（18）：4600－4610．

［11］王明波，王瑞和，陳煒卿．單個磨料顆粒沖擊巖石過程的數值模擬研究［J］．石油鉆探技術，2009，37（5）：34－38．WANG Ming-bo，WANG Rui-he，CHEN Wei-qing．Numerical simulation study of rock breaking mechanism and process under abrasive water jet［J］．Petroleum Drilling Techniques，2009，37（5）：34－38．

［12］徐依吉，趙紅香，孫偉良，等．鋼粒沖擊巖石破巖效果數值分析［J］．中國石油大學學報（自然科學版），2009，33（5）：68－71．XU Yi-ji，ZHAO Hong-xiang，SUN Wei-liang，et al．Numerical analysis on rock breaking effect of steel particles impact rock［J］．Journal of China University of Petroleum，2009，33（5）：68－71．

［13］Anwar S，Axinte D A，Becker A A．Finite elementmodelling of overlapping abrasive water jetmilled footprints［J］．Wear，2013，303（1－2）：426－436．

［14］Wang JM，Gao N，Gong W J．Abrasive water-jetmachining simulation by coupling smoothed particle hydrodynamics／finite element method［J］．Chinese Journal of Mechanical Engineering，2010，23（5）：568－573．

［15］劉飛宏，王建明，余豐，等．基于SPH耦合有限元法的噴丸殘余應力場數值模擬［J］．山東大學學報（工學版），2010，40（6）：67－71．LIU Fei-hong，WANG Jian-ming，YU Feng，et al．Numerical simulation for compressive residual stress of shot-peening based on SPH coupled FEM［J］．Journal of Shandong University（Engineering Science），2010，40（6）：67－71．

［16］呂東喜，黃燕華，唐永健，等．基于SPH算法的磨粒沖擊工件表面過程數值模擬［J］．振動與沖擊，2013，32（7）：169－174．LUDong-xi，HUANG Yan-hua，TANG Yong-jian，et al．Simulating process of abrasive impacting a workpiece surface based on SPH method［J］．Journal of Vibration and Shock，2013，32（7）：169－174．

［17］紀沖，龍源，方向．基于FEM-SPH耦合法的彈丸侵徹鋼纖維混凝土數值模擬［J］．振動與沖擊，2010，29（7）：69－74．JIChong，LONGYuan，FANG Xiang．Numerical simulation for projectile penetrating steel fiber reinforced concrete with FEM-SPH coupling alforithm［J］．Journal of Vibration and Shock，2010，29（7）：69－74．

［18］ MA Li，BAO Rong-hao，GUO Yi-mu．Waterjet penetration simulation by hybrid code of SPH and FEA［J］．International Journal of Impact Engineering，2008，35（9）：1035－1042．

［19］張運．高壓水射流切割原理及其應用［J］．武漢工業大學學報，1994，16（4）：13－18．ZHANGYun-qi．High pressure waterjet cutting principle with application［J］．Journal of Wuhan University of Technology，1994，16（4）：13－18．

［20］張鳳國，李恩征．混凝土撞擊損傷模型參數的確定方法［J］．彈道學報，2001，13（4）：12－16．ZHANGFeng-guo，LI En-zheng．A method to determine the parameters of themodel for concrete impact and damage［J］．Journal of Ballistics，2001，13（4）：12－16．

［21］巫緒濤，李耀，李和平．混凝土HJC本構模型參數的研究［J］．應用力學學報，2010，27（2）：340－344．WUXu-tao，LIYao，LI He-ping．Research on the material constants of the HJC dynamic constitutive model for concrete［J］．Chinese Journal of Applied Mechanics，2010，27（2）：340－344．

［22］王建明，宮文軍，高娜．基于ALE法的磨料水射流加工數值模擬［J］．山東大學學報（工學版），2010，40（1）：48－52．WANGJian-ming，GONG Wen-jun，GAO Na． Numerical simulation for the abrasive water jet machiningbased on the ALE algorithm ［J］．Journal of Shandong University（Engineering Science），2010，40（1）：48－52．

［23］Lawn B R，Evans A G，Marshall D B．Elastic／plastic indentation damage in ceramics：the median／radial crack system［J］．Journal of the American Ceramic Society．2006，63（9－10）：574－581．