页岩毫米粗糙裂隙内支撑剂运移及裂隙损伤规律研究

Study on proppants transportation and fracture damage in shale millimeter⁃scale rough fractures

摘要

页岩储层压裂缝网中毫米裂隙是页岩气渗流的重要通道，探究毫米裂隙内支撑剂运移沉降特征及支撑剂沉降后对裂隙的支撑性能具有重要意义.基于计算流体力学与离散元、有限差分与离散元耦合方法，提出了一种综合探究支撑剂在裂隙内运移及其沉降后对裂隙支撑作用的数值研究思路.首先，通过轮廓扫描获得劈裂页岩裂隙面几何模型，建立了毫米开度裂隙内压裂液携带支撑剂运移的流固耦合数值模型，分析了不同参量下支撑剂的运移特征并获得支撑剂沉降分布；进一步，建立了毫米开度裂隙页岩数值模型，并将获得的支撑剂沉降分布导入页岩模型裂隙内，探究了应力作用下的裂隙开度演化、支撑剂受力特征以及裂隙面损伤规律.研究结果表明：支撑剂粒径越小、压裂液流速越大、压裂液黏度越高，越有利于支撑剂运移；应用混合粒径支撑剂时，大粒径支撑剂沉降会阻碍小粒径支撑剂运移，从而形成沉降区；主要承载应力的支撑剂占比较小且主要分布于裂隙面倾斜程度较高的位置；裂隙开度压缩量与裂隙面支撑剂覆盖率无明显对应关系，但受应力大小和支撑剂尺寸影响显著；应力作用下，支撑剂周围产生应力集中造成裂隙面发生拉伸和剪切破坏，拉伸破坏是裂隙面的主要损伤形式，且破坏体积随着支撑剂粒径的增加而增加.

Abstract

Millimeter⁃scale fractures in the hydraulic fracture network of shale reservoirs are important pathways for shale gas seepage. It is of great significance to investigate the characteristics of proppant migration and settlement in these millimeter⁃scale fractures and the supporting performance of proppants after settlement. Based on the coupled methods of computational fluid dynamics （CFD） with discrete element method （DEM）， and finite difference method （FDM） with DEM， a numerical research approach was proposed to comprehensively explore proppant migration in fractures and their supporting effect on fractures after settlement. First， the geometric model of the fracture surface for split shale was obtained through laser scanning， and a CFD⁃DEM model for the proppants transport carried by fracturing fluid within millimeter⁃aperture fracture was established. The transport characteristics of proppants under different parameters were analyzed and the settlement distribution of proppants was obtained. Then， a shale model with millimeter⁃aperture fracture was established， and the proppants obtained in CFD⁃DEM simulation were imported into the fracture of the shale model. The evolution of fracture aperture under stress， the stress⁃bearing characteristics of proppants， and the damage law of fracture surfaces were investigated. The research results indicate that smaller proppant particle size， higher fracturing fluid flow rate， and higher fracturing fluid viscosity are more conducive to proppant migration. When mixed proppants with different particle sizes are used， the settlement of large⁃sized proppants hinders the migration of small⁃sized proppants， thereby forming settlement zones. The proportion of proppants that mainly bear stress is small， and they are mainly distributed in areas where the fracture surfaces have large inclination. There is no obvious correlation between the compression amount of fracture aperture and the proppant coverage rate on fracture surfaces， but it is significantly affected by stress magnitude and proppant size. Under the action of stress， stress concentration occurs around proppants， leading to tensile and shear damage to fracture surfaces. Tensile failure is the main damage form of fracture surfaces， and the damage volume increases with the increase of proppant particle size.

Graphical abstract

关键词

支撑剂 / 毫米粗糙裂隙 / 运移沉降 / 裂隙压缩 / 裂隙损伤 / 耦合数值模拟

Key words

proppant / millimeter rough fracture / transport and settlement / fracture compression / fracture damage / coupled numerical method

2023年我国页岩气产量达到250亿m³，水力压裂技术无疑发挥了重要作用^［1］.如图1所示，高压水驱动下储层最大主应力方向上形成双翼对称主裂隙，并进一步发育次级裂隙.主裂隙、次级裂隙与天然裂隙相互导通形成复杂的页岩气渗流网络^［2］.为了防止压裂液反排后裂隙在应力作用下闭合，往往在压裂液中添加支撑剂^［3］.因此，支撑剂在裂隙内的运移及其沉降后对裂隙的支撑性是水力压裂技术中的重要研究内容.

为了提高对支撑剂运移的认识，文献［4］基于欧拉双流体模型建立了支撑剂与压裂液相互作用的流固耦合模型，探究了平滑垂直主裂隙内支撑剂的沉降特点；文献［5‑7］搭建了平滑垂直主裂隙试验平台，通过试验分析了支撑剂在垂直裂隙内的沉积特征及影响因素.文献［8‑9］进一步基于搭建的含多级分支的垂直裂隙模型，揭示了支撑剂在多级裂缝中的运移和沉降规律；文献［10‑11］则利用数值模拟方法分析了支撑剂在多级垂直裂隙内的沉积特征.文献［12］搭建了迂曲裂隙模型，揭示了不同施工因素下支撑剂的沉降规律；文献［13］建立了倾斜裂隙模型，获得了裂隙倾角与支撑剂沉降的关系.然而，上述研究虽获得了支撑剂在不同条件下的运移规律，但均假定了裂隙面为平滑面，难以反映支撑剂在真实裂隙中的运移.为此，文献［14‑15］搭建了支撑剂在粗糙垂直裂隙内运移的试验平台，得出了粗糙裂隙内的支撑剂沉降特征；文献［16‑17］基于数值计算方法获得了粗糙垂直主裂隙内的支撑剂运移特征.尽管上述研究考虑了裂隙面粗糙度，但仍聚焦于垂直主裂隙内的支撑剂运移.文献［18］对Barnett页岩气田的研究表明，与主裂隙相比，页岩气井的产量更受控于裂隙网络的复杂程度，其中毫米裂隙的闭合状态对页岩气产量的影响尤为显著.然而，目前对毫米裂隙内支撑剂运移规律的研究较少.

另一方面，支撑剂对裂隙的支撑性最终决定了裂隙网络的导流能力.为了探究支撑剂对裂隙的支撑性能，文献［19］开发了支撑剂支撑下裂缝闭合的数学模型，评估了应力作用下的裂隙残余开度和渗透率；文献［20］基于格子玻尔兹曼方法（LBM）分析了应力作用下裂隙的渗流能力，获得了支撑剂粒径与渗透率的关系.上述研究均假定裂隙被支撑剂完全充填，实际上支撑剂在裂隙内的充填状态受多方面因素影响.在试验方面，文献［21］制备了含支撑剂的页岩试样，利用MTS压力机探究了静力压缩过程中裂隙开度的演化特征；文献［22］搭建了压力加载台，得出了人工单层支撑剂填充裂隙开度随压力的变化规律；文献［23］自主设计了环形钢套，探究了支撑剂在单轴压缩和循环加载过程中的强度变化特征；文献［24］基于低场核磁共振成像，分析了应力加载条件下，支撑剂充填裂隙的应力敏感性和渗流能力.上述研究提高了对支撑剂在裂缝中力学行为的认识，但受限于设备水平和研究手段，裂隙中的支撑剂是人工填充，这与压裂液携带作用下的支撑剂分布存在明显差异.支撑剂在裂隙内的分布决定了支撑剂对裂隙的支撑能力，因此研究压裂液携带作用下沉降支撑剂对裂隙的支撑性更具意义.

