耦合计算流体力学与离散元方法的CFD-DEM方法是模拟裂缝封堵层结构形成过程的有效手段。根据塔里木盆地深层裂缝性储集层裂缝几何特征, 选用楔形裂缝作为裂缝模型的几何形态, 裂缝入口和出口宽度分别设定为3 mm和1 mm, 裂缝长度设定为50 mm, 接近室内实验常用裂缝模块尺寸。裂缝几何模型如图1所示。

图1

	Figure Option ViewDownloadNew Window
	图1 裂缝几何模型

本文CFD-DEM模拟中, 堵漏材料相关参数依据裂缝性储集层最常用的方解石类堵漏材料进行设置, 堵漏材料摩擦系数采用COF-1型堵漏材料表面摩擦系数测量装置^[5]测得。流体相关参数依据塔里木盆地深井钻井液性能参数和深井裂缝性地层钻井液漏失速率确定。堵漏材料和钻井液基础参数如表1所示。

表1

表1

表1 堵漏材料和钻井液基础参数 参数 数值 参数 数值 颗粒密度 2 700 kg/m3 流体密度 1 700 kg/m3 颗粒等效直径与 裂缝出口宽度比 0.4~1.0 回弹系数 0.5 颗粒入射速度 0.5 m/s 弹性模量 1 GPa 流体初始速度 0.5 m/s 泊松比 0.3 流体动力黏度 50 mPa· s 摩擦系数 0.6, 0.8, 1.0, 1.2

表1 堵漏材料和钻井液基础参数

耦合CFD-DEM模拟方法分为解析方法和非解析方法, 解析方法的颗粒尺寸明显大于流体计算单元（见图2）, 运用动态网格使局部网格细化（见图2）, 进而可对颗粒周围流场和颗粒受力进行更详细的求解, 精确计算颗粒运动行为^[15]。本文研究裂缝封堵层结构形成机制, 需特别关注堵漏材料的微观架桥和封堵行为, 因此采用解析方法开展模拟。解析CFD-DEM方法的计算域由流体域和颗粒域组成^[16]。

图2

	Figure Option ViewDownloadNew Window
	图2 解析CFD-DEM方法颗粒及流体网格示意图

对于流体域部分, 解析CFD-DEM模拟方法的流体控制方程（Navier-Stokes方程）、边界条件、初始条件以及流体-颗粒耦合条件如下^[17]。

流体连续性方程：

$\nabla \cdot {{u}_{\text{f}}}=0$ （1）

流体运动方程：

${{\rho }_{\text{f}}}\frac{\partial {{u}_{\text{f}}}}{\partial t}+{{\rho }_{\text{f}}}\left( {{u}_{\text{f}}}\cdot \nabla \right){{u}_{\text{f}}}=-\nabla p+{{\mu }_{\text{f}}}{{\nabla }^{2}}{{u}_{\text{f}}}$ （2）

流体域边界条件：

${{u}_{\text{f}}}={{u}_{\text{b}}}$ （3）

流体域初始条件：

${{u}_{\text{f}}}\left( x, t=0 \right)={{u}_{0}}\left( x \right)$ （4）

流体-颗粒界面耦合条件：

