王雪飞，王素玲，侯 峰，王 明，李雪梅，孙丹丹

(1.东北石油大学，黑龙江 大庆 163318；2.中国石油新疆油田分公司，新疆 克拉玛依 834002；3.中科院高能物理研究所，北京 100000；4.散裂中子源科学中心，广东 东莞 523000；5.齐齐哈尔大学，黑龙江 齐齐哈尔 161006；6.中国石油大庆油田有限责任公司，黑龙江 大庆 163514)

0 引 言

在压裂中，压力梯度驱动的流动以及支撑剂和周围携砂流体的密度差导致的对流，二者的共同作用是支撑剂在裂缝中运移的主要原因[1-8]。国内外学者通过实验和数值模拟方法，对支撑剂运移规律进行了大量研究。侯磊等[9]将复杂裂缝系统抽象成单元化物理模型，开展了分流量和转向条件研究。乔继彤等[10]考虑支撑剂沉降因素，建立了支撑剂二维输送方程以及含砂液的二维运动方程，应用有限元方法模拟了水力压裂中的支撑剂输送过程。Kern等[11]基于室内实验方法，对砂浆混合物中砂粒运移机理进行了研究。Shokir等[12]采用小玻璃模型模拟水力裂缝，研究注入速率、支撑剂用量和聚合物类型等参数对支撑剂运移的影响。Chun等[13]研究了裂缝倾角、支撑剂粒径和注入点位置对支撑剂运移的影响。张涛等[14]建立Euler-Euler流固耦合双流体模型，模拟单缝内水输送支撑剂的行为。徐加祥等[15]利用COMSOL软件中的自由与多孔介质流动耦合模型，对压裂液在微裂缝中的流动进行模拟。Zeng等[16]采用代表性颗粒模型(RPM)对支撑剂颗粒进行放大，研究支撑剂在裂缝中的运移机理。Siddhamshetty等[17-18]采用三维多相颗粒单元(MP-PIC)模型，利用欧拉-拉格朗日方法，模拟多尺寸支撑剂的运移。Rivas等[19]利用有限体积法建立支撑剂通过水力裂缝运移的数值模型，模拟KGD裂缝中支撑剂的运移。Kou等[20]使用CFD-DEM方法研究支撑剂在原生裂缝和次生裂缝中运移。Zeng等[21]基于浸没边界法和解析CFD-DEM方法模拟水-颗粒耦合机制的支撑剂运移过程。Wu等[22]利用CFD-DEM方法，模拟水平井射孔套管中支撑剂的运移过程。Wang等[23]采用CFD-DEM方法研究影响支撑剂在裂缝及裂缝网络中分布规律的因素，研究裂缝闭合程度、注入速度和入口开度对单个裂缝中支撑剂分布的影响。

上述对支撑剂运移和铺置的研究中均将裂缝简化为平直裂缝，忽略了迂曲裂缝对支撑剂运移和铺置的影响。现场微地震监测和实验室模拟压裂裂缝均显示压裂裂缝具有迂曲特点[24-27]。因此，研究支撑剂铺置规律时考虑裂缝迂曲情况十分必要。基于CFD-DEM方法，研究了裂缝迂曲度、支撑剂注入速度、压裂液砂比、支撑剂粒径对迂曲裂缝中支撑剂运移和铺置的影响，为现场压裂作业提供理论依据。

1 颗粒与流体耦合计算方法

目前颗粒与流体的耦合计算方法主要有欧拉-欧拉法[28-29]和欧拉-拉格朗日法[30]。此次研究采用欧拉-拉格朗日法，流体采用Navier-Stokes方程求解，而颗粒采用离散元方法，利用牛顿运动定律预测颗粒运动轨迹，并判断颗粒相互碰撞情况。欧拉-拉格朗日法结合了流体相的连续描述和基于牛顿运动方程的颗粒相拉格朗日表征，这种模拟两相流的组合方法可提供两相耦合微观行为的详细信息。

