A Comparative Study of Closure Relations for CFD Modelling of Bubbly Flow in a Vertical Pipe

Commercial code CFX was used to examine the performance of a two-fluid model to predict the details of upward isothermal bubbly flow of air and water in a vertical pipe. The model equations are volume-averaged Navier-Stokes equations that require closure models for interfacial forces and bubble-induced turbulence effects. Two-equation SST and k-epsilon RANS turbulence models were also used. A parametric study of closure models included both standard options in CFX and previously published novel closure models that were implemented with user-defined functions. The CFD simulations were compared with two cases from the MTLoop experiments by Lucas et al. at the Helmholtz-Zentrum Dresden Rossendorf: one with wall-peak void fraction profile (MT039), and another with a core-peak void fraction profile (MT118). The effect of changing the drag force closures was not significant for the set examined. Poor predictions were found when the lift force and wall lubrication models were incompatible in magnitude. There was no significant effect of changing the liquid phase turbulence model. Changing the bubble-induced turbulence models, however, had a significant impact on the radial void fraction profile. The novel wall force from Lubchenko et al. at the Massachusetts Institute of Technology significantly improved the prediction of the near wall void fraction in the wall peak profile.


Introduction
Two-phase gas-liquid flows occur in many important industrial processes. It is challenging to gain a strong understanding of such flows because of the complex physical interactions between the phases that often occur. Many different com-plex flow structures can be found, depending on the flow geometry, orientation, fluid properties, and the mass flow rates of gas and liquid. The capability to predict the two-phase flow behaviour would be a valuable tool in the design of important industrial equipment, so the development of suitable computational fluid dynamics (CFD) models is an active area of ongoing research.
In the area of isothermal bubbly flow (adiabatic flow of a dispersed gas in a continuous liquid phase), a common modelling approach for practical applications is the two-fluid (Euler-Euler) framework [1]. In this framework, the governing equations are volume-averaged, which makes CFD simulations on the large length scales possible, but the details of the small-scale interfacial interactions must then be supplied by auxiliary closure relations. The key interfacial forces that are provided in the closure relations are drag, lift, wall lubrication, turbulence dispersion, and virtual mass. Another key effect that is modelled is the effect of the bubbles on the turbulent flow of the liquid. Reynolds-averaged Navier-Stokes (RANS) models are most often used to model turbulence quantities in the continuous phase. In addition, bubbly flows are generally polydisperse (have a variety of bubble sizes), so there are complex bubble breakup and coalescence processes occurring that influence the bubble size distribution.
A significant number of studies using numerical models of bubbly flow in vertical pipes have also been carried out. Khan et al. [22] provide an excellent over- There are an enormous number of possible combinations of these relations and no consensus has been reached on the set to select. Chuang and Hibiki [23] also discuss interfacial force in gas-liquid CFD modelling, with a focus on a bubble-wall collision force. Lucas et al. [24] [25] and Liao et al. [26] proposed approaches for development of improved closure models. Because there has been no consensus on the set of closure relations to use for a two-fluid model for bubbly flow in a vertical pipe, many earlier works explored the effects of using different closure relations and validation of codes using the experimental data previously mentioned. Examples of studies using polydisperse multi-dimensional two-fluid models are in [27]- [35]. A significant number of studies using monodisperse multi-dimensional two-fluid models have originated from the group at HZDR (Helmholtz-Zentrum Dresden-Rossendorf); some examples are in [36]- [41]. Examples of other studies in the same area are those of Yamoah et al. [42], Wang and Yao [43], Lote et al. [44], and Colombo et al. [45].
All of these studies examined the effects of the choice of closure relations on the predictions of the model.
The present work is restricted to monodisperse two-fluid modelling of isothermal bubbly flow in a vertical pipe using commercial CFD code ANSYS CFX. The focus of the work is comparing conventional built-in closure models with newer closure relations that have not been examined elsewhere. Those newer closure relations were implemented with user-defined functions. The performance of the closure set combinations is examined by comparison with two experimental data sets from the HZDR MTLOOP experiments [15]. The comparison is made for data in the fully developed region. One data set has a wall-peak void fraction distribution (MT039) and the other has a core-peak void fraction distribution (MT118).

