# A Second-Order Turbulence Model Based on a Reynolds Stress Approach for Two-Phase Flow—Part I: Adiabatic Cases

Authors: S. Mimouni; F. Archambeau; M. Boucker; J. Laviéville; C. Morel
Journal: Science and Technology of Nuclear Installations
 (200)

---

## Abstract

In our work in 2008, we evaluated the aptitude of the code Neptune_CFD to reproduce the incidence of a structure topped by vanes on a boiling layer, within the framework of the Neptune project. The objective was to reproduce the main effects of the spacer grids. The turbulence of the liquid phase was modeled by a first-orderK-ε model. We show in this paper that this model is unable to describe the turbulence of
rotating flows, in accordance with the theory. The objective of this paper is to improve the turbulence
modeling of the liquid phase by a second turbulence model based on a Rij-ε approach. Results obtained on typical single-phase cases highlight the improvement of the prediction for all computed values. We tested the turbulence model Rij-ε implemented in the code versus typical adiabatic two-phase flow experiments. We check that the simulations with the Reynolds stress transport model (RSTM) give satisfactory results in a simple geometry as compared to a K-ε model: this point is crucial before calculating rod bundle geometries where the K-ε model may fail.


---

## Body

## 1. Introduction

In
a pressurized water nuclear reactor (PWR), an optimum heat removal from the
surface of the nuclear fuel elements (rod bundle with spacer grids) is very
important for thermal margin and safety.One important goal is to carry out sensitivity
analyses on the angle of the vanes of the fuel assembly spacer grids. In [1], the critical heat flux (CHF) experiment on the effect of the angle and
of the position of mixing vanes was performed in a 2×2 rod bundle. The authors
show that the mixing vanes increase the value of the CHF and the result is
correlated to the magnitude of the swirl generated by the mixing vanes. If the
angle of the mixing vanes is relatively small, the magnitude of the swirling
flow is smaller because the rotating force created by the mixing vanes is weak.
If the angle of the mixing vanes is relatively large, the mixing vanes play the
role of flow obstacle under the departure from nucleate boiling (DNB)
condition. Therefore, it is important that the turbulence modeling deals with
rotation effects.There have been several studies on flow mixing and
heat transfer enhancement caused by a mixing-vane spacer grid in rod bundle
geometry. Lee and Choi [2] simulate the flow field and heat transfer in a
single-phase flow for a 17×17 rod bundle with eight spans of mixing vanes. The
FLUENT commercial code is employed and a Reynolds stress transport model (RSTM)
is used for turbulence. According to the authors, RSTM is helpful. Ikeda et al.
[3] study an assembly consisting of a 5×5 heater rod bundle and eight
specific mixing vane grids. For Ikeda et al., it might be insufficient to apply
a standard K−ε model to
swirl-mixing flow and narrow-channel flow conditions that include nonisotropic
effects. Moreover, In et al. [4] have performed a series of CFD single-phase
flow simulations to analyze the heat transfer enhancement in a fully heated rod
bundle with mixing-vane spacers. For future work, In et al. recommend that a
refined computational fluid dynamic (CFD) model be developed to include details
of the grid structure and a higher-order turbulence model be employed to
improve the accuracy of such simulations.Additional two-phase effects like accumulation of
bubbles in the center of subchannel or pockets of bubbles on the rods should be
taken into account to improve the simulation of flow close to DNB. Indeed,
single-phase simulations remain insufficient and boiling flows simulations are
required. In [5], the authors describe CFD approaches to subcooled
boiling and investigate their capability to contribute to fuel assembly design.
A large part of their work is dedicated to the modeling of boiling flows and to
forces acting on the bubbles. The authors note that the size of bubbles in the
bulk is correlated to the local subcooling which is an important parameter (see
[1]).Considering flow in a PWR core in conditions close to
nominal, when boiling occurs, a high-velocity steady flow takes place with very
small times scales associated to the passage of bubbles (10-4 second −10-3 second) and with quite small
bubble diameters (10-5 m to 10-3 m) compared to the
hydraulic diameter (about 10-2 m). According to the synthesis of the
work performed in WP2.2 of the NURESIM project [6], these are
perfect conditions to use a time average or ensemble average of equations as
usually done in the RANS approach. All turbulent fluctuations and two-phase
intermittency scales can be filtered since they are significantly smaller than
the scales of the mean flow.The large-eddy simulation is also a possible approach.
In the context of the NURESIM project, several studies have been carried out
with a large-eddy simulation to study the axial development of air-water bubbly
flows in a pipe. But in the synthesis of the work performed in WP2.2, Bestion
[6] notes that several open modeling and numerical issues still
remain. So, we will focus on the RANS approach in this paper.According to all these recommendations, a better
understanding of the detailed structure of a flow mixing and heat transfer
downstream of a mixing-vane spacer in a nuclear fuel rod bundle has to be
investigated with an RSTM.Moreover, the overwhelming majority of industrial CFD applications today
are still conducted with two-equation eddy viscosity model, especially the
standardK−ε model, while RSTM remains exceptional.As an example of RSTM development, an RSTM model adapted to bubbly flows
is studied in [7] and used to perform simulations of three basic
bubbly flows (grid, uniform shear, and bubbly wake). The authors decomposed the
Reynolds stress tensor of the liquid into two independent parts: a turbulent
part produced by the mean velocity gradient that also contains the turbulence
of the bubble wakes and a pseudoturbulent part induced by bubble displacements;
each part is predicted from a transport equation. This model is interesting but
has not been selected here for the following reasons. Firstly, the computation
effort is doubled (turbulent and pseudoturbulent parts). Secondly, considering
the flow close to nominal PWR core conditions, when boiling occurs, a high-velocity
steady flow takes place and the bubble diameter is quite small (10-5 m to 10-3 m), therefore the bubbles follow the liquid streamlines
and so the modeling of the pseudoturbulent part induced by bubble displacements
can be omitted. Thirdly, the two-phase flow modeling proposed in [7]
does not tend to a single-phase flow formulation when the void fraction tends
to zero. These three arguments have imposed the choice of the higher-order turbulence
model described in the paper.A second-order moment turbulence model for simulating a bubble column is
also proposed in [8]. The authors defined a Reynolds tensor for each
phase. Furthermore, following a similar method used in deriving and closing the
Reynolds stress equations, a modeled transport equation of two-phase velocity
correlation is also solved. For the same reasons as those mentioned above for
the model proposed in [7], we have not adopted the model proposed in
[8]. The turbulence modeling described in the present paper takes into
account the Reynolds tensor for the liquid only, while a more basic modeling is
used for the vapor phase.However, to the authors’ knowledge, no industrial CFD
approach for boiling flows with an RSTM approach is available in the context of
the fuel assembly design. This may be due to the fact that numerical problems
may occur when using an RSTM approach without caution. Furthermore, the
turbulence modeling of boiling flows is not straightforward. The use of RSTM
also requires finer meshes than eddy
viscosity models (EVMs) and RSTM
may therefore be more time and storage consuming. Developing an industrial RSTM
approach is a quite challenging task but it is worth working at it: RSTM and EVM
results will probably differ and it is important to determine what consequences
it may have on thermal margin and safety of
reactors.In the framework of an R&D program carried out in
the Neptune project (EDF, CEA, AREVA-NP,
IRSN), the following strategy has been adopted:(1)
validation of the
Neptune_CFD code with an RSTM approach on single-phase flow with mixing
vanes and on more academic cases of air-water adiabatic bubbly flows in a pipe;(2)
validation of the
Neptune_CFD code with an RSTM approach on boiling flows in a pipe and
sensitivity to the angle of the vanes for fuel assembly spacer grids
performed in a2×2 rod bundle in a boiling flow configuration;(3)
validation of the
Neptune_CFD code with an RSTM approach for a5×5 rod bundle with mixing
vanes currently used for commercial nuclear fuel.
Step (1) is described in this paper. The second step
should be finalized by the end of 2008. The step (3) is not yet started. Our main objective in this paper is
to check that the simulation with the RSTM gives satisfactory results in a
simple geometry as compared to an EVM: this point is crucial before calculating
rod bundle geometries where the EVM model may fail.This
paper is organized as follows. In Section2.2, the general model we use for adiabatic bubbly flow
simulations is presented in details. In Section 2.3, we underline the weaknesses of the EVM models. In Sections
3.1 and 3.2, the second-moment closure model for high Reynolds
number two-phase flows is presented. In Section 2.4, we give some examples of typical PWR problems that
EVM models fail to represent. In Sections 4.1 and 4.2, respectively, the Liu and Bankoff case and the
sudden expansion experiment are briefly described. The comparison of the
results of Neptune_CFD calculations and the experimental data are presented.
The sensitivity of the numerical results to the turbulence model for the fluid
and to the most important models is studied.Finally,
conclusions are drawn about our current capabilities to simulate bubbly flows
with an RSTM model and perspectives for future work are given.


## 2. The Neptune_CFD Solver and Physical Modeling

### 2.1. Introduction

Neptune_CFD is a three-dimensional two-fluid code
developed more especially for nuclear reactor applications. This local three-dimensional module is based on the
classical two-fluid one pressure approach, including mass, momentum, and energy
balances for each phase.The Neptune_CFD
solver, based on a pressure correction approach, is able to simulate
multicomponent multiphase flows by solving a set of three balance equations for
each field (fluid component and/or phase) [9–12]. These fields can represent many kinds of multiphase flows: distinct
physical components (e.g., gas, liquid, and solid particles), thermodynamic
phases of the same component (e.g., liquid water and its vapor), distinct
physical components, some of which split into different groups (e.g., water and
several groups of different diameter bubbles), and different forms of the same
physical components (e.g., a continuous liquid field, a dispersed liquid field,
a continuous vapor field, and a dispersed vapor field). The solver is based on
a finite volume discretization, together with a collocated arrangement for all
variables. The data structure is totally face-based, which allows the use of
arbitrary shaped cells (tetraedra, hexahedra, prisms, pyramids, etc.) including
nonconforming meshes (meshes with hanging nodes).


### 2.2. Governing Equations and Physical Modeling

The CFD
module of the Neptune software platform is
based on the two-fluid approach [13, 14]. In this
approach, a set of local balance equations for mass, momentum, and energy is
written for each phase. These balance equations are obtained by ensemble
averaging of the local instantaneous balance equations written for the two
phases. When the averaging operation is performed, the major part of the
information about the interfacial configuration and the microphysics governing
the different types of exchanges are lost. As a consequence, a certain number
of closure relations (also called constitutive relations) must be supplied for
the total number of equations (the balance equations and the closure relations)
to be equal to the number of unknown fields. We can distinguish three different
types of closure relations: those which express the interphase exchanges
(interfacial transfer terms), those which express the intraphase exchanges
(molecular and turbulent transfer terms), and those which express the
interactions between each phase and the walls (wall transfer terms). The
balance equations of the two-fluid model we use for adiabatic bubbly flows and
their closure relations are described in the following subsections.


#### 2.2.1. Main Set of Balance Equations

The two-fluid model we use for our adiabatic bubbly
flow calculations consists of the following balance equations.Two mass balance equations(1)∂αkρk∂t+∇⋅(αkρkV¯k)=0,k=l,g, where t is the time, αk,ρk,V¯k denote the time fraction of phase k, its averaged density and velocity.
The phase index k takes the values l for the liquid phase and g for gas bubbles.Two momentum balance equations(2)∂αkρkV¯k∂t+∇⋅(αkρkV¯kV¯k)=−αk∇¯p+M¯k+αkρkg¯+∇⋅[αk(Σ¯¯k+R¯¯k)],k=l,g,
where p is the pressure, g¯ is the
gravity acceleration, M¯k is the interfacial momentum
transfer per unit volume and unit time, and Σ¯¯k and R¯¯k denote the molecular and turbulent stress
tensors, the latter being also called the Reynolds stress tensor.The interfacial transfer of momentumM¯k appearing in the RHS of (2) is assumed to be the sum of four forces: (3)M¯k=M¯kD+M¯kAM+M¯kL+M¯kTD.The
four terms are the averaged drag, added mass, lift, and turbulent dispersion
forces per unit volume. Now we will give the expressions we use for these
forces and for their coefficients.(i)
Drag force(4)M¯gD=−M¯lD=−18AiρlCD|V¯g−V¯l|(V¯g−V¯l),where CD is the drag coefficients
for bubbles which can be determined experimentally.(ii)
Added mass force(5)M¯gAM=−M¯lAM=−CAlg1+2αg1−αgαgρl[(∂V¯g∂t+V¯g⋅∇¯¯V¯g)−(∂V¯l∂t+V¯l⋅∇¯¯V¯l)], where CAlg is the added mass coefficient which is equal
to 1/2 for a spherical bubble and the factor (1+2α)/(1−α) takes into account the effect of the bubbles concentration [15, 16].(iii)
Lift force(6)M¯gL=−M¯lL=−CLαgρl(V¯g−V¯l)⋀(∇¯⋀V¯l), where CL is the lift coefficient.
This coefficient is equal to 1/2 in the particular case of a weakly rotational
flow around a spherical bubble in the limit of infinite Reynolds number [17].(iv)
Turbulent dispersion force(7)M¯gTD=−M¯lTD=−CTDρlKl∇¯αg, where Kl is the liquid turbulent
kinetic energy and CTD is a numerical constant of order 1. This
expression was proposed by Lance and Lopez de Bertodano [18].
An alternative approach is proposed by [19, 20] to model the turbulence induced by bubbles: an algebraic model developed
in the framework of Tchen’s theory, where the turbulent kinetic energy for the
dispersed phase and the covariance are calculated from the turbulent kinetic
energy of the continuous phase.For the dispersed phase, the Reynolds stress tensor is
closed using a Boussinesq-like hypothesis:(8)R¯¯g=ρgνgT(∇¯¯V¯g+∇¯¯TV¯g)−23I¯¯(ρgKg+ρgνgT∇⋅V¯g), where I¯¯ is the identity tensor, with νgT=(1/3)qlg(τlgt−(CAlgτlgF/((ρg/ρl)+CAlg))(τlgF/2))+(1/3)KgτlgF, the turbulent viscosity for the dispersed
phase,  Kg=Kl((b2+ηr)/(1+ηr)), the gas turbulent kinetic energy,  qlg=Kl((b+ηr)/(1+ηr)), the covariance of the dispersed phase,  b=(1+CAlg)/((ρg/ρl)+CAlg),ηr=τlgt/τlgF,  the ratio between the time scale of the
continuous phase turbulence viewed by the dispersed phase (takes into account
crossing trajectories effect) and the characteristic time scale of the momentum
transfer rate between the liquid and dispersed phases: (9)τlgt=τltσα(1+Cβξr2)−1/2,τlgF=CAlg+ρgFDlg, with σα is the turbulent Schmidt or Prandtl turbulent
for the continuous phase, Cβ is the crossing trajectories coefficient taken
equal to 1.8, and CAlg is the added mass coefficient; (10)ξr=〈|Vr→|〉g(1/3)Kl,τlt=32CμKlεl. Turbulent dispersion force and Tchen’s model will be
compared below for bubbly flows in a straight pipe and in a sudden expansion.


#### 2.2.2. Turbulent Transfer Terms

TheK−ε model describes energy
processes in terms of production and dissipation, as well as transport through
the mean flow or by turbulent diffusion. The Kolmogorov spectral equilibrium
hypothesis also enables one to predict a large eddy length-scale. On the other
hand, the anisotropy of the stresses is quite crudely modeled. First of all the
EVM model assumes the Reynolds stress tensor is aligned with the strain rate
tensor (Boussinesq approximation): (11)R¯¯l=ρlνlT(∇¯¯V¯l+∇¯¯TV¯l)−23I¯¯(ρlKl+ρlνlT∇⋅V¯l), where I¯¯ is the identity tensor, Kl is the liquid turbulent kinetic energy, and νlT is the liquid turbulent eddy viscosity. The
liquid turbulent eddy viscosity is expressed by the following relation: (12)νlT=CμKl2εl, where Cμ=0.09. The turbulent kinetic energy Kl and its dissipation
rate εl are calculated by using the
two-equation K−ε approach.


