兩個氣泡相互作用的某些動力特性研究

李章銳 ， 宗 智 ，b， 董 婧 ， 孫 雷 ，b

（大連理工大學a.船舶工程學院；b.工業裝備結構分析國家重點實驗室，遼寧 大連 116023）

兩個氣泡相互作用的某些動力特性研究

李章銳a， 宗 智a，b， 董 婧a， 孫 雷a，b

（大連理工大學a.船舶工程學院；b.工業裝備結構分析國家重點實驗室，遼寧 大連 116023）

文章采用邊界積分方法模擬三維水下氣泡的動力特性，詳細闡述氣泡計算的數學模型和數值實施過程，探究了三種不同浮力參數情況下（無重力、弱浮力、強浮力）氣泡的演變情況。計算結果表明，在不同的浮力參數下，氣泡形狀和射流特點有著顯著的區別，射流形成和發展與重力有著密切的聯系。文中計算了兩氣泡在不同距離和不同強度參數下的相互作用的特點，為氣泡動力特性的研究提供了參考。

邊界積分法；氣泡動力特性；浮力參數；射流

1 引 言

氣泡為我們日常生活所熟知，并在物理、化學、醫學和技術等領域占據著重要地位，然而它的行為特點有時很令人吃驚和難以預測，而且在很多情況下不為人們所知。由于氣泡在流體系統的眾多領域存在，因此氣泡在當今的科學和技術領域占據著十分重要的地位。舉幾個簡單例子：在石油的生產和傳輸過程中，人為地注入氣泡可以把重油從地底輸送到地表；在能源生產中通過沸騰使水變成蒸汽從而推動渦輪轉動的整個過程中，氣泡起著十分重要的作用；化工領域，液氣反應堆里依靠氣泡來增加不同相之間的接觸面積；在噴墨打印中，泡沫室產生的氣泡可以用來標記電離子的運動軌跡。

最早研究氣泡潰滅的研究可以追溯到Rayleigh（1917）[1]，他應用能量方程計算了球形空化氣泡在無粘無旋流體中壓力情況。后來Kornfeld和Suvarov等人（1944）[2]提出射流的形成可能是空化損傷的一個原因。Benjamin和Ellis（1966）[3]為氣泡射流的機理做了進一步的研究，他們用氣泡沖量的概念來解釋氣泡潰滅，射流形成，以及環形氣泡形成等現象。他們的實驗為氣泡由于重力作用或受到附近固壁效應的作用而形成射流提供了很好的驗證。Blake和Gibson（1981）[4]提出了利用近似積分方程來模擬近自由液面附近空化氣泡的生長和潰滅過程，觀察到了氣泡潰滅時射流的加速形成過程和自由液面的明顯水冢現象。Wilkerson（1990）[5]利用邊界積分方法計算了水下爆炸氣泡在剛性固壁附近的潰滅情況。近年來Wang，Zhang等人[8-9]進一步在Wilkerson等人的基礎上，采用邊界積分法（Boundary Integral Method,BIM）模擬了結構物附近三維非對稱水下爆炸氣泡的運動情況，他們提出的光順計算[6]和彈性網格技術[7]使得三維氣泡運動模擬過程變得更加穩定可靠。在國內，張阿漫等人[16]基于勢流理論用高階邊界元法對氣泡的運動特性進行了研究，得到了一些有價值的結論。當前對單個氣泡特性的研究相對較多，但對于多個氣泡研究還不夠深入。本文在前人的研究成果的基礎上，數值模擬兩個氣泡在不同浮力和距離、強度情況下，兩個氣泡的相互作用情況，試圖揭示水下氣泡的某些動態特性。

2 計算模型

假設流體是無粘、不可壓縮的，并且引起的運動是無旋的，因此在流體域內存在勢函數φ x,y,（）z，滿足Laplace方程：

采用Green第二定理，控制方程Laplace方程的邊界積分方程的表達形式為：

其中，V0，P0為初始不可凝結氣體體積和壓強，λ為熱比率。對蒸汽氣泡λ=1.4，對TNT炸藥產生的氣泡，λ=1.25。

為便于計算，所有的變量采用無因次的表達形式。無量綱參數參見表1所示。

圖1 計算坐標系示意圖Fig.1 Schematic plan of coordinate system for calculation

表1 無量綱參數Tab.1 Dimensionless parameters

上表中Rm代表Rayleigh[1]氣泡半徑，相當于在一個無限流場，不考慮重力情況下，球對稱氣泡的最大半徑，d為空間中兩物體相距的實際距離。ΔP=P∞-Pc定義為壓強標尺，即與氣泡等深無窮遠處流體靜水壓力P∞與氣泡內的飽和蒸汽壓力Pc的壓力差，ρ為流體密度，g為重力加速度，P0為氣泡內的初始壓強。

流體域內對應該控制方程的邊界條件為：

運動學條件：

無窮遠邊界條件：

動力學條件：非定常伯努利方程滿足動力學邊界條件：

結合氣體狀態方程，氣泡運動的動力學邊界條件可改寫為：

結合無量綱參數表1，以上動力學邊界條件可以化簡為：

氣泡剛開始，假定在t=0時刻為初速度為0、半徑為R0的高壓氣團，在整個流域內，初始時刻φ=0在氣泡初期運動中，體積很小，浮力相對影響微弱，氣泡的初始狀態可以采用Rayleigh-Plesset[15]（無因次形式）方程描述：

依據該方程可以推導出[10]

其中R為氣泡的無量綱半徑，當給定ε和λ，方程（10）采用非線性方程的牛頓迭代法可以求解出氣泡的初始半徑 R0。 這樣以上方程（2）、（4）、（8）和（10）構成了氣泡運動求解的完整方程。

3 數值方法

對氣泡的數值模擬，本文對其進行三維邊界單元模擬，采用三角單元模擬球形氣泡表面。氣泡表面的初始離散是基于20面體來完成的。20面體形狀由20個等邊三角形和位于球形邊界上的12個節點組成。這種劃分可以通過不斷地把原來的等邊三角形細分成更小等邊三角形來實現球形表面的逼近。每個三角形被細分的次數稱為級，1級精度三角單元即為二十面體，2級精度的三角形單元有42個節點和80個三角形單元（如圖2所示）。隨著劃分級別的不斷提高，原始二十面體的形狀的粗糙度降低，整個多面體的表面逐漸逼近于球面。對于n級精度來說，球面三角形由10n2+2個節點和20n2個單元構成。通常情況下氣泡求解采用6級以上精度能準確描述氣泡的運動狀態。

圖2 不同離散級別的初始氣泡表面示意圖Fig.2 Schematic plan of discretized initial bubble surface with different levels

3.1 離散氣泡邊界的積分計算

（1）氣泡邊界參數的表達

當完成氣泡表面單元劃分后，邊界積分方程（2）可以表示成為所有邊界單元上的積分求和。因此需要引入以正交斜坐標系統來表示氣泡表面上一任意的三角單元ABC將氣泡表面離散化，如圖3（a）所示，氣泡每一個單元的計算，采用線性單元局部坐標系（ξ,η ），把任意的三角形單元ABC映射到ξη平面內的正三角形上。如圖3（b）所示，分別把BC和AC映射到坐標軸η和ξ上，則氣泡的任意一個三角形單元的三個頂點A、B和C分別落在（1,0）、（0,1）和（0,0 ）上，則線性插值表示的幾何位置函數 r，速度勢 φ 和法向速度 ψ=?φ/?n 如下：

圖3 離散的氣泡邊界及單元映射Fig.3 The discrete bubble surface and the element coordinate mapping

將（11）、（12）式和（13）式代入方程（2）右端項，則第一部分在三角形ABC上的積分表達為：

其余部分表達為：

即Ai（i=1,3）和Bi（i=1,3 ）為影響系數；n是三角形ABC的外法向單位向量；（14）式和（15）式中各參數由（16）式至（29）式給出。

系數矩陣A和B為影響系數是整個求解的關鍵，也是耗費計算機CPU運算的主要部分。系數矩陣A和B的非對角元素采用7點高斯積分進行計算。而對于系數矩陣A和B矩陣的對角元素的積分存在奇異性。由于A矩陣對角元素是弱奇異，通過將直角坐標系求解轉化到極坐標條件下求解即可以消除積分的奇異性[11]。B矩陣的對角元素具有強奇異，可采用4π法則求解,詳細推導可參見Klaseboer文獻 [12]，即：

當邊界的初始位置和邊界上初始速度勢（對于氣泡表面，初始時刻的速度勢φ=0）給定后，通過求解Laplace方程（2）可以得到速度勢法向速度?φ/?n，然后通過有限差分法[6]求解得到氣泡的切向速度，當邊界上的▽φ求得后，利用公式（4）和（8）可以得到下一時刻的氣泡節點的位置和速度勢，求解過程由于會產生數值不穩定，需采用光順技術[6]和EMT[7]技術。

4 計算模型驗證

4.1 與普通非爆炸氣泡的對比

Lawson[13]曾采用質點影像測速技術（簡稱PIV技術）在實驗中模擬過水下普通氣泡破裂過程。他們在尺寸為600×600×800（mm）水箱中給一橡膠氣球充氣至300 cm3的大小，以模擬真實氣泡，然后用一長釘刺穿氣球，使之產生初始運動，再由PIV設備記錄氣泡在水下運動各時刻的情況。為模擬真實情況，本文數值計算中參數依照實驗設置如下，如表2所示。

表2 氣泡初始參數設置Tab.2 Parameters of initial bubble

圖4 50，60，70，80ms時刻計算得到氣泡形狀與試驗對比（虛線代表實驗數據，實線代表數值計算）Fig.4 The calculated bubble shapes compared with those of experiment at 50,60,70,80 ms

數值計算采用6級精度，初始氣泡表面離散為362個節點和720個三角單元。氣泡運動的實驗與仿真結果如圖4所示。圖中依次給出了50 ms、60 ms、70 ms和80 ms時刻的氣泡運動剖面圖。實線代表數值仿真的結果，虛線代表實驗測得的結果，從圖中可以看到，仿真過程中考慮了重力的影響，與實驗條件一致。整個過程中氣泡都在豎直向上運動，初始時刻運動較為緩慢，隨著時間的增長，開始加速運動，由球形變為非球形，在t=50 ms時，氣泡底部出現明顯凹陷，在t=60 ms和70 ms時刻，氣泡繼續上浮，其底部的凹陷進一步加深，逐漸形成向上的射流，在t=80 ms時刻，射流即將沖破氣泡對面的壁面，在實驗中氣泡在90 ms被射流沖破，而數值仿真中是92 ms，結果十分接近真實值。整個數值模擬與實驗基本吻合，說明邊界積分方法能夠較為真實地模擬出氣泡運動情況。數值仿真中我們也發現數值模擬與實際觀察仍然存在一些差別，實驗中橡膠氣球中氣體有5%到10%的壓縮而實際仿真中暫沒有考慮。此外，數值仿真中存在著截斷誤差和累計誤差（在下個實驗對比中可以看到），這些都是誤差產生的重要原因，但整體說來，邊界積分法能夠較為準確地模擬出氣泡運動的真實狀況。

4.2 與Rayleigh 氣泡半徑的對比

在不考慮重力影響的情況下，在無限流域（忽略自由液面和固壁效應）中的一個振動氣泡的運動狀況可以用Rayleigh-Plesset方程描述，如（9）式所示。該方程可以通過相應數值算法，如四階Runge-Kutta或五階Runge-Kutta-Felhberg[14]（RKF），本文采用五階RKF求解。在BIM數值仿真中氣泡采用7級精度，初始氣泡表面離散為492個節點和980個三角單元，如圖5所示。求解參數如下。

