流體黏性及表面張力對(duì)氣泡運(yùn)動(dòng)特性的影響?

艾旭鵬 倪寶玉

(哈爾濱工程大學(xué)船舶工程學(xué)院,哈爾濱 150001)

流體黏性及表面張力對(duì)氣泡運(yùn)動(dòng)特性的影響?

艾旭鵬 倪寶玉?

(哈爾濱工程大學(xué)船舶工程學(xué)院,哈爾濱 150001)

氣泡,邊界層,黏性,表面張力

1 引 言

黏性流體中多氣泡的運(yùn)動(dòng)和耦合作用在生物、化學(xué)、海洋工程中應(yīng)用十分廣泛,比如組織和細(xì)胞損傷[1]、水下爆炸氣泡[2]、微氣泡減阻[3]等.從數(shù)學(xué)模型和數(shù)值模擬角度而言,這一問(wèn)題主要涉及兩大有趣但又具有挑戰(zhàn)性的難題,一是氣泡與氣泡間的強(qiáng)非線(xiàn)性耦合效應(yīng),二是流體黏性效應(yīng)對(duì)于氣泡耦合效應(yīng)的影響.對(duì)于第一個(gè)問(wèn)題,目前較常用的有效解決方法之一是采用完全非線(xiàn)性邊界元(BEM)方法.文獻(xiàn)[4—8]針對(duì)這一問(wèn)題進(jìn)行了大量的研究,取得了較大的進(jìn)展.對(duì)于第二個(gè)問(wèn)題,一方面可以采用計(jì)算流體力學(xué)(CFD)方法近似求解Navier-Stokes(N-S)方程[9,10],但是計(jì)算量和計(jì)算精度都有待進(jìn)一步優(yōu)化;另一方面,采用“邊界層”理論,將黏性效應(yīng)限制于很薄的流體層內(nèi),而層外則依然可以采用勢(shì)流理論進(jìn)行計(jì)算.邊界層理論的思想最早源于普朗特的平板阻力,在船舶領(lǐng)域應(yīng)用廣泛.但是有關(guān)氣泡邊界層理論,發(fā)展和研究都相對(duì)較晚.

Miksis等最早在氣泡動(dòng)力學(xué)中引入邊界層理論的思想[11].Miksis等認(rèn)為氣泡周?chē)嬖谝槐∵吔鐚?邊界層外流體滿(mǎn)足無(wú)黏不可壓縮勢(shì)流理論,邊界層內(nèi)考慮流體黏性效應(yīng),邊界層兩側(cè)滿(mǎn)足法向應(yīng)力連續(xù),但不滿(mǎn)足切向應(yīng)力連續(xù).由于此邊界層很薄,忽略邊界層內(nèi)外速度的變化,這樣便可采用邊界積分法數(shù)值模擬計(jì)及弱黏性效應(yīng)的氣泡動(dòng)態(tài)變化.此后,研究人員從不同角度改進(jìn)Miksis等的氣泡邊界層理論.Lundgren和Mansour[12]進(jìn)一步分析了邊界層方程,獲得了邊界層兩側(cè)壓力差和法向速度差的表達(dá)式.Boulton-Stone和Blake[1,13]又進(jìn)一步擴(kuò)展了Lundgren和Mansour的方法,使得邊界層兩側(cè)速度的切向分量也連續(xù),在此基礎(chǔ)上首次成功地?cái)?shù)值模擬了單個(gè)軸對(duì)稱(chēng)氣泡在自由面處的破裂現(xiàn)象.Georgescu等[14]也采用Miksis的氣泡邊界層理論模擬氣泡在自由面處的破裂現(xiàn)象,在邊界積分和時(shí)間步進(jìn)上均采用二階單元,以獲得更好的精度.Klaseboer等[15]采用Joseph和Wang[16]的思想修正邊界層理論,同時(shí)考慮了法向和切向應(yīng)力的連續(xù),應(yīng)用邊界積分法數(shù)值模擬了高雷諾數(shù)下軸對(duì)稱(chēng)氣泡在流體中的脈動(dòng)上浮現(xiàn)象.Lind和Phillips[17,18]分別研究了考慮法向應(yīng)力連續(xù)的軸對(duì)稱(chēng)氣泡在固壁面和自由液面附近的脈動(dòng).倪寶玉等[19]聯(lián)合模型實(shí)驗(yàn)和數(shù)值模擬研究了常壓小氣泡在自由表面的破裂,數(shù)值模擬中考慮了法向應(yīng)力連續(xù)的軸對(duì)稱(chēng)模型,主要分析了射流斷裂后的水滴形態(tài)以及水滴撕裂現(xiàn)象.Zhang和Ni[20]研究了軸對(duì)稱(chēng)和三維氣泡在法向和切向應(yīng)力同時(shí)連續(xù)的情況下在無(wú)界重力場(chǎng)中的脈動(dòng),得到氣泡射流速度隨流體黏性增加而衰減的結(jié)論.Li和Ni[21]在此基礎(chǔ)上又進(jìn)一步研究了氣泡與自由面相互作用過(guò)程中流體黏性效應(yīng)的影響.

