李延銘，吳 睿，崔國民，程 樹

（上海理工大學 能源與動力工程學院新能源所,上海 200093）

環(huán)境與化工

非均質多孔介質內超臨界二氧化碳驅替鹽水過程的孔隙網(wǎng)絡模擬

李延銘，吳 睿，崔國民，程 樹

（上海理工大學 能源與動力工程學院新能源所,上海 200093）

針對CO2在非均勻多孔介質內的儲存過程,建立二維動態(tài)孔隙網(wǎng)絡模型,對含鹽水的介質中超臨界CO2驅替鹽水的過程進行模擬計算。提出用非均勻性來衡量二維動態(tài)孔隙網(wǎng)絡模型中喉道的分布特性,通過建立3種非均勻性不同的二維動態(tài)孔隙網(wǎng)絡模型,分析非均勻性對CO2封存的影響。模擬計算結果表明,非均勻性越大的二維動態(tài)孔隙網(wǎng)絡模型,驅替過程中超臨界CO2的前進端越不穩(wěn)定,且最后突破時的超臨界CO2飽和度越低；通過分析穩(wěn)定驅替下不同時刻的CO2分布,揭示非均勻性的影響機理；非均勻性越大,陷住的含鹽水孔隙更多是導致突破時超臨界CO2飽和度越低的主要原因。

二氧化碳封存；二維動態(tài)孔隙網(wǎng)絡模型；超臨界二氧化碳；非均勻性

CO2捕集與儲存技術（CCS）是當今減緩大氣中CO2濃度過高和溫室效應的有效措施［1］。CCS是指將大型CO2排放源排放出的CO2收集并進行長時間儲存,使其與大氣相隔絕的一種技術。地下深層多孔的地質結構,例如廢棄的油氣田和沉積含鹽水地層,可以為CO2提供足夠的儲存空間。沉積含鹽水地層是全球分布范圍較廣的一種理想儲存空間,其中的部分砂層、砂巖、石灰?guī)r具有較高的孔隙度和滲透率,可以為CO2儲存提供儲存空間,而頁巖和泥巖雖然有較高的孔隙度,但滲透率很低,這種低滲透率的巖層稱之為巖蓋［2］,可以很好地對CO2實現(xiàn)隔離。

實際操作中,CO2封存是一個非常復雜的過程,其封存過程主要受到黏性力、浮力和毛細壓力的影響［3-5］。封存在地下的CO2部分溶解可能會導致地下水酸化形成離子沉淀［6］,造成地質結構發(fā)生改變［7］。因此地質結構的孔隙度、滲透率、穩(wěn)定性和地下水鹽度等都是CO2封存所必須考慮的因素。

高效安全地實現(xiàn)CO2封存是目前研究的重點,這方面的研究主要分為可視化實驗［8-12］以及數(shù)值計算［13-16］兩大部分。目前,在數(shù)值計算方面主要采用建立孔隙網(wǎng)絡模型研究CO2的封存特性?？紫毒W(wǎng)絡模型從微觀上揭示了多孔介質中兩相流動的原理,為了解巖石微觀結構對CO2封存的影響提供了理論依據(jù)。

圖1 二維動態(tài)孔隙網(wǎng)絡模型Fig.1 Two-dimensional dynamic pore network model.

本工作建立二維動態(tài)孔隙網(wǎng)絡模型,模擬地質條件下儲存CO2的過程,分析不同的孔隙結構對CO2驅替鹽水過程的影響。

1 模型描述

將無規(guī)則的巖石微觀結構抽象為理想的幾何形狀,建立二維動態(tài)孔隙網(wǎng)絡模型（見圖1）。二維動態(tài)孔隙網(wǎng)絡模型用來描述一個含潤濕相（原流體,如鹽水）和非潤濕相（侵入流體,如超臨界CO2）的正方形晶格。由圖1可見,二維動態(tài)孔隙網(wǎng)絡模型由圓柱形喉道和圓形孔隙構成,每個孔隙有4條喉道連接形成格子型網(wǎng)絡。其中,喉道代表實際介質中較小的空間,孔隙代表較大的空間。二維動態(tài)孔隙網(wǎng)絡模型中孔隙半徑（rpore）和喉道長度（Lthroat）均為定值,喉道半徑大小介于rmin和rmax之間,非均勻性用式（1）表示。

2 計算過程

二維動態(tài)孔隙網(wǎng)絡模型所有孔隙中一開始充滿原流體,侵入流體自下向上侵入。邊界條件為：二維動態(tài)孔隙網(wǎng)絡模型左右兩側為封閉的邊界；下端為體積流量恒定,其大小給定；上端為壓力恒定,其大小為0。

