Numerical Investigation of the Impacts of Large Particles on the Turbulent Flow and Surface Wear in Series-Connected Bends

Abstract

The deep sea harbors abundant mineral, oil, and gas resources, making it highly valuable for commercial development, including the extraction of minerals. Due to the relatively large particle size of these minerals, how they interact with fluids is significantly different from that of small particles. However, there has been limited simulation research on the impacts of large particles (the diameter of particles is at the level of centimeters) on the flow and wear characteristics in bends, because the simulation of the particles at such a size is difficult. Additionally, in the field of deep-sea mining, multiple bends are simultaneously connected in series, and the wear in such bends has garnered increasing attention. Based on an improved CFD-DEM model, this article solved the issue that traditional unresolved CFD-DEM methods cannot accurately simulate large particles in a hydraulic conveying pipe and bend. After validating the accuracy of this model against classical experiments, the paper comprehensively analyzes the modulation effect of large particles on turbulence, and the effects of different particle diameters, particle transport concentrations, and transport velocities on the wear of bends connected serially. Finally, the bends connected serially in various configurations are simulated to study the wear on the bent interior surfaces. Results indicate a pronounced modulation effect of large particles on turbulence at higher transport concentrations; the wear rate in the combined bends does not exhibit a linear correlation with the collision frequency of particles on the wall surface. Furthermore, different configurations of serially connected bends exhibit significant differences in the wear morphology of the second bend.

1. Introduction

With the depletion of terrestrial resources, the deep-sea environment has emerged as a vital successor to resource extraction on land in the 21st century. The deep sea is highly valued for its abundant reserves of minerals, oil, gas, and biological resources, making it part of a significant strategic objective for nations worldwide [1]. Despite its immense commercial potential, exploring and developing deep-sea resources present significant challenges due to the difficulty of exploration and high technological requirements. As such, it has become a focal point of strategic interest for many countries [2]. Currently, the transportation of raw materials relies mainly on pipeline lifting, with pipeline systems and centrifugal pumps serving as the core power devices for the entire deep-sea mining transportation system. The mixed transportation system is mainly used to transport solid–liquid two-phase flow mixtures containing ore particles, posing a significant technical problem due to the wear caused by solid particles on the pipeline walls. The continuous impact and collision of particles with the pipeline walls directly erode the wall, weakening the reliability of the pipeline, increasing hydraulic losses on the wall surfaces, and ultimately reducing the lifespan of relevant mining equipment. This poses substantial safety hazards to the transport system [1]. To address these practical challenges, researchers have investigated the interaction between solid and liquid phases and the wear patterns induced by particle impact on pipeline bends using both experimental and numerical simulation methods.

Since the 1970s, many scholars have conducted extensive experimental research on surface wear in bends that convey particle-laden flows. Vieira et al. [3,4] investigated the wear rate’s influence using particle image velocimetry and slip velocity between air and sand. Zeng et al. [5], using array electrode technology at various locations within a 90° bend conveying solid–liquid two-phase flow, measured pure corrosion rates, wear-enhanced corrosion rates, and corrosion-enhanced wear rates. They also identified the percentages of these four mechanisms, revealing the inherent differences in wear–corrosion synergy at different positions in the bend. These experimental results served as a reference for validating the numerically modeled solid–liquid biphasic flow developed in this study. Wood et al. [6] used non-destructive ultrasonic testing, weight measurements, and micrometer measurements to study the wear characteristics of straight and bent pipes. The results showed significantly greater wear rates on the extrados of the bend and conspicuous wear marks at the bottom of the straight pipe. Sedrez et al. [7] used ultrasonic sensors to measure wall thickness losses at 63 positions along the outer diameter of the bend in the water–sand biphasic flow and dispersed bubble–sand multiphase flow. Bilal et al. [8] conducted a study on the maximum wear positions and wear patterns of bends with different bending angles and curvature radii using the paint method under the conditions of water–sand and water–air–sand flows. However, these studies’ particles range in size from μm to mm, and there has been limited research on the wear of bends caused by ultra-large particles.

Experimental studies offer valuable insights into the relationship between material wear and fluid-particle conditions, but they struggle to capture detailed information about two-phase flows and particle–wall collisions. This makes it difficult to analyze the individual impact of various factors on wear. Fortunately, computational fluid dynamics (CFD) has become an invaluable tool for investigating complex flows thanks to advancements in computer technology. Numerical models for solid–liquid biphasic flows can be categorized as either continuous or discrete, depending on the particle treatment approach. The Two-Fluid Model (TFM) is a popular model of the continuous approach, widely used to study the hydraulic transport of small particles [9,10,11]. However, there is currently no effective TFM model for simulating the hydraulic transport of large particles. Moreover, since the TFM model views particles as a continuous medium, it cannot accurately simulate the interaction between large particles and fluids, or the wear of large particles on walls. Therefore, TFM models are unsuitable for simulating large particle movement in hydraulic conveying systems.

The Eulerian–Lagrangian simulation method, based on the Eulerian approach for describing the fluid phase and the Lagrangian method for describing the particle phase, provides a new avenue for addressing the limitations associated with the Two-Fluid Model (TFM). Common approaches within this framework include the Discrete Phase Model (DPM) and the Discrete Element Method (DEM). The DPM model is capable of tracking the center of mass motion of particles but neglects the interaction forces between particles and the fluid, and it is only applicable to dilute phase flows [12]. Furthermore, although the CFD-DPM coupling method can establish two-way coupling, it fails to consider the influence of fluid volume fractions in the fluid control equations. Wang and Shirazi [13], based on the CFD-DPM model, proposed a wear prediction model for a 90° bend. Chen et al. [14] used the CFD-DPM model to investigate the relative wear severity between blockage in plugged tees and bends in dilute gas–solid two-phase flows. They discussed the impact of gravel loading and phase density on the wear in both the three-way pipe and standard bends. In these studies, the CFD-DPM model was primarily employed for the dilute phase flows with fine particles, typically requiring particle volume fractions of less than 10%. However, in practical processes such as deep-sea mining, where large particles like ores are in the size of centimeters and local particle volume fractions can be significant, considerations of the particle volume fraction and particle–particle collisions become imperative. In such cases, the CFD-DPM model becomes less applicable.

Lately, there has been increased focus on the CFD-DEM coupling method among researchers. Chu and Yu [15,16] have successfully integrated discrete element codes with computational fluid dynamics (CFD) software to simulate complex three-dimensional flow in solid–liquid biphasic flow systems. Zeng et al. [17] utilized the CFD-DEM numerical simulation method to investigate the formation of rare scouring wear patterns, analyzing the influence of secondary flow, vortex intensity, and particle trajectories. Their findings suggested that direct collisions and sliding friction were responsible for the V-shaped wear pattern. Similarly, Zhou et al. [18] conducted an initial investigation into the effects of particle shape and vortex intensity on wear in pneumatic conveying bends using CFD-DEM numerical simulations. Uzi et al. [19] made predictions regarding corrosion in industrial conveying pipelines, highlighting that particle collision velocity and collision angle were crucial factors contributing to wear in bends. Zhao et al. [20] simulated the effects of particle mass load and various interphase forces on the wear of a 90° bend based on the CFD-DEM model, finding that particle mass loading was the most significant parameter causing erosion in the bend, while the lift force exerted on the particle by the fluid was insignificant. However, the pressure gradient force and the virtual mass force should be taken into account in liquid–solid flow and the prediction of bend erosion.

