海洋天然气水合物降压开采地层井壁力学稳定性分析1)

时间：2024-05-22

袁益龙许天福辛欣夏盈莉李冰

*(吉林大学地下水资源与环境教育部重点实验室,长春 130021)

†(吉林大学新能源与环境学院,长春 130021)

**(吉林大学建设工程学院,长春 130026)

引言

海洋天然气水合物(简称水合物)通常形成于深水浅覆盖层的松散沉积物中,其成岩性差、抗剪强度低.水合物以固体形态赋存于这些沉积层的孔隙中起到了有效的胶结作用,改变了地层的结构特性[1-2].降压开采水合物使其分解导致地层强度和承载能力大幅度降低,进而弓起生产井周围地层变形,甚至弓发井壁失稳、海底滑坡、海底面沉降等地质灾害[3-4].由此可见,水合物降压开采涉及到松散沉积层系统中热对流传导(T)、水动力(H)、力学(M)过程以及三者之间的相互耦合作用.然而,以往有关水合物开采的研究主要集中于产能评价,往往忽略对水合物地层和井壁稳定性的评估[5-8].

目前国内外有关水合物开采地层井壁稳定的模型研究尚处于起步阶段,能够同时考虑水合物开采THM 耦合过程的数值计算软件不多.Rutqvist 等[3]通过耦合TOUGH+Hydrate 和力学分析软件FLAC3D形成水合物开采THM 耦合软件,并成功应用到冻土区[3]和海洋[2]水合物开采地层力学稳定性分析.Kimoto 等[9]基于化学-热力学-岩土力学耦合分析,建立了预测甲烷水合物分解弓发的地层沉降模型.Lin等[10]基于印度国家水合物项目第二航次(NGHP-02)16 站位的储层资料,分析了砂泥互层结构水合物储层降压开采的地层力学稳定性.Qiu 等[11]通过MH21-HYDRES 和3D FE 地质力学模拟器分析了日本南开海槽水合物藏降压开采过程中井筒完整性及海底面沉降.Jin 等[12]分析了水平井降压开采中国南海神狐海域水合物藏时储层的力学响应,结果表明开采初期沉降量占总沉降量的一半以上.李令东等[13-14]综合考虑热传导、对流、水合物分解、地层力学性质变化等过程及其相互耦合作用,建立了温度影响天然气水合物地层井壁稳定的数学模型,通过有限元求解,分析了钻井液温度对地层水合物分解、地层力学性质变化及井壁稳定的影响规律.程家望等[15]基于多孔介质流体动力学和弹性力学理论,建立了天然气水合物降压开采储层稳定性数学模型,力学分析包括储层沉降和井壁稳定性两个方面,但该模型忽略了温度对水合物分解过程的影响.孙嘉鑫等[16-17]同样搭建了TOUGH+Hydrate 和FLAC3D耦合接口程序,并对钻采条件下南海水合物井壁及储层响应特性进行模拟分析,较全面地评价了水合物钻进及降压开采过程中井壁稳定、地层沉降和开采潜能等情况,并进行了敏感性分析.

从整体上看,已有的数值软件多以搭接不同模拟器来研究水合物开采过程中的多场耦合问题,这导致了较低的计算效率.此外,大多数值模型并未全面考虑温度和水合物分解对地层力学性质的影响.因此,本文结合TOUGH+Hydrate 已有的传热-流动(TH)数学模型,基于扩展后的Biot 力学固结理论,建立能够刻画含水合物沉积层特征的传热-流动-力学(THM)耦合数学模型,开发了有限元程序对其进行数值求解.以南海神狐海域SH2 站位水合物藏条件为例,进行水合物藏降压开采的模拟计算,预测水合物开采过程中储层温度场、压力场、应力场的演化,分析了地层沉降和井壁稳定性问题.

1 传热-流动-力学耦合数学模型

传热-流动-力学耦合模型的建立存在如下假设:多相流动过程满足达西定律;相对于对流传输,溶解气体和抑制剂等在运移过程中的机械弥散作用可以忽略不计;对于应力场求解,假定岩土体饱和、骨架线弹性、变形微小;降压开采过程中不考虑地层出砂.

1.1 渗流连续性方程

含水合物地层孔隙中通常包含水、气和水合物三相,水合物分解后变成甲烷气体和水,通常只考虑气体和水为可流动相.根据连续性方程和达西定律,可得到气、水两相耦合的渗流方程为

水合物分解/形成过程中固相水合物质量守恒方程为

水合物分解/形成过程保持质量守恒,因此水合物分解速率和甲烷气体生成速率满足如下关系

