孔隙尺度下超临界CO2驱水两相流数值模拟

姜水生，赵万东，张莹，李培生，王昭太，钟源

（南昌大学机电工程学院，江西 南昌 330031）

孔隙尺度下超临界CO2驱水两相流数值模拟

姜水生，赵万东，张莹，李培生，王昭太，钟源

（南昌大学机电工程学院，江西 南昌 330031）

以咸水层封存二氧化碳（CO2）为研究背景，利用流体体积函数（Volume of Fluid，VOF）多相流模型建立孔隙尺度多孔介质计算模型，研究了超临界二氧化碳 （Sc-CO2）注入到含水多孔介质中的两相运移规律。对比分析了毛细管数、地质封存压力、Sc-CO2注射温度、两相表面张力系数、接触角等因素对两相运移速率以及驱替效率的影响，同时将不同毛细管数下的驱替效率与实验数据进行了对比。研究结果表明，随着毛细管数的增大，驱替效率先减小然后趋于稳定，数值模拟与实验数据吻合良好；在同一孔隙率下，在壁面表现为亲水性时，壁面润湿性越好，驱替速率越快，但驱替效率有所下降；同时毛细管数越小、地质封存压力越低、注射温度越高、张力系数越小驱替速率越快，且驱替效率越高。

超临界二氧化碳；咸水层封存；孔隙尺度；数值模拟；两相流；毛细管数；接触角

随着温室效应问题的日益严重，温室气体的排放与控制成为当今学者们研究的热点。碳元素的捕集与封存（carbon capture and storage，简称CCS）[1-2]是减少温室气体的重要措施之一，而CCS中地质封存被认为是中短期内减少大气中CO2含量的最有效方法之一。地质封存分为深部咸水层封存（简称盐水层封存）[3-4]、枯竭油气藏封存[5]以及不可开采的煤层封存，CO2封存示意图如图1所示。其中咸水层封存因为含地下水以及可渗透的多孔介质，且广泛分布于世界各地，故被认为是最具潜力的封存地点。超临界二氧化碳（Sc-CO2）扩散性好且动力黏度小，故将Sc-CO2注入到咸水层中成为封存CO2研究的热点。

图1 CO2地质封存示意图

近年来国内外学者对Sc-CO2咸水层封存运移机理进行了不同程度的研究，SUEKANE等[6]将Sc-CO2注入到使用玻璃珠堆积形成的人工多孔介质中，使用核磁共振（MRI）技术进行了可视化实验研究；WANG等[7]通过实验，使用连续和不连续的注射速率将Sc-CO2注射到二维均质咸水多孔介质中，研究了毛细管数对驱替过程中Sc-CO2饱和度的影响；TREVISAN等[8-9]通过实验构建均质和非均质多孔介质，进行岩心驱替实验，并使用X射线衰减分析流域饱和度，确定毛细阻力的影响因素；SASAKI等[10]通过直接数值模拟方法对不同压力下的Sc-CO2注入到地下岩石进行了研究；CHOI等[11]为了研究CO2封存过程中的安全性，建立了孔隙尺度下的多孔介质模型，使用CFD软件STAR-CCM研究了不同注射温度下的CO2驱替过程中的流动和传热过程。姜培学课题组[12]采用核磁共振设备，通过实验研究了Sc-CO2和盐水流过玻璃珠过程中浮升力和孔隙度的影响，同时使用可视化技术研究了岩心中超临界压力下CO2-水相对渗透率与饱和度的关系；XU等[13]在孔隙尺度条件下直接求解N-S方程方法对比研究了Sc-CO2驱水和空气驱水过程，但只模拟了一种方案下的驱替过程；罗庶等[14]使用数值模拟方法研究了地质封存Sc-CO2突破压力梯度的问题，但对多种因素影响下的驱替过程中两相界面变化、突破过程、Sc-CO2饱和度等研究较少；蒋兰兰[15]通过MRI技术以及CT成像技术，对玻璃砂多孔介质进行了单相和多相驱替实验，并使用瞬态和稳态法研究了相对渗透率饱和曲线。

