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

Geoffrey S. Gray, Scott J. Ormiston^*

Department of Mechanical Engineering, University of Manitoba, Winnipeg, Canada.

DOI: 10.4236/ojfd.2021.112007                

Abstract

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.

Keywords

Two-Phase Flow, Bubbly, Numerical Analysis, Vertical Pipe

1. 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 complex 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.

In support of improved understanding of the flow phenomena, there have been numerous detailed experimental studies of isothermal bubbly flow in vertical pipes. The studies that are most often used for validation of CFD simulations are those of Serizawa et al. [2] [3] [4], Fu [5], Liu [6] [7], Liu and Bankoff [8] [9], Hibiki and Ishii [10], Hibiki et al. [11], Hosokawa and Tomiyama [12], Shawkat et al. [13], Sun et al. [14], Lucas et al. [15], Wang et al. [16], bin Mohd Akbar et al. [17], Monrós-Andreu et al. [18], Ohnuki and Akimoto [19], Shen et al. [20], and Prasser et al. [21]. Radial profiles of volume fraction and velocity are the most commonly reported quantities. Experimental measurement of turbulence quantities is sparse; the work of Liu [7] is the most commonly cited source. More recent experimental works have focused on the measurement of bubble size, and bubble size distribution [15] [21].

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 overview of CFD modelling of bubbly flow and a summary of the many possible closure relations available for the interfacial forces. Khan et al. listed 21 possible drag force relations, 18 possible lift force relations, 9 possible wall lubrication relations, 8 possible turbulence dispersion relations, and 13 possible virtual mass relations. 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).

2. Model Description

2.1. 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

Figure 1. Computational domain: (a) side view, and (b) top view (not to scale).

by a no slip boundary for both phases. The outlet is a zero relative pressure boundary at the top of the pipe.

2.2. 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.

2.2.1. Continuity Equations

The mass conservation equations for both phases are:

$\frac{\partial \left({\alpha }_L{\rho }_L\right)}{\partial t}+\nabla \cdot \left({\alpha }_L{\rho }_L{U}_L\right)=0$ (1)

$\frac{\partial \left({\alpha }_G{\rho }_G\right)}{\partial t}+\nabla \cdot \left({\alpha }_G{\rho }_G{U}_G\right)=0$ (2)

The volume fraction constraint enforces ${\alpha }_G+{\alpha }_L=1$ in each control volume. Both phases share the same pressure field.

2.2.2. Momentum Equations

The momentum equations for both phases are:

$\frac{\partial \left({\alpha }_L{\rho }_L{U}_L\right)}{\partial t}+\nabla \cdot \left({\alpha }_L{\rho }_L{U}_L{U}_L\right)=-{\alpha }_L\nabla P+\nabla \cdot \left({\alpha }_L{T}_L\right)+{\alpha }_L{\rho }_{L}g+F_L^{\text{inter}}$ (3)

$\frac{\partial \left({\alpha }_G{\rho }_G{U}_G\right)}{\partial t}+\nabla \cdot \left({\alpha }_G{\rho }_G{U}_G{U}_G\right)=-{\alpha }_G\nabla P+\nabla \cdot \left({\alpha }_G{T}_G\right)+{\alpha }_G{\rho }_{G}g+F_G^{\text{inter}}$ (4)

The inter-phase momentum transfer consists of:

$F_G^{\text{inter}}=-F_L^{\text{inter}}=F^D+F^L+F^W+F^{TD}$ (5)

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:

$F^D=-\frac{3}{4d_b}{C}_D{\rho }_L{\alpha }_G\Vert U_R\Vert \left(U_R\right)$ (6)

where $U_R$ is the relative slip velocity, defined as:

$U_R=U_G-U_L$ (7)

The lateral lift force is expressed as [46]:

$F^L=-C_L{\rho }_L{\alpha }_G{U}_R\times \left(\nabla \times U_L\right)$ (8)

The wall lubrication force is generally expressed as [47]:

$F^W=C_{W\text{}L}{\rho }_L{\alpha }_G{\Vert U_{R\left|\right|}\Vert }^2\cdot {\hat{n}}_{\text{wall}}$ (9)

where ${\hat{n}}_{\text{wall}}$ is the wall normal unit vector, and $U_{R\left|\right|}$ is the component of relative velocity orthogonal to ${\hat{n}}_{\text{wall}}$.

The turbulent dispersion force is calculated using the Favre-averaged drag (FAD) model [48]:

$F^{TD}=-\frac{3}{4}{C}_D\frac{{\alpha }_G}{d_b}\Vert U_R\Vert \frac{{\mu }_{t\mathrm{,}L}}{{\sigma }_{TD}}\left(\frac{1}{{\alpha }_G}+\frac{1}{{\alpha }_L}\right)\nabla {\alpha }_G$ (10)

where ${\sigma }_{TD}=0.9$, except when using the Magolan and Baglietto bubble-induced turbulence model [49], as described in Section 2.6. Coefficients $C_D$, $C_L$, and $C_{W\text{}L}$ are given with the details of the interfacial force closure models in Section 2.4.

The stress tensor, $T$, includes viscous and turbulent stresses:

$T_L={\mu }_{\text{eff}\mathrm{,}L}\left(\nabla U_L+{\left(\nabla U_L\right)}^{\text{T}}\right)$ (11)

$T_G={\mu }_{\text{eff}\mathrm{,}G}\left(\nabla U_G+{\left(\nabla U_G\right)}^{\text{T}}\right)$ (12)

where the effective viscosities are ${\mu }_{\text{eff}\mathrm{,}L}={\mu }_L+{\mu }_{t\mathrm{,}L}$ and ${\mu }_{\text{eff}\mathrm{,}G}={\mu }_G+{\mu }_{t\mathrm{,}G}$.

2.3. 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):

$\begin{array}{l}\frac{\partial \left({\alpha }_L{\rho }_L{k}_L\right)}{\partial t}+\nabla \cdot \left({\alpha }_L{\rho }_L{U}_L{k}_L\right)\\ =\nabla \cdot \left({\alpha }_L\left({\mu }_L+\frac{{\mu }_{t\mathrm{,}L}}{{\sigma }_k}\right)\nabla k_L\right)+{\alpha }_L\left(P_{k\mathrm{,}L}-{\rho }_L{\epsilon }_L\right)+S_{k\mathrm{,}L}\end{array}$ (13)

$\begin{array}{l}\frac{\partial \left({\alpha }_L{\rho }_L{\epsilon }_L\right)}{\partial t}+\nabla \cdot \left({\alpha }_L{\rho }_L{U}_L{\epsilon }_L\right)\\ =\nabla \cdot \left({\alpha }_L\left({\mu }_L+\frac{{\mu }_{t\mathrm{,}L}}{{\sigma }_{\epsilon }}\right)\nabla {\epsilon }_L\right)+{\alpha }_L\frac{{\epsilon }_L}{k_L}\left(C_{\epsilon 1}{P}_{k\mathrm{,}L}-C_{\epsilon 2}{\rho }_L{\epsilon }_L\right)+S_{\epsilon \mathrm{,}L}\end{array}$ (14)

$P_{k\mathrm{,}L}=2{\mu }_{t\mathrm{,}L}{S}_L{S}_L$ (15)

where $S_L$ is the strain rate tensor of the liquid phase defined as:

$S_L=\frac{1}{2}\left(\nabla U_L+{\left(\nabla U_L\right)}^{\text{T}}\right)$ (16)

${\mu }_{t\mathrm{,}L}=C_{\mu }{C}_{\mu B}{\rho }_L\frac{k_L^2}{{\epsilon }_L}$ (17)

Source terms $S_{k\mathrm{,}L}$, $S_{\epsilon \mathrm{,}L}$, and $S_{\omega \mathrm{,}L}$ (appearing in Equation (19) later) account for BIT, and they are specified in Section 2.6. The $C_{\mu B}$ coefficient is specified later in Section 2.6. Standard model constants (as for a single phase flow) are used, as shown in Table 1.

• Shear-Stress Transport (SST):

$\begin{array}{l}\frac{\partial \left({\alpha }_L{\rho }_L{k}_L\right)}{\partial t}+\nabla \cdot \left({\alpha }_L{\rho }_L{U}_L{k}_L\right)\\ =\nabla \cdot \left({\alpha }_L\left({\mu }_L+\frac{{\mu }_{t\mathrm{,}L}}{{\sigma }_k}\right)\nabla k_L\right)+{\alpha }_L\left(P_{k\mathrm{,}L}-C_{\mu }{\rho }_L{\omega }_L{k}_L\right)+S_{k\mathrm{,}L}\end{array}$ (18)

$\begin{array}{l}\frac{\partial \left({\alpha }_L{\rho }_L{\omega }_L\right)}{\partial t}+\nabla \cdot \left({\alpha }_L{\rho }_L{U}_L{\omega }_L\right)\\ =\nabla \cdot \left({\alpha }_L\left({\mu }_L+\frac{{\mu }_{t\mathrm{,}L}}{{\sigma }_{\omega }}\right)\nabla {\omega }_L\right)+{\alpha }_L\left(C_{\omega P}\frac{{\rho }_L{P}_k}{{\mu }_t}-C_{\omega D}{\rho }_L{\omega }_L^2\right)\\ +{\alpha }_L{\xi }_{\omega }\left(1-F_1\right)+S_{\omega \mathrm{,}L}\end{array}$ (19)

where

${\xi }_{\omega }=\frac{2{\rho }_L}{{\sigma }_{\omega }{}_2}\frac{\nabla k_L\cdot \nabla {\omega }_L}{{\omega }_L}$ (20)

The blending function $F_1$ varies between 1 for k-ω at the wall, and 0 for k-ε away from the wall:

$F_1=\tanh \left({\left(\min \left[\max \left[\frac{\sqrt{k_L}}{C_{\mu }{\omega }_L{y}_{\text{wall}}}\mathrm{,}\frac{500\mu }{{\rho }_L{\omega }_L{y}_{\text{wall}}^2}\right]\mathrm{,}\frac{\left(\frac{1}{y_{\text{wall}}^2}\right)4{\rho }_L\frac{k}{{\sigma }_{\omega }{}_2}}{\max \left[{\xi }_{\omega }\mathrm{,}{\xi }_{\text{min}}\right]}\right]\right)}^4\right)$ (21)

where $\begin{array}{l}{\xi }_{\text{min}}=1\times {10}^{-10}\hfill \end{array}$ and $y_{\text{wall}}$ is the distance to the nearest wall. The turbulent viscosity is then calculated using:

${\mu }_{t\mathrm{,}L}=\frac{{\rho }_L{k}_L}{\max \left[{\omega }_L\mathrm{,}{F}_2\frac{\sqrt{2S_L{S}_L}}{C_{\omega P}}\right]}$ (22)

where

$F_2=\tanh \left({\left(\max \left[\frac{2\sqrt{k_L}}{C_{\mu }{\omega }_L{y}_{\text{wall}}}\mathrm{,}\frac{500{\mu }_L}{{\rho }_L{\omega }_L{y}_{\text{wall}}^2}\right]\right)}^2\right)$ (23)

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: $\varphi =F_1{\varphi }_1+\left(1-F_1\right){\varphi }_2$, where $\varphi$ is ${\sigma }_k$, ${\sigma }_{\omega }$, $C_{\omega P}$, or $C_{\omega D}$.

Table 1. Turbulence model constants.

Table 2. Turbulence model constants.

The turbulence in the dispersed phase is approximated from ${\mu }_{t\mathrm{,}L}$ using a zero equation model:

${\mu }_{t\mathrm{,}G}=\frac{{\rho }_G}{{\rho }_L}\frac{{\mu }_{t\mathrm{,}L}}{{\sigma }_t}$ (24)

where ${\sigma }_t$ is the ratio of liquid and gas eddy viscosities; the default value of 1 was used in this work.

2.4. 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.

2.4.1. 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 Ishii and Zuber (as implemented in CFX [50] )

$C_D=\max \left(\min \left(C_{D\mathrm{,}\text{ellipse}}\mathrm{,}{C}_{D\mathrm{,}\text{cap}}\right)\mathrm{,}{C}_{D\mathrm{,}\text{sphere}}\right)$ (25)

CFX applies the following swarm corrections to the original IZ model

$C_{D\mathrm{,}\text{sphere}}=\frac{24}{{\text{Re}}_E}\left(1+0.1{\mathrm{Re}}_E^{0.75}\right)$ (26)

$C_{D\mathrm{,}\text{cap}}={\left(1-{\alpha }_G\right)}^2{C}_{D\mathrm{,}\text{cap}}^{\infty }$ (27)

$C_{D\mathrm{,}\text{ellipse}}={\text{E}}_{\alpha }{C}_{D\mathrm{,}\text{ellipse}}^{\infty }$ (28)

${\text{Re}}_E=\frac{{\rho }_L{d}_b\Vert U_R\Vert }{{\mu }_E}$ (29)

${\text{E}}_{\alpha }={\left(\frac{1+17.67{\left(\frac{{\mu }_L}{{\mu }_E}\sqrt{{\alpha }_L}\right)}^{\frac{6}{7}}}{18.67\left(\frac{{\mu }_L}{{\mu }_E}\sqrt{{\alpha }_L}\right)}\right)}^2$ (30)

${\mu }_E={\mu }_L{\left(1-\frac{{\alpha }_G}{{\alpha }_{\text{crit}}}\right)}^{-2.5{\alpha }_{\text{crit}}{\mu }_{\ast }}$ (31)

${\mu }_{\ast }=\frac{{\mu }_G+0.4{\mu }_L}{{\mu }_G+{\mu }_L}$ (32)

where ${\alpha }_{\text{crit}}$ is the maximum packing; the default value of 1 was used.

The drag correlations for the cap and ellipse sparsely distributed regimes from IZ [52] are:

$C_{D\mathrm{,}\text{cap}}^{\infty }=\frac{8}{3}$ (33)

$C_{D\mathrm{,}\text{ellipse}}^{\infty }=\frac{2}{3}\sqrt{\text{Eo}}$ (34)

where

$\text{Eo}=\frac{g{d}_b^2\left({\rho }_L-{\rho }_G\right)}{\sigma }$ (35)

TKZS Tomiyama et al. [53]

$C_D=\max \left(\min \left(\frac{24}{{\text{Re}}_{\text{b}}}\left(1+0.15{\mathrm{Re}}_{\text{b}}^{0.687}\right)\mathrm{,}\frac{72}{{\text{Re}}_{\text{b}}}\right)\mathrm{,}\frac{8}{3}\frac{\text{Eo}}{\text{Eo}+4}\right)$ (36)

where

${\text{Re}}_{\text{b}}=\frac{{\rho }_L{d}_b\Vert U_R\Vert }{{\mu }_L}$ (37)

BVRD Buffo et al. [54]

$C_D=C_{D\mathrm{,0}}{f}_B$ (38)

$f_B=(\begin{array}{ll}{\left(1-{\alpha }_G\right)}^{C_A}\hfill & {\alpha }_G\le 0.8\hfill \\ 1\hfill & {\alpha }_G>0.8\hfill \end{array}$ (39)

$C_{D\mathrm{,0}}=\max \left(\frac{24}{{\text{Re}}_{\text{b,mod}}}\left(1+0.15{\mathrm{Re}}_{\text{b,mod}}^{0.687}\right)\mathrm{,}\frac{8}{3}\frac{\text{Eo}}{\text{Eo}+4}\right)$ (40)

${\text{Re}}_{\text{b,mod}}=\frac{{\rho }_L{d}_b\Vert U_R\Vert }{{\mu }_L+C_B{\rho }_L\frac{k_L^2}{{\epsilon }_L}}$ (41)

where $C_A=-1.3$, and $C_B=0.002$.

2.4.2. 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 Tomiyama [58] (T) is a standard lift model built into CFX that accounts for lift force reversal. The recent correlation by Ziegenhein et al. [59] (ZTL) is not currently in CFX and was implemented by a user-defined function in CFX. Although the two models have the same functional dependencies, it will be demonstrated they predict significantly different C_L under some circumstances.

T Tomiyama [58] (as implemented in CFX)

$C_L=(\begin{array}{ll}\min \left(0.288\tanh \left(0.121{\text{Re}}_{\text{b}}\right)\mathrm{,}{f}_T\right)\hfill & {\text{Eo}}_{\text{}\perp }\le 4\hfill \\ f_T\hfill & 4<{\text{Eo}}_{\text{}\perp }\le 10\hfill \\ -0.27\hfill & {\text{Eo}}_{\text{}\perp }>10\hfill \end{array}$ (42)

where:

$f_T=0.00105{\text{Eo}}_{\text{}\perp }^3-0.0159{\text{Eo}}_{\text{}\perp }^2-0.0204{\text{Eo}}_{\text{}\perp }+0.474$ (43)

and

${\text{Eo}}_{\text{}\perp }=\frac{g{d}_{\perp }^2\left({\rho }_L-{\rho }_G\right)}{\sigma }$ (44)

The lateral bubble diameter, $d_{\perp }$, 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.

$d_{\perp }=d_b\sqrt[3]{1+A{\text{Eo}}^B}$ (45)

using $A=0.163$ and $B=0.757$.

ZTL Ziegenhein et al. [59]

$C_L=(\begin{array}{ll}0.002{\text{Eo}}_{\text{}\perp }^2-0.1{\text{Eo}}_{\text{}\perp }+0.5\hfill & 1.2<{\text{Eo}}_{\text{}\perp }\le 10.5\hfill \\ -0.3295\hfill & {\text{Eo}}_{\text{}\perp }>10.5\hfill \end{array}$ (46)

${\text{Eo}}_{\text{}\perp }$ is defined by Equation (44) and $d_{\perp }$ is calculated using Equation (45) with $A=0.7$ and $B=0.7$.

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.

2.4.3. 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]

$C_{W\text{}L}=C_W\frac{d_b}{2}\left(\frac{1}{y_{\text{wall}}^2}-\frac{1}{{\left(D-y_{\text{wall}}\right)}^2}\right)$ (47)

$C_W=(\begin{array}{ll}0.47\hfill & \text{Eo}<1\hfill \\ {\text{e}}^{-0.933\text{Eo}+0.179}\hfill & 1\le \text{Eo}\le 5\hfill \\ 0.00599\text{Eo}-0.0187\hfill & \mathrm{5\; <mo>\; <\; </mo>}\text{Eo}\le 33\hfill \\ 0.179\hfill & \text{Eo}>33\hfill \end{array}$ (48)

FZKPL Frank et al. [28]

$C_{W\text{}L}=C_W\max \left(\mathrm{0,}\frac{1}{C_{wd}}\frac{\left(\mathrm{<mo\; stretchy="false">\; (\; </mo>1}-\frac{y_{\text{wall}}}{C_{wc}{d}_b}\right)}{y_{\text{wall}}{\left(\frac{y_{\text{wall}}}{C_{wc}{d}_b}\right)}^{p-1}}\right)$ (49)

where $C_W$ is calculated using Equation (48). The default model constants are $C_{wc}=10$, $C_{wd}=6.8$, and $p=1.7$.

RKL Rzehak et al. [36]

$C_{W\text{}L}=\frac{2}{d_b}0.0217\text{Eo}{\left(\frac{d_b}{y_{\text{wall}}}\right)}^2$ (50)

2.4.4. 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]

$C_L=(\begin{array}{ll}0\hfill & y_{\text{wall}}\le \frac{d_b}{2}\hfill \\ C_{L\mathrm{,0}}\left(3{\left(2\left(\frac{y_{\text{wall}}}{d_b}\right)-1\right)}^2-2{\left(2\left(\frac{y_{\text{wall}}}{d_b}\right)-1\right)}^3\right)\hfill & \frac{d_b}{2}<y_{\text{wall}}\le d_b\hfill \\ C_{L\mathrm{,0}}\hfill & y_{\text{wall}}>d_b\hfill \end{array}$ (51)

where $C_{L\mathrm{,0}}$ is a nominal lift coefficient. In the present study the lift coefficient predicted by the ZTL correlation (Equation (46)) was used.

LMSB Lubchenko et al. [62]

$F^W=(\begin{array}{ll}\frac{3}{4}{C}_D\Vert U_R\Vert \frac{{\mu }_t}{{\sigma }_{TD}}\left(1+\frac{{\alpha }_G}{1-{\alpha }_G}\right){\alpha }_G\frac{1}{y_{\text{wall}}}\frac{d_b-2y_{\text{wall}}}{d_b-y_{\text{wall}}}\cdot {\hat{n}}_{\text{wall}}\hfill & y_{\text{wall}}\le \frac{d_b}{2}\hfill \\ 0\hfill & y_{\text{wall}}>\frac{d_b}{2}\hfill \end{array}$ (52)

Equation (52) assumes the use of FAD turbulent dispersion model (Equation (10)), but other turbulent dispersion models could be used.

2.5. Summary of Interfacial Force Closure Models

Table 3 summarizes the interfacial force closure models used in the present study. The table indicates the acronym (Name) used for future reference and the functional dependency.

As mentioned in Section 2.2.2, the Favre-averaged drag (FAD) model by Burns et al. [48] was used for the turbulent dispersion force.

Table 3. Summary of interfacial force closure relations selected and their dependencies.

2.6. 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 $k_L$, and ${\epsilon }_L$ or ${\omega }_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:

$S_{k\mathrm{,}L}=K_{BI}{F}^D\cdot U_R$ (53)

$S_{\epsilon \mathrm{,}L}=C_{\epsilon B}\frac{S_{k\mathrm{,}L}}{\tau }$ (54)

$S_{\omega \mathrm{,}L}=\frac{1}{C_{\mu }{k}_L}{S}_{\epsilon \mathrm{,}L}-\frac{{\omega }_L}{k_L}{S}_{k\mathrm{,}L}$ (55)

The relation presented for $S_{\omega \mathrm{,}L}$ (Equation (55)) is simply a conversion of the $\epsilon$ source (Equation (54)) into an equivalent $\omega$ source for use in the SST model. The specifics of the BIT models are summarized in Table 4.

The Sato Sekoguchi (SS) eddy viscosity modifier approach [65] does not alter the turbulence transport equations via source terms. Instead, an extra term ( ${\mu }_t^{BI}$ ) is added to the standard definition of the eddy viscosity. The new modified eddy viscosity ( ${\tilde{\mu }}_{t\mathrm{,}L}$ ) becomes:

${\tilde{\mu }}_{t\mathrm{,}L}={\mu }_{t\mathrm{,}L}+{\mu }_t^{BI}$ (56)

Table 4. Summary of BIT model coefficients.

where ${\mu }_{t\mathrm{,}L}$ is the standard shear-induced eddy viscosity defined by (Equation (17) or Equation (22)) and ${\mu }_t^{BI}$ is calculated as:

${\mu }_t^{BI}=C_{\mu }^{BI}{\rho }_L{\alpha }_G{d}_b\Vert U_R\Vert$ (57)

using $C_{\mu }^{BI}=0.6$. The effective viscosity becomes ${\mu }_{\text{eff}\mathrm{,}L}={\mu }_L+{\tilde{\mu }}_{t\mathrm{,}L}.$

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 ${\sigma }_{TD}$ to produce more realistic results. The values of ${\sigma }_{TD}$ used in this study are given in Section 3.

2.7. 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.

3. 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

Table 5. Selected MTLoop test conditions.

flow rates at the bottom of a pipe with a diameter of 51.2 mm. A series of measurement locations along the length of the pipe were recorded. Because the modelling in this study was limited to fully developed flow, only the final measurement location at $z_m=3.03\text{m}$ was used. A wire-mesh sensor (WMS) developed by Prasser et al. [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 was 3.3 m. Constant values of density and dynamic viscosity were used for air and water: ${\rho }_G=1.185\left[\text{kg}\cdot {\text{m}}^{-3}\right]$, ${\mu }_G=1.831\times {10}^{-5}\left[\text{kg}\cdot {\text{m}}^{-1}\cdot {\text{s}}^{-1}\right]$, ${\rho }_L=997\left[\text{kg}\cdot {\text{m}}^{-3}\right]$, and ${\mu }_L=8.899\times {10}^{-4}\left[\text{kg}\cdot {\text{m}}^{-1}\cdot {\text{s}}^{-1}\right]$. The surface tension between air and water was specified as 0.072 [N∙m⁻¹]. 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 $\frac{L}{D}\approx 60$ 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 $C_{\mu B}$ coefficient to the eddy viscosity for the SST model is not as straightforward as with the k-ε model. Another constraint on the combination of interfacial forces is the novel lift-wall closures are not intermixed with the conventional lift-wall approach. For example, combining the SP lift force modification with the T wall model is not a valid configuration. 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

Table 6. Summary of parametric study cases.

^aThe turbulent dispersion model was modified.

different linearization for the user-defined approach (compared to built-in) and a resulting less robust scheme that could not be made to converge for the larger near-wall body force for MT118. Cases 25 and 26 were only deemed necessary for MT118. Because of its frequent use in the open literature, the IZ drag model is used as the baseline drag model. A value of ${\sigma }_{TD}=3.0$ was used for cases 13, 19, and 23 (for both MT039 and MT118); for Case 26 (only for MT118), it was set to 1.8, as discussed later.

4. Numerical Methods

4.1. 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 turbulence equations. The transient terms in the governing equations are retained in the solution procedure as a means of obtaining relaxation toward a steady-state solution. This pseudo-transient approach uses the time step size as a means of controlling convergence behaviour. At each time step, the linearized coupled set of mass, momentum, and volume fraction equations is solved iteratively using additive correction multigrid acceleration of an incomplete lower upper (ILU) factorization.

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⁻³ to 1 × 10⁻⁶, negligible difference was observed in the results. Radial void profiles had a 0.04% root mean square (RMS) difference between 1 × 10⁻³ and 1 × 10⁻⁶ 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⁻⁵ 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 residual 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.

4.2. 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 (Δ₁) 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.

4.3. 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

Table 7. Details of grids used for mesh sensitivity tests.

Figure 2. Cross-sectional view of working mesh.

models were implemented in CFX as user-defined momentum sources. The expression for y_wall was limited to a minimum value of 1 × 10⁻⁴ 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 $C_{\mu B}$ 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 ${\sigma }_{TD}$ cannot be adjusted directly. Instead, the turbulent dispersion coefficient was specified to obtain the desired value of ${\sigma }_{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.

5. 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 ( $z_m=3.0\text{3}\text{m}$ ). 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.

Figure 3. Comparison with Rzehak and Kriebitzsch [39] for MT039 radial profiles of: (a) void fraction, (b) axial gas velocity, and (c) liquid eddy viscosity. Closures are drag: IZ; lift: T; wall: RKL; BIT: RK; turbulence: SST; dispersion: FAD.

The lateral lift, and wall lubrication force, along with the turbulent dispersion force act normal to the streamwise direction of the flow. These lateral redistribution forces strongly influence the shape of the radial void profile and have little effect on the mean void fraction or velocity. Although the lateral lift and wall lubrication models can be varied independently, it is difficult to isolate the physical accuracy of either closure model. Because both lateral lift and wall lubrication force act normal to the wall and their effective regions overlap, the net effect of the two is what influences the prediction of the radial void profile. Consequently, any conclusions drawn with respect to the physical accuracy of a lateral lift model must consider the wall lubrication model it is paired with, and vice versa.

5.1. 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 Figure 5. The BVRD swarm modifiers are functions of local void fraction, and turbulence (proportional to the eddy viscosity) and act to increase the value of C_D. The void fraction and turbulence fields in the MT039 case are relatively small, to the extent that the BVRD drag model shows a negligible difference from the IZ or TKZS drag models in all cases considered. For the MT118 case, the BVRD swarm modifiers have a larger effect on C_D because of the increased void fraction and intensity of turbulence. Although no significant effect is observed in Figure 5 where the SS BIT model is used, the BVRD drag model does have a noticeable effect when other BIT models are used, as seen later in Section 5.2. In instances where the BVRD drag model and MSZLF BIT model are used together, the accuracy of the radial void fraction profile improved with a minor reduction in accuracy of the radial gas velocity profile. A more detailed discussion of the sensitivity of the BVRD drag model to the BIT model is provided in Section 5.2.

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

Figure 4. Effect of drag and lift models on MT039 radial profiles of: (a) void fraction and (b) axial gas velocity. Legend in (b) applies to (a). Other closures are: wall: T; BIT: SS; turbulence: SST.

Figure 5. Effect of drag and lift models on MT118 radial profiles of: (a) void fraction and (b) axial gas velocity. Legend in (b) applies to (a). Other closures are: wall: T; BIT: SS; turbulence: SST.

predict similar lift coefficients for an equivalent value of ${\text{Eo}}_{\text{}\perp }$. The lift models, however, predict different values of ${\text{Eo}}_{\text{}\perp }$ for the same d_b because different correlations for $d_{\perp }$ are used. Consequently, the ZTL model predicts a smaller lift force for the MT039 case, and a smaller critical diameter (d_b where lift force reversal occurs), as seen in Figure 6. The difference in critical diameter between lift models does not affect the simulation of the MT039 and MT118 cases, as d_b is either well below or above the critical diameter range. Hypothetically, in other flow conditions where $d_b\approx \text{5}.\text{5}\text{mm}$, the T and ZTL lift models would predict completely different radial void profiles. For the MT118 case, the ZTL model predicts a larger magnitude of C_L than the T model.

Considering the effect of different C_L on the radial void profile, the difference in C_L has a diminished effect on the radial void profile of the MT118 case (Figure 5) compared to the MT039 case (Figure 4) for two reasons. First, the difference in C_L between the T and ZTL models is approximately three times greater for MT039 than for MT118. Second, because the lift and wall forces do not oppose each other in the core peak profile like they do in the wall peak, it is not critical that the magnitude of the lift and wall lubrication forces be compatible. In the wall peak conditions of MT039, the smaller ZTL lift force is overpowered by the T wall lubrication model, resulting in a physically inconsistent void fraction profile. The results of Cases 04 to 06 do not necessarily indicate that ZTL model inaccurately predicts the lift force. This issue is discussed further in Section 5.3, when the novel lift-wall approach is used. Provided the lift force and wall lubrication force are balanced appropriately, the choice of lift model does not appreciably affect the velocity profile.

As indicated in the discussion of lift models, the wall lubrication model force does not play a significant role in the core peak profile. Because of the lift force reversal, both the lift force and the wall lubrication force are acting to skew the void fraction profile toward the tube centre. The effect of the wall lubrication force on the results for the T and FZKPL models for MT118 was negligible, so those cases (07, 09) are not plotted.

Figure 6. Lift coefficient as a function of spherical bubble diameter.

Figure 7 shows that the choice of a conventional wall lubrication model does have a noticeable effect on the radial void profile in MT039. All three wall lubrication models considered are compatible with the magnitude of the T lift force. The agreement of the void profile predicted by the wall lubrication model is dependent on the BIT model used. To illustrate this, Case 24 (from Figure 3) is also plotted in Figure 7. The centre-line void fraction of Case 24 and Case 08 differ by 10% between the RK and SS BIT models. The choice of wall lubrication model does have some impact on the velocity profile due to the change in void fraction distribution.

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 $1/y_{\text{wall}}^2$ 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.

5.2. 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 ${\mu }_t$. Considering Case 01 for MT039, the SS eddy viscosity modifier produces a large value of ${\mu }_t$ relative to other BIT models, leading to the much flatter void profile seen in Figure 8. In the case of MT118, the ${\mu }_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

Figure 7. Effect of drag and wall models on MT039 radial profiles of: (a) void fraction and (b) axial gas velocity. Legend in (b) applies to (a). Other closures are: lift: T; BIT: SS; turbulence: SST.

Figure 8. Effect of drag and BIT models on MT039 radial profiles of: (a) void fraction and (b) axial gas velocity. Legend in (b) applies to (a). Other closures are: lift: T; wall: T; turbulence: SST (except keps for Cases 13 and 19).

Figure 9. Effect of drag and BIT models on MT118 radial profiles of: (a) void fraction and (b) axial gas velocity. Legend in (b) applies to (a). Other closures are: lift: T; wall: T; turbulence: SST.

turbulent dispersion force was re-calibrated for the MB BIT model by adjusting ${\sigma }_{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 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 ${\mu }_{t\mathrm{,}L}$ 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 opposite trend for Cases 01 and 16 between MT039 (Figure 10) and MT118 (Figure 13) is not contradictory; rather, it is a result of the lift force reversal which alters the behaviour of turbulent dispersion force.

The MB BIT model is the only BIT model to modify the eddy viscosity. The use of $C_{\mu B}$ increases the eddy viscosity with increase in local volume fraction. In the regions of zero void fraction, $C_{\mu B}=1$, so the standard single phase form

Figure 10. Effect of BIT and primary phase turbulence models on MT039 radial profiles of: (a) void fraction and (b) axial gas velocity. Legend in (b) applies to (a). Other closures are: drag: IZ; lift: T; wall: T.

Figure 11. Effect of BIT and primary phase turbulence models on MT039 radial profiles of: (a) eddy viscosity, (b) turbulent kinetic energy, and (c) dissipation rate. Legend in (c) applies to (a) and (b). Other closures are: drag: IZ; lift: T; wall: T.

Figure 12. Effect of BIT and primary phase turbulence models on MT118 radial profiles of: (a) eddy viscosity, (b) turbulent kinetic energy, and (c) dissipation rate. Legend in (c) applies to (a) and (b). Other closures are: drag: IZ; lift: T; wall: T.

Figure 13. Effect of BIT and primary phase turbulence models on MT118 radial profiles of: (a) void fraction and (b) axial gas velocity. Legend in (b) applies to (a). Other closures are: drag: IZ; lift: T; wall: T.

is regained. The effects of $C_{\mu B}$ are clear in Figure 12 where the MB and MSZLF model results have similar curves for k and ε, yet substantially different profiles of ${\mu }_t$.

As previously stated, the BVRD drag modifier is a function of local turbulence and void fraction, which causes a complex interaction between momentum and turbulence equations. Turbulence and void fraction is much lower in the MT039 case than in the MT118 case. Consequently, the BVRD drag modifiers have negligible effect for the MT039 case. The SS eddy viscosity modifier does not alter the k and ε or ω equations, so k²/ε is relatively small despite the value of ${\mu }_{t\mathrm{,}L}$ seen in Figure 12. Thus, the BVRD drag modifier has negligible effect in Case 02, as seen in Figure 9. When the BVRD drag model is used in conjunction with the MSZLF BIT model, however, the agreement of the radial profile and mean void fraction is improved, with negligible effect on velocity profile. Feedback between the drag force and $S_k$ occurs when any of the source term style BIT models is used. The MSZLF closure is the only BIT model to have an additional feedback between the drag model and $S_{\epsilon }$. The couplings between drag, turbulence, and void fraction result in an increased void fraction and decreased eddy viscosity at the centre line, seen in Figure 12. Because this effect appears to be unique to the MSZLF dissipation source term and no measurements of turbulence quantities are available for comparison, it remains uncertain if the improvements made with combination of BVRD and MSZLF have physical grounds.

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 $k_L$ and ${\epsilon }_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.

5.3. Effect of Novel Lift-Wall Approach

The difference between the conventional and LMSB lift-wall approach is highlighted by plotting the radial distribution of lift and wall forces, as seen in Figure 14. Theoretically, the conventional lift and wall forces are at a maximum at the wall. In Cases 01 and 04, however, the forces go to zero at the wall because the solution produced ${\alpha }_G=0$ at the wall. The effect of the SP model near-wall modification of the lift coefficient is demonstrated by the lift force curves of Case 04 (standard) and Case 24 (modified). It is noted that the lateral lift force in Case 24 goes to zero due to the SP model modification in the region where the conventional lift force is at a maximum. Because the lift force near the wall is reduced 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.

Figure 14. Radial profiles of lift and wall force for select cases.

Despite the difference in magnitude of the individual forces between the two lift-wall approaches, the corresponding radial void profiles are not that dissimilar. This similarity occurs because it is the net difference, not the magnitude of lift or wall force alone, that influences the radial void profile.

Figure 15 shows that the LMSB wall treatment produces a non-zero void fraction in the region immediately adjacent to the wall in MT039, which is physically more realistic. The reduction of the lift force near the wall, however, does produce some undesirable inflections in the void profile for MT118, as seen in Figure 16. The LMSB lift-wall approach has a more noticeable impact than the conventional wall models on MT118 because of the near wall lift force suppression. The lift force of Case 20 is zero for $r/R\ge 0.85$, where it is at a maximum in Case 01. Consequently for core-peak flows, the wall force becomes the dominant mechanism controlling the void fraction near the wall when the LMSB wall model is used whereas the lift force is the dominant mechanism when the conventional approach is used. Furthermore, the effective distance away from the wall increases with increased bubble diameter. The prediction of the velocity profile is marginally improved with the LMSB approach for both MT039 and MT118. This is a promising result because the agreement of the void and velocity profile is often inversely related.

Figure 15. Effect of BIT models with novel lift-wall approach on MT039 radial profiles of: (a) void fraction and (b) axial gas velocity. Legend in (b) applies to (a). Other closures are: drag: IZ; lift: SP; wall: LMSB; turbulence: SST.

Figure 16. Effect of BIT models with novel lift-wall approach on MT118 radial profiles of: (a) void fraction and (b) axial gas velocity. Legend in (b) applies to (a). Other closures are: drag: IZ; lift: SP; wall: LMSB; turbulence: SST.

The LMSB wall force is proportional to ${\mu }_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 ${\sigma }_{TD}$ for Case 26. Magolan [71] suggested a variable ${\sigma }_{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, ${\sigma }_{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.

5.4. 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(a). In general, the void fraction profile is

Figure 17. Comparison of cases with experimental data: (a) RMSE for radial void fraction, and (b) mean void fraction values compared to experimental values.

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 dispersion 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.

6. 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 23 was selected for comparison against the benchmark solutions published by Colombo et al. [45]. The prediction of the present CFD model for MT039 was on par with the published results.

The key findings are as follows:

• The novel LMSB lift-wall closure approach significantly improved the prediction of radial void fraction profile for wall peak conditions. The LMSB wall force is the only model tested that predicts non-zero void fraction at the wall.

• The MB BIT model gave better predictions when used with the LMSB lift-wall approach.

• The ZTL lift model is overpowered by the T wall lubrication model in wall peak conditions.

• The BVRD drag modifier combined with the MSZLF BIT model improved the agreement of the void profile for core peak conditions. Otherwise there was no significant difference observed between drag models.

• Significant difference between k-ε and SST turbulence models was only observed when the SS eddy viscosity modifier was used.

• Of all the closure parameters, the BIT model has the greatest effect on the mean void fraction and gas velocity.

• Agreement of radial distribution and mean void fraction were inversely related with change in turbulence modelling.

The behaviour of the LMSB lift-wall approach should be examined for a polydispersed simulation of MT118. Because the mean diameter would be smaller at the wall and increase towards the centre of the pipe, the active distance of the lift force reduction and LMSB wall force would vary. This is likely to reduce the undesirable inflections observed in the radial void profile of MT118 when the novel lift-wall approach was used.

Acknowledgements

The financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC) is gratefully acknowledged. The authors also recognize the support of the University of Manitoba for the use of high performance computers. Special thanks to Dr. Hassan M. Soliman for his comments on the manuscript.

Nomenclature

Symbols

A, B coefficients in lateral bubble diameter equation

C coefficient in various equations

D pipe diameter, [m]

d_b spherical bubble diameter, [m]

$d_{\perp }$ bubble lateral diameter, [m]

Eo Eötvös number

${\text{E}}_{\alpha }$ coefficient in IZ dragclosure, [kg∙m⁻³∙s⁻²]

$F$ Force per unit volume, [N∙m⁻³]

$f_B,f_T$ functions in drag and lift closure relations

F₁ blending function in SST turbulence model

g gravitational acceleration [m∙s⁻²]

J superficial velocity, [m∙s⁻¹]

k turbulence kinetic energy, [m²∙s⁻²]

K_BI coefficient in BIT source term for k

L tube length, [m]

${\hat{n}}_{\text{wall}}$ wall normal unit vector

P pressure, [N∙m⁻²]

P_k turbulence production, [N∙m⁻²∙s⁻¹]

R tube radius, [m]

Re Reynolds number

$S$ strain rate, [s⁻¹]

t time, [s]

$T$ stress tensor, [N∙m⁻²]

U velocity, [m∙s⁻¹]

N number of nodes

w axial velocity, [m∙s⁻¹]

x, y, z Cartesian coordinates

$\langle \alpha \rangle$ area average of $\alpha$

Greek symbols

$\alpha$ volume fraction

$\epsilon$ turbulence dissipation rate, [m²∙s⁻³]

$\mu$ dynamic viscosity, [kg∙m⁻¹∙s⁻¹]

${\mu }_t$ eddy viscosity, [kg∙m⁻¹∙s⁻¹]

$\rho$ density, [kg∙m⁻³]

$\sigma$ surface tension, [N∙m⁻¹]

${\sigma }_{TD}$ coefficient in turbulence dispersion force closure

$\omega$ turbulence specific dissipation rate, [s⁻¹]

Subscripts

b bubble

eff effective

G gas phase

k k equation

L liquid phase

r radial direction

R relative

t turbulent

TD Turbulence Dispersion

$\perp$ lateral

$\epsilon$ $\epsilon$ equation

$\omega$ $\omega$ equation

Superscripts

inter interfacial

D drag

L lift

W wall lubrication

TD turbulence dispersion

BI bubble induced turbulence

$\infty$ sparse bubble limit

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ishii, M. and Hibiki, T. (2011) Thermo-Fluid Dynamics of Two-Phase Flow. 2nd Edition, Springer Verlag, New York. https://doi.org/10.1007/978-1-4419-7985-8
[2]	Serizawa, A., Kataoka, I. and Michiyoshi, I. (1975) Turbulence Structure of Air-Water Bubbly Flow—I. Measuring Techniques. International Journal of Multiphase Flow, 2, 221-233. https://doi.org/10.1016/0301-9322(75)90011-7
[3]	Serizawa, A., Kataoka, I. and Michiyoshi, I. (1975) Turbulence Structure of Air-Water Bubbly Flow—II. Local Properties. International Journal of Multiphase Flow, 2, 235-246. https://doi.org/10.1016/0301-9322(75)90012-9
[4]	Serizawa, A., Kataoka, I. and Michiyoshi, I. (1975) Turbulence Structure of Air-Water Bubbly Flow—III. Transport Properties. International Journal of Multiphase Flow, 2, 247-259. https://doi.org/10.1016/0301-9322(75)90013-0
[5]	Fu, X. (2001) Interfacial Area Measurement and Transport Modeling in Air-Water Two-Phase Flow. Ph.D. Thesis, Purdue University, West Lafayette, IN, USA.
[6]	Liu, T. (1993) Bubble Size and Entrance Length Effects on Void Development in a Vertical Channel. International Journal of Multiphase Flow, 19, 99-113. https://doi.org/10.1016/0301-9322(93)90026-Q
[7]	Liu, T. (1998) The Role of Bubble Size on Liquid Phase Turbulent Structure in Two-Phase Bubbly Flow. Proceedings of the 3rd International Conference on Multiphase Flow, 8-12 June 1998, Lyon, 8-12.
[8]	Liu, T. and Bankoff, S. (1993) Structure of Air-Water Bubbly Flow in a Vertical Pipe—I. Liquid Mean Velocity and Turbulence Measurements. International Journal of Heat and Mass Transfer, 36, 1049-1060. https://doi.org/10.1016/S0017-9310(05)80289-3
[9]	Liu, T. and Bankoff, S. (1993) Structure of Air-Water Bubbly Flow in a Vertical Pipe—II. Void Fraction, Bubble Velocity and Bubble Size Distribution. International Journal of Heat and Mass Transfer, 36, 1061-1072. https://doi.org/10.1016/S0017-9310(05)80290-X
[10]	Hibiki, T. and Ishii, M. (1999) Experimental Study on Interfacial Area Transport in Bubbly Two-Phase Flows. International Journal of Heat and Mass Transfer, 42, 3019-3035. https://doi.org/10.1016/S0017-9310(99)00014-9
[11]	Hibiki, T., Ishii, M. and Xiao, Z. (2001) Axial Interfacial Area Transport of Vertical Bubbly Flows. International Journal of Heat and Mass Transfer, 44, 1869-1888. https://doi.org/10.1016/S0017-9310(00)00232-5
[12]	Hosokawa, S. and Tomiyama, A. (2009) Multi-Fluid Simulation of Turbulent Bubbly Pipe Flows. Chemical Engineering Science, 64, 5308-5318. https://doi.org/10.1016/j.ces.2009.09.017
[13]	Shawkat, M., Ching, C. and Shoukri, M. (2008) Bubble and Liquid Turbulence Characteristics of Bubbly Flow in a Large Diameter Vertical Pipe. International Journal of Multiphase Flow, 34, 767-785. https://doi.org/10.1016/j.ijmultiphaseflow.2008.01.007
[14]	Sun, X., Paranjape, S., Kim, S., Ozar, B. and Ishii, M. (2004) Liquid Velocity in Upward and Downward Air-Water Flows. Annals of Nuclear Energy, 31, 357-373. https://doi.org/10.1016/j.anucene.2003.08.002
[15]	Lucas, D., Krepper, E. and Prasser, H.M. (2005) Development of Co-Current Air-Water Flow in a Vertical Pipe. International Journal of Multiphase Flow, 31, 1304-1328. https://doi.org/10.1016/j.ijmultiphaseflow.2005.07.004
[16]	Wang, S., Lee, S., Jones, O. and Lahey, R. (1987) 3-D Turbulence Structure and Phase Distribution Measurements in Bubbly Two-Phase Flows. International Journal of Multiphase Flow, 13, 327-343. https://doi.org/10.1016/0301-9322(87)90052-8
[17]	bin Mohd Akbar, M.H., Hayashi, K., Hosokawa, S. and Tomiyama, A. (2012) Bubble Tracking Simulation of Bubble-Induced Pseudo Turbulence. Multiphase Science and Technology, 24, 197-222. https://doi.org/10.1615/MultScienTechn.v24.i3.20
[18]	Monrós-Andreu, G., Chiva, S., Martínez-Cuenca, R., Torró, S., Juliá, J.E., Hernández, L. and Mondragón, R. (2013) Water Temperature Effect on Upward Air-Water Flow in a Vertical Pipe: Local Measurements Database Using Four-Sensor Conductivity Probes and LDA. EPJ Web of Conferences, 45, Article No. 01105. https://doi.org/10.1051/epjconf/20134501105
[19]	Ohnuki, A. and Akimoto, H. (2000) Experimental Study on Transition of Flow Pattern and Phase Distribution in Upward Air-Water Two-Phase Flow along a Large Vertical Pipe. International Journal of Multiphase Flow, 26, 367-386. https://doi.org/10.1016/S0301-9322(99)00024-5
[20]	Shen, X., Mishima, K. and Nakamura, H. (2005) Two-Phase Phase Distribution in a Vertical Large Diameter Pipe. International Journal of Heat and Mass Transfer, 48, 211-225. https://doi.org/10.1016/j.ijheatmasstransfer.2004.06.034
[21]	Prasser, H.M., Beyer, M., Carl, H., Gregor, S., Lucas, D., Pietruske, H., Schutz, P. and Weiss, F.P. (2007) Evolution of the Structure of a Gas-Liquid Two-Phase Flow in a Large Vertical Pipe. Nuclear Engineering and Design, 237, 1848-1861. https://doi.org/10.1016/j.nucengdes.2007.02.018
[22]	Khan, I., Wang, M., Zhang, Y., Tian, W., Su, G. and Qiu, S. (2020) Two-Phase Bubbly Flow Simulation Using CFD Method: A Review of Models for Interfacial Forces. Progress in Nuclear Energy, 125, Article ID: 103360. https://doi.org/10.1016/j.pnucene.2020.103360
[23]	Chuang, T.J. and Hibiki, T. (2017) Interfacial Forces Used in Two-Phase Flow Numerical Simulation. International Journal of Heat and Mass Transfer, 113, 741-754. https://doi.org/10.1016/j.ijheatmasstransfer.2017.05.062
[24]	Lucas, D., Rzehak, R., Krepper, E., Ziegenhein, T., Liao, Y., Kriebitzsch, S. and Apanasevich, P. (2016) A Strategy for the Qualification of Multi-Fluid Approaches for Nuclear Reactor Safety. Nuclear Engineering and Design, 299, 2-11. https://doi.org/10.1016/j.nucengdes.2015.07.007
[25]	Lucas, D., Krepper, E., Liao, Y., Rzehak, R. and Ziegenhein, T. (2020) General Guideline for Closure Model Development for Gas-Liquid Flows in the Multi-Fluid Framework. Nuclear Engineering and Design, 357, Article ID: 110396. https://doi.org/10.1016/j.nucengdes.2019.110396
[26]	Liao, Y., Ma, T., Liu, L., Ziegenhein, T., Krepper, E. and Lucas, D. (2018) Eulerian Modelling of Turbulent Bubbly Flow Based on a Baseline Closure Concept. Nuclear Engineering and Design, 337, 450-459. https://doi.org/10.1016/j.nucengdes.2018.07.021
[27]	Lucas, D., Krepper, E. and Prasser, H.M. (2007) Use of Models for Lift, Wall and Turbulent Dispersion Forces Acting on Bubbles for Poly-Disperse Flows. Chemical Engineering Science, 62, 4146-4157. https://doi.org/10.1016/j.ces.2007.04.035
[28]	Frank, T., Zwart, P., Krepper, E., Prasser, H.M. and Lucas, D. (2008) Validation of CFD Models for Mono- and Polydisperse Air-Water Two-Phase Flows in Pipes. Nuclear Engineering and Design, 238, 647-659. https://doi.org/10.1016/j.nucengdes.2007.02.056
[29]	Duan, X., Cheung, S., Yeoh, G., Tu, J., Krepper, E. and Lucas, D. (2011) Gas-Liquid Flows in Medium and Large Vertical Pipes. Chemical Engineering Science, 66, 872-883. https://doi.org/10.1016/j.ces.2010.11.031
[30]	Colombo, M. and Fairweather, M. (2016) RANS Simulation of Bubble Coalescence and Break-up in Bubbly Two-Phase Flows. Chemical Engineering Science, 146, 207-225. https://doi.org/10.1016/j.ces.2016.02.034
[31]	Pena-Monferrer, C., Passalacqua, A., Chiva, S. and Munoz-Cobo, J. (2016) CFD Modelling and Validation of Upward Bubbly Flow in an Adiabatic Vertical Pipe Using the Quadrature Method of Moments. Nuclear Engineering and Design, 301, 320-332. https://doi.org/10.1016/j.nucengdes.2016.03.006
[32]	Liao, Y., Lucas, D., Krepper, E. and Schmidtke, M. (2011) Development of a Generalized Coalescence and Breakup Closure for the Inhomogeneous MUSIG Model. Nuclear Engineering and Design, 241, 1024-1033. https://doi.org/10.1016/j.nucengdes.2010.04.025
[33]	Liao, Y., Rzehak, R., Lucas, D. and Krepper, E. (2015) Baseline Closure Model for Dispersed Bubbly Flow: Bubble Coalescence and Breakup. Chemical Engineering Science, 122, 336-349. https://doi.org/10.1016/j.ces.2014.09.042
[34]	Rzehak, R. and Krepper, E. (2015) Bubbly Flows with Fixed Polydispersity: Validation of a Baseline Closure Model. Nuclear Engineering and Design, 287, 108-118. https://doi.org/10.1016/j.nucengdes.2015.03.005
[35]	Krepper, E., Lucas, D., Frank, T., Prasser, H.M. and Zwart, P.J. (2008) The Inhomogeneous MUSIG Model for the Simulation of Polydispersed Flows. Nuclear Engineering and Design, 238, 1690-1702. https://doi.org/10.1016/j.nucengdes.2008.01.004
[36]	Rzehak, R., Krepper, E. and Lifante, C. (2012) Comparative Study of Wallforce Models for the Simulation of Bubbly Flows. Nuclear Engineering and Design, 253, 41-49. https://doi.org/10.1016/j.nucengdes.2012.07.009
[37]	Rzehak, R. and Krepper, E. (2013) Closure Models for Turbulent Bubbly Flows: A CFD Study. Nuclear Engineering and Design, 265, 701-711. https://doi.org/10.1016/j.nucengdes.2013.09.003
[38]	Rzehak, R. and Krepper, E. (2013) CFD Modeling of Bubble-Induced Turbulence. International Journal of Multiphase Flow, 55, 138-155. https://doi.org/10.1016/j.ijmultiphaseflow.2013.04.007
[39]	Rzehak, R. and Kriebitzsch, S. (2015) Multiphase CFD-Simulation of Bubbly Pipe Flow: A Code Comparison. International Journal of Multiphase Flow, 68, 135-152. https://doi.org/10.1016/j.ijmultiphaseflow.2014.09.005
[40]	Rzehak, R., Ziegenhein, T., Kriebitzsch, S., Krepper, E. and Lucas, D. (2017) Unified Modeling of Bubbly Flows in Pipes, Bubble Columns, and Airlift Columns. Chemical Engineering Science, 157, 147-158. https://doi.org/10.1016/j.ces.2016.04.056
[41]	Krepper, E., Rzehak, R. and Lucas, D. (2018) Validation of a Closure Model Framework for Turbulent Bubbly Two-Phase Flow in Different Flow Situations. Nuclear Engineering and Design, 340, 388-404. https://doi.org/10.1016/j.nucengdes.2018.09.042
[42]	Yamoah, S., Martínez-Cuenca, R., Monrós, G., Chiva, S. and Macián-Juan, R. (2015) Numerical Investigation of Models for Drag, Lift, Wall Lubrication and Turbulent Dispersion Forces for the Simulation of Gas-Liquid Two-Phase Flow. Chemical Engineering Research and Design, 98, 17-35. https://doi.org/10.1016/j.cherd.2015.04.007
[43]	Wang, Q. and Yao, W. (2016) Computation and Validation of the Interphase Force Models for Bubbly Flow. International Journal of Heat and Mass Transfer, 98, 799-813. https://doi.org/10.1016/j.ijheatmasstransfer.2016.03.064
[44]	Lote, D.A., Vinod, V. and Patwardhan, A.W. (2018) Computational Fluid Dynamics Simulations of the Air-Water Two-Phase Vertically Upward Bubbly Flow in Pipes. Industrial & Engineering Chemistry Research, 57, 10609-10627. https://doi.org/10.1021/acs.iecr.8b01579
[45]	Colombo, M., Rzehak, R., Fairweather, M., Liao, Y. and Lucas, D. (2021) Benchmarking of Computational Fluid Dynamic Models for Bubbly Flows. Nuclear Engineering and Design, 375, Article ID: 111075. https://doi.org/10.1016/j.nucengdes.2021.111075
[46]	Zun, I. (1980) The Transverse Migration of Bubbles Influenced by Walls in Vertical Bubbly Flow. International Journal of Multiphase Flow, 6, 583-588. https://doi.org/10.1016/0301-9322(80)90053-1
[47]	Antal, S., Lahey, R. and Flaherty, J. (1991) Analysis of Phase Distribution in Fully Developed Laminar Bubbly Two-Phase Flow. International Journal of Multiphase Flow, 17, 635-652. https://doi.org/10.1016/0301-9322(91)90029-3
[48]	Burns, A., Frank, T., Hamill, I. and Shi, J.M. (2004) The Favre Averaged Drag Model for Turbulent Dispersion in Eulerian Multi-Phase Flows. Proceedings of the 5th International Conference on Multiphase Flow, Yokohama, 30 May-4 June 2004, Paper No. 392.
[49]	Magolan, B. and Baglietto, E. (2019) Assembling a Bubble-Induced Turbulence Model Incorporating Physical Understanding from DNS. International Journal of Multiphase Flow, 116, 185-202. https://doi.org/10.1016/j.ijmultiphaseflow.2019.04.009
[50]	Ansys, Inc. (2020) ANSYS CFX-Solver Theory Guide Release 2020-R1. ANSYS, Inc., Canonsburg, PA, USA.
[51]	Schiller, L. and Naumann, A. (1935) A Drag Coefficient Correlation. Zeitschrift des Vereins Deutscher Ingenieure, 77, 318-320.
[52]	Ishii, M. and Zuber, N. (1979) Drag Coefficient and Relative Velocity in Bubbly, Droplet or Particulate Flows. AIChE Journal, 25, 843-855. https://doi.org/10.1002/aic.690250513
[53]	Tomiyama, A., Kataoka, I., Zun, I. and Sakaguchi, T. (1998) Drag Coefficients of Single Bubbles under Normal and Micro Gravity Conditions. JSME International Journal Series B, 41, 472-479. https://doi.org/10.1299/jsmeb.41.472
[54]	Buffo, A., Vanni, M., Renze, P. and Marchisio, D. (2016) Empirical Drag Closure for Polydisperse Gas-Liquid Systems in Bubbly Flow Regime: Bubble Swarm and Micro-Scale Turbulence. Chemical Engineering Research and Design, 113, 284-303. https://doi.org/10.1016/j.cherd.2016.08.004
[55]	Lucas, D. and Tomiyama, A. (2011) On the Role of the Lateral Lift Force in Poly-Dispersed Bubbly Flows. International Journal of Multiphase Flow, 37, 1178-1190. https://doi.org/10.1016/j.ijmultiphaseflow.2011.05.009
[56]	Legendre, D. and Magnaudet, J. (1998) The Lift Force on a Spherical Bubble in a Viscous Linear Shear Flow. Journal of Fluid Mechanics, 368, 81-126. https://doi.org/10.1017/S0022112098001621
[57]	Mei, R. and Klausner, J. (1994) Shear Lift Force on Spherical Bubbles. International Journal of Heat and Fluid Flow, 15, 62-65. https://doi.org/10.1016/0142-727X(94)90031-0
[58]	Tomiyama, A. (1998) Struggle with Computational Bubble Dynamics. Multiphase Science and Technology, 10, 369-405. https://doi.org/10.1615/MultScienTechn.v10.i4.40
[59]	Ziegenhein, T., Tomiyama, A. and Lucas, D. (2018) A New Measuring Concept to Determine the Lift Force for Distorted Bubbles in Low Morton Number System: Results for Air/Water. International Journal of Multiphase Flow, 108, 11-24. https://doi.org/10.1016/j.ijmultiphaseflow.2018.06.012
[60]	Wellek, R., Agrawal, A. and Skelland, A. (1966) Shape of Liquid Drops Moving in Liquid Media. AIChE Journal, 12, 854-862. https://doi.org/10.1002/aic.690120506
[61]	Tomiyama, A., Sou, A., Zun, I., Kanami, N. and Sakaguchi, T. (1995) Effects of Eotvos Number and Dimensionless Liquid Volumetric Flux on Lateral Motion of a Bubble in a Laminar Duct Flow. In: Serizawa, A., Fukano, T. and Bataille, J., Eds., Multiphase Flow 1995, Elsevier Science, Amsterdam, 3-15. https://doi.org/10.1016/B978-0-444-81811-9.50005-6
[62]	Lubchenko, N., Magolan, B., Sugrue, R. and Baglietto, E. (2018) A More Fundamental Wall Lubrication Force from Turbulent Dispersion Regularization for Multiphase CFD Applications. International Journal of Multiphase Flow, 98, 36-44. https://doi.org/10.1016/j.ijmultiphaseflow.2017.09.003
[63]	Shaver, D.R. and Podowski, M.Z. (2015) Modeling of Interfacial Forces for Bubbly Flows in Subcooled Boiling Conditions. Transactions of the American Nuclear Society, 113, 1368-1371.
[64]	Rzehak, R. and Krepper, E. (2013) Bubble-Induced Turbulence: Comparison of CFD Models. Nuclear Engineering and Design, 258, 57-65. https://doi.org/10.1016/j.nucengdes.2013.02.008
[65]	Sato, Y. and Sekoguchi, K. (1975) Liquid Velocity Distribution in Two-Phase Bubble Flow. International Journal of Multiphase Flow, 2, 79-95. https://doi.org/10.1016/0301-9322(75)90030-0
[66]	Ma, T., Santarelli, C., Ziegenhein, T., Lucas, D. and Frohlich, J. (2017) Direct Numerical Simulation-Based Reynolds-Averaged Closure for Bubble-Induced Turbulence. Physical Review Fluids, 2, Article ID: 034301. https://doi.org/10.1103/PhysRevFluids.2.034301
[67]	Prasser, H.M., Scholz, D. and Zippe, C. (2001) Bubble Size Measurement Using Wire-Mesh Sensors. Flow Measurement and Instrumentation, 12, 299-312. https://doi.org/10.1016/S0955-5986(00)00046-7
[68]	Prasser, H.M., Misawa, M. and Tiseanu, I. (2005) Comparison between Wire-Mesh Sensor and Ultra-Fast X-Ray Tomography for an Air-Water Flow in a Vertical Pipe. Flow Measurement and Instrumentation, 16, 73-83. https://doi.org/10.1016/j.flowmeasinst.2005.02.003
[69]	Magolan, B., Lubchenko, N. and Baglietto, E. (2019) A Quantitative and Generalized Assessment of Bubble-Induced Turbulence Models for Gas-Liquid Systems. Chemical Engineering Science: X, 2, Article ID: 100009. https://doi.org/10.1016/j.cesx.2019.100009
[70]	Hosokawa, S., Tomiyama, A., Misaki, S. and Hamada, T. (2002) Lateral Migration of Single Bubbles Due to the Presence of Wall. Proceedings of the ASME Fluids Engineering Division Conference, Vol 1: Fora, Parts A and B, 14-18 July 2002. Montreal, 855-860. https://doi.org/10.1115/FEDSM2002-31148
[71]	Magolan, B. (2018) Extending Bubble-Induced Turbulence Modeling Applicability in CFD through Incorporation of DNS Understanding. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge.