Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 12.4.2 - Direct Numerical Simulation of a Liquid-Solid Fluidized Bed
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00306
 Material Information
Title: 12.4.2 - Direct Numerical Simulation of a Liquid-Solid Fluidized Bed Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Corre, C.
Estivalezes, J.L.
Vincent, S.
Simonin, O.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: direct numerical simulation
fluidized bed
Eulerian-Lagrangian method
particle tracking
 Notes
Abstract: The characterization of fluidized beds still requires specific investigation for understanding and modeling the local coupling between the dispersed phase and the carrier fluid. The aim of this work is to simulate this type of unsteady particle laden flows via Direct Numerical Simulations in order to provide a local and instantaneous description of particle flow interactions and to extract statistical parameters useful for large scale models. With the advantage of powerful supercomputers and the development of new numerical methods, direct numerical simulation of real fluidized beds involving thousands of particles becomes possible with an Eulerian treatment of the particulate phase. In this work, a fictitious domain model and an Eulerian description of the two-phase flow is investigated, the particle being larger than the local space step of the mesh. The originality of the numerical methods is based on a one-fluid formulation of the incompressible Navier Stokes equations in which the presence of the particles is tackle with a penalty method to ensure the solid behavior (Randrianarivelo 2005) and (Randrianarivelo 2007). The particle tracking is obtained via an original hybrid Eulerian-Lagrangian Volume Of Fluid (VOF) method, previously presented by Randrianarivelo (2005) and Corre (2010). Particles are managed in a Lagrangian manner whereas the local solid fraction is described in an Eulerian way. In addition, collisions between particles are treated via a soft sphere model of Cundall (1979). At each time step, the particle characteristics are stored and statistics are performed on the particle-particle and fluid-particle interactions. A structured Cartesian grid is used and the cylindrical tank in which the fluidized bed occurs is represented by a penalty method. The simulated cylinder tank corresponds to the apparatus used in the experimental study of Aguilar Corona (2008) but at a half scale. Numerical results are compared with available experimental results for validation of the numerical modeling. Instantaneous and average flow characteristics are reported, as well as the pair distribution function, the streamlines of the mean particle velocity field, the lagrangian particle auto-correlation coefficients as well as the auto-diffusion coefficents. The results agree with the literature, notably concerning the maximum of the pair distribution function, good trends are also observed concerning the order of magnitude of the lagrangian auto-correlation coefficients, the variance of particle movement and the auto-diffusion coefficients.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00306
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1242-Corre-ICMF10.pdf

Full Text



7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Direct Numerical Simulation of a Liquid-Solid Fluidized Bed


C. Corre*, J.L. Estivalezes*, S. Vincentti and 0. SimoninT

ONERA/DMAE, 2 avenue Edouard Belin, BP 74055, 31055 Toulouse, France
t TREFLE-ENSCPB, UMR 8508 CNRS, 16 avenue Pey-Berland, 33607 Pessac Cedex, France
University de Toulouse ; INPT, UPS ; IMFT ; All6e Camille Soula, F-31400 Toulouse, France
SCNRS; TREFLE ; F-33407 Pessac, France
CNRS; IMFT ; F-31400 Toulouse, France
cedric.corre@onera.fr, Jean-Luc.Estivalezes@onera.fr, vincent@enscpb.fr and simonin@imft.fr
Keywords: direct numerical simulation, fluidized bed, Eulerian-Lagrangian method, particle tracking




Abstract

The characterization of fluidized beds still requires specific investigation for understanding and modeling the local
coupling between the dispersed phase and the carrier fluid. The aim of this work is to simulate this type of unsteady
particle laden flows via Direct Numerical Simulations in order to provide a local and instantaneous description of
particle flow interactions and to extract statistical parameters useful for large scale models. With the advantage
of powerful supercomputers and the development of new numerical methods, direct numerical simulation of real
fluidized beds involving thousands of particles becomes possible with an Eulerian treatment of the particulate
phase. In this work, a fictitious domain model and an Eulerian description of the two-phase flow is investigated,
the particle being larger than the local space step of the mesh. The originality of the numerical methods is based
on a one-fluid formulation of the incompressible Navier Stokes equations in which the presence of the particles is
tackle with a penalty method to ensure the solid behavior (Randrianarivelo 2005) and (Randrianarivelo 2007). The
particle tracking is obtained via an original hybrid Eulerian-Lagrangian Volume Of Fluid (VOF) method, previously
presented by Randrianarivelo (2005) and Corre (2010). Particles are managed in a Lagrangian manner whereas the
local solid fraction is described in an Eulerian way. In addition, collisions between particles are treated via a soft
sphere model of Cundall (1979). At each time step, the particle characteristics are stored and statistics are performed
on the particle-particle and fluid-particle interactions. A structured Cartesian grid is used and the cylindrical tank
in which the fluidized bed occurs is represented by a penalty method. The simulated cylinder tank corresponds to
the apparatus used in the experimental study of Aguilar Corona (2008) but at a half scale. Numerical results are
compared with available experimental results for validation of the numerical modeling. Instantaneous and average
flow characteristics are reported, as well as the pair distribution function, the streamlines of the mean particle velocity
field, the lagrangian particle auto-correlation coefficients as well as the auto-diffusion coefficents. The results agree
with the literature, notably concerning the maximum of the pair distribution function, good trends are also observed
concerning the order of magnitude of the lagrangian auto-correlation coefficients, the variance of particle movement
and the auto-diffusion coefficients.


