Benchmarking of computational fluid dynamic models for bubbly flows Nuclear

Eulerian-Eulerian computational fluid dynamic (CFD) models allow the prediction of complex and large-scale industrial multiphase gas – liquid bubbly flows with a relatively limited computational load. However, the interfacial transfer processes are entirely modelled, with closure relations that often dictate the accuracy of the entire model. Numerous sets of closures have been developed, often optimized over few experimental data sets and achieving remarkable accuracy that, however, becomes difficult to replicate outside of the range of the selected data. This makes a reliable comparison of available model capabilities difficult and obstructs their further development. In this paper, the CFD models developed at the University of Leeds and the Helmholtz- Zentrum Dresden-Rossendorf are benchmarked against a large database of bubbly flows in vertical pipes. The research groups adopt a similar modelling strategy, aimed at identifying a single universal set of widely appli- cable closures. The main focus of the paper is interfacial momentum transfer, which essentially governs the void fraction distribution in the flow, and turbulence modelling closures. To focus on these aspects, the validation database is limited to experiments with a monodispersed bubble diameter distribution. Overall, the models prove to be reliable and robust and can be applied with confidence over the range of parameters tested. Areas are identified where further development is needed, such as the modelling of bubble-induced turbulence and the near-wall region, as well as the best features of both models to be combined in a future harmonized model. A benchmark is also established and is available for the testing of other models. Similar exercises are encouraged to support the confident application of multiphase CFD models, together with the definition of a set of experiments accepted community-wide for model benchmarking.


Introduction
Bubbly flows are widespread in a multitude of multiphase flow engineering applications and industrial fields, such as nuclear thermal hydraulics and chemical and process engineering equipment. In bubbly boiling flows, extremely high heat transfer coefficients are reached and the dispersion of small bubbles in a background liquid phase is often employed when high rates of heat and mass transfer are needed between two or more fluids (Risso, 2018). On the other hand, the interaction between the dispersed bubbles and the continuous liquid at the interface between them makes the hydrodynamics of bubbly flows complex and challenging. Bubbly flows are normally polydispersed, with multiple bubbles of sometimes largely different sizes that breakup by interacting with the surrounding liquid and, when the bubble concentration increases, frequently collide and coalesce with their neighbours (Lucas et al., 2010). As a consequence, the density of the interfacial area, the major driver of the desired heat and mass transfer processes, is continuously altered. Empirical correlations and simplified onedimensional models are not equipped to capture such physics that occurs at the bubble scale, as they can only correlate with values of averaged or bulk parameters (Woldesemayat and Ghajar, 2007;Vasavada et al., 2009). Therefore, they have limited accuracy and normally struggle when applied outside the specific range of parameters for which they were developed. In view of this, research has more recently focused on three-dimensional, time-dependent computational fluid dynamic (CFD) methods (Yao and Morel, 2004;Yeoh and Tu, 2006;Hosokawa and Tomiyama, 2009;Dabiri and Tryggvason, 2015;Colombo and Fairweather, 2016b;Santarelli and Fröhlich, 2016;Mimouni et al., 2017;Feng and Bolotnov, 2018;Liao et al., 2018;Lubchenko et al., 2018). These are best equipped to account for the many local phenomena at the bubble scale that impact the macroscopic behaviour of the flow, and provide reliable numerical tools that are much needed to underpin improved bubbly flow understanding as well as efficient industrial equipment design and process optimization.
For bubbly flows of practical relevance, where large-scale, complex geometries with hundreds of thousands or millions of bubbles are involved, multifluid Eulerian-Eulerian models are the preferred choice (Yao and Morel, 2004;Yeoh and Tu, 2006;Hosokawa and Tomiyama, 2009;Liao et al., 2015Liao et al., , 2018Rzehak et al., , 2017Fairweather, 2016a, 2016b;Mimouni et al., 2017;Lubchenko et al., 2018). Although the continuous increase in computational resources has favoured the development of interface resolving approaches where all the interfacial and flow scales are resolved, these remain mainly constrained to much fewer bubbles in simplified flow conditions (Dabiri and Tryggvason, 2015;Santarelli and Frohlich, 2015;Santarelli and Fröhlich, 2016;Bolotnov, 2017, 2018). In contrast, in multifluid models, physical processes at the bubble scale are not resolved and only the spatial and temporal distribution of the averaged (over small volumes) phase concentration is known. A set of conservation equations is solved for each phase and coupling between the phases is achieved with closure relations that model the unresolved exchanges of mass, momentum and energy at the interface (Ishii and Hibiki, 2006;Prosperetti and Tryggvason, 2007;Yeoh and Tu, 2010). It is not surprising that most of the research undertaken using these approaches has focused on the development of more accurate and physically based versions of these closures, given their impact on the accuracy of the overall model.
Of considerable importance are the interfacial forces used to model the momentum interfacial exchange and how continuous liquid and bubble velocities and spatial distributions mutually interact (Colombo and Fairweather, 2015;Rzehak et al., 2017;Liao et al., 2018;Lubchenko et al., 2018). This interaction is modelled with a number of forces that reproduce different physical effects. The drag force models the opposition of the surrounding liquid to bubble motion by interfacial shear. A well-known effect in closed ducts, modelled with the lift force, is the force on the bubble in the direction perpendicular to a wall, and the main fluid motion, induced by the gradient in the same direction in the fluid velocity (Auton, 1987;Tomiyama et al., 2002b). Lift forces considerably alter the bubble spatial distribution, by pushing small, relatively spherical bubbles towards regions of higher relative velocity, e.g. to the wall in upward vertical flows. For larger, deformed bubbles, driven by the altered fluid circulation around the bubble surface, the lift force acts in the opposite direction (Lucas and Tomiyama, 2011). In air-water bubbly flows, experimental evidence suggests that at atmospheric conditions this change in direction occurs for bubble diameters between 5 and 6 mm (Tomiyama et al., 2002b) and, with small bubbles, vertical pipe flows exhibit a peculiar wall-peaked void fraction distribution (Liu andBankoff, 1993a, 1993b;Lucas et al., 2005;Hosokawa and Tomiyama, 2009). In multifluid models, this peak has generally been predicted with a linear superposition of the lift force and an additional repulsive wall force that, at a sufficiently small distance from the wall, prevents bubbles moving closer to it (Antal et al., 1991;Hosokawa et al., 2002;Rzehak et al., 2012). Over the years, numerous lift and wall force formulations have been developed, often optimized over a limited amount of data (Hibiki and Ishii, 2007). Although good predictive accuracy is often achieved over the range of the data selected, extension to other conditions has proven difficult. For example, agreement with data has been reported for values of the lift coefficient ranging from 0.01 (Wang et al., 1987;Yeoh and Tu, 2006) to 0.5 (Mimouni et al., 2010). This, and the multiple combinations of different closure models, clearly impacts the general applicability of multifluid models and complicates any genuine assessment of their overall accuracy (Lucas et al., 2016;Podowski, 2018). At the same time, it obstructs the community from reaching agreement over the best model available and the most pressing developments needed.
An additional open and active area of research is the development of modelling closures for multiphase turbulence in Reynolds-averaged Navier-Stokes CFD approaches. In this regard, modelling still often relies on the eddy viscosity assumption (Yao and Morel, 2004;Rzehak et al., 2017;Sugrue et al., 2017;Liao et al., 2018). More recently, however, driven by the desire to move beyond the limitations of Boussinesq's assumption, well-documented for single-phase flows (Benhamadouche, 2018), progress has been made in the development and application of second-moment Reynolds stress closures (Lopez de Bertodano et al., 1990;Fairweather, 2015, 2020;Mimouni et al., 2017;Parekh and Rzehak, 2018). Recent studies suggest that such closures can account for additional influences of turbulence, and its modelling, on the dispersed phase distribution (Ullrich et al., 2014;Santarelli and Frohlich, 2015;Colombo and Fairweather, 2019). Specific to bubbly flows, and the subject of numerous studies, is the modelling of the bubble-induced contribution to the continuous phase turbulence. This is often modelled based on the conversion of energy from drag to turbulence kinetic energy in bubble wakes (Troshko and Hassan, 2001;Rzehak and Krepper, 2013;Colombo and Fairweather, 2015). Continuous advances are being achieved in this area, including specific implementations for Reynolds stress closures (Ma et al., 2017(Ma et al., , 2020Parekh and Rzehak, 2018;Magolan and Baglietto, 2019;Colombo and Fairweather, 2020).
In response to the multiplication of modelling closures, researchers at the Helmholtz -Zentrum Dresden -Rossendorf (HZDR) have proposed their baseline closure strategy Liao et al., 2018;Lucas et al., 2020). A default modelling setup is established, including closure relations for all the relevant physical processes in bubbly flows. The model is then systematically validated over a continuously increasing experimental database, and a specific closure is modified only if it improves the overall prediction of the entire database, to the benefit of the robustness of the model, even if at the expense of a slight decrease in accuracy over specific experiments (Lucas et al., 2016). This strategy, including validation over large databases of experimental measurements, has been embraced at the University of Leeds (UoL), where the same principles were adopted in CFD multifluid model development Fairweather, 2015, 2020).
In this paper, the models from HZDR and UoL are systematically benchmarked blindly against each other over a large range of adiabatic bubbly pipe flow experiments. The paper is specifically focused on the interfacial momentum closure framework, and the multiphase and bubbleinduced turbulence modelling employed. If present, a population balance model influences the model's overall accuracy and the assessment of other closures, the accuracy of which will also depend on the accuracy of the average bubble diameter prediction from the population balance. For this reason, experiments where a monodispersed bubble diameter distribution was measured are selected for comparison purposes. In these flows, the bubble population can be effectively approximated with a fixed average bubble diameter taken from the experimental measurements, without the need for a population balance. The work aims at systematically identifying the level of confidence and overall accuracy that can be expected when applying these models over the range of flows tested. Equivalent exercises, extending over a similarly large database, are difficult to find in the literature, and the present work also establishes a benchmark that is available for other models to be tested against. The strengths and weaknesses of each model are identified, with the aim of establishing a path towards a future harmonized best possible model as well as pointing out areas where further joint developments will be beneficial.

