基于SPH-FEM耦合算法的后混合磨料水射流冲击破岩数值模拟研究

2021-09-27米建宇李树清王荣荣

振动与冲击 2021年16期

米建宇，黄飞,2,3，李树清,2,3，王荣荣，李丹

(1.湖南科技大学资源环境与安全工程学院，湖南湘潭 411201；2.湖南科技大学南方煤矿顶板及煤与瓦斯突出灾害预防控制应急管理部重点实验室，湖南湘潭 411201；3.湖南科技大学南方煤矿瓦斯与顶板灾害预防控制安全生产重点实验室，湖南湘潭 411201)

磨料水射流是一种将高速水射流与磨料相互混合而成的一种固、液两相射流，目前被广泛运用于破碎岩石、切割金属等方面[1]。根据磨料掺合方式的不同，磨料水射流可分为前混合式和后混合式两种。后混合磨料射流的磨料是在射流形成后通过独立软管在喷嘴内部完成混合，有效避免了在前混合磨料射流中高压水管堵塞和高压水管管路磨损严重等问题。后混合磨料是在高速水射流的推动牵引下获得一定的加速，鉴于磨料与水的物理差异性，磨料粒子很难与水充分混合。为了获得较高速度的冲击磨料，通常采用提高系统泵压、改变喷嘴结构等手段。然而，目前对于后混合磨料在喷嘴内的加速情况以及对破岩效果的影响规律方面的研究仍较为棘手。因此，开展后混合磨料水射流的磨料加速及其对破岩影响的数值模拟研究，对提高后混合磨料射流工作效率具有重要的工程意义。

磨料水射流破岩作为一个涉及到诸多因素的非线性冲击动力学问题，其具有的瞬时强值动荷载、大变形以及高应变率等特点使得磨料水射流破岩难以从传统理论以及实验手段研究其冲蚀损伤及破岩机理。然而随着计算机技术和工程商用软件的发展，传统理论和实验手段难以解决的问题可以通过相关软件的数值模拟对上述问题进行研究分析，这一方法已被大量学者运用并解决了诸多通过传统理论手段所不能解决的问题。林晓东等[2]通过SPH-FEM(smoothed particle hydrodynamics with finite element method)耦合方法对前混合磨料水射流中磨料粒子加速进行了研究，发现了磨料粒子在喷嘴内部不同阶段的加速特征及喷嘴结构对加速效果的影响规律。Li等[3]通过SPH-FEM耦合方法，从射流角度和施加围压方面出发，研究了射流倾角和围压对岩石破坏的影响。Wu等[4]通过对LS-DYNA程序的二次开发，建立了SPH-FEM/DEM数字高程模型模型，采用Voronoi细分技术研究了了在微观结构和细观力学性质下水射流破碎岩石的性能。表明了水射流冲击岩石的破碎性能受颗粒大小和不规则性、延性、微观强度、微观参数的非均匀性等因素的影响较大，而受接触刚度的影响较小。Ma等[5]通过研究基于光滑粒子流体动力学(smoothed particle hydrodynamics，SPH)法磨料水射流破碎岩石的动态过程，揭示了提高磨料浓度和射流速度可以大大提高磨料水射流破岩效果。Jiang等[6-7]通过SPH/FEA(smoothed particle hydrodynamics/finite element analysis)方法建立了水射流破岩的数值模型，并发现了岩石的裂纹和损伤的产生机制。刘佳亮等[8-9]通过对高压水射流冲击混凝土的研究，揭示了混凝土的致裂机理及裂纹扩展规律。文献[10-11]通过对水射流中喷嘴流场的研究，揭示了喷嘴最佳收敛角度并对喷嘴进行了优化设计。

以往研究主要对前混合磨料射流加速以及其破岩规律开展了研究，对于后混合磨料射流破岩的数值模拟通常将磨料和水的速度设为相同，与实际情况不符。本文拟采用SPH-FEM耦合算法并限定SPH领域搜索范围的方法来模拟后混合磨料射流的形成，进而研究射流速度、磨料浓度与围压大小等因素对后混合磨料水射流破碎岩石的影响规律。

