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

This study proposes a theoretical model for predicting forced convective critical heat flux (CHF) in rod bundles based on the bubble crowding phenomenon. The theoretical model applied to the rod bundle is based on Weisman-Pei's basic tube model. In order to make it suitable for rod bundles in Pressurized Water Reactors (PWR), the flow and heat transfer characteristics of rod bundles are considered, including velocity distribution, flow patterns, the grid effect on CHF, and turbulence intensity. The theoretical model is applied together with the subchannel code ATHAS to assess it against uniformly heated CHF data obtained from a 5 × 5 rod bundle, and good agreement is 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.

FIGURE 1

Figure 1. Bubble crowding in vertical flow.

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:

$\begin{array}{ll}{q}_{CHF}/(h_{fg}{G}_3)=(x_2-x_1)(\frac{h_f-h_{ld}}{h_l-h_{ld}})& (1)\end{array}$

where q_CHF is the predicted value of CHF calculated by the model, h_fg is the latent heat of evaporation, x₂ is the average quality at the bubble layer and x₁ 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₃ is the lateral mass velocity from the core to the bubble layer due to turbulence, which is determined by:

$\begin{array}{ll}{G}_3=G{i}_b\psi & (2)\end{array}$

in which G₃ 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:

$\begin{array}{ll}\psi =\left\{\frac{1}{\sqrt{2\pi }}\text{exp}\left[-\frac{1}{2}{(\frac{{\text{v}}_{\text{1}l}}{{\sigma }_{v^{\prime }}})}^{\text{2}}\right]-\frac{1}{2}(\frac{{\text{v}}_{\text{1}l}}{{\sigma }_{v^{\prime }}})erfc(\frac{\text{1}}{\sqrt{\text{2}}}\frac{{\text{v}}_{\text{1}l}}{{\sigma }_{v^{\prime }}})\right\}& (3)\end{array}$

where v' is the radial fluctuating velocity, ${\sigma }_{{\text{v}}^{\text{′}}}$ is the standard deviation of v, and v₁l is the radial velocity created by vapor generation.

$\begin{array}{ll}{v}_{1l}=\frac{q_b}{h_{fg}{\rho }_g}& (4)\end{array}$

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:

$\begin{array}{ll}{q}_b=q\left[\frac{h_l-h_{ld}}{h_f-h_{ld}}\right]& (5)\end{array}$

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$\left[\sqrt{\overline{{(v^{\prime })}^2}}/U_{\tau }\right]/(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:

$\begin{array}{ll}\frac{\sqrt{\overline{{(v^{\prime })}^2}}}{U_{\tau }}/(\frac{l_e}{r_0})=2.318{\left(\frac{r_0}{y}\right)}^{0.538}& (6)\end{array}$

where U_τ is the frictional velocity, l_e is the Prandtl mixing length, r₀ 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):

$\begin{array}{ll}\frac{\sqrt{\overline{{(v^{\prime })}^2}}}{U_{\tau }}=2.318F_1{\left(\frac{r_0}{y}\right)}^{0.538}\left(\frac{l_e}{r_0}\right)& (7)\end{array}$

FIGURE 2

Figure 2. Comparisons of radial turbulence fluctuations and fitted representations.

And the general relationship of l_e is:

$\begin{array}{ll}{l}_e=0.4y& (8)\end{array}$

The frictional velocity U_τ is defined as:

$\begin{array}{ll}{U}_{\tau }\text{=}{\left(\frac{f}{2}\right)}^{1/2}\frac{G}{\rho }& (9)\end{array}$

where f in the turbulent region is:

$\begin{array}{ll}f=0.046{\mathrm{Re}}^{-0.2}& (10)\end{array}$

From Equations (6) to (10), we have:

$\begin{array}{ll}\frac{\sqrt{\overline{{(v^{\prime })}^2}}}{G}=\frac{0.176}{\rho }{F}_1{\mathrm{Re}}^{-0.1}{\left(\frac{y}{r_0}\right)}^{0.538}& (11)\end{array}$

Weisman and Pei (1983) assumed the distance from the wall at which the bubbly layer–core interface occurs, y_c, is:

$\begin{array}{ll}{y}_c=\left(k{D}_p\right)/(0.4F_2)& (12)\end{array}$

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:

$\begin{array}{ll}\sqrt{\overline{{(v^{\prime })}^2}}\left(\frac{\rho }{G}\right)=0.2962\left(\frac{F_1}{F_2^{0.6}}\right){(k)}^{0.462}{\mathrm{Re}}^{-0.1}\times {\left(\frac{D_p}{D}\right)}^{0.462}=i_b& (13)\end{array}$

As indicated by Weisman and Pei, the two-phase factor F₁ represents the effect of bubble motion on turbulence intensity, and the two-phase factor F₂ accounts for the effect of bubble motion on the fluctuation velocity. Weisman and Pei assumed the ratio F₁/F₂ to be in the form:

$\begin{array}{ll}(F_1/F_2^{0.6})=\left[1+a\left(\frac{{\rho }_l-{\rho }_g}{{\rho }_g}\right)\right]& (14)\end{array}$

The physical meaning of Equation (14) is the effect of two-phase flow on turbulence intensity.

Constitutive Relations

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:

$\begin{array}{ll}{h}_{ld}=h_f-\Delta T_{ld}{C}_{pf}& (15)\end{array}$

$\begin{array}{ll}\{\begin{array}{l}0\le y_b^{+}\le 5.0;\Delta T_{ld}=\frac{q}{H_{db}}-\frac{q}{G{\left(\frac{f}{8}\right)}^{1/2}}\mathrm{Pr}{y}_b^{+}\hfill \\ 5.0\le y_b^{+}\le 30.0;\Delta T_{ld}=\frac{q}{H_{db}}-5.0\frac{q}{G{\left(\frac{f}{8}\right)}^{1/2}}\left\{\mathrm{Pr}+In\left[1+\mathrm{Pr}(y_b^{+}/5.0-1.0)\right]\right\}\hfill \\ y_b^{+}\ge 30.0;\Delta T_{ld}=\frac{q}{H_{db}}-5.0\frac{q}{G{\left(\frac{f}{8}\right)}^{1/2}}\left\{\mathrm{Pr}+In(1.0+5.0\mathrm{Pr})+0.5In(y_b^{+}/30.0)\right\}\hfill \end{array}& (16)\end{array}$

where,

$\begin{array}{ll}{y}_b^{+}=\frac{0.015}{{\mu }_f}\sqrt{\frac{\sigma g_{c}D}{{\nu }_f}}& (17)\end{array}$

and H_db is the single-phase heat transfer coefficient given by the Dittus-Boelter equation:

$\begin{array}{ll}{H}_{db}=0.023{\mathrm{Re}}^{0.8}{\mathrm{Pr}}^{0.4}{k}_f/D& (18)\end{array}$

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:

$\begin{array}{ll}{D}_p=0.015\left(\frac{8\rho \sigma D}{f{G}^2}\right)& (19)\end{array}$

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.

FIGURE 3

Figure 3. Enhancement effect of grids on CHF.

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 hydraulic-equivalent diameter from the grid. Nagayoshi and Nishida (1998) expressed the change in radial turbulence intensity factor as:

$\begin{array}{ll}{\left(\frac{\sqrt{\overline{{(v^{\prime })}^2}}}{U_{\tau }}\right)}_{grid}/\left(\frac{\sqrt{\overline{{(v^{\prime })}^2}}}{U_{\tau }}\right)=1+6.5{\epsilon }^2{e}^{-0.27\frac{x}{D}}& (20)\end{array}$

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:

$\begin{array}{ll}\{\begin{array}{l}{i}_{b-use}=i_b{F}_{grid}\hfill \\ F_{grid}=1+a_2{\epsilon }^2{\text{e}}^{-b_2\frac{x}{D}}\hfill \end{array}& (21)\end{array}$

The coefficients a₂ and b₂ depend on the type of grid. The experimental data of Yao et al. (1982) indicated that the heat-transfer 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₂ has been established as−0.13. A value of 6.5 has been derived for the coefficient a₂ 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₁/F₂), the physical meaning of which is the effect of boiling or two phases on turbulence intensity, as shown in Equation (14), rewritten here as:

$\begin{array}{ll}(F_1/F_2^{0.6})=\left[1+a\left(\frac{{\rho }_l-{\rho }_g}{{\rho }_g}\right)\right]& (22)\end{array}$

In Weisman and Pei's study, the coefficient “a” to account for the enhancement in turbulence intensity for two-phase flow, is expressed as:

$\begin{array}{ll}a=0.135{\left(\frac{G}{G_{cr}}\right)}^{a_1}& (23)\end{array}$

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:

$\begin{array}{ll}(F_1/F_2^{0.6})=1+a_3{\left(\frac{G}{G_{cr}}\right)}^{a_1}{\alpha }^{b_1}& (24)\end{array}$

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₁ equals 1.7 and a₃ equals 197.9 when α is greater than or equal to 0.52, and b₁ equals −0.06 and a₃ equals 69.65 when α is smaller than 0.52. For coefficient a₁, 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 α₂ 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:

$\begin{array}{ll}{\rho }_2=(1-{\alpha }_2){\rho }_f+{\alpha }_2{\rho }_g& (25)\end{array}$

$\begin{array}{ll}{\rho }_1=\frac{r^2}{{(r-s)}^2}\rho -{\rho }_2\frac{2(r-\frac{s}{2})s}{{(r-s)}^2}& (26)\end{array}$

$\begin{array}{ll}{\alpha }_1=\frac{{\rho }_l-{\rho }_1}{{\rho }_l-{\rho }_g}& (27)\end{array}$

$\begin{array}{ll}{x}_1=\frac{{\alpha }_1\times {\rho }_g}{{\rho }_1}& (28)\end{array}$

$\begin{array}{ll}{x}_2=\frac{{\alpha }_2\times {\rho }_g}{{\rho }_2}& (29)\end{array}$

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.

Assessment and Analysis of the Theoretical CHF Model

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.

FIGURE 4

Figure 4. Test bundle setup in the NPIC 5 × 5 rod-bundle experiment (rods setup).

FIGURE 5

Figure 5. Test bundle setup in the NPIC 5 × 5 rod-bundle experiment (grids setup).

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.

TABLE 1

Table 1. Range of CHF data selected in the assessment.

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:

$\begin{array}{ll}DNBR=\frac{q_{pre}}{q_{\exp }}& (30)\end{array}$

$\begin{array}{ll}{R}_{av}=\frac{1}{N}{\displaystyle \sum _{i=1}^{N}MDNB{R}_i}& (31)\end{array}$

$\begin{array}{ll}SD={\left[\frac{1}{N}{\displaystyle \sum _{i=1}^N{(MDNB{R}_i-R_{av})}^2}\right]}^{1/2}& (32)\end{array}$

$\begin{array}{ll}RMS={\left[\frac{1}{N}{\displaystyle \sum _{i=1}^N{(MDNB{R}_i-1)}^2}\right]}^{1/2}& (33)\end{array}$

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.

FIGURE 6

Figure 6. Comparison of predicted and experimental CHF values by 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).

