# Euler-Euler Large Eddy Simulation of a Square Cross-Sectional Bubble Column Using the Neptune_CFD Code

Authors: B. Ničeno; M. Boucker; B. L. Smith
Journal: Science and Technology of Nuclear Installations
 (2008-02-13)

---

## Abstract

In this work, we report on Euler-Euler large eddy simulation (EELES) of dispersed bubbly flow in a square cross-sectional bubble column. Simulations are performed using the Neptune_CFD package, and results processed using the SALOME platform. The motivation to undertake this study is to check our implementation of the Smagorinsky subgrid-scale (SGS) model into NEPTUNE_CFD. We outline all the physical models used, and we present instantaneous realizations of velocity and void fraction fields in order to illustrate the structure of the turbulence field, and long-time averaged results, to compare with analogous simulations performed using the CFX-4 code and experimental data. The same physical models and constants have been used in both the CFX-4 and NEPTUNE_CFD codes, except the SGS model, which is Smagorinsky in case of NEPTUNE_CFD and a one-equation model in CFX-4. The results obtained with EELES compare reasonably well with experiment, meaning in particular that the implementations have been successful. Some perspectives on the further use of EELES are also given.


---

## Body

## 1. Introduction

The aim of the NURESIM
[1] integrated project is
to provide the initial step towards a common European
software platform for analysing, modelling, and exchanging data for nuclear
reactor simulations. One of the main goals of the project is the integration of
state-of-the-art physical models in a common, open-software platform, including
the latest advances in reactor core physics, thermal-hydraulics, and coupled multiphysics
modelling capabilities. The objective of the thermal-hydraulic subproject within NURESIM is the
improvement of the predictive capabilities of the simulation tools for key
two-phase-flow, thermal-hydraulic processes that can occur in nuclear reactors,
focusing on two high-priority
issues: critical heat flux (CHF) [2] and pressurized thermal shock (PTS) [3]. As a part of
this activity, the Paul Scherrer Institute is investigating the potential of
the EELES approach for simulating bubbly flows, which are important for understanding
CHF and PTS, using the Neptune_CFD code [4], the CFD component of the NURESIM
platform. In view of this, it was thought desirable
to carry out EELES in a three-dimensional, square-cross-sectional bubble
column, based on the experiments of Deen et al. [5], to validate the
implementation of the SGS model in the code. In the present work, for reasons
of simplicity, the performance of Neptune_CFD with the EELES paradigm is
assessed based on the Smagorinsky SGS model. These results, and those obtained
with the CFX-4 code, are compared against experimental findings [5]. We point
to the potential weaknesses of the current approach, and propose directions for future research and development in the area.Previous work on EELES (often
referred to as two-fluid LES) for bubbly flows include Milelli et al. [6], who
used Smagorinsky and dynamic SGS models to study the motion of a bubble plume
in a cylindrical test section 50 cm
in diameter, 40 cm high, for a maximum superficial
gas velocity of 6.11 mm/s. A similar approach was taken by Deen et al. [7] to
simulate a bubble plume in a test section with a square cross-section of 15 cm × 15 cm,
a height of 45 cm, and superficial gas velocity of 4.9 mm/s. The authors
compare the EELES prediction with those obtained using a K−ε model and conclude that EELES results are
closer to experimental measurements. Turbulent bubbly shear flows have also
been analyzed [8] for a square channel 30 cm wide, 4 cm deep, and 60 cm high. A
splitter plane was used to divide the inlet channel in two equal parts. A mixture
of water and air, with different velocities and void fractions, was introduced
at the inlet sections. More
recently, an LES calculation incorporating the Smagorinsky SGS model has been
reported for a rectangular column [9]. LES predictions were compared with transient
results obtained using a mixing length model and the K−ε model; all
three models produced similar results.


### 1.1. The Deen Bubble Column Experiment

In the Deen series
of experiments [5], the hydrodynamics of the bubble plume is studied in the
same apparatus, but for different inlet conditions. In the test to be
investigated here, bubbles are created by injecting air with a superficial gas velocity
of 4.9 mm/s through a central sparger at the bottom of the vessel of dimension
5 cm × 5 cm.
The experiment operates at atmospheric pressure. A sketch of the Deen
experiment is given in Figure 1.Figure 1
Geometry of the Deen bubble column
experiment.


## 1.1. The Deen Bubble Column Experiment

