近壁聲空泡潰滅微射流沖擊流固耦合模型及蝕坑反演分析*

葉林征，祝錫晶，王建青

（中北大學機械工程學院，山西 太原 030051）

空化是由于液體環(huán)境壓力變化而產(chǎn)生的結(jié)合非線性動力聲學和極端物理學的動態(tài)現(xiàn)象，以空泡的反復形成和強烈潰滅為主要特征。19 世紀末，Parsons 和Barnaby 首次提出空化效應(yīng)，并指出空泡潰滅使螺旋槳葉片使用效率下降。隨后，大量學者在空化領(lǐng)域研究并取得了很多突破性進展，對空化效應(yīng)的研究重點也逐漸由抑制其破壞性轉(zhuǎn)移至利用其破壞性，使空化起到積極作用，如超聲清洗、石油鉆孔、海洋開發(fā)等領(lǐng)域[1]。

依據(jù)產(chǎn)生條件，空化主要分為超聲空化和水力空化，兩者分別通過聲壓載荷施加和特殊流動系統(tǒng)幾何結(jié)構(gòu)改變液體環(huán)境壓力，進而產(chǎn)生空化現(xiàn)象。水力空化在艦船領(lǐng)域十分常見，諸多學者對水下爆炸氣泡運動特性進行了廣泛研究。張阿漫等[2-3]運用高階三維邊界元法模擬了水下爆炸氣泡的動態(tài)特性，自主開發(fā)程序系統(tǒng)地研究了近自由面水下爆炸氣泡的動態(tài)特性。李健等[4]基于VOF 方法對水下爆炸氣泡運動進行了數(shù)值模擬，研究了距離參數(shù)對氣泡坍塌、射流及環(huán)狀氣泡的影響。郭婭等[5]考慮了水下近場爆炸沖擊波及射流的聯(lián)合作用，真實地模擬了水下近場爆炸過程。沈壯志等[6]研究了聲場中水力空化泡的動力學特性，認為超聲將水力空泡運動調(diào)制為穩(wěn)態(tài)空化，增強了空化效果。目前，超聲空化的研究重點集中在兩個方面：空泡動力學理論研究以及超聲空化對材料作用的試驗研究，前者重點集中在單泡、雙泡以及泡群動力學分析[7]，后者則重點研究超聲空化作用下材料的行為變化及試驗分析[8-9]，而聯(lián)系兩者間的紐帶—空泡潰滅微射流沖擊的研究卻少有報導。

眾所周知，超聲空化對材料的作用中，空泡潰滅微射流沖擊占主導地位。1944 年，Kornfeld 等[10]首先提出微射流理論，指出當空泡受壓力梯度作用或在邊壁附近潰滅時，空泡變形為扁平形或元寶形，最后分裂、潰滅，并在潰滅前瞬間，產(chǎn)生一束直徑幾微米的微型射流沖向壁面。隨后，研究者們陸續(xù)從理論、試驗驗證了此觀點，并對微射流進行了數(shù)學分析。Gregor?i?等[11]分別通過高速攝像機觀察到了近壁空泡潰滅時產(chǎn)生的沖向壁面的微射流。Plesset 等[12]通過數(shù)值計算，指出附壁空泡潰滅時微射流速度約為130 m/s，近壁空泡潰滅時速度約為170 m/s。Dular 等[13]研究了不同溫度下的空化腐蝕，結(jié)果表明：微射流速度為幾百米每秒，且應(yīng)用微射流方法能很好地預測空蝕趨勢。

在微射流對材料具體作用的研究方面，Tzanakis 等[14]認為潰滅微射流是蝕坑形成的主要原因，用干涉顯微鏡觀察了聲空化下鋼試件表面的點坑，并通過逆向工程及數(shù)理統(tǒng)計預測：近壁面聲空泡潰滅產(chǎn)生0.4～1 GPa 的沖擊壓力，且對應(yīng)的微射流速度為200～700 m/s。Reuter 等[15]通過納秒激光脈沖誘導產(chǎn)生單空泡和高速攝影記錄等手段，研究了近壁單空泡潰滅微射流對受微粒污染玻璃板的清洗作用，對應(yīng)不同的無量綱空泡壁面距離，提出了三種不同的清洗機理。Soyama 等[16]則利用空泡潰滅產(chǎn)生的沖擊作用處理硬鋁合金板表面，試驗結(jié)果表明，材料的疲勞壽命提高超過10 倍?？梢?，雖然上述空蝕分析、單空泡清洗以及空化表面強化各領(lǐng)域均以空泡潰滅微射流沖擊作用為主導，但均未對微射流沖擊進行詳細理論分析，而此理論分析對于深入揭示超聲空化機理、聯(lián)系空泡理論與作用材料行為以及進一步超聲空化工業(yè)應(yīng)用均具有重要意義。

