含孔板微通道內隨流氣泡運動的兩相格子Boltzmann模擬

丁弘毅, 王治云, 趙敬帥, 王 楠, 婁 欽

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

0 引 言

隨著能源利用的發展,含氣泡運動的兩相流系統在自然界和工程應用中廣泛存在,如艦船廢氣隱蔽排放、頁巖油開發、生物柴油制備、CO2地質封存等[1-2],這些過程都存在氣泡通過孔板的現象．撞擊孔板過程中,氣泡在表面張力和慣性力的作用下發生形變甚至破碎現象,涉及復雜的界面演化問題[3]．隨著計算科學的快速發展,此類流動的數值模擬能揭示氣泡破裂聚合的機理,為提高兩相流傳熱傳質效率提供理論依據,是實驗研究的重要補充[4],具有重要的工程實踐和學術研究意義[5]．

許多研究者都對微通道內的氣泡流動問題進行過實驗研究．Hallmark等[6]研究了含孔板矩形通道內,氣泡在高黏度Newton流體中上升時的形狀變化,發現氣泡穿越孔板后演化為“新月形”．Dawson等[7]研究了體積流量恒定條件下氣泡在狹窄微通道內傳播變形的過程,發現該過程由黏性剪切力主導,從而影響氣泡前后部的曲率變化規律．Losi等[8]對氣泡在水平管道中的運動進行了深入實驗,分析了黏度對氣泡在水平管道中漂移速度和形態的影響．Wu等[9]研究了孔喉微結構中的氣泡破碎現象,實驗表明隨著孔徑比、喉長、母泡速度和毛細數的增大,子泡大小呈線性或指數減小趨勢．

過去三十余年,作為一種介觀數值計算方法,格子Boltzmann方法(lattice Boltzmann method, LBM)在模擬多相流時具有算法簡單、并行效率高和復雜邊界易于實現等顯著優勢[10],因此開發能夠模擬大密度比流場的多相流模型在LBM研究者中是一個很有吸引力的課題[1]．Inamuro等[11]提出了第一個基于自由能方法的LB模型,實現了密度比高達1 000的液滴與氣泡模擬計算．隨后Yi和Xing[12]使用其改進模型,模擬了光滑毛細管和窄喉管道中的氣泡-水流動,系統分析了煤的潤濕性對切割兩相流動的影響．Sankaranarayanan等[13]運用多組分偽勢LB方法建立了氣泡流動的LB基準模型,成功模擬了單個氣泡和規則氣泡群的上升過程,其結果與實驗數據一致．Yuan和Schaefer[14]將偽勢模型推廣到大密度比的情況,但應用于動態問題時會受到限制．部分研究者也嘗試建立基于相場理論的LB模型[15-16],該模型可以很好地進行界面捕捉,在多相流模擬計算中越來越流行[17]．Lou等[18]基于大密度比LB模型,對氣泡在不同類型多孔介質中的上升過程進行數值模擬,研究了不同Eotvos數、黏度比、孔隙形狀、氣泡數量及不同多孔介質排列方式下的氣泡和流體動力學．

以上研究揭示了氣泡在復雜結構微通道內的運動機理,但主要集中于氣泡通過自身浮力豎直上升方面,對水平通道內受液相強制流動下驅動的氣泡運動行為關注較少,而此類兩相流動問題在工程應用中十分常見(如潛艇水下排氣系統水平通入廢氣進行降溫降噪)．本文使用Liang等[19]提出的基于Allen-Cahn方程的相場LB模型,對受液體推動的氣泡在含孔板結構的水平微通道內的運動形態進行了模擬研究;通過考察氣泡變形、分裂、合并、質量損失和速度變化等動力學行為,分析了不同氣泡狀態參數和表面潤濕性對氣泡穿過孔板過程及其后續演化的影響．

1 數 值 方 法

在氣液兩相流動的數值計算中,兩相界面的捕捉計算是一個復雜的問題,為準確描述整個流體系統中多種非混溶流體共存的狀態以及它們在相互作用下的運動,相場模型將任意位置上非混溶流體所占的體積分數定義為該點的序參數,氣液相界面由一個連續的薄層表示．

Liang等[19]提出的相場LB模型采用兩個分布函數分別用于捕獲相界面和描述速度/壓力場,其中序參數分布函數滿足保守Allen-Cahn方程,如下所示：

(1)

