劉高潔 邵子宇 婁欽

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

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

采用格子Boltzmann 方法研究了孔隙尺度下多孔介質(zhì)內(nèi)含流固溶解反應(yīng)的互溶驅(qū)替過程,重點(diǎn)研究了被驅(qū)替流體與驅(qū)替流體黏性差異較大的情況下,溶解反應(yīng)引起的多孔介質(zhì)內(nèi)部結(jié)構(gòu)變化對(duì)驅(qū)替過程的影響;定量分析了不同達(dá)姆科勒數(shù)及佩克萊數(shù)下多孔介質(zhì)孔隙率和驅(qū)替過程驅(qū)替效率隨時(shí)間的演變.研究結(jié)果表明:達(dá)姆科勒數(shù)較大時(shí),溶解反應(yīng)的發(fā)生會(huì)在多孔介質(zhì)內(nèi)部生成蟲洞,導(dǎo)致一部分被驅(qū)替流體不能被波及,驅(qū)替流體沿蟲洞離開多孔介質(zhì),造成驅(qū)替效率的減少.在此基礎(chǔ)上,隨著達(dá)姆科勒數(shù)的增大,孔隙率變化越大,生成的蟲洞越寬,最終驅(qū)替效率變大,但仍小于無溶解反應(yīng)時(shí)的驅(qū)替效率;隨著佩克萊數(shù)的增大,指進(jìn)增長(zhǎng)速度越快,孔隙率變化越小,驅(qū)替效率越小.

1 引言

多孔介質(zhì)中含流固溶解反應(yīng)的互溶驅(qū)替過程在地下水污染治理[1]、二氧化碳埋存[2]及溫室氣體強(qiáng)化采油[3]等領(lǐng)域廣泛存在,是能源、環(huán)境等學(xué)科共同關(guān)注的典型問題.如在地下水污染治理中,污水會(huì)不斷侵蝕和驅(qū)替地下水,同時(shí)污水中的酸根離子會(huì)與地層中的巖石(如方解石、白云石等)發(fā)生溶解反應(yīng),使巖石結(jié)構(gòu)發(fā)生顯著變化,進(jìn)而影響流體流動(dòng).在此驅(qū)替過程中,流動(dòng)、擴(kuò)散以及化學(xué)反應(yīng)三者同時(shí)存在,相互耦合.如何有效地控制含溶解反應(yīng)下的驅(qū)替是相關(guān)領(lǐng)域的關(guān)鍵技術(shù)之一.

關(guān)于互溶驅(qū)替問題的研究由來已久.當(dāng)驅(qū)替流體黏性小于被驅(qū)替流體時(shí),可能產(chǎn)生黏性指進(jìn)現(xiàn)象.Saffman 和Taylor[4]采用Hele-Shaw 裝置對(duì)多孔介質(zhì)內(nèi)的驅(qū)替問題進(jìn)行實(shí)驗(yàn)研究,并提供了黏性指進(jìn)指端發(fā)展的數(shù)學(xué)模型;Zimmerman 和Homsy[5]通過數(shù)值模擬研究了二維非均質(zhì)多孔介質(zhì)中互溶驅(qū)替過程;Wit 和Homsy[6]研究了多孔介質(zhì)中流體間化學(xué)反應(yīng)與黏性指進(jìn)的非線性相互作用,發(fā)現(xiàn)反應(yīng)會(huì)使指進(jìn)形態(tài)特征和傳播機(jī)制發(fā)生顯著的變化;之后,Nagatsu 等[7]通過實(shí)驗(yàn)及理論方法研究了流體間沉淀反應(yīng)對(duì)指進(jìn)的影響,并討論了滲透率變化與指進(jìn)發(fā)展之間的關(guān)系.需要指出的是,上述研究涉及的化學(xué)反應(yīng)均是流體與流體之間的化學(xué)反應(yīng),少有對(duì)流體與多孔介質(zhì)骨架之間存在化學(xué)反應(yīng)的互溶驅(qū)替現(xiàn)象的研究.

另一方面,研究人員采用不同方法研究了多孔介質(zhì)中含化學(xué)反應(yīng)的各種問題.Békri 等[8]率先將隨機(jī)模型和矩量法相結(jié)合以計(jì)算固體壁面處的反應(yīng)速率;Luo 等[9]和Oltéan 等[10]通過arbitrary-Lagrangian-Eulerian (ALE)方法求解歐拉網(wǎng)格上的物理問題(如溶質(zhì)擴(kuò)散、界面反應(yīng)等),該方法需要及時(shí)更新歐拉網(wǎng)格;Soulaine 等[11]提出了一種基于Dancy-Brinkman-Stokes (DBS)方程的方法來模擬單相流動(dòng)下固體的溶解反應(yīng),并通過實(shí)驗(yàn)驗(yàn)證了該方法的有效性.由于孔隙尺度下多孔介質(zhì)內(nèi)固體骨架形狀復(fù)雜,尤其當(dāng)固體骨架隨著反應(yīng)發(fā)生而不斷變化時(shí),傳統(tǒng)數(shù)值方法處理起來難度較大.近三十年來,格子Boltzamnn 方法(lattice Boltzmann method,LBM)作為一種介觀的模擬方法,可以更方便地處理復(fù)雜微觀孔隙結(jié)構(gòu)內(nèi)流體與固體之間、不同流體組分之間復(fù)雜的相互作用[12-14],對(duì)于多孔介質(zhì)內(nèi)流動(dòng)與反應(yīng)耦合的問題,也受到了學(xué)者們的廣泛關(guān)注.如He 等[15]將流體-壁面間的表面反應(yīng)與擴(kuò)散機(jī)理相耦合,通過邊界條件實(shí)現(xiàn)了表面反應(yīng)的產(chǎn)生;Kang 等[16]在He 模型的基礎(chǔ)上建立了一種格子Boltzmann 模型,模擬了多孔介質(zhì)中兩種水合物的流動(dòng)與反應(yīng)過程,并考慮了多孔介質(zhì)幾何形狀的演化,該模擬結(jié)果與實(shí)驗(yàn)觀察的現(xiàn)象相符合;之后又進(jìn)一步研究了多組分溶液在多孔介質(zhì)中的反應(yīng)性傳輸[17];張婷等[18]詳細(xì)分析了各無量綱參數(shù)對(duì)多孔介質(zhì)內(nèi)溶解和沉淀過程的影響;Ju 等[19]提出了一種適用于不規(guī)則表面的反應(yīng)邊界格式,并成功將其拓展至三維模型;Meng 等[20]從多尺度分析了非均相反應(yīng)對(duì)密度驅(qū)動(dòng)流的影響,發(fā)現(xiàn)界面的不穩(wěn)定性由流體密度差、反應(yīng)速率和孔隙率變化三者共同決定.

