多孔介質(zhì)內(nèi)氣泡Ostwald 熟化特性三維孔網(wǎng)數(shù)值模擬*

張沐安 王進卿? 吳睿 馮致 詹明秀 徐旭 池作和

1) (中國計量大學計量測試工程學院,杭州 310018)

2) (上海交通大學機械與動力工程學院,上海 200240)

多孔介質(zhì)內(nèi)氣泡的Ostwald 熟化行為廣泛存在于CO2 地質(zhì)封存、多孔材料制備、燃料電池等領域.為探究孔隙尺度下多孔介質(zhì)內(nèi)氣泡的熟化特性,建立了基于濃度耦合計算的三維孔隙網(wǎng)絡模型,該模型考慮了氣泡形態(tài)、多孔介質(zhì)結(jié)構以及氣液之間的傳質(zhì),通過求解三維孔隙網(wǎng)絡內(nèi)各孔體的氣相濃度,得到各氣泡的演化過程,并采用四孔隙結(jié)構微流體芯片可視化實驗驗證了模型的可靠性.為分析多孔介質(zhì)非均質(zhì)性對氣泡熟化過程影響,構建了兩種不同孔隙尺寸的三維孔網(wǎng)結(jié)構,對兩個區(qū)域內(nèi)的氣泡熟化過程進行了數(shù)值模擬.結(jié)果表明: 氣泡初始分布會對熟化過程產(chǎn)生影響,當氣泡非均勻分布時,氣泡從小孔隙區(qū)域傳輸至大孔隙區(qū)域的同時,也會各自向自身區(qū)域的大氣泡區(qū)域傳質(zhì);氣泡初始尺寸的差異,會加速熟化進程,使熟化時間明顯短于均勻分布狀態(tài).孔隙數(shù)的選取對平均毛細力、飽和度等連續(xù)尺度等效參數(shù)具有顯著影響,孔隙數(shù)增加時毛細力與飽和度呈現(xiàn)更具規(guī)律性的非線性變化.該模型的建立可以預測地質(zhì)封存過程中CO2 的演化過程,為CO2長期封存過程中非均質(zhì)性的影響機制研究提供指導.

1 引言

氣泡Ostwald 熟化效應是一種在氣液兩相流體系中,受氣泡間不同毛細壓力的影響,小氣泡通過潤濕相分子擴散將自身質(zhì)量傳遞給附近大氣泡的現(xiàn)象.毛細力決定氣泡間濃度分布,濃度梯度驅(qū)動的擴散效應是傳質(zhì)的主導機制.該過程普遍存在于自然界及能源、化工、生化和環(huán)境保護等領域,如超疏水表面自清潔技術[1,2]、油水分離工藝[3]、巖漿脫氣[4,5]、多孔材料制備[6,7]和CO2地質(zhì)封存[8,9]等.最典型的是在CO2地質(zhì)封存中,這種熟化效應會引起CO2氣體在局部區(qū)域的重新聚集,增加了地質(zhì)儲存泄漏的風險[10,11].因此,對多孔介質(zhì)內(nèi)氣泡熟化特性的研究至關重要.

Ostwald 熟化效應可根據(jù)空間狀態(tài)分為自由流體空間與受限空間(多孔介質(zhì))熟化兩種形式.早在20 世紀眾多國內(nèi)外學者已經(jīng)對氣泡在溶液中的熟化特性進行研究,提出并改進了氣泡簇半徑與時間的穩(wěn)態(tài)生長理論[12-15].在自由流體中,氣泡熟化可以不受空間限制,呈現(xiàn)標準的Ostwald 熟化現(xiàn)象.而在多孔介質(zhì)中,空間結(jié)構限制會增加Ostwald熟化過程的復雜性.目前,國內(nèi)外學者已對多孔介質(zhì)環(huán)境下的熟化效應進行了廣泛的實驗研究.Xu 等[16]在實驗過程中成功捕捉到小氣泡與大氣泡的熟化過程,驗證了CO2地質(zhì)封存過程中Ostwald 熟化機制的存在;Xu 等[17]通過2.5D 玻璃微流體芯片實驗發(fā)現(xiàn)氣泡在均質(zhì)多孔介質(zhì)內(nèi)的熟化受限于孔喉的幾何結(jié)構,最終氣泡尺寸趨于一致,發(fā)生反Ostwald 熟化效應.Ostwald 熟化是長時間尺度的過程,在毫米尺度空間域,CO2達到毛細平衡需要數(shù)月至數(shù)年的時間,而對于米級系統(tǒng),熟化則需要一萬年甚至更久[18].可見實驗僅能對短時間內(nèi)的熟化效應進行研究,采用數(shù)值模擬方法對長時間尺度氣泡熟化過程進行預測是有效手段.

