曲面固壁附近空化泡潰滅動力學的流體體積數值研究*

單鳴雷 楊 雨 邱順俐 束方勇 韓慶邦

(1 河海大學常州市傳感網與環境感知重點實驗室并江蘇省輸配電裝備技術重點實驗室 常州 213022)

(2 江蘇省“世界水谷”與水生態文明協同創新中心 南京 211100)

0 引言

空化是指液體中局部壓強降低，其內部發生汽化并形成微小空穴(空化泡)的過程。流體中的空化泡在遇到周圍液體壓強增大的情況下，體積將急劇縮小或潰滅。空化泡潰滅時間很短(微秒級)，潰滅時會在局部產生極高的瞬時壓強。當潰滅發生在固體表面附近時，流體中不斷潰滅的空化泡所產生的極高壓強反復作用于固體表面，對固壁面造成破壞，這種現象稱為空蝕，也稱氣蝕[1-2]。一方面，空化效應會產生如水輪機的葉片損傷[3]、高速軸承表面材料呈點狀脫落及螺旋槳的空蝕噪聲[4-5]等不利影響。另一方面，人們發現合理利用空化效應會給生產生活帶來極大的便利，如空化清洗技術用于裂縫表面的清潔[6-8]、使用空化泡作為藥物的運輸工具從而將藥物送達病變組織表面[9-10]。因此，無論是避免還是利用空化效應，都有必要對固壁附近空化泡潰滅過程進行深入研究。

目前空化泡在固壁附近潰滅的動力學研究已相當豐富，模型結構也從原始的平直固壁附近空化泡潰滅模型進一步發展為任意幾何形狀固壁附近空化泡潰滅模型。Tomita 等[11]從理論和實驗的角度給出了曲面邊界影響下的氣泡動力學特性，研究表明，如果邊界稍微凹陷，最顯著的遷移發生在第一次氣泡潰滅的時候。液體射流沖擊氣泡表面遠側的速度隨著邊界曲率的增加而增加。Ma等[12]利用霍普金森桿(the Hopkinson bar,HPB)測量技術，對不同曲率邊界和不同爆炸距離下的電火花氣泡動力學和載荷特性進行了一系列實驗研究，并通過實驗結果總結了一些物理量的變化規律。Brujan等[13]對超聲誘導產生的空化泡進行了實驗研究，比較了兩種尺寸的空化泡近壁潰滅時的射流速度與沖擊壓力的大小，發現尺寸較小的空化泡其射流速度與沖擊壓力都要小于大尺寸空化泡。Li 等[14]利用流體體積(Volume of fluid,VOF)法研究了空化泡在錐形固壁附近潰滅的過程，提出了錐形固壁角與空化泡潰滅形態、潰滅時間以及液體射流的關聯，計算了近錐形固壁空化泡潰滅時間的理論值，并通過數值模擬驗證了此理論值的有效性。魏夢舉等[15]基于合理的簡化模型，重點分析了空泡潰滅過程對球體的影響，研究表明，百微米級氣泡能夠通過射流方式驅動同尺度固體顆粒向前顯著運動，而氣泡受限的程度對射流驅動的效果具有重要影響，顆粒半徑越小、氣泡與顆粒間距越近，空泡潰滅時對顆粒的推動效果越明顯。

