天然气水合物降压开采数值模拟研究

王文博， 刘晓， 崔伟， 肖加奇

齐鲁工业大学(山东省科学院), 济南 250353

0 引言

天然气水合物是在低温、高压条件下形成的一种冰晶状化合物(刘军等, 2015; 刘建军等, 2017),具有分布广、存储量大等优点,被认为是21世纪最具有开发潜力的能源之一.

天然气水合物通常以固态的形式填充于多孔介质中,不具备流动性,可行的开采方法是使水合物分解产生天然气和水,然后将天然气采集(Konno et al.，2010; 李刚等, 2011; Chong et al., 2016; 吴传芝等, 2016).由于水合物成为了岩石结构的一部分,其分解会使储层胶结强度、孔隙度、渗透率和地质结构等发生改变(Delli and Grozic, 2013; 李栋梁等, 2019),对天然气水合物开采研究必须关注水合物分解过程及其对储层参数的影响.

目前,天然气水合物开采的方法主要有降压法、注热法、化学试剂法以及CO2置换法等.其中降压法被认为是最经济有效的天然气水合物开采方式,该方法将压力降低至水合物平衡压力之下,促使水合物分解(梁海峰和宋永臣, 2008; 朱自浩, 2015; 刘昌岭等, 2019; 陈强等, 2020).在加拿大Mallik冻土区(栾锡武等, 2007)、日本南海海槽(Fujii et al., 2013; 赵克斌, 2019)以及中国神狐海域(万义钊等, 2018; Li et al., 2018;叶建良等, 2020)等水合物试采均采用了降压法并取得一定的成功,证明了降压法开采水合物的可行性.在降压开采过程中,水合物分解区和尚未分解区间存在分解前沿,随着开采时间的推进,分解前沿发生移动,因此,对分解前沿形状及其移动规律的研究成为了降压开采研究的核心之一.

近年来,国内外许多学者对天然气水合物的分解过程进行了研究.Verigin等(1980)、Ji等(2003)将水合物的分解看成是移动边界消融的过程.Selim和Sloan(1985)、Ullerich等(1987)认为水合物的分解只发生在分解界面处.Tsypkin(1999)、Yang等(2012)、Vasil'ev等(2006)、 杜庆军等(2007)、周锡堂等(2007)、 张旭辉等( 2012)通过降压、注热分解实验或数值模拟认为水合物分解存在一个较大的分解区.Selim和Sloan(1985)、 Jamaluddin等(1989)、李淑霞等(2005)研究认为水合物分解前沿与开采时间呈线性关系.Yousif等(1990)、 Yousif等(2013)、张旭辉等(2012)研究表明水合物分解前沿随着开采时间的增加移动速度逐渐变慢.喻西崇等(2004)、刘乐乐等(2013)、赵振伟(2013)通过数值模拟得出水合物分解前沿与开采时间的平方根成正比.Sung等(2000)研究表明在一维模型中,分解前沿与时间呈线性关系,而在三维模型中,分解前沿移动速度逐渐变慢.Sun等(2005)在建立的数学模型中得出小尺度条件下分解前沿立即到达介质末端,而在大尺度条件下分解前沿逐渐推进且推进速度变慢.李淑霞等(2019)利用CMG-STARS模拟软件得出水合物的降压分解存在一个较大的分解区,分解过渡带外沿与开采时间近似线性关系,分解过渡带内沿与开采时间近似指数关系.

对水合物分解进行研究时,许多学者在直角坐标系下模拟了一个较小的、封闭的水合物藏区域,因此得出了分解前沿移动随开采时间而越来越快的假象.在实际的降压开采过程中,水合物藏是一个较大的区域,在压力驱动力的作用下,水合物分解区域沿径向扩散,因此,在柱坐标系下建立数学模型更为合适,且模型外边界必须足够大.基于降压开采机理,本研究建立了柱坐标系下的降压开采数学模型,应用有限差分法,对模型进行了求解并以此为工具,分析了水合物降压开采过程的压力、水合物饱和度、渗透率的变化规律,以及水合物分解过渡带的移动规律,研究了不同降压开采参数对水合物分解过程的控制作用.

1 物理模型及其基本假设条件

研究中天然气水合物藏物理模型和假设条件如下:

(1) 在均质的水合物储层中心有一半径r0的直井,储层半径re(图1).整个水合物层采用等比级数网格划分方法(李淑霞和谷建伟,2008),网格序号和径向距离的关系见图2.

图1 水合物降压分解模型示意图Fig.1 Schematic of depressurization decomposition of natural gas hydrate

图2 径向网格划分示意图Fig.2 Radical grid distribution

(2)完全分解区水合物全部分解,水合物饱和度为0;未分解区水合物未发生分解,水合物饱和度为初始值;完全分解区到未分解区之间的区域为分解过渡带,内沿半径r1,该处水合物饱和度刚开始大于0;分解过渡带外沿半径r2,该处水合物饱和度刚好低于初始值.

(3)将水合物分解过程视为是一个移动边界问题,r=r(t)表示移动边界的位置:完全分解区(r0≤r≤r1)、分解过渡带(r1≤r≤r2)和未分解区(r2≤r≤re).

(4)完全分解区只含有气、水两相,储层的渗透率为储层绝对渗透率;分解过渡带含有水合物、气和水三相,水合物不参与流动, 储层的渗透率为小于绝对渗透率而大于初始渗透率,水合物饱和度、渗透率随着时间变化;未分解区储层的渗透率为初始渗透率,水合物饱和度为初始值.

(5)不考虑岩石的压缩性,且流体流动满足达西定律,忽略开采过程中水合物的二次生成,忽略储层温度的变化,忽略重力的影响.

2 数学模型

天然气水合物的分解涉及到三相(水相、水合物相、气相)、三组分(水、甲烷、水合物),应用质量守恒方程、分解动力学方程来具体描述水合物的分解过程.

2.1 天然气水合物分解动力学方程

天然气水合物化学反应式为:CH4·6H2O⟺CH4+6H2O(Ahmadi et al., 2007).当水合物稳定存在的条件被打破后,水合物开始分解,其相平衡关系式为(Soloan, 1998)

(1)

式中Pe为水合物平衡压力(Pa),T为温度(K).水合物相平衡示意图(图3),当储层压力低于水合物相平衡曲线,水合物分解产生甲烷和水;当储层压力高于水合物相平衡曲线,水合物不分解以固体形式存在.

图3 水合物相平衡曲线图Fig.3 Hydrate phase equilibrium curve

(2)

式中Mg为甲烷的摩尔质量,取值0.016 kg·mol-1;Pe为水合物平衡压力(Pa);Pg为甲烷气体压力(Pa).kd为反应动力学常数(Pa·s);As为反应比表面积(1/m2).

其中,反应动力学常数(kd)计算公式为

(3)

采用Masuda(2002)提出的反应比表面积(As)计算公式:

As=φShAsi,

(4)

式中Asi为常数,取值3.75×105,Sh为水合物饱和度,φ为孔隙度.

(5)

(6)

式中Mw为水的摩尔质量,取值0.018 kg·mol-1,n为水合数,取值6.

2.2 气、水、水合物质量守恒方程

气质量守恒方程:

(7)

水质量守恒方程:

(8)

水合物质量守恒方程:

(9)

式中t为时间,d(天数);ρg为气体密度,kg·m-3;ρw为水的密度,kg·m-3;ρh为水合物的密度,kg·m-3;Sg为气体饱和度;Sw为水饱和度;Sh为水合物饱和度;vg为气体渗流速度,m·s-1;vw为水渗流速度,m·s-1.其中,气体密度(ρg)采用理想状态方程进行计算:

(10)

根据多相流达西定律公式可计算气体渗流速度(vg)和水渗流速度(vw):

(11)

(12)

式中μg为气黏度,MPa·s;μw为水相黏度,MPa·s,K为含水合物时储层渗透率(md),Kri(i=g、w)为流体相对渗透率(10-3μm2)

采用Masuda(2002)提出的渗透率(K)计算公式:

K=K0(1-Sh)N,

(13)

式中K0为绝对渗透率(md),N为渗透率衰减指数,取值为4.

采用Genuchten(1980)提出的相对渗透率(Kri(i=g、w))计算模型,其公式为:

(14)

(15)

气相压力(Pg)与水相压力(Pw)之间的压差为毛细管力,采用(Genuchten, 1980)提出的毛细管力(Pc)计算公式:

Pc=Pg-Pw,

(16)

(17)

储层孔隙由水合物相、水相、气相共同占据,各相饱和度关系为

Sh+Sw+Sg=1.

(18)

2.3 初始条件和边界条件

Sg=Sg0,Sh=Sh0,Sw=Sw0,t=0

式中Sg0为气相初始饱和度;Sh0为水合物初始饱和度;Sw0为水相初始饱和度;Pout为开采井压力,MPa;l为外边界半径(m).

上述数学模型经过求解后可以得到任意时刻储层压力、水合物饱和度、渗透率的变化规律,具体流程图如图4所示.

图4 模型求解流程图Fig.4 Flow chart of the numerical simulation

3 神狐海域水合物降压开采实例分析以及模型验证

2017年5月10日起中国地质调查局在神狐海域水深1266 m、海底以下203 m至277 m开采出天然气,并成功进行点火(吴时国和王吉亮, 2018; 乔玲茜等, 2018; 程聪等, 2019).至5月18日,神狐海域水合物藏连续产气18 d,累积产气超过12×104m3,其中甲烷含量高达99.5%.至5月26日,连续产气16 d,平均日产气超过1×104m3.至6月10日,连续产气31 d,累积产气超过21×104m3,6月21日,已连续试采42 d,累积产气量超过23.5×104m3.截至7月9日,神狐海域水合物藏连续产气60 d后试采结束,累积产气超过30.9×104m3(李淑霞等, 2018; 刘佳丽, 2018).

应用开发的程序对神狐海域水合物藏降压开采进行了模拟计算,模型中:初始压力6 MPa,水合物初始饱和度40%,初始渗透率为10 md,孔隙度为38%.模拟的产气速率和累积产气量如图5所示.其中,蓝色点线为产气速率,黑色虚线为累积产气量,红色实线为神狐海域试开采过程中累积产气量.模拟结果与实际试采累积产气量变化趋势基本一致.由于模型中没有考虑水合物层厚度以及储层出砂等因素,水合物累积产气量模拟结果略高于实际的产气量.

图5 神狐海域水合物试采模拟结果Fig.5 Simulation results of hydrate pre-production test in Shenhu area

4 数值模拟结果分析

4.1 水合物降压开采动态参数变化规律

模型采用了开放型外边界条件,半径为5000 m;水合物储层的物性参数和初始条件如表1所示.针对这一模型，对降压开采过程中压力、水合物饱和度以及渗透率的变化规律进行分析.

图6为储层压力与开采时间的平面关系图.降压开始后,井眼周围形成压降漏斗,随着开采时间的增加,压降辐射区域逐步变大,在压降漏斗未涉及的区域,压力保持为储层初始压力.

表1 水合物藏模型物性参数Table 1 Physical properties of the hydrate reservoir model

图7为不同开采时间水合物饱和度随径向距离的变化关系.图中可以看出,开采初期(头30天内),水合物饱和度下降的范围较小;到约第90天,井眼附近才出现了水合物饱和度为0的区域(即完全分解区);第90天以后(如第240天),径向上有一个区域其水合物饱和度从0变化到初始值,这个区域即是典型的分解过渡带;随着开采时间的推进,分解过渡带沿半径方向向外移动.

图8为不同开采时间渗透率随径向距离的变化关系.储层绝对渗透率为100×10-3μm2,储层初始渗透率为9.2×10-3μm2;储层渗透率是水合物饱和度的函数,成相反关系,即水合物的分解使储层渗透率增大.水合物分解以井眼为中心随时间而向外推进,渗透率增大的区域也与之相应.

图6 储层压力变化平面图Fig.6 Contour plot of reservoir pressure change

