State of the Art on Two-Phase Non-Miscible Liquid/Gas Flow Transport Analysis in Radial Centrifugal Pumps Part C: CFD Approaches with Emphasis on Improved Models

: Predicting pump performance and ensuring operational reliability under two-phase conditions is a major goal of three-dimensional (3D) computational ﬂuid dynamics (CFD) analysis of liquid/gas radial centrifugal pump ﬂows. Hence, 3D CFD methods are increasingly applied to such ﬂows in academia and industry. The CFD analysis of liquid/gas pump ﬂows demands careful selection of sub-models from several ﬁelds in CFD, such as two-phase and turbulence modeling, as well as high-quality meshing of complex geometries. This paper presents an overview of current CFD simulation strategies, and recent progress in two-phase modeling is outlined. Particular focus is given to different approaches for dispersed bubbly ﬂow and coherent gas accumulations. For dispersed bubbly ﬂow regions, Euler–Euler Two-Fluid models are discussed, including population balance and bubble interaction models. For coherent gas pocket ﬂow, essentially interface-capturing Volume-of-Fluid methods are applied. A hybrid model is suggested, i.e., a combination of an Euler– Euler Two-Fluid model with interface-capturing properties, predicting bubbly ﬂow regimes as well as regimes with coherent gas pockets. The importance of considering scale-resolving turbulence models for highly-unsteady two-phase ﬂow regions is emphasized.


