Numerical Investigation of Turbulence Production Under Two-Phase Flow in a Simulation of Electrical Submersible Pump Performance

Centrifugal pumps applications range from domestic applications to large industrial applications. Certain industrial applications, such as pressurized water reactors, require liquid-gas mixture as their primary working fluid; however, these pumps are impractical for two-phase flow. On the other hand, an electrical submersible pump (ESP) is an extended centrifugal pump takes liquid directly below the discharge point. ESPs are primarily needed in the petroleum industry to boost production from wells that contain gas residue in an oil-gas mixture flow [1]. ESPs are positive displacement pumps installed downhole to transfer kinetic energy to the fluid via a downhole motor driving a multistage centrifugal pump [2]. Relevant studies can be found in Refs. [3,4].

The formation of gases in the process is destructive to the pump resulting in the loss of pump performance, which in turn causes safety issues, economic losses, and significant damage to the pump [1]. Higher gas formation near the impeller blades blocks the flow channels and opposes fluid flow [5]; moreover, gaseous flow conditions and pressure surging significantly reduce the head. ESPs currently on the market are capable of handling two-phase flows effectively. However, the exact mechanisms that lead to deterioration of the head under gaseous flow and surging initiation are still unclear [6]. Variables such as bubble coalescence and breakup mechanisms are difficult to visualize due to the complex geometry of multistage ESPs, as illustrated in Fig. 1.

According to Kundu et al. [7], turbulence involves fluctuations that are unpredictable and difficult to capture by statistical analysis. A turbulent flow instantaneously satisfies the Navier–Stokes equation, and it is practically impossible to predict turbulent flows in detail due to the large length and time scales involved, especially at higher Reynolds numbers. Situations where the duration of the flow exceeds the time scale of the turbulent-fluctuation are not uncommon. The tensor $−ρ0ui¯uj¯$⁠, which relates to the Reynolds stress, shows that the Reynolds stresses are usually larger than viscous stresses. However, near a solid surface, the fluctuations approach zero and the mean flow gradients reach an all-time high.

According to Dol et al. [8] and Kundu et al. [7], a particle starts at y and is displaced with a positive vertical velocity by a positive velocity change v upward toward y + dy. Assuming that the fluid particle moves upward due to a slight velocity change (v > 0), it tends to keep its horizontal velocity unchanged during the displacement. However, once it reaches y + dy, the horizontal velocity acting on the particle increases, resulting in a slower horizontal velocity (u < 0) compared to the neighboring particles. On the other hand, fluid particles moving downward due to fluctuations (v < 0) experience a positive horizontal velocity (u > 0) when they reach y − dy. In summary, a positive vertical velocity v correlates with a negative horizontal velocity u and a negative vertical velocity v correlates with a positive horizontal velocity u always with a negative $uv¯$⁠, where according to Fig. 2, du/dy > 0. The Reynolds stress terms are obtained from the nonlinear advection term U_j(∂u_i/∂x_j) in the momentum equation and the average stress due to turbulence at mean flow. Figure 2 illustrates the evolution of the nonzero Reynolds shear stress in a simple shear flow. Another view of the Reynolds stress involves the rate of mean momentum transfer due to turbulent fluctuations

Current data on wall shear stress in two-phase flows involving bubbly flow are inadequate, although they have been applied extensively in recent technology, including nuclear reactors, chemical and petroleum industries, and heat exchangers for electronic cooling. The data are illuminated by realizing the bubble behavior in the impeller channels as shown in Fig. 3. The wall region of the impeller in a detailed controlled flow environment has not been extensively observed and requires further investigation [9].

Reynolds shear stress in bubbly flow contributes to the frictional pressure drop between the liquid and gas phases, which increases turbulence in the impeller channels and decreases pump performance as more gas flows through the liquid. The study of Reynolds shear stress along with turbulence in a bubbly flow has not been thoroughly investigated.

This research work investigates the effect of turbulence on pump performance in two-phase flow by numerical simulation in ANSYS of a single-stage ESP. GDIWA ISO 9906-Annex ESP provides a rational perspective of the internal flow profile for a descriptive analysis of flow structures. The numerical simulation considers the two-phase parameters, including gas volume fraction (GVF) and impeller rotational speeds [10]. The boundary conditions include varying gas volume fractions associated with the liquid flowrate, which applies the mixing model (simplified Eulerian model).

This research work creates a comprehensive and descriptive application of fluid behavior using a comprehensive understanding of the fundamental structures of incompressible fluid flows. This research work is intended to accomplish the following objectives:

1. Simulation of two-phase flow in ANSYS Fluent under 1%, 6%, and 10% inlet gas volume fraction (IGVF).

2. Observation of turbulence characteristics through vortices, pressure, and turbulent kinetic energy contours between impeller channels.

3. Investigation of turbulent flow while maintaining the concepts of Reynolds stress and wall shear stresses at a fixed RPM of 1500.

Current research in turbulence addresses two discrete aspects of fluid flow, namely, overlooking near-wall turbulent structures and skin friction on the wall; both aspects are combined to detect near-wall turbulent structures [13]. The turbulent flow takes into account an additional diffusion mechanism that serves to transport the fluid momentum toward the boundary. Above the mean shear, which is parallel to the base, smaller macroscopic eddies are present that tend to balance the velocity distribution of the momentum near the boundary layer. According to Newton’s second law, the rate of momentum transport through the turbulent motion is identical to the shear stress in the entire plane. This is called turbulent shear stress, or Reynolds stress and has the same influence as the frictional force, but between two fluid layers and not between fluid and solid [14].

The ESP model used in the numerical simulation is replicated with attention to detail, as opposed to a detailed geometric model of the GDIWA A series ESP. Pump replication suitable for numerical simulation starts with an accurate and precise model of the pumps, including exclusive components, as shown in Fig. 4. Then, ANSYS forms a fluid zone to proceed with the computation. Water is chosen as the typical fluid material, and the developed fluid zone is shown in Fig. 5.

Table 1 illustrates the impeller geometry, which is the most challenging pump geometry due to the dimensions used to design the impeller.

The mesh coverage of the fluid zone is based on the fluid flow through the zone, and the mesh elements must interlock between the interfaces. The mesh encloses the entire fluid zone with tetrahedral elements, the most appropriate element shape used. The refined mesh for the fluid zone is shown in Fig. 6, where the vanes and diffuser represent the smallest element size. The mesh is refined according to the complex geometry of the impeller blades and diffuser, while the volute and inlet are covered by a less refined mesh due to their simple shape.

The initial mesh model is simulated as a steady-state fluent model with single-phase fluid flow and is gradually scaled up to find the best-fitting model, as shown in Fig. 7. In this way, the optimal mesh coverage can be selected for the available computing power. Figure 7 also shows the simulation results in accordance with different mesh element sizes for a comprehensive comparison. For a network coverage of 3 million elements, the simulation consumed about 20 min per simulation with a rational uncertainty, which is reasonable given the available computing power and is selected for further simulations. A mesh coverage of more than 3 million elements would not be feasible in terms of the tradeoff with error reduction. The Y+ mesh value is less than 1 for the shear stress transport (SST) turbulence model due to the lack of a wall function, and the need to investigate areas of low Reynolds number and areas contributing to pressure drop.

A key difference between the k − ω and the k − ɛ SST turbulence models is the calculation of turbulence viscosity, which accounts for turbulent shear stress transport. The shear stress turbulence model contributes to a constitution term, shown in Eqs. (6) and (7), which allows accurate calculations of the near-wall shear stresses. The software built-in functions trigger the k − ω turbulence model when the flow reaches a near-wall region and trigger the k − ɛ model in regions away from surfaces. Compared to the traditional k − ω and k − ɛ models [17], this allows mimicking accurate results for a range of different flow structures. Further explanation can be found in Ref. [16]

The impeller rotation speed is set at 1500 rpm as this provides a good balance between rotation speed and GVF fluctuations. The liquid flowrate is kept constant at 250 L/min, which is the Best Efficiency Point (BEP) in the initial single-phase tests of ESP. The GVF (λ_g) is varied by 1%, 6%, and 10%.

The fluid in the pump is liquid water, while the gas is air. The transition occurs at the inlet, followed by the upper volute, which forms an interface between the upper volute and the impeller, leading to a lower volute that exits at the diffuser. Second-order coupled equations are used for the simulation, ignoring all relaxation terms and adding the patch to all components of the pump. The rotational speed of the impeller is assigned in terms of the flowrate, which is assigned as the velocity inlet along with other components of the pump, such as the rear shroud, impeller vanes, and impeller hub. The outlet of the pump is assigned as the pressure outlet.

Gas–liquid flow is a complicated shear stress mechanism with difficult correlations. Gas–liquid mixture involves a complicated shear stress mechanism with problematic calculations. Shear stress models are based on both assumptions and empirical closure equations. Moreover, turbulence models aim to represent the Reynolds stresses TRANS according to Reynolds averaging. Turbulence models are characterized as first-order models that refer to a single equation (Spalart–Allmaras) or a dual equation model (k − ɛ and k − ω) [18].

The compiled shear stress due to the shear plane in the flow sums the turbulent shear stress that occurs due to the macroscopic diffusion of the fluid momentum along with viscous shear stresses (molecular diffusion of the fluid momentum along with attractive molecular forces). The macroscopic nature associated with turbulent shear stress can be assigned to a point of shear stress, even if it has a real physical meaning from one point to another. It is derived as an average over the area of the shear plane, since the turbulent flow is associated with both the magnitude and orientation of the shear, which varies (continuously and on the scale of eddies) from point to point [14].

Since the turbulent shear stresses are larger than the viscous shear stresses (except near the boundary), the differences in time-averaged velocity from one layer to another in a turbulent flow structure spread over a large portion of the flow depth compared to a laminar flow structure. This spread is associated with a simple velocity gradient du/dy over most of the flow depth in a turbulent flow compared to a laminar flow; however, this results in a higher velocity gradient profile near the lower boundary where viscous effects dominate. This occurs because the shear generated by the transition from a still large velocity near the boundary to a velocity of zero at the boundary is compressed in a thin layer that immediately alternates with the boundary [14].

After several simulations for the k − ω SST for single-phase flow, the results show a reasonable comparison within the uncertainty limits, as shown in Fig. 8. Shojaeefard et al. [19] performed computational fluid dynamics (CFD) simulations using the k − ɛ turbulence SST model to investigate the effect of viscous fluids on the centrifugal pump. The use of the k − ω SST model ensures the use of both the k − ω and k − ɛ turbulence models [20]. The percentage deviation of the simulated pressure head at the flowrates of interest is compared to the actual pressure head results in a minimal deviation $(μ<5%)$⁠, which confirms the reliability of the model.

In comparison with the ESP simulation performed by Zhu and Zhang [21], this research model is verified by comparing the pressure contours, as can be seen in Fig. 9, which shows a similar pressure distribution at 1500 rpm with a slight difference of 1% in the GVF.

Figure 10(a) shows the pressure contour on the impeller blades for 1500 rpm at 1% GVF and 250 L/min. The contour represents a straightforward pressure distribution along the pump, with pressures reaching up to 1.494 × 10⁴ Pa on the suction side of the impeller and dropping to 1.458 × 10⁴ Pa on the pressure side, reaching maximum pressure of 3.603 × 10⁴ Pa at the impeller tips. Figure 10(b) shows the contours of static pressure at the impeller blades with the highest peak pressure of 3.134 × 10⁴ Pa at the impeller tips due to increased wall shear stresses resulting from interaction with gas and liquid molecules. Pressures up to 7.217 × 10³ Pa prevail on the suction side of the impeller. Detailed information can be found in Ref. [10]. It is imperative to create a lower pressure or negative pressure at the suction side because the water is displaced from the suction side and this lower pressure helps to draw fresh water back into the system due to the pressure differential [22]. Figure 10(c) illustrates the pressure contours at the impeller blades when the GVF is increased to 10%. The maximum pressure at the impeller tips remains at 3.435 × 10⁴ Pa, while the pressure at the suction side of the impeller decreases to 2.250 × 10³ Pa. The maximum pressure at the impeller tips increases due to the higher pressure required to pump air bubbles along with the liquid, which consumes energy.