為深入探究聲空泡潰滅微射流沖擊機理，首先估算空泡潰滅微射流速度，隨后依據(jù)流體力學、沖擊動力學建立率相關(guān)的微射流沖擊流固耦合三維數(shù)學模型，并用耦合歐拉拉格朗日方法分析，探究微射流沖擊形成凹坑的幾何尺寸及壁面壓強的變化，最后通過超聲空化試驗及反演分析驗證理論的可靠性并給出微射流沖擊強度、速度的合理預測。

1 微射流沖擊流固耦合數(shù)學模型

1.1 微射流速度

當傳入液體中的超聲波聲壓幅值達到空化閾值時，液體中產(chǎn)生空化效應(yīng)，在超聲波正負壓相的交替作用下，液體中的微小泡核歷經(jīng)生長、膨脹、壓縮、潰滅等一系列行為?？张菰诮瘫诿嫜莼瘯r發(fā)生不對稱潰滅，在潰滅歷程早期，空泡遠離壁面一側(cè)就產(chǎn)生指向壁面的微射流，空泡由初始的近似球形變?yōu)槭蓍L形再趨近于扁平形，最終微射流突破空泡近壁邊側(cè)沖向壁面，造成壁面空蝕損傷?？张轁邕^程極其復雜且充滿隨機性，潰滅微射流速度很難精確計算，文獻[17]得到了微射流速度估算公式：

式中：H 為空泡中心距壁面距離，R0為空化泡初始半徑，p∞為環(huán)境壓力，pV為泡內(nèi)飽和蒸氣壓，ρ2為液體密度。由于本文液體環(huán)境中附加了超聲，故p∞=p0+p2sin (ωt)，p0為大氣壓，p2為超聲聲壓幅值，ω 為超聲角頻率。故針對超聲空化這一物理現(xiàn)象，潰滅微射流的速度為：

從中可以看出，雖然空泡壁面距離越大，潰滅微射流速度越大，但微射流在沖向壁面過程中速度會衰減，到達壁面時速度反而不高，沖擊效果不如空泡壁面距離近時強烈，為便于分析，排除空泡壁面距離的影響，僅考慮H=R0的情況，故在本文具體研究環(huán)境中，微射流速度范圍為0～479 m/s。

1.2 數(shù)學模型

近壁空泡潰滅時產(chǎn)生一束沖向壁面的微射流，會造成材料表面空蝕損傷。微射流以較高速度沖擊壁面，在接觸壁面一瞬間速度突躍降低，同時產(chǎn)生激波傳入液體。激波波陣面后會發(fā)生密度、壓力突躍變化，稱之為受擾動區(qū)，此處液體壓縮性起主要作用，不可忽略；而波陣面前各物理量基本未變，稱之為未受擾動區(qū)，此處液體壓縮性可忽略。激波離體之后，液體壓力降低，側(cè)向射流形成。微射流沖擊壁面屬于強非線性流固耦合問題，整個過程強烈、復雜且短暫，最有研究意義的是微射流受擾動區(qū)的物理行為變化及材料壁面行為變化。首先依據(jù)流體力學，建立微射流沖擊三維數(shù)學模型。

以靜止壁面為參考系，微射流以初始速度v0沖擊壁面，產(chǎn)生的激波以速度c-v0傳入液體，波陣面前后的密度、壓力分別為ρ0、p0和ρ、p，微射流沖擊示意圖見圖1。

波陣面后壓強突躍增強，通常稱為水錘壓力，可達到幾百甚至上千兆帕，液體的黏性、表面張力及重力與水錘壓力相比極小，其影響可忽略，并且微射流沖擊強烈且短暫，可忽略外界因素影響。采用歐拉法描述液體行為變化，其控制方程為：

