倪寶玉，張阿漫，張忠宇，李 帥

（哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001）

復(fù)雜流場下氣泡界面不穩(wěn)定現(xiàn)象的數(shù)值模擬方法綜述

倪寶玉，張阿漫，張忠宇，李 帥

（哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001）

氣泡界面不穩(wěn)定現(xiàn)象一直是國內(nèi)外關(guān)注和研究的熱點。近年來，邊界元方法（Boundary Element Method，BEM）被廣泛應(yīng)用于氣泡動力學(xué)模擬，文章則側(cè)重綜述邊界元法及其相關(guān)技術(shù)在模擬氣泡界面不穩(wěn)定現(xiàn)象中的應(yīng)用與發(fā)展。首先，簡要回顧不同邊界不同背景流場下氣泡潰滅坍塌的邊界元法研究進(jìn)展；其次，分別回顧和討論氣泡撕裂、融合和在自由面破裂三種典型界面不穩(wěn)定現(xiàn)象的研究進(jìn)展；再次，回顧和討論了氣泡模擬過程中可壓縮性、表面張力和粘性三種因素的影響；最后，在上述基礎(chǔ)上提出了一些尚需進(jìn)一步解決的問題。該文旨在向國內(nèi)同行學(xué)者介紹目前氣泡界面不穩(wěn)定現(xiàn)象的最新研究進(jìn)展，希望能為國內(nèi)學(xué)者研究氣泡動力學(xué)提供一定參考。

氣泡；界面不穩(wěn)定；邊界元法；撕裂；融合

1 引 言

氣泡在流體中的運動學(xué)形態(tài)和動力學(xué)特性一直是流體力學(xué)、環(huán)境工程、醫(yī)學(xué)和船舶與海洋工程等多種領(lǐng)域的研究熱點之一，在自然界中隨處可見流體中的氣泡現(xiàn)象，在各工程領(lǐng)域的應(yīng)用更不勝枚舉。僅以船舶與海洋工程領(lǐng)域為例，較著名的氣泡現(xiàn)象就有高速運轉(zhuǎn)的螺旋槳周圍形成的多種形態(tài)螺旋槳空泡[1]，有水下攻擊武器爆轟后形成的水下爆炸氣泡[2]，也有近年來綠色造船引入的氣泡減阻技術(shù)中的微氣泡或氣穴[3]等。

氣泡作為氣液兩相的自由交界面，在數(shù)學(xué)模型上屬于非線性初值—邊值問題。在數(shù)值模擬上，對于處理氣泡這種帶有自由交界面的問題，可以采用有限元法（FEM）或者邊界元法（BEM）等方法。有限元法求解通?？煞譃閮纱箢惙椒ǎ航缑孀粉櫡ǎ↖nterface Tracking Methods）[4]和界面捕捉法（Interface Capturing Methods）[5]，二者各有利弊。界面追蹤法可以精確地定義流體界面，但是需要非常精密而復(fù)雜的網(wǎng)格重建工作，尤其當(dāng)界面發(fā)生破損或融合等強非線性運動時。此外，計算過程中通常存在質(zhì)量和能量損失。相反，以流體體積法（VOF）[6]、水平集法（Level-Set）[7]以及任意拉格朗日-歐拉法（ALE）[8]為代表的界面捕捉法在流域中弱化界面的概念，而是通過附加流體變量如質(zhì)量分?jǐn)?shù)、空隙度或比熱率等來區(qū)分不同的流體。然而若想精確捕捉界面可能需要很多網(wǎng)格單元，同時在仿真過程中易出現(xiàn)界面模糊效應(yīng)。

不同于有限元法，邊界元法僅對流域的邊界（如氣泡表面）進(jìn)行劃分，而無需對整個流域內(nèi)部及氣泡內(nèi)部進(jìn)行劃分，從而將問題的維度降低一維。BEM的最大優(yōu)點在于高效性，網(wǎng)格和節(jié)點數(shù)的減少使BEM所消耗的時間相對于有限元法至少縮短幾個量級，通常有限元法需要十幾個小時或幾天才能運行完成的結(jié)果應(yīng)用邊界元法幾個小時就可完成。雖然邊界元也有自己的不足，例如通常假設(shè)流體為勢流，對于求解粘性有很大局限性，或無法計及物面內(nèi)外的質(zhì)量和能量交換等，但是由于其避免了重構(gòu)交界面兩側(cè)流體網(wǎng)格的優(yōu)勢，在研究氣泡這種瞬態(tài)自由面大變形問題方面，還是獲得了巨大的應(yīng)用和發(fā)展。

通常而言，氣泡并不會穩(wěn)定地存在于流體中，而是在流體環(huán)境或周圍邊界的復(fù)雜影響下表現(xiàn)出極不穩(wěn)定現(xiàn)象，如劇烈潰滅、撕裂以及融合等。氣泡界面的這種不穩(wěn)定現(xiàn)象一方面會改變氣泡自身的動力學(xué)行為，另一方面也會改變其對周圍結(jié)構(gòu)的載荷特性，故關(guān)于氣泡界面不穩(wěn)定現(xiàn)象的研究一直是國內(nèi)外的熱點之一。本文將回顧并討論氣泡不穩(wěn)定現(xiàn)象的研究進(jìn)展和技術(shù)難點，并主要側(cè)重于邊界元法（BEM）研究氣泡不穩(wěn)定現(xiàn)象的相關(guān)數(shù)值處理和技術(shù)手段。

2 氣泡潰滅坍塌

基于勢流理論的邊界元法（BEM）或稱邊界積分法（BIM）在瞬態(tài)自由交界面上最早應(yīng)用，可能是Longuet-Higgins和Cokelet（1976）[9]。他們應(yīng)用第二格林函數(shù)將拉普拉斯方程轉(zhuǎn)化為邊界積分方程，并采用橢圓積分?jǐn)?shù)值求解積分方程，從而成功地解決了非線性瞬態(tài)水波問題。此后，邊界積分法迅速發(fā)展，在氣泡動力學(xué)上應(yīng)用廣泛（Lenoir（1979）[10]，Guerri等人（1981）[11]，Blake 等人（1986，1987）[12-13]，Chahine（1990）[14]，Best（1993）[15]，Zhang 和 Duncan（1994）[16]，Brujan 等 人 （2002）[17]，Lee，Klaseboer 和Khoo（2007）[18]）。邊界元法求解氣泡動力學(xué)問題的基本思路如下：已知氣泡表面的速度勢，應(yīng)用格林第二函數(shù)求解可獲得氣泡表面的法向速度?φ/?n，同時應(yīng)用速度勢差分法可獲得表面切向速度?φ/?τ，進(jìn)而可合成節(jié)點的全速度▽φ。在時間步進(jìn)中，采用拉格朗日法追蹤氣泡表面節(jié)點，并應(yīng)用拉格朗日坐標(biāo)求解氣泡表面伯努利方程，獲得每一時間步內(nèi)速度勢增量Δφ。在獲得新的速度勢和新的氣泡節(jié)點坐標(biāo)后，即可進(jìn)入循環(huán)，獲得氣泡隨時間的推進(jìn)演化歷程。

