Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Discrete Particle Model for simulating liquidsolid fluidization
Micheline Abbas a, Martin Van der Hoef b and Hans Kuipers b
University de Toulouse, INPT UPS, Laboratoire de Genie Chimique, Toulouse, France
bFaculty of Science and Technology, University of Twente, The Netherlands
Email: Micheline.Abbas@ensiacet.fr
Keywords: Liquidsolid fluidization, granular temperature, Discrete Particle Model
Abstract
We have used the discrete particle model combined with computational fluid dynamics for the continuous phase
(DPM) to simulate the fluidization of 0.165 mm glass particles in water, and compare with the experiments of Zivkovic et al.
(2' "i r The added mass and lubrication forces are included in the model, whereas they are generally neglected in the
simulation of gassolid and also most of the liquidsolid systems. The statistical quantities related to the discrete phase were
calculated using the space and / or time averages (pressure drop, bed height, stress developed in the system ...). The
expansion curve of the fluidized bed, obtained using the DPM, has a similar trend as the experimental one. This expansion
curve depends on the choice of the drag force which is the key element in the twoway coupling between the solid and liquid
phases. The granular temperature of the solid phase is calculated from the fluctuating velocities of the particles. As in the
experiments, the granular temperature is nearly equal to the square root of the superficial liquid velocity, and reaches a ceiling
at high velocities when an equilibrium between the kinetic and collisional momentum transfer is reached. The particle
velocity fluctuations are anisotropic. Small and large scale fluctuations are observed and their influence on the local granular
temperature is shown explicitly. In contrast to the gassolid fluidization, the choice of the restitution coefficient has a minor
influence on the statistical behavior of the liquid fluidized bed. The inclusion in the model of the added mass force and
viscous lubrication forces between colliding particles does not have a significant impact on the macroscopic bed properties.
Introduction
Liquidsolid fluidized beds have many applications such as
mineral processing operations (Mukherjee et al., 2009),
heat exchangers (Pronk et al., 2009), water treatment and
others. Nevertheless, this topic has been investigated much
less intensely than gasfluidization which is involved in
processes that are economically more significant like Fluid
Catalytic Cracking for instance.
Models based on Lagrangian tracking of particles
combined with computational fluid dynamics for the
continuous phase (called DPM in this paper) have become
stateoftheart for simulating gassolid flows, especially in
fluidization processes. This type of models falls in between
the twofluid model (TFM) used for simulations of large
industrial processes (Kuipers et al. (1993), Gobin et al.
(2 11')), and the direct numerical simulations (DNS)
applicable only for small scale systems (Pan et al. 2002).
The TFM model has the disadvantage of relying heavily on
closures based on the kinetic theory of granular media. The
DNS is computationally expensive and is limited to the
small scale systems. The clear advantage of a model based
on the tracking of the particles compared to the Eulerian
models used for liquidsolid fluidized beds (Lettieri et al.,
2006) is that the first method gives information about the
behavior of the individual particles in parallel with the
continuous phases. This allows one to obtain qualitative
and quantitative distributions of the particle properties over
the different bed regions (velocity fluctuation distribution
for example). Moreover particleparticle interactions can
be modeled much more accurately in a discrete particle
approach.
The flows encountered in the gassolid fluidized beds are
usually of agitated nature (bubbly flows, fast fluidization
...). By contrast, most of the liquidsolid fluidized beds are
homogeneous and quiescent systems due to the high
density and viscosity of the liquid compared to the gas.
Hence, for liquidsolid systems, models which do not fully
resolve the fluid flow between the particles is more critical
than for gassolid mixtures, since even slight imperfections
in the model can induce instabilities in the system. Studies
of liquidsolid fluidized beds based on a DPM approach
are scarce. Recently Di Renzo and Di Maio (21" i) have
performed a study on the expansion and stability of
liquidsolid fluidized bed. However their analysis was only
twodimensional.
In this work, we use the DPM as described in Hoomans et
al. (1996) and test its validity for studying a particular type
of liquidsolid fluidization, namely that of 0.165 mm glass
particles in water following the system studied by Zivkovic
et al. (2'" 'r '. We are particularly interested in studying the
solid phase agitation, often called "granular temperature",
induced by the fluidization process. This agitation arises
from the particle velocity fluctuations which are induced
by the particleparticle and particlefluid interactions. It is
important to study the particle velocity fluctuations since
they are closely related to the particle diffusion which was
a focus for many studies since the early sixties (Handley et
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
al. (1966), Carlos and Richardson (1968), Kang et al.
(1990), Aguilar (2'" ')). Note that particle diffusion is a
key property since it is directly related to the efficiency for
physical and chemical phenomena (chemical or biological
reactions, heat exchange, and others).
Note that additional forces, which are usually more
significant for the simulation of liquidsolid than gassolid
mixtures, are included in our DPM model. The
implementation of the added mass and lubrication forces
allows for taking into consideration the respective effects
of the accelerating fluid around the particles and of the
interstitial fluid between two particles. These forces are
generally neglected in the simulation of gassolid and also
most of the unresolved liquidsolid systems. In this study,
we test their influence on the behavior of the particular
fluidized bed considered here, which to our knowledge has
not been performed before, in these type of models.
The paper is organized as following. First a brief
description of the numerical scheme is given, including the
equations of the continuous and discrete phases. Next we
summarize the simulation parameters of the modeled
system. The computation of the macroscopic properties of
the bed is then shown: bed expansion, velocity fluctuation
distribution and granular temperature, where we focus in
particular on the dependence on the superficial fluidization
velocity, and hence on the solid volume fraction. We end
the paper with some concluding remarks.
Nomenclature
a Particle radius (m)
dp Particle diameter (m)
dt Time step (s)
e Restitution coefficient ()
F Force (N)
g Gravitational constant im qi
h Gap lenth (m)
k Spring stiffness (N/m)
m Mass (kg)
Np Particle number
p Pressure field (NIi )
r Position in the fluid field (m)
Re Reynolds number ()
t Time (s)
u Fluid velocity field (m/s)
Uo Superficial fluidization velocity (m/s)
Ut Particle terminal settling velocity (m/s)
v Particle Velocity (m/s)
V Volume (m3)
Greek letters
p Interphase drag coefficient ()
e Volume fraction ()
p Density (kg/m3)
u Viscosity (Pa.s)
1 Friction coefficient
Subscripts
a Individual particle
cut Cutoff barrier
D Drag
f Fluid phase
G Gravity
Lubrication
Normal direction
Pressure
Particle phase
Particleparticle
Virtual mass
Tangential direction
Numerical Scheme
Fluid equations
The flow field of the liquid phase (pressure p and flow
velocity u) is computed from the volumeaveraged
NavierStokes equations given by
Pf + PfV (u) = 0
at
f(Efu)+ (f)
Pf a +pfV (efuu)
at
eVp V (eT,) Sp +fpfg
The viscous stress tensor ; is assumed to obey the general
form for a Newtonian fluid (Bird et al. 1960). p is the fluid
pressure field. The effect of the excluded volume of the
particles appears in the fluid equations through the fluid
volume fraction Ef. It is calculated using the local volume
occupied by the particles in the different fluid cells (see eq.
(6)). Twoway coupling is achieved via the momentum
sink/source term S, which includes the fluidparticle drag
force and the added mass force contributions.
The source term S is defined as
1
S v= cell
Ven acen
Va u a)FvM,a 8(rra)
1Ef j J
where Vce,, represents the volume of a computational cell,
Va and v, the volume and the linear velocity of particle a
respectively. The function ensures that the drag force acts
as a point force at the central position of this particle. The
source term contains, according to the principle of
action/reaction, the opposite of the drag force and of the
virtual mass force Fvm,a (defined further in eq. (10)), which
are both applied by the fluid on particle a.
Fixed Eulerian mesh is used to solve the NavierStokes
equations of the fluid. Eq (1) is valid provided that the
particle volume fraction is less than 1 (sp=l indicates that
the cell is fully occupied by particle). In our study, a fluid
cell size is nearly 3 times larger than the particle diameter.
Note that for liquidsolid flows, the fluid phase should be
treated as incompressible, which requires a solution
procedure that is slightly different from what is described
in Hoomans et al. (1996) and Van der Hoef et al. (2006).
In eq. (2), f represents the interphase drag coefficient.
Three correlations for the drag coefficient are considered
in the present paper. The first one is the combined Ergun
(1952) Wen & Yu (1966) correlation suggested by
Gidaspow (1994) (eq. (3)). This correlation is deduced
from static fluidized beds and has a slight discontinuity at
Ep=0.2.
Paper No
Paper No
[3 PfEfEpluvJ 265
CD dp f
15o0 +1.75 lv
rdp d,
p0 0.2
E >0.2
where p=1g and d, are solid volume fraction and particle
diameter, respectively. Eq. (3) is the most common choice
for fluidized bed simulations, with the coefficient CD
calculated from the Schiller & Naumann (1933)
expression,
24 (1+0.15Re0687)
CD = Rep
0.44,
Rep <1000
Rep >1000
Re =EfPfd
Lf
(5)
Note that in DPM, the porosity f (also called "voidage") in
a cell is calculated by using fLVa which is the fractional
volume of particle a residing in the cell under
consideration
=_V1
Vce accell
(6)
The second drag correlation used in our study originates
from Lattice Boltzmann simulations and tested for
Rep<1000 and 0=0.1 0.65, in static bed of monodisperse
particles (see Beetstra et al., 2007)
fe2 3 e (+1.5 p
PBeets =180 p+18 d2
+0.31 gp Rep f1 +3EfEp +8.4Re 0343]
efd [1+10, Re 052]
(7)
The third one is developed by Di Felice (1994) and based
on liquid fluidization experiments which makes this
correlation particularly popular in sedimentation and liquid
fluidization studies
3 PfEplUvo X
PFehce 4 CD dp 
(8)
4.8 2
CD= 0.63+ 4 '
2
(1.5log,(,Rep,)2
x =2.7 0.65 exp L 2 j
Particle equations:
The motion of each spherical particle is calculated from the
second Newton law of motion
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
dva
ma = F, +FD + F +Fp + F +Flu
(9)
All the relevant forces applied on a solid particle of mass
m, in a liquid phase are taken into account. The forces
appearing on the right handside are due respectively to
pressure, drag, gravity, particleparticle collision, added
mass and lubrication effects. The lubrication force is based
on the Stokes interaction between a pair of particles that
are close together (Nguyen and Ladd, 2002).
Particleparticle collision dynamics are described by
collision laws, which account for energy dissipation due to
the nonideal particle interaction by means of the empirical
coefficients of normal and tangential restitution, and the
coefficient of friction. The particle collisions are
represented by a spring/dashpot system following the
softsphere model (see Van der Hoef et al., 2006). Small
overlapping between two particles is allowed (less than 2%
of the average particle radius).
The Virtual Mass force is a resistance which rises due to an
"added mass" of liquid that has to be accelerated when a
particle accelerates. According to Auton (1983), the
resulting added or virtual mass force can be modeled using
eq. (10).
F. +i1 Vu
(10)
where I=CvMpfVp(vu). CVM is the virtual mass coefficient
related to the shape of the particle, and it is equal to 0.5 for
nondeformable spherical particles. The material derivative
in this equation for the virtual mass force should be the
derivative pertaining to the particle.
This force consists of the sum of two contributions. The
first part mappdv/dt is added directly to the particle
acceleration on the left hand side of the equation of motion
(eq. (9)) and contributes only to increase the particle inertia
by mapp=CvMpfVp, the apparent particle mass. The second
contribution mpp (du /dt + (u v)Vu) is added as a
supplementary contribution to the particle acceleration.
The lighter the particles are the more relevant the
contribution of the added mass to the total forces becomes.
Note that switching on the added mass force in the DPM
simulations slows down the computations by a factor 1.5,
at least for the simulations described in this paper.
Our unresolved discrete particle model does not solve for
the fluid flow perturbation induced by the particle presence
at a scale below that of the particle diameter. However, in a
dense suspension, qualitatively important physics occurs at
small particle separations, such as pairwise lubrication
forces. Hence these forces, which are repulsive (resp.
attractive) for approaching (resp. depleting) particles, have
to be included explicitly in the model, using the expression
of the lubrication force based on the Stokes flow (the
complete development of the lubrication effects is
described in Kim and Karilla, 1991). The Stokes
approximation is fair in our case since the Reynolds number
at the scale of the particles is smaller than one. Only the
normal component of the lubrication force is considered
Paper No
and it is written to the first order following Nguyen and
Ladd (2'" '2):
F 67cwt [(1 l v,)n] n h < hu,
Flub,. 6 (a+a 2 h h
0 h>h,t
(11)
where n is the normal vector connecting the particle
centers, and a, and a2 are the respective radii of nearby
particles. Eq. (11) is written as a lubrication correction of
the form of a difference between the lubrication force at gap
h and the force at some cut off distance h,,t. The lubrication
force applied on a particle close to the wall is
67f ,a2 [ 1 I [v1.n]n h
Flub,1= h11 hutJ
0 h >hut
(12)
System and Simulation Parameters
We used the discrete particle model (DPM) as described in
Hoomans et al. (1996) to simulate the fluidization of 0.165
mm glass particles in water. The simulated fluidized bed is
rectangular (the width is 10 times larger than the thickness)
similarly to the experiments of Zivkovic et al. (2" 1.r A
comparison at a onetoone scale with the experimental
case is out of reach of our computational resources due to
the high particle number; hence we used a smaller system
for numerical simulations. Different superficial velocities
Uo were used to study the system (2
actual fluidparticle mixture is characterized by low inertia
at the scale of the particles, i.e. low Reynolds Rep and
Stokes Stp dimensionless numbers defined at the scale of
the particles (0.04
Reynolds number is defined according to eq. (5). The
Stokes number characterizes the time that a particle needs
to respond to a local fluid velocity perturbation and it is
defined as
St = PP 8
p Pf 3C
(13)
where CD is the drag coefficient defined by the equation of
SchillerNeuman (1933) (eq. (4)), depending on both the
solid volume fraction and the particle Reynolds numbers.
The fluidized bed properties are shown in Table 1. This
table also contains the parameters settings for most
simulations in this paper (restitution coefficient, bed size
and particle number ...). However the influence of varying
these parameters is considered in the section "Influence of
the Simulation Parameters".
X 17 NX 30
Fluidized bed Y 1.7 mm NY 3 meshes
Z 40 NZ 80
Particles dp=165pm N,= 44300 pp=2.5g/cm3
Collisions en= kn16kt=
e=l 100 N/m
Fluid pf=lg/cm3 Uo: variable
Time steps dt,=2.5x106s dtf =2.5x105s
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Table 1
Properties and simulation parameters of the fluidized beds.
Figure 1
Front view of the 3D fluidized bed captured from a simulation
where the superficial velocity is U,,=:linii. and the Di Felice
(1994) drag correlation is used. The snapshots a), b), c) and d) are
consecutively taken at t = 0, 5.9, 12.2 and 18.7s.
The fluidized beds were typically simulated for a period of
30 seconds. Each simulation took around 30 computation
days on the Huygens cluster at SARA, Amsterdam (IBM
Power6, 4.7 GHz, with 256 GB per computation node).
The time steps used to solve the fluid and particle
equations (dty and dt, respectively) were chosen small
enough to guarantee convergence of the iterative
procedures, and also keep the particleparticle overlap
below 1% of the diameter. The fluid equations are solved
using freeslip boundary conditions on the wall cells,
Paper No
otherwise the fluidized system would be highly damped in
the depth direction (Y/X=0.1), and the particle velocity
fluctuations would be considerably reduced. Constant
velocity and atmospheric pressure are respectively
imposed at the inlet and the outlet of the channel. The
particles are initially stacked at the bottom of the bed in a
regular array. When a constant fluid velocity is imposed at
the inlet of the system, the particle bed expands until it
reaches a steady state. The numerical model gives access to
all information related to the discrete and continuous
phases. It is hence possible to calculate all the statistical
quantities from the space and / or time averages (pressure
drop, bed height, stress developed in the system ...). These
macroscopic quantities are calculated after the system
reaches its steady state.
Bed Expansion
The bed height is computed from the position of the
interface between fluidized particles and clear fluid at the
top of the bed. Because of the relatively small density
difference between the particles and the liquid (when
compared to gas fluidized beds), the transient time during
which the bed expands before reaching a steady state is
very long in liquidsolid fluidization, unlike the gassolid
systems. Figure 2 shows a typical bed expansion curve.
0 5 10 15 20 25 30 35
time
Figure 2
Evolution of the bed height in time. The drag correlation of Di
Felice (1994) is used and the superficial velocity is U,,=',n.. i
In order to check the homogeneity of the solid phase
distribution, we plotted the crossaveraged solid volume
fraction along the fluidized bed, where the solid volume
fraction e, 1 fis calculated for each cell using eq. (6). In
this, we checked that there is no significant variation of the
solid volume fraction in the cross section direction. Figure
3 shows that, for each superficial fluidization velocity Uo,
the solid volume fraction is nearly constant along the bed,
and its value increases when Uo is decreased. The most
striking observation though is the fluctuation in the solid
volume fraction, which becomes more pronounced near the
bed surface. Since the wavelength of these fluctuations
coincide with the length of a CFD cell, it is clear that this
structure is an artifact of the discretization method, and
inherent to the implementation of the fluidsolid interaction
in the DPM model. In particular, for calculating the drag
force on a particle, a Eulerian Lagrangian mapping is
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
applied, which may introduce a (slight) artificial
dependency of the drag force on the precise location of a
particle with respect to the Eulerian grid. This problem
seems to be hardly relevant for gassolid fluidized beds
(see for example Ye et al. 2005) where the particles are
much more agitated. However for these liquidfluidized
beds, where the particles are nearly static due the almost
exact balance of the drag force and the gravity force, only
a slight bias in the forces due to the underlying lattice may
produce a "tipping of the scale", and create a local
structure.
So our simulation reveal that the implementation of the
twoway coupling in the standard DPM model needs to be
further improved, in particular when applied to more
quiescent systems like liquidfluidized beds.
Liquidfluidized beds are far more critical system for the
EulerianLagrangian mapping,
Nevertheless, we checked that global properties (so
averaged over the bed domain) like the mean bed voidage
or the mean granular temperature, are hardly influenced by
these structures, so that it is still meaningful to perform a
statistical study of these properties.
0.1 0.15 0.2 0.25
Solid Volume Fraction
0.35
Figure 3
Variation of the crossaveraged solid volume fraction with height.
The curves correspond from right to left to Uo=4, 5.37, 6, 7 and 8
mm/s respectively. The corresponding simulations are run using
the drag correlation ofDi Felice (1994).
The mean bed voidage , i.e. the average fluid volume
fraction in the fluidized part, can be determined either from
the bed height H following eq. (14) where A is the
crosssection area and V, is the particle volume,
N, .V
=1 Np.Vp
HA
(14)
or from the void fraction (eq. (6)) averaged over all the
cells containing particles ( 1= where =). The
brackets <.> denote the time average over the ensemble of
particles after the system has reached its steady state. Both
methods lead to equal results for the bed height.
Series of simulations using different drag correlations were
run to find the relationship between fluidization superficial
velocity Uo, and overall bed voidage . Figure 4 shows
the numerical results obtained using different drag
correlations (eqs. (3), (7) and (8)) compared to the
experimental results of Zivkovic et al. ( 'i ,).
Paper No
10
8 V A
E 6A
4 3
2 A
8.4 0.5 0.6 0.7 0.8 0.9
Average Fluid Volume Fraction
Figure 4
Variation of averaged bed voidage with superficial fluid
velocity Uo. The results of the numerical simulations are obtained
using the drag correlations of Di Felice (1994) (0), Beetstra et al.
(2007) (A) and Ergun(1952) Wen & Yu(1966) (0). The
simulation results are compared to the experimental results (*) of
Zivkovic et al. (2008) fitted with a solid line.
The solid line in Figure 4 represents the correlation of
Richardson and Zaki (1954)
U =(
Ut
(15)
that best fits the experimental results with an exponent
n=4.040.03 and a particle terminal settling velocity of
Ut=16.10.2mm/s. These values are in agreement with the
prediction of Khan & Richardson (1990). Figure 4 shows
that the expansion curves obtained by the DPM
simulations follow the right trend. For a given superficial
velocity Uo, the solid volume fraction is in general slightly
underpredicted by the simulations. Nevertheless, the
expansion curves depend on the drag correlation
implemented for the coupling between the fluid and the
particle equations. The values of the exponent n and the
settling velocity Ut of each case are given in Table 2. The
expansion coefficient n is overpredicted by the
simulations. As for the settling velocity, it is only
underestimated by the drag correlation of Beetstra et al.
(21" 1) and overestimated by the others simulations.
Drag correlation n Ut (mm/s) R2
Di Felice (1994) 5.47 18.07 0.9993
ErgunWen&Yu (1966) 5.1 19.47 0.9974
Beetstra et al. (2007) 4.81 13.14 0.9993
Exp. Zivkovic et al. (2008) 4.04 16.1 0.9987
Table 2
Exponent n and particle terminal settling velocity U, obtained by
fitting eq. (15) with the numerical simulations for different drag
correlations. R2 is the regression coefficient.
Granular Temperature
The granular temperature or the agitation of the solid phase,
results from the different interactions in the system
consisting of the particles, the fluid and the wall. Similarly
to the Ogawa granular temperature used in Zivkovic et al.
('I I ), we define the granular temperature from the
particle velocity fluctuations v'a about the average of the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
particle velocities
V'a = Va < V >
T ((V' VJ K') Ka))
(16)
where a {x, y, z} which represent respectively the
directions of the width, depth and height of the bed. In the
simulations we have checked that the particle averaged
velocity is negligible and that the products
are at least 104 times less than the crossproduct .
Hence we can simply write the granular temperature as
T=II^ 'v1'
3 ,
Before exploring the granular temperature, we checked the
state of the particle velocity fluctuations in the fluidized
bed. The normalized velocity fluctuation distribution
functions calculated for the particular case of Uo=6mm/s
are shown in Figure 5 (the fluctuations in the width v and
flow v' directions are noted u and w in Figure 5). These
velocity fluctuations follow nearly a Gaussian distribution.
However the velocity fluctuations close to zero are slightly
underpredicted by the Gaussian function. The peak at
small velocity fluctuations suggests the existence of some
regions in the fluidized bed where the particle agitation is
very small. The velocity fluctuations v' are not explored
here since they are very small compared to the other
components due to geometry limitations in the depth
direction (Y/dp 10 whereas X/dp,100).
10
10
10 221 0 1 2 3
3 2 1 0 1 2 3
u/GJ,
2 1 0 1 2
wlc
Figure 5
Probability density function (PDF) of the velocity fluctuations in
the width direction (left) and in the flow direction (right). The
velocity fluctuations are normalized by their respective standard
deviation. The plus signs are the numerical simulations (obtained
with Uo=6mm/s and using the drag correlation of Di Felice
(1994)) and the dashed curves are the normalized Gaussian
functions.
The evolution of the squareroot of the granular temperature
in the width direction is shown in Figure 6 for different
heights. It can be seen that close to the bottom of the bed,
the particles have low agitation level with a minimum at the
center of the fluidized bed. However at higher levels, all the
granular temperature curves are similar. They exhibit local
maximum at the center of the fluidized bed and local
minima at both sides of the bed. Close to the wall, the
particle agitation is the highest, suggesting that a flow
bypass is taking place there, due to the fluid freeslip
boundary conditions at the walls. The clear apparition of
local peaks in the curve of the granular temperature is an
evidence of a large scale motion. Indeed, the movies of
these simulations show a clear upward motion of the
particles in the center of the bed, a downward motion on
both sides and an upward bypass close to the walls. Such a
Paper No
large scale motion is a phenomena commonly observed in
liquidsolid fluidization (Aguilar(_' i1 ,,), Latif & Richardson
(1972), Wilus (1970), Handley et al. (1966)).
n
0.005 0.01
Width(m)
0.015 0.02
Figure 6
Variation of the squareroot of the granular temperature in the bed
width direction, for superficial velocity U,,=,iiii The dotted,
solid, dashed and dasheddotted lines are computed respectively at
9, 27, 45 and 64% of the height of the liquidsolid mixture.
In Figure 7 we show results for the average granular
temperature as a function of the superficial velocity Uo.
With "average" we mean that it is averaged over all
particles in the domain, and averaged over time for the
period that the bed is in steady state (roughly from 15 to 25
seconds). The experimental results of Zivkovic et al.
(2" "I ,) are plotted on the same figure. Although they do not
match exactly the experimental values, the numerical
results follow the right trend. The values of the granular
temperature depend strongly on the choice of the drag
correlation, especially at high superficial velocities. Most
of the simulations were run using the correlation of Di
Felice (1994) which has been constructed for dynamic
fluidized beds. Similarly to the experimental results, the
rate of increase of the granular temperature is much higher
at low superficial velocities. A kind of saturation in particle
agitation is reached at high superficial velocities.
10
0
8
*""
6
E A
o
Uo (mm/s)
superficial velocity Uo. The + signs linked with a dotted line are
points are from the numerical simulations results using the drag
A p
2
6'
2 4 6 8 10
U0 (mm/s)
Figure 7
Evolution of the squareroot of granular temperature T1/2 with
superficial velocity Uo. The + signs linked with a dotted line are
the experimental results of Zivkovic et al. (2008). The remaining
points are from the numerical simulations results using the drag
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
correlations of Di Felice (1994) (0), Beetstra et al. (2007) (A) and
Ergun (1952) Wen & Yu (1966) (0).
Gevrin et al. (2" I) tried to model the granular
temperature as a function of the solid volume fraction
based on a global energy balance in the liquidsolid
fluidized bed. The authors separated between large and
small scale velocity fluctuations. They found that the
particle velocity fluctuations at the small scale 1/2 are
proportional to Ut, (CD))13 and (1_/m)n"m, where n and
U, are the same exponent and settling particle velocity
appearing in the correlation of Richardson and Zaki (1954),
Om is the maximum packing fraction (set to 0.64) and CD is
the drag coefficient given by Schiller & Neumann equation
(eq. (4)). The coefficients of proportionality written in
Gevrin et al. (2'" I,) depend on the nature of the fluidized
particles, and were found by fitting the expression of the
granular pressure (depending on /2) with the
experimental results. Later on, the authors have revisited
their correlation, and they fitted the expression of l/2
with the experimental results carried out for different
particle types and sizes. The expression of the small scale
velocity fluctuations for low particle inertia is given by
Aguilar (2, 1 ,,):
Figure 8 shows the numerical results of T1/2 compared to
the theoretical prediction of the velocity fluctuations which
presents a maximum value for 0 close to 0.1. The existence
of the peak was also observed in the work of Zivkovic et al.
(2" ,). Gevrin et al. (i ,) explained that this maximum
is related to the production of the particle fluctuating
motion, in the absence of turbulence, from the local
instantaneous variations of the solid volume fraction which
induce local velocity gradients (to satisfy the continuity
equation) and hence kinetic energy fluctuations of the solid
phase. The increase of the solid volume fraction in "dilute
regime" leads to an increase of the granular temperature.
At high concentration, the fluctuations of the solid fraction
are limited by the increase of the collision frequency,
reducing the mean free path of the particles.
10
8
6 n@
n E
E3D
n s
0.1 0.2 0.3
Solid volume fraction
0.4 0.5
Figure 8
Evolution of the square root of the granular temperature with the
average solid volume fraction. The numerical simulations using the
drag correlations ofDi Felice (1994) (0) are compared to eq. (17).
q< p>= 0.45U, (CD)" (1 / m)
Paper No
Figure 8 shows an acceptable agreement between the
simulations and the prediction of eq. (17). Nevertheless the
following differences should be kept in mind. On the one
hand, the solid line is the prediction of 1/2 which
represents the granular temperature based only on the
small scale fluctuations. However the velocity fluctuations
T1/2 deduced from our simulations are based on the global
velocity fluctuations (including all scales). On the other
hand, Aguilar (2'" i ") has found equal velocity fluctuations
in the crosssection planes whereas the velocity
fluctuations in the depth direction practically vanish in our
simulations. It seems from Figure 8 that these differences
cancel each other, in such a way that the simulation results
are not far from the theory prediction.
The advantage of the numerical simulation over the
experiments such as those of Zivkovic et al. .2 I ), is that
in simulation one can also measure the average
components of the velocity fluctuations in different
directions, i.e. Txx=, Tyy= and Tzz=,
where we remind that x, y and z are the width, depth and
height directions, respectively. In a homogeneous fluidized
bed, the transverse velocity fluctuations (Txx)1/2 and (Tyy)1/2
are generally found nearly equal. However in our
simulations the velocity fluctuations in the depth direction
(Tyy)1/2 are negligible compared to other components
((T/Tyy,)/2and (Tzz/Tyy)1/2being larger than 10). Hence in
the actual system most of the energy transmitted to the
particles is redistributed as a fluctuating motion in the
plane formed by the flow and the bed width directions.
Figure 9 shows the anisotropy coefficient, i.e. (Tz/Txx)1/2
obtained from the DPM simulations. The value of this
coefficient fluctuates around the value of 1.89. This
constant is of the same order of magnitude of the
anisotropy coefficient found experimentally by Aguilar
2' 1').
2.5
2
x
N
1.5
1
n
.0.1
0.15 0.2 0.25 0.3
Solid volume fraction
0.35 0.4
Figure 9
Variation of the anisotropy coefficient (Tzz/Txx)12 with the solid
volume fraction, where z and x are the height and width
directions repsectively. The dashed line is the average of the
anisotropy coefficient obtained from the DPM simulations using
the drag correlation ofDi Felice (1994) (0).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Influence of Simulation Parameters
Added mass and lubrication forces
Although it is common practice to neglect the added mass
and lubrication forces in liquidsolid mixtures, the validity
of this assumption was never tested for the liquidsolid
fluidized beds particularly. In order to do so, we added
both forces to the numerical code, and looked at the
differences when they are switched on or off. We found
that the influence of adding the virtual mass and
lubrication forces on the fluidized bed behavior is not
significant: for instance the variation of the bed height is
less than 8% and the granular temperature less than 5% for
Uo=5.37mm/s. Note that the values of the relative error are
only indicative due to the uncertainty in the numerical
measurements mainly related to the DPM formulation.
Restitution c ,iffic ient
In the gassolid systems, the energy provided to the
fluidized bed via the superficial velocity is dissipated
mainly via the particle collisions. The energy dissipation is
usually expressed in terms of a restitution coefficient
which is the ratio of the velocity after and before particle
surface contact. However in liquidsolid fluidization, the
contribution of particleparticle interaction to dissipation is
far less relevant than the contribution of the hydrodynamic
interactions (fluidparticle interaction, i.e. the drag
dissipation). We considered the effect of the restitution
coefficient on the macroscopic behavior of the fluidized
bed. Note that the restitution coefficient in our study does
not take into consideration the lubrication effects (as in
Legendre et al. (2006), Abbas et al. (2010)). These effects
are taken into account in the lubrication forces. The
numerical simulations show that by changing the
restitution coefficient from 0 to 0.97, the granular
temperature and the bed height are influenced by less than
1% which is in agreement with our a priori expectation.
Size of the bed
As already mentioned, the particle velocity fluctuations are
only partly induced by the particleparticle interactions
(small scale fluctuations). The remaining part is due to the
bed confinement effect which induces large scale
fluctuations inside the fluidized bed. The movies resulting
from DPM simulations show clearly the presence of large
scale motion. The quasi2D shape of the fluidized bed
enhances the existence of recirculation zones. The scale
separation of velocity fluctuations are not studied in detail
in this work. Nevertheless we checked if the wall has an
effect on the macroscopic bed behavior.
In order to get an approximate idea about the influence of
the bed size on the macroscopic behavior of the fluidized
bed, a larger system is considered ((34x3.4x70) mm3 with
354000 particles), while the properties of the system are
kept the same. The bed expansion observed is comparable
to the results with smaller bed size, with an exponent
n=6.06 and a terminal settling velocity Ut=29.54mm/s. For
the granular temperature, Figure 10 shows a better
agreement with the experimental results except for the case
with Uo=5.37mm/s. We believe that this difference is due
a a
DD nl
n] [n n
FL
'4
Paper No
to the large celltoparticle size ratio used for this case:
AX/dp=5 for Uo=5.37mm/s whereas AX/dp=3.5 for other
cases (Uo=2.06mm/s and Uo=7.87mm/s).
10
0
8o
6
E
+'* V
I
2 ,
'4.
0 mrm
0 2 4 6
Uo (mm/s)
8 10
Figure 10
Influence of the bed size on the computation of the granular
temperature using the drag correlation of Ergun (1952)Wen & Yu
(1966). The + symbols related by dotted line are as in Figure 7.
The other symbols are obtained with different
simulationtoexperiment bed size ratio: 1/6 (V) and 1/12 (0)
(same as in Figure 7).
Conclusions
This paper shows the validity of using a DPM approach
(discrete particle model combined with computational fluid
dynamics for continuous phase) to simulate liquidsolid
fluidization. The bed expansion was reproduced by the
simulations and compared to the experiments of Zivkovic
et al. (2'" 1) for different drag correlations.
The particle velocity fluctuations were calculated. The
velocity fluctuations are anisotropic, the strongest
fluctuation being in the flow direction. Similarly to
experimental observations, both small and large scale
fluctuations were found in the simulations. The granular
temperature was obtained by averaging the variance of the
velocity fluctuations over all the particles. The variation of
the granular with the superficial velocity shows a similar
trend as the experiments despite the different system sizes:
the numerical simulated bed is 12 times smaller than the
experimental bed. The numerical results of granular
temperature show good agreement with the theoretical
prediction of Aguilar (2i' 11,) and Gevrin et al. (2'" ".
The added mass and lubrication forces were included in the
model and no significant impact on the macroscopic
behavior of the fluidized bed was observed. Also the
particle restitution coefficient was not found to have a
significant impact on the results, in contrast to
gasfluidized beds.
The DPM has shown some weakness like the oscillations
observed in the solid volume fractions curves (Figure 3).
Such symptoms can be avoided by improving the twoway
coupling between the fluid and particle equations. It would
then be possible to measure accurately the local solid
volume fraction and pressure fluctuations, which are
fundamental to understand some phenomena that are
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
behind the local agitation and instability inside the
fluidized bed systems.
Acknowledgements
M. A. would like to thank O. Masbernat for fruitful
discussions on the theoretical background of agitation in
fluidized beds.
References
Abbas, M., Climent, E., Parmentier, J.F. & Simonin, O.
(2010). Flow of particles suspended in a sheared viscous
fluid : effects of finite inertia and inelastic collisions, AIChE,
under press.
Aguilar, A. (2'" 1"). Agitation des particles dans un lit
fluidise liquid: etude experimental, Ph.D. thesis of
University of Toulouse, Toulouse, France.
Auton, T. R., (1983). The dynamics of bubbles, drops and
particles in motion in liquids, Ph.D. thesis, University of
Cambridge, Cambridge, United Kingdom.
Beetstra R., van der Hoef M.A. & Kuipers J.A.M. i ').
Drag force of intermediate Reynolds number flow past
mono and bidisperse arrays of spheres, A.I.Ch.E Journal 53,
489501.
Carlos, C.R. & Richardson J.F. (1968). Solids movements in
liquid fluidized beds I Particle velocity distribution, Chem.
Eng. Sci., 23, 813824.
Di Felice R. (1994). The voidage function for fluidparticle
interaction systems, International Journal of Multiphase
Flow, 20, 153159.
Di Renzo A. & Di Maio, F.P. (2. i'). Homogeneous and
bubbling fluidization regimes in DEMCFD simulations:
hydrodynamic stability of gas and liquid fluidized beds.
Chem. Eng. Sci., 62, 116130.
Ergun, S. (1952). Flow through Packed Columns, Chem.
Eng. Prog., 48, 89.
Gevrin, F., Masbernat, O., & Simonin, O. (2i'"). Granular
pressure and particle velocity fluctuations prediction in
liquid fluidized beds, Chem. Eng. Sci., 63, 24502464.
Handley D., Doraisamy, A., Butcher, K.L. & Franklin, N.L.
(1966). A study f the fluid and particle mechanics in
liquidfluidized beds, Trans. Inst. Chem. Eng., 44,
T260T273.
Hoomans B.P.B., Kuipers J.A.M., Briels W.J. & Van Swaaij
W.P.M. (1996). Simulation of bubble and slug formation in
a twodimensional gasfluidized bed: a hardsphere
approach, Chem. Eng. Sci., 51:99118.
Kang, Y, Nah, J.B., Min, B.T. & Kim, S.D. (1990).
Dispersion and fluctuation of fluidized particles in a
liquidsolid fluidized bed, Chem. Eng. Commun., 97,
Paper No
197208.
Khan, A.R. & Richardson, J.F. (1990). Pressure gradient and
friction factor for sedimentation and fluidisation of uniform
spheres in liquids. Chem. Eng. Sci. 45 (1), 255265.
Kuipers, J.A.M., Van Duin, K.J., Van Beckum, P.H., Van
Swaaij, W.P.M. (1993). Computer simulation of the
hydrodynamics of a twodimensional gasfluidized bed,
Computers Chem. Eng., 17, 839258.
Latif, B.A.J., & Richardson, J.F. (1972). Circulation patterns
and velocity distribution for particles in a liquid fluidized
bed, Chem. Eng. Sci., 27, 19331949.
Lettieri, P., Di Felice, R., Pacciani, R. & Owoyemi O.
(2006). CFD modeling of liquid fluidized beds in slugging
mode, Pow. Tech., 167, 94103.
Legendre, D., Zenit, R., Daniel, C. & Guiraud (2006), A
note on the modelling of the bouncing of spherical drops or
solid spheres on a wall in viscous fluid, Chem. Eng. Sci., 61,
35433549.
Mukherjee, A.K., Mishra, B.K. & Ran Vijay Kumar (2i 1 r'.
Application of liquid/solid fluidization technique in
beneficiation of fines, Int. J. Mineral Processing, 92, 6773.
Nguyen N.Q. & Ladd A. J. C. (2'"'2). Lubrication
corrections for latticeBoltzmann simulations of particle
suspensions, Physical Review E, 66, 0467081:12.
Pan, T.W., Joseph, D.D., Bai, R., Glowinski, R., Sarin, V
('"'2). Fluidization of 1204 spheres: simulation and
experiment, J. Fluid Mech., 451, 169191.
Pronk, P, Infante Ferreira, C.A., Witkamp, G.J.
J. Heat and Mass Transfer 52, 38573868.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Zenit, R., Hunt, M.L., Brennen, C.E. (1997). Collisional
particle pressure measurements in solidliquid flows, J.
Fluid Mech. 353, 261283.
Zivkovic, V, Biggs, M.J., Glass, D., Pagliai, P. & Buts, A.
(2'" i'). Granular temperature in a liquid fluidized bed as
revealed by diffusing wave spectroscopy, Chem. Eng. Sci.,
64, 1102 1110.
(2009). Int.
Richardson, J.F. & Zaki, W.N. (1954). Sedimentation and
fluidization: part I. Transactions of the Institute of Chemical
Engineers, 32, 3553.
Schiller, L., & Neumann, A. (1933).Uber die grundlegenden
Berechungen beider Schwer kraftaufbereitung," Verein
Deutscher Ingenieure 77, 318.
Van der Hoef, M.A., Ye, M., Van Sint Annaland, M.,
Andrews, A.T., Sundaresan S. & Kuipers J.A.M. (2006).
Multiscale modelling of gasfluidized beds, Adv. Chem.
Eng., 31, 65149.
Wen, C. Y. & Yu, Y H. (1966). Mechanics of Fluidization,
Chem. Eng. Prog. Symp. Ser., 62, 100.
Wilus, C.A. (1970). An experimental investigation of
particle motion in liquid fluidized bed. PhD thesis of
California Insitute of Technology.
Ye, M., Van der Hoef, M.A., Kuipers, J.A.M. (2'" ). The
effects of particle and gas properties on the fluidization of
Geldart A particles, Chem. Eng. Sci., 60, 45674580.