在Ostwald 熟化數(shù)值模型的建立方面,相關的研究工作主要集中在孔網(wǎng)模型方面.孔網(wǎng)模型的核心思想是將孔隙結(jié)構簡化為孔隙和孔喉的組成[19-21],從而對其中內(nèi)部平衡問題和傳質(zhì)計算進行簡化.與直接模擬法(DNS)相比具有計算范圍廣、計算效率高的優(yōu)點,能夠較好地反映多孔介質(zhì)內(nèi)孔隙結(jié)構對擴散過程的影響,揭示傳質(zhì)過程孔隙尺度機理.Xu 等[22]通過雙氣泡孔網(wǎng)模型證明重力會誘導CO2向上遷移影響多孔介質(zhì)內(nèi)Ostwald 熟化,并放大孔隙尺度得到一維連續(xù)介質(zhì)模型以評估CO2地質(zhì)封存安全性;Xu 等[23]基于實驗建立了能準確預測氣泡演變的均質(zhì)孔網(wǎng)模型,并指出多孔介質(zhì)內(nèi)孔喉的結(jié)構會影響Ostwald 熟化方向;Mehmani 等[24]通過孔網(wǎng)模型模擬混相氣泡的熟化過程,提出一種完全隱式算法求解多孔介質(zhì)內(nèi)被捕獲氣泡的演化過程,并開發(fā)了一種預測氣泡穩(wěn)定平衡分布狀態(tài)理論,證明氣泡在非均質(zhì)多孔介質(zhì)中的平衡比均勻介質(zhì)中的平衡要慢得多.以上文獻主要在一維和二維條件下開展數(shù)值研究,而自然界和工業(yè)上的問題多為三維多孔介質(zhì).Chalendar 等[25]建立了一種三維錐形孔喉數(shù)值模型,模擬地質(zhì)封存時CO2氣泡演化,通過算法尋找每個有效擴散路徑并對每條路徑的傳質(zhì)進行疊加,實現(xiàn)了對簡單連續(xù)氣泡群演化時間尺度的數(shù)值模擬.但該模型隨著熟化過程中氣泡數(shù)量減少,可搜索的有效擴散路徑數(shù)量大幅增加,算法復雜度將呈指數(shù)增長.

基于質(zhì)量守恒與氣液傳質(zhì)方程,本文構建一種通過耦合計算每個孔隙中氣相濃度場來獲得氣泡生長速率的孔隙網(wǎng)絡數(shù)值模型.該模型考慮氣泡形態(tài)、多孔介質(zhì)結(jié)構以及氣液之間的傳質(zhì),并通過四孔隙結(jié)構微流體芯片可視化實驗驗證模型的可靠性.最后以CO2地質(zhì)封存為研究背景,采用該模型對包含兩種不同孔隙尺度多孔介質(zhì)內(nèi)的氣泡熟化過程進行三維數(shù)值模擬,分析氣泡初始分布和孔隙數(shù)對熟化的影響特性,并探究孔隙尺度平均毛細力與飽和度的關系,為后續(xù)的連續(xù)尺度熟化機制研究提供指導.

2 數(shù)值計算模型

2.1 基本方程

氣泡熟化過程主要受3 個方程控制: 拉普拉斯方程、亨利定律和菲克定律:

式中,pn為氣泡內(nèi)部毛細壓力,ps為溶液壓力,Ra和Rb分別為氣泡的平面和層面曲率半徑,C為溶質(zhì)的濃度,A為兩氣泡間擴散面積,L為兩氣泡間擴散距離.本文計算時選用的擴散面積與擴散距離如圖1 所示,其中A1為氣泡1 與喉道切線對應的弧面面積,A2為喉道截面積,A3為氣泡2 與喉道切線對應的弧面面積;L1為氣泡1 到孔體的距離,L2為喉道長度,L3為氣泡2 到孔體的距離.因此,擴散面積A與擴散距離L之比為

圖1 雙氣泡系統(tǒng)傳質(zhì)過程Fig.1.Two-bubble system mass transfer process.

2.2 計算模型

2.2.1 模型框架

建立的孔網(wǎng)模型假設熟化過程是準靜態(tài)的,只考慮氣體擴散作用.在整個孔隙網(wǎng)絡中,孔體為圓柱體結(jié)構,喉道為長方體結(jié)構,當氣泡直徑小于孔體深度時,氣泡為球形氣泡,其余情況氣泡為圓柱形氣泡.同時做出以下假設: 氣泡與液相均處于局部熱力學平衡狀態(tài),即溶解與生長都是瞬時的;一個孔體內(nèi)只能存在一個氣泡且不可移動,在整個溶解或生長過程中,氣泡質(zhì)量中心始終與孔體中心重合;當氣泡尺寸達到孔體尺寸后,氣泡即停止生長,不侵入喉道.

如前所述,國內(nèi)外學者已對Ostwald 熟化問題開展了數(shù)值研究,不同于他們建立的數(shù)值模型,本模型通過計算孔體內(nèi)的氣相濃度分布得到氣泡生長速率,使整個計算過程更加簡便,提高了計算效率.如圖2 所示,該模型會依次判斷各孔體內(nèi)是否存在氣泡以及該孔體周圍充滿液相的孔體個數(shù),并根據(jù)孔體內(nèi)是否存在氣泡將(3)式中的濃度計算分為兩種情況: 當孔體內(nèi)為氣泡時,根據(jù)拉普拉斯方程和亨利定律求解氣泡濃度;當孔體內(nèi)為液相時,基于質(zhì)量守恒方程建立線性方程組求解孔體內(nèi)氣相濃度.以下對兩種情況的濃度計算進行介紹.

圖2 計算程序算法流程圖Fig.2.Flow chart for the calculation algorithm.

2.2.2 孔體內(nèi)為氣泡

如圖3(a)所示,孔體i內(nèi)存在氣泡i,其相鄰六個方向的孔體j內(nèi)均存在氣泡j.根據(jù)這些氣泡的初始半徑,代入(1)式和(2)式可以分別得到i和j孔體中氣相濃度Ci與Cj,即

圖3 數(shù)值計算模型 (a)孔體i 內(nèi)為氣泡;(b)孔體i 內(nèi)為液相Fig.3.Numerical calculation model: (a) Pore body i with a bubble;(b) pore body i filled with liquid phase.

2.2.3 孔體內(nèi)為液相

當孔體i內(nèi)充滿液相時,無法依據(jù)拉普拉斯方程與亨利定律求解氣相濃度Ci,需基于質(zhì)量守恒與氣液傳質(zhì)方程建立線性方程組求解.構建的孔網(wǎng)模型考慮了孔體連接方式的結(jié)構化,即每個孔體周圍都有6 個孔隙連接.如圖3(b)所示,孔體i中充滿液相,相鄰6 個方向的孔體j內(nèi)均存在氣泡j,根據(jù)氣泡j的初始半徑,代入(5)式可分別得到相鄰氣泡的氣相濃度Cj,則充滿液相的孔體i中氣相濃度Ci可基于傳質(zhì)方程(6)求解:

在整個三維孔網(wǎng)模型中,孔體內(nèi)氣泡與液相隨機分布,因此對于孔體i內(nèi)充滿液相的情況,需要依次判斷周圍孔體j是否充滿液相,建立關于濃度的線性方程組,并采用BiCGSTAB 方法[26]求解出充滿液相的孔體內(nèi)氣相濃度.根據(jù)各孔體的氣相濃度,可以得到氣泡的生長速率.氣泡i通過熟化效應在單位時間變化的質(zhì)量為該氣泡與其相鄰氣泡傳遞質(zhì)量的和,見(7)式.當熟化系統(tǒng)中所有孔隙內(nèi)的氣相濃度相等時,熟化結(jié)束,停止計算.