FIGURE 7

Figure 7. Variation of predicted MDNBRs with pressure.

FIGURE 8

Figure 8. Variation of predicted MDNBRs with mass flux.

FIGURE 9

Figure 9. Variation of Predicted MDNBRs with quality.

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).

Conflict of Interest

QY was employed by the company Guangxi Fangchenggang Nuclear Power Co., Ltd. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Bruder, M., Bloch, G., and Sattelmayer, T. (2015). Critical heat flux in flow boiling-review of the current understanding and experimental approaches. Heat Transf. Eng. 38, 347–360. doi: 10.1080/01457632.2016.1189274

CrossRef Full Text | Google Scholar

Celata, G. P., Cumo, M., Mariani, A., Simoncini, M., and Zummo, G. (1994). Rationalization of existing mechanistic models for the prediction of water subcooled flow boiling critical heat flux. Int. J. Heat Mass Transf. 37, 347–360. doi: 10.1016/0017-9310(94)90035-3

CrossRef Full Text | Google Scholar

Chang, S. H., and Lee, K. W. (1989). A critical heat flux model based on mass, energy, and momentum balance for upflow boiling at low qualities. Nucl. Eng. Des. 113, 35–50. doi: 10.1016/0029-5493(89)90294-X

CrossRef Full Text | Google Scholar

Kinoshita, H., Nariai, H., and Inasaka, F. (2001). Study of critical heat flux mechanism in flow boiling using bubble crowding model. JSME Int. J. 44, 81–89. doi: 10.1299/jsmeb.44.81

CrossRef Full Text | Google Scholar

Kodama, S., and Kataoka, I. (2002). Study on analytical prediction of forced convective CHF in the wide range of quality. Int. Conf. Nucl. Eng. 62, 147–154. doi: 10.1115/ICONE10-22128