為了簡化計算過程,對二維動態(tài)孔隙網(wǎng)絡模型及流體流動過程做以下假設：1） 超臨界CO2與固體表面間的接觸角為定值（165°）；2） 雖然喉道決定導流率的大小,仍認為喉道的體積遠小于孔隙的體積；3） 對于每一個喉道,只有一種流體可以滯留；4） 喉道中流體的流動狀態(tài)假設為層流并遵從哈根-泊肅葉定律；5） 忽略孔隙的流動阻力；6）兩種流體均為不可壓縮流體。

整個計算過程可分為4步：1）壓力場及流場的預處理；2）壓力場及流場的求解；3）時間步長的確定；4）彎液面的更新。驅替過程從下邊界有侵入流體流入開始,到上邊界有侵入流體流出（突破）為止。

兩相流體間的接觸面稱為彎液面（見圖1）,其穩(wěn)定性由式（2）確定。

若孔隙中充滿侵入流體且Δp≤0,則彎液面是不可流動的,即處于穩(wěn)定狀態(tài)；若Δp＞0,則彎液面是可流動的,即處于不穩(wěn)定狀態(tài)；若孔隙只有一部分被侵入流體占據(jù),由于孔隙中流動阻力大小不計,所以此時彎液面仍處于不穩(wěn)定狀態(tài)。

首先,根據(jù)式（2a）判斷彎液面的穩(wěn)定性。若有不穩(wěn)定的彎液面,則不需要預處理,若所有彎液面均為穩(wěn)定,則進行預處理。當所有彎液面均處于穩(wěn)定狀態(tài)時,為了保證驅替過程的連續(xù),假設最大孔道內的彎液面處于不穩(wěn)定狀態(tài),且侵入流體占據(jù)此孔道。侵入流體在喉道內的流量由式（3）計算。

式（3b）的意義為當相連接兩孔隙間的壓力差小于門闕壓力時（Δp≤0）,喉道內侵入流體的流量為0。原流體在喉道內的流量由式（4）計算。

根據(jù)質量守恒定律,并聯(lián)立式（2）～（4）可獲得壓力場和流場。壓力場及流場確定后,需要一個時間步長推動彎液面前進。先求得每個與彎液面相連接的孔隙填滿所需的時間,然后取其中的最小值確定為時間步長,見式（5）。

時間步長確定后,由式（6）迭代更新彎液面。

陷住的孔隙示意圖見圖2。

圖2 陷住的孔隙示意圖Fig.2 Trapped pores.

由圖2可見,若出現(xiàn)侵入流體將原流體包圍的情況,則被困住的原流體不能再移動［17-18］,即孔隙被陷住,陷住的孔隙可能是一個也可能是多個。

3 結果與討論

采用兩個無量綱參數(shù)：毛細數(shù)（Ca）和黏度比（M）,表征流體的流動狀態(tài)。Ca是黏性力與毛細力的比值,見式（7）。M為兩液相黏度的比值,見式（8）。

通過改變超臨界CO2的流速和黏度來改變Ca和M。

二維動態(tài)孔隙網(wǎng)絡模型的參數(shù)見表1。

表1 二維動態(tài)孔隙網(wǎng)絡模型的參數(shù)Table 1 Parameters of the two-dimensional dynamic pore network model

3種二維動態(tài)孔隙網(wǎng)絡模型結構的非均勻性分別為1.8,1.2,0.4,對應的喉道半徑分別為0.5～9.5,2.0～8.0,4.0～6.0 μm。采用的M分別為0.1,1.0；Ca分別為10-1,10-3,10-5。以超臨界CO2飽和度衡量超臨界CO2驅替鹽水的效率。

不同參數(shù)下突破時超臨界CO2在二維動態(tài)孔隙網(wǎng)絡模型中的分布見圖3。

圖3 不同參數(shù)下突破時超臨界CO2在二維動態(tài)孔隙網(wǎng)絡模型中的分布Fig.3 Distribution of supercritical CO2in the two-dimensional dynamic pore network model when breakthrough occurred with different parameters.Ca：capillary number；M：viscosity ratio；χ：heterogeneity. The dark part in the f gure expressed supercritical CO2.

