7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Eulerian models and threedimensional numerical simulation of polydisperse
sprays
L. Fr6ret*, S. de Chaisemartint, J. Reveillont, F. Laurent and M. MassotO
Laboratoire d'Imagerie Parametrique, Universit6 Paris VI, France
t IFP, Reuil Malmaison, France
CORIA, Saint Etienne du Rouvray, France
SLaboratoire EM2C, Ecole Centrale Paris, ChatenayMalabry, France
lucie.freret@upmc.fr, stephane.dechaisemartin@ifp.fr, julien.reveillon@coria.fr, frederique.laurent@em2c.ecp.fr,
marc.massot@em2c.ecp.fr
Keywords: Liquid Sprays, multifluid models, numerical analysis, scientific and parallel computing,
EulerEuler/EulerLagrange comparisons.
Abstract
Providing accurate simulations of polydisperse evaporating sprays dynamics in unsteady gaseous flows with large
scale vortical structures is both a crucial issue for industrial applications and a challenge for modeling and scientific
computing. The usual Lagrangian approaches developed in polydisperse unsteady configurations require tremendous
computational costs and may lead to a low level of resolution if not enough numerical parcels are used. Besides,
they induce coupling issues due to the different kind of description of the two phases that are involved. A large
range of Eulerian models have been recently developed to describe the dispersed liquid phase with a lower cost
and an easier coupling with a carrier gaseous phase. Among these models, the multifluid model allows a detailed
description of polydispersity and size/velocity correlations of droplets of various sizes. It has been studied in depth
from a mathematical and numerical point of view (see Laurent & al (2004); Laurent (2006); Massot & al (2009))
and validated through comparisons versus Lagrangian simulations in de Chaisemartin (2009) and experimental
measurements in Fr6ret & al (2008) in 2D and 2Daxisymmetrical configurations. However, the validation in
threedimensional unsteady configurations still remains to be done. In this work, we study the nonevaporating
droplet segregation in threedimensional Homogeneous Isotropic Turbulence (HIT) using a reference Lagrangian
spray model versus the Eulerian multifluid model. A spectral Direct Numerical Simulation solver is used to describe
the evolution of the turbulent carrier phase, whose characteristic properties remain statistically stationary due to a
semideterministic forcing scheme. We focus on the optimization via a parallel implementation of the multifluid
model and dedicated numerical methods which demonstrates the ability of the Eulerian DNS model to be used in high
performance computing for academic threedimensional configurations. We provide qualitative comparisons between
the EulerLagrange and the EulerEuler descriptions for two different values of the Stokes number based on the initial
fluid Kolmogorov time scale, St = 0.17 and 1.05. A very good agreement is found between the mesoscopic Eulerian
and Lagrangian predictions. We go further with first quantitative comparisons of the segregation effect of the vortices
on the spray mass density distribution showing the accuracy and ability of the multifluid model to be used in 3D
configurations from the tracer limit (St 0) to unity.
Introduction character of the droplet size distribution can signifi
cantly affect flame structures. Size distribution effects
In many industrial combustion applications such as are also encountered in a crucial way in solid propel
Diesel engines, fuel is stocked in condensed form and lant rocket boosters, where the cloud of alumina parti
burnt as a dispersed liquid phase carried by a gaseous cles experiences coalescence and becomes polydisperse
flow. Two phase effects as well as the polydisperse in size, thus determining their global dynamical behav
ior, see Doisneau & al (2010). Consequently, it is of
interest to have reliable models and numerical methods
able to describe precisely the twophase flows physics
where the dispersed phase is constituted of a cloud of
particles of various sizes that can evaporate, coalesce
or aggregate, breakup and also have their own inertia
and sizeconditioned dynamics. In a "mesoscopic" de
scription of the liquid phase, droplets are considered as
a cloud of point particles for which the exchanges of
mass, momentum and heat are described using a statis
tical point of view, with eventual correlations, and the
details of the interface behavior, angular momentum of
droplets, detailed internal temperature distribution in
side the droplet, etc., are not predicted. Instead, a finite
set of global properties such as size of spherical droplets,
velocity of the center of mass, temperature are modeled.
Since it is the only one which provides numerical simu
lations at the scale of a combustion chamber or in a free
jet, this mesoscopic point of view will be adopted in the
present paper.
The main physical processes that must be accounted
for are (1) transport in real space, (2) acceleration of
droplets due to drag, (3) droplet heating and evaporation,
and (4) coalescence and breakup of droplets leading to
polydispersivity. Spray models have a common basis at
the mesoscopic level under the form of a number den
sity function (NDF) satisfying a Boltzmann type equa
tion, the socalled Williams equation (Williams (1958)).
The internal variables characterizing one droplet are the
size, the velocity and the temperature, so that the to
tal phase space is usually highdimensional. Such a
transport equation describes the evolution of the NDF
of the spray due to evaporation, to the drag force of the
gaseous phase, to the heating of the droplets by the gas
and finally to the dropletdroplet interactions, such as
coalescence and breakup phenomena. There are sev
eral strategies in order to solve the liquid phase and the
major challenge in numerical simulations is to account
for the strong coupling between all the involved pro
cesses. A first choice is to approximate the NDF by a
sample of discrete numerical parcels of particles of var
ious sizes through a LagrangianMonteCarlo approach
(see O'Rourke (1981)). It is called Direct Simulation
MonteCarlo method (DSMC) by Bird (1994) and gen
erally considered to be the most accurate for solving
Williams equation; it is especially suited for DNS since
it does not introduce any numerical diffusion, the par
ticle trajectories being exactly solved. This approach
has been widely used and has been shown to be effi
cient in numerous cases. Its main drawback is the cou
pling of an Eulerian description for the gaseous phase
to a Lagrangian description of the dispersed phase, thus
encountering difficulties of vectorization/parallelization
and implicitation. Besides, it brings another issue asso
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
ciated with the repartition of the evaporated mass at the
droplet location onto the Eulerian grid for the gas de
scription. Moreover for unsteady computations of poly
disperse sprays, a large number of parcels in each cell
of the computational domain is generally needed, thus
yielding large memory requirement and CPU cost. This
drawback makes attractive the use of a Eulerian formu
lation for the description of the disperse phase, at least
as a complementary tool for Lagrangian solvers.
The Eulerian MultiFluid model, extended by Laurent &
Massot (2001) from the ideas of Greenberg & al (1993),
allows to describe polydispersivity of a spray in size
and the associated sizeconditioned dynamics. This ap
proach relies on the derivation of a semikinetic model
from the Williams equation using a moment method for
velocity conditioned by droplet size while keeping the
continuous size distribution function. This distribution
function is then discretized using a finite volume ap
proach in the size phase space that yields conservation
equations for mass, momentum (and eventually other
properties such as temperature) of droplets in fixed size
intervals. This MultiFluid model is developed in the
framework of a DNS for laminar flows. However, it
is specific of the difficulties one will encounter in the
development of Large Eddy Simulation (LES) tools for
turbulent flows (see Boileau & al (2010)).
In the present work, we consider monodisperse non
evaporating sprays in large scale vortical structures of an
unsteady forced HIT gaseous field. As we do not want
to cope with the difficulties of the twoway coupling but
solely compare two descriptions of the dispersed liquid
phase, we restrict the simulation to oneway coupling
and thus isolate the behavior of each method. The ex
tension to polydisperse evaporating droplets has already
been discussed in Laurent & al (2004); de Chaisemartin
& al (2008); Reveillon & Demoulin (2007), crossing
phenomena in Kah & al (2010) and was not taken into
account here since we tend to validate the multifluid
model in 3D unsteady flows. Simulations were car
ried out thanks to the coupling of a solver dedicated to
the Eulerian spray description (Muses3D) with an other
one which solves the Lagrangian liquid phase and the
gas phase on a Eulerian grid (Asphodele). Muses3D
solver has been developed by S. de Chaisemartin dur
ing his thesis ( de Chaisemartin (2009)) and L. Fr6ret
at EM2C laboratory and Asphodele by J. R6veillon and
coworkers at CORIA laboratory. The coupling of these
two solvers allows the simultaneous numerical simula
tion with the use of two models Eulerian and Lagrangian
for the liquid phase coupled to an unique gas carrier
phase.
The paper is organized as follows. We first recall briefly
the modelling of the liquid phase, both from a La
grangian and Eulerian descriptions, as detailed in de
(I miisc'.iiiiii (2009). We also discuss the numerical
methods for the two different descriptions as well as for
the gaseous phase. Next, we detail the optimization via
MPI parallel implementation done in Muses3D connect
ing with numerical methods and the way we achieved
the coupling between parallel Muses3D and sequential
Asphodele. This allows to carry out simultaneous sim
ulations of two descriptions of a spray dispersion in an
unsteady threedimensional HIT configuration with two
different Stokes, from the tracer limit St 0.17 to
St 1.05 which leads to maximal segregation effects.
Qualitative comparisons are performed considering the
Eulerian liquid density, the Lagrangian droplet positions
and particle number density. Moreover, statistical prop
erties of the dispersed phase such as particlefluid ve
locity correlation and Eulerian droplet segregation are
computed from both Lagrangian and Eulerian simula
tions. Presented results account for our first quantitative
outcome on droplet segregation in an unsteady 3D con
figuration.
Governing equations and modelling
In this study, we restrict ourselves to physical processes
such as (1) transport in real space and (2) acceleration of
droplets due to drag.
In order to introduce the nondimensional equations, we
define the reference velocity Uo and length xo based on
the macroscopic characteristics of the computational do
main allowing to define a reference time scale for the
gaz: T7 xo /Uo. These quantities, along with the phys
ical constants for a reference physical mixture, po Poo
are taken to define the dimensionless system. To derive
the dimensionless equations, we define a normalization
Reynolds number based on the reference quantities
Re0 = pX (1)
The Stokes number is given by St = Tp/r where the
drag relaxation time is defined by Tp piS/(187p,),
the liquid density is pi and p, is the gas viscosity and T,
is the Kolmogorov time scale.
Dimensionless variables are given such as:
u* u/Uo, x* = x/xo, t* = t/9g, S* = S/So.
A number density function (NDF) of the spray f' is
introduced, the quantity f(t, x, S, u)dtdxdSdu be
ing the probable number of droplets with a position
in [x, x + dx], a surface in [S, S + dS], a velocity
in [u, u + du] and at time t. The NDF satisfies a simpli
fied WilliamsBoltzmann equation:
otf + x, (uf) + (Ff) 0,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where F is the dimensionless Stokes law drag force
given by F = (ug u)/St with u, the gas velocity.
Lagrangian liquid phase. To solve the kinetic equa
tion (2) of the spray, we can use Lagrangian Monte
Carlo methods. This leads to EulerLagrange numerical
methods, commonly used for the calculation of poly
disperse sprays in various application fields (generally,
the gas phase is computed using a deterministic Eule
rian solver, while the disperse phase is treated in a La
grangian way). In our study, the Lagrangian reference
is not taken as a converged DSMC method, in order to
be closer to industrial concern. For the given gas DNS
configuration, we perform a Discrete Particle Simulation
(DPS). Indeed there is no need to use Stochastic Parcel
method, since all the droplets contained in the computa
tional domain can be tracked. One can note that in the
infinite Knudsen limit meaning that there is no droplet
interaction, the DSMC computation is equivalent to an
ensemble of DPS, each numerical particle representing
one droplet with a weight equal to one.
The physical processes are then described by the classi
cal following nondimensional equations:
dvk
dt
dxk
dt
1
St(S (x, t) Vk) ,
St(Sk)
where vk and xk denote the dimensionless velocity and
position vectors of each droplet k. The vector ug rep
resents the gas velocity at the droplet position xk. The
right handside term of the first equation in (3) stands for
a drag force applied to the droplet of size Sk.
The Lagrangian numerical particles ODE systems (3)
are solved with an explicit third order Runge Kutta
solver.
Eulerian liquid phase. The alternative to Lagrangian
particle tracking is the resolution of spray Eulerian
global quantities, as number or mass density and mo
mentum. These Eulerian methods can be seen as mo
ment methods derived from the kinetic equation (2).
The formalism and the associated assumptions needed
to derive the Eulerian multifluid model are introduced
in Laurent & Massot (2001). Two steps are to be real
ized in order to obtain the Eulerian multifluid model's
equations. In a first step, the size of the phase space is re
duced by considering zero and one order moments with
respect to the velocity variable at a given time t, position
x and droplet size S : n = fdu and u = ufdu/n
which depend on (t, x, S). The closure of the system is
obtained through the following assumptions:
[H1] For a given droplet size, at a given point (t, x),
there is only one characteristic averaged velocity
u(t, x, S).
[H2] The velocity dispersion around the averaged veloc
ity u(t, x, S) is zero in each direction, whatever the
point (t, x, S) is.
It is equivalent to presume the following NDF condi
tioned by droplet size:
f(t, x, S, u) = n(t, x, S)6(u u(t, x, S)), (4)
that is to reduce the support of the NDF to a one dimen
sional submanifold parametrized by droplet size.
Such an assumption leads to a closed system of con
servation equations called the semikinetic model on
the moments of order zero and one in velocity. It is
given by two partial differential equations in the vari
ables n(t, x, S) and u(t, x, S) which express the con
servation of the number density of droplets and their
momentum, respectively, at a given location x and for
a given size S:
0 9in + 9x (nu) 0,
(5)
t dt(nu) +x (nu iu) nF 0,
where F(t, x, S) is the Stokes's drag force taken at
u u.
The second step consists in choosing a discretization
0 = S(1) < S(2) < ... < S(p) < ... < S(Ns + 1) for
the droplet size phase space and to average the obtained
system of conservation laws over each fixed size inter
vals [Sp, S+1 [, called section. The set of droplets in
one section can be seen as a "fluid" for which conserva
tion equations are written. The sections exchange mass
and momentum. To close the system, the following as
sumptions are introduced:
[H3] In one section, the characteristic averaged velocity
does not depend on the size of the droplets.
[H4] The form of n as a function of S is supposed to
be independent of t and x in a given section, thus
decoupling the evolution of the mass concentration
of droplets in a section from the repartition in terms
of sizes.
The conservation equations for the kth section then read,
in our simplified case:
,1 ,,.' + 9x. (mfu ) = 0,
S (6)
St(,,' i') + dOx (mkuk i uk) = mkFk,
where m is the mass concentration of droplets in the
kth section.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Eulerian multifluid numerical methods
Phenomena involved in our problem are of two differ
ent types: transport induces an evolution in the physical
space without any interaction between the sections, and
by contrast, the transport in internal coordinate space
(velocity) are encountered through drag. This induces
an evolution without any coupling with the spatial co
ordinates. It is then interesting to separate these two
transport types using an operatorsplitting method and
to treat efficiently the different difficulties of the multi
fluid system: a complex transport term (for the physi
cal space) and stiff source terms (for the phase space).
The multifluid system (6) is then split into two systems
that we solve alternatively. We choose a Strang splitting
which is second order in time provided that all the steps
are second order in time, see (Descombes and Massot
2004). The scheme then takes the form:
HAt o GAto GAt o Ht Ht o At o H at, (7)
where G' (respectively Hs) denotes the physical trans
port (respectively the phase space transport) during time
a. This Strang splitting is composed of two Lie split
ting step (H' o G0) of length a At and it alternates
the order in which they are performed (left part of equa
tion (7)). It is equivalent to one Strang splitting step of
length 2At (right part of equation (7)). The splitting
approach has the great advantage to preserve the proper
ties of the schemes we use for the different contributions
such as maximum principle on the velocity or positivity
of density.
In physical space, the system that we get from the
operator splitting is weakly hyperbolic and can gen
erate 6shock and vacuum zones. As precised in de
Chaisemartin (2009), we use second order kinetic
schemes which are finite volume schemes based on the
equivalence between a macroscopic and a microscopic
level of description for the pressureless gas equations
(see Bouchut & al (2003)). These schemes preserve the
positivity of mass density and reproduce a discrete max
imum principle on the velocity. In addition, we use a
dimensional Strang type splitting. This allows to use
schemes in 1D configuration and it preserves the sec
ond order of the method (see LeVeque (2002)). In the
3D case, the scheme then takes the form:
GAt At o o At o GAt o GAt o At
Soy z z oy x
where G' denotes the physical transport in the direction
d (which can be x, y, or z) during time a. This form of
splitting does not lead to CFL reduction in a direction
as a classical Strang splitting algorithm would. Indeed
all the transport substeps are computed for a timestep
At. This splitting allows the same numerical diffusion
in each direction. The global scheme we solve is then
obtained by replacing G At by the expression (8) in (7).
The system obtained in the phase space from the
splitting of system (6) is solved using an efficient ODE
solver based on an implicit RungeKutta 5th order
method RadauIIA (see Hairer & Wanner (1996)).
Eulerian multifluid optimization. As for standard
Eulerian computational fluid dynamics, domain decom
position appears, for multifluid computation, as a very
interesting way to achieve parallel computation. Indeed
it offers the ability to use an arbitrary high number
of points, the number of processes used allowing to
have a subdomain on each process with a reasonable
number of points, leading to reasonable computational
time and memory requirement. The difficulty of such
parallelization lies in the data communications due both
to multifluid peculiarities and dedicated numerical
methods.
The main issue of the multifluid, when dealing with
domain decomposition, is the size discretization lead
ing to an extra dimension of the problem. As this
dimension typically contains five to twenty sections
for "classical" multifluid model, a 3D computational
domain leads to a 4D computation. Furthermore, the
operator splitting used for the numerical scheme makes
two blocks appear with different properties, as far as
domain decomposition is concerned. On the one hand
the physical space transport is local in size and would
naturally lead to a domain decomposition in size, each
process realizing the transport for one size section. In
this case no communication would be necessary in the
physical space. On the other hand, the phase space
transport is local in physical space and would then lead
to a decomposition in space, a process treating a space
subdomain with all the size sections on it. Here, no
communication needs to take place on the phase space.
The different decomposition strategies available in
this context have been evaluated (see de ( I.ic'iii.iiiiii
(2009)). The first strategy is limited by the number
of sections used and will have to be coupled with
a partially spatial decomposition, in order to use an
important number of parallel processes ( 100). An
hybrid method, with size decomposition in the physical
transport, and space decomposition in the size phase
space transport lead to an array reorganization and
consequently to a high amount of communications.
Finally, the domain decomposition in space is a good
compromise between data calculations and MPI com
munications. For this strategy, an efficacity close to one
has been obtained up to 128 cores that was the limit
of our cluster. The extensibility to a thousand core is
eventually foreseen and is the subject of our current
investigations.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Gas solver and turbulence forcing. The gaseous
carrierphase is solved using a solver for Low Mach or
incompressible flows. This solver is based on finite
difference (FD) prediction correction method for the ve
locity evolution, as introduced in Chorin (1968). As far
as numerical methods are concerned, the time resolu
tion is provided by a third order explicit Runge Kutta
scheme. Spatial evolution is done with a FD scheme,
the derivatives being computed with a Pade 6th order
scheme (see Lele (1992)).
At each time step, we perform a turbulence forcing
method to generate an isotropic homogeneous turbu
lence. To maintain the major properties (energy, dis
sipation rate, integral scale) of the spectral turbulence
close to constant values, a controlled amount of en
ergy must be transferred into the spectral simulation
through a forcing procedure. There are various ways to
achieve the forcing of isotropic homogeneous turbulence
in a spectral DNS. In this work, we use the fully con
trolled deterministic forcing scheme (FCDFS) devel
opped by Guichard & al (2004). Inspired from Overholt
& Pope (1998)'s deterministic scheme, FCDFS scheme
has an efficient convergence rate and reduces drastically
the fluctuations of the prescribed properties. Turbulence
is forced by adding a linear source term to the balance
equation for the velocity field a in wavenumber space:
ia =+ u,
at Tf
where a represents the classical NavierStokes contribu
tions for an incompressible flow. The forcing function
fk is real and depends on both time t and wavenum
ber magnitude K. The value Tf is the characteristic re
laxation delay of the simulated spectrum E, towards a
model spectrum Em. The principle of FCDFS model
is to relax the simulated spectrum E, towards a model
one Em only for a given range of low wavenumbers
(K < KF). The interested reader is referred to Guichard
& al (2004) for further information. The vorticity of the
field extracted at dimensionless time t = 20 from this
computation is plotted in Fig 1(a).
3D DNS configuration
A uniformly monodispersed nonevaporating spray
with a zero initial velocity is distributed initially in the
gaseous field in order to study droplet ejection from
the core of the vortices. The drag force sets particles
in motion. Owing to the multifluid size distribution
description, the evaporating case is easy to carry out and
will be addressed in a near future, the nonevaporating
case was a necessary first step.
Droplet dispersion and preferential segregation have
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
25 )4C
3
2 2
x 1
(a) Gas vorticity norm, obtained with
t = 20
x 1 1
spectral resolution at time
N 2
1.5
1
3
(a) Normalized particle number density (/o measured from the La
grangian simulation for St = 0.17 at time t 20
(b) Eulerian spray number density for St
0.17 at time = 20
x 1 1 y
Y
(c) Eulerian spray number density for St = 1.05 at time t = 20
Figure 1: (a) Gas vorticity norm, obtained on a 1293
grid with spectral resolution at time t = 20, slice planes
inside the domain (x = 2, y = 2, z = 1).
Eulerian spray number density obtained with the multi
fluid method on a 1293 Cartesian grid at time t 20.
Slice planes at x = 2, y = 2, z = 1, for two Stokes
number: (a) St = 0.17, (b) St = 1.05.
1 1
Y
(b) Normalized particle number density (/o measured from the La
grangian simulation for St = 1.05 at time t 20
Figure 2: Eulerian spray number density /to obtained
from a Lagrangian simulation of 1283 particles by con
sidering the droplets accumulated around each numeri
cal grid node, at time t = 20. Slice planes at x = 2,
y = 2, z 1 for two Stokes number: (a) St = 0.17, (b)
St 1.05.
been analyzed from a Eulerian point of view. Instanta
neous fields of density are plotted for St = 0.17 in
Fig l(b) and for St = 1.05 in Fig l(c). This latter
case is a case of a Stokes number close to unity leading
to maximal segregation effects. The corresponding
vorticity field is presented in Fig l(a). These three fields
have been captured at exactly the same dimensionless
time t = 20. Even without any quantitative analysis, it
is possible to see the dramatic impact on the particles'
inertia on their dispersion properties. Particles tend to
leave the vortex cores and segregate in weak vorticity
areas. For St = 1.05, the normalized liquid density
varies between 0 (no droplets) and 4 (four times the
mean density).
Numerical setup. We use a cubic computational
domain with periodic boundary conditions. As far as
spatial discretization is concerned, the same grid of
1293 points is used for the gas solver and the Eulerian
spray solver. One has to note that the space domain
decomposition for the multifluid method, implemented
using the MPI library, allows to use for the spray com
putation a refined Eulerian grid at a low computational
time cost. This refinement is already available through
a trilinear interpolation procedure but was not used in
this study since our objectives were the validation of the
Eulerian model in a 3D unsteady and realistic context.
As the gas field is solved sequentially, it is globally
known on each processor and this allows to distribute
uniformly Lagrangian particles on each core. This
reduces significantly the number N1 of Lagrangian
particles on each core keeping the total particles number
N very high. We have N1 N/ni where np is the
core number. As the evaporation phenomena is not
taken into account and there is no interaction between
Lagrangian droplets, numerical particles are indepen
dent from one to each other and no information goes
from droplets to the gas. Hence, there is no need to
consider MPI communications between nodes. No
particles have to be moved from one core to another
as it has to be done when the gas field is divided into
distinct subdomains, see Garcia (2009). An other
point of view is to consider that we are doing np DPS
simulations with few numerical particles. The parallel
efficiency of such method is thus optimal and equal to
one. The simulations corresponding to the results that
are presented were performed on a cluster made of 8
nodes with 2 processors AMD Opteron 64 bits dual core
with speed 2.4 GHz, the 32 cores being connected by
an infinyband gigabit network. A computational time of
twenty hours was necessary to obtain the Eulerian and
the Lagrangian simulation up to a dimensionless time
t 20.
In our simulations, N = 1283 monodispersed
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
nonevaporating particles are randomly embedded
throughout the computational domain with a zero
initial velocity. Hence, thanks to the use of 32 cores,
N1 65530 droplets are tracked through the unsteady
gas field on each core. Droplet dispersion is usually
characterized by the Stokes number St = Tp/Tk, Tk
being the Kolmogorov length scale, which indicates
the ability of droplets to capture local variations of
the carrierphase velocity. Turbulence properties being
fixed, simulations were carried out by varying the Tp
parameter. Two Stokes numbers have been considered,
0.17 and 1.05, corresponding to droplet diameters of
15 and 45 /m. The Reynolds number defined in (1)
is equal to 1000. We eventually provide dimensional
quantities for illustration purposes, based on an esti
mated velocity of Uo 1.5m/s and xo 0.01m, as
well as a typical value of 1.5 x 10 m 2 / s for /po.
In addition, we will let do 0.001 (for St 0.17)
and do = 1i iI.:,, (for St 1.05), where do is the
diameter corresponding to the droplet surface area So.
The computational domain has a size L3 with L = 2,
which corresponds to 8 cm3 in dimensional values.
We take as a reference solution for the liquid phase
the Lagrangian Discrete Particle Simulation with N
particles in the computational domain. We provide
comparisons between this Lagrangian reference and
the Eulerian multifluid monokinetic computations
by plotting the Lagrangian particle positions and the
particle number density measured from the Lagrangian
simulation versus the Eulerian number density.
EulerianLagrangian comparisons
Qualitative liquid dispersion comparisons. The Eu
lerian multifluid description of the spray dynamics are
presented in this section for two Stokes number, based
on the Kolmogorov length scale:
* St = 0.17, corresponding to droplet with diameter
d = 15/m, see Fig l(b).
* St = 1.05, corresponding to droplet with diameter
d = 45/m, see Fig l(c).
These two different inertia allow to study a spray ejected
from the center core and segregated in weak vorticity
areas. They are thus well suited for robustness eval
uation of the multifluid method. Indeed high den
sity regions, as well as vacuum, are created, that rep
resent a challenging issue for a Eulerian method. Higher
Stokes number are not tackled here since it was shown
in Reveillon & Demoulin (2007) that, for Stokes num
ber greater than unity, effects of maximal segregation in
turbulent flow occur. The droplets are inertial enough
to be ejected from a vortex and not follow the fluid
particle like a tracer. Their velocities become decorre
lated from the gaseous carrier phase velocity. In this
case droplet trajectory crossings have a strong impact on
the spray repartition and the monokinetic assumption
of the multifluid method might not allow to describe
it. Recently, Kah & al (2010) extended the multifluid
method to higher order moment method to take into ac
count these droplet crossings and it would be interesting
to enrich the presented comparisons with Eulerian re
sults provided by the extended multifluid multivelocity
model.
To assess the multifluid description of the size
conditioned dynamics, Eulerian density fields are com
pared to Lagrangian droplet positions computed for the
same Stokes numbers at the same time t = 20, see
Fig 2(a) and Fig 2(b). Qualitative comparisons are ap
plied in the planes x = 2 and y = 2. In the chosen iner
tial range, the spray is ejected from the vortex cores and
accumulated in low vorticity areas. In order to link the
spray dispersion given by both methods to the gas vortic
ity structure, the square norm of the gas vorticity is given
in Fig 3(d) and Fig 4(d), for the planes x = 2 and y 2,
respectively. Qualitative comparisons between both ap
proach can be done focusing on the vacuum zones de
scription. These zones correspond to the gas vortex
cores, that can be identified from the vorticity represen
tation. The repartition of these vacuum zones obtained
by the classical Lagrangian method is very precisely re
produced by the multifluid on the different planes (see
Fig 3(a) and Fig 3(e) for the plane x = 2 and Fig 4(a)
and Fig 4(e) for the plane y = 2). Furthermore, the
evolution of droplet repartition with inertia is very well
captured by the multifluid. Indeed, the Eulerian density
fields for higher Stokes number still a present very good
agreement with the Lagrangian droplet repartitions, see
Fig 3(b), Fig 3(c) and Fig 3(f) for the plane = 2 and
Fig 4(b), Fig 4(c) and Fig 4(f) for the plane y =2.
Quantitative spray equilibrium comparisons. As
studied in Reveillon & Demoulin (2007), the equilib
rium of the spray with its carrier phase is detected
through the statistics of the slip velocity
WE(X, ) (U U)(X,t),
where u. and u are the Eulerian gas velocity field and
the liquid velocity field both taken at the numerical grid
point x. From a Lagrangian point of view, the slip ve
locity of a droplet k is given by
wL(k,t)(= (u,(x) V)(t),
where vk is its velocity and u,(xk) is the gas ve
locity taken at the droplet position Xk. The mean
value of the Lagrangian slip velocity is defined by
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
w(t) = 1/N E L(k, t) with N the to
tal particle number while its Eulerian equivalent is
w(t) L3 f wE(x(, t)dx. In the following, an over
bar stands for a Eulerian quantity whereas a tilde denotes
a Lagrangian one. Because of the homogeneous nature
of the turbulence and the dispersion, the mean value
of the slip velocity remains equal to zero. However,
the slip velocity root mean square respectively defined
I ( 2 1/2
by c (t) l/N ( (WL(k,t) ())) and
'(t) L3 (fWE(X t) (t))2 dx) 12 evolves
towards a stationary value co which corresponds to the
equilibrium of the spray with the carrier gas phase. At
initial time t = 0, the liquid spray is initially distributed
in the computational domain with a zero velocity, so
wu (0)/Ti = '(O)/u = 1 with ug the initial
gas velocity r.m.s. As prescribed previously, the drag
force set the liquid spray in motion and, depending on
the Stokes number, the mean value of the slip velocity
reaches a steady state at t = 0.2 for a small Stokes
number (St 0.17) whereas a longer time is needed
for the largest value of the Stokes number (St 1.05).
The final mean stationary value of the slip velocity stan
dard deviation wc is close to zero when the droplets are
small enough to follow all the velocity fluctuations of
the flow (St = 0.17) and it increases to reach 0.21 cor
responding to droplets of St 1.05. This result accords
with those presented in Reveillon & Demoulin (2007)
and are plotted in Fig 5. Values of wZ and W'o show
a very good agreement for St 0.17, whereas wZU is
smaller than c', and equal to 0.16 for St 1.05.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
25 3
15 2
Y
(a) Eulerian multifluid number density
for low inertia droplets St 0.17
(b) Lagrangian droplet positions
for low inertia droplets St 0.17
(c) Normalized particle number density
(/o measured from the Lagrangian
simulation for low inertia droplets St
(d) Instantaneous Vorticity field
16~
15 2 25 3
(e) Eulerian multifluid number density
for higher inertia droplets St 1.05
(f) Lagrangian droplet positions
for higher inertia droplets St 1.05
(g) Normalized particle number density
(/o measured from the Lagrangian
simulation for low inertia droplets St =
Figure 3: Eulerian Lagrangian comparisons of the liquid phase for two different Stokes St = 0.17 and St = 1.05, in
the (y z) plane at x = 2 at time t = 20.
15 2
Y
15 2
Y
0.17
15 2
Y
 i"
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
2
25 15
05
15 2 25 3 1 15 2 25
Y x
(a) Eulerian multifluid number density
for low inertia droplets St 0.17
(b) Lagrangian droplet positions
for low inertia droplets St 0.17
(c) Normalized particle number density
(/(o measured from the Lagrangian
simulation for low inertia droplets St
(d) Instantaneous Vorticity field
1 15
25 3
(e) Eulerian multifluid number density
for higher inertia droplets St 1.05
(f) Lagrangian droplet positions
for higher inertia droplets St 1.05
(g) Normalized particle number density
(/(o measured from the Lagrangian
simulation for higher inertia droplets St
Figure 4: Eulerian Lagrangian comparisons of the liquid phase for two different Stokes St = 0.17 and St = 1.05, in
the (x: z) plane at y = 2 at time t = 20.
0.17
Eulerian rms. St=0.17
Lagrangian rms. St=i 17
Euleran rlls. SL=1.05
Lagranmoan rms. St=l )5
0 02 0 04 U 06 018
time (en ms)
0.1 0 12
Figure 5: Time evolution of the Eulerian (w (u)/u )
and Lagrangian (_' (t)/ug) slip velocity r.m.s. for
St 0.17 and St 1.05.
Conclusions and perspectives
We have presented comparisons between Lagrangian
and Eulerian liquid spray computations in a context
of threedimensional unsteady forced HIT carrierphase
gas. Results were presented for two different Stokes
values showing the impact of inertia on droplet seg
regation as well as the robustness of the multifluid
model since Stokes values from tracer limit to unity were
used. Qualitative comparisons were conducted consid
ering two slice planes and a very good agreement has
been obtained up a Eulerian density prediction quanti
tatively close to DPS results. Preferential concentration
was quantitatively studied through the comparison of the
slip velocity r.m.s. These results are a first step in spray
equilibrium comparisons between the Eulerian and La
grangian approach and will be completed in a near fu
ture in both improved quantitative comparisons and an
extension to evaporating case.
Acknowledgements
This research was supported by an ANR05JCJC0013
Young Investigator Award (PI : M. Massot, 20052009)
which provided a postdoctoral grant for L. Freret
at EM2C Laboratory (20072008) and a DGA/CNRS
Ph.D. grant for S. de ( Ii.iinc'sii.llii (Mathematics and
Engineering Departments of CNRS 20052009). The
authors also acknowledge the support from IDRIS
CNRS (Institut de Developpement et de Ressources
en Informatique Scientifique, Centre National de la
Recherche Scientifique) where some of the computa
tions were performed.
S05V
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
References
Bird, G.A., Molecular gas dynamics and the direct simu
lation of gas flows. Oxford science publications, vol 42,
1994.
Boileau, M., Chalons, C., Bourgouin, J.F, Terrier, C.,
Laurent, F., de Chaisemartin, S., Massot, M., Robust
numerical schemes for Eulerian spray DNS and LES
in twophase turbulence flows, 7th International Confer
ence on Multiphase Flow, 2010.
Bouchut, F., Jin, S., Li, X., Numerical approximations
of pressureless and isothermal gas dynamics. SIAM J.
Num. Anal, vol 41, pp 135158, 2003.
Chorin, A.J., Numerical solution of the NavierStokes
equations. Math. comp. 22, vol 217, pp 745762, 1968.
de Chaisemartin, S., Polydisperse evaporating spray tur
bulent dispersion: Eulerian model and numerical sim
ulation. PhD thesis, Ecole Centrale Paris, available
in english at http://tel.archivesouvertes.fr/tel00443982/en/,
2009.
de Chaisemartin, S. Fr6ret, L., Kah, D., Laurent, F, Fox,
R.O., Reveillon, J., Massot, M., Eulerian models for tur
bulent spray combustion with polydispersity and droplet
crossing. C. R. Mecanique vol 337(67), pp 438448,
2009.
Descombes, S. Massot, M. Operator splitting for non
linear reactiondiffusion systems with an entropic struc
ture: singular perturbation and order reduction. Numer.
Math. 97(4), pp 667698, 2004.
Doinseau, F, Laurent, F, Murrone, A., Dupays, J., Mas
sot, M., Optimal Eulerian model for the simulation of
dynamics and coalescence of alumina particles in solid
propellant combustion, 7th International Conference on
Multiphase Flow, 2010.
Fr6ret, L., Lacour, C., de Chaisemartin, S., Laurent, F.,
Massot, M., Birbaud, A.L., Ducruix, S., Durox, D., Pul
sated free jets with polydisperse spray injection: exper
iments and numerical simulations. Proceedings of the
combustion institute, vol 32(2), pp 22152229, 2009.
Garcia, M., Development and validation of the Euler
Lagrange formulation on a parallel and unstructured
solver for largeeddy simulation. PhD thesis, Institut Na
tional Polytechnique de Toulouse, available in english at
http://tel.archivesouvertes.fr/tel00414067/en/, 2009.
Greenberg, J.B., Silverman, I. Tambour, Y., On the ori
gin of spray sectional conservation equations. Combus
tion and Flame, 93: 9096, 1993.
Guichard, L. R6veillon, J., Hauguel, R., Direct numer
ical simulation of statistically stationary one and two
phase turbulent combustion: a turbulent injection pro
cedure, Flow Turbulence Combustion, Vol. 73, pp 133
167, 2004.
Hairer, E., Wanner, G., Solving ordinary differen
tial equations. II. Berlin: SpringerVerlag. Stiff and
differentialalgebraic problems, second revised edition,
1996.
Kah, D., Laurent, F., Fr6ret, L., de Chaisemartin,
S., Fox, R.O., R6veillon, J., Massot, M., Eule
rian quadraturebased moment models for polydisperse
evaporating sprays, accepted in Flow, Turbulence and
Combustion, available at http://hal.archivesouvertes.fr/hal
00449866/en/, 2010.
Kaufmann, A., Moreau, M., Simonin, O., Helie, J.,
Comparison between Lagrangian and mesoscopic Eu
lerian modelling approahces for inertial particles sus
pended in decaying isotropic turbulence, Journal of
Comp. Phys. vol 227, pp:66486472, 2008.
Laurent, F. and Massot, M., Multifluid modeling of
laminar polydispersed spray flames: origin, assump
tions and comparison of the sectional and sampling
methods, Combustion theory and Modelling, vol 5, is
sue 4, pp 537572, 2001.
Laurent, F., Massot, M. Villedieu, P., Eulerian multi
fluid modeling for the numerical simulation of coales
cence in polydisperse dense liquid sprays. Journal of
computational physics, vol 194, pp 505543, 2004.
Laurent F., Numerical analysis of Eulerian mutlifluid
models in the context of kinetic formulations for dilute
evaporating sprays, M2AN Math. Model. Numer. Anal.
40(3): 431468, 2006.
Lele, S.K., Compact finite difference schemes with spec
tral like resolution. J. Comput. Phys. 103, pp 1642,
1992.
LeVeque, R.J., Finite volume methods for hyperno
lic problems. Cambridge texts in applied Mathematics.
Cambridge University press, Cambridge, 2002.
Massot, M., de Chaisemartin, S., Fr6ret, L., Kah, D.,
Laurent, F., Eulerian multifluid models: modeling
and numerical methods. In Modelling and Computation
of Nanoparticles in Fluid Flows, Lectures of the von
Karman Institute. NATO RTO AVT 169, available at
http://hal.archivesouvertes.fr/hal00423031/en/, 2009.
O'Rourke, P.J., Collective drop effects on vaporizing liq
uid sprays. PhD thesis, Princeton University, 1981.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Overholt, M.R. and Pope, S.B., A deterministic forcing
scheme for direct numerical simulation of turbulence.
Computer fluids, vol 27, pp 1128, 1998.
Poinsot, T. and Lele, S.K., Boundary conditions for di
rect simulations of compressible viscous flows. J. Com
put. Phys. 101(1), pp 104129, 1992.
Reveillon, J., P6ra, C., Massot, M., Knikker, R., Eulerian
analysis of the dispersion of evaporating polydispersed
sprays in a statistically stationary turbulent flow. Journal
of turbulence, 5(1), pp 127, 2004.
Reveillon, J. and Demoulin, EX., Effects of the prefer
ential segregation of droplets on evaporation and turbu
lent mixing. Journal of fluids mechanics, 583, pp 273
302, 2007.
Reveillon, J. and Demoulin, F.X., Evaporating droplets
in turbulent reacting flows, Proceedings of the combus
tion institute, vol 31, issue 2, pp 23192326, 2007.
Simonin, O., Fevrier, P., Lavi6ville, J., on the spatial dis
tribution of heavyparticle velocities in turbulent flow:
from continuous field to particulate chaos. J Turbulence,
vol 118,pp 97118, 1993.
Squires, K.D. and Eaton, K.J., Presential concentration
of particles by turbulence. Phys. Fluides, vol 3, pp 1969
1178, 1991.
Williams, F.A., Spray combustion and atomization.
Physics of fluids, vol 1, pp 541545, 1958.