以上所述的數(shù)值研究中對(duì)于黏性效應(yīng)的處理大部分側(cè)重于法向應(yīng)力連續(xù),對(duì)于切向應(yīng)力連續(xù)的考慮較少;數(shù)值模型中以軸對(duì)稱(chēng)為主,三維模型則十分少見(jiàn);而對(duì)于考慮黏性效應(yīng)和表面張力下多氣泡耦合作用的三維模型公開(kāi)發(fā)表的文獻(xiàn)還很少.鑒于此,本文在前人基礎(chǔ)之上,首先將建立的數(shù)值模型與Rayleigh-Plesset的解析解進(jìn)行對(duì)比,驗(yàn)證數(shù)值模型的有效性;其次研究不同雷諾數(shù)下雙氣泡耦合過(guò)程中出現(xiàn)的各類(lèi)物理現(xiàn)象及變化規(guī)律;最后研究雷諾數(shù)和韋伯?dāng)?shù)對(duì)于單個(gè)球狀脈動(dòng)氣泡的影響規(guī)律.

2 理論與方法

圖1給出了雙氣泡耦合作用的示意圖,這里假定初始位置較低的氣泡為“氣泡1”,位置較高的為“氣泡2”.定義笛卡爾坐標(biāo)O-xyz,坐標(biāo)系的原點(diǎn)處于氣泡1的初始中心處,z軸豎直向上.氣泡1與氣泡2的初始中心距離為h.

圖1 雙氣泡耦合作用示意圖Fig.1.Sketch of interaction between two bubbles.

2.1 控制方程

本文主要考慮不可壓縮、動(dòng)力黏性系數(shù)恒定的牛頓流體,根據(jù)亥姆霍茲速度分解定理[22],流場(chǎng)中任意一點(diǎn)速度U均可以分解為無(wú)旋的速度u=?φ和有旋的速度v=?×a.對(duì)于不可壓縮流體的質(zhì)量守恒方程則為

式中?是哈密爾頓算子,φ是速度勢(shì),a是一個(gè)非零矢量勢(shì).對(duì)于不可壓縮流體的動(dòng)量守恒方程為N-S方程:

式中t是時(shí)間,g是重力加速度,P是壓力,ρ是流體密度,ν是流體的運(yùn)動(dòng)黏性系數(shù).將u=?φ代入到方程(2)并考慮到無(wú)旋運(yùn)動(dòng)的特點(diǎn),方程(2)可簡(jiǎn)化為

式中P∞是無(wú)窮遠(yuǎn)處參考?jí)毫?模型中不考慮背景聲壓,氣泡的背景壓力只有靜水壓力.Pvc為黏性壓力,滿(mǎn)足

用(4)式直接求解Pvc十分困難,本文采用間接求解方法.

2.2 邊界層理論

對(duì)于法向應(yīng)力,在第i個(gè)氣泡的界面處,根據(jù)楊氏-拉普拉斯方程,可知法向應(yīng)力平衡條件為

式中下標(biāo)i表示第i個(gè)氣泡,i可取1或者2;Pl是界面處液體壓力,Pb是界面處氣體壓力,對(duì)于氣泡有Pb=Pc+P0(V0/V)ι,其中Pc是飽和蒸汽壓,P0和V0分別是氣泡初始形成時(shí)的壓力和體積,V是氣泡體積,ι是氣泡內(nèi)氣體的比熱比;τn是法向黏性應(yīng)力;σ是表面張力系數(shù);κ是界面的局部曲率.當(dāng)流場(chǎng)的無(wú)旋速度?φ在總速度中占主導(dǎo)地位時(shí),在氣泡附近很薄的邊界層內(nèi)法向黏性應(yīng)力可近似表達(dá)為滿(mǎn)足

對(duì)于切向應(yīng)力,由于氣泡內(nèi)氣體的動(dòng)力黏性系數(shù)相對(duì)于流體十分小,故可認(rèn)為在氣泡表面切向黏性應(yīng)力為0.但是,無(wú)旋速度?φ在界面處誘導(dǎo)的切向應(yīng)力卻不為零,即τs=2ρυ(?φτ/?n)/=0. 如果直接采用切向應(yīng)力為零的切向應(yīng)力連續(xù)條件,則無(wú)旋速度誘導(dǎo)的切向應(yīng)力在邊界處的做功則被人為忽略了.為了彌補(bǔ)這一非零的切向應(yīng)力的貢獻(xiàn),根據(jù)Joseph和Wang[16]的推導(dǎo),假設(shè)在氣泡附近很薄的邊界層內(nèi),切向應(yīng)力τs的切向做功與黏性壓力的法向做功相等:

式中t和n分別是邊界的切向單位矢量和法向單位矢量,S是流域的所有邊界,包含氣泡1和氣泡2的邊界Sb1和Sb2.

如前所述,采用(4)式直接求解黏性壓力Pvc在數(shù)值上具有極大困難,這使得我們不得不尋找求解Pvc的其他方法.這里引用Joseph和Wang[16]的一個(gè)經(jīng)驗(yàn)性假設(shè)條件,認(rèn)為黏性壓力與無(wú)旋速度誘導(dǎo)的法向應(yīng)力成正比,即

其中C是未知的比例系數(shù),考慮到(6)式中的τs在三維問(wèn)題中難以直接求解,這里采用能量守恒法進(jìn)行求解.對(duì)于三維工況,首先需要將?2φ/?n2轉(zhuǎn)化到笛卡爾坐標(biāo)系中,有

根據(jù)Lamb[22]的公式,對(duì)于無(wú)旋流場(chǎng),黏性能量耗散率D可表達(dá)為

式中(u,v,w)是速度矢量u在笛卡爾坐標(biāo)系中的三個(gè)分量.另一方面,通過(guò)功能關(guān)系,黏性能量耗散率D與單位時(shí)間黏性應(yīng)力做功W相等:

將方程(6),(7),(8)代入方程(10),再與方程(9)聯(lián)立,最終可獲得未知系數(shù)C的表達(dá)式為

2.3 邊界條件

將方程(5)和(7)代入方程(3)中,并考慮到物質(zhì)導(dǎo)數(shù)與局部導(dǎo)數(shù)間的關(guān)系?,則可以獲得在第i個(gè)氣泡表面的全非線(xiàn)性動(dòng)力學(xué)邊界條件:

拉格朗日系統(tǒng)中第i個(gè)氣泡的全非線(xiàn)性運(yùn)動(dòng)學(xué)邊界條件為

式中x=(x,y,z)是氣泡上任一個(gè)質(zhì)點(diǎn)的矢量坐標(biāo).最后,在無(wú)窮遠(yuǎn)處擾動(dòng)為零,則無(wú)窮遠(yuǎn)邊界條件為

2.4 三維邊界元方法和無(wú)量綱化

采用格林第三公式和格林函數(shù)G,方程(1)可轉(zhuǎn)化為如下的邊界積分方程:

式中S是包含氣泡表面Sb1和Sb2以及無(wú)窮遠(yuǎn)邊界S∞在內(nèi)的所有邊界;p和q分別是場(chǎng)點(diǎn)和源點(diǎn),ε(p)是在p點(diǎn)處觀(guān)察流場(chǎng)的立體角;三維格林函數(shù)G(p,q)=1/|R?r|,其中R和r分別是p和q的位置矢量.

對(duì)于y和z方向有類(lèi)似的公式.對(duì)于偏微分方程(16)可以參照方程(15)的方法進(jìn)行求解,將求解的結(jié)果代入方程(8)和(11),即可得到和C.

其中?_i是第i個(gè)氣泡的強(qiáng)度參數(shù),不同氣泡的強(qiáng)度參數(shù)可以一樣也不可以不一樣;Rmax為所有氣泡在無(wú)界自由場(chǎng)中球狀運(yùn)動(dòng)時(shí)能達(dá)到的最大半徑,初始強(qiáng)度參數(shù)越大的氣泡Rmax也越大.后文中無(wú)量綱量均用“一拔”表示.

3 結(jié)果與討論

3.1 有效性驗(yàn)證

為了校核本文數(shù)值模型的有效性,這里采用單個(gè)球狀氣泡為對(duì)象,忽略重力的影響(Fr→∞),分別采用解析解和數(shù)值模型求解單個(gè)球狀氣泡在黏性流體中的脈動(dòng).當(dāng)氣泡是球狀時(shí),氣泡半徑的脈動(dòng)方程可由Rayleigh-Plesset方程[24]求解:

另一方面,可采用本文建立的邊界元模型進(jìn)行求解,且考慮到球狀氣泡為軸對(duì)稱(chēng)模型,可分別采用軸對(duì)稱(chēng)邊界元和三維邊界元獲得數(shù)值結(jié)果.將數(shù)值結(jié)果與解析解進(jìn)行對(duì)比分析,驗(yàn)證數(shù)值程序.

選取初始無(wú)量綱參數(shù)為?=100,We=通過(guò)初始條件(19)式可以計(jì)算得到無(wú)量綱初始半徑ˉR0=0.1485.對(duì)應(yīng)的軸對(duì)稱(chēng)模型母線(xiàn)離散為33個(gè)節(jié)點(diǎn),32個(gè)線(xiàn)性單元;三維模型表面離散為362個(gè)節(jié)點(diǎn),720個(gè)三角形單元.圖2給出了球狀氣泡半徑隨時(shí)間變化的軸對(duì)稱(chēng)、三維數(shù)值結(jié)果與精確解的對(duì)比曲線(xiàn).圖中圖標(biāo)Re(∞)表示雷諾數(shù)趨于無(wú)窮大,即不考慮黏性效應(yīng)的工況;Re(102)表示采用Rayleigh-Plesset方程(18)直接求解的精確解;axi表示對(duì)應(yīng)的軸對(duì)稱(chēng)邊界元模型的結(jié)果;3d表示對(duì)應(yīng)的三維邊界元模型的結(jié)果.由圖2可見(jiàn),不考慮黏性效應(yīng)的情況下,氣泡在每次膨脹到最大時(shí)刻均可精確達(dá)到無(wú)量綱最大半徑1,且無(wú)量綱運(yùn)動(dòng)周期維持在1.99不變.此時(shí)可將氣泡視為一個(gè)無(wú)能量損耗的彈簧振子,在內(nèi)外壓力差和慣性作用下不斷往復(fù)地在平衡位置振動(dòng).然而,在考慮了黏性效應(yīng)之后,氣泡每次膨脹到最大的半徑將無(wú)法達(dá)到1,且隨著脈動(dòng)而逐漸減小,前三周期最大半徑依次為0.966,0.914,0.873;相應(yīng)地,收縮到最小的半徑也大于初始半徑ˉR0=0.1485,且隨脈動(dòng)而逐漸增大.同時(shí),氣泡脈動(dòng)的周期逐漸減小,首次和二次依次為1.95和1.89.

圖2 (網(wǎng)刊彩色)黏性球狀氣泡半徑時(shí)歷曲線(xiàn)Fig.2.(color online)Time histories of bubble radius with viscous effect of fluid.

對(duì)于軸對(duì)稱(chēng)邊界元模型,在前三周期氣泡半徑最大時(shí)刻與精確解的相對(duì)誤差分別為0.76%,1.01%和1.13%,對(duì)于三維邊界元模型,在前兩周期氣泡最大半徑時(shí)的相對(duì)誤差分別為0.19%和0.36%.可見(jiàn)軸對(duì)稱(chēng)和三維模型的數(shù)值解與Rayleigh-Plesset方程精確解符合度很高,證明了本文建立的邊界元模型是有效可行的.由于軸對(duì)稱(chēng)模型與三維模型的網(wǎng)格構(gòu)造方法不同:對(duì)于軸對(duì)稱(chēng)模型,只需在母線(xiàn)上劃分若干段線(xiàn)段網(wǎng)格;對(duì)于三維模型,需要在整個(gè)球面上劃分若干個(gè)三角形網(wǎng)格.數(shù)值計(jì)算中網(wǎng)格劃分方法會(huì)對(duì)計(jì)算結(jié)果有一定影響,所以在數(shù)值積分時(shí),二者結(jié)果會(huì)有一定差異.從圖中可見(jiàn)兩種模型的差異十分微小,在數(shù)值計(jì)算中是可以理解和接受的.此時(shí)氣泡可視為計(jì)入結(jié)構(gòu)內(nèi)黏性的彈簧振子模型,在往復(fù)振動(dòng)中不斷有黏性能量消耗.Zhang等[25]給出了考慮流場(chǎng)可壓縮性后氣泡半徑脈動(dòng)曲線(xiàn),也觀(guān)察到了類(lèi)似的半徑衰減的現(xiàn)象,然而二者產(chǎn)生的原因是截然不同的,計(jì)及黏性效應(yīng)時(shí)初始?xì)馀輧?nèi)能不斷被黏性耗散,而可壓縮性是將初始內(nèi)能不斷傳遞并儲(chǔ)存到流場(chǎng)中.

圖5給出了球狀氣泡無(wú)量綱黏性耗散率的變化曲線(xiàn).

圖3 (網(wǎng)刊彩色)黏性球狀氣泡徑向速度時(shí)歷曲線(xiàn)Fig.3.(color online)Time histories of radial velocity of the bubble with viscous effect of fluid.

圖4 (網(wǎng)刊彩色)黏性球狀氣泡內(nèi)壓時(shí)歷曲線(xiàn)Fig.4.(color online)Time histories of inner pressure of the bubble with viscous effect of fluid.

圖5 (網(wǎng)刊彩色)球狀氣泡黏性耗散率時(shí)歷曲線(xiàn)Fig.5.(color online)Time histories of viscous dissipation rate of the bubble

3.2 典型結(jié)果分析

在驗(yàn)證數(shù)值模型有效性的基礎(chǔ)上,研究?jī)蓚€(gè)串列的氣泡的相互作用.假定兩個(gè)氣泡初始均靜止且二者內(nèi)壓相同,即強(qiáng)度參數(shù)一致.這里給定初始參數(shù)分別為?_1=?_2=100,We=1.37×106,Re=100,Fr=5和λ=2.5.零時(shí)刻氣泡表面的初值條件為和

圖6 (網(wǎng)刊彩色)重力場(chǎng)中考慮流體黏性效應(yīng)的兩氣泡演化過(guò)程(a)ˉt≈0.0241;(b)ˉt≈1.2336;(c)ˉt≈1.9911;(d) ˉt≈ 2.2025Fig.6.(color online)The evolution of the bubbles with viscous effect of fluid under gravity:(a)ˉt≈0.0241;(b)ˉt≈1.2336;(c)ˉt≈1.9911;(d)ˉt≈2.2025.

圖6給出了考慮流體黏性效應(yīng)的兩個(gè)氣泡耦合作用的演化過(guò)程,其中云圖表示速度勢(shì)大小(紅色表示高、藍(lán)色表示低).圖6(a)對(duì)應(yīng)氣泡初始狀態(tài),兩個(gè)球狀氣泡在各自高壓的作用下向外膨脹,由于氣泡間的耦合作用,相互臨近的表面的速度勢(shì)和速度較低,而相互遠(yuǎn)離的表面的速度勢(shì)和速度較高.圖6(b)對(duì)應(yīng)上氣泡體積達(dá)到最大的時(shí)刻,此時(shí)下氣泡的下表面在浮力的作用下已經(jīng)開(kāi)始進(jìn)入收縮階段.從圖6(c)中可清晰地觀(guān)察到下氣泡的下表面已經(jīng)變得平坦,而上氣泡在下氣泡的吸引作用下呈現(xiàn)“鴨蛋狀”.圖6(d)對(duì)應(yīng)下氣泡坍塌結(jié)束時(shí)刻,此時(shí)可見(jiàn)在浮力、流體黏性、氣泡之間的吸引力(又稱(chēng)為Bjerknes力)等多種耦合作用的影響下,下氣泡的下表面內(nèi)卷形成典型的氣泡射流,射流即將穿透下氣泡而形成環(huán)狀射流.同時(shí)上氣泡在多種力的作用下,仍然以“鴨蛋狀”較慢地收縮.

由于缺乏可對(duì)比的實(shí)驗(yàn)數(shù)據(jù),這里分別采用軸對(duì)稱(chēng)邊界元模型和三維邊界元模型計(jì)算圖6的工況,并將計(jì)算得到的射流速度曲線(xiàn)進(jìn)行對(duì)比,如圖7所示,其中viscous axi和viscous 3d分別是考慮黏性作用的軸對(duì)稱(chēng)模型和三維模型.本文將射流速度定義為氣泡下表面內(nèi)凹最高點(diǎn)處流體質(zhì)點(diǎn)的法向速度.從圖7的對(duì)比可見(jiàn)二者計(jì)算的結(jié)果精確符合,再度驗(yàn)證了數(shù)值模型的有效性.此外,為了突出流體黏性效應(yīng)的作用,再采用三維邊界元模型計(jì)算不考慮黏性作用的圖6工況,即令Re無(wú)窮大,而保持其他參數(shù)不變.不考慮黏性效應(yīng)的計(jì)算結(jié)果與考慮黏性效應(yīng)的計(jì)算結(jié)果對(duì)比顯示,考慮黏性作用后,射流尖端的速度降低,射流砰擊氣泡表面的時(shí)間推遲.表明流體的黏性效應(yīng)將抑制氣泡射流的發(fā)展,降低氣泡射流的速度和能量.

圖7 (網(wǎng)刊彩色)兩個(gè)氣泡耦合作用中下氣泡射流速度曲線(xiàn)對(duì)比圖Fig.7.(color online)Variation of jet velocity of the lower bubble during interaction of two bubbles.

在進(jìn)一步討論黏性效應(yīng)對(duì)能量的影響之前,這里給出流場(chǎng)動(dòng)能、勢(shì)能和黏性耗散能及總能量的定義.首先,整個(gè)流場(chǎng)的無(wú)量綱動(dòng)能如下:

式中v為流場(chǎng)的體積,S為包含氣泡表面在內(nèi)的流域所有邊界面.氣泡的無(wú)量綱勢(shì)能[26]定義如下:

式中無(wú)量綱黏性耗散率ˉD可通過(guò)(9)或(10)式數(shù)值求得.

氣泡在脈動(dòng)過(guò)程中,守恒的總能量Etotal除了動(dòng)能、勢(shì)能相互轉(zhuǎn)化的機(jī)械能Em外,還將存在一部分黏性耗散能Ed,即

(23)式的最后一個(gè)等號(hào)是因?yàn)槌跏紩r(shí)刻流場(chǎng)動(dòng)能和黏性耗散能均為0.

圖8給出了圖6對(duì)應(yīng)的工況中兩個(gè)氣泡相互作用過(guò)程中,流場(chǎng)中能量相互轉(zhuǎn)換的曲線(xiàn).從圖中可見(jiàn),隨著兩個(gè)氣泡的運(yùn)動(dòng),流場(chǎng)能量由最開(kāi)始的氣泡勢(shì)能轉(zhuǎn)換為兩個(gè)氣泡的動(dòng)能以及流場(chǎng)的黏性耗散能.對(duì)于動(dòng)能,明顯地隨著氣泡的膨脹、坍塌等階段呈現(xiàn)波峰波谷的特性;黏性耗散能則不同,它隨著氣泡的運(yùn)動(dòng)而不斷上升,尤其在氣泡射流坍塌即將發(fā)生階段,上升得更快,說(shuō)明黏性摩擦的作用在氣泡射流速度高的時(shí)候表現(xiàn)得更明顯.總體而言,總能量在整個(gè)變化過(guò)程中數(shù)值誤差允許范圍內(nèi)是守恒的.可見(jiàn)黏性效應(yīng)在氣泡射流速度、流場(chǎng)能量的轉(zhuǎn)化中是很重要的.

圖8 (網(wǎng)刊彩色)兩個(gè)氣泡耦合作用中流場(chǎng)無(wú)量綱能量的轉(zhuǎn)化過(guò)程Fig.8.(color online)Variation of dimensionless energy during interaction of two bubbles.

3.3 影響因素分析

本節(jié)分析雷諾數(shù)和韋伯?dāng)?shù)對(duì)于氣泡運(yùn)動(dòng)過(guò)程的影響規(guī)律.為了清晰地觀(guān)察二者的影響,忽略重力影響(Fr→∞),采用簡(jiǎn)單的球狀氣泡為例,計(jì)算不同雷諾數(shù)和韋伯?dāng)?shù)下氣泡脈動(dòng)特性,同時(shí)在考慮其中任一參數(shù)(稱(chēng)為主要參數(shù))時(shí),令另外一參數(shù)(稱(chēng)為次要參數(shù))足夠大,以降低次要參數(shù)對(duì)主要參數(shù)的影響,從而突出考察黏性效應(yīng)和表面張力效應(yīng).

圖9給出了不同雷諾數(shù)下氣泡半徑變化曲線(xiàn),其中Re(∞)仍表示對(duì)應(yīng)不考慮黏性效應(yīng)的工況.從圖中可見(jiàn),氣泡半徑在高雷諾數(shù)下變化很不敏感,雷諾數(shù)Re在104時(shí),氣泡半徑基本與不考慮黏性效應(yīng)是重合的;Re在104—103之間變化時(shí),氣泡半徑損失十分微弱,在前兩個(gè)周期內(nèi),可以忽略黏性效應(yīng);僅當(dāng)Re在102以下時(shí),氣泡半徑損失明顯,黏性效應(yīng)必須考慮.

圖9 (網(wǎng)刊彩色)不同雷諾數(shù)下球狀氣泡半徑的變化Fig.9.(color online)Time histories of spherical bubble radius at different Re number.

圖10給出了當(dāng)雷諾數(shù)Re分別取102—103間某幾個(gè)典型數(shù)值時(shí)黏性耗散率的變化曲線(xiàn).從圖中可見(jiàn)隨著Re減小,黏性耗散率逐漸增大,且隨著時(shí)間推進(jìn),峰值衰減更快.數(shù)值計(jì)算表明,當(dāng)Re大于103之后,黏性耗散率的數(shù)值已經(jīng)很小,故沒(méi)有在圖中繪制.

圖10 (網(wǎng)刊彩色)不同雷諾數(shù)下球狀氣泡黏性耗散率的變化Fig.10.(color online)Time histories of spherical bubble radius at different Reynolds number.

圖11給出了不同韋伯?dāng)?shù)We下氣泡半徑變化曲線(xiàn).從圖中可見(jiàn),無(wú)論We如何變化,氣泡的無(wú)量綱最大半徑一直是1,即韋伯?dāng)?shù)不會(huì)影響氣泡脈動(dòng)的幅值,但是會(huì)影響脈動(dòng)的周期.而且當(dāng)We在高于106量級(jí)時(shí),氣泡的半徑幾乎重合,可以忽略表面張力效應(yīng);當(dāng)We處于106—102之間時(shí),氣泡半徑和周期變化不敏感;當(dāng)We取101量級(jí)及以下時(shí),表面張力效應(yīng)顯著,氣泡脈動(dòng)周期顯著減小,即氣泡振動(dòng)更快.

圖11 (網(wǎng)刊彩色)不同韋伯?dāng)?shù)下球狀氣泡半徑的變化Fig.11.(color online)Time histories of spherical bubble radius at different Weber number.

圖12給出了韋伯?dāng)?shù)分別取137和13.7兩種工況下流場(chǎng)動(dòng)能和勢(shì)能以及總能量的變化曲線(xiàn).由于此時(shí)雷諾數(shù)Re=1010,黏性耗散能已經(jīng)十分微小,可忽略不計(jì).從圖中的對(duì)比曲線(xiàn)可見(jiàn),韋伯?dāng)?shù)的增加對(duì)動(dòng)能影響較小,對(duì)勢(shì)能影響很大,這是因?yàn)閯?shì)能中含有表面張力效應(yīng)項(xiàng)

圖12 (網(wǎng)刊彩色)不同韋伯?dāng)?shù)下球狀氣泡能量的變化Fig.12. (color online)Time histories of energy of spherical bubble at different Weber number.

總體而言,黏性效應(yīng)會(huì)使得氣泡振動(dòng)幅值衰減,同時(shí)氣泡周期縮短;而表面張力會(huì)使得氣泡周期縮短但幅值不變,同時(shí)勢(shì)能增加.然而,氣泡半徑對(duì)高雷諾數(shù)和高韋伯?dāng)?shù)變化均不敏感,此時(shí)基本可以忽略黏性效應(yīng)和表面張力效應(yīng).

4 結(jié) 論

基于氣泡邊界層理論,修正現(xiàn)有的勢(shì)流理論,建立考慮流體弱黏性效應(yīng)的三維邊界元模型,研究考慮流體黏性效應(yīng)和表面張力效應(yīng)下氣泡的脈動(dòng)和氣泡間的耦合作用.首先,采用Rayleigh-Plesset方程求解球狀氣泡,將邊界元模型的計(jì)算結(jié)果與解析解進(jìn)行對(duì)比,驗(yàn)證數(shù)值模型的有效性.在此基礎(chǔ)上,選擇兩個(gè)氣泡相互耦合作為典型算例,分析考慮黏性作用下氣泡耦合、氣泡射流和流場(chǎng)能量等的變化.最后研究雷諾數(shù)和韋伯?dāng)?shù)兩個(gè)參數(shù)對(duì)氣泡脈動(dòng)的影響規(guī)律.得到以下主要結(jié)論:

1)采用黏性壓力法向做功與切向應(yīng)力切向做功等效的方法,可巧妙地解決無(wú)旋速度?φ在界面處誘導(dǎo)的切向應(yīng)力不為零,但勢(shì)流理論中選取切向應(yīng)力為零作為邊界條件的矛盾,從而彌補(bǔ)了經(jīng)典勢(shì)流理論中忽略了切向應(yīng)力在邊界處做功的缺失;

2)就非球狀氣泡而言,當(dāng)考慮氣/液交界面處法向和切向黏性效應(yīng)后,氣泡射流在黏性流體中速度減慢,氣泡射流能量降低,流場(chǎng)損失的動(dòng)能和勢(shì)能逐漸轉(zhuǎn)化為黏性耗散能,且在射流速度較大的坍塌后期,黏性耗散能增長(zhǎng)得越快;

3)就球狀氣泡而言,黏性效應(yīng)會(huì)使得氣泡振動(dòng)幅值衰減,同時(shí)氣泡周期縮短,黏性耗散能不斷增加;表面張力會(huì)使得氣泡周期縮短但脈動(dòng)幅值不變,同時(shí)氣泡整體勢(shì)能增加;當(dāng)雷諾數(shù)和韋伯?dāng)?shù)較大后,可忽略黏性和表面張力的影響.

[1]Boulton-Stone J M,Blake J R 1993J.Fluid Mech.254 437

[2]Cui P,Zhang A M,Wang S P 2016Phys.Fluids28 117103

[3]Xue Y Z,Cui B,Ni B Y 2016Ocean Eng.118 58

[4]Sato K,Tomita Y,Shima A 1994J.Acoust.Soc.Am.95 2416

[5]Rungsiyaphornrat S,Klaseboer E,Khoo B C,Yeo K S 2003Comput.Fluids32 1049

[6]Chew L W,Klaseboer E,Ohl S W,Khoo B C 2011Phys.Rev.E84 0663078

[7]Han R,Li S,Zhang A M,Wang Q X 2016Phys.Fluids28 062104

[8]Guo X,Cai C,Xu G,Yang Y,Tu J,Huang P,Zhang D 2017Ultrason.Sonoche.39 863

[9]Zhang Y,Zhang Y,Li S 2017Ultrason.Sonoche.35 431

[10]Liu L T,Yao X L,Zhang A M,Chen Y Y 2017Phys.Fluids29 012105

[11]Miksis M J,Vanden-Broeck J M,Keller J B 1982J.Fluid Mech.123 31

[12]Lundgren S,Mansour N 1988J.Fluid Mech.194 479

[13]Boulton-Stone J M 1995J.Fluid Mech.302 231

[14]Georgescu S C,Achard J L,Canot E 2002Euro.J.Mech.B-Fluids21 265

