Effect of interfacial drag force model on code prediction for upward adiabatic two-phase bubbly flow in vertical channels

Accurate modeling of the interfacial drag force is one of the keys to predicting thermo-fluid parameters using one-dimensional nuclear thermal-hydraulic system analysis code architected through the two-fluid model. The interfacial drag force appears in the interfacial momentum transfer term and governs the velocity slip or the relative velocity between gas and liquid phases. The most straightforward method to model the interfacial drag force is to model the force through the drag law (drag law approach). A drag coefficient and interfacial area concentration should be given to close the interfacial drag force model. Among them, the modeling of the interfacial area concentration has been one of the weakest links in the interfacial drag force modeling due to the lack of reliable data covering a wide test condition including prototypic nuclear reactor conditions and lack of physically sound interfacial area model. To avoid a considerable uncertainty in the prediction of the interfacial area concentration, Andersen and Chu (1982) proposed the interfacial drag force model using the drift-flux parameters (Andersen approach). The Andersen approach is practical for the simulation of a slow transient flow and a steady flow. Major system analysis codes such as USNRC TRACE have adopted the Andersen approach in the interfacial drag force modeling. Some attempts to improve the code performance have been considered using the drag law approach with the interfacial area transport equation. The dynamic modeling of the interfacial area concentration has the potential to improve the prediction accuracy of the interfacial area concentration in a transient flow and developing flow. Due to the importance of the improved interfacial drag force modeling, the implementation and evaluation of the interfacial area transport equation in USNRC TRACE code has been performed by Talley et al. (2011, 2013). The study claimed that the introduction of the interfacial area transport equation into the TRACE code improved the code performance in an adiabatic bubbly flow analysis significantly. The present study assessed the code calculation made by Talley et al. (2011) and identified several issues in the code calculation results. The present study analytically demonstrated that the drag law approach became identical with the Andersen approach for the distorted particle regime (or a major bubble shape regime in bubbly flow) due to the balancing-out of the interfacial area concentration (or bubble size) in the numerator and denominator of the interfacial drag force formulation. The code calculation using TRAC code endorsed the analytical assessment of the insignificant or no merit of the interfacial area transport equation in the code performance of the adiabatic bubbly flow analysis. The present study also pointed out the inconsistency of the code calculation made by Talley et al. (2011).

considered one of the most effective two-phase flow formulations (Ishii and Hibiki, 2011). The two-fluid model is composed of six equations, i.e., mass, momentum, and energy conservation equations for gas and liquid phases. Since the two-fluid model can treat gas and liquid phases separately, kinematic and thermal non-equilibrium between gas and liquid phases can be simulated using the two-fluid model with proper interfacial transfer terms.
Most of major one-dimensional nuclear thermal-hydraulic system analysis codes such as TRACE, RELAP5, and TRAC-BF1 codes (Borkowski et al., 1992;Ransom et al., 2001;NRC US, 2008) have adopted the two-fluid model to architect the simulation code. One of the key models in the code is the interfacial drag force, which governs the kinematic non-equilibrium between gas and liquid phases or velocity slip, and the accurate modeling of the interfacial drag force is indispensable to predict void fraction accurately. The most straightforward method of the interfacial drag force modeling is to formulate the interfacial drag force using a drag law. In this approach (hereafter, drag law approach), the interfacial drag force is represented by a function of the interfacial area concentration, drag coefficient depending on bubble shape regime, liquid density and relative velocity between two phases. The area-averaged value of the relative velocity is formulated with the aid of the distribution parameter in the drift-flux model, and reliable drag coefficient model is available (Ishii and Hibiki, 2011). However, the accurate modeling of the interfacial area concentration has been one of the weakest links in the interfacial drag force modeling. Due to its importance of the interfacial area concentration modeling in the two-phase flow simulation, analytical, experimental, and computational studies have been performed and some review articles are available (Hibiki and Ishii, 2009;Lin and Hibiki, 2014;Chuang and Hibiki, 2015).
Due to the lack of a reliable and accurate model of the interfacial area concentration during the code development period, alternative modeling method of the interfacial drag force was proposed to avoid a significant uncertainty in the prediction of the interfacial area concentration. Andersen and Chu (1981) formulated the interfacial drag force with the aid of the drift velocity in the drift-flux model (hereafter, Andersen approach or drift-velocity approach). The Andersen approach is practical for a slow transient flow and steady flow. Major one-dimensional nuclear thermal-hydraulic system analysis codes have adopted the Andersen approach to calculating the interfacial drag force. Mortensen (1995) and Kelly (1997) identified several unsolved issues in the current modeling approach based on "static" flow regime transition criteria such as (1) lack of dynamic predictive capability in the current flow regime transition criteria, (2) compound errors due to two-step methods (errors in constitutive equations and flow regime transition criteria), and (3) insufficient validation of existing flow regime dependent correlations and criteria under prototypic nuclear reactor conditions. In order to partly solve the issues, the introduction of the interfacial area transport equation has been considered in the middle of the 1990s.
The interfacial area transport equation is considered as a simplified population balance approach often used in chemical engineering. Considering a computational load for additional gas momentum equations, a two-group approach to treat bubbles in two groups is adopted. The group-one bubbles are spherical and distorted bubbles, whereas the group-two bubbles are cap, slug, and churn-turbulent bubbles. Since the characteristic bubble transport velocity is different between the group-one and group-two bubbles, two gas momentum equations for the group-one and group-two bubbles should be solved in a code. Two-group interfacial area transport equation is reduced to one-group interfacial area transport equation in bubbly flow regime. Some review articles are available for the progress of the interfacial area transport equation development (Hibiki and Ishii, 2009;Ishii and Hibiki, 2011). In the current development of the interfacial area transport equation, some critical defects exist. The prediction accuracy of the interfacial area transport equation is susceptible to the initial (or inlet) conditions. Unless correct initial conditions are specified, the interfacial area transport equation cannot predict the interfacial area concentration accurately. Most of the current interfacial area transport equation has not been validated under prototypic nuclear reactor conditions, and its applicability to the prototypic conditions has not been confirmed. Many source and sink terms of the interfacial area transport with complicated functional forms and many empirical constants hamper the validation of each term. Additional gas momentum equation may cause unexpected numerical instability in the code calculation.
Limited study has been performed to implement the interfacial area transport equation in a nuclear thermalhydraulic system analysis code and to evaluate the effect of the interfacial area transport equation on the code performance. Talley et al. (2011Talley et al. ( , 2013 implemented one-group interfacial area transport equation into USNRC TRACE code and claimed that "the inclusion of the one-group interfacial area transport equation in TRACE yields significant improvements to prediction results when compared with those made by the flow regime based algebraic closure relations." However, a brief review of their code calculation results identifies some critical defects and inconsistencies in their modeling. In view of these, the present study is aiming at the revisit of the implementation and evaluation study made by Talley et al. (2011). The present study will identify the role of the interfacial area concentration in adiabatic bubbly two-phase flow simulation. T. Hibiki,T. Ozaki 214 2 Assessment of TRACE code calculation done by Talley et al. Talley et al. (2011) implemented a one-group interfacial area transport equation into the USNRC TRACE code and evaluated the performance of the code modified with the interfacial area transport equation. Talley et al. (2011) utilized TRACE 4.291b architected based on the three-field two-fluid model, namely two gas momentum equations for the twogroup interfacial area transport equation and one liquid momentum equation. However, TRACE 4.291b applied only to one-group case, namely bubbly flow in the work done by Talley et al. (2011). Talley et al. (2011) gave the gas momentum equation as

Implementation
and θ are the void fraction weighted velocity of the group-one bubble, time, one-dimensional flow direction, interfacial drag coefficient between continuous liquid and dispersed gas, void fraction of the group-one bubble, gas density, liquid velocity, pressure, gravitational acceleration, and inclination angle from vertical, respectively. It should be noted here that Eq. (1) is the gas momentum equation given in the original paper by Talley et al. (2011). It seems that Talley et al. (2011) misused a vectorial notation in Eq. (1) for the one-dimensional twophase flow formulation. In what follows, the equations that appeared in the original paper will be given, but the subscript of "1" for the group-one bubble is dropped for simplicity. Talley et al. (2011) showed the interfacial drag "coefficient" used in the conventional TRACE code as where g j v  , Δρ , 0 C , and s P are the drift velocity, density difference between two phases, distribution parameter, and profile slip factor, respectively. The interfacial drag in the conventional TRACE code is formulated based on Andersen approach (Andersen and Chu, 1981).
According to Talley et al. (2011), "in the revised implementation, this is replaced using an interfacial drag coefficient for the group one bubbles to be more consistent with the transport equation approach. Hence, the interfacial drag is given by" where D C , f ρ , and Sm D are the drag coefficient of the groupone bubble, liquid density, and Sauter mean diameter of the group-one bubble, respectively. It should be noted that 3/4 and the profile slip factor are missing in the right-hand side of Eq. (3) in the original paper (Talley et al., 2011). The correct form of Eq. (3) is as follows: Re C Re where the mixture Reynolds number is given by In Eq. (6), f μ is the liquid viscosity. It should be noted that Eq. (5) is the drag coefficient for the viscous regime (or low mixture Reynolds number).
The applicability of Eq. (5) is limited by the following maximum spherical bubble diameter or threshold bubble diameter between the viscous and distorted particle regimes as (Hibiki and Ishii, 2002a): where σ is the surface tension and the viscous number is defined by When the bubble diameter exceeds max1 D , the drag coefficient for the distorted particle regime should be utilized. When the bubble diameter further increases and exceeds max 2 D , the drag coefficient for the cap bubble regime should be utilized. The threshold diameter between the distorted particle and cap bubble regimes, max 2 D , is given by The drag coefficient for the distorted particle regime is given as (Ishii and Hibiki, 2011): Effect of interfacial drag force model on code prediction for upward adiabatic two-phase bubbly flow in vertical channels

215
where Re ¥ is the Reynolds number for a single particle. The threshold diameters, max1 D and max 2 D , for air-water flow at atmospheric pressure and 20 °C are estimated to be 2 mm and 10 mm, respectively. Talley et al. (2011) compared the calculated interfacial area concentration with the measured interfacial area concentration. " Figure 5" in the paper by Talley et al. (2011) showed the comparison for Run No. 2-7. Talley et al. (2011) indicated that "the experimentally measured size for this condition is only 2.4 mm". Talley et al. (2011) misused Eq. (5) in calculating the drag coefficient for the distorted particle regime and should have used Eq. (10) in their code calculation for Run No. 2-7.
The conventional TRACE code utilizes the following relationship in calculating the interfacial area concentration in the bubbly flow regime. Since the conventional TRACE code adopts the Andersen approach to calculate the interfacial drag force, the correlation of the interfacial area concentration is used in calculating the interfacial heat transfer.
i Sm where the bubble Sauter mean diameter is calculated by  (4) and (2) were utilized in TRACE-T and TRACE-NT codes, respectively.

Evaluation
In the evaluation process, Talley et al. (2011) utilized two datasets taken for atmospheric air-water bubbly flow in round pipes with the diameters of 25.4 mm and 48.3 mm (Hibiki and Ishii, 1999;Fu, 2001). The flow conditions are listed in Table 1. Some inconsistency of the flow condition in Run No. 2-7 is identified. The flow condition labeled by Run No. 2-7 ( g j =0.138 m/s and f j =2.49 m/s) does not exist in Fu (2001). The flow condition closest to Run No. 2-7 is g j =0.138 m/s and f j =2.34 m/s. In the present study, the superficial gas velocity, g j , and superficial liquid velocity,  Talley et al. (2011). In the figures, the predictions made by TRACE-T and TRACE-NT are indicated by solid and broken lines, respectively, and measured data are represented by an open symbol. The error bar indicates the ±10% range. Figure 1(a) shows that the gas velocity predictions made by TRACE-NT (Andersen approach) tend to overestimate the experimental data compared to those made by the TRACE-T (drag law approach using the interfacial area transport equation) which yields better results. However, the results shown in Fig. 1(a) are inconsistent with the results shown in Figs. 2(a) and 3(a). As shown in Fig. 2(a), the predicted local pressure using TRACE-NT is the same as that using TRACE-T, which means that the predicted local superficial gas velocity using TRACE-NT should be the same as that using TRACE-T. Figure 3(a) indicates that the predicted local area-averaged void fraction using TRACE-NT is similar to that using TRACE-T. The relationship among the gas velocity, superficial gas velocity, and void fraction is given by The computational results shown in Figs. 2(a) and 3(a) indicate that the predicted local gas velocity using TRACE-NT should be similar to that using TRACE-T for Run No. 2-5, which contradicts the results shown in Fig. 1(a). Figure 4(a) compares the predicted interfacial area concentrations with the measured data. Talley et al. (2011) claimed that " i a predictions made by TRACE-T agree very well with the experimental data, within ±10% difference, while those made by TRACE-NT can show significant disagreement." However, Talley et al. (2011) also recognized that "(t)he main difference observed from Fig. 4(a) (Fig. 5(a) in the original paper), however, is that TRACE-T and TRACE-NT do not have the same initial value. For TRACE-T the initial experimental value is specified, whereas for TRACE-NT the initial value is calculated based on Eqs. (11) and (12) (Eqs. (26) and (27) in the original paper). Using these relations, the initial bubble size of TRACE-NT is approximately 5.4 mm, while the experimentally measured size for this condition is only 2.4 mm" and "if an improved bubble size relation is employed in TRACE-NT, a prediction similar to TRACE-T can be obtained." The argument by Talley et al. (2011) indicates that the interfacial area transport equation may not provide accurate prediction unless the initial bubble size is accurately given. This is one of the weakest links in the approach using the interfacial area transport equation (Chuang and Hibiki, 2015;Hibiki et al., 2018).
As discussed in Section 2.1, since the bubble size measured for Run No. 2-7 is between max1 D and max 2 D , the drag coefficient for the distorted particle regime, Eq. (10), should have been used in TRACE-T calculation instead of Eq. (5). If the correct drag coefficient, Eq. (10), is used, the drag coefficient for TRACE-T, Eqs. (4) and (10), becomes identical with the drag coefficient for TRACE-NT, Eq. (2). This indicates that the drag coefficient is independent on the bubble size or the interfacial area concentration. In other words, the implementation of the one-group interfacial area transport equation has no impact on the code prediction in adiabatic bubbly flow, which is completely different from the conclusions given by Talley et al. (2011). The detailed derivation is given below.
As demonstrated in Fig. 5, the functional value of ( ) f α ¢ can be approximated to be unity in the void fraction range between 0 and 0.3. Equation (17) is reduced to Eq. (2) resulting in Talley et al. (2011) emphasized the lack of modeling dynamic bubble interactions in TRACE-NT using the comparison between calculated and measured interfacial area concentrations for Run No. 1-11 as shown in Fig. 6(a). Talley et al. (2011) claimed that "the prediction made by TRACE-NT is linear, while the data show a non-linear increase in i a " and "(s)uch a trend in i a development can occur in practice when the bubble break-up due to turbulence impact." However, the non-linear increase in i a is insignificant, and the "non-linear" increase identified by Talley et al. (2011) can be approximated by the linear increase within the measurement uncertainty. The argument by Talley et al. (2011) over-magnifies the importance of the interfacial area transport equation for such an adiabatic quasi-fully developed two-phase bubbly flow. Talley et al. (2011) also emphasized the necessity of the two-group approach using the comparison between calculated  T. Hibiki,T. Ozaki 218 and measured interfacial area concentration for Run No. 1-12 as shown in Fig. 7(a). According to Talley et al. (2011), "the area-averaged void fraction at the third measurement port is greater than 0.15, implying the presence of group two bubbles" and "the presence of group two bubbles in this condition is evidenced by the sudden reduction in i a at the final measurement port." They concluded that "in order for TRACE-T to be applicable to comprehensive dispersed two-phase flow conditions, implementation of the two-group interfacial area transport equation is necessary." The percentage of the void fraction for the two-group bubbles is about 15% of the total void fraction at the third measurement port (Hibiki and Ishii, 1999). The reduction of the interfacial area concentration at the third measurement port up to 15% can be explained by the presence of the group-two bubbles. However, the deviation of the interfacial area concentration calculated by TRACE-T from the measured interfacial area concentration more than 15% is attributed to the insufficient prediction capability of the one-group interfacial area transport equation used in the TRACE-T calculation (Wu et al., 1998;Kim, 1999).

Identified issues
Talley et al. (2011) claimed that "the inclusion of the one-group interfacial area transport equation in TRACE yields significant improvements to prediction results when compared with those made by the flow regime based algebraic closure relations." However, as discussed in Sections 2.1 and 2.2, some issues in their implementation study have been identified as follows: 1) The predictions of the local pressure (or local superficial gas velocity) and local void fraction made by TRACE-T and TRACE-NT are nearly identical and agree very well with the experimental data. Nevertheless, some differences in the predicted local gas velocity between the two codes are observed. This is the violation of the identity between gas velocity and superficial gas velocity divided by void fraction if the same boundary conditions are used for TRACE-T and TRACE-NT calculations.
2) The drag coefficient for the viscous regime is given in the implementation study made by Talley et al. (2011). The dependence of the bubble size (or interfacial area  Effect of interfacial drag force model on code prediction for upward adiabatic two-phase bubbly flow in vertical channels 219 concentration) on the interfacial drag force should disappear in the distorted particle regime. The implementation of the one-group interfacial area transport equation into TRACE code should not improve the prediction results for the adiabatic two-phase flow in the distorted particle regime.
3) The superiority of the one-group interfacial area transport equation to the algebraic correlation is due to a priori specified initial bubble size in the one-group interfacial area transport equation. The implementation study made by Talley et al. (2011) overemphasizes the code improvement with the implementation of the one-group interfacial area transport equation in the adiabatic bubbly flow.
In addition to the above issues, 4) In the code calculation, the following correlations of the distribution parameter and drift velocity are utilized.
1 4 Equation (20) can predict the distribution parameter for slug and churn flow regimes where the void fraction distribution has a core peak. However, Eq. (20) overestimates the distribution parameter for the adiabatic bubbly flow regime where the void fraction distribution has a wall peak (Hibiki and Ishii, 2002a). An accurate correlation to predict the distribution parameter in the bubbly flow regime (Hibiki and Ishii, 2002a) is recommended for a code calculation in a future study.
5) The effect of the covariance due to void fraction distribution on the interfacial drag force has been identified (Hibiki and Ozaki, 2017). The inclusion of the covariance effect on Eqs. (2), (4), and (19) is recommended for a code calculation in a future study.
In what follows, a replication study of the code improvement using the one-group interfacial area transport equation will be conducted to address the issues 1)-3).

Implementation
In the present implementation study, TRAC-BF1 code (Borkowski et al., 1992) is utilized due to the unavailability of the TRACE code. The one-dimensional interfacial drag force, ig M , is given by where i C is the drag force coefficient. In the distorted particle regime, the drag force coefficient is given by Eq. (23) using Eq. (10) .
The relative velocity is calculated by Equation (20) is utilized for calculating the distribution parameter.
It should be noted here that the dependence of the drag force coefficient on the interfacial area concentration disappears in the distorted particle regime, see Eq. (23). If the drag force coefficient is to be formulated in the viscous regime, the drag coefficient given by Eq. (5) Here, H D and ε are the hydraulic equivalent diameter of a flow channel and energy dissipation rate per unit mass, respectively. Hibiki-Ishii correlation (2002b), Eq. (26), has been validated by extensive database taken in adiabatic bubble columns and forced convective bubbly flows under fully developed conditions. The databases used for the model development include 459 datasets taken in various loop and flow conditions such as channel geometry (circular and rectangular channels), hydraulic equivalent diameter (9.0-5500 mm), flow direction (vertical and horizontal directions), superficial gas velocity (0.000788-4.87 m/s), and superficial liquid velocity (0-6.55 m/s). The fluid systems in the databases cover an extensive range of physical properties such as liquid density (684-1594 kg/m 3 ), liquid viscosity (0.410-21.1 mPa·s), and surface tension (20.0-75.0 mN/m). The prediction accuracy of Hibiki-Ishii correlation is estimated to be ±22.0%. If 30% prediction accuracy is accepted, Hibiki-Ishii correlation applies to a developing flow at H z D larger than 10. It is expected that Hibiki-Ishii correlation predicts the interfacial area concentration much accurate than Eqs. (11) and (12).
The drag force coefficient in Andersen approach is given by Equation (16) is utilized for calculating the drift velocity.
In what follows, two code versions, TRAC-CD and TRAC-VGJ are utilized in the code calculation. TRAC-CD code calculates the interfacial drag force based on drag law approach with Eqs. (4) and (23). TRAC-VGJ code calculates the interfacial drag force based on Andersen approach with Eqs. (2) and (31). Other than the interfacial drag force formulation, the constitutive equations used in two codes are identical.

Evaluation
The code calculations using TRAC-CD and TRAC-VGJ codes are performed using the flow conditions the same as those shown in Figs. 1(a), 2(a), 3(a), 4(a), 6(a), and 7(a). The corresponding calculation results are shown in Figs. 1(b), 2(b), 3(b), 4(b), 6(b), and 7(b), respectively. In the figures, the predictions made by TRAC-CD and TRAC-VGJ are indicated by solid and broken lines, respectively, and measured data are represented by an open symbol. The error bar indicates the ±10% range.
As compared with the figures (a) and (b), the predictions made by TRAC-VGJ in the figure (b) are almost identical with those made by TRACE-NT in the figure (a) due to similar constitutive equations used in the two codes. The excellent agreement in the predictions between TRAC-VGJ and TRACE-NT indicates that the two codes have similar prediction accuracy. Figure 2(b) compares the local pressures predicted by TRAC-CD with the measured local pressures. The TRAC-CD code prediction shows an excellent prediction capability of the local pressure. The local pressures predicted by TRAC-CD also agree with those predicted by TRACE-T, TRACE-NT, and TRAC-VGJ. Figure 3(b) compares the local void fractions predicted by TRAC-CD with the measured local void fractions. The local void fractions predicted by TRAC-CD code agree with the measured local void fractions reasonably well. The local pressures predicted by TRAC-CD also agree with these predicted by TRACE-T, TRACE-NT, and TRAC-VGJ. As can be seen from Figs. 2(a), 2(b), 3(a), and 3(b), the code predictions for local pressure and void fraction made by four codes, TRAC-CD, TRAC-VGJ, TRACE-T, and TRACE-NT, are nearly identical. Figure 1(b) indicates that the gas velocity predictions made by TRAC-CD and TRAC-VGJ are nearly identical, which is consistent with the results shown in Figs. 2(b) and  3(b). The similar predictions of the local pressure (or local superficial gas velocity) and the local void fraction should result in the similar prediction of the local gas velocity. However, the gas velocity predictions made by TRACE-T is much lower than those made by other three codes, TRAC-CD, TRAC-VGJ, and TRACE-NT in spite of the similar predictions of the local pressure and void fraction made by all four codes. The local gas velocity predictions made by TRACE-T are not consistent with those made by TRAC-CD, TRAC-VGJ, and TRACE-NT. As discussed in Section 2.2, if the correct drag coefficient is used, the gas velocity predictions based on drag law approach (TRAC-CD and TRACE-T) should be the same as those based on Andersen approach (TRAC-VGJ and TRACE-NT). The interfacial area concentration should not appear in the formulation of the interfacial drag force, and thus the code prediction results of the gas velocity should be independent on the interfacial area concentration. It is recommended that the implementation of the constitutive equations in TRACE-T should be verified and the consistency of the boundary conditions between TRACE-T and TRACE-NT should be carefully checked. Figure 4(b) compares the local interfacial area concentrations predicted by Eq. (26), with the measured local interfacial area concentrations. The predictions labeled by TRAC-CD and TRAC-VGJ are the interfacial area concentration predicted by Hibiki-Ishii correlation with the input of the flow parameters predicted by TRAC-CD and TRAC-VGJ, respectively. Figure 4 (b) shows that the predicted interfacial area concentrations are lower than the measured interfacial area concentration of the order of 15%. This discrepancy may be due to the prediction accuracy of Hibiki-Ishii correlation, ±22.0%, or insufficient accuracy of the input flow parameters predicted by TRAC-CD and TRAC-VGJ. Considering the prediction accuracy of Hibiki-Ishii correlation, ±22.0%, it can be concluded that the predicted interfacial area concentration agrees with the measured interfacial area concentration within its prediction accuracy of the correlation. Figure 6(b) shows an excellent agreement between the predicted and measured interfacial area concentrations. The prediction accuracy of Hibiki-Ishii correlation is similar to that of the one-group interfacial area transport equation with a given inlet condition. Figure 7 (b) shows that the predicted interfacial area concentrations agree with the measured interfacial area concentrations at the second and third measurement ports where the bubbly flow is fully developed. The prediction accuracy of Hibiki-Ishii correlation is deteriorated at the first measurement port because the prediction accuracy of Hibiki-Ishii correlation is deteriorated for a developing flow.

Sensitivity analysis
In what follows, some sensitivity analyses will be performed to identify the critical parameter which affects the code calculation. In the sensitivity analyses, gas and liquid momentum equations for steady-state fully developed twophase flow are considered. Eliminating the pressure gradient term in the gas and liquid momentum equations yields: Equation (32) indicates that the buoyancy force is balanced with the interfacial drag force given by Eq. (22). In the following discussion, the right-hand-side of Eq. (32) is expressed by b M , whereas the right-hand-side of Eq. (22) is expressed by ig M . Figure 8 shows the dependence of calculated b M and ig M on the void fraction. In Fig. 8(a), the solid line indicates the calculated b M , and the broken, chain, and chain double-dashed lines indicate the calculated ig M with 0.5 i C , 1.0 i C , and 1.5 i C , respectively. The sensitivity analysis on ig M with varied i C practically means the sensitivity analysis on ig M with varied i a . In Fig. 8(b), the solid line indicates the calculated b M , and the broken, chain and chain double-dashed lines indicate the calculated ig M with C ¥ =1.2, 1.05, and 0.90, respectively. Here, C ¥ is the asymptotic value of the distribution parameter which appears in the following equation.
( ) The sensitivity analysis on ig M with varied C ¥ (or 0 C ) practically means the sensitivity analysis on ig M with varied r v , see Eq. (25). The void fraction satisfying the gas and liquid momentum equations is the one at the intersection point of ig M and b M . Figure 8(a) indicates that the characteristic curves of ig M are insensitive to the value of i C for the flow condition at the third measurement port in Run No. 2-5. Figure 8(b) shows that the value of the distribution parameter affects the characteristic curve of ig M . In conclusions, the drag force coefficient including the interfacial area concentration does not affect the calculated void fraction, whereas the distribution parameter used for calculating the relative velocity affects the calculated void fraction. This conclusion is consistent with that obtained at the void fraction of about 0.5 for steam-water flow under the pressure of 7 MPa . Figure 9 shows the sensitivity analysis for the lower liquid velocity case . The result similar to the sensitivity analysis for Run No. 2-5 is obtained. Figures 10(a), 10(b), and 10(c) show the sensitivity analyses for local gas velocity, pressure, and void fraction with varied distribution parameter, respectively. Table 2 shows the flow parameters at three measurement ports in Run No. 2-5. The distribution parameter is estimated by As can be seen in Table 2, the distribution parameter changes along the flow direction, and the flatter void fraction profile becomes the void fraction profile with more enhanced core peak as the flow is developed. As can be seen in Fig. 10(c), the change of the asymptotic distribution parameter from 0.90 to 1.2 reduces the calculated void fraction of 0.02. In other words, the change of the distribution parameter of the order of 0.1 changes the void fraction of 0.01 in Run No. 2-5. Figure 10(b) indicates that the calculated pressure    (or calculated superficial gas velocity) is insensitive to the asymptotic distribution parameter. Figure 10(a) shows the change of the asymptotic distribution parameter from 0.90 to 1.2 increases the gas velocity from 5.2 to 6.8 m/s, which corresponds to about 30% increase. The distribution parameters measured at the first and second measurement ports are 1.08 and 1.03, respectively, and agree with the calculation with the asymptotic distribution parameter of 1.05. The distribution parameter measured at the third measurement port is 1.19 and agrees with the calculation with the asymptotic distribution parameter of 1.20. Talley et al. (2011) attributes the discrepancy in the gas velocity between the calculation and measurement to the estimated interfacial area concentration. The above sensitivity analyses indicate that the calculated gas velocity is mainly dependent on the distribution parameter.

Conclusions
The implementation of the interfacial area transport equation into a one-dimensional nuclear thermal-hydraulic system analysis code is being considered to improve the code prediction accuracy. Talley et al. (2011) claimed that the introduction of the one-group interfacial area transport equation into USNRC TRACE code improved the prediction accuracy significantly. The present study first assessed the code calculation made by Talley et al. (2011) and identified the following issues in the code calculation made by Talley et al. (2011).
1) The predictions of the local pressure (or local superficial gas velocity) and local void fraction made by TRACE-T and TRACE-NT were nearly identical and agreed very well with the experimental data. Nevertheless, some differences in the predicted local gas velocity between the two codes were observed. This was the violation of the identity between gas velocity and superficial gas velocity divided by void fraction if the same boundary conditions were used for TRACE-T and TRACE-NT calculations.
2) The drag coefficient for the viscous regime was given in the implementation study made by Talley et al. (2011). The dependence of the bubble size (or interfacial area concentration) on the interfacial drag force should disappear in the distorted particle regime. The implementation of the one-group interfacial area transport equation into TRACE code should not improve the prediction results for the adiabatic two-phase flow in the distorted particle regime.
3) The superiority of the one-group interfacial area transport equation to the algebraic correlation was due to a priori specified initial bubble size in the one-group interfacial area transport equation. The implementation study made by Talley et al. (2011) overemphasized the code improvement with the implementation of the one-group interfacial area transport equation in the adiabatic bubbly flow.
The present study analytically demonstrated that the functional form of the interfacial drag force formulated based on the drag coefficient (drag law approach) became identical with that formulated based on the drift velocity (Andersen approach or drift velocity approach) in the distorted particle regime. The present study developed two code versions such as TRAC-CD based on drag law approach and TRAC-VGJ based on Andersen approach. As expected, two codes' predictions were nearly identical in the adiabatic bubbly flow, and the interfacial area concentration did not have any effects on the code predictions in the bubbly flow regime.
The present study performed the sensitivity analysis for the interfacial drag force with varied drag force coefficient including the interfacial area concentration and the distribution parameter in the relative velocity. The sensitivity analysis demonstrated that the interfacial area concentration was insensitive to the calculated gas velocity and the distribution parameter affected the calculated gas velocity.
The above conclusions do not mean that the interfacial area concentration is not critical in the analysis of the twophase bubbly flow with phase change. An accurate prediction of the bubble interfacial area concentration is crucial for predicting heat transfer between two phases accurately. It is worth mentioning that a study to implement two-group interfacial area correlation into a one-dimensional nuclear thermal-hydraulic system analysis code has been initiated  because some fairly reliable two-group interfacial area correlations are available (Hibiki and Ishii, 2002b;Hibiki et al., 2006;Ozar et al., 2012;Schlegel and Hibiki, 2015;Shen and Hibiki, 2015).