0.   引言

页岩油是一种在中国存储量较丰富且仅次于石油的油气资源^[1]。为了满足中国日益增长的能源需求，大力发展页岩产业和页岩开采技术是十分必要的。射流破岩技术能够大面积且有导向地碎裂储层岩体，方便页岩气和煤层气等油气资源的解吸与渗流，已被广泛应用于页岩的开采领域^[2-3]。

目前，国内外学者已对磨料水射流破岩问题展开了大量的研究^[4-6]，但由于射流和岩体具有不同的物理属性和力学性能，尤其是岩体中还存在着大量的孔隙^[7]，再加上磨料颗粒对岩体的冲蚀作用，使得岩体的变形破损机制变得极其复杂^[8]。因此，开展磨料水射流冲击破岩的研究，并对含有孔隙结构岩体的破损过程进行分析，能够为射流破岩开采技术的优化和完善提供一定的理论依据。

随着计算机技术的发展，数值模拟方法已经成为研究射流破岩问题的有效手段之一^[9-10]。黄中伟等^[11]基于有限元法模拟了液氮磨料射流冲击岩体的过程，发现冲击是岩体破损的主要因素。庄欠伟等^[12]应用有限元方法模拟研究了磨料水射流冲击钢筋和混凝土的过程。穆朝民等^[13]使用网格法构建了岩体，并分析了在受到磨料水射流冲击后岩体的损伤演化过程。由于网格类的数值方法在模拟大变形等问题时存在着网格畸变的情况，因此，近年来无网格方法得到了大量的关注并广泛应用于射流破岩问题的研究。

光滑粒子流体动力学（smoothed particle hydro- dynamics，SPH）方法作为一种典型的无网格粒子方法，在处理大变形、应力局部化等问题中具有一定的优势，被广泛采纳并应用于冲击大变形和流固耦合等领域。刘勇等^[14]运用SPH方法模拟分析了磨料空气射流破岩的冲蚀过程。李世杰等^[15]采用SPH方法模研究了不同直径的水射流对土壤冲蚀坑的影响。赵健等^[16]将SPH方法和有限元方法进行耦合，模拟分析了磨料水射流冲击下岩体的损伤演化机制。虽然目前有许多学者已采用SPH方法模拟研究了不同的射流破岩问题，但针对射流冲击下含孔隙岩体的破损过程以及孔隙结构对破岩效果影响的研究相对较少。鉴于实际岩体中含有孔隙结构，本文在岩体模型的内部引入了孔隙，模拟研究了磨料水射流冲击含孔隙岩体的破损过程，并分析了多种孔隙特性对射流破岩效果的影响。

本文依据SPH方法的基本理论，分别引入描述射流和岩体力学特性的NULL型和H-J-C型本构模型，建立磨料水冲击含孔隙岩体的数值模型。通过模拟磨料水射流破岩问题，验证了方法和模型的有效性。在此基础上，加入随机混合的孔隙，构建出磨料水射流冲击含孔隙结构岩体的动力学模型，对射流破岩的动态过程进行模拟研究，分析了岩体的孔隙特性（孔隙尺寸、孔隙度和孔隙填充物）以及射流参数（射流速度和磨料浓度）对岩体破损效果的影响，为深入理解射流开采中岩体的破损机理和提高开采效率提供了一定的理论依据。

1.   模型的建立

1.1   SPH方法的基础理论

SPH方法是一种通过利用离散点来构造近似函数的纯拉格朗日方法。对于磨料水射流冲击含孔隙岩体的问题而言，可将含磨料颗粒的水射流和具有孔隙结构的岩体离散为一系列的粒子，通过差值求解非稳态N-S控制方程即可得到粒子的物理属性。因此，宏观变量可通过计算域内离散点的积分插值得到，质点的近似函数为^[17]

