陈志强， 李凤霞， 潘林华， 周彤

(中国石化石油勘探开发研究院， 北京 100083)

21世纪以来，页岩油、页岩气进入快速发展期，成为全世界油气产量的重要组成部分[1-3]，其中复杂裂缝压裂(体积压裂)技术是实现页岩油气商业开发的关键技术之一[4-7]。室内裂缝扩展研究[8]以及现场微地震信号监测[9-10]显示，页岩水力压裂形成的裂缝不是简单的双翼型，而是呈现网状分布，主裂缝周围存在大量次级分支裂缝。压裂液和支撑剂在主/次裂缝交汇处发生分流，如何保障支撑剂顺利进入次级分支裂缝并进行有效支撑是实现体积压裂的关键。

以往体积压裂相关研究主要侧重于如何形成复杂裂缝，对于复杂裂缝形成后“支撑剂去哪了，能否进入分支缝？”的问题还不明晰[4，11]。数值模拟、室内实验提供了部分定性认识，但缺乏归纳总结形成定量规律，特别是还没有工程可用的数学模型。设计人员在进行压裂液排量、支撑剂粒径等参数选择时主要依赖现场经验，缺乏理论依据，导致在施工过程中出现主裂缝加砂困难、分支缝支撑不足等问题，最终影响体积压裂改造效果，因此亟需开展相关理论研究[4]。

支撑剂沉降、运移实验研究目前主要集中在单一裂缝[12-13]，相关学者分析了加砂浓度、压裂液黏度以及速率等对于支撑剂运移与展布规律的影响。针对复杂裂缝，虽有相关报道[5，7，14]，但受限于场地与实验装置，主要考虑一条主裂缝与一条分支缝构成的二级裂缝系统。然而，页岩水力压裂形成的复杂裂缝网络往往呈现多级、多尺度特征，主裂缝(100 m)与最小尺度裂缝(100 μm)跨越4～5个量级并通过不同级别的分支缝沟通互联，且小尺度缝对于页岩等致密储层同样发挥重要作用，但是这种多级、多尺度特征实验研究往往难以直接考虑。

支撑剂输运数值模拟方法主要包括传统的双流体模型[15]以及离散元模型[16]。双流体模型将离散的支撑剂颗粒简化为连续的流体相进行处理，这种方法整体计算量小，在模拟简单裂缝支撑剂运移时计算精度较高，但缺点是难以准确描述支撑剂在主/次裂缝交汇处的输运行为，且无法考虑支撑剂粒径的影响。为了更准确模拟支撑剂在复杂裂缝中的输运行为，有学者提出计算流体力学(computational fluid dynamics， CFD)与离散单元法(discrete element method，DEM)耦合的方法[16]，该算法采用CFD模拟压裂液流动过程，并计算支撑剂所受的流体作用力，采用DEM模拟该作用力下支撑剂的输运行为。CFD与DEM耦合算法的优势是不需要将支撑剂简化为连续的流体相，而是直接模拟支撑剂颗粒在流体中的运动，因此计算精度更高，但缺点是计算量相对较大。

综上所述，现提出工程可用数学模型，具体模型包括两个方面。首先，主裂缝以及各级分支缝内压裂液流速对于支撑剂输运行为具有重要影响，因此，本文研究关注的第一个问题是在注入排量(速率)已知的情况下，压裂液如何在多级裂缝系统中进行液量(速率)分配(液量分配模型)。第二个关键问题是，在已知主裂缝与各级分支缝压裂液排量的情况下，支撑剂能否在裂缝交汇处发生转向(支撑剂转向输运条件)，以及有多少支撑剂可以发生转向(砂量分配模型)。旨在通过理论分析的方法对复杂裂缝支撑剂输运规律进行研究，为现场施工提供参考。

1 复杂裂缝液量分配模型

1.1 理论模型