建立考虑颗粒-颗粒、颗粒-壁面和流体-颗粒相互作用的固液两相流耦合模型。对模型做以下假设：流体是不可压缩的，流体流变性和温度在模拟过程中不发生变化，流体与壁面间不滑移，支撑剂是规则的球形颗粒。由于欧拉-拉格朗日模型计算耗时较长，因此，采用实验室规模的裂缝模型。

1.1 流体控制方程

流体质量守恒方程：

(1)

流体动量守恒方程：

(2)

两相流与单向流的区别在于颗粒在计算单元内占据的流体体积。改变流体体积分数 将影响流体运动状态，流体运动状态影响颗粒受到的流体作用力，从而改变流体与颗粒动量交换力源项 。

(3)

(4)

1.2 颗粒受力及动力学方程

颗粒-流体的流动系统中，颗粒进行的平移和旋转运动符合牛顿第二定律[31]，动力学方程为：

(5)

(6)

颗粒受到的流体作用力的方程为：

Ff=Fd+Fvm+Fpg+Fsa+Fb+Fm

(7)

式中：Fd为颗粒受到的流体阻力(曳力)，N；Fvm为虚拟质量力，N；Fpg为压力梯度力，N；Fsa为Saffman升力，N；Fb为Basset力，N；Fm为Magnus力，N。

已有研究[32]认为，作用于颗粒上的力主要是流体阻力，其他力比较小可忽略不计。因此，颗粒受到流体作用力只考虑流体阻力。

当黏性流体绕流颗粒时，流体对颗粒产生与流动方向一致的作用力，该作用力与流体和颗粒相对速度(uf-up)有关。当流体的流速小于颗粒的运动速度时，作用于颗粒上的力表现为阻力；当流体的流速大于颗粒运动速度时，作用在颗粒上的力表现为曳力，计算公式为：

(8)

(9)

(10)

式中：Cd为曳力系数；Rep为颗粒雷诺数；Dp为颗粒直径，m；μf为流体表观黏度，Pa·s。

1.3 颗粒碰撞模型

如图1所示，在外力作用下颗粒i与颗粒j在c点碰撞，虚线表示开始碰撞时颗粒j的位置，随着两颗粒相对运动，产生法向重叠量γ和切向位移δ(单位均为m)。通过引入弹簧、阻尼器、滑动器和耦合器等计算颗粒间产生碰撞力，如图2所示。

图1 两颗粒碰撞示意图Fig.1 The schematic diagram of collision between two particles

图2 弹簧-阻尼系统接触力模型Fig.2 The contact force model of spring-damping system

对于三维坐标中的球体颗粒，Fσij为作用于颗粒i的碰撞力的法向分量，可等效为弹簧和阻尼器产生的力之和：

(11)

式中：uij为颗粒i与颗粒j的相对速度矢量，uij=ui-uj；a为颗粒i的中心至颗粒j的中心直线方向上的单位矢量；kσ、cσ分别为颗粒i的法向弹性系数和法向阻尼系数。

碰撞力的切向分量Fτij为：

Fτij=-kτδ-cτucτ

(12)

ucτ=uij-(uij·a)a+riωi×a+rjωj×a

(13)

式中：kτ、cτ分别为切向弹性系数和切向阻尼系数；ucτ为接触点的滑移速度，m/s；ri、rj分别为颗粒i和颗粒j的半径，mm；ωi、ωj分别为颗粒i和颗粒j的角速度，rad/s；

综上所述，作用在颗粒i上的合力和合力矩分别为：

(14)

式中：Fij为作用在颗粒i上的合力，N；Tij为作用在颗粒i上的合力矩，N·m；。

当多个颗粒同时与颗粒i接触时，作用在颗粒i上的总力和总力矩为：

(15)

式中：Fi为作用在颗粒i上的总力，N；Ti为作用在颗粒i上的总力矩，N·m；。