式中$f(x)$为三维坐标向量x的函数；$\mathit{\Omega }$为x的支持域；$x - x'$为粒子间距离；h为SPH粒子光滑长度，光滑长度随时间和空间变化，如图 1所示；$W(x - x', h)$为核函数，用辅助函数$\theta (x - x')$定义：

图  1  粒子近似法

Figure  1.  Particle approximation

下载: 全尺寸图片 幻灯片

式中，d为空间维数。用粒子近似法将连续形式积分方程换算为离散形式的方程，即

式中，ρ_i为粒子i的密度，m_i为粒子i的质量。

SPH方法下的N-S方程中质量守恒方程为

动量守恒方程为

能量守恒方程为

式中，v为速度，e为热能，σ为应力张量，P为压力，m为质量，ρ为密度，μ为动力黏度系数，ε为应变率，N为质点总数，x为坐标轴符号，α，β代表不同坐标轴，i，j代表不同质点。

为了保证计算中的精度和稳定性，提高计算效率，研究采用了应用较广的三次B样条核函数^[17]。

1.2   磨料和水的本构模型

研究采用了LS-DYNA中的NULL本构模型描述水和磨料，塑性材料是一种在外力的作用下，能够产生显著变形而不被破坏的材料，由于在射流破岩问题中，磨料和水较符合塑性材料的特性^[2]，因此本文在仿真模拟中将两者都视为完全塑性材料并赋予Gruneisen状态方程^[18]，即

式中$C$为冲击波速度与质点速度变化曲线截距；${S_1}$，${S_2}$，${S_3}$为冲击波与质点速度变化曲线的斜率系数；${\gamma _0}$为Gruneisen常数；$a$为${\gamma _0}$和$\mu = \rho /{\rho _0} - 1$的一阶体积修正量。水和磨料颗粒的本构模型参数^[16]如表 1所示。为了实现磨料粒子在流场中的真实分布情况，从射流粒子中进行随机抽样，并赋予磨料的属性，实现射流中磨料颗粒的随机混合。

表  1  水和磨料本构模型的参数

Table  1.  Parameters of constitutive model for water and abrasive

名称	ρ0 /(g·cm-3)	C/(m·s-1)	S1	S2	S3
水	1.0	1480	2.56	-1.986	0.228
磨料	3.5	4569	1.49	0	0.460

下载: 导出CSV 

| 显示表格

1.3   岩体的本构模型

在射流的冲击作用下，不同性质的岩体会呈现出不同程度的破损。研究选用页岩材料，并使用H-J-C本构模型描述，页岩本构模型的参数^[19]如表 2所示。H-J-C本构模型综合考虑了岩石材料损伤、应变率和静水压力对屈服力的影响，是一种能够在岩体大变形、高应变率和高压等情况下应用的模型，其等效强度通过应力、应变速率和损伤状态来描述，屈服面的方程为^[20]

表  2  页岩本构模型的参数

Table  2.  Parameters of constitutive model for shale

岩石	ρ/(g·cm-3)	G/GPa	fc/GPa	T/GPa	μcrush
页岩	2.44	13.6	0.0048	0.004	0.001

下载: 导出CSV 

| 显示表格

式中${\sigma ^*}$为岩体冲击载荷作用下等效应力与静态屈服强度之比；${P^*}$为岩体实际压力与静态屈服强度之比；${\overline \varepsilon ^*}$为冲击载荷作用下岩体应变率与静态应变率之比；$A$，$B$，$N$和c为材料的强度参数。

损伤因子$D$通过等效塑性应变和塑性体积应变累计得到

式中$\Delta {\varepsilon _{\text{p}}}$和$\Delta {\mu _{\text{p}}}$分别为一个计算循环内的等效塑性应变和塑性体积应变；${D_1}$与${D_2}$为损伤常数；${T^*} = T/{f_{\text{e}}}$，$T$为材料的最大拉伸强度，${f_{\text{e}}}$为材料静态抗压强度。

1.4   孔隙的构建

为了构建与天然岩体更接近并且能够进行数值计算的孔隙结构，研究中参考最大球算法^[21]建立孔隙的过程和随机堆积球模型^[22-24]，将孔隙离散为球体，允许孔隙相互重叠，重叠程度与孔隙自身相比较小，如图 2所示。

图  2  球体重叠程度

Figure  2.  Sphere overlap

下载: 全尺寸图片 幻灯片

在使用SPH方法将岩体与孔隙的离散后，空间中的任意一个代表孔隙中心位置的粒子，与最近岩体粒子的欧氏距离，为球心与岩体表面相切的最大球的半径，计算公式如下：

式中，$({x_1}, {y_1}, {z_1})$为孔隙中心点的坐标，$({x_2}, {y_2}, {z_2})$为岩体点的坐标。

根据岩体与孔隙的质量比α能够确定出页岩的孔隙度β，之后通过β即可确定出孔隙的个数：

式中m_a，m_w分别为岩体与孔隙的质量；V_a，V_w分别为岩体与孔隙的体积；ρ_a，ρ₀分别为岩体与孔隙的密度；n_a为页岩模型中孔隙的个数；d_w为岩体的边长，R为孔隙最大内切球的半径。

页岩等沉积岩的形成过程十分复杂，因此在建立孔隙结构需要有以下特征^[25]：①随机性。在岩石的沉积过程中，孔隙在岩体中随机生成，且在任何点都有均等形成的机会。②独立性。孔隙生成在岩体中达到平衡且静止的状态，不会受到其他孔隙的影响。③长周期性。在岩体沉积、形成的过程中，孔隙生成的方向为重力势能梯度最大的方向，选取的孔隙的位置需要符合此规律。研究中采用了蒙特卡洛法随机生成孔隙，对于空间内任意一个SPH粒子，具有相同的概率被采样，能够使需要生成的n_a个孔隙均匀的分布在岩体内部，被选取的单个孔隙位置表示为${P_i}(x, y, z)$，其描述如下：

式中，Random(L)是能够满足孔隙分布的随机函数，x，y，z为坐标信息。

在随机得到的${P_i}(x, y, z)$位置上生成一个SPH粒子之后，沿着X，Y，Z三个方向继续增加粒子，最终生成单个孔隙，演化过程如图 3所示。

图  3  孔隙生成过程

Figure  3.  Formation process of pore

下载: 全尺寸图片 幻灯片

本文离散球体的规则与SILIN等^[21]和文献[22，23]一致。为了比较不同尺寸的孔隙之间的差异，建立球体大小相等的孔隙模型，如图 4所示。孔隙的本构模型及状态方程与水相同

图  4  球体的离散与分布

Figure  4.  Dispersion and distribution of sphere

下载: 全尺寸图片 幻灯片

经过对比，发现在相同参数下的构建的岩体仿真模型与真实岩体的孔隙分布^[26]基本相同，并且通过三轴压缩试验，得到了在载荷作用下孔隙岩体模型的应力应变曲线，如图 5所示。经计算得知，该模型的弹性模量与文献[27]中真实孔隙岩体相差1%，其他力学参数基本一致，说明了构建的孔隙岩体模型能够非常接近真实岩体孔隙结构。

图  5  孔隙岩体模型的应力应变曲线

Figure  5.  Stress-strain curve of model for porous rock

下载: 全尺寸图片 幻灯片

1.5   磨料水射流冲击孔隙岩体模型

针对所研究的磨料水射流冲击含孔隙岩体的问题，进行如下的简化假设：①磨料水射流破岩过程中忽略气相环境的影响，只包含水射流、磨料颗粒、岩体和孔隙4种物质；②磨料颗粒与水粒子为相同直径的球体。

为了保证计算效率，采用了1/4的对称模型，对称面设置为SPH_SYMMETRY_PLANE对称约束，并在对称面的基础上增加SPH_NON_REFLECTING约束来模拟岩石无限大的情况。对称后岩体为50 mm×50 mm×25 mm的立方体，孔隙的直径介于0.3~2 mm，射流喷距为1 mm，模型共由424395个光滑粒子构成，其中2394个为射流粒子，包含了水与磨料。磨料水射流冲击孔隙岩体的几何模型如图 6所示。

图  6  磨料水射流冲击孔隙岩体模型

Figure  6.  Model for abrasive water jet impacting pore rock

下载: 全尺寸图片 幻灯片

研究中选取了影响岩体破损效果的以下因素进行分析：孔隙特性（孔隙尺寸ϕ、孔隙度$\beta$和孔隙填充物）以及射流参数（射流速度v和磨料浓度δ）。

1.6   模型的验证

为了验证所建模型的有效性，对李井慧^[28]所研究的磨料水射流冲击大理石过程进行了模拟，各项参数均与文献[28]保持一致，图 7为不同射流压力下破损坑截面演化图，由图 7可以看出，随着射流压力的增加，大理石的侵彻深度在逐渐的增大。

图  7  不同射流压力下破损坑演化图

Figure  7.  Evolution of damaged pits under different jet pressures

下载: 全尺寸图片 幻灯片

图 8为试验与模拟结果对比图。由图 8可知，由于真实试验的冲击时间大于仿真模拟，且试验中的磨料水射流能够持续冲击岩体，导致试验切割深度整体大于模拟结果，但岩体的切割深度随射流压力的变化趋势基本相同，均随着射流压力的增加而增大，验证了本文所构建的磨料水射流模型的有效性。

图  8  试验与模拟结果对比

Figure  8.  Comparison between experimental and simulation results

下载: 全尺寸图片 幻灯片

本文还对林晓东等^[29]所研究的磨料水射流冲击灰岩问题进行了模拟，对岩体破损形状进行分析和验证。射流直径为4 mm，射流高度为75 mm，磨料浓度为30%，均与文献[29]保持一致。射流冲击灰岩破损坑截面演化过程如图 9所示，当射流撞击岩体的瞬间，受到冲击的区域出现了破损坑，随着冲击时间的增加，破损坑逐渐扩大，最终接近子弹型。由于林晓东等^[29]采用网格法描述岩体，受到冲蚀作用的岩石网格单元破坏失效后删除，而本文采用SPH方法对岩体进行了离散，使得破损坑的形状略有差异，但冲孔截面形态基本相同，与试验的冲孔扫描结果^[29]吻合，说明使用SPH方法替代网格法构建岩体是可行的。

图  9  灰岩破损坑对比图

Figure  9.  Comparison of damaged limestone pits

下载: 全尺寸图片 幻灯片

图 10给出了射流速度模拟结果对比，由图 10可知，当射流速度小于200 m/s时，冲击速度对破损坑半径的影响较小；当射流速度超过200 m/s时，破损坑半径随着射流速度的增加急剧增大，不同模拟方法所得的岩体冲孔半径随冲击速度的变化趋势较为相似，进一步验证了本研究基于SPH方法建立的数值模型的有效性。后续关于磨料水射流冲击含孔隙岩体的研究以此模型为基础，为了更好的观察破损坑形态的变化，射流高度调整为25 mm。

图  10  射流速度模拟结果对比

Figure  10.  Comparison of simulated results of jet velocity

下载: 全尺寸图片 幻灯片

2.   破岩指标参数定义

为了全面分析岩体的破损情况，对不同工况下的破岩过程进行了模拟，并提取了岩体破损图像的特征，研究表明破损坑深和破损面积的变化规律基本一致，因此本文选取了表面坑径（r/R）和破损面积（s/S）为岩体破损指标，并进行了无量纲处理，如图 11所示，其中r，s分别为岩体破损坑的表面坑径和纵截面积，R，S分别是岩体的宽度和纵截面积。

图  11  破岩指标参数示意图

Figure  11.  Diagram of rock breaking index parameters

下载: 全尺寸图片 幻灯片

3.   Von-Mises应力分析

选取孔隙度β=10%，孔隙尺寸为1 mm³，孔隙填充物为水，射流速度v=500 m/s，磨料浓度δ=5%的数值模型进行Von-Mises应力分析，图 12为射流冲击过程的Von-Mises应力演化图。

图  12  Von-Mises应力演化图

Figure  12.  Evolution diagram of Von Mises stress

下载: 全尺寸图片 幻灯片

由图 12可知，射流撞击到岩体的瞬间，受到冲击的岩体区域出现了较大的应力集中现象，当Von-Mises应力超过岩体的强度极限时，会导致岩体与射流的接触区域破损。由于岩体部分采用了脆性较高的页岩，当t=0.03 ms时，应力波开始呈现出裂纹状；距离射流较近的孔隙边缘出现了应力集中现象，并且孔隙的形状发生了改变，因孔隙填充物为水，导致孔隙内部应力较小。随着冲击时间的增加，冲击区域两侧的孔隙持续受到岩体挤压。当t≥0.06 ms时，孔隙被完全压碎，导致无法看出孔隙具体的形状和位置，岩体的破损程度不断加剧，应力波在岩体内部向四周传播并且以裂纹状继续扩散。

4.   孔隙特性及射流参数的影响

4.1   孔隙尺寸

孔隙的尺寸会影响射流破岩的效率和岩体的破损程度。研究中取孔隙尺寸ϕ=1~6（以ϕ=3为例说明，孔隙直径离散为3个粒子），如图 13所示。孔隙度β=15%，孔隙填充物为水，磨料体积浓度δ=5%，射流速度v=500 m/s。

图  13  孔隙尺寸示意图

Figure  13.  Diagram of pore size

下载: 全尺寸图片 幻灯片

图 14给出了不同孔隙尺寸下t=0.1 ms时岩体破损区域与破损指标随冲击时间的变化曲线。由图 14（a）可看出，在相同的孔隙度下，当孔隙尺寸较小时(ϕ≤4)，岩体内孔隙数量较多，能够均匀的分布在岩体中，随着孔隙体积的增大，岩体破损加剧；当孔隙尺寸ϕ≥5时，岩体的破损坑的形状变化较大。由图 14（b）可以看出，随着射流冲击时间的增加，岩体破损坑的表面坑径和破损面积都逐渐增大。当孔隙较小时(ϕ≤4)，随着孔隙尺寸的增加，破损面积不断增大，但表面坑径之间的差异较小；当孔隙尺寸ϕ≥5时，表面坑径之间出现了较大的差异。

图  14  不同孔隙尺寸

Figure  14.  Different pore sizes

下载: 全尺寸图片 幻灯片

为了探究岩体的破损是否受到孔隙分布位置的影响，研究中对孔隙尺寸ϕ≤4的4种情况进行了多次模拟并进行了对比，发现在相同孔隙尺寸下，岩体的破损情况受到孔隙分布位置的影响较小。在孔隙尺寸较大的两种情况下（ϕ≥5），分别增加3种不同的孔隙随机分布。图 15为不同的孔隙分布位置下t=0.1 ms时岩体破损区域与破损指标随冲击时间的变化曲线。结合图 14，15可以得知，当孔隙较大并且在相同的孔隙度下，不同孔隙分布位置的破损坑形状有显著的差别。当孔隙尺寸ϕ≥5时，不同的孔隙分布位置下破损坑的表面坑径以及破损面积差异较大，说明了在磨料水射流冲击的过程中含大孔隙岩体的破损受到孔隙随机分布位置的影响，因此后续的研究中孔隙尺寸选取为1 mm³（ϕ=3）。

图  15  不同孔隙分布位置

Figure  15.  Distribution location of different pores

下载: 全尺寸图片 幻灯片

4.2   孔隙度

孔隙度是影响岩体破损的重要原因。研究中取孔隙度β=0~25%，孔隙尺寸为1 mm³（ϕ=3），孔隙填充物为水，射流速度v=500 m/s，磨料体积浓度δ=5%。图 16给出了不同孔隙度下t=0.1 ms时岩体破损区域与破损指标随冲击时间的变化曲线。从图 16可以看出，随着岩体内部孔隙的增多，在磨料水射流的冲击作用下岩体破损坑的表面坑径和破损面积均增大，尤其是破损面积急剧增加，说明孔隙率越高的岩体越容易破碎。当岩体内不含孔隙（β=0%）时，岩体的破损面积略大于含孔隙较少的岩体（0%＜β≤10%），这也说明少量的孔隙及其填充物可能会对磨料水射流产生干扰并起到一定的缓冲作用，降低了射流的冲击破岩能力。

图  16  不同孔隙度

Figure  16.  Different porosities

下载: 全尺寸图片 幻灯片

4.3   孔隙填充物

不同区域的页岩孔隙中蕴藏着不同的填充物质。研究选取了水、轻油、重油和空气4种孔隙填充物，由于空气的密度与液体相比非常小，可近似忽略，因此将孔隙粒子删除视为空气，填充物的参数^[30]如表 3所示，孔隙尺寸为1 mm³（ϕ=3），孔隙度β=15%，孔隙的分布位置完全相同；射流速度v=500 m/s，磨料体积浓度δ=5%。

表  3  孔隙填充物的参数

Table  3.  Parameters of pore fillers

材料	ρ /(g·cm-3)	动力黏度系数/(pa·s)
水	1.00	1.009×10-3
重油	1.10	2.5×10-3
轻油	0.85	2.5×10-3

下载: 导出CSV 

| 显示表格

图 17给出了不同填充物下t=0.1 ms时岩体破损区域与破损指标随冲击时间的变化曲线。由图 17可以看出，当填充物为液体时，岩体破损坑的形状较为相似，表面坑径以及破损面积差异较小；当孔隙填充物为空气时，岩体破损坑的面积最大，但表面坑径最小，这也说明与其他填充物相比，含气体的岩体更容易沿着磨料水射流的冲击方向发生破损。

图  17  不同孔隙填充物

Figure  17.  Different pore fillers

下载: 全尺寸图片 幻灯片

4.4   射流速度

射流速度是决定开采破损效率的重要因素。射流速度取v=100~500 m/s，磨料浓度δ=5%，孔隙尺寸为1 mm³（ϕ=3），孔隙度β=15%，孔隙填充物为水。图 18给出了不同射流速度下t=0.1 ms时岩体的破损区域与破损指标随冲击时间的变化曲线。由图 18可以看出，岩体破损坑的表面坑径和破损面积随着射流速度的增大而增加。当射流速度v＞200 m/s时，表面坑径和破损面积急剧增加，说明磨料水射流速度越高，破损孔隙岩体能力越强。

图  18  不同射流速度

Figure  18.  Different jet velocities

下载: 全尺寸图片 幻灯片

4.5   磨料浓度

磨料浓度也是影响射流破岩效率和破损精度的重要因素之一。根据实际磨料水射流破岩过程中的磨料颗粒体积占比，磨料浓度δ=0~25%，射流速度v=500 m/s，孔隙尺寸为1 mm³（ϕ=3），孔隙度β=15%，孔隙填充物为水。图 19给出了不同磨料浓度下t=0.1 ms时岩体的破损区域与破损指标随冲击时间的变化曲线。

图  19  不同磨料浓度

Figure  19.  Different abrasive concentrations

下载: 全尺寸图片 幻灯片

由图 19可知，随着磨料浓度的增加，岩体破损坑的纵截面形状发生了明显的改变。当水射流中添加磨料颗粒后，岩体破损坑的表面坑径急剧减小，磨料浓度越大，岩体破损坑的表面坑径和破损面积越小；当磨料浓度增大到10%以上时，各浓度之间的表面坑径以及破损面积差异较小。

为了深入分析磨料浓度对岩体破损效果的影响，给出了底部坑径（l/R）指标，如图 20所示。

图  20  底部坑径示意图

Figure  20.  Diagram of diameter of bottom pit

下载: 全尺寸图片 幻灯片

图 21为不同磨料浓度下底部坑径随冲击时间的变化曲线，由图 21可知，当水射流中混有磨料颗粒时，岩体破损坑的底部坑径急剧增加，并且随着磨料浓度的增加而增大。结合图 19（a）可知，磨料颗粒主要集中在破损坑四周的表面上并产生冲蚀作用，使破损坑的表面坑径和底部坑径相差较小，说明添加磨料后能够有效的提高射流破岩的精准性。

图  21  不同磨料浓度下底部坑径随冲击时间的变化曲线

Figure  21.  Variation curves of diameter of bottom pit with impact time under different abrasive volume concentrations

下载: 全尺寸图片 幻灯片

5.   结论

（1）模拟分析了射流冲击岩体过程的破损形态和指标，通过与已有研究结果进行对比，验证了本研究所建模型的有效性。

（2）文中采用PYTHON语言编写离散后的球体孔隙，并实现随机分布，采用该方法可以对数值模拟中的孔隙特性实现精准的控制。

（3）模拟研究了孔隙特性和射流参数对岩体破损的影响。孔隙尺寸和孔隙分布情况会显著影响岩体的破损形状；当孔隙度介于0~25%时，相同的孔隙尺寸下，岩体孔隙度越大，岩体破损坑的表面坑径和破损面积越大；当孔隙填充物之间的密度以及动力黏度系数相对接近时，填充物对岩体破损的影响相对较小，其中当填充物为气体时，岩体则更容易沿着射流冲击方向发生破损变形。射流的初始速度和磨料浓度也是影响岩体破损情况的主要因素，射流初速度越大，射流的冲击能力越强，岩体破损坑的表面坑径和破损面积也随之增大；对于多孔岩体而言，冲击过程中的磨料颗粒主要集中在破损坑四周的壁面上，随着磨料浓度的增大，磨料颗粒对破损坑壁面的冲蚀效果越强，坑底部的直径明显增大。