### 2.3. EVM Weaknesses: Theoretical Approach

Flows
encountered in vertical pipe are of great interest to validate the most
important heat, mass, and momentum closure relations. However, some negligible
effects in simple geometry sometimes become preponderant in complex geometries.
For example, the modeling of two-phase flow in water-cooled nuclear reactors
needs to take into account swirls and stagnation points. Applications in
complex geometries also need to take into account the complex features of the
secondary motions which are observed experimentally. These requirements
highlight the need for meticulous turbulence modeling.A reason for the persistent widespread use of low-level turbulence
modeling in two-phase CFD is perhaps the fact that the use of two-phase CFD in
complex industrial geometries is only starting. Moreover, many studies merely
require “order of magnitude” or “good tendencies” answers.However, extensive testing and application over the past three decades
have revealed a number of shortcomings and deficiencies in EVM models, and
among them theK−ε model, such as(i)
limitation to linear algebraic stress-strain relationship (poor
performances wherever the stress transport is important, e.g., non equilibrium,
fast evolving, separating, and buoyant flows),(ii)
insensitivity to the orientation of turbulence structure and stress
anisotropy (poor performances where normal stresses play an important role,
e.g., stress-driven secondary flows in noncircular ducts),(iii)
inability to account for extra strain (streamline curvature, skewing,
rotation),(iv)
poor prediction particularly of flows with strong adverse pressure
gradients and in reattachment regions.
In a plane strain situation, such as upstream of a
stagnation point on a bluff body, the exact (as obtained by an RTSM) production
and that obtained from an EVM are, respectively, [21] (13)Pexact=−(ux′2¯−uy′2¯)∂Vl,x∂x,PEVM=4νl,t(∂Vl,x∂x)2. The difference between the normal stresses actually
grows slowly on the short-time scale needed for the flow to travel around the
stagnation point, so the production remains moderate, and in any case is
bounded, whereas PEVM usually yields a severe over-prediction when
the strain is high.The simulation of swirling flow generated by the
mixing vanes is our main goal since it plays an important role for the
prediction of the CHF for the fuel assemblies. For this reason, the rotation
effects are more specifically addressed hereafter.It can be easily shown [22] that in the
presence of an initially anisotropy turbulence, rotation will cause a
redistribution of energy between normal components without affecting the value
of this quantity. In fact, the angular velocity Θ does not appear explicitly in the K-equation, obtained
by adding the normal stresses: (14)dKdt=−ε. Thus, the K−ε model is
totally blind to rotation effects. The swirling flows can be regarded as a
special case of fluid rotation with the axis usually aligned with the mean flow
direction so that the Coriolis force is zero. This aspect is crucial for the
simulation of hot channel of a fuel assembly. In fact, mixing vanes at the
spacer grids generate a swirl in the coolant water to enhance the heat transfer
from the rods to the coolants in the hot channels and to limit boiling.In the following section, we present some examples of
large-scale industrial applications, performed using eddy-viscosity models, and
subsequently discuss areas of weakness of the models, highlighting some
improvements that can be obtained through the use of more advances stress
transport closures.


### 2.4. EVM Weaknesses: Illustration on the AGATE-Mixing and DEBORA-Mixing Experiment

Keeping in mind the long-term objective (two-phase CFD
calculations validated under typical pressurized water reactor (PWR) geometries
and thermalhydraulic conditions), we started very recently to evaluate
Neptune_CFD against spacer grid type experiments. An experimental device
representing three mixing blades (Figure2) was introduced in a heated tube (diameter =
19.2 mm) and used for two different programs.(i)AGATE-mixing experiment [23]. Single-phase
liquid water tests, with laser-doppler liquid velocity measurements
upstream and downstream the mixing blades (for each of the 15 horizontal
planes, the liquid velocity is measured along 12 different diameters and
there are 12 points for each radius); the velocity at inlet is 3 m/s and
the pressure is 2 bar.(ii)DEBORA-mixing experiment [24]. Boiling R-12
freon tests on the same geometry
but the total length of the calculation domain is 3.5 m;
the tube is heated and the
uniform wall heat flux is 109300 W/m2 which gives about 2% of
vapor at outlet; the outlet pressure is 26.2 bar; the inlet liquid
temperature is 63.3°C; the inlet liquid mass is 0.873 kg/s. The main physical phenomena to reproduce
are wall boiling, entrainment of bubbles in the wakes, and recondensation.
So, the prediction of the swirl is crucial.For the mixing blades part (60 mm),
77000 cells are needed. This grid is considered as a reasonable compromise
between the numerical accuracy and the computational effort. Figure1 compares computed and experimental orthoradial
(circular component in a horizontal plane) liquid velocity downstream the
mixing vane (AGATE test). One can notice that the rotating flow is
qualitatively well reproduced by Neptune_CFD although the velocity is
underestimated. This is mainly due to the turbulence model (standard K−ε here) which is not optimum for this type of
geometry.Figure 1
Orthoradial liquid velocity downstream the mixing vane (Agate-mixing experiment).Figure 2
View of the mixing device.TheK−ε model
underestimates orthoradial velocities downstream the blades, but the results
remain qualitatively satisfactory. The Rij−ε model gives
satisfactory results (Figure 1).In the following section, we propose a second-moment closure model to
take into account the liquid turbulence in order to validate in the long-term
calculations in typical pressurized water reactor (PWR) geometries and
thermalhydraulic conditions.In
the present paper, we suppose that RSTM is well known in single-phase flow
[21]. Now, our objective is to test and if possible to improve our
RSTM model adapted to bubbly flows as compared to
experimental data and K−ε results. Indeed, we are interested by two-phase high
Reynolds numbers flows, but beforehand, the mechanical models implemented in
the Neptune_CFD code must be tested on the simpler cases of air-water adiabatic
bubbly flows.


## 2.1. Introduction

Neptune_CFD is a three-dimensional two-fluid code
developed more especially for nuclear reactor applications. This local three-dimensional module is based on the
classical two-fluid one pressure approach, including mass, momentum, and energy
balances for each phase.The Neptune_CFD
solver, based on a pressure correction approach, is able to simulate
multicomponent multiphase flows by solving a set of three balance equations for
each field (fluid component and/or phase) [9–12]. These fields can represent many kinds of multiphase flows: distinct
physical components (e.g., gas, liquid, and solid particles), thermodynamic
phases of the same component (e.g., liquid water and its vapor), distinct
physical components, some of which split into different groups (e.g., water and
several groups of different diameter bubbles), and different forms of the same
physical components (e.g., a continuous liquid field, a dispersed liquid field,
a continuous vapor field, and a dispersed vapor field). The solver is based on
a finite volume discretization, together with a collocated arrangement for all
variables. The data structure is totally face-based, which allows the use of
arbitrary shaped cells (tetraedra, hexahedra, prisms, pyramids, etc.) including
nonconforming meshes (meshes with hanging nodes).


## 2.2. Governing Equations and Physical Modeling

The CFD
module of the Neptune software platform is
based on the two-fluid approach [13, 14]. In this
approach, a set of local balance equations for mass, momentum, and energy is
written for each phase. These balance equations are obtained by ensemble
averaging of the local instantaneous balance equations written for the two
phases. When the averaging operation is performed, the major part of the
information about the interfacial configuration and the microphysics governing
the different types of exchanges are lost. As a consequence, a certain number
of closure relations (also called constitutive relations) must be supplied for
the total number of equations (the balance equations and the closure relations)
to be equal to the number of unknown fields. We can distinguish three different
types of closure relations: those which express the interphase exchanges
(interfacial transfer terms), those which express the intraphase exchanges
(molecular and turbulent transfer terms), and those which express the
interactions between each phase and the walls (wall transfer terms). The
balance equations of the two-fluid model we use for adiabatic bubbly flows and
their closure relations are described in the following subsections.


### 2.2.1. Main Set of Balance Equations