$\left\{ \begin{align} & {{u}_{\text{f}}}={{u}_{\text{p}}} \\ & \sigma \hat{n}={{t}_{\text{p}}} \\ \end{align} \right.$ （5）

对于颗粒域部分, 在考虑每个颗粒与周围其他颗粒及固体边界间的相互作用力的条件下计算每个颗粒的轨迹。颗粒的运动包括平移和转动两部分, 颗粒的平移加速度和角加速度基于相应的动量平衡计算得到。根据牛顿第二定律, 颗粒i的平移控制方程为：

${{m}_{i}}\frac{\text{d}{{v}_{i}}}{\text{d}t}=\sum\limits_{j=1}^{n}{{{f}_{\text{pp}}}_{, ij}}+{{f}_{\text{pf}, i}}+{{m}_{i}}g$ （6）

颗粒j作用于颗粒i的转矩由切向力矩和滚动摩擦力矩两部分组成。根据欧拉第二运动定律, 转动惯量为I_i的颗粒i的转动控制方程表示为：

${{I}_{i}}\frac{\text{d}{{\omega }_{i}}}{\text{d}t}=\sum\limits_{j=1}^{n}{\left( {{M}_{\text{t}, ij}}+{{M}_{\text{r}, ij}} \right)}+{{M}_{\text{pf}, i}}$ （7）

耦合计算过程中, 颗粒所占据的流体计算域需要细致的网格分辨率, 为取得较精确的模拟结果, 流体计算单元特征长度与颗粒直径之比需不超过1/10^[18]。为提高计算效率, 应用动态网格对颗粒所占据的网格进行细化。当颗粒移动后, 之前所占据的网格又恢复粗化。CFD和DEM的耦合通过流体和颗粒之间的质量、动量和能量交换来完成。CFD-DEM耦合的具体流程如图3所示。

图3

	Figure Option ViewDownloadNew Window
	图3 CFD-DEM耦合流程

为了验证CFD-DEM耦合模型的可靠性, 利用所建立的模拟方法进行了一系列颗粒沉降问题模拟, 沉降区域的大小是25d_p× 25d_p× 250d_p。球体初始位置位于顶面几何中心, 高度约240d_p, 共模拟了10种不同粒径的颗粒。将颗粒沉降数值模拟结果与Concha^[19]推导出的颗粒沉降速度方程计算结果及Lapple和Shepherd^[20]的实验数据进行了对比, 验证了本文模拟方法的可靠性（见图4）。

图4

	Figure Option ViewDownloadNew Window
	图4 颗粒沉降模拟结果与文献中方程计算结果及实验结果对比

1.3.1 材料粒径及加量对裂缝封堵效率的影响

堵漏材料架桥是裂缝封堵的关键, 架桥发生后后续材料逐渐填充堆积, 进而形成裂缝封堵层。架桥材料加量和粒度分布是决定裂缝内架桥和堵漏配方设计的重要因素。图5表明, 堵漏材料加量一定时, 材料粒径与裂缝宽度的比值（R值）直接决定了裂缝封堵效率。在堵漏材料加量为5%（体积分数）, 摩擦系数为0.8的条件下, R值取0.5, 0.6和0.7时, 模拟过程中裂缝出口不断有颗粒流出, 裂缝内未能形成有效封堵; 随着R值提高到0.8, 一段时间后裂缝出口端无颗粒流出, 裂缝内形成封堵层。R值越大, 形成封堵需要的时间越短, 裂缝封堵效率越高。

图5

	Figure Option ViewDownloadNew Window
	图5 材料粒径对裂缝封堵效率的影响（堵漏材料加量为5%, 摩擦系数为0.8）

堵漏材料加量是堵漏配方设计的重要参数, 模拟中用体积分数表示, 即材料体积与堵漏浆体积之比。加量主要影响颗粒之间相互作用, 随着材料加量增加, 颗粒之间相互作用增强, 架桥发生的可能性增加。颗粒的架桥行为是一个概率性事件, 由粒径与裂缝宽度比值、加量等因素共同决定^{[21, 22]}。裂缝封堵效率体现为颗粒形成架桥所需的时间, 本质上取决于架桥发生的概率, 即架桥概率。图6表明, 在R值为0.7、摩擦系数为0.8的条件下, 当材料加量为5%时, 裂缝出口不断有颗粒流出, 裂缝内无法发生架桥; 当材料加量增加到10%及以上时, 颗粒之间相互作用增强, 架桥概率增加。材料加量越大, 形成封堵需要的时间越短, 裂缝封堵效率越高。

图6

	Figure Option ViewDownloadNew Window
	图6 材料加量对裂缝封堵效率的影响（堵漏材料粒径与裂缝宽度的比值为0.7, 摩擦系数为0.8）

1.3.2 临界与绝对架桥加量

堵漏材料在裂缝内的架桥行为具有随机性, 架桥概率用于表征裂缝封堵过程中堵漏材料在裂缝中发生架桥封堵的难易程度。基于CFD-DEM模拟方法, 架桥概率定义为：

${{P}_{\text{b}}}\text{=}\frac{{{N}_{\text{b}}}}{N}\times 100\text{ }\!\!%\!\!\text{ }$ （8）

基于架桥概率与材料加量的关系提出了临界架桥加量与绝对架桥加量概念（见图7）。当材料加量较低时, 无法发生架桥, 即架桥概率为零; 当材料加量增加到某一临界值后, 开始有一定概率发生架桥, 此加量即为临界架桥加量; 随着材料加量继续增加, 架桥概率也逐渐上升, 当架桥概率达到100%时, 对应的加量即为绝对架桥加量; 当材料加量高于绝对架桥加量时, 裂缝内一定发生架桥, 即架桥概率为100%。

图7

	Figure Option ViewDownloadNew Window
	图7 堵漏材料临界和绝对架桥加量示意图

R值是影响临界与绝对架桥加量的关键因素。图8表明, 在摩擦系数为0.8的条件下, 随着R值的增加, 临界与绝对架桥加量均呈现下降趋势。基于临界架桥加量和绝对架桥加量可将图8分成3个区域。当材料加量和R值位于红色区域时, P_b=0, 架桥不会发生; 当材料加量和R值位于黄色区域时, 0< P_b< 100%, 有一定概率发生架桥; 当材料加量和R值位于蓝色区域时, P_b=100%, 架桥必然发生。临界和绝对架桥加量指标为架桥材料加量设计提供了理论依据, 综合考虑裂缝封堵效率和堵漏材料经济性, 架桥材料加量下限应大于临界架桥加量, 上限应略大于绝对架桥加量（考虑到架桥材料粒度降级等额外消耗）。

图8

	Figure Option ViewDownloadNew Window
	图8 不同材料粒径与裂缝宽度之比、不同加量下架桥概率图版（摩擦系数为0.8）

1.3.3 材料形状、粒径与绝对架桥加量的关系

形状是表征材料外形的重要参数, 裂缝封堵过程中形状影响材料摩擦系数从而影响裂缝封堵效率^[23]。Xie等^[24]通过改变材料摩擦系数来分析非球形颗粒在狭小通道处的架桥行为。图9表明, 在不同R值条件下, 绝对架桥加量均随材料摩擦系数的增加而降低; 当R值大于0.8时, 摩擦系数对绝对架桥加量和材料架桥行为的影响较小; 随着R值的降低, 摩擦系数的影响越来越显著。通过增加材料的不规则程度和滑动/滚动摩擦系数, 可有效降低绝对架桥加量, 提高裂缝封堵效率。

图9

	Figure Option ViewDownloadNew Window
	图9 不同摩擦系数下绝对架桥加量随材料粒径与裂缝宽度之比的变化关系

1.3.4 裂缝封堵层结构形成机制

架桥是裂缝封堵层结构形成的第一步, 也是关键的一步, 架桥形成后, 材料在架桥的基础上进一步堆积填充, 封堵层致密性和强度逐渐提高。本文通过分析材料加量、R值、形状/摩擦系数等因素的影响, 揭示了楔形裂缝内材料架桥机理, 如图10所示。当R≥ 1时, 发生单粒架桥。当0.7≤ R< 1时, 靠前的两颗粒首先呈一定角度在裂缝内架桥, 后续颗粒又进一步滞留形成稳定架桥, 即顺序双粒架桥。单粒架桥和顺序双粒架桥的发生主要受材料粒径主导。随着粒径进一步减小, 当0.5≤ R< 0.7时, 架桥以平行双粒架桥和多粒架桥为主。当R< 0.5时, 裂缝内只发生多粒架桥。平行双粒架桥和多粒架桥均受材料粒径和摩擦力共同主导。随着R值的降低, 材料绝对架桥加量升高, 架桥由粒径主导向粒径-摩擦力共同主导转化。在设计堵漏配方时, 优先选取能够发生单粒架桥和顺序双粒架桥的R值。但是在应对毫米— 厘米级的裂缝宽度时, 架桥材料粒径过大往往会造成堵漏浆泵送困难或井下工具无法使用的问题。在这种情况下, 只能首先确定可满足堵漏浆泵送要求和井下工具正常使用的R值, 然后根据R值确定绝对架桥加量, 通过增加材料加量或提高材料摩擦系数的方法实现高效封堵。

图10

	Figure Option ViewDownloadNew Window
	图10 楔形裂缝内材料架桥机理

堵漏材料进入裂缝后, 经历运移、架桥、填充而形成裂缝封堵层, 封堵层几何结构形成的快慢由架桥主导的裂缝封堵效率决定。而裂缝封堵层形成后的强度与承压演化机制则由封堵层力学结构, 即细观力链网络所决定（见图11）。本文将通过开展光弹实验来分析裂缝封堵层结构演化机制。

图11

	Figure Option ViewDownloadNew Window
	图11 裂缝封堵层几何结构与力学结构

光弹实验法通过光学干涉原理反映颗粒体系中接触力的分布, 在颗粒物质力学、岩土力学、岩石力学等领域用于研究材料体系内部应力分布问题, 是表征裂缝封堵层细观力链网络的有效方法^[26]。剪切失稳是裂缝封堵层最普遍的结构失稳方式^[27]。

基于光弹实验原理, 采用自主研制的表征裂缝封堵层细观力链网络的光弹实验系统（见图12）, 研究裂缝封堵层形成后在垂向和剪切载荷作用下裂缝封堵层细观力链网络结构演化行为。该系统裂缝面为透明有机玻璃板, 裂缝开度为5 mm, 缝面尺寸为260 mm× 260 mm, 水平剪切带宽度为20 mm, 最大剪切深度为10 mm。该系统能施加垂向载荷、水平剪切载荷, 可模拟封堵层剪切失稳模式。实验材料为聚碳酸酯, 其具有光学灵敏度和透明度高、室温下蠕变小等优势^[28]。

图12

	Figure Option ViewDownloadNew Window
	图12 表征封堵层细观力链网络的光弹实验系统[28]

2.2.1 裂缝封堵层结构承压演化机制

裂缝封堵层剪切过程中压力变化如图13所示, 选取A、B、C、D这4个加压点进行分析, 各压力点对应的细观力链网络光弹图像如图14所示。图14中越亮部分表示接触力越大, 可以看出, 力链网络结构变化是宏观封堵层受力变化的细观响应。根据光弹图像像素点的平均平方灰度梯度可计算各颗粒之间接触力大小^[29]。将大于平均接触力的力链定义为强力链。光弹实验结果表明强力链占比与封堵层受力大小呈现很好的正相关性（见图13）。因此, 强力链占比是评价细观力链网络结构稳定性的重要指标。

图13

	Figure Option ViewDownloadNew Window
	图13 裂缝封堵层剪切过程中压力变化

图14

	Figure Option ViewDownloadNew Window
	图14 4个压力点下封堵层细观力链网络光弹图像

2.2.2 裂缝封堵层细观结构优化

进一步对比了裂缝封堵层剪切过程中不同材料形状、摩擦与流体环境、材料类型、封堵层孔隙度对应的细观力链网络中强力链占比与封堵层最高承压能力（见表2）。结果表明, 不同参数条件下较高的强力链占比均对应较高的封堵层承压能力。通过降低材料圆球度, 增加摩擦系数, 加入弹性材料, 提高封堵层致密性, 可有效优化裂缝封堵层细观力链网络结构, 提高封堵层强度与承压稳定性。采用该光弹实验系统, 以强力链占比为指标, 可形成堵漏材料评价新方法。

表2

表2

表2 不同参数对应的裂缝封堵层细观力链网络结构 参数 参数值 强力链 占比/% 最大承压/ MPa 最高承压点对应的 力链网络图像 材料形状 球形 11.64 0.22 球形+ 矩形 16.20 0.35 摩擦与流体环境 干燥 16.61 0.53 油浸润 14.64 0.44 材料类型 刚性+ 弹性 13.86 0.46 刚性+ 纤维 13.78 0.41 封堵层孔隙度 8.95% 16.61 0.53 9.95% 15.50 0.48

表2 不同参数对应的裂缝封堵层细观力链网络结构