Introduction
In various industrial and engineering applications, a mixture of liquid and noncondensable gas, such as air, is encountered. Centrifugal pumps, which have been designed for single-phase liquid transport, are frequently required to convey two-phase mixtures. Twophase pump applications can be found in nuclear [1] and geothermal [2] power stations, in the pulp and paper industry [3] or in the petroleum industry for artificial lifting [4]. Even a small load of gas may considerably disturb conveyance, and may even cause a break-down of the pump head [5]. As discussed in detail in Part A of this review paper series, a drop in the pump head is associated with a separation of air and water in the blade channels [3,[6][7][8].
In pumps with high specific speeds, in terms of axial and mixed flow pumps, the action of the Coriolis force is directed from the shroud to the hub [9]. Thereby, secondary flow is promoted in the blade channels and this, to some extent, counteracts phase separation. In circumferential pumps, the Coriolis force acts in the direction of the blade pressure side. As a consequence, pressure gradients in the cross-flow direction develop. Since the densities of liquid and gas considerably deviate, pressure gradients lead to a slip between phases, and the gas tends to accumulate on the blade suction side, especially in part-load conditions, where the flow incidence at the blade trailing edge is high [10].
Although there is no difference in principle in the interaction of body forces in the liquid and the pressure gradients, which finally leads to phase segregation, the effect on segregation is considerably more pronounced for circumferential pumps, due to the pockets are present in the transition zone between bubbly and pocket flow regimes. Even inherently steady adherent void regions in the pocket flow regime show unsteady wakes, which will be discussed in Section 4. This observation should be taken into account by reasonable turbulence models. For highly unsteady flow, the validity of Reynoldsaveraging might be challenged, and turbulence-scale resolving models might be more accurate. Therefore, a brief review of turbulence-scale resolving simulations of pump flow is provided in Section 2.5.   [27,39]. Some operational points are marked with an "X", and are slightly pulled apart for clarity due to their proximity. They are picked up in the flow simulations further below.

Euler-Euler Two-Fluid (EE2F) Approach
Most available studies on radial centrifugal pump flow have been performed using the EE2F model family. An Eulerian approach for the dispersed phase is preferred to the Lagrangian approach, since more moderate grid dependence and better statistical convergence of the Eulerian approach may be assumed [45]. Early simulations of two-phase flow in centrifugal pumps by Pak and Lee [46] and Minemura and Uchiyama [47] pointed out that bubbles move to the blade tip and accumulate on the blade's suction side and that the head drop is associated with phase separation and large air accumulation zones. Even with a simplified monodisperse approach and single-channel simulations, Müller et al. [48] obtained good agreement with the measured head drop at low air loading, while, with a rise in the inlet gas volume fraction, even a qualitative prediction could not be obtained, confirming the conclusion of several previous studies, e.g., [49][50][51]. A significant influence of the computational grid quality on the location and size of bubble accumulation was found by Müller et al. [52]. Si et al. [40][41][42] and Wang et al. [53] presented validated numerical results of a 3D-bladed centrifugal pump in two-phase air-water flow and compared different impeller designs. Zhou et al. [54] investigated transient, i.e., time-resolved, characteristics of the radial pump force, and Zhong et al. [55] analyzed the transient pump characteristics that start when air is injected in the suction pipe. The EE2F model was also adopted in our previous studies, comparing a closed and a semi-open impeller [56][57][58]. Several authors have utilized polydisperse bubble size distribution, together with population balance modeling (PBM) within the EE2F framework, enabling temporal and spatial variation of the bubble size. He et al. [59] highlighted the importance of considering the change of bubble diameter for volute-type centrifugal pumps and compared these results to the results of a monodisperse EE2F simulation. An example result obtained by He et al. [59], in terms of the distribution of the air volume fraction α a , is given in Figure 2. Stronger air aggregation was observed with the EE2F-PBM model than with the monodisperse approach, supposedly because the chosen monodisperse bubble diameter was too small. The benefit of using a PBM in combination with the EE2F approach is emphasized by having a better match to experimental data.
Air volume fraction α a a) b) Figure 2. Development of air accumulations, calculated by the coupled EE2F-PBM model (left) and the monodisperse EE2F model (right) for ε = 2.21% (a) and ε = 4.86% (b) at n = 1500 1/min and Q = 7.7 m 3 /h. The figure is adopted from He et al. [59].
Yan et al. [60] confirmed this conclusion. In Figure 3, the pump's pressure coefficient versus inlet gas volume fraction is given for the simulation results obtained by a coupled EE2F-PBM and by a monodosperse EE2F. For the latter, the bubble diameter d was varied. By an adjustment of d, the experimental data could be approached. Thus, d can be considered a kind of calibration parameter in the monodisperse EE2F variant, which is dispensable for the polydisperse EE2F-PBM. With this coupled EE2F-PBM model, Yan et al. [60] analyzed the bubble size distribution in a multistage centrifugal pump.
It is important to note that the term homogeneous (abbreviated "Hom." in Figure 3) used by the authors means that a common velocity field for water and air is assumed. This should not be confused with the homogeneous mixture approach, which means a treatment of bubbles on the subgrid scale and is inherent to any EE2F variant. The mismatch to data by using a common velocity field (result termed "Hom." in Figure 3) underlines the necessity of distinct velocity fields for water and air, to enable phase segregation.
In the simulations by Stel et al. [61,62], air accumulation zones within the impeller observed in the experiments of the same research group from Ref. [37] could be reproduced. Zhang et al. [63] observed large bubble sizes in regions of high air loading within the pump impeller. Chen et al. [64] analyzed the head drop of the first stage of an electrical submersible pump (ESP). Recently, Si et al. [65] reproduced the head drop observed in measurements of a five-stage ESP. They showed that the most significant portion of the head drop occurs in the first and the second stages under two-phase flow conditions, which is underlined by Figure 4, in terms of the stage head coefficient Ψ. For single-phase, i.e., pure water flow, the first stage does most of the work, according to Figure 4a. For two-phase conditions and ε > 2%, Ψ drops significantly for stages one and two, according to Figure 4b.  The EE2F model inherently presumes a dispersed bubbly flow [58], i.e., a homogeneous mixture of water and air within each computational cell. Thus, it can be assumed that in the EE2F approach, air accumulations are represented by clustering of densely dispersed bubbles. However, in the experimental investigations on a planar diffuser flow [66], a sharp interface between the water and the air accumulations was observed, which points out a coherent attached air structure, rather than a clustering of separated bubbles. It may be concluded that such large interface structures should be resolved by a VoF rather than an EE2F approach [67], which is outlined next.

Volume-of-Fluid (VoF) Approach
By utilizing a VoF approach, Parikh et al. [68] optimized an inducer design by employing a multi-objective optimization and assessed the influence of the inducer design on the pump flow. The same authors investigated the air distribution for different semi-open impellers [69], and compared these simulation results with the experimental results from the same research group [24,25,39]. For example, the air distribution is shown in Figure 5 for a variation of the inlet gas volume fraction. Parikh et al. [69] concluded that the instantaneous distribution of void regions observed in the experiments by a scattered-light technique could be reproduced well by the void fraction distribution in the simulation. However, steady air accumulations formed in the impeller channels for ε ≥ 3%, according to the experiments, which were not observed in the VoF simulation results. Pineda et al. [70] utilized the VoF method for flow simulations in ESPs and experienced up to 25% deviation from the results of the experimentally measured pump head. Zhu et al. [71] concluded that the VoF model is inferior to the EE2F model for the simulation of ESP flows. Mansour et al. [72] figured out the important physical mechanisms of attenuating pump performance and the beneficial effect of an inducer. However, it remains unclear if the spatial resolution limit allows an adequate resolution of void structures within the VoF simulations. According to De Santis et al. [73], VoF methods are only appropriate for capturing large interfacial structures, provided that a high spatial resolution is applied. Hundshagen et al. [74] showed that the grid resolution, which is feasible in centrifugal pumps, is by far too coarse to resolve all small scales encountered in bubbly, dispersed flow regimes. It is worth noting that, of course, like VoF methods, level set methods also provide an interface resolution. However, to the best knowledge of the authors, there are no centrifugal pump flow studies employing level set methods.

