# CFD Code Validation against Stratified Air-Water Flow Experimental Data

Authors: F. Terzuoli; M. C. Galassi; D. Mazzini; F. D'Auria
Journal: Science and Technology of Nuclear Installations
 (2008)

---

## Abstract

Pressurized thermal shock (PTS) modelling has been identified as one of the most important industrial needs related to nuclear reactor safety. A severe PTS scenario limiting the reactor pressure vessel (RPV) lifetime is the cold water emergency core cooling (ECC) injection into the cold leg during a loss of coolant accident (LOCA). Since it represents a big challenge for numerical simulations, this scenario was selected within the European Platform for Nuclear Reactor Simulations (NURESIM) Integrated Project as a reference two-phase problem for computational fluid dynamics (CFDs) code validation. This paper presents a CFD analysis of a stratified air-water flow experimental investigation performed at the Institut de Mécanique des Fluides de Toulouse in 1985, which shares some common physical features with the ECC injection in PWR cold leg. Numerical simulations have been carried out with two commercial codes (Fluent and Ansys CFX), and a research code (NEPTUNE CFD). The aim of this work, carried out at the University of Pisa within the NURESIM IP, is to validate the free surface flow model implemented in the codes against experimental data, and to perform code-to-code benchmarking. Obtained results suggest the relevance of three-dimensional effects and stress the importance of a suitable interface drag modelling.


---

## Body

## 1. Introduction

The  European
Platform for Nuclear Reactor Simulations (NURESIM) Integrated Project aims at
developing a  common European Standard Software Platform for modelling,
recording, and recovering computer simulation data for current and future
nuclear reactor systems [1]. NEPTUNE [2] is the thermal-hydraulics tool of NURESIM, and is
designed to simulate two-phase flow in all situations encountered in nuclear
reactor power plants. The present work is aimed at the validation and
benchmarking of NEPTUNE_CFD, the two-phase CFD tool of the NEPTUNE platform [3, 4].Since PTS has been
identified as one of the most important aspects related to nuclear reactor
safety, some relevant PTS scenarios were chosen as reference test cases for CFD
code validation within the NURESIM IP [5, 6]. A severe PTS scenario limiting the reactor pressure
vessel (RPV) lifetime is the cold water emergency core cooling (ECC) injection
into the cold leg during a loss of coolant accident (LOCA). Complex phenomena
take place during this scenario, such as turbulent mixing in the jet region and
downstream of the impingement zone, stratified two-phase flows, phase change at
the steam water interface. This paper deals with the study of a stratified
air-water flow experiment performed at the Institut de Mécanique des Fluides de
Toulouse in 1985 [7, 8]; this flow configuration is likely to share common
physical features with the chosen PTS scenario.To validate the
two-phase models implemented in NEPTUNE_CFD V1.0.6, numerical simulation
results were compared with both experimental data and predictions from two
commercial CFD codes, namely, ANSYS CFX 10.0 and FLUENT 6.1.


## 2. Experimental Facility and Tests

The experimental
facility (see Figure1) consists of a quasihorizontal channel with
rectangular cross-section (0.1 m high, 0.2 m wide, and 13.0 m long), with an inclination of 0.1° with respect
to the horizontal plane. The rectangular channel is connected to the water and
air inlet and outlet volumes. Desired mass flow rates are granted by a
recirculation pump and all the facility walls are adiabatic.Figure 1
Experimental facility—conceptual scheme.The facility is
equipped with sensors located at 7.05 m, 9.10 m, and 11.10 m from the inlet section, which provide
the measurements of mass flow rates, local instantaneous water height (by
capacitance wire sensors), local mean and fluctuating values of horizontal and
vertical velocity components, as well as Reynolds stress tensor (by a laser
Doppler anemometer in water and by hot wire anemometer in air). Documentation
is available [7] for three test cases conducted at ambient pressure
and temperature, characterized by constant water mass flow rate and different
air mass flow rates. This work deals only with one of these tests, namely, the
“T250” experiment, in which water and air bulk velocities are 0.395 m/s and
3.66 m/s, respectively. Under this condition, mean water depth is measured to
be 38 mm.


## 3. Experimental Test Simulation with Ansys Cfx and Fluent

