0 引言

大规模水力压裂是实现非常规油气藏商业开发的主要手段。然而大规模水力压裂会带来储层和地下水污染及水资源浪费等问题^[1-3]。针对水力压裂存在的局限性，中国工程院沈忠厚院士率先提出利用超临界CO₂进行储层压裂改造^[4-6]。超临界CO₂压裂不仅对储层和地下水无污染，而且能够使储层产生多而复杂的裂缝网络，大大提高油气采收率，同时还可以实现CO₂埋存。然而超临界CO₂压裂液与其他低黏压裂液相比，黏度更低，因此其在裂缝内支撑剂携带能力和规律不得而知。在超临界CO₂携带支撑剂能力研究方面，L.HOU等^[7]利用试验和模型推导的方法，建立了支撑剂在超临界CO₂中的沉降末速计算模型，能够预测支撑剂在裂缝中沉降特性。郭斌^[8]通过数值模拟证实了超临界CO₂在连续管水平段和裂缝中携带支撑剂的可行性。侯磊等^[9]在BBO方程的基础上，建立了支撑剂在超临界CO₂中的跟随性计算模型，并通过试验验证了模型的准确性。

此前的研究主要关注超临界CO₂在裂缝中携带支撑剂的可行性，而对超临界CO₂在裂缝中携带支撑剂填充和铺设规律研究较少；此外，上述研究还忽略了超临界CO₂压裂液黏度变化的影响，且所取砂比不符合超临界CO₂压裂的特点。鉴于此，笔者采用计算流体力学(CFD)的方法，考虑支撑剂运移中流固和固固间耦合作用的影响，建立了求解裂缝中超临界CO₂和支撑剂两相流动模型，并进行了关键参数敏感性分析，进而得到超临界CO₂在裂缝中携带支撑剂的最优参数，可以为超临界CO₂压裂参数优化及技术适应性研究提供借鉴和参考。

1 基本模型和模拟方案

基于Eulerian-Eulerian的多相流模型主要包含VOF(Volume of Fluid)模型、Mixture(混合)模型和Eulerian(欧拉)模型。其中，VOF模型主要用来求解分层流和自由表面流，欧拉模型求解多相流，常将动量方程和连续性方程分开求解，计算量较大，而混合模型是一种简化的欧拉模型，适用于支撑剂运移和沉降问题，并且混合模型计算量较小，精度高。故笔者采用Mixture模型求解超临界CO₂和支撑剂两相流动。为了考虑颗粒之间的相互作用，采用了颗粒动力学中固相的质量、动量及颗粒温度的输运方程^[10-11]。

1.1 控制方程

式(1)~式(3)分别给出了携砂液的连续性方程、动量守恒方程以及能量守恒方程。

连续性方程为：

(1)

式中：ρ_m为混合物密度，kg/m³；为混合物的质量平均速度，m/s; ▽为算子符号；t为时间，s。

动量方程为：

(2)

式中：ρ_k为k相的密度，kg/m³；p为压力，Pa；为重力加速度， =9.81 m/s²；α_k为k相的体积分数；为k相颗粒的漂移速度，m/s；μ_m为混合物的黏度，Pa·s； →为体力，Pa。

能量守恒方程为：

(3)

式中：E_k为k相的能量，J；为k相的速度，m/s；k_eff为有效的传率，W/(m·K)；T表示温度，K。

根据颗粒流动模型能量理论可知，固相颗粒温度正比于颗粒随机脉动动能，其输运方程为：

(4)

式中：Θ_p为颗粒拟温度，m²/s²；κ_p为颗粒能量传导系数，kg/(m·s)；τ_p为颗粒相的应变张量；为颗粒相的速度张量；J_p为非弹性碰撞导致的能量耗散，kg/(m·s³)；Π_Θ为因流体黏性阻尼而产生的能量耗散，kg/(m·s³)。

1.2 本构方程

求解上述控制方程需要添加本构方程加以封闭，其中对于固相与液相之间的拖曳力计算，在裂缝内往往采用Gidaspow模型^[12-14]。固液两相间的动量交换系数β计算式如下。

当α₁≥0.8时，有：

(5)

当α₁＜0.8时，有：

