A New Method for Predicting Erosion Damage of Suddenly Contracted Pipe Impacted by Particle Cluster via CFD-DEM

A numerical study on the erosion of particle clusters in an abrupt pipe was conducted by means of the combined computational fluid dynamics (CFD) and discrete element methods (DEM). Furthermore, a particle-wall extrusion model and a criterion for judging particle collision interference were developed to classify and calculate the erosion rate caused by different interparticle collision mechanisms in a cluster. Meanwhile, a full-scale pipe flow experiment was conducted to confirm the effect of a particle cluster on the erosion rate and to verify the calculated results. The reducing wall was made of super 13Cr stainless steel materials and the round ceramsite as an impact particle was 0.65 mm in diameter and 1850 kg/m3 in density. The results included an erosion depth, particle-wall contact parameters, and a velocity decay rate of colliding particles along the radial direction at the target surface. Subsequently, the effect of interparticle collision mechanisms on particle cluster erosion was discussed. The calculated results demonstrate that collision interference between particles during one cluster impact was more likely to appear on the surface with large particle impact angles. This collision process between the rebounded particles and the following particles not only consumed the kinetic energy but also changed the impact angle of the following particles.


Introduction
Erosion prediction is usually carried out using metal walls to evaluate the safety of process equipment, such as sudden contraction pipes. During prediction, numerous parameters, such as fluid flow characteristics, particle characteristics, particle impingement information, and target surface properties, jointly affect the efficiency and accuracy of erosion calculations, which have been widely studied and have accounted for a large number of prediction models [1]. However, almost all of these models are based on parameters or statistics of individual particles and the erosion behavior of particle clusters has received limited attention. When a particle cluster impacts the wall, the particles that initially make contact with the wall bounce back and collide with subsequent particles in a phenomenon called collision interference. This oncoming collision can decrease the impact energy and change the impact angle of the following particles, thereby causing a large error in erosion prediction.
Erosion prediction in complex geometries adopts the computational fluid dynamics (CFD) method. CFD generally consists of three steps: Flow modeling, particle tracking, and erosion calculation. A flow field calculation is typically solved based on Navier-Stokes equations and a set of conditions, such as transient or steady state, compressible or incompressible, and laminar or turbulent flow. When the

Solution Methodology
DEM takes into account the normal deformation, tangential slip, and rotation of each collision particle in solving the motion equation by using a hard-sphere or a soft-sphere model. In this study, we present the governing equations of fluid flow and particle motion; discuss the determination of particle interference for near-wall particle motion and provide a brief description of the erosion model.

Particle-Particle Interaction
The first step for getting the evolution in time of the (translational and rotational) velocity and position of the particles is to calculate the particle translation and rotational motions according to Newton's second law [18]. The total force and torque acting on a solid particle can be expressed by [19].
The total force acting on a particle is calculated as a sum of the total contact force (F ij ), fluid interaction force (F l,i ), gravitational force (F g ) and buoyancy (F b ). The total contact force is split into two components: Normal (F n,ij ) and tangential (F t,ij ). The total torque T i results from a vector summation of the torque at each particle-particle contact [20].
For particle-particle contact, we adopted the Hertz-Mindlin (no slip) contact model to calculate the inter-particle contact force and damping force. A possible contact model in the system is given in Figure 1. Each force and moment can be treated as a spring or a damper and then, the model can be described by four elastic-damping equations. In this elastic-damping model, the normal force component is based on Hertzian contact theory and the tangential force model is based on the work of Mindlin-Deresiewicz [21]. This model has been used for contact parameter estimation for the collision of iron ore balls [22], the collision of non-spherical particles in a vibrating bed [23], and collision of particles in liquid-solid flow [6]. The results showed that this model has great calculation precision and is applicable to the calculation of particle collisions in a liquid. The total force between two particles, which contain a contact force and damping force, can be expressed by where the tangential contact force (F t ) is limited by Coulomb's law of friction and meets the relationship F t ≤ −µ s F n,ij δ t /|δ t | in a system without sliding [18,22]. Therefore, the normal and tangential force based on the normal and tangential overlap can be represented as −S n δ n and −S t δ t , respectively. The normal and tangential damping forces are described in Tsuji's research [24] and can be expressed as C n v rel n and C t v rel t . Therefore, the total force is given by F n,ij + F t,ij = −S n δ n + C n v rel n + S t δ t + C t v rel t (4) where the stiffness constants (S n and S t ), damping coefficients (C n and C t ), and other parameters are defined in Table 1. In addition to the contact force between the particles, there may be some non-contact forces between the particles, which can determine the particle movement state. Van der Waals forces, capillary forces and electrostatic forces are the major non-contact forces between the particles. For particles with a micrometer or nanometer diameter, Van der Waals forces may have a great effect on the movement of the particles but the effect on the particles (d p > 0.6 mm) in this study was ignored [25]. The reason for neglecting the other forces won't be discussed here owing to the limited space.

Particle-Particle Interaction
The first step for getting the evolution in time of the (translational and rotational) velocity and position of the particles is to calculate the particle translation and rotational motions according to Newton's second law [18]. The total force and torque acting on a solid particle can be expressed by [19]. The total force acting on a particle is calculated as a sum of the total contact force (Fij), fluid interaction force (Fl,i), gravitational force (Fg) and buoyancy (Fb). The total contact force is split into two components: Normal (Fn,ij) and tangential (Ft,ij). The total torque Ti results from a vector summation of the torque at each particle-particle contact [20]. For particle-particle contact, we adopted the Hertz-Mindlin (no slip) contact model to calculate the inter-particle contact force and damping force. A possible contact model in the system is given in Figure 1. Each force and moment can be treated as a spring or a damper and then, the model can be described by four elastic-damping equations. In this elastic-damping model, the normal force component is based on Hertzian contact theory and the tangential force model is based on the work