The two-fluid model we use for our adiabatic bubbly
flow calculations consists of the following balance equations.Two mass balance equations(1)∂αkρk∂t+∇⋅(αkρkV¯k)=0,k=l,g, where t is the time, αk,ρk,V¯k denote the time fraction of phase k, its averaged density and velocity.
The phase index k takes the values l for the liquid phase and g for gas bubbles.Two momentum balance equations(2)∂αkρkV¯k∂t+∇⋅(αkρkV¯kV¯k)=−αk∇¯p+M¯k+αkρkg¯+∇⋅[αk(Σ¯¯k+R¯¯k)],k=l,g,
where p is the pressure, g¯ is the
gravity acceleration, M¯k is the interfacial momentum
transfer per unit volume and unit time, and Σ¯¯k and R¯¯k denote the molecular and turbulent stress
tensors, the latter being also called the Reynolds stress tensor.The interfacial transfer of momentumM¯k appearing in the RHS of (2) is assumed to be the sum of four forces: (3)M¯k=M¯kD+M¯kAM+M¯kL+M¯kTD.The
four terms are the averaged drag, added mass, lift, and turbulent dispersion
forces per unit volume. Now we will give the expressions we use for these
forces and for their coefficients.(i)
Drag force(4)M¯gD=−M¯lD=−18AiρlCD|V¯g−V¯l|(V¯g−V¯l),where CD is the drag coefficients
for bubbles which can be determined experimentally.(ii)
Added mass force(5)M¯gAM=−M¯lAM=−CAlg1+2αg1−αgαgρl[(∂V¯g∂t+V¯g⋅∇¯¯V¯g)−(∂V¯l∂t+V¯l⋅∇¯¯V¯l)], where CAlg is the added mass coefficient which is equal
to 1/2 for a spherical bubble and the factor (1+2α)/(1−α) takes into account the effect of the bubbles concentration [15, 16].(iii)
Lift force(6)M¯gL=−M¯lL=−CLαgρl(V¯g−V¯l)⋀(∇¯⋀V¯l), where CL is the lift coefficient.
This coefficient is equal to 1/2 in the particular case of a weakly rotational
flow around a spherical bubble in the limit of infinite Reynolds number [17].(iv)
Turbulent dispersion force(7)M¯gTD=−M¯lTD=−CTDρlKl∇¯αg, where Kl is the liquid turbulent
kinetic energy and CTD is a numerical constant of order 1. This
expression was proposed by Lance and Lopez de Bertodano [18].
An alternative approach is proposed by [19, 20] to model the turbulence induced by bubbles: an algebraic model developed
in the framework of Tchen’s theory, where the turbulent kinetic energy for the
dispersed phase and the covariance are calculated from the turbulent kinetic
energy of the continuous phase.For the dispersed phase, the Reynolds stress tensor is
closed using a Boussinesq-like hypothesis:(8)R¯¯g=ρgνgT(∇¯¯V¯g+∇¯¯TV¯g)−23I¯¯(ρgKg+ρgνgT∇⋅V¯g), where I¯¯ is the identity tensor, with νgT=(1/3)qlg(τlgt−(CAlgτlgF/((ρg/ρl)+CAlg))(τlgF/2))+(1/3)KgτlgF, the turbulent viscosity for the dispersed
phase,  Kg=Kl((b2+ηr)/(1+ηr)), the gas turbulent kinetic energy,  qlg=Kl((b+ηr)/(1+ηr)), the covariance of the dispersed phase,  b=(1+CAlg)/((ρg/ρl)+CAlg),ηr=τlgt/τlgF,  the ratio between the time scale of the
continuous phase turbulence viewed by the dispersed phase (takes into account
crossing trajectories effect) and the characteristic time scale of the momentum
transfer rate between the liquid and dispersed phases: (9)τlgt=τltσα(1+Cβξr2)−1/2,τlgF=CAlg+ρgFDlg, with σα is the turbulent Schmidt or Prandtl turbulent
for the continuous phase, Cβ is the crossing trajectories coefficient taken
equal to 1.8, and CAlg is the added mass coefficient; (10)ξr=〈|Vr→|〉g(1/3)Kl,τlt=32CμKlεl. Turbulent dispersion force and Tchen’s model will be
compared below for bubbly flows in a straight pipe and in a sudden expansion.


### 2.2.2. Turbulent Transfer Terms

TheK−ε model describes energy
processes in terms of production and dissipation, as well as transport through
the mean flow or by turbulent diffusion. The Kolmogorov spectral equilibrium
hypothesis also enables one to predict a large eddy length-scale. On the other
hand, the anisotropy of the stresses is quite crudely modeled. First of all the
EVM model assumes the Reynolds stress tensor is aligned with the strain rate
tensor (Boussinesq approximation): (11)R¯¯l=ρlνlT(∇¯¯V¯l+∇¯¯TV¯l)−23I¯¯(ρlKl+ρlνlT∇⋅V¯l), where I¯¯ is the identity tensor, Kl is the liquid turbulent kinetic energy, and νlT is the liquid turbulent eddy viscosity. The
liquid turbulent eddy viscosity is expressed by the following relation: (12)νlT=CμKl2εl, where Cμ=0.09. The turbulent kinetic energy Kl and its dissipation
rate εl are calculated by using the
two-equation K−ε approach.


## 2.2.1. Main Set of Balance Equations

The two-fluid model we use for our adiabatic bubbly
flow calculations consists of the following balance equations.Two mass balance equations(1)∂αkρk∂t+∇⋅(αkρkV¯k)=0,k=l,g, where t is the time, αk,ρk,V¯k denote the time fraction of phase k, its averaged density and velocity.
The phase index k takes the values l for the liquid phase and g for gas bubbles.Two momentum balance equations(2)∂αkρkV¯k∂t+∇⋅(αkρkV¯kV¯k)=−αk∇¯p+M¯k+αkρkg¯+∇⋅[αk(Σ¯¯k+R¯¯k)],k=l,g,
where p is the pressure, g¯ is the
gravity acceleration, M¯k is the interfacial momentum
transfer per unit volume and unit time, and Σ¯¯k and R¯¯k denote the molecular and turbulent stress
tensors, the latter being also called the Reynolds stress tensor.The interfacial transfer of momentumM¯k appearing in the RHS of (2) is assumed to be the sum of four forces: (3)M¯k=M¯kD+M¯kAM+M¯kL+M¯kTD.The
four terms are the averaged drag, added mass, lift, and turbulent dispersion
forces per unit volume. Now we will give the expressions we use for these
forces and for their coefficients.(i)
Drag force(4)M¯gD=−M¯lD=−18AiρlCD|V¯g−V¯l|(V¯g−V¯l),where CD is the drag coefficients
for bubbles which can be determined experimentally.(ii)
Added mass force(5)M¯gAM=−M¯lAM=−CAlg1+2αg1−αgαgρl[(∂V¯g∂t+V¯g⋅∇¯¯V¯g)−(∂V¯l∂t+V¯l⋅∇¯¯V¯l)], where CAlg is the added mass coefficient which is equal
to 1/2 for a spherical bubble and the factor (1+2α)/(1−α) takes into account the effect of the bubbles concentration [15, 16].(iii)
Lift force(6)M¯gL=−M¯lL=−CLαgρl(V¯g−V¯l)⋀(∇¯⋀V¯l), where CL is the lift coefficient.
This coefficient is equal to 1/2 in the particular case of a weakly rotational
flow around a spherical bubble in the limit of infinite Reynolds number [17].(iv)
Turbulent dispersion force(7)M¯gTD=−M¯lTD=−CTDρlKl∇¯αg, where Kl is the liquid turbulent
kinetic energy and CTD is a numerical constant of order 1. This
expression was proposed by Lance and Lopez de Bertodano [18].
An alternative approach is proposed by [19, 20] to model the turbulence induced by bubbles: an algebraic model developed
in the framework of Tchen’s theory, where the turbulent kinetic energy for the
dispersed phase and the covariance are calculated from the turbulent kinetic
energy of the continuous phase.For the dispersed phase, the Reynolds stress tensor is
closed using a Boussinesq-like hypothesis:(8)R¯¯g=ρgνgT(∇¯¯V¯g+∇¯¯TV¯g)−23I¯¯(ρgKg+ρgνgT∇⋅V¯g), where I¯¯ is the identity tensor, with νgT=(1/3)qlg(τlgt−(CAlgτlgF/((ρg/ρl)+CAlg))(τlgF/2))+(1/3)KgτlgF, the turbulent viscosity for the dispersed
phase,  Kg=Kl((b2+ηr)/(1+ηr)), the gas turbulent kinetic energy,  qlg=Kl((b+ηr)/(1+ηr)), the covariance of the dispersed phase,  b=(1+CAlg)/((ρg/ρl)+CAlg),ηr=τlgt/τlgF,  the ratio between the time scale of the
continuous phase turbulence viewed by the dispersed phase (takes into account
crossing trajectories effect) and the characteristic time scale of the momentum
transfer rate between the liquid and dispersed phases: (9)τlgt=τltσα(1+Cβξr2)−1/2,τlgF=CAlg+ρgFDlg, with σα is the turbulent Schmidt or Prandtl turbulent
for the continuous phase, Cβ is the crossing trajectories coefficient taken
equal to 1.8, and CAlg is the added mass coefficient; (10)ξr=〈|Vr→|〉g(1/3)Kl,τlt=32CμKlεl. Turbulent dispersion force and Tchen’s model will be
compared below for bubbly flows in a straight pipe and in a sudden expansion.