为此，本文以毫米裂隙为关注点，提出了一种综合探究支撑剂在毫米裂隙内运移及其沉降后对裂隙支撑作用的数值研究思路.首先，通过轮廓扫描获得人工劈裂页岩裂隙面的几何模型，并基于计算流体力学和离散元（CFD⁃DEM）耦合方法，建立毫米开度裂隙内压裂液携带支撑剂运移的流固耦合模型，分析了不同支撑剂粒径及配比、压裂液流速和黏度条件下支撑剂的运移特征并获得了支撑剂沉降分布规律.然后，基于有限差分与离散元（FDM⁃DEM）耦合方法，建立了毫米裂隙页岩模型，并将获得的沉降支撑剂导入页岩模型裂隙内，分析了应力作用下裂隙开度演化过程、支撑剂受力特征以及裂隙面损伤规律.研究结果有助于提高对支撑剂在毫米级水平裂隙内的运移沉降特征以及应力作用下沉降支撑剂支撑裂隙机理的认识.

1 研究方法与流程

本文首先基于CFD⁃DEM耦合方法，借助OpenFOAM和Liggghts软件探究了支撑剂在裂隙中的运移特征，然后基于FDM⁃DEM耦合方法，借助FLAC3D和PFC3D软件探究了支撑剂对裂隙的支撑作用.

1.1 CFD⁃DEM耦合方法

流体流动由Navier⁃Stokes方程（N⁃S方程）控制.对于不可压缩流体，在含颗粒的流场中，质量守恒方程和动量方程如式（1）~（2）所示^{［9， 25］}.

\frac{\partial {\alpha }_{\mathrm{f}}}{\partial t}+Δ·（α_fu_s）=0，（1）

\frac{\partial {\alpha }_{\mathrm{f}}{\rho }_{\mathrm{f}}{u}_{\mathrm{f}}}{\partial t}+\mathrm{\Delta }·({\alpha }_{\mathrm{f}}{\rho }_{\mathrm{f}}{u}_{\mathrm{f}})=-{\alpha }_{\mathrm{f}}\mathrm{\Delta }·p-K_{\mathrm{s}\mathrm{f}}(u_{\mathrm{f}}-u_{\mathrm{s}})+

Δ·（α_fτ）+α_fρ_fg，（2）

式中：α_f为流体网格中的流体体积分数；t为流体运移时间，s；Δ为Hamilton算子；ρ_f为流体密度，kg/m³；u_f为流体速度，m/s；u_s为颗粒速度，m/s；p为压力，Pa；τ为流体应力，Pa；K_sf为固体颗粒与流体单元的动量交换系数；g为重力加速度m/s². K_sf根据文献［26‑27］提出的拽力模型计算，在流体单元的体积分数超过0.8时，采用Gidspow模型，在不超过0.8时采用Ergun模型，如式（3）所示.

