基于嵌入黏聚单元法的页岩储层压裂缝网扩展规律

位云生 林铁军 于 浩 齐亚东 王军磊 金亦秋 朱汉卿

1.中国石油勘探开发研究院 2.“油气藏地质及开发工程”国家重点实验室·西南石油大学 3.中国石油勘探与生产分公司

0 引言

体积压裂技术已经成为有效开发页岩气的主要手段[1]，在相同有效改造体积范围的储层中形成的裂缝数量越多、缝网越复杂，改造效果就越好，产能也就越高[2-5]。缝网形态及其复杂程度与岩石力学性质、天然裂缝空间展布、压裂施工参数等都有着密切的关系。因此，体积压裂过程中，对天然裂缝发育形态的认识与人工裂缝的控制是压裂改造的关键和难点问题。

目前国内外学者对不同储层条件或不同工况条件下水力裂缝的扩展开展了大量研究。室内试验方法可以直观地获取压裂裂缝的宏观和微观扩展规律[6-10]，然而当岩石离开地下原位环境，其结构及力学属性可能发生变化[11-14]，导致试验结果偏离真实情况[15-17]。且由于较大尺度的试验成本高昂，运用数值方法对压裂过程中裂缝扩展机理及规律进行模拟成为了主要手段。目前国内外学者针对水力压裂中裂缝扩展的数值模拟研究方法主要包括线弹性断裂力学法（LEFM）、边界元法（BEM）、离散元法（DEM）、扩展有限元法（XFEM）以及黏聚单元法（CZM）。Yue等[18]和Marco等[19]使用LEFM法对脆性岩石水力裂缝扩展中断裂力学行为进行表征。然而，LEFM法忽略了裂缝尖端的奇异性，不适用于准脆性或非均质性较强地层的裂缝扩展过程模拟。Zou等[20]、周彤等[21]、李玉梅等[22]、Gordeliy等[23]利用离散元的方法（DEM）将裂缝的扩展简化成线性弹簧的断裂失效行为，从而分析了不同层理弱面发育密度、强度、裂缝间距、力学参数以及压裂施工参数对水力裂缝展布规律的影响。但是基于DEM法的裂缝扩展只能沿着刚性块体的边界而不能很好地处理连续体扩展问题，同时计算消耗量极大。Olson[24]和胥云等[25-27]认为BEM法可以通过在定义域边界和裂尖划分、加密单元并对边界进行插值离散来避免尖端奇异性，在裂缝扩展模拟上具有独特优势，但是BEM法不能充分考虑裂缝内流体流动而只能假设裂缝内为均匀孔压[26-27]。Dahi-Taleghani等[28]、盛广龙等[29]提出了基于XFEM法的复杂水力裂缝模式扩展模型作为设计工具，可用于在复杂扩展条件下优化处理参数。XFEM法的优点是裂缝扩展不需要重置网格和预置路径，可以有效模拟水力裂缝沿任意路径的扩展，但是无法构建原始地层中的天然裂缝。方修君等[30-32]证明了CZM法模拟沿任意路径的裂纹扩展问题的可行性。Guo等[33]采用流动—应力—渗流耦合的CZM方法，考虑了盖层—储层的多层结构，对页岩气压裂过程水力裂缝扩展进行了模拟研究。Dahi-Taleghani和Yu等[34-37]提出不同性质的黏聚单元对不同方位的天然裂缝进行表征。但是，CZM法必须预置裂缝的扩展路径，难以模拟人工水力裂缝在基岩中沿任意方向起裂扩展的过程。目前多数成果主要针对的是单条水力裂缝的起裂与扩展，或者两条裂缝的交汇形式的研究，缺少了对缝网模拟及其复杂程度的探索与研究。因此需要一种既可以准确表征天然裂缝展布，而且允许水力裂缝沿任意方向扩展的方法，以实现模拟体积压裂缝网扩展演化的目的。