圖 1 微射流沖擊示意圖Fig. 1 Schematic diagram of micro-jet impact

上述兩式分別為液體的連續(xù)性方程和動量方程，v 為速度向量， ?為微分算符，且代表散度。上述方程組展開后有4 個方程、5 個未知量，方程組不閉合，需補充液體狀態(tài)方程，通常情況下液體密度僅是液體壓力的函數(shù)，如下式所示[18]：

式中：p0、ρ0分別為大氣壓及大氣壓下水的密度，B 是常數(shù)，可取為3 000，n 為無量綱常數(shù)，可取為7。由此，方程組閉合，理論上可得到唯一解，但上述式子本質(zhì)為流體N-S 方程，得到其一般精確解極其困難，因此，簡化上述方程。

以微射流為參考系，看作壁面以速度v0撞擊微射流，微射流初始速度、密度及壓力為0、ρ0和p0。沖擊過程中受擾動區(qū)液體的速度、密度和壓強可表示為液體初始值與相關(guān)攝動值的和，此時物理量均為空間坐標和時間的函數(shù)，如下式所示：

式(8)本質(zhì)為波動方程，沖擊較強烈時，聲速c 不再是一常數(shù)，有[19]：=?φ，

式中：k=7.15。此時式(8)成為非線性偏微分方程。此外，對式(6)中動量方程進行積分計算且引入速度勢函數(shù)，可得：

此時，微射流沖擊壁面液體區(qū)域的物理量已全部得出，控制方程如下：

在固體區(qū)域，使用拉格朗日法描述材料行為變化，設(shè)壁面質(zhì)點位移為U(ux, uy, uz)，其運動方程為：

式中：ρs為材料密度，σ 為質(zhì)點應(yīng)力張量。由于微射流為短時歷高強度沖擊載荷，材料變形過程中應(yīng)變率可高達104～106s-1量級，此時，材料應(yīng)變率效應(yīng)不能忽略，故選用典型的率相關(guān)本構(gòu)模型—J-C 模型[20]，此模型考慮了材料的應(yīng)變硬化效應(yīng)、應(yīng)變率強化效應(yīng)及熱軟化效應(yīng)，不過本研究中忽略溫度影響，則JC 方程可表示為：

式中：σe、σ0、εe分別為等效應(yīng)力、屈服強度及等效應(yīng)變， 為等效應(yīng)變率， 為參考應(yīng)變率，B、n1為應(yīng)變硬化參數(shù)，C 為應(yīng)變率強化參數(shù)。本研究壁面材料選用Al 1060，依據(jù)文獻[21]的有限元模擬結(jié)果，可得到Al 1060 的材料參數(shù)，見表1。

至此，液固區(qū)域的控制方程均已得到，加之初始條件及邊界條件，即可得到微射流沖擊流固耦合三維數(shù)學模型。

液體區(qū)域：

表 1 Al 1060 材料參數(shù)Table 1 Material parameters of Al1060

固體區(qū)域：

數(shù)學模型是一個偏微分方程組，得到其解析解非常困難，故借用有限元模擬的手段對其進行分析，利用耦合歐拉拉格朗日方法模擬分析微射流沖擊壁面的流固耦合特性。假設(shè)微射流為頭部半球形的圓柱，直徑為6 μm，初始速度分別為50、100、200、300、400、479 m/s，分析步長為1 ns。

2 結(jié)果與討論

2.1 凹坑幾何形貌分析

記微射流半徑為r0，r 為壁面上沿x 軸坐標，壁面無量綱距離即為r/r0。微射流沖擊下材料發(fā)生塑性變形，圖2 給出了不同沖擊速度下壁面變形深度h 的分布圖，結(jié)果表明：微射流沖擊下材料表面出現(xiàn)凹坑，凹坑兩側(cè)有材料塑性隆起，速度越大的微射流沖擊形成的凹坑深度越大；相比于凹坑深度，凹坑直徑受微射流速度影響較小。為進一步分析凹坑幾何形貌與微射流速度的關(guān)系，繪制凹坑直徑dp、深度hp、徑深比dp/hp以及凹坑直徑與微射流直徑之比dp/dj隨微射流速度的變化曲線，如圖3 所示?？芍S微射流速度增大，凹坑直徑和凹坑直徑與微射流直徑之比雖略有增大但變化不明顯，dp/dj處于0.95～1.15，見圖3(a)、(d)，故凹坑直徑與微射流速度無明顯相關(guān)性。由圖3(b)可知：微射流速度小于300 m/s 時，雖然凹坑深度隨微射流速度增大而緩慢增大，但其數(shù)值很小，形成的凹坑極不明顯，故可忽略；當微射流速度大于300 m/s 時，凹坑深度則隨微射流速度增大而明顯非線性增大，通過比較可知凹坑的直徑遠大于深度。由圖3(c)可見，凹坑徑深比與微射流速度有明顯的相關(guān)性，且凹坑徑深比隨微射流速度增大而減小，在形成明顯凹坑的前提下(即v0≥300 m/s 時)，凹坑的徑深比處于15～80 之間，這與其他學者的研究結(jié)論[22]相符。

圖 2 不同速度壁面變形深度分布曲線Fig. 2 Depth distribution of wall deformation at different velocities

圖 3 dp、 hp 、dp/hp、dp/dj 隨微射流速度變化曲線Fig. 3 Curves of dp , hp , dp/hp , dp/dj with different velocities

為進一步分析凹坑幾何形貌與微射流直徑的相關(guān)性，定義微射流初始速度為300 m/s，微射流直徑分別為4、6、8、10、12 μm，s 為壁面距離，從而進行仿真分析，不同微射流直徑下壁面變形深度分布如圖4 所示?？梢?，微射流直徑增大，形成的凹坑幾何尺寸也相應(yīng)增大。同樣，為對凹坑幾何尺寸進行詳細分析，提取圖4 中數(shù)據(jù)，繪制凹坑直徑dp、深度hp、徑深比dp/hp以及凹坑直徑與微射流直徑之比dp/dj隨微射流直徑的變化曲線，如圖5 所示。隨微射流直徑增大，凹坑直徑近似線性增大，而凹坑直徑與微射流直徑之比則變化不大，故凹坑直徑主要由微射流直徑?jīng)Q定，且dp/dj處于0.95～1.2，見圖5(a)、(d)。凹坑深度及凹坑徑深比隨微射流直徑增大而略有增加，見圖5(b)、(c)，但相比于微射流速度，凹坑徑深比受微射流直徑影響較小。

圖 4 不同微射流直徑壁面變形深度分布曲線Fig. 4 Depth distribution of wall deformation at differentmicro-jet diameters

圖 5 dp 、hp 、dp/hp 、dp/dj 隨微射流直徑變化曲線Fig. 5 Curves of dp , hp , dp/hp , dp/dj with different micro-jet diameters

由上可知，凹坑深度由微射流速度及微射流直徑共同決定；凹坑徑深比主要由微射流沖擊速度決定，且沖擊形成明顯凹坑的徑深比約為15～80；而凹坑直徑則主要由微射流直徑?jīng)Q定，且dp/dj≈0.95～1.2。

2.2 壁面壓強分析

圖 6 不同時刻壁面壓強分布曲線Fig. 6 Wall pressure distribution at different times

微射流沖擊下不同時刻壁面壓強分布，如圖6 所示。微射流沖擊壁面瞬間產(chǎn)生極高壓強，并在1～3 ns 內(nèi)達到最大，見圖6(a)、(b)，這是由于微射流撞擊壁面時速度急劇降低，液體積壓以致產(chǎn)生的激波波陣面后液體壓縮性不容忽視，從而密度壓強突躍增大，隨后壓力釋放，積壓液體以側(cè)向射流形式流出，壓力降低至一相對穩(wěn)定狀態(tài)，如圖6(c)、(d)、(e)、(f)。圖6 中可明顯看出壁面壓強隨微射流沖擊速度增大而增大，壁面壓強基本呈對稱分布且最大壓強出現(xiàn)在微射流沖擊邊緣，這是由微射流中反向傳播的球形激波首先由沖擊中心產(chǎn)生并逐漸擴散至兩邊而導致的。整理不同速度下壁面最大壓強，并由式(4)、(9)得到不同微射流速度下液體的最大聲速，見圖7，結(jié)果表明：在沖擊速度為479 m/s 時，壁面壓強可達682 MPa，此時對應(yīng)液體中聲速即激波速度可達2 435 m/s，這一方面體現(xiàn)了液體壓縮性對結(jié)果的重要性，另一方面也體現(xiàn)了空泡潰滅微射流沖擊的劇烈程度。隨后，為驗證數(shù)學模型及結(jié)果的準確性，進行超聲空化試驗并進行反演分析。