### 3.1. Preliminary Results of Two-Phase Calculations

Since the width of
the duct is large compared to the height, a two-dimensional velocity field has
been assumed in order to perform preliminary calculations. These analyses, carried
out with the ANSYS CFX 10.0 code [9], adopted a 1-mm wide one-cell thick grid since it
does not allow assuming a real two-dimensional domain. A spatial discretization
of the channel has been created using ANSYS ICEM CFD 10.0 code, counting ~20 k
hexahedral elements (~40 k computational nodes). Elements refinement has been
provided near the walls and at the expected height of interface between fluids;
anyway it is worth noting that in more realistic applications, the interface is
moving in the domain, so that such grid refinement could be obtained only with
dynamic meshing adaptation.The “inhomogeneous”
two-phase flow model was selected, since some interface instability has been
observed in all preliminary simulations. This model solves one velocity field for
each phase (resulting in two different coexisting fields in the domain); while
the “homogeneous” setting has been adopted for turbulence together with thek−ω
shear-stress transport (SST) model, providing only one field shared by the
phases. Moreover, the code default interface treatment model has been used, the
so-called “standard free surface” [9].A uniform velocity
for both air and water has been assumed at the inlet. Air constant pressure and
water hydrostatic profile have been imposed at the outlet section according to measured
water height. Finally, adiabatic upper and lower walls with no slip and symmetric
lateral faces of the domain have been selected.Figure2 shows the calculated velocity in the test section
compared with experimental data. Air-water stratification has been correctly
predicted, but relevant differences can be observed: water level is calculated
about 24% lower, while the maximum air velocity is reduced by 10% and it is no
longer placed between the top wall and the interface (66 mm), but closer to the
wall (81 mm,
~20% higher). These relevant mismatches suggest that the modelling of the
frictional drag between the phases is overestimated, and a deeper understanding
of the experimental data is needed. In order to investigate these problems, single-effect
analyses have been carried out as described in the following Sections, together with some
sensitivity analyses on the most relevant parameter.Preliminary results—velocity profiles: (a) water; (b) air.
(a)(b)


### 3.2. Experimental Data Understanding

All performed
experimental tests assume the same value of water flow rate. Different water
heights are measured into the channel depending on the different equilibrium
conditions between the forces acting on the fluids, that is, the drag force
between air and water, the longitudinal component of the gravitational force
due to the inclination of the duct, and the friction forces acting on walls. Except
for the gravitational force, the others depend on fluid velocity, thus changing
their values flowing into the channel: the drag force increases while the
interface friction force decreases up to reaching the equilibrium condition. An
incorrect prediction of one of them can justify the mismatch between calculated
profiles and experimental data.Furthermore, it is
worth noting that water average velocity resulting from measured data at probe
section (calculated by means of the trapezes rule, thus underestimating the
real value) is >4% greater than the value provided in [7] for the water bulk velocity. The same occurs for the
air. Since the flow rate is constant along the channel, this mismatch could be
due to the development of a 3D profile. In fact, in a real 3D geometry, the
velocity profile measured at the symmetry plane (as shown in all the plots) is
the maximum profile of a developed 3D flow (see Figure 3). It is possible to conclude that considering a 2D
computational domain implies a loss of information related to the flow
development.Figure 3
Three-dimensional water domain.


### 3.3. Single-Phase Analysis

In this analysis,
the computational domain was splitted into two separate parts along the experimental
interface level. Spatial discretizations have been created for each phase
channel in both two- and three-dimensional configurations. Different node
distributions have been employed to evaluate the grid requirements for CFD
simulation to correctly reproduce near-wall effect and flow developing. The
most relevant grid details are reported in Table1 for both FLUENT 6.1 [10] and CFX 10.0 [9] codes. These characteristics have been established
following the main findings of preliminary mesh sensitivity investigations.Table 1
Details of grids for single-phase analysis.
CodeDomain2D3DFLUENT 6.1Air320×35 and 640×70320×35×40Water320×27 and 320×54320×27×40CFX
10.0Air520×34520×34×38Water520×23520×23×44In  single-phase
calculations, the interface has been modelled as a wall moving at the expected
velocity of the free surface. Since this value is not available, it has been
imposed in the range 0.50–0.65 m/s,
according to the top measured water velocity, 0.502 m/s, which is the available
data closer to the interface. The boundary conditions have been imposed according
to the preliminary calculation documented in Section3.1.Direct numerical
simulation (DNS) and large eddy simulation (LES) of stratified flows [11] were used in the framework of the NURESIM Integrated
Project to derive some closure laws for interfacial momentum, turbulence, and
heat transfer. Future work is still necessary to implement these laws in CFD
codes as well as to compare predictions with DNS-LES studies on the same flow
conditions [12] and to validate them against experimental data. Anyway,
this subject is beyond the aim of the present article.


#### 3.3.1. Single-Phase Water Flow

In Figure4, the obtained results are shown in terms of longitudinal
water velocity profile at the test section located at 9.1 mm from the inlet. The
experimental profile is correctly predicted from a qualitative point of view by
both CFX and FLUENT codes. No relevant improvements are obtained varying the
interface velocity except for the region closer to the moving wall in both two-
and three-dimensional calculations. However, two-dimensional calculations
underestimate the velocity values by about 10% than three-dimensional ones.
Thus, domain three dimensionality has a great relevance on water velocity
profile and cannot be neglected in the simulations. Moreover, grid sensitivity
analysis has shown that limited improvement is obtained by increasing the
number of cells.Water velocity: (a) 2D grid with
different interface velocities; (b) 3D grid with different interface
velocities; (c) 3D versus 2D grids.
(a)(b)(c)Figure5 shows the comparison between calculated turbulent
kinetic energy and experimental data both for two- and three-dimensional
simulations. The third dimension has great relevance on near-wall values,
increasing them by about 20%. Calculated values in the region near the
interface have a different behavior from experimental data owing to the
presence of a solid wall. Turbulence produced by the contact between fluids at
the interface is not included in the model. Although these differences in shape
and local values, predicted turbulent kinetic energy is in good agreement with
measurements, especially for three-dimensional calculations.Water turbulent kinetic energy: (a)
two-dimensional simulations;
(b) three-dimensional simulations.
(a)(b)


#### 3.3.2. Single-Phase Air Flow

Figure6(a) shows the transversal air velocity profiles at the
test section in three-dimensional simulations calculated with CFX, which is
predicted with relevant differences on both shape and values. Negligible
effects on the results are obtained varying the interface velocity.
Two-dimensional calculations show the same behavior as well as results obtained
using FLUENT. Moreover, as shown in Figure 6(b), the experimental data provide a greater maximum
velocity with respect to both two-dimensional (15.8% higher), and
three-dimensional simulations (11.4% higher). This is a further confirmation of
the three-dimensional flow structure supposed in Section 3.2. However, the systematic underestimation of the
velocity needs further investigations.Air velocity: (a) 3D grid with different interface velocities; (b) 3D versus 2D.
(a)(b)Finally, from
sensitivity analysis performed with FLUENT, no relevant advantages are obtained
by increasing the cell number or changing the turbulence near-wall treatment.


## 3.1. Preliminary Results of Two-Phase Calculations

Since the width of
the duct is large compared to the height, a two-dimensional velocity field has
been assumed in order to perform preliminary calculations. These analyses, carried
out with the ANSYS CFX 10.0 code [9], adopted a 1-mm wide one-cell thick grid since it
does not allow assuming a real two-dimensional domain. A spatial discretization
of the channel has been created using ANSYS ICEM CFD 10.0 code, counting ~20 k
hexahedral elements (~40 k computational nodes). Elements refinement has been
provided near the walls and at the expected height of interface between fluids;
anyway it is worth noting that in more realistic applications, the interface is
moving in the domain, so that such grid refinement could be obtained only with
dynamic meshing adaptation.The “inhomogeneous”
two-phase flow model was selected, since some interface instability has been
observed in all preliminary simulations. This model solves one velocity field for
each phase (resulting in two different coexisting fields in the domain); while
the “homogeneous” setting has been adopted for turbulence together with thek−ω
shear-stress transport (SST) model, providing only one field shared by the
phases. Moreover, the code default interface treatment model has been used, the
so-called “standard free surface” [9].A uniform velocity
for both air and water has been assumed at the inlet. Air constant pressure and
water hydrostatic profile have been imposed at the outlet section according to measured
water height. Finally, adiabatic upper and lower walls with no slip and symmetric
lateral faces of the domain have been selected.Figure2 shows the calculated velocity in the test section
compared with experimental data. Air-water stratification has been correctly
predicted, but relevant differences can be observed: water level is calculated
about 24% lower, while the maximum air velocity is reduced by 10% and it is no
longer placed between the top wall and the interface (66 mm), but closer to the
wall (81 mm,
~20% higher). These relevant mismatches suggest that the modelling of the
frictional drag between the phases is overestimated, and a deeper understanding
of the experimental data is needed. In order to investigate these problems, single-effect
analyses have been carried out as described in the following Sections, together with some
sensitivity analyses on the most relevant parameter.Preliminary results—velocity profiles: (a) water; (b) air.
(a)(b)


## 3.2. Experimental Data Understanding

All performed
experimental tests assume the same value of water flow rate. Different water
heights are measured into the channel depending on the different equilibrium
conditions between the forces acting on the fluids, that is, the drag force
between air and water, the longitudinal component of the gravitational force
due to the inclination of the duct, and the friction forces acting on walls. Except
for the gravitational force, the others depend on fluid velocity, thus changing
their values flowing into the channel: the drag force increases while the
interface friction force decreases up to reaching the equilibrium condition. An
incorrect prediction of one of them can justify the mismatch between calculated
profiles and experimental data.Furthermore, it is
worth noting that water average velocity resulting from measured data at probe
section (calculated by means of the trapezes rule, thus underestimating the
real value) is >4% greater than the value provided in [7] for the water bulk velocity. The same occurs for the
air. Since the flow rate is constant along the channel, this mismatch could be
due to the development of a 3D profile. In fact, in a real 3D geometry, the
velocity profile measured at the symmetry plane (as shown in all the plots) is
the maximum profile of a developed 3D flow (see Figure 3). It is possible to conclude that considering a 2D
computational domain implies a loss of information related to the flow
development.Figure 3
Three-dimensional water domain.


## 3.3. Single-Phase Analysis

In this analysis,
the computational domain was splitted into two separate parts along the experimental
interface level. Spatial discretizations have been created for each phase
channel in both two- and three-dimensional configurations. Different node
distributions have been employed to evaluate the grid requirements for CFD
simulation to correctly reproduce near-wall effect and flow developing. The
most relevant grid details are reported in Table1 for both FLUENT 6.1 [10] and CFX 10.0 [9] codes. These characteristics have been established
following the main findings of preliminary mesh sensitivity investigations.Table 1
Details of grids for single-phase analysis.
CodeDomain2D3DFLUENT 6.1Air320×35 and 640×70320×35×40Water320×27 and 320×54320×27×40CFX
10.0Air520×34520×34×38Water520×23520×23×44In  single-phase
calculations, the interface has been modelled as a wall moving at the expected
velocity of the free surface. Since this value is not available, it has been
imposed in the range 0.50–0.65 m/s,
according to the top measured water velocity, 0.502 m/s, which is the available
data closer to the interface. The boundary conditions have been imposed according
to the preliminary calculation documented in Section3.1.Direct numerical
simulation (DNS) and large eddy simulation (LES) of stratified flows [11] were used in the framework of the NURESIM Integrated
Project to derive some closure laws for interfacial momentum, turbulence, and
heat transfer. Future work is still necessary to implement these laws in CFD
codes as well as to compare predictions with DNS-LES studies on the same flow
conditions [12] and to validate them against experimental data. Anyway,
this subject is beyond the aim of the present article.


### 3.3.1. Single-Phase Water Flow

In Figure4, the obtained results are shown in terms of longitudinal
water velocity profile at the test section located at 9.1 mm from the inlet. The
experimental profile is correctly predicted from a qualitative point of view by
both CFX and FLUENT codes. No relevant improvements are obtained varying the
interface velocity except for the region closer to the moving wall in both two-
and three-dimensional calculations. However, two-dimensional calculations
underestimate the velocity values by about 10% than three-dimensional ones.
Thus, domain three dimensionality has a great relevance on water velocity
profile and cannot be neglected in the simulations. Moreover, grid sensitivity
analysis has shown that limited improvement is obtained by increasing the
number of cells.Water velocity: (a) 2D grid with
different interface velocities; (b) 3D grid with different interface
velocities; (c) 3D versus 2D grids.
(a)(b)(c)Figure5 shows the comparison between calculated turbulent
kinetic energy and experimental data both for two- and three-dimensional
simulations. The third dimension has great relevance on near-wall values,
increasing them by about 20%. Calculated values in the region near the
interface have a different behavior from experimental data owing to the
presence of a solid wall. Turbulence produced by the contact between fluids at
the interface is not included in the model. Although these differences in shape
and local values, predicted turbulent kinetic energy is in good agreement with
measurements, especially for three-dimensional calculations.Water turbulent kinetic energy: (a)
two-dimensional simulations;
(b) three-dimensional simulations.
(a)(b)


### 3.3.2. Single-Phase Air Flow

Figure6(a) shows the transversal air velocity profiles at the
test section in three-dimensional simulations calculated with CFX, which is
predicted with relevant differences on both shape and values. Negligible
effects on the results are obtained varying the interface velocity.
Two-dimensional calculations show the same behavior as well as results obtained
using FLUENT. Moreover, as shown in Figure 6(b), the experimental data provide a greater maximum
velocity with respect to both two-dimensional (15.8% higher), and
three-dimensional simulations (11.4% higher). This is a further confirmation of
the three-dimensional flow structure supposed in Section 3.2. However, the systematic underestimation of the
velocity needs further investigations.Air velocity: (a) 3D grid with different interface velocities; (b) 3D versus 2D.
(a)(b)Finally, from
sensitivity analysis performed with FLUENT, no relevant advantages are obtained
by increasing the cell number or changing the turbulence near-wall treatment.


## 3.3.1. Single-Phase Water Flow

In Figure4, the obtained results are shown in terms of longitudinal
water velocity profile at the test section located at 9.1 mm from the inlet. The
experimental profile is correctly predicted from a qualitative point of view by
both CFX and FLUENT codes. No relevant improvements are obtained varying the
interface velocity except for the region closer to the moving wall in both two-
and three-dimensional calculations. However, two-dimensional calculations
underestimate the velocity values by about 10% than three-dimensional ones.
Thus, domain three dimensionality has a great relevance on water velocity
profile and cannot be neglected in the simulations. Moreover, grid sensitivity
analysis has shown that limited improvement is obtained by increasing the
number of cells.Water velocity: (a) 2D grid with
different interface velocities; (b) 3D grid with different interface
velocities; (c) 3D versus 2D grids.
(a)(b)(c)Figure5 shows the comparison between calculated turbulent
kinetic energy and experimental data both for two- and three-dimensional
simulations. The third dimension has great relevance on near-wall values,
increasing them by about 20%. Calculated values in the region near the
interface have a different behavior from experimental data owing to the
presence of a solid wall. Turbulence produced by the contact between fluids at
the interface is not included in the model. Although these differences in shape
and local values, predicted turbulent kinetic energy is in good agreement with
measurements, especially for three-dimensional calculations.Water turbulent kinetic energy: (a)
two-dimensional simulations;
(b) three-dimensional simulations.
(a)(b)


## 3.3.2. Single-Phase Air Flow

Figure6(a) shows the transversal air velocity profiles at the
test section in three-dimensional simulations calculated with CFX, which is
predicted with relevant differences on both shape and values. Negligible
effects on the results are obtained varying the interface velocity.
Two-dimensional calculations show the same behavior as well as results obtained
using FLUENT. Moreover, as shown in Figure 6(b), the experimental data provide a greater maximum
velocity with respect to both two-dimensional (15.8% higher), and
three-dimensional simulations (11.4% higher). This is a further confirmation of
the three-dimensional flow structure supposed in Section 3.2. However, the systematic underestimation of the
velocity needs further investigations.Air velocity: (a) 3D grid with different interface velocities; (b) 3D versus 2D.
(a)(b)Finally, from
sensitivity analysis performed with FLUENT, no relevant advantages are obtained
by increasing the cell number or changing the turbulence near-wall treatment.


## 4. T250 Experimental Test Simulation with Neptune_cfd

In the hypothesis
of planar symmetry, the computational domain is constituted by a 2D section of
the channel; it has been modelled with three successively refined 2D grids built
up with the ANSYS ICEM code: 75×30, 100×40 and 150×60 cells, respectively.
All imposed boundary conditions were as for CFX calculations (Section 3.1); initially a stratified air/water flow was
established and parabolic velocity profiles were imposed for both water and air
flows at inlet. Calculations were run with NEPTUNE_CFD V1.0.6 by means of an
input deck kindly provided by Mr. Pierre COSTE (CEA/Grenoble). The k−ε model was considered for both air
and water turbulence together with interface turbulence production, “Pierre Coste Large Interface Model” [13]
was selected
for the drag coefficient.As Figure7 shows, water velocity profile is quite well predicted
in all three cases (maximum error ~7%), while air velocity profile is
appreciably underestimated, especially in the bulk region between wall and free
surface (maximum error ~12% for the coarser grid and ~10% for the finer one).
It seems that mesh refinement does not produce important improvements, except for
the air velocity profile in the region near the interface.Velocity profiles: (a) water; (b) air.
(a)(b)Figures8(a) and 8(b) show the turbulent kinetic energy profile at probe
location for water and air flow, respectively. As in the previous case, the
profile is qualitatively well predicted for the water flow, while for the air
flow, it is significantly underestimated in the bottom region and overestimated
in the upper one. It can be observed that calculated values with refined grids
better match experimental data in the upper region of the water stream, with
some underestimation near the wall (maximum error ~45%). It is worth noting
that the code is able to catch the increase of water turbulence near the free
surface due to the air stream, but not due to the presence of solid walls. On
the contrary, air turbulent kinetic energy profile does not get significantly
better with mesh refinement; results are underestimated near the interface
(maximum error ~66%) and slightly overestimated elsewhere.Turbulent kinetic energy: (a) water; (b) air.
(a)(b)Calculations were
also run considering the “separated phases model” for the interphase drag
coefficient. The resulting velocity profiles seem to be very similar to that
predicted by ANSYS CFX and shown in Figure2. In both cases, the interface level is underestimated
and the maximum air velocity is reached in the region near the upper wall
instead of the air stream core. This could be due to an incorrect modelling of
the drag coefficient. Taking into account results presented in Section 3.3, a 3D simulation was also set up considering a 100×40×20 grid and the “Pierre Coste Large Interface Model” for the drag coefficient.
Unfortunately, only 29 seconds were calculated after one-week run on two
processors, but preliminary results were encouraging.


## 5. Conclusions

A Computational fluid
dynamic analysis of a stratified air-water flow experimental investigation
performed at the Institut de Mécanique des Fluides de Toulouse in 1985 [7]
was performed. The aim was the comprehension of the
experimental data and of the role played by some fundamental parameters. The
simulation has been performed by means of three different CFD codes:
NEPTUNE_CFD V1.0.6, FLUENT 6.1, and CFX 10.0. The spatial discretizations have
been modelled with GAMBIT 1.0 and ANSYS ICEM 10.0 softwares.Preliminary
results of two-phase CFD calculations with a two-dimensional domain suggested
that three-dimensional effects are not negligible, so that 2D simulations are
not suitable to correctly predict this stratified fluid flow. To better
understand the physics of the problem, single-phase analyses were conduced with
FLUENT and ANSYS CFX comparing 2D and 3D simulations for both air and water
single-phase domain. As a result, relevant improvements of both water and air
velocity profile were achieved with 3D simulations. It is worth noting that in
such single-phase analyses, the water level was not
calculated but fixed according to experimental data.Two-phase
simulations by means of NEPTUNE_CFD code, despite taking into consideration a
2D domain, showed better agreement with measured data when considering the new “Pierre
Coste Large Interface Model” for the drag coefficient: water level was
correctly predicted and error in velocity profiles decreased, even if some
underestimation of the air velocity is still present. Moreover, CFX and NEPTUNE_CFD
standard models gave similar results, putting in evidence the fundamental role
played by the drag coefficient modelling. Nevertheless, a systematic
underestimation of the air medium velocity suggests that further information on
the experiment and boundary conditions is needed.



---

*Source: 434212-2008-11-24.xml*