Prediction of Critical Heat Flux for Subcooled Flow Boiling in Annulus and Transient Surface Temperature Change at CHF

: The ability to predict critical heat ﬂux (CHF) is of considerable interest for high-heat equipment, including nuclear reactors. CHF prediction from a mechanistic model for subcooled ﬂow boiling in rod bundles still remains unsolved. In this paper, we try to predict the CHF in an annulus, which is the most basic ﬂow geometry simpliﬁed from a fuel bundle, using a liquid sublayer dryout model. The prediction is validated with both water and R113 data, showing an accuracy within ± 30%. After the CHF in an annulus is calculated successfully, a near-wall vapor–liquid structure is proposed on the basis of the liquid sublayer dryout model. Modeling of heat transfer modes over the heating surface at CHF is performed, and predictions of the changes in liquid sublayer thickness and heater surface temperature at the CHF occurrence point are carried out by solving the heat conduction equation in cylindrical coordinates with a convective boundary condition, which changes with the change in ﬂow pattern over the heating surface. Transient changes in the liquid sublayer thickness and surface temperature at the CHF occurrence point are reported.


Introduction
The thermal output of high-heat equipment using forced flow subcooled boiling as a cooling method is limited by its cooling limit, the so-called critical heat flux (CHF).For subcooled flow boiling, CHF is closely related to the departure from nucleate boiling (DNB).DNB is the starting point of transition boiling and film boiling, which causes a significant increase in the surface temperature.The ability to predict CHF is of considerable interest for high-heat equipment, including nuclear reactors.
CHF prediction methods can be categorized into three types: CHF correlation, CHF look-up table (LUT), and mechanistic model.For fuel assembly geometries that represent a typical PWR, specific CHF correlations have been developed.The correlations [1,2] can predict CHF with high accuracy as they are developed by fitting experimental data.However, the applicable ranges of correlations are limited by the experimental database.
The LUTs are developed for easy use in CHF prediction.Standard tables of the CHF for pressure, mass flux, and subcooling were developed by the USSR Academy of Science [3] and Chalk River Nuclear Laboratories [4][5][6][7].The tables present CHF values at discrete ranges of pressure, mass flux, and quality for 8 mm tubes.A correlation is used to account for the diameter effect and to extend the application to other tube diameters.In conditions where data are scarce or unavailable, CHF values are obtained by extrapolation using empirical correlations of available data.
Mechanistic models have the advantage, with respect to the correlations, of being able to characterize not only the existing and developing database, but also being able to be used to predict the CHF beyond the validated database.In future thermal design of nuclear reactors, it is expected that each boiling heat transfer mode, including DNB, can be modeled on the basis of their physical phenomena, and the CHF can be predicted when an excursion of wall temperature is detected.Therefore, it is needed to understand the CHF triggering mechanism, construct the CHF prediction method on the basis of the triggering mechanism, and predict the transient wall temperature change at the CHF.
So far, CHF predictions from mechanistic models for subcooled flow boiling have mainly been focused on circular tubes, and the predictions in rod bundles of nuclear reactors still remain unsolved.As the first step for the CHF prediction in rod bundles, in this study, we try to predict the CHF in an annulus, which is the most basic flow geometry simplified from a fuel bundle, using a liquid sublayer dryout model.After the CHF in the annulus is predicted successfully from the from mechanistic model, a near-wall vaporliquid structure at the CHF is proposed, and modeling of heat transfer modes over the heating surface is performed.Then, predictions of the transient surface temperature change at CHF are carried out by solving the heat conduction equation.The transient changes of liquid sublayer thickness and surface temperature at CHF occurrence are reported.

