臧晨強(qiáng) 婁欽

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

復(fù)雜微通道內(nèi)非混相驅(qū)替過程的格子Boltzm ann方法?

臧晨強(qiáng) 婁欽?

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

(2017年1月10日收到;2017年5月4日收到修改稿)

本文采用改進(jìn)的基于偽勢(shì)模型的格子Boltzmann方法研究復(fù)雜微通道內(nèi)的非混相驅(qū)替問題.這種方法克服了原始偽勢(shì)模型中計(jì)算結(jié)果對(duì)網(wǎng)格步長的依賴.首先用Lap lace定律驗(yàn)證模型的正確性,然后用該方法研究壁面潤濕性、粗糙結(jié)構(gòu)、黏性比以及距離對(duì)非混相驅(qū)替過程的影響.模擬結(jié)果表明：與壁面粗糙結(jié)構(gòu)和黏性比相比,壁面潤濕性的影響是決定性的因素.隨著接觸角的增加,驅(qū)替效率增加,當(dāng)接觸角大于某一值后,驅(qū)替效率不再變化;隨著黏性比的增加,驅(qū)替效率增加;而壁面粗糙性對(duì)驅(qū)替過程的影響較復(fù)雜,只有凸起半圓的半徑在一定范圍內(nèi)增加時(shí),驅(qū)替效率增加;距離較小時(shí)將促進(jìn)驅(qū)替過程.

非混相驅(qū)替,格子Boltzmann方法,驅(qū)替效率,驅(qū)替時(shí)間

1 引 言

非混相驅(qū)替過程是一個(gè)重要的研究課題,廣泛存在于日常生活和工業(yè)生產(chǎn)中,例如土壤中的非水相液體的運(yùn)輸、聚合物電解質(zhì)膜燃料電池的氣體擴(kuò)散層中水的輸運(yùn)、石油的采集等.該問題是一個(gè)典型的接觸線問題,涉及復(fù)雜的流體-流體間相互作用力和流體-固體間相互作用力.因此,盡管很多學(xué)者關(guān)注到了這類問題[1-21],但是目前的研究尚不完善.

一些研究者采用宏觀方法和實(shí)驗(yàn)方法研究非混相驅(qū)替問題.Tecklenburg等[1]采用雙連續(xù)模型研究水平裂隙巖體中非混相驅(qū)替過程.Zhu等[2]采用流體體積(VOF)方法研究聚合物電解質(zhì)燃料電池氣體通道中液態(tài)水的驅(qū)替過程,他們研究了壁面的潤濕性、空氣入口速度、注水速度和孔隙的大小對(duì)驅(qū)替過程的影響.Yang等[3]用實(shí)驗(yàn)的方法研究了毛細(xì)管力驅(qū)動(dòng)的液-液驅(qū)替問題,重點(diǎn)研究了相對(duì)黏度對(duì)驅(qū)替效率的影響.Islam et等[4]研究了親水性顆粒懸浮液中水的驅(qū)替問題.在他們的研究中,主要考察不同體積的水滴對(duì)驅(qū)替過程的影響.李維仲等[5]將實(shí)驗(yàn)方法和格子Boltzmann方法結(jié)合,研究了單個(gè)氣泡在具有三個(gè)半圓形喉部的復(fù)雜流道內(nèi)的上升過程.

盡管用宏觀方法和實(shí)驗(yàn)方法研究非混相驅(qū)替過程可以得到流動(dòng)的宏觀參數(shù),如驅(qū)替效率等,但得不到流動(dòng)的細(xì)節(jié).為了彌補(bǔ)該不足,一些研究者用微觀方法研究驅(qū)替問題[6,7].Primkulov[6]利用分子動(dòng)力學(xué)理論研究了微米級(jí)球形玻璃表面水滴的驅(qū)替問題.在他們的研究中,考察了毛細(xì)數(shù)、流體黏度和界面特性對(duì)驅(qū)替過程的影響.數(shù)值計(jì)算結(jié)果表明,分子動(dòng)力學(xué)模型很好地?cái)M合了毛細(xì)數(shù)大于10-4的實(shí)驗(yàn)數(shù)據(jù).利用分子動(dòng)力學(xué)方法,Kop lik等[7]研究了低雷諾數(shù)下的非混相驅(qū)替過程.研究表明,在低雷諾數(shù)的情況下,無滑移條件在接觸線附近不適用.雖然流動(dòng)的細(xì)節(jié)可以通過微觀方法獲得,但由于微觀方法的計(jì)算量較大,因此僅適用于較小的計(jì)算區(qū)域.

