郭旭洋，金 衍，林伯韬，卢运虎，訾敬玉

(中国石油大学(北京)石油工程学院，北京 102249)

天然气水合物储量巨大，并具有清洁高效的特点，在中国国南海已经开展2次试采，其中第2次为水平井试采[1]。水平井降压开采过程中的工程扰动会诱发水合物分解，导致井周沉积物力学性质劣化，井周应力时空演化特征复杂，存在出砂和海底滑坡的可能性[2]。因此明确试采过程中的沉积物变形破坏机制和流-固-热-化多场耦合机制对于南海水合物安全、高效钻采具有意义。水合物沉积物胶结性较差，水合物的分解对其力学性质影响显著。研究人员采用人工合成样品和现场保压取芯等手段开展了大量三轴力学试验，确定多种水合物沉积物试样的强度参数和变形特性，为现场尺度建模分析提供关键力学参数。水合物和沉积物基质共同决定水合物沉积物的力学性质，但是基质占据主导作用。增加水合物饱和度或减小有效围压均能够使沉积物应力-应变关系由应变硬化转向应变软化。沉积物内聚力受水合物饱和度影响作用明显，而内摩擦角仅在水合物饱和度较高时受水合物饱和度影响较大[3-9]。针对中国南海北部深水浅层天然气水合物沉积物的力学试验显示，有效围压增加的情况下，沉积物性质会呈塑性，有效围压和水合物饱和度综合影响沉积物力学性质，导致水合物沉积物钻采过程中地层失稳风险较大[10]。力学试验获取的水合物沉积物力学性质为更大尺度的水合物沉积物钻采过程建模分析提供了基础数据支持，结合水合物分解动力学、质量守恒和能量守恒等，能够对水合物沉积物钻采过程中的多种物理场演化开展定量建模分析。通过有限元建模分析发现，较小的井眼半径更易导致水合物水平井发生井壁失稳，这种影响随钻井时间加剧[11]。水合物相变、传热传质和地层变形的耦合建模有助于确定水合物分解前缘和解释出砂和稳产时间短等现象，为商业化开发提供理论参考，而生产压力、渗透率、初始饱和度和温度等都与地层变形有一定联系[12-13]。基于TOUGH等平台，能够定量表征水合物开采过程中的传热传质问题和水合物分解规律[14]。在此基础上通过孔隙弹性力学理论进行流固耦合研究，认为固体力学场的耦合方法对模拟结果会产生显著影响。全耦合对收敛性和准确性的提升有较大帮助，而如果采用顺序耦合则需要更严密的耦合条件以确保数值求解的稳定性和准确性[15-17]。笔者根据南海水合物沉积物力学特征，建立水平井降压开采过程的流-固-热-化多场耦合有限元模型，分析沉积物内的变形、破坏和海底沉降规律。

1 数学模型建立

以水合物沉积物为对象构建多物理场耦合模型，表征沉积物内水合物分解、热传递、气水渗流和变形破坏过程。根据质量守恒，多孔介质内的流动过程可以表示为

(1)

式中，φ为孔隙度；Si为流体饱和度；ρi为流体密度，kg/m3；vi为流体速度，m/s；si为汇/源质量速度，kg/(m3·s)；下角i代表相态类型。

沉积物中流体渗流速度低，根据达西定律[18]，多孔介质内流体流速表示为

(2)

式中，k为绝对渗透率，m2；Kri为相对渗透率；μi为黏度，Pa·s；pi为孔隙流体压力，Pa；g为重力加速度，m/s2。

水合物降压开采过程中的传热问题需要考虑沉积物固相骨架内的热传导、多孔介质内流体的热对流和水合物分解产生的吸热过程[19]。该过程中的传热问题可以表示为

(hC+hE)=Q.

(3)

式中，E为比内能，J/kg；S为饱和度；h为对流传热项；Q为热汇/源，J/(m3·s)；下角L、G、H分别代表液相、气相和水合物相3种相态。

降压开采过程中，水平井井底压力降低，导致天然气水合物稳定区域平衡被破坏，水合物开始分解[20-21]。根据水合物分解动力学模型，降压过程中水合物分解速率表示为