Problem Description
Vertical flow in a pipe is modelled assuming axisymmetry. Therefore, similar to the approach used by Rzehak [39], a 5-degree wedge geometry was used. The computational domain is a pipe with length, L and radius R, as shown schematically in Figure 1.
Because the focus of the present study is limited to fully developed flow, it is not necessary to resolve the air injection apparatus. Instead, the flow enters the bottom of the domain with uniform velocity and mass flow rates of air and water specified to be consistent with the fully developed condition. Symmetry boundary conditions are used for the sides of the wedge. The pipe wall is represented

Two-Fluid Model
Incompressible flow of air and water is modelled using the two-fluid framework [1]. In the theory manual for CFX, this is referred to as the inhomogeneous two-phase model approach.

Continuity Equations
The mass conservation equations for both phases are: The volume fraction constraint enforces 1 G L α α + = in each control volume.
Both phases share the same pressure field.

Momentum Equations
The momentum equations for both phases are: The inter-phase momentum transfer consists of: Because the flow is steady and the axial development is not considered, the virtual mass force is neglected in the present study.
The drag force is calculated as: where R U is the relative slip velocity, defined as: The lateral lift force is expressed as [46]: The wall lubrication force is generally expressed as [47]: where wall n is the wall normal unit vector, and || R U is the component of relative velocity orthogonal to wall n .
The turbulent dispersion force is calculated using the Favre-averaged drag (FAD) model [48]: Open Journal of Fluid Dynamics , 3 1 where 0.9 TD σ = , except when using the Magolan and Baglietto bubble-induced turbulence model [49], as described in Section 2.6. Coefficients D C , L C , and WL C are given with the details of the interfacial force closure models in Section 2.4.
The stress tensor, T , includes viscous and turbulent stresses: where the effective viscosities are eff ,

Turbulence Model
Although the SST turbulence model has been the predominant choice amongst previous numerical works, both the SST and k-ε turbulence models were included in this study. The two-equation turbulence models are applied to the continuous liquid phase. Due to the large density difference and small void fraction, the turbulence within the dispersed gas phase is not significant. The effect of the dispersed phase on the turbulence field of the continuous phase, however, is significant. This phenomenon, commonly referred to as bubble-induced turbulence (BIT) is accounted for by additional terms in the liquid phase turbulence model, as described in Section 2.6. The governing equations for the turbulence models documented in the CFX theory manual [50] are as follows:  k-ε (keps): , , where L S is the strain rate tensor of the liquid phase defined as: The blending function 1 where 10 min 1 10 ξ − = × and wall y is the distance to the nearest wall. The turbulent viscosity is then calculated using: where 2 2 2 wall wall 2 500 tanh max , Standard turbulence model constants were used, as listed in Table 2. In the SST turbulence model, the constants are interpolated between the values for the k-ε and k-ω models using the blending function:  The turbulence in the dispersed phase is approximated from ,L t µ using a zero equation model: where t σ is the ratio of liquid and gas eddy viscosities; the default value of 1 was used in this work.

Interfacial Force Closures
The primary goal of the present study is to examine the behaviour and performance of combinations of closure relationships. Two types of closure models are examined: 1) standard closure relations that were previously applied to the experiments of interest and 2) novel closure models that have not been tested on the dataset selected. Only closure models that are applicable for both wall peak and core peak profiles were considered. Furthermore, closure models found by other authors to give poor predictions of the flow were not considered.

Drag Force
The drag force acts in line with the direction of flow and is typically the largest of the interfacial forces. It represents the resistance to relative motion the continuous fluid imposes on the dispersed phase (bubble or particle). As the mean bubble diameter increases, the bubbles deform, and the analytically derived expression for the drag force on rigid spheres is no longer valid. Consequently, drag models such as Schiller and Naumann [51] are not used in the present work. Many drag models accounting for bubble deformation are available and most have little effective difference on C D as demonstrated by Wang and Yao [43].
The Ishii-Zuber [52] (IZ) drag model was selected as the standard because it is included by default in CFX. Although the Tomiyama et al. [53] (TKZS) drag model is not in CFX and is very similar to the IZ drag model, it was included because it is the basis of the novel drag model by Buffo et al. [54] (BVRD). The BVRD model was selected because its bubble swarm modifier has an additional functional dependency on local turbulence. The bubble swarm modifier accounts for the interactions between closely packed bubbles (a swarm). The swarm correction is important in situations where the flow conditions differ greatly from the low bubble concentration studies from which drag models have been derived [54]. IZ where crit α is the maximum packing; the default value of 1 was used.