1.4 CFD-DEM耦合计算方法

CFD-DEM耦合求解流程：首先建立连续相与离散相数值模型，设置边界条件，计算流体体积分数αf和动量交换力源项Sf，启动CFD求解器，计算首次迭代的动量离散方程中的系数和常数项，求解动量离散方程，对速度和压力进行修正，并判断是否满足收敛准则。若不满足，则重新求解动量离散方程的系数和常数，再次修正速度和压力，直至满足收敛准则；若满足收敛准则，流体域推进至下一个计算时间步。启动DEM求解器，根据插值函数获得颗粒中心位置流体速度uf和颗粒受到的流体作用力Ff，使用质点坐标方法追踪颗粒并更新颗粒位置和颗粒运动速度up，判断颗粒碰撞情况并计算颗粒碰撞力，进行颗粒动力学方程求解，更新颗粒运动速度up，再次计算体积分数αf和动量交换力源项Sf，开始下一次流体求解，直至耦合求解结束。

2 边界条件及验证

2.1 边界条件建立

在压裂过程中，先将前置液注入地层，打开裂缝，然后注入均匀混合的压裂液。为了准确模拟这一过程，先在CFD中计算稳态流场直至收敛，将稳态流场作为初始条件，建立裂缝物理模型，对流体域进行六面体网格划分，并对边界层进行网格加密处理。采用SIMPLEC方法和瞬态解算器求解流体控制方程。设置速度入口边界和压力出口边界，颗粒与流体初始速度相同。

颗粒选用石英砂材质，密度为2 650 kg/m3，弹性模量为5.5×1010Pa，泊松比为0.13。裂缝边界为岩石，密度为3 000 kg/m3，弹性模量为6.5×1010Pa，泊松比为0.20。颗粒采用Hertz-Mindlin接触模型，颗粒接触参数采用虚拟实验进行标定，恢复系数为0.40，静摩擦系数为0.68，动摩擦系数为0.01。

2.2 数值模型验证

基于建立的CFD-DEM耦合数值模型，对文献[33]中“一”型裂缝砂丘铺置形态实验过程进行数值模拟，对比结果如图3所示。由图3可知：“一”型裂缝中砂丘铺置高度和铺置的距离、铺置形态基本一致，验证了建立的数值模型的准确性。

图3 “一”型裂缝中支撑剂铺置砂丘形态对比Fig.3 The comparison of proppant displacement shape in linear fracture

3 迂曲裂缝中支撑剂运移影响因素分析

裂缝迂曲度是裂缝中心线长度与裂缝首尾直线距离的比值。数值模拟中，流体和支撑剂从同一入口流入，初始速度相同。

3.1 迂曲度对支撑剂运移的影响

模拟条件：裂缝迂曲度为1.0、1.2、1.5、2.0，如图4所示；流体注入速度为1.0 m/s，模拟时间为0.5 s。压裂液在不同迂曲裂缝中的流速如图5所示，流速损失率随进入裂缝距离变化如图6所示，支撑剂在裂缝内分布如图7所示。

图4 裂缝迂曲度示意图Fig.4 The schematic diagram of crack tortuosity

图5 裂缝中流速分布Fig.5 The flow rate distribution in fracture

图6 流速损失率Fig.6 The loss ratio of flow rate

图7 迂曲度对裂缝中支撑剂分布的影响Fig.7 The effect of tortuosity on proppant distribution in fractures

由图5可知：压裂液沿裂缝流动时，流速随着进入裂缝的距离增大而减小。其原因为压裂液具有一定的黏度，在流动过程中与裂缝壁面间的摩阻力和裂缝形状引起的局部阻力造成流体能量损失。迂曲度为1.0的裂缝内流速呈线性下降，迂曲度为1.2、1.5、2.0的裂缝内流速呈非线性下降。由图6可知，随着流体进入裂缝距离的增大，流体流速损失率增大。迂曲度越大流速损失越快，进入裂缝20 mm处，与迂曲度为1.0的裂缝中的流速损失率相比，迂曲度为1.2、1.5、2.0的裂缝的流速损失率分别增大2.2、3.7、4.5倍。其原因为流体进入裂缝距离相同时，迂曲度越大实际流过距离越长，流体在迂曲裂缝内还会形成涡旋，流体产生的沿程和局部压力损失增大。进入裂缝初期流速损失率增长最快，随着进入裂缝距离的增大，流速损失率逐渐减缓。

由图7a可知：迂曲度为1.0时，支撑剂向前平稳运移，形成的砂丘坡度较缓，铺置距离较长，支撑裂缝面积大。砂丘底层颗粒受上部颗粒的压实作用和周围颗粒的阻力停止运动，砂丘中部颗粒受周围颗粒的阻力和压裂液的携带，以较小速度向前滚动，大部分上层颗粒以较大速度和流体一起向前运移。

由图7b、c、d可知：迂曲度由1.2增至2.0时，迂曲裂缝阻力较大，流体流动状态复杂，流体流动速度不均匀；颗粒与壁面以及颗粒与颗粒之间碰撞加剧，导致流体和颗粒剧烈减速，颗粒在裂缝内发生沉降，形成堵塞，造成大面积支撑剂沉降，只有少部分上层颗粒以较小速度向前滚动；砂丘铺置距离较短，坡度较大，颗粒在裂缝入口堆积，无法有效支撑裂缝。

3.2 注入速度对支撑剂运移的影响

模拟条件：支撑剂注入速度为0.1、0.5、1.0、1.5 m/s，裂缝迂曲度为1.0、1.2、1.5、2.0，流体和支撑剂的初始速度相同，压裂液砂比为10%。支撑剂在裂缝中分布如图8所示(以迂曲度为1.2的裂缝为例，下同)；各迂曲度裂缝中注入速度对铺置距离的影响如图9所示。

图8 注入速度对裂缝中支撑剂分布的影响Fig.8 The effect of injection rate on proppant distribution in fractures

图9 注入速度对铺置距离的影响Fig.9 The analysis of the effect of injection rate on displacement distance

由图8a可知：注入速度为0.1 m/s时，砂丘堆积在裂缝入口处，继续注入支撑剂砂丘高度增长迅速，最终堵塞裂缝入口，后续支撑剂将无法进入裂缝，而砂丘长度增长甚微，说明后续支撑剂起不到支撑裂缝作用。

由图8b可知：注入速度为0.5 m/s时，迂曲裂缝内的颗粒与颗粒、颗粒与壁面剧烈碰撞，流体局部压力损失增大，大部分支撑剂在裂缝迂曲处沉降堆积成砂丘，继续注入支撑剂，在裂缝迂曲狭窄处继续堆积形成堵塞，后续注入支撑剂不能进入裂缝深处，导致裂缝不能得到有效支撑。

由图8c可知：注入速度为1.0 m/s时，支撑剂在迂曲裂缝中减速明显，但仍以一定速度边向前运移边沉降形成砂丘，相比注入速度低的情况，高注入速度形成的砂丘形状较平缓，砂丘铺置面积增大，继续注入支撑剂能够延长裂缝支撑长度。

由图8d可知，注入速度为1.5 m/s时，支撑剂克服裂缝迂曲带来的阻力，在裂缝中铺置距离更长，更均匀。

由图9可知：支撑剂的铺置距离随注入速度的增加而增大，但不同迂曲度的裂缝中支撑剂对注入速度的敏感度不同，注入速度从0.1 m/s增至1.5 m/s，迂曲度为1.0、1.2、1.5、2.0的裂缝中支撑剂铺置距离分别延长1.85、1.72、1.52、1.30倍。

