天然气水合物降压分解诱发储层变形破坏正交数值模拟实验研究

2022-03-10孙可明

特种油气藏 2022年1期

翟诚，孙可明

(1.山西工程技术学院，山西阳泉 045000；2.辽宁工业大学，辽宁锦州 121001；3.青岛理工大学，山东青岛 266520)

0 引言

天然气水合物因具有储量丰富、清洁、高效等优点，被认为是极具开发价值和良好发展前景的新型替代能源之一。勘探成果表明，中国的天然气水合物资源主要蕴藏在南海北部的神狐海域、青藏高原陆地冻土区以及漠河盆地陆地冻土区[1-3]。降压法因具有操作简单、经济、高效等优点，被认为是未来工业化开采水合物的重要方法之一。但如果防控措施不当，水合物降压分解引起储层失稳破坏后，会导致包括海底边坡失稳等一系列灾害效应的发生[4-5]。目前，很多国内外学者运用数值模拟方法分析天然气水合物分解引起的储层变形破坏问题，并取得了一些具有借鉴意义的成果[6-14]。但在降压分解条件下，缺少影响水合物沉积层近井储层变形破坏因素的敏感性分析。鉴于此，以中国南海神狐海域天然气水合物降压开采为工程背景，采用数值模拟与正交实验相结合的方法，基于ABAQUS有限元软件的USDFLD子程序进行二次开发，探讨初始水合物饱和度、井底压力和水合物藏所处地应力状态对水合物沉积层近井储层变形破坏影响的敏感程度。

1 模型的建立

建立模型时，做出如下假设：①将水合物沉积层看作含有固、液、气三相的多孔介质，其中，水合物与沉积物骨架胶结构成多孔介质的固体骨架，甲烷气体和水则充填在固体骨架围成的孔隙中；②固体骨架为各向同性的弹塑性材料，其变形满足小变形假设，且在变形过程中沉积物骨架的密度和水合物的密度不发生变化；③不考虑水合物降压分解过程中储层温度的变化，即温度对水合物沉积层有效应力和孔隙流体性质的影响忽略不计；④甲烷气体为理想气体，水合物沉积层中水、气的渗流满足Darcy定律。

1.1 水合物分解动力学方程

水合物分解的产气速率采用Kim等[15]建立的分解动力学方程来计算，即：

(1)

As=φShAhs

(2)

lgpe=0.0342(T-273.15)+0.0005(T-273.15)2+6.4804

(3)

根据式(1)和化学反应方程式CH4·5.75H2O=CH4+5.75H2O，可得单位体积水合物的分解速率和产水速率分别为：

(4)

1.2 变形场方程

变形场方程包括固体骨架的静力平衡方程、几何方程和弹塑性本构方程3个部分。根据有效应力原理和弹性力学理论，水合物沉积层固体骨架的静力平衡方程可用张量形式表示为：

(5)

(6)

孔隙中气相压力与水相压力之间的关系满足如下的方程：

pg-pw=pc

(7)

式中：pc为毛管力，Pa。

采用VG模型计算毛管力，其表达式为：

(8)

式中：Swr为残余水饱和度；Swmax为最大含水饱和度；a、m、n为VG模型参数，且m=1-1/n。

几何方程为：

(9)

式中：εij为应变张量；u为位移，m；ui，j、uj，i为采用张量表示的不同方向位移的偏微分形式。

固体骨架增量形式的弹塑性本构方程为：

[dσ′]=[De][dεe]=
[De]{[dε]-[dεP]-[dεpl]}

(10)

式中：[dσ′]为有效应力增量矩阵；[De]为弹性模量矩阵；[dεe]为弹性应变增量矩阵；[dε]为总应变增量矩阵；[dεP]为孔隙压力产生的压应变增量矩阵；[dεpl]为塑性应变增量矩阵。

孔隙压力变化产生的单向应变增量为：

(11)

式中：Ksh为固体骨架的体积模量，Pa；vsh为固体骨架的泊松比；Esh为固体骨架的弹性模量，Pa。

塑性应变增量为：

(12)

(13)

1.3 渗流场方程

考虑水合物分解效应和固体骨架变形对孔隙流体渗流的影响，根据质量守恒方程和达西定律，可推导得到气、水两相的流固耦合渗流场方程。

气相：

(14)

孔隙中气相相对渗透率为[18]：

(15)

水相：

(16)

孔隙中水相的相对渗透率为[19]：

(17)

水合物相的质量守恒方程：

(18)

将上述数学模型与水合物沉积层弹性模量和渗透率随有效应力及水合物饱和度变化的关系模型结合[20]，再给定具体的初始条件和边界条件，即可构成水合物降压分解诱发储层变形破坏的流固耦合弹塑性模型。

2 正交数值模拟实验

2.1 有限元模型、边界条件及模拟参数

