非均质多孔介质内超临界二氧化碳驱替盐水过程的孔隙网络模拟

李延铭，吴 睿，崔国民，程 树

（上海理工大学 能源与动力工程学院新能源所,上海 200093）

环境与化工

非均质多孔介质内超临界二氧化碳驱替盐水过程的孔隙网络模拟

李延铭，吴 睿，崔国民，程 树

（上海理工大学 能源与动力工程学院新能源所,上海 200093）

针对CO2在非均匀多孔介质内的储存过程,建立二维动态孔隙网络模型,对含盐水的介质中超临界CO2驱替盐水的过程进行模拟计算。提出用非均匀性来衡量二维动态孔隙网络模型中喉道的分布特性,通过建立3种非均匀性不同的二维动态孔隙网络模型,分析非均匀性对CO2封存的影响。模拟计算结果表明,非均匀性越大的二维动态孔隙网络模型,驱替过程中超临界CO2的前进端越不稳定,且最后突破时的超临界CO2饱和度越低；通过分析稳定驱替下不同时刻的CO2分布,揭示非均匀性的影响机理；非均匀性越大,陷住的含盐水孔隙更多是导致突破时超临界CO2饱和度越低的主要原因。

二氧化碳封存；二维动态孔隙网络模型；超临界二氧化碳；非均匀性

CO2捕集与储存技术（CCS）是当今减缓大气中CO2浓度过高和温室效应的有效措施［1］。CCS是指将大型CO2排放源排放出的CO2收集并进行长时间储存,使其与大气相隔绝的一种技术。地下深层多孔的地质结构,例如废弃的油气田和沉积含盐水地层,可以为CO2提供足够的储存空间。沉积含盐水地层是全球分布范围较广的一种理想储存空间,其中的部分砂层、砂岩、石灰岩具有较高的孔隙度和渗透率,可以为CO2储存提供储存空间,而页岩和泥岩虽然有较高的孔隙度,但渗透率很低,这种低渗透率的岩层称之为岩盖［2］,可以很好地对CO2实现隔离。

实际操作中,CO2封存是一个非常复杂的过程,其封存过程主要受到黏性力、浮力和毛细压力的影响［3-5］。封存在地下的CO2部分溶解可能会导致地下水酸化形成离子沉淀［6］,造成地质结构发生改变［7］。因此地质结构的孔隙度、渗透率、稳定性和地下水盐度等都是CO2封存所必须考虑的因素。

高效安全地实现CO2封存是目前研究的重点,这方面的研究主要分为可视化实验［8-12］以及数值计算［13-16］两大部分。目前,在数值计算方面主要采用建立孔隙网络模型研究CO2的封存特性。孔隙网络模型从微观上揭示了多孔介质中两相流动的原理,为了解岩石微观结构对CO2封存的影响提供了理论依据。

图1 二维动态孔隙网络模型Fig.1 Two-dimensional dynamic pore network model.

本工作建立二维动态孔隙网络模型,模拟地质条件下储存CO2的过程,分析不同的孔隙结构对CO2驱替盐水过程的影响。

1 模型描述

将无规则的岩石微观结构抽象为理想的几何形状,建立二维动态孔隙网络模型（见图1）。二维动态孔隙网络模型用来描述一个含润湿相（原流体,如盐水）和非润湿相（侵入流体,如超临界CO2）的正方形晶格。由图1可见,二维动态孔隙网络模型由圆柱形喉道和圆形孔隙构成,每个孔隙有4条喉道连接形成格子型网络。其中,喉道代表实际介质中较小的空间,孔隙代表较大的空间。二维动态孔隙网络模型中孔隙半径（rpore）和喉道长度（Lthroat）均为定值,喉道半径大小介于rmin和rmax之间,非均匀性用式（1）表示。

2 计算过程