2.2.4 時間步長

在該程序中,一個時間步長為整個孔隙空間內(nèi)所有氣泡中半徑增大或減小指定長度R'所需要的最短時間.因此,時間步長會隨著熟化程序的迭代進行自動優(yōu)化.在熟化初期,氣泡間尺寸半徑差異大,熟化速率快,時間步長短;而熟化后期,熟化速率減慢,時間步長增大.本文計算模型選取的R'值為0.1 μm.

(8)式為適用于球形與圓柱形氣泡的體積公式,計算程序中單次輸出的熟化時間步長為

2.3 計算對象

在實際多孔介質(zhì)環(huán)境中,多孔介質(zhì)結(jié)構通常呈現(xiàn)非均質(zhì)性,為了探究孔隙結(jié)構非均質(zhì)性對熟化過程的影響,選取了由兩種不同孔隙尺度區(qū)域構成的三維孔網(wǎng)結(jié)構作為計算對象,對兩個區(qū)域內(nèi)的氣泡熟化過程進行數(shù)值模擬.如圖4 所示,小孔隙區(qū)域、大孔隙區(qū)域孔體半徑分別為20 和50 μm,兩個區(qū)域深度均為30 μm;所有喉道半徑均為5 μm,深度為5 μm.

圖4 多孔隙結(jié)構熟化模型Fig.4.Multi bubble ripening model.