虽然目前石油行业对于非常规储层开发中复杂裂缝的重要性已达成共识，但是针对复杂裂缝系统液量分配模型的研究还相对较少。侯磊等[4]参照限流压裂思路提出了一个简化模型，该模型假设上级裂缝中液体流量都分配到了下级裂缝，分流比与下级裂缝截面积2/3次方成正比，并采用Fluent模拟结果对于模型进行了验证，结果显示该模型预测的流量分布与Fluent模拟结果定性上趋势一致，但定量上误差较大。因此有必要进一步分析前人模型误差来源，提出更为准确的液量分配模型。

在模拟多孔介质渗流时，一个经典的方法是孔隙网络模型(pore network)[17]，该方法将多孔介质简化为一系列“孔”与“喉道”的集合，“孔”之间通过“喉道”相连，并采用解析方法计算每一个“喉道”的流动阻力，最后利用基尔霍夫定律(流入节点和流出节点的流量和为零)联立求解，得到多孔介质内流场分布。该方法既考虑了多孔介质微观细节，通过“孔”“喉道”简化又提高了计算速率，因此在渗流领域广泛应用。

借鉴孔隙网络模型的思路，如图1所示，将已知形态的复杂裂缝进行网络结构提取，其中裂缝端点与裂缝交汇点作为“节点”，节点之间由“裂隙”相连。将压裂液流动与热流、电流进行类比，节点储存压强信息，流体在节点之间流动时会受到阻力，进而产生压降，该压降由经典的哈根泊肃叶公式得到

(1)

(2)

式中：q为裂隙内的液体流量；ΔP为节点间的压力降；Rh为节点间裂隙流动阻力；μ为液体黏度；L为裂隙长度；Dr为裂隙截面积水力半径；w和d分别为截面的宽和高。将裂隙截面简化为矩形，那么水力半径Dr的表达式为

qij为流进j节点的流量；qjk、qjn为流出j节点的流量图1 四级多尺度复杂裂缝示意图Fig.1 Geometric diagram of multi-scale complex fracture

(3)

利用基尔霍夫定律，每个节点流入和流出的液量之和为0，以节点j为例，可得

(4)

式(4)中：qij为流进j节点的流量；qjm、qjn为流出j节点的流量。

设裂缝入口处为恒定流量边界条件，各个出口处压力为定值，均为原始地层压力，通过联立式(1)～式(4)则可得到压裂液在主裂缝以及各级分支缝中的分配情况。

1.2 模型验证

为了验证当前模型的有效性，采用Fluent对不同复杂裂缝系统液量分配情况进行模拟，并将模拟结果与本文液量分配模型[式(1)～式(4)]预测结果进行对比，进而说明当前模型在预测复杂裂缝液量分配方面的有效性。为了更全面地对模型进行检验，模拟的复杂裂缝系统包括单分支二级裂缝、多分支二级裂缝、多分支四级裂缝，同时为了说明当前模型对于主/次裂缝斜交情况的适用性，第4个算例考虑了主裂缝与分支裂缝呈30°、45°、60°以及120°斜交的情况。

1.2.1 单分支二级裂缝

单分支二级裂缝包括一条主裂缝与一条分支裂缝，主裂缝与分支裂缝夹角为90°，如图2(a)所示。主裂缝缝长、缝宽、缝高分别为20、1、4 cm，分支裂缝位于距主裂缝左侧入口12 cm处，缝长、缝宽、缝高分别为8、1、4 cm，计算域内采用均匀网格，网格尺寸为0.05 cm。主裂缝左侧为入口，主裂缝右侧与分支裂缝为两个出口，压裂液选用清水。模拟开始时，压裂液从主裂缝左侧流入，从主裂缝右侧与分支裂缝流出。入口采用恒定流量边界条件，出口1与出口2为恒定压力边界条件，模拟中考虑了三种不同的入口流量。

达到稳态后，裂缝系统内速度分布如图2(b)所示，并统计了不同入口流量下，出口1与出口2的流量分布。同时为了检验本文液量分配模型有效性，采用式(1)～式(4)计算了相同入口流量下，出口1与出口2的流量，图3为当前预测模型与Fluent模拟结果的对比，两者吻合较好，误差在3%之内。

