大密度比雙氣泡在孔板結構微通道內(nèi)上升行為的格子Boltzmann 方法模擬

嚴 裕，婁 欽，*

(1. 上海理工大學 能源與動力工程學院，上海 200093；2. 上海市動力工程多相流動與傳熱重點實驗室，上海 200093)

0 引 言

氣泡在復雜微通道內(nèi)的流動現(xiàn)象廣泛存在于自然界和工業(yè)生產(chǎn)過程，如沸騰現(xiàn)象、頁巖氣開采、燃料電池的使用等。例如在燃料電池中氣泡間的碰撞和合并時常發(fā)生，氣泡在液體中的運動和相互作用直接影響相間的接觸時間和接觸面積，并改變氣液間的傳熱傳質(zhì)[1]，進而影響其使用性能[2-4]。因此對氣泡進行微通道內(nèi)的流動研究以及工業(yè)生產(chǎn)過程中的優(yōu)化具有十分重要的意義。

學者們對復雜微通道內(nèi)的氣泡流動問題進行了大量的實驗研究。Chen 等[5]研究了高黏度流體中氣泡在孔板通道內(nèi)上升時的形狀變化，發(fā)現(xiàn)氣泡穿過孔板通道后變?yōu)椤霸卵罓睢?。Corapcioglu 等[6]提出采用定量方法預測氣泡在粗糙多孔介質(zhì)中的上升速度，其預測結果表明氣泡的上升速度經(jīng)過很短的時間后會達到平衡。Dawson 等[7]研究了在恒定體積流量中孔板微通道內(nèi)氣泡的流動過程，發(fā)現(xiàn)隨著流速的增加，氣泡的后尖端速度發(fā)散且氣泡尾部變薄向中間收縮回主體。

以上實驗結果得到了氣泡運動過程中的直觀現(xiàn)象和宏觀規(guī)律，為人們認識氣泡在通道內(nèi)的運動規(guī)律提供了理論基礎。然而由于實驗方法難以捕獲微孔內(nèi)的流動細節(jié)，采用數(shù)值模擬來研究兩相流動問題也逐漸成為人們研究氣泡運動的一種常用方法[8-9]。Sussman 等[10]應用Level-set 方法，研究了單個氣泡上升過程中的界面形變以及破裂現(xiàn)象，發(fā)現(xiàn)隨著表面張力和黏度比的增加，氣泡越不容易發(fā)生變形和破裂現(xiàn)象；而隨著雷諾數(shù)的增加氣泡的變形會越來越劇烈。Sun 等[11]利用格子Boltzmann 方法模擬了單個氣泡在含三個半圓形障礙物的通道內(nèi)上升的問題，發(fā)現(xiàn)了障礙物的存在可以影響氣泡的運動形態(tài)，使其運動軌跡發(fā)生偏轉。Lou 等[12]通過兩相流格子Boltzmann模型，模擬研究了氣泡在浮力作用下上升并通過復雜微通道的問題。研究發(fā)現(xiàn)障礙物表面潤濕性越強，氣泡變形越劇烈、氣泡終端速度越?。磺译S著障礙物半徑的增大，氣泡的終端速度呈線性減小。

以上研究揭示了一些氣泡在復雜微通道內(nèi)的運動機理，但主要集中于單氣泡的研究，較少涉及多氣泡在復雜微通道內(nèi)的運動行為，且大部分研究只能考慮較小的氣液密度比。因此，本文選用Liang 等[13]基于Allen-Cahn 相場理論提出的大密度比格子Boltzmann 模型，對兩個上下排列的氣泡在含孔板結構的微通道內(nèi)的運動形態(tài)進行研究；主要調(diào)查了Eo數(shù)、氣泡之間的相對大小、兩氣泡之間的距離以及氣泡和孔板之間的距離等因素對兩個氣泡變形、合并、分裂的動力學行為以及氣泡穿過孔板的質(zhì)量和通過時間的影響。

1 數(shù)值方法

2 物理問題

圖1 D2Q9 速度模型示意圖Fig. 1 Discrete velocities for D2Q9 model

圖2 模型示意圖Fig. 2 Model sketch

3 模型驗證

3.1 拉普拉斯驗證

在本文驗證中，初始氣泡半徑選20、25、30、35、40 這五種情況，并考慮不同表面張力σ= 0.2、0.4、0.6、0.8、1.0，模擬結果如圖3 顯示。氣泡半徑的倒數(shù)與表面張力成正比，模擬結果與Laplace 定律吻合，證明了本文所使用模型的準確性。

圖3 Laplace 驗證Fig. 3 Numerical validation

3.2 氣泡在光滑通道內(nèi)上升行為

氣泡在光滑壁面通道內(nèi)受浮力上升的問題與本文所要研究的問題相似, 且已有學者對該問題進行了大量研究, 因此本小節(jié)將采用該算例驗證模型的正確性。模擬計算區(qū)域為Lx×Ly=80×300，液相密度ρl=1 000.0， 氣泡密度 ρg=1.0， 氣泡半徑r=20，界面厚度W=3.0 ，初始時刻氣泡圓心設在（Lx/2，40）處。邊界條件設置：左右為無滑移固壁，上下為周期邊界。

考慮Eo=5、10、20、40 四種工況，Re數(shù)固定為17.7。表1 列出了不同工況下，本文所使用的LBM方法和其他數(shù)值方法得到的氣泡穩(wěn)定形狀對比，可以看出隨著Eo數(shù)的增大，氣泡從球形變?yōu)樯习嫉那蛎毙?。本文得到的結果與Takada 等[16]使用雙流體LBM 模型得到的模擬結果以及Ren 等[17]使用相場LBM 法模擬小密度比氣泡得到的模擬結果高度相似，且在大密度比工況下得到了數(shù)值穩(wěn)定的結果。之后統(tǒng)計氣泡終端速度，如表2 所示，并與Alizadeh 等[18]和VOF 方法[19]進行對比。本文結果與兩者結果一致，故本模型適用于研究氣泡動態(tài)行為。

表1 不同 Eo數(shù)下氣泡穩(wěn)定形態(tài)Table 1 Bubble shapes at steady states under different E o numbers

表2 不同 E o 數(shù)下氣泡終端速度Table 2 Terminal velocities of bubble under different E o numbers.

4 數(shù)值結果與分析

本文中主要研究E?tv?s數(shù)（Eo）、氣泡之間的相對大小和氣泡的初始相對位置對氣泡在孔板微通道內(nèi)流動過程的影響。

4.1 E o數(shù)的影響

的中間以外，上方氣泡與下方氣泡在兩側也開始合并，導致上下方氣泡之間包裹兩個極小液泡（t?=7.10）；隨后在表面張力作用下液泡從兩側被排出，此時通道內(nèi)的兩個氣泡在孔板附近完全合并形成一個完整的大氣泡（t?=9.98），大氣泡通過孔板通道后在孔板正下方殘留兩個小氣泡，而其余部分通過孔板后在孔板上方近似恢復為橢球形（t?=18.13）。

圖4（b）所示為Eo=10的工況，此時氣泡表面張力相對于浮力減小，因此兩個氣泡在上升的過程中都有較大的形變。上方氣泡下表面向上凹陷，在兩端形成尾翼，下方氣泡水平拉長且上表面向上凸起與上方氣泡下表面近似貼合；t?=7.09時刻下方氣泡的頂點與上方氣泡中間接觸，接著上方氣泡兩端向下彎曲，下方氣泡上表面兩側與其接觸，氣泡間形成兩個液泡，因為表面張力相對浮力減小，所以此時形成的液泡比Eo=5的工況下生成的大。由于孔口的阻力小于兩側，下方氣泡質(zhì)量主要從中間的接觸點并入上氣泡，t?=9.19時刻上氣泡與兩端斷裂，其兩端并入下氣泡兩端，液泡也被氣泡從斷口排出。隨后在t?=11.04時刻兩個氣泡合成了一個主氣泡（能夠排出的氣泡），在部分質(zhì)量穿過孔口后，主氣泡運動速度增大與孔板正下方部分質(zhì)量再次斷裂，這些質(zhì)量最終殘留在孔板下方。主氣泡在穿過孔口后由于拉伸再次破裂出一個小氣泡，該小氣泡尾隨主氣泡上升（t?=20.76）。對比圖4（a）和圖4（b）可以發(fā)現(xiàn)，隨著Eo數(shù)增加，氣泡在穿越孔板前后的形變更明顯。

圖4（c）為Eo=20得到的氣泡運動過程。從圖中也發(fā)現(xiàn)隨著Eo數(shù)的持續(xù)增大導致氣泡變形和斷裂更嚴重。例如在氣泡的上升過程中，上方氣泡尾端形成了更長的尾翼，下方氣泡也更扁長。而更多的接觸面積和更小的表面張力也讓上下兩個氣泡間形成更多的接觸點（t?=7.36），兩個氣泡中間存在著更多更大的液泡，盡管此時氣泡在過孔板通道前依舊可以將液泡擠出，但也導致主氣泡和其兩端的連接容易斷裂（t?=11.30），殘余在孔板下的氣泡質(zhì)量增加。主氣泡在通過孔板后也分裂出一個更大的尾氣泡（t?=20.76）。而當Eo數(shù)增加到40 時（圖4（d）），上下氣泡發(fā)生了多次接觸、合并與破裂的過程。在主氣泡通過孔板通道后，氣泡內(nèi)部還夾帶著一個小液泡與其上升了一段距離才被排出（t?=21.55）。從圖中可以看出四種工況下氣泡首次接觸的時間分別為t?=6.71、7.09、7.36、7.61，接觸時上氣泡都已進入孔道，主氣泡到達計算區(qū)域上方的時間如圖5 所示。因此隨著Eo數(shù)的增大，氣泡接觸時間越晚，氣泡間接觸面積增大，接觸點隨之增多；接觸點的增多并不能促進兩氣泡的合并，反而導致中間生成更多液泡，加劇氣泡分裂，氣泡通過孔板通道需要的時間也越長。氣泡到達頂部的所需時間隨著Eo數(shù) 的增大而增多，并且當Eo數(shù)大于10 后呈線性上升的趨勢。

圖4 不同 Eo數(shù)下氣泡運動過程Fig. 4 The temporal evolution of bubbles under different E o numbers

上述四種工況下，Eo數(shù)對兩個氣泡通過孔板通道的剩余質(zhì)量比Ms有影響，為了分析這種影響，本文對不同Eo數(shù)下的氣泡剩余質(zhì)量進行統(tǒng)計。由圖5 對Ms曲線的擬合發(fā)現(xiàn)，氣泡剩余質(zhì)量的變化趨勢隨Eo數(shù)的增加呈負指數(shù)級地減小；減小是因為隨著Eo數(shù)的增大，氣泡表面張力相對于浮力減小，氣泡難以維持初始形狀，其在上升和變形的過程中更容易包裹液泡，并更容易發(fā)生破裂，導致殘留在孔板下方。而隨著Eo數(shù)的持續(xù)增大，氣泡的最終剩余質(zhì)量趨于穩(wěn)定；這是因為盡管氣泡被拉伸，但氣泡初始位置與孔板間的距離導致氣泡還未充分形變便到達了孔板位置，有一半多的質(zhì)量處在孔口正下方且最終通過孔板。

圖5 不同 Eo數(shù)下剩余質(zhì)量比和氣泡到達頂部時間Fig. 5 The variation of residual mass and the time for bubble to reach the top for different E o number

4.2 氣泡相對大小的影響

本小節(jié)研究了氣泡相對大小對氣泡通過孔板微通道的影響。數(shù)值模擬中參數(shù)設置為 σ =0.624；為了方便比較，選擇孔口間距w作為參照尺寸，上下兩氣氣 泡 大 小d1:d2分 別 設 為 2.4w:1.5w、 2.0w:2.0w、1.5w:2.4w；三種氣泡方案的質(zhì)量和近似；兩氣泡初始位置分別設為（x1c，y1c） =（Lx/2 ， 100）和（x2c，y2c）=（Lx/2，40）。

圖6（a）展示了d1:d2=2.4w:1.5w，即上方氣泡大于下方氣泡時氣泡穿過孔板的運動過程。如圖所示，在t?=9.99時下方氣泡接觸到上方氣泡，此時上方氣泡部分已經(jīng)穿過孔板通道，下方氣泡與上方氣泡尾翼接觸并將一個直徑約為1 .0w寬的液泡包裹在內(nèi)。液泡隨氣泡的上升與其一起穿過孔道。當液泡處于孔板通道內(nèi)時氣泡被夾斷分成上下兩個氣泡（t?=11.56），之后上氣泡雖然已經(jīng)穿過孔板通道，但由于形變較嚴重，維持了較低的運動速度。當液泡受到下氣泡對其的推動力后從孔口被排出，下氣泡加速向上運動與上氣泡合并（t?=17.35），之后尾部再次破裂產(chǎn)生一個小氣泡，并在表面張力的作用下，逐漸恢復為球帽狀。

圖6（b）給出了上下兩氣泡直徑相等時（d1:d2=2.0w:2.0w）氣泡的上升過程。如圖所示，上氣泡靠近孔板后拉伸為長條狀，而下氣泡也被拉長，上氣泡下方和下氣泡上方貼合，在中間和兩端產(chǎn)生接觸點（t?=7.36），氣泡內(nèi)部形成兩個液泡。上氣泡與其兩端斷開，接著液泡從斷口排出（t?=11.04），上氣泡端部匯入下氣泡的端部。但此時下方氣泡質(zhì)量主要在中部并且已并入上方氣泡，其與兩端的連接也迅速斷開。主氣泡加速通過孔板通道并繼續(xù)上升（t?=15.51），之后在上升過程中氣泡尾部又斷裂出一個小氣泡。

上方氣泡直徑小于下方氣泡（d1:d2=1.5w:2.4w）的結果如圖6（c）所示。在t?=6.05時上方氣泡靠近孔板后變扁平狀中間凸起。下方氣泡從中間和兩側同時接觸到上方氣泡，形成兩個小液泡。與氣泡直徑相等的情況不同的是，下氣泡主要從兩端匯入上氣泡，再從兩端向孔板通道中間收縮（t?=11.04），過程中兩個小液泡被擠壓合并為一個直徑為w的液泡。液泡在孔板通道下方主要受到氣泡對其的表面力，方向從向上轉為向中間，液泡幾乎處于靜止狀態(tài)，未跟隨氣泡上升，隨后液泡被主氣泡從下表面排出。主氣泡穿過孔板通道繼續(xù)上升。

圖6 不同氣泡直徑下氣泡運動過程Fig. 6 The temporal evolution of bubbles with different diameters

為了研究氣泡相對大小對氣泡通過孔板通道后剩余氣泡質(zhì)量的影響，增加了幾組氣泡直徑比不同的工況以便比較。圖7 給出了三組相同變化比例但不同初始質(zhì)量和的情況，統(tǒng)計了氣泡穿過孔板后的剩余質(zhì)量以及氣泡運動到計算區(qū)域上方的時刻。從圖7中間組可以看出，當上氣泡質(zhì)量大于下氣泡質(zhì)量時，上氣泡位于孔板正下方部分質(zhì)量接觸到下氣泡后在表面張力作用下向中間收縮，最終有較多的剩余質(zhì)量可以達到頂部。而當下方氣泡大于上方氣泡時，兩氣泡發(fā)生合并的時間和位置提前，盡管下方氣泡是從兩側匯入上氣泡，但兩個氣泡形變程度較小，合并的氣泡質(zhì)量向孔口方向移動，故最終穿過孔板通道的質(zhì)量比前者情況少。從中間組可以看出上下氣泡直徑對比度越大，剩余質(zhì)量越多；而當上下兩個氣泡大小相等，穿過孔板到達通道上半部分的剩余質(zhì)量最少，為0.82。在相反的直徑比例下，質(zhì)量大的氣泡處在上方的情況最終的剩余質(zhì)量始終比處在下方的情況多。另外，隨著下方氣泡相對上方氣泡直徑增大，氣泡通過孔板到達頂部所需時間變短。結合圖7 上下兩組可以看出氣泡尺寸相對孔徑越小，其通過得越完整，即剩余質(zhì)量越多；而較小的質(zhì)量也導致其受到的浮力相對較小，到達頂部時間變慢，在不同直徑比下的剩余質(zhì)量和到達時間也與中間組的變化趨勢相似。

圖7 不同氣泡直徑下的剩余質(zhì)量比和達到時間Fig. 7 Residual mass and the time for bubbles to reach the top for different diameters

4.3 兩氣泡之間的距離與氣泡與孔板的距離的影響

本小節(jié)研究兩氣泡之間的距離以及氣泡與孔板的距離對雙氣泡運動行為的影響。首先研究不同氣泡之間距離的情況，為了方便比較，此算例也以孔板間距w作為參照尺寸，l1和l2都 表示為w的倍數(shù)。在數(shù)值模擬中設置了三組工況，首先使上氣泡的位置固定不變，此時上方氣泡圓心與孔板間下表面距離l1=3.0w為固定值，接著改變下氣泡的位置，使上下氣泡圓心之間的距離l2分 別為2.5w、 3.0w和 3.5w三種情況。此外，Eo=10，d1=d2=20。

當兩個氣泡圓心距離l2=2.5w時，所得到的結果如圖8（a）所示。上氣泡下表面和下氣泡上表面的最近點距離為 0.5w，t?=3.68時下氣泡向上浮動從中間接觸并入上氣泡，并在上氣泡到達孔板之前完成合并（t?=6.04），氣泡間未有液泡形成。合并完成后的主氣泡橫向拉抻，氣泡中間部分質(zhì)量在浮力作用下穿過孔口；在表面張力作用下氣泡兩側向中間收縮穿越孔口，最快地穿過了孔板通道（t?=12.35），穿過后恢復為球帽形并保持勻速上升運動。

當兩個氣泡圓心距離l2=3.0w時，上氣泡下表面和下氣泡上表面的最近距離為1 .0w，此時下氣泡并入上方氣泡的同時上方氣泡準備穿過孔板通道；當t?=8.41時下氣泡先后從中間和兩端接觸到了上方氣泡，中間夾著液泡。液泡把上氣泡與其兩端擠斷，接著液泡從斷口排出（t?=10.77），上方氣泡兩端與下方氣泡的兩側合并；但此時下方氣泡多數(shù)質(zhì)量集中在中間并已并入上方氣泡。下方氣泡與兩側的連接斷開，繼續(xù)上升。

圖8（c）給出了兩個氣泡圓心距離l2=3.5w時的結果。從圖中可以看出，t?=9.20時刻在兩個氣泡合并前，上氣泡中間部分已經(jīng)與兩側脫離后先行穿過孔板通道，穿過孔板通道時速度上升，但其所占整體質(zhì)量的比重小，平均上升速度提升不明顯；兩側殘余的質(zhì)量與下氣泡合并，t?=12.88時上氣泡離開孔板通道的同時下方氣泡進入孔板通道，氣泡之間存在液泡，上下兩氣泡運動速度均下降；當下方氣泡頭部穿過孔板通道并把液泡從孔板通道擠出后，其速度上升并加速和上氣泡合并成一個氣泡。

圖8 氣泡不同初始位置情況下氣泡運動過程Fig. 8 The temporal evolution of bubbles under different initial positions

經(jīng)模擬計算發(fā)現(xiàn)，當l1=3.0w、l2=3.3w時，幾乎在上方氣泡與兩側脫離的同時下方氣泡與其兩側接觸，而上氣泡主體已開始加速上升。相比l2=3.25w的情況，兩氣泡主體接觸的位置從孔口下方約 0.5w處躍遷到了距離孔口上方1 .5w的位置；氣泡主體第一次接觸的 時 間 也 從t?≈6.0推 遲 到 了t?≈14.2。 在l1從 2.5w增加至 3.5w的 過程中，臨界的l2從 3 .1w變 化為 3 .45w。

圖9 展示了兩個氣泡從不同位置開始上升時統(tǒng)計得到的速度變化曲線。從圖9 可以看出，在氣泡未到達障礙物通道前，氣泡的上升速度變化趨勢基本相同；但當上氣泡與孔板的距離l1不 變時，在l2=2.5w氣泡間距較小的情況下，初始運動的速度提升最大。這是因為兩氣泡在通過孔板通道前完成了合并，合并后的氣泡受到的浮力增大；在完成合并后通過孔板通道時速度變化劇烈，作為整體通過孔道時減速更明顯，通過孔道后的加速也越劇烈；之后因表面張力作用與孔板下方質(zhì)量發(fā)生斷裂時, 速度略微下降, 最后再度趨于平穩(wěn)。而當兩氣泡間距增大（如l2=3.0w），兩氣泡在孔板下方發(fā)生接觸但未完全合并，這也使其速度變化與間距較小時相同，但變化幅度更平和。如果兩氣泡間距進一步增大，則上氣泡部分質(zhì)量先進入孔道，另外部分質(zhì)量并入下氣泡，穿越孔道的過程分成兩次，加速過程隨之變?yōu)閮纱?；其中第二次在t?=17.35時，出現(xiàn)轉折是在下氣泡通過孔道后接觸到上氣泡并匯入的時刻。

圖9 不同初始位置下氣泡上升的平均速度 v?與 時間 t?的關系Fig. 9 Average bubble velocity versus time for different initial positions

接下來研究上方氣泡與孔板間的距離對雙氣泡運動行為的影響，這里固定兩氣泡間距，通過改變上氣泡的位置來改變氣泡與孔板的間距l(xiāng)1。圖9 的實線和虛線分別表示l1為 3.0w和 3.5w的情況。從圖中虛線和實線的對比可以發(fā)現(xiàn)，氣泡離孔板越遠（l1越大），氣泡的運動相對越遲滯；初始時表現(xiàn)為速度上升得更高，之后在穿越孔道前后時的速度變化規(guī)律相同；氣泡穿過孔板通道后又恢復到勻速運動。

為了分析氣泡初始位置和初始間距對氣泡通過孔板通道后剩余氣泡質(zhì)量的影響，本文在上述三種工況的基礎上，通過改變l1的大小增加了幾組工況。在不同l1和l2的情況下，氣泡穿過孔板后到達計算區(qū)域上方的剩余質(zhì)量進行了統(tǒng)計和比較，如圖10 所示。在上氣泡位置不變的情況下，下氣泡所處位置越高（氣泡之間間距l(xiāng)2越短），兩氣泡在越早的時刻合并，可以穿過孔板到達頂部的質(zhì)量也越多。在氣泡間距相同的情況下，上氣泡與孔板的距離l1越大，通過的剩余質(zhì)量反而越少。對于剩余質(zhì)量，氣泡間的距離比氣泡與孔板間的距離更重要。結合圖10 可以發(fā)現(xiàn)終端速度和剩余質(zhì)量也存在著聯(lián)系：剩余質(zhì)量越多，氣泡受到的浮力越大，其終端速度也越大。

圖10 氣泡初始位置與剩余質(zhì)量的關系Fig. 10 The relationship between the residual mass and the position of the bubble

5 結 論

本文基于大密度比兩相流LBM 模型，對雙氣泡在含孔板微通道內(nèi)的運動過程進行了研究。分別考慮了Eo數(shù)、氣泡之間的相對大小和相對位置等因素對氣泡在孔板通道內(nèi)流動的影響，得到如下結論：

（1）隨著Eo數(shù)的增大，氣泡在通過孔板通道時所包裹的液泡越多，氣泡容易發(fā)生變形、合并和破裂；另外當氣泡在相同的初始位置上升時，隨著Eo數(shù)的增大，氣泡剩余質(zhì)量Ms呈負指數(shù)級地減小。

（2）大氣泡在上、小氣泡在下時，運動過程中氣泡受到的阻力更大，表現(xiàn)為上方氣泡通過孔板的過程中形變更明顯，兩個氣泡在孔口附近合并得更充分，通道內(nèi)的殘留質(zhì)量較少；但到達通道頂部需要的時間更長。大氣泡在下、小氣泡在上時，氣泡到達孔板附近的時間更快，此種情況兩氣泡的形變最小，氣泡到達頂部所需要的時間最短。而當上下兩氣泡相等時，穿過孔板到達通道上半部分的剩余質(zhì)量相對于其他兩種情況最少。

（3）當兩個氣泡初始距離相近時，兩氣泡之間的合并通常發(fā)生在上方氣泡到達孔板通道之前，此時通過的剩余質(zhì)量最多；隨著兩個氣泡之間的距離拉大，兩氣泡之間的合并更容易發(fā)生在上方氣泡經(jīng)過孔板通道時，通過的剩余質(zhì)量次之；當兩氣泡之間的距離增加到一定程度，上方氣泡穿過孔板通道并發(fā)生斷裂之后，下方氣泡才追上并與之合并，到達計算區(qū)域上方的剩余質(zhì)量最小。另外，上下氣泡初始位置同時升高也可以提高通過孔道的剩余質(zhì)量。