近年來,作為連接微觀方法和宏觀方法橋梁的介觀方法開始被諸多學(xué)者采用,格子Boltzmann方法作為介觀方法的一種開始受到廣泛關(guān)注.該方法的優(yōu)點(diǎn)是其有著宏觀方法的本質(zhì)和微觀方法的特性,因此它能夠直接處理流體與流體間的作用力及流體與固壁間的作用力.目前基于格子Boltzm ann方法的非混相驅(qū)替問題的研究越來越多.Kang等[12,13]在二維和三維通道里研究了不混溶液滴受重力作用下的驅(qū)替流動(dòng),包括壁面潤濕性、Bond數(shù)、液滴大小和密度比對(duì)驅(qū)替流體的影響.之后,他們又研究了驅(qū)替過程中的指進(jìn)現(xiàn)象并得到了一些重要的結(jié)果[14].然而他們的研究主要集中在沒有孔隙介質(zhì)的單通道的驅(qū)替問題.Huang等[15]采用基于顏色梯度的多相格子Boltzmann方法,研究了毛細(xì)數(shù)和黏性比對(duì)非均相多孔介質(zhì)中非混相驅(qū)替過程的影響.Dong等[16,17]采用格子Boltzmann方法研究了微通道和多孔介質(zhì)中毛細(xì)數(shù)、Bond數(shù)和黏性比對(duì)非混相驅(qū)替過程的影響.李維仲等[18]采用格子Boltzm ann方法研究了水平通道內(nèi)流體黏度比以及壁面潤濕性對(duì)非混相驅(qū)替過程的影響.李娟等[19]采用多相多組分格子Boltzm ann模型中的偽勢(shì)模型研究了單通道中壁面上液滴驅(qū)替過程.彭本利等[20]采用自由能格子Boltzmann方法研究了在不同蒸汽速度剪切作用下,液滴在具有不同潤濕性固體表面上的變形和運(yùn)動(dòng)過程.以上研究推動(dòng)并揭示了一些驅(qū)替機(jī)理,但對(duì)于一些典型的多孔結(jié)構(gòu),如孔喉結(jié)構(gòu),其流動(dòng)細(xì)節(jié)卻不清楚.為了得到孔喉結(jié)構(gòu)中的非混相驅(qū)替問題的詳細(xì)信息,Liang等[21]利用數(shù)值模擬方法研究了含腔微通道內(nèi)的非混相驅(qū)替問題.在該工作中,主要研究了壁面潤濕性、毛細(xì)數(shù)和密度比對(duì)驅(qū)替效率的影響,然而沒有考慮微通道的粗糙性.本文在Liang等的工作基礎(chǔ)上,研究了含粗糙壁面的微通道內(nèi)的驅(qū)替問題,其中微通道的粗糙結(jié)構(gòu)由一個(gè)凸起的半圓來體現(xiàn),在研究中采用了格子Boltzmann方法,并主要考察了壁面潤濕性、粗糙度、黏性比和粗糙表面與含液相腔間的距離對(duì)非混相驅(qū)替過程的影響.

目前有很多常用的多相格子Boltzmann模型被相繼提出,例如顏色梯度模型[22]、偽勢(shì)模型[23,24]、自由能模型[25,26]、動(dòng)力學(xué)模型[27-29].在這些方法中,偽勢(shì)模型由于其簡(jiǎn)單性成為最常用的方法之一.然而它也有一些局限性,例如Yu等[30]發(fā)現(xiàn)流體性質(zhì)對(duì)網(wǎng)格步長有依賴性,并提出了一個(gè)改進(jìn)的格子Boltzmann方法來克服原模型的缺陷.本文采用了這種改進(jìn)的格子Boltzmann方法來研究復(fù)雜通道內(nèi)的非混相驅(qū)替問題.

2 數(shù)值方法

本文采用了Yu等[30]提出的改進(jìn)的偽勢(shì)模型進(jìn)行數(shù)值模擬研究.這種模型具有簡(jiǎn)單性和可靠性.由于該模型在文獻(xiàn)[30]中已經(jīng)被詳細(xì)描述,因此在這里我們只給出一個(gè)簡(jiǎn)短的描述.在這個(gè)模型中粒子的分布函數(shù)fi的演化方程為

其中i=0,1,2,···,b-1,b是離散速度的個(gè)數(shù);x和t分別表示是位置和時(shí)間;ci表示離散速度;τ是松弛時(shí)間,它與運(yùn)動(dòng)黏度ν的關(guān)系為(δt表示時(shí)間步長),是一個(gè)模型常數(shù),其中c=δx/δt(δx代表網(wǎng)格步長).是平衡態(tài)分布函數(shù),且與局部密度ρ和速度u有關(guān)：

其中,wi是權(quán)重系數(shù).方程(1)中的Fi是力項(xiàng),它的表達(dá)式為[31]

其中F是總的相互作用力,它包含三部分：流體間相互作用力(Fint)、流固作用力(Fads)和外力(G).在本文使用的模型中,導(dǎo)致相分離的流體間作用力F in t為[24,30]