(6)

其中：

(7)

(8)

式中：α₁为液相的体积分数；v_l为液相速度，m/s；ρ_l为液相的密度，kg/m³；d_p为固相颗粒粒径，m；C_D为固液两相间曳力系数；μ_l为液相的黏度，Pa·s；Re_p为固液两相间滑脱速度定义的雷诺数。

在混合物体系中，颗粒黏度对混合物有效黏度的计算至关重要，因此在体积加权平均黏度的计算中需要考虑颗粒动能引起的剪切黏度。固体剪切黏度μ_s由碰撞黏度部分μ_{s, col}和动能黏度部分μ_{s, kin}组成。

(9)

(10)

(11)

(12)

式中：e为弹性恢复系数；g₀为径向分布函数，是当固体相稠密时用于修正颗粒之间碰撞概率的校正系数；s为颗粒之间的距离，m。

1.3 超临界CO₂状态方程

CO₂的密度和黏度是影响其携带支撑剂性能的因素，考虑CO₂黏度和密度受温度影响，利用Fluent中的Aungier-Redlich-Kwong真实气体模型及Vesovic模型分别计算CO₂的密度和黏度^[15-16]。

(13)

(14)

(15)

(16)

(17)

式中：R为气体常数，J/(mol·K)，k；T_r为对比温度，k；T_c为临界温度，K；p_c为临界压力，Pa；V_c为临界体积，m³/kg；V为对比体积，m³/kg；ω为偏差因子。

1.4 几何模型和边界条件

考虑超临界CO₂压裂裂缝特点，确定采用与文献模型相似的尺寸，模拟裂缝的长、宽、高分别为3 000、10和400 mm，如图 1所示^[17-21]。为了提高计算的准确性和收敛速度，网格划分采用结构化划分方式，且网格节点数72万个。入口的边界条件为质量流量入口，出口边界条件为压力出口。模型的其他部分边界条件为壁面，且假设壁面为无滑脱壁面，温度恒定。考虑到非常规油藏渗透率较低，故忽略超临界CO₂的滤失效应。使用可实现k-ε(Realizable k-ε)模型描述支撑剂在裂缝内流动，同时为了提高计算的准确性，设置入口处的湍流强度和水力当量直径。此外，使用SIMPLE算法求解压力和速度耦合场，压力采用标准的离散格式，其他均采用一阶迎风的离散格式。

1.5 模拟方案

为了探究支撑剂密度、粒径、携砂液的排量以及砂比等参数对超临界CO₂在裂缝中携带支撑剂流动规律的影响，在参考滑溜水等低黏压裂液参数的基础上^[22-23]，制定了数值模拟方案(见表 1)。其中基准的模拟条件为支撑剂密度1 540 kg/m³，粒径为0.25 mm，砂比为0.3，携砂液排量为2 kg/s。

数值模拟分4组进行，对每组模拟结果分析时，均取缝宽方向位于中心平面上的数据。

2 结果与讨论 2.1 超临界CO₂携带支撑剂在裂缝中分布规律

图 2为基准条件下不同时刻裂缝垂直剖面上支撑剂体积分数分布云图。由图 2a可知，当超临界CO₂和支撑剂注入后，支撑剂由于重力作用在裂缝底部沉降形成砂丘；同时由于支撑剂密度较小，加上射流作用，其易被超临界CO₂携带到深部裂缝，而在靠近入口位置形成凹陷。图 2b显示，随着注入时间持续，裂缝中的砂丘会影响携砂液射流速度，导致支撑剂因速度减小发生沉降，从而在距入口一定距离处形成堤峰。图 2c表明，由于堤峰的存在，致使裂缝中过流截面减小，流经裂缝中的支撑剂运移速度会加快，促使支撑剂向砂丘的背部方向堆积和运移，此时沙丘的高度不发生变化，长度不断向出口延伸。上述即是超临界CO₂携带支撑剂在裂缝中沉降和运移的3个阶段，与文献[18]介绍的滑溜水携带支撑剂在裂缝中运移和铺设过程大致相同，说明超临界CO₂同滑溜水具有相似的携带支撑剂的性能。

