Multi-Fluid Computational Fluid Dynamic Predictions of Turbulent Bubbly Flows Using an Elliptic-Blending Reynolds Stress Turbulence Closure

The accurate prediction of bubbly flows is critical to many areas of nuclear reactor thermal hydraulics, mainly, but not only, in relation to the key role bubble behavior plays in boiling flows. Large scale computations of flows with hundreds of thousands of bubbles are possible at a reasonable computational cost using computational fluid dynamic, multi-fluid Eulerian-Eulerian models. The main limitation of these models is the need to entirely model interfacial transfer processes with proper closure relations. Here, the capabilities and advantages provided by a model that includes an elliptic-blending Reynolds stress turbulence closure (EB-RSM), allowing fine resolution of the velocity field in the near-wall region, are tested over a large database. This database includes mostly monodispersed bubbly flows over a wide range of operating conditions and geometrical parameters, including upward and downward pipe flows, large diameter pipes and a square duct. The model shows encouraging accuracy and robustness, with good agreement over most void fraction distributions and accurate prediction of the magnitude and position of the near-wall void fraction peak. The model does not include any wall force, avoiding all the related uncertainties, and the prediction of the void fraction peak relies on the fine resolution of the near-wall pressure gradient induced by the turbulence field. Overall, the EB-RSM allows accurate resolution of the velocity and turbulence field near the wall, and the transition to this and similar turbulence closures is of value in assisting the ongoing quest for thermal hydraulic models that are accurate and of general applicability. Additional modifications to the near-wall modeling approach, which is still based on its single-phase counterpart, may be required to deal with high void fraction conditions and, in the overall model, additional improvements to momentum and, most importantly, bubble-induced turbulence closures are desirable.


INTRODUCTION
Bubbly flows are frequent in nature and in numerous industrial and engineering applications. Bubbles dispersed in a continuous liquid promote mixing in the liquid phase and large interfacial areas in industrial processes that cause high rates of heat and mass transfer between the phases or the different fluids/components (Lehr et al., 2002;Risso, 2018). In nuclear thermal hydraulics, bubbly flows have a critical role in boiling, and therefore impact the operation, safety limits and accident response of water-cooled reactors. Detailed knowledge of the behavior and size distribution of bubbles is essential for the accurate prediction of boiling, three dimensional void fraction distribution and evolution of the flow regime in boiling flows. Given its impact on nuclear reactor thermal hydraulics, a significant amount of recent research has been focused on improving our understanding of, and predictive capability for, gas-liquid bubbly flows (Mimouni et al., 2010(Mimouni et al., , 2015Hosokawa et al., 2014;Colombo and Fairweather, 2015;Hazuku et al., 2016;Lucas et al., 2016;Sugrue et al., 2017;Liao et al., 2018;Lubchenko et al., 2018).
The presence of two or more phases complicates the physics of the flow and their multiple interactions make modeling particularly challenging. The liquid phase alters the motion and distribution of bubbles, which mutually impact the continuous phase flow in multiple ways (Liu and Bankoff, 1993a,b;Feng and Bolotnov, 2017). Bubble size distribution continuously evolves, driven by collision and coalescence between bubbles, which can also breakup following interactions with the continuous fluid phase (Liao et al., 2015;Colombo and Fairweather, 2016;Liu and Hibiki, 2018). Most of the time, these processes that impact the large scale fluid behavior are governed by phenomena at much smaller scales (Prince and Blanch, 1990;Legendre and Magnaudet, 1997;Martinez-Bazan et al., 1999;Lucas, 2009, 2010;Bolotnov, 2017, 2018). This multiscale complexity has limited the accuracy achievable with empirical or one-dimensional modeling approaches, and is driving the development of computational fluid dynamic (CFD) models given their increased ability to account for small scale physical effects on the three-dimensional fluid motion (Bestion, 2014;Liao et al., 2018;Podowski, 2018).
Recently, a significant research effort has been dedicated to the development of more advanced CFD models for bubbly flows. Significant improvements have been made in interface tracking/resolving methods which are able to resolve all the interfacial details of each individual bubble and that are helping to improve understanding of many still poorly known physical details (Santarelli and Fröhlich, 2016;Feng and Bolotnov, 2017;Magolan et al., 2017;Mehrabani et al., 2017;Chen et al., 2019;du Cluzeau et al., 2019). However, these methods are computationally expensive and remain prohibitive when more than a few hundred bubbles are present. Therefore, for industrial flows and applications where hundreds of thousands of bubbles may be present, Eulerian-Eulerian multi-fluid models remain the preferred choice (Hosokawa and Tomiyama, 2009;Rzehak and Krepper, 2013;Colombo and Fairweather, 2015;Mimouni et al., 2017;Liao et al., 2018). In these approaches, field equations are averaged and interfacial transfer processes entirely modeled by means of often empirically-based closure relations that clearly impact the overall accuracy and applicability of the models. At present, no general agreement has been found on the best closure models available, and the use in these of multiple adjustable constants, eventually optimized on a case-by-case basis over limited experimental databases, remains a major constraint (Lucas et al., 2016;Podowski, 2018).
Momentum transfer governs the dynamic interaction between bubbles and the liquid phase and is modeled by introducing interfacial force terms. In bubbly flows, the lift force has received special attention, with it being the major driver of the bubbles' transverse motion and the accumulation of near spherical bubbles near the wall, and large cap bubbles in the center, of closed ducts. Near the wall, the addition of a wall force has been used to prevent bubbles moving below a certain distance from the wall. In the literature, general agreement on a unified formulation has not been reached and numerous liftwall formulations exist that often differ only in the value of some model coefficients (Colombo and Fairweather, 2019). In addition to this uncertainty, not only the validity of some of the modeling assumptions, but the entire existence of the wall lubrication force, essentially introduced to predict the correct wall-peaked nearwall void fraction profiles in closed ducts, has been recently questioned. Rzehak et al. (2012) have questioned the linear decrease of the wall force with the distance from the wall, and the accuracy of the widely used model of Antal et al. (1991). More recently, Lubchenko et al. (2018), have predicted the wallpeaked void profile in bubbly flows with a modified expression for the turbulence dispersion force, but without the direct action of any wall lubrication force. The authors' work started from the observation, supported by measurements (Hassan, 2014) and interface tracking results (Lu and Tryggvason, 2013), that the liquid film that, in wall lubrication theory, remains between the bubble and the wall, is either negligibly small or absent.
Another open area in multi-fluid modeling is the development of multiphase turbulence closures. At the present time, the most often adopted strategy consists in adding specific source terms to the turbulence model equations to account for bubble-induced turbulence. Although significant advances have been made in recent years, particularly in modeling based on the conversion of energy from drag to turbulence kinetic energy in bubble wakes (Troshko and Hassan, 2001;Rzehak and Krepper, 2013;Ma et al., 2017;Magolan et al., 2017), turbulence modeling still often relies on the eddy viscosity assumption (Yao and Morel, 2004;Rzehak and Krepper, 2013;Sugrue et al., 2017;Liao et al., 2018). In contrast, second moment closures have only been applied in a few studies (Lopez de Bertodano et al., 1990;Lahey et al., 1993). The necessity for moving beyond the limiting assumptions of eddy viscosity based approaches is welldocumented in single-phase flows (Benhamadouche, 2018), and similar arguments can be made for multiphase flows. Using Reynolds stress based closures, it has been demonstrated that turbulence and its modeling can have an additional impact on the distribution of the dispersed phase (Ullrich et al., 2014;Santarelli and Frohlich, 2015) and, as an example in nuclear thermal hydraulics, in the prediction of flow rotation effects induced by mixing blades in pressurized water reactor coolant channels (Mimouni et al., 2017). Advances in bubble-induced turbulence modeling for implementation in Reynolds stress closures have also recently been achieved (Colombo and Fairweather, 2015;Parekh and Rzehak, 2018).
In a recent paper, we made a preliminary assessment of the additional impact that flow turbulence has on the bubble distribution through the pressure gradient in the fluid phase generated by the anisotropic turbulence field (Colombo and Fairweather, 2019). This was captured using an elliptic-blending Reynolds stress model (EB-RSM) implemented in the STAR-CCM+ CFD code that, to the best of the authors' knowledge, represented the first application of near-wall modeling to bubbly flows. Improved near-wall predictions of multiple quantities were achieved and, through the action of the pressure gradient induced by the turbulence field, accurate prediction of the usual near-wall peak in the void fraction distribution was found without the need for any additional wall force. In the present work, the accuracy of the same multi-fluid model coupled with the EB-RSM is tested against a large database of air-water bubbly flows and its assessment extended over a significantly wider range of conditions. The EB-RSM is coupled with a set of momentum transfer interfacial closures previously assessed with other turbulence models and a specific source for bubble-induced turbulence modeling and, initially, this model is compared against a high-Reynolds number RSM that was previously tested against a similarly large database (Colombo and Fairweather, 2015). Given the already mentioned abundant availability of slightly different and not extensively validated models (Lucas et al., 2016), assessment against a large database including a wide range of conditions is particularly important. The databases used in this paper therefore include, in addition to commonly used upward flows in pipes, other flows less frequently or never tested such as downward flows, large pipes and a square duct. To focus on momentum transfer and turbulence closure effects without adding the additional complication of a population balance model, some effort is taken to limit as much as possible the database to monodispersed bubble diameter distributions that can be characterized reasonably well by a single value of the average bubble diameter. The advantages of RSM closures and near-wall modeling in comparison to more standard methodologies in the critical area of reactor thermal hydraulics are highlighted, and areas for further improvements identified and discussed.

COMPUTATIONAL FLUID DYNAMIC MODEL
The CFD model employed is a two-fluid Eulerian-Eulerian model, where a set of averaged conservation equations is solved for each phase. Given the focus on air-water bubbly flows, only continuity and momentum balances are solved, with the phases treated as incompressible with constant properties: In the above equations, U is the velocity and p the pressure. τ and τ Re are the laminar and turbulent stress tensors, respectively, and g is the gravitational acceleration. ρ k and α k are the density and the volume fraction of phase k although, for simplicity, in the following we will refer to α only to identify the void fraction of the gas phase. The last term on the right hand side of Equation (2) is the interfacial momentum transfer source M k and includes drag, lift and the turbulent dispersion force. Wall lubrication, following wall-peaked void fraction profiles obtained even in the absence of any wall force when a Reynolds stress turbulence model is used (Ullrich et al., 2014;Colombo and Fairweather, 2019), and the reported uncertainties in its theoretical foundation (Lubchenko et al., 2018), is neglected. Also, the virtual mass force is neglected in view to its negligible effect in the considered multiphase flow conditions (Politano et al., 2003;Rzehak and Krepper, 2013;Colombo and Fairweather, 2015). The drag force expresses the resistance opposed to bubble motion by the surrounding liquid and is modeled using the correlation of Tomiyama et al. (2002a): The model is a function of the Eötvös number (Eo = ρgd B / σ , where σ is the surface tension and d B the diameter of the bubbles) and accounts for the effect of the bubble aspect ratio E on the drag coefficient C D : Experimental evidence shows that the bubble aspect ratio and drag coefficient increases near the wall, causing a reduction in the relative velocity between the bubbles and the fluid (Hosokawa and Tomiyama, 2009). Consistently with this, in Equation (4) the aspect ratio tends to a value of 1 (spherical bubble) at the wall and reduces with distance from the wall, y w , to approach a reference value E 0 calculated from the correlation of Welleck et al. (1966). In Equation (3), F is also a function of the bubble aspect ratio (Tomiyama et al., 2002a). In a shear flow, each bubble experiences a lift force perpendicular to its direction of motion that impact the lateral movement of the bubble and the void fraction distribution. Overall, the lift coefficient is positive for spherical bubbles and pushes them in the direction of lower liquid, and higher relative, velocity, i.e., toward the wall in upflow, resulting in typical wall-peaked void fraction profiles, and toward the center of the duct in downflow. For larger and more deformed non-spherical bubbles, the lift force reverses its direction, producing a shift in void fraction profiles from wall-peaked to core peaked. as observed in numerous upflow experiments (Tomiyama et al., 2002b;Lucas et al., 2005Lucas et al., , 2010. Over the years, agreement with experiments has been reported using values of the lift coefficient ranging from 0.1 (Wang et al., 1987;Yeoh and Tu, 2006) to 0.5 (Mimouni et al., 2010), and there remains a lot of uncertainty on the best possible model available. Numerous authors have used the correlation of Tomiyama et al. (2002b) that also predicts the change in sign of the lift coefficient. In this work, a constant value C L = 0.1 is adopted, following good agreement over a wide range of experimental conditions found by Colombo and Fairweather (2015) and in the work of other authors (Lopez de Bertodano et al., 1994;Lahey and Drew, 2001). Overall, the choice of a constant positive value of the lift coefficient will limit the applicability of the model to polydispersed bubbly flows, but is acceptable in the present work where only experimental conditions exhibiting a positive lift force coefficient were selected. In addition, as a consequence of the much more refined resolution near the wall required by the EB-RSM, very high lift values would be predicted in the small cells adjacent to the wall and at a distance much smaller than the bubble diameter. Therefore, and in the absence of a physically based approach, the lift force is decreased at a distance from the wall lower than the bubble diameter, to approach zero at the wall (Shaver and Podowski, 2015): The model of Burns et al. (2004) is used to model the turbulent dispersion force: In Equation (6), U r is the relative velocity between the phases, ν t,c the turbulent kinematic viscosity of the continuous phase and σ α the turbulent Prandtl number for the volume fraction, assumed equal to 1.0.

Multiphase Turbulence Modeling
Turbulence is resolved in the continuous phase only using the EB-RSM (Manceau and Hanjalic, 2002;Manceau, 2015), based on an extension of the single-phase transport equations for the Reynolds stresses R ij = τ Re i,j /ρ c (CD-adapco, 2016): (7) Here, P ij is the turbulence production due to shear and D ij the diffusion of the turbulent stresses, modeled following Daly and Harlow (1970). The turbulence energy dissipation rate ε ij is modeled following the isotropic hypothesis. Φ ij is the pressurestrain correlation and accounts for the redistribution of the turbulence kinetic energy between the stress components. Away from the wall, the pressure-strain is modeled using the "SSG model" model (Speziale et al., 1991), which is quadratically nonlinear in the turbulence anisotropy tensor a ij : In Equation (8), S ij and W ij are the strain and the rotation rate tensors. In the elliptic-blending formulation (Manceau and Hanjalic, 2002;Manceau, 2015), the SSG model is blended with a near-wall model that predicts the correct asymptotic behavior of the turbulent stresses near the wall. In this way, the flow can be entirely resolved in the near-wall region, without the need of imposing the velocity via a wall function. The near-wall model for the pressure-strain reads: ε k u i u k n j n k +u j u k n i n k − 1 2 u k u l n k n l n i n j +δ ij (9) In Equation (9), n are the components of the wall-normal vector. Blending between the near-wall and the SSG model is achieved using the elliptic relaxation function α EB , calculated by solving an elliptic relaxation equation with the α EB = 0 wall boundary condition: Here, L is the turbulence length scale given by L = C l max (Cην 3/4 ε −1/4 ,k 3/2 ε −1 ). Near-wall blending, in addition to the pressure-strain, is also enforced for the turbulence energy dissipation rate: The bubble-induced turbulence contribution to the continuous phase turbulence field is included by assuming that the energy lost by the bubbles to drag is converted into turbulence kinetic energy in the bubble wakes (Troshko and Hassan, 2001;Rzehak and Krepper, 2013). Source terms can accordingly defined for the turbulence kinetic energy and dissipation rate transport equations: Here, F d is the drag force and τ BI the timescale of the bubbleinduced turbulence. This is modeled, following Rzehak and Krepper (2013), from the turbulence velocity scale and the bubble length scale. The constant K BI is introduced to account for the modulation of the turbulence source, with its value optimized to 0.25 following comparison with a large database of bubbly flows (Colombo and Fairweather, 2015). The bubble-induced turbulence source in Equation (13) needs to be partitioned amongst the normal Reynolds stress components for use with the EB-RSM approach. In a previous work, we have obtained good agreement with experimental data by apportioning a higher fraction of the bubble-induced contribution to the streamwise direction, and the same assumption is used here (Colombo and Fairweather, 2015): In view of the low value of the density ratio in air-water bubbly flows, and the consequent much smaller values of the Reynolds stresses in the gas-phase, turbulence is resolved only in the continuous phase, while it is derived from this in the dispersed gas phase (Gosman et al., 1992;Behzadi et al., 2004). Values of the model coefficients used can be found in Table 1.

Numerical Settings
Numerical simulations were performed using the STAR-CCM+ code (CD-adapco, 2016). For pipe flows, a radial slice of the geometry was considered, whereas a quarter section of the domain was employed for the square duct. At the inlet, constant phase velocities and void fractions were imposed. An isotropic inlet turbulence profile with 2% turbulence intensity was also imposed at the inlet, to facilitate the development toward fullydeveloped conditions. At the outlet, a fixed pressure boundary condition, and zero gradient conditions on all the other variables, were employed. At the wall, no-slip boundary conditions were enforced on the velocity, with a zero value imposed for the turbulence stresses and the asymptotic limit ε = 2ν (k / y w ) yw→0 used for the turbulence energy dissipation rate. At the wall, the elliptic relaxation function α EB was set to zero, while a constant value equal to one was imposed at the inlet. On the lateral surfaces of the section, symmetry boundary conditions were imposed on all variables. Convective terms were discretized using second order upwind schemes and the pressure-velocity coupling was solved using a multiphase extension of the SIMPLE algorithm. Simulations were advanced implicitly in time using a second order scheme and, after an inlet development region, fully developed steady-state conditions were reached before recording the results. Strict convergence of residuals as ensured and the mass balance was checked to ensure an error always <0.1 % for both phases. Grid sensitivity studies ensured that solutions independent of the mesh resolution were reached. Given that significant resolution is needed to properly resolve the near-wall region, the first near-wall cell in the mesh employed was always located at a non-dimensional distance from the wall of y + ∼ 1. A sensitivity study is reported in Figure 1 for one of the experiments from Hosokawa and Tomiyama (2009) (experiment H22 in Table 2).
Specifically, the figure shows the radial profiles of the water mean velocity, the air void fraction, the radial turbulent stress (responsible for the radial pressure gradient and its impact on the bubble distribution) and the turbulent shear stress. Values are plotted as a function of the non-dimensional radial co-ordinate, which equals 0 in the center of the pipe and 1 at the pipe wall. Three different meshes were tested, having a number of mesh elements equal to 256 × 250 (on the surface perpendicular to the main flow motion and the axial direction, respectively), 420 × 400 and 975 × 750. Both velocity ( Figure 1A) and void fraction ( Figure 1B) do not show any appreciable difference even from the low to the medium refinement. For the radial normal stress and the turbulent shear stress in Figures 1C,D, some changes are visible from the low to the medium refinement, while no additional differences are found from the medium to the high refinement. Therefore, the mesh with medium refinement (168,000 elements) was used in the simulations. Comparable refinements were used for the other experiments. Given the fact that Hosokawa and Tomiyama (2009) employed a pipe of small diameter, a significantly higher number of elements was necessary for the larger diameter pipes, and 1,280,000 elements were employed for the square duct of Sun et al. (2014) (more details on the specific experiments are provided in the following section).

Experimental Database
The CFD model will be assessed against a database built using measurements from a number of literature datasets. These includes measurements in vertical upward pipe flows from the works of Liu (1998), Hosokawa and Tomiyama (2009), and Liu and Bankoff (1993a,b) in downward pipe flows from Kashinsky and Randin (1999), in large vertical pipes using the TOPFLOW facility (Shawkat et al., 2008;Lucas et al., 2010), and in a vertical square duct flows (Sun et al., 2014).
The data include a large variety of liquid and gas flow rates, void fractions, bubble diameters and geometries. Since the work described focused on the modeling of multiphase turbulence and interfacial momentum transfer, experiments having as much as possible a monodispersed bubble size distribution were selected as these can be effectively simulated using a constant value of the bubble diameter. Since the present database extends to cases at a high void fraction (for a bubbly flow), it is difficult to assume a monodispersed bubble size distribution in some instances. However, data that showed a marked wall-peaked (in upflow) and core-peaked (in downflow) void fraction profile were selected, allowing the use of a constant lift coefficient. Liquid and gas superficial velocities, and averaged values of the void fraction and the bubble diameter, were used to obtain the correct inlet conditions in the experiments. For the experiments of Liu (1998) and Shawkat et al. (2008), corrections to the nominal gas velocity and averaged void fractions were made from integration of measured radial profiles.
Hosokawa and Tomiyama (2009) measured liquid and gas velocities, void fraction, turbulence in the liquid phase, and bubble number and shape in a vertical upward pipe flow of  inside diameter 25 mm. Liu (1998) studied vertical upward flows of air bubbles in water in a pipe of inside diameter 57.2 mm, measuring the liquid velocity, turbulence and the void fraction. In the present simulations, the liquid superficial velocity and the averaged void fraction were imposed, and the latter was used to adjust the gas velocity to achieve the correct flux of air at the measurement position. Liu and Bankoff (1993a) measured liquid and gas velocities, turbulence levels, void fraction and bubble diameter in an upward air-water bubbly flow in a pipe of inside diameter 38 mm. Kashinsky and Randin   Table 2. Given that some of the experiments have been used multiple times in previous works [see, for example, Rzehak and Krepper (2015), Rzehak et al. (2017), Parekh and Rzehak (2018), and Magolan and Baglietto (2019)], to reduce confusion and favor consistency we have tried as much as possible to maintain the same naming convention.

Comparison With High-Reynolds Number RSM
In this section, using a first set of simulation results, the EB-RSM predictions are compared with existing results from a previously assessed high-Reynolds number model based on the SSG Reynolds stress model. Details of the model can be found in one of our previous publications (Colombo and Fairweather, 2015). Results are compared against experiments LB18 from Liu and Bankoff (1993a) (Figure 2), H22 from Hosokawa and Tomiyama (2009) (Figure 3) and L22A from Liu (1998) (Figure 4). On the ordinate, the plots show the void fraction, mean liquid (and gas where available) velocity and radial profiles of the turbulence kinetic energy, and the r.m.s. of the velocity fluctuations for Liu and Bankoff (1993a) and Hosokawa and Tomiyama (2009). This first comparison serves as an assessment of the capabilities of the EB-RSM compared with a more tested model. The models show similar capabilities and good accuracy. Void fraction profiles (Figures 2A, 3A, 4A) have the typical and expected wall-peaked features and are in very good agreement with the experiments. Most importantly, the magnitude and location of the wall peak are well-predicted by the EB-RSM without the addition of any wall force, and the relative uncertainties connected to a formulation that is, unavoidably, at least partially empirical. This is achieved because, as a consequence of the radial turbulent stress, a radial pressure gradient is generated in the liquid phase and bubbles are pushed toward the pressure minimum in the near-wall region. The radial behavior of the pressure for the three experiments considered is shown in Figure 5. Although the variation in pressure is only a few Pascal, it acts over a few millimeters and is sufficient to affect the bubble distribution and, in the near-wall region,  to produce the wall-peaked void fraction profile. Clearly, this pressure gradient, being related to the radial turbulent stress, can only be predicted with accuracy if the anisotropic structure of the turbulence field is accounted for through a Reynolds stress turbulence model, and the near-wall region properly resolved without resorting to wall functions.
Velocity profiles are also in good agreement with experiments and very similar predictions are found between the two models ( Figures 2B, 3B, 4B), although the EB-RSM provides a much more refined resolution in the near-wall region. Turbulence predictions are in agreement with data for Hosokawa and Tomiyama (2009) (Figure 3C) and Liu (1998) (Figure 4C), although both models show discrepancies for Liu and Bankoff (1993a) (Figures 2C,D), for which turbulence levels and the Reynolds shear stress are underpredicted. Additional discussion of these discrepancies will be provided later when predictions   with additional data from the same database are presented. It is worth noting here that, while specifically the streamwise, but also the radial, r.m.s. velocities are underpredicted in experiment LB18, much better agreement is found in all three co-ordinate directions for experiment H22 (Figure 3D) from Hosokawa and Tomiyama (2009), proving the models' ability to predict the anisotropic turbulence structure. Away from the wall, predictions from the two models are similar. Near the wall, in contrast, the EB-RSM's higher resolution provides a much improved prediction of the turbulence kinetic energy peak, at least for the cases where turbulence measurements near the wall are available (Figures 3C,D, 4C). In contrast, the SSG model employs a high-Reynolds formulation and wall functions near the wall. For this reason, the first near-wall cell is located at a nondimensional distance from the wall of y + ∼ 30, necessary for the wall function to be valid, and the formulation is not able, as expected, to predict the near-wall peak. Overall, the superior capabilities of the EB-RSM in the near-wall region, and their impact on model predictions, in particular, for the void fraction and turbulence in the fluid phase, suggests that models with nearwall modeling capabilities should be adopted wherever possible. For this reason, only predictions from the EB-RSM will be addressed in the following.

Upward Pipe Flows
In this section, results from the comparisons against the remaining upward pipe flows (excluding large diameter pipes that will be addressed separately) in the database are discussed, and these include the cases from Tomiyama (2009), Liu andBankoff (1993a) and Liu (1998). Void fraction profiles were almost always found to be in good agreement with experiments, and Figure 6 provides some examples. Discrepancies are only found in experiment H12 (Figure 6A) in the center of the pipe, and these can be attributed to the presence of some larger bubbles that cannot be predicted using a constant positive value of the lift coefficient. From Hosokawa and Tomiyama (2009), measurements show how a tail of bubbles in the range 6-8 mm is present in the bubble size distribution of this experiment only. The predictions remain accurate at large void fractions, approaching or even exceeding 0.4, such as in experiments LB3 and LB21 (Figure 6B), although for larger liquid velocities (LB33, Figure 6C), the model has a tendency to overpredict the void fraction peak. Overall, the accuracy of the model, even without the addition of any wall force contribution, is encouraging. Clearly, improvements to the lift model, allowing it to predict the change in sign of the lift coefficient, are necessary to extend the overall model's applicability to polydispersed bubbly flows.
Velocity profiles are in general in good agreement with data, and well-predicted in the very near-wall region (Figure 7A), for cases where the void fraction is not excessively high. At high void fraction, the velocity, driven by bubble buoyancy, tends to peak at the wall (Figures 7B,C), and specific modifications to the drag and near-wall turbulence models might be necessary. Overall, the velocity profiles are almost flat and the slight decrease toward the center of the pipe causes the discrepancies in the prediction of the turbulent shear stress mentioned previously ( Figure 2D). It is worth noting that in Liu and Bankoff (1993a) measurements were taken at the shortest distance from the inlet and the turbulent shear stress behavior is not entirely consistent with the close-to-flat velocity profile also observed in the experiments (Figures 2C,D). Additional model testing is therefore desirable, although detailed experimental characterization of the turbulence structure in pipe flows is not frequently found in the literature. The relative velocity, mainly governed by the drag model, tends to be underpredicted on occasion (Figure 7C), contributing to the peak in the liquid velocity profile.
Overall, the area in need of most improvement is multiphase turbulence modeling, and the model for bubble-induced turbulence specifically. Although accurate in many cases (Figures 8A, 3C, 4C) discrepancies are found, with the data from Liu and Bankoff (1993a) generally underpredicted ( Figure 8C). It is therefore necessary to progress beyond the simple constant values used in the bubble-induced source and the partitioning of this between the different stresses [Equations (13,15)]. Nevertheless, available measurements confirm how the nearwall resolution allows accurate prediction of the peak in the turbulence kinetic energy at the wall (Figures 8A,B).

Downward Pipe Flows
This section extends the previous comparisons, limited to upflow cases, to downward flow conditions. In these flows, the bubbles travel at a lower velocity with respect to the liquid phase, and the same lift force pushes them toward the higher relative velocity region that is now in the center of the pipe. This produces core-peaked void fraction profiles, such as those shown in Figures 9A,D for the experiments of Kashinsky and Randin (1999). This behavior is predicted by the model which maintains consistent predictions of the void fraction and shows very good agreement in the near-wall region.
Velocity profiles are in very good agreement (Figures 9B,E), and predictions of the mean velocity and turbulence (Figures 9C,F) are remarkably good in the near-wall region. For these experiments only, the r.m.s. of the streamwise velocity fluctuations is shown as a function of the non-dimensional distance from the wall, in agreement with the way they were originally provided by Kashinsky and Randin (1999). The turbulence intensity is underpredicted by 20-30 % in the center of the pipe. However, in these experiments an additional complication is included since the bubbles have a much smaller diameter than in previous cases, for which the contribution to turbulence from their wakes may become negligible with respect to that due to their random motion. The latter is not properly captured by the type of model used for bubbleinduced turbulence in the present work, it being based on the conversion of energy from drag to turbulence kinetic energy in the bubble wakes.

Vertical Pipes of Large Diameter
In previous comparisons, measurements were taken in pipes with diameters of a few centimeters. In the database, two cases with much larger diameter pipes are also included; TL12-041 from the TOPFLOW facility  and S1 from Shawkat et al. (2008). These are included to extend the model assessment as much as possible. The pipes considered have diameters of tens of centimeters and their hydrodynamics can be considered to have features similar to those of bubble columns. Wall effects are still present, but their impact extends into the body of the flow much less than in smaller diameter pipes, with uniform velocity and void fraction distributions dominated by mixing found in the bulk of the pipe cross-section. This is clearly evident in the results of Figure 10, where wall effects are confined to the very near-wall region. The model maintains satisfactory agreement with data in these conditions, with the overprediction of the void fraction peak at the wall being the major discrepancy. The air velocity, which has shown some discrepancies in a few previous cases, is instead well-predicted in TL12-041 ( Figure 10B). It must also be pointed out that measurements of air velocity and turbulence levels are also available for case S1. However, large and not entirely explicable discrepancies with model predictions have been observed in other studies, for example Parekh and Rzehak (2018), and, for this reason, these data are not included here.

Square Duct
Finally, predictions are compared against data in a square channel having a relatively large cross sectional area (Sun et al., 2014). In Figures 11, 12, results are shown for three distinct locations, along the diagonal D of the duct, with the abscissa showing the non-dimensional distance from the duct center d / D, and along two lines parallel to the duct wall, one passing through the center of the duct and the other in the near-wall region. For the latter two locations, the abscissa identifies the non-dimensional distance along the line from the center of the duct, x / L. Void fraction behavior has similarities with pipe flows,    Figures 11C,F, it may be noted that the latter are higher, given that they essentially detect the near-wall peak along the wall. The highest values are found in the duct corners, as shown by the peaks along the diagonal in Figures 11A,D, and the line parallel to the near-wall region in Figures 11C,F. The void fraction is well-reproduced by the model for both experiments considered in Figure 11, particularly along the diagonal and the line through the duct center. Some underestimation is found on the line in the near-wall region, although the peak in the duct corner remains well-predicted.
Lastly, Figure 12 shows mean velocity profiles for experiment SN1 and turbulence kinetic energy, normalized with the square of the mean velocity, for experiment SN2. The qualitative behavior of the velocity, peaking at the walls and corners and slightly decreasing toward the duct center, are well-captured, although the wall-peak is underpredicted. Turbulence levels, in line with previous results, are underestimated in the center of the duct, again confirming remaining challenges in obtaining accurate predictions. In this case, the cross-sectional area is large and the void fraction relatively high, and additional effects promoted by the interaction between bubbles may add to the contribution from their wakes, explaining the underestimation. Peaks in the turbulence kinetic energy are visible at the walls, although the absence of measurements taken in the very near-wall region prevents any comprehensive assessment. Overall, additional model improvements are necessary, but the model performance remains robust even in the square duct geometry.

CONCLUSIONS
The performance of an Eulerian-Eulerian multi-fluid CFD model coupled with the EB-RSM turbulence closure, allowing fine resolution of the near-wall region, has been assessed against a large database of air-water bubbly flows. The overall model has demonstrated robust applicability and encouraging accuracy over a wide range of operating conditions and geometrical parameters, with the experimental database including upward and downward pipe flows, large diameter pipe and square duct flows, some of which have rarely, if ever, been tested against.
Good agreement with data for the void fraction distribution was obtained for all experimental cases, with accurate predictions of the magnitude and positon of the near-wall peak found in most experiments. The model does not account for any wall lubrication force, and all the associated uncertainties in modeling this force, and the wall peak is obtained from the action of the radial pressure gradient induced by the turbulence field and its accurate prediction near the wall. In this regard, these results confirm and support recent findings where wall peaked void fraction distributions were also obtained without any wall lubrication effect, as a consequence of which the theoretical basis of the wall lubrication force was questioned. Further evidence will be required, as well as extension to the modeling to handle laminar conditions where turbulent effects are not present, such as with the addition of the turbulent dispersion regularization model of Lubchenko et al. (2018). Overall, and although the lift force remains dominant, these mentioned additional turbulent effects should be taken into account as accurately as possible in any model that aims to be general applicability. For this reason, the transition to models based on Reynolds stress turbulence model formulations, able to capture the anisotropy of the turbulence field, with near-wall resolution capabilities is to be preferred when possible. Such approaches allow fine resolution of the velocity and turbulence field near the wall, although modifications seem necessary at high void fraction given that the elliptic-blending model is still based on its single-phase counterpart. For this to be achieved, the availability of detailed measurements of the turbulent stresses in the near-wall region for further model validation is desirable.
Overall, multiphase turbulence models and the modeling of bubble-induced turbulence are areas in need of major improvement. Over the wide range of conditions tested, multiple physical processes coexist which all contribute in determining a flows' characteristics, which makes predictions particularly challenging. In the short term, more physically based modeling improvements of the bubble-induced source and the attribution of a larger portion of that source to the streamwise direction are necessary. In this regard, some advances have started to appear in the very recent literature (du Cluzeau et al., 2019;Liao et al., 2019).
Finally, it is worth stressing the role of assessments of the kind described for the future development of computational models of bubbly flows. When focusing on only a few experiments for model validation purposes, these can show trends not observable in other cases (see, for example, the differences between Figures 7C, 10B for the air velocity, and Figures 2C, 3D for the turbulence field), and results and modeling improvements should be accepted only when consistent over multiple datasets.

DATA AVAILABILITY STATEMENT
Datasets obtained from the CFD simulations with the EB-RSM for the experiments included in this paper are available at https:// doi.org/10.5518/800.

AUTHOR CONTRIBUTIONS
MC and MF conceived and planned the research, wrote and finalized the manuscript. MC developed the model and performed the simulations.