其中g(shù)是一個(gè)常數(shù),表示流體間作用力強(qiáng)度.ψ(x)是“有效質(zhì)量”,它是與局部密度有關(guān)的一個(gè)函數(shù)：

流固作用力Fads的表達(dá)式為[32]

其中ρw是壁面密度;s(x)是一個(gè)開關(guān)函數(shù),當(dāng)x是固體點(diǎn)時(shí),s(x)取值為1,而當(dāng)x為流體點(diǎn)時(shí),s(x)取值為0.宏觀物理量的表達(dá)式為：

本文使用D2Q9模型進(jìn)行數(shù)值模擬研究.其中權(quán)重系數(shù)wi的參數(shù)設(shè)置如下：w0=4/9;當(dāng)i=1-4時(shí),wi=1/9;當(dāng)i=5-8時(shí),wi=1/36. c i為

3 模型驗(yàn)證

與其他研究者的模型驗(yàn)證方法一樣[21,33-36],本文用Laplace定律驗(yàn)證程序的正確性.初始時(shí),在N x×N y的格子區(qū)域中心放置一個(gè)半徑為r、密度為ρl的靜止液滴,其余區(qū)域充滿著蒸汽,密度為ρv.當(dāng)系統(tǒng)達(dá)到平衡狀態(tài)時(shí),液滴呈圓形.根據(jù)Laplace定律,液滴內(nèi)、外壓力差(Δp)與表面張力σ和液滴半徑r滿足關(guān)系式

本文模擬不同的液滴半徑下的系統(tǒng)狀態(tài)以驗(yàn)證Lap lace定律,半徑的變化為r=18,20,22,24, 26,28,30,32,35,40,45,50,并考慮了不同黏性比(M)下的系統(tǒng)狀態(tài),以確?？尚哦?其取值分別為M=3.78,6.19,7.57,9.10.表1列出了由Maxwell重構(gòu)得到的具體參數(shù)值.在以上所有的模擬中N x×N y=128×128,x方向和y方向均為周期邊界條件.圖1給出了不同黏性比M 下的液滴內(nèi)外壓力差(pi-po)和液滴半徑1/r的關(guān)系.從圖中可以得到數(shù)值模擬結(jié)果和Lap lace定律很符合.同時(shí),根據(jù)方程(11)可知,圖中直線的斜率即為對(duì)應(yīng)的表面張力,在本文中黏性比M=3.78,6.19,7.57,9.10下的表面張力分別為 0.0176,0.0445,0.0605和0.0782.

圖1 不同黏性比M下的(pi-po)和1/r的關(guān)系Fig.1.Relationship between(pi-po)and 1/r at different viscosity ratios M.

表1 不同黏性比下的參數(shù)值Tab le 1.Param eter values of different viscosity ratios.

4 模擬結(jié)果與分析

圖2給出了微通道結(jié)構(gòu)示意圖,灰色區(qū)域代表寬度為R的固體,由于其宏觀參數(shù)不發(fā)生變化因此不參與計(jì)算;白色區(qū)域?yàn)榱黧w流動(dòng)的區(qū)域,數(shù)值模擬時(shí)只需要模擬白色區(qū)域.模型包含一個(gè)圓心O1半徑為R1的凸起半圓(obstacle A),和一個(gè)圓心O2半徑為R2潤濕性為θ的半圓腔(cavity B),兩者之間的距離為D,凸起A代表壁面的粗糙度.初始時(shí)刻,密度為ρl的液體放在腔B中,其他計(jì)算區(qū)域?yàn)槊芏葹棣裿的氣體.氣體從左側(cè)進(jìn)口流入,右側(cè)出口流出.在模擬中,整個(gè)系統(tǒng)的網(wǎng)格為N×L,其中N=R+W是豎直方向上的網(wǎng)格數(shù),L為水平方向上的網(wǎng)格數(shù),在流體流動(dòng)的x方向上還施加有一個(gè)質(zhì)量力G.不失一般性,本文忽略了重力的影響,因?yàn)橹亓υ诖怪狈较蛏系淖饔煤苄?而且流體是在毛細(xì)力和黏性力的作用下流動(dòng)的[21].固體壁面上選用無滑移邊界條件,進(jìn)出口選用周期性邊界條件以保證質(zhì)量守恒.

在驅(qū)替過程中,流體流動(dòng)的無量綱參數(shù)設(shè)置如下：毛細(xì)數(shù)Ca=Uρvν/σ,其中U是在質(zhì)量力作用下流體的最大速度,表達(dá)式為U=W2G/(8ν).在模擬中,主要是通過兩個(gè)參數(shù)衡量驅(qū)替效果：驅(qū)替效率De和驅(qū)替時(shí)間t.其中De是被驅(qū)出腔B的液體質(zhì)量與原始腔B液體總質(zhì)量的比值;t是液體從初始時(shí)刻開始到被驅(qū)到微通道出口的時(shí)間,這里取W/U為特征時(shí)間,在文中出現(xiàn)的時(shí)間為無量綱時(shí)間.

由于不同的接觸角對(duì)應(yīng)不同的潤濕性,在這里計(jì)算了接觸角以滿足后續(xù)研究的需要.潤濕性有三種形式：親水性、疏水性、和中性.當(dāng)液滴在材料表面接觸角大于90°時(shí),材料表面具有疏水性;當(dāng)液滴在材料表面接觸角小于90°時(shí),材料表面具有親水性;當(dāng)液滴在材料表面接觸角等于90°時(shí),材料表面具有中性.在本文使用的模型中,通過改變壁面密度ρw來控制接觸角的大小[37,38].

我們模擬了一個(gè)液滴在壁面上的展開情況來得到壁面密度和接觸角的關(guān)系.初始時(shí),一個(gè)半徑為r密度為ρl的液滴放在一個(gè)矩形計(jì)算區(qū)域底部的中心,其他區(qū)域充滿著密度為ρv的氣體.在所有的模擬中,計(jì)算區(qū)域網(wǎng)格數(shù)為Nx×Ny=150×300.液體密度ρl=3.0404,氣體密度ρv=0.3342, g=-10.5.x方向?yàn)橹芷谛赃吔鐥l件,y方向?yàn)闊o滑移邊界條件.圖3給出了不同壁面密度(ρw)下分別對(duì)應(yīng)的接觸角(θ).從圖中可以看出當(dāng)壁面密度從ρl減小到ρv時(shí),對(duì)應(yīng)的接觸角從0°增加到180°.

圖2 微通道結(jié)構(gòu)示意圖Fig.2.Schem atic illustration of the channel structure.

圖3 (網(wǎng)刊彩色)不同壁面密度下的接觸角Fig.3.(color on line)Contact angle at d iff erent su rface densities.

4.1 潤濕性的影響

固體壁面的潤濕性是驅(qū)替過程的關(guān)鍵因素之一,本節(jié)研究了潤濕性對(duì)驅(qū)替過程的影響.眾所周知,固體壁面潤濕性可以量化為氣液固三相的接觸角,因此本節(jié)研究了不同的接觸角θ對(duì)驅(qū)替過程的影響.在研究中,分別考慮了以下接觸角θ=20°, 40°,50°,60°,75°,90°,105°,130°和160°.在數(shù)值模擬中,其他參數(shù)設(shè)置如下：N×L=100×1000, ρl=3.0404,ρv=0.3342,g=-10.5,ν=0.1667, R=W=R2=2R1=2D=50,Ca=0.05,圓心O2的坐標(biāo)為(200,50).

表2給出了驅(qū)替過程在不同接觸角下的驅(qū)替效率(D e)和驅(qū)替時(shí)間(t).從表中可以清晰地得到壁面的潤濕性對(duì)驅(qū)替過程有很重要的影響：當(dāng)接觸角小于某一角度時(shí)(這里是75°),驅(qū)替效率D e隨著θ的增加而增加,當(dāng)接觸角增加到某值后,De不再變化;當(dāng)接觸角很小時(shí),例如θ=20°,驅(qū)替效率D e為0.6552,在這種情況下腔B液體有大部分未被驅(qū)出,隨著接觸角的增加,驅(qū)替效率增加;當(dāng)θ=60°, De增加到0.8951,這說明腔B中液體大部分已經(jīng)被驅(qū)出;當(dāng)θ大于75°時(shí),D e等于1.0,即腔B液體全部被驅(qū)出.上述數(shù)值模擬的結(jié)果與理論預(yù)測(cè)的一致.事實(shí)上,固體壁面的潤濕性受黏附力(adhesive forces)和凝聚力(cohesive forces)的影響.接觸角較小對(duì)應(yīng)著較大的黏附力,液體不容易與壁面分離,因此較難從腔B中被驅(qū)出;隨著接觸角的增加,黏附力變小,部分液體容易被驅(qū)出腔B,因而De增加;當(dāng)凝聚力大于黏性力時(shí),腔B中的液體全部被驅(qū)出.從表2還可以看出,接觸角越大,驅(qū)替時(shí)間(t)越短.綜上所述,增加接觸角對(duì)驅(qū)替過程有利.

為了更直觀地揭示接觸角對(duì)驅(qū)替過程的影響,圖4給出了不同接觸角下的瞬時(shí)流動(dòng)情況.當(dāng)驅(qū)替過程在接觸角為20°和60°時(shí),初始時(shí)一部分液滴流出腔B,并鋪展在壁面上,這部分液體上受到來自進(jìn)口流體的拖拽力,使得液體在壁面上緩慢地向出口流動(dòng).同時(shí)在固體壁面上的液體和腔內(nèi)剩余液體之間的連接處,表面張力和拖拽力同時(shí)施加在此處液體上.隨著越來越多的液體流出,拖拽力會(huì)克服表面張力,使得液體在此處斷裂.隨后壁面上的液體流向出口,腔內(nèi)液體依然留在腔內(nèi).同時(shí)接觸角20°驅(qū)替過程的驅(qū)替時(shí)間大于接觸角為60°驅(qū)替過程的驅(qū)替時(shí)間.另外,接觸角大于75°時(shí)會(huì)略有不同,主要體現(xiàn)在液體不會(huì)發(fā)生斷裂,并全部被驅(qū)出微通道,這是由于接觸角增大,腔B中液體受到的黏附力變小,因此作用在液體上的驅(qū)力更容易將液體驅(qū)出.

表2 不同接觸角下的驅(qū)替效率和時(shí)間Tab le 2.D isp lacem ent efficiency and times at d iff erent contact angle.

圖4 (網(wǎng)刊彩色)不同接觸角下的流動(dòng)過程 (a)θ=20°,從左到右t=0,20,40,50;(b)θ=60°,從左到右t=0,20,34,44;(c)θ=90°,從左到右t=0,20,24,36Fig.4.(color on line)Flow process at d iff erent contact angle:(a)θ=20°,from left to right t=0,20,40, 50;(b)θ=60°,from left to right t=0,20,34,44;(c)θ=90°,from left to right t=0,20,24,36.

4.2 壁面粗糙度的影響

粗糙的壁面將影響流場(chǎng)的變化,因而本節(jié)研究壁面粗糙度對(duì)驅(qū)替過程的影響.在本節(jié)中,通過改變凸起A的尺寸(R1),研究不同粗糙度的壁面對(duì)驅(qū)替過程的影響.其中凸起A半徑的參數(shù)設(shè)置如下：R1=0.0R2,0.1R2,0.2R2,0.3R2,0.4R2,0.5R2, 0.6R2,0.7R2,0.8R2和0.9R2.其他參數(shù)和邊界條件與4.1節(jié)中設(shè)置相同.在上一節(jié)中的研究已經(jīng)說明固體壁面的潤濕性對(duì)驅(qū)替過程有重要的影響,因此為了模擬結(jié)果的一般性,在這里也考慮了不同接觸角下(θ=40°,50°,60°,75°和90°)的驅(qū)替過程.

不同尺寸的凸起A在不同接觸角下對(duì)驅(qū)替效率(D e)和驅(qū)替時(shí)間(t)的影響如圖5所示,其中對(duì)應(yīng)θ=40°和θ=50°的垂直向上的折線代表t=∞.從圖中可以看出凸起A對(duì)驅(qū)替過程有著不可忽略的影響,并且R1對(duì)驅(qū)替過程的影響依賴于壁面潤濕性(θ).當(dāng)接觸角(θ)小于90°并大于60°時(shí),D e隨著R1的增加而增加,當(dāng)De增加到1.0后便保持不變;然而當(dāng)接觸角等于或小于50°時(shí),D e首先隨著R1的增加而增加,當(dāng)R1大于某一值后De急劇減小.當(dāng)θ大于等于90°,由于效率均為1.0,曲線重合因而未給出.情況下所有的液體都可以從腔B中驅(qū)出,而且驅(qū)替時(shí)間(t)幾乎不隨著凸起A的尺寸(R1)改變而改變.這說明對(duì)于液體可以全部被驅(qū)出的情況,凸起A的尺寸R1對(duì)驅(qū)替過程的影響并不大,從這個(gè)意義上說,與粗糙度的大小相比,潤濕性對(duì)驅(qū)替過程影響更大.相應(yīng)地,t的變化也是一個(gè)復(fù)雜的過程.對(duì)于所有的接觸角來說,當(dāng)R1小于0.7R2時(shí),t隨著R1的增加而減小;當(dāng)R1大于0.7R2時(shí),t隨著R1的增加而增加,液體不能被驅(qū)出腔B,驅(qū)替時(shí)間無窮大.這種現(xiàn)象出現(xiàn)的原因可以從驅(qū)替過程的流場(chǎng)分析,如圖6所示,隨著R1的增加,腔B中x方向的速度減小,y方向的速度增加,這使得液體在腔B中更容易被驅(qū)出：驅(qū)替效率增加,驅(qū)替時(shí)間減小.但當(dāng)R1大于某一值后,通過凸起的速度會(huì)變得很小,因此會(huì)使驅(qū)替過程變難：驅(qū)替時(shí)間增加(如θ=60°和θ=75°的情況);或者液體不能從腔B中驅(qū)出,使得驅(qū)替時(shí)間無窮大(如θ=50°和θ=40°的情況).

圖5 不同尺寸的凸起A在不同接觸角下對(duì)驅(qū)替效率和驅(qū)替時(shí)間的影響Fig.5.Dependency of D e and t on the d iff erent size of obstacle A at different contact angle.

圖6 接觸角為50°時(shí)的流動(dòng)結(jié)果:從上到下:R1/R2= 0.1,0.6,0.8Fig.6.F low resu lts of d iff erent sizes of obstacle at contact angle 50°:from top to bottom:R1/R2=0.1, 0.6,0.8.

4.3 黏性比的影響

在本節(jié)中我們研究驅(qū)替過程中另一個(gè)重要的參數(shù)黏性比(M).研究中考慮了七組不同的黏性比,分別是M=3.78,4.23,4.94,5.68,6.19,7.57, 9.10(具體的參數(shù)如表1所列).Ca=0.1,其他參數(shù)值和邊界條件與上節(jié)相同.

圖7給出了不同黏性比(M)對(duì)驅(qū)替效率(D e)和驅(qū)替時(shí)間(t)的影響,圖中虛線代表t=∞.從圖中可以很清楚地看出：如圖中虛線所示,在小黏性比時(shí),液體不能被驅(qū)出微通道,驅(qū)替時(shí)間無窮大.隨著黏性比(M)的增加,驅(qū)替效率De增加,驅(qū)替時(shí)間t減小,即在黏性比較大時(shí),腔B液體驅(qū)出質(zhì)量較多并且驅(qū)替過程更快.為了使結(jié)果更直觀,圖8展示了不同黏性比(M)下的驅(qū)替結(jié)果.如圖所示,小黏性比的情況下,液體不能從腔B中被驅(qū)出;增加黏性比會(huì)使得腔B中液體被驅(qū)出,同時(shí)被驅(qū)出液體的質(zhì)量隨著黏性比的增加而增加.從液體鋪在壁面上的長度可以看出,黏性比越大,長度越長,說明大黏性比下的液體受到的驅(qū)力較大.這也同時(shí)解釋了黏性比(M)越大,驅(qū)替效率(D e)越高,驅(qū)替時(shí)間(t)越小.

圖7 不同黏性比(M)對(duì)驅(qū)替效率(D e)和驅(qū)替時(shí)間(t)的影響Fig.7.Dependency of D e and t on the d iff erent viscosity ratio(M).

圖8 (網(wǎng)刊彩色)不同黏性比(M)下的驅(qū)替結(jié)果Fig.8.(color on line)Disp lacem ent results at different density ratios(M).

4.4 距離的影響

本節(jié)研究凸起A和腔B間的距離(D)對(duì)驅(qū)替過程的影響. 在數(shù)值模擬中D分別取 0.0R2, 0.25R2,0.5R2,0.75R2,1.0R2,1.5R2,2.0R2,2.5R2, 3.0R2,3.5R2,4.0R2,4.5R2,5.0R2,6.0R2,7.0R2,其中R2=50.數(shù)值計(jì)算中用到的其他參數(shù)如下：腔B的圓心為(500,50),R1=0.2R2;Ca=0.08; ρl=3.0404,ρv=0.3342,g=-10.5;θ=60°.

圖9給出了在不同距離(D)下的驅(qū)替效率(D e)和驅(qū)替時(shí)間(t).如圖所示,當(dāng)D=0R2即凸起A和腔B的距離最小時(shí),驅(qū)替效率最大,此時(shí)驅(qū)替時(shí)間也最短.隨著距離D的增大,驅(qū)替效率D e減小,同時(shí)驅(qū)替時(shí)間t增加.當(dāng)D增加到0.75R2時(shí)驅(qū)替效率達(dá)到最小值,對(duì)應(yīng)的驅(qū)替時(shí)間t達(dá)到最大值.然后隨著D繼續(xù)增大,驅(qū)替效率增大繼而保持不變,對(duì)應(yīng)的驅(qū)替時(shí)間減小并保持不變.進(jìn)一步研究發(fā)現(xiàn),最后不再變化的驅(qū)替效率和驅(qū)替時(shí)間對(duì)應(yīng)沒有凸起結(jié)構(gòu)時(shí)的驅(qū)替效率和驅(qū)替時(shí)間.此模擬結(jié)果說明有凸起結(jié)構(gòu)的情況與沒有凸起結(jié)構(gòu)的情況相比,當(dāng)凸起結(jié)構(gòu)和腔B之間的距離較小時(shí)該凸起結(jié)構(gòu)對(duì)驅(qū)替結(jié)果是有利的,當(dāng)凸起結(jié)構(gòu)和腔B之間的距離很大時(shí),凸起結(jié)構(gòu)對(duì)驅(qū)替沒有影響,當(dāng)凸起結(jié)構(gòu)和腔B之間的距離在兩個(gè)極值之間時(shí),凸起結(jié)構(gòu)將對(duì)驅(qū)替有阻礙作用.

圖9 不同距離(D)對(duì)驅(qū)替效率(D e)和驅(qū)替時(shí)間(t)的影響Fig.9.Dependency of D e and t on the d iff erent distances(D).

5 結(jié) 論

本文采用改進(jìn)的偽勢(shì)格子Boltzmann方法研究復(fù)雜微通道內(nèi)的兩相流驅(qū)替問題,主要研究了壁面潤濕性、粗糙度、黏性比以及距離對(duì)驅(qū)替過程的影響.研究表明：1)隨著接觸角的增加,驅(qū)替效率增加,驅(qū)替時(shí)間減小;當(dāng)接觸角增加使得驅(qū)替效率等于1.0以后,驅(qū)替效率保持不變,驅(qū)替時(shí)間仍然隨著接觸角的增加而減小;2)壁面的粗糙度對(duì)驅(qū)替過程也有很重要的影響,本文用一個(gè)凸起的半圓凸起A代表壁面的粗糙性,其半徑的大小代表粗糙度的大小,數(shù)值模擬結(jié)果表明,隨著凸起半徑A的增加,驅(qū)替效率增加,驅(qū)替時(shí)間減小;然而當(dāng)其半徑增加到某一值后,腔B液體不能被驅(qū)出;3)黏性比同樣影響驅(qū)替過程,增加流體的黏性比會(huì)使得驅(qū)替效率增加,驅(qū)替時(shí)間減小;4)距離D為0時(shí)驅(qū)替效率最大,隨著D的增加,驅(qū)替效率減小然后增加最后趨于常數(shù).

[1]Teck lenburg J,Neuweiler I,Dentz M,Carrera J,Geiger S,Ab ram ow ski C,Silva O 2013 Adv.Water Res.62 475

[2]Zhu X,Sui P C,D jilali N 2008 J.Power Sources 181 101

[3]Yang D,K rasowska M,Priest C,Ralston J 2014 Phys. Chem.Chem.Phys.16 24473

[4]Islam S F,Sundara R V,W hitehouse S,Althaus T O, Palzer S,Hounslow M J,Salm an A D 2016 Chem.Eng. Res.Des.110 160

[5]Li W Z,Sun H M,Dong B 2013 Chin.J.Com putat. M ech.1 106(in Chinese)[李維仲,孫紅梅,董波2013計(jì)算力學(xué)學(xué)報(bào)1 106]

[6]Prim ku lov B K,Lin F,Xu Z 2016 Co lloids Surf.A: Physicochem.Eng.Aspects 497 336

[7]Kop lik J,Banavar J R,W illem sen J F 1988 Phys.Rev. Lett.60 1282

[8]Zhou G,Chen Z,Ge W,Li J 2010 Chem.Eng.Sci.65 3363

[9]Jam aloei B Y,K harrat R 2010 Transp.Porous M ed.81 1

[10]Pram anik S,M ishra M 2016 Phys.Rev.E 94 043106

[11]Yang K,Guo Z 2016 Com put.F luids 124 157

[12]K ang Q,Zhang D,Chen S 2002 Phys.F luids 14 3203

[13]Kang Q,Zhang D,Chen S 2005 J.Fluid M ech.545 41

[14]K ang Q,Zhang D,Chen S 2004 Adv.W ater Res.27 13

[15]Huang H,Huang J J,Lu X Y 2014 Com put.F luids 93 164

[16]Dong B,Yan Y Y,LiW,Song Y 2010 Com put.F luids 39 768

[17]Dong B,Yan Y Y,LiW Z 2011 Transp.Porous.M ed. 88 293

[18]LiW Z,Dong B,Song Y C 2012 J.Dalian Univ.Techno l.3 343(in Chinese)[李維仲,董波,宋永臣2012大連理工大學(xué)學(xué)報(bào)3 343]

[19]Li J,Song Y C,LiW Z 2009 J.Therm al Sci.Technol. 4 284(in Chinese)[李娟,宋永臣,李維仲2009熱科學(xué)與技術(shù)4 284]

[20]Peng B L,Xu W,W en R F,Lan Z,Bai T,M a X H 2015 J.Eng.Therm ophys.4 820(in Chinese)[彭本利,徐威,溫榮福,蘭忠,白濤,馬學(xué)虎2015工程熱物理學(xué)報(bào)4 820]

[21]Liang H,Chai Z,Shi B,Guo Z,Li Q 2015 In t.J.M od. Phys.C 26 1550074

[22]Gunstensen A K,Rothm an D H,Zaleski S,Zanetti G 1991 Phys.Rev.A 43 4320

[23]Shan X,Chen H 1993 Phys.Rev.E 47 1815

[24]Shan X,Chen H 1994 Phys.Rev.E 49 2941