Fluent进行液量分配模拟时计算量大，涉及的网格数众多，如引言部分所述，页岩裂缝系统具有多尺度特征，主裂缝(100 m)与分支裂缝(100 μm)相差4～5个数量级，若采用Fluent对该系统进行液量分配模拟时，网格数将达到计算机无法承受的量级。而本文提出的液量分配模型采用解析求解的方法，不需要进行网格划分，只需计算不同节点间裂隙流动阻力即可，在保证与Fluent具有相同精度的情况下，计算量与计算时间几乎可以忽略不计。

图2 单分支二级裂缝系统示意图Fig.2 Schematic of single branch secondary fracture system

图3 单分支二级裂缝本文预测结果与Fluent模拟结果对比Fig.3 Comparison between prediction results and fluent simulation results

1.2.2 多分支二级裂缝

准确的液量分配模型是研究复杂裂缝支撑剂运移的关键与基础，文献[4]参照限流压裂法提出复杂裂缝液量分配模型，但与Fluent模拟结果差异较大。主要原因是前人模型在计算流动阻力时只考虑了裂缝截面积大小而忽略了裂缝长度的影响，本文模型综合考虑了裂缝截面积与长度的影响，因此计算精度更高。

图4为本文模型与前人模型预测结果的对比，可以发现对于多分支二级裂缝，本文模型大大提高了分支缝流量的预测精度，与Fluent模拟结果吻合更好，弥补前人模型没有考虑沿程流动阻力的不足。

1.2.3 多分支四级裂缝

上述两个算例只考虑了一级主裂缝与二级分支缝，实际压裂形成的复杂裂缝可能存在多级特征，因此本节第三个算例考虑了具有四级分支缝的复杂裂缝。如图5所示，主裂缝左侧为压裂液入口，其余各级分支裂缝以及主裂缝右侧为出口，主裂缝、二级裂缝、三级裂缝、四级裂缝缝长分别为20、7、5、5 cm，缝宽分别为1、0.8、0.6、0.4 cm，缝高均为4 cm。Fluent模拟得到的各级分支缝流量如图5所示，可见当前模型对于多分支多级裂缝系统也具有较好的预测效果。

图4 多分支二级裂缝系统本文预测结果与Fluent结果对比Fig.4 Comparison between prediction results in this paper and fluent simulation results

图5 多分支四级裂缝系统本文预测结果与Fluent结果对比Fig.5 Comparison between prediction results in this paper and fluent simulation results

1.2.4 分支裂缝与主裂缝斜交系统

前3个算例中主裂缝与分支裂缝夹角均为90°，然而受到最小主应力方向、天然裂缝分布影响，主裂缝与分支裂缝可能会发生斜交，此时两者之间存在一定的夹角。侯腾飞[6]基于参考文献[4]的模型进行了角度修正，在流量分配关系式中引入了角度影响f(θ)，f(θ)的具体表达式通过参数拟合得到。为了明晰分支裂缝与主裂缝夹角对于流量分配影响的物理机制，并验证当前模型在斜交情况下的适用性，本节最后一个算例考虑不同角度的斜交裂缝，包括30°、45°、60°、120° 4种情况，如图6所示。

w为分支裂缝实际缝宽；L为主裂缝出口宽度；θ为分支裂缝与主裂缝夹角图6 单分支二级裂缝本文预测结果与Fluent结果对比Fig.6 Comparison between prediction results in this paper and fluent simulation results

