0 引言
很多学者对压裂裂缝中支撑剂的流动和输运进行了大量的实验和数值模拟研究[7⇓⇓-10]，但主要基于可视化光滑平板裂缝模型开展研究，与真实粗糙裂缝内支撑剂的运移具有较大差距。仅有少数学者开展了支撑剂在粗糙裂缝内的运移实验研究，Raimbay等[3]采用不同岩石劈裂复刻而成的粗糙裂缝模型开展一系列水平粗糙裂缝输砂实验，证实裂缝表面粗糙度对缝内携砂液的流动和支撑剂的输送具有控制作用，支撑剂的运移和分布取决于压裂液的输送能力和表面粗糙度，其在粗糙裂缝中更容易形成非均匀分布。Huang等对粗糙垂直裂缝中支撑剂的运移进行了定量和可视化实验研究[4⇓-6]，模型制作工艺与Raimbay一致，实验发现不同类型的岩石具有不同的表面特征，导致支撑剂运移、沉降和充填裂缝的方式有所差异。徐加祥等[11]基于COMSOL Multiphysics采用数值模拟研究了支撑剂在二维迂曲微裂缝中的分布；Zhang等[12]通过对砂岩扫描获得壁面CT图像，建立3.90 cm×2.19 cm尺寸的粗糙裂缝数值三维模型，采用CFD-DEM（计算流体力学-离散单元法）耦合方法模拟了稀疏支撑剂流在裂缝中的运移，结果表明高黏度压裂液能够防止颗粒在裂缝内聚集并携带支撑剂深入裂缝。
受到岩石劈裂及复刻工艺的限制，物理模拟实验制作裂缝模型较为复杂且模拟案例单一，而数值模拟中的裂缝与真实裂缝形态存在一定的差距：目前这两种模拟的研究结果均具有一定的局限性。因此，本文基于分形插值理论，建立了不规则真实复杂粗糙壁面裂缝的生成方法，考虑了颗粒-流体、颗粒-颗粒及颗粒-壁面之间的相互作用，形成了支撑剂-压裂液固液双向耦合流动计算模型，并结合相关实验数据完成了模型的可靠性验证，在此基础上开展了不同粗糙度裂缝的输砂模拟，分析了裂缝壁面粗糙度对支撑剂起堆过程的影响，揭示了粗糙裂缝内支撑剂的运移沉降微观机理，可为现场粗糙裂缝输砂施工参数优化提供参考。
1 基于分形理论的粗糙裂缝生成方法
压裂裂缝复杂的壁面粗糙结构是导致实际输砂过程中支撑剂铺置异于理想平板模型的主要原因。研究表明岩石断裂表面几何要素具有分形特征[13⇓⇓-16]，故可采用分形理论[14]进行描述。孙洪泉[14]提出了改进自仿射分形插值方法，可生成比较精确的岩石断裂粗糙表面。利用文献[14]中的岩石表面数据，建立200 mm× 200 mm的单侧粗糙裂缝面，然后对裂缝结构面以局部邻域网格为1进行重建，可得图1a 所示粗糙裂缝结构面，进一步增加局部插值邻域网格，可得图1b 所示插值分形裂缝结构面。为使模拟更贴近实际，基于Huang等输砂实验中所用粗糙裂缝物理模型的相关数据[4⇓-6]，重构分形粗糙裂缝结构面，如图1c 、图1d 所示。
以生成的粗糙裂缝结构面为基础面，将上下表面节理匹配合成几何数值裂缝，考虑网格划分及计算收敛性等问题，建立的粗糙裂缝模型上下粗糙曲面无互相接触，且整个裂缝面上的缝宽处处相等（见图2 ）。具体生成方法为：将插值计算得到的粗糙裂缝结构面数据矩阵导入建模程序，以粗糙裂缝结构面的底面为基准面确定裂缝面上各点的垂向坐标，再根据设定的缝宽，得到裂缝的另一面坐标值，依次生成点、线、面，最后成体，随后通过计算将所有的边界融合，上下裂缝结构面之间形成缝宽均匀的流动通道，完成上下粗糙裂缝结构面的叠合，形成节理匹配粗糙裂缝模型（见图3 ）。
利用上述生成方法，以图1 中的4个粗糙裂缝结构面为基础面，生成长度、宽度、缝宽分别为200，200，1 mm的粗糙裂缝模型，如图4 所示。
2 粗糙裂缝内支撑剂-压裂液两相耦合数值模型
支撑剂在粗糙裂缝内压裂液中的运移属于固液两相流，其中压裂液以连续相形式存在，支撑剂颗粒以离散相形式在压裂液中分散。利用连续性方程、动量方程及湍流方程进行压裂液流体流场的表征，利用固相方程进行支撑剂颗粒运动属性的分析。充分考虑支撑剂颗粒间、颗粒与粗糙壁面间、压裂液与支撑剂颗粒的相互作用，将流体的流场与固相颗粒的运动进行双向耦合，从而实现支撑剂颗粒在粗糙裂缝中动态充填的数值模拟。
2.1 固相方程
压裂液与支撑剂在裂缝模型内自由流动，其中单个支撑剂颗粒的运动可分解为平动和转动两种[17]。
颗粒间的接触可通过非线性黏弹性Hertz模型表征，颗粒法向接触力是法向重叠量δn的函数[17]：
颗粒法向弹性阻尼表达式[17]为：
其中
切向接触力取决于切向重叠量和切向刚度，可由Mindlin-Deresiewicz接触理论[17]计算：
颗粒切向弹性阻尼表达式[17]为：
由此可得，第j颗支撑剂颗粒作用在第i颗支撑剂颗粒上的接触力可表示为[17]：
在计算颗粒与壁面之间的接触力时，可假定壁面为一个无限大的静止球形，仍采用上述作用力模型进行计算。
2.2 固-液两相耦合方程
压裂液携带支撑剂在裂缝内流动，压裂液对支撑剂颗粒的作用力主要为曳力和浮力两种。支撑剂颗粒随流体流动时受到的曳力由压差阻力和摩擦阻力组成，是流体、颗粒两相发生动量交换的主要原因，其中曳力可以统一表述为[9]：
根据支撑剂颗粒“稠密”与“稀疏”的分布情况，常采用Gidaspow经验模型[18]来计算曳力的大小。
颗粒受到流体的曳力，要实现CFD-DEM耦合，必须计算颗粒与流体的动量交换源项，动量交换源项代表颗粒对流体的阻力（相当于流体对颗粒曳力的反作用力），单位体积内的阻力表达式为[19]：
2.3 数值模型验证
考虑到计算的复杂性，为提高计算效率，对相关模拟参数进行适当简化。模拟时将孔眼注入口简化为一个宽0.001 m、高0.038 m的近似矩形入口，以3#粗糙裂缝模型为例（见图7 ），注入口在裂缝模型的左侧壁面（距模型底面10 mm处），出口位于整个模型的右壁面，同时建立了与实验模型完全等尺寸的光滑平板模型。模型验证基础数据如表1 所示。
表1 数值模型验证模拟方案 |
| 模型 | 裂缝缝长/ mm | 裂缝缝宽/ mm | 裂缝缝高/ mm | 粒径/ mm | 颗粒密度/ （kg•m−3） | 支撑剂质量 分数/% | 模拟颗粒数量/ （颗•s−1） | 黏度/ （mPa•s） | 实验排量/ （mL•s−1） | 模拟注入速度/ （m•s−1） | 裂缝 类型 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 实验1 | 200 | 1 | 200 | 0.375 | 2 630 | 10 | 2.96 | 9.84 | 4# | ||
| 模拟1 | 200 | 1 | 200 | 0.375 | 2 630 | 14 700 | 2.96 | 0.25 | 4# | ||
| 实验2 | 200 | 1 | 200 | 0.375 | 2 630 | 10 | 2.96 | 9.84 | 3# | ||
| 模拟2 | 200 | 1 | 200 | 0.375 | 2 630 | 14 700 | 2.96 | 0.25 | 3# | ||
| 实验3 | 200 | 1 | 200 | 0.375 | 2 630 | 10 | 2.96 | 9.84 | 光滑平板 | ||
| 模拟3 | 200 | 1 | 200 | 0.375 | 2 630 | 14 700 | 2.96 | 0.25 | 光滑平板 |
物理实验与数值模拟结果对比如图8 所示，可以看到，裂缝模型内支撑剂铺置形态与实验结果整体相似。通过对粗糙裂缝内支撑剂铺置区域与非铺置区域的前景色（黑）、背景色（白）进行填充，提取特征参数计算“支撑剂覆盖率”[4]（指支撑剂铺置平面面积占裂缝平面面积的比值），如模拟3中支撑剂填充部分前景色色阶占比为22.49%，表征支撑剂在缝内覆盖率为22.49%（见图9 ）。对比发现，物理实验中裂缝支撑剂覆盖率与数值模拟结果的平均误差仅为6.49%，证实无论是支撑剂颗粒堆积形态还是覆盖率，数值模拟结果均与实验结果具有良好的一致性，这说明本文建立的粗糙裂缝内支撑剂运移数值模型能够较为准确地模拟支撑剂在裂缝内的铺置过程。
3 裂缝内支撑剂的堆积及分布形态
3.1 单口注入时裂缝内支撑剂的运移沉降
以图4 中3#、4#粗糙裂缝模型为基础分析支撑剂的运移沉降过程，为提高计算效率，沿左侧壁面底部设置高为15 mm的注入口，注入口距底边界10 mm，当缝内支撑剂铺置不再变化时模拟停止。模拟过程中，光滑平板、3#粗糙裂缝、4#粗糙裂缝模型均采用相同的模拟参数，基本参数如表2 所示。
表2 模拟基础参数 |
| 参数名称 | 数值 | 参数名称 | 数值 |
|---|---|---|---|
| 压裂液黏度 | 1.0 mPa•s | 颗粒产生速率 | 6 000颗/s |
| 缝宽 | 1.0 mm | 注入速度 | 0.5 m/s |
| 支撑剂粒径 | 0.35 mm | 支撑剂密度 | 2 630 kg/m3 |
表3 不同粗糙裂缝内支撑剂沉降堆积过程 |
| 注入时间/s | 光滑平板裂缝 | 3#粗糙裂缝 | 4#粗糙裂缝 |
|---|---|---|---|
| 2 |  |  |  |
| 14 |  |  |  |
| 22 |  |  | 进口堵塞 模拟停止 |
通过对平板、粗糙裂缝内标记的支撑剂颗粒、颗粒群沿运移方向的平均速度统计（见图12 ），可进一步说明平板、粗糙缝内砂团的起堆过程：①光滑平板裂缝中，注入口上中下位置的3号、2号、1号颗粒同一时间注入，其中1号、2号颗粒经历加速—减速，最后速度变为0，运移迹线显示不再运动，沿堆积砂堤表面堆积；3号颗粒沿射流轨迹在堆积砂堤深处沉降，经历加速—减速—冲刷向前—减速，最后颗粒速度在16.75 s降至0，运移迹线停止，颗粒堆积；平板裂缝不同位置注入的砂粒均匀，砂粒主要按照射流轨迹分布，最终平板缝内支撑剂堆积直达稳态。②3#粗糙裂缝模型中，不同时间段注入的1、2、3号颗粒群速度变化趋势相似，最先注入的2号颗粒群0.7 s时在裂缝中部位置堆积后速度减为0；0.4 s时1号颗粒群注入，受2号颗粒群起堆砂团影响，1.3 s时在入口及砂团之间位置堆积，速度减为0；随着空隙中砂团高度的增加，后续注入的3号颗粒群运移轨迹抬升跃过砂团在背面沉降，4 s时速度基本减为0，后续注入颗粒重复该过程直至砂团形态稳定。
垂直粗糙裂缝中支撑剂的起堆过程可用图14 示意。流体携带支撑剂颗粒在狭窄粗糙通道中流动，可分解为水平运动和垂直运动，水平运动受注入压力和流动阻力控制，垂直运动受重力、浮力与下沉阻力控制。当颗粒遇到粗糙程度较大的位置（绿色阴影部分）时，水平运动被裂缝壁影响，运动速度减慢，此时颗粒在重力作用下沉降起堆，形成砂团。一旦砂团形成，将严重影响底部砂团附近的流场，阻碍颗粒的水平运动，大量支撑剂在入口与起堆砂团之间沉降，颗粒沉降越来越多，最终形成堵塞。粗糙程度越高，砂团越容易在入口处垂直堆积，直至填满整个入口。图中粗糙裂缝位置绿色区域处的粗糙程度远大于蓝色、黄色区域，此处颗粒撞击裂缝壁面几率最大，卡槽沉降可能性也最大[4]，因此砂粒在绿色区域处极易形成堵塞。同时在这个过程中，砂团将部分携砂液向上分流，携带一定量的颗粒越过砂团，最终在重力作用下，在砂团背面沉降堆积。可以看到，裂缝模型的粗糙程度越高，裂缝入口附近的支撑剂颗粒沉降速度越快，其水平运移距离越短，越倾向于在裂缝入口附近堆积，并在较短时间内形成缝内砂堵。
3.2 多口注入时裂缝内支撑剂的运移沉降
以上述研究为基础，采用光滑平板模型和图4 中的1#粗糙裂缝模型和2#粗糙裂缝模型，重新规划注入位置，在裂缝模型左侧沿壁面均匀设置3个高为15 mm的注入口，上下注入口各自距边界10 mm。光滑平板模型、1#粗糙裂缝模型、2#粗糙裂缝模型的模拟均采用相同基础参数，压裂液黏度为2.0 mPa·s，每个注入口的颗粒产生速率为6 000颗/s，其余参数借用表2 中的数据。
利用改进的立方体覆盖法[16]计算粗糙裂缝结构面的分形维数，可实现对复杂表面粗糙度的数学描述，从而间接量化裂缝数值模型粗糙度。经计算1#、2#粗糙裂缝结构面的分形维数值分别为2.084 5和2.218 9。分形维数值越大，裂缝粗糙结构表面凸起变化频率越大，粗糙程度越大；分形维数值越小，裂缝粗糙结构表面越平缓，粗糙程度越小。
表4 支撑剂沉降堆积过程 |
| 注入时间/s | 光滑平板裂缝模型 | 1#粗糙裂缝模型 | 2#粗糙裂缝模型 |
|---|---|---|---|
| 2 |  |  |  |
| 14 |  |  |  |
| 18 | 砂堤形态不再变化 |  |  |
| 22 |  | 进口堵塞 模拟停止 |
与平板裂缝模型相比，粗糙裂缝模型裂缝壁面粗糙程度高，对缝内流场的扰动大，支撑剂运移轨迹不稳定，颗粒之间的聚集、成簇趋势更为剧烈。图15 为支撑剂颗粒沿裂缝剖面的填充数与注入时间的关系曲线，分析可知，裂缝粗糙程度越高，支撑剂在裂缝进口附近堆积越快，裂缝内部填充越稀疏，并在较短时间内堵塞裂缝进口。裂缝壁面粗糙度在一定程度上控制了流体的运移路径，改变了支撑剂填充裂缝的方式。从表4 和图15 可以看出：①2#粗糙裂缝模型裂缝内支撑剂的覆盖面积要大于1#粗糙裂缝模型，裂缝壁面粗糙度对提高支撑剂覆盖面积起有利作用，壁面越粗糙，支撑剂覆盖面积增长越迅速、铺置高度增加越快；②壁面更粗糙的2#粗糙裂缝在注入口发生了砂堵，输砂困难。
光滑平板模型、1#粗糙裂缝模型、2#粗糙裂缝模型的最终支撑剂覆盖率分别为24.68%、10.57%、43.57%。结合表4 裂缝内支撑剂起堆过程可以发现，2#粗糙裂缝模型壁面粗糙程度最高，支撑剂颗粒沉降快，支撑剂覆盖率最大；光滑平板模型没有壁面微凸的影响，注入开始时主要受重力影响直接沉降，支撑剂覆盖率居中；1#裂缝模型壁面微凸抬升了支撑剂颗粒流轨迹，导致大部分颗粒随携砂液流出裂缝模型，支撑剂覆盖率小于光滑平板模型。在粗糙裂缝中支撑剂的沉降并不总是遵循重力方向，粗糙壁面使携砂液在凸起接触点附近发生转向流动，在一定程度上粗糙裂缝中支撑剂分布覆盖范围更大，尤其沿缝高堆积现象明显。因此，裂缝壁面粗糙度对支撑剂运移沉降的影响不一，粗糙裂缝输砂施工过程中，需要根据实际情况合理优化参数，趋利避害。
4 结论
采用分形插值构建真实粗糙壁面形态裂缝，同时考虑颗粒-颗粒、颗粒-壁面、颗粒-流体的相互作用所建立的CFD-DEM耦合固液两相流动模型，经实验数据验证可以较好地匹配粗糙裂缝内支撑剂的运移情况及堆积过程。
支撑剂在粗糙裂缝内输送，壁面粗糙凸起会显著改变支撑剂水平运移轨迹，裂缝模型的粗糙程度越高，裂缝入口附近的支撑剂颗粒沉降速度越快，其水平运移距离越短，越倾向于在裂缝入口附近堆积，并在较短时间内形成缝内砂堵。
裂缝壁面粗糙度在一定程度上可控制流体的运移路径，改变支撑剂填充裂缝的方式。一方面壁面凸起抬升了支撑剂运移轨迹、使支撑剂流出裂缝，导致覆盖率减小。另一方面携砂液易在粗糙壁面凸起接触点附近发生转向流动，可在一定程度上扩大支撑剂分布覆盖范围。
符号注释：
Cd——曳力系数，无因次；E——碰撞恢复系数，无因次；Fd——流体对固相颗粒施加的总曳力，N；Fd,i——流体对颗粒i施加的曳力，N；Fn——颗粒法向接触力，N；Fnd——颗粒法向弹性阻尼，N；Ftd——颗粒切向弹性阻尼，N；Fc,ij——第j颗支撑剂颗粒作用在第i颗支撑剂颗粒上的接触力，N；Fpf,i——流体对颗粒i的作用力，N；Ft,ij，Ftd,ij，Fn,ij，Fnd,ij——支撑剂颗粒i和j之间的切向弹性力、切向阻尼、法向弹性力、法向阻尼，N；Fp——单位体积内固相颗粒对流体的阻力，N/m3；Ft——切向接触力，N；g——重力加速度，m/s2；G——当量剪切模量，Pa；i，j——颗粒的编号；Ii——颗粒i的转动惯量，kg·m2；mi，mj——颗粒i、颗粒j的质量，kg；m*——等效质量，kg；Mt,ij、Mr,ij——由切向力和滚动摩擦引起的力矩，N·m；n——计算网格中包含的颗粒数目；r——粒子半径，m；R——有效颗粒半径，m；sg——颗粒密度与4 ℃水的密度之比，无因次；sn——法向刚度，N/m；st——切向刚度，N/m；t——时间，s；uf——颗粒位置的流体速度，m/s；up——粒子的速度，m/s；vi——颗粒i的线速度，m/s；vn，vt——颗粒相对速度的法向分量、切向分量，m/s；ΔV——计算单元格体积，m3；Y——颗粒弹性模量，Pa；z——颗粒i接触其他颗粒、固体壁面的总次数；β——中间变量，由恢复系数决定，无因次；μpi，μpj——颗粒i、j的泊松比，无因次；δn——法向间距，m；δt——切向间距，m；ωi——颗粒i的角速度，rad/s。