In the Deen series
of experiments [5], the hydrodynamics of the bubble plume is studied in the
same apparatus, but for different inlet conditions. In the test to be
investigated here, bubbles are created by injecting air with a superficial gas velocity
of 4.9 mm/s through a central sparger at the bottom of the vessel of dimension
5 cm × 5 cm.
The experiment operates at atmospheric pressure. A sketch of the Deen
experiment is given in Figure 1.Figure 1
Geometry of the Deen bubble column
experiment.


## 2. Physical Models

### 2.1. Euler-Euler Approach

The governing
equation for conservation of mass, with no mass exchange between the phases,
takes the form:(1)∂∂t(rρ)α+∇⋅(rρu)α=0,
in which t is time, r is the volume fraction, ρ is the
phase density, u the phase velocity,
and α is the phase indicator. Conservation of momentum is governed
by the relation:
(2)∂(rρu)α∂t+∇⋅(rρuu)α=−∇⋅(rτ)α−rα∇p+(rρ)αg+MF,α,
where  τ is shear
stress, p is pressure, and g is the gravity vector. The last
term on the right-hand
side represents the momentum exchange between the phases due to interface
forces and requires modelling. The models
we have used are given in Section 2.3.


### 2.2. LES of Turbulence

Turbulence is
modelled in the liquid phase of the flow. The first step in the derivation of
the governing equations for LES is the decomposition of the velocity field into
resolved (grid scale) and unresolved (SGS) parts. For each component, we write(3)uL=u˜L−uL′,
where the
resolved velocity on the left side is equal to the difference between the true
instantaneous velocity (first term on the right) and the unresolved velocity (second
term on the right). The filtering of the nonlinear
advection term on the left-hand side of (2)
leads to an additional stress-like term, which has to be modelled. Hence, we
write
(4)τL=−2(μS)L+τSGS,L.The first term
represents the laminar stress, with the strain rate tensor defined as
(5)SL=12[∇uL+(∇uL)T],and the second
term derives from the turbulence effects and needs closure. In this work, we
have used two closure models: Smagorinsky model (in Neptune_CFD), and a one-equation
model (in CFX-4).


### 2.3. Smagorinsky Model

To complete (4), we use a formulation [10] for
the turbulent stress:(6)τSGS,L−13(τkkδij)L=−2(μtS)L.The second term
on the left-hand
side, the trace of the SGS stress tensor, is implicitly added to the pressure,
while the deviatoric part is modelled by the expression on the right-hand
side, that is, by the Boussinesq eddy viscosity (μt) concept. The Smagorinsky model is an
algebraic one, with the eddy viscosity calculated from:
(7)μt,L=ρL(CSΔ)2|SL|.In (7), Δ is a
filter width, associated with the cell size. In the present work, it was
estimated as the cube root of the cell volume:
(8)Δ=V1/3,although other
definitions are possible. CS is the Smagorinsky
constant, and was set here to 0.12 since the
flow is dominated by the large-scale structures generated by the buoyancy force,
as pointed out in [11]. The eddy viscosity, defined by (7),
is used in the stress tensor appearing in the momentum (4),
giving
(9)τL=−2[(μ+μt)S]L.


### 2.4. One-Equation SGS Model

A disadvantage of the Smagorinsky model is the loss of information resulting
from use of the deviatoric part of the SGS stress tensor only. In CFX-4, we
have implemented a transport equation for SGS kinetic energy(KSGS), defined by
(10)kSGS,L=12ρLτkk,L,of the following transport form [12]:
(11)∂(rLρLkSGS)∂t+∇⋅(rLρLuLkSGS)=∇[rL(μL+μt,L)∇kSGS]=+rLPkSGS−rLρLCεkSGS3/2Δ,where PkSGS is the turbulence production term defined as:
(12)Pk,LSGS=(μtS:S)L.
Using this model, the eddy viscosity is then calculated from:
(13)μt,L=CkρLΔkSGS,L1/2,
with model constants Ck=0.07 and Cε=1.05. A detailed description of the one-equation SGS model for dispersed bubbly flows
can be found in [12].


### 2.5. Interfacial Forces