## 2.2.2. Turbulent Transfer Terms

TheK−ε model describes energy
processes in terms of production and dissipation, as well as transport through
the mean flow or by turbulent diffusion. The Kolmogorov spectral equilibrium
hypothesis also enables one to predict a large eddy length-scale. On the other
hand, the anisotropy of the stresses is quite crudely modeled. First of all the
EVM model assumes the Reynolds stress tensor is aligned with the strain rate
tensor (Boussinesq approximation): (11)R¯¯l=ρlνlT(∇¯¯V¯l+∇¯¯TV¯l)−23I¯¯(ρlKl+ρlνlT∇⋅V¯l), where I¯¯ is the identity tensor, Kl is the liquid turbulent kinetic energy, and νlT is the liquid turbulent eddy viscosity. The
liquid turbulent eddy viscosity is expressed by the following relation: (12)νlT=CμKl2εl, where Cμ=0.09. The turbulent kinetic energy Kl and its dissipation
rate εl are calculated by using the
two-equation K−ε approach.


## 2.3. EVM Weaknesses: Theoretical Approach

Flows
encountered in vertical pipe are of great interest to validate the most
important heat, mass, and momentum closure relations. However, some negligible
effects in simple geometry sometimes become preponderant in complex geometries.
For example, the modeling of two-phase flow in water-cooled nuclear reactors
needs to take into account swirls and stagnation points. Applications in
complex geometries also need to take into account the complex features of the
secondary motions which are observed experimentally. These requirements
highlight the need for meticulous turbulence modeling.A reason for the persistent widespread use of low-level turbulence
modeling in two-phase CFD is perhaps the fact that the use of two-phase CFD in
complex industrial geometries is only starting. Moreover, many studies merely
require “order of magnitude” or “good tendencies” answers.However, extensive testing and application over the past three decades
have revealed a number of shortcomings and deficiencies in EVM models, and
among them theK−ε model, such as(i)
limitation to linear algebraic stress-strain relationship (poor
performances wherever the stress transport is important, e.g., non equilibrium,
fast evolving, separating, and buoyant flows),(ii)
insensitivity to the orientation of turbulence structure and stress
anisotropy (poor performances where normal stresses play an important role,
e.g., stress-driven secondary flows in noncircular ducts),(iii)
inability to account for extra strain (streamline curvature, skewing,
rotation),(iv)
poor prediction particularly of flows with strong adverse pressure
gradients and in reattachment regions.
In a plane strain situation, such as upstream of a
stagnation point on a bluff body, the exact (as obtained by an RTSM) production
and that obtained from an EVM are, respectively, [21] (13)Pexact=−(ux′2¯−uy′2¯)∂Vl,x∂x,PEVM=4νl,t(∂Vl,x∂x)2. The difference between the normal stresses actually
grows slowly on the short-time scale needed for the flow to travel around the
stagnation point, so the production remains moderate, and in any case is
bounded, whereas PEVM usually yields a severe over-prediction when
the strain is high.The simulation of swirling flow generated by the
mixing vanes is our main goal since it plays an important role for the
prediction of the CHF for the fuel assemblies. For this reason, the rotation
effects are more specifically addressed hereafter.It can be easily shown [22] that in the
presence of an initially anisotropy turbulence, rotation will cause a
redistribution of energy between normal components without affecting the value
of this quantity. In fact, the angular velocity Θ does not appear explicitly in the K-equation, obtained
by adding the normal stresses: (14)dKdt=−ε. Thus, the K−ε model is
totally blind to rotation effects. The swirling flows can be regarded as a
special case of fluid rotation with the axis usually aligned with the mean flow
direction so that the Coriolis force is zero. This aspect is crucial for the
simulation of hot channel of a fuel assembly. In fact, mixing vanes at the
spacer grids generate a swirl in the coolant water to enhance the heat transfer
from the rods to the coolants in the hot channels and to limit boiling.In the following section, we present some examples of
large-scale industrial applications, performed using eddy-viscosity models, and
subsequently discuss areas of weakness of the models, highlighting some
improvements that can be obtained through the use of more advances stress
transport closures.


## 2.4. EVM Weaknesses: Illustration on the AGATE-Mixing and DEBORA-Mixing Experiment

Keeping in mind the long-term objective (two-phase CFD
calculations validated under typical pressurized water reactor (PWR) geometries
and thermalhydraulic conditions), we started very recently to evaluate
Neptune_CFD against spacer grid type experiments. An experimental device
representing three mixing blades (Figure2) was introduced in a heated tube (diameter =
19.2 mm) and used for two different programs.(i)AGATE-mixing experiment [23]. Single-phase
liquid water tests, with laser-doppler liquid velocity measurements
upstream and downstream the mixing blades (for each of the 15 horizontal
planes, the liquid velocity is measured along 12 different diameters and
there are 12 points for each radius); the velocity at inlet is 3 m/s and
the pressure is 2 bar.(ii)DEBORA-mixing experiment [24]. Boiling R-12
freon tests on the same geometry
but the total length of the calculation domain is 3.5 m;
the tube is heated and the
uniform wall heat flux is 109300 W/m2 which gives about 2% of
vapor at outlet; the outlet pressure is 26.2 bar; the inlet liquid
temperature is 63.3°C; the inlet liquid mass is 0.873 kg/s. The main physical phenomena to reproduce
are wall boiling, entrainment of bubbles in the wakes, and recondensation.
So, the prediction of the swirl is crucial.For the mixing blades part (60 mm),
77000 cells are needed. This grid is considered as a reasonable compromise
between the numerical accuracy and the computational effort. Figure1 compares computed and experimental orthoradial
(circular component in a horizontal plane) liquid velocity downstream the
mixing vane (AGATE test). One can notice that the rotating flow is
qualitatively well reproduced by Neptune_CFD although the velocity is
underestimated. This is mainly due to the turbulence model (standard K−ε here) which is not optimum for this type of
geometry.Figure 1
Orthoradial liquid velocity downstream the mixing vane (Agate-mixing experiment).Figure 2
View of the mixing device.TheK−ε model
underestimates orthoradial velocities downstream the blades, but the results
remain qualitatively satisfactory. The Rij−ε model gives
satisfactory results (Figure 1).In the following section, we propose a second-moment closure model to
take into account the liquid turbulence in order to validate in the long-term
calculations in typical pressurized water reactor (PWR) geometries and
thermalhydraulic conditions.In
the present paper, we suppose that RSTM is well known in single-phase flow
[21]. Now, our objective is to test and if possible to improve our
RSTM model adapted to bubbly flows as compared to
experimental data and K−ε results. Indeed, we are interested by two-phase high
Reynolds numbers flows, but beforehand, the mechanical models implemented in
the Neptune_CFD code must be tested on the simpler cases of air-water adiabatic
bubbly flows.


## 3. The Second-Moment Closure Model for High Reynolds Number Flows Dedicated to the Continuous Phase (Liquid)

In this section, we omit the subscript “l” for the
liquid and “α” is the void fraction for the
sake of simplicity.


### 3.1. Equation onRij

We have(15)(1−α)DρRijDt=∂∂xk{(ρν+ρCsKεRij)∂∂xk((1−α)⋅Rij)}+(1−α)(Pij+Gij+Φij+εij). In this model, the Reynolds stress tensor of the
continuous phase is split into two parts, a turbulent dissipative part produced
by the gradient of mean velocity and by the wakes of the bubbles and a
pseudoturbulent nondissipative part induced by the displacements of the
bubbles. The displacements of the bubbles should be taken into account in
experiments, where air is injected at the bottom of a water pool creating a
large, axisymmetric bubble plume with a large-scale recirculation flow around
the plume. But swirling flows and high Reynolds number characterize our
industrial applications.Hence, we neglected, in our approach and in first
analysis, the nondissipative component called “pseudoturbulent.” We
consider only the “turbulent” dissipative part. Within this
framework, the term of production by the bubbles interfaces is written as
[7] (16)−(pρui′nj+pρuj′ni)δI+ν(∂∂xkui′uj′)nkδI¯, where n indicates the normal to the interface
and δI a Dirac
function on the interface. It was omitted in [7]. Indeed, according
to [7], dissipation in the wakes is balanced by the interfacial
production: the equation of transport of the Reynolds stress tensor has the
same form as in the single-phase case and is given by (15). When the void
fraction is vanishing, the two-phase flow modeling naturally degenerates to the
single-phase flow modeling.Some terms of the equation of transport of the
Reynolds stress tensor cannot be computed directly and must be modeled. A
modeling resulting from [21] is proposed below.A common way to model the viscous destruction of
stresses for high Reynolds number flows is(17)εij=23εδij. The
turbulent diffusion is of diffusive nature and the most popular model is the
generalized gradient diffusion: (18)Dijt=∂∂xk(CsKεuk′ul′¯∂ui′uj′¯∂xl). Pressure fluctuations tend to
disrupt the turbulent structures and to redistribute the energy to make turbulence
more isotropic: Φij=Φij,1+Φij,2+Φij,3+Φij,1ω+Φij,2ω+Φij,3ω with(i)(19)Φij,1=−C1ε(ui′uj′¯K−23δij)withC1=1.8,(ii)(20)Φij,2=−C2(Pij−23Pδij),withP=12Pkk,C2=0.6,Φij,3=−C3(Gij−23Gδij),withG=12Gkk,C3=0.55,(iii)(21)Φij,1ω=C1ωεK(uk′um′¯nk¯nm¯δij−32ui′uk′¯nk¯nj¯−32uk′uj′¯nk¯ni¯)⋅fω,(iv)(22)Φij,2ω=C2ω(Φkm,2nk¯nm¯δij−32Φik,2nk¯nj¯−32Φjk,2nk¯ni¯)⋅fω,withC1ω=0.5,C2ω=0.3,fω=0.4⋅K3/2/ε⋅xn, where xn is the distance to the
wall and nk¯ the base vector normal to the wall.(v)G is the production by body force.


### 3.2. Equation onε

In the RTSM closures, the same basic form of model
equation forε is used as in the K−ε model,
except that now (uk′ul′¯) is available, which has the following
implications.The production of kinetic energy (P and G) in the source term of ε is treated in exact form.The generalized gradient
hypothesis is used to model turbulent diffusion.Hence, the model equation forε has the form (23)(1−α)DεDt=∂∂xk(CεKεuk′ul′¯∂(1−α)ε∂xl)+(Cε1P+Cε3G+Cε4K∂Vk∂xk−Cε2ε)⋅(1−α)εK. The coefficients of the Rij−ε  model are shown in Table 2.


## 3.1. Equation onRij

We have(15)(1−α)DρRijDt=∂∂xk{(ρν+ρCsKεRij)∂∂xk((1−α)⋅Rij)}+(1−α)(Pij+Gij+Φij+εij). In this model, the Reynolds stress tensor of the
continuous phase is split into two parts, a turbulent dissipative part produced
by the gradient of mean velocity and by the wakes of the bubbles and a
pseudoturbulent nondissipative part induced by the displacements of the
bubbles. The displacements of the bubbles should be taken into account in
experiments, where air is injected at the bottom of a water pool creating a
large, axisymmetric bubble plume with a large-scale recirculation flow around
the plume. But swirling flows and high Reynolds number characterize our
industrial applications.Hence, we neglected, in our approach and in first
analysis, the nondissipative component called “pseudoturbulent.” We
consider only the “turbulent” dissipative part. Within this
framework, the term of production by the bubbles interfaces is written as
[7] (16)−(pρui′nj+pρuj′ni)δI+ν(∂∂xkui′uj′)nkδI¯, where n indicates the normal to the interface
and δI a Dirac
function on the interface. It was omitted in [7]. Indeed, according
to [7], dissipation in the wakes is balanced by the interfacial
production: the equation of transport of the Reynolds stress tensor has the
same form as in the single-phase case and is given by (15). When the void
fraction is vanishing, the two-phase flow modeling naturally degenerates to the
single-phase flow modeling.Some terms of the equation of transport of the
Reynolds stress tensor cannot be computed directly and must be modeled. A
modeling resulting from [21] is proposed below.A common way to model the viscous destruction of
stresses for high Reynolds number flows is(17)εij=23εδij. The
turbulent diffusion is of diffusive nature and the most popular model is the
generalized gradient diffusion: (18)Dijt=∂∂xk(CsKεuk′ul′¯∂ui′uj′¯∂xl). Pressure fluctuations tend to
disrupt the turbulent structures and to redistribute the energy to make turbulence
more isotropic: Φij=Φij,1+Φij,2+Φij,3+Φij,1ω+Φij,2ω+Φij,3ω with(i)(19)Φij,1=−C1ε(ui′uj′¯K−23δij)withC1=1.8,(ii)(20)Φij,2=−C2(Pij−23Pδij),withP=12Pkk,C2=0.6,Φij,3=−C3(Gij−23Gδij),withG=12Gkk,C3=0.55,(iii)(21)Φij,1ω=C1ωεK(uk′um′¯nk¯nm¯δij−32ui′uk′¯nk¯nj¯−32uk′uj′¯nk¯ni¯)⋅fω,(iv)(22)Φij,2ω=C2ω(Φkm,2nk¯nm¯δij−32Φik,2nk¯nj¯−32Φjk,2nk¯ni¯)⋅fω,withC1ω=0.5,C2ω=0.3,fω=0.4⋅K3/2/ε⋅xn, where xn is the distance to the
wall and nk¯ the base vector normal to the wall.(v)G is the production by body force.


## 3.2. Equation onε

In the RTSM closures, the same basic form of model
equation forε is used as in the K−ε model,
except that now (uk′ul′¯) is available, which has the following
implications.The production of kinetic energy (P and G) in the source term of ε is treated in exact form.The generalized gradient
hypothesis is used to model turbulent diffusion.Hence, the model equation forε has the form (23)(1−α)DεDt=∂∂xk(CεKεuk′ul′¯∂(1−α)ε∂xl)+(Cε1P+Cε3G+Cε4K∂Vk∂xk−Cε2ε)⋅(1−α)εK. The coefficients of the Rij−ε  model are shown in Table 2.


## 4. Validation on Adiabatic Bubbly Flow Cases

### 4.1. The LIU&BANKOFF Case [25]

