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

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

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

0 引言

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

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

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

1 理論分析

1.1 氣泡初始條件

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

式中:V為氣泡體積;ρ為流體密度;Pc=;mc、ac分別為藥包的藥量和半徑;K1、K2、A、B 為與炸藥有關的材料參數。對于密度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時間段內，氣泡半徑及徑向速度的時程曲線。

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

炸藥引爆后，壓力瞬間增加，故初始時刻氣泡徑向速度變化率非常明顯，之后的沖擊波傳播過程中氣泡膨脹速度逐漸平緩。針對這一規律，由圖1 可見，當t≥0.8 ms 時，氣泡徑向速度的變化已經很小了，即可認為此時初始氣泡形成。本文，t 取0.8 ms時對應的氣泡膨脹速度為174.5 m/s，初始半徑為0.303 m.

1.2 氣泡運動的經驗公式

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

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

2 數值模擬

2.1 有限元建模

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

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

為在有限區域內真實反映流場特性，需定義流場的靜水壓力以及與之對應的邊界條件。由于MSC. DYTRAN 有限元軟件并不能通過直接定義重力加速度實現流場靜水壓力的定義。因此，需根據每個網格所處的位置，通過EXINIT 子程序將靜水壓力值逐一賦到對應的網格上。同樣，為了能夠在有限的計算區域真實模擬無限水域流場特點，即當流場邊界兩側出現壓力差時，流體會從高壓區通過邊界流向低壓區，可利用EXFLOW2 子程序將不同的靜水壓力值賦到流場邊界上，最終實現三維流場靜水壓力與邊界條件的定義。

2.3 狀態方程

文中忽略氣泡運動對氣體壓力的影響，認為氣體的壓力僅和氣泡的初始狀態及其體積有關，即氣泡內的壓力p 與氣泡體積V 的關系:

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

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

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

3 計算結果及討論

3.1 單個氣泡數值計算結果

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

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

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

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

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

3.2 兩氣泡數值計算結果

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

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

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

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 較小時，氣泡之間的Bjerknes 力影響顯著，明顯強于浮力的影響，表現為坍塌發生于每個氣泡的外側，氣泡A 形成由左下指向右上方向的射流，氣泡B 產生右下指向左上方向的射流，最終兩股射流均會擊穿氣泡并相互作用，而氣泡也會由初始的球形變成環狀氣泡。

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

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

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

3.2.2 垂直距離的影響

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

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

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

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

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

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

3.2.3 初始角度的影響

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

圖8 d=1.0 時氣泡A 射流角度、速度與初始角度之間關系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 射流角度、速度與初始角度之間關系Fig.9 Relationship among the velocity and angle of jet and initial angle of bubble B for d=1.0

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

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

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

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

4 結論

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

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

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

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

References)

［1］金輝，張慶明，高春生，等.不同邊界條件水下爆炸氣泡脈動對比的試驗研究［J］. 兵工學報，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］. 兵工學報，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］. 船舶力學，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］.力學學報，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］. 兵工學報，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］. 北京:國防工業出版社，1960:135 -239.Cole. Underwater explosions［M］. Beijing:National Defense Industry Press，1960:135 -239. (in Chinese)

［15］李健，榮吉利，楊榮杰，等.水中爆炸沖擊波傳播與氣泡脈動的實驗及數值模擬［J］. 兵工學報，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)