在不同邊界或不同背景流場條件下，氣泡在收縮階段可能會發(fā)生不穩(wěn)定坍塌潰滅，并形成高速射流。最典型的兩類邊界是剛性壁面和自由液面。對于剛性壁面附近氣泡的求解中，對于氣泡表面，速度勢已知，滿足狄利克雷（Dirichelet）邊界條件；對于剛性壁面，法向速度已知，滿足紐曼（Neumann）條件，僅需對邊界積分方程得到的系數(shù)矩陣重組即可獲得所需未知量。研究發(fā)現(xiàn)氣泡在壁面附近膨脹階段會輕微遠(yuǎn)離壁面，而收縮階段會被壁面強烈吸引，形成指向壁面的射流。壁面對氣泡的這種效應(yīng)稱為“Bjerknes效應(yīng)”（Bjerknes，1966[19]）。對于自由液面附近的氣泡，自由面和氣泡面均滿足狄利克雷邊界條件，在時間步進(jìn)中也采用拉格朗日法，這樣即可采用與求解氣泡類似的方法追蹤自由面的變形。自由面附近的氣泡動態(tài)研究基本上是與壁面附近氣泡同步展開的，然而與壁面不同，在不考慮浮力影響下，氣泡在坍塌階段將遠(yuǎn)離自由面，并形成指向背離自由面的射流，而自由面將形成明顯的水?，F(xiàn)象（Blake和Gibson，1987[20]）。無論是剛性壁面還是自由液面，氣泡在潰滅后期，射流都將穿透氣泡另一側(cè)壁面形成環(huán)狀氣泡。環(huán)狀氣泡屬于雙連通域，現(xiàn)在應(yīng)用的較廣泛的求解方法有“渦面模型”[17]和“渦環(huán)模型”[21]。關(guān)于邊界附近氣泡不穩(wěn)定潰滅，尤其是環(huán)狀氣泡的研究方法和進(jìn)展，請參見文獻(xiàn)[22]和[23]。

對于不同背景流場，本文選擇船舶與海洋工程中經(jīng)常遭遇的兩類背景流場為例，一是高速螺旋槳誘導(dǎo)的漩渦流場，一是水下爆炸引發(fā)的沖擊波流場。

對于漩渦流場，經(jīng)研究表明，螺旋槳誘導(dǎo)流場中氣泡將隨著渦流做螺旋線運動，在運動過程中自身還有可能撕裂為多個子氣泡（Choi，2006[24]）。最早采用邊界積分法研究漩渦流場中氣泡的變形和運動的當(dāng)屬 Chahine（1983）[25]。Duraiswami和 Chahine（1992）[26]采用漸進(jìn)展開法求解勢流中某線渦附近氣泡的變形，其中假設(shè)氣泡本身的尺度相對于氣泡與線渦的距離是個小量。數(shù)值計算表明氣泡將比正常流中以更快的速度向下游運動，同時將形成指向渦心的斜射流。Hsiao和Pauley（1999）[27]應(yīng)用雷諾平均NS方程模擬翼梢渦流場，并建立球狀氣泡模型研究空化初生現(xiàn)象。之后，Hsiao和Chahine（2003）[28]又進(jìn)一步擴(kuò)展球狀氣泡模型，研究Rankine線渦內(nèi)空化氣泡的尺度效應(yīng)。然而，這種算法的缺點在于無法考慮氣泡與流場的相互作用，也無法計及非球狀效應(yīng)。為了克服這個缺點，Hsiao和Chahine等人（2004，2006）[29-30]基于勢流理論開發(fā)了非球狀模型來研究空化初生、氣泡坍塌潰滅和輻射噪聲等，尤其是研究了氣體擴(kuò)散和溶解濃度的影響，獲得了許多有益的結(jié)論。圖1給出了文獻(xiàn)[31]采用直接計算法（DNS）和文獻(xiàn)[32]采用粒子追蹤法（PTM）計算得到的螺旋槳流場中氣泡的捕獲過程。在國內(nèi)，王獻(xiàn)孚（2009）[33]基于修正的RP方程總結(jié)了螺旋槳流場中球狀空泡初生判別方法。韓寶玉等人（2011）[34]應(yīng)用Hsiao等人開發(fā)的球狀氣泡平均壓力模型（SAP）預(yù)報了不同初始半徑的氣泡的空泡初生數(shù)，并預(yù)報了相應(yīng)的輻射噪聲。然而，國內(nèi)直接應(yīng)用邊界積分法求解渦流場中氣泡變形和撕裂等現(xiàn)象的文獻(xiàn)尚十分缺乏。

圖1 直接計算法（a）[31]和粒子追蹤法（b）[32]計算氣泡在漩渦流場的運動軌跡Fig.1 Bubble’s trajectory in a vortex field calculated by DNS(a)[31]and PTM(b)[32]

對于沖擊波流場，除了重力場中流體的靜壓力梯度和邊界效應(yīng)（壁面或自由面）外，當(dāng)氣泡處于沖擊波流場時，強大的沖擊波壓力也會促使氣泡形成射流，且射流速度可能達(dá)到數(shù)千米每秒的量級（Bourne和Field，1992[35]）。Blake等人（1999）[36]應(yīng)用邊界元法建立了低頻振動的剛性邊界附近氣泡的動態(tài)響應(yīng)，他們在伯努利方程的重力項中附加了一項考慮了邊界正弦振動引起的附加壓力，此項附加壓力是時間的函數(shù)但不是空間的函數(shù)，以此模擬不同振動頻率和振幅強度下氣泡的變形和射流等。Klaseboer等人（2006）[37]借鑒Blake等人的思想，在伯努利方程中引入隨時空均變化的沖擊參考壓力，進(jìn)而應(yīng)用邊界元研究了某一氣泡與矩形脈沖的相互作用。在模擬中假設(shè)可壓縮性僅在氣泡坍塌后期影響明顯，在大部分氣泡運動周期內(nèi)慣性項占主導(dǎo)作用。Klaseboer等人（2007）[38]又進(jìn)一步考慮了某一碎石沖擊波與氣泡的相互作用，在數(shù)值模擬中考慮了球狀氣泡初始形態(tài)（如膨脹、靜止或收縮）對于射流速度等影響。國內(nèi)，倪寶玉等人（2011）[39]在Klaseboer等人的基礎(chǔ)上，考慮爆炸沖擊波在船底和水底形成的反射波對水下爆炸氣泡的影響，將邊界元理論計算應(yīng)用到工程領(lǐng)域。圖2即給出了文獻(xiàn)[39]采用邊界元法（BEM）與文獻(xiàn)[40]采用自由拉格朗日法（FLM）數(shù)值模擬沖擊波流場中氣泡坍塌射流的演化過程。經(jīng)過計算發(fā)現(xiàn)，反射波對氣泡射流載荷的影響較弱，在工程計算中可將沖擊波載荷和氣泡射流載荷分別計算，無需考慮二者間的耦合作用。當(dāng)水下武器多發(fā)攻擊目標(biāo)時，異相爆炸間也會存在沖擊波與氣泡載荷的相互作用，尚需進(jìn)一步研究。

圖2 自由拉格朗日法（上）[40]和邊界元法（下）[39]計算沖擊波下氣泡坍塌對比圖Fig.2 Contrast of bubble collapse under shockwave calculated by FLM(top)[40]and BEM(bottom)[39]

3 氣泡撕裂和融合

除了上文闡述的潰滅現(xiàn)象和環(huán)狀氣泡現(xiàn)象，在日常生活中氣泡不穩(wěn)定界面中更為常見的有氣泡的撕裂、融合，以及自由面處破裂等。本章則在上述的基礎(chǔ)上進(jìn)一步分析氣泡撕裂和融合的動力學(xué)行為，尤其關(guān)注邊界元法在處理撕裂和融合問題的處理手段。

3.1 撕裂

