Mathematical Model Predicting the Critical Heat Flux of Nuclear Reactors

Boiling heat transfer system keeps a nuclear power plant safe without getting over-heated. Crisis will occur if the dissipated heat flux exceeds the critical heat flux value. This study assumes the flow boiling system at high heat flux is characterized by the existence of a very thin liquid layer, known as the “sublayer”, which is trapped between the heated surface and the vapor blankets. In the present study, it is hypothesized that the heat transfer through the liquid sublayer is dominated by the heat conduction and the sublayer is dried out due to occurrence of Helmholtz instability as the relative velocity of the vapor blanket to the local liquid in the sublayer reaches a critical value. By recognizing this hypothesis, a theoretical model for low-quality flow is developed to predict boiling heat transfer and Critical Heat Flux (CHF). To verify the validity of the present model, the predictions are compared with the experimental data of flow boiling heat transfer in the simulation of Pressurized Water Reactor (PWR) and Boiling Water Reactor (BWR) conditions. For the PWR low-quality flow, the comparison demonstrates that the Helmholtz instability is the trigger condition for the onset of CHF.


INTRODUCTION
Subcooled or low-quality forced convective boiling system is utilized in the design of nuclear reactor cores, because its high heat fluxes can be dissipated. Therefore, to protect the heated surface from overheating or burnout, a reliable estimation of the boiling heat transfer and its limitation Critical Heat Flux (CHF) is extremely required in the flow boiling system. In general, the boiling heat transfer mechanisms are different at the low heat flux (q<0.6 CHF) and the high heat flux (q>0.6 CHF) as shown in Fig. 1. At the former condition, most of the heat flux is provided for nucleation of bubbles (Fig. 1a); while at the latter, the boiling heat transfer is dominated by the heat conduction in the liquid sublayer trapped between the heated surface and the vapor blanket (Fig. 1b).
In the review of high heat flux boiling mechanisms, many studies have found the existence of liquid sublayer between the heated surface and the bubble layer at high heat flux conditions. Based on the formation of the liquid sublayer, Katto and Yokoya (1968) developed a hydrodynamic model near the CHF according to the consumption of liquid film.
Through a measurement of the transient variation of heated surface temperature during nucleate pool boiling of water, Chi-Liang and Mesler (1977) confirmed the existence of a liquid film beneath the growing vapor which is paramount in transferring heat. Bhat et al. (1983) hypothesized that the heat transfer, in the high heat flux region between 0.6 CHF and CHF, takes place mainly due to the heat conduction through the liquid layer formed on the heated surface. In a subsequent paper of Bhat et al. (1986), experimental results of sublayer thickness and frequency of vapor mass show good agreement with their previous predictions. Later, many high heat flux models based on heat conduction in the liquid layer have been reported (Chinm et al., 1989;Jairajpuri and Saini, 1991;Rajvanshi et al., 1992). Although all the available models are developed specifically for pool boiling, Chinm et al. (1989) expected that the same concept of heat conduction through the sublayer may also be applied to flow boiling conditions. The sublayer concept with heat balance rather than heat conduction has also been applied in the research of flow boiling for predicting the CHF. For example, Haramura and Katto (1983), as ones of the pioneers of the sublayer dry out theory, originally derived a high heat flux model, based on the heat balance of liquid film located between the heated surface and a vapor blanket, to predict the flow boiling CHF on flat plates. Lee and Mudawwar (1988) postulated that CHF occurs when the liquid film underneath the vapor blanket dries out due to the oscillation motion of the blanket for an intearnal flow boiling system. They emphasized that the onset of sublayer dry out was triggered by Helmhotz instability in the sublayer-vapor blanket interface whiale as the length of the vapor blanket is equal to the Helmholtz critical wavelength. Later, by replacing the single-phase fluid properties with the two-phase homogeneous equilibrium fluid properties, Lin et al. (1988a;1988b) made an improved sublayer dry out model to extend the applicable range from the subcooled flow boiling to the low-quality saturated flow boiling. By accounting for the passage time of the blanket, Katto (1990a;1990b;1992) used the similar idea yet different approach to evaluate the DHF required to vaporize the liquid film underneath the blanket. Celata (1991) suggested that this mechanism may also be applied in the thermal-hydraulics studies of high heat flux, high mass flow rate fusion reactor.
Although thickness of liquid sublayers calculated by Lee and Mudawwar (1988) and Katto (1990a;1990b;1992) models are quite different, it is worth to note that the magnitude of sublayer thickness is extremely tiny (order of 10 −7~1 0 −3 m). It is thus appropriate to hypothesize that the convective boiling heat transfer is dominated by the heat conduction through the liquid sublayer at superheated conditions. Therefore, based on heat conduction in the sublayer, Lin et al. (1994) Science Publications JCS developed a theoretical model for subcooled flow boiling heat transfer at high heat flux conditions. Lee and Lin (1993) conducted the flow transient CHF experiments in the simulation of Pressurized Water Reactor (PWR) conditions. The experimental data were compared with the predictions of improved model by Lin et al. (1988a). Their results reveal that the sublayer dry out model is appropriate. All of the afore-mentioned theoretical approaches for predicting the flow boiling CHF (Lee and Mudawwar, 1988;Lin et al., 1988a;1988b;Katto, 1990a;1990b;1992) acknowledged that the sublayer dry out is triggered by Helmoholtz instability at the sublayer-vapor blanket interface. However, so far there is no theoretical model based on the Helmholtz instability concept proposed for prediction of occurrence of CHF. In the present study, a critical condition for onset of CHF is derived based on the Helmholtz instability at interface of two streams. The relative velocity of two streams is calculated using the heat transfer model by Lin et al. (1994). The predicted transient time to the onset of CHF provided reasonable agreement with the Lee and Lin (1993) experimental data.