Lateral Lift Force
The lateral lift force acts perpendicular to the streamwise direction of flow. It arises from the rotation and deformation of bubbles due to velocity gradient in the continuous phase. For relatively small bubbles, the lift force pushes the gas phase towards the wall, thereby creating the wall peak void profile. Larger bubbles experience greater deformation which causes the direction of lift force to reverse, concentrating the void fraction at the centre of the pipe, thereby creating the core peak profile [55].
Lift force reversal is the primary factor to consider when selecting a broadly applicable lift model. For this reason, common lift models such as Legendre-Magnaudet [56], and Saffman-Mei-Klausner [57] were not considered. The where: and The lateral bubble diameter, d ⊥ , is the maximum horizontal dimension of the bubble measured along the wall-normal direction. It serves as an indicator of bubble deformation. An empirically derived expression by Wellek et al. [60] is used.
Eo ⊥ is defined by Equation (44) and d ⊥ is calculated using Equation (45) with 0.7 A = and 0.7 B = .
An additional lift force closure, the novel Shaver and Podowski (SP) model, was also implemented. It is discussed with its companion wall force model in Section 2.4.4.

Wall Lubrication Force
The wall lubrication force acts normal to the wall, repelling the bubbles and thus moving the void fraction peak away from the wall. In the conventional lift-wall modelling approach, the lift force reaches a maximum at the wall because the velocity gradient increases near the wall. The wall lubrication force counteracts the lift force in the near-wall region.
The shortcoming of many wall lubrication models is the need for tuning parameters. For example, the commonly used wall lubrication model by Antal et al. [47] was excluded because it is very sensitive to its fitting coefficients and a wide range of values have been cited by various authors.
In the set of wall models selected, the Tomiyama et al. [61] (T) and Frank et al. [28] (FZKPL) wall models are built-in to CFX, whereas the Rzehak et al. [36] (RKL) model was implemented with a user-defined function.
Recent studies have questioned the physical mechanisms of the conventional lift-wall modelling approach [62]. A novel lift-wall modelling approach based on the wall force model by Lubchenko et al. [62] (LMSB) will be discussed in Section 2.4.4.
T Tomiyama et al. [61] ( ) where W C is calculated using Equation (48). The default model constants are RKL Rzehak et al. [36] Open Journal of Fluid Dynamics

Novel Lift-Wall Model
The Shaver and Podowski [63] (SP) approach is not a conventional lift model that predicts a lift coefficient based on a correlation. Rather, it is a near-wall modification of the lift coefficient. Specifically, the lift coefficient is decreased to zero near the wall to prevent the non-physically large lift force at the wall previously mentioned, thus eliminating the need for the conventional wall lubrication force.
Lubchenko et al. [62] (LMSB) proposed a fundamentally different formulation for the wall force, treating it as a renormalization of turbulent dispersion [62].
The use of the LMSB wall lubrication force requires the SP lift force modification.
The combination of the Shaver and Podowski lift model and the Lubchenko et al.
wall model was implemented in CFX using user-defined functions. The novel lift-wall approach has an added advantage of being more numerically stable when using a fine mesh at the wall because the lift force is reduced in wall-adjacent cells where the velocity gradient is very large.

SP Shaver and Podowski [63]
is a nominal lift coefficient. In the present study the lift coefficient predicted by the ZTL correlation (Equation (46)) was used.

