Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 14.4.2 - Robust numerical schemes for Eulerian spray DNS and LES in two-phase turbulent flows
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00353
 Material Information
Title: 14.4.2 - Robust numerical schemes for Eulerian spray DNS and LES in two-phase turbulent flows Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Boileau, M.
Chalons, C.
Bourgouin, J.-F.
Terrier, C.
Laurent, F.
de Chaisemartin, S.
Massot, M.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: liquid sprays
multi-fluid models
spray equation
LES
relaxation schemes
 Notes
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); Fréret 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.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00353
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1442-Boileau-ICMF2010.pdf

Full Text



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 two-phase
turbulent flows


M. Boileau*, C. Chalonsi*, J.-F. Bourgouin*, C. Terrier*, F. Laurent*,

S. de Chaisemartin~ and M. Massot*

Laboratoire EM2C-UPR CNRS 288, Ecole Centrale Paris, 92295 Chiltenay-Malabry, France
t DEN/DANS/DM2S/SFME/LETR CEA-Saclay, 91191 Gif-sur-Yvette, 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.de-chaisemartin~ifp.fr and mare.massot~em2c.eep.fr (corresponding author)
Keywords: Liquid Sprays; Multi-Fluid 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 two-phase 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 non-zero dispersion in the framework of ensemble
averages and modeling non-resolved 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 higher-inertia particles induces a
random-uncorrelated 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 random-uncorrelated 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
semi-kinetic 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 so-called 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 high-dimensional. 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 droplet-droplet 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 Lagrangian-Monte-Carlo approach
(see O'Rourke (1981)). It is called Direct Simulation
Monte-Carlo 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 Multi-Fluid solver capable of de-
scribing the polydispersity of a spray in size and the as-
sociated size-conditioned dynamics. This approach re-
lies on the derivation of a semi-kinetic 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 test-cases 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 dynamics-type system of con-
servation equations. The hybrid relaxation scheme is de-
rived in the second section. Finally, results for relevant
test-cases, 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 du-neighborhood 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 multi-fluid 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 multi-fluid 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 delta-shocks 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 Euler-type system of
equations in regions of high Stokes number or in re-
gions where sub-grid 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
non-linearity 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 non-linearity 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 mono-kinetic 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 mono-kinetic 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
mono-kinetic 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) = nF-u
2 n69 + add. terms

where the total energy reads I u-u/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 so-called 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 3-1







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 so-called sub-characteristic condition a > pc
where c denotes the sound speed.

The first-order system extracted from (15) is hyperbolic
and admits the following three eigenvalues,

a a

with second-order 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 Godunov-type 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 self-similar 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 second-order
accuracy

So far, we focused ourselves on the mono-dimensional
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 two-dimensional
governing equations into a pair of quasi one-dimensional
equations, and then to solve the underlying sequence of
two one-dimensional 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 quasi-1D 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 quasi-1D 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 = (Y-1)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 second-order
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 second-order ex-
tension, we used an usual Runge-Kutta 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/S-shock interaction

The robustness of the second-order relaxation scheme
is tested in a configuration where a shock propagates
through a pressureless region and meets a 6-shock in
density. The 6-shock 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 6-shock. The trace of this interaction on the
velocity profile is an n-wave downstream the shock posi-
tion (Fig. 3.b). A corresponding density n-wave appears
on Fig. 3.a. After having interacted with the 6-shock
(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/6-shock:
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 Taylor-Green 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 test-cases.
Pressureless transport of a non-uniform 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 Taylor-Green vortices


Figure 4.a shows the velocity field U = (U, V) of the
carrier phase corresponding to the four contra-rotating
Taylor-Green 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 low-inertia
particles dynamics described by PGD and, on the
Other hand, to account for the effects of high-inertia par-
ticles in the turbulent regions of the flow falling under
the general gas dynamics framework. The zero-density
is also explicitly handled, which is a key feature for
simulating spray injection. In terms of accuracy, one
and two-dimensional 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 test-case. The high
robustness of the method is demonstrated, in particular
in a shock/6-shock interaction. Two-dimensional simu-
lations in the framework of Taylor-Green 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 gas-type 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):484-505, 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 6-shock, 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 high-inertia 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 6-shock 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 large-eddy simulation. PhD thesis, Institut
National Polytechnique de Toulouse, available in english
at http://tel.archives-ouvertes.fr/tel-0041407e/ 2009.