氣泡撕裂是自然界中常見現(xiàn)象之一，導(dǎo)致氣泡撕裂的真實物理原因可能十分復(fù)雜，這里總結(jié)以下幾個方面：首先，較容易理解的穩(wěn)態(tài)氣泡遭受外力或外物干擾，Kim和Kim（1990）[41]應(yīng)用不同截面形狀的物體主動切割某一穩(wěn)態(tài)上浮氣泡，總結(jié)浮體的形狀、速度等因素對氣泡撕裂的影響；其次，較常見的是氣體體積增大，當(dāng)大于某一臨界值后，氣泡在微小擾動下即自發(fā)破裂，Wilkinson等人（1993）[42]研究了不同氣體密度和不同液體特性對于氣泡碎裂的影響；再次，在湍流漩渦作用下，由于粘性中心半徑、速度環(huán)量和壓力等在渦流場內(nèi)變化劇烈，氣泡十分容易撕裂，Martinez-Bazan等人（1999）[43]即研究了湍流場中氣泡的運動與撕裂現(xiàn)象；最后，不同邊界的作用也會導(dǎo)致氣泡撕裂，例如在兩個壁面的Bjerknes力的相互吸引下，一個氣泡在收縮階段即可能分裂為兩個子氣泡。

為了揭示氣泡撕裂現(xiàn)象及其載荷特性，人們進(jìn)行了大量的研究。在實驗方面，Ishida等（2001）[44]應(yīng)用水箱模擬氣泡在兩水平壁面間的對稱撕裂與非對稱撕裂，并通過壓力傳感器測得壁面與流場的壓力變化；Choi和Chahine（2003）[45]應(yīng)用高速攝影拍攝了氣泡在兩傾斜壁面和兩垂直壁面間的撕裂過程。圖3給出了采用電火花氣泡模擬兩變化距離剛性壁面內(nèi)氣泡的撕裂現(xiàn)象。

圖3 兩平行壁面間氣泡撕裂現(xiàn)象實驗Fig.3 Experiment on breakup of a bubble between two horizontal walls

數(shù)值模擬上，早期 Ellis（1956）[46]，Lauterborn 和 Hentschel（1985）[47]等人都曾研究過兩剛性邊界間單個氣泡的運動，預(yù)測氣泡收縮后期可能形成沙漏狀；Mohammad（2006）[48]等系統(tǒng)地研究了軸對稱氣泡在兩平行壁面間的頸縮狀態(tài)，對氣泡撕裂前的行為做了詳細(xì)的工作，但沒有模擬出氣泡撕裂及以后的過程；Choi和Chahine等（2003）[45]提出了軸對稱氣泡的撕裂準(zhǔn)則，建立了軸對稱撕裂模型并模擬了對稱撕裂與非對稱撕裂；張阿漫和倪寶玉等人（2010）[49]在Choi和Chahine研究的基礎(chǔ)上提出了三維氣泡撕裂準(zhǔn)則，研究了兩剛性壁面間單個氣泡的對稱和非對稱撕裂現(xiàn)象，給出了不同工況下氣泡形態(tài)、射流速度、開爾文沖量等物理量的變化。

圖4 氣泡撕裂數(shù)值過程示意圖（來自于文獻(xiàn)[49]）Fig.4 Process of the split of a 3D bubble(data from Ref.[49])

不直接模擬氣泡撕裂的真實過程，在數(shù)值上人工干預(yù)，重新分布滿足分裂準(zhǔn)則的氣泡節(jié)點，使得一個氣泡分裂為兩個子氣泡，圖4給出了邊界元法處理氣泡撕裂的過程，具體實行過程如下：

（a）確定局部“對稱軸”：如圖4（a）所示，氣泡在撕裂前期會發(fā)生較明顯的“頸縮”現(xiàn)象，在局部呈現(xiàn)筒狀結(jié)構(gòu)，選擇筒狀結(jié)構(gòu)的中心線為局部“對稱軸”。

（b）計算最小距離；在每個時間步內(nèi)，計算頸縮區(qū)節(jié)點到局部“對稱軸”的距離，并將最小距離設(shè)定為Δs。當(dāng)無量綱距離Δs≤0.03后，認(rèn)為滿足撕裂準(zhǔn)則。否則繼續(xù)下一步計算，直至滿足撕裂準(zhǔn)則為止，這個過程通過程序自動完成，無需人工干涉。

（c）確定“破裂線”；當(dāng)Δs滿足撕裂準(zhǔn)則后，設(shè)定對應(yīng)此Δs的節(jié)點編號為I1，與I1在同一圈的節(jié)點I2，…，Im則被定義為“破裂線”?！捌屏丫€”上所有節(jié)點將全部移除，并將I1，I2，…，Im節(jié)點坐標(biāo)的算術(shù)平均值作為新的一點J的空間坐標(biāo)。

（d）撕裂中心點：將中心點J沿局部“對稱軸”軸向上下各移動微小距離0.001，得到兩個新節(jié)點J1和J2，并將原來斷開的節(jié)點分別與J1和J2連接上，重新編號節(jié)點號和單元號，從而實現(xiàn)了一個氣泡到兩個子氣泡的分裂，如圖4（b）所示。

（e）計算新節(jié)點物理量；J1和J2點的速度勢均取“破裂線”上點I1，I2，…，Im的算術(shù)平均值即可更新子氣泡的速度勢和表面位置。

3.2 融合

與撕裂相對，融合現(xiàn)象是氣泡動力學(xué)研究的另一熱點，在許多領(lǐng)域如氣-液/泥漿反應(yīng)堆、氧化/氫化反應(yīng)、發(fā)酵、污水處理以及最近的細(xì)胞工程等均有廣泛應(yīng)用。通常而言，對于氣泡融合可以化為三個階段（Rungsiyaphornrat等人，2003[50]）：第一階段，子氣泡相互吸引，相對運動至較接近的位置，足以保證一段時間的演化后能達(dá)到融合條件；第二階段，兩個氣泡在繼續(xù)吸引下表面開始變平坦，其間的液體逐漸被擠壓出去，很長一段時間內(nèi)液膜越來越薄；當(dāng)然，如果此時氣泡在過大的碰撞力下分開而反向運動的話，則不會發(fā)生融合。第三階段，當(dāng)液膜變得很薄，進(jìn)入某一臨界值后，分子力等內(nèi)力加之微小外界擾動則將促使薄膜破裂，氣泡融為一體。Chen等人（2005）[51]給出了氣泡融合三個階段的示意圖，如圖5所示。

圖5 氣泡融合示意圖（圖形改自于文獻(xiàn)[51]）Fig.5 The sketch of bubble coalescence(figure revised from reference[51])

與液滴一樣，在兩個氣泡融合前的小段時間，其間形成的液膜厚度非常小，可達(dá)到微米數(shù)量級。處理液滴融合的薄膜問題時，人們常采用流體潤滑理論模型（Chevaillier等人，2006[52]），這樣可以很好地解決液體薄膜的計算。在模擬氣泡融合時，起初也試圖引入流體潤滑理論模型，但經(jīng)過對兩個氣泡的數(shù)值模擬，發(fā)現(xiàn)氣泡融合前的液膜運動規(guī)律并不能用流體潤滑理論控制。這是因為液膜的雷諾數(shù)[50]至少在100以上，有的甚至達(dá)到1萬以上，遠(yuǎn)遠(yuǎn)不能滿足流體潤滑理論應(yīng)用時雷諾數(shù)小于1的要求[53]。關(guān)于氣泡融合的實驗數(shù)據(jù)相對較少，Klaseboer等人（2000）[53]實驗研究了勻速運動接近的兩氣泡的融合現(xiàn)象，而后，Rungsiyaphornrat等人（2003）[50]應(yīng)用軸對稱勢流理論模型求解了兩個水下爆炸氣泡同相和異相融合，數(shù)值結(jié)果與Klaseboer等人的實驗進(jìn)行對比，吻合度良好，發(fā)現(xiàn)同相氣泡融合后將形成相互吸引的對射流。張阿漫等人（2010）[54]在Rungsiyaphornrat等人的基礎(chǔ)上，采用邊界元法數(shù)值建立了同相氣泡三維融合模型，數(shù)值模擬了兩同相爆炸氣泡水中融合現(xiàn)象，并對比融合氣泡和單個氣泡引起的滯后流和壓力的差異。然而，上述模型均計算到對射流相互接觸時刻為止，沒有考慮對射流碰撞后的環(huán)狀氣泡和環(huán)狀射流階段。

