#### Abstract

To investigate the rock fragmentation and its influence factors under the impact load of water jet, a numerical method which coupled finite element method (FEM) with smoothed particle hydrodynamics (SPH) was adopted to simulate the rock fragmentation process by water jet. Linear and shock equations of state were applied to describe the dynamic characteristics of rock and water, respectively, while the maximum principal stress criterion was used for the rock failure detection. The dynamic stresses at the selected element containing points in rock are computed as a function of time under the impact load of water jet. The influences of the factors of boundary condition, impact velocity, confining pressure, and structure plane on rock dynamic fragmentation are discussed.

#### 1. Introduction

Rock fragmentation by the water jet technology has been widely used in mining, petroleum drilling, civil construction, gas drainage, and cleaning operations [1–4]. 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 [5]. The rock fragmentation by water jet impact not only is a difficult problem in mechanics, but also presents a new challenge for physics [6], 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 [7]. Meanwhile, they studied the fracture of brittle material and noted that shear and tear occurred on the rock surface due to the high speed liquid, and the so-called water-hammer pressure can cause substantial damage to brittle material [8]. 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 different stress waves [9]. Momber investigated the rock erosion due to the liquid jet impact and observed that the failure was mainly attributed to the lateral jetting [10]. 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 [11, 12]. 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 [13]. In recent years, many researchers adopted the numerical method to investigate the failure of rock or rock-like materials by water jet impact, and their research works were mainly focused on the influences of the impact velocity, diameter, standoff distance, and incidence angle of water jet, as well as the confining pressure and type of rock on the erosion depth, diameter, damage, and stress distributions of rock [14–18]. However, there is less investigation on the rock fragmentation mechanism by the compression, shear, or tensile stress. Although many significant results have been published, they are far from complete for the numerical study of rock fragmentation by water jet impact.

The present paper is aiming at revealing the rock fragmentation mechanism and explaining the reasons for crushing zone formation, crack initiation, and propagation under the impact load of water jet. To achieve this goal, 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 effects of boundary condition, impact velocity, confining pressure, and structure plane on the rock failure and crack initiation and propagation.

#### 2. Numerical Model

##### 2.1. Coupled SPH/FEM Method

As we all know, the large deformation of water exists in the rock fragmentation process by water jet impact. The calculation would easily terminate due to the mesh distortion while adopting the conventional Lagrangian FEM. The coupled Lagrangian/Eulerian method can effectively avoid the mesh distortion but would increase the computational cost and reduce the computational efficiency.

SPH is one of the mesh-free particle methods in Lagrangian frame, which has been widely used in various fields since it is first invented to solve the astrophysical problems. In the SPH method, a series of particles with some physical quantity, such as mass and velocity, are used to express the continuous material. The mesh distortion can be well avoided due to no structural mesh among these particles. Therefore, it is a powerful method for solving the multiphysics flow and large deformation problems. Meanwhile, the FEM is suitable for the simulation of mechanical characteristics of solid materials. Nowadays, the coupled SPH/FEM method has been verified to solve the fluid-solid interaction problem commendably, and it could overcome the disadvantages associated with the Lagrangian/Eulerian method [19, 20]. Therefore, the coupled SPH/FEM method was adopted in this paper to simulate the rock fragmentation process by the water jet impact. The coupling program is shown in Figure 1: the left 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 [21], where the slave part was defined with SPH particles and the master part was defined with finite elements.

##### 2.2. Theory of SPH

Compared with the finite element method, a kernel approximation is used in SPH based on 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 has a compact supporting domain with a radius of . The approximation and its discretized differential form at point can be respectively expressed as [19, 20] where is the smooth kernel function with B-spline type, is the smooth length which defines the supporting domain of the particle, and are the location vectors in different position, and is the total number of the particles, including the particle within the supporting domain of the given particle , represents those influenced particles nearby the particle , is the mass of particle , is the density of particle .

The equations of conservation governing the evolution of mechanical variables can be expressed as follows:

conservation of mass:

conservation of momentum:

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 coefficient of fluid, represents the relative velocity between two particles in direction.

##### 2.3. Material Models

###### 2.3.1. Equation of State for Water Jet

As mentioned above, there is a large deformation of water in the rock fragmentation process. Therefore, the shock equation of state is adopted to describe the mechanical characteristics of water
where is the water pressure, is the density of water jet, represents the compression ratio of water jet, represents that water is in compression, and, on the contrary of it, is in expansion; is the internal energy of water, and , , , , , , and are the material constants of water medium. In this numerical study, = 1000 kg/m^{3}, = 2.2 GPa, = 9.54 GPa, = 14.6 GPa, = 2.2 GPa, = 0, = 0.28, and = 0.28 [22].

###### 2.3.2. Equation of State and Failure Criterion for Rock

The pressure or deformation of the rock is relative small under the dynamic load by water jet impact, and their variations have less influence on the thermodynamic entropy. Therefore, the pressure variation could be deemed to be only related to the density and volume variations of the rock element. Consequently, the linear equation of state is adopted to describe the rock mechanical properties, which is very suitable for solving the dynamic problem with small deformation and pressure. The equation can be expressed as [21] where is the rock pressure, is bulk elastic modulus of the rock, and 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. The maximum principal stress criterion can be expressed as
where and are the maximum tensile principal stress of rock element, represents the rock dynamic tensile strength, is the maximum shear principal stress, and is rock dynamic shear strength. However, it is known that dynamic fracture simulation under high strain rate requires dynamic strength, and this value should increase with increasing strain rate. But the relationship between dynamic strength and strain rate for the target rock is not available. Therefore, we have to select a dynamic strength in this simulation and the impact target with the following mechanical properties: reference density = 2200 kg/m^{3}, elasticity modulus = 58.9 GPa, dynamic tensile strength = 57 MPa, dynamic shear strength = 192 MPa, and poisson ratio = 0.22.

##### 2.4. Geometric Model and Boundary Conditions

The geometric model of rock fragmentation under dynamic load by water jet impact is shown in Figure 2. The water jet is simplified as a rectangle (2 × 20 mm). The smooth particle size is 0.25 mm, and there are 640 smooth particles of water jet. Similarly, the rock is also simplified as a rectangle (50 × 30 mm). The quadrilateral element is adopted to mesh rock with the side length of 0.25 mm, and there are 24000 elements of rock. The rock bottom is set as the free boundary, and its side is set as the no-reflection boundary. The numerical model in this paper all adopts the above geometric model and boundary conditions unless otherwise stated.

#### 3. Simulation Results

##### 3.1. Rock Fragmentation Process

Figure 3 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-reflection boundary, respectively, and the impact velocity of water jet is 800 m/s. Due to the water jet impact, the rock element containing point 1 (0 mm, 29 mm) sustains an extremely high pressure of 1.49 GPa, and the pressure versus time at point 1 is shown in Figure 4. In the initial stage, the high-density shear stress field forms nearby the impact point under the effect 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 effect of water jet impact load, propagating disturbance of stress wave, compressive energy release, and so on. The formation of crushing zone nearby the impact point is conducted by the high-density shear stress field as shown in Figure 3(a). Meanwhile, radial crack initiation around the crushed zone forms the crack initiation zone. Then, as shown in Figure 3(b), the crack propagation zone is formed at the time of 3.1 *μ*s. It is worth noting that a centerline crack is formed perpendicular to the bottom in rock interior, and the radial crack is generally symmetrical about the impact direction. These phenomena appear as a result of the rock continuity and isotropy which have no effect on the rock failure behavior and mechanism. As shown in Figure 3(c), the spall crack appears nearby the rock bottom because the stress wave reflects at the bottom which is set as the free boundary, and several radial cracks nearly propagate to the rock side at 6.2 *μ*s. At the time of 20 *μ*s, several spall cracks formed as a result of the stress wave propagation, reflection, and superposition in rock bottom is shown in Figure 3(d), and the radial and spall cracks initiate and propagate alternately within this zone. We can notice that the depth of the crushing zone is 4.2 mm at the initial stage of 0.12 *μ*s. However, it increases by 0.2 mm at 30 *μ*s. The 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. The reason is that the walls of crushing zone will suffer 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 reflection 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 3(e) and 3(f), the rock failure modes of simulation and experiment (conducted by Lu et al.) at the upper and lower part are basically consistent [18]. Hence, the developed numerical model in this paper can well reappear in the rock fragmentation process by the water jet impact.

