Simulation of Boiling Flow Experiments Close to CHF with the Neptune_CFD Code

Abstract

A three-dimensional two-fluid code Neptune_CFD has been validated against the Arizona State University (ASU) and DEBORA boiling flow experiments. Two-phase flow processes in the subcooled flow boiling regime have been studied on ASU experiments. Within this scope a new wall function has been implemented in the Neptune_CFD code aiming to improve the prediction of flow parameters in the near-wall region. The capability of the code to predict the boiling flow regime close to critical heat flux (CHF) conditions has been verified on selected DEBORA experiments. To predict the onset of CHF regime, a simplified model based on the near-wall values of gas volume fraction was used. The results have shown that the code is able to predict the wall temperature increase and the sharp void fraction peak near the heated wall, which are characteristic phenomena for CHF conditions.

1. Introduction

In the case of offnormal operation of the pressurized water reactor (increase of coolant inlet temperature at full power and full pressure), a nucleate flow boiling on the surface of the fuel rods may occur, which influences the operating and accident conditions in several ways. The increase of vapour content due to the boiling in narrow flow passages between the fuel rods increases the pressure drop through the reactor core and reduces the moderation ability of the coolant (effect on the core reactivity). When the applied heat flux reaches the CHF value, the departure from nucleate boiling (DNB) may be observed, resulting in rapid reduction of the heat removal ability. Overheating and consequently damaging of the fuel rods may occur. A better understanding of this phenomenon is a basis for optimization of fuel elements and for planning of appropriate safety measures in the case of fuel rods overheating.

In this study, the capability of the three-dimensional code Neptune_CFD to simulate the local boiling flow processes over a wide range of operating conditions, including those close to CHF, has been assessed. This work has been carried out within the scope of NURESIM project (NUclear REactor SIMulations, 6th EURATOM FP). The main objective of the NURESIM activity related to CHF is to use two-phase computational fluid dynamics (CFD) as a tool for understanding boiling flow processes. Numerical simulations of boiling flows up to CHF conditions should help new fuel assembly design and contribute to better local predictions of CHF phenomenon [1].

The three-dimensional two-fluid code Neptune_CFD [2] is developed within the framework of the Neptune project, financially supported by Commissariat à l'Énergie Atomique (CEA), Électricité de France (EDF), Institut de Radioprotection et de Sûreté Nucléaire (IRSN) and AREVA-NP. It has been more specifically designed for simulation of transients in nuclear power plants and it is currently used and further developed within the NURESIM project. To validate the code, relevant benchmark experiments have to be used, which provide useful local information about distributions of flow parameters (e.g., cross-sectional profiles of void fraction, phase velocities, temperature, bubbles size, etc.). Nucleate boiling processes in the subcooled boiling regime were studied on ASU experiments. Boiling flow close to CHF conditions was investigated on the selected experiments performed in DEBORA facility located at Commissariat à l'Énergie Atomique (CEA), Grenoble.

2. Experimental Data

The ASU and DEBORA boiling flow experiments were used to validate the boiling model of the Neptune_CFD code. Considered experimental facilities differ in the geometry of the test channel (annulus, pipe) and also in working fluids. In ASU experiments the refrigerant R-113 is used, while DEBORA facility uses the refrigerant R-12 as a working fluid. Density ratios and other scaling numbers are therefore different.

The measurement section of the ASU experimental facility [3] consists of a vertical annular channel with the heated inner tube with outer diameter of 15.8 mm and insulated outer tube with inner diameter of 38.02 mm. The total length of the annulus is 3.66 m and the 2.75 m long upper part of the inner tube is heated by the direct current. The 0.91 m long lower part of the annulus is not heated. The local measurements of transversal profiles of void fraction, phase velocities, velocity fluctuations, and liquid temperature were performed at a single axial location located 1.99 m downstream from the beginning of the heated section. The measurement probes and measurement techniques used in ASU experiments are described in the original paper of Roy et al. [3].

The DEBORA experiments (performed at CEA [4, 5]) were selected to analyze boiling processes close to CHF. The refrigerant R-12 at the pressure conditions of about 30 bar has been used as the working fluid to simulate steam-water flow at pressurized water reactor (PWR) conditions. The test section consists of a vertical pipe with internal diameter of 19.2 mm, divided into three axial parts: the adiabatic inlet section (1 m length), the heated section (3.5 m length), and the adiabatic outlet section (≈0.5 m length). At the end of the heated section, the radial profiles of void fraction, gas velocity, and bubble diameter were measured by the optical probe. ASU and DEBORA experimental conditions are presented in Table 1.

Table 1 

ASU [3] and DEBORA [4, 5] experimental conditions.

3. Physical Modelling in the Neptune_CFD Code

The basic model of Neptune_CFD is the classical six-equation two-fluid model together with 𝑘-𝜀 transport equations used for modelling of the liquid phase turbulence [2]. The version V.1.0.6 has been used to perform numerical simulations. Only some most relevant models and new model improvements used for the computations presented in this study are described.

3.1. Turbulence Modelling

The turbulent stress tensor for the liquid phase is modelled using the Boussinesq approximation𝜏^turb_𝑙=−𝜌_𝑙⟨^⇀𝑢^′_𝑙^⇀𝑢^′_𝑙⟩=𝜇^turb_𝑙(∇^⇀𝑢_𝑙+(∇^⇀𝑢_𝑙)^𝑇2)−3𝐼(𝜌_𝑙𝑘_𝑙+𝜇^turb_𝑙∇⋅^⇀𝑢_𝑙), (1) where ^⇀𝑢^′_𝑙 is the fluctuating part of liquid velocity, ^⇀𝑢_𝑙 is the liquid velocity in axial direction, 𝑘_𝑙 is the turbulent kinetic energy of liquid phase 𝑘_𝑙=0.5⟨^⇀𝑢^′_𝑙^⇀𝑢^′_𝑙⟩, 𝐼 is identity tensor, and 𝜇^turb_𝑙 is the eddy viscosity𝜇^turb_𝑙=𝐶_𝜇𝜌_𝑙𝑘²_𝑙𝜀_𝑙. (2) Eddy viscosity in (2) is defined by the turbulent kinetic energy 𝑘_𝑙 and its dissipation rate 𝜀_𝑙, both calculated from the two-equation 𝑘-𝜀 model. Parameter 𝐶_𝜇 is set to 0.09. The effect of wakes behind the bubbles on the liquid turbulence is taken into account by additional terms in 𝑘-𝜀 transport equations. These additional source terms represent the turbulent contribution of the gas phase on liquid and are modelled as follows [4]:𝑆^𝑘_𝑙=−(𝑀^⇀_𝐷+𝑀^⇀_𝐴𝑀)⋅(^⇀𝑢_𝑔−^⇀𝑢_𝑙𝑆),_𝜀𝑙=𝐶_𝜀3𝑆^𝑘_𝑙𝜏, (3) where 𝑀^⇀_𝐷 and 𝑀^⇀_𝐴𝑀 are interfacial drag and added mass volumetric forces and 𝜏 is a characteristic time for bubble induced turbulence, which depends on bubble departure diameter 𝑑_𝑏 and dissipation rate 𝜀_𝑙:𝑑𝜏=(^𝑏2𝜀_𝑙)^1/3. (4)

3.2. Interfacial Transfer Terms

The interfacial transfer of momentum is modelled by interfacial forces per unit volume, which include drag force 𝑀^⇀_𝐷, added mass force 𝑀^⇀_𝐴𝑀, lift force 𝑀^⇀_𝐿, and turbulent dispersion force 𝑀^⇀_TD. The interfacial drag force is calculated according to the Ishii correlation [6] and the added mass force is calculated by Zuber model [7]. Turbulent dispersion force, which considers the diffusion of the vapour phase due to the liquid phase turbulence, is calculated as^⇀𝐹_TD=−𝐶_TD𝜌_𝑙𝑘_𝑙∇𝛼, (5) where 𝐶_TD is the turbulent dispersion coefficient.

The bubble diameter in the bulk 𝑑_𝑏 determines the interfacial momentum transfer and interfacial heat and mass transfer. In the present work, constant values for bubble diameter have been used. The approximate values (1.2 mm for ASU and 0.3 mm for DEBORA experiments) were estimated from the available experimental data. The interfacial heat and mass transfer due to condensation in the subcooled bulk flow is modelled by Ranz-Marshall correlation [8]. Description of mass and heat transfer models is given in the code manual [2].

3.3. Wall to Fluid Transfer Model at Conditions Close to CHF

To take into account the phenomenon of temperature excursion at CHF conditions, the standard heat flux partitioning model of Kurul and Podowski [9] is extended by additional heat transfer to the gas phase. To distribute the wall heat flux between different phases, a phenomenological function 𝑓_𝛼1 is introduced. The function 𝑓_𝛼1 depends on liquid volume fraction 𝛼₁ and takes care of the numerically smooth transition between nucleate boiling regime and DNB regime. In the Neptune_CFD, the transition to DNB regime is determined by an arbitrary value for liquid volume fraction in the near-wall cell, set to 𝛼_1𝑐𝑟=0.2𝑓_𝛼1=1−0.5⋅exp[−20(𝛼₁−𝛼_1𝑐𝑟)]for𝛼₁≥𝛼_1𝑐𝑟,𝑓_𝛼1𝛼=0.5⋅(¹𝛼_1𝑐𝑟)^{20⋅𝛼1𝑐𝑟}for𝛼₁<𝛼_1𝑐𝑟. (6) The wall heat flux is then split into four different componentsΦ_𝑤=𝑓_𝛼1(Φ_𝐶1+Φ_𝑄+Φ_𝐸)+(1−𝑓_𝛼1)Φ_𝐶2, (7) where Φ_𝐶1 denotes the single-phase convection heat flux to the liquid, Φ_𝑄 denotes quenching heat flux that transfers cold liquid from the bulk flow to the wall periodically, Φ_𝐸 is the heat flux component needed to generate vapor bubbles, and Φ_𝐶2 denotes the heat flux used to preheat the vapor phase in the DNB regime. The first three heat flux components are extensively described in the code manual [2] and in our previous references [10, 11]. The convective heat flux to the vapor phase Φ_𝐶2 is modelled by a single-phase heat transfer at the wall:Φ_𝐶2=ℎ^𝑔_log(𝑇_𝑤−𝑇_𝑔), (8) where ℎ^𝑔_log is the convective heat transfer coefficient in the laminar regime [2], 𝑇_𝑤 is the wall temperature, and 𝑇_𝑔 is the temperature of the vapor phase. In the case of standard nucleate boiling model the function 𝑓_𝛼1 is equal to one, leading to zero value of Φ_𝐶2. The proposed approach does not have the ambition to predict the CHF triggering mechanism, but it is a promising attempt how to model the transition to DNB conditions. Local CHF triggering mechanisms are much more complex and depend on the microscale conditions, which cannot be resolved by CFD codes. Recently, Le Corre and Yao [12] have proposed a high resolution CHF model based on local wall hot spot mechanism, which has a potential to bridge the gap between local modelling of CHF mechanisms and implementation in the CFD codes.

3.4. Two-Phase Wall Law for Boiling Flows

At nucleate flow boiling, the liquid velocity profile in the boiling boundary layer is significantly disturbed by the bubble formation and detachment mechanisms on the heated wall. The use of single-phase log law for boiling flow calculations may lead to significant overprediction of liquid and gas velocities in the boundary region near the heated wall [10, 13, 14]. In our previous work, [10] a modified wall law following the formulation of Ramstorfer et al. [15] has been proposed. The main idea of the new wall function is that nucleating bubbles on the wall disturb the boundary layer flow in a similar way as the surface roughness. As a basis, a logarithmic law for turbulent flows over rough walls is used:𝑢⁺=1𝜅ln(𝑦⁺)+𝐵−Δ𝑢⁺, (9) where velocity 𝑢⁺=𝑢_𝑡/𝑢_𝑤 and distance from the wall 𝑦⁺=𝜌_𝑙𝑢_𝑤Δ𝑦/𝜇_𝑙 are written in nondimensional wall units scaled by wall friction velocity 𝑢_𝑤=𝜏_𝑤/𝜌_𝑙 (𝜏_𝑤 is the wall shear stress). Here 𝑢_𝑡 is the known velocity tangential to the wall and Δ𝑦 is the distance from the wall. Coefficients 𝜅 and 𝐵 are standard single-phase constants with the values of 0.41 and 5.3, respectively. The last term in (9) represents the offset of 𝑢⁺ due to the wall roughness:Δ𝑢⁺=1𝜅ln(1+𝐶_𝑘𝑟𝑘⁺_𝑟),𝑘⁺_𝑟>11.3,0,𝑘⁺_𝑟≤11.3, (10) where 𝐶_𝑘𝑟 is a roughness constant, which depends on the type of roughness (𝐶_𝑟=0.5 for sand-grain roughness) and 𝑘⁺_𝑟 is the roughness Reynolds number:𝑘⁺_𝑟=𝜌_𝑙𝑘_𝑟𝑢_𝜏𝜇_𝑙. (11) The quantity 𝑘_𝑟 represents the physical roughness height of the surface. For 𝑘⁺_𝑟>11.3, the wall is considered to be smooth, otherwise the wall is rough. Although Ramstorfer et al. [15] have studied the flow boiling in a horizontal channel this type of log-law may be applied to all boiling flows where the flow motion along the wall is dominant. The model assumes that the roughness height can be represented by a functional dependence on the bubble departure diameter 𝑑_𝑏𝑤 and by the contribution of nucleate boiling heat flux 𝑞_𝑛𝑏 to the total heat flux 𝑞_𝑤: 𝑘_𝑟=𝜂𝑑_𝑏𝑤(𝑞_𝑛𝑏𝑞_𝑤)^𝜍=𝜂𝑑_𝑏𝑤𝑞(1−^1𝜙𝑞_𝑤)^𝜍. (12)

In this study, bubble departure diameter is calculated according to extended Unal model [14]. The ratio of the nucleate boiling component to the total heat flux Φ_𝑛𝑏/Φ takes into account the thickening of the boiling boundary layer with increasing boiling activity. The coefficients 𝜂 and 𝜁 in (12) are empirical parameters set to the values 𝜂=0.5 and 𝜁=0.174 for the considered experimental cases. The proposed “boiling” law of the wall is implemented in the Neptune_CFD code in the form of blended linear-logarithmic wall function as follows:𝑢⁺=min(𝑦⁺,1𝜅ln(𝑦⁺)+𝐵−Δ𝑢⁺). (13) From this wall law the wall friction velocity is computed, which is used as a “near-wall” boundary condition for the liquid momentum equations.

4. Results

Experiments presented in Table 1 were simulated. The calculations for ASU and DEBORA experiments were performed on 2D numerical meshes (19×220 for ASU and 20×220 for DEBORA), where the lower number denotes the number of cells in radial direction. These meshes were selected as a reasonable compromise between the numerical accuracy and the computational effort [4, 10].

The heat transfer regime in ASU experiments is subcooled nucleate boiling, significantly below the CHF conditions. In this regime local nucleate boiling processes governing the lateral distribution of flow parameters were studied. In Figures 1 to 5 two different Neptune_CFD calculations are validated against ASU experiments. The calculations differ in wall function models. The “base” calculation uses standard single-phase law of the wall, whereas “wall func” calculation uses modified wall function model as described in Section 3.4. Other models, described in Section 3, are the same for both calculations.

(a)

(b)

(c)

Figure 1 

ASU experiments: void fraction profile.

Radial void fraction profiles are shown in Figure 1. The “wall func” simulation predicts somewhat wider boiling region comparing to the base calculation. The calculation of the axial gas velocity depends on the model for interfacial drag and interfacial area density (e.g., bubble size) whereas the axial liquid velocity profile in the wall boundary layer mainly depends on the wall friction, determined by the velocity wall function. Other influencing parameters are nondrag forces. The liquid and gas axial velocities are compared in Figures 2 and 3. The “base” calculation significantly overpredicts measured phase velocities for all experimental cases. In the calculation with the new wall function model, the liquid velocity adjacent to the wall is significantly decreased and is closer to the measured data, but overprediction somewhat away from the wall is still notable. Due to the coupling through the interfacial drag, a similar trend of decreased velocity near the heated wall may be observed also for the gas phase (Figure 3). In this case, the agreement between the “wall func” calculation and experiments is significantly improved over the entire gas velocity profile.

(a)

(b)

(c)

Figure 2 

ASU experiments: liquid velocity profile.

(a)

(b)

(c)

Figure 3 

ASU experiments: gas velocity profile.

Turbulent kinetic energy profiles are presented in Figure 4. The measurements showed that turbulent kinetic energy is the highest in the boiling region near the inner wall and then rapidly decreases towards centre of the channel. This trend was not adequately reproduced by the “base” calculation, whose 𝑘_𝑙 profile tends to be more gradual. The “base” calculation also significantly overpredicts the 𝑘_𝑙 values for ASU3 experiment. The “wall func” calculation gives improved prediction of the 𝑘_𝑙 profile for the experiment ASU3 over the entire channel cross-section but still overestimates measured data for the other two experiments. Due to the new wall function model, the calculated values of 𝑘_𝑙 are very high near the wall. In spite of discrepancies between the measured and simulated results, “wall func” calculations and measurements show similar trend-most of the turbulent kinetic energy production occurs in the boiling region close to the inner wall.

(a)

(b)

(c)

Figure 4 

ASU experiments: turbulent kinetic energy profile.

(a)

(b)

(c)

Figure 5 

ASU experiments: liquid temperature.

The liquid temperature profiles are compared in Figure 5. For both calculations, a reasonable agreement with experiments may be observed for all cases. The temperature profile is somewhat steeper in the case of “wall func” calculation.

The simulation results for DEBORA experiments are presented in Figures 6, 7, 8, and 9. The calculation denoted as “base” includes the models described in Section 3; the wall function is standard single-phase log law. In the calculation “𝐶_TD=2,” the turbulent dispersion coefficient 𝐶_TD in (4) is increased to the value of 2. Three different experimental cases (see Table 1) with significantly different operating conditions are simulated. According to the measured void fraction profile, the experiment DEB2 is supposed to be close to DNB, since a sharp increase of void fraction is observed near the heated wall. Radial void fraction profiles in Figure 6 are best predicted by the modified calculation with turbulent dispersion coefficient, though the profiles are smoother and underpredict the near-wall values. Especially in the case of DEB1 experiment, the “base” calculation strongly overpredicts measured void fraction on the heated wall, whereas the calculated void fraction in the centre of the pipe is too low. Both calculations are able to predict a void fraction jump near the wall at DEB2 experiment.

(a)

(b)

(c)

Figure 6 

DEBORA experiments: void fraction profile.

(a)

(b)

(c)

Figure 7 

DEBORA experiments: gas velocity profile.

(a)

(b)

(c)

Figure 8 

DEBORA experiments: liquid temperature profile.

Figure 9 

DEBORA experiments: distribution of heat flux components.

Gas velocity profiles are presented in Figure 7. In general, both calculations tend to overpredict the measured gas velocities in the near-wall region. A better agreement can be observed in for DEB1 experiment, which has a higher mass flow rate and heat flux than the other two experiments. The experiments DEB2 and DEB3 differ only in the inlet temperature. The simulations with increased 𝐶_TD give profiles closer to measured data. The measured liquid temperature profile is available only for DEB1 experiment. Like in the case of other flow parameters, the calculation with increased 𝐶_TD predicts smoother temperature profiles closer to experimental data.

Distribution of heat flux components for DEB2 case is presented in Figure 9. Heat fluxes are multiplied by function 𝑓_𝛼1 according to (6). The liquid single-phase convection component Φ_𝐶1 decreases with increasing evaporation heat flux Φ_𝐸. Towards the end of the heated channel the convective heat flux to the gas phase Φ_𝐶2 is activated, which denotes the onset of CHF conditions. According to (6), the incipience of CHF is defined by function 𝑓_𝛼1 and predefined near-wall void fraction value 0.8 (Figure 10). Although the wall temperature shows fluctuating behavior along the wall after the onset of CHF, the temperature excursion at CHF is fair predicted (Figure 10).

Figure 10 

DEBORA experiments: axial distribution of near-wall void fraction, 𝑓𝛼1, wall temperature, saturation temperature.

5. Conclusions

The set of ASU and DEBORA boiling flow experiments have been calculated with the version V1.0.6 of the Neptune_CFD code. The distribution of flow parameters in the subcooled boiling regime was studied on ASU experiments. Within this scope, a new wall function model based on the surface roughness analogy has been implemented, leading to improved agreement of flow parameters with measured data.

The capability of the code to predict boiling flow phenomena close to CHF conditions was tested on the selected DEBORA experiment. Though the simplified model for the onset of CHF was used, it was demonstrated that the code is able to predict the wall temperature increase and the sharp void fraction peak characteristic for CHF conditions. However, it should be emphasized that more generic CHF criteria based on physical mechanisms need to be developed. Further investigations of critical heat flux mechanisms are therefore necessary, both experimentally and numerically.

Acknowledgments

This work was financially supported in part by the Slovenian Research Agency and partly through the NURESIM project (6th Euratom Framework Program).