Three-Dimensional Flow Simulation by a Hybrid Two-Phase Solver for the Assessment of Liquid/Gas Transport in a Volute-Type Centrifugal Pump with Twisted Blades

: A hybrid two-phase ﬂow solver is proposed, based on an Euler–Euler two-ﬂuid model with continuous blending of a Volume-of-Fluid method when phase interfaces of coherent gas pockets are to be resolved. In a preceding study on a two-dimensional bladed research pump with reduced rotational speed, the transition from bubbly ﬂow to coherent steady gas pockets observed in optical experiments with liquid/gas ﬂow could be well captured by the hybrid solver. In the present study, the experiments and solver validation are extended to an industrial-scale centrifugal pump with twisted three-dimensional blades and elevated design rotational speed. The solver is combined with a population balance model, and a scale-adaptive turbulence model is employed. Compared to the two-dimensional bladed pump, the transition from agglomerated bubbles ﬂow to attached gas pockets is shifted to larger gas loading, which is well captured by the simulation. The pump head drop with increasing gas load is also reproduced, showing the hybrid solver’s validity for realistic pump operation conditions.


Introduction
Radial volute-type centrifugal pumps are frequently required in industrial processes to transport multi-phase flows consisting of liquid and non-condensable gas.However, they are usually designed to convey pure liquids and, hence, even a small load of gas, i.e., an inlet gas volume fraction ε ≈ 3% [1], may cause a considerable disturbance in conveying.In the worst case, even a full breakdown of the pump head may occur [2], and flow instabilities may destroy the pump's components.
The resolution of experimental measurements is inherently limited [3,4].As pointed out by, for example, Yin et al. [5] and De Santis et al. [6], computational fluid dynamics (CFD) simulation of liquid/gas flow provides a high amount of local flow information and may complement experimental results.In a preceding study [7], we presented a hybrid twophase (H2P) solver for the 3D simulation of a water/air pump flow.A transparent research pump with a 2D cylindrical blade design facilitated detailed optical accessibility and a unique validation database.However, the impeller speed was reduced to n = 650 min −1 to protect the damageable transparent impeller material.It was shown that the transition of flow morphology from dispersed bubbly to segregated pocket flow was captured by the H2P solver.In the current study, the experiments and H2P solver validation are extended to an industrial-scale volute-type centrifugal pump with twisted 3D blades and a design rotational speed n = 1450 min −1 , which is considerably higher than the corresponding design speed of the 2D bladed pump used by Hundshagen et al. [7].
Experimental validation data is essential for assessing CFD models, and optical data of the two-phase flow field offer a sound validation basis.For example, Stel et al. [8,9], He et al. [10], Zhang et al. [11], Parikh et al. [12] and, also, the present authors [7] used optical measurement data to access the prediction of gas distribution within 2D bladed impellers in CFD studies.Recent high-speed imaging experiments [1,[13][14][15][16][17][18] confirmed the existence of gas pockets inside the pump impeller channels, which block a part of the blade channel, making the flow guidance inefficient and causing the head drop.Based on optical measurements, together with the head drop characteristics, several authors [1,8,[13][14][15][19][20][21][22][23][24] provided performance maps of the conveyance of a mixture.Subject to the volumetric flow rate and the inlet gas volume fraction, some of the present authors [1,[13][14][15] observed a transition from bubbly flow to steady gas accumulation in the blade channel and introduced a categorization in terms of flow regime maps.Experimental studies by Mansour et al. [1,[13][14][15], and the corresponding simulations with the H2P solver by Hundshagen et al. [7], were, however, confined to a 2D bladed impeller at reduced speed, as mentioned above.In studies where 3D twisted blades were used, the assessment of CFD models was restricted to integral pump characteristics, such as those of the pump head, or inner efficiency [25,26].In this study, optical measurements of a 3D bladed impeller, having an elevated design speed, were performed and used for the validation of the H2P solver, which is a major novelty of this study.
The optical experimental results provided by Mansour et al. [15] and Parikh et al. [12] also pointed out that there is a significant size variation of gas bubbles within the pump.Chen et al. [27] even concluded that coalescence is one of the primary mechanisms for the formation of gas pockets.Inspired by the observations made in experiments, the H2P solver is extended by a population balance model (PBM) to take into account bubble polydispersity, which is a further novelty of this study.The class method of Lo [28] is adopted, together with coalescence and breakup models.
A high unsteadiness of particular void regions was observed in experiments by Mansour et al. [15] and Liao et al. [29], which is also confirmed by an inspection of video sequences in this study.Alternating unsteady pockets are present in the transition zone between bubbly and pocket flow regime.Even inherently steady adherent void regions in the pocket flow regime show an unsteady wake [7].Statistical turbulence models may be inappropriate to resolve the highly unsteady flow regions observed in experiments, so that a turbulence-scale resolving approach is preferred.The scale-adaptive simulation (SAS) approach of Egorov & Menter [30], together with the hybrid H2P solver, has already been applied in our preceding study [7].The SAS treatment of turbulence is also adopted in this study.The liquid/gas flow is investigated in terms of a water/air system in the following.This manuscript corresponds to our paper published in the proceedings of the 15th European Turbomachinery Conference [31].

Experimental Set-Up and Pump Characteristics
A 3D view of the pump and a photograph of the test rig are provided in Figure 1.The pump casing and the impeller are made of transparent acrylic glass and transparent polycarbonate, respectively, to allow flow visualization.The two-phase flow patterns are captured using two high-speed cameras with a 2016 × 2016 pixel resolution.The cameras are located in axial (camera 1) and isometric (camera 2) views, as depicted in Figure 1b.
Acrylic glass, polycarbonate and water have different optical refractive indices.In this study, only images obtained by camera 1 are used for direct comparison with simulation results, while images from camera 2 serve as a means for additional visual inspection of air accumulation zones.The arrangement of camera 1 is orthogonal to the casing wall.Therefore, it is assumed that the effect of refraction is minor.For a more quantitative comparison, the actual location and size of bubbles may be calibrated and double-checked, based on the geometrical pump dimensions, such as impeller diameter and blade thickness.However, this is out of scope of this study.Six LED lamps are distributed around the pump casing to illuminate the flow, as shown in Figure 1a, making the water/air interfaces visible.A dashed rectangle in Figure 1a indicates the window where the two-phase flow patterns are observed.A rotational speed of n = 1450 min −1 is applied, resulting in a specific speed of about n q = 24 min −1 , according to the definition of Gülich [32], and corresponding to a non-dimensional specific speed of N s = 0.44.A semi-open impeller is applied, and its meridional view is sketched in  It is found that the void patterns are not significantly affected by tip clearance size.Therefore, this study focuses on the smaller clearance gap.The volumetric flow rates of air and water, Q a and Q w , are separately measured to calculate the inlet air volume fraction ε.The static pressure difference across the pump is recorded to calculate the pump head.No bubble size spectra are measured.Therefore, we confine our investigations to a qualitative comparison between experimental and simulation results.Details of the experimental set-up, measurement techniques, and data evaluation are provided elsewhere [1,7,15,33].

Flow Regimes
The flow regime map observed in the experiments, in terms of pump head H versus total flow rate Q t is depicted in Figure 3.The bubbly flow, agglomerated bubbles flow, alternating pocket flow, and pocket flow regime are discernible.The basic structure of the regimes essentially corresponds to the one of the 2D bladed pump from our previous study [7].However, there are distinctive shifts in individual regimes.A detailed analysis of the effect of blade geometry and impeller design speed on the flow regimes is out of scope of this study.We confine our discussion to the effects which are immediately helpful for solver validation.For example, the bubbly flow regime is shifted to overload conditions, while the agglomerated bubbles flow regime is near the nominal flow rate, presumably because a rise in speed enhances Coriolis forces, so that the pressure gradient in the cross-flow direction and the tendency of bubble agglomeration are also enhanced.Furthermore, the pocket flow regime is present at rather large values of the inlet air volume fraction, ε ≥ 9%.A wide surging range between 3% < ε < 9% is discernible, which means that the pump jumps back and forth between two operation points.Alternating pockets occur only at part-load and overload, and at moderate air loading 3% < ε < 7%.Ref. [15] showed that, by increasing the rotational speed, the water/air mixing and, therefore, the two-phase pump performance is improved.Hence, it is assumed that the shift of pocket flow to higher levels of ε can be traced back to a rise in design speed.It is interesting to note that the Reynolds number evaluated by the impeller diameter and circumferential impeller outlet velocity is the same for the present pump and the 2D bladed pump of our previous investigations [7], so that the shift of flow regime maps to higher values of ε was not a Reynolds number effect.A more in-depth analysis of experimental data for assessing flow regime maps is left to further studies.In this study, the validation of the H2P method for a 3D bladed pump with elevated design speed is the focus.Therefore, three operation points are considered with ε = 1%, 3%, and 9%, additionally marked in Figure 3. ε = 1% and 3% correspond to the agglomerated bubbles flow regime, while ε = 9% corresponds to the pocket flow regime.

Flow Solver
The H2P solver from our preceding studies [7,34] is adopted, and its main features are summarised here.It is an extension of the Euler-Euler two-fluid (EE2F) model with the ability to continuously activate Volume-of-Fluid (VOF) features for a proper interface resolution of coherent void structures.According to Wardle & Weller [35], the balance equations for mass and momentum for each isothermal phase ϕ (water and air) read: 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 the marking of Reynolds-averaged quantities is omitted for convenience.Both phases share the same static pressure field p. Water and air are treated as incompressible fluids, which enhances solver stability [36].The interfacial momentum exchange terms, i.e., drag and virtual mass force, are symbolized by M D,ϕ and M vm,ϕ .For the former, the drag model of Schiller & Naumann [37], and, for the latter, a constant virtual mass coefficient of 0.5 [38], are applied, respectively.The force M s,ϕ is the surface tension force, estimated as a continuum force, according to Brackbill et al. [39].
The hybrid multi-scale character of the solver is facilitated by two blending functions for the interface compression term in Equation ( 1) and the drag force M D,ϕ , respectively.By following Hänsch et al. [40], the blending function for the interface compression C α reads: We set the constants a B , a min , and a max to 35, 0.1, and 0.9, respectively.If C α is set to 0, the original EE2F formulation is retained, while with C α = 1, a VOF-like interface resolution is achieved, which counteracts the numerical diffusion at the phase interface by including the compression velocity c c in Equation (1).Details of the c c evaluation can be found in Wardle & Weller [35].A second blending function C β blends the drag force between the formulation of [37] for spherical bubbles and a treatment of large coherent air structures where the drag is set to a low finite value of 10 −4 .The C β function is presented in [7] and is not repeated here.Note that applying two different blending functions for the interface compression and the drag force is attributed to the complexity of the pump flow, since it would be more consistent to use the same blending function for both sub-models.However, using Equation (3) for the drag force as well has produced stability issues [7].
The limitations of a monodisperse approach have been pointed out by Hundshagen et al. [7].Therefore, a polydisperse bubble size distribution is considered here.An additional transport equation for the specific number density n B (d B ) of bubbles is solved in each computational cell and timestep.Beyond the transport of n B (d B ), birth and death rates, due to coalescence and breakup, are treated by source and sink terms.Thereby, a variable bubble diameter d B is facilitated, e.g., in the drag force M D,ϕ .By adopting the approach of Kumar & Ramkrishna [41], the population balance equation is discretized in size groups, leading to the class method of Lo [28].Twelve size classes, with a duplication of the bubble volume between consecutive classes, are adopted, resulting in a spectrum of bubble diameters for each size class, ranging from roughly 0.16 mm < d B < 2.00 mm.This bubble diameter range corresponds to the one observed by Barrios & Prado [20] in electrical submersible pump flow for speed variation.Each bubble size class comprises the same velocity, and coalescence and breakup kernels are evaluated for each class.The coalescence kernel of Prince & Blanch [42], and the breakup kernel of Luo & Svendsen [43] are chosen.The coalescence kernel is used with an extension for maximum bubble packing, as proposed by [44], resulting in a formulation similar to the one suggested by Stel et al. [9] for a centrifugal rotor flow.Although the breakup and coalescence kernels are simply adopted from entirely different applications, e.g., bubble columns or mixer vessels, they have also been adopted in other two-phase pump studies, such as those by He et al. [10] and Yan et al. [25].Hence, this is also done in this study.From the resulting bubble size distribution, the Sauter mean diameter d 32 is evaluated, corresponding to a time and space distribution of d B .
The SAS approach by Egorov & Menter [30] is adopted, which is based on the k-ω shear stress transport turbulence model of Menter [45].A separate turbulence field, in terms of turbulence kinetic energy k ϕ and specific dissipation rate ω ϕ , is calculated for air and water.The wall functions of Kalitzin et al. [46] are employed for both phases.Hence, c ϕ , k ϕ , and ω ϕ are directly integrated to the wall or calculated by the logarithmic wall law depending on the value of the non-dimensional wall distance y + w,ϕ .The multiphaseEulerFoam solver of OpenFOAM-9 (https://github.com/OpenFOAM/OpenFOAM-9 accessed on 1 May 2022) was used, customized by means of several model extensions, as described above.The PIMPLE solution algorithm with pressure relaxation was chosen, which is a combination of the pressure-implicit with splitting of operators (PISO) [47] method and the semi-implicit method for pressure-linked equation (SIMPLE) [48] algorithms.The timestep is set in such a way that the Courant-Friedrich-Levy number does not exceed 0.75 at the maximum, to achieve numerical stability.Three outer PIMPLE iteration loops turned out to yield a drop of non-linear residuals below 5 • 10 −4 in each time step.Additionally, a drop of residuals of linear solvers below 10 −8 is applied.Different discretization schemes for the respective balance equations are used, which are summarized in Table 1.
Table 1.Discretization schemes of the simulations.

Simulation Set-Up
The simulation set-up consists of the suction and discharge pipe, volute casing, side chamber, and 360 • impeller to reflect the entire pump geometry.An explosion view of the computational domain is shown in Figure 4. Constant Dirichlet inlet boundary conditions are chosen for c ϕ , α ϕ , d B , k ϕ , and ω ϕ .The phase velocity is set according to the measured flow rate and water and air have the same inlet velocity.A constant bubble diameter of 0.5 mm is set at the inlet, following the estimation by Barrios & Prado [20].At the suction side cross-section, i.e., the inlet boundary of the computational domain, a constant value of α ϕ is prescribed.A Neumann boundary condition with zero gradient is set for pressure at the inlet, while a constant static pressure is defined at the outlet.At the outlet, Neumann zero gradient boundary conditions are set for c ϕ , α ϕ , d B , k ϕ , and ω ϕ .
Numeca AutoGrid® is used for the impeller and side chamber grid generation.For the suction and discharge pipes and the volute casing, Numeca HexPress® is employed, and local grid refinement by hanging nodes is utilized.The grids are optimized in terms of non-orthogonality, following the OpenFOAM nomenclature of Jasak [47].In this way, grids with reasonably high quality are generated, taking into account the complexity of the 3D-twisted blade design.A grid study, by means of refining grid level G1 to G2, is performed, and some grid parameters are summarized in Table 2.The single-phase simulation results, regarding the pump head and the inner efficiency for the same impeller geometry, but with a different volute casing, were investigated elsewhere [53,54].Therefore, the single-phase results are not again discussed here in detail.Preliminary simulation results showed that the change in single-phase pump head between G1 and G2 simulations approaches at maximum 1.5% in over-load conditions and less than 0.5% in nominal and part-load conditions.Regarding the two-phase flow, by means of a grid refinement, successively finer void structures are resolved by the H2P solver, which transitions from EE2F mode towards a VOF-like scheme with phase-interface resolving capabilities.Thus, an inherent grid dependence of the results cannot be avoided.The same holds for the SAS, since a successively larger portion of the turbulence spectrum is directly resolved by the computational grid.Details of a two-phase grid study are omitted here for brevity, but were presented elsewhere [7], and it is ensured that the conclusions drawn from the simulation results in the following are not affected by grid dependence.The G1 results are described below.

Results and Discussion
In the 2D bladed pump of our prior study [7], air accumulations extended over the entire blade height, which enabled a direct comparison between experimental results in the front view and contour plots of air volume fraction in the midspan of the impeller.This assumption does not hold any more for 3D-twisted blades.Figure 5 presents plots over lines located near the hub, midspan, and shroud in the middle of a blade-to-blade surface in an example blade channel.The aligned coordinate along the line is denoted by l, and l max is the line length.It can be seen that significantly less air is located near the shroud region than near the midspan or hub regions.Therefore, the α a distribution changes over the blade height, and a direct comparison of the air distribution from simulation results with experiments is facilitated by virtual Schlieren imaging.This visualization technique of the simulation results is based on spatial gradients of α a , and described elsewhere by Nguyen et al. [34].Virtual Schlieren images are shown in Figure 6, together with the grey-scale images obtained by the experimental scattered-light technique.White areas in the virtual Schlieren images and the experimental results represent regions with high dispersed air loading.Dark regions indicate continuous phase regions of either water or air.Example snapshots of the flow field are shown, and coloured arrows are depicted to highlight similarities between simulation and measurement results derived from an intensive inspection of video sequences.For ε = 1% (Figure 6a), a maximum of dispersed air loading is visible near the suction side of the blade (red arrow), while a minimum is discernible on the pressure side (blue arrow).For ε = 3% (Figure 6b), regions of dispersed air loading increase in size.Thus, the tendency of water and air to separate increases with a rising value of ε.For ε = 9% (Figure 6c), a vast region with continuous air phase is discernible with maxima of dispersed loading downstream of this continuous region.Hence, the local air distribution obtained by the virtual Schlieren results from the simulations agrees well with the experimental results for the investigated operation range from ε = 1% to ε = 9%.
Contour plots of the time-averaged air volume fraction ᾱa and its temporal standard deviation α a,RMS are shown in the impeller midspan in Figure 7a and b, respectively.Even for ε = 1%, a slightly elevated level of ᾱa is discernible near the suction side of each blade, which means that a slight tendency for phase separation is present.This tendency has already been observed in Figure 6, and is confirmed here.The highest values of α a,RMS correlate with the highest values of ᾱa .These temporal and spatial characteristics of the air distribution indicate the presence of the agglomerated bubbles flow regime, according to Figure 3.For ε = 3%, the peaks of both ᾱa and α a,RMS increase significantly and span over a greater region within each blade channel.For ε = 9%, the α a,RMS field shows minima in the region of the highest values of ᾱa , which means that the correlation between ᾱa and α a,RMS observed for low levels of ε completely reverses.Downstream of these stationary regions, highly unsteady wakes develop.Hence, a steady air pocket with an unsteady wake is observable.These characteristics can be assigned to the pocket flow regime, according to Figure 3.The analysis of transition between flow regimes has been extensively discussed by Hundshagen et al. [7] for a 2D bladed pump with reduced speed.From our brief discussion here, it is clear that, essentially, the same transition processes are present for the 3D bladed pump with elevated design speed, albeit on higher levels of ε.This transition is well captured by the H2P simulation approach.
The instantaneous distribution of the Sauter mean diameter d 32 at impeller midspan is presented in Figure 8. d 32 rises with rising ε, which indicates enhanced coalescence activity, according to the findings of Chen et al. [27].A cross-check to Figure 7a reveals that this is most pronounced for ε = 9%.The location of large bubble diameter correlate with the air accumulation zones, i.e., preferably, large bubbles accumulate, which is in agreement with similar findings by Zhang et al. [11].Small bubbles are present downstream of the accumulation, indicating that the breakup mostly occurs in the wake of the accumulation zone.Due to these spatial characteristics, the d 32 distribution appears to be reasonable.The time-averaged pump head versus inlet air volume fraction is shown in Figure 9.The experimental trend of decreasing head with increasing ε and the head drop between ε = 3% and ε = 9% are well reproduced.The transition of flow regimes, which is captured by the new H2P approach, is also reflected in a good reproduction of the measured pump head characteristics.

Conclusions
A hybrid two-phase (H2P) solver for the simulation of liquid/gas transport has been presented and assessed by experimental validation data of an industrial-scale centrifugal pump flow with 3D-twisted blades and increased impeller design speed.As a continuation of our previous study [7], it is the elevated speed that raises the pump resistance against air loading [15], probably due to its more substantial mixing effect, which shifts the pocket flow regime to higher values of inlet air volume fraction ε.It could be shown that the transition from bubbly flow to adherent pockets is also reproduced by the H2P method for a more industry-oriented pump design, i.e., twisted 3D blading and elevated design speed, which is a major novelty of this study.
A further novelty of the simulation method is the inclusion of a PBM to capture a variation of bubble sizes.It is interesting to note that the transition between flow regimes can also be captured by a monodisperse approach (not shown here), i.e., constant diameter, which, however, has to be prescribed and is, thus, not a result of the simulation.By the polydisperse approach employed here, a reasonable bubble size distribution is obtained.The hybrid approach, in combination with a scale-adaptive turbulence treatment and a PBM, enables the description of fluid dynamical processes of gas accumulation on a sound physical basis.This is a basic requirement for more predictive CFD tools and a knowledge-based design and optimization of centrifugal pumps for gas-laden liquid transport.Coalescence and breakup kernels have been, so far, adopted from entirely different applications, such as bubble columns or mixer vessels.For a sound validation of bubble spectra and the development of enhanced kernels, more in-depth experiments are indispensable and should be the focus of future studies.

Figure 2 .
Figure 2. Sketch of the meridional view of the semi-open impeller.

Figure 3 .
Figure 3. Pump head performance map and flow regimes observed in the experiments.The operation points investigated in the simulations are marked with an "X".

Figure 4 .
Figure 4. Explosion view of the computational domain.

Figure 5 .
Figure 5. Plots of time-averaged air volume fraction ᾱa along arbitrarily chosen lines in the middle of a blade-to-blade surface near the hub, midspan, and shroud for the ε = 3% operation point.

Figure 6 .
Figure 6.Comparison between virtually-derived Schlieren images (left) and experimental results of a grey-scale analysis from the front view (middle) and isometric view (right) for ε = 1% (a), ε = 3% (b), and ε = 9% (c).Similarities between simulation results and measurement data are marked by red and blue arrows.

Author Contributions:
Conceptualization, M.H. and R.S.; methodology, M.H., K.R. and R.S.; software, M.H. and K.R.; investigation, M.H., M.M., D.T. and R.S.; writing-original draft preparation, M.H. and R.S.; writing-review and editing, M.H. and R.S.; supervision, D.T. and R.S.; funding acquisition, D.T. and R.S.All authors have read and agreed to the published version of the manuscript.Funding: The research project was carried out in the framework of the industrial collective research program (IGF no.20638 BG).It was supported by the Federal Ministry for Economic Affairs and Climate Action (BMWK) through the AiF (German Federation of Industrial Research Associations eV), based on a decision taken by the German Bundestag.Simulations were performed with computing resources granted by RWTH Aachen University under project bund0013.