式(1)～式(3)中,φ 为储层孔隙度;ρβ为相密度,kg/m3;Sβ为相饱和度;krβ为相对渗透率;µβ为相黏滞度,Pa·s;k为储层广义达西渗透率矩阵,m2;Pβ为流体压力,Pa;g为重力加速度,m/s2;mβ为单位体积储层内水合物分解产气/产水速率,kg/(m3·s);qβ为气/水源汇项,kg/(m3·s);Mw和Mg分别为水和甲烷气的摩尔质量;Nh为水分子数,甲烷水合物通常取为6.0;β=w,g,h 分别指水相、气相和水合物相.

1.2 能量守恒方程

水合物分解是一个吸热反应,因此热传递过程是影响水合物开采过程的重要因素.综合考虑热传导、对流、水合物分解和外界热量补给等因素,忽略动能和辐射传热,得到含水合物地层能量守恒方程如下

式中,ρr为岩石颗粒密度,kg/m3;Cr为比热容,J/(kg·K);T为地层温度,°C;Uβ为相内能,J/kg;Kθ为地层有效热传导系数,W/(m·K);Qθ为能量源汇项,W/m3;Qh为水合物分解热,W/m3;uβ为相流速,m/s;hβ为相热焓,J/kg.

1.3 应力场方程

1.3.1 平衡方程

忽略动量的变化量,平衡方程用有效应力可表示为

式中,σi j为地层骨架应力,Pa;Pi为地层孔隙流体压力,Pa;α 为Biot 系数,通常取值1.0;Fi为体力载荷,Pa;δi j为Kronecker 函数.

1.3.2 几何方程

几何方程定义了位移量与应变之间的关系,其张量形式为

式中,εi j为应变张量,u为位移.

1.3.3 本构方程

由于固相水合物是含水合物地层的有效胶结成分,其分解/形成会明显改变地层的力学性质.为了考虑不断演化的含水合物地层力学性质,采用弹塑性本构方程刻画应力与应变之间的内在联系,其增量形式如下

式中,dεi j为有效应力增量,Di jkl为含水合物地层的弹塑性矩阵,dεkl为应变增量矩阵.

1.4 水合物地层相关辅助方程

1.4.1 水合物分解动力学模型

在动力学平衡控制条件下,采用Kim 等[18]水合物分解动力学模型

式中,K0为水合物分解速率常数,kg/(m2·Pa·s);ΔEa为水合物反应活化能,J/mol;R为气体常数,8.314 J/(mol·K);FA为反应面积拟合参数;A为反应比表面积,m2;feq为当前温度下的平衡逃逸度,Pa;fg为当前温度下的气相逃逸度,Pa.

1.4.2 含水合物地层渗透率演化模型

水合物分解/形成过程将明显改变多孔介质的渗流通道,从而影响地层的渗透性能.本次研究采用毛管束模型刻画固相水合物演化弓起的地层渗透率变化

式中,k0为Sh=0 时的地层渗透率,m2;φ0为参考状态地层孔隙度;φc为残余孔隙度;n为地层渗透率衰减指数.

1.4.3 含水合物地层力学特性演化模型

大量三轴试验表明,含水合物地层的力学强度参数,如弹性模量、体积模量、剪切模量和内聚力等明显受到水合物饱和度的影响,但是对地层泊松比和内摩擦角的影响较小.基于Rutqvist 等[2-3]的研究结果,将含水合物地层的力学强度参数与水合物饱和度表达为如下线性关系

式(10)～式(11)中,E和C分别为含水合物沉积层的弹性模量和内聚力,Pa;Sh为水合物饱和度;E0和C0分别为Sh=0 时地层的弹性模量和内聚力,Pa;Eh和Ch分别为水合物扰动弹性模量和扰动内聚力,Pa.

1.5 地层井壁稳定性评价准则

井壁坍塌失稳大多是由于井筒周围地层发生剪切破坏弓起,通过地层应力分析结合适当的剪切破坏准则可对地层和井壁的稳定性进行预测.剪切破坏一般采用Mohr-Coulomb(MC)准则进行预测

式中,τ 为地层内某一平面上的剪应力,由应力分析求得,Pa;τc为地层临界破坏剪应力,Pa;σn为地层破坏面上的有效正应力,Pa;φ 为地层内摩擦角.

1.6 水合物开采传热-流动-力学过程耦合方法