綜上所述,雖然目前的研究對(duì)于多孔介質(zhì)互溶驅(qū)替過程和多孔介質(zhì)中含化學(xué)反應(yīng)的問題均有充足的認(rèn)識(shí),但對(duì)于多孔介質(zhì)內(nèi)含溶解反應(yīng)的互溶驅(qū)替過程研究尚不充分,此過程中的驅(qū)替規(guī)律尚待研究.本文采用LBM 對(duì)孔隙尺度下多孔介質(zhì)內(nèi)含溶解反應(yīng)的互溶驅(qū)替過程進(jìn)行研究,著重分析驅(qū)替過程中溶解反應(yīng)對(duì)多孔介質(zhì)結(jié)構(gòu)和驅(qū)替效率造成的影響,尋求控制流體驅(qū)替的方法,為地下水污染治理,二氧化碳混相驅(qū)油等技術(shù)提供思路.

2 問題描述

2.1 物理問題

本文研究孔隙尺度下多孔介質(zhì)中含溶解反應(yīng)的互溶流體驅(qū)替過程.其中,多孔介質(zhì)固體骨架結(jié)構(gòu)是通過對(duì)CT 掃描得到的多孔介質(zhì)灰度圖進(jìn)行二值化處理而獲得的,多孔介質(zhì)區(qū)域?qū)挾群透叨确謩e為lx和ly,lx=ly.多孔介質(zhì)孔隙率ε的定義為

其中Vpore為孔隙所占體積,Vtotal為多孔介質(zhì)的總體積.圖1 所示的多孔介質(zhì)中,藍(lán)黑色不規(guī)則部分代表的是多孔介質(zhì)固體骨架B,該多孔介質(zhì)初始孔隙率ε0為0.864,下標(biāo)“0”表示初始時(shí)刻.

圖1 系統(tǒng)初始狀態(tài)示意圖Fig.1.Schematic illustration of the initial state of the system.

系統(tǒng)設(shè)置如圖1 所示.在初始時(shí)刻,0.05lx≤x≤lx區(qū)域的多孔介質(zhì)中充滿了黏性為μ2的被驅(qū)替流體(Fluid 2),該流體處于靜止?fàn)顟B(tài).剩余部分的多孔介質(zhì)中充滿了黏性為μ1的驅(qū)替流體(Fluid 1),然后從多孔介質(zhì)左側(cè)以恒定速度u0持續(xù)注入該驅(qū)替流體,兩流體等溫,可互溶,且均為不可壓縮流體.在驅(qū)替流體中含有溶質(zhì)A,流體無量綱濃度為C1=1;被驅(qū)替流體中不含該溶質(zhì),流體無量綱濃度C2=0.在驅(qū)替過程中溶質(zhì)A與固體骨架B發(fā)生溶解反應(yīng),A+B→P,反應(yīng)產(chǎn)物P為溶于流體的酸根離子.

2.2 數(shù)學(xué)描述

在孔隙尺度下,多孔介質(zhì)內(nèi)流體的流動(dòng)可用如下方程組來刻畫:

其中(2)式為不可壓流體的連續(xù)性方程,(3)式為不可壓Navier-Stokes 方程.u為流體流過多孔介質(zhì)孔隙時(shí)的速度,t為時(shí)間,p為壓力,ν為流體的運(yùn)動(dòng)黏度.

多孔介質(zhì)內(nèi)兩種流體間溶質(zhì)的擴(kuò)散采用對(duì)流擴(kuò)散方程來描述.本文考慮的是流體與固體之間的溶解反應(yīng),不包括流體之間的化學(xué)反應(yīng).因此,對(duì)流擴(kuò)散方程中不包含源項(xiàng),即

其中,C為流體中溶質(zhì)的濃度,D為溶質(zhì)間的擴(kuò)散系數(shù).

系統(tǒng)邊界條件如下:當(dāng)x=0時(shí),u=(u0,0),C=C1;當(dāng)x=lx時(shí),上下邊界的速度場(chǎng)滿足無滑移條件,對(duì)于上下邊界的濃度場(chǎng),則有?C=0.

此外,對(duì)于多孔介質(zhì)固體骨架與溶質(zhì)間的溶解反應(yīng),也通過邊界條件實(shí)現(xiàn).考慮到實(shí)際工程中溶質(zhì)與固體間的溶解反應(yīng)機(jī)理過于復(fù)雜,本文對(duì)該反應(yīng)進(jìn)行簡(jiǎn)化處理,采用如下的線性反應(yīng)邊界條件[18]來描述:

其中,kr為溶解反應(yīng)的反應(yīng)速率常數(shù),n為固體壁面指向流體的法線方向,Cw為固體壁面處流體的濃度,Ceq代表平衡濃度,本文中溶解反應(yīng)是不可逆的,Ceq設(shè)為0.

該系統(tǒng)有4 組無量綱參數(shù),分別為雷諾數(shù)(Reynolds number,Re),Re=u0lx/ν2;佩克萊數(shù)(Pèlcet number,Pe),Pe=u0lx/D;黏性比(viscosity ratio,M)定義為兩流體的黏性之比,M=μ2/μ1;達(dá)姆科勒數(shù)(Damkohler number,Da),Da=krly/D.在驅(qū)替過程中,由于兩種流體互溶,流體黏性與濃度相關(guān),即

其中R=lnM.需要注意的是,當(dāng)驅(qū)替流體黏性遠(yuǎn)小于被驅(qū)替流體時(shí),會(huì)產(chǎn)生明顯的指進(jìn)現(xiàn)象[21,22].本文研究的驅(qū)替問題考慮的均是兩互溶流體黏性差異較大的情況,黏性比設(shè)定為M=100.

3 格子Boltzmann 方法

本文采用多孔介質(zhì)中流體流動(dòng)及擴(kuò)散的耦合格子Boltzmann (lattice Boltzmann,LB)模型[23]來求解(2)式—(4)式組成的方程組.該模型為雙分布模型,用兩組演化方程分別描述速度場(chǎng)和濃度場(chǎng).

3.1 耦合LB 模型

速度場(chǎng)可用如下演化方程來描述:

其中q是離散速度的個(gè)數(shù),fi(x,t) 是粒子在t時(shí)刻x處離散速度為ci的分布函數(shù),δt為時(shí)間步長(zhǎng),Ωi是碰撞矩陣.流體的密度和速度可以表示為

已有研究證明多松弛(multiple-relaxtiontime,MRT)模型在模擬孔隙尺度下多孔介質(zhì)內(nèi)流動(dòng)問題時(shí),相比于Bhatnagar-Gross-Krook (BGK)模型,具有更好的數(shù)值穩(wěn)定性[24].故本文采用MRT模型求解流動(dòng)方程.

MRT 模型的碰撞算子可表示為

其中m和m(eq)分別是矩向量和矩空間中的平衡態(tài),M是將分布函數(shù)映射到矩空間的q×q的變換矩陣,m=M ·f,S=diag(s0,s1,···,sq-1) 為松弛矩陣.

本文考慮的是二維問題,采用二維九速(D2Q9)離散速度模型,該模型中離散速度ci定義為

其中c=δx/δt,δx和δt分別為網(wǎng)格步長(zhǎng)和時(shí)間步長(zhǎng),本文取c=1.

變換矩陣M為

其對(duì)應(yīng)的矩向量為

其中T 表示轉(zhuǎn)置變換.(12)式中的各個(gè)矩的物理含義如下:ρ是流體密度;e和ε分別是總能和能量的平方;jx和jy分別為x和y分量上的動(dòng)量,jx=ρux,jy=ρuy;qx和qy則是x和y分量上的能量通量;pxx和pxy分別為應(yīng)力張量的對(duì)角線和非對(duì)角線上的元素.對(duì)于不可壓流體,流體密度近似為常數(shù),記作ρ0,而密度波動(dòng)為 δρ,從而有ρ=ρ0+δρ.相應(yīng)矩空間中的平衡態(tài)如下:

其中,流體密度ρ和動(dòng)量j=ρ0u=(jx,jy) 是碰撞守恒量.而其他非守恒的平衡態(tài)是守恒量的函數(shù),具有如下形式:

此模型的松弛參數(shù)為

注意到流體密度和動(dòng)量在碰撞過程中守恒,這些矩對(duì)應(yīng)的松弛率sρ和sj可取任意值.其他松弛率為se=sε=sν=sν=1/τ,sq=8(2-sν)/(8-sν)[24],其中τ與流體動(dòng)力黏性ν相關(guān),

濃度場(chǎng)則用另一組演化方程來描述:

其中g(shù)i(x,t) 為濃度C的分布函數(shù),離散速度和上文中用于描述速度場(chǎng)的MRT 模型中的離散速度一致;平衡態(tài)分布函數(shù)為

其中ωi為權(quán)系數(shù),ω0=4/9,ω1-4=1/9,ω5-8=1/36,而=c2/3.參數(shù)A與擴(kuò)散系數(shù)D相關(guān):

對(duì)于擴(kuò)散系數(shù)較小(或Pe數(shù)較大)的情況,可通過調(diào)整參數(shù)A來提高數(shù)值穩(wěn)定性.

此模型中,濃度定義如下:

(23)式中的濃度梯度可由該點(diǎn)非平衡態(tài)函數(shù)的一階矩得到,計(jì)算如下:

3.2 邊界處理

對(duì)于不同的邊界條件,均需進(jìn)行邊界處理.其中左右邊界均采用非平衡外推格式進(jìn)行處理[25].對(duì)于上下邊界和多孔介質(zhì)內(nèi)部的固體骨架,速度滿足無滑移條件,采用半反彈格式來處理[26].標(biāo)準(zhǔn)反彈格式固體邊界位于格點(diǎn)上,而半反彈格式的固體邊界位于格點(diǎn)的中間,即 (xf+xb)/2處,其中xb是固體格點(diǎn),xf是臨近固體壁面的流體格點(diǎn).