Liquid Sublayer Dryout Model and Its Application to an Annulus
Although modeling of the CHF for subcooled flow boiling can be categorized into six groups [8], only the bubble crowding [9][10][11] and liquid sublayer dryout models [12][13][14][15][16] are currently receiving attention.The bubble crowding model assumes a bubbly layer existing adjacent to the heating wall and CHF occurring when the void fraction of the bubbly layer reaches a critical value.Although the model has been reported with some good prediction results, it is quite empirical because it employs the determination of the turbulent exchange in the bubbly layer and liquid bulk region.
As shown in Figure 1, the liquid sublayer dryout model assumes that an elongated vapor clot, which is formed as a consequence of coalescence of small bubbles rising along the near-wall region, is overlying a very thin liquid sublayer adjacent to heater wall.The CHF is assumed to occur at the complete dryout of liquid sublayer during the passage time of the vapor clot.The liquid sublayer dryout model seems to be promising from a number of observational experiments [17][18][19].As a result, CHF is described as where u B , L B , and δ 0 are the vapor clot velocity, vapor clot length, and initial thickness of the liquid sublayer, respectively.L B /u B represents the passage time of the vapor clot, ρ l is the liquid density, and h lg is the latent heat.
Conventionally, an instability wave at the interface between the liquid sublayer and the vapor clot is considered, and the length of the vapor clot is assumed to be determined by the Helmholz instability wavelength, which is inversely proportional to u 2 B .
Then, the CHF is written as Therefore, the liquid sublayer dryout model relies on the calculations of u B and δ 0 .Different models employ different approaches.C is a coefficient which is only a function of the physical properties.
Liu et al. [16] considered the Helmholz instability waves at the interface between the liquid sublayer and the vapor clot, as well as at the interface between the vapor clot and liquid bulk; they proved that the two instability wavelengths at the two interfaces were equal to each other.The Helmholz instability wavelengths at the two interfaces are expressed in Equations ( 4) and (5), respectively.With the assumption that the length of the vapor clot is determined by the Helmholz instability wavelength, we get Equation (6).
At the interface between the liquid sublayer and the vapor clot, At the interface between the vapor clot and main flow, where λ 1 and λ 2 are the Helmholz instability wavelengths at the interface between the liquid sublayer and the vapor clot and at the interface between the vapor clot and the liquid bulk, respectively.u ls is the velocity of the liquid sublayer.As the initial liquid film thickness δ 0 is on the order of several to ten microns, the liquid sublayer velocity u ls can be approximated to be zero.ρ v is the vapor density, σ is the surface tension, and ρ c and u c are the average density and velocity of liquid bulk, respectively.From Equation ( 6), the vapor clot velocity u B can be obtained using Equation (7).
where G is the mass flow rate, and α is the average void fraction at CHF elevation.Liu [16] recommended to use Ahmad model [20] is to calculate the void fraction.For the case that water is used as a coolant, the Levy model [21] is recommended to calculate the net vapor generation point (NVG), which is necessary in the calculation of the void fraction.As shown in Figure 2, by considering the force balance in the flow direction of the vapor clot, i.e., where the drag force F D is balanced by the buoyancy force F B , the liquid velocity at the centerline of the vapor clot is calculated using Equation (9).
where u BL is the liquid velocity at the centerline of the vapor clot, and C D is drag coefficient, which is calculated from the Harmathy model [22].After u BL is calculated, the distance from the wall to the centerline of the vapor clot can be calculated from the velocity profile in the liquid phase.For the case of an annulus, the velocity distribution is calculated from the Nouri correlation [23] developed for an annulus, which is shown in Equation (11).
where u l is local liquid velocity, U τ is the friction velocity, µ l is the liquid viscosity, y + and u + l are the dimensionless distance from the tube wall and the velocity, respectively, τ w is the wall shear stress, f is the friction factor, and r 2 is the inner diameter of the annulus flow channel.The friction factor f for an annulus is calculated from Equation ( 13), as recommended by Nouri [23].
where R e is the Reynolds number, G is the mass flux, and D is the equivalent diameter of the test section.
By substituting the calculated u BL into Equation ( 11), the distance from the wall to the centerline of the vapor clot, y, can be calculated.Then, the initial thickness of the liquid sublayer δ 0 , as shown in Figure 2, is obtained as follows: where D B is the thickness of the vapor clot, which is calculated from Levy model [21] for the case of water in the Liu model [16].