表3 Rayleigh氣泡初始參數設置Tab.3 Initial parameters for Rayleigh bubble

從計算結果看，實線代表采用RKF求解得到的氣泡運動半徑，星號代表的是用邊界積分法（BIM）數值仿真得到的氣泡半徑，綠色（Err）虛線代表兩者的誤差，可以看到兩者吻合良好，說明邊界積分法能夠很好地模擬無限流場的氣泡運動狀況，從誤差比較可以看到，在氣泡膨脹期間兩者的誤差較小，在氣泡收縮時兩者誤差較大，而且隨著求解周期的增加，從圖中可以看到，到第三個周期，誤差增長已十分迅速。采用BIM計算得到的結果誤差會逐漸增大。但通常情況下數值仿真求解氣泡周期為1-2個（第二個周期后的氣泡運動狀態有待進一步探究），因此不會對求解結果造成太大誤差。

圖5 采用RKF方法和BIM方法求得瑞利氣泡半徑隨時間的變化Fig.5 The radius variation and error accumulation of Rayleigh bubble obtained by RKF and BIM

5 兩個氣泡的某些特性研究

5.1 浮力參數對氣泡運動特性的影響

當浮力參數δ比較小的時候，氣泡所受的浮力的影響可以忽略不計，這種近似適用于大壓力（深水情況下）的體積較小的氣泡。為此我們計算了在不同浮力參數下（δ為零，δ較小，δ較大情況下）氣泡的運動特點。這里浮力參數分別選取δ2=0，δ2=0.05和δ2=0.25計算條件，計算氣泡的運動及其射流情況，以下參數均為無因次參數。

（a） 無重力（δ2=0）

首先考慮100 m水深下的兩氣泡，由于水深較深，自由液面的影響可以忽略，初始半徑均為0.139 1，最大半徑為1，氣泡中心相距2，強度參數ε=120，浮力參數為δ=0（即零重力），氣泡精度為12級（即1 442個節點，2 880個單元 ），光順技術每10實施一次，并采用EMT技術更新氣泡位置和速度勢，圖6給出了幾個典型時刻氣泡的計算結果圖。

初始時刻兩氣泡為一高溫高壓球體，由于氣泡內部壓力遠遠高于外界水壓力，氣泡開始膨脹，在t=0.210 6時刻看到氣泡體積明顯增長，但依舊保持球形，在t=1.133 4時刻，氣泡體積達到最大，在兩氣泡相鄰的部分氣泡已不再保持球形，微微呈現出扁平形狀，說明氣泡間發生了相互排斥。隨后氣泡體積開始縮小，在t=1.917 8時刻，中部出現明顯的扁平，排斥作用進一步加強，在t=2.184 0時刻兩氣泡外側出現水平凹陷，開始形成水平射流，在t=2.253 8時刻，氣泡體積繼續縮小，射流方向依舊保持水平并進一步發展，即將要沖向氣泡的對面壁面。

（b） 弱浮力（δ2=0.05）

在弱浮力即δ2=0.05條件下，在初始時刻及膨脹的初期（t=0.089 7時刻）氣泡保持球形，在t=1.270 7時刻氣泡體積達到最大，但中部仍然出現扁平，兩氣泡發生排斥作用。隨后氣泡體積開始減小，在t=1.941 7時刻，氣泡已不再保持球形，遠離兩氣泡的部分的下側出現凹陷，在t=2.194 1時刻遠離兩氣泡的部分的下側出現斜向上的射流，在t=2.237 4斜向上的射流進一步發展，沖向氣泡壁面。射流發展的方向是重力和兩氣泡之間的吸引力共同作用的結果，但由于重力作用較小，寬度較窄，規模較小，如圖7所示。

圖6 無重力狀態下兩氣泡的演化過程Fig.6 The evolution of two bubbles without gravity

圖7 弱浮力狀態下兩氣泡的演化過程Fig.7 The evolution of two bubbles with weak buoyancy effect

（c） 強浮力（δ2=0.25）