Parameter Expression
Normal stiffness constant (S n ) The relative velocities of the centers of the spheres before and after a collision (∆v) The relative velocity of the contact points v rel = ∆v + d p After we set the numerous parameters in the DEM model, the particle collision system can be well-established without considering the fluid action. However, for the particle clusters movement in an actual structure, fluid force actions on the particles play a key role in particle motion and collision. The fluid force F l,i often includes the drag force (F d ), virtual mass force (F v ), pressure gradient force (F p ), Saffman lift force (F s ), Basset force (F ba ), and Magnus lift force (F m ), respectively. The drag force plays a major role in the force acting on the particles by a fluid, which has been proven by literature [26]. The drag force model can be expressed by where the following equation for the drag coefficient (C d ) is obtained where the particle Reynolds number Re p is defined as In the Lagrangian description, the Saffman lift force and Magnus lift force are not considered because the particle diameter was close to 0.6 mm rather than sub-micrometer. In addition, the pressure varied very slightly over a distance of one particle diameter due to the reasonably small particles. For this reason, we neglected the pressure gradient force. Virtual mass force is not considered due to the small relative velocity between the viscous fluid and particles.
Particle collisions involve energy absorption and dissipation through momentum exchange, friction and impact restitution. It is an important factor when calculating the erosion rate using the energy method. Figure 2a shows the process of collision between two particles. The subscripts (1) and (2) denote the pre-collisional and post-collisional velocities, respectively. The difference in momentum (v) and angular momentum (ω) between two instants of time is equal to the impulsive force acting on the particle, which can be expressed by [27]

Momentum Exchange and Energy Dissipation
Materials 2018, 11, x FOR PEER REVIEW 5 of 24 Figure 2. The sketch of impingement between a particle and the wall, (a) the particle-particle collision process; (b) the particle-wall contact process.
Equation (8a,b) are not sufficient to determine the relationship between the pre-and postcollisional velocities. Auxiliary equations, which contain the restitution coefficient or friction coefficient, are necessary to close a set of equations. Walton [28] has proposed a hard sphere analogy of the soft-sphere model, which can be used if only two particles are in contact at one time. Jenkins [29] followed Walton and distinguished between two kinds of binary collisions: Sliding and sticking collisions. A nearly elastic collision was studied by using the constant coefficient of normal restitution and the tangential restitution to relate the translational and rotational velocities. Here we assume that all particles have the same geometry and physical properties and ignore the liquid forces. The momentum exchange between two particles can be expressed as for the tangential direction. (9b) When the relative velocity of the contact points is substituted into the formula, the condition then becomes According to Walton's research, the post-collisional translational velocities are Figure 2. The sketch of impingement between a particle and the wall, (a) the particle-particle collision process; (b) the particle-wall contact process. Equation (8a,b) are not sufficient to determine the relationship between the pre-and post-collisional velocities. Auxiliary equations, which contain the restitution coefficient or friction coefficient, are necessary to close a set of equations. Walton [28] has proposed a hard sphere analogy of the soft-sphere model, which can be used if only two particles are in contact at one time. Jenkins [29] followed Walton and distinguished between two kinds of binary collisions: Sliding and sticking collisions. A nearly elastic collision was studied by using the constant coefficient of normal restitution and the tangential restitution to relate the translational and rotational velocities. Here we assume that all particles have the same geometry and physical properties and ignore the liquid forces. The momentum exchange between two particles can be expressed as When the relative velocity of the contact points is substituted into the formula, the condition then becomes v rel(2) n · r = −e n v rel(1) n · r (10a) According to Walton's research, the post-collisional translational velocities are v and the rotational velocities are where the expressions of parameter ∆v, v rel and K are shown in Table 1. There is an angle between v and r that can determine whether a collision is sliding or sticking. When this angle is greater than or equal to a critical angle θ*, a sliding collision occurs. Otherwise, a sticking collision takes place. The tangential coefficient of restitution can be expressed by [30] This hard sphere analogy of the soft-sphere collision dynamics greatly reduces the computational time for a soft-sphere model. When the particle cluster impacts on the wall, many particles will bounce back from the wall and collide with others. Therefore, the kinetic energy of particle j may be reduced by the collision of particle i before it impacts the wall. We define the ratio of pre-collisional velocity to post-collisional velocity as the decay coefficient, which includes the normal momentum decay coefficient (λ n ), the tangential momentum decay coefficient (λ t ) and the energy decay coefficient (λ e ). The following definitions have been used in general.
When the particle-particle collision is considered as a nearly elastic one, according to the law of conservation of energy, the exchange of kinetic energy between two particles can be expressed by the following relation 1 2 mv where the change of translational and rotational kinetic energy have the forms of [18]  If the energy exchange is expressed in the form of a reduced velocity, i.e., 1 2 m p (∆v ) 2 = ∆E n + ∆E t , the decrease in velocity of particle j can be expressed by an energy decay coefficient, which is

Particle-Wall Interaction
In order to determine whether there is a collision between two particles near the wall, the particle-wall contact is taken into account as shown in Figure 2b. Unlike particle-particle contact models, particles are treated as fully rigid bodies when they come into contact with the wall. Comparing the rebound time, t 1 , for particle i and the moving time t 2 for particle j, if t 1 > t 2 , the particle j will collide with particle i before it leaves the wall region, on the contrary, it will not. The rebound process for particle i consists of two steps including a contact step and a release step. Therefore, the rebound time can be expressed as t 1 = t e + t r . When the effect of elastic deformation is ignored, the depth of the extruding lips is given by The strain caused by the normal stress of indentation can be calculated by According to Hooke's law, the average force acting on a metal surface during the indentation process is where the maximum depth of indentation, H, can be described as [31] According to the mutual principle of forces, the normal force for a single particle is equal to the force of a metal surface. When the indentation depth reaches a maximum, the particle's normal velocity, v i,z , closes to 0 m/s. Therefore, the contact time with the metal surface, based on the theorems of particle momentum, can be expressed by the following equation.
In the process of contact between the particle i and the wall, the moving distance of the particle j along the centerline is According to Habib's research [32,33], particles move toward the centerline of the pipe as they approach the sudden contraction section. As shown in Figure 2b, the movement of two particles can be treated as a meeting problem in the Z direction as well as a catching-up problem in the Y direction after the rebound of particle i. Therefore, the release time (t r ) should meet the following relationship in order to satisfy the time adaptive condition of particle collision. (24) If the distance between the two particles is far greater than the particle diameter and if the acceleration of particles is ignored, the critical inter-particle distance in the release process will be given by where i,z ) Then, the critical inter-particle distance in the process of a particle-particle collision near the wall can be expressed by where the release velocity for particle i can be obtained using the velocity component based on the 3 The velocity loss rate for a single particle in the process of particle-wall contact is given by

Governing Equations of Fluid Flow
For the liquid phase, the governing Equations comply with the law of conservation of mass and momentum in terms of local-average variables. Taking the volume fraction of the continuous phase and the particle-liquid interaction force into account, the continuous carrier liquid flow equations are based on Model A [35] and are given by where the particle force acting on the liquid (F p,i ) should be equal to the force of the liquid phase (F l,i ) acting on the solid phase in the opposite direction. In the above equations, the fluid volume fraction is obtained from the relation α l 1 − N ∑ j=1 V p /∆V l and the viscous stress tensor τ is defined as In the process of the coupling calculation, the motion of a particles phase is obtained by DEM, which applies Newton's laws of motion at the individual particle level, while the flow of continuum fluid is described by the local averaged Navier-Stokes equations that can be solved using the CFD approach.