Prediction of CHF in an Annulus
For given geometric condition (equivalent diameter of the test section is D and heating length of the test section is L) and thermal hydraulic conditions (mass flux G, system pressure P, and inlet subcooling ∆T in ), the CHF can be predicted by an iterative procedure through the abovementioned equations.Figure 3 shows the CHF calculation flow chart.More detailed information for the calculation can be found in [16].The above model applied to an annulus was validated using Fiori [24] and Hino [25] data.The Fiori experiment used water as a coolant.The inner diameter of the annulus was 0.312 inches, the pitch between the inner and outer walls was 0.095 to 0.112 inches, and the heating length was 10 inches.A total of 16 CHF data points were acquired at pressures of 0.2-0.6MPa, mass flow rates of 800-1700 kg/(m 2 s), and inlet subcooling of 64-120 K.The Hino experiment used R113 as a coolant.As shown in Figure 4, the Hino test channel was an annulus with an inner heater rod whose outer diameter was 8 mm and an outer pipe whose inner diameter was 18 mm.The inner heater rod was a circular annulus with an inner diameter of 7 mm and a thickness of 0.5 mm.The heating length was 400 mm, and the heater material was stainless steel.A total of 6 CHF data points were obtained at a pressure of 0.147 MPa, mass flow rates of 1239 and 512 kg/(m 2 s), and inlet subcooling rates of 30, 20, and 10 K.The prediction results for the experimental data are shown in Figures 5 and 6, respectively.Note that the Saha-Zuber model [26], which was recommended by Hino [25], is used to calculate the NVG point for Hino's data.The prediction accuracy, which is defined as q CHF, cal /q CHF, exp , revealed a mean of 1.06 and a standard deviation of 14.7% for Fiori's data.For Hino's data, relatively good results were obtained at a high flow rate and high subcooling degree, but the prediction accuracy decreased for a low flow rate and low subcooling degree, with CHF generally predicted in the range of ±30%.From Figure 6, it can be observed that the present model showed better prediction results for comparatively high mass flux and high subcooling degree, where CHF is nearer to the DNB type.The present model was developed for a DNB-type CHF.Although the model itself creates no empirical constants, it uses correlations in the CHF calculation process for the friction factor, liquid velocity distribution, drag coefficient, NVG point, void fraction, etc.The correlations developed for high mass flux and high subcooling degree were selected in this paper.Therefore, the present model is more adaptable to high mass flux and high subcooling conditions.The present model was only validated using the Fiori [24] and Hino [25] data.More validations are needed in the future to check its adaptability range.16), where k l , ρ l , and T l are the thermal conductivity, density, and temperature of the liquid sublayer, respectively.Because the liquid sublayer is very thin and adjacent to the heater wall, T l can be taken as the saturation temperature.The CHF occurs when the liquid sublayer completely dries out.As a result, in the region where the liquid sublayer is present (L ls ), heat is transferred through the thin liquid film, whereas, in the region where the liquid film has been completely evaporated, L dry , heat is transferred through steam vapor.Assuming that the bubble distribution across the flow channel is uniform, Podowski et al. [27] calculated the nucleate boiling region length L NB from Equation (17), where α is the average void fraction.
Figure 7. Near-wall vapor-liquid structure at the CHF based on liquid sublayer dryout model.

