Numerical Investigation on Bubble Distribution of a Multistage Centrifugal Pump Based on a Population Balance Model

: This work aimed to study the bubble distribution in a multiphase pump. A Euler-Euler inhomogeneous two-phase ﬂow model coupled with a discrete particle population balance model (PBM) was used to simulate the whole ﬂow channel of a three-stage gas-liquid two-phase centrifugal pump. Comparison of the computational ﬂuid dynamic (CFD) simulation results with experimental data shows that the model can accurately predict the performance of the pump under various operating conditions. In addition, the liquid phase velocity distribution, gas-phase distribution, and pressure distribution of the second stage impeller at a 0.5 span of blade height under three typical working conditions were compared. Results show that the region with high local gas volume fraction (LGVF) mainly appears on the suction surface (SS) of the blade. With the increase in inlet gas volume fraction (IGVF), vortices and low velocity recirculation regions are generated at the impeller outlet and SS of the blade, the area with high LGVF increases, and gas–liquid separation occurs at the SS of the blade. The liquid phase ﬂows out of the impeller at high velocity along the pressure surface of the blade, and the limited pressurization of ﬂuid mainly happens at the impeller outlet. The average bubble size at the impeller outlet is the smallest while that at the impeller inlet is the largest. Under low IGVF conditions, bubbles tend to break into smaller ones, and the broken bubbles mainly concentrate at the blade pressure surface (PS) and the impeller outlet. Bubbles tend to coalesce into larger ones under high IGVF conditions. With the increase in IGVF, the bubble aggregation zone di ﬀ uses from the blade SS to the PS. the higher the bubble size fraction of the coalesced bubbles, which mainly accumulate near the SS of the blade leading edge and di ﬀ use toward the blade PS and the middle-downstream of the ﬂow channel with the increase in IGVF. This study provides a reference for the numerical simulation of gas-liquid two-phase ﬂow in rotating machinery, and it helps for the design and optimization of gas-liquid two-phase ﬂow centrifugal pumps through the research of bubble distribution.


Introduction
Multistage centrifugal pumps have been widely used in the field of oil gas mixed transportation due to their advantages of high head, high efficiency, and simple structure [1,2]. The multi-stage centrifugal pump needs to operate under different inlet gas volume fraction (IGVF) conditions, as the volume fraction of oil and gas always changes during crude oil production and transportation. At low IGVF conditions, the pump has a higher head and efficiency, but under high IGVF conditions, its pressurization capacity drops sharply or it even stops working, which seriously endangers the continuity and safety of industrial production [3]. Therefore, it is of great significance to carry out an investigation on gas-liquid two-phase flow in a multistage centrifugal pump and improve its performance at IGVF conditions. Many scholars have studied the flow characteristics in gas-liquid two-phase flow pumps. Caridad [4] pointed out that the accumulation of gas-phase on the blade surface caused the performance

Experimental Test Rig
The gas-liquid two-phase experimental system of a three-stage centrifugal pump is shown in Figure 1. It is a branch system of a large-scale deep-sea gathering and transmission experimental system of State Key Laboratory of Multiphase Flow in Power Engineering, Xi'an Jiaotong University.

Experimental Test Rig
The gas-liquid two-phase experimental system of a three-stage centrifugal pump is shown in Figure 1. It is a branch system of a large-scale deep-sea gathering and transmission experimental system of State Key Laboratory of Multiphase Flow in Power Engineering, Xi'an Jiaotong University. This system is composed of a gas-phase line, a liquid-phase line, and a gas-liquid mixture line. In the gas-phase line, air passes through the compressor and enters the gas buffer tank to stabilize the pressure. After passing through the gas mass flowmeter, it finally enters the gas-liquid mixing section. In the liquid-phase line, the water from the water tank is pressurized by the plunger pump, passes through the liquid mass flowmeter, and finally enters the gas-liquid mixing section to mix with the air. In the gas-liquid mixture line, the gas-liquid phase mixed by the mixer enters the experimental pump and finally returns to the gas-liquid separator. The test rig of the experimental section is shown in Figure 2. The RHM30 and RHM06 high-precision mass flow meters were used to measure the flow rate of air and water. The gas and liquid phase measurement accuracy is 0.5% and 0.15%, respectively. The inlet pressure, outlet pressure, and pressure differentials between pump inlet and outlet of the experimental pump were measured by Rosemount 3051 series pressure and differential pressure sensors with a measurement accuracy of 0.075%. The HY-5005 torque tachometer was used to measure the torque and rotational speed of the pump, and its accuracy is 0.1%. The experimental test maintains the pump rotational speed, inlet pressure, and liquid phase mass flow rate under different IGVFs at an identical value to investigate the effect on the pump performance by adjusting the air volume fraction λg at the pump inlet. The experimental parameters are as follows: pump rotational speed is 3500 rpm, high pressure inlet, water mass flow rate is 7.2 kg/s, and the IGVF λg of the pump ranges from 0% to 25.7%. The expression of the IGVF λg is shown in Equation (1).
where, Qg and Ql represent the gas and liquid volume flow rate, respectively, in m 3 /h. λl represents the liquid volume fraction, in %.
This system is composed of a gas-phase line, a liquid-phase line, and a gas-liquid mixture line. In the gas-phase line, air passes through the compressor and enters the gas buffer tank to stabilize the pressure. After passing through the gas mass flowmeter, it finally enters the gas-liquid mixing section. In the liquid-phase line, the water from the water tank is pressurized by the plunger pump, passes through the liquid mass flowmeter, and finally enters the gas-liquid mixing section to mix with the air. In the gas-liquid mixture line, the gas-liquid phase mixed by the mixer enters the experimental pump and finally returns to the gas-liquid separator. The test rig of the experimental section is shown in Figure 2. The RHM30 and RHM06 high-precision mass flow meters were used to measure the flow rate of air and water. The gas and liquid phase measurement accuracy is 0.5% and 0.15%, respectively. The inlet pressure, outlet pressure, and pressure differentials between pump inlet and outlet of the experimental pump were measured by Rosemount 3051 series pressure and differential pressure sensors with a measurement accuracy of 0.075%. The HY-5005 torque tachometer was used to measure the torque and rotational speed of the pump, and its accuracy is 0.1%. The experimental test maintains the pump rotational speed, inlet pressure, and liquid phase mass flow rate under different IGVFs at an identical value to investigate the effect on the pump performance by adjusting the air volume fraction λ g at the pump inlet. The experimental parameters are as follows: pump rotational speed is 3500 rpm, high pressure inlet, water mass flow rate is 7.2 kg/s, and the IGVF λ g of the pump ranges from 0% to 25.7%. The expression of the IGVF λ g is shown in Equation (1).
where, Q g and Q l represent the gas and liquid volume flow rate, respectively, in m 3 /h. λ l represents the liquid volume fraction, in %.

Computing Domain and Mesh Generation
The gas-liquid two-phase flow pump used in the experiment is a three-stage centrifugal pump with a specific speed of 26.6. The pump design volume flow rate Qd is equal to 26 m 3 /h, design head Hd is equal to 25 m, and design rotational speed nd is equal to 3500 r/min. The number of impeller blades is 7, and the impeller outlet diameter is 127 mm. Figure 3 shows the computational domain of the pump model, including the inlet pipe, the first-stage pump, the second-stage pump, the finalstage pump, and the outlet pipe. Each stage of the pump has the same structure and is composed of an impeller, a cavity, and a diffuser.

Computing Domain and Mesh Generation
The gas-liquid two-phase flow pump used in the experiment is a three-stage centrifugal pump with a specific speed of 26.6. The pump design volume flow rate Q d is equal to 26 m 3 /h, design head H d is equal to 25 m, and design rotational speed n d is equal to 3500 r/min. The number of impeller blades is 7, and the impeller outlet diameter is 127 mm. Figure 3 shows the computational domain of the pump model, including the inlet pipe, the first-stage pump, the second-stage pump, the final-stage pump, and the outlet pipe. Each stage of the pump has the same structure and is composed of an impeller, a cavity, and a diffuser.

Computing Domain and Mesh Generation
The gas-liquid two-phase flow pump used in the experiment is a three-stage centrifugal pump with a specific speed of 26.6. The pump design volume flow rate Qd is equal to 26 m 3 /h, design head Hd is equal to 25 m, and design rotational speed nd is equal to 3500 r/min. The number of impeller blades is 7, and the impeller outlet diameter is 127 mm. Figure 3 shows the computational domain of the pump model, including the inlet pipe, the first-stage pump, the second-stage pump, the finalstage pump, and the outlet pipe. Each stage of the pump has the same structure and is composed of an impeller, a cavity, and a diffuser.

Computing Domain and Mesh Generation
The gas-liquid two-phase flow pump used in the experiment is a three-stage centrifugal pump with a specific speed of 26.6. The pump design volume flow rate Qd is equal to 26 m 3 /h, design head Hd is equal to 25 m, and design rotational speed nd is equal to 3500 r/min. The number of impeller blades is 7, and the impeller outlet diameter is 127 mm. Figure 3 shows the computational domain of the pump model, including the inlet pipe, the first-stage pump, the second-stage pump, the finalstage pump, and the outlet pipe. Each stage of the pump has the same structure and is composed of an impeller, a cavity, and a diffuser.    Five grid numbers are created to validate the grid independence of the pump model. The details are listed in Table 1. Where n represents the grid scheme (n = 1,2,3,4,5). H n /H 1 and η n /η 1 represent the ratio of the head and efficiency of the pump model with grid n to the grid 1. As shown in Table 1, when the grid total number of the pump model reaches 7.3223 million, the values of H n /H 1 and η n /η 1 slightly change with increasing grid number. Therefore, after considering the calculation resources and the accuracy of the numerical calculation results comprehensively, grid 3 was finally selected for numerical simulation research.

Two-Fluid Model of Gas-Liquid Two-Phase Flow
In the two-phase flow, each phase diffuses randomly in space and time, and there is dynamic interaction at the same time. Scholars have put forward a variety of mathematical models to solve this complex problem. The homogeneous model does not consider the difference of the two phases at all. The inhomogeneous model considers the slip velocity and the interphase force between the gas and liquid phases, which is closer to the actual two-phase flow. Therefore, the Euler-Euler two-fluid model was used to simulate the gas-liquid two-phase flow in the centrifugal pump. This model assumes that the liquid phase is a continuous fluid medium and the gas-phase is a dispersed fluid medium, and all the media fills the entire flow field. The mass and momentum conservation equations for the liquid and gas-phases are listed separately, and the two sets of equations are coupled together through the interaction between the phase interfaces. Due to the strong compression characteristics of the gas-phase, the concentration gradient of the gas in the flow field causes mass diffusion. Therefore, the Navier-Stokes equation [19], considering the medium mass diffusion, was adopted. The mass conservation equation is shown in Equation (2), and the momentum conservation equation is shown in Equation (3).
where i represents liquid phase l and gas-phase g. In this paper, the gas-liquid two-phase working medium is air and tap water, respectively. ρ i , α i are the density and volume fraction of the ith phase, respectively. u i represents the velocity component in the three directions of XYZ. The source terms S i and S Mi represent the transmission of the gas-phase mass and momentum during bubble breakup and coalescence. M i represents the sum of the forces acting between the gas and liquid phase, and τ i represents the strain force tensor of the ith phase.

Discrete Particle Population Balance Model (PBM)
Discrete PBM was used to simulate the distribution of the number and size of bubbles in the flow field. This model considers the breakup and aggregation of bubbles due to the forces between the gas and liquid phases. The bubble breakup model used was the Luo and Svendsen model [22], and the bubble coalescence model used was the Prince and Blanch model [23]. The mass transfer S i of the gas-phase can be obtained from Equation (4): where B B , D B , B C , and D C , respectively, represent the bubble birth rates due to the breakup of larger particles, the bubble death rate due to the breakup of larger particles, the bubble birth rate due to coalescence of smaller particles, and the bubble death rate due to the coalescence of smaller particles. In this paper, 10 groups of bubbles with different sizes were selected. All bubbles in the same group have the same diameter. The distribution of the size and number of these bubbles obeys the balance of discrete particle population. It is known that the bubble diameter in the pump is mainly distributed between 0.1 and 10mm according to the experimental research of scholars [6,8]. In order to improve the calculation efficiency, this bubble size range was used in this work. The minimum bubble diameter d0 is given as 0.1 mm, and the maximum bubble diameter d9 is given as 10 mm. The bubble size [24] calculation of each group is shown in Equations (5) and (6).
The bubble sizes of the 10 groups are shown in Table 2. At the pump inlet, the bubble size fraction of the groups with medium-diameter bubbles [21] d 4 and d 5 are set to 0.5, respectively, and the other bubble size groups are set to 0. In the numerical calculation, the bubble size fraction in the flow field reaches a balanced state through the bubble breakup and coalescence. The bubble size fraction refers to the ratio of the volume flow rate of the group to the total volume flow rate of the gas-phase. The sum of the bubble size fraction of 10 bubble groups is 1.00 under each working condition.

Basic Assumptions and Parameter Settings
The basic assumptions for the numerical simulation of the gas-liquid two-phase pump model are as follows: (1) The gas-liquid two-phase flow pattern at the pump inlet is bubble flow. (2) The pump inlet bubble is spherical with diameter d.  Table 3.

Experimental Validation
The pressure increment coefficient C p is defined as the ratio of the pump pressure increment value ∆P and the experimental value ∆P 0 under pure water conditions (λ g = 0%), as shown in Equation (7).
The homogeneous flow model and inhomogeneous flow model coupled with the particle size model (single particle size model and PBM model) are used for numerical simulations in this work. The comparison between the numerical simulation results and the experimental results is shown in Figure 5. The error analysis is shown in Figure 6. The calculation results using the inhomogeneous flow coupled with PBM agree the best with the experimental values, followed by the inhomogeneous flow model, and the homogeneous flow model that has the largest errors. The numerical simulation errors under different conditions are shown in Table 4.
The homogeneous flow model does not consider the slip velocity [26] between the gas and liquid phases. As the gas volume fraction increases, the slip velocity between the gas and liquid phases also increases. Therefore, the homogeneous model is applicable to the low gas volume fraction conditions of bubble flow or mist flow [27,28]. The numerical simulation results of the inhomogeneous flow coupled with the bubble size of 0.1 mm are larger than the experimental values, whereas the numerical simulation results of the inhomogeneous flow coupled with the bubble size of 0.2mm are lower than the experimental values at C2 and C3, but higher than the experimental values at C4, C5, C1-C6 refers to six operating conditions with different IGVFs. From C1 to C6, the IGVF λ g is 0%, 3.5% (±0.0175%), 10% (±0.05%), 13.5% (±0.0675%), 18.53% (±0.09265%), and 25.7% (±0.1285%), respectively. It can be seen from Figure 5 that the uncertainty, the experimental values, and the numerical simulation results using the PBM model are within the range of the measurement uncertainty. The experimental uncertainties correspond to stationary measurements. The pressure increment curve of the pump decreases with the increase of the IGVF under the condition that the liquid phase flow Energies 2020, 13, 908 8 of 15 rate, inlet pressure, and rotation speed were the same. The curve descending process was mainly divided into three stages. The first stage is from C1 to C3, the second stage is from C3 to C4, and the third stage is from C4 to C6. The three curves are all approximately linear functions. The slope of the second curve is the largest, which indicates that the pressure increment of the pump suddenly decreases, and the sudden drop value is 22.3% of the pressure increment value in the C3 operating condition. This phenomenon of pressure increment drop is called surging [25].
The error analysis is shown in Figure 6. The calculation results using the inhomogeneous flow coupled with PBM agree the best with the experimental values, followed by the inhomogeneous flow model, and the homogeneous flow model that has the largest errors. The numerical simulation errors under different conditions are shown in Table 4. The error analysis is shown in Figure 6. The calculation results using the inhomogeneous flow coupled with PBM agree the best with the experimental values, followed by the inhomogeneous flow model, and the homogeneous flow model that has the largest errors. The numerical simulation errors under different conditions are shown in Table 4.
The homogeneous flow model does not consider the slip velocity [26] between the gas and liquid phases. As the gas volume fraction increases, the slip velocity between the gas and liquid phases also increases. Therefore, the homogeneous model is applicable to the low gas volume fraction conditions of bubble flow or mist flow [27,28]. The numerical simulation results of the inhomogeneous flow coupled with the bubble size of 0.1 mm are larger than the experimental values, whereas the numerical simulation results of the inhomogeneous flow coupled with the bubble size of 0.2mm are lower than the experimental values at C2 and C3, but higher than the experimental values at C4, C5, and C6.   The homogeneous flow model does not consider the slip velocity [26] between the gas and liquid phases. As the gas volume fraction increases, the slip velocity between the gas and liquid phases also increases. Therefore, the homogeneous model is applicable to the low gas volume fraction conditions of bubble flow or mist flow [27,28]. The numerical simulation results of the inhomogeneous flow coupled with the bubble size of 0.1 mm are larger than the experimental values, whereas the numerical simulation results of the inhomogeneous flow coupled with the bubble size of 0.2mm are lower than the experimental values at C2 and C3, but higher than the experimental values at C4, C5, and C6.
The coalesced bubbles in the flow field reduce the effective flow area of the liquid phase in the flow channels, and the large gas pocket blocks the flow channel, resulting in decreased performance of Energies 2020, 13, 908 9 of 15 the pump. The probability of bubble coalescence into large bubbles increases with increasing bubble diameter and gas volume fraction. Therefore, when using a single-bubble size model, the initial bubble size d needs to be reset as the operating conditions change [29][30][31].
In the gas-liquid two-phase flow field, the bubble size distribution in the gas-liquid two-phase flow changes with the occurrence of multiphase system reactions and transfer phenomena, and the coalescence and fragmentation between the bubbles cause the particle size to change. The PBM can make up for the shortcomings of the single-bubble size model well. There are three main characteristics of the PBM: (1) Multiple groups of bubble diameters are set, as shown in Table 2, and the mass transfer, break up, and coalescence between bubbles are considered. (2) Various bubble sizes coexist and the numbers are balanced with each other. (3) Under small gas volume fraction conditions, large bubbles are broken into small ones, which then coalesce into large particles under large gas volume fraction conditions until the particle group reaches a balanced state.