Erosion Model
A detailed literature survey by Meng and Ludema [36] revealed that more than 30 erosion models exist for particle impact erosion to date. All prediction models basically include the particle impact velocity and impact angle, such as the model proposed by Chen [11], McLaury [12] and Ahlert [13]. According to Chen's research, erosion can be calculated by using Equation (30) where X and Y represent the percentage of grouped particles and independent particles in the total impact particles, respectively. Therefore, the total erosion rate is controlled by multiple impacts of different kinds of particles. The detailed erosion models for each impact are given by [11] where λ = λ n ·λ t ·λ e ; F s = 1.0 for sharp, 0.53 for semi-rounded, or 0.2 for fully rounded sand particles.
The impact angle function f (α) is given by The erosion properties of the target wall material (super 13Cr stainless steel) have been studied by a jet flow system at the Xi'an Shiyou University [37] and at a critical impact angle, θ = 70 • . The new erosion sub-model parameters used in the erosion calculation are shown in Table 2. In the post-processor, the erosion was visualized as the local wall thickness loss rate on different radial positions of the reducing wall.

Numerical Solving Step
In the current study, erosion prediction involved the calculation of liquid-particle, particle-particle, and particle-wall interactions, which were completed via CFD coupled with DEM. The CFD code ANSYS FLUENT [38] was employed to model the continuous phase flow. Meanwhile, discrete particle motion modeling was accomplished using the DEM code EDEM [39] as shown in Figure 3. The governing equations (Equation (28a,b)) were discretized in a finite volume form and then, the discretized equations were solved using a semi-implicit method for pressure linked equations (SIMPLE), which is described by Ferziger [40] in detail. The explicit time integration method was used to calculate the motion of the particles. When a stable fluid region was obtained, the drag and contract forces between the fluid and particle were calculated using the DEM time step. Subsequently, new particle positions were transferred to the coupling module between CFD and DEM. The flow field was updated until the results satisfied the accuracy requirement. Two calculations were performed independently but were coupled at regular intervals, commonly with multiple DEM time steps for a single CFD time step.

Boundary Conditions
The detailed parameters for fluid, particle, and geometry are summarized in Table 3. As shown in Figure 4, the particles were generated on two layers (Particle Layers I and J) of virtual quadrilateral planes, with the two planes radially distributed and spaced from each other by 0.05 mm. Each virtual plane was10 mm in length, 3 mm in width, perpendicular to the pipe axis and close to the pipe wall. One of the purposes of this setting was to induce the maximum number of particles to impact the reducing wall because only particles near the wall of the upstream pipe can impact the wall. Another purpose was to distribute the particles as far as possible in the radial direction of the target wall.
Particle Layer I, which was in front of Particle Layer J, represented particles that initially impacted the wall and were not subject to particle collisions. By contrast, the particles in Layer J may suffer from particle-particle collision interference near the wall. In each particle layer, particles were randomly positioned, thus satisfying the multiple possibilities of particles impacting the wall. Therefore, the number of particles mentioned in Table 3 is an average value and may change with different particle-producing positions. Meanwhile, a particle-moving calculation at each velocity was repeated 10 times to count the particle impact parameters on different radial walls because not all surfaces in the calculation process were impacted. Moreover, the inducing wall was considered to be under a no-slip condition, where ur|r = ux|r = 0 = 0 and әur/әx = 0 and a fully developed incompressible flow was considered at the velocity inlet and exit sections. In the coupling process, space and time steps of CFD and DEM must satisfy a certain correspondence. Although two phases were created in Fluent using the Eulerian-Eulerian model, conservation equations for the solid phase were not solved. Therefore, it was necessary to calculate the volume fraction of the CFD mesh cell occupied by the particles. In order to obtain the position of the particles before and after the collision, the CFD mesh cell occupied by the particles was usually divided into a number of small spaces. The particle volume fraction within a particular mesh cell was, therefore, the percentage of the points as given by The grid size, ∆x, which can be expressed by ∆x = (V p /N p ) 1/3 , determined the number of CFD points within the mesh cell of a particle. It was usually necessary to ensure that N p Materials 2018, 11, x FOR PEER REVIEW 10 of 24 the particles before and after the collision, the CFD mesh cell occupied by the particles was usually divided into a number of small spaces. The particle volume fraction within a particular mesh cell was, therefore, the percentage of the points as given by The grid size, ∆x, which can be expressed by ∆x = (Vp/Np) 1/3 , determined the number of CFD points within the mesh cell of a particle. It was usually necessary to ensure that Np ≧ 1 and Np < Nl. Meanwhile, a grid size that was too small would increase the amount of calculation and a grid size that was too large size may ignore collision details. Therefore, the grid size was set to 0.2 mm for the CFD mesh in this work.
In general, the CFD calculation time size and the DEM calculation time size must satisfy the following relationship: Firstly, the minimum time step size needed to satisfy the CFD calculation convergence condition. Secondly, the ratio of the DEM time step to the Rayleigh time step should be in the range of 5-30%. Thirdly, the ratio of the CFD time step to the DEM time step must be an integer. The Rayleigh time step between two particles is given by [41] 877 Studies have shown that the DEM time step is proportional to the Rayleigh time step and the ratio is within the range of 0.1-0.5 [42][43][44]. Here, the DEM time step was equal to one-fifth of the Rayleigh time step, that is, ΔtD = 0.2 × ∆tRa and the shear modulus can be expressed as G = Y/2(2 + ν)(1 − ν). In order to obtain the contact parameters between the particles and the wall, the DEM time step needs less than the contact time (Equation (22)). Therefore, the time step was set to be 2 × 10 −6 s for DEM and 4 × 10 −5 s for CFD.
When the particles move closer to the reducing wall, they will change direction and accelerate 1 and N p < N l . Meanwhile, a grid size that was too small would increase the amount of calculation and a grid size that was too large size may ignore collision details. Therefore, the grid size was set to 0.2 mm for the CFD mesh in this work.
In general, the CFD calculation time size and the DEM calculation time size must satisfy the following relationship: Firstly, the minimum time step size needed to satisfy the CFD calculation convergence condition. Secondly, the ratio of the DEM time step to the Rayleigh time step should be in the range of 5-30%. Thirdly, the ratio of the CFD time step to the DEM time step must be an integer. The Rayleigh time step between two particles is given by [41] ∆t Ra = πr p ρ p /G 0.163ν + 0.877 (34) Studies have shown that the DEM time step is proportional to the Rayleigh time step and the ratio is within the range of 0.1-0.5 [42][43][44]. Here, the DEM time step was equal to one-fifth of the Rayleigh time step, that is, ∆t D = 0.2 × ∆t Ra and the shear modulus can be expressed as G = Y/2(2 + ν)(1 − ν). In order to obtain the contact parameters between the particles and the wall, the DEM time step needs less than the contact time (Equation (22)). Therefore, the time step was set to be 2 × 10 −6 s for DEM and 4 × 10 −5 s for CFD.
When the particles move closer to the reducing wall, they will change direction and accelerate under the action of the liquid force [32,33], resulting in an increase in the distance between the particles. Therefore, whether the particles in a cluster collide with each other near the target and thus affect the erosion prediction is determined by the following steps: Calculation of the critical interparticle distance: The distance between adjacent particles was calculated when the front particle made contact with the wall. If L < L p , the subsequent impact on the wall was treated as a clustered particle impact; otherwise, the impact was considered as an independent particle impact.
Acquisition of clustered particles: The distance determination method was used to divide the form of particle impingement on the wall and impact parameters were calculated by classification.
Determination of the decay coefficient: For the clustered particles, the next step was to obtain the decay coefficient λ and the new impact angle for each colliding particle after a particle-particle collision.

