天然气水合物藏井眼蠕变缩径研究

李 阳， 程远方, 闫传梁， 周晓辉

(中国石油大学(华东)石油工程学院， 青岛 266580)

天然气水合物具有资源量大、能量密度高和低污染等优点[1]，被公认为是具有良好前景的重要后续清洁能源[2]，全球天然气水合物的碳储量相当于现已探明的石油、煤和天然气等常规化石能源总含碳量的2倍[3-5]。受水合物相平衡条件和地温梯度的限制，天然气水合物主要分布在深海和永久冻土区，埋存深度浅，沉积物胶结程度差[6-8]。各国学者对于水合物的力学性质展开了大量研究[9-11]，发现水合物储层具有极强的蠕变特性[12]。由于水合物在储层中起到了很强的胶结作用，因此水合物饱和度的不同会对其蠕变特性产生非常大的影响。对于不同的水合物藏，其水合物饱和度必然不同，同时钻井和开采等工程作业也会导致水合物分解，影响其饱和度和蠕变特性[13]。

由于钻井工程的扰动，井眼附近会产生应力集中，对于蠕变性地层，应力的增大会使其产生蠕变变形。蠕变对钻井最大的影响是会造成井眼缩径，影响起下钻作业。目前对于蠕变地层井壁稳定的研究主要在盐岩地层和软泥岩地层，近年也开始关注脆性岩石蠕变。殷有泉等[14]采用理论计算和数值模拟的方法，分析了蠕变地层中套管载荷，认为套管载荷的黏弹性解最终趋近于弹性解。Willson等[15]、Poiate等[16]建立了二维有限元或有限差分模型，分析了水平地应力作用下盐膏层中套管载荷。Machay等[17]、Wang等[18]等建立了三维有限元或有限差分模型，分析了三维地应力作用下盐膏层中套管力学状态。Machay等[19]又提出了分步模拟方法，对钻井和固井等过程井眼的力学状态进行了模拟分析。Lao等[20]建立了盐膏层中斜井井筒的有限元模型，分析了井斜角和方位角对井眼缩径和套管力学状态的影响规律。Li等[21]建立了热流固三场耦合的冻土地层井眼收缩力学模型，并分析了钻井液温度、地层温度及井眼打开时间对井眼状态的影响。

这些不同地层井眼缩径的研究对于水合物地层钻井作业过程中井眼的变形分析有很大的借鉴意义，但水合物地层的蠕变特性及其影响因素都与其他地层有很大的不同，目前关于水合物储层蠕变特性对井眼缩径规律的影响还鲜有研究。因此，现从应力的角度对钻穿水合物储层时井眼的蠕变缩径规律进行研究，并分析水合物储层井眼蠕变缩径的影响因素及其作用规律。

1 数学模型

1.1 流固耦合模型

地层的应力场与渗流场是相互作用、互相影响的。应力场作用使得岩石骨架发生变形，改变地层的渗流特征，同时渗流场的变化又会改变地层的孔隙压力分布从而影响岩石骨架承受的有效应力，孔隙压力对地层有效应力影响的计算式为

σ′=σ-αpp

(1)

式(1)中：σ′ 为有效应力；σ为地层中任一点的总应力；α为有效应力系数；pp为孔隙压力。

由于实际的多孔介质微观组成复杂，很难用精确的数学关系将流体固体分开，因此所建立的模型为理想的连续系统，忽略岩石孔隙的微观结构变化，但通过一些特殊变量的变化可以表征孔隙流体流动的过程，以此模拟流固耦合的过程。

对于固相材料的应力变化采用虚功原理表征[22]，即

(2)

式(2)中：V为计算区域体积；Sa为牵引力作用下的表面积；σ为总应力矩阵；δε为虚应变率矩阵；ts为表面牵引力矢量；f为体力矢量；δv为虚速度矢量；：为双点积。

通过拉格朗日法将固相的虚功方程离散到网格之中，是流体流经这些网格，流体满足连续性方程，如式(3)所示，流体的渗流过程遵循Darcy定律。

(3)

式(3)中：pf为孔隙流体密度；φ为多孔介质孔隙度；vfp为渗流速率矩阵；n为外表面外法线方向向量。

1.2 水合物储层蠕变模型

蠕变分析是用于表征模型受力发生弹性变形后基于其黏塑性属性而发生的一种后效反应。单轴蠕变实验中，当材料处于弹性阶段进行蠕变实验时，其总应变可分解为弹性应变和蠕变应变之和，即

ε(t)=εe+εc

(4)

式(4)中：ε(t)为总应变；εe为弹性应变；εc为蠕变应变，是温度T、时间t以及应力σ的函数，即

εc=f(σ,t,T)=f1(σ)f2(t)f3(T)

