Analysis of bubble distribution in a multiphase rotodynamic pump

ABSTRACT Recent researches show that the pressure increment of a multiphase pump is affected by bubble size and distribution. In order to study the bubble distribution characteristics in such pumps, a novel approach describing the variable bubble size in the pump is proposed. The bubble number density equation, which has taken into account the phenomena of break-up and coalescence, is introduced into the flow simulation, and the drag coefficient is revised because of the interaction of multiple bubbles. The reliability of the approach is verified by comparison with the experiment. It was established that the bubbles move to the impeller hub due to the difference in centrifugal force between gas and liquid. Despite the high collision rate near the hub, bubble size changes little with the stirring action of the impeller. The mixture flows in a disorderly way and the bubble diameter increases due to the rotor–stator interaction. Owing to the increasing flow area in the diffuser, bubbles move to the mainstream region, and bubble size reaches its maximum owing to the flow separation near the hub. The distribution of bubbles is also analyzed under a different inlet gas volume fraction (IGVF) and inlet bubble diameter (d0). Bigger IGVF brings about a higher collision rate of bubbles, while smaller d0 makes the diffusion of bubbles easier.


Introduction
The exploitation of offshore oil and gas has been strengthened by increased energy consumption . Multiphase transportation technology is widely used in the exploitation since it can prolong the life of the oilfield, reduce investment, and promote production efficiency. As the core equipment of this technology, the multiphase rotodynamic pump can convert the rotational kinetic energy to the hydraulic pressure head, and bring about remarkable economic benefits. In order to take measures to improve pump performance, research about flow characteristics should be conducted. Murakami and Minemura (1983) investigated the movement of bubbles in the pump impeller. Their results showed that the trajectory of bubbles deflects from the path of liquid, and the extent of deviation relates to the bubble diameter and liquid flow rate. Barrios (2007), Barrios, Prado, and Kenyery (2009), and Barrios and Prado (2011) found that the enlargement of bubbles corresponds to the pump's poorer performance when the pressure increment in an electrical submersible pump (ESP) was studied. It is obvious that the flow characteristics and pump performance are directly affected by bubble size, but few studies allow for variable bubble size in the pump. Kim et al. (2015) applied the orthogonality design method for the simulation to improve the hydraulic performance of an axial multiphase pump, but the simulation was performed with a constant bubble size. Shi, Zhu, Zhang, Zhang, and Zhao (2018) studied the flow characteristics in a multi-stage multiphase pump based on the Eulerian-Eulerian inhomogeneous model, but the bubble size was also considered as constant. Zhu and Zhang (2017) proposed a multiple-bubble model and employed it in the pressure increment calculation for a three-stage ESP, where a bubble distribution function was provided. However, the key coefficients were determined by the experimental data, which makes the scope of its application limited. Therefore, bubble size prediction is still needed for flow simulation in the multiphase rotodynamic pump.
Methods for predicting bubble size can be found in studies on bubble columns and stirred tanks. The population balance equation was established by Kocamustafaogullari and Ishii (1995), where the number density function was used, and the change mechanism of the bubble size comprehensively described, while also making considerable demands on computational resources. Lehr and Mewes (2001) simplified the equation and successfully applied it in predicting bubble size in a bubble column, whereas most coefficients in the source term of the equation are dependent on the measurement data. Bakker (1992) and Bakker and Van den Akker (1994) replaced the number density function with the bubble number density, so the bubble size could be calculated in a convenient way. On the basis of Bakker's study, Lane, Schwarz, and Evans (2005) proposed the bubble number density equation. In a stirred tank, the equation was introduced into the calculation of the flow and a good prediction on the bubble diameter and bubble distribution was obtained. Similar to the case of a stirred tank, a rotational flow field also exists in a rotodynamic pump, thus the bubble diameter in such pumps could be predicted by using the bubble number density equation.
In this paper, a novel numerical method is presented to analyze the two-phase flow in the pump. First, secondary development technology is used to bring in the bubble number density equation, and the phenomena of break-up and coalescence are taken into account. Then, the numerical method is verified through the comparison between the simulation and the experiment. Based on the above, the distribution characteristics of bubbles in the flow passage are explored, and the influence of IGVF and d 0 are analyzed before conclusions are drawn in the final section.