Internal Flow Field Analysis
The first and last stages of the three-stage pump were affected by the inflow of the inlet and outlet. Therefore, the internal flow field was analyzed on the second-stage impeller to avoid the influence of other inflows on the internal flow field of the pump. The IGVF is the main factor that affects the performance of the pump. Therefore, three typical working conditions of the pump model (at the identical liquid flow rates) were numerically simulated. The IGVF λ g for the three operating conditions is 3.5%, 13.5%, and 25.7%, as shown in Figure 7. Figure 7a-c show the local gas-phase volume fraction (LGVF) α g distribution, liquid phase streamline distribution, and static pressure increment distribution at a 0.5 span of the impeller in the second-stage pump, respectively.  It shows that with the increase of IGVF, the flow field distribution in the impeller changes drastically, especially from the IGVF of 3.5% to 13.5%. The flow loss is the largest at this stage, and the external characteristics of the pump drop suddenly (as shown in Figure 5).
At high IGVF conditions, the area with high LGVF near the suction surface also increases, indicating that the gas-phase has accumulated and stuck on the suction surface. At the same time, it can be seen that the accumulation of gas-phase makes the liquid flow space be severely squeezed, which causes the relative velocity of the liquid phase near the pressure surface (PS) of the blade to increase significantly, and leading to large flow losses in the pump.
As shown in Figure 7a, the gas-phase in the impeller flow channel is mainly distributed near the suction surface (SS) of the blade leading edge under low IGVF conditions, and it diffuses to the PS of the blade and downstream of the flow channel as the IGVF increases. There are three main reasons for the gas-phase separation and why it is stuck near the SS of the impeller blade: (1) The gas and liquid exhibit flow velocity separation at the blade inlet, as shown in Figure 7b. The flow direction of the liquid phase at the inlet of the blade is mainly from the SS to the PS of the blade. Due to the low density of gas-phase, it obtains less energy from the impeller blades, resulting in its velocity being lower than that of the liquid phase, which causes the gas-phase to stay at the leading edge of the blade SS. (2) The high-pressure and high velocity fluid near the outlet of the blade PS flows towards the SS of the same blade generating a jet region and a vortex region in an anticlockwise direction. The jet flows towards the outlet of the impeller while taking away the bubbles near the outlet of the impeller. This vortex is located inside of the jet region, as shown by the black arrow line in Figure 7b. The vortex region blocks the fluid near the suction surface, forms a backflow and vortex, and prevents the bubbles from flowing out of the impeller. (3) Comparing Figure 7a and 7b, it can be seen that under the condition of high IGVF, the region with high LGVF near the suction surface corresponds to the liquid vortex and low velocity recirculation zone, as shown by the red arrow line in Figure 7b. In this region, the relative velocity direction of the liquid phase is directed from the impeller outlet to the inlet, and the drag force of the bubble is directed to the impeller inlet along the liquid streamline. Therefore, it is difficult for the bubbles to be carried away from the impeller channel by the liquid, causing the gas to aggregate at the liquid phase vortex and gas-liquid separation.
As shown in Figure 7c, the pressure increment value of static pressure increases uniformly in the radial direction from the impeller inlet to the impeller outlet under low IGVF conditions. Under high IGVF conditions, a reverse pressure gradient zone is generated near the blade PS of the downstream impeller flow channel. The pressure increment value is almost unchanged within 0.7 times the relative radius of the impeller and the value increases slightly near the outlet of the impeller. The main reason for the uneven pressure distribution is that there is a large amount of gas accumulation in the flow channel under high IGVF conditions, and the energy applied by the impeller to the fluid is dissipated by the backflow. Therefore, the pressurized area is mainly concentrated downstream of the impeller flow channel, and the value of pressure increment decreases rapidly with the increase of the IGVF.