圖6 氣泡融合數(shù)值過程示意圖Fig.6 Process of the mergence of two 3D bubbles

與撕裂類似，并不真實模擬氣泡融合過程，僅在數(shù)值計算過程中，當(dāng)氣泡網(wǎng)格節(jié)點距離滿足融合準(zhǔn)則后，進(jìn)行人工干預(yù)，重新組合網(wǎng)格而形成一個新的氣泡。圖6給出了邊界元法處理融合過程示意圖，具體實行過程如下：

（a）計算最小距離：設(shè)氣泡表面節(jié)點個數(shù)分別為N1和N2，則每個時間步內(nèi)均計算兩個氣泡各節(jié)點間的距離，得到N1×N2個距離，并假定其中最小距離為Δs，對應(yīng)的兩個節(jié)點分別標(biāo)注為I和I′，如圖6（a）所示。當(dāng)無量綱距離Δs≤0.03后，則認(rèn)為滿足融合準(zhǔn)則，否則繼續(xù)下一步計算，直到滿足融合準(zhǔn)則為止，這個過程通過程序自動完成，無需人工干涉。

（b）確定“融合線”：當(dāng)獲得了滿足融合準(zhǔn)則的節(jié)點I和I′后，以節(jié)點I為例，將與I相連接的第一圈節(jié)點I1，I2，…，Im定義為“融合線”。在實際計算中，如果網(wǎng)格尺寸較小，也可以將融合線擴(kuò)展到I周圍的第二圈第三圈節(jié)點。類似地，在I′附近也會找到對應(yīng)的一條“融合線”。

（c） “融合線”縫合：將融合線上的節(jié)點 I1，I2，…，Im和 I1′，I2′，…，In′人為去除，并取對應(yīng)節(jié)點坐標(biāo)的算數(shù)平均值而得到新的一圈節(jié)點坐標(biāo)J1，J2，…，Jk，其中k=max m,（）n 。假若m≠n，需在認(rèn)為去除融合線上的節(jié)點之前，先自身插值，確保在取平均值時，所有節(jié)點均可尋找對應(yīng)的節(jié)點。

（d）重新編號：分別將J1，J2，…，Jk與原來分離的節(jié)點連接上，并將節(jié)點號和單元號重新編號，至此兩個子氣泡就融合為一個氣泡了，如圖6（b）所示。

（e）計算新節(jié)點物理量；對于J1，J2，…，Jk點的速度勢均取移除點的算術(shù)平均值即可更新子氣泡的速度勢和表面位置。

3.3 自由面處破裂

上浮氣泡或接近自由面附近水下爆炸氣泡均會產(chǎn)生在自由面的破裂現(xiàn)象。對于上浮氣泡而言，實驗和數(shù)值均發(fā)現(xiàn)氣泡破裂后，會形成向上的小水滴，小水滴對氣液兩相間的質(zhì)量交換起到重要作用，在氣泡體積譜分析、通風(fēng)、沸騰、去氣、蒸餾以及蒸發(fā)方面均有重要應(yīng)用。與氣泡融合現(xiàn)象類似，通常認(rèn)為氣泡破裂的全部過程也可劃分為三個階段（Georgescu等人，2002[55]）：第一階段，接近階段。此階段氣泡加速上浮，越來越接近自由液面，同時自由面微微隆起。有趣的是，實驗觀察（Malysa等人，2005[56]）到氣泡不是立即與自由面融合，而是在自由面附近微幅彈跳，直至振幅衰減，氣泡在自由面下處于平衡位置；第二階段，融合階段。此階段氣泡上表面與自由面間形成的水膜逐漸變薄，直至厚度小于某個臨界值，則水膜在微小擾動下將破碎；第三階段：坍塌階段。此階段氣泡內(nèi)部將和大氣連通，氣泡將坍塌。此時，表面波會匯聚于氣泡坍塌軸心并產(chǎn)生一股向上的高速射流。并形成指向上方的不穩(wěn)定射流。圖7即給出了氣泡在自由面處破裂的示意圖。在這三個階段中，將可能形成兩種類型的水滴[53]：其一發(fā)生在第二階段水膜破碎時刻，破裂的水膜碎片將四處飛濺，形成所謂的“水膜液滴”；其二發(fā)生在坍塌階段末期，此時指向上方的高速射流是不穩(wěn)定的，在射流中上部可能發(fā)生連續(xù)的破碎，并產(chǎn)生若干“射流水滴”。射流水滴可能攜帶較大的能量，故學(xué)者們相對更關(guān)心射流水滴的變化。

圖7 氣泡在自由面處破裂示意圖Fig.7 The sketch of bubble bursting at free surface

由于射流水滴可能攜帶較高的能量，本文主要關(guān)注“射流水滴”的相關(guān)進(jìn)展。在實驗上，Kientzler等人（1954）[57]即開始應(yīng)用高速攝影儀觀察氣泡在自由面的破裂現(xiàn)象，并指出第一滴射流水滴的體積約為原始?xì)馀蒹w積的十分之一。Garner等人（1954）[58]描述了蒸餾水和酒精內(nèi)射流水滴的速度和數(shù)量，指出原始?xì)馀蒹w積由幾十微米增加到幾毫米時，第一滴射流水滴的速度將由幾十米每秒迅速降到幾米每秒，同時形成的射流水滴的個數(shù)將減少，最多的射流水滴數(shù)目也不會超過10個。Resch和Afeti（1991，1992）[59-60]實驗得到，對于海水而言，當(dāng)原始?xì)馀莅霃酱笥? mm后，則無法形成射流水滴現(xiàn)象。在數(shù)值模擬上，由于自由面邊界條件的強非線性，氣泡在自由面處的破裂現(xiàn)象一直是數(shù)值建模方面的一大難點。最早的數(shù)值模擬研究應(yīng)屬Boulton-Stone和Blake（1993，1995）[61-62]，他們應(yīng)用邊界元法數(shù)值模擬氣泡破裂過程，并引入邊界層考慮弱粘性效應(yīng)。研究表明，表面張力對于氣泡在自由面處破裂起到十分重要的作用，而粘性對于前期氣泡運動影響很弱，僅在射流形成后起到較小的耗散效應(yīng)，阻礙射流水滴的形成。Sussman和Smereka（1997）[63]采用Level set法求解NS方程，數(shù)值模擬初始時刻軸對稱球狀氣泡在水平自由面處破裂，但氣泡半徑均選取大于4 mm的工況，與之前的0.5-3 mm小氣泡沒有對比性。Georgescu等人（2002）[55]采用帶有粘性效應(yīng)修正的邊界元法模擬氣泡破裂過程的射流水滴現(xiàn)象，在時間步進(jìn)上采用二階泰勒展開。數(shù)值結(jié)果與實驗值進(jìn)行對比，表明此方法在模擬破裂前期精度很好，在射流水滴形成的后期誤差有所增大，故得到結(jié)論是僅在整個氣泡演化過程的后期NS方程才可能需要采用。

