Numerical Simulation of Erosion Wear for Continuous Elbows in Different Directions

Abstract

The purpose of the present study is to simulate the continuous bend erosion process in different directions, using the dense discrete particle model (DDPM). The influence of the length of the straight pipe in the middle of the continuous bend is investigated. The Rosin–Rammler method is introduced to define the diameter distribution of erosion particles, which is theoretically closer to the actual engineering erosion situation. The numerical model is based on the Euler–Lagrange method, in which the continuous phase and the particle phase are established on a fixed Euler grid. The Lagrange model is used to track the particles, and the interaction between particles is simulated by particle flow mechanics theory. The velocity field distribution, pressure variation, and turbulent kinetic energy of gas–solid two-phase flow, composed of natural gas and gravel in the pipeline, are studied. The simulation results, using the one-way coupled DPM and the four-way coupled DDPM, are compared and analyzed. The results show that the DDPM has good accuracy in predicting the distribution of the continuous bend erosion processes in different directions. The erosion rates of particles with an average distribution size of 50 μm are significantly increased (8.32 times), compared with that of 10 μm, at the same gas transmission rate. It is also indicated that it is important to consider the impact between particles and the coupling between fluid and particles in the erosion simulation of the continuous elbow when using the CFD method.

1. Introduction

Natural gas is clean and efficient energy, with methane as the main component, and its combustion products are mostly water and carbon dioxide, which plays an important role in today’s world energy field. However, the process of natural gas exploitation and transportation still faces many challenges. In natural gas exploitation, although some filtering measures are taken to filter solid particles in produced gas, they still cannot completely prevent the passage of some particles with very small diameters [1]. Erosion caused by these small particles is one of the main problems faced by oil and gas production systems [2]. Over time, the damage caused by this erosion could lead to catastrophic energy security and environmental problems. Therefore, the study and analysis of erosion damage are very important for oil and gas production, the energy field, and other industries.

Due to the importance of erosion, many researchers have carried out various experimental and simulation studies to determine the influencing factors of erosion severity. Hong et al. [3] analyzed the influence of gas flow rate, solid mass flow rate, particle diameter, and other factors on the maximum erosion rate. It was found that the pipe diameter and elbow curvature radius were negatively correlated with the maximum erosion rate. They proposed a new dimensionless correlation to predict the maximum corrosion rate of the elbow. Mamoo et al. [4] studied the effects of inlet flow rate, temperature, particle density, and diameter on tube wall erosion, through numerical simulation by using the bidirectional coupling Eulerian–Lagrangian method and the Oka erosion and Grant and Tabakoff particle-wall rebound model methods. It is found that the increase in fluid velocity and temperature, particle mass flow rate, and particle density will increase the erosion rate. Akrami et al. [5] revealed the evolution process of pipe size in the pipeline process and the erosion mechanism of BEP in vertical layered sand, by measuring and giving the depth, width, and the ratio of depth to width at the tip of pipe in the critical erosion stage. They developed the idea that the three-dimensional study of pipe is more meaningful than the two-dimensional study.

Agrawal et al. [6] proposed an algorithm that coupled the erosion rate calculation with a moving deformation grid, which is a solid wall geometry and computational grid, based on dynamic deformation of the local erosion rate, within a certain time interval. The results show that computational fluid dynamics (CFD) are a powerful tool for erosion prediction and can provide detailed erosion patterns with complex geometric shapes. High local erosion rates can lead to significant changes in pipe geometry. Singh et al. [7] used Fluent to study the erosion wear of 90-degree bends. The k-epsilon turbulence model was used to track solid particles, and the erosion rate of the continuous flow field was calculated by the erosion model. The effects of velocity, particle size, and concentration on solid–liquid interaction were studied. It was found that erosion wear increases exponentially with velocity, particle size, and concentration, and the maximum value and location of erosion wear are more serious in the elbows than in the straight pipes. Nam et al. [8] used computational fluid dynamics (CFD) to simulate the erosion rate of the pipe bend, rationalized the geometry and grid structure, and used a discrete phase model (DPM) and corrosion model to calculate the failure cycle of pipe, by using the limit state function of the erosion rate. Xu et al. [9] established a finite element model of erosion for a gas–solid two-phase flow pipeline, according to the structural characteristics and working conditions of a gas–solid two-phase flow pipeline in a gas transmission station. They studied the erosion characteristics of gas–solid two-phase flow pipelines with different pressures, solid content, throttle valve opening, and pipe diameter. The prediction equation of the maximum erosion rate is presented. Pei et al. [10] used the CFD method to investigate the relationship between flow field, particle trajectory, maximum erosion zone, and influencing factors in curves. It was found that the Stokes number of particles moving in the bend has no decisive influence on the position of the maximum erosion zone, and the erosion zone is directly related to the particle diameter. The law that increasing the radius of curvature will change the flow field in the pipeline and then change the position of the maximum erosion zone was discovered.

Zhang et al. [11] studied the effects of the turbulence model and particle rebound model and they found that the grid configuration has a significant impact on fine-grained sand erosion prediction, due to the limitations of existing CFD codes, and proposed a diffusion and migration mesh strategy that captures fine sand, which provides a significant guide for the meshing of this study.

Some researchers have improved the particle phase motion model of computational fluid dynamics. Cloete et al. [12] used a refined analytical gas–solid flow simulation using an improved Lagrange method, known as the dense discrete phase model (DDPM). Compared with the standard Euler fluid method, the DDPM takes the forces between particles into account, greatly improves the mesh independent behavior, and easily contains particle size distributions. Pouraria et al. [13] studied the effects of particle loads on the erosion rate and flow field in the two-phase flow of gas sand, transported by pipe bending, based on the dense discrete particle model. The results obtained by the DDPM and one-way coupled DPM were compared with experimental data, and it was proved that the DDPM has better accuracy in predicting erosion distribution.

However, the previous research mainly focused L-pipes, U-pipes, and Tee-pipes [14,15,16,17]. The study on the erosion of continuous bend in different directions is still lacking, especially the effect of the length of straight pipe in the middle, on the erosion of bend. In the transportation of natural gas, some terrain factors or user demands are often encountered, and the pipeline has to be bent continuously to meet the demand. Therefore, relevant research is urgently needed to provide references for the energy field.

