婁欽 黃一帆 李凌

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

2) (上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室,上海 200093)

基于不可壓格子玻爾茲曼氣?液兩相流模型,建立了一個(gè)新的非牛頓冪律流體氣?液兩相流模型,并采用該模型研究了多孔介質(zhì)內(nèi)牛頓氣體驅(qū)替非牛頓冪律流體液體的驅(qū)替問(wèn)題,主要探究了Ca數(shù)、動(dòng)力黏度比M、固體表面潤(rùn)濕性θ、多孔結(jié)構(gòu)幾何類型及冪律指數(shù)n對(duì)驅(qū)替過(guò)程的影響.研究發(fā)現(xiàn)： 不論被驅(qū)替液體為剪切變稀流體、牛頓流體還是剪切變稠流體都有隨著Ca數(shù)增加,驅(qū)替速度越快,指進(jìn)現(xiàn)象越明顯,驅(qū)替效率越低.然而對(duì)于不同的冪率指數(shù)n,驅(qū)替效率隨Ca數(shù)的增加而減小的速率不同： n越大驅(qū)替效率隨著Ca數(shù)增加而減小的速率越慢.另一方面,隨著黏性比M增加,驅(qū)替效率減小,且冪律指數(shù)n越小,M對(duì)驅(qū)替效率的影響越大.此外,固體表面接觸角θ對(duì)驅(qū)替過(guò)程的影響也和被驅(qū)替流體的冪率指數(shù)n相關(guān),雖然對(duì)于n ＞ 1和n ＜ 1的情況都有隨著多孔介質(zhì)固體表面接觸角θ增加,驅(qū)替過(guò)程受到的阻力越小,指進(jìn)越來(lái)越不明顯,驅(qū)替完成時(shí)間和驅(qū)替效率增加,然而當(dāng)n ＞ 1時(shí),隨著n的增加,接觸角對(duì)驅(qū)替過(guò)程的影響越來(lái)越小.還研究了孔隙率相同的情況下,孔隙幾何類型不同時(shí)的驅(qū)替過(guò)程,從數(shù)值結(jié)果可以發(fā)現(xiàn)與多孔結(jié)構(gòu)為圓形和方形障礙物相比,當(dāng)多孔結(jié)構(gòu)的幾何類型為三角形時(shí),指進(jìn)現(xiàn)象最明顯,驅(qū)替效率最低.

1 引 言

剪切應(yīng)力與剪切應(yīng)變率之間不呈線性關(guān)系的非牛頓流體廣泛存在于我們的日常生活中,如牛奶、蛋青、血液、建筑中的泥漿、果醬等.此外,非牛頓流體對(duì)應(yīng)的氣?液兩相流動(dòng)在廢水處理、石油開采、塑料泡沫加工以及環(huán)境保護(hù)等大量工程應(yīng)用和科學(xué)研究中有著重要的作用,例如,在石油開采過(guò)程中[1,2],原油作為一種典型的非牛頓流體,其與注入地下深處儲(chǔ)油層的水或二氧化碳構(gòu)成了非牛頓多相滲流問(wèn)題.隨著能源短缺日益嚴(yán)峻,開展多相非牛頓滲流相關(guān)問(wèn)題的研究具有重要的科學(xué)和現(xiàn)實(shí)意義.

近幾十年來(lái),眾多學(xué)者對(duì)微通道內(nèi)的非牛頓兩相流流動(dòng)機(jī)理進(jìn)行了一系列研究,且大多數(shù)研究都基于冪率流體.Xu等[3]通過(guò)實(shí)驗(yàn)研究了單一氣泡在剪切變稀冪律流體中的上升行為,發(fā)現(xiàn)剪切變稀冪律流體的流變性質(zhì)對(duì)氣泡上升的影響非常顯著,剪切變稀冪律流體的冪律指數(shù)與羧甲基纖維素(CMC)溶液的質(zhì)量百分?jǐn)?shù)有關(guān),其溶液質(zhì)量百分?jǐn)?shù)越高,對(duì)應(yīng)的冪律指數(shù)越低;在質(zhì)量百分?jǐn)?shù)較低的CMC溶液中,氣泡上升過(guò)程根據(jù)其形狀變化可分為三個(gè)階段： 第一階段,氣泡脫離噴嘴后,形狀基本為球形,并逐漸轉(zhuǎn)化成扁橢球形;第二階段,氣泡形狀變形不連續(xù);最后階段,氣泡變形連續(xù);隨著CMC溶液質(zhì)量百分?jǐn)?shù)的增加,過(guò)渡階段逐漸消失;在質(zhì)量百分?jǐn)?shù)最大的CMC溶液中,氣泡首先是平坦階段的連續(xù)形狀,然后是形狀基本不變以恒定速度沿直線上升.Du等[4]采用實(shí)驗(yàn)方法研究了在流體聚焦裝置中分散螺紋對(duì)剪切變稀冪律流體中液滴破裂行為的影響,發(fā)現(xiàn)剪切變稀冪律流體的流變特性對(duì)液滴的破裂影響較小.Fu等[5]通過(guò)實(shí)驗(yàn)研究了微通道內(nèi)氣泡在冪律流體中的融合過(guò)程.在研究中觀察到了兩種氣泡融合機(jī)理： 第一種是兩個(gè)氣泡融合成一個(gè)氣泡,第二種是一個(gè)氣泡先分裂成兩個(gè)氣泡,然后再融合成一個(gè)氣泡,該工作還分析了氣泡融合的概率、氣泡的位置隨時(shí)間變化的關(guān)系以及氣泡融合時(shí)周圍的流場(chǎng)變化.采用實(shí)驗(yàn)方法,Salehi等[6]研究了含有聚合物的牛頓流體對(duì)牛頓液滴和非牛頓冪律流體液滴形成過(guò)程中的表面差異的影響,發(fā)現(xiàn)在液滴增長(zhǎng)的初始階段,牛頓流體和非牛頓冪律流體的行為相似,然而在液滴頸縮過(guò)程中,隨著變形速率的變化,牛頓流體和冪律流體之間存在明顯的差異.

以上學(xué)者運(yùn)用實(shí)驗(yàn)方法對(duì)非牛頓流體兩相流進(jìn)行了研究.隨著計(jì)算機(jī)科技和數(shù)值方法的迅速發(fā)展,采用數(shù)值方法研究復(fù)雜系統(tǒng)流動(dòng)問(wèn)題被廣大學(xué)者認(rèn)可.近年來(lái),有大量學(xué)者運(yùn)用數(shù)值方法研究了非牛頓流體兩相流問(wèn)題.Sontti等[7]運(yùn)用計(jì)算流體動(dòng)力學(xué)(computational fluid dynamics,CFD)方法對(duì)T型微通道內(nèi)牛頓液滴在非牛頓冪律流體氣體中的形成過(guò)程進(jìn)行了研究,主要分析了冪律流體冪律指數(shù)、牛頓流體和冪律流體流量變化以及兩相表面張力對(duì)液滴形成尺寸的影響,發(fā)現(xiàn)了冪律流體冪律指數(shù)越大,產(chǎn)生的液滴尺寸越小;表面張力越大,產(chǎn)生的液滴尺寸越小;并發(fā)現(xiàn)了牛頓相和非牛頓冪律流體相的流量比對(duì)產(chǎn)生的液滴尺寸有顯著影響.相比于傳統(tǒng)CFD方法,格子Boltzmann方法(LBM)作為近幾十年發(fā)展并不斷完善的一種介觀數(shù)值算法,因其易于處理流體之間以及流體和固體壁面之間的相互作用、計(jì)算效率高以及不需要追蹤界面等優(yōu)點(diǎn),被廣泛用來(lái)研究多相流問(wèn)題.目前已有眾多學(xué)者用LBM研究了牛頓兩相流問(wèn)題[8?12]和非牛頓流體多相流問(wèn)題[13?23].謝馳宇等[13]基于兩相流自由能模型[14],構(gòu)建了非牛頓流體兩相流模型,并用構(gòu)建的模型研究了非牛頓流體驅(qū)替牛頓流體越過(guò)簡(jiǎn)單障礙物和多孔介質(zhì)時(shí)的流動(dòng)機(jī)理.Shi和Tang[15]基于相場(chǎng)多相流模型[16],分析了光滑微通道內(nèi)牛頓流體驅(qū)替非牛頓冪律流體的指進(jìn)現(xiàn)象,他們分析了被驅(qū)替相的冪律指數(shù)、固壁潤(rùn)濕性與Bo數(shù)對(duì)驅(qū)替指進(jìn)現(xiàn)象的影響;接下來(lái),Shi和Tang[17]又分析了在含多孔介質(zhì)的微通道內(nèi)牛頓液體驅(qū)替冪律流體液體的流動(dòng)機(jī)理,重點(diǎn)研究了冪律流體的冪律指數(shù)、毛細(xì)數(shù)、黏性比、障礙物表面潤(rùn)濕性以及多孔介質(zhì)的排列順序?qū)χ高M(jìn)現(xiàn)象的影響.Ba等[18]基于顏色模型[19,20],構(gòu)建了一個(gè)正規(guī)則格子Boltzmann (RLB)的顏色梯度模型,并用該模型模擬兩相冪律流體在兩個(gè)平板之間的流動(dòng)問(wèn)題以及牛頓液滴在冪律流體中的變形和破裂問(wèn)題.此外,閔琪等[21]基于Shan?Chen (S?C)模型[22,23]構(gòu)建了冪律流體兩相流格子Boltzmann模型,模擬了牛頓液滴以及冪律流體液滴在固壁面的鋪展情況.

以上研究工作提出了一些非牛頓氣液兩相流模型,并用提出的模型研究了微通道內(nèi)非牛頓兩相流的一些運(yùn)動(dòng)機(jī)理.但現(xiàn)有模型存在一些不足： 由于自由能模型不滿足伽利略不變性[24,25],因此基于該模型提出的非牛頓兩相流模型[14]在數(shù)值模擬過(guò)程中會(huì)產(chǎn)生一些非物理效應(yīng);基于顏色模型和基于相場(chǎng)模型改進(jìn)的非牛頓兩相流模型只能模擬密度比為1∶1的情況[15,18];而在S?C模型中,平衡態(tài)性質(zhì)和網(wǎng)格相關(guān)[26];此外,現(xiàn)有的非牛頓兩相流模型中流體壓力的演化與流體密度直接相關(guān),容易產(chǎn)生數(shù)值不穩(wěn)定.而He等[27]基于Enskong理論提出的不可壓氣?液兩相流模型相對(duì)于其他兩相流模型具有物理機(jī)理更清晰,壓力的演化和密度的演化相獨(dú)立,處理復(fù)雜界面變化時(shí)數(shù)值穩(wěn)定性更好等優(yōu)點(diǎn)[25,28?30].鑒于此,本文基于冪率流體,在He等[27]提出的針對(duì)牛頓流體的不可壓多相流模型基礎(chǔ)上,提出了不可壓非牛頓冪律流體氣?液兩相流模型.此外目前對(duì)于非牛頓冪律流體兩相流問(wèn)題的研究大都局限在簡(jiǎn)單通道內(nèi),而對(duì)于多孔介質(zhì)內(nèi)牛頓氣體驅(qū)替非牛頓冪律流體液體的機(jī)理研究尚不充分.因此,本文又用發(fā)展的非牛頓冪律流體模型研究了多孔介質(zhì)內(nèi)牛頓氣體驅(qū)替非牛頓冪律流體液體的問(wèn)題.

2 不可壓冪律流體氣?液兩相流格子Boltzmann模型

He等直接從不可壓多相流的Boltzmann方程出發(fā),提出了一種穩(wěn)定性好的適用于牛頓流體的不可壓多相流模型[27],簡(jiǎn)稱為HCZ[31?35]模型,下面我們將該模型推廣到冪律流體.

在HCZ模型中分別采用兩個(gè)分布函數(shù)fα和gα描述指標(biāo)參數(shù)和速度/壓力的演化過(guò)程,其對(duì)應(yīng)的分布函數(shù)fα和gα分別為

式中α=0,1,2,···,b-1,b為離散速度方向個(gè)數(shù);x和t分別表示粒子運(yùn)動(dòng)的位置和時(shí)間;δt代表時(shí)間步長(zhǎng);τ為松弛時(shí)間,其與運(yùn)動(dòng)黏度ν相關(guān);?,ρ,u分別代表指標(biāo)函數(shù)、流體密度和流體速度;κ為決定表面張力大小的參數(shù);ψ(ρ)=p-ρcs2,其中p為流體壓力,演化方程中ψ(?) 和狀態(tài)方程相關(guān);在方程(1)中函數(shù)Γα(u) 表達(dá)式為

