Assessment of a Theoretical Model for Predicting Forced Convective Critical Heat Flux in Rod Bundles

A theoretical model for predicting forced convective critical heat flux (CHF) in rod bundles was proposed, which is based on the bubble crowding phenomenon. The theoretical model applied to the rod bundle is based on the Weisman-Pei's basic tube model. In order to make it suitable for Pressurized Water Reactor (PWR) in rod bundles, the flow and heat transfer characteristics of rod bundles are considered, including velocity distribution, flow patterns, grid effect on CHF and turbulent intensity. The theoretical model has been applied together with the subchannel code ATHAS for assessing against uniformly heated CHF data obtained through a 5x5 rod bundle, and good agreement has been observed.


INTRODUCTION
Effective and accurate prediction of critical heat flux (CHF) is essential to ensure the safe operation of forced convective equipment. It is of particular concern in the operation of water-cooled nuclear reactors, where cladding and core integrities must be maintained and a safe operating power envelope and margin must be established. Moreover, it is closely related to the economy of the reactors.
Theoretical models (Weisman, 1991;Celata et al., 1994;Bruder et al., 2015) have received the attention of some scholars because of their incorporation of physical mechanisms, accurate parameter trends, and wide application range, which provide flexibility in predicting CHF for new bundle concepts or configurations, such as those proposed for small modular reactors (SMRs). The bubble crowding model (Weisman, 1991) is a successful theoretical model for predicting forced convective critical heat flux in a tube. The model was first proposed by Weisman and Pei (1983) based on bubble crowding. The main idea is that the bubbles generated on the heating wall hinder the radial flow of liquid from the bulk flow area to the near-wall bubble layer area. Ying and Weisman (1986) modified the model to accommodate the non-uniform void profile in a tube and extended the CHF prediction for void fractions up to 0.8. In addition, they revised the calculation for bubble diameters and included the slip ratio in the bubble layer to improve CHF prediction at low mass flow rates. Weisman and Illeslamlou (1988), on the other hand, extended the model to high subcooling conditions based on the energy balance at the outer edge of the bubble layer in round tubes. Chang and Lee (1989) revised the calculation of the lateral mass velocity from the core to the bubble layer at low qualities in uniformly heated tubes in their CHF model. Kwon and Chang (1999) introduced a drag force due to the roughness of wall-attached bubbles in the momentum balance, which determines the limiting transverse interchange of mass flux crossing the interface of the bubbly wall layer and the core. Furthermore, the bubbly layer was assumed to be a single layer of wall-attached bubbles that acts as an equivalent of surface roughness. The critical void fraction at the bubble layer was represented with an exponential function related to the quality, which was determined from CHF data for round tubes. Kodama and Kataoka (2002) expressed the critical void fraction at the bubble layer in terms of the channel-average void fraction, which was also determined from CHF data for round tubes. Kinoshita et al. (2001) considered spherical bubbles in the bubble layer at high velocity and subcooling conditions. These bubbles were assumed to contact each other in a cubic configuration. The interference of bubbles was taken at the ratio of bubble diameter to bubble distance, equaling 0.5. The critical void fraction at the bubble layer was determined to be π/12.
Previous studies on bubble crowding models are based on a tube, which is not suitable for rod bundle structures. The objective of this study is to develop a generalized theoretical CHF model in rod bundles based on the bubble crowding concept (Weisman and Pei, 1983). The radial velocity distribution, grid enhancement effect, and flow pattern change in rod bundles are considered.

DESCRIPTION OF THE THEORETICAL CHF MODEL
It is generally believed that the mechanism of CHF under subcooling and low quality is different from that under high quality. To distinguish between the two mechanisms, subcooled and low-quality CHF is called "departure from nucleate boiling" (DNB), and the other is termed dry-out.