**(a) 1.2 μs**

**(b) 3.1 μs**

**(c) 6.2 μs**

**(d) 20 μs**

**(e) 50 μs**

**(f) Experimental result**

##### 3.2. Formation Mechanism of the Crushed Zone

To investigate the formation mechanism of the crushing zone by water jet impact, the stresses of the element containing point 2 (0 mm, 27 mm) in the crushing zone as a function of time are recorded as shown in Figure 5. Because of a certain distance between point 2 and the impact point, the stress wave propagates to point 2 at 0.15 *μ*s; then, the element stresses in both *x* and *y* directions increase with time at the compressive status. But the stress in *y* direction increases faster than that in *x* direction because the impact load is perpendicular to the rock bottom at the beginning. At 0.36 *μ*s, the maximum shear principal stress reaches up to 196 MPa, which is larger than the rock dynamic shear strength = 192 MPa). Then the element containing point 2 fails in shear, and the shear stress declines to 0 in the next time step (0.365 *μ*s). Because the maximum principal stress criterion is based on the shear stress equal to 0, the element stresses in *x* and *y* directions are equal to each other once the element fails. There 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 different from each other before failure, the failure behaviors are similar. Therefore, it is unnecessary to present the failure behaviors of the other elements in this paper; (ii) the element containing point 2 fails while the maximum shear principal stress reaches 196 MPa but not 192 MPa in theory, which is caused by the time step effect 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 efficiency 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.

##### 3.3. Crack Formation

In Figure 3, the element failure is used to simulate the crack initiation and propagation approximatively. The element containing point 3 (0 mm, 25 mm) in the crack initiation zone is selected to investigate the crack initiation because a radial crack propagation crosses it. The stresses as a function of time at point 3 are recorded shown in Figure 6. The stress in *y* direction is always compressive stress before failure, so the crack through the point 3 is a radial crack. The maximum principal stress alternates between compression and tension status twice due to the coupling effect between the stresses and element strain energy release; then it reaches 60.3 MPa which is larger than the rock dynamic tensile strength MPa). In the next time step, the element containing point 3 fails, resulting in the radial crack. After the element fails, the shear stress decreases to 0, and the normal stresses , , , and are equal. Although the element fails because the maximum principal stress exceeds the rock dynamic tensile strength, its maximum shear principal stress (177.9 MPa) is very close to the rock dynamic shear strength at this time. Thus it can be seen that the elements nearby point 3 may fail in tension or shear, in which the radial crack initiated nearby point 3 can also be confirmed as shown in Figure 3.

The rock bottom is set as the free boundary in the simulation, which results in the appearance of spall crack at rock bottom. The element containing point 4 (6 mm, 7 mm) is selected to investigate the mechanism of spall crack and its stresses as a function of time are recorded as shown in Figure 7. The stress wave propagates to point 4 at 3.5 *μ*s; then, the element stresses in both *x* and *y* direction increase with time at the compressive status. Because the distance between point 4 and the free boundary is short, the element stresses are affected by the reflected wave at the first time after 0.75 *μ*s. Furthermore, the maximum principal stress reaches 57.7 MPa at 6.42 *μ*s, and the element fails because the maximum principal stress exceeds the rock dynamic tensile strength. In particular, the stress in *y* direction is larger than that in *x *direction at tension status before failure, which indicates that the crack through point 4 is a spall crack. Furthermore, the stress variations in *x* and *y* direction both change with time, which shows that the element stresses around the point 4 are very complex. The crack through point 4 may propagate along *x* direction or *y* direction due to the coupling effect between the dynamic stress wave and reflected stress wave, which confirms that the initiation and propagation of the radial and spall crack are decussate just as shown in Figure 3.