Computational fluid dynamics model
The CFD models are based on the multifluid Eulerian-Eulerian method (Prosperetti and Tryggvason, 2007;Yeoh and Tu, 2010). The flows considered are adiabatic, and a set of averaged continuity and momentum equations is solved for each phase: In the above equations, p is the pressure, common to both phases, and U k , α k and ρ k are the velocity, volume fraction, and density of phase k, respectively. In the following, we will use the indices c and d to denote continuous liquid and dispersed gas, but we will often refer only to α to identify the void fraction of the gas phase. Indices i and j denote Cartesian coordinates. g is the gravitational acceleration and τ k and τ Re k the laminar and turbulent stress tensors, respectively. The term M k is the interfacial momentum transfer source and models the interfacial momentum transfer between the phases with a set of closure relations that account for the different forces that act on the bubbles. The closures employed are commonly used in the modelling of bubbly flows, although in a multitude of combinations and often with modified coefficients, and both the HZDR and UoL momentum transfer modelling frameworks have been systematically validated in numerous recent publications (Colombo and Fairweather, 2015Liao et al., 2018;Lucas et al., 2020).

Interfacial forces
In the HZDR model the drag force, lift force, wall force, turbulent dispersion force and virtual mass are all considered: In the UoL model, only the drag, lift and turbulent dispersion forces are modelled. The wall force is neglected since wall-peaked void fraction profiles have been predicted without it when a Reynolds stress turbulence model with wall resolution capabilities is employed Fairweather, 2019, 2020), and due to recently reported drawbacks in the theoretical foundations of some wall force models (Lubchenko et al., 2018). The absence of the wall force is one of the two major differences between the two overall models, together with the multiphase turbulence modelling approach employed (more details are given in the turbulence modelling section). Also the virtual mass force is neglected in the UoL model but, for the steady fully-developed flows considered in this work, no significant impact is expected from its neglecting.
The drag force models the resistance exerted by the surrounding liquid on the bubble motion, and the corresponding momentum source to the liquid phase is given by: In Eq. (4), U r = U d − U c is the relative velocity between the bubbles and the liquid and d B the average bubble diameter. The drag coefficient C D needs to be calculated with a specific model and, at HZDR, the model of Ishii and Zuber (1979) is employed, where C D is expressed as a function of the bubble Reynolds number Re (Re = |U r |d B /ν mol c , where ν mol c is the liquid kinematic viscosity) and the Eötvös number Eo (Eo = Δρgd 2 B /σ, where Δρ is the density difference and σ the surface tension): At UoL, the drag coefficient is instead calculated from the model of Tomiyama et al. (2002a), which also accounts for the effect of the bubble aspect ratio E: In the above equation, F is also a function of the bubble aspect ratio (Tomiyama et al., 2002a). E is determined from the following expression (Hosokawa and Tomiyama, 2009): In Eq. (8), derived to reproduce experimental evidence of an aspect ratio close to 1 near solid walls (Hosokawa and Tomiyama, 2009), y w is the distance from the wall and the reference aspect ratio E 0 is calculated from the model of Welleck et al. (1966). The presence of neighbour bubbles altering the velocity field and the drag coefficient experienced by a bubble at high bubble concentration is tentatively accounted for with an additional correction factor (Hosokawa and Tomiyama, 2009): As mentioned in the introduction, in a shear flow the lift force expresses the force experienced by the bubble in the direction perpendicular to the main fluid motion (Auton, 1987): For a spherical bubble, the lift coefficient C L is positive and the bubble travels in the direction of lower liquid velocity. Instead, for large deformed bubbles, C L becomes negative. The HZDR model employs the correlation of Tomiyama et al. (2002b), derived from the experimental observation of single air bubbles rising in a glycerol-water solution: where Eo ⊥ is the Eötvös number based on the maximum horizontal dimension of the bubble as the characteristic length and This maximum dimension is calculated from the Welleck et al. (1966) aspect ratio correlation. The correlation of Tomiyama et al. (2002b) predicts the change of sign in the lift coefficient for air bubbles in water at d B ~ 6 mm in atmospheric conditions. In the UoL model, a constant value of C L0 = 0.10 is assumed, after improved accuracy was obtained with the model over a wide range of experiments with wall-peaked void profiles Fairweather, 2015, 2020). However, the UoL approach employs a near-wall turbulence model that requires fine mesh resolution near the wall (more details on the turbulence model are provided in the following section). Therefore, at a distance from the wall smaller than the bubble diameter, the lift coefficient is decreased to approach zero at the wall and avoid very high, unphysical values of the lift force in the very small cells adjacent to the wall (Shaver and Podowski, 2015): When a bubble approaches a solid wall, it experiences an additional wall lift force, driven by the modification of the flow field around the bubble that prevents the bubble from moving further towards the wall. This force has been often modelled with an additional lateral wall force (Rzehak et al., 2012): where n w is the vector normal to the wall pointing into the fluid. This force is only included in the HZDR model, where the wall coefficient is calculated using the model of Hosokawa et al. (2002), derived based on the trajectories of single bubbles in the range 2.2 ≤ Eo ≤ 22: where The turbulent dispersion force accounts for the effect of the turbulent fluctuations in the continuous phase on a bubble. In both models, the formulation of Burns et al. (2004), based on the Favre averaging of the drag force, is employed: In Eq. (16), ν turb c is the turbulent kinematic viscosity of the continuous phase and σ α the turbulent Prandtl number for the void fraction, taken equal to 1 in the UoL and 0.9 in the HZDR models. In the HZDR model, the virtual mass force F vm is also accounted for, using a fixed virtual mass coefficient C VM equal to 0.5 (Rzehak et al., 2017).