Introduction excellent heat and mass transfer characteristics. Al-
though rather simple in its conception, the application
Particles and processes involving particles are of of a fluidized process still faces some challenges: a
paramount importance in the chemical and applied sound understanding of the mechanisms governing the
industries. Fluidized beds are widely employed in complex flow phenomena involved in a fluidized bed
industrial operations, ranging from the pharmaceutical still remains an open technical and scientific issue.
and food industry, to processes such as catalytic crack-
ing of petroleum, combustion and biomass gasification. Several modeling issues can be investigated to
The success of fluidized bed processes is due to their simulate dense particle flows as fluidized beds. On a











macroscopic point of view, Eulerian-Eulerian models
based on the kinetic theory of granular flows KTGF can
be used to simulate the large scales of the flow, the par-
ticle being assumed to be of small size compared to the
characteristic scales of the flow. They are represented
through a diffuse volume fraction. The work of Simonin
(1990) in this field has allowed to recover good results
concerning granular pressure and concentration waves.
The main drawback of these models is associated to the
requirement of defining a priori macroscopic particle
parameters such as the drag force. The models have
been mainly used in two-dimensions. Under the same
scale separation assumption concerning the size of
the particles, Eulerian-Lagrangian models also exist
(Boivin 1998) in which the particle are considered as
force points transported by flow. One or two-way cou-
pling model exists in order to account of the interacting
between particles and the surrounding motion. This
type of modeling is generally applied to gas particle
fluidized beds. As for Eulerian-Eulerian approaches,
it requires the a priori definition of interaction terms
concerning the drag force or particle collisions.


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Numerical modeling

Similarly to existing numerical models allowing to sim-
ulate particle flows at the particle scale (ten Cate 2004)
(Glowinski 2001), our approach is based on a fixed con-
trol volume kernel and Eulerian flow characteristics, in
order to avoid the difficulty of adapting meshes at each
time step when particles move (Maury 1999). Our goal
is to propose an original two-way coupling model in
which the particles are submitted to the flow motion
through Eulerian particle tracking models and at the
same time, the effect of particle collisions with walls or
neighboring particles is directly taken into account in the
momentum conservation equations through a particular
force. These last effects are taken into account explic-
itly with a time splitting modeling by Pan (2002) and
Derksen (2007), which induces a specific modeling er-
ror when dense particle flows are under consideration.

Governing equations

The modeling of incompressible two-phase flows is clas-
sically achieved on fixed elementary control volume ker-
nels by means of a 1-fluid formulation of the mass and
momentum conservation equations (Kataoka 1986) as
follows:


The constant increase of the power brought by super-
computers and the development of new numerical meth-
ods have allowed to carry out in rare research works
direct numerical simulation DNS of real fluidized beds
with thousands of particles. To the acknowledge of the
authors, only three publications exists in the field of
DNS of mono-dispersed particles fluidized by a liquid,
all being investigated on Eulerian fixed grids: the Dis-
tributed Lagrangian method DLM of Pan (2002) which
is validated and compared to dedicated experiments,
the lattice Boltzmann method LBM of Derksen (2007)
which is used to simulate wave instabilities in fluidized
beds and the present work which uses fictitious domains
approaches and pl'c.llii nilli', d, (Randrianarivelo 2007)
(Corre 2010). The interest of these unsteady and small
scale simulations, i.e. the grid space scale is smaller than
the particle diameter, is to enable to model both the fluid
and particle motions in a coupled way without consider-
ing any assumption on the particle size nor on any space
scale separation. The domain of interest of the DNS is
the particle motions in which the size of the particle lies
in the range of the small flow scale, typically the Kol-
mogorov scale for turbulent flows, and the integral scale.
The present work aims at presenting for the first time
the analysis of a real DNS of a fluidized bed in terms
of instantaneous local flow motion, statistical variables
required in KTGF based models, time spectra of particle
velocities or solid fractions and averaged quantities.


P u + u Vu)
a t


-VP + V 7 + Fp


V u 0 (2)
where u is the velocity, P the pressure, t the time, p
the density and 7 p= (Vu + VTu) the viscous stress
tensor for a Newtonian fluid, p being the dynamic vis-
cosity. The source term Fp is due to the implicita-
tion of the particle forces involved by particle-particle
or particle-wall collisions (two-way coupling). The one
fluid model (1-2) is built from a convolution and a fil-
tering of the Navier-Stokes equations in each phase in
which the jump equations valid at the interface (velocity
and stress tensor continuity across fluid-solid interfaces
Delhaye (1974)) are used to close the specific terms due
to space filtering and summation over phases. The lo-
cation of fluid and solid phases is achieved thanks to a
local solid fraction C which is equal to 1 in the parti-
cles and 0 elsewhere. The particle fluid interface is de-
fined as C 0.5. The motion of the particles under the
surrounding flow constraints is described by an Eulerian
material advection equation on C:
aC
+ u VC 0 (3)

The one-fluid model (1-3) is interesting as it applies in
the whole calculation domain, the solid phase being con-
sidered as a fluid with specific physical properties. The











local characteristics of the equivalent fluid, valid in both
the fluid and solid phases, are build with the following
mixing laws:

P Cp + (1 C)pi (4)
t = ts + (1 C) t (5)
where indexes s and I refer to the properties of solid and
fluid media respectively.

Fictitious domain approach and
penalty method

As explained in the previous section, the presence of
particles in the two-phase flow is modeled by a ficti-
tious domain method, i.e. the mesh does not conform
to the spherical particle shape, the motion of all phases
being described by an Eulerian one-fluid approach. Fic-
titious domain approaches were initially introduced by
Caltagirone (1994) and Khadra (2000) to represent fixed
obstacles interacting with incompressible flows through
penalty terms. These penalty methods consist in the in-
troduction of specific physical control terms in the mo-
mentum equations as follows

p ( +u-Vu + u+ B (u u)
Sat ) K
-VP + V 7 + Fp (6)

