Development of Bubble Size Correlation for Adiabatic Forced Convective Bubbly Flow in Low Pressure Condition Using CFD Code

: In a multidimensional two-phase ﬂow analysis, bubble size signiﬁcantly a ﬀ ects interfacial transfer terms such as mass, momentum, and energy. With regard to bubbly ﬂow, the application of a simple correlation-type bubble size model presents certain advantages, including short calculation times and ease of usage. In this study, we propose a semi-theoretical correlation developed from a steady state bubble number density transport equation for predicting the distribution of local bubble size using a computational ﬂuid dynamics (CFD) code. The coe ﬃ cients of the new correlation were determined using the local bubble parameters obtained on the basis of three existing vertical air-water experiments. Finally, these were implemented in commercial CFD code and evaluated against experimental data, which showed that the proposed correlation exhibits good prediction capability for forced convective air-water bubbly ﬂows under low pressure conditions.


Introduction
Two-phase flows are commonly employed in multiphase contactors and reactors in the chemical, nuclear, biochemical, and petrochemical sectors, among several others.Gas-liquid two-phase flows in vertical pipes can have different flow regimes, such as bubbly, slug, churn, and annular flows, depending on the flow conditions and flow channel geometry.Bubbly flow is the most common flow regime in industrial applications and has a significant advantage in terms of interfacial mass transfer between two phases because of its large interfacial area.To evaluate or increase the efficiency of systems such as steam generators and heat exchangers with such a regime, it is necessary to accurately predict bubble behavior.
For 20 years, the application of computational fluid dynamics (CFD) codes has been extended to the analysis of multidimensional two-phase flows.Many researchers have analyzed bubbly flow with a strong focus on the distribution of void fraction, interfacial forces, and turbulent model.Wang and Sun [1] applied an interfacial area transport equation (IATE) proposed by Ishii et al. [2] to Fluent to investigate their effects on two-phase flow calculations under two types of boundary conditions at the inlet: uniform and non-uniform void fraction profiles, measured via experiments.They reported similar results at the developed region under adiabatic bubbly flow conditions, regardless of the type of boundary condition.Furthermore, the lift force is essential to obtain an accurate radial distribution.In this regard, Magnaudet and Legendre [3] and Chahed and Masbernat [4] performed studies on lift force for a wide range of Reynolds numbers and interfacial forces near the wall, respectively.
In addition, studies and evaluations of forces acting on the interface between the two phases have been actively conducted [5][6][7][8][9][10].Gemello et al. [11], through CFD analysis, proposed an appropriate set of breakup and coalescence kernels for a population balance model for bubble columns with a heterogeneous flow regime.Marfaing et al. [12] studied the prediction of an adiabatic bubbly flow in a vertical pipe using the NEPTUNE code, where constant coefficients are applied to the interfacial transfer terms, such as drag, lift, and wall lubrication forces.Their results showed that the distribution of the void fraction is mainly affected by the balance between the interfacial forces.Particularly, in the case of the bubble lift force, it depends on the direction of the liquid velocity.Further, it was noted that the influences of dispersion and wall lubrication forces were dominant in the vicinity of the wall.Yamoah et al. [8] evaluated the models for drag, lift, wall lubrication, and turbulent dispersion forces using Ansys-CFX for a vertical gas-liquid flow in a pipe.They suggested baseline models for the interfacial forces based on the evaluation results of radial distributions of void fraction, interfacial area concentration (IAC), and phase velocity against the experimental data.
As revealed in the aforementioned applications, bubble size is a common parameter in the model of interfacial forces, and it significantly affects the prediction accuracy of multidimensional bubble behaviors in two-phase flows because it determines the mass, momentum, and energy transfer between phases.In particular, the interfacial transfer terms are proportional to the driving forces multiplied by the IAC, which has the following relation for the spherical bubble expected in a bubby flow.
where a i is the IAC, D b is the bubble size, and α is the void fraction.Many mechanistic bubble size models have been proposed to predict bubble size (or IAC) in a flow channel accurately.Ishii et al. [2] developed a one-dimensional interfacial area transport equation including sink and source terms which dynamically predicts the axial variations of interfacial area concentration.Wu et al. [13], Hibiki and Ishii [14], as well as Yao and Morel [15] put forth constitutive models for the source and sink terms that compose the IATE.Subsequently, Yeoh and Tu [16] evaluated the performance of the bubble number density transport equation with the constitutive models in the simulation of adiabatic air-water flows.Lo and Zhang [17] and Yun et al. [18] attempted to apply the S γ model, which is defined as a generalized moment parameter applicable to air-water and subcooled boiling flows.Additionally, Duan et al. [19] evaluated two types of population balance models, namely, the inhomogeneous multi-size group and average bubble number density models, for air-water flows.However, these mechanistic models require various additional and complicated models for the source and sink terms in the transport equations, and a large amount of calculation time in the CFD application.
Meanwhile, previous researchers have presented correlation-type bubble size models that assume bubble size as one group [20][21][22][23].Furthermore, a two-group model consisting of spherical and slug bubbles has been proposed [24][25][26][27].Because these models have been developed on the basis of the average two-phase parameters and properties in the cross-section of the flow channel, they are only applicable to one-dimensional analysis.
In a bubbly flow, we can assume a one-group bubble size.A correlation-type bubble size model provides advantages for practical CFD applications because of its simplicity.In this study, we propose a simple semi-theoretical correlation generated on the basis of a steady-state fully developed bubble number density equation to precisely predict the local bubble diameter for a CFD application.The coefficients of this correlation were determined based on three air-water bubbly flow experiments performed in a vertical pipe under low pressure conditions.For the development of the new correlation, data not provided for the experiments, such as turbulent dissipation rate and relative velocity, were obtained from the CFD calculation.Finally, the proposed correlation and the existing correlations were evaluated against four air-water experiments under low pressure conditions.

Existing One Group Bubble Size Correlations
Hinze's model [28], which employs Kolmogorov's theory, is an empirical correlation that predicts bubble size in a liquid-liquid system.This model includes the liquid properties such as density and surface tension and energy dissipation rate.After that, Akita and Yoshida [20] presented an experimental bubble size correlation using 58 bubble column data under air-water flow.However, this correlation does not properly predict bubble size if it deviates from the bubble size used to develop the correlation [23].The mechanism of bubble breakup can be expressed as a balance between external stress and surface tension, which can be defined as the Weber number.Kocamustafaogullari et al. [21] developed a theoretical bubble size correlation on the basis of the critical Weber number.Because the error increases as the energy dissipation rate decreases [23], it is difficult to apply the correlation to flow conditions with a low energy dissipation rate.Hibiki and Ishii [22] derived the bubble size correlation applicable to bubbly flow from the IATE and validated it with 204 datasets measured under the high convective air-water flow.Subsequently, Hibiki and Ishii [23] developed another bubble size correlation based on 459 data from various fluid conditions.This correlation was validated for a wide range of adiabatic bubbly flows, including various channel sizes, a forced convection flow, and a two-phase flow system in the tank.All the correlations for the gas-liquid system described above were developed using the cross-sectional averaged properties and two-phase parameters of the flow channel under atmospheric pressure conditions.However, when the models were applied to the CFD codes for their evaluation in this study, the calculations were performed by considering local properties and parameters.The correlations with their application conditions are summarized in Table 1.

Two-Phase Flow Model for CFD Calculation
A two-fluid model based on the Eulerian equation and a standard k-ε model were used for the liquid phase [29,30].The turbulent response model of Issa [31] was utilized for the gas phase.Constitutive models for the interfacial drag, turbulence models, wall lift force, wall lubrication force, etc., are required for the Eulerian two-phase flow simulation of bubbly flow.In the present calculation, the drag force proposed by Tomiyama et al. [32], lift force proposed by Auton et al. [33], turbulent dispersion force proposed by Lahey et al. [34], particle-induced turbulence force proposed by Troshko and Hassan [35], wall lubrication force proposed by Antal et al. [36], as well as the virtual mass force proposed by Auton et al. [33] were applied.All the numerical simulations were conducted using the commercial CFD code STAR CCM+ ver 13.06, Munich, Germany [31].

Turbulence Model
The standard k-ε model was used as the turbulent model, and a high-y+ wall treatment (applicable to 30 ≤ y+ ≤ 100) was applied for a continuous liquid phase in this work.This model is widely used in CFD calculations because it exhibits practical applicability for the two-phase flow conditions, in which bubble-induced turbulence plays a dominant role.The transport equations of the turbulent kinetic energy, k, and the turbulent dissipation rate, ε, for a two-phase flow are, respectively, calculated from the following two equations: where µ is the viscosity, ρ C is the density of the continuous liquid phase, S k and S ε are the source terms, and C µ , C ε1 , C ε2 , σ k , and σ ε are constant model coefficients.In a bubbly flow, the liquid turbulence caused by the bubbles is considered as the source term, S k .Rzehak and Krepper [9] evaluated the bubble-induced turbulence (BIT) models proposed by Troshko and Hassan [35], Morel [37], and Politano et al. [38] with different time scales on the radial distribution of the void fraction, liquid velocity, and turbulent kinetic energy for a dispersed gas-liquid flow.Furthermore, the BIT model introduced by Sato et al. [39], including the addition of bubble-induced eddy diffusivity, was evaluated as a reference model.All these models, except for those of Sato et al., showed similar calculation results [9].Therefore, the BIT model proposed by Troshko and Hassan was chosen for the present study.The source term (S k ) and dissipation term (S ε ) of Troshko and Hassan's BIT model applied to the simulation are as follows: where C 3 is the model coefficient, C D is the drag coefficient, ω b is the bubble pseudo turbulence dissipation frequency proposed by Lopez de Bertodano et al. [40], and U r is the relative velocity of the phases.Troshko and Hassan suggested a value of 0.45 for C 3 ; however, this tends to be overestimated at high gas velocities.Therefore, in this study, parametric analysis was performed to obtain the optimum C 3 by changing its value via the Hibiki experiment.When the coefficient is 0.45, the phase velocity is overpredicted, as in Troshko and Hassan, as shown in Figure 1.In this calculation, the coefficient was determined to be 0.2, which exhibited the best prediction with the experimental results.
(a) (b) The turbulent response model was applied as the turbulent model for the gas phase.The model is based on the assumption that the Reynolds stress in the gas phase is directly related to the stress in the liquid phase by the turbulent response coefficient, Ct, which is defined as the ratio of the velocity fluctuation of gas phase to that of the liquid phase as follows: The turbulent eddy viscosity of the gas phase is as follows: The turbulent response model was applied as the turbulent model for the gas phase.The model is based on the assumption that the Reynolds stress in the gas phase is directly related to the stress in the liquid phase by the turbulent response coefficient, C t , which is defined as the ratio of the velocity fluctuation of gas phase to that of the liquid phase as follows: The turbulent eddy viscosity of the gas phase is as follows: where ρ d is the density of the dispersed gas phase.The turbulent kinetic energy (k d ) and turbulent dissipation rate (ε d ) of the gas phase are directly related to those of the liquid phase and are defined below.
Issa clarifies the turbulence response coefficient (C t ) model for bubbly flow [31].The model, derived by considering the Lagrangian equation of motion for a single fluid element encountering a turbulence eddy of the liquid phase, is described using the following equations: ) where A D is the drag correlated coefficient, l e is the eddy length scale, Re t is the turbulent Reynolds number, ν c is the kinematic viscosity of the liquid phase, and U C is the velocity fluctuation of the liquid phase.Issa's turbulence response model is defined as a correlation of C t with a void fraction correction [41].This estimates C t to be up to 3 for a significantly small void fraction of the gas phase but rapidly decreases to unity for dispersed phase volume fractions greater than approximately 6%.

Two-Phase Interaction Models
Drag force refers to a force acting opposite to the relative motion of any object with respect to the surrounding fluid.The drag force ( ) is given by: There are several models available for the drag coefficient.In this study, the drag coefficient model proposed by Tomiyama [32], which covers a wide range of bubbly flow conditions and considers different bubble shapes via the aspect ratio, was adopted as follows: where Re is the Reynolds number, Eo is the Eotvos number, and σ is the surface tension.This model accounts for the deformation effect of a spherical bubble.
The lift force is a component perpendicular to the oncoming flow direction when a fluid flows past the surface of the body.This force was derived by Auton et al. [33] as follows: where U c is the velocity of the liquid phase and C L is the lift coefficient, which has a constant value.
For the bubble size distribution, the coefficient in the present study was adjusted to represent the experimental void fraction distribution.
In two-phase flow, a significantly low void fraction is often observed in the vicinity of the wall.This phenomenon is caused by the force pushing the bubbles from the wall, which is known as the wall lubrication force.This force, as proposed by Antal et al. [36], is expressed as the distance from the wall (y w ), as follows: C WL (y w ) = max −0.01 + 0.05 Here, C WL is a function of the inverse length and is activated only near the wall region.In Equation ( 23), → n is the outward-facing unit normal to the nearest point on the wall such that the force is pointing inwards to prevent bubbles from contacting the wall.
Turbulence and lateral forces play an important role in the lateral movement of bubbles.The effect of turbulence on the redistribution of the phase concentration is modeled by a turbulent dispersion force in the phase momentum equations.When a Reynolds stress tensor is modeled for the liquid phase, the force of Lahey et al. [34] is expressed as: where C µ is the model coefficient of the standard k-ε model, Pr t is the turbulent Prandtl number, k c and ε c are the turbulent kinetic energy and turbulent dissipation rate of the continuous phase, and τ c is the Reynolds stress tensor for the continuous phase.The inertia of the surrounding fluid influences the acceleration of a particle that is submerged in a flow.By including this virtual mass force, the inviscid flow theory can provide a realistic representation of the accelerating flows.This acceleration by the bubble is taken into account through the virtual mass force proposed by Auton et al. [33], provided below: Here, C VM is the virtual mass coefficient depending on the bubble shape and has a value of 0.5 for the individual sperical bubbles, and → a r is the relative acceleration.

Development of a New Bubble Size Model
For the development of a bubble size correlation reflecting the bubble breakup and coalescence mechanisms, the major governing parameters for stable bubble size in the flowing channel are obtained from the one-dimensional bubble number density equation below.
where n, t, υ p , and z are the bubble number density, time, average bubble velocity weighted by the bubble number, and axial position along the flow direction, respectively.φ B and φ PN are the increased rates of bubble number density by bubble breakup and bubble nucleation, respectively.φ C and φ PC are the decreased rates of bubble number density due to bubble coalescence and bubble condensation, respectively.The steady-state fully developed flow condition is defined as a constant bubble number density along the flow direction.Because the bubble number density in the flow direction remains constant, the second term on the left-hand side of Equation ( 27) is zero.Furthermore, the bubble nucleation (φ PN ) and condensation (φ PC ) terms are zero because of the adiabatic flow.Here, one of the mechanisms for increasing the bubble number density is through the turbulent impact between bubbles and turbulent eddies (φ TI ), causing the bubbles to break.Conversely, the mechanism for reducing bubble number density is through collapsing of the bubbles, owing to the random collision between the bubbles (φ RC ).Another bubble coalescence mechanism is wake entrainment (φ WE ), which dominates in small-diameter tubes.When large and small bubbles coexist in the flow channel, the small bubbles are accelerated and merge into large bubbles, owing to the reduced drag in the wake region behind the large bubbles.Therefore, assuming the steady-state fully developed condition of an adiabatic flow, the bubble number density transport equation is simplified as follows: The modeling for each term in Equation ( 28) is summarized in Table 2. Finally, the bubble size of a bubbly flow under steady-state fully developed conditions can be obtained from the modeling of bubble breakup and coalescence as: where α max is the maximum void fraction.The remaining variables (K B , K C , γ B , γ RC , γ WE ) are the constants of bubble breakup and coalescence terms, which are adjustable parameters to be determined experimentally.

Mechanism Model Reference
Turbulent impact Hibiki and Ishii [42] Random collision Hibiki and Ishii [42] Wake entrainment Hibiki et al. [43] However, implementation of the above equation directly into the CFD code is still complicated.Considering the dependence of variables appearing in Equation ( 29), the equation can be simplified as: In an air-water flow under low pressure conditions, the density of the gas is significantly low compared to that of the liquid, thus, (ρ c − ρ d )/ρ c can be assumed to be 1.Subsequently, Equation ( 30) can be non-dimensionalized from the turbulent Reynolds number defined by Hibiki and Ishii [14] (Re b = ε 1/3 Lo 4/3 /ν c ), Laplace length (Lo = σ/g∆ρ), and Weber number defined by relative velocity (We r = 2ρ c U r Lo/σ) as follows: Finally, considering the dependency of variables, we proposed the following bubble size correlation for a bubbly flow under a steady-state fully developed condition by applying Equation (31) with the help of the iterative CFD calculation against the experimental database by following the procedure shown in Figure 2.
The above correlation is a semi-theoretical correlation whose constant and exponents were obtained via multiple linear regression against the experiments by DEDALE [44], Hibiki et al. [45], and Kalkach-Navarro et al. [46].The turbulent dissipation rate and relative velocity, which were not measured in the experiments, were obtained through CFD analysis.
Equation (32) satisfies the precondition that the bubble size should converge to zero when the void fraction tends toward zero.As the amount of void fraction increases, the random collision and the wake entrainment occur actively compared to the bubble breakup mechanism, and bubble size increases, and thus, the exponent of the void fraction is positive.As the turbulent dissipation rate and relative velocity increase, the turbulent impact for bubble breakup increases, and thus, bubble size decreases.Therefore, the exponent of the turbulent Reynolds numbers has a negative value.

Experimental Database for the Simulation
We performed a CFD analysis of four forced convective air-water experiments under lowpressure conditions, as presented in Table 3, to evaluate the bubble size correlations.Table 4 presents the experimental conditions.All experiments were conducted in vertical pipes under the up-flow condition and obtained a radial distribution of local bubble parameters.The bubble shapes of the DEDALE [44], Hibiki [45], and Kalkach-Navarro [46] experiments are estimated as spherical, ellipsoid, and wobbling.However, we assumed a spherical bubble shape in the present study because the three experimental datasets and existing one group IATEs were obtained by applying a spherical The above correlation is a semi-theoretical correlation whose constant and exponents were obtained via multiple linear regression against the experiments by DEDALE [44], Hibiki et al. [45], and Kalkach-Navarro et al. [46].The turbulent dissipation rate and relative velocity, which were not measured in the experiments, were obtained through CFD analysis.
Equation ( 32) satisfies the precondition that the bubble size should converge to zero when the void fraction tends toward zero.As the amount of void fraction increases, the random collision and the wake entrainment occur actively compared to the bubble breakup mechanism, and bubble size increases, and thus, the exponent of the void fraction is positive.As the turbulent dissipation rate and relative velocity increase, the turbulent impact for bubble breakup increases, and thus, bubble size decreases.Therefore, the exponent of the turbulent Reynolds numbers has a negative value.

Experimental Database for the Simulation
We performed a CFD analysis of four forced convective air-water experiments under low-pressure conditions, as presented in Table 3, to evaluate the bubble size correlations.Table 4 presents the experimental conditions.All experiments were conducted in vertical pipes under the up-flow condition and obtained a radial distribution of local bubble parameters.The bubble shapes of the DEDALE [44], Hibiki [45], and Kalkach-Navarro [46] experiments are estimated as spherical, ellipsoid, and wobbling.However, we assumed a spherical bubble shape in the present study because the three experimental datasets and existing one group IATEs were obtained by applying a spherical bubble assumption.In this study, the forced convective air-water bubbly flow is considered to be a steady, adiabatic, and fully developed condition on the measuring plane [47].

Grid Sensitivity Test
The fully developed forced convective bubbly flow in the vertical pipe was axisymmetric along the azimuthal direction; hence, we employed a quasi-2D cylindrical geometry for the simulation.This geometry is a narrow cylindrical sector with a central angle of 10 • , which has symmetry boundary conditions on the side faces.The grid was divided into equal intervals along the axial direction.The flow domain and cross-sectional view of the grid are shown in Figure 3.The size in the axial direction was fixed at 5 mm, and the grid sensitivity test was conducted by changing the number of grids along the radial direction.The grid test was carried out with four meshes.Furthermore, we performed a grid geometry study with the full and half circular shapes for the simulation; however, there was no substantial difference in the calculation results.Here, the physical models introduced in Section 3 were used.Figure 4 compares the results obtained with different meshes in the Hibiki experiments.For the Hibiki experiment, 25 × 612 meshes revealed good independence on mesh and were, therefore, selected for the following calculation.A similar grid test was also performed for other experiments.The grids utilized in each of the calculations are summarized in Table 5.

Evaluation of the New Bubble Size Model
The proposed bubble size correlation was applied to the benchmark calculation of the CFD co nst experiments.

Evaluation of the New Bubble Size Model
The proposed bubble size correlation was applied to the benchmark calculation of the CFD code against experiments.Figures 5-8 show a comparison of the local two-phase flow parameters, such as the void fraction, bubble diameter, and phase velocities predicted with the proposed bubble size correlation and experimental data.Of the four experiments used to evaluate the correlation, the first three experiments listed in Table 5 were used to determine the coefficients of the correlation when developing the correlation; however, VAWL data were not used.The existing correlations investigated in Section 2 were also implemented in the CFD code and evaluated.As illustrated in Figure 5, the numerical calculation result for the DEDALE experiment shows a flat void fraction distribution at the center of the pipe and a rapid increase in the vicinity of the wall with reasonable accuracy.Moreover, the bubble size and phase velocity are well predicted against the experimental data.For the four cases of Hibiki experiments, the peak void fraction in the channel center was well predicted, as depicted in Figure 6a.Additionally, the prediction of the phase velocities shows good agreement, as shown in Figure 6c,d.In Case 3 of the Hibiki experiment, the local void fraction for r/R < 0.7 exceeds 0.25 and, thus, the appearance of a slug bubble is expected [49].Therfore, the relatively large bubbles at the center of the pipe could not be accurately predicted, even though their profiles closely followed the tendency of the bubble size, as shown in Figure 6b.For the Kalkach-Navarro experiment, the data on the velocity of the liquid and gas phase are not provided; only the void fraction and the bubble size are shown in Figure 7. Simulation results show that large bubbles, which appeared in the center in Cases 1-4, were not well predicted, similar to those in the Hibiki experiment.Furthermore, the increase in bubble size near the wall was not predicted in Cases 5-8.Thus, a second group bubble size model for large bubbles is required for better prediction.Nevertheless, the overall value and tendencies of the void fraction and bubble size were predicted with reasonable accuracy.
Finally, the new bubble size model was evaluated for the VAWL experiment, which was not used for determining the correlation coefficients.Figure 8 shows a comparison of the simulation results and experimental data for the local two-phase flow parameters, such as the void fraction, bubble diameter, and phase velocities.The VAWL data show both the core and wall peaks, where bubbles are concentrated at the center and near the wall, respectively, according to the experimental conditions.These phenomena were well predicted by the calculations.In V01 and V03, the void fraction in the nearby wall was overpredicted, and hence, bubble size was also overpredicted near the wall.The velocity of each phase was predicted to be slightly higher at the center.Nevertheless, the evaluation confirms that the proposed bubble size correlation predicts the overall value and tendencies of bubble size and phase velocity with reasonable accuracy.
Figure 9 shows the performances of the existing bubble size correlations and the new bubble size correlation against the four experiments.The Akita and Yoshida correlation as well as the Kocamustafaogullari correlation significantly overpredicted the local bubble size and were unsuitable for the CFD code for a forced convective bubbly flow in a vertical pipe because the two correlations were developed under bubble column and horizontal tube conditions, respectively.Although the two Hibiki and Ishii correlations were developed using the average bubble size data in a flow channel, these correlations yielded acceptable predictions for the experimental data.The new bubble size correlation shows the best performance among all correlations.For quantification of the prediction accuracy of the existing and new correlations for the local bubble diameter, the average error and root-mean-square error (RMSE) are defined as follows: where N is the number of experimental datasets, D b, meas. is the local bubble size obtained through experiments, and D b, calc. is the calculated local bubble size with the correlation.The average error and RMSE are summarized in Table 6.Through the assessment of the new correlation with other correlations, the new correlation has an excellent agreement within an RMSE of ± 11% in the vertical adiabatic bubbly flow under the low pressure condition.

Conclusions
In this study, we developed a semi-theoretical bubble size correlation based on the bubble number density transport equation.The new correlation considers the random collision and wake entrainment for bubble coalescence and turbulent impact for bubble breakup, respectively.To determine the coefficient and exponent of the correlation, we utilized three different experimental databases obtained in a vertical adiabatic bubbly flow under low-pressure conditions.
To validate the performance of the new correlation, three-dimensional numerical simulations for four air-water experiments were carried out using commercial CFD code.The numerical results indicated that the correlation can be used to accurately predict air-water bubbly flow behavior.Through the assessment of the new correlation together with other existing one-dimensional correlations, the new correlation is in excellent agreement with the experimental data.The applicable condition of the new correlation covers forced convective bubbly flows, which are in the ranges of 38.1 mm ≤ hydraulic diameter (D H ) ≤ 80 mm, 0.05 m/s ≤ superficial gas velocity (j g ) ≤ 0.624 m/s, 0.5 m/s ≤ superficial liquid velocity (j f ) ≤ 5 m/s, 1.32 mm ≤ bubble size (D b ) ≤ 5.84 mm, 0.47 ≤ Eo ≤ 2.65, 171 ≤ Re ≤ 8785, and 1 bar ≤ absolute pressure (P) ≤ 3.5 bar.

Figure 1 .
Figure 1.Prediction of radial distribution of (a) void fraction, and (b) liquid velocity by the different coefficients of Troshko and Hassan's BIT model.

Figure 1 .
Figure 1.Prediction of radial distribution of (a) void fraction, and (b) liquid velocity by the different coefficients of Troshko and Hassan's BIT model.

Figure 2 .
Figure 2. Iterative calculation procedure for the development of the bubble size model.

Figure 2 .
Figure 2. Iterative calculation procedure for the development of the bubble size model.

Figure 3 .
Figure 3. (a) Flow domain with boundary conditions, (b) grid near the outlet, and (c) cross-section of the grid.

Figure 3 . 4 Figure 4 .
Figure 3. (a) Flow domain with boundary conditions, (b) grid near the outlet, and (c) cross-section of the grid.

Figures 5 -Figure 4 .
Figure 4. Results of grid sensitivity study for the Hibiki experiment.

Table 1 .
Existing one group bubble size correlations for the gas-liquid system.

Table 2 .
Models for bubble breakup and coalescence.

Table 3 .
Air-water bubbly flow experiments used to evaluate the bubble size correlation.

Table 4 .
Experimental conditions of each experiment used to simulation.

Table 5 .
Grids used in the benchmark calculation of each experiment.

Table 5 .
Grids used in the benchmark calculation of each experiment.

Table 6 .
Average error and RMSE of the bubble size correlation for each experiment.