Figure 11(a) shows the pressure contour in plane view of the pump, which has high pressures of 4.036 × 10⁴ Pa at the edges of the volute casing and 1.052 × 10⁴ Pa on the suction side of the impeller. Figure 11(b) at 6% GVF shows the lower pressure on the suction side of the impeller with increasing fluid pressure as it exits the impeller channels. Prado and Pessoa [23] found that the pressure increase in the first stage tends to be positive at higher GVF, while it is negative at low GVF. However, further investigation is needed to determine whether the pressure drop occurs at the pump inlet or in the first stage where the mixing of water and gas occurs.

The impeller eye indicates a higher maximum pressure of 3.628 × 10⁴ Pa because it averages the pressure at the edges of the volute, which is not accounted for in the pressure contours of the impeller blades. The lower pressure on the suction side is 6.524 × 10² Pa. The view of the impeller suction shown in Fig. 11(c) shows a detailed picture at 10% GVF; however, the maximum pressure is 3.730 × 10⁴ Pa at the volute wall, while a pressure of 2.075 × 10³ Pa covers the inlet of the impeller channel closer to the pressure side of the impeller. There is higher pressure in the impeller channels than at 6% GVF because additional kinetic energy is consumed to generate the same pressure that occurs without the presence of gas.

The shear stress is a product of the dynamic viscosity μ and the shear rate. Since the mixture is incompressible, the dynamic viscosity remains constant, which means that the shear rate is directly determined by the wall shear stresses. Figure 12(a) below shows the wall shear stresses at the rear impeller shroud at 1% GVF. It shows the peak wall shear stress at the edge of the back shroud through which the fluid flows downstream toward the lower volute and eventually exits through the diffuser. The edge or rim of the shroud experiences a τ_w of 588.7 Pa since it is a steep slope. Shear stresses of up to 264.2 Pa occur in the impeller channels, illustrating the enhanced interactions between air bubbles and fluid particles. The shear stresses in Fig. 12(b) at the edge of the shroud increase to 534.4 Pa, while they decrease to 324.9 Pa on the blade pressure side, corresponding to a 3.39% decrease for a 5% increase in GVF. According to Fig. 12(c), the wall shear stresses at the shrouds increase to 500.6 Pa, while they are 303.38 Pa in the impeller channels, a decrease of 6.49% at 6% GVF.

Figure 13(a) shows shear stresses of 495.3 Pa at the impeller tips and much lower shear stresses of 126.3 Pa at the suction side of the blades. The shear stresses increase by 369 Pa as the mixture moves downstream through the impeller channels. The impeller blades in Fig. 13(b) have maximum shear stresses at the impeller tips ranging at 512.6 Pa, while the impeller vanes on the immediate suction and discharge sides have 96.24 Pa, showing a decrease of 0.775% with the increase in GVF. The wall shear stresses in Fig. 13(c) on the impeller blades show that the maximum shear stress at the impeller tips is 502.7 Pa, while on the impeller suction and pressure side experience 94.46 Pa.

Figure 14 shows a graphical approximation of the shear stresses on the path from the hub to the shroud edge. The reduction in shear stresses at the walls is the result of increasing bubble coalescence in the impeller channel, which is accompanied by a reduction in viscosity, causing the impeller to experience lower flow resistance (as seen in Figs. 13(b) and 13(c)), resulting in lower shear stresses. The total shear on the walls decreased by 9.41% from 1% to 10% gas volume fraction. Because of the direct relationship between wall shear stress and shear rate, a decrease in wall shear stress results in a decrease in shear rate and thus a smaller decrease in energy gain from turbulence. This reduction must favor an increase in pump pressure, which improves performance. However, the contribution of turbulence to the degradation of pump performance is minimal; other factors such as gas pockets, impeller inlet blockage due to bubble coalescence, pressure gradient, and bubble size have a much greater impact on performance [13].

The degree of circulation in centrifugal pumps can be assessed from the contours of the air volume fraction as well as the velocity vectors, which provides better insight into the circulation and vortices contributing to the shear rate and ultimately the turbulent kinetic energy. Figure 15(a) shows the volume fraction of air at 1%. The figure shows a significantly higher air volume fraction on the suction side of the impeller compared to the pressure side. In the impeller channels, bubble coalescence and interaction of the gas pockets with the fluid and the impeller walls occur. Figure 15(b) illustrates the air volume fraction at 6% GVF and shows the formation of large gas pockets on the pressure side, known as agglomerated bubble flow, with slight gas coverage at the impeller inlet.

The more bubbles form in the impeller channels, the smaller the empty spaces between the molecules become, leading to a stronger interaction between the molecules. Caridad et al. [24] have illustrated that the bubble diameter increases steadily with increasing no-slip GVF, leading to the formation of gas pockets along the impeller channels, which further degrade the gas-handling capability of an ESP in two-phase flow.

The enhanced interactions create much larger bubbles that disperse smaller bubbles that completely cover the impeller channels. The drag force is not sufficient to pull the bubbles downstream, so the bubbles circulate in the flow and create an unfavorable pressure gradient in the channels [1]. 1% and 6% GVF showed stronger gas pocket formation on the impeller suction side, while further increasing the GVF to 10% showed bubble coalescence at the blade pressure side. An experiment by Zhu et al. [25] showed that the bubbles accumulate near the leading edge and on the pressure side of the impeller blades (see Fig. 15(c)), which is explained by the detection of vortices at similar locations. Another observation comparable to this simulation was the recirculation around the trailing edges of the impeller blades, as shown in Fig. 15(b).

The figure clearly shows how the bubbles recirculate toward the impeller inlet. Figure 15(c) illustrates the volume fraction of air at 10% GVF. Larger bubbles with irregular shapes form at the impeller inlet; however, in this case, they merge with the bubbles in the impeller channels, suggesting that stationary bubbles reduce the kinetic energy transferred to the fluid so that larger bubbles do not decay [26]. Paternost et al. [27] agreed that the presence of larger bubbles at the pump inlet has an impact on pump performance, with bubble coalescence playing a crucial role. They further stated that large bubbles stall the impeller and lead to surge conditions.

Figure 16(a) shows the velocity vectors, which in conjunction with Fig. 15(a) illustrate the ability of turbulence to maintain its shape by absorbing the kinetic energy of the fluid. Figure 16(a) illustrates the key features such as downwash, circulation zones, and the mean flow directions. Figure 16(b) shows the vortex formation and circulation using velocity vectors. Compared to Fig. 16(a), Fig. 16(b) shows a highly organized gas–liquid flow with circulation extending along the pressure side of the impeller. The fluid has a velocity of 2.886 m/s at the center of the circulation, which gradually increases as it moves away from the center. Compared to 1% GVF, the velocity increased by 14.25% with a 5% increase in GVF, suggesting higher turbulence in the impeller channels. Figure 16(c) shows large circulating vortices in the impeller channel with a velocity of 2.898 m/s moving near the blade suction side. The figure shows a greater frequency of vortices compared to lower GVFs. The occurrence of recirculation at the edges of the impeller channels is relatively low. The recirculation constitutes negatively the turbulent kinetic energy and results in a lower pressure head at the gas passage. The velocity vectors agree well with those of Barrios et al. [28], with high velocities observed at the channel inlet where air generally flows in the circumferential direction from the suction side to the pressure side.