Turbulence modelling
To solve Eq. (2), the turbulent stress tensor τ Re k needs to be obtained from a turbulence model. Both HZDR and UoL model turbulence only in the continuous phase, HZDR using a two-equation model based on the eddy viscosity assumption: 2 is the strain rate tensor, μ turb c the turbulent dynamic viscosity and k c the turbulence kinetic energy of the continuous phase. UoL in contrast employs a Reynolds stress model that directly models the individual turbulent stress components by solving a transport equation for each: Both models are detailed below. Consideration of only the continuous phase turbulence is justified since in bubbly flows the dispersed phase has a much lower density and turbulent stresses are much higher in the continuous phase (Gosman et al., 1992;Rzehak and Krepper, 2013;Colombo and Fairweather, 2015). Therefore, in the UoL model the dispersed phase turbulence is derived from the continuous phase turbulence field, directly relating the turbulent viscosities of the two phases: where C t is a constant assumed equal to 1 (Behzadi et al., 2004). In the HZDR model, instead, a laminar flow is assumed for the dispersed phase, i.e. τ Re ij,d = 0 with negligible effects expected on the results (Rzehak and Krepper, 2013).

HZDR turbulence model
The HZDR model adopts the two-equation k-ω SST model, which combines the advantages of the k-ω model near the wall and the k-ε formulation away from it (Menter, 2009). The model solves an equation for the turbulence kinetic energy k and the turbulence frequencyω: together with standard k-ε model equations, where ε is the rate of dissipation of turbulence kinetic energy. The combination of the k-ω and k-ε models is achieved by a blending function, which is explicitly prescribed in terms of the wall distance as: This blending function is also used to interpolate the model constants ω between the corresponding values of the k-ω model (index '1') and k-ε model (index '2'). The usual values of the above constants for single-phase flows are applied as summarized in Table 1 at the end of the section. The production term for the shearinduced turbulence: includes a limiter to prevent the build-up of turbulence kinetic energy in stagnation zones. Since bubble-induced turbulence effects are included in k and ω due to the respective source terms S BI k and S BI ω discussed below, the turbulent viscosity is evaluated from the standard relation of the SST model: which includes a limiter with a second blending function: and a further model constant C γ = 1/0.31. In addition, a turbulent wall function is applied, as described in detail, for example, by .