In this section, we evaluate the modeling capabilities
to simulate an upward bubbly flow in adiabatic conditions. We do not
specifically optimize the coefficients and modeling of the momentum transfer
terms to get results as good as possible. Our main objective is to validate theRij−ε turbulence model for the fluid against the K−ε one.The test section was aZ=2800 mm long, vertical
smooth acrylic tubing, with inside diameter of D=38 mm.The set of physical properties is the following:(24)ρliquid=994.9kg⋅m−3,ρgas=1.6kg⋅m−3,μliquid=7.97⋅10−4N⋅s⋅m−2,μgas=1.748⋅10−5N⋅s⋅m−2,P0=101500Pa,g=9.81m⋅s−2,dbubble=2.5⋅10−3m.  A uniform axial liquid profile is imposed at the inlet
and is equal to 1.138 m/s. A uniform axial gas profile is imposed at the inlet
and is equal to 1.333 m/s. The void fraction at the inlet is 0.045.The interfacial momentum transfer term is assumed to
be the sum of four different forces. The three first ones are simplified
averaged expressions of the classical drag, added mass, and lift force. The
fourth is the turbulent dispersion force.The flow is assumed to be axisymmetric therefore a
two-dimensional axisymmetric meshing is used. Computations have been performed
on two kinds of meshing: a coarse grid (20 cells in the radial direction and 50
cells in the axial direction) and a fine grid (30 cells in the radial direction
and 100 cells in the axial direction). Results are similar and computations are
performed on the first grid.At the measuring station (Z/D=36), we compare numerical results against
experimental data for the axial liquid velocity and the void fraction.As
recommended in [26], the lift coefficient is taken to be equal to
0.1. We have also tested the turbulent dispersion force of Davidson model
[26] written as (25)FDT=−34CDdbμt(1−α)|Vg−Vl|∂α∂xiwithμt=0.09ρlKl2ε. This expression gives values
negligible with respect to the Lopez de Bertodano expression [18] (the ratio is about 1000) and the void fraction profile is
observed to be similar to the profile calculated without any turbulence model
for the dispersed phase. If we take the bubble fluctuation into account with
the Hinze-Tchen algebraic model of bubble turbulence, we get the same results
with the turbulent dispersion force.In [27], Grossetête
considers that bubbles are deformed near the wall. To take into account this
effect, the author proposes to put a negative lift coefficient near the wall
and otherwise a positive one. Nevertheless, calculations (not presented in this
paper) show a peak of void fraction at the wall. We have also tested the
Tomiyama lift force [28, 29], but results are not improved.
Moreover, a wall lubrification force [30] can push the bubbles away
from the wall and improve the results.Figures3 and 4 show that the simulation results are in quite good
agreement with the experimental data. In [31], computations performed
with a K−ε turbulence model for the liquid produced comparable results. Our main
objective in this paper is to check that the simulation with the Rij−ε
turbulence model gives satisfactory results in a simple geometry, which is
crucial before calculating industrial geometries, where the K−ε turbulence model is susceptible to fail.Figure 3
Radial profile of the liquid velocity, comparisons of turbulent dispersion forces withCL=0.1.Figure 4
Radial profile of void fraction, comparisons of turbulent dispersion forces withCL=0.1.But improvement of the modeling of the interfacial
forces exerted on bubbles by the surrounding liquid is required. A strong
sensitivity to the lift coefficient has been found in our calculations. Other
forces, like the turbulent dispersion force, have also a crucial effect on the
void fraction distribution. These forces depend on uncontrolled parameters like
the bubble shape, the liquid turbulence, the bubbles collective effects, and so
on [31].


### 4.2. Sudden Expansion Experiment

Bel Fdhila [32] investigated
experimentally several upwards bubbly flows in a vertical pipe with a sudden
expansion. The total length of the pipe was equal to 14 meters. The bottom part
of the tube had a length equal to 9 meters and an internal diameter equal to 50 mm, the top part of the tube (5 m length) having an internal diameter equal to
100 mm. The fluids used were water and air under atmospheric pressure and
ambient temperature. Six measuring sections were located upstream and
downstream of the singularity. The first measuring section was located two
centimetres before the singularity, the other five were located at 7, 13, 18,
25, and 32 cm above the singularity. In each measuring section, the radial
profile of the void fraction was measured by means of a single optical probe,
and two components of the liquid velocity were measured by means of a hot film
anemometer. The time-averaged components of this velocity field and three
components of the liquid Reynolds stress tensor were deduced (the flow being
assumed axisymmetric). The bubble size was not measured in the experiment.
According to the author, the observed bubble diameter was equal to a few
millimetres, the bubbles remaining relatively small due to the strong
turbulence existing in the liquid phase.In our calculations, only a small part of the tube,
containing the singularity and the six measuring sections, was reproduced. The
radial profiles of the void fraction and the liquid mean and fluctuating
velocities measured upstream of the singularity were used as inlet conditions.
The length of the calculation domain is equal to 38 cm. The flow is assumed to
be axisymmetric therefore a two-dimensional axisymmetric meshing is used.
Several calculations have been done in order to test the sensitivity to the
axial and radial grids. Four different grids have been tested in [33],
the number of radial meshes in the largest section multiplied by the number of
axial meshes being 10*38, 20*76, 20*152, and 40*152, respectively. The
comparison of different calculations of the same two-phase flow, realized on
these four different grids show that the calculations performed with the finest
grid can be considered as converged. All the calculations presented here have
been done on the finest 40*152 calculation grid. The flow studied here is
characterized by the liquid and gas superficial velocities and the
area-averaged void fraction in the two sections given in Table 1. It can be
noted that the averaged void fraction has important values for this test (12%).Table 1
Simulated test case.
I.D.* (mm)JL (m/s)JG (m/s)〈α〉2501.570.30.121000.390.0750.0903*internal diameter.Table 2CsC1C2C1ωC2ωCεCε1Cε2Cε3Cε40.21.80.60.50.30.181.441.921.440.33Cokljat [34] performed calculations of the
sudden expansion experiment with the FLUENT code. Predictions were obtained
using the standard K-ε model as
well as an RSTM for the continuous phase, while the turbulence closure for the
dispersed phase is achieved by the algebraic model of Tchen [19, 20]. With this approach, similar as ours, the authors indicate that both
models produce similar results for the axial velocity but void fraction results
are improved with the RSTM model.We only consider the classical drag, added mass, and dispersion
turbulent force. The dispersion force coefficient is equal to 2 in the
computations. The bubble diameter is equal to 2 mm. Following [18, 33, 34], we underline the necessity to discard the
lift force. In fact, the effect of the lift force is
to produce sharp peaks near the wall because the classical modeling of the lift
force seems not well adapted to this case.The simulation results are in reasonable agreement with the experimental
data for the void fraction profiles (Figure5) and have been improved as compared to
[33, 34]. But the profiles at z=7 cm and
z=13 cm show that the void fraction is underestimated which mean that a
better understanding of the physical mechanisms is still needed.Figure 5
Void fraction profiles.Especially for the axial and radial mean liquid velocity profiles(Figure6,7), we
have obtained a good agreement which means that the recirculation zone is well
captured.Figure 6
Liquid mean radial velocity profiles.Figure 7
Liquid mean axial velocity profiles.We have obtained only qualitatively good results for the RMS quantities (Figure8)
because the turbulence mechanisms in a bubbly flow are far from being fully
understood [18]. But the turbulence modeling of the dispersed phase in
a PWR core in conditions close to nominal is less crucial than the liquid
turbulence modeling.Figure 8
Liquid r.m.s. radial (u′) and axial (w′) velocity profiles (Rij−ε).Finally, results with
theRij−ε turbulence model for the fluid are similar to the K−ε one,
which is our main objective in this case, before calculating rod bundle
geometries, where the K−ε model may fail.


## 4.1. The LIU&BANKOFF Case [25]