Physical Mechanism and Basic Equations
The derived prediction procedure of Weisman and Pei's tube CHF theoretical model (Weisman and Pei, 1983) is based on the bubble crowding concept. At the heating wall, bubbles will be generated. CHF will occur when there are too many bubbles to wet the heating wall effectively with cold liquid from the bulk flow area, as shown in Figure 1.
They divide the tube area into the bulk flow area and the near-wall bubble layer area. For these two areas, mass, and energy conservation equations are established to derive the CHF expression as follows: where q CHF is the predicted value of CHF calculated by the model, h fg is the latent heat of evaporation, x 2 is the average quality at the bubble layer and x 1 is the average quality at the core layer, h f is the saturated liquid enthalpy, h l is the enthalpy of liquid, and h ld is the enthalpy at the point of bubble detachment. G 3 is the lateral mass velocity from the core to the bubble layer due to turbulence, which is determined by: in which G 3 is expressed in terms of the total axial mass velocity, G, the turbulence intensity at the interface between the bubbly layer and the core, i b , and a miscellaneous function, ψ, which represents the share of liquid reaching the wall in the bulk flow. Weisman and Pei (1983) expressed ψ as follows: where v' is the radial fluctuating velocity, σ v ′ is the standard deviation of v, and v 1 l is the radial velocity created by vapor generation.
According to Lahey and Moody (1977), q b , the portion of the total heat flux effective in generating vapor, has the following relationship with total heat flux, q: The following is the derivation of turbulence intensity i b in this study. Trupp and Azad's (1973) measurements of turbulent radial velocity fluctuations for a P/D of 1.35 (which is close to that of the PWR fuel assembly) are applied in establishing the radial turbulence fluctuation in the subchannel. Lee and Durst (1980) pointed out that the ratio (v ′ ) 2 /U τ /(l e /r 0 ) does not depend on the Reynolds number and can be considered to be only a function of (r/r0). Figure 2 compares the measurements of Trupp and Azad (1973) and those of a tube (Laufer, 1953). In addition, parabolic fits of all of and the near-wall (up to radius ratios, r/r0, of 0.5) measurements of Trupp and Azad (1973) are also shown.
In view of the CHF occurrence at the wall, the near-wall velocity fluctuation (pink line) is of the most interest and is expressed as: where U τ is the frictional velocity, l e is the Prandtl mixing length, r 0 is the outer radius of the tube, and y is the radial distance from the wall. With the assumption that the ratio of two-phase to single-phase turbulence intensity is independent of radial position, we have Equation (7): And the general relationship of l e is: The frictional velocity U τ is defined as: where f in the turbulent region is: From Equations (6) to (10), we have: Weisman and Pei (1983) assumed the distance from the wall at which the bubbly layer-core interface occurs, y c , is: where k is an undetermined coefficient and D P is the bubble detachment diameter. By replacing 'y' in Equation (11) with 'y c ' in Equation (12), we have: As indicated by Weisman and Pei, the two-phase factor F 1 represents the effect of bubble motion on turbulence intensity, and the two-phase factor F 2 accounts for the effect of bubble motion on the fluctuation velocity. Weisman and Pei assumed the ratio F 1 /F 2 to be in the form: The physical meaning of Equation (14) is the effect of two-phase flow on turbulence intensity.

Bubble Detachment Point and Detachment Diameter
Reliable predictions of the onset of bubble departure and detached bubble diameter are essential for modeling in Equations (1) and (13). In a review in which an extensive amount of data was compared, Lee et al. (1992) found the Levy model (Levy, 1966) to achieve the best fit with the whole set of data of the many analytical models. Levy proposed the bubble departure enthalpy, which is defined as the liquid enthalpy at which vapor begins to break away from a heated surface as: where, and H db is the single-phase heat transfer coefficient given by the Dittus-Boelter equation: Meanwhile, the bubble detachment diameter, D P , which is determined from a balance of the fluid forces acting on the bubble and the surface force, was given by Levy (1966) as: Grid Enhancement on Turbulence Intensity i b In rod bundles, grids enhance the turbulence intensity in the flow stream and, accordingly, the heat transfer. This enhancement effect is captured in the turbulence intensity factor, i b , of the present model. Figure 3 schematically illustrates the enhancement effect of grids on CHF in a bundle with axially uniform heating. The local CHF is reduced monotonically along the bare bundle (i.e., without grids) due to the increase in flow quality and reaches a minimum at the end of the heated length. In the presence of grids, the local CHF is enhanced at the grid location but decays gradually with increasing distance downstream from the grid. The CHF decays until encountering the next grid, where the enhancement recovers, or it returns back to that of a bundle without grids if the distance between neighboring grids is long.
The enhancement effect of the grid on CHF is proportional to the enhancement in turbulence intensity at locations downstream of the grid. Nagayoshi and Nishida (1998) studied the radial turbulent velocity distribution of the straight-type grid and observed a decay in the radial turbulence intensity factor downstream of the grid. The normalized velocity fluctuation reached a maximum value at the grid and decayed with increasing distance from the grid. It approached the level for the bare bundle after a distance of approximately 10 times the hydraulicequivalent diameter from the grid. Nagayoshi and Nishida (1998) expressed the change in radial turbulence intensity factor as: where x is the axial distance downstream from the grid, D is the hydraulics-equivalent, and ε is the blockage-area ratio of the grid.
A grid correction factor, F grid , was introduced into the radial turbulence intensity factor i b (Equation 13) in this study to include the enhancement effect, which is expressed as: The coefficients a 2 and b 2 depend on the type of grid. The experimental data of Yao et al. (1982) indicated that the heattransfer enhancement due to a grid falls down to the reference value (i.e., a bundle without grids) at a distance of 25-30 times the hydraulics-equivalent diameter. Based on these data, the coefficient b 2 has been established as−0.13. A value of 6.5 has been derived for the coefficient a 2 from the data of Nagayoshi and Nishida (1998).

Boiling Effect on Turbulence Intensity i b
The turbulence intensity factor, i b , expressed in Equation (13), is based primarily on single-phase parameters. Weisman and Pei (1983) introduced the ratio (F 1 /F 2 ), the physical meaning of which is the effect of boiling or two phases on turbulence intensity, as shown in Equation (14), rewritten here as: In Weisman and Pei's study, the coefficient "a" to account for the enhancement in turbulence intensity for two-phase flow, is expressed as: where Gcr is the reference mass flow rate at 970,000 kg/m2/h. Lim (1988) indicated a slight pressure dependency for the coefficient "a" based on their experimental data. The pressure effect is captured by associating it with the void fraction in this study. Therefore, we consider a function with a void fraction dependent variable to describe the effect of pressure. Through the analysis of the bundle CHF data in section assessment and analysis of the theoretical CHF model, we find that the two-phase factor has a closer relationship with the void fraction, and the influence of the void fraction on the two-phase factor is obviously different between the high void fraction and the low void fraction, which is due to the change in flow pattern. In the square rod bundle flow pattern experiments by Venkateswararao et al. (1982), they divided the flow pattern into bubbly flow, dispersed bubble flow, slug flow, churn flow, and annular flow when the ratio of the rod bundle spacing to the diameter (S/D) is about 1.38. They also defined the boundary between dispersed bubbly flow and churn flow as where void fraction "α" equals 0.52.
By comparing the experimental data and the flow pattern map given by Venkateswararao et al. (1982), we find that all of the bundle CHF data in section assessment and analysis of the theoretical CHF model are in the region of dispersed bubble flow and churn flow. Therefore, we can conclude that the boundary line due to the different flow patterns leads to different effects of the void fraction on the two-phase effect. The expression of Equation (22) can be simplified to: According to the above boundary conditions, the boundary void fraction 'α' equals 0.52 between dispersed bubble flow and churn flow. We used the CHF data in section assessment and analysis of the theoretical CHF model below to optimize the coefficients and found that b 1 equals 1.7 and a 3 equals 197.9 when α is greater than or equal to 0.52, and b 1 equals −0.06 and a 3 equals 69.65 when α is smaller than 0.52. For coefficient a 1 , we still use the form of Weisman and Pei (1983).

Other Relations
According to the assumption that there is homogeneous flow in both the bubble layer and the bulk flow region, this study follows Weisman and Pei's (1983) original idea that when CHF occurs, the void fraction in the bubble layer α 2 equals 0.82, which was based on the maximum packing density of ellipsoids with an axis ratio of three to one. The two-phase flow parameters are calculated as follows: x 1 = α 1 × ρ g ρ 1 (28) Subscript 1 denotes the parameters in the bulk flow, and subscript 2 denotes the parameters in the bubble layer. "r" is the channel radius and "s" is the bubbly layer thickness, which equals 5.5 times D P , as recommended by Weisman and Pei (1983). The detailed calculation process of CHF model is shown in the Appendix A.

CHF Data Bank
The experiments were performed at the Nuclear Power Institute of China (NPIC) with water flow over a 5 × 5 rod bundle simulating a PWR fuel assembly with flat mixing vanes. Three test bundles were constructed for the experiments. Two of these bundles simulated a PWR fuel assembly with nine hot rods at the central locations and 16 cold rods at the peripheral locations inside a square frame 66.1 mm in width. Each rod had an outer diameter of 9.5 mm. The spacing between heated rods and between the cold rods and the frame was 3.1 mm (i.e., the rod pitch was 12.6 mm). One of these bundles was heated over a length of 2438 mm (or 8 ft) and the other over a length of 3657 mm (or 12 ft). Figure 4 illustrates the rod configuration of the 5 × 5 rod bundle experiment. Eight grids were installed along the axial length of the short bundle [four of these grids were simple support grids without mixing vanes (SS), while four had mixing vanes installed in the grid (MVG)]. Thirteen grids were installed along the axial length of the long bundle at the locations shown in the figure (six SS grids and seven MVG grids). Figure 5 illustrates the grid configurations of the long and short bundles.
The span between grids with mixing vanes was 0.56 m. A set of thermocouples was installed at the location 56 mm upstream of the end of the heated length (10 mm upstream of the last grid with mixing vanes). The coolant traveled vertically upward from the bottom to the top.
The third test bundle in the experiments was equipped with a guide tube of 12.45 mm in outer diameter that replaced the hot rod at the central position (see Figure 4). The spacing between the hot rods and the guide tube was 1.625 mm. Thirteen grids were installed over the heated length of 3657 mm. The other geometric configurations of this bundle were the same as those of the long bundle described above. Six sets of data have been selected from the database to assess the prediction accuracy of the bubble-crowding-based mechanistic CHF model for subchannels in the test bundles. These data covered bundles of both short and long heated lengths with and without the guide tube (GT). Table 1 lists the overall flow conditions covered by these data, which are 417 data under the PWR operations.

Prediction Results of Bundle CHF Data
CHF is predicted using the bubble-crowding-based mechanistic CHF model with the local flow conditions at each subchannel along the axial distance of the bundle evaluated using the subchannel code ATHAS, which was developed by the Nuclear Safety and Operation Laboratory at Xi'an Jiao Tong University (Liu et al., 2014). The ATHAS code is based on the drift-flow basic model and can be used for core thermal-hydraulic analysis of PWRs and BWRs. Local flow conditions in the subchannel are directly applied in the model to evaluate the CHF. This is referred to as the Direct Substitution Method (DSM) for determining CHF, which is more commonly used in PWRs. The DSM method uses the experimental inlet mass rate flux, inlet temperature, and CHF value as the input of the subchannel code and calculates the real local parameters as the inputs of the   CHF model. A turbulent mixing coefficient of 0.066 is applied for the type of mixing grid installed in the tested bundles in ATHAS code. The predicted CHF value is compared against the experimental heat flux to determine the DNBR in the subchannel. The location with the minimum DNBR (or MDNBR) is considered the initial CHF point. The prediction accuracy of the model is assessed from the average value of predicted MDNBRs, R av , the standard deviation on MDNBR, SD, and the root-mean-square error on MDNBR, RMS, which are defined as: where MDNBRi is the MDNBR for each data point and N is the total number of data points. Figure 6 compares the predicted and experimental CHF values for 417 data points. The average value of MDNBRs, R av , is 1.0025, with a standard deviation, SD, of 10.51%, and a root-mean-square error, RMS, of 10.50%. Almost all of the data are predicted within the ±30% error band. Applying the original model of Weisman and Pei (1983) leads to a R av of 0.88 with a standard deviation of 22% for the same set of data points. This demonstrates the improvement in prediction accuracy achieved by the current model.
Parametric trends of predicted MDNBR in the bundles are examined against the experimental conditions (i.e., pressure, mass flux, and quality) in Figures 7-9. There are no apparent trends of predicted MDNBRs with pressure (Figure 7), mass flux (Figure 8), and quality (Figure 9).

CONCLUSIONS
A new mechanistic CHF model for subchannels in a bundle has been developed. It is based on the bubble-crowding concept for tubes but considers the effects of subchannel geometry on axial and radial velocity distributions and the effects of boiling characteristics and spacer grids on CHF. In the model, the turbulence intensity factor has been revised to capture the change in velocity profile between a tube and a subchannel. It also includes the enhancement effect due to a grid. Moreover, the influence of flow pattern change on CHF is also considered. This mechanistic CHF model has been assessed against experimental data obtained with water flow through 5 × 5 rod bundles. Local flow conditions in subchannels of the bundle were calculated using the subchannel code "ATHAS." Experimental CHF values were predicted with an average value of 1.0025 and a standard deviation of 10.51% for 417 data points at PWR conditions of interest. The predicted MDNBRs appear independent of pressure, mass flux, and quality. More detailed local parameter distributions and evidence for CHF visualization need to be acquired for future optimization of this model.

DATA AVAILABILITY STATEMENT
The datasets analyzed in this manuscript are not publicly available. Requests to access the datasets should be directed to The Nuclear Power Institute of China [Contact mail: liuwei0958@126.com].

AUTHOR CONTRIBUTIONS
YL was responsible for the specific work of this paper. QY carried out some of the simulation work. JS guided the work of this article. BZ was responsible for the embedding of the model in the subchannel code. WL provided experimental data. LL optimized the structure and tone of this article.

FUNDING
This work was financially supported by the National Key R&D Program of China (No. 2018YFB1900402).