Distribution of Bubbles with Different Particle Sizes in the Pump
The numerical simulation results of the discrete particle population balance model were analyzed to understand the distribution of bubbles with different particle sizes in the pump under different operating conditions. The bubble size fraction of 10 bubble groups (Table 2) at a 0.5 span of the pump impeller and guide vane is shown in Figure 8.
The curve of bubble size fraction is normally distributed with the bubble size. When λ g is 3.5%, 10%, 13.5%, 18.53%, and 25.7%, the maximum bubble size fraction corresponds to d 6 , d 5 , d 5 , d 4 , and d 4 , respectively. The bubble diameter corresponding to the maximum bubble size fraction is larger as the IGVF increases. The bubble size fraction of the bubble group with a diameter of d 9 (0.1 mm) is the largest at the 3.5% operating condition, which is about four times that of the 25.7% operating conditions. The bubble size fraction of the bubble group with diameter d 0 (10 mm) under 3.5%, 10%, and 13.5% operating conditions approaches 0, and it reaches the maximum at the 25.7% operating condition.
Considering that the bubble size fraction of both the d 4 and d 5 bubble groups at the pump inlet is 0.5, the bubbles with the size of d 4 and d 5 in the flow field are assumed to be neither broken nor coalesced. The bubble size fraction with neither broken nor coalesced bubbles is shown by the curve between the two black dotted lines in Figure 8. The curve on the left side of the black dotted line indicates the bubble size fraction of broken bubbles with diameters smaller than d 5 , and the curve on the right side of the black dashed line represents the bubble size fraction of the coalesced bubbles with diameters larger than d 4 . Under different IGVF conditions, the value of bubble size fraction at the areas of broken bubbles and coalescence bubbles is shown in Table 5. It shows that the bubbles tend to break from large to small bubbles under small IGVFs while coalesce from small to large bubbles under large IGVFs.
10%, 13.5%, 18.53%, and 25.7%, the maximum bubble size fraction corresponds to d6, d5, d5, d4, and d4, respectively. The bubble diameter corresponding to the maximum bubble size fraction is larger as the IGVF increases. The bubble size fraction of the bubble group with a diameter of d9 (0.1 mm) is the largest at the 3.5% operating condition, which is about four times that of the 25.7% operating conditions. The bubble size fraction of the bubble group with diameter d0 (10 mm) under 3.5%, 10%, and 13.5% operating conditions approaches 0, and it reaches the maximum at the 25.7% operating condition.
Considering that the bubble size fraction of both the d4 and d5 bubble groups at the pump inlet is 0.5, the bubbles with the size of d4 and d5 in the flow field are assumed to be neither broken nor coalesced. The bubble size fraction with neither broken nor coalesced bubbles is shown by the curve between the two black dotted lines in Figure 8. The curve on the left side of the black dotted line indicates the bubble size fraction of broken bubbles with diameters smaller than d5, and the curve on the right side of the black dashed line represents the bubble size fraction of the coalesced bubbles with diameters larger than d4. Under different IGVF conditions, the value of bubble size fraction at the areas of broken bubbles and coalescence bubbles is shown in Table 5. It shows that the bubbles tend to break from large to small bubbles under small IGVFs while coalesce from small to large bubbles under large IGVFs.    The bubble size fractions of the second stage pump at the impeller inlet, impeller outlet, cavity outlet, and diffuser outlet under 3.5% and 25.7% operating conditions were analyzed, as shown in Figure 9. The bubble size fraction curve at the impeller inlet is on the right end, and that at the impeller outlet is on the left end. The rate of bubble breakup at the impeller outlet is greater than the rate of coalescence due to the shear stress, and the bubble size fraction curve moves in the direction of a small diameter. However, the bubble coalescence rate at the cavity outlet is greater than the breaking rate, which causes the bubble size fraction of the larger bubbles to increase, and the curve moves toward the direction of large bubbles. Similarly, the bubble size fraction curve at the diffuser outlet moves to the right of the horizontal axis. The bubble size fraction under the IGVF 25.7% operating condition moves to the direction of the large bubble size compared to the 3.5% operating condition. The reason for this phenomenon is that as the IGVF increases, the distance between the bubbles decreases, and the probability of bubble collision increases. It is easy to aggregate from small bubbles into large bubbles. Therefore, the bubble size fraction of the large bubble group increases, and the curve moves to the right of the horizontal axis. outlet moves to the right of the horizontal axis. The bubble size fraction under the IGVF 25.7% operating condition moves to the direction of the large bubble size compared to the 3.5% operating condition. The reason for this phenomenon is that as the IGVF increases, the distance between the bubbles decreases, and the probability of bubble collision increases. It is easy to aggregate from small bubbles into large bubbles. Therefore, the bubble size fraction of the large bubble group increases, and the curve moves to the right of the horizontal axis. In this paper, the isosurface where the sum of the bubble size fraction of the coalesced bubbles is greater than or equal to 0.9 was used to represent the bubble accumulation zone, and the isosurface where the sum of the bubble size fraction of broken bubbles is greater than or equal to 0.9 was used to represent the broken bubble accumulation zone. Figure 10a,b shows the distribution of the accumulation area and the accumulation area of the coalesced bubbles in the second-stage impeller. In this paper, the isosurface where the sum of the bubble size fraction of the coalesced bubbles is greater than or equal to 0.9 was used to represent the bubble accumulation zone, and the isosurface where the sum of the bubble size fraction of broken bubbles is greater than or equal to 0.9 was used to represent the broken bubble accumulation zone. Figure 10a,b shows the distribution of the accumulation area and the accumulation area of the coalesced bubbles in the second-stage impeller. outlet moves to the right of the horizontal axis. The bubble size fraction under the IGVF 25.7% operating condition moves to the direction of the large bubble size compared to the 3.5% operating condition. The reason for this phenomenon is that as the IGVF increases, the distance between the bubbles decreases, and the probability of bubble collision increases. It is easy to aggregate from small bubbles into large bubbles. Therefore, the bubble size fraction of the large bubble group increases, and the curve moves to the right of the horizontal axis. In this paper, the isosurface where the sum of the bubble size fraction of the coalesced bubbles is greater than or equal to 0.9 was used to represent the bubble accumulation zone, and the isosurface where the sum of the bubble size fraction of broken bubbles is greater than or equal to 0.9 was used to represent the broken bubble accumulation zone. Figure 10a,b shows the distribution of the accumulation area and the accumulation area of the coalesced bubbles in the second-stage impeller.  Figure 10a shows almost no coalesced bubbles in the impeller at lower IGVFs (λ g = 3.5%). The coalesced bubbles begin to gather near the SS of the blades leading edge at the 13.5% operating condition. With the further increase in IGVF, the accumulation area of the coalesced bubbles diffuses to the PS of the blade and the middle-downstream of the flow channel at the 25.7% operating condition. The reason is that the gas-liquid separation causes a large number of coalesced bubbles at the leading edge of the SS forming gas pocket. Therefore, the gas-phase volume fraction in this region is higher, shown as the black circle part in Figure 7a.
As shown in Figure 10b, the accumulation area of broken bubbles is distributed throughout the impeller at the 3.5% operating condition. As IVGF increases, the broken bubbles begin to gather near the blade PS and the impeller outlet, and almost no coalesced bubbles appear in this area. The reason for the accumulation of broken bubbles near the blade PS and the impeller outlet is that the secondary flow and jet phenomenon cause the turbulence intensity of the flow field to increase, and the shear stress on the bubble increases. As a result, the bubbles easily break into small bubbles and do not easily coalesce into larger ones. The location of the broken bubble distribution area corresponds to the secondary flow and jet area. The small bubbles tend to flow out of the impeller due to the drag force of water, and the volume fraction of the gas-phase in this area is lower, shown as the red circle part in Figure 7a.