实验中考虑主裂缝与分支缝夹角影响时，通常固定L改变θ，通过这种方法改变裂缝夹角时导致分支缝实际缝宽变窄W=Lsinθ，而前人模型直接将L作为分支缝宽进行计算，因此需要修正。如果将实际缝宽W=Lsinθ代入本文模型计算流动阻力，预测结果与Fluent模拟结果吻合良好，不需要再次进行修正。因此，本文研究澄清了分支裂缝夹角对于液量分配影响的本质是对于分支裂缝实际宽缝的影响，只要计算时将实际缝宽带入预测模型就不需要再对角度进行修正，同时也说明本文模型不仅适用于正交情况，对于主缝与分支缝斜交的情况也同样适用。

2 支撑剂转向输运模型

2.1 理论模型

影响支撑剂在主/次裂缝交汇处输运行为的因素主要有以下几点：主裂缝与分支裂缝几何形态、压裂液物性参数与流动速率、支撑剂物性参数等。本文支撑剂转向输运模型包括两部分：第一部分为在已知主裂缝与分支裂缝压裂液流速的情况下，支撑剂能否发生转向；第二部分为如果可以发生转向输运，有多少支撑剂能够进入下一级分支裂缝。

2.1.1 支撑剂转向条件

支撑剂被压裂液携带进入主裂缝后会随着压裂液一起运移，运移一段距离后发生沉降，形成砂堤，随着砂堤高度逐渐增加，压裂液过流面积减小，流速开始增大，当压裂液流速大于支撑剂颗粒水平临界启动速率时，沉降后的支撑剂再次被压裂液携带运移。

Sahai等[18]实验研究了支撑剂在主/次裂缝交汇处的输运行为，发现主裂缝内支撑剂进入分支缝主要通过两种模式：第一种为砂堤中处于沉降状态的支撑剂依靠重力滚落进入分支缝；第二种为处于悬浮状态的支撑剂跟随压裂液一起进入分支缝。第一种模式容易在裂缝交汇处发生堵塞，进入分支缝的支撑剂大多集中在入口处且数量有限，所以只有通过第二种模式才能将大量支撑剂有效输运到分支缝。因此，支撑剂能够进入分支缝的前提条件是主裂缝内压裂液流速大于支撑剂水平临界启动速率(水平临界携砂流速)，保证支撑剂在主裂缝内处于悬浮状态而不发生沉降。

目前中外学者针对垂直井的临界携砂流速研究较多，而关于水平方向固体颗粒的临界启动速率关注较少。曾思睿等[19]通过实验分析了不同角度对于临界携砂流速的影响，发现水平方向临界启动速率是垂直方向的1.71倍。

垂直方向临界启动速率可以通过以下受力分析得到，固体颗粒在液体中受到的拖曳力FR与浮力Ff为动力，自身重力Fg为阻力，分别表示为

(5)

(6)

(7)

式中：ds为颗粒直径；ρs、ρl分别为固体颗粒密度、流体密度；CD为阻力系数取决于颗粒雷诺数(Rep)；Vs、Vl分别为颗粒速度、流体速度；g为重力加速度。

(8)

式(8)中：

(9)

颗粒受力平衡后，求得

(10)

令Vs=0得到垂直方向临界启动流速Vl，即

(11)

根据参考文献[19]实验结果，得到水平方向临界启动速率，即

(12)

Sahai等[18]通过实验发现，主裂缝内流速越大，支撑剂越容易在裂缝交汇处发生转向，且对于一条主裂缝与一条分支缝构成的裂缝系统，颗粒发生转向输运的临界流速介于0.1～0.23，平均值为0.165 m/s。将实验中采用的支撑剂、压裂液物性参数代入式(11)得到本文模型预测的临界速率位于0.14～0.29 m/s，平均为0.2 m/s，两者吻合较好，说明了当前临界启动模型的有效性。

值得注意的是，主裂缝内压裂液流速大于水平临界流速Vmin[式(12)]只是支撑剂发生转向的必要条件，而非充分条件。当支撑剂大小与分支裂缝缝宽尺寸相当时，支撑剂能否发生转向还必须考虑支撑剂粒径影响。如图7所示，假设主裂缝内压裂液流速大于颗粒A与颗粒B的水平临界启动速率，但显然两个颗粒在转向处的输运行为不同，只有颗粒B可以完成转向，颗粒A由于尺寸较大，无法进入分支缝，继续在主裂缝内运动。

