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

張沐安 王進(jìn)卿? 吳睿 馮致 詹明秀 徐旭 池作和

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

2) (上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院,上海 200240)

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

1 引言

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

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

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

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

2 數(shù)值計(jì)算模型

2.1 基本方程

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

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

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

2.2 計(jì)算模型

2.2.1 模型框架

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

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

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

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

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

圖3 數(shù)值計(jì)算模型 (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)為液相

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

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

2.2.4 時(shí)間步長

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

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

2.3 計(jì)算對象

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

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

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

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

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

3 實(shí)驗(yàn)方法

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

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

4 討論部分

4.1 模型驗(yàn)證

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

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

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

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

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

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

圖8 工況1#模擬結(jié)果 (a)氣泡演化過程,其中白色區(qū)域?yàn)镃O2 氣泡,藍(lán)色區(qū)域?yàn)槿ルx子水,灰色區(qū)域?yàn)楣桀w粒;(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é)構(gòu)的氣泡熟化過程進(jìn)行說明(圖9).一維孔隙結(jié)構(gòu)由3 個(gè)小孔體和1 個(gè)大孔體構(gòu)成,每個(gè)孔體間通過喉道相連.小孔體內(nèi)的氣泡大小相同,且小于大孔體內(nèi)氣泡.氣泡間傳質(zhì)受毛細(xì)力決定的濃度梯度與擴(kuò)散距離影響.因此,熟化開始時(shí),與大氣泡4 毗鄰的小氣泡3 最先發(fā)生熟化,導(dǎo)致氣泡3 半徑減小,此時(shí)氣泡4、氣泡3 和氣泡2 的半徑大小關(guān)系為R4＞R2＞R3,故氣泡3 同時(shí)向氣泡4 和氣泡2 傳質(zhì)(見圖9 箭頭方向).氣泡2 獲得傳質(zhì)后半徑增大,使得R2＞R1,引起氣泡1 向氣泡2 傳質(zhì).上述傳質(zhì)現(xiàn)象的存在,使得氣泡發(fā)生規(guī)則性的間隔消失.在多孔隙結(jié)構(gòu)(圖8(a))中,這種間隔變化使小孔隙區(qū)域內(nèi)氣泡質(zhì)量往大孔隙區(qū)域傳輸?shù)耐瑫r(shí),還會往中間區(qū)域積累,因此在小孔隙區(qū)域內(nèi)位于中間列的氣泡最慢消失.

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

工況2#中,由于孔體結(jié)構(gòu)和氣泡初始尺寸的對稱性,第1 層和第5 層的熟化情況相同,第2 層與第4 層的熟化情況相同.工況2#中氣泡熟化過程如圖10 所示.觀察圖10 發(fā)現(xiàn),由于氣泡非均勻分布,熟化過程與工況1#存在以下區(qū)別: 1)工況1#中由于氣泡均勻分布,不存在向不同深度的中間層區(qū)域傳質(zhì)的現(xiàn)象,而工況2#中氣泡從小孔隙區(qū)域傳輸至大孔隙區(qū)域的同時(shí),大、小孔隙區(qū)域也會各自向自身區(qū)域的大氣泡區(qū)域傳質(zhì),因此整個(gè)熟化過程會向中間層區(qū)域逐層推進(jìn);2)如圖10(b),(d)所示,由于氣泡初始分布的差異,小孔隙區(qū)域內(nèi)氣泡熟化進(jìn)行更快,整個(gè)孔隙空間內(nèi)所有氣泡完成熟化的時(shí)間為1.7×1018h,與工況1#中完成熟化的時(shí)間量級上較為相近.進(jìn)一步證明了氣泡熟化是一個(gè)長時(shí)間尺度的過程,CO2地質(zhì)封存中達(dá)到平衡需要數(shù)百萬年時(shí)間[30],時(shí)間尺度的研究可以確定封存過程中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泄漏路徑和風(fēng)險(xiǎn),為CO2地質(zhì)封存技術(shù)提供重要的指導(dǎo)意義.

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

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

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

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

5 結(jié)論

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

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

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

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