Conclusions
Experiments and numerical simulations of a gas-liquid two-phase three-stage centrifugal pump were conducted. The following conclusions were obtained: (1) The Euler-Euler inhomogeneous flow coupled discrete particle population balance model was applied to a gas-liquid two-phase three-stage centrifugal pump. The maximal discrepancy between the numerical simulation results and the experimental values is 5%, which verified the numerical model.
(2) The area with high LGVF diffuses from the suction of the blade leading edge toward the blade PS and downstream of the flow channel as the IGVF increases. The jet flow at the impeller outlet causes a liquid-phase vortex rotating counterclockwise. This vortex blocks the fluid on the suction surface, which results in the local vortices and backflow. The liquid phase vortex rotating clockwise and the backflow near the SS trap the gas in the flow channel. Under high IGVF conditions, the pressurized area is mainly located at the downstream region of the impeller, and the pressurized value is very small.
(3) The average bubble size at the impeller outlet is the smallest, followed by the cavity outlet and then the diffuser outlet, and the average bubble size at the impeller inlet is the largest.
(4) Bubbles tend to break into smaller bubbles under small IGVF conditions, while bubbles tend to aggregate into larger bubbles under large IGVF conditions. The lower the IGVF, the higher the bubble size fraction of the broken bubbles, which mainly accumulate near the blade pressure surface and the impeller outlet. And the higher the IGVF, the higher the bubble size fraction of the coalesced bubbles, which mainly accumulate near the SS of the blade leading edge and diffuse toward the blade PS and the middle-downstream of the flow channel with the increase in IGVF.
This study provides a reference for the numerical simulation of gas-liquid two-phase flow in rotating machinery, and it helps for the design and optimization of gas-liquid two-phase flow centrifugal pumps through the research of bubble distribution.
Author Contributions: Numerical analysis and design of experiments were conducted by S.Y., S.S., X.L., S.C., C.L. and J.F. Analysis of the results of numerical simulation was conducted by S.Y. and S.S. S.Y. wrote the manuscript with revision by S.S. and X.L. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.