Boundary Conditions
The detailed parameters for fluid, particle, and geometry are summarized in Table 3. As shown in Figure 4, the particles were generated on two layers (Particle Layers I and J) of virtual quadrilateral planes, with the two planes radially distributed and spaced from each other by 0.05 mm. Each virtual plane was10 mm in length, 3 mm in width, perpendicular to the pipe axis and close to the pipe wall. One of the purposes of this setting was to induce the maximum number of particles to impact the reducing wall because only particles near the wall of the upstream pipe can impact the wall. Another purpose was to distribute the particles as far as possible in the radial direction of the target wall.  (c) the location of a particle cluster.

Full-Scale Experiment
A full-scale pipe flow experiment was used to study particle clusters during erosion and to verify the simulation results. As shown in Figure 5, the experimental loop mainly comprised a contraction section, a screw pump (with a flow range of 5-60 m 3 /h), an electric heating agitator, a temperature and pressure sensor, a magnetic flowmeter (8712HR, Emerson, Rosemount. Co., USA), four gate values, a control cabinet, a computer, and connecting pipes. After the test fluid containing the particles was mixed, stirring and electric heating units were turned on until the temperature reached the required value and then the pump was turned on. The gate valve in the test pipe was opened when the flow reached stability and the related data, including a particle motion image, flow pressure, and temperature were recorded.
The test section (Figure 6a) consisted of a large-diameter pipe (D = 50 mm), a small-diameter pipe (D = 25 mm), and two annular samples (super 13Cr stainless steel), which were fixed on the reducing wall by four screws (Figure 6b). The reducing wall was divided into two parts according to Particle Layer I, which was in front of Particle Layer J, represented particles that initially impacted the wall and were not subject to particle collisions. By contrast, the particles in Layer J may suffer from particle-particle collision interference near the wall. In each particle layer, particles were randomly positioned, thus satisfying the multiple possibilities of particles impacting the wall. Therefore, the number of particles mentioned in Table 3 is an average value and may change with different particle-producing positions. Meanwhile, a particle-moving calculation at each velocity was repeated 10 times to count the particle impact parameters on different radial walls because not all surfaces in the calculation process were impacted. Moreover, the inducing wall was considered to be under a no-slip condition, where u r | r = u x | r = 0 = 0 and

Boundary Conditions
The detailed parameters for fluid, particle, and geometry are summarized in Table 3. As shown Figure 4, the particles were generated on two layers (Particle Layers I and J) of virtual quadrilateral lanes, with the two planes radially distributed and spaced from each other by 0.05 mm. Each virtual lane was10 mm in length, 3 mm in width, perpendicular to the pipe axis and close to the pipe wall. ne of the purposes of this setting was to induce the maximum number of particles to impact the ducing wall because only particles near the wall of the upstream pipe can impact the wall. Another urpose was to distribute the particles as far as possible in the radial direction of the target wall.
Particle Layer I, which was in front of Particle Layer J, represented particles that initially pacted the wall and were not subject to particle collisions. By contrast, the particles in Layer J may ffer from particle-particle collision interference near the wall. In each particle layer, particles were ndomly positioned, thus satisfying the multiple possibilities of particles impacting the wall. herefore, the number of particles mentioned in Table 3 is an average value and may change with ifferent particle-producing positions. Meanwhile, a particle-moving calculation at each velocity was peated 10 times to count the particle impact parameters on different radial walls because not all rfaces in the calculation process were impacted. Moreover, the inducing wall was considered to be nder a no-slip condition, where ur|r = ux|r = 0 = 0 and әur/ ә x = 0 and a fully developed compressible flow was considered at the velocity inlet and exit sections.

Boundary Conditions
The detailed parameters for fluid, particle, and geometry are summarized in Table 3. As shown in Figure 4, the particles were generated on two layers (Particle Layers I and J) of virtual quadrilateral planes, with the two planes radially distributed and spaced from each other by 0.05 mm. Each virtual plane was10 mm in length, 3 mm in width, perpendicular to the pipe axis and close to the pipe wall. One of the purposes of this setting was to induce the maximum number of particles to impact the reducing wall because only particles near the wall of the upstream pipe can impact the wall. Another purpose was to distribute the particles as far as possible in the radial direction of the target wall.
Particle Layer I, which was in front of Particle Layer J, represented particles that initially impacted the wall and were not subject to particle collisions. By contrast, the particles in Layer J may suffer from particle-particle collision interference near the wall. In each particle layer, particles were randomly positioned, thus satisfying the multiple possibilities of particles impacting the wall. Therefore, the number of particles mentioned in Table 3 is an average value and may change with different particle-producing positions. Meanwhile, a particle-moving calculation at each velocity was repeated 10 times to count the particle impact parameters on different radial walls because not all surfaces in the calculation process were impacted. Moreover, the inducing wall was considered to be under a no-slip condition, where ur|r = ux|r = 0 = 0 and әur/ ә x = 0 and a fully developed incompressible flow was considered at the velocity inlet and exit sections.
x = 0 and a fully developed incompressible flow was considered at the velocity inlet and exit sections.

