Impact characteristics and stagnation formation on a solid surface by a supersonic abrasive waterjet

A computational fluid dynamics (CFD) study of the impact characteristics and stagnation formation on a solid target surface by an abrasive waterjet at supersonic velocities is presented to understand the impact process. A CFD model is developed and verified by experimental water and particle velocities and then used to simulate the jet impact process. The trends of the stagnation formation and its effect on the jet flow with respect to the jetting and impacting parameters are amply discussed. It is found that stagnation formation at the impact site increases with an increase in the impact time, nozzle standoff distance and nozzle diameter, while the initial peak velocity at the nozzle exit has little effect on the size of the stagnation zone. It is shown that stagnation markedly changes the water and particle flow direction, so that the particle impact angle is varied and the jet impact area is enlarged. The jet structure may be classified to have a free jet flow region, a jet deflection region with a stagnation zone and a wall jet region. Furthermore, the stagnation affects significantly the waterjet and particle energy transferred to the target surface. The average particle velocity across the jet is reduced by approximately one third due to the damping effect of the stagnation under the conditions considered in this study.

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. V avg,w , V avg,p average velocities of water and particle V w , V p velocities of water and particle z axial distance from the nozzle exit α p particle impact angle ρ fluid density ρ w , ρ a density of water and air ρ p particle density μ dynamic viscosity μ t turbulent viscosity ω specific turbulent dissipation θ jet impact angles