以上研究都是關于非平直壁面附近或其他復雜邊界附近空化泡潰滅動力學問題的相關研究工作，這些工作并沒有涉及空化泡潰滅過程中的熱力學特性的研究。然而在空化泡潰滅的最后階段，熱效應發揮了極其重要的作用，隨著空化泡的潰滅收縮，泡內溫度急劇上升，潰滅瞬間溫度擴散到周圍流體中，對泡外流體或固壁面都將產生重要影響。為更全面地了解空化泡潰滅過程，近年來，專家學者們對空化泡潰滅階段的熱力學問題開展了一系列實驗。Vokurka 等[16]對電火花誘導生成的空化泡熱力學特性進行了實驗研究，空化泡表面的溫度通過光學傳感器和聲傳感器進行測量，實驗結果顯示，空化泡表面溫度可以達到5900 K。但實驗存在過多的不可控因素，且空化泡潰滅過程較為復雜，其溫度測量存在較大的難度。因此，采用數值方法[17-21]研究空化泡熱力學問題顯得十分必要。Yang等[22]采用格子玻爾茲曼方法研究了空化泡潰滅的熱力學效應，描述了氣液相的溫度演化規律。Liu 等[23-24]使用VOF 模型討論了單壁與平行雙壁附近空化泡潰滅對固壁表面的熱影響，研究結果表明，在空化泡被射流穿透時，射流速度達到最大值，由于射流沖擊，壁面溫度明顯降低。但以上研究的物理模型均為平直固壁，未探討其他固壁模型對空化泡潰滅的影響。

本文采用VOF 數值模擬方法研究曲面固壁附近空化泡潰滅的過程，并對此模型下空化泡形態演化進行實驗觀測，驗證數值模擬方法的有效性。首先，分析了空化泡近壁潰滅過程中相場、壓力場、溫度場的演化，討論了空化泡近壁潰滅階段的熱動力學作用機理。其次，研究了曲面固壁尺寸、位置參數對空化泡潰滅時間、射流速度的影響。最后，詳細探討了不同條件下，曲面固壁溫度的變化情況，深入分析了高溫損傷對固壁表面的作用。

1 數值方法

數值計算固壁附近空化泡運動過程時，假設空化泡內為水蒸氣，滿足理想氣體狀態方程。整個空化泡潰滅過程不考慮重力因素的影響[25]，忽略氣體與液體之間的質量傳遞。存在相界面的流動中，直接求解流場中的物理量會使其在界面處發生間斷，導致計算的發散或計算結果的虛假，為保證界面處所要求解物理量值的合理性，對這類流動進行數值研究時，需要進行足夠精確的界面捕捉。VOF方法可以較好地保持流體質量守恒并且可以處理界面拓撲結構的變化[26]，因此本文采用VOF 方法來計算空化泡運動過程。

1.1 控制方程

空化泡流場中的流動方程由連續性方程、動量方程及能量方程組成，表達式分別為[27]

其中，ρ表示流體的密度，μ表示流體的黏度，U表示流場速度，p表示流場壓力，σ表示表面張力系數，κ表示表面曲率，n表示兩相流界面中的單位法向量，T表示溫度，E表示總能量，k表示流體的傳熱系數。

VOF 方法通過一個標量的輸運方程來描述各相的體積分數，表達式為[28]

其中，α表示構成相的體積分數大小，當α= 0 代表該網格處于空化泡內部，當0＜α ＜1 代表該網格處于空化泡外即液體內部，當α= 1 代表該網格處于兩相交界面上。

各相的密度以及黏度可分別表示為

其中，α1+α2= 1，α1代表液體的體積分數，α2代表氣體的體積分數。

1.2 計算模型

本文采用VOF 方法研究空化泡鄰近曲面固壁潰滅的過程，幾何模型如圖1 所示。曲面固壁為球形，半徑設為r，球形空化泡的初始半徑設為R0。為方便描述空化泡與固壁間的距離，引入無量綱位置參數λ=d/R0，其中d為空化泡中心到曲面固壁右側一點的距離，在數值模擬過程中，改變曲面固壁的半徑時始終保持曲面固壁右側一點固定。為減少壓力邊界條件對潰滅過程的影響，將計算域大小設為30R0×30R0。計算域的四周均為壓力出口邊界，曲面固壁表面為無滑移固壁邊界條件。

圖1 計算模型Fig.1 Computational model

在數值模擬軟件中導入網格文件，由于Gambit 中構建模型結構時默認的長度單位為m，因此要根據實際目標需求對網格的尺寸進行重新定義并檢查網格質量。采用基于壓力的瞬態求解器，選用多相流模型中的VOF 模型進行計算，開啟能量方程。壓力-速度耦合方式運用PISO 算法，空間離散方法中，梯度項使用Least Squares Cell Based 格式，壓力項使用PRESTO!格式，體積分數項使用Geo-Reconstruct格式，密度項、能量項、動量項均采用Second Order Upwind 格式[29]。計算的時間步長設置為10-8s，開啟殘差曲線監控即可開始迭代計算，觀察計算結果是否收斂。

1.3 網格無關性驗證

不同的網格數量會產生不同的計算結果，為使計算結果誤差較小，需要選取恰當的網格數。在比較不同網格數的計算結果時，采用如下初始計算條件：曲面固壁半徑r= 0.2 mm，空化泡半徑R0=0.2 mm，無量綱位置參數λ=1.1，液體中壓強pl=3×106Pa，空化泡內壓強pv=1.01×105Pa。

圖2 為空化泡橫軸上射流速度的變化情況，其中橫軸表示空化泡軸上的位置坐標，縱軸表示軸上的射流速度。觀察發現射流速度在某一位置達到最大值，且不同網格數下的空化泡軸上射流速度差異較小，幾乎可以忽略不計。綜合考慮計算時間、計算資源等的限制，最終選取網格數為200000的工況進行數值模擬。

圖2 網格無關性驗證Fig.2 Verification of the grid independence

1.4 實驗驗證

為驗證數值模擬結果的有效性，本章節將曲面固壁附近空化泡演化的實驗研究與數值仿真結果進行對比，如圖3 所示。圖3(a)為高壓放電致泡實驗下空化泡近曲面固壁潰滅的形態變化圖，圖3(b)為數值模擬仿真得到的空化泡近曲面固壁潰滅的相圖。

圖3 空化泡形態演化實驗結果圖和仿真結果圖Fig.3 Comparison of cavitation bubble profiles between the experimental results and the simulation results

將實驗結果與仿真結果進行對比，由于實驗研究和數值仿真中空化泡尺寸的差異，其潰滅時間不能直接比較，需要進行歸一化處理。比較兩圖可以發現空化泡的形態演化過程一致，且歸一化后得到的無量綱時間相同。空化泡形態演化均為從初始的球狀逐漸被壓縮成類似橢球形狀，可以明顯地看出空化泡上下泡壁收縮速度大于左右泡壁收縮速度，左側泡壁收縮的速度與右側泡壁收縮的速度相比較慢。因此，在t*=0.5789時刻，左側泡壁出現尖端現象，隨后繼續潰滅。在這過程中空化泡沒有出現類似于近平直固壁潰滅時產生的凹陷或者被射流穿透現象，而是維持一個整體潰滅直至空化泡表面完全消失。由于空化泡潰滅形態演變的實驗結果與數值仿真結果的吻合性較好，故進一步驗證了數值模擬曲面固壁附近空化泡潰滅過程的正確性，為接下來對數值模擬結果的分析奠定基礎。

2 結果分析

2.1 壓力場、速度場、溫度場演化分析

為了進一步探究空化泡鄰近曲面固壁潰滅機理，對曲面固壁附近空化泡潰滅過程中壓力場分布情況進行了分析。圖4為曲面固壁半徑r=0.3 mm，空化泡初始半徑R0=0.2 mm，位置參數λ=1.4的工況下的計算域壓力場演化圖。

圖4(a)～圖4(d)為球形空化泡形變為長空化泡的壓力場演化圖。由于固壁的阻滯效應，空化泡縱向方向受到的擠壓大于橫向方向受到的擠壓，空化泡呈現出被左右拉長的狀態。受到液體回彈效應的影響，空化泡右側形成錐形高壓區域。圖4(e)～圖4(g)顯示了空化泡潰滅末期壓力場演化的過程，空化泡右側的高壓使得泡壁表面出現輕微凹陷現象，但是由于垂直于曲面固壁右側的沖擊壓力沒有達到能夠擊穿泡壁左側表面的大小，因此空化泡壁沒有被完全穿透，而是維持此狀態繼續潰滅直至泡壁表面完全消失，最后空化泡完全消失的壓力場圖如圖4(g)所示。圖4(h)和圖4(i)為空化泡潰滅后一段時間內的壓力場分布圖，通過觀察可以發現曲面固壁右側表面形成圓形高壓區，這是造成固壁損傷不可忽略的原因。

圖4 曲面固壁附近空化泡潰滅階段壓力場演化Fig.4 Pressure field evolution of cavitation bubble collapse near curved wall

對曲面固壁半徑r= 0.3 mm，空化泡初始半徑R0= 0.2 mm，位置參數λ= 1.4 的工況進行空化泡潰滅時的速度場演化分析，如圖5 所示，其中t=6.89 μs 時刻空化泡完全潰滅。計算域中速度最高的位置用紅色箭頭表示，速度最低的位置用藍色箭頭表示，箭頭的方向即表示流場速度的方向。

觀察圖5(a)和圖5(b)可以發現在空化泡潰滅初期，空化泡上下泡壁收縮的速度較快，而由于曲面固壁的阻滯作用使得靠近曲面固壁一側的空化泡壁收縮的速度較慢，空化泡右側泡壁收縮速度明顯高于空化泡左側收縮速度。圖5(c)和圖5(d)顯示了空化泡潰滅末期流場速度的分布情況，圖中指向曲面固壁的速度值較大，隨著空化泡的潰滅，會產生對空化泡附近曲面固壁的射流沖擊，對曲面固壁造成損傷。圖5(e)和圖5(f)為空化泡完全潰滅之后流場速度分布圖，圖中指向曲面固壁的速度仍然存在，然而隨著時間的推移，指向曲面固壁的速度值越來越小，并且壁面附近的速度方向發生改變，沿著空化泡軸線處形成上下對稱的渦流。

圖5 曲面固壁附近空化泡潰滅階段速度場分布Fig.5 Velocity field evolution of cavitation bubble collapse near curved wall

空化泡潰滅時會產生高溫高壓效應，其會帶來空化泡周圍溫度的變化。圖6 為空化泡潰滅過程中溫度場的演化情況，其中圖6(a)～圖6(f)為空化泡潰滅階段溫度場的演化過程，圖6(g)和圖6(h)為空化泡完全潰滅后一段時間溫度場的演化圖。初始時刻，曲面固壁半徑r= 0.3 mm，空化泡半徑R0= 0.2 mm，位置參數λ= 1.4，空化泡內溫度與泡外溫度均為300 K。

從圖6(a)～圖6(c)中可以看出空化泡潰滅初期，泡內外溫度差異較小，泡內溫度上升緩慢。隨著空化泡被壓縮，泡內溫度逐漸升高，泡內外溫度差異逐漸增大，空化泡潰滅后期，泡內溫度上升速度與空化泡潰滅初期相比較快，如圖6(d)到圖6(f)所示。然而在此過程中，泡外溫度值變化較小，與泡內溫度的變化形成鮮明對比。這是由于空化泡體積小，周圍流體計算域范圍大，當空化泡收縮產生高溫效應時，周圍流體對其冷卻速度較快，導致空化泡周圍流體總體溫度上升較少。觀察圖6(g)和圖6(h)可以發現空化泡壁完全消失后，隨著流場中流體的運動，空化泡潰滅產生的高溫區域將向曲面固壁一側偏移。這將導致曲面固壁受到空化泡潰滅的影響，從而對曲面固壁表面的材料造成高溫損傷或其他不可忽略的影響。

圖6 曲面固壁附近空化泡潰滅階段溫度場圖Fig.6 Temperature field evolution of cavitation bubble collapse near curved wall

2.2 潰滅時間分析

為了進一步考察曲面固壁下空化泡潰滅時間的變化情況，畫出了空化泡初始半徑R0=0.2 mm，曲面固壁半徑r分別為0.1 mm、0.2 mm、0.3 mm時，空化泡潰滅時間隨位置參數變化的折線圖，如圖7所示。圖中橫軸表示位置參數，縱軸表示潰滅時間，3 條折線分別表示不同曲面固壁半徑與空化泡初始半徑的比值下的空化泡潰滅時間。

從圖7 中可以看出，曲面固壁半徑與空化泡初始半徑的比值相同時，空化泡潰滅時間隨著位置參數的增大而減小，空化泡與固壁間的距離越遠，潰滅速度越快。保持位置參數不變，曲面固壁半徑與空化泡初始半徑的比值越大，曲面固壁對空化泡潰滅過程的阻滯效應越強，空化泡潰滅時間越長。隨著位置參數的增大，空化泡逐漸遠離固壁面，大尺寸固壁下空化泡潰滅時間下降幅度較大，不同曲面固壁半徑與空化泡初始半徑的比值所帶來的空化泡潰滅時間差異在逐漸減小。

圖7 不同曲面固壁半徑下潰滅時間隨位置參數的變化Fig.7 The change of collapse time with position parameter under different curved wall radius

2.3 射流速度分析

空化泡近壁潰滅產生的射流沖擊會對材料表面造成不可忽略的影響，選取空化泡初始半徑R0= 0.2 mm，不同曲面固壁半徑r，位置參數λ分別為1.1、1.4、1.7、2.0的工況進行空化泡軸上射流速度的討論，如圖8 所示。圖中橫軸表示空化泡軸線上的位置，縱軸表示空泡軸線上的射流速度。

圖8 不同曲面固壁半徑下射流速度隨位置參數的變化Fig.8 The change of jet velocity with position parameter under different curved wall radius

從圖8 中可以發現，射流速度大致變化情況為隨著橫軸值的增大，空化泡軸上射流速度逐漸增大再逐漸減小，在軸上某一位置會迎來射流速度的峰值，且相同位置參數下，射流速度峰值的位置隨著曲面固壁半徑與空化泡初始半徑比值的增大而減小。出現以上現象的原因為，曲面固壁半徑越大，空化泡潰滅點與固壁間的距離越小，射流速度峰值的位置與固壁間的距離也越小。隨著位置參數的增大，空化泡與曲面固壁間的距離越遠，固壁對空化泡潰滅過程中產生的液體射流的影響越小，空化泡軸上射流速度的峰值越小。另外，從λ= 1.4、λ= 1.7、λ= 2.0 三幅圖中可以看出，同一位置參數下，曲面固壁半徑與空化泡初始半徑的比值越大，固壁對空化泡潰滅產生的阻滯效應越強，空化泡潰滅階段其軸上射流速度的峰值越大。

2.4 固壁溫度分析

本部分主要對空化泡潰滅過程中熱力學效應進行討論，由于實驗研究中測量溫度變化情況較為困難，且測得數據受較多因素的影響，所得結果與實際結果相比誤差較大。而數值模擬可以便捷地調整參數大小，且數值仿真成本較低，故采用數值模擬的方式考察空化泡潰滅的熱力學特性。根據仿真結果，重點研究了不同條件對壁面溫度的影響。由于模型為軸對稱結構，觀察了曲面固壁上表面3 個位置點的溫度隨位置參數的變化情況來探究曲面固壁上的溫度變化情況，其中A、B、C三個位置點分布如圖9所示。

為探究不同位置參數下曲面固壁上的溫度分布情況，畫出了固定曲面固壁半徑與空化泡初始半徑的比值為1，位置參數λ分別為1.1、1.4、1.7、2.0時，曲面固壁上A、B、C三點的溫度隨時間變化的情況，如圖10 所示。其中A、B、C三點在曲面固壁上的具體位置分布如圖9 所示，圖10 中橫軸表示空化泡潰滅的時間，縱軸表示固壁上的溫度值。

圖9 A、B、C 點位置分布示意圖Fig.9 Schematic diagram of location distribution of points A、B and C

圖10 不同位置參數下曲面固壁溫度分布Fig.10 Temperature distribution of curved wall under different position parameters

