水下爆炸兩氣泡相互作用的數(shù)值計算研究

李健，榮吉利，林賢坤，項大林

(1.廣西科技大學(xué) 汽車與交通學(xué)院，廣西 柳州545006;2.北京理工大學(xué) 宇航學(xué)院，北京100081;3.廣西汽車零部件與整車技術(shù)重點實驗室，廣西 柳州545006)

0 引言

兩氣泡或多氣泡融合現(xiàn)象在物理、化學(xué)、醫(yī)學(xué)及工業(yè)生產(chǎn)中廣泛存在并占據(jù)著重要的地位，如化學(xué)領(lǐng)域中的氣體與液體、液體與液體之間的相互反應(yīng);醫(yī)學(xué)領(lǐng)域中利用多氣泡擊碎結(jié)石現(xiàn)象;能源生產(chǎn)中利用沸騰使水變成蒸汽從而推動渦輪轉(zhuǎn)動的過程等。整個氣泡融合過程主要包括氣泡間相互吸引而引起的運動、氣泡壁間液體的擠壓、失穩(wěn)、射流形成及最終的融合等階段。影響氣泡的融合現(xiàn)象涉及的因素較多，主要有質(zhì)量交換、氣泡表面張力、范德華力、布朗運動等。致使氣泡融合機(jī)理還沒有被完全揭示，很多行為特點還難以預(yù)測，因此，研究兩氣泡或多氣泡之間的相互作用能加深人們對這些領(lǐng)域中相關(guān)物理現(xiàn)象的理解，對許多工程應(yīng)用也有著重要的意義，有必要對氣泡運動、水射流、氣泡融合問題進(jìn)行深入的研究。

氣泡融合問題可通過實驗與仿真的方法進(jìn)行研究，實驗方法主要是利用氣泡發(fā)生器在液體中產(chǎn)生大量氣泡，通過相應(yīng)的測試設(shè)備研究氣泡運動與融合的行為規(guī)律［1］，此外也可利用水下爆炸這一簡單而有效方法進(jìn)行氣泡融合力學(xué)特性的研究。數(shù)值仿真計算方面目前大都采用邊界元或有限體積法［2-5］。多氣泡數(shù)值計算方面，Chew 等［6］利用邊界積分方法研究了近剛性壁面兩氣泡的運動特性，結(jié)合實驗測試驗證方法的正確性。Rungsiyaphornrat等［7］利用邊界積分方法研究了水下爆炸相互融合問題。李章銳等［8］采用邊界積分方法對水下三維氣泡的動力學(xué)特性進(jìn)行模擬，研究了浮力對氣泡運動特性與射流形成之間的聯(lián)系。王詩平等［9］通過自行研制的水中氣泡發(fā)生裝置，通過實驗測試方法了解了兩氣泡相互作用過程中出現(xiàn)的融合、射流等現(xiàn)象，通過實驗提出了氣泡之間的無量綱距離、無量綱周期等參數(shù)來描述氣泡耦合特性。

邊界元方法只在邊界離散，大大降低了計算成本，但是在射流產(chǎn)生、氣泡破碎及環(huán)形氣泡形成方面有一定的局限性，此方面目前研究成果較少，而采用有限元方法進(jìn)行氣泡特性的研究相對更少。為了對兩氣泡運動數(shù)值模擬研究提供有益的補充，本文對單個氣泡在水下運動特性進(jìn)行仿真分析，通過將氣泡最大半徑、脈動周期與經(jīng)驗公式的對比，驗證有限元模型與計算方法的正確性，以此為基礎(chǔ)，研究兩氣泡水平距離、垂直距離、初始角度等參數(shù)對氣泡運動特性的影響規(guī)律。

1 理論分析

1.1 氣泡初始條件

目前關(guān)于氣泡初始狀態(tài)求解的方法較多，主要包括Rayleigh 模型［10-11］和文獻(xiàn)［12 -14］提出的氣泡體積加速度模型，后者主要是分別建立了藥包參數(shù)和流場壓力以及氣泡體積加速度與流場壓力之間的關(guān)系，并利用壓力等效關(guān)系最終確定氣泡的初始半徑與初始膨脹速度，其體積與體積膨脹速度的表達(dá)式如(1)式～(2)式所示:

式中:V為氣泡體積;ρ為流體密度;Pc=;mc、ac分別為藥包的藥量和半徑;K1、K2、A、B 為與炸藥有關(guān)的材料參數(shù)。對于密度1 600 kg/m3的TNT 炸藥，K1=52.1，K2=0.090，A = 0.18，B =0.185. 考慮r0=即可求出不同時刻氣泡的初始半徑及徑向速度。

如圖1 所示為1.0 kg 球形TNT 炸藥在0 ～1.2 ms時間段內(nèi)，氣泡半徑及徑向速度的時程曲線。

圖1 氣泡半徑及徑向速度時間歷程曲線Fig.1 Time-history curves of radius and velocity of bubble

炸藥引爆后，壓力瞬間增加，故初始時刻氣泡徑向速度變化率非常明顯，之后的沖擊波傳播過程中氣泡膨脹速度逐漸平緩。針對這一規(guī)律，由圖1 可見，當(dāng)t≥0.8 ms 時，氣泡徑向速度的變化已經(jīng)很小了，即可認(rèn)為此時初始?xì)馀菪纬伞１疚模瑃 取0.8 ms時對應(yīng)的氣泡膨脹速度為174.5 m/s，初始半徑為0.303 m.

1.2 氣泡運動的經(jīng)驗公式

研究人員通過大量的實驗獲取足夠的數(shù)據(jù)，結(jié)合理論推導(dǎo)、曲線擬合等得到了水下爆炸中沖擊波壓力、氣泡壓力、氣泡半徑及脈動周期等經(jīng)驗公式，為計算結(jié)果的準(zhǔn)確性提供有力的參考。氣泡脈動最大半徑與周期的經(jīng)驗公式分別如(3)式和(4)式［14］所示:

式中:K 為與炸藥性質(zhì)有關(guān)的系數(shù)，對TNT 炸藥，K取1.6;p0為周圍水的靜壓力(kg/cm2);W 為氣泡的裝藥量(kg);z0為氣泡所處位置的流體靜壓的等效水深(m)。

2 數(shù)值模擬

2.1 有限元建模

本節(jié)首先建立流場區(qū)域的三維有限元模型，流場區(qū)域為15 m×15 m ×15 m，流場全部用六面體歐拉網(wǎng)格劃分，單元數(shù)量為135×135×135=2 460 375 個，其中每個單元的幾何尺寸相同，且每個方向上的長度均為0.111 m，約為氣泡初始半徑的1/3，能夠滿足計算精度的要求［15］。通過對MSC.DYTRAN 有限元軟件的二次開發(fā)，實現(xiàn)流場初始與邊界條件的定義，并將數(shù)值計算結(jié)果與經(jīng)驗公式對比，驗證有限元模型與數(shù)值計算方法的正確性。

2.2 初始條件及邊界條件定義

為在有限區(qū)域內(nèi)真實反映流場特性，需定義流場的靜水壓力以及與之對應(yīng)的邊界條件。由于MSC. DYTRAN 有限元軟件并不能通過直接定義重力加速度實現(xiàn)流場靜水壓力的定義。因此，需根據(jù)每個網(wǎng)格所處的位置，通過EXINIT 子程序?qū)㈧o水壓力值逐一賦到對應(yīng)的網(wǎng)格上。同樣，為了能夠在有限的計算區(qū)域真實模擬無限水域流場特點，即當(dāng)流場邊界兩側(cè)出現(xiàn)壓力差時，流體會從高壓區(qū)通過邊界流向低壓區(qū)，可利用EXFLOW2 子程序?qū)⒉煌撵o水壓力值賦到流場邊界上，最終實現(xiàn)三維流場靜水壓力與邊界條件的定義。

2.3 狀態(tài)方程

