Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 2.4.1 - A high order moment method with mesh movement for the description of a polydisperse evaporating spray
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00053
 Material Information
Title: 2.4.1 - A high order moment method with mesh movement for the description of a polydisperse evaporating spray Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Kah, D.
Massot, M.
Tran, Q.H.
Jay, S.
Laurent, F.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: Eulerian models
polydispersivity
high order moment method
ALE formalism
 Notes
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 multi-fluid 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 multi-fluid model. The objective is to write a polydisperse spray model in a formulation close to the usual two-fluid 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 multi-fluid approach with a single section would be too coarse, standard two-fluid 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, IFP-C3D (Bohbot et al. (2009)), which is in progress.
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: VID00053
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 241-Kah-ICMF2010.pdf

Full Text
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 Bois-Preau, 92852 Rueil-Malmaison, FRANCE
Laboratoire EM2C, CNRS, Grande Voie des Vignes, 92295 Chatenay-Malabry, FRANCE
damien.kah@ifp.fr, marc.massot@em2c.ecp.fr, q-huy.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 multi-fluid 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 multi-fluid model. The objective is to write a polydisperse spray model in a formulation
close to the usual two-fluid 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 multi-fluid
approach with a single section would be too coarse, standard two-fluid 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, IFP-C3D (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 multi-phase 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 two-fluid approaches
whereas, even if the two-fluid approach can also describe
the disperse phase, the particulate aspect of the disperse
phase makes kinetic-based 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 multi-fluid 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 Williams-Fokker-Planck 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 particle-particle interaction terms such as col-
lision, coalescence. However, its numerical solution in
multi-dimensional 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 multi-fluid 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 log-normal
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 multi-fluid 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 Lagrangian-Eulerian (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 IFP-C3D (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-
Fokker-Planck based equation:

a9f+Vd(vf)-Gs(Rf) = V,[((v-u,)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 one-way 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 Williams-Boltzmann 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 semi-kinetic 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 semi-kinetic 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-
Fokker-Planck equation reads:
1
9t f + V(vf)-) s(Kf) V= V [((v g) f + aV-f)]
(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 Chapman-Enskog 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


-NKmNi-I


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 one-to-one 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 k-1 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) (C-1)

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(ck-1)
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 pressure-less gas dynamics
equation. The properties of the canonical moments
provides a second order formulation of the fluxes which
preserves the moment space. For a two-dimensional
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 Runge-Kutta 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 2-D 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 one-way 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 2-D Cartesian low Mach number
Navier-Stokes 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 free-jet 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 free-jet
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 non-dimensional 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 multi-fluid
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
Lagrangian-Eulerian (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(v-ux)-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
one-dimensional 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 one-dimensional 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 p-systems, 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 time-step 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 600ms-1. 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 time-step 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 p-system, 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 shock-wave 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 P-Bt 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,i-1/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 non-linear 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 cell-defined 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,i-1/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 off-centered 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
cell-defined 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 cell-defined 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 cell-defined 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/2-Fi-1/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 ['i-1/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
]Xj-1/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 non-negativity 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 6-shock 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
IFP-C3D

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 IFP-C3D
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
non-conventional Diesel engines using for instance HCCI
combustion mode, to direct injection hydrogen internal
combustion engines. Large super-scalar computers up to
1 000 processors are being widely used and IFP-C3D has
been optimized for running on these Cluster machines.
IFP-C3D 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 Test-cases

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
one-dimensional 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 300ms-1. 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 6-shock 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., IFP-C3D: 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
135-139, 2003

Bouchut F., On zero pressure gas dynamics, Advances
in kinetic theory and computing, World Sci. Publishing,
Pages 171-190, 1994

Coquel F., Nguyen Q. L., Postel M. and Tran, Q. H.,
Entropy-satisfying 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


0025-5718(10)02339-2, 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 265-276, 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.archives-ouvertes.fr/hal-
00443982/en/

de Chaisemartin S., Laurent, F., Massot, M. and
Reveillon, J., Evaluation of Eulerian multi-fluid ver-
sus Lagrangian methods for the ejection of polydis-
perse evaporating sprays by vortices, European project re-
port, TIMECOP-AE Project (2007). Available on HAL:
http://hal.archives-ouvertes.fr/hal-00169721/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 Non-Uniform Gases, Cambridge University Press,
1970

Despr6 B. and Mazeran C., Lagrangian gas dynamics in
two-dimensions and Lagrangian systems, Arch. Ration.
Mech. Anal., Volume 178, number 3, Pages 327-372,
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 Rodriguez-Ferran
A., Arbitrary Lagrangian-Eulerian methods, Encyclope-
dia of Computational Mechanics, Pages 413-437, 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 90-96, 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 multi-fluid method,
Journal of Computational Physics, Volume 227, Pages
3058-3088, 2008

Freret L., de Chaisemartin S., Reveillon J., Laurent F.
and Massot M., Eulerian models and three-dimensional
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 271-306, 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 Lagrangian-Eulerian computing method for all flow
speeds, J. Comput. Phys., Volume 14, Pages 227-253,
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, North-Holland, 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 652-665, 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:
multi-fluid 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.archives-ouvertes.fr/hal-00332423/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
43-73, 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
euler-euler 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 541-545, 1958

Wright D. L., Numerical advection of moments of the
particule size distribution in Eulerian models, Journal of
Aerosol Science, Pages 352-369, 2007




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 - - mvs