In this paper, the erosion mechanism and law of continuous elbows with different lengths, with the transportation of natural gas-containing sand, are studied by using the CFD method, based on the gas–solid two-phase flow particle erosion theory, according to the actual working conditions of natural gas transmission pipeline. The Rosin–Rammler distribution function [18] is used to set the grain diameter of sand, so that the diameter of erosion particles can be randomly generated, according to certain rules, which are closer to reality. We use the DDPM to consider a collision between particles and the influence of the coupling of the fluid and particles., etc., and to more accurately reveal the straight pipe length between successive bend to the maximum erosion rate and erosion morphology, particle trajectories and continue the influence law of flow field. The influences of the length of straight tube between continuous bends on the maximum erosion rate, erosion morphology, particle trajectory and flow field in continuous media are revealed more accurately by the method. The present paper provides a theoretical basis for gas pipeline improvement, selection, installation, and use specification revision.

2. Numerical Model

2.1. Mathematical Model for Continuous Phase

The dense discrete particle model builds the continuous phase on a fixed Euler grid and uses the Lagrange model to track the particles, based on the Euler–Lagrange method. In this model, the void ratio and collision of particles are considered, and the calculation of collision is modeled instead of the calculation of collision process by the soft sphere model. The collision between particles is described by particle flow mechanics theory, and the conservative equation of continuous phase is as follows:

(1)$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g\right)+\nabla ·\left({\alpha }_g{\rho }_g\overrightarrow{v_g}\right)=0$

(2)$\frac{\partial }{\partial t}\left({\alpha }_g{\rho }_g\overrightarrow{v_g}\right)+\nabla ·\left({\alpha }_g{\rho }_g\overrightarrow{v_g}\overrightarrow{v_g}\right)=\nabla ·\left[{\alpha }_g{\mu }_g\left(\nabla \overrightarrow{v_g}+\nabla {\overrightarrow{v_g}}^T\right)\right]+{\alpha }_g{\rho }_g\overrightarrow{\mathrm{g}}-{\alpha }_g\nabla p+\overrightarrow{F_{ex}}$

The above two equations are a continuity equation and momentum equation, respectively. The subscript $g$ denotes the natural gas fluid, $\alpha$ is the volume fraction and $\rho$ is the density of the continuous phase, kg/m³; $v$ represents the fluid velocity, m/s; $p$ indicates the fluid pressure, Pa. $F_{ex}$ denotes the exchange forces between two phases, which can be calculated based on the employed forces in DPM method.

The realizable k-ε model adds a new transmission equation to the dissipation rate, which can better complete the calculation and have higher calculation accuracy with adequate grid quality. The turbulent kinetic energy $k$ and the dissipation rate $S$ are calculated using the following differential equation [19]:

(3)$\rho \frac{Dk}{Dt}=\frac{\partial }{\partial X_n}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial K}{\partial X_n}\right]+G_k+G_{b }-\rho \epsilon -Y_M$

(4)$\rho \frac{D\epsilon }{Dt}=\frac{\partial }{\partial X_n}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_t}\right)\frac{\partial \epsilon }{\partial X_n}\right]+\rho \epsilon S{C}_{1 }+G_b{C}_{1\epsilon }{C}_{3\epsilon }\frac{\epsilon }{K}-\rho C_2\frac{{\epsilon }^2}{K+\sqrt{\gamma \epsilon }}$

Here, $\gamma$ is the viscosity coefficient of molecular motion and ${\mu }_t$ is the turbulent viscosity coefficient. $C_{1\epsilon }$, $C_2$,$C_{3\epsilon }$, ${\sigma }_k$, and ${\sigma }_t$ all take the default values. $G_k$ is the turbulent kinetic energy effect caused by the average velocity gradient. $G_b$ is the turbulent kinetic energy effect caused by buoyancy. $Y_M$ represents the influence of turbulent pulsating expansion on the total dissipation rate. The formula for $C_1$ is shown as follows. $\eta$ is the apparent viscosity of the fluid.

(5)$C_1=max\left[0.43, \frac{\eta }{\eta +5}\right]$

2.2. Dense Discrete Particle Motion Model

The trajectory of a particle is calculated by the equilibrium equation of the force acting on the particle. The conservation equation of particle phase is not solved, which is contrary to the standard multi-fluid method [13]. The Lagrange method calculates the trajectory on each parcel by integrating the force balance on each parcel. The volume fraction and the velocity field of the particle phase are obtained directly from the Lagrange method. The trajectory of a particle is calculated by the equilibrium equation of the force acting on the particle. The motion equation of a particle in the Cartesian coordinate system is given by the following equation:

(6)$m_p\frac{d\overrightarrow{u_p}}{dt}=m_p\frac{{\rho }_g}{{\rho }_p}\overrightarrow{v_g}\nabla \overrightarrow{v_g}+\frac{3}{4}{C}_d{m}_p\frac{\left|\overrightarrow{u_p}-\overrightarrow{v_g}\right|}{d_p}\left(\overrightarrow{v_g}-\overrightarrow{u_p}\right)+m_p\overrightarrow{\mathrm{g}}\frac{\left|{\rho }_p-{\rho }_g\right|}{{\rho }_p}+{\overrightarrow{F}}_{KTGF}$

In the above equation, the first term on the right side of the equation is pressure gradient force, the second term is the drag force, and the third term is buoyancy. The first three items are common in the DPM method. The last item is added to the DDPM. $m_p$ is the mass flow rate of the sand particles, kg/s; $C_d$ is the drag coefficient; $\overrightarrow{u_p}$ is the particle velocity, m/s. Subscript $g$ represents natural gas fluid and subscript $p$ represents solid particles. Solid particles in natural gas are equivalent to spherical particles with random diameter, following Rosin–Rammler size distribution. The KTGF model is used to predict the stresses and forces on particles, caused by collisions between particles and spherical particle translation. The KTGF forces are modeled from the solid phase stress tensor. The normal viscous stress, shear viscous stress, and normal pressure are added together to form the solid phase stress tensor.

(7)${\overrightarrow{F}}_{KTGF}=-\frac{m_p}{{\rho }_p}\nabla \stackrel{=}{{\tau }_s}$

The calculation method of $\stackrel{=}{{\tau }_s}$ in the above equation is as follows [13]. $p_p$ denotes the sand particle pressure, which can be calculated by the coefficient of restitution for particle collisions ${\vartheta }_p$, the granular temperature $T_p$, and the radial distribution function $g_0$.

(8)$\stackrel{=}{{\tau }_s}={\alpha }_p{\mu }_p\left(\nabla \overrightarrow{{\mu }_p}+\nabla {\overrightarrow{{\mu }_p}}^T\right)+{\alpha }_p\left({\lambda }_p-\frac{2}{3}{\mu }_p\right)\nabla \overrightarrow{v_p}\stackrel{=}{I}-p_p\stackrel{=}{I}$