Governing equations
As a widely recognized flow pattern in the multiphase rotodynamic pump, bubbly flow is assumed in this simulation. Eulerian-Lagrange and Eulerian-Eulerian are the two main frameworks used in the simulation of twophase flow. Compared to the Eulerian-Lagrange method, it is not necessary to solve the balance equations for every single bubble in the Eulerian-Eulerian method, thus the computational resource is saved (Hanafizadeh, Saidi, Karimi, & Zamiri, 2010;Kim et al., 2015;Sokolichin, Eigenberger, Lapin, & Lübert, 1997;Zhang, Li, Vafai, & Zhang, 2018). On the basis of the Eulerian-Eulerian method, the two-fluid model is applied. The continuity equation and momentum equation are given by Continuity equation Where w i is the relative velocity, M i is the interphase force, ω represents the rotational angular velocity, and τ i stands for viscous stress tensor. −2α i ρ i ω and α i ρ i ω 2 r i respectively denote the coriolis force and the centrifugal force.
The turbulence k-model is often used because of its simplicity, strong stability, and high efficiency (Akbarian et al., 2018;Mou, He, Zhao, & Chau, 2017). However, due to the strong rotating flow in the pump, the shear stress transport k-ω model is more suitable since it can give an accurate prediction for the flow with a big area of separation.

Interphase force
For two-phase flow, the drag force is the most important, while the impact of other forces such as lift and added mass are slightly smaller. The turbulent dispersion force should be considered when the spreading of bubbles is affected by eddies. Therefore, the interphase momentum transfer term (Yu, Zhu, & Cao, 2015) is given by Drag force D i , whose magnitude is the largest among these forces, denotes the resistance suffered by bubbles, calculated by (Tabib & Schwarz, 2011). Taking into account the bubble interaction and gas volume fraction effect, the drag model is modified based on the research of Ishii and Zuber (1979), then the following equations are given: The added mass force A i is an inertia force acting on the accelerating bubble, while the lift force L i may generate when liquid flows past the surface of bubbles. These two forces are expressed as Here C A and C L are the coefficients for calculating added mass and lift forces. Generally, C A and C L are both set as 0.5 (Yu, Zhu, Cao, & Liu, 2014).
Turbulent dispersion force T i has a significant effect of "spreading" gas out due to the interaction between bubbles and eddies. For the calculation of this force, the Lopez de Bertodano model is used as follows, where the dispersion coefficient C td is set as 0.1 (Zhang, Yu, Zahid, & Li, 2018).

Calculation of bubble size
As a simplification of the equation proposed by Bakker (1992) and Bakker and Van den Akker (1994), the bubble number density equation can describe the change mechanism of bubble size, shown as Here, φ br and φ co are respectively the break-up rate and the coalescence rate, while n denotes the variable bubble number density, calculated by According to Lane et al. (2005), the break-up of bubbles is related to turbulent fluctuations at a length scale of the bubbles. The critical Weber number, whose value is taken as 1.5 (Lehr & Mewes, 2001), is used as the break-up criterion. The break-up rate can be expressed as The coalescence of bubbles is assumed as a binary process, which consists of two steps: two bubbles collide and liquid film forms between them; then the film ruptures and coalescence occurs. For the two steps, the binary collisions rate and the coalescence efficiency between bubbles are respectively calculated. Finally, the coalescence rate is as follows By transformation, the bubble number density equa tion is expressed as a generic form of transport equation. Then with the secondary development technology of CFX, i.e. user CEL routine and function (ANSYS, Inc., 2013), the equation can be brought into the numerical simulation.

Geometry model
As is shown in Figure 1, a multiphase rotodynamic pump includes four parts. The gas-liquid mixture flows into the pump through the inlet pipe and obtains kinetic energy from the rotational impeller, then the kinetic energy is converted to the pressure energy when the mixture flows through the stationary diffuser. Finally, the gas-liquid mixture flows out of the outlet pipe.