Modeling of Heat Transfer Modes over Heated Surface near the CHF
The heater temperature T can be calculated by solving the thermal conduction equation of the heater in cylindrical coordinates, which is expressed in Equation ( 18), where q is the volumetric heating density, while ρ w , C pw , and k w are the density, specific heat, and heat conduction of the heater.
The initial condition inside the heater is written in Equation ( 19), which is derived by solving the thermal conduction equation (Equation ( 18)) when ∂T ∂t = 0, with a surface wall temperature calculated from the Jens-Lottes correlation (Equation ( 20)) for a fully developed nucleate boiling region.q w and P in Equation ( 20) are the wall surface heat flux and system pressure, respectively.The unit for the system pressure P in Equation ( 20) is MPa.∆T nucl is the wall surface superheat.
The boundary condition for the heater inner surface is expressed as For the heater outer surface, the boundary condition changes with the change in heat transfer mode.
In the nucleate boiling region (L NB ), the boundary condition at the heater outer surface is expressed in Equation ( 22).q w is calculated from the Shaver [28] heat transfer model for subcooled flow boiling, which is shown in Equation (23).
where q b and q 1φ are the heat fluxes transferred by the boiling and liquid single phases, respectively, A b is the ratio of the area where boiling occurs on the heating surface, ∆h l and ∆T l denote subcooling of the liquid phase, h 1p is the liquid single-phase heat transfer coefficient, which is calculated from the Dittus-Bolter equation, ∆T w and ∆T 0 are the surface superheat and the minimum wall superheat necessary for boiling, respectively, ∆T nucl is the wall superheat in the fully developed nucleate boiling region (Equation( 20)), υ lg is the difference in the specific volume between the vapor phase and liquid phase, and σ is the surface tension.
In the region under the vapor clot where the liquid sublayer is present (L ls ), the boundary condition at the heater outer surface is expressed in Equation ( 24), where T w is the heater surface temperature, and h l is the liquid heat transfer coefficient.The calculation to h l is found in Equation (24), where k l and δ are the thermal conductivity and the thickness of the liquid sublayer, respectively.δ decreases due to evaporation and is calculated from Equation ( 16).T l is the liquid sublayer temperature and is taken as the saturation temperature.
In the region under the vapor clot where no liquid film exists (L dry ), heat is transferred through the steam vapor film, and the heat transfer coefficient h g can be calculated from Equation (25), where k g is the thermal conductivity of steam vapor, and D B is the thickness of the vapor clot.T g is the temperature of the vapor clot and is taken as the saturation temperature.

Analysis Case
One case from the Hino data [25] was used as the analysis case, which is shown in Table 1.A cross-section view of the heater and flow area in the Hino CHF test is shown in Figure 4. Table 2 shows the predicted CHF for the analysis case, which was reported in Section 3, as well as detailed information including the initial liquid sublayer thickness, as well as the length and velocity of the vapor clot, calculated from the liquid sublayer dryout model.As shown in Figure 7, around the CHF occurrence point (generally channel exit), the heater surface is periodically covered with a nucleate boiling region (L NB ) and an elongated vapor clot hovered region (L B ).The vapor clot hovering region can be further divided into the liquid sublayer existing region and liquid sublayer exhausted region.The liquid sublayer exhausted region is the CHF location.If we assume that the nucleate boiling region flows as the same speed as the vapor clot, a period for one cycle of the near-wall phase structure is 3.731 ms, of which the vapor clot passage time is 2.709 ms and the nucleate boiling passage time is 1.022 ms.The conditions shown in Table 2 were used as the calculation conditions for the transient heater surface temperature change at the CHF.In order to obtain a liquid sublayer dryout area, the volumetric heating density was set to a value corresponding to 1.01 times the calculated CHF.Table 3 shows some other conditions needed in the calculation.Figure 8 shows the calculated liquid sublayer thickness change under the elongated vapor clot at CHF.If we assume that the time at which the elongated vapor clot first arrives at the channel exit (CHF occurrence point) is 0 ms, the liquid sublayer completely evaporates after 1.739 ms, and a dry spot forms, which results in the vapor clot contacting the heater surface directly.For simplicity, in this study, the velocity of the vapor clot u B remains unchanged even after the dryout of the liquid sublayer.Therefore, the elongated vapor clot flows away, and the heater surface is recovered by the nucleate flow boiling, even after the vapor clot contacts the heater surface.
Figure 9 shows the transient change in heater surface temperature.When the liquid sublayer thickness is evaporated and approaching 0, the heater surface temperature decreases to nearly the saturation temperature.Thereafter, as the liquid sublayer is completely dried out, the heater surface is covered by vapor steam and the temperature rises; then, the temperature drops due to the heater surface being recovered with the nucleate boiling flow.
In the presented calculation results, there is no temperature excursion observed.The main reason is the modeling of the speed of the elongated vapor clot after the dryout.In the presented calculation, the speed of elongated vapor clot after the dryout was calculated to be the same speed as that before the dryout, without considering the drag from the heater wall, even though the vapor clot was in direct contact with the wall.The drag from the heater wall results in a longer contact time of the vapor clot with the heater surface, leading to a higher wall surface temperature excursion.In order to reproduce the temperature excursion observed in experiments, it is necessary to model the vapor clot dynamics more accurately after the dryout of the liquid sublayer, i.e., when the vapor clot is in direct contact with the wall.

