Impact of inter-particle collision on elbow erosion based on DSMC-CFD method

In this work, computational fluid dynamics (CFD) is used to study elbow erosion due to a gas–solid two-phase flow. In particular, the direct simulation Monte Carlo (DSMC) method is used to study the impact of inter-particle collision on the erosion behavior. The two-way coupled Euler–Lagrange method is used to solve the gas–solid flow, and the DSMC method is used to consider the collision behavior between particles. The effects of key factors, such as the particle concentration distribution and inter-particle collision, on the erosion ratio are evaluated and discussed. The effectiveness of the method is verified from experimental data. The results show that the inter-particle collision significantly influences the particle movement path and erosion ratio. When the inter-particle collision is considered, the maximum erosion position is offset. The erosion model proposed by Oka et al., who used the DSMC method, agrees best with the experimental data, and the average percentage error decreases from 39.2 to 27.4%.

(0) Velocity before collision, m/s e Coefficient of restitution (dimensionless) n Normal unit vector directed from particle i to particle j (dimensionless) t Unit vector in the tangential direction (dimensionless) G (0) Relative velocity of the particles before collision, m/s G ct (0) Tangential component of (

Introduction
Pipeline erosion is one of the main problems in pneumatic conveying systems, particularly in the field of oil and gas. Erosion may lead to equipment failure or even plant failure, limiting the efficiency and safety of production (Zhang et al. 2007). It is important to develop an effective method to predict elbow erosion distribution and build a high-precision erosion prediction model for reducing maintenance time and saving resources (Pereira et al. 2014). The fundamental cause of erosion is the interaction between particles and wall, resulting in the removal of wall materials. The velocity and impact angle of the particles play a major role in the erosion phenomenon (Karimi et al. 2017). Many erosion prediction models have been proposed to reveal the interaction mechanism between particles and wall. Finnie (1960) put forward the theory of micro-cutting of plastic materials for the first time. Bitter (1963a, b) classified erosion wear of surface materials into deformation wear and cutting wear and proposed an erosion model. Neilson and Gilchrist (1968), Hutchings (1981), Zhang et al. (2007), and Oka et al. (2005) proposed empirical and semi-empirical erosion prediction models based on experimental results. Computational fluid dynamics (CFD)-based erosion prediction technology has been effectively applied to predict erosion under specific conditions. Zhang et al. (2012) studied the effects of slurry velocity, bending direction, and angle of an elbow on erosion by conducting a CFD-based numerical simulation. The maximum erosion position shifted downstream with the increase in the flow velocity. Duarte et al. (2015) numerically studied the effect of particle mass concentration on elbow erosion and found that the erosion rate decreases with the increase in the particle mass concentration. The reason given was the impact of inter-particle collision; however, the mechanism of inter-particle collision was not discussed. Peng and Cao (2016) studied the erosion of an elbow due to a liquid-solid two-phase flow using the two-way coupled Euler-Lagrange method. Three mechanisms of particle collision were proposed to explain the effect of Stokes number on erosion morphology. Zhang et al. (2018a) studied the influence of mesh structure on the erosion prediction of small particles and proposed a CFD mesh generation strategy to better match the experimental results. Zeng et al. (2018) used the CFD-discrete element method (DEM) coupling method to study the effect of particle shape on the erosion behavior in a dense particle flow. The particle sphericity was found to affect the contribution of key parameters, such as the impact frequency, velocity, and elbow angle, to the erosion rate. Li et al. (2017) used the nonlinear finite element method (FEM) to study deformation behavior and different failure modes of unpressurized and pressurized pipes under bending conditions. Laín and Sommerfeld (2019) calculated erosion depth for two-and four-way coupling and for mono-sized spherical glass beads as well as a size distribution of particles. When particle mass loading is increased, bend erosion is reduced due to modifications of particle impact velocity and angle, although wall collision frequency grows. Azimian and Bart (2016) proposed a model based on the Euler-Euler approach; the effects of input flow velocity, solid concentration, and the particle size on erosion depth and erosion spatial distribution are investigated. Zhu and Li (2018) used CFD-DPM Eulerian-Lagrangian approach to evaluate the anti-erosion effect. Particles are more uniformly distributed in a lowspeed flow, rebound effect attenuates the direct collisions against the same spot. Farokhipour et al. (2019) used the CFD-DEM method to study the effects of mass loading and one-, two-and four-way coupling regimes on the model prediction on erosion rate. The simulation results indicated that, as the particle mass loading increases, the effectiveness of plugged tee compared to the standard elbow increases, and influence of particles interaction with flow and particle collision cannot be neglected. Zhao et al. (2020) used dense grid by one-way coupling of solid phase to simulate the single-phase impinging jet due to its dilute distribution, and the simulation results agreed well with experiments. These studies have promoted the development of CFD-based erosion prediction technology; however, there is no research on the impact of inter-particle collision on the erosion behavior. The main reason is that the Euler-Lagrange particle trajectory model is used to track the trajectories of a large number of real particles in the flow field and to find collision pairs from the trajectories in a deterministic manner. If the number of particles increases, huge computational resources and computational time are required, making it impossible to perform the calculation based on individual particles.
In the field of rarefied gas dynamics, the direct simulation Monte Carlo (DSMC) method is widely used to simulate particle motion in gas-particle two-phase flows, as it can effectively deal with particle motion and collision. Xiong et al. (2017) used the grand canonical Monte Carlo simulation method to study the methane adsorption behaviors in slit-like chlorite nanopores, and the influences of the pore sizes, temperatures, water, and compositions on methane adsorption on chlorite were discussed. In the DSMC method, real particles are replaced with a small number of sampled particles, and the inter-particle collisions are determined randomly, rather than deterministically as in the case of the discrete phase model (DPM) or DEM. This helps save computational power.
The purpose of this work was to study the impact of interparticle collision on the erosion of a 90° elbow. We compared the simulation and experimental results to verify the effectiveness of the DSMC-CFD erosion prediction method.

Numerical model
In this study, the Eulerian-Lagrangian method is used to solve the gas phase as a continuous phase based on the Navier-Stokes equation, and the solid particles are treated as discrete phases and are solved using Newton's second law. The numerical simulation involves four models: a continuous phase model, a discrete phase model, an inter-particle collision model, and an erosion model.

Continuous phase model
The Navier-Stokes equation is used to model the continuous phase. The general equations of continuity and momentum are as follows: Here ρ is the continuous phase density, kg/m 3 ; ⃗ u is the instantaneous velocity vector, m/s; p is the pressure, Pa; is the stress tensor; ⃗ g is the acceleration due to gravity, m/ s 2 ; and � ⃗ S D is the additional source term for the interaction of particles with the continuous phase.
The standard k-ε turbulence model is used for turbulence modeling in the following form: Here t is the turbulence viscosity: G k is the generation term for the turbulent kinetic energy due to the velocity gradient; u i is the velocity component in the i direction, m/s; k is the turbulent kinetic energy, J; is the turbulent dissipation rate, W/m 3 ; x i and x j are the space coordinates, m, i ≠ j ; k (= 1) and ε (= 1.3) are the turbulent Prandtl numbers for k and , respectively; S k and S ε are source terms; C 1 = 1.44, C 2 = 1.92, and C μ = 0.09.

Discrete phase model
In the DPM, the translational motion of the particles in the gas phase conforms to Newton's second law of motion, and the trajectory of the solid particles can be obtained by integrating the motion equation of the particles in Lagrangian coordinates. The governing equation of particle motion is as follows: Here m i and � ⃗ v i , respectively, represent the mass and velocity of the particle i , and ∑ F i is the sum of all the forces acting on the particles, including the drag force ⃗ F d,i , gravity ⃗ F g,i , pressure gradient force ⃗ F p,i , and virtual mass force ⃗ F vm,i . The drag force is the main force acting on the solid particles and is defined as follows: Here ⃗ v i is the velocity vector of the solid particles, m/s; d i is the diameter of the particles, m; i is the density of the solid particles, kg/m 3 ; and Re s is the Reynolds number of the solid particles: C d is the drag coefficient: Here a 1 , a 2 , and a 3 are constants for the smooth spherical particles given by Zhang et al. (2012).
The fluid affects the dispersed phase through drag and turbulence, while the particles affect the fluid because of the reduction in average momentum and turbulence, resulting in momentum exchange between the gas and dispersed phases. In this study, the two-way coupling method is used to solve the interaction between the particles and the fluid.
Re 2 s

Inter-particle collision model
When the rotation of the particles is neglected, the movement of the sampled particles can be divided into translation motion and collision process. The translation motion of the sampled particles obeys the DPM model, and the collision process follows the inter-particle collision model (Zhang et al. 2018b).

Particle collision dynamics model
Assuming that the particles are rigid spheres, the collision process between the particles is instantaneous and is controlled by an impulse equation, which is not affected by other forces. If a collision occurs between particles i and j, the velocity after collision can be calculated using the impulse equation, as follows: Here J is the impulse vector applied to particle i during collision, and (0) is the velocity before collision. Through mathematical deduction, the velocity of the two particles after collision is as follows: Here e is the coefficient of restitution, is the coefficient of friction, is the normal unit vector directed from particle i to particle j in collision, is the unit vector in the tangential direction, G (0) is the relative velocity of the particles before collision, and G (0) ct is the tangential component of G (0) .

Collision probability
Assuming that the particles are randomly distributed in a certain spatial range, particle i may collide with any particle in this range. The collision probability depends on the relative velocity, size, and local concentration of the particles, as shown in Fig. 1.
The collision frequency of particle i can be calculated using the following formula: Here D i and D j are the diameters of particles i and j, respectively, is the mode of relative velocity of particles i and j, V 0 is the volume of the search area, the search radius is Δt p , N is the actual number of particles in the search range, and n is the actual number of particles represented by the sampled particles. Pawar et al. (2014) showed that an unbiased estimation of the average collision frequency can be made when at least eight samples are taken from each search area volume. Therefore, we take N n ≥ 8. The average free path of the particles is the average distance between particles. According to the collision frequency, the average free path i of particle i with velocity v i can be calculated as follows: The particle time step Δt p is calculated as follows: Here Δt g is the time step of the continuous phase.
Therefore, the collision probability of particles i and j during the particle time step can be obtained from the following formula: The modified Nanbu method is used to determine whether collision occurs (Zhang et al. 2018b). First, a random number R with a uniform distribution ranging from 0 to 1 is extracted from a random number generator, and the serial number of the collision particle pair j is generated using Eq. (16).
Here int[R × N] is an integer operation of R and N.
The collision probability of the particles is calculated using Eq. (15). The collision occurrence criterion is expressed as follows: When Eq. (17) is satisfied, the particles i and j will collide in time step Δt p . The velocity after collision is calculated using Eqs. (9)-(11). Figure 2 shows a flowchart of the sample particle calculation process in the DSMC method.

Erosion prediction model
Many models and methods have been proposed for predicting the erosion rate. In the current erosion prediction models, the erosion rate is modeled as follows: Here ER is the ratio of the mass loss of the target material to the mass of the eroded particles, F s is the shape coefficient of the particles, BH is the hardness of the target material, K and n 1 are empirical constants based on the properties of the materials, u p is the particle velocity, n 2 is the velocity correlation index based on experience, and F( ) is the collision angle function. The erosion models have a strong empirical correlation, and the parameters of the. different models vary significantly. The erosion prediction models considered in this study include those proposed by Oka et al. (2005), DNV model (Haugen et al. 1995), Zhang et al. (2007), Neilson andGilchrist (1968) model, andVieira et al. (2016). The specific expression formulae and detailed parameters of each model can be found in their papers.
The particle-wall collision rebound model is used to calculate the change in velocity after particle-wall collision. The random particle-wall collision rebound model proposed by Grant and Tabakoff (1975) is used. The expressions are as follows: Here e n and e t are the restitution coefficients in the normal and tangential directions, respectively.
(20) e t = 0.988 − 1.66 + 2.11 2 − 0.67 3 The calculation process for particle i The calculation for the next particle i + 1 Calculate mean free path λ i , using Eq.(13) Calculate particle time step ∆t p ,using Eq. (14) Generate a random number R j = int[R•N] + 1 Calculate collision probability P ij , using Eq. (15) Particles i and j collide

Geometric model and mesh generation
The subject of this study is a 90° elbow with a diameter of 76.2 mm, a curvature radius of 1.5D (114.3 mm), an upstream pipe length of 15D (1143 mm), and a downstream pipe length of 10D (762 mm), as shown in Fig. 3. The near-wall meshing strategy in the CFD model is very important to the stability, and accuracy of the simulation results in this article refers to Zhang et al. (2018b) research. The geometry of the runner is simple, hexahedron elements are used for the structural mesh, and using enhanced wall function near-wall modeling approach. The first layer thickness is set to 1.1 mm, the maximum Y + value is 80, the boundary layer expansion coefficient is 1.2, and the number of boundary layers is 10. Grid independence analyzed the number of grids from 100,000 to 15 million. The results show that the erosion ratio changes only slightly with the increase in the number of meshes, and no significant change occurs when the number of meshes is more than 3.66 million. Therefore, the number of grids in this study is set to 3.66 million.

Simulation conditions
Solid particles and air are injected uniformly through the inlet of the elbow at the same speed. An inlet velocity boundary condition is applied to the inlet, and a free outflow boundary condition is applied to the outlet. The pipe wall is the wall boundary, the wall roughness constant is set to 0.5, and the turbulence intensity is set to 5%. Previous studies have shown that the erosion ratio has no effect on the number of particles when the number of particles is more than 20,000. To ensure a sufficient number of particles, 50,000 particles were injected into the entrance of the simulation model.
The calculation conditions are consistent with the experimental conditions set by Vieira et al. (2016). The gas velocity is in the range of 11-27 m/s, the mass flow rate of the solid particles is in the range of 103-452 kg/day, the mass loading (the ratio of solid mass flow rate to gas mass flow rate) is in the range of 0.015-0.059, and the particle diameters are 150 and 300 µm. The air density is 1.2 kg/m 3 , the wall density is 7990 kg/m 3 , and the solid particle density is 2650 kg/m 3 . Table 1 lists 11 sets of experimental conditions obtained from the study by Vieira et al. Figure 4 shows the fluid velocity distribution at the XY section of the elbow under condition 3. There is an obvious velocity gradient between the center of the elbow and the wall, and the fluid velocity at the center of the elbow is generally greater than that near the wall. When considering  inter-particle collision, the velocity gradient near the wall is closer, and the boundary layer is thinner. This is because the particles deviate from the main streamline due to collisions and consequently impact the outer wall and transfer their kinetic energy to the fluid near the wall. Because of the coupling mechanism involved in the exchange of kinetic energy between the particles and the fluid, the velocity of the fluid at the inlet extrados and outlet intrados of the 90° elbow is higher. When the fluid enters the elbow, the gas velocity near the extrados of the elbow decreases, whereas that near the intrados increases, resulting in a significant change in the pressure. The extrados pressure of the 90° elbow increases, and the intrados pressure decreases. The influence of particle collision on the fluid pressure distribution can be neglected. Figure 5 shows the distribution of the fluid pressure field.

Particle concentration distribution
The mass loading ratio is defined as the ratio of the maximum local particle mass concentration to the average inlet particle mass concentration in the entire flow process. Figure 6 shows the mass loading ratio during the erosion simulation process. The closer the color to the blue area in the figure, the closer the local maximum concentration to the uniformly distributed mass loading. The closer the color to the red area, the greater the local maximum concentration and particle concentration in this area. The transparent areas indicate that very few particles pass this region. The particles are distributed evenly in the straight section of the entrance, and the local concentration of the particles varies dramatically after entering the 90° elbow. Because of inertia, the particles in the straight section accumulate at the extrados of the elbow when passing through the elbow, and at the same time, a particle free zone is formed at the intrados side of the elbow because of the shadow effect. The particles rebounding at the extrados contact the wall partially, forming a "V-shaped" high-concentration particle  area under the action of secondary flow. The maximum local particle concentration occurs at the intersection point of the V-shaped wings and is 19.86 times the average particle mass concentration at the entrance. Figure 7 shows the particle mass loading ratio and the velocity vector at the cross sections perpendicular to the axis of the elbow. The effect of secondary flow on the local concentration of the particles is more obvious, as shown in Fig. 7. Initially, at the cross section (a), the directions of the velocity vectors are the same, and the particle concentration distribution is more uniform. At cross sections (b) and (c), the directions of the velocity vectors change, and the fluid near the extrados is separated from the main stream, and a large number of particles accumulate at the extrados. At cross section (d), the directions of the velocity vectors are consistent again, and the particles diffuse to the periphery of the aggregation zone, as observed from the particle concentration distribution. Figure 8 shows the particle trajectories with and without considering inter-particle collision under condition 3. As shown in Fig. 8, when the collision between the particles is not considered, the particle trajectories are concentrated in a V-shaped region, similar to that shown in Fig. 6, and a wider particle free zone is generated because of the particle motion mechanism. When considering the collision between particles, only a small number of particles collide in the straight section of the entrance, with most of the inter-particle collisions occurring in the high concentration area of the elbow. Some particles collide outside of the V-shaped region, resulting in a narrower particle zone. The mass loading at the inlet under condition 3 is 0.059, and the volume fraction is 2.67 × 10 −5 . Moreover, even the particles in the dilute particle flow may have 1-2 times or more collisions at the elbow.  Figures 9 and 10 show the relationship between the particle mass flow rate, particle diameter, and particle collision location distribution. The inter-particle collisions are mainly distributed at the 90° elbow. With the increase in the particle mass flow rate, the proportion of the number of collisions in the straight pipe section to the number of collisions in the entire pipe increases; however, it is still far less than that at the elbow. The inflection point, shown in Fig. 9, is due to the change in the flow velocity, and the number of interparticle collisions decreases with the increase in the inlet velocity. Figure 10 shows the effect of particle diameter on the inter-particle collisions. The total number of collisions decreases with the increase in the particle diameter. This is because with the increase in the particle size, the number of particles decreases, and the collision pairing of the particles becomes more difficult. Figure 11 shows the erosion simulation results with and without considering particle collision under condition  3. Because of the special particle trajectory at the elbow extrados, the "V-shaped" erosion morphology is formed on the outer wall of the elbow. The maximum erosion ratio considering inter-particle collisions is less than that without considering them. Notably, the erosion ratio decreases at the maximum erosion position when particle collision is not considered and increases around and outside the intersection of the two wings of the V-shaped erosion scar. Therefore, the collisions between the particles affect the erosion ratio distribution of gas-solid flows. Figure 12 shows a comparison of the distribution of the erosion ratio, particle trajectory, particle concentration, and collision location. The particle trajectory is the most dense, the local concentration is the highest, and the collision frequency is the highest where the erosion ratio is maximum. Therefore, the impact of particle collision should be considered when predicting the maximum erosion location.

Effect of inter-particle collision on erosion behavior
By comparing the three key parameters affecting the erosion, we can have a better understanding of the change in the erosion mechanism after particle collision. Figure 13 shows a comparison of the variations in the number of particle hits to the wall, the average impact angle, and the average impact velocity. The trajectory of particle hits to the wall forms a V-shaped high-frequency impact area. The number of particle hits to the wall is highest at the end of the "V-shaped" area, and it decreases after considering interparticle collision. The average impact angle and velocity do not exhibit any obvious change in this area, because most of the collisions are first collisions between the upstream particles and the inner wall of the elbow. In the outer flanks of the V-shaped region, the number of hits, the average impact angle, and the average impact velocity of the particles on the wall are significantly increased. The difference is due to the change in the particle trajectory caused by particle collision. Figure 14a and b shows the difference between the CFD and DSMC-CFD results in predicting the maximum and total erosion ratios, respectively. In Fig. 14a axis represents the percentage bias of the maximum erosion ratio predicted using the CFD and maximum erosion ratio predicted using the DSMC-CFD, and the transverse axis represents the total number of collisions between the particles. When the total number of collisions is less, the difference is more. With the increase in the total number of collisions between the particles, the percentage error of the maximum erosion ratio predicted using the CFD and  Fig. 12 Relationship between erosion ratio, particle trajectory, particle concentration, and collision location  . 13 Effects of inter-particle collisions on particle parameters: a Number of particle hits, b Average impact velocity, c Average impact angle DSMC-CFD is approximately in the range of 5-15%, indicating that the maximum erosion ratio predicted using the DSMC-CFD is 5-15% lower than that predicted using the CFD. Figure 14b shows the total erosion ratio obtained using the two methods. The total erosion ratios calculated using the CFD and DSMC-CFD are the same, implying that the total metal loss due to the particles calculated using the two methods is similar. Figure 15 shows the distribution of the erosion ratio along the center line of the elbow, predicted using the models proposed by Oka et al., DNV model, Zhang et al., Neilson and Gilchrist, and Vieira et   collision, the maximum erosion location is approximately 45° on the elbow axis. After considering the inter-particle collision, the maximum erosion position shifts downstream by 5°, i.e., to approximately 50° on the elbow axis. The erosion ratio at 45° decreases significantly, from 15.2 to 17.3%. The maximum erosion ratio considering particle collision is less than that without considering it. The numerical simulations are carried out considering the 11 working conditions, five erosion models, and two cases (with and without inter-particle collision). Figure 16 shows a comparison of the erosion prediction results with the experimental values. The transverse axis represents the maximum erosion ratio measured by experiment, and the longitudinal axis represents the maximum erosion ratio predicted by numerical simulation. Figure 16 shows that the average error between the predicted maximum erosion ratio and the experimental values for the different erosion models varies when particle collision is considered in the modeling. The models proposed by Oka et al., Zhang et al. and Vieira et al. overestimate the actual erosion ratio. The error is less when using the DSMC method wherein the particle collision is considered; the percentage error is reduced by 11.8, 14.1, and 27.4%, respectively. The models proposed by DNV, Neilson and Gilchrist underestimate the actual erosion ratio; the error when using the DSMC method is increased by 3.6 and 5.1%, respectively. The model proposed by Oka et al. with the DSMC method agrees best with the experimental data, and the average percentage error between the model and experimental data decreases from 39.2 to 27.4%.

Conclusions
The DSMC method was used to study the impact of interparticle collision on elbow erosion due to a gas-solid twophase flow. To verify the models and methods analyzed in this study, the simulation results are compared with the experimental results. The following conclusions can be drawn from this work: (1) Under the action of inertia and secondary flow, the particles form a high concentration region, and the maximum local concentration is 19.86 times the average concentration. The particles in the dilute particle flow may collide 1-2 times or more in the high concentration region. The particle mass flow rate, inlet velocity, and particle size affect the frequency of particle collisions.
(2) When inter-particle collision is considered in the DSMC method, the lateral erosion ratio of the V-shape erosion trace increases, and the predicted maximum erosion ratio decreases by 5-15%. This is mainly because the inter-particle collision results in a change in the location and number of hits.
(3) Considering inter-particle collision, the maximum erosion position is offset by 5° downstream, i.e., to approximately 50° on the elbow axis. The erosion ratio at 45° decreases from 15.2 to 17.3%. (4) When inter-particle collision is considered through the DSMC method, the errors are reduced for the erosion model that overestimates the actual erosion ratio, whereas the errors are increased for the model that underestimates the actual erosion ratio. The model proposed by Oka et al., wherein the DSMC method is used, agrees best with the experimental data, and the average error decreases from 39.2 to 27.4%.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.