图7 水合物分解过程饱和度变化曲线Fig.7 Saturation variation during hydrate dissociation

图8 水合物分解过程渗透率变化曲线Fig.8 Permeability variation during hydrate dissociation

模型外边界条件的选择对模拟结果影响很大.当模型采用了封闭型外边界条件、半径为225 m时,得到的水合物饱和度随径向距离的变化关系如图9所示.图中可以看出,当分解外沿达到了封闭边界后,外边界附近的水合物饱和度(饱和度最大值)随开采时间推进而开始下降,且最终(1500天时)降至零,即整个封闭区域的水合物全部分解.图9所示的水合物饱和度随径向距离的变化关系在文献中多次出现,但实际水合物开采时外边界条件应该是开放的,因此需要小心使用参考文献所得到的这一计算结果.

图9 外边界封闭时，模拟出的水合物分解饱和度变化曲线Fig.9 Simulated saturation variation during hydrate decomposition when the outer boundary is assumed as closed

4.2 分解过渡带移动规律分析

为进一步表明水合物饱和度随时间的变化关系,绘制了水合物饱和度随时间变化的平面图，如图10所示.图中红色区域代表未分解区,黄色区域代表完全分解区,最外侧黑色圆圈为水合物分解过渡带外沿,最内侧黑色圆圈为水合物分解过渡带内沿.随着开采时间的增加,分解过渡带以井为中心沿径向扩散,经过很长一段开采时间后,井眼周围的水合物才全部分解, 而不是打破水合物相平衡条件后就立即全部分解.

基于上述水合物饱和度的分析,绘制了两种模型外边界条件下分解过渡带随开采时间的变化关系图. 图11为分解过渡带移动距离与开采时间的关系,其中模型外边界为开放型,且半径为5000 m.可以看出:在压力驱动力的作用下,水合物发生分解,分解过渡带沿着半径方向向外移动;分解过渡带宽度随着开采时间逐渐增大,到一定天数后趋于稳定;分解过渡带外沿、内沿移动速度不同步,开采前期,水合物分解过渡带外沿移动速度远远大于内沿移动速度,开采后期,由于水合物分解驱动力的减小,分解过渡带外沿、内沿移动速度都变慢.

图12为模型外边界为封闭型,且半径为225 m.图中可以看出:随着开采时间的推进,分解过渡带内沿和外沿呈指数形式移动,即内沿和外沿移动速度越来越快.分解过渡带外沿推进至模型外边界处之后,过渡带外沿停止在外边界处(225 m),过渡带内沿继续移动,整个封闭区域水合物总量减少,表现出过渡带内沿移动速度快速增大.这一移动规律看似水合物在降压开采过程中会加快地无限分解下去,实际上是因为模拟计算中外边界设置为封闭型,且外边界半径较小所致.类似这种关系图在文献中多次出现,但不符合实际的降压开采过程,因此需要小心使用文献中所得到的这一计算结果,避免误解.

图10 降压开采过程中水合物饱和度随时间的变化Fig.10 Evolution of gas hydrate saturation during depressurization process

图11 水合物分解过渡带移动规律Fig.11 The evolution of the hydrate dissociation transition zone

图12 外边界设置为封闭型且半径较小时，模拟得出的水合物分解过渡带移动规律Fig.12 Simulated evolution of the hydrate dissociation transition zone when the outer boundary is assumed as closed and small in radius

4.3 降压开采参数敏感性分析

目前,降压法是水合物试采工程中最为常用的开采方法,但试采工程中不同降压开采参数对水合物分解过程的影响仍不够明确.本节计算了不同开采井压力、水合物初始饱和度、储层绝对渗透率以及分解动力学常数对水合物分解过程的影响.表2给出了模拟计算分析中各开采参数的取值.

4.3.1 开采井压力的影响

图13是开采井压力分别为5 MPa、7 MPa、9 MPa分解过渡带随开采时间的变化规律.可以看出:井底压力增大时,分解过渡带沿径向方向推进距离变短,过渡带宽度变窄,分解外沿、内沿移动速度变慢.这主要是因为单一的降压开采,水合物分解的驱动力来自于储层和开采井之间的压差.当储层压力一定时,井底压力越低,储层压力和开采井之间的压力梯度越大,水合物分解越快,分解过渡带推进距离越远.在实际的降压开采过程中,建议选用较低的开采井压力以取得良好地降压开采效果.