建立如图1a所示的有限元模型进行正交数值模拟实验，并对井眼附近区域进行网格细化，坐标如图1b所示。整体模型的几何尺寸为：边长L=10 m，井眼半径r=0.2 m。地质条件参考中国南海神狐海域水合物藏进行选取：海水深度为1 235 m，水合物藏距离海床的深度为182 m，储层沉积物骨架主要由粉砂岩和黏土构成。模拟所用基本参数如表1所示[7，10，13]。

变形场边界条件如下。

(1)σBC=σhmax，σCD=σhmin，(σBC、σCD分别为作用在BC、CD边的应力，Pa；σhmax、σhmin分别为最大和最小水平地应力，Pa)。

(2)UAB=Uy=0，UDE=Ux=0，AE为固定端约束(UAB、UDE为AB、DE边位移，m；Ux、Uy为x、y方向位移，m)。

渗流场边界条件：

(2)AB、DE为渗流封闭边界。

图1 有限元模型示意图

表1 数值模拟基本参数

2.2 正交数值模拟实验因素的确定与实验结果

表2 正交数值模拟实验因素水平

表3 正交数值模拟实验方案及结果

3 正交数值模拟实验结果分析

3.1 实验结果的极差分析

根据表3的实验结果可分别计算得到以PEEQmax和塑性区范围为评价指标的各因素水平的均值(表4、5)，并进一步采用直观分析法对表4、5进行极差分析，从而确定每个因素对近井储层变形破坏影响的敏感程度。由表6、7的极差数据分析可知：各因素对近井储层的影响程度由大到小依次为井底压力、有效主应力差、初始水合物饱和度，说明井底压力对近井储层变形破坏的影响最大，有效主应力差的影响次之，而初始水合物饱和度的影响最小；各因素对近井储层塑性区范围的影响程度由大到小依次为有效主应力差、井底压力、初始水合物饱和度，说明有效主应力差对近井储层塑性区范围的影响最大，井底压力的影响次之，而初始水合物饱和度的影响最小。

表4 以PEEQmax为评价指标的各因素水平的实验均值

表5 以塑性区范围为评价指标的各因素水平的实验均值

表6 PEEQmax的直观分析

表7 塑性区范围的直观分析

为了进一步说明各因素对近井储层变形破坏和塑性区范围影响的变化趋势，根据表4、5的分析结果，分别作出各影响因素对近井储层影响的直观分析图(图2、3)。

图2 各因素对近井储层PEEQmax的影响规律

图3 各因素对近井储层塑性区范围的影响规律

由表6、7和图2、3可得到如下结论。

(1) 近井储层发生变形破坏的程度主要受井底压力和有效主应力差的影响。其中，井底压力对PEEQmax的影响最为显著，而塑性区范围则受有效主应力差的影响最为明显。相比较于前2个参数，初始水合物饱和度对近井储层变形破坏的影响相对较小。

(2) 随着初始水合物饱和度的增大，近井储层PEEQmax和塑性区范围逐渐变小，即发生变形破坏的程度逐渐减弱。这是因为，初始水合物的含量越多，其与沉积物骨架胶结后形成的储层固体骨架的力学性质越好，抵抗变形破坏的能力也就越强。

(3) 井底压力越低，近井储层PEEQmax和塑性区的范围越大。这是因为，井底压力越低，水合物降压分解的效果越明显，近井储层水合物分解区的范围以及分解区储层的弹性模量、黏聚力以及摩擦角等力学参数下降的幅度越大，其发生变形破坏的程度也就越严重。虽然井底压力的降低，会使近井储层的有效应力增大，在一定程度上抑制和削弱水合物饱和度降低引起的储层力学参数劣化的程度，但这与水合物分解带来的负面效应相比，后者占据主导地位。因此，近井储层的抗剪强度总体仍表现出下降趋势，并最终发生剪切破坏。

(4) 随着有效主应力差的增大，近井储层PEEQmax和塑性区的范围越来越大。根据Mohr-Coulomb准则，在其他条件相同的前提下，近井储层的有效主应力差越大，发生剪切破坏的程度越严重。另外，受水平地应力非均匀性的影响，近井储层塑性区范围的分布表现出明显的各向异性，主要沿最大水平地应力方向扩展。图4为选取的其中3组实验的近井储层塑性区分布等值线图以及沿储层AB方向PEEQ分布的对比分析图，由图4可明显看出塑性区沿最大水平地应力方向的扩展更为显著。

图4 近井储层塑性区分布等值线及对比分析

3.2 实验结果的方差分析

比较实验因素水平变化和误差波动对实验结果的差异，将表3中的“空白列”作为误差列对数值模拟的实验结果进行方差分析，结果如表8、9所示。由表8、9可知：影响近井储层PEEQmax和塑性区范围的3个因素中，井底压力和有效主应力差2个影响因素对应的P值皆小于0.05，说明它们是影响近井储层变形破坏的显著因素。根据F值的大小进行排序，3个因素对近井储层PEEQmax的影响程度由大到小依次为井底压力、有效主应力差、初始水合物饱和度；对塑性区范围的影响程度由大到小依次为有效主应力差、井底压力、初始水合物饱和度。这与极差分析所得结论是一致的。