文中忽略氣泡運動對氣體壓力的影響，認(rèn)為氣體的壓力僅和氣泡的初始狀態(tài)及其體積有關(guān)，即氣泡內(nèi)的壓力p 與氣泡體積V 的關(guān)系:

式中:pcon為可冷凝氣體的飽和蒸汽壓，量級與大氣壓相當(dāng)，一般可忽略;p0和V0分別為氣泡形成時的初始壓力和初始體積;γ 為氣體的比熱比，對于TNT 水中爆炸的爆轟產(chǎn)物，γ 取1.25.

水的狀態(tài)方程利用MSC.DYTRAN 軟件中已有的多項式狀態(tài)方程來描述，其形式為

式中:μ=(ρ/ρ0)-1;ρ0為材料參考密度;e 為比內(nèi)能;對于水介質(zhì)，取a1=2.134 ×109Pa，a2=6.561 ×109Pa，a3=1.126 ×1010Pa.

3 計算結(jié)果及討論

3.1 單個氣泡數(shù)值計算結(jié)果

對于單一氣泡，計算工況為裝藥量為1.0 kg 的TNT 炸藥在水下10 ～45 m 處爆炸，并將氣泡脈動最大半徑與周期與經(jīng)驗公式(3)式、(4)式進(jìn)行對比，驗證數(shù)值計算的正確性。如圖2 所示，不同工況下氣泡脈動半徑的時間歷程曲線。

圖2 氣泡脈動半徑時間歷程曲線Fig.2 Time-history curves of radius of bubbles detonated at different depth

由圖2 可看出，氣泡脈動的初始半徑為0.303 m，隨著爆心處水深的增加，氣泡最大半徑及與之對應(yīng)的脈動周期也會逐漸減小。

如表1 所示，1 kg 球形TNT 炸藥在10 ～45 m 等不同水深處運動，其最大半徑與脈動周期與經(jīng)驗公式的對比。可看到，氣泡脈動最大半徑的仿真結(jié)果略大于經(jīng)驗公式，而脈動周期略小于經(jīng)驗公式結(jié)果，但從整體上而言，數(shù)值計算結(jié)果與經(jīng)驗公式吻合得較好，氣泡最大半徑的誤差能夠控制在8%以內(nèi)，而脈動周期的誤差能夠控制在7%以內(nèi)。通過與經(jīng)驗公式的對比分析，驗證了水下爆炸氣泡脈動模型建立和數(shù)值計算方法的正確性。

表1 不同工況下氣泡脈動最大半徑與周期仿真結(jié)果與經(jīng)驗公式的對比Tab.1 Comparison computational results and experimental results of the radius and period of bubble pulsation

3.2 兩氣泡數(shù)值計算結(jié)果

為研究兩氣泡水中相互作用的運動規(guī)律，以單一氣泡數(shù)值計算模型為基礎(chǔ)，建立水下爆炸兩氣泡脈動的數(shù)值計算模型，兩氣泡的狀態(tài)方程完全相同。由于氣泡坍塌、射流方向、射流速度等參數(shù)與爆心初始位置、氣泡脈動周期、氣泡最大半徑等關(guān)鍵參數(shù)有密切關(guān)系，因此本節(jié)將研究距離參數(shù)、初始角度對氣泡坍塌、射流角度、射流速度的影響。計算模型示意圖如圖3 所示。

圖3 兩氣泡初始位置示意圖Fig.3 Schemaic diagram of initial locations of two bubbles

當(dāng)兩氣泡連線與水平線夾角α 為0°時，表示兩氣泡沿水平方向排列;當(dāng)夾角α 為90°時，表示兩氣泡沿豎直方向排列;當(dāng)兩氣不處于同一水平面時，相同條件下的氣泡由于所處水深的不同會導(dǎo)致其最大半徑和脈動周期的不同。為方便討論，以兩氣泡連線中點處的參考?xì)馀葑畲蟀霃絉max為特征長度，中點處靜水壓力pamb為特征壓力，則特征時間為t =Rmax(ρ/pamb)1/2，速度的無量綱表達(dá)式為v/(Rmax/t)，兩氣泡中心初始位置之間的距離定義為d =2L/Rmax，即將參考?xì)馀葑畲蟀霃降? 倍作為特征長度，分別研究水平、垂直距離及夾角α 對氣泡運動特性的影響規(guī)律。

3.2.1 水平距離對氣泡運動特性的影響

如圖4 所示為2 個裝藥量相同的TNT 炸藥在相同水深且初始水平距離d =1.0 時其相互作用的演化過程。

圖4 水平距離d 為1.0 時兩氣泡運動演化過程Fig.4 Evolution processes of interaction of two bubbles for horizontal standoff distance d=1.0

由圖4 可看出，當(dāng)初始水平距離d 較小時，氣泡之間的Bjerknes 力影響顯著，明顯強于浮力的影響，表現(xiàn)為坍塌發(fā)生于每個氣泡的外側(cè)，氣泡A 形成由左下指向右上方向的射流，氣泡B 產(chǎn)生右下指向左上方向的射流，最終兩股射流均會擊穿氣泡并相互作用，而氣泡也會由初始的球形變成環(huán)狀氣泡。

由于兩氣泡水平放置時其運動具有對稱性，因此僅研究左側(cè)氣泡的運動特性。如圖5 所示，描述的是左側(cè)氣泡射流角度、速度與氣泡水平距離之間的關(guān)系。

圖5 氣泡水平排列時射流角度、速度和距離之間關(guān)系Fig.5 The relationship among angle and velocity of jet and horizontal distance

由圖5 可看出，氣泡射流角度會隨著水平距離的增加而增加，而射流速度會隨著水平距離的增加而減小。表現(xiàn)為當(dāng)兩氣泡初始距離d 取1.0，即2 倍最大氣泡半徑時，氣泡之間的Bjerknes 力影響顯著，所形成射流的角度較小，約為7°左右，且射流速度也較大，其無量綱速度為5.25;隨著氣泡間水平距離的增加，其氣泡射流與水平線的夾角隨之增加，且增幅明顯，表現(xiàn)為當(dāng)兩氣泡初始距離d 取3.0，即6 倍最大氣泡半徑時，所形成射流角度已為76°，而當(dāng)d 取4.0 時，射流角度約為86°，此時可認(rèn)為氣泡間的Bjerknes 力較弱，浮力占主導(dǎo)，在其作用下形成豎直向上的射流，且與之對應(yīng)的射流速度相對較低，其無量綱量約為1.4.

3.2.2 垂直距離的影響

如圖6 所示為兩氣泡在初始垂直距離d =1.0，即垂直距離為2 倍氣泡最大半徑時氣泡脈動的演化過程。

由圖6 可看出，當(dāng)初始水平距離d 較小時，氣泡之間的Bjerknes 力影響顯著，明顯強于浮力的影響，表現(xiàn)為氣泡B 的坍塌發(fā)生于氣泡的頂部，形成由豎直向下的射流;而氣泡A 在Bjerknes 力浮力共同作用下，坍塌產(chǎn)生于氣泡底部，形成豎直向上的射流，最終射流擊穿氣泡并相互作用，氣泡也會由初始的球形變成環(huán)狀氣泡。

圖6 豎直排列初始距離d 為1.0 時兩氣泡運動演化過程Fig.6 Evolution processes of interaction of two bubbles for vertical standoff distance d=1.0

由于兩氣泡沿垂直方向分布，因此其坍塌與射流方向均會沿豎直方向，為此，主要研究垂直距離對氣泡射流速度的影響。如圖7 所示，描述的是兩氣泡射流速度與垂直距離的變化規(guī)律。

圖7 氣泡豎直排列各氣泡射流速度與距離之間的關(guān)系Fig.7 Relationship between jet velocity and vertical distance

由圖7 可看出，隨著兩氣泡垂直距離的增加，氣泡之間的Bjerknes 力影響逐漸減弱，表現(xiàn)為上下氣泡所產(chǎn)生的射流速度會隨著垂直的距離的增加而減小，特別是上氣泡。當(dāng)垂直距離超過6 倍氣泡最大半徑(d=3.0)時，氣泡頂部的運動速度已接近于氣泡收縮時的速度，此時可認(rèn)為上氣泡頂部不會產(chǎn)生坍塌;而與之對應(yīng)的下氣泡，由于Bjerknes 力與浮力對氣泡作用效果相同，導(dǎo)致下氣泡底部發(fā)生坍塌并形成豎直向上的射流，因此，雖然隨著初始距離的增加下氣泡形成向上射流的速度會降低，但其射流速度會略高于上氣泡頂部的射流速度。

3.2.3 初始角度的影響

與氣泡豎直排列時類似，仍以兩氣泡連線中心點處氣泡為參考，對兩氣泡各參數(shù)進(jìn)行無量綱化。如圖8 和圖9 所示，分別為兩氣泡初始距離為d =1.0 時其最終射流角度和速度與兩氣泡初始夾角之間的關(guān)系曲線。

圖8 d=1.0 時氣泡A 射流角度、速度與初始角度之間關(guān)系Fig.8 Relationship among the velocity and angle of jet and the initial angle of bubble A for d=1.0

圖9 d=1.0 時氣泡B 射流角度、速度與初始角度之間關(guān)系Fig.9 Relationship among the velocity and angle of jet and initial angle of bubble B for d=1.0

由圖8 和圖9 可看出，由于初始距離較小，此刻Bjerknes 力影響顯著，表現(xiàn)為氣泡射流方向基本與原初始角度相同，而射流的速度變化不大且相對較高。由圖5、圖7 距離對射流速度影響規(guī)律可看出:氣泡射流速度隨著距離的增加降低明顯，為此，當(dāng)氣泡初始距離d=3.0 時，氣泡間距離相對較遠(yuǎn)，射流速度相對較低，此時主要考慮初始角度與射流角度之間的關(guān)系。

如圖10 所示為氣泡A 與氣泡B 初始距離d =3.0 時，最終射流角度與氣泡初始角度之間的關(guān)系曲線。

圖10 d=3.0 時氣泡射流角度與初始角度之間關(guān)系Fig.10 Relationship among the velocity and angle of jet and the initial angles of two bubbles for d=3.0

由圖10 可看出，當(dāng)初始距離增大時，氣泡間的Bjerknes 力對氣泡運動影響效果降低，浮力的影響逐漸顯現(xiàn)。隨著氣泡連線與水平線夾角的增加，氣泡A 最終射流的角度由42°逐漸增加至90°，氣泡B最終射流角度由132°逐漸減小至90°，即隨著兩氣泡連線與水平夾角的增加，兩氣泡坍塌所形成的射流方向會逐漸沿豎直方向移動。

4 結(jié)論

1)以氣泡體積加速度模型作為氣泡初始條件定義，結(jié)合自行開發(fā)的定義流場初始條件和邊界條件的子程序，可較準(zhǔn)確地計算出水下爆炸氣泡脈動的全物理過程，且滿足計算精度要求。

2)兩氣泡水平排列時，隨著氣泡間距離的增加，Bjerknes力影響效果逐漸減弱，浮力隨之占主導(dǎo)，表現(xiàn)為隨著氣泡射流角度隨著氣泡間距離的增加而增大，而射流速度會隨著水平距離的增加而減小。

3)兩氣泡豎直排列時，隨著氣泡間距離的增加，Bjerknes力影響效果逐漸減弱，氣泡射流速度逐漸減小，上氣泡會先后形成豎直向下射流、豎直向下與向上的對射流及豎直向上射流，而下氣泡豎直向上的射流速度會逐漸減小。

4)兩氣泡初始時存在夾角時，隨著距離的增加，兩氣泡最終的射流方向會逐漸向豎直方向移動，并最終與豎直方向重合，射流速度會隨著距離的增加而減小。

References)