[15]Klaseboer E,Manica R,Chan D Y C,Khoo B C 2011Eng.Anal.Bound.Elem.35 489

[16]Joseph D D,Wang J 2004J.Fluid Mech.505 365

[17]Lind S J,Phillips T N 2012Theor.Comp.Fluid Dyn.26 245

[18]Lind S J,Phillips T N 2013Phys.Fluids25 022014

[19]Ni B Y,Li S,Zhang A M 2013Acta Phys.Sin.62 124704(in Chinese)[倪寶玉,李帥,張阿漫 2013物理學(xué)報(bào) 62 124704]

[20]Zhang A M,Ni B Y 2014Comput.Fluids92 22

[21]Li S,NI B Y 2016Eng.Anal.Bound.Elem.68 63

[22]Lamb H 1932Hydrodynamics(sixth Ed.)(Cambridge:Cambridge University Press)pp580–581

[23]Wu G X 1991Appl.Ocean Res.13 317

[24]Rayleigh J W 1917Philos.Mag.34 94

[25]Zhang A M,Wang S P,Wu G X 2013Eng.Anal.Bound.Elem.37 1179

[26]Best J P 1993J.Fluid Mech.251 79

In fluence of viscosity and surface tension of fluid on the motion of bubbles?

Ai Xu-Peng Ni Bao-Yu?

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

18 July 2017;revised manuscript

12 August 2017)

Boundary integral simulation has been conducted to study the motion and deformation of bubbles with weak viscous and surface tension effects in fluid.Both normal and tangential stress boundary conditions are satisfied and the weak viscous effects are con fined to the thin boundary layers around bubble surfaces,which is also known as boundary layer theory of bubble.By using this method,the influence of viscosity and surface tension of fluid on the motion of bubbles has been studied.Both axisymmetric and three-dimensional numerical results are compared with analytical results of Rayleigh-Plesset equation.Good agreement between them is achieved,which validates the numerical model.On this basis,interaction model between two vertically placed bubbles is established,by taking the surface tension,gravity,and viscous effects into consideration.Variations of physical quantities including bubble deformation,jet velocity,and energy of fluid are studied.Last but not least,the influence of viscosity and surface tension on the motion of a spherical bubble is investigated.It is found that viscous effects of fluid depress the pulsation of bubble and part of fluid energy is transformed into viscous dissipation energy.As a result,the development of bubble jet,the radius of the bubble,and the jet velocity are reduced gradually.On the other hand,the surface tension of fluid does not change the range of the bubble pulsation but reduces the period of the bubble pulsation and enhances the potential energy of the bubble.This model and numerical results aim to provide some references for bubble dynamics in bioengineering,chemical engineering,naval architecture,and ocean engineering,etc.

bubble,boundary layer,viscosity,surface tension

PACS:47.55.db,47.55.dr,47.55.df,47.11.HjDOI:10.7498/aps.66.234702

*Supported by the National Natural Science Foundation of China(Grant Nos.51639004,51579054,11472088),the Fundamental Research Funds for the Central Universities,China(Grant Nos.HEUCFM170110,HEUCFP201701,HEUCFP201777),and the 111 Project of Harbin Engineering University,China.

?Corresponding author.E-mail:nibaoyu@hrbeu.edu.cn

(2017年7月18日收到;2017年8月12日收到修改稿)

基于氣泡邊界層理論,引入黏性修正,采用邊界積分法,考慮黏性效應(yīng)和表面張力在單氣泡以及雙氣泡耦合作用過(guò)程中的影響.首先將建立的數(shù)值模型與Rayleigh-Plesset的解析解進(jìn)行對(duì)比,發(fā)現(xiàn)二者符合良好,驗(yàn)證了數(shù)值模型的有效性;在此基礎(chǔ)上,建立考慮流體弱黏性效應(yīng)的雙氣泡耦合模型,研究流體黏性和表面張力作用下,氣泡表面變形、射流速度、流場(chǎng)能量轉(zhuǎn)換等物理量的變化規(guī)律;最后研究雷諾數(shù)和韋伯?dāng)?shù)對(duì)于氣泡脈動(dòng)特性的影響規(guī)律.結(jié)果表明,流體黏性會(huì)抑制氣泡脈動(dòng)和氣泡射流發(fā)展,降低氣泡半徑和射流速度;表面張力不改變氣泡脈動(dòng)幅值,但縮短了脈動(dòng)周期,提升氣泡勢(shì)能.

10.7498/aps.66.234702

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):51639004,51579054,11472088)、中央高?；究蒲袠I(yè)務(wù)基金(批準(zhǔn)號(hào):HEUCFM170110,HEUCFP201701,HEUCFP201777)和哈爾濱工程大學(xué)學(xué)科創(chuàng)新引智計(jì)劃資助的課題.

?通信作者.E-mail:nibaoyu@hrbeu.edu.cn