Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 6.4.3 - Quadrature-based Third Order Moment Models for the Simulation of Particle Laden Flows
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00156
 Material Information
Title: 6.4.3 - Quadrature-based Third Order Moment Models for the Simulation of Particle Laden Flows Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Le Lostec, N.
Villedieu, P.
Simonin, O.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: quadrature
method of moments
fluid-particle flow
 Notes
Abstract: We present here a new method of moments for the numerical simulation of particle-laden flows. The closure needed in Eulerian methods relies on writing the kinetic descriptor, the velocity destribution function, as a sum of delta-functions instead of the one-delta-function or close-to-Maxwellian assumption in existing methods. The closure velocity distribution function parameters are computed from the transported moments using a quadrature method. 2D simulation results are compared to those of a one-delta-function-based Eulerian method and those of a reference Lagrangian simulation, considering only transport and drag of monodisperse solid particles in a Taylor-Green 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 quadrature-based method, compared to the one-delta-function method, allows to consider particle inertia effects in the Taylor-Green 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.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00156
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 643-LeLostec-ICMF2010.pdf

Full Text



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


Quadrature-based Third Order Moment Models for the Simulation of
Particle-Laden 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, fluide-particle flow




Abstract

We present here a new method of moments for the numerical simulation of particle-laden flows. The closure needed in
Eulerian methods relies on writing the kinetic descriptor, the velocity distribution function, as a sum of delta-functions
instead of the one-delta-function or close-to-Maxwellian 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 one-delta-function-based Eulerian method and those of a reference
Lagrangian simulation, considering only transport and drag of monodisperse solid particles in a Taylor-Green 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 quadrature-based method, compared to the one-delta-function method, allows to consider particle
inertia effects in the Taylor-Green 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 fluid-droplet or fluid-particle 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 fluid-particle flows, more precisely on
the treatment of particle velocity dispersion, for solid
particles. Leaving aside polydispersity effects, let us
consider only simple interactions: Stokes-shaped drag
force and elastic particle-particle 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 Stokes-shaped 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
higher-order 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 close-to-Maxwellian (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
quadrature-based method simulation results for aca-
demic unsteady particle-laden flows: a solid-particle-
spreaded Taylor-Green 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 particle-particle 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,
particle-particle and fluid-particle 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 quadrature-based methods presented here re-
quire, in 3-D, 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, particle-particle interactions, particle evapora-
tion and so on, the computational cost increase due to
the quadrature method would be marginal.
Quadrature-based methods evaluations in academic
cases were previously performed considering jet cross-
ing, jet-wall rebound, Taylor-Green 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 quadrature-based velocity dis-
persion with polydisperse modeling was also explored
considering on the one hand particle-particle 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 Taylor-Green 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 (positions-Y 0 velocities-v) (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 1-D 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.
2-D moment equations
Considering the general moment transport equation
given above, in two dimensions and for moments up to
third order ( = {1, vi,. . }), the 2-D 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


Qijk-l vi'j v'dv


An Eulerian model is characterized by its order, which
corresponds to the maximum order of the transported
moments:

first-order model n (number density) and vi
(mean velocity) are transported (equivalent to
pressure-less gas dynamics),

reduced second-order model n, v, and E> Rii
(kinetic energy) are transported (equivalent to
Navier-Stokes equation),

full second-order model n, vi, Rij (kinetic
stresses) are transported (equivalent to Grad's mo-
ment method for rarefied gases).

The most basic Eulerian model (first-order model), also
called bifluid method, assumes that the velocity distri-
bution function is a mean-velocity 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 "sticky-particle" 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
gradient-diffusion 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.

Quadrature-based method

Transported moments
In the quadrature-based moment method we consider
the ten 2-D 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 quadrature-based 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 non-negative 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 4-node quadrature reconstruction of f. Note that
the set of quadrature nodes contains 12 unknowns (i.e.,
four weights and four two-component 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 1-D quadratures and to couple them, with respect
to particle velocity covariance, to form a 2-D 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
two-node 1-D 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 1-D weights and X(j)i the 1-D ab-
scissas, for the node i in j-direction.


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


The 2-D quadrature discrete velocities are those of the
1-D quadrature, combined using a tensor-product 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 1-D quadrature weights a, and 3, that yield the
correct value of the second-order cross-moment M2,.
Defining p = R12/RR22, a straightforward com-
putation gives the 2-D 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 four-node 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 third-order moments to set the
weight equal to zero, and recompute the other weights.
Moreover as the off-diagonal 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 frame-change-based 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
zero-mean-velocity 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 2-D closure velocity distribution in the ro-
tated frame F* is taken as the product of two two-nodes
1-D 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 1-D 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 1-direction 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.

Quadrature-based method: schemes and
algorithm

A first-order finite-volume 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 1-D
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 Runge-Kutta 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

Taylor-Green flow configuration


0 0.2 04 0.6 0.8 1


Figure 1: Taylor-Green fluid flow velocity field.


The Taylor-Green 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 3-D HIT sim-
ulation on a 1283 mesh then the 2-D fluid velocity field
is not divergence free. The 3-D HIT was generated ac-
cording to a Passot-Pouquet 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 2-D 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 Taylor-Green 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 Taylor-Green case and
5122 cells mesh for HIT case.


1 Y I 05
0 0.2 0.4 0.6 0.8
X


Figure 3: Taylor-Green 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 Taylor-Green 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, Taylor-Green 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 Taylor-Green 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 Taylor-Green
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 Taylor-Green 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 Taylor-Green 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 Taylor-Green 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 Taylor-Green 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 second-order 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 non-uniform 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 multi-fluid 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 265-276,
2009.

O. Desjardins, RO Fox, and P. Villedieu. A quadrature-
based moment method for dilute fluid-particle flows.
Journal of Computational Physics, 227(4):2514-2539,
2008.

P. F6vrier, O. Simonin, and K.D. Squires. Partitioning
of particle velocities in gas-solid turbulent flows into
a continuous field and a spatially uncorrelated random
distribution: theoretical formalism and numerical study.
Journal of Fluid Mechanics, 533:1-46, 2005.


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


R. O. Fox. A quadrature-based third-order moment
method for dilute gas-particle flows. Journal of Com-
putational Physics, 227(12):6313-6350, 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 277-288, 2009.

J. He and O. Simonin. Non-equilibrium prediction of the
particle-phase stress tensor in vertical pneumatic con-
veying. ASME-FED, 166:253-263, 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):6448-6472, 2008.

N. Le Lostec, R. O. Fox, O. Simonin, and P. Villedieu.
Numerical description of dilute particle-laden flows by
a quadrature-based moment method. In Proceedings of
the Center for Turbulence Research Summer Program
2008, pages 209-221, 2009a.

N. Le Lostec, P. Villedieu, and O. Simonin. Comparison
between Grad's and quadrature-based 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 non-equilibrium effects in dilute granular
flows. In Proceedings of the Twenty-First International
Symposium on Rarefied Gas Dynamics, volume 1, pages
287-294. C6padubs Editions, July 1999.

L. Schneider, N. Le Lostec, P. Villedieu, and A. Sadiki.
A moment method for polydisperse sprays. Comptes
rendus-Mathimatique, 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:261-272, 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: Taylor-Green 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: Taylor-Green 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: Taylor-Green 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.




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

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