Introduction
Abrasive waterjet (AWJ) finds wide applications in the manufacturing industry, particularly for the processing of difficult-to-machine materials, such as composites and ceramics, due to its various advantages [1]. In a high pressure AWJ, the jet is formed by a fine nozzle attached to a mixing chamber, before entering into an atmosphere as a free jet to impinge a target workpiece. The erosion process in AWJ machining is mainly affected by particle impact velocity and impact angle in addition to the target material properties. A deeper understanding of the waterjet and the particle dynamic characteristics in the AWJ is therefore essential to further improve the AWJ machining performance. Numerical analysis is an effective method to study the AWJ dynamic characteristics since it can directly visualise and record the local results of the AWJ, such as water and particle velocities, which are in general very difficult to measure through experiments. The numerical studies reported on the AWJ processing of materials can be divided into two main streams, i.e. the impact mechanics of particles on a target material [2][3][4][5][6][7] and the AWJ flow characteristics [8][9][10][11]. While the effort to understand the AWJ flow characteristics is small, the knowledge of the impact phenomenon for waterjet and AWJ on a target is even in more dearth [12][13][14][15][16][17][18][19].
Computational fluid dynamic (CFD) has been proven to be a powerful numerical method to study the AWJ flow characteristics [9,18,19]. Liu et al [9] developed a CFD model to study water and particle velocities and trajectories in an AWJ after exiting a nozzle. The CFD model provided an important insight into the jet flow characteristics, such as the jet velocity variation across and along the jet flow. Data from the CFD model [9] was further utilised to develop a mathematical model for predicting the particle velocity inside AWJs [20].
Since the material removal in AWJ cutting is caused by the impact of the water stream and particles, it is important to understand the jet impact characteristics. For this purpose, a number of important investigations using numerical approaches, such as CFD [14,[17][18][19], have been reported. The study of the AWJ impact phenomenon for a polishing process was conducted by Lv et al [14]. The effect of the jet impact angle was observed under the condition of low water pressure and high standoff distance. It was found that an elliptical contour of AWJ impact pressure occurred at non-perpendicular impact conditions, with the peak pressure located at the vicinity of the centre area. Ahmed et al [17] developed a model to investigate the characteristics of particle impact on a target workpiece by considering the effect of kerf wall radius corresponding to the depth of cut. It was ascertained that particle impact angles depended on the kerf wall radius, and an increase in the radius of curvature resulted in a decrease in the particle impact angles. A numerical-empirical technique was used by Kowsari et al [18] and Nouraei et al [19] to simulate the impact and material removal process for both brittle (sintered ceramics) and ductile materials by AWJs. It was found that jet impact angle affected the erosion rate. Recently, Schwartzentruber et al [21,22] employed a fluidstructure interaction (FSI) technique coupling the finite element method (FEM) and CFD, to investigate the effects of AWJ cutting process on anisotropic materials. They found that the hydraulic AWJ shock impact at the initial cutting stage caused ply debonding and delamination [21].
From the previous investigations, a damping flow or stagnation zone has been observed at the impact site on the target surface [3,15,19]. Jafar et al [23] found that stagnation suddenly decreased particle velocities in their study of a micro-AWJ impingement on a brittle material. It was claimed [24] that stagnant formation may affect the machined kerf formation because the direction of jet and particle flow may be affected by the stagnant layer. Pang [25] investigated the micro-AWJ drilling process for amorphous glass and found that the bottom surface of the non-through drilled holes exhibited a 'W' shape because of the stagnation effect. A bulge was also created at the bottom of machined channels by AWJ milling in the case of high nozzle traverse speeds [26]. It has been further claimed that an oblique jet impact causes the stagnation zone to shift away from the jet centre [23], while an increase in the depth of cut in AWJ cutting generates a stronger stagnation [18]. Matsumura et al [13] studied a lowpressure AWJ polishing process on a glass workpiece, with and without a tapered mask. It was found that the tapered mask has a significant influence on the stagnation formation, whereby the mask provided a larger stagnation area than without the mask.
It is apparent from the previous investigations that stagnation does exists and affects the AWJ impact process and material removal. For the low-pressure or micro-AWJ cutting condition, it has been shown that stagnation has a significant effect on the jet impact energy and the highest reduction of the impact energy was found at the centre of the jet impact area [25]. However, there is little understanding of stagnation formation and its effect on the jet impact process, particularly under ultrahigh pressure, ultrahigh pressure conditions relevant to AWJ machining. This work is aimed to gain an insight into the AWJ impact characteristics. For this purpose, a CFD analysis is employed to develop an AWJ impact model, which after experimental verifications, is used to study the formation process of the stagnation zone, the effects of jetting and impact parameters on the stagnation formation, and the effect of stagnation on the jet and particle impact characteristics.

CFD model formulation
In AWJ processing, a high pressure AWJ is formed by a fine nozzle attached to a mixing chamber before entering into an atmosphere as a free jet to impinge a target surface. The abrasive particles are mixed with the waterjet in the mixing chamber by an assistance of air entrainment during the abrasive feeding process. The air entrainment creates a suction effect that the particles are accelerated and axially oriented before exiting the nozzle [27,28]. As it is a normal expectation in AWJ machining that the velocity of abrasives at the nozzle exit is close to that of their surrounding water through appropriate cutting head design [29]. Thus, in this study, the CFD model for the AWJ starts where the jet exits a nozzle and enters into a computational domain. The initial particle velocity at the nozzle exit or inlet boundary of the AWJ model is set to be equal to the velocity of the surrounding water. It should be noted that in the entrainment AWJ system commonly used for ultrahigh pressure AWJ, the jet inside the nozzle is not anything like flow in pipe. Thus, modelling the flow inside the nozzle is rather difficult and the reported models in the literature mostly incorrectly represent the actual flow.
This CFD study is conducted using ANSYS Fluent 16.2. The volume of fluid (VOF) model is applied for two-phase flow of air and water, whereas the discrete phase method (DPM) model is employed for the abrasive particles. In the VOF model, the water phase is treated as a primary phase. Surface tension at the air-water interface is set at the constant value of 0.072 N m −1 [9]. In the DPM model, garnet particles of 80 mesh (average diameter of 0.18 mm and density of 4.1 g cm −3 ) are used. The AWJ model is simulated under a transient state and turbulent flow, and the fluid is considered incompressible. In addition, a two-way coupled approach is adopted to accommodate the interaction between water and abrasives, in which the water phase is treated as a continuum and the abrasives as a discrete phase.

Geometry and boundary conditions
In this study, a three-dimensional (3D) computational domain of a cylindrical shape is chosen, as shown in figure 1(a), where D and H are the domain diameter and height, respectively. For the boundary conditions, the inner diameter of nozzle (Nozzle ID or d n ) is used as the inlet boundary. To perform the calculation, a suitable velocity profile at the inlet boundary has to be applied. According to Liu [30], a oneseventh power law distribution has deemed to be appropriate and is used in this study. It has been reported [31] that before the jet exits the nozzle, there is a small proportion of air mixture inside the jet, as the air entrined into the mixing chamber may flow outside the jet at the junction between the core jet surface and inner nozzle wall. Therefore, the air contamination inside the jet at the inlet boundary (i.e. nozzle exit) was neglected by setting the volume fraction of air as zero. The atmospheric pressure is applied at the pressure outlet regions, i.e. the top and side surfaces. In evaluating the top surface, the nozzle wall thickness is considered. The volume fraction of the air phase is set as unity at the outlet boundaries since only air enters into the computational domain. At the bottom boundary, a no-slip wall condition is selected to represent a target surface. The domain height of the model is selected with reference to the standoff distances commonly used in AWJ machining [29], that is, from 3 to 10 mm, which will be used in the simulation study later in this paper. The one-seventh power law inlet velocity profile is coded in user-defined function, and is given by where V (0,r) is the velocity profile at the nozzle exit and the other symbols are as defined in the Nomenclature. The jet diameter (d j ) is assumed to be equal to the nozzle inner diameter.
A hexahedral mesh with non-uniform spacing is applied for the mesh generation of the AWJ model, as shown in figure 1(b) where a more concentrated mesh is adopted at close to the wall or target surface in order to better understand the impact characteristics in that region.  velocity profiles in the jet radial direction at 1 mm downstream distance from the nozzle exit under different mesh sizes. It can be noticed that the difference between the water velocities when the element size changes from 0.025 to 0.050 mm is less than 5%. Thus, the element size of 0.050 mm is selected for shorter computation time. From about 0.5 mm to the impact wall, the mesh size was reduced by 5% after each grid toward the wall in z-direction with reference to other studies [32] and was found to give good convergence of the calculations. From figure 2(b), the velocity profiles have no discernible difference for the three domain diameters considered with the mesh size of 0.050 mm. Nevertheless, the domain diameter of 12 mm, rather than 8 mm, is selected to enable to observe such impact characteristics as particle rebounding in upward and side directions. Furthermore, as the mesh size is smaller than the particle size, a node-basedaveraging function is employed to investigate the interaction between the waterjet and particle phase by evaluating the effect of the particles among the neighbouring nodes of each mesh [33].

Governing equations
The governing equations for modelling the two-phase flow of water and air are based on the conservation of mass, momentum and energy. The equations can be expressed as follows [34]: The continuity equation: The momentum equation: The energy conservation equation: where α is the volume fraction in mesh, F is the external force,  vv is the Reynolds stress, E is energy, T is temperature, k eff is the effective thermal conductivity, and S h is the volumetric heat sources. As the AWJ model assumed no heat generation, the variable of S h is zero. The effect of abrasive particles on the fluid flow (e.g. drag force) is taken into account through the force variable (F) in equation (3).
For the two-phase flow of water and air, the mixture fluid density, ρ, is determined based on the volume of fraction, i.e.
where α w and α a are, respectively, the volume fraction of water and air in a mesh. For the turbulent flow, the shear-stress transport (SST) k-ω model is applied for the AWJ impact model to allow an additional wall function to bridge the viscosity-affected region between the wall and the fully turbulent region of the fluid flow near the wall [35]. It has been shown [18,19] that the k-ω turbulent model is efficient in capturing the turbulence of flow in the region close to the wall and results in better calculation accuracy. The SST k-ω model is given by where σ k and σ ω are the turbulent transport constants, G k and G ω represent, respectively, the generation of turbulent kinetic energy and the specific turbulent dissipation, Y k and Y ω represent the dissipation of k and ω, respectively, and S k and S ω are the user-defined source terms. The turbulent viscosity, μ t , can be determined based on Menter's model [35], which accounts for the transport of the turbulence shear stress in the definition of the turbulent viscosity. The viscosity model and its factors are given by [35]: where τ is the magnitude of strain rate, α * is the damp coefficient of the turbulent viscosity, α 1 is a constant, and y is the distance to the next surface. For high-Reynolds number, α * is approximately one and α 1 is 0.31 [35].
For the discrete particle phase, particle velocities and trajectories are considered by the influence of the surrounding fluid flow, including the gravity force (m p g), drag force (  F D ), virtual mass force (  F V ) and pressure gradient force (  F P ). Particle trajectories and velocities are calculated by integrating the force balance on each particle based on Newton's second law, so that the governing equation for the discrete phase becomes [34] ( ) The drag force per unit particle mass is given by [36] ( ) From equation (12), the drag coefficient (C D ) is dependent on the particle sphericity, which can be divided into spherical particles and non-spherical particles, respectively [37], and are given as: The drag force coefficient for a spherical particle: The drag force coefficient for a non-spherical particle: where a 1 and a 2 are the constants for a spherical particle over a range of particle Reynolds number, R esp , and the relevant data can be found from Morsi and Alexander [38]. The constants b 1 , b 2 , b 3 and b 4 can be calculated from: From equations (15)-(18), the particle sphericity (S p ) is defined by the ratio of the surface area of a sphere having the same volume as a particle divided by the actual surface area of the particle. Thus, the sphericity of one represents a round particle, while a particle having the sphericity less than one is a non-spherical particle.
The effects of the virtual mass and pressure gradient are taken into account in this model since the density ratio of water to particles is above 0.1 [39].
The virtual mass force is calculated from: V v m p p p and the pressure gradient force is given by: where C vm is the virtual mass factor defined as a constant with the value of 0.5 for a spherical particle [39].
The effect of particle-particle collisions is neglected in this study due to the relatively small number of particles inside the AWJ [30]. The abrasive particles are assumed as a rigid body with a uniform size and shape distribution.
In the CFD model, the exchange in heat and mass is negligible because of a very small heat generation and no particle phase transformation in the AWJ system [30]. Thus, only the momentum exchange is taken into account for the waterjet-particle interaction, and the exchanged force can be computed from: where F m is the exchanged force and Δt is the time step size. At the target surface, the impacting particles transfer energy to the target during the particle collision process, in which the particle impact velocity and impact angle as well as the material properties of the particle and target play a key role. The transferred energy may be consumed by a frictional loss, residual energy storage in the target after the impact and kinetic energy of the removed material [4]. In Fluent, the energy loss of the particles after colliding a target surface can be determined by using a coefficient of restitution, e. A coefficient of restitution less than unity implies energy loss during an impact, while equal to unity indicates no energy loss.
To investigate the particle impact characteristics in this study, cases with and without considering energy loss are conducted. However, the effect of erosion on the target surface has been ignored as the target is assumed to be a rigid wall boundary. This is considered valid for pure waterjet impact or the initial stage of AWJ impact on a solid surface. Further, the surface roughness in a micro/nano-meter scale is not considered to cause any discernible effect at such high pressures and is therefore neglected. In the case of ignoring energy loss, the coefficient of restitution is assumed to be unity, yielding the rebounding particles to retain the same velocity magnitude after colliding a target. By contrast, with the consideration of energy loss, the coefficients of restitution proposed by Grant and Tabakoff [40], which have been popularly accepted in simulating particle-wall interaction [41][42][43] are applied, which are given by [40]: where e t and e n are the normal and tangential coefficients of restitutions, respectively. From equations (22) and (23), the particle impact angle (α p ) is defined as the angle between the particle moving direction immediately before impacting a target surface and the surface of the target. In this study, the particle moving direction is determined based on the particle velocity components in the jet axial (v pz ) and radial (v pr ) directions, so that the particle impact angle is obtained by the relation of ( The finite volume method is used to convert the above governing equations into algebraic equations that can be solved numerically in the computational domain [44]. In this study, a coupled scheme for pressure-velocity is selected [45], while a second order upwind differencing scheme and an implicit scheme are employed for the space and time discretisation, respectively [46]. For numerical calculations, the time step of the fluid phase and the discrete particle phase may be set as small as possible to overcome any divergence problem. The fluid time step size can be estimated from the mesh size and the inlet velocity, while the time step for the particle phase is estimated corresponding to the fluid time step. Approximately 20 steps of particle phase calculation are recommended for each step of fluid phase calculation to ensure that the exchange between particles and water achieves the momentum conservation [47]. Thus, three levels of fluid time step sizes, 1×10 −6 , 1×10 −7 and 1×10 −8 s, and three levels of particle time step sizes, 5×10 −8 , 5×10 −9 and 5×10 −10 s, are considered to assess the effect of the fluid and particle time step sizes. It is found that the water velocity is independent of the fluid and particle time steps, and the variation of particle velocity is negligible when the fluid time step varies from 1×10 −7 to 1×10 −8 s. Nevertheless, in order to gain enough and accurate results of particle velocity and trajectory, especially where close to the target surface, the time step sizes for fluid and particle phases are set at 1×10 −8 s and 5×10 −10 s, respectively. The parallel version of the Fluent solver is used over 4-16 computer nodes on a multi-processor (Trentino cluster, UNSW) to reduce the calculation time of the CFD model.

Model verification
The assessment of the CFD model is made by comparing the model calculations with experimentally measured water and particle flow velocities, as the dynamic characteristics during an impact are difficult to measure, if not impossible. The water velocity at selected downstream locations was obtained by measuring the impact force using a dynamometer, while the technique of particle image velocimetry (PIV) and laser induced fluorescence (LIF) was employed for measuring particle velocities. It should be noted that the velocity of the jet flow at close to the target surface is expected to vary significantly. However, due to small jet size (<1 mm in diameter) and short travel distance at this region (about 0.5 mm), the measurement of velocity variation is almost impossible particularly considering the water rebounding and spreading from the target surface. Thus, the understanding of velocity variation in this region may be achieved by numerical studies.

Waterjet velocity
The average waterjet velocities were obtained from the jet impact force on the target surface measured by a dynamometer similar to [30,48]. The experimental setup is shown in figure 3. The system included a thin plat specimen, a jetting system, a Kistler type 9257A 3-component dynamometer, signal amplifiers, analogue-to-digital converters and LabVIEW data acquisition software. The dynamometer was covered for waterproof and the specimen was mounted on the top of it. In this experiment, four levels of water pressure (P), 200, 250, 300 and 350 MPa, were considered. For each pressure, the impact force was measured at five downstream locations (S d ) of 5, 15, 25, 35 and 45 mm. The nozzle used for all the tests has 0.76 mm in diameter and 76.2 mm in length. Each test was repeated for three repetitions and the average force was taken as the final reading.
The experimentally measured force were then transferred into jet impact velocity, V avg,w(exp.) , according to the equation below: where F j is the measured impact force and A is the measured impact area that is equal to A=π/4d j 2 . To verify the model-calculated velocities with the experimental data under the corresponding conditions, the initial peak velocity of water (U) used for simulation (i.e. in equation (1)) needs to be determined according to the experimental parameters. The initial peak water velocity can be expressed as The average water velocity, V avg,w , at the nozzle exit in equation (25) can be found using the Bernoulli equation, i.e.

( )
where χ d is a discharge factor to account for velocity loss during jet formation. The water pressure (P) used was the same as that in the experiment and the discharge factor χ d is taken as 0.8 according to Sheikh-Ahmad [49]. A comparison of the waterjet velocities obtained from the CFD model and the experiment is shown figure 4. It can be seen that the calculated velocities agree well with the corresponding experimental data both qualitatively and quantitatively. It is shown that the velocity decreases very marginally with an increase in downstream distance from the nozzle exit and increases with an increase in water pressure. These trends are consistent with previous findings of Liu [30] and Hou et al [48]. A quantitative assessment is made based on the percentage deviations of the model predictions from the experimental data under the corresponding conditions. The percentage deviation of the jet velocity is found to be 6.80% on average with a standard deviation of 4.15%.

Particle velocity
In order to measure the particle velocities inside the AWJ flow, the technique combining the use of PIV and LIF was applied and the measurement was made in the jet normal flow zone away from the impact region for the reasons discussed earlier. This technique is capable of obtaining particle trajectory and velocity profile in an AWJ [50]. In the PIV system, a laser light sheet penetrates the whole jet to excite the abrasive particles inside the jet. Due to the effect of water droplets and water mist, it is difficult to capture the particles inside the waterjet, so that the LIF technique was employed [50][51][52]. In the LIF technique, the abrasive particles are coated with fluorescent powder, Rhodamine B, so that the coated particles can absorb the laser light and emit higher light wavelengths. A high-pass filter was attached to the camera lens to eliminate laser light and transmit only the light from the coated particles for clearer images. Figure 5 shows the experimental setup. The PIV system used consisted of a PIV synchronizer (from ILA GmbH), a Gemini Nd:YAG laser and a CCD (charge coupled device) SensiCam camera. In this experiment, an in-house developed water jetting system that could generate a water pressure of up to 67 MPa was used to form the AWJ and the nozzle diameter was 0.762 mm. Three levels of water pressures, 10, 15 and 20 MPa, were considered. Garnet particles of 80 mesh (average diameter of 0.18 mm and density of 4.1 g cm −3 ) were used with a mass flow rate of 0.3 g s −1 through an abrasive supply tube. The measurements of particle velocities were taken at up to 40 mm downstream from the nozzle exit. The experiment for each water pressure was repeated a number of times so that at least 200 particle image pairs were taken.
When measuring particle velocities using a PIV, the laser emits two laser pulses with a very short time interval, while the camera acquires two images corresponding to the two laser pulses for the same particle that moves in the flow. Figure 6(a) shows some of the image pairs (called the 1st and 2nd images) in this experiment. By knowing the time interval of the two laser pulses and the distance between the two particle locations in the images based on the calibrated scaling factor of the camera, the particle velocity can be evaluated. In this study, the time interval between the two laser pulses was set at 20 μs. A detailed description of the PIV and LIF measurement technique can be found in [53]. The particle displacement increases with an increase in the water pressure as shown in figure 6(a). As shown in figure 6(b), particles also rotated while moving downstream, which might be induced by the frictional force from the surrounding water and the irregular shapes of the particles.
For this model, the particle velocity profile at the inlet boundary (nozzle exit) is determined from the waterjet velocity profile (equation (1)), where the peak particle velocity, U, is defined by equation (25). From equation (25), the average particle velocity, V avg,p , at the nozzle exit can be found using the momentum transfer equation, i.e. where ψ is the momentum transfer coefficient. The coefficients of ψ and χ d are taken as 0.8 and 0.9, respectively, according to Sheikh-Ahmad [49]. The abrasive mass flow rate,  m a , used has the same value as in the experiment work presented in the last section. Figure 7 shows a comparison of the model-calculated and experimental particle velocities, in which the symbol represents experimental data and the trend line shows the calculated results. The calculated velocities from the CFD model agree well with the experimental results both at the jet centerline (r=0) and the jet edge (r=0.38 mm). With an increase in the downstream distance from the nozzle exit, the particle velocity at the jet centerline decreases ( figure 7(a)), but at the jet edge increases ( figure 7(b)). This trend follows what was reported by Liu et al [9] due to momentum exchange within the jet. The percentage deviation of the calculated particle velocity on the jet centerline from the corresponding experimental data is found to be 3.70% on average with a standard deviation of 2.78%, while the percentage deviation on the jet edge is 8.77% on average with a standard deviation of 9.90%. Although velocity measurements under ultrahigh pressure conditions were not possible due to the institution's strict health and safety policy which did not allow us to use PIV in the waterjet laboratory, it may be stated that the CFD model is adequate to evaluate particle velocities under the tested conditions.

Jet impact characteristics
The phenomenon when a waterjet impacts a target surface at a 3 mm standoff distance is shown in figure 8 with respect to the time of computation, in which both the water volume fraction and water velocity contours are given. At an early stage (t=2 μs), a jet core is generated with a jet tip at the end of the jet flow, which is induced by the resistance force from the atmospheric pressure ( figure 8(a)). In addition, water droplets are developed at the edge of the jet tip due to the surface tension between the water and air. As the calculation proceeds, the jet flows further downstream and impinges the target surface at about t=6 μs. After the jet reaches the target surface, it rebounds and results in water droplets and water hydraulic jumps as shown in figure 8(a) when t=8 and 10 μs. For the velocity contours illustrated in figure 8(b), a stagnant layer appears to have been developed at the impact site on the target surface and its size increases with an increase in the impinging duration (from t=6 to 10 μs). This finding is consistent with the simulation results by Schwartzentruber et al [21].
Stagnation zone thickness is determined with reference to the jet velocity profile at the downstream location where the jet velocities reduced sharply and significantly with a small variation across the jet. There appears to be a clear separation between the jet flow velocity profile and the jet velocity profile at stagnation zone that enables to identify the stagnation zone. For the case shown in figure 9, the stagnation thickness is about 0.4 mm above the target surface. As the stagnation formation  affects the jet impact force transferred to the target, the finding that stagnation zone dimension increases with the impinging duration may be used to explain why the jet impingement at the initial stage has a higher impact energy [24]. As a wall jet has been fully developed on the target surface by the impact time of t = 100 μs, the jet flow is steady and, hence, a stable stagnation formation, as shown in figure 8.
A comparison of the jet flow and jet impact on a target surface at a 3 mm standoff distance is shown in figure 9(a). In the impact case, the jet structure may be classified to have three main regions; namely a free jet flow region, a jet deflection region with a stagnation zone and a wall jet region. In the free jet region, a power law velocity profile across the jet applies, and a flat cap profile evolves as jet flows downstream [9]. It is also found that the jet flow does not spread significantly in the free jet region within the small standoff distance considered. In the deflection or stagnation region, the jet rebounds from the wall (or target surface), resulting in the formation of a stagnant layer at the impact site ( figure 9(a)). As such, the jet flow expands and the jet velocity decreases sharply when it enters the stagnation zone. The axial flow velocity drops immediately due to the reflection force from the target surface, while the radial velocity increases along the target surface due to the change of flow direction ( figure 9(b)). In the wall jet region, the flow is forced to travel parallel to the wall (or target) surface, in which the velocity component normal to the surface is much less than the tangential component.
Stagnation formation may be explained by the compression waves reflecting from a flat surface [24]. It is formed  according to the Newton's third law when the jet hits a target surface [44]. The portion of reflection force after impacting a target surface is decreased due to the free phase flow of waterjet that is forced to change direction by the rigid wall, leading to a decrease in the axial velocity, but an increase in the radial velocity. Stagnation formation is also affected significantly by the shear stress on the wall and the fluid viscosity [54]. According to the assumption of no-slip wall, the shear stress or friction force is imposed on the wall, resulting in a damping flow next to the stagnation zone before the jet flows parallel to the wall surface ( figure 9(a)). This phenomenon has been explained with the logarithmic law of the wall [54]. Figure 10(a) shows the velocity profiles in the jet axial direction at different radial locations for a 3 mm standoff distance. It is found that the jet velocity decreases as the jet downstream distance increases, and the decreasing rate reduces as the flow is closer to the jet rim (e.g. r=0.3 mm). Interestingly, after flowing for about 2.6 or 0.4 mm from the target surface, the jet velocities across the jet merge into almost the same value and decrease sharply ( figure 10(a)). The variation of jet velocity along the jet radial direction at different jet downstream locations is shown in figure 10(b). In order to identify and study the damping effect of the stagnation zone, the velocity profile during an AWJ impact process is compared with the velocity profile calculated by the AWJ flow model at the respective jet downstream locations. It should be noted that these velocity profiles may go outside the jet flow field and include the surrounding air that is entrained by the jet to gain a velocity. Thus, the water flow referred to here may be a combination of water and air flow.
From figure 10(b), it is found that the greatest jet reflection is generated at close to the target surface (i.e. z= 3 mm), where the jet velocity is close to zero at the jet centreline. As shown at z=2.7 in figure 10(b), the stagnation remains influencing the jet impingement, so that the minimum velocity occurs at the jet centre and the velocity increases with an increase in radial distance. The velocity across the jet field is almost constant at z=2.6 mm, and the velocity profile from the AWJ impact model is similar to that from the AWJ flow model at z=2.3 and 2.5 mm ( figure 10(b)). As the jet flows further downstream, the velocity profile evolves to a flatter cap shape. However, when the jet arrives at the target surface at z=3.0 mm, the jet bounces back and has a totally different profile. From figure 10(b), stagnation may start to affect the jet flow at about 0.4 mm above the target surface, which may be used to define the thickness of the stagnation zone under this condition.

Effects of impact parameters on stagnation formation
As discussed above, stagnation affects the jet energy transferring to the target and greatly changes the jet flow direction when the jet approaches the target surface. Figure 11 shows the effects of impact parameters, i.e. peak inlet velocity, standoff distance, nozzle diameter and jet impact angle, on the stagnation formations. It should be noted that the velocity profiles in the case of perpendicular jet impact have an axisymmetric distribution. Thus, the characteristics of stagnation formation is considered only in the r-z plane. Figure 11(a) shows the effect of peak inlet velocity (U=500, 600, 700 and 800 m s −1 ) at the nozzle diameter d n =0.76 mm and standoff S d =3 mm. It is found that the difference in the stagnation thickness under the four peak inlet velocities is insignificant. Three standoff distances (S d =3, 5 and 10 mm) and three nozzle diameters (d n =0.76, 1.00 and 1.20 mm) are considered at the condition of U=800 m s −1 , as shown in figures 11(b) and (c). It can be noticed that a higher standoff distance or a larger nozzle diameter causes a larger stagnation zone. This finding agrees with the study of Hammad and Milanovic [55]. Figure 11(d) shows the velocity contour with the variation of jet impact angle, where jet impact angle is defined as the angle between the initial jet flowing direction and the target surface. Thus, the perpendicular impact angle is when θ=90°. In this study, three jet impact angles (θ=70°, 80°and 90°) are considered. It is shown in figure 11(d) that the stagnation zone formed is asymmetric to the jet axis under oblique impact conditions, where the stagnation centre is away from the jet centreline. A decrease in jet impact angle causes the stagnation centre to move further away from the jet centreline and a stronger tangential flow along the radial direction in the downhill side parallel to the target surface. Further, at the downhill side, the stagnation layer is smaller and vice versa ( figure 11(d)).
According to Bernoulli's equation, the jet velocity approaches zero when it impacts at a solid surface, resulting in only static pressure on the target surface [44]. As the static pressure represents the potential energy of the jet, it can also describe the jet impact energy. Figure 12 shows the effect of impact parameters on the static pressure on the target surface at different jet radial locations. It is found that the static pressure increases with an increase in the peak inlet velocity and nozzle diameter (figures 12(a) and (c)), but decreases with an increase in standoff distances ( figure 12(b)). As discussed earlier, a larger stagnation zone is formed at a larger standoff distance, resulting in less static pressure. A larger nozzle diameter increases the static pressure because a larger jet impact area and larger jet flow volume. The effect of jet impact angle on the static pressure is shown in figure 12(d). Similar to the velocity profile, it is found that the static pressure profile in the jet radial direction is symmetric about the jet axis at a normal impact (θ=90°), but not so for an oblique impact (θ<90°) where the maximum pressure is not at the jet centreline. In addition, a reduction in the impact angle decreases the maximum static pressure, probably due to the inclined target surface that facilitates the fluid flow and hence generates a smaller static pressure.
For a perpendicular impact, a centrosymmetric impingement area is generated, where a stagnation centre at the jet centreline (S s =0), as seen in figure 13(a). By contrast, an oblique impact causes an ellipse impingement area on the target surface, with a major axis of 2r 0 /sinf and a minor axis of 2r 0 , as shown in figure 13(b), in which the stagnation centre is no longer coincident with the jet centre, but is away at a distance of S s =r 0 cot(θ) on the major axis.
From figure 13, the effective jet impingement area can be determined based on the polar radius and the stagnation centre as [56]: where r s is the impingement radius for both perpendicular and oblique impacts, a j is the jet expansion coefficient, r 0 is the jet radius at the nozzle exit and f is the angle of an arc in a circle or an ellipse, as shown in figure 13. Figure 14 shows the particle impact characteristics during a normal impact. It can be seen that particles flow initially in their trajectory in the free jet region and are then decelerated when going through the stagnation zone with an increased decay rate as approaching the target surface. It has been noticed that particles remain at almost the same radial location when they are in the free jet flow region. This agrees with the experimental results of Balz's work [57] in which they found that the radial velocity component of the particles was approximately 0%-4%, of the axial velocity. It is also shown that particles move further away from the jet centreline as they approach the target surface because of the jet divergence induced by the stagnation zone, resulting in a larger jet impact area and decreased particle impact angles (α p ). After an impact on the target surface, as denoted by p a , the particle rebounds and may move outside of the jet domain or inside the main jet flow field to impact the target again. However, the subsequent impact (denoted by p b ) may be associated with a much lower velocity not to cause any considerable effect even if they may regain some energy from the subsequent flow.

Particle impact characteristics
From figures 14(a) and (b), it is noticed that the characteristics of the rebounding particles are different between with and without considering energy loss during an impact. In the case that the energy loss is not considered ( figure 14(a)), some particles may turn back for another impact at the edge of the stagnation zone or outside the main jet impact site. By contrast, when considering energy loss ( figure 14(b)), most of the particles follow the jet streamline along the wall (target surface), rather than flying out of the water flow field, because a small particle velocity is retained after colliding the target surface. It should be noted that in this study, the target surface is assumed to be absolutely solid with no deformation by any impact. A further study may be undertaken to examine the particle impact characteristics when erosion or material removal takes places during particle impacts. Figure 15(a) shows the variation of water and particle velocities along the jet axial direction at the radial locations of r=0 and 0.3 mm. It is found that the particle velocity decreases slightly with an increase in downstream distance, but drops abruptly when approaching the target surface. The water and particle velocities have a similar trend, but the particle velocities are associated with a lower decay rate due to their greater mass and momentum [30]. The variation of particle velocities crossing the jet flow field at different downstream locations is further studied, as shown in figure 15(b). It is apparent that particle velocity has the highest value at the jet centreline and the velocity decreases with an increase in radial distance from the jet centre. The particle velocity is decelerated due to the stagnation effect, and the decay rate of the particle velocity at close to the jet centre is higher than that at the jet edge because of a higher  damping effect of the stagnation zone at the jet centreline. As a result, the particle velocity across the jet domain becomes flatten out as the jet downstream distance is increased ( figure 15(b)).
In the following sections, the effect of impact parameters on the particle impact velocity and impact angle is presented. Because the particle velocity on the target surface is actually zero, as shown in figure 15(a), the particle velocity immediately before an impact on the target surface (at 0.002 mm above the target) is used to represent the particle impact velocity in the model calculations. In addition, the total velocity (i.e. the summation of the axial and radial velocities) is used as the particle impact velocity. 4.3.1. Effects of impact parameters on particle impact velocity. Figure 16 shows the variation of particle impact velocity along the jet radial direction under different peak inlet velocities, standoff distances, nozzle diameters and jet impact angles. It is found that particle impact velocity increases with an increase in peak inlet velocity and has a similar profile for all peak inlet velocities considered ( figure 16(a)). The particle impact velocity profiles under different standoff distances or nozzle diameters are similar, especially at close to the jet centreline, as shown in figures 16(b) and (c). As mentioned earlier, a greater standoff distance or nozzle diameter induces a stronger stagnation effect, leading to a lower particle impact velocity. However, the impact velocity under different standoff  distances has no discernible difference at the jet edge ( figure 16(b)).
As shown in figure 16(c), a larger nozzle diameter causes a more flatten velocity profile, because of the momentum changes as the jet flows downstream. From figure 16(d), an asymmetrical particle impact velocity profile along the jet radial direction is found under the oblique jet impact condition, where the particle velocity at the target downhill side increases with a decrease in the impact angle, while a reverse trend is seen at the target uphill side for the reasons discussed earlier. It is expected that the material removed at the downhill side is greater than that at the uphill side. After impacting the target surface, most particles travel in the downhill direction due to the fluid drag as well as the gravity force, as seen in figure 17.
In order to assess the effect of stagnation zone on the particle impact energy, the ratio of the particle impact velocity to the particle flow velocity without impact at the same downstream location (V p,imp /V p,flow ) is used. The particle impact velocity ratio is plotted against the dimensionless radial distance (or location) with respect to the jet (or nozzle) radius (r/R n ), as shown in figure 18. It can be seen that the particle impact velocity ratio increases when the dimensionless particle location in the jet radial direction increases. It is apparent that the damping effect of the stagnation has its maximum at the jet centre and reduces as the dimensionless radial location of particles increases.
Considering the whole jet cross section and the conditions used for figure 18, the average particle impact velocities is about 67% of the particle velocity if the particles flow to the same location without impact or without stagnation effect (or the stagnation decreases the particle velocities by approximately 33%). From figure 18(a), the particle velocity ratio at close to the jet centreline has no discernible difference under different dimensionless standoff distances. However, the velocity ratio profile becomes flatten out with an increase in dimensionless standoff distance. However, the velocity ratio at close to the jet centerline decreases with an increase in nozzle diameter because of the increase of the stagnation zone size at the jet centre, but the difference of the velocity ratio for all nozzle diameters at the jet edge is not significant, as shown in figure 18(b).
As discussed earlier, it has been found that particles rotate while moving downstream from the nozzle exit, due to the irregular particle shape and jet divergence. The effects of these parameters are now explained in which the particle mass is represented by particle diameter and the particle shape is represented by particle sphericity. Four particle sphericities (0.7, 0.8, 0.9 and 1) and five particle diameters (0.076, 0.12, 0.16, 0.18 and 0.20 mm) are considered, where a round particle has the sphericity, S p , equal to one and a non-round particle has the sphericity smaller than one. Figures 19(a) and (b) shows the variation of particle impact velocity with respect to particle sphericity and particle  diameter, respectively. It is found that a greater particle sphericity rises the particle impact velocity, whereby a round particle (S p =1) has the highest velocity profile. This is due to the fact that a round particle is associated with a lower drag force than non-spherical particles [58], and hence retaining a higher velocity. From figure 19(b), particle impact velocity increases with an increase in particle diameter for smaller particles (d p =0.076-0.120 mm), but it does not show any discernible difference for larger particles (d p =0.16, 0.18 and 0.20 mm), possibly due to their large momentum. Nevertheless, the difference of particle velocities at the close the jet rim is very small.

4.3.2.
Effect of impact parameters on particle impact angle.
As discussed before, particles deviate significantly from its original trajectory as they approach the target surface because of the stagnation effect. As such, the particle impact angle is not equal to the jet impact angle. Figure 20 shows the variation of particle impact angle along the jet radial direction under the variation of other parameters. It is found that under the normal impact mode (figures 20(a)-(c)), particle impact angles exhibit an axisymmetric profile and decrease with an increase in the jet radial distance from the jet centreline. Particle impact angle slightly varies with the change of the standoff distance or peak inlet velocity at close to the jet rim. From figure 20(c), a larger  nozzle diameter yields a larger particle impact angle, hence a larger impact area. It is also noticed that particle impact angle at close to the jet rim changes significantly, with about 30°r eduction from the normal impact angle (90°). For the oblique impacts shown in figure 20(d), there is an asymmetric profile for the particle impact angle along the jet radial direction, in which the peak impact angle moves towards the uphill side (with negative radial distance) as the jet impact angle decreases. From figure 20(d), the decreasing rate of particle impact angle for oblique impacts is smaller than the normal impact. This is because the tangential flow along the target surface changes the stagnation shape and induces the stagnation centre to move away from the jet centre ( figure 11(d)), leading to a weaker stagnation effect. Figure 21 shows the effect of particle diameter (d p ) and sphericity (S p ) on particle impact angle. It can be noticed that particle impact angle decreases with an increase of the particle sphericity or particle diameter, and the decreasing rate increases along the jet radius outwards. This is due to the fact that the particles with a smaller sphericity or diameter are expected to more easily deviated from their trajectory by jet divergence, resulting in a reduction in the impact angle. Particles at the jet centreline is expected to received a balanced force in the jet radial direction so that their flow trajectory is not affected by particle size or sphericity.

Conclusions
A study of the impact characteristics of ultrahigh pressure AWJ on a solid surface has been presented. A CFD model has been developed and verified by the experimentally acquired velocities of water and particles. It has been found that the structure of an impacting waterjet can be divided into three major regions; namely a free jet flow region, a stagnation region, and a wall jet region. The stagnantion formed has a major effect on the jet and particle kinetic energy that is transferred to the target, and the strongest effect is found at the centre of impact area for normal impacts. The impact parameters have been found to affect stagnant formation, whereby an increase in the standoff distance or nozzle diameter results in an increase in the dimension of the stagnation zone, hence a reduction in the jet impact velocity. An oblique jet impact yields an asymmetry stagnation zone, where the stagnation centre is away from jet centreline. Using the static pressure to represent the jet impact energy, it has been found that the static pressure increases as the peak inlet velocity or nozzle diameter increases, or the standoff distance decreases. An oblique impact causes the highest static pressure to move away from the jet centreline. Particle velocities decrease abruptly when going through the stagnation zone. The decay rate at close to the jet centre is higher than that at the jet edge because of the higher damping effect of the stagnation zone at the jet centre. It is found that the average particle velocity across the jet flow is reduced to about two third of its value after going through the stagnation zone. The particles also diverge from the jet centreline when they approach the target surface, resulting in a change in the actual impact position and angle. It has also been noticed that a larger particle diameter or a greater particle sphericity yields an increase in particle impact velocity, particularly at locations close to the jet centreline, and in particle impact angles.