D Kah

M Massot
email: marc.massot@em2c.ecp.fr

Q H Tran

S Jay

F Laurent
email: laurent@em2c.ecp.fr

A high order moment method with mesh movement for the description of a polydisperse evaporating spray

Keywords: Eulerian models, polydispersivity, high order moment method, ALE formalism

)), which is in progress.

Introduction

In order to provide an accurate numerical tool for the prediction of engines performance, the simulation of liquid 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 separate 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 models [START_REF] Greenberg | On the origin of spray sectional conservation equations[END_REF], de Chaisemartin ( 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) solution of a Williams-Fokker-Planck equation. This equation 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 collision, coalescence. However, its numerical solution in multi-dimensional systems by finite volume schemes is intractable due to the high phase space dimension. Consequently, Lagrangian approaches have been widely used and have been shown to be efficient in numerous cases, for example in [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF]) and references herein. But its main drawback is the coupling of an Eulerian description 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 formulation 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 kinetic equation using a moment method. The use of moment 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 considers 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 [START_REF] Freret | Eulerian models and three-dimensional numerical simulation of polydisperse sprays[END_REF]). 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 latter. Consequently the possibility of high order moment method considering only one size section is very attractive. Some high order moment methods have already been designed. The first one consists in solving the evolution of moments of a presumed NDF (assumed as a log-normal law) [START_REF] Mossa | Extension polydisperse pour la description euler-euler des écoulements diphasiques réactifs[END_REF]). But this approach leads to serious numerical instabilities in the treatment of evaporation. The second method uses Direct Quadrature Method of Moment (DQMOM) [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF]) where equations on nodes and abscissas of the distribution function are directly written [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF]). Contrary to the first one, this method is stable, but its accuracy can still be improved. A new high order moment method, explained in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]), provides a stable numerical scheme for evaporation, and reaches unequivalent levels of accuracy. In [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF]), 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 polydispersivity, and have proved to be a serious alternative to the multi-fluid model, without the constrain of discretization 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 Chaisemartin ( 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-sides, for engine simulations however, fuel injection in combustion chamber is done in the more general context 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 [START_REF] Donea | Arbitrary Lagrangian-Eulerian methods[END_REF]) has been adopted in the design of the numerical scheme. Indeed the ALE formalism is adapted to fluidstructure interactions as it combines advantages from the Lagrangian and the Eulerian representations. These methods are enforced in the context of engine simulation in the unstructured CFD code IFP-C3D [START_REF] Bohbot | IFP-C3D: an Unstructured Parallel Solver for Reactive Compressible Gas flow with Spray[END_REF]). In order to introduce the high order moment method described 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 stability 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 polydisperse 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 developments are validated in the context of the Riemann problem for aerosols and advection with discontinuous velocity 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 piston for aerosols, and the Riemann problem for sprays.

Moment equation and property of the moment space

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 = (x 1 , . . . , x d ) t , where d is the dimension of the physical space, with size S, and velocity 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 d 2 law, but we remark that the treatment of more complex laws is done in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]). 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:

∂ t f +∇ x •(vf )-∂ S (Rf ) = β∇ v •[((v-u g )f + Q β ∇ v •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, β -1 is the relaxation time of the particles to the gas velocity u g , Q is the matrix of the temporal correlation of forces due to brownian motion. If

F = (F 1 , F 2 , F 3 ) t is the force due to Brownian motion in three dimensions, then Q ij (t) = ∞ 0 F i (t)F j (t+τ ) dτ .
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 ′ = t τ g , x ′ = x L 0 , S ′ = S S 0 , v ′ = v U 0 , u ′ g = u g U 0 , f ′ = f U 3 0 L 3 0 S 0 n 0 , τ g = L 0 U 0 , St = U 0 β -1 L 0 , K = RL 0 S 0 U 0 , β -1 σ = Q U 2 0
where L 0 is a reference length, U 0 is a reference velocity, S 0 is the maximum droplet size, n 0 is a reference number density. The quantity τ g is the characteristic time of the particles, St is the Stokes number, K is the dimensionless evaporation rate, σ is the velocity dispersion of the particles due to Brownian motion.

Written in terms of dimensionless quantities and omitting the prime sign for the sake of legibility, the Williams-Fokker-Planck equation reads:

∂ t f + ∇ x •(vf ) -∂ S (Kf ) = 1 St ∇ v •[((v -u g )f + σ∇ v •f )]
(2) The low Stokes number of particles belonging to the aerosol class makes the rhs appear as a singular perturbation for Eq. ( 1). The NDF is then solved through a Chapman-Enskog development [START_REF] Cowling | The Mathematical Theory of Non-Uniform Gases[END_REF], [START_REF] Kaper | Mathematical theory of transport processes in gases[END_REF]). After some algebra, we find that the number density n(t, x, S) = R 3 f dv verifies the Smoluchowski equation [START_REF] Chandrasekhar | [END_REF]):

∂ t n + ∇ x •(u g n) -∂ S (Kn) = St[∇ x •(σ∇ x n) -nD t u g ], (3)
where

D t . = ∂ t . + u g ∇ x •..
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/β∇ v • f is negligible in comparison with the term (vu g )f . Thus f verifies the following kinetic equation, called Williams-Boltzmann equation [START_REF] Williams | Spray Combustion and Atomization[END_REF]):

∂ t f + ∇ x •(vf ) -∂ S (Kf ) + ∇ v •(β(u g -v)f ) = 0 (4)
Contrary to the previous case a hypothesis on the velocity 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 multifluid 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)δ(v -u(t, x, S))
, where u is the mean velocity of the particles. The semi-kinetic equation is obtained by taking the velocity moments of order 0 and 1 of Eq. (4):

∂ t n+ ∇ x •(nu)-∂ S (Kn) = 0 ∂ t nu+∇ x •(nu ⊗ u)-∂ S (Knu)- 1 St n(u g -u) = 0 (5)
In both cases, we write the resulting system of equations on the variables of interest, which are the size moments 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,

m k = 1 0 S k n(t, x, S) dS (6)
Taking the N + 1 first size moments of system (5) leads to the independant equations:

∂ t (m 0 )+ ∇ x •(m 0 u) = -Kn(t, x, 0) ∂ t (m 1 )+ ∇ x •(m 1 u) = -Km 0 ∂ t (m 2 )+ ∇ x •(m 2 u) = -2Km 1 . . . ∂ t (m N )+ ∇ x •(m N u) = -N Km N -1 ∂ t (m 1 u)+∇ x •(m 1 u ⊗ u) = m 0 m 1 (u g -u) θ (7)
We consider m 1 u as the momentum and not m 0 u, because, as St depends on S, the quantity (ug-u) St would have not been integrable. In this context, θ is such as m 1 m 0 (t = 0) θ = St. The subsystem corresponding to the 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 = 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 [START_REF]th International Conference on Multiphase Flow[END_REF], the gas phase is solution of the Euler system for compressible gas:

∂ t (ρ g )+ ∇ x •(ρ g u g ) = 0 ∂ t (ρ g u g )+∇ x •(ρ g u g ⊗ u g ) = -∇ x P g ∂ t (ρ g E g )+ ∇ x •(ρ g E g u g ) =-∇ x •(P g u g ) (8)
where ρ g , P g , E g 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 [START_REF]th International Conference on Multiphase Flow[END_REF] 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 [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]) 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.

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 M 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 m k , 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 [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]). Thus, checking if a N + 1 component vector belongs to the moment space becomes quite tedious when N ≥ 2 because of the complex geometry of the moment space. Indeed, the Hankel determinant for m k depends on all the moments of lower order. Geometrically, the interval of admissibility for m k 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 of the canonical moments, called the canonical moment space, has the very convenient property of fully lying in the set [0, 1] k , 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 c j , P (c k-1 ), where c k-1 = (c 0 , . . . , c k-1 ) t , where c k = m k m 0 are the normalized moments. Dette and Studden have given the definition of the canonical moments of order k(k ≥ 1):