Full-Scale Experiment
A full-scale pipe flow experiment was used to study particle clusters during erosion and to verify the simulation results. As shown in Figure 5, the experimental loop mainly comprised a contraction section, a screw pump (with a flow range of 5-60 m 3 /h), an electric heating agitator, a temperature and pressure sensor, a magnetic flowmeter (8712HR, Emerson, Rosemount. Co., Boulder, CO, USA), four gate values, a control cabinet, a computer, and connecting pipes. After the test fluid containing the particles was mixed, stirring and electric heating units were turned on until the temperature reached the required value and then the pump was turned on. The gate valve in the test pipe was opened when the flow reached stability and the related data, including a particle motion image, flow pressure, and temperature were recorded.  Figure  6c). The mass concentration of the particles in the suspension was 50 kg/m 3 . The flow velocities of the experimental liquid were set to 1.5, 2.5, and 3.5 m/s, respectively. The exposed surface was sealed with epoxy resin and ground using SiC emery paper grade 1200 prior to installation. A MEMRECAM SP-614 high-speed camera (NAC. Co., Ltd., Tokyo, Japan) was used to observe the particle clusters with a spatial resolution of 0.01 μm and the sample surface profiles were verified using an H1200WIDE confocal scanning laser microscope (Lasertec. Co., Ltd., Tokyo, Japan) with a scanning rate of 120 fps.  Figure 8. The results show notably different erosion depths on the inner edge surface and the erosion difference gradually decreased along the radial direction. The maximum differences in erosion depth were 18.37, 11.96, and 8.43 μm, at the flow velocities of 1.5, 2.5, and 3.5 m/s, respectively. By ignoring the effects of metal corrosion, fatigue damage, and surface hardening, we identified that the critical cause of the erosion difference was particle collision interference, which diminished the impact energy of a fraction of the particles on the wall. As shown by the measurement results, the effect of particle interference was magnified by more particle impacts on the inner edge of the sample, whereas extremely small effects were observed in areas that were rarely impacted by particles (for a detailed description of the division of the erosion area on the reducing wall, we refer the reader to our previous work [45]).  The test section (Figure 6a) consisted of a large-diameter pipe (D = 50 mm), a small-diameter pipe (D = 25 mm), and two annular samples (super 13Cr stainless steel), which were fixed on the reducing wall by four screws (Figure 6b). The reducing wall was divided into two parts according to the particle impact frequency. The inner edge of the sample surface withstood multiple overlapping impacts with the formation of craters, platelets, and extruding lips. The outer circumference surface mainly suffered from an independent particle impingement. The test fluid was added with 0.4 wt.% hydroxypropyl guar gum as a thickener and 0.3 wt.% inorganic boron as a cross-linking agent. The viscosities of the liquid before and after crosslinking were 26 and 375 mPa. Round ceramsite, with a density of 1850 kg/m 3 and an average diameter of 0.65 mm, was used as an erosion abrasive (Figure 6c). The mass concentration of the particles in the suspension was 50 kg/m 3 . The flow velocities of the experimental liquid were set to 1.5, 2.5, and 3.5 m/s, respectively. The exposed surface was sealed with epoxy resin and ground using SiC emery paper grade 1200 prior to installation. A MEMRECAM SP-614 high-speed camera (NAC. Co., Ltd., Tokyo, Japan) was used to observe the particle clusters with a spatial resolution of 0.01 µm and the sample surface profiles were verified using an H1200WIDE confocal scanning laser microscope (Lasertec. Co., Ltd., Tokyo, Japan) with a scanning rate of 120 fps.      Figure 8. The results show notably different erosion depths on the inner edge surface and the erosion difference gradually decreased along the radial direction. The maximum differences in erosion depth were 18.37, 11.96, and 8.43 µm, at the flow velocities of 1.5, 2.5, and 3.5 m/s, respectively. By ignoring the effects of metal corrosion, fatigue damage, and surface hardening, we identified that the critical cause of the erosion difference was particle collision interference, which diminished the impact energy of a fraction of the particles on the wall. As shown by the measurement results, the effect of particle interference was magnified by more particle impacts on the inner edge of the sample, whereas extremely small effects were observed in areas that were rarely impacted by particles (for a detailed description of the division of the erosion area on the reducing wall, we refer the reader to our previous work [45]).    Figure 7b, most of the particles formed many particle clusters and gathered together across the upstream section. When the particle cluster impacted the wall, the particle-wall contact process may involve three steps: Initial contact between the front particle and the wall, rebound of the front particle, and collision interference between particles. The initial contact occurred between particle i and the metal surface (Figure 9a). If particles were too close, particle j in the cluster would collide with rebounding particle i (Figure 9b). Afterward, numerous subsequent particles may collide with each other and even form an unstable accumulation layer (Figure 9c).     Figure 7b, most of the particles formed many particle clusters and gathered together across the upstream section. When the particle cluster impacted the wall, the particle-wall contact process may involve three steps: Initial contact between the front particle and the wall, rebound of the front particle, and collision interference between particles. The initial contact occurred between particle i and the metal surface (Figure 9a). If particles were too close, particle j in the cluster would collide with rebounding particle i (Figure 9b). Afterward, numerous subsequent particles may collide with each other and even form an unstable accumulation layer (Figure 9c).  Figure 7b, most of the particles formed many particle clusters and gathered together across the upstream section. When the particle cluster impacted the wall, the particle-wall contact process may involve three steps: Initial contact between the front particle and the wall, rebound of the front particle, and collision interference between particles. The initial contact occurred between particle i and the metal surface (Figure 9a). If particles were too close, particle j in the cluster would collide with rebounding particle i (Figure 9b). Afterward, numerous subsequent particles may collide with each other and even form an unstable accumulation layer (Figure 9c). . Schematic diagrams of the particle-wall impact and particle-particle collision in a cluster, (a) distribution of the particles in pipe flow; (b) initial contact between the front particle and the wall; (c) the rebound of the front particle; (d) mutual interference of the particles.

As shown in
In the current study, the impact parameters between Particle Layer I and the reducing wall were used to obtain the distribution of the erosion crater depths, the particle-wall contact time, and the particle velocity attenuation. Meanwhile, the particles in Particle Layer J were used to determine the particle interference and calculate particle cluster erosion.

