Numerical Simulation of Rock Fragmentation under the Impact Load of Water Jet

Toinvestigatetherockfragmentationanditsinfluencefactorsundertheimpactloadofwaterjet,anumericalmethodwhichcoupled finiteelementmethod(FEM)withsmoothedparticlehydrodynamics(SPH)wasadoptedtosimulatetherockfragmentation processbywaterjet.Linearandshockequationsofstatewereappliedtodescribethedynamiccharacteristicsofrockandwater, respectively,whilethemaximumprincipalstresscriterionwasusedfortherockfailuredetection.Thedynamicstressesatthe selectedelementcontainingpointsinrockarecomputedasafunctionoftimeundertheimpactloadofwaterjet.Theinfluences ofthefactorsofboundarycondition,impactvelocity,confiningpressure,andstructureplaneonrockdynamicfragmentationare discussed.


Introduction
Rock fragmentation by the water jet technology has been widely used in mining, petroleum drilling, civil construction, gas drainage, and cleaning operations [ -]. However, the rock fragmentation mechanism is still unclear because of the opacity and damage instantaneity of the rock. Nowadays, there are extensive arguments on the fragmentation patterns and processes [ ].
e rock fragmentation by water jet impact not only is a di cult problem in mechanics, but also presents a new challenge for physics [ ], whereas there are few investigations of the rock deformation and fragmentation by the water jet impact. Bowden and Brunton took note of the failure pattern on the brittle material surface under the high speed liquid impact and found that the discrete cracks formed regularly surrounding the impact area [ ]. Meanwhile, they studied the fracture of brittle material and n o t edth a ts h ea ra n dt ea roc cu rr edo nth er oc ks u rf a c ed u e to the high speed liquid, and the so-called water-hammer pressure can cause substantial damage to brittle material [ ]. Bourne et al. studied the fracture of brittle material (PMMA) under liquid jet impact using high-speed photography and schlieren visualization, and they found that the failure resulted from the interactions of di erent stress waves [ ]. Momber investigated the rock erosion due to the liquid jet impact and observed that the failure was mainly attributed to the lateral jetting [ ]. Ni et al. systematically studied the rock fragmentation under high-pressure water jet drilling through experimental, theoretical, and numerical method, respectively, and they considered that the stress wave played the leading role in rock fragmentation [ , ]. Li and Liao investigated the micromechanism of rock failure by water jet impact using scanning electron microscope and observed that the fracture was mainly caused by two types of transgranular fracture and shear dislocations [ ]. In recent years, many researchers adopted the numerical method to investigate the failureofrockorrock-likematerialsbywaterjetimpact,and their research works were mainly focused on the in uences of the impact velocity, diameter, stando distance, and incidence angle of water jet, as well as the con ning pressure and type of rock on the erosion depth, diameter, damage, and stress distributions of rock [ -]. However, there is less investigation on the rock fragmentation mechanism by the compression, shear, or tensile stress. Although many signicant results have been published, they are far from complete for the numerical study of rock fragmentation by water jet impact. e present paper is aiming at revealing the rock fragmentation mechanism and explaining the reasons for crushing zone formation, crack initiation, and propagation u n d e rt h ei m p a c tl o a do fw a t e rj e t .T oa c h i e v et h i sg o a l , a numerical method coupled FEM with SPH was adopted and the following topics were investigated: (a) the fracture process of rock sample under water jet impact; (b) the mechanism of crushing zone formation, crack initiation, and propagation; (c) the e ects of boundary condition, impact velocity, con ning pressure, and structure plane on the rock failure and crack initiation and propagation.

Numerical Model
. . Coupled SPH/FEM Method. As we all know, the large deformation of water exists in the rock fragmentation process b yw a t e rj e ti m p a c t . ec a l c u l a t i o nw o u l de a s i l yt e r m i n a t e due to the mesh distortion while adopting the conventional Lagrangian FEM. e coupled Lagrangian/Eulerian method cane ectivelyavoidthemeshdistortionbutwouldincrease the computational cost and reduce the computational eciency.
S P Hi so n eo ft h em e s h -f r e ep a r t i c l em e t h o d si n Lagrangian frame, which has been widely used in various elds since it is rst invented to solve the astrophysical problems. In the SPH method, a series of particles with some physicalquantity,suchasmassandvelocity,areusedtoexpressthe continuousmaterial. emeshdistortioncanbewellavoided due to no structural mesh among these particles. erefore, it is a powerful method for solving the multiphysics ow and large deformation problems. Meanwhile, the FEM is suitable for the simulation of mechanical characteristics of s o l i dm a t e r i a l s .N o w a d a y s ,t h ec o u p l e dS P H / F E Mm e t h o d has been veri ed to solve the uid-solid interaction problem commendably, and it could overcome the disadvantages associated with the Lagrangian/Eulerian method [ , ]. erefore, the coupled SPH/FEM method was adopted in this paper to simulate the rock fragmentation process by the water jet impact. e coupling program is shown in Figure : the le part shows the computational process of SPH while the right part indicates the FEM process. Meanwhile, the SPH and FEM parts are combined by the node to surface contact algorithm [ ], where the slave part was de ned with SPH particles and the master part was de ned with nite elements.
. . eory of SPH. Compared with the nite element m e t h o d ,ak e r n e la p p r o x i m a t i o ni su s e di nS P Hb a s e do n randomly distributed interpolation points with no assumptions about which points are neighbors to calculate spatial derivatives. Considering a problem domain Ω discretized by a group of particles by collection of particles and assuming kernel function W has a compact supporting domain with a radius of ℎ. ea p p r o x i m a t i o n (x) and its discretized di erential form ∇ ( ) at point can be respectively expressed as [ , ] ( where is the smooth kernel function with B-spline type, ℎ is the smooth length which de nes the supporting domain of the particle, x and x ὔ are the location vectors in di erent position, and is the total number of the particles, including the particle within the supporting domain of the given particle , represents those in uenced particles nearby the particle , is the mass of particle , is the density of particle .
e equations of conservation governing the evolution of mechanical variables can be expressed as follows: conservation of mass: ( ) conservation of energy: where is the pressure, is the coordinate of particle in direction, and are the stress and strain tensor of particle ,and and are the contravariant indexes, is the viscosity coe cient of uid, V represents the relative velocity between two particles in direction.
. . Material Models . . . Equation of State for Water Jet. As mentioned above, there is a large deformation of water in the rock fragmentation process. erefore, the shock equation of state is adopted to describe the mechanical characteristics of water = 1 + 2 2 + 3 3 + 0 + 1 ≥ 0 where is the water pressure, is the density of water jet, represents the compression ratio of water jet, ≥0represents that water is in compression, and, on the contrary of it, is in expansion; is the internal energy of water, and 1 , 2 , 3 , 1 , 2 , 0 ,and 1 are the material constants of water medium. In this numerical study, = kg/m 3 , 1 = . GPa, 2 = . GPa, 3 = . GPa, 1 = . GPa, 2 = , 0 = . ,and 1 = . [ ]. . . . Equation of State and Failure Criterion for Rock. e pressure or deformation of the rock is relative small under the dynamic load by water jet impact, and their variations have less in uence on the thermodynamic entropy. erefore, the pressure variation could be deemed to be only related to the density and volume variations of the rock element. Consequently,thelinearequationofstateisadoptedtodescribethe rock mechanical properties, which is very suitable for solving the dynamic problem with small deformation and pressure. e equation can be expressed as [ ] where is the rock pressure, is bulk elastic modulus of the rock, and 0 represent real density and reference density of the rock, respectively. In order to investigate the mechanism of rock fragmentation under the dynamic load by water jet impact, the maximum principal stress criterion is adopted to describe the rock failure behavior. Once the maximum tensile or shear principal stresses exceed the rock dynamic tensile or shear strength, the rock element fails. e maximum principal stress criterion can be expressed as jet impact is shown in Figure . e water jet is simpli ed as arectangle( × mm). esmoothparticlesizeis . mm, and there are smooth particles of water jet. Similarly, the rock is also simpli ed as a rectangle ( × mm). e quadrilateral element is adopted to mesh rock with the side length of . mm, and there are elements of rock. e rock bottom is set as the free boundary, and its side is set as the no-re ection boundary. e numerical model in this paper all adopts the above geometric model and boundary conditions unless otherwise stated.

Simulation Results
. . Rock Fragmentation Process. Figure shows the rock statuses as a function of time by water jet impact. In this simulation, the rock bottom and side are set as the free boundary and no-re ection boundary, respectively, and the impact velocity of water jet is m/s. Due to the water jet impact, the rock element containing point ( mm, mm) sustains an extremely high pressure of . GPa, and the pressure versus time at point is shown in Figure .I n the initial stage, the high-density shear stress eld forms nearby the impact point under the e ect of the extremely high pressure. Complex variation of the pressure nearby the impact point is induced by the intricate stress condition due to the combined e ect of water jet impact load, propagating disturbance of stress wave, compressive energy release, and so on. e formation of crushing zone nearby the impact point is conducted by the high-density shear stress eld as shown in Figure  failure behavior and mechanism. As shown in Figure (c),the spall crack appears nearby the rock bottom because the stress wave re ects at the bottom which is set as the free boundary, and several radial cracks nearly propagate to the rock side at . s. At the time of s, several spall cracks formed as a result of the stress wave propagation, re ection, and superposition in rock bottom is shown in Figure (d),a n d the radial and spall cracks initiate and propagate alternately within this zone. We can notice that the depth of the crushing zone is . mm at the initial stage of .
s. However, it increases by . mm at s. e phenomenon indicates that the formation process of crushing zone is extremely transitory. Although the increase in the depth of crushing zone is inconspicuous, the zone width increases obviously. e reason is that the walls of crushing zone will su er the erosion damage under the combined action of compression and shear due to the return water jet. As mentioned above, the spall crack initiation and propagation are caused by the stress wave re ection at the free boundary. However, the phenomena are related to the load magnitude of water jet impact and attenuation characteristics of stress wave, as well as rock volume and mechanical properties. For example, the spall cracks will not appear under a certain impact load while the rock is very big. As compared in Figures (e)  both and directions increase with time at the compressive status. But the stress in y direction increases faster than that in directionbecausetheimpactloadisperpendiculartothe rockbottomatthebeginning.At . s, the maximum shear principal stress 12 reaches up to MPa, which is larger than the rock dynamic shear strength ( = MPa). en the element containing point fails in shear, and the shear stress declines to in the next time step ( . s). Because the maximum principal stress criterion is based on the shear stress equal to , the element stresses in x and y directions are equal to each other once the element fails. ere are two highlights of the simulation results: (i) the element failures in the crushing zone are all due to that the maximum shear principal stress reaches the rock dynamic shear strength. Although element stresses in crushing zone are di erent from each other before failure, the failure behaviors are similar. erefore, it is unnecessary to present the failure behaviors of the other elements in this paper; (ii) the element containing point fails while the maximum shear principal stress reaches MPa but not MPa in theory, which is caused by the time step e ect in the simulation. Of course, the phenomenon can be avoided if the time step is set to be small enough. In order to improve the computational e ciency of the developed numerical model, the time step is automatically calculated according to the mesh size, and the similar phenomena will not be explained anymore below.
. . Crack Formation. In Figure , the element failure is used to simulate the crack initiation and propagation approximatively. e element containing point ( mm, mm) in the crack initiation zone is selected to investigate the crack initiation because a radial crack propagation crosses it. e stresses as a function of time at point are recorded shown in Figure . es t r e s si ny direction is always compressive stress before failure, so the crack through the point is a radial crack. e maximum principal stress 1 alternates between compression and tension status twice due to the coupling e ect between the stresses and element strain energy release; then it reaches . MPa which is larger than Shock and Vibration the rock dynamic tensile strength ( = 5 7MPa). In the next time step, the element containing point fails, resulting in the radial crack. A er the element fails, the shear stress decreases to , and the normal stresses , , 1 ,a n d 2 areequal.Althoughtheelementfailsbecausethemaximum principal stress 1 exceeds the rock dynamic tensile strength, its maximum shear principal stress 12 ( . MPa) is very close to the rock dynamic shear strength at this time. us i tc a nb es e e nt h a tt h ee l e m e n t sn e a r b yp o i n t m a yf a i li n tension or shear, in which the radial crack initiated nearby point can also be con rmed as shown in Figure . e rock bottom is set as the free boundary in the simulation, which results in the appearance of spall crack at rock bottom. e element containing point ( mm, mm) is selected to investigate the mechanism of spall crack and its stresses as a function of time are recorded as shown in Figure . es t r e s sw a v ep r o p a g a t e st op o i n t a t . s; then, the element stresses in both and direction increase with time at the compressive status. Because the distance betweenpoint andthefreeboundaryisshort,theelement s, and the element fails because the maximum principal stress exceeds the rock dynamic tensile strength. In particular, the stress in direction is larger than that in direction at tension status before failure, which indicates that the crack through point is a spall crack. Furthermore, the stress variations in and direction both change with time, which shows that the element stresses around the point are very complex. e crack through point may propagate along x direction or y direction due t otheco u p lin ge ectbetw eenthed yna mics tr es swa v ea nd re ected stress wave, which con rms that the initiation and propagation of the radial and spall crack are decussate just as shown in Figure . As shown in Figure , there are transverse cracks in the upper and lower parts of the rock, and the transverse crack in the lower part due to the re ected stress wave at the free boundary has been investigated above. To study the mechanism of transverse crack in the upper part, the stress change of the element containing point ( mm, mm) is recorded as shown in Figure . en o r m a ls t r e s s e sa r e basically coincident with the maximum principal stresses because the shear stress variation is very small, and the stress in direction and the maximum principal stress 2 are both negative before failure. us, the maximum principal stress 1 and the stress in direction reach . MPa, and . MPa respectively, and they conduce to the crack through point to open along y direction and propagate along direction. e cracks through point and are similar, but the failure o fp o i n t l a g sb e h i n dt h a to fp o i n t b ya b o u t . sa si t canbeseenfromFigure . Furthermore, the re ected stress wave action on point happens before that on point , which indicates that the crack through point is actually a radial crack but not a spall crack.

. . Stress Wave Propagation and Attenuation.
e element stresses in di erent positions (shown in Figure )  s when the pressure falls from the maximum value to the rock dynamic compressive strength, and it shows that action of the stress wave on the rock elements is a load/unload process. erefore, the rock fragmentation by water jet impact can also be regarded as a load/unload process, which can provide a theoretical basis for the formation of the crushing zone at the initial stage.

Influence Factors on Crack Initiation and Propagation
Many factors in uence the rock fragmentation process by water jet impact, such as boundary conditions, water jet velocity and diameter, rock con ning pressure, and structure plane. Based the above-mentioned numerical model, the in uence of these factors on rock fragmentation will be investigated in the following sections.
. . Boundary Condition. e rock bottom is set as the free boundary and no-re ection boundary, respectively; then the rock fragmentation with the two di erent boundary conditions is investigated based on the developed numerical model. Figure shows the rock failure status at the time of s. e free boundary means that the stress wave will be re ected when it reaches the rock bottom; however, it will go through with the no-re ection boundary which is usually usedtosimulatethewavepropagationinthein nitemedium. B a s e do nt h ea b o v ea n a l y s i s ,t h ef r e eb o u n d a r yc a nb r i n g about the spall cracks in the lower part of the rock. In contrast tothecasewithfreeboundary ,thereisnospallcracknearby the rock bottom because the stress wave goes through the interface. Due to the e ect of the re ected stress wave, the extent of rock fragmentation with the free boundary is better than that with the no-re ection boundary, which indicates that the free boundary can improve the rock fragmentation ability of water jet.
e velocity of water jet determinates the impact energy; therefore, it has a direct e ect on rock fragmentation. ree di erent velocities of m/s, m/s, and m/swereselectedinthesimulation,androckfailure behaviors at s were obtained as shown in Figures (a), (b),a n d (d), respectively. According to the simulation results, the extent of rock fragmentation increases with impactvelocity.Whilethevelocityis m/s,thestresswave strengthistooweaksothattheradialcrackisunabletopropagate to rock bottom or side. us, the rock fragmentation is mainly concentrated around the impact point. When the velocity is m/s, the radial crack length is longer than that of m/s, but there is few spall cracks at rock bottom. As for the velocity of m/s, the re ected stress wave from free boundary can bring out the formation of spall crack due to the large amount of impact energy, which in return increases the extent of rock fragmentation. As a consequence, it can be seen from the above analysis that the surface erosion of rock is primary at low impact velocity and the actual failure (such as radial and spall cracks) will occur only when impact velocity reaches up to a certain value.
. . Con ning Pressure. In the deep geotechnical engineering, the con ning pressure greatly in uences the rock fragmentation and the crack initiation and propagation. In order toavoidthee ectofthere ectedstresswaveonthenumerical simulation results, the rock bottom is set as the no-re ection boundary, and both two rock sides were applied with the conning pressure. Four con ning pressures of MPa, MPa, MPa, and MPa were separately applied to both two rock sides, and the rock fragmentation statuses at swere shown in Figure . e radial crack can propagate to the rock bottom when the con ning pressure is MPa or MPa, but it cannot propagate to the bottom with the con ning pressure of MPaor MPa. ereasonisthatthecon ningpressures can e ectively suppress the tensile stress perpendicular to the impact direction, and the inbibitional e ect increases with the con ning pressure. e crack length is . mm (as shown in Figure (c)) with the con ning pressure of MPa, but it decreases to . mm (as shown in Figure (d))when the con ning pressure is MPa. is phenomenon indicates that the con ning pressure can not only inhibit the tensile stress perpendicular to the impact direction, but also restrain the cracks propagation rate along the pressure direction. As a result, the stress concentration area is formed nearby the impactpointunderthecombinede ectofthestresswaveand the con ning pressure. us, the severe rock fragmentation and high-density tensile crack nearby the impact point are caused by the stress concentration.
. . Structure Plane. e relative position between the impact point and the structure plane also in uences the rock fragmentation and crack initiation and propagation. e sa n d s t o n ewi thawid tho f . mmwa su sedt osim ula t eth e structure plane, and the distances of mm, mm, and mm as well as the angles of ∘ and ∘ were set in the following simulations. e rock bottom and sides were a l ls e ta st h en o -r e e c t i o nb o u n d a r yt oa v o i dt h ee e c to f there ectedstresswa veonthesim ula tionresults. erock fragmentation statuses at s with di erent structure planes were obtained as shown in Figure . Similar with the free boundary, part of the stress wave was re ected back and the other part would go through the structure plane. us, the stress wave strength will be weakened greatly. e stress wave through the structure plane can hardly conduce to the rock fragmentation outside the structure plane. But many spall cracks were initiated inside the structure plane because the re ected compressive stress wave changed into the tensile stress wave. Moreover, the extent of rock fragmentation inside thestructureplaneincreaseswiththedecreaseinthedistance between the impact point and the structure plane. As shown in Figures (e) and (f), the spall cracks would propagate along the structure plane direction due to the stress wave re ection. According to the de nition of the speci c energy (the energy required to crush a unit volume of rock) [ ], the rock fragmentation e ect is better with a larger crushing zone under the same impact energy. e rock is overcrushed when the distance is only mm, resulting in a high speci c energy. e rock cannot be crushed e ectively while the distance is mm, whereas it can be mainly crushed to moderate lumpiness with a distance of mm. As a consequence, a n o p t i m a ld i s t a n c eb e t w e e nt h ei m p a c tp o i n ta n dt h er o c k structural plane can be obtained once the other conditions are determined.

Conclusions
In the present study, the coupled SPH/FEM method is implemented to simulate the rock fragmentation under the impact load of water jet. e in uence factors on the rock fragmentation and crack initiation and propagation are extensively investigated. From the numerical simulation results, we come to the following conclusions.
( ) e rock nearby the impact point is crushed severely due to the extremely high local pressure at the initial impact stage. en, the pressure attenuates sharply when it propagates in the rock medium due to the energy dissipation. erefore, the rock fragmentat i o nb yw a t e rj e ti m p a c tc a na l s ob er e g a r d e da sa load/unload process.
( ) e rock fragmentation modes at the upper and lower parts from the numerical simulation are consistent with that in the experiment. e rock fragmentation by water jet impact is due to the combined action of shear and tensile failure. e crushing zone nearby the impact point is mainly caused by the shear failure as a result of the high compressive stress, while the radial crack initiation and propagation is aroused by the tensile failure. However, spall crack nearby the rock bottom is caused by the re ected stress wave. At the same time, the micromechanism of the cracks initiation and propagation is investigated by analyzing the element stresses as a function of time in di erent positions ( Figures -).
( ) ee e c to ft h ef r e eb o u n d a r yo nt h er o c kf r a gmentation is very signi cant, and the spall crack can improve the fragmentation extent. e surface erosionofrockisprimaryatlowimpactvelocityand t h ea c t u a lf a i l u r e( s u c ha sr a d i a la n ds p a l lc r a c k s ) will occur only when impact velocity reaches up to a certain value. e con ning pressure can restrain the stress wave and cracks propagation, but it will conduce to a severe rock fragmentation as a result of the stress concentration closed to the impact point. e e ect of structure plane on rock fragmentation i ss i m i l a rt ot h a to ft h ef r e eb o u n d a r y .A no p t i m a l distance between the impact point and the rock structural plane can be obtained once the other conditions are determined. Meanwhile, the structure plane can lead the spall crack to propagate along the plane direction. References [ ] X.F .Y a n g,X.H.Li,a ndY .Y .L u," W ea rcha racteristicso fthe cemented carbide blades in drilling limestone with water jet, " International Journal of Refractory Metals and Hard Materials, vol. , no. , pp.