基于cohesive单元海域天然气水合物储层水力压裂模拟

蒋亚峰，田英英，李小洋，韩泽龙，张欣，赵明，刘聃，尤志伟，梁金强，黄伟

（1.中国地质科学院勘探技术研究所，河北 廊坊 065000；2.清华大学，北京 100083；3.中国地质调查局广州海洋地质调查局，广东 广州 511466）

0 引言

天然气水合物是由水和甲烷气体在低温高压条件下组成的结晶化合物［1］，普遍分布在深海海底沉积物和陆地永久冻土中［2］，据初步预测，赋存于水合物中的天然气资源量达到3×1015m³，其中可开采量达到3×1014m³，被视为是21世纪一种重要的潜在的清洁替代能源［3-5］。由于海域天然气水合物富集程度较低，分布范围广，因此，与其他常规化石能源相比，从天然气水合物储层中开采天然气难度极大。目前，世界上仅有中国和日本2个国家进行了海域天然气水合物试开采。2013年，日本首次在Nankai海槽开展了海域天然气水合物试采作业，经过6 d作业，产气量达11.95×104m³，平均日产气量2×104m³［6-8］。2017年在相同海域，根据第一次试采的相关经验，再次开展了第二次试开采。其中，第一口井由于井底严重出砂只开采了12 d，产气4.1×104m³；第二口试采井持续开采了24 d，累计产气22.25×104m³［9］。2017年，我国在南海神狐海域成功完成了天然气水合物试采工程，实现了60 d内持续产气30.9×104m³［10］。上述3次海域水合物试采工程均采用单垂直井降压法开采。2020年，我国首次采用水平井降压法对南海神狐海域天然气水合物进行第二次开采，累计产气达149.86×104m³，日均产气3.57×104m³［11］。虽然通过多次尝试在总产气量和产气速率上取得了极大的进步，并验证了相关钻井技术、开采工艺和相关装备的可行性，但是与商业化开采标准还有一段距离。

与日本水合物试采区的砂质地层对比，我国南海海域水合物主要赋存于泥质粉砂地层中，其赋存环境更加复杂，泥质含量高导致渗透率低，严重制约了降压开采的效果，从而使得常规的降压开采难以获得理想的产气量［12］。其实，很早之前便已有研究人员提出将水力压裂应用于海域水合物储层改造中，利用水力压裂形成的高导流裂缝增加了降压效果和水合物分解面积，从而提高产气量［13］。张永勤等［14］提出利用旋转喷射钻井将改造剂注入到水合物储层中形成蜂窝状支柱，从而起到提供流动通道和支撑地层的作用。彭赛宇等［15］提出通过对储层灌注特殊的泡沫水玻璃浆液改造剂，形成方状骨架结构，从而达到提高储层渗透率和增大产气量的效果。KONNO等［16］针对含水合物砂质沉积物，开展了不同侧重的水力压裂试验。含甲烷水合物砂质沉积物水力裂缝表现出拉伸破坏特征，扩展方向垂直于最小主应力方向，与完整固结岩石压裂行为相似。虽然相关理论提出了很多，但是，针对水合物储层的水力压裂裂缝的起裂机理、扩展规律方面的研究较少。因此，研究水合物储层水力压力过程中裂缝的起裂、扩展和形态，对后续分析储层水合物分解速率以及推动海域天然气水合物的安全高效开采具有重要意义。基于此，本研究在前人理论的基础上，利用abaqus软件建立基于cohesive单元的二维水力压裂数值模型，研究压裂过程中的储层基质变形、裂缝起裂和扩展、裂缝形态等参数，以验证采用cohesive单元法分析水力压裂及水合物储层采用水力压裂进行改造的可行性。

1 水力压裂基本方程

1.1 cohesive单元简介

在Abaqus软件中，通过在原始单元内插入cohesive单元，将其分为上下两面，其中上表面和下表面的节点距离表示裂缝宽度，中间面的节点表示压裂切向流动的状态，如图1所示。该方法可以模拟压裂过程中的储层基质变形、裂缝起裂和扩张、裂缝形态等参数，并且在页岩气和煤层气等非常规能源储层水力压裂模型分析中得到了成功运用［17-19］，为后续研究其他特殊地层水力压裂裂缝扩展规律提供了基础［20］。

图1 cohesive单元示意Fig.1 Schematic diagram of a cohesive unit

1.2 水力压裂裂缝扩展方程

CHEN［20］和ZOU等［21］提出，采用牵引扩展定律描述水合物储层发生破坏前的线性弹性行为、裂缝起裂和扩展过程，如图2所示。随着流体的注入，储层介质牵引力不断增大，压裂储层发生线性弹性行为，裂缝尖端开始不断延伸，当注入流体达到内部弹性模量Tmax时，扩展达到临界值δ0。当其大于δ0时，牵引力随着裂缝的扩展而减小即发生损伤演化行为，直到分离达到临界值δf（压裂的裂纹尖端），此时牵引力或注入流体强度消失，裂纹尖端完全停止延伸［22］。