在孔隙尺度,氣泡群初始分布會顯著影響熟化過程,氣泡群非均勻分布時的熟化規(guī)律相較均勻分布有明顯差別[23],因此選取了兩種不同氣泡群初始分布工況(工況1#、2#)進行分析(表1).與工況1#相比,工況2#僅改變第3 層氣泡群初始半徑.同時,孔隙數(shù)的選取對平均毛細力、飽和度等連續(xù)尺度等效參數(shù)的獲取具有顯著影響,為了分析孔隙數(shù)對熟化過程的影響,本文構建了孔隙數(shù)分別為325 個(工況1#)和10800 個(工況3#)兩種多孔介質(zhì),具體參數(shù)見表1.工況3#氣泡群初始分布與工況1#相同.模擬在常溫常壓條件下進行,表2 整理了氣泡熟化模擬時選用的參數(shù).

表1 模擬工況參數(shù)Table 1.Simulation parameters of three conditions.

表2 四孔隙及多孔隙結(jié)構氣泡熟化模擬參數(shù)Table 2.Simulation parameters of four bubble and muti bubble ripening system.

3 實驗方法

本實驗選用硅基微流體芯片,以去離子水作為液相,CO2作為氣相.芯片分別對硅基底進行兩次刻蝕制備出孔隙與喉道圖案,得到具有不同深度流道的2.5D 結(jié)構,能對芯片內(nèi)的氣泡進行夾斷,使每個孔隙里只有一個氣泡[17].芯片觀察部分為四孔隙結(jié)構,具體結(jié)構如圖5 所示,其中孔隙深度為30 μm,喉道深度為10 μm.本文采用氣泡熟化可視化實驗裝置對氣泡在四孔隙結(jié)構內(nèi)熟化過程進行研究,具體的實驗裝置與方法見參考文獻[29].

圖5 微流控芯片內(nèi)部結(jié)構Fig.5.Internal structure of microfluidic chip.

4 討論部分

4.1 模型驗證

對四孔隙結(jié)構氣泡熟化系統(tǒng)進行實驗及模擬研究,結(jié)果如圖6 所示.熟化過程中,雖然氣泡1 的半徑小于氣泡3,但熟化結(jié)束后氣泡3 消失而氣泡1 并未消失.這是由于氣泡1 受到了氣泡2 的質(zhì)量傳輸,而氣泡3 周圍存在半徑更大的氣泡4,使得氣泡3 向氣泡4 傳質(zhì),最終導致上述現(xiàn)象的發(fā)生.可見氣泡的初始分布不同,會導致不同熟化過程的發(fā)生.

圖6 四孔隙結(jié)構氣泡熟化過程 (a)實驗過程;(b)模擬過程Fig.6.Four bubble ripening system: (a) Experiments process;(b) simulations process.

圖7 為四孔隙結(jié)構熟化過程中氣泡半徑變化的模擬和實驗結(jié)果對比,模擬和實驗結(jié)果符合較好.在熟化時間上,除氣泡3 存在差別,其余氣泡達到穩(wěn)定的模擬時間與實驗較為吻合.造成這種差別的原因在于模擬過程中氣泡位置始終位于孔體正中心,而在實際實驗中,氣泡的位置會與孔體正中心存在一定偏差.隨著熟化的進行,模擬和實驗結(jié)果偏差會隨之變大.

圖7 四孔隙結(jié)構氣泡半徑隨時間變化Fig.7.Evolution of curvature radius in four bubble ripening system.

4.2 氣泡初始分布對熟化過程影響

由于工況1#各層孔隙結(jié)構和氣泡初始半徑均相同,因此三維空間內(nèi)不同深度方向之間不存在傳質(zhì),故選用單層二維平面數(shù)據(jù)進行討論.工況1#中氣泡熟化過程如圖8 所示.觀察圖8 發(fā)現(xiàn):1)由于小孔隙區(qū)域氣泡半徑小,毛細壓力大,因此在熟化初期小孔隙區(qū)域氣泡半徑顯著減小,氣體從該區(qū)域傳輸至大孔隙區(qū)域;2)熟化過程中,大、小孔隙區(qū)域氣泡的消失順序呈現(xiàn)間隔性,如熟化進行至14 h 時,小孔隙區(qū)域殘留氣泡呈現(xiàn)規(guī)則性的間隔分布,4.4×1018h 的大孔隙區(qū)域殘留氣泡分布也呈現(xiàn)相同規(guī)律;3)在熟化末期,小孔區(qū)域殘留氣泡集中在中間區(qū)域,即靠近及遠離大孔隙區(qū)域的氣泡最先消失,該現(xiàn)象對連續(xù)尺度模型的建立具有指導意義,目前文獻認為小孔隙區(qū)域氣泡消失的順序應從靠近大孔隙區(qū)域開始逐列向遠離該區(qū)域的方向推進.

圖8 工況1#模擬結(jié)果 (a)氣泡演化過程,其中白色區(qū)域為CO2 氣泡,藍色區(qū)域為去離子水,灰色區(qū)域為硅顆粒;(b) CO2 氣泡曲率半徑變化圖Fig.8.Simulation results of condition 1#: (a) Bubble evolution process,where the white regions are CO2 bubble,blue regions are filled with DI water,and the gray regions are silicon grains;(b) curvature radius variation diagram of each CO2 bubble.

為了解釋上述現(xiàn)象,采用一維孔隙結(jié)構的氣泡熟化過程進行說明(圖9).一維孔隙結(jié)構由3 個小孔體和1 個大孔體構成,每個孔體間通過喉道相連.小孔體內(nèi)的氣泡大小相同,且小于大孔體內(nèi)氣泡.氣泡間傳質(zhì)受毛細力決定的濃度梯度與擴散距離影響.因此,熟化開始時,與大氣泡4 毗鄰的小氣泡3 最先發(fā)生熟化,導致氣泡3 半徑減小,此時氣泡4、氣泡3 和氣泡2 的半徑大小關系為R4＞R2＞R3,故氣泡3 同時向氣泡4 和氣泡2 傳質(zhì)(見圖9 箭頭方向).氣泡2 獲得傳質(zhì)后半徑增大,使得R2＞R1,引起氣泡1 向氣泡2 傳質(zhì).上述傳質(zhì)現(xiàn)象的存在,使得氣泡發(fā)生規(guī)則性的間隔消失.在多孔隙結(jié)構(圖8(a))中,這種間隔變化使小孔隙區(qū)域內(nèi)氣泡質(zhì)量往大孔隙區(qū)域傳輸?shù)耐瑫r,還會往中間區(qū)域積累,因此在小孔隙區(qū)域內(nèi)位于中間列的氣泡最慢消失.