Hybrid Two-Phase (H2P) Approach
When using the EE2F and VOF approaches, two limiting morphological situations, encountered in radial centrifugal pumps, are addressed. Regarding the former, a homogeneous mixture is assumed, which means that a multitude of bubbles is dispersed within each computational cell and that the bubble diameter is much smaller than the cell size. This assumption is particularly valid in the bubbly flow region, while it is challenged for adherent air pockets. In particular, the validity of the EE2F model is challenged when the bubble size exceeds the grid scale, as discussed by, for example, Marschall [75]. Hundshagen et al. [56][57][58] showed that solver convergence issues may arise for the EE2F model, either in large coherent air regions or if the computational grid is very fine, i.e., if the bubble size locally exceeds the size of the computational cell. Regarding VoF schemes, the resolution of all small scales encountered in the bubbly flow region is not feasible with available computational resources. In Ref. [76], it was shown that, by adopting a hybrid ansatz to pump flow simulations, a transition of bubbly flow to coherent air accumulations could be accounted for even without the tremendous spatial grid resolution required for interface capturing by pure VoF schemes. In the hybrid models by Hänsch et al. [77,78], Meller et al. [79], Yin et al. [80], and Frederix et al. [81], air is treated either as a continuous or a disperse phase. Since CFD simulations of centrifugal pumps comprise complex grids and are computationally expensive, a preferable simple approach of a hybrid model has been adopted for this first application, which means that air is treated as a single phase, irrespective of its morphology. Thus, a local blending of the EE2F and VoF method, in terms of the interface compression, as proposed by, for example, Wardle and Weller [82], Shonibare and Wardle [83] and Mathur et al. [84], seems to be an appropriate choice. This hybrid approach comprises only one single momentum and volume fraction equation set for the air phase, which means that the solution of additional momentum and volume fraction equations for the air phase is omitted. We adopted the hybrid two-phase model of Wardle and Weller [82], with in-house extensions (called the H2P model in the following), including the interface compression concept, for the simulation of two-phase centrifugal pump flows. The capability of the H2P model was demonstrated on a research pump where optical experimental data had been measured by Mansour et al. [26,27]. Based on this data, together with head measurements, the transition from dispersed bubbly flow to air accumulations, termed the bubbly and pocket flow regimes, respectively, was categorized [26,27], as described above in Section 2.1 and summarized in Figure 1. Thus, a unique validation database is available for the H2P approach, which is described in more detail in Section 3 below.

Turbulence-Scale Resolving Approach and Scale-Adaptive Simulations (SASs)
We assume that inherently unsteady void structures may only inaccurately be resolved by statistical, i.e., unsteady Reynolds-averaged Navier Stokes (URANS) turbulence models. Even for single-phase flow, there are several studies that show the limitations of statistical models and the benefits of scale-resolving models in terms of large-eddy simulations (LES). Examples are the prediction of an inhomogeneous flow distribution at part load [85], head instabilities [86], part load instabilities and flow separation [87][88][89][90][91][92][93][94], tip vortices [95], pressure and velocity fluctuations due to rotor-stator interaction [96][97][98][99][100][101][102][103][104][105], as well as entropy production [106]. According to Refs. [107,108], 80% of the spectral energy should be resolved, so an adequately resolved LES demands an extensively high number of computational cells. Therefore, highly-efficient numerical schemes are employed, e.g., in terms of immersed boundary methods, as proposed by Posa et al. [91,92], Kim and Choi [109] and Kye et al. [101], or small Reynolds numbers are investigated, e.g., [103,106]. Posa et al. [97][98][99][100] and Kye et al. [101] utilized cylindrical computational grids with up to 500 Mio cells. By a finite-element method, with body-fitted overset grids, Pacot et al. [102] performed simulations at a reduced Reynolds number and estimated a number of 600 · 10 9 computational cells required for adequate resolution at realistic Reynolds numbers. These examples illustrate the high demands of the simulation method and the spatial resolution required for an adequate LES.
One main issue concerning a wall-resolved LES is the poor prediction of the near wall behavior of resolved quantities by a subgrid scale model, which can be overcome by hybrid URANS-LES turbulence models. Several hybrid methods in terms of detached eddy simulation (DES) [110,111], delayed DES (DDES) [112], Very LES [113], partially-averaged Navier-Stokes (PANS) [114,115] or partially-integrated transport model (PITM) [116] have been proposed. Hybrid turbulence models and, in particular, DES and DDES are increasingly utilized for pump flow simulations as their computational requirements are lower than those of LES [117,118]. In hybrid approaches and in near-wall regions with attached boundary layers, an URANS model is utilized, while, far away from the wall, an LES approach is employed. The local grid density is utilized as a transition criterion between LES and URANS [119]. One drawback of DES may be mesh-induced separation zones [120]. Beyond DDES, this issue has been addressed, e.g., by shielding functions [121,122]. A differ-ent approach for scale resolution is based on the exact length-scale equation by Rotta [123]. This scale-adaptive class of turbulence models can be understood as a class of enhanced statistical models [124], which may resolve the turbulent structures down to the grid limit. This is achieved by a local reduction of eddy-viscosity based on the von Kárman and integral length scale. According to Jakirlić et al. [125], the local grid resolution may be considered a further model parameter in the hybrid approaches listed above. Related to this assessment, it is interesting to note that, for the SAS, the filter which controls grid resolution is not immediately governed by the local grid resolution but rather by the local von Kármán length scale. A rather broadly used variant of this model is the k-ω-SST-SAS model [126], which can be understood as an extension of the URANS k-ω-SST model [120,127] towards scale resolution. A fall-back to an URANS solution in stationary flow regions or regions of low spatial and temporal resolution is an especially convenient way to avoid uncertainties in spectral energy resolution [128]. Several studies have shown the potential of scale-adaptive simulations (SAS) in highly transient pump flows compared to conventional statistical turbulence models, in terms of velocity distributions, pressure fluctuations, and integral characteristics [89,93,[129][130][131][132]. The present authors have also demonstrated the benefit of scale-adaptive simulations in single-phase flow in a centrifugal pump [133] and a pump mixer [134]. It was concluded that the SAS yields a sound prediction of the local flow field at only moderately enhanced computational effort compared to URANS. This observation was confirmed for multiphase flow in a centrifugal pump [74,76] and an impeller channel [58].

Summary of State-of-the-Art Method Algorithms
As summarized in Section 2, the EE2F and VoF schemes are utilized for the simulation of multiphase flow transport in circumferential pumps. The present authors have also used EE2F [48,52,[56][57][58]135] and VoF [74] approaches in their preceding studies, and these algorithms are briefly summarized here, before we proceed to recent enhancements, in terms of the H2P approach, in Section 3.2. Regarding the EE2F model, an Eulerian approach for the disperse phase is used. A dispersed bubbly air phase is assumed to be mixed within a continuous water phase, which means that each computational cell contains a homogeneous mixture of dispersed bubbles. Separate mass and momentum balance equations for water and air are solved, yielding a separate velocity field for each phase, according to Wardle and Weller [82]: Using separate velocity fields is often termed an inhomogeneous approach in the literature, e.g., by Yan et al. [60], and should not be confused with the assumption of a homogeneous mixture of two phases in each computational cell. Using separate velocity fields is indispensable to predict a phase separation and air accumulations, and its omission leads to an underestimation of head drop, as has been demonstrated by Yan et al. [60] and shown in Figure 3. Mass transfer between the phases is neglected in Equation (1). Note that the interface compression term in Equation (1) is also omitted in the pure EE2F model, but is noted here to facilitate the summary of the VoF algorithm further below. The diffusive mass transport is neglected, an assumption that is justified by a Peclet number higher than 100 for the investigated pump flow. The volume fraction, density, and velocity of phase ϕ are symbolized by α ϕ , ρ ϕ , and c ϕ , respectively, and R eff ϕ represents the stress tensor combining Reynolds (turbulent) and viscous stresses. Note that indicating Reynolds averaged quantities is omitted here for convenience. The sets of momentum equations are coupled by interfacial momentum transfer forces. The interfacial momentum transfer terms, i.e., drag and virtual mass force, are represented by M ϕ , while M s,ϕ is the surface tension force. Regarding M ϕ , the drag force M a,D and the virtual mass force M a,vm dominate all other interfacial forces [37], so that only these forces are considered here.
In the pure EE2F model, M s,ϕ is set to zero. The drag coefficient is evaluated according to Schiller and Naumann [136], and the virtual mass force coefficient is set constant to 0.5, according to Refs. [137,138]. A separate turbulence field is applied for both phases, i.e., balance equations for turbulence kinetic energy k ϕ and turbulence specific dissipation ω ϕ are solved for both phases. Both phases share the same pressure field, which means that a single mixture continuity equation, in terms of a pressure Poisson equation, is solved.
Thus, in the EE2F approach, and also in the H2P approach (the latter being presented further below), two sets of conservation equations are solved, which means one set for each phase ϕ. In the VoF approach, typically, one single set of conservation equations is solved, formally obtained by summing up the momentum equations of both phases (Equation (2)) and setting ϕ to mixture properties, assuming a common velocity and pressure field for both phases in the entire flow domain. A mass balance equation is solved for water, while the volume fraction of the air phase is calculated by α a = 1 − α w . The water/air interface is sharply resolved by an interface sharpening process [82] (interface compression term in Equation (1)) that counteracts the numerical diffusion and which is explicated in more detail in Section 3.2. For VoF, the surface tension force M s,ϕ is considered as a continuum force, according to Brackbill et al. [139], with a constant surface tension factor between water and air of 0.07275 (see also Equation (5) in Section 3.2), and M ϕ is set to zero.
For both EE2F and VoF schemes, both phases are treated as incompressible fluids with a density of ρ w = 998 kg/m 3 and ρ a = 1.185 kg/m 3 , which was justified, for example, by Hundshagen et al. [56][57][58]74,76]. The results hardly deviate when the air phase is treated compressible or incompressible, which was shown by Müller et al. [48] for the EE2F approach, and which is also assumed for the VoF and H2P approaches. The flow is assumed to be isothermal, so that a constant temperature of 25 • C was chosen, omitting the energy conservation equation. Both EE2F and VoF schemes were implemented in OpenFOAM in different solvers, which were customized, where applicable, for example, for moving grid capabilities [74]. The EE2F and VoF implementations are based on OpenFOAM version 6 within the reactingMultiphaseEulerFoam and interFoam solver, respectively.

Hybrid Two-Phase (H2P) Approach
Recently, we presented an enhanced EE2F model with the capability of local phase interface resolution [76,140], so it locally works like a VoF scheme. The functionality of this hybrid scheme is best demonstrated by the interface compression term in Equation (1). To enable an interface resolution, according to Wardle and Weller [82], a correction flux at the phase interface is introduced in terms of interface compression. This correction flux is added in Equation (1) in the interface compression term by introducing the compression velocity c c and the volume fraction term α ϕ 1 − α ϕ . The compression velocity c c counteracts the numerical diffusion at the phase interface to prohibit interface smearing and to maintain a sharp interface. It is estimated by the blending function C α , and the magnitude of the mixture velocity |c m | yielding: The ratio ∇α/|∇α| affects that the interface compression term acts orthogonally to the phase interface, while α ϕ 1 − α ϕ in Equation (1) ensures that the interface compression term essentially acts in proximity of the phase interface. The mixture velocity c m is calculated by a volume fraction-average of water and air velocity. For the evaluation of C α in Equation (3), a blending function adopted from Hänsch et al. [77] is implemented: In the limiting case C α = 0, the original EE2F formulation is retained, while with C α = 1, a VOF-like interface resolution is applied. At this point, it is important to mention that for pure VoF schemes, as described in Section 3.1, no blending for C α is applied. Hence, a constant value of C α = 1 is set [141]. A second difference between pure VoF and H2P schemes is that, for VoF, only one set of mass and momentum equations is solved, while for H2P, one mass and momentum equation is solved for each phase.
The surface tension force M s,ϕ in Equation (2) is considered as a continuum force, according to Brackbill et al. [139], and reads: The surface tension factor for a water/air system and the interface curvature are σ = 0.07275 and κ ϕ = −∇ ∇α ϕ |∇α ϕ | [82]. M s,ϕ is also blended with C α according to Equation (5), so that the surface tension force is only active in regions where C α differs from the value zero, which means that the phase interface resolving mode is activated.
In the dispersed phase (subscript ϕ = a), dispersed air bubbles are mixed within a continuous water phase (subscript ϕ = w). The momentum equations (Equation (2)) for each of the two phases ϕ are coupled by the forces M ϕ . As already noted for the EE2F model in Section 3.1, the drag force M a,D and the virtual mass force M a,vm dominate all other interfacial forces [62], so only these forces were considered here. While the expression for M a,vm is presented elsewhere [76], we note, for the drag force: d B is the bubble diameter and is explained further below in Section 3.3. In Equation (6), a second blending function C β is used, adopted from Ref. [142] and customized for the H2P solver in Ref. [76]. C β reads: This drag blending affects that the drag formulation is only active in the dispersed part of the flow and turned off within air accumulations, which is equivalent to the assumption that, inside air accumulations, only a continuous air phase is present while there is no dispersed water in the form of droplets. We set the constants to b min = 10 −4 and b max = 0.8. Thus, air accumulations were assumed to be present in a range of air volume fraction between 0.8 and 1.0.
In terms of Equations (4) and (7), different blending functions C α and C β were used for interface compression and drag. The particular functional and the constants of C α and C β were chosen and implemented to achieve stability for the complex centrifugal pump flow. Using the same blending function for interface compression and drag might be more consistent, as has been proposed, e.g., by De Santis et al. [73,143] and Mathur et al. [84].
To be more precise, we used two different blending functions in terms of C α and C β since we encountered stability issues when using C α for drag as well. This is traced back to the complex grid and flow field of the centrifugal pump. Therefore, for the sake of simplicity, we limited the dispersed drag force to a low finite value by clipping C β to 10 −4 in regions of air accumulations, i.e., α a > 0.8. A further optimization of C α and C β is left for future work.
The hybrid model described above was implemented in OpenFOAM versions 6 and 9 as an extension of the EE2F model within the reactingMultiphaseEulerFoam and multi-phaseEulerFoam solver. It was the first time it was adopted for centrifugal pump flow and certainly leaves potential for further improvements. In the literature, more sophisticated hybrid multiphase approaches have been reported, such as the generalized multiphase modeling approach (GEMMA) of De Santis et al. and Colombo et al. [73,143,144] and the hybrid multi-scale model of Meller et al. [79]. The GEMMA model fundamentally deviates from the one described here by including an interfacial drag formulation for large-scale interfaces, which was first proposed by Marschall [75] within the EE2F framework. Meller et al. [79] coupled two continuous phases by applying a drag model of Štrubelj anď Tiselj [145]. The GEMMA model also introduced an interface resolution criterion, which was first described in Ref. [83] and is based on the ratio of local cell size V c to the Sauter mean diameter of the dispersed phase d 32 . The interphase compression term is activated when the criterion d 32 /V c > Γ is fulfilled. The constant Γ is set to a value of ten [73,83,143], which may be case-dependent. Thus, the GEMMA model deviates from our approach also by including the droplet size in the blending function formulation. Both the GEMMA model and the ansatz of Meller et al. [79] should be tested for centrifugal pump flow in future studies.

Population Balance Model (PBM)
Optical measurement data by 69], and described in Part B of this review paper series, suggest a polydisperse bubble size distribution within the pump, which means that dispersed bubbles have a spectrum of properties. As summarized in Section 2.2, several authors [59][60][61][62][63][64][65] coupled an EE2F method with a population balance model (PBM) [146,147]. Thereby, an additional transport equation for the specific number density n B (d B ) of bubbles is solved. Beyond the transport of n B (d B ), birth and death rates due to coalescence and breakup of bubbles with diameter d B are treated by source and sink terms. We adopted the OpenFOAM implementation of the PBM of Lehnigk et al. [148] and merge it into the dispersed part of the hybrid H2P solver. By adopting the approach of Kumar & Ramkrishna [149], the population balance equation is discretized in definite size groups, leading to the class method of Lo [146], also referred to as the multi-size group (MUSIG) model [147,150]. The coalescence kernel of Prince & Blanch [151] and the breakup kernel of Luo & Svendsen [152] were chosen. The coalescence model was extended by following Liao et al. [153], resulting in a formulation similar to the one suggested by Stel et al. [61] for the simulation of a liquid/gas centrifugal rotor flow. The kernels comprise empirical parameters, which were not adapted but retained at the original values suggested by Refs. [151,152]. From the resulting bubble size distribution, a statistical mean, in terms of d 32 , is evaluated and enters the disperse part of the H2P, as, for example, in Equation (6).

Simulation Method
Preferably, customized versions of OpenFOAM were used as a software platform for our CFD simulations of multiphase flow in centrifugal pumps or pump mixers [74,76,154]. The implementation employed the Pressure Implicit with Splitting of Operators (PISO) algorithm of Jasak [155]. The PISO scheme works like a quasi-explicit solver and has turned out to be particularly efficient for small time steps, corresponding to a convective Courant-Friedrich-Lewy number not exceeding 0.5 [155]. For larger time steps, non-linear iteration loops and an under-relaxation were performed, resulting in the well-established PIMPLE algorithm, a combination of the PISO and Semi-Implicit Pressure Method for Pressure Linked Equations (SIMPLE) scheme, introduced by [156]. Discretization of convective terms was performed by second order Total Variation Diminshing (TVD) schemes. Particular care was taken when turbulence scale-resolving models, such as SAS, were used. Customized hybrid discretization schemes are used for adequate scale resolution, as described elsewhere [76,133].

Test Case and Simulation Setup
A research pump with a semi-open impeller and 2D-blading, i.e., six cylindrical (nontwisted) blades, was investigated. The pump head H = ∆p/(ρ w g) + ∆z was evaluated by the static pressure and height difference ∆p and ∆z, respectively, between suction and pressure port. Beyond pump head and flow rate measurements, resulting in performance maps, according to Figure 1, high-speed visualization of the two-phase flow structures was performed by means of the scattered-light technique. Therefore, the im-peller and the casing were manufactured from transparent acrylic glass to obtain optical access. Experiments were performed by Mansour et al. [26,27] at reduced rotational speed n = 650 1/min to avoid damage of acrylic components. Details of the experiments are provided elsewhere [24][25][26][27]39,56,58,76], as well as in the accompanying experimental Part B of this review paper series.
The simulation domain is depicted in Figure 6a. A preferably realistic image of the real geometry was aimed at, so that the domain includes the impeller, suction and pressure pipe, volute casing, and side chamber. To capture the rotor-stator interaction and to account for any unsteady effects, unsteady pump simulations in the absolute frame of reference and with moving impeller mesh were performed, and the arbitrary mesh interface of Farrell and Maddison [157] was used. The SAS turbulence model, as described in Section 2.5, was applied. Velocity and pressure were prescribed as Dirichlet boundary conditions at the inlet and outlet of the computational domain, respectively. Dispersed air bubbles were assumed at the inlet, which were prescribed in terms of a Dirichlet condition for α ϕ and d B = 0.5 mm. The operation points considered in the simulation are indicated above in Figure 1. By having this range of operation conditions, all flow regimes from bubbly flow via agglomerated bubbles and alternating pocket flow to pure pocket flow were covered.
Particular care was taken for computational grid generation. By a grid refinement, first of all, the accuracy of the discretization scheme rose. At the same time, grid refinement also affected the turbulence and two-phase models. On the one hand, the turbulence scale-resolution capability of the SAS was enhanced, since turbulent fluctuations could be resolved down to the local grid resolution. On the other hand, successively finer void structures were resolved by the H2P solver, which transitioned towards a VoF-like scheme with phase-interface resolving capabilities. Thus, an inherent grid dependence of the results could not be avoided. What is more, the complex geometry of circumferential pumps required special care with regard to grid quality. For example, the meshing of round trailing edges or narrow gaps in semi-open impellers is particularly demanding [56][57][58]. Despite the considerable effort for grid generation, we prefer hexahedral cell grids wherever applicable due to their higher accuracy than, for example, triangular cells. Peculiarities of grid generation are provided in Refs. [56][57][58]76] and are not all repeated here. A grid study was performed, and the preferred computational grid, for which the results are presented further below, had about 16 Mio cells. A coarser version of the grid is shown in Figure 6b-d. Visually, smoothness and orthogonality appear to be poor. However, non-dimensional grid quality criteria, in terms of maximum non-orthogonality and minimum face volume ratio, were adopted from the OpenFOAM nomenclature [155] and indicated a fairly high quality, taking into account the complexity of the pump geometry. It was ensured that the conclusions drawn from the simulation results were not affected by grid dependence. Details of the simulation setup and the grid study are presented elsewhere [76].

Selected Results
Selected results obtained by the H2P model in combination with the PBM are presented here. In Figure 7, the pump head drop by a rise of the inlet gas volume fraction ε is demonstrated. Albeit there were remaining differences to the experimental result, the trend was captured well at the first attempt. This is remarkable since, by means of a monodisperse ansatz, the predicted head level in the simulation can immediately be tuned to the experimental data by the prescribed bubble size d B [56][57][58]. In fact, for the applied grid resolution, a value of d B = 2.0 mm also yielded a reasonable head drop curve. Obviously, by using the PBM, a reasonable bubble size distribution was obtained without the need for tuning, which affected proper head drop characteristics. What is more, solver robustness was enhanced by the PBM, compared to the use of a constant value of d B , which we trace back to the fact that, in flow regions with small cell size, shear and, thus, bubble breakup were enhanced, and bubbles larger than the computational cells occurred less often.
The air accumulation is visualized in Figure 8a by means of the ensemble-averaged experimental results of a gray-scale analysis. With rising ε, air successively accumulated at the blade suction side, which, in turn, caused a head drop. Air accumulation was observable for ε = 1% neither in the experimental nor in the simulation result. For ε = 3%, the air started to agglomerate at the suction side of the blade. For ε = 4%, the tendency of water and air to separate within the blade channel further increased. Air accumulations formed at the blade suction side near the leading edge. It is interesting that the accumulations were observable at every second blade in the simulation result. This behavior, in fact, reflects the alternating pocket flow regime, as discussed in detail by Hundshagen et al. [76]. For ε = 5%, the air accumulations were attached at each blade and increased in size, which was, again, in good agreement with the experimentally observed trend. Summarizing, the size of air accumulations was captured well by the simulations, according to Figure 8b. Again, it should be pointed out that these accumulations developed without further tuning of the simulation parameters in terms of, for example, coalescence and breakup parameters. Thus, from ε = 1% to ε = 5%, a transition from the bubbly flow regime over the agglomerated bubbles flow regime to the pocket flow regime was observed. In Figure 8c, the standard deviation of the air volume fraction, in terms of α a,RMS , is depicted. In the agglomerated bubbles flow regime for ε = 3%, the bubble clustering was highly unsteady since the highest values of α a,RMS were observable in regions having the highest values of α a . In the alternating pocket flow regime at ε = 4%, the air pocket, attached only at every second blade, was steady according to the low α a,RMS level within the attachment zone. At the other blades, local maxima of α a,RMS were present at the same blade location (as in the agglomerated bubbles flow regime at every blade), which, in fact, reflected the unsteady behavior of the alternating pocket flow regime. With a rising value of ε from ε = 3% to ε = 5%, the values of α a,RMS successively decreased in the region of air accumulations, which meant that the accumulations became more and more steady. The highest values of α a,RMS were present in the wake of the accumulations, which corresponded to the experimental observations of an unsteady wake and underlined the benefit of using a scale-resolving turbulence model.
In Figure 9a, the Sauter mean diameter d 32 of the bubbles is shown, and is, again, evaluated from the simulation results at mid-span. A cross-check to Figure 8b reveals that the location of large bubble diameter correlated with air accumulation zones, which meant that preferably large bubbles accumulated. A more detailed look at the coalescence kernel outcome (not shown here) showed high kernel activity in these regions, so that a self-energizing effect was present, leading to enforced growth of the accumulation. A look at the blending function C β , according to Equation (7), which is illustrated in Figure 9b, by an isosurface C β = 0.8 revealed that the H2P solver worked in a VoF-like mode in the accumulation regions. It is important to point out that by the transition of the solver towards VoF within the accumulation zones, the dispersed nature of the mixture had locally faded out, and the physical meaning of d B and, thus, d 32 degenerated. In fact, the diameter mutated into a purely numerical variable with a diminishing feedback to the local flow solution. It is interesting to note that no special means had to be applied to enhance solver stability in these flow situations. Thus, the air accumulations were treated in a reasonable way by a combination of an EE2F and VoF approach in terms of the H2P solver.
Bubbly flow Agglomareted bubbles flow Alternating pocket flow Pocket flow

Conclusions
We have pointed out the limitations of recent approaches for the 3D simulation of liquid/gas flow in centrifugal pumps. On the one hand, the resolution of all phase interfaces by Volume-of-Fluid (VoF) methods is not feasible, so that this kind of approach is confined to the resolution of large coherent air accumulations. On the other hand, Euler-Euler Two-Fluid (EE2F) models treat dispersed bubbles on a subgrid level. Inspired by CFD approaches resulting from chemical and process engineering research [73,[77][78][79][82][83][84]143,144,158,159], a hybrid method (H2P) is suggested, which is based on an EE2F scheme in most regions of the flow field, and transitions to VoF-like resolution, when air accumulations are detected. The H2P is combined with a population balance model (PBM) to capture bubble size distribution in flow regions with dispersed bubbles, i.e., where the H2P works in the pure EE2F mode. It was shown that coalescence activities correlate with air accumulations, which means that, in accumulation zones, bubbles grow by coalescence, leading to large air pockets. By means of this approach, the transition of flow regimes from bubbly to pocket flow detected in the experiments could, for the first time, be reproduced in 3D flow simulations. The importance of a scale-resolving turbulence approach to capture unsteady wakes downstream of the air pockets was shown. In spite of this quite promising first application of the H2P scheme to centrifugal pumps [76], and rather generic diffusor flows [140], improvements, for example in regard to blending functions, are suggested for future work. The PBM, particularly the breakup and coalescence kernels, have, so far, simply been adopted from completely different applications, such as bubble columns, and should be tailored for centrifugal pump flow applications in the future. In future studies, it should also be assessed if the H2P method, in fact, converges to a VoF solution with excessive grid refinement.

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

Nomenclature and Abbreviations
The following nomenclature is used in this manuscript: Roman characters a (min,max,B) Model constants for C α (-) b (min,max) Model constants for C β (-) c Velocity vector (m/s)