海域天然气水合物降压开采诱发储层力学性质劣化及沉降规律建模研究

摘要

降压法是海域天然气水合物储层开采的一种常见方法。降压开采会导致近井储层出现复杂的多物理场耦合响应，诱发压力变化、温度变化、水合物分解、储层力学性质劣化及地层沉降。本研究通过一种全耦合流固热化数值模型分析水平井筒降压导致的海域天然气水合物储层力学性质劣化及沉降特征，表征水平井筒及井周储层的多场耦合响应规律，明确储层力学性质劣化区域及沉降的影响因素。模拟结果显示：储层压力和温度变化的波及区域远大于水合物分解前缘，有效正应力的分布在不同方向差异明显，降压开采诱发的内聚力劣化区域与塑性区和水合物分解区关联程度高，水平井筒以浅区域和以深区域的沉降呈现出不同特征，沉降程度可相差5 mm以上。模拟结果对水平井降压开采海域天然气水合物的储层稳定性分析具有参考意义。

关键词

天然气水合物; 降压开采; 储层力学性质; 数值模拟; 弹塑性变形; 地层沉降; 水平井

0 引言

天然气水合物资源量大，常存在于陆域冻土和海底沉积层内^［1］。海域天然气水合物主要通过降压法、热激法和注抑制剂法进行开采^［2-5］。天然气水合物试采现场常通过降压法降低天然气水合物储层压力，引起储层内的天然气水合物相变^［6-7］。然而，降压法开采诱发天然气水合物分解的同时也会导致储层力学性质发生变化。由于海域水合物储层多为弱胶结沉积物，导致储层骨架出现变形、破坏、地层沉降等现象，影响储层的稳定性，严重时可能进一步诱发地质灾害，严重影响天然气水合物降压法开采的安全性和效率，需要开展针对性研究^［8-10］。

多场耦合理论被广泛运用于海域天然气水合物钻采诱发的储层岩土力学与岩石力学响应研究中。Rutqvist和Moridis^［11］提出了一种联合TOUGH+HYDRATE和FLAC3D软件的模拟方法并结合Mohr-Coulomb判据表征了水合物开采诱发的储层稳定性问题。这种方法能够分别发挥TOUGH-HYDRATE在水合物分解、气水流动、热传递方面的优势和FLAC3D在计算固体骨架变形方面的优势。该方法未将岩石固体骨架变形的求解与水合物相变分解放在同一数值系统求解。在此基础上，J. Kim等^［12］提出了一种双向流固耦合算法，提高了孔隙压力和骨架变形的求解精度。Sanchez等^［13］提出了一种全耦合方法，联合室内物理实验和天然气水合物开采现场尺度模拟，表征了降压法开采天然气水合物过程中的流固热化耦合规律。通过流固热化多场耦合模拟，能够预测降压过程中的孔隙压力和温度传递时空演化规律，表征水合物分解过程中储层岩石骨架的应力释放和体积变化。改变降压方式和工程参数会对水合物储层的变形产生明显的影响^［14］。在天然气水合物储层流固热化耦合模型的控制方程推导和模型构建的过程中，储层骨架力学性质对模型稳定性和准确性的影响较大，内聚力和刚度等典型的强度和变形参数的取值需要进行准确的判断。在多场耦合条件下，这些力学参数会影响储层骨架的应力应变特征、剪切破坏规律和水合物分解及二次生产^［15］。当天然气水合物储层的压力和有效应力发生显著变化后，塑性变形和塑性破坏等现象出现，并可能诱发储层出砂，影响天然气水合物储层的流动特性和井筒的稳定性^［16-17］。水合物储层出砂过程可以通过计算流体力学方法和离散元法建立数值模型表征，水合物颗粒和砂颗粒的运移过程中，颗粒受力情况复杂，受到流动方向的影响十分显著^［18］。水平井降压开采诱发的储层流固热化响应与直井开采具有差异，水平井的布井方式和布井数量均会对多场耦合响应造成影响，导致井周储层内的岩石骨架变形破坏规律出现不同^［19］。李阳等^［20］提出了一种适用于南海神狐海域水合物储层的多场耦合井壁稳定性分析方法，为钻井过程中的钻井液密度、钻井液温度和钻井时间的优化提供了理论参考。孙金等^［21］利用Duncan-Chang模型描述水合物储层的变形破坏特征，并结合渗流场、温度场、变形场进行了耦合建模，发现降压开采的井底压力、储层厚度、弹性模量、储层渗透率等因素均控制降压开采诱发海床沉降的程度。由于水合物开采过程中储层沉降可能导致井筒出现位移和变形，井筒-土壤耦合模型能够同时表征水合物储层变形和井筒的应力应变规律，为安全试采提供工程参数的参考^［22］。天然气水合物储层的流固热化建模过程中，空间建模常采用二维或三维模型。二维模型能够表征特定平面内的储层响应，而且数值模型求解效率高。三维建模则能够表征近井地带多个方向的压力、应力应变和温度的变化，能够比二维模型考虑更多的渗流流型，但也会显著提高算力需求，导致数值模型的求解速度降低^［23］。