Concluding Remarks
The ability to predict the CHF is of considerable interest for high-heat equipment, including nuclear reactors.As the first step for the CHF prediction from the mechanistic model in rod bundles, in this study, we predicted the CHF in the annulus, which is the most basic flow geometry simplified from a fuel bundle, using a liquid sublayer dryout model.The Nouri liquid-phase velocity distribution correlation was employed to apply the model to an annulus.The prediction was validated using water and R113 data.The results show that the CHF in an annulus can be predicted with an accuracy of about ±30%.
For the future thermal design of nuclear reactors, it is expected that each boiling heat transfer mode, including the DNB, can be modeled on the basis of their physical phenomena; then, the CHF can be predicted when an excursion of wall temperature is calculated.Therefore, after the CHF in an annulus was successfully calculated from the mechanistic model, a near-wall vapor-liquid structure at CHF was proposed on the basis of the liquid sublayer dryout model.Then, the heat transfer modes over the heating surface were modeled.The predictions of the transient surface temperature at the CHF were carried out by solving the heat conduction equations in cylindrical coordinates with a convective boundary condition, which changed with the change in heat transfer modes over the heated surface.The transient changes in liquid sublayer thickness and wall surface temperature at the CHF were reported.In the presented calculation, no temperature excursion leading to heater burnout was observed.One of the reasons was the modeling of the speed of the elongated vapor clot after the dryout.In the presented calculation, the speed of the elongated vapor clot after the dryout was calculated as the same speed as that before the dryout, without considering the drag from the heater wall, even though the vapor clot was in direct contact with the wall.The drag from the heater wall resulted in a longer contact time of the vapor clot with the heater surface, leading to a higher wall surface temperature.In order to reproduce the temperature excursion observed in experiments, it is necessary to model the vapor clot dynamics more accurately after dryout of the liquid sublayer, i.e., when the vapor clot is in direct contact with the wall.

Figure 4 .
Figure 4. Cross-section view of the heater and flow area in the Hino CHF test [25].

Figure 5 .
Figure 5. CHF prediction in an annulus using Fiori data (Levy model for NVG point).

Figure 6 .
Figure 6.CHF prediction in an annulus for Hino data (Saha−Zuber model for NVG point).

4 .
Predictions of the Changes in Liquid Sublayer Thickness and Heater Surface Temperature at the CHF 4.1.Modeling of Near Wall Vapor-Liquid Structure at the CHF On the basis of the liquid sublayer dryout model, a near-wall vapor-liquid structure at the CHF is proposed, which is shown in Figure 7. Around the CHF, the heated surface is periodically covered with a nucleate boiling region (L NB ) and an elongated vapor clot hovered region (L B ).The liquid sublayer trapped between the vapor clot and the heated surface has an initial thickness δ 0 and becomes thin due to evaporation over the passage time of the vapor clot.The change in thickness is calculated from Equation (

Figure 8 .
Figure 8. Calculated liquid sublayer thickness change at the channel exit (CHF point).

Figure 9 .
Figure 9. Calculated heater surface temperature change at channel exit (CHF occurrence point).

Table 2 .
Predicted CHF results for the analysis case.Passage time of the vapor clot τ B = L B /u B 2.709 ms Passage time of the boiling region τ NB = L NB /u B 1.022 ms 4.4.Calculation of Transient Heater Surface Temperature Change at CHF

Table 3 .
Other information used in the calculation of surface temperature change at the CHF.Volumetric heating density q corresponding to 1.01 q CHF,cal 7.224 × 10 5 kW/m 3 BL liquid velocity at the centerline of the vapor clot, [m/s] u l velocity profile inside liquid phase, [m/s] u ls velocity of the liquid sublayer, [m/s] u + u