3D simulation of gas-laden liquid flows in centrifugal pumps and the assessment of two-fluid CFD methods

An assessment of a two-fluid model assuming a continuous liquid and a dispersed gas phase for 3D computational fluid dynamics (CFD) simulations of gas/liquid flow in a centrifugal research pump is performed. A monodisperse two-fluid model, in conjunction with a statistical eddy-viscosity turbulence model, is utilized. By a comprehensive measurement database, a thorough assessment of model inaccuracies is enabled. The results on a horizontal diffuser flow reveal that the turbulence model is one main limitation of simulation accuracy for gas/liquid flows. Regarding pump flows, distinctions of single-phase and two-phase flow in a closed and semi-open impeller are figured out. Even single-phase flow simulations reveal challenging requirements on a high spatial resolution, e.g., of the rounded blade trailing edge and the tip clearance gap flow. In two-phase pump operation, gas accumulations lead to coherent gas pockets that are predicted partly at wrong locations within the blade channel. At best, a qualitative prediction of gas accumulations and the head drop towards increasing inlet gas volume fractions (IGVF) can be obtained. One main limitation of two-fluid methods for pump flow is figured out in terms of the violation of the dilute, disperse phase assumption due to locally high disperse phase loading within coherent gas accumulations. In these circumstances, bubble population models do not appear beneficial compared to a monodisperse bubble distribution. Volume-of-Fluid (VOF) methods may be utilized to capture the phase interface at large accumulated gas cavities, requiring a high spatial resolution. Thus, a hybrid model, i.e., a dispersed phase two-fluid model including polydispersity for flow regions with a dilute gas phase, should be combined with an interphase capturing model, e.g., in terms of VOF. This hybrid model, together with scale-resolving turbulence models, seems to be indispensable for a quantitative two-phase pump performance prediction.


Introduction
Centrifugal pumps are employed in many technical applications, especially when reliable and flexible pump performance is needed. Often, their conveying task is not limited to pure liquids, since also mixtures, e.g., liquid and gas, must be delivered. This is the case when economically expensive phase separation processes are avoided, e.g., in oil, see, e.g., Caridad et al. (2008), or wastewater industry, see, e.g., Cappelino et al. (1992). However, the transport of liquidgas mixtures may cause a significant drop of pump head even at very low inlet gas volume fractions (IGVF) of about 1% (Murakami et al., 1971). This effect is usually more dominant in part-load and overload than in optimal load (Cappelino et al., 1992).
Several experimental studies pointed out that the drop of the pump head is caused by the separation of air and water in the blade channel (Murakami et al., 1971;Furukawa et al., 1988;Cappelino et al., 1992;Tillack, 1998). This process is forced by the Coriolis force inside the impeller of centrifugal pumps, the relative motion between water and air, the pressure gradients in cross-flow direction, and the density difference between the phases (Sato et al., 1996). Monte Verde et al. (2017), Stel et al. (2020a), andV Vol. 3, No. 3, 2021, 186-207 Experimental andComputational Multiphase Flow https://doi.org/10.1007/s42757-020-0080-4 3D simulation of gas-laden liquid flows in centrifugal pumps and the assessment of two-fluid CFD methods  Mansour et al. (2018bMansour et al. ( , 2018c categorized the two-phase air-water flow patterns in centrifugal pumps in five groups, i.e., bubbly flow, agglomerated bubble flow, pocket flow, alternating pocket flow, and segregated flow. These groups are ordered by increasing IGVF and bubble interaction, particularly coalescence effects. Air separation processes become stronger with increasing air bubble size. Thus, a reliable two-phase flow pump performance requires small gas bubbles inside the impeller, which can be forced by employing semi-open impellers. The leakage flow over the blade tip enforces secondary flow, which breaks-up big bubbles and flushes the air out of the blade channel (Mansour et al., 2018c). A review of the two-phase performance of various centrifugal pump designs is presented by Jiang et al. (2019), mainly based on experimental results. High sophisticated measurement methods in terms of X-ray and gamma-ray tomography are utilized by Schäfer et al. (2015Schäfer et al. ( , 2020 and Neumann et al. (2016) to characterize both the time-averaged and time-dependent void structures in a transparent research pump. Although the measurement of the spatially and temporally resolved gas-liquid mixture flow field demands sophisticated measurement techniques and is thus associated with a high effort, the level of detail of measurements is limited, e.g., by the flow accessibility and the temporal and spatial resolution of the probes and data acquisition systems. On the other hand, 3D simulations provide a high amount of local flow information. In Zhu and Zhang (2018), a recent literature review concerning gas-liquid flows in centrifugal pumps is presented that also summarizes the progress of recent computational fluid dynamics (CFD) efforts. Mostly, two-fluid simulation methods in terms of Euler-Euler methods are employed. Two of the first numerical investigations of two-phase flow in centrifugal pumps are carried out by Pak and Lee (1998) and Minemura and Uchiyama (1993). Both studies pointed out that bubbles move to the blade tip and accumulate on the blade's suction side. It was proven that the head drop at increasing IGVF is caused by phase separation and large gas accumulation zones. Müller et al. (2015) employed a monodisperse Eulerian-Eulerian two-phase model to evaluate the head drop in a single-channel simulation domain of a centrifugal research pump with IGVF variation. A good agreement with the measured head drop was obtained up to IGVF = 3%. At higher IGVF, even a qualitative prediction could not be obtained, confirming the conclusion of several previous studies, e.g., Minemura and Uchiyama (1993), Caridad et al. (2008), and Yu et al. (2012). A significant influence of the computational grid quality on the location and size of bubble accumulation was found by Müller et al. (2016). Numerical simulation studies on electric submersible pumps (ESP) were carried out to provide a better understanding of flow parameters and phase interaction processes (Rutter et al., 2017;Zhu and Zhang, 2017). Both studies pointed out that an adjustment of the gaseous phase bubble size at high IGVF has to be done to capture a realistic flow behavior with increasing IGVF. Dupoiron (2018) reveals that the maximum bubble diameter directly influences the relative magnitude of the drag force in ESP. Si et al. (2017Si et al. ( , 2018Si et al. ( , 2020 presented experimental and numerical results of a 3D-bladed centrifugal pump in two-phase air/water flow. A further numerical investigation of this topic was undertaken by Zhang et al. (2017) to enable a pump impeller design optimization.
The performance of Euler-Euler methods in two-phase flow simulations is very sensitive to the modeling of the inter-phase momentum transfer and the bubble size of the dispersed phase. On the one hand, significant deviations may occur between experimental data and simulation results when these models are inappropriately chosen, as, e.g., shown by Caridad and Kenyery (2004) for ESP, or Yamoah et al. (2015) for vertical pipe flows. On the other hand, Stel et al. (2020b) recently experienced a satisfactory agreement between experimentally and numerically predicted gas accumulation zones in a centrifugal rotor by utilizing a population balance method for the air bubble diameter. These numerical studies emphasize that there is up to now no standard modeling of inter-phase momentum transfer or bubble size, which can be applied in general for different two-phase flows (Wang and Yao, 2016).
In these exemplary cited simulation studies, two-fluid simulation methods in terms of an Euler-Euler method were employed, i.e., a dispersed gas phase is assumed to be mixed within a continuous liquid phase. This assumption is strictly valid only in bubbly flow where a homogeneous mixture of distinct small bubbles within each computational cell is present, and bubbles are significantly smaller than the cell size, which is equivalent to the assumption that the computational grid does not resolve bubble-liquid interfaces. Distinctive velocity fields are considered for each phase by the solution of phase momentum equations, which enables a de-mixing of phases. An explicit limitation of multi-fluid methods is a situation with a locally too high loading of the dispersed phase within the continuous phase that may occur either in flow regions with high gas accumulation or for fine computational grids, which will be discussed later in Section 4. In the studies cited above, where multi-fluid methods are employed for centrifugal pump gas/liquid flow, these limitations of the dispersed multi-fluid assumptions have not been foregrounded.
Another family of multiphase methods in terms of Volume-of-Fluid (VOF) methods aims at a direct resolution of the liquid-gas interface and void structures down to the available grid limit with the assumption of a common velocity field at the phase interface. Pineda et al. (2016) and Zhu et al. (2019) investigated the two-phase flow patterns in ESP using VOF methods and compared them to experimental data. Mansour et al. (2020aMansour et al. ( , 2020b, Kopparthy et al. (2020), andParikh et al. (2020) employed VOF methods on the same test cases, i.e., diffuser and research pump flow, investigated in the present study, figured out important physical mechanisms of attenuating pump performance, and the beneficial effect of an inducer (Mansour et al., 2020a(Mansour et al., , 2020b. However, it remains unclear if the spatial resolution limit allows a sufficient resolution of void structures within the VOF simulations.
The aim of this paper is the assessment of two-fluid simulation methods for centrifugal pump mixture flow by direct comparison with corresponding experiments. Although bubble interaction, i.e., break-up and coalescence, is assumed to influence the air accumulation at the impeller blades and thus the two-phase flow performance of pumps, the present numerical study is essentially confined to a monodisperse two-phase model to assess the bubble diameter impact separately. A further central assumption is the incompressibility of gas, which is justified by the low head of the research pump and was verified in preliminary simulations (Müller et al., 2015). To investigate the influence of different impellers on the two-phase flow, a semi-open (OI) and a closed (CI) impeller, whose blades are geometrically identical, are investigated. Pump characteristics such as pump head in single-and two-phase flow are evaluated and compared. Gas accumulation within the impellers predicted by the simulations is compared to corresponding optical measurements. The experimental results presented in this paper are adopted from Mansour et al. (2018aMansour et al. ( , 2018bMansour et al. ( , 2018c and Hundshagen et al. (2019aHundshagen et al. ( , 2019b. Since the planar pump design has been optimized for optical access of the two-phase flow, it enables a unique validation database. Before we study the pump flow, an investigation of a simpler test case in terms of diffuser flow is presented.
The paper is organized as follows: the test cases in terms of diffuser and pump flow, and the experimental data are presented in Section 2. Although this paper focuses on the flow simulation by CFD methods, we provide a thorough review of measurement methods and results to enable a critical assessment of the simulation. In Section 3, a description of the simulation method is provided. The results of single-phase, as well as gas-laden liquid diffuser and pump flow simulations, are presented in Section 4, together with an assessment of the limitations of two-fluid CFD methods. The paper is closed with our conclusions in Section 5.

Horizontal planar diffuser flow
As a first test case to study the transition from dispersed bubbly flow to a continuous air accumulation, a diffuser flow is considered. The test case and measurement data have been presented in detail by Mansour et al. (2018a) and are briefly summarized here.
A turbulent two-phase air/water flow through a horizontal, diverging channel has been investigated. The mixture flows from a rectangular channel of 40 mm × 44 mm through a diffuser to a rectangular channel of 100 mm × 44 mm, corresponding to a hydraulic diameter (dhu) ratio of 1.45.
Further geometric details, as well as the position of the pressure sensors P1 to P8, are illustrated in Fig. 1.
Due to an increasing opening angle in the flow direction, flow separation is enforced in the diffuser within the flow conditions under consideration, which is associated with a coherent air accumulation in the separation zone. The superficial Reynolds number of water and air phases has been varied in the range of Re w = 50130-87730, and Re a = 3-18.5, respectively. The results show large gas accumulation even at a meager air volume fraction of 0.05%, which significantly affects the velocity field and the pressure recovery of the diffuser. The impact of Reynolds numbers on the gas pocket size is studied. An increase of Re a is accompanied by a continuously enhanced gas accumulation, while for high Re w values, it attenuated again, what is associated with the break-up and shedding of gas pockets due to enhanced turbulence that is transported downstream. Albeit much more straightforward, the diffuser flow is similar to a centrifugal pump flow in terms of a positive pressure gradient and the occurrence of separation, e.g., at off-design operation. The Coriolis force perpendicular to the mean flow direction is mimicked by gravitation. However, since no rotating parts are involved, the diffuser configuration allows a far more accurate and complete experimental characterization. Therefore, this study provides a first step toward understanding the complex flow patterns occurring into centrifugal pumps transporting gas/liquid mixtures and provides a database for CFD methods that is picked-up here for the validation of the two-fluid simulation model. Selected measurement results are presented together with simulation results in Section 4.1.

Centrifugal pump
As for the diffuser flow, only a brief summary of the experimental setup is presented here, while details can be found in Mansour et al. (2018c). Whereas in Mansour et al. (2018c), mainly characteristic curves and associated flow patterns are discussed, in Hundshagen et al. (2019b) further details of flow visualizations have been presented that are also summarized here.

Experimental setup
A research centrifugal pump has been installed at the Lab. of Fluid Dynamics and Technical Flows at the University of Magdeburg. The entire pump casing and the impellers are manufactured from acrylic glass to allow a precise flow visualization. By employing a high-speed camera, LED lamps and a scattered light technique, the gas-liquid interfaces appear brighter than the single-phase liquid due to light reflections, making the bubbles and the accumulated gas cavities clearly observable. To obtain the time-averaged size of the accumulated gas zones in the impeller passages, ensemble-averaging, i.e., cyclic image recording, is used with the same frequency of impeller rotation, always acquiring one image per rotation. Afterwards, an average image is deduced from 220 instantaneous images of the flow, and gas accumulations are captured and colored by a grayscale analysis of the average image. A similar time-averaging technique has been already employed previously by Mansour et al. (2018a) to analyze the size of the accumulated gas in two-phase air/water flows in the plane diffuser, as described in Section 2.1 above.
The pump front view, the OI, and the CI are shown in Figs. 2(a), 2(b), and 2(c), respectively, while Fig. 3 shows a 3D view of the test-rig. Pump data are listed in Table 1. Six cylindrical, non-twisted blades, which are geometrically identical in both impellers, are applied to maximize the flow visualization ability. A rotational speed of n = 650 min − 1 -corresponding to a specific speed of approximately n q = 21 min − 1 -was kept for all investigations since at higher rotational speeds, strong system vibrations occurred in the experiments, which could destroy the acrylic parts. The recording of pump characteristics is described in detail in Mansour et al. (2018bMansour et al. ( , 2018c. The air is injected in the center of the suction pipe through a nozzle (as marked by "Air injection" in Fig. 3). The pressure is measured across the pump by two pressure sensors and is used to calculate the air volume flow rate and volume fraction. The temperature is measured near the injected air. A frequency controller is used to set the rotational speed. A torque The pump efficiency is calculated according to

Measurement results
The key parameters of the experiment have been presented in Section 2.2.1. Some key findings are presented here, and more details can be found in Mansour et al. (2018c) and Hundshagen et al. (2019b). Further experimental results are presented together with simulation results in Section 4.2. Figure 4 shows the single-phase characteristics for both impellers in terms of pump head and pump efficiency, respectively. The values are normalized by the optimal pump head and nominal flow rate of the CI. Due to the tipleakage flow, the pump head and efficiency of the OI are slightly lower than those of the CI, particularly in overload operation. Notably, the CI efficiency is approximately 1%-3% higher than that of the OI.
In Fig. 5, the head degradation for a gradual increase of IGVF is presented. For the CI, a continuous drop occurs up to IGVF of about 6%, which is slightly more pronounced at  overload, i.e., Q w /Q opt = 1.25. For the OI, up to an IGVF of 3%, only a slight head drop is discernible so that H does hardly deviate from single-phase flow. For higher values of IGVF, a rather abrupt head drop occurs due to massive air accumulations along the blades, preventing the pump from a proper transport of the mixture. Therefore, for IGVF = 1%-3%, the OI seems beneficial, while the CI is more effective for IGVF = 4%-6%.
Based on flow visualization illustrated in Fig. 6, the two-phase flow pattern in the impeller passages is classified into five regimes:  Bubbly flow-Bubbles are dispersed everywhere in the impeller without significant interaction between them.  Agglomerated flow-Bubbles interaction starts by coalescence, forming bigger bubbles near the blades.  Alternating pocket flow-A gas pocket near the blade entrance part with unsteady properties appears (strong oscillations, appearing/ disappearing).  Pocket flow-A big gas pocket (cavity) can steadily stand near the blades.  Segregated flow-The gas pocket is connected throughout the blades till the outer impeller diameter. The bubbly flow regime (see Fig. 6(a)) occurs for both impellers mainly at overload operation and low inlet gas volume fraction, IGVF < 2%, and extends to a broader range of flow conditions for the OI, as a result of the higher shear rates exerted on the blades because of the secondary flow. The alternating pocket flow regimes (see Fig. 6(c)) appear at part-load as well as overload operation and is smaller in the OI, corresponding to lower flow instabilities compared to the CI. In the CI, bubble-bubble interaction is stronger; hence, the agglomerated flow regime (see Fig.  6(b)) was found to appear more likely than in the OI, at IGVF < 2% in nominal and part load. The pocket flow regime (see Fig. 6(d)) appears in the CI at lower IGVF, starting from IGVF ≈ 3% and showing no remarkable   change in the pump performance. On the other hand, as soon as the pocket appears in the OI (starting from IGVF ≈ 4%), an abrupt reduction in the pump performance can be seen in particular at nominal and for overload operation (see Fig. 5(b)). As the IGVF is increased from 4% to 7%, the OI performance shows several discontinuities. When surging occurs in either of the impellers, the pump fluctuates between two different operational points. Therefore, the size of the accumulated gas is strongly unstable as a result of the oscillating pumping behavior. In the CI, the gas cavity can be large enough to cover the whole length of the blade, reaching down to the outer diameter of the CI. This two-phase flow pattern occurs for IGVF > 9% and corresponds to a segregated regime ( Fig. 6(e)). This regime is not observed in the OI because of the intense bubble break-up that occurs near the impeller outlet. This exemplary discussed image series provides physical reasoning for the characteristics shown in Fig. 5, later in Fig.  18 and for the performance maps that are presented in detail in Mansour et al. (2018c). Here we do not revisit the performance maps but focus our attention on particular operation points that are proposed in Section 4.2 for validation of CFD. The particular two-fluid CFD method is presented next.

Flow solver
The numerical method is described in detail in Hundshagen et al. (2019b) and is briefly summarized here. The commercial CFD solver ANSYS CFX 18.0 ① is utilized for the solution of the unsteady Reynolds-averaged Navier-Stokes equations. If not otherwise mentioned, the statistical k-ω shear stress transport (SST) turbulence model by Menter (1994) and Menter et al. (2003)-termed SST model in what follows-in conjunction with an automatic wall treatment of Esch and others (Menter and Esch, 2001;Vieser et al., 2002) is applied. In Sections 4.1.2 and 4.2.3.2, a turbulence model variation is performed and described therein. An inhomogeneous monodisperse Eulerian-Eulerian two-fluid model, i.e., a dispersed bubbly gas phase is assumed to be mixed within a continuous liquid phase, and a constant bubble diameter is chosen. Thus, bubble interaction and bubble size distributions are neglected. We term this model class disperse two-fluid model in what follows. Note that we use the term two-fluid model instead of multi-fluid model since only one gas and one liquid phase are present in our particular gas/liquid application. An Eulerian approach for the dispersed phase is preferred to the Lagrangian approach since a more moderate grid dependence, and better statistical convergence of the Eulerian approach is expected. Preliminary tests showed that the results hardly deviate when using compressible or incompressible modeling of the gaseous phase. Therefore, both phases are treated as incompressible, which leads to a significant stabilization of the CFD solver (Müller et al., 2015). To enable phase separation, a separate momentum balance equation for each phase, water and air, is solved and yields a separate velocity field for each phase. Preliminary simulations showed that a homogenous treatment of turbulence and pressure fields, i.e., both phases share the same turbulence and pressure fields, is sufficiently accurate and, at the same time, numerically stable. Double-precision accuracy of floating-point numbers and second-order discretization in space and time is applied. Momentum and continuity equations are solved with a coupled procedure. The volume fraction equation is solved in a segregated manner for both phases, due to a better convergence than with coupled volume fraction solution procedure, see Müller et al. (2015Müller et al. ( , 2016. Because the drag force is assumed to dominate all other interfacial forces in centrifugal pumps, all non-drag forces are neglected. The drag between the dispersed air and the continuous water is calculated by the non-dimensional drag coefficient c D , which is modeled by the Schiller-Naumann drag law (Schiller and Naumann, 1933). This drag model is valid for perfectly spherical bubbles, an assumption which is assumed to be valid for any bubble and void cavity size in this study.
Albeit most of the simulations are performed with a monodisperse bubble distribution and a prescribed bubble diameter, on the diffuser flow, exemplary tests with a polydisperse model are conducted. From investigations on bubble columns, e.g., Domgin et al. (1999), Krepper et al. (2007), Cheung et al. (2008), Duan et al. (2011), andXiang et al. (2011), it can be concluded that polydispersity is associated with strongly different bubble velocities, so that small bubbles move towards the wall and large bubbles accumulate within the main flow, which has a substantial impact on bubble break-up and coalescence. Thus, it can be assumed that a polydisperse bubble field may also be associated with the accumulation of coherent gas clouds, e.g., within the blade channel of centrifugal pumps and, therefore, with the head drop. A population balance model based on separate bubble property classes with inhomogeneous velocity treatment of the dispersed phase (iMUSIG = inhomogeneous Multiple Size Group) by Krepper et al. (2007) is utilized, where on the one hand, the size population of bubbles, and on the other hand, the change of bubbles to different bubble size classes is enabled. The population balance model is applied together with bubble interaction models in terms of break-up (Luo and Svendsen, 1996) and coalescence (Prince and Blanch, 1990). A variation of the particular interaction model parameters is performed and presented in the result section 4.1. All computational grids are generated with the commercial software ANSYS ICEM CFD ① .

Diffuser
The computational domain is adopted from the experimental ① https://www.ansys.com setup and is sketched in Fig. 7(a). Since a high grid quality is particularly crucial for two-phase flow, as shown in Müller et al. (2015Müller et al. ( , 2016, hexahedral grids with virtually orthogonal cells are applied. Cell quality and further key simulation parameters are summarized in Table 2. The cross-section geometry is rectangular and, therefore, symmetric around the horizontal as well as vertical axis. Nonetheless, in Kopparthy et al. (2020), the entire cross-section was included in the computational domain, and a significantly asymmetric flow field within the cross-section was obtained even for steady-state simulations despite a perfectly symmetric grid topology and sufficient solver convergence. Thus, to circumvent artificial flow field asymmetries, only one quadrant of the computational domain is considered for single-phase flow simulations in the present study, applying symmetry boundary conditions. For two-phase flow simulations, to allow a vertical gas phase distribution due to buoyancy, a half cross-section is considered by cutting in half the domain and employing a vertical symmetry plane. Note that in Fig. 7(a), the computational domain with the entire cross-section is depicted. The inlet of the computational domain is chosen in the square inflow duct at a distance far upstream of the diffuser so that at the diffuser entrance, a Homogeneous turbulence field, Shear Stress Transport (SST) (Menter, 1994;Menter et al., 2003) with automatic wall treatment (Menter and Esch, 2001;Vieser et al., 2002), k-ω and k-ε based explicit algebraic Reynolds-stress models (BSLEARSM, keEARSM) (Wallin and Johansson, 2000), k−ω-based differential Reynolds-stress model (Launder et al., 1975) (BSLReyStr) with automatic wall treatment (for ω-based) (Menter and Esch, 2001;Vieser et al., 2002), and scalable wall function (for ε-based) (Grotjans and Menter, 1998)   fully-developed flow profile is obtained. A homogeneous distribution of gas bubbles is assumed within the duct inlet cross-section. Albeit in the experiment, the bubble distribution might not be homogeneous immediately after injection, at a certain distance between duct inlet and diffuser entrance, bubbles accumulate at the top of the duct due to gravity, so that in fact there is a fully-developed flow with a gas layer at the duct top wall entering the diffuser entrance.
A monodisperse bubble distribution is chosen, i.e., the bubble diameter d B is constant within each simulation run. A variation of prescribed d B in the range of 1-8 mm is performed. It should be noted that all single-phase simulation results are obtained by a steady-state solution in terms of stationary solver set-up, whereas two-phase simulations need to be performed in a time-accurate way in terms of a timemarching scheme with a Courant-Friedrich-Lewy number around one to achieve a sufficient residuum drop what is in accordance with our experience from former pump singlechannel simulations (Müller et al., 2015(Müller et al., , 2016. Despite the time-accurate run, the flow variables only marginally vary in time after a certain initial transient phase, so that for post-processing, a time-average is inessential, and the last recorded time step may be used. We perform a grid study by a successively refined grid that is summarized in Table 3. For the grid study, water and air flow rate in terms of Q w = 11 m 3 /h and Q a = 60 L/h are chosen that are considered representative for the operation range of the diffuser. Pressure at the bottom wall and GVF at top wall are evaluated at the measurement probe positions that are depicted in Fig. 7(b). Exemplary, the difference pressure p i -p 2 at the bottom wall is illustrated in Fig. 8(a) and shows a remaining grid dependence even on G3, which further reduces with a further grid refinement (not shown here). In Fig. 8(b), the GVF is shown exemplary at the top wall and is essentially grid-independent on grid G3. Thus, for the results presented in Section 4.1, grid G3 is utilized.
It should be pointed out that bubble diameters larger than d B ≥ 4 mm exceed the average cell size in terms of the third root of cell volume on the finest diffuser grid G3. We will discuss in detail in Section 4.2.2 this possible limitation of the disperse two-fluid model. Here, we confine ourselves with the fact that we have verified the validity of the conclusions drawn from our simulations with d B = 4 mm by comparative simulations with smaller bubble sizes what will be demonstrated in the result section 4.1.2.

Test pump
The simulation set-up is adopted from Hundshagen et al.  3D simulation of gas-laden liquid flows in centrifugal pumps and the assessment of two-fluid CFD methods

195
(2019a, 2019b) and is briefly summarized here. Since no rotation-symmetrical impeller discharge is present (see Fig.  9(a)), the computational domain includes the entire pump geometry, i.e., pressure and suction pipe, side chambers, 360°-impeller, and volute casing. It is noteworthy that the ratio of the bubble diameter to the outer impeller diameter d B /D 2 in the simulation is in the order of 3 3 10 -. A pure hexahedral computational grid with a one-to-one grid connection within each reference frame domain and fulfilling grid quality parameters according to Müller et al. (2016), e.g., minimum cell face angle, is used. Three different grid refinement levels are applied: a coarse grid (named G1), a medium grid (named G2) with bisected node distance in all three spatial directions compared to G1, and a fine grid (named G3). G3 is based on the grid G2 and is again refined in the close-to-blade region by bisected node distances. The coarse and fine side of G3 is connected via a general grid interface (GGI) in the simulations. Table 4 lists the main quality parameters of the three grids. The results of the grid study are presented in the result section 4.2.2.
In the impeller geometry, rounded blade trailing edges are ending flush at the circular impeller outlet area. Hence, two concentric circles appear in the blocking of the impeller geometry, which would produce zero-angled computational cells (see detailed view in Fig. 9(a)). Therefore, the trailing edge must be cut artificially by a linear approximation, which leads to a slight geometry simplification (see Figs. 9(b)-9(d)). We assume that the influence of this simplification on the simulation results has a minor effect on gas accumulations that are primarily expected upstream of the trailing edge. Nevertheless, in the result section 4.2.2, we provide a discussion in the single-phase trailing edge flow. In further studies, a verification of the slight geometry simplification, e.g., by unstructured tetrahedral grids, is envisaged.
All two-phase simulations are started from scratch, i.e., with velocity-fields and GVF-field initialized by zero, respectively, to emulate a pump starting process and Core-hours per run (single-/two-phase) 400/6000 6400/-24000/prohibit influences of the initial solution. Time-averaged values of pump head and inner efficiency of the last revolution are evaluated after it has been ensured that the time-average of pump head and impeller momentum deviate by less than 1% over the last three revolutions. This procedure leads to approximately 20-25 revolutions for each simulation run, with an approximate consumption of 6000 core-hours per two-phase flow simulation run on G1. For single-phase flow simulations, due to a better convergence per time step and a requirement of fewer revolutions to obtain statistical convergence, only 400 corehours per run on G1 are consumed, which increases to 6400 and 24,000 core-hours for G2 and G3, respectively. We assume that the combination of a sufficient residual dropresidual root mean square less than 10 − 5 -with this temporal convergence assessment represents a criterion to ensure an utmost quality of the numerical solution.

Bubble size variation
Simulations are performed on grid G3, as discussed in Section 3.2.1, utilizing the SST turbulence model and a monodisperse bubble size distribution. The simulation results are assessed by the existence, size, and shape of air cavities that are observed in the measurements (Mansour et Fig. 9 Explosion-view of the computational domain for the CI with a detailed view of the geometry at blade's trailing edge (a) and grid at the blade trailing edge for three grid refinement levels (b-d).
al., 2018a) and whose time-average is depicted in Fig. 10(a) for a selected operation point Q w = 11 m 3 /h and Q a = 42.42 L/h, entitled OP 1 in what follows. In Fig. 10(b), the simulation result for d B = 4 mm in terms of the gas volume fraction (GVF) in the vertical symmetry plane, i.e., mid-channel, is presented. Only a tiny gas cavity is observed in contrast to the measurement where a quite extensive gas accumulation is present according to Fig. 10(a).
An inspection of other planes than the mid-channel plane shows the same result. It should be pointed out that we performed preliminary simulations with the full computational domain, i.e., no symmetry boundary conditions, which result in the same finding. A variation of the bubble size in the range d B = 1-8 mm yields essentially the same result that no extensive gas accumulation as in the measurement is obtained. Also a variation of the operation point in the range Q w = 8-14 m 3 /h and Q a = 7.07-42.42 L/h-an operation range where each combination of water and air flow rate is associated with an extensive cavity at the upper diffuser wall according to the experiment (Mansour et al., 2018a)-does not significantly change the simulation result, i.e., no or a minimal gas accumulation is found, in contrast to measurement data. Only for an artificially high air flow rate Q a = 120 L/h, termed OP 2 in Fig. 10(c), and thus with a significant deviation from experimental conditions, we find a gas accumulation that is comparable in size and location to the measurement. To conclude, the measured gas accumulation is reproduced by the simulation only with the prescription of a strongly increased gas flow rate.
To further explore the influence of the bubble size, we repeat the simulations for selected operation points with the polydisperse iMUSIG model taking into account either break-up or coalescence. Fifteen bubble classes are applied in the same range as the monodisperse bubble diameter variation, i.e., d B = 1-8 mm. This number of classes is sufficient, which we have verified in preliminary simulations. For the dispersed phase, two sets of momentum equations are solved, corresponding to the iMUSIG ansatz (Krepper et al., 2007): one for bubble classes with a bubble diameter d B ≤ 4 mm, and one for bubble classes with bubble diameter d B > 4 mm. At the diffuser entrance, again, a fully-developed two-phase flow profile in combination with a constant bubble diameter d B = 4 mm is prescribed. The break-up and coalescence model parameters are each varied from their default value 1 towards 0.1 and 10. In Fig. 11, the Sauter mean diameter (d 32 ) at the top diffuser wall is depicted.  For an increase of coalescence parameter, d 32 rises, and for an increased break-up, it decreases as expected. Albeit the bubble interaction parameter range is chosen quite arbitrarily, it spans a wide range, so that we see a significant impact on bubble size. However, we do still not find any hint that a large gas cavity at the diffuser top wall may form. Thus, our finding from the monodisperse simulation is confirmed even with a polydisperse bubble distribution, and we can conclude that a monodisperse bubble distribution alone is not the origin of the misprediction of gas cavity formation. Thus, we perform a variation of the turbulence model, which is presented in the subsequent section.

Turbulence model variation
In Section 4.1.1, the k-ω-based eddy-viscosity SST turbulence model has been utilized, which is based on the assumption of isotropic turbulence. In a square duct flow, however, significant non-isotropic turbulence statistics are present that are associated with a distinct secondary flow (Huser and Biringen, 1993;Zhu et al., 2009). We thus utilize explicit algebraic k-ε-based and k-ω-based (baselineframework) Reynolds-stress models adopted from (Wallin and Johansson, 2000), termed keEARSM and BSLEARSM, respectively, as well as a k-ω-based differential Reynoldsstress model (BSLReyStr), which was adopted from (Launder et al., 1975) and reformulated for the k-ω-(baseline) framework. The specific implementation of the models can be found in (Ansys, 2017).
All models are combined with a flexible wall treatment in terms of automatic wall functions, i.e., the k-ε-based model is coupled with a scalable wall function (Grotjans and Menter, 1998) and the k-ω-based models utilize a blending between logarithmic wall function and a low-Reynolds wall treatment according to Esch and others (Menter and Esch, 2001;Vieser et al., 2002).
In Fig. 12, exemplary results are shown in terms of wall pressure distribution and velocity profile and compared to measurement data. The SST model completely fails to predict the separation that is also reflected by a lower pressure recovery compared to the non-isotropy resolving models. Any model underestimates the peak velocity, and BSLReyStr best reproduces the local backflow at the top and bottom wall. Since with BSLReyStr no solver convergence could be obtained for the two-phase flow, the BSLEARSM is utilized in the following in combination with the monodisperse two-fluid model and d B = 4 mm.
Results in terms of GVF and secondary velocity vectors within the duct cross-section at the diffuser entrance are presented in Fig. 13. According to the BSLEARSM, secondary vortices are present that move gas towards the channel center, which is schematically sketched by large arrows. The accumulation towards the channel center can also be confirmed by visual inspection of the measurement images and is completely missing in the SST simulation results, where a thin gas film at the top wall without any local thickening is predicted. This gas accumulation at the duct top wall center upstream from the diffuser entrance is associated with a completely different cavity formation within the diffuser than predicted by the SST model ( Fig.  10(b)).
In Fig. 14, it can be seen that the size of the gas accumulation is now predicted in good agreement to data, albeit its beginning is slightly more upstream, and its   extension is somewhat overestimated, which might be traced back to the monodisperse assumption. This finding could also be confirmed for further operation points, which are not shown here.
At the end of Section 3.2.1, we have pointed out that the bubble diameter d B = 4 mm exceeds the cell size in terms of the third root of cell volume on the finest diffuser grid G3 what might put the validity of the disperse two-fluid model into question. Thus, we have repeated the simulations with d B = 2 mm, which does not exceed the average grid size limit and obtained essentially the same results that are presented in terms of secondary flow pattern in Fig. 13 and extension of the gas cavity in Fig. 14. On the pump flow in Section 4.2.3.1, we will demonstrate that a bubble size exceeding the grid cell size may cause numerical issues in terms of solver convergence and grid dependence of results. It is interesting to note that we have not observed such numerical issues for the diffuser flow case, even for large values of bubble diameter, maybe since we utilize a diffuser grid with essentially homogeneous grid cell size distribution, and the ratio of d B to the grid size is still moderate.
Summarizing, the resolution of secondary flow in the duct, as well as the prediction of diffuser separation by an adequate turbulence model, are decisive for the prediction of the diffuser two-phase flow field. Thus, the turbulence model and particular consideration of un-isotropy rather than the bubble diameter or polydispersity seems to be the key to reproduce the main flow pattern in terms of gas cavity size. It is interesting to note that Kopparthy et al. (2020) obtained a distinctive gas accumulation even by utilizing the linear k-ε model for the simulation of the same diffuser test case. A VOF method rather than a disperse multi-fluid model has been utilized by Kopparthy et al. (2020) so that it can be assumed that the turbulence model performance is closely associated with the particular multi-phase model it is combined with.

Choice of operation points
In the measurements, a multitude of operation points in terms of water and air flow rate have been investigated and presented by performance maps in Mansour et al. (2018c). For significantly less operation points, detailed flow visualizations that have been presented in Section 2.2.2 are available. The operation points that we have chosen for the simulation study are marked in Fig. 15. The pump head and Fig. 15 Extract from the measured performance maps for CI (a) and OI with tip clearance gap s/b2 = 2.5% (b) (Mansour et al., 2018c). The operation points that are selected for simulation are marked. 3D simulation of gas-laden liquid flows in centrifugal pumps and the assessment of two-fluid CFD methods 199 flow rate are normalized by the best efficiency conditions of the CI in the experiments (H opt and Q opt ), i.e., at CI's nominal flow rate. These two-phase flow operation points are close to the nominal flow rate so that the flow is not dominated by incidence effects and separation that are an additional challenge for turbulence modeling and are avoided here.

Grid study
Firstly, we perform a grid study for single-phase flow that is extended to two-phase flow further below. Despite the better performance of the non-isotropy resolving BSLEASRM for the duct and diffuser flow, for the pump flow, no better prediction is obtained compared to the SST model, neither for single nor for two-phase flow, as will be detailed further below in Section 4.2.3.2. Thus, the SST model is applied in the following.
Simulation results are compared to experimental data in Fig. 16 in terms of pump head for both impeller types and three different grid refinement levels. For the CI (see Fig.  16(a)), a supposedly good agreement to the experimental data in terms of single-phase pump head is observed on G3. The results between the grids G1 and G2 hardly deviate, slightly more pronounced towards overload operation. However, between the grids G2 and G3, the deviation increases, so that the results on the fine grid G3 cannot yet be considered grid-independent. One main reason for the distinctive grid dependence and thus the remaining deviation from measurement data is assumed to be due to the simplified rounded blade trailing edge, which is associated with a geometrically un-defined flow separation location. This separation location is strongly grid dependent, as a detailed inspection of the trailing edge separation (not shown here) figures out. As will be outlined later in this section, a further grid refinement is not feasible for the two-phase flow investigation since the limitations of the disperse multi-phase model prohibits this.
In the OI simulations (see Fig. 16(b)), the measured head curve is only qualitatively captured, since the experimental data is numerically under-predicted. This underestimation is at first glance unexpected since several authors experienced an over-prediction of pump head in the simulation of centrifugal pumps, e.g., Barrio et al. (2010), Limbach and Skoda (2017), and Melzer et al. (2019), which is assumed to be attributed to the neglect of roughness effects in the simulation or the simulation's limitation of resolving secondary flows, e.g., flow separation. However, in other studies also an underestimation of experimental data is observed, which is assumed to occur mainly due to limitations of both the numerical grid and the turbulence model, e.g., Zhang et al. (2018) and Perissinotto et al. (2019). In our case, we assume that besides the rounded trailing edge, particular the relatively coarse and complex tip clearance gap meshing of the OI contributes to the mismatch of experimental and numerical data, which is discussed next in the context of Fig. 17.
For a more detailed tip clearance gap investigation, we consider a single-channel of the impeller and the side gaps and replace the volute casing by a flow contraction for flow acceleration that avoids separation in the impeller outflow region and thus enhances the stability of the simulation. The simulations are performed in the rotating frame of reference, and thus in a steady-state mode. In addition to the tip clearance of s/b2 = 2.5%, a larger clearance of s/b 2 = 5% is investigated. A further grid refinement is performed on the single-channel model by the use of four grids that are summarized in Table 5. The resolution of SC G1 and SC G2 corresponds to the resolution of the full impeller grids G1 and G2. Grid SC G3 is obtained by bisection of the respective coarser grid SC G2, and SC G4 is obtained by a local refinement of SC G3. Note that SC G3 and SC G4 are not feasible for the full pump simulation due to computational effort so that they are restricted to the single-channel-e.g., the grid resolution of SC G4, transferred to the full impeller   Fig. 17, the head for the three impellers and the four grid levels is depicted at CI's nominal flow rate. Note that the evaluated pump head is generally higher with the single-channel model compared to the full pump model for the distinct impeller types since losses due to rotor-stator interaction are neglected. As expected, the head of the CI is higher than for the OIs and attenuates with larger tip clearance gap. Both observations are in accordance with the measurements (not shown here). For the CI, while still no final grid-independent solution is obtained, the head at least flattens towards the finest grid SC G4 and seems to approach a grid-independent value. A closer inspection of the local flow field (not shown here) reveals that the location and extension of separation at the rounded trailing edge are significantly grid-dependent. The different separation pattern is immediately associated with the grid-dependence of pump head, which demonstrates the particular challenge for the simulation in terms of a rounded trailing edge.
For the OIs, the head still significantly rises between SC G3 and SC G4. We assume that this significant griddependence of the OIs is attributed to the spatial gap resolution that is still insufficient, which will be demonstrated next.
Based on the grid SC G2, the number of cells is locally increased in the gap by inserting a separate grid domain in the environment of the gap, which is connected by a GGI to the otherwise unchanged blade channel resolution of SC G2. In Table 6, the gap refinement levels of the gap grid domain, termed G2A, G2B, and G2C, are listed. For G2A, the number of cells in the blade tip direction is quintupled, Fig. 17 Grid study for the single-channel simulation at CI's nominal flow rate. resulting in 20 cells. For G2B, it is increased tenfold, resulting in 40 cells in blade tip direction. G2A and G2B thus correspond to a 1D refinement in blade tip direction. For grid G2C, the grid number within the gap is quintupled in each direction, resulting in 20 cells in blade tip direction and in 125 times as many gap cells as SC G2, corresponding to a local 3D refinement of the gap grid domain. In Table 6, a sketch of the 1D and 3D refinement procedure is included. The results of the local gap grid refinement are also depicted in Fig. 17 in terms of G2A to G2C. The gap grid refinement has essentially the same effect for both gap widths, i.e., the head rises with refinement but does not approach a grid-independent solution. The grid dependence is more pronounced for s/b 2 = 5%, presumably due to the lower spatial resolution than for s/b 2 = 2.5%, at the same number of cells. The significant difference between G2A and G2C shows that it is not only the number of cells in blade tip direction but also in the other directions that have a significant impact on the gap flow and thus on the head. A local inspection of the gap flow (not shown here) reveals that a vortex is generated in the gap close to the leading edge that moves downstream and hits the pressure side of the neighbor blade of the channel so that a strong interaction between the gap flow and the flow field around the neighbor blade exists. This vortex is significantly better resolved by gap grid refinement and associated with the strong gap grid-dependence of results. While the adequate resolution of the tip clearance gap is an important issue even for single-phase flow, it might be of particular importance for two-phase flow where the gap flow is associated with the flushing of gas accumulations and the improvement of twophase pump performances as we reviewed in the introduction section 1.
Summarizing, we have identified two critical zones that are prone to a grid dependence of the numerical flow solution. For the separation at the round trailing edge, on the one hand, the vortex generated at the tip clearance gap of the OI; on the other hand, it is challenging to obtain a grid-independent result. Since this issue is not restricted to two-phase flow but is inherent to any pump flow simulation, we confine ourselves here by pointing out that particular difficulty and proceed to the two-phase flow study in the pump.
Also, for two-phase flows, a grid study is performed together with a bubble diameter variation in the d B range  Hundshagen et al. (2019b), and only the main conclusions are summarized here. It has been found that the deviation between two grid refinement levels increases when either the IGVF or the bubble diameter rises. Besides, the solver convergence gets worse with greater bubble diameter, e.g., the residuals in simulations with a bubble diameter of 2.0 mm are approximately two times higher than in simulations with a bubble diameter of 0.5 mm. The higher grid sensitivity and lower solver stability with larger d B are assumed to be attributed to the particular properties of the disperse two-fluid model. The main assumptions of the Eulerian−Eulerian disperse two-fluid model are the existence of a dilute disperse phase and a bubble size which does not exceed the grid-scale, see, e.g., Marschall (2011). A violation of these constraints is associated with a degeneration of the drag force model, which, in fact, is based on the assumption of dispersed bubble distribution within each computational cell. We assume that the rising grid-dependence with increasing d B is associated with the fact that the bubble diameter becomes taller than the surrounding cell size. Table 7 shows the ratio between the average cell volume of the OI's grid and the volume of three different bubble diameters. We assume that the ratio must be < 1 for a valid disperse model. In the investigated bubble size range, it is clear that a diameter d B = 2.0 mm significantly exceeds the grid-scale. A further grid refinement with bisected node distances is, therefore, for the same reason not reasonable. There is a second disperse two-fluid model constraint that is violated and may be associated with higher griddependence with enlarged IGVF: according to Marschall (2011) andHänsch et al. (2012), the morphology change from dispersed to continuous gas phase occurs in an IGVF range between 0.25 and 0.35 and, therefore, the assumption of the disperse multi-fluid model is violated for IGVF > 0.35. Within gas accumulations that are evidently present in the impeller channels as shown by the measurements in Section 2.2.2, no dispersed gas phase is present any more since the gas accumulates to a continuous gas zone so that the assumption of dispersed bubbly flow is violated.
To conclude the grid study, we have to cope with two presumed limitations of the disperse two-fluid model: firstly, the bubble diameter exceeds the grid size, and secondly, the degeneration of the dispersed phase towards a coherent gas accumulation. It is particularly the first constraint that prohibits a further grid refinement so that the utilization of the relatively coarse grid G1 for further investigations presented in this study is preferred. It should be pointed out that for the diffuser flow discussed in Section 4.1.1, we have run in the same issue with the finest diffuser grid G3. However, by a bubble diameter variation, we have not found a fundamental change of flow pattern in terms of gas cavity accumulation in the diffuser, so that we can conclude that the main conclusion from the horizontal diffuser investigations-i.e., the dominating impact of non-isotropic turbulence modeling-is substantially unaffected by the presumed violation of two-fluid model assumptions. Subsequent studies will focus on the improvement of the two-phase model to combine a dispersed two-fluid with an interphase-resolving (Volume-of-Fluid) method as proposed, e.g., in Hänsch et al. (2012), Wardle and Weller (2013), Hänsch et al. (2014), and Shonibare and Wardle (2015) to enable two-phase simulations even on fine grids. In the present study, coarse grid results in terms of grid G1 are presented.

Impact of monodisperse bubble diameter
In Fig. 18, the pump head is presented for a variation of IGVF for simulation and measurement. In the simulation, the bubble diameter is varied in the range d B = 0.5−2.0 mm. The bubble diameter has got a significant influence on the IGVF level, where the head breaks down. For bubble diameter of d B = 0.5 mm in the simulation, which is a good guess for centrifugal pumps according to Minemura et al. (1985), for both impellers, only a small drop of pump head is observed in the investigated IGVF range up to 5%. Thus, a lower head drop than in the measurements is predicted. Increasing the air bubble diameter to 1.0 mm leads to a head drop for the CI at IGVF = 5%, which does not occur for the OI and is associated with a different gas accumulation size between both impellers in this operating point and will be later discussed in this subsection. When d B is further increased to d B = 2.0 mm, for both impellers, a significant head drop is discernible even for IGVF = 3%. In the measurement, however, a gradual head decrease for the CI is observed. The measured OI head drops rapidly at IGVF = 5% since large gas pockets start to form quickly along the blades, reducing the ability of the pump to convey the mixture (Mansour et al., 2018b). It can be concluded that the head attenuation observed in the simulation is strongly dependent on the prescribed bubble diameter, and the distinct measured characteristics of CI and OI are, at best, only qualitatively captured by the simulation.
These results show that the monodisperse bubble diameter has a significant influence on the pump head. For a more detailed analysis, in Fig. 19, 2D contours of Fig. 18 Normalized pump head H/Hopt for two-phase flow with a variation of IGVF and dB for CI (a) and OI with tip clearance gap s/b2= 2.5% (b).

Fig. 19
Instantaneous GVF-contour plot at CI's and OI's mid-span calculated on G1 with different air bubble diameters compared to time-averaged experimental data in front view. 3D simulation of gas-laden liquid flows in centrifugal pumps and the assessment of two-fluid CFD methods

203
instantaneous GVF at CI's and OI's mid-span together with ensemble-averaged images from measurements are shown. The usage of instantaneous images for the simulation is reasonable since, due to only low temporal fluctuations of the gaseous zones, they hardly deviate from the timeaveraged solution, a fact that has been verified in preliminary tests. In the experiments, the air accumulation size is observed in front view by scattered light measurements. It is visually ensured that the accumulation is spread over the whole blade height. Therefore, we assume that a comparison of mid-span results from the simulation with experimental scattered light results is justified. Regions in the experimental data, where higher bubble densities are observable but are not highlighted in red, represent zones of unsteady dense bubble motions, see, e.g., Fig. 19(a). Only stable, i.e., essentially steady gas cavities, are compared between experimental data and simulation results.
A first general observation is that in the simulation, gas accumulations grow with increasing d B and increasing IGVF (see Figs. 19(a)−19(d)). We start our discussion of simulation results with the CI, i.e., Figs. 19(a)−19(c). For IGVF = 1% (Fig. 19(a)), the gas accumulation is predicted at the pressure side close to the leading edge as in the experiment. Regarding the size of the accumulation region, the simulation result of dB = 0.5 mm seems to fit the experimental data best. For higher IGVF = 3%, the gas accumulation switches from the pressure to the suction side according to the measurement, due to a transition from agglomerated flow regime-i.e., start of bubble coalescence and formation of bigger bubbles near the blades-to a pocket flow regime-i.e., a large gas cavity is attached to the blade-as detailed in Mansour et al. (2018c). This gas cavity attached to the suction side grows according to the experiment when IGVF is further increased from 3% ( Fig.  19(b)) to 5% (Fig. 19(c)). Regarding the simulation of the CI, even a qualitative mismatch of the gas accumulation to experimental data is observed. At IGVF = 3% and IGVF = 5% (Figs. 19(b) and 19(c)), the gas cavity has not switched from the pressure to suction side according to the simulation results, and the simulation predicts a quite pronounced gas accumulation still at the pressure side in contrast to the measured result. Obviously, none of the bubble diameters fit well the experimental data. Albeit there is also a suction side gas accumulation present in the simulation, the erroneous pressure side gas accumulation is even more pronounced, so that there is even a qualitative mismatch to the measurement.
For the OI, the gas accumulation is located at the suction side according to the measurement, which is exemplarily shown in Fig. 19(d) for IGVF = 5%. For the simulation results, the accumulation region is essentially at the correct location, i.e., suction side, and its size fits to measurement data best for a bubble diameter of about dB = 1.0 mm. A more detailed analysis of the OI simulation results reveals that gas zones at the pressure side are flushed away by the secondary tip clearance gap flow of the OI, so that only gas accumulation at the suction side is present, in accordance to the measurement data, see Fig. 19(d).
Summarizing, at best qualitative trends can be captured regarding a successive head drop with increasing IGVF. The choice of the bubble diameter in the monodisperse approach has a decisive effect on the head drop. Different bubble diameters need to be chosen for different operation points and different impeller variants-CI or OI. This observation is in agreement with findings by Rutter et al. (2017), Zhu and Zhang (2017), and Dupoiron (2018) who observed that the best choice of the bubble diameter depends on the operation point of ESP. An important observation is that indeed the bubble diameter has a significant influence on the gas accumulations size but has only a subordinate influence on the location of gas accumulation. Regarding the location of gas accumulation, even a qualitative mismatch to measurement data is observed.
The results reveal two potential improvements for two-phase flow simulations of centrifugal pumps, which will be subject of future studies. Firstly, more precise modeling of bubble diameter as, e.g., supposed by Stel et al. (2020b), seems to be promising also for volute-type centrifugal pumps. Secondly, due to the formation of coherent gas structures and the corresponding flow transition between dispersed and continuous gas phase, the need for a kind of hybrid approach is emphasized. This hybrid approach aims at capturing both, the dispersed dilute bubble clusters by a multi-fluid model on the one hand, and sharp interfaces between liquid and coherent gas accumulations on the other hand, e.g., by a VOF-method. A further potential of flow prediction improvement may be turbulence modeling, which is discussed next.

Impact of turbulence model
In the preceding section, we have seen that the simulation fails to predict the switch of gas accumulation from the pressure to the suction side. This switch is observed in the measurements for the CI and occurs when IGVF is increased from 1% to 3%. Again, this switch is illustrated in Fig. 20(d) and is associated with a transition of agglomerated flow to pocket flow. So far, the SST turbulence model has been employed for pump flow simulations. In Section 4.1.2, we have demonstrated on the diffuser flow, that the choice of an un-isotropic turbulence model and the reproduction of secondary flow structures may have a decisive impact on the prediction of gas accumulation by the simulation. Thus, it may be assumed that the misprediction of gas accumulation location observed in the pump flow, particularly in the CI, may be attributed to the isotropic SST model, so that we utilize the same turbulence models already employed on diffuser flow in terms of BSLEARSM, keEARSM, and BSLReyStr on the single-channel set-up introduced in Section 4.2.2. We consider the CI operation point IGVF = 3% corresponding to Fig. 19(b) since there the mismatch of gas location has been particularly apparent, and choose d B = 0.5 mm for the simulation. The simulations are performed on the grid SC G2, which fulfills the presumed constraints of the dispersed two-fluid model, as discussed in Section 4.2.2, for a bubble diameter of 0.5 mm. First, it should be noted that with the single-channel domain and the SST model (see Fig. 20(a)), essentially the same location of the gas cavity at the pressure side as with the full model (see Fig.  19(b) most left image) is observed, what leads to the conclusion that the single-channel model reflects the main flow features. Regarding the other turbulence models, BSLEARSM, keEARSM, and BSLReyStr (results not shown here), although an inspection of the local flow field in terms of GVF shows a slight change of the cavity size, any further turbulence model predicts the gas accumulation virtually at the same pressure side location as the SST model, in contrast to measurement where the gas cavity is attached at the suction side. This result is confirmed by the impeller head evaluation of the single-channel simulations (not shown here) that does show only a marginal change in dependence on the particular turbulence model. Thus, the observed misprediction of gas accumulation location in the CI cannot immediately be attributed to the isotropic assumption of the SST turbulence model. It should be pointed out that all turbulence models tested so far are statistical (URANS) models that, due to a-priori Reynoldsaveraging of the governing equations, may suppress any unsteadiness of the flow field. We, therefore, utilize the SST-SAS model by Egorov and Menter (2008) that directly resolves turbulent structures down to the available grid limit and, hence, may show characteristics of a large-eddy simulation. GVF results are depicted in Figs. 20(b) and 20(c) in terms of a snapshot and a time-average, respectively. A detailed analysis of the single-phase flow that is not shown here reveals that turbulent vortex structures are, in fact, resolved on the employed grid SC G2, which is also reflected in the two-phase result where the instantaneous image ( Fig.  20(b)) significantly deviates from the time-average (Fig.  20(c)).
Albeit for the SAS result, the gas cavity is still located at the pressure side and, therefore, does not reproduce the experimental result (where it is at the suction side), it is at least located further downstream than the SST solution and therefore closer to the measurement result at the lower IGVF = 1%. Of course, it remains speculative to conclude from that observations a better performance of the SAS model, but in any way, the downstream-shifted gas accumulation location obtained by the SAS model again emphasizes the importance of an appropriate turbulence modeling and maybe even the need for scale-resolving turbulence models.

Conclusions
The performance of a multi-fluid model is assessed for gas/liquid two-phase flow in a centrifugal research pump with closed and semi-open impeller. Due to the excellent measurement database, an assessment of model inaccuracies is performed. In a preliminary test case in terms of a horizontal diffuser, inaccuracies of the turbulence model are revealed to have a significant impact on the flow prediction. Furthermore, the limitations of the validity of the disperse two-fluid ansatz are shown either for fine grids or large bubbles. One main limitation is figured out in terms of the violation of the dilute, disperse phase assumption due to locally high disperse phase loading within coherent gas accumulations. In these circumstances, bubble population models taking into account coalescence effects and large gas accumulations, in their present form, are assumed to be of limited use for conditions-e.g., high forces that are associated with coherent gas accumulations-encountered in centrifugal pump flow simulations. In pump flow, even single-phase flow simulations reveal essential requirements on the spatial resolution of the rounded trailing edge and the tip clearance gap flow. The two-phase pump simulations reveal that, at best, a qualitative agreement of gas accumulations and the head drop towards increasing IGVF can 3D simulation of gas-laden liquid flows in centrifugal pumps and the assessment of two-fluid CFD methods 205 be obtained. Besides more advanced turbulence models in terms of scale-resolving capability, VOF methods should be utilized to capture the phase interface at large accumulated gas cavities, requiring a high spatial resolution. Hybrid approaches, which aim at capturing both, dispersed dilute bubble clusters by a multi-fluid model as well as sharp phase interfaces by a VOF method, have been proposed, e.g., by Wardle and Weller (2013) and Shonibare and Wardle (2015) or Hänsch et al. (2012Hänsch et al. ( , 2014 and will be adopted in our further studies for the simulation of gas/liquid two-phase flow in centrifugal pumps. It has to be ensured that the dispersed multi-fluid part of the hybrid model converges towards the VOF model either with grid refinement or with an increased gas load of computational cells. In this hybrid framework, it can be expected that the benefits of population and bubble interaction models are revealed in flow regions where a dilute disperse gas phase is present.