7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Quadraturebased Third Order Moment Models for the Simulation of
ParticleLaden Flows
N. Le Lostec*, P. Villedieu* and 0. Simonint
Office National d'Etudes et de Recherches A6rospatiales (ONERA),
2 Avenue Edouard Belin 31055 Toulouse, France
t University de Toulouse; INPT, UPS; Institut de M6canique des Fluides de Toulouse (IMFT),
All6e du Professeur Camille Soula, 31400 Toulouse, France
CNRS; Institut de M6canique des Fluides de Toulouse (IMFT)
nechtan.lelostec@onera.fr, villedieu@onera.fr and simonin@imft.fr
Keywords: quadrature, method of moments, fluideparticle flow
Abstract
We present here a new method of moments for the numerical simulation of particleladen flows. The closure needed in
Eulerian methods relies on writing the kinetic descriptor, the velocity distribution function, as a sum of deltafunctions
instead of the onedeltafunction or closetoMaxwellian assumption in existing methods. The closure velocity distri
bution function parameters are computed from the transported moments using a quadrature method.
2D simulation results are compared to those of a onedeltafunctionbased Eulerian method and those of a reference
Lagrangian simulation, considering only transport and drag of monodisperse solid particles in a TaylorGreen fluid
flow and in a frozen homogeneous isotropic turbulence fluid flow. For a range of Stokes number the particle velocity
distribution function is far from equilibrium and particle trajectory crossing is an important feature.
We find that the quadraturebased method, compared to the onedeltafunction method, allows to consider particle
inertia effects in the TaylorGreen case and slightly extends the particle Stokes number range that can be simulated
with good agreement with the lagrangian reference particle density map in the homogeneous isotropic turbulence case.
Introduction
Unsteady fluiddroplet or fluidparticle flows are in
volved in both industrial systems like fuel combustors
in piston and jet engines, coal boilers, chemical reactors
or rocket engines and in natural phenomenon like dust
storms or volcanic plumes. In particular due to energy
and pollution issues, multiphase flow modeling and sim
ulation are a growing subject of applied and fundamental
studies.
In this work, we focus on particle kinematic descrip
tion in turbulent fluidparticle flows, more precisely on
the treatment of particle velocity dispersion, for solid
particles. Leaving aside polydispersity effects, let us
consider only simple interactions: Stokesshaped drag
force and elastic particleparticle collisions. Then the
three relevant adimentionalized numbers are particle
Reynolds number, dynamic Stokes number build on
fluid large eddy characteristic time and collision Stokes
number defined in terms of mean collision time. As the
particle Reynolds number only acts as a correction in
the Stokesshaped drag, the two Stokes numbers are the
most important. For a dynamic Stokes number much
smaller than collision Stokes number, the particles ve
locity will stick to the local fluid velocity and the par
ticle velocity distribution will be unimodal, because of
drag force damping particle agitation faster than it is
produced in the flow and isotropized by collisions. At
the opposite, when dynamic particle Stokes number is
much greater than collision Stokes number, because of
predominant collisions the particle velocity distribution
is Maxwellian. Between these two limits the velocity
distribution function can have any shape.
While in Lagrangian simulations the trajectory of
each numerical particle is computed, the Eulerian ap
proach characterizes the particles on a statistical ba
sis. The kinetic equation gives the evolution of the ve
locity distribution function whose moments are the de
sired macroscopic quantities. As the kinetic equation
is not solved directly but a finite set of moment trans
port equations are calculated, each of which involving
higherorder moments and unknown source terms, clo
sures have to be made for the unknown required mo
ments and source terms.
Two existing Eulerian methods give closures appro
priate for unimodal (standard Eulerian method, referred
to as bifluid method) and closetoMaxwellian (Grad's
method) velocity distribution function. But they obvi
ously fail in intermediary cases between predominant
drag force and predominant collisions. Quadrature
based methods introduce a new kind of closure, done
directly on the particle velocity distribution function: it
consists in writing the velocity distribution function as a
weighted sum of delta functions whose weights and ab
scissas are computed from the known moments. Then
the closure distribution function is used to compute the
unknown moments and source terms. Furthermore the
method can a priori be extended to any arbitrary order,
the closure problem being written as a quadrature inver
sion problem.
Due to the discrete nature of particles in Lagrangian
methods, no limitation arise on the shape of the particle
velocity distribution. However the number of particles
needed to reduce the statistical noise in Eulerian statis
tics increases with the particle load and the unsteadiness
of the flow and so does the computational cost. More
over the coupling between the continuous phase and the
particulate phase is difficult in Lagrangian methods. At
the opposite, Eulerian methods computational cost does
not depend on particle density and remain constant for
any load as the particles are seen as a continuous me
dia. Furthermore the coupling between the two phases
is more easily done. Consequently Eulerian methods are
better suited to the simulation of unsteady dense fluid
particle flows.
The object of our work is to compare bifluid (one
delta velocity distribution function) method and a
quadraturebased method simulation results for aca
demic unsteady particleladen flows: a solidparticle
spreaded TaylorGreen fluid flow and a frozen homo
geneous isotropic turbulence fluid flow. The Taylor
Green flow case was introduced in de ( Ii.,isi.i.i li et al.
(2007) as a test case for polydisperse sprays dynamic
simulation and is of much interest because the fluid flow
seen by the particles is unsteady and it presents particle
interaction with fluid vortices which are the basic struc
tures that turbulence is made of. Homogeneous isotropic
turbulence fluid flow is closer to real turbulent flows then
its study is a step towards real cases simulation. Refer
ence results are provided by a Lagrangian simulation.
We focus on particle velocity dispersion, in particular
on the ability of Eulerian methods to handle far from
equilibrium particle velocity distribution. Then in the re
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
mainder the particulate phase is monodisperse and we do
not consider neither particleparticle collisions nor cou
pling between the two phases. Note that these choices of
modelisation are not limits of Eulerian modeling, no dif
ficulty arise adding for example polydispersity effects,
particleparticle and fluidparticle interactions to the fol
lowing kinetic models. The dynamic Stokes number will
range in [0.1; 10] where noticeable patterns in the parti
cle density have been observed in Lagrangian simulation
with important particle trajectory crossing involved.
The quadraturebased methods presented here re
quire, in 3D, to solve three times as much equation
as the bifluid method. Because these two methods in
volve computations of the same complexity, a three
times higher computational cost can be expected. Note
that this cost only concerns the inertial part of the par
ticle movement. Then in a more general simulation, in
cluding solving the fluid flow, coupling between the two
phases, particleparticle interactions, particle evapora
tion and so on, the computational cost increase due to
the quadrature method would be marginal.
Quadraturebased methods evaluations in academic
cases were previously performed considering jet cross
ing, jetwall rebound, TaylorGreen flow, Riemann prob
lem, impinging particle planes and vertical channel flow
in Desjardins et al. (2008); Fox (2008); Le Lostec et al.
(2009a,b). Coupling of quadraturebased velocity dis
persion with polydisperse modeling was also explored
considering on the one hand particleparticle collisions
and evaporation in de Chaisemartin et al. (2009); Fr6ret
et al. (2009), on the other hand evaporation and particle
wall interaction (splashing) in Schneider et al. (2009,
2010). In this work we present results obtained with
a modified quadrature that overcomes some drawbacks
of the quadrature method presented in Le Lostec et al.
(2009a,b) and for an extended range of Stokes number
in the TaylorGreen fluid flow test case.
Moment methods
Kinetic description
The kinetic descriptor is the velocity distribution
function f, which is the particle number density in phase
space (positionsY 0 velocitiesv) (see Chapman and
Cowling (1970)). Its transport equation is a Boltzmann
like equation:
8f(t, x, v)
(t, x + v gradx f(t, x, v)
Ot
+div, (f(t,x,v) a) 0
where a is the acceleration the particles experience (we
will provide its expression later).
The Eulerian macroscopicc) descriptors are the mo
ments of the velocity distribution function, defined, for
a 1D velocity space, as
J J= vqf(t, x,v)dv
q 0, particle number density (n),
q = 1, n times particle mean velocity,
q = 2, n times particle mean square velocity.
Eulerian methods solve transport equations for such mo
ments.
Moment transport equations
The moment transport equations come from integrat
ing over velocity space the kinetic equation multiplied
by powers of the velocities components, denoted y. Its
general form is
] a fdv
t [j f dvj + x, [j \if dvj j f dv
Note that the moment transport equations can be writ
ten for either the moments or the central moments corre
sponding to integral of total velocity v or velocity fluc
tuation v' v v with respect to the mean velocity
v = f vf dv. The equations for moments and central
moments are strictly equivalent. Any total moment can
be expressed in terms of lower order central moments.
The model for drag acceleration is the Stokes drag
with fixed dynamic response time. Thus, a
T 1(Vf v) where Tp = 2ppr2/(9pf), Vf is the lo
cal fluid velocity, v the particle velocity, pp the particle
mass density, r the particle radius, and pf the fluid dy
namic viscosity.
2D moment equations
Considering the general moment transport equation
given above, in two dimensions and for moments up to
third order ( = {1, vi,. . }), the 2D moment
equations are
93\ 3M1
/,+ V 0
at Ox,
3M1 39'.
3t 3Xj
kt r 2
+ Aij
at O,
k + j k A 3
9t 9xi ijk
where acceleration terms A, A? A3k can be expressed
in terms of the moments up to third order (i.e., they are
closed).
Existing moment methods
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The existing Eulerian methods developed in the
frame of classical kinetic theory of gases are using trans
port equations for particle number density:
n f dv
and central velocity moments:
Rij 1 'f dv
n vi j
S ijk 1f '.'fdv
Sijk a vi vj
Qijkl vi'j v'dv
An Eulerian model is characterized by its order, which
corresponds to the maximum order of the transported
moments:
firstorder model n (number density) and vi
(mean velocity) are transported (equivalent to
pressureless gas dynamics),
reduced secondorder model n, v, and E> Rii
(kinetic energy) are transported (equivalent to
NavierStokes equation),
full secondorder model n, vi, Rij (kinetic
stresses) are transported (equivalent to Grad's mo
ment method for rarefied gases).
The most basic Eulerian model (firstorder model), also
called bifluid method, assumes that the velocity distri
bution function is a meanvelocity delta function: f
n 6 (v v). Thus all particles at the same location have
the same velocity, equal to the mean local velocity and
there is no particle agitation. For this reason, it is also
referred to as the "stickyparticle" model.
For existing Eulerian models of second or higher or
der, the velocity distribution is assumed to be close to
an equilibrium distribution corresponding to the local
Maxwellian. Then, Enskog's or Grad's method may
be used to close the possible collision terms and derive
gradientdiffusion model for the unknown higher order
moments appearing in the computed transport equations.
See He and Simonin (1993), Sakiz and Simonin (1999),
Simonin (2000) and Kaufmann et al. (2008) for the de
tails of these models derivation.
Quadraturebased method
Transported moments
In the quadraturebased moment method we consider
the ten 2D moments up to third order:
(;l", Mf{, M2, M1Af1, M122, ,
1f M 12 ).
where
I" = / (t,x,v) dv, I = / v f(t,x,v) dv,
3 [ .. f(t, x, v) dv,
Mi= J f(t, x,v) dv
Quadrature closure
In quadraturebased closures, the particle velocity dis
tribution function f is reconstructed from the moments
in the form
N
f(v) a(v U,)
al
where no(x, t) are nonnegative weights and Ua(x, t)
are the velocity abscissas. Let [(na, Ua)] with a E
(1,..., 4) denote the set of weights and abscissas for
the 4node quadrature reconstruction of f. Note that
the set of quadrature nodes contains 12 unknowns (i.e.,
four weights and four twocomponent velocity vectors).
To find these unknowns, we work with the velocity mo
ments up to third order defined above, which are related
to the quadrature weights and abscissas by
1" Z=n, Mi:
a=l
4
E ns Uai
a=1
'. = noV Uai U Mjk n UniUa Uaks
a=1 a1=
(2)
Quadrature inversion is the algorithm for computing
weights and abscissas from moments. The inverse op
eration, finding moments from weights and abscissas,
is given in Eq. (2). In general, it will not be possible
to represent all possible moment sets using weights and
abscissas.
Direct inversion
The principle behind direct inversion is to compute
two 1D quadratures and to couple them, with respect
to particle velocity covariance, to form a 2D quadra
ture, as reported in Le Lostec et al. (2009a,b). In
this inversion method only the diagonal third order mo
ments A311 and M222 are used. We start by defin
ing two subsets of moments: (II", M1, M21, MA311)
and ( 2", M1, Mff, M22). Each subset yields a
twonode 1D quadrature approximation using the
procedure described in Desjardins et al. (2008):
[(al, X(1)1), (02, X(1)2)] and [(/3i, X(2)1), (/32, X(2)2)].
ai and 3, being the 1D weights and X(j)i the 1D ab
scissas, for the node i in jdirection.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The 2D quadrature discrete velocities are those of the
1D quadrature, combined using a tensorproduct so that
we have only two discrete velocities per direction but
four velocity vectors:
X( )
X(2)2 '
(X(1)2)
X(2)2
U (X(1)1 ), U2
Xv(2)1 )
U3 ( 1)2) U4
X(2)1 '
The coupling is done by computing the weights n,
from the 1D quadrature weights a, and 3, that yield the
correct value of the secondorder crossmoment M2,.
Defining p = R12/RR22, a straightforward com
putation gives the 2D quadrature weights:
S ,/3 aa23P+1 32
aO )
2 (
1i/ 2 ai02/31/32
P+ a ',I"
r J2 131
a131 J
S= P+ \  "
Then the fournode quadrature approximation is defined
as
[(ni, U1), (n2, U2), (n3, U3), (n4, U4)].
Yet this inversion method does not ensure positive
weights. Thus, if one weight happens to be negative, we
change the value of the thirdorder moments to set the
weight equal to zero, and recompute the other weights.
Moreover as the offdiagonal third order moments are
not used, their values are not given back by the quadra
ture, the resulting quadrature method is not tensorially
consistent.
Frame change inversion
For the framechangebased inversion method the
central moments up to order three are considered:
(n, v, v, R11, R22, R12,
S111, S222, S112, S122)
The inversion procedure consists in computing the
quadrature's weights and abscissas in a frame F* where
quadrature computation is simplified, following the
method introduced in Fox (2008). However the inver
sion procedure proposed here is different in that the
third order moments used are the contracted moments
T S111 + S122 and T2 = S112 + S222 instead
of the whole set of third order moments (S111, S112,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
S122 and S222) then the number of quadrature param
eters (weights and abscissas) is equal to the number of
transported moments.
The first step is the normalization of central moments
setting n 1 and the mean velocities vi equal to zero,
which can be seen as a translation of the moments to the
zeromeanvelocity point.
Then the moments are rotated into the frame F*
where the second order moments tensor R* is diagonal,
by computing the eigenvalues and eigenvectors of R and
applying the transformation to the third order contracted
moments vector T. Under the Ii\i s1l, c\i, that the veloc
ity distributions in x: and x* directions are independent
(which is reasonable considering that the covariance R*2
is zero), the 2D closure velocity distribution in the ro
tated frame F* is taken as the product of two twonodes
1D distributions:
f* (a 8(vi X1)1) + (1 a) 6(v X()2))
x (/3 8(2 X(2)) + (1 /3) 8(V2 X2)2))
This closure function implies that the third order off
diagonal moments S~22 and S112 are zero, then in the
rotated frame F* the third order contracted moments are
equal to the third order diagonal moments: Ty = S1i
and T* = S2,,. Thus the 1D quadrature inversion
method of Desjardins et al. (2008) can be used to com
putes the weights a and 3 and the abscissas X* from
the two sets of central normalized rotated moments:
{n = 1, O O, RX, T, Sii=1,2.
The final weights and abscissas are obtained devel
opping the expression of the closure velocity distribu
tion function, rotating the velocity vectors back to the
initial cartesian frame and adding the initial mean ve
locity components. Compared to the direct inversion
method, the coupling between 1direction quadratures is
done through velocities instead of weights. Denoting 0
the angle of the rotation from the initial cartesian frame
to the frame F*, the closure velocity distribution func
tion weights are:
nj n a 3 n2 = n a(1 /3)
ns = n (1 a) 3 n4 = n (1 a) (1 3)
and the abscissas:
= (v + cos 0 X() sin 0 X)
Ssin0 X2 +cosOX
2 + sin 0 X1) + cos 0 X2)1
71 + cos e x* sin e x*
U2 (1)1 (2)2
72 + sin 0 X) + cos 0 X2)
7 1 + cos 0 X)2 sin 0 X*
S(1)2 (2)1
o2 + sin 0 X*()2 + cos 0 X*&)
U4 (1) (2)2
72 + sin 2 X)21) + COS 0 X(*2)
Note that this inversion method gives positive
weights, due to a and 3 taking their values in [0; 1], and
is tensorially consistent, all the third order moments be
ing used.
Quadraturebased method: schemes and
algorithm
A firstorder finitevolume method gives the spatial dis
retisation while fluxes are written according to a kinetic
scheme. A fractional step algorithm allows to treat sepa
rately transport and drag terms, similarly to what is done
in Desjardins et al. (2008). Time evolution is explicitly
computed.
The moment equations to solve, considering a 1D
case for the sake of clarity, can be written as:
aW aH(W)
t ax A
Ot ox
where
M1 ,H M2
W M2 ,H(W)= M3 ,A
M3 M4
In the first step of the algorithm, transport only is con
sidered, and gives intermediary moments W* from the
current iteration initial moments Win in the cell denoted
W Wn A [G (W[ ni W 1) (W 1 W]
where At and Ax are time and space steps, G being the
numerical flux function. Having
K(,) '2
\ 3
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
such that
W = K(v) f dv and H(W) = vK(v) fdv
we can write the numerical flux function as:
G (W, W)= l (v + v) K(v) fw, dv
+ (v v) K(v) fw, dv
Then the numerical flux are expressed in terms of the
closure velocity distribution function (see Eqn. (1), fw,
for cell I and fw, for cell r), whose weights nji and
abscissas Ui (j being the quadrature node subscript)
are computed from the moments Wi.
The second step of the algorithm is the drag effect
computation. We proceed in the same way as Desjardins
et al. (2008), the drag force results in a relaxation of the
quadrature abscissas to the fluid velocity, the weights not
being changed. The relaxation applies to quadrature ab
scissas Ujp* computed from the intermediary moments
Ww:
nj, =nj,i
{uj 1 (e At AtM
U ~ U y T) +V1,4 1 exp y
The final moment for the current iteration W["1 are
then computed from the weights nf 1 and abscissas
Uj1 according to Eqn. (2).
Lagrangian method equations and scheme
Lagrangian method consists in solving the following set
of equations for the position Xi and velocity Vi of each
numerical particle denoted i:
dXi Vi
di= Vi
dt
dVi
dt
/i Vf(Xi)
T_
The scheme is an explicit fourth order RungeKutta with
a time step At = T/10. A volume average gives the
Eulerian statistics maps according to the moments defi
nitions from the Lagrangian particle velocity distribution
function:
flag = wi (x Xi)
Flow configuration and results
TaylorGreen flow configuration
0 0.2 04 0.6 0.8 1
Figure 1: TaylorGreen fluid flow velocity field.
The TaylorGreen fluid flow is made of periodic vor
tices with alternate rotating direction (Fig. 1). The fluid
velocity is defined as follows:
Vf, = sin (27x) cos (27y)
Vf,y = sin (27y) cos (27x)
with x, y E [0; 1].
The fluid characteristic time 7f is the ratio of the cell
size to the maximum velocity of the fluid then f = s.
Initially the particle density is homogeneous and the
particles velocity is set to zero.
Frozen homogeneous isotropic turbulence flow
configuration
1
06 .' . .. .. ~',.;;: : . .
*" IL I ' *
0.4
0.....
0 02 04 0.6 0.8 1
Figure 2: HIT fluid flow velocity field.
The fluid flow for frozen homogeneous isotropic tur
bulence (HIT), see Fig. 2, is a slice of a 3D HIT sim
ulation on a 1283 mesh then the 2D fluid velocity field
is not divergence free. The 3D HIT was generated ac
cording to a PassotPouquet spectrum, the most ener
getic eddy adimentionalized length is L, = 2, the tur
bulent Reynolds number Ret = 100, the integral length
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
scale based Reynolds number Re, 48.5, the turbu
lent kinetic energy k 1.5.10 2 and the turbulent dis
sipation rate 6.09.10 4. A linear interpolation gives
the fluid velocity in each cell of greater than 1282 cells
meshes for Eulerian simulations or at particle position in
Lagrangian simulations.
The fluid characteristic time is the ratio of the greatest
vortex size L present in the flow, around 0.25m here, to
the 2D rms fluctuating velocity on the whole domain,
then Tf L/ ((v V)2 2.5s.
Initially the particle density is homogeneous and each
particle velocity is equal to the local fluid velocity.
Simulation parameters
Simulations results were analyzed for times t = f,
t 2 Tf and t = 5 Tf. Because fluid velocity field struc
tures have a limited life time of order Tf, whe choose
the fluid characteristic time Tf as time scale and re
stricted simulations to times up to t 5 Tf.
The particle Stokes number St is the inertia parameter:
St = where Tp is the particle dynamic response time
defined above. Results are presented for Stokes ranging
in [0.1; 10].
For Lagrangian simulations the numerical particles
number is 108 and the corresponding moments are com
puted on a 2002 cells mesh for TaylorGreen case and
on a 5122 cells mesh for HIT. Numerical convergence is
achieved.
Note that because we consider here no particle
particle interactions and no influence of the particles on
the fluid flow, the kinetic equation is linear (this is the
case, for example, of dilute flows). Then as the fluid
flow is given a priori, it is equivalent in terms of the
particle velocity distribution function to compute an en
semble average of the volume average in each cell over
Ns simulations with Np numerical particles or to do one
simulation for Ns x Np numerical particles and volume
averaging in each cell. Indeed as the particle velocity
distribution function is a probability density function it
can represent any number of real particles. However, to
compute it, it is necessary to average over enough nu
merical particles to guarantee the numerical convergence,
i.e. a sufficient number of numerical particles in each
cell to reduce the statistical noise to an acceptable level.
For Eulerian simulations a convergence study shows
that increasing the mesh size above 5002 cells only re
duces numerical diffusion then sharpens the structures
observed on particle number density map, as illustrated
on Fig. 3. Then as no more information is obtained
refining the mesh from 5002 cells size while the compu
tational cost quickly increases, all Eulerian simulation
are based on 5002 cells mesh for TaylorGreen case and
5122 cells mesh for HIT case.
1 Y I 05
0 0.2 0.4 0.6 0.8
X
Figure 3: TaylorGreen fluid flow. Particle number den
sity maps. Frame change quadrature. Stokes St 10,
time t 5 Tf. Top: 5002 cells mesh, bottom: 10002
cells mesh.
Results
Only the most illustrative results are presented here.
We compare the particle number density maps given by
the different methods (Lagrangian, bifluid, direct inver
sion quadrature and frame change inversion quadrature)
which are a good iIheL',i of each method performance
because the errors on all moments induce error on the
particle number density.
For TaylorGreen flow case with small inertia parti
cles (St = 0.1), due to their small inertia, all the parti
cles velocity stick to the local fluid velocity then when
reaching the vortex edges, all the particles remains in the
vortex delimiting cells. In this particular case, the bi
fluid method hip' 1lc ii of only one velocity for all the
particles in each cell is verified then the bifluid method
particle density map is similar to the Lagrangian ref
erence. Quadrature based methods also perform well
for all times. The particle density maps for all meth
ods show two lines of high particle density similar to the
bifluid method result in Fig. 4, top right.
In Fig. 4, TaylorGreen flow case with small inertia
particles (St 0.5) at long time (t 5 Tf), we observe
stripes of particles moving towards the same direction.
These stripes cross each other, involving multimodal
particle velocity distribution. Particles have small but
sufficient inertia not to follow the fluid streamlines then
the bifluid method fails while direct quadrature method
particle number density map only have some features
of the lagrangian reference such as the empty squares
centered on the vortices. The frame change quadra
ture method performs much better, the particle stripes
are present but widened compared to the lagrangian due
to numerical diffusion and even the background lower
density structures are evoqued. At smaller times the
two quadrature methods results are very similar to La
grangian.
For higher Stokes number in the TaylorGreen case,
as the particles velocity is different from the fluid ve
locity, the bifluid method always fails and produces the
same pattern of particles concentrated at the edges of the
vortices then we will not discus its performance.
Fig. 5 shows the simulation results for TaylorGreen
flow case with larger inertia particles (St 5) at long
time (t 5 Tf). Lagrangian reference particle number
density map has a complex structure that is not matched
by the quadrature methods except around the vortex cen
ters for lines of high density on the direct quadrature
method density map and for some similarities of the
square structures and high density lines on the frame
change quadrature map. For times smaller than t 2 Tf,
quadrature methods particle density maps are very close
to the lagrangian.
The highest inertia particles TaylorGreen case anal
ysed, St=10, at long time t 5 Tf, presented on Fig.
6 exhibits a very complex structure which is quite well
predicted by the quadrature methods, with only small
distortions. The same conclusions are valid for smaller
times.
Note that for Stokes greater than a value between
1 5, the particle number density map of frame change
quadrature method show a small amplitude noise pertur
bation which we think is related to a numerical instabil
ity. This point is the object of current work.
HIT fluid flow with small inertia (St = 0.1) particle
and medium time (t 2 Tf) results are presented on Fig.
7. We can see that bifluid and quadrature methods par
ticle density maps match the Lagrangian reference: the
particles concentrate on high density thin lines which are
widened due to numerical diffusion in Eulerian methods
results.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
For Stokes below 1, even at long times t < 5 Tf, the
bifluid method provides particle density maps in good
agreement with the Lagrangian reference. However for
Stokes above 1 this method is no more able to give re
sults comparable to the Lagrangian reference, as can be
seen on Fig. 8. At the contrary, quadrature methods
predict well the complicated structures visible on the
Lagrangian reference particle density map with small
perturbations appearing on the frame change quadrature
map due to the supposed numerical instability.
High inertia particles (St > 5) in HIT fluid flow do not
concentrate in some noteworthy pattern, see Fig. 9. Eu
lerian methods particle density map show some lines of
high density which stand for nothing on the Lagrangian
reference map. Note the noise on frame change quadra
ture result. Even for small times, none of the three Eule
rian methods is able to predict a density map similar to
the Lagrangian reference.
Conclusions
The comparison of quadrature based methods perfor
mance to those of the widespreaded bifluid and of ref
erence Lagrangian is promising. The quadrature based
methods extends the range of particle inertia that can be
simulated with good agreement with the reference, by
a great step in the case of TaylorGreen flow and by a
smaller one for the HIT fluid flow case.
The frame change based quadrature method presented
here overall performs better than the direct inversion
quadrature method, particularly for medium inertia par
ticles in the TaylorGreen fluid flow case where the
stripes are well predicted. It also gives a density map
much closer to the Lagrangian reference than the frame
change inversion method of Fox (2008) for Stokes equal
to 1 at long times (t 10Tf), result not presented here
(see Le Lostec et al. 21' I' ., for illustration). There we
see that the tensorial consistency, that is to say consider
ing and giving back the right value of all the third order
moments, under their contracted form here, instead of
only the diagonal third order moments, brings in apre
ciable improvement in the particle density map.
A secondary conclusion is that the TaylorGreen fluid
flow is a very stressing test case for particle velocity
dispersion models, involving complex patterns in par
ticle number density map with far from equilibrium par
ticle velocity distribution. Surprisingly, compared to
the HIT simulation results presented in Desjardins et al.
(2008), the bifluid method performance for HIT fluid
flow is quite good for small to medium particle inertia.
Moreover, at the opposite of expectations expressed in
Le Lostec et al. (2009b), quadrature based methods do
not bring noticeably better results in the HIT fluid flow
case for high inertia particles compared to the bifluid
method.
Further perspectives to this work can be splitted in two
complementary directions. On the one hand doing com
parisons with more evolved Eulerian methods, for other
cases including coupling between the two phases, for ex
ample comparisons with a full secondorder model, spe
cially for large inertia particles with a velocity distribu
tion closing to the equilibrium. On the other hand more
detailed comparison and analysis, examining all the
transported moments and terms of the transport equa
tions versus Lagrangian reference, similarly to what was
done in Kaufmann et al. (2008), would provide valuable
informations on the quality of the different moments and
terms modelisation and allow to determine which are the
most important for an accurate modeling of the particu
late phase.
Acknowledgements
The authors are grateful to E. Masi for providing the HIT
fluid flow velocity field and its characteristics.
References
S. Chapman and TG Cowling. The mathematical theory
of nonuniform gases: an account of the kinetic theory
of viscosity, thermal conduction and diffithion in gases.
Cambridge University Press, 1970.
S. de Chaisemartin, F. Laurent, M. Massot, and
J. Reveillon. Evaluation of Eulerian multifluid ver
sus Lagrangian methods for the ejection of polydisperse
evaporating sprays by vortices. Journal of Computa
tional Physics, submitted for publication, 2007.
S. de Chaisemartin, L. Fr6ret, D. Kah, F. Laurent, Fox R.
O., J. Reveillon, and M. Massot. Turbulent combustion
of polydisperse evaporating sprays with droplet cross
ing: Eulerian modeling and validation in the infinite
knudsen limit. In Proceedings of the Center for Turbu
lence Research Summer Program 2008, pages 265276,
2009.
O. Desjardins, RO Fox, and P. Villedieu. A quadrature
based moment method for dilute fluidparticle flows.
Journal of Computational Physics, 227(4):25142539,
2008.
P. F6vrier, O. Simonin, and K.D. Squires. Partitioning
of particle velocities in gassolid turbulent flows into
a continuous field and a spatially uncorrelated random
distribution: theoretical formalism and numerical study.
Journal of Fluid Mechanics, 533:146, 2005.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
R. O. Fox. A quadraturebased thirdorder moment
method for dilute gasparticle flows. Journal of Com
putational Physics, 227(12):63136350, 2008.
L. Fr6ret, F. Laurent, S. de Chaisemartin, D. Kah, R. O.
Fox, P. Vedula, J. Reveillon, O. Thomine, and M. Mas
sot. Turbulent combustion of polydisperse evaporating
sprays with droplet crossing: Eulerian modeling of col
lisions at finite knudsen and validation. In Proceedings
of the Center for Turbulence Research Summer Program
2008, pages 277288, 2009.
J. He and O. Simonin. Nonequilibrium prediction of the
particlephase stress tensor in vertical pneumatic con
veying. ASMEFED, 166:253263, 1993.
A. Kaufmann, M. Moreau, O. Simonin, and J. Helie.
Comparison between Lagrangian and mesoscopic Eu
lerian modelling appraoches for inertial particles sus
pended in decaying isotropic turbulence. Journal of
Computational Physics, 227(13):64486472, 2008.
N. Le Lostec, R. O. Fox, O. Simonin, and P. Villedieu.
Numerical description of dilute particleladen flows by
a quadraturebased moment method. In Proceedings of
the Center for Turbulence Research Summer Program
2008, pages 209221, 2009a.
N. Le Lostec, P. Villedieu, and O. Simonin. Comparison
between Grad's and quadraturebased methods of mo
ments for the numerical simulation of unsteady particle
laden flows. In Proceedings ofASME 2009 Fluid Engi
neering Division Summer Meeting, 2009b.
M. Sakiz and O. Simonin. Numerical experiments and
modelling of nonequilibrium effects in dilute granular
flows. In Proceedings of the TwentyFirst International
Symposium on Rarefied Gas Dynamics, volume 1, pages
287294. C6padubs Editions, July 1999.
L. Schneider, N. Le Lostec, P. Villedieu, and A. Sadiki.
A moment method for polydisperse sprays. Comptes
rendusMathimatique, 2009.
L. Schneider, N. Le Lostec, P. Villedieu, and A. Sadiki.
A moment method for splashing and evaporation pro
cesses of polydisperse sprays. International Journal of
Multiphase Flow, 36:261272, 2010.
O. Simonin. Statistical and continuum modelling of
turbulent reactive particulate flows, part 1: Theoretical
derivation of dispersed Eulerian modelling from prob
ability density function kinetic equation. Theoretical
and Experimental Modeling of Particulate Flows, Lec
ture Series, von Karman Institute, 6, 2000.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
2
0 0.2 0.4 0.6 0.8
2 .
0.4
Figure 4: TaylorGreen fluid flow particle density map. Stokes 0.5, time 5 Tf. Top left: lagrangian, top right: bifluid
method, bottom left: direct quadrature method, bottom right: frame change quadrature method.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
2
0 0.2 0.4 0.6 0.8
Figure 5: TaylorGreen fluid flow particle density map. Stokes 5, time 5 Tf. Top left: lagrangian, top right: bifluid
method, bottom left: direct quadrature method, bottom right: frame change quadrature method.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.2 0.4 0.6 0.8
0.2 0.4 0.6 0.8
2
0 0.2 0.4 0.6 0.8
Figure 6: TaylorGreen fluid flow particle density map. Stokes 10, time 5 Tf. Top left: lagrangian, top right: bifluid
method, bottom left: direct quadrature method, bottom right: frame change quadrature method.
n
e4N.. I
)0
I I ____ L
r
 i
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 7: THI fluid flow particle density map. Stokes 0.1, time 2 Tf. Top left: lagrangian, top right: bifluid method,
bottom left: direct quadrature method, bottom right: frame change quadrature method.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
4
3.5
2.5
2
1.5
1
0.5
Figure 8: THI fluid flow particle density map. Stokes 1, time 5 Tf. Top left: lagrangian, top right: bifluid method,
bottom left: direct quadrature method, bottom right: frame change quadrature method.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
4
35
2.5
2
1.5
1
0.5
4
3.5
2.5
2
1.5
1
0.5
Figure 9: THI fluid flow particle density map. Stokes 5, time 2 Tf. Top left: lagrangian, top right: bifluid method,
bottom left: direct quadrature method, bottom right: frame change quadrature method.