其中ωα代表權(quán)重系數(shù).演化方程中(x,t) 和(x,t)是分布函數(shù)對(duì)應(yīng)的平衡態(tài),其形式如下：

宏觀量指標(biāo)參數(shù)?、壓力p以及速度u的統(tǒng)計(jì)由下面方程給出：

流體密度ρ(?) 和運(yùn)動(dòng)黏度ν(?) 可由指標(biāo)參數(shù)?計(jì)算得到

其中ρg和ρ1分別代表氣相流體和液相流體密度,?h和?l為指標(biāo)參數(shù)的最大值和最小值,可根據(jù)狀態(tài)方程由Maxwell重構(gòu)得到.

以上模型只適用于剪切應(yīng)力與剪切應(yīng)變率之間呈線性關(guān)系的牛頓流體的多相流動(dòng)問(wèn)題,為了將該模型推廣到非牛頓流體,本文用冪率流體模型來(lái)體現(xiàn)流體的非牛頓特性.對(duì)于冪率流體,動(dòng)力黏度η的表達(dá)式為

其中γ為剪切速率,η0為稠度系數(shù),

為應(yīng)變張量,n為非牛頓流體冪律指數(shù).當(dāng)n＜1時(shí),流體為剪切變稀流體,即其動(dòng)力黏度η隨著剪切速率γ的增大而減小;當(dāng)n＞ 1時(shí),流體為剪切變稠流體,其動(dòng)力黏度η隨著剪切速率γ的增大而增大;當(dāng)n=1時(shí),流體就是通常的牛頓流體,其動(dòng)力黏度η保持為一個(gè)定值,此時(shí)η=η0=ρv.因此,根據(jù)n的取值來(lái)區(qū)分流體是牛頓流體還是非牛頓冪律流體.

通過(guò)對(duì)分布函數(shù)gα進(jìn)行Chapman?Enskog分析可以得到應(yīng)變張量Sαβ與分布函數(shù)之間有如下關(guān)系：

從(7)式可以看出,在該冪律流體模型中,應(yīng)變張量Sαβ可直接用分布函數(shù)和局部點(diǎn)宏觀量信息計(jì)算得到,避免了差分計(jì)算,可提高數(shù)值穩(wěn)定性.通過(guò)Chapman?Enskog分析還可以得到方程(1)對(duì)應(yīng)的宏觀動(dòng)力學(xué)方程為

其中Π是黏性應(yīng)力張量,Π=ρν(?u+u?).運(yùn)動(dòng)黏度系數(shù)和松弛時(shí)間τ之間關(guān)系為ν=(τ-0.5)cs2δt,這里是與格子速度c=dx/dt相關(guān)的常數(shù).

本文使用二維九數(shù)(D2Q9)模型來(lái)進(jìn)行數(shù)值模擬研究,雖然本文是開展二維研究,但已有文獻(xiàn)[36,37]表明,盡管二維得到的結(jié)果和實(shí)際三維的結(jié)果定量上可能會(huì)有所不同,但是定性上趨勢(shì)一致.在D2Q9模型中權(quán)重系數(shù)ωα設(shè)置為： 當(dāng)α=0時(shí)當(dāng)α=5－8,為離散速度,其表達(dá)式為[38]

演化方程中梯度和拉普拉斯算子的離散方法均采用二階中心各向同性方法(ICS)[39].

演化方程中ψ(?) 由狀態(tài)方程決定,在本文中采用Carnahan?Starling狀態(tài)方程[27],其對(duì)應(yīng)的ψ(?) 可寫為如下形式：

其中a決定分子間相互吸引力強(qiáng)度,R為氣體常數(shù),T為流體溫度.

3 模型驗(yàn)證

3.1 拉普拉斯定律驗(yàn)證

本小節(jié)采用Laplace定律來(lái)驗(yàn)證模型的正確性.初始時(shí)在網(wǎng)格數(shù)為 128×128 區(qū)域中心內(nèi)放置半徑r,密度ρl=0.5 的靜止冪律流體圓形液滴,其余區(qū)域充滿著密度為ρg=0.1 的牛頓氣體.計(jì)算域四周的邊界條件均為周期性邊界條件.根據(jù)Laplace定律可知,當(dāng)系統(tǒng)達(dá)到穩(wěn)定時(shí)表面張力σ恒定,且液滴內(nèi)外的壓力差Pi-Po與半徑的倒數(shù)1/r呈線性關(guān)系,關(guān)系式如下：

為了驗(yàn)證Laplace定律,在數(shù)值模擬中分別考慮了r=20,22,24,26,28和30六種情況.為了保證數(shù)值模擬結(jié)果具有一般性,對(duì)于每一種情況都研究了n=0.7,n=1.0與n=1.3三種情形.圖1給出了在不同的冪律指數(shù)下得到的液滴內(nèi)外的壓力差Pi-Po與半徑的倒數(shù) 1/r的關(guān)系,可知,對(duì)于所有的n(流體分別為剪切變稀流體、牛頓流體和剪切變稠流體),模擬結(jié)果與Laplace定律均一致.

圖1 液滴內(nèi)外壓力差 Pi-Po 和半徑倒數(shù)1/r之間的關(guān)系Fig.1.Relationship between pressure jump across the droplet interface Pi-Po and inverse of droplet radius 1/r.

3.2 潤(rùn)濕性固壁面靜態(tài)接觸角驗(yàn)證

潤(rùn)濕現(xiàn)象在自然界中廣泛存在,如水滴在玻璃板上將迅速鋪展開,而水銀滴在玻璃板上將呈現(xiàn)為球滴狀,這就是潤(rùn)濕性不同程度的結(jié)果.潤(rùn)濕性反映流體和固體之間相互作用力的強(qiáng)度.在復(fù)雜微通道內(nèi),其是影響氣?液?固或液?液?固三相動(dòng)態(tài)行為的重要指標(biāo)參數(shù).在格子Boltzmann方法(LBM)中,通過(guò)潤(rùn)濕性邊界條件來(lái)描述壁面的潤(rùn)濕性.本文考慮復(fù)雜微通道內(nèi)固體表面潤(rùn)濕性對(duì)牛頓氣體驅(qū)替液相非牛頓冪律流體的驅(qū)替動(dòng)態(tài)影響,潤(rùn)濕性邊界條件采用Davies等[40]提出的方法,在該方法中壁面的潤(rùn)濕強(qiáng)度采用表面親和性αs來(lái)刻畫,并把表面親和性與指標(biāo)參數(shù)?聯(lián)系起來(lái),其關(guān)系式為