通过大量文献调研发现，Yamada等[20]在研究微流动通道细胞分流时，考虑了不同尺寸细胞在转向处的输运行为，并提出细胞能够发生转向的条件，其核心思想如图7所示。

设主裂缝入口处压裂液流量为Q且速度V沿缝宽度方向服从抛物线分布，w1为主裂缝宽，w2为分支裂缝缝宽，S1、S2为速度沿主裂缝方向积分的面积，则有以下关系：

图7 颗粒转向输运示意图Fig.7 Schematic of particle diversion transportation

(13)

流量Q在裂缝交汇处发生分流，根据流线分布可知，距壁面距离小于H处的压裂液将进入分支缝(图7)流量为Q2，其余压裂液继续在主缝内流动流量为Q1(Q1=Q-Q2)。由于支撑剂跟随压裂液做随体运动，因此只有半径小于H的支撑剂颗粒才可能跟随压裂液Q2一起进入分支缝，否则支撑剂只能在流量Q1的携带下继续在主裂缝内运动。特征长度H主要取决于Q1与Q2的相对大小：

(14)

如图7所示，Q1与Q2可以表示为

(15)

(16)

由1.1节给出的液量分配模型与式(14)～式(15)联立求解则可以得到H。根据当前模型可知，如果分支裂缝较窄时，流动阻力大，进入分支缝的压裂液比例较低(Q2/Q)，因此得到的特征尺度H较小，只有尺寸更小的支撑剂颗粒才能够进入。

综上所述，复杂裂缝系统颗粒转向输运条件为：①主裂缝压裂液流速大于支撑剂颗粒水平临界启动速度；②支撑剂颗粒半径小于特征长度H。

(17)

(18)

2.1.2 砂量分配模型

当支撑剂满足转向输运条件后，另一个重要问题是有多少支撑剂进入了分支缝，即砂量如何在主裂缝与分支缝之间进行分配[21-22]。假设处于悬浮状态的支撑剂在压裂液中均匀分布，并在惯性力作用下跟随压裂液一起做随体运动，此时砂量在裂缝分叉处的分配等同于液量分配。以二级分支裂缝为例(图7)，由液量分配模型[式(1)～式(3)]可得分支裂缝与主裂缝砂量比为

(19)

2.2 模型验证

针对支撑剂临界转向速率模型已在2.1节进行了验证，本节主要对砂量分配准则[式(19)]进行验证，目前定量的实验结果较少，固采用数值方法对支撑剂在复杂裂缝内的输运行为进行模拟，并通过统计主裂缝与分支缝内支撑剂的个数确定砂量比，并将模拟结果与砂量分配模型[式(19)]进行对比，以验证模型的正确性。

采用计算流体力学(CFD)与离散单元法(DEM)耦合模拟支撑剂在复杂裂缝内的输运过程，该耦合算法中，采用CFD计算复杂裂缝内流场信息以及支撑剂受到的流体作用力，采用DEM计算支撑剂在压裂液中的输运行为。该方法相较于传统的双流体模型求解精度更高，可以更准确地描述支撑剂颗粒在裂缝交叉处的输运行为，但缺点是计算时间较长，因此，主要考虑一个主裂缝与一个分支裂缝构成的二级裂缝系统，如图8所示共有5组不同的缝长、缝宽、夹角组合。

模拟结果如图9所示，支撑剂通过裂缝交汇处后，一部分进入了分支缝，另一部分仍在主裂缝内运移，通过统计进入分支裂缝以及留在主裂缝内的支撑剂数目，可以得到不同二级裂缝系统对应的砂量分配情况。图9为进入分支裂缝的支撑剂数目随时间的变化，可以看出在0.32 s时，支撑剂开始进入分支裂缝，并随着时间的增加，进入的支撑剂数量逐渐增加，其斜率代表单位时间进入分支裂缝内的支撑剂数量。同样可以算出残留在主裂缝内的支撑剂数量，进而得到分支裂缝与主裂缝之间的砂量分配情况。

