7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Direct Numerical Simulation of a LiquidSolid Fluidized Bed
C. Corre*, J.L. Estivalezes*, S. Vincentti and 0. SimoninT
ONERA/DMAE, 2 avenue Edouard Belin, BP 74055, 31055 Toulouse, France
t TREFLEENSCPB, UMR 8508 CNRS, 16 avenue PeyBerland, 33607 Pessac Cedex, France
University de Toulouse ; INPT, UPS ; IMFT ; All6e Camille Soula, F31400 Toulouse, France
SCNRS; TREFLE ; F33407 Pessac, France
CNRS; IMFT ; F31400 Toulouse, France
cedric.corre@onera.fr, JeanLuc.Estivalezes@onera.fr, vincent@enscpb.fr and simonin@imft.fr
Keywords: direct numerical simulation, fluidized bed, EulerianLagrangian 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 twophase 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 onefluid 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 EulerianLagrangian 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 particleparticle and fluidparticle 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 autocorrelation coefficients as well as the autodiffusion 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 autocorrelation coefficients, the variance of particle movement
and the autodiffusion 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, EulerianEulerian 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 twodimensions. Under the same
scale separation assumption concerning the size of
the particles, EulerianLagrangian models also exist
(Boivin 1998) in which the particle are considered as
force points transported by flow. One or twoway 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 EulerianEulerian 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 twoway 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 twophase flows is clas
sically achieved on fixed elementary control volume ker
nels by means of a 1fluid 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 monodispersed 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 particleparticle
or particlewall collisions (twoway coupling). The one
fluid model (12) is built from a convolution and a fil
tering of the NavierStokes equations in each phase in
which the jump equations valid at the interface (velocity
and stress tensor continuity across fluidsolid 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 onefluid model (13) 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 twophase 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 onefluid 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 ( +uVu + 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 twodimensional 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 onefluid 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 NavierStokes
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 NavierStokes 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 twoway 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 onefluid 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 twophase mixture laws presented in the pre
vious section. Knowing the characteristics of the two
phase medium, the onefluid NavierStokes 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 onefluid 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 fluidsolid
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 fluidsolid 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 fluidsolid interface. The
point M corresponds to the center of the con
trol volume
Collision model
To complete the twoway 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 twoway 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 onefluid NavierStokes 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 twoway 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 onefluid 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 NavierStokes 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 BiCGStab 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.m3 and a kinematic viscosity / equal to
3.8 103 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
RichardsonZaki 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 RichardsonZaki
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 autocorrelation 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 autocorrelation
function and autodiffusion 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 autocorrelation function Rii (r)
is obtained as the average of Ri' (7) over all the times to.
The evolution of the autocorrelation 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 autocorrelation 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.56e4 7.78e 4
(u' 2 ) (m2.s2) 7.79e4 7.74e4
(u'/,z) (m2.s 2) 1.30e3 2.15e 3
D11 (m2.s 1) 2.7es 5.4e
Dyy (m2.s 1) 2.7e5 5.1e5
Dz; (m2.s 1) 1.4e 4 6.8e 4
Table 1: Autodiffusion 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 autocorrelation
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 autodiffusion 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 autodiffusion
coefficients. However, the simulation clearly underesti
mate the axial and vertical autodiffusion coefficients.
From table 1, it can be deduced that these differences
between the autodiffusion 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 autocorrelation
functions. However, the treatment of available results
provide values of D, Dy, 6e5 m2 s 1 and
Dzz le3 m2 s which are closer to expected
experimental measurements.
The autodiffusion coefficients can also be defined ac
cording to the variance Mi, of the particle movement in
M, along xand yradius =1 cm
,e06 Slope along adius = 1cm
along x andyradius 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
loplog 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 loglog 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 autodiffusion coefficient can be refor
mulated as
M (t to)
Dii = limtto+oo o) (24)
The use of equation (24) to estimate the autodiffusion
coefficient requires the estimation of the limit val
ues of Mi(tto) In each direction, they are re
(t to)
spectively equal to 4e 5 m2.s 1, 4e 5 m2.s 1 and
2.5e4 2. 1. It is observed that DDn and D22 are
".i larger than the coefficients estimated thanks to the
Lagrangian particle autocorrelation 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 liquidsolid
fluidized bed has been investigated with a onefluid
model, a fictitious domain approach, a mixed Eulerian
Lagrangian VOF technique for tracking particle
evolutions and an augmented Lagrangian method for
handling the pressurevelocity 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
autocorrelation functions and autodiffusion 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 autodiffusion 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 autocorrelation 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 256processor 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 MidiPyrenees 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 fluidporous 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 NavierStokes equations, C.R.
Acad. Sci., s6rie IIb, vol. 329, pp. 607613, 2001.
Carlos C.R. and Richardson J.F., Solids movements
in liquid fluidised bedsI Particle velocity distribution,
Chem. Eng. Sci., vol. 23(8), pp. 813824, 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. 4765, 1979.
Delhaye J.M., Jump conditions and entropy sources in
twophase systems. Local instant formulation, Int. J.
Multiphase Flow, vol. 1, pp. 395409, 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. 303336,
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. 363426, 2001.
Handley D., Doraisamy A., Butcher K.L. and Franklin
N.L., A study of the fluid and particle mechanics in
liquidfluidised beds, Trans. Inst. Chem. Eng., vol. 44,
pp. 260273, 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. 21822189,
1965
Kataoka I., Local instant formulation of twophase flow,
Int. J. Multiphase Flows, vol. 12(5), pp. 745758, 1986.
Khadra K., Angot P., Parneix S. and Caltagirone J.P.,
Fictious domain approach for numerical modelling of
navierstokes equations, Int. J. Numer. Meth. Fluids, vol.
34,pp.651684,2000.
Maury B., Direct simulations of 2D fluidparticle flows
in biperiodic domains, J. Comput. Phys., vol. 156(2), pp.
325351, 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. 19331949, 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. 1544, 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. 169191, 2002
Qi D., LatticeBoltzmann simulations of particles in
nonzeroReynoldsnumber flows, J. Fluid Mech., vol.
385, pp. 4162, 1999.
Randrianarivelo T.N., Numerical study of fluidsolid 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 (1011), pp. 12451251, 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.
3553, 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. 10671090,
1999.
Simonin 0. and Viollet PL., Modelling of turbulent two
phase jets loaded with discrete particles, PhaseInterface
Phenomena Mult. Flow, pp. 259269, 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. 233271, 2004
Vincent S., Caltagirone J.P, Lubin P and Randrianariv
elo T.N., An adaptative augmented Lagrangian method
for threedimensional multimaterial flows, Comput.
Fluids, vol. 33, pp. 12731289, 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. 902913, 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.
