婁欽 黃一帆 李凌

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

2) (上海市動力工程多相流動與傳熱重點實驗室,上海 200093)

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

1 引 言

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

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

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

以上研究工作提出了一些非牛頓氣液兩相流模型,并用提出的模型研究了微通道內(nèi)非牛頓兩相流的一些運動機理.但現(xiàn)有模型存在一些不足： 由于自由能模型不滿足伽利略不變性[24,25],因此基于該模型提出的非牛頓兩相流模型[14]在數(shù)值模擬過程中會產(chǎ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理論提出的不可壓氣?液兩相流模型相對于其他兩相流模型具有物理機理更清晰,壓力的演化和密度的演化相獨立,處理復雜界面變化時數(shù)值穩(wěn)定性更好等優(yōu)點[25,28?30].鑒于此,本文基于冪率流體,在He等[27]提出的針對牛頓流體的不可壓多相流模型基礎(chǔ)上,提出了不可壓非牛頓冪律流體氣?液兩相流模型.此外目前對于非牛頓冪律流體兩相流問題的研究大都局限在簡單通道內(nèi),而對于多孔介質(zhì)內(nèi)牛頓氣體驅(qū)替非牛頓冪律流體液體的機理研究尚不充分.因此,本文又用發(fā)展的非牛頓冪律流體模型研究了多孔介質(zhì)內(nèi)牛頓氣體驅(qū)替非牛頓冪律流體液體的問題.

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

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

在HCZ模型中分別采用兩個分布函數(shù)fα和gα描述指標參數(shù)和速度/壓力的演化過程,其對應的分布函數(shù)fα和gα分別為

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

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

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

流體密度ρ(?) 和運動黏度ν(?) 可由指標參數(shù)?計算得到

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

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

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

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

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

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

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

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

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

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

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

3 模型驗證

3.1 拉普拉斯定律驗證

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

為了驗證Laplace定律,在數(shù)值模擬中分別考慮了r=20,22,24,26,28和30六種情況.為了保證數(shù)值模擬結(jié)果具有一般性,對于每一種情況都研究了n=0.7,n=1.0與n=1.3三種情形.圖1給出了在不同的冪律指數(shù)下得到的液滴內(nèi)外的壓力差Pi-Po與半徑的倒數(shù) 1/r的關(guān)系,可知,對于所有的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 潤濕性固壁面靜態(tài)接觸角驗證

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

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

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

圖2 不同初始靜態(tài)接觸角 θeq 時得到的穩(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)接觸角 θ 與指標參數(shù) ?wall 的線性關(guān)系Fig.3.Linear relationship between steady state contact angle θ and the order parameter of a solid wall ?wall.

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

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

本小節(jié)驗證在不同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ù)相與分散相的進口是速度進口邊界,連續(xù)相出口是對流邊界條件[30],管徑的固壁面都采用無滑移邊界條件.圖5給出了在不同的Ca數(shù)下,分離的牛頓流體液滴(黃色).圖6給出了在不同的Ca數(shù)下,分離的牛頓流體液滴尺寸,并與數(shù)值結(jié)果[41]進行對比,得到了一致的結(jié)果.

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

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

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

4.1 Ca數(shù)的影響

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

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

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

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

4.2 動力黏度 M 的影響

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

圖9 驅(qū)替完成時,不同冪律指數(shù)情況下得到的氣液兩相動力黏度示意圖 (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對冪律流體驅(qū)替效率的影響Fig.10.Effects of Ca and power?law exponent n on power?law fluid displacement efficiency.

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

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

圖11 不同的動力黏性比M下,被驅(qū)替液為剪切變稀、牛頓與剪切變稠流體時得到的指進形態(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.

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

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

4.3 潤濕性的影響

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

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

圖13 不同的?潤濕性角度 θ 下,被?驅(qū)替液為剪切變稀、牛頓與剪切變稠流體時得到的指進形態(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對驅(qū)替過程的影響,圖14給出了不同固體表面潤濕性以及冪律指數(shù)n下得到的驅(qū)替效率.從圖14可以看出無論剪切變稀、牛頓以及剪切變稠流體,都有固體表面接觸角越小,驅(qū)替效率De越低.具體而言,n=0.7時,固體表面接觸角θ從45°到135°時驅(qū)替效率De的取值為0.611到0.838,增加了37.15%;n=1.0時,固體表面接觸角θ從45°到135°時驅(qū)替效率De的取值為0.530到0.703,增加了32.64%;n=1.3時,固體表面接觸角θ從45°到135°時驅(qū)替效率De的取值為0.471到0.620,增加了31.63%.同時觀察到,固體表面潤濕性固定時,冪律指數(shù)n越大,驅(qū)替效率De越低.

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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