7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Robust numerical schemes for Eulerian spray DNS and LES in twophase
turbulent flows
M. Boileau*, C. Chalonsi*, J.F. Bourgouin*, C. Terrier*, F. Laurent*,
S. de Chaisemartin~ and M. Massot*
Laboratoire EM2CUPR CNRS 288, Ecole Centrale Paris, 92295 ChiltenayMalabry, France
t DEN/DANS/DM2S/SFME/LETR CEASaclay, 91191 GifsurYvette, France
i:Institut Frangais du petrole, 92852 Rueil Malmaison, France
matthieu.boileau ~em2c.eep.fr, chalons ~math~jussieu.fr, frederique.1aurent ~em2c.eep.fr,
stephane.dechaisemartin~ifp.fr and mare.massot~em2c.eep.fr (corresponding author)
Keywords: Liquid Sprays; MultiFluid models; Spray Equation; Large Eddy Simulation; Relaxation schemes
Abstract
Large Eddy Simulation (LES) and Direct numerical Simulation (DNS) of polydisperse evaporating sprays with
Eulerian models are very promising simulation tools for combustion applications for high performance computing
in order to describe turbulent dispersion and evaporation and properly predict combustion regimes. However, the
resulting systems of conservation equations have a convective part which is either similar to gas dynamics Euler
equations with a real gas type state law or to the pressureless gas dynamics (PGD), depending on the local flow regime
and droplet Stokes number; so, they usually involve singularities due to model closure assumptions and require
dedicated numerical schemes. Besides, it is desirable to cope with exactly zero droplet density in some zones of the
flow, especially near the injection zone, where droplets are injected in only some spatial locations. Even if the issue
has been successfully tackled in de Chaisemartin (2009); Frdret et al. (2010) in the framework of PGD with the use of
kinetic schemes with very low diffusion levels, it can not be directly extended to general gas dynamics. The purpose
of the present contribution is to introduce a new generation of numerical methods based on relaxation schemes which
are able to treat both PGD and general gas dynamics, as well as to cope in a robust manner with vacuum zones and
natural singularities of the resulting system of conservation equations. The proposed hybrid relaxation scheme and
algorithms are validated through comparisons with analytical solutions and other numerical strategies on 1D and 2D
configurations. They exhibit a very robust behavior and are a very promising candidate for more complex applications
since they provide solutions to key numerical issues of the actual Eulerian spray DNS and LES models.
Introduction
Many industrial devices involve turbulent combustion of
a liquid fuel. The transportation sector, rocket, aircraft
or car engines are almost exclusively based on storage
and injection of a liquid phase, which is sprayed into
a combustion chamber. It is of primary importance to
understand and control the physical process as a whole,
from the injection into the chamber up to the combus
tion phenomena. Numerical simulation is now a stan
dard industrial tool to optimize the turbulent combus
tion process in such devices (Duchaine et al. cano ' ,
Thanks to LES, unsteady phenomena such as jet igmi
tion (Lacaze et al. (2009)) or combustion instabilities
(Selle et al. (2006); Roux et al. (2008)) can now be
accurately predicted in simplified configurations where
purely gaseous flames are encountered. Nevertheless,
the liquid fuel injection needs special attention in or
der to properly predict the combustion regimes and con
sists in two parts. The first one is related to the at
omization process near the injector and requires dedi
cated models and methods. On the other hand, the sec
ond part is related to the spray dynamics once the liquid
has reached the structure of a cloud of of polydisperse
droplets; some promising advances have been performed
in the field of spray combustion in real devices (Boileau
et al. (2008a,b); Vid et al. (2010)). However, the re
liable prediction of such complex twophase reacting
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
in order to capture strongly non equilibrium velocity
distributions, has been developed in order to capture
droplet crossing trajectories and will not be treated in
the present contribution. The second one, adopted in
the present paper, is based on equilibrium velocity dis
tributions either with zero dispersion around the mean
in the framework of DNS such as in de Chaisemartin
(2009); Freret et al. (2010) and references therein or
with nonzero dispersion in the framework of ensemble
averages and modeling nonresolved scales such as in
de Chaisemartin (2009) and Vi6 et al. (2010). In the
context of spray dynamics, a zero pressure assumption
means that the probability density function of particle
velocity is a Dirac in the velocity space, i.e. that no
dispersion in the local instantaneous particle velocity is
considered. For this assumption to be true, the relax
ation time of particles must not exceed the timescale
of the fluid turbulence. Otherwise, the effect of trajec
tory crossing due to higherinertia particles induces a
randomuncorrelated component of particle motion that
must be taken into account. Simonin et al. (2002) and
Kaufmann et al. (2008) proposed a formalism that ac
count for this randomuncorrelated motion in the con
text of DNS. Moreau et al. (2005) and Riber et al. (2005,
2009) have extended this approach for LES. They used a
spatial filtering of Kaufmann's system of Eulerian con
servation equations and proposed a model for the result
ing particulate subgrid stresses. Another way to use the
LES concept for the spray equations is to apply the LES
filtering directly to the Williams equation for the PDF of
particle velocity (see Pandya and Mashayek (2002) and
Zaichik et al. (2009)). The resulting kinetic equation for
the filtered PDF has the same form as the statistical PDF
equation initially derived by Reeks (1992) in the con
text of RANS and used in Massot et al. (2004). What
ever the approach for turbulence modelling (DNS, LES
or RANS) and the level of corresponding filtering (on
the kinetic equation or on the moment equations at the
semikinetic level), the local velocity dispersion of par
ticles introduces stress, and more specifically a pressure
like term, in the spray conservation equations. Finally,
the Eulerian equations for inertial particles dynamics
are similar to the gas dynamics equations featuring, in
particular, they include a real gas type state law which
can eventually degenerate in some parts of the flow to
a zero pressure term leading to the peculiar PGD. Let
us emphasize that the size distribution function is then
discretized using a finite volume approach in the size
phase space that yields conservation equations for mass,
momentum (and eventually other properties such as en
thalpy) of droplets in fixed size intervals which have the
same mathematical structure. In the present paper, we
will consider a monodisperse spray so that the semi
kinetic model is sufficient, keeping in mind that all the
flows requires further work in the modeling of the triple
spray/turbulence/combustion interaction. In particular,
the description of the turbulent spray dispersion remains
a challenging issue. 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 total
phase space is usually highdimensional. Such a trans
port equation describes the evolution of the NDF of the
spray due to convection, heating and evaporation, drag
force from the gaseous phase and dropletdroplet inter
actions. Different strategies can be used to solve the dis
persed phase dynamics. A first choice is to approximate
the NDF by a sample of discrete numerical parcels of
particles 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 del
icate coupling between the Lagrangian description of
the dispersed phase and the Eulerian description for the
gaseous phase. From a computational point of view, a
Lagrangian solver is difficult to efficiently parallelize us
ing the domain decomposition used by the gas solver.
This is particularly true in massively parallel calcula
tions where only a few parallel blocks may contain most
of the Lagrangian particles and requires dedicated algo
rithms as in Garcia (2009). Finally, unsteady computa
tions of polydisperse sprays requires a large number of
parcels in each Eulerian cell, leading to large memory
needs and high CPU costs. As a consequence, an Eu
lerian formulation for the description of the dispersed
phase, as long as it is able to describe the essential fea
ture of polydispersity, is more attractive for massively
parallel simulations of industrial configurations.
Based on the ideas of Greenberg et al. (1993) and
following the work of Laurent and Massot (2001),
de Chaisemartin (2009) have developed a multi
dimensional Eulerian MultiFluid solver capable of de
scribing the polydispersity of a spray in size and the as
sociated sizeconditioned dynamics. This approach re
lies on the derivation of a semikinetic model from the
Williams equation using a moment method for veloc
ity conditioned by droplet size while keeping the con
tinuous size distribution function. The key issue then
relies in this velocity moment closure. Two strategies
can be chosen. The first one is based on quadrature
method (see Kah et al. (2010) and references therein)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
port step, as well as a strong relaxation step related to
a singular perturbation parameter. A large literature on
the subject has shown the impact of such seminal ideas
(Chalons and Coquel (2005), Bouchut (2004), Chalons
and Coulombel (2008)).
In this contribution, we conduct three new steps : 1)
we propose a scheme for PGD based on successive en
ergy and pressure relaxation which can deal with vac
uum and extends the work of Berthon et al. (2006), 2)
based on this new scheme, we introduce a hybrid nu
merical method which is able to deal with vacuum and
can treat both regions with and without pressure and still
remain accurate and robust, 3) we finally prove the po
tential of the proposed schemes by comparing them on
several tough testcases to standard approaches in both
1D and 2D configurations. Since relaxation methods
are able to treat arbitrary state law, we only provide the
schemes in the framework of ideal gas law; besides we
focus on the purely convective part of the system of con
servation laws and do not treat the potential stress ten
sors which can be dealt with using standard schemes.
The anticle is organized as follows. The first section de
scribes how the Eulerian description of turbulent spray
dynamics leads to a gas dynamicstype system of con
servation equations. The hybrid relaxation scheme is de
rived in the second section. Finally, results for relevant
testcases, in 1D/2D and pres sureles s/pres sure/hybrid
configurations, are presented and discussed.
1 Eulerian modelling of turbulent spray
dynam ics
Conservation equations on moments of the particle
number density function
At the inesoscopic level, spray models have a common
basis called the kinetic model by analogy with kinetic
theory of gases. The spray is described as a statisti
cal cloud of point particles experimenting exchanges of
mass, momentum and heat with the carrier phase. This
kinetic model is described by a Boltzmann type equa
tion (equation 1) for the number density function (NDF)
f of the spray, where f(t, x, u)dxdu denotes the proba
ble number of particles at time t, in a volume of size dx
around x, with a velocity in a duneighborhood of u. As
mentioned in the introduction, other physical properties
like the particle size and temperature can be introduced
in the NDF for a finer description of the spray in the
framework of the multifluid model introduced in Lau
rent and Massot (2001) and for which references are to
be found in de ( Ilnlc'sllliallion (2009). For sake of sim
plicity, constant particle size (monodisperse spray) and
temperature are considered here so they will not appear
in the equations.
developed tools can be easily extended to polydispersity
in the framework of the multifluid method de Chaise
martin (2009); Vid et al. (2010).
The main difficulty of the resulting system of conser
vation equations is related to transport in physical space,
that is the convective part of the system, which is ei
ther hyperbolic of weakly hyperbolic, and thus leads to
singularity formation. In the framework of the PGD
system, de ( Ilnlc'sllliallion (2009) has solved the prob
lem by using the second order kinetic scheme proposed
by Bouchut et al. (2003) and by proposing a numerical
strategy which leads to second order in space and time
schemes with very little diffusionde ( Insivllusil~lllli et al.
(2009). This numerical scheme makes it possible to cap
ture accurately the deltashocks in density and vacuum
states which naturally emerge from the weakly hyper
bolic system. However, this strategy can not be extended
in a natural manner to gas dynamics either with perfect
gas law or with real gas state laws.
Thus, the numerical method we are looking for must
have the ability 1) to handle a Eulertype system of
equations in regions of high Stokes number or in re
gions where subgrid scales induce \ignlilk~.lll pressure
effects, 2) to degenerate to the PGD system in regions
of Stokes number below the critical value for particular
trajectory crossing or in regions where the subgrid scales
do not play any role in particle velocity dispersion, 3)
to treat exact vacuum regions for both pressureless and
full gas dynamics systems in various regions of the flow.
Moreover, this method must feature the same properties
of robustness with singularities and vacuum treatment
as the Bouchut's kinetic scheme for PGD proposed in
de ( Ilnlc'sllliallion (2009); Massot et al. (2009) where the
level of numerical diffusion is kept very low even if the
scheme is second order in space and time. Finally, since
the pressure law can bear some real gas effects, the nu
merical method has to handle such cases while keeping
a high level of accuracy as required by the DNS/LES
approach.
In that context, the purpose of the present paper is to
introduce a new numerical method based on relaxation
schemes which has the ability to match all the previous
requirements. Relaxation methods, introduced in Jin and
Xin (1995), and further developed in Suliciu (1998) and
Coquel and Perthame (1998), have a common basis and
introduce auxiliary variables in the framework of Go
dunov schemes in order to treat more easily the strong
nonlinearity due to the treatment of pressure and state
law. They allow to avoid the resolution of complex non
linear Riemann solvers or of their approximated versions
which can have a very high computational cost with non
standard pressure laws. The nonlinearity treatment is
replaced by a splitting like strategy in the framework
of a linear or linearly degenerate version of the trans
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The evolution of the spray NDF is given by the
Williams transport equation Williams (1958):
8, f + u 8xf + a, (F f) = 0, (1)
where F is the drag force due to the velocity difference
with the gaseous phase and given by the Stokes law:
U(t, x) u P d2
F(t, x, u) = 7 with 7>> = 18, (2)
where U is the gas velocity at the particle location, p,
its dynamic viscosity, pi is the liquid density and d is the
droplet diameter.
The first possibility is to write conservation equations
for the zero and first order moments with respect to the
velocity variable at a given time t and position x:
n(t, x) = f (tx)du, (3)
n(t,x) U(t, x) = ftxd (4)
where n is the particle density and u is the particle mean
velocity. At this stage, there are two different ways
of deriving the conservation equations for these two
moments according to the value of the particle Stokes
number St, defined by: St viap/7K, where 7K is the
Kolmogorov time microscale.
Pressureless gas system
For Stokes number low enough, particles have a low iner
tia and do not experiment any trajectory crossings. Ac
cordingly, the velocity dispersion around the averaged
velocity u(t, x, S) is assumed to be zero in each direc
tion the spray is called monokinetic and the NDF
writes:
f (t, x, u) = n(t, x)6(u u(t, x)). (5)
Such an assumption leads to a closed system of conser
vation equations given by two partial differential equa
tions in the variables n(t, x) and u(t, x) which express
the conservation of the number density of droplets and
their momentum respectively:
Gas dynamics system
As pointed out in the introduction, the monokinetic as
sumption is not verified for larger Stokes number, i.e. for
particle relaxation times greater than the Kolmogorov
time scale, where the effects of particles trajectory cross
ings yields the need for additional higher order moment
modeling. In particular, these crossings are expected
to reduce the particle segregation induced by inertia ef
fects. A way to account for the uncorrelated motion of
inertial particles is to the used the mesoscopic formalism
proposed by Fevrier et al. (2005), starting from the fol
lowing decomposition: u u(t, x) + 6u, where 6u is
Called the random uncorrelated component of the parti
cle velocity. The equations in system (6), obtained for a
monokinetic spray, now becomes (see Kaufmann et al.
(2008)):
8,n + 8 (nu) = 0
8t(nu) + 8 (nu <.u +P) = nF + 867
(7)
dt(nS)+8x(nlu+Pu) = nFu
2 n69 + add. terms
where the total energy reads I uu/2+69, with 68 the
random uncorrelated energy defined as half the trace of
the random uncorrelated stress tensor and where P is
called the random uncorrelated pressure which is linked
to the random uncorrelated energy through the following
equation of state:
p = n bO (8)
In system (7), by is the deviatoric part of the random
uncorrelated motion tensor and it can be modeled by a
viscosity assumption. These equations correspond to the
case where the gas flow is entirely resolved and no mod
eling of the gas turbulence is used (DNS approach). In
the context of statistical (RANS) filtering (Reeks (1992);
Zaichik (1999); Massot et al. (2004); de Chaisemartin
(2009)) or LES filtering (Moreau et al. (2005); Riber
et al. (2005); Pandya and Mashayek (2002); Zaichik
et al. (2009)), the pressure law becomes more com
plicated than the Eq. (8), involving contributions from
turbulent or subgrid motion respectively. The modeled
scales involve real gas effects through a modification of
the state law, as well as source terms in the rhs of sys
tem (7) for the random uncorrelated energy.
The simplified general form of the system of conser
vation equation finally considered in the following is
then system (7) with by 0 without the additional
terms in the energy equation but with a potential source
term, that is 2so69 is replaced by 2n (68 ~t), where it is
the energy source term due to subgrid turbulence agita
tion. From a numerical point of view, we thus isolate the
difficulties to solve system (7), whicharise mainly from
8 ,(n u) + 812x (nd (u~u u)
where the Stokes drag F is taken at u u. Equation (6)
is similar to the PGD system with an additional velocity
relaxation source term.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
At least formally, we observe from the last equation in
(10) that the relaxation internal energy pe tends to zero
as the relaxation parameter A > 0 goes to infinity. By
(11), the solutions of the relaxation system (10) are thus
expected to provide a good approximation of the solu
tions of the PGD for large values of A. Note that if we
define the temperature T and the mathematical entropy
S according to the second principle of thermodynamics
TdS = de pd7, 7 = 1/p,
easy calculations lead to the expected entropy inequality
the pressure terms. They requires numerical method for
highly compressible flows. The additional source terms
and second order derivatives usually do not lead to nu
merical difficulties and can be treated through operator
splitting, whereas the main difficulties arise from the
convective first order part involving the pressure effects.
Therefore, in the next subsection, we will focus on the
the lhs and forget temporarily the rhs in order to build
the numerical schemes. The strategy adopted is to fo
cus on the perfect gas state law in the following, but us
ing relaxation methods, which can be easily extended to
any state law with real gas effects usually encountered
in LES. Finally, we also need to be able to treat cases
where the random uncorrelated energy can be zero and
the previous system degenerates toward the PGD.
2 A hybrid relaxation scheme for
gas/pressureless gas dynamics problems
Our objective in this section is to describe a global nu
merical strategy in 1D, able to deal with both gas dynam
ics and PGD at the same time, and to handle vacuum. It
is based on the concept of relaxation approximation for
systems of conservation laws. The basic idea is to pro
pose an enlarged system with a stiff relaxation source
term, the solutions of which are expected to converge
to the solutions of the initial system in the asymptotic
limit. In the following, the notations of the previous sec
tion are abandoned and replaced by the usual notations
for hyperbolic systems of conservation laws.
Towards a unified treatment of gas and pressureless
gas dynamics
We first propose to write the pressureless gas system (6)
under the equivalent form
dt 8,+ 8, (per) =0,
8, (per) + 8, (pU2 + p) =0, (9
with p = 0. Then, following the general idea of Co
quel and Perthame (1998), we propose to approximate
the solutions of this system by the ones of the energy
relaxation system
dt pS)+ dpg
Xp~< O
It is important to note that the zero internal energy equi
librium manifold, that is also the zero pressure manifold
in the limit of infinite A, is stable in the sense that for
initial data with zero pressure, the dynamics remain with
zero pressure.
The numerical procedure we are going to propose in
order to approximate the solutions of the pressureless
gas dynamics system (6) is very classical in the context
of relaxation approximations. It is based on an operator
splitting for (10) and is made of two steps that we now
briefly describe.
First step: We solve the convective part of the model:
dtP d pu)=0
dt (p)+d p"2 p) =0
8, (pE) 8 (pEtt + per) 0,
which is nothing but the classical gas dynamics system.
We will use the condensed form
8,U + 8, F(U) 0
for (12) with clear definitions for U and F(U).
Second step: In the second step, the contribution of the
stiff relaxation source term is accounted for by solving
the ODE system
i, p = O
dt(pr) = 0
dt (pE) =
i, p = O
d t(pr)= O
87 = 
8, (pt) + (pU2 + p) = 0,
8, (pE) 8 (pEtt + per)
in the asymptotic regime A 00. This clearly amounts
to keep p and pts unchanged and to set e 0 that is
pE ~pu2 and p 0.
A pressure relaxation model for gas dynamics
In this paragraph, we propose a pressure relaxation sys
tem in order to approximate the solutions of the gas dy
namics system (12). Motivated by the seminal work of
Jin and Xin (1995) and by Suliciu (1998), we propose to
relax the nonlinearities associated with the pressure law
Ape,
where the socalled relaxation internal and total energies
are related by1
[InsposuIhI.I1 the pressure p here no longer equals zero
but obeys for instance a perfect gas equation of state
~~ (11)
P 31
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
p only, and to retain the other ones for the sake of accu
racy. With this in mind, we introduce the following non
linear first order system with singular perturbation:
dep + d, (pu) = 0,
a, (pu) + 8, (pU2 + I) (15)
a, (pE) + 8, (pEu + Inu) 0 ,
a, (pil) + Ox (plnu + aku) = pp(p II),
that we write for shortness
amv + a,g(v) = ex(v/).
As p goes to infinity, we observe at least formally that
the relaxation pressure II tends to p so that the equilib
rium system (12) is recovered in this asymptotic regime.
The additional equation associated with II is easily seen
to be equivalent to
8~ti + u 8il, +  i p I).
This equation is then very similar to the one associated
with the exact pressure p given by
8tp + u xp + pcxu = 0.
The choice of the parameter a > 0 is crucial for the
stability of the relaxation procedure and is determined
by the socalled subcharacteristic condition a > pc
where c denotes the sound speed.
The firstorder system extracted from (15) is hyperbolic
and admits the following three eigenvalues,
a a
with secondorder multiplicity for Ag. We note that At
and As approximate the characteristic speeds n c and
u + c of (12). Insposuih.I11 these eigenvalues are now
associated with linearly degenerate characteristic fields.
This implies that the Riemann problem associated with
(15) (with p 0) can be explicitly solved, unlike the
one associated with (12). Riemann solutions being the
key ingredient to devise Godunovtype methods, this
mathematical property justifies the introduction of the
relaxation model (15).
Here again, the proposed numerical procedure to
approximate the solutions of the gas dynamics system
(12) is based on an operator splitting for (15) and is
made of two steps:
First step: We solve the convective part of the pressure
relaxation model taking p 0 in (15):
dtp + d, (pu) = 0,
at (pu) + 8, (pU2 + I) ,,
8t (pE) + 8,(pEu + Inu) 0,
8t (pil) + Ox (plnu + aku) = 0,
or equivalently
at? + a,0(v) = 0.
In practice, we will use in this step a Godunov method
based on the exact Riemann solution of (16).
Second step: We then solve
iSep = 0,
at (pu) = 0,
8,(pE) = 0,
in the asymptotic regime p 00o. The conservative
variables p, pu and pE are thus constant, while II is set
to be equal to p in this step.
For the sake of completeness, we now give the
Riemann solution associated with (16). We propose
to take a nonconstant in the Riemann solution and we
choose to solve
dta + uda =0.
The diagonal form of (16)(17) is given by
Be (1 + an) + (u + ar)8, (11 + an)~a+u 3a=0
at(II au) + (u ar)8x(II au)
at (I + a 7) + u 8, (II + a2) = 0,
8ta+ u a = ~ 0. 28
In other words, the quantities (II f an), respectively
(II+ a7), (E ~i) and a, are (strong) Riemann invari
ants for the eigenvalues vi ar, resp. u. The calculations
are left to the reader.
Let be given jL (ZAL, (pll)L) and 1a (Zda, (pll)n)
two constant states and let at and aR be two values
for a. The selfsimilar Riemann solution (x, t) H
l/(x/t; IjL, VR; at, aR) associated with (16) and initial
data i <0
,,.x, t 0 VRV ifi x > 0,:
is made of four constant states VL, VI,~ VA and VR, sep
arated by three contact discontinuities associated with
Xk =X(1/), k 1 2, 3and propagating with speeds
denoted by X(1J, VIj), X(VLI, VA) and X(1/, VR). More
precisely, we have
yJx vt if x(vL, vt)< "
t v~~P if X(VA, Va) < X (,V)
The intermediate states 1I, 1/A, as well as the speeds
of propagation, are determined using for all k 1 2, 3
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
of propagation X1(1/) and Xs(VR) remain finite. In par
ticular, discrete entropy inequalities as well as maximum
principles can be proved. These results are pretty tech
nical to establish and are not presented in this paper. We
refer the reader to Bouchut (2004) for the details.
In the case of PGD, these formulas are to be considered
with pt pR 0 and ct cR 0 We then observe
that the threshold cmin allows to guarantee (18) when
UL < UR and then to avoid the resonance phenomenon.
A relaxation scheme for the gas and pressureless
gas dynamics
In this paragraph, we present a relaxation scheme for ap
proximating the solutions of the gas or pressureless gas
dynamics equations (13) and (10) separately. The case
of mixed computations involving both the gas and pres
sureless gas dynamics at the same time will be consid
ered in the next paragraph. It is important to notice that
the same formalism will be used for both systems. Just
note that in the pressureless case, E must be understood
as a function of the unknowns p and pu, namely
(p)2
E=2p
but not as an unknown with evolution given by the pas
sive transport equation
the continuity of the (strong) Riemann invariants for
As across the contact discontinuity associated with At,
I / k. We get after easy calculations A(Vl, 7~)
As (1J) = UL aL7L, Aff L, VA) = U*, A(VA, VR) =
AS(1VR) = UR + agr7 and
n, up= u =areL + anuR + IIL II
at + aR
Sanle + aLIn araR(uR uL)
II1 = II, =
at + aR
1 1_ aR(uR UL) + IIL IIR
1 1_ a (UR UL) I ng IIL
Pk PR aR L + R)
2 L
tR = tR +, 2a, ,
R~ 2a 2,R
At this stage, the initial states IjL and VR and more pre
cisely the free parameters at and aR are implicitly as
sumed to be such that the waves in the Riemann solu
tions are ordered as they should, namely
us  < S(IR) un+ .~ (18)
PL PR
As (1J)
dtpE + 8,(pEu) = 0.
Initial condition is denoted
Following Bouchut (2004), we define at = aL(VL) and
aR aR(VR) as follows:
if pR > pt
U(x, 0) = Mo(x),
a max(cL, cmin) +a us
PL PRcR
a ma(cR, cmin) + a( + Us
pR aL
ifpR < pm
a ma(cR, cmin) + a( + Us
PR PLcL
max(cL, cmin) + ( +us
PL aR
UR) ,
withEo () = in the case of PGD.
We first set some notations. Let Ax: and at be two
constant steps for space and time discretizations. Let
(xy)jz be a sequence of equidistributed points in RW:
xy I xy ax:. For all j EZ and all ne N we
define
and consider the following discretization of the compu
tational domainlW Rx x R
R:, x R =U c;
jeZn)0
with
C"n [xj/2 Zl/[X FI tn [~
On the one hand and as usual in the context of finite vol
ume methods, the approximate solution Mat,az(x, t) of
UR) ,
with a (*y + 1)/2, cmin > 0 and where ptL,R =
PL,R(ML, a), cL,R cL,R(ML,R) are the values of the
pressures and sound speeds evaluated on UL and us.
This choice has several advantages. First, it is shown to
fullfil (18) and to give the positive of the intermediate
densities pZ and p*,. Then, it complies with the sub
characteristic condition a > pc. At last, it guarantees the
nonlinear stability of the underlying relaxation scheme
that will be described in the next paragraph, and the pos
sibility of handling vacuum in the sense that the speeds
(13) or (10) with initial data uo is sought as a piecewise
constant function on each slab Of:
unat,az(x, t) = u for (x:, t)E 6 O.
At time t = 0, we set
1 5 1/2
Of=
Ax 1/2
IUo(x)dx, jE Z.
On the other hand, we define from ust~ax the piecewise
constant approximate solution l/at,ax by
This solution is set to be at equilibrium, that is
(pH),n = p(up3), jZ E
for the gas dynamics and
(pll), = 0, jZ E
for the PGD.
Let us assume that the solution unt~ax (x, t") at time tn
is known. In order to advance it to the next time level
t" l, we now describe the two steps of the method in
details.
First step: evolution in time (t" t"l)
In this step, we solve (16) with lVat,(zx, tn) as initial
data and for times t E [0, at]. Under the CFL condition
at 1
Smax(Ag(V), i = 1, 2, 3) < (19)
ax: v 2
where the maximum is taken over all the 1/ under con
sideration, the solution is obtained by solving a sequence
of non interacting Riemann problems set at each cell in
terface x, 17. It is explicitly known by the previous
paragraph and we have
for (:, t) 6[xy~x, zj~]x]0, at], jEZ.
We then get back a piecewise constant function in : E
[x 1ye, x 17,] by means of a classical L2 projection'
that is
V~z~1 5 1/2
V~x, ) = (x, t)dx,
A 1/2 '
for (:, t) 6[x 2,x, zj7,~x]0, at], jE Z'
and we set
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Of course, this first step is nothing but the celebrated
Godunov method applied to (16). As a consequence, the
update formula (20) can easily be given the following
conservation form
je Z n '> 0,
where the numerical flux function writes for all j eZ
(21)
Let us recall that the numerical flux (21) is here explic
itly known.
Second step: relaxation (t"+ t" l)
We now project the solution l/at,ax (x, t"+ ) obtained
at the end of the previous step on the equilibrium
manifold p = +oo. More precisely, we set for all j eZ
' ( 3 1 (22)
with
& 1= M and (pll) = ~p(M ~)
in the case of the gas dynamics equations, and
U1"+] (pU)2
i =(p, pu, 2p ) and (pll) = 0
in the pressureless case. This is equivalent to solve in
the asymptotic regime p = +oo
~p(p 11),
in the case of the gas dynamics equations, and
8t ( ) = 0,
at( E) = ApE, (4
at (Pn) = ~p(P 11),
in the case of PGD.
In agreement with the description of these two
Steps, the approximate solution unt~ax is then updated
according to the following consistent finite volume
method:
at
Vn+1U"+13n+l
)
at
(p);+l= (~U)3r nLCnfPU(~S~I:~4~1):
St p = 0,
dt (pu) = 0,
dt (pE) = 0
dt (pTI)
together with
(pE);+1=(pE); A fpE 3 3 1
in the case of gas dynamics and
(pE) = ) +
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
we consider that the cell Cl should be considered with
pressure in the update formula. Thus, we propose to use
At (29)
and
(pE) = (pE)3 Ax nf(ZU, M (30)
where for k = j 1, j, j + 1, M = u" if Y~ = 1 and
in the case of PGD. Here of course,
(f fpu" fpE ) ( pn 13n ) denote the first three
components of g(lj 1/3n z) and
for a = p, pu, pE. Ohrie
Extension to 2D configurations and to secondorder
accuracy
So far, we focused ourselves on the monodimensional
case. In order to perform the 2D computations presented
in the next section on cartesian meshes, we used a very
classical dimensional splitting method. We briefly re
call that it first consists in splitting the twodimensional
governing equations into a pair of quasi onedimensional
equations, and then to solve the underlying sequence of
two onedimensional problems with the proposed nu
merical strategy. Recall that if we denote (u, v) the two
components of the velocity field, v being associated with
the additional space dimension, the governing equation
fOr v in the quasi1D system reads
8t (pv) + Ox (pv) = 0. (31)
This equation means that v is simply passively trans
ported with the flow. From a numerical point of view,
a natural discretisation of (31) is given by
with
a f" (&,9, 4) = fP" (14 4 z)l f 1(up 7l, up)
and
f Pv (upn, 4 n)
S fP(I4,IM) v" if fP(43,M3l )> 0,
fP(3, 3~)V~~ if fP (UPn, M~) < 0
This formula has been first introduced in Larrouturou
(1991) and complies with the exact resolution of the Rie
mann problem for the quasi1D relaxation model
dtP + d,(pu) 0,
88 (pu) + 8,(pU2 + I) = 0,
at (pE) + 8, (pEu + Inu) 0 ,
88 (#11) + d, (plnu + a u) = 0,
dt (pv) + Ox (pv) = 0.
Coupling the gas and pressureless gas dynamics
In order to perform computations involving both the gas
and PGD at the same time, we have to describe how
to couple the relaxation schemes we have developed for
both systems. Recall that the conservative unknowns are
p, pu and pE for the gas dynamics and p and pu for the
PGD. The main difference then clearly lies in the treat
ment of the energy equation.
For the sake of clarity, we begin by introducing a color
function Y such that Y 1 for gas dynamics and Y 0 O
for PGD. From a numerical point of view, a given cell
Cf is said to be pressureless, or equivalently such that
Y," = 0, if the intemal energy Ey = (pE. P 3)3 is
less than a given threshold Emi, and with pressure, that
is Y" 1, otherwise. Introducing the threshold Em "
is a convenient way to switch from one algorithm to the
other. In agreement with the threshold cmi, already in
troduced for the sound speed in the definition of at and
tmin = mi" (28)
Recall indeed that for perfect gas equations of state we
have c2 = (Y1)E. We thus distinguish between zones
with PGD where the internal energy is exactly zero and
zones where the energy level is above the defined small
threshold, a property which is preserved by the pure con
vective part of the evolution.
Let us consider a given cell Cy. Four different situa
tions must be distinguished, depending in particular on
whether Y Y" Y" or not.
The case Yn_ / YP n P In this case, we
simply use (25) and (27) without any modification.
The case Yn I/ / Y"Y In this case, we
simply use (25) and (26) without any modification.
The case Y" _, f Y," and/or YP1 f YP. In this case,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The calculations are left to the reader. The secondorder
extension in space we used in the numerical experiments
is based on a very classical MUSCL reconstruction tech
nique on the primitive variables p, u and =, using a min
mod slope limiter. Regarding the time secondorder ex
tension, we used an usual RungeKutta method coupled
with a Strang splitting.
3 Results and discussion
1D Bouchut test
To evaluate the performance of the present relaxation
scheme in the PGD configuration, the first numerical
test of Bouchut et al. (2003) is performed. In this test,
the initial solution is designed to create a vacuum state
and a mass accumulation. Figure 1 compares the results
between Bouchut's kinetic scheme and the relaxation
scheme for first and second orders. Figure 1.a shows
that vacuum is properly captured by the first and second
order relaxation schemes. As noticed by Bouchut et al.
(2003), the first order scheme forms an artificial density
peak. This problem does not appear for both second or
der schemes. Compared to Bouchut's scheme, the sec
ond order relaxation scheme is slightly more diffusive
and the density overshoots created in the zone of nega
tive velocity divergence are a little bit stronger. All the
schemes capture perfectly the discontinuity of the veloc
ity in the vacuum region (see Fig. 1.b). In the mass ac
cumulation zone, the most accurate results are obtained
with Bouchut's scheme which is closely followed by the
second order relaxation scheme.
1D Sod shock tube
In order to evaluate the hybrid PGD/gas dynamics relax
ation method, the Sod shock tube test is performed with
the following initial conditions:
Figure 1: Profiles of density (a) and velocity (b) for
Bouchut numerical test at time t= 0.5: exact
solution (), 2"d order Bouchut scheme (0),
1'st order relaxation scheme (0), 2"d order
relaxation scheme (+) (80 nodes, CFL = 0.5).
1D Shock/Sshock interaction
The robustness of the secondorder relaxation scheme
is tested in a configuration where a shock propagates
through a pressureless region and meets a 6shock in
density. The 6shock is created by an initial velocity
perturbation located in the pressureless region (see black
line in Fig. 3.b). Figure 3 shows two instants of the cal
culation: before (t = 0.1) and after (t = 0.5) the shock
meets the 6shock. The trace of this interaction on the
velocity profile is an nwave downstream the shock posi
tion (Fig. 3.b). A corresponding density nwave appears
on Fig. 3.a. After having interacted with the 6shock
(t = 0.5), the shock has a higher density ratio than be
fore (t = 0.1) while its velocity jump stays unchanged.
This test demonstrates the high robustness of the present
relaxation scheme.
n 0Pn
v n
1, 1" = 1.1 if .r < 0.5
0.125, " = 0 if .r > 0.5
At the initial time, Jr > 0 corresponds to a zero pres
sure field computed with the pressureless gas algorithm
while Jr < 0 is computed with the gas dynamics algo
rithm. In this test case as in all other coupled method
calculations, en,;, 10" and is evaluated through
Eq. (28). Figures 2 shows the density and pressure pro
files at time t = 0.1644 for the first and second order
relaxation schemes. The interface between pressure and
pressureless regions do not present any numerical arte
fact. Due to the poor discretization of the surface dis
continuity and the shock, the density solution is smeared
by the numerical diffusion. The second order scheme
presents a significantly better accuracy than the first or
der one.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
o
o O
o
0 0.2 0.4 0.6 0.8 1
0 0.25 0.5 0.75
1 1.25 1.5
~O
" 0.4
0
0.4
0.2 0.4 0.6 0.8 1
0 0.25 0.5 0.75
1 1.25 1.5
Figure 2: Profiles of density (a) and pressure (b) for the
hybrid relaxation scheme in Sod numerical
test at time t 0.1644: 1'st order relaxation
scheme (pressure region: 0, pressureless re
gion: +), 2'"d order relaxation scheme (pres
sure region: O, pressureless region: x) (80
nodes, CFL = 0.5).
Figure 3: Second order hybrid relaxation scheme calcu
lation of the interaction of a shock/6shock:
profiles of density (a) and velocity (b) at time
t 0 (), t 0.1 (before the interaction,
pressure region: 0, pressureless region: *)
and t = 0.5 (after the interaction, pressure re
gion: O, pressureless region: +) (120 nodes,
CFL= 0.5).
regimes. For St < St,, the particles cannot escape from
the TaylorGreen vortices while, for St > St,, they are
ejected out of their original vortices. Therefore, the fol
lowing tests consider two values of St in order to cover
these two regimes: St 0.9St, and St 13Stc. From
a numerical point of view, the drag source term F is ap
plied via operator splitting through an analytical expres
sion of the exponential relaxation (Eq. 10). The initial
spray velocity is uniformly zero for all testcases.
Pressureless transport of a nonuniform initial dis
tribution at supercritical Stokes number. In order to
test the capability of the method to treat multidimen
sional transport of inertial particles, the Stokes number
is fixed at a supercritical value St = 13Stc. Figure 4.a
2D TaylorGreen vortices
Figure 4.a shows the velocity field U = (U, V) of the
carrier phase corresponding to the four contrarotating
TaylorGreen vortices used in the following numerical
tests:
( U(.r, y)
sin(2x~r) cos(2ixy)
~~~~~~ co(iz i(i
The spray dynamics is coupled to the gaseous flow field
through a Stokes drag source term in the momentum
equation, which amounts to relaxing the spray velocity
field toward the gaseous one at a rate set by the Stokes
number St, i.e. the non dimensional relaxation time.
From de ( Ilndisonl.elllion (2009) we know that there ex
ists a critical value St, = 1/8x7 which separates two
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
d). Most of the density is concentrated close to y = 0
in the PGD case, while a smoother density distribution
is observed in the hybrid case. Note that with the hybrid
scheme, the maximum of density is not located on the
y = 0 line because the pressure gradient resulting from
the energy source term prevents particles from accumu
lating there.
Conclusions
A novel hybrid numerical method for solving Eulerian
models for spray dynamics has been proposed. Based
on the relaxation method, it can deal with both PGD
and general gas dynamics system of equations in vari
ous zones of the same configuration. Therefore, it has
the ability, on the one hand, to compute the lowinertia
particles dynamics described by PGD and, on the
Other hand, to account for the effects of highinertia par
ticles in the turbulent regions of the flow falling under
the general gas dynamics framework. The zerodensity
is also explicitly handled, which is a key feature for
simulating spray injection. In terms of accuracy, one
and twodimensional tests in PGD configurations show
that the scheme matches the kinetic scheme of Bouchut
previously used and thus validate the approach. Beside,
the hybrid PGD/gas dynamics approach predicts accu
rate results in the 1D shock tube testcase. The high
robustness of the method is demonstrated, in particular
in a shock/6shock interaction. Twodimensional simu
lations in the framework of TaylorGreen vortices with
eventually localized turbulent subgrid energy source al
low to exhibit the potential of the method. Besides, the
relaxation framework makes it possible to handle arbi
trary pressure law such as the real gastype behaviour
of turbulent sprays. Therefore the present investigation
shows that this method has the ingredients needed to
simulate turbulent sprays in a DNS/LES framework.
Acknowledgements
Quang Huy Tran from IFP is gratefully acknowledged
for his precious help to the present work.
References
C. Berthon, M. Breu8, and M.O. Titeux. A relaxation
scheme for the approximation of the pressureless Eu
ler equations. Nwner. Methods for Partial Dliffeiva,rial
Equations, 22(2):484505, 2006.
G. A. Bird. Molecular gas dynamics and the direct siln
ulation of gas flows, volume 42 of Oxford Engineering
Science Series. Oxford University Press, 1994.
shows the initial density distribution provided by a cardi
nal sinus function. To allow comparison with Bouchut's
scheme, the pressureless relaxation scheme is used. Fig
ure 4.b and 4.c show the results at time t = 0.8 for
the second order Bouchut scheme and the second or
der relaxation scheme. Both schemes predict very sim
ilar density fields featuring a 6shock, as expected (see
de Chaisemartin (2009)).
Pressure vs pressureless transport of an uniform
initial distribution at subcritical Stokes number. Fig
ures 5.a and 5.b compares the density fields obtained
with the second order Bouchut scheme and the second
order relaxation scheme respectively, starting from a
uniform density distribution in subcritical Stokes num
ber conditions St 0.9Stc. As for the supercritical
case, both scheme predicts very similar results, with
mass concentrating in the high strain regions because
of the ejection of particles from the center of vortices.
The hybrid second order scheme still being under de
velopment, calculations of the standard gas dynamics
relaxation scheme are presented only for the first order
scheme. In order to limit the effect of different numeri
cal diffusion according to the scheme order, a spatial res
olution of two times larger is used for the gas dynamics
first order scheme (400 nodes vs. 200 nodes for the sec
ond order calculations). Figure 5.c shows that pressure
effects limits the segregation of particles (density con
centration here) because. Physically, this can be inter
preted as a mixing effect from the fluid turbulence. Here,
this turbulence is very simply modeled by a uniform re
laxation term of internal energy (target value E, 1 ).
Hybrid pressure/pressureless transport of a non
uniform initial distribution at supercritical Stokes
number. The hybrid scheme is evaluated in a config
uration where two parcels of highinertia particles are
ejected from their initial vortex and collides together.
The density distribution is given by a cardinal sinus
function whose center is (0.125, 0.375) and radius
0.125 (see Fig. 6). A PGD calculation and a gas dy
namics calculation are performed. For the gas dynamics
case, a relaxation term of internal energy is imposed us
ing a cardinal sinus function centered on y 0 with a
radius of 0.125 and a maximum value of Et = 0.5 (see
Fig. 6.b). This energy source term simulates the effect of
a local turbulence region of the carrier flow on the panti
cles transport. As expected, Fig. 6 shows that, for both
PGD and hybrid schemes, each particles parcel is ej ected
from its vortex and start to interact with its mirror image
at t 0.75. In the PGD case (Fig. 6.a), this interaction
forms a 6shock at the meeting line y 0 On the other
hand, the hybrid case (Fig. 6.b) features only a small in
crease in density at y = 0 because pressure effects limit
the particles concentration. Moments later, the behavior
of both schemes are even more different (see Fig. 6.c and
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
M. Garcia. Development and validation of the Euler
Lagrange foundation on a parallel and unstructured
solver for largeeddy simulation. PhD thesis, Institut
National Polytechnique de Toulouse, available in english
at http://tel.archivesouvertes.fr/tel0041407e/ 2009.
J. B. Greenberg, I. Silverman, and Y. Tambour. On the
origin of spray sectional conservation equations. Corn
bust. Flame, 93:9096, 1993.
S. Jin and Z. Xin. The relaxation schemes for systems of
conservation laws in arbitrary space dimensions. Coini.
Pure Appl. Math., 48(3), 1995.
D. Kah, F. Laurent, L. Frdret, S. de Chaisemartin, R.O.
Fox, J. Reveillon, and M. Massot. Eulerian quadrature
based moment models for dilute polydisperse evap
orating sprays. Flow Turbul. Cornbust., accepted
for publication, available on HAL http://hal.archives
ouvertes.fr/hal00449866/en/ 2010.
A. Kaufmann, M. Moreau, O. Simonin, and J. Helie.
Comparison between Lagrangian and mesoscopic Eu
lerian modelling approaches for inertial particles sus
pended in decaying isotropic turbulence. J. Comput.
Phys., 227(13):6448 6472, 2008.
G. Lacaze, B. Cuenot, T. Poinsot, and M. Oschwald.
Large eddy simulation of laser ignition and compress
ible reacting flow in a rocketlike configuration. Corn
bust. Flame, 156(6):11661180, 2009.
B. Larrouturou. How to preserve the mass fractions pos
itivity when computing compressible multicomponent
flows. J. Comput. Phys., 95:5984, 1991.
F. Laurent and M. Massot. Multifluid modeling of
laminar polydispersed spray flames: origin, assump
tions and comparison of sectional and sampling meth
ods. Comb. Theory and Modelling, 5:537572, 2001.
M. Massot, R. Knikker, C. P~ra, and J. Reveillon. La
grangian/Eulerian analysis of the dispersion of evaporat
ing sprays in nonhomogeneous turbulent flows. In Proc.
5th Int. Conf Multiphase Flow, 2004.
M. Massot, F. Laurent, S. de Chaisemartin, L. Frdret,
and D. Kah. Eulerian multifluid models: modeling
and numerical methods. In Modelling and Computa
tion of Nanoparticles in Fluid Flows, Lecture Notes of
the von Karman Institute. NATO RTO AVT169, 2009.
In Press, Available at http://hal.archivesouvertes.fr/hal
00423031/en/.
M. Moreau, B. Bedat, and O. Simonin. A priori test
ing of subgrid stress models for eulereuler twophase
LES from eulerlagrange simulations of gasparticle tur
bulent flow. In 18th Ann. Conf on Liquid Atornization
and Spray Systems. ILASS Americas, 2005.
M. Boileau, S. Pascaud, E. Riber, B. Cuenot, L.Y.M.
Gicquel, T. Poinsot, and M. Cazalens. Investigation of
twofluid methods for Large Eddy Simulation of spray
combustion in Gas Turbines. Flow Turbul. Cornbust., 80
(3):291321, 2008a.
M. Boileau, G. Staffelbach, B. Cuenot, T. Poinsot, and
C. B~rat. LES of an ignition sequence in a gas turbine
engine. Cornbust. Flame, 154(12):222, 2008b.
F. Bouchut. Nonlinear stability offinite volume methods
for hyperbolic conservation laws, and wellbalanced
schemes for sources. Frontiers in Mathematics series,
2004.
F. Bouchut, S. Jin, and X. Li. Numerical approximations
of pressureless and isothermal gas dynamics. SIAM J.
Nwner Anal., 41(1):135158, 2003.
C. Chalons and F. Coquel. Navierstokes equations with
several independent pressure laws and explicit predictor
corrector schemes. Numer Math., 101(3):451478,
2005.
C. Chalons and J.F. Coulombel. Relaxation approxima
tion of the euler equations. J. Math. Anal. Appl., 348(2):
872893, 2008.
F. Coquel and B. Perthame. Relaxation of energy and
approximate riemann solvers for general pressure laws
in fluid dynamics. SIAM J. Numer Anal., 35(6):2223
2249, 1998.
S. de Chaisemartin. Modules enle'riens et simula
tion nuindrique de la dispersion turbulente de brouil
lards qui s'dvaporent. PhD thesis, Ecole Centrale
Paris, France, available in english at http://tel.archives
ouvertes.fr/tel00443982/en/, 2009.
S. de Chaisemartin, L. Frdret, D. Kah, F. Laurent, R.O.
Fox, J. Reveillon, and M. Massot. Eulerian models
for turbulent spray combustion with polydispersity and
droplet crossing. C R. Mi'canique, 337:438448, 2009.
F. Duchaine, T. Morel, and L.Y.M. Gicquel.
ComputationalFluidDynamicsBased Kriging Opti
mization Tool for Aeronautical Combustion Chambers.
AIAA J., 47(3), 2009.
P. Fevrier, O. Simonin, and K. Squires. Partitioning
of particle velocities in gassolid turbulent flows into
a continuous field and a spatially uncorrelated random
distribution theoretical formalism and numerical study.
J.Fluid Mech, 533:146, 2005.
L. Frdret, S. de Chaisemartin, J. Reveillon, F. Laurent,
and M. Massot. Eulerian models and threedimensional
numerical simulation of polydisperse sprays. In Proc.
7th Int. Conf Multiphase Flow, Tampa, FL, USA, 2010.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
P. J. O'Rourke. Collective drop aT i, r\ on vaporizing
liquid sprays. PhD thesis, Princeton University, 1981.
R.V.R. Pandya and F. Mashayek. Twofluid largeeddy
simulation approach for particleladen turbulent flows.
Int. J. Heat Mass Transfer, 45(24):47534759, 2002.
M.W. Reeks. On the continuum equations for dispersed
particles in nonuniform flows. Phys. Fluids A, 4(6):
12901303, 1992.
E. Riber, M. Moreau, O. Simonin, and B. Cuenot. To
wards large eddy simulation of nonhomogeneous par
ticle laden turbulent gas flows using eulereuler ap
proach. In 11th Workshop on TwoPhase Flow Predic
tions, Merseburg, Germany, 2005.
E. Riber, V. Moureau, M. Garcia, T. Poinsot, and O. Si
monin. Evaluation of numerical strategies for large eddy
simulation of particulate twophase recirculating flows.
J. Comput. Phys., 228(2):539564, 2009.
A. Roux, L.Y.M. Gicquel, Y. Sommerer, and T. Poinsot.
Large eddy simulation of mean and oscillating flow in
sidedump ramjet combustor. Combust. Flame, 152(1
2):154176, 2008.
L. Selle, L. Benoit, T. Poinsot, F. Nicoud, and W. Krebs.
Joint use of compressible largeeddy simulation and
helmholtz solvers for the analysis of rotating modes in
an industrial swirled burner. Combust. Flame, 145(12):
194205, 2006.
O. Simonin, P. Fevrier, and J. Lavieville. On the spa
tial distribution of heavyparticle velocities in turbulent
flow: from continuous field to particulate chaos. J. Tur
bul.,3(40):118, 2002.
I. Suliciu. Energy estimates in ratetype thermo
viscoplasticity. Int. J. Plast., 14:227244, 1998.
A. Vi6, M. Sanjos6, S. Jay, C. Angelberger, B. Cuenot,
and M. Massot. Evaluation of a multifluid mesoscopic
eulerian formalism on the Large Eddy Simulation of an
aeronauticaltype configuration. In Proc. 7th Int. Conf
Multiphase Flow, Tampa, FL, USA, 2010.
F. A. Williams. Spray combustion and atomization.
Phys. Fluids, 1:541545, 1958.
L.I. Zaichik. A statistical model of particle transport and
heat transfer in turbulent shear flows. Phys. Fluids, 11:
1521, 1999.
L.I. Zaichik, O. Simonin, and V.M. Alipchenkov. An
Eulerian approach for large eddy simulation of parti
cle transport in turbulent flows. J. Turbul., 10(4):121,
2009.
1
0.6
0.4 z.
0.2 l
0.2 .4 06 0.
0.8.406 .
0.2 0.4 0.6 0.8
X
0.2 0.4 0.6 0.8
X
Figure 4: Carrier phase velocity field (TaylorGreen pe
riodic vortices) and initial density contours
(a). Snapshots of the density distribution at
time t 0.8 for Stokes number St 13St,
(200 nodes, CFL 1 ): 2nd order Bouchut
scheme (b) and 2"d order pressureless relax
ation scheme (ck.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
a) 0.4
>0
0.2 0.4 0.6 0.8
X
0.2
>
0.2 0.4 0.6 0.8
X
>
0.2
...... ..:...j i
0.2 0.4 0.6 0.8
X
0.4 0.2 0
x
0.2 0.4
Figure 5: Snapshots of the density distribution at time
t 0.33 for Stokes number St 0.9Ste: 2""
order Bouchut scheme (a), 2"d order pres
sureless relaxation scheme (b) and 2"d or
der relaxation scheme with pressure (c). (a,
b): 200 nodes, CFL = 0.5. (c): 400 nodes,
CFL = 1.
Figure 6: Snapshots of the density for Stokes number
St = 13St, (1rst order relaxation scheme, 200
nodes, CFL 1 ): pressureless scheme (a, c)
and hybrid scheme (b, d). Time t =0.75 (a, b)
and t = 1.1 (c, d). Dotted circles and dashed
lines are the limits of the cardinal sinus func
tion of the initial density distribution and the
energy source term (b, d) respectively.