CrossRef Full Text | Google Scholar

Kwon, Y. M., and Chang, S. H. (1999). A mechanistic critical heat flux model for wide range of subcooled and low-quality flow boiling. Nucl. Eng. Des. 188, 27–47. doi: 10.1016/S0029-5493(99)00025-4

CrossRef Full Text | Google Scholar

Lahey, R. T., and Moody, F. (1977). The Thermal Hydraulics of a Boiling Water Reactor. La Grange Park: American Nuclear Society.

Google Scholar

Laufer, J. (1953). The structure of turbulence in fully developed pipe flow. NACA Technical Note No. 2954.

Google Scholar

Lee, S. C., Dorra, H., and Bankoff, S. G. (1992). A critical review of predictive models for the onset of significant void in forced-convection subcooled boiling, HTD, fundamentals of Subcooled Flow Boiling, ASME 217, 33–39. doi: 10.2172/10104504

CrossRef Full Text

Lee, S. L., and Durst, F. (1980). On the motions of particles in turbulent flow, U.S. Nuclear Regulatory Commission Report NUREG/CR-1556.

Google Scholar

Levy, S. (1966). Forced convection subcooled boiling - prediction of Vapor volumetric fraction. Int.J. Heat Mass Transf. 10, 91–965. doi: 10.1016/0017-9310(67)90071-3

Lim, J. C. (1988). Improvements in the theoretical prediction of the departure from nucleate boiling (Ph.D. thesis). University of Cincinnati, Cincinnati, OH, United States.

Google Scholar

Liu, W., Bai, N., Zhu, Y-B., Shan, J-Q., Zhang, B., Gou, J.-L, et al. (2014). Study on reactor thermal-hydraulic sub-channel analysis code ATHAS. Nucl. Sci. Eng. 34, 59–66.

Google Scholar

Nagayoshi, T., and Nishida, K. (1998). Spacer effect model for subchannel analysis. J. Nucl. Sci. Technol. 35, 399–405. doi: 10.1080/18811248.1998.9733881

CrossRef Full Text | Google Scholar

Trupp, A. C., and Azad, R. S. (1973). The structure of turbulent flow in triangular array rod bundles. Nucl. Eng. Des. 32, 47–84. doi: 10.1016/0029-5493(75)90090-4

CrossRef Full Text | Google Scholar

Venkateswararao, P., Semiat, R., and Dukler, A. E. (1982). Flow pattern transition for gas-liquid flow in a vertical rod bundle. Int. J. Multiphase Flow 8, 509–524. doi: 10.1016/0301-9322(82)90021-0

CrossRef Full Text | Google Scholar

Weisman, J. (1991). The current status of theoretically based approaches to the prediction of the critical heat flux in flow boiling. Nucl. Technol. 9, 1–21. doi: 10.13182/NT92-A34699

CrossRef Full Text | Google Scholar

Weisman, J., and Illeslamlou, S. (1988). Phenomenological model for prediction of critical heat flux under highly subcooled conditions. Fusion Technol. 13, 654–659. doi: 10.13182/FST88-A25140

CrossRef Full Text | Google Scholar

Weisman, J., and Pei, B. (1983). Prediction of critical heat flux in flow boiling at low qualities. Int. J. Heat Mass Transf. 10, 1463–1477. doi: 10.1016/S0017-9310(83)80047-7

CrossRef Full Text | Google Scholar

Yao, S. C., Hochreiter, L. E., and Leech, W. J. (1982). Heat-transfer augmentation in rod bundles near grid spacers. Int. J. Heat Mass Transf. 104, 76–81. doi: 10.1115/1.3245071

CrossRef Full Text

Ying, S. H., and Weisman, J. (1986). Prediction of critical heat flux in flow boiling at intermediate qualities. Int. J. Heat Mass Transf. 29, 1639–1648. doi: 10.1016/0017-9310(86)90105-5

CrossRef Full Text | Google Scholar

Nomenclature

Appendix A: Evaluation Procedure Of The Mechanistic Bubble-Crowding Model

The calculation procedure is as follows:

$q_{pre}=G{h}_{fg}(x_2-x_1)\psi i_{b-use}\frac{h_f-h_{ld}}{h_l-h_{ld}}$

Calculate the bubble departure enthalpy h_ld