MATERIALS AND METHODS
Based on the boiling configuration of subcooled flow through a vertical tube at the high heat flux conditions as shown in Fig. 1b, a theoretical boiling heat transfer model was proposed with the following assumptions: • In the subcooled flow through a vertical tube, the vapor blankets are formed from small bubbles pilling up as vertical distorted vapor cylinders. All the vapor blankets glide up parallel to the heated wall and shield the heated wall from the cooling of bulk flow. As a result, a liquid sublayer forms between the vapor blankets and the heated wall • For high heat flux, the passing period of the vapor blanket becomes sufficiently long as that reported in Hino and Ueda (1985) investigation. Thus, the sublayer seems always exist between the heated wall and the vapor blanket • Since the thickness of the sublayer is very thin with low flow rate, it is reasonable to assume that, in the sublayer, the bubble generation is ceased and the heat convection can be ignored due to small velocity, the heat transfer across the liquid sublayer is dominated by heat conduction Figure 2 shows the temperature distribution in the sublayer and the force balance on the vapor blanket. As mentioned above, the q heat flux can be calculated by heat conduction through the thickness of sublayer δ with the temperature difference of heated wall T w and the vapor blanket T sat . i.e.: In this equation k f is the conductivity of liquid. By newton's law of cooling, the heat flux also can be expressed as: From Eq. 1 and 2, the subcooled boiling heat transfer coefficient h can be obtained as Eq. 3: Thus to determine the heat transfer coefficient, it is needed to calculate the sublayer thickness δ which can be obtained by force balance on the vapor blanket in the axial and the radial directions (Fig. 2). For the two-phase flow, the effective fluid properties such as density ρ and viscosity µ need to be modified first.

Effective Homogeneous Fluid Properties
The homogeneous two-phase flow model is assumed to be suitable for the present subcooled or low-quality flow boiling conditions.

Magnitude of True Quality
The true quality x can be evaluated by using Saha and Zuber (1974) formula Eq. 4: where, x e is the thermodynamic equilibrium quality and x d is the thermodynamic quality at the point of bubble detachment from that heated wall. The value of x d is Eq. 5: where, c pf the specific heat capacity of liquid, D the inner diameter of the tube, H fg is the latent heat of vaporization, G the mass velocity and P ef the Peclet number of liquid. Note that, as x e ≤x d the true quality is identically zero for the single-phase flow.

ρ and µ for Homogeneous Flow
For the homogeneous two-phase flow of true quality x, the fluid density ρ is generally given by Eq. 6: In this equation ρ g and ρ f respectively, denote the density of vapor and liquid, while the mean two-phase viscosity µ will be evaluated by using the formula of Dukler et al. (1964) Eq. 7: where, µ g and µ f represent the viscosity of vapor and liquid; and v g and v f are specific volume of vapor and liquid, respectively.

Relative Velocity of Vapor Blanket to Liquid
At high heat flux condition, the vapor blanket is formed by the coalescence of small bubbles. The vapor blanket is assumed to be a distorted cylinder of length L b and D b diameter, which forms a flat interface near the wall. Consider the force balance in the axial direction, the relative velocity of the blanket with respect to the liquid will be determined by a balance between the buoyancy force F b and the drag force F d Eq. 8 and 9: Where: And: in which ∆ρ = ρ f -ρ g is the density difference between the two phases, C D the drag coefficient and U bL the relative velocity of the vapor blanket with respect to the liquid at the position corresponding to the centerline of the blanket. The negative sign in Eq. 10 indicates that the direction of the drag force is opposite to the flow direction. Chan and Prince (1965) proposed an expression of the drag coefficient for a small deformed bubble Eq. 11: Since the vapor blanket is formed initially by the coalescence a vertical column of small bubbles, the Science Publications JCS equivalent diameter D eb and D b are assumed to be equal to the bubble departure diameter (Lee and Mudawwar, 1988) and the latter can be obtained from the correlation of Cole and Rohsenow (1968) Eq. 12: where, σ is the surface tension. Weisman and Pei (1983) indicated that the bubbles are approximately ellipsoidal with the ratio of long to short ellipse axis being about 3:1 in the case of high heat flux. So the length of the vapor blanket L b is assumed to be three times of the bubble diameter D b . Combining Eq. 8 through (12) with the relations D eb =D b and L b =3D b gives the relative velocity of the bubble with respect to the liquid as Eq. 13:

Liquid Velocity Gradient
The local liquid velocity gradient can be evaluated by the force balance of a vapor blanket in the radial direction and the velocity profile distribution.

Force Balance in Radial Direction
The force balance of a vapor blanket in the radial direction is shown as Fig. 2. Vapro generation due to sublayer evaporation creates a rate of change of momentum F 1 which pushes the vapor blanket away from the wall. However, the lateral motion of the blanket is resisted by a lateral force F R caused by the vapor blanket rotation, which is resulted due to the velocity gradient associated with the liquid boundary layer in the tube. The inertial force F 1 is given by Eq. 14: where, V b is the vapor velocity due to evaporation of the sublayerand can be expressed as Eq. 15: Beyerlein et al. (1985) derived an expression for the lateral force on a bubble in turbulent two-phase flow in a vertical tube. The lateral force on the vapor blanket is determined by the relative velocity of the blanket and the gradient of the liquid velocity profile, i.e Eq. 16: where, U L is the local liquid velocity and C is a parameter which accounts for the effects of turbulent fluctuations and local bubble concentration on the rotation of the vapor blanket. From a liquid-gas two-phase flow experiment in the adiabatic boundary condition, Beyerlein et al. (1985) found that C is a function of the average void fraction and the liquid Reynolds number. In the present model for vapor-liquid system with wall heating, the dependences of two-phase Reynolds number and the boiling number in the parameter C are taken into account and is modeled by Eq. 17: In which a 1 , a 2 and a 3 are empirical constants, Re = GD/µ is the effective Reynolds number for homogeneous flow and Bo * is a modified boiling number defined as Eq. 18: In which Bo = q/(GH fg ) is the conventional boiling number and α is the void fraction Eq. 19: Upon combining the above equations, the liquid velocity gradient ∂U L /∂ y can be evaluated as the values of a 1 , a 2 and a 3 are treated known quantities, it is Eq. 20:

Velocity Profile Distribution
Applying the Karman's three-layer structure of turbulent flow in a tube, the friction velocity U + can be derived by the non-dimensional wall distance y + , that is Eq. 21a-d: The effective value of he velocity gradient causing vapor blanket rotation can be approximated as velocity gradient at the radial position of δ+ D b /2 from the wall. Lee and Mudawwar (1988) found that the liquid velocity profile around the vapor blanket locates in the buffer region and the effective velocity gradient can be calculated from Eq. (21a-d):