其中σ12為氣?液表面張力;σs1,σs2分別代表固?液表面張力與固?氣表面張力.

下面對(duì)冪律流體靜態(tài)液滴的靜態(tài)接觸角進(jìn)行驗(yàn)證.數(shù)值模擬中網(wǎng)格數(shù)為128 × 65,在計(jì)算區(qū)域下邊界中心處放置半徑r=20 ,密度ρl=0.5 以及冪律指數(shù)n=0.7的非牛頓冪律流體半圓液滴,液滴周圍充滿了ρg=0.1 的牛頓氣體,初始兩相的運(yùn)動(dòng)黏度均為νg=νl=1/6.計(jì)算域的邊界條件設(shè)置為： 上下壁面是無(wú)滑移邊界條件,左右是周期邊界條件.

圖2 不同初始靜態(tài)接觸角 θeq 時(shí)得到的穩(wěn)態(tài)接觸角θ(a) θeq=60o ;(b) θeq=90o;(c)θeq=120oFig.2.Steady state contact angles θ obtained with the dif?ferent values of static contact angles θeq ： (a) θeq=60o ;(b) θeq=90o;(c) θeq=120o.

圖3 穩(wěn)態(tài)接觸角 θ 與指標(biāo)參數(shù) ?wall 的線性關(guān)系Fig.3.Linear relationship between steady state contact angle θ and the order parameter of a solid wall ?wall.

圖2展示了當(dāng)數(shù)值模擬達(dá)到穩(wěn)定時(shí),壁面靜態(tài)接觸角θeq分別為 60°,90°與 120°時(shí),模擬所得到的穩(wěn)態(tài)接觸角θ,其結(jié)果分別為57.6°,87.3°與117.7°.模擬結(jié)果和理論值之間的相對(duì)誤差分別為4.0%,3.0%與1.9%.圖3展示了數(shù)值模擬得到的壁面穩(wěn)態(tài)接觸角θ與固壁面上的指標(biāo)參數(shù)?wall的關(guān)系,結(jié)果與(13)和(14)式符合得較好.

3.3 T型通道內(nèi)液滴分離

本小節(jié)驗(yàn)證在不同Ca數(shù)下,T型微通道內(nèi)形成液滴的大小.T型微管道結(jié)構(gòu)如圖4所示,其中W0=W1=30,L=520,Y1=75,Y2=120.初始分散相的子管道充滿牛頓液相流體,連續(xù)相主通道充滿冪律流體,其冪律指數(shù)n=0.4.邊界條件設(shè)置為： 連續(xù)相與分散相的進(jìn)口是速度進(jìn)口邊界,連續(xù)相出口是對(duì)流邊界條件[30],管徑的固壁面都采用無(wú)滑移邊界條件.圖5給出了在不同的Ca數(shù)下,分離的牛頓流體液滴(黃色).圖6給出了在不同的Ca數(shù)下,分離的牛頓流體液滴尺寸,并與數(shù)值結(jié)果[41]進(jìn)行對(duì)比,得到了一致的結(jié)果.

圖4 T型通道問(wèn)題物理模型Fig.4.Physical model for the case of T shape channel.

圖5 不同Ca數(shù)對(duì)應(yīng)的液滴形態(tài) (a) Ca=0.06370;(b) Ca=0.06835;(c) Ca=0.07300;(d) Ca=0.07750;(e) Ca=0.0820;(f) Ca=0.08650;(g) Ca=0.0910Fig.5.Droplet morphology obtained under various values of Ca： (a) Ca=0.06370;(b) Ca=0.06835;(c) Ca=0.07300;(d) Ca=0.07750;(e) Ca=0.0820;(f) Ca=0.08650;(g) Ca=0.0910.

圖6 在剪切變稀冪律流體中,不同的 Ca 數(shù)下形成液滴的無(wú)量綱直徑(其中D是形成的液滴的直徑,H是管徑的直徑)Fig.6.Droplet dimensionless diameters at different values of Ca in shear thinning power?law fluid.D is diameters of the droplet and H is width of the main channel.

4 多孔介質(zhì)內(nèi)牛頓流體驅(qū)替冪律流體的數(shù)值結(jié)果

本節(jié)研究多孔介質(zhì)內(nèi)冪律流體氣液兩相流驅(qū)替問(wèn)題,所得到的結(jié)論對(duì)石油開采、二氧化碳地質(zhì)埋存等相關(guān)問(wèn)題均有一定的指導(dǎo)意義,由于本文主要關(guān)注牛頓流體驅(qū)替非牛頓冪律流體的機(jī)理研究,如最近文獻(xiàn)[42?46]中的處理方法一樣,將實(shí)際地質(zhì)結(jié)構(gòu)簡(jiǎn)化為一種理想的多孔結(jié)構(gòu),其示意圖見圖7.計(jì)算域的網(wǎng)格數(shù)為M×N,多孔障礙物的尺寸分別為A與B,兩障礙物中心線之間的水平距離為X,垂直距離為Y,最靠近進(jìn)口的障礙物中心離進(jìn)口的水平距離為S1.初始時(shí)刻在x＜S處充滿密度ρg=0.1 ,運(yùn)動(dòng)黏度νg=1/6 ,動(dòng)力黏度ηg=0.01667的牛頓流體(藍(lán)色),在x＞S處充滿被驅(qū)替液相(黃色)非牛頓冪律流體,其密度為ρl=0.5.邊界條件設(shè)置為： 左邊入口處采用速度邊界條件,右邊出口處采用出口邊界條件[47],多孔表面與上下固壁面均采用無(wú)滑移邊界條件.

