2.1 Computational geometry and mesh generation

A 3-D geometry was created in ANSYS DesignModeler using the same computational geometry dimensions of an elbow pipe used in the study conducted by Wang et al. [2]. The model geometry is a 90° standard elbow pipe (R/D = 1.5) with an internal diameter of 40 mm and a curvature radius of 60 mm. In order to achieve a fully developed flow, the lengths of 400 mm (10D) for both horizontal pipe upstream and vertical pipe downstream of the elbow were used. The general flow conditions are shown in Table 1.

The grid was defined to provide a reduced solving time and increase the accuracy of the numerical solution. A finer grid scheme was created near the pipe wall with five (5) inflation-layers to obtain an accurate flow field near the pipe wall and provide a better erosion prediction [16]. The entire computational domain was discretised using hexahedral grid type. Hexahedral mesh has the advantage of providing high-quality solution with less number of computational cells and it is generally more aligned with the flow. The total number of nodes in the model was 108 315, with 306 109 faces and 103 912 cells.

2.2 Flow field modelling

The liquid is treated as a continuous phase and modelled by solving the Navier-Stokes conservation of mass and momentum equations for fully developed incompressible pipe flow. Taking into account the volume fraction of the liquid phases, the conservation equations are expressed as follows [1, 17].

Where α is the volume fraction, ρ is the density, u⃗ is the velocity vector, P is the static pressure, τ= is the stress tensor, and ρg⃗ is the gravitational body force. The subscript f represents the liquid phase: water (W) and oil (O).

The stress tensor τ= is calculated as follows:

Where μ is the molecular fluid viscosity, I is the unit tensor, and the second term on the right hand side is the effect of volume dilation.

The realizable k-ε two-equation model is adopted to describe the kinetic energy encountered in the continuous phase and dissipated in the dispersed phase. And the standard wall function approach is used to model the flow in near-wall region. The modelled transport equations for k and ε in the realizable k-ε model are as follows:

Turbulent kinetic energy transport equation:

Dissipation of turbulent kinetic energy transport equation:

The turbulent eddy viscosity is computed from:

In the above equations, G_k represents the generation of turbulence kinetic energy due to the mean velocity gradients, C₂ and C₁ are constants, σk and σε are the turbulent Prandtl numbers for k and ε, Sk and Sε are the user-defined source terms.Turbulence model coefficients are given in the Table 2.

2.3 Particle tracking

The motion of the particles in discrete phase is described by the Newton's laws of motion which is the sum of all fluid-solid interaction forces. For the purpose of this study, the motion of a particle in fluid is mostly affected by the drag force. The governing equation of particle motion proposed by Newton's second law is:

In the expression, the first term on the right hand side represents the drag force, FD→, caused by the velocity difference of the fluid and the particles, and the other terms Fg→, FP→ and FVM→ represent the gravity force, the pressure gradient force and the Virtual mass force, respectively.

The drag force is the most dominant term computed by the following relationship:

The drag coefficient,CD, in this study is expressed based on spherical drag law for smooth sand particles as the following equation:

Where a1, a2 and a3 are constants that are applied over a wide range of Re given by Morsi and Alexander [18].

Re is the particle Reynolds number calculated from the following relationship:

Where up is the particle velocity, dp is the particle diameter, ρp is the density of the particles, ρf and μ are the density and viscosity of the fluid, respectively.

The gravity force is expressed as:

Which is computed by the difference in sand particle density (ρP) and the density of the fluid (ρf)

The pressure gradient is defined as:

The virtual mass force can be expressed as:

2.4 Particle-wall collision and erosion model

During particle-wall collision, particles hit the wall surface of the material and then rebound back to the fluid domain. Due to the effects of the impact, particles lose some kinetic energy; the rebound velocity is always lower than the incident velocity. This impact is characterized by the change in particle momentum through the coefficients of restitution, en and et for the normal and tangential velocity components, respectively [16, 19]. en and et are written as:

Where upn and upt are the normal and tangential velocity components of the particle, respectively. Subscript 1 refers to the condition prior collision and subscript 2 is after collision.

In this study, the restitution coefficients in the normal and tangential directions to the wall are the default correlations by Ansys Fluent, as follows:

Where α is the particle impact angle.

When the particles impact the wall material, the impact information such as location, velocity and angle of impingement is recorded in each CFD cell next to the pipe wall. These data are then introduced to the selected erosion equation monitored at the wall boundaries. The erosion equation model used in this study is the default equation provided by Fluent, defined as follows:

Where:

1. Cdp : particle mass flow rate

2. α : impact angle of particle path with the wall face

3. fα : function of impact angle

4. v : relative particle velocity

5. bv : function of relative particle velocity

6. Aface : area of cell face at the wall subject to erosion

C, f and b are the constant functions defined as part of the wall boundary conditions and their numerical values based on the pipe material used. The values of the functions C, f and b are 1.8E-09, 1, and 0, respectively, as specified in Ansys Fluent guide.

A piecewise-linear profile in Fluent is used to define the impact angle function found in the literature, as given in Table 3. The diameter function is defined at the value of 1.8E-9 and the velocity exponent is set at the constant value of 2.6 because it is consistent with the value in the open literature for sand.

5.1 Model verification

A verification case is carried out to validate the numerical model presented in this study under the same experimental conditions analysed by Blanchard et al. [7]. They conducted this experimental analysis of erosion process of a 90-degree pipe elbow with various sand particle sizes ranging from 95 μm to 605 μm in diameter. They found out that the maximum point, at which erosion occurs, is located to be at an average angle of 85-degree, measured from the starting point of the curve. Furthermore, the results indicated that for an elbow whose R/D ratio is 1.5, a consistency was observed in the location point of the maximum wear.

The verification case for erosion prediction involves two phase flow particle-water, simulated using CFD code in a 90-degree pipe elbow with R/D ratio equal to 1.5. The pipe diameter is 25.4 mm, and the curvature radius is 38.1 mm. The fluid with an inlet velocity of 5 m/s, enters from the vertical downstream pipe and discharges horizontally with lengths straight pipe of 10D each (D is the inner diameter of the pipe). The solid particles of 605 μm in size with density of 2500 kg/m³, are injected uniformly at the pipe inlet with the same velocity as the fluid and a mass flow rate of 0.1 kg/s [2, 7]. The numerical solution was obtained after 200 iterations.

The comparison of the erosion profile between the predicted numerical results and the experimental data along the wall of the elbow can be seen in Figure 1. It can certainly be shown that there is a similar tendency between the predicted results and the experimental data. In addition, the maximum wear location occurs at the exit of the elbow (pick point of the graphs) at an average angle ϕ = 86°, measured from the starting point of the bend wall. Comparably, the maximum erosion depth occurs at the average angle ϕ = 85° of the bend wall in the experiment. However, discrepancies can be observed (under prediction) between CFD predictions and actual measured values of erosion depths. This may be due to the dissimilarities in wear coefficient and mostly the velocity exponents between Finnie erosion equation used for the experiment and the erosion equation provided by Fluent used in the CFD modelling at the wall boundary. Additionally, the authors indicated that the 2-D experimental model failed to account for the secondary flow effects which were present, whereas the present study in 3-D successfully predicted these flows. Therefore, these reasons are suspected to be the causes of the lower CFD erosion depths prediction by two orders of magnitude compared with the experimental erosion depths.

Figure 1:

Comparison of erosion depth profile between simulation results with experimental data of Blanchard et al. [7].

5.2 Flow analysis of oil-water in elbow pipe

Figure 2(a-d) below shows the contours of volume fraction of oil and water at different oil viscosity conditions. It can commonly be observed in each case that the horizontal upstream section of the pipe presents a constant green colour for the water volume fraction and a constant yellow colour for oil volume fraction. Since the water volume fraction in that section is about 65 %, thus the remaining 35 % of the flow is oil when we are looking at the contours. Therefore, oil and water are fully mixed. The mixed flow is developed before and after the bend with insignificant change in volume fraction of each phase. In addition, Figure 3 illustrates the velocity contours of oil and water which indicate that the phases are flowing at the same velocity. However, oil and water are two immiscible fluids due to the difference in physical properties. The flow may be suspected to form emulsions favour by the high velocity, in which tiny oil droplets could be dispersed in water since in this case water occupies 65 % of the total volume in the pipe. Therefore, this result comes into agreement with Nädler and Mewes [22], that the flow characteristic of oil and water in pipes highly depends on the volume fraction and droplet distribution of the dispersed phase because of the density difference between the two phase.

