CFD Analysis of Subcooled Flow Boiling in 4 × 4 Rod Bundle

: Rod bundle ﬂow is an important research ﬁeld related to reactor cooling in nuclear power plants. Owing to the rapid development of computerized performance assessments, interest in coolant ﬂow analysis using computational ﬂuid dynamics has garnered research interest. Rod bundle ﬂow research data compared with experimental results under various conditions are thus needed. To address this, a boiling model veriﬁcation study was conducted with reference to experiments. This study adopts the Reynolds-averaged Navier–Stokes equation, a practical analysis method compared to direct numerical simulation and large eddy simulation, including turbulence modeling, to predict the ﬂow of coolant inside a rod bundle. This study also investigates void behavior in low-pressure subcooled ﬂow boiling using a Eulerian approach (two-ﬂuid model). The rod bundle has a length of 0.59 m and a hydraulic diameter of approximately 14.01 mm. At the cross-section at a height of 0.58 m, near the exit, numerical results were compared with the experimental values of the volume fraction of vapor and interfacial area concentration. The simulation results showed good agreement with the experimental data for six di ﬀ erent initial conditions with constant density. W.-G.P.; formal analysis, Y.-B.S. and S.S.P.; funding acquisition, W.-G.P.; investigation, Y.-B.S. and S.S.P.; methodology, W.-G.P.; project administration, W.-G.P. and B.Y.; Supervision, W.-G.P. and B.Y.; writing—original draft, Y.-B.S. and J.-Y.B.; writing—review & editing, Y.-B.S., W.-G.P. and B.Y.;