However, the traditional unresolved CFD-DEM method requires the fluid mesh size to be three times or more than the particle size [21]. To tackle this challenge, alternative methods have been proposed such as the overlapping grid method [22], the cubic averaging method [23], and the artificial diffusion method [24]. However, the overlapping grid method can be difficult to implement and the cubic averaging method significantly increases computation time in simulations with a large number of particles [25]. The artificial diffusion method, on the other hand, is equivalent to the Gaussian kernel-based coarse-grained method. By solving a transient diffusion equation of the physical quantity, the excessively large and unreasonable physical quantity is diffused to a larger range in the fluid grid. This, in turn, smooths the quantity in the fluid grid and improves the stability of the algorithm. The artificial diffusion method is easier to implement in commercial solvers and parallel computing compared to the Gaussian kernel-based coarse-grained method. Sun and Xiao [26] successfully used this method to numerically simulate a fluidized bed and obtained convergent results of the coarse particle distribution in refined grids with size ratios ranging from 1/2 to 2. Similarly, Zhou et al. [27] utilized this method to calculate liquid–solid two-phase flow in a vertical pipe transporting coarse particles, and while the results showed convergent results in the refined mesh, significant differences were found between the simulation results and experimental values.

Based on a comprehensive analysis of relevant literature, it is clear that the CFD-DEM approach has shown promise in addressing computational obstacles that arise with the flow of solid–liquid biphasic systems and the wear on bent surfaces in engineering applications. However, there are still some unresolved issues in the study of bend wear. Firstly, due to limitations in experimental conditions, research on the wear of bends caused by large particles in solid–liquid biphasic hydraulic transport systems is limited. Additionally, the traditional unresolved CFD-DEM method is unable to simulate the wear of centimeter-sized particles on elbows and their modulation effect on turbulence. As a result, it is necessary to conduct further research on the wear patterns caused by large particles on curved surfaces and the relationship between particle turbulence interactions using an improved CFD-DEM model. Secondly, understanding the wear characteristics of differently configured serial bends is crucial for designing hydraulic conveying pipelines for deep-sea mining engineering since serial bends are installed in different orientations based on the spatial and structural requirements of the equipment.

This study presents an improved fully coupled CFD-DEM model that solves the issue of traditional unresolved CFD-DEM methods’ inability to accurately simulate the interaction between large particles and fluids. The authenticity of the numerical model is validated by comparing it with typical experimental cases. Further, the study analyzes the interaction between centimeter-sized large particles and turbulence within differently configured serial bends connected with pipes of a diameter of 100 mm. The mathematical principles of the numerical model are detailed in Section 2, while Section 3 outlines the numerical model settings and provides comparisons with typical experiments. Section 4 systematically analyzes the influence of various factors on bend wear, as well as their effects on the turbulent kinetic energy and turbulent dissipation rate within the pipeline. The study’s main conclusions are summarized in the last section.

2. Methodology

2.1. Control Equation for Liquid Phase

The fluid is treated as an incompressible continuous medium, and the Reynolds-averaged Navier–Stokes equations are employed for solving. The mass and momentum conservation equations for the fluid are respectively expressed as

$$\frac{\partial }{\partial t}\left({\alpha }_f{\rho }_f\right)+\nabla \cdot \left({\alpha }_f{\rho }_f{\overline{\mathbf{v}}}_f\right)=0$$

(1)

$$\begin{array}{c}\frac{\partial }{\partial t}\left({\alpha }_f{\rho }_f{\overline{v}}_f\right)+\nabla \cdot \left({\alpha }_f{\rho }_f{\overline{v}}_f{\overline{v}}_f\right)=-\nabla p+\nabla \cdot \left[{\alpha }_f({\mu }_f+{\mu }_t)(\nabla {\overline{v}}_f+\nabla {\overline{v}}_f{}^T)\right]\\ +{\alpha }_f{\rho }_f\mathbf{g}+F_{p-f}\end{array}$$

(2)

$${\mathbf{F}}_{p-f}=-\frac{1}{V_{cell}}{\displaystyle {\sum }_{i=1}^n{\mathbf{F}}_{f-p}}$$

(3)

In the equation, ρ_f is the density of the fluid, ${\overline{v}}_f$ is the average velocity of the fluid, p is the pressure, μ_f and μ_t are the viscosity and turbulent viscosity of the fluid, α_f represents the volume fraction of the liquid phase. In Equation (3), n is the number of particles in the calculated CFD cell and $V_{cell}$ is the volume of the CFD cell. F_f−p is the total force that is enforced in the fluid by the particles and F_f−p is the fluid-induced force exerted on each particle.

The particles are divided into several segments to compute the particle volume fraction within the fluid grid [28,29], expressed by the following equation:

$${\alpha }_f=1-\frac{{\displaystyle {\sum }_{i=1}^k{V}_{i,part}}}{V_{cell}}$$

(4)

In the equation, $V_{i,part}$ represents the total volume of the discrete portion of particle i in the computational grid, and k denotes the number of particles overlapping with the computational grid. Through a discretization approach, the volume of large particles is distributed proportionally into the CFD grid occupied by the particles. Subsequently, the artificial diffusion method, employing an equivalent statistical Gaussian kernel function, is utilized to compute the diffusion smoothing field for the initial values of the particle volume fraction calculated using Equation (4). In comparison to traditional unresolved methods and statistical Gaussian kernel-based coarse-graining approaches, this method applies to various grid sizes, features simplified boundary condition handling, and is conducive to parallel computing. The principle is shown in Figure 1. The equation for the diffusion algorithm is as follows:

$$\frac{{{\displaystyle \partial \phi }}_{}}{{{\displaystyle \partial t}}_{}}={{\displaystyle \nabla }}^2{{\displaystyle \phi }}_{}\hspace{1em}X\in \mathsf{\Omega }, t>0$$

(5)

$$\phi (X,t=0)={\phi }_0(X)$$

(6)

where φ is the diffused variable, which could be the solid volume fraction α_p, the momentum source term F_f−p, and the source terms of turbulence S_k and S_ε. Ω denotes the boundary of the solved domain. φ₀ is the initial value of the diffused variable. The non-flux boundary condition is defined at the domain boundary to ensure the variable conservation. The pseudo-time τ controls the scale of the diffusion and computes the diffused variables until t = τ [26].

The focus of this study is on the elbow in solid–liquid two-phase flow transport. The elbow induces secondary flow. Considering that the RNG k-ε turbulence model corrects for vortices in turbulence, this study employs the RNG k-ε turbulence model. The transport equations for turbulent kinetic energy and turbulence dissipation are expressed as follows:

$$\frac{\partial }{\partial t}\left({\alpha }_f{\rho }_f{k}_f\right)+\frac{\partial }{\partial x_i}\left({\alpha }_f{\rho }_f{k}_f{u}_i\right)=\frac{\partial }{\partial x_j}\left[{\alpha }_f{\alpha }_k{\mu }_{eff}\frac{\partial k_f}{\partial x_j}\right]+G_k+G_b-{\alpha }_f{\rho }_f{\epsilon }_f+S_k$$

(7)

$$\begin{array}{c}\frac{\partial }{\partial t}\left({\alpha }_f{\rho }_f{\epsilon }_f\right)+\frac{\partial }{\partial x_i}\left({\alpha }_f{\rho }_f{\epsilon }_f{u}_i\right)=\frac{\partial }{\partial x_j}\left[{\alpha }_f{\alpha }_{\epsilon }{\mu }_{eff}\frac{\partial {\epsilon }_f}{\partial x_j}\right]+{\alpha }_f{C}_{1\epsilon }\frac{{\epsilon }_f}{k_f}\left(G_k+C_{3\epsilon }{G}_b\right)\\ -{\alpha }_f{C}_{2\epsilon }{\rho }_f\frac{{\epsilon }_f^2}{k_f}+S_{\epsilon }\end{array}$$

(8)

where k_f and ε_f are the turbulent kinetic energy and dissipation rate, respectively. G_k is the production of turbulence kinetic energy by the velocity gradient. G_b is related to the buoyancy caused by the temperature gradient and it is zero in the isothermal condition of this work. σ_k = 1.0 and σ_ε = 1.3 are the turbulent Prandtl numbers for k_f and ε_f, respectively. C_1ε = 1.44 and C_2ε = 1.92 are two constants. S_k and S_ε represent the source terms of the turbulence kinetic energy and turbulence dissipation rate caused by the presence of particles in the fluid. According to the work of Crowe et al. [30], the influence of particles on the turbulence kinetic energy and turbulence dissipation rate can be represented as

$$S_k=\frac{1}{V_{cell}}{{\displaystyle \sum _{n}3\pi {\mu }_f{d}_p{f}_k\left|v_f-v_p\right|}}^2$$

(9)

$$S_{\epsilon }={\rho }_f{C}_{\epsilon 3}{}^{\text{'}}\frac{{\nu }^2}{V_{cell}}{\displaystyle \sum _n{f}_k}\frac{{\left|v_f-v_p\right|}^2}{d_p}$$

(10)

where v_f and v_p represent the momentary velocities of the fluid and particles, respectively. The drag factor is represented by f_k, which is defined as $f_k=1+0.15{\mathrm{Re}}_p{}^{0.687}$. $C_{\epsilon 3}{}^{\text{'}}=C_{\epsilon 2}^{\text{'}}{C}_{\epsilon p}{\mathrm{Re}}_p{}^{1.416}$, $C_{\epsilon 2}^{\text{'}}\approx C_{\epsilon 2}=1.92$ and $C_{\epsilon p}=0.058$. With the help of a user-defined function in Fluent, the source terms were defined in the interface code and linked to the turbulence equations.

2.2. Control Equation for Solid Phase

The Discrete Element Method (DEM), introduced by Cundall and Starck [31], describes particles as discrete entities and their motion in two types—translational and rotational. Newton’s second law governs both types of motion, while forces acting on each particle and their trajectories are tracked in Lagrangian coordinates. The motion equation for each particle is expressed as follows:

$$m_p\frac{\mathrm{d}{v}_p}{dt}=m_p\mathbf{g}+{\mathbf{F}}_{drag}+{\mathbf{F}}_{lift}+{\mathbf{F}}_p+{\mathbf{F}}_{vm}+{\mathbf{F}}_b+{\displaystyle \sum {\mathbf{F}}_{c,ij}}$$

(11)

In the equation, m_p is the mass of individual particles, and g is the gravitational acceleration. F_drag, F_lift, F_p, and F_vm are the drag force, the lift force, the pressure gradient force, and the virtual mass force, respectively. F_b is the buoyant force, which is calculated as ${\mathbf{F}}_b=-{\rho }_f{V}_p\mathbf{g}$. The contact force between particles (or particles and wall) is represented by F_c,ij, which can be obtained from the particle contact model.

A natural aspect of particle motion is particle rotation, particularly noticeable for heavier or larger particles with high moments of inertia. Thus, in this model, we consider the impacts of rotation caused by tangential forces, rolling friction, and fluid-induced rotation in the particle torque equation. The equation for the torque of an individual particle is expressed as follows:

$$I_p\frac{\mathrm{d}{\mathsf{\omega }}_p}{\mathrm{dt}}={\displaystyle \sum T_c}+T_f$$

(12)

where ω_p is particle angular velocity. The moment of inertia for a sphere is represented by I_p, which is defined as $I_p=0.1m_p{d}_p^2$. T_c refers to the contact torque between particles, while T_f refers to the torque caused by the fluid.

(1)

Drag force model

The force of drag is a crucial factor in the interaction between particles and fluid at the interphase. In situations where fluid and multiphase particles are coupled, this study utilizes the drag force model put forth by Di Felice. This model is calculated by multiplying the drag force experienced by particles in the absence of fluid obstruction under identical volume flow rates with the void fraction function [32]. The formula is expressed as follows:

$$F_{drag}=\frac{1}{2}{\rho }_f{C}_D\frac{\pi d_p{}^2}{4}\left|{\overline{v}}_f+v_f^{\text{'}}-v_p\right|\left({\overline{v}}_f+v_f^{\text{'}}-v_p\right)\left({\alpha }_f^{1-\gamma }\right)$$

(13)

(2)

Lift force model

Currently, the most commonly utilized lift force models are the Saffman and Magnus models, which originated from gas–solid two-phase flow studies [33,34]. However, recent research by Zhao et al. [35] has revealed drawbacks of using Saffman and Magnus lift models in solid–liquid two-phase flow. In this study, we build upon the findings of Loth et al. [36] and modify the lift coefficient induced by rotation and shear through superimposition. Additionally, we take into account the contributions of fluid vorticity and particle rotation to lift, which can be expressed as follows:

$$F_{lift}=\frac{1}{8}\pi d_p^2{\rho }_f{C}_L\left|{\overline{v}}_f+v_f^{\text{'}}-v_p\right|\left[\left({\overline{v}}_f+v_f^{\text{'}}-v_p\right)\times \frac{{\omega }_i}{\left|{\omega }_i\right|}\right]$$

(14)

In the equation, C_L represents the lift coefficient, and its expression is as follows:

$$C_L=J^{\ast }\frac{12.92}{\pi }\sqrt{\frac{{\omega }^{\ast }}{R{e}_p}}+{\mathsf{\Omega }}_{p,eq}^{\ast }{C}_{L,\mathsf{\Omega }}^{\ast }$$

(15)

The remaining expressions containing parameters are as follows:

$$J^{\ast }=0.3\left\{1+\tanh \left[\frac{5}{2}\left({\log }_{10}\sqrt{\frac{{\omega }^{\ast }}{R{e}_p}}+0.191\right)\right]\right\}\left\{\frac{2}{3}+\tanh \left[6\sqrt{\frac{{\omega }^{\ast }}{R{e}_p}}-1.92\right]\right\}$$

(16)

$$C_{L,\mathsf{\Omega }}^{\ast }=1-\left\{0.675+0.15\left\{1+\tanh \left[0.28\left({\mathsf{\Omega }}_p^{\ast }-2\right)\right]\right\}\right\}\tanh \left(0.18R{e}_p^{1/2}\right)$$

(17)

$${\mathsf{\Omega }}_{p,eq}^{\ast }=\frac{{\omega }^{\ast }}{2}\left(1-0.0075R{e}_{\omega }\right)\left(1-0.062R{e}_p^{1/2}-0.001R{e}_p\right)$$

(18)

$${\omega }^{\ast }=\frac{\left|{\mathsf{\omega }}_f\right|d_p}{\left|v_f-v_p\right|}$$

(19)

$${\mathsf{\Omega }}_p^{\ast }=\frac{\left|{\mathsf{\Omega }}_p\right|d_p}{\left|v_f-v_p\right|}$$

(20)

$$R{e}_{\omega }=\frac{{\rho }_f{d}_p^2\left|{\mathsf{\omega }}_f\right|}{{\mu }_f}$$

(21)

(3)

Models of virtual mass and pressure gradient forces

Virtual mass force is the force required to accelerate the fluid around a particle. The virtual mass of a fluid is the mass of the fluid surrounding a particle calculated based on its acceleration. The expression for the virtual mass force is as follows [37]:

$$F_{vm}=C_{vm}{m}_p\frac{{\rho }_f}{{\rho }_p}\left[v_p\nabla \left({\overline{v}}_f+v_f^{\text{'}}\right)-\frac{d{v}_p}{dt}\right]$$

(22)

Pressure gradient force refers to the additional force generated by particles in a fluid with a pressure gradient, which is opposite to the direction of the fluid pressure gradient [37]. The expression is as follows:

$$F_p=-V_p\nabla p$$

(23)

(4)

Torque model

In Equation (10), according to the research by Rubinow and Keller [38], the expression for the rotational torque T_f exerted by the fluid on the particle is

$$T_f=\frac{{\rho }_f}{2}{\left(\frac{d_p}{2}\right)}^5{C}_R\left|{\mathsf{\omega }}_r\right|{\mathsf{\omega }}_r$$

(24)

In the equation, ${\mathsf{\omega }}_r=0.5\nabla \times v_f-{\mathsf{\Omega }}_p$, C_R is the rotational torque coefficient, and based on the research by Rubinow and Keller [38] as well as the direct numerical simulation study by Dennis et al. [39], the expression for the rotational torque coefficient C_R is as follows:

$$C_R=\left\{\begin{array}{ll}\frac{64\pi }{R{e}_r}\hfill & \hfill R{e}_r\le 32\hfill \\ \frac{12.9}{\sqrt{R{e}_r}}+\frac{128.4}{R{e}_r}\hfill & \hfill 32<R{e}_r<1000\hfill \end{array}\right.$$

(25)

The contact torque T_c between two particles, i and j, has two components. One is the torque T_t,ij that is generated by the tangential contact force, and the other is the rolling friction torque T_r,ij. The formula to calculate the torque generated by the tangential contact force is as follows:

$${\mathbf{T}}_{t,ij}=\mathbf{r}\times F_{ct,ij}$$

(26)

The formula used to calculate the rolling friction torque is as follows [40]:

$${\mathbf{T}}_{r,ij}=-{\mu }_r{F}_{cn,ij}\left|\mathbf{r}\right|\frac{\mathsf{\omega }}{\left|\mathsf{\omega }\right|}$$

(27)

where μ_r is the rolling friction coefficient, and $F_{cn,ij}$ is the normal contact force, which can be found in Equation (28).

(5)

Particle contact model

The DEM software EDEM utilizes its default model to directly simulate collisions. The resulting contact force, Fc, takes into account both normal and tangential components between two particles or between a particle and a wall. This is achieved through the implementation of the soft-sphere contact model, originally developed by Cundall and Strack [31]. The contact force is further composed of both normal and tangential forces, which are based on Hertzian contact theory and Middlin-Deresiewicz’s research [41], respectively. Study [42] is referenced to accurately calculate the normal and tangential damping forces.

$$F_{c,ij}=F_{cn,ij}+F_{ct,ij}$$

(28)

$$F_{cn,ij}=-k_{n,ij}{\delta }_{n,ij}{\mathbf{n}}_i-{\gamma }_{n,ij}({\mathbf{v}}_r\cdot {\mathbf{n}}_i){\mathbf{n}}_i$$

(29)

$$F_{ct,ij}=-k_{t,ij}{\delta }_{t,ij}{\mathbf{t}}_i-{\gamma }_{t,ij}\left[({\mathbf{v}}_r\cdot {\mathbf{t}}_i){\mathbf{t}}_i+({\mathsf{\omega }}_i\times {\mathbf{r}}_i-{\mathsf{\omega }}_j\times {\mathbf{r}}_j)\right]$$

(30)

The equation for calculating the normal and tangential unit vectors is represented by n and t, respectively. The relative velocity vector between particle i and particle j during their collision is denoted by v_r = v_i − v_j. The elastic stiffness, damping coefficient, and overlapping amount of the particle contact are represented by k, γ, and δ, respectively. This model can also be extended to the collision between a particle and a wall. When rolling friction occurs, the Coulomb friction model can be used to calculate the tangential contact force, which can be expressed as follows:

$$F_{ct,ij}=-{\mu }_s\left|F_{cn,ij}\right|t_i$$

(31)

Within the equations, μ refers to the sliding friction coefficient. Additionally, k_n,ij denotes the normal spring stiffness, γ_n,ij represents the normal damping coefficient, k_t,ij signifies the tangential spring stiffness, γ_t,ij denotes the tangential damping coefficient, δ_n,ij is the amount of overlap in the normal direction, δ_t,ij is the amount of overlap in the tangential direction, R^* denotes the equivalent radius of the particle, Y^* signifies the equivalent Young’s modulus of the particle, G^* denotes the equivalent shear modulus of the particle, and S_n,ij and S_t,ij denote the normal and tangential stiffness, respectively. The specific expressions for these pertinent parameters are displayed in

Table 1. The formula for collision model variables.

$k_{n,ij}=\frac{4}{3}{Y}^{*}\sqrt{R^{*}{\delta }_{n,ij}}$		$k_{t,ij}=8G^{*}\sqrt{R^{*}{\delta }_{t,ij}}$
${\gamma }_{n,ij}=2\sqrt{\frac{5}{6}\beta }\sqrt{S_{n,ij}m*}$		${\gamma }_{t,ij}=2\sqrt{\frac{5}{6}\beta }\sqrt{S_{t,ij}m*}$
$\frac{1}{Y^{*}}=\frac{(1-{\nu }_i^2)}{Y_i}+\frac{(1-{\nu }_j^2)}{Y_j}$		$\frac{1}{G^{*}}=\frac{2(2+{\nu }_i)(1-{\nu }_j)}{Y_i}+\frac{2(2+{\nu }_i)(1-{\nu }_j)}{Y_j}$
$S_{n,ij}=2Y^{*}\sqrt{R^{*}{\delta }_{t,ij}}$		$S_{t,ij}=8G^{*}\sqrt{R^{*}{\delta }_{t,ij}}$
$\frac{1}{R^{*}}=\frac{1}{R_i}+\frac{1}{R_j}$	$\frac{1}{m^{*}}=\frac{1}{m_i}+\frac{1}{m_j}$		$\beta =\frac{ln(e)}{\sqrt{l{n}^2(e)+{\pi }^2}}$

.

2.3. Discrete Random Walk Model

To incorporate turbulent velocity fluctuations into interphase force models and turbulent modulation source terms, the discrete random walk (DRW) model [43] is utilized. This model uses a numerical approach and considers the instantaneous velocity of the fluid $v_f$ within the vortex, which includes both mean velocity ${\overline{v}}_f$ and velocity fluctuation $v_f^{\text{'}}$:

$$v_f={\overline{v}}_f+v_f^{\text{'}}$$

(32)

The fluctuation of turbulent pulsation velocity follows a Gaussian distribution, which is obtained using the following formula:

$$\left\{\begin{array}{c}u\text{'}={\xi }_i\sqrt{\overline{u{\text{'}}^2}}\\ v\text{'}={\xi }_j\sqrt{\overline{v{\text{'}}^2}}\\ w\text{'}={\xi }_k\sqrt{\overline{w{\text{'}}^2}}\end{array}\right.$$

(33)

where ${\xi }_i, {\xi }_j, {\xi }_k$ are zero-mean Gaussian random variables with a standard deviation of 1 and the remaining portion on the right represents local root mean square (RMS) velocity fluctuations, represented as

$$\sqrt{\overline{u{\text{'}}^2}}=\sqrt{\overline{v{\text{'}}^2}}=\sqrt{\overline{w{\text{'}}^2}}=\sqrt{2\frac{k_f}{3}}$$

(34)

where k_f represents the fluid turbulent kinetic energy, solved according to the standard k-ε equation.

2.4. Erosion Prediction Model

Currently, numerous scholars have created various wear models, usually developed through a regression analysis based on experiments. Each wear rate in classical wear models has different physical definitions, encompassing various influencing factors. In solid–liquid two-phase flow, numerous factors affect wear, such as material properties, particle characteristics, and liquid-phase factors. Zhang et al. [44] conducted experimental research and numerical simulations on particle velocity near the target material and wear at the elbow of a 90° bend in water and airflow. The study results indicate that the E/CRC wear model and Oka’s wear model were originally developed for particles in high-speed airflow (>50 m/s). However, their application to wall wear in elbow hydraulic transport under low-flow conditions demonstrates their suitability and accurate prediction of wall wear. Peng and Cao [45] combined five wear models and two particle–wall rebound models to predict the wear rate in solid–liquid two-phase flow bends. The research results suggest that the E/CRC wear model is closer to experimental data. Therefore, this section will mainly introduce the E/CRC wear model, and the next section will discuss the validation of the E/CRC wear model. The E/CRC wear model is defined in the particle–wall contact model based on the DEM solver’s Application Programming Interface (API) development functionality, with the specific expression as follows:

$$ER=C{\left(BH\right)}^{-0.59}{F}_s{v}_p^{n}F\left(\theta \right)$$

(35)

$$F\left(\theta \right)={\displaystyle \sum _{i=1}^5{A}_i{\theta }^i}$$

(36)

In this equation, ER denotes the wear rate, which is defined as the mass lost by the material due to particle impact divided by the mass of the impacting particle. BH represents the material’s Brinell hardness, while $F_s$ is the particle shape factor. For circular, semi-circular, and sharp-edged particles, $F_s$ takes values of 0.2, 0.53, and 1, respectively. The particle impact velocity, measured in m/s, is denoted by v_p, while θ represents the collision angle of the particles in radians. The constants used in the equation are n = 2.41 and C = 2.17 × 10⁻⁷, and the values of Ai (i = 1–5) are 5.3983, −10.1068, 10.9327, −6.3283, and 1.4234, respectively.

2.5. Geometric Model and Boundary Condition Setting and Numerical Method

For this study, a standard elbow with an inner diameter (D) of 100 mm was selected. To ensure optimal fluid and particle development, the inlet straight pipe and outlet pipe were set to 20D in length. The curvature radius to diameter ratio (R/D) was set at 1.5, and the elbow bend angles ranged from 0° to 90°, evenly distributed. The geometric model is depicted in Figure 2, with the first elbow arranged from vertically upwards to horizontally, the second elbow in case 1 arranged downwards from horizontal to vertical, the second elbow in case 2 arranged from horizontal to horizontal, and the second elbow in case 3 arranged upwards from horizontal to vertical. The gravitational direction follows the negative Z-axis and the flow direction is indicated by the arrows in the diagram. To divide the mesh, the O-type grid method was used with a boundary layer mesh y+ around 30, which ensured the optimal application of the standard wall function. In EDEM, particles are randomly generated at any position at the inlet cross-section of the geometric model through the particle factory. The volume flow rate (mass flow rate) of the hydraulic system per unit time multiplied by the corresponding particle concentration determines the volume (mass) of particles generated per unit time. The particles are immediately eliminated upon exiting the pipe outlet. The initial velocity of the particles matches the uniform inlet velocity of the fluid. The simulated long pipe allows for the full development of the mixed flow. The inlet velocities were set at 5 m/s and 10 m/s, with the outlet set as a pressure outlet at standard atmospheric pressure. Non-slip walls were assumed for the entire computational domain. A Hertz–Mindlin (no-slip) model was used for particle–particle contacts, and a Hertz–Mindlin (no-slip) model with an added wear model was used for particle–wall contacts. Refer to

Table 2. Particles and wall parameters.

	Parameter	Unit	Value
Particle	Density	kg/m3	2650
	Diameter	mm	10, 15, 20
	Poisson’s Ratio		0.3
	Young’s Modulus	Pa	1 × 108
	Incident Velocity	m/s	5, 10
	Coefficient of Restitution		0.8
	Coefficient of Static Friction		0.2
	Coefficient of Rolling Friction		0.01
Wall	Density	kg/m3	7800
	Poisson’s Ratio		0.3
	Young’s Modulus	Pa	2 × 1011
	Coefficient of Restitution		0.5
	Coefficient of Static Friction		0.2
	Coefficient of Rolling Friction		0.01

for the characteristics of the particles and walls.

The ANSYS Fluent software uses the finite volume method to solve the liquid phase’s governing equations, while the EDEM software solves the solid phase equations based on DEM. To model turbulence in the bend, we selected the RNG k-ε turbulence model due to its correction for swirling flow in the model. Our team has developed an in-house interface code that couples the calculations between the two tools. This interface program includes the diffusion equation and ensures the movement of large particles in the hydraulic system by smoothing the solid volume fraction and other parameters. Our code also considers the influence of turbulence on particle motion, with the instantaneous velocity of the fluid composed of the mean velocity from the in-house interface code and turbulent velocity fluctuation. The fluid forces acting on the particle are calculated using the instantaneous fluid velocity, to express the influence of turbulence on particle motion. The influence of particles on fluid turbulence is considered by adding appropriate source terms defined in Equations (7) and (8) in the transport equations of the turbulence. Through the above methods, we achieved bidirectional coupling between particles and turbulence. The data transmission between both software is roughly as follows:

The fluid field is solved using ANSYS Fluent, and the resulting information, including grid data, fluid velocity, and pressure, is then input into EDEM to calculate the fluid forces acting on each particle. Solve each particle using Newton’s second law. To simulate particle motion, we follow the acted forces and record the new particle data for the next calculation step. In the interface code, we survey the particles simulated in EDEM to locate their positions in the CFD grids. When looping over the particles, we calculate the interaction forces between the particle and fluid and sum the forces from each particle to the located CFD cell as the force source. The same scheme is applied for S_k and S_ε, and all these terms are treated as the body source in the calculation.

In the CFD solver, all variables are saved in the center of the cells since it employs the finite volume method. When integrating convective and diffusive terms, we interpolate the variable at the interfaces between two fluid cells and use it to calculate fluxes. The coupling of velocity and pressure terms is accomplished using the SIMPLE algorithm, and a second-order upwind scheme is employed to discretize the momentum, turbulent kinetic energy, and turbulent dissipation rate. The variable values in the center of the cell are also used to calculate source terms, such as the volumetric interaction force and turbulence sources.

Fluent and EDEM utilize different time steps for transient simulation, with the former at 5 × 10⁻⁴ s and the latter at 5 × 10⁻⁵ s. The simulation convergence default criterion for each time step is 10⁻⁴ for mass, momentum, and turbulent parameters, with a maximum iteration number of 60 for each.

2.6. Verification of Grid Independence and Statistical Time Independence

For the case geometric model, this paper considers three mesh schemes to grid the pipes: fine mesh (296,176 cells), medium mesh (106,623 cells), and coarse mesh (62,197 cells). Figure 3a displays the time-average particle impact frequency on the wall as predicted using the different mesh schemes. The medium grid scheme is selected to ensure accuracy and efficiency in all subsequent simulations. Figure 3b demonstrates the influence of statistical time on the time average result, with data after 2.0 s being used for statistics. It is worth noting that the simulation’s physical time was 6s to ensure that the liquid–solid flow was established for a statistical analysis in each simulation. All results are time averages over the simulation time (excluding the initial 2 s). After 4 s, the parameters stabilize, and there is minimal difference in the average collision frequency of particles on the curved pipe wall. Consequently, the following simulation’s duration is set at 4 s.

3. Model Validation with the Experiments

3.1. Validation of Solid–Liquid Two-Phase Flow Characteristics

To guarantee the precision of the simulation, it is crucial to validate the numerical model through experimental data. For this particular study, the experimental data sourced from Vlasák et al. [46] and Zeng et al. [5] were selected as the reference points for comparison, thereby confirming the fidelity of the numerical simulation.

3.1.1. Validating the Flow Properties of 11 mm Particles inside a 4000 mm Horizontal Pipe

Vlasák conducted experiments on the transportation of 11mm particles with an initial velocity of 4.1 m/s and initial concentrations of 8.2% and 13.4% in a 4000-mm-long horizontal pipe with a 100 mm inner diameter. We replicated the experiment’s conditions in our study, and

Table 3. Particles and wall parameters used to verify the flow characteristics of solid-liquid two-phase flow in horizontal and vertical pipes.

	Parameter	Unit	Value
Particle	Density	kg/m3	2650
	Diameter	mm	11
	Poisson’s Ratio		0.3
	Young’s Modulus	Pa	1 × 108
	Incident Velocity	m/s	4.1, 2.05
	Coefficient of Restitution		0.8
	Coefficient of Static Friction		0.2
	Coefficient of Rolling Friction		0.01
Wall	Density	kg/m3	7800
	Poisson’s Ratio		0.3
	Young’s Modulus	Pa	2 × 1011
	Coefficient of Restitution		0.5
	Coefficient of Static Friction		0.2
	Coefficient of Rolling Friction		0.01

presents the simulation parameters. Figure 4 shows the horizontal pipe’s geometric model and the cross-section mesh. Cross-section data processing occurred at a location 3500 mm from the pipe inlet. As shown in Figure 4c, the values of the parameters of the fluid and particle were extracted from 50 sampling points uniformly distributed along a horizontal line. A total of 10 horizontal lines with an equivalent interval were created along the height of the pipe. Therefore, the parameters’ values were averaged from the 50 sampling points at each horizontal line. Figure 5 compares the simulation and experimental results, demonstrating excellent agreement in particle concentration distribution for different initial volume fractions. The maximum error is within 20%.

3.1.2. Validating the Flow Properties of 11 mm Particles inside a 6000 mm Vertical Pipe

In this study, the movement of 11mm particles in a vertical pipe was examined and compared to experimental data. The pipe’s length and diameter are 6000 mm and 100 mm, respectively, as shown in Figure 6. The data processing for the cross-section was conducted at a location 5500 mm from the pipe inlet. As shown in Figure 6c, the values of the parameters of the fluid and particle were extracted from 50 sampling points uniformly distributed along the pipe radius and 50 sampling points uniformly distributed along the azimuthal direction of the cross-section. The radial distribution of the parameter was obtained by averaging the 50 sampling points along the azimuthal direction at each radial position. The transport velocity was 2.05 m/s, and the total particle volume fraction was 4.5%. Table 3 provides the simulation parameters.

The experiment referenced here presents the distributions of the particle volume fraction on the cross-section of data processing, as shown in Figure 7, alongside the simulation outcomes. The observations reveal that the particles are distributed circumferentially across the cross-section, with the highest particle volume fraction not found at the center, but rather more likely concentrated in an area of an annular ring shape. The simulation results also demonstrate this feature, with particles aggregating in an annulus. The maximum value observed in the simulation is 10%, which is consistent with the experimental results. Furthermore, the radial distributions of the particle volume fraction in both the experiment and simulation are depicted in Figure 8, showing that the radial distribution of the particle volume fraction in the experimental results is accurately reproduced by the simulation, thus validating its effectiveness.

3.1.3. Validation of the Prediction of Bend Erosion

In a study conducted by Zeng et al. [5], various experimental analyses were performed to examine the pure corrosion rate, pure wear rate, wear rate under corrosion enhancement, and corrosion rate under wear enhancement in a 90° elbow. The elbow had an internal diameter of D = 50 mm and a curvature radius to diameter ratio of 1.5. To create a well-developed two-phase flow velocity field, a 20 D horizontal pipe was installed in front of the bend. The downstream pipe length was set to 10 D to eliminate the influence of secondary flow in the bend. The geometric model of the simulated bend and the computational domain grid are illustrated in Figure 9. Using a fully coupled CFD-DEM numerical model, the wear process of spherical particles with a diameter of 450 μm and a concentration of 1.2% in the bend was numerically simulated. The inlet flow velocity of the bend was set to 4 m/s, with a hydraulic radius of 50 mm, and the outlet of the bend was set as a pressure outlet boundary condition. To account for the generation of vortices in the bend, the RNG k-ε turbulence model was employed. The DEM solver utilized the Hertz–Mindlin (no slip) contact model, and the E/CRC wear model was utilized to calculate the wall wear rate. The simulation parameters are listed in

Table 4. Particles and wall parameters used to verify wear of bent pipes.

	Parameter	Unit	Value
	Density	kg/m3	2650
	Diameter	mm	0.45
	Poisson’s Ratio		0.23
Particle	Young’s Modulus	Pa	5.9 × 1010
	Mass Flow Rate	kg/s	0.235
	Density	kg/m3	8200
	Poisson’s Ratio		0.3
Wall	Young’s Modulus	Pa	2.07 × 1011
	Hardness	GP	3.43

.

In Figure 10, the wear rates on the outer convex inner wall of the elbow are shown. The curve at an azimuth angle of 180° corresponds to the outermost inner wall of the elbow’s outer convex side. To determine the accuracy of the predicted wear rates, they were compared to the experimental results of Zeng et al. [5]. The research results confirm that the E/CRC wear model is precise in predicting the wear rate on the inner side of the outer convex of the elbow. The wear rate gradually increases as the bending angle of the elbow increases, with the maximum wear occurring near the exit of the elbow. As a result, the E/CRC wear model was adopted in this study to predict the wear rate on the elbow’s inner wall.

4. Results and Discussion

4.1. Influence of Large Particles on Turbulence in Solid–Liquid Two-Phase Interactions

To comprehend the impact of particles on turbulence modulation, one can analyze the dispersion of the turbulence kinetic energy (k) and turbulence dissipation rate (e) across a sequence of cross-sections established in the elbow. The figure below illustrates a comparison of wear maps, turbulence kinetic energy, and turbulence dissipation rate maps under varying operating conditions, with legend values confined to the same range.

4.1.1. Fifteen-Millimeter Particles at 5% and 5 m/s

Comparing the wear maps of 15 mm particles at a 5% transport concentration and a transport velocity of 5 m/s with and without considering turbulent modulation, as is shown in Figure 11, it can be observed that the wear morphology of the first elbow is relatively concentrated when turbulent modulation is considered, whereas without turbulent modulation, the maximum wear occurs near the outlet. For the second elbow, it is evident that the wear rate of particles on the wall is higher when turbulent modulation is not considered.

In the case of 15 mm particles at a transport concentration of 5% and a transport velocity of 5 m/s, the maximum values of k occur near the outer wall, while the maximum values of e occur near the inner and outer walls. As is shown in Figure 12 and Figure 13, the simulation considering turbulent modulation shows higher values of the turbulent kinetic energy (k) and turbulent dissipation rate (e) on the outer wall and near the pipe center compared to the case where turbulent modulation is neglected. This phenomenon is similar to the findings of Hetsroni [47], who proposed that particles with lower particle Reynolds numbers weaken turbulence. When particles have higher particle Reynolds numbers, the shedding of vortices from the particle tails enhances turbulence. In this study, considering the centimeter-scale particle size with a significant particle Reynolds number, the modulation of turbulence by large particles results in a noticeable increase in the turbulent kinetic energy and turbulent dissipation rate in the outer wall region and near the center of the pipe.

4.1.2. Twenty-Millimeter Particles at 10% and 5 m/s

Comparing the wear maps of 20 mm particles at a concentration of 10% and a conveying velocity of 5 m/s with and without considering turbulent modulation, as is shown in Figure 14, it is evident that, when considering turbulent modulation, the wear concentration in the first bend occurs around 60° from the inlet. In contrast, without considering turbulent modulation, wear in the first bend is more pronounced around 60° from the inlet and near the outlet. Regarding the second bend, it is noticeable that without considering turbulent modulation, the wear rate of particles on the bend surface is significantly higher than when turbulent modulation is considered.

In the case of 20 mm particles at a concentration of 10% and a conveying velocity of 5 m/s, the modulation effect of particles on the turbulent kinetic energy and turbulent dissipation rate becomes more pronounced, as shown in Figure 15 and Figure 16. This observation is similar to the findings of Rashidi et al. [48], who experimentally studied the interaction between particles and turbulence in the wall boundary layer under different particle sizes, particle densities, particle mass fractions, and fluid Reynolds numbers. Their research indicated that large-scale polystyrene particles increase the number of wall-induced vortices, thereby enhancing turbulence and Reynolds stress. On the other hand, small-scale polystyrene particles decrease the number of wall-induced vortices, leading to a weakening of turbulence and Reynolds stress. Moreover, this enhancement or attenuation effect is strengthened with an increase in the particle mass fraction.

4.1.3. Twenty-Millimeter Particles at 5% and 10 m/s

The comparison of the wear maps for 20 mm particles at a concentration of 5% and a conveying velocity of 10 m/s, considering and neglecting turbulent modulation, as is shown in Figure 17, reveals noticeable differences in the wear patterns of the first elbow. When considering turbulent modulation, the region with the maximum wear in the first elbow is more concentrated compared to the case where turbulent modulation is ignored. However, the wear pattern in the second elbow remains almost the same in both cases.

For the condition of 20 mm particles at a concentration of 5% and a conveying velocity of 10 m/s, as is shown in Figure 18 and Figure 19, the simulation results considering and neglecting turbulent modulation follow patterns similar to the two conditions mentioned above. These patterns are not reiterated here. Considering the importance of the particle trajectories, particle distribution, and influence of particles on turbulent modulation throughout the entire elbow, the turbulent modulation source term is considered in the subsequent numerical simulations. Because the influence of particles on the turbulence dissipation rate is roughly similar to that of particles on turbulence kinetic energy, the following sections only introduce the influence of particles on turbulence kinetic energy under different conditions.

4.2. The Influence of Different Particle Diameters

This section explores how particle size affects wear rates, wear patterns, and turbulent kinetic energy in the series elbow at a conveying velocity of 5 m/s and a concentration of 5%. As is shown in Figure 20, we can see that when conveying 10 mm particles, the first elbow experiences maximum wear at approximately 85° near the elbow outlet, with a distinct strip-like distribution. When conveying 15 mm particles, the maximum wear position of the first elbow occurs between 75° and 90°, with a wear pattern that combines strip-like and point-like distributions. With 20 mm particles, the maximum wear position of the first elbow is between 60° and 90°, and the wear pattern is notably different from that of 10 mm particles, resembling a point-like distribution. The variations in wear patterns are primarily due to the different Stokes numbers of the three particle sizes.

The Stokes number is the ratio of the particle relaxation time to the characteristic time of the fluid and reflects the magnitude of particle inertial effects relative to drag. It is a dimensionless number used to characterize the inertia of particle motion. A smaller Stokes number indicates less particle inertia, making it easier for particles to follow fluid motion. Conversely, a larger Stokes number means less drag force on the particles, resulting in less momentum exchange between particles and the fluid, making it more difficult for particles to follow fluid motion, as shown in Figure 21, the formula is defined as $St=\frac{{\rho }_p{d}_p^{2}u}{18\mu D}$, where D is the characteristic length (hydraulic diameter in this case). Under similar conditions, compared to 20 mm particles, 10 mm particles have smaller diameters, and the fluid forces they experience are relatively larger. As a result, 10 mm particles can better follow the streamline motion. Consequently, most particles gather near the outlet of the first elbow, creating a strip-like distribution, as shown in the figure. When the particle diameter is larger, particle inertia plays a dominant role. Twenty-millimeter particles are more likely to deviate from the streamlines, carrying more energy and directly colliding with the wall along their previous trajectories, resulting in a point-like wear distribution as shown for 20 mm particles.

As particles exit the first elbow, they continue to move within a horizontal pipe that is 3D in length. However, because the pipe is relatively short, the secondary flow from the first elbow does not have enough time to fully develop by the time it reaches the second elbow. Furthermore, due to a flow velocity of 5 m/s, the gravitational effect on the particles is not significant, and inertial effects still dominate. The particles tend to gather around 30° from the outlet of the second elbow until they exit. The primary wear occurs from around 30° from the outlet of the second elbow to the elbow outlet. When comparing the wear rates of the elbows for different particle sizes (10 mm, 15 mm, and 20 mm), it is observed that the wear rates initially increase and then decrease. Analyzing the collision frequency histogram, as shown in Figure 22, it is discovered that the wear rate does not exhibit a linear relationship with particle size or the collision frequency of particles with the elbow. This could be because larger particles, although with fewer particles at the same conveying concentration and a lower collision frequency with the wall, carry more energy. Thus, larger particles directly impact the outer surface of the elbow, causing more significant wear. Therefore, this study concludes that the wear rate of the elbow is not simply linearly correlated with the collision frequency of particles with the wall or the energy carried by individual particles.

As shown in Figure 23, the distribution of turbulent kinetic energy (TKE) within the first elbow when conveying particles of different diameters. It is observed that changes in particle diameter do not cause significant variation in the TKE within the first elbow.

When conveying particles of different diameters, the distribution of turbulent kinetic energy (TKE) within the second elbow is shown in the Figure 24. For the conveyance of 10 mm particles, the maximum turbulent kinetic energy is mostly near the outer wall surface of the section at 30° from the exit. With 15 mm particles, the maximum turbulent kinetic energy is at the outer wall surface of the section at 60° from the exit. The turbulent distribution for 20 mm particles is similar to that of 15 mm particles, but the turbulent kinetic energy on the outer wall surface and at the center of the pipe increases with particle diameter, consistent with Section 4.1 and the study by Hetsroni et al., and further details are not reiterated here.

4.3. Effect of Different Particle Concentrations

This section delves into the effects of conveying speed and concentration on the wear morphology of series-connected elbows and turbulent kinetic energy. The Figure 25 illustrate the impact of 20 mm particles at a conveying speed of 5 m/s but in different conveying concentrations. As shown in Figure 26, at 1% concentration, the first elbow experiences maximum wear between 70° and the exit, while at 10%, the maximum wear location is near 45° to 60°. At 5% and 10% concentrations, the wear morphology takes on a strip-like appearance. However, at 1% concentration, the wear morphology of the second elbow resembles a “Y” shape. This is because, with fewer particles, most of the conveyed particles collide with the outer wall surface near 30° from the exit upon entering the second elbow. Subsequently, they slide along the outer wall surface, and the influence of the secondary flow within the elbow widens the collision range with the outer wall surface near 30° from the exit, forming a “Y”-shaped wear morphology. As shown in Figure 27, with the increase in particle transport concentration, the collision frequency between particles and the wall surface significantly increases.

The figure below depicts the effects of 20 mm particles at a conveying speed of 5 m/s, varying in conveying concentrations, on turbulent kinetic energy. As shown in Figure 28 and Figure 29, the turbulent kinetic energy increases in both elbows as the particle conveying concentration rises. This trend is especially apparent on the outer wall surfaces of the elbows and in their central regions. These observations are consistent with the results presented in Section 4.1 and corroborated by Rashidi’s research. Further elaboration on this topic is not reiterated here.

4.4. Effect of Different Conveying Speeds

The distribution of 20 mm particles at a conveying concentration of 5% but different conveying speeds in the pipeline, their impact on the wear patterns of the elbows are shown in the Figure 30 and Figure 31. It can be observed that, compared to a conveying speed of 5 m/s, at a conveying speed of 10 m/s, the trajectories of particles are closer to the outer wall surface of the elbows. As analyzed earlier regarding the Stokes number of particles, under other identical conditions, particles at higher flow speeds have a larger Stokes number, where the dominant effect is inertia, making them less likely to follow fluid motion.

Comparing the wear patterns of particles on the elbows at different conveying speeds, it is evident that particles cause greater wear on the elbow wall surface at higher speeds. Compared to 5 m/s, at a conveying speed of 10 m/s, particles concentrate and cause severe wear on the first elbow near the inlet at 60° to 70°. This can be attributed to particles carrying greater energy and causing more direct impacts in this region. Comparing the wear on the second elbow, particles at a conveying speed of 10 m/s also cause more severe wear near the inlet compared to a conveying speed of 5 m/s.

The Figure 32 and Figure 33 shows the impact of 20 mm particles at a conveying concentration of 5% but with different conveying speeds on turbulence. Turbulent kinetic energy increases with the flow velocity.

4.5. Effect of Series Elbows with Different Structures

As shown in Figure 34, the wear patterns of 20 mm particles at a conveying concentration of 5% and a conveying speed of 5 m/s are presented for three different series-connected elbows. The first elbow is consistently arranged vertically upward to horizontally, and its wear characteristics with changing conditions have been analyzed in the preceding sections, and will not be reiterated here.

Based on the wear maps, it is apparent that the second elbow experiences the most wear in the range of 0–30° from the exit in all three structures. Further examination of the collision frequency between particles and the respective wall surfaces indicates that, as shown in Figure 35, when all other conditions are constant, the wear rate is closely linearly linked to the collision frequency.

In case 2, the orientation of the elbow changes from vertical to horizontal. When the particles flow out of the first elbow, their inertia causes wear to occur in the upper part of the second elbow. As the particles settle during the flow process, a significant area of wear appears on the second elbow, extending from top to bottom at an inclined angle.

In case 3, when the elbow moves from a vertically upward position to a horizontal position and then back to a vertically upward position, the flow direction of particles in the second elbow becomes opposite to the direction of gravity. This results in a lower wear rate of the second elbow as compared to the first two arrangements. However, the position with the highest wear rate is still found near the exit of the elbow.

The turbulent kinetic energy distribution inside the second bend of the three different series-connected bend structures is shown in the Figure 36.

5. Conclusions

This study utilized an advanced CFD-DEM numerical model to simulate the flow of solid–liquid two-phase large particles in a series of connected elbows and their impact on bent surfaces. The accuracy of the model was verified by conducting numerical simulations of the two-phase flow of large particles in horizontal and vertical pipes and elbows, and the simulation outcomes are consistent with experimental results. Additionally, the study examined the characteristics of solid–liquid two-phase flow in a series elbow, as well as the impact of large particles on the turbulent kinetic energy and turbulent dissipation rate inside the series elbow. Based on these findings, we can draw the following conclusions:

In this study, the turbulence kinetic energy and turbulence dissipation rate, when neglecting the modulation effect of particles on turbulence, were significantly lower compared to the results obtained when considering the modulation effect. Moreover, there is a noticeable difference in the wall wear rate between the simulations without considering turbulence modulation and those considering turbulence modulation. The modulation effect of centimeter-sized large particles on turbulence cannot be ignored.

When particles of varying diameters were conveyed at 5 m/s and 5% concentration, a study on the wear rates of elbow components showed various patterns and the most severely worn positions were also different. Initially, the wear rate increased and then decreased for 10 mm, 15 mm, and 20 mm particles. This suggests that the correlation between the wear rate and particle collision frequency and energy is not linear. Additionally, the study revealed that larger particles had a more pronounced effect on the turbulence kinetic energy and turbulence dissipation rate in the elbow and pipe center.

When conveying 20 mm particles at a velocity of 5 m/s, the wear on the curved pipe’s extrados increases as the particle concentration increases. Particle concentration is a significant factor that causes erosion in the bend. As particle concentrations increase, the modulation effect on turbulent kinetic energy becomes more pronounced, emphasizing the importance of considering the modulation effect of large particles on turbulence when dealing with high conveying concentrations. However, at low concentrations of only 1%, the wear pattern of large particles on the second bend exhibits a Y-shaped distribution instead of a belt-like distribution typically seen in moderate to high concentrations.

At a conveying concentration of 5%, 20 mm particles moving at high conveying velocities tend to have higher Stokes numbers. This means that the dominant force acting on them is inertia, which makes them less likely to follow fluid motion and more likely to directly impact the first elbow. As a result, the location of maximum wear on the first elbow changes, and the wear pattern becomes more concentrated. By observing the particle distribution in the pipeline, it can be seen that at higher conveying velocities, particles tend to move closer to the outer wall of the pipe.

The wear pattern and wear rate of the second elbow are affected by its arrangement, particularly the direction of particle flow. In case 3, where the flow direction is opposite to gravity, the collision frequency of particles with curved pipes is the lowest. As a result, the wear rate of the second elbow in case 3 is minimal compared to the other two arrangements. It is recommended to use series elbows similar to case 3 in engineering to reduce wear. Additionally, it has been observed that the wear amount of structurally connected elbows is linearly correlated with the frequency of particle collisions, provided the particle size, conveying velocity, and conveying concentration remain consistent.