由圖3可知,M較小、Ca較大時（M=0.1；Ca=10-1,10-3）,3種二維動態(tài)孔隙網(wǎng)絡模型結構的驅替過程均呈現(xiàn)出沿前進方向不穩(wěn)定擴散的黏性指進；M較大、Ca較大時（M=1.0；Ca=10-1,10-3）, 3種二維動態(tài)孔隙網(wǎng)絡模型結構的驅替過程均較為穩(wěn)定；Ca較小時（Ca=10-5；M=0.1,1.0）,3種二維動態(tài)孔隙網(wǎng)絡模型結構的驅替過程均呈現(xiàn)出沿任意方向不穩(wěn)定擴散的毛細指進。

黏性指進主要是由于黏性力占優(yōu)時兩種流體黏度差異造成的一種不穩(wěn)定驅替現(xiàn)象,而毛細指進是由于毛細力占優(yōu)時喉道隨機分布造成的一種不穩(wěn)定驅替現(xiàn)象［19-20］。在黏性力占優(yōu)的情況下,兩種流體黏度越接近,侵入流體的前進端越穩(wěn)定。隨非均勻性的增大,陷住的含鹽水孔隙越來越多甚至出現(xiàn)局部大面積陷住的情況,而且超臨界CO2的前進端越不穩(wěn)定,尤其在穩(wěn)定驅替情況下,孔隙由單個到多個、從少量到大量的陷住更為明顯。

不同參數(shù)下非均勻性對超臨界CO2飽和度的影響見圖4。由圖4可知,超臨界CO2飽和度隨非均勻性的增大而減小。

圖4 不同參數(shù)下非均勻性對超臨界CO2飽和度的影響Fig.4 Effect of χ on the supercritical CO2saturation with the different parameters.Ca=10-1,M=0.1；Ca=10-1,M=1.0；Ca=10-3,M=0.1；Ca=10-3, M=1.0；Ca=10-5,M=0.1；Ca=10-5,M=1.0

為了解釋非均勻性導致的超臨界CO2飽和度差異,將穩(wěn)定驅替情況下χ=1.8和χ=0.4時二維動態(tài)孔隙網(wǎng)絡模型抽象為兩種簡單的局部示意圖,見圖5。

圖5 不同二維動態(tài)孔隙網(wǎng)絡模型結構的侵入過程Fig.5 Displacement process in different two-dimensional network models.

由式（2b）可知,超臨界CO2從網(wǎng)絡模型下端選擇最大喉道開始侵入（見圖5（a））；隨驅替過程的進行,為保證流量的恒定,底部壓力會有一定的波動,此時對于χ=0.4的結構若底部壓力超過右側較小喉道的門闕壓力時則較小喉道開始有超臨界CO2侵入,對于χ=1.8的結構較小喉道的門闕壓力過大時仍難以開始侵入（見圖5（b））；隨驅替過程的進一步進行,χ=0.4的結構中超臨界CO2的體積與χ=1.8的結構中超臨界CO2的體積差距增大（見圖5（c））；χ=1.8的結構右下角出現(xiàn)陷住,最終導致在不同的非均勻性下超臨界CO2飽和度存在差異（見圖5（d））。

穩(wěn)定驅替下超臨界CO2前進端的形態(tài)變化見圖6。由圖6（a）可知,驅替過程開始階段,不同的二維動態(tài)孔隙網(wǎng)絡模型結構中超臨界CO2的前進端就開始呈現(xiàn)了不同的現(xiàn)象。非均勻性越大的二維動態(tài)孔隙網(wǎng)絡模型結構中超臨界CO2的前進端越不穩(wěn)定,隨驅替過程的進行,這種前進端穩(wěn)定性的差異越來越明顯,見圖6（b）～（d）。

結合圖5和圖6可知,非均勻性越大的二維動態(tài)孔隙網(wǎng)絡模型結構中陷住的含鹽水孔隙越多,這是導致突破時超臨界CO2飽和度越低的主要原因。

圖6 穩(wěn)定驅替下超臨界CO2前進端的形態(tài)變化Fig.6 Patterns of the CO2front in stable displacement.t*： ratio of total time that current state spent.Parameters：Ca=10-1,M=0.1.

雖然在實際CO2封存過程中孔隙率占主要作用［21］,但通過模擬計算結果發(fā)現(xiàn),地質微觀結構的非均勻性也對CO2封存有重要影響。非均勻性較大的結構孔隙率更高,但封存過程中的超臨界CO2飽和度較低。相反非均勻性較小的結構有較低的孔隙率,但封存過程中的超臨界CO2飽和度較高。因此在實際CO2封存過程中應該考慮地質微觀結構非均勻性的影響。

4 結論

1） 為了解決地質微觀結構內喉道分布非均勻特性對超臨界CO2封存的影響,建立了3種非均勻性不同的二維動態(tài)孔隙網(wǎng)絡模型。通過對二維動態(tài)孔隙網(wǎng)絡模型喉道分布非均勻特性的數(shù)值模擬研究,發(fā)現(xiàn)不同的喉道分布模型對CO2的驅替效率有重要影響。

2）非均勻性越大的孔隙網(wǎng)絡模型中超臨界CO2的前進端越不穩(wěn)定,導致最后突破時的超臨界CO2飽和度越低。基于穩(wěn)定驅替下的分步計算結果,發(fā)現(xiàn)非均勻性較大的二維動態(tài)孔隙網(wǎng)絡模型中門闕壓力分布不均勻導致了超臨界CO2前進端的不穩(wěn)定,而更多陷住的含鹽水孔隙是導致突破時超臨界CO2飽和度低的主要原因。

符 號 說 明

Ca 毛細數(shù)

M 黏度比

pi第i個孔隙的壓力,Pa

pj第j個孔隙的壓力,Pa

pt門闕壓力,Pa

Δp 兩個孔隙之間的壓差,Pa

rmin最小喉道半徑,m

rmax最大喉道半徑,m

rave喉道半徑平均值,m

rpore孔 隙半徑,m

rthroat喉道半徑,m

Δt 時間步長,s

Vi第i個孔隙的體積,m3

v 流體的流速,m/s

γ 潤濕相與非潤濕相的流體界面張力,N/m

θ 超臨界CO2與固體表面的接觸角,°

μnw非潤濕相的流體動力黏度,Pa?s

μw潤濕相的流體動力黏度,Pa?s s/kg

χ 非均勻性

上角標

nw 非潤濕相

w 潤濕相

下角標

i 孔隙的編號

j 孔隙的編號

nw 非潤濕相

［n］ 迭代步數(shù)

w 潤濕相

［1］ Metz B,Davidson O,De Coninck H C,et al. Special Report on Carbon Dioxide Capture and Storage［R］. Working Group III of the Intergovernmental Panel on Climate Change,2005.

［2］ Blunt M. Carbon Dioxide Storage［J］. Imperial College London Grantham Institute Climate Change Briefing Paper,2010（4）：1 - 14.

［3］ Ide T S,Jessen K,Orr Jr F M. Storage of CO2in Saline Aquifers：Effects of Gravity,Viscous,and Capillary Forces on Amount and Timing of Trapping［J］. Int J Greenh Gas Con, 2007,1（4）：481 - 491.

［4］ Bryant S L,Lakshminarasimhan S,Pope G A. Buoyancy-Dominated Multiphase Flow and Its Effect on Geological Sequestration of CO2［J］. SPE J,2008,13（4）：447 - 454.

［5］ Mo S,Akervoll I. Modeling Long-Term CO2Storage in Aquifer with a Black-Oil Reservoir Simulator［C］//SPE/EPA/DOE Exploration and Production Environmental Conference. Society of Petroleum Engineers,2005： SPE - 93951 - MS.

［6］ Duan Zhenhao,Sun Rui,Zhu Chen,et al. An Improved Model for the Calculation of CO2Solubility in Aqueous Solutions Containing Na+,K+,Ca2+, Mg2+,Cl-, and S［J］. Marine Chem,2006,98（2）：131 - 139.

［7］ Pruess K,Müller N. Formation Dry-Out from CO2Injection into Saline Aquifers：1. Effects of Solids Precipitation and Their Mitigation［J］. Water Resour Res,2009,45（3）：W03402.

［8］ Perrin J C,Benson S. An Experimental Study on the Inf uence of Sub-Core Scale Heterogeneities on CO2Distribution in Reservoir Rocks［J］. Trans Porous Media,2010,82（1）：93 -109.

［9］ Zhao Yuechao,Song Yongchen,Liu Yu,et al. Visualization and Measurement of CO2Flooding in Porous Media Using MRI［J］. Ind Eng Chem Res,2011,50（8）：4707 - 4715.

［10］ Zhang C,Oostrom M,Grate J W,et al. Liquid CO2Displacement of Water in a Dual-Permeability Pore Network Micromodel［J］. Environ Sci Technol,2011,45（17）： 7581 -7588.

［11］ Müller N. Supercritical CO2-Brine Relative Permeability Experiments in Reservoir Rocks：Literature Review and Recommendations［J］. Trans Porous Media,2011,87（2）：367 - 383.

［12］ Cinar Y,Riaz A,Tchelepi H. Experimental Study of CO2Injection into Saline Formations［J］. SPE J,2009,14（4）：588 - 594.

［13］ Xi Jiang. A Review of Physical Modelling and Numerical Simulation of Long-Term Geological Storage of CO2［J］. Appl Energy,2011,88：3557 - 3566.

［14］ Ellis J S,Bazylak A. Investigation of Contact Angle Heterogeneity on CO2Saturation in Brine-Filled Porous Media Using 3D Pore Network Models［J］. Energy Convers Manage, 2013, 68：253 - 259.

［15］ Chaczykowski M,Osiadacz A J. Dynamic Simulation of Pipelines Containing Dense Phase/Supercritical CO2-Rich Mixtures for Carbon Capture and Storage［J］. Int J Greenh Gas Con, 2012,9：446 - 456.

［16］ Audigane P,Gaus I,Czernichowski-Lauriol I,et al. Two-Dimensional Reactive Transport Modeling of CO2Injection in a Saline Aquifer at the Sleipner Site, North Sea［J］. Am J Sci, 2007,307（7）：974 - 1008.

［17］ Al Mansoori S K,Itsekiri E,Iglauer S,et al. Measurements of Non-Wetting Phase Trapping Applied to Carbon Dioxide Storage［J］. Int J Greenh Gas Con,2010,4（2）：283 - 288.

［18］ Yoo Seungyoul,Mito Yoshitada,Ueda Akira,et al. Geochemical Clogging in Fracture and Porous Rock for CO2Mineral Trapping［J］. Energy Procedia,2013,37：5612 -5619.

［19］ Lenormand R. Liquids in Porous Media［J］. J Physics: Condensed Matter,1990,2（S）：SA79.

［20］ Dullien F A L. Porous Media： Fluid Transport and Pore Structure［J］. AIChE J,1991,38（8）：1303-1304.

［21］ Ellis J S,Bazylak A. Dynamic Pore Network Model of Surface Heterogeneity in Brine-Filled Porous Media for Carbon Sequestration［J］. Phys Chem Chem Phys,2012,14（23）：8382 - 8390.

（編輯 李治泉）

·最新專利文摘·

用于C—O鍵氫解和加氫脫氧的過渡金屬催化劑

該專利公開了一種用于C—O鍵氫解和加氫脫氧的亞磷酰胺-金屬催化劑。該催化劑含有過渡金屬（如鎳、鈷和鐵）,過渡金屬與亞磷酰胺陰離子的配比為1∶1。與工業(yè)加氫脫氧用的普通催化劑相比,該金屬催化劑可以在較低溫度和壓力下,使含氧有機化合物的C—O鍵氫解。（Governors of the University of Alberta）/US 20140179954 A1, 2014-06-26

Simulation of Pore Network in Heterogeneous Brine-Filled Porous Media for Displacement of the Brine with Supercritical CO2

Li Yanming,Wu Rui,Cui Guomin,Cheng Shu
（Institute of New Energy Science and Engineering,School of Energy and Power Engineering, University of Shanghai for Science and Technology,Shanghai 200093,China）

A two-dimensional dynamic pore network model was established to simulate the displacement of brine in porous media with supercritical CO2for the CO2storage. The effect of heterogeneity on the displacement was investigated by comparing the supercritical CO2dynamic f ow process in three different two-dimensional network models. The simulation results indicated that, with the increase of the heterogeneity,the more unstable the invasion front of CO2,the lower the saturation of CO2at the breakthrough moment. The effective mechanism of the heterogeneity was revealed by investigating the CO2distribution in the stable displacement at different time. It was indicated that the decrease of the CO2saturation was mainly due to the increase of the trapping pore number in the network.

carbon dioxide storage；two-dimensional dynamic pore network model；supercritical carbon dioxide；heterogeneity

1000 - 8144（2014）08 - 0954 - 06

TQ 021.4

A

2014 - 03 - 15；［修改稿日期］ 2014 - 05 - 15。

李延銘（1988—）,男,山東省濰坊市人,碩士生,電話 18817849064,電郵 liyanming_n@163.com。聯(lián)系人：吳睿,電郵 ruiwu1986@gmail.com。

國家自然科學基金項目（51306124）；上海市青年科學基金項目（13ZR1458300）；低品位能源利用技術及系統(tǒng)教育部重點實驗室基金項目（LLEUTS - 201305）；上海大學生創(chuàng)新創(chuàng)業(yè)訓練計劃（SH2013007）；滬江基金研究基地專項（D14001）。