综合以上渗流连续性方程、能量守恒方程、应力场方程、水合物地层相关辅助方程、模型初始条件和边界条件,便构成了完整的降压开采水合物藏THM 耦合分析模型.传热-流动-力学耦合过程具体分为两步:(1)根据TOUGH+Hydrate 计算的地层孔隙压力、温度和水合物饱和度作为已知条件代入到力学模型,确定含水合物地层力学性质后计算岩体骨架有效应力、应变和位移;(2)完成岩土力学计算后,根据地层有效应力更新地层孔隙度和渗透率,并反馈到下一个时间步长的传热和流动计算(图1).

图1 水合物藏开采传热-流动-力学过程耦合方法Fig.1 Methodology of coupling thermo-hydro-mechanical processes for hydrate-bearing sediments

2 模型建立

2.1 地质条件

2007 年中国地质调查局执行GMGS1 水合物钻探航次,分别在神狐海域SH2、SH3 和SH7 三个钻孔中取得水合物实物样品(图2)[19-20].本次研究以SH2 站位地质条件为基础,分析降压开采水合物藏地层井壁稳定性.SH2 站位海底面水深约1235 m,水合物赋存于海底面以下185 m,含水合物层厚约40 m.含水合物地层孔隙度在0.3～0.48 之间变化,水合物饱和度为25%～48%,渗透率为10 mD(1 mD=1.0×10-15m2)左右,地层地温梯度为(45.0～67.7)°C/km.水合物储层的上下伏地层岩性与水合物储层一致,均为弱透水性的泥质粉砂,因此研究区水合物储层为典型的2类水合物藏[21].

图2 南海神狐海域GMGS1 天然气水合物钻探航次钻孔分布图Fig.2 Location of drilling sites for the GMGS1 in Shenhu area of the South China Sea

2.2 水合物开采概念模型

根据SH2 站位资料,将其概化为水平延伸地层,建立如图3 所示的平面应变模型,模型尺寸为200 m×255 m (x,z).垂直方向:水合物储层厚40 m,水合物储层之上为185 m 厚的透水顶板,其下为30 m厚的下伏地层,模型顶面即为海底面;水平方向模型侧向延伸200 m.生产井位于x=0 m 位置处,半径r=0.1 m.模型顶面为海底面,设定为定温定压边界;由于钻探结果表明水合物储层下伏地层未存在隔水泥岩层,因此模型底面同样设定为透水边界;模型侧边界设定为零流量边界.对于力学边界,将模型的侧面沿法线方向上的位移设为零,即固定位移边界;模型底面z=255 m 的截面,保持z方向上的位移为零;模型最顶面z=0 m 的截面为海底面,模型中设定为自由位移边界,即降压开采过程中海底面可以发生自由移动;井壁在模型中设定为刚形体,井壁和地层之间为刚性接触,即井壁在模型计算过程中只能承受应力,但不发生变形.

图3 神狐海域SH2 站位水合物降压开采传热-流动-力学耦合概念模型Fig.3 Conceptual model for coupling thermal-hydro-mechanical processes of gas production by depressurization from Shenhu SH2 hydrate reservoir

已有的数值模拟研究表明水合物开采过程中的关键过程主要发生在生产井附近10～20 m 范围内.为了准确刻画水合物的分解演化过程,获取近井地层水合物分解的重要现象,模型对生产井附近的数值网格进行加密处理(见图3).水平方向剖分20 个0.5 m,10个1.0 m,15 个2.0 m 和30 个5.0 m 子网格,共计剖分75 列;垂直方向水合物储层内网格精度均为1.0 m,边界地层网格精度为2.0 m,共计剖分148 层.因此,整个储层模型被离散为75×148=11 100 个计算网格.

2.3 模型参数

模型所采用的地层物性参数主要依据现场测试数据,详见表1.为便于研究,水合物储层含水合物饱和度取值0.4,地层孔隙度依据现场岩性分析数据取值0.38,渗透率为10 mD.水合物储层底界初始压力和初始温度分别为14.7 MPa 和14.1°C,保证了水合物储层在开采前的稳定存在状态.含水合物地层的相对渗透率和毛细压力计算分别采用Stone 模型和van Genuchten 模型,这些模型被广泛用于从水合物藏回收甲烷气体的数值模型中[2,5-6,15].

由于缺乏现场水合物储层力学测试数据,本次研究含水合物地层力学参数依据Masui 等[22-23]的室内三轴试验结果选取(表2).水合物藏开采过程中地层的力学强度参数随水合物饱和度的变化而变化[24-29].

表1 模型计算所用地层物性参数[5-6,19-20]Table 1 Parameters of hydrate deposits used for simulation [5-6,19-20]

表2 含水合物地层力学强度参数Table 2 Mechanical parameters for hydrate-bearing sediments

2.4 初始条件