從圖10中可以看出，固壁溫度在空化泡潰滅前期基本保持不變，空化泡潰滅末期，固壁溫度急劇上升，在潰滅時刻達到溫度最大值，隨后固壁溫度迅速下降。這是因為潰滅前期流場溫度變化較小，對固壁溫度產生的影響幾乎可以忽略不計。隨著空化泡潰滅時刻的到來，空化泡內的高溫在潰滅瞬間向外迅速擴散，且潰滅瞬間對壁面產生高壓沖擊，會對空化泡附近的固壁面造成較大的影響。固壁面在短暫的高溫沖擊之后，被周圍流體迅速冷卻，使溫度急劇下降。當位置參數相同時，曲面固壁上距離空化泡越近的點，其溫度峰值越大，空化泡的潰滅對其溫度造成的影響越大。比較不同位置參數下壁面溫度變化可得，位置參數越大，曲面固壁上溫度峰值越小，即空化泡與曲面固壁間的距離越遠，空化泡潰滅時刻到達固壁上的溫度峰值越小，空化泡潰滅對固壁造成的高溫損傷越小。

由于曲面固壁上C點受空化泡潰滅的影響最大，本文討論了相同位置參數下曲面固壁右側C點的溫度變化情況，如圖11 所示。圖中橫軸表示空化泡潰滅時間，縱軸表示固壁上的溫度值。

圖11 不同曲面固壁半徑及泡內外壓差下曲面固壁溫度的變化Fig.11 The change of the temperature of curved wall under different curved wall radius and pressure difference between inside and outside the bubble

由圖11可知，曲面固壁溫度峰值隨著曲面固壁半徑與空化泡初始半徑比值的增大而增大，且曲面固壁半徑與空化泡初始半徑的比值越大，固壁溫度到達峰值所需的時間越長。固壁與空化泡間的距離相同時，曲面固壁的尺寸越大，空化泡軸上的射流速度越快，其溫度到達固壁面時越高。而曲面固壁半徑越大，空化泡潰滅時間越長，熱量到達固壁所需的時間越長。泡內壓強越大，泡內外壓差越小，空化泡潰滅的時間越長，導致曲面固壁上溫度到達最大值所需的時間越久，壁面溫度峰值越小。

3 結論

本文采用VOF 方法對空化泡鄰近曲面固壁潰滅的熱動力學行為進行研究，將曲面固壁附近高壓放電致泡實驗結果與數值仿真結果進行對比，從壓力場、速度場、溫度場、空化泡潰滅時間、空化泡軸上射流速度等方面分析不同物理參數對空化泡潰滅過程的影響。實驗所得結果與數值模擬所得空化泡潰滅形態變化相同，驗證了數值模擬空化泡鄰近曲面固壁潰滅過程的正確性。通過改變曲面固壁尺寸與空化泡初始半徑的比值及空化泡與曲面固壁間的距離可以發現，空化泡潰滅時間隨著位置參數的增大而減小，曲面固壁半徑較小時可以加速空化泡的潰滅。空化泡軸上射流速度的峰值隨位置參數的增大而減小，而曲面固壁半徑與空化泡初始半徑的比值越大，射流速度的峰值越大。

空化泡潰滅時產生的高溫在壓力和流體速度的影響下，熱量沿著射流方向轉移到壁面上。位置參數越大，壁面溫度越小，壁面上距離空化泡越近的點，空化泡潰滅時對其溫度造成的影響越大。曲面固壁的尺寸越大，固壁上溫度峰值越大，到達溫度峰值所需的時間越長。固壁上的溫度峰值隨著泡內壓強的增大而降低，泡內壓強增大，空化泡潰滅的時間變長，隨之帶來曲面固壁溫度到達峰值所需的時間變長。

本文對空化泡近曲面固壁潰滅的熱動力學研究，有利于人們深入了解空化泡近壁潰滅機理，并對工程中不同形狀固壁附近空化效應的認知提供一些幫助。同時揭示了空化泡近壁潰滅的高溫影響，對工程實際中空化在材料表面的應用或規避具有重要作用，為未來研究復雜幾何固壁附近空化效應提供參考。