天然气水合物试采井筒-土壤三维非线性耦合模型研究

2022-04-13畅元江黄帅王康孙宝江李昊孙焕钊陈国明

中南大学学报(自然科学版) 2022年3期

畅元江，黄帅，王康，孙宝江，李昊，孙焕钊，陈国明

(1.中国石油大学(华东)海洋油气装备与安全技术研究中心，山东青岛，266580；2.中国石油大学(华东)石油工程学院，山东青岛，266580)

天然气水合物井筒稳定性是井筒抵抗结构性破坏、维持井筒功能的重要属性，是海洋水合物勘探开发过程中井下安全的保证。由于水合物具有埋藏浅、未成岩、强度低以及欠压实等特点[1-2]，试采时，水合物分解的最直接效应是近井界面地层力学特性弱化，造成井筒失稳甚至套管损坏，因此，有必要考虑水合物分解引起的储层物性参数在近井界面的时变特征，研究天然气水合物三维井筒-土壤非线性耦合模型，并针对水合物储层井筒的稳定性进行模拟分析。近年来，国内外学者对深水水合物开发过程中的井筒、井壁稳定性进行了大量研究。FREIJ-AYOUB等[3]根据摩尔-库仑失效准则，利用单向达西流动模型分析水合物地层孔隙压力变化，考虑孔隙压力变化对应力场的影响，建立了水合物地层钻井井壁稳定性模型；RUTQVIST 等[4-5]考虑天然气水合物降压开采过程中的渗流场、应力场和温度场的耦合作用，对地层以及井筒的稳定性进行了分析，提出垂直井射孔段是最易出砂和应力集中的位置；KIM等[6]利用FLAC3D 软件进行了油藏规模的热-液-力耦合模拟，考虑井筒组件的尺寸和沉积物与井筒相互作用，分析了Ulleung 盆地天然气水合物在降压开采过程中高压储层的力学响应和井筒稳定性；NING等[7-9]研究了钻井液侵入水合物地层条件下钻井液参数(即温度、密度、盐度)和初始储层条件(即渗透率和饱和度)对井筒稳定性的影响。程远方等[10-11]将水合物分解效应融合到渗流场与岩土变形场的耦合作用中，建立了天然气水合物藏气、水两相非等温流固耦合数学模型，并开发了有限元程序；张若昕[12]从裸眼段井筒及套管段井筒2个方面对天然气水合物开采过程中的井壁稳定性进行了研究；毕圆圆[13]建立了二维水合物地层-水泥环-套管力学模型，分析了孔隙压力和温度耦合作用对井筒应力的影响。上述研究为开展水合物试采井筒稳定性分析提供了理论依据，然而，这些研究多以二维为主平面模型为主，不能反映由于水合物分解导致的地层沉降对井筒的影响，且没有考虑多场耦合下管土相互作用的时变特征。为此，本文作者基于海洋水合物井筒面临的复杂环境载荷，考虑水合物分解效应，以井筒-储层为研究对象建立水合物试采井的三维井筒-土壤非线性耦合模型，研究井筒-储层之间随时间变化的相互作用，揭示深海水合物储层渗流场、应力场和井筒内外压力场的耦合作用对井筒的力学响应规律，以便为我国海域天然气水合物资源安全高效试采提供参考。

1 水合物分解多场耦合模型

1.1 Mohr-Coulomb 土壤本构模型

Mohr-Coulomb 模型可以较好地描述土体在弹性或塑性变形上的破坏机理[14]。土体在弹性阶段的变形屈服方程如下：

式中：F为土体剪切力，MPa；p为平均主应力；MPa；q为偏应力，MPa；Rmc为度量π 平面上屈服面形状的参数；φ为土体摩擦角，即q-p面上Mohr-Coulomb屈服面的倾斜角；C为土体黏聚力，MPa。

当土体位于塑性变形阶段时，在应力空间子午线上，流动势G为双曲线函数，可表示为

式中：Ψ为土体剪胀角；C0为土体塑性变形为0时的初始黏聚力，MPa；ε为子午线平面上的偏心率；Rmw为流动势G在π平面上的形状参数。