RMH=-kdMMH(pe-pG)As.

(4)

式中，RMH为分解速率；kd为反应速率，Pa·s；MMH为水合物摩尔质量，kg/mol；pe为相平衡压力，Pa；pG为气相压力，Pa；As为反应比表面积，1/m2。

由于温度和压力的变化会导致多孔介质及多孔介质内流体发生体积变化，具体可以表示为

(5)

(6)

式中，χ为压缩系数，1/Pa；β为膨胀系数，1/K；T为温度，K。

由于降压，沉积物内水合物分解，水合物饱和度下降，导致沉积物力学性质出现劣化[22]。

水合物分解会改变内聚力，c可以表示为水合物饱和度的函数

(7)

式中，c0为不考虑水合物部分的沉积物内聚力，Pa；θ内摩擦角，(°)；SMH为水合物饱和度；α和β为系数。

为表征降压开采过程中的沉积物力学响应，需要根据动量平衡对应力张量建立表达式

(8)

式中，σ为总应力张量；ρb为体积密度，kg/m3。

根据无穷小变换，应变张量ε和位移矢量u的关系可以表示为

(9)

在多孔介质中，本构关系表示为

δσ=C:δ(ε-εp)-αδpI.

(10)

式中，C为弹性张量；εp为塑性应变；α为耦合系数；I为二阶单位张量。

利用Mohr-Coulomb和Drucker-Prager模型表征弹塑性变形下的剪切破坏。屈服函数f和Drucker-Prager模型的塑性势函数g可分别表示为

(11)

(12)

式中，I1为有效应力第一不变量；J2为有效偏应力第二不变量；βf、κf、βg和κg为破坏包络线系数[23-25]。

Mohr-Coulomb模型可表示为

(13)

(14)

(15)

(16)

2 水平井降压开采模型

根据上述控制方程，建立水合物沉积物降压开采模型。由于水平井钻采是南海试采工作的一项重要内容，也是提高产气效率的关键手段，因此建立单个水平井模型和双水平井模型，其中双水平井的水平段深度不同。图1为单个水平井筒降压开采模型和双水平井筒降压开采模型。由于水平井筒沿x方向贯穿模型，井筒长度较长，可以采用平面应变假设简化数值运算，仅分析y-z截面内的二维流-固-热-化耦合问题。对于单井模型，井筒位于y-z截面正中(25 m，25 m)；对于双井模型，井筒分别位于y-z截面的(25 m，15 m)和(25 m，35 m)。对于二维y-z平面模型，y方向为SHmax，z方向为Sv，x方向施加Shmin。在厚度为50 m的沉积物层内，假设地层压力梯度和地应力梯度的影响较小，地层压力为14 MPa，垂向主应力为16.6 MPa，水平最大主应力为14.76 MPa，不随深度而增大。沉积物密度为2 300 kg/m3，初始饱和度为40%，孔隙度为0.15，渗透率为1.8×10-2μm2，地层温度为284.15 K，干湿导热系数分别为1和3 W/(m·K-1)，降压生产压力为3 MPa。

图1 单个水平井筒和双水平井筒降压开采模型Fig.1 Sketch of models for depressurization in a single horizontal well and two horizontal wells

3 结果分析

根据图1所示模型，开展水平井筒降压开采过程中的数值模拟，分别针对单个水平井筒和双水平井筒开展分析，模拟降压开采时长为15 d。降压开采模型能够表征该过程中井周沉积物内的水合物分解范围、分解吸热引发温度变化、孔隙压力降低和弹塑性演化，并能够预测水平井产能。在二维y-z平面内，由于数值解关于z轴对称，仅展示y轴正半轴数值结果。

图2为采用单个水平井和2口水平井降压开采5和15 d时的井周孔隙流体压力分布。近井区域压力降最明显，这是由井筒降压开采使用的较低生产压力造成的。随着距离向远场移动，压力越接近初始压力。随着降压开采时间的增加，压降前缘的波及范围更广。对于单井模型，降压开采5 d时远场压力仍等于初始压力14 MPa，但开采15 d后远场压力也开始下降。相对应地，采用双井模型能够增加压降前缘的波及范围。对比开采5 d时的单井模型和双井模型的孔隙流体压力结果，发现双井模型井周压降范围更大，压降更加剧烈。开采15 d后双井模型的远场区域压力普遍降至低于11 MPa，单井模型原厂区域压力仍可达13 MPa，说明采用两口水平井降压开采能够更有效地建立地层内的压降，为甲烷流动提供促进作用。