本文以四川盆地南部奥陶系五峰组—志留系龙马溪组某页岩气井为研究对象，提出了“批量嵌入黏聚单元”的建模方法，利用Python编程实现每个基岩单元的外边界批量插入零厚度的黏聚单元，用于表征人工水力裂缝的潜在扩展路径以及天然裂缝展布。因此，该方法既利用了黏聚单元描述天然裂缝准确性的优势，又保持了人工水力裂缝在基岩中任意路径扩展的随机性。鉴于此，建立了考虑水力裂缝与天然裂缝交汇的水力压裂储层—裂缝有限元模型，模拟了大型水力压裂人工水力裂缝与天然裂缝竞争起裂、交错扩展而形成缝网扩展的演化过程。获取了水力裂缝缝网形态为“星形”，并根据缝网形态拟合了压裂体积的计算公式，模拟分析了压裂施工排量对缝网扩展演化的影响，并通过对压裂后实际产量对比，验证了压裂体积计算的可靠性，为准确认识页岩压后缝网形态提供了重要的理论依据。

1 黏聚单元法模拟裂缝扩展原理

裂缝的起裂和扩展可以视为黏聚材料的渐近脱胶过程。在水力压裂过程中，压裂液流体进入胶结面中抵消原地应力后并超过了岩石材料的胶结强度时，对于水力裂缝起裂、扩展及与天然裂缝交汇行为也可以采用黏聚方法进行表征，水力裂缝扩展的潜在路径和天然胶结裂缝等特征也可以视为胶结在一起的两个面，其中胶结厚度很小，接近于零厚，胶结面就开始按照牵引—分离损伤理论进行损伤演化直至完全打开。胶结界面的打开和分离对应着水力裂缝的起裂和扩展以及与天然裂缝的交汇作用[1]。

为了解决模拟人工水力裂缝在基岩中沿任意方向起裂扩展的过程这一难题，实现对人工水力裂缝随机起裂扩展行为的模拟，提出了一种基于Abaqus有限元软件，利用Python二次开发技术的批量嵌入黏聚单元方法（BCZM，Batched Cohesive Zone Method）。该方法可以根据基础网格坐标，利用Python编程在每个基岩单元的外边界都批量嵌入黏聚单元，用于表征水力主裂缝的潜在扩展路径，如图1-a所示。其中蓝色与紫色线条分别表示储层中天然存在的不同方向和性质的天然裂缝，而所指出的红色黏聚单元即为储层中人工水力裂缝可能扩展的路径，其余绿色单元均为其潜在的扩展路径，能有效表征大型水力压裂人工水力裂缝扩展的随机性。红色和绿色单元，均为相同单元，并赋予了相同的参数设置，红色部分仅为裂缝可能的扩展路径的一个示意，所有绿色单元都可以是潜在的扩展路径，裂缝实际的扩展情况还与注入泵压、原始地应力场等多种因素有关。因此，该方法既利用了黏聚单元的描述天然裂缝准确性的优势，又保持了人工水力裂缝在基岩中扩展的任意随机性。

图1-b展示了基于黏聚单元方法的不同黏聚单元交汇扩展方法示意图。黏聚单元具有位移节点和孔隙压力节点，红色填充圆为位移节点，黄色填充圆为水平方向压裂的压力节点，灰褐色填充圆为垂直方向压裂的压力节点，黄色方块为水平方向的黏聚单元，绿色方块为垂直方向的黏聚单元。在水力压裂过程中，压裂液的孔隙压力通过孔隙压力节点进行传递，当压力节点值超过岩石基岩强度或天然裂缝胶结强度时，位移节点就会逐渐张开，以此来表征裂缝起裂过程。在压力节点传递过程中，黏聚单元逐渐被打开，就形成了裂缝扩展的过程，并且可以通过对不同的黏聚单元赋值不同强度参数分别表征岩石基岩及不同性质的天然裂缝，进而对页岩的各向异性特征进行表征。

图1 批量嵌入黏聚单元模拟人工水力裂缝起裂扩展方法示意图

牵引—分离损伤理论（Traction-Separation Law，TSL）是用于描述物体损伤及其断裂行为的方法[38-39]，如图2所示。可以采用TSL对基岩和不同天然裂缝的黏聚强度进行力学表征。当黏聚单元没有损伤时，其可以被视为一个双面胶结的整体，随着黏聚单元内部压力逐渐升高直至达到黏聚单元的初始损伤强度，黏聚单元就开始逐渐的出现损伤，其损伤系数与分离距离呈线性关系，当黏聚单元的开启位移超过其失效位移时，黏聚单元变成完全失效并自动删除，两个胶结的面被完全分离开，变成两个孤立的壁面，几何上表现为裂缝的起裂。

图2 牵引—分离损伤理论示意图

牵引—分离模型首先假设了线弹性行为，然后是损伤的起始和发展。弹性行为以弹性本构矩阵的形式表示，该弹性本构矩阵将界面上的名义应力与名义应变联系起来。名义应力是在每个积分点上的力分量除以原始面积。名义应变是在每个积分点上除以原始厚度的分离值，并对横向剪切分量施加一些平均值。名义牵引应力矢量t由三个分量组成（二维问题中有两个分量）。牵引—分离损伤准则定义了在裂缝尖端黏结层的黏结界面之间的本构关系。可以通过胡克定律简单描述沿内聚区的弹性行为：

式中Eij(i,j=n,s,t)表示材料在对应三个方向的黏聚刚度，Pa/m；tn、ts和tt分别表示法向和两个切向的牵引，Pa；δn、δs和δt分别表示对应的分离量，m。

当黏聚单元损伤后，弹性参数开始退化，退化程度用损伤参数（D）表示（0～1）。标量损伤参数表示在损伤开始后进一步加载时，岩石的整体损伤从0单调演变为1。黏聚单元牵引—分离损伤的应力分量受以下因素的影响：

随着孔隙压力增加到一定程度，多孔岩石中的破坏开始并发展，从而导致液压驱动的裂缝扩散。在裂缝内部，假定流动状态遵循泊松流动（层流）。裂缝内流体流包括两个分量：裂缝内的切向流和垂直于裂缝表面往地层的法向滤失。

对于牛顿流体，在其进入黏聚单元时的切向流动方程为：

式中q表示沿裂缝的质量流量，N/(m3·s)；d表示裂缝开度，m；kt表示切向渗流系数，m/s；∇p表示压力梯度，Pa/m；μ表示流体黏度，Pa·s。

法向流动方程为：

式中qt、qb分别表示裂缝上壁面和下壁面的法向流速，m/s；ct、cb分别表示从裂缝中流体通过裂缝上壁面和下壁面的滤失系数，Pa·s；pt、pb分别表示裂缝上壁和下壁单元表面的孔隙压力，Pa；pi表示裂缝中的流体压力，Pa。

然后，岩石基体的控制方程包括流体流动和岩石变形的耦合，表示为：

式 中σij、分别表示总应力和初始应力的分量，Pa；E表示弹性模量，Pa；ν表示泊松比；表示应变张量的分量；表示 Kronecker's delta函数 ；α表示Biot系数；pw、分别表示孔隙压力与初始孔隙压力，Pa。

在运用黏聚单元法模拟水力压裂原理进行理论计算或建模时，黏聚单元损伤前的流动方式即是在地层中的渗流，且黏聚单元损伤后的流动方式分为切向上的平板流动和法向上往地层中的滤失。

2 水力压裂缝网演化模型建立

在对页岩储层进行压裂改造的过程中，人工裂缝扩展不仅受到页岩层理、天然裂缝等先天地质条件的影响，而且还受到井口压裂注入排量等作业工况的影响。解析方法无法耦合多种因素定量刻画水力压裂缝网演化过程与形态，故需采用基于黏聚单元的有限元方法，对薄层状海相页岩储层人工裂缝的动态扩展进行模拟分析。

2.1 储层基本情况

以四川盆地南部五峰组—龙马溪组某区块页岩气井为研究对象，该区块目的层埋深2 500～3 000 m，从下至上有效储层包括五峰组、龙一11、龙一12、龙一13、龙一14小层，层理和天然微裂缝发育，尤其是龙一11小层，如图3所示。水平井靶体和轨迹位于龙一11小层中部，各小层储层参数见表1所示。

表1 五峰组—龙马溪组各小层储层参数表

图3 川南龙一11小层页岩层理与天然微裂缝分布图

水平井体积压裂作业参数如表2所示。室内压裂大型物理模拟实验表明，页理导致水力裂缝垂向穿透受限，平面上离射孔点越近，页理缝越长；现场微地震监测数据表明，纵向水力裂缝高度介于35～40 m，裂缝向上沟通了龙一14小层；但非放射性示踪剂测井解释结果表明，支撑裂缝仅9～12 m，且由于滑溜水黏度低、携砂能力弱，支撑剂集中在水力裂缝下部；现有主流软件采用矩形裂缝形态，拟合生产动态数据，动用高度介于10～20 m，向上仅沟通了龙一13小层的底部。

表2 压裂施工参数表

2.2 黏聚单元有限元模型建立

基于实际储层中层理、天然裂缝分布情况及各小层储层参数数据，采用黏聚单元法，建立储层—裂缝平面应变有限元模型。由中深层致密砂岩储层的压裂实践可知，天然裂缝不发育的块状储层中水平井分段压裂形成的是沿着水平最大主应力方向延伸的多条近似平行的主裂缝。海相页岩储层不同的是，发育层理和天然微裂缝多种弱面，在动态裂缝延伸的过程中，当泵入液体的压力超过弱面破裂压力时，弱面会打开，并沿弱面延伸一段距离，不断减弱压裂液传递的能量，故远离射孔点的水力裂缝长度和高度会快速减小。模型中水力裂缝、层理、天然微裂缝的分布如图3所示。图3-c为水力裂缝与层理、天然微裂缝交汇示意图以及天然微裂缝与层理的接触关系示意图。综合考虑以上因素，所建立的模型横向范围为300 m，纵向范围为40 m，如图4所示。

图4 页岩储集层黏聚单元地层—压裂缝网模型示意图

对图4中模型的A、B、C、D四个边界进行约束，作为模型中地层的远场边界条件。然后再在模型中添加应力场，以还原地层中原始地应力情况。最后再在图中注入点位置进行压裂液注入模拟。

3 水力压裂缝网演化模拟研究

为了准确描述页岩储层水力压裂缝网的演化过程及形态，分析井口压裂液注入排量对压裂缝网形态的影响，运用已建立的黏聚单元有限元模型进行数值模拟计算，并对计算结果进行分析。

3.1 水力压裂缝网演化过程及形态分析

水力压裂过程中，从井口不断向地层中高压注入压裂液，使得地层人工裂缝中流体压力增大，裂缝得以在地层中扩展延伸，形成如图5-a所示的孔隙压力云图。通过图5-a可以看出，井眼周围孔压区域随着注入压裂液时间而不断增长，最大孔压覆盖区域形状也在不断演变。

从图5-a可得，层理和天然裂缝开启对压裂液能量起到分流作用，由于水平层理薄弱面起裂压力低，沿层理面的裂缝横向扩展速度最快，对纵向缝高方向上分流作用明显，故压裂液能量衰竭迅速，水力裂缝高度受限，且远离射孔点，缝高急剧降低。随着井周缝网形态逐渐演化，最终形成一个主体区域为星形的缝网区，缝网区域外的边界线为星形包络线。

图5 水力压裂缝网扩展演化过程孔压及缝网形态云图及压裂缝网长度随时间变化曲线图

水平井靶体位置和轨迹在龙马溪组最底部的龙一11小层，上部地层为龙马溪组地层，下部地层为间隔较薄的五峰组页岩，与宝塔组灰岩接触。由于五峰组地层破裂压力比页岩高得多，人工裂缝无法或难以开启与扩展，故裂缝形态呈现出上半部分面积大于下半部分面积的特点。

通过对缝网纵横向长度测量，当注入排量为12 m3/min时，缝网的最大纵向高度逐渐从零增加到24.1 m，缝网的最大横向长度逐渐增加，这表明，压裂过程中水力主裂缝不断扩展，同时层理也在不断开启和扩展。压裂液进入地层后，首先沿着层理与人工主裂缝方向快速扩展，即在横向和纵向两个方向快速延伸。

另外，以水平井眼位置为圆心，做内切于星形区域上半部分边缘的内切圆，令其半径为压裂缝网区域有效半径。由图5-a可知，缝网的有效半径达到8.82 m。

从模拟压裂作业过程时间来看，不同时刻形成的压裂缝网横向长度、纵向高度以及压裂缝网区域有效半径随时间变化存在较大差异，沿层理方向的横向长度随着时间增长速度远大于纵向，有效半径在一定压裂作用时间后趋于稳定，如图5-b所示。

水力主裂缝与水平层理交汇扩展演化缝网在水力主裂缝扩展延伸过程中，与人工主裂缝交汇的天然裂缝被开启并延伸，裂缝之间不断沟通和交错，地层逐渐形成压裂缝网。

由图5中压裂缝网结构演化过程可以看出，由于靶体层位龙一11小层的水平层理发育密度最大，故形成的水力裂缝密度最大，其他小层的水力裂缝密度沿着远离水平井筒的方向逐渐减小。同时由于龙一11小层的岩石强度相较于其他层位更低，故横向上水力裂缝扩展长度也是最长的。当水力主裂缝的横向长度达到93.6 m时，其缝内流体压力不足以打开其他层理或沟通更多的天然裂缝，呈现单缝扩展。由于在水力压裂后期裂缝呈现单缝扩展，所形成的裂缝对于页岩地层改造增产效果不大，故在缝网面积计算时不考虑，缝网区域的有效横向长度忽略该段区域。通过测量计算，有效改造缝网区域长轴长度为48.36 m，并且关于井眼所在位置的垂线呈现出横向两侧对称关系，两侧横向长度基本相等。

模拟所得该工况下水力压裂最终形成的缝网改造区域横截面如图6所示，其中有效长轴（横向）长度48.36 m，有效短轴（纵向）长度24.1 m，星形横截面面积 421.85 m2。

通过拟合水力压裂最终缝网横截面轮廓曲线，可以定量计算星形形状缝网面积，围成区域为最右端位置为A，最左端位置为－A，有效短轴最上端位置为B，最下端位置为C。

根据图6缝网区域轮廓拟合曲线，可以得出缝网上、下区域星形包络线轮廓边线拟合曲线表达式：

图6 缝网横截面示意图

式中mi，ki，ni表示与实际地层情况以及施工工况有关的参数。mi，ni与A，B，C有以下关系：

根据拟合曲线公式可以进一步推导出星形轮廓面积计算公式：

例如：根据图6缝网轮廓线，可得该模型注入排量为12 m3/min的星形缝网拟合轮廓边线表达参数为m1=17.02，k1=0.087，n1=－2.80，m2=－7.52，k2=0.086，n2=1.56，计算拟合星形轮廓面积为426.93 m2，与模拟所得压裂改造缝网区域面积421.85 m2基本相符。

3.2 注入排量对缝网演化的影响分析

为研究井口压裂液注入排量对页岩储集层水力压裂缝网演化的影响，进一步分别模拟了注入排量为 6 m3/min、18 m3/min、24 m3/min 三种工况下裂缝扩展演化的情况。

在设置井口压裂液注入排量为6 m3/min时，人工水力压裂缝形成改造缝网区域如图7-a所示。从图中可以看出，最终形成的星形压裂缝网区域最大横向长度为65.49 m。压裂缝网区域有效横向长度为41.8 m，有效纵向长度为15.18 m，模拟所得星形裂缝网区域面积为339.26 m2。相比较12 m3/min排量工况，水力压裂缝网面积减少了，缝网区域有效横向长度和有效纵向长度都减少了，说明注入排量对缝网最终形态影响很大。

图7 不同排量下人工水力压裂缝形成改造缝网区域及压裂缝网星形面积随井口注入排量变化关系曲线图

当压裂液注入排量提高到18 m3/min时，人工水力压裂缝形成改造缝网区域面积示意图如图7-b所示。从图中可以看出，最终形成的星形压裂缝网区域最大横向长度为97.2 m，已经贯通模型边界并逃逸。而井眼周围压裂缝网区域有效横向长度为54.55 m，有效纵向长度为32.81 m，模拟所得星形裂缝网区域面积为493.43 m2，持续注入压裂液不会继续增大缝网有效纵横向长度。较之于12 m3/min，压裂缝网面积进一步增大，缝网密度与复杂程度进一步增高，说明压裂增产效果更好。

当井口压裂液注入排量为24 m3/min时，人工水力压裂缝形成改造缝网区域面积示意图如图7-c所示。从图中可以看出，最终形成的星形压裂缝网区域最大横向长度为98.6 m，压裂液更早贯通模型边界而被天然层理所捕获。压裂缝网区域有效横向长度为50.64 m，有效纵向长度为30.79 m，模拟所得星形裂缝网区域面积为478.64 m2。较之于12 m3/min，压裂缝网面积增加，但较之于18 m3/min有一定的下降趋势。

井口压裂液注入排量设置为24 m3/min，已经大于正常压裂作业时使用的压裂液注入排量大小，所以在此工况下，较大量的压裂液进入地层，较大的压力使得压裂液快速地沿人工水力主裂缝穿过地层而被天然层理捕获，如图7-b所示。故此时压裂液不能有效地沟通与开启天然裂缝，从而使得此工况下的压裂缝网面积以及压裂效果不及18 m3/min排量工况。

综合各个排量工况的模拟结果，可以得到压裂缝网星形面积以及拟合轮廓曲线所围成面积随井口注入排量变化关系曲线如图7所示。可以看出，在井口压裂液注入排量较小时，实际压裂星形面积即压裂改造区域面积也较小；随着注入排量的增加，星形面积不断增加，压裂改造效果不断提高；当注入排量为18 m3/min时，模拟压裂星形面积达到最大值493.43 m2，拟合星形面积达到最大值 499.37 m2，之后随排量的增加实际压裂星形面积呈现出下降的趋势。

对于提出的运用拟合轮廓曲线以计算压裂缝网区域面积的方法，在不同井口压裂液注入排量工况下，对应的压裂区域顶点位置A、B、C的数值如表3所示，对应的拟合曲线参数如表4所示。

表3 不同工况下拟合轮廓顶点位置数值统计表

表4 不同工况下拟合轮廓曲线参数及拟合面积统计表

根据对页岩压裂缝网演化过程与形态的模拟研究，进而得出图7所呈现的井口压裂液注入排量与实际压裂改造缝网面积之间关系。随着压裂液注入排量的提高，更多的天然裂缝在压裂液压力的作用下与人工水力裂缝相互沟通，开启起裂并得以扩展延伸，水力压裂改造区域面积随之增大，缝网密度也更为密集。但当压裂液注入排量超过某一临界值时，压裂液会快速形成井周缝网，并沿强度最低的层理面形成人工主裂缝且不断延伸，由此便不能有效开启天然裂缝，使得压裂缝网区域面积减小，压裂改造效果变差。

如表4所示压裂缝网区域拟合轮廓公式在不同工况下的参数值，其中实际地层中压裂缝网有效横向长度与有效纵向能够延伸的长度决定了mi、ni取值，ki和实际地层情况与压裂作业参数有关，可以看出，压裂作业时井口压裂液注入排量越大，压裂缝网轮廓拟合曲线中ki的值越大，对应的曲线曲率半径越小。在实际压裂作业时，需要根据地层情况，合理调整压裂施工参数，使得压裂缝网轮廓达到预期缝网演化形态及有效缝网面积。

4 结论

1）压裂液进入地层后会沿人工主裂缝方向进行延伸，由于地层对裂缝扩展的阻力以及压裂液压力的共同作用，与人工主裂缝交汇的天然裂缝被开启起裂并延伸，地层逐渐形成压裂缝网，最终井眼横截面压裂缝网区域形状呈星形，并在到达一定的有效半径后趋于稳定。

2）根据模拟所得压裂缝网结构图，可以测量计算得到星形缝网区域的面积和轮廓边线，提出拟合轮廓曲线表达式，可作为评价水力压裂改造效果的一个新方法，为更准确估算压裂后产量提供理论依据。

3）压裂液注入排量对压裂缝网的演化形态具有重要影响，随着注入排量增加缝网面积先增加后趋于降低，故在实际压裂施工中，可根据地层情况获得合理压裂液注入排量。