在強浮力即δ2=0.25條件下，氣泡在膨脹期間與無浮力和弱浮力條件類似，依舊是球形膨脹，如t=0.202 4時刻所示。在t=1.073 9時刻，氣泡繼續膨脹，接近氣泡的最大體積，當兩氣泡距離足夠近時氣泡相鄰部分出現排斥。在t=1.6434時刻氣泡底部出現明顯凹陷，在t=1.954 2時刻氣泡底部出現斜向上射流，與弱浮力條件相比，射流的寬度明顯增加，規模增大。在t=2.227 9時刻，射流沖向對面的氣泡壁面。與弱浮力條件下相比，強浮力狀況下射流形成的時間明顯提前，如圖8所示。

圖8 強浮力狀態下兩氣泡的演化過程Fig.8 The evolution of two bubbles with strong buoyancy effect

（d） 參數分析

圖9給出了不同浮力參數下兩氣泡中心橫向位置的變化，從圖中可以看到無論在無重力還是弱浮力或強浮力狀況下，兩氣泡的距離都是先增大后減小，也就是在膨脹期間，兩氣泡間發生排斥，氣泡之間距離增大，在氣泡潰滅期間兩氣泡發生吸引，氣泡之間距離減小。圖10給出了不同浮力參數下氣泡中心垂向位置的變化，可以看到隨著浮力參數的增加，氣泡上升的速度加快。圖11給出了不同浮力參數下體積隨時間的變化。氣泡在膨脹過程中體積增長基本一致，但潰滅時隨著浮力參數增加，氣泡體積收縮的速度有所減慢。圖12給出了氣泡在不同浮力下動能（Kinetic）、勢能（Potential）和總能量（Total）隨時間的變化，在氣泡膨脹時，氣泡的動能逐漸降低而勢能逐漸升高，而氣泡潰滅時，氣泡的動能逐漸升高而勢能逐漸降低，氣泡的總能量在整個氣泡運動過程中是守恒的。

圖9 氣泡中心橫向距離的變化Fig.9 The variation of transverse distance of two bubbles

圖10 氣泡中心垂向位置的變化Fig.10 The variation of vertical position of bubble

圖11 氣泡體積隨時間的變化Fig.11 The bubble volumes as a function of time

圖12 氣泡能量隨時間的變化Fig.12 The bubble energy as a function of time

圖13給出了不同浮力參數下射流點速度的變化，可以看到初始階段氣泡射流點速度隨浮力參數的增大而增大，對于無重力或弱浮力 （δ2= 0,0.05 ）情況，隨著時間的增長，氣泡射流點速度逐步增大，會達到一個較高峰值。對于浮力參數較大 （δ2=0.25,0.3,0.4）時射流點速度卻經歷了先增大后減小的過程，射流點速度最大峰值相對于無重力或弱浮力情況小，也就是說，浮力較小的情況下，射流的寬度較窄，速度可能很大，浮力較大時，射流的寬度較大，但速度峰值卻要小一些。

5.2 橫向距離參數和強度參數對氣泡運動特性的影響

（a）不同橫向距離下的比較

為進一步考察氣泡橫向距離對氣泡運動特性的影響。在無重力情況下，我們計算了氣泡中心橫向（x方向

）距離為兩倍、三倍、五倍和七倍最大半徑情況下氣泡的運動情況，計算的其他參數和5.1節（δ2=0）中相同。計算結果表明兩氣泡在不同橫向距離的演化形態與兩倍條件下類似。不同的是，從圖14可以看到隨著兩氣泡距離的增大，兩氣泡中心變化的幅值越來越小，也就是說氣泡間作用排斥和吸引越來越小，在五倍以上的最大半徑距離時，氣泡的作用已比較微弱了。圖15給出了不同距離下氣泡體積隨時間的變化，隨著距離的增大氣泡的收縮速率加快，距離越大，氣泡射流點后期體積越小。氣泡體積變化率（體積對時間求導）如圖16所示，從圖中可以看到氣泡膨脹初期和潰滅末期氣泡變化速率都較為劇烈。圖17給出了不同距離下氣泡射流點速度隨時間的變化，射流點初期在不同距離下無明顯差別，后期的速度隨距離的增大而增大。