1.2 地层多孔介质渗流-变形模型

储层岩石的受力主要为岩石骨架承受的有效应力和储层孔隙流体的孔隙压力。根据有效应力原理，土体内某一微元的静力学平衡微分方程为

式中：τzx，τxy和τyz分别为xz，xy和yz平面的切应力；P为孔隙压力；分别为x，y和z轴方向的单位孔隙压力；σ′x，σ′y和σ′z分别为x，y和z轴方向的有效应力；γ为土的容重。

只考虑弹性变形情况时，土体形状变化和有效应力改变需遵守胡克定律要求，其表达式如下：

式中：E为土体骨架的弹性模量；v为泊松比；G为剪切模量，G=E/[2(1+v)]；εx，εy和εz分别为x，y和z方向的应变分量。

根据小变形理论，应变与位移的表达式为

式中：u，v和w分别为x，y和z方向上的位移分量。

在水合物开采过程中，孔隙中水和气体的压力远大于毛细管压力，在此忽略毛细管作用，认为空隙中水和气体压力相等，考虑水合物分解对地层渗透系数的影响，采用比奥固结理论求解渗流-变形耦合模型：

单元体内水量的变化率等于土体积的变化率，由达西定律可得补充方程：

式中：γw为水的容重；k为渗透系数；εV=εx+εy+εz，为体积应变量。

1.3 水合物储层物性参数动态变化模型

由于水合物分解效应，储层物性参数将随着水合物的开采发生变化，根据水合物储层的渗透率、孔隙度与水合物的饱和度和土的体积应变建立函数关系[15-16]，得到渗透系数与孔隙度的动态计算模型：

式中：SH为水合物储层的饱和度；K0为绝对渗透率，即储层SH=0 时的渗透率，10-3μm2；ρ为流体密度，kg/m3；g为重力加速度；η为动力黏滞性系数；φ0为初始孔隙度；N为地层渗透率递减指数，由水合物赋存形式决定，通常取2～15。

在水合物开采过程中，水合物储层的强度和胶结特性等会伴随水合物的分解逐步降低和减弱，表现为水合物储层弹性模量减小以及黏聚力呈线性软化[17]，在钻采过程中，水合物储层弹性模量、黏聚力与水合物饱和度相关的计算模型可以表示为

式中：E0为水合物未分解时的弹性模量，GPa；SH0为水合物储层的初始饱和度；Δφ=φ0(SH0-SH)，为水合物分解引起的地层孔隙度变化量；ξ为实验拟合参数，通常取9.5；ψ取为1.2。

对于我国南海海域广泛存在的胶结或骨架支撑型水合物，在实际模拟过程中需要不断地更新调整相应的力学参数和渗流参数。式(1)和式(2)描述的Mohr-Coulomb 本构模型可作为水合物地层的土壤模型；式(7)可为Mohr-Coulomb土壤模型计算提供劲度矩阵，用来控制三维渗流过程中地层的变形；式(8)～(11)可为水合物开采模拟中储层物性参数的变化提供计算依据。

2 三维井筒-土壤非线性耦合模型建模方法

2.1 水合物试采井井身结构

由于海洋水合物储层埋藏浅、未固结，井身结构采用如图1所示的三开方式，储层段采用套管不固井完井[18-19]，通过降低储层压力，进行水合物降压开采。

图1 天然气水合物试采井井身结构Fig.1 Well structure of gas hydrate production test well

本文选取受水合物分解影响最显著的储层段井筒作为研究对象，其力学模型如图2所示。通过有限元软件建立三维实体模型，研究水合物在试采过程中井筒与地层随时间的变化。为此，编写了水合物分解多场耦合模拟程序，进而对试采过程中水合物分解进行模拟。

图2 储层段井筒力学模型Fig.2 Wellbore mechanics model of reservoir section

2.2 水合物分解多场耦合模型的建立

