## 1 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. [PERSON] and [PERSON] [1] applied an interfacial area transport equation (IATE) proposed by [PERSON] 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, [PERSON] and [PERSON] [3] and [PERSON] and [PERSON] [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]. [PERSON] 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. [PERSON] 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. [PERSON] 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. \[a_{i}=\frac{6\alpha}{D_{b}} \tag{1}\]  
where \(a_{i}\) is the IAC, \(D_{b}\) is the bubble size, and \(\alpha\) is the void fraction. Many mechanistic bubble size models have been proposed to predict bubble size (or IAC) in a flow channel accurately. [PERSON] 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. [PERSON] et al.