Two-dimensional axisymmetric simulations of pressurized thermal shock (PTS) phenomena through Neptune_CFD module are presented aiming at two-phase models validation against experimental data. Because of PTS complexity, only some thermal-hydraulic aspects were considered. Two different flow configurations were studied, occurring when emergency core cooling (ECC) water is injected in an uncovered cold leg of a pressurized water reactor (PWR)—a plunging water jet entering a free surface, and a stratified steam-water flow. Some standard and new implemented models were tested: modified turbulent - models with turbulence production induced by interfacial friction, models for the drag coefficient, and interfacial heat transfer models. Quite good agreement with experimental data was achieved with best performing models for both test cases, even if a further improvement in phase change modelling would be suitable for nuclear technology applications.

1. Introduction

The Integrated Project, European Platform for Nuclear Reactor Simulations (NURESIM), financially supported by the European Commission, aims at developing a common European standard software platform for modelling, recording, and recovering computer simulation data for current and future nuclear reactor system [13]. Neptune_CFD [46] is the thermal-hydraulic two-phase CFD tool of NURESIM and is designed to simulate most of two-phase flow configurations encountered in nuclear reactor power plants. Neptune_CFD is developed within the framework of the NEPTUNE project, financially supported by Commissariat à Énergie Atomique (CEA), Électricité de France (EDF), Institut de Radioprotection et de Sûreté Nucléaire (IRSN), and AREVA-NP.

One task of the NURESIM Project [7] is the prediction of PTS phenomena through computational fluid dynamics’ (CFDs) codes, in order to improve the operational safety and remnant life assessment of the PWRs.

A PTS scenario limiting to the reactor pressure vessel lifetime is the cold water ECC injection into the cold leg during a small-break loss of coolant accident (SBLOCA) [8]; rapid wall cooling can lead to strong thermal gradients and consequently to high stresses in the pressurized components, while local reduction of fracture toughness occurs due to temperature decrease. Complex phenomena take place when cold water is transported from injection nozzle to the downcomer, such as

(i)turbulent mixing of momentum and heat in the jet region and downstream of the impingement zone,(ii)stratified two-phase flow with condensation at the free surface. The two flow configurations considered in this paper are likely to share common physical features with these scenarios and represent challenging cases for multiphase models validation. As a matter of fact, they were identified as relevant PTS scenarios [9] and selected as test cases for CFD codes validation [10] within the ECORA Project.

The first concerns a jet flow impinging on a free surface, with air carry under and subsequent bubble dispersion in the water bath. Since turbulence strongly influences the condensation and the diffusion of heat within the water, it also influences the temperature field and plays a major role in the severity of the PTS Scenario. CFD models have to capture this turbulence production by jet kinetic energy in order to provide reliable predictions, and the Iguchi tests [11] are used as separate effect tests for code validation.

The second one relates to a turbulent stratified steam-water flow with interfacial heat and mass transfers. Simulating such a problem involves many critical aspects for the successful validation of CFD models: turbulence should be accurately predicted near the interface, and turbulence models should account for anisotropy effects; all interphase transport source terms have to be accurately represented in the solved equations, needing a numerically stable and robust solver. Numerical results are compared with experimental data from the cocurrent LAOKOON test case at high Reynolds number of steam [1214].

2. Modelling

In this paper, a local 3D two-fluid approach [15] for turbulent flows with/without condensation is presented. In this approach, a set of local balance equations of mass, momentum, and energy is written for each phase. These balance equations are obtained by time averaging (or ensemble averaging) the local instantaneous balance equations written for the two phases. When the averaging operation is performed, most of the information about the interfacial configuration and exchanges is lost. As a consequence, a certain number of constitutive relations are needed for the closure of the equations system.

Three different types of closure relations can be identified: those which express the intraphase exchanges (molecular and turbulent transfer terms), those which express the interphase exchanges (interfacial transfer terms), and those which express the interaction between each phase and the walls (wall transfer terms) [16].

Together with Neptune_CFD standard models, various modified models developed at CEA/Grenoble are tested regarding the turbulence production induced by interfacial friction, the drag coefficient, and the interfacial heat transfer [17]. In the following will be presented the balance equations together with the closure laws of the most important terms used to simulate the considered two-phase problems.

2.1. Mass, Momentum, and Energy Averaged Balance Equations

The Neptune_CFD code is based on the classical two-fluid model, which consists of the following six balance equations.

(i) Two mass-balance equations are where the quantities and are the averaged fraction of presence, the averaged density, and the averaged centre of mass velocity for phase , with the phase index k being equal to L for the liquid phase and to G for the gaseous phase. The right-hand side (RHS) of (1), denoted by , is the rate of phase change (evaporation or condensation) per unit volume of mixture.

(ii) Two momentum balance equations are where p is the averaged pressure, common to the two phases, and g is the gravity acceleration vector. The two tensors are the averaged molecular and turbulent Reynolds stress tensors, respectively, and the vector is the averaged interfacial momentum transfer between phases.

(iii) Two total enthalpy balance equations are where is the bulk-averaged enthalpy of phase k and is the interfacial-weighted averaged enthalpy. The two vectors are the molecular and turbulent heat fluxes. The term is the heat flux exchanged between phase k and the interface per unit volume, where is the interfacial area concentration (or interfacial area per unit volume), and the term is a possible heat exchange term between phase k and the wall.

2.2. Turbulent Transfer Terms

In the Neptune_CFD code [5], for each phase k, the Reynolds stress tensor is closed using a Boussinesq-like hypothesis [18].

A two-equation model for the calculation of the turbulent eddy viscosity is used, which is an extension to multiphase flows of the classical model used in single phase flows. The two equations for the turbulent kinetic energy and the turbulent dissipation rate are written in the following nonconservative form: where the turbulent viscosity is given by with .

represents the (positive) production due to the mean velocity gradients, and is a stratification attenuation term modelling the correlation between fluctuating densities and velocity (more details are available in [5]). and are the two classical constants taken from the single-phase model.

By the models being primarily valid only for turbulent core flows, the near-wall region will be modelled with the standard wall function approach [18].

The two terms and take into account the additional turbulence production (or destruction) due to the influence of each phase on the other one.

2.2.1. Various Options for Bubbly Flows Tested in the Jet Experiment Calculations

“Ke Liq + TRC”
TRC stands for turbulence reverse coupling. For the liquid phase of a bubbly flow, the two terms and are given by where and are, respectively, the averaged drag and added mass forces, is a constant parameter equal to 0.6 in our calculations, and is the characteristic time defined as a function of the imposed bubble diameter and the turbulence dissipation rate . It represents the time scale of liquid eddies having the same size as the bubble diameter .

“Ke EDF”
The difference to notice with the previous one is that the dissipation production term is

“Ke Liq”
Means that the influence of these interfacial turbulence effects of (6a), (6b) was not considered in the calculations.

Means in the jet experiment calculations that the transport equations were not used for the gas phase.

2.3. Interfacial Transfer Terms
2.3.1. Heat and Mass Transfers

If the mechanical terms are neglected in comparison to the thermal terms in the averaged form of the energy jump condition, it reduces to [15] This important relation (together with the mass jump condition ) allows to compute the mass-transfer term as a function of the two volumetric heat fluxes and the interfacial-averaged enthalpies . We assume that the interfacial-averaged enthalpies are identical to the phase-averaged enthalpies . Each interfacial heat transfer term is the product of the interfacial heat flux density which is expressed as where is a heat transfer coefficient, and are, respectively, the averaged temperature of phase k and the saturation temperature, and α is the void fraction (or vapour void fraction).

The interfacial area concentration in bubbly flow is In stratified flow, with the “continuous approach,” it is In stratified flow, with the “discrete approach,” it is where is the cell size equal to the length of the segment which includes the gravity centre of the cell and which has the direction normal to the free surface.

A two-phase liquid-vapour flow is considered, where the liquid is identified by index “” and the vapour by index “”. model is not relevant for tests presented here: air-water or saturated vapour. is calculated with Two models for the heat exchange coefficient , as follows, are used in the present calculations. They are selected amongst the different choices available in the standard version of the code because they are dedicated for PTS applications. Both models can be applied either with a discrete approach, where the heat transfer is calculated only in interface cells and is zero elsewhere, using (9c), or with a continuous approach, where the heat transfer is calculated in all the grid cells of the domain, using (9b) at each time step. In practice, the liquid volume fraction gradient which gives the interfacial area is nonzero only near the interface, and consequently the calculated heat transfer tends rapidly to zero elsewhere.

Neptune_CFD 2004 [19] (HD1)
with Since this model takes into account only condensation effects at the water-steam interface, it has been completed by a residual droplet contribution (in the upper zone where ) taken as a “return to saturation” term, with a constant time scale arbitrary equal to 1 second, and the weighting function .