To close the
momentum equation set for the two phases, the various interphase exchange
terms have to be modelled. The  interphase exchange
terms relevant for
the configuration considered in this work aredrag, buoyancy, virtual mass, and lift forces. In the framework of the EELES presented in this work,
interfacial forces are applied only to the resolved field, while their SGS
counterparts are neglected.Thedrag force is modelled from the resolved
velocity field as follows [13]:
(14)MD,L=−MD,G=34rGρLCDdG|uG−uL|(uG−uL),
where CD is the drag coefficient, defined here
in terms of the Eötvös number (EO=gΔρdG2/σ):(15)CD=23E01/2.
In the present
application, a bubble size of 4.0 mm is specified (the value reported in [5]), and
the coefficient of surface tension σ=0.072 N/m, giving EO=2.2 and CD=1.0.For thevirtual mass force, we use the model of Drew and Lahey Jr. [14]:(16)MVM,L=−MVM,G=rGρLCVM(DuGDt−DuLDt),
where CVM is the
virtual mass force coefficient, here set to 0.5.Thelift force accounts for the transverse
migration of bubbles under the influence of the liquid shear field. Following
earlier work [14, 15], it is modelled here according to(17)ML,L=−ML,G=ρLCL(uG−uL)×∇×uL,
in which the
lift coefficient, CL, is
set to 0.5.


## 2.1. Euler-Euler Approach

The governing
equation for conservation of mass, with no mass exchange between the phases,
takes the form:(1)∂∂t(rρ)α+∇⋅(rρu)α=0,
in which t is time, r is the volume fraction, ρ is the
phase density, u the phase velocity,
and α is the phase indicator. Conservation of momentum is governed
by the relation:
(2)∂(rρu)α∂t+∇⋅(rρuu)α=−∇⋅(rτ)α−rα∇p+(rρ)αg+MF,α,
where  τ is shear
stress, p is pressure, and g is the gravity vector. The last
term on the right-hand
side represents the momentum exchange between the phases due to interface
forces and requires modelling. The models
we have used are given in Section 2.3.


## 2.2. LES of Turbulence

Turbulence is
modelled in the liquid phase of the flow. The first step in the derivation of
the governing equations for LES is the decomposition of the velocity field into
resolved (grid scale) and unresolved (SGS) parts. For each component, we write(3)uL=u˜L−uL′,
where the
resolved velocity on the left side is equal to the difference between the true
instantaneous velocity (first term on the right) and the unresolved velocity (second
term on the right). The filtering of the nonlinear
advection term on the left-hand side of (2)
leads to an additional stress-like term, which has to be modelled. Hence, we
write
(4)τL=−2(μS)L+τSGS,L.The first term
represents the laminar stress, with the strain rate tensor defined as
(5)SL=12[∇uL+(∇uL)T],and the second
term derives from the turbulence effects and needs closure. In this work, we
have used two closure models: Smagorinsky model (in Neptune_CFD), and a one-equation
model (in CFX-4).


## 2.3. Smagorinsky Model

To complete (4), we use a formulation [10] for
the turbulent stress:(6)τSGS,L−13(τkkδij)L=−2(μtS)L.The second term
on the left-hand
side, the trace of the SGS stress tensor, is implicitly added to the pressure,
while the deviatoric part is modelled by the expression on the right-hand
side, that is, by the Boussinesq eddy viscosity (μt) concept. The Smagorinsky model is an
algebraic one, with the eddy viscosity calculated from:
(7)μt,L=ρL(CSΔ)2|SL|.In (7), Δ is a
filter width, associated with the cell size. In the present work, it was
estimated as the cube root of the cell volume:
(8)Δ=V1/3,although other
definitions are possible. CS is the Smagorinsky
constant, and was set here to 0.12 since the
flow is dominated by the large-scale structures generated by the buoyancy force,
as pointed out in [11]. The eddy viscosity, defined by (7),
is used in the stress tensor appearing in the momentum (4),
giving
(9)τL=−2[(μ+μt)S]L.


## 2.4. One-Equation SGS Model

A disadvantage of the Smagorinsky model is the loss of information resulting
from use of the deviatoric part of the SGS stress tensor only. In CFX-4, we
have implemented a transport equation for SGS kinetic energy(KSGS), defined by
(10)kSGS,L=12ρLτkk,L,of the following transport form [12]:
(11)∂(rLρLkSGS)∂t+∇⋅(rLρLuLkSGS)=∇[rL(μL+μt,L)∇kSGS]=+rLPkSGS−rLρLCεkSGS3/2Δ,where PkSGS is the turbulence production term defined as:
(12)Pk,LSGS=(μtS:S)L.
Using this model, the eddy viscosity is then calculated from:
(13)μt,L=CkρLΔkSGS,L1/2,
with model constants Ck=0.07 and Cε=1.05. A detailed description of the one-equation SGS model for dispersed bubbly flows
can be found in [12].


## 2.5. Interfacial Forces

To close the
momentum equation set for the two phases, the various interphase exchange
terms have to be modelled. The  interphase exchange
terms relevant for
the configuration considered in this work aredrag, buoyancy, virtual mass, and lift forces. In the framework of the EELES presented in this work,
interfacial forces are applied only to the resolved field, while their SGS
counterparts are neglected.Thedrag force is modelled from the resolved
velocity field as follows [13]:
(14)MD,L=−MD,G=34rGρLCDdG|uG−uL|(uG−uL),
where CD is the drag coefficient, defined here
in terms of the Eötvös number (EO=gΔρdG2/σ):(15)CD=23E01/2.
In the present
application, a bubble size of 4.0 mm is specified (the value reported in [5]), and
the coefficient of surface tension σ=0.072 N/m, giving EO=2.2 and CD=1.0.For thevirtual mass force, we use the model of Drew and Lahey Jr. [14]:(16)MVM,L=−MVM,G=rGρLCVM(DuGDt−DuLDt),
where CVM is the
virtual mass force coefficient, here set to 0.5.Thelift force accounts for the transverse
migration of bubbles under the influence of the liquid shear field. Following
earlier work [14, 15], it is modelled here according to(17)ML,L=−ML,G=ρLCL(uG−uL)×∇×uL,
in which the
lift coefficient, CL, is
set to 0.5.


## 3. Requirements for the Numerical Grid

In the present
work, we couple the Euler-Euler approach for multiphase flow with LES, and
therefore have to consider the resolution requirements of both techniques simultaneously
in order to choose a satisfactory grid. A basic requirement is that the control
volume size should be large enough to encompass all the interface details. This
is the intrinsic assumption in the derivation of the Euler-Euler model
equations, and strictly has to be satisfied at the discrete level as well.In LES, the SGS
model is often very simple, and only drains energy from the resolved field without feed-back.
Therefore, our goal in LES is to resolve as much of the flow field as possible,
and to have as fine a grid as feasible for the available computer hardware. Since
the Euler-Euler approach specifies the minimum control volume size, whereas for
LES we are invariantly seeking as fine a grid as possible, the requirements for
the numerical grid may sometimes be in conflict. The point is illustrated in
Figures2 and 3, which show the turbulence spectra under different modelling conditions. For
a successful LES, we must have a filter width (Δ) in the inertial subrange
region, and all scales of motion larger than that (left of Δ in Figures 2 and 3), must be accurately resolved on the
numerical grid. If, however, we have a bubble diameter (dG) larger than Δ (see Figure 2), it is obvious that they would
induce some large-scale motions which are not properly accounted for by the
EELES, since there is no information on the interface details or their
influence on the resolved large-scale motions. If in addition we use a model
for bubble-induced turbulence, it would drain the energy from the resolved
field, further deteriorating the accuracy of the resolved field. This is
illustrated by the saw-like shading in Figure 2. This influence on the resolved part of the spectra is not acceptable for
LES. This is not just a conceptual consideration, as discussed in [6], and it is explained in more detail in Section 3.1
. The situation which is safe for EELES is shown in Figure 3, where dG is smaller than Δ, and all bubble-induced scales which cannot be
calculated using EELES approach fall into the SGS part of the spectra.Figure 2
Bubble size larger than filter
width.Figure 3
Bubble size smaller than filter width.


### 3.1. The Milelli Condition

The grid size
considerations discussed above, are not new. In the present work, we have indicated why the bubbles must be
smaller than the cell size from the point of view of the energy spectra and
modelling closures for the interfacial forces. A systematic a posteriori analysis of the minimum
ratio of the bubble and cell sizes for LES modelling of free bubble plumes is
reported by Milelli et al. [6], which resulted in the following criterion:(18)hdG≥1.5.
Equation (18)
states that the cell size must be at least 50% larger than the bubble diameter
for accurate LES. The situation is illustrated in Figure 4. In the Deen experiment, the mean bubble size is 4 mm. Since the flow is
dominated by the energetic, large-scale structures in the core of the flow,
with wall effects having a smaller impact on the overall flow field, we have
created a uniform grid with 30 × 30 × 100 cubic cells. This results in a cell
size of 5 mm, hence, (19)hdG=1.2,which is below
the optimum value. A coarser mesh which satisfies the Milelli condition has also
been used, but with no significant change in the computed results [12].Figure 4
The Milelli condition.It might be argued at this point that the grid used in
the present simulations is too coarse for LES, at least from the point of view
of capturing all the relevant (large) scales. Judging from the flow patterns
that can be expected in a square column such as this, with the bubble plume
meandering from one direction to the other, it is quite conceivable that the
largest, and therefore the most energetic, eddies will be of the size of the
domain cross-section. Therefore, we are confident that the grid resolution we
employ (30× 30 cells in the cross-section) is sufficient to resolve these
large structures. Should the topology of the flow be different, and should the
flow be dominated by small, near-wall scales of motion, this resolution would
not be adequate. The adequacy of the grid resolution used is proved in the
following section.


## 3.1. The Milelli Condition

The grid size
considerations discussed above, are not new. In the present work, we have indicated why the bubbles must be
smaller than the cell size from the point of view of the energy spectra and
modelling closures for the interfacial forces. A systematic a posteriori analysis of the minimum
ratio of the bubble and cell sizes for LES modelling of free bubble plumes is
reported by Milelli et al. [6], which resulted in the following criterion:(18)hdG≥1.5.
Equation (18)
states that the cell size must be at least 50% larger than the bubble diameter
for accurate LES. The situation is illustrated in Figure 4. In the Deen experiment, the mean bubble size is 4 mm. Since the flow is
dominated by the energetic, large-scale structures in the core of the flow,
with wall effects having a smaller impact on the overall flow field, we have
created a uniform grid with 30 × 30 × 100 cubic cells. This results in a cell
size of 5 mm, hence, (19)hdG=1.2,which is below
the optimum value. A coarser mesh which satisfies the Milelli condition has also
been used, but with no significant change in the computed results [12].Figure 4
The Milelli condition.It might be argued at this point that the grid used in
the present simulations is too coarse for LES, at least from the point of view
of capturing all the relevant (large) scales. Judging from the flow patterns
that can be expected in a square column such as this, with the bubble plume
meandering from one direction to the other, it is quite conceivable that the
largest, and therefore the most energetic, eddies will be of the size of the
domain cross-section. Therefore, we are confident that the grid resolution we
employ (30× 30 cells in the cross-section) is sufficient to resolve these
large structures. Should the topology of the flow be different, and should the
flow be dominated by small, near-wall scales of motion, this resolution would
not be adequate. The adequacy of the grid resolution used is proved in the
following section.


## 4. Simulation Details

In this
section, we briefly summarize all the relevant simulation details.Inlet Superficial Velocities:
for liquid:JL = 0.0 m/s;
for gas:JG = 0.0049 m/s.Spatial Discretization:
number of cells in computational domain: 90 000;
advection scheme: bounded central scheme for Neptune_CFD;
quick scheme for
CFX-4.Temporal
Discretization:
the linear backward time differencing for Neptune_CFD;
quadratic backward differencing for CFX-4;
CFL = 1.0, variable time step,
approximately equal to 0.05 second for Neptune_CFD constant time step of 0.05 second for CFX-4, leading to CFL≈  1.0;
integral time of simulation: 400 seconds;
simple time-average is used in both, not
weighted by volume fraction.


## 5. Results

In this
section, the results for our EELES with Neptune_CFD are presented. We start
with the instantaneous velocity and void fraction fields to show the flow
patterns in the bubble column flow, and continue with time-averaged velocity and fluctuating liquid velocity
profiles. Time-averaged profiles obtained with Neptune_CFD and CFX-4 are
compared with the experimental measurements of Deen et al. [5].


### 5.1. Instantaneous Results