迂曲度为1.0的裂缝中，支撑剂运移时所受到阻力较小，随着注入速度增加铺置距离增长较快。随着裂缝迂曲度的增大，支撑剂向前运移阻力越大，受力复杂，随着注入速度增加铺置距离增长变慢。当注入速度小于1.0 m/s时，支撑剂和流体的动能较小，支撑剂在裂缝内运移受阻，发生快速沉降，铺置距离较短。而注入速度大于1.0 m/s时，支撑剂和流体具有足够的动能抵抗迂曲裂缝的形状阻力，使得支撑剂越过迂曲地带后仍然具有一定动能，运移一段距离才发生沉降。因此，迂曲度越大的裂缝，需要更大的注入速度才能获得较理想的铺置效果，注入速度较小时，支撑剂会阻塞裂缝入口，使后续支撑剂无法运移至裂缝深处，支撑裂缝效果较差。

3.3 压裂液砂比对支撑剂运移的影响

模拟条件：压裂液砂比为10%、20%、30%、40%，注入速度为1.0 m/s，支撑剂粒径为0.6 mm。支撑剂在裂缝中分布如图10所示。

由图10可知：压裂液砂比由10%增至40%时，支撑剂铺置距离缩短。其原因为裂缝中注入高砂比压裂液时，颗粒与颗粒、颗粒与壁面碰撞增多，颗粒动能损失增大，颗粒发生快速沉降，在裂缝入口附近发生堵塞，阻碍后续颗粒的进入，无法获得长距离支撑效果。而低砂比压裂液运移时，颗粒碰撞次数减少，颗粒能够跟随压裂液进入裂缝深处，铺置距离增大。由此可见，砂比越小支撑剂运移过程阻力损失越小，铺置越均匀，但减小砂比会增加泵注时间。

图10 压裂液砂比对裂缝中支撑剂分布的影响Fig.10 The effect of fracturing fluid-proppant ratio on proppant distribution in fractures

3.4 支撑剂粒径对支撑剂运移的影响

模拟条件：支撑剂粒径为0.2、0.4、0.6、0.8 mm，注入速度为1.0 m/s，压裂液砂比为10%。撑剂在裂缝中分布如图11所示。

由图11可知：支撑剂粒径越大，压裂液对支撑剂的携带性越弱，在迂曲裂缝内通过性越差，与壁面发生多次碰撞后发生快速沉降，支撑剂铺置距离短，且不均匀。支撑剂粒径越小，压裂液对其携带能力越强，支撑剂通过迂曲裂缝数量增多，能够形成铺置距离长且均匀的裂缝。

图11 支撑剂粒径比对裂缝中支撑剂分布的影响Fig.11 The effect of proppant particle size ratio on proppant distribution in fractures

4 结 论

(1) 压裂液在裂缝中流动，入口流速相同条件下，随着裂缝迂曲度的增大，压裂液速度损失越大，进入裂缝初期流速衰减越快，随着进入裂缝距离的增大，流速衰减变缓。裂缝迂曲度越小，支撑剂铺置越均匀，铺置距离越长，迂曲度过大会导致支撑剂在裂缝入口处堵塞，支撑剂运移距离缩短。

(2) 随着支撑剂入口速度的增大，不同迂曲裂缝中支撑剂运移距离均有所增加，迂曲度较大的裂缝，支撑剂运移距离随着速度的增加缓慢增加，迂曲度小的裂缝支撑剂运移距离随着速度增加明显增大。支撑剂注入速度大于1.0 m/s时，不同迂曲裂缝内支撑剂均能铺置一定距离，对裂缝形成有效支撑，速度越大支撑剂铺置距离越长分布越均匀，对裂缝的支撑效果越好。

(3) 支撑剂在迂曲裂缝中运移距离随压裂液砂比的增加而减小。支撑剂粒径越大，流体对其携带能力越差，支撑剂铺置距离越短，且铺置不均匀。