J. B. Greenberg, I. Silverman, and Y. Tambour. On the
origin of spray sectional conservation equations. Corn-
bust. Flame, 93:90-96, 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/hal-00449866/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 rocket-like configuration. Corn-
bust. Flame, 156(6):1166-1180, 2009.

B. Larrouturou. How to preserve the mass fractions pos-
itivity when computing compressible multi-component
flows. J. Comput. Phys., 95:59-84, 1991.

F. Laurent and M. Massot. Multi-fluid modeling of
laminar poly-dispersed spray flames: origin, assump-
tions and comparison of sectional and sampling meth-
ods. Comb. Theory and Modelling, 5:537-572, 2001.

M. Massot, R. Knikker, C. P~ra, and J. Reveillon. La-
grangian/Eulerian analysis of the dispersion of evaporat-
ing sprays in non-homogeneous turbulent flows. In Proc.
5th Int. Conf Multiphase Flow, 2004.

M. Massot, F. Laurent, S. de Chaisemartin, L. Frdret,
and D. Kah. Eulerian multi-fluid models: modeling
and numerical methods. In Modelling and Computa-
tion of Nanoparticles in Fluid Flows, Lecture Notes of
the von Karman Institute. NATO RTO AVT-169, 2009.
In Press, Available at http://hal.archives-ouvertes.fr/hal-
00423031/en/.

M. Moreau, B. Bedat, and O. Simonin. A priori test-
ing of subgrid stress models for euler-euler two-phase
LES from euler-lagrange simulations of gas-particle 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
two-fluid methods for Large Eddy Simulation of spray
combustion in Gas Turbines. Flow Turbul. Cornbust., 80
(3):291-321, 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(1-2):2-22, 2008b.

F. Bouchut. Nonlinear stability offinite volume methods
for hyperbolic conservation laws, and well-balanced
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):135-158, 2003.

C. Chalons and F. Coquel. Navier-stokes equations with
several independent pressure laws and explicit predictor-
corrector schemes. Numer Math., 101(3):451-478,
2005.

C. Chalons and J.-F. Coulombel. Relaxation approxima-
tion of the euler equations. J. Math. Anal. Appl., 348(2):
872-893, 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/tel-00443982/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:438-448, 2009.

F. Duchaine, T. Morel, and L.Y.M. Gicquel.
Computational-Fluid-Dynamics-Based 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 gas-solid turbulent flows into
a continuous field and a spatially uncorrelated random
distribution theoretical formalism and numerical study.
J.Fluid Mech, 533:1-46, 2005.

L. Frdret, S. de Chaisemartin, J. Reveillon, F. Laurent,
and M. Massot. Eulerian models and three-dimensional
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. Two-fluid large-eddy
simulation approach for particle-laden turbulent flows.
Int. J. Heat Mass Transfer, 45(24):4753-4759, 2002.

M.W. Reeks. On the continuum equations for dispersed
particles in nonuniform flows. Phys. Fluids A, 4(6):
1290-1303, 1992.

E. Riber, M. Moreau, O. Simonin, and B. Cuenot. To-
wards large eddy simulation of non-homogeneous par-
ticle laden turbulent gas flows using euler-euler ap-
proach. In 11th Workshop on Two-Phase 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 two-phase recirculating flows.
J. Comput. Phys., 228(2):539-564, 2009.

A. Roux, L.Y.M. Gicquel, Y. Sommerer, and T. Poinsot.
Large eddy simulation of mean and oscillating flow in
side-dump ramjet combustor. Combust. Flame, 152(1-
2):154-176, 2008.

L. Selle, L. Benoit, T. Poinsot, F. Nicoud, and W. Krebs.
Joint use of compressible large-eddy simulation and
helmholtz solvers for the analysis of rotating modes in
an industrial swirled burner. Combust. Flame, 145(1-2):
194-205, 2006.

O. Simonin, P. Fevrier, and J. Lavieville. On the spa-
tial distribution of heavy-particle velocities in turbulent
flow: from continuous field to particulate chaos. J. Tur-
bul.,3(40):1-18, 2002.

I. Suliciu. Energy estimates in rate-type thermo-
viscoplasticity. Int. J. Plast., 14:227-244, 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
aeronautical-type configuration. In Proc. 7th Int. Conf
Multiphase Flow, Tampa, FL, USA, 2010.

F. A. Williams. Spray combustion and atomization.
Phys. Fluids, 1:541-545, 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):1-21,
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 (Taylor-Green 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.




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

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