(9)$p_p=2{\rho }_p{\alpha }_p{}^2{g}_0{T}_p+2{\rho }_p{\vartheta }_p{\alpha }_p{}^2{g}_0{T}_p+{\alpha }_p{\rho }_p{T}_p$

Here, $\stackrel{=}{I}$ is the unit stress tensor; ${\mu }_p$ represents the shear forces resulting from the momentum exchange of particles due to translation and collision; ${\lambda }_p$ is the particle bulk viscosity of the particle. The radial distribution function $g_0$ can be calculated as the following equation [20], where ${\alpha }_{p,max}$ is the packing limit. In the present study, the coefficient of restitution for particle collisions ${\vartheta }_p$ is set to 0.90 and the ${\alpha }_{p,max}$ is set to 0.63.

(10)$g_0={\left[1-{\left(\frac{{\alpha }_p}{{\alpha }_{p,max}}\right)}^{\frac{1}{3}}\right]}^{-1}$

The Rosin–Rammler distribution functions have long been used to describe the size distribution of particles of various types and sizes. This function is especially suitable for representing particles with different diameters in natural gas pipelines in the actual environment. Compared with the particle model, whose particle diameter is assumed to be constant, the method presented in this paper is closer to the real situation of engineering application [21]. The governing equation of Rosin–Rammler particle size distribution can be expressed as follows.

(11)$W_p=100e^{{\left[-\left(\frac{D}{De}\right)\right]}^n}$

Here, $W_p$ is the retained weight fraction, %; $n$ is the Rosin–Rammler skewness parameter; $D$ is the particle size, m; $De$ is the equivalent mean diameter, m. In the present study, the range of equivalent mean particle diameter $De$ based on Rosin–Rammler distribution is 10~50 μm, and the volume fraction is less than the packing limit [13], so the influence of friction viscosity is ignored. The calculation method [22] of ${\mu }_p$ in the above equation is as follows.

(12)${\mu }_{p,c}=\frac{4}{5}{\alpha }_p{}^2{\rho }_{p}De{g}_0\left(1+{\vartheta }_p\right){\left(\frac{T_p}{\pi }\right)}^{\frac{1}{2}}$

(13)${\mu }_{p,k }=\frac{{\alpha }_{p}De{\rho }_p\sqrt{\pi T_p}}{6\left(3-{\vartheta }_p\right)} \left[1+\frac{2g_0}{5}\left(3{\vartheta }_{p }-1\right)\left(1+{\vartheta }_p\right){\alpha }_p\right]$

(14)${\mu }_p={\mu }_{p,c}+ {\mu }_{p,k}$

The calculation method of the particle bulk viscosity ${\lambda }_p$ in the above equation is as follows [23].

(15)${\lambda }_p=\frac{5}{3}{\mu }_{p,c}$

The rotation of a particle has a significant effect on its trajectory in a fluid. In this study, particle flow through two consecutive bends will inevitably produce rotation. In this case, if the particle rotation is ignored in the numerical calculation, the particle flow trajectory and erosion results will be significantly different from the actual results. Here, we solve an ordinary differential equation for angular momentum to take the rotation of the particle into account.

(16)$I_p\frac{d{\overrightarrow{\omega }}_p}{dt}=\frac{{\rho }_g}{64}D{e}^5{C}_w\left|\overrightarrow{\Omega }\right|·\overrightarrow{\Omega }=\overrightarrow{T_o}$

Here, $I_p$ denotes the moment of inertia of the particle, kg·m²; ${\overrightarrow{\omega }}_p$ is the angular velocity of the particle, rad/s; ${\rho }_g$ indicates the density of natural gas fluid, kg/m³; $C_w$ is the rotational resistance coefficient; $\overrightarrow{T_o}$ is the torque applied to the sand particle, N·m; $\overrightarrow{\Omega }$ is the relative particle–fluid angular velocity, which can be calculated as follows:

(17)$\overrightarrow{\Omega }=\frac{1}{2}\nabla \times v_g-{\overrightarrow{\omega }}_p$

Particles in this paper are equivalent to spherical particles that obey Rosin–Rammler particle size function distribution, so the moment of inertia $I_p$ can be calculated as follows:

(18)$I_p=\frac{\pi }{60}{\rho }_{p}D{e}^5$

As described, the four-way coupled DDPM considers the theory of particle flow me-chanics on the basis of the DPM coupling of buoyancy, pressure gradient force, and drag force, adding coupling of collision between particles, translation, and particle and fluid coupling.

2.3. Erosion Model

The general expression of erosion prediction models is as follows [24]:

(19)$E_R={\displaystyle \sum }_{p=1}^{N_p}\frac{m_{p}C\left(d_p\right)f\left(\vartheta \right)v^{n\left(v\right)}}{A_f}$

Here, $E_R$ denotes the erosion ratio, which is defined as the rate of the mass loss rate at the wall to the particle mass flow rate per unit area per unit time, kg/(m²s); $A_f$ indicates the unit superficial area of the wall; $m_p$ is the mass flow rate of the sand particles, kg/s; $N_p$ is the number of the particles that impact the wall; $C\left(d_p\right)$ is the particle diameter function and the value is 1.8 × 10⁻⁹. $f\left(\vartheta \right)$ can be expressed as a piecewise linear function [25], as shown in Table 1, that can be obtained through normalizing the erosion data. $n\left(v\right)$ is the particle impact velocity function and the exponent value is set to 2.6 [26].

When particles impact the fitting wall, energy is lost and the reflection velocity is less than the incident velocity. According to previous research methods, the recovery coefficient (the ratio of spring back to incident velocity) is usually used to describe this effect, as follows [24]:

(20)$e_n=0.988-0.78\vartheta +0.19{\vartheta }^2-0.02{\vartheta }^3+0.027{\vartheta }^4$

(21)$e_{t }=1-0.78\vartheta +0.85{\vartheta }^2-0.21{\vartheta }^3+0.028{\vartheta }^4-0.022{\vartheta }^5$

Here, $e_n$ represents the normal restitution coefficients and $e_t$ represents the parallel restitution coefficients. Simultaneously, the boundary condition of the wall for the discrete phase is set as “reflect”. $\vartheta$ is the impact angle of the particles.

3. Physical Model and Method Validation

3.1. Geometric Model and Grid Independence Validation