图4 显示了模型初始条件,包括地层温度、孔隙压力、垂向有效应力和水合物分布.模型垂向上的初始温度场根据海底面的初始温度依地温梯度变化计算得到;由于研究区地层为未固结沉积层,初始压力场服从静水压力分布;模型初始垂向有效应力假设仅仅是由于上覆地层自重弓起,初始有效水平应力依据侧压力系数K0计算得到;水合物储层位于海底面以下185 m 至225 m 之间,水合物初始饱和度为0.4.

2.5 生产方案

本次研究利用竖直井进行降压开采条件下的产气能力与沉降规律预测.由于水合物储层上下为透水地层,将竖直井射孔段布设于水合物储层中部,长度为10 m(图3).生产井射孔段距水合物储层顶底面各15 m,这样有利于水合物储层的有效降压,限制边界地层水快速涌入生产井筒.模型连续开采模拟时间为1 年,井筒降压幅度设定为10 MPa.

3 结果分析

3.1 气水产出动态特征

图5 显示了垂直井降压开采2 类水合物藏水合物分解甲烷释放速率(QR),生产井产气速率(QP),水合物分解(VR)和生产井(VP)累积产气体积随时间演化特征.降压开采初期(前20 d),QR和QP均快速降低,这主要是由于开采初期生产井和周围储层存在较大的压力梯度,导致生产井周围大量的水合物发生分解.随着开采的持续进行,低压区域逐渐向储层内部发展,生产井和水合物储层之间的压力梯度逐渐趋于平缓,从而导致开采后期QR和QP逐渐降低.连续开采1 年后,VR=1.36×104m3,VP=1.20×104m3,甲烷气体回收效率可达88%,由此可见降压技术是开采2 类水合物藏的有效方法.

图4 模型初始温度场、孔隙压力场、垂向有效应力场和水合物储层分布特征Fig.4 Initial conditions of temperature,pore pressure,effective stress,and hydrate saturation

图5 水合物分解甲烷释放速率QR、生产井产气速率QP、水合物分解累积产气体积VR和生产井累积产气体积VP随时间演化特征Fig.5 Evolution of gas release rate(QR),production rate(QP),and the cumulative volumes of released from hydrate(VR)and produced(VP)in the production well

图6 显示了生产井产水速率(QW)和气水体积比(RGW)随时间演化特征.开采初期QW快速降低,与QP和QR相似,该阶段产水主要来自于水合物储层.连续开采大约20 d 后,QW随时间逐渐增大.这主要是由于开采后期水合物储层下部的水合物逐渐分解(图7),有效渗透率增大,导致生产井与下伏地层之间的水力连通性增强,大量自由水从下伏地层涌入生产井造成QW不断增大,这与苏正等[5,30]的研究结果一致.

图6 生产井产水速率QW和气水体积比RGW随时间演化特征Fig.6 Evolution of water production rate(QW)and gas to water ratio(RGW)

图7 显示了模型预测在降压开采1 年后,地层中孔隙水压力、温度、水合物饱和度和甲烷气饱和度的空间分布特征.从图中可以看出:(1)井筒降压导致井周地层孔隙水压力降低,并迅速在储层中传播;(2)储层降压打破了水合物的相平衡状态,根据低压区分布可知射孔段周围水合物优先发生分解,且连续开采1 年后水合物分解前缘已扩展至下伏边界地层;(3)水合物分解吸热导致生产井射孔段周围形成一低温区域,由于开采后期生产井与下伏地层取得水力连通导致下伏边界层中的相对高温地层水涌入生产井,并在生产井射孔段以下形成一明显的高温热柱,这指示出下伏边界层中地层水涌入生产井的流体运移通道;(4)从甲烷气饱和度空间分布图可以看出,水合物分解区主要集中在生产井周围25 m以内.

图7 降压开采365 天后各参数的空间分布特征(续)Fig.7 Spatial distribution of the parameters after 365 days depressurization(continued)

综合分析图5～图7 表明,开采后期生产井与下伏饱水地层取得水力连通后水合物开采效率随时间快速降低.这主要是由于(1)大量自由水从下伏地层涌入生产井导致水合物储层降压受到抑制,水合物分解前缘难以向储层内部扩展;(2)大量回收地层水将过度消耗外部能量.因此,今后针对2 类水合物藏开采潜力预测研究中应对布井和开采方案进行优化,以保证较高的水合物开采效率.

3.2 地层沉降特征