Thickness of Sublayer
Comparing Eq. 20 and 24, the thickness of sublayer δ is obtained by: Thus, the heat transfer coefficient h can be predicted by substituting Eq. 25 into Eq. 3.

Velocity of Vapor Blanket
To determine the onset of CHF, it is needed to calculate the velocity of vapor blanket U b which can be approximated as the superposition of the local liquid velocity at the radial position of δ+D b /2 from the wall and the relative vapor blanket velocity. Therefore, the velocity of the vapor blanket can be calculated from Eq. 21b and 13, it is Eq. 26: Figure 3 represents the flow boiling configuration of subcooled or low-quality flow immediately just before and occurrence of CHF. This phenomenon is called Helmholt instability which causes a wavy motion of the sublayer-vapor blanket interface. Based on the concept of Helmholtz instability as the relative velocity of the two streams, i.e., the liquid in the sublayer and the vapor in the vapor blanket, reaches a critical value, a wavy motion of the interface is induced. Due to the instability nature, the amplitude of the wavy motion can be amplified. The critical value will be evaluated based on the following assumptions:

Critical Velocity for Oneset of CHF
• The length of the blanket changes suddenly to the Helmholtz critical wavelength when the vapor blanket velocity reaches the critical valueand the sublayer also adjust its thickness • The vapor blanket touches the heated surface as result of the Helmholtz instability, the dry patch persists and spreads very quickly, which induces a sudden rise in wall temperature • CHF occurs when the rate of sublayer mass loss by evaporation exceeds the mass flow rate of the liquid entering the sublayer from the core region

