A Review of Experiments and Modeling of Gas-Liquid Flow in Electrical Submersible Pumps

As the second most widely used artificial lift method in petroleum production (and first in produced amount), electrical submersible pump (ESP) maintains or increases flow rate by converting kinetic energy to hydraulic pressure of hydrocarbon fluids. To facilitate its optimal working conditions, an ESP has to be operated within a narrow application window. Issues like gas involvement, changing production rate and high oil viscosity, greatly impede ESP boosting pressure. Previous experimental studies showed that the presence of gas would cause ESP hydraulic head degradation. The flow behaviors inside ESPs under gassy conditions, such as pressure surging and gas pockets, further deteriorate ESP pressure boosting ability. Therefore, it is important to know what parameters govern the gas-liquid flow structure inside a rotating ESP and how it can be modeled. This paper presents a comprehensive review on the key factors that affect ESP performance under gassy flow conditions. Furthermore, the empirical and mechanistic models for predicting ESP pressure increment are discussed. The computational fluid dynamics (CFD)-based modeling approach for studying the multiphase flow in a rotating ESP is explained as well. The closure relationships that are critical to both mechanistic and numerical models are reviewed, which are helpful for further development of more accurate models for predicting ESP gas-liquid flow behaviors.


Introduction
The electrical submersible pump (ESP), as a high-efficiency downhole equipment for converting kinetic energy to hydraulic pressure head, has improved significantly since it was invented in the 1910's by the Russian engineer Arutunoff.Globally, there are more than 100,000 ESP installations, making ESP the second most widely used artificial lift method in oil production, but the first based on the produced amount [1,2].ESPs excel in producing crude oils at very high flow rate, but they have to be operated within a narrow application window.Gas involvement, changing production rate and high oil viscosity can greatly affect the ESP performance [3,4].
Previous studies showed that the presence of gas would cause ESP hydraulic head degradation.The flow behaviors inside ESPs under gassy conditions, such as pressure surging and gas pockets, further deteriorate ESP boosting pressure.Surging may result in vibrations and short service life, and gas pockets can severely limit liquid production rates [5].Although handling gas-liquid mixture has gradually become common for ESPs, the physical mechanism of two-phase flow affecting ESP hydraulic performance is not well understood.The gas bubble formation, coalescence and breakup mechanisms inside ESPs, which affect the two-phase flow characteristics, are still unclear.Due to the compact and complex geometries of multistage ESPs, the visualization of internal flow structures and bubble movement is very difficult.The main objective of this paper is to provide a comprehensive review of literature concerning experimental studies and modeling approaches for predicting ESP pressure increment under gas-liquid flow conditions.The review helps better understand the multiphase flow characteristics inside a rotating ESP, in order to develop mechanistic model to predict ESP boosting pressure under gassy flow conditions.Section 2 presents the fundamental concepts and definitions of a multistage centrifugal pump.Section 3 discusses the previous knowledge of ESP hydraulic performance and flow structures under gassy conditions observed from experiments.Section 4 is divided into four subsections.It begins with the empirical correlations, followed by one-dimensional two-fluid modeling approach.Then, the CFD-based (Computational Fluid Dynamics) modeling and mechanistic approaches are discussed.Finally, the closure relationships to make models solvable are evaluated in Section 5.

Boosting Pressure of Centrifugal Pump
ESPs are widely used in petroleum industry to increase the hydrocarbon fluid production rates, especially for off-shore deep-water oil fields.Normally, multistage ESP is assembled in series (Figure 1a) with each stage comprising the rotating impeller and stationary diffuser (Figure 1b,c).In oil fields, the ESP system usually contains hundreds of stages to meet the boosting pressure requirement due to the great depth of reservoir.The impeller is the dominant part of ESP since it forces the fluid flow and adds kinetic energy to the fluids by spinning the blades.At the outlet of impeller, the main part of kinetic energy in the fluid is converted to pressure potential by diffuser vanes, which are fixed and work as the guiding channels for the impeller of the next stage.
Energies 2018, 11, 180 2 of 41 compact and complex geometries of multistage ESPs, the visualization of internal flow structures and bubble movement is very difficult.
The main objective of this paper is to provide a comprehensive review of literature concerning experimental studies and modeling approaches for predicting ESP pressure increment under gas-liquid flow conditions.The review helps better understand the multiphase flow characteristics inside a rotating ESP, in order to develop mechanistic model to predict ESP boosting pressure under gassy flow conditions.Section 2 presents the fundamental concepts and definitions of a multistage centrifugal pump.Section 3 discusses the previous knowledge of ESP hydraulic performance and flow structures under gassy conditions observed from experiments.Section 4 is divided into four subsections.It begins with the empirical correlations, followed by one-dimensional two-fluid modeling approach.Then, the CFD-based (Computational Fluid Dynamics) modeling and mechanistic approaches are discussed.Finally, the closure relationships to make models solvable are evaluated in Section 5.

Boosting Pressure of Centrifugal Pump
ESPs are widely used in petroleum industry to increase the hydrocarbon fluid production rates, especially for off-shore deep-water oil fields.Normally, multistage ESP is assembled in series (Figure 1a) with each stage comprising the rotating impeller and stationary diffuser (Figure 1b,c).In oil fields, the ESP system usually contains hundreds of stages to meet the boosting pressure requirement due to the great depth of reservoir.The impeller is the dominant part of ESP since it forces the fluid flow and adds kinetic energy to the fluids by spinning the blades.At the outlet of impeller, the main part of kinetic energy in the fluid is converted to pressure potential by diffuser vanes, which are fixed and work as the guiding channels for the impeller of the next stage.As a type of centrifugal pumps, ESPs can be classified into three different categories: radial, axial and mixed types based on the dimensionless specific speed (N S ).The non-dimensional N S is given by: (gH) 3/4  (1 where Ω is the rotational speed (rad/s), Q is liquid flow rate (m 3 /s), g is local gravitational acceleration (m/s 2 ), H is pump head (m).The pump industry uses a more practical expression for the specific speed N S as below: where N, q and h are rotational speed (rpm), flow rate (gpm) and pump head (ft), respectively.Based on Equation ( 2), centrifugal pumps are categorized into radial, mixed and axial types.The radial pumps usually fall in the range 500 < N S < 1800, while the mixed pumps can reach a maximum N S = 4500.For a centrifugal pump, there are three important variables to characterize its hydraulic and mechanical performance, namely: pumping head (H), efficiency (η) and brake horsepower (BHP).Figure 2 below shows the typical pump performance curves provided in the product brochure, which are obtained experimentally using tap-water.
Energies 2018, 11,180 3 of 41 As a type of centrifugal pumps, ESPs can be classified into three different categories: radial, axial and mixed types based on the dimensionless specific speed (NS).The non-dimensional NS is given by: where Ω is the rotational speed (rad/s), Q is liquid flow rate (m 3 /s), g is local gravitational acceleration (m/s 2 ), H is pump head (m).The pump industry uses a more practical expression for the specific speed NS as below: where N, q and h are rotational speed (rpm), flow rate (gpm) and pump head (ft), respectively.Based on Equation (2), centrifugal pumps are categorized into radial, mixed and axial types.The radial pumps usually fall in the range 500 < NS < 1800, while the mixed pumps can reach a maximum NS = 4500.
For a centrifugal pump, there are three important variables to characterize its hydraulic and mechanical performance, namely: pumping head (H), efficiency (η) and brake horsepower (BHP).Figure 2 below shows the typical pump performance curves provided in the product brochure, which are obtained experimentally using tap-water.With the liquid flow rate increase, the ESP performance curves exhibit different trends.The horsepower and pumping head render a monotonic increasing and decreasing trend versus flow rate, respectively, whereas the efficiency is in a semi-quadratic relationship with the pump capacity.The best efficiency point (BEP) presented with the dashed line, corresponds to the highest efficiency of 68.9% at Q = 2700 bpd, N = 3500 rpm in Figure 2. It is an important parameter in characterizing the ESP overall performance.In this study, the main objective is to review different models in predicting ESP hydraulic boosting pressure under gas-liquid flow conditions.

Euler Head
For an ESP impeller, the velocity triangles at the intake and discharge are shown in Figure 3.All the variables are in International System of units (SI).The absolute velocity C can be decoupled into two components: relative velocity W and peripheral velocity U, which is calculated by = Ω × .With the liquid flow rate increase, the ESP performance curves exhibit different trends.The horsepower and pumping head render a monotonic increasing and decreasing trend versus flow rate, respectively, whereas the efficiency is in a semi-quadratic relationship with the pump capacity.The best efficiency point (BEP) presented with the dashed line, corresponds to the highest efficiency of 68.9% at Q = 2700 bpd, N = 3500 rpm in Figure 2. It is an important parameter in characterizing the ESP overall performance.In this study, the main objective is to review different models in predicting ESP hydraulic boosting pressure under gas-liquid flow conditions.

Euler Head
For an ESP impeller, the velocity triangles at the intake and discharge are shown in Figure 3.All the variables are in International System of units (SI).The absolute velocity C can be decoupled into The subscripts 1 and 2 denote the inlet and outlet of impeller, respectively.Here, W is relative to the impeller and C is relative to the global coordinate and is equal to the vector summation of U and W, e.g., 3c, C' 2 is an ideal absolute velocity assuming infinite number of impeller blades, while C 2 is the real absolute velocity at the impeller outlet.The subscripts 1 and 2 denote the inlet and outlet of impeller, respectively.Here, W is relative to the impeller and C is relative to the global coordinate and is equal to the vector summation of U and W, e.g., = + .In Figure 3c, C'2 is an ideal absolute velocity assuming infinite number of impeller blades, while C2 is the real absolute velocity at the impeller outlet.According to the ideal conservation law of angular momentum in rotating centrifugal pump, the Euler head (HE) can be expressed as [6]: where C1U and C2U are the projection of absolute velocities at the impeller inlet and outlet to the direction of peripheral velocities.Applying the velocity triangle relationship in Figure 3, one can obtain: where r is the radius of impeller, h is the channel height, β is the blade angle.If the fluids enter the impeller without pre-rotation, Equation ( 5) can be written as:

Head Loss Mechanisms
As shown in Equation ( 6), the ideal Euler head HE is in linear relationship with liquid flow rate.In reality, pressure losses in the impeller, the diffuser, and losses from the interaction between them must be deducted from the ideal Euler head, including friction, shock and recirculation losses.Therefore, the actual pump head can be calculated by: Figure 4a shows the change of the ideal Euler head with the outlet blade angle β. Figure 4b shows the losses to the ideal head and the final H-Q curve of a centrifugal pump according to Equation (7).Friction loss becomes prominent at high flow rates.In contrast, leakage loss is more at relatively low According to the ideal conservation law of angular momentum in rotating centrifugal pump, the Euler head (H E ) can be expressed as [6]: where C 1U and C 2U are the projection of absolute velocities at the impeller inlet and outlet to the direction of peripheral velocities.Applying the velocity triangle relationship in Figure 3, one can obtain: and: where r is the radius of impeller, h is the channel height, β is the blade angle.If the fluids enter the impeller without pre-rotation, Equation ( 5) can be written as:

Head Loss Mechanisms
As shown in Equation ( 6), the ideal Euler head H E is in linear relationship with liquid flow rate.In reality, pressure losses in the impeller, the diffuser, and losses from the interaction between them must be deducted from the ideal Euler head, including friction, shock and recirculation losses.Therefore, the actual pump head can be calculated by: Energies 2018, 11, 180 5 of 41 Figure 4a shows the change of the ideal Euler head with the outlet blade angle β. Figure 4b shows the losses to the ideal head and the final H-Q curve of a centrifugal pump according to Equation (7).Friction loss becomes prominent at high flow rates.In contrast, leakage loss is more at relatively low flow rates.Shock loss takes place when the liquid flow rate differs from the designed flow rate.Tables 1-5 below discuss each pressure loss mechanism by presenting existing models that are available in literature.Table 1 summarizes the studies on the friction losses inside a centrifugal pump impeller.Takacs [3] pointed out that friction losses progressively increase with liquid rate and are due to fluid friction in the impeller.

Reference
Model Coefficient Ito [8], Jones [9], Churchill [10], Shah [11], Sun [ The mechanism of shock losses in centrifugal pump is still not well studied.Thus, only empirical correlations are available in the literature.Table 2 lists the existing prediction models of shock loss as implemented in the prediction models.Shock losses are negligible at the best efficiency point (BEP) of the pump, but increase at lower and higher liquid flow rates.They occur at the entrance and the exit of the impeller and are caused by sudden changes in the direction of flow [3].

Reference
Model Coefficient Stepanoff [6], Amaral [18] Thin et al.Table 1 summarizes the studies on the friction losses inside a centrifugal pump impeller.Takacs [3] pointed out that friction losses progressively increase with liquid rate and are due to fluid friction in the impeller.

Reference Model Coefficient
Ito [8], Jones [9], Churchill [10], Shah [11], Sun [12] h f riction = Wiesner [13], Sun and Prado [14], The mechanism of shock losses in centrifugal pump is still not well studied.Thus, only empirical correlations are available in the literature.Table 2 lists the existing prediction models of shock loss as implemented in the prediction models.Shock losses are negligible at the best efficiency point (BEP) of the pump, but increase at lower and higher liquid flow rates.They occur at the entrance and the exit of the impeller and are caused by sudden changes in the direction of flow [3].

Reference Model Coefficient
Stepanoff [6], Amaral [18] Thin et al. [15] Wiesner [13] Sun and Prado [14] Thin et al. [15]  As suggested by Tackas [3], the leakage losses represent the losses through the clearances between the rotating and stationary parts of the pump stage (at the impeller eye, through balancing holes, etc.), which decrease with increase of liquid rates.Table 3 below shows the calculation models available in literature.

Bing et al. [17]
In a rotating centrifugal pump, recirculation loss is interpreted as the dissipated fluid energy which is continuously obtained from the shaft due to adverse pressure gradient from the impeller inlet to outlet.As liquid flow rate decreases, the adverse pressure gradient increases, resulting more severe head loss due to fluid recirculation.Currently, the agreement has not been reached as for the definition of recirculation, which in turn a fully understanding in terms of recirculation loss mechanism is unavailable by far.However, the empirical correlations have been proposed in literature, as summarized in Table 4 below.
C 2E -effective outlet velocity Table 5. Diffuser loss models in literature.

Reference Model Coefficient
Ito [8], Jones [9], Churchill [10], Shah [11], Sun [12] Sun and Prado [12,14,22] claimed that the diffuser loss is caused mainly by the friction on the diffuser walls.Table 5 presents the calculation formulae of diffuser losses that were implemented in the prediction models.Most equations are similar to the friction loss formulae in Table 1.
Disk friction losses, as Stepanoff [6] defined, are due to the contact between a rotating disk and fluid.Thus, additional power is consumed to keep the disk rotating since the viscous shear forces act on the disk surfaces [7].Table 6 shows some of existing models of calculating disk friction losses in literature.
By deducting the pressure losses from Euler head in Equation ( 7), the real pressure increment of a centrifugal pump is obtained for single-phase flow.As can be seen in the Tables 1-6, different pressure loss models are available in literature.A proper combination of these models may output an optimized estimation of pump performance.Yet, due to the empirical nature of these models, the universal validity of such combination to deduce the most appropriate prediction model is questioned.Thus, more and more researchers resort to numerical approaches to calculate the hydraulic performance and understand the complicated flow structures inside rotating ESPs.

Experimental Studies
The ESP boosting pressure depends on many hydraulic factors including pump geometry, fluid properties and multi-phase flow conditions.In recent years, with more and more installations of ESP in oil production systems, the effects of fluid viscosity and gas entrainment on ESP's pressure boosting ability have been the focal areas of research interests.

Single-Phase Tests
For single-phase flow, ESP performance is sensitive to the fluid viscosity.With the increase of oil viscosity, ESP boosting pressure becomes lower corresponding to the same flow rate.By the same token, the flow rate through the ESP decreases corresponding to the same boosting pressure in ESPs.Ippen [25] conducted over 200 performance tests for oil viscosities up to 10,000 Saybolt Second Universal (SSU) on four variants of centrifugal pumps.The experimental results were summarized by plotting the ratio of oil head to water head, brake horsepower and efficiency against a Reynolds-type dimensionless number, based on which the general correction factors for specific speeds from 800 to 2200 were proposed to correlate pump's boosting pressure under viscous fluid flow.
Hydraulic Institute provided a typical empirical approach with correction factors to estimate the conventional centrifugal pump boosting pressure for viscous liquid flow if the water performance were known.However, the accuracy of this approach was questioned by Gülich [19,20] and Li [26] since the experiments carried out by Hydraulic Institute were within a narrow range of the pump specific speeds.Unreasonable errors were found if extrapolation was beyond that range.
Stepanoff [6] proposed a similar Reynolds-type number by using only one correction factor to get the new H-Q curves if the water performance were known.A more general model based on the evaluation of viscous dissipation for disk and hydraulic frictions to predict the boosting pressure of centrifugal pumps was proposed by Gülich.The friction losses on disk and in flow passage, as the author claimed, were dominating factors impairing centrifugal pump's ability to handle high viscosity fluids.Compared with available data, Gülich also pointed out that friction losses were affected significantly by pump geometry, fluid properties and thermal conditions.More recent experimental studies conducted by Amaral et al. [27] and Solano [28] further revealed that the Hydraulic Institute charts and empirical correlations available in literature were unable to give appropriate correction factors to predict the ESP boosting pressure for viscous oil flow.The discrepancies between experimental measurements and predictions by existing correlations or charts are discussed in Barrios et al. [29] and Banjar et al. [30].

Two-Phase Tests
The pioneering experimental studies conducted by Murakami and Minemura [31,32] investigated centrifugal pump performance with gas entrainment.By employing a semi-opened impeller pump with a transparent casing, they experimentally observed the behavior of entrained air bubbles.The decreasing total head of the pump caused by air admission and the work consumed for air delivery were reported.Since then, investigators have conducted more experimental measurements and mathematical modeling on centrifugal pump performance under gassy flow conditions.Experimental researches of ESP performance under gas-liquid flow conditions have been carried out extensively at the Tulsa University Artificial Lift Project (TUALP) in the past decades.Table 7 lists a summary of these studies.As can been seen in Table 7, the previous experimental studies, which were conducted by Cirilo [33], Pessoa [35], Beltur [36], Duran [37], Zapata [38], Barrios [39], Gamboa [40], and Salehi [42] among others, covers a wide range of flow conditions.Cirilo built the experimental flow loops for testing ESP performance at the TUALP.Using water and air as the working fluids, he measured both the water and air-water performances of three different types of ESPs as a function of GVFs (the volumetric fraction of gas phase at the ESP inlet), intake pressure and rotational speeds.The obtained data indicated that the mixed type pumps were able to handle much higher GVFs (>30%) than radial type pumps (<10%).With necessary modifications of the same testing flow loops, Pessoa conducted experimental investigations of two-phase performance of ESP using a 22-stage GC6100 pump.By monitoring the stage-by-stage pressure increment, his results revealed that the ESP average behavior was significantly different from that observed per stage.Also, phenomena like surging and gas locking were observed and their boundaries were mapped.Additionally, a second region after pressure surging was observed in mapping test curves, where the slope of pressure increment versus flow rate changed again.
Using the same experimental flow loop of Pessoa, Beltur, Duran, Zapata and Salehi conducted extensive experimental measurements of ESP performance under both liquid-and gas-liquid flow conditions.Beltur focused on ESP performance deterioration with the presence of gas for varying intake GVFs and pressure.Data analysis revealed that the intake GVF is the most important factor affecting ESP boosting pressure under gassy flow.A higher deterioration of pumping head occurs at GVFs above 6%.Duran and Zapata developed empirical correlations for predicting the pressure increment across the stage and the flow regime boundaries.Zapata also carried out further measurements with a wide range of rotational speed to study its effects on the average efficiency of ESP.Salehi did similar measurements using a 14-stage TE2700 ESP, on which the effects of the stage number, intake pressure and fluid properties were investigated.It was found that ESP two-phase boosting pressure varied stage by stage only when the GVF exceeded a certain value, below which the deterioration was mild and independent of stage number.However, the degradation of pumping head became more prominent and affected by the stage number if the GVF reached a critical value, at which the pump boosting pressure of the downstream stage is better than upstream ones.
A typical testing result of ESP two-phase performance with pure water-air flow is shown in Figure 5 [40], including surging (Figure 5a) and mapping tests (Figure 5b).For surging test, the liquid flow rate is fixed, but gas flow rate or inlet gas volumetric fraction (GVF) varies.On the contrary, the liquid flow rate changes while the gas flow rate is fixed in mapping test.In Figure 5a, different curves correspond to different pump intake pressures.Obviously, the pump pressure increment declines non-linearly with the increase of inlet GVF.For mapping test results in Figure 5b, each curve corresponds to a different gas flow rate, which increases from the top to bottom curves.With a fixed gas flow rate, a sudden drop of pump stage pressure increment is observed when flow rate is reduced to a certain value.
Energies 2018, 11, 180 9 of 41 intake GVFs and pressure.Data analysis revealed that the intake GVF is the most important factor affecting ESP boosting pressure under gassy flow.A higher deterioration of pumping head occurs at GVFs above 6%.Duran and Zapata developed empirical correlations for predicting the pressure increment across the stage and the flow regime boundaries.Zapata also carried out further measurements with a wide range of rotational speed to study its effects on the average efficiency of ESP.Salehi did similar measurements using a 14-stage TE2700 ESP, on which the effects of the stage number, intake pressure and fluid properties were investigated.It was found that ESP two-phase boosting pressure varied stage by stage only when the GVF exceeded a certain value, below which the deterioration was mild and independent of stage number.However, the degradation of pumping head became more prominent and affected by the stage number if the GVF reached a critical value, at which the pump boosting pressure of the downstream stage is better than upstream ones.
A typical testing result of ESP two-phase performance with pure water-air flow is shown in Figure 5 [40], including surging (Figure 5a) and mapping tests (Figure 5b).For surging test, the liquid flow rate is fixed, but gas flow rate or inlet gas volumetric fraction (GVF) varies.On the contrary, the liquid flow rate changes while the gas flow rate is fixed in mapping test.In Figure 5a, different curves correspond to different pump intake pressures.Obviously, the pump pressure increment declines non-linearly with the increase of inlet GVF.For mapping test results in Figure 5b, each curve corresponds to a different gas flow rate, which increases from the top to bottom curves.With a fixed gas flow rate, a sudden drop of pump stage pressure increment is observed when flow rate is reduced to a certain value.A more recent experimental study regarding surfactant effect on ESP gas-handling ability and boosting pressure under gassy flow was conducted by Zhu et al. [44].Surfactants are molecules with a hydrophobic and a hydrophilic part, and therefore they preferentially adsorb at the interface of continuous/dispersed phases.In the process, they reduce the interfacial tension [45].Studies of surfactant effect on bubble/droplet formation in air/water or water/oil binary immiscible two-phase flows were carried out by Hu et al. [46], Omer and Pal [47].Surfactant effect on liquid loading in gas well has been studied by van Nimwegen et al. [48,49], Ajani et al. [50,51].Both reported a lower critical gas velocity at which liquid loading occurs, meaning that the surfactant presence helps prevent liquid loading in wells.
However, the studies on the surfactant effect on centrifugal pump performance under two-phase flow conditions are few.Zhu et al. [44] measured the ESP two-phase pressure increment with/without surfactant injection in the flow loop.The surfactant used in their study is isopropanol alcohol (IPA).The decline trend of normalize pressure (Np) in Figure 6a versus the intake GVF improves with surfactant injections, and the pressure surging phenomenon depicted by a sudden drop of ESP boosting pressure disappears.Figure 6b compares ESP H-Q performance curves under air-water flow with different IPA volumetric concentrations.A clear difference can be seen at low QL in terms of the stage pressure increment.Without IPA injection, the pressure increment drops to zero if QL becomes very low.But, the ESP performance improves significantly with surfactant presence.A more recent experimental study regarding surfactant effect on ESP gas-handling ability and boosting pressure under gassy flow was conducted by Zhu et al. [44].Surfactants are molecules with a hydrophobic and a hydrophilic part, and therefore they preferentially adsorb at the interface of continuous/dispersed phases.In the process, they reduce the interfacial tension [45].Studies of surfactant effect on bubble/droplet formation in air/water or water/oil binary immiscible two-phase flows were carried out by Hu et al. [46], Omer and Pal [47].Surfactant effect on liquid loading in gas well has been studied by van Nimwegen et al. [48,49], Ajani et al. [50,51].Both reported a lower critical gas velocity at which liquid loading occurs, meaning that the surfactant presence helps prevent liquid loading in wells.
However, the studies on the surfactant effect on centrifugal pump performance under two-phase flow conditions are few.Zhu et al. [44] measured the ESP two-phase pressure increment with/without surfactant injection in the flow loop.The surfactant used in their study is isopropanol alcohol (IPA).The decline trend of normalize pressure (N p ) in Figure 6a versus the intake GVF improves with surfactant injections, and the pressure surging phenomenon depicted by a sudden drop of ESP boosting pressure disappears.Figure 6b compares ESP H-Q performance curves under air-water flow with different IPA volumetric concentrations.A clear difference can be seen at low Q L in terms of the stage pressure increment.Without IPA injection, the pressure increment drops to zero if Q L becomes very low.But, the ESP performance improves significantly with surfactant presence.

Flow Pattern Recognition
In multiphase pipe flow, flow pattern varies due to phase interaction and interfacial momentum transfer changes, which depend on fluid properties, flow rates and pipe geometry.Many studies on the recognition of flow patterns and their transition boundaries are available in literature [52][53][54].Flow pattern is important to accurately predict the liquid holdup and pressure gradient in pipe flow.Take the upward gas/liquid flow in a vertical pipe for example, the main flow patterns are bubble, slug, churn and annular flows as illustrated in Figure 7a.Correspondingly, the flow pattern map showing the transition boundaries between different flow patterns is plotted in Figure 7b with the superficial velocities of gas and liquid as its coordinates.
The dispersed bubble flow prevails if the gas flow rate is low and the liquid flow rate is high, as shown in Figure 7a.For this flow pattern, the turbulence forces are high enough to breakup gas-phase into small bubbles and disperse them in the continuous liquid-phase [52].
When both gas flow rate and liquid flow rate are low, gas bubbles become larger and move upward in a zigzag path.The bubble flow regime may or may not exist depending on the pipe diameter.Confined by maximum lattice packing, the in-situ gas void fraction of bubble flow should be below 0.52 [55].However, for a practical application, the maximum gas void fraction below 0.25 is commonly used for vertical co-current pipe flow.

Flow Pattern Recognition
In multiphase pipe flow, flow pattern varies due to phase interaction and interfacial momentum transfer changes, which depend on fluid properties, flow rates and pipe geometry.Many studies on the recognition of flow patterns and their transition boundaries are available in literature [52][53][54].Flow pattern is important to accurately predict the liquid holdup and pressure gradient in pipe flow.Take the upward gas/liquid flow in a vertical pipe for example, the main flow patterns are bubble, slug, churn and annular flows as illustrated in Figure 7a.Correspondingly, the flow pattern map showing the transition boundaries between different flow patterns is plotted in Figure 7b with the superficial velocities of gas and liquid as its coordinates.
The dispersed bubble flow prevails if the gas flow rate is low and the liquid flow rate is high, as shown in Figure 7a.For this flow pattern, the turbulence forces are high enough to breakup gas-phase into small bubbles and disperse them in the continuous liquid-phase [52].
When both gas flow rate and liquid flow rate are low, gas bubbles become larger and move upward in a zigzag path.The bubble flow regime may or may not exist depending on the pipe diameter.Confined by maximum lattice packing, the in-situ gas void fraction of bubble flow should be below 0.52 [55].However, for a practical application, the maximum gas void fraction below 0.25 is commonly used for vertical co-current pipe flow.

Flow Pattern Recognition
In multiphase pipe flow, flow pattern varies due to phase interaction and interfacial momentum transfer changes, which depend on fluid properties, flow rates and pipe geometry.Many studies on the recognition of flow patterns and their transition boundaries are available in literature [52][53][54].Flow pattern is important to accurately predict the liquid holdup and pressure gradient in pipe flow.Take the upward gas/liquid flow in a vertical pipe for example, the main flow patterns are bubble, slug, churn and annular flows as illustrated in Figure 7a.Correspondingly, the flow pattern map showing the transition boundaries between different flow patterns is plotted in Figure 7b with the superficial velocities of gas and liquid as its coordinates.
The dispersed bubble flow prevails if the gas flow rate is low and the liquid flow rate is high, as shown in Figure 7a.For this flow pattern, the turbulence forces are high enough to breakup gas-phase into small bubbles and disperse them in the continuous liquid-phase [52].
When both gas flow rate and liquid flow rate are low, gas bubbles become larger and move upward in a zigzag path.The bubble flow regime may or may not exist depending on the pipe diameter.Confined by maximum lattice packing, the in-situ gas void fraction of bubble flow should be below 0.52 [55].However, for a practical application, the maximum gas void fraction below 0.25 is commonly used for vertical co-current pipe flow.Slug flow regime corresponds to the steady alternating flow of gas pockets (Taylor bubble) and liquid slugs [53] at even higher gas flow rates (see Figure 7b), where Taylor bubbles occupy the entire cross section of the pipe except for a thin liquid film on the wall.Slug flow can be considered as steady-state flow pattern since its slug length, frequency, translational velocity etc. are time-independent.However, a further increase of gas flow rates triggers the collapse of slugs and lead to an unstable flow that is described as the churn flow regime due to the large degree of turbulence in the flow [54].
Annular flow is characterized by extremely high gas flow rates and continuous gas core flowing upward with thin liquid film on the pipe wall.Due to high interfacial shear and drag force effects, the liquid droplets from the wavy interface between gas core and liquid film can be swept into the gas core, and pushed upward by the high-speed gas phase.
Several visualization studies have been conducted to observe multiphase flow structures/patterns inside a rotating centrifugal pump.Table 8 below lists the visualization experiments on flow patterns inside a centrifugal pump in literature.Most of these studies used transparent casing to visually observe the pump internal flow patterns.Starting with the original volute-type centrifugal pumps [31,32,[56][57][58][59][60][61][62][63][64][65][66][67][68][69][70], they cover several research topics including flow pattern recognition, bubble movement visualization, mapping transition boundaries etc. Murakami and Minemura [31,32] experimental study first investigated the gas entrainment effect on the hydraulic pump head of a centrifugal pump visually, and associated the degradation of pump performance with gas-liquid flow pattern inside the rotor.Using a transparent casing, Estevam [62] conducted the first visualization experiment of flow pattern recognition inside a rotating ESP impeller.Later visualization experimental work, conducted by Gamboa [40] presented similar observations as Estevam's study, based on which the first flow pattern map with transition boundaries between different flow patterns inside an ESP was proposed.His work offers significant insight regarding how to characterize the flow patterns in an ESP quantitatively.
Among previous visualization experiments, two recent studies from Schäfer et al. [68] and Neumann et al. [69] are different from the others.They both used a new non-intrusive technique, the so called high resolution computed tomography (HireCT) to accurately measure the distribution of in-situ gas void fractions.Compared to the intrusive technique, the new measurement approach, as the authors claimed, can provide better image resolution of the gas void distribution and thus a higher measurement accuracy.
By visualizing the ESP internal flow, Barrios [39], Barrios and Prado [67] observed that bubbles enlarged when inlet GVFs increased and pump rotational speeds decreased.Such enlargement corresponding to the poorer pump performance indicated that bubble behaviors played a significant role in ESP's ability of handling gas-liquid mixtures.In addition, visualization experiments also showed different flow patterns prevailing inside ESP channels at higher GVFs.
In Figure 8, λ G denotes the intake GVF.It is evident that flow behaviors inside ESP impeller altered significantly with flow conditions.From Figure 8a,b, the GVF increase caused formation of larger bubbles and gas pockets, which in turn choked the flow passage for liquids and decreased the pump hydraulic head.Moreover, the flow patterns prevailing in ESP impeller at a specific value of λ G , are comparable to that of gas-liquid two-phase flow in pipelines.[69] are different from the others.They both used a new non-intrusive technique, the so called high resolution computed tomography (HireCT) to accurately measure the distribution of in-situ gas void fractions.Compared to the intrusive technique, the new measurement approach, as the authors claimed, can provide better image resolution of the gas void distribution and thus a higher measurement accuracy.
By visualizing the ESP internal flow, Barrios [39], Barrios and Prado [67] observed that bubbles enlarged when inlet GVFs increased and pump rotational speeds decreased.Such enlargement corresponding to the poorer pump performance indicated that bubble behaviors played a significant role in ESP's ability of handling gas-liquid mixtures.In addition, visualization experiments also showed different flow patterns prevailing inside ESP channels at higher GVFs.
In Figure 8, λG denotes the intake GVF.It is evident that flow behaviors inside ESP impeller altered significantly with flow conditions.From Figure 8a,b, the GVF increase caused formation of larger bubbles and gas pockets, which in turn choked the flow passage for liquids and decreased the pump hydraulic head.Moreover, the flow patterns prevailing in ESP impeller at a specific value of λG, are comparable to that of gas-liquid two-phase flow in pipelines.Similar as the gas/liquid pipe flow, Estevam [64] and Estevam et al. [71] categorized the flow patterns in a rotating ESP impeller as bubbly flow, transition flow, and elongated-bubble flow.A more thorough experimental study on flow pattern recognition of two-phase flow in ESP impeller was carried out by Gamboa [40], Gamboa and Prado [72].Under a given flow condition, the typical flow pattern map from their visualization experiments is shown in Figure 9, where the curves denoted by different symbols represent the transition boundaries between flow patterns, including homogenous flow, bubbly flow, gas-pocket formation, and segregated flow.Figure 9 is important to understand the two-phase flow mechanisms in a rotating ESP.From the perspective of mechanistic modeling, each flow regime corresponds to the specific governing equations for flow characteristics, such as bubble size (d b ), in-situ gas void fraction (α G ), slippage velocity (V slip ) between gas and liquid phases.With the flow pattern prevailing inside a rotating ESP determined, the governing equations based on mass and momentum conservations can be simplified.
Energies 2018, 11, 180 13 of 41 Similar as the gas/liquid pipe flow, Estevam [64] and Estevam et al. [71] categorized the flow patterns in a rotating ESP impeller as bubbly flow, transition flow, and elongated-bubble flow.A more thorough experimental study on flow pattern recognition of two-phase flow in ESP impeller was carried out by Gamboa [40], Gamboa and Prado [72].Under a given flow condition, the typical flow pattern map from their visualization experiments is shown in Figure 9, where the curves denoted by different symbols represent the transition boundaries between flow patterns, including homogenous flow, bubbly flow, gas-pocket formation, and segregated flow.Figure 9 is important to understand the two-phase flow mechanisms in a rotating ESP.From the perspective of mechanistic modeling, each flow regime corresponds to the specific governing equations for flow characteristics, such as bubble size (db), in-situ gas void fraction (αG), slippage velocity (Vslip) between gas and liquid phases.With the flow pattern prevailing inside a rotating ESP determined, the governing equations based on mass and momentum conservations can be simplified.Verde et al. [70] conducted visualization experiments on flow pattern recognition inside a rotating ESP impeller using high speed and resolution imaging technique.As shown in Figure 10, four different flow patterns were classified, including bubble flow (Figure 10a), agglomerated bubble flow (Figure 10b), gas pocket flow (Figure 10c) and segregated flow (Figure 10d).They observed that the intensity of pump performance degradation is directly influenced by the flow pattern within the impeller.The occurrence of the gas pocket flow pattern is linked to the intensification of the deterioration of pump performance and the appearance of operation instabilities.Moreover, the segregated flow patterns correspond to the severe performance degradation which makes the pump incapable of generating pressure.
Summarized by Verde et al. [70], the schematic representations of each flow pattern are shown in Figure 11 with intake GVF increasing from left to right.Due to the relatively small intake GVF, the homogenous flow regime is featured by tiny and evenly-dispersed bubbles inside impeller channels, as shown in Figure 11a.Under this regime, the bubbles are deemed to move together with liquid phase.Slippage between gas and liquid is small, meaning the in-situ αG is almost the same as λG.
As the intake GVF increases, the tiny bubbles are prone to collide and aggregate to form bigger ones, resulting in bubbly flow regime.In contrast to homogenous flow regime, the phase slippage between gas and liquid, shown by Figure 11b can no longer be neglected.Thus, depending on the intake GVF, the in-situ αG under bubbly flow becomes higher than λG.A further increase of GVF induces more severe collision and aggregation of bubbles so that the large gas pocket forms.This  Verde et al. [70] conducted visualization experiments on flow pattern recognition inside a rotating ESP impeller using high speed and resolution imaging technique.As shown in Figure 10, four different flow patterns were classified, including bubble flow (Figure 10a), agglomerated bubble flow (Figure 10b), gas pocket flow (Figure 10c) and segregated flow (Figure 10d).They observed that the intensity of pump performance degradation is directly influenced by the flow pattern within the impeller.The occurrence of the gas pocket flow pattern is linked to the intensification of the deterioration of pump performance and the appearance of operation instabilities.Moreover, the segregated flow patterns correspond to the severe performance degradation which makes the pump incapable of generating pressure.
Summarized by Verde et al. [70], the schematic representations of each flow pattern are shown in Figure 11 with intake GVF increasing from left to right.Due to the relatively small intake GVF, the homogenous flow regime is featured by tiny and evenly-dispersed bubbles inside impeller channels, as shown in Figure 11a.Under this regime, the bubbles are deemed to move together with liquid phase.Slippage between gas and liquid is small, meaning the in-situ α G is almost the same as λ G .
As the intake GVF increases, the tiny bubbles are prone to collide and aggregate to form bigger ones, resulting in bubbly flow regime.In contrast to homogenous flow regime, the phase slippage Energies 2018, 11, 180 14 of 41 between gas and liquid, shown by Figure 11b can no longer be neglected.Thus, depending on the intake GVF, the in-situ α G under bubbly flow becomes higher than λ G .A further increase of GVF induces more severe collision and aggregation of bubbles so that the large gas pocket forms.This flow pattern is similar as slug/churn flow patterns in pipelines, which are featured by a gas core/pocket followed by a liquid slug.As shown in Figure 11c, the Taylor-bubble-like gas pocket forms near the suction side of ESP impeller, which occupies a significant portion of the impeller flow passage.
Energies 2018, 11, 180 14 of 41 flow pattern is similar as slug/churn flow patterns in pipelines, which are featured by a gas core/pocket followed by a liquid slug.As shown in Figure 11c, the Taylor-bubble-like gas pocket forms near the suction side of ESP impeller, which occupies a significant portion of the impeller flow passage.A relatively high gas flow rate and low liquid flow rate lead to the segregated flow pattern, similar to the annular/stratified flows in pipe.As it can be seen in Figure 11d, the elongated bubble expands and occupies the entire impeller length.For practical applications, the segregated flow pattern corresponds to gas-locking, an operational failure which gives null pump head and null flow rate [70].
As a direct observation of flow patterns inside the ESP impeller, visualization experiments can help reveal gas-liquid flow behaviors intuitively.However, the experimental facility needs special designs associated with necessary modifications on pump geometries, such as the removal of A relatively high gas flow rate and low liquid flow rate lead to the segregated flow pattern, similar to the annular/stratified flows in pipe.As it can be seen in Figure 11d, the elongated bubble expands and occupies the entire impeller length.For practical applications, the segregated flow pattern corresponds to gas-locking, an operational failure which gives null pump head and null flow rate [70].
As a direct observation of flow patterns inside the ESP impeller, visualization experiments can help reveal gas-liquid flow behaviors intuitively.However, the experimental facility needs special designs associated with necessary modifications on pump geometries, such as the removal of A relatively high gas flow rate and low liquid flow rate lead to the segregated flow pattern, similar to the annular/stratified flows in pipe.As it can be seen in Figure 11d, the elongated bubble expands and occupies the entire impeller length.For practical applications, the segregated flow pattern corresponds to gas-locking, an operational failure which gives null pump head and null flow rate [70].
As a direct observation of flow patterns inside the ESP impeller, visualization experiments can help reveal gas-liquid flow behaviors intuitively.However, the experimental facility needs special designs associated with necessary modifications on pump geometries, such as the removal of impeller hub and attaching the Pyrex glass on its top for visualized observation.Barrios [39] pointed out that the visualization of two-phase flow structures in a multistage ESP assembled in series was much more difficult.Thus, how to characterize the multiphase flow in ESPs has become a challenging topic.Although several technologies [68,69,[73][74][75] to visualize the internal flows inside a centrifugal pump were discussed, they required some modifications on pump geometry and the implementation was difficult to carry out.It required mounting HireCT into the apparatus so that the internal flow structures could be visualized regardless of the opaque pump casing or volute.In addition, the data processing involved time-averaged rotation-synchronized CT scanning techniques, adding further complexity in analyzing the obtained experimental results.
Figure 12 shows the measurements using the HireCT technique by Schäfer et al. [68].The horizontal axis denotes the inlet GVF and the vertical axis is the volumetric averaged in-situ α G in the rotating pump impeller.Clearly, there is a sharp α G increase at a GVF corresponding to the severe gas-pocket formation and pumping head degradation.The step change is between inlet GVF = 2.5% and 3. impeller hub and attaching the Pyrex glass on its top for visualized observation.Barrios [39] pointed out that the visualization of two-phase flow structures in a multistage ESP assembled in series was much more difficult.Thus, how to characterize the multiphase flow in ESPs has become a challenging topic.Although several technologies [68,69,[73][74][75] to visualize the internal flows inside a centrifugal pump were discussed, they required some modifications on pump geometry and the implementation was difficult to carry out.It required mounting HireCT into the apparatus so that the internal flow structures could be visualized regardless of the opaque pump casing or volute.In addition, the data processing involved time-averaged rotation-synchronized CT scanning techniques, adding further complexity in analyzing the obtained experimental results.
Figure 12 shows the measurements using the HireCT technique by Schäfer et al. [68].The horizontal axis denotes the inlet GVF and the vertical axis is the volumetric averaged in-situ αG in the rotating pump impeller.Clearly, there is a sharp αG increase at a GVF corresponding to the severe gas-pocket formation and pumping head degradation.The step change is between inlet GVF = 2.5% and 3.The aforementioned experimental studies on ESP gas-liquid performance used tap water as the working fluid, while the gas phase was compressed air or nitrogen.Several recent experimental studies focusing on ESP gas-handling ability under viscous fluid flow were conducted by Trevisan [41], Trevisan and Prado [76], Banjar et al. [30], and Paternost et al. [77].Using a visualization prototype built from original ESP components and with minimal geometrical modifications, Trevisan [41], Trevisan and Prado [76] conducted experiments to investigate the viscosity effect on liquid/gas two-phase flow through ESP.The authors identified four liquid/air flow patterns inside the impeller channels: agglomerated bubble, gas pocket, segregated gas and intermittent gas flows.It was concluded that the agglomerated bubble flow is responsible for pressure surging phenomenon, which is the initiation of pump head deterioration due to gas entrainment.The authors also observed that the increase in viscosity caused surging to occur at relatively lower inlet GVFs.Similar experimental observations were made by Banjar et al. [30].Paternost et al. [77] further investigated the performance of a centrifugal pump handling single-phase viscous liquids and analyzed the impact of free gas entrainment.They concluded that the degradation of pump head was due to the stagnation of large gas-pocket formation, which became worse with the increase of liquid viscosity.An empirical correlation similar to the dimension analysis procedure in Solano [28] was proposed to correlate ESP pressure increment under two-phase flow conditions accounting the inlet GVF and fluid viscosity.The aforementioned experimental studies on ESP gas-liquid performance used tap water as the working fluid, while the gas phase was compressed air or nitrogen.Several recent experimental studies focusing on ESP gas-handling ability under viscous fluid flow were conducted by Trevisan [41], Trevisan and Prado [76], Banjar et al. [30], and Paternost et al. [77].Using a visualization prototype built from original ESP components and with minimal geometrical modifications, Trevisan [41], Trevisan and Prado [76] conducted experiments to investigate the viscosity effect on liquid/gas two-phase flow through ESP.The authors identified four liquid/air flow patterns inside the impeller channels: agglomerated bubble, gas pocket, segregated gas and intermittent gas flows.It was concluded that the agglomerated bubble flow is responsible for pressure surging phenomenon, which is the initiation of pump head deterioration due to gas entrainment.The authors also observed that the increase in viscosity caused surging to occur at relatively lower inlet GVFs.Similar experimental observations were made by Banjar et al. [30].Paternost et al. [77] further investigated the performance of a centrifugal pump handling single-phase viscous liquids and analyzed the impact of free gas entrainment.They concluded that the degradation of pump head was due to the stagnation of large gas-pocket formation, which became worse with the increase of liquid viscosity.An empirical correlation similar to the dimension analysis procedure in Solano [28] was proposed to correlate ESP pressure increment under two-phase flow conditions accounting the inlet GVF and fluid viscosity.

Modeling Approaches
In this section, several empirical correlations for calculating centrifugal pump performance under gas/liquid two-phase flow are reviewed.Although empirical correlations are much easier to implement compared to other numerical or mechanistic models, the application range is rather limited.With this consideration, the more sophisticated one-dimensional two fluid models coupled with interfacial momentum transfer mechanisms are discussed.The 3D CFD simulations, based on the conservation laws of mass and momentum, can provide detailed numerical results in terms of flow structure, phase distribution and slippage etc.However, CFD simulations for multiphase flow, which usually demand high computational cost, are not as reliable or robust as for single-phase flow.Furthermore, both one-dimensional two-fluid model and CFD-based numerical model are extremely sensitive to the closure relationships, such as bubble size, drag force model and other interphase momentum transfer terms, etc.Thus, the mechanistic models rooted in physics, including flow pattern, mass and momentum conservation equations, have become a new direction for developing the best prediction model of boosting pressure in ESPs with gas involvement.

Emperical Correlations
The first empirical correlation based on experimental data of ESP pressure degradation due to free gas entrainment was proposed by Turpin et al. [78]: where, H m is the head of ESP with gas/liquid flow (ft); H is the single-phase pump head (ft); q g and q l are gas and liquid flow rates (bpd), respectively; p in is pump intake pressure (psi); C 1 = 0.145 is a unit conversion constant.As the authors claimed, a critical constraint of stable flow condition inside ESPs should be applied to this model so as to achieve a reasonable prediction accuracy.Thus, another empirical correlation to distinguish the stable/unstable flow boundaries are given as below: The stable flow condition corresponds to φ < 1, which is the recommended application range of Equation (8).
Romero [34] established a multiphase head model for the type of mixed flow pumps based on Cirilo's [33] experimental data: where H max is the shut-in pump head in ESP; q ld is the dimensionless liquid flow rate, a ratio of the intake liquid flow rate to the open flow rate (OFR).The application range for Equation (10), as pointed out by the authors, is dispersed bubble flow or low degradation of ESP performance: To extend the application range of above two empirical correlations, Duran and Prado [79] presented models for the mild and severe head degradations in ESPs, which correspond to bubbly and elongated bubble flow regimes: where ∆p m is the pressure increment of single stage ESP with gas/liquid flow (psi), α is the in-situ gas void fraction.Equation ( 13) is applied to bubbly flow and elongated bubble flow, respectively.A closure relationship is presented below to solve for α in above equations: q max 1.622 0.435 q l (1−α) 1.6213q max (14) Similarly, Equation ( 14) corresponds to Equation ( 13), which can be used to calculate the in-situ gas void fraction and then the pressure increment in ESP with gas presence can be obtained accordingly.Although the application range of Equations ( 13) has been extended compared to Equation ( 10), its formulation is only based on the experimental data from a single ESP geometry without any validation against different pump models.Thus, its general validity is questioned due to its empirical nature.

One-Dimension Two-Fluid Modeling
Zakem [80] first developed a mathematical model using a one-dimensional control volume method to analyze gas bubbles and liquid interaction for straight blade impellers.Furuya [81] developed a similar analytical model by incorporating the pump geometry, void fraction, flow slippage, and the flow regime but neglecting the compressibility and condensation effects.The basic formula in Furuya study is given by: .
. m l and .m g are the mass flow rates of liquid and gas.s is the streamline coordinate.n is the coordinate normal to the streamline coordinate s. β is the angle between relative flow velocity and circumferential direction.γ is the angle between the radial direction and the stream surface.C d is the drag coefficient.r b is the bubble radius.Comparing with experimental data in literature, the model predictions were within the relative average error band of ±30% for GVF <20%, and ±50% for GVF >30%.
Sachdeva et al. [82,83] conducted a comprehensive investigation of two-phase flow in ESPs with air/water and diesel/CO2 mixtures.A dynamic five-equation, one-dimensional, two-fluid model accounting for pump geometry, intake pressure and GVF, fluid properties, was developed and verified to predict ESP boosting pressure.The basic equations of Sachdeva model are presented below.
The mass balance equations are: where r is radius in ESP, W G and W L are the mass flow rates of gas and liquid, respectively.
The momentum balance equations are: where V R,G and V R,L are the radial components of gas and liquid absolute velocities, respectively.F W,G and F W,L are the friction forces of gas and liquid against the channel walls per unit volume of fluids.F i is the interfacial friction force between gas and liquid per unit volume.F v is the virtual mass force per unit volume.Minemura et al. [84] also studied the performance of centrifugal pumps in the nuclear industry under air-water two phase flow conditions with a low inlet GVF (<10%).Based on energy change in the flow from the rotating impeller to the stationary volute, a 1D, two-fluid model considering fluid viscosity and air-phase compressibility in a rotating impeller was proposed.This model can be solved numerically with a prediction error of ±20% of the related flow capacity.However, both Sachdeva et al. and Minemura et al. models are valid only under a narrow application range or specific experimental conditions.Compared to Sachdeva model, the major difference in Minemura et al. model is that the momentum balance equations are based on the relative velocities of gas/liquid rather than the radial components of the absolute velocities, which are given as: where the subscript p is g or l for gas or liquid phase, respectively.And M p,s is the interfacial momentum transfer term.Eliminating the pressure increment term at the left hand side in Equation ( 24), the combined momentum balance equation can be expressed as Applying a finite differencing algorithm, a preferred numerical scheme [85], to solving this model, a good agreement of pump performance curve, α G distribution, surging and gas lock conditions against correspondent experimental measurements was obtained.Although many studies for modeling ESP performance under gassy conditions were conducted, mechanistic modeling is still needed due to the over-simplification and assumptions or narrow application ranges of the existing models.

Three-Dimension Computational Fluid Dynamics (CFD)
With the advances of computer technology, computational fluid dynamics (CFD) becomes a more and more powerful tool to study centrifugal-pump performance under single-phase and multiphase flow conditions.Due to complicated ESP geometries, it is difficult to investigate the internal velocity and pressure fields experimentally.However, CFD offers an alternative way to simulate the complex internal flow structures.A centrifugal pump consists of an impeller rotating at a set angular velocity and a volute which is stationary.For an ESP, the rotating and stationary parts are the impeller and diffuser, which are accommodated in the rotating and stationary computational domains in CFD, respectively.
Using a 3D CFD code with the frozen-rotor interface model, the internal flow inside centrifugal pumps can be simulated, including velocity and pressure fields [86][87][88] as well as flow recirculation and separation [89,90].The frozen-rotor model is considered as steady-state simulation because it holds the rotating and stationary parts in two separate reference frames.Some transient simulations were conducted using the sliding-mesh technique to investigate the dynamic flow structures in centrifugal pumps [91][92][93][94].Gonzalez et al. [91] performed numerical simulations of unsteady flow in a single-phase centrifugal pump, considering impeller-volute interaction.By solving viscous, incompressible Navier-Stokes (N-S) equations with the sliding mesh technique, the unsteady flow behavior inside a centrifugal pump due to impeller-volute interaction was captured.A relationship between the global variables, such as torque, as a function of impeller relative position, secondary flow in volute, etc., was obtained numerically.This approach was successful in predicting the dynamic interaction between impeller and volute.Huang et al. [93] studied unsteady flow and pressure fluctuations due to interaction between impeller and diffuser vanes by the sliding mesh technique.His study confirmed that the global variables are primarily affected by impeller blade passing frequency.
In addition to designing turbomachinery, CFD can help engineers study the viscosity effects on centrifugal pump performance.Shojaeefard et al. [95,96] conducted both experimental study and numerical simulation of a centrifugal pump handling viscous fluids.The authors stated that a good agreement between simulation and experimental data was obtained by solving the steady state Reynolds-averaged Navier-Stokes (RANS) equations with Shear Stress Transport (SST) k-ω turbulence model.Using the same pump geometry, Sirino et al. [97] and Stel et al. [98] performed numerical investigations of viscosity effects on single-stage and three-stage ESPs, respectively.Similar numerical methodologies were used in their work including SST turbulence model with transient rotor-stator techniques.Both simulation results matched experiments well under a wide range of fluid viscosity.Stel et al. [98,99] pointed out that CFD simulated boosting pressure with multistage ESP geometry agreed with experimental results better than that with a single-stage geometry.The phenomenon of rising head with moderate increase of fluid viscosity was studied by Li [100].By implementing the standard k-ε turbulence model and non-equilibrium wall function into RANS equations, the author confirmed that the rising pump head was due to the transition from hydraulically rough regime to hydraulically smooth regime.Zhu et al. [84] solved a set of 3D, steady-state RANS equations with standard SST turbulence model using the commercial CFD solver ANSYS CFX by employing the frozen-rotor technique.Their results matched correspondent experimental data well.
CFD has also been used to simulate pump performance with gas involvement, such as cavitation phenomenon [101,102], and free-gas entrainment flow [103][104][105].Unlike single-phase simulations, the two-phase simulation requires the solution of conservation equations of mass, momentum and energy for the continuous and dispersed phases.Meanwhile, another set of constitutive equations, for describing the interphase interactions, such as interfacial momentum/heat transfer, need to be solved simultaneously.
Minemura and Uchiyama [106] solved 3D momentum equations with a finite-element method to predict gas/liquid flow behavior in a rotating centrifugal pump impeller.Tremante et al. [107] numerically studied gas/liquid flow through a cascade axial pump by CFD simulation.Coupled with a modified k-ε turbulence model by considering the viscosity of the liquid phase and the compressibility of the gas phase, the gas/liquid distributions versus different attack angles were obtained.Caridad et al. [108,109] studied ESP impeller performance handling water/air mixtures using CFD simulation.Applying two fluid models in 3D CFD simulations, the pressure and velocity fields, as well as the gas phase distributions (Figure 13) were obtained.The gas pocket near the impeller blade was also identified and compared with experimental observations.Minemura and Uchiyama [106] solved 3D momentum equations with a finite-element method to predict gas/liquid flow behavior in a rotating centrifugal pump impeller.Tremante et al. [107] numerically studied gas/liquid flow through a cascade axial pump by CFD simulation.Coupled with a modified k-ε turbulence model by considering the viscosity of the liquid phase and the compressibility of the gas phase, the gas/liquid distributions versus different attack angles were obtained.Caridad et al. [108,109] studied ESP impeller performance handling water/air mixtures using CFD simulation.Applying two fluid models in 3D CFD simulations, the pressure and velocity fields, as well as the gas phase distributions (Figure 13) were obtained.The gas pocket near the impeller blade was also identified and compared with experimental observations.The sensitivity analysis on GVF and bubble diameter indicated that ESP performance deteriorated with GVF and bubble diameter increase.Barrios [39], Barrios et al. [103] conducted multiphase CFD simulations on a single-stage ESP impeller with the new models of bubble size and drag coefficient predictions.Their simulations agreed with laboratory visualization images of streamlines and gas-accumulation zones.Qi et al. [110] designed ESPs for geothermal application with high-temperature gas-liquid two-phase flow.Using CFD simulations, the designed mix-type centrifugal impeller and diffuser were optimized for better gas-handling and higher efficiency within a wide range of production rate.Zhu and Zhang [104] conducted multiphase CFD simulations on a three-stage ESP with each stage comprising of an impeller and a diffuser.Comparing to experimental measurements, their work revealed that bubble size was a critical factor affecting ESP performance under gassy conditions.A new bubble size prediction model was later proposed in Zhu and Zhang [105] based on CFD simulations.They also predicted the in-situ gas void fraction (αG) with a mechanistic model and validated the results with numerically simulated values (Zhu and Zhang [90]).A typical CFD simulation result of gas/liquid flow inside a rotating ESP from their studies is shown in Figure 14 below.
Multiphase flow phenomena in ESP are transient in nature, such as fluctuations of local pressure, gas void fractions, and breakup or coalescence of bubbles.To better simulate the hydrodynamics of gas-liquid two-phase flow in a rotating centrifugal pump, the unsteady CFD simulation code coupled with multiphase flow model and transient rotor-stator algorithm to account for the interactions between impeller and diffuser should be adopted.However, the computational cost for transient CFD simulations is far more than that for steady-state simulations.Marsis et al. [111] performed transient two-phase CFD simulations on eight multi-vane ESP designs.The predicted pump performance was confirmed by experimental results.The final design was achieved with the stage The sensitivity analysis on GVF and bubble diameter indicated that ESP performance deteriorated with GVF and bubble diameter increase.Barrios [39], Barrios et al. [103] conducted multiphase CFD simulations on a single-stage ESP impeller with the new models of bubble size and drag coefficient predictions.Their simulations agreed with laboratory visualization images of streamlines and gas-accumulation zones.Qi et al. [110] designed ESPs for geothermal application with high-temperature gas-liquid two-phase flow.Using CFD simulations, the designed mix-type centrifugal impeller and diffuser were optimized for better gas-handling and higher efficiency within a wide range of production rate.Zhu and Zhang [104] conducted multiphase CFD simulations on a three-stage ESP with each stage comprising of an impeller and a diffuser.Comparing to experimental measurements, their work revealed that bubble size was a critical factor affecting ESP performance under gassy conditions.A new bubble size prediction model was later proposed in Zhu and Zhang [105] based on CFD simulations.They also predicted the in-situ gas void fraction (α G ) with a mechanistic model and validated the results with numerically simulated values (Zhu and Zhang [90]).A typical CFD simulation result of gas/liquid flow inside a rotating ESP from their studies is shown in Figure 14 below.
Multiphase flow phenomena in ESP are transient in nature, such as fluctuations of local pressure, gas void fractions, and breakup or coalescence of bubbles.To better simulate the hydrodynamics of gas-liquid two-phase flow in a rotating centrifugal pump, the unsteady CFD simulation code coupled with multiphase flow model and transient rotor-stator algorithm to account for the interactions between impeller and diffuser should be adopted.However, the computational cost for transient CFD simulations is far more than that for steady-state simulations.Marsis et al. [111] performed transient two-phase CFD simulations on eight multi-vane ESP designs.The predicted pump performance was confirmed by experimental results.The final design was achieved with the stage pressure increased by 4% for single-phase flow and 23% for two-phase flow at inlet GVF = 20% by optimizing the meridional profile and number of blades.
Yu et al. [112,113] conducted unsteady numerical simulation on gas-liquid flow in a multiphase centrifugal pump.Considering multiple interfacial momentum transfer components, such as drag force, lift force, virtual mass force (VMF) and turbulent dispersion force (TDF), they concluded that the two-fluid multiphase model was able to capture the transport process more accurately than the homogeneous model.Compared to turbulent dispersion force, the drag force plays a more dominant role, as shown in Figure 15a.In Figure 15, z is in the streamwise direction, while L is the maximum streamline as the fluid flows from the intake to discharge of the computational domain.As bubble size (D b ) may vary in a rotating centrifugal pump, its effect on the local drag force is limited, which can be verified in Figure 15b.
optimizing the meridional profile and number of blades.
Yu et al. [112,113] conducted unsteady numerical simulation on gas-liquid flow in a multiphase centrifugal pump.Considering multiple interfacial momentum transfer components, such as drag force, lift force, virtual mass force (VMF) and turbulent dispersion force (TDF), they concluded that the two-fluid multiphase model was able to capture the transport process more accurately than the homogeneous model.Compared to turbulent dispersion force, the drag force plays a more dominant role, as shown in Figure 15a.In Figure 15, z is in the streamwise direction, while L is the maximum streamline as the fluid flows from the intake to discharge of the computational domain.As bubble size (Db) may vary in a rotating centrifugal pump, its effect on the local drag force is limited, which can be verified in Figure 15b.pressure increased by 4% for single-phase flow and 23% for two-phase flow at inlet GVF = 20% by optimizing the meridional profile and number of blades.Yu et al. [112,113] conducted unsteady numerical simulation on gas-liquid flow in a multiphase centrifugal pump.Considering multiple interfacial momentum transfer components, such as drag force, lift force, virtual mass force (VMF) and turbulent dispersion force (TDF), they concluded that the two-fluid multiphase model was able to capture the transport process more accurately than the homogeneous model.Compared to turbulent dispersion force, the drag force plays a more dominant role, as shown in Figure 15a.In Figure 15, z is in the streamwise direction, while L is the maximum streamline as the fluid flows from the intake to discharge of the computational domain.As bubble size (Db) may vary in a rotating centrifugal pump, its effect on the local drag force is limited, which can be verified in Figure 15b.unsteady CFD simulation results of multiphase flow patterns in a three-stage centrifugal pump with the corresponding visualization experiments.Good agreement on the positions and shapes of the gas pocket was observed.Ye et al. [117] combined the advanced transient CFD multiphase simulation with finite-element-analysis (FEA) to optimize the 3D metal printing process of hybrid stage prototype and ESPs for high-gas wells.A design-validation tool was developed and a prototype ESP was manufactured which can pump up to 75% GVF gas/liquid mixture without gas locking.

Mechanistic Modeling
Unlike multiphase pipe flow, whose flow structures and characteristics are easier to visualize and quantify, the mechanistic modeling of ESP performance under gas/liquid flow is more difficult due to its highly twisted flow channel and compact assembly.At present, only limited studies are available in the literature related to this topic, such as Barrios [39], Gamboa [40] and Zhu [1].Similar to the mechanistic model of multiphase pipe flow, the flow patterns in ESP two-phase flow should be identified first, based on which the corresponding flow models can be developed.
Based on Sachdeva [82] and Sun [12] momentum equations in Equation ( 24) along the streamline inside a rotating centrifugal pump, Zhu [1] derived the combined momentum equation of gas and liquid phases by assuming the Taylor bubble flow pattern prevailing in rotating ESP two-phase flow: where the subscripts L, F, S, C, LF, and I denote liquid, film, slug, gas core, liquid film and interface, respectively.v T is the translational velocity.Compared to the combined momentum equation for slug flow in pipe [115,116,[118][119][120], the difference in Equation ( 24) is the body force term.All the velocities are the relative velocities to the ESP channel.By applying the flow pattern and necessary assumptions, Equation ( 24) can be reduced to the simplified flow model, under which the flow characteristics and pump performance can be solved.

Flow Pattern Map and Transition Boundary
As aforementioned, the visualization of flow patterns inside ESP with gas/liquid flow either involves difficult experimentation to install transparent glasses so as to enable the visualization of internal flow fields [39,70] or requires the complicated and expensive non-intrusive instrumentation [66].However, the identification of flow patterns and transition boundaries inside a rotating ESP impeller is still not well understood.Gamboa and Prado [72] proposed an alternative approach to recognize the flow patterns of ESP multiphase flow by relating the inflection points on pump performance curves.Figure 16 below illustrates the identification process schematically.
Gamboa and Prado [72] ascribed the inflections on each performance curve to the flow pattern transitions inside an ESP with gas/liquid flow.As the authors pointed out that such analogy can merely achieve a rough accordance, it did offer significant convenience to obtain the transition boundary among flow patterns.In Figure 16, the horizontal axis is normalized liquid flow rate, and the vertical axis is the measured pump pressure increment.Compared to single-phase performance curve with no inflections, it is found that the two-phase performance curve suffers from significant pressure degradation with inflections at varying loci as q ld drops.Different markers representing the inflection points on each mapping test performance curve correspond to different flow pattern transitions.The white circles denote the transition boundary between dispersed bubble flow and bubbly flow; green diamond shows the transition from bubbly flow to intermittent flow; green square corresponds to the transition from intermittent flow to segregated flow.Gamboa and Prado [72] ascribed the inflections on each performance curve to the flow pattern transitions inside an ESP with gas/liquid flow.As the authors pointed out that such analogy can merely achieve a rough accordance, it did offer significant convenience to obtain the transition boundary among flow patterns.In Figure 16, the horizontal axis is normalized liquid flow rate, and the vertical axis is the measured pump pressure increment.Compared to single-phase performance curve with no inflections, it is found that the two-phase performance curve suffers from significant pressure degradation with inflections at varying loci as qld drops.Different markers representing the inflection points on each mapping test performance curve correspond to different flow pattern transitions.The white circles denote the transition boundary between dispersed bubble flow and bubbly flow; green diamond shows the transition from bubbly flow to intermittent flow; green square corresponds to the transition from intermittent flow to segregated flow.

Transition from Dispersed Bubble Flow to Bubbly Flow (1st Transition Boundary)
As discussed above, the dispersed bubble flow is considered as no-slip flow.On the contrary, bubbly flow, featured by phase slippage, corresponds to the uneven dispersion of gas bubbles in liquids.For a rotating ESP, the initiation of pressure surging is closely related to the 1st transition boundary [44,72].Due to the lack of a valid mechanistic model to predict the initiation of pressure surging inside ESPs, the empirical correlations have been developed from multi-variable regression of experimental data.Table 9 below summarizes several existing pertinent studies.
Figure 17 shows the comparison of surging initiation models in ESPs.The same flow conditions with rotational speed N = 3500 rpm and separator pressure at Psep = 1034 kpa are applied to all calculations.As can be seen, the predictions of ESP surging initiation from existing models vary significantly.Turpin et al. [78] and Cirilo [33] suggested that pressure surging should be independent on liquid flow rates (QL).Other studies [37,38,72] claimed that the critical gas volumetric fraction (λC) at which ESP pressure surging initiates is in a close relationship with QL.Duran [37] and Zapata [38] correlations predict λC in monotonic-increase trend with respect to QL. Gamboa and Prado correlated λC and QL with a concave quadratic function, and claimed that this model was validated for ⁄ > 0.2, beyond which the prediction error might occur.
Cirilo [33] 4342 .0 0187 .0 Transition from Dispersed Bubble Flow to Bubbly Flow (1st Transition Boundary) As discussed above, the dispersed bubble flow is considered as no-slip flow.On the contrary, bubbly flow, featured by phase slippage, corresponds to the uneven dispersion of gas bubbles in liquids.For a rotating ESP, the initiation of pressure surging is closely related to the 1st transition boundary [44,72].Due to the lack of a valid mechanistic model to predict the initiation of pressure surging inside ESPs, the empirical correlations have been developed from multi-variable regression of experimental data.Table 9 below summarizes several existing pertinent studies.
Figure 17 shows the comparison of surging initiation models in ESPs.The same flow conditions with rotational speed N = 3500 rpm and separator pressure at P sep = 1034 kpa are applied to all calculations.As can be seen, the predictions of ESP surging initiation from existing models vary significantly.Turpin et al. [78] and Cirilo [33] suggested that pressure surging should be independent on liquid flow rates (Q L ).Other studies [37,38,72] claimed that the critical gas volumetric fraction (λ C ) at which ESP pressure surging initiates is in a close relationship with Q L .Duran [37] and Zapata [38] correlations predict λ C in monotonic-increase trend with respect to Q L .Gamboa and Prado correlated λ C and Q L with a concave quadratic function, and claimed that this model was validated for Q L /Q max > 0.2, beyond which the prediction error might occur.
In Figure 17, the red-asteroid curve is obtained based on the mechanistic model predictions with a dome-shape with the maximum λ C at Q L close to BEP (grey dashed line).At low Q L , turbulent kinetic energy in ESP impeller is low.It is easy for bubbles to coalesce and form large gas pocket.The flow pattern transition and pressure surging will occur earlier at smaller λ C .At high Q L , the hydraulic head of ESP is close to zero due to open flow conditions.Thus, the turbulent kinetic energy of fluids is also low according to Padron [121] postulation ε = k∆PQ L /(ρV), where k is a constant value obtained from experiments.Then, bubbles are more likely to collide and coalesce to generate bigger ones, leading to further decrease of pump boosting pressure if the gas entrainment rate increases.Thus, the dome-shape of λ C versus Q L is reasonable.
Zhu et al. [44]    17, the red-asteroid curve is obtained based on the mechanistic model predictions with a dome-shape with the maximum λC at QL close to BEP (grey dashed line).At low QL, turbulent kinetic energy in ESP impeller is low.It is easy for bubbles to coalesce and form large gas pocket.The flow pattern transition and pressure surging will occur earlier at smaller λC.At high QL, the hydraulic head of ESP is close to zero due to open flow conditions.Thus, the turbulent kinetic energy of fluids is also low according to Padron [121] , where k is a constant value obtained from experiments.Then, bubbles are more likely to collide and coalesce to generate bigger ones, leading to further decrease of pump boosting pressure if the gas entrainment rate increases.Thus, the dome-shape of λC versus QL is reasonable.For this flow pattern transition, it associates with the severe collision, breakup and coalescence of free bubbles, which triggers the formation of large gas pocket and unstable flow characteristics.Unfortunately, quite few modeling studies are available in literature on this topic.Similar to slug flow in pipes, the intermittent flow inside a rotating ESP involves transient flow coupled with complex pump geometry, which brings in great difficulties to the quantitative analysis.Zhu [1] claimed a critical in-situ gas void fraction (α Crit ) based on the maximum packing theory in a 3D cubic that initiates the transition of bubbly flow to intermittent flow.In multiphase pipe flow, α Crit = 0.25 is mostly used as the critical value.To account for the rotation effect in ESP, a new value was proposed by Zhu [1] as below: where N REF is the rotational speed at best efficiency point, which is 3500 rpm for most ESPs.If N = 0 rpm (no-rotation), Equation ( 25) is reduced to α Crit = 0.25, the same value as used in prediction of the transition boundary from bubbly flow to slug flow in two-phase pipe flow.If N = +∞, the equation reduces to α Crit = 0.52, corresponding the maximum packing of bubbles.The index n (n > 0, usually bounded by 4.0) is an empirical constant determined by experimental data [1].When the transition from slug flow to annular flow occurs, the momentum exchange term can be neglected.Given the superficial gas velocity v SG , and making a guess for the superficial liquid velocity v SL , the critical liquid holdup of the film can be obtained by iterations.
The combined momentum balance equation (Equation ( 26)) can be reduced to: With the calculated H LF (liquid holdup in film), v F (film velocity), H LC (liquid holdup in gas core), and v C (gas core velocity), Equation ( 26) can be solved to get a new v F .Several iterations are required to reach convergence.Finally, the superficial gas/liquid velocity can be calculated.

Flow Model
With the analytical understanding of flow patterns prevailing inside a rotating ESP and the transition boundaries, the assumptions that simplify the governing equation (Equation ( 26)) can be made.Then the flow characteristics of each flow pattern can be solved by iterations.However, due to the limited studies on ESP flow models with gas/liquid flow in literature, the material is insufficient to conduct the corresponding review.Thus, this part is not discussed in detail.

Mechanistic Model Calculation
For the mechanistic models of two-phase flow in a rotating ESP discussed above, a calculation flow chart is given below in Figure 19.
From our previous studies [1,90,122], the inputs of the mechanistic model contain pump geometrical parameters such as radius, blade angles, channel volume and cross sectional area etc., as

Flow Model
With the analytical understanding of flow patterns prevailing inside a rotating ESP and the transition boundaries, the assumptions that simplify the governing equation (Equation ( 26)) can be made.Then the flow characteristics of each flow pattern can be solved by iterations.However, due to the limited studies on ESP flow models with gas/liquid flow in literature, the material is insufficient to conduct the corresponding review.Thus, this part is not discussed in detail.

Mechanistic Model Calculation
For the mechanistic models of two-phase flow in a rotating ESP discussed above, a calculation flow chart is given below in Figure 19.

Closure Relationship
In mechanistic or numerical models, the closure relationships are needed on top of the conservation equations.The closure relationships in modeling centrifugal pump two-phase performance include bubble size, drag force coefficient, in-situ gas void fractions, friction factors etc.

Bubble Size
The bubble size prediction is a critical closure relationship in mechanistic modeling of ESP performance under gassy conditions.However, a generally validated mechanistic model for predicting bubble size in centrifugal pumps is not available.Several proposed bubble size models for centrifugal pump are either empirical or semi-empirical.They were verified with specific pumps and air/water as working fluids [31,32,39,40].The generality of these models is questionable.The bubble characteristics inside the rotating impeller are affected by many factors, including fluid properties (density, viscosity, surface tension), pump geometry and operating parameters (rotational speed, flow rates).
By photographing the bubble dispersion in a pump with a transparent Plexiglas casing, Murakami and Minemura [31,32] (27) Although Equation ( 27) is based on rotational speed (N) and inlet GVF, the empirical nature limited its applicability (λ < 8%, 1000 < N < 1500 rpm).Another model of bubble size in a centrifugal pump was proposed by Estevam [64] based on analogy to gas-liquid two-phase pipe flow.Applying Hinze [123] theory for droplet breakup mechanism in turbulent flow to bubble size prediction in centrifugal pump, Estevam obtained the following equation to calculate the maximum dispersed bubble size: From our previous studies [1,90,122], the inputs of the mechanistic model contain pump geometrical parameters such as radius, blade angles, channel volume and cross sectional area etc., as well as the fluid properties, e.g., gas/liquid density, viscosities, surface tension.The ESP operation conditions (rotational speed, flow rates etc.) are also required.The flow pattern transition boundaries are similar to two-phase pipe flow.The flow map can be divided into four different flow regimes, namely dispersed bubble flow, bubbly flow, intermittent flow, and segregated flow.Based on the limited visualization results, the intermittent and segregated flows in a rotating ESP are slug flow and co-current annular flow, respectively.
With the flow pattern determined, the respective flow model can be called to solve the governing equations to obtain the flow structure and other hydraulic parameters, among which the in-situ gas void fraction (α G ) inside the rotating impeller is the most important parameter to estimate ESP boosting pressure under gassy flow conditions.

Closure Relationship
In mechanistic or numerical models, the closure relationships are needed on top of the conservation equations.The closure relationships in modeling centrifugal pump two-phase performance include bubble size, drag force coefficient, in-situ gas void fractions, friction factors etc.

Bubble Size
The bubble size prediction is a critical closure relationship in mechanistic modeling of ESP performance under gassy conditions.However, a generally validated mechanistic model for predicting bubble size in centrifugal pumps is not available.Several proposed bubble size models for centrifugal pump are either empirical or semi-empirical.They were verified with specific pumps and air/water as working fluids [31,32,39,40].The generality of these models is questionable.The bubble characteristics inside the rotating impeller are affected by many factors, including fluid properties (density, viscosity, surface tension), pump geometry and operating parameters (rotational speed, flow rates).
By photographing the bubble dispersion in a pump with a transparent Plexiglas casing, Murakami and Minemura [31,32] correlated the observed bubble sizes with a linear relationship of Sauter Mean Diameter (d 32 ) versus GVF (Equation ( 29)): Although Equation ( 27) is based on rotational speed (N) and inlet GVF, the empirical nature limited its applicability (λ < 8%, 1000 < N < 1500 rpm).Another model of bubble size in a centrifugal pump was proposed by Estevam [64] based on analogy to gas-liquid two-phase pipe flow.Applying Hinze [123] theory for droplet breakup mechanism in turbulent flow to bubble size prediction in centrifugal pump, Estevam obtained the following equation to calculate the maximum dispersed bubble size: where the constant of 1.17 is a parameter accounting for the curvature of impeller flow passages; σ, ρ l and D H are interfacial tension, liquid density and hydraulic diameter, respectively; the friction factor (f β,ω ) is determined by analogy to fluid flow in pipeline and considering rotational speed effect.Following a similar methodology of modeling bubble size inside a centrifugal pump, Barrios [39] proposed a bubble size prediction model based on the visualization experimental data inside ESP taken by high-speed camera.Barrios pointed out that the experimental results were necessary to determine the relationship of the maximum bubble size, the critical bubble size and inlet GVF.The critical Weber number (We crit ) is a parameter dominating gas bubble breakup.Hence, Equation (30) explicitly relates the bubble size with the rotational speed and liquid properties: where r 1 is the impeller radius.Gamboa [40] employed the Levich [124] model for maximum stable bubble size in pipe flow and proposed an alternative way of modeling bubble size inside an ESP.Gamboa improved the bubble size model by introducing dispersed gas phase density and We crit based on Kouba [125] droplet breakup studies: where ρ g is in-situ gas density, v is liquid kinetic viscosity and Ω, D are impeller rotational speed and diameter, respectively.A new correlation for computing the representative bubble sizes inside or rotating ESP impeller was recently proposed by Zhu and Zhang [105], as given by Equation (31).The model is based on CFD simulations of pump pressure increment under gassy flows, where the bubble size plays a dominant role to make the numerically simulated ESP heads comparable to the corresponding experimental data.The proposed correlation was verified by further comparisons of the CFD simulation results with the incorporation of Equation ( 33) against the experimental pump pressure increment, which shows a good match: where d 32 is the Sauter mean diameter; λ G is the inlet gas volumetric fraction; subscripts of c and d denote continuous and dispersed phases, respectively.∆P is the pressure increment of ESP under single-phase flow; V is the volume of the entire impeller flow domain. Energies The comparison of existing bubble size models in literature is shown in Figure 20, where the horizontal axis corresponds to λ G and the vertical axis shows the calculated bubble diameters.As can be seen, the predicted bubble sizes in the rotating ESP impeller by different models vary significantly.The reason might be due to the derivation of prediction model that is based on either empirical constant or partial regression of experimental data.The complex physical reality of multiphase flow through a rotating ESP impeller further contributes to the difficulty of reaching a generally validated prediction model of bubble size.Meanwhile, the mutual validation among different models seems to be unavailable due to the lack of relevant experimental measurement of in-situ bubble sizes as well as geometrical details (channel volume etc.).A similar trend of d b versus λ G is used by Barrios [39] and Gamboa [40] to correlate the experimental measurement of bubble diameters.However, the linear relationship of d b versus λ G is applied by Murakami and Minemura [31,32], and Zhu and Zhang [105] studies, where a good agreement between the model predicted d b against either visualization experiments or CFD simulation results can be reached.Although several models are available in literature for bubble size prediction in a centrifugal pump with rotating turbulent flow, their empirical/semi-empirical natures limited the range of applications.As discussed by Gamboa [40], the challenge in modeling bubble size inside a rotating ESP is the mechanism that dominates bubble formation, including coalescence and breakup.Therefore, investigation is needed to better understand the bubble dispersion mechanism and develop better model for bubble size prediction.

Drag Coefficient
In gas-liquid two-phase flow, the drag force is the interfacial momentum transfer due to velocity difference between gas and liquid phases [102]: where U and v are velocities of liquid and gas phases, respectively.The drag coefficient (CD) for bubbles in no-rotating flow fields without shearing was studied by Schiller and Naumann [126], Clift et al. [127], Ishii and Zuber [128], Mei et al. [129]   Although several models are available in literature for bubble size prediction in a centrifugal pump with rotating turbulent flow, their empirical/semi-empirical natures limited the range of applications.As discussed by Gamboa [40], the challenge in modeling bubble size inside a rotating ESP is the mechanism that dominates bubble formation, including coalescence and breakup.Therefore, investigation is needed to better understand the bubble dispersion mechanism and develop better model for bubble size prediction.

Drag Coefficient
In gas-liquid two-phase flow, the drag force is the interfacial momentum transfer due to velocity difference between gas and liquid phases [102]: where U and v are velocities of liquid and gas phases, respectively.
Ishii and Zuber incorporated a correlation term (1 + α d ) −γ to account for effects of bubble volumetric fraction and flow regime, where γ is an empirical constant determined by fluid properties and particle shapes.Mei et al. [129] studied the behavior of clean bubbles (no contaminations or surfactants involved) in a uniform flow and proposed an empirical drag C D for 0.1 < Re < 100: For shear-induced flow, the viscous drag force exerted on bubbles is increased by broadening the near wake [131].Thus, the dimensionless shear rate, Strouhal number (Sr) as an indicator of shear strength, was employed by Legendre and Magnaudet [132] to calculate C D in shear flow for moderate-to-large Re (≥50): and: where C D,0 is the drag coefficient calculated without shear effect [131].Rastello et al. [133,134] revised Equation (36) to get a better fitting of their experimental data for low, moderate and high Reynolds numbers with a broad range of Sr: Barrios [39] measured bubble sizes inside a single-stage ESP with a visualization experimental system and calculated the drag coefficients on stagnant bubbles in rotating flow field.Then, the drag coefficients were correlated by: where Y and f (Re,Y) are given by: and: Figure 21 shows the comparison of drag coefficients predicted by the empirical models discussed above.As can be seen, the evident discrepancy is observed under relatively low Reynolds number with/without considering the rotating effects in flow fields.At higher Reynolds number, the similar drag coefficients, which are close to 0.44, can be obtained from the existing empirical models.The drag coefficients predicted by Clift et al. [127] and Mei et al. [128] is close to each other due to neglecting rotation flow.However, when accounting for fluid rotating effect, the model predictions of C D vary significantly.Barrios [39] correlation outputs C D that is several times larger than that calculated by Legendre and Magnaudet [132] model.A possible explanation is that Barrios model is specifically on basis of the bubble flow inside a rotating ESP impeller, while the Legendre and Magnaudet model is regressed from a single bubble flowing inside a rotating cylinder.

In-Situ Gas Void Fraction (αG)
The in-situ gas void fraction (αG) is an important variable in two-phase flow related to the velocity slippage between two immiscible phases and the gas accumulation.Moreover, αG is needed to determine the in-situ gas/liquid mixture density by = + (1 − ) , which in turn affects the pump boosting pressure significantly according to ∆ = .ρM is the mixture density, and HM is the hydraulic head of ESP with gas/liquid flow.However, due to the complicated pump geometries and fluid flow dynamics in the ESP impeller, it is challenging to develop a mechanistic model to predict αG with general validity.Very few studies on mechanistic modeling of local gas void fraction in ESPs can be found in literature so far.
The simplest model for predicting αG in multiphase centrifugal pump flows is the homogeneous model, which assumes no slippage between gas and liquid phases: The homogenous model is valid for very low inlet GVFs when the slippage between gas and liquid is minimal.Errors will result from applying the homogeneous model to flow conditions with high inlet GVFs when the slippage is not negligible.Accounting for the phase slippage, several empirical correlations were proposed by Chisely [62], Estevam [64], Zapata [38] and Pineda et al. [114], among others.
Chisely studied loss of coolant accident with a volute-type centrifugal pump in the nuclear industry.The pressure distribution and flow regimes were determined with the experimental measurement and high-speed photographing.A model for predicting the centrifugal pump pressure increment under two-phase flow conditions was proposed.To make this model solvable, the in-situ αG was correlated as: where χ is the gas quality, μG and μL are the gas and liquid viscosities.
Zapata experimentally studied the rotational speed effect on ESP two-phase performance.Using least-square regression, a new correlation of αG was presented to predict the pump boosting pressure.The empirical correlation as function of gas and liquid flow rates and the rotational speed is given as:

In-Situ Gas Void Fraction (α G )
The in-situ gas void fraction (α G ) is an important variable in two-phase flow related to the velocity slippage between two immiscible phases and the gas accumulation.Moreover, α G is needed to determine the in-situ gas/liquid mixture density by ρ M = α G ρ G + (1 − α G )ρ L , which in turn affects the pump boosting pressure significantly according to ∆P M = ρ M gH M .ρ M is the mixture density, and H M is the hydraulic head of ESP with gas/liquid flow.However, due to the complicated pump geometries and fluid flow dynamics in the ESP impeller, it is challenging to develop a mechanistic model to predict α G with general validity.Very few studies on mechanistic modeling of local gas void fraction in ESPs can be found in literature so far.
The simplest model for predicting α G in multiphase centrifugal pump flows is the homogeneous model, which assumes no slippage between gas and liquid phases: The homogenous model is valid for very low inlet GVFs when the slippage between gas and liquid is minimal.Errors will result from applying the homogeneous model to flow conditions with high inlet GVFs when the slippage is not negligible.Accounting for the phase slippage, several empirical correlations were proposed by Chisely [62], Estevam [64], Zapata [38] and Pineda et al. [114], among others.
Chisely studied loss of coolant accident with a volute-type centrifugal pump in the nuclear industry.The pressure distribution and flow regimes were determined with the experimental measurement and high-speed photographing.A model for predicting the centrifugal pump pressure increment under two-phase flow conditions was proposed.To make this model solvable, the in-situ α G was correlated as: where χ is the gas quality, µ G and µ L are the gas and liquid viscosities.Zapata experimentally studied the rotational speed effect on ESP two-phase performance.Using least-square regression, a new correlation of α G was presented to predict the pump boosting pressure.The empirical correlation as function of gas and liquid flow rates and the rotational speed is given as: α G = q G q max (q L /q max ) −1.277−0.034N where q G and q L are the gas and liquid flow rates, N d N d is normalized rotational speed, and q max is the maximum single-phase liquid flow rate.Pineda et al. [114] proposed an empirical correlation based on the Lockhart-Martinelli type parameter X tt by non-linear regression of CFD simulated values of in-situ α G : α G = 7.119X −0.8778 tt − 0.002138 (44) where: Zhu and Zhang [90] proposed a mechanistic model to compute the in-situ gas void fraction.The model prediction was validated by the 3D CFD simulation.A better match over the existing empirical models was achieved.The model is based on the balance of centrifugal buoyancy force and drag force on a stable gas bubble in ESP, which incorporates the gas-liquid interaction mechanism and is applicable to a wide flow conditions.The prediction model is given as, where: V SR is the slippage velocity in radial direction, R I is the representative impeller radius, Z I is the impeller number, T B is the blade thickness, Y I is the impeller channel height, Q and Q LK is the liquid flow rate and its corresponding leakage flow rate.The drag coefficient is calculated by C D,0 1 + 0.55Sr 2 Re > 50 C D,0 1 + 0.3Sr 2.5 Re ≤ 50 (48) where Sr is Strouhal number and defined by Sr = d B Ω/|U − v|.U and v are phase velocities.C D,0 is the drag coefficient proposed by Clift et al. [127] without consideration of shear effect.Figure 22 shows the comparison of existing models in predicting in-situ gas void fraction in a rotating pump impeller.As can be seen, Zhu and Zhang [90] model predicts in-situ α G with an error less than 25%, and most predictions are bounded by 10% error line.Nevertheless, α G predicted by previous empirical correlations deviates from CFD simulated α G increasingly as inlet λ G increases.Since most of the correlations were based on homogeneous assumption rather than velocity slippage between gas and liquid phases, the predictions from empirical correlations are acceptable only at low λ G .For relatively higher λ G , the phase slippage dominates bubble behaviors, which results in the failure of empirical correlations.

Conclusions and Future Work
In this paper, the experimental studies and modeling approaches for investigating gas/liquid two-phase flow inside a rotating electrical submersible pump have been reviewed in detail.Different empirical and mechanistic models for predicting ESP boosting pressure under gassy flows were discussed.The state-of-art CFD-based numerical simulation as a complemental method for ESP twophase performance prediction was elaborated.
The review of the published experimental studies for gas-liquid flow in centrifugal pumps has highlighted the complications behind the accurate prediction of two-phase performance inside a rotating ESP.The turbulent flow with transient behaviors make it difficult to visualize the internal flow structures and flow patterns.This is not helped by the twisted blade/vane geometries and multistage assembly.Due to these observation difficulties, the multiphase flow mechanisms that dominate the ESP head degradation were not well understood.This in turn delayed the development of mechanistic models.Although CFD simulation serves as an alternative way to study multiphase flow, its validity and reliability, especially coupled with the strong shearing effect and phase interactions are questioned.
The ESP performance under gassy flow is affected by many factors.There is room for further improvements on the mechanistic modeling and prediction:

Conclusions and Future Work
In this paper, the experimental studies and modeling approaches for investigating gas/liquid two-phase flow inside a rotating electrical submersible pump have been reviewed in detail.Different empirical and mechanistic models for predicting ESP boosting pressure under gassy flows were discussed.The state-of-art CFD-based numerical simulation as a complemental method for ESP two-phase performance prediction was elaborated.
The review of the published experimental studies for gas-liquid flow in centrifugal pumps has highlighted the complications behind the accurate prediction of two-phase performance inside a rotating ESP.The turbulent flow with transient behaviors make it difficult to visualize the internal flow structures and flow patterns.This is not helped by the twisted blade/vane geometries and multistage assembly.Due to these observation difficulties, the multiphase flow mechanisms that dominate the ESP head degradation were not well understood.This in turn delayed the development of mechanistic models.Although CFD simulation serves as an alternative way to study multiphase flow, its validity and reliability, especially coupled with the strong shearing effect and phase interactions are questioned.
The ESP performance under gassy flow is affected by many factors.There is room for further improvements on the mechanistic modeling and prediction:

•
Experimental visualization of internal two-phase flow structures in a rotating ESP is inadequate, which restricts the flow pattern assessment and the formulations of governing equations.Moreover, the measurement of in-situ gas phase is insufficient, resulting in further difficulty to characterize pump performance quantitatively.Future experimental work should take these challenges into consideration, thus providing in-depth knowledge of multiphase flow mechanisms in ESPs.

•
The computational fluid dynamics techniques need to be further developed in order to capture the complex behaviors and movement of gas bubbles inside ESPs and thus better understand the interactions between the gas and liquid phases in different flow patterns.

•
The mechanistic models should focus more on the physics of the two-phase flow in rotating ESPs to propose better closure relationships.From the one-dimensional two-fluid model, the simplified mass and momentum conservation equations are solvable with the help of the improved closure relationships in future, ending up with a more accurate and reliable mechanistic model to predict ESP two-phase flow performance.

Figure 2 .
Figure 2. Typical ESP pump performance curves (figure courtesy of Wood Group ESP, Inc., Oklahoma City, OK, USA).

Figure 2 .
Figure 2. Typical ESP pump performance curves (figure courtesy of Wood Group ESP, Inc., Oklahoma City, OK, USA).
relative velocity W and peripheral velocity U, which is calculated by →

Energies 2018 ,
11, 180 5 of 41 flow rates.Shock loss takes place when the liquid flow rate differs from the designed flow rate.Tables 1-5 below discuss each pressure loss mechanism by presenting existing models that are available in literature.(a) (b)

Figure 6 .
Figure 6.Comparison of ESP two-phase pressure increment with/without surfactant injection [44].(a) Constant liquid flow rate testing; (b) constant gas flow rate testing.

Figure 7 .
Figure 7.Typical flow patterns and transition boundaries in vertical pipe, (a) main observed flow patterns [54]; (b) transition boundaries and flow pattern map [52].

Figure 6 .
Figure 6.Comparison of ESP two-phase pressure increment with/without surfactant injection [44].(a) Constant liquid flow rate testing; (b) constant gas flow rate testing.

Figure 6 .
Figure 6.Comparison of ESP two-phase pressure increment with/without surfactant injection [44].(a) Constant liquid flow rate testing; (b) constant gas flow rate testing.

Figure 7 .
Figure 7.Typical flow patterns and transition boundaries in vertical pipe, (a) main observed flow patterns [54]; (b) transition boundaries and flow pattern map [52].Figure 7. Typical flow patterns and transition boundaries in vertical pipe, (a) main observed flow patterns [54]; (b) transition boundaries and flow pattern map [52].

Figure 7 .
Figure 7.Typical flow patterns and transition boundaries in vertical pipe, (a) main observed flow patterns [54]; (b) transition boundaries and flow pattern map [52].Figure 7. Typical flow patterns and transition boundaries in vertical pipe, (a) main observed flow patterns [54]; (b) transition boundaries and flow pattern map [52].

Energies 2018 ,
11, 180 12 of 41 Later visualization experimental work, conducted by Gamboa [40] presented similar observations as Estevam's study, based on which the first flow pattern map with transition boundaries between different flow patterns inside an ESP was proposed.His work offers significant insight regarding how to characterize the flow patterns in an ESP quantitatively.Among previous visualization experiments, two recent studies from Schäfer et al. [68] and Neumann et al.
0%.As confirmed by Schäfer et al. observation, such change is due to the rapid flow pattern change from bubbly flow to intermittent flow.Energies 2018, 11, 180 15 of 41 0%.As confirmed by Schäfer et al. observation, such change is due to the rapid flow pattern change from bubbly flow to intermittent flow.
Based on Sachdeva et al. and Minemura et al. one-dimensional two-fluid models, Sun[12,22] developed a new two-phase model including a set of one-dimensional mass and momentum balance equations for predicting ESP performance.He also improved analytical models for wall frictional losses and shock loss, as well as new correlations for the drag coefficient.The general momentum balance equation along the streamlines is given by: dp dr streamline = −ρ p W p dW p dr + ρ p Ω 2 r + dp ds f ,p

Figure 14 .Figure 15 .
Figure 14.CFD predicted gas void fraction distribution and the corresponding pressure contour [90], (a) gas void fraction color map; (b) pressure contour on half impeller span.

Figure 14 .
Figure 14.CFD predicted gas void fraction distribution and the corresponding pressure contour [90], (a) gas void fraction color map; (b) pressure contour on half impeller span.

Figure 14 .Figure 15 .
Figure 14.CFD predicted gas void fraction distribution and the corresponding pressure contour [90], (a) gas void fraction color map; (b) pressure contour on half impeller span.

Figure 15 .
Figure 15.Magnitude analysis of interface forces by CFD simulations [112], (a) magnitude ratios of non-drag forces to drag force; (b) bubble size effect on drag force magnitude.Figure 15.Magnitude analysis of interface forces by CFD simulations [112], (a) magnitude ratios of non-drag forces to drag force; (b) bubble size effect on drag force magnitude.
shows a typical calculation result for predicting the transition boundary, where the regimes I, II, III and IV correspond to the dispersed bubble flow, bubbly flow, intermittent flow and segregated flow, respectively.The legend labels of DB, BF, SF, and AF denote dispersed bubble flow, bubbly flow, slug flow and annual flow.The colored markers represent the flow patterns from experimental results.As can be seen, the flow patterns predicted by mechanistic model are comparable to that detected from experimental performance curves.Compared to Gamboa [40] visualization experiments, similar flow pattern transition boundaries are well predicted.At low gas flow rate (Q G ), the dispersed bubble flow encompassed by the blue line or bubbly flow bounded by the red line prevails.With the increase of gas flow rate or decrease of liquid flow rate (Q L ), slug flow takes place with the boundary of the black dot-dash line.The segregated flow (regime IV) corresponds to extremely low Q L .Energies 2018, 11, 180 26 of 41

Figure 18 .
Figure 18.Validation of mechanistic model for predicting the flow pattern transition boundary in ESP under gas/liquid flow [1].

Figure 18 .
Figure 18.Validation of mechanistic model for predicting the flow pattern transition boundary in ESP under gas/liquid flow [1].

Figure 19 .
Figure 19.Flow chart for gas-liquid flow calculation inside a rotating ESP impeller.
among others.Schiller and Naumann proposed an empirical correlation of CD for 0.1 < Re < 800.Clift et al. developed a more accurate expression which is valid for higher Re up to 3 × 10 5[130]: incorporated a correlation term (1 + αd) −γ to account for effects of bubble

Figure 20 .
Figure 20.Comparison of bubble size prediction models.

Table 1 .
Friction loss models in literature.

Table 2 .
Shock loss models in literature.

Table 1 .
Friction loss models in literature.

Table 2 .
Shock loss models in literature.

Table 3 .
Leakage loss models in literature.

Table 4 .
Recirculation loss models in literature.

Table 6 .
Disk loss models in literature.

Table 7 .
Experimental studies on ESP two-phase performance conducted at TUALP.

Table 8 .
Experimental studies on two-phase flow patterns inside a centrifugal pump.
[60]32,56]bservationMurakami and Minemura[31,32,56]Transparent casing 1st study that associated pump perfor-mance with gas-liquid flow pattern inside an impellerPatel & Runstadler[57]Transparent casing Small bubble flow and stationary large bubble coinciding with significant reduction of pump head Sekoguchi et al.[58]Transparent casing electric resistivity probeThe evolvement of dispersed bubble to slug flow and a large gas pocket Kim et al.[59]Transparent casing Phase slippage leads to bubble agglomeration Sato et al.[60]Transparent casing Bubble flow to separated flow as gas flow rate increases [39] pattern map for GC-6100 ESP at 2400 rpm and 150 psig[39].
d N d 0.598 + 0.223N d N d , 1 bpd ≈ 1.8942 × 10 −6 m 3 /s psi pound per square inch, 1 psi ≈ 6894.76 pa psia pound per square inch for absolute pressure psig pound per square inch for gauge pressure