Particle-Wall Impingement
Target wall deformation included the extrusion and recovery processes. The time spanning all Figure 9. Schematic diagrams of the particle-wall impact and particle-particle collision in a cluster, (a) distribution of the particles in pipe flow; (b) initial contact between the front particle and the wall; (c) the rebound of the front particle; (d) mutual interference of the particles.
In the current study, the impact parameters between Particle Layer I and the reducing wall were used to obtain the distribution of the erosion crater depths, the particle-wall contact time, and the particle velocity attenuation. Meanwhile, the particles in Particle Layer J were used to determine the particle interference and calculate particle cluster erosion.

Particle-Wall Impingement
Target wall deformation included the extrusion and recovery processes. The time spanning all processes was dominated by the particle impact velocity and angle. Figure 10a,b show the normal and tangential velocity components before and after the particles impacted the wall, respectively. The particle tangential velocities were markedly higher than that of the normal velocities because of the steering effect of the particles [32,33], that is, the small particle impact angle increased the particle tangential velocity. According to the relationship between the velocity decay ratio and impact angle, the normal velocity loss rate changes from 60 to 90% along the radial surface (Figure 10a). Most of the impact energy was dissipated in the plastic work and elastic wave energies and the loss rate reached 90%, as confirmed by Hutchings [46,47]. Meanwhile, the tangential velocity loss rates, which were almost half of the normal velocity loss rate, approached 30% on the inner surface and 40% on the outer circumference ( Figure 10b).   The depth of an erosion crater can be approximated using Equation (31) when the particle impact velocity and angle were obtained. However, a significant empirical coefficient P n , which indicates the average pressure between a particle and the metal surface, needs to be obtained experimentally. The geometric dimensions of the 20 similar craters (Figure 11a) on the outer circumference of a sample surface were measured to determine the coefficient value. The obtained P n was equal to 2.1 × 10 6 Pa for super 13Cr stainless steel. Figure 11b shows the calculated erosion depth along the radial surface. The depth of an independent erosion crater gradually increased from the inner edge to the outer circumference, with maximum and minimum erosion depths of 17.8 and 2.9 µm, respectively. Hence, the calculated contact time using Equation (32) stabilized in the microsecond time scale and was notably lower than the particle release time (a difference of 1000 times) so that t 1 was approximately equal to t r . Figure 12 shows the release time for Particle I along the radial direction. The minimum release time for a particle that impacts the inner surface was 1.8 × 10 −4 s and the maximum release time was 3.3 × 10 −3 s. If a particle collision occurred after particle i rebounded from the wall, then a lower release time indicates that a shorter interparticle distance is required to allow two adjacent particles to collide.

Particle Cluster Erosion
According to the values of particle velocity, particle release time, and the force acting on the particle, critical interparticle distances were calculated along the radial direction and are shown in Figure 13. The distances increased gradually from the inner edge to the outer circumference, reflecting that more particles may collide with each other near the outer circumference owing to a Figure 11. Depths for an erosion crater obtained by experimental measurement and numerical calculation, (a) SEM micrograph of the single impact crater on the outer circumference surface; (b) erosion depth caused by independent particle impacts and cluster impacts versus experimental results.
(a) (b) Figure 10. Velocity components of a particle in Particle Layer I along the radial surface. (a) Normal impact velocities; (b) Tangential impact velocities.

Particle Cluster Erosion
According to the values of particle velocity, particle release time, and the force acting on the particle, critical interparticle distances were calculated along the radial direction and are shown in Figure 13. The distances increased gradually from the inner edge to the outer circumference, reflecting that more particles may collide with each other near the outer circumference owing to a

Particle Cluster Erosion
According to the values of particle velocity, particle release time, and the force acting on the particle, critical interparticle distances were calculated along the radial direction and are shown in Figure 13. The distances increased gradually from the inner edge to the outer circumference, reflecting that more particles may collide with each other near the outer circumference owing to a loose distance limit. By repeatedly calculating the distance between the two-layer particles when Particle I made contact with the wall, the clustered particle percentages of total impact particles were calculated and are exhibited in Figure 14. Calculations and statistical results showed that only 10% of the total impacted particles on the inner edge and 35% on the outer circumference were affected by collision interference before impacting the wall. Therefore, the particles impacting on the wall as a cluster were less than half of the total number of particles. Meanwhile, particles that were closer to the inner edge showed more independent particle impacts, thereby demonstrating the diminished effect of a high velocity gradient on clustered particle collision. calculated and are exhibited in Figure 14. Calculations and statistical results showed that only 10% of the total impacted particles on the inner edge and 35% on the outer circumference were affected by collision interference before impacting the wall. Therefore, the particles impacting on the wall as a cluster were less than half of the total number of particles. Meanwhile, particles that were closer to the inner edge showed more independent particle impacts, thereby demonstrating the diminished effect of a high velocity gradient on clustered particle collision.   Figure 15 shows the velocity decay coefficient for each particle in Layer J along the radial surface. The coefficient values approached 0.9 on the inner edge and 0.55-0.7 on the outer circumference, indicating that the kinetic energy consumption of Particle J on the outer circumference exceeded that on the inner edge. When the inlet flow velocity was altered, the velocity gradient and the angle of velocity component near the reducing wall partly changed and then the critical interparticle distance reached a new range. Figure 16 shows that the critical interparticle distance decreased with increasing flow velocity, reflecting that particles must be closer to each other for a collision to occur at a high flow velocity. In addition to its influence on the critical interparticle distance, an increasing flow velocity changed the decay coefficient of the particle velocity, as shown in Figure 17. Owing to the growth of the opposite collision velocity between two particles with an increasing flow velocity, the decay coefficient increased on average by 11% and approached 1 with an increased flow velocity by 1 m/s. Therefore, the effect of a cluster impact on the erosion calculation at a high flow velocity was less than that at a low flow velocity on the reducing wall.
The calculated results of the erosion depth by particle cluster and by independent particle impact in the radial direction are shown in Figure 11b for comparison. Compared with the calculated results of the total impacted particles on the inner edge and 35% on the outer circumference were affected by collision interference before impacting the wall. Therefore, the particles impacting on the wall as a cluster were less than half of the total number of particles. Meanwhile, particles that were closer to the inner edge showed more independent particle impacts, thereby demonstrating the diminished effect of a high velocity gradient on clustered particle collision.   Figure 15 shows the velocity decay coefficient for each particle in Layer J along the radial surface. The coefficient values approached 0.9 on the inner edge and 0.55-0.7 on the outer circumference, indicating that the kinetic energy consumption of Particle J on the outer circumference exceeded that on the inner edge. When the inlet flow velocity was altered, the velocity gradient and the angle of velocity component near the reducing wall partly changed and then the critical interparticle distance reached a new range. Figure 16 shows that the critical interparticle distance decreased with increasing flow velocity, reflecting that particles must be closer to each other for a collision to occur at a high flow velocity. In addition to its influence on the critical interparticle distance, an increasing flow velocity changed the decay coefficient of the particle velocity, as shown in Figure 17. Owing to the growth of the opposite collision velocity between two particles with an increasing flow velocity, the decay coefficient increased on average by 11% and approached 1 with an increased flow velocity by 1 m/s. Therefore, the effect of a cluster impact on the erosion calculation at a high flow velocity was less than that at a low flow velocity on the reducing wall.
The calculated results of the erosion depth by particle cluster and by independent particle impact in the radial direction are shown in Figure 11b for comparison. Compared with the calculated results  Figure 15 shows the velocity decay coefficient for each particle in Layer J along the radial surface. The coefficient values approached 0.9 on the inner edge and 0.55-0.7 on the outer circumference, indicating that the kinetic energy consumption of Particle J on the outer circumference exceeded that on the inner edge. When the inlet flow velocity was altered, the velocity gradient and the angle of velocity component near the reducing wall partly changed and then the critical interparticle distance reached a new range. Figure 16 shows that the critical interparticle distance decreased with increasing flow velocity, reflecting that particles must be closer to each other for a collision to occur at a high flow velocity. In addition to its influence on the critical interparticle distance, an increasing flow velocity changed the decay coefficient of the particle velocity, as shown in Figure 17. Owing to the growth of the opposite collision velocity between two particles with an increasing flow velocity, the decay coefficient increased on average by 11% and approached 1 with an increased flow velocity by 1 m/s. Therefore, the effect of a cluster impact on the erosion calculation at a high flow velocity was less than that at a low flow velocity on the reducing wall.
The calculated results of the erosion depth by particle cluster and by independent particle impact in the radial direction are shown in Figure 11b for comparison. Compared with the calculated results of erosion depth caused by an independent particle impact, the erosion depths caused by particle cluster were also determined in the radial direction. The value of the cluster impact erosion depth approximates that of the independent particle impact on the inner edge but shows a significant difference on the outer circumference. With a change in the fluid flow velocity, as shown in Figure 18, the main erosion position was on the inner surface and the erosion depth markedly increased on the inner edge with an increasing flow velocity. More particles impacted the wall at a high flow velocity because the particles possessed a higher inertial force to hit the wall head-on. Moreover, the particle impact areas also increased with the increased flow velocity.
approximates that of the independent particle impact on the inner edge but shows a significant difference on the outer circumference. With a change in the fluid flow velocity, as shown in Figure  18, the main erosion position was on the inner surface and the erosion depth markedly increased on the inner edge with an increasing flow velocity. More particles impacted the wall at a high flow velocity because the particles possessed a higher inertial force to hit the wall head-on. Moreover, the particle impact areas also increased with the increased flow velocity.   difference on the outer circumference. With a change in the fluid flow velocity, as shown in Figure  18, the main erosion position was on the inner surface and the erosion depth markedly increased on the inner edge with an increasing flow velocity. More particles impacted the wall at a high flow velocity because the particles possessed a higher inertial force to hit the wall head-on. Moreover, the particle impact areas also increased with the increased flow velocity.   18, the main erosion position was on the inner surface and the erosion depth markedly increased on the inner edge with an increasing flow velocity. More particles impacted the wall at a high flow velocity because the particles possessed a higher inertial force to hit the wall head-on. Moreover, the particle impact areas also increased with the increased flow velocity.

Discussion
The above results exhibit the effect of particle collision prior to impacting a wall on erosion and the following chapters discuss the changes in particle impact velocity under different particle collision forms. Based on the trajectory and collision form of two particles, the clustered particle

Discussion
The above results exhibit the effect of particle collision prior to impacting a wall on erosion and the following chapters discuss the changes in particle impact velocity under different particle collision forms. Based on the trajectory and collision form of two particles, the clustered particle erosion can be divided into three patterns, including IIE, interferential particle erosion (IPE), and stacked particle erosion (SPE). When the trajectories of the two particles are changed, the different particle-wall impact forms as shown in Figure 19. . Schematic diagrams of particle-wall and particle-particle collisions between two particles. Figure 20. The second particle impact velocity on the wall for interferential particle erosion (IPE) and stacked particle erosion (SPE). The second impact particle is particle j in IPE as well as is particle i in SPE.

Conclusions
The erosion characteristics of a reducing wall in a sudden contraction were analyzed by CFD-DEM coupling methods. A novel identification approach of particle cluster erosion was applied in the particle-wall impact calculation and the results provide the distribution of particle impact Figure 19. Schematic diagrams of particle-wall and particle-particle collisions between two particles. Independent Impact Erosion (IIE): When particle j moves on the outer side of particle i (Figure 19a), the particle i initially impacts the wall and then rebounds back to the other side. Meanwhile, particle j follows particle i and does not collide with the latter. The impact between the two particles and the wall is independent of each other.
Interferential Particle Erosion (IPE): When the trajectories of two particles are extremely close, as shown in Figure 19b, particle i impacts the wall and then collides with particle j. Particle j may impact the wall again after the collision, but the new impact energy of particle j should be attenuated during the particle-particle collision.
Stacked Particle Erosion (SPE): When the motion path of particle j is inside particle i (Figure 19c), particle j may be subjected to a head-on collision after particle i rebounds off the reducing wall. Therefore, the second impact between a particle and the wall may be attributed to particle i instead of particle j. This result causes the actual velocity in the second wall impingement to decrease twice.
IIE has been investigated in numerous aspects and is thus no longer discussed. As for IPE and SPE, different collision positions during momentum transfer between two particles result in different scattering angles and velocity decay rates for particle j. According to the relationship between the scattering angle and the new impact angle, that is, α 1 = α − ϕ (Figure 2a), the particle impact velocities in the second particle-wall impingement at different new impact angles are displayed in Figure 20. Given that the particle during the second impingement suffers twice the momentum loss in SPE and only once in IPE, the particle impact velocity in IPE is almost three times the velocity in SPE on the inner edge and almost 10 times on the outer circumference. Therefore, the effect of SPE on decreasing erosion extends far beyond the IPE and both erosion types exert more influence than IIE. Figure 19. Schematic diagrams of particle-wall and particle-particle collisions between two particles. Figure 20. The second particle impact velocity on the wall for interferential particle erosion (IPE) and stacked particle erosion (SPE). The second impact particle is particle j in IPE as well as is particle i in SPE.

