多孔介質中含溶解反應的互溶驅替過程格子Boltzmann 研究*

劉高潔 邵子宇 婁欽

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

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

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

1 引言

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

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

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

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

2 問題描述

2.1 物理問題

本文研究孔隙尺度下多孔介質中含溶解反應的互溶流體驅替過程.其中,多孔介質固體骨架結構是通過對CT 掃描得到的多孔介質灰度圖進行二值化處理而獲得的,多孔介質區域寬度和高度分別為lx和ly,lx=ly.多孔介質孔隙率ε的定義為

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

圖1 系統初始狀態示意圖Fig.1.Schematic illustration of the initial state of the system.

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

2.2 數學描述

在孔隙尺度下,多孔介質內流體的流動可用如下方程組來刻畫:

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

多孔介質內兩種流體間溶質的擴散采用對流擴散方程來描述.本文考慮的是流體與固體之間的溶解反應,不包括流體之間的化學反應.因此,對流擴散方程中不包含源項,即

其中,C為流體中溶質的濃度,D為溶質間的擴散系數.

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

此外,對于多孔介質固體骨架與溶質間的溶解反應,也通過邊界條件實現.考慮到實際工程中溶質與固體間的溶解反應機理過于復雜,本文對該反應進行簡化處理,采用如下的線性反應邊界條件[18]來描述:

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

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

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

3 格子Boltzmann 方法

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

3.1 耦合LB 模型

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

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

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

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

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

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

其中c=δx/δt,δx和δt分別為網格步長和時間步長,本文取c=1.

變換矩陣M為

其對應的矩向量為

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

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

此模型的松弛參數為

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

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

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

其中ωi為權系數,ω0=4/9,ω1-4=1/9,ω5-8=1/36,而=c2/3.參數A與擴散系數D相關:

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

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

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

3.2 邊界處理

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

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

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

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

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

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

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

其中,V表示多孔介質固體骨架的體積,Vm表示其摩爾體積,as表示化學計量數,此處即反應表面積.每個固體節點表示為一個控制體積,初始時其大小為 1×1 的單位網格單元.在發生反應時,固體節點的體積在每個時刻依據以下式子更新:

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

3.3 模型驗證

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

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

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

4 數值結果與分析

4.1 網格無關性驗證

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

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

4.2 達姆科勒數的影響

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

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

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

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

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

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

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

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

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

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

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

為了定量描述驅替效率,定義驅替效率

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

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

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

綜上所述,當流場內未生成明顯蟲洞時,Da數越大,驅替效率越小;當流場內生成明顯蟲洞時,隨著時間發展,Da數越大,驅替效率會越大,但是仍遠小于無反應時的驅替效率.

4.3 佩克萊數的影響

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

圖8 t=0.2441時不同Pe 數下的濃度場瞬時圖 (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 為相應情形下多孔介質內水平方向上的平均濃度圖.可以看出,Pe=5 時,平均濃度的變化十分平滑.隨著Pe數的增大(Pe=65),濃度曲線有了細微的波動,說明此時已有指進生成.當Pe數為524 時,濃度曲線上出現了明顯的波動,說明高Pe數下驅替過程中產生大量指進,流體界面不穩定性增加.

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

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

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

4.4 多孔介質結構的影響

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

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

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

5 結論

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

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

2)隨著Da數的增大,反應速率變大,孔隙率變化越大,所生成的蟲洞越寬,驅替效率逐漸增大.

3)隨著Pe數的增大,指進增長的速率越快,反應的量越小,孔隙率變化越小,驅替效率也隨之變小.

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

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