多孔介質(zhì)中含溶解反應(yīng)的互溶驅(qū)替過程格子Boltzmann 研究*

劉高潔 邵子宇 婁欽

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

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

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

1 引言

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

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

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

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

2 問題描述

2.1 物理問題

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

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

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

系統(tǒng)設(shè)置如圖1 所示.在初始時刻,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)流體的流動可用如下方程組來刻畫:

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

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

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

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

此外,對于多孔介質(zhì)固體骨架與溶質(zhì)間的溶解反應(yīng),也通過邊界條件實現(xiàn).考慮到實際工程中溶質(zhì)與固體間的溶解反應(yīng)機理過于復(fù)雜,本文對該反應(yīng)進行簡化處理,采用如下的線性反應(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;達姆科勒數(shù)(Damkohler number,Da),Da=krly/D.在驅(qū)替過程中,由于兩種流體互溶,流體黏性與濃度相關(guān),即

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

3 格子Boltzmann 方法

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

3.1 耦合LB 模型

速度場可用如下演化方程來描述:

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

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

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)格步長和時間步長,本文取c=1.

變換矩陣M為

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

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

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

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

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

濃度場則用另一組演化方程來描述:

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

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

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

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

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

3.2 邊界處理

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

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

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

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

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

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

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

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

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

3.3 模型驗證

本節(jié)通過模擬兩平板間的互溶驅(qū)替過程,對所使用的耦合LB 模型進行驗證.如圖2 所示,計算區(qū)域為 0≤x≤lx,0≤y≤ly;初始時刻,計算區(qū)域內(nèi)0.1lx≤x≤lx,0≤y≤ly的部分充滿了濃度為C=C2=0,動力黏度為μ=μ2的被驅(qū)替流體,剩余部分充滿了濃度為C=C1=1,動力黏度為μ=μ1的驅(qū)替流體;之后,從計算區(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時,u=(u0,0),c=c1;當(dāng)x=lx時,上下邊界為壁面,速度滿足無滑移條件,對于濃度場,?C=0.本節(jié)模擬了文獻[21]中的4 種工況:1)M=1,Pe=5;2)M=1,Pe=262;3)M=100,Pe=5;4)M=100,Pe=262.模擬結(jié)果如圖3 所示,圖中還給出了縱截面的平均濃度分布圖,通過對比可以發(fā)現(xiàn)與文獻中結(jié)果完全一致,即驗證了該耦合LB 模型模擬互溶驅(qū)替過程的準確性.

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

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

為確保選取網(wǎng)格的計算精度和計算效率,在數(shù)值模擬之前要進行網(wǎng)格無關(guān)性的驗證.計算區(qū)域的長和寬為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 處的橫截面作為觀測面,圖4 為觀測面上流體的平均濃度和平均速度隨時間的演變圖.需要注意的是,在本文的模擬中,特征時間τ=ly/u0,無量綱時間t=δtNstep/τ,其中Nstep為模擬所用時間步,后文中t均為無量綱時間.如圖4 所示,在512×512和 1024×1024 這兩種網(wǎng)格下演化過程符合良好,因此選取 512×512 的網(wǎng)格即可滿足網(wǎng)格無關(guān).

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

4.2 達姆科勒數(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ù)對驅(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 種工況下不同時刻流體的濃度場圖.需要指出的是,濃度場圖中藍黑色不規(guī)則部分代表的是多孔介質(zhì)固體骨架B,隨著溶解反應(yīng)的發(fā)生,部分固體骨架被溶解,多孔介質(zhì)結(jié)構(gòu)隨之發(fā)生改變.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖8 t=0.2441時不同Pe 數(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 時,平均濃度的變化十分平滑.隨著Pe數(shù)的增大(Pe=65),濃度曲線有了細微的波動,說明此時已有指進生成.當(dāng)Pe數(shù)為524 時,濃度曲線上出現(xiàn)了明顯的波動,說明高Pe數(shù)下驅(qū)替過程中產(chǎ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.

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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

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