其中 δt為時(shí)間步長(zhǎng),c?=-ci(ci指向壁面),(xf,t)為fi(xf,t) 碰撞后的值.

對(duì)于上下邊界的濃度場(chǎng),?C=0.可用一種類反彈格式來實(shí)現(xiàn)[27]:

多孔介質(zhì)固體骨架與溶質(zhì)間的溶解反應(yīng),用(5)式來表示,采用修正半反彈格式[28]來處理:

其中Cw可由有限差分格式得出:

其中,Cf為固體壁面附近流體節(jié)點(diǎn)的濃度.

溶解反應(yīng)在固相上的表現(xiàn)為固體體積的不斷變化,同時(shí)假定固體骨架主要由易溶鹽礦物組成,能夠全部溶解.節(jié)點(diǎn)處固體體積隨時(shí)間的變化規(guī)律滿足公式[17]:

其中,V表示多孔介質(zhì)固體骨架的體積,Vm表示其摩爾體積,as表示化學(xué)計(jì)量數(shù),此處即反應(yīng)表面積.每個(gè)固體節(jié)點(diǎn)表示為一個(gè)控制體積,初始時(shí)其大小為 1×1 的單位網(wǎng)格單元.在發(fā)生反應(yīng)時(shí),固體節(jié)點(diǎn)的體積在每個(gè)時(shí)刻依據(jù)以下式子更新:

當(dāng)固體節(jié)點(diǎn)的體積V=0 時(shí),說明該部分固體已被完全溶解,節(jié)點(diǎn)被更新為流體節(jié)點(diǎn),其濃度由附近流體點(diǎn)濃度的平均值求得,而速度則取為0.

3.3 模型驗(yàn)證

本節(jié)通過模擬兩平板間的互溶驅(qū)替過程,對(duì)所使用的耦合LB 模型進(jìn)行驗(yàn)證.如圖2 所示,計(jì)算區(qū)域?yàn)?0≤x≤lx,0≤y≤ly;初始時(shí)刻,計(jì)算區(qū)域內(nèi)0.1lx≤x≤lx,0≤y≤ly的部分充滿了濃度為C=C2=0,動(dòng)力黏度為μ=μ2的被驅(qū)替流體,剩余部分充滿了濃度為C=C1=1,動(dòng)力黏度為μ=μ1的驅(qū)替流體;之后,從計(jì)算區(qū)域左邊界以速度u=(0.005,0)持續(xù)注入該驅(qū)替流體.本節(jié)采用的網(wǎng)格大小為 256×32.

圖2 兩平板間的互溶驅(qū)替問題初始狀態(tài)示意圖Fig.2.Schematic illustration of the initial status of the miscible displacement of two parallel plates.

在此驅(qū)替過程中,系統(tǒng)的控制方程為(2)式—(4)式組成的方程組,本節(jié)的邊界條件給定如下:當(dāng)x=0時(shí),u=(u0,0),c=c1;當(dāng)x=lx時(shí),上下邊界為壁面,速度滿足無滑移條件,對(duì)于濃度場(chǎng),?C=0.本節(jié)模擬了文獻(xiàn)[21]中的4 種工況:1)M=1,Pe=5;2)M=1,Pe=262;3)M=100,Pe=5;4)M=100,Pe=262.模擬結(jié)果如圖3 所示,圖中還給出了縱截面的平均濃度分布圖,通過對(duì)比可以發(fā)現(xiàn)與文獻(xiàn)中結(jié)果完全一致,即驗(yàn)證了該耦合LB 模型模擬互溶驅(qū)替過程的準(zhǔn)確性.

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

4.1 網(wǎng)格無關(guān)性驗(yàn)證

為確保選取網(wǎng)格的計(jì)算精度和計(jì)算效率,在數(shù)值模擬之前要進(jìn)行網(wǎng)格無關(guān)性的驗(yàn)證.計(jì)算區(qū)域的長(zhǎng)和寬為lx=ly=32,分別選擇 512×512 和1024×1024的網(wǎng)格尺寸來描述圖1 中的多孔介質(zhì)系統(tǒng).在本次模擬中,設(shè)定u0=(0.005,0),Re=32,Pe=262及Da=0.為量化不同網(wǎng)格間的差異,選取y=16 處的橫截面作為觀測(cè)面,圖4 為觀測(cè)面上流體的平均濃度和平均速度隨時(shí)間的演變圖.需要注意的是,在本文的模擬中,特征時(shí)間τ=ly/u0,無量綱時(shí)間t=δtNstep/τ,其中Nstep為模擬所用時(shí)間步,后文中t均為無量綱時(shí)間.如圖4 所示,在512×512和 1024×1024 這兩種網(wǎng)格下演化過程符合良好,因此選取 512×512 的網(wǎng)格即可滿足網(wǎng)格無關(guān).

圖4 不同網(wǎng)格下橫界面(a)平均濃度和(b)平均速度隨時(shí)間的演變圖Fig.4.Time evolution of the transverse-average concentration (a) and transverse-average velocity (b) with different grids.

4.2 達(dá)姆科勒數(shù)的影響

圖34 種工況下的互溶驅(qū)替現(xiàn)象和平均濃度分布圖Fig.3.Miscible displacement phenomena and the average concentration distributions of four different cases.