需要特別指出的是下文中需要用到如下兩個(gè)變量： 無(wú)量綱毛細(xì)數(shù)Ca,其定義為Ca=uηg/σ,其中u是驅(qū)替相的速度,ηg是驅(qū)替相動(dòng)力黏度,σ是氣液兩相表面張力;驅(qū)替效率De,其定義為當(dāng)驅(qū)替流體到達(dá)出口時(shí)通道內(nèi)剩余的被驅(qū)替相體積與初始被驅(qū)替相體積的比值.

4.1 Ca數(shù)的影響

本小節(jié)研究Ca數(shù)對(duì)不混溶冪律流體驅(qū)替過(guò)程的影響.在數(shù)值模型中M×N=240 × 600,A=30,B=20,S=6,S1=40,X=60,Y=60,孔隙率ξ=0.8507.初始動(dòng)力黏度比M=5.0,即被驅(qū)替液初始運(yùn)動(dòng)黏度νl=1/6 ,動(dòng)力黏度ηl=0.08333 ,σ=0.0056.固體表面均為中性潤(rùn)濕(θ=90o).

圖8給出了在不同Ca數(shù)下被驅(qū)替液為剪切變稀、牛頓以及剪切變稠三種流體時(shí)得到的驅(qū)替完成時(shí)指進(jìn)形態(tài)圖.從圖8可以看出,不論被驅(qū)替液是剪切變稀(圖8(a)－(c)n=0.7)、牛頓(圖8(d)－(f)n=1.0)還是剪切變稠(圖8(g)－(i)n=1.3)流體,都有Ca數(shù)越大,指進(jìn)現(xiàn)象越明顯[17],驅(qū)替完成時(shí)花費(fèi)的時(shí)間越少.具體而言,當(dāng)n=0.7時(shí)(圖8(a)－(c)),Ca=0.0298對(duì)應(yīng)的驅(qū)替完成時(shí)間t=36.2,當(dāng)Ca數(shù)增加到0.0877時(shí),驅(qū)替時(shí)間減少到了t=11.1,減少了69.3%;而當(dāng)n=1.0時(shí)(圖8(d)－(f))和n=1.3 (圖8(g)－(i))時(shí),Ca數(shù)從0.0298增加到0.0877,驅(qū)替時(shí)間分別減少了68.8%和67.7%.發(fā)生以上現(xiàn)象的原因是由于Ca數(shù)是表征黏性力與表面張力比值的無(wú)量綱參數(shù),Ca數(shù)越大說(shuō)明表面黏性力越大而表面張力越小,而隨著黏性力的增加,驅(qū)替過(guò)程受到的阻力增加,隨著表面張力減小,氣液界面更容易發(fā)生變形,因此Ca數(shù)越大指進(jìn)現(xiàn)象越明顯.Shiri等[48]以及Liu等[49]在研究多孔介質(zhì)內(nèi)的驅(qū)替問(wèn)題時(shí)也得到了相似的結(jié)論.另一方面,從圖8還可以發(fā)現(xiàn)當(dāng)Ca數(shù)相同時(shí),當(dāng)被驅(qū)替相的冪率指數(shù)n越大時(shí),指進(jìn)現(xiàn)象越明顯,驅(qū)替完成所需時(shí)間越短[17].例如當(dāng)Ca=0.0298時(shí)(圖8(a)、圖8(d)、圖8(g)),n=0.7,1.0,1.3對(duì)應(yīng)的驅(qū)替完成時(shí)間分別為36.2,32.4,29.1;Ca=0.0595時(shí)(圖8(b)、圖8(e)、圖8(h)),n=0.7,1.0,1.3對(duì)應(yīng)的驅(qū)替完成時(shí)間分別為16.9,15.3,14.2;而Ca=0.0877時(shí)(圖8(c)、圖8(f)、圖8(i)),n=0.7,1.0,1.3對(duì)應(yīng)的驅(qū)替完成時(shí)間分別為11.1,10.1,9.4.因此對(duì)應(yīng)這三種Ca數(shù)的情況,隨著冪率指數(shù)n的增加驅(qū)替時(shí)間分別減小了19.6%,16.0%和15.3%.也就是說(shuō),冪率指數(shù)n越大,Ca數(shù)的增加導(dǎo)致的驅(qū)替時(shí)間減少的速率越來(lái)越慢.導(dǎo)致該現(xiàn)象的原因是對(duì)于剪切變稠流體隨著驅(qū)替過(guò)程的進(jìn)行其動(dòng)力黏度會(huì)大于初始時(shí)刻的動(dòng)力黏度,剪切變稀流體的動(dòng)力黏度會(huì)小于初始時(shí)刻的動(dòng)力黏度,而牛頓流體的動(dòng)力黏度會(huì)保持不變.為了說(shuō)明這一現(xiàn)象,圖9給出了Ca=0.0298時(shí)對(duì)應(yīng)驅(qū)替完成時(shí)刻的兩相流體的動(dòng)力黏度分布.從圖9可以看出,對(duì)于研究的所有n的值,當(dāng)驅(qū)替過(guò)程完成時(shí),驅(qū)替氣相的動(dòng)力黏度一直保持在初始值0.01667附近,而對(duì)應(yīng)冪律指數(shù)n=0.7 (圖9(a))、1.0 (圖9(b))以及1.3 (圖9(c))的被驅(qū)替液得到的動(dòng)力黏度分別為0.04409,0.08333與0.1579.即剪切變稀、牛頓以及剪切變稠流體被驅(qū)替過(guò)程中,其分別對(duì)應(yīng)的兩相動(dòng)力黏度比M小于5、等于5以及大于5.而兩相動(dòng)力黏度比越大,兩相流體間黏性力影響變大,即所受黏性阻力越大,因此指進(jìn)現(xiàn)象越明顯,驅(qū)替越快[15,17].

圖7 多孔介質(zhì)驅(qū)替模型Fig.7.The model for porous media displacement problem.

圖8 不同的 Ca 數(shù)下,被驅(qū)替液為剪切變稀、牛頓與剪切變稠流體時(shí)得到的指進(jìn)形態(tài)圖 (a)?(c) n=0.7;(d)?(f) n=1.0;(g)?(i) n=1.3Fig.8.Final finger patterns obtained under different values of Ca for shear thinning,Newtonian and shear thickening fluids： (a)?(c) n=0.7;(d)?(f) n=1.0;(g)?(i) n=1.3.