p k = c k -c - k (c k-1 ) c + k (c k-1 ) -c - k (c k-1 ) (9) c + k (c k-1 ) = max µ∈P (c k-1 ) c k (µ) and c - k (c k-1 ) = min µ∈P (c k-1 )
c k (µ) are respectively the upper and lower boundary of the admissible interval for the vector c k = (1, c 1 , . . . , c k ) t to be in the moment space. The canonical moments have two major intrinsic properties 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, p k = 0 or 1, then ∀j > 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 ∀k ≥ 1, p k (f ) = p k (f SminSmax ), where f SminSmax denotes the distribution induced by the linear transformation S = S min + (S max -S min )x of [0, 1] onto [S min , S max ] (refer to [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]) 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 [START_REF] Mead | Maximum entropy in the problem of moments[END_REF]). Consequently, the value of the NDF appearing in the first equation of system [START_REF]th International Conference on Multiphase Flow[END_REF] is not exactly n, but a reconstructed distribution, we will note ñME . This closure is proven to preserve the moment space during the whole dynamical process. Moreover [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]) 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, [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]) 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 property given from the convective system extracted from system ( 7) will be very useful when designing the numerical scheme. The canonical moments are transported quantities, which means that they verify the transport equation:

∂ t p k + u∂ x p k = 0.
The demonstration can be found in [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF]). These properties of the canonical mo-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 issue discusted in [START_REF] Mcgraw | Numerical advection of correlated tracers: Preserving particle size/composition moment sequences during transport of aerosol mixtures[END_REF], [START_REF] Wright | Numerical advection of moments of the particule size distribution in Eulerian models[END_REF]).

Numerical validation of the Eulerian high order moment method

The aim of this section is to show that numerical tools have been succesfully 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 [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF]) 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 ∆t/2 the transport in phase space, then for ∆t the transport in physical space, and then for ∆t/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 [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF]) 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 Chaisemartin ( 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 injected in the jet core with either a lognormal size NDF, corresponding to the beginning of a typical experimental distribution [START_REF] Laurent | Accurate treatment of size distribution effects in polydispersed spray diffusion flames: multi-fluid modeling, computations and experiments[END_REF]). The simulations are conducted with an academic solver, coupling the AS-PHODELE solver, developed at CORIA by Julien Reveillon and collaborators [START_REF] Reveillon | DNS of Spray Combustion, Dispersion Evaporation and Combustion[END_REF]), with the multifluid solver MUSES3D developed by Stéphane de Chaisemartin and Lucie Fréret at EM2C Laboratory (de Chaisemartin ( 2009)). The ASPHODELE solver couples a Eulerian description of the gas phase with a Lagrangian description of the spray. One of the key feature of this tool is to allow, in the framework of one-way coupling, the simultaneous 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 computed on a 400 × 200 uniformly spaced grid. To destabilize the jet, we inject turbulence using the Klein method with 10% fluctuations [START_REF] Klein | A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations[END_REF]). The Reynolds number based on U 0 , ν 0 and L 0 is 1, 000. The gas vortic- Figure (2) shows the results for the 0 th order moment m 0 and the 3 rd order moment m 3 of the injected evaporating spray. The corresponding non-dimensional evaporation coefficient is K = 0.14. The mean Stokes number of the droplets is St = 1.88. These results have been validated by comparison with resuls given by the multi-fluid model [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF]).

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.

Derivation of the system in ALE formalism

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 [START_REF] Donea | Arbitrary Lagrangian-Eulerian methods[END_REF]). In the ALE description, the nodes of the computational mesh may be moved with the continuum in normal Lagrangian 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 follows 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 R X . The motion of the material points relates the material coordinates X, to the spatial ones x. It is defined by an application φ such that:

φ : R X × [t 0 , t f inal [→ R x × [t 0 , t f inal [ (X, T ) → φ(X, T ) = (x, t) (10)
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 ) = ∂ t x| X (11)
which explicitly states the particular nature of φ. 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 R X nor the spatial one R x is taken as a reference. Thus, a third domain is needed : the referential configuration R χ , where the reference coordinates, χ are introduced to identify the grid points. The referential domain R χ is mapped into the spatial domain by Φ. The mapping Φ from the referential domain to the spatial domain, which can be understood as the motion of the grid points in the spatial domain, is represented by:

Φ : R χ × [t 0 , t f inal [→ R x × [t 0 , t f inal [ (χ, t) → Φ(χ, t) = (x, t) (12)
and its gradient is defined by

∂ (χ,t) Φ = ∇ χ x u χ 0 T 1 , u χ (χ, t) = ∂ t x| χ (13)
where u χ is the mesh velocity. We define also the quantity J(χ, t) = det(∇ χ x), which is the dilatation rate of a volume d 3 χ with time. In order to close the final system, we need the time evolution of J, given by the transport theorem:

∂ t J = J∇ x •u χ (14)

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 (Després (2005)), [START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF]).

Let us consider the term α, where α can refer either to the gas density, ρ, the gas momentum, ρu g , the gas total energy, ρE g , m k,k=1..N , and m 0 u p . The equations are written over the control volume, written V (Φ(χ, t), t), whose envelope velocity is u χ . Applying the dynamic theorem leads to :

d dt V (Φ(χ,t),t) α d 3 x + ∂V (Φ(χ,t),t) α (v -u χ ) d 2 x = ∂V (Φ(χ,t),t) σ.n d 2 x + V (Φ(χ,t),t) S d 3 x ( 15
)
where σ is the generalized stress tensor and S the volumic production rate, for α.

After doing the change of variable (x, t) = Φ(χ, t) leading to the change of elementary volume d 3 x = Jd 3 χ, we obtain:

V (χ,0) ( ∂ t Jα| χ +J∇ x •α(v-u χ )-J∇ x •σ-JS) d 3 χ = 0 (16)
where V (χ, 0) is the controle volume at time t = 0.

The final model writes:

∂ t (J)| χ - J ∇ x •(u χ ) =0, ∂ t (Jρ)| χ + J ∇ x (dρw g ) =0, ∂ t (Jρu g )| χ + J ∇ x •(ρu g w g )+ J ∇ x (P )=0, ∂ t (JρE g )| χ + J ∇ x •(ρE g w g )+ J ∇ x (P u g )=0, ∂ t (Jm 0 )| χ + J ∇ x •(m 0 w p )- J K ñME =0, ∂ t (Jm k )| χ + J ∇ x •(m k w p )+J K(k -1)m k-1 =0, ∂ t (Jm 1 u p )| χ +J ∇ x •(m 1 u p w p )-J m 0 m 1 (u g -u) θ =0, (17)
where the coefficient k ranges in [1, N ], and where w g = u pu χ and w p = u gu χ .

3 Numerical scheme

Dimensional and operator splitting

In this section, we restrict ourselves to cartersian meshes to show some results for the design of a second 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:

∂ t (J)| χ - ∂ χ (u χ ) =0, ∂ t (Jρ)| χ + ∂ χ (ρw g ) =0, ∂ t (Jρu g )| χ + ∂ χ (ρu g w g )+ ∂ χ (P )=0, ∂ t (JρE g )| χ + ∂ χ (ρE g w g )+ ∂ χ (P u g )=0, ∂ t (Jm 0 )| χ + ∂ χ (m 0 w p )- J K ñME =0, ∂ t (Jm k )| χ + ∂ χ (m k w p )+J (k -1)m k-1 =0, ∂ t (Jm 1 u p )| χ +∂ χ (m 1 u p w p )-J m 0 m 1 (u g -u) θ =0, (18)
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 ∆t. This approach is fully accurate when the kinematic waves speed is negligible compared with the accoustic waves speed. This hypothesis 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 decomposed 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:

∂ t (J g )| Xg -∂ χ (u g )=0, ∂ t (J g ρ)| Xg =0, ∂ t (J g ρu g )| Xg + ∂ χ (P )=0, ∂ t (J g ρE g )| Xg +∂ χ (P u g )=0, ∂ t (J p )| Xp -∂ χ (u p )=0, ∂ t (J p m k )| Xp =0, k = 0..N ∂ t (J p m 1 u p )| Xp =0, (19)
J g denotes the dilatation due to the gaseous velocity field and J p the dilatation due to the liquid velocity field. This 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 u χ . The quantities are expressed in this final mesh, which is the reasion why this step is also called the projection step.

We divide the domain, [0, Z] into N cells [x i-1/2 , x i+1/2 ], each cell having its own size ∆x i . The innner cells are numbered from 1 to N . We also define two ghost cells labeled 0 and N + 1. The thermodynamical quantities are defined on the cells, so for example ρ i is defines as:

ρ i = 1 ∆x i x i+1/2
x i-1/2 ρ(x) dx. The velocities are defined on the nodes, we write them: u g,i-1/2 , u g,i+1/2 , u p,i-1/2 , u p,i+1/2 . Moreover, with U g = (ρ, ρu g , ρE g ), U p = (m k,,k=1..N , m 0 u p ), and the subscript m either refering to the gas (g) or the particles (p), we denote by

(J n m , U n m ), (J A m , U A m ), (J B m , U B m ), (J C m , U C m ) = (J n+1 m , U n+1 m
) the values of (J m , U m ) respectively at time t = n∆t, 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)∆t.

Resolution of the Lagrange step

Because of the different nature of the waves introduced in the equation systems for the gas and the particles in system (17), the Lagrange step is solved with dedicated solvers for each phase.

Resolution for the gas with an implicit scheme

For the gas, the Lagrangian step corresponds to the resolution 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 essential requirement.

From the mass conservation equation, we first use that the mass in a material volume is conserved: (J g ρ) B = (J g ρ) A . Extracting this relation in the three other equations leads to:

∂ t τ - 1 ρ ∂ χ u g =0, where τ = 1/ρ (20) ∂ t u g + 1 ρ ∂ χ P =0 (21)
∂ t E g + 1 ρ ∂ χ P u=0 (22) (23)
Before discretizing these equations, let us remark that the total energy E g is not well defined on the grid, as the internal energy is defined on a cell, and the kinetic energy on a node. We will then work with the gas internal energy e g , the equation of which is deduced from a combinaison of the equations of system (20). We assume that all the quantities are C 1 , and so that no shock-wave occurs. This assumption is relevant as we consider the spray downstream the injector, with a substantialy lower velocity 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 standar acoustic Godunov scheme [START_REF] Godunov | A difference scheme for numerical computation of discontinuous solution of hydrodynamic equations[END_REF]), which can be interpreted as a relaxation scheme [START_REF] Coquel | Entropy-satisfying relaxation method with large timesteps for Euler IBVPS[END_REF]). The discretized equations read:

τ B i -τ A i - ∆t(u * g,i+1/2 -u * g,i-1/2 ) ρ A i ∆x A i =0 (24) u B g,i+1/2 -u A g,i-1/2 + ∆t(P * i+1 -P * i ) ρ A i ∆x A i =0 (25) e B g,i -e A g,i +P B i (τ B i -τ A i ) =0 (26)
where

u * g,i+1/2 =u B g,i+1/2 - P B i+1 -P B i 2a (27)
P * i = P B i -a u B g,i+1/2 -u B g,i-1/2 2 (28) (29)
The parameter a is homogeneous to ρc, the Lagrangian speed of sound in the gas. However, it can shown that a must be subjected to the subcharacteristic condition

> max i (ρc) i .
The resulting non-linear matrix equation is inverted by the Newton method.

Explicit resolution for the liquid

As for particles, there is no pressure term in the momentum equation. This implies that the dilatation coefficient J p is modified only by the divergence of the particle 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 τ k,i = 1/m k,i , k = 1...N , one gets:

τ B k,i = τ A k,i (1 - ∆t ∆x i (u p,i+1/2 -u p,i-1/2 )) (30)
For τ k,i to remain positive, we impose ∆t ≤ u p,i+1/2 -u p,i-1/2 ∆x A i .

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

∂ t J m | χ +∂ χ w m =0, ∂ t J m U m | χ +∂ χ U m w m =0, (31)
Combining the equations of system (31), we obtain the componentwise advection equation, one gets:

∂ t U m + w m J m ∂ χ U m = 0, w m = u m -u χ (32)
At the discrete level, the transport velocity for the gas should be taken equal to u * g -u χ , where u * g is the output of phase B. This ensures that after phase A and B, the mesh has moved exactly at velocity u χ The projection step is merely a remap of the variables contained in U m , at the velocity relative to the referential domain.

Let ψ refer to a cell-defined quantity. As ψ is a solution of a conservative equation, a finite volume discretization is appropriate for its discretization. We denote by ∆x B m = J m ∆x n the cell volume at the end of the Lagrange step, and by ∆x n+1 the common final cell volume for both phase. Integrating the second equation of system (31) over the i th cell, we get the discretized equation:

ψ n+1 ∆x n+1 =ψ B ∆x B -∆t((ψ B w m )(x i+1/2 ) -(ψ B w m )(x i-1/2 )) (33)
The values (ψ B w m )(x i+1/2 ) and (ψ B w m )(x i-1/2 ) 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 ρu g and m 0 u p . This is not a priori obvious as ρ, m 0 , and u g , u 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 ρ and m 0 on each side of the node. For the gas, phase, this writes:

(ρug) B i+1/2 = u B g,i+1/2 ρ B i ∆x B i /2 + ρ B i+1 ∆x B i+1 /2 ∆x B i /2 + ∆x B i+1 /2 (34)
In the same time, we define the dual interfacial velocities as the average of the velocities of each regular node:

w m,i = w m,i-1/2 + w m,i+1/2 2 .
The discretized equation write, for the gas:

(ρu g ) n+1 i+1/2 (∆x n+1 i /2 + ∆x n+1 i+1 /2) = (ρu g ) B i+1/2 (∆x B i /2 + ∆x B i+1 /2)- ∆t((ψ B u g w m )(x i+1 ) -(ψ B u g w m )(x i )) (35)
As before we use a first order upwind scheme to compute the fluxes.

We may also want to have a scheme with a higher order 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 quantites. 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 quantites. 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 enforced 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.

Second order in space and time and preservation of the moment space

In [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF]), 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 [START_REF] Wright | Numerical advection of moments of the particule size distribution in Eulerian models[END_REF]) and [START_REF] Mcgraw | Numerical advection of correlated tracers: Preserving particle size/composition moment sequences during transport of aerosol mixtures[END_REF]) 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

m B k,i = 1 ∆x B x B i+1/2 x B i-1/2 m B k (x) dx, k = 0, ..N (36)
The discretization of the convection part in the equations on moments in system (17) leads to:

M n+1 i ∆x n+1 i = M B i ∆x B i -∆t(F i+1/2 -F i-1/2 ) (37)
where F i+1/2 , the moment flux, writes:

F i+1/2 = ∆t 0 R v      m 0 m 1 . . . m N      (x i+1/2 ) dvdt (38)
To design a second order scheme in space, one chooses to reconstruct the moments linearly in the i th 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 moment space is not a vectorial space.

Both difficulties are solved when considering the kinetic 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:

∂ t (f ) + v∂ χ (f ) = 0, in R χ × R v × R S ×]t n , t n+1 [ (39) so that f (t, x, v, S) = f (0, x -vt, v, S).
Thus, expressing the moments according to f leads to the following expression for F i+1/2

F i+1/2 = 1 ∆t ∆t Z 0 Z R 1 Z 0 0 B B B @ 1 S . . . S N 1 C C C A vf (t, x i+1/2 , S, 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 i+1/2 = F + i+1/2 + F - i+1/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 + i+1/2 can be expressed as a space integral of the moments over the i th cell. Provided that the reconstruction of f ensures that the scheme is conservative and that the realisability condition is satisfied ∀x ∈ [x i-1/2 , x i+1/2 ], the updated vector M n+1 i belongs to the moment space. In fact, in that case, when considering positive velocities with no loss of generality, the quantity M B i -

∆t ∆x B i F + i+1/2 is a moment vector.
To achieve a reconstruction meeting the two above requirements, we reconstruct the canonical moments, in ]x j-1/2 , x j+1/2 [ taking benefit from the fact that they are transported and they lie in [0, 1]:

         m 0 (x) = m 0,j + D m0j (x -x j ) p 1 (x) = p 1j + D p1j (x -x j ) . . . p N (x) = p N j + D p N j (x -x j )
Where D p1,i and D p N,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 conservative. Indeed, given T a transported quantity and C the corresponding conserved quantity, the scheme is conservative if:

C(x) = m 0 (x)T (x) C i = T i m 0,i = 1 ∆x x B i+1/2 x B i-1/2 m 0 (x)T (x) dx
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 p 2 , and 6 for p 3 ). The quantity p ki writes p ki = a k,i + b k,i D p k,i . The calculation of the coefficients a 2j , a 3j and b 2j , b 3j 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:

r ki ≤ p k (x) ≤ R ki , x ∈ (x i-1/2 , x i+1/2 ), where r ij = min(p k,i-1 , p k,i , p k,i+1 ) and R ij = max(p k,i-1 , p k,i , p k,i+1 )
We must have:

         r ki ≤ a k,i + b k,i D p k,i + ∆x B i 12 D p k,i ≤ R ki r ki ≤ a k,i + b k,i D p k,i - ∆x B i 12 D p k,i ≤ R ki
The slopes must then verify:

           D p k,i ≤ min R ki -a k,i b k,i + ∆x B i /2 , a k,i -m ki ∆x B i /2 -b k,i D p k,i ≥ min r ki -a k,i b k,i + ∆x B i /2 , a k,i -M ki ∆x B i /2 -b k,i
In practice, we use the following slope limiter to satisfy all the conditions:

D p k,i = 1 2 (sgn(p k,i+1 -p k,i ) + sgn(p k,i -p k,i-1 ))× min |p k,i+1 -a k,i | 2( ∆x B i 2 + b k,i ) , |a k,i -p k,i+1 | 2( ∆x B i 2 + b k,i )
The same limiter is used for the slope on the droplet number m 0 . In order to ensure the non-negativity of m 0 , one adds the following condition:D m0,i ∆x B i /2 ≤ m 0,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:

F + i+1/2 = 1 ∆t x i+1/2 x L i+1/2      m n 0 (x) m n 0 p1(x) m n 0 (p1((1 -p1)p2 + p1))(x) m n 0 (p1((1 -p1)p2(1 -p2)p3 +(p1 + p2(1 -p1)) 2 ))(x)      dx where x L i+1/2 = x i+1/2 -u m,i+1/2 ∆t
, where u m,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.

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 possibiltity to transport and evaporate and aerosol. The interest for a second order in space scheme is then 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 dimensionless. One writes P = P ′ P 0 , T = T ′ T 0 , ρ = ρ ′ ρ 0 , where P 0 , T 0 , ρ 0 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 ′ l = 3, ρ ′ l = 1, and for the right state P ′ r = 1, ρ ′ r = 0.125. The NDF for the polydisperse spray is set as n(x, S) = ρ ′ , S ∈ [0, 1] so that m 0 = ρ ′ . Physically, m 0 and ρ have no reason to be equal. m 0 is the droplet number density and ρ is the gas density. m 0 and ρ are initialized with the same value only for a numerical reason. That way, the profil of m 0 can be easily compared to the profil of ρ ′ . The particles are evaporated with the coefficient K = 2. K is such as ∆S = K∆t, so that all the droplets must be evaporated at t = ∆S/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, m 0 , and the analytical solution for this case at t = 0.1. Focusing on the resul for m 0 , it is everywhere equal to 0.8ρ ′ . This is coherent with the initial NDF, and the evaporation rate K = 2, implying that 20% of the total droplet number must have been evaporated. We have only displayed the results for m 0 and m 1 for the sake of legibility, but m 2 and m 3 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 m 0 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 displacememt of the grid follows the law : x(t) = 0.2 cos( 2π 0.1 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 profil of m 0 is much more diffused with the first order scheme than with the second order scheme.

The particles are now considered to have their own dynamics. The NDF is more complex than before, in order to show that the evaporation solver does not presume any The particles velocity is initiated by:

u(x) = 0.5 if x ≤ 0.25 2 if x > 0.25
We set periodic boundary conditions on the space interval [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 analytical 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 moment 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 0 th order moment, and the lowest curve stands for the 3 rd 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 m 0 and m 1 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 δ-shock is therefore engendered when two droplet clouds with different velocity meet [START_REF] Bouchut | On zero pressure gas dynamics, Advances in kinetic theory and computing[END_REF]). 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 Chaisemartin ( 2009)) and references therein for details on this matter. 4 Polydispersivity in a CFD code

Introduction of polydispersivity in IFP-C3D

An Unstructured Parallel Solver for Reactive Compressible 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

Test-cases

The implementation of the moment method is validated on two academic test cases. The first test case considers an aerosol in a homogeneous high pressure cell with a moving piston. The second test case considers the Riemann 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 bondary 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 stationnary. 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. [START_REF]th International Conference on Multiphase Flow[END_REF].

The size distribution is taken constant like in the one-dimensional tests with aerosol. the results are displayed for the number density m 0 and the surface density m 1 with a 1200(30 × 40) cell grid. 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 hypothesis made in the ALE formalism assuming that the acoustic time scale is negligible in comparison to the convective time scale. This hypothesis 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. 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.

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 [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF]) 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 [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]), 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 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.

Figure 1 :

 1 Figure 1: Gas vorticity field of the free-jet configuration at final time (t = 20), on a 400 × 200 grid

Figure 3 :

 3 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 represented by the solid black line, the dashed maroon line represents the numerical solution for the gas, the blue curve with stars markers and the red curve with up triangle markers respectively stands for m 0 and m 1 . The computation is carried out on a 200 cell grid

Figure 4 :

 4 Figure 4: Focus on the interfacial area in Fig.(3) andcomparison of the first and second order solution for m 0 , with and without mesh movement. The solid blue line corresponds to the second order solution without mesh movement, the dashed red line to the second order solution with mesh movement, the dotted violet line to the first order solution without mesh movement, and the green line with square markers to the first order solution with mesh movement.The computation is carried out on a 200 cell grid

Figure 5 :

 5 Figure 5: Evolution of a spray in a discontinuous velocity field calculated with the moment method, compared 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 0 th order m 0 (blue curve with circles) to the 3 rd m 3 (violet curve with diamonds). The analytical solution is represented by marquers. The computation is carried out on a 200 grid.

Figure 6 :

 6 Figure 6: Evolution of a spray in a discontinuous velocity field calculated with the moment method, compared to the analytical solution of the problem, at time t = 0.45. For legibility, only m 0 (blue curve) and m 1 (red curve) are represented. 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

Figure 7 :

 7 Figure 7: Considered piston movement during the computation. The computation starts at cad = -180 and ends at cad = 180.

  Figures (8)-(top left) displays the initial conditions for the moments, and the other charts of 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    (m 0 , m 1 , m 2 , m 3 ) = (1, 0.5, 0.333, 0.25), ∀x ∈ [0, 1] v = 1000, x ∈ [0.4, 0.6] v = 500 otherwise (41)

Figure 8 :

 8 Figure 8: Solution for aerosol particles in the case of piston movement, where m 0 and m 1 are represented. (top left) Initial condition for the homogeneous 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 × 40 grid

Figure 9 :

 9 Figure 9: Solution for spray particles in the Riemann problem. (Top) Solution for m 0 and m 1 . (Bottom) Solution for the velocity field. (left) solution at time t = 10 -4 . (Right) Solution at time t = 10 -3 . The computation is carried out on a 100 grid

th International Conference on Multiphase Flow, ICMF 2010, Tampa, FL, May 30 -June 4, 2010

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.