As shown in Figure 3, there are transverse cracks in the upper and lower parts of the rock, and the transverse crack in the lower part due to the reflected 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 5 (10 mm, 28 mm) is recorded as shown in Figure 8. The normal stresses are basically coincident with the maximum principal stresses because the shear stress variation is very small, and the stress in *x* direction and the maximum principal stress are both negative before failure. Thus, the maximum principal stress and the stress in *y *direction reach 58.25 MPa, and 58.23 MPa respectively, and they conduce to the crack through point 5 to open along *y* direction and propagate along *x* direction. The cracks through point 4 and 5 are similar, but the failure of point 4 lags behind that of point 5 by about 3.4 *μ*s as it can be seen from Figure 9. Furthermore, the reflected stress wave action on point 4 happens before that on point 5, which indicates that the crack through point 5 is actually a radial crack but not a spall crack.

##### 3.4. Stress Wave Propagation and Attenuation

The element stresses in different positions (shown in Figure 3) as a function of time have been recorded as shown in Figure 9. After the water jet starts to impact the rock, the pressure of the element containing point 1 achieves the maximum value of 1.49 GPa at 0.285 *μ*s. After 1.35 *μ*s, the maximum pressure value of point 3 falls down to 202 MPa due to energy dissipation associated with the stress wave propagation. And the maximum pressure is only 65.2 MPa when the stress wave propagates to point 6, which indicates that the stress wave intensity weakens sharply during the propagation in the rock medium. The stress wave is not able to crush rock (shear failure) when the peak pressure is smaller than the rock dynamic compressive strength. At this time, they can only conduce to the radial and spall cracks. The time interval is only about 1.35 *μ*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. Therefore, 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.

#### 4. Influence Factors on Crack Initiation and Propagation

Many factors influence the rock fragmentation process by water jet impact, such as boundary conditions, water jet velocity and diameter, rock confining pressure, and structure plane. Based the above-mentioned numerical model, the influence of these factors on rock fragmentation will be investigated in the following sections.

##### 4.1. Boundary Condition

The rock bottom is set as the free boundary and no-reflection boundary, respectively; then the rock fragmentation with the two different boundary conditions is investigated based on the developed numerical model. Figure 10 shows the rock failure status at the time of 20 *μ*s. The free boundary means that the stress wave will be reflected when it reaches the rock bottom; however, it will go through with the no-reflection boundary which is usually used to simulate the wave propagation in the infinite medium. Based on the above analysis, the free boundary can bring about the spall cracks in the lower part of the rock. In contrast to the case with free boundary, there is no spall crack nearby the rock bottom because the stress wave goes through the interface. Due to the effect of the reflected stress wave, the extent of rock fragmentation with the free boundary is better than that with the no-reflection boundary, which indicates that the free boundary can improve the rock fragmentation ability of water jet.

**(a) Free boundary**

**(b) No-reflection boundary**

##### 4.2. Impact Velocity

The velocity of water jet determinates the impact energy; therefore, it has a direct effect on rock fragmentation. Three different velocities of 300 m/s, 500 m/s, and 800 m/s were selected in the simulation, and rock failure behaviors at 20 *μ*s were obtained as shown in Figures 11(a), 11(b), and 3(d), respectively. According to the simulation results, the extent of rock fragmentation increases with impact velocity. While the velocity is 300 m/s, the stress wave strength is too weak so that the radial crack is unable to propagate to rock bottom or side. Thus, the rock fragmentation is mainly concentrated around the impact point. When the velocity is 500 m/s, the radial crack length is longer than that of 300 m/s, but there is few spall cracks at rock bottom. As for the velocity of 800 m/s, the reflected 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.

**(a) 300 m/s**

**(b) 500 m/s**

##### 4.3. Confining Pressure