国内外学者对Sc-CO2驱水过程进行了大量的实验研究，然而在实验过程中，MRI在压力的选定上存在限制，且实验周期长，成本大；在驱替实验中，由于其构造的多孔介质模型单一，表面张力以及接触角等因素较难改变；同时由于出口末端效应和入口段效应的存在[15]，不能精确地测量多孔介质两端的压力，因此存在一定的实验误差。而数值模拟方法在选取不同的孔隙率、毛细管数、注射压力、表面张力、接触角等参数较容易。故本文使用流体体积函数（VOF）方法对Sc-CO2驱水进行了孔隙尺度单相驱替研究，研究不同注射压力和温度等因素对两相界面变化、突破过程、驱替速率、驱替效率以及残余水饱和度的影响，将有助于CO2地质封存场地的选取以及Sc-CO2注射条件的选取。

1 物理模型、控制方程和数值求解方法

1.1 物理模型

描述真实多孔介质三维模型十分困难，而通过孔隙尺度下的有限微小单元可以描述不同孔隙率的多孔介质；假设固体颗粒大小均匀，同时颗粒排布均匀，如实验中人工填充的玻璃珠模型。改变两个简单的不同直径大小的圆球来实现流体域的建立，例如简单立方体（simple cube）、体积中心立方体（body-centered cube）以及面中心立方体（face-centered cube）方式。XU和JIANG[16]已经通过数值模拟研究了不同排列方式下，考虑壁面的剪切速度以及从实验方面对比3种模型的误差，其中体积中心立方体模型的误差小于5%。因此本文使用体积中心立方体模型来描述孔隙尺度模型，将孔隙结构中固体域分离即可得出流体计算域，计算域孔隙尺度单元以及尺寸见图2与表1。求解域离散成CFD计算网格，并将流体域分为质量流量入口、压力出口、壁面、对称边界。网格边界层局部放大图以及边界条件见图3，流体计算域离散网格使用ANSYS生成。

图2 孔隙尺度结构单元以及体积中心立方体模型

表1 单元尺寸以及孔隙度

图3 离散网格边界条件及横截面局部网格放大图

1.2 控制方程

Sc-CO2驱替含水多孔介质是一个复杂的多相流过程，涉及到相与相之间的界面作用力，如表面张力的作用，同时还有Sc-CO2、水以及壁面之间的三相接触角的问题。VOF方法已经广泛应用于多相流仿真，如模拟溃坝、沸腾与冷凝、气泡的运动、自由界面等模拟[17-18]。VOF基于追踪技术应用于固定的欧拉网格上的自由界面，求解模型中引入一个体积比函数α，假设αq为其中主相的体积分数，也即该相占整个网格体积的体积分数[19]。通过此方法能较好地实现对自由面以及两相流内部运动交界面的追踪，因此可以实现Sc-CO2对驱替水过程中的内部运移界面的追踪。

多相界面的追踪通过求解连续性方程主相来完成，对于其中的某一相连续性方程表达如式(1)[19]。

式中，mpq为从p到q的质量转移，相反mpq是p到q的质量转移，同时首相的体积分数满足如下式(2)定值。

其中体积分数方程通过显示时间离散求解，在显示求解过程中，通过有限差分插值格式应用于体积分数计算的前一个时间步，计算格式如式(3)。

式(3)中，n+1为当前的计算时间步；n代表上一个时间计算步。

通过一个动量方程求解整个计算域，得到的速度、温度场应用于整个流场，式(4)是动量方程的表达式。

在控制体的输运方程中密度和黏度由混合相确定，在两相求解模型中，密度和黏度如式(5)、式(6)定义。

求解过程中添加能量项、能量方程的计算如式(7)、式(8)。

式(7)中，ρ为式(5)计算得到的平均密度；keff为有效热导率。

求解过程引入表面张力，也即是CSF模型，该模型把表面张力看作求解过程中的体积力添加在动量方程中，引入的体积力如式(9)。

式(9)中，σqp为两相表面张力系数；ρ为式(5)计算得到的平均密度。

三相界面接触角采用固定接触角，接触角定义如式(10)。

式(10)中，θw为壁面接触角；为分别为壁面法向和切向单位向量。

1.3 数值求解方法与边界条件

