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

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

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

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

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

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

氣泡界面不穩(wěn)定現(xiàn)象一直是國內(nèi)外關(guān)注和研究的熱點(diǎn)。近年來，邊界元方法（Boundary Element Method，BEM）被廣泛應(yīng)用于氣泡動(dò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é)者研究氣泡動(dòng)力學(xué)提供一定參考。

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

1 引 言

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

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

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

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

2 氣泡潰滅坍塌

基于勢(shì)流理論的邊界元法（BEM）或稱邊界積分法（BIM）在瞬態(tài)自由交界面上最早應(yīng)用，可能是Longuet-Higgins和Cokelet（1976）[9]。他們應(yīng)用第二格林函數(shù)將拉普拉斯方程轉(zhuǎn)化為邊界積分方程，并采用橢圓積分?jǐn)?shù)值求解積分方程，從而成功地解決了非線性瞬態(tài)水波問題。此后，邊界積分法迅速發(fā)展，在氣泡動(dòng)力學(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]）。邊界元法求解氣泡動(dòng)力學(xué)問題的基本思路如下：已知?dú)馀荼砻娴乃俣葎?shì)，應(yīng)用格林第二函數(shù)求解可獲得氣泡表面的法向速度?φ/?n，同時(shí)應(yīng)用速度勢(shì)差分法可獲得表面切向速度?φ/?τ，進(jìn)而可合成節(jié)點(diǎn)的全速度▽?duì)?。在時(shí)間步進(jìn)中，采用拉格朗日法追蹤氣泡表面節(jié)點(diǎn)，并應(yīng)用拉格朗日坐標(biāo)求解氣泡表面伯努利方程，獲得每一時(shí)間步內(nèi)速度勢(shì)增量Δφ。在獲得新的速度勢(shì)和新的氣泡節(jié)點(diǎn)坐標(biāo)后，即可進(jìn)入循環(huán)，獲得氣泡隨時(shí)間的推進(jìn)演化歷程。

在不同邊界或不同背景流場條件下，氣泡在收縮階段可能會(huì)發(fā)生不穩(wěn)定坍塌潰滅，并形成高速射流。最典型的兩類邊界是剛性壁面和自由液面。對(duì)于剛性壁面附近氣泡的求解中，對(duì)于氣泡表面，速度勢(shì)已知，滿足狄利克雷（Dirichelet）邊界條件；對(duì)于剛性壁面，法向速度已知，滿足紐曼（Neumann）條件，僅需對(duì)邊界積分方程得到的系數(shù)矩陣重組即可獲得所需未知量。研究發(fā)現(xiàn)氣泡在壁面附近膨脹階段會(huì)輕微遠(yuǎn)離壁面，而收縮階段會(huì)被壁面強(qiáng)烈吸引，形成指向壁面的射流。壁面對(duì)氣泡的這種效應(yīng)稱為“Bjerknes效應(yīng)”（Bjerknes，1966[19]）。對(duì)于自由液面附近的氣泡，自由面和氣泡面均滿足狄利克雷邊界條件，在時(shí)間步進(jìn)中也采用拉格朗日法，這樣即可采用與求解氣泡類似的方法追蹤自由面的變形。自由面附近的氣泡動(dòng)態(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)展，請(qǐng)參見文獻(xiàn)[22]和[23]。

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

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

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

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

圖2 自由拉格朗日法（上）[40]和邊界元法（下）[39]計(jì)算沖擊波下氣泡坍塌對(duì)比圖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)一步分析氣泡撕裂和融合的動(dòng)力學(xué)行為，尤其關(guān)注邊界元法在處理撕裂和融合問題的處理手段。

3.1 撕裂

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

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

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

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

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

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

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

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

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

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

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

3.2 融合

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

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

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

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

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

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

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

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

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

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

3.3 自由面處破裂

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

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

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

4 幾類主要影響因素

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

4.1 可壓縮性

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

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

4.2 表面張力

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

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

4.3 粘性

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

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

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

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

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

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

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

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

（4）目前業(yè)內(nèi)基本已公認(rèn)氣泡射流是造成附近結(jié)構(gòu)毀傷的主要載荷之一，然而氣泡射流對(duì)結(jié)構(gòu)的毀傷模式，以及不同爆距等參數(shù)如何影響氣泡射流對(duì)結(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]張阿漫,姚熊亮.基于邊界積分法的氣泡動(dòng)態(tài)特性綜述[J].力學(xué)進(jìn)展,2008,38(5):561-570.

[23]張阿漫,汪 玉,聞學(xué)友,倪寶玉,等.水下爆炸氣泡動(dòng)態(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)孚.空化泡和超空化泡流動(dòng)理論及應(yīng)用[M].北京:國防工業(yè)出版社,2009.

[34]韓寶玉,熊 鷹,葉金銘.梢渦空泡初生預(yù)報(bào)方法[J].海軍工程大學(xué)學(xué)報(bào),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]倪寶玉,張阿漫,朱 楓,吳巧瑞.艦船近場爆炸反射波對(duì)氣泡射流特性影響[J].哈爾濱工程大學(xué)學(xué)報(bào),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]張阿漫,曾令玉,王詩平,楊樹濤,姚熊亮.水下爆炸氣泡融合動(dòng)態(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)物附近爆炸氣泡運(yùn)動(dòng)特性研究[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]黃繼湯,陳嘉范,丁 彤,田立言.表面張力對(duì)單空泡運(yùn)動(dòng)特性的影響[J].水利學(xué)報(bào),1996,12:1-7.

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

[82]張振宇,王起棣,張慧生.表面張力對(duì)近固壁二空化泡影響的數(shù)值研究[J].力學(xué)學(xué)報(bào),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.靜止水中單個(gè)上升氣泡的直接數(shù)值模擬[J].工程熱物理學(xué)報(bào),2005,26(6):980-982.

[93]李彥鵬,關(guān)衛(wèi)省.應(yīng)用Level Set方法模擬三維氣泡與液滴的運(yùn)動(dòng)[J].地球科學(xué)與環(huán)境學(xué)報(bào),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é)基金重點(diǎn)項(xiàng)目（50939002）；國家自然科學(xué)基金委員會(huì)—中國工程物理研究院聯(lián)合基金資助項(xiàng)目（10976008）；國防基礎(chǔ)科研計(jì)劃（B2420110011）；國家自然科學(xué)基金青年項(xiàng)目（51009035）；黑龍江省自然科學(xué)基金（E201047,A200901）資助課題；中央高校基本科研業(yè)務(wù)費(fèi)專向資金資助（HEUCF140116）

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

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