由图 2可知，支撑剂体积分数在裂缝垂直剖面上具有明显的差异，且由下到上依次可以分为4个区域：静止砂丘层(图中1)、接触层(图中2)、跃移层(图中3)和悬浮层(图中4)^[17-18]。图 3为裂缝半长处支撑剂体积分数与无因次缝高关系。由图可知，在静止砂丘层支撑剂的体积分数维持在大于0.60的值保持不变，接触层体积分数在介于0.25~0.60之间缓慢变化，跃移层支撑剂体积分数在介于0.20~0.25范围内波动剧烈，而悬浮层的支撑剂体积分数则在介于0.0~0.2之间缓慢变化。此外，在10~30 s之间，随着时间的增加，静止砂丘层、接触层以及悬浮层在所占的面积比例不断增大，跃移层所占面积比例不断减小。而在40 s之后，裂缝中各层所占面积比例基本保持不变，即达到了平衡状态。

为了获得支撑剂在裂缝内铺设和填充规律，选取裂缝半长处的砂丘高度和支撑剂体积分数进行分析，并利用无因次平衡高度(平衡高度/裂缝高度)和无因次砂丘高度(当前砂丘堆积高度/裂缝高度)对影响超临界CO₂携带支撑剂性能的参数进行优选。

2.2 支撑剂密度的影响

图 4为无因次平衡高度和平衡时间随支撑剂密度变化的曲线。由图可知，随着支撑剂密度的增大，砂丘的无因次平衡高度逐渐增加，而平衡时间逐渐减小。其中支撑剂密度为1 540 kg/m³时，在裂缝中支撑剂形成的砂丘高度最小，约为裂缝高度的40%，平衡时间是50 s。支撑剂密度为2 540 kg/m³时，支撑剂形成的砂丘高度最大，约为裂缝高度的65%，平衡时间缩短为30 s。

图 5为不同支撑剂密度下砂丘无因次高度随时间变化曲线。由图可知，不同密度的支撑剂所形成的砂丘高度均是先迅速增加，然后逐渐趋向稳定。这是因为当支撑剂和超临界CO₂受到射流加速作用时，可以在射流的滞止点附近较快沉降。然而随着砂丘高度增加，携砂液的射流速度会受到影响，导致支撑剂沉降速度减慢。此外，由于携砂液的过流截面变窄而速度增大，支撑剂易被流化运移到砂丘背面。当支撑剂颗粒沉降数目和砂丘中支撑剂颗粒被流化数目相等时，则达到了平衡状态，砂丘高度不发生变化。

图 5还表明，当支撑剂密度为2 540 kg/m³时，砂丘的高度随时间变化曲线在0~30 s之间的斜率最大，砂丘高度增加较快；而支撑剂密度为1 540 kg/m³时，曲线斜率变化相对较小，砂丘高度变化较慢。因此采用低密度支撑剂有利于其在裂缝中形成更加均匀的砂丘剖面，避免发生砂堵。

2.3 支撑剂粒径的影响

图 6为无因次平衡高度和平衡时间随着支撑剂粒径变化曲线。由图可知，随着支撑剂粒径的增加，砂丘的无因次平衡高度不断增加，而平衡时间逐渐缩短。这是因为支撑剂的粒径越大，支撑剂颗粒越易受到裂缝壁面的影响，故支撑剂易于发生沉降。支撑剂的粒径为0.425 mm时，砂丘的无因次高度最大，约为裂缝高度的52.5%，平衡的时间约为30 s。而支撑剂粒径为0.250 mm时，砂丘的无因次高度最小，约为裂缝高度的40%；平衡的时间接近50 s。

图 7是在不同的支撑剂粒径条件下，砂丘的无因次高度随时间变化曲线。由图可知，不同的支撑剂粒径，无因次砂丘高度随时间变化的趋势均是先迅速增加，之后趋于稳定。其中支撑剂粒径越大，其前期沉降速度越大，无因次砂丘高度变化得越快；并且当支撑剂沉降达到平衡状态以后，支撑剂在裂缝中形成的砂丘高度也较高。因此，支撑剂粒径越大，越容易发生砂堵，不利于支撑剂在深部裂缝形成有效的支撑。故使用超临界CO₂加砂压裂，建议选用高于20/40目数的支撑剂有利于支撑裂缝。