圖9 一維孔隙結(jié)構氣泡熟化過程Fig.9.Bubble evolution process of one-dimensional structure.

工況2#中,由于孔體結(jié)構和氣泡初始尺寸的對稱性,第1 層和第5 層的熟化情況相同,第2 層與第4 層的熟化情況相同.工況2#中氣泡熟化過程如圖10 所示.觀察圖10 發(fā)現(xiàn),由于氣泡非均勻分布,熟化過程與工況1#存在以下區(qū)別: 1)工況1#中由于氣泡均勻分布,不存在向不同深度的中間層區(qū)域傳質(zhì)的現(xiàn)象,而工況2#中氣泡從小孔隙區(qū)域傳輸至大孔隙區(qū)域的同時,大、小孔隙區(qū)域也會各自向自身區(qū)域的大氣泡區(qū)域傳質(zhì),因此整個熟化過程會向中間層區(qū)域逐層推進;2)如圖10(b),(d)所示,由于氣泡初始分布的差異,小孔隙區(qū)域內(nèi)氣泡熟化進行更快,整個孔隙空間內(nèi)所有氣泡完成熟化的時間為1.7×1018h,與工況1#中完成熟化的時間量級上較為相近.進一步證明了氣泡熟化是一個長時間尺度的過程,CO2地質(zhì)封存中達到平衡需要數(shù)百萬年時間[30],時間尺度的研究可以確定封存過程中CO2的長期儲存能力,評估地質(zhì)封存系統(tǒng)的穩(wěn)定性和長期效果;3)對比工況1# (圖8(b))與工況2# (圖10(b),(d),(f))中氣泡曲率半徑變化發(fā)現(xiàn),由于工況2#中自身熟化現(xiàn)象更顯著,在熟化末期大孔隙區(qū)域內(nèi)大部分氣泡都會充滿孔隙.該模型可以模擬CO2氣泡在地下儲層中的傳輸和分布過程,確定潛在的CO2泄漏路徑和風險,為CO2地質(zhì)封存技術提供重要的指導意義.

圖10 工況2#模擬結(jié)果 (a)第1 層和第5 層氣泡演化過程;(b)第1 層和第5 層CO2 氣泡曲率半徑變化;(c)第2 層和第4 層氣泡演化過程;(d)第2 層和第4 層CO2 氣泡曲率半徑變化;(e)第3 層氣泡演化過程;(f)第3 層CO2 氣泡曲率半徑變化Fig.10.Simulation results of condition 2#: (a) Bubble evolution process in layers 1 and 5;(b) curvature radius variation of each CO2 bubble in layers 1 and 5;(c) bubble evolution process in layers 2 and 4;(d) curvature radius variation of each CO2 bubble in layers 2 and 4;(e) bubble evolution process in layer 3;(f) curvature radius variation of each CO2 bubble in layer 3.

4.3 孔隙數(shù)對熟化過程影響

熟化過程連續(xù)尺度模型參數(shù)的構建依賴于孔隙尺度熟化特性的研究,可基于孔隙尺度模擬獲得連續(xù)尺度等效參數(shù),如平均毛細力、飽和度等參數(shù).孔隙尺度下多孔介質(zhì)孔隙數(shù)的選取對連續(xù)尺度等效參數(shù)的獲取具有顯著影響.為此,本文對工況1#和工況3#熟化過程中CO2氣泡受到的毛細力與飽和度變化進行統(tǒng)計.熟化過程中氣泡受到的平均毛細力如圖11(a)所示,當孔隙數(shù)較少時,個別氣泡半徑的突變也會對整個孔隙空間產(chǎn)生較大的影響,使毛細力曲線走勢振蕩;當孔隙數(shù)較多時,毛細力呈下降趨勢,曲線波動較少.熟化過程中飽和度隨時間的變化如圖11(b)所示,熟化初期小孔隙區(qū)域內(nèi)的氣泡向大孔隙區(qū)域傳質(zhì),導致小孔隙區(qū)域飽和度持續(xù)下降,而在熟化后期大孔隙區(qū)域內(nèi)的氣泡之間也會發(fā)生熟化使大孔隙區(qū)域整體飽和度降低.孔隙數(shù)的增加提供了更多的氣體擴散路徑,使大孔隙區(qū)域內(nèi)氣泡自身熟化的情況更復雜.

圖11 CO2 氣泡各參數(shù)變化圖 (a)毛細力;(b) CO2 飽和度Fig.11.Variation diagram of CO2 bubble parameters: (a) Capillary pressure;(b) CO2 saturation.

在構建連續(xù)尺度熟化模型時,毛細力與飽和度關系式的選取直接影響熟化模型精度[31].目前國內(nèi)外學者在建立多孔介質(zhì)內(nèi)CO2熟化連續(xù)尺度模型時,認為熟化過程中CO2氣泡毛細力與飽和度呈線性關系,如Li 等[18]提出的連續(xù)尺度數(shù)值模型中選用了“Pc=a+bSg”的毛細力模型.本文分別對兩種工況模擬過程中氣泡毛細力與飽和度進行統(tǒng)計,其變化關系如圖12 所示,兩種工況下大、小孔隙區(qū)域內(nèi)的毛細力與飽和度的變化均呈現(xiàn)非線性關系.對比圖12(a)和圖12(b)發(fā)現(xiàn),當孔隙數(shù)增加時,大孔隙區(qū)域和小孔隙區(qū)域內(nèi)毛細力與飽和度關系的規(guī)律性顯著提升.這是由于熟化過程中氣泡的變化情況非常復雜,氣泡所受的毛細力與飽和度都經(jīng)歷了大幅改變.工況3#中毛細力與飽和度關系的變化具有更明顯的規(guī)律性,故對該工況下的數(shù)據(jù)點進行非線性擬合,得到小孔隙區(qū)域內(nèi)毛細力與飽和度的關系式為大孔隙區(qū)域內(nèi)毛細力與飽和度的關系式為Pc=.研究毛細力和飽和度的關系,可以提升連續(xù)尺度熟化模型的可靠性和預測能力,為CO2地質(zhì)封存過程提供指導和決策支持.

圖12 氣泡毛細力與CO2 飽和度關系對比圖 (a)小孔隙區(qū)域;(b)大孔隙區(qū)域Fig.12.Comparison of capillary pressure with CO2 saturation curves: (a) Small pore region;(b) large pore region.

5 結(jié)論

本文構建了基于濃度耦合計算的三維孔隙網(wǎng)絡模型,并采用四孔可視化實驗對模型準確度進行了驗證,隨后對三維條件下不同初始分布和不同孔隙數(shù)的氣泡開展數(shù)值模擬研究,結(jié)論如下.

1)本文建立的模型考慮了氣泡形態(tài)、多孔介質(zhì)結(jié)構以及氣液之間的傳質(zhì),通過求解三維孔網(wǎng)內(nèi)各孔體的氣相濃度,可以得到各氣泡的演化過程.該模型數(shù)值模擬結(jié)果與四孔隙結(jié)構實驗吻合較好.

2)數(shù)值模擬結(jié)果表明氣泡初始分布和孔隙數(shù)均會對熟化過程產(chǎn)生影響,當氣泡非均勻分布時,由于大、小孔隙區(qū)域會各自向自身區(qū)域的大氣泡傳質(zhì),熟化過程中殘留氣泡不再呈現(xiàn)規(guī)律性的間隔分布,并且熟化進程明顯加速;孔隙數(shù)增加時,大、小孔隙區(qū)域的毛細力與飽和度均呈現(xiàn)更具規(guī)律性的非線性變化.

3)本文得到大、小孔隙區(qū)域毛細力與飽和度的關系式分別為:,與現(xiàn)有文獻的假設有所不同,該結(jié)果對連續(xù)尺度熟化模型的構建具有重要指導意義.但該關系式并未考慮不同氣泡初始分布、孔隙數(shù)、孔喉比等參數(shù)的影響,后續(xù)將在這些方面開展孔隙尺度模擬.