7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A Priori Test of Effective Drag Mlodeling for Filtered TwoFluid Mlodel Simulation of
Circulating and Dense GasSolid Fluidized Beds
A. iizel", J.F. Parmentier ', O. Simonind and P. Fededi
University de Toulouse: INPT, UPS: IMFT: 31400 Toulouse, France
SCNRS: Institut de M~canique des Fluides de Toulouse: 31400 Toulouse, France
ozeld imft.fr, parmend imft.fr, simonin a uft.fr and feded imft.fr
Keywords: Twofluid model, LES, fluidized bed, drag force
Abstract
Eulerian twofluid approach is generally used to simulate gassolid flows in industrial fluidized beds. Because of
limitation of computational resources, simulations of large vessels are usually performed using too coarse mesh
to capture the influence of the fine flow scales which can play an important role in the dynamic behaviour of the
beds. In particular, neglecting the particle segregation effect at small scale leads to an inadequate modelling of the
mean interfacial momentum transfer between phases. Then, an appropriate modelling approach which accounts for
influences of unresolved structures has to be proposed for "coarse simulations". For this purpose, computational grids
are refined to get meshindependent results for a dense and a periodic circulating fluidized beds in which statistical
quantities do not change with further mesh refinement. These meshindependent results are filtered by volume
averaging and then used to perform a priori analysis on the filtered drag term. Results show that filtered momentum
equation can be computed on "coarse simulations" but must take into account the particle to fluid drift velocity due
to the sub grid correlation between the local fluid velocity and the local particle volume fraction. In the present paper
we propose a model, for subgrid the drift velocity, written in terms of the difference between the averaged of gas
velocity weighted by solid volume fraction and the averaged of gas velocity weighted by gas volume fraction. We use
a systematic procedure to provide constitutive closures for the subgrid drift velocity and the closure depends of both
the filtered solid volume fraction and a characteristic filter size.
Nomenclature
117,; sub grid drift velocity
Subscripty
<; Gas
p Particle
Symbols
Si
6;y
As
gravity
mean volume fraction of phase k
Kronecker delta
bulk viscosity of phase k
viscosity of phase k
Introduction
Numerical analysis of fluidized beds by Eulerian two
fluid approach appear to be a powerful tool to improve
design and performance of industrial facilities. How
ever, in most industrial applications involving large de
vices, because of limited computational resources, two
fluid model equations for unsteady gasparticle flows
in dense and circulating fluidized bed risers are rou
tinely simulated over too coarse spatial grids which can
not resolve all the finescale structures. Comparisons
with finegrid simulation results showed that unresolved
structures of gasparticle fluidized bed can have a dras
<1, particle diameter
e, elasticity coefficient during interparticle collision
Isinterphase momentum exchange
P, gas pressure
q~i disperse phase fluctuations
Re particle Reynolds number
7,, characteristic time scale of
gasparticle momentum transfer
us,; velocity ofphase k
Ux, = ux,) Amean velocity of phase k
11; mean relative velocity
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
tic influence on the flow dynamic due to the inadequate
modeling of the resolved drag term coupling the gas and
particle resolved momentum equations (Andrews et al.
2005). To account for the effect of unresolved struc
tures on the macroscopic behavior in coarse grid simula
tions, suitable sub grid model can be proposed to model
accurately the effective drag term. In this study, a pri
ori analysis of the effective drag modeling is performed
by spatial filtering of highly resolved 2D and 3D nu
merical simulations of both dense and periodic circu
lating fluidized beds. The numerical simulations have
been carried out using an unsteady Eulerian multifluid
approach implemented in the unstructured parallelized
code NEPTUNE CFD V1.ll' &Thei NEPTUNE CFD
is a multiphase flow software developed in the frame
work of the NEPTUNE project, financially supported by
CEA (Commissariat g l'Energie Atomique), EDF (Elec
tricit6 De France), IRSN (Institut de Radioprotection et
de Stiret6 Nucl~aire) and AREVANP.
as follows:
i3 i3
( i)P ?)~i
at  p ;+ A;+ R
with P, the mean pressure, 9; acceleration due to
gravity, Ex,;y is the effective stress tensor, IA~i is the
mean transfer momentum transfer rate after substraction
of the mean gas pressure effect.
Interfacial Momentum Transfer
The term It s accounts for momentum transfer rate
between the gas phase (carrier or continuous phase) and
particle (disperse) phase. This term can be modeled by
taking into account only the drag force between phases:
Qppp
7
Qp
Mathematical Modeling Approach
The modelling approach is based on the twofluid
model formalism that involves mean separate transport
equations of mass, momentum and energy for each
phases. Interactions between phases are coupled
through interphase transfers. The transport equation for
disperse phase fluctuations, q~ developed in the frame
of kinetic theory of granular media supplemented by
the interstitial fluid effect and the interaction with the
turbulence (Balzer et al. (1995), Gobin et al. (2003)),
are resolved with taking into account interparticle
collisions on the dispersed phase hydrodynamic. The
effect of the fluctuations of the gas velocity at small
scales is neglected. Concerning the transfers between
the phases with nonreactive isothermal flow, drag
force was only taken into account for the transfer of
momentum. Neither terms of added mass nor lift force
were considered.
Transport Equation
The twofluid mass balance equation for the phase k is
written:
(asPk)+ (asPaLT<,;)= 0 (1)
at dir
with asi, p e, ET the volume fraction, the density and the
mean velocity of phase k (when subscript k = g, we
refer to the gas and k = p to the particle phase). The
momentum balance equation for the phase k is defined
with the particle relaxation time scale, 7,'
aMI, di V  and (.) the ensemnble aver
ag~e orator over the particle (Simonin 1996).
The mean drag coefficient of a single particle,
(C D) can be written as function of a particle
Reynolds number and defined by Wen & Yu (1966),
Dp)i 24/(Rr p) 1+ 0.15(Reg~6i Qr) a with
the definition of the mean particle Reynolds number,
(Re p) = ad,(. \p/v where d, is the particle
diameter and v, is the molecular viscosity of gas,
respectively. In this study, only the Wen&Yu correlation
is used while the combination of Ergun's and Wen&Yu
correlations is generally employed in fluidized beds
(Gobin et al. 2003). The term, (vr )p, in (3) represents
the local instantaneous relative velocity, v,,;, and is
equal to difference between the local particle velocity,
undisturbed by presence of the particle at the particle
position. V,,; is the averaged of the local relative
velocity and equal to (il,,; up,;) p. It can be expressed
in terms of the averaged velocity between phase,
Effective Solid Stress Tensor
The effective stress tensor for particle phase has two
contributions. The first contribution is the kinetic stress
tensor which represents the transport of the momentum
by the particle velocity fluctuations. The second one
is the collisional stress tensor which accounts for
destruction and exchange of the momentum during
interparticle collisions. The constitutive equation of the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Free outlet
effective stress tensor is:
1F: 17, iif 6
i~"' i
OWall
g Periodicity
H=0 22 m
1
~sij~] (4)
Walls
.oslmii'~
The constitutive relations for viscosity and diffusivity
are derived in frame of the extension of kinetic theory
of dry granular flow and given by Balzer et al. (1995).
The collisional pressure and the bulk viscosity are writ
ten according to Balzer et al. (1995) as follow:
PP = 3apPPq ~[1+299,,(1+ec)]
Xp n pylppun(1 + ec) 2q (5)
3 V 3 x
with the random kinetic energy of particles, q, the resti
tution coefficient, e 0, that determines energy loss dur
ing interparticle collisions and go is the pair correla
tion function. The shear viscosity is the linear combi
nation of the collisional and the kinetic stress: ps =
a~ppp [va~n, + y<>i] where
21~p
c 2
g it Aag (1+ec v'" + cld q (6)
5 ; 3 x
with ac = (1+ e ) (3 e ) and Ae =
(1 c)(30c 1).Such closure laws were also
used by Agrawal et al. (2001) in the case of 2Dperiodic
flows.
Flow Configurations
Gasparticle flows were simulated for two configura
tions: a twodimensional dense and a threedimensional
periodic circulating fluidized bed. Typical FCC parti
cles, <1, 75 pm, pp 1500 ky/m3, are interact
ing with the ambient gas (p, 1.186 ky/m3, p,
1.8 x 10") for both cases. The computational domains
of cases are shown in Figure 1. The dense fluidized bed
is initialized by the homogenous distribution with the
particle volume fraction equal to 0.55 and the fluidiza
tion velocity, Uf, is set to 0.2 at/s. The periodic cir
culating fluidized bed is initialized by the homogenous
distribution of the particle volume fraction equal to 0.05
and very small fluctuations applied on the particle vol
ume fraction to destabilize the bed. The flow in the pe
riodic fluidized bed is driven by the pressure gradient
Uniform gaz velocity
Figure 1: The computational domain of two
dimensional dense and threedimensional periodic cir
culating fluidized bed.
due to the total mass which is defined as the opposite di
rection to the gravity. For both cases, noslip condition
were imposed at the walls.
Mesh Independent Results
Agrawal et al. (2001) stated that statistics quantities over
the whole domain is strongly dependent on the mesh size
and they became meshindependent when mesh sizes
are the order of few particle diameters. In this study,
the mesh refinement studies were carried out to be com
pletely sure that the mesh resolution is sufficient and all
spatial and temporal scales of particle and gas phases
are captured. Figure 2 shows different instantaneous
particle volume fraction fields in the dense and periodic
circulating fluidized bed obtained by numerical simula
tions using different mesh sizes. As the resolution of
mesh increases, inhomogenous structures are better re
solved. These structures have a drastic influence on the
height of dense bed and the solid holdup of periodic bed.
The influence of mesh size on the time averaged bed
height of dense fluidized bed and the vertical solid flux,
# ff a, Up clS where the integral is performed on a
horizontal surface, in the periodic circulating fluidized
bed is shown in Fig. 3. Dramatic changes of macro
scopic behaviors of beds can be attributed to poor pre
diction of drag force. These converged results can be
called "DNS" results by analogy with single phase flows.
They are statistically sufficient and not needed to have
further mesh refinement, then used to perform a priori
analysis on the drag term filtered by volume averaging
to propose subgrid model for the modeling of the effec
tive drag term.
'7=0 0275 m
L=0 0275 m
.
0 0. .
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 2: Instantaneous particle volume fraction field in the dense (le f t) and periodic circulating fluidized bed (rig ht)
for different grid mesh sizes. >From left to right, the mesh resolution increases for both cases. White color corresponds
to ap = 0. Black color corresponds to ap = am = 0.64.
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Figure 3: The influence of grid size on the height of dense fluidized bed (le ft) and the vertical solid flux in the
periodic circulating fluidized bed brightt), 0: h~omogenlous case, conv~:: converged case, A/(ist)ag: nondimensional
mesh size with characteristic mesh size, (AxAyAz)l/ and Stokes' relaxation time, 7,st
Filtered TwoFluid Model Equations and
The twofluid model equations are spatially averaged
over some chosen filter length scales. Let a,(x, t) de
notes the particle volume fraction at location, x, and
time, t, obtained by solving the twofluid equations. We
can define the filtered particle volume fraction 0, (x, t)
ap,(x, t) Glx Ciu) a(ut) du (7)
where G is a weight function which satisfies
fff G(u) du 1 Filtered phase velocities
are defined according to
U,(xi) =X G u) a,(u, t) U (u, t) du
as . .(9)
Applying such a filter to the continuity equations of
phases, one can obtain
8agpa pkcasUlc4
0 (10)
Upx t) = GJ )a(,t pu )d Repeating this filtering process to momentum balances,
(8) the filtered momentum balance for particle and gas
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
as filtered particle volume fraction, gas and particle ve
locities at several time instants. We carry out correlation
analysis between the filtered drag term and the filtered
variables. Finally, we perform conditional average of the
filtered drag term for given values of the filtered particle
volume fraction. We propose to decompose the filtered
in two contributions: the difference between filtered gas
and particle velocity Up,4 Up,s and a sub grid drift ve
locity Vd,i defined by
phase are
i3i3
dP,
agc (P + arPpagg + Irc g +
8xy 8xyprcc~rar~ij (11)
Additional terms arise in (11) due to the filtering pro
cess and require closure models. The term, cpi, repre
sents the contribution between fluctuations of the vol
ume fractions of phases and gas pressure and defined as:
dP, dP,
(Pi =r an 8 ag 8xy (12)
The modeling of this term was proposed by De Wilde
(2005) with the introduction of the global added mass
coefficient. A Reynolds stresslike contribution com
ing from the gas or particle phase velocity fluctuations,
arc~ij, is defined by the following equation:
Ercarc,ij = ag Uc,iUlc,j ag Elc,i1p,j (13)
This term can be modeled by using Boussinesq approx
imation with introducing the subgrid viscosity, psy,,~,
for the gas phase:
appp pppy ~ p.
~ ~
U + Va
gi 
This decomposition allows us to identify physically the
meaning of sub grid drift velocity, and then to understand
the origin of the difference between the filtered and re
solved drag force. A correlative analysis has shown that
the filtered drag force can be approximated by the fol
lowing epxression.
Qpp' V apV17
7p ,i pr,i (7
Both sides of (17) are correlated with more than 99 %,
even for large filter size. It shows that the subgrid mo
mentum transfer is occurred by the filtered relative ve
locity averaged by particle volume fraction, a,V,,s. Us
ing (16) and (17), the subgrid drift velocity is approxi
mated by:
3 8xm
V4,4 g~p~i g(i8)
where Ugep~i =pUg,~i / op iS the filtered gas velocity
seen by the particle phase.
where 6ij is the Kronecker delta. The subgrid viscosity
can be modelled by subgrid scale model proposed by
Deardoff (1971):
Psas,9 = a~~Pp(ctA)ay22/SjSi (15)
where a = (AxAyAz)l/ and Sy,j,
4 [+i~dzj1'. This term for disperse phase
was investigated by Moreau (2005) with mesoscopic
Eulerian approach. The term, I,, is the filtered drag
term and we focus on this term in this study. To define
variables contributed to the subgrid drag term, correla
tion analysis can be performed between the filtered drag
term and computed variables. To fill up the main goal,
suitable model accounted for main contributions of the
subgrid drag term will be proposed by the definition
of the effective drag term. We present a priori test on
the filtered drag term and briefly show primary results
obtained by the underdeveloping model.
A Priori Analysis on the Filtered Drag term
We apply top hat volume filters with various sizes on lo
cal instantaneous drag term and computed variables such
Modeling of the subgrid drift velocity
Heynderick (2l1***4, Andrews et al. (2005) and Igci
(2007) introduce an effective drag coefficient Pe to
express the filtered drag force term as
appp ri=# 7~ ~
Heynderick (21 *4) and Andrews et al. (2005) write the
effective drag coefficient as a function of the filtered par
ticle volume fraction while Igci (2007) suggest that this
coefficient is a function of the filter size. We propose
to write the filtered drag force by modeling the subgrid
drift velocity as:
Vas = g(A, a,) ET,,s
rET,,4)
with a function g that will be determined by the condi
tional averaging procedure applied on the high resolu
tion simulations. Then, the effective drag term can be
written as follows:
appy pp
ET,, (21)
agge,ij = psys,,g 8 Us,jdz 8Ug,idz
0
~0.25
0 0.2 0.
Fitrdsld ouofato
Figure 4: The function "y" for dense fluidized bed
casefor different filter sizes: 0 : 11A/ADNS, O :
15A/ADNS, o2 2 DDN.
The function "9"
The function "g" can be calculated by the ratio between
the sub grid drift velocity and the difference between the
filtered gas and particle velocities which are condition
ally averaged by the filtered volume fraction of the dis
persed phase. This allows us to derive the following re
lation:
(Mid.
9 a, Gr)
The computed values of the function g are shown for
dense and periodic circulating fluidized bed in Figure
4 and Figure 5. The function "g" can be written as a
product of two functions: h (o,) (for the filtered vol
ume fraction dependance) and f (A) (for the filter size
effect).
U(a, Oap) = f (a) h(op) (23)
>From our database, we propose the following form of
the function h:
h(K4) = [iK," '(ap,,, a (24)
with contents m and I which are equal to 1.2 and 1.6 ,
respectively. The function h is shown for different filter
sizes in Figure 6. The function f is modeled as the
following equation:
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.1
0.2
0.3
0.4
0.5
S0.6. x
0.7 x
0.8 x
0 0.05 0.1 0.15 0.2 0.25
Figure 5: The function "g" for periodic circulating flu
idized bed case for different filter sizes; : 3A/ADNS,
x :5A/ADNS,* /DNS 09/DNS.
1.8
1.6 *
1.4 
1.2
0.8
0.6
0.4
0.2 f(A)
0 1, Proposed f(A)
0 0.2 0.4 0.6 0.8 1 1.2 1.4
Figure 6: The function f". Symbols stand for the
measured values in highly resolved simulation and the
dashed line is (25).
A Posteriori Test
To validate the model, a posteriori test has been done
with twodimensional dense fluidized bed. The non
dimensional bed heights with and without model for
coarses meshes are shown in Figure 7 and results are
promising. Herein, we present results briefly because the
nondimensional mesh size, A/(TS') 9, is questionable
parameter and more detailed analysis have to be done by
different type of particles and different bed widths to be
completely sure that this parameter represents physical
interpretations.
IP
( S) 9
Conclusions
f (a)= 1.75
The global dynamics of fludized beds are strongly de
pendent on the mesh size and simulations conducted by
coarse meshes can not predict accurately the statistics
n a
.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
References
Andrews A.T., Loezos P.N. and Sundaresan S., 3D
Coarse Grid Simulation of GasParticle Flows in Ver
tical Risers, Ind. Eng. Chem. Res., Vol. 44, pp. 6022
6037, 2005
Balzer G., Boelle A., Simonin O., Eulerian GasSolid
Flow Modelling of Dense Fluidized Bed, FLUIDIZA
TION VIII, Proc. International Symposium of the Engi
neering Foundation, pp 409418, 1995.
K. Agrawal, P. N. Loezos, M. Syamlal, and S. Sundare
san, The Role of Mesoscales Structures in Rapid Gas
solid Flows, J. Fluid Mech. 445, 151, 2001.
J. De Wilde, Reformulating and Quantifying the Gener
alized Added Mass in Filtered Gassolid Flow Models,
Plws. Fluids 17, 113304, 2005.
O. Simonin, Continuum Modelling of Dispersed Two
Phase Flows, in Combustion and Turbulence in Two
Phase Flows, Lecture Series 199602, von Karman In
stitute for Fluid Dynamics, Rhode Saint Gendse (Bel
gium), 1996.
W. Wang, J. Li, Simulation of Gassolid Two Phase Flow
a Multiscale CFD Approach Extension of the EMMS
model to the subgrid level, Chem. Eng. Sci. 62, 208,
2007.
G. J. Heynderickx, A. K. Das, J. De Wilde, and G. B.
Marin, Effect of Clustering on Gassolid Drag in Dilute
Twophase Flow, Ind. Eng. Chem. Res. 43, 4635, 2004.
Y. Igci, A. T. Andrews IV, S. Sundaresan, S. Pannala,
T. O'Brien, Filtered Twofluid Models for Fluidized Gas
Particle Suspensions, AIChE J. 54, 1431, 2008.
C. Y. Wen, Y. H.Yu, Mechanics of Fluidization, Chem.
Eng. Prog. Sym. Ser., 62, 100113, 1966.
J. W.Deardoff, On the Magnitude of the Subgrid Scale
Eddy Coefficient, J. Comp. Phy. 7, 120133, 1971.
Moreau, M., Simonin, O., Bedat, B., Development
of GasParticle EulerEuler LES Approach: A Priori
Analysis of Particle Subgrid Models in Homogeneous
Isotropic Turbulence, Flow Turbulence Combust. 84,
295324, 2010.
Gobin A., Neau H., Simonin O., Llinas J.R., Reiling
V, Selo J.L., Fluid Dynamic Numerical Simulation of
a Gas Phase Polymerisation Reactor, Int. J. for Numeri
cal Methods in Fluids 43, 11991220, 2003.
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0.5 0.6 0.7 0.8
1 1.1
Figure 7: The influence of grid size on the height of
dense fluidized bed: with model a, without model O.
quantities in the beds. This poor predictions of gassolid
flows in beds can be attributed to lack of an appropri
ate model of the drag term for coarse mesh simulations.
For this propose, the mesh refinement studies were car
ried out for twodimensional dense fludized and three
dimensional circulating fludized bed to get mesh inde
pendent results. Then, these results are used to perform
a priori analysis on drag term and the subgrid contribu
tion of drag is presented with plwsical explanation. It is
stated that the filtered momentum equation can be com
puted on "coarse simulations" but must take into account
the particle to fluid drift velocity due to the subgrid cor
relation between the local fluid velocity and the local
particle volume fraction. We propose a model for the
subgrid drift velocity which can be written the differ
ence between the averaged of gas velocity weighted by
solid volume fraction and the averaged of gas velocity
weighted by gas volume fraction. This model is a func
tion of the filtered particle volume fraction and the non
dimensionless mesh size. A posteriori test is performed
on twodimensional dense fluidized bed and results have
good agreement. However, the nondimensionless mesh
size is the questionable part of the approach and it has to
be investigated by different flow configurations. Addi
tion that, closures are still missing for terms in (12) and
(13) and it will be the object of further work.
Acknowledgements
This work was granted access to the HPC resources
of CINES under the allocation 2010026012 made by
GENCI (Grand Equipement National de Calcul Inten
sif) and of CALMIP under the allocation P0 111 (Calcul
en MidiPyr~ndes ).