Thickness of Sublayer
Since CHF is postulated to occur as a result of Helmholtz instability, the length of the sublayer and the vapor blanket are assumed to be equal to the Helmholtz critical wavelength L bH , as shown in Fig. 3a, i.e Eq. 27: In which U bH , is the critical velocity of the vapor blanket to onset of CHF, U m the liquid velocity in the sublayer. Since U bH is always much higher than U m , the above expression can be reduced to: Since the occurrence of CHF is a result of local sublayer dry out, the minimum critical heat flux q c necessary to evaporate the mass flux in the sublayer (see Fig. 3a is: where, δ H is the thickness of the sublayer to onset of CHF, T m the temperature of the liquid entering the sublayer, G m the liquid mass flux flowing into the sublayer, can be express as Eq. 30: Since H fg = C pf and in general, T sat = T m the second term in RHS of Eq. 29 can be neglected. Therefore, combining Eq. 28 through (30) gives δ H the thickness of the sublayer as Eq. 31:

Relative Velocity of Capor Blanket to Liquid
Near the CHF condition, the relative velocity of the vapor blanket to local liquid U bLH can be determined by a balance between the buoyancy force F B and the drag force F D Eq. 33 and 34: Where: And: In which C D is the drag coefficient, can be evaluated through the Chan and Prince (1965) Eq. 35: Combining Eq. 28 and 32 through 35 gives the relative velocity of the vapor blanket:

Critical Velocity of Vapor Blanket
The critical velocity of the vapor blanket is derived by the supersession of local liquid velocity at the radial position of δ H +D b /2 from the wall and the blanket relative velocity U bH which can be obtained from Eq. 21b and 36. it is Eq. 37: Thus, the critical blanket velocity U bH can be predicted by the numerical iteration.

Verification of Heat Transfer Model
By a regression analysis of 321 experimental data in the simulation of Light Water Reactors (LWRs) conditions (P = 6.9~15.5MPa), that is conducted by Lin et al. (1994), empirical constants a 1 , a 2 and a 3 are found to be 0.50, 1.20 and 2.21, respectively, Fig. 4 provides a comparison of predicted and experimental heat transfer coefficient for subcooled flow boiling at high heat flux (0.6 CHF < q<CHF). The predictions agree well with the experimental dataand the majority of data points fall within the error band of 40%.
Gungor and Winterton (1986) evaluated the various correlations for subcooled flow boiling, by comparing with measured data, they reported that the correlations by Moles and Shaw (1972); Shah (1977) and Gungor and Winterton (1986) are the best ones among others. Therefore, the above mentioned three correlations are selected in the study of subcooled flow boiling heat transfer coefficients. The comparison of these correlations with the experimental data is in Table 1. The correlations by Moles and Shaw (1972) and Gungor and Winterton (1986) that performed well with the high pressure (6.9MPa~15.5 MPa) boiling data give mean deviation of 32.7% (Moles and Shaw, 1972) and 46.1% (Gungor and Winterton, 1986), this error is closed to the report for subcooled flow biling data at the pressure from 13.2 to002020.3 MPa.
The correlation by Moles and Shaw present a better agreement with the experimental data within average deviation of 23.9% and mean deviation of 32.7%. However, this correlation is derived by direct dimensionless analysisand no significant physical meanings. Since the present model is developed based on the high heat flux mechanisms of vapor blanket and heat conduction in sublayer, it is superior to other correlations.

Verification of Steady-State CHF
To verify the trigger condition for onset of CHF, it is needed to examine the comparison between the critical value U bH and the vapor blanket velocity U b as CHF occurs.  Therefore, an analysis of relation between U bH and U b will be study based on a total of 362 data points of subcooled and low-quality saturated flow boiling water included in the tabular CHF data of the USSR Academy of Sciences (Collier, 1981), at the conditions of P = 6.9∼17.6MPa, G = 1000~5000 kg m −2 s and α = 0∼.07 . The analysis provided the Average Deviation (AD) of 5.3% and mean deviation of 22.8%. Thus, based on the Helmholtz instability at the sublayer-vapor blanket interface is the trigger condition for the onset of CHF as the blanket velocity reaches a critical value, this hypothesis is appropriate. Lee and Lin (1993) conducted an experiment flow transient CHF in the simulation of PWR at the linear mass flow decay rate from 0.1 to 30%/s. Figure 5 shows the variation of the inflow mass velocity and the wall temperature with the transient time at the flow decay rate of 0.1%/s. The onset of CHF is determined at the exit of the test section while the wall temperature excursion occurred, since the abrupt drop and the followed jump in wall temperature implies rewetting and dry out process of liquid film underneath the vapor blanket.   The experimental data will be compared with present model predictions. In present study, the predicted time to onset of CHF is determined as the blanket velocity is reached the critical velocity. Figure 6 show the variation a vapor blanket velocity and sublayer thickness with the transient time at the conditions of T in = 290°C, P = 15.5MPa, G initial = 3800 kg m −2 sec and flow decay rate =0.1%/. The CHF is occurred when U b is equal to U bH , it is worth to note that the thickness of sublayer is very near zero during the occurrence of CHF, this is a good evidence for sublayer dry out model. In Fig. 7, as the flow decay rate increase from 0.1-30%/s, it is shown that the experimental and the predicted transient time to CHF is decreased. Figure 8 provides a comparison of the predicted and measured time to ones of CHF. The prediction agrees well with the experimental dataand the majority of data points fall within the error bank of 30%

Verification of Tr0061nsient CHF
The well-known CHF correction of BandW-2 (Gellerstedy, 1969), EPR-1 (Fighetti and Reddy, 1983) and Lin et al. (1988a) are used to predict the onset of CHF in comparison with the present model. The comparison is shown in Table 2. Only the BandW-2 correlation with 13.9% of mean deviation is better than the prediction of present model. But the process from boiling heat transfer to the onset of CHF only can be demonstrated by the present model, this is why the present model is superior to the other correlations.

CONCLUSION
A theoretical model based on the Helmholtz instability has been developed for the evaluation of heat transfer performance and occurrence of CHF in lowquality flow boiling problem. For prediction of the heat transfer coefficient and the critical heat flux, the validity of the present model has been demonstrated by comparing with the existing experimental data and empirical corrections. However, it is worthy to note that the present model is developed based on the flow and heat transfer mechanisms. Therefore, unlike most of the previous models, the present one is of more physical meanings. It is believed that the present theoretical model is also useful to the analysis of the cooling technology related to the flow boiling heat transfer.