已有研究说明，多场耦合理论能够为天然气水合物储层降压法开采诱发的储层压力、变形、温度、破坏等响应提供数学建模依据，为海域天然气水合物安全、高效钻采提供工程参数优化的参考。本文提出了一种海域天然气水合物储层全耦合流固热化数值模型，模型考虑了水合物分解、热传导和储层骨架弹塑性变形与破坏本构。基于全耦合模型，提高了各物理场数值求解时的稳定性和准确性，分析典型天然气水合物储层和降压参数条件下的近井储层力学响应规律，量化相关因素对海域天然气水合物储层降压开采诱发的多场耦合响应及储层沉降的影响特征。

1 控制方程

本研究通过连续介质表征天然气水合物储层，考虑水合物分解过程储层内的质量守恒、能量守恒、静力平衡和水合物分解动力学^［19，23-24］。根据质量守恒和渗流理论，储层中气和水的渗流表示为：

$$\frac{\partial }{\partial t}\left(\phi S_{\mathrm{w}}{\rho }_{\mathrm{w}}\right)+\nabla ·({\rho }_{\mathrm{w}}{v}_{\mathrm{w}})=s_{\mathrm{w}}$$

（1）

$$\frac{\partial }{\partial t}\left(\phi S_{\mathrm{g}}{\rho }_{\mathrm{g}}\right)+\nabla ·({\rho }_{\mathrm{g}}{v}_{\mathrm{g}})=s_{\mathrm{g}}$$

（2）

式中：$\phi$——储层孔隙度；$S_{\mathrm{w}}$、$S_{\mathrm{g}}$——各相的饱和度；${\rho }_{\mathrm{w}}$、${\rho }_{\mathrm{g}}$——各相的密度；$v_{\mathrm{w}}$、$v_{\mathrm{g}}$——各相的速度；$s_{\mathrm{w}}$、$s_{\mathrm{g}}$——各相的流入或流出质量速度。

降压法开采过程中，水合物分解速度通过动力学模型表示：

$$R_{\mathrm{M}\mathrm{H}}=-k_{\mathrm{d}}{M}_{\mathrm{M}\mathrm{H}}\left(p_{\mathrm{e}}-p_{\mathrm{g}}\right)A_{\mathrm{s}}$$

（3）

式中：$R_{\mathrm{M}\mathrm{H}}$——水合物分解速率；$k_{\mathrm{d}}$——反应速率；$M_{\mathrm{M}\mathrm{H}}$——摩尔质量；$p_{\mathrm{e}}$——相平衡压力；$p_{\mathrm{g}}$——气相压力；A_s——反应比表面积。

在多场耦合过程中，储层内温度和压力的变化会导致储层固体骨架和孔隙流体发生体积变化，该过程表示为：