在含溶解反應(yīng)的驅(qū)替過程中,Da數(shù)為表征溶解反應(yīng)速率的無量綱參數(shù),本節(jié)將研究Da數(shù)對(duì)驅(qū)替過程的影響.首先模擬了3 種不同工況下的互溶驅(qū)替過程,工況一為無溶解反應(yīng)(Da=0)的互溶驅(qū)替過程;工況二為含有溶解反應(yīng)且反應(yīng)速率較小(Da=4)的互溶驅(qū)替;工況三中溶解反應(yīng)速率較大(Da=32).雷諾數(shù)設(shè)定為Re=32,佩克萊數(shù)設(shè)定為Pe=262,黏性比設(shè)定為M=100.圖5 給出了3 種工況下不同時(shí)刻流體的濃度場(chǎng)圖.需要指出的是,濃度場(chǎng)圖中藍(lán)黑色不規(guī)則部分代表的是多孔介質(zhì)固體骨架B,隨著溶解反應(yīng)的發(fā)生,部分固體骨架被溶解,多孔介質(zhì)結(jié)構(gòu)隨之發(fā)生改變.

如圖5(a)所示,在驅(qū)替過程早期(t=0.0488),驅(qū)替流體剛進(jìn)入多孔介質(zhì),與多孔介質(zhì)固體骨架的接觸面積較小,溶解反應(yīng)的量也較小,因此3 種工況下流體的濃度場(chǎng)之間無明顯差異.

隨著驅(qū)替過程的進(jìn)行(t=0.2441,見圖5(b)),受驅(qū)替流體和被驅(qū)替液流體間的黏性差異以及多孔介質(zhì)孔隙結(jié)構(gòu)兩種因素的影響,在兩流體界面上產(chǎn)生黏性指進(jìn)現(xiàn)象.此外,在工況二及工況三中,驅(qū)替流體與多孔介質(zhì)固體骨架發(fā)生流固溶解反應(yīng),部分固體骨架逐漸被溶解,擴(kuò)大了多孔介質(zhì)的孔隙率.但由于溶解反應(yīng)的發(fā)生消耗了一部分溶質(zhì),因此Da數(shù)較大時(shí)指進(jìn)發(fā)展的速度相對(duì)較慢.總體上,此時(shí)3 種工況下流體的濃度場(chǎng)差異仍不大.

隨后(t=0.4883,見圖5(c)),隨著驅(qū)替流體不斷向前推進(jìn),工況二和工況三中溶解反應(yīng)沿指進(jìn)發(fā)展路徑不斷發(fā)生,從而在多孔介質(zhì)內(nèi)部形成裂縫,使后續(xù)驅(qū)替流體沿裂縫前進(jìn).與此同時(shí),在工況三中,驅(qū)替流體與多孔介質(zhì)固體骨架間持續(xù)發(fā)生溶解反應(yīng),促進(jìn)了裂縫的擴(kuò)展,形成兩條明顯的通道,即蟲洞現(xiàn)象[29-31],而工況二中由于溶解反應(yīng)速率較小,沒有出現(xiàn)明顯的蟲洞.由此可以推測(cè),存在一個(gè)臨界Dacr數(shù),當(dāng)Da數(shù)達(dá)到臨界值時(shí),互溶驅(qū)替過程中由于流固溶解反應(yīng)的發(fā)生會(huì)使多孔介質(zhì)中出現(xiàn)蟲洞.工況二中由于Da數(shù)未達(dá)到臨界值,所以驅(qū)替過程和工況一類似.本文將在下文中對(duì)臨界Dacr數(shù)作進(jìn)一步的討論.此外,隨著驅(qū)替的不斷進(jìn)行,工況一和工況二中驅(qū)替流體仍以指狀不斷前進(jìn),且隨著溶質(zhì)的擴(kuò)散,驅(qū)替流體與被驅(qū)替流體之間的界面逐漸增厚,與工況三形成明顯差異.另外要注意的是,工況二和工況三中流固溶解反應(yīng)的發(fā)生同時(shí)會(huì)降低驅(qū)替流體的濃度,減小驅(qū)替流體與被驅(qū)替流體間的黏性差異,從而減少黏性指進(jìn)的持續(xù)發(fā)展.因此,工況二中指進(jìn)的發(fā)展較工況一略有減少.

圖53 種工況下不同時(shí)刻多孔介質(zhì)內(nèi)流體濃度場(chǎng)瞬時(shí)圖Fig.5.Snapshots of the concentration field in porous media at different times under three cases.

最后(t=0.7324,見圖5(d)),工況一和工況二中驅(qū)替流體繼續(xù)沿著指進(jìn)的路徑不斷驅(qū)離多孔介質(zhì)區(qū)域,且工況二中指進(jìn)數(shù)量明顯少于工況一;而工況三中驅(qū)替流體則沿著蟲洞通道驅(qū)離.在3 種工況中,隨著溶質(zhì)間的擴(kuò)散仍不斷進(jìn)行,指進(jìn)和蟲洞均被不斷拓寬.

為研究臨界Dacr數(shù)的存在,圖6 給出了t=0.7324時(shí)刻Da=4 和8 的濃度場(chǎng)圖.可以看出,Da=4時(shí)由于反應(yīng)速率較小,驅(qū)替過程中無明顯的蟲洞現(xiàn)象;通過實(shí)驗(yàn)發(fā)現(xiàn),當(dāng)Da=8 時(shí),驅(qū)替過程中可見蟲洞的形成.由此可推出,存在臨界達(dá)姆科勒數(shù)Dacr,當(dāng)Da≥Dacr時(shí),驅(qū)替過程中會(huì)有蟲洞現(xiàn)象生成.由于驅(qū)替現(xiàn)象會(huì)隨著多孔介質(zhì)結(jié)構(gòu)的改變而呈現(xiàn)出不同的細(xì)節(jié),不同多孔介質(zhì)下Dacr的具體數(shù)值會(huì)有所差異.本節(jié)各工況中的多孔介質(zhì)結(jié)構(gòu)見圖1,對(duì)應(yīng)的Dacr≈8.