图2 牵引扩展定律Fig.2 Law of traction propagation

在单一类型荷载下，GONG等［23］也指出牵引扩展定律在分层过程中具有线性荷载和线性退化两大特征：

式中：k0——cohesive单元表面刚度；δ0、δf——分别为单元界面的初始形变和完全破坏的相对位移；D——截面破坏损伤因子，其取值在0和1之间，当D为0时，表示材料完好无损，当D为1时，表示材料已完全破坏。

应用二次应力准则则可以判断水力压裂裂缝是否开裂，即当cohesive网格单元任一方向的三向应力与其对应的强度极限之比等于1时，水合物储层开始破裂：

水合物储层在完全起裂后，牵引扩展定律中的损伤演化模型如下：

同时通过线性位移扩展准则，可以推出损伤因子D的公式：

1.3 流体流动方程

图3为水合物储层压裂流体的切向流动和法向流动示意图。由图3可知，根据水合物储层压裂的流固耦合关系模型，可以看出压裂液主要由切向流动和法向流动2部分组成，分别导致裂缝的扩展和压裂液滤失。

图3 流体的切向流动和法向流动示意Fig.3 Tangential and normal flow patterns of fluid

流体切向流动方程：

流体法向流动方程：

式中：q——切向流量；Δp——单元长度方向压力梯度；ω——裂缝宽度；μ——压裂液粘度；qt、qb——分别为单元上下表面的压裂液流速；ct、cb——分别为上下表面的渗滤系数；pt、pf、pb——分别为上、中、下表面的流体压力。

1.4 孔隙介质固体颗粒的控制方程

水力压裂过程会改变地层中岩石的受力场，对于岩石这类的多孔介质，其应力平衡方程可以应用虚功原理描述：固定时间内岩石的虚功，等于作用在整个岩石上的体力与面力产生的虚功的叠加［24］。

岩体介质的平衡方程为［25］：

式中：δε——虚应变率矩阵；t——面力矩阵；δv——虚速度矩阵；f——体力矩阵；dV——单元体。

2 水力压裂有限元模型

本文的研究对象为水合物储层，与常规油气储层相比岩性较软，其特殊的赋存形态形成了复杂的地层结构，加之不稳定性、沉积物和应力状态等多种因素影响，水力压裂过程极其复杂，目前相关的研究均处于起步阶段。因此，本文借助泥岩、页岩水力压裂相关模型基础，并只选取双向单缝为研究目标，简化模型。

2.1 数值模拟参数选取

2.1.1 储层弹性模量和泊松比的选取

Lijith等［26］等借助大量的拟合数据结果，提出了水合物沉积物的弹性模量计算公式：

式中：E——水合物储层弹性模量；——有效围压；SH——水合物饱和度；a——代表原始含天然气水合物沉积物的初始刚度的常数，大小取决于沉积物的组成及地层压力；b——代表受水合物形态、水合物形成条件及饱和度影响的常数。

WILLIAM等［27］通过单轴压裂试验测得水合物储层等温杨氏模量7800 MPa，泊松比0.317；DVORKIN等［28］测得天然气水合物储层等温杨氏模量为7900 MPa。通过对比，本次选取弹性模量7800 MPa，泊松比0.317。

2.1.2 渗透率参数的选取

DAI等［29］基于Kozeny-Carman模型提出了水合物储层的渗透率模型，该模型与Kumar等［30］、Johnson等［31］和Kleinberg等［32］的实验结果取得较好的吻合，其公式为：

式中：K′—水合物储层的渗透率；K0——水合物中沉积物的渗透率；SH——水合物饱和度。

2.1.3 孔隙比的计算

水合物储层的有效孔隙度受水合物饱和度的影响，具体公式为：

式中：φ′——水合物储层的有效孔隙度；φ——无水合物状态下的孔隙度。

孔隙比是指水合物储层中有效孔隙占水合物与沉积物体积之和的百分比，其计算公式为：

式中：e——水合物储层孔隙比。

2.2 二维水力压裂有限元模型

为了更好地分析储层在水力压裂过程中的裂缝扩展及应力、应变变化规律，对二维水力压裂数值模拟做出如下假设：

（1）水合物储层为均质地层；

（2）孔隙压力为3.3 MPa；

（3）储层上下、左右边界固定，不产生位移。

同时，分别选用20 m×20 m、100 m×100 m两种尺度的模型进行模拟，其中20 m×20 m仿真模型主要用于研究裂缝及周边的应力、应变，100 m×100 m模型主要用于研究裂缝扩展方向、延伸机理方面的研究。此外还设置了水平方向和竖直方向的最大地应力σ=5 MPa与最小地应力σz=4 MPa，模型中心位置设置为压裂液的注入点，cohesive预设单元为裂缝发展方向，对称分布在井筒两端，如图4所示。2个模型均以（x，z）为坐标的格式离散划分成100×100=10000个cohesive网格块。模型选取的主要参数，如表1所示。