為了定量分析Ca數(shù)以及冪律指數(shù)n對(duì)驅(qū)替過(guò)程的影響,圖10給出了不同Ca數(shù)以及冪律指數(shù)n下得到的驅(qū)替效率.從圖10可以看出,不論被驅(qū)替液是剪切變稀、牛頓還是剪切變稠流體,都有Ca數(shù)越大,驅(qū)替效率De越低.具體而言,n=0.7時(shí),Ca數(shù)從0.0298到0.0877時(shí)驅(qū)替效率De的值從0.744減小到0.688,減小了7.52%;n=1.0時(shí),Ca數(shù)從0.0298到0.0877時(shí)驅(qū)替效率De的值從0.663減小到0.617,減小了6.93%;n=1.3時(shí),Ca數(shù)從0.0298到0.0877時(shí)驅(qū)替效率De的值從0.590減小到0.550,減小了6.78%.從以上分析以及圖10中曲線變化可以發(fā)現(xiàn),不論被驅(qū)替相是牛頓流體還是冪律流體,驅(qū)替效率隨著Ca數(shù)的增加而減小,然而驅(qū)替效率減小的速率隨著n的增加而減小.另外,Ca數(shù)一定時(shí),冪律指數(shù)n越大,驅(qū)替效率De越低.

4.2 動(dòng)力黏度 M 的影響

本小節(jié)研究動(dòng)力黏度比M對(duì)不混溶冪律流體驅(qū)替過(guò)程的影響,這里黏性比的增加是通過(guò)改變被驅(qū)替液的黏性實(shí)現(xiàn)的.在本小節(jié)Ca數(shù)設(shè)置為0.0446,其余參數(shù)與4.1節(jié)相同.

圖9 驅(qū)替完成時(shí),不同冪律指數(shù)情況下得到的氣液兩相動(dòng)力黏度示意圖 (a)n=0.7 ;(b)n=1.0 ;(c)n=1.3Fig.9.Schematic diagram of gas?liquid two phase dynamics viscosity obtained under different values of power?law exponent：(a)n=0.7;(b)n=1.0 ;(c)n=1.3.

圖10 Ca 數(shù)和冪律指數(shù)n對(duì)冪律流體驅(qū)替效率的影響Fig.10.Effects of Ca and power?law exponent n on power?law fluid displacement efficiency.

圖11給出了被驅(qū)替液為剪切變稀、牛頓以及剪切變稠三種流體在不同初始動(dòng)力黏度比M的情況下,驅(qū)替完成時(shí)指進(jìn)形態(tài)圖.從圖11可以發(fā)現(xiàn),對(duì)于被驅(qū)替液是剪切變稀流體的情況(圖11(a)－(c)n=0.7),黏度比從2.5增加到12.5,驅(qū)替時(shí)間從26.5減少到19.7,減少了25.7%.對(duì)于被驅(qū)替液為牛頓流體(圖11(d)－(f)n=1.0)和剪切變稠(圖11(g)－(i)n=1.3)流體的情況,動(dòng)力黏度比從2.5增加到12.5時(shí)對(duì)應(yīng)的驅(qū)替時(shí)間分別減少了19.6%和9.3%.因此對(duì)于所有n的情況都有初始動(dòng)力黏度比M越大,指進(jìn)現(xiàn)象越明顯,驅(qū)替完成時(shí)所花費(fèi)的時(shí)間越少[15,17].這是因?yàn)轲ば员仍酱笳f(shuō)明被驅(qū)替液黏性越大,而一個(gè)輕流體(黏性較小)驅(qū)替一個(gè)重流體(黏性較大)是非常困難的,因此指進(jìn)現(xiàn)象會(huì)更明顯.其他學(xué)者[48?51]也發(fā)現(xiàn)了類似的現(xiàn)象.從以上數(shù)據(jù)以及指進(jìn)形態(tài)圖還可以發(fā)現(xiàn)驅(qū)替相的冪律指數(shù)n越大,黏性比的增加對(duì)驅(qū)替過(guò)程的影響越小.另一方面從圖11可以觀察到當(dāng)初始動(dòng)力黏度比M相同時(shí),被驅(qū)替相的冪律指數(shù)n越大,指進(jìn)現(xiàn)象越明顯,對(duì)應(yīng)的驅(qū)替完成所花費(fèi)的時(shí)間越少.這與4.1節(jié)流體越黏稠,流體越難被驅(qū)替結(jié)論一致[15,17].如M=2.5 (圖11(a)、圖11(d)、圖11(g))時(shí),對(duì)應(yīng)n=0.7,1.0和1.3的驅(qū)替完成時(shí)間分別為26.5,23.5和20.3,此時(shí)隨著冪率指數(shù)n的增加,驅(qū)替完成時(shí)間減小了23.4%.

M=7.5 (圖11(b)、圖11(e)、圖11(h))時(shí),對(duì)應(yīng)以上三個(gè)冪率指數(shù)的驅(qū)替完成時(shí)間分別為21.4,19.6和18.5.M=12.5 (圖11(c)、圖11(f)、圖11(i))時(shí)對(duì)應(yīng)的驅(qū)替完成時(shí)間分別為19.7,18.9和18.4.因此對(duì)應(yīng)這兩種黏性比的情況,隨著冪率指數(shù)增加驅(qū)替時(shí)間分別減小了13.6%和6.6%.從以上數(shù)據(jù)分析可以發(fā)現(xiàn)隨著黏性比的增加,驅(qū)替速率減少,且當(dāng)黏性比較小時(shí),冪律指數(shù)n越大,流體越難被驅(qū)替,而隨著黏性比的增加,被驅(qū)替相是牛頓流體和冪律流體的驅(qū)替結(jié)果之間的差異越來(lái)越小.

圖11 不同的動(dòng)力黏性比M下,被驅(qū)替液為剪切變稀、牛頓與剪切變稠流體時(shí)得到的指進(jìn)形態(tài)圖 (a)?(c) n=0.7;(d)?(f) n=1.0;(g)?(i) n=1.3Fig.11.Final finger patterns obtained under different values of viscosity ratios M for shear thinning,Newtonian and shear thicken?ing fluids： (a)?(c) n=0.7;(d)?(f) n=1.0;(g)?(i) n=1.3.

為了進(jìn)一步分析初始動(dòng)力黏度比M以及冪律指數(shù)n對(duì)驅(qū)替過(guò)程的影響,圖12給出了不同初始動(dòng)力黏度比M以及冪律指數(shù)n下得到的驅(qū)替效率.從圖12可以看出,當(dāng)不論被驅(qū)替液是剪切變稀、牛頓還是剪切變稠流體,都有動(dòng)力黏度比M越大,驅(qū)替效率De越低,而且隨著n的增加得到的驅(qū)替效率曲線的斜率越來(lái)越平緩,說(shuō)明黏性比M對(duì)于驅(qū)替效率的影響隨著n的增加而減小.具體而言,n=0.7時(shí),M從2.5到12.5時(shí)驅(qū)替效率De的值從0.827減小到0.596,減小了27.93%;n=1.0時(shí),M從2.5到12.5時(shí)驅(qū)替效率De的值從0.730減小到0.562,減小了23.01%;n=1.3時(shí),M從2.5到12.5時(shí)驅(qū)替效率De的值從0.626減小到0.535,減小了14.54%.因此,不論被驅(qū)替相是牛頓流體還是冪律流體,驅(qū)替效率隨著M的增加而減小,然而驅(qū)替效率減小的速率隨著n的增加而減小.另外,M一定時(shí),冪律指數(shù)n越大,驅(qū)替效率De越低.

圖12 動(dòng)力黏度比M和冪律指數(shù)n對(duì)冪律流體驅(qū)替效率的影響Fig.12.Effects of viscosity ratio M and power?law expo?nent n on power?law fluid displacement efficiency.

4.3 潤(rùn)濕性的影響

本小節(jié)研究固體表面潤(rùn)濕性對(duì)冪律流體驅(qū)替效率的影響.在數(shù)值模擬中動(dòng)力黏度比M=5.0 ,毛細(xì)數(shù)Ca=0.0446 ,其余參數(shù)與4.2節(jié)相同.

圖13給出了被驅(qū)替液為剪切變稀、牛頓以及剪切變稠三種流體在四種不同潤(rùn)濕性(θ=45o,60o,120o與135°)情況下,驅(qū)替完成時(shí)得到的指進(jìn)形態(tài)圖.從圖13可以發(fā)現(xiàn),對(duì)于所有的n都有隨著固體表面接觸角θ增加,指進(jìn)現(xiàn)象越不明顯.這是因?yàn)楫?dāng)固體面的接觸角越大,固體表面的疏水性越強(qiáng),則被驅(qū)替液受到固壁的黏附力越小,所以被驅(qū)替液更容易驅(qū)替,則指進(jìn)現(xiàn)象越不明顯.同時(shí)由于接觸角越大被驅(qū)替流體越容易被驅(qū)替,表現(xiàn)為指進(jìn)長(zhǎng)度變短,所以被驅(qū)替液到出口的時(shí)間,也即驅(qū)替完成的時(shí)間增加.因此驅(qū)替流體到達(dá)出口的時(shí)間越長(zhǎng).其他學(xué)者[49,52]也發(fā)現(xiàn)了類似的現(xiàn)象.具體而言,對(duì)于所研究的四種接觸角的情況(θ=45°,60°,120°與135°),當(dāng)n=0.7 (圖13(a)、圖13(d)、圖13(g)、圖13(j))對(duì)應(yīng)的驅(qū)替完成時(shí)間分別為23.1,23.5,25.5和26.8;n=1.0 (圖13(b)、圖13(e)、圖13(h)、圖13(k))對(duì)應(yīng)的驅(qū)替完成時(shí)間分別為20.7,20.9,22.2,22.6;n=1.3 (圖13(c)、圖13(f)、圖13(i)、圖13(l))對(duì)應(yīng)的驅(qū)替時(shí)間分別為19.0,19.1,19.9和20.3.通過(guò)計(jì)算可以得到接觸角θ從45o增加到 135o時(shí),n=0.7,1.0和1.3對(duì)應(yīng)的驅(qū)替時(shí)間分別增加了16.0%,9.2%和6.8%,說(shuō)明隨著冪率指數(shù)n的增加,固體表面潤(rùn)濕性對(duì)驅(qū)替完成時(shí)間段的影響逐漸減小.另一方面,從圖13可以發(fā)現(xiàn)當(dāng)固體表面潤(rùn)濕性θ相同時(shí),被驅(qū)替相的冪律指數(shù)n越大,指進(jìn)越明顯,又一次說(shuō)明了被驅(qū)替相流體越黏稠越不容易被驅(qū)替.

圖13 不同的?潤(rùn)濕性角度 θ 下,被?驅(qū)替液為剪切變稀、牛頓與剪切變稠流體時(shí)得到的指進(jìn)形態(tài)圖 (a)?(c) θ=45° ;(d)?(f) θ=60° ;(g)(i) θ=120° ;(j)(l)θ=135°Fig.13.Final finger patterns obtained under different values of contact angles θ for shear thinning,Newtonian and shear thicken?ing fluids： (a)?(c) θ=45° ;(d)-(f) θ=60° ;(g)?(i) θ=120° ;(j)?(l) θ=135°.

為了定量分析固體表面接觸角θ以及冪律指數(shù)n對(duì)驅(qū)替過(guò)程的影響,圖14給出了不同固體表面潤(rùn)濕性以及冪律指數(shù)n下得到的驅(qū)替效率.從圖14可以看出無(wú)論剪切變稀、牛頓以及剪切變稠流體,都有固體表面接觸角越小,驅(qū)替效率De越低.具體而言,n=0.7時(shí),固體表面接觸角θ從45°到135°時(shí)驅(qū)替效率De的取值為0.611到0.838,增加了37.15%;n=1.0時(shí),固體表面接觸角θ從45°到135°時(shí)驅(qū)替效率De的取值為0.530到0.703,增加了32.64%;n=1.3時(shí),固體表面接觸角θ從45°到135°時(shí)驅(qū)替效率De的取值為0.471到0.620,增加了31.63%.同時(shí)觀察到,固體表面潤(rùn)濕性固定時(shí),冪律指數(shù)n越大,驅(qū)替效率De越低.

圖14 潤(rùn)濕性 θ 和冪律指數(shù)n對(duì)冪律流體驅(qū)替效率的影響Fig.14.Effects of contact angles θ and power?law expo?nent n on power?law fluid displacement efficiency.

4.4 多孔介質(zhì)幾何類型的影響