Instantaneous
velocity fields are shown in Figure5. Velocity vectors are coloured
with the velocity magnitude. Two things are immediately apparent. First, the
velocity field does not show any signs of the growing instabilities to which
LES computations are very often prone. This is not surprising,
since the advection scheme is bounded. Another thing
worth noting is that the velocity field pattern is changing dramatically in
time. Bubbles are moving up through the large turbulent structures, driven by
their buoyancy, resulting in different bubble paths at each time instant. This
behaviour was also reported for bubble plume
experiments [16]. Figure 6 shows two isosurfaces of
constant void fraction. The blue one corresponds to the low-void fraction of r=.05
and the red to r=.5. The threshold for the blue isosurface gives a
good illustration of the bubble plume shape, whereas the red one indicates the
position of the free surface.Instantaneous velocity vectors, calculated using Neptune_CFD, at
instants: 60 seconds, 80 seconds, and 120 seconds.
(a)(b)(c)Isosurfaces of constant void fraction (0.05 and 0.5), calculated using Neptune_CFD, at instants: 60 seconds, 80 seconds, and 120 seconds.
(a)(b)(c)


### 5.2. Time-Averaged Results

In this
section, we compare the simulated profiles of liquid and gas vertical
velocities and liquid turbulence intensities against the reported experimental
data [5]. The measurements were taken along the centreline of the horizontal
plane at height 25 cm.Figure7 shows
comparisons of vertical liquid and gas vertical velocity profiles. The liquid
velocity profile is somewhat under-predicted by both codes. In contrast, the
gas velocity prediction compares very well with experiment in both cases. From
these time-averaged results, it might be concluded that the drag force is
underestimated in both codes, probably due to a wrong assumption for the drag
coefficient (CD), or for
the relative velocity between the two phases [16].Comparison of (a) liquid and (b) gas velocity profiles obtained with CFX-4 using a one-equation model, Neptune_CFD with the Smagorinsky model, and experimental data.
(a)(b)The vertical
component of the resolved vertical fluctuating liquid velocity is plotted in
Figure8(a).
The qualitative comparison is encouraging for both codes, since the twin-peaked
shape seen in the experiments is reproduced. Quantitatively, however, both codes under-predict the magnitude to some extent. Figure 8(b)shows the resolved lateral fluctuating liquid velocity, which is predicted very
well. The total turbulent kinetic energy, plotted in Figure 9, in spite of the under-prediction
of vertical fluctuating liquid component, is well predicted by both codes in
the middle of the column and somewhat over-predicted close to the walls. This
over-estimation of the near-wall kinetic energy might be attributed to
insufficient resolution for the near-wall structures.Comparison of (a) vertical and (b) lateral fluctuating liquid velocities, obtained using CFX-4 with a one-equation model, Neptune_CFD with the Smagorinsky model, and experimental data.
(a)(b)Figure 9
Comparison of liquid
turbulent kinetic energy obtained with CFX-4 using a one-equation model,
Neptune_CFD with the Smagorinsky model, and experimental data. The blue
dashed line is the resolved, the blue continuous line is the total (resolved
plus SGS) kinetic energy.


## 5.1. Instantaneous Results

Instantaneous
velocity fields are shown in Figure5. Velocity vectors are coloured
with the velocity magnitude. Two things are immediately apparent. First, the
velocity field does not show any signs of the growing instabilities to which
LES computations are very often prone. This is not surprising,
since the advection scheme is bounded. Another thing
worth noting is that the velocity field pattern is changing dramatically in
time. Bubbles are moving up through the large turbulent structures, driven by
their buoyancy, resulting in different bubble paths at each time instant. This
behaviour was also reported for bubble plume
experiments [16]. Figure 6 shows two isosurfaces of
constant void fraction. The blue one corresponds to the low-void fraction of r=.05
and the red to r=.5. The threshold for the blue isosurface gives a
good illustration of the bubble plume shape, whereas the red one indicates the
position of the free surface.Instantaneous velocity vectors, calculated using Neptune_CFD, at
instants: 60 seconds, 80 seconds, and 120 seconds.
(a)(b)(c)Isosurfaces of constant void fraction (0.05 and 0.5), calculated using Neptune_CFD, at instants: 60 seconds, 80 seconds, and 120 seconds.
(a)(b)(c)


## 5.2. Time-Averaged Results

In this
section, we compare the simulated profiles of liquid and gas vertical
velocities and liquid turbulence intensities against the reported experimental
data [5]. The measurements were taken along the centreline of the horizontal
plane at height 25 cm.Figure7 shows
comparisons of vertical liquid and gas vertical velocity profiles. The liquid
velocity profile is somewhat under-predicted by both codes. In contrast, the
gas velocity prediction compares very well with experiment in both cases. From
these time-averaged results, it might be concluded that the drag force is
underestimated in both codes, probably due to a wrong assumption for the drag
coefficient (CD), or for
the relative velocity between the two phases [16].Comparison of (a) liquid and (b) gas velocity profiles obtained with CFX-4 using a one-equation model, Neptune_CFD with the Smagorinsky model, and experimental data.
(a)(b)The vertical
component of the resolved vertical fluctuating liquid velocity is plotted in
Figure8(a).
The qualitative comparison is encouraging for both codes, since the twin-peaked
shape seen in the experiments is reproduced. Quantitatively, however, both codes under-predict the magnitude to some extent. Figure 8(b)shows the resolved lateral fluctuating liquid velocity, which is predicted very
well. The total turbulent kinetic energy, plotted in Figure 9, in spite of the under-prediction
of vertical fluctuating liquid component, is well predicted by both codes in
the middle of the column and somewhat over-predicted close to the walls. This
over-estimation of the near-wall kinetic energy might be attributed to
insufficient resolution for the near-wall structures.Comparison of (a) vertical and (b) lateral fluctuating liquid velocities, obtained using CFX-4 with a one-equation model, Neptune_CFD with the Smagorinsky model, and experimental data.
(a)(b)Figure 9
Comparison of liquid
turbulent kinetic energy obtained with CFX-4 using a one-equation model,
Neptune_CFD with the Smagorinsky model, and experimental data. The blue
dashed line is the resolved, the blue continuous line is the total (resolved
plus SGS) kinetic energy.


## 6. Conclusions and Perspectives

Implementation
of the Smagorinsky SGS model in the Neptune_CFD code has been validated using
data from the Deen bubble column experiment. Results have also been compared
with predictions obtained from simulations performed using the code CFX-4.
Generally, liquid velocity profiles are under-predicted compared with measured
data, whereas the gas velocity profiles are predicted very well. This can only
be attributed to a wrong assumption for the drag coefficient. The liquid velocity
fluctuations predicted by both codes compare well with experiments, except for
some under-prediction of the vertical component. Overall, the results obtained
with the two codes are consistent, regardless of the fact that different SGS
models have been
used. This is a clear indication that the SGS model influences results only to
a very small extent.When envisaging
the potential of LES for modelling multiphase flows in general, one should
clearly distinguish between the scales at which LES might be used, since this
implies the level of detail of accurate interface resolution to the degree of modelling
that can be tolerated. Simulations at
mesoscales imply an
Euler-Euler description of the interface between the phases, such as the one
described in this work. However, if LES for multiphase flows is applied at microscales,
explicit interface tracking procedures would be needed. This has already been
proposed in [17], and is generally referred to as large-scale
simulation (LSS).The principal
advantage of EELES over LSS is that since the interface details are not
calculated explicitly, the simulations may be carried out at lower cost. The
principal disadvantage of EELES is that the most influential interfacial forces
(lift and drag) are modelled for the large-scale field, meaning that the
question of how to model these forces remains as open for EELES as it is with RANS. LSS, on the
other hand, explicitly resolves the large-scale part of the interfacial forces,
leaving the modelling at the SGS level, where the effects are smaller and hence
less influential on the accuracy of the results. LSS modelling comes at a heavy
price, due to the need to resolve all relevant interface details, imposing huge
demands on computing power. Applying LSS to industrial-scale problems is beyond
the current state of computing resources.Another
disadvantage of EELES stems from the different resolution requirements imposed
by the Euler-Euler description of the two fluids and that of the LES approach itself,
expressed in [6]. Since for an Euler-Euler description, the minimum cell size must be larger than the interface detail, accurate LES
requires as fine a grid as possible. This requirement is particularly stringent
in the near-wall regions, where the turbulent structures (e.g., streaks) are very small, and a very
fine grid is needed to capture all relevant details. The disadvantages of EELES
just outlined are not so pronounced for flows in which the bubbles are small
compared to the grid size, such as for the bubble column studied in this work.
Generally, one should use EELES if the flow is likely to be unsteady, with
large coherent structures generated either by buoyancy or obstacles in the
flow. EELES, as LES itself, should be avoided for wall-dominated flows. In
other words, large
eddy simulation should be used in the cases in which large eddies are present.
Flows with density stratification in large channels (relevant to the PTS
issue), in addition to bubble column devices, are examples where EELES has the
potential to give accurate insight into the phenomena taking place and the
ability to quantify them.



---

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