Simulation of Boiling Flow Experiments Close to CHF with the Neptune CFD Code

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.


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.

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.

PHYSICAL MODELLING IN THE NEPTUNE CFD CODE
The basic model of Neptune CFD is the classical sixequation two-fluid model together with k-ε 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.

Turbulence modelling
The turbulent stress tensor for the liquid phase is modelled using the Boussinesq approximation where u l is the fluctuating part of liquid velocity, u l is the liquid velocity in axial direction, k l is the turbulent kinetic energy of liquid phase k l = 0.5 u l u l , I is identity tensor, and μ turb l is the eddy viscosity Eddy viscosity in (2) is defined by the turbulent kinetic energy k l and its dissipation rate ε l , both calculated from the two-equation k-ε model.Parameter C μ is set to 0.09.The effect of wakes behind the bubbles on the liquid turbulence is taken into account by additional terms in k-ε transport equations.These additional source terms represent the turbulent contribution of the gas phase on liquid and are modelled as follows [4]: where MD and MAM are interfacial drag and added mass volumetric forces and τ is a characteristic time for bubble induced turbulence, which depends on bubble departure diameter d b and dissipation rate ε l :

Interfacial transfer terms
The interfacial transfer of momentum is modelled by interfacial forces per unit volume, which include drag force MD, added mass force MAM, lift force ML, and turbulent dispersion force MTD.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 where C TD is the turbulent dispersion coefficient.The bubble diameter in the bulk d b 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].

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 f α1 is introduced.The function f α1 depends on liquid volume fraction α 1 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 α 1cr = 0.2 The wall heat flux is then split into four different components where Φ C1 denotes the single-phase convection heat flux to the liquid, Φ Q denotes quenching heat flux that transfers cold liquid from the bulk flow to the wall periodically, Φ E is the heat flux component needed to generate vapor bubbles, and Φ C2 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 Φ C2 is modelled by a single-phase heat transfer at the wall: where h g log is the convective heat transfer coefficient in the laminar regime [2], T w is the wall temperature, and T g is the temperature of the vapor phase.In the case of standard nucleate boiling model the function f α1 is equal to one, leading to zero value of Φ C2 .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.

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: where velocity u + = u t /u w and distance from the wall y + = ρ l u w Δy/μ l are written in nondimensional wall units scaled by wall friction velocity u w = τ w /ρ l (τ w is the wall shear stress).Here u t is the known velocity tangential to the wall and Δy is the distance from the wall.Coefficients κ and B are standard single-phase constants with the values of 0.41 and 5.3, respectively.The last term in (9) represents the offset of u + due to the wall roughness: where C kr is a roughness constant, which depends on the type of roughness (C r = 0.5 for sand-grain roughness) and k + r is the roughness Reynolds number: The quantity k r represents the physical roughness height of the surface.For k + r > 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 d bw and by the contribution of nucleate boiling heat flux q nb to the total heat flux q w : 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 Φ nb /Φ 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 linearlogarithmic wall function as follows: From this wall law the wall friction velocity is computed, which is used as a "near-wall" boundary condition for the liquid momentum equations.

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.
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.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 k l profile tends to be more gradual.The "base" calculation also significantly overpredicts the k l values for ASU3 experiment.The "wall func" calculation gives improved prediction of the k l 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 k l 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.
The liquid temperature profiles are compared in Figure 5.For both calculations, a reasonable agreement with exper-iments 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 "C TD = 2," the turbulent dispersion coefficient C 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.
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 C 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 C 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 f α1 according to (6).The liquid single-phase convection component Φ C1 decreases with increasing evaporation heat flux Φ E .Towards the end of the heated channel the convective heat flux to the gas phase Φ C2 is activated, which denotes the onset of CHF conditions.According to (6), the incipience of CHF is defined by function f α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).

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