圖 7 壁面最大壓強及激波速度隨微射流速度變化曲線Fig. 7 Maximum wall pressure and shock wave velocityvaried with the increase of micro-jet velocity

3 試驗及反演分析

為驗證前述理論及結(jié)果的準確性，進行超聲諧振系統(tǒng)空化試驗，并針對材料出現(xiàn)的微小凹坑進行反演分析。為保證試驗的穩(wěn)定性及準確性，試驗材料選取厚約1 mm 的1 060 鋁薄板。超聲諧振系統(tǒng)由超聲波發(fā)生器、換能器、變幅桿和輻射桿組成，超聲波發(fā)生器會根據(jù)負載變化而自動跟蹤頻率，從而達到共振。本試驗中頻率范圍為20～20.6 kHz，工具頭端面振幅為8 μm，可通過調(diào)節(jié)振幅百分比控制空化強弱。除超聲諧振系統(tǒng)外，主要試驗儀器物品有VHX-600ESO 數(shù)碼顯微鏡、直徑為160 mm 的燒杯、酒精杯、1 060 鋁薄板等，放置好超聲諧振系統(tǒng)后，將鋁薄板固定在鋼塊上放入燒杯，設(shè)置好試驗參數(shù)后進行試驗，圖8 為試驗示意圖。為避免試驗后材料表面凹坑重疊，保證單凹坑由單束微射流沖擊造成，合理選取試驗參數(shù)：工具頭與材料表面距離3 mm，振幅為8 μm，試驗時間為1 min。試驗后將薄鋁板放入酒精杯中去除殘余水分，取出后待鋁薄板上酒精完全揮發(fā)，將其放置在VHX-600ESO 數(shù)碼顯微鏡下觀察，試樣試驗前后表面形貌圖如圖9 所示，挑選圖中有代表性的5 個凹坑并分別記錄其直徑及深度，見表2。

可見，試驗后試樣出現(xiàn)近似圓形的微小凹坑，即空蝕坑，其直徑大致在4～14 μm，對應(yīng)深度大致在0.06～0.88 μm，徑深比處于16～68，與前述理論結(jié)果基本吻合。空蝕坑幾何特征符合空泡潰滅微射流沖擊結(jié)果，可認為空蝕坑主要由微射流沖擊導致，這也和空蝕損傷主要由微射流引起的論述相符。由于空化微射流沖擊下材料形成的空蝕坑與球形壓痕試驗中基底材料的變形有很高的相似性，故應(yīng)用球形壓痕試驗理論，加之空蝕坑幾何特征的測量，可預測空蝕坑應(yīng)力及空化沖擊載荷，進而反演推出空化微射流沖擊強

圖 8 試驗示意圖Fig. 8 Test schematic diagram

表 2 凹坑幾何參數(shù)Table 2 Pit geometry parameters

圖 9 試驗前、后材料表面形貌圖Fig. 9 Surface topography of material before and after the test

度及速度。假設(shè)單空蝕坑由單束微射流沖擊導致空蝕坑為球冠幾何形狀，且其幾何特征可由直徑和深度完全描述，依據(jù)Tabor 壓痕試驗的開創(chuàng)性工作，可推導出空蝕坑球冠接觸邊緣的等效應(yīng)變εp為[23]：

式中：dp和hp分別為空蝕坑的直徑和深度。由于本試驗中凹坑直徑遠大于深度，故式(16)可變?yōu)椋?/p>

由式(17)可知，空蝕坑等效應(yīng)變僅與徑深比負相關(guān)。由前述J-C 方程可得空蝕坑等效應(yīng)力σp為：

空化沖擊載荷Lp為：