由于多孔介质结构的复杂性，为了表征边界层网格，因此离散网格采用四面体网格。使用Fluent16.0中的VOF方法模拟Sc-CO2驱替含水多孔介质过程，求解模型设定为瞬态VOF模型，考虑体积力打开implicit body force选项；由于驱替过程雷诺数较小，因此求解模型选择层流模型；模拟过程中有能量的传递，故加入能量方程；考虑到瞬态求解的收敛性，科拉数设定为0.5，压力与速度的耦合采用PISO算法，动量方程采用二阶迎风格式，同时采用几何重构的方法获取Sc-CO2驱替的界面形状；残差收敛设定为1×10–4，求解过程中设定时间步长为1×10–6s。入口边界条件设定为Sc-CO2质量流量入口，质量流量为5×10–4g/s，出口设定为压力出口，静压为0，壁面选择无滑移边界条件，环境变量根据选取的地质封存深度来决定，含水多孔介质初始温度设定为333K，其他边界设定为对称边界条件。

VOF材料的添加选取水和Sc-CO2。对于Sc-CO2材料的添加，可通过界面操作命令来实现，界面操作命名中输入NIST Real Gas Models，选取单组分Sc-CO2材料，并在界面命名中创建计算节点表。本次数值模拟方案以及Sc-CO2、水的物性参数分别见表2和图4、图5。

1.4 数据处理方法

对于评定Sc-CO2驱替过程参数有渗透率、绝对渗透率、相对渗透率等。对于服从达西定律的多孔介质流动，渗透率可以计算出，而绝对渗透率是在多孔介质内的单相流动评定参数。对多孔介质内的绝对渗透率定义如式(11)[15]。

表2 数值模拟方案及相应参数的选取

图4 Sc-CO2物性参数随温度和压力变化关系

式中，Q为流量；L为长度；μ为流体动力黏度；ΔP为压降。多孔介质渗透率只与多孔介质本身有关，而与流体属性无关，本文采取稳定流动获取多孔介质渗透率，计算出的渗透率与实验数据吻合[15]。

计算过程中采用量纲为1毛细管数来量化不同方案，毛细管数定义如式(12)[15]。

式中，μ为Sc-CO2动力黏度，Pa·s；V为速度，m/s；σ为两相表面张力系数，N/m。

计算过程中水和Sc-CO2的饱和度定义如式(13)、式(14)。

2 结果与讨论

2.1 验证性计算

2.1.1 网格无关性验证

数值模拟过程中，网格的数量决定了计算的精度，因此在进行计算前，本文先验证了网格数量对计算精度的影响。分别选取了网格数量225671、446242、682542进行计算精度验证，表3给出了3种网格数量下，计算稳定时水的饱和度。计算结果表明，网格数量较小时，计算误差较大，但网格数量提高大约3倍时，以网格数量最多的计算结果作为参考，中等网格数量相对误差为4.81%。因此从计算精度和计算资源考虑，本次选中等网格数量，即446242。

表3 网格无关性验证

2.1.2 VOF模型可行性验证

为了验证VOF模型计算Sc-CO2的可行性。分别使用单相层流模型和VOF层流模型（Sc-CO2体积分数为1）进行验证，入口使用质量流量入口，入口流量大小都为0.0005g/s，出口使用压力出口[14]。两种模型对应的进出口压降相差为3.62×10–5Pa，由此验证VOF模型可用于计算Sc-CO2驱替含水多孔介质。

2.2 Sc-CO2驱替

2.2.1 Sc-CO2驱替过程两相界面分析

整个驱替过程中，Sc-CO2以恒定质量流量从左侧入口进入，初始温度设定为333K，环境压力为15MPa，表面张力系数设定为0.045N/m，图6给出了两相界面与壁面接触角分别为30°和60°时整个驱替过程中不同时刻的两相分布图。计算结果表明，在t=0.0005s时刻，Sc-CO2运移到第一个对称颗粒，两相界面对称分布；在t=0.001s，Sc-CO2运移到第一个上下非对称的小颗粒处，30°与60°接触角出现了不同结果，30°接触角选择小孔隙吼道进行突破，而60°接触角选择大孔隙吼道进行突破。这是因为当吼道较小时，60°接触角中吼道较小的地方形成毛细压力，故对于非润湿相驱替润湿相，毛细压力变成阻力，而相比于较小的接触角时，毛细压力成为驱替动力，因此驱替过程中出现不同现象。随着Sc-CO2连续的注入到多孔介质中，在壁面形成的润湿相液膜也逐渐被驱替，在t=0.0027s时刻通过最后一个小颗粒。

2.2.2 毛细管数的影响

图6 不同时刻驱水两相分布图

由于多孔介质中含有不均匀的孔隙，因此为了探讨孔隙大小对驱替过程的影响，本文通过改变Sc-CO2的注入流量来获取不同的毛细管数。图7给出了数值模拟结果与WANG等[7]在2012年通过实验获取毛细管数影响Sc-CO2饱和度数据图。计算结果表明，采用数值模拟下的不同毛细管数的Sc-CO2饱和度能很好的与实验数据吻合；随着注入流量的增加，Sc-CO2饱和度下降，但随着流量的再次加大，毛细管数达到Sc-CO2的饱和度变化不再明显，最后趋于平稳。

图7 Sc-CO2饱和度与毛细管数关系

2.2.3 Sc-CO2注射压力的影响

咸水层地质封存深度大约为1000～2500m，不同深度下的封存过程中的环境压力会有所不同。因此本文通过改变模拟中的环境压力来实现不同封存压力，压力方案A、B、C分别选取10MPa、15MPa、20MPa。图8给出了在t=0.0012s时刻Sc-CO2驱水过程两相分布图；3种方案残余水饱和度随驱替时间发生的变化关系如图9所示。数值模拟结果表明，水的饱和度在10%以下时，3种方案所花费的时间大致为0.0015s、0.0026s、0.0033s，因此封存压力越低，驱替速率越快，其次由最后驱替达到稳定状态结果表明，环境压力越低，残余水饱和度越低。这是因为地质封存选取深度越高，即驱替环境压力越高，这时Sc-CO2的密度越大，在同样的质量流量下，体积流量反而越小，因此Sc-CO2样能运移到更远处；同时驱替所遇到的阻力越大，导致驱替速率下降，但是驱替过程所花费的时间并不与压力变化成正比，压力越高驱替速率变化减小。

图8 t=0.0012s时刻两相分布图

图9 H2O饱和度与驱替时间关系

2.2.4 Sc-CO2注射温度的影响

由于地质封存选取地层深度约为1000～2500m，该地层的环境温度大致为313～383K，本文研究选取含水多孔介质温度为313K，采用不同注射温度下的Sc-CO2来驱替，方案D、E、F温度分别为313K、323K、333K。图10给出了Sc-CO2在3种温度条件、t=0.002s时刻的两相分布图；图11给出了3种温度条件下残余水饱和度随驱替过程的变化图。计算结果表明，Sc-CO2所选取的温度越高，驱替速率越快，在驱替中期，驱替速率大致和Sc-CO2温度成线性关系，这是因为随着温度增加，Sc-CO2黏度越小，流动性强，即在相同时刻，方案F中Sc-CO2运移到更远处，但是随着时间的推移，驱替稳定后，驱替效率大致相等。

2.2.5 表面张力的影响

图10 t=0.002s时刻两相分布图

图11 H2O饱和度与驱替时间关系

图12 t=0.0018s时刻两相分布图

图13 H2O饱和度与驱替时间关系

为了研究表面张力对两相驱替过程的影响，方案G、I分别选取表面张力为0.03N/m、0.06N/m来进行研究。图12给出了两种方案在t=0.0018s时刻两相分布图，图13给出了残余水饱和度随驱替时间发生变化关系图。计算结果表明，随着表面张力的增加，驱替速率下降，同时驱替效率也下降。这是因为随着表面张力的减小，Sc-CO2在驱替过程中两相毛细压力越小，因此在相同的时间内，图12方案G中Sc-CO2能运移到更远处；同时由于表面张力的增加，润湿相更易黏附在固体壁面，导致驱替效率有所下降，故残余水饱和度更大。

2.2.6 接触角的影响

不同的地质条件下，岩石的浸润性会有所不同，也即表现为三相（岩石、水、Sc-CO2）接触角不一样，因此本文以此为研究出发点，选取不同的接触角（方案J、K、M的接触角分别为30°、60°、90°），在其他条件不变的情况下进行数值模拟计算。3种不同接触角方案在t=0.0016s时两相分布见图14；图15给出了驱替过程中残余水饱和度随时间变化的过程。数值模拟结果表明，接触角越小，驱替速率越快，但驱替效率有所下降。这是因为当岩石壁面表现为亲水性时，随着接触角的减小，润湿相润湿性越好，这样会在壁面形成液膜，使得非润湿相Sc-CO2在驱替过程中阻力减小，也即表现为方案J中Sc-CO2运移到更远处；但是接触角越小，润湿相在壁面形成液膜越厚，Sc-CO2驱替效率下降。

图14 t=0.0016s时刻两相分布图

图15 H2O饱和度与驱替时间关系

3 结论

本文建立了孔隙尺度下各向同性的三维多孔介质模型，使用数值模拟方法研究了Sc-CO2注入含水多孔介质的两相驱替过程。分析了毛细管数、地质封存压力、Sc-CO2注入温度、两相表面张力、接触角影响因素对Sc-CO2驱替含水多孔介质中的运移规律，得到如下结论。

毛细管数越大，Sc-CO2饱和度先下降，但随着毛细管数对数进一步加大到–5.5时，Sc-CO2的饱和度趋于平稳；地质封存地点越深，达到同一残余水饱和度驱替所花费的时间越长，但花费的时间并不与封存压力成线性关系，当压力大于15MPa时，驱替速率大大下降；Sc-CO2注射温度越高，驱替速率越快；随着两相表面张力系数的增加，驱替速率减小，同时驱替效率也越低；在壁面表现为润湿性条件下，接触角越小，驱替速率越快，同时较小的接触角，将导致过大的毛细压力，这可能成为驱替的动力。

[1] BODE S，JUNG M. Carbon dioxide capture and storage—liability for non-permanence under the UNFCCC[J]. International Environmental Agreements: Politics，Law and Economics，2006，6（2）：173-186.

[2] BENTHAM M，KIRBY M. CO2storage in saline aquifers[J]. Oil &Gas Science and Technology，2005，60（3）：559-567.

[3] CHADWICK R，ARTS R J，BERNSTONE C，et al. Best practice for the storage of CO2in saline aquifers-observations and guidelines form the SACS and CO2STORE projects[M]. British Geological Survey，2008.

[4] THOMAS S. Enhanced oil recovery- an overview[J]. Oil & Gas Science and Technology-Revue de l'IFP，2008，63（1）：9-19.

[5] BACHU S，SHAW J C，PEARSON R M. Estimation of oil recovery and CO2storage capacity in CO2EOR incorporating the effect of underlying aquifers[C]//Society of Petroleum Engineers，2004.

[6] SUEKANE T，SOUKAWA S，IWATANI S，et al. Behavior of supercritical CO2injected into porous media containing water[J].Energy，2005，30：2370–2382.

[7] WANG Y，ZHANG C Y，WEI N，et al. Experimental study of crossover from capillary to viscous fingering for supercritical CO2-water displacement in a homogeneous pore network[J].Environmental Science & Technology，2013，47：212-218.

[8] TREVISAN L，PINI R，CIHAN A，et al. Experimental investigation of supercritical CO2trapping mechanisms at the intermediate laboratory scale in well-defined heterogeneous porous media [J].Energy Procedia，2014，63：5646–5653.

[9] TREVISAN L，CIHAN A，FAGERLUND F，et al. Investigation of mechanisms of supercritical CO2trapping in deep saline reservoirs using surrogate fluids at ambient laboratory conditions[J].International Journal of Greenhouse Gas Control，2014，29：35–49.

[10] SASAKI K，FUJI T，NIIBORI Y，et al. Numerical simulation of supercritical CO2injection into subsurface rock masses[J]. Energy Conversion and Management，2008，49：54–61.

[11] CHOI H S，PARK H C，HUH C，et al. Numerical simulation of fluid flow and heat transfer of supercritical CO2in micro-porous media[J].Energy Procedia，2011，4：3786–3793.

[12] 马瑾，胥蕊娜，罗庶，等. 超临界压力CO2在深部咸水层中运移规律研究[J]. 工程热物理学报，2012，33（11）：1971-1975.MA J，XU R，LUO S，et al. Core-scale experimental study on supercritical-pressure CO2migration mechanism during CO2geological storage in deep saline aquifers[J]. Journal of Engineering Thermophysics，2012，33（11）：1971-1975.

[13] XU R，LUO S，JIANG P X. Pore scale numerical simulation of supercritical CO2injection into porous media containing water[J].Energy Procedia，2011，4：4418-4424.

[14] 罗庶，胥蕊娜，姜培学. CO2地质封存突破压力梯度的孔隙尺度数值模拟[J]. 工程热物理学报，2011，32（5）：819-823.LUO S，XU R N，JIANG P X. Pore-sclae numerical simulation on breakthrough pressure gradient during CO2geological storage[J].Journal of Engineering Thermophysics，2011，32（5）：819-823.

[15] 蒋兰兰. CO2地质封存多孔介质内气液两相渗流特性研究[D]. 大连：大连理工大学，2014.JIANG L L. Study on gas/liquid flow characteristic in porous media during CO2geological storage[D]. Dalian: Dalian University of Technology，2014.

[16] XU R N，JIANG P X. Numerical simulation of fluid flow in microporous media[J]. International Journal of Heat and Fluid Flow 2008，29（5）：1447-1455.

[17] POURYOUSSEFI S M，ZHANG Y W. Numerical investigation of chaotic flow in a 2D closed-loop pulsating heat pipe[J]. Applied Thermal Engineering，2016，98：617–627.

[18] CARCIOF B A M，PRAT M，LAURINDO J B. Homogeneous volume-of-fluid（VOF）model for simulating the imbibition in porous media saturated by gas[J]. Energy & Fuels，2011，25（5）：2267-2273.

[19] HIRT C W，NICHOLS B D. Volume of fluid （VOF） method for the dynamics of free boundaries[J]. Journal of Computational Physics，1981，39：201-225.

Pore-scale two-phase numerical simulation of supercritical carbon dioxide displacement of water

JIANG Shuisheng，ZHAO Wandong，ZHANG Ying，LI Peisheng，WANG Zhaotai，ZHONG Yuan
（School of Mechanical and Electrical Engineering，Nanchang University，Nanchang 330031，Jiangxi，China）

The computational method of pore-scale porous media is established by using Volume of Fluid（VOF）numerical simulation method based on the background of CO2storage in Saline aquifer，and the migration mechanism of supercritical carbon dioxide（Sc-CO2）injection into porous media containing water is studied. The effects of capillary number，geological storage pressure，Sc-CO2injection temperature，two-phase surface tension coefficient and contact angle on the two-phase migration rate and the displacement efficiency were analyzed. At the same time，displacement efficiency was compared with experimental data under different capillary numbers. The result showed that with the increase of capillary number，displacement efficiency decreases first and then becomes stable，and the numerical simulation of displacement efficiency agrees well with experimental data under different capillary numbers. Under the same porosity，when wall surface is hydrophilic，wall surface wettability is better. The faster the rate of displacement，while the displacement efficiency decreased. At the same time，the lower geological storage pressure，the higher the injection temperature.Small surface tension coefficient has a higher displacement rate and efficiency.

supercritical carbon dioxide；saline aquifer storage；pore-scale；numerical simulation；two-phase flow；capillary number；contact angle

TK124；TE122

A

1000–6613（2017）11–3955–08

10.16085/j.issn.1000-6613.2017-0567

2017-03-31；修改稿日期2017-06-22。

国家自然科学基金（51566012，11562011）、江西省科技厅支撑项目（2009BGA01800）及江西省研究生创新基金（YC2017－S056）项目。

姜水生（1957—），男，硕士，教授。联系人张莹，博士，教授，研究内容为复杂热流场模拟。E-mail：yzhan2033@163.com。