基于SPH-FEM耦合算法的后混合磨料水射流沖擊破巖數值模擬研究

米建宇，黃 飛,2,3，李樹清,2,3，王榮榮，李 丹

(1.湖南科技大學 資源環境與安全工程學院，湖南 湘潭 411201；2.湖南科技大學 南方煤礦頂板及煤與瓦斯突出災害預防控制應急管理部重點實驗室，湖南 湘潭 411201；3.湖南科技大學 南方煤礦瓦斯與頂板災害預防控制安全生產重點實驗室，湖南 湘潭 411201)

磨料水射流是一種將高速水射流與磨料相互混合而成的一種固、液兩相射流，目前被廣泛運用于破碎巖石、切割金屬等方面[1]。根據磨料摻合方式的不同，磨料水射流可分為前混合式和后混合式兩種。后混合磨料射流的磨料是在射流形成后通過獨立軟管在噴嘴內部完成混合，有效避免了在前混合磨料射流中高壓水管堵塞和高壓水管管路磨損嚴重等問題。后混合磨料是在高速水射流的推動牽引下獲得一定的加速，鑒于磨料與水的物理差異性，磨料粒子很難與水充分混合。為了獲得較高速度的沖擊磨料，通常采用提高系統泵壓、改變噴嘴結構等手段。然而，目前對于后混合磨料在噴嘴內的加速情況以及對破巖效果的影響規律方面的研究仍較為棘手。因此，開展后混合磨料水射流的磨料加速及其對破巖影響的數值模擬研究，對提高后混合磨料射流工作效率具有重要的工程意義。

磨料水射流破巖作為一個涉及到諸多因素的非線性沖擊動力學問題，其具有的瞬時強值動荷載、大變形以及高應變率等特點使得磨料水射流破巖難以從傳統理論以及實驗手段研究其沖蝕損傷及破巖機理。然而隨著計算機技術和工程商用軟件的發展，傳統理論和實驗手段難以解決的問題可以通過相關軟件的數值模擬對上述問題進行研究分析，這一方法已被大量學者運用并解決了諸多通過傳統理論手段所不能解決的問題。林曉東等[2]通過SPH-FEM(smoothed particle hydrodynamics with finite element method)耦合方法對前混合磨料水射流中磨料粒子加速進行了研究，發現了磨料粒子在噴嘴內部不同階段的加速特征及噴嘴結構對加速效果的影響規律。Li等[3]通過SPH-FEM耦合方法，從射流角度和施加圍壓方面出發，研究了射流傾角和圍壓對巖石破壞的影響。Wu等[4]通過對LS-DYNA程序的二次開發，建立了SPH-FEM/DEM數字高程模型模型，采用Voronoi細分技術研究了了在微觀結構和細觀力學性質下水射流破碎巖石的性能。表明了水射流沖擊巖石的破碎性能受顆粒大小和不規則性、延性、微觀強度、微觀參數的非均勻性等因素的影響較大，而受接觸剛度的影響較小。Ma等[5]通過研究基于光滑粒子流體動力學(smoothed particle hydrodynamics，SPH)法磨料水射流破碎巖石的動態過程，揭示了提高磨料濃度和射流速度可以大大提高磨料水射流破巖效果。Jiang等[6-7]通過SPH/FEA(smoothed particle hydrodynamics/finite element analysis)方法建立了水射流破巖的數值模型，并發現了巖石的裂紋和損傷的產生機制。劉佳亮等[8-9]通過對高壓水射流沖擊混凝土的研究，揭示了混凝土的致裂機理及裂紋擴展規律。文獻[10-11]通過對水射流中噴嘴流場的研究，揭示了噴嘴最佳收斂角度并對噴嘴進行了優化設計。

以往研究主要對前混合磨料射流加速以及其破巖規律開展了研究，對于后混合磨料射流破巖的數值模擬通常將磨料和水的速度設為相同，與實際情況不符。本文擬采用SPH-FEM耦合算法并限定SPH領域搜索范圍的方法來模擬后混合磨料射流的形成，進而研究射流速度、磨料濃度與圍壓大小等因素對后混合磨料水射流破碎巖石的影響規律。

1 基本算法描述

1.1 SPH基本原理

SPH法是無網格法的一種，其核心理論為差值理論[12-14]。對問題域采用離散化的粒子進行描述，有效避免了在非線性計算中單元畸變所導致的計算終止。其核近似函數為

