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

Authors: Maria Cristina Galassi; Pierre Coste; Christophe Morel; Fabio Moretti
Journal: Science and Technology of Nuclear Installations
 (2009)

---

## 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 turbulentk-ε 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.


---

## Body

## 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)∂αkρk∂t+∇⋅(αkρkV¯k)=Γk,k=L,G, where the quantities αk,ρk and V¯k are the averaged fraction of
presence, the averaged density, and the averaged centre of mass velocity for
phase k, 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 Γk, is the rate of phase change (evaporation or condensation) per unit
volume of mixture.(ii) Two momentum balance equations are(2)∂αkρkV¯k∂t+∇⋅(αkρkV¯kV¯k)=−αk∇¯p+M¯k+αkρkg¯+∇⋅[αk(τ¯¯k+τ¯¯kT)],k=L,G, where p is the averaged pressure, common to the two phases, and g is the gravity acceleration vector. The two tensors τ¯¯kandτ¯¯kT are the averaged molecular and turbulent
Reynolds stress tensors, respectively, and the vector M¯k is the averaged interfacial momentum transfer between phases.(iii) Two total enthalpy balance equations are(3)∂∂t[αkρk(hk+Vk22)]+∇⋅(αkρk(hk+Vk22)V¯k)=αk∂p∂t+αkρkg¯⋅V¯k+Γk(hki+Vk22)+qki′′ai+qwk′′′−∇⋅[αk(q¯k+q¯kT)],k=L,G, where hk is the bulk-averaged enthalpy of phase k and hki is
the interfacial-weighted averaged enthalpy. The two vectors q¯kandq¯kT are the molecular and turbulent heat fluxes.
The term qki′′ai is the heat flux exchanged between phase k and the interface per unit volume, where ai is the
interfacial area concentration (or interfacial area per unit volume), and the
term qwk′′′ 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-equationk−ε 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: (4)ρk[∂Kk∂t+Vk,i∂Kk∂xi]=1αk∂∂xj[αkμkTσK∂Kk∂xj]+ρk(PRODk+Gk−εk)+PKi,(5)ρk[∂εk∂t+Vk,i∂εk∂xi]=1αk∂∂xj[αkμkTσε∂εk∂xj]+ρkεkKk(Cε1PRODk+Cε1max(Gk,0)−Cε2εk)Pεi, where the turbulent viscosity is given by μkT=Cμρk(Kk2/εk) with Cμ=0.09.PRODk  represents the (positive) production due to the mean velocity gradients,
and Gk is a stratification
attenuation term modelling the correlation between fluctuating densities and
velocity (more details are available in [5]). Cε1=1.44 and Cε2=1.92 are the two classical
constants taken from the single-phase model.By thek-ε 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 termsPKi and Pεi 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 termsPKi and Pεi are given by
(6a)PKi=−(M¯GD+M¯GAM)⋅(V¯G−V¯L),Pεi=Cε3PKiτ,τ=(d322εL)1/3, where M¯GD and M¯GAM are, respectively, the averaged drag and added
mass forces, Cε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 d32 and the turbulence dissipation rate εL.
It represents the time scale of liquid eddies having the same size as the
bubble diameter d32.“Ke EDF”
The difference to notice with the previous one is that the dissipation
production term is(6b)Pεi=Cε4εkKkPKi.“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 thek−ε 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] (7)∑k(Γkhki+qki′′ai)≈0. This important relation (together with the mass jump condition ΓG=−ΓL) allows to compute the mass-transfer term as a
function of the two volumetric heat fluxes qki′′ai and the interfacial-averaged enthalpies hki.
We assume that the interfacial-averaged enthalpies are identical to the phase-averaged
enthalpies (hki=hk). Each interfacial heat transfer term qki′′ai is the product of the interfacial heat flux
density which is expressed as (8)qki′′=Cki(Tsat(p)−Tk), where Cki is a heat transfer coefficient, Tk and Tsat(p) 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(9a)ai=6αGd32.
In stratified flow, with the “continuous approach,” it is (9b)ai=|∇αG|. In stratified flow, with the “discrete approach,”  it is (9c)ai=1hn→,
where hn→ is the cell size equal to the length of
the segment which includes the gravity centre of the cell and which has the
direction n→ normal to the free surface.A two-phase liquid-vapour flow is considered, where the liquid is
identified by index  “L” and the vapour by index “G”. qGi′′ model is not relevant for tests
presented here: air-water or  saturated
vapour. qLi′′ is calculated with (10)qLi′′=CLi(Tsat−TL). Two models for the heat exchange coefficient CLi, 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 ai is nonzero only near the interface, and
consequently the calculated heat transfer tends rapidly to zero elsewhere.Neptune_CFD 2004 [19] (HD1)
(11)CLi=2πλLμLρLVL*Pr, with (12)VL*=max(VL,0.001),VL=min(|V¯L|,Cμ0.25KL). 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 αL<0.1) taken as a
“return to saturation” term, with a constant time scale τL arbitrary equal to 1 second, and the weighting function f(αL)=1−αL=αG.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
(13a)Pr(Ret)7/16≫1. The heat transfer coefficient is (13b)CLi=ρLCpLaL(εLνL)1/4(δ4εLνL3)1/24,aL=λLρLCpL,
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 momentumM¯k appearing in (2) is
assumed to be the sum of five forces: (14)M¯k=ΓkV¯k+M¯kD+M¯kAM+M¯kL+M¯kTD, where we have assumed that the
interfacial-averaged velocity is equal to the phase-averaged velocity (V¯ki=V¯k), 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
(15)M¯GD=−M¯LD=−18aiρLCD|V¯G−V¯L|(V¯G−V¯L), where CD is the
nondimensional drag coefficient.
Concerning the drag coefficient closure law, theseparated 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:(16)CD=23d32g|ρG−ρL|σ{1+17.67|f(αG)|6/718.67f(αG)}2,f(αG)=(1−αG)1.5, where σ is the surface tension. Equation
(16) assumes that we are in the distorted bubble regime.
Added
Mass Force is(17)M¯GAM=−M¯LAM=−CAM1+2αG1−αGαGρL[(∂V¯G∂t+V¯G.∇¯¯V¯G)−(∂V¯L∂t+V¯L.∇¯¯V¯L)], where CAM is the added mass coefficient which is equal
to 0.5 in the case of spherical bubbles.
Lift
Force is(18)M¯GL=−M¯LL=−CLαGρL(V¯G−V¯L)⋀(∇¯⋀V¯L), where CL 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(19)M¯GTD=−M¯LTD=−CTDρLKL∇¯αG, where CTD 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αG<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 αG. 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 (αG<0.7) by a factor of 10−4. 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.