4 幾類主要影響因素

氣泡動力學(xué)進(jìn)展中，學(xué)者們十分關(guān)心可壓縮性、表面張力和粘性等對于氣泡不穩(wěn)定潰滅坍塌，撕裂融合以及相關(guān)載荷的影響規(guī)律，并試圖揭示在什么樣的條件下這些因素起主要作用必須考慮，或起微弱作用可以忽略。故在本節(jié)中，分別闡述國內(nèi)外對這幾類影響因素的研究進(jìn)展。

4.1 可壓縮性

Prosperetti和Lezzi（1986）[64]應(yīng)用漸近匹配展開法對氣泡壁馬赫數(shù)展開到一階項，從而研究在可壓縮流體中氣泡的徑向運動。他們獲得了描述弱可壓縮流體內(nèi)球狀氣泡運動的單一參數(shù)方程組，此方程組包含了著名的 Herring-Keller方程（Herring，1941[65]；Keller和 Kolodner，1956[66]）。 Prosperetti和 Lezzi（1987）[67]又進(jìn)一步將馬赫數(shù)展開到二階項，并指出對于不可壓縮流體的一階項修正已經(jīng)在很大程度上考慮了可壓縮效應(yīng)，二階項僅有較弱的影響。Pozrikidis（2001）[68]采用邊界積分法數(shù)值模擬某無粘二維可壓縮氣泡在Stokes流中的運動，發(fā)現(xiàn)在純剪切流或純拉伸流中，初始圓環(huán)狀的氣泡在收縮后最終均將達(dá)到一穩(wěn)定態(tài)，且穩(wěn)定態(tài)的面積比初始圓環(huán)狀面積大。Darren和Crowdy（2003）[69]擴(kuò)展了Pozrikidis的邊界積分法，將原本的積分方程轉(zhuǎn)化為兩個耦合的非線性常微分方程組，此耦合方程組的計算更直接更穩(wěn)定，且需要更小的計算機(jī)內(nèi)存。數(shù)值模擬還發(fā)現(xiàn)絕熱氣泡較等溫氣泡更容易膨脹。Wang和Blake（2010，2011）[70-71]也采用漸近匹配法對氣泡壁馬赫數(shù)分別進(jìn)行了一階和二階展開，研究可壓縮流場中非球狀氣泡的運動。與Prosperetti和Lezzi不同的是，Wang和Blake將流場劃分為內(nèi)外兩部分，對于接近氣泡的內(nèi)流域認(rèn)為一階和二階均滿足不可壓縮的拉普拉斯方程，對于遠(yuǎn)離氣泡的外流域直接采用一階和二階線性波動方程求解，在內(nèi)外域交接面進(jìn)行展開匹配。通過混合拉格朗日—歐拉法和修正邊界元法求解氣泡的變形，此方法的主要優(yōu)點在于保證精度的同時提高了全局求解波動方程的效率。

Geers和Hunter（2002）[72]采用雙漸近法研究了可壓縮流場中水下爆炸氣泡的球狀運動，給出了氣泡運動的經(jīng)驗公式。王詩平（2011）[73]在Geers和Hunter模型的基礎(chǔ)上，采用雙漸近法求解了可壓縮流場非球狀氣泡的運動。Geers和Hunter采用頻域內(nèi)局部和全局近似，高低頻匹配后再通過拉普拉斯變換到時域問題，而王詩平（2011）[73]采用時域內(nèi)前期逼近和后期逼近，在時域內(nèi)匹配求解，在數(shù)值實現(xiàn)上更為方便。在此基礎(chǔ)上，王詩平研究了自由場和邊界附近可壓縮氣泡的運動，并指出對于水下爆炸氣泡第一周期的運動，可壓縮性影響較小，工程應(yīng)用中可以忽略。

4.2 表面張力

Beeching（1942）[74]首次在Rayleigh球狀氣泡運動中考慮了表面張力2σ/R的影響。此后，Akhatov等人（1997）[75]數(shù)值研究了表面張力和氣體擴(kuò)散對于聲致發(fā)光氣泡穩(wěn)定性的影響，指出表面張力對于小氣泡運動十分重要，可能是聲致發(fā)光氣泡穩(wěn)定脈動的主要原因之一。在數(shù)值模擬上，表面張力對于氣泡運動影響的難點在于局部曲率的求解，因為對于給定的流體和氣體，表面張力系數(shù)是已知的。對于二維或者軸對稱氣泡而言，曲率有解析公式，可參見Bronshtein和Semendyayev（1997）[76]；對于三維氣泡，曲率的數(shù)值求解是有一定難度的。Brackbill等人（1991）[77]開發(fā)連續(xù)譜法求解曲率，核心思想為將氣液交界面認(rèn)為是有限厚度的過渡區(qū)域，在過渡區(qū)域內(nèi)定義與局部曲率成比例的力密度函數(shù)。連續(xù)譜法消除了對交界面重構(gòu)的需要，簡化了表面張力的計算方法，可精確求解二維或三維曲率，而無需受限于模型的網(wǎng)格數(shù)目、復(fù)雜度以及動態(tài)演化等因素。Zinchenko等人（1997）[78]闡述了一種更簡便的求三角形單元表面曲率和法向量的迭代算法，該算法是在Rallison（1981）[79]思想上的進(jìn)展，主要是通過最小二乘法擬合局部最優(yōu)拋物面，數(shù)值計算表明該算法魯棒性較強。Georgescu等人（2002）[55]實驗研究了氣泡在不同表面張力液體自由表面處的破裂，發(fā)現(xiàn)表面張力大的液體形成的射流水滴體積大于表面張力小的液體，且原始?xì)馀輰?yīng)的臨界半徑（可形成射流水滴的最大氣泡半徑）也越大。

在國內(nèi)，黃繼湯等人（1996）[80]研究了表面張力對于單個氣泡運動特性的影響，研究表明表面張力在氣泡膨脹階段起到延緩作用，在收縮階段則起到了加速作用。并通過高速攝影實驗對五種表面張力不同的液體內(nèi)氣泡運動進(jìn)行拍攝分析，分析得到表面張力越小的液體中，氣泡壁面速度達(dá)到的最大值越大。王起棣和張慧生等人（2003，2005）[81-82]分別研究了表面張力對于剛性壁面附近單個和兩個氣泡運動的影響，數(shù)值模擬采用軸對稱邊界元法，結(jié)果表明表面張力系數(shù)越大的流體中氣泡潰滅速度越快，射流隨之變寬變短。而對于兩個氣泡情況而言，表面張力的影響與壁面和上氣泡對于下氣泡的Bjerknes力的大小有關(guān)。

4.3 粘性