Observing the circulation at a sufficient shear rate to restore turbulent kinetic energy leads to an increase in energy consumption. Overall, Fig. (16) shows a pattern of increasing circulation in the impeller channels with increasing GVF. Higher vortices occupy the impeller inlet, restricting mass flow and consuming energy to maintain their turbulent shape. Most importantly, the rotational speed of the impeller is unable to eliminate the vortices when the gas enters the flow rapidly. The flow toward the impeller outlet forms a vortex in the liquid phase, which blocks the fluid at the suction surface, resulting in local vortices and backflows [29].

Figure 17(a) illustrates the turbulence intensity at different sections of the impeller channels. As expected, the turbulence intensity is slightly lower for two-phase flow than for single-phase flow, with maximum values of 2.819 m²/s² on the impeller suction side. The figure illustrates the extreme turbulent kinetic energy on the pressure sides of the blades as the fluids move through the impeller toward the volute. Figure 17(b) shows the turbulent kinetic energy at a higher GVF. It illustrates a highly ordered perspective on the turbulence with an intensity of 2.873 m²/s² at the impeller pressure side. Figure 17(b) correlates with Fig. 17(a), with the circulation extending along the pressure side and flowing downstream. The turbulence intensity at the impeller channel shows the vorticial positions contributing to the turbulent kinetic energy. Figure 17(c) shows an illustrative representation of the turbulence in the impeller channel, which reaches 2.765 m²/s², corresponding to a 3.76% reduction in turbulent kinetic energy. The clear trend shows a reduction in turbulent kinetic energy with an increase in GVF. As both terms in the equation for shear production diminish, the turbulence condition reduces, thus implying that turbulence consumes less energy as more and more gas occupies the impeller channels. As more gas pockets form and bubbles coalesce, the fluid forms less turbulence.

Regarding the reduction in the head, the pump head decreases dramatically by 26.53% when the GVF increases from 1% to 10%. The increasing coalescence of the bubbles and interaction with the fluid consumes a lot of energy and reduces the ability of the pumps to pump successfully. Larger bubbles tend to accumulate much more easily compared to bubbles with smaller diameters [30]. The size of the bubbles and the forces acting on the impeller are the two major causes of surging.

Current analysis of ESP multiphase flow (flow visualization and turbulence analysis) has helped identify the difficulties and complexities that lead to significant degradation of pump performance. Bubble coalescence and gas pockets play an important role in the head loss, as they can drain energy through turbulence in the impeller channels. The result shows that when the GVF increases from 1% to 10%, the shear rate due to wall shear stress decreases by 30.2%. Comparison of the different contours in the study illustrates that as GVF increases, the circulation begins to show a pattern of having an organized outline that contributes negatively to turbulence. Both terms reduce shear production, which in turn reduces the energy gain from turbulence. With increasing GVF in the impeller channels, the interaction between molecules increases along with the increase in turbulence. Compared to 1% GVF, the velocity increased by 14.25% when the GVF was increased by 5%, indicating higher turbulence in the impeller channels. The lower the GVF, the higher the ESP performance; lower kinetic energy consumption would result in higher pressure head and lower losses. The proposed flow structure as well as wall shear stress and circulation phenomena, such as vorticity and shear layers, can be further investigated in the context of two-phase flows and abrasive and erosive wear.

• This research was funded by Abu Dhabi Award for Research Excellence (AARE).

There are no conflicts of interest. This article does not include research in which human participants were involved. Informed consent is not applicable. This article does not include any research in which animal participants were involved.

The data sets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper. Data provided by a third party are listed in Acknowledgment.

m =

index for mixture

k =

Index for specific phase k

u =

velocity

S =

coefficient in turbulence model

α =

volume fraction

β =

coefficient in turbulence model

ρ =

density

τ =

viscous shear stress

τ_D =

diffusion stress

ω =

specific turbulent dissipation