图2 单个水平井筒和双水平井筒降压开采5和15 d孔隙流体压力分布Fig.2 Pore pressure distribution after 5 and 15 d of depressurization for single- and two-well scenarios

图3为单井模型和双井模型内降压开采5和15 d时的水合物饱和度分布特征，水合物饱和度为0的区域表示水合物完全分解。结果显示水合物分解区域近似呈圆形，这是因为y-z平面内的井周压降前缘也近似呈圆形，压降导致水合物稳定区的相平衡遭到破坏而出现水合物分解，分解区域受到压降前缘影响较为直接。根据图3数据，5 d时单井水合物分解前缘距水平井筒扩展1.01 m，双井为0.96 m；15 d时水合物分解前缘扩展的距离在单井模型中为1.74 m，在双井模型中则为1.58 m。说明增加水平井数量和增加降压开采时间均会对水合物分解范围产生影响。随着降压开采时间的增加，水合物分解前缘向远场扩展，15 d内能够扩展1.58～1.74 m。

图3 单个水平井筒和双水平井筒降压开采5 d和15 d水合物饱和度分布Fig.3 Hydrate saturation distribution after 5 and 15 d of depressurization for single- and two-well scenarios

图4为单井模型和双井模型降压开采5和15 d后的温度分布。由于水合物分解是吸热过程，导致分解区域及邻近沉积物内温度降低。降温前缘也近似呈圆形，与压降前缘和水合物分解前缘形状类似。造成这一现象的原因除了压降与水合物分解直接相关外，也与式(1)、(3)有关，两者均为抛物型微分方程，呈现的传质传热特征具有类似性。图4表明，单井模型中的降温前缘在5 d时为1.25 m，在15 d时为2.71 m；双井模型中的降温前缘在5 d时为1.21 m，在15 d时为2.67 m。这说明开采时间仍是影响降温范围的主要因素，水平井数量对降温范围的影响相对较小。

图4 单个水平井筒和双水平井筒降压开采5和15 d温度分布Fig.4 Temperature distribution after 5 and 15 d of depressurization for single- and two-well scenarios

图5 水平井累积产气对比Fig.5 Comparison of cumulative gas production between studied scenarios

根据降压开采多物理场耦合模拟可以进行水平井产气预测(图5)。对比了单井模型水平井产量、双井模型井1(z=35 m)产量、双井模型井2(z=15 m)产量和双井模型两口水平井平均产量。总体累产曲线显示单井累产气差异不明显，这是因为两口水平井井距较远，尽管井周压降区域互相波及，但水合物分解区域和温度降低区域仍相隔一段距离，而水平井产气主要是由水合物分解区释放的甲烷提供，因此在井周水合物分解范围相似的情况下，获得的水平井累产气结果也较为近似。比较单井模型和双井模型内两口井的累积产气曲线，发现15 d内累产气由高到低排序为双井模型井1>单井模型>双井模型井2。分析发现这一关系能够和三口水平井井位深度建立起联系，累产气最高得双井模型井1深度最小，单井模型其次，双井模型井2深度最大。由于在y-z平面内多场耦合模型考虑了重力的作用，导致降压开采后深度更大的区域压力相对更高、压降相对更弱。据此，在深度更大的区域进行水平井降压开采形成的压降差异会略小于在深度较小的区域进行降压开采，一定程度上减缓水合物分解和甲烷释放，减小了对应水平井累产气。此外，还计算了双井模型两口井累产的平均值，发现其数值略高于单井模型。这说明在本研究设定的条件下，增加水平井数量能够适当提升单井累产效果。4条累产曲线均呈现出开采初期产气效率很高、开采一段时间后产量急剧下降的特点，与前期国内外部分试采工作结果具有一定的一致性。