The physical model is composed of five parts. The parts of the computing domain are Straight-1, Elbow-1, Straight-2, Elbow-2, and Straight-3, from the inlet to the outlet, as shown in Figure 1. The length of the straight pipe between the two elbows is an important variable, which increases from 1 m to 7 m, with a span of 1 m. A total of seven 3D physical models are established and the dimensions are listed in Table 2. The calculation domain is meshed with hexagonal structured units, and the number of the boundary layer is 10. Local mesh encryption is carried out at the elbows to enhance the accuracy of results. The y+ number is tested by calculation and the value basically satisfies 30 < y+ < 300 in all cases. The scalable wall functions are used to prevent the deterioration of numerical results and calculation errors of unbounded shear forces. The minimum value is 0.261, the maximum value is 1, and the mean value is 0.802, when using the Element as the evaluation criterion. As the Skewness is used as the criterion, the average value is 8.566 × 10⁻². Good mesh quality can further ensure the reliability of the CFD results.

The physical model of L = 4 m was selected to examine the variation of the maximum erosion rates of the three sections (Elbow-1, Elbow-2, and Straight-2) with the number of grids. The calculation of grid refers to the study of Zhang et al. [12], and the data results of grid independence verification are shown in Figure 2. The similarity is that we also use structured grids and divided boundary layers. The thickness of the first layer of the wall grids is strictly controlled to avoid an inaccurate prediction. The difference is that the size transition of radial mesh is gentler. Although the calculation amount is increased to a certain extent, the poor mesh quality caused by excessive aspect ratio is prevented and the convergence is improved. When the mesh unit number reaches 1,307,132, the influence of increasing the mesh density on the calculation results is negligible. Considering the accuracy and efficiency of the calculation, the method and size of this mesh are used to complete the mesh drawing of other physical models.

3.2. Boundary Conditions and Calculating Method

The boundary condition setting includes velocity inlet (12 m/s) and pressure outlet (gauge pressure 0). The particle size of erosion particle flow follows the Rosin–Rammler function and is divided into three groups. The particle size of the first group ranged from 8 μm to 12 μm, with an equivalent average particle size of 10 μm. The particle size of the second group ranged from 16 μm to 24 μm, with an equivalent average particle size of 20 μm. The particle size of the third group ranged from 40 μm to 60 μm, with an equivalent average particle size of 50 μm. The inlet and outlet are set to “escape” and the wall is set to “reflect”. The pressure–velocity double precision solver is used to solve the boundary conditions using the pseudo-transient method and the Coupled algorithm. The PRESTO! pressure solution method is selected to obtain a better convergence effect and maintain high accuracy. The second order upwind algorithm is used to solve the momentum and turbulent kinetic energy, and the convergence standard of the residual in the monitor is set as 10⁻⁶. The Warped-Face Gradient Correction (WFGC) method is used to solve the problem of gradient accuracy reduction, caused by the extremely high aspect ratio of the grid at boundary layer, or non-flat surface elements and highly deformed elements, whose centroid of formal elements is outside the control volume.

3.3. Model and Calculation Method Validation

The numerical model calculation method, described above, is used to conduct a simulation and verification calculation, with the same conditions as the experimental study of Mazumder et al. [27]. The medium is air and the inlet velocity is 15.24 m/s. The equivalent average diameter of particles following the Rosin–Rammler distribution function is set at 300 μm. The results are compared to verify the simulation modeling and method of erosion wear, as seen in Table 3. The model studied in the experiment is the inner wall erosion area of the S-shaped bend. The physical model studied in this paper is somewhat different from the experimental model mentioned above, so only the erosion area range at Elbow-1 is compared. The results show that the error between the proposed method and their experimental results is less than 10%. Compared with the DPM, the DDPM has smaller errors, and the simulation results are in good agreement with the experimental results.

The calculation model used in this paper was used to set the conditions, consistent with Pouraria [13], with a velocity of 11 m/s and a particle size of 150 μm, for the second data verification. The results show the proposed method in this paper and the results of the Pouraria method have a high degree of coincidence, as shown in Figure 3. In addition, Ou et al. [25] used a similar method to assist in the proof. It is verified that the calculation method used in this paper can reliably predict the erosion wear of elbows, and can be used in engineering calculations, through repeated parameter adjustment and repeated calculation.

4. Results and Discussion

All the residuals of cases are less than 10⁻⁶. The difference between the inlet mass flow rate and the outlet mass flow rate is less than 0.5%. Based on the above analysis, the convergence standard is judged to be reached. The results and discussion of the cases will be below.

4.1. Pressure Distribution and Variation along the Path

Some iso-surfaces at different positions are established in the computational domain to clearly observe the pressure field distribution and changes along with the flow, when passing through the two elbows, as shown in Figure 4. The pressure gradient of the fluid of Elbow-1 is less than that of Elbow-2. There is a narrow area of low pressure inside Elbow-1 that increases with the deflection angle. The area of low pressure inside Elbow-2 is much larger than Elbow-1. The shape of the low-pressure region changes from b–0° to b–90°. As L decreases from 7 m to 1 m, the pressure gradient at Elbow-2 increases. When L is long enough, the fluid will get more adequate flow development in Straight-2, after passing through Elbow-1. The distance makes the influence of the flow field brought by the first bend gradually decrease. The gas–solid two-phase flow will be affected by centrifugal force when passing through the elbows, resulting in local low pressure. As the length of L decreases, the two-phase flow through Elbow-2 is still greatly affected by the centrifugal deflection of Elbow-1, which further enlarges the influences on the pressure gradient at Elbow-2.

Figure 5 shows the data of the maximum facet total differential pressure of the iso-surfaces and the outlet, which can represent the magnitude of the pressure gradient. When the length of Straight-2 is 5 m, the maximum facet total differential pressure keeps the minimum in the range of a–0° to a–90° of Elbow-1 and b–0° to b–45° of Elbow-2. The pressure gradient of these surfaces is maximum at L = 1 m in the range b–0° to b–60° of Elbow-2. The maximum facet total differential pressure occurred at the a-90° position of Elbow-1, reaching 194.8 Pa, corresponding to L = 7 m.

Pressure loss is caused by collision, friction, and other factors in the process of fluid flow. The pressure change curve along the flow is shown in Figure 6. The pressure curves coincide almost exactly when passing flow through Elbow-1, indicating that the flow field at Elbow-1 is not affected by the length of Straight-2. As the value of L increases gradually, the total distance of the fluid in the computational domain will also increase, and the pressure drop loss will be higher, up to 65.07 Pa.

4.2. Velocity and Particle Trajectory

In gas–solid two-phase flow, the fluid movement has a direct impact on the elbow erosion rate. The velocity field distribution and the vector are shown in Figure 7. The fluid entering the elbow is decelerated by the impact and friction of the elbow wall. The fluid velocity at the inner side of Elbow-1 is high. By comparing the pressure contours in Figure 4, it can be seen that the cause of high velocity may be the acceleration in pressure difference, and the fluid in the high-pressure area flows to the low-pressure area with the action of pressure. The flow rate increases to 20.01 m/s in the process, which exceeds the inlet velocity at 12 m/s. The phenomenon is also present in Elbow-2 and there is a low-speed vortex in the center of the bend that is smaller than the main flow velocity. The flow entering the elbow rotates to the outer side of the upper wall at low speed and collides with the wall, causing erosion and reflecting to change the movement track. These particles join the surrounding particle stream and spiral down the outer wall together, as shown in the velocity vector diagram and the particle track diagram in Figure 8. The particle flow in the process will cause a severe impact on the elbow and its adjacent wall.

Particle trajectories subject to the Rosin–Rammler particle size distribution function are tracked, as shown in Figure 8, which shows 80,169 incoming particles of the particle flow, numbered and distinguished by different colors. It is divided into three columns, according to the average particle size. Particle trajectories vary from the length of Straight-2. The reason for the difference seems to be that the gas–solid two-phase flow passing through Straight-2 needs to overcome gravity and friction, so the particle energy decreases gradually. In addition, the centrifugal force through the elbows also affects the particle trajectory. Particle trajectories using the DPM are compared with those calculated by the DDPM, as shown in Figure 9. The interactions between particles in the DPM are ignored, so there are deviations from particle trajectories calculated by the model used in this paper.

4.3. Turbulent Kinetic Energy

Turbulent kinetic energy is half of the product of turbulence velocity fluctuation variance and fluid mass, which is a significant index to measure turbulence intensity. It determines the ability of a flow to remain turbulent or become turbulent and indicates the stability of the fluid flow. The value of the maximum turbulent kinetic energy of the gas–solid two-phase fluid increases gradually, as it passes through the two elbows. The slope of increase is greater at Elbow-2 than at Elbow-1, as shown in Figure 10a. The minimum turbulent kinetic energy of Elbow-1 is almost identical when the length of Straight-2 is different, as per Figure 10b. However, the values of the minimum turbulent kinetic energy in Elbow-2 vary and tend to increase gradually. It can be found that the turbulent kinetic energy increases gradually, and the flow stability decreases during the process of fluid passing through the elbows. The flow process in Straight-2 will gradually decrease after exiting from Elbow-1. The phenomenon indicates that the straight pipe flow contributes to improved flow stability. It should be noted that the flow cannot be fully developed as the length of the straight pipe is 1 m.

The numerical values of the maximum and minimum turbulent kinetic energy solved by the DPM are compared with the data results solved by the DDPM used in the present paper, and the results are shown in Figure 10c,d. It can be found that the calculation results of the maximum turbulent kinetic energy in the elbows are different, and the small deviations are from a–0° to a–60° of Elbow-1 and from b–0° to b–60° of Elbow-2. It should be noted that the simulation results of the DPM for a–90° and b–90° differ from the DDPM. The reason may be that the turbulent kinetic energy of the fluid at the positions is more affected by the interaction between particles, which is improper to ignore.

4.4. Accretion Rate and Erosion

When the natural gas fluid with particles is transported and flowing in the pipeline, the particles are affected by various factors and may be diffused in space or deposited on the wall. It is of great significance for the erosion investigation to study the deposition of fine particles on the wall surface, to find the law of the deposition position of fine particles on the wall surface, and to explore the influence of the different particle sizes and lengths of Straight-2 on the deposition of fine particles on the wall surface. Figure 11 shows the variation in the maximum accretion rate in different computational domains. The spline curves are interpolated by the Akima method [28]. When the average particle size is 10 μm, the gas–solid two-phase flow does not cause the maximum deposition when passing through Elbow-1, due to the small average particle diameter and a good following of natural gas fluid. The maximum deposition rate is located behind Elbow-1, appearing in Straight-2. When the length of Straight-2 is 3 m, the max DPM accretion rate curve shows two maximal values at the elbows.

The reason for the low max accretion rate at L = 1 m may be that the two-phase flow continuously passes through two elbows, within a very short distance, resulting in extremely high turbulence, which promotes the particles to be carried away by the gas flow at a very fast speed. The accretion capacity is reduced. As the length of Straight-2 increases, the max accretion rate remains low. The reason seems to be that the particles need to overcome gravity to function as they move vertically upwards. The potential energy increases and the kinetic energy decreases gradually, which weakens the impact of the wall deposition. When the length is 7 m, a maximum occurs at Straight-2. The reason may be that the velocity of particle flow decreases gradually, due to gravitational factors, when passing through a long vertical pipe, resulting in flow arrest, to a certain extent. The phenomenon results in secondary reflux of the particles that had passed Elbow-1, exacerbating sedimentation. When the average particle size is 20 μm, the max accretion rate curve generally shows two maximal values at the elbows, as shown in Figure 11b. It can be seen from Figure 8 that when the length of Straight-2 is 2 m, the particle trajectory of gas–solid two-phase flow in the continuous bend and intermediate fluid domain is smooth. The integrity of particles is stronger and the occurrence of the max accretion rate lags. When the average particle size is 50 μm, the max accretion rate increases by an order of magnitude compared to that at 20 μm, as shown in Figure 11c. The centrifugal force of the gas–solid two-phase flow through the bend has a significant effect on the trajectory of fine particles, which is the main reason for the max accretion rate on the outer wall surface. According to Ibrahim et al. [29], the maximum deposition rate is related to the flow velocity, residence time, secondary flow, electrostatic force, blocking rate, and other factors of gas–solid two-phase flow. The max accretion rate at Elbow-1 peaks at 0.00035412729 kg/(m²s), when Straight-2 is 5 m in length. It is speculated that the current calculation conditions are most favorable for the deposition of the elbow wall surface.

In addition, the max accretion rate calculated by the four-way coupled DDPM with the data calculated by the one-way coupled DPM are compared, as shown in Figure 12. It is found that the results for the max accretion rate calculated in the Straight-1 region are almost identical. Because the turbulence of the two-phase flow, flowing in a straight direction from the inlet is small, the force between particles has little influence on the motion trajectory. When the particle flow passes through Elbow-1, the flow path will change, becoming complicated. Ignoring the interaction forces can lead to bias. In particular, it can have a huge impact on Elbow-2 and the further calculations.