圖13 不同浮力參數下射流點速度的變化Fig.13 The jet-tip velocity of bubble under different buoyancy parameters

圖14 氣泡中心橫向距離的變化Fig.14 The variation of transverse distance of two bubbles

圖15 氣泡體積隨時間的變化Fig.15 The bubble volumes as a function of time

圖16 氣泡體積變化率隨時間的變化Fig.16 The bubble volume change rate as a function of time

圖17 不同距離下射流點速度的變化Fig.17 The maximum velocity of bubble under different distances

（b）不同強度參數下的比較

為進一步考察氣泡強度參數對氣泡運動特性的影響，在無重力情況下，我們計算不同強度ε=20，60，120，180參數情況下氣泡的運動情況，計算的其他參數和5.1節（δ2=0）中相同。氣泡演化的形態與ε=120的類似。從圖18可以看到，氣泡運動過程中，兩氣泡中心橫向距離都有先緩慢增大后迅速減小的規律，但隨著兩氣泡強度參數的增大，兩氣泡中心橫向距離變化的幅值增加，也就是說氣泡間排斥和吸引作用增強。圖19給出了不同強度參數下氣泡體積隨時間的變化，氣泡強度增大，氣泡的運動周期明顯增大，氣泡膨脹到最大體積的時間會越晚，圖20給出氣泡體積變化率隨時間的變化，氣泡在膨脹初期和潰滅末期的運動都十分劇烈。圖21給出了不同強度參數下氣泡射流點速度隨時間的變化，射流點在初期速度變化較為平緩，在末期射流速度急劇增加；隨著強度參數增加，射流發生穿透的時間越晚。

圖18 氣泡中心橫向距離的變化Fig.18 The variation of transverse distance of two bubbles

圖19 氣泡體積隨時間的變化Fig.19 The bubble volumes as a function of time

圖20 氣泡體積變化率隨時間的變化Fig.20 The bubble volume change rate as a function of time

圖21 不同強度參數下射流點速度的變化Fig.21 The maximum velocity of bubble under different strength parameters

6 結論與討論

本論文中采用三維邊界積分模擬水下兩個氣泡相互作用的問題，詳細論述了采用邊界積分法求解氣泡運動的算法和數值實施過程，并總結出以下結論：

（1）邊界積分法（BIM）能夠較好地模擬實際氣泡的運動狀況，無論是非爆炸氣泡還是爆炸氣泡都能準確地描述和模擬出氣泡的運動狀態和各個時刻氣泡特征參數變化。

（2）無論是無重力、弱浮力還是強浮力條件下，氣泡兩者之間主要表現為膨脹時相互排斥，潰滅時相互吸引。與無重力和弱浮力比較，強浮力作用下射流的出現時間明顯提前，并且射流的寬度和規模都更大，也就是說，重力對射流形成和發展有著重要的影響。無重力和弱浮力條件下，射流速度隨時間逐步增加并在末期急劇上升，強浮力作用下，射流速度經歷了由大變小的過程，強浮力條件下的射流速度峰值小于無重力和弱浮力條件下的。

（3）氣泡在膨脹初期和潰滅末期的變化都十分劇烈，兩氣泡的作用隨著距離參數的減小和強度參數的增大越來越明顯，在五倍以上的距離參數條件下，兩氣泡的作用已比較微弱了。無重力條件下，在不同距離下射流點速度在初期無明顯差別，后期的速度隨距離的增大而增大；不同強度參數下，氣泡射流點在初期速度變化較為平緩，在末期射流速度急劇增加；隨著強度參數增加，射流發生穿透的時間越晚。

[1]Rayleigh L.On the pressure developed in a liquid during the collapse of a spherical cavity[J].Phil.Mag.,1917,34:94-98.