圖6 Da=4,8 的濃度場(chǎng)圖(t=0.7324)Fig.6.Concentration contours for Da=4,8 at t=0.7324.

下文以Dacr為界,討論不同Da數(shù)下,多孔介質(zhì)孔隙率ε和驅(qū)替效率η隨時(shí)間的變化情況,其余無量綱數(shù)與4.1 節(jié)保持一致.圖7(a)為Da=0,2,4,8時(shí)多孔介質(zhì)孔隙率隨時(shí)間的變化曲線.可以看出,當(dāng)Da=0,即無溶解反應(yīng)時(shí),多孔介質(zhì)孔隙率保持不變;當(dāng)Da=2,4,8,即溶解反應(yīng)較小時(shí),由于反應(yīng)的發(fā)生,孔隙率隨時(shí)間不斷增大.且隨著Da數(shù)的增大,同時(shí)刻多孔介質(zhì)的孔隙率越大.圖7(b)為Da=16,32,160,320 時(shí)多孔介質(zhì)孔隙率隨時(shí)間的變化曲線,此時(shí)溶解反應(yīng)會(huì)生成明顯的蟲洞.從圖7(b)可以看出,與Da數(shù)較小時(shí)不同的是,孔隙率隨時(shí)間的變化趨勢(shì)為s 型曲線.初始階段,驅(qū)替流體與多孔介質(zhì)固體骨架接觸面積小,溶解多孔介質(zhì)骨架的量小,多孔介質(zhì)孔隙率緩慢增長(zhǎng);伴隨驅(qū)替流體的逐漸進(jìn)入,驅(qū)替流體與多孔介質(zhì)接觸面積變大,溶解多孔介質(zhì)骨架的量逐漸變多,多孔介質(zhì)孔隙率越來越大,且孔隙率的增長(zhǎng)也逐漸加快,在t=0.35左右孔隙率的變化速率達(dá)到頂峰,此時(shí)為驅(qū)替流體初次到達(dá)流場(chǎng)出口的時(shí)間;之后,由于大量固體骨架已經(jīng)被溶解,驅(qū)替流體只能持續(xù)性拓寬蟲洞周圍,因此孔隙率變化較為緩慢,這在圖5(d)中也可體現(xiàn).此外,隨著Da數(shù)的增大,同時(shí)刻多孔介質(zhì)的孔隙率越大,這與Da數(shù)較小時(shí)情況相同.這是因?yàn)楫?dāng)Re數(shù),Pe數(shù)和黏性比M保持不變時(shí),增大Da數(shù)等價(jià)于增大反應(yīng)速率常數(shù)kr.反應(yīng)速率越大,相同時(shí)間內(nèi)溶解的多孔介質(zhì)骨架越多,因此孔隙率越大.具體地,當(dāng)t=0.7324,Da=4 時(shí)孔隙率為0.9122,較初始孔隙率提高了5.49%;當(dāng)t=0.7324,Da=32 時(shí)孔隙率為0.9421,較初始時(shí)刻孔隙率提高了9.03%;當(dāng)t=0.7324,Da=320 時(shí)孔隙率為0.9524,較初始時(shí)刻孔隙率提高10.22%,說明此時(shí)Da數(shù)的影響已經(jīng)減弱.

圖7 不同Da 數(shù)下,(a),(b)多孔介質(zhì)孔隙率和(c),(d)驅(qū)替效率隨時(shí)間的變化圖Fig.7.Time evolutions of (a),(b) the porosity and (c),(d) displacement efficiency for different Damkohler numbers.

為了定量描述驅(qū)替效率,定義驅(qū)替效率

其中V(t)表示在t時(shí)刻驅(qū)替流體所占的體積,Vtotal為多孔介質(zhì)的總體積.需要說明的是,驅(qū)替效率應(yīng)為驅(qū)替流體所占體積與流體所占總體積之比.本文中溶解反應(yīng)的發(fā)生會(huì)改變多孔介質(zhì)結(jié)構(gòu),從而改變流體所占的總體積.為使驅(qū)替效率中所對(duì)比的量保持一致,本文設(shè)定驅(qū)替效率為驅(qū)替流體所占的體積與多孔介質(zhì)總體積之比.此外,由于溶質(zhì)間會(huì)發(fā)生擴(kuò)散,此處定義溶質(zhì)濃度C≥0.5 時(shí)為驅(qū)替流體.

圖7(c)為Da=0,2,4,8 時(shí)驅(qū)替效率隨時(shí)間的變化曲線.該過程可以分前后兩個(gè)部分進(jìn)行討論.在t=0—0.35 時(shí),驅(qū)替流體并未流出多孔介質(zhì)區(qū)域,驅(qū)替效率增長(zhǎng)較快.由于此時(shí)Da數(shù)較小,不同Da數(shù)下驅(qū)替效率相差不大.當(dāng)t ＞0.35 時(shí),部分驅(qū)替流體已流出多孔介質(zhì)區(qū)域,驅(qū)替效率在微降后緩慢增長(zhǎng).此階段,隨著Da數(shù)的增大,驅(qū)替效率越小.這是因?yàn)镈a數(shù)較大時(shí),指進(jìn)波及范圍減少,蟲洞現(xiàn)象逐漸明顯(如圖6 所示),驅(qū)替效率越低.