表1 水合物储层水力压裂模型主要参数Table 1 Main parameters of the hydrate reservoir hydraulic fracturing model

3 结果与分析

本文分别对20 m×20 m、100 m×100 m两个模型尺寸下的裂缝扩展进行了数值模拟，得出了水合物储层水力压力过程中裂缝扩展的最大半长及最大宽度，通过2种尺寸下应力、应变、孔压的分布状态分析最优模拟结果，最后对比不同时间段裂缝扩展的速率，得出最优压裂时长，为后续的水合物水力压裂模拟提供指导。

3.1 裂缝扩展尺寸

模拟结果，如表2所示，2种模型尺寸下裂缝的半长均为6 m，大尺寸模型裂缝最大宽度为5.8 mm，小尺寸模型裂缝最大宽度5.5 mm。2组模拟结果基本吻合，验证了cohesive单元法应用于水合物水力压裂研究的可行性，为研究水合物储层裂缝变化、扩展机理提供了良好的基础。

表2 模拟结果Table 2 Simulation results

3.2 仿真结果分析

模型建立的过程中，我们分别设立了（x，z）2个方向的边界值，即水平方向和竖直方向的最大地应力σx=5 MPa与最小地应力σz=4 MPa。在模拟的过程中，随着压裂液的不断注入，裂缝不断扩张的同时，其内部的压裂液沿法向切力的方向作用于裂缝内表面，并不断向外辐射。如果应力位移传递量设置过小，很容易导致应力不能有效传递，在边界位置产生变形，影响模拟结果的准确性。

图5为2种模型尺寸裂缝扩展模拟结果，从图5可以明显看到20 m×20 m裂缝扩展应变图中应力在传播的过程中在边界位置立即出现变形，而100 m×100 m模型中应力呈扇形均匀的扩散，模拟效果更加真实；同时，可以明显看出20 m×20 m裂缝扩展尖端应力和孔压图中受位移传递量的限制存在发散现象，而100 m×100 m模型中裂缝延展尖端的应力和孔压变化不受边缘的限制，模拟的区域图形呈现规律性和对称性。所以对比模拟结果后裂缝最大宽度在100 m×100 m模型下的5.8 mm更加准确。

图5 裂缝扩展模拟结果Fig.5 Simulation results of crack propagation

3.3 裂缝宽度随时间的变化规律

图6为注入点裂缝宽度随注入时间与压力的变化特征，从图中可以看出裂缝的发生形态与注入压力的关系，随着注入压力的急剧增加，储层开始起裂，裂缝瞬间张开，此时注入压力最大，达到了25 MPa，然后再经历一段波动后趋于稳定，注入压力约为18 MPa。裂缝宽度在起裂后，以小幅度波动的方式缓慢变大。原因在于，当开始向储层注入压裂液时，在注入点附近很容易形成憋压导致裂缝瞬间增大。起裂发生后，裂缝瞬间沿预设单元扩展方向裂开。注入压力出现微降，在水平地应力的作用下，裂缝宽度出现短暂的收缩，后续随着注入液的补充，开始逐步趋于缓慢的增长状态，形成这种“阶梯式”的裂缝扩展形态。

图6 裂缝宽度随时间与压力的变化规律Fig.6 Variation law of crack width vs time and injection pressure

4 结论

本文通过利用abaqus软件建立基于cohesive单元的二维水力压裂数值模型，研究压裂过程中裂缝起裂和扩展、裂缝形态等参数，得出了水合物储层水力压裂裂缝扩展的最大半长及最大宽度。具体结论如下：

（1）采用100 m×100 m和20 m×20 m大、小两尺寸模型分别建立相同参数下的水合物二维水力压裂裂缝扩展模型，最终得出2组基本相同的数据结果：裂缝长度均为6 m，最大宽度分别为5.8、5.5 mm，论证了cohesive单元法分析水合物水力压裂方案的可行性。

（2）在二维水合物储层水力压裂模拟中，通过对2种尺寸模型数值模拟结果的分析，在利用cohesive单元法分析水合物水力压裂模型中应尽量构建尺寸较大的模型，这样更利于应力的传递及裂缝扩展的模拟，得出的实验结果更加准确。

（3）裂缝在扩展的过程中会伴随着阶梯性延展，一方面，由于刚开始注压时，较高的排量很容易在注入点附近形成憋压导致裂缝瞬间增大，并出现了短暂的波动。同时，随着预设单元扩展方向迅速裂开，此时裂缝中的压裂液内压会出现微降，在水平地应力的作用下，裂缝宽度出现短暂的收缩的现象，后续随着注入液的补充，开始逐步趋于缓慢的增长。这种急速扩张—收缩的现象形成了水合物储层“阶梯式”的水力压裂扩展规律。该研究能够为海域天然气水合物储层水力压裂提供一定的理论指导，提高水合物储层的渗流能力。