[25]Sw ift M R,Osborn W R,Yeom ans JM 1995 Phys.Rev. Lett.75 830

[26]Sw ift M R,O rland ini E,Osborn W R,Yeom ans J M 1996 Phys.Rev.E 54 5041

[27]Luo L S 1998 Phys.Rev.Lett.81 1618

[28]He X,Chen S,Zhang R 1999 J.Com put.Phys.152 642

[29]Guo Z,Zhao T S 2005 Phys.Rev.E 71 026701

[30]Yu Z,Fan L S 2009 J.Com put.Phys.228 6456

[31]Guo Z,Zheng C,Shi B 2002 Phys.Rev.E 65 046308

[32]M artys N S,Chen H 1996 Phys.Rev.E 53 743

[33]Zhang R,He X,Chen S 2000 Com pu t.Phys.Comm un. 129 121

[34]Fakhari A,Rahim ian M H 2009 Comm un.Non linear Sci.Num er.Sim u lat.14 3046

[35]Zheng H W,Shu C,Chew Y T 2006 J.Com put.Phys. 218 353

[36]Hao L,Cheng P 2000 J.Power Sources 190 435

[37]Huang H,Lu X 2009 Phys.Fluids 21 092104

[38]Li Q,Luo K H,Kang Q J,Chen Q 2014 Phys.Rev.E 90 053301

(Received 10 January 2017;revised manuscript received 4 May 2017)

Lattice Boltzmann simulation of immiscible displacement in the complex micro-channel?

Zang Chen-Qiang Lou Qin?

(School of Energy and Power Engineering,University of Shanghai for Science and Technology,Shanghai 200093,China)

The imm iscible disp lacem ent process inmicro-channel,which w idely existes in daily life and industrialp roduction,is an im portant research sub ject.This sub ject is a typical contact line prob lem involving com p licated fl uid-fluid interactions and fl uid-solid interactions which have attracted the interest ofmany scholars.Although the imm iscib le disp lacement in micro-channels has been studied by some researches,the problem is still not fully understood because them echanism of the imm iscib le disp lacement is very com p lex.In order to further explain the physical mechanism of imm iscible disp lacement process in micro-channels,detailed numerical simu lations are carried out in a com p lex micro-channel containinga semicircular cavity and a sem icircular by bulge using an im p roved pseudo-potential lattice Boltzm ann method(LBM).Thismodel overcomes the drawback of the dependence of the fl uid properties on the grid size,which exists in the original pseudo-potential LBM.Initially,the cavity is fi lled with the liquid and the rest of the area is fi lled with its vapour.The sem icircular bulge rep resents the roughness of the micro-channel.The approach is fi rst validated by the Lap lace law.The results show that the numerical resu lts are in good agreement with the theoretical predictions.Then themodel is em p loyed to study the imm iscible disp lacem ent process in themicro-channel.The effects of the surface wettability,the surface roughness,the viscosity ratio between the liquid phase and the gas phase,and the distance between the sem icircular cavity and the sem icircu lar bu lge are studied.The simulation resu lts show that the in fluence of the surface wettability on the disp lacem ent process is a decisive factor com pared with other factors. W ith the increase of the contact angle,the disp lacem ent efficiency increases and the disp lacem ent time decreases.W hen the contact angle is larger than a certain value,all of the liquid can be disp laced from the cavity.At that time,the disp lacement efficiency is equal to 1.The above resu lts are consistent with the theoretical prediction that with the increase of the contact angle,the liquid is easily driven out of the cavity because the adhesion force of the liquid in the cavity decreases.On the other hand,the influence of the surface roughness on the disp lacem ent process ism ore com p lex. The disp lacem ent efficiency increases with the radius of the sem icircle bulge increasing in a certain range.W hen the radius is larger than a certain value,the liquid cannot be ejected from the cavity due to the velocity around the cavity is too sm all.Furtherm ore,the liquid cannot be disp laced from the cavity at a sm all viscosity ratio.As the viscosity ratio increases,the disp lacem ent efficiency increases and the disp lacem ent time decreases.As for the distance between the sem icircular bulge and the sem icircular cavity,it promotes the disp lacement process at an early stage.W hen the distance exceeds a certain value,it has little effect on the disp lacem ent process.

imm iscible disp lacement,lattice Boltzmann method,disp lacement efficiency,disp lacem ent time

PACS:47.11.-j,47.55.Ca,47.56.+r DO I:10.7498/aps.66.134701

?國家自然科學(xué)基金(批準(zhǔn)號(hào):51406120)資助的課題.

?通信作者.E-m ail:louqin560916@163.com

PACS:47.11.-j,47.55.Ca,47.56.+r DO I:10.7498/aps.66.134701

*Pro ject supported by the National Natural Science Foundation of China(G rant No.51406120).

?Corresponding author.E-m ail:louqin560916@163.com