Boundary conditions and convergence settings
Since the mass flow rate of a gas or liquid is relatively easy to measure, the inlet boundary condition is set as bulk mass flow rate. The mixture finally flows into the water tank, so the average static pressure is set for the outlet boundary condition, and the relative pressure is given as zero. No-slip condition is applied for the walls, and the multiple frames of reference are adopted in the calculation. For the data transmission between the impeller and  the diffuser, the frozen rotor model is used (Zhu & Zhang, 2014). Air and water are selected as the medium. For the domain inlet, the IGVF is from 5% to 10% (Zhang, Zhu, Wei, & Ying, 2011), while d 0 is from 0.4 mm to 1.0 mm (Falcimaigne, Brac, Charron, Pagnier, & Vilagines, 2002;Yu et al., 2015). Meanwhile, the rotational speed N is kept as 2950r/min and the total flow rate Q holds the value of 50 m 3 /h.

Verification of the numerical simulation
The test bench is shown as Figure 2. Air and water, which are respectively transported by a compressor and a pump, are uniformly mixed in the buffer tank. After that, the mixture is transported to the water tank through the multiphase rotodynamic pump. Pressure sensors are installed to acquire the pressure increment of the pump, and the high-speed camera is used to record the instantaneous flow details.
The pump head (Vieira, Siqueira, Bueno, Morales, & Estevam, 2015) obtained from the simulation and the experiment are compared in different IGVF when d 0 = 0.4 mm. As is shown in Table 1, the relative errors between the simulation and the experiment are below 5%.
The comparison of flow characteristics between the simulation and the experiment is shown in Figure 3. The distribution of gas volume fraction taken by a high-speed camera is presented in Figure 3 (a), while the overlap of nine isosurfaces of gas volume fraction from 0.1 to 0.9 at an interval of 0.1 is given in Figure 3 (b). From Figure 3 (a), a great value of gas volume fraction is observed in the interaction region. Moreover, low gas volume fraction is found near the suction side of the diffuser blade at the leading edge. It is obvious that the flow characteristics obtained from the simulation agree well with the experiment, and the rationality of the simulation is verified.

Distribution characteristics of bubbles in the flow passage
The gas distribution in gas-liquid flow is usually characterized by the variable gas volume fraction. According to Equation (13), the gas volume fraction is closely related to the bubble number density n and the bubble size d. For example, in a condition of constant d, bigger n results in a bigger gas volume fraction. In this analysis, n and d are directly used to present the distribution characteristics of bubbles, while the distribution of the gas volume fraction can be obtained through the distribution analysis of the other two.
Since some common characteristics of bubble distribution are shared under different conditions of IGVF and d 0 , the movement of bubbles is analyzed in the case of Q = 50m 3 /h, IGVF = 5% and N = 2950r/min.
As is shown in Figure 4, bubbles gradually accumulate near the impeller hub, which leads to an increased collision rate of bubbles. Meanwhile, big bubbles shown in Figure 5 break up into small bubbles in the impeller. These phenomena can be explained by the influence of the rotational impeller. The difference between gas and liquid in centrifugal force leads to the gathering of bubbles near the hub. Despite a high collision rate in the hub, bubble diameter changes little with the stirring action of the impeller, thus d is uniform in the impeller.
The gas volume fraction in the impeller is also presented. As the variation of d in the impeller is not very large under this condition, the value of gas volume fraction depends mainly on the value of n. The gas volume fraction is large near the hub yet small near the shroud of the impeller, this is the same as the case of n.
With the decrease of the cross-sectional area of the flow passage, n reaches a big value ( Figure 6). In addition, the existence of rotor-stator interaction makes the mixture flows disorderly, so the bubble diameter is increased (Figure 4).
It can be observed from Figure 7 that n decreases along the streamwise direction, and that the largest bubble is formed near the diffuser hub. Bubbles move to the mainstream region gradually due to the diffusion of the diffuser structure. Because of the rising pressure and the falling air velocity along the streamwise direction ( Figure  8), flow separation occurred near the diffuser hub, and the collision rate of bubbles increased.
It can be observed in Figure 9 that the bubble diameter near the suction surface (SS) is bigger than that near the pressure surface (PS) in the diffuser, owing to the formation of large vortexes near the SS. As is shown in the distribution of water velocity, the PS in the diffuser is directly flushed by the mixture, and vortexes formed near the SS in the diffuser. Due to the effect of the vortexes, the collision phenomena among bubbles increase and big bubbles appear.

Sensitivity of bubble distribution to IGVF
Two-phase flow characteristics will change under conditions of different IGVF in a multiphase rotodynamic pump (Caridad, Asuaje, Kenyery, Tremante, & Aguillón, 2008;Pirouzpanah, Gudigopuram, & Morrison, 2017). In this work, the distribution of bubbles is analyzed when IGVF is respectively 5%, 10%, and 15%, while d 0 is kept as 0.4 mm. Figure 10 presents the bubble distribution in the inlet pipe at different IGVF. In the condition of IGVF = 15%, n is highest at the pipe inlet while lowest at the pipe outlet. This indicates that d increases with a greatest growth rate in the condition of IGVF = 15%. According to Equation (13), higher IGVF results in larger n but a smaller distance between bubbles; this brings about a bigger interaction and collision rate of bubbles.
From Figures 11 and 12, it can be seen that the big bubble area and the vortex formed near the SS in the diffuser enlarge with an increased IGVF. Higher IGVF means more obvious flow separation, which is characterized by a larger vortex. Therefore, the big bubble area is largest when IGVF = 15%.

Sensitivity of bubble distribution to d 0
The distribution of bubbles is analyzed when d 0 is respectively 0.4 mm, 0.7 mm, and 1.0 mm, while IGVF is kept as 10%. In these cases, both the inlet bubble diameter   and inlet bubble number density are different. To make the analysis reasonable, the dimensionless relative bubble diameter d r and relative bubble number density n r are defined as follows In Figure 13, the relative bubble diameter decreases with an increased d 0 . This means that the collision rate of bubbles is high at a small d 0 . According to Equation (13), small d 0 indicates big n 0 under the condition of a constant IGVF, which will lead to a smaller distance between bubbles and more frequent bubble interaction, as the bubble diameter obviously increases.
In Figure 14, the distribution of relative bubble number density when d 0 is 0.4 mm is uniform at the diffuser outlet, while nonuniform when d 0 is 0.7 mm or 1.0 mm. This means that the diffusion of bubbles is easier at a smaller d 0 . In addition, bubbles tend to accumulate near the hub region with an increased d 0 , and the accumulation area locates mainly at the diffuser inlet. Therefore, for a better operation of the multiphase pump,  the control of bubble diameter to a smaller value is recommended, and the hydraulic design of the hub region should be emphasized.

Conclusions and remarks
With the consideration of changeable d in the simulation, the distributions of n and d in the multiphase rotodynamic pump have been analyzed. Conclusions are as follows: (1) The bubble distribution in the impeller is mainly influenced by the rotation of the impeller. On the one hand, the stirring action caused by the impeller leads to the break-up of bubbles; on the other hand, the difference between gas and liquid in centrifugal force makes bubbles move to the hub. In interaction regions, the bubble diameter increases due to the rotor-stator interaction in this region, and n increases obviously owing to the decreased crosssectional area. Bubbles move to the mainstream region gradually in the diffuser, and d reaches the largest value, since the flow separation occurred near the diffuser hub.
(2) Higher IGVF results in a larger n but a smaller distance between bubbles; this makes the interaction of bubbles in the inlet pipe more violent. The most obvious flow separation which is characterized by a largest vortex is observed when IGVF = 15%, and this results in a largest area of big bubbles in the diffuser.
(3) Small d 0 with a constant IGVF indicates small distance but large interaction between bubbles, so the bubble diameter increases obviously. The diffusion of bubbles is easier at a smaller d 0 , and bubbles tend to accumulate near the hub region with an increased d 0 .
The flow mechanism is not totally the same between a stirred tank and a pump, so some adjustments are needed when the bubble number density equation is applied in this work. For example, the fitted constants about the rates of break-up and coalescence may be improved to reflect more precisely the flow mechanism as well as the rotor-stator interaction. Additionally, based on the study of the single-stage multiphase pump, simulations and experiments on the multi-stage pump should be conducted in the next step.