where Ku is a Darcy term imposing a solid behavior
when the permeability K tends to 0. In the same way,
the volume penalty term B (u u,) imposes the refer-
ence velocity uo when the components of the tensorial
penalty parameter B tends to infinity. The Darcy and
volume penalty approaches have been demonstrated to
be easy to implemented and are particularly adapted
to the resolution of conservation equations on fixed
Eulerian grid. The extension of penalty methods to
moving objects whose motion is known a priori was
first realized with success by Ritz (1999) concerning
the two-dimensional sedimentation of particles in tanks.
The previously cited Darcy or viscosity penalty method
are of first order convergence in space as they consider
the Eulerian volume projection C of the real object
shape on the structured grid to build K or B according
to the phase function C.

An original method was proposed by Caltagirone
(2001) for handling particle flows with a higher order
of convergence in space, based on the one-fluid set of
equations, the carrier fluid and the solid particles still
being located by means of a color function C. The
idea is to split the stress tensor in the Navier-Stokes
equations into four parts, associated to incompressibil-
ity, elongation, pure shearing and rotation respectively


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


(Caltagirone 2001). Thus, the viscous stress tensor can
be written as:

= (-P+AV.u)Id+KF(u)+CO(u) -q (u) (7)

Taking, A K 2, C 2p and r = 2p, the
standard viscous stress tensor for a Newtonian fluid is re-
covered. Thus, in the fluid zones, the artificial viscosity
A will be chosen such as the term AV u will be pre-
dominant in order to impose incompressibility through
an augmented Lagrangian approach. In this case, aug-
mented Lagrangian parameter A is automatically esti-
mated by means of the algebraic approach proposed by
Vincent (2010). In the solid zones, the different artifi-
cial viscosities K, ( and r are chosen in order to impose
a solid behavior. The term KF(u) + (C(u) q?(u)
is made predominant in the Navier-Stokes equations for
the cells belonging to a particle by choosing K -- +oo,
( +00 and r - +oo if C > 0.5. Consequently, with
these different choices of the artificial viscosities, the in-
compressibility of the flow is guaranteed in the fluid and
the non deformation of particles is also maintained with
a second order.

Tracking of particles


Tracking particles is obtained via a hybrid Eulerian-
Lagrangian method (Ritz 1999) (Randrianarivelo 2005).
This numerical algorithm corresponds to a first step in
the two-way coupling between the flow motion and the
particle dynamics: on the one hand, the fluid dynamics
impose the solid character with the penalty method and
provide the particle velocities, and on the other hand,
the movements of the particles modify the local values
of the solid fraction and so involve an unsteady and
local evolution of the local viscosities and densities
used in the one-fluid model.

The particle are characterized by two complementary
quantities: a color function C which is used on the Eu-
lerian grid to locate the solid medium and furthermore
the coordinates of the particle mass center which are uti-
lized to locate the particles in a Lagrangian manner in
order to maintain at each time step their spherical shape.
The color function C is defined over the computational
domain as being the solid fraction, allowing to build ps
and ps in each cell of the Eulerian grid by using nu-
merical two-phase mixture laws presented in the pre-
vious section. Knowing the characteristics of the two
phase medium, the one-fluid Navier-Stokes equations
are solved and velocity vectors are calculated in each
cell. From the Eulerian local velocity field so obtained
in both fluid and solid media according to the penalty
methods, each particle velocity is built by averaging the







7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


velocities in the cells belonging to the particle, i.e. cells
verifying C > 0.5. Finally, each particle can be ad-
vected by a standard Lagrangian advection equation as
follows:


Xn+1 = Xn + V,At


where X' is the position vector of particle i at the
instant nAt and V, is the corresponding translation
velocity vector of particle i. In this way, the new
position of each particle is known. Then, the an-
alytical spherical shape of the particles is projected
onto the Eulerian mesh to build C in a statistical manner.


Let M be the center of a given grid cell. If, for each
particle i, the following relation is verified


. LM >R+ Ax2 Ay2 + Az2
M 2


then this cell of center M is full of fluid and so the color
function C in this cell is equal to 0. Similarly, if, for
each particle i, the following relation is verified


X. 1M11 | 2M l -j (10)

then the Eulerian cell is full of solid and so the color
function C in this cell is equal to 1. The previous re-
lations implicitly use the fact that the particles have a
spherical shape of radius R. This shape is consequently
perfectly maintained at each time step without defor-
mation, even if the global one-fluid modeling is Eule-
rian, which is a very important property when hundred
of thousand time steps are simulated to obtain average
quantities and representative physical times. It is as-
sumed that Ax, Ay and Az are the grid space steps in
each direction respectively. Concerning the grid cells
which contain solid and fluid parts, i.e. the fluid-solid
interface intersects the Eulerian cell, the calculation of
C is obtained on a statistical point of view by testing
the belonging of some points, distributed regularly in the
grid cell, to the interior or the exterior of the particles, as
is represented in figure 1, and by locating the points in-
side or outside the particle. In this figure, the point M
corresponds to the center of the plotted control volume
at which the solid fraction C is estimated. Finally, C in
a given fluid-solid mixed cell is equal to the number of
points inside the particle divided by the total number of
test points. In this way, the solid fraction is known over
the entire computational domain and the time can be in-
cremented to the next iteration. In the calculations, ten
test points by direction are used to evaluate C.


M


Figure 1: Definition sketch in two dimensions of the
seeding of a grid cell by test points in order to
evaluate the local solid fraction C. The black
line represents the fluid-solid interface. The
point M corresponds to the center of the con-
trol volume


Collision model


To complete the two-way coupling modeling, an explicit
model, called collision model, devoted to taking into
account interactions between solid particles, has to be
implemented as soon as the lubricating length and time
scales are not resolved explicitly, i.e. the space grid cell
and the numerical time step are larger than the typical
scales of the lubrication process.

In order to avoid overlapping of the Eulerian zones
resulting from the projection of the particles, a collision
model is implemented. To maintain a reasonable time
step, a soft sphere model is adopted. This model is ap-
plied via an implicit force, called Fp, which is similar to
a restoring force for a spring. Indeed, this force between
two particles i and j is expressed by:


Fp,ij = -Kp a x i (11)

where a is the sum of the artificial deformations due to
numerical overlapping when particles are advected, as
represented in figure 2, and Kp is a constant of stiffness
whose value has to be tuned numerically.

Particles being non deformable, avoiding the particle
overlapping requires to define virtual particles for han-
dling the collision force model, as represented in figure
3. In order to reproduce correctly the rebound between
two particles, some numerical tests have shown that the


* 6 6 6 6 6 6 6 6
* S 0 S 0 0 S 0 S
* 6 0 S 0 6 S 0 S
* 0 0 0 0 0 0 0 0
* S 6 S 0 6 S 0 S
* S 6 S 0 6 S 0 S
S
* S 0 S 0 S 0
* S 0 6 S 6 S 6 S
* 6 0 6 S 6 6 6 S
* S 0 0 S 0 0I S 0 0 S
* S 0 6 S 0 6 S 0 6 S







7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Figure 2: Apparition of an overlapping region between
to particles after the Lagrangian advection
step definition of coefficient a used in the
soft sphere collision model.


selection of a virtual radius equal to R + Ax and a con-
stant of stiffness equal to 2.5 kg.s 2 give good results.
Consequently, the zone of deformation a is equal to:


Figure 3: Definition of virtual spheres for handling the
collision force model the real particle shapes
are plotted with a black line whereas their vir-
tual shape is denoted by a red dotted line.


a -XXj XI 2 x (R + Ax)


The collision force is activated only when a is negative.
The key point in the two-way coupling algorithm is
associated to a correct introduction of this force in the
conservation equations.

At the end of each time iteration, the velocity V' of
each particle is interpolated according to u"+1. This
tentative Lagrangian velocity is used to predict at the
next time iteration the likely position of each particle
by carrying out a tentative advection to obtain tentative
particle positions X'. Then, tests on the possible col-
lisions between every pairs of particle are realized with
theses likely positions. If a collision is detected, then
the collision force Fp,ij is calculated and added in the
right hand of the one-fluid Navier-Stokes equations for
all the cells belonging to the particles i and j concerned
by the collision. For each pair of particles, the force is
different because the value of a depends on the magni-
tude of the overlapping during each collision. Finally,
the two-way coupling is realized in an implicit manner
by adding Fp Fp,ij directly in the momentum con-
servation equations. The velocity field is calculated and
converged at each time step by taking into account the
collision forces as well as the inertia, viscous dissipation
and divergence free constraints at the same time. Such
an implicit model has never been proposed before in the
literature. All the authors (Pan 2002) (Derksen 2007)
advect the particles according to an interpolation of the
Eulerian velocity field and modify afterwards the parti-
cle position by adding a Lagrangian force to the trajec-
tory equation of the particles. The position X' are not
conserved to track the evolutions of the particles over


time These tentative positions are only used to build the
collision force model.


Discretization and solvers


The one-fluid model generalized for modeling dense
particle flows has been implemented in the computa-
tional fluid library CFD named Th6tis, which is devel-
oped at the Trefle laboratory, in collaboration with IMFT
laboratory. Orthogonal structured staggered grids (Har-
low 1965) are used to discretized physical space. The
decoupling between pressure and velocity is fulfilled
via an augmented Lagrangian approach (Fortin 1982),
adapted for handling multiphase flows (Vincent 2004)
(Vincent 2007). The discretization is carried out in an
implicit way with finite volumes. This Eulerian nu-
merical approach is commonly used for its efficiency in
tracking particulate flows and for its easy programming.
The space derivatives of the Navier-Stokes equations are
discretized with a second order centered scheme and
time derivatives are approximated by an Euler scheme
of first order. Concerning the linear algebra, the non-
symmetric discretization matrix resulting from the ap-
proximation of the momentum equations is inverted with
an iterative Bi-CGStab II solver, preconditionned under
a Modified and Incomplete LU method. All these nu-
merical techniques have still been published and vali-
dated in previous publications concerning particle flows.







7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Direct Numerical Simulation of a fluidized bed

Presentation of the simulated fluidized
bed

In the experiments of Aguilar Corona (2008) the flu-
idized bed consists of a cylinder of radius R,6c 4cm
and height of 64 cm filled with 2133 solid particles of
Pyrex, with a density p, equal to 2230 kg.m 3. The
particle diameter is Dp 6mm. The fluid is a solution
of potassium dliii L.\ i.il KSCN with a density pi equal
to 1400 kg.m-3 and a kinematic viscosity / equal to
3.8 10-3 Pa.s. In the present study, in order to save
CPU time, we have performed the simulations on a half
scale configuration but kept the same particle diameter
and the same initial solid fraction. Thus, the flow occurs
in a cylinder of radius R,6y 2 cm and of height
32 cm. A fluidization velocity equal to 0.12 m.s 1 is
imposed. Initially, particles are randomly distributed
into the cylinder on a height of 24 cm. For a particle
diameter of DP 6mm the terminal velocity is
equal to 0.226 m.s 1 and the particle Reynolds number
based on this velocity is equal to 530. For this case, a
parallelepiped box of dimensions 4 cm x 4 cm x 32 cm
is used with a grid mesh containing 75 x 75 x 800 cells.
With such grid resolution, nearly 12 Ax per particle
diameter is achieved. The cylinder is represented
numerically with a Darcy penalty method by taking a
permeability tending to infinity in the cells being outside
the cylinder. The rebound of particles on the walls of
the cylindrical tank is treated as a rebound between two
particles by inserting a virtual size of the cylindrical
tank whose radius is R,6y Ax, as is done for the
collision model presented previously. The simulation
are realized on 64 processors during a total physical
time of 42s with a numerical time step of At = 0.001s.

A global view of the particles inside the fluidized bed
at a given time is provided in figure 4. Six horizon-
tal planes are selected along the vertical axis in order
to carry out flow analysis in the next sections. These
slices, called S1 to S6, correspond to vertical coordi-
nates z = 0.04 m, z = 0.09 m and z = 0.14 m inside
the bed and z = 0.19 m, z = 0.24 m and z = 0.29 m
outside of it.


Macroscopic properties of the fluidized
bed
Figure 4: View of an instantaneous position of the par-
In the figure 5, the time evolution of the bed height and tides inside the fluidized bed six horizontal
the solid fraction are presented. It is observed that af- planes are defined for future analysis.
ter 5 s, a stabilized particle flow behavior is reached.
All the averaging and statistical analysis provided in the
next sections will be carried out after this stabilization




















021 -

0.2 -



-
g017
m
0.16

0.15

0.14 L -L
5 10 15 20 25 30
Time (s)
0 16










0.1504 -
0 014 -
"o
V 013









0.04



S004-


C 303
0 0 0 2 A

0


-0.02


5 10 15
Time (s)


20 25


Figure 5: Time evolution of bed height, averaged solid
fraction and fluctuations of averaged solid
fraction.


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


time. The averaged solid fraction ap is equal to 0.135
and the corresponding bed height is hb 0.18 m. The
solid fraction is deduced from the bed height by using
the following equation:

N,4/37R3
ap(t) = (13)
wR91.. (t)

where Np is the total number of particles and Rc the
radius of the cylinder. Concerning the solid fraction
fluctuations, they are oscillating between 1'. of the
average solid fraction. If we try to compare the aver-
aged solid fraction obtained from the simulation to the
Richardson-Zaki correlation (Richardson 1954) for the
same fluidization velocity given by


p 1 ( )1/2.41 (14)
Ut
where Uo is the fluidization velocity and
Ut 0.226m/s is the terminal velocity, we find a
discrepancy of nearly I, on the averaged solid frac-
tion. This calculation of ap for the simulation have been
achieved using equation (13) with a value of the particle
radius of 6mm. However, the effective numerical radius
due to the collision model is Re Rp + Ax. By
using this effective radius, an averaged solid fraction of
ap 0.21 is found, to be compared to Richardson-Zaki
value which is aprz 0.23. A good agreement is now
obtained with the reference correlation (Richardson
1954) and also with Aguilar's experiments (Aguilar
Corona 2008).

The fluid properties are analyzed in slice planes
S1 to S6 in figures 6, 7 and 8. It is observed that
the maximum fluid velocity inside the bed (slices
S1 to S3) is located near the cylinder wall, with a
value of 0.2m.s 1 which is I;'. larger than the mean
fluidization velocity, while the maximum fluctuations
of the fluid velocities are measured inside the bed
due to particle motion and collisions. These fluctu-
ations are of the order of the fluidization velocity,
reflecting the fact that particles induce blocking flow
effects and strong accelerations alternatively. To finish
with, the averaged fluid fraction demonstrates that
the fluid flow path is mainly located near the cylinder
walls, with localized cylinder fluid pathes inside the bed.

Another macroscopic behavior of the fluidized bed
is described in figure 9 concerning the stream lines of
the mean particle velocity in the bottom of the bed, for
z < 0.15m. This mean particle velocity field has been
obtained by first defining a new Eulerian grid, much
coarser than the fluid simulation mesh. On this new grid
called Gp, a time averaging has been done by checking,
for each cell of Gp, the presence of a solid particle and








7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Alphe flud

095

06
045
o 55

oos
035


Vz
C1
C18
C1 4
C 02
C12
C1
08
C 08
C02


Figure 6: View of time averaged vertical fluid velocities
in slice planes S1, S2 and S3 from top to bot-
tom in the left column and planes S4, S5 and
S6 from top to bottom in the right column.


Figure 7: View of vertical fluid velocity fluctuations in
slice planes S1, S2 and S3 from top to bottom
in the left column and planes S4, S5 and S6
from top to bottom in the right column.


Figure 8: View of averaged fluid fractions in slice
planes S1, S2 and S3 from top to bottom in
the left column and planes S4, S5 and S6 from
top to bottom in the right column.



by storing the particle velocity if a particle lies inside
the related cell of the new Eulerian grid. Time averag-
ing has been performed between t = 5s and t = 40s.
Moreover, axisymetry of the mean field has been sup-
posed. It is demonstrated that recirculation loops ex-
ist in the first 0.05m of the vertical cylinder, the parti-
cles going upward in the near wall zones and coming
down in the central zone. This toroidal motion of the
particles has still been observed in the experiments of
Aguilar Corona (2008), Handley (1966), Carlos (1968),
Latif (1972), where the characteristic size of the recircu-
lation loop in the inferior part of the fluidized bed was
measured to be of the order of the cylinder radius, as in
our simulations.



Particle trajectories

The trajectories of three different particles are plotted in
figure 10 with a top view and in figure 11 with a side
view in half a cylinder diameter. It is observed that dur-
ing the same time, some particles travel in all the cylin-
der section and across its whole height while other parti-
cle trajectories are restricted to local horizontal and ver-
tical bed zones. Similar behaviors have been noticed by
Aguilar Corona (2008) in her measurements. On a statis-
tical point of view, the motion of the ensemble of spheres
occupies the whole cylinder volume, in the zone where
the fluidized bed is active, i.e. where z < hb.








7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


0.04


0.03 -


>, 0.02


0.01 -


Figure 9: Streamlines of the mean particle velocity


Pair distribution function g(r)

This is the conditional probability density of finding a
particle at a distance r, given that there is a particle at a
chosen origin coordinate r = 0. Thus, g(r) provides a
measure of local spatial ordering in a particulate flow. It
is calculated numerically as follows:

NdV, Vtot (15)
g (r) (15)
N, dV N,

where NdV, is the number of particles belonging to
the volume delimited by r = 0, i.e. the center of the
reference particle, and r = R + Ax + dr, dr is the
tentative distance used for building the probability, Np
is the total number of particles in the fluidized bed,
Vto, is the volume of the cylinder which contains the
particles and dV, is the volume of the spherical shape
associated to the tentative radius dr included in the flow
cylinder.

As can be seen in figure 12, the maximum of g(r) is
obtained for a distance r equal to 2.5 Rp, which corre-
sponds to two times the radius of the virtual spheres, in
other words, the distance for which the particles touch.
The value of the maximum is important because this
value is used in the statistical modelings to evaluate col-
lision terms. This value, named go, depends on the solid
fraction, as provided by the correlation of Lun (1986):


go (1


a -2.5a_
max 7


where amax = 0.64 is the solid fraction due to packing
effects for spherical particles. In our numerical case, for


0.03


>, 0.02



0.01


0.03


>, 0.02



0.01



0


) 0.01 0.02 0.03 0.(
x



















) 0.01 0.02 0.03 0.<
X
I I I I










S0.01 0.02 0.03 0.(
X
I I I


-


-/











o 0.01 0.02 0.03 0.
x


Figure 10: Examples of time evolution of horizontal co-
ordinates of three different particles trajec-
tories are observed from the top of the flu-
idized bed.


S' I I I







7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


0.2 ----


0.2 -




0.15-




N 0.1-




0.05-




0 00
0 0.02


Figure 13: Lagrangian particle auto-correlation func-
tion estimated in a centered cylinder of ra-
dius 2 cm.


0.02


Figure 11: Examples of time evolution of vertical co-
ordinates of three different particles trajec-
tories are observed from the side of the flu-
idized bed in half a cylinder section.



a 0.21, the value of go is 2.12 whereas Lun (1986)
give a value of 1.88 for the same solid fraction.


r /(R+Ax)


Figure 12: Pair distribution function g(r) according to
dimensionless radius r/(R + Ax).




Lagrangian particle auto-correlation
function and auto-diffusion coefficient

For given time to, the Lagrangian particle auto-
correlation function Rt, (r) is expressed as


Ui!/n( ) (17)
( (to)) u, (to + T)

where i y, z is the considered direction, Pn
is the number of the ,' particle, the bracket operator
(.) is an averaging operator over all the particles, and
u'pn,' (to), is the particle velocity fluctuation defined as
follows


u'prni (to) = Upn,, (to) (upni (to))


and u'J,, (to) is the square of the particle velocity
fluctuation given by

U i (to) p ,i (to) (u p,i (to)) (19)

Finally, the Lagrangian auto-correlation function Rii (r)
is obtained as the average of Ri' (7) over all the times to.

The evolution of the auto-correlation functions Ri (r)
is plotted in figure 13 for the three space directions. It is
observed that the magnitude of the Rii(r) is greater in
the vertical direction than in the horizontal directions,
indicating large scale movements in the axial direction.
Moreover, the values of Rx and Ry are superimposed,
as in the experiments of Aguilar Corona (2008).

From those auto-correlation functions, the auto-
diffusion coefficients Dii can be calculated in the three
directions. They are of importance in the interpreta-
tion and characterization of the mixing in fluidized beds.
They are defined as the product between the time and


Time (en s)


0.15-




N 0.1-




0.05-


0 --0.0
0 0.02







7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Simulation Experiment
T W (s) 0.0357266 0.0694
Ty, (s) 0.0346635 0.0659
T. (s) 0.1076160 0.3163
(u4'~ (m'.s -) 7.56e-4 7.78e 4
(u' 2 ) (m2.s2) 7.79e-4 7.74e-4
(u'/,z) (m2.s 2) 1.30e-3 2.15e 3
D11 (m2.s 1) 2.7e-s 5.4e
Dyy (m2.s 1) 2.7e-5 5.1e-5
Dz; (m2.s 1) 1.4e 4 6.8e 4

Table 1: Auto-diffusion coefficients in each space direc-
tion for the simulation of the reduced scale flu-
idized bed and corresponding experiments of
Aguilar Corona (2008).


space average of the square of the particle velocity fluc-
tuation times the time integral of the auto-correlation
function as follows :

Di ( u'/ Th (20)

where the Lagrangian time scale of the particles is given
by

Tii = Ri(T)d(7) (21)

The comparisons of numerical auto-diffusion co-
efficients, given in table 1, with experimental values
show that a correct order of magnitude is obtained
concerning the Lagrangian time scales, the variance of
the particle fluctuation velocity and the auto-diffusion
coefficients. However, the simulation clearly underesti-
mate the axial and vertical auto-diffusion coefficients.
From table 1, it can be deduced that these differences
between the auto-diffusion coefficients are first due
to the estimate of the Lagrangian time scales Tu, for
which the simulations are two to three times lower
than the experimental values. The gap between nu-
merical results and measurements are also caused by
the underestimate of the particle velocity fluctuations
in the vertical direction. It should be recalled that the
simulations are performed on a half scale fluidized bed
which could explain the observed differences. Direct
numerical simulations of the real fluidized bed are
under investigation on a 100 x 100 x 800 grid. These
simulations are not completely converged concern-
ing statistical quantities such as the auto-correlation
functions. However, the treatment of available results
provide values of D, Dy, 6e5 m2 s 1 and
Dzz le-3 m2 s which are closer to expected
experimental measurements.

The auto-diffusion coefficients can also be defined ac-
cording to the variance Mi, of the particle movement in


M,-- along xand y-radius =1 cm
,-e06- Slope along adius = 1cm
along x andy-radius 2 cm
along- rdius 2cm




0,001 0.01 OA 1
Time (s)

Figure 14: Variance of the particle movement with time
in centered cylinders of radius 1 and 2 cm
in the three space directions to = 5s a
lop-log scale is used.


each direction. This variance is the time average of the
particle relative position in direction i at time t compared
to an initial time to. It reads :

Mii ((xi(t) Xi(to)]2) (22)

with to 5 s is the time for which the height of the bed
is stabilized. The diffusion coefficient in direction i is
then
Dii -1 Mi(t to) (23)
2 dt
The values obtained for the Di calculated in two cylin-
ders of radii 1 cm and 2 cm, centered in the fluidized
bed, are provided in figure 14 in log-log scale. It is ob-
served that the x and y direction are superimposed and
that the variance in the vertical direction is larger than in
the other directions. In addition, the slope of all Di is
equal to 2 until time t to 1 s. These observations
were also reported by Aguilar Corona (2008) in her ex-
periments. The saturation of the variance for the largest
times t to is classically obtained due to confinement.
From the time evolution of the variance, it can be de-
2M, (t)
duced that Mu(t) at2 and M'(t) 2at 2M t
In this way, the auto-diffusion coefficient can be refor-
mulated as

M (t to)
Dii = limt-to-+oo o) (24)

The use of equation (24) to estimate the auto-diffusion
coefficient requires the estimation of the limit val-
ues of Mi(t-to) In each direction, they are re-
(t to)
spectively equal to 4e 5 m2.s 1, 4e 5 m2.s 1 and
2.5e-4 2. -1. It is observed that DDn and D22 are
".i larger than the coefficients estimated thanks to the











Lagrangian particle auto-correlation function Ri while
D33 is twice larger than the previously calculated value.
No clear explanation can be given on this result.


Conclusions

The direct numerical simulation of a liquid-solid
fluidized bed has been investigated with a one-fluid
model, a fictitious domain approach, a mixed Eulerian-
Lagrangian VOF technique for tracking particle
evolutions and an augmented Lagrangian method for
handling the pressure-velocity coupling in the fluid. The
solid behavior is ensured in the particles by a penalty
approach based on a decomposition of the viscous
stress tensor and the collision between the particles is
explicitly accounted for with a soft sphere model.

Based on the parameters of the experiments of
Aguilar Corona (2008), the simulations were carried
out on a fluidized bed whose dimensions were half
the real scales while the size of the particles were
keep unchanged. A reduced calculation time cost
were so involved. The macroscopic characteristics
measured numerically have been favorably compared
to experiments and existing laws of the literature in
terms of fluidization velocity, solid concentration after
a stabilized state is reached or toroidal recirculation
loop at the inferior part of the bed. Concerning particle
trajectories, it has been demonstrated that the whole bed
volume is statically concerned by the particle motion.

On a statistical point of view, the pair distribution
function has been found equal to 2.12, which is in
good agreement with the value of 1.88 given by the
law of Lun (1986). Concerning Lagrangian particle
auto-correlation functions and auto-diffusion coef-
ficients, they have been found to be in the order of
magnitude than that measured by Aguilar Corona
(2008), the numerical values being 2 two lower than in
the measurements. The auto-diffusion coefficients have
been observed to be equal while the time evolution of
the variance of the particle movement has been found to
grow as t2, as expected.

The simulation of the fluidized bed of Aguilar Corona
(2008), with the real dimensions, is under investigation
on 256 processors with a time step of 0.001 s. Obtaining
converged statistical quantities such as the pair distri-
bution or the auto-correlation functions require to sim-
ulate a physical time of 40 s, corresponding to 40000
iterations. This point is one of the main remaining dif-
ficulties, even on HPC computers, for leading real DNS
comparable to "numerical experiments". Concerning the
analysis of the particle flow motion, the spectral power


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


of the fluctuations of the particle concentrations and par-
ticle velocities will be extracted in order to estimate the
inertial character of the fluidized bed. To finish with, the
effect of the collision model on the behavior of the bed
will have to be studied and solid fraction superior to 0.21
will also be considered.


Acknowledgements

We acknowledge the Aquitaine Regional Council for the
financial support dedicated to a 256-processor cluster in-
vestment, located in the TREFLE laboratory. This work
was also granted access to the HPC resources of CCRT,
CINES and IDRIS by GENCI (Grand Equipement
National de Calcul Intensif) under reference number
x2009026115 and x2010026115. Part of C. Corre PhD
thesis has been funded by a Midi-Pyrenees Regional
Council grant.


References

Aguilar Corona A., Fluctuations of particles in a flu-
idized bed. Experimental study, PhD thesis, 2008, INPT.

Boivin M., Simonin O. and Squires K. D., Direct numer-
ical simulation of turbulence modulation by particles in
isotropic turbulence, J. Fluid Mech., vol. 375, pp. 235-
263, 1998.

Caltagirone J.-P., On the fluid-porous interaction: appli-
cation to the calculation of forces exerted by fluids on an
obstacle, C. R. Acad. Sci, S6rie IIb, vol. 318, pp. 571-
577, 1994.

Caltagirone J.-P and Vincent S., A tensorial penalty
method for solving the Navier-Stokes equations, C.-R.
Acad. Sci., s6rie IIb, vol. 329, pp. 607-613, 2001.

Carlos C.R. and Richardson J.F., Solids movements
in liquid fluidised beds-I Particle velocity distribution,
Chem. Eng. Sci., vol. 23(8), pp. 813-824, 1968.

Corre C., Estivalezes J.L., Vincent S. and Simonin O.,
Simulation of a fluidized bed using a hybrid Eulerian-
Lagrangian method for particle tracking, accepted in
"Notes on Numerical Fluid Mechanics and Multidisci-
plinary Design", 2010

Cundall PA. and Strack O.D.L., A discrete numeri-
cal model for granular assemblies. G6otechnique, vol.
29(1), pp. 47-65, 1979.

Delhaye J.M., Jump conditions and entropy sources in
two-phase systems. Local instant formulation, Int. J.
Multiphase Flow, vol. 1, pp. 395-409, 1974.











Derksen J.J. and Sundaresan S., Direct numerical simu-
lations of dense suspensions: wave instabilities in liquid-
fluidized beds, J. Fluid Mech., vol. 587, pp. 303-336,
2007.

Fortin M. and Glowinski R., M6thodes de lagrangien
augment. Application a la resolution num6rique de
problems aux limits. Collection "M6thodes math6ma-
tiques de l'informatique", Dunod, Paris, 1982.

Glowinski R., Pan T.W., Hesla T.I., Joseph D.D. and
P6riaux J., A fictitious domain approach to the direct nu-
merical simulation of incompressible viscous flow past
moving rigid bodies: application to particulate flow, J.
Comput. Phys., vol. 169, pp. 363-426, 2001.

Handley D., Doraisamy A., Butcher K.L. and Franklin
N.L., A study of the fluid and particle mechanics in
liquid-fluidised beds, Trans. Inst. Chem. Eng., vol. 44,
pp. 260-273, 1966.

Harlow F.H. and Welch J.E., Numerical calculations
of time dependent viscous incompressible flow of fluid
with a free surface, Phys. Fluids, vol. 8, pp. 2182-2189,
1965

Kataoka I., Local instant formulation of two-phase flow,
Int. J. Multiphase Flows, vol. 12(5), pp. 745-758, 1986.

Khadra K., Angot P., Parneix S. and Caltagirone J.-P.,
Fictious domain approach for numerical modelling of
navier-stokes equations, Int. J. Numer. Meth. Fluids, vol.
34,pp.651-684,2000.

Maury B., Direct simulations of 2D fluid-particle flows
in biperiodic domains, J. Comput. Phys., vol. 156(2), pp.
325-351, 1999.

Latif B.A.J. and Richardson J.F., Circulation patterns
and velocity distribution for particles in a liquid fluidised
bed, Chem. Eng. Sci., vol. 27, pp. 1933-1949, 1972.

Lun C.K.K. and Savage S.B., The effects of an impact
velocity dependent coefficient of restitution on stresses
developed by sheared granular materials, Acta Mechan-
ica, vol. 63(1), pp. 15-44, 1986.

Pan T.-W., Joseph D.D., Bai R., Glowinski R. and Sarin
V., Fluidization of 1204 spheres: simulation and experi-
ment, vol. 451, pp. 169-191, 2002

Qi D., Lattice-Boltzmann simulations of particles in
non-zero-Reynolds-number flows, J. Fluid Mech., vol.
385, pp. 41-62, 1999.

Randrianarivelo T.N., Numerical study of fluid-solid hy-
drodynamic interactions: application to fluidized beds,
PhD thesis, 2005, Bordeaux 1 University.


7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010


Randrianarivelo T.N., Pianet G., Vincent S. and Calta-
girone J.-P, Numerical modelling of the solid particle
motion using a new penalty method, Int. J. Numer. Meth.
Fluids, vol. 47 (10-11), pp. 1245-1251, 2005.

Randrianarivelo T.N., Vincent S., Simonin 0. and Calt-
agirone J.-P, A DNS approach dedicated to the analysis
of fluidized beds, Fluid Mech. Appl., vol. 81, pp. 207-
214, 2007.

Richardson J.F. and Zaki W.N., Sedimentation and flu-
idization. Part 1, Trans. Inst. Chem. Eng., vol. 32, pp.
35-53, 1964.

Ritz J.B. and Caltagirone J.-P., A numerical continuous
model for the hydrodynamics of fluid particle systems,
Int. J. Numer. Meth. Fluids, vol. 30, pp. 1067-1090,
1999.

Simonin 0. and Viollet PL., Modelling of turbulent two-
phase jets loaded with discrete particles, Phase-Interface
Phenomena Mult. Flow, pp. 259-269, 1990.

ten Cate A., Derksen J.J., Portela L.M. and Van den
Akker H.E.A., Fully resolved simulations of colliding
monodisperse spheres in forced isotropic turbulence. J.
Fluid Mech., vol. 519, pp. 233-271, 2004

Vincent S., Caltagirone J.-P, Lubin P and Randrianariv-
elo T.N., An adaptative augmented Lagrangian method
for three-dimensional multi-material flows, Comput.
Fluids, vol. 33, pp. 1273-1289, 2004.

Vincent S., Randrianarivelo T.N., Pianet G. and Calta-
girone J.-P, Local penalty methods for flow interacting
with moving solids at high Reynolds numbers, Comput.
Fluids, vol. 36, pp. 902-913, 2007.

Vincent S., Sarthou A., Caltagirone J.-P, Sonilhac F.,
Fevrier P, Mignot C. and Pianet G., Augmented La-
grangian and penalty methods for the simulation of two-
phase flows interacting with moving solids. Application
to hydroplaning flows interacting with real tire tread pat-
terns, under correction in J. comput. Phys., 2010.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs