婁欽 李濤 楊茉

(上海理工大學(xué)能源與動力工程學(xué)院,上海 200093)

(2018年7月6日收到;2018年9月2日收到修改稿)

1 引 言

氣液兩相流中的氣泡在浮力作用下的上升行為廣泛存在于自然界和生物、化學(xué)、工業(yè)加工等工程應(yīng)用中,如:沸騰現(xiàn)象,氣泡在流化床系統(tǒng)內(nèi)的形成、變形及運(yùn)動現(xiàn)象[1],醫(yī)學(xué)領(lǐng)域中微氣泡所產(chǎn)生的動力學(xué)行為導(dǎo)致的聲致穿孔現(xiàn)象[2],在冶金領(lǐng)域因化學(xué)反應(yīng)而導(dǎo)致的氣泡形成、長大、運(yùn)動現(xiàn)象[3],微生物驅(qū)油過程中氣泡的產(chǎn)生、變形和移動現(xiàn)象[4],熱交換器中氣泡和液滴的傳熱傳質(zhì)過程[5]等.因此,開展微通道內(nèi)氣泡上升行為的研究具有重要的現(xiàn)實(shí)意義.

近幾十年來,許多學(xué)者對微通道內(nèi)的氣泡運(yùn)動特性進(jìn)行了實(shí)驗(yàn)以及數(shù)值研究.早期的研究大都關(guān)注氣泡在光滑壁面微通道內(nèi)的運(yùn)動行為.Bhagat等[6]通過實(shí)驗(yàn)考察了高莫頓數(shù)下氣泡在黏性流體中上升時形狀和終端速度的變化,并發(fā)現(xiàn)當(dāng)雷諾數(shù)小于110時氣泡后部的流線為封閉形態(tài),而當(dāng)雷諾數(shù)大于110時氣泡的運(yùn)動呈現(xiàn)不穩(wěn)定狀態(tài),后部的流線出現(xiàn)開放形態(tài).Weber等[7]在此基礎(chǔ)上,進(jìn)一步研究了不同雷諾數(shù)下上升氣泡的尾部流線軌跡.為了進(jìn)一步刻畫氣泡上升過程中的運(yùn)動特性,Hua等[8]應(yīng)用波前追蹤法(front tracking)研究了高雷諾數(shù)、高密度比和高黏性比情況下,氣泡在黏性流體中的平均上升速度、形狀等特性,并將氣泡的終端速度和終端形狀與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比,發(fā)現(xiàn)氣泡上升速度近似呈對數(shù)形式增長.艾旭鵬等[9]基于邊界層理論,采用邊界積分法,研究了流體黏性效應(yīng)和表面張力對氣泡運(yùn)動特性的影響,發(fā)現(xiàn)流體黏性會使氣泡振動幅值衰減,表面張力的增加會縮短氣泡脈動周期,并提高勢能.Maxworthyt等[10]考查了多個氣泡在浮力作用下的運(yùn)動特性,并給出了不同莫頓數(shù)下氣泡群在上升過程中的阻力系數(shù)和雷諾數(shù)之間的關(guān)系.Rensen等[11]調(diào)查了氣泡群在不同上升速度時的形狀變化,發(fā)現(xiàn)水箱內(nèi)水的深度基本不影響氣泡群的整體形狀.Takada等[12]采用格子Boltzmann自由能模型研究了初始垂直放置的兩個氣泡的變形情況和內(nèi)部流線圖.除了管道內(nèi)的流線軌跡以及氣泡形狀、速度等宏觀特性,還有學(xué)者研究了氣泡運(yùn)動過程中界面的動力學(xué)行為.Sussman等[13]應(yīng)用Level set方法,研究了單個氣泡上升過程中的界面形變以及破裂現(xiàn)象,發(fā)現(xiàn)隨著表面張力和黏性比的增加,氣泡越不容易發(fā)生變形和破裂現(xiàn)象,隨著雷諾數(shù)的增加,氣泡的變形會越來越劇烈.Baltussen等[14]采用流體體積函數(shù)(volume of f l uid,VOF)方法,考察了不同愛特威數(shù)下氣泡上升過程中界面的變化,發(fā)現(xiàn)當(dāng)愛特威數(shù)小于1時,氣泡上升過程中基本不發(fā)生形變.以上文獻(xiàn)主要研究了氣泡初始位置在管道中間時的情形.Fakhari等[15]基于格子Boltzmann自由能模型,研究了氣泡初始緊貼壁面,并沿壁面上升的情況,發(fā)現(xiàn)在高愛特威數(shù)下,氣泡變形加劇,并且在足夠高的阿基米德數(shù)和低莫頓數(shù)時會出現(xiàn)氣泡破裂的情況.

對豎直光滑通道內(nèi)氣泡上升過程中的分裂、合并等現(xiàn)象及其運(yùn)動路徑、速度流場已經(jīng)有了廣泛的研究,為了研究氣泡在復(fù)雜管道內(nèi)的運(yùn)動特性,Uchiyama等[16]通過實(shí)驗(yàn)考察了氣泡群繞豎直通道中心單個圓柱障礙物上升的運(yùn)動狀況,得到了氣泡群演化過程中的速度場和繞圓柱障礙物后的運(yùn)動形狀.Salcedo等[17]研究了通道內(nèi)有兩個相同的圓形障礙物時雷諾數(shù)對氣泡形狀和速度的影響.Li等[18]研究了單個氣泡繞圓形障礙物上升時表面張力和浮力的變化對氣泡形狀的影響.Alizadeh等[19]考察了通道內(nèi)圓形障礙物尺寸和數(shù)量對單個氣泡運(yùn)動特性的影響.除了氣泡在浮力作用下的上升行為,Yi等[20]采用格子Boltzmann自由能模型研究了壓力驅(qū)動下氣泡在含方形障礙物的微通道內(nèi)的上升問題,主要考察了三種特定的壁面潤濕性(接觸角θ=60?,90?,120?)以及不同的氣泡初始位置對氣泡形狀和終端速度的影響.

以上研究工作推動并揭示了氣泡運(yùn)動過程中的一些運(yùn)動機(jī)理,然而還不夠深入,例如微通道表面潤濕性對氣泡的動力學(xué)行為影響較大,但現(xiàn)有文獻(xiàn)對壁面潤濕性的范圍考慮不夠全面.其次,如工業(yè)上換熱器微氣泡強(qiáng)化換熱技術(shù)[21]中通道的表面都不是光滑的,而目前針對壁面上障礙物對氣泡運(yùn)動特性影響的研究非常少.此外,盡管已有研究考慮了氣泡運(yùn)動過程中的變形、破裂現(xiàn)象,但其內(nèi)部機(jī)理尚不清晰,并且氣泡和固體壁面的相互作用以及氣泡穿過通道后的剩余質(zhì)量也鮮有報(bào)道.鑒于此,本文采用格子Boltzmann方法(lattice Boltzmann method,LBM),主要針對氣泡經(jīng)過換熱器圓形管束的上升問題,研究氣泡在含有半圓形障礙物微通道內(nèi)的界面動力學(xué)行為以及宏觀運(yùn)動特性,主要考察障礙物表面潤濕性、浮力大小、障礙物尺寸和氣泡初始位置對氣泡在復(fù)雜微通道內(nèi)受浮力上升問題的影響,探究氣泡變形、分裂等動力學(xué)行為以及氣泡運(yùn)動的剩余質(zhì)量、上升速度和終端速度.

2 數(shù)值方法

2.1 格子Boltzmann模型

LBM在模擬氣液兩相流時有一些獨(dú)特優(yōu)勢,例如可以直觀描述流體和流體之間以及流體和壁面之間的相互作用、不需要追蹤相界面、邊界條件易于處理,因此該方法近幾年受到學(xué)者們的廣泛關(guān)注并得到迅速發(fā)展.目前,基于不同相互作用力的處理方式,已提出多種氣液兩相流模型,主要有顏色模型[22]、偽勢模型[23,24]、自由能模型[25,26]、基于Enskog理論和相場理論的模型[27,28].以上模型都成功地用于研究氣液兩相流動問題.在本文中采用He等[29]基于Enskog理論提出的LBM研究復(fù)雜微通道內(nèi)的氣泡上升問題.

He等[29]提出的兩相流格子Boltzmann模型采用兩個分布函數(shù),fα和gα分別描述指標(biāo)參數(shù)和速度/壓力的演化過程,可以提高數(shù)值穩(wěn)定性.由于該模型已在文獻(xiàn)[29]中有了詳細(xì)的介紹,故本文只做簡要的描述.在此模型中,分布函數(shù)fα和gα分別為

式中 α=0,1,2,···b?1,b為離散速度方向個數(shù);x和t分別表示位置和時間;τ代表松弛時間,與運(yùn)動黏度ν相對應(yīng)的關(guān)系為代表時間步長,是與格子速度c=dx/dt相關(guān)的常數(shù));?,ρ,u分別代表指標(biāo)函數(shù)、流體密度和流體速度;κ代表表面張力強(qiáng)度系數(shù);G為體積力;其中p為流體壓力,演化方程中ψ(?)依賴于模型中用的狀態(tài)方程,本文中采用Carnahan-Starling狀態(tài)方程[30],對應(yīng)的ψ(?)為

其中a決定分子間相互吸引力強(qiáng)度,R為氣體常數(shù),T為流體溫度.在方程(1)中,函數(shù)Γα(u)表達(dá)式為

其中ωα代表權(quán)重系數(shù).演化方程中(x,t)和(x,t)是分布函數(shù)對應(yīng)的平衡態(tài),其形式如下:

宏觀量指標(biāo)參數(shù)?,壓力p以及速度u的統(tǒng)計(jì)由下面方程給出:

流體密度ρ(?)可由指標(biāo)參數(shù)?計(jì)算得到,

其中 ρg和ρl分別代表氣相流體和液相流體密度;?h和?l為指標(biāo)參數(shù)的最大值和最小值,可由Maxwell重構(gòu)得到.采用Chapmann-Enskog多尺度展開可以得到方程(1)對應(yīng)的宏觀動力學(xué)方程[31]:

本文使用D2Q9模型來進(jìn)行數(shù)值模擬研究,權(quán)重系數(shù)ωα設(shè)置為:當(dāng)α=0時,當(dāng)α =1—4時,當(dāng)α =5—8時,eα為離散速度,表達(dá)式為

關(guān)于梯度和拉普拉斯算子的離散方法采用二階中心各向同性方法(ICS)[32]:

2.2 潤濕性處理

壁面潤濕性反映流體和固體之間相互作用力的強(qiáng)度.在微通道內(nèi),潤濕性是影響氣液兩相動態(tài)行為的重要參數(shù).在LBM中,壁面的潤濕性通過邊界條件來描述.本文將考慮潤濕性對氣泡通過障礙物通道的動態(tài)影響,潤濕性邊界條件采用Davies等提出的方法[33],在該方法中采用表面親和性αs刻畫壁面的潤濕強(qiáng)度,并把表面親和性與指標(biāo)參數(shù)?聯(lián)系起來,其關(guān)系式為

其中σs1,σs2分別代表固-液表面張力和固-氣表面張力;σ12為氣-液表面張力.

3 模型驗(yàn)證

3.1 拉普拉斯定律

首先采用Laplace定律來驗(yàn)證模型的正確性.初始時在Lx×Ly區(qū)域中心內(nèi)放置半徑為R,密度為ρl的靜止圓形液滴,其余區(qū)域充滿密度為ρg的氣體.根據(jù)Laplace定律可知,當(dāng)系統(tǒng)達(dá)到穩(wěn)定時表面張力σ恒定,且液滴內(nèi)外的壓力差?P與半徑的倒數(shù)1/R呈線性關(guān)系,關(guān)系式為

本文的驗(yàn)證中,初始液滴半徑R分別取24,28,32,36,40五種情況,并考慮不同表面張力強(qiáng)度κ=0.1,0.15,0.2,以保證驗(yàn)證可信度.在數(shù)值模擬中網(wǎng)格數(shù)均取為128×128,狀態(tài)方程中分子間強(qiáng)度a為4.0,RT的取值為1/3,對應(yīng)的?h和?l分別為0.250291和0.022838,其余參數(shù)設(shè)置為:ρl=0.5,ρg=0.1,ν=1/6.計(jì)算域四周的邊界條件皆為周期性邊界條件.數(shù)值模擬達(dá)到穩(wěn)態(tài)的條件為

圖1 液滴內(nèi)外壓力差pi?po和半徑倒數(shù)1/R之間的關(guān)系Fig.1 . Relationship between pressure jump across the droplet interface pi?p0and inverse of droplet radius 1/R.

3.2 氣泡在光滑壁面管道內(nèi)的上升行為

氣泡在光滑壁面通道內(nèi)受浮力上升的問題與本文所要研究的問題較貼近,且已有學(xué)者對該問題進(jìn)行了大量研究,因此本章節(jié)將采用該算例驗(yàn)證模型的正確性.對于該問題,初始時在Lx×Ly的計(jì)算區(qū)域內(nèi)放置半徑為R密度為ρg的靜止圓形氣泡,氣泡的圓心設(shè)為(xc,yc),其余區(qū)域充滿密度為ρl的液體,在豎直方向施加浮力G=(ρ?ρl)g,其中g(shù)為重力加速度.在數(shù)值模擬中,Lx×Ly=80×300,R=10,ρl=1.48,ρg=0.52,κ =0.15,xc為微通道水平方向的中心位置,yc為通道高度的1/4,計(jì)算域的邊界條件設(shè)置為:左右固體壁面采用無滑移邊界條件,進(jìn)出口選用周期性邊界條件.

該問題有兩個重要無量綱數(shù),E?tv?s數(shù)(Eo)和Morton數(shù)(Mo),其中Eo=g(ρl?ρg)d2/σ代表浮力與表面張力的比值;Mo=g(ρl? ρg)μ4l/ρ2lσ3代表黏性力和表面張力的比值,式中d為氣泡直徑,μl為液體的動力黏度.表1給出了不同Eo數(shù)和Mo數(shù)時氣泡的終端速度Ut,并與LBM方法中基于Enskog理論和相場理論的模型[19]以及VOF方法[34]進(jìn)行對比,發(fā)現(xiàn)結(jié)果基本一致.

表1 不同模型對應(yīng)的氣泡終端速度Table 1 .Terminal velocities of the bubble obtained by dif f erent models.

圖2給出了表1中當(dāng)Eo=10,Mo=0.4535以及Eo=40,Mo=1.813(對應(yīng)(b),(d)兩種情況)時對應(yīng)的氣泡上升的速度變化趨勢,并與Alizadeh等[19]的數(shù)據(jù)進(jìn)行對比.從圖2中可以看出,本文得到的速度變化趨勢與Alizadeh等得到的結(jié)果基本一致.具體地說,當(dāng)Eo數(shù)較小時(Eo=10),演化初期氣泡速度在浮力作用下迅速增加,隨著浮力、表面張力以及黏性力達(dá)到平衡,氣泡速度基本穩(wěn)定.當(dāng)Eo數(shù)較大時(Eo=40),演化初期氣泡上升速度的增加速率更快,達(dá)到的極值更大.當(dāng)氣泡上升速度達(dá)到最大值時,由于氣泡變形更加明顯,因此會在表面張力的作用下,出現(xiàn)速度值略微下降的現(xiàn)象[19],此后各項(xiàng)力之間達(dá)到平衡,氣泡上升速度近似不變.

圖2 由本文和Alizadeh等(文獻(xiàn)[19])得到的氣泡上升速度 (a)Eo=10;(b)Eo=40Fig.2 .Bubble rising velocities obtained by the present study and by Alizadeh et al.(Ref.[19]):(a)Eo=10;(b)Eo=40.

圖3 不同條件下氣泡到達(dá)微通道終端時的形狀 (a)Eo=5;(b)Eo=10;(c)Eo=20;(d)Eo=40Fig.3 .Bubble shapes at the micro-channel terminal under dif f erent conditions:(a)Eo=5;(b)Eo=10;(c)Eo=20;(d)Eo=40.

為了進(jìn)一步驗(yàn)證本文模型,圖3展示了不同Eo數(shù)下,氣泡到達(dá)通道終端時的形狀.從圖3中可以看出,當(dāng)Eo=5時,氣泡上升過程中基本保持初始圓形,而隨著Eo數(shù)的增加,氣泡在上升過程中會出現(xiàn)更加明顯的形變.此現(xiàn)象與文獻(xiàn)[12,35]得到的結(jié)果一致.

4 物理問題描述

本文研究的物理問題如圖4所示.微通道寬度和高度分別為Lx和Ly,微通道內(nèi)黑色半圓區(qū)域代表障礙物,其半徑為R1,藍(lán)色圓形區(qū)域?yàn)槌跏荚谕ǖ纼?nèi)的氣泡,其密度為ρg,半徑為R2,圓心為(xc,yc),其中xc為微通道水平方向的中心位置,yc為通道高度的五分之一,其余白色區(qū)域充滿密度為ρl的液體.在流體流動的y方向施加浮力G=(ρ?ρl)g.該問題采用的邊界條件與3.2節(jié)相同.

圖4 模擬問題示意圖Fig.4 . illustration of simulation problem.

在下文分析中,為了說明氣泡的運(yùn)動特性,引入兩個參數(shù):剩余質(zhì)量百分比Rq和t時刻氣泡在y方向的上升速度vb(t).剩余質(zhì)量百分比為到達(dá)微通道終端的氣泡質(zhì)量占初始?xì)怏w質(zhì)量的百分比,這里需要指出的是到達(dá)微通道終端的氣泡質(zhì)量即剩余氣泡質(zhì)量可以由初始?xì)馀菘傎|(zhì)量減去殘留在障礙物表面的氣相質(zhì)量得到,也可由氣泡穿過障礙物到達(dá)微通道終端時所占據(jù)的氣相區(qū)密度求和得到,其中氣相區(qū)為密度ρ小于0.5·(ρl+ρg)的區(qū)域;而t時刻氣泡上升√速度定義為氣泡豎直方向的上升速度其中x代表t時刻被氣泡占據(jù)的格點(diǎn)位置,v(x,t)代表t時刻氣泡在x點(diǎn)垂直方向上的速度,N代表此時刻被氣泡占據(jù)的所有格點(diǎn)數(shù).在本文中取為特征時間對時間進(jìn)行無量綱化,在下文中無量綱時間用T表示.

5 數(shù)值結(jié)果與分析

5.1 潤濕性的影響

本小節(jié)主要研究障礙物表面潤濕性不同時的氣泡運(yùn)動特性.數(shù)值模擬中分別考慮接觸角θ=30?,40?,50?,60?,70?,80?,82?,85?,90?,100?,110?,120?,130?,140?,150?的情況,其他參數(shù)設(shè)置如下:Lx和Ly分別設(shè)置為160和600格子單位,Eo=10,Mo=0.4535,R1=60,R2=32,ρl=0.5,ρg=0.1,μl=0.05.

圖5給出了接觸角θ=60?時氣泡上升速度演化圖和對應(yīng)速度極值的氣泡運(yùn)動瞬時圖.從圖5中可以看出,在氣泡的運(yùn)動過程中,速度呈波浪形變化,依次出現(xiàn)增加-減小-增加-減小-增加-減小-基本不變幾個階段(與文獻(xiàn)[20]中氣泡穿過方形障礙物通道時得到的結(jié)果類似),對應(yīng)兩個極小值和三個極大值,這些速度極值點(diǎn)均發(fā)生在氣泡穿過通道較窄部分的過程中.具體而言,演化過程初期,在浮力作用下,氣泡速度逐漸增加,在T=1.40時,氣泡靠近障礙物表面,通道截面積減小,氣泡運(yùn)動受到阻礙,其速度開始減小.速度減小趨勢一直持續(xù)到T=3.76,此時氣泡頂部剛剛穿過通道最窄部分,隨后氣泡受到的阻礙力減小,氣泡上升速度開始增加.在T=6.86時,氣泡底部也穿過了通道最窄部分,隨著通道截面積的增加,氣泡在表面張力作用下變形,氣泡上升速度開始減小,直至氣泡前端接近球形(此時T=7.82),之后氣泡形狀變化不明顯,其速度在浮力作用下增加,而在T=8.70時,氣泡與障礙物脫離(速度對應(yīng)第三個極大值).由于氣泡脫離障礙物時,與障礙物表面接觸部分生成兩個較長的“小觸須”,發(fā)生嚴(yán)重變形,導(dǎo)致氣泡與障礙物脫離之后在表面張力作用下收縮,其速度出現(xiàn)下降趨勢.隨著氣泡形狀變化越來越小,氣泡上升速度趨于穩(wěn)定(T=11.43之后),此情況與文獻(xiàn)[8]中氣泡受浮力上升、速度最終會趨于穩(wěn)定的情況類似.

圖5 接觸角θ=60?時對應(yīng)的氣泡上升速度演化圖和運(yùn)動瞬時圖Fig.5 .Time histories of the bubble velocity and the snapshots of the bubble at extreme positions under θ =60?.

與氣泡的速度軌跡對應(yīng),其形狀也發(fā)生了巨大的變化.由于浮力的作用,初始圓形氣泡在上升速度出現(xiàn)第一次極大值時變成了下凹的球冠形.此時氣泡除了受到向上的浮力作用外,還受到障礙物的擠壓力,因此氣泡發(fā)生變形以穿過通道較窄部分.氣泡在穿過通道較窄部分過程中,由于障礙物表面親水而疏氣,與障礙物表面始終存在明顯的間隔,并且在穿過通道時與障礙物表面接觸面積較少.一旦大部分氣泡穿過障礙物之間的通道,即快速與半圓形障礙物脫離,恢復(fù)成下凹的球冠形穩(wěn)定上升.

為了進(jìn)一步說明氣泡通道障礙物表面潤濕性對氣泡動力學(xué)行為的影響,圖6給出了障礙物表面接觸角為120?的氣泡上升速度演化圖和對應(yīng)速度極值的運(yùn)動瞬時圖.對比圖5和圖6可以發(fā)現(xiàn),當(dāng)障礙物表面潤濕性不同時,氣泡上升速度都呈波浪形變化,速度極值點(diǎn)均出現(xiàn)在氣泡通過障礙物的過程中,而且在氣泡接觸障礙物之前,不同時刻氣泡上升速度的值基本一致.因此,這兩種情況下,氣泡都是在T=3.76時刻擠入通道最窄處.然而,隨著障礙物表面接觸角增加,其對氣泡的黏附力增加,氣泡通過障礙物通道的上升速度減小,所需要的時間增加.例如,障礙物表面接觸角θ=120?的工況下,氣泡從開始擠入通道最窄處到穿過障礙物通道的時間為T=3.76至T=8.55,比θ=60?時增加了17.98%;氣泡由障礙物上部到完全脫離障礙物表面的時間為T=8.55到T=11.06,比θ=60?時增加了185.22%.完全脫離后的氣泡上升速度變化趨勢與θ=60?的情況類似,都是先因氣泡形狀發(fā)生較大的變化,對應(yīng)的速度有所下降,再隨著氣泡形狀變化越來越小,氣泡上升速度趨于穩(wěn)定.氣泡到達(dá)微通道終端的總時間T=17.70,相比于θ=60?的情況增加了14.34%.

圖6 接觸角θ=120?時對應(yīng)的氣泡上升速度演化圖和運(yùn)動瞬時圖Fig.6 .Time histories of the bubble velocity and the snapshots of the bubble at extreme positions under θ =120?.

從以上結(jié)果可以看出,隨著接觸角的增加,障礙物對氣泡的黏附力增加,于是氣泡穿過障礙物通道所需要的時間更長,困難性增加,有可能造成到達(dá)通道終端的速度減小.為了驗(yàn)證這一理論預(yù)測,圖7給出了氣泡終端速度和通道內(nèi)障礙物表面潤濕性之間的關(guān)系.從圖7中可以看出,隨著障礙物表面接觸角增加,氣泡的終端速度減小,該結(jié)果驗(yàn)證了理論預(yù)測.

圖7 不同接觸角下氣泡終端速度Fig.7 .Terminal velocities of the bubble under dif f erent values of contact angle.

此外,對比圖5和圖6還可以發(fā)現(xiàn),當(dāng)接觸角較大時,由于障礙物對氣泡的黏附力增加,氣泡在通過障礙物之間的通道時,與障礙物表面基本貼合,氣泡形狀變化出現(xiàn)較大不同.例如,在障礙物表面接觸角θ=120?時,氣泡底部通過通道最窄部分時變形嚴(yán)重,底部兩側(cè)分別出現(xiàn)小圓孔;氣泡整體通過障礙物通道并將要脫離障礙物表面時,其底部被拉伸,出現(xiàn)“拱門”形狀.為了更清晰地說明障礙物表面潤濕性對氣泡形狀的影響,圖8給出了障礙物表面接觸角為150?時的氣泡運(yùn)動瞬時圖.從圖8中可以看出,由于接觸角的進(jìn)一步增加,氣泡受到的黏附力增大,在通過障礙物之間的通道時,其兩側(cè)與障礙物表面完全接觸,孔隙消失;在氣泡脫離障礙物過程中,其“拱門”形狀向兩邊拉伸更明顯.對比圖5、圖6和圖8可以發(fā)現(xiàn),隨著障礙物表面接觸角的增加,氣泡完全脫離障礙物后剩余的質(zhì)量更少,說明有部分氣泡殘留在通道障礙物表面上.

圖8 接觸角θ=150?時對應(yīng)的氣泡運(yùn)動瞬時圖(從左到右T=3.61,6.49,8.92,12.09)Fig.8 .Snapshots of the bubble at θ =150? (from left to right T=3.61,6.49,8.92,12.09).

為了定性刻畫障礙物表面潤濕性對氣泡運(yùn)動剩余質(zhì)量的影響,圖9給出了不同接觸角時氣泡剩余質(zhì)量百分比Rq.從圖9中可以看出,當(dāng)障礙物表面接觸角θ小于或等于82?時,氣泡的剩余質(zhì)量百分比Rq均為100%,說明氣泡上升過程中雖然變形嚴(yán)重,但由于障礙物表面疏氣,氣泡不會黏附在障礙物表面上,而是完整地到達(dá)微通道終端.而當(dāng)障礙物表面接觸角θ等于85?時,如圖9中氣泡運(yùn)動瞬時圖所示,在氣泡通過障礙物通道時,與障礙物表面兩側(cè)緊密接觸,并受到黏附力作用,部分殘留在障礙物表面上,氣泡剩余質(zhì)量百分比為87.37%.隨著障礙物表面接觸角越來越大,其對氣泡的黏附力增大,則殘留在障礙物表面上的氣泡質(zhì)量越多,氣泡剩余質(zhì)量百分比Rq明顯下降.

圖9 不同接觸角下氣泡剩余質(zhì)量百分比Fig.9 .Percentage of the bubble residual mass under dif f erent values of contact angle.

5.2 浮力的影響

浮力作為驅(qū)使氣泡在微通道內(nèi)上升的作用力,具有重要的研究意義,其中上文提到的Eo數(shù)表征浮力與表面張力的相對大小,本節(jié)在維持初始表面張力不變的情況下,通過調(diào)節(jié)Eo數(shù)的值來得到不同的氣泡運(yùn)動狀態(tài),即研究浮力對氣泡在微通道內(nèi)上升的影響.在本文模擬中,分別考慮了Eo數(shù)為10,15,20,25,30,35和40七種情況.在所有情況中,θ=90?,其他參數(shù)設(shè)置與5.1小節(jié)中相同.

圖10為Eo=10時,氣泡上升速度演化圖和對應(yīng)速度極值的運(yùn)動瞬時圖.從圖10中可以看出,在演化初期,氣泡上升速度在浮力的作用下逐漸增加,當(dāng)氣泡靠近障礙物時,受到障礙物的阻礙作用,其上升速度開始減小.隨后氣泡頂部擠過通道最窄部分其速度又逐漸增加,當(dāng)所有氣泡都穿過障礙物通道最窄部分時,由于通道截面的增加,其在表面張力作用下開始變形,對應(yīng)的上升速度減小.隨著氣泡的形狀趨近于圓形,表面張力的影響減弱,上升速度再次增加,氣泡發(fā)生分裂現(xiàn)象,一部分殘留在障礙物表面上,另一部分繼續(xù)上升且形狀逐漸恢復(fù)成下凹的球冠形,出現(xiàn)速度下降現(xiàn)象,此后氣泡形狀變化越來越小,氣泡上升速度趨于穩(wěn)定.該現(xiàn)象與上一小節(jié)現(xiàn)象一致.

圖10 Eo=10時對應(yīng)的氣泡上升速度演化圖和運(yùn)動瞬時圖Fig.10 .Time histories of the bubble velocity and the snapshots of the bubble at extreme positions under Eo=10.

圖11 Eo=30時對應(yīng)的氣泡上升速度演化圖和運(yùn)動瞬時圖Fig.11 .Time histories of the bubble velocity and the snapshots of the bubble at extreme positions under Eo=30.

圖11為Eo=30時,氣泡上升速度演化圖和對應(yīng)速度極值的運(yùn)動瞬時圖.對比圖10和圖11可以發(fā)現(xiàn),Eo數(shù)的增大造成氣泡上升速度變化趨勢以及界面動力學(xué)行為出現(xiàn)顯著的不同.一方面,當(dāng)Eo=10時,氣泡的上升速度相對較低,在運(yùn)動過程中出現(xiàn)三個極大值、兩個極小值的情況,而當(dāng)Eo=30時,浮力相對于表面張力增加,氣泡穿過障礙物通道的過程中與障礙物有明顯的間隙,更容易通過障礙物通道,而且由于表面張力作用減弱,氣泡通過障礙物通道后不再出現(xiàn)因?yàn)樽冃味鴮?dǎo)致速度減小的趨勢,在運(yùn)動過程中速度只出現(xiàn)兩個極大值,一個極小值.Eo數(shù)的增大使得氣泡在運(yùn)動過程中整體上升速度增加,到達(dá)微通道終端所需時間減少,當(dāng)Eo=30時,氣泡會在T=13.89時刻到達(dá)通道終端,而對應(yīng)Eo=10時,氣泡會在T=16.52時刻到達(dá)通道終端.另一方面,Eo數(shù)不同時氣泡的界面動力學(xué)行為也有較大不同.當(dāng)Eo=10時,氣泡緊貼障礙物壁面上升,在大部分氣泡通過通道最窄處時,氣泡底部形成兩個“小圓孔”(T=7.01),此后氣泡發(fā)生嚴(yán)重變形,逐漸形成一個類似“拱門”的形狀.隨后,在浮力作用下氣泡發(fā)生分裂行為(T=9.59),一部分殘留在障礙物表面上,另一部分則繼續(xù)上升,完全脫離后的氣泡在浮力作用下開始變形并逐漸恢復(fù)成下凹的球冠形.當(dāng)Eo=30,氣泡逐漸擠入障礙物通道時(T=3.56),受到的浮力更大,氣泡下半部分變形更加劇烈,變形的弧度更明顯,而且在氣泡通過障礙物通道的過程中,上半部分基本不會與障礙物有接觸,而下半部分在較大的浮力和一定的障礙物擠壓力下,生成兩個較長的“小觸須”.隨著時間演化,有一部分“小觸須”穿過障礙物表面,另一部分殘留在障礙物表面上,穿過的“小觸須”在表面張力的作用下,重新形成球冠狀隨大氣泡一起上升,初始兩者并未融合,隨著時間演化由“小觸須”形成的那部分追上了上部分大氣泡,兩氣泡合并后一起上升.

圖12給出了不同Eo數(shù)下,氣泡經(jīng)障礙物通道到達(dá)微通道終端的剩余質(zhì)量百分比Rq.從圖12中可以看出,剩余質(zhì)量百分比Rq隨Eo數(shù)增大基本呈對數(shù)形式增加.當(dāng)Eo=10時,剩余質(zhì)量百分比Rq=73.26%,在此情況下,由于浮力相對于表面張力給予氣泡的作用力并不明顯,氣泡在經(jīng)過障礙物表面時,在浮力、表面張力和障礙物阻力共同作用下,發(fā)生破裂,破裂的氣泡會殘留在障礙物表面上.隨著浮力逐漸增加,開始在整個力場占據(jù)主導(dǎo)地位,氣泡能相對順利地通過障礙物通道,從而破裂后殘留在障礙物表面的氣泡越來越少,其剩余質(zhì)量百分比Rq增加.當(dāng)浮力增加到一定值(Eo=30)時,氣泡整體基本能完整地通過障礙物通道,氣泡剩余質(zhì)量百分比Rq基本保持不變.圖13給出了不同Eo數(shù)對氣泡終端速度vb的影響.從13圖中可以看出,隨著Eo數(shù)的增加,氣泡終端速度vb的變化也基本呈現(xiàn)對數(shù)形式增加.當(dāng)Eo數(shù)較小時,氣泡到達(dá)微通道終端的終端速度vb隨著Eo數(shù)增加得較快,當(dāng)增加到一定值(Eo=30)時,增長速度逐漸放緩.

圖12 不同Eo數(shù)下氣泡剩余質(zhì)量百分比Fig.12 .Percentage of the bubble residual mass obtained under dif f erent Eo numbers.

圖13 不同Eo數(shù)下氣泡的終端速度Fig.13 .Terminal velocities of the bubble under different Eo numbers.

5.3 障礙物尺寸的影響

障礙物作為微通道結(jié)構(gòu)中重要的一環(huán),其尺寸的改變會對氣泡運(yùn)動造成重要的影響,本節(jié)分析不同障礙物半徑大小(R1=50,55,60,65,70)對氣泡運(yùn)動的影響.數(shù)值模擬中Eo=15,Mo=0.6802,θ=90?,其余參數(shù)設(shè)置與5.1小節(jié)相同.

