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

（大連理工大學(xué)a.船舶工程學(xué)院；b.工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023）

兩個(gè)氣泡相互作用的某些動(dòng)力特性研究

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

（大連理工大學(xué)a.船舶工程學(xué)院；b.工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023）

文章采用邊界積分方法模擬三維水下氣泡的動(dòng)力特性，詳細(xì)闡述氣泡計(jì)算的數(shù)學(xué)模型和數(shù)值實(shí)施過(guò)程，探究了三種不同浮力參數(shù)情況下（無(wú)重力、弱浮力、強(qiáng)浮力）氣泡的演變情況。計(jì)算結(jié)果表明，在不同的浮力參數(shù)下，氣泡形狀和射流特點(diǎn)有著顯著的區(qū)別，射流形成和發(fā)展與重力有著密切的聯(lián)系。文中計(jì)算了兩氣泡在不同距離和不同強(qiáng)度參數(shù)下的相互作用的特點(diǎn)，為氣泡動(dòng)力特性的研究提供了參考。

邊界積分法；氣泡動(dòng)力特性；浮力參數(shù)；射流

1 引 言

氣泡為我們?nèi)粘Ｉ钏熘⒃谖锢?、化學(xué)、醫(yī)學(xué)和技術(shù)等領(lǐng)域占據(jù)著重要地位，然而它的行為特點(diǎn)有時(shí)很令人吃驚和難以預(yù)測(cè)，而且在很多情況下不為人們所知。由于氣泡在流體系統(tǒng)的眾多領(lǐng)域存在，因此氣泡在當(dāng)今的科學(xué)和技術(shù)領(lǐng)域占據(jù)著十分重要的地位。舉幾個(gè)簡(jiǎn)單例子：在石油的生產(chǎn)和傳輸過(guò)程中，人為地注入氣泡可以把重油從地底輸送到地表；在能源生產(chǎn)中通過(guò)沸騰使水變成蒸汽從而推動(dòng)渦輪轉(zhuǎn)動(dòng)的整個(gè)過(guò)程中，氣泡起著十分重要的作用；化工領(lǐng)域，液氣反應(yīng)堆里依靠氣泡來(lái)增加不同相之間的接觸面積；在噴墨打印中，泡沫室產(chǎn)生的氣泡可以用來(lái)標(biāo)記電離子的運(yùn)動(dòng)軌跡。

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

2 計(jì)算模型

假設(shè)流體是無(wú)粘、不可壓縮的，并且引起的運(yùn)動(dòng)是無(wú)旋的，因此在流體域內(nèi)存在勢(shì)函數(shù)φ x,y,（）z，滿足Laplace方程：

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

其中，V0，P0為初始不可凝結(jié)氣體體積和壓強(qiáng)，λ為熱比率。對(duì)蒸汽氣泡λ=1.4，對(duì)TNT炸藥產(chǎn)生的氣泡，λ=1.25。

為便于計(jì)算，所有的變量采用無(wú)因次的表達(dá)形式。無(wú)量綱參數(shù)參見(jiàn)表1所示。

圖1 計(jì)算坐標(biāo)系示意圖Fig.1 Schematic plan of coordinate system for calculation

表1 無(wú)量綱參數(shù)Tab.1 Dimensionless parameters

上表中Rm代表Rayleigh[1]氣泡半徑，相當(dāng)于在一個(gè)無(wú)限流場(chǎng)，不考慮重力情況下，球?qū)ΨQ氣泡的最大半徑，d為空間中兩物體相距的實(shí)際距離。ΔP=P∞-Pc定義為壓強(qiáng)標(biāo)尺，即與氣泡等深無(wú)窮遠(yuǎn)處流體靜水壓力P∞與氣泡內(nèi)的飽和蒸汽壓力Pc的壓力差，ρ為流體密度，g為重力加速度，P0為氣泡內(nèi)的初始?jí)簭?qiáng)。

流體域內(nèi)對(duì)應(yīng)該控制方程的邊界條件為：

運(yùn)動(dòng)學(xué)條件：

無(wú)窮遠(yuǎn)邊界條件：

動(dòng)力學(xué)條件：非定常伯努利方程滿足動(dòng)力學(xué)邊界條件：

結(jié)合氣體狀態(tài)方程，氣泡運(yùn)動(dòng)的動(dòng)力學(xué)邊界條件可改寫(xiě)為：

結(jié)合無(wú)量綱參數(shù)表1，以上動(dòng)力學(xué)邊界條件可以化簡(jiǎn)為：

氣泡剛開(kāi)始，假定在t=0時(shí)刻為初速度為0、半徑為R0的高壓氣團(tuán)，在整個(gè)流域內(nèi)，初始時(shí)刻φ=0在氣泡初期運(yùn)動(dòng)中，體積很小，浮力相對(duì)影響微弱，氣泡的初始狀態(tài)可以采用Rayleigh-Plesset[15]（無(wú)因次形式）方程描述：

依據(jù)該方程可以推導(dǎo)出[10]

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

3 數(shù)值方法

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

圖2 不同離散級(jí)別的初始?xì)馀荼砻媸疽鈭DFig.2 Schematic plan of discretized initial bubble surface with different levels

3.1 離散氣泡邊界的積分計(jì)算

（1）氣泡邊界參數(shù)的表達(dá)

當(dāng)完成氣泡表面單元?jiǎng)澐趾螅吔绶e分方程（2）可以表示成為所有邊界單元上的積分求和。因此需要引入以正交斜坐標(biāo)系統(tǒng)來(lái)表示氣泡表面上一任意的三角單元ABC將氣泡表面離散化，如圖3（a）所示，氣泡每一個(gè)單元的計(jì)算，采用線性單元局部坐標(biāo)系（ξ,η ），把任意的三角形單元ABC映射到ξη平面內(nèi)的正三角形上。如圖3（b）所示，分別把BC和AC映射到坐標(biāo)軸η和ξ上，則氣泡的任意一個(gè)三角形單元的三個(gè)頂點(diǎn)A、B和C分別落在（1,0）、（0,1）和（0,0 ）上，則線性插值表示的幾何位置函數(shù) r，速度勢(shì) φ 和法向速度 ψ=?φ/?n 如下：

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

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

其余部分表達(dá)為：

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

系數(shù)矩陣A和B為影響系數(shù)是整個(gè)求解的關(guān)鍵，也是耗費(fèi)計(jì)算機(jī)CPU運(yùn)算的主要部分。系數(shù)矩陣A和B的非對(duì)角元素采用7點(diǎn)高斯積分進(jìn)行計(jì)算。而對(duì)于系數(shù)矩陣A和B矩陣的對(duì)角元素的積分存在奇異性。由于A矩陣對(duì)角元素是弱奇異，通過(guò)將直角坐標(biāo)系求解轉(zhuǎn)化到極坐標(biāo)條件下求解即可以消除積分的奇異性[11]。B矩陣的對(duì)角元素具有強(qiáng)奇異，可采用4π法則求解,詳細(xì)推導(dǎo)可參見(jiàn)Klaseboer文獻(xiàn) [12]，即：

當(dāng)邊界的初始位置和邊界上初始速度勢(shì)（對(duì)于氣泡表面，初始時(shí)刻的速度勢(shì)φ=0）給定后，通過(guò)求解Laplace方程（2）可以得到速度勢(shì)法向速度?φ/?n，然后通過(guò)有限差分法[6]求解得到氣泡的切向速度，當(dāng)邊界上的▽?duì)涨蟮煤?，利用公式?）和（8）可以得到下一時(shí)刻的氣泡節(jié)點(diǎn)的位置和速度勢(shì)，求解過(guò)程由于會(huì)產(chǎn)生數(shù)值不穩(wěn)定，需采用光順技術(shù)[6]和EMT[7]技術(shù)。

4 計(jì)算模型驗(yàn)證

4.1 與普通非爆炸氣泡的對(duì)比

Lawson[13]曾采用質(zhì)點(diǎn)影像測(cè)速技術(shù)（簡(jiǎn)稱PIV技術(shù)）在實(shí)驗(yàn)中模擬過(guò)水下普通氣泡破裂過(guò)程。他們?cè)诔叽鐬?00×600×800（mm）水箱中給一橡膠氣球充氣至300 cm3的大小，以模擬真實(shí)氣泡，然后用一長(zhǎng)釘刺穿氣球，使之產(chǎn)生初始運(yùn)動(dòng)，再由PIV設(shè)備記錄氣泡在水下運(yùn)動(dòng)各時(shí)刻的情況。為模擬真實(shí)情況，本文數(shù)值計(jì)算中參數(shù)依照實(shí)驗(yàn)設(shè)置如下，如表2所示。

表2 氣泡初始參數(shù)設(shè)置Tab.2 Parameters of initial bubble

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

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

4.2 與Rayleigh 氣泡半徑的對(duì)比