Introduction
Computational fluid dynamics (CFD) is being increasingly used in the analysis of three-dimensional two-phase flow for the safe operation of a nuclear power plant. In recent years, subcooled boiling in the rod bundle is one of the important phenomena observed in small break loss of coolant accidents (SBLOCA) for pressurized water reactor (PWR) and is also expected in the loss of heat sink accident of nuclear spent fuel pool (SFP) as a safety issue. Thus, experiments and measurements on multiphase turbulent flow in pipes or rod bundles have been conducted, and numerous fluid properties have been studied. In recent years, various simulations have been conducted in various papers [1][2][3][4][5]. Numerical analysis of subcooled two-phase flow at high pressure has confirmed that the Eulerian multiphase model including a Reynolds stress turbulence model is suitable for predicting the actual phenomenon [6][7][8][9][10]. Single-phase convective heat transfer and vapor-water two-phase boiling heat transfer experiments have been conducted in a vertical seven-rod bundle at a low flow rate near atmospheric pressure. In this case, bubble generation is dominantly determined by the enhancement factor of forced convection by the inlet mass flux and the suppression factor of nucleate boiling. Thus, it was necessary to react sensitively to the fluid velocity, and in the next step, excessive bubble prediction is prevented through a nuclear boiling suppression factor. By taking in account the enhancement factor of forced convection and the suppression factor of nucleate boiling, a new correlation in terms of boiling number, Reynolds number, and Martinelli number was found to be effective in predicting two-phase-flow boiling heat transfer coefficients [11]. Analytical studies using various turbulence models for rod bundles in various shapes have also been conducted. A 3 × 3 rod bundle flow was analyzed using the k-ε model, shear stress transport k-ω model, and Reynolds stress model; and the pressure rate, flow velocity, and temperature were predicted [12]. For supercritical pressure conditions, the importance of the turbulence mixing factor was studied thermal hydraulics characteristics of 2 × 2 rod bundle by numerical analyzing [13]. In total of 2000 mm length of the 4 × 4 rod bundle with a mixing grid installed at 2/3 height, turbulence was processed in a Reynolds-averaged Navier-Stokes method using a nonlinear and viscous model to simulate CFD [14]. A study on the effect of spacer grid distance considered the influence of spacer grid blockage ratio and spacer grid spacing on heat transfer in a vertical rod bundle [15]. Reynolds number of the working fluid in this paper is 24,800 to 30,100 to investigate the flow effect on heat transfer. As such, this study was carried out with reference to other studies on the two-phase model, turbulence model, and various rod bundle shapes and grids.
The present study addresses the local flow inside the rod bundle. As the inside of a bundle of nuclear fuel rods has a very narrow flow path, a two-phase flow study using three-dimensional heat transfer analysis and a boiling model is required for reactor safety design and analysis. For two-phase flow analysis, the bubble size is a very important factor for accurately predicting bubble behavior in the flow path. In line with this focus, this study analyzes the results of the numerical analysis of the volume fraction of vapor (VF) and interfacial area concentration (IAC) compared with experimental values. VF is the volume of vapor in volume rate and IAC is defined by VF and interaction length scale. It focuses on the analysis of coolant boiling in a narrow flow path inside a nuclear fuel reactor 4 × 4 rod bundle. This phenomenon is solved by the Eulerian multiphase model as a segregated flow. The important elements to accurately predict the behavior of vapor are the interfacial momentum transfer force and wall boiling model of heat and mass transfer force. The usability of the detailed model used in this study has been evaluated by researchers as follows. Tomiyama (1998) studied simple but reliable correlations for a drag coefficient by comparing experimental data. The correlations were obtained for the drag coefficient of single bubbles under a wide range of fluid properties, bubble diameters, and acceleration of gravity. These were developed according to the balance of forces acting on the bubbles in stagnant liquid, and the available empirical correlations of the terminal rising velocities of single bubbles [16].
It has been shown that in many flows, the inertial force can simply be added to the generalization of the lift force introduced by Auton (1987) [17]. The lift force of Auton (1988) consists of the rotational force and inertial force or additional mass [18].
Owing to the neglect of the wall area force, the predicted pore profile peaked instead of approaching zero at the wall. Therefore, Achard & Cartellier (1985) concluded that a different kind of transverse force is needed to break the bubbles away from the wall [19]. Later, Antal (1991) developed a two-fluid model of multidimensional laminar bubble two-phase flow to analyze vertical pipe flow. In this study, the wall lift force for the dispersed phase (i.e., bubbles) and the wall lubrication force about the wall force were calculated. The results showed good agreement with the data when the considered model was used [20]. Lahey et al. (1993) used C TD = 1 to predict the lateral phase distribution of a three-dimensional water/air bubble upflow or downflow in a pipe using the k-ε turbulence model and an algebraic stress law. It can be confirmed [21]. These results are very interesting because a properly formulated multidimensional two-fluid model has the inherent ability to predict the lateral phase distribution phenomena in simple and complex geometric conduits for bubbles of size 1 mm < D b < 6 mm.
The model of Kurul and Podowski (1990) was used to determine the influence of bubble wall area fraction properties [22]. This model can estimate the fraction of the wall area affected by the sweeping of the liquid inflow under the bubble departure.
In the two-fluid model, active nucleation site density is an important factor for estimating the IAC. Thus, the model by Hibiki-Ishii (2003) proposed that the active nucleation site density is a function of the critical cavity size and contact angle [23]. This model clearly demonstrated the dependence of active nucleation site density on the superheated wall by various studies. Kocamustafaogullari (1983) reported that the bubble departure diameter was associated with water over a wide range of pressures [24]. Using the correlation proposed by Kocamustafaogullari, the mean deviation of the experimental data is approximately ±33%. It accurately simulated the available experimental water data, considering various free surfaces. Cole (1960) demonstrated through experimental measurements that the main forces acting on bubbles leaving the surface at high heat flows are the buoyancy and drag forces [25]. The dimensionless relationship was developed by the bubble velocity, bubble diameter, and contact angle at breakoff.
There are many other studies. However, the boiling phenomenon of the subchannel is still a field that requires various studies depending on the total length, hydraulic diameter, flow rate, and heat transfer rate. This study analyzes the coolant flow in 4 × 4 rod bundles with applied symmetrical conditions. The diameter, pitch, and hydraulic diameter of the rods are 16, 21, and 14.01 mm respectively, and the total length is 590 mm. The cross-sectional area of the coolant is 0.004008 m 2 . The grid size has three prism layers total of 0.8 mm near the wall, and 1 mm in other areas. Since the fluid belongs to high-Re, High-y+ wall treatment (applicable to 30 ≤ y+ ≤ 100) was applied. This simulation is calculated with steady state and constant density. Initial condition covered a mass flux ranging from 250 to 310 kg/m 2 s, heat flux ranging from 148 to 200 kW/m 2 , and inlet subcooled temperature ranging from 20.2 • C to 25.2 • C. Numerical analysis is performed by STAR-CCM+ ver. 12.06.011 considering three-dimensional turbulence.

