Two-Phase Flow Simulations for PTS Investigation by Means of Neptune_CFD Code

Abstract

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 [1–3]. Neptune_CFD [4–6] 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 à l^′É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 [12–14].

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 𝜕𝛼_𝑘𝜌_𝑘𝜕𝑡+∇⋅(𝛼_𝑘𝜌_𝑘𝑉_𝑘)=𝛤_𝑘,𝑘=𝐿,𝐺, (1) 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 𝜕𝛼_𝑘𝜌_𝑘𝑉_𝑘𝜕𝑡+∇⋅(𝛼_𝑘𝜌_𝑘𝑉_𝑘𝑉_𝑘)=−𝛼_𝑘∇𝑝+𝑀_𝑘+𝛼_𝑘𝜌_𝑘𝑔+∇⋅[𝛼_𝑘(𝜏_𝑘+𝜏^𝑇_𝑘)],𝑘=𝐿,𝐺, (2) where p is the averaged pressure, common to the two phases, and g is the gravity acceleration vector. The two tensors 𝜏_𝑘and𝜏^𝑇_𝑘 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 𝜕𝜕𝑡[𝛼_𝑘𝜌_𝑘(ℎ_𝑘+𝑉²_𝑘2)]+∇⋅(𝛼_𝑘𝜌_𝑘(ℎ_𝑘+𝑉²_𝑘2)𝑉_𝑘)=𝛼_𝑘𝜕𝑝𝜕𝑡+𝛼_𝑘𝜌_𝑘𝑔⋅𝑉_𝑘+𝛤_𝑘(ℎ_𝑘𝑖+𝑉²_𝑘2)+𝑞^′′_𝑘𝑖𝑎_𝑖+𝑞^′′′_𝑤𝑘−∇⋅[𝛼_𝑘(𝑞_𝑘+𝑞^𝑇_𝑘)],𝑘=𝐿,𝐺, (3) where ℎ_𝑘 is the bulk-averaged enthalpy of phase k and ℎ_𝑘𝑖 is the interfacial-weighted averaged enthalpy. The two vectors 𝑞_𝑘and𝑞^𝑇_𝑘 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: 𝜌_𝑘[𝜕𝐾_𝑘𝜕𝑡+𝑉_𝑘,𝑖𝜕𝐾_𝑘𝜕𝑥_𝑖]=1𝛼_𝑘𝜕𝜕𝑥_𝑗[𝛼_𝑘𝜇^𝑇_𝑘𝜎_𝐾𝜕𝐾_𝑘𝜕𝑥_𝑗]+𝜌_𝑘(PROD_𝑘+𝐺_𝑘−𝜀_𝑘)+𝑃^𝑖_𝐾,𝜌_𝑘[𝜕𝜀_𝑘𝜕𝑡+𝑉_𝑘,𝑖𝜕𝜀_𝑘𝜕𝑥_𝑖]=1𝛼_𝑘𝜕𝜕𝑥_𝑗[𝛼_𝑘𝜇^𝑇_𝑘𝜎_𝜀𝜕𝜀_𝑘𝜕𝑥_𝑗]+𝜌_𝑘𝜀_𝑘𝐾_𝑘(𝐶_𝜀1PROD_𝑘+𝐶_𝜀1max(𝐺_𝑘,0)−𝐶_𝜀2𝜀_𝑘)𝑃^𝑖_𝜀, (4)(5) where the turbulent viscosity is given by 𝜇^𝑇_𝑘=𝐶_𝜇𝜌_𝑘(𝐾²_𝑘/𝜀_𝑘) with 𝐶_𝜇=0.09.

PROD_𝑘 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]). 𝐶_𝜀1=1.44 and 𝐶_𝜀2=1.92 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 𝑃^𝑖_𝐾=−(𝑀^𝐷_𝐺+𝑀^AM_𝐺)⋅(𝑉_𝐺−𝑉_𝐿𝑃),_𝑖𝜀=𝐶_𝜀3𝑃^𝑖_𝐾𝜏𝑑,𝜏=(²³²𝜀_𝐿)^1/3, (6a) where 𝑀^𝐷_𝐺 and 𝑀^AM_𝐺 are, respectively, the averaged drag and added mass forces, 𝐶_𝜀3 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 𝑃^𝑖_𝜀=𝐶_𝜀4𝜀_𝑘𝐾_𝑘𝑃^𝑖_𝐾. (6b)

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

“Laminar”
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] _𝑘(𝛤_𝑘ℎ_𝑘𝑖+𝑞^′′_𝑘𝑖𝑎_𝑖)≈0. (7) 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 𝑞^′′_𝑘𝑖=𝐶_𝑘𝑖(𝑇_sat(𝑝)−𝑇_𝑘), (8) where 𝐶_𝑘𝑖 is a heat transfer coefficient, 𝑇_𝑘 and 𝑇_sat(𝑝) 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 𝑎_𝑖=6𝛼_𝐺𝑑₃₂. (9a) In stratified flow, with the “continuous approach,” it is 𝑎_𝑖=|∇𝛼_𝐺|. (9b) In stratified flow, with the “discrete approach,” it is 𝑎_𝑖=1ℎ^→_𝑛, (9c) 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 𝑞^′′_𝐿𝑖=𝐶_𝐿𝑖(𝑇_sat−𝑇_𝐿). (10) 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)
𝐶_𝐿𝑖=2𝜋𝜆_𝐿𝜇_𝐿𝜌_𝐿𝑉^∗_𝐿√𝑃𝑟, (11) with 𝑉^∗_𝐿=max(𝑉_𝐿𝑉,0.001),_𝐿=min(|𝑉_𝐿|,𝐶_𝜇^0.25𝐾_𝐿). (12) 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 𝛼_𝐿<0.1) taken as a “return to saturation” term, with a constant time scale 𝜏_𝐿 arbitrary equal to 1 second, and the weighting function 𝑓(𝛼_𝐿)=1−𝛼_𝐿=𝛼_𝐺.

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 𝑃𝑟(Re_𝑡)^7/16≫1. (13a) The heat transfer coefficient is 𝐶_𝐿𝑖=𝜌_𝐿𝐶_𝑝𝐿√𝑎_𝐿(𝜀_𝐿𝜈_𝐿)^1/4(𝛿⁴𝜀_𝐿𝜈³_𝐿)^1/24,𝑎_𝐿=𝜆_𝐿𝜌_𝐿𝐶_𝑝𝐿, (13b) 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: 𝑀_𝑘=𝛤_𝑘𝑉_𝑘+𝑀^𝐷_𝑘+𝑀^AM_𝑘+𝑀^𝐿_𝑘+𝑀^TD_𝑘, (14) 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
𝑀^𝐷_𝐺=−𝑀^𝐷_𝐿1=−8𝑎_𝑖𝜌_𝐿𝐶_𝐷|𝑉_𝐺−𝑉_𝐿|(𝑉_𝐺−𝑉_𝐿), (15) 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: 𝐶_𝐷=23𝑑₃₂𝑔|𝜌_𝐺−𝜌_𝐿|𝜎{1+17.67|𝑓(𝛼_𝐺)|^6/718.67𝑓(𝛼_𝐺)}²,𝑓(𝛼_𝐺)=(1−𝛼_𝐺)^1.5, (16) where 𝜎 is the surface tension. Equation (16) assumes that we are in the distorted bubble regime.
Added Mass Force is𝑀^AM_𝐺=−𝑀^AM_𝐿=−𝐶_AM1+2𝛼_𝐺1−𝛼_𝐺𝛼_𝐺𝜌_𝐿[(𝜕𝑉_𝐺𝜕𝑡+𝑉_𝐺.∇𝑉_𝐺)−(𝜕𝑉_𝐿𝜕𝑡+𝑉_𝐿.∇𝑉_𝐿)], (17) where 𝐶_AM is the added mass coefficient which is equal to 0.5 in the case of spherical bubbles.
Lift Force is𝑀^𝐿_𝐺=−𝑀^𝐿_𝐿=−𝐶_𝐿𝛼_𝐺𝜌_𝐿(𝑉_𝐺−𝑉_𝐿)(∇𝑉_𝐿), (18) 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𝑀^TD_𝐺=−𝑀^TD_𝐿=−𝐶_TD𝜌_𝐿𝐾_𝐿∇𝛼_𝐺, (19) where 𝐶_TD 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 𝛼_𝐺<0.7 (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 (𝛼_𝐺<0.7) by a factor of 10⁻⁴. 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.

Figure 1 

Iguchi’s experiment: fluid domain.

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

Table 1 

Geometry and flow conditions.

From the experiments, the authors found that in this configuration the jet produces a significant amount of small bubbles (𝑑<0.001m 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 74×214 cells in the vessel and 5×102 in the pipe).

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

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 10cm<𝑧<30cm), with some underestimation (see Figure 2).

(a)

(b)

Figure 2 

Calculated and measured values on the centreline.

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 𝑢^′_rms and 𝑣^′_rms, failing to catch the anisotropy of the problem (see Figure 3).

(a)

(b)

Figure 3 

Radial distribution of the root-mean-square values of turbulence components.

This is a classical feature of the 𝑘-𝜀 model; important differences are observed especially near the bath surface (𝑧=10cm) for the 𝑢^′_rms prediction. Nevertheless, moving downstream (𝑧>20cm), 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 𝐶_𝜀3, 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)

(b)

Figure 4 

Radial distribution of the root-mean-square values of the axial turbulence component (TRC term sensitivity analysis).

A sensitivity analysis to mesh refinement was conducted for the best performing turbulence model, considering three successively refined grids (CASE1, 104×212 cells; CASE2, 52×106 cells; CASE3, 26×53 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.

(a)

(b)

(c)

Figure 5 

Mesh sensitivity study (CASE1 = finer mesh, CASE2 = mean mesh, CASE3 = coarser mesh): radial velocity profiles considering the modified 𝑘-𝜀 model together with (a) separated phases model or (b) the Ishii’s correlation for 𝐶𝐷, and (c) adding the nondrag forces (AM, 𝐿, TD) contribution.

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

(a)

(b)

Figure 6 

Radial distribution of rms values of the turbulence components-comparison between single- and two-phase results.

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.

Figure 7 

2D simulation domain.

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

Table 2 

Flow regime parameters and measured temperature profiles.

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 𝑡=30 seconds and 𝑡=40 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)).

(a)

(b)

Figure 8 

Water temperature profiles for (a) standard (“St”) and (b) modified condensation models (“cont” = continue approach, “rt” = relaxation time, “discr” = discrete approach).

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

Figure 9 

Water temperature profiles for HD2 model (discreet approach) mesh sensitivity study.

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 (∼70%), while increasing mesh refinement in height direction caused a reduction (∼15%); for mesh refinement in both directions the sensitivity was reduced (∼30%).

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.

(a)

(b)

Figure 10 

Water temperature profiles: inlet turbulence sensitivity study.

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.

Acknowledgments

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.