在不考慮重力影響的情況下，在無(wú)限流域（忽略自由液面和固壁效應(yīng)）中的一個(gè)振動(dòng)氣泡的運(yùn)動(dòng)狀況可以用Rayleigh-Plesset方程描述，如（9）式所示。該方程可以通過(guò)相應(yīng)數(shù)值算法，如四階Runge-Kutta或五階Runge-Kutta-Felhberg[14]（RKF），本文采用五階RKF求解。在BIM數(shù)值仿真中氣泡采用7級(jí)精度，初始?xì)馀荼砻骐x散為492個(gè)節(jié)點(diǎn)和980個(gè)三角單元，如圖5所示。求解參數(shù)如下。

表3 Rayleigh氣泡初始參數(shù)設(shè)置Tab.3 Initial parameters for Rayleigh bubble

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

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

5 兩個(gè)氣泡的某些特性研究

5.1 浮力參數(shù)對(duì)氣泡運(yùn)動(dòng)特性的影響

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

（a） 無(wú)重力（δ2=0）

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

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

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

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

圖6 無(wú)重力狀態(tài)下兩氣泡的演化過(guò)程Fig.6 The evolution of two bubbles without gravity

圖7 弱浮力狀態(tài)下兩氣泡的演化過(guò)程Fig.7 The evolution of two bubbles with weak buoyancy effect

（c） 強(qiáng)浮力（δ2=0.25）

在強(qiáng)浮力即δ2=0.25條件下，氣泡在膨脹期間與無(wú)浮力和弱浮力條件類似，依舊是球形膨脹，如t=0.202 4時(shí)刻所示。在t=1.073 9時(shí)刻，氣泡繼續(xù)膨脹，接近氣泡的最大體積，當(dāng)兩氣泡距離足夠近時(shí)氣泡相鄰部分出現(xiàn)排斥。在t=1.6434時(shí)刻氣泡底部出現(xiàn)明顯凹陷，在t=1.954 2時(shí)刻氣泡底部出現(xiàn)斜向上射流，與弱浮力條件相比，射流的寬度明顯增加，規(guī)模增大。在t=2.227 9時(shí)刻，射流沖向?qū)γ娴臍馀荼诿妗Ｅc弱浮力條件下相比，強(qiáng)浮力狀況下射流形成的時(shí)間明顯提前，如圖8所示。

圖8 強(qiáng)浮力狀態(tài)下兩氣泡的演化過(guò)程Fig.8 The evolution of two bubbles with strong buoyancy effect

（d） 參數(shù)分析

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

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

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

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

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

圖13給出了不同浮力參數(shù)下射流點(diǎn)速度的變化，可以看到初始階段氣泡射流點(diǎn)速度隨浮力參數(shù)的增大而增大，對(duì)于無(wú)重力或弱浮力 （δ2= 0,0.05 ）情況，隨著時(shí)間的增長(zhǎng)，氣泡射流點(diǎn)速度逐步增大，會(huì)達(dá)到一個(gè)較高峰值。對(duì)于浮力參數(shù)較大 （δ2=0.25,0.3,0.4）時(shí)射流點(diǎn)速度卻經(jīng)歷了先增大后減小的過(guò)程，射流點(diǎn)速度最大峰值相對(duì)于無(wú)重力或弱浮力情況小，也就是說(shuō)，浮力較小的情況下，射流的寬度較窄，速度可能很大，浮力較大時(shí)，射流的寬度較大，但速度峰值卻要小一些。

5.2 橫向距離參數(shù)和強(qiáng)度參數(shù)對(duì)氣泡運(yùn)動(dòng)特性的影響

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

為進(jìn)一步考察氣泡橫向距離對(duì)氣泡運(yùn)動(dòng)特性的影響。在無(wú)重力情況下，我們計(jì)算了氣泡中心橫向（x方向

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

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

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

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

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

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

（b）不同強(qiáng)度參數(shù)下的比較

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

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

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

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

圖21 不同強(qiáng)度參數(shù)下射流點(diǎn)速度的變化Fig.21 The maximum velocity of bubble under different strength parameters

6 結(jié)論與討論

本論文中采用三維邊界積分模擬水下兩個(gè)氣泡相互作用的問(wèn)題，詳細(xì)論述了采用邊界積分法求解氣泡運(yùn)動(dòng)的算法和數(shù)值實(shí)施過(guò)程，并總結(jié)出以下結(jié)論：

（1）邊界積分法（BIM）能夠較好地模擬實(shí)際氣泡的運(yùn)動(dòng)狀況，無(wú)論是非爆炸氣泡還是爆炸氣泡都能準(zhǔn)確地描述和模擬出氣泡的運(yùn)動(dòng)狀態(tài)和各個(gè)時(shí)刻氣泡特征參數(shù)變化。

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

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

[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

創(chuàng)新研究群體科學(xué)基金(50921001)；國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃項(xiàng)目（2010CB83270）

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