在水合物试采时，储层孔隙流体流失，孔隙压力水平降低，使水合物相变分解，导致岩石骨架的应力状态发生变化，相应的储层渗透率、孔隙度等物性参数随之改变，这些由孔隙流体流失造成的储层变化又会对孔隙流体的渗流过程产生影响，这是一个复杂的非线性耦合过程。本文利用建立的渗流-变形耦合模型与水合物储层物性参数动态变化模型，对ABAQUS 进行二次开发，基于Fortran 语言对模型积分点的场变量进行重新定义，开发编写了水合物模拟求解器。考虑水合物分解引起的储层应力场、渗流场和井筒内外压力场的特征参数在近井界面区域的时变特征，对水合物开采过程中多场耦合作用下的井筒受力进行模拟，计算流程如图3所示。

图3 多场耦合分析计算流程Fig.3 Flow chart of multi-field coupling analysis

进行多场耦合分析时，首先定义储层物性参数随场变量的变化；然后，对压力场与渗流场进行计算，通过调用场变量子程序USDFLD 对水合物分解区域的储层物性参数进行实时更新，以此实现管土相互作用边界的动态变化；根据更新的边界对下一个增量步进行求解，直至达到设定的模拟分析时长为止。

2.3 三维井筒-土壤非线性耦合模型实现

按以下步骤建立三维井筒-土壤非线性耦合模型(见图4)：

图4 三维井筒-土壤非线性耦合模型建模流程Fig.4 Modeling flow chart of 3D nonlinear coupling wellbore-soil model

1)确定研究位置的水合物地层物性参数和井筒结构参数，建立有限元几何模型。

2)建立渗流-变形耦合模型和水合物储层物性参数动态变化模型。

3)根据建立的水合物分解多场耦合数学模型编写场变量子程序，模拟水合物的分解。

4)对模型赋予初始预定义场(地应力场和孔隙压力场)以及管土相互作用边界。

由于材料刚度存在差异，井筒与水合物地层之间的界面处产生摩擦力，这些摩擦力导致井筒产生轴向应力。在定义管土相互作用时，采用接触面法，套管表面为主控表面，土体表面为从属表面；考虑管土接触界面之间存在黏聚力，管土间法向接触定义为修正的硬接触，切向接触定义为具有各向同性的库仑摩擦，采用罚函数法求解。

3 算例分析

3.1 三维非线性耦合模型基本参数

本研究采用的数值模拟参数来源于我国神狐海域某水合物试采井，水深为1 230 m，上覆土层厚度为200 m，基本参数如表1所示[20-21]。

表1 天然气水合物藏地层物性参数Table 1 Formation physical parameters of gas hydrate reservoir

在建立管土耦合模型时，为了避免模型端部边界效应的影响，模型参数选取如下[22]：地层设置为长×宽×高为5 m×10 m×17 m 的长方体，其中上覆层厚度为5 m，水合物储层厚度为10 m，下伏层厚度为2 m，水泥环的长度为5 m，套管贯穿整个地层的长度取17 m。井筒基本参数如表2所示。

表2 井筒基本参数Table 2 Basic parameters of wellbore

3.2 基于三维非线性耦合模型的井筒稳定性分析

利用图4所示流程图进行建模分析。由于套管和水合物地层的结构及载荷呈对称性，且在井筒与地层胶结良好的情况下，建模型时可沿加载方向取模型的1/2 进行分析，如图5所示。在模型的上表面施加海水及上覆岩层等效压力15 MPa，储层土体受重力的作用，对土体上、下以及四周表面赋予初始渗透压力，在套管内壁施加生产压力；在模型底部对土体和套管施加U3方向的位移约束，在土体外表面施加U1和U2方向的位移约束，以模拟周围地层对模型的约束作用。给模型对称面施加对称约束。对土体赋予地应力、初始孔隙压力、孔隙比等预定义场后，进行地应力平衡计算，结果见图6。由图6可知储层的位移数量级为10-7m，远小于地应力平衡要求的位移数量级10-5m，表明符合地应力平衡条件，可以进行下一步计算。

图5 载荷和边界条件Fig.5 Loading and boundary conditions

图6 地应力平衡结果Fig.6 Result of geo stress balance