The relations between the mass flow rate of the particle flow and average particle diameter are shown in Table 4. The Erosion contours of different particle diameters in Elbow-1 and Straight-2 computational domains are shown in Figure 13. Particles enter from the bottom of the diagram and exit from the left. When the average particle size is 50 μm, the erosion area is concentrated on the lateral side of Elbow-1 and the two sides of Straight-2, resulting in a “Y” shape. When the average particle size is 10 μm, the erosions at the Elbow-1 site are concentrated in the connection between Elbow-1 and Straight-2. The two sides of Straight-2 are almost covered by the erosion zones. With the increase in average particle size, the area of maximum erosion is gradually closer to the elbow. This phenomenon also occurs in two computing domains, Elbow-2, and Straight-3, as shown in Figure 14. The reason for this phenomenon may be that when the particle size is small enough, the particles are more easily carried by the fluid. In other words, provided that the particle density is constant, the smaller the particle size is, the smaller its mass is, and the easier it is to be affected by the viscous force of the fluid. Under the influence of airflow, particles with smaller diameters erode the wall surface in a wider area, and the damaged area tends to be farther back.

When straight-2 is 5 m in length, the erosion area at the Elbow-2 location is small, as shown in Figure 14. The movement trajectory of particles is relatively smooth, as per Figure 8, which makes the movement of particle flow more stable and the collision area between particles and the wall surface more concentrated. The erosion area is smaller and the local maximum erosion rate is higher. As the average particle diameter is 50 μm and the Straight-2 length is 5 m, the maximum erosion rate reaches 1.1144613 × 10⁻¹¹ kg/(m²s) at Elbow-2, as shown in Figure 15. As the length of Straight-2 is 6 m, the erosion area is large, corresponding to the erosion contours with the average particle size of 10 μm and 20 μm. The obstruction occurred in the fluid domain in the straight pipe section. The particles overcome the gravity and work to reduce the velocity of the upper fluid. The lower fluid and the upper fluid are squeezed together, and the particles collide with each other to form a reflux zone and vortex zone, which aggravate the erosion and damage of the straight pipe area.

Most particles move in straight lines in the direction of flow. Only a few particles collide with other particles to change their trajectory and bounce off the wall. Therefore, the erosion rate of Straight-1 remains low, with different particle sizes, and the corrosion degree of the straight pipe wall is slight. When the average particle size is 10 μm, the maximum erosion rate corresponding to the L = 3 m is the highest.

As the average particle size is 20 μm, the overall average value of the maximum erosion rate, corresponding to the length of Straight-2, is 5 m, which is the highest. The gas–solid two-phase flow loses little kinetic energy when it works against gravity in the vertical pipe, as the length of Straight-2 is 1 m. The continuous passage of the two elbows makes the Straight-2 area highly subject to centrifugal force. The particles are thrown heavily towards the straight tube wall with the action of centrifugal force, resulting in serious erosion. When the average particle size is 20 μm, the max erosion rate reaches a peak.

When the average particle size is 50 μm, the erosion degree at the elbow sharply increases. The max erosion rate with the average particle size of 50 μm is about 8.32 times that of the average particle size of 10 μm. The simulation data for Straight-2 lengths less than 6 m demonstrate that the maximum erosion rates are all located at Elbow-2.

As the length of Straight-2 continues to grow, the kinetic energy of the particle flow to Elbow-2 is reduced due to gravity, resulting in lower erosion rates. The reason why the maximum erosion rate of the 7 m length is higher than that of the 6 m length is that the complex reflux in the latter straight pipe further dissipates the erosion energy.

The max erosion rate data and contours calculated by the one-way coupled DPM and the four-way coupled DDPM are compared. There are some differences in the erosion areas calculated by the two models, as shown in Figure 16. The DPM lacks many of the erosion details because it does not consider the interaction between particles, and also fails to capture the smooth and concentrated particle flow as L = 5 m, resulting in a large erosion area.

The difference between the maximum erosion rate calculated by the two models should not be overlooked. The values calculated by the DPM are close to those calculated by the DDPM, only in the average particle diameter of 10 μm and in Straight-1 and Elbow-1 domains, as shown in Figure 17, because the interaction between particles in the above two computational domains is still not obvious. When the gas–solid two-phase flow passes through Elbow-1, the particle trajectory changes dramatically, and the collision probability between particles increases. Ignoring the interactions between particles may lead to undesirable predictions. The results also confirm that it is necessary to use the four-way coupled DDPM when simulating the particle erosion of continuous elbows in different directions.

5. Conclusions

In this paper, the double-precision steady-state numerical simulation is carried out on seven continuous elbows in different directions, to investigate the flow field characteristics and particle erosion law of the straight pipe lengths and particle sizes as variables. The DDPM-KTGF method and the pseudo-transient calculation strategy are used. An integrated CFD method is established, including the discrete phase model, conforming to the Rosin–Rammler particle size function distribution law, the realizable k-ε turbulence model, and the erosion prediction model. The calculation results of the one-way coupled DPM and the four-way coupled DDPM are compared. The most important conclusions can be drawn as follows:

• (1). The erosion degree of the two elbows is different and the erosion rate of Elbow-2 is generally higher than that of Elbow-1. When the particle size ranges from 8 μm to 60 μm, the erosion degree increases with the increase in particle size, and the larger the particle mass, the stronger the destructive force to the tube wall.

• (2). In general, the degree of erosion of the elbow is larger than that of the straight pipe, but when the distance L between the two elbows is small, the velocity separation and secondary reflux phenomenon are obvious, and the local pressure gradient is large. Due to the high energy of the particles hitting the wall, the centrifugal force throws the particle stream towards the middle straight section, where the erosion degree is higher than that of the elbows.

• (3). As the distance L between the two elbows increases, the particle flow overcomes gravity to do its work and consumes kinetic energy, and the difference in the erosion between the two elbows weakens on the whole. As the length of L is 5 m, the particle flow trajectory is concentrated and smooth, resulting in a smaller, but more destructive, erosion area. As the length of L is 6 m, the phenomenon of flow blocking and an obvious whirlpool appear in the vertical pipe, which will aggravate the erosion of the central straight pipe. When continuous elbows are encountered in practical engineering, attention should be paid to the design of the wall thickness of the two elbows and the middle straight pipe, and a sufficient erosion margin should be considered. Erosion tests should be carried out regularly in this area to avoid potential safety hazards.