表2 不同降压开采参数取值Table 2 Parameter values of different depressurization production

图13 不同开采井压力下水合物分解过渡带移动随开采时间变化曲线图Fig.13 Evolution of hydrate dissociation transition zone at different bottom hole pressure

4.3.2 储层绝对渗透率的影响

图14是储层绝对渗透率取值为50 md、75 md、100 md时分解过渡带随开采时间的变化规律.可以看出:储层绝对渗透率越大,井中低压在储层中传播速度越快,水合物的分解越快,表现在水合物分解过渡带推进距离和外沿、内沿移动速度越大.当储层绝对渗透率较小时,可以先通过压裂技术增大储层绝对渗透性,再结合降压法开采水合物.

图14 不同渗透率下水合物分解过渡带移动随开采时间变化曲线图Fig.14 Evolution of hydrate dissociation transition zone at different permeability

图15 不同初始饱和度下水合物分解过渡带移动随开采时间变化曲线图Fig.15 Evolution of hydrate dissociation transition zone at different initial hydrate saturation

4.3.3 水合物初始饱和度的影响

图15是水合物初始饱和度分别为0.35、0.45、0.55时,分解过渡带随开采时间的移动规律.可以看出:随着水合物初始饱和度的增大,分解过渡带沿径向推进距离越近,外沿、内沿移动速度越慢.这主要是因为水合物饱和度较高时,储层中更多的流动空间被固体水合物占据,储层渗流作用减弱,压降传播速度降低,水合物分解速度变慢.

4.3.4 反应动力学常数的影响

天然气水合物的开采与常规能源的开采不同,水合物藏的开采主要受气、水流动特性和水合物分解动力学特性的影响.图16是反应动力学常数分别为3.5×10-11、5.5×10-11、9.5×10-11时分解过渡带随开采时间移动规律.可以看出:分解外沿几乎不受分解动力学常数的影响;分解动力学常数对分解过渡带内沿影响较大,分解动力学常数越大,分解过渡带内沿移动距离越远.反应动力学常数与温度正相关,降压开采可以结合注热法进一步促进水合物的分解.

图16 不同分解动力学常数下水合物分解过渡带移动随开采时间变化曲线图Fig.16 Evolution of hydrate dissociation transition zone at different decomposition kinetic constant

5 结论

为揭示天然气水合物降压开采过程中水合物分解规律,将水合物储层分为完全分解区、分解过渡带和未分解区,建立了柱坐标系下水合物降压开采的数学模型,应用有限差分法对其进行求解,并以神狐海域试采数据进行了验证,进而分析了水合物分解过程中压力、水合物饱和度、渗透率的变化规律,以及水合物分解过渡带移动规律,研究了不同降压开采参数对水合物分解过程的影响.得出以下结论:

(1)降压开始后,开采井周围迅速形成压降漏斗,随着开采时间的增加,压降漏斗向储层远处扩散,在压降漏斗未涉及的区域,压力保持为储层初始压力;

(2)井眼处的水合物最先分解,饱和度随着开采时间的推进逐渐降低,储层渗透率逐渐增大,饱和度降低的区域沿径向向外扩散,渗透率增大的区域也与之相应;

(3)水合物分解过渡带沿着半径方向向外移动,分解过渡带外沿、内沿移动速度不同步,开采前期,外沿移动速度远远大于内沿移动速度,开采后期,外沿、内沿移动速度都变慢,分解过渡带宽度随着开采时间逐渐增大,到一定天数后趋于稳定;

(4)水合物降压开采的主要控制参数包括开采井压力、水合物初始饱和度、储层绝对渗透率、水合物分解动力学常数等;

(5)模拟水合物降压开采时,如果选择封闭型边界且半径较小,则所得出的模拟结果与实际开采情况会有较大的差别.