超临界CO2压裂缝内支撑剂运移规律

郭兴，孙晓，穆景福，乔红军，罗攀，李恒

（1.陕西延长石油(集团)有限责任公司研究院,西安 710065；2.陕西省二氧化碳封存与提高采收率重点实验室,西安 710065；3.北京大学工程学院,北京 100871）

0 引言

压裂作为一种有效的储层改造增产技术，在油气藏勘探开发过程中应用广泛，但其浪费大量的水资源、对环境及储层造成较大的污染等弊端和局限性也逐渐暴露出来[1–2]。如何在当前“双碳”背景下，“清洁、低碳、环保”的进行油气田勘探开发，保障我国能源安全，是目前亟需解决的问题。

超临界CO2由于其高密度、低黏度、高扩散性等特殊性质[3]，非常适合页岩油气等非常规油资源的勘探开发[4–5]。国内外进行了诸多理论研究和现场试验以及规模化CO2压裂技术现场应用[6–7]，其可以降低破裂压力，实现低渗致密储层造复杂缝网，不仅节省大量水资源，对环境和储层无污染，还能在油气井增产提效的同时实现二氧化碳埋存[8–9]。超临界CO2压裂过程中支撑剂在裂缝中的有效运移和沉降是压裂施工中极其重要的环节之一。

目前对压裂过程中支撑剂运移规律研究主要采用室内实验和数值模拟2种方式[10]。通过实验研究可以很直观地看到支撑架在运移过程中的运动规律，分析得到小尺度（相比于储层裂缝尺度）下描述颗粒群在裂缝中的沉降规律[11]。而对于超临界CO2压裂，在实验室很难通过物理手段模拟其压裂过程中支撑剂在CO2流体流动过程中的运移情况，具有极大的限制性，难以真正反映支撑剂在真实裂缝中的运移情况，如研究的时间和空间尺度与实际生产情况相差太大，边界条件设置与实际情况相差甚远，监测统计流体、颗粒两相速度的手段不够完善等[12]。通过利用数值计算手段则可以弥补实验研究的不足，结合实验研究总结得到的校正模型，能够精确地描述大尺度下支撑剂的运移情况[13]。因此笔者通过数值计算研究不同参数对CO2压裂液携砂效果的影响，为超临界CO2压裂设计优化及现场施工提供重要的技术支撑。

1 数学模型

由于欧拉-欧拉型方法将颗粒相整体视作单一连续相的假设与实际情况相差较大，而欧拉-拉格朗日型粒子类方法中的多相质点网格法（MP-PIC）则是通过颗粒相压力梯度力替代颗粒之间的碰撞，从而表征颗粒间的相互作用力；将具有相同动力学特性的颗粒群作为计算单元，以便达到减小计算量的目的，对于大规模粒子运动模拟计算具有独特的优势[14]。MP-PIC方法中颗粒相的运动不存在因差分格式导致在空间上产生的非物理扩散，因此其运动描述更为准确[15]。考虑到CO2黏度低，裂缝中流体运动复杂，存在湍流效应，以及由于CO2携砂能力差，导致颗粒之间的作用以及颗粒-固壁作用不可忽略，因此采用多相质点网格模型（MPPIC）方法模拟颗粒团在CO2压裂液中的运移沉降过程。

1.1 流体控制方程

基于Patankar and Joseph[16]扩展的二维MPPIC方法，考虑大量支撑剂颗粒在裂缝扩展中的运移，流体相的控制方程可以通过Favre平均法推导得出，鉴于篇幅限制，具体过程不再赘述，具体可参考文献[17]。流体控制方程为体积平均纳维-斯托克斯方程（VANS），包括连续性方程和动量方程[15,18]：

式中，νf为流体运动黏度，其与局部颗粒浓度有关，m2/s。

由于颗粒的存在，携砂液的有效黏度高于纯流体，因此该部分使用Barree和Conway（1994）的黏度模型：

式中，v0是纯流体的运动黏度，m2/s；φcp是紧密堆积颗粒体积分数，%。

1.2 颗粒动力学方程

颗粒群的运动描述采用代表质点法的思想，即一个质点代表多个颗粒，而这些颗粒具有相同的性质如密度、大小等，它们具有相似的运动行为，因此可以看作一个粒子团，运移过程中质量保持不变。质点运动由牛顿第二定律控制，其控制方程[20–21]可表示为：

采用Snider[22]的粒子法向应力模型来解决粒子碰撞和紧密堆积问题，τp可表示为：

式中，ps和 β 为模型参数，Pa和无量纲；φcp是紧密堆积时的颗粒体积分数，%。这里ps和 β分别设置为4 Pa和2[22]。

采用基于BGK近似模型的一个粒子分布函数输运方程来对描述碰撞阻尼机制[21]，为了考虑壁摩擦效应，这里只考虑颗粒紧密堆积区域。为粒子相的质量平均速度，表示为可以表示为：

式中，α为颗粒速度分布的质量加权标准差，无量纲；η为模型系数，无量纲；r32为Sauter平均半径，m，可表示为[15]：

式中，rp是粒子半径，m；在颗粒均匀情况下，Sauter平均半径是颗粒相的平均半径。

流体-颗粒之间存在相互作用，其中曳力最为重要。这里采用Wen-Yu模型[23]，则曳力系数可表示为：

式中，r是颗粒半径，m；为颗粒与流体密度比，无量纲；Cd为无量纲的曳力系数，表示为：

1.3 求解方法

流体部分求解采用SIMPLE压力校正迭代方法进行求解，时间项离散采用欧拉格式，对流项采用显式一阶迎风格式。颗粒部分与颗粒速度相关的项均采用隐式离散，其余项如压力，颗粒应力等都由流场数据显式给出。由于MP-PIC为混合型方法，即流体描述为欧拉型方法，而颗粒描述采用的是拉格朗日型方法，二者之间的数据交换需要通过插值得到，采用二维双线性插值。

数值求解过程中裂缝扩展、颗粒团运动和流体运动3个过程是通过顺序和显示耦合实现的，即假设流体和颗粒运动对裂缝扩展没有影响，首先将裂缝宽度嵌入到控制流方程中，以考虑裂缝传播效应；然后，流体-颗粒相互作用是完全耦合的，为简化计算，只考虑曳力和流体-颗粒相互作用力等。具体计算程序如图1所示，其中上面的箭头表示3种物理之间的数据传输。上标“n”和“n+1”分别表示第n个和（n+1）个时间步；q滤失是 流体滤失量；w是 裂缝宽度；φp是颗粒的体积分数；分别是流体和颗粒速度。

图1 计算程序流程图

2 实验对比验证

通过与实验结果对比验证方法的可行性，随后利用到大规模问题当中，由简单边界到复杂边界条件出发，逐步考虑在以上不同因素组合下支撑剂运移的最终分布，并作相关参数敏感性分析，找到影响支撑剂运移效果的关键因素。

2.1 物模实验

物模实验采用可视化平面裂缝模拟实验系统，实验系统主体为可视化平面裂缝模型，为透明钢化玻璃，两端设置有模拟井筒，配合流体泵和混砂罐，用于模拟支撑剂在垂直裂缝中运移及铺置。

2.1.1 实验方案

以滑溜水压裂液为实验对象，使用40/70目陶粒支撑剂（视密度为2800 kg/m3），在排量为8 m3/min，并设定砂比为6%的条件下进行物理模型实验，实验过程中保证均匀加砂，通过记录某些支撑剂的运移轨迹和时间，获得支撑剂运动速度，并对砂堤高度进行记录。

2.1.2 实验结果

在裂缝模型中出现支撑剂时开始计时，每隔可知，随着携砂液不断流入裂缝，不同时间所形成砂堤的形状大体一致，但砂堤高度随着时间逐渐增大。通过图2不难看出，裂缝模型总长L=4 m，裂缝各部位的沉砂量也差异较大，约95%的支撑剂均铺置于裂缝长度L=0～3 m的裂缝区间位置处，即近缝口位置的3/4裂缝长度内，仅有约5%的支撑剂铺置于近裂缝尖端1/4裂缝长度内。

图2 实验系统示意图

2.2 数模对比验证

为便于分析，将支撑剂运移在垂直裂缝剖面上的砂堤高度分为2类，即支撑剂沉降区和支撑剂悬浮区，分别取支撑剂体积分数为40%和10%为分界线[12]，即体积分数大于40%时支撑剂完全沉降，认为其为纯固相；体积分数在10%～40%范围内支撑剂处于悬浮状态，认为其为固液混合项；支撑剂体积分数小于10%时可以认为其为纯液相。

数值计算中参数均与物模实验参数一致，唯一不同的是在处理粒径上采取40/70目的平均粒径作为代表。时间节点分别为0.5，1，1.5，2 min的数值模拟结果如图3所示。

从图3可以看出，在砂堤长度上，大部分支撑剂均处于距入口3 m以内，这与实验结果图3也比较吻合。图4为不同时刻室内实验和数值计算的砂堤峰值高度对比，不难看出其结果相差不大。综合比较二者结果，差异较小，说明了该模型及数值计算的正确性和可靠性。

图3 砂比为6%时实验砂堤形态分布图

图4 砂比为6%时砂堤形态数值计算结果

3 超临界CO2压裂缝内支撑剂运移规律分析

二氧化碳物性参数随温度和压力的变化而变化，其相变情况复杂，在CO2压裂过程中其地面管线流动、井筒流动等温压变化较大的过程，研究多基于CO2物性模型进行耦合计算分析[24]。而本文研究的为CO2压裂过程中其进入地层后在裂缝扩展过程中的支撑剂运移规律，考虑到CO2在裂缝中的平均流速较小，且可以认为二氧化碳的密度和压力与裂缝中运移过程中的地层应力相比变化不大，为便于计算，本文计算假设CO2是弱可压缩性，主要考虑二氧化碳的平均黏度和密度对支撑剂运移的影响。

3.1 考虑流体滤失的简单矩形裂缝区域内的支撑剂运移分析

为了便于对比分析，考虑CO2流体滤失，不考虑裂缝扩展和裂缝宽度变化，考察区域简化为一个矩形区域，长160 m，高20 m，裂缝宽度设为3 mm，如图5所示。从左侧均匀速率注入，右侧定压，假设整个矩形裂缝区域内流体滤失速率一致，最终使得右侧出口流体速率由于滤失作用接近于零，计算分析支撑剂在一矩形裂缝区域内的运移情况。

图5 不同时刻实验和数值计算的砂堤峰值高度对比图

假设CO2在携砂过程中密度为600 kg/m3，初始黏度为0.1 mPa·s，假设增黏后黏度为10 mPa·s，支撑剂密度为2650 kg/m3和1000 kg/m32种，支撑剂尺寸为20/40目和40/70目两种，排量5 m3/min，即入口速率约为0.2 m/s；因此本算例共设计支撑剂密度、尺寸和携砂液黏度3个参数分别各2个取值正交共8组数值实验。计算时间共30 min，计算结果如图6所示，支撑剂沉降区和支撑剂悬浮区分别取支撑剂体积分数为40%和10%为分界线。从图6可以看出，当流体黏度较低时，即黏度为0.1 mPa·s时，均会在近井区域形成砂床。由于未经增黏，液态二氧化碳的黏度比清水还要低很多，颗粒雷诺数非常大，导致颗粒快速在近井区域沉降堆积形成砂床。由于黏度极低，优化其他参数对结果影响不大。这些近井区域支撑剂堆积预示着着实际压裂施工中的砂堵现象，泵压会急剧升高，说明携砂效果极差。

图6 简单矩形裂缝示意图

将二氧化碳增黏100倍后，携砂效果有明显改善。但是仅增黏，效果也不是很理想。当二氧化碳黏度为10 mPa·s时，相比于支撑剂直径为0.6 mm和密度为2650 kg/m3时，减小支撑剂尺寸，即直径为0.2 mm时，或者减小支撑剂密度，即密度为1000 kg/m3时，均可以提高携砂效果。但优化支撑剂密度对携砂效果的提高更为明显。

不难发现，将二氧化碳增黏，同时减小支撑剂密度和尺寸，即二氧化碳黏度为10 mPa·s，支撑剂密度为1000 kg/m3，支撑剂直径为0.2 mm时，其携砂效果最理想，仅在裂缝前端沉降堆积，且支撑剂悬浮区大面积覆盖了裂缝。

通过分析可以发现，颗粒密度会影响收尾速度，改变流场，一般密度越大颗粒收尾速度大沉降越快，然而在颗粒尺寸小的配置下，导致一方面沉降速度小而改变流场能力凸显出来，低黏条件下易产生异重流；颗粒粒径会影响收尾速度，粒径越大沉降越快，易形成砂床；携砂液黏度会影响颗粒收尾速度及流场受颗粒群扰动能力，黏度越大，颗粒沉降越慢，同时流场受到壁面黏性力增大，即使是大密度颗粒也难以影响流场。

3.2 考虑流体滤失和裂缝扩展的PKN型裂缝中的支撑剂运移分析

考虑CO2流体滤失，基于PKN型裂缝扩展模型[15,25,26]，计算分析CO2压裂裂缝中的支撑剂的运移情况，充分考虑了实际复杂边界条件。由于篇幅所限，本部分计算所用PKN型裂缝扩展模型不做详述，具体参考详见文献[15,25,26]。在PKN型裂缝的假设中，裂缝是等高的，裂缝断面为椭圆形，最大缝宽在裂缝中部，前缘不断向前扩展，如图7所示。

图7 不同模拟条件下支撑剂运移结果

3.2.1 不同支撑剂及流体黏度影响分析

考虑一定排量下支撑剂粒径、密度以及流体黏度影响，只考虑裂缝扩展，暂不考虑流体滤失的情况下。假设CO2在携砂过程中密度为600 kg/m3，初始黏度为0.1 mPa·s，假设增黏后黏度为10 mPa·s，支撑剂密度为2650 kg/m3和1200 kg/m32种，支撑剂尺寸为20/40目和40/70目2种，排量5 m3/min，即入口速率约为0.2 m/s；因此本算例共设计支撑剂密度、尺寸和携砂液黏度3个参数分别各两个取值正交共8组数值实验。计算时间共30 min，计算结果如图8所示，支撑剂沉降区和支撑剂悬浮区分别取支撑剂体积分数为40%和10%为分界线。

图8 PKN型裂缝示意图

由图8可以，看出支撑剂在扩展裂缝中的运移行为大致与在固定边界的矩形裂缝中的结果类似，但支撑剂在竖直方向上的砂堤浓度分布比简单裂缝中更为复杂，尤其支撑剂悬浮区形状在竖直方向上缝前端部更高，也说明了考虑裂缝扩展后，在裂缝高度方向上中部的速度较大。增黏后支撑剂沉降区域面积减少，说明支撑剂没有太多沉降，当支撑剂密度为1200 kg/m3，支撑剂直径为0.2 mm，二氧化碳黏度为10 mPa·s时，支撑剂大部分处于悬浮区，携砂效果是最好的，这与固定边界的矩形裂缝中的最优结果也是一致的。

3.2.2 不同排量影响分析

考虑到实际增黏效果达不到10 mPa·s这个量级，根据实际情况，假设二氧化碳增黏后黏度为2.5 mPa·s，颗粒尺寸为40/70目，密度取1200 kg/m3，排量分别为1、3、5 m3/min的3组算例以作比较。不同排量下支撑剂运移的数值计算结果如图9所示。由图9可以发现，增加黏度至2.5 mPa·s后，采用超轻支撑剂与细尺寸颗粒组合，携砂效果仍然十分理想，与10 mPa·s情况相差不大，能够满足工程需求。此外，值得注意的是，排量为3 m3/min与5 m3/min的结果从颗粒分布上看差异不大。这种情况说明在该排量范围下，支撑剂的平面铺置效果类似，而当排量过小时，支撑剂的平面铺置效果亦显然不如高排量的效果。

图9 PKN裂缝中的支撑剂运移结果

3.2.3 流体滤失影响分析

为了考虑流体滤失的影响，这里分别针对小尺寸低密度支撑剂在黏度分别为0.1 mPa·s和2.5 mPa·s 2种情形下进行模拟。参考压裂方案设计模拟计算结果，当黏度为0.1 mPa·s时流体滤失率为20.53%，黏度为2.5 mPa·s时流体滤失率为6.26%。计算结果如图10所示。

图10 增黏至2.5 mPa·s时不同排量下支撑剂运移结果

由图10可以发现，在当前滤失率条件下，支撑剂运移情况与不考虑流体滤失情况相差很小。黏度为2.5 mPa·s时，由于滤失率本身很小，为6%左右，因此运移结果十分接近；黏度为0.1 mPa·s时，尽管流体滤失20%以上，但是由于该情况下支撑剂主要在近井范围堆积，流体滤失对颗粒运移的效果在近井周围的影响并不明显，因此依然与不考虑滤失情况下支撑剂运移结果相似。

4 结论与建议

考虑CO2压裂过程中携砂液中支撑剂颗粒之间相互作用，采用欧拉-拉格朗日方法中的多相质点网格方法，建立超临界CO2压裂支撑剂运移数学模型，通过室内支撑剂运移物模实验验证模型准确性，最后对超临界CO2压裂支撑剂运移规律计算和分析，并给出相应的建议。

1.CO2压裂液的增黏是不可或缺的。未增黏CO2由于黏度低，携砂效果极差，优化其他参数对携砂效果影响不大；CO2黏度增加到2.5 mPa·s可有效提高携砂效果；采用超轻支撑剂与细尺寸颗粒组合，携砂效果与增黏到10 mPa·s效果相差不大。

2.减小支撑剂密度和尺寸可以提高CO2压裂液携砂效果。优化支撑剂密度比尺寸对携砂效果的影响更为明显；CO2黏度较低时，减小支撑剂密度和粒径可以提高携砂效果，在携砂液充分增黏的前提下可以选择目数稍大的以有效支撑裂缝。在施工成本允许的前提下可以利用一部分超轻支撑剂，与常规支撑剂交替注入，从而同时保证裂缝上端与低端同时支撑，提高最终的支撑面积。

3.增大排量可以提高携砂效果，但排量继续增大，其携砂效果变化较小；工程上不需要只为了提高携砂效果而过大的增加排量，导致摩阻增大，造成井口超压和成本升高。

4.流体滤失对CO2携砂效果影响不大；只要满足携砂和成本要求，不需要通过对CO2过度增黏来减少流体滤失。从CO2压裂与地层埋存角度出发，CO2流体在裂缝扩展过程中的滤失是有益的。