Coste ICMF’2004 (HD2)
This model was implemented by CEA/Grenoble. Like many others, it is based on the old concept of surface renewal (Higbie, 1935 [20]). It differs in the definition of the characteristic renewal frequency scale. As discussed in [17, 21], the assumption of renewal by Kolmogorov eddies gives rise to a theoretical contradiction when the Prandtl number approaches unity, that is generally true for water. An alternative was then proposed [21], called HD2 hereafter, where the frequency is calculated with the Kolmogorov eddies length scale and velocity fluctuation due to turbulence. The validity domain of the surface renewal model framework is then The heat transfer coefficient is where represents the large eddies length. This model has been validated with Simmer and Neptune_CFD codes on about twenty test cases of COSI experiment. From this point of view, its validity domain is given by the characteristic nondimensional numbers of this experiment.

2.3.2. Momentum Transfer

In the case of a bubbly flow, the interfacial transfer of momentum appearing in (2) is assumed to be the sum of five forces: where we have assumed that the interfacial-averaged velocity is equal to the phase-averaged velocity (), the first term on the right-hand side of (14) represents the interfacial transfer of momentum associated to the interfacial transfer of mass. The other terms are, respectively, the averaged drag, added mass, lift and turbulent dispersion forces per unit volume.

Standard Models

Drag Force is
where is the nondimensional drag coefficient.
Concerning the drag coefficient closure law, the separated phases’ model and the Ishii correlation were considered. The separated phases’ model, used for liquid-gas separated flows, is a Simmer-like model [22] which considers either dispersed gas bubbles in a continuous liquid flow, or dispersed liquid droplets in a continuous gas flow with regard to the volumetric fraction.
The Ishii’s empirical correlation, used in case of bubbly flow, provides the automatic calculation of the drag coefficient based on the local regime: where is the surface tension. Equation (16) assumes that we are in the distorted bubble regime.
Added Mass Force is where is the added mass coefficient which is equal to 0.5 in the case of spherical bubbles.
Lift Force is where is the lift coefficient. This coefficient is equal to 0.5 in the particular case of a weakly rotational flow around a spherical bubble in the limit of infinite Reynolds number.
Turbulent Dispersion Force is where is a numerical constant of order 1.

New Implemented Models

Drag Force
This model too was inspired by the Simmer code but differs from the separated phases in the definition of the drag coefficient for (mixing case, defined previously). In fact, the separated phases’ model in this case considers that in every cell is present different percentage of liquid with bubbles inclusion and gas with drops inclusion, depending on the value of . In free surface flows with a flat surface, bubbles and droplets are not present, and such a model will be not coherent with physical reality. In fact, this model will lead to an overestimation of the friction due to high bubble drag coefficient, which depends on the fluid density. In order to account roughly for this, as a first step towards the adaptation of drag force closure law to the free surface case, Dev model [21, 23] multiplies the bubble drag coefficient for bubbly flows by a factor of . A second step is to use a wall law type approach but it has been implemented too recently [24] to be tested in the present work.

3. Calculations Discussion and Results

Neptune_CFD is based on a fully unstructured finite volume meshing, together with a collocated arrangement for all flow variables [16]. The solver, based on an elliptic oriented fractional step approach, is able to simulate multifield and multiphase flows. The nonlinear behaviour between pressure and the phase fractions and the symmetric treatment of the fields are taken into account in an iterative procedure, within the time step.

3.1. Prediction of Turbulence Distribution below a Plunging Jet

Experimental data of Iguchi et al. [11] for a plunging water jet entering a free surface (see Figure 1) are considered to evaluate prediction capabilities of the two-phase models implemented in the code.

The considered experimental conditions and measurements are summarised in Table 1.

From the experiments, the authors found that in this configuration the jet produces a significant amount of small bubbles (m or 1 mm) in the entire water bath.

In all calculations presented in this section, the fluid domain was represented by nonuniform 2D spatial grids, adequately refined in the near wall regions. All boundary conditions were taken as suggested by the NURESIM Programme specifications [8].

Preliminary calculations testing different grids showed that the predicted flow reaches a steady-state configuration in a physical time of 4 seconds, and resulted in the selection of a reference spatial meshing (with  cells in the vessel and in the pipe).

First results showed that the numerical simulation effectively accounts for air entrainment near the bath surface (cm) and bubble dispersion in the whole water vessel ( and 30 cm), but with very small void fraction values .

Various two-phase turbulence models described in Section 2.2 were tested together with the separated phases model for the drag coefficient also evaluating the influence of interfacial turbulence effects. Predicted results [25] for the axial mean velocity and the root-mean-square value of axial fluctuation on the centreline along the vertical direction were generally in quite good agreement with experimental data (available for ), with some underestimation (see Figure 2).