二维动态孔隙网络模型所有孔隙中一开始充满原流体,侵入流体自下向上侵入。边界条件为：二维动态孔隙网络模型左右两侧为封闭的边界；下端为体积流量恒定,其大小给定；上端为压力恒定,其大小为0。

为了简化计算过程,对二维动态孔隙网络模型及流体流动过程做以下假设：1） 超临界CO2与固体表面间的接触角为定值（165°）；2） 虽然喉道决定导流率的大小,仍认为喉道的体积远小于孔隙的体积；3） 对于每一个喉道,只有一种流体可以滞留；4） 喉道中流体的流动状态假设为层流并遵从哈根-泊肃叶定律；5） 忽略孔隙的流动阻力；6）两种流体均为不可压缩流体。

整个计算过程可分为4步：1）压力场及流场的预处理；2）压力场及流场的求解；3）时间步长的确定；4）弯液面的更新。驱替过程从下边界有侵入流体流入开始,到上边界有侵入流体流出（突破）为止。

两相流体间的接触面称为弯液面（见图1）,其稳定性由式（2）确定。

若孔隙中充满侵入流体且Δp≤0,则弯液面是不可流动的,即处于稳定状态；若Δp＞0,则弯液面是可流动的,即处于不稳定状态；若孔隙只有一部分被侵入流体占据,由于孔隙中流动阻力大小不计,所以此时弯液面仍处于不稳定状态。

首先,根据式（2a）判断弯液面的稳定性。若有不稳定的弯液面,则不需要预处理,若所有弯液面均为稳定,则进行预处理。当所有弯液面均处于稳定状态时,为了保证驱替过程的连续,假设最大孔道内的弯液面处于不稳定状态,且侵入流体占据此孔道。侵入流体在喉道内的流量由式（3）计算。

式（3b）的意义为当相连接两孔隙间的压力差小于门阙压力时（Δp≤0）,喉道内侵入流体的流量为0。原流体在喉道内的流量由式（4）计算。

根据质量守恒定律,并联立式（2）～（4）可获得压力场和流场。压力场及流场确定后,需要一个时间步长推动弯液面前进。先求得每个与弯液面相连接的孔隙填满所需的时间,然后取其中的最小值确定为时间步长,见式（5）。

时间步长确定后,由式（6）迭代更新弯液面。

陷住的孔隙示意图见图2。

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

由图2可见,若出现侵入流体将原流体包围的情况,则被困住的原流体不能再移动［17-18］,即孔隙被陷住,陷住的孔隙可能是一个也可能是多个。

3 结果与讨论

采用两个无量纲参数：毛细数（Ca）和黏度比（M）,表征流体的流动状态。Ca是黏性力与毛细力的比值,见式（7）。M为两液相黏度的比值,见式（8）。

通过改变超临界CO2的流速和黏度来改变Ca和M。

二维动态孔隙网络模型的参数见表1。

表1 二维动态孔隙网络模型的参数Table 1 Parameters of the two-dimensional dynamic pore network model

3种二维动态孔隙网络模型结构的非均匀性分别为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驱替盐水的效率。

不同参数下突破时超临界CO2在二维动态孔隙网络模型中的分布见图3。

图3 不同参数下突破时超临界CO2在二维动态孔隙网络模型中的分布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种二维动态孔隙网络模型结构的驱替过程均呈现出沿前进方向不稳定扩散的黏性指进；M较大、Ca较大时（M=1.0；Ca=10-1,10-3）, 3种二维动态孔隙网络模型结构的驱替过程均较为稳定；Ca较小时（Ca=10-5；M=0.1,1.0）,3种二维动态孔隙网络模型结构的驱替过程均呈现出沿任意方向不稳定扩散的毛细指进。

黏性指进主要是由于黏性力占优时两种流体黏度差异造成的一种不稳定驱替现象,而毛细指进是由于毛细力占优时喉道随机分布造成的一种不稳定驱替现象［19-20］。在黏性力占优的情况下,两种流体黏度越接近,侵入流体的前进端越稳定。随非均匀性的增大,陷住的含盐水孔隙越来越多甚至出现局部大面积陷住的情况,而且超临界CO2的前进端越不稳定,尤其在稳定驱替情况下,孔隙由单个到多个、从少量到大量的陷住更为明显。

不同参数下非均匀性对超临界CO2饱和度的影响见图4。由图4可知,超临界CO2饱和度随非均匀性的增大而减小。

图4 不同参数下非均匀性对超临界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饱和度差异,将稳定驱替情况下χ=1.8和χ=0.4时二维动态孔隙网络模型抽象为两种简单的局部示意图,见图5。

图5 不同二维动态孔隙网络模型结构的侵入过程Fig.5 Displacement process in different two-dimensional network models.

由式（2b）可知,超临界CO2从网络模型下端选择最大喉道开始侵入（见图5（a））；随驱替过程的进行,为保证流量的恒定,底部压力会有一定的波动,此时对于χ=0.4的结构若底部压力超过右侧较小喉道的门阙压力时则较小喉道开始有超临界CO2侵入,对于χ=1.8的结构较小喉道的门阙压力过大时仍难以开始侵入（见图5（b））；随驱替过程的进一步进行,χ=0.4的结构中超临界CO2的体积与χ=1.8的结构中超临界CO2的体积差距增大（见图5（c））；χ=1.8的结构右下角出现陷住,最终导致在不同的非均匀性下超临界CO2饱和度存在差异（见图5（d））。

稳定驱替下超临界CO2前进端的形态变化见图6。由图6（a）可知,驱替过程开始阶段,不同的二维动态孔隙网络模型结构中超临界CO2的前进端就开始呈现了不同的现象。非均匀性越大的二维动态孔隙网络模型结构中超临界CO2的前进端越不稳定,随驱替过程的进行,这种前进端稳定性的差异越来越明显,见图6（b）～（d）。

结合图5和图6可知,非均匀性越大的二维动态孔隙网络模型结构中陷住的含盐水孔隙越多,这是导致突破时超临界CO2饱和度越低的主要原因。

图6 稳定驱替下超临界CO2前进端的形态变化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］,但通过模拟计算结果发现,地质微观结构的非均匀性也对CO2封存有重要影响。非均匀性较大的结构孔隙率更高,但封存过程中的超临界CO2饱和度较低。相反非均匀性较小的结构有较低的孔隙率,但封存过程中的超临界CO2饱和度较高。因此在实际CO2封存过程中应该考虑地质微观结构非均匀性的影响。

4 结论

1） 为了解决地质微观结构内喉道分布非均匀特性对超临界CO2封存的影响,建立了3种非均匀性不同的二维动态孔隙网络模型。通过对二维动态孔隙网络模型喉道分布非均匀特性的数值模拟研究,发现不同的喉道分布模型对CO2的驱替效率有重要影响。

2）非均匀性越大的孔隙网络模型中超临界CO2的前进端越不稳定,导致最后突破时的超临界CO2饱和度越低。基于稳定驱替下的分步计算结果,发现非均匀性较大的二维动态孔隙网络模型中门阙压力分布不均匀导致了超临界CO2前进端的不稳定,而更多陷住的含盐水孔隙是导致突破时超临界CO2饱和度低的主要原因。

符 号 说 明

Ca 毛细数

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］ 迭代步数

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。与工业加氢脱氧用的普通催化剂相比,该金属催化剂可以在较低温度和压力下,使含氧有机化合物的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。联系人：吴睿,电邮 ruiwu1986@gmail.com。

国家自然科学基金项目（51306124）；上海市青年科学基金项目（13ZR1458300）；低品位能源利用技术及系统教育部重点实验室基金项目（LLEUTS - 201305）；上海大学生创新创业训练计划（SH2013007）；沪江基金研究基地专项（D14001）。