Conclusions
The erosion characteristics of a reducing wall in a sudden contraction were analyzed by CFD-DEM coupling methods. A novel identification approach of particle cluster erosion was applied in the particle-wall impact calculation and the results provide the distribution of particle impact velocity, impact angle, erosion crater depth, critical interparticle distance, and velocity decay coefficient along the radial surface. The influence of cluster collision on the erosion rate was verified by a full-scale experiment. Meanwhile, the simulation results of particle cluster erosion were compared with the experimental results with different decay velocities. The conclusions derived from the present study can be summarized as follows: Figure 20. The second particle impact velocity on the wall for interferential particle erosion (IPE) and stacked particle erosion (SPE). The second impact particle is particle j in IPE as well as is particle i in SPE.

Conclusions
The erosion characteristics of a reducing wall in a sudden contraction were analyzed by CFD-DEM coupling methods. A novel identification approach of particle cluster erosion was applied in the particle-wall impact calculation and the results provide the distribution of particle impact velocity, impact angle, erosion crater depth, critical interparticle distance, and velocity decay coefficient along the radial surface. The influence of cluster collision on the erosion rate was verified by a full-scale experiment. Meanwhile, the simulation results of particle cluster erosion were compared with the experimental results with different decay velocities. The conclusions derived from the present study can be summarized as follows: (1) Formation of particle clusters results in interparticle collisions when certain particles bounce back from the wall and the succeeding particle erosion rate decreases with the change in particle impact velocity. (2) During particle impact and rebound, the contact time of the target surface is less than that of the particle rebound time before the particle collides with another particle. The decay rate of a normal impact velocity changes from 60 to 90% along the radial surface and the decay rate of the tangential velocity changes in the range of 30 to 40%. (3) A smaller critical interparticle distance, which indicates that particles must be closer to each other for a collision to exist, results in a lower clustered particle percentage of the total impact particles. Therefore, the probability and frequency of a particle cluster impact on the outer circumference are larger than those on the inner edge. Meanwhile, by updating the particle impact velocity and angle, the relative error between the calculation results and the experimental results was reduced by almost 11% compared with that of the complete independent impact setting. (4) The second impact velocity of the particle in IPE is almost three times the velocity in SPE on the inner edge and almost 10 times on the outer circumference because the particle velocity decays once in IPE and twice in SPE. The growth of flow velocity not only decreases the critical interparticle distance but also the decay coefficient in particle collision, indicating that the effect of a cluster impact on erosion can be weakened by an increasing flow velocity.
Author Contributions: J.C., Z.L. and N.Z. designed the experiments and supervised experimental work; J.C. and Z.W. wrote the program; Y.D. and N.Z. verified and changed the manuscript; Y.D. provided experimental instruments.
Acknowledgments: This work was supported by the National Natural Science Foundation of China (grant no. 51674199), and it was made possible by Study on Working Behavior of Subtending Pipe Column Supporting Tools (grant no. 290017039).

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
A p particle cross-sectional area, m 2 t e contact time, s C d drag coefficient, dimensionless t r release time, s C n normal damping coefficient, kg·s −1 t 1 rebound time, s C t tangential damping coefficient, kg·s −1 t 2 moving time for particle j, s D inside diameter of circular pipe, m ∆t D DEM time step, s d p diameter of particle, m ∆t Ra Rayleigh time step, s e n restitution coefficient in the normal direction, dimensionless T r rolling friction torque, N·m e t restitution coefficient in the tangential direction, dimensionless T t torque generated by tangential forces, N·m E elastic modulus, Pa u fluid flow velocity, m E R total erosion rate, kg/kg ∆v the relative velocities of the centers of the spheres before and after a collision, m·s −1 E RA erosion rate caused by independent particle impact, kg/kg v n particle normal velocity, m·s −1 E RI erosion rate caused by grouped particle impact, kg/kg v t particle tangential velocity, m·s −1 E* equivalent elastic modulus, Pa v ref the relative velocity of the contact points, m·s −1 ∆E n normal energy loss of a particle caused by particle-particle collision, J V p volume of particle, m 3 ∆E t tangential energy loss of a particle caused by particle-particle collision, J V l volume of liquid cell, m 3 F b buoyancy, N ∆x space step, m F d drag force, N Y Young's modulus, Pa F g gravitational force, N a empirical constants in angle functions, dimensionless F l,i the fluid force acting on a particle, N b F n inter-particle contact force in the normal direction, N w F p,i the particle force acting on the liquid, N x F ij contact force between particles, N y F t inter-particle contact force in the tangential direction, N z F s particle shape coefficient, dimensionless Greek symbols F n d inter-particle damping force in the normal direction, N ω i unit angular velocity vector of the object, rad·s −1 F t d inter-particle damping force in the tangential direction, N λ decay coefficient in the collision, dimensionless G shear modulus, Pa κ particle velocity loss rate in process of particle-wall impact, % G* equivalent shear modulus, Pa ν poisson ratio, dimensionless H maximum depth of indentation, m δ tangential displacement, m ∆h depth of the extruding lips, m µ dynamic viscosity, Pa·s I moment of inertia, m 4 β The angle between the particle velocity and the centerline of the two particles, • J particle momentum, kg·m·s −1 δ k identity tensor K dimensionless coefficient, dimensionless δ n normal overlap, m L distance from the inner edge of the sample, m δ t tangential overlap, m L p critical inter-particle distance of particle, m ρ l liquid density, kg·m −3 L p,e the moving distance of the particle j along the centerline, m µ s friction coefficient, dimensionless L p,r the critical inter-particle distance in the release process, m µ r rolling friction coefficient, dimensionless L y particle spacing in the Y direction, m ϕ scattering angle for particle j, • L z particle spacing in the Z direction, m ψ scattering angle for particle i, • m p single particle mass, kg α particle impact angle, • m* equivalent mass, kg α 1 new particle impact angle, • N total particle number, dimensionless α l fluid volume fraction, % N p the number of sample points contained within the mesh cell of the particle, dimensionless θ critical angle of particle impact, • N 1 the number of particles in the particle group, dimensionless θ* critical angle between particles, • N t total number of sample points of the particle, dimensionless ε strain caused by the normal stress of indentation, m P n function of properties of the material, Pa Superscripts r normal unit vector, dimensionless n flow behavior index r p radius of the erosion crater, m m impact velocity power law coefficient R radius of particle, m (1) pre-collision R i maximum radius of deformation, m (2) post-collision R* equivalent radius, m Subscripts R l inside radius of a circular pipe, m i particle i S n normal stiffness, N·m −1 j particle j S t tangential stiffness, N·m −1 z normal direction t tangential unit vector, dimensionless y tangential direction