UoL turbulence model
In the UoL model, turbulence is modelled with an elliptic-blending Reynolds stress model (EB-RSM) (Manceau and Hanjalic, 2002;Manceau, 2015) that has near-wall resolution capabilities. The model solves an equation for each turbulent stress − u i u j = τ Re ij /ρ c and the turbulence energy dissipation rate ε: Here, P ij is the production of turbulence due to shear, with . The turbulence dissipation rate ε is assumed isotropic in the bulk of the flow away from the wall, but in the near wall region it becomes a tensor, ε ij , as defined below. Φ ij is the pressure-strain correlation, which mainly redistributes the turbulence kinetic energy between the normal stress components. It is modelled following the SSG model of Speziale et al. (1991). The turbulent timescale T is equal to (18) and (19), S BI are the source terms for the bubble-induced contribution to the continuous phase turbulence.
In the elliptic-blending model, the correct near-wall behaviour of the turbulent stresses is achieved by blending a near-wall formulation with the SSG model in the bulk flow region, avoiding the need for any wall function (Manceau and Hanjalic, 2002;Manceau, 2015). Near the wall, the pressure-strain is modelled as: In the previous equation, n i are the components of the wall-normal unit vector. Blending from the near-wall behaviour to the bulk flow region model for the pressure-strain and the turbulence dissipation rate is achieved with a relaxation function α EB that is obtained by solving an elliptic relaxation equation (Manceau, 2015): In Eq. (29), L is the turbulence length scale given by L = . When required, such as in Eq. (19) to estimate the turbulence in the dispersed phase, the turbulent viscosity is calculated with the usual single-phase relation: Values of the numerous constants employed in the model are summarized in Table 1.

Bubble-induced turbulence modelling
A significant challenge in multiphase turbulence modelling is how to properly model the portion of the continuous phase turbulence generated by the bubbles, normally referred to as bubble-induced turbulence. Both the HZDR and UoL models include source terms for bubble-induced turbulence kinetic energy and dissipation rate in the turbulence model equations.
The bubble-induced turbulence kinetic energy production is derived from the approximation that the energy lost by the bubbles due to the drag force, F drag , is converted into turbulence kinetic energy in their wakes (Rzehak and Krepper, 2013): The numerical factor C BI k is set to unity for the HZDR model, following the original proposal of Rzehak and Krepper (2013). In the UoL model, a lower coefficient C BI k = 0.25 is introduced after the improved agreement obtained with the model over a wide range of bubbly flows (Colombo and Fairweather, 2015). Moreover in the UoL approach, for use with the Reynolds stress model, the source, once calculated, is divided between the three normal stresses, with a larger portion assigned to the axial stress along the direction of the mean flow: The source of turbulence energy dissipation rate is obtained from the turbulence kinetic energy source divided by a bubble-induced turbulence timescale τ BI : The bubble-induced turbulence timescale is modelled as in Rzehak and Krepper (2013) using the length scale of the bubble and the velocity scale of the continuous phase turbulence τ BI = dB k 0.5 . The C BI ε coefficient is taken equal to 1. For use with the ω -equation, the turbulence dissipation rate source is converted into an equivalent ω-source in the HZDR model.
A summary of the modelling closures employed in both models is provided in Table 2.

Numerical solution method
The HZDR model was solved in ANSYS CFX (ANSYS, 2019). In ANSYS CFX, the Navier-Stokes equations are solved with a control volume based finite-element discretization. In the present work, the advection terms are discretized using the high resolution scheme proposed in Barth and Jespersen (1989), while the solution is advanced in time with a second-order backward Euler scheme. The gas fraction coupling was achieved using the coupled solver option, and other details regarding the discretization of the diffusion and pressure gradient terms as well as the solution strategy are detailed in ANSYS (2019). Simulations were run in time until steady-state conditions were reached and this was evaluated by checking that values of velocity, void fraction and turbulence quantities showed variations in time of under 1% with respect to their mean values. In the experiments, measurements were taken at a sufficient distance from the inlet to avoid any flow development or inlet effects. Similarly, in the simulations results were recorded at the same distance from the inlet, sufficient for the velocity and void distributions to reach fully-developed conditions. At the wall, the no-slip boundary condition was imposed on the liquid phase, with the velocity in the first near-wall cell imposed using the single-phase wall law for a smooth wall. The free-slip boundary condition was instead imposed at the wall for the gas phase. At the inlet, the velocity of the phases and the void fraction were imposed based on the experimental measurements (adjusted if required, as will be discussed in the experimental data section). The average bubble diameter was kept fixed and was also taken from the experimental measurements. In this way, the development of the bubble diameter distribution before the measurement point, and any effect on it of the bubble injection method, can be neglected. A fixed pressure was imposed at the outlet section. Only a narrow axisymmetric section of each pipe was simulated and a mesh independence study ensured that grid independent solutions were achieved, and the distance from the wall of the first grid point was sufficient to ensure the validity of the law of the wall.
The UoL model was solved with STAR-CCM+ (CD-adapco, 2016). In STAR-CCM+, conservation equations are solved using a finite volume discretization. In this work, the second-order upwind scheme was used to discretize velocity, volume fraction, turbulence stresses and turbulence dissipation rate convective terms. The time derivative was discretized with a second-order implicit scheme and a multiphase extension of the SIMPLE algorithm (Patankar and Spalding, 1972) was used to solve the pressure-velocity coupling. Boundary conditions exactly matched those employed in the HZDR model, expect that the no-slip condition at the wall was also imposed on the gas velocity. For the EB-RSM model, at the wall a zero value was imposed on the turbulent stresses and the relaxation function α EB, and for the turbulence dissipation rate the asymptotic limit ε = 2ν mol c ( k/y w ) yw→0 was employed (Manceau, 2015). Clearly, inlet values of velocities and void fraction, and the fixed bubble diameters, were equal to those used for the HZDR model. Water and air properties were taken at a temperature of 25 • C and a pressure of 1 bar. A 1/4 section of each pipe was simulated and a sensitivity study ensured that mesh independent solutions were achieved.
For both models, sensitivity studies were performed by looking at changes in the water and air velocity, void fraction, turbulence kinetic energy, Reynolds stresses and turbulence frequency (or dissipation rate) radial profiles as a function of the mesh refinement. Mesh independence was considered achieved when negligible changes (of the order of 1-2% or lower) were observed with a further refinement of the mesh. The meshes employed are summarized in Table 3, where the total number of elements, the mesh elements in the radial and axial directions and their respective refinements are included. Meshes of the order of 10 4 were sufficient for the HZDR model (6800-31,520), while at least 10 5 elements were necessary for the UoL model (220,800-2,553,600). Even though the larger number of elements was partially due to the quarter pipe geometry employed, the UoL model still requires 5-10 times more elements for the same geometry, with an associated increase of computational time, due to the wall refinement requirements of the EB-RSM.

