7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A high order moment method with mesh movement for the description of a
polydisperse evaporating spray
D. Kah*t, M. Massott, Q.H. Tran *, S. Jay *, F. Laurent t
IFP, 1 et 4 Avenue de BoisPreau, 92852 RueilMalmaison, FRANCE
Laboratoire EM2C, CNRS, Grande Voie des Vignes, 92295 ChatenayMalabry, FRANCE
damien.kah@ifp.fr, marc.massot@em2c.ecp.fr, qhuy.tran@ifp.fr, stephane.jay@ifp.fr, laurent@em2c.ecp.fr,
Keywords: Eulerian models, polydispersivity, high order moment method, ALE formalism
Abstract
This work investigates the issue of describing polydispersivity in an Eulerian framework of a disperse spray
with potentially mesh movement. In this perspective, the multifluid model and and associated numerical schemes
provide robust tools (de Chaisemartin (2009)). However, there is a substantial interest for the development of a
method able to account for the size distribution of particles without discretizing the size phase space into sections
such as in the classical multifluid model. The objective is to write a polydisperse spray model in a formulation
close to the usual twofluid one. Eulerian high order size moment methods are well suited for such problems, but
the design of numerical algorithms usually faces two difficulties: accuracy and stability. Whereas a multifluid
approach with a single section would be too coarse, standard twofluid models only provide an approximated mean
diameter which is also too coarse an assumption in order to capture polydispersivity. This issue has been successfully
tackled considering a high order size moment method. The vector of successive size moments of a number density
function (NDF) over a given size interval belongs to what is called the moment space. But its complex geometry
makes difficult to ensure that a moment vector stays in the moment space during the computation. The achievement
of this work consists in designing numerical schemes to address these issues. These methods are displayed in the
context of evaporation (Massot et al. (2010)) and advection (Kah et al. (2010)) of a dispersed droplet flow. These
results are the first step, for our application of interest, of the simulation of an injected polydisperse spray into a
combustion chamber. Indeed, downstream the injector, after primary and secondary atomization, the spray consists
in a polydisperse droplet population. An accurate description of the droplet size distribution will improve precision
on the fuel fraction evaporated in the gaseous phase. These injection cases involve a mobile mesh to deal with the
engine operating conditions. Arbitrary Lagrangian Eulerian (ALE) methods provide numerical schemes adapted for
problems involving interactions between fluids and a moving structure. In this paper, we present a new formalism
which allows to keep the conservation properties of the numerical scheme for the evolution of the moments in the
ALE formalism, in order to preserve the moment space. This new developed tool is validated on academic test cases,
showing its feasability in the perspective of its implementation in an unstructured flow solver, IFPC3D (Bohbot et al.
(2009)), which is in progress.
1 Introduction
In order to provide an accurate numerical tool for the
prediction of engines performance, the simulation of liq
uid fuel jets in combustion chambers is necessary. In this
perspective, the ability to simulate multiphase flow is of
crucial importance. Indeed fuel is stored as a liquid and
injected at high pressure (up to 2000 bars) in a chamber
filled with gas. Right before ignition, we can distinguish
two different regimes. The fluid structure near the injector
is called separate phase, whereas downstream the injector,
the fuel has the shape of a disperse droplet cloud. Thus
different models are used to describe each area. The sep
arate phase can only be modeled by twofluid approaches
whereas, even if the twofluid approach can also describe
the disperse phase, the particulate aspect of the disperse
phase makes kineticbased approaches more relevant.
This paper focuses on the disperse phase and displays
a model and a numerical tool to access polydispersivity
without having to cope with a costful discretization of
size phase space associated to classical multifluid mod
els (Greenberg et al. (1993), de ( ii.mcini.iii (2009)).
This model can be applied to the class of sprays or the
class of aerosols. The cloud is described by a statistical
model involving a number density function (NDF) solu
tion of a WilliamsFokkerPlanck equation. This equa
tion contains all the transport terms in the physical and
the phase space (containing the droplet velocities, sizes,
...), in order to describe convection, drag, evaporation,
and also particleparticle interaction terms such as col
lision, coalescence. However, its numerical solution in
multidimensional systems by finite volume schemes is
intractable due to the high phase space dimension. Con
sequently, Lagrangian approaches have been widely used
and have been shown to be efficient in numerous cases, for
example in (Fox et al. (2008)) and references herein. But
its main drawback is the coupling of an Eulerian descrip
tion for the gaseous phase to a Lagrangian description of
the disperse phase, thus offering limited possibilities of
vectorization/parallelization. Besides, as in any statistical
approach, Lagrangian methods require a relatively large
number of parcels to control statistical noise, and thus are
computationally expensive. The computational cost gets
even higher when unsteady flows are considered.
This drawback makes the use of an Eulerian formula
tion for the description of the disperse phase attractive.
This method consists in solving transport equations not
for the NDF directly, but for selected moments of the ki
netic equation using a moment method. The use of mo
ment methods leads to the loss of some information but
the cost of such methods is usually much lower than the
Lagrangian ones for two reasons. The first is related to the
fact that the number of unknowns we solve for is limited;
and the second is related to the high level of optimization
one can reach when the two phases are both described by
the Eulerian model. The multifluid model only consid
ers one size moment which accounts for the liquid mass
on small intervals of the size phase space called sections.
This has been shown to yield simple transport algorithms
for transport in physical space in (de Chaisemartin (2009))
easily implemented on parallel architectures (Freret et al.
(2010)). However, the cost of the discretization in size
phase space is high. Moreover, its implementation in an
existing CFD code with an architecture not adapted for a
sectional method requires to completely rethink this lat
ter. Consequently the possibility of high order moment
method considering only one size section is very attrac
tive.
Some high order moment methods have already been de
signed. The first one consists in solving the evolution of
moments of a presumed NDF (assumed as a lognormal
law) (Mossa (2005)). But this approach leads to seri
ous numerical instabilities in the treatment of evapora
tion. The second method uses Direct Quadrature Method
of Moment (DQMOM) (Marchisio et al. (2005)) where
equations on nodes and abscissas of the distribution func
tion are directly written (Fox et al. (2008)). Contrary to
the first one, this method is stable, but its accuracy can
still be improved. A new high order moment method, ex
plained in (Massot et al. (2010)), provides a stable nu
merical scheme for evaporation, and reaches unequivalent
levels of accuracy. In (Kah et al. (2010)), a stable scheme
for the transport in physical space of these moments is
added.
These methods, by solving theoretical problems, have
opened new perspectives in the simulation of polydis
persivity, and have proved to be a serious alternative to
the multifluid model, without the constrain of discretiza
tion in size sections. A simulation of an evaporating
free jet is presented in order to assess the feasability of
this high order size moment method. The computation is
lead with the code MUSES3D (de ( Ii.iii.iiiiilll (2009))
where this method has been implemented. But so far,
they have been designed on fixed Eulerian grids and only
for academic test cases. First, the applicability of such
a formalism to CFD still needs to be investigated. Be
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
sides, for engine simulations however, fuel injection in
combustion chamber is done in the more general con
text of moving geometries, involving a moving piston.
To deal with the interaction between the piston and the
fluid, and the computational area change with time, the
Arbitrary LagrangianEulerian (ALE) formalism (Donea
(2004)) has been adopted in the design of the numerical
scheme. Indeed the ALE formalism is adapted to fluid
structure interactions as it combines advantages from the
Lagrangian and the Eulerian representations. These meth
ods are enforced in the context of engine simulation in the
unstructured CFD code IFPC3D (Bohbot et al. (2009)).
In order to introduce the high order moment method de
scribed above in C3D, two challenges have to be taken up.
The first one is to adapt a high order moment method on
a compact subset in the ALE formalism. Given the sta
bility problems of high order methods on Eulerian grids
already, this is an important task, and is first done in a
one dimensional grid in order to isolate the theoretical
problems. Therefore, in this paper, is provided a stable
numerical scheme able to cope with a moving mesh, for
evaporation and convection of the particles.
The paper is organized as follows. In the next section,
the equation system on moments is derived in an Eulerian
framework and the key properties of the moment space
are highlighted to design a robust numerical scheme for
the advection of moments. The resolution of a polydis
perse evaporating free jet in 2D with the high order size
moment method is displayed. Then the equation system
for the moments and for the gas phase are written in a
moving frame. In section 3, the ALE numerical scheme
is displayed as well as the second order in space and time
scheme for the advection of the moments. These develop
ments are validated in the context of the Riemann prob
lem for aerosols and advection with discontinuous veloc
ity fields for sprays. Finally some test cases validate the
implementation of polydispersivity in C3D, which is still
in progress: compression of a high pressure cell by a pis
ton for aerosols, and the Riemann problem for sprays.
2 Moment equation and property
of the moment space
2.1 Kinetic model and Derivation of the
moment equations
We take the general point of view of a dilute
droplet flow described by a NDF f(t, x, S, v), such that
f(t, x, S, v)dxdSdv represents the probable number of
particles located in x (xz,..., a)t, where d is the di
mension of the physical space, with size S, and veloc
ity v. As the particles are assumed to be spherical, we
choose the particle surface to describe the size. We could
have also chosen to work with their radius, r or their
volume, V, all linked by the relation f(t, x, r, v)dr
f(t, x, S, v)dS = f(t, x, V, v)dV. In the applications
we study in this paper, particles undergo evaporation and
collision with the molecules of gas. For the purpose of
this paper, evaporation is solved through a d2 law, but we
remark that the treatment of more complex laws is done
in (Massot et al. (2010)). The collision term with the gas
molecules is expressed with a Stokes drag.
The kinetic equation verified by the NDF is a Williams
FokkerPlanck based equation:
a9f+Vd(vf)Gs(Rf) = V,[((vu,)f + Q f)] (1)
The second term of lhs represents transport in physical
space whereas the third term of the lhs accounts for
the evaporation of the droplets. The term of the rhs
stands for the drag force (for which the expression given
by Stokes is considered) and Brownian motion due to
the interaction of the droplets with the gas molecules.
The evaporation rate, R, is the rate of the decrease of
the droplet surface, 3 1 is the relaxation time of the
particles to the gas velocity ug, Q is the matrix of the
temporal correlation of forces due to brownian motion. If
F (FI, F2, F3)t is the force due to Brownian motion
in three dimensions, then Qij(t) = o Fi(F(t)F(t+T) dT.
In parallel, the gas is described by the Euler system
for compressible gas dynamics.
We consider only a oneway coupling, that is to say the
influence from the gas on the particles. We would like
to emphasize that this is not a limitation of the model,
but we focus only on the description of the liquid droplet
cloud.
In what follows we use dimensionless variables:
t XI S' V
ST Lo St
f,=f O Lo St
no Uo
RLo 1 Q
K SoUo' u2
V Ug
Uo 1
Lo'
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where Dr. = Ot. + UgV..
Our work focuses on the transport and the evaporation
terms. The diffusion term will not be solved in this paper,
but one remarks that it can be treated with a dedicated
solver in the context of a splitting algorithm. Moreover
our purpose is to design a robust moment method for the
particle. So, in a first approach, the gas acceleration term
can be omitted.
In the case of spray particles, their Stokes number is
too important for their trajectory to be influenced by the
small scale fluctuations of the Brownian motion, but are
still coupled with the gas through drag. We assume that
the term Q//3, f is negligible in comparison with the
term (v ug)f. Thus f verifies the following kinetic
equation, called WilliamsBoltzmann equation (Williams
(1958)):
9tf + V,(vf) 9s(Kf) + Vv.(,(u v)f) = 0 (4)
Contrary to the previous case a hl p1 ,Ili, on the veloc
ity distribution has to be made in order to obtain a closed
equation system on the total number density n(t, x, S)
and the mean particle velocity u(t, x, S)) conditioned
by size. This system is called the semikinetic system.
Following the example of what is done in the multi
fluid model (de Chaisemartin (2009)), we project f on
a distribution with a single velocity conditioned by size:
f(t, x, S, v) = n(t, x, v)6(v u(t, x, S)), where u is
the mean velocity of the particles.
The semikinetic equation is obtained by taking the veloc
ity moments of order 0 and 1 of Eq. (4):
1 ,
9tn+ V,(nu) 9s(Kn)
ctnu+Vct(nnu (9) u)9s(Knu, r(ug
St
) =0
(5)
where Lo is a reference length, Uo is a reference velocity,
So is the maximum droplet size, no is a reference number
density. The quantity T7 is the characteristic time of the
particles, St is the Stokes number, K is the dimensionless
evaporation rate, o is the velocity dispersion of the
particles due to Brownian motion.
Written in terms of dimensionless quantities and omit
ting the prime sign for the sake of legibility, the Williams
FokkerPlanck equation reads:
1
9t f + V(vf)) s(Kf) V= V [((v g) f + aVf)]
(2)
The low Stokes number of particles belonging to the
aerosol class makes the rhs appear as a singular pertur
bation for Eq. (1). The NDF is then solved through
a ChapmanEnskog development (Cowling et al. (1970),
Kaper et al. (1972)). After some algebra, we find that the
number density n(t, x, S) = jR f dv verifies the Smolu
chowski equation (Chandrasekhar (1943)):
9tn + ,V(i .) 9s(Kn) St[V, (aVn)
nDtug],
In both cases, we write the resulting system of equa
tions on the variables of interest, which are the size mo
ments of order 0 ( total number density), 1 (giving the
mean size), 2 ( giving dispersion around the mean size)
up to N on the non dimensional interval [0, 1] are defined
as follows,
mk j Skn(t,x,S)dS
Taking the N + 1 first size moments of system (5) leads
to the independent equations:
9t(mTo)+
at(mn)+
'9t(Tn2)+
Vt(monu)
Vt (mn i)
Ve' (mn21)
9t(mN)+ Vt (mNNu)
Kn(t, x,0)
Kmo
2Kmi
NKmNiI
mo (^o u)
t (miu)+V (miu u) =o(9
m\ u
We consider m1u as the momentum and not monu,
because, as St depends on S, the quantity gt would
have not been integrable. In this context, 0 is such as
m (t 0) 0 St. The subsystem corresponding to the
mo
first N + 1 equations of system (7) is the equation system
for aerosols. It is also given by taking the first N + 1
size moments of Eq. (3) without the diffusion term and
the gas acceleration term. Indeed in the case of aerosol,
u = g, and so there is no equation on the momentum
but diffusion of the droplet number density and higher
order moments which we do not take into account here.
In parallel to system (7), the gas phase is solution of the
Euler system for compressible gas:
at (p,)+ V(pgUg)
dt(PgU9)+V(pgU i9 9 Ug)
9t(pE,)+ V(p9E9u,)
~,P9
V, (P iU9)
where pg, P9, E, are respectively the gas density,
pressure, and total energy. The gas is assumed to be a
perfect gas and so follows the corresponding equation of
state.
System (7) contains an unclosed term, n(t, x, S 0).
In practical, Kn(t, x, 0) represents the disapearrance flux
of droplets through evaporation. The challenge here is
to find continuous values of n(t, x,.) from the data of
its first moments. This problem is the finite Hausdorff
problem (Dette et al. (1997)) for the set [0,1]. Such a
reconstruction is possible if and only if the numerical
scheme allows to stay in the moment space. This requires
to know the key properties of the moment space detailed
in the next paragraph.
2.2 Property of the moment space, and
closure of the evaporation term
In this part, we focus only on the size phase space.
Therefore, we consider a homogeneous and stationary
flow. Another way to say that AM belongs to the (size)
moment space for the compact interval [0, 1] is to be sure
that at least one (size) NDF exists on [0, 1] the moment
of order k of which is equal to mrn, k ranging from 0 to
N. We call it the realizability condition. Necessary and
sufficient conditions for existence of a (non unique) n(S)
are non negative Hankel determinants from the theory of
canonical moments (Dette et al. (1997)). Thus, checking
if a N +1 component vector belongs to the moment space
becomes quite tedious when N > 2 because of the com
plex geometry of the moment space. Indeed, the Hankel
determinant for mk depends on all the moments of lower
order. Geometrically, the interval of admissibility for mn
recursively depends on all its predecessors. Nevertheless,
it is quite easy to see that the moment space is convex, but
that it is not a vectorial space.
However, from the Hankel determinants, it is possible
to derive quantities called canonical moments, linked
with a onetoone mapping to the moment space. The set
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
of the canonical moments, called the canonical moment
space, has the very convenient property of fully lying
in the set [0, 1]h, and not in a subset of it. Therefore, in
order to check if the moment space is preserved, it is
more practical to work with the canonical moments.
Let us introduce the set of all probability measures on
[0, 1] whose moments of order up to k1 are cj, P(c_ 1),
where Ck 1 (c, ...,
mo
normalized moments. Dette and Studden have given the
definition of the canonical moments of order k(k > 1):
P C Ck (k 1) (9)
C k k <' (c)
1()1) (C1)
C((Ck 1) max c (/t) and ck (k 1)
P P(ck 1)
min Cek(p) are respectively the upper and lower
P P(ck1)
boundary of the admissible interval for the vector
c = (1, Cl,..., ck)t to be in the moment space.
The canonical moments have two major intrinsic prop
erties which make them attractive to work with. First,
according to the definition (9), and as it as already
been stated, each canonical moment independently
lies in the interval [0, 1]. If, for a certain value of k,
pk 0 or 1, then Vj > k, the canonical moments
are not defined, and the distribution is a sum of Dirac
distributions. Secondly, the canonical moments remain
invariant under linear transformation of the distribution,
i.e Vk > 1,pk(f) p= (fs,_s....), where fs_ ,s,,.
denotes the distribution induced by the linear transfor
mation S = S,in + (Sm,,, S,min,) of [0, 1] onto
[Sn,;i, S,,,] (refer to (Dette et al. (1997)) for proof).
That is the reason why we can work on the size interval
[0, 1] without loss of generality.
The closure of the evaporation term is done using
an Entropy Maximization (Mead et al. (1984)). Con
sequently, the value of the NDF appearing in the first
equation of system (7) is not exactly n, but a recon
structed distribution, we will note n ME. This closure is
proven to preserve the moment space during the whole
dynamical process. Moreover (Massot et al. (2010))
have designed an adapted solver based on a kinetic
scheme because it provides a dedicated tools which
relies on an integral formulation of the fluxes making use
of the underlying kinetic equation Eq. (2). Moreover,
(Massot et al. (2010)) shows that combining this kinetic
scheme with a DQMOM approach allows to guarantee
the preservation of the moment space and to solve an old
issue in the literature. Once evaporation is dealt with, we
have to preserve the moment space through convection.
In addition to the two previous properties, another prop
erty given from the convective system extracted from sys
tem (7) will be very useful when designing the numerical
scheme. The canonical moments are transported quanti
ties, which means that they verify the transport equation:
Otpk + uOpk =1 0. The demonstration can be found in
(Kah et al. (2010)). These properties of the canonical mo
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
ments have enabled to design a high order in space and
time numerical scheme for the transport of the moments,
meeting the realizability condition per se, solving the is
sue discussed in (McGraw (2007), Wright et al. (2007)).
2.3 Numerical validation of the Eulerian
high order moment method
The aim of this section is to show that numerical
tools have been successfully designed to solve system (7)
on a Eulerian cartesian mesh. This paragraph gives an
outline of the validations in this context, and we refer to
(Kah et al. (2010)) for more details.
Because of the transport in physical space and the
transport in phase space due to evaporation and drag have
different structures, we use a Strang splitting algorithm
(de Chaisemartin (2007)). We first solve for At/2 the
transport in phase space, then for At the transport in
physical space, and then for At/2 the transport in phase
space. The transport in physical space obeys a system of
weakly hyperbolic conservation laws and relies on kinetic
finite volume schemes as introduced in (Bouchut et al.
(2003)) in order to solve the pressureless gas dynamics
equation. The properties of the canonical moments
provides a second order formulation of the fluxes which
preserves the moment space. For a twodimensional
space, we further use a dimensional Strang splitting of
the 1D scheme previously described in (de ( Il.icii.iiilii
(2007)). For the transport in phase space through drag
the model equations reduce to a system of ODE's,
which can be stiff, for each point of the domain. Each
system is solved using an implicit RungeKutta Radau
IIA method of order 5 with adaptive time steps. The
transport in the phase space through evaporation is solved
using on the principles highlighted is the previous section.
In order to assess the Eulerian methods we focus on
a 2D Cartesian free jet. A polydisperse spray is in
jected in the jet core with either a lognormal size NDF,
corresponding to the beginning of a typical experimen
tal distribution (Laurent et al. (2004)). The simulations
are conducted with an academic solver, coupling the AS
PHODELE solver, developed at CORIA by Julien Reveil
lon and collaborators (Reveillon (2007)), with the multi
fluid solver MUSES3D developed by St6phane de Chaise
martin and Lucie Fr6ret at EM2C Laboratory (de Chaise
martin (2009)). The ASPHODELE solver couples a Eu
lerian description of the gas phase with a Lagrangian de
scription of the spray. One of the key feature of this tool
is to allow, in the framework of oneway coupling, the si
multaneous computation of the gaseous phase as well as
both spray descriptions within the same code run.
As far as the gas phase is concerned, we do not solve
system (8) but we use a 2D Cartesian low Mach number
NavierStokes compressible solver. The gas jet is com
puted on a 400 x 200 uniformly spaced grid. To destabi
lize the jet, we inject turbulence using the Klein method
with 10% fluctuations (Klein et al. (2003)). The Reynolds
number based on Uo, vo and Lo is 1, 000. The gas vortic
Figure 1: Gas vorticity field of the freejet configuration
at final time (t 20), on a 400 x 200 grid
0.8
0.6
0.4
0.2
8 10 12
Figure 2: Polydisperse evaporating spray of the freejet
configuration at final time (t 20) for droplets
with Stokes 1.88. (Top) Results for the num
ber density mo. (Bottom) Results for the third
order moment ms. The computation is carried
out on a 400 x 200 grid
ity is presented in Fig. (1).
Figure (2) shows the results for the 0th order moment
mo and the 3rd order moment m3 of the injected evapo
rating spray. The corresponding nondimensional evapo
ration coefficient is K = 0.14. The mean Stokes number
of the droplets is St = 1.88. These results have been val
idated by comparison with results given by the multifluid
model (Kah et al. (2010)).
These schemes have been designed on fixed eulerian
grid in the context of finite volume schemes. However,
the problem we focuses on involves a moving mesh to
take into account the displacement of the piston. Our task
is then to adapt the proposed moment method in order to
describe polydispersivity in a moving grid formalism.
2.4 Derivation of the system in ALE
formalism
2.4.1 Description of the motion
Because of the shortcomings of purely Lagrangian
and purely Eulerian descriptions, a technique has been
developed that succeeds to a certain extent in combining
the best features of both the Lagrangian and the Eulerian
approaches. Such a technique is known as the Arbitrary
LagrangianEulerian (ALE) description (Donea (2004)).
In the ALE description, the nodes of the computational
mesh may be moved with the continuum in normal La
grangian fashion, or be held fixed in the Eulerian manner,
or be moved in some arbitrarily specified way to give a
continuous rezoning capability. Because of this freedom
in moving the computational mesh offered by the ALE
description, greater distortions of the continuum can be
handled than would be allowed by a purely Lagrangian
method, with more resolution than is afforded by a purely
Eulerian approach.
The Lagrangian viewpoint consists in following the
material particles of the continuum in their motion. To
this end, one introduces a computational grid which fol
lows the continuum in its motion, the grid nodes being
permanently connected to the same material points. The
material coordinates X, allow us to identify the material
configuration Rx. The motion of the material points re
lates the material coordinates X, to the spatial ones x. It
is defined by an application p such that:
9: Rx x [to,ttfinal[ R. x [to,jtfinal[
(10)
(X, T) O(X, T) = (x, t)
which allows to link X and x by the law of motion and to
deduce the material velocity, namely:
x = x(X, T), t = T, v(X, T) = dXlx (11)
which explicitly states the particular nature of p. v can
refer either to the liquid or gaseous velocity. The spatial
coordinates x, depend both on the material particle X,
and time t, and, second, physical time is measured by the
same variable t in both material and spatial domains.
However, in the ALE description of motion, neither the
material configuration Rx nor the spatial one R, is taken
as a reference. Thus, a third domain is needed : the refer
ential configuration Rx, where the reference coordinates,
X are introduced to identify the grid points. The referen
tial domain Rx is mapped into the spatial domain by b.
The mapping D from the referential domain to the spa
tial domain, which can be understood as the motion of the
grid points in the spatial domain, is represented by:
D : Rx x [to,tfinal[j R. x [to,tfinal[
(x,t) (X,,t) (x,t)
and its gradient is defined by
O(x,t ) ( 0oT
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where ux is the mesh velocity. We define also the quan
tity J(X, t) det(Vxx), which is the dilatation rate of a
volume d3X with time. In order to close the final system,
we need the time evolution of J, given by the transport
theorem:
it J J=J ux (14)
2.4.2 Equations in a moving frame
From the Eulerian system (7) and (8) for the moments and
the gas, we can now define the final system of equations
in the ALE formalism. The details of the algebra can be
found in (Despr6s (2005)), (Hirt (1974)).
Let us consider the term a, where a can refer either to
the gas density, p, the gas momentum, pu,, the gas total
energy,, mnk,k 1..N, and moup.
The equations are written over the control volume, written
V(((X, t), t), whose envelope velocity is ux.
Applying the dynamic theorem leads to :
d
dt V(4'(xt),
t) a 3 av(+(x,t),t)
r.n d2x + /
Jv7(,(xf)t)
a (v u) d2x
Sd3x
(15)
where a is the generalized stress tensor and S the
volumic production rate, for a.
After doing the change of variable (x, t) = D(X, t)
leading to the change of elementary volume d3 = Jd3X,
we obtain:
v() (OtJa x+JVa(vux)JVo'JS) d3x 0
(16)
where V(X, 0) is the control volume at time t = 0.
The final model writes:
t(J)lx J%(ux)
dt(Jp)lx+ JV (dpw,)
dt(Jpug,)l+ JV.(,. w I.)+
dt(JpEg)l+ JV.(pE wi.)+
at(Jmo)lx+ J%(mowp)
dt(Jinm)I+ JVt.(mnkwp)+J
t (Ji~miUp)x+J% (mi7UpWp)
0,
=0,
JV (P)=0,
JV1(Pu,)=0,
J KnME O,
K(k l)mk=0,
mo (u u)
mi 0
where the coefficient k ranges in [1, N], and where wg
up ux and wp = ug uX.
(12) 3 Numerical scheme
3.1 Dimensional and operator splitting
(13) In this section, we restrict ourselves to cartersian
meshes to show some results for the design of a second
x x(X,t)= datlx
1 )
order scheme in space and time in the ALE formalism
for the transport of moments. We can first work on a
onedimensional scheme, and extent the dimension of the
domain easily through a dimensional splitting algorithm.
Indeed, designing a high order scheme in space with an
unstructured grid brings additional difficulties we have
not tackled yet. For now on, all the quantities are defined
on the onedimensional grid.
The system writes now more easily in the referential
configuration:
at(Jp) l
ot(Jp)l+
Oa (a )
ax(pw9)
(J' ). Ox(,' .. . )+
at(JpE )l + (,. .. .)+
at(Jmo)l+ a(mowp)
9t(Jmk)l + 9X(mkwp)+J (k
O (P)
9 (P., ):
JKniME
 l)mT 1
9t(Jm1up) +x(mlUpWp)J m
X X? 1 0
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
step corresponds to the resolution of two psystems, one
for each the gaseous and the liquid phase. As the mesh
nodes moves at the speed of each phase, to virtual meshes
coexist at the end of this step.
In the last step (phase C), also called the rezoning step,
the convective terms are solved in a Eulerian formalism.
The two virtual meshes, output of the Lagrange phase,
are moved to their common final location corresponding
to a displacement with the velocity ua. The quantities are
expressed in this final mesh, which is the reason why this
step is also called the projection step.
0, We divide the domain, [0, Z] into N cells
0 [xi1/2,xi+1/2], each cell having its own size Axi.
The inner cells are numbered from 1 to N. We also
0, define two ghost cells labeled 0 and N + 1. The ther
0, modynamical quantities are defined on the cells, so for
0, example pi is defines as: pi J1/2 p(x) dx.
The velocities are defined on the nodes, we write them:
, 1/+2 +2p,i l/2, up,i+1l/2
The formulation (17) turns out to separate acoustic
waves of the gas from kinematic waves. Therefore, the
basic idea of ALE approaches is to perform a splitting
of (17), within a timestep At. This approach is fully
accurate when the kinematic waves speed is negligible
compared with the accoustic waves speed. This hypoth
esis is not obvious in the context of injection in Diesel
engines, where the highest speeds of injection can now
reach more than 600ms1. Nevertheless the fact that this
method enables to take into account a mesh movement
without the drawbacks of the purely Lagrangian methods
is a prevailing argument.
We use a Lie fractional timestep algorithm decom
posed in three steps. The first step (phase A) corresponds
to the resolution of the source terms of the system, which
in our case are the evaporation and the drag term.
The second step (phase B) is called the Lagrange
step, taking into account only acoustic effects due to the
gaseous pressure, which is the only source of accoustic
effects as the spray is treated as a pressureless gas. The
system solved in during that step reads:
dt( Jp) Ix
at(JgP) x ,
a9.. .)=0,
=0,
at(J,',, )lx + 9X(P)= ,
O(Jp)lxp +x(p)'O,
at(Jpinm)x, =0, k = 0..N
at(Jpml p)x =0,
Jg denotes the dilatation due to the gaseous velocity field
and Jp the dilatation due to the liquid velocity field. This
Moreover, with Ug = (p,,.. pE,), p
("'' 1..N, moUp), and the subscript m either referring
to the gas (g) or the particles (p), we denote by (J,', U ),
(JAU), (J ,U(), (JCU ) (Jg 1,U 1) the
values of (Jm, U,) respectively at time t nAt, at
the end of phase A, at the end of phase B, and at the
end of phase C which are also the updated values at time
t (n + 1)At.
3.2 Resolution of the Lagrange step
Because of the different nature of the waves intro
duced in the equation systems for the gas and the particles
in system (17), the Lagrange step is solved with dedicated
solvers for each phase.
3.2.1 Resolution for the gas with an implicit
scheme
For the gas, the Lagrangian step corresponds to the reso
lution of a psystem, with the treatment of the accoustic
waves. In order to avoid small time steps leading to high
computational costs, an implict time integration is an es
sential requirement.
From the mass conservation equation, we first use that
the mass in a material volume is conserved: (JYp)B
(Jgp)A. Extracting this relation in the three other equa
tions leads to:
1
Otr .. =0, where
P
+ ., +aP =0
P
p
1_
atE9+9 Pu=0
p
Before discretizing these equations, let us remark that the
total energy E, is not well defined on the grid, as the in
ternal energy is defined on a cell, and the kinetic energy
on a node. We will then work with the gas internal en
ergy eg, the equation of which is deduced from a combi
naison of the equations of system (20). We assume that
all the quantities are C1, and so that no shockwave oc
curs. This assumption is relevant as we consider the spray
downstream the injector, with a substantially lower veloc
ity in comparison to its injection velocity. Besides, we
are not interested in the description of the shock waves in
the gas. The numerical scheme used to solve system (20)
is a simplified version of the standard acoustic Godunov
scheme (Godunov (1959)), which can be interpreted as a
relaxation scheme (Coquel et al. (2010)). The discretized
equations read:
At(U* u*
S TA g,i 1/2 g,
TA iA
B A At(P+1 P)
1 2 12 A A
PiB A
CB A +P (T _ )
"3,t "3,t PBt 7'z j4l)
0 (24)
0 (25)
0 (26)
where
P PB
U9,i1/2 U9,it1/2
P1* P
pB pB
2a
UB
9,i 1/22
a
U B
9,i1/2
The parameter a is homogeneous to pc, the Lagrangian
speed of sound in the gas. However, it can shown that
a must be subjected to the subcharacteristic condition
a > max(pc)i.
The resulting nonlinear matrix equation is inverted by
the Newton method.
3.2.2 Explicit resolution for the liquid
As for particles, there is no pressure term in the mo
mentum equation. This implies that the dilatation coef
ficient Jp is modified only by the divergence of the par
ticle velocity field which is unchanged during this phase.
We can thus use an explicit method for the resolution of
the system concerning the moments. If we write T7,
1/mk,i, k 1...N, one gets:
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
3.3 Resolution of the convective terms
The outcome of the Lagrange step is now the input
data for the projection step. The latter amounts to solving
dtJmlx+Oxwm =0,
atJmU'lx +OxU. =0, (31)
Combining the equations of system (31), we obtain the
componentwise advection equation, one gets:
W+m
OtUm + OxUJ = 0, wM = uM ux
'"ri,
At the discrete level, the transport velocity for the gas
should be taken equal to ... ua, where .* is the out
put of phase B. This ensures that after phase A and B,
the mesh has moved exactly at velocity ux The projection
step is merely a remap of the variables contained in Um,
at the velocity relative to the referential domain.
Let y refer to a celldefined quantity. As y is a solu
tion of a conservative equation, a finite volume discretiza
tion is appropriate for its discretization. We denote by
Ax JmAx' the cell volume at the end of the La
grange step, and by Axz+1 the common final cell volume
for both phase. Integrating the second equation of system
(31) over the ith cell, we get the discretized equation:
'*A +1 p BAXZB At((yBw,)(nX+l/2)
( Wr 1/(i 2))
The values ( BWrm)(xi+l/2) and ( BWm)(zi 12) are
computed using a first order upwind scheme.
The resolution of the momentum equations are more
tedious, as the first thing is to define the quantity ,.. .
and moup. This is not a priori obvious as p, mo, and
p are not defined on the same domain. We define
the momentum on the nodes. This means that we solve
the momentum on an off centered mesh or a dual mesh
with respect to the regular one. The dual cell centers are
the nodes of the regular cells, and the dual cell interfaces
are the regular cell centers. Using a weighted average by
p and mo on each side of the node. For the gas, phase,
this writes:
nB B pf Aa /2 + p B+iAzx+1/2
(P),+1/2 B +1/2 Ax/2 + AxZ,/2 (34)
In the same time, we define the dual interfacial veloci
ties as the average of the velocities of each regular node:
Wm,i1/2 + Wm,i+1/2
Wm,i =2
The discretized equation write, for the gas:
(, .)n 2(A+ 1/2 + A /2/2)
At(( 1 ..... )( +i1) ( 3 ^ .... )(3i))
B A At
tk,i Tk,i( (l X Up,i+l/2 Up,i 1/2))
For Tk i to remain positive, we impose At <
Up,i 1/2 Up,i1/2
Ax zA
As before we use a first order upwind scheme to compute
the fluxes.
We may also want to have a scheme with a higher or
der in space and time. We will focus in this paper on the
increase of the spatial order. The approach in the context
of offcentered mesh is different from the approaches on
fixed Eulerian grids as the velocities as defined on nodes
and therefore cannot be reconstructed the same way as the
celldefined quantities. We present a second order scheme
in space.
For the gas phase, a classical linear reconstruction with a
minmod limiter is used for the celldefined quantities. The
momentum, expressed on each node from the velocity and
the density, are reconstructed on the dual mesh, with the
same type of reconstruction as for the celldefined quan
tites. This scheme guarantees the sine qua non condition:
positivity of the density and internal energy. But when it
comes to the particles, the size moments must not only
stay positive, but also verify the realizability condition, as
explained in (2.2). However the last property is not en
forced by a scheme where the moments are independently
reconstructed. Therefore, we present in the next section a
scheme which enables to verify the realizabiliy condition.
3.4 Second order in space and time and
preservation of the moment space
In (Kah et al. (2010)), a scheme has been proposed,
relying on the reconstruction of canonical moments,
in order to preserve the realisability condition directly
during the convection process. Let us recall that this issue
has been already tackled before in (Wright et al. (2007))
and (McGraw (2007)) for example, but these authors
have to project the moments back into the moment space
after the convection step.
Defining the discrete values of the moments over the
mesh at the end of the Lagrangian step
B 1 +1/2
mk,i AxB j
L/2
0, ..N (36)
The discretization of the convection part in the equa
tions on moments in system (17) leads to:
M 1A 1 MfAf At(Fi 1/2Fi1/2) (37)
where Fi+1/2, the moment flux, writes:
F 1/2 j j V
S(X.i1/2) dvdt
/
To design a second order scheme in space, one chooses
to reconstruct the moments linearly in the ith cell. Two
difficulties are encountered. The first one is that, given
the complexity of the moment space, we are not sure
that, the reconstructed vector at the point X.i+1/2 is still a
moment vector. Moreover, assuming that we can be sure
of that, there is no guarantee that, after adding the fluxes,
the updated moment is still a moment vector. Indeed, the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
moment space is not a vectorial space.
Both difficulties are solved when considering the ki
netic level of description and looking at the macroscopic
equations as derived from the kinetic equation (2). The
NDF f is the solution of the transport equation:
dt(f)+vd(f) =0, in T. x R, xRsx]t,,, +i[ (39)
so that f(t, x, v, S) = (0, x vt, v, S).
Thus, expressing the moments according to f leads to the
following expression for F,+ 1/2
At 1
0 R 0
1
s
vf(t, Xi+1/2S, v, ) dSdvdt
(40)
We can compute the integrals over R+ and ]R separately,
and then the numerical flux can be written in the flux
vector splitting form: F,+1/2 = F+1/2 + Fi1/2, where
the distinction is made between positive and negative
velocities.
Using the fact that f is transported at the velocity v
enables to make a change of variable in the expression
(40), so that F+,2 can be expressed as a space integral
of the moments over the ith cell. Provided that the
reconstruction of f ensures that the scheme is conser
vative and that the realisability condition is satisfied
VX E ['i1/2,Xi+1/2], the updated vector M '1 be
longs to the moment space. In fact, in that case, when
considering positive velocities with no loss of gener
a t
ality, the quantity Mf F 1/2 is a moment vector.
To achieve a reconstruction meeting the two above re
quirements, we reconstruct the canonical moments, in
]Xj1/2, j+ 1/2 [ taking benefit from the fact that they are
transported and they lie in [0, 1]:
Smo(X)
{l(ar)
< 9 ( '
".. i + Dmo ( Xj )
Plj + Dp1, (x xi)
pNj + DpNi (x X j)
PNJD^PNjX x)
Where Dpi' and DpN,i are the slopes and the quantities
with bars, called the modified averages, are different from
their real mean value in order for the scheme to be conser
vative.
Indeed, given T a transported quantity and C the corre
sponding conserved quantity, the scheme is conservative
if:
SC(x) = mo(x)T(x)
Ci = imo,i =  I L2 mo(x)T(x) dx
ax t/2
The recursive dependance of higher order canonical
moments makes their modified average more difficult
to express. High order polynomial must be integrated,
(up to order 4 for p2, and 6 for 2p). The quantity "
writes pki = ak,i + bk,iDpk, The calculation of the
coefficients .. and b2j, bj is achieved using Maple
(Maplesoft, a division of Waterloo Maple, Inc 2007),
and implemented in the code written in Fortran. Their
expression is quite heavy, but, as it is just an algebraic
relation, the corresponding CPU cost is low.
The slopes for the canonical moments are determined
using limiters in order to satisfy maximum principles:
rki < Pk(x) < Rki, x E (xi, 1/2, Xi+ 1/2),
where rij min(pk,i 1,pk,i,Pk,i+i) and Rij
max(pk,i 1,Pk,i,Pk,i+l)
We must have:
rki < ak,i + bk,iDp,k +
rci < ah,i + bc,iDP,i 
AxB
Dp,,,i < Rkci
AxB
 D,,,i < Rhi
The slopes must then verify:
SDk < in Rki ak,i ak,i mi
xbh,i + Ax/2' Ax /2 b, i
DPk ki a,i 
Kbh, + Axf/2' Axf /2 b,
In practice, we use the following slope limiter to satisfy
all the conditions:
Dp,, = (s.gn(ph,i+ ,i) + sgn((Pk,i Pk,i 1))x
i (IPk,i+1 ak,i lak,i Pk,i+l 1
(2( + bk,i) 2( + bk,i)
The same limiter is used for the slope on the droplet
number mo. In order to ensure the nonnegativity of mo,
one adds the following condition:Dm,,o Ax. /2 < mo,i.
In practice, we set N 3, so we work with the four
first size moments. The final expression of the size flux
containing the canonical moments is:
1 7i 1/2 **L(x)
F1/2 ,, pi(( p)p2 + p))() dx
1/2 J '.P((1 )P2) 1 P2)P3
Xi+1/2 \+(pi +p(1 p))2))(X) /
where x +1/2 X+1/2 Um,i+1/2At, where am,i+1/2
either represents the gas velocity in the case of aerosols,
or the particle velocity in the case of sprays. This scheme
is second order in space and time for the moments.
3.5 Validations
Some computations have been performed with a code
dedicated to validate the use of moment methods in the
ALE formalism in one dimension. A first result shows
the possibility to transport and evaporate and aerosol.
The interest for a second order in space scheme is then
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
highlighted when the grid moves. Finally a test case
involving an evaporating spray validates our approach for
sprays as well.
For the first test case, the particles are considered as
tracers for the gas phase. We remark that, as for the
particles, the quantities describing the gas are dimen
sionless. One writes P P'P, T T'T, p p'po,
where Po, To, po are respectively a referential pressure,
temperature and density. The initial conditions for the
gas are the Riemann problem's conditions: the left
state is set as P' 3, p 1, and for the right state
P' 1, p' 0.125. The NDF for the polydisperse
spray is set as n(x, S) p', S E [0, 1] so that mo p'.
Physically, mo and p have no reason to be equal. mo
is the droplet number density and p is the gas density.
mo and p are initialized with the same value only for
a numerical reason. That way, the profile of mo can be
easily compared to the profile of p'. The particles are
evaporated with the coefficient K = 2. K is such as
AS KAt, so that all the droplets must be evaporated
at t = AS/K. Neumann boundary conditions are set on
the space interval [0,1], and the CFL is set as 1, on the
grid with 200 cells.
Figure (3) shows the results for the gas density, mo,
and the analytical solution for this case at t 0.1.
Focusing on the resul for mo, it is everywhere equal to
0.8p'. This is coherent with the initial NDF, and the
evaporation rate K 2, implying that _'' r. of the total
droplet number must have been evaporated. We have
only displayed the results for mo and mi for the sake
of legibility, but m2 and ms are also solved. The first
conclusion to draw is that it is possible to transport an
evaporating aerosol in the context of the ALE formalism.
But this is not a challenging issue as the scheme used for
that case is first order in space.
Figure (4) displays the results for mo of the same
test case, but comparing a first order scheme with the
second order scheme explained in (3.4). The accuracy of
these schemes is compared on a fixed grid but also in a
moving grid. The displacement of the grid follows the
27r
law : x(t) 0.2cos(0 t). The results focused on the
interface region. When no grid movement is involved, the
second order scheme is already more accurate than the
first order. This conclusion validates the implementation
of our high order moment method. But the interest of
a second order scheme becomes obvious when the grid
moves. In our case, the high grid velocity leads to small
CFL numbers for the fluids. Therefore the profile of mo is
much more diffused with the first order scheme than with
the second order scheme.
The particles are now considered to have their own dy
namics. The NDF is more complex than before, in order
to show that the evaporation solver does not presume any
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
type of NDF:
SA({x) sin(7s) + (1
n[ 0
Figure 3: Solution of the Riemann problem at t 0.1 for
the gas and the aerosol, with the first order in
space scheme. The analytical solution is repre
sented by the solid black line, the dashed ma
roon line represents the numerical solution for
the gas, the blue curve with stars markers and
the red curve with up triangle markers respec
tively stands for mo and mi. The computation
is carried out on a 200 cell grid
Figure 4: Focus on the interfacial area in Fig.(3) and
comparison of the first and second order so
lution for mo, with and without mesh move
ment. The solid blue line corresponds to the
second order solution without mesh movement,
the dashed red line to the second order solu
tion with mesh movement, the dotted violet line
to the first order solution without mesh move
ment, and the green line with square mark
ers to the first order solution with mesh move
ment.The computation is carried out on a 200
cell grid
A(z))exp(s) ifx < 0.5
otherwise
where A(x) = 4(0.5 X)2.
The particles velocity is initiated by:
a(x) f 0.5 if x < 0.25
2 if x > 0.25
We set periodic boundary conditions on the space in
terval [0, 1]. The analytical solution is the translation of
the two parts of the density profile corresponding to each
value of the velocity. Figure (5) displays the four ana
lytical size moments, and the size moments given by the
calculation at time t 0.2 on the left, and at the time
t 0.45 on the right. One can notice first that the mo
ment space is preserved, at any time.
At t 0.2, Fig.(5), the initial distribution breaks
into two parts. Vacuum is created at the initial velocity
discontinuity. The numerical solutions are represented
by solid lines with decreasing ordering in terms of value,
meaning that the highest curve stands for the 0th order
moment, and the lowest curve stands for the 3rd order
moment. The analytical solutions for the corresponding
moments are represented by marquers. It can be seen
that the numerical solution very accurately matches the
analytical one.
At t = 0.45, Fig.(6), because of periodic boundary
conditions, the fastest portion catches up the slower one.
In this figure, the numerical solutions are represented by
dashed line, the numercial solution by solid lines, and
only mo and mi are represented for the sake of legibilty.
The two types of solution are very different. Indeed,
the models corresponding to each type of solution are
different. The numerical solution corresponds to the
resolution of system (17), written in the pressureless gas
formalism. In this context, as the pressure is considered to
be null, nothing prevents the particles from accumulating.
A 6shock is therefore engendered when two droplet
clouds with different velocity meet (Bouchut (1994)). On
the other hand, the analytical solution of this problem in
the infinite Knudsen limit for the particles corresponds to
the solution at the kinetic level, where the droplets cross.
One can notice though that in the region where only one
droplet cloud is considered (in the interval [0.4, 0.5]), the
numerical and analytical solution match. However our
purpose is to show that system (17) can be implemented
in the ALE formalism. With respect to this objective,
Fig.(6) validates our numerical approach for sprays as
the dynamic of the shock is well captured. Nevertheless
simulating jet crossing is an issue in Eulerian models and
methods and has been recently resolved in the literature.
We refer to (de (I.iiciii.iiiiii (2009)) and references
therein for details on this matter.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 5: Evolution of a spray in a discontinuous velocity
field calculated with the moment method, com
pared to the analytical solution of the problem,
at time t = 0.2. The plain curves represent,
with decreasing ordering in terms of value, the
four first moments from the 0th order mo (blue
curve with circles) to the 3rd m3 (violet curve
with diamonds). The analytical solution is rep
resented by marquers. The computation is car
ried out on a 200 grid.
4 Polydispersivity in a CFD code
4.1 Introduction of polydispersivity in
IFPC3D
An Unstructured Parallel Solver for Reactive Com
pressible Gas Flow with Spray is used to integrate the
previous developments. It is a hexahedral unstructured
parallel solver dedicated to multiphysics calculation
being developed to compute internal combustion engines.
Original algorithms and models such as the conditional
temporal interpolation methodology for moving grids,
the remapping algorithm for transferring quantities on
different meshes during the computation enable IFPC3D
to deal with complex moving geometries with large
volume deformation induced by all moving geometrical
parts (intake/exhaust valve, piston). The Van Leer and
Superbee slop limiters are used for advective fluxes
and the wall law for the heat transfer model. Physical
models developed at IFP for combustion, for ignition
and for spray modelling enable the simulation of a
large variety of innovative engine configurations from
nonconventional Diesel engines using for instance HCCI
combustion mode, to direct injection hydrogen internal
combustion engines. Large superscalar computers up to
1 000 processors are being widely used and IFPC3D has
been optimized for running on these Cluster machines.
IFPC3D is parallelized using the Message Passing
Interface(MPI) library to distribute calculation over a
large number of processors.
The moment method is being implemented according to
Figure 6: Evolution of a spray in a discontinuous veloc
ity field calculated with the moment method,
compared to the analytical solution of the prob
lem, at time t = 0.45. For legibility, only
mo (blue curve) and mi (red curve) are rep
resented. The numerical solution is represented
by dashed lines, the analytical solution by solid
line. (Top) Solution in the space interval [0, 1].
(Bottom) Focus on the interest area from the top
figure.The computation is carried out on a 200
grid
the structure of the code, to enable the code to describe a
polydisperse spray with mesh movement.
4.2 Testcases
The implementation of the moment method is vali
dated on two academic test cases. The first test case con
siders an aerosol in a homogeneous high pressure cell with
a moving piston. The second test case considers the Rie
mann problem for a spray.
The objective of the first test case is to ensure that the
implemented model is stable with mesh movement. The
evolution of homogeneous fields of liquid and gas in a
closed high pressure cell are considered. The bottom
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 7: Considered piston movement during the com
putation. The computation starts at cad
180 and ends at cad = 180.
boundary of this cell corresponds to a moving piston
being at the bottom dead center. The gas is taken as air.
and the particles are initially stationary. No ignition
occurs, and no thermic effect is considered. Also, no
special treatment of the boundary is considered. The
computation ends after a revolution of the crank, the
crank angle degree (cad) ranging in [180, 180]. The
high pressure cell and the movement of the piston are
described in Fig. (7).
The size distribution is taken constant like in the
onedimensional tests with aerosol. All the results are
displayed for the number density mo and the surface
density mi with a 1200(30 x 40) cell grid.
Figures (8)(top left) displays the initial conditions for
the moments, and the other charts of Fig. (8) show
the results for cad 100, 30, 140. In the different
graphs, the distance where the value of the moments
is null is the distance traveled by the piston. One first
notice than the flow stays homogenous during the whole
computation. This is a consequence of the Ilpllhsis,
made in the ALE formalism assuming that the acoustic
time scale is negligible in comparison to the convective
time scale. This hi\plslli, is verified in our case. The
value of the speed of sound exceeds 300ms1. In the
same time, with a rating of 1200trmn 1, and stroke
of 10cm, the piston velocity and so the velocity of the
dragged fluid is much smaller than the speed of sound in
the gas. Mass conservation impose that the gas density
and so the particle number increase as the piston heads
to the top dead center, because the volume of the high
pressure cell decreases (Fig.(8)(top right),(8)(bottom
left)) The moments recover their initial values at the end
of the computation.
The second test case involves a spray, solution of the
Riemann problem for the following initial conditions
{ (mo, mi, mn2, m3) (1, 0.5, 0.333, 0.25), Vx E [0, 1]
v 1000, x E [0.4,0.6]
v 500 otherwise
Figure 8: Solution for aerosol particles in the case of pis
ton movement, where mo and mi are repre
sented. (top left) Initial condition for the homo
geneous chamber. (top right) Solution at cad
140. (bottom left) Solution at cad 30.
(bottom right) Solution for cad 140. The
computation is carried out on a 30 x 40 grid
Periodic boundary conditions are imposed on a 100 grid
and the final time t 104 is chosen such as a particle at
constant velocity 1000 comes back to its initial position.
Figure (9) displays the results of this test case. In Fig
(9)left a 6shock occurs at the location of the velocity
discontinuity corresponding to the imping of particles.
At the same time, a rarefaction wave occurs on the other
velocity discontinuity. At time t 103, the dynamics
of the shock and the strucutre of the rarefaction wave are
established and one can verify that this corresponds to the
solution of this Riemann problem.
These tests validate the developments done so far in
C3D. This is a first step towards a full implementation
and multi dimensional aerosol and spray computations.
5 Conclusion
This paper shows the feasability to design a stable
scheme for high order moment methods in an ALE
formulation. The main difficulty when working with
a moment vector is to preserve the moment space, or
equivalently, to ensure that an NDF can be reconstructed
from the updated moment vector. In order to achieve
that, we have adapted, in the context of convection, the
kinetic scheme designed in (Kah et al. (2010)) during
the rezoning phase, despite the stability problems of
these kinds of method. A second order scheme in space
and time has been designed for aerosols and sprays.
Combined with the use of high order moment method
for evaporation developed in (Massot et al. (2010)), this
has lead to the capability of transporting an evaporating
spray with a second order scheme in a moving mesh,
which takes up the second challenge highlighted in the
3 10
2.5 8 
6 
1.5
0.5 
0 0
0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1
X X
1000 1000
900 900
800 800
700 700
600 600
.... I500
SI 0 0.2 0.4 0.6 0.8 1
x x
Figure 9: Solution for spray particles in the Riemann
problem. (Top) Solution for mo and mi. (Bot
tom) Solution for the velocity field. (left) solu
tion at time t 10 4. (Right) Solution at time
t 10 The computation is carried out on a
100 grid
introduction and is the main point of this paper.
The third and second main task of our contribution is
being addressed with the implementation of these high
order moment method in a unstructured CFD code. The
validation cases shown in this paper are the first step
towards the objective of the multidimensional simulation
of a polydisperse jet in combustion chambers.
Acknowledgements
This research was supported by a Ph.D. grant CIFRE from
IFP and EM2C laboratory, Ecole Central Paris. We also
would like to thank Antony Velghe and Julien Bohbot
(IFP) for their help.
References
Bohbot, J., Gillet N. and Benkenida A., IFPC3D: an Un
structured Parallel Solver for Reactive Compressible Gas
flow with Spray, Oil and Gas Science and Technology,
Volume 24, 2009
Bouchut F., Jin S. and Li, X., Numerical approximations
of pressureless and isothermal gas dynamics, SIAM Jour
nal on Numerical Analysis, Volume 41, number 1, Pages
135139, 2003
Bouchut F., On zero pressure gas dynamics, Advances
in kinetic theory and computing, World Sci. Publishing,
Pages 171190, 1994
Coquel F., Nguyen Q. L., Postel M. and Tran, Q. H.,
Entropysatisfying relaxation method with large time
steps for Euler IBVPS, Mathematics of Computation, S
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
00255718(10)023392, Article electronically published,
2010
de ( Ii.icisNni.iiiii S., Freret L., Kah D., Laurent F., Fox R.
O., Reveillon J. and Massot M., Eulerian models for tur
bulent spray combustion with polydispersity and droplet
crossing : modeling and validation, Proceedings of the
Summer Program 2008, Center for Turbulence Research,
Stanford University, Pages 265276, 2009
de Chaisemartin S., Polydisperse evaporating spray tur
bulent dispersion: Eulerian model and numerical simu
lation, Ph.D. thesis, Ecole Centrale Paris, 2009, avail
able in English on TEL http://hal.archivesouvertes.fr/hal
00443982/en/
de Chaisemartin S., Laurent, F., Massot, M. and
Reveillon, J., Evaluation of Eulerian multifluid ver
sus Lagrangian methods for the ejection of polydis
perse evaporating sprays by vortices, European project re
port, TIMECOPAE Project (2007). Available on HAL:
http://hal.archivesouvertes.fr/hal00169721/en/
Chandrasekhar S., Stochastic Problems in Physics and
Astronomy, Reviews of Moder Physics, Volume 15,
number 1, 1943
Cowling T. G. and Chapman S., The Mathematical The
ory of NonUniform Gases, Cambridge University Press,
1970
Despr6 B. and Mazeran C., Lagrangian gas dynamics in
twodimensions and Lagrangian systems, Arch. Ration.
Mech. Anal., Volume 178, number 3, Pages 327372,
2005
Dette H. and Studden W. J., The theory of canonical mo
ments with applications in statistics, probability, and anal
ysis, Wiley Series in Probability and Statistics: Applied
Probability and Statistics, John Wiley & Sons Inc., 1997
Donea J., Huerta A., Ponthot J.Ph. and RodriguezFerran
A., Arbitrary LagrangianEulerian methods, Encyclope
dia of Computational Mechanics, Pages 413437, 2004
Greenberg J.B., Silverman I. and Tambour Y., On the ori
gin of spray sectional conservation equations, Combus
tion and Flame, Volume 93, Pages 9096, 1993
Fox R. O., Laurent F. and Massot M., Numerical simula
tion of spray coalescence in an Eulerian framework: direct
quadrature method of moments and multifluid method,
Journal of Computational Physics, Volume 227, Pages
30583088, 2008
Freret L., de Chaisemartin S., Reveillon J., Laurent F.
and Massot M., Eulerian models and threedimensional
numerical simulation of polydisperse sprays, 7th Interna
tional Conference on Multiphase Flow, ICMF 2010
Godunov S., A difference scheme for numerical compu
tation of discontinuous solution of hydrodynamic equa
tions, Math. Sib., Volume 47, Pages 271306, 1959
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Hirt C. W., Amsden A. A. and Cook J. L., An arbi
trary LagrangianEulerian computing method for all flow
speeds, J. Comput. Phys., Volume 14, Pages 227253,
1974
Kah D., Laurent F., Massot, M. and Jay, S., A high order
moment method simulating evaporation and advection of
a polydisperse spray, to be submitted to Journal of Com
putational Physics, 2010
Kaper H. G. and Ferziger J. H., Mathematical theory of
transport processes in gases, NorthHolland, Amsterdam,
1972
Klein M., Sadiki A. and Janicka, J., A digital filter based
generation of inflow data for spatially developing direct
numerical or large eddy simulations, J. Comput. Phys.,
Volume 186, Pages 652665, 2003
Laurent F., Santoro V., Noskov M., Gomez A., Smooke
M.D. and Massot M., Accurate treatment of size distri
bution effects in polydispersed spray diffusion flames:
multifluid modeling, computations and experiments,
Combust. Theory and Modelling, Volume 8, Pages 385
412, 2004
Massot M., Laurent F., Kah D. and de Chaisemartin,
S., A robust moment method for evaluation of the
disappearance rate of evaporating sprays, SIAM Jour
nal on Applied Mathematics, 2010, available on HAL:
http://hal.archivesouvertes.fr/hal00332423/en/
Marchisio D. L. and Fox R. O., Solution of population
balance equations using the direct quadrature method of
moments, Journal of Aerosol Science, Volume 36, Pages
4373, 2005
McGraw R., Numerical advection of correlated tracers:
Preserving particle size/composition moment sequences
during transport of aerosol mixtures, Journal of Physics:
Conference Series, Volume 78, 2007
Mead L. R. and Papanicolaou N., Maximum entropy in
the problem of moments, J. Math. Phys., Pages 2404
2417, 1984
Mossa J.B., Extension polydisperse pour la description
eulereuler des 6coulements diphasiques r6actifs, Ph.D.
thesis, Institut National Polytechnique de Toulouse, 2005
Reveillon J., DNS of Spray Combustion, Dispersion
Evaporation and Combustion, In: Computational Models
for Turbulent Multiphase Reacting Flows, CISM Courses
and Lectures, Volume 442, Pages 229, SpringerWien
NewYork, Vienna, 2007, Editors Marchisio D. L. and Fox
R. O., Udine, July 2006
Williams F. A., Spray Combustion and Atomization,
Physics of Fluid, Volume 1, Pages 541545, 1958
Wright D. L., Numerical advection of moments of the
particule size distribution in Eulerian models, Journal of
Aerosol Science, Pages 352369, 2007