降压开采诱发的沉积物变形与破坏也是需要研究的重点之一。图6为采用单井模型和双井模型降压开采5和15 d后的沿垂直方向的总应力分布情况。由于近井区域受到降压开采影响，孔隙流体压力下降明显，导致垂直方向总应力在井筒和近井区域处于较低值。在同一时间，增加水平井数量进一步加剧了垂直方向总应力的降低。随着开采时间增加，垂直方向总应力降低的区域逐渐增大，而且远场区域也出现一定的垂向总应力降低。通过与图2对比，发现垂直方向总应力的演化规律与孔隙流体压力额演化规律具有一定关联性，只是因为在孔隙弹性力学理论中，孔隙流体压力是总应力中各正应力分量的重要组成部分。在15 d时，双井模型内远场垂向总应力比单井模型更低，也是由双井模型远场孔隙流体压力更低造成的。

图6 单个水平井筒和双水平井筒降压开采5 d和15 d垂向应力分布Fig.6 Vertical stress distribution after 5 and 15 d of depressurization for single- and two-well scenarios

表征降压开采诱发的弹塑性变形能够帮助预测气井出砂和海底沉降。图7为单井模型水平井、双井模型中井1和井2同深度的塑性变形时间演化规律。由图7(a)看出，在降压开采起始时，仅在井壁处存在塑性应变，储层沉积物均无塑性应变发生。降压开采0.5～1 d后，由井壁至沉积物层内出现随距离逐渐减小的塑性应变，这是由于近井区域压降剧烈，沉积物固体骨架变形程度更大。0.5～1 d时，由于水合物分解尚不明显，水合物饱和度对塑性变形的影响不显著。降压开采5 d后，由于水合物分解前缘的扩展，分解区向未分解区过渡的塑性变形呈现出突变特征，分解区展现的塑性应变高于未分解区。从0.5～15 d的结果显示，水解物分解前缘和塑性应变突变点具有较强相关性，分解前缘外的未分解储层内的塑性应变则由近及远平滑降低。该结果说明水合物分解区塑性变形特点明显，更易出砂和诱发地层失稳。

图7 不同模型井塑性变形时间演化Fig.7 Temporal evolution of plastic deformation in different model wells

由图7(b)、(c)看出双井模型中井1和井2各时间点的塑性变形分布特征与单井模型结果比较接近，均呈现出水合物分解前缘与塑性应变突变点吻合的情况。5、10和15 d的塑性应变在10 m以外出现小幅度升高的现象，这是受到了模型的边界条件的影响。在y-z平面的右侧边界设置为固定位移边界，双井模型的压降和变形程度和范围更大，导致近井区域压缩变形程度较高，由于右侧边界为固定边界，在靠近右侧边界的部分区域出现拉伸变形。说明在本研究中模型尺寸和边界条件对塑性变形求解过程有影响。

4 结 论

(1)由于降压开采中的传热和传质过程控制方程均具有抛物型微分方程的特征，沉积物层内的孔隙压力、温度和水合物分解前缘的时空演化规律具有较强的相关性。

(2)沉积物层内的总应力时空演化规律受孔隙压力演化的影响明显，近井区域的总应力数值小于远场受压降影响较小区域的总应力数值。这一现象能够通过孔隙弹性力学理论进行解释。

(3)在水平井降压开采的多场耦合研究时，重力作用对水平井产能、力学响应和孔隙压力分布的影响较为明显。降压开采后，重力作用导致深层水平井压降程度低于浅层水平井，一定程度上抑制水合物分解前缘的扩展。本研究数模显示，浅层水平井累积产气量更高。

(4)水合物分解会明显地劣化沉积物力学性质，诱发较强的塑性应变。在降压开采初期水合物分解前缘未形成时，沉积物层塑性应变不明显；随着降压开采时间的增加，水合物分解区力学性质劣化显著，塑性应变明显高于未分解区,内聚力等强度参数变低,在水合物分解前缘附近的塑性应变会出现明显的阶梯式下降特征。此外，建模分析时的模型尺寸和边界条件设置均在一定程度上可能影响塑性变形求解过程。