In the deep geotechnical engineering, the confining pressure greatly influences the rock fragmentation and the crack initiation and propagation. In order to avoid the effect of the reflected stress wave on the numerical simulation results, the rock bottom is set as the no-reflection boundary, and both two rock sides were applied with the confining pressure. Four confining pressures of 10 MPa, 20 MPa, 40 MPa, and 60 MPa were separately applied to both two rock sides, and the rock fragmentation statuses at 20 *μ*s were shown in Figure 12. The radial crack can propagate to the rock bottom when the confining pressure is 10 MPa or 20 MPa, but it cannot propagate to the bottom with the confining pressure of 40 MPa or 60 MPa. The reason is that the confining pressures can effectively suppress the tensile stress perpendicular to the impact direction, and the inbibitional effect increases with the confining pressure. The crack length is 17.3 mm (as shown in Figure 12(c)) with the confining pressure of 40 MPa, but it decreases to 15.1 mm (as shown in Figure 12(d)) when the confining pressure is 60 MPa. This phenomenon indicates that the confining 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 impact point under the combined effect of the stress wave and the confining pressure. Thus, the severe rock fragmentation and high-density tensile crack nearby the impact point are caused by the stress concentration.

**(a) 10 MPa**

**(b) 20 MPa**

**(c) 40 MPa**

**(d) 60 MPa**

##### 4.4. Structure Plane

The relative position between the impact point and the structure plane also influences the rock fragmentation and crack initiation and propagation. The sandstone with a width of 0.5 mm was used to simulate the structure plane, and the distances of 10 mm, 15 mm, and 20 mm as well as the angles of 30° and 60° were set in the following simulations. The rock bottom and sides were all set as the no-reflection boundary to avoid the effect of the reflected stress wave on the simulation results. The rock fragmentation statuses at 20 *μ*s with different structure planes were obtained as shown in Figure 13. Similar with the free boundary, part of the stress wave was reflected back and the other part would go through the structure plane. Thus, the stress wave strength will be weakened greatly. The 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 reflected compressive stress wave changed into the tensile stress wave. Moreover, the extent of rock fragmentation inside the structure plane increases with the decrease in the distance between the impact point and the structure plane. As shown in Figures 13(e) and 13(f), the spall cracks would propagate along the structure plane direction due to the stress wave reflection. According to the definition of the specific energy (the energy required to crush a unit volume of rock) [23], the rock fragmentation effect is better with a larger crushing zone under the same impact energy. The rock is overcrushed when the distance is only 10 mm, resulting in a high specific energy. The rock cannot be crushed effectively while the distance is 20 mm, whereas it can be mainly crushed to moderate lumpiness with a distance of 15 mm. As a consequence, anoptimal distance between the impact point and the rock structural plane can be obtained once the other conditions are determined.

**(a) Structural plane**

**(b) 10 mm**

**(c) 15 mm**

**(d) 20 mm**

**(e) 30°**

**(f) 60°**

#### 5. Conclusions

In the present study, the coupled SPH/FEM method is implemented to simulate the rock fragmentation under the impact load of water jet. The influence factors on the rock fragmentation and crack initiation and propagation are extensively investigated. From the numerical simulation results, we come to the following conclusions.(1)The rock nearby the impact point is crushed severely due to the extremely high local pressure at the initial impact stage. Then, the pressure attenuates sharply when it propagates in the rock medium due to the energy dissipation. Therefore, the rock fragmentation by water jet impact can also be regarded as a load/unload process.(2)The rock fragmentation modes at the upper and lower parts from the numerical simulation are consistent with that in the experiment. The rock fragmentation by water jet impact is due to the combined action of shear and tensile failure. The 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 reflected 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 different positions (Figures 5–8).(3)The effect of the free boundary on the rock fragmentation is very significant, and the spall crack can improve the fragmentation extent. 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. The confining 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. The effect of structure plane on rock fragmentation is similar to that of the free boundary. An optimal 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.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

The authors would like to acknowledge the Foundation of National 863 Plan of China (2012AA062104), the National Natural Science Foundation of China (51375478), the project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (SZBF2011-6-B35), and the Graduate Education Innovation Project of Jiangsu Province (CXLX12_0948).