2.4 砂比的影响

图 8为无因次平衡高度和平衡时间随着砂比变化曲线。由图可知，随着砂比的增加，无因次平衡高度逐渐增加，而平衡时间先在低砂比时保持不变之后逐渐缩短。这是因为随着砂比的增加，支撑剂颗粒之间及其与裂缝壁面的作用加剧，促使更易于发生沉降。其中支撑剂的砂比为8%，砂丘的无因次平衡高度最小，约为裂缝高度的10%，平衡时间约为60 s。支撑剂的砂比为3%，砂丘的无因次平衡高度最大，约为裂缝高度的40%，平衡时间接近50 s。

图 9是在不同的砂比条件下，砂丘的无因次高度随时间变化曲线。由图可知，不同的砂比条件下，无因次砂丘高度变化趋势先迅速增加，之后趋于稳定。同时随着砂比的增加，支撑剂颗粒之间的对流沉降作用加强，平衡时间逐渐缩短。此外，在达到平衡状态前，砂比越大，无因次砂丘高度随时间变化曲线斜率越大，支撑剂沉降越快；在达到平衡状态之后，砂比越大，无因次砂丘高度也越大。因此利用低黏的超临界CO₂进行压裂时，适合采用低砂比(小于8%)的携砂液进行裂缝填充。

2.5 质量流量的影响

图 10为无因次平衡高度和平衡时间随着携砂液排量变化曲线。由图可知，随着携砂液排量的增加，无因次平衡高度降低，平衡时间缩短。这是因为入口质量流量越大，单位时间内注入的支撑剂质量越大，形成砂丘所需的时间就会缩短；同时注入的质量流量越大，注入支撑剂的速度较快，故悬浮区和跃移区的面积将变大，而静止砂丘的区域面积会变小。图中携砂液的质量流量为1 kg/s时，无因次平衡高度最大，约为裂缝高度的60%，平衡时间接近70 s；而携砂液的质量流量为4 kg/s时，无因次平衡高度最小，约为裂缝高度的5%，时间接近10 s。

图 11是在不同的携砂液排量条件下，砂丘的无因次高度随时间变化的曲线。由图可知，随着时间的增加，无因次砂丘高度变化趋势是先增加之后趋于稳定。随着单位时间入口排量的增加，单位时间注入的支撑剂质量会增加，则平衡时间会逐渐减小。同时排量越小，达到平衡状态的时间越长，平衡高度越大，形成砂堵的风险越大。因此在实际的超临界CO₂压裂作业中，应根据裂缝特点，选择合理的携砂液排量，避免发生砂堵。

3 结论

(1) 超临界CO₂携带支撑剂在裂缝中铺设时大致分为3个阶段：首先是支撑剂在射流的滞止点附近因重力发生沉降，在裂缝底部形成砂丘；其次砂丘随着支撑剂颗粒沉降不断抬高，在距离入口一定位置处形成凸起；最后，砂丘达到平衡高度后，支撑剂床不断向出口方向延伸。

(2) 使用超临界CO₂携带支撑剂填充支撑裂缝时，支撑剂体积分数在裂缝垂直剖面上差异明显；静止沙丘层的支撑剂体积分数随注入时间增加基本保持不变，接触层和悬浮层体积分数变化较为平缓，而跃移层的体积分数变化剧烈。

(3) 当超临界CO₂最初在裂缝中携带支撑剂流动时，静止砂丘层、接触层以及悬浮层所占的面积比例随着注入时间增加不断增大，而跃移层所占面积比例不断减小。而当达到一定时间后，裂缝中各层所占面积比例基本保持不变。

(4) 支撑剂的密度、粒径、携砂液的砂比以及排量会对超临界CO₂携带支撑剂在裂缝中铺设规律和填充效果产生影响。在模拟条件下，建议选用密度小于1 540 kg/m³、粒径大于20/40目的支撑剂，并且携砂液的砂比要低于8%，排量要高于4 kg/s，从而有益于超临界CO₂携带支撑剂在裂缝中形成均匀的砂丘形态，避免发生砂堵。