本小節(jié)研究多孔介質(zhì)幾何類型對(duì)不混溶冪律流體驅(qū)替過(guò)程的影響.在數(shù)值模擬中M×N=240 × 600,S=6,S1=40,X=60,Y=60,初始動(dòng)力黏度比M為5.0,Ca=0.0446,固體表面都是中性潤(rùn)濕(θ=90o).障礙物的幾何類型分別為正方形、圓形與等邊三角形三種情況,當(dāng)障礙物的幾何類型為正方形時(shí),A=B=30;當(dāng)障礙物的幾何類型為圓形時(shí),障礙物都是半徑為16.93的圓形;當(dāng)障礙物為等邊三角形時(shí),其邊長(zhǎng)為49.96.保證對(duì)于所研究的所有的情況孔隙率都為ξ=0.78125.

圖15給出了多孔介質(zhì)幾何類型和被驅(qū)替液為剪切變稀、牛頓以及剪切變稠三種流體驅(qū)替完成時(shí)指進(jìn)形態(tài)圖.從圖15可以發(fā)現(xiàn),對(duì)于所有的冪率指數(shù)n的情況,都有當(dāng)多孔介質(zhì)類型為三角形時(shí)指進(jìn)現(xiàn)象最明顯,因此驅(qū)替完成時(shí)所需的時(shí)間最短.當(dāng)n=0.4時(shí)(圖15(a)－(c)),對(duì)應(yīng)多孔介質(zhì)結(jié)構(gòu)為方形、圓形和三角形的驅(qū)替完成時(shí)間分別為25.0,24.5和20.9;當(dāng)n=0.7時(shí)(圖15(d)－(f)),對(duì)應(yīng)上述多孔介質(zhì)結(jié)構(gòu)的驅(qū)替完成時(shí)間分別為22.1,22.0和19.0;n=1.0時(shí)(圖15(g)－(i)),對(duì)應(yīng)上述多孔介質(zhì)結(jié)構(gòu)的驅(qū)替完成時(shí)間分別為20.7,20.3和17.7;n=1.3時(shí)(圖15(j)－(l)),對(duì)應(yīng)上述多孔介質(zhì)結(jié)構(gòu)的驅(qū)替完成時(shí)間分別為19.8,19.1和17.1;n=1.6時(shí)(圖15(m)－(o)),對(duì)應(yīng)上述多孔介質(zhì)結(jié)構(gòu)的驅(qū)替完成時(shí)間分別為19.3,18.5和16.9.從上述數(shù)據(jù)可以發(fā)現(xiàn)相對(duì)于其他兩種情況,多孔介質(zhì)結(jié)構(gòu)為三角形時(shí)驅(qū)替完成時(shí)間明顯減少,這是因?yàn)樵谙嗤目紫堵是闆r下,相對(duì)于其他兩種多孔結(jié)構(gòu),流體通過(guò)三角形結(jié)構(gòu)的多孔介質(zhì)結(jié)構(gòu)時(shí)所對(duì)應(yīng)的通道最小,因此驅(qū)替過(guò)程受到的阻力最大.從圖15還可以發(fā)現(xiàn)對(duì)于相同的多孔介質(zhì)結(jié)構(gòu)被驅(qū)替相的冪律指數(shù)n越大,指進(jìn)現(xiàn)象越明顯,這是因?yàn)楸或?qū)替相的冪率指數(shù)n越大,驅(qū)替過(guò)程的阻力就越大,因此指進(jìn)越明顯,對(duì)應(yīng)的驅(qū)替完成時(shí)間越短.

為進(jìn)一步定性分析障礙物幾何類型以及冪律指數(shù)n對(duì)驅(qū)替過(guò)程的影響,圖16給出了在多孔介質(zhì)幾何類型不一樣時(shí)不同冪律指數(shù)n下得到的驅(qū)替效率.從圖16可以看出對(duì)于所研究的多孔介質(zhì)結(jié)構(gòu),驅(qū)替效率都隨著n的增加而減小,該結(jié)論與圖15的指進(jìn)現(xiàn)象得到對(duì)應(yīng)：n越大指進(jìn)現(xiàn)象越明顯.這是因?yàn)殡S著冪率指數(shù)的增加,被驅(qū)替流體更黏稠,不利于驅(qū)替過(guò)程的進(jìn)行.從圖16還可以看出,對(duì)于所有的冪率指數(shù)n的情況,圓形多孔結(jié)構(gòu)和方形多孔結(jié)構(gòu)得到的驅(qū)替效率比較接近(與圖15中的驅(qū)替時(shí)間對(duì)應(yīng)),而三角形多孔介質(zhì)結(jié)構(gòu)得到的驅(qū)替效率明顯減小.

圖16 障礙物幾何類型和冪律指數(shù)n對(duì)冪律流體驅(qū)替效率的影響Fig.16.Effects of geometric type and power?law exponent n on power?law fluid displacement efficiency.

5 結(jié) 論

基于HCZ兩相流模型,建立了一個(gè)不可壓非牛頓冪律流體兩相流模型,并用該模型研究了在含多孔介質(zhì)的微通道內(nèi)牛頓氣體驅(qū)替非牛頓冪律流體液體的流動(dòng)機(jī)理,主要分析了非牛頓流體的冪律指數(shù)n、Ca數(shù)、初始動(dòng)力黏度比M、多孔介質(zhì)表面潤(rùn)濕性θ以及多孔介質(zhì)障礙物幾何類型對(duì)驅(qū)替過(guò)程的影響,得出以下結(jié)論：

1) 冪律指數(shù)n越大,驅(qū)替過(guò)程的指進(jìn)現(xiàn)象越明顯,驅(qū)替效率越低,且冪率指數(shù)n越大,其驅(qū)替過(guò)程受Ca數(shù)的影響,黏性比的影響越小;

2)對(duì)于剪切變稀,牛頓與剪切變稠流體,都有隨Ca數(shù)和黏性比的增加指進(jìn)現(xiàn)象越明顯,驅(qū)替效率越低;

3)固體表面接觸角越小,其表面殘留的液體就越多,驅(qū)替效率越低;

4)在孔隙率ξ相同的情況下,相對(duì)于多孔介質(zhì)內(nèi)障礙物為圓形和方形的情況,障礙物為等邊三角形時(shí),指進(jìn)現(xiàn)象最明顯,驅(qū)替效率最低.