$$x_i=\frac{\mathrm{1}}{{\rho }_i}{\left(\frac{\partial {\rho }_i}{\partial p_i}\right)}_T$$

（4）

$${\beta }_i=-\frac{\mathrm{1}}{{\rho }_i}{\left(\frac{\partial {\rho }_i}{\partial T_i}\right)}_p$$

（5）

式中：x——与压力变化相关的压缩系数；$\beta$——与温度变化相关的膨胀系数；$p$——压力；$T$——温度。

在表征储层变形过程中，通过应力张量和静力学平衡的控制方程进行计算：

$$\nabla ·\sigma =\mathrm{0}$$

（6）

式中：$\sigma$——应力张量。

方程未考虑重力对应力变化的作用。

由于本研究关注近井地带的沉积物损伤与力学性质劣化，对于储层固体骨架采用弹塑性本构：

$$\delta \sigma =C:\delta \left(\mathit{\epsilon }-{\mathit{\epsilon }}_{\mathrm{p}}\right)-\alpha \delta pI$$

（7）

式中：$C$——弹性张量；$\epsilon$——应变；${\mathit{\epsilon }}_{\mathrm{p}}$——塑性应变；$\alpha$——有效应力系数；$I$——二阶单位张量。

近井地带的储层固体骨架破坏判据通过Mohr-Coulomb模型和Drucker-Prager模型表示，判断降压法开采不同阶段、不同位置的弹、塑性状态及塑性破坏程度。

通过有限元方法对上述方程进行空间离散，通过后向欧拉法进行时间离散，构建数值系统进行求解。求解时采用全耦合法，使用稀疏矩阵直接求解器进行矩阵求解。

2 数值模型建立

基于上述控制方程，建立海域天然气水合物储层降压法开采二维模型，针对近井短水平井降压开采的问题进行模拟和分析。图1所示为近井地带40 m长度和40 m厚度区域，水平井筒长度2 m位于模型中心。渗流场中，顶部、底部、左侧、右侧边界均为封闭边界，无流体流入和流出；温度场中，4个边界均为隔热边界，无热量流入和流出；应力场中，左侧、右侧、底部边界均为辊支撑边界，顶部为自由边界以表征沉降作用。由于全耦合模型对数值稳定性要求高，如图1所示采用加密网格处理。储层固体骨架杨氏模量为40 GPa，泊松比为0.2，骨架密度为2600 kg/m³，渗透率为2 mD，孔隙度为0.15，初始天然气水合物饱和度为40%，初始孔隙压力为14 MPa，井底生产压力为3 MPa。储层强度参数中，初始内聚力为0.861 MPa，内摩擦角为12°，降压过程中，内聚力变化可以表示为：

图1  二维模型及边界条件示意

Fig.1  Sketch of the 2D model and the boundary conditions

$$c=c_{\mathrm{0}}+\frac{\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}\phi }\alpha S_{\mathrm{M}\mathrm{H}}^{\beta }$$

（8）

式中：$c_{\mathrm{0}}$——初始内聚力；$\alpha$、$\beta$——系数。

3 结果与讨论

根据建模参数，模拟降压法20 d内诱发的储层压力、温度、有效应力、水合物饱和度、沉降等参数变化规律，为海域天然气水合物降压开采过程中井周地层的多场耦合响应规律提供理论参考。

3.1 压力、温度和水合物分解规律

图2所示为降压开采0.5、5、10、20 d时储层内的压力分布情况。结果显示水平井筒降压能够较快波及到近井地带，降压开采初期（例如0.5和5 d），由于水平井筒几何形状的影响，渗流形态为椭圆形流动，导致压降区域的扩展不是严格的径向流。当降压开采时间进一步增加，流动形态逐渐偏向于径向流，压降均匀向储层边界波及。这说明采用水平井筒进行天然气水合物降压开采时，储层压降规律不同于点状压降诱发的径向流，而是前期为椭圆流动，后期为近似径向流的状态。5 d的压力分布结果显示，水平井筒降压开采能够较快降压前缘波及至近井地带，至20 d时储层压力已经显著下降。这说明水平井筒降压能够建立近井地带压差，促进水合物分解和向井内流动。在模拟时间内，压降能够有效波及至井周20 m范围内。