1 基本算法描述

1.1 SPH基本原理

SPH法是无网格法的一种，其核心理论为差值理论[12-14]。对问题域采用离散化的粒子进行描述，有效避免了在非线性计算中单元畸变所导致的计算终止。其核近似函数为

(1)

式中：f(x)为核近似函数；W(x-x′,h)为光滑核函数；x-x′为粒子间距；h为粒子的光滑长度；Ω为点x处光滑核函数的支持域。由粒子近似可将连续SPH积分形式写成离散化粒子近似式

(2)

Wij=W(xi-xj,h)

(3)

式中：f(xi)为粒子i处的近似函数值；mj为粒子j的质量(j=1,2,…,N)；N为在粒子j的支持域内粒子的总量；pj为粒子j的密度。

式(2)说明粒子i在任一函数值都可通过光滑核函数对其紧支域内所有粒子相对应的函数值进行加权平均近似。

1.2 SPH邻域搜索

虽然无网格粒子法可以解决传统的有限元网格大变形问题，但是其计算精度与计算效率较低等问题仍不容忽视[15-16]。当SPH磨料粒子在喷嘴内运动至收敛段、直线段时，截面粒子数量增多、粒子间距缩小，如图1(a)所示。这将导致SPH邻域搜索范围内的粒子数急剧增多，进而急剧增加领域搜索时间。随着柱塞的持续推进，SPH磨料粒子持续被压缩，整个算例的计算时间将会无限延长，最终导致计算失效。现有的SPH数值模拟技术难以从算法上来提高模型的计算效率，为了节省计算时间同时不影响SPH粒子在喷嘴内的相互作用，本文将采用图1(b)所示的方法对SPH粒子的邻域搜索范围进行限定，从而节省计算时间，提高计算精度。

图1 邻域搜索示意图Fig.1 Schematic diagram of neighborhood searching

2 计算模型建立

2.1 磨料水射流模型

为了更好的模拟真实工况，首先将水射流假定为一段圆柱进行几何建模，并在其中建立若干随机分布的直径为0.6 mm(30目)的磨料粒子。为防止水射流和磨料初始穿透，对水射流几何模型和磨料粒子几何模型进行布尔运算后分别进行有限元网格划分，最后将水射流模型网格转化为SPH粒子，保留磨料粒子的Lagrange有限元网格属性。其中：水采用NULL材料模型，并赋予其Mie-Grueisen状态方程；磨料粒子采用刚体模型。水与磨料的材料参数如表1与表2所示。其中：ρ0为材料密度；E为弹性模量；v为泊松比；φ为直径。Mie-Grueisen状态方程为

表1 水的本构模型参数Tab.1 The constitutive parameters for water

表2 磨料的模型参数Tab.2 Parameters of state for abrasive

(4)

式中：E为单位体积内能；C为vs-vp曲线截距；S1，S2，S3为vs-vp曲线斜率系数；γ0为Grueisen常数；a为一介体积修正量。

2.2 柱塞、喷嘴模型

将后混合磨料水射流考虑成通过柱塞以恒定的速度来推动纯水，然后由运动的水来推动磨料粒子并形成初步混合，最后经过收敛喷嘴加速从而喷射出高速磨料水射流。本文中假设柱塞与喷嘴均不发生形变，故不考虑它们的厚度变化。因此，本文采用普通钢制材料模型来模拟柱塞，而采用刚体材料模型来模拟喷嘴。柱塞与喷嘴均采用壳单元，并用Lagrange法进行有限元网格划分。

2.3 岩石本构模型

岩石采用HJC模型，该模型综合考虑了应变率效应、损伤演化效应、围压效应和压碎、压实效应影响，可以较好的描述混凝土、岩石材料在大变形、高应变率和高静水压力下的力学行为。该模型形式简单、参数意义明确，现已被广泛应用于冲击爆炸等强动载作用下混凝土、岩石类材料的分析中[17-18]。该模型强度采用规范化等效应力描述