[2]Kornfeld M,Suvarov L.On the destructive action of cavitation[J].Journal of Applied Physics,1944,15:495-506.

[3]Benjamin T B,Eillis A T.The collapse of cavitation bubbles and the pressures thereby produced against solid boundaries[J].Phil.Trans.Roy.Soc.,London,Ser.A,1966,260:221-240.

[4]Blake J R,Gibson D C.Growth and collapse of a vapour cavity near a free surface[J].J Fluid Mach.,1981,111:123-140.

[5]Wilkerson S A.A boundary integral approach to three-dimensional underwater explosion bubble dynamics[D].PHD.Thesis.The Johns Hopkins University,1990.

[6]Zhang Y L,Yeo K S,Khoo B C,Wang C.3D jet impact and toroidal bubbles[J].Journal of Computational Physics,2001,166:336-360.

[7]Wang C,Khoo B C,Yeo K S.Elastic mesh technique for 3D BIM simulation with an application to underwater explosion bubble dynamics[J].Computers and Fluids,2003,32:1195-1212.

[8]Wang Q X.The evolution of gas bubble near an inclined wall[J].Theoretical and Computational Fluid Dynamics,1998,12:29-51.

[9]Zhang Y L,Yeo K S,Khoo B C,Chong W K.Three-dimensional bubbles near a free surface[J].Comput.Phy.,1998,146:105-123.

[10]Klaseboer E,Khoo B C.A modified Rayleigh Plesset model for a non-spherically symmetric oscillating bubble with applications to boundary integral methods[J].Eng.Anal.Bound.Elem.,2006,30:59-71.

[11]Wilkerson S A.Boundary integral technique for explosion bubble collapse analysis[C]//Proc.ASME Energy-Source Technology Conference and Exhibition.Houston,Texas,1989.

[12]Klaseboer E,Fernandez C R,Khoo B C.A note on true desingularisation of boundary integral methods for three-dimensional potential problems[J].Engineering Analysis with Boundary Elements,2009,33,6:796-801.

[13]Lawson N J,Rudman M,Guerra A,Liow J L.Experimental and numerical comparisons of a large bubble[J].Experiments in Fluids,1999,26:524-534.

[14]Alehossein H,Qin Zongyi.Numerical analysis of Rayleigh-Plesset equation for cavitating water jets[J].Int.J Numer.Meth.Engng.,2007,72:780-807.

[15]Best J P,Kucera A.A numerical investigation of non-spherical rebounding bubbles[J].Fluid Mech,1992,245:137-154.

[16]Zhang A M,Yao X L,Yu X B.The dynamics of three-dimensional underwater explosion bubble[J].Journal of Sound and Vibration,2008,311(3-5):1196-1212.

Some dynamic characteristics of the interactions of two bubbles

LI Zhang-ruia,ZONG Zhia,b,DONG Jinga,SUN Leia,b

(a.School of Naval Architecture;b.State Key Laboratory of Structural Analysis for Industrial Equipment,Dalian University of Technology,Dalian 116023,China)

Dynamic characteristics of three-dimensional underwater bubbles are simulated using boundary integral method.The mathematical model and numerical implementation are presented in detail.The evolution of the two bubbles under three different buoyancy parameters(no gravity,weak buoyancy parameter and strong buoyancy parameter)was investigated.The calculated results indicate that the bubble shapes and jet characteristics are very different under different buoyancy parameters.The jet formation and development have very close relationship with gravity effect.The interactions of two bubbles under different distances and strength parameters are also simulated.This study could provide a reference to the study of bubble dynamics.

boundary integral method;bubble dynamics;buoyancy parameter;jet formation

O351

A

1007-7294（2012）07-0717-13

2011-11-15 修改日期：2012-02-26

創新研究群體科學基金(50921001)；國家重點基礎研究發展計劃項目（2010CB83270）

李章銳（1985-），男，大連理工大學船舶工程學院博士研究生，E-mail:lizhangruix@yahoo.com.cn；宗 智（1964-），男，教授，博士生導師。