式中：ψ 為Francis[24]提出的約束因子，當材料為完全彈性行為時ψ=1.11，在彈塑性轉(zhuǎn)變過程中逐漸增加且當材料變?yōu)橥耆苄孕袨闀r達到最大值2.87。本研究中，假設(shè)空蝕坑發(fā)生了完全塑性變形，故約束因子取為2.87，Ap為空蝕坑投影面積，因空蝕坑直徑遠大于深度，故可寫為Ap=πdp2/4。因空化侵蝕主要由微射流沖擊導致，假設(shè)空蝕坑完全由微射流沖擊形成，即空化沖擊載荷全部由微射流提供，由前述結(jié)論可知凹坑直徑主要與微射流直徑相關(guān)，此處假設(shè)dp=dj，則可計算微射流沖擊強度σj：

依據(jù)水錘方程，即可得到微射流沖擊速度vj：

由上述分析可知，等效應(yīng)變、等效應(yīng)力、微射流沖擊強度及速度均與凹坑徑深比密切相關(guān)。依據(jù)試驗結(jié)果反演分析，可預測微射流沖擊強度及速度隨徑深比變化曲線，如圖10 所示。結(jié)果表明：由球形壓痕試驗理論反演預測的微射流沖擊強度及速度與前述理論結(jié)果匹配良好，證實了選用率相關(guān)J-C 模型的合理性，在16～68 的徑深比下，微射流沖擊強度為420～500 MPa，對應(yīng)的微射流速度為310～370 m/s。

本文所得到的空泡潰滅微射流的速度與相關(guān)文獻的結(jié)果相符。一些有關(guān)預測空泡潰滅微射流速度的例子在表3 中列出。Grant 等[25]記錄了水動力空泡內(nèi)爆產(chǎn)生的微射流的速度，結(jié)果為112～500 m/s。Chen 等[26]分別針對空泡距壁面10 μm 和1 μm 的情況，計算得到了對應(yīng)產(chǎn)生的微射流速度分別為192 m/s和755 m/s。Tzanakis 等[14]依據(jù)逆向工程理論，預測聲空泡潰滅產(chǎn)生的微射流的速度為200～700 m/s。Petkovsek 等[27]應(yīng)用激光誘導產(chǎn)生空泡，并觀測到微射流速度可達到780 m/s。此外，針對不同的試驗條件，Roberto 等[28]和Futakawa[29]分別提出微射流的速度為200～400 m/s 和250～300 m/s。通過與上述提到的相關(guān)文獻結(jié)果的對比，證實了本文的結(jié)果是合理的，選用的方法也是可靠的。

圖 10 微射流沖擊強度及速度隨凹坑徑深比變化曲線Fig. 10 Impact pressure and velocity of micro-jet with theincrease of dp/hp

表 3 空泡潰滅微射流速度預測的案例Table 3 Examples of bubble collapse micro-jet velocityprediction

4 結(jié) 論

本文利用流體力學及沖擊動力學建立了近壁聲空泡潰滅微射流沖擊的率相關(guān)三維流固耦合模型，對其沖擊特性進行了分析，并進行了超聲空化試驗驗證及反演分析。結(jié)果表明：近壁聲空泡潰滅產(chǎn)生的微射流速度最大可達479 m/s；微射流沖擊下材料表面會出現(xiàn)微型凹坑，凹坑深度由微射流速度和微射流直徑共同決定且隨其增大而增大，凹坑直徑主要由微射流直徑?jīng)Q定且dp/dj≈0.95～1.2，而凹坑徑深比則主要由微射流速度決定且其比值為15～80；壁面壓強隨微射流速度增大而增大，基本呈對稱分布且最大壓強出現(xiàn)在微射流沖擊邊緣，壓強增大可導致液體中激波速度增大，在沖擊速度為479 m/s 時，壁面壓強及激波速度最大分別可達682 MPa 及2 435 m/s；超聲空化試驗驗證了微射流沖擊下材料表面出現(xiàn)的微小凹坑，基于球形壓痕理論的反演分析法表明凹坑等效應(yīng)變、等效應(yīng)力、微射流沖擊強度及速度均與凹坑徑深比密切相關(guān)，在16～68 的徑深比下，微射流沖擊強度為420～500 MPa，對應(yīng)的微射流速度為310～370 m/s；反演分析結(jié)果與理論分析結(jié)果相符，驗證了選用J-C 材料模型及反演分析方法的合理性及準確性，也為后續(xù)工程應(yīng)用中控制空化強度、微射流速度等提供了理論參考。