图8 5种不通过缝长、缝宽、夹角的二级裂缝系统Fig.8 Five fracture systems with different fracture lengths， widths and angles

图9 分支裂缝中支撑剂累计个数随时间变化曲线Fig.9 Variation curve of cumulative number of proppant in branch fracture with time

图10显示了采用CFD耦合DEM模拟得到的五种二级裂缝系统(不同缝宽、缝长、夹角)进入分支裂缝与残留在主裂缝内的支撑剂比值，同时采用式(19)预测了不同情况对应的砂量比，结果显示两者吻合较好，最大误差小于10%。

3 结果与分析

前文介绍了本文提出的液量分配模型与支撑剂转向输运条件，本节主要应用该模型定量分析天然裂缝开启程度对于压裂液流量分配以及支撑剂输运情况的影响。考虑三级分支裂缝系统(图11)，其中一级裂缝(缝长100 m、缝宽10 mm、缝高10 m)与二级裂缝(缝长30 m、缝宽 5 mm、缝高10 m)为水力压裂形成的人工缝，三级裂缝为人工缝沟通的天然裂缝(缝高10 m、缝长10 m)。压裂液与支撑剂从入口处注入，出口1、2、3为恒定压力边界条件。

如图12所示，考虑5种天然裂缝开启程度(1、2、3、4、5 mm)，随着天然裂缝宽度的增加，进入天然裂缝的压裂液流量(以入口流量进行归一化)逐渐加大。当天然裂缝宽度为1 mm时，进入天然裂缝流量仅为0.001，当宽度增加为5 mm时，流量增加为0.064，提高了64倍。

图10 不同结构对应的进入分支裂缝与残留在主裂缝内的支撑剂比值Fig.10 Ratio of proppant number of entering branch crack and remaining in main crack corresponding to different structures

如前文所述，压裂液注入排量越大，支撑剂越容易发生转向，也更容易被压裂液携带进入下一级裂缝。利用本文提出的颗粒转向模型，可以预测不同排量下，能够进入图11所示的二级裂缝(缝宽5 mm)与天然裂缝(缝宽2.5 mm)对应的最大支撑剂粒径(表1)。

如表1显示，虽然二级裂缝缝宽5 mm，但只有直径小于1.6 mm的支撑剂方可进入，对于缝宽2.5 mm的天然裂缝，随着排量的增加，可进入天然裂缝的支撑剂直径由0.4 mm增加到1.5 mm。因此为了有效提高复杂裂缝的支撑效果，必须增加压裂液注入排量，以保证支撑剂顺利进入开度较小的次级分支裂缝。

图11 三级裂缝系统示意图Fig.11 Schematic diagram of tertiary fracture system

图12 不同天然裂缝缝宽对应的归一化流量Fig.12 Normalized flow corresponding to different natural fracture widths

表1 不同排量进入二级裂缝、天然裂缝最大支撑剂直径Table 1 Maximum proppant can enter the branch and natural cracks with different velocity

4 结论

提出了一种多级裂缝系统支撑剂转向输运条件以及液量、砂量分配计算模型，该模型物理图像清晰、计算简单、工程可用，可以有效预测支撑剂在复杂裂缝中的分布情况，同时与前人模型相比大幅度提高了预测精度。本文模型采用解析法进行预测，不需要进行网格划分以及大规模数值模拟，因此克服了Fluent等商业软件由于计算量大，无法考虑实际复杂裂缝系统多级、多尺度特性的不足。同时澄清了分支裂缝夹角对于液量分配的影响本质是对于分支裂缝缝宽的影响。最后应用该模型定量研究了天然裂缝开启程度对于压裂液分配情况以及支撑剂输运效果的影响，有望为现场压裂设计支撑剂粒径、压裂液排量选取提供理论依据。