關(guān)于氣泡粘性效應(yīng)的考慮，最直接的做法是求解NS方程，然而由于數(shù)值直接求解NS方程的難度，學(xué)者們通常從高雷諾數(shù)和低雷諾數(shù)兩方面入手簡化問題。對于高雷諾數(shù)的情況，最經(jīng)典的當(dāng)屬Miksis等人（1982）[83]引入的氣泡邊界層理論，Miksis等人認(rèn)為氣泡周圍存在一薄邊界層，邊界層外流體滿足無粘不可壓縮勢流理論，邊界層內(nèi)考慮流體粘性效應(yīng)，邊界層兩側(cè)滿足法向應(yīng)力連續(xù)，但不滿足切向應(yīng)力連續(xù)。由于此邊界層很薄，忽略邊界層內(nèi)外速度的變化，這樣便可采用邊界積分法數(shù)值模擬計及弱粘性效應(yīng)的氣泡動態(tài)變化。此后，研究人員從不同角度改進(jìn)Miksis等人的氣泡邊界層理論。Lundgren和Mansour（1988）[84]進(jìn)一步分析邊界層方程，獲得了邊界層兩側(cè)壓力差和法向速度差的表達(dá)式。Boulton-Stone和Blake（1993，1995）[61-62]又進(jìn)一步擴(kuò)展了Lundgren和Mansour的方法，使得邊界層兩側(cè)速度的切向分量也連續(xù)，在此基礎(chǔ)上首次成功地數(shù)值模擬了單個氣泡在自由面處的破裂現(xiàn)象。Georgescu等人（2002）[55]也采用Miksis的氣泡邊界層理論模擬氣泡在自由面處的破裂現(xiàn)象，在邊界積分和時間步進(jìn)上均采用二階單元，以獲得更好的精度。Klaseboer等人（2011）[85]采用Joseph和Wang（2004）[86]的思想修正邊界層理論，同時考慮了法向和切向應(yīng)力的連續(xù)，應(yīng)用邊界積分法數(shù)值模擬了高雷諾數(shù)下氣泡在流體中的脈動上浮現(xiàn)象。

對于低雷諾數(shù)問題，當(dāng)雷諾數(shù)很低時，流場可簡化為Stokes流。Youngren和Acrivos（1976）[87]首先應(yīng)用邊界積分方程求解了Stokes流內(nèi)單個對稱氣泡的變形，對于表面張力較大的氣泡形狀，數(shù)值計算結(jié)果與Barthes-Biesel和Acrivosde（1973）[88]的解析解吻合良好；對于表面張力較小的氣泡形狀，與Buckmaster（1972）[89]應(yīng)用細(xì)長體理論獲得的表達(dá)式基本一致。 近期，Eshpuniyani等人（2008，2010）[90-91]采用邊界元法研究了Stokes流中貼附于壁面上的微氣泡的變形和滑移，旨在為強粘性微血管循環(huán)流中氣泡的運動形態(tài)提供參考。

目前在國內(nèi)，關(guān)于考慮粘性的氣泡運動，如粘性流中氣泡上升問題等，多采用直接數(shù)值模擬法如Front Tracking（陳斌等人，2005[92]），Level Set（李彥鵬等人，2007[93]）等，很少有對邊界積分法進(jìn)行粘性修正的。關(guān)于邊界積分法考慮粘性問題還有待于進(jìn)一步研究。

5 需進(jìn)一步研究的問題

國內(nèi)外關(guān)于氣泡的不穩(wěn)定界面現(xiàn)象的研究日益豐富，實驗技術(shù)和數(shù)值模擬手段相應(yīng)地完善成熟。尤其對于邊界元法，在處理氣泡這類具有大變形氣液交界面的問題上具有較大的優(yōu)勢，獲得了研究人員的青睞，取得了較多有益的進(jìn)展與結(jié)論。不過也存在以下幾方面不足：

（1）目前建立的氣泡邊界元模型中，通常忽略了粘性效應(yīng)的影響，而在氣泡減阻技術(shù)或氣泡在漩渦流場中的運動等情況，粘性效應(yīng)可能會很大地影響氣泡的運動，并改變氣泡的載荷形式，關(guān)于如何在邊界元模型中考慮粘性修正，仍需要加強研究；

（2）目前關(guān)于氣泡的撕裂和融合現(xiàn)象主要集中在軸對稱模型，關(guān)于三維氣泡的撕裂和融合數(shù)值模型和準(zhǔn)則的研究相對薄弱，尚需加強。對于融合問題，同向融合研究成果較多，異相融合的相關(guān)數(shù)值模擬尚需深入；

（3）氣泡和自由面的相互作用研究比較豐富，但關(guān)于水下爆炸氣泡在接近自由面附近的運動，包括氣泡在自由面的破裂和自由面自身的飛濺等研究尚需加強；

（4）目前業(yè)內(nèi)基本已公認(rèn)氣泡射流是造成附近結(jié)構(gòu)毀傷的主要載荷之一，然而氣泡射流對結(jié)構(gòu)的毀傷模式，以及不同爆距等參數(shù)如何影響氣泡射流對結(jié)構(gòu)的毀傷等問題，仍不十分明朗，需要進(jìn)一步深入研究。

[1]Hsiao C T,Pauley L L.Numerical computation of the tip vortex flow generated by a marine propeller[J].Journal of Fluids Engineering,1999,121(3):638-645.

[2]van Wijingaarden L.On the collective collapse of a large number of gas bubbles in water[C]//Proceeding of 11th International Congress on Applied Mechanics.Munich,1964:854-861.

[3]Li Y B,Wu X Y,Ma Y,Wang J G.A method based on potential theory for calculating air cavity formation of an air cavity resistance reduction ship[J].Journal of Marine Science Application,2008,7:98-101.

[4]Unverdi S O,Tryggvason G.A front-tracking method for viscous,incompressible,multi-fluid flows[J].Journal of Computational Physics,1992,100:25-37.

[5]Bonometti T,Magnaudet J.An interface-capturing method for incompressible two-phase flows.Validation and application to bubble dynamics[J].International Journal of Multiphase Flow,2007,33:109-133.

[6]Hirt C W,Nichols B D.Volume of Fluid(VOF)method for the dynamics of free boundaries[J].Journal of Computational Physics,1981,39:201-225.

[7]Huang J T,Zhang H S.Level set method for numerical simulation of a cavitation bubble,its growth,collapse and rebound near a rigid wall[J].Acta Mechanica Sinica,2007,23:645-653.

[8]Ding Z,Gracewsi S M.The behaviour of a gas caviyy impacted by a weak or strong shock wave[J].Journal of Fluid Mechanics,1996,309:183-209.

[9]Longuet-Higgins M S,Cokelet E D.The deformation of steep surface waves on water.I.A numerical method of computation[J].Philosophical Transactions of the Royal Society of London,1976,A 350:1-26.

[10]Lenoir M.A calculation of the parameters of the high-speed jet formed in the collapse of a bubble[J].Journal of Applied Mechanical Technology Physics,1979,20(3):333-337.

[11]Guerri L,Lucca G,Prosperetti A.A numerical method for the dynamics of nonspherical cavitation bubbles[C]//In Proceedings of the Second International Colloquium on Drops and Bubbles.California,USA,1981:175.

[12]Blake J R,Taib B B,Doherty G.Transient cavities near boundaries.Part I:Rigid boundary[J].Journal of Fluid Mechanics,1986,170:479-497.

[13]Blake J R,Taib B B,Doherty G.Transient cavities near boundaries.Part 2:Free surface[J].Journal of Fluid Mechanics,1987,181:197-212.

[14]Chahine G L,Numerical modeling of the dynamic behavior of bubble in nonuniform flow field[C]//In ASME symposium on numerical methods for multiphase flows.Toronto,1990.

[15]Best J P.The formulation of toroidal bubbles upon collapse of transient cavities[J].Journal of Fluid Mechanics,1993,251:79-107.

[16]Zhang S,Duncan J H.On the nonspherical collapse and rebound of a cavitation bubble[J].Physics of Fluids,1994,6(7):2352-2362.

[17]Brujan E A,Keen G S,Vogel A,et al.The final stage of the collapse of a cavitation bubble close to a rigid boundary[J].Physical of Fluids,2002,14(1):85.

[18]Lee M,Klaseboer E,Khoo B C.On the boundary integral method for the rebounding bubble[J].Journal of Fluid Mechanics,2007,570:407-429.

[19]Bjerknes.Fields of force[M].Columbia:Columbia University Press,1966:45-47.

[20]Blake J R,Gibson D C.Cavitation bubbles near boundaries[J].Journal of Fluid Mechanics,1987,19:99-123.

[21]Wang Q X,Yeo K S,Khoo B C,Lam K Y.Vortex ring modelling of toroidal bubbles[J].Theoretical Computational Fluid Dynamics,2005,00:1-15.

[22]張阿漫,姚熊亮.基于邊界積分法的氣泡動態(tài)特性綜述[J].力學(xué)進(jìn)展,2008,38(5):561-570.

[23]張阿漫,汪 玉,聞學(xué)友,倪寶玉,等.水下爆炸氣泡動態(tài)特性研究綜述[J].船舶力學(xué),2009,13(5):828-840.

Zhang Aman,Wang Yu,Wen Xueyou,Ni Baoyu,et al.Review of the dynamics of the underwater explosion bubble[J].Journal of Ship Mechanics,2009,13(5):828-840.

[24]Choi J H.Dynamics and noise emission of vortex cavitation bubbles[R].University of Michigan,USA,2006.

[25]Chahine G L,Genoux Ph F.Collapse of a cavitating vortex ring[J].Trans.ASME,Journal of Fluids Engineering,1983,105:400.

[26]Duraiswami R,Chahine G L.Analytical study of a gas bubble in the flow field of a line vortex[J].Cavitation and Multiflow Forum,ASME,1992.

[27]Hsiao C T,Pauley L L.Study of tip vortex cavitation inception using Navier-Stokes computation and bubble dynamics model[J].Journal of Fluids Engineering,1999,l21:198-204.

[28]Hsiao C T,Chahine G L.Scaling effect on prediction of cavitation inception in a line vortex flow[J].Journal of Fluid Engineering,2003,125:53-60.

[29]Hsiao C T,Chahine G L.Prediction of tip vortex cavitation inception using coupled spherical and nonspherical bubble models and Navier-Stokes computations[J].Journal of Marine Science Technology,2004,8:99-108.

[30]Hsiao C T,Jain A,Chahine G L.Effect of gas diffusion on bubble entrainment and dynamics around a propeller[C]//26th Symposium on Naval Hydrodynamics,September,17-22,2006.Rome,Italy,2006.

[31]Oweis G F,Hout I E,Iyer C,Tryggvason G,Ceccioa S L.Capture and inception of bubbles near line vortices[J].Physics of Fluids,2005,17:1-14.

[32]Ni B Y,Zhang A M,Yao X L,Wang B.Numerical simulation of trajectory and deformation of bubble in tip vortex[J].Applied Mathematics and Mechanics(English edition),2012,33(6):1-16.

[33]王獻(xiàn)孚.空化泡和超空化泡流動理論及應(yīng)用[M].北京:國防工業(yè)出版社,2009.

[34]韓寶玉,熊 鷹,葉金銘.梢渦空泡初生預(yù)報方法[J].海軍工程大學(xué)學(xué)報,2011,23(4):27-32.

[35]Bourne N K,Field J E.Shock-induced collapse of single cavities in liquids[J].Journal of Fluid Mechanics,1992,244:225-240.

[36]Blake J R,Keen G S,Tong R P,Wilson M.Acoustic cavitation:the fluid dynamics of non-spherical bubbles[J].Philosophy Transaction of Royal Society,London,1999,A 357:251-267.

[37]Klaseboer E,Turangan C K,Fong S W,et al.Simulations of pressure pulse-bubble interaction using boundary element method[J].Computer Methods in Applied Mechanics and Engineering,2006,195:4287-4302.

[38]Klaseboer E,Fong S W,Turangan C,et al.Interaction of lithotripter shockwaves with single inertial cavitation bubbles[J].Journal of Fluid Mechanics,2007,593:33-56.

[39]倪寶玉,張阿漫,朱 楓,吳巧瑞.艦船近場爆炸反射波對氣泡射流特性影響[J].哈爾濱工程大學(xué)學(xué)報,2011,32(8):976-983.

[40]Jamaluddin A R,Ball G J,Leighton T G.Free-lagrange simulations of shock/bubble interaction in shock wave lithotripsy[C].Proceedings of the Second International Conference on Computational Fluid Dynamics,2002.

[41]Kim J O,Kim S D.Bubble breakage phenomena,phase holdups and mass transfer in three-phase fluidized beds with floating bubble breakers[J].Chemical Engineering Process,1990,28:101-111.

[42]Wilkinson P M,Schayk A V,Spronken J P M.The influence of gas density and liquid properties on bubble breakup[J].Chemical Engineering Science,1993,48(7):1213-1226.

[43]Martinez-Bazan C,Montanes J L,Lasheras J C.On the breakup of an air bubble injected into a fully developed turbulent flow.Part 1.Breakup frequency[J].Journal of Fluid Mechanics,1999,401:157-182.

[44]Ishida H,Nuntadusit C,Kimoto H.Cavitation bubble behavior near solid boundaries[J].CAV,2001,A5-003.

[45]Choi J K,Chahine G L.Noise due to extreme bubble deformation near inception of tip vortex cavitation[J].FEDSM,2003:45313.

[46]Ellis A T.Techniques for pressure pulse measurements and high-speed photography in ultrasonic cavitation,cavitation in hydrodynamics[C].HMSO,London,Paper 8,1956.

[47]Lauterborn W,Hentschel W.Cavitation bubble dynamics studied by high speed photography and holography.Part I:Ultrasonics[M].In press,1985:258-293.

[48]Mohammad T S,Amir A,Mohammad R S.Numerical study on the dynamics of an electrical discharge generated bubble in EDM[J].Engineering Analysis with Boundary Elements,2006,30:503-514.

[49]張阿漫,倪寶玉,宋炳月,姚熊亮.狹窄流域內(nèi)氣泡破裂現(xiàn)象數(shù)值模擬[J].應(yīng)用數(shù)學(xué)和力學(xué),2010,31(4):420-432.

[50]Rungsiyaphornrat S,Klaseboer E,Khoo B C,Yeo K S.The merging of two gaseous bubbles with an application to underwater explosions[J].Computers&Fluids,2003,32:1049-1074.

[51]Chen P,Sanyal J,Dudukovi M P.Numerical simulation of bubble columns flows:effect of different breakup and coalescence closures[J].Chemical Engineering Science,2005,60:1085-1101.

[52]Chevaillier J P,Klaseboer E,Masbernat O,Gourdon C.Effect of mass transfer on the film drainage between colliding drops[J].Colloid and Interface Science,2006,299:472-485.

[53]Klaseboer E,Chevaillier J P,Gourdon C,Masbernat O.Film drainage between colliding drops at constant approach velocity:Experiments and modeling[J].Colloid and Interface Science,2000,229:274-285.

[54]張阿漫,曾令玉,王詩平,楊樹濤,姚熊亮.水下爆炸氣泡融合動態(tài)特性研究[J].應(yīng)用數(shù)學(xué)和力學(xué),2010:31(2):163-170.

[55]Georgescu S C,Achard J L,Canot E.Jet drops ejection in bursting gas bubble processes[J].European Journal of Mechanics B/Fluids,2002,21:265-280.

[56]Malysa K,Krasowska,Krzan M.Influence of surface active substances on bubble motion and collision with various interfaces Advances in Colloid and Interface Science,2005,114-115:205-225.

[57]Kientzler C F,Arons A B,Blanchard D C,Woodcock A H.Photographic investigation of the projection of droplets by bubbles bursting at a water surface[J].Tellus,1954,6(1):1-7.

[58]Garner F H,Ellis S R M,Lacey J A.The size distribution and entrainment of droplets[J].Transaction of Institute Chemistry Engineering,1954,32:222-235.

[59]Resch F J,Afeti G M.Gas-liquid conversion at a free interface and its environmental implications[J].Fluid Mechanics and Thermodynamics,1991:1066-1076.

[60]Afeti G M,Resch F J.From breaking waves to breaking bubbles:Production characteristics of the marine droplet aerosol[J].Trends Geophysics Research,1992,1:113-125.

[61]Boulton-Stone J M,Blake J R.Gas bubbles bursting at a free surface[J].Journal of Fluid Mechanics,1993,254:437-466.

[62]Boulton-Stone J M.The effect of surfactant on bursting gas bubbles[J].Journal of Fluid Mechanics,1995,302:231-257.

[63]Sussman M,Smereka P.Axisymmetric free boundary problems[J].Journal of Fluid Mechanics,1997,341:269-294.

[64]Prosperetti A,Lezzi A.Bubble dynamics in a compressible liquid.Part I.First-order theory[J].Journal of Fluid Mechanics,1986,168:457-478.

[65]Herring C.Theory of the pulsations of the gas bubble produced by an underwater explosion[R].Office of Science Research and Development Report 236(NDRC Report C-4-sr 10-010,Columbia University),1941.

[66]Keller J B,Kolodner I I.Damping of underwater explosion bubble oscillations[J].Journal of Applied Physics,1956,27:1152-1161.

[67]Lezzi A,Prosperetti A.Bubble dynamics in a compressible liquid.Part II.Second-order Theory[J].Journal of Fluid Mechanics,1987,185:289-321.

[68]Pozrikidis C.Expansion of a compressible gas bubble in Stokes flow[J].Journal of Fluid Mechanics,2001,442:171-189.

[69]Darren G,Crowdy.Compressible bubbles in Stokes flow[J].Journal of Fluid Mechanics,2003,476:345-356.

[70]Wang Q X,Blake J R.Non-spherical bubble dynamics in a compressible liquid.Part 1:Travelling acoustic wave[J].Journal of Fluid Mechanics,2010,659:191-224.

[71]Wang Q X,Blake J R.Non-spherical bubble dynamics in a compressible liquid.Part 2:Acoustic standing wave[J].Journal of Fluid Mechanics,2011,679:559-581.

[72]Geers T L,Hunter K S.An integrated wave-effects model for an underwater explosion bubble[J].Acoustic Society of America,2002,111(4):1584-1601.

[73]王詩平.水中結(jié)構(gòu)物附近爆炸氣泡運動特性研究[J].哈爾濱:哈爾濱工程大學(xué),2011.

[74]Beeching R.Resistance to cavitation erosion[J].Transaction Institute of Engineering and Shipbuilders,1942,85:273-277.

[75]Akhatov I,Gumerov N,Ohl C D,Parlitz U,Lauterborn W.The role of surface tension in stable single-bubble sonoluminescence[J].Physical Review Letters,1997,78(2):227-230.

[76]Bronshtein I N,Semendyayev K A.Handbook of Mathematics[M].3rd edition.Springer Publication,Berlin,1997.

[77]Brackbill J U,Kothe D B,Zemach C.A continuum method for modeling surface tension[J].Journal of Computational Physics,1992,100:335-354.

[78]Zinchenko A Z,Rother A M,Davis R H,A novel boundary-integral algorithm for viscous interaction of deformable drops[J].Physics Fluids,1997,9(6):1493-1511.

[79]Rallison J M.A numerical study of the deformation and burst of a viscous drop in general shear flows[J].Journal of Fluid Mechanics,1981,109:465.

[80]黃繼湯,陳嘉范,丁 彤,田立言.表面張力對單空泡運動特性的影響[J].水利學(xué)報,1996,12:1-7.

[81]王起棣,張慧生.表面張力對固壁近旁空化氣泡潰滅特性的影響[J].復(fù)旦學(xué)報（自然科學(xué)版）,2003,42(2):201-207.

[82]張振宇,王起棣,張慧生.表面張力對近固壁二空化泡影響的數(shù)值研究[J].力學(xué)學(xué)報,2005,37(1):100-104.

[83]Miksis M J,Vanden-Broeck J M,Keller J B.Rising bubbles[J].Journal of Fluid Mechanics,1982,123:31-41.

[84]Lundgren S,Mansour N.Oscillations of drops in zero gravity with weak viscous effects[J].Journal of Fluid Mechanics,1988,194:479-510.

[85]Klaseboer E,Manica R,Chan D Y C,Khoo B C.BEM simulations of potential flow with viscous effects as applied to a rising bubble[J].Engineering Analysis with Boundary Elements,2011,35:489-494.

[86]Joseph D D,Wang J.The dissipation approximation and viscous potential flow[J].Journal of Fluid Mechanics,2004,505:365-377.

[87]Youngren G K,Acrivos A.On the shape of a gas bubble in a viscous extensional flow[J].Journal of Fluid Mechanics,1976,76(3):433-442.

[88]Barthes-Biesel D,Acrivosde A.Deformation and burst of a liquid droplet freely suspended in a linear shear field[J].Journal of Fluid Mechanics,1973,61:1.

[89]Buckmaster D.Pointed bubbles in slow viscous flow[J].Journal of Fluid Mechanics,1973,55:385.

[90]Eshpuniyani B,Fowlkes J B,Bull J L.A boundary element model of microbubble sticking and sliding in the microcirculation[J].International Journal of Heat and Mass Transfer,2008,51:5700-5711.

[91]Calderon A J,Eshpuniyani B,Fowlkes J B,Bull J L.A boundary element model of the transport of a semi-infinite bubble through a microvessel bifurcation[J].Physics of Fluids,2010,22,061092(1-11).

[92]陳 斌,Kawamura T,Kodama Y.靜止水中單個上升氣泡的直接數(shù)值模擬[J].工程熱物理學(xué)報,2005,26(6):980-982.

[93]李彥鵬,關(guān)衛(wèi)省.應(yīng)用Level Set方法模擬三維氣泡與液滴的運動[J].地球科學(xué)與環(huán)境學(xué)報,2007,29(2):217-220.

Review on the numerical simulation methods of the instability of bubble interface under complex flow field

NI Bao-yu,ZHANG A-man,ZHANG Zhong-yu,LI Shuai

(School of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China)

The study on the instability of the bubble interface is a hot topic.Among the various numerical methods,Boundary element method(BEM)is widely adopted to simulate the bubble dynamics due to its high efficiency.The development and application of BEM and relative numerical treatments on bubble dynamics are reviewed in this paper.The collapse of bubble near various boundaries under various ambient flows is summarized;The splitting,mergence and breakup of the bubble near free surface are concluded and discussed;the influences of three significant factors,namely compressibility,surface tension and viscosity,are reviewed and analyzed;and,the problems needed to be studied further are put forward.

bubble;interface instability;BEM;split;mergence

O35

A

10.3969/j.issn.1007-7294.2014.04.015

1007-7294（2014）04-0470-14

2013-08-03

國家自然科學(xué)基金重點項目（50939002）；國家自然科學(xué)基金委員會—中國工程物理研究院聯(lián)合基金資助項目（10976008）；國防基礎(chǔ)科研計劃（B2420110011）；國家自然科學(xué)基金青年項目（51009035）；黑龍江省自然科學(xué)基金（E201047,A200901）資助課題；中央高?；究蒲袠I(yè)務(wù)費專向資金資助（HEUCF140116）

倪寶玉（1986-），男，博士，哈爾濱工程大學(xué)船舶工程學(xué)院講師，E-mail:baoyuni@gmail.com；

張阿漫（1981-），男，教授，博士生導(dǎo)師。