In this section, we evaluate the modeling capabilities
to simulate an upward bubbly flow in adiabatic conditions. We do not
specifically optimize the coefficients and modeling of the momentum transfer
terms to get results as good as possible. Our main objective is to validate theRij−ε turbulence model for the fluid against the K−ε one.The test section was aZ=2800 mm long, vertical
smooth acrylic tubing, with inside diameter of D=38 mm.The set of physical properties is the following:(24)ρliquid=994.9kg⋅m−3,ρgas=1.6kg⋅m−3,μliquid=7.97⋅10−4N⋅s⋅m−2,μgas=1.748⋅10−5N⋅s⋅m−2,P0=101500Pa,g=9.81m⋅s−2,dbubble=2.5⋅10−3m.  A uniform axial liquid profile is imposed at the inlet
and is equal to 1.138 m/s. A uniform axial gas profile is imposed at the inlet
and is equal to 1.333 m/s. The void fraction at the inlet is 0.045.The interfacial momentum transfer term is assumed to
be the sum of four different forces. The three first ones are simplified
averaged expressions of the classical drag, added mass, and lift force. The
fourth is the turbulent dispersion force.The flow is assumed to be axisymmetric therefore a
two-dimensional axisymmetric meshing is used. Computations have been performed
on two kinds of meshing: a coarse grid (20 cells in the radial direction and 50
cells in the axial direction) and a fine grid (30 cells in the radial direction
and 100 cells in the axial direction). Results are similar and computations are
performed on the first grid.At the measuring station (Z/D=36), we compare numerical results against
experimental data for the axial liquid velocity and the void fraction.As
recommended in [26], the lift coefficient is taken to be equal to
0.1. We have also tested the turbulent dispersion force of Davidson model
[26] written as (25)FDT=−34CDdbμt(1−α)|Vg−Vl|∂α∂xiwithμt=0.09ρlKl2ε. This expression gives values
negligible with respect to the Lopez de Bertodano expression [18] (the ratio is about 1000) and the void fraction profile is
observed to be similar to the profile calculated without any turbulence model
for the dispersed phase. If we take the bubble fluctuation into account with
the Hinze-Tchen algebraic model of bubble turbulence, we get the same results
with the turbulent dispersion force.In [27], Grossetête
considers that bubbles are deformed near the wall. To take into account this
effect, the author proposes to put a negative lift coefficient near the wall
and otherwise a positive one. Nevertheless, calculations (not presented in this
paper) show a peak of void fraction at the wall. We have also tested the
Tomiyama lift force [28, 29], but results are not improved.
Moreover, a wall lubrification force [30] can push the bubbles away
from the wall and improve the results.Figures3 and 4 show that the simulation results are in quite good
agreement with the experimental data. In [31], computations performed
with a K−ε turbulence model for the liquid produced comparable results. Our main
objective in this paper is to check that the simulation with the Rij−ε
turbulence model gives satisfactory results in a simple geometry, which is
crucial before calculating industrial geometries, where the K−ε turbulence model is susceptible to fail.Figure 3
Radial profile of the liquid velocity, comparisons of turbulent dispersion forces withCL=0.1.Figure 4
Radial profile of void fraction, comparisons of turbulent dispersion forces withCL=0.1.But improvement of the modeling of the interfacial
forces exerted on bubbles by the surrounding liquid is required. A strong
sensitivity to the lift coefficient has been found in our calculations. Other
forces, like the turbulent dispersion force, have also a crucial effect on the
void fraction distribution. These forces depend on uncontrolled parameters like
the bubble shape, the liquid turbulence, the bubbles collective effects, and so
on [31].


## 4.2. Sudden Expansion Experiment

Bel Fdhila [32] investigated
experimentally several upwards bubbly flows in a vertical pipe with a sudden
expansion. The total length of the pipe was equal to 14 meters. The bottom part
of the tube had a length equal to 9 meters and an internal diameter equal to 50 mm, the top part of the tube (5 m length) having an internal diameter equal to
100 mm. The fluids used were water and air under atmospheric pressure and
ambient temperature. Six measuring sections were located upstream and
downstream of the singularity. The first measuring section was located two
centimetres before the singularity, the other five were located at 7, 13, 18,
25, and 32 cm above the singularity. In each measuring section, the radial
profile of the void fraction was measured by means of a single optical probe,
and two components of the liquid velocity were measured by means of a hot film
anemometer. The time-averaged components of this velocity field and three
components of the liquid Reynolds stress tensor were deduced (the flow being
assumed axisymmetric). The bubble size was not measured in the experiment.
According to the author, the observed bubble diameter was equal to a few
millimetres, the bubbles remaining relatively small due to the strong
turbulence existing in the liquid phase.In our calculations, only a small part of the tube,
containing the singularity and the six measuring sections, was reproduced. The
radial profiles of the void fraction and the liquid mean and fluctuating
velocities measured upstream of the singularity were used as inlet conditions.
The length of the calculation domain is equal to 38 cm. The flow is assumed to
be axisymmetric therefore a two-dimensional axisymmetric meshing is used.
Several calculations have been done in order to test the sensitivity to the
axial and radial grids. Four different grids have been tested in [33],
the number of radial meshes in the largest section multiplied by the number of
axial meshes being 10*38, 20*76, 20*152, and 40*152, respectively. The
comparison of different calculations of the same two-phase flow, realized on
these four different grids show that the calculations performed with the finest
grid can be considered as converged. All the calculations presented here have
been done on the finest 40*152 calculation grid. The flow studied here is
characterized by the liquid and gas superficial velocities and the
area-averaged void fraction in the two sections given in Table 1. It can be
noted that the averaged void fraction has important values for this test (12%).Table 1
Simulated test case.
I.D.* (mm)JL (m/s)JG (m/s)〈α〉2501.570.30.121000.390.0750.0903*internal diameter.Table 2CsC1C2C1ωC2ωCεCε1Cε2Cε3Cε40.21.80.60.50.30.181.441.921.440.33Cokljat [34] performed calculations of the
sudden expansion experiment with the FLUENT code. Predictions were obtained
using the standard K-ε model as
well as an RSTM for the continuous phase, while the turbulence closure for the
dispersed phase is achieved by the algebraic model of Tchen [19, 20]. With this approach, similar as ours, the authors indicate that both
models produce similar results for the axial velocity but void fraction results
are improved with the RSTM model.We only consider the classical drag, added mass, and dispersion
turbulent force. The dispersion force coefficient is equal to 2 in the
computations. The bubble diameter is equal to 2 mm. Following [18, 33, 34], we underline the necessity to discard the
lift force. In fact, the effect of the lift force is
to produce sharp peaks near the wall because the classical modeling of the lift
force seems not well adapted to this case.The simulation results are in reasonable agreement with the experimental
data for the void fraction profiles (Figure5) and have been improved as compared to
[33, 34]. But the profiles at z=7 cm and
z=13 cm show that the void fraction is underestimated which mean that a
better understanding of the physical mechanisms is still needed.Figure 5
Void fraction profiles.Especially for the axial and radial mean liquid velocity profiles(Figure6,7), we
have obtained a good agreement which means that the recirculation zone is well
captured.Figure 6
Liquid mean radial velocity profiles.Figure 7
Liquid mean axial velocity profiles.We have obtained only qualitatively good results for the RMS quantities (Figure8)
because the turbulence mechanisms in a bubbly flow are far from being fully
understood [18]. But the turbulence modeling of the dispersed phase in
a PWR core in conditions close to nominal is less crucial than the liquid
turbulence modeling.Figure 8
Liquid r.m.s. radial (u′) and axial (w′) velocity profiles (Rij−ε).Finally, results with
theRij−ε turbulence model for the fluid are similar to the K−ε one,
which is our main objective in this case, before calculating rod bundle
geometries, where the K−ε model may fail.


## 5. Conclusion and Perspectives

An analysis of turbulence modeling for two-phase flows
has been proposed. Indeed, the use of eddy viscosity models is widespread and
may be sufficient for parallel flows in vertical pipes, but that type of model
does not account for effects that are preponderant in complex geometries,
especially when swirling flows are involved, for example, in pressurized water
reactor cores downstream of mixing vanes and spacer grids. In accordance with
the theory, it is demonstrated in the case of a flow downstream of a mixing
vane that using a Reynolds stress model is an efficient way to improve the
simulation of such complex flows. To demonstrate that the use of a Reynolds stress
model is not bound to deteriorate the classical results obtained with an eddy
viscosity model, a validation step on more analytical experiments is detailed
(bubbly flows in a straight pipe and in a sudden expansion): the study shows
that the Reynolds Stress model implemented in the multiphase 3D code
Neptune_CFD satisfactorily reproduces the results obtained with the standard
eddy viscosity model and both compare reasonably well with the experiments.As concern the computational cost, we note that in the
case of the DEBORA-mixing test which is under process, the time required by
iteration is, respectively, 3.09 seconds and 2.81 seconds for RSTM and EVM. The
time step is, respectively, 5 milliseconds and 5.4 milliseconds with a CFL
equal to 1. In this particular case, the RSTM over-cost is about 18.8%.Moreover, among the developments planned in the medium
term, we have identified the need for a polydispersion model.Besides,
the Neptune project has set up a medium and
long-term experimental program to acquire detailed measurements in simplified
and real geometries, both in adiabatic and real conditions [9, 10].



---

*Source: 792395-2008-10-14.xml*