使用材料点方法研究砂岩在磨蚀性水射流作用下的矿物破坏行为及侵蚀过程

《Journal of Rock Mechanics and Geotechnical Engineering》：Mineral failure behavior and erosion in sandstone by abrasive waterjet using material point method

　　张思松|夏斌伟|周延敏|卢怡宇|周雷
中国重庆市重庆大学煤矿灾害动力学与控制国家重点实验室，邮编400044

摘要
磨料水射流（AWJ）在地下工程中具有显著的应用潜力。然而，关于磨料和高压水对砂岩矿物破坏行为的报道很少，岩石破坏的过程和机制仍不清楚。本研究讨论了现有AWJ模拟方法的局限性以及材料点方法（MPM）的优势。构建了两个真实的异质砂岩模型，并进行了AWJ实验和数值模拟。阐明了AWJ的应力传递特性和矿物的破坏行为，解释了AWJ的岩石破坏机制。研究结果如下：
(1) 磨料的影响导致应力集中和矿物断裂，从而引发剪切破坏。应力波在弹性区域传播能量，当材料发生塑性破坏时能量被耗散。
(2) 材料和结构的异质性显著影响了应力分布。机械强度较高的矿物对应力反应更强烈，岩石中的缺陷（孔洞、弱面）在加载过程中成为应力集中区。
(3) 水楔效应导致因脆弱性引起的破坏。裂纹沿矿物界面传播或穿透受冲击表面的较弱矿物。
(4) 在施加持续拉力时，水楔使矿物松动并在受冲击面上形成不规则的断裂面，从而磨料颗粒有效地清除松动的材料。

1. 引言
地下开发不断向更深的岩层推进，深地质工程（如隧道建设和采矿）取得了显著进展（Munoz等人，2016；Li等人，2017）。然而，深部地下工程在破碎硬岩方面面临诸多挑战。传统的开挖方法主要包括钻孔、爆破和机械挖掘（Wang等人，2019；Sun等人，2023）。然而，钻孔和爆破存在诸多局限性，包括对周围岩层的严重扰动、过度或不足的挖掘以及在城市地区的适用性有限（Guo等人，2023）。机械挖掘的岩石破碎效率低、工具磨损严重且成本高昂（Cheng等人，2021）。因此，绿色高效的岩石破碎技术对深部地下工程的改进至关重要。
学者们提出了高压磨料水射流（AWJ）来克服传统岩石破碎方法的局限性并提高深部岩石的挖掘效率（Liu等人，2020；Banan等人，2021）。该方法利用高压水和水磨料颗粒冲击和破碎岩体，其优势在于能量集中、对周围岩层的扰动小、截面平整以及环境适应性强，被认为是破碎硬岩的有前景的技术（Selvan等人，2019；Singh等人，2021）。因此，分析岩石破坏过程和机制对于提高高压AWJ的效率和能量利用率以及发挥该技术在地下工程中的潜力至关重要。
AWJ岩石破碎通常具有高压和高速度。许多学者使用实验和数值模拟方法研究了这项技术。Liu等人（2019）开发了一个三维（3D）数值模型，结合光滑粒子流体动力学（SPH）和有限元方法（FEM），研究了围压和射流压力对AWJ切割工具性能的影响。Zai等人（2025）使用实验和数值模拟研究了标准水射流（WJ）与AWJ在破碎人工砂岩方面的差异。Chen等人（2025）通过数值和实验方法研究了AWJ的能量密度和最佳磨料体积分数。Hu等人（2025）使用SPH和离散元方法（FDEM）分析了磨料浓度和射流速度对裂纹和能量消耗的影响，比较了反射边界和非反射边界之间的破坏差异。Zhao等人（2024）使用SPH-FEM方法研究了AWS在围压条件下的岩石切割。Ge等人（2023a，2023b）进行了在三轴应力下AWJ破碎花岗岩的实验研究，确定了不同应力水平下的切割特征以解释AWJ的岩石破坏机制。Zhou等人（2025）使用SPH-FEM研究了AWJ侵蚀下层状砂岩的破坏模式和机制。
高速硬质磨料颗粒与高压水的耦合作用是高效岩石破碎的关键。由于使用实验设备观察岩石破坏过程具有挑战性，大多数研究采用数值模拟。然而，这些模拟方法存在以下缺点：
(1) 在传统模型中，岩石通常被视为均匀或随机结构材料，无法捕捉强相（如石英、长石）与弱界面之间的关键差异，从而导致对破坏机制的解释过于概括，缺乏矿物级别的力学分析。
(2) 岩石和水相的离散化表示过于简化，对磨料颗粒使用了不现实的固定几何参数，导致断裂行为表征不准确。这是因为SPH中使用粒子模拟水射流（Monaghan，2005）。一些粒子被离散化为磨料，因此磨料和水元素具有相似的体积且呈球形。此外，由于SPH和DEM元素中邻居检测的复杂性（Karekal等人，2011），计算效率非常低，导致建模中使用的元素数量较少，水的离散化不足，计算精度低。
(3) 基于网格的方法（FEM、FDEM）在模拟断裂和大幅变形方面存在固有局限性，不适合进行岩石破碎的完整周期模拟（Fukuda等人，2021）。数值模拟中的岩石破坏和大幅变形可能导致网格元素高度变形，可能出现负体积元素，从而导致计算终止（Sulsky等人，1994）。因此，受损目标元素被删除而不是随回流一起移除，这是不现实的。
(4) 计算时间过短。由于计算效率低，计算时间通常只有几百微秒。一些研究通过降低目标材料的性质来达到所需的破坏效果。AWJ岩石破碎涉及形成岩石的矿物断裂，磨料与矿物之间以及它们的体积、质量和材料性质之间存在复杂的能量相互作用。具有不同材料性质的矿物对应力的响应不同，因此一个真实的异质目标体对于理解岩石破坏过程和机制至关重要。
我们使用了AWJ进行了岩石破碎的实验和数值模拟。使用TESCAN集成矿物分析仪（TIMA）确定了数值模拟中使用的两种砂岩的矿物组成和分布，并根据实验条件建立了一个真实的异质砂岩模型。基于材料点方法（MPM）对这两种异质砂岩的AWJ侵蚀进行了模拟。对这两种砂岩类型分别进行了WJ和AWJ的实验，并使用微纳计算机断层扫描（CT）重建了它们的侵蚀形态。在侵蚀坑的中部和底部进行了扫描电子显微镜（SEM）观察，分析了射流冲击应力波的传递、矿物应力响应以及岩石和形成岩石的矿物的破坏演变。通过评估不同位置的微观损伤形态来阐明矿物的破坏行为，以研究AWJ的侵蚀机制。

2. 数值模拟方法
2.1. 控制方程
AWJ岩石破碎是一个涉及流固耦合的非线性动力学问题。砂岩和磨料的机械行为包括弹性变形和塑性破坏，而流体的行为基于高速射流的流变过程。流体和固体接触和碰撞过程中发生动量交换。弹性理论用于描述砂岩和磨料的动态行为，该理论使用运动方程（式（1））、几何方程（式（2）和弹性本构方程（式（3））（Tang，1997）：
$$
\sigma, \rho, b, v, u, K, G, \epsilon, \varepsilon, \varepsilon_v, \varepsilon_{dev}
$$
其中σ是应力张量（MPa）；ρ是密度（kg/m3）；b是体积力（N/m3）；v是速度（m/s）；u是位移（m）；K是体积模量（MPa）；G是剪切模量（MPa）；ε是应变张量；ε_v和ε_{dev分别是体积应变张量和偏应变张量。
我们使用应变软化本构模型来模拟岩石损伤的演变及断裂的起始和传播。莫尔-库仑准则（式（4）和拉伸破坏准则（式（5）分别用于描述剪切和拉伸屈服强度：
$$
\sigma_f, \sigma_t
$$
其中fs是剪切屈服强度（MPa）；ft是拉伸屈服强度（MPa）；σ?是最大主应力（MPa）；σ?是最小主应力（MPa）；c是内聚力（MPa）；Nφ = (1 + \sinφ)/(1 - \sinφ)；φ是内摩擦角（°）；σt是拉伸强度（MPa）。
当应力超过屈服强度时，材料发生塑性变形。塑性应变的方向和大小由流动规则确定：
$$
\epsilon = \phi \cdot \left( \frac{\sigma - \sigma_f}{\sigma_2 - \sigma_1} \right)
$$
其中φ是塑性修正系数，g是势函数。
本研究采用基于连续介质力学的应变软化模型（式（7）来描述高应力下岩石的损伤：
$$
\epsilon_{\text{plastic}} = \phi \cdot \left( \frac{\sigma - \sigma_f}{\sigma_2 - \sigma_1} \right)^n
$$
当塑性应变达到临界阈值时，由于可见裂纹的形成，材料发生分离。值得注意的是，此过程中没有材料点被删除。
$$
\epsilon_{\text{plastic}} = \phi \cdot \left( \frac{\sigma - \sigma_f}{\sigma_2 - \sigma_1} \right)^n
$$
其中n是指数。
本研究使用牛顿流体模型来模拟流体介质。牛顿流体抵抗压缩，但不抵抗拉伸或剪切。流体的本构方程为式（8）（Janna，2020）：
$$
\tau = \mu \cdot \left( \frac{\partial \varepsilon}{\partial \varepsilon_x} \right)
$$
其中τ是粘性剪切应力（MPa）。

2.2. 材料点方法（MPM）
本研究使用MPM数值求解控制方程，并采用对流粒子域插值（CPDI）方法来处理元素间误差和数值裂纹（Sadeghirad等人，2011）。MPM将物体离散为材料点，并使用线性函数的背景网格。运动方程在欧拉背景网格上求解，拉格朗日粒子用于跟踪物理状态（Ma等人，2009）。MPM的原理和数值计算步骤由Sulsky等人（1994）描述。图1提供了所提模型的数值解示意图。
下载：下载高分辨率图像（685KB）
下载：下载完整尺寸图像
图1. 射流冲击模型的材料点离散化和CPDI方法。r1和r2是配置向量；p是角点的位置。
使用单值映射函数描述MPM元素与背景网格之间的动量交换，确保对象之间不发生穿透。因此，使用平均动量计算粘性接触：
$$
\vec{m} = \sum_{i = 1}^{n} \vec{m}_i \cdot \vec{n}
$$
其中mI,1和mI,2表示背景网格节点处的固体和流体质量（kg）；vI,1和vI,2表示背景网格节点处的固体和流体速度（m/s）；n是时间步长。
基于 control equations 和数值解原理，使用C++编程语言创建了程序以提高计算效率。我们利用Nvidia的CUDA库实现了GPU并行计算（Liu等人，2021）。计算代码用于模拟磨料和水的流动。与SPH和DEM不同，MPM不需要复杂的邻居搜索，显著提高了模型效率。在动态过程的模拟中观察到了以下优势：（1）可以使用大量更多的计算元素，从而充分表征复杂的异质结构和断裂模式，提高模型精度；（2）该模型能够对低粘度流体（如水）进行离散化，提高液体与固体相之间相互作用的建模精度；（3）计算时间显著增加（增加了10个数量级），反映了冲击断裂的不同阶段。MPM方法能够区分材料界面（由点定义）和计算网格（背景网格），消除了使用元素边界来划分复杂异质材料的需要，这是传统基于网格的方法所必需的。离散点包含了物理性质的信息，因此MPM在模拟异质材料和极端变形方面具有显著优势。

3. 模型与验证
3.1. 异质砂岩模型的构建
为了表征AWJ对砂岩矿物造成的损伤，需要一个真实的异质砂岩目标模型。这两个砂岩样本取自中国山西省大同市塔山煤矿8101工作面。#1砂岩的密度为2.37 g/cm3，#2砂岩的密度为2.68 g/cm3。建立异质模型的步骤如下：使用线切割技术将砂岩切割成长度25毫米、宽度25毫米、高度25毫米的尺寸，并将其表面打磨光滑。砂岩样本如图2所示。

在已发表的研究中，我们分析了矿物组成与岩石尺寸之间的关系，证实当边缘长度超过4毫米时，代表体积元素（RVE）可以达到误差收敛（Zhou等人，2025年）。因此，在本研究中，我们使用能量分散光谱（EDS）和背散射电子（BSE）成像技术对样本中8毫米×8毫米区域的矿物进行了定量分析（Liu等人，2024年）。TIMA3系统库用于统计分析、比较和拟合测试点数据（Tang等人，2022年）。该系统采用晶格数据采集方式，每个点记录1200次X射线计数。BSE图像分析的像素间距为3微米，EDS的像素间距为9微米。由此获得了矿物类型、矿物丰度和元素分布的信息。两种砂岩表面的矿物组成和分布如图3所示。

在#1砂岩中识别出六种主要矿物成分，占总量的96.1%：石英占40.7%，钠长石占32.5%，正长石占11.9%，孔隙和其他未鉴定的次要矿物占3.9%。在#2砂岩中识别出三种主要矿物成分，占总量的95.6%：石英占64.9%，白云母占22.7%，孔隙和其他未鉴定的次要矿物占4.4%。我们获取了矿物分布图像中像素的RGB值，并为不同类别分配了不同的材料属性和层间界面，以建立矿物分布的异质模型。上述建模方法的详细信息已在之前的研究中描述（Zhou等人，2025年）。

3.2. 数值验证
我们确定了并验证了构成岩石的矿物的力学参数。我们利用μTS介观尺度单轴压缩实验和MPM单轴压缩模拟来验证这些矿物力学参数及数值模型的可靠性。图2中所示的砂岩样本被重新切割成8毫米×8毫米×4毫米的尺寸。在四川大学教育部深部地下科学与工程重点实验室中使用μTS介观尺度测试系统进行了单轴压缩实验（图4a）。该实验系统能够对小型岩石样本进行高精度力学性能测试（位移分辨率为5纳米）。单轴压缩试验以1微米/秒的加载速率进行。试验后，将断裂的样本用透明环氧树脂固化，并使用微纳CT扫描（图4b）以获取其断裂形态。

我们根据以往的研究和文献（Du等人，2025年；Liu等人，2025年；Xi等人，2025年）获得了矿物的力学参数。构成岩石的矿物具有很高的机械强度，但岩石的强度受矿物之间弱界面的影响。然而，获取和准确表征矿物之间的界面参数极其困难。因此，我们根据微观力学实验结果调整了界面力学参数，以获得与实验数据一致的单轴压缩强度。矿物的力学参数和界面参数列在表1中。

表1. 砂岩模型中矿物和界面的力学参数
| 微粒 | 杨氏模量（GPa） | 泊松比 | 内聚力（MPa） | 摩擦角（°） |
|-----------|--------------|------------|--------------|---------------------|
| 石英 | 104 | 0.09 | 448 | 45 |
| 白云母 | 51 | 10.2 | 291 | 238 |
| 霉云母 | 75 | 0.26 | 226 | 38 |
| 褐帘石 | 36.7 | 0.32 | 043 | |
| 方解石 | 73 | 0.31 | 484 | |
| 正长石 | 79 | 0.29 | 346 | 22 |
| #1界面 | 10.3 | 0.25 | 152 | 25 |
| #2界面 | 21.5 | 0.25 | 363 | 30 |

我们使用表1中的参数和MPM，在与μTS实验相同的条件下对异质砂岩模型进行了单轴压缩模拟。图5显示了单轴压缩实验和模拟结果的失效模式对比。与实验类似，模拟中的样本主要经历了剪切失效，伴有剪切断裂和少量拉伸裂纹。数值模拟表明，裂纹形成是由于材料点的损伤、软化和分离。我们没有使用断裂算法或移除材料点。通过比较损伤模式，验证了数值模型描述失效行为的能力。

杨氏模量和单轴压缩强度被用作岩石宏观力学性能的指标。如图6所示，数值模拟与实验结果在峰值应力点和弹性阶段的斜率方面表现出高度一致性。这一发现证实了所提出的矿物尺度力学模型能够准确模拟岩石的宏观性能。关于失效行为的表征和力学性能的表示，异质砂岩的模拟结果与实验数据非常吻合。

4. 结果与讨论
4.1. 应力载荷特性
4.1.1. 应力载荷响应
材料和结构的异质性显著影响应力分布。图8显示了在不同20 MPa原位应力下不同岩石样本的应力分布。#1和#2砂岩的应力分布表现出明显的异质性。由于长石和石英等矿物的强力学性能，它们的应力响应较高。此外，岩石中的缺陷（孔隙、弱平面）也会显著影响应力分布，因为在加载过程中应力通常会在缺陷附近集中（Zhang等人，2024年）。

WJ对岩石的冲击包括水锤压力阶段和停滞压力阶段（Heymann，1969年）。当高压水流撞击固体表面时，在接触点会产生冲击波。我们通过时间步长细化并放大了这一中间过程。图9和图10展示了AWJ冲击下#1和#2砂岩的应力波传播过程，包括四个阶段：（a）水锤应力的产生和传播；（b）应力波的反射；（c）应力波的衰减；（d）停滞压力的集中。

岩石在水射流冲击下的失效是由于动态应力波和准静态停滞压力的共同作用。这一过程分为两个阶段：（1）动态应力波阶段。初始冲击产生了显著的水锤压力，导致应力波（纵波和剪切波），#1砂岩的峰值应力为157 MPa，#2砂岩的峰值应力为189 MPa，持续约3微秒。这些波通过岩石传播，并由于边界反射形成了复杂的叠加场，导致次级峰值应力分别为75.7 MPa（#1砂岩）和86.9 MPa（#2砂岩），持续时间为2微秒。这一阶段的主要失效模式包括靠近冲击区的剪切失效，主要由应力波的剪切分量引发；以及波反射、叠加和收敛区域（如边界）处的拉伸失效。（2）停滞压力阶段。初始应力波能量衰减后，来自水射流的持续应力显著降低，峰值应力分别为30.2 MPa（#1砂岩）和42.3 MPa（#2砂岩）。射流的停滞压力影响了冲击区域，导致了塑性失效。

在应力传输过程中，当磨料颗粒接触岩石时观察到了显著的应力集中（图10b），#1砂岩的峰值应力为116 MPa，#2砂岩的峰值应力为171 MPa。磨料冲击所产生的应力集中远高于水冲击产生的停滞压力，主要导致了剪切失效。由于计算时间步长取决于材料属性，强度更高的#2砂岩表现出更长的时间步长。在相同的冲击条件下，#2号砂岩在冲击区域表现出更显著的应力集中（见图10b和c），并且在应力传递和耗散过程中产生的拉伸损伤较少。应力波传递的基本机制是通过连续（弹性）介质传递能量。当材料发生塑性变形时，能量会被耗散。这可以归因于两个主要因素：（1）材料进入塑性阶段时消耗能量；（2）塑性应变会产生裂纹或不连续区域，应力波可能会在这些区域反射或传递。这些原因解释了在静压阶段应力波传递明显低于初始水锤压力阶段的现象。连续冲击会产生塑性区域，塑性应变软化会降低材料的机械强度。因此，材料在应力作用下吸收更多能量，形成受损和松动的区域。这些区域吸收冲击能量，降低了传递效率。

4.1.3 应力状态的演变
我们分析了整个周期内岩石元素的应力状态和损伤状态，以了解AWJ破碎岩石的过程，并揭示AWJ冲击岩石的损伤机制。在#1号和#2号砂岩的冲击表面的径向和周向选择了八个元素，并跟踪了它们的应力状态和损伤状态。原位应力为10 MPa，这两个方向的元素覆盖了目标材料的所有部分。图11显示了#1号和#2号砂岩中选定的岩石元素及其矿物。

我们分析了元素的最大主应力σ1、最小主应力σ3、冯·米塞斯应力σvm、剪切损伤程度DS和拉伸损伤程度DT。剪切损伤程度和拉伸损伤程度分别由剪切塑性应变和拉伸塑性应变的累积产生。图12和图13分别展示了#1号砂岩径向方向上元素#1-1至#1-4和#1-5至#1-8的应力状态和损伤状态。岩石元素在受到AWJ冲击时经历了四个阶段：（1）弹性阶段，（2）塑性应变累积阶段（塑性阶段），（3）损伤阶段，（4）剥落阶段。

图12显示了#1号砂岩径向方向上元素#1-1至#1-4的应力状态和损伤状态。
图13显示了#1号砂岩径向方向上元素#1-5至#1-8的应力状态和损伤状态。

在弹性阶段，径向方向上元素的应力状态是稳定的，没有剪切和拉伸塑性损伤。随着岩石破碎深度的增加，元素上的应力也随之增加。随后，元素进入塑性阶段，在此期间塑性应变开始累积，表明矿物达到了其拉伸/剪切屈服强度。塑性阶段是由磨料和水的共同作用引起的。元素的塑性应变通过塑性流动缓慢累积，直到达到临界值，此时元素进入损伤阶段。当达到岩石破碎深度时，由于磨料冲击产生的应力集中，元素迅速受损。进入损伤阶段后，元素能够承受压缩应力和剪切应力，并形成了松动或受损的区域。在AWJ的持续冲击下，这些区域发生剥落，并被侵蚀坑的回流带走。

基于网格的方法（FEM、FDEM）不太适合模拟损伤和剥落阶段，因为当元素达到临界塑性应变时会被删除，以防止计算错误和终止。

径向跟踪的8个元素的损伤演变表明，剪切损伤的累积始终高于拉伸损伤。这一发现表明，矿物元素的破坏主要是由于剪切失效，伴随着拉伸失效的累积，无论应力是由直接水流还是磨料冲击引起的。此外，并非所有岩石元素都经历了所有四个阶段。表面的元素在受到冲击后进入了塑性阶段，例如元素#1-1。位于岩石中间或底部的元素由于岩石破碎深度不足而保持在弹性阶段，因此应力载荷没有达到屈服强度，例如元素#1-8。另一方面，有些元素在受损后没有发生剥落，例如元素#1-7。

图14和图15分别展示了#1号砂岩周向方向上元素#2-1至#2-4和#2-5至#2-8的应力状态和损伤状态。周向元素位于喷射轴的两端，两个元素的位置几乎对称，例如元素#2-1和#2-8。与径向元素类似，周向岩石元素也经历了相同的四个阶段。然而，随着岩石的侵蚀，周向元素位于侵蚀坑的不同位置。除了#2-4和#2-5之外的元素没有受到水流和磨料的直接冲击，主要受到回流和侧向应力波的影响。因此，在弹性阶段出现了应力波动，不同位置的矿物元素表现出不同的损伤模式。在侵蚀坑和坑底，剪切失效占主导地位；而在坑边缘，观察到了拉伸失效，例如元素#2-1和#2-8。

图16和图17分别展示了#2号砂岩径向方向上元素#1-1至#1-4和#1-5至#1-8的应力状态和损伤状态。#2号砂岩的矿物组成包括石英、白云母和长石，其中石英占比为64.9%。因此，#2号砂岩的强度显著高于#1号砂岩。冯·米塞斯等效应应力和剪切塑性应变的累积速率也高于#1号砂岩。如4.1.2节所解释的，应力波的传递是弹性区域能量传递的主要原因。较高的弹性模量导致更高的传递效率，使矿物在弹性阶段的应力更高。

图18和图19分别展示了#2号砂岩周向方向上元素#2-1至#2-4和#2-5至#2-8的应力状态和损伤状态。这些元素没有经历显著的拉伸失效，主要表现出剪切失效，并显示出应力集中。高应力集中是由于磨料失效造成的，表明回流的冲击没有达到矿物元素的拉伸屈服强度。此外，侵蚀坑壁的矿物受到回流、剪切波和磨料作用的影响。水流冲击导致剪切失效并产生了水楔形拉伸失效，而磨料从各个角度冲击岩石，导致剪切失效。

图20展示了AWJ对异质砂岩的侵蚀过程。磨料和水在挖掘坑阶段影响岩石，回流与水平面形成一定角度，不断侵蚀坑缘。随着冲击时间的增加，岩石裂缝深度增加，回流几乎变为垂直。随着回流方向的变化，磨料和水的作用从坑缘转移到坑壁，使侵蚀坑宽度增加。AWJ将岩石元素磨损并剥落。

分析表明，受到AWJ冲击的矿物元素遭受了损伤，形成了受损或松动的区域。图21显示了不同原位应力水平下#1号砂岩的损伤形态。损伤区域是从AWJ到塑性区域的过渡区，覆盖了整个冲击表面。随着冲击时间的增加，损伤区域的厚度增加。此外，通常在界面和周围矿物晶体上出现明显的水楔形裂纹，形成晶间断裂。损伤区域的扩展方向与岩石的断裂方向一致。在4000 μs时，由于砂岩尺寸较小，横截面应力较高，且在0原位应力下发生了显著的横向拉伸断裂（I-e）。原位应力的增加提高了矿物强度，减少了裂纹，并降低了损伤区域的厚度。#1砂岩在不同侵蚀时间下的损伤区演变：(a) 冲击时间 = 808 μs；(b) 冲击时间 = 1616 μs；(c) 冲击时间 = 2424 μs；(d) 冲击时间 = 3232 μs；(e) 冲击时间 = 4040 μs；(I) 原位应力 = 0 MPa；(II) 原位应力 = 10 MPa；(III) 原位应力 = 20 MPa。图22展示了在不同侵蚀时间和不同原位应力水平下#2砂岩的损伤形态。与#1砂岩类似，#2砂岩也出现了水楔裂纹和拉伸破坏现象。两种砂岩在水楔损伤中都表现出明显的强度减弱特征，裂纹沿矿物界面一致发展或穿透冲击表面最弱的造岩矿物。随着原位应力的增加，裂纹数量和损伤区的厚度略有减小。由于#2砂岩的强度更高，原位应力对岩石和矿物的强化效果不如#1砂岩显著。在坑缘或坑壁上，整个矿物或多个矿物被剥离。这种损伤模式是在拉伸应力作用下的脆性断裂，这解释了坑缘和坑壁的扩大。此外，由于矿物的临界拉伸塑性应变较低，#2砂岩更容易发生拉伸破坏，因此随着冲击时间的增加，拉伸裂纹更多，坑缘也变得更大。通过观察损伤模式来分析AWJ对砂岩的影响。

图22. #2砂岩在不同侵蚀时间下的损伤区演变：(a) 冲击时间 = 840 μs；(b) 冲击时间 = 1680 μs；(c) 冲击时间 = 2520 μs；(d) 冲击时间 = 3360 μs；(e) 冲击时间 = 4200 μs；(I) 原位应力 = 0 MPa；(II) 原位应力 = 10 MPa；(III) 原位应力 = 20 MPa。

4.2.3 损伤模式和均匀性
受到AWJ冲击的矿物元素经历了塑性阶段，在此阶段塑性应变的积累降低了它们的机械强度。随后，它们进入损伤阶段，被喷射流迅速剥离并通过回流从侵蚀坑中移除。图23展示了在10 MPa原位应力下不同侵蚀时间内#1砂岩的塑性区形态。I表示剪切损伤，II表示拉伸损伤。塑性区是损伤区和弹性区之间的过渡区，代表了裂纹和损伤的增长方向。#1砂岩的剪切损伤发生在冲击表面附近，其扩展方向与应力集中方向相同。相比之下，拉伸损伤主要发生在边界附近。如4.1.2节所述，这些位置是应力波的反射和重叠区域。

图24展示了在10 MPa原位应力下不同侵蚀时间内#2砂岩的塑性区形态。与#1砂岩相比，#2砂岩的主裂纹生长方向更为明显，塑性区较小，包括剪切和拉伸损伤。然而，由于缺陷（孔隙、弱层）的生长，观察到了更多的微孔损伤。如4.1.1节所述，岩石中的缺陷在载荷作用下会导致应力集中并加速其扩展。由于#2砂岩的强度较高，弹性波在传播过程中大部分被缺陷吸收，从而形成了孔洞。相比之下，#1砂岩中的弹性波被低强度矿物和界面吸收，增加了拉伸塑性应变。

图25显示了原位应力对非均质砂岩中受损元素数量和侵蚀坑面积的影响。受损元素数量是剥离元素和损伤区内元素的总和。#1砂岩对原位应力更敏感，在应力梯度上受损元素数量和侵蚀坑面积的变化更大。这一现象主要归因于两种机制：(1) 围压抑制了拉伸破坏。由于#1砂岩的矿物强度和界面参数较低，在围压增加时缺陷较少。(2) #2砂岩的密度比#1砂岩高13%，其孔隙率更低，缺陷也更少，导致矿物间的压实程度较低。因此，强度的提升不如#1砂岩明显。随着侵蚀时间的增加，单位时间内受损元素数量和侵蚀坑面积减小。

图25. 原位应力对受损元素数量和侵蚀坑面积的影响：(a) #1砂岩中的受损元素数量；(b) #2砂岩中的受损元素数量；(c) #1砂岩中的侵蚀坑面积；(d) #2砂岩中的侵蚀坑面积。

4.3 侵蚀坑的宏观-微观形态
我们使用AWJ和WJ对#1和#2砂岩进行了侵蚀实验，以研究磨料和高压水对造岩矿物的破坏行为。实验装置在第3.3节中有描述。侵蚀岩石样品的尺寸为25毫米（长度）× 25毫米（宽度）× 50毫米（高度）。喷射压力为100 MPa，磨料浓度约为6%，与数值模拟结果一致。所有实验的侵蚀时间为4秒，无原位应力加载。

图26显示了WJ和AWJ侵蚀实验后#1砂岩的坑形态。左侧显示横截面，右侧显示图像重建结果。WJ侵蚀后#1砂岩的坑深度为5.94毫米，坑宽度为1.96毫米（图26a）。AWJ侵蚀后坑深度为14.55毫米，坑宽度为2.35毫米（图26b）。模拟的腔体宽度为2.5毫米，误差为6%。AWJ侵蚀后坑深度和宽度显著增加，表明磨料由于回流作用影响了坑壁。WJ引起的坑具有更多分支裂纹，坑壁上存在更多的水楔损伤，而AWJ作用的坑壁相对光滑。

图27展示了WJ和AWJ侵蚀实验中#2砂岩的坑形态。WJ侵蚀后#2砂岩的坑深度为4.87毫米，坑宽度为4.08毫米（图27a）。AWJ侵蚀后坑深度为10.84毫米，坑宽度为3.05毫米（图27b）。模拟的腔体宽度为3.17毫米，误差为3.7%。实验和数值结果证明了在拉伸载荷下发生了广泛的脆性断裂，导致#2砂岩的坑宽度明显大于#1砂岩。这种现象间接验证了模型的准确性。生成的裂纹几乎穿透了整个岩石样品，坑周围的区域成块剥落，导致坑尺寸更大。在第4.2.3节的模拟结果中也观察到了类似现象，并给出了解释，这是脆性岩石中拉伸破坏的特征。我们分离出由喷射流在#2砂岩中形成的主坑并进行了重建（右侧）。与#1砂岩不同，#2砂岩中没有发生分支裂纹。WJ和AWJ侵蚀后坑壁相对光滑，这与石英、云母和长石的断裂或节理形态密切相关。

图28展示了WJ和AWJ穿孔在#1砂岩中部的微观特征。损伤特征位于喷射穿孔的凹槽内。WJ侵蚀导致大量矿物松动和水楔裂纹，这些是晶界断裂。如4.2.2节所述，水楔显著削弱了材料并破坏了矿物间的弱界面，导致矿物在冲击过程中脱落。AWJ侵蚀引起的损伤区域没有观察到松动的矿物，但存在水楔裂纹。产生了大量粗糙表面，粗糙表面上有许多撞击坑和犁沟痕迹。犁沟是由磨料以低角度撞击矿物表面造成的切割效应。两种侵蚀方法都在节理面上产生了晶界断裂，其中石英表现为弯曲断裂，长石表现为条纹状或层状断层结构。

图29展示了WJ和AWJ穿孔在#1砂岩底部的微观特征。WJ导致大量矿物松动和水楔裂纹。AWJ造成的犁沟痕迹较少，撞击坑是主要特征。水流中的磨料接触到了坑底或坑壁的松散小矿物，使其脱落并产生粗糙表面。大颗粒矿物（如石英）上形成了不规则的晶界断裂表面。撞击坑或切割坑是由于后续的磨料冲击形成的。整个过程中都发生了水楔损伤，对冲击表面产生了拉伸应力，增强了磨料的切割和冲击破碎效果。

图30展示了WJ和AWJ穿孔在#2砂岩中部的微观特征。由于拉伸断裂导致的广泛矿物剥落，#2号砂岩在WJ侵蚀过程中没有显示出规则的侵蚀孔。此外，#2号砂岩密度非常高，在侵蚀坑的中间观察到了少量松散的矿物。WJ仅通过晶间裂纹去除矿物，导致岩石断裂并出现大面积剥落。这种效应在细粒矿物之间形成了粗糙的断裂表面和晶间裂纹，而在粗粒矿物之间形成了平滑的节理表面和穿晶裂纹。磨料增加了节理和断裂表面的粗糙度。粗糙的表面上出现了冲击坑、切割坑以及由磨料冲击形成的微裂纹。

图30. #2号砂岩中水射流穿孔区域的显微特征：(I) 水射流，(II) AWJ。(a) 50倍放大，(b) 150倍放大。

图31显示了#2号砂岩中WJ和AWJ穿孔底部的显微特征。WJ形成了侵蚀坑，#2号砂岩的底部出现了穿透岩石样品的拉伸断裂，导致矿物松散、节理表面光滑且呈穿晶状，而晶间裂纹表面粗糙。暴露于AWJ的坑底也显示出粗糙的表面以及冲击和切割坑。总之，磨料的冲击和切割作用有效地去除了由水楔效应引起的松散矿物和不规则的断裂表面，显著提高了破岩效率。

4.4. AWJ的破岩机制
AWJ去除材料的机制包括切割、推开和挤压断裂（Desale等人，2006年；Tarodiya和Levy，2021年）。材料磨损研究对此机制提供了很好的理解。此外，不同颗粒的冲击角度可能会产生不同的破坏效果（Al-Bukhaiti等人，2017年）。本研究使用数值模拟和实验结果来阐明AWJ影响砂岩的侵蚀机制，同时考虑了高压水射流与磨料对造岩矿物之间的耦合效应。

图32展示了AWJ在初始阶段和坑形成阶段影响砂岩的机制。数值模拟和样品的显微形态表明，冲击表面由松散的矿物和断裂区组成。这些区域通常无法承受高应力，容易发生剥落，因此是损伤区。损伤区和完整岩石之间存在一个过渡区，在这个区域中，虽然应力很高，但材料不会松散或断裂。这一区域可以被视为塑性区。完整岩石区域可以定义为弹性区。

图32. AWJ对砂岩的侵蚀损伤机制：(a) 裂纹起始；(b) 坑形成。
在裂纹起始阶段(a)，水锤压力和磨料的冲击迅速在冲击区产生了塑性状态。磨料引起的应力集中显著高于水压造成的应力。在这一阶段，磨料主要在高角度下冲击，主要导致冲击引起的破碎。塑性区域内的矿物受到损坏并被带走，伴随着应力波的传播和反射。这一阶段是侵蚀坑形成的关键时期。

随着岩石破碎深度的增加，进入了坑形成阶段(b)。磨料与冲击表面之间的冲击角度发生变化，射流的水锤效应转变为水楔效应。高压水产生的水楔效应对侵蚀坑施加了拉应力，导致岩石变弱。它破坏了矿物之间的胶结作用，增加了岩石的缺陷，使矿物松散并生成水楔裂纹。脆性岩石经历了明显的矿物剥落。由于冲击 angle 的变化，磨料在侵蚀坑处引起了冲击破碎和切割。在坑底，射流的直接冲击导致了冲击破碎，而坑壁则因回流效应而发生切割。

需要强调的是，磨料的切割和冲击破碎作用都属于剪切损伤，两者都是由于水楔效应引起的。因此，AWJ的破岩机制包括磨料在水楔效应的拉应力作用下对造岩矿物的剪切破坏。此外，磨料有效地去除了由水楔效应引起的松散矿物和不规则的断裂表面。磨料与高压水对岩石的耦合效应是高压AWJ高效破碎硬岩的原因。

在工程实践中，可以根据储层岩石类型调整参数以优化AWJ的破岩效果。关于岩石侵蚀机制的研究成果为开发理论模型提供了指导。这些模型可用于预测应用AWJ的最佳深度并优化其参数。根据周围岩石特性可以预测岩石破碎深度，为地下工程设计及AWJ操作的成本优化提供关键支持。

5. 结论
我们进行了AWJ实验和数值模拟。在数值模拟中构建了两个异质砂岩目标体。AWJ侵蚀过程基于相同的条件下的MPM进行模拟。实验后的侵蚀形态利用微纳CT进行了重建，并进行了SEM观察。模拟和实验结果用于分析AWJ对矿物的损害情况，并解释了侵蚀机制。得出以下主要结论：
(1) 材料和结构的异质性显著影响了应力分布。长石和石英具有更高的强度和更明显的应力响应，在损伤过程中通常伴随着更高的应力或更长的塑性阶段。岩石中的缺陷（如孔隙和弱面）在加载时成为应力集中区。
(2) 不能忽视磨料的体积效应。磨料的冲击产生了显著的应力集中，超过了水锤应力，导致剪切破坏。造岩矿物的冲击破碎是AWJ测试中主要的坑挖掘方式。能量传递发生在弹性区域内的应力波传播过程中，当材料进入塑性损伤状态时能量被耗散。应力波主要在冲击区附近引起剪切破坏，而在耗散区发生拉伸破坏。
(3) 高压水产生的水楔效应导致了由应力驱动的破坏。裂纹总是沿着矿物之间的界面或通过冲击表面上强度较低的矿物传播，使单个或多个矿物脱离，从而在拉应力下发生脆性断裂。这种破坏模式在石英含量较高的脆性岩石中更为明显，导致在模拟和实验中Sandstone #2的坑宽度始终大于Sandstone #1。
(4) 持续的水楔效应在冲击表面上产生了拉应力，使磨料能够在拉应力下有效去除松散的矿物和不规则的断裂表面。磨料与高压水的耦合效应使得高压AWJ能够高效破碎硬岩。

研究展望包括：
(1) 三维（3D）异质建模。利用训练有素的深度学习算法处理矿物分布图像，将现有的2D模型扩展到3D框架中。使用高分辨率CT数据进行3D模型重建，从而准确模拟复杂岩层中的AWJ侵蚀过程。
(2) 基于机制的深度预测建模。开发基于岩石侵蚀机制的理论深度预测模型，而不是依赖算法或经验拟合。

作者贡献声明：
Lei Zhou：验证、监督、软件开发。
Yiyu Lu：监督、资金获取、概念构思。
Yanmin Zhou：软件开发、方法论、正式分析、数据管理。
Binwei Xia：资源准备、资金获取、正式分析。
Sisong Zhang：写作——审稿与编辑、初稿撰写、数据管理、概念构思。

未引用的参考文献：
Zhou等人，2025年。