Subcooled Boiling Experiment in 4 × 4 Rod Bundle
In this study, an experiment was conducted in the subcooled boiling water flow with a 4 × 4 rod bundle heater. Instruments such as Coriolis flow meters, static pressure gauges, differential pressure gauges, thermal couples, and power meters were equipped in the facility, and control valves for precisely controlling the flow rate and preheater to maintain the target inlet temperature were installed. The measurement instrument used in the experiment and the uncertainty are listed in Table 1. The pressure could be adjusted with a pressurizer connected to the water storage tank consisting of a test loop. A separator for separating steam and water was installed at the outlet of the test section, and a condenser and heat exchanger for heat removal from steam and water, respectively, were installed at the rear end of the water separator. An additional heat exchanger was installed at the outlet of the test section to prevent the flashing of high-temperature water. To minimize heat loss from the flow test channel during the experiment, insulating material was applied to the outside wall of the flow channel. The test section was designed as a square channel of which the length of the sides was 85 mm. Each heater rod was 16 mm in outer diameter and 590 mm in heated length (Z h ). The pitch-to-diameter ratio was 1.3 in the bundle arrangement; therefore, the hydraulic diameter (D H ) of the whole subchannel was 14 mm. A two-sensor optical fiber probe (2S-OFP) consisting of one front and one rear sensor was applied to measure local two-phase flow parameters in the rod bundle. The local two-phase flow parameters measured in the experiment were the void fraction (α t ), interfacial area concentration (a i t , IAC), interfacial area per unit volume, is a parameter of prime importance in terms of mass, momentum, and energy transfers between phases. In this study, the interfacial velocity is the velocity of the actual velocity of the bubbles along the axial direction of a flowing channel. The detailed equations for each parameter are as follows.
where T is the total measurement time, N b is the number of bubbles in contact with local optical fiber sensor, t F2 is the time when the bubble completely passes through the front sensor, t F1 is the time when the bubble starts to pass through the front sensor, t R1 is the time when the bubble starts to pass through the rear sensor, I φ j,max is the angle between the mean flow direction vector and the velocity vector of the interface, proposed by Hibiki et al. (1998) [26], and superscript t is the time average. The bubble diameter is six times the void fraction over the interfacial area concentration. Shen and Nakamura (2013) [27] showed a measurement error of ±7.8% when measuring the void fraction of small bubbles with an optical fiber probe technique. The interfacial velocity for small bubbles estimated with error less than ±2.3% when comparing the photographic and double sensor probe method used in our laboratory.
The measurement was performed on a total of 26 points in the area corresponding to 1/8 of the center subchannel at Z h /D H = 41.8 downstream from the inlet of the test section. The flow condition covered a mass flux ranging from 250 to 310 kg/m 2 s, heat flux ranging from 148 to 200 kW/m 2 , and an inlet subcooled temperature ranging from 20.2 • C to 25.2 • C at a 200 kPa inlet pressure.

Numerical Methods
For the analysis performed in this study, a two-fluid analysis method based on the Eulerian equation and a standard k-ε turbulence model were used [28,29]. In particular, this study considers the drag force of Tomiyama et al. (1998) [16], the lift force of Auton et al. (1988) [18], the wall lubrication force of Antal et al. (1991) [20], the turbulent dispersion force of Lahey et al. (1993) [21], and the virtual mass force of Auton et al. (1988). [18] The wall heat distribution model of Kurul and Podowski (1990) [22] was used to calculate the boiling flow, considering the interfacial length scale from Bak (2015) [30] of the bubble diameter model. Hibiki-Ishii (2003) [23] nucleation density, Kocamustafaogullari (1983) [24] bubble departure diameter, and the Cole (1960) [25] bubble departure frequency model were used to calculate the wall bubble generation.

Turbulence Model
Prior to the numerical analysis, the standard k-ε and realized k-ε turbulence models were compared to confirm the turbulence model sensitivity [31]. Figure 1 shows that there is no significant difference between the results of the two models. However, when calculations were made using the standard k-ε model under the same conditions, the results were closer to the experimental values. Therefore, the standard k-ε model was adopted for all cases.