Both the models considered for water (the standard and modified versions) underestimated the turbulence production and predicted almost the same radial profiles for the turbulence components and , failing to catch the anisotropy of the problem (see Figure 3).

This is a classical feature of the model; important differences are observed especially near the bath surface for the prediction. Nevertheless, moving downstream , the flow decelerates and the calculated values better match the experimental data.

Considering the “TRC” contribution (see (6a)) in the “ke liq” model brings worst results; contrary to what expected, it generally caused a further reduction in turbulence estimation, especially far from the surface (see Figure 4). Looking at (4), it can be observed that the additional turbulence production terms (the “TRC” contribution), and , act on both turbulence kinetic energy and turbulence dissipation equation, so that the overall effect can be either an increase or a decrease of turbulence. Trying to modify the relative importance of production and dissipation terms in the turbulent kinetic energy equation (4), by means of changing the values of bubble diameter, , and the parameter , did not lead to better results. Probably, a more accurate study is needed to clarify the relation between these parameters and TRC term effects.

A sensitivity analysis to mesh refinement was conducted for the best performing turbulence model, considering three successively refined grids (CASE1,  cells; CASE2,  cells; CASE3,  cells). At first, the separated phases’ model was adopted for the drag coefficient (see Figure 5(a)); then the Ishii’s correlation given in (16) was tested (see Figure 5(b)), and finally the added mass (AM), lift , and turbulent dispersion (TD) force contributions (see Figure 5(c)) described in Section 2.3.2 above were taken into account. The dependence of numerical results from grid refinement is strongly reduced but the convergence on spatial meshing is not reached. Moreover, the void fraction prediction always seems to change without coherence considering the different grids; this fact probably depends on the air entrainment modelling at the free surface.

Iguchi [11] presented a comparison between experimental measurements for the two-phase case and theoretical values for the single-phase free jet, showing that the mean velocity and turbulence characteristics were not affected by the bubbles, and agreed well with those for the single-phase case.

An analogous single-phase study was carried out finding that turbulence prediction was not considerably influenced by bubble entrainment, according to experimental evidence (see Figure 6).

Even if there is some confidence in the predicting capabilities of models on water jet impact effect, some uncertainties are left regarding the prediction of turbulent parameters near the free surface (important to determine the interfacial transfers) because experimental data are not available. Moreover, the effects of water jet on gas entrainment need to be clarified.

3.2. Prediction of Direct Contact Condensation in Turbulent Steam-Water Stratified Flow

This test case concerns a horizontal stratified flow of subcooled water and saturated dry steam along a rectangular straight channel with adiabatic walls. Available experimental data have been measured in the Technical University of Munich using the LAOKOON test facility [14]. The correspondent two-dimensional geometry recommended in [9] for CFD simulations is presented in Figure 7.

One experimental test was simulated, namely, the cocurrent case at high Reynolds number of steam. The regime parameters and the water temperature profile at the measurement section are listed in Table 2. The fluid domain was represented by uniform 2D grids, with different refinement. All boundary conditions were taken as suggested by the NURESIM Program specifications in [8].

Calculations were generally run for 40 seconds, since steady-state conditions were typically reached after 30 to 40 seconds. In all CFD simulations, small surface waves were observed, causing values of flow characteristics calculated near the interface to oscillate, so that presented results are time averaged between seconds and seconds.

Simulations were run testing at first standard two-phase models (see Figure 8(a)) and then the new implemented models (see Figure 8(b)).

Results showed that modified turbulence model and drag coefficient Dev model perform better than the standard versions, so they were chosen for the sensitivity study of the interface transfers models.

The two models HD1 and HD2 were compared. Predicted temperature profiles were qualitatively correct, and calculated condensation rates were all around 40%, according to experimental data (the dimensionless condensation rate (CR) is given by the ratio between the overall calculated mass transfer [kg/s] and the vapour mass-flow rate entering at inlet [kg/s]).

The HD2 (discrete approach) model resulted as the best performing model and was chosen as reference to develop a mesh sensitivity study (see Figure 9) (the discreet approach represents the free surface as a sharp discontinuity and allows to consider condensation only in surface cells; in the continuous approach, the free surface is smeared for few grid cells, and the flow variables are distributed continuously with high gradients at the phase separation surfaces).