Bubble-Induced Turbulence (BIT) Closure
Based on the findings of Rzehak and Krepper [64], the Sato and Sekoguchi [65] (SS) and Rzehak-Krepper [38] (RK) BIT models were selected as commonly used models for the present study. The novel BIT models selected are the Ma et al. [66] (MSZLF) and Magolan and Baglietto [49] (MB) models that were assembled from recent DNS studies. The SS eddy viscosity modifier is the only BIT model built-in to CFX.
The source terms in the L k , and L ε or L ω equations account for bubble-induced contributions to the primary phase turbulence, in accordance with the BIT model used. The source terms share the same general form between all the BIT models used here except for the one by Sato and Sekoguchi, as discussed shortly. In general, the source terms are given by the following expressions: (53) , , The relation presented for ,L S ω (Equation (55)) is simply a conversion of the ε source (Equation (54)) into an equivalent ω source for use in the SST model. The specifics of the BIT models are summarized in Table 4.
where ,L t µ is the standard shear-induced eddy viscosity defined by (Equation (17) or Equation (22)) and BI t µ is calculated as: The effective viscosity becomes eff , , .
The turbulent dispersion force was modified when the MB BIT model was used. When developing their BIT model, Magolan and Baglietto decoupled the feedback between interfacial forces and turbulence by replacing the eddy viscosity in the turbulent dispersion force with the Sato eddy viscosity modifier (Equation (57)). When re-coupling the closure relations with the new BIT model, it was necessary to re-calibrate the turbulent dispersion force. This was accomplished by adjusting TD σ to produce more realistic results. The values of TD σ used in this study are given in Section 3.

Boundary Conditions
At the bottom of the domain, a per-phase mass flow rate inlet boundary condition was applied. The mass flow rates of air and water were calculated based on the inlet area and superficial velocities listed in Table 5 and were uniformly distributed over the inlet. The average void fraction at the inlet was specified to the experimentally measured value at the fully-developed location in Table 5. The velocity of each phase was calculated internally by CFX based on the density and volume fraction. At the inlet, the turbulence intensity was set to 5% and the ratio of eddy viscosity to molecular dynamic viscosity was 10. Symmetry boundary conditions were used for the sides of the wedge. The arc-shaped pipe wall was prescribed by a no-slip and impermeable boundary for both phases. An average relative pressure of zero was specific over the outlet area at the top of the pipe.
For the SST model, the automatic wall function consistent with the k-ω model was used. For the k-ε model, the scalable wall function was used.

Definition of Test Cases
The MTLoop dataset by Lucas et al. [15] was selected for comparisons in the present work. In the experiments, air and water entered with controlled mass  [67] was used to measure the spatial variation of void fraction and bubble size at a given cross-section. The accuracy of the WMS and its effect on the flow were studied by Prasser et al. [68] by comparing against the established benchmark of X-ray tomography. The agreement between the WMS measurements and the benchmark was within 4% for the range of bubbly and slug flows tested. Because the wire-mesh sensor is an invasive technique, each axial location was measured one at a time in separate trials. From the approximately 100 combinations of air and water flow rates recorded in the MTLoop dataset, a configuration with a wall peak void profile (MT039) and a configuration with a core peak void profile (MT118) were selected to give a representative sample of the range of physics associated with bubbly flow. The details of the selected test conditions are given in Table 5.
An extra length of approximately 10% was added to the length of the computational domain to ensure the outlet boundary condition did not affect the results at z m . Therefore, for this work the length, L, of the computational domain shown in Figure 1  The assumption of constant air density is consistent with previous work (e.g., [28] [39] [44]). The experimentally measured mean bubble sizes reported in Table 5 correspond to the 60 L D ≈ measurement location. The measured bubble diameters are used to specify the bubble diameter in the simulations in the present work.
Exhaustively testing every possible combinations of the selected closure relations would require 288 cases. Not only would that method be unnecessarily time-consuming and inefficient, there are practical constraints on the combination of interfacial forces that could be used. For example, the MB BIT model is only used with the k-ε turbulence model in this work because that was the primary-phase turbulence model used in its development and adding the Recognizing that some parameters are more strongly connected than others, a more feasible set of cases can be assembled. For example, the lift model will have the same relative effect regardless of the primary-phase turbulence model used. The BIT model, however, is likely to have a different interaction with different primary-phase turbulence models. Table 6 summarizes the closure relations used in the parametric study cases in this work. Solutions for cases 08, 10 and 24, were not obtained for MT118 because of solution divergence. It is thought that the divergence is a result of the  T  T  RK  SST   12  BVRD  T  T  MSZLF  SST   13  BVRD  T  T  MB  keps   14  IZ  T  T  RK  SST   15  IZ  T  T  MSZLF  SST   16  IZ  T  T  SS  keps   17  IZ  T  T  RK  keps   18  IZ  T  T  MSZLF  keps   19  IZ  T  T

Methodology
Commercial CFD code CFX (Version 2020-R1) was used to solve the governing equations presented in Section 2.2. In CFX, an element-based finite volume approach is used. Standard finite element gradient approximations were used for diffusion terms and the high resolution scheme was used for advection terms for the momentum equations. A first order upwind scheme was used for the turbu- The effect of convergence criteria on the solution accuracy was studied for a representative case (Case 01 for the MT039 data set discussed earlier). As the maximum residual target was progressively decreased from 1 × 10 −3 to 1 × 10 −6 , negligible difference was observed in the results. Radial void profiles had a 0.04% root mean square (RMS) difference between 1 × 10 −3 and 1 × 10 −6 maximum residual targets. The solution (run) time for the initial case was relatively low (under 1 hour using 4 Intel 2.40 GHz i7 cores), and the run time of subsequent cases was reduced by using a similar prior case for initialization. Because of the relatively low computational effort needed, a maximum residual target of 1 × 10 −5 was used unless otherwise noted. A fixed time step of 0.01 s was used for the majority of cases. In some instances with multiple user-defined source terms or for the larger bubble diameter of MT118 cases, it was necessary to reduce the time step size to 0.001 s in order to achieve convergence.
Because of the small time step size needed for numerical stability, monitoring global conservation is important to ensure a steady-state solution is reached.
Special attention must be paid to the air mass imbalance. By default, all phase mass flows are normalized by the largest mass flow of any phase. Because the mass flow of water is much larger than that for air, the percentage domain imbalance of air reported by CFX can be misleading in this type of problem. To address this, a user-defined variable equal to the air mass imbalance normalized by the air inlet mass flow rate was monitored. In addition to the maximum resi-Open Journal of Fluid Dynamics dual criterion, the maximum allowable domain imbalances were 0.001%. The average volume fraction at the outlet was also monitored to confirm that a steady-state solution had been obtained.

Mesh Sensitivity Study
When using the multiphase particle model, the inter-phase momentum transfer models limit how fine the grid can be made before there are numerical stability issues [39]. As a general rule, the larger the particle diameter, the larger the grid must be for numerical stability. A compromise was reached between a finer grid to better resolve the velocity gradients and turbulence quantities, and a coarser grid to improve numerical stability and convergence behaviour.
A mesh sensitivity study was performed for Case 01 for the MT039 data set using with three grids listed in Table 7. The RMS difference between the void fraction profiles at z m for a mesh and the next finest mesh were compared. The finest grid (M3) required reducing the time step size by a factor of ten to prevent non-convergence because of oscillating residuals. In order to ensure good convergence behaviour with an acceptable change relative to the next finer grid, M1 was selected as the working mesh. This mesh has similar resolution to that used by Rzehak [39] and is shown schematically in Figure 2. The size of the first cell at the wall (Δ 1 ) influences the turbulence model through the y + , as well as the wall lubrication force. The first grid spacing near the wall (and its corresponding y + value of 16) for M1 is within the range commonly used for this problem in previous studies [39] [43] [44] [69]. The first node distance is within the active range for all wall lubrication forces considered.

User-Defined Closure Models
User-defined models were implemented in CFX using CEL command language.
Non-standard drag and lift models were implemented by specifying a user-defined expression for C D and C L , respectively. It is not as straightforward to implement a custom wall force, as it is not possible to specify a user-defined expression for C WL . Furthermore, in the case of the LMSB model, a completely different expression for wall force is required. Consequently, the RKL and LMSB models were implemented in CFX as user-defined momentum sources. The expression for y wall was limited to a minimum value of 1 × 10 −4 to prevent the 1/y wall term in the wall force from becoming too large and causing convergence problems with the user-defined momentum source. The BIT models were implemented with user-defined sources in the appropriate turbulence equations. To add the B C µ coefficient for the MB BIT model a user-defined expression was supplied for the eddy viscosity. To adjust the turbulent dispersion force, CFX gives the option of a user-specified turbulent dispersion coefficient (in the numerator of Equation (10)) but TD σ cannot be adjusted directly. Instead, the turbulent dispersion coefficient was specified to obtain the desired value of TD σ .
In the implementation the BVRD drag model, a CEL expression could not be given for C D because CFX could not evaluate the local turbulence quantities in that expression. As a work around, a nominal drag coefficient was specified and then a user-defined momentum source was created to subtract the nominal drag force and add the drag force calculated from the BVRD model. It was not possible to verify all user-defined closures by comparisons with previous numerical results. It was possible, however, to verify the combination of the RKL wall lubrication model and the RK BIT model by comparison with the work of Rzehak and Kriebitzsch [39]. The closure set of Case 24 (IZ/T/RKL/RK/SST/FAD for drag/lift/wall/BIT/turbulence/dispersion) was used in the present study to compare with the results for gas void fraction, gas axial velocity and eddy viscosity with those for MT039 from Rzehak and Kriebitzsch [39], who also used CFX. Figure 3 shows excellent agreement between the two results, which demonstrates that the user-defined source terms correctly implemented the RKL wall lubrication and RK BIT models.

Results and Discussion
The effects of the closure models on the CFD solution are mainly examined in terms of agreement of radial profiles of void fraction and gas axial velocity with the MTLoop data from the fully developed measurement location ( 3.03 m m z = ). An effort was made to limit the focus to a single parameter at a time, but this is not always possible due to the highly non-linear interaction between equations.

Effect of Drag and Conventional Lift-Wall Closures
Considering the effect of drag force, there is little difference between the three drag models, as seen in Figure 4 and To better understand the discrepancy between the ZTL and T lift models in Figure 4, a deeper examination of the lift coefficient is required. When the assumption of monodispersed and constant density is made, Eo does not vary spatially. Therefore, the lift coefficient is constant. The T and ZTL lift models    The FZKPL and T wall lubrication models share the same C W (Equation (48)) but have different means of determining the wall force C WL . The T model uses the pipe diameter in determining the wall force, limiting its application to circular pipes. The FZKPL model removes the dependence on a pipe diameter, but introduces cut-off coefficients. The FZKPL model is more universal in terms of geometry, but has a less physical basis. The RK wall model was derived from the data of Hosokawa et al. [70] and shares the same correlated functional dependency on Eo but the effective distance from the wall differs. The RK wall model decreases at a rate of 2 wall 1 y with no dependence on pipe diameter or cutoff coefficients, whereas the original Hosokawa et al. work follows the approach of either Tomiyama et al. [61] or Antal et al. [47]. The FZKPL and RKL wall models are compatible with the magnitude of the T lift force, and give reasonable results.

Effect of Drag and Turbulence Closures
The BIT model influences the radial void fraction profile primarily through the turbulent dispersion force which is proportional to t µ . Considering Case 01 for MT039, the SS eddy viscosity modifier produces a large value of t µ relative to other BIT models, leading to the much flatter void profile seen in Figure 8. In the case of MT118, the t µ predicted by the RK BIT model is too small, which decreases the smoothing of the void fraction profile and permits the non-physical inflections seen in Figure 9. It is important to remember the   turbulent dispersion force was re-calibrated for the MB BIT model by adjusting TD σ as suggested in [49].
Some of the source term style BIT models do not predict well the void fraction distribution for the MT118 case, as seen in Figure 9. It is important to recall that a monodispersed assumption is used to model the flow but the MT118 flow has polydispersity. The mechanisms of BIT and bubble size distribution are highly intertwined [33], notably because d b appears in the BIT model equations. The Open Journal of Fluid Dynamics RK and MB BIT models exhibit physically inconsistent inflection of the radial void profile but they better predict the mean void fraction, as discussed later in Section 5.4.
The k-ε model has been reported by a number of authors to overpredict the near wall void fraction, relative to the SST model [42]. This trend is observed in Figure 10 with the use of the Sato eddy viscosity modifier. Interestingly, there is little difference between k-ε and SST turbulence models when the source term BIT models are used. The BIT source terms added to the k and ε equations are larger in magnitude than the difference between the unmodified k-ε and SST formulations, as demonstrated in Figure 11 and Figure 12. The difference in ,L t µ between Cases 01 and 16 is far less in MT118 than MT039, but the effect on the radial void profile is the opposite.
In Figure 13, Cases 01 and 16 indicate that the void fraction predicted by the k-ε model is lower at the wall and higher in the core than for the SST model. The    The SS BIT model performed well overall, but the approach has disadvantages for modelling additional physics that depend on local k or ε. A common example is the modelling of breakup and coalescence rate for polydispersed simulations [33]. Of the source term style BIT models considered, the MSZLF and MB BIT models performed reasonably well in terms of void fraction and velocity for both MT039 and MT118, yet predict very different turbulent quantities. The eddy viscosity predicted by the MB model is larger than by the MSZLF model, whereas the L k and L ε values are similar between the two. Generally, the MB BIT model provided better prediction of the velocity profile than other BIT models. Without experimental measurements of turbulence quantities it is uncertain which BIT model is more accurate.

Effect of Novel Lift-Wall Approach
The difference between the conventional and LMSB lift-wall approach is hig- duced by the SP modification, the LMSB wall force must be much smaller than the conventional wall lubrication force to yield a physically reasonable void profile. Not only is the LMSB wall force much smaller in magnitude than the conventional approach, it is also active over a much shorter distance. Figure 14 clearly demonstrates the disproportionate magnitudes of the lateral redistribution forces. The T lift force from Case 01 is approximately eight times larger in magnitude than Case 04 at their respective peaks.   The LMSB wall force is proportional to t µ so it is necessary to re-examine the effects of BIT models. As expected, the MB BIT model performs better with LMSB wall model (for which it was derived) than conventional lift-wall models. The same trends of radial void fraction profiles are observed between BIT models, regardless of the wall model used. The various BIT models have a somewhat lesser effect on the radial void profile when the LMSB lift-wall approach is used (Figure 15), than when a conventional lift-wall approach is used (Figure 8). It is noted that the peaks of the void fraction profile are closer together in Figure  15(a) than in Figure 8(a).
The radial profiles of turbulence quantities vary slightly with the change in lift-wall approach but the general shape of the curves remains the same. The differences in magnitude can be explained by the difference in void fraction distribution. There does not appear to be a high degree of feedback between the turbulence models and the LMSB wall force. The mean void fraction and gas velocity were not affected significantly by the change in the lift wall approach.
Although the purpose of the present work is not to optimize the agreement of a closure set for a particular flow condition, it is not unreasonable to adjust the TD σ for Case 26. Magolan [71] suggested a variable TD σ correlated to void fraction and local turbulence for use with the MB BIT model. The complete formulation by Magolan [71] was not used in the present work due to technical limitations with CFX. Instead, TD σ was simply decreased from 3 to 1.8 for the MT118 case. This change increased the turbulent dispersion force and improved the shape and quantitative agreement of the radial void profile at the expense of the mean value, as discussed in Section 5.4. Further optimization of the turbulent dispersion force was outside the scope of the present study. Combining the BVRD drag model with the MSZLF BIT model improved the agreement of the void fraction profile. Both Cases 25 and 26 improved the agreement of the void profile in Figure 16 but marginally worsened the prediction of velocity.

Agreement with Experimental Data
It has been demonstrated how changes in interfacial momentum and turbulence closures affect the radial profiles of void fraction and velocity. With the exception of the few instances previously discussed, the average void fraction and gas velocity do not change significantly. The CFD simulations under-predict the mean void fraction by approximately 10% -15%. The gas velocity is over-predicted in MT039, and under predicted in MT118 by about 5%, and 10%, respectively. Whether it is the radial distribution or mean value, the gas velocity is far less sensitive to a given change in closure model than the void fraction. Consequently, only the comparative agreement of void fraction is shown in Figure 17.
In order to access the physical accuracy of the present simulations, the root mean square error (RMSE) is calculated between the simulations and the appropriate experimental profile. The RMSE is normalized by the cross-sectional average of the corresponding experimental quantity. The agreement of the simulated cases is summarized in Figure 17  better predicted in MT118 than MT039. It was noted, but not plotted here, that the velocity profile is better predicted in MT039 than MT118. The agreement of the velocity profile is better than the void profile for both MT039 and MT118. The LMSB lift-wall model improves the agreement of the radial void profile for MT039 by 10% -20% from the conventional lift-wall approach, with negligible change to the mean void fraction. For MT118 cases, the quantitative agreement is similar between LMSB and conventional lift-wall models for otherwise equivalent configurations. The LMSB lift-wall model appears to be a promising advancement in closure modelling.
The results for all the cases were also compared in terms of mean void fraction against the experiments. The mean void fraction for cases MT039 and MT118 are 0.02126 and 0.189, respectively, and are shown as dashed lines in Figure  17(b). The comparison of mean void fraction and gas velocity between primary-phase turbulence and BIT models is consistent with the findings stated in Section 5.2. There was only a significant difference observed between k-ε model and SST model results when the SS BIT model was used. In that case, the k-ε model predicted a mean void fraction 3% higher in MT039 and 2.5% lower in MT118. Thus, it is concluded the source term style BIT models dominate the effect of primary phase turbulence model.
An inverse relationship between the agreement of the radial distribution and mean void fraction occurred. This was most pronounced for MT118 Cases 11,14,17, where the RK BIT model produced excellent agreement of mean void fraction but physically unreasonable void distribution and a large RSME. The same trend can be observed comparing MT039 Cases 01, 07, 08, where the wall lubrication model had a significant effect on the results. When the turbulent disper-sion force was altered between Cases 23 and 25, the RMSE improved by 10%, but the mean agreement worsened by 8%. The inverse relationship between the agreement of radial distribution and mean quantities suggests further development of the lateral redistribution force closures is required. This contradictory relationship speaks to the complex physics being modeled through averaged momentum, void fraction, and closure relations.
For the parameters considered, the turbulence closure had the most significant effect on the mean quantities. For MT039, SST + SS (Case-01) best predicts the mean void fraction, and MB is second best. For MT118, the RK and MB models best predict the mean void fraction, but only the RK model best predicts the mean gas velocity.
Although there is no case that ranks best in all comparisons, Cases 20 to 23 are considered to have performed well because of the improved prediction of near wall void fraction in the wall peak profile. With some exceptions, the MB BIT model better predicts the gas velocity than other turbulence closures tested.
For these reasons, Case 23 was selected for comparison against other published closure sets.
The best performing case from the present work was compared against the benchmark solutions for MT039 in Colombo et al. [45] presented as work from the University of Leeds (UofL) and the Helmholtz-Zentrum Dresden Rossendorf (HZDR). The UofL model and Case 23 predict the radial void profile 10% more accurately than the HZDR model. The UofL model predicts the gas velocity profile within 5%, whereas the HZDR and Case 23 are within 6.5%. The primary distinction between the UofL model, and HZDR and Case 23 is the turbulence modeling. The UofL model uses an elliptic-blending Reynolds stress transport model for primary-phase turbulence which leads to further differences in the BIT modelling. Because of this difference in turbulence modeling, the UofL model underpredicts the mean velocity and both the HZDR and the present work overpredict it. Furthermore, the mean void fraction predicted by the UofL model is 11% under the measured value, 4% better than HZDR and Case 23.
Case 23 is the only model of the three to predict a non-zero void fraction on the wall due to the LMSB lift-wall approach.

Summary and Conclusions
Vertical bubbly flow in a vertical pipe was modelled using commercial CFD code ANSYS CFX. Two experimental data sets (wall peak and core peak) from the MTLoop experiments were used for comparison. The accuracy of the two-fluid CFD model relies heavily on the accuracy of closure models. Many bubbly flow closure models new to CFX were implemented through user-defined CEL and source terms, including recently proposed closure models not yet applied to the MTLoop data set. In total, 26 closure configurations were considered, and 47 cases were simulated and compared for wall and core peak profiles. Although there is no clear best closure combination, Case