σ*=[A(1-D)+BP*N][1+Clnε*]

(5)

D=∑[(Δεp+Δμp)/D1(p*+T*)D2]

(6)

式中：σ*=σ/fc；σ*≤SMAX,SMAX为岩石所能承受最大强度；P*=P/fc；ε*=ε/ε0，ε为应变率；ε0为参考应变率，ε0=1.0 s-1；T*=T/fc；σ为等效应力；P为单元内静压；T为材料最大拉伸强度；A，B，C，N，D1，D2为岩石材料常数；D为损伤度，0≤D≤1.0，且D1(p*+T*)D2≥EFMIN,EFMIN为岩石最小断裂应变；Δεp和Δμp分别为在一个积分步长内单元的等效塑性应变和塑型体积应变。在单元的变形中分为抗压和抗拉两种变形情况。fc为材料的抗压强度，其余的岩石本构模型参数如表3所示。

表3 岩石本构模型参数Tab.3 Parameters of state for rock

2.4 接触耦合处理

本模型中采用的独立单元较多，有采用SPH粒子模拟的纯水介质单元，还有采用Lagrange网格模拟的喷嘴、柱塞、磨料与岩石单元。为了实现柱塞推动水，水带动磨料粒子经过收敛段加速后喷射，最终作用于岩石进行破岩的全过程，本模型中对水与磨料、水与喷嘴、水与柱塞均采用点对面单向接触；涉及到磨料粒子与喷嘴及岩石的接触碰撞，由于磨料整体较小，也采用点对面单向接触；涉及到水和磨料整体对岩石的冲蚀，采用带有侵蚀算法点对面单向接触；磨料粒子之间的接触采用自接触定义。

2.5 几何模型及网格划分

根据以上建模方法建立的几何模型，如图2所示。图中：柱塞直径为14 mm并限制其自由度使其只在Z方向上运动；纯水段直径为14 mm，长度为50 mm；混合段直径为14 mm，长度为50 mm；收敛段长度为20 mm；直线段长度为10mm；喷嘴出口直径为3 mm；岩石为40 mm×40 mm×30 mm的长方体，根据林晓东和刘佳亮等的研究对其模型的设定，借鉴文献中所述边界条件及求解策略，对岩石周围四面与底面施加非反射边界，同时约束岩石的竖向移动。磨料粒子采用0粒、25粒、50粒、75粒、100粒5种浓度分别计算并分析对岩石的冲击损伤效果。

图2 后混合磨料水射流模型示意图 (mm)Fig.2 The model of post-mixed abrasive water jet (mm)

3 后混合磨料水射流加速与破岩数值计算结果与分析

3.1 磨料粒子加速分析

根据2.1节～2.5节建模方法建立模型后进行初步计算，获得水与磨料各阶段的加速过程示意图，如图 3 所示。根据图示，纯水在柱塞的推动下带动磨料向喷嘴收敛段流动，最终从喷嘴射出形成磨料水射流。由于本模型中所建立的磨料粒子并非单独的一粒，且在磨料加速过程中涉及到磨料与喷嘴的碰撞和磨料间自身的碰撞，为了尽可能完整的提取到磨料粒子的加速过程，选取其中碰撞较少的一颗磨料粒子并同时选中与磨料处于同一断面的SPH 粒子，提取两者的速度并且绘制速度-时间曲线，如图3所示。柱塞以一定速度作用于水，使得水在喷嘴混合腔内瞬间获得一个较高的速度，与此同时高速水射流与混合于其中的静止的磨料发生相互作用，强烈的速度差导致磨料在混合腔内发生剧烈的加速效果[19-20]；当磨料运动到喷嘴收敛段时，收敛段积聚的能量被部分释放，使得磨料在收敛段被继续加速；由于磨料粒子的质量较大，提速较为缓慢，在进入直线段后仍与水有较大的速度差，且直线段过流断面较小，使得磨料粒子继续加速并逐渐趋近于水，最终在喷嘴出口磨料的速度加速至纯水速度的 80%左右，并与水一起作用于岩石；当水与磨料分别撞击到岩石后，动能用于破岩能量消耗，各自的速度急剧减小。

图3 水与磨料的各阶段加速过程与速度变化规律Fig.3 The accelerating process of water and abrasive and their velocities’changing law

3.2 磨料水射流破碎岩石损伤规律

柱塞速度为300 m/s、磨料浓度为75粒的岩石冲蚀截面形状演化图，如图4(a)所示。其原理是：柱塞作用于水射流，给予水射流一定速度，水射流带动混合与其中的磨料经过收敛喷嘴从而喷射出高速磨料水射流并与岩石发生接触；首当其冲的岩石中心区域单元被迅速破坏失效并删除，形成初始孔径。随着冲蚀的进行，磨料水射流的冲击使得孔径逐渐加深，与此同时与底部接触的粒子反弹而后向周围飞溅对周围孔壁进行冲蚀并扩大孔径，由于反弹飞溅粒子能量较低，随着冲蚀的进行，孔径大小逐渐稳定并且深度不断增加，形成了如图4所示的“V”型剖面。图4(b)与图4 (c)为相同工况下磨料水射流冲击岩石的宏观破碎形态与细观CT扫描图像。通过对比可知：本论文的数值模拟现象与冲击破岩实验中的岩石破碎形态几乎一致，且与文献[21-23]中的结果相符。由此间接表明本文所建立的后混合磨料水射流冲击岩石模型恰当，可用于后续的后混合磨料水射流破岩的损伤效果研究。

图4 冲蚀坑对比图Fig.4 The comparison of erosion pit

3.3 磨料浓度对岩石破碎损伤的影响

为了探究磨料浓度对岩石损伤的影响，有限元法(finite element method，FEM)磨料采取0粒、25粒、50粒、75粒、100粒5种不同的浓度以柱塞300 m/s的速度分别进行计算，并提取在相同时间内不同磨料浓度的岩石损伤情况并绘制关于磨料浓度的破碎曲线图，如图5、图6所示。

图5 不同磨料浓度冲蚀坑截面演化图Fig.5 Evolution of erosion pits of different abrasive concentration

图6 破碎效果与磨料浓度的关系曲线Fig.6 The relationship of broken effect vs.abrasive concentration

在3.1节～3.2节的分析中，磨料在喷嘴中加速后以趋近于水的速度和水射流一起作用于岩石，鉴于磨料颗粒大，惯性大，动能高，在磨料冲击碰撞岩石后，产生的冲击压力与磨削力超过了岩石的抗压强度，使得撞击点附近岩石瞬间产生压剪破碎，冲击压力所引发的环向拉应力及应力波反射拉应力使得岩石进一步破裂[24]。磨料在撞击岩石后速度衰减并向四周反弹，逐渐形成冲蚀漏斗，随着冲蚀的进行，破碎深度逐渐增加。从图5可知：在没有磨料的纯水射流中，岩石的损伤程度较小，在加入25粒磨料后，冲蚀深度和宽度都有一个大幅度的提升，随着磨料的增多，第一时间内作用于岩石的磨料也随之增加，岩石的破坏深度也随之加深[25]。不过在图6中不难看出，当磨料在加入70粒以上后，破碎深度却几乎不再增加，而破碎宽度反而有一定减小。分析表明：磨料粒子的质量大、惯性大，在水带动其加速的过程中存在较大的能量交换。随着磨料的增加，水与磨料交换的能量增多，从而导致射流喷出时的总能量减弱，破岩的效果反而降低。因此可以推测，磨料的浓度需要控制在一个适宜的范围以内，以保证磨料射流喷出时拥有更大的能量。在本数值模型中，磨料浓度在70～80粒内最佳。

3.4 射流速度对岩石破碎损伤的影响

在了解了磨料浓度对岩石破碎损伤影响之后，采用75粒磨料浓度的模型依次进行柱塞速度为100 m/s，150 m/s，200 m/s，250 m/s，300 m/s，350 m/s的模拟，并获取在相同时间内不同速度下岩石损伤的效果并提取岩石的破碎深度及宽度并绘制关于速度的破碎情况曲线图，如图7、图8所示。

图7 不同速度冲击下冲蚀坑演化图Fig.7 Evolution of erosion pits under different velocity

根据文献[26-27]中的阐述，在水射流冲击岩石中，会产生水锤压力和滞止压力，其产生的水锤压力是导致岩石破坏的主要因素。在距离冲击中心的一定范围内，有效应力随径向距离的增大而减小，最终趋于稳定。图8所示的曲线表明：岩石的破碎深度与射流速度呈现近似线性增加的趋势。这是由于随着磨料水射流速度的增加，冲击的能量增大，岩石的破碎深度也随之增加。然而，岩石的破碎宽度随着射流速度增加有一定程度的增大，随后逐渐趋于稳定。分析表明：射流冲击压力随着射流速度的增加而增加，岩石受冲击点附近的应力也随之增加，破碎宽度有一定程度增加。继续增加射流速度，冲击中心的应力增加了，但是射流的直径并未有变化，中心处所产生的冲击应力影响的范围有限，破碎的宽度逐渐趋于稳定。由此推测，对于相同直径喷嘴的磨料水射流，单一增加射流速度工况下，岩石的破碎宽度存在一个最大值，本次数值模拟的最大破碎宽度约为13 mm。

图8 破碎效果与射流速度的关系曲线Fig.8 The relationship of broken effect vs.jet velocity

3.5 围压作用对岩石破碎损伤的影响

由于磨料水射流在应用于煤矿、石油钻探等工程领域时，岩石大都处于高围压状态，为了更加贴合真实情况，故在数值模拟计算时对岩石施加围压约束，通过模拟后磨料水射流冲击破碎高围压岩石的过程，分析磨料水射流冲击作用下高围压岩石的损伤演化过程，具有重要意义。为了更好的模拟真实工况，对岩石周围4个面施加压力以模拟在真实工况中有围压作用下的岩石破碎情况。依旧采用柱塞速度为300 m/s，磨料浓度为75粒的模型对岩石施加2 MPa，4 MPa，6 MPa，8 MPa，10 MPa的围压分别进行模拟计算，并绘制围压作用下岩石破碎情况曲线图，如图9与图10所示。

图9 不同围压作用下冲蚀演化图Fig.9 Erosion evolution map under different confining pressures

图10 破碎效果与围压关系曲线Fig.10 The relationship of broken effect vs.confining pressure

根据图示可知：随着围压的增加，岩石的破碎损伤效果出现减弱的趋势。其中，破碎深度随围压的增加呈现线性的减小趋势，而破碎宽度随围压先减小后趋于稳定；另外，随着围压的增加，岩石损伤破碎口的形态由“V”形向“U”形转变，这种变化趋势与文献[28]所得结果基本一致。经过分析，这种变化趋势可能是由多种因素影响所导致的。首先，在磨料水射流的冲击下，高速的磨料射流与岩石内部的孔隙流体之间形成巨大的压力差，在压力差的作用下流体向岩石内部渗透，从而引发“水楔作用”，该作用有利于岩石的进一步损伤破碎。在围压的作用下，磨料水射流的水楔作用会被减弱，且减弱的程度随围压的增加而加强，不利于岩石的损伤破碎。其次，根据岩石强度理论原理可知，岩石的力学强度会随着围压的增加而增大，从而不利于磨料水射流冲击下岩石的损伤发展。再次，根据文献[29-30]可知，围压会影响磨料水射流冲击下岩石内应力波的传播，从而减弱应力波对岩石的损伤作用，并增加岩石损伤所需的时间。