图2  不同降压时间的储层压力分布

Fig.2  Pressure distribution in the reservoir at different depressurization time steps

图3所示为温度随降压开采时间变化。与压力时空演化规律类似，温度变化在初期也呈现近似椭圆形扩展，后期近似径向扩展。这是由于气水渗流和温度在储层的传递耦合程度更高导致的。对比图2的压力分布和图3的温度分布发现，压降波及程度远大于降温范围，这是由于降压诱发水合物分解，水合物分解吸热、储层降温、温度梯度导致热传递等过程随后发生，因此储层压降比温度降低波及更远。

图3  不同降压时间的温度分布

Fig.3  Temperature in the reservoir at different depressurization time steps

图4所示为天然气水合物饱和度随着降压开采的进行的变化规律，水合物分解区域随着降压时间的增加而增加，且分解区形状近似椭圆形，椭圆形长轴方向与水平井筒走向一致。这一结果与图2的压降规律存在关系：压降导致相平衡破坏，水合物开始分解，由于降压初期孔隙压力降低近似椭圆形流动，在长轴方向更加明显，导致水合物分解在椭圆长轴（水平井筒方向）更加明显。随着降压开采时间增加，水合物分解区逐渐扩大，降压20 d后，水合物分解前缘最远扩展至井筒外3 m左右。通过分析色标梯度发现，降压0.5、5和10 d后，水合物分解前缘明显，由完全分解至未分解过渡很快；降压20 d后，水合物分解前缘相对平缓，水合物饱和度由0过渡至40%的距离略微变长。

图4  不同降压时间的天然气水合物饱和度分布

Fig.4  Natural gas hydrate saturation in the reservoir at different depressurization time steps

通过图2~4所示不同时刻的储层压力、温度、水合物饱和度分布发现，在降压开采后的特定时间内，压力场、温度场和水合物分解场呈现出相关性，水平井筒降压开采导致的多场耦合响应应当区别于其他井型开采诱发的储层响应。在温度、压力、水合物饱和度变化规律中，压降波及最远，温度降低波及的程度次之，水合物分解前缘推进的速度最慢。这一现象是由于各物理场不同的控制机理导致的。

3.2 井周储层力学响应规律

降压开采诱发近井地带储层的多场耦合响应，在已有温度、压力、水合物分解等特征的基础上，进一步分析近井储层力学响应规律。

图5所示为降压开采0.5 d和20 d的储层${\mathrm{\sigma }}_{\mathrm{x}}\text{'}$和${\sigma }_{\mathrm{z}}\text{'}$分布情况，分别表示沿水平井筒方向和垂向的有效正应力。此应力由降压开采诱发，能直接反映各方向储层骨架受到的压缩状态应力分布。结果显示，由于水平井筒直接施加压降，水平井筒区域有效应力更高，且最大值随着压降时间的增加而上升。${\sigma }_{\mathrm{z}}\text{'}$在水平井筒两端更易出现应力集中现象。结果显示，在近井地带较小范围内建模的条件下，模型边界均会由于压力变化而出现有效应力上升，且有效应力上升的特征在水平方向和垂直方向均有不同。整体上，由于海域天然气水合物常位于深水浅层，上覆地层应力相对于水平方向地应力较高，也会导致${\sigma }_{\mathrm{x}}\text{'}$和${\sigma }_{\mathrm{z}}\text{'}$分布出现差异。

图5  不同降压时间的水平有效正应力和垂向有效正应力分布

Fig.5  Effective stresses in the horizontal and the vertical directions in the reservoir at different depressurization time steps

图6所示为近井地带储层破坏规律云图，通过塑性应变（体积应变）表示。由于水平井筒两端存在应力集中现象，井筒两端出现较明显的塑性应变。由于降压开采引入较大的生产压差，0.5 d就诱发了明显的井周塑性变化，塑性区主要出现在井筒和井周10 m范围内，模型边界塑性变形不明显。随着降压开采时间的增加（20 d），井周塑性区及塑性变形进一步增加，但波及范围并未扩展至模型边界，塑性破坏仍以近井区域为主。

图6  不同降压时间的近井地带塑性破坏区分布

Fig.6  Plastic regions in the near‑well area at different depressurization time steps

图7所示为降压开采20 d后近井地带内聚力和沉降分布情况。由于近井地带水合物分解及塑性演化，储层内聚力出现明显降低，强度发生明显劣化，内聚力劣化的区域与水合物分解区域高度重合，体现了天然气水合物对储层力学强度的主控作用。水平井筒对应区域的上覆地层出现垂向沉降，水平井筒正下方的下伏地层沉降程度略低于水平井筒以外区域，这一现象受到降压诱发的多场耦合作用控制，在水平井筒上方储层骨架压缩明显促进沉降，在水平井筒下方降压诱发的储层骨架压缩减缓地层沉降的程度。

图7  降压开采20 d后的近井区域内聚力分布及沉降特征

Fig.7  Plastic regions in the near‑well area at different depressurization time steps

图8所示为水平井筒中心（0 m，0 m）及其以浅2 m处（0 m，2 m）的内聚力劣化及储层沉降程度随降压开采时间的关系。在该结果中，展示沉降程度的绝对值以表征降压开采对地层沉降的影响大小，正值代表地层沉降的距离。在水平井筒及其以浅区域，地层沉降均随着时间单调递增，且水平井筒以浅的观测点沉降值更高。内聚力劣化结果显示，井筒中心由于在0时刻就施加压降载荷，导致其初期就出现明显的水合物分解、破坏及内聚力劣化，随后，由于水合物完全分解，并已进0入塑性区，内聚力变化不明显。在井筒中心以前2 m处，内聚力在降压开采5 d后开始出现下降，表明水合物分解和塑性区演化影响到该处，随后内聚力在2 d内持续下降，直到水合物完全分解，其数值趋于稳定。

图8  水平井筒中心及其以浅2 m处的内聚力劣化及沉降演化规律

Fig.8  The evolution patterns of cohesion and subsidence at the center of the horizontal wellbore and at 2 m above it

图9所示为降压开采过程中的压力和水合物饱和度分布情况。结果显示在20 d内，压降有效波及至模型边界，建立近井沉积物层内压差；水合物分解前缘扩展相对缓慢，20 d内扩展至短井筒以外5 m。

图9  不同时刻近井区域压力和水合物饱和度分布

Fig.9  The distribution of pressure and hydrate saturation near the wellbore at different time steps

4 结论

本研究通过海域天然气水合物储层流固热化耦合建模，形成了数值稳定性较好的全耦合模型，分析了水平井筒降压开采诱发的近井地带压力、温度、变形、破坏区变化规律，考虑了内聚力劣化及塑性破坏对全耦合过程的影响。得到以下结论：

（1）降压法开采过程中，近井地带流动形态受到水平井筒的几何尺寸影响，前期呈现近似椭圆形，后期呈现径向特征；压降波及程度较温度降更广，而水合物分解区扩展速度最慢。

（2）降压诱发的有效正应力在水平井筒方向和垂向呈现不同的响应特征，水平井筒两端的位置出现应力集中现象。这种响应特征的差异化也与深水浅层的上覆地层应力和水平地应力场有关。

（3）降压法诱发近井地带的水合物分解和塑性区域演化对地层的强度劣化及沉降均有直接影响。内聚力的劣化与水合物分解和塑性区演化关联程度高。在水平井筒以浅处，地层沉降效果被水平井筒降压增强；在水平井筒以深处，地层沉降一定程度上被水平井筒降压所缓解，沉降程度可以降低约5 mm。