(1)

式中：f(x)為核近似函數；W(x-x′,h)為光滑核函數；x-x′為粒子間距；h為粒子的光滑長度；Ω為點x處光滑核函數的支持域。由粒子近似可將連續SPH積分形式寫成離散化粒子近似式

(2)

Wij=W(xi-xj,h)

(3)

式中：f(xi)為粒子i處的近似函數值；mj為粒子j的質量(j=1,2,…,N)；N為在粒子j的支持域內粒子的總量；pj為粒子j的密度。

式(2)說明粒子i在任一函數值都可通過光滑核函數對其緊支域內所有粒子相對應的函數值進行加權平均近似。

1.2 SPH鄰域搜索

雖然無網格粒子法可以解決傳統的有限元網格大變形問題，但是其計算精度與計算效率較低等問題仍不容忽視[15-16]。當SPH磨料粒子在噴嘴內運動至收斂段、直線段時，截面粒子數量增多、粒子間距縮小，如圖1(a)所示。這將導致SPH鄰域搜索范圍內的粒子數急劇增多，進而急劇增加領域搜索時間。隨著柱塞的持續推進，SPH磨料粒子持續被壓縮，整個算例的計算時間將會無限延長，最終導致計算失效。現有的SPH數值模擬技術難以從算法上來提高模型的計算效率，為了節省計算時間同時不影響SPH粒子在噴嘴內的相互作用，本文將采用圖1(b)所示的方法對SPH粒子的鄰域搜索范圍進行限定，從而節省計算時間，提高計算精度。

圖1 鄰域搜索示意圖Fig.1 Schematic diagram of neighborhood searching

2 計算模型建立

2.1 磨料水射流模型

為了更好的模擬真實工況，首先將水射流假定為一段圓柱進行幾何建模，并在其中建立若干隨機分布的直徑為0.6 mm(30目)的磨料粒子。為防止水射流和磨料初始穿透，對水射流幾何模型和磨料粒子幾何模型進行布爾運算后分別進行有限元網格劃分，最后將水射流模型網格轉化為SPH粒子，保留磨料粒子的Lagrange有限元網格屬性。其中：水采用NULL材料模型，并賦予其Mie-Grueisen狀態方程；磨料粒子采用剛體模型。水與磨料的材料參數如表1與表2所示。其中：ρ0為材料密度；E為彈性模量；v為泊松比；φ為直徑。Mie-Grueisen狀態方程為

表1 水的本構模型參數Tab.1 The constitutive parameters for water

表2 磨料的模型參數Tab.2 Parameters of state for abrasive

(4)

式中：E為單位體積內能；C為vs-vp曲線截距；S1，S2，S3為vs-vp曲線斜率系數；γ0為Grueisen常數；a為一介體積修正量。

2.2 柱塞、噴嘴模型

將后混合磨料水射流考慮成通過柱塞以恒定的速度來推動純水，然后由運動的水來推動磨料粒子并形成初步混合，最后經過收斂噴嘴加速從而噴射出高速磨料水射流。本文中假設柱塞與噴嘴均不發生形變，故不考慮它們的厚度變化。因此，本文采用普通鋼制材料模型來模擬柱塞，而采用剛體材料模型來模擬噴嘴。柱塞與噴嘴均采用殼單元，并用Lagrange法進行有限元網格劃分。

2.3 巖石本構模型

巖石采用HJC模型，該模型綜合考慮了應變率效應、損傷演化效應、圍壓效應和壓碎、壓實效應影響，可以較好的描述混凝土、巖石材料在大變形、高應變率和高靜水壓力下的力學行為。該模型形式簡單、參數意義明確，現已被廣泛應用于沖擊爆炸等強動載作用下混凝土、巖石類材料的分析中[17-18]。該模型強度采用規范化等效應力描述

σ*=[A(1-D)+BP*N][1+Clnε*]

(5)

D=∑[(Δεp+Δμp)/D1(p*+T*)D2]

(6)