• (4). When the particle size is 50 μm, the erosion area is more concentrated, showing the shape of “Y” at the elbow. When the particle size is reduced to 10–20 um, the erosion area is more dispersed. The solid particles have a good concomitant with natural gas, and a large range of weak erosion phenomena will occur.

• (5). It is necessary to consider the interaction between particles in the erosion simulation of continuous bend pipes. It is recommended that readers use the four-way coupled DDPM to perform similar CFD simulation calculations, as computing resources are available.

Author Contributions

B.L. designed the study, conducted the numerical simulation, analyzed the data and wrote the paper, M.Z. and Q.W. reviewed and revised the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

The work is supported by the National Key Research and Development Program of China (2018YFE0108900).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

This paper is supported by HPC Platform, Xi’an Jiaotong University.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Schematic diagram of 3D geometric model and grid element.

Figure 2. Max erosion rate of Elbow-1, Elbow-2, and Straight-2 with different grid numbers.

Figure 3. Comparative data of the effect of particle volume fraction on erosion rate by using the two models.

Figure 4. Position diagram of iso-surfaces and pressure distribution contours.

Figure 5. Max Facet Total Differential Pressure.

Figure 6. Area-Weighted Average Total Pressure along the flow direction.

Figure 7. Position diagram of iso-surfaces and velocity distribution contours with vectors.

Figure 8. Particle trajectory diagram with average particle diameters and the length of Straight-2.

Figure 9. Particle trajectories calculated by the DDPM and the DPM.

Figure 10. The turbulent kinetic energy varies with the location of different iso-surfaces: (a) The maximum turbulent kinetic energy varies with the length of Straight-2; (b) The minimum turbulent kinetic energy varies with the length of Straight-2; (c) Comparison of the maximum turbulent kinetic energy obtained by the two models; (d) Comparison of the minimum turbulent kinetic energy obtained by the two models.

Figure 11. The max accretion rate varies with the location of different computation domains: (a) The average particle size is 10 μm; (b) The average particle size is 20 μm; (c) The average particle size is 50 μm.

Figure 12. Comparison of the maximum accretion rate calculated by the DDPM and the DPM: (a) The average particle size is 10 μm; (b) The average particle size is 50 μm.

Figure 13. The erosion contours at Elbow-1 and Straight-2 with the different Straight-2 lengths and different average diameters: (a) The length of Straight-2 is from 1 m to 4 m; (b) the length of Straight-2 is from 5 m to 7 m.

Figure 14. The erosion contours at Elbow-2 and Straight-3 with the different Straight-2 lengths and different average diameters.

Figure 15. The max erosion rate at all computation domains with the different Straight-2 lengths and different average diameters: (a) The average particle size is 10 μm; (b) the average particle size is 20 μm; (c) the average particle size is 50 μm.

Figure 16. Erosion contours calculated by the DDPM and the DPM. Each pair of rendered colors corresponds to the same erosion range.

Figure 17. The max erosion rate calculated by the DDPM and the DPM: (a) Erosion in the elbows; (b) erosion in the straight pipes.

Table 1

Parameters of the impact angle function.

Point	1	2	3	4	5	6	7
Angle	0°	15°	30°	45°	60°	75°	90°
Value	0	1	0.75	0.51	0.18	0.12	0.07

Table 2

The dimensions and units of the three-dimensional geometric models.

Parameters	Value (m)
H1	5.00
H2	5.00
L	1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 7.00
Di	1.00
Do	1.08
Ri	0.50
Ro	1.54

Table 3

Research Data for the location of the maximum erosion zone at Elbow-1.

Research Data	Mazumder et al.	Simulation-DDPM	Simulation-DPM
Max Erosion Location	19°–69°	20°–71°	21°–74°

Table 4

The relations between the mass flow rate of particle flow and average particle diameter.

Average Particle Diameter (μm)	Mass Flow Rate (kg/s)
10	1 × 10−5
20	2 × 10−5
50	5 × 10−5

References