## 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)∂αkρk∂t+∇⋅(αkρkV¯k)=Γk,k=L,G, where the quantities αk,ρk and V¯k are the averaged fraction of
presence, the averaged density, and the averaged centre of mass velocity for
phase k, 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 Γk, is the rate of phase change (evaporation or condensation) per unit
volume of mixture.(ii) Two momentum balance equations are(2)∂αkρkV¯k∂t+∇⋅(αkρkV¯kV¯k)=−αk∇¯p+M¯k+αkρkg¯+∇⋅[αk(τ¯¯k+τ¯¯kT)],k=L,G, where p is the averaged pressure, common to the two phases, and g is the gravity acceleration vector. The two tensors τ¯¯kandτ¯¯kT are the averaged molecular and turbulent
Reynolds stress tensors, respectively, and the vector M¯k is the averaged interfacial momentum transfer between phases.(iii) Two total enthalpy balance equations are(3)∂∂t[αkρk(hk+Vk22)]+∇⋅(αkρk(hk+Vk22)V¯k)=αk∂p∂t+αkρkg¯⋅V¯k+Γk(hki+Vk22)+qki′′ai+qwk′′′−∇⋅[αk(q¯k+q¯kT)],k=L,G, where hk is the bulk-averaged enthalpy of phase k and hki is
the interfacial-weighted averaged enthalpy. The two vectors q¯kandq¯kT are the molecular and turbulent heat fluxes.
The term qki′′ai is the heat flux exchanged between phase k and the interface per unit volume, where ai is the
interfacial area concentration (or interfacial area per unit volume), and the
term qwk′′′ 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-equationk−ε 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: (4)ρk[∂Kk∂t+Vk,i∂Kk∂xi]=1αk∂∂xj[αkμkTσK∂Kk∂xj]+ρk(PRODk+Gk−εk)+PKi,(5)ρk[∂εk∂t+Vk,i∂εk∂xi]=1αk∂∂xj[αkμkTσε∂εk∂xj]+ρkεkKk(Cε1PRODk+Cε1max(Gk,0)−Cε2εk)Pεi, where the turbulent viscosity is given by μkT=Cμρk(Kk2/εk) with Cμ=0.09.PRODk  represents the (positive) production due to the mean velocity gradients,
and Gk is a stratification
attenuation term modelling the correlation between fluctuating densities and
velocity (more details are available in [5]). Cε1=1.44 and Cε2=1.92 are the two classical
constants taken from the single-phase model.By thek-ε 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 termsPKi and Pεi 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 termsPKi and Pεi are given by
(6a)PKi=−(M¯GD+M¯GAM)⋅(V¯G−V¯L),Pεi=Cε3PKiτ,τ=(d322εL)1/3, where M¯GD and M¯GAM are, respectively, the averaged drag and added
mass forces, Cε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 d32 and the turbulence dissipation rate εL.
It represents the time scale of liquid eddies having the same size as the
bubble diameter d32.“Ke EDF”
The difference to notice with the previous one is that the dissipation
production term is(6b)Pεi=Cε4εkKkPKi.“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 thek−ε transport equations
were not used for the gas phase.


## 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 termsPKi and Pεi are given by
(6a)PKi=−(M¯GD+M¯GAM)⋅(V¯G−V¯L),Pεi=Cε3PKiτ,τ=(d322εL)1/3, where M¯GD and M¯GAM are, respectively, the averaged drag and added
mass forces, Cε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 d32 and the turbulence dissipation rate εL.
It represents the time scale of liquid eddies having the same size as the
bubble diameter d32.“Ke EDF”
The difference to notice with the previous one is that the dissipation
production term is(6b)Pεi=Cε4εkKkPKi.“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 thek−ε 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] (7)∑k(Γkhki+qki′′ai)≈0. This important relation (together with the mass jump condition ΓG=−ΓL) allows to compute the mass-transfer term as a
function of the two volumetric heat fluxes qki′′ai and the interfacial-averaged enthalpies hki.
We assume that the interfacial-averaged enthalpies are identical to the phase-averaged
enthalpies (hki=hk). Each interfacial heat transfer term qki′′ai is the product of the interfacial heat flux
density which is expressed as (8)qki′′=Cki(Tsat(p)−Tk), where Cki is a heat transfer coefficient, Tk and Tsat(p) 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(9a)ai=6αGd32.
In stratified flow, with the “continuous approach,” it is (9b)ai=|∇αG|. In stratified flow, with the “discrete approach,”  it is (9c)ai=1hn→,
where hn→ is the cell size equal to the length of
the segment which includes the gravity centre of the cell and which has the
direction n→ normal to the free surface.A two-phase liquid-vapour flow is considered, where the liquid is
identified by index  “L” and the vapour by index “G”. qGi′′ model is not relevant for tests
presented here: air-water or  saturated
vapour. qLi′′ is calculated with (10)qLi′′=CLi(Tsat−TL). Two models for the heat exchange coefficient CLi, 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 ai is nonzero only near the interface, and
consequently the calculated heat transfer tends rapidly to zero elsewhere.Neptune_CFD 2004 [19] (HD1)
(11)CLi=2πλLμLρLVL*Pr, with (12)VL*=max(VL,0.001),VL=min(|V¯L|,Cμ0.25KL). 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 αL<0.1) taken as a
“return to saturation” term, with a constant time scale τL arbitrary equal to 1 second, and the weighting function f(αL)=1−αL=αG.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
(13a)Pr(Ret)7/16≫1. The heat transfer coefficient is (13b)CLi=ρLCpLaL(εLνL)1/4(δ4εLνL3)1/24,aL=λLρLCpL,
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 momentumM¯k appearing in (2) is
assumed to be the sum of five forces: (14)M¯k=ΓkV¯k+M¯kD+M¯kAM+M¯kL+M¯kTD, where we have assumed that the
interfacial-averaged velocity is equal to the phase-averaged velocity (V¯ki=V¯k), 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
(15)M¯GD=−M¯LD=−18aiρLCD|V¯G−V¯L|(V¯G−V¯L), where CD is the
nondimensional drag coefficient.
Concerning the drag coefficient closure law, theseparated 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:(16)CD=23d32g|ρG−ρL|σ{1+17.67|f(αG)|6/718.67f(αG)}2,f(αG)=(1−αG)1.5, where σ is the surface tension. Equation
(16) assumes that we are in the distorted bubble regime.
Added
Mass Force is(17)M¯GAM=−M¯LAM=−CAM1+2αG1−αGαGρL[(∂V¯G∂t+V¯G.∇¯¯V¯G)−(∂V¯L∂t+V¯L.∇¯¯V¯L)], where CAM is the added mass coefficient which is equal
to 0.5 in the case of spherical bubbles.
Lift
Force is(18)M¯GL=−M¯LL=−CLαGρL(V¯G−V¯L)⋀(∇¯⋀V¯L), where CL 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(19)M¯GTD=−M¯LTD=−CTDρLKL∇¯αG, where CTD 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αG<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 αG. 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 (αG<0.7) by a factor of 10−4. 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.


## 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] (7)∑k(Γkhki+qki′′ai)≈0. This important relation (together with the mass jump condition ΓG=−ΓL) allows to compute the mass-transfer term as a
function of the two volumetric heat fluxes qki′′ai and the interfacial-averaged enthalpies hki.
We assume that the interfacial-averaged enthalpies are identical to the phase-averaged
enthalpies (hki=hk). Each interfacial heat transfer term qki′′ai is the product of the interfacial heat flux
density which is expressed as (8)qki′′=Cki(Tsat(p)−Tk), where Cki is a heat transfer coefficient, Tk and Tsat(p) 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(9a)ai=6αGd32.
In stratified flow, with the “continuous approach,” it is (9b)ai=|∇αG|. In stratified flow, with the “discrete approach,”  it is (9c)ai=1hn→,
where hn→ is the cell size equal to the length of
the segment which includes the gravity centre of the cell and which has the
direction n→ normal to the free surface.A two-phase liquid-vapour flow is considered, where the liquid is
identified by index  “L” and the vapour by index “G”. qGi′′ model is not relevant for tests
presented here: air-water or  saturated
vapour. qLi′′ is calculated with (10)qLi′′=CLi(Tsat−TL). Two models for the heat exchange coefficient CLi, 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 ai is nonzero only near the interface, and
consequently the calculated heat transfer tends rapidly to zero elsewhere.Neptune_CFD 2004 [19] (HD1)
(11)CLi=2πλLμLρLVL*Pr, with (12)VL*=max(VL,0.001),VL=min(|V¯L|,Cμ0.25KL). 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 αL<0.1) taken as a
“return to saturation” term, with a constant time scale τL arbitrary equal to 1 second, and the weighting function f(αL)=1−αL=αG.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
(13a)Pr(Ret)7/16≫1. The heat transfer coefficient is (13b)CLi=ρLCpLaL(εLνL)1/4(δ4εLνL3)1/24,aL=λLρLCpL,
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 momentumM¯k appearing in (2) is
assumed to be the sum of five forces: (14)M¯k=ΓkV¯k+M¯kD+M¯kAM+M¯kL+M¯kTD, where we have assumed that the
interfacial-averaged velocity is equal to the phase-averaged velocity (V¯ki=V¯k), 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
(15)M¯GD=−M¯LD=−18aiρLCD|V¯G−V¯L|(V¯G−V¯L), where CD is the
nondimensional drag coefficient.
Concerning the drag coefficient closure law, theseparated 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:(16)CD=23d32g|ρG−ρL|σ{1+17.67|f(αG)|6/718.67f(αG)}2,f(αG)=(1−αG)1.5, where σ is the surface tension. Equation
(16) assumes that we are in the distorted bubble regime.
Added
Mass Force is(17)M¯GAM=−M¯LAM=−CAM1+2αG1−αGαGρL[(∂V¯G∂t+V¯G.∇¯¯V¯G)−(∂V¯L∂t+V¯L.∇¯¯V¯L)], where CAM is the added mass coefficient which is equal
to 0.5 in the case of spherical bubbles.
Lift
Force is(18)M¯GL=−M¯LL=−CLαGρL(V¯G−V¯L)⋀(∇¯⋀V¯L), where CL 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(19)M¯GTD=−M¯LTD=−CTDρLKL∇¯αG, where CTD 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αG<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 αG. 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 (αG<0.7) by a factor of 10−4. 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 Table1.Table 1
Geometry and flow conditions.
(a)
Geometry of test section
Circular injector, vertical downward flowVessel diameterD = 200 mmVessel heightH = 390 mmInjector diameterd = 5 mmHeight of the injector above free surfaceh = 2 mm(b)
Flow conditions
Turbulent nonfragmented jet entering a free surfaceFluidsAir-waterHeatingnoPressure1 barJet velocity2.54 m/sJet turbulence intensity5–15%From the experiments, the authors found that in this configuration the jet
produces a significant amount of small bubbles (d<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 (with74×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 (z=10cm) and bubble dispersion
in the whole water vessel (z=20 and 30 cm), but with very small void
fraction values (α∼10−3).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<z<30cm), with some
underestimation (see Figure 2).Calculated and measured values on
the centreline.
(a)(b)Both thek-ε models considered
for water (the standard and modified versions) underestimated the turbulence
production and predicted almost the same radial profiles for the turbulence
components urms′ and vrms′, failing to catch the
anisotropy of the problem (see Figure 3).Radial distribution of the
root-mean-square values of turbulence components.
(a)(b)This is a classical feature of thek-ε model; important differences are observed especially near the
bath surface (z=10cm) for the urms′ prediction.
Nevertheless, moving downstream (z>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), PKi and Pεi,
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, d32, and the parameter Cε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.Radial
distribution of the root-mean-square values of the axial turbulence component
(TRC term sensitivity analysis).
(a)(b)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 (L), 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.Mesh sensitivity
study (CASE1 = finer mesh, CASE2 = mean mesh, CASE3 = coarser mesh): radial
velocity profiles considering  the
modifiedk-ε model together with (a)
separated phases model or (b) the Ishii’s correlation for CD,
and  (c) adding the nondrag forces (AM,
L, TD) contribution.
(a)(b)(c)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 Figure6).Radial distribution of rms values
of the turbulence components-comparison
between single- and two-phase results.
(a)(b)Even
if there is some confidence in the predicting capabilities ofk−ε 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 Table2. 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.
(a)Flow regime parametersInlet flow rate velocity of water0.28Inlet flow rate velocity of steam3.2 m/sInlet water temperature300.2 KMeasured mean water temperature324.43 KPressure0.697 MPaHeight of water layer0.031 mMeasured steam velocity1.94 m/sCondensation rate∼40%(b)y [mm]T [K]4310.948310.4812311.3116313.1420316.4224322.2626331.7928354.7330411.9731435.6734434.7736435.74Calculations 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 betweent=30 seconds and t=40 seconds.Simulations were run testing at first standard
two-phase models (see Figure8(a)) and then the new implemented models (see Figure
8(b)).Water temperature profiles for
(a) standard (“St”) and (b) modified condensation models (“cont” = continue
approach, “rt” = relaxation time, “discr” = discrete approach).
(a)(b)Results showed that modifiedk−ε 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
modelsHD1 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]).TheHD2 (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 forHD2 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 Figure10(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.Water
temperature profiles: inlet turbulence sensitivity study.
(a)(b)


## 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 Table1.Table 1
Geometry and flow conditions.
(a)
Geometry of test section
Circular injector, vertical downward flowVessel diameterD = 200 mmVessel heightH = 390 mmInjector diameterd = 5 mmHeight of the injector above free surfaceh = 2 mm(b)
Flow conditions
Turbulent nonfragmented jet entering a free surfaceFluidsAir-waterHeatingnoPressure1 barJet velocity2.54 m/sJet turbulence intensity5–15%From the experiments, the authors found that in this configuration the jet
produces a significant amount of small bubbles (d<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 (with74×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 (z=10cm) and bubble dispersion
in the whole water vessel (z=20 and 30 cm), but with very small void
fraction values (α∼10−3).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<z<30cm), with some
underestimation (see Figure 2).Calculated and measured values on
the centreline.
(a)(b)Both thek-ε models considered
for water (the standard and modified versions) underestimated the turbulence
production and predicted almost the same radial profiles for the turbulence
components urms′ and vrms′, failing to catch the
anisotropy of the problem (see Figure 3).Radial distribution of the
root-mean-square values of turbulence components.
(a)(b)This is a classical feature of thek-ε model; important differences are observed especially near the
bath surface (z=10cm) for the urms′ prediction.
Nevertheless, moving downstream (z>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), PKi and Pεi,
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, d32, and the parameter Cε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.Radial
distribution of the root-mean-square values of the axial turbulence component
(TRC term sensitivity analysis).
(a)(b)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 (L), 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.Mesh sensitivity
study (CASE1 = finer mesh, CASE2 = mean mesh, CASE3 = coarser mesh): radial
velocity profiles considering  the
modifiedk-ε model together with (a)
separated phases model or (b) the Ishii’s correlation for CD,
and  (c) adding the nondrag forces (AM,
L, TD) contribution.
(a)(b)(c)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 Figure6).Radial distribution of rms values
of the turbulence components-comparison
between single- and two-phase results.
(a)(b)Even
if there is some confidence in the predicting capabilities ofk−ε 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 Table2. 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.
(a)Flow regime parametersInlet flow rate velocity of water0.28Inlet flow rate velocity of steam3.2 m/sInlet water temperature300.2 KMeasured mean water temperature324.43 KPressure0.697 MPaHeight of water layer0.031 mMeasured steam velocity1.94 m/sCondensation rate∼40%(b)y [mm]T [K]4310.948310.4812311.3116313.1420316.4224322.2626331.7928354.7330411.9731435.6734434.7736435.74Calculations 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 betweent=30 seconds and t=40 seconds.Simulations were run testing at first standard
two-phase models (see Figure8(a)) and then the new implemented models (see Figure
8(b)).Water temperature profiles for
(a) standard (“St”) and (b) modified condensation models (“cont” = continue
approach, “rt” = relaxation time, “discr” = discrete approach).
(a)(b)Results showed that modifiedk−ε 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
modelsHD1 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]).TheHD2 (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 forHD2 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 Figure10(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.Water
temperature profiles: inlet turbulence sensitivity study.
(a)(b)


## 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 ak−ε 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 thatk−ε 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.



---

*Source: 950536-2008-11-06.xml*