-
天然气水合物开采过程是一个多相多组分非等温的复杂过程,涉及水合物相变分解、流体运移、热量传递等物理过程。为便于分析计算,理论模型作以下假设:(1)仅考虑三相:气相、液相以及水合物固相,水合物分解产生的气体为甲烷气体,且不考虑甲烷气体在水中的溶解;(2)多孔介质基质中为气、液两相渗流,流动规律符合达西定律且多孔介质具有均质性、各向同性;(3)裂隙处水合物分解瞬间转换为气体和水。
-
天然气水合物生成和分解的化学过程:
$$ {{\mathrm{CH}}_4}\cdot n{{\mathrm{H}}_2}{\mathrm{O}} \leftrightharpoons {{\mathrm{CH}}_4} + n{{\mathrm{H}}_2}{\mathrm{O}} $$ (1) 式中:
n——水合数取6。
天然气水合物分解生成甲烷气体和水,论文中下标h、g、w分别表示水合物固相、气相和液相。分解动力学方程采用Kim-Bishnoi模型[13],气体和水的生成量由化学反应方程控制,水合物动力学反应速率表示为:
$$ m_{\mathrm{h}} = k_{\mathrm{d}}M_{\mathrm{h}}A_{\mathrm{s}}(P_{\mathrm{e}} - P_{\mathrm{g}}) $$ (2) 式中:
kd ——动力学反应速率(Pa·s);
Mh ——水合物摩尔质量(kg/mol);
Pe ——水合物相平衡压力(kPa),表示为${P_{\mathrm{e}}} = \exp (a - {b \mathord{\left/ {\vphantom {b T}} \right. } T})$,a和b为回归系数取值为39.08和8 533,T为不同时刻下的温度(K)。
$$ {k_{\mathrm{d}}} = k_{\mathrm{d}}^0{{\mathrm{exp}}{\left( { - \tfrac{{\Delta {E_a}}}{{{\mathrm{R}}T}}} \right)}} $$ (3) R ——理想气体常数,一般为8.314 J/(mol·K)。
$$ A_{\mathrm{s}} = \varphi S_h\sqrt {\frac{{{\varphi ^3}{{(1 - S_{\mathrm{h}})}^3}}}{{2k}}} $$ (4) 式中:
As ——水合物反应比表面积(1/m),比表面积的大小与水合物储层的孔隙度、绝对渗透率,水合物饱和度有关[14-15];
${S_{\mathrm{h}}}$ ——水合物饱和度;
k ——绝对渗透率,遵从Masuda[16]渗透率下降模型。
$$ k = {k_0}{(1 - {S_{\mathrm{h}}})^N} $$ (5) 式中:
$ {k}_{0} $ ——无水合物多孔介质初始绝对渗透率(m2);
N ——渗透率下降指数。
水合物分解过程中的质量守恒可表示为:
$$ \frac{\mathrm{\partial}(\varphi\rho_{\mathrm{h}}S_{\mathrm{h}})}{\mathrm{\partial}t}=m_{\mathrm{h}} $$ (6) 式中:
φ ——孔隙度;
$ {\rho }_{h} $ ——水合物密度。
-
水合物分解过程气体和水的渗流可以表示为[17]:
$$ {\rho }_{{\mathrm{w}}}{C}_{{\mathrm{p}},{\mathrm{w}}}\frac{\partial {P}_{{\mathrm{w}}}}{\partial t}+\nabla \cdot{\rho }_{w} \boldsymbol{v}_{{\mathrm{w}}}={m}_{{\mathrm{w}}}+{\rho }_{{\mathrm{w}}}{C}_{{\mathrm{p}},{\mathrm{w}}}\frac{\partial {P}_{{\mathrm{g}}}}{\partial t} $$ (7) $$ {\rho }_{{\mathrm{g}}}S\frac{\partial {P}_{{\mathrm{g}}}}{\partial t}+\nabla \cdot{\rho }_{{\mathrm{w}}} \boldsymbol {v}_{{\mathrm{g}}}={m}_{{\mathrm{g}}}+{\rho }_{{\mathrm{g}}}{C}_{{\mathrm{p}},{\mathrm{w}}}\frac{\partial {P}_{{\mathrm{w}}}}{\partial t} $$ (8) $$ {C_{{\mathrm{p}},{\mathrm{w}}}} = - {\varphi _{{\mathrm{wg}}}}\frac{{{\mathrm{d}}{S_{\mathrm{w}}}}}{{{\mathrm{d}}{P_{\mathrm{c}}}}} $$ (9) $$ {\varphi _{{\mathrm{wg}}}} = \varphi (1 - {S_{\mathrm{h}}}) $$ (10) $$ {P_{\mathrm{c}}} = {P_{\mathrm{g}}} - {P_{\mathrm{w}}} $$ (11) $$ {\chi _{\mathrm{f}}} = \frac{1}{{{\rho _{\mathrm{g}}}}}\frac{{\partial {\rho _{\mathrm{g}}}}}{{\partial {P_{\mathrm{g}}}}} $$ (12) $$ S = {\varphi _{{\mathrm{wg}}}}{S_{\mathrm{g}}}{\chi _f} + {C_{{\mathrm{p}},{\mathrm{w}}}} $$ (13) 式中:
$ {C}_{{\mathrm{p}},{\mathrm{w}}} $ ——毛细管容量系数[18]。
$ {\chi }_{{\mathrm{f}}} $ ——甲烷气体的压缩系数。
S ——储水系数。
式中气体和水流速由达西定律表示。
$$ {v_{\mathrm{i}}} = - \frac{{k{k_{{\mathrm{r}},{\mathrm{i}}}}}}{{{\upsilon _{\mathrm{i}}}}}\nabla {\rho_{\mathrm{i}}} $$ (14) 式中:
i ——气相(g)和水相(w);
$ {\nu }_{\mathrm{i}} $ ——i相的速度;
$ {\upsilon}_\mathrm{i} $ ——i相的动力粘度(Pa·s);
$ {P}_\mathrm{i} $ ——i相的压力(Pa);
$ {k}_{\mathrm{ri}} $ ——i相的相对渗透率;
$ {\rho }_\mathrm{i} $ ——i相的密度。
-
气体和水的相对渗透率采用Corey模型[19]来描述:
$$ {k_{{\mathrm{rw}}}} = {\left( {\dfrac{{\dfrac{{{S_{\mathrm{w}}}}}{{{S_{\mathrm{g}}} + {S_{\mathrm{w}}}}} - {S_{{\mathrm{wr}}}}}}{{1 - {S_{{\mathrm{wr}}}} - {S_{{\mathrm{gr}}}}}}} \right)^{{n_{\mathrm{w}}}}} $$ (15) $$ {k_{{\mathrm{rg}}}} = {\left( {\dfrac{{\dfrac{{{S_{\mathrm{g}}}}}{{{S_{\mathrm{g}}} + {S_{\mathrm{w}}}}} - {S_{{\mathrm{gr}}}}}}{{1 - {S_{{\mathrm{wr}}}} - {S_{{\mathrm{gr}}}}}}} \right)^{{n_{\mathrm{g}}}}} $$ (16) $$ {P_{\mathrm{c}}} = {P_{\mathrm{g}}} - {P_{\mathrm{w}}} = {P_{{\mathrm{ce}}}}{\left( {\dfrac{{\dfrac{{{S_{\mathrm{w}}}}}{{{S_{\mathrm{g}}} + {S_{\mathrm{w}}}}} - {S_{{\mathrm{wr}}}}}}{{1 - {S_{{\mathrm{wr}}}} - {S_{{\mathrm{gr}}}}}}} \right)^n} $$ (17) 式中:
$ {n}_{{\mathrm{w}}} $、$ {n}_{{\mathrm{g}}} $ ——模型参数,分别取2、4;
n ——模型参数取2;
$ {S}_{{\mathrm{gr}}}、{S}_{{\mathrm{wr}}} $ ——束缚气饱和度、束缚水饱和度;
$ {S}_{{\mathrm{g}}} $、$ {S}_{{\mathrm{w}}} $ ——气体饱和度、水饱和度,满足${S_{\mathrm{g}}} + {S_{\mathrm{w}}} +{S_{\mathrm{h}}} = 1$。
-
固相水合物分解生成气体和水的过程需要不断吸收热量,水合物储层中热量交换包括反应热、热对流与热传导,水合物储层基质与裂隙流体的热量交换通过热对流和热传导进行。水合物分解过程的能量转移采用能量方程进行描述。
$$ {C}_{{\mathrm{eq}}}\frac{\partial T}{\partial t}+\left({\rho }_{{\mathrm{g}}}{C}_{{\mathrm{g}}}\boldsymbol{v}_{{\mathrm{g}}}+{\rho }_{{\mathrm{w}}}{C}_{{\mathrm{w}}}\boldsymbol{v}_{{\mathrm{w}}}\right)\cdot \nabla T+\nabla \cdot \left(-{k}_{{\mathrm{eq}}}\nabla T\right)={Q}_{{\mathrm{h}}}+q $$ (18) $$ {C_{{\mathrm{eq}}}} = (\varphi ({\rho _{\mathrm{w}}}{S_{\mathrm{w}}}{C_{\mathrm{w}}} + {\rho _{\mathrm{g}}}{S_{\mathrm{g}}}{C_{\mathrm{g}}} + {\rho _{\mathrm{h}}}{S_{\mathrm{h}}}{C_{\mathrm{h}}}) + (1 - \varphi ){\rho _{\mathrm{s}}}{C_{\mathrm{s}}}) $$ (19) $$ {k_{{\mathrm{eq}}}} = \varphi {S_{\mathrm{g}}}{\lambda _{\mathrm{g}}} + \varphi {S_{\mathrm{w}}}{\lambda _{\mathrm{w}}} + \varphi {S_{\mathrm{h}}}{\lambda _{\mathrm{h}}} + (1 - \varphi ){\lambda _{\mathrm{s}}} $$ (20) $$ {Q_{\mathrm{h}}} = \frac{{{m_{\mathrm{h}}}}}{{{M_{\mathrm{h}}}}}\left( {{C_1} + {C_2}T} \right) $$ (21) $$ q = h({T_{\mathrm{r}}} - {T_{\mathrm{w}}}) $$ (22) 式中:
$ {\rho }_{i} $ ——气体(g)、水(w)、水合物(h)、固体骨架颗粒(s)的密度(kg/m3);
$ {S}_{i} $ ——气体、水、水合物饱和度;
$ {C}_{i} $ ——气体、水、水合物、固体骨架颗粒的比热容[J/(kg·K)];
$ {k}_{{\mathrm{eq}}} $ ——等效导热系数[W/(m·K)];
$ {\lambda }_{i} $ ——气体、水、水合物、固体骨架颗粒的导热系数[W/(m·K)];
h ——对流换热系数[W/(m2·K)];
$ {Q}_{{\mathrm{h}}} $ ——水合物相变潜热(J);
$ {C}_{1} $、$ {C}_{2} $ ——回归系数;
q ——裂隙流体与水合物储层接触界面热通量(W/m2);
Tr、Tw ——水合物储层温度、裂隙流体温度(K)。
-
基于Comsol Multiphysics数值模拟软件PDE模块建立以水合物饱和度、孔隙气压力、孔隙水压力、储层温度为场变量的THC多物理场耦合数值模型并进行求解,求解过程如图1所示。为验证理论模型与数值模型的正确性,分别与Masuda[16]砂岩岩芯降压分解实验结果、Wilder[20]天然气水合物降压分解模型模拟结果进行对比验证。
Masuda砂岩岩芯降压分解实验中,圆柱形岩芯样品直径为5.1 cm,长度为30 cm,A、B、C三个监测点分别距离降压边界0.375 cm、15 cm、22.5 cm。模型如图2所示,左侧边界为降压边界,P=2.84 MPa,其他边界为恒温和无流动边界,模型初始条件如表1所示。
参数 数值 初始气体压力Pg0/MPa 3.75 初始水合物饱和度Sh0 0.443 初始水饱和度Sw0 0.206 水合物储层初始温度T/K 275.45 孔隙度φ 0.182 无水合物多孔介质绝对渗透率k0/mD 98 Table 1. Parameters for Masuda hydrate decomposition by depressurization experiment
图3为测点温度变化和产气量的结果对比图。由图3(a)可以看出,模拟和实验测点的温度基本吻合,整体上处于先降低后升高的趋势,距离降压边界近的测点温度先降至最低点然后温度回升。这是由于靠近降压边界处的储层压力先降低至水合物相平衡压力之下,水合物大量分解并消耗大量热量,外界补充热量小于水合物分解消耗的热量,储层温度降低。随后水合物分解速度减慢,外界补充热量大于水合物分解消耗的热量,储层温度逐渐升高,上升至环境温度后保持不变。由图3(b)可知,Masuda实验中最终产气量为
8985.7 cm3,产气初期产气速率较高,随着运行时间的变化产气速率逐渐减小进入稳定产气阶段直至产气结束,累积产气量趋于稳定。模型设置左侧降压边界处压力瞬间降低至2.84 MPa,而实验中通过抽取水合物储层流体致使储层压力缓慢降低至2.84 MPa,文章模型中压降传递速率更快,水合物分解速率加快,因此,本文模型中产气速率高于Masuda实验结果。模型数值模拟结果的最终产气量为8990 cm3与Masuda实验最终产气量基本相同。图4为Wilder水合物降压分解模型,模型尺寸为L×W=1.5 m×1 m,左端为降压端P=2.8 MPa,并且保持温度T=6 ℃,其他边界为绝热和无流动边界,左侧降压驱动力作用下水合物分解。模型初始条件及基本参数如表2所示。
参数 数值 初始气体压力Pg0/MPa 8 初始水合物饱和度Sh0 0.5 初始水饱和度Sw0 0.5 水合物储层初始温度T/K 279.15 孔隙度φ 0.2 无水合物多孔介质绝对渗透率k0/mD 300 残余水饱和度Swr 0.1 残余气饱和度Sgr 0.01 Table 2. Parameters for Wilder decomposition by depressurization model
图5为数值模拟结果与其他数值分析结果对比图,可以看出,所建立模型的计算结果与其他数值分析结果一致。降压3 d后,近降压端水合物饱和度由0.5降低至0,远离降压端处水合物饱和度由0.5降至0.4后保持不变,这是由于近降压端模拟域温度高,给水合物分解提供足够的热量,水合物分解较快,离降压端较远处模拟域温度较低,水合物分解所需热量供给不足,水合物停止分解,如图5(b)所示。
计算结果与Masuda砂岩岩芯水合物降压分解试验累积产气量、测点温度以及Wilder水合物降压分解的模拟结果进行对比分析,所有计算结果吻合较好,验证了所建模型的有效性和可靠性,可用于进一步研究。
-
压裂使低渗透率水合物储层形成裂缝网络,随后通过井筒向水合物储层注入高温水,高温水通过裂缝网络到达储层深处为水合物分解提供热量,这是一种提高储层产气速率的方法,方案示意图如图6所示。单裂隙是组成裂缝网络的基本单元,探究单裂隙水合物储层注热开采产气速率的影响机理为后续研究压裂后含复杂裂隙水合物储层注热开采产气速率的影响机理奠定基础。
基于验证后Wilder水合物储层降压开采数值模型构建单裂隙水合物储层注热开采数值模型,如图7所示,模型尺寸为L×W=1.5 m×1 m,横向贯穿一条开度为e的裂隙,裂隙流体流场与温度场分别用达西定律与流体传热物理场表示。模拟域的初始温度和初始压力分别为279.15 K和8 MPa,模型左侧边界压力为定值,大小与水合物储层初始压力相同,其余三边为热绝缘边界和无流动边界。裂隙注入温度为Tin,裂隙上下边界为热通量边界,裂隙入口和出口处为固定的压强差ΔP。数值模拟具体参数如表3所示。
参数 参考值 参数 参考值 孔隙度φ 0.2 裂缝孔隙度Φf 0.2 初始温度T0/K 279.15 裂隙初始温度Tf0/K 279.15 初始气体压力Pg0/MPa 8 裂缝初始压力Pf0/atm 1 无水合物多孔介质绝对渗透率k0/m2 2.96×10-13 裂缝绝对渗透率kf0/m2 10-10 初始水合物饱和度Sh0 0.5 渗透率下降指数N 9 初始水饱和度Sw0 0.5 相对渗透率指数nw、ng 2、4 进气值参考压力Pce/MPa 0.1 对流换热系数h/[W·(m2·K)-1] 30 Table 3. Calculation parameters for numerical model of heat injection extraction in hydrate reservoirs with a single fracture
-
水合物分解是吸热反应,高温水沿裂缝到达储层深处,一方面高温流体与水合物储层间进行对流换热,为水合物分解提供充足的热量;另一方面,沿裂缝分布的水合物在高温水作用下分解,分解生成的甲烷气体和水在相邻结构单元压力差作用下流出,沿裂缝处储层渗透率增大。图8为注热温度20 ℃,裂隙开度e=1 mm,运行2、10、60、120 h后的水合物储层水合物饱和度分布图。由图8可知,靠近裂缝区域的水合物先分解,沿裂缝形成无水合物锥形区域,随着运行时间的增加,无水合物区域逐渐增大,在整个运行过程中,无水合物区域沿着裂隙中心线对称。这是因为裂隙注热温度与水合物储层温度温差较大,受热水影响区域裂隙流体与水合物储层之间产生剧烈的对流换热效应,靠近裂隙区域水合物储层温度迅速升高,水合物储层压力低于升温后水合物相平衡压力,水合物先分解。
图9为不同注热温度对单裂隙水合物储层产气速率和产气量的影响。从图9(a)可以看出,裂隙注热时产气速率分解可分为3个阶段:(1)第1阶段:产气速率骤升后快速回落到较低水平,产生第一个产气速率峰值。这是因为靠近注热裂隙的水合物储层与裂隙流体对流换热剧烈,靠近注热裂隙部分水合物储层升温快,水合物相平衡条件破坏,水合物快速分解。但由于靠近裂隙部分水合物完全分解,离裂隙较远处水合物未达到分解条件,产气速率回落到较低水平;(2)第2阶段:产气速率逐渐上升,达到第二个峰值产气速率后逐渐下降。在裂隙对流换热和水合物储层内部热传导的作用下,离裂隙较远处水合物储层温度逐渐升高,该区域水合物大量分解,产气速率逐渐上升。由于水合物大量分解吸收热量多,水合物储层温度上升速度减慢,水合物分解速率随之逐渐降低;(3)第3阶段:储层内水合物饱和度不断减小,产气速率骤降至0后保持不变。
Figure 9. Variation of gas production rate and gas production volume with time at different heat injection temperatures
由图9(a)可知,注热温度较高时(T=24、26、28、30 ℃),第2阶段产气速率上升时间较短,出现峰值产气速率的时间更早,且第2阶段峰值产气速率小于第1阶段峰值产气速率。注热温度较低时(T=16、18、20、22 ℃),第2阶段产气速率上升时间较长且第2阶段峰值产气速率大于第1阶段峰值产气速率。这是因为水合物储层初始温度相同,注水温度越高,水合物储层与裂隙流体单位时间交换的热量越多,靠近裂隙部分水合物储层温度单位时间内上升幅度越大,第1阶段峰值产气速率越大。随后,由于水合物储层内部的热传导作用,远离注热裂隙区域水合物储层温度在更短时间内提高,水合物大量分解,产气速率迅速增加。由于温度较高时,第1阶段的对流换热作用大于第2阶段水合物储层热传导作用,第1阶段峰值产气速率大于第2阶段峰值产气速率。随着运行时间的增加,水合物分解产气量不断增大,水合物分解完成后产气量保持不变。并且,随着裂隙注热温度的提高,产气量曲线出现拐点的时间不断提前,越早达到最终产气量,如图9(b)所示。这是因为随着裂隙注热温度的提高,裂隙流体与水合物储层温差越大,单位时间交换的热量越多(见式(21)),水合物分解得越快。
-
图10为注热温度为20 ℃、注热压差为0.05 MPa时不同裂隙开度对水合物储层产气速率和产气量的影响。如图10(a)所示,不同裂隙开度下产气速率随运行时间的变化分为3个阶段,与图7(a)相同,但裂隙开度为1 mm、2 mm时,第1阶段并未出现峰值产气速率。这是因为注热温度较低、裂隙开度较小时,对流换热效果不显著,靠近裂隙部分水合物分解较慢,随着运行时间的增加,水合物储层升温区域逐渐扩大,产气速率逐渐增大。不同裂隙开度下,第1阶段峰值产气速率均低于第2阶段峰值产气速率,由前文可知,这主要由注热温度起主导作用。随着裂隙开度增加,第1阶段峰值产气速率与第2阶段峰值产气速率均有所增加并且峰值产气速率出现的时间不断提前,如图11所示。并且裂隙开度越大,产气量曲线越先出现拐点,越早达到最终产气量,如图10(b)所示。这是因为裂隙开度越大,可用于对流换热的流体越多,水合物储层温度升温更快,水合物分解越快。但是裂隙开度对产气速率的影响有一个阈值,超过该阈值时变化不大,如图11所示,峰值产气速率及峰值产气速率出现的时间随裂隙开度的变化逐渐趋于平缓。
-
图12为不同注热速度下水合物储层产气速率及产气量随运行时间的变化。由图12(a)可知,产气速率随运行时间的变化分为3个阶段,这与图9(a)、图8(a)产气速率随运行时间变化的规律相同。注热速度较大时,初始注热阶段产气速率短暂升高后降低,这是由于裂隙注热速度越大,单位时间内裂隙流体与水合物储层交换的热量越多,靠近裂隙区域水合物储层温度迅速升高,水合物分解较快,但由于水合物分解吸收热量较多,离裂隙较远处水合物储层温度升温速度减慢,产气速率降低。图12(b)为不同注热速度下产气量随运行时间的变化,由图12(b)可以看出,产气量随运行时间逐渐增加,水合物完全分解后,产气量保持不变。
-
基于正交试验设计分析,可确定注热温度、裂隙开度、注热速度对水合物储层峰值产气速率的影响程度,将3个影响因素的水平数均设置为3,正交设计所用具体参数如表4所示。
水平 1 2 3 Tin/℃ 20 24 28 e/mm 1 2 3 vin/(m·s−1) 0.009 26 0.027 78 0.055 56 Table 4. Level of heat injection parameters
基于单裂隙水合物储层注热开采影响因素水平数,设计L9(33)正交表,该正交表共包含9个方案,设置峰值产气速率(qmax)为综合评价指标,对正交试验设计计算结果进行极差分析,如表5所示。
方案 Tin/℃ e/mm vin/(m·s−1) qmax/(ml·min−1) 1 20 1 0.009 26 36.674 2 20 2 0.027 78 42.865 3 20 3 0.055 56 42.264 4 24 1 0.027 78 48.105 5 24 2 0.055 56 56.277 6 24 3 0.009 26 51.993 7 28 1 0.055 56 71.670 8 28 2 0.009 26 62.262 9 28 3 0.027 78 70.434 平均数1 40.601 52.150 50.310 - 平均数2 52.125 53.801 53.801 - 平均数3 68.122 54.897 56.737 - 极差 27.521 2.747 6.427 - 较优水平 28 3 0.055 56 - 因素主次 1 3 2 - Table 5. Orthogonal test scheme for heat injection parameters
极差分析是考虑某个因素对结果的影响时,其他因素对结果的影响是均衡的。以峰值产气速率为指标评价注热参数的影响程度,由表5可以看出,极差最大的因素是注热温度,其最佳水平为28 ℃;注热速度次之,最佳水平为
0.05556 m/s;裂隙开度的极差最小,最佳水平为3 mm。通过以上分析的选取参数范围内的最优组合为:注热温度28 ℃,注热速度0.05556 m/s,裂隙开度3 mm,该组合下的峰值产气速率为71.670 ml/min。
Numerical Simulation of Heat Injection Extraction in Fractured Natural Gas Hydrate Reservoir
doi: 10.16516/j.ceec.2024-120
- Received Date: 2024-04-24
- Accepted Date: 2024-05-20
- Rev Recd Date: 2024-05-15
-
Key words:
- natural gas hydrate /
- heat injection extraction /
- fracture /
- parameter sensitivity analysis /
- gas production rate
Abstract:
Citation: | XUE Luanluan, PENG Rui. Numerical Simulation of Heat Injection Extraction in Fractured Natural Gas Hydrate Reservoir[J]. SOUTHERN ENERGY CONSTRUCTION. doi: 10.16516/j.ceec.2024-120 |