1. Parsi, M.; Najmi, K.; Najafifard, F.; Hassani, S.; McLaury, B.S.; Shirazi, S.A. A comprehensive review of solid particle erosion modeling for oil and gas wells and pipelines applications. J. Nat. Gas Sci. Eng.; 2014; 21, pp. 850-873. [DOI: https://dx.doi.org/10.1016/j.jngse.2014.10.001]

2. Darihaki, F.; Zhang, J.; Vieira, R.E.; Shirazi, S.A. The near-wall treatment for solid particle erosion calculations with CFD under gas and liquid flow conditions in elbows. Adv. Powder Technol.; 2021; 32, pp. 1663-1676. [DOI: https://dx.doi.org/10.1016/j.apt.2021.03.020]

3. Hong, B.; Li, X.; Li, Y.; Li, Y.; Yu, Y.; Wang, Y.; Gong, J.; Ai, D. Numerical Simulation of Elbow Erosion in Shale Gas Fields under Gas-Solid Two-Phase Flow. Energies; 2021; 14, 3804. [DOI: https://dx.doi.org/10.3390/en14133804]

4. Mamoo, M.S.; Goharrizi, A.S.; Abolpour, B. CFD Simulation Of Erosion by Particle Collision in U-Bend and Helical Type Pipes. J. Serbian Soc. Comput. Mech.; 2020; 14, pp. 1-13. [DOI: https://dx.doi.org/10.24874/jsscm.2020.14.02.01]

5. Akrami, S.; Bezuijen, A.; van Beek, V.; Terwindt, J. Influence of a permeable sand layer on the mechanism of backward erosion piping using 3D pipe depth measurements. Acta Geotech.; 2021; 16, pp. 3881-3898. [DOI: https://dx.doi.org/10.1007/s11440-021-01350-2]

6. Agrawal, M.; Khanna, S.; Kopliku, A.; Lockett, T. Prediction of sand erosion in CFD with dynamically deforming pipe geometry and implementing proper treatment of turbulence dispersion in particle tracking. Wear; 2019; 426–427, pp. 596-604. [DOI: https://dx.doi.org/10.1016/j.wear.2019.01.018]

7. Singh, J.; Kumar, S.; Singh, J.P.; Kumar, P.; Mohapatra, S.K. CFD modeling of erosion wear in pipe bend for the flow of bottom ash suspension. Part. Sci. Technol.; 2018; 37, pp. 275-285. [DOI: https://dx.doi.org/10.1080/02726351.2017.1364816]

8. Nam, C.; Lee, Y.; Park, G.; Lee, G.; Lee, W.B. Analysis of Pipe Failure Period Using Pipe Elbow Erosion Model by Computational Fluid Dynamics (CFD). Korean Chem. Eng. Res.; 2018; 56, pp. 133-138. [DOI: https://dx.doi.org/10.9713/kcer.2018.56.1.133]

9. Xu, J.; Lian, Z.; Hu, J.; Luo, M. Prediction of the Maximum Erosion Rate of Gas–Solid Two-Phase Flow Pipelines. Energies; 2018; 11, 2773. [DOI: https://dx.doi.org/10.3390/en11102773]

10. Pei, J.; Lui, A.; Zhang, Q.; Xiong, T.; Jiang, P.; Wei, W. Numerical investigation of the maximum erosion zone in elbows for liquid-particle flow. Powder Technol.; 2018; 333, pp. 47-59. [DOI: https://dx.doi.org/10.1016/j.powtec.2018.04.001]

11. Zhang, J.; McLaury, B.S.; Shirazi, S.A. Modeling sand fines erosion in elbows mounted in series. Wear; 2018; 402–403, pp. 196-206. [DOI: https://dx.doi.org/10.1016/j.wear.2018.02.009]

12. Cloete, S.; Johansen, S.T.; Amini, S. Performance evaluation of a complete Lagrangian KTGF approach for dilute granular flow modelling. Powder Technol.; 2012; 226, pp. 43-52. [DOI: https://dx.doi.org/10.1016/j.powtec.2012.04.010]

13. Pouraria, H.; Darihaki, F.; Park, K.H.; Shirazi, S.A.; Seo, Y. CFD modelling of the influence of particle loading on erosion using dense discrete particle model. Wear; 2020; 460–461, [DOI: https://dx.doi.org/10.1016/j.wear.2020.203450]

14. Sedrez, T.A.; Shirazi, S.A.; Rajkumar, Y.R.; Sambath, K.; Subramani, H.J. Experiments and CFD simulations of erosion of a 90 degrees elbow in liquid-dominated liquid-solid and dispersed-bubble-solid flows. Wear; 2019; 426, pp. 570-580. [DOI: https://dx.doi.org/10.1016/j.wear.2019.01.015]

15. Elemuren, R.; Evitts, R.; Oguocha, I.; Kennell, G.; Gerspacher, R.; Odeshi, A. Slurry erosion-corrosion of 90 degrees AISI 1018 steel elbow in saturated potash brine containing abrasive silica particles. Wear; 2018; 410, pp. 149-155. [DOI: https://dx.doi.org/10.1016/j.wear.2018.06.010]

16. Zhu, H.; Qi, Y. Numerical investigation of flow erosion of sand-laden oil flow in a U-bend. Process Saf. Environ. Prot.; 2019; 131, pp. 16-27. [DOI: https://dx.doi.org/10.1016/j.psep.2019.08.033]

17. Jamwal, A.S. Comparative Study of Particulate Erosion Phenomena Between Elbow and Target Tee in Pipe Flow. Trends in Manufacturing Processes; Springer Science and Business Media LLC: Berlin/Heidelberg, Germany, 2020; pp. 153-159.

18. Vesilind, P. The Rosin-Rammler particle size distribution. Resour. Recover. Conserv.; 1980; 5, pp. 275-277. [DOI: https://dx.doi.org/10.1016/0304-3967(80)90007-4]

19. El-Amin, M.F.; Sun, S.; Heidemann, W.; Muller-Steinhagen, H. Analysis of a turbulent buoyant confined jet modeled using realizable k-E > model. Heat Mass Transf.; 2010; 46, pp. 943-963. [DOI: https://dx.doi.org/10.1007/s00231-010-0625-3]

20. Ogawa, S.; Umemura, A.; Oshima, N. On the equations of fully fluidized granular materials. Zeitschrift für angewandte Mathematik und Physik; 1980; 31, pp. 483-493. [DOI: https://dx.doi.org/10.1007/BF01590859]

21. Djamarani, K.; Clark, I.M. Characterization of particle size based on fine and coarse fractions. Powder Technol.; 1997; 93, pp. 101-108. [DOI: https://dx.doi.org/10.1016/S0032-5910(97)03233-6]

22. Syamlal, M.; Rogers, W.; O’Brien, T.J. MFIX Documentation: Volumel, Theory Guide; National Technical Information Service: Springfield, VA, USA, 1993; DOE/METC-9411004, NTIS/DE9400087, 1993

23. Lun, C.K.K. Kinetic theory for granular flow of dense, slightly inelastic, slightly rough spheres. J. Fluid Mech.; 1991; 233, pp. 539-559. [DOI: https://dx.doi.org/10.1017/S0022112091000599]

24. Forder, A.; Thew, M.; Harrison, D. A numerical investigation of solid particle erosion experienced within oilfield control valves. Wear; 1998; 216, pp. 184-193. [DOI: https://dx.doi.org/10.1016/S0043-1648(97)00217-2]

25. Ou, G.; Bie, K.; Zheng, Z.; Shu, G.; Wang, C.; Cheng, B. Numerical simulation on the erosion wear of a multiphase flow pipeline. Int. J. Adv. Manuf. Technol.; 2017; 96, pp. 1705-1713. [DOI: https://dx.doi.org/10.1007/s00170-017-0834-8]

26. Vieira, R.E.; Mansouri, A.; McLaury, B.S.; Shirazi, S.A. Experimental and computational study of erosion in elbows due to sand particles in air flow. Powder Technol.; 2016; 288, pp. 339-353. [DOI: https://dx.doi.org/10.1016/j.powtec.2015.11.028]

27. Mazumder, Q.H.; Ahmed, K.; Zhao, S. Experimental Investigation of Solid Particle Erosion in S-Bend. J. Fluids Eng.; 2016; 138, 044501. [DOI: https://dx.doi.org/10.1115/1.4031685]

28. Fried, J.; Zietz, S. Curve fitting by spline and akima methods: Possibility of interpolation error and its suppression. Phys. Med. Biol.; 1973; 18, pp. 550-558. [DOI: https://dx.doi.org/10.1088/0031-9155/18/4/306]

29. Ibrahim, H.M.; Awad, M.; Al-Farraj, A.S.; Al-Turki, A.M. Effect of Flow Rate and Particle Concentration on the Transport and Deposition of Bare and Stabilized Zero-Valent Iron Nanoparticles in Sandy Soil. Sustainability; 2019; 11, 6608. [DOI: https://dx.doi.org/10.3390/su11236608]