復雜微通道內非混相驅替過程的格子Boltzm ann方法?

臧晨強 婁欽

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

復雜微通道內非混相驅替過程的格子Boltzm ann方法?

臧晨強 婁欽?

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

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

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

非混相驅替,格子Boltzmann方法,驅替效率,驅替時間

1 引 言

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

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

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

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

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

2 數值方法

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

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

其中,wi是權重系數.方程(1)中的Fi是力項,它的表達式為[31]

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

其中g是一個常數,表示流體間作用力強度.ψ(x)是“有效質量”,它是與局部密度有關的一個函數：

流固作用力Fads的表達式為[32]

其中ρw是壁面密度;s(x)是一個開關函數,當x是固體點時,s(x)取值為1,而當x為流體點時,s(x)取值為0.宏觀物理量的表達式為：

本文使用D2Q9模型進行數值模擬研究.其中權重系數wi的參數設置如下：w0=4/9;當i=1-4時,wi=1/9;當i=5-8時,wi=1/36. c i為

3 模型驗證

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

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

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

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

4 模擬結果與分析

圖2給出了微通道結構示意圖,灰色區域代表寬度為R的固體,由于其宏觀參數不發生變化因此不參與計算;白色區域為流體流動的區域,數值模擬時只需要模擬白色區域.模型包含一個圓心O1半徑為R1的凸起半圓(obstacle A),和一個圓心O2半徑為R2潤濕性為θ的半圓腔(cavity B),兩者之間的距離為D,凸起A代表壁面的粗糙度.初始時刻,密度為ρl的液體放在腔B中,其他計算區域為密度為ρv的氣體.氣體從左側進口流入,右側出口流出.在模擬中,整個系統的網格為N×L,其中N=R+W是豎直方向上的網格數,L為水平方向上的網格數,在流體流動的x方向上還施加有一個質量力G.不失一般性,本文忽略了重力的影響,因為重力在垂直方向上的作用很小,而且流體是在毛細力和黏性力的作用下流動的[21].固體壁面上選用無滑移邊界條件,進出口選用周期性邊界條件以保證質量守恒.

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

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

我們模擬了一個液滴在壁面上的展開情況來得到壁面密度和接觸角的關系.初始時,一個半徑為r密度為ρl的液滴放在一個矩形計算區域底部的中心,其他區域充滿著密度為ρv的氣體.在所有的模擬中,計算區域網格數為Nx×Ny=150×300.液體密度ρl=3.0404,氣體密度ρv=0.3342, g=-10.5.x方向為周期性邊界條件,y方向為無滑移邊界條件.圖3給出了不同壁面密度(ρw)下分別對應的接觸角(θ).從圖中可以看出當壁面密度從ρl減小到ρv時,對應的接觸角從0°增加到180°.

圖2 微通道結構示意圖Fig.2.Schem atic illustration of the channel structure.

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

4.1 潤濕性的影響

固體壁面的潤濕性是驅替過程的關鍵因素之一,本節研究了潤濕性對驅替過程的影響.眾所周知,固體壁面潤濕性可以量化為氣液固三相的接觸角,因此本節研究了不同的接觸角θ對驅替過程的影響.在研究中,分別考慮了以下接觸角θ=20°, 40°,50°,60°,75°,90°,105°,130°和160°.在數值模擬中,其他參數設置如下：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的坐標為(200,50).

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

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

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

圖4 (網刊彩色)不同接觸角下的流動過程 (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 壁面粗糙度的影響

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

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

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

圖6 接觸角為50°時的流動結果:從上到下: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 黏性比的影響

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

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

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

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

4.4 距離的影響

本節研究凸起A和腔B間的距離(D)對驅替過程的影響. 在數值模擬中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.數值計算中用到的其他參數如下：腔B的圓心為(500,50),R1=0.2R2;Ca=0.08; ρl=3.0404,ρv=0.3342,g=-10.5;θ=60°.

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

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

5 結 論

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

[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計算力學學報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大連理工大學學報3 343]

[19]Li J,Song Y C,LiW Z 2009 J.Therm al Sci.Technol. 4 284(in Chinese)[李娟,宋永臣,李維仲2009熱科學與技術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)[彭本利,徐威,溫榮福,蘭忠,白濤,馬學虎2015工程熱物理學報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

?國家自然科學基金(批準號: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