Two-Phase Interaction Models
The drag force and coefficient of Tomiyama et al. (1998) [16] are given by In Equation (5), is the Reynolds number and is the Eotvos number, expressed as This method is suitable for liquid-gas bubble systems with a low Reynolds number and low Morton number flow; for example, small bubbles in an air-water system. It is only available when the dispersed phase is a gas. The lift force and coefficient of Auton et al. (1988) [18] are given by The wall lubrication force that prevents bubbles from touching the wall, and the coefficient of Antal et al. (1991) [20] are given by In Equation (11), and are constant. The turbulence dispersion force and coefficient of Lahey et al. (1993) [21] are given by

Two-Phase Interaction Models
The drag force and coefficient of Tomiyama et al. (1998) [16] are given by In Equation (5), Re is the Reynolds number and Eo is the Eotvos number, expressed as This method is suitable for liquid-gas bubble systems with a low Reynolds number and low Morton number flow; for example, small bubbles in an air-water system. It is only available when the dispersed phase is a gas. The lift force and coefficient of Auton et al. (1988) [18] are given by The wall lubrication force that prevents bubbles from touching the wall, and the coefficient of Antal et al. (1991) [20] are given by In Equation (11), C 1 and C 2 are constant. The turbulence dispersion force and coefficient of Lahey et al. (1993) [21] are given by Appl. Sci. 2020, 10, 4559 6 of 17 The virtual mass force and coefficient of Auton et al. (1988) [18] are given by The nucleation site number density determines the number of locations on the heated surface where bubbles form, per unit area. The nucleation site number density of Hibiki-Ishii (2003) [23] was used as follows where d 1 is 2.64 × 10 −5 m/deg, and θ is 41.36 • for the active nucleation site number density as a function of critical cavity size and contact angle. It is designed for use with the Kocamustafaogullari (1983) [24] model for bubble departure diameter as follows The bubble departure frequency model of Cole (1960) [25] was used as follows

Bubble Diameter Model
In this three-dimensional two-phase flow analysis, the bubble size is a very important factor for accurate prediction of bubble behavior in the flow path. Yao and Morel (2004) and Yeoh and Tu (2005) developed a kinematic model based on the IAC and bubble number density transport equations for accurate bubble size prediction. The model has been developed and requires subsidiary equations [32,33].
The bubble size is known to be related to the hydraulic diameter of the flow path, Laplace length scale, volume fraction, and Reynolds number of the bubble through the previously developed models of Zeitoun et al. (1994) [34], Zeitoun and Shoukri (1996) [35], and Hibiki et al. (2006) [36]. Thus, bubble size model is formed by where α is the volume fraction. The bubble Reynolds number (N Reb ) to consider the turbulent dissipation term (ε) and the Laplace length scale (L O ) are given by In Equation (20), the subscript f means liquid phase. In Equation (21), σ is the surface tension, g is the gravitational acceleration, and ∆ρ is the density difference between phases.
Based on Equation (19) (20) is derived by adding the density ratio and excluding the hydraulic diameter included in the bubble Reynolds number [37]. Bak derived a bubble diameter model from Yun by applying DEDALE, Hibiki, and DEBORA experimental data for both low-and high-pressure conditions. The bubble size model of Bak (2015) [30] used in this study is given by where the density ratio (N ρ ) to reflect pressure is given by and the subscripts f and g are the liquid phase and gas phase, respectively.

Operating Conditions
The numerical analysis of the rod bundle has six initial conditions. It also proves that the grid test was conducted before showing the numerical analysis results.

Initial and Boundary Conditions
The total length of the rod bundles was 590 mm. The diameter, pitch, and hydraulic diameter are 16, 21, and 14.02 mm, respectively. The cross-sectional area through which the coolant flows is approximately 0.00401 m 2 . The experimentally measured data on Lines 1-5 comparing the numerical analysis results at the cross-section of 580 mm in height are shown in Figure 2. . .

/ (22)
A simple form of Equation (20) is derived by adding the density ratio and excluding the hydraulic diameter included in the bubble Reynolds number [37]. Bak derived a bubble diameter model from Yun by applying DEDALE, Hibiki, and DEBORA experimental data for both low-and high-pressure conditions. The bubble size model of Bak (2015) [30] used in this study is given by where the density ratio ( ) to reflect pressure is given by = (24) and the subscripts f and g are the liquid phase and gas phase, respectively.

Operating Conditions
The numerical analysis of the rod bundle has six initial conditions. It also proves that the grid test was conducted before showing the numerical analysis results.

Initial and Boundary Conditions
The total length of the rod bundles was 590 mm. The diameter, pitch, and hydraulic diameter are 16, 21, and 14.02 mm, respectively. The cross-sectional area through which the coolant flows is approximately 0.00401 m². The experimentally measured data on Lines 1-5 comparing the numerical analysis results at the cross-section of 580 mm in height are shown in Figure 2. The initial conditions for the six experimental cases used in the calculations are listed in Table 2. At the outlet of all cases, the fluid has a pressure of 194.45 kPa. This shows a decrease of about 5 kPa from the inlet. Calculation was performed with constant density.  The initial conditions for the six experimental cases used in the calculations are listed in Table 2. At the outlet of all cases, the fluid has a pressure of 194.45 kPa. This shows a decrease of about 5 kPa from the inlet. Calculation was performed with constant density. In Table 2, the different conditions are highlighted in gray. Cases 2 and 3 have different inlet temperatures compared with Case 1. Cases 4 and 5 have different heat flux compared with Case 1. Case 6 has a different inlet velocity from that of Case 3.

Grid Test
Given that the fluid does not have different results at a centerline of 4 × 4 and 2 × 2, the numerical analysis of the rod bundle was modeled by constructing a shape with one 2 × 2 side and by applying a symmetrical condition twice. The geometry part is the gray area in Figure 3b.
Appl. Sci. 2020, 10, x 8 of 17 In Table 2, the different conditions are highlighted in gray. Cases 2 and 3 have different inlet temperatures compared with Case 1. Cases 4 and 5 have different heat flux compared with Case 1. Case 6 has a different inlet velocity from that of Case 3.

Grid Test
Given that the fluid does not have different results at a centerline of 4 × 4 and 2 × 2, the numerical analysis of the rod bundle was modeled by constructing a shape with one 2 × 2 side and by applying a symmetrical condition twice. The geometry part is the gray area in Figure 3b.  The analysis results are shown in Figure 4, with approximately one million, two million, and four million grids. A comparison of the analysis values was made between Line 1 and Line 4.  The grid test showed no significant difference between the approximately one million, two million, and four million grids. Therefore, numerical analysis was performed with one million grids, The analysis results are shown in Figure 4, with approximately one million, two million, and four million grids. A comparison of the analysis values was made between Line 1 and Line 4.
Appl. Sci. 2020, 10, x 8 of 17 In Table 2, the different conditions are highlighted in gray. Cases 2 and 3 have different inlet temperatures compared with Case 1. Cases 4 and 5 have different heat flux compared with Case 1. Case 6 has a different inlet velocity from that of Case 3.

Grid Test
Given that the fluid does not have different results at a centerline of 4 × 4 and 2 × 2, the numerical analysis of the rod bundle was modeled by constructing a shape with one 2 × 2 side and by applying a symmetrical condition twice. The geometry part is the gray area in Figure 3b.  The analysis results are shown in Figure 4, with approximately one million, two million, and four million grids. A comparison of the analysis values was made between Line 1 and Line 4.  The grid test showed no significant difference between the approximately one million, two million, and four million grids. Therefore, numerical analysis was performed with one million grids,  The grid test showed no significant difference between the approximately one million, two million, and four million grids. Therefore, numerical analysis was performed with one million grids, which had the shortest calculation time. The cross-sectional shape of the one million grid is shown in Figure 5.
Appl. Sci. 2020, 10, x 9 of 17 which had the shortest calculation time. The cross-sectional shape of the one million grid is shown in Figure 5. The one million grid consists of hexahedron grids. By applying three prism layers to the wall surface, the spacing is approximately 0.2667 mm near the wall, for a total of 0.8 mm. In other areas, the grid spacing was 1 mm. This grid method calculates three important measures of mesh quality at the start of a run: aspect ratio-maximum 0.99: good; skewness angle-maximum 4.79: good; volume change-maximum 0.98: good. Mean aspect ratio is 0.81, mean skewness is 4.75, mean volume change is 0.8. In terms of overall mesh quality minimum quality is 0.53, maximum 0.99, and mean 0.97.

Results
Before the results of the six cases are shown, Figure 6 illustrates a stable flow falling by approximately 5.53 kPa of pressure from about 200 kPa from the inlet to the outlet. The pressure distribution is shown in Figure 6. For the parameters of the results, the meanings of the volume fraction and IAC are as follows.
In Equation (26), the interaction length scale is defined as the bubble diameter. The one million grid consists of hexahedron grids. By applying three prism layers to the wall surface, the spacing is approximately 0.2667 mm near the wall, for a total of 0.8 mm. In other areas, the grid spacing was 1 mm. This grid method calculates three important measures of mesh quality at the start of a run: aspect ratio-maximum 0.99: good; skewness angle-maximum 4.79: good; volume change-maximum 0.98: good. Mean aspect ratio is 0.81, mean skewness is 4.75, mean volume change is 0.8. In terms of overall mesh quality minimum quality is 0.53, maximum 0.99, and mean 0.97.

Results
Before the results of the six cases are shown, Figure 6 illustrates a stable flow falling by approximately 5.53 kPa of pressure from about 200 kPa from the inlet to the outlet. The pressure distribution is shown in Figure 6.
Appl. Sci. 2020, 10, x 9 of 17 which had the shortest calculation time. The cross-sectional shape of the one million grid is shown in Figure 5. The one million grid consists of hexahedron grids. By applying three prism layers to the wall surface, the spacing is approximately 0.2667 mm near the wall, for a total of 0.8 mm. In other areas, the grid spacing was 1 mm. This grid method calculates three important measures of mesh quality at the start of a run: aspect ratio-maximum 0.99: good; skewness angle-maximum 4.79: good; volume change-maximum 0.98: good. Mean aspect ratio is 0.81, mean skewness is 4.75, mean volume change is 0.8. In terms of overall mesh quality minimum quality is 0.53, maximum 0.99, and mean 0.97.

Results
Before the results of the six cases are shown, Figure 6 illustrates a stable flow falling by approximately 5.53 kPa of pressure from about 200 kPa from the inlet to the outlet. The pressure distribution is shown in Figure 6. For the parameters of the results, the meanings of the volume fraction and IAC are as follows.
In Equation (26), the interaction length scale is defined as the bubble diameter. For the parameters of the results, the meanings of the volume fraction and IAC are as follows.
Volume f raction = Volume o f vapor Volume (25) Volume f raction Interaction length scale (26) In Equation (26), the interaction length scale is defined as the bubble diameter. Figures 7-12 show the comparison of the results of Lines 1-5 in Cases 1-6. In each figure, (a) shows the volume fraction and (b) shows the IAC.
In Case 1, it is considered that the VF and IAC were accurately predicted. Cases 2 and 3 (subcooled temperature of 20.2 • C or 25.2 • C) have different inlet temperature from Case 1 (subcooled temperature of 22.7 • C). The axis shows numerical values similar to experimental values, but has smaller values to the rod surface. There was no difficulty in interpreting the same tendency, such as the slight decrease in numerical results near the rod surface, as the experimental results. The result of Case 4 shows accurate prediction. Case 4 also indicates low VF and IAC, the same as the experiments, because of the low heat flux of 148 kw/m 2 . In Case 5, which had a high heat flux of 200 kw/m 2 , the IAC was low near the rod surface, but it can be said that it measures some reasonable value in the VF. Lastly, Case 6, which had a low flow rate, was under-interpreted near the rod surface except on Line 1, but it shows that it interprets the appropriate values near the axis and across Line 1. In Case 1, it is considered that the VF and IAC were accurately predicted. Cases 2 and 3 (subcooled temperature of 20.2 °C or 25.2 °C) have different inlet temperature from Case 1 (subcooled temperature of 22.7 °C). The axis shows numerical values similar to experimental values, but has smaller values to the rod surface. There was no difficulty in interpreting the same tendency, such as the slight decrease in numerical results near the rod surface, as the experimental results. The result of Case 4 shows accurate prediction. Case 4 also indicates low VF and IAC, the same as the experiments, because of the low heat flux of 148 kw/m². In Case 5, which had a high heat flux of 200 kw/m², the IAC was low near the rod surface, but it can be said that it measures some reasonable value in the VF. Lastly, Case 6, which had a low flow rate, was under-interpreted near the rod surface except on Line 1, but it shows that it interprets the appropriate values near the axis and across Line 1.  In Case 1, it is considered that the VF and IAC were accurately predicted. Cases 2 and 3 (subcooled temperature of 20.2 °C or 25.2 °C) have different inlet temperature from Case 1 (subcooled temperature of 22.7 °C). The axis shows numerical values similar to experimental values, but has smaller values to the rod surface. There was no difficulty in interpreting the same tendency, such as the slight decrease in numerical results near the rod surface, as the experimental results. The result of Case 4 shows accurate prediction. Case 4 also indicates low VF and IAC, the same as the experiments, because of the low heat flux of 148 kw/m². In Case 5, which had a high heat flux of 200 kw/m², the IAC was low near the rod surface, but it can be said that it measures some reasonable value in the VF. Lastly, Case 6, which had a low flow rate, was under-interpreted near the rod surface except on Line 1, but it shows that it interprets the appropriate values near the axis and across Line 1.      Figure 13 shows the numerical analysis results for the volume fraction distribution contour of all cases for the cross-section with a height of 580 mm. In Figure 13c,d, the volume fraction distributions are characterized by a low value between the rod surface and the wall. However, in Figure 13a,b,e,f, the volume fraction distributions are characterized by a high value between the rod surface and the wall.  Figure 13 shows the numerical analysis results for the volume fraction distribution contour of all cases for the cross-section with a height of 580 mm. In Figure 13c,d, the volume fraction distributions are characterized by a low value between the rod surface and the wall. However, in Figure 13a,b,e,f, the volume fraction distributions are characterized by a high value between the rod surface and the wall.  Figure 13 shows the numerical analysis results for the volume fraction distribution contour of all cases for the cross-section with a height of 580 mm. In Figure 13c,d, the volume fraction distributions are characterized by a low value between the rod surface and the wall. However, in Figure 13a,b,e,f, the volume fraction distributions are characterized by a high value between the rod surface and the wall.  Figure 14a, the liquid velocity decreased more as it approached the outlet in the wall edge. Therefore, the liquid temperature increased because of heat transfer. However, because Case 4 was a case with very low heat flux, the temperature at the wall edge excluding the rod surface did not rise sufficiently to the boiling point, as shown in Figure 14b. Therefore, as shown in Figure 14c, air bubbles were not generated at the corners. Because no bubbles were generated, the velocity of the bubbles remained low, as shown in Figure 14d.   Figure 14a, the liquid velocity decreased more as it approached the outlet in the wall edge. Therefore, the liquid temperature increased because of heat transfer. However, because Case 4 was a case with very low heat flux, the temperature at the wall edge excluding the rod surface did not rise sufficiently to the boiling point, as shown in Figure 14b. Therefore, as shown in Figure 14c, air bubbles were not generated at the corners. Because no bubbles were generated, the velocity of the bubbles remained low, as shown in Figure 14d.  Figure 14a, the liquid velocity decreased more as it approached the outlet in the wall edge. Therefore, the liquid temperature increased because of heat transfer. However, because Case 4 was a case with very low heat flux, the temperature at the wall edge excluding the rod surface did not rise sufficiently to the boiling point, as shown in Figure 14b. Therefore, as shown in Figure 14c, air bubbles were not generated at the corners. Because no bubbles were generated, the velocity of the bubbles remained low, as shown in Figure 14d.   Figure 15a, however, shows that the liquid flowed with a low velocity at the wall edge from the inlet to the middle height of the rod bundle. Therefore, as the liquid velocity decreased, the temperature increased enough to reach the boiling point from the middle height, as shown in Figure  15b. Consequently, more bubbles were generated, as shown in Figure 15c. Figure 15d shows that the velocity of the liquid increased more as the generated vapor velocity increased.

Model Validation
Various initial conditions were used for the numerical analysis of the six cases, with inlet fluid velocity of approximately 0.260417 or 0.322917 m/s, heat flux of 148-200 kw/m², and inlet subcooled temperature of 20.2-25.2 °C. The graphs on the cross-section near the outlet from the rod can confirm the validity of the model in  In particular, this model predicts the reasonable value on center of rod bundles in a wide flow path and the high inlet temperature or the large heat flux; the sub-cooled temperature is 20.2 °C or the heat flux is 200 kw/m². All the volume fractions and IACs in Figures 7-12 tend to be similar to the experimental values, which indicates that this study will be useful for understanding the vapor behavior of the entire steady-state flow in the simulations for predicting the boiling phenomenon. The results shown in Figures 14 and 15 are expected to aid in predicting the correlation between liquid and volume fraction, which is the most important factor considered in bubbly flow CFD research. These CFD results generally reproduced the trend of experimental values. However, this simulation has the same kind of error in all cases. Thus, we must describe the error. Numerical analysis results in higher errors as the mass flux is higher and the heat flux and temperature are lower. Moreover, it can be confirmed that the VF is smaller than the experimental value. The analysis of this cause from the model aspect shows that is due to the use of a constant value for each coefficient of the interface momentum transfer force models in this paper. In each model, the coefficient was inputted as a certain and reasonable constant with reference to the papers, but it is not a value derived from the exact same experiment and the same working fluid, so  Figure 15a, however, shows that the liquid flowed with a low velocity at the wall edge from the inlet to the middle height of the rod bundle. Therefore, as the liquid velocity decreased, the temperature increased enough to reach the boiling point from the middle height, as shown in Figure 15b. Consequently, more bubbles were generated, as shown in Figure 15c. Figure 15d shows that the velocity of the liquid increased more as the generated vapor velocity increased.

Model Validation
Various initial conditions were used for the numerical analysis of the six cases, with inlet fluid velocity of approximately 0.260417 or 0.322917 m/s, heat flux of 148-200 kw/m 2 , and inlet subcooled temperature of 20.2-25.2 • C. The graphs on the cross-section near the outlet from the rod can confirm the validity of the model in Figures 7-12. In particular, this model predicts the reasonable value on center of rod bundles in a wide flow path and the high inlet temperature or the large heat flux; the sub-cooled temperature is 20.2 • C or the heat flux is 200 kw/m 2 . All the volume fractions and IACs in Figures 7-12 tend to be similar to the experimental values, which indicates that this study will be useful for understanding the vapor behavior of the entire steady-state flow in the simulations for predicting the boiling phenomenon. The results shown in Figures 14 and 15 are expected to aid in predicting the correlation between liquid and volume fraction, which is the most important factor considered in bubbly flow CFD research. These CFD results generally reproduced the trend of experimental values. However, this simulation has the same kind of error in all cases. Thus, we must describe the error. Numerical analysis results in higher errors as the mass flux is higher and the heat flux and temperature are lower. Moreover, it can be confirmed that the VF is smaller than the experimental value. The analysis of this cause from the model aspect shows that is due to the use of a constant value for each coefficient of the interface momentum transfer force models in this paper.
In each model, the coefficient was inputted as a certain and reasonable constant with reference to the papers, but it is not a value derived from the exact same experiment and the same working fluid, so an error occurs. In the interfacial momentum transfer force model that predicts the movement of bubbles, the most dominant force is the drag force, and then the lift force. When looking at the horizontal cross section of a vertical pipe, the lift force, which determines the location of the bubble, is one of the important factors in determining how well the bubble moves after generation and departure from the heating surface [38][39][40].

Conclusions
In this study, numerical analysis of two-phase flow boiling by heat transfer was conducted.
The experimental values and CFD results showed good agreement in all six cases considered. The analysis was a steady state of 3D and the calculated model was performed at constant density. The fluid was calculated using the Eulerian multiphase model including a Reynolds stress turbulence model. The interfacial momentum transfer force model and wall boiling model are important models in two-phase flow. The momentum transfer force was calculated by considering the drag force, lift force, wall lubrication force, turbulent dispersion force, and virtual mass force. The nucleation site number density, bubble departure diameter, and bubble departure frequency of the wall boiling model were considered in the heat and mass transfer force. This means that the model considered the behavioral characteristics of the fluid to be as realistic as possible. Thus, it is possible to sufficiently interpret the code. In particular, it shows stable analysis results in all cases with different conditions. The main contributions of this study are as follows.

1.
It enables various water-air flow analyses of low-and high-pressure conditions of the model used in numerical analysis.

2.
When the error rate is very small in a wide flow path, and the inlet temperature is high or the heat flux is large, the flow boiling phenomenon is properly predicted.

3.
Through an analysis of vertical plane flow, which cannot be measured by experiments, the cause of the high-volume fraction at the edge was investigated by using the liquid velocity and correlation.