式中：σ*=σ/fc；σ*≤SMAX,SMAX為巖石所能承受最大強度；P*=P/fc；ε*=ε/ε0，ε為應變率；ε0為參考應變率，ε0=1.0 s-1；T*=T/fc；σ為等效應力；P為單元內靜壓；T為材料最大拉伸強度；A，B，C，N，D1，D2為巖石材料常數；D為損傷度，0≤D≤1.0，且D1(p*+T*)D2≥EFMIN,EFMIN為巖石最小斷裂應變；Δεp和Δμp分別為在一個積分步長內單元的等效塑性應變和塑型體積應變。在單元的變形中分為抗壓和抗拉兩種變形情況。fc為材料的抗壓強度，其余的巖石本構模型參數如表3所示。

表3 巖石本構模型參數Tab.3 Parameters of state for rock

2.4 接觸耦合處理

本模型中采用的獨立單元較多，有采用SPH粒子模擬的純水介質單元，還有采用Lagrange網格模擬的噴嘴、柱塞、磨料與巖石單元。為了實現柱塞推動水，水帶動磨料粒子經過收斂段加速后噴射，最終作用于巖石進行破巖的全過程，本模型中對水與磨料、水與噴嘴、水與柱塞均采用點對面單向接觸；涉及到磨料粒子與噴嘴及巖石的接觸碰撞，由于磨料整體較小，也采用點對面單向接觸；涉及到水和磨料整體對巖石的沖蝕，采用帶有侵蝕算法點對面單向接觸；磨料粒子之間的接觸采用自接觸定義。

2.5 幾何模型及網格劃分

根據以上建模方法建立的幾何模型，如圖2所示。圖中：柱塞直徑為14 mm并限制其自由度使其只在Z方向上運動；純水段直徑為14 mm，長度為50 mm；混合段直徑為14 mm，長度為50 mm；收斂段長度為20 mm；直線段長度為10mm；噴嘴出口直徑為3 mm；巖石為40 mm×40 mm×30 mm的長方體，根據林曉東和劉佳亮等的研究對其模型的設定，借鑒文獻中所述邊界條件及求解策略，對巖石周圍四面與底面施加非反射邊界，同時約束巖石的豎向移動。磨料粒子采用0粒、25粒、50粒、75粒、100粒5種濃度分別計算并分析對巖石的沖擊損傷效果。

圖2 后混合磨料水射流模型示意圖 (mm)Fig.2 The model of post-mixed abrasive water jet (mm)

3 后混合磨料水射流加速與破巖數值計算結果與分析

3.1 磨料粒子加速分析

根據2.1節～2.5節建模方法建立模型后進行初步計算，獲得水與磨料各階段的加速過程示意圖，如圖 3 所示。根據圖示，純水在柱塞的推動下帶動磨料向噴嘴收斂段流動，最終從噴嘴射出形成磨料水射流。由于本模型中所建立的磨料粒子并非單獨的一粒，且在磨料加速過程中涉及到磨料與噴嘴的碰撞和磨料間自身的碰撞，為了盡可能完整的提取到磨料粒子的加速過程，選取其中碰撞較少的一顆磨料粒子并同時選中與磨料處于同一斷面的SPH 粒子，提取兩者的速度并且繪制速度-時間曲線，如圖3所示。柱塞以一定速度作用于水，使得水在噴嘴混合腔內瞬間獲得一個較高的速度，與此同時高速水射流與混合于其中的靜止的磨料發生相互作用，強烈的速度差導致磨料在混合腔內發生劇烈的加速效果[19-20]；當磨料運動到噴嘴收斂段時，收斂段積聚的能量被部分釋放，使得磨料在收斂段被繼續加速；由于磨料粒子的質量較大，提速較為緩慢，在進入直線段后仍與水有較大的速度差，且直線段過流斷面較小，使得磨料粒子繼續加速并逐漸趨近于水，最終在噴嘴出口磨料的速度加速至純水速度的 80%左右，并與水一起作用于巖石；當水與磨料分別撞擊到巖石后，動能用于破巖能量消耗，各自的速度急劇減小。

圖3 水與磨料的各階段加速過程與速度變化規律Fig.3 The accelerating process of water and abrasive and their velocities’changing law

3.2 磨料水射流破碎巖石損傷規律

柱塞速度為300 m/s、磨料濃度為75粒的巖石沖蝕截面形狀演化圖，如圖4(a)所示。其原理是：柱塞作用于水射流，給予水射流一定速度，水射流帶動混合與其中的磨料經過收斂噴嘴從而噴射出高速磨料水射流并與巖石發生接觸；首當其沖的巖石中心區域單元被迅速破壞失效并刪除，形成初始孔徑。隨著沖蝕的進行，磨料水射流的沖擊使得孔徑逐漸加深，與此同時與底部接觸的粒子反彈而后向周圍飛濺對周圍孔壁進行沖蝕并擴大孔徑，由于反彈飛濺粒子能量較低，隨著沖蝕的進行，孔徑大小逐漸穩定并且深度不斷增加，形成了如圖4所示的“V”型剖面。圖4(b)與圖4 (c)為相同工況下磨料水射流沖擊巖石的宏觀破碎形態與細觀CT掃描圖像。通過對比可知：本論文的數值模擬現象與沖擊破巖實驗中的巖石破碎形態幾乎一致，且與文獻[21-23]中的結果相符。由此間接表明本文所建立的后混合磨料水射流沖擊巖石模型恰當，可用于后續的后混合磨料水射流破巖的損傷效果研究。

圖4 沖蝕坑對比圖Fig.4 The comparison of erosion pit

3.3 磨料濃度對巖石破碎損傷的影響

為了探究磨料濃度對巖石損傷的影響，有限元法(finite element method，FEM)磨料采取0粒、25粒、50粒、75粒、100粒5種不同的濃度以柱塞300 m/s的速度分別進行計算，并提取在相同時間內不同磨料濃度的巖石損傷情況并繪制關于磨料濃度的破碎曲線圖，如圖5、圖6所示。

圖5 不同磨料濃度沖蝕坑截面演化圖Fig.5 Evolution of erosion pits of different abrasive concentration

圖6 破碎效果與磨料濃度的關系曲線Fig.6 The relationship of broken effect vs.abrasive concentration

在3.1節～3.2節的分析中，磨料在噴嘴中加速后以趨近于水的速度和水射流一起作用于巖石，鑒于磨料顆粒大，慣性大，動能高，在磨料沖擊碰撞巖石后，產生的沖擊壓力與磨削力超過了巖石的抗壓強度，使得撞擊點附近巖石瞬間產生壓剪破碎，沖擊壓力所引發的環向拉應力及應力波反射拉應力使得巖石進一步破裂[24]。磨料在撞擊巖石后速度衰減并向四周反彈，逐漸形成沖蝕漏斗，隨著沖蝕的進行，破碎深度逐漸增加。從圖5可知：在沒有磨料的純水射流中，巖石的損傷程度較小，在加入25粒磨料后，沖蝕深度和寬度都有一個大幅度的提升，隨著磨料的增多，第一時間內作用于巖石的磨料也隨之增加，巖石的破壞深度也隨之加深[25]。不過在圖6中不難看出，當磨料在加入70粒以上后，破碎深度卻幾乎不再增加，而破碎寬度反而有一定減小。分析表明：磨料粒子的質量大、慣性大，在水帶動其加速的過程中存在較大的能量交換。隨著磨料的增加，水與磨料交換的能量增多，從而導致射流噴出時的總能量減弱，破巖的效果反而降低。因此可以推測，磨料的濃度需要控制在一個適宜的范圍以內，以保證磨料射流噴出時擁有更大的能量。在本數值模型中，磨料濃度在70～80粒內最佳。

3.4 射流速度對巖石破碎損傷的影響

在了解了磨料濃度對巖石破碎損傷影響之后，采用75粒磨料濃度的模型依次進行柱塞速度為100 m/s，150 m/s，200 m/s，250 m/s，300 m/s，350 m/s的模擬，并獲取在相同時間內不同速度下巖石損傷的效果并提取巖石的破碎深度及寬度并繪制關于速度的破碎情況曲線圖，如圖7、圖8所示。

圖7 不同速度沖擊下沖蝕坑演化圖Fig.7 Evolution of erosion pits under different velocity

根據文獻[26-27]中的闡述，在水射流沖擊巖石中，會產生水錘壓力和滯止壓力，其產生的水錘壓力是導致巖石破壞的主要因素。在距離沖擊中心的一定范圍內，有效應力隨徑向距離的增大而減小，最終趨于穩定。圖8所示的曲線表明：巖石的破碎深度與射流速度呈現近似線性增加的趨勢。這是由于隨著磨料水射流速度的增加，沖擊的能量增大，巖石的破碎深度也隨之增加。然而，巖石的破碎寬度隨著射流速度增加有一定程度的增大，隨后逐漸趨于穩定。分析表明：射流沖擊壓力隨著射流速度的增加而增加，巖石受沖擊點附近的應力也隨之增加，破碎寬度有一定程度增加。繼續增加射流速度，沖擊中心的應力增加了，但是射流的直徑并未有變化，中心處所產生的沖擊應力影響的范圍有限，破碎的寬度逐漸趨于穩定。由此推測，對于相同直徑噴嘴的磨料水射流，單一增加射流速度工況下，巖石的破碎寬度存在一個最大值，本次數值模擬的最大破碎寬度約為13 mm。

圖8 破碎效果與射流速度的關系曲線Fig.8 The relationship of broken effect vs.jet velocity

3.5 圍壓作用對巖石破碎損傷的影響

由于磨料水射流在應用于煤礦、石油鉆探等工程領域時，巖石大都處于高圍壓狀態，為了更加貼合真實情況，故在數值模擬計算時對巖石施加圍壓約束，通過模擬后磨料水射流沖擊破碎高圍壓巖石的過程，分析磨料水射流沖擊作用下高圍壓巖石的損傷演化過程，具有重要意義。為了更好的模擬真實工況，對巖石周圍4個面施加壓力以模擬在真實工況中有圍壓作用下的巖石破碎情況。依舊采用柱塞速度為300 m/s，磨料濃度為75粒的模型對巖石施加2 MPa，4 MPa，6 MPa，8 MPa，10 MPa的圍壓分別進行模擬計算，并繪制圍壓作用下巖石破碎情況曲線圖，如圖9與圖10所示。

圖9 不同圍壓作用下沖蝕演化圖Fig.9 Erosion evolution map under different confining pressures

圖10 破碎效果與圍壓關系曲線Fig.10 The relationship of broken effect vs.confining pressure

根據圖示可知：隨著圍壓的增加，巖石的破碎損傷效果出現減弱的趨勢。其中，破碎深度隨圍壓的增加呈現線性的減小趨勢，而破碎寬度隨圍壓先減小后趨于穩定；另外，隨著圍壓的增加，巖石損傷破碎口的形態由“V”形向“U”形轉變，這種變化趨勢與文獻[28]所得結果基本一致。經過分析，這種變化趨勢可能是由多種因素影響所導致的。首先，在磨料水射流的沖擊下，高速的磨料射流與巖石內部的孔隙流體之間形成巨大的壓力差，在壓力差的作用下流體向巖石內部滲透，從而引發“水楔作用”，該作用有利于巖石的進一步損傷破碎。在圍壓的作用下，磨料水射流的水楔作用會被減弱，且減弱的程度隨圍壓的增加而加強，不利于巖石的損傷破碎。其次，根據巖石強度理論原理可知，巖石的力學強度會隨著圍壓的增加而增大，從而不利于磨料水射流沖擊下巖石的損傷發展。再次，根據文獻[29-30]可知，圍壓會影響磨料水射流沖擊下巖石內應力波的傳播，從而減弱應力波對巖石的損傷作用，并增加巖石損傷所需的時間。

4 結 論

本文基于SPH-FEM 耦合算法，建立了后混合磨料水射流沖擊破碎巖石的有限元模型，并通過限定鄰域搜索來完成磨料水射流的加速，最終分析了磨料濃度、射流速度以及巖石圍壓對巖石損傷破碎效果的影響規律，研究結果表明：

(1)在柱塞的推動下，純水帶動磨料粒子分別在 混合腔、收斂段與直線段獲得加速后噴出形成磨料水射流，最終磨料噴出速度約為純水的 80%。

(2)巖石的損傷破碎程度隨磨料濃度的增加呈現先變大后減小的趨勢，在相同工況下存在一個最佳破巖的磨料濃度范圍，本數值模擬中最佳磨料濃度在70～80粒。

(3)巖石的破碎深度隨射流速度近似線性增加，而破碎寬度隨射流速度的增加呈現先增大后穩定的趨勢，最大破碎寬度約為13 mm。

(4)巖石的破碎深度隨圍壓的增加呈現線性的減 小趨勢，而破碎寬度隨圍壓先減小后趨于穩定；相對于無圍壓狀態，高圍壓狀態下巖石沖擊深度銳減5 mm，破碎深度減少1 mm，巖石損傷破碎口的形態由“V”形向“U”形轉變。