［1］金輝，張慶明，高春生，等.不同邊界條件水下爆炸氣泡脈動對比的試驗研究［J］. 兵工學(xué)報，2009，12，30(增刊2):213-217.JIN Hui，ZHANG Qing-ming，GAO Chun-shen，et al. Comparison experimental study of underwater explosion bubble pulse among the different boundaries［J］. Acta Armamentarii，2009，12，30(S2):213 -217. (in Chinese)

［2］Klaseboer E，Manica R，Chan D Y C，et al. BEM simulations of potential flow with viscous effects as applied to a rising bubble［J］.Engineering Analysis with Boundary Elements，2011，35(3):489-494.

［3］S Shervani-Tabar M T，Seyed-Sadjadi M H，Shabgard M R.Numerical study on the splitting of a vapor bubble in the ultrasonic assisted EDM process with the curved tool and workpiece［J］.Ultrasonics，2013，53(1):203 -210.

［4］Menon S，Lal M. On the dynamics and instability of bubbles formed during underwater explosions［J］. Experimental thermal and fluid science，1998，16(4):305 -321.

［5］榮吉利，李健. 基于DYTRAN 軟件的三維水下爆炸氣泡運動研究［J］. 兵工學(xué)報，2008，29(3):331 -336.Rong Ji-li，LI Jian. Research on the motion of three-dimensional underwater explosion bubble with DYTRAN software ［J］. Acta Armamentarii，2008，29(3):331 -336. (in Chinese)

［6］Chew L W，Klaseboer E，Ohl S W，et al. Interaction of two oscillating bubbles near a rigid boundary［J］. Experimental Thermal and Fluid Science，2013，44:108 -113.

［7］Rungsiyaphornrat S，Klaseboer E，Khoo B C，et al. The merging of two gaseous bubbles with an application to underwater explosions［J］.Computers ＆ Fluids，2003，32(8):1049 -1074.

［8］李章銳，宗智，董婧，等. 兩個氣泡相互作用的某些動力特性研究［J］. 船舶力學(xué)，2012，16(7):717 -729 LI Zhang-rui，ZONG Zhi，DONG Jing，et al. Some dynamic characteristics of the interactions of two bubbles［J］. Journal of Ship Mechanics，2012，16(7):717 -729. (in Chinese)

［9］王詩平，張阿漫，劉云龍，等. 同相氣泡耦合特性實驗研究［J］.力學(xué)學(xué)報，2012，44(1):56 -64.WANG Shi-ping， ZHANG A-man， LIU Yun-long， et al.Experimental study on interaction of in phase bubbles［J］. Chinese Journal of Theoretical and Applied Mechanics，2012，44(1):56-64. (in Chinese)

［10］Rayleigh L. On the pressure developed in a liquid during the collapse of a spherical cavity［J］. The London，Edinburgh，and Dublin Philosophical Magazine and Journal of Science，1917，34(200):94 -98.

［11］Price R S. Similitude equations for explosives fired underwater，TR 80-299［R］. NSWC:Naval Surface War-fare Center，1979.

［12］Hunter K S. Global-shape-function models of an underwater explosion bubble［D］. Colorado:University of Colorado，2001.

［13］周霖，謝中元，陳勇. 炸藥水下爆炸氣泡脈動周期工程計算方法［J］. 兵工學(xué)報，2009，30(9):1202 -1205.ZHOU Lin，XIE Zhong-yuan，CHEN Yong. An engineering calculation method on the bubble pulse period of underwater explosion［J］. Acta Armamentarii，2009，30(9):1202 -1205.(in Chinese)

［14］庫爾. 水下爆炸［M］. 北京:國防工業(yè)出版社，1960:135 -239.Cole. Underwater explosions［M］. Beijing:National Defense Industry Press，1960:135 -239. (in Chinese)

［15］李健，榮吉利，楊榮杰，等.水中爆炸沖擊波傳播與氣泡脈動的實驗及數(shù)值模擬［J］. 兵工學(xué)報，2008，29(12):1437 -1443.LI Jian，RONG Ji-li，YANG Rong-jie，et al. Experiment and numerical simulation of shock wave propagation and bubble impulse of underwater explosion［J］. Acta Armamentarii，2008，29(12):1437 -1443. (in Chinese)