式中,φ為序參數,φ取0和1分別表示氣相和液相,當0<φ<1時為氣液間界面;M為遷移系數;t表示時間;n=?φ/|?φ|,n是垂直于界面單位向量;λ為序參數φ的函數,

(2)

式中,W為界面厚度．

流體運動速度u和壓力p等宏觀量還需滿足不可壓縮的Navier-Stokes方程．對于Allen-Cahn方程和Navier-Stokes方程,具有BGK碰撞算子的LB演化過程可分別表示為

(3)

(4)

在兩相流系統中,由于存在不同物質,流體動力黏度μ并非定值,動力黏度由序參數和氣液兩相各自的動力黏度線性差值得到,計算公式如下：

μ=φ(μl-μg)+μg．

(5)

(6)

(7)

(8)

(9)

在式(3)、(4)中,外力項Gi(x,t)和Fi(x,t)的計算公式如下：

(10)

(11)

其中偏導項?t(φu)可由顯式Euler法計算：

(12)

為描述氣泡與黏性液體的相互作用,力項F由表面張力Fs與外力Fg相加而得,Fs=μφ?φ,μφ=4βφ(φ-1)(φ-0.5)-κ?2φ是化學勢．式中,κ和β是常數,由界面厚度W和表面張力σ決定,即κ=1.5σW,β=12σ/W．

通過對分布函數求和,可以得到序參數、壓力和速度等宏觀統計量：

(13)

(14)

(15)

數值迭代中,對模型內導數項進行差分離散,梯度項和Laplace算子采用二階各向同性中心格式計算：

(16)

(17)

其中,離散速度ei的表達式為

(18)

式中,c為格子速度,計算方式為c=δx/δt,δx和δt分別表示格子長度和格子時間,本文中取δx=δt=1．

而對流固界面處的潤濕性參數,本文采取虛擬密度法進行構建[20],以微通道底部邊界為例,對于固體邊界處的格點,虛擬節點上的序參數計算式如下：

(19)

其中導數項采用二階中心差分法計算：

(20)

2 物 理 模 型

本文所研究水平管道的幾何模型如圖1所示,在計算域為Nx×Ny=840×120的微通道內,包含一塊垂直放置的孔板結構,孔板的左表面距離計算區域左邊界Lqx=360,孔板厚度wd=15,孔板的開孔半徑rd=6,板上四個開孔的中心位置縱向坐標分別為15,45,75,105．

圖1 模型示意圖

初始時刻,在孔板前方放置一個半徑為r,動力黏度μg=1的圓形氣泡,其圓心距計算區域左沿距離cx=60,距計算區域上沿距離cy=60,周圍區域為動力黏度μl=100的液相流體．氣相密度ρg設置為1.0,液相密度ρl為1 000,其他參數固定為界面厚度W=4.0,遷移率M=0.1．氣泡在水平管道流動的體積力推動下(體積力Fg=gx×ρ)經過孔板并流出計算區域．本文并未考慮豎直方向的浮力．對微通道的上下壁面采用半步反彈邊界,左右邊界采用周期邊界．如無特殊說明,文中所涉及的仿真參數和計算結果均為格子單位．

3 模 型 驗 證

3．1 Laplace驗證

與其他研究一樣,本文首先使用Laplace定律驗證了模型對于氣液分離的正確性．Laplace定律表明,對處于穩定狀態的氣泡,氣泡的內外壓差ΔP與氣泡半徑r的倒數呈線性關系,即滿足以下關系式：

ΔP=|Pin-Pout|=σ/r,

(21)

式中,Pin為氣泡內部流場的平均壓力,Pout為氣泡外部流場的平均壓力．

設置計算區域Nx×Ny=200×200,四周均采用周期邊界．液相和氣相密度分別為ρl=1 000和ρg=1．在初始時刻,計算區域被液相填充,一個圓形氣泡被置于區域中心,通過改變氣泡半徑大小得到不同壓差．半徑分別選取r=15,20,25,30,35,并分別在三種不同表面張力σ=0.1,0.2,0.3下進行模擬．如圖2所示,數值結果與式(21)吻合較好,ΔP與1/r呈線性關系,符合Laplace定律．

圖2 Laplace驗證

3．2 潤濕性驗證

為了驗證當前潤濕邊界方案在實現非穩態運動過程中接觸角方面的能力,本小節以毛細管的侵入為典型潤濕現象進行模擬．在尺寸為L×H=400×20的毛細管中,初始時刻由未潤濕性氣體填充,然后外部潤濕性液體侵入管內．總計算區域設為Lx×Ly=800×36,其中毛細管占據區域為200≤x≤600,氣體組分的位置為240≤x≤750．液氣密度比為ρl/ρg=1 000,黏度比設置為μl/μg=100(其中μl=100,μg=1)．其他參數取σ=0.2,M=0.1,W=4．對所有管壁施加無滑移邊界條件,對毛細管管壁邊界采用潤濕性邊界方案,其余邊界設置為周期邊界．如果忽略重力和慣性效應,相界面的演化可以描述為

(22)

如圖3所示,本小節計算了在三種不同的規定接觸角θ=30°, 45°和60°下,使用潤濕邊界方案模擬得到的毛細管侵入距離和公式計算結果的比較．需要注意,由于式(22)中未考慮慣性效應、初始條件以及進氣道效應等因素,本模型數值結果會小于式(22)的計算結果．且隨著接觸角的增加,由σcosθ決定的毛細管力減小,在較大接觸角下,二者之間的偏差還會增大．圖3中還將本文模擬的毛細管的侵入距離與采用類似潤濕性方法的模擬結果[20]進行對比,兩者最終結果在三個不同接觸角下的偏差分別為1.5%,0.9%,0.4%,故可認為該潤濕性方法是有效的．

(a) θ=30°

(b) θ=45°(c) θ=60°圖3 潤濕邊界對毛細管侵入的模擬圖

3．3 網格獨立性驗證

圖4給出了在四個不同密度的網格下,氣泡在進入及離開孔板時,其界面輪廓的局部圖．可以看出Nx×Ny=840×120和1 120×160網格得到的氣泡輪廓幾乎相同．而由圖5顯示,氣泡通過孔板時,隨著網格密度的增大,氣泡的氣體平均速度vg和通過孔板時間t*的相對偏差ε逐漸減小,并且在Nx×Ny=840×120和1 120×160兩種網格密度下與密集網格的相對偏差ε均小于1%,這樣的偏差是可以接受的．綜合考慮計算精確度需要和運算成本,本文采用840×120網格進行計算．

圖4 氣泡進入孔板和離開孔板時,不同網格下的氣泡輪廓

圖5 網格獨立性驗證

4 計算結果與分析

4．1 We數對氣泡通過孔板的影響

氣泡在液體中的變形和破碎主要受外力和表面張力之比的影響,這一關系可以用We數進行描述,本小節主要討論不同We數下氣泡在水平管道的運動及撞擊孔板表面的動態過程．

圖6為不同We數(We=0.588,1.176,4.704)下氣泡通過孔板的動態過程,此時Re=4.38,氣泡半徑r為32個格子單位．可以看出, 隨著We數的增加, 氣泡在通過孔板時呈現出不同的形態．如圖6(a)所示,We=0.588時,氣泡在孔板結構前受到連續流體推動向前運動,在表面張力和黏性力作用下逐漸達到穩定,呈現“子彈形”．當氣泡逐漸向右運動并接近孔板時,氣泡受中間障礙產生的兩個前緣呈現曲率半徑呈減小趨勢,當前緣到達孔板位置時達到最小值．然后,隨著前緣穿過孔喉,其曲率半徑回升增大,這一趨勢與實驗中氣泡穿過孔隙結構的現象相吻合[21]．由于孔板部分管道寬度快速變窄,氣泡平均速度劇烈上升,留在孔板前端的“尾部”氣膜變薄,氣泡受到強烈的黏性切斷力[7]產生變形和破裂．由于氣泡內部存在對稱的黏性應力,氣泡破裂為兩個大小幾乎相同的子氣泡．受孔板后方回流影響,破裂子氣泡在通過孔板后速度隨即下降,快速互相靠近發生二次聚合．聚合形成的新氣泡在表面張力作用下,排出內部小液滴逐漸恢復穩定,再度發展至子彈形．如圖6(b)所示,在We=1.176工況下,氣泡表面張力減小,故在運動過程中形成子彈形后,會繼續發展為“彈弓形”,處在管道中間氣膜的厚度減小而兩側體積增大．因此,當氣泡在通過孔板結構時,部分氣體會從靠近管道壁面的開孔通過,整體被分割為四個主要子氣泡．中間兩個子氣泡如同剛才一樣重新匯聚后發展為稍小的子彈形氣泡,由于破碎后結構的不對稱性,很多情況下該破碎氣泡會被流動推向管道一側(t*> 3.5時)．靠近兩側的子氣泡由于速度差距較大而獨立發展,形成扁平狀氣泡并在靠近壁面的位置運動．繼續增大We數,如圖6(c)所示,當We=4.704時,由于表面張力進一步縮小,氣泡在運動到孔板之前已經產生巨大形變,直接被液體流動完全分解(t*=1.991),大部分殘余氣體被擠壓到兩側的低速流動區,在通過孔板時被進一步撕裂,在管道兩側形成連續的扁平泡狀流．而“頭部”的少量氣體在通過孔板后重新匯集成箭頭形向前流動,運動速度遠高于兩側的分裂氣泡．

(a) We=0.588 (b) We=1.176 (c) We=4.704圖6 d/D=0.533 33的氣泡通過孔板的動態行為

從圖7的氣體與液體流速隨時間變化圖可以觀察到,氣泡在低We數下通過孔板是一個連續過程,氣體與液體速度存在一個共同的波峰．隨著We數的增加,氣泡變形程度提高,部分被分割的子氣泡在孔板靠外的位置通過,從而使氣體平均速度產生了兩個存在時間錯位的峰值．由于中部匯聚氣泡在一段時間后被推離管道中心,氣體平均速度在后續時間(t*> 3.5時)相比通過孔板前會更低．而We數更高時,氣泡經過前孔板已經被撕裂,但位于管道中部的“頭部”殘余子氣泡也因此有了更高的速度,氣泡的“頭部”部分更早運動到孔板位置后流過,使氣體與液體流動速度均有提前的小幅加速．而兩側氣泡偏向長條的形狀,使孔板結構持續會有小型子氣泡分解流過,故氣液平均速度提高的時間會持續得比之前兩種情況更久．

(a) 液體速度 (b) 氣體速度(a) The liquid velocity (b) The gas velocity圖7 氣液速度隨無量綱時間t*的變化

4．2 直徑比對氣泡融合和脫落過程的影響

本小節討論在Re=4.38條件下,氣泡大小對其通過橫向管道孔板的影響．為方便參考,選擇氣泡直徑與管道直徑之比d/D為氣泡大小的參照．在模擬中分別設置了d=0.266 7D至d=0.533 3D等多種不同大小的氣泡．

圖8給出了在不同直徑比d/D和We數下,氣泡流經孔板的相圖．不同氣泡在通過孔板時,會表現出三種形態：先分裂再聚合、直接分裂和先撕裂再分裂．在直徑比d/D和We數均較小時,氣泡主體會在被孔板分裂為兩部分后再次聚合為一個氣泡運動．當氣泡處于較高We數時,其直徑比d/D會出現兩個臨界值．當d/D大于第一個臨界值而小于第二個臨界值時,氣泡會被孔板分裂為三個主要子氣泡運動,而當d/D超過了第二個臨界值時,氣泡將在孔板前直接被流動撕裂,并被孔板進一步分裂為眾多小氣泡．這兩個臨界直徑比的大小都會隨著We數的增加而降低．

圖8 不同1/We和d/D下氣泡通過孔板的動態行為

4．3 潤濕性影響

孔板孔喉部分毛細壓力受表面潤濕性的影響,會改變非潤濕相通過孔喉的能力．本小節為研究孔板表面潤濕性對水平微通道內氣泡-水流動特性的影響,分別對不同接觸角情況進行了一系列數值模擬,包括θ=30°,60°,90°,120°,150°．氣泡參數如下：d/D=0.4,Re=4.28,W=4.0．

其中三種不同接觸角下氣泡的運動過程如圖10所示,孔板表面潤濕性變化對氣泡在橫向微通道內的動態行為存在較大影響．如圖10(a)所示,當θ=60°時,孔板表面處在疏氣親水工況,在橫向驅動力、表面張力和孔板表面阻力等共同作用下,氣泡前端進入孔板之中并與孔板表面發生接觸．之后,氣泡以極快速度通過孔板,在加速和減速過程中受到孔板表面阻力和表面張力的共同作用,有小部分氣體被撕裂獨立發展(t*=2.756 3),主體氣泡重新恢復為子彈狀(t*=3.368 8)．

圖9 氣泡運動過程中在不同d/D下幾種運動參數隨We數的變化

(a) θ=60°

(b) θ=120°

(c) θ=150°圖10 氣泡通過不同潤濕性孔板動態行為

當接觸角θ=120°時,如圖10(b)所示,相比于θ=60°工況,由于孔板障礙表面親水性減弱,導致氣泡在通過孔板過程中與表面的接觸面積增大．在通過孔板結構的速度變化過程中,氣泡同樣出現了破裂現象,但在親氣表面產生的破碎小氣泡會傾向于黏在孔板后方(t*=2.756 3)．主氣泡與親水工況一樣演化為子彈形狀(t*=3.368 8)．如圖10(c)所示,隨著孔板表面親氣性進一步加強,接觸角增大為θ=150°,氣泡幾乎緊貼著孔板表面穿過孔喉(t*=2.327 5),受到結構表面毛細力影響,氣泡會整體被孔板結構捕獲,被捕獲后的氣泡滯留在孔板結構后方(t*=2.753 6)．氣泡前端會產生兩個先行分離的小氣泡,氣泡的主體隨后在液流驅動下離開孔板,繼續發展演化．

(a) εg

圖11 氣泡運動參數隨孔板潤濕性變化

為定量分析孔板潤濕性改變對氣泡運動的影響,本文統計了不同接觸角下不同We數氣泡的運動數據和殘余質量εg．氣泡殘余率如圖11(a)所示,當孔板表面親氣性較弱時,氣泡不會被結構吸附,幾乎可以整體穿過孔板(θ=30°,60°)．當θ>60°后,隨著表面親氣性的加強,氣泡殘余質量逐漸減小,但整體質量損失不大(均小于10%)．觀察圖11(b)不難發現,當θ=60°,90°,120°,150°時,氣泡的最大平均速度比分別為3.92,3.93,4.04,4.65(不同We數的氣泡取平均),即當孔板表面開始成為親氣表面時,接觸角以相同幅度增大,但氣泡運動中的最大平均速度分別提高了0.3%,2.8%和15.2%．因此隨著接觸角的增大,通過孔板結構的平均氣泡速度也逐漸增大,此外當接觸角大于120°后,平均最大氣泡速度隨著接觸角增大的程度更加劇烈,特別對We數較低的氣泡尤為明顯．圖11(c)顯示,在接觸角增大的過程中,周圍液體的最大平均速度也出現了和氣體速度類似的現象,不過在大接觸角下的增加程度沒有氣體速度顯著．

5 結 語

本文基于Allen-Cahn方程的大密度比相場LB模型,對氣泡在含孔板微通道內的運動過程進行了數值模擬研究．分別考慮了We數、氣泡當量直徑和孔板表面潤濕性等因素對氣泡平均速度和破碎性等運動學特性的影響,得到了如下結論：

1) 隨著We數的增大,氣泡在孔板前的形變會更加劇烈,表現為氣泡邊緣到兩側管壁的距離逐漸縮短,氣泡在通過孔板時更容易發生破裂、匯聚等形變;氣泡通過孔板過程中,當氣泡大小和流動強度均相同時,We數的增大會使速度變化的時間更長．

2) 在本文所研究的參數范圍內存在兩個臨界直徑比,這兩個臨界直徑比將氣泡通過孔板的運動分成了“分裂后聚合為一個氣泡”“被孔板分裂”和“在孔板前被撕裂”三種形態．而氣泡直徑比d/D的增加會使得氣泡運動過程中出現最高速度時所需要的We數逐漸變小,最高氣泡速度、單點速度和液體速度峰值則會隨之增加．

3) 隨著接觸角的增大,孔板結構親氣性越強,孔板表面對氣泡的吸附作用越明顯．此時氣泡與孔板的接觸面積增大,被吸附的氣體也越多,導致穿過孔板的氣泡體積逐漸減少．且隨著接觸角的增大,氣泡的最大平均速度增加,氣泡通過孔板的時間減少,這在接觸角大于120°后尤為明顯．

實際微通道中的氣泡運動是一個三維問題,本文采取二維模型計算,后續應通過優化并行效率計算三維情況,使結果更符合真實情況,尤其是對于計算含浮力水平通道內流動的問題[22]．我們希望這一數值研究有助于理解氣泡動力學的基本物理原理,以及在未來的微結構流體應用中設計更為接近現實的復雜流動．