K_sf=\left\{\begin{array}{ll}\frac{\mathrm{3}}{\mathrm{4}}{C}_{\mathrm{d}}\frac{{\alpha }_{\mathrm{f}}(\mathrm{1}-{\alpha }_{\mathrm{f}})|u_{\mathrm{f}}-u_{\mathrm{s}}|}{D_{\mathrm{s}}}{\alpha }_{\mathrm{f}}^{-\mathrm{2.65}},& {\alpha }_{\mathrm{f}}>\mathrm{0.8},\\ \mathrm{150}\frac{(\mathrm{1}-{\alpha }_{\mathrm{f}}{)}^{\mathrm{2}}{\nu }_{\mathrm{f}}}{{\alpha }_{\mathrm{f}}{D}_{\mathrm{s}}^{\mathrm{2}}}+\mathrm{1.75}\frac{(\mathrm{1}-{\alpha }_{\mathrm{f}})|u_{\mathrm{f}}-u_{\mathrm{s}}|}{D_{\mathrm{s}}},& {\alpha }_{\mathrm{f}}\le \mathrm{0.8},\end{array}\right.（3）

式中：D_s为固体颗粒的直径，m；ν_f为流体的动力黏度，Pa·s；C_d为颗粒的拖拽力系数，按式（4）计算^［25］.

C_d=\frac{\mathrm{24}}{{\alpha }_{\mathrm{f}}R{e}_{\mathrm{s}}}［1+0.15（α_fRe_s）^0.687］,（4）

式中 Re_s为固体颗粒的雷诺数，Re_s按式（5）计算^［25］.

Re_s=\frac{|u_{\mathrm{f}}-u_{\mathrm{s}}|D_{\mathrm{s}}}{{\nu }_{\mathrm{f}}}.（5）

CFD和DEM的双向耦合计算过程，如图2所示，颗粒与流体单元的一个耦合过程可以简单分为4步：1）生成颗粒模型并计算颗粒的速度和位置；2）颗粒的信息通过耦合程序发送至流体求解器；3）流体求解器计算流体网格的体积分数，然后根据式（1）~（2）计算流场；4）流场数据通过耦合程序将流场信息发送至DEM求解器.耦合计算过程中CFD和DEM具有不同的时间步长，一般来说，DEM的时间步长要小于CFD的时间步长.因此，为了符合物理过程，整个系统的时间步长必须统一为两者中的较小值.同时，如果每次迭代均要求实现双向耦合，则计算时间消耗巨大.为了降低计算成本，通常将CFD时间步长设置为DEM时间步长的整数倍^［28］.此外，需保证颗粒尺寸要小于流体网格尺寸^［29］.

1.2 FDM⁃DEM耦合方法

FDM和DEM耦合是基于DEM方法中颗粒与墙体的相互作用原理.颗粒无法直接与有限差分网格产生相互作用，故首先在有限元网格面边界创建墙体单元，墙体与网格面有相似的几何形状^［28］，并且依附于网格的墙面随网格的变化而变化.如果网格面为四边形，则会分为2个三角形.图3展示了颗粒与墙面之间作用力等效至网格节点的原理.墙面受到的力按照式（6）~（8）^［30］等效至每个网格面节点x_i （i=1， 2， 3），从而实现FDM⁃DEM耦合计算.

S={\sum }_{i=\mathrm{1}}^{\mathrm{3}}S_i，（6）

w_i =\frac{S_i}{S}，（7）

| F_i |=w_i | F |，（8）

式中： F 为颗粒与墙面之间的作用力矢量，N； F_i 为顶点x_i 处的作用力矢量，N；S为墙面的总面积，m²；S_i 为划分后的三角形区域面积，m²；w_i 为权重因子.

1.3 模型建立

从我国四川长宁地区选取新鲜露头的五峰组⁃龙马溪组页岩，页岩无风化现象，较致密，完整性良好.对该页岩样品进行了成分测试，X射线衍射试验结果表明该页岩的主要成分为石英、白云石和方解石，其质量分数分别为36%，29%和22%，黏土矿物的质量分数约为12%，还含有微量的其他矿物.如图4a所示，将页岩加工成尺寸为⌀50 mm×50 mm的圆柱试样，利用MTS压力机将试样劈裂.然后通过三维轮廓扫描仪获取裂隙面的三维几何模型.试验所用轮廓扫描仪型号为基恩士VR⁃5000，扫描精度为0.001 mm.仪器配套软件可实现裂隙面几何参数的提取和粗糙度计算，结果表明该裂隙面的高度范围为0~2.08 mm，裂隙面的最低处位于裂隙面中心偏右上的位置，三维高度一阶导数均方根为0.3，裂隙面面积为1 611.8 cm³.基于获取的裂隙面建立流体模型，模型上下面几何相同.模型在X和Y方向上的长度均为40 mm，裂隙开度为1 mm.模型左侧为速度入口，流体按恒定的流速流入.同时，在流体模型速度入口一侧随机生成支撑剂颗粒，颗粒与流体具有相同的初速度，随后在流体作用下移动.

在流固耦合模型计算完成后，获得了支撑剂在裂隙内的分布，进一步开展支撑剂对裂隙的支撑作用研究.如图4b所示，仍然基于扫描获取的裂隙面建立裂隙页岩模型.将流固耦合模型计算得到的支撑剂分布导入裂隙内，建立含支撑剂的裂隙页岩模型.模型高度为80 mm，顶部设置为应力边界，底部设置为固定边界，四周设置为滚支边界.

制备尺寸为⌀50 mm×100 mm和⌀50 mm×25 mm的页岩试样，开展单轴压缩试验和巴西劈裂试验；同时建立相应的数值模型，进行模型参数标定.岩石模型采用应变⁃软化本构模型，经过参数标定后的结果如图5所示.模型参数设置见表1.为了更好地对比试验结果与模拟结果，将模拟曲线沿横坐标轴平移，结果表明该页岩的单轴抗压强度约为140 MPa，巴西劈裂强度约为5 MPa.

1.4 研究方案与流程

在支撑剂运移特征的研究中，考虑了支撑剂粒径d、流体速度u_f、流体黏度ν_f和支撑剂配比4个影响因素.工程中常用的支撑剂粒径范围为0.2~0.8 mm^［9，31］，在此区间内选取了4种不同的支撑剂粒径，分别为0.25，0.27，0.38和0.55 mm.常规滑溜水压裂液的黏度范围一般为1~10 mPa·s，故选择3，5，7和9 mPa·s这4种黏度.同时，选择了0.8，1.0，1.6和2.0 cm/s共4种不同的压裂液流速.具体方案设置如表2所示，方案A（A1~A4）、方案B（B1~B3）及方案A3、方案C（C1~C3）及方案A3分别探究了支撑剂粒径、流体速度和流体黏度对支撑剂运移的影响.

此外，工程应用中通常应用混合粒径的支撑剂^{［32‑33］}，因此设置3组混合粒径支撑剂方案，如表3所示.方案D（D1~D3）中流体黏度和流体速度保持一致；D1~D3方案中均含有4种不同粒径的支撑剂，表中配比指不同粒径支撑剂质量占总支撑剂质量的比例.以方案D1为例：粒径分别为0.25，0.27，0.38和0.55 mm的支撑剂质量分别占支撑剂总质量的20%，20%，30%和30%.在模拟支撑剂运移的过程中，不同粒径支撑剂在裂隙入口处按照上述质量比例生成.墙面和支撑剂接触的切向和法向刚度设置为10 MPa，阻尼系数为0.7，摩擦系数为0.7.支撑剂的密度设置为2 500 kg/m³.DEM时间步长设置为0.1 μs，CFD时间步长设置为5 μs. CFD时间步长是DEM时间步长的50倍，也即DEM每计算50次，进行1次CFD计算.本文的研究过程假定支撑剂为刚性球形颗粒^{［10， 16， 34］}，假定压裂液流动对裂隙形态无影响.

完成支撑剂运移研究后获得支撑剂分布结果，进一步开展沉降支撑剂对裂隙的支撑作用研究.首先通过FISH语言编写脚本文件，将沉降支撑剂导入页岩模型的裂隙内.导入裂隙内后，支撑剂与裂隙面之间仍设置为线性接触模型，模型参数与前文所述的支撑剂与墙面之间的参数设置保持一致.模型顶面施加垂直应力，应力由0增加至50 MPa，以1 MPa为增量，每次增加应力后模型计算至平衡状态.为了方便叙述，根据A，B，C和D这4组支撑剂沉降结果而开展裂隙压缩研究，仍沿用原来设置的组别编号.

2 结果及分析

2.1 支撑剂在裂隙内运移及沉降特征

图6展示了不同条件下支撑剂在裂隙内的沉降过程.图6a表明相同条件下，压裂液输运4 s后，只有粒径为0.25 mm的支撑剂未在裂隙内沉降；粒径为0.27，0.37和0.55 mm的支撑剂表现出不同的沉降特点：粒径为0.27 mm的支撑剂在裂隙面边缘处产生沉降，粒径为0.38 mm的支撑剂在裂隙面的最低处产生沉降，而粒径为0.55 mm的支撑剂同时在裂隙面边缘和裂隙面最低处产生沉降.输运16 s后，粒径为0.25 mm支撑剂仅在边界处产生沉降区；粒径为0.27和0.38 mm的支撑剂均在边界和中心低凹处产生2处沉降区；而粒径为0.55 mm的支撑剂则产生3处沉降区.沉降结果表明相同裂隙面条件下，不同粒径支撑剂产生的沉降区位置具有相似的形态，且支撑剂粒径越大，越容易沉降.

图6b展示了不同压裂液流速下的支撑剂沉降过程，结果表明，压裂液流速对支撑剂的沉降有重要影响.压裂液流速为0.8 cm/s时，支撑剂在裂隙面边缘和中心低凹处均产生大量沉降.在1.6 cm/s的流速下，支撑剂仅在裂隙面边缘处产生沉降，而当压裂液流速达到2.0 cm/s时，支撑剂不发生沉降.图6c分别展示了压裂液黏度为3，7和9 mPa·s时的支撑剂沉降过程.黏度为3 mPa·s时，压裂液输运4 s后裂隙内已形成多处沉降；16 s后，支撑剂覆盖了裂隙内较大区域.观察沉降区的演化可知，沉降区在压裂液流动的反方向上扩展速度更大.黏度为7 mPa·s时，裂隙中心低凹处在支撑剂输运8 s后产生了沉降区，但在12 s时消失，表明中心沉降区并不稳定，难以沉降.而当压裂液黏度达到9 mPa·s时，支撑剂不发生沉降.

通过上述支撑剂运移特征可以总结，工程中常用的0.2~0.8 mm粒径支撑剂、压裂液泵送速率和1~10 mPa·s的常规滑溜水压裂液均对支撑剂沉降有显著影响.支撑剂粒径越小、压裂液黏度越大、泵送速率越高，越有助于支撑剂在毫米裂隙内运移.由此可得：在压裂作业初期，应用较高黏度压裂液和较高的泵入速度，同时使用更小粒径的支撑剂，有利于支撑剂在毫米裂隙内的运移.

为了定量化分析支撑剂粒径、压裂液黏度和压裂液流速对支撑剂沉降的影响，将沉降速率低于0.02 cm/s的支撑剂定义为沉降支撑剂.图7a统计了支撑剂沉降量（沉降支撑剂体积）随时间的变化，可以看出不同条件下支撑剂沉降量与时间均呈较好的线性分布规律.通过线性拟合得到支撑剂的沉降速率，并统计在图7b~7d中.支撑剂粒径为0.25，0.27，0.38和0.55 mm时的沉降速率分别为0.458，0.710，2.157和2.893 mm³/s.粒径越大支撑剂沉降量越大，粒径为0.25和0.55 mm的支撑剂的沉降速率相差6.3倍，并且当支撑剂粒径由0.27 mm增加至0.38 mm时，沉降速率增加明显.压裂液黏度为3 mPa·s时，支撑剂沉降速率为4.585 mm³/s，在所有方案中沉降速率最大，而黏度增加至9 mPa·s时支撑剂不发生沉降，说明压裂液黏度的增加会极大提高支撑剂在裂隙内的运移.压裂液流速为0.8，1.6和2 cm/s时支撑剂的沉降速率分别为2.52，0.38和0 mm³/s，表明压裂液流速对支撑剂的运移也有重要影响，压裂液流速超过2 cm/s时，粒径为0.38 mm的支撑剂难以在裂隙内沉降.

图8以方案D1为例展示了混合粒径支撑剂在3个运移时刻的俯视图，图8a中颗粒颜色代表颗粒运移速度，颗粒的大小反应了支撑剂粒径，黑色箭头的大小和方向代表压裂液流速大小和流动方向.压裂液输运3.24 s时，粒径为0.55 mm的支撑剂首先产生沉降，之后，在3.84 s时，支撑剂沉降区开始形成.观察沉降颗粒的大小及分布可知，大粒径支撑剂的沉降阻碍了小粒径支撑剂的运移，从而导致小粒径支撑剂也随之沉降.与图6a进行对比，粒径为0.25和0.27 mm的支撑剂在压裂液输运4 s后并未产生沉降，进一步印证了上述结论.分析大粒径支撑剂沉降导致小粒径支撑剂沉降的原因：从3.24 s后沉降支撑剂周围的压裂液流速和方向可知，支撑剂沉降导致周围压裂液流速降低，而图6b表明流速越低支撑剂越容易沉降，因此小粒径支撑剂沉降是由于大粒径支撑剂沉降改变了周围压裂液流速引起的.压裂液输运16 s后，裂隙面低凹处已经形成大范围沉降区.观察图8b所示的沉降区不同粒径支撑剂三维分布可以发现，沉降区右侧边缘几乎都为粒径是0.55 mm的大颗粒支撑剂.这表明混合粒径支撑剂中较大粒径支撑剂更容易在裂隙面低凹处沉降，从而阻碍了小粒径支撑剂的运移.

2.2 应力作用下支撑剂受力特征

2.2.1 支撑剂接触比变化

将上节中不同方案获得的支撑剂分布导入页岩模型的裂隙中，探究应力作用下支撑剂对裂隙的支撑作用.由于裂隙面具有粗糙度，裂隙压缩过程中不同位置处的支撑剂受力不同，因此有必要研究支撑剂的应力承载特点.定义支撑剂接触比为承载应力的支撑剂（与上、下2个裂隙面均接触的支撑剂）数量与裂隙内支撑剂总数量的比值，记为k_c，%.图9a为方案A1~A4应力加载过程中接触比随应力（σ）的变化.可以看出对于粒径为0.55，0.38，0.27和0.25 mm的支撑剂，当应力从5 MPa增加到50 MPa时，接触比分别从26.9%，15.5%，8.5%和7.5%增加到97.4%，83.3%，67.8%和63.4%，接触比随应力的增加呈非线性增加.应力增加初始阶段，支持剂接触比增加更快，这是因为初始阶段承载应力的支撑剂数量较少，承载能力较弱.而随着应力增加，承载应力的支撑剂数量不断增加，承载能力增强导致新增加承载应力的支撑剂数量增长速度降低，表现为支撑剂接触比增加幅度降低.相同应力条件下，支撑剂粒径与接触比具有正相关关系.这是因为本研究中不同方案支撑剂输运总质量保持一致，因此总质量相同条件下，方案中支撑剂粒径越大，支撑剂数量越少.因此，在相同的压裂液输运时间内，方案A4中的支撑剂数量最少，因此接触比相对更大.此外，不同粒径支撑剂接触比均未达到100%，表明并非所有支撑剂都能发挥承受应力的作用，而是存在某些支撑剂不承载应力.

图9b展示了混合粒径方案中支撑剂接触比的变化情况，图中虚线以上部分表示粒径为0.55 mm的支撑剂接触比随应力的变化情况.可以发现不同方案中0.55 mm粒径支撑剂的接触比有相似的演化过程.5 MPa应力下，3种方案中0.55 mm粒径支撑剂的接触比为75%~85%，并且在应力加载至20 MPa时均达到100%，表明20 MPa应力条件下，粒径为0.55 mm的支撑剂全部发挥应力承载作用.此时，其他粒径的支撑剂接触比几乎为0，说明20 MPa应力条件下，几乎只有粒径为0.55 mm的支撑剂发挥了应力承载作用.当应力达到25 MPa时，粒径为0.38 mm的支撑剂的接触比开始增加，表明粒径为0.38 mm的支撑剂开始发挥承载作用.随着应力由25 MPa加载至50 MPa，3种混合粒径方案中粒径为0.38 mm的支撑剂的接触比也分别达到38.1%，43%和24.1%.而不同方案中粒径为0.27和0.25 mm的支撑剂的接触比几乎为0，表明粒径为0.27和0.25 mm的支撑剂尚未发挥应力承载作用.结果表明混合粒径支撑剂在裂隙内的应力承载过程表现出逐级承载的特点，即大粒径支撑剂率先承载应力.

由于不同条件下裂隙内的支撑剂数量不同，定义裂隙内支撑剂的总投影面积与裂隙面的投影面积之比为裂隙面覆盖率，记为k_Δs，%.图9c统计了支撑剂粒径为0.38 mm的方案（A3，B1~B3，C1~C3）中k_Δs 和k_c之间的关系.计算得出7组方案中的裂隙面覆盖率介于11.8%~23.5%之间.由图9c可知，相同应力条件下，随着裂隙面覆盖率增加，支撑剂接触比呈非线性降低趋势.统计相同应力条件下k_Δs =23.5%时和k_Δs =11.8%时对应的k_c的差值，以明确k_c对k_Δs 的敏感性变化情况.应力加载至10，20，30，40和50 MPa时，对应的k_c差值分别为22.5%，32.6%，30.9%，27.9%和23.4%，表明随着应力逐渐增加，k_c对k_Δs 的敏感性表现出先增加后降低的趋势，也就是应力为20 MPa时裂隙面覆盖率对支撑剂接触比的影响更大.此外，可以发现不同条件下k_c始终未达到100%，表明裂隙内总是存在不发挥应力承载作用的支撑剂.

2.2.2 不同位置支撑剂受力分布

前文表明不同条件下支撑剂接触比始终未达到100%，因此进一步探究裂隙面不同位置处支撑剂的应力承载特点.图10a和10b以不同压裂液黏度方案（方案C及方案A3）为例展示了支撑剂受力分布.图10a中黑色颗粒点代表垂直应力超过100 MPa的支撑剂，灰色颗粒点代表应力低于100 MPa的支撑剂.统计发现受力超过100 MPa的支撑剂颗粒仅有114个，而受力低于100 MPa的支撑剂颗粒则有3 206个（占比96.6%），这表明承载较高应力的支撑剂仅占较小部分.将方案C及方案A3中受力超过100 MPa的支撑剂统计在图10b中.方案A3，C2和C3中应力超过100 MPa的支撑剂数量分别为181，145和142个，同样说明发挥主要应力承载作用的支撑剂仅占较少部分.结合图中的高度等值线分布可以发现，主要承载应力的支撑剂分布表现出明显规律性，即主要承载应力的支撑剂主要分布在裂隙面倾斜程度较大的位置，如图10b中黑色椭圆线框所示.

由图10a中的局部放大图可以看出，主要承载应力的支撑剂分布在裂隙面倾斜较大的位置.分析其原因，如图10c所示，对于开度为h的裂隙，裂隙面倾斜会导致裂隙间的垂直距离（h'）减小，并且裂隙倾斜程度越大，垂直距离越小.当支撑剂滞留在裂隙的倾斜位置时，在应力作用下，上部裂隙面率先与倾斜位置的支撑剂接触（如图中绿色支撑剂所示），而在水平裂隙中的支撑剂并未起到支撑作用.因此，在裂隙面倾斜程度较大的位置，支撑剂更容易起到应力承载作用.

2.3 应力作用下裂隙开度分布及演化

裂隙内不同位置处支撑剂受力不同，因此裂隙局部变形不同，导致裂隙开度分布不同.裂隙开度可以反映裂隙空间的压缩特点，进而反映支撑剂的支撑性能.图11以方案A为例统计了应力作用下裂隙开度分布及演化过程.可以看出，支撑剂粒径越小，裂隙开度越低.随着应力增加，裂隙开度逐渐降低，不同区域表现出不同的变化特点.局部开度较大的区域分布在裂隙面低凹处边缘，图10c可以解释上述现象：在裂隙面低凹处边缘区域的倾斜程度较大，因此边缘处支撑剂率先支撑裂隙，限制了裂隙压缩；而相对平整区域的支撑剂未与裂隙上表面接触，因此相同应力条件下，该区域的裂隙进一步压缩.A1和A2方案中支撑剂粒径相差较小，裂隙开度演化情况基本一致.A3和A4方案中支撑剂粒径较大，从裂隙开度演化情况可以看出，50 MPa应力条件下，支撑剂沉降区裂隙开度仍保持较大值，表明支撑剂沉降区能起到更好的裂隙支撑作用.

定量化探究裂隙开度变化，将裂隙的平均开度随应力的变化情况统计在图12a中，可以看出应力与裂隙平均开度之间存在线性负相关关系.支撑剂粒径越大，裂隙平均开度越大，且随应力增加裂隙平均开度的降低量也越大.定义拟合直线的斜率为裂隙平均开度的降低速率，并统计在图12b中.

由图12b可知，粒径为0.55，0.38，0.27和0.25 mm的支撑剂充填裂隙的平均开度降低速率分别为1.630，0.453，0.258和0.229 μm/MPa，支撑剂粒径与裂隙平均开度降低速率表现出非线性关系.图10表明并非裂隙内所有支撑剂均发挥承载作用，因此也统计了方案B和C的结果，如图12b中虚线框中内容所示.裂隙面覆盖率为11.8%~23.5%时，裂隙平均开度降低速率基本一致，表明支撑剂数量对裂隙平均开度的降低速率影响较小，从而再一次说明关键位置处的支撑剂发挥了主要的应力承载作用，而大部分支撑剂并未起到支撑作用.

2.4 裂隙面应力演化及损伤特征

2.4.1 裂隙面拉应力分布

一般情况下页岩的抗拉强度低于剪切和压缩强度，因此静力作用下岩石主要发生拉伸破坏.从软件中提取了支撑剂颗粒与裂隙面之间的接触力和裂隙面的拉应力分布，分析裂隙面的拉应力分布特征.图13以50 MPa应力作用下的A4方案为例展示了裂隙下表面的拉应力分布和接触力分布情况，图中圆点位置代表支撑剂与裂隙面的接触位置.

如图13中①号区域所示，对于单个支撑剂支撑裂隙的局部区域，在接触点处的岩体不受拉应力，而在支撑剂周围产生压应力.当发挥裂隙支撑作用的2个支撑剂间隔较远时，两个支撑剂之间的裂隙面上形成拉应力，如图13中②号区域所示.当2个或多个发挥支撑作用的支撑剂间距小于临界距离时，2个或多个支撑剂即可视为1个支撑单元，共同发挥裂隙支撑作用，在共同支撑单元范围内，裂隙面不受拉应力，如图13中③号区域所示.可以看出在沉降区中，接触点密集，形成一个共同支撑裂隙的支撑单元，从而使得沉降区内岩体不处于受拉应力状态.通过接触力的大小可以看出，在沉降区边缘接触力更大，而沉降区内部的接触力相对较小，结合图10c可知，这是因为沉降区边缘处的支撑剂先发挥支撑作用，导致相同应力条件下边缘处支撑剂产生的局部应力更大.

图14展示了不同粒径支撑剂充填条件下页岩模型下部裂隙面拉应力分布及其随应力增加的演化过程，图中灰色区域为受压区域.由图14可知，不同粒径支撑剂充填的裂隙，应力作用下的拉应力较大区域均主要分布在沉降区边缘，与图13所得结论一致.随着应力由10 MPa增加至50 MPa，不同粒径支撑剂充填的裂隙表现出不同的拉应力演化规律.对于粒径为0.25 mm的支撑剂充填的裂隙（方案A1），拉应力区域的分布范围随应力增加而迅速收缩，这是因为裂隙内小粒径支撑剂数量多，应力加载初期接触点分布较为分散，在接触点之间的裂隙面上形成如图13中②号区域所示的拉应力.而随应力增加，接触点位增加，裂隙面的应力分布由图13中的②号区域向③号区域转变，支撑单元逐步扩大合并、数量减少.因此，粒径为0.25 mm的支撑剂充填的裂隙拉应力分布区域随应力增加而迅速缩小.对于粒径为0.27 mm的支撑剂充填的裂隙，其裂隙内的支撑剂分布与方案A1相似，因此裂隙面的拉应力演化也与A1方案结果相似；但粒径增加导致裂隙内支撑剂数量减少，意味着相同应力条件下形成的支撑单元数量更多，因此与相同应力条件下的A1方案相比，拉应力分布范围有所增加.

对于粒径为0.38 mm的支撑剂充填的裂隙，裂隙面的拉应力范围同样随应力增加而缩小，但缩小的范围明显减小，且拉应力较高区域分布范围增加，如图14c中黑线圈出区域所示，这一现象对于粒径为0.55 mm的支撑剂充填的裂隙（方案A4）更加明显.结合图13分析原因：与方案A，B中小粒径支撑剂不同，方案C，D中支撑剂粒径更大、数量更少；在形成图13中②号区域所示的应力状态后，尽管加载的应力增加，但由于支撑剂数量少，在两支撑剂之间未形成新的接触位点，也就无法形成新的支撑单元，导致支撑单元之间无法有效合并，造成裂隙面拉应力分布范围随加载应力增加的缩减范围较小，且2个接触点之间的拉应力随加载应力的增加而不断累积.

2.4.2 裂隙面塑性破坏区分布

在应力作用于充填支撑剂的裂隙时，由于支撑剂与裂隙面近似为点接触，因此地应力恢复过程中裂隙面会产生局部应力集中.在岩石弹性阶段变形量小，裂隙开度变化量很小，而裂隙面局部应力集中超过岩石强度后发生的塑性破坏是导致裂隙闭合的主要原因.因此，分析裂隙面的塑性破坏过程具有重要意义.

图15展示了A4和A1方案中应力作用下裂隙面塑性破坏区域的分布及演化，图中黑色线框内为支撑剂沉降区.应力作用下，裂隙发生拉伸、剪切和拉剪混合破坏，拉伸破坏是裂隙面损伤的主要形式.方案A4中，裂隙面拉伸损伤区随着应力增加而显著增加，A1方案中的拉伸损伤区则不明显，这与图13中的拉应力分布情况相对应.观察A4方案中沉降区内裂隙损伤范围可以发现，裂隙损伤范围随应力增加变化较小.这是因为应力增加使得沉降区形成一个支撑单元，沉降区范围内的裂隙面不受拉应力，因此沉降区内的裂隙损伤较小.由于其他位置处支撑剂分布稀疏，支撑单元之间产生较高的拉应力，形成较大的拉伸破坏范围.A1方案中支撑剂数量较多，裂隙内支撑剂形成更大范围的支撑单元，因此在应力加载至30 MPa时，裂隙面损伤分布基本达到稳定.在50 MPa应力下的裂隙面损伤图中显示了支撑剂分布位置，可以看出在每个支撑剂与裂隙面的接触点均产生了剪切或拉剪混合破坏.

2.4.3 塑性破坏体积统计

图16统计了岩石模型中的塑性破坏体积.图16a为不同支撑剂粒径条件下的岩石破坏体积，可以看出低应力时，4种不同粒径支撑剂充填裂隙的塑性破坏体积介于120~200 mm³之间，差异较小.随着闭合应力的增加，塑性破坏体积表现出不同的变化特征.粒径为0.25和0.27 mm的支撑剂充填的裂隙具有相似的塑性破坏体积变化规律：随着应力增加，破坏体积的增长速度逐渐降低，在应力超过25 MPa时，破坏体积开始趋于稳定，这与图15中A1方案裂隙面损伤演化情况一致.对于粒径为0.38 mm的支撑剂充填的裂隙，其破坏体积随着应力增加呈线性增加趋势；而对于粒径为0.55 mm的支撑剂充填的裂隙，其塑性破坏体积在15 MPa后增加迅速.对比50 MPa应力下的破坏体积大小可知，支撑剂粒径越大，裂隙面损伤程度越大.

图16b统计了剪切破坏体积的占比.可以看出应力加载过程中，剪切破坏程度明显低于拉伸破坏程度.随着支撑剂粒径增加，剪切破坏占比增加；随应力增大，剪切破坏占比均表现出先减小后增加的趋势.例如方案A4中剪切破坏占比由5 MPa应力时的0.22降低至25 MPa应力时的0.16，再增加至50 MPa时的0.23.这是因为应力加载初期，裂隙内支撑单元数量较多，支撑单元之间形成较高的拉应力分布区，随着加载应力增大拉应力区的破坏增加；而当应力进一步增加后，支撑单元相互合并形成更大的支撑单元，裂隙面拉应力分布区域减少，同时支撑剂颗粒周围产生的剪应力达到岩体剪切破坏强度，支撑剂颗粒处的岩体发生剪切破坏，因此剪切破坏占比开始增加.

图16c为不同加载应力条件下粒径为0.38 mm的支撑剂充填裂隙时裂隙面k_Δs 与塑性破坏体积的关系.塑性破坏体积与支撑剂的充填量并未表现出严格的对应关系，仅大致呈现出裂隙面覆盖率越大，破坏体积越大的趋势.在50 MPa应力条件下，裂隙面覆盖率最高（23.5%）和最低（11.8%）时的塑性破坏体积均非最值，而覆盖率为19.2%时的塑性破坏体积最大，表明支撑剂的分布位置是影响裂隙面破坏程度的重要因素.不同裂隙面覆盖率下的塑性破坏体积增长特征保持一致，说明支撑剂粒径是影响裂隙面塑性破坏体积变化趋势的主要原因.

3 结　论

1）不同粒径支撑剂在相同裂隙面条件下沉降区位置形态相似，粒径越大越易沉降；流速、压裂液黏度对支撑剂沉降影响显著，小粒径、高黏度、高泵送速率利于支撑剂在毫米裂隙内运移.应用混合粒径支撑剂时，大粒径支撑剂的沉降，阻碍了小粒径支撑剂运移，不利于支撑剂在毫米裂隙内输运.

2）应力加载过程中，支撑剂接触比随应力非线性增加且与粒径正相关.混合粒径支撑充填裂隙，大粒径支撑剂先发挥应力承载作用.相同应力下，裂隙面覆盖率增加，支撑剂接触比非线性降低.发挥主要应力承载作用的支撑剂仅占裂隙内总支撑剂数量的较少部分，且主要分布在裂隙面倾斜程度较大的位置.应力作用下，裂隙平均开度与应力呈线性负相关关系，平均开度降低速度与粒径非线性相关.相同粒径、不同裂隙面覆盖率时，裂隙平均开度降低速度基本一致.

3）支撑剂在裂隙内的不均匀分布使裂隙面不同区域受力不同，拉应力主要分布于沉降区边缘.裂隙面与支撑剂接触点处不受拉应力，而在接触点周围产生拉应力.2个支撑剂间隔较远时，其间形成拉应力；间距小于临界距离时，成为一个支撑单元共同支撑裂隙，支撑单元内岩体不受拉应力.沉降区接触点密集形成支撑单元，岩体不受拉应力，故应力增加裂隙面拉应力分布区域范围减小.

4）应力作用下裂隙面主要发生拉伸破坏，剪切和拉剪混合破坏主要发生在裂隙面与支撑剂接触点位置.小粒径支撑剂充填裂隙时，破坏体积增长速度随应力增加逐渐降低，裂隙面损伤趋于稳定；大粒径支撑剂充填裂隙，裂隙面损伤程度随应力增加而增加.支撑剂粒径增加，裂隙面剪切损伤占比增加.随着应力增加，岩体剪切损伤占比先减小后增加.裂隙面覆盖率与塑性破坏体积大致正相关，但非严格关系.裂隙面关键位置处的支撑剂分布和支撑剂粒径是影响裂隙面损伤的重要因素.

参考文献

[1] 李　亭，张金发，管英柱，等.水平井分段压裂各段产能评价技术研究进展［J］.科学技术与工程，2023，23（21）：8916⁃8927.

[2] LI Ting， ZHANG Jinfa， GUAN Yingzhu，et al.Progress on each section productivity evaluation technology of horizontal well staged fracturing ［J］.Science Technology and Engineering，2023，23（21）：8916⁃8927.

[3] BARBATI A C， DESROCHES J， ROBISSON A，et al.Complex fluids and hydraulic fracturing ［J］.Annual Review of Chemical and Biomolecular Engineering，2016，7：415⁃453.

[4] 李昊擂，刘怀东，刘长友，等.巷道顶板水力裂缝卸压的相场模型及数值模拟研究［J］.中国矿业大学学报，2025，54（5）：971⁃983.

[5] LI Haolei， LIU Huaidong， LIU Changyou，et al.Phase field model and numerical simulation study of hydraulic fractures for roadway roof pressure relief ［J］.Journal of China University of Mining & Technology，2025，54（5）：971⁃983.

[6] LI J T， WEI J O， ZHOU X F，et al.Influence of proppant physical properties on sand accumulation in hydraulic fractures［J］.Journal of Petroleum Exploration and Production Technology，2021，12：625⁃632.

[7] 桑　宇，周小金，郭兴午，等.水平井压裂楔形缝内支撑剂输运特征数值模拟［J］.科学技术与工程，2023，23（23）：9911⁃9917.

[8] SANG Yu， ZHOU Xiaojin， GUO Xingwu，et al.Numerical simulation of proppant transport characteristics in horizontal well fracturing wedge fracture ［J］.Science Technology and Engineering，2023，23（23）：9911⁃9917.

[9] WU Z Y， WU C F， ZHOU L B.Experimental study of proppant placement characteristics in curving fractures ［J］.Energies，2022，15（19）：7169.

[10] 苏现波，范　渐，王　然，等.煤储层水力压裂裂缝内支撑剂运移控制因素实验研究［J］.煤田地质与勘探，2023，51（6）：62⁃73.

[11] SU Xianbo， FAN Jian， WANG Ran，et al.An experimental study on factors controlling the proppant transport in hydraulic fractures of coal reservoirs ［J］.Coal Geology & Exploration，2023，51（6）：62⁃73.

[12] 郭建春，唐　堂，张　涛，等.深层页岩压裂多级裂缝内支撑剂运移与分布规律［J］.天然气工业，2024，44（7）：1⁃11.

[13] GUO Jianchun， TANG Tang， ZHANG Tao，et al.Migration and distribution laws of proppant in multi⁃scale fractures during deep shale fracturing ［J］.Natural Gas Industry，2024，44（7）：1⁃11.

[14] 汪　杰，赵康佳，付珊珊，等.非常规储层压裂支撑剂多级裂缝运移和沉降规律［J］.天然气工业，2024，44（7）：109⁃119.

[15] WANG Jie， ZHAO Kangjia， FU Shanshan，et al.Migration and settlement laws of proppant in multi⁃scale fractures during the fracturing in the unconventional reservoirs ［J］.Natural Gas Industry，2024，44（7）：109⁃119.

[16] 郭天魁，宫远志，刘晓强，等.复杂裂缝中支撑剂运移铺置规律数值模拟［J］.中国石油大学学报（自然科学版），2022，46（3）：89⁃95.

[17] GUO Tiankui， GONG Yuanzhi， LIU Xiaoqiang，et al.Numerical simulation of proppant migration and distribution in complex fractures ［J］.Journal of China University of Petroleum（Edition of Natural Science），2022，46（3）：89⁃95.

[18] ROOSTAEI M， NOURI A， FATTAHPOUR V，et al.Numerical simulation of proppant transport in hydraulic fractures ［J］.Journal of Petroleum Science and Engineering，2018，163：119⁃138.

[19] 汪　杰，赵康佳，龚润璞，等.压裂支撑剂在迂曲粗糙裂缝中运移沉降规律研究［J］.长江大学学报（自然科学版），2025，22（1）：81⁃90.

[20] WANG Jie， ZHAO Kangjia， GONG Runpu，et al.Study on migration and settlement rules of fracturing proppant in tortuous and rough fractures ［J］.Journal of Yangtze University （Natural Science Edition），2025，22（1）：81⁃90.

[21] 陈　裕.含倾斜支撑裂缝的湘西北页岩导流能力研究［D］.长沙：中南大学，2022.

[22] CHEN Yu.Research on the flow conductivity of shale with inclined propped fractures in Northwestern Hunan ［D］.Changsha：Central South University，2022.

[23] HUANG H， BABADAGLI T， LI H，et al.A visual experimental study on proppants transport in rough vertical fractures ［J］.International Journal of Rock Mechanics and Mining Sciences，2020，134：14446.

[24] 魏东亚，陈省身，王树森，等.支撑剂在粗糙裂缝内的沉降规律的实验研究［C］//中国石油新疆油田分公司（新疆砾岩油藏实验室），西安石油大学，陕西省石油学会.2022油气田勘探与开发国际会议论文集Ⅳ.北京：《中国学术期刊（光盘版）》电子杂志有限公司，2022：231⁃237.

[25] WEI Yadong， CHEN Xingshen， WANG Shusen，et al.The experiment research about proppant settlement in the rough fracture ［C］//PetroChina Xinjiang Oilfield Company （Xinjiang Conglomerate Reservoir Laboratory），Xi'an Shiyou University，Shaanxi Petroleum Society.Proceedings of the 2022 International Conference on Oil and Gas Field Exploration and Development，Volume Ⅳ.Beijing：China Academic Journals （CD Edition） Electronic Publishing Co.，Ltd.，2022：231⁃237.

[26] 刘善勇，尹　彪，楼一珊，等.粗糙裂缝内支撑剂运移与展布规律数值模拟［J］.石油钻探技术，2024，52（4）：104⁃109.

[27] LIU Shanyong， YIN Biao， LOU Yishan，et al.Numerical simulation of migration and placement law of proppants in rough fractures ［J］.Petroleum Drilling Techniques，2024，52（4）：104⁃109.

[28] 吕显先.基于DEM⁃CFD耦合的水力压裂与支撑剂运移机理研究［D］.青岛：中国石油大学（华东），2022.

[29] Xianxian LYU.Study on mechanism of hydraulic fracturing and proppant migration based on DEM⁃CFD coupling ［D］.Qingdao：China University of Petroleum （East China），2022.

[30] FISHER M K， WRIGHT C A， DAVIDSON B M，et al.Integrating fracture mapping technologies to optimize stimulations in the barnett shale ［J］.SPE Prodution & Facilities，2005，20（2）：85⁃93.

[31] WANG H， SHARMA M M.Modeling of hydraulic fracture closure on proppants with proppant settling ［J］.Journal of Petroleum Science and Engineering，2018，171：636⁃45.

[32] 邓守春，左　鸿，李海波，等.地应力作用下水力裂缝中支撑剂导流能力的数值模拟研究［J］.煤炭学报，2017，42（增2）：434⁃440.

[33] DENG Shouchun， ZUO Hong， LI Haibo，et al.Numerical investigation of flow conductivity of propping agent in hydraulic fracture under in⁃situ stresses ［J］.Journal of China Coal Society，2017，42（Sup 2）：434⁃440.

[34] 黄　婷，翟　成，刘　厅，等.自支撑条件页岩裂缝的压缩特性分布及各向异性演化［J］.煤炭学报，2023，48（12）：4336⁃4351.

[35] HUANG Ting， ZHAI Cheng， LIU Ting，et al.Compression characteristic distribution and anisotropic evolution of shale fracture with shale self⁃supporting particles ［J］.Journal of China Coal Society，2023，48（12）：4336⁃4351.

[36] SHEKHAWAT D S， PATHAK K.Proppant's performance with reservoir rock under variable closure pressure：Results of experiments with a newly developed experimental set⁃up ［J］.Journal of Petroleum Science and Engineering，2016，147：34⁃46.

[37] BANDARA K M， RANJITH P G， RATHNAWEERA T D，et al.Crushing and embedment of proppant packs under cyclic loading：An insight to enhanced unconventional oil/gas recovery ［J］.Geoscience Frontiers，2021，12（6）：100970.

[38] 刘　厅，翟　成，赵　洋，等.基于LF⁃NMR的页岩多尺度孔裂隙应力敏感性评价［J］.煤炭学报，2021，46（增2），887⁃897.

[39] LIU Ting， ZHAI Cheng， ZHAO Yang，et al.Evaluation on stress sensitivity of multiscale pore and fracture in shale based on LF⁃NMR ［J］.Journal of China Coal Society，2021，46（Sup 2）：887⁃897.

[40] GONIVA C， KLOSS C， HAGER A，et al.An open source CFD⁃DEM perspective ［C］//5th Openfoam Workshop，Chalmers，Gothenburg，Sweden，June 21⁃24，2010：1-10.

[41] GIDASPOW D， BEZBURUAH R， DING J.Hydrodynamics of circulating fluidized beds：Kinetic theory approach ［C］//7th International Conference on Fluidization.Chicago：Illinois Institute of Technology Press，1991：75⁃82.

[42] ERGUN S， ORNING A A.Fluid flow through randomly packed columns and fluidized beds ［J］.Industrial & Engineering Chemistry，1949，41：1179⁃1184.

[43] 周　游，陈根发，陈明洪，等.全耦合CFD⁃DEM颗粒水力输送模型及其在倾斜管中的应用［J］.清华大学学报（自然科学版），2024，64（11）：1987⁃1996.

[44] ZHOU You， CHEN Genfa， CHEN Minghong，et al.Fully coupled CFD⁃DEM model for hydraulic transport of dense particles and its application in inclined pipe ［J］.Journal of Tsinghua University （Science & technology），2024，64（11）：1987⁃1996.

[45] 陈闻潇，石　崇，单治钢，等.基于OpenFOAM与PFC耦合方法的水下滑坡数值模拟研究［J］.工程地质学报，2021，29（6）：1823⁃1830.

[46] CHEN Wenxiao， SHI Chong， SHAN Zhigang，et al.Numerical simulation of subaqueous landslide based on OpenFOAM and PFC coupling method ［J］.Journal of Engineering Geology，2021，29（6）：1823⁃1830.

[47] TANG W， ZHAI C， XU J Z，et al.Numerical simulation of expansion process of soundless cracking demolition agents by coupling finite difference and discrete element methods ［J］.Computers and Geotechnics，2022，146：104699.

[48] 郭建春，路千里，何佑伟.页岩气压裂的几个关键问题与探索［J］.天然气工业，2022，42（8）：148⁃161.

[49] GUO Jianchun， LU Qianli， HE Youwei.Key issues and explorations in shale gas fracturing ［J］.Natural Gas Industry，2022，42（8）：148⁃161.

[50] 王俊亭，许冬进，陶振强，等.非常规油气支撑剂粒度粒形对其破碎率及裂缝导流能力的影响［J］.科学技术与工程，2025，25（9）：3672⁃3679.

[51] WANG Junting， XU Dongjin， TAO Zhenqiang，et al.Effects of particle size and shape of unconventional oil and gas proppants on breakage rate and fracture conductivity ［J］.Science Technology and Engineering，2025，25（9）：3672⁃3679.

[52] 陈思源，刘　浩，金　衍，等.压裂支撑剂发展综述与展望［J］.石油科学通报，2023，8（3）：330⁃346.

[53] CHEN Siyuan， LIU Hao， JIN Yan，et al.Review and prospect of fracturing proppant development ［J］.Petroleum Science Bulletin，2023，8（3）：330⁃346.

[54] 尤元元，杨仁树，左进京，等.基于相似模型和PBM⁃FEM⁃DEM耦合算法的大结构采场爆破参数优化［J］.中国矿业大学学报，2025，54（4）：796⁃810.

[55] YOU Yuanyuan， YANG Renshu， ZUO Jinjing，et al.Optimization of blasting parameters in large structure mining sites based on similarity models and PBM⁃FEM⁃DEM coupling algorithm ［J］.Journal of China University of Mining & Technology，2025，54（4）：796⁃810.

基金资助

国家重点研发计划项目(2020YFA0711800)

中国博士后科学基金项目(2025M771795)