圖7(d)為Da=16,32,160,320 時(shí)驅(qū)替效率隨時(shí)間的變化曲線.當(dāng)t=0—0.35 時(shí),驅(qū)替效率同樣快速增長(zhǎng).Da數(shù)越大,驅(qū)替效率越小.當(dāng)t ＞0.35時(shí)驅(qū)替效率會(huì)有一個(gè)小幅下降的過程.驅(qū)替效率出現(xiàn)短時(shí)下降一方面是因?yàn)橐徊糠烛?qū)替流體已經(jīng)離開了流場(chǎng),另一方面是因?yàn)榇藭r(shí)驅(qū)替流體波及范圍較廣,溶解反應(yīng)發(fā)生部位較多,反應(yīng)消耗的驅(qū)替流體較多,導(dǎo)致流入多孔介質(zhì)驅(qū)替流體的量少于流出的量加上反應(yīng)的量,驅(qū)替效率下降.隨著時(shí)間的發(fā)展,如圖5 所示,驅(qū)替流體波及范圍逐漸增大,驅(qū)替效率會(huì)逐漸變大.另外,在驅(qū)替過程后期,隨著Da數(shù)的增大,驅(qū)替效率也越大(如圖7(d)中插圖所示).這是因?yàn)镈a數(shù)越大,生成的蟲洞相對(duì)越寬,驅(qū)替流體波及范圍相對(duì)越廣,驅(qū)替效率越大.具體地,在t=0.7324時(shí)刻,Da=32時(shí)驅(qū)替效率為0.3478,而Da=320 的驅(qū)替效率為0.3710.可見雖然Da數(shù)越大驅(qū)替效率相應(yīng)越大,但驅(qū)替效率增加幅度不明顯.

綜上所述,當(dāng)流場(chǎng)內(nèi)未生成明顯蟲洞時(shí),Da數(shù)越大,驅(qū)替效率越小;當(dāng)流場(chǎng)內(nèi)生成明顯蟲洞時(shí),隨著時(shí)間發(fā)展,Da數(shù)越大,驅(qū)替效率會(huì)越大,但是仍遠(yuǎn)小于無反應(yīng)時(shí)的驅(qū)替效率.

4.3 佩克萊數(shù)的影響

Pe數(shù)表征驅(qū)替過程中對(duì)流速率與擴(kuò)散速率的相對(duì)大小,具有重要研究意義.模擬中設(shè)定Da=32,通過改變Pe數(shù)來觀察其對(duì)驅(qū)替過程的影響.圖8 給出了Pe=5,65,524,t=0.2441 時(shí)的濃度場(chǎng)瞬時(shí)圖.可以看出,在Pe數(shù)很小時(shí)(Pe=5,圖8(a)),流體的擴(kuò)散能力強(qiáng),使得溶質(zhì)可充分?jǐn)U散,在兩流體間的界面上無法產(chǎn)生指進(jìn)現(xiàn)象.此時(shí)驅(qū)替流體以平界面形式沿驅(qū)替方向溶解固體骨架,圖中左半側(cè)的多孔介質(zhì)骨架已被溶解.隨著Pe數(shù)的逐步增大(Pe=65,圖8(b)),流體間溶質(zhì)的擴(kuò)散減弱,對(duì)流作用增強(qiáng),開始產(chǎn)生指進(jìn)現(xiàn)象,但此時(shí)指進(jìn)個(gè)數(shù)較少,多孔介質(zhì)骨架被溶解的量也隨之減少.而當(dāng)Pe數(shù)較大時(shí)(Pe=524,圖8(c)),出現(xiàn)了明顯的指進(jìn)現(xiàn)象,一些指進(jìn)被多孔介質(zhì)骨架分裂開來,形成一系列更小的指進(jìn),即發(fā)生了指進(jìn)分裂現(xiàn)象.其他研究人員也對(duì)此類現(xiàn)象進(jìn)行了研究[32].此外,注意到此工況下多孔介質(zhì)骨架被溶解的量最少.

圖8 t=0.2441時(shí)不同Pe 數(shù)下的濃度場(chǎng)瞬時(shí)圖 (a)Pe=5;(b) Pe=65;(c)Pe=524Fig.8.Snapshots of the concentration field at t=0.2441 :(a) Pe=5;(b) Pe=65;(c) Pe=524.

圖9 為相應(yīng)情形下多孔介質(zhì)內(nèi)水平方向上的平均濃度圖.可以看出,Pe=5 時(shí),平均濃度的變化十分平滑.隨著Pe數(shù)的增大(Pe=65),濃度曲線有了細(xì)微的波動(dòng),說明此時(shí)已有指進(jìn)生成.當(dāng)Pe數(shù)為524 時(shí),濃度曲線上出現(xiàn)了明顯的波動(dòng),說明高Pe數(shù)下驅(qū)替過程中產(chǎn)生大量指進(jìn),流體界面不穩(wěn)定性增加.

圖9 不同Pe 數(shù)下多孔介質(zhì)內(nèi)水平方向上的平均濃度Fig.9.Average concentration along the horizontal direction of the porous media for different Pèclet numbers.

接下來選取指進(jìn)現(xiàn)象較為明顯的情況作進(jìn)一步討論.圖10 為Pe=131,262,524 時(shí),多孔介質(zhì)的孔隙率和驅(qū)替效率隨時(shí)間的變化情況.如圖10(a)所示,Pe數(shù)越大,同一時(shí)刻下孔隙率越小.這是因?yàn)镻e數(shù)越大,擴(kuò)散的作用越小,指進(jìn)現(xiàn)象越明顯,指進(jìn)的增長(zhǎng)速度也越快,溶解的多孔介質(zhì)骨架越少,孔隙率越小.此外,不同Pe數(shù)下孔隙率的變化趨勢(shì)相同,均為s 型曲線.如圖10(b)所示,當(dāng)Pe數(shù)越大時(shí),初始階段驅(qū)替效率越大.這是因?yàn)镻e數(shù)越大,同時(shí)刻下指進(jìn)數(shù)量越多,驅(qū)替流體波及范圍越廣(如圖10(b)中插圖①和②所示),驅(qū)替效率相應(yīng)越大.但隨著時(shí)間的發(fā)展,Pe數(shù)越小,蟲洞越寬,且入口和蟲洞附近溶解的固體骨架越多(如圖10(b)中插圖③和④所示),驅(qū)替流體占據(jù)更多的空間,驅(qū)替效率越大.

圖10 不同Pe 數(shù)下多孔介質(zhì)的(a)孔隙率和(b)驅(qū)替效率隨時(shí)間的變化Fig.10.Time evolutions of (a) porosity and (b) displacement efficiency for different Pèlcet numbers.

4.4 多孔介質(zhì)結(jié)構(gòu)的影響

在上文的模擬中,所選取的多孔介質(zhì)均如圖1所示.為驗(yàn)證上述結(jié)果具有普適性,本節(jié)選取同一多孔介質(zhì)CT 掃描灰質(zhì)圖中的其他部分進(jìn)行二值化處理,所獲取的多孔介質(zhì)初始孔隙率分別為ε0=0.8和0.75,并分別模擬了兩種多孔介質(zhì)中含溶解反應(yīng)的驅(qū)替過程.圖11 為兩種多孔介質(zhì)下孔隙率和驅(qū)替效率隨時(shí)間的變化圖.可以看出,孔隙率的變化雖然會(huì)顯著影響孔隙率和驅(qū)替效率的數(shù)值大小,但并不會(huì)改變整體的變化趨勢(shì),即孔隙率變化依然呈s 型,驅(qū)替效率隨時(shí)間先增大然后減小最后增大.這表明上文研究的驅(qū)替過程的基本規(guī)律在不同的孔隙率下仍然適用,研究結(jié)果不失一般性.

圖11 不同初始孔隙率下,多孔介質(zhì)(a)孔隙率和(b)驅(qū)替效率隨時(shí)間的變化情況Fig.11.Time evolutions of (a) porosity and (b) displacement efficiency for different initial porosities.

此外,如圖11(b)所示,隨著多孔介質(zhì)初始孔隙率的增加,驅(qū)替效率也相應(yīng)增加.這是由于同一多孔介質(zhì)中,初始孔隙率越大的部分,流體越容易通過多孔介質(zhì),為驅(qū)替的進(jìn)行創(chuàng)造更好的條件,從而為溶解反應(yīng)在流固界面的進(jìn)行也提供了更多的機(jī)會(huì),使得驅(qū)替效率隨之增加.因此不同多孔介質(zhì)結(jié)構(gòu)對(duì)含溶解反應(yīng)的驅(qū)替有著顯著影響,初始孔隙率較大的多孔介質(zhì)中驅(qū)替效率也越高.

5 結(jié)論

本文采用LBM 模擬了孔隙尺度下多孔介質(zhì)內(nèi)含有溶解反應(yīng)的互溶驅(qū)替過程,觀察了含流固溶解反應(yīng)時(shí)多孔介質(zhì)內(nèi)流體濃度場(chǎng)隨時(shí)間的演化情況,并與不含溶解反應(yīng)的工況進(jìn)行對(duì)比,研究溶解反應(yīng)對(duì)驅(qū)替過程的影響;之后通過改變Da數(shù)和Pe數(shù)研究孔隙率和驅(qū)替效率的變化情況,得到如下結(jié)論.

1)在兩流體黏性差異較大的前提下,在驅(qū)替過程早期,溶解反應(yīng)會(huì)消耗驅(qū)替流體的濃度,抑制指進(jìn)的繼續(xù)發(fā)展;隨著時(shí)間的增加,溶解反應(yīng)會(huì)不斷溶解多孔介質(zhì)骨架,在多孔介質(zhì)內(nèi)部生成蟲洞,驅(qū)替流體會(huì)直接沿著該蟲洞離開多孔介質(zhì),減小了驅(qū)替效率.

2)隨著Da數(shù)的增大,反應(yīng)速率變大,孔隙率變化越大,所生成的蟲洞越寬,驅(qū)替效率逐漸增大.

3)隨著Pe數(shù)的增大,指進(jìn)增長(zhǎng)的速率越快,反應(yīng)的量越小,孔隙率變化越小,驅(qū)替效率也隨之變小.

4)不同的多孔介質(zhì)結(jié)構(gòu)下獲得的孔隙率和驅(qū)替效率隨時(shí)間的變化趨勢(shì)一致,說明以上結(jié)論具有普適性.此外,初始孔隙率較大的多孔介質(zhì)中驅(qū)替效率也越高.

在互溶驅(qū)替過程中,若兩種流體黏性差異較大且無法避免流固溶解反應(yīng)發(fā)生的情況下,為獲取較高的驅(qū)替效率,應(yīng)選擇與多孔介質(zhì)發(fā)生溶解反應(yīng)速率高的驅(qū)替流體,但無需一味追求反應(yīng)速率.此外,系統(tǒng)的Pe數(shù)越小驅(qū)替效率越高.為減緩驅(qū)替,可采取相反的策略.