功能性多孔介質(zhì)中氣泡輸運調(diào)控的格子Boltzmann分析

王恒博，蘭忠，馬學虎，宋天一，董曉強

（1 大連理工大學化工學院，遼寧大連116024；2 隆華科技集團（洛陽）股份有限公司，河南洛陽471132）

多孔介質(zhì)中的兩相流動廣泛存在于自然界及工業(yè)生產(chǎn)中的各個領域，其中氣泡的輸運特性是相關流動及傳熱過程的關鍵問題。Dai等[1]制備的微孔膜表面如圖1所示，在毛細蒸發(fā)海水淡化技術中，首先利用多孔絕熱層將主體海水隔絕，海水通過微孔膜的毛細作用上升至熱表層從而吸熱蒸發(fā)實現(xiàn)海水的淡化。這一過程中，由于海水中存在一定量的不凝氣，或者在微孔膜層中已經(jīng)開始部分蒸發(fā)，則在微孔道中將析出氣體形成氣泡。由于多孔介質(zhì)骨架結(jié)構(gòu)對氣泡產(chǎn)生的額外阻力[2]，阻礙氣泡的繼續(xù)上升，甚至使氣泡阻塞在通道中，導致海水無法通過被堵塞的微通道。隨著氣泡堵塞率的增大，總體蒸發(fā)效率將大幅下降。那能否僅僅通過加大孔隙率來消除這個問題呢？在絕熱和毛細輸運過程中，需要滿足一定的孔隙率以保持毛細壓力，還要滿足一定的骨架成分以降低熱導率。因此，如何在骨架結(jié)構(gòu)排列較密集的情況下，使氣泡更容易穿過多孔介質(zhì)，從而提高兩相流輸運效率是一個核心的技術問題。

圖1 海水薄膜蒸發(fā)

在多數(shù)多孔材料和微流控設備中，氣泡的運動往往受到流道結(jié)構(gòu)的限制與約束，較多文獻對氣泡在不同幾何結(jié)構(gòu)約束下的運動特征進行了研究。例如，Qin 等[3]通過格子Boltzmann 方法發(fā)現(xiàn)氣泡的合并有助于氣泡更快地穿過多孔介質(zhì)孔隙，同時壁面接觸角的增大不利于氣泡的脫離。Qin 等[4]基于泡沫金屬的SEM 圖像，提出了泡沫金屬二維重構(gòu)模型。采用格子Boltzmann 方法研究了泡沫金屬內(nèi)部的沸騰傳熱機理，表明在低孔隙率下，由于氣泡未產(chǎn)生“清掃”效應而導致在加熱表面形成蒸氣薄膜，從而降低沸騰傳熱性能。Sattari 等[5]通過格子Boltzmann 方法研究了Eo 數(shù)及壁面潤濕性對氣泡通過多孔介質(zhì)的影響，表明低Eo 數(shù)時，氣泡主要受表面張力的影響，形變較小并無法通過多孔介質(zhì)。隨著Eo 數(shù)的增加，氣泡形變增大并能順利通過多孔介質(zhì)。Sun 等[6]通過格子Boltzmann 方法模擬了三維空間加熱壁面上的氣泡動力學，數(shù)值結(jié)果與實驗結(jié)果相吻合，發(fā)現(xiàn)氣泡的合并可以獲得更高的傳熱效率。婁欽等[7]構(gòu)建了非牛頓流體兩相流模型，模擬了多孔介質(zhì)內(nèi)牛頓氣體驅(qū)替非牛頓流體的過程，研究發(fā)現(xiàn)無論被驅(qū)替流體為何種流體，隨著Ca 數(shù)增加，驅(qū)替速度加快，指進現(xiàn)象越明顯，驅(qū)替效率越低。同時在相同孔隙率下，多孔介質(zhì)的形狀對驅(qū)替過程也有較大的影響。婁欽等[8]模擬了復雜微通道內(nèi)氣泡的上升過程，研究了浮力大小、障礙物尺寸、氣泡初始位置等對氣泡變形、合并、分裂的動力學行為，以及對氣泡上升速度、剩余質(zhì)量的運動特性的影響。Xu 等[9]通過格子Boltzmann 方法模擬了土壤孔隙中的兩相流過程，表明Shan-Chen模型可以很好地用于模擬土壤中的水氣分布過程，并且模擬了不同孔隙率、壁面潤濕性下土壤中的水氣分布特征。Yi 等[10]通過格子Boltzmann 方法研究了氣泡通過狹窄孔道的情況，表明氣泡的尺寸以及壁面的潤濕性對液體及氣泡的流動通量影響很大，尤其在通道狹窄處，氣泡流動通量急劇減小。Uchiyama等[11]實驗觀測了氣泡上升過程繞圓柱流動的情況，得到了其通過圓柱障礙物后的形貌變化。Li等[12]通過數(shù)值模擬研究了單氣泡通過圓形障礙物時的情況。研究發(fā)現(xiàn)當氣泡遇見障礙物時，氣泡的形貌會發(fā)生較大的變化。同時，隨著Eo 數(shù)的增加，氣泡離開圓形障礙物的脫離角度也增大。Salcedo等[13]通過粒子成像測速技術（PIV）研究了通道內(nèi)有兩個水平的圓形障礙物對氣泡流動的影響。這些文獻所研究領域涉及微通道、多孔介質(zhì)以及較小的氣泡，在這個尺度上，格子Boltzmann 方法是一個較為理想的數(shù)值模擬工具。目前，研究者更多關注氣泡運動特性本身，并針對較為簡單的通道結(jié)構(gòu)或繞流物時兩相流流動開展研究，而對多孔介質(zhì)模型中，基于多孔介質(zhì)結(jié)構(gòu)和潤濕特性的調(diào)控對氣泡輸運的影響研究較少。因此，本文構(gòu)建了多孔介質(zhì)參數(shù)化模型，通過格子Boltzmann 方法研究了多孔介質(zhì)中氣泡的輸運特性，分別從多孔介質(zhì)孔隙率、孔道潤濕性、多孔介質(zhì)排布以及氣泡水平方向初速度等因素探究氣泡通過多孔介質(zhì)的速度、位移、形貌、與壁面接觸面積和孔隙毛細力變化，為設計合理的孔道結(jié)構(gòu)和分布提供一定的理論依據(jù)。

1 數(shù)值方法

模型采用Qian等[14]提出的偽勢模型。該模型通過對微觀相互作用力的描述，從而反映多相流體動力學的物理本質(zhì)。偽勢模型可以自動追蹤相界面的運動，在模擬復雜的多相流系統(tǒng)中具有較大的優(yōu)勢。在單組分兩相流格子Boltzmann 方法中，流體通過介觀粒子分布函數(shù)進行表征和描述，流體粒子分布函數(shù)在特定結(jié)構(gòu)的格子上進行演化及碰撞。本文采用單松弛模型的流體粒子分布函數(shù)演化方程可以用式(1)表示[15]。

式中，fi(x,t)為密度分布函數(shù)，即t 時刻位于x處的速度為ei的粒子分布函數(shù)；τ 為松弛時間；Δfi(x,t)代表體積力項，包括流體粒子間作用力、壁面黏滯力、重力等。粒子的平衡態(tài)分布函數(shù)fieq(x,t)由式(2)給出。

式中，ωi為權系數(shù)；cs為格子聲速；ρ和u分別為流體的密度和速度。

針對二維模擬，較為成熟的格子模型為D2Q9模型。對D2Q9 格子來說，權系數(shù)ωi分別為4/9（i=0）、1/9（i=1～4）、1/36（i=5～8），離散速度由式(3)給出。

式中，c=δx/δt為格子速度；δx為x方向的格子步長，δt為時間步長，在LBM中兩者通常取1.0。

體積力項Δfi(x,t)采用精確差分格式[16]引入，該方法數(shù)值精度高同時易于實現(xiàn)，可將體積力項直接添加在流體粒子分布函數(shù)演化方程中，見式(4)。

其中Δu=F·δt/ρ 是在時間步δt內(nèi)由于力的作用而產(chǎn)生的速度變化。

在D2Q9 格子中，cs2=c2/3。運動黏度與松弛時間相關，并通過式(5)得到。

流體的密度和速度可通過分布函數(shù)的統(tǒng)計得到［式(6)］。

真實的流體速度U定義為碰撞前后速度的平均值，由式(7)表示[17]。

作用力F是流體所受到作用力的總和，包括粒子間作用力Fint、流體與壁面間作用力Fads和浮力Fg等，可用式(8)表示。

粒子間作用力Fint采用Gong 等[18]提出的綜合加權形式，該模型通過調(diào)節(jié)權重系數(shù)β，可以在不同的狀態(tài)方程下均獲得較為準確的模擬結(jié)果，在本文中β選取為0.886［式(9)］。

式中G表示粒子間作用力的大小。為了使模擬結(jié)果具有各向同性[19]，取值為［式(10)］。

式中，g0為粒子間作用強度系數(shù)。

ψ(x)用來模擬流體粒子間相互作用的偽勢函數(shù)，可用式(11)表示。

式中壓力p 采用Shan-Chen 狀態(tài)方程［式(12)］[20]。

式中，ρ0=1.0。在D2Q9 模型中c0=6.0。令壓力與密度的一階和二階導數(shù)為零，可得到狀態(tài)方程的臨界參數(shù)。對Shan-Chen 狀態(tài)方程，臨界參數(shù)見式(13)[21]。

在Shan-Chen 狀態(tài)方程中，溫度定義T=-1/g0，臨界溫度即為Tcr=4.5ρ0。

固體壁面的潤濕性在LBM 偽勢兩相流模型中較易于實施，是通過引入固體壁面與流體之間作用力Fads來實現(xiàn)的［式(14)］[22]。

式中，Gads用來描述流體與固體壁面間作用強度，用于調(diào)控壁面潤濕性；s(x)為指示函數(shù)：x為固體格點時，s(x)為1；x為流體格點時，s(x)為0。

浮力Fg可用式(15)表示。

2 模型驗證

根據(jù)本文所涉及的多孔介質(zhì)孔道內(nèi)氣泡輸運特性的研究，相關的特征物理規(guī)律包括表面張力、潤濕特性以及氣泡運動特性等。因此，首先在這三個方面對構(gòu)建的格子玻爾茲曼方法進行模型驗證。

2.1 Laplace定律

首先采用Laplace 定律來驗證模型對表面張力模擬的正確性。根據(jù)Laplace 定律可知，當液滴與周圍氣體處于熱平衡和氣平衡的條件下，表面張力σ恒定，且液滴內(nèi)外壓差Δp與半徑的倒數(shù)1/R呈線性關系［式(16)］。

定義一個200×200的計算區(qū)域。初始時刻在中心處設置一定半徑的液滴，其余區(qū)域充滿了T0=0.80Tcr溫度下對應的飽和氣體，初始密度為該溫度下的飽和密度。四周采用周期性邊界條件。經(jīng)過30000 步迭代后計算達到穩(wěn)定。圖2 給出了不同半徑下液滴內(nèi)外壓力差與半徑的關系。擬合方程的斜率為0.0687（Adj.R2=0.996），可以得出計算流體的表面張力為0.0687。

圖2 液滴內(nèi)外壓力差與半徑的關系

2.2 接觸角驗證

為了考察固體壁面接觸角對氣泡運動行為的影響，需要引入表面與流體之間的相互作用。在偽勢模型中，通過式(14)引入的固體壁面與流體之間的作用力Fads來表示流體的潤濕性，其中改變流體與固體壁面間作用力強度Gads可以實現(xiàn)不同接觸角的模擬。

定義一個200×200的計算區(qū)域，初始時刻在下壁面中心（100,25）處給定一個半徑R0=25的液滴，其余區(qū)域充滿了T0=0.80Tcr溫度下的飽和氣體。計算流體的表面張力σ=0.0687，液相密度ρl=1.90，氣相密度ρg=0.14。計算域上下邊界采用半步長反彈格式，左右邊界為周期性邊界。液滴在固體壁面與流體間的作用下不斷演化，直至達到平衡，形成固定的接觸角。圖3 顯示了不同固液作用參數(shù)下，液滴穩(wěn)定時的接觸角。擬合得到方程：θ=212.23-456.07Gads（Adj.R2=0.999），該結(jié)果與文獻[23]中固液作用參數(shù)與接觸角呈線性相關的規(guī)律相符。

圖3 平衡時接觸角與固液作用參數(shù)的關系

2.3 單氣泡上升過程的形貌特征

定義一個240×400的計算區(qū)域，初始時刻在中心（120,70）處設置一個半徑R0=28 的氣泡，其余區(qū)域充滿了T0=0.80Tcr溫度下對應的飽和液體，初始密度為該溫度下的飽和密度。計算流體的表面張力σ=0.0687，液相密度ρl=1.90，氣相密度ρg=0.14。計算域四周采用周期性邊界。在浮力作用下，氣泡開始上浮運動。隨著氣泡的上浮過程，上浮速度和形貌特征不斷變化，最終在浮力以及表面張力的作用下趨于穩(wěn)定。其中，氣泡形貌主要與莫頓數(shù)Mo以及厄特沃什數(shù)Eo有關[24]，其中Mo=g(ρl-ρg)μ4l/ρ2lσ3代表黏性力和表面張力的比值，Eo=g(ρl-ρg)d2/σ 代表浮力與表面張力的比值。

圖4展現(xiàn)了3組不同Eo數(shù)及Mo數(shù)下，氣泡穩(wěn)定后的形貌特征。模擬結(jié)果與文獻[25-26]結(jié)果相吻合，再次驗證了此模型可以正確模擬氣泡的輸運特性。

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

3.1 物理模型

圖4 不同條件下氣泡達到穩(wěn)定狀態(tài)的最終形貌模擬結(jié)果（左）與文獻［25］（中）、文獻［26］（右）的對比

基于毛細蒸發(fā)過程傳遞特性優(yōu)化背景，絕熱多孔層需要滿足一定的孔隙率以保持毛細壓力。另一方面，應調(diào)控多孔介質(zhì)骨架結(jié)構(gòu)參數(shù)以減小氣泡所受阻力。將多孔介質(zhì)簡化為由柱狀骨架間隔形成的通道[27]，模型如圖5 所示，并以結(jié)構(gòu)參數(shù)來調(diào)整孔隙率，見表1。采用這種簡化物理模型，便于模擬計算以及結(jié)果分析。

為了驗證網(wǎng)格無關性，以單氣泡上升過程速度為參考（Eo=10, Mo=0.4535），分別考察網(wǎng)格數(shù)為170×280、240×400、340×560、480×800 時氣泡速度的變化。如圖6 所示，當計算域為170×280 時，網(wǎng)格數(shù)目較少，此時氣泡速度誤差較大，誤差存在于氣泡上升過程受表面張力和浮力的作用下速度逐漸減小至穩(wěn)定的階段。隨著網(wǎng)格數(shù)的加密，對氣泡運動過程的描述更加準確，誤差逐漸減小。當網(wǎng)格加密至240×400 時，與480×800 網(wǎng)格數(shù)計算結(jié)果最大誤差為2.72%，平均誤差為0.72%。因此，為了消除網(wǎng)格數(shù)對計算的影響，同時基于計算資源的考慮，本文選取網(wǎng)格數(shù)240×400進行計算。

圖5 物理模型

表1 多孔介質(zhì)結(jié)構(gòu)參數(shù)

圖6 網(wǎng)格無關性檢驗

計算域為240×400 的矩形區(qū)域，下部中心（120, 70）處給定一個半徑R0=28 的氣泡，其余區(qū)域充滿了T0=0.80Tcr溫度下的飽和液體。計算流體的表面張力σ=0.0687，液相密度ρl=1.90，氣相密度ρg=0.14。氣泡量綱為1數(shù)Eo數(shù)及Mo數(shù)分別選取為Eo=10、Mo=0.4535，可計算得到重力加速度g=0.00007。計算域四周采用周期性邊界條件，多孔介質(zhì)壁面處采用半步長反彈格式，以保證整個流場在計算過程中保持質(zhì)量守恒。松弛因子τ=1.0，根據(jù)式(5)計算得到流體的動力學黏度ν=1/6。

在下文數(shù)據(jù)分析中，氣泡區(qū)域的密度ρ定義為小于(ρl+ρg)/2 的區(qū)域。氣泡在上浮過程中沿著豎直方向上的瞬態(tài)速度定義見式(17)。

式中，V 代表整個氣泡區(qū)域。量綱為1 時間定義為t*=t(g/2R0)1/2。其中t為格子時間，g為重力加速度，R0為氣泡初始時的半徑。

氣泡在多孔介質(zhì)間隙的運動過程中，所受毛細力Fcap用式(18)表示。

圖7 多孔介質(zhì)孔隙率ε=0.96時氣泡輸運特性

式中，σ為液體的表面張力；θ為壁面接觸角；r為多孔介質(zhì)的孔道間隙。

3.2 孔隙率的影響

本小節(jié)主要研究多孔介質(zhì)孔隙率對氣泡輸運特性的影響，分別考察了孔隙率ε=0.99、0.96、0.93、0.91、0.88、0.84的情況。其他模擬參數(shù)設置如下：Eo=10，Mo=0.4535，接觸角θ=30°。

圖7給出了多孔介質(zhì)孔隙率ε=0.96時氣泡上升過程的速度v 及形貌變化。其中藍色區(qū)域為氣泡，紅色區(qū)域表示該溫度下的飽和液體，灰色區(qū)域為多孔介質(zhì)骨架結(jié)構(gòu)。氣泡穿過多孔介質(zhì)孔道時的速度呈現(xiàn)出周期性的振蕩，并且速度的極小值均出現(xiàn)在穿過孔道最狹窄的結(jié)構(gòu)時，與文獻[8]中氣泡穿過障礙物時得到的結(jié)果類似。具體過程如下：演化初期階段，在浮力的作用下，氣泡速度逐漸增大。在t*=1.5 時，氣泡上表面接觸多孔介質(zhì)，由于孔道截面積變窄，導致氣泡運動受到阻礙，其速度減小。在t*=2.6 時，氣泡上表面剛通過第1 層孔道，氣泡所受阻力減小，速度增大出現(xiàn)一個極大值。在t*=4.3 時，整個氣泡完全進入多孔介質(zhì)中，此時氣泡速度最小。當氣泡依次穿過孔道時，氣泡速度呈現(xiàn)出周期性的振蕩，同時氣泡形貌在豎直方向上拉長。在t*=14.0時，氣泡上部分穿出多孔介質(zhì)區(qū)域，所受阻力減小，使氣泡速度逐漸增大。在t*=16.3時，氣泡完全穿出多孔介質(zhì)，此時氣泡速度為整個過程中最大值。隨后在表面張力和浮力的作用下，氣泡開始收縮變得扁平，速度逐漸減小至平穩(wěn)，穩(wěn)定后形貌與圖4中結(jié)果相吻合。

圖8為不同多孔介質(zhì)孔隙率時氣泡上浮過程速度變化。隨著孔隙率的減小，孔道截面積減小，導致氣泡所受阻力增大，因此氣泡穿過多孔介質(zhì)所用時間更長。同時，由于通道截面積變窄，使得氣泡在上升過程中形變程度變大，反映在圖中即氣泡速度的周期性振蕩幅度增大。

圖8 不同多孔介質(zhì)孔隙率時氣泡上浮速度變化

圖9 不同多孔介質(zhì)孔隙率時流場分布

圖9為不同多孔介質(zhì)孔隙率時的流場分布。在氣泡上升過程中，氣泡的兩側(cè)產(chǎn)生對稱的渦流，此處流體運動較為劇烈。在遠離氣泡的部位流場分布較為平緩。隨著氣泡上浮過程中速度越來越大，氣泡的形變越大，渦旋結(jié)構(gòu)的尾跡也越長。而由于多孔介質(zhì)的存在，在靠近多孔介質(zhì)的部分兩相流流動發(fā)生改變，流線發(fā)生彎曲，從而影響到氣泡的形貌及速度變化。隨著孔隙率的減小，孔道截面積減小，導致多孔介質(zhì)壁面附近流線彎曲程度更大。氣泡兩側(cè)渦流結(jié)構(gòu)也隨孔隙率的減小而變小，尾跡逐漸變短，從而導致氣泡的上浮速度減小，同時產(chǎn)生更大的形變。

圖10 為不同多孔介質(zhì)孔隙率時氣泡通過多孔介質(zhì)過程中與壁面的平均接觸面積。隨著孔隙率的減小，多孔介質(zhì)骨架尺寸增大，使得氣泡與壁面的接觸面積增大。在毛細輸運過程中，需要降低熱導率以防止氣泡吸收熱量后體積變大，從而導致阻塞孔道。因此，多孔介質(zhì)孔隙率越小，氣泡與壁面?zhèn)鳠崦娣e越小，越有利于毛細蒸發(fā)海水淡化過程。

圖10 氣泡通過不同多孔介質(zhì)孔隙率時與壁面平均接觸面積

圖11 為不同多孔介質(zhì)孔隙率時氣泡上浮位移隨時間的變化。當孔隙率較大時，即ε=0.99、0.96，氣泡在多孔介質(zhì)中上升位移接近勻速運動，達到多孔介質(zhì)頂端所用時間較少。而隨著孔隙率的減小，相同時間內(nèi)氣泡上升位移也越小，達到頂端所用時間增多。

圖11 不同多孔介質(zhì)孔隙率時氣泡上浮位移隨時間變化

圖12 為氣泡通過多孔介質(zhì)所用時間與氣泡所受毛細力隨不同孔隙率的變化。隨著孔隙率的減小，氣泡通過多孔介質(zhì)所用時間更長。另一方面，隨著孔隙率的減小，由于孔道間隙r減小，毛細力逐漸增大。因此，為使氣泡更快脫離多孔介質(zhì)，同時使多孔介質(zhì)絕熱層保持一定的毛細力，多孔介質(zhì)孔隙率應在0.88～0.96范圍內(nèi)。

圖12 不同孔隙率時氣泡通過多孔介質(zhì)所用時間與所受毛細力的變化

3.3 壁面潤濕性的影響

圖13 不同潤濕性壁面中氣泡輸運特性

壁面潤濕性對兩相流流動有著十分重要的影響。圖13 分別為壁面接觸角θ=30°、90°、120°時氣泡在多孔介質(zhì)中的輸運特性。在圖13(a)中，當氣泡通過親水壁面時，由于孔道壁面呈現(xiàn)親水疏氣的效應，氣泡與壁面接觸面積較小，因此氣泡可以較為快速地穿過多孔介質(zhì)。在圖13(b)中，壁面接觸角增加至90°，在t*=2.2 時，氣泡上表面接觸多孔介質(zhì)壁面，壁面呈現(xiàn)出對氣泡的黏附效應，使氣泡上表面與壁面接觸部分出現(xiàn)凸起形變。在t*=13.0時，由于氣泡與壁面接觸面積增大，壁面的黏附效應更明顯，使氣泡在豎直方向上產(chǎn)生較大的形變。壁面的黏附效應增加了氣泡運動過程的阻力，因此氣泡通過多孔介質(zhì)所用的時間則更長。在圖13(c)中，壁面接觸角增加至120°，壁面呈現(xiàn)疏水親氣的效應。此時氣泡將部分壁面包裹至其中，因此壁面對氣泡的黏附效應進一步增強，使得氣泡完全被多孔介質(zhì)所“捕獲”，無法脫離多孔介質(zhì)。

圖14 為不同潤濕性壁面時多孔介質(zhì)內(nèi)流場分布。當壁面接觸角θ從30°增加至60°時，壁面對氣泡的黏附效應增強，導致氣泡通過多孔介質(zhì)時的形變更大，從而使其渦流結(jié)構(gòu)也變大。當接觸角增加至90°時，壁面開始呈現(xiàn)疏水親氣效應，增大了氣泡上升過程中的阻力，使形變進一步增大，渦流尾跡也變長，對流場影響較大。當接觸角增加至120°時，氣泡無法克服壁面的黏附效應，被捕獲在多孔介質(zhì)間，此時渦流集中在氣泡兩側(cè)，對流場影響較小。

圖15 為氣泡通過不同潤濕性壁面時與壁面平均接觸面積。隨著壁面潤濕性增強，壁面對氣泡的黏附效應增強，使得氣泡與壁面的接觸面積增大。因此，在毛細蒸發(fā)海水淡化過程中，親水壁面有助于降低熱導率以使氣泡快速通過。

圖14 不同潤濕性壁面時多孔介質(zhì)內(nèi)流場分布

圖15 氣泡通過不同潤濕性壁面時與壁面平均接觸面積

圖16 不同潤濕性壁面中氣泡上升位移隨時間的變化

圖16 為不同潤濕性壁面中氣泡上升位移隨時間的變化。當θ=30°、60°時，壁面呈現(xiàn)親水疏氣的效應，此時氣泡在多孔介質(zhì)中的上浮過程接近勻速運動。當接觸角θ=80°、90°時，壁面的黏附效應逐漸加強，增加了氣泡上浮過程的阻力，使得氣泡通過多孔介質(zhì)所需增長。隨著接觸角的進一步增大，當θ=100°、120°時，由于壁面的疏水親氣效應進一步增強，對氣泡的黏附效應也逐漸增強，最終使得氣泡無法通過多孔介質(zhì)。

圖17 為氣泡通過多孔介質(zhì)時間與所受毛細力隨壁面接觸角的變化。研究發(fā)現(xiàn)隨著壁面接觸角的增大，由于壁面的黏附效應，氣泡通過多孔介質(zhì)所用時間增加。另一方面，隨著接觸角的增大，孔道毛細壓力隨之減小。因此，選取壁面接觸角θ=30°，可以使氣泡在最快脫離多孔介質(zhì)的同時，使多孔介質(zhì)孔道的毛細壓力也保持最大，從而最有利于毛細蒸發(fā)過程的強化。

圖17 不同潤濕性壁面中氣泡通過多孔介質(zhì)所需時間與所受毛細力的變化

3.4 孔道排布的影響

探究了多孔介質(zhì)孔道排布對氣泡輸運特性的影響。對比圖18(a)、(b)，當多孔介質(zhì)孔隙率和接觸角相同時，多孔介質(zhì)孔道排布方式對氣泡的輸運行為有著較大的影響。當多孔介質(zhì)孔道間錯排列時，多孔介質(zhì)骨架結(jié)構(gòu)阻礙在氣泡的上浮路徑中，使得氣泡必須穿過部分骨架結(jié)構(gòu)，同時氣泡與壁面的接觸面積也更大。因此與順次排列的結(jié)構(gòu)相比，氣泡在間錯排列的多孔介質(zhì)中所受到的阻力更大，導致氣泡穿過多孔介質(zhì)所用時間更長。在圖18(c)中，當接觸角增大至60°時，由于氣泡與多孔介質(zhì)的接觸面積更大，壁面對氣泡的黏附效應也更加顯著，因此導致氣泡在豎直方向上產(chǎn)生了更大的形變。同時由于壁面對氣泡的阻力進一步增強，氣泡被“捕獲”在孔道間隙中。在圖18(d)中，當接觸角進一步增大至90°時，黏附效應更加顯著。而在3.3 的模擬中，只有在接觸角大于100°時，氣泡才會被“捕獲”在多孔介質(zhì)中，無法繼續(xù)上升。這是因為當多孔介質(zhì)孔道間錯排列時，由于氣泡與壁面的接觸部位更多，使得壁面的黏附效應更為顯著。因此說明孔道排布方式對氣泡輸運特性的影響最為顯著，孔道結(jié)構(gòu)越復雜，壁面的黏附效應越顯著，氣泡則越難以穿過多孔介質(zhì)。

圖18 多孔介質(zhì)孔道排布方式對氣泡輸運特性的影響

圖19為不同孔道排布時多孔介質(zhì)內(nèi)流場分布。相比于順次排列，圖19(b)間錯排列的多孔介質(zhì)孔道更加復雜，因此在靠近多孔介質(zhì)的地方，兩相流的流場所受限制更大，流線也更加彎曲。同時，氣泡通過間錯排列多孔介質(zhì)所受阻力更大，因此導致其速度減小，渦流結(jié)構(gòu)也減小。在圖19(c)中，壁面接觸角增加至60°，此時壁面對氣泡的黏附效應增強，使形變進一步增大，渦流結(jié)構(gòu)變大，與圖14壁面潤濕性中流場分布規(guī)律一致。

圖19 不同孔道排布時多孔介質(zhì)內(nèi)流場分布

相同孔隙率下，當孔道排布方式為順次排列時，氣泡通過多孔介質(zhì)過程中與壁面的平均接觸面積為35.43。當孔道排布為間錯排列時，平均接觸面積為35.73。由此表明，孔道排布方式對氣泡與壁面的平均接觸面積影響較小。

圖20 為不同多孔介質(zhì)孔道排布方式時氣泡上升位移隨時間變化。在相同孔隙率和壁面潤濕性的情況下，氣泡通過間錯排列的多孔介質(zhì)所需時間更長。同時，由于氣泡上升過程中需要穿過多孔介質(zhì)骨架結(jié)構(gòu)，因此其位移變化呈現(xiàn)出幾個不同的階段。

圖20 不同多孔介質(zhì)孔道排布方式氣泡上升位移隨時間變化

3.5 水平方向初速度的影響

在流動沸騰過程中，氣泡生成時存在一定的水平方向初速度。本節(jié)模擬了不同水平方向初速度v0=0.01、0.02、0.04、0.06、0.08、0.10、0.12 時，氣泡在多孔介質(zhì)中的輸運特性。其他模擬參數(shù)設置如下：Eo=10，Mo=0.4535，接觸角θ=30°，孔隙率ε=0.99.

圖21(a)為水平方向初速度v0=0.02 時，氣泡在多孔介質(zhì)中輸運特性。由于存在水平方向的速度，導致氣泡在遇到金屬骨架的阻礙時，改變其運動路徑，從而可以從孔道間隙穿過。因此不需要直接克服骨架的阻力，使得氣泡所受阻力減小。與圖18(b)相比，在相同孔隙率和接觸角下，氣泡穿過多孔介質(zhì)所用時間也更短。圖21(b)為水平方向初速度v0=0.06 時的情況。氣泡的水平方向速度更大，因此水平方向上的位移更大，使氣泡更容易在孔道間隙運動，減小了需要克服骨架的阻力，氣泡穿過多孔介質(zhì)所用時間也更短。當初速度增加至v0=0.12時，由于水平方向初速度過大，導致氣泡的實際運動路線接近斜向方向上的運動。氣泡需要克服骨架的阻力，阻礙了氣泡的運動。

圖22 為不同水平方向初速度時多孔介質(zhì)內(nèi)流場分布。當存在水平方向初速度時，氣泡運動方向發(fā)生改變，氣泡兩側(cè)渦流不再對稱，呈現(xiàn)出一定的傾角。同時，渦流形狀也不再相同，左側(cè)渦流尾跡短于右側(cè)渦流尾跡。并且隨著水平方向初速度增大，兩側(cè)渦流差別越明顯，所呈現(xiàn)出的傾角也越大。

圖23 為不同水平方向初速度的氣泡與壁面平均接觸面積。結(jié)果表明隨著水平方向初速度的增大，氣泡與壁面平均接觸面積增加。因此，在毛細蒸發(fā)海水淡化過程中，為了降低熱導率以使氣泡快速通過，氣泡水平方向初速度應處在較低的范圍內(nèi)。

圖24 為不同水平方向初速度時氣泡上浮速度隨時間的變化。當水平方向初速度v0≤0.10 時，隨著水平方向初速度的增大，氣泡穿過多孔介質(zhì)所用時間更短，但時間縮短的程度越來越小。當水平方向初速度v0=0.12 時，由于氣泡運動方向的變化，使得氣泡需要克服更大的阻力，因此穿過多孔介質(zhì)時間增長。因此，一定的水平方向初速度有助于改變氣泡運動軌跡，減小其運動過程所受阻力，縮短氣泡穿過多孔介質(zhì)所用時間，從而有利于毛細蒸發(fā)過程。而當水平方向初速度超過一定值后，反而不利于氣泡在多孔介質(zhì)中的運動。

圖21 不同水平方向初速度對氣泡輸運特性的影響

圖22 不同水平方向初速度時多孔介質(zhì)內(nèi)流場分布

圖23 不同水平方向初速度的氣泡與壁面平均接觸面積

圖24 不同水平方向初速度時氣泡上浮速度隨時間的變化

圖25 為不同水平方向初速度時氣泡上浮位移隨時間變化。隨著水平方向上初速度增大，氣泡在相同時間內(nèi)上升位移增大。當水平方向上速度增大至0.12時，由于明顯改變了氣泡的運動軌跡，因此氣泡所受阻力增大，通過多孔介質(zhì)所用時間更長。將不同水平方向上氣泡通過多孔介質(zhì)所用時間繪制于圖26，研究發(fā)現(xiàn)，隨著水平方向上初速度增大，氣泡所用時間減小，但減小幅度越來越弱。當速度增大至0.12時，所用時間增大。

圖25 不同水平方向初速度時氣泡上浮位移隨時間的變化

4 結(jié)論

在毛細蒸發(fā)過程中，為維持一定的毛細壓力，多孔介質(zhì)孔隙率不宜過大。另一方面則需要調(diào)控多孔介質(zhì)骨架結(jié)構(gòu)使氣泡更快速地通過。本文構(gòu)建了功能性多孔介質(zhì)結(jié)構(gòu)，通過格子Boltzmann 模型對氣泡輸運調(diào)控提出策略，得出以下結(jié)論。

圖26 不同水平方向初速度時氣泡通過多孔介質(zhì)所用時間

（1）隨著多孔介質(zhì)孔隙率的減小，孔道截面積減小，氣泡兩側(cè)渦流結(jié)構(gòu)也隨孔隙率的減小而變小，氣泡穿過多孔介質(zhì)所用時間更長。同時，多孔介質(zhì)孔隙間毛細力增大，氣泡與壁面?zhèn)鳠崦娣e減小，有益于毛細蒸發(fā)過程。因此，為使氣泡更快脫離多孔介質(zhì)，同時保持一定的毛細力，多孔介質(zhì)孔隙率應在0.88～0.96范圍內(nèi)。

（2）氣泡通過親水壁面多孔介質(zhì)時所用時間較少，與壁面接觸面積也較小。同時，多孔介質(zhì)孔隙間毛細力較大。隨著壁面接觸角的增大，壁面對氣泡的黏附效應顯著使形變進一步增大，渦流尾跡也變長，對流場影響較大。當壁面接觸角大于100°時，氣泡則被“捕獲”在多孔介質(zhì)中。因此，選取接觸角θ=30°最適宜毛細蒸發(fā)過程。

（3）間錯排列的孔道結(jié)構(gòu)會顯著加劇壁面的黏附效應，氣泡越容易被多孔介質(zhì)所“捕獲”。在靠近多孔介質(zhì)的地方，流線也更加彎曲，渦流結(jié)構(gòu)也減小。因此，對于多孔膜設計時，宜考慮有序孔道的要求。

（4）當氣泡存在一定的水平方向初速度時，氣泡的運動路徑發(fā)生改變，可以更快地通過多孔介質(zhì)結(jié)構(gòu)。同時，氣泡兩側(cè)渦流不再對稱，渦流形狀也不再相同。但當水平方向初速度增加到一定值后，氣泡所受阻力又增大。本文分析的條件范圍顯示，對于自由上升過程的氣泡而言，水平初速在v0=0.06～0.10為宜。

符號說明

c—— 格子速度，m/s

cs—— 格子聲速，m/s

d0—— 氣泡直徑，m

e—— 格子速度矢量，m/s

F—— 體積力項，N

Fads—— 固體壁面與流體間作用力，N

Fcap—— 毛細壓力，N

Fg—— 浮力，N

Fint—— 粒子間作用力，N

Gads—— 流體與固體壁面作用系數(shù)

g—— 重力加速度，m/s2

g0,gc—— 粒子間作用強度系數(shù)

PPI—— 孔密度

p—— 壓力，N

R0—— 初始氣泡直徑，m

r—— 孔道間隙，m

S—— 氣泡與壁面平均接觸面積

T0—— 初始溫度，K

Tcr—— 臨界溫度，K

t—— 格子時間，s

t*—— 量綱為1時間

x—— 流體格點

β—— 權重系數(shù)

δt—— 時間步長，s

ε—— 孔隙率

θ—— 液滴表觀接觸角，(°)

ρ,ρ0,ρc—— 流體密度，kg/m3

σ—— 表面張力，N

τ—— 松弛因子

υ—— 流體動力學黏度，N·s/m2

Ψ—— 勢函數(shù)

角標

0—— 初始值

c—— 臨界參數(shù)

l—— 液體

v—— 氣體

*—— 該物理量經(jīng)過量綱為1化處理