As a result, temperature profiles seemed not very sensitive to mesh refinement in the bottom part of the channel, while some differences are observed near the interface. Concerning the condensation rate, calculated values were strongly influenced by mesh refinement: increasing axial refinement caused a rise in condensation , while increasing mesh refinement in height direction caused a reduction ; for mesh refinement in both directions the sensitivity was reduced .

The configuration is subject to the Kelvin-Helmotz instability. In physical reality, the surface tension is opposed to it at length scales which can be compared with the typical cells size of our meshes. In the calculations, no surface tension is taken into account on interfaces larger than cells size. Then, the instabilities tend to be overestimated with small cells. This is what happens in present calculations; when refining axially the mesh, the sensitivity on the condensation rate is clear. Such instabilities generate a higher microscale turbulence which in turn generates a higher condensation rate, so it is not possible to quantify how much can be attributed to one or the other. With large enough cells, the instabilities growth cannot be simulated, and therefore the calculated large interfaces are more stable. When refining radially, as the axial length is kept the same and remains quite large, the effect on the calculation is much less coming from an increase of the instabilities than from the modelling sensitivity to the mesh. This sensitivity rather due to the models is not satisfactory, and future work will be dedicated to it.

Moreover, since no exhaustive information is available about the water inlet turbulence, numerical experiments were carried out to investigate the results sensibility to turbulence boundary conditions.

At first only the water turbulence intensity was changed (see Figure 10(a)). Then, the equation which governs the relationship between the turbulence kinetic energy k and the turbulent dissipation rate was also modified (see Figure 10(b)). In both cases, water temperature profiles and condensation rate showed a strong dependence on turbulence boundary conditions, so that it seems necessary to carry on a close examination on this point.

4. Conclusions

This study presents a validation of Neptune_CFD against plunging jet data of Iguchi [11] and DCC on stratified steam-water flow for the LAOKOON experiment [14], with an extensive performance comparison of different two-phase models in predicting flow characteristics in PTS scenarios.

One test of Iguchi’s experiment is simulated, where small bubbles are entrained below the free surface. Numerical simulations effectively account for air entrainment and bubble dispersion but with very small void fraction values (less than 0.5%), which is qualitatively consistent with experimental observations (quantitative comparison is not possible since void fraction measurement is not available). Predictions of the mean velocity field were always in a rather good agreement with experimental data. Calculated turbulence was generally not bad but significant underestimation is obtained far from the jet axis region. Big differences are also observed in the prediction of turbulent velocities components near the free surface. In fact, none of the used turbulence models succeeded in predicting the anisotropy of the problem; however the main objective of the turbulence model is here to predict turbulent diffusion of heat in the liquid layer and a model may be sufficient. Considering the interfacial turbulence effects, in contrast to what expected, brought a further reduction in turbulence estimation. Normally, due to the small size of bubbles, there should not be a significant influence of bubbles on the turbulence as mentioned by Iguchi. It is then recommended to further consider the formulation of the interfacial production and dissipation terms in case of small bubbles. Anyway, an under prediction of turbulence would, in real PTS applications, lead to an overestimation of the scenario severity, since steam condensation in the bulk liquid is strongly enhanced by the turbulence level.

For the considered case of LAOKOON experiment, calculations were run at first with standard models for both drag coefficient and interface transfers. Recent models for free surfaces developed at CEA/Grenoble were also tested—a method for interfacial friction and two methods for interfacial heat transfers. In all CFD simulations, small surface waves were observed, causing flow characteristics calculated near the free surface to oscillate. Calculated condensation rates were often higher than measured value, probably because of a turbulence overestimation at the interface. Most of numerical predictions gave correct qualitative water temperature profiles at probe location, while some performance dissimilarities were found in the near surface region, and temperatures are underestimated in the bottom part of the channel. Considering recent free surface models allowed calculated values to better match experimental data, but as well as for standard models, numerical predictions were found to be mesh dependent.

Improving considered models, reducing grid sensitivity, and increasing accuracy in calculating turbulence and interface transfers would then be advantageous.

As a conclusion, we can state that the present work contributed to the assessment of CFD code applicability to PTS scenarios; presented simulations of Iguchi’s jet test demonstrated that models could predict reasonably well the jet induced turbulence and that it is important to consider also the coupling of the turbulence fields. LAOKOON study showed that the two-fluid approach is appropriate to study a stratified steam-liquid flow with condensation, even if further improvement on heat transfer modelling is required.


Understanding the code and clearing up the encountered problems would not be possible without the precious collaboration of CEA and University of Pisa colleagues. Special thanks are also due to the NEPTUNE Project team for their continuous support.