圖14給出了障礙物半徑為60和70時,氣泡上升速度演化圖和對應(yīng)速度極值的運(yùn)動瞬時圖.從圖14中可以看出,在氣泡未到達(dá)障礙物通道前,氣泡的上升速度變化趨勢基本相同;在進(jìn)入障礙物通道后,由于障礙物尺寸不同,氣泡上升速度的值出現(xiàn)明顯的不同.例如當(dāng)R1=60時,氣泡在通過障礙物通道過程中,受到障礙物的阻礙作用,上升速度由0.0403下降到0.0342;當(dāng)R1=70時,氣泡受到的阻礙更大,氣泡上升速度由0.0278下降到0.0108.在氣泡脫離障礙物表面并上升到微通道終端的過程中,不同障礙物半徑對應(yīng)的氣泡上升速度變化趨勢也不相同,比如當(dāng)R1=60,氣泡整體脫離障礙物時,氣泡上升速度略微下降后出現(xiàn)振蕩,之后因形狀劇烈變化造成速度持續(xù)下降,直到形狀穩(wěn)定后(圖14(a)T=10.38),速度逐漸趨于穩(wěn)定;當(dāng)R1=70時,由于氣泡通過障礙物通道時受到的阻力較大,則氣泡與障礙物脫離時上升速度較小,隨后在浮力的作用下,速度略有上升,再因表面張力作用產(chǎn)生劇烈變形,變成下凹的球冠形,速度略微下降,再度趨于平穩(wěn)(圖14(b)中T=16.56時).障礙物半徑的增加使得氣泡通過障礙物通道更加困難,到達(dá)微通道終端的時間增加.當(dāng)半徑R1分別為60和70時,氣泡到達(dá)微通道終端的時間分別為T=16.17和T=24.84.

圖14 障礙物半徑不同時得到的氣泡上升速度演化圖和運(yùn)動瞬時圖 (a)R1=60;(b)R1=70Fig.14 .Time histories of the bubble velocity and the snapshots of the bubble under dif f erent values of radius:(a)R1=60;(b)R1=70.

相對應(yīng)的,氣泡形狀的變化也會因障礙物尺寸不同而出現(xiàn)變化.在氣泡擠入障礙物之間的通道后,半徑為60的情況中,氣泡在通過障礙物通道時,與障礙物表面之間有一定的間隙,氣泡能較完整地通過障礙物通道,并以半圓形狀繼續(xù)上升;而半徑為70的情況中,氣泡完全緊貼障礙物表面上升,隨后會出現(xiàn)“拱門”形狀,且當(dāng)氣泡通過障礙物通道后,有部分氣泡殘留在障礙物表面,該情形下脫離后的氣泡明顯更小.

為進(jìn)一步說明障礙物尺寸對氣泡運(yùn)動的影響,圖15和圖16展示了不同障礙物半徑下,氣泡剩余質(zhì)量百分比Rq以及氣泡終端速度vb的變化趨勢.從圖15中可以看出,氣泡剩余質(zhì)量百分比Rq的變化趨勢表現(xiàn)為初期下降較慢而后期下降較快,這是因?yàn)檎系K物半徑較小時,大部分氣泡能較順利地通過障礙物通道,而隨著障礙物半徑增加到60以上時,氣泡運(yùn)動過程中受到的阻力增加,有相當(dāng)一部分氣泡破裂后殘留在障礙物表面上,剩余質(zhì)量百分比Rq隨障礙物半徑增加而迅速下降.從圖16中可以看出,氣泡終端速度隨著障礙物半徑的增加而近似呈線性減小,這是由于障礙物半徑的增加,使得氣泡通過障礙物通道更加困難,產(chǎn)生的變形更嚴(yán)重,從而終端速度降低.

圖15 不同障礙物尺寸下氣泡剩余質(zhì)量百分比Fig.15 .Percentage of the bubble residual mass obtained for dif f erent sizes of the obstacle.

圖16 不同障礙物尺寸下氣泡終端速度Fig.16 .Terminal velocities of the bubble for dif f erent sizes of the obstacle.

5.4 氣泡初始位置的影響

本小節(jié)研究氣泡初始位置對氣泡運(yùn)動情況的影響.如圖17所示,初始?xì)馀輬A心位置為(xd,yc),xd為距離左壁面的R2位置,yc則仍為微通道高度的五分之一,在數(shù)值模擬中考慮不同的Eo數(shù)情況:Eo=10,15,20,25,30,35和40.在所有情況中,θ=90?,其他參數(shù)以及邊界條件設(shè)置同5.1小節(jié).

圖17 模擬問題對比示意圖Fig.17 .Comparison of simulation problem.

圖18給出了氣泡初始位置圓心為(xd,yc)情況下的氣泡上升速度變化趨勢以及對應(yīng)速度極值的運(yùn)動瞬時圖(對應(yīng)Eo=10).與圖10相比,可以發(fā)現(xiàn)氣泡初始放置在一側(cè)時對應(yīng)的上升速度變化趨勢與氣泡放置在管道中間時基本一致,同樣出現(xiàn)三個極大值點(diǎn)和兩個極小值點(diǎn),但整體的上升速度均小于初始放置在通道中間的氣泡上升速度,此情況與文獻(xiàn)[20]給出的不同氣泡初始位置的速度變化趨勢類似.其中,氣泡初始受浮力上升并即將擠入障礙物之間的通道時,初始?xì)馀莘胖迷谝粋?cè)的情況中,氣泡受到障礙物(尤其是左側(cè)障礙物)的阻力更大,因此上升速度相對較低,氣泡放置在一側(cè)對應(yīng)的上升速度值為0.0116,而氣泡放置在管道中間對應(yīng)的上升速度值為0.0165.相應(yīng)地,在通過障礙物通道的過程中,初始?xì)馀莘胖迷谝粋?cè)的情況下,對應(yīng)氣泡受到的阻力更大,上升速度由0.0209下降到0.0073,而氣泡放置在管道中間的情況中,上升速度由0.0285下降到0.0182.兩者到達(dá)微通道終端的終端速度也存在明顯差異,氣泡放置在一側(cè)對應(yīng)的終端速度為0.0128,氣泡放置在中間對應(yīng)的終端速度為0.0179.

從圖18中還可以看出,演化初期由于氣泡受力不平衡,在進(jìn)入通道內(nèi)最狹窄部分時,氣泡已經(jīng)發(fā)生了破裂現(xiàn)象,且在T=8.41時有少部分氣泡會殘留在微通道的左側(cè)障礙物表面上.由于氣泡進(jìn)入障礙物通道時,左側(cè)的障礙物與氣泡的接觸面積始終較大,在通過障礙物通道的過程中,左側(cè)障礙物的阻力作用更明顯,形成了整體偏左的“拱門”形狀(此時T=11.72).隨后在浮力和表面張力的共同作用下,氣泡逐漸與障礙物脫離,有相當(dāng)一部分氣泡殘留在障礙物表面上(此時T=15.19),而脫離的氣泡繼續(xù)上升,其形狀慢慢恢復(fù)成為向左下方傾斜的球冠形,最終移動到微通道水平中心線位置并沿著水平中心線繼續(xù)上升,但由于脫離氣泡在擠出障礙物時,形狀不對稱,因此氣泡到達(dá)微通道終端時仍會出現(xiàn)左右不對稱現(xiàn)象.由此情況可以看到氣泡與障礙物表面接觸更普遍,障礙物的阻力作用更加明顯,因此氣泡變形更嚴(yán)重,上升過程更加困難,破裂并殘留在障礙物表面的氣泡更多.

圖18 氣泡初始放置在左側(cè)時氣泡上升速度演化圖和運(yùn)動瞬時圖(Eo=10)Fig.18 .Time histories of the bubble velocity and the snapshots of the bubble at extreme positions when the bubble is initially set at the left side of the channel(Eo=10).

圖19 不同的氣泡初始位置下氣泡剩余質(zhì)量百分比Fig.19 .Percentage of the bubble residual mass for dif f erent initial positions of the bubble.

圖19給出了兩種氣泡初始位置下,剩余質(zhì)量百分比Rq隨Eo數(shù)的變化趨勢.從圖19中可以發(fā)現(xiàn),初始?xì)馀菸恢貌煌瑫r,其剩余質(zhì)量百分比Rq整體趨勢類似,Rq首先隨Eo數(shù)增大而增加,然后再趨于穩(wěn)定,但當(dāng)氣泡初始放置在通道一側(cè)時,剩余質(zhì)量百分比Rq相對較小.發(fā)生該現(xiàn)象是由于氣泡初始放置在一側(cè)時,上升過程中受力嚴(yán)重不平衡,因此變形現(xiàn)象更嚴(yán)重,表現(xiàn)為氣泡上升過程中基本需完全接觸左側(cè)半圓形障礙,受到障礙物的阻力更明顯,殘留更多的氣泡在障礙物表面,因此氣泡的剩余質(zhì)量百分比Rq相對較小.圖20給出了兩種不同氣泡初始位置下,氣泡終端速度vb隨Eo數(shù)的變化趨勢.從圖20中可以清晰地看到,氣泡初始位置不同時氣泡終端速度vb的趨勢基本一致,但氣泡初始放置在一側(cè)的情況中,氣泡在上升過程中受到障礙物阻力更大,因此終端速度vb相對較小,與文獻(xiàn)[20]中得到的結(jié)論一致.

圖20 不同的氣泡初始位置下氣泡終端速度Fig.20 .Terminal velocities of the bubble for dif f erent initial positions of the bubble.

6 結(jié) 論

采用LBM研究了復(fù)雜微通道內(nèi)氣泡在浮力作用下上升的運(yùn)動特性,主要研究了障礙物表面潤濕性、浮力大小、障礙物尺寸、初始?xì)馀菸恢脤馀萆仙^程中的變形、分裂、合并等動力學(xué)行為以及氣泡瞬時速度、終端速度以及氣泡剩余質(zhì)量等宏觀特性的影響,得出以下主要結(jié)論.

1)障礙物表面潤濕性表現(xiàn)為親水時,氣泡基本可以完整地通過障礙物通道;障礙物表面潤濕性表現(xiàn)為中性或者親氣時,氣泡通過障礙物通道時會發(fā)生分裂行為并有部分氣泡殘留在障礙物表面上,且隨著障礙物表面潤濕性的增加,氣泡通過障礙物通道時變形越嚴(yán)重,殘留氣泡質(zhì)量越大,同時氣泡的上升速度和終端速度越小.

2)隨著浮力增加,通過微通道的氣泡質(zhì)量和氣泡終端速度均近似呈對數(shù)形式增長.

3)隨著障礙物半徑增加,氣泡殘留在其表面上的質(zhì)量增加,趨勢表現(xiàn)為初期增加較慢而后期增加較快,氣泡終端速度隨著障礙物半徑的增加而近似線性減小.

4)當(dāng)氣泡的初始位置偏離管道中心時,氣泡在上升過程中受力嚴(yán)重不平衡,因此氣泡變形更明顯,且此時當(dāng)氣泡穿過障礙物時,殘留在障礙物表面以及氣泡與障礙物表面接觸面上的質(zhì)量更多,氣泡上升速度以及終端速度更小.