(5)

采用时间硬化法则时，σ为有效应力，f1(σ)采用Norton幂律法则，则其初始蠕变和第二蠕变阶段可表达[23]为

(6)

式(6)中：A、m、n为蠕变参数，由蠕变实验获得。

将式(5)对时间微分并化简即可获得时间硬化法则，即

(7)

1.3 水合物储层蠕变模型

实验研究表明，低温水合物具有较强的蠕变特性，其蠕变曲线与富冰冻土在温度接近0 ℃时蠕变曲线接近[24-25]，其蠕变本构模型为

(8)

式(8)中：ε为应变；δ为轴向差应力；θ为负温；t为时间；A、B、C、D为实验拟合参数。拟合结果如图1所示。

图1 冻土水合物层蠕变模型拟合结果[24]Fig.1 Creep curves of the sediments with gas hydrates[24]

2 数学模型

采用ABAQUS有限元软件对所建立的数学模型进行计算。ABAQUS软件是一个功能强大的非线性有限元软件，广泛应用于岩土力学的相关计算[26-27]。同时ABAQUS拥有丰富材料数据库可满足大多数工程需要，UMAT子程序自主定义任意材料属性。

2.1 水合物储层井眼模型

水合物储层开采流固耦合力学模型如图2所示，储层半径为2.5 m，AB边为井眼，井眼半径为0.25 m，分别在AB边、CD边和DE边设置孔压边界，CD、DE同时设置位移边界；BC边和AE边为轴对称边界。模型单元类型为CPE4P，即4节点平面应变双线性位移、双线性孔压单元，模型总网格数为3 200。

图2 流固耦合力学模型示意图Fig.2 Schematic diagram of fluid-solid coupling mechanics model

2.2 水合物储层井眼模型

假设该水合物层所在深度为500 m。模拟所需数据如表1[28-30]所示。

表1 地层基本参数[28-30]Table 1 Basic parameters of the model[28-30]

2.3 地应力平衡

地应力平衡是用于模拟地层在地应力作用下的状态，在上覆岩层压力、最大水平地应力、最小水平地应力作用下地层处于应力平衡状态，当使用ABAQUS模拟地应力平衡不会出现明显的位移和应变时，说明所确定的地应力参数是合理的。地应力平衡通常是使用ABAQUS模拟岩土问题时的第一步。图3显示了模型模拟的结果，土体内的最大位移为10-18m的量级，这表明在Geostatic步骤中给出的所有初始设置在随后的分析中是有效的。

图3 地应力平衡Fig.3 Initial in-situ stress equilibrium

3 水合物层井眼缩径规律

3.1 井眼钻开时初始变形

使用钻头钻开地层的实质是地应力的释放过程，在ABAQUS模拟岩土开挖过程最常用的方法是“生死单元法”，通过赋予土体单元“死”的特性来表示外部作用使这一部分土体离开地层。同时对井眼处施加钻井液液柱压力，在地应力、钻井液液柱压力、孔隙压力作用下，地层会出现一个初始的弹性变形，其应力分布如图4所示，位移分布如图5所示。同时由于钻井液液柱压力大于地层压力，钻井液会向地层发生渗流，使井周孔隙压力发生改变。

图4 井眼钻开后初始米塞斯应力分布Fig.4 Initial Mises stress distribution after drilling

图5 井眼钻开后地层初始位移Fig.5 Initial displacement of the formation after drilling

由地层初始应力分布和位移分布可以看出，地层的井口会出现一个明显的应力集中，且在最小水平地应力方向上应力最大，在最大水平地应力方向上位移最大，在最小水平地应力方向上位移最小，此时井眼将呈现为椭圆形。

3.2 水合物储层蠕变特征

对于非蠕变地层来说，井壁稳定的判断是一个瞬时或者短时效过程。但对于蠕变地层来说，井眼的变形是一个持续的过程，随着蠕变的进行井眼将发生缩径，同时井眼的形状也将不断改变。使用ABAQUS模拟蠕变过程时使用VSICO分析步进行分析，井眼打开240 h后其应力分布与位移分布分别如图6、图7所示。

图6 蠕变后地层米塞斯应力分布Fig.6 Mises stress distribution after creep

图7 蠕变后地层位移Fig.7 Displacement after creep

由图6可以看出，经240 h蠕变后，地层近井地带应力分布与井眼刚打开时的应力分布出现很大变化，选择井壁处1/4圆周路径上的点的应力绘制在同一坐标轴，如图8所示，起点为最小水平地应力方向上的点，既0°方向为最小水平地应力方向，90°为最大水平地应力方向。

图8 井壁应力分布图Fig.8 Stress distribution on the wellbore wall

由图8可以看出，井眼打开初期只发生弹性变形阶段，井壁处出现应力集中。井周一周在0°方向应力最大，在90°方向应力最小，这也符合非均匀地应力下直井井眼应力应力分布情况。此时在0°方向应力大小为11.7 MPa，90°方向应力大小为10.6 MPa，最大应力与最小应力比值约为1.1，这与初始设置的地应力比相同。当井眼钻开240 h以后应力分布发生了大的改变，总体上与初始时刻相比井壁处应力有所减小，这是由于井眼缩径所导致，井眼缩径将减轻应力集中程度。在最0°方向应力减小值最大，减小为为初始时刻的87%。同时井周应力差减小，最大应力处任为0°处但最小应力分布处为45°处，二者应力比值约为1.03小于初始时刻1.1，90°与0°方向应力比为1.005 5，井周应力分布已近似相等。

对于蠕变地层来说应力分布的变化必将导致蠕变速率的变化，从而影响井眼的变形。由图5和图7可以看出，蠕变发生前井眼受非均匀地应力作用在最大水平地应力出发生最大变形，在最小水平地应力处变形最小，选择蠕变发生前后1/4圆周地层变形如图9所示。

由图9可以看出，蠕变前后井壁位移发生变化，但总体趋势相同。0°方向(最小水平地应力方向)最小，90°(最大水平地应力方向)最大，但最小位移与最大位移比值却发生很大变化。蠕变发生前，井壁最大位移为1.43 cm，最小位移为1.25 cm，最大位移与最小位移的比值约为1.14；蠕变发生后井壁最大位移为6.79 cm,最小位移为6.7 cm，最大位移与最小位移的比值为1.01。这表明蠕变时井眼不同方位会以不同的蠕变速率进行缩径，图10所示为最大水平地应力方向和最小水平地应力方向井眼缩径随时间变化曲线。

图9 蠕变前后井周位移分布Fig.9 Wellbore displacement distribution before and after creep

图10 井眼缩径曲线Fig.10 Wellbore shrinkage curve

图11为最小水平地应力方向与最大水平地应力方向蠕变速率曲线，从图11中可以看出，不同方位井壁上的点初期缩径速率很大并随着时间增加而下降，这是由于水合物衰减蠕变的特征决定。图12所示为不同时刻最小水平地应力方向井眼缩径速率与最大水平地应力方向缩径速率比值，从图12中可以看出，初始时刻最小地应力方向缩径速率与最大主应力方向缩径速率的比值约为1.17，略大于初始最大水平地应力与最小水平地应力比值。

图11 蠕变速率曲线Fig.11 Wellbore shrinkage rate

图12 0°方向与90°方向缩径速率比值Fig.12 The ratio of shrinkage rate at the direction of 0° to 90°

这主要是受井眼的应力分布及井眼的初始位移所决定。从图8看出，弹性阶段结束后井壁出现应力集中，在最小水平地应力方向最大，最大水平地应力方向最小，二者比值约为1.1。水合物沉积物的蠕变速率受应力影响，在井眼不同方位地层在初始时刻二者所受应力比值为1.1，故二者的应变速率比值也约为1.1。同时由图9可知，初始时刻最大水平地应力方向井壁由于弹性变形产生初始缩径量大于最小水平地应力方向，因此二者的缩径速率比值将略大于1.1。

由图12可以看出，随着蠕变的进行，0°与90°方向的缩径速率比值不断减小逐渐趋近于1。这是由于初始时刻0°缩径量小于90°方向，但其缩径速率大于90°方向，这将导致0°与90°方向井眼缩径差距将不断减小，井眼形状的改变将影响经井眼周围应力分布，减小两个方向的应力差，使得缩径速率的比值随时间减小并趋近于1。这也解释了井眼打开240 h后，井眼周围最大应力与最小应力的比值远小于蠕变发生前并趋近于1。由此对于水合物地层的井眼形状变化如图13所示。

井眼打开初期在地应力作用下井眼由于弹性变形将呈现为椭圆形，长轴为最小水平地应力方向，短轴为最大水平地应力方向。随着蠕变的进行，井筒趋近于圆。井周应力分布初始时刻在最小水平地应力方向与最大水平地应力方向应力集中程度差异越大，随着蠕变的进行，应力分布将趋于均匀。

3.3 水合物储层井眼缩径影响因素分析

由于实际地层中水合物层的深度、地应力的比值水合物饱和度都影响水合物层的井周应力分布，因此建立一系列模型模拟不同因素对于水合物层井眼缩径规律的影响。

3.3.1 水合物层深度对井眼缩径的影响

选择地应力比为1∶1.1、钻井液密度为1.05 g/cm3、井眼直径为0.25 m模拟不同深度下水合物层井眼缩径。图14所示为不同深度下蠕变240 h后的井眼最大缩径量。

由图14可以看出，井眼的缩径量随水合物层所在深度的增加而增大，既水合物层所在深度越深，最终的井眼尺寸越小。这是由于地层深度越大，地应力的值越大，虽然地应力的比值相同，但应力差值增大；同时地应力与井眼液柱压力的差值越大从而导致井眼的缩径量随深度的增加而增加。

图14 不同深度井眼缩径量Fig.14 Wellbore shrinkage of the hydrate reservoir in different depth

3.3.2 地应力非均匀性对井眼缩径的影响

选择地层深度为500 m、钻井液密度为1.05 g/cm3、井眼直径为0.25 m模拟不同地应力比条件下水合物层的井眼缩径。图15(a)和图15(b) 分别为不同地应力比值下0°和90°蠕变240 h前后的井眼变形量。

图15 蠕变前后井眼缩径量Fig.15 Wellbore shrinkage before and after creep

由图15可以看出，地应力比值的增大会增大井眼的缩径程度。对井眼形状来说，蠕变发生前井筒微小变形是由地应力的挤压作用而造成的弹性变形，地应力比值越大，造成井口的不同方向变形不同程度越大，井眼向椭圆变化。但地层蠕变的发生后井眼的变形将不再是持续原本的弹性变形方向，而是从椭圆向圆发展，地应力差值越大蠕变同一段时间后地层的椭圆程度越大。

图16所示为初始时刻0°方向井眼缩径速率与90°方向比值。从图16可以看出，地应力比值越大，井眼不同方向的缩径速率比值越大，既蠕变初期的井眼收缩的不均匀性越大。一方面是由于地应力比值越大，初始时刻最小主应力方向应力集中程度较最大主应力方向越大；同时初始位移差也越大使得蠕变速率的不一致性随地应力的增大而增大。

图16 不同地应力比下的缩径速率比值Fig.16 Ratio of wellbore shrinkage rate at different in-situ stress ratios

3.3.3 水合物饱和度对眼缩径的影响

由于水合物饱和度与地层蠕变目前研究较少，但水合物饱和度对地层的其他力学性质的影响已有一个较为清晰的研究。因此，从弹性模量的角度研究水合物饱和度对地层井眼缩径的影响。实验得出水合物沉积物弹性模量和饱和度为0的沉积物弹性模量比值与其饱和度曲线，如图17所示[31]。

图17 水合物饱和度与弹性模量比值关系[31]Fig.17 Relationship between hydrate saturation and elastic modulus ratio[31]

拟合曲线可用三次多项式表示为

E/E0=8.85Sh3+1.54Sh2+1.249 4Sh+1

(9)

式(9)中：E0为水合物饱和度为0时沉积物弹性模量；E为某一饱和度下沉积物弹性模量。

建立水合物饱和度与水合物沉积物弹性模量的关系，由此模拟水合物层井眼缩径规律。图18所示为水合物饱和度与井眼缩径量关系。

图18 水合物饱和度与井眼缩径Fig.18 Hydrate saturation and wellbore shrinkage

随着水合物饱和度的降低井眼的缩径量逐渐增加，当水合物饱和度由30%变为10%时，蠕变缩径量增加约14%。这是由于饱和度的降低水合物层的弹性模量降低，初始弹性变形阶段弹性位移增大；其次弹性模量减小会增大井眼的变形能力，使井眼的缩径速率更快，最终使井眼的缩径量增加。

4 结论

(1)井眼打开时，水合物蠕变地层会在非均匀地应力作用下发生弹性变形，在最大水平地应力方向位移最大，最小水平地应力方向最小。同时井眼应力集中出现不均匀现象，在最小水平地应力方向应力最大，最大水平地应力方向最小，二者比值约为地应力比值。

(2)井眼发生蠕变缩径时，不同方位的缩径速率不同，最小水平地应力方向缩径速率最大，最大水平地应力方向缩径速率最小。

(3)随着蠕变的发展，井眼的缩径速率随时间发生指数衰减，当井眼钻开时间较长时，最大缩径速率与最小缩径速率比直接近1，最大缩径量与最小缩径量也接近1，井眼逐渐呈现为圆形等速缩径。

(4)蠕变发生一定时间后，井眼的应力分布将从不均匀分布向均匀分布发展，最终将消除地应力的不均匀对井眼应力分布的影响，从岩石受力角度井壁趋于稳定。

(5)随着水合物层所在深度的增加，井眼的蠕变量将增大；井眼的缩径量将随着地应力比值的增大而增大，同一时刻井眼的缩径速率、不均匀性也随地应力比值的增大而增大；随着水合物饱和度的降低，井眼的缩径量增加。