Experimental data
Over the years, numerous experimental studies have addressed the behaviour of bubbly flows in pipes. The database built for this work includes 16 experiments taken from the studies of Tomiyama (2009), Liu (1998), Lucas et al. (2005) in the MTLoop facility, built and operated at HZDR over the last few decades, and Liu and Bankoff (1993a). As discussed previously in the introduction, experiments were selected with wall-peaked void fraction distributions that could be sufficiently-well predicted using a monodispersed bubble distribution with a fixed value of the average bubble diameter. This allowed the study to focus exclusively on interfacial momentum transfer and multiphase turbulence closures, without considering changes in the bubble diameter distribution induced by breakup and coalescence. A summary of the experimental conditions and the averaged values employed in the CFD simulations is provided in Table 4.
In Hosokawa and Tomiyama (2009), measurements were taken in a vertical upward air-water bubbly flow at atmospheric pressure and temperature in a pipe of inside diameter 25 mm. Radial profiles of gas volume fraction, liquid and gas velocity and liquid turbulence kinetic energy were measured using laser Doppler velocimetry and shadowgraphy at an axial location L/D = 68. Bubble concentration, size and shape were reconstructed from stereoscopic images obtained with two high-speed cameras. The average measured bubble diameter was used in the CFD simulations, and the superficial velocities and averaged values of the void fraction over the pipe cross-section were imposed as the inlet conditions. Liu (1998) studied vertical upward air-water bubbly flows in a pipe of inside diameter 57.2 mm, at atmospheric pressure and a temperature of 26 • C. A dual resistivity probe and a single hot film anemometry probe were used to measure radial profiles of the liquid velocity and turbulence intensity, the gas volume fraction and the average bubble diameter Table 3 Parameters of the meshes employed for the simulations, including the total number of mesh elements N, the number of elements in the axial and radial directions N z and N r and the refinement in the axial and radial direction n z and n r . For HZDR, total number refers to a narrow axisymmetric section. For UoL, total numbers refer to a 1/4 section of the pipe, and the range is provided as a minmax range. at L/D = 60. From integration of the void fraction radial profiles, averaged values were obtained and imposed at the inlet in the CFD simulations. The averaged void fractions were also used to correct the value of the air superficial velocity to achieve the correct flux of air through the pipe cross-section at the measurement position. The MTLoop facility (Lucas et al., 2005) was built at HZDR and employed to study the development of upward vertical flows of air and water in a pipe of inside diameter 51.2 mm using the wire-mesh sensor technique. Radial profiles of the gas average velocity and volume fraction, and the bubble size distribution, were measured at different heights from the inlet up to L/D = 60, at atmospheric pressure and 30 C temperature. Measurements with a bubble size distribution almost constant along the axial direction were selected, and average bubble diameter and void fraction from the last measuring station used to setup the CFD simulations. The average void fraction was also used to adjust the nominal value of the gas superficial velocity. Liu and Bankoff (1993a) studied upward air-water bubbly flows in a vertical pipe of inside diameter 38 mm at atmospheric pressure and temperature conditions. Measurements were taken at L/D = 36, and the liquid velocity was measured using one-and two-dimensional hot-film anemometer probes, while void fraction and bubble velocity and frequency were obtained using an electrical resistivity probe. The measurements cover a large range of flow conditions and include radial profiles of liquid and gas velocities, turbulence levels, void fraction and bubble diameter. Provided values of superficial velocities and average void fraction, and bubble diameter, were used to setup the CFD simulations.

Results and discussion
4.1. Hosokawa and Tomiyama (2009) Predictions of the Hosokawa and Tomiyama (2009) experiments are summarized in Fig. 1 (for cases H11 and H12) and Fig. 2 (for cases H21  and H22). For these experiments, measurements are available for liquid velocity, relative velocity (between the bubbles and the liquid), void fraction and turbulence kinetic energy (although not shown here, data are also available for the individual normal turbulent stresses). Here, and in all the following figures, comparisons are made against radial profiles of the physical quantities measured as a function of the nondimensionalized (by the pipe radius) radial distance from the pipe centreline, with 0 being the pipe centreline and 1 the pipe wall. Given that all the experiments considered are for air bubbles in water, the subscripts w and a will be used in the following to identify the two phases.
Good agreement is achieved by both models for the liquid mean velocity profiles as shown in parts (a) and (e) of Figs. 1 and 2. The UoL model, with the near-wall refinement required by the EB-RSM model, matches almost exactly the liquid velocity decrease near the wall. Away from the wall, in contrast, it is the HZDR model that provides the best prediction, showing a remarkable accuracy in the centre of the pipe for all four experiments, with an average relative error of 3% and always lower than 4.5%. The UoL model consistently overestimates the liquid velocity away from the wall, even though the discrepancy is always less than 10% and 7.5% on average. The reason for this is found in the relative velocity plots, in Figs. 1(b) and (f) and 2(b) and (f). Near the wall, the UoL model predicts well the decrease in the relative velocity, induced by the higher drag of the more spherically shaped bubbles in this region, which is only partially captured by the HZDR model. However, outside the near-wall region, and despite the low spatial resolution of the measurements in the centre of the pipe, the UoL model tends to underpredict the relative velocity, with the HZDR approach found to be in better agreement with data. A lower relative velocity is induced by a higher drag coefficient. Therefore, the mentioned overestimation of the liquid velocity by the UoL model can be explained with the excessive drag from the bubbles to the liquid predicted by the drag model employed.
Both models provide robust predictions of the void fraction distribution shown in parts (c) and (g) of Figs. 1 and 2, with marked wallpeaked radial void fraction profiles and a lower void fraction concentration in the centre of the pipe. Notable discrepancies with data are found in the pipe centre in experiment H12 which has the lowest liquid velocity and the highest void fraction. In these conditions, larger bubbles (H12 has indeed the largest average bubble diameter of the four experiments) may form and migrate towards the pipe centre, increasing the void fraction there. The UoL model, which uses a constant positive value of the lift coefficient, is unable to capture this behaviour. In the HZDR model, the Tomiyama et al. (2002b) correlation predicts the change in the sign of the lift coefficient. However, this happens at d B ≈ 5.8 mm and, until d B ≈ 4.25 mm (the measured averaged bubble diameter was d B = 4.1 mm), the model returns an almost constant positive lift coefficient equal to 0.28. Therefore, the HZDR model, in the present "monodispersed" configuration, is also unable to entirely capture the void fraction profile in H12.
The UoL model, despite not using any wall force, and avoiding all the related uncertainties, shows good predictions of the void peak position and magnitude, the latter being predicted with an average relative error of 20%. In a turbulent flow, the radial turbulent stress is not constant and induces a radial pressure gradient that shows a minimum around the location of the void peak. This gradient in the stress, which is properly resolved by the Reynolds stress turbulence model, contributes to the lateral void fraction distribution by pushing bubbles towards the minimum pressure region. From this minimum, the pressure increases again towards the wall and this increase, predicted by the EB-RSM near-wall   model, is sufficient to predict the void peak in the simulations without any additional wall force Fairweather, 2019, 2020). Reasonable accuracy for the void fraction peak is also shown by the HZDR model, although the model tends to predict an excessive bubble accumulation near the wall, visible in the somewhat overpredicted peak magnitude in Figs. 1 and 2. Nevertheless, it has also to be pointed out that the relative error on the peak may be affected by uncertainty related to the discrete nature of experimental measurements. Therefore, the highest value measured may not be exactly the peak value, and this can also contribute to explaining why the models tend, here and in the following experiments, to predict a higher peak. Lastly, Figs. 1(d) and (h) and 2(d) and (h) show the turbulence kinetic energy. Very good and similar agreement is found for cases H21 and H22 (Fig. 2(d) and (h)), where the maximum relative error on the centreline is 25% for the UoL model in H22. These cases have the highest liquid velocity and a very low void fraction concentration in the pipe centre where, therefore, the turbulence production is mainly sheardriven. In contrast, cases H11 and H12, with a smaller decrease of the turbulence kinetic energy away from the wall, show a more significant contribution of the bubble-induced turbulence and more evident differences between the models and the experiments. The HZDR model better predicts case H11 while the UoL approach is superior for case H12, although neither is particularly accurate for the latter, with errors as high as 50%. In the wall region, the EB-RSM in the UoL model reproduces well the behaviour of the turbulence kinetic energy and its near-wall peak.

Liu (1998)
Compared to Hosokawa and Tomiyama (2009), data from Liu (1998) were measured in a larger pipe and at significantly higher void fractions. Measurements of liquid velocity and void fraction are available and comparisons with CFD predictions are provided in Fig. 3, together with the turbulence kinetic energy, estimated from the axial turbulent normal stress measurements assuming its value is two times that of the radial and angular normal stresses, as always observed in pipe flows outside of the near-wall region (Rzehak and Krepper, 2013;Colombo and Fairweather, 2015).
Velocity profiles, as a consequence of the higher void fraction and the still marked near-wall peak ( Fig. 3(b), (e), (h) and (k)), are considerably flatter than observed previously. Overall, good predictions are still obtained, with average relative errors on the centreline of 5.6% for the UoL and 4.2% for the HZDR models. The CFD models have a tendency to predict flatter profiles with respect to the experiments, in particular for cases L11A (Fig. 3(a)) and L21C ( Fig. 3(g)). This is likely due, at the high void fractions considered, to larger bubbles travelling in the centre of the pipe that are not captured by the models in the present configuration. The UoL model even predicts a slight peak at the wall in the velocity profile, which is likely caused by the excessive drag from the air bubbles predicted with the Tomiyama et al. (2002a) model. Near the wall, the calculated velocity gradients are much steeper than the measured ones. Both models remain in reasonable agreement with data, although the HZDR model predicts the water velocity slightly better. It is possible that, despite the finer resolution of the EB-RSM, the model, still based on a single-phase formulation, does not capture some two-phase effects induced by the high void fraction. Another possible reason is the freeslip condition imposed on the air bubbles in the HZDR model, which may trigger higher air and water velocities near the wall. For a more precise assessment, and any future developments, availability of detailed measurements in the near-wall boundary layer is a priority.
Predictions of the void fraction remain robust, although both models have a tendency to overpredict the near-wall peak. The UoL model has the best agreement with data, and maintains good accuracy for both the peak magnitude and position. On the other hand, the HZDR model once again shows a tendency to predict an excessive accumulation of void near the wall. The wall force model from Hosokawa et al. (2002) is employed in the latter, which has been proven to have a greater impact than other formulations that assume a linear decrease of the wall force with distance from the wall (Rzehak et al., 2012). However, its effect can still be too weak, contributing to the overestimated accumulation of bubbles near the wall. In the pipe centre, both models demonstrate remarkable accuracy and maximum deviations from the data are limited to a few percent, with the exception of L21C in Fig. 3(h). This case, similarly to H12 for Hosokawa and Tomiyama (2009), has the highest average bubble diameter in the group, and is the only one where this is greater than 4 mm (d B = 4.22 mm). Therefore, it is again plausible that in this case larger bubbles flow in the centre of the pipe which are not resolved in the simulations. For the other three experiments, the void fraction on the centreline is predicted with an average relative error of 10.8% by the UoL and 7.2% by the HZDR models.
The largest discrepancies are again found in the turbulence kinetic energy comparisons, confirming the complexity of predicting turbulence in flows that contain a contribution from the bubbles. The two models differ only by a coefficient, introduced in the UoL model to limit the turbulence kinetic energy source. Therefore, the HZDR model always returns the highest turbulence kinetic energy between the two models. With respect to the experiments, mixed results are obtained, with the HZDR model better in L21C (Fig. 3(i)), the UoL model in L21B (Fig. 3(f)) and neither able to properly predict L11A (Fig. 3(c)) and L22A (Fig. 3  (l)). On the centreline, the relative error varies from less than to 2% for the UoL model in L21B to values as high as 50-100%. This suggests further developments are needed, specifically improving on the constant coefficients employed in the bubble-induced source. In the near-wall region, the EB-RSM is more accurate, although less clearly than in the case of the Hosokawa and Tomiyama (2009) experiments. In L11A (Fig. 3(c)) and L22A (Fig. 3(l)), the UoL model is closer to the experimental peak in the turbulence kinetic energy. However, in one of the two other cases where the HZDR model performs better in this respect (L21B, Fig. 3(f)), this seems to be more a consequence of a general overprediction of k across the entire pipe. Still, as observed for the water velocity, the more limited improvement at high void fractions suggests that there are relevant two-phase effects that the present EB-RSM singlephase based formulation is still not able to capture.

MTLoop
With the MTLoop experiment (Lucas et al., 2005), the focus is back to low void fraction cases, but in a larger pipe than used by Hosokawa and Tomiyama (2009). Comparisons against the four experiments for air velocity and void fraction profiles can be found in Fig. 4.
Largely, predictions of the air velocity are in very good agreement with the experiments (Fig. 4(a), (c), (e) and (g)), with average relative errors on the centreline lower than 2.5% for both models. Results from the UoL model are always lower than for the HZDR model, and always on the lower side of the measurements, confirming the excessive drag (slightly in this case) that results in lower relative velocities. Near the wall, the HZDR model is in line with the experiments, while the UoL model bubble velocity reduces excessively approaching the wall. Again, the HZDR approach employs a free slip boundary condition, whilst no slip is imposed in the UoL model. Therefore, the free slip boundary condition appears to be most appropriate for the gas phase.
Good agreement is found for the void fraction, both in terms of the peak near the wall and in the pipe centre ( Fig. 4(b), (d), (f) and (h)). The UoL model predictions confirm their previously observed accuracy, in particular for the void peak, which is predicted with an average relative error of 20%. For the HZDR model, the only experiment where the previously noted tendency to overpredict bubble accumulation at the wall is of a noticeable extent is MT041 (Fig. 4(d)). The underprediction of data in the pipe centre for case MT063 (Fig. 4(h)) is most probably due to the already discussed presence of larger bubbles in the central region, with the average measured bubble diameter being 5.2 mm for this case.     Liu and Bankoff (1993a) Last to be discussed are the experiments from Liu and Bankoff (1993a). These are most similar to those of Liu (1998) in terms of void fraction, but were performed in a smaller diameter pipe. Results for the liquid and air velocity, as well as void fraction, are provided in Fig. 5, and those for turbulence quantities in Fig. 6. Fig. 5 further confirms previous findings, particularly those from the cases of Liu (1998). Calculated velocity profiles, although generally in agreement (average relative errors on the centerline between 5% and 6% for the water and 7% for the air), are flatter than their experimental counterparts when the void fraction in the pipe centre is high, such as for experiments LB16 and LB17 (Fig. 5(a) and (d)). Drag in the UoL model is excessive, given that the air velocity is lower than for the HZDR model and the experiments, even when the corresponding water velocity is higher (LB30 and LB31, Fig. 5(h) and (k)); in these experiments, where the void fraction in the pipe centre is very low, flat velocity profiles are not obtained. In experiment LB17 (Fig. 5(d)), this causes the water velocity to peak near the wall, sometimes also found in the results of Fig. 3. In experiments LB30 and LB31, both models underpredict the air velocity in the pipe centre, although the experimental relative velocity values are substantially higher (greater than0.3 ms − 1 ) than those observed in all other experiments considered. Similar high values have only been reported for bubbles of smaller diameter rising in ultrapurified water which points to some undetected uncertainties in the measurements as a plausible cause for the noted discrepancies (Kriebitzsch and . Near the wall, it is confirmed that some of the advantages of the EB-RSM, due to a combination of two-phase effects in the liquid phase turbulent wall function and the no-slip condition for the gas phase, are lost at high void fraction, with the HZDR model performing better for this dataset, although some spikes in the air velocity near the wall require further verification.

4.4.
The void fraction shows an overall good agreement. Between the two models, the UoL model still shows the best agreement. For cases LB30 and LB31 (Fig. 5(i) and (l)), the UoL model still overpredicts the peak magnitude, but always returns the correct peak location. In the pipe centre, no significant differences with experimental data are apparent.
Turbulence measurements are available for the streamwise and the radial root mean square (r.m.s.) of the velocity fluctuations and the Reynolds shear stress ( Fig. 6(a), (c), (e) and (g)). Having been calculated from the isotropic assumption, the r.m.s. values predicted by the HZDR model are always between the two sets of experimental values. For LB16 and LB17, where bubble-induced turbulence is significant, the HZDR model performs best, while the UoL model underpredicts the turbulence levels. This trend is more marked than in previous cases, where inconclusive results were found. The UoL model also shows smaller differences between the streamwise and radial components than found in the experiments, possibly confirming some recent findings that the major contribution from the bubble-induced turbulence source should act in the streamwise direction (du Cluzeau et al., 2019;Ma et al., 2020). In LB30 and LB31, where the void fraction in the pipe centre is low and turbulence dominated by the shear contribution, predictions are in close agreement with measurements.
At high void fraction, both models are unable to predict the Reynolds shear stress except for in the very near-wall region. However, it has to be remarked that, although the velocity profiles for LB16 and LB17 look flatter than for typical single-phase flow behaviour, the shear stresses instead are far from flat and appear more similar to their single-phase counterparts. Therefore, additional testing is recommended in this area. Unfortunately, however, Reynolds shear stress measurements are rarely available from bubbly flow experiments, and additional data is highly desirable. In cases LB30 and LB31 (Fig. 6(f) and (h)), the UoL model is in reasonably good agreement with data, and predicts well measurements from r/R = 0.5-0.6 towards the pipe centre. The eddy viscosity assumption employed in the HZDR model, in contrast, remains inaccurate in this region. In these experiments, the closest measurements are still at a considerable distance from the wall. However, it is worth noting that the UoL model shows in all cases good agreement with the measurement point closest to the wall.

Conclusions
The accuracy of modelling closures for interfacial momentum transfer and multiphase turbulence employed and developed at the Helmholtz-Zentrum Dresden-Rossendorf (HZDR) and the University of Leeds (UoL), embodied in overall Eulerian-Eulerian computational fluid dynamic models of two-phase flows, has been assessed against experimental data for monodispersed bubbly pipe flows. Overall, both models demonstrate robustness and accuracy over the wide range of operating conditions considered, and can be employed with a considerable degree of confidence up to averaged void fractions of 10% or even higher, provided that the flow still exhibits a distinctive monodispersed behaviour. In the range 10%-20% void fraction, additional modelling such as a population balance model is necessary when the flow starts exhibiting polydispersed features with larger bubbles flowing in the centre of the pipe.
Velocity profiles were generally well-predicted, with average relative errors on the centreline of 6.5% (UoL) and 4.3% (HZDR) for the water and 4.8% for the air, and more than 90% of the data were predicted with a relative error lower than 10%. The HZDR drag model (Ishii and Zuber, 1979) is slightly more accurate than that of UoL (Tomiyama et al., 2002a), which underpredicts the relative velocity, impacting both liquid and gas velocity profiles. However, the capability of the latter to account for the bubble aspect ratio increase near the wall should be maintained in future models. Comparisons also show that the free slip boundary condition at the wall has to be imposed on the gas phase rather than the no slip condition.
The wall-peaked void fraction profiles, typical of these flows, were successfully and consistently predicted. However, the lift-wall force combination of the HZDR model tends to overpredict bubble accumulation at the wall, while the UoL model shows consistency in predicting the peak position and magnitude, even without a wall force term. Overall, the average relative error on the peak magnitude is 25%. Additional wall effects, and a non-constant lift coefficient, will however certainly be necessary to predict polydispersed flows or laminar conditions. On the centreline, it is more difficult to provide accurate figures for the void fraction, given that some values of the void fraction are extremely low (minimum deviations result in high relative errors) and some experiments were affected by the presence of large bubbles. However, considering only the experiment with a value of the void fraction higher than 0.005, and excluding those not showing monodispersed features (H12, L21C and MT63), the average relative error is 26% for the UoL and 29% for the HZDR models.
Bubble-induced turbulence is the area in need of major development, with none of the models employed being consistently accurate. On the centreline, the relative error for the turbulence kinetic energy can be as low as 2%, but is below 50% for only half of the data. It is also the area where more new findings and additional physical understanding were made available in recent publications (Ma et al., 2017(Ma et al., , 2020Magolan et al., 2017;du Cluzeau et al., 2019;Magolan and Baglietto, 2019). Improvement of the current formulations, which are limited to a drag force contribution multiplied by a constant modulation coefficient, is the most urgent need, with HZDR already having done relevant work in this area (Ma et al., 2017;Liao et al., 2019).
Overall, the Reynolds stress turbulence closure introduces the additional effect of a radial pressure gradient, from which models that are of general applicability should benefit. The development and use of nearwall treatments should be encouraged, given the improved accuracy achievable for velocity and turbulence quantities. At high void fraction, however, two-phase specific models, in place of the mostly single-phase formulations currently available, are required, in the absence of which any advantage over high-Reynolds number treatments is not guaranteed.
Building on the present results, future efforts will also aim at a harmonized best possible model that will include near-wall and turbulence modelling from UoL and the HZDR momentum exchange closure framework. The latter demonstrated better overall performance and is already equipped to predict the change of sign in the lift force and extend the model to polydispersed flows. The coupled model accuracy will be re-evaluated and sensitivity studies to quantify the impact on the overall accuracy of any change in each closure will be necessary, to focus better future development efforts.
Finally, the availability of more detailed experimental data, such as measurements with high resolution in the near-wall region, individual components of the Reynolds stress tensor, and resolved polydispersity of the bubble size distribution, needs to be improved to support the further development of these models. In addition, further works of this kind, as well as the identification of proven sets of measurements accepted community-wide, over which models can be quantitatively benchmarked and judged, will support the confidence in, and increased applicability of, multiphase CFD models.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.