復雜微通道內氣泡在浮力作用下上升行為的格子Boltzmann方法模擬?

婁欽 李濤 楊茉

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

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

1 引 言

氣液兩相流中的氣泡在浮力作用下的上升行為廣泛存在于自然界和生物、化學、工業加工等工程應用中,如:沸騰現象,氣泡在流化床系統內的形成、變形及運動現象[1],醫學領域中微氣泡所產生的動力學行為導致的聲致穿孔現象[2],在冶金領域因化學反應而導致的氣泡形成、長大、運動現象[3],微生物驅油過程中氣泡的產生、變形和移動現象[4],熱交換器中氣泡和液滴的傳熱傳質過程[5]等.因此,開展微通道內氣泡上升行為的研究具有重要的現實意義.

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

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

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

2 數值方法

2.1 格子Boltzmann模型

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

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

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

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

其中ωα代表權重系數.演化方程中(x,t)和(x,t)是分布函數對應的平衡態,其形式如下:

宏觀量指標參數?,壓力p以及速度u的統計由下面方程給出:

流體密度ρ(?)可由指標參數?計算得到,

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

本文使用D2Q9模型來進行數值模擬研究,權重系數ωα設置為:當α=0時,當α =1—4時,當α =5—8時,eα為離散速度,表達式為

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

2.2 潤濕性處理

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

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

3 模型驗證

3.1 拉普拉斯定律

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

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

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

3.2 氣泡在光滑壁面管道內的上升行為

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

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

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

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

圖2 由本文和Alizadeh等(文獻[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 不同條件下氣泡到達微通道終端時的形狀 (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.

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

4 物理問題描述

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

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

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

5 數值結果與分析

5.1 潤濕性的影響

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

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

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

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

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

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

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

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

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

圖8 接觸角θ=150?時對應的氣泡運動瞬時圖(從左到右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).

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

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

5.2 浮力的影響

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

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

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

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

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

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

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

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

5.3 障礙物尺寸的影響

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

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

圖14 障礙物半徑不同時得到的氣泡上升速度演化圖和運動瞬時圖 (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.

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

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

圖15 不同障礙物尺寸下氣泡剩余質量百分比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 氣泡初始位置的影響

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

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

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

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

圖18 氣泡初始放置在左側時氣泡上升速度演化圖和運動瞬時圖(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 不同的氣泡初始位置下氣泡剩余質量百分比Fig.19 .Percentage of the bubble residual mass for dif f erent initial positions of the bubble.

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

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

6 結 論

采用LBM研究了復雜微通道內氣泡在浮力作用下上升的運動特性,主要研究了障礙物表面潤濕性、浮力大小、障礙物尺寸、初始氣泡位置對氣泡上升過程中的變形、分裂、合并等動力學行為以及氣泡瞬時速度、終端速度以及氣泡剩余質量等宏觀特性的影響,得出以下主要結論.

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

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

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

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