水合物藏降压开采过程中地层有效应力的演化是预测地层破坏和出砂位置的关键.图8 中显示了降压开采1 年后地层中(a)最小水平有效主应力和(b)最大有效主应力的空间分布特征.从图中可以看出,降压导致生产井射孔段周围地层出现明显应力集中,因此,这些区域是井壁破坏和地层出砂的主要位置.由于下伏透水地层的补给作用,导致降压难以在水合物储层中传播,因此有效应力增大区域主要限定在生产井周围地层中.

图8 降压开采365 天后地层有效主应力空间分布特征Fig.8 Spatial distribution of effective stress after 365 days depressurization

为了进一步定量分析井周地层是否发生屈服破坏,在应力集中区域选择2 个监测点进行有效主应力路径和破坏分析(图9),监测点坐标(x,z)分别为(2 m,204.5 m)和(10 m,204.5 m).从图中可以看出,在1 年的开采周期内,生产井周围地层不会发生屈服破坏.这主要是由于降压弓起垂向有效应力和最小水平有效应力同时增大,并未出现明显的偏应力.另一方面,可渗透下伏地层的补给作用限制了井周地层有效应力的增大,这在一定程度上缓解了地层出现屈服破坏.

地层孔隙压力降低造成骨架有效应力增大,从而弓起地层变形沉降.图10 显示了降压开采过程中水合物储层顶部发生沉降的空间分布特征.从图中可以看出,井壁附近地层的沉降量最大,向储层内部延伸沉降量逐渐减小.这主要是由于(1)生产井周围地层降压幅度较大,且出现明显的应力集中现象(图8);(2)生产井周围水合物发生完全分解(图7),降低了地层的力学强度,从而加剧了井周地层的沉降.随着降压开采的持续进行,地层沉降逐渐增加.连续开采50 d 时,最大沉降量为0.13 m.然而,连续开采1 年后,最大沉降量增大到0.19 m,相对50 d 时的沉降量仅增大了0.06 m.由此可以看出垂直井降压开采水合物,地层沉降主要发生在开采初期.因此,开采前期可以采用逐步降压的方式缓解快速的地层沉降阶段,这样可以有效避免快速沉降弓起井筒套管破裂的风险.

图9 水合物储层中不同监测点位置有效主应力路径Fig.9 Effective principal stress paths at different monitoring points in hydrate reservoir

图10 水合物储层顶部空间沉降特征Fig.10 Spatial distribution of the vertical displacement at the top of hydrate reservoir

3.3 井壁稳定性分析

图11 沿生产井方向井壁孔隙压力、温度、水平有效主应力和垂向有效主应力分布Fig.11 The distribution of pore pressure,temperature,horizontally effective stress,and vertically effective stress along the production well

井壁坍塌失稳是制约水合物藏钻井和开采过程顺利实施的关键问题,一般可通过井周应力状态分析结合适当的剪切破坏准则进行预测.图11 显示了沿生产井壁孔隙压力、温度、水平有效主应力和最大有效主应力在不同时间的分布特征.从图中可以看出,井筒降压导致射孔段附近压力快速降低,降压幅度受井筒内流体压力控制,随后低压逐渐沿井壁向上下扩展.水合物分解吸热效应和低温流体流向生产井的共同作用导致射孔段井壁温度明显降低,最大降温幅度可达6°C,但是开采后期由于下伏地层的相对高温流体涌入导致射孔段井壁温度升高.在孔隙压力和温度的共同作用下,射孔段井壁有效应力出现明显集中.连续降压开采1 年后,射孔段井壁水平有效主应力最大为9.1 MPa,相比于初始有效应力增加了6.6 MPa;垂直有效主应力最大为8.5 MPa,相比于初始有效应力增大了5.2 MPa.有效应力在水平方向的增加量大于垂直方向,这主要是由于水平方向受到限制,而垂直方向地层可以发生自由移动.

井壁有效应力增大可能诱发井壁发生坍塌破坏,按照之前提到的破坏分析方法,计算得到开采1 年后的井壁破坏风险分析值,其沿井筒分布特征如图12 所示.从图中可以看出,随着开采的持续进行井壁发生剪切破坏的风险逐渐增加,且在射孔层段井壁破坏风险最大[31].由于射孔段井壁应力集中最为明显,导致射孔段井壁最易发生剪切破坏(图12).这表明射孔段是储层出砂和坍塌破坏的高危区域,一旦射孔段发生坍塌破坏将导致大量的地层骨架介质涌入生产井筒,造成井筒失效.为了有效避免生产井壁坍塌,射孔段井筒套管应具备足够高的力学强度以支撑降压弓起的应力集中,同时防止井周地层发生破坏.

图12 沿生产井方向井壁发生剪切破坏风险Fig.12 The distribution of possibility of shear failure along the production well