$\begin{array}{l}{h}_{ld}=h_f-\Delta T_{ld}{C}_{pf}\\ \{\begin{array}{l}0\le y_b^{+}\le 5.0;\Delta T_{ld}=\frac{q}{H_{db}}-\frac{q}{G{(\frac{f}{8})}^{1/2}}\mathrm{Pr}{y}_b^{+}\hfill \\ 5.0\le y_b^{+}\le 30.0;\Delta T_{ld}=\frac{q}{H_{db}}-5.0\frac{q}{G{(\frac{f}{8})}^{1/2}}\left\{\mathrm{Pr}+In\left[1+\mathrm{Pr}(y_b^{+}/5.0-1.0)\right]\right\}\hfill \\ y_b^{+}\ge 30.0;\Delta T_{ld}=\frac{q}{H_{db}}-5.0\frac{q}{G{(\frac{f}{8})}^{1/2}}\left\{\mathrm{Pr}+In(1.0+5.0\mathrm{Pr})+0.5In(y_b^{+}/30.0)\right\}\hfill \end{array}\end{array}$

where,

$y_b^{+}=\frac{0.015}{{\mu }_f}\sqrt{\frac{\sigma g_{c}D}{{\nu }_f}}$

$H_{db}=0.023{\mathrm{Re}}^{0.8}{\mathrm{Pr}}^{0.4}{k}_f/D$

$\mathrm{Re}=\frac{GD}{{\eta }_f},\mathrm{Pr}=\frac{c_{pf}{\eta }_f}{k_f}$

The turbulence intensity factor i_b is calculated as

$i_b=0.2962\times {(k)}^{0.462}{\mathrm{Re}}^{-0.1}\times {(\frac{D_p}{D})}^{0.462}[1+a_3{(\frac{G}{G_{\text{cr}}})}^{a_1}{\alpha }^{b_1}]$

where k equals 0.58 and

$\begin{array}{l}{a}_1=0.3,G>G_{\text{c}};a_1=-0.6,G<G_c;b_1=1.7,\\ =a_{3}197.9,\alpha \ge 0.52;b_1=-0.06,a_3=69.65,\alpha <0.52;\end{array}$

$\begin{array}{l}{F}_{grid}=1+a_2{\epsilon }^2{e}^{-b_2\frac{x}{D}}\\ i_{b-use}=i_b{F}_{grid}\end{array}$

where a₂ is 6.5, and b₂ is −0.13.

$D_p=0.015\left(\frac{8\rho \sigma D}{f{G}^2}\right)$

Calculation of function ψ:

$\psi =\left\{\frac{1}{\sqrt{2\pi }}\exp \left[-\frac{1}{2}{(\frac{v_{1l}}{{\sigma }_{v^{\prime }}})}^2\right]-\frac{1}{2}(\frac{v_{1l}}{{\sigma }_{v^{\prime }}})erfc(\frac{1}{\sqrt{2}}\frac{v_{1l}}{{\sigma }_{v^{\prime }}})\right\}$

$\begin{array}{l}{v}_{1l}=\frac{q_b}{{\rho }_g{h}_{fg}}\\ q_b=\frac{h_l-h_{ld}}{h_f-h_{ld}}{q}_{\exp }\\ {\sigma }_{v^{\prime }}=\frac{i_{b-use}G}{\overline{\rho }}=i_{b-use}\left[\frac{G}{{\rho }_l}+(\frac{G}{{\rho }_g}-\frac{G}{{\rho }_l})x_{avg}\right]\end{array}$

To calculate the x₁ and ρ₁ based on α₂ = 0.82:

$\begin{array}{l}s=5.5D_p\\ {\rho }_2=(1-{\alpha }_2){\rho }_f+{\alpha }_2{\rho }_g\\ {\rho }_1=\frac{r^2}{{(r-s)}^2}{\rho }_{avg}-{\rho }_2\frac{2(r-\frac{s}{2})s}{{(r-s)}^2}\\ {\alpha }_1=\frac{{\rho }_l-{\rho }_1}{{\rho }_l-{\rho }_g}\\ x_1=\frac{{\alpha }_1\times {\rho }_g}{{\rho }_1},x_2=\frac{{\alpha }_2\times {\rho }_g}{{\rho }_2}\end{array}$

The predicted value is the result of the calculation on the left side of the first equation in Appendix A.