Figure 2:

Contours of volume fraction of oil and water at different oil viscosity conditions in symmetric plan: (a) 0.0572 Pa.s, (b) 0.0911 Pa.s, (c) 0.155 Pa.s, (d) 0.287 Pa.s.

Figure 3:

Velocity distribution of oil and water in symmetric plan at different oil viscosity conditions: (a) 0.0572 Pa.s, (b) 0.0911 Pa.s, (c) 0.155 Pa.s, and (d) 0.287 Pa.s.

Furthermore, due to the effects of oil viscosity increase in the model, the shear stress (τ) is developed by the friction between the fluids particles and consequently induce pressure drop in pipe. Because of the viscous stress created along the pipe, according to Newton's law of viscosity, the fluid (oil-water) molecules interact and exchange momentum with the discrete phase sand particle by the influence of viscosity and velocity gradient, through the linear relationship in eq. (20). Additionally, high oil viscosity might produce the presence of oil film with the tendency to stick on the pipe wall, thus reduce the kinetic energy of the solid particle to collide on the wall bend [23].

Where, μ is the dynamic viscosity (coefficient of viscosity), ux is the flow velocity along the boundary, and y is the height above the boundary.

5.3 Velocity distribution along the elbow

The results of velocity distribution for different oil viscosity conditions are shown in Figure 3(a-d) as symmetry plane. The magnitude of the velocity of the continuous phase is clearly indicated by the colour. It can be observed that the maximum velocity zone is developed at the inner wall (intrados) of the bend and a lower velocity zone is developed at the outer wall (extrados), with a separation of a clear boundary layer. The maximum velocities at the intrados are 40.23 m/s, 40.43 m/s, 40.76 m/s and 41.01 m/s for (a), (b), (c) and (c), respectively and tend to slightly increase with the oil viscosity. This numerical profile is in qualitative good agreement with the experimental results by Enayet et al. [24] for flow behaviour in a pipe bend. The high velocity near the intrados of the bend can be initially explained by the sudden change in flow direction; hence the flow becomes more turbulent in that zone than in the horizontal upstream and vertical downstream zones. In the circular cross-section the pressure distribution becomes uneven and the high velocity is directing the flow from intrados to extrados of the pipe due to the action of centrifugal force. Additionally, the presence of secondary flow is expected in that region of the pipe, as a result of high velocity pattern close to the intrados of the bend [25].

5.4 Pressure distribution along the elbow

The static pressure distribution along the pipe domain can be shown in Figure 4(a-d) for different viscosity conditions. The contours show that the maximum pressure at the inlet of the pipe elbow is generally increased with the increase in oil viscosity. The common trend observed in each case is that the pressure drops gradually along the flow direction. At the bend section the direction of the fluid is changed and this in manifested by the migration of the fluid and solid particle from the inner to the outer side of the bend. Thus, the motion of the fluid is under the effect of centrifugal force as a result of inertia. This makes the static pressure on the extrados higher than that near the intrados, as is indicated in Figure 4 [26]. The local static pressures located at the extrados of the bend are 3.95E + 05 Pa, 4.14E + 05 Pa, 5.10E + 06 Pa, and 6.52E + 05 Pa with the respective increase of oil viscosity. In addition, when the fluid exits the bend the static pressure along the vertical downstream straight pipe will not restore its original state due to inertia effect. Therefore, the colour in the contours indicates that the static pressure is gradually decreasing along the flow direction as a result of friction [27].

Figure 4:

Pressure distribution in the elbow for different oil viscosity conditions: (a) 0.0572 Pa.s, (b) 0.0911 Pa.s, (c) 0.155 Pa.s, and (d) 0.287 Pa.s.

5.5 Effect of oil viscosity on erosion rates

The increase in oil viscosity clearly demonstrates that it has a significant influence on the erosion rates of solid material. It can be shown in Figure 5 below that, erosion rate decreases with the increase in oil viscosity. It is worth noticing that when oil viscosity is increased from 0.0572 Pa.s to 0.09106 Pa.s, the maximum erosion rate on the elbow is slightly decreased from 3.24E-05 kg/m².s to 3.035E-05 kg/m².s and then drastically decreases to 2.154E-05 kg/m².s when the viscosity is increase to 0.155 Pa.s. This result indicates that the viscous fluid reduces the potentiality of the sand particle to strike on the pipe wall, resulting in a decrease of material failure [28]. Furthermore, the motion of the sand particle is given by the surrounding liquid phase through the drag force due to their relative velocity difference. This may imply that, as the fluid becomes more viscous it restricts the motion of the particles and reduces particle kinetic energy, hence their ability to erode [29]. In addition, a presence of liquid film may be expected to form at the wall surface of the elbow with a finer grid generation. As the carrier liquid mixture becomes viscous, the liquid film also becomes thicker. As a result, the liquid film slows the particle's impact velocity [30].

Figure 5:

Effect of oil viscosity on erosion rate.

5.6 Distribution of erosion rate and location of maximum wear point

Erosion rate distributions on the walls of the elbow are presented as contours in Figure 6. The contours clearly show that, erosion is distributed on the extrados and on side walls of the downstream straight pipe near the bend exit for all simulations. Erosion is negligible on the upstream and downstream straight pipes, since the solid particles are flowing along the fluid streamlines parallel to the pipe axis, whereas the bend wall is more vulnerable due to the change in flow direction. In addition, the erosion rates per unit area reach their maximum values of 3.24E-05 kg/m².s, 3.035E-05 kg/m².s, 2.154E-05 kg/m².s and 1.895E-05 kg/m².s, respectively at the extrados near the bend exit, approximately at an average angle of 87-degree measured from starting point of the bend, as further illustrated in Figure 6. Therefore, the maximum erosive locations occur at an average angle of 87-degree, which is in reasonable agreement with the experimental data of Blanchard et al. [7].

Figure 6:

Maximum erosion location on the elbow.

A sample particle trajectory is displayed in Figure 7 in order to analyse the motion of the particle in the flow field for the erosion phenomenon. The trajectories of the particles are mostly governed by the continuous phase through the drag force. When approaching the bend, the flow is fully developed, and the solid particles have sufficient momentum and also the sudden change in flow direction causes them to deviate from the fluid streamlines. Therefore, the solid particles impinge directly on the extrados near the bend exit and cause severe erosion in location A as illustrated in Figure 7. Moreover, the secondary flows due to centrifugal effects inside the bend section drive a number of particles in a circular motion. Hence, the effects of secondary flow convey the particles to the side walls of the downstream straight pipe near the exit and as a result, causing erosion in location B [2].

Figure 7:

Sample particle velocity trajectories in elbow.

5.7 Effect of secondary flow in elbow

The three-dimensional simulations successfully predicted the presence of strong secondary flow effects at the exit section of the elbow. Figure 8(a-d) illustrates the secondary flow paths by the velocity vectors at the exit section of the elbow for respective oil viscosity conditions. As discussed earlier on, when the fluid flow through the elbow, the velocity field near the intrados is greater than that near the extrados. The circumferential motion of the fluid at the intrados as shown in Figure 8 is due to the strong effect of centrifugal forces as compared to the extrados.

Figure 8:

Velocity vectors in the cross-section exit of the bend for different oil viscosity conditions: (a) 0.0572 Pa.s, (b) 0.0911 Pa.s, (c) 0.155 Pa.s, and (d) 0.287 Pa.s.

In addition, the pressure gradient at the extrados is larger than that at the intrados. As an outcome, the flow is directed back from the extrados region to the intrados region near the pipe walls, thus the secondary flow occurs. Therefore, the secondary flow vortices direct a number of solid particles by the surrounding fluid from the intrados region to the extrados region [27]. As a result, severe erosion is observed on the side walls of the downstream straight pipe near the bend exit. According to the velocity vectors, it can clearly be confirmed that the fully developed secondary flows force the particle to impinge on the side wall of a 90-degree bend pipe. Based on this result, the distribution of erosion on the side walls near the bend exit can be justified.