天然气水合物的开采与常规油气的最主要区别在于开采过程存在水合物的分解效应。图7和图8所示分别为以10 MPa 的生产压差试采60 d 后孔隙压力和水合物饱和度分布情况。从图7和图8可以看出：随着开采生产过程进行，储层中的孔隙压力从生产区域向周围地层逐渐增加，越靠近井眼位置，孔隙压力越小，孔隙压力变化梯度越大；井筒周围储层的水合物饱和度随着开采逐渐降低，且水合物完全分解区域与未分解区域的过度区域较窄。

图7 孔隙压力分布云图Fig.7 Nephogram of pore pressure distribution

图8 水合物饱和度分布云图Fig.8 Nephogram of hydrate saturation distribution

在试采井系统中，井筒是否稳定关系到试采作业能否安全、顺利进行，影响井筒稳定性的关键因素为试采井井筒中套管的受力。储层段井筒应力情况如图9所示，图中正值表示套管受到的是拉应力，负值表示套管受到的是压应力。从图9可以看出：由于水合物的分解，水合物储层下沉和下伏地层向上隆起，使得储层段井筒产生较大的轴向压应力，且套管的最大轴向压应力位于储层段套管中部即降压开采的位置；套管最大环向应力同样位于储层段套管中部降压开采位置，套管受到的径向应力可忽略不计。套管Mises应力的变化主要受轴向应力影响，其变化趋势与轴向应力变化趋势基本一致，表明轴向应力可以作为井筒稳定性的判据。

图9 储层段井筒套管应力分布Fig.9 Stress distribution of wellbore casing in reservoir section

3.3 井筒稳定性影响因素分析

3.3.1 初始水合物饱和度对井筒稳定性的影响

以生产压差为13 MPa，试采时间为30 d为例，分别模拟水合物初始饱和度为20%，40%和60%时井筒的稳定性，套管轴向应力情况如图10所示。

由图10可知：在水合物初始饱和度为20%，40% 和60% 时，套管的最大轴向应力分别为-472.76，-417.00 和-381.10 MPa(负号表示轴向应力方向向下)，表明在相同的试采时间和温压条件下，储层初始水合物饱和度越小，储层段套管受到的轴向应力就越大，井筒的稳定性越差。这主要是因为水合物初始饱和度越高，储层渗透率相对越低，从而降低了孔隙压力的降压速率，减慢了水合物分解的进程。

图10 不同初始水合物饱和度对套管轴向应力分布的影响Fig.10 Influence of different initial hydrate saturations on the axial stress distribution of casing

3.3.2 生产压差对井筒稳定性的影响

设定储层初始水合物饱和度为40%，试采生产时间为10 d，分析生产压差为5，8，10，12 和和13 MPa 共5 种试采条件下井筒的稳定性，套管轴向应力如图11所示。

由图11可知：套管的轴向应力随着生产压差的增加而增加，其中，当生产压差由5 MPa增加到8 MPa 时，套管的最大轴向应力由154.1 MPa 增加到246.9 MPa，增大了60.22%，可见生产压差越大，储层水合物的分解速率越快，范围越广，相应的储层段套管因地层有效应力变化受到的压力也更加大。这与常规疏松砂岩油气藏开采所得到的结果类似。

图11 不同生产压差对套管轴向应力分布的影响Fig.11 Influence of different production pressures on axial stress distribution of casing

3.3.3 试采时间对井筒稳定性的影响

设定储层初始水合物饱和度为40%，生产压差压为10 MPa，在试采时间分别为1，5，10，30和60 d 时，套管轴向应力如图12所示。参照API标准，考虑到水合物储层生产段套管存在射孔情况，选取N-80 套管屈服强度552 MPa 的60%即331.2 MPa 作为参照[23]，分析不同生产压差下套管轴向应力与生产时间的关系，如图13所示。

由图12和图13可知：随着生产压差提高和试采时间延长，井筒的轴向应力逐渐增加，且轴向应力在试采初期增加较明显，随着试采进行逐渐趋于平稳；当生产压差为13 MPa与12 MPa时，最大轴向应力在试采初期就超过了套管的屈服强度(331.2 MPa)，若生产压差为10 MPa，则安全试采仅能维持15 d；若生产压差不超过8 MPa，则安全试采可达60 d以上。