Mohamed Essadki
email: mohamed.essadki@centralesupelec.fr

Stephane De Chaisemartin
email: stephane.de-chaisemartin@ifpen.fr

F R Éd

Érique Laurent
email: frederique.laurent@centralesupelec.fr-marc.massot@centralesupelec.fr

Marc Massot

HIGH ORDER MOMENT MODEL FOR POLYDISPERSE EVAPORATING SPRAYS TOWARDS INTERFACIAL GEOMETRY DESCRIPTION

Keywords: High order moment method, moment space, realizable high order numerical scheme, polydisperse spray, evaporation, entropy maximization, interface geometry AMS subject classifications. 76T10, 35Q35, 65M08, 65M12, 65M99, 65D99

In this paper we propose a new Eulerian model and related accurate and robust numerical methods, describing polydisperse evaporating sprays, based on high order moment methods in size. The main novelty of this model relies on the use of fractional droplet surface moments and their ability to predict some geometrical variables of the droplet-gas interface, by analogy with the liquid-gas interface in interfacial flows. Evaporation is evaluated by using a Maximum Entropy (ME) reconstruction. The use of fractional moments introduces some theoretical and numerical difficulties. First, relying on a study of the moment space, we extend the ME reconstruction to the case of fractional moments. Then, we propose a new high order and robust algorithm to solve the moment evolution due to evaporation, which preserves the structure of the moment space. It involves some negative order fractional moments for which a novel treatment is introduced. The present model and numerical schemes yield an accurate and stable evaluation of the moment dynamics with minimal number of variables, as well as computational cost, but also provides an additional capacity of coupling with diffuse interface model and transport equation of averaged geometrical interface variables, which are essential in order to describe atomization.

1. Introduction. Modeling ans simulating the fuel injection in automotive engines and particularly Diesel engine, where the fuel is stored as a liquid phase and injected at high pressure in the combustion chamber, faces major challenges, because of the strongly multi-scale character of this two-phase flow problem. It is however essential, since it has a direct impact on the combustion regime and pollutant emissions and because of the difficulties of getting accurate measurements inside an engine and the high cost of experiments [START_REF] Pickett | Engine combustion network[END_REF].

The "Direct Numerical Simulations" approach aims at resolving the monophasic Navier-Stokes equation as well as the interface dynamics and the whole range of scales and relies on several classical methods ranging from Volume of Fluid to level set methods and combinations of both, thus constituting a whole branch of the two-phase literature (see [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. application to atomization process[END_REF][START_REF] Arienti | An embedded level set method for sharp-interface multiphase simulations of Diesel injectors[END_REF] and references therein). However, it has a hard time in realistic configurations of high Reynolds and Weber numbers, where it cannot resolve the scale spectrum of both the separated-phases two-phase flow close to the injector and the disperse phase, after atomization leading to a polydisperse evaporating spray in the downstream region, where evaporation and combustion are taking place. Therefore, reduced-order models have to be considered [START_REF] Andreini | Development of a turbulent liquid flux model for eulerian-eulerian multiphase flow simulations[END_REF], since a full resolution is out of reach for both academic and industrial applications.

Two reduced-order model classes have been designed depending on the flow re-gion. In the dense core region, diffuse interface models (see [START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF][START_REF] Touze | Coupling between separated and dispersed two-phase flow models for the simulation of primary atomization in cryogenic combustion[END_REF][START_REF] Saurel | Diffuse-interface capturing methods for compressible two-phase flows[END_REF] for a literature review) can be used to simulate interfacial flows with lower computational cost, even if the artificial fluid mixing and mesh size lead to some level of diffusion of the interface, and we lose important information about the interface geometry. To gain more precision in capturing atomization, some models add a transport equation of the interface area density [START_REF] Vallet | Development of a eulerian model for the "atomization" of a liquid jet[END_REF][START_REF] Jay | Combined surface density concepts for dense spray combustion[END_REF][START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF][START_REF] Devassy | Atomization modelling of liquid jets using a two-surface-density approach[END_REF], which is not not sufficient to describe a key property for applications: the polydispersion in the disperse phase region.

The second reduced-order model class relies on a kinetic approach, which has been widely used to describe accurately a population of particle at mesoscopic level. In this approach, a cloud of droplets is modeled by a number density function (NDF), which satisfies a Williams-Boltzmann Equation [START_REF] Williams | Spray combustion and atomization[END_REF]. The internal variables of the NDF provide a statistical description of some relevant physical properties such as the droplet size, velocity and temperature. In the following, we focus on size and velocity but temperature can easily be introduced. The numerical resolution of the Williams-Boltzmann Equation can be achieved by the Lagrangian Monte-Carlo approach. This method is considered to be the most accurate for solving this equation, but leads to a high computational cost for unsteady flows and requires complex load-balancing and coupling algorithms in parallel computations. An alternative method consists in deriving an Eulerian moment model from the Williams-Boltzmann Equation. In this approach, a differential system of a finite set of moments of the NDF is closed by expressing the velocity NDF conditioned in size as a function of the known velocity moments [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF][START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Levermore | The Gaussian moment closure for gas dynamics[END_REF][START_REF] Grad | Principles of the kinetic theory of gases[END_REF][START_REF] Sabat | Eulerian modeling and numerical methods for the description of turbulent polydisperse sprays[END_REF]. For the modeling of the size distribution, three possible approaches can be used: 1-Multi-fluid models also called sectional methods (see [START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays: modeling and numerical issues[END_REF][START_REF] Laurent | Two-size moment Eulerian multi-fluid model: a robust and high-fidelity description of polydisperse moderately dense evaporating sprays[END_REF] and references therein), where the size phase space is discretized into intervals called sections using the conservation of size moments up to two. 2-Method of moments with interpolative closure [START_REF] Frenklach | Method of moments with interpolative closure[END_REF], essentially used for soot modeling and simulation, provides a closure of negative as well as fractional moments from the integer ones through a logarithm Lagrangian interpolation. However, such an approach suffers from several issues such as handling multi-modal distributions, preserving the moment space [START_REF] Mueller | Hybrid method of moments for modeling soot formation and growth[END_REF], as well as important inaccuracy due to the way the negative order moments are approximated. 3-The third approach involves high order moments on the whole size range, which will be a compact interval in the present study, as well as a procedure in order to reconstruct the NDF. It allows to limit the computational cost and the number of systems of partial differential equations to be solved. We will not consider in the present study quadrature methods such as QMOM or EQMOM [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments for aggregation-breakage processes[END_REF][START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF][START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF][START_REF] Pollack | Bivariate extensions of the extended quadrature method of moments (eqmom) to describe coupled droplet evaporation and heat-up[END_REF] since they lead to either difficulties in evaluating the evaporating flux at zero size or can not reach the integrality of the moment space, but rather focus on Maximum Entropy (ME) reconstruction such as with Eulerian Multi Size Moment (EMSM) and Correlated Size Velocity Moment models in [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays: modeling and numerical issues[END_REF][START_REF] Gumprich | Assessment of the direct quadrature-based sectional method of moments for the simulation of evaporating polydisperse sprays[END_REF]. They cover the whole moment space and show a great capacity in modeling the polydispersion and evaporation with a minimal number of variables.

So far, the existing models do not provide a unified description of the two regions. In the present contribution, we propose a new high order moment model for the disperse phase with the capacity of describing the interface geometry statistics [START_REF] Drew | Evolution of geometric statistics[END_REF]. Our strategy is to resolve the polydispersion by using a set of variables, which can be identified as averages of the gas-liquid interface geometry. We show that some geometrical variables can be expressed as fractional surface-moments of NDF in the disperse zone. We introduce the mathematics fundamentals of the model and show that we can preserve the advantages of previous methods [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] in terms of both accuracy, realizability and robustness, and low computational cost but with a much higher potential in terms of coupling with a diffuse interface model. Evaporation flux are evaluated using ME reconstruction based on fractional moments1 . The existence and uniqueness of this convex optimization problem under constraints is given in [START_REF] Mead | Maximum entropy in the problem of moments[END_REF] in the case of integer moment. While some elements of proof are to be found in the fractional moment case in [START_REF] Kapur | Entropy optimization principles with applications[END_REF], we propose a generalization of the result in the case of a special set of fractional moments. Moreover, we generalize useful properties of the fractional moment space such as canonical moments as well as lower principal representation and quadrature [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]8,[START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. These properties are relevant ingredients to design high order realizable schemes and algorithms to solve a high order moment system. A new realizable algorithm to solve the evolution of fractional moments due to the evaporation is proposed, which involves negative order of moments and requires an original strategy compared to the integer moment problem. The accuracy and robustness of the proposed strategy is then assessed by a careful investigation of the numerical errors as well as a detailed comparison with the original approach in 0D, 1D and 2D academic configurations. A companion paper [START_REF] Essadki | High order moment methods and adaptive mesh refinement for polydisperse evaporating spray simulations[END_REF] aims at implementing the proposed model and numerical strategy into a massively parallel code with adaptive mesh refinement and showing its potential towards realistic engine simulations.

The paper is organized as follows. First, the two-phase flow modeling of polydisperse evaporating sprays in a carrier gaseous flow field as well as the original high order moment modeling are introduced; we also recall the classical averaged geometrical description of interfacial flows in order to identify the relevant geometrical averaged variables. Section 3 is devoted to the introduction of the new geometrical high order moments for a polydisperse spray as well as the resulting system of partial differential equations on moments derived from the Williams-Boltzmann Equation. The closure of the system through ME is then presented and its mathematical properties detailed. Once the proposed system is closed, Section 4 and 5 are dedicated to its numerical resolution. While Section 4 is devoted to the transport in phase space, Section 5 focuses on the transport in physical space. Section 6 is eventually concerned with the numerical verification and results in 0D, 1D and 2D, thus assessing the proposed modeling and numerical strategy, before concluding.

Two phase flows modeling.

2.1. Kinetic modeling of polydisperse spray and semi-kinetic model. The spray consists in a cloud of polydisperse droplets, which can be described statistically with the number density function (NDF) f φ (t, x, c, φ, T ). This function represents the probable number of droplets located at position x, travelling with velocity c and having temperature T and size φ. In general, the size φ of a spherical droplet can be given by its volume V , its surface S or its radius R of the droplet. By considering a spherical form, these three geometrical variables are equivalent

f V dV = f S dS = f R dR.
In the following, we use the surface S as the size variable. The NDF will be simply noted by f . In the following, dimensionless variables are considered, in such a way that S ∈ [0, 1].

For the sake of simplicity and the clarity of the presentation, the derivation of the model is done here in a simplified context: we consider a dilute spray at high Knudsen number and small and spherical droplets at low Weber number, obtained after the atomization. Under these assumptions, secondary breakup, coalescence and collision can be neglected. We also assume that thermal transfer can be neglected, in such a way that the variable T is no longer considered, that the only force acting on the droplets is the drag force, modeled by the Stokes law and we will mainly focus on one-way coupling. Then, the NDF satisfies the following simplified and dimensionless form of the Williams-Boltzmann Equation [START_REF] Williams | Spray combustion and atomization[END_REF]:

(1)

∂ t f + ∂ x • (cf ) + ∂ c • U g -c St(S) f + ∂ S (R S f ) = 0,
where St(S) = θS is the Stokes number and R S is the dimensionless evaporation rate, assumed to be constant: R S (S) = -K. We refer the readers to the following articles and references therein [START_REF] Laurent | Two-size moment Eulerian multi-fluid model: a robust and high-fidelity description of polydisperse moderately dense evaporating sprays[END_REF][START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF][START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF], showing that such a mesoscopic approach is capable of describing droplet heating, coalescence and break-up and two-way coupling.

The high dimension of the phase space makes its direct discretization too costly for complex applications. Since the resolution of the NDF is not required and only macroscopic quantities of the flows are needed for applications, an Eulerian moment method can be used to reduce the complexity of the problem. The closure of the velocity distribution requires a specific treatment, especially when we are concerned with modeling particle trajectory crossing at high Knudsen numbers and moderate to large inertia. For accurate modeling of such phenomenon, one can use high order velocity-moment based on classical equilibrium approaches in kinetic theory, and we rather rely on the Levermore hierarchy [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Levermore | The Gaussian moment closure for gas dynamics[END_REF][START_REF] Sabat | Eulerian modeling and numerical methods for the description of turbulent polydisperse sprays[END_REF]. In the present work, we are not concerned with these modeling issues, and we will consider a monokinetic velocity distribution [START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF][START_REF] Sabat | Eulerian modeling and numerical methods for the description of turbulent polydisperse sprays[END_REF]:

(2) f (t, x, c, S) = n(t, x, S)δ(c -U (t, x, S)).
This closure does not take into account particle trajectory crossing, since only one velocity is defined per position and time. This assumption is valid for low inertia droplets, when the droplet velocities are rapidly relaxed to the local gas velocity [START_REF] Sabat | Eulerian modeling and numerical methods for the description of turbulent polydisperse sprays[END_REF].

We then derive the following semi-kinetic system from equation (1) by considering moments of order 0 and 1 in velocity:

(3)

∂ t n + ∂ x • (nU ) = K ∂ S n, ∂ t nU + ∂ x • (nU ⊗ U ) = K ∂ S (nU ) + n U g -U St(S) .
2.2. High order size-moment model: closure and moment space. In the present contribution, we adopt the high order size-moment method with a continuous reconstruction of the size distribution to model a polydisperse spray. This approach, based mainly on [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF][START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF], consists in deriving a dynamical system on a finite set of size-moments of the NDF. The integer size-moments are defined as follows, with N ≥ 1 for a non dimensional size interval [0, 1], m = (m 0 , . . . , m N ) t , with m k = 1 0 S k n(t, x, S)dS. For the so-called Eulerian Multi Size Moment (EMSM) model, the droplet velocity is assumed to be independent on the size: U (t, x). Vié at al. [START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays: modeling and numerical issues[END_REF] developed a model able to deal with this dependence, but it is beyond the scope of this paper. The system of the EMSM model obtained from an integration of the semi-kinetic system (3) over S ∈ [0, 1] multiplied by S k reads:

(4)

∂ t m 0 + ∂ x .(m 0 U ) = -Kn(t, x, S = 0) , ∂ t m k + ∂ x .(m k U ) = -kK m k-1 , ∂ t (m 1 U ) + ∂ x .(m 1 U ⊗ U ) = -K m 0 U +m 0 U g -U θ .
where k = 1 . . . N . The unclosed term -Kn(t, x, S = 0) expresses the pointwise disappearance flux of droplets through evaporation. Even though, this term is only involved in the first equation, it contributes in the evolution of the other moments in the same way through the terms -kKm k-1 [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]. Before proposing a closure of this model, let us recall the definition of the moment space and some useful properties. We denote by P the set of all probability density measures of the interval [0, 1]. Then the Nth "normalized" moment space M N is defined as follows:

(5)

M N = {c N (µ), µ ∈ P } , c N = (c 0 (µ), . . . , c N (µ)), c k (µ) = 1 0 x k dµ(x).
We have c 0 = 1, since we use probability density measures. We can also normalize by m 0 the moment vector (m 0 , . . . , m N ) t and obtain (c 0 , . . . , c N ) t ∈ M N , where c k = m k /m 0 . The Nth "normalized" moment space is a convex and bounded space. Definition 1. The N th moment space is defined as the set of vectors m, for which the vector normalized by m 0 belongs to the N th normalized moment space. Considering this definition and some results from [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]: if (m 0 , . . . , m N ) t is in the interior of the moment space, there exists an infinity of size distributions, which represent this moment vector. In other words, there exists an infinity of size distributions n(S), which are the solution of the following finite Hausdorff moment problem:

(6) m k = 1 0 S k n(S)dS, k = 0 . . . N.
Massot et al [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] proposed to use a continuous reconstruction of the size distribution through the maximization of Shannon entropy:

(7) H(n) = - 1 0 n(S) ln(n(S))dS.
The existence and uniqueness of a size distribution n M E (S), which maximizes the Shannon entropy [START_REF] Tagliani | Hausdorff moment problem and maximum entropy: a unified approach[END_REF] and is the solution of the finite Hausdorff moment problem [START_REF] Descombes | Analysis of operator splitting in the non-asymptotic regime for nonlinear reaction-diffusion equations. application to the dynamics of premixed flames[END_REF] was proved in [START_REF] Mead | Maximum entropy in the problem of moments[END_REF] for the moments of integer order, when the moment vector belongs to the interior of the moment space and the solution is shown to have the form:

(8) n(S) = exp -(λ 0 + λ 1 S + . . . + λ N S N ) ,
where coefficients λ k , are determined from [START_REF] Descombes | Analysis of operator splitting in the non-asymptotic regime for nonlinear reaction-diffusion equations. application to the dynamics of premixed flames[END_REF]. The resolution of this nonlinear problem can be achieved by using Newton-Raphson method. The limitation of this algorithm when the moment vector is close to the moment space boundary, or equivalently when the ME reconstructed size-distribution degenerates to a sum of Dirac delta functions, is discussed in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]. Vié et al [START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays: modeling and numerical issues[END_REF] proposed some more advanced solutions to cope with this problem, by tabulating the coefficients depending on the moments and by using an adaptive support for the integral calculation, which enables an accurate computation of the integral moments when the NDFs are nearly singular.

2.3. Limitation of such an approach to a disperse phase model. Even though this high order moment formalism provides some key informations about polydispersion using only one size-section, it is restricted to the disperse phase zone. Coupling such an approach with a separated-phases model requires some complementary informations, which the usual approaches of diffuse interface models can not provide. Indeed, diffuse interface models [START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF][START_REF] Saurel | Diffuse-interface capturing methods for compressible two-phase flows[END_REF] consider the interface as a smooth transition layer, where we have lost important informations about the interface geometry. The first step, would be to enrich the diffuse interface models [START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF][START_REF] Essadki | Statistical modeling of the gas-liquid interface using geometrical variables: toward a unified description of the disperse and separated phase flows[END_REF] in order to transport averaged geometrical variables to gain accuracy. Nevertheless, even in the framework of such an enriched diffuse interface model, the coupling of two very different models is usually cumbersome and relies on parameters the described physics will depend on.

Consequently, we adopt an original strategy and build a high order moment model for the disperse phase, which possesses the same key properties as the EMSM model in terms of accuracy, robustness and computational cost, but involves a different set of variables that are describing the averaged interface geometry, so that we end up with a set of common variables in the two zones and can potentially help in building a single unified model able to capture the proper physics in both zones. In order to introduce the new set of variables, we first have to recall the natural geometrical variables in the separated-phases zone, before extending this description to the disperse phase.

2.4. Geometrical description of interfacial flows. In many two-phase flow applications, the exact location of each phase is difficult to determine precisely because of the different unpredictable phenomena such as turbulence, interface instabilities and other small scale phenomena. Fortunately, in industrial applications, we are more concerned with the averaged features than to the small details of the flow. Therefore, we can use Diffuse Interface Models [START_REF] Drew | Mathematical modeling of two-phase flow[END_REF][START_REF] Ishii | Thermo-fluid dynamics of two-phase flow[END_REF][START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF][START_REF] Saurel | Diffuse-interface capturing methods for compressible two-phase flows[END_REF] to describe the interface location in terms of probability and averaged quantities based on averaged operators (ensemble averaged, time averaged or volume averaged). In the following, we define some averaged geometrical variables, which can be used to model the interface in separated-phases for a complementary geometrical description. Their definitions are based here on the volume averaged operator following the derivation of Drew [START_REF] Drew | Evolution of geometric statistics[END_REF]. First, we define the phase function χ k (t, x) for a given phase k by χ k (t, x) = 1 if x ∈ k and χ k (t, x) = 0 otherwise. The volume-averaged operator is defined by: (

•)(t, x) = 1 |V | V (•)χ k (t, x )dx
, where V ∈ R 3 is a macroscopic space around the position x, and |V | is the occupied volume, with η 3 << |V | << L 3 , where η (resp. L) is the fluctuations (resp. system) length scale.

The Diffuse Interface Models can be obtained by averaging the monophasic fluid equations. The obtained equations involve the volume fraction, which allows to locate the interface up to the averaging scale and is then a first piece of information about the interface geometry:

(9) α k (t, x) = 1 |V | V χ k (t, x )dx .
The second variable treated by Drew [START_REF] Drew | Evolution of geometric statistics[END_REF] and used also in other two-phase flow models [START_REF] Vallet | Development of a eulerian model for the "atomization" of a liquid jet[END_REF][START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF] is the interfacial area density. The importance of this variable relies mainly on the modeling of exchange terms (evaporation, thermal transfer and drag force) as well as modeling the primary atomization. This variable is interpreted as the ratio of the surface area of an interface contained in a macroscopic volume and this volume.

(10) Σ(t, x) = 1 |V | V ||∇χ k (t, x )||dx .
So far, the interface modeling is still incomplete, since no information on the interface shape is being given. In fact, the small details of the interface can not be modeled accurately using only two geometrical variables. Drew proposed to introduce the mean H = 1 2 (k 1 + k 2 ) and Gauss G = k 1 k 2 curvatures of the interface in his model, where k 1 and k 2 are the two principal curvatures. These variables are defined only at the interface. Therefore, we need a specific averaging for these interfacial variables. So, we introduce the interfacial averaging operator (•), defined as follows: Σ

(•)(t, x) = 1 |V | V (•)||∇χ k (t, x )||dx .
The interfacial averaged Gauss and mean curvatures weighted by the interface are defined as follows:

(11) Σ H = 1 |V | V H||∇χ k (t, x )||dx , Σ G = 1 |V | V G||∇χ k (t, x )||dx .
Drew [START_REF] Drew | Evolution of geometric statistics[END_REF] derived conservative equations for these variables with source terms, which describe the stretch and the wrinkling of the interface. Its derivation is based on a kinematic evolution of an interface, when the interface velocity is given. In real application, the interface velocity can be determined from the diffuse interface model.

In the following, we express these geometrical variables in the disperse phase as sizemoments of the NDF, and derive a new high order moment model using such variables.

3. Geometrical high order moment model.

3.1. Interfacial geometrical variables for the disperse phase. Let us consider a population of spherical droplets represented by their size distribution n(t, x, S). Then, by analogy with the separated phases, we express the averaged geometrical variables: volume fraction, interface area density, Gauss curvature and mean curvature in the disperse phase. The definition of these geometrical variables is based on the phase function χ k . This function contains all the microscopic information about the interface. In the disperse phase, we use the statistical information about the spherical droplet distribution, which is given by the size distribution n(t, x, S). Considering this function, we define the different geometrical variables in the context of a polydisperse spray as follows:

A. The volume fraction α d is the sum of the volume of each droplet divided by the contained volume at a given position:

(12) α d = 1 0 V (S)n(t, x, S)dS, V (S) = S 3/2 6 √ π .
B. The interfacial area density Σ d is the sum of the surface of each droplet divided by the contained volume at a given position:

(13) Σ d = 1 0 Sn(t, x, S)dS.
C. The two local curvatures are equal for a spherical droplet k

1 = k 2 = 2 √ π √ S .
But since we use the mean and Gauss curvatures, we can define two different averaged quantities. Let us notice that, in the case of separated phases, the average mean and Gauss curvatures were defined as an average over a volume and weighed by the interfacial area. In the disperse phase case, this becomes:

(14) Σ d H d = 1 0 H(S)Sn(t, x, S)dS, Σ d G d = 1 0 G(S)Sn(t, x, S)dS.
These four geometrical variables are expressed as fractional moments of the size distribution m k/2 = 1 0 S k/2 n(S)dS:

(15) Σ d G d = 4πm 0/2 , Σ d H d = 2 √ πm 1/2 , Σ d = m 2/2 , α d = 1 6 √ π m 3/2 .
These moments can be expressed as integer moments by simple variable substitution x = √ S. However, we prefer to hold the droplet surface as the size variable, since we consider a d 2 evaporation law, where the evaporation rate R S is constant.

3.2. The governing moment equation. In this section, we derive from the kinetic equation ( 1) a high order fractional moment model. This model gives the evolution of the mean geometrical interfacial variables due to transport, evaporation and drag force and reads:

(16)          ∂ t m 0/2 + ∂ x • m 0/2 U = -Kn(t, x, S = 0), ∂ t m k/2 + ∂ x • m k/2 U = - kK 2 m (k-2)/2 , ∂ t m 2/2 U + ∂ x • m 2/2 U ⊗ U = -Km 0/2 U +m 0/2 U g -U θ ,
where k ∈ {1, 2, 3} and -Kn(t, x, S = 0) represents the pointwise disappearance flux, and the moment of negative order, m -1/2 = 1 0 S -1/2 n(t, S)dS, naturally appears in the system after integrating by part the evaporation term in the Williams-Boltzmann Equation. In the following, the source terms will be closed by a smooth reconstruction of the size distribution through entropy maximisation.

The use of fractional moments introduces a new mathematical framework as well as some numerical difficulties, which require a specific treatment. Some useful mathematical properties of the fractional moments space are discussed in the Appendix B. These results will be used to design realizable numerical schemes, i.e. schemes that ensure the preservation of the moment vector in the moment space.

3.3. Maximum Entropy reconstruction. NDF reconstruction through the maximum entropy (ME) provides a smooth reconstruction to close the moment system [START_REF] Frenklach | Method of moments with interpolative closure[END_REF] as it was done in the EMSM model. The ME reconstruction consists in maximizing the Shannon entropy [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF], under the condition that the first N + 1 (here N = 3) moments of the size distribution are equal to the computed moments:

(17) m k/2 = 1 0 S k/2 n(S)dS, k = 0 . . . N.
In the Appendix C.1, the existence and uniqueness of the ME distribution function is shown, as presented in the following Theorem. The proof is based on the method of Lagrange multipliers, using some ideas of Mead and Papanicolaou [START_REF] Mead | Maximum entropy in the problem of moments[END_REF], who already showed this result in the case of integer moments, but in a different and simplified way. Indeed, in our case, we only use the definition of the fractional moment space, given in the Appendix B. Theorem 2. If the vector m N = (m 0 , m 1/2 , .., m N/2 ) belongs to the interior of the Nth fractional moment space, then there exists a unique NDF n M E maximizing the Shannon entropy defined by [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] under the constraints [START_REF] Gordon | Error bounds in equilibrium statistical mechanics[END_REF]. Moreover, it has the following form:

(18) n M E (S) = exp -λ 0 - N i=1 λ i S i/2 .
Moreover, in the Appendix C.2, an algorithm is given to compute the parameters λ i from the fractional moments.

4. Evaporation source term. In this section, we focus on the numerical resolution of the evaporation source terms. The kinetic equation then reads:

(19) ∂ t n -∂ S (Kn) = 0.
The integration of this equation multiplied by the vector (1, S 1/2 , S, S 3/2 ) t yields to a system of ordinary differential equations for the moments defined on [0, 1]:

(20) d t m 0/2 = -n M E | S=0 , d t m 1/2 = -K 2 m -1/2 , d t m 2/2 = -Km 0/2 , d t m 3/2 = -3K 2 m 1/2 ,
where the

n M E | S=0 (m 0/2 , m 1/2 , m 2/2 , m 3/2
) is obtained by maximum entropy reconstruction. Solving this system using classical integrators such as Euler or Runge-Kutta methods does not ensure the preservation of the moments in the moment space [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]. This property is essential for robustness and accuracy to reconstruct a positive NDF.

4.1. Exact kinetic solution through the method of characteristics. The exact solution of the NDF evolution, in the case of d 2 evaporation law, can be easily obtained by solving analytically the kinetic equation ( 19): [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow[END_REF] n(t, S) = n(0, S + Kt).

For more general evaporation law, when the evaporation rate R S (S) is a smooth function of the size, the kinetic equation can be solved by using the method of characteristics. Indeed, by multiplying the kinetic equation by R S (S), we obtain the following equation in the function Γ(t, S) = R S (S)n(t, S):

(22) ∂ t Γ(t, S) + R S (S)∂ S (Γ(t, S)) = 0.
Then, for a given initial time t 0 and size S o , we define the one variable function g(t) = Γ(t, S(t; t 0 , S o )), such that S(t; t 0 , S o ) is the characteristic curve verifying:

(23)    d S(t; t 0 , S o ) dt = R S ( S(t; t 0 , S o )), S(t 0 ; t 0 , S o ) = S o .
It is a simple fact that the derivative of g(t) vanishes, which implies g(t) = g(t 0 ) and then Γ(t, S(t; t 0 , S o )) = Γ(t 0 , S o ). Finally, we obtain the exact solution of the size distribution as follows:

(24) n(t, S) = R S ( S(t 0 ; t, S)) R S (S) n(t 0 , S(t 0 ; t, S)).
4.2. Inefficiency of the original evaporation scheme. It is possible to use the exact solution [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow[END_REF] or ( 24) of [START_REF] Gumprich | Assessment of the direct quadrature-based sectional method of moments for the simulation of evaporating polydisperse sprays[END_REF], at the kinetic level, to compute directly the moments. In this case, a Gauss-Legendre quadrature can approximate the momentintegrals. In [START_REF] Mead | Maximum entropy in the problem of moments[END_REF], the authors showed that 24 quadrature points are needed to have an accurate estimation of such integrals. However, when the evaporation rate R S is not a constant, we will also need to solve the ODE [START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF] for each quadrature points. For this reason, we consider this direct method in computing the moments will be very costly in terms of CPU time for realistic simulations.

Massot et al [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] proposed a realizable algorithm to solve the evaporation moment system in the case of integer moments. This algorithm is based on evaluating the disappearance evaporation flux, that is the droplets that will completely evaporate during one time step, then use a simplified quadrature to evaluate the kinetic evolution of the moments. This algorithm reads in the case of fractional moments:

• From the moment vector m(t n ) = (m 0/2 , m 1/2 , m 2/2 , m 3/2 ) t , we reconstruct the NDF by using the ME algorithm. Then, we compute the disappearance fluxes of the droplets which will be totally evaporated at the next step:

(25) Φ -(t n ) = K∆t 0 n M E (t n , s)     1 s 1/2 s s 3/2     ds.
• Using the Product-Difference algorithm [START_REF] Gordon | Error bounds in equilibrium statistical mechanics[END_REF], we compute the abscissas S j ≥ K∆t and the weights w j ≥ 0 of the quadrature corresponding to of the mo-

ments 2 m [K∆t,1] (t n ) = m(t n ) -Φ -(t n ).
This corresponds to the resolution:

(26) m [K∆t,1] k/2 (t n ) = 2 j=1 w j S k/2 j , k ∈ {0, 1, 2, 3},
where m [K∆t,1] (t n ) is the moment vector on [K∆t, 1]. • Finally, we calculate the moments at the next step:

(27) m k/2 (t n + ∆t) = 2 j=1 w j (S j -K∆t) k/2 .
However, in the case of the fractional moments, the present algorithm can not predict the correct kinetic evolution of the NDF. We illustrate this problem in Figure 1, where we compare the maximum entropy reconstructed NDF of the moment computed with the above algorithm and the exact kinetic solution of the NDF. It is clear that the solution obtained by the algorithm diverges from the exact solution. In the following, we will provide the reader with an explanation of the inefficiency of this algorithm in the case of fractional moments and propose a new specific solution.

4.3. Involvement of negative moment order. In order to understand the inefficiency of the previous algorithm proposed in section 4.2, we propose to consider again a constant evaporation rate R S (S) = -K. The moment solution expressed directly from the exact kinetic solution reads:

(28) m k/2 (t+∆t)= 1 0 S k/2 n(t+∆t, S)dS = 1 0 S k/2 n(t, S +K∆t)dS = 1+K∆t K∆t (S -K∆t) k/2 n(t, S)dS = 1 K∆t (S -K∆t) k/2 n(t, S)dS.
In the previous equality, we consider that n(t, S) = 0 for S ≥ 1. In the following, we show that the updated moments can be written as a function of an infinite set of moments on the support [K∆t, 1].

Lemma 3. For all positive integer k and for all x ∈ [-1, 1] the function f : x → (1 -x) k/2 admits a power series which converges normally to the function f :

(29) (1 -x) k/2 = +∞ n=0 a k/2 n x n ,
a k/2 n = (-1) n (m + 1/2)(m -1/2) . . . (m + 1/2 -n + 1) n! , n > 0.
By using the Stirling's approximation, for n ≥ m + 3/2, we can show the following equivalence relation when n tends to infinity.

(

) |a k/2 n | = (-1) m+1 (2(n -m -1))!(2m + 1)! (n -m -1)!m!2 2n-1 n! ∼ (2m + 1)! 4 m √ πm! 1 n 3/2+m . 31

Therefore the series |a

k/2 n | is convergent for any integer k ≥ 0, thus concluding the proof.

We deduce that ∞ n=0 a k/2 n (K∆t) n s k/2-n converges normally to (s -K∆t) k/2 for s ≥ K∆t. Thus, we can invert the sum and the integral in the moment expression:

(32) m k/2 (t + ∆t) = 1 K∆t (s -K∆t) k/2 n(t, s)ds, = ∞ n=0 a k/2 n (K∆t) n 1 K∆t s k/2-n n(t, s)ds, = ∞ n=0 a k/2 n (K∆t) n m [K∆t,1] k/2-n (t).
where m [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF] shows that the fractional moments at t + ∆t depend on the four fractional moments (m where l = -1/2, -2/2, . . . -∞). In the case of the EMSM model, where only integer moments are used, the same expansion of the exact kinetic solution of the integer moments involves only the moments of the support [K∆t, 1] and of order l ∈ {0, 1, 2, 3}. For this reason, the evolution of the integer moments can be evaluated exactly by translating the abscissas [START_REF] Laurent | Two-size moment Eulerian multi-fluid model: a robust and high-fidelity description of polydisperse moderately dense evaporating sprays[END_REF] of the quadrature [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. The same strategy for fractional moments leads to the divergence of the method observed previously. Since we aim at describing the evolution of the moments through the evolution of a minimal number of quadrature nodes, we propose a new quadrature, suitable for fractional moments, such that the updated moments in [START_REF] Laurent | Two-size moment Eulerian multi-fluid model: a robust and high-fidelity description of polydisperse moderately dense evaporating sprays[END_REF] can approximate accurately the exact kinetic evolution (32).

[K∆t,1] l/2 = 1 K∆t s k/2-n n(t, s)ds is a fractional moment of support [K∆t, 1]. Equation

4.4.

A specific quadrature for negative moment order. Our objective is to find an adequate quadrature with the lowest possible quadrature number n q , such that the following approximation:

(33) m k/2 (t + ∆t) ≈ nq j=1 w j (S j -K∆t) k/2 ,
is accurate. More precisely, we would like to find a quadrature such that the difference,

for k ∈ {0, 1, 2, 3} (34)
k/2 (∆t) = m k/2 (t + ∆t) - nq j=1 w j (S j -K∆t) k/2 = ∞ n=0 a k/2 n (K∆t) n m [K∆t,1] k/2-n (t) - nq j=1 w j S k/2-n j
, is at least o(∆t) to ensure the convergence of the numerical scheme. It is a difficult task to prove such a result, but we provide here a first result is this direction. Our strategy consists in using a quadrature, which cancels a finite set of the first terms in the sum [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments for aggregation-breakage processes[END_REF], by enforcing

E p/2 = m [K∆t,1] p/2 (t)- nq j=1 w j S p/2 j
to be zero for p ∈ {-2n - q , . . . , 2, 3}.

Thus 2n q fractional moments on [K∆t, 1] are considered, with n q = n - q +2, from order -n - q to 3/2. As soon as n - q ≥ 0, this is sufficient to have 0/2 = 0 and 2/2 = 0. The following Lemma shows the existence and uniqueness of such quadrature and gives the corresponding value of the error terms k/2 (∆t).

Lemma 4. Let m [K∆t,1] p/2 p∈{-2n - q ,...,3}
be a set of fractional moments on [K∆t, 1] in the interior of the moment space and n - q a non negative integer. There exists a unique set of abscissas (S j ) j∈{1,...,nq} and weights (w j ) j∈{1,...,nq} , with n q = n - q + 2 such that [START_REF] Mead | Maximum entropy in the problem of moments[END_REF] ∀k ∈ {-2n - q , . . . , 3} m

[K∆t] k/2 = nq j=1 w j S k/2 j .
Moreover one have w j = w j r 2n - q j and S j = r 2 j , where r j and w j are the abscissas and weights of the Gauss quadrature corresponding to the moments

m k = m [K∆t,1] k/2+n - q for k ∈ {0, . . . , 2n q -1}.
Using this quadrature, the error terms defined by (34) can be written: . Proof. To be able to compute the weights w j and abscissas S j of the quadrature [START_REF] Mead | Maximum entropy in the problem of moments[END_REF], the fractional moments are transformed to integer ones, with non negative orders. Indeed, using the change of variable r = √ S, E p/2+n - q can be written, for k ∈ {0, . . . , 2n q -1}:

(36) 0/2 (∆t) = 0, 1/2 (∆t) = +∞ n=nq-1 a 1/2 n Ē(2n-1)/2 (K∆t) 1/2 , 2/2 (∆t) = 0, 3/2 (∆t) = +∞ n=nq a 3/2 n Ē(2n-3)/2 (K∆t)
(37) E p/2+n - q = 1 K∆t s k/2-n - q n(s)ds- nq j=1 w j S k/2-n - q j = 1 √ K∆t r k µ(r)dr- nq j=1 w j r k j , with µ(r) = (2rn(r 2 ))/(r 2n - q ) for r ∈ [ √ K∆t, 1], r j = S j ∈ [ √ K∆t, 1] and w j = w j r -2n - q j
≥ 0. Then w j and r j are defined as the unique weights and abscissas of the Gauss quadrature of order n q corresponding to moments [START_REF] Andreini | Development of a turbulent liquid flux model for eulerian-eulerian multiphase flow simulations[END_REF] for k ∈ {0, . . . , 2n q -1}. By using this quadrature in equation ( 34), we cancel the k + n - q first terms in the series, thus obtaining [START_REF] Mueller | Hybrid method of moments for modeling soot formation and growth[END_REF]. The quantities Ēp/2 are bounded, since S j ≥ K∆t and

m k = m [K∆t,1] k/2+n - q of µ on [ √ K∆t,
| Ēp/2 | ≤ max    (K∆t) p/2 m [K∆t,1] -p/2 , (K∆t) p/2 nq j=1 w j S -p/2 j    ≤ m [K∆t,1] 0/2 .
However, it is still a rough bound, since this only leads to 1/2 = O(∆t 1/2 ), when ∆t goes to zero, which is not sufficient to prove the convergence. An accurate estimation of the Gauss quadrature errors is needed to improve the bound, which should be investigated further in future works.

It seems now clear that the time evolution of the fractional moments depends on other moments than the ones involved in the moment system [START_REF] Frenklach | Method of moments with interpolative closure[END_REF] and we cannot restrict the quadrature to the mere moments of positive order, which also corresponds to taking n - q = 0, if we want to be able to provide a good approximation of the moment evolution. This is thus the key issue with fractional moment evolution, since we need at least to consider some moments of negative order in the quadrature, i.e. n - q > 0. After a series of test-cases, we would be tempted to conjecture that for n - q = 1, where two supplementary moments of negative order (m

[K∆t,1] -1/2 and m [K∆t,1]
-2/2 ) are also represented by the quadrature, the solution approximates accurately the exact kinetic solution. We were not able to complete the proof of such a conjecture and some of these results are presented in Section 6. It seems that the terms Ēl/2 where l = 1, 2 in equation ( 36) cause the divergence of the standard algorithm proposed in section 4.2. This is an important result, since we need only a total of three quadrature nodes n q = 3 to cancel the terms E l/2 for l = 1, 2 and to capture correctly the kinetic evolution. Besides, we can use more quadrature points to increase the precision by choosing n - q ≥ 2. Let us underline that the proper approximation and closure of the negative moments is here a key issue.

New algorithm:. according to the last results, we propose a 4-steps algorithm. This algorithm is named NEMO (Negative Moments) algorithm and described below: B. We calculate the negative order moments in the interval [K∆t, 1]

(38) m [K∆t,1] -a/2 = 1 K∆t s -a/2 n M E (s)ds,
for a = 1, . . . , 2n - q , where 2n - q ≥ 2 is the number of additional moments of negative order used in this algorithm and chosen by the user. The other moments of positive order are computed using the disappearance flux: [START_REF] Pollack | Bivariate extensions of the extended quadrature method of moments (eqmom) to describe coupled droplet evaporation and heat-up[END_REF] m

[K∆t,1] k/2 = m k/2 (0) -Φ -,k/2 (t n ), k = 0, . . . , 3,
where Φ is the disappearance flux given in [START_REF] Kapur | Entropy optimization principles with applications[END_REF].

C. The abscissas S j = r 2 j ∈ [K∆t, 1] and the weights w j = w j r 2n - q j of the quadrature corresponding to the moments m

[K∆t,1] p/2
for p = -2n - q , . . . , 3 are computed using the Product-Difference Algorithm [START_REF] Gordon | Error bounds in equilibrium statistical mechanics[END_REF] corresponding to the moments

m k = m [K∆t,1] k/2+n - q for k ∈ {0, . . . , 2n q -1}, leading to (40) m [K∆t,1] p/2 = nq j=1 w j S p/2 j , p = -2n - q , . . . , 3
D. Finally, we calculate the updated moments as follows:

(41) m k/2 (t n + ∆t) = nq j=1 w j (S j -K∆t) k/2 .
The singularities of the negative moment-integrals, when ∆t is very small, limits the use of high values of n - q . But in practice K∆t > 1.e -4 and we will show that the choice of n - q = 1 or n - q = 2 are sufficient to obtain an accurate solution. In these cases, the integral computation of the negative order moments can be achieved correctly with 24 Gauss-Legendre quadrature points. For more complex evaporation laws, the algorithm can be straightforwardly generalized by computing the Lagrangian evolution of the abscissas. In other words, equation [START_REF] Saurel | Diffuse-interface capturing methods for compressible two-phase flows[END_REF] becomes: [START_REF] Tagliani | Hausdorff moment problem and maximum entropy: a unified approach[END_REF] m k/2 (t n + ∆t) = nq j=1 w j S(t n + ∆t; t n , S j ) k/2 , and the integral of negative moments in [START_REF] Pickett | Engine combustion network[END_REF] becomes:

(43) m [K∆t,1] -a/2 = 1 S(tn;tn+∆t,0)
s -a/2 n M E (s)ds.

5. Numerical resolution of the moment governing equations. In this section, we present the global numerical resolution of system of equations [START_REF] Frenklach | Method of moments with interpolative closure[END_REF]. To resolve this system, we adopt operator splitting techniques [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF][START_REF] Descombes | Analysis of operator splitting in the non-asymptotic regime for nonlinear reaction-diffusion equations. application to the dynamics of premixed flames[END_REF] to separate the resolution of the convective transport part and the source terms. For the transport part, we extend the kinetic finite volume scheme [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF] used for the resolution of the EMSM model to the present case of fractional moments, while ensuring robust and realizable resolution. The obtained kinetic finite volume scheme is detailed in D. The source terms consists of two parts: evaporation and drag force. In this paragraph, we present a coupled solver for the spray evolution under evaporation and drag force. The corresponding system of equations reads:

(44)                    d t m 0/2 = -Kn(S = 0), d t m 1/2 = - K 2 m -1/2 , d t m 2/2 = -Km 0/2 , d t m 3/2 = - 3K 2 m 1/2 , d t m 2/2 U = -Km 0/2 U + m 0/2 U g -U θ .
Since the first four equations do not depend on the last one, the updated moments can be computed using NEMO algorithm. For the last equation, we use the method developed in [START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays: modeling and numerical issues[END_REF] to solve the velocity evolution due to the drag force and evaporation. The momentum evolution is conducted in two steps: first, we remove the part of the droplets, which will completely evaporate during the time interval [t n , t n+1 ], by evaluating the disappearance fluxes of the moments and momentum:

(45) U [K∆t,1] = U - Φ - Φ -,2/2 U ,
where U = (m 0/2 , . . . , m 3/2 , m 2/2 U t ) t and Φ -is the disappearance flux vector of the moments [START_REF] Kapur | Entropy optimization principles with applications[END_REF]. Then, the quadrature abscissas S i (t n ) and weights w i (t n ) are computed from [START_REF] Sabat | Eulerian modeling and numerical methods for the description of turbulent polydisperse sprays[END_REF] from the moments on the support [K∆t, 1]. The computation of the moments is achieved with [START_REF] Saurel | Diffuse-interface capturing methods for compressible two-phase flows[END_REF]. To update the momentum in this last step, we consider a size-velocity correlation: at time t n , we attribute the velocity c i (t n ) = U (t n ) to each abscissa S i (t n ). Then, the evolution on the time step [t n , t n+1 ] of c i and S i are computed by solving the following system:

(46)

dc i dt = - U g -c i θS i , dS i dt = -K.
The final momentum is computed at t = t n+1 as follows:

(47) (m 2/2 U )(t n+1 ) = nq i=1 w i S i (t n+1 )c i (t n+1 ).
The method can be generalized to more complex evaporation law by replacing -K in equation ( 46) by a general evaporation rate R S (S).

6. Numerical results. This section is dedicated to some representative testcases and analysis of numerical results, to verify the robustness and the accuracy of the proposed numerical schemes. In the first part, we test the new evaporation algorithm in the case of d 2 law with R S (S) = -K = -1 for two different initial conditions. A complementary evaporation test in the case of R S (S) = -(1 + 0.5S) is also presented to test the algorithm accuracy with a non constant evaporation rate. The second part focuses on two cases of transport in 1D. First, we test the accuracy of the kinetic finite volume schemes dedicated to the transport. Second, a critical case of a δ-shock is performed to evaluate the robustness of these schemes. Finally, a 2D case of an evaporating spray in the presence of a steady gaseous flow field, given by Taylor-Green vortices, is presented, in order to qualify the robustness and accuracy of the method compared to EMSM with its original numerical schemes [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF]. In this section, we compare the solution obtained with NEMO algorithm (n - q = 1), using the time step ∆t = 0.002, with the exact kinetic solution. One can see from Figure 2, where we plot the reconstructed NDF of NEMO solution and the exact kinetic distribution at t = 0.1 and t = 0.2, that the ME-reconstructed distribution follows accurately the exact distribution. It is important to notice that we are able to capture an accurate solution of the size distribution by using only four moments. To go further in quantitative comparison, Figure 3 shows the evolution of the relative error of the four fractional moments. The moment errors relatively to their initial values do not exceed 0.2%. NEMO algorithm then allows to obtain an accurate solution of the moments by using only three quadrature points [START_REF] Saurel | Diffuse-interface capturing methods for compressible two-phase flows[END_REF]. This avoids Fig. 5: Solutions of ME-reconstructed NDF at t = 0.3 (left) and t = 0.6 (right) using: NEMO with n - q = 1 (cross), exact kinetic solution (solid line) and the initial distribution (dashed line).

to compute directly the integral in ( 28) by using a large number of Gauss-Legendre quadrature points to approximate the moment-integrals.

6.1.2. d 2 law evaporation for a discontinuous initial NDF. In this second case, we test the new algorithm NEMO in the case of a discontinuous initial NDF:

(49) n 0 (S) = 1, if S ∈ [0.1, 0.6] 0, otherwise
The initial NDF defined in [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF] and the initial ME-reconstructed NDF are plotted in Figure 4. In Figure 5, we present the reconstructed NDF computed using NEMO algorithm with n - q = 1 at two different times. In the same figures, we compare the obtained results with the exact kinetic solution 3 . As in the previous case, NEMO algorithm shows an accurate prediction of the exact kinetic solution. Furthermore, Figure 6-left shows that the relative error of the moments are less than 1.%, in the case of n - q = 1 and ∆t = 6.e -3, which is an accurate result. We can improve the accuracy of NEMO algorithm by using a smaller time step ∆t = 6.e -4 or by using n - q = 2 and keeping the same time step as the previous case. In figure 6-right, we present the relative moment error using NEMO algorithm with n - q = 2 and ∆t = 6.e -3. In this case, the relative moment errors are less than 0.3%. 6.2. Accuracy of NEMO algorithm for non constant evaporation laws. NEMO scheme has been developed under the assumption of a d 2 law, but as it was explained before, the algorithm can be generalized for more complex laws by solving the Lagrangian equation [START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF] for each abscissas S j given in the third step of the algorithm. In this section, we propose to evaluate the accuracy of the algorithm in the case where the evaporation rate depends linearly on the size:

(50) R s (S) = -(a + bS).
The exact kinetic solution can be computed according to equation [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF]. In the following, we set a = 0.5 and b = 1. Figure 7 presents the ME-reconstructed NDF computed by NEMO algorithm, and the exact kinetic solution of the NDF at t = 0.3 Fig. 6: Evolution of the relative moment errors using NEMO with n - q = 1 (left) and n - q = 2 (right): m 0 (solid line), m 1/2 (dash-dotted line), m 1 (cross) and m 3/2 (circle). Fig. 7: The evolution of the NDF in the case of a linear evaporation rate: initial ME reconstructed solution (dashed line), NEMO algorithm using n - q = 1 (cross), and exact kinetic solution (solid line), at times t = 0.3 and t = 0.6. and t = 0.6. The relative errors of the moments are given in Figure 8. We can see from these results, the accuracy of the generalized NEMO algorithm to predict the kinetic evolution.

Transport in 1D simulation.

In this part, we investigate the accuracy and robustness of the transport schemes, the first order and the second order detailed in the Appendix D. The evaporation and the drag force are not considered in this study. First, we investigate the order of accuracy of the transport schemes. Then, we test the two numerical schemes in a critical case, where we generate a δ-shock singularity.

6.3.1. Accuracy order study. We consider an initial size distribution with a form depending on the coordinate space x. In fact, we need a non trivial initial sizespace distribution, such that the slopes used in the second order scheme do not vanish all the time. The chosen initial size distribution has the following profile:

(51) n(t, S, x) = 10 exp - (x -x c ) 2 σ 2 x exp - ( √ S -(1 -x)/2) 2 σ R ,
where x c = 0.25, σ x = 0.1 and σ R = 0.3. The initial velocity field is initiated as follows:

(52) u(t = 0, x) = 0.5 -x x < 0.5, u(t = 0, x) = 0.

x ≥ 0.5.

Figure 9 shows the spatial distribution of the moment m 0 obtained with the first and second order schemes at t = 0.5. Comparing the two schemes, it is clear that the results obtained with the second order scheme are more accurate than the ones obtained with the first order. In Figure 10, we show the grid convergence for the first and second order to the analytic solution at t = 0.8. The second order scheme converges faster than the first order scheme to the solution. In the extremum of the moment spatial profile, we can see that the convergence to solution is slow compared to the other points. In fact, the slope limitation in this point is activated to ensure a non oscillating solution, but, it introduces some numerical diffusion.

We also compute the L 1 -error for each numerical scheme depending on the grid size. In Figure 11, we display the order curve for the two schemes. As it can be seen from these curves, the order of the first order scheme is around 0.6 and of the second order, around 1.5. We mention that the order was computed in a critical case. Indeed, the characteristic curves cross at t = 1.0 and the order of accuracy was evaluated at t = 0.8. At this time the solution starts to have high gradient variations, thus the observed decrease of the order of accuracy. 6.3.2. Robustness and capacity of capturing δ-shocks. The monokinetic assumption is not a valid hypothesis when droplets cross. In such an event, the monokinetic model generates a δ-shock. Despite the non-physical solution, the two kinetic schemes should be able to provide a solution in all critical situations and even resolve the created singularities. In this section, we test the robustness of the numerical schemes in a case of a strong crossing. Initially, we consider two spatial Gaussian distributions centered at two positions and traveling at opposite velocity as shown in Figure 12. Figure 13 shows the spatial profile of the moment m 0 at two different times. We can see the shock generation when the two packets cross. At the end, the droplets accumulate at the center and the corresponding spatial distribution has the form of a Dirac delta measure. The computation is still robust in this case and realizability preserved.

6.4. 2D simulation: Transport, evaporation and the drag force. After the model and numerical schemes have been tested in 0D and 1D simulations, we propose in this part to compare in a classical 2D configuration the present model and numerical schemes with the EMSM model and its original numerical schemes [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF], in a case where we consider transport, evaporation and drag. The simulations are performed using the CanoP code, developed within the collaboration of Maison de la Simulation, IFPEn and EM2C Laboratory. It is based on the p4est library [START_REF] Burstedde | p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF], a library providing Adaptive Mesh Refinement capability highly scalable in massively parallel computations [START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF][START_REF] Essadki | High order moment methods and adaptive mesh refinement for polydisperse evaporating spray simulations[END_REF]. In the present simulation, we use only a uniform grid, since we are not concerned with Adaptive Mesh Refinement in the present work. We consider an evaporating spray in the presence of Taylor-Green vortices for the gas, which is a steady solution of the inviscid incompressible Euler equations. The Fig. 13: The spatial profile of m 0 in the crossing case at t = 0.6 (left) and t = 1.2 (right), using 128 cells: first order scheme (cross) and second order scheme (circle). non-dimensional velocity field of the gas reads: We have used second order scheme for the transport resolution in both cases (see [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF] for EMSM and D for the fractional moment model) and the original evaporation algorithm [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] to solve evaporation in the EMSM model and NEMO algorithm with n - q = 1 for fractional moments. For the EMSM model, the volume fraction is not resolved but calculated through ME reconstruction of the size distribution (1/6 √ π)

1 0 S 3/2 n M E (S)dS. Instead, for the new model, the volume fraction is directly calculated as (1/6 √ π)m 3/2 . The results of the two computations are closely similar, and the L 1 -norm difference relatively to the initial volume fraction field is less than 3% at t = 1.

Figure 16 provides a quantitative comparison, where we plot the volume fraction along the segments: (x = 0.3, y ∈ [0, 0.5]) and (x ∈ [0, 0.5], y = 0.3) at t = 0.5. Let us underline that the EMSM model in this configuration was validated through a detailed comparison with the Multi-fluid model in [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF] and shown to be able to predict the evaporation and the mean dynamics of the spray4 . Besides, the computational time was shown to be much shorter for the EMSM model, by a factor almost related to the ratio of the number of PDE solved in the model 5 .

Thus we can conclude on two points: 1-the comparison of the two models validates the results of the new model and assesses its ability to capture the dynamics and Even if we do not claim this number is generic, the new model is however going to be even cheaper in terms of computational time. Let us also emphasize that in our numerical strategy, even if a large amount of time can be spent in evaluating the source terms (including maximization of entropy) and evolution in phase space, operator splitting yields an embarrassingly parallel problem for which work stealing and task-based optimization can be extremely efficient.

Conclusions.

In the present paper, we have proposed a novel Eulerian high order size moment model and numerical methods in order to describe polydisperse evaporating sprays, with the same accuracy and robustness based on built-in realizability preserving numerical schemes as previously developed methods. Its main originality in terms of modeling is to involve a set of variables, which can be interpreted physically as interface geometry quantities and are thus easier to couple to other zones of the flow, where the spray is denser or even in the separated-phases region. Since the new model involves fractional order moments, we have tackled the various mathematical difficulties related to ME reconstruction and moment space theory [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF][START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] and proposed a new precise and robust numerical method, within the framework of operator splitting techniques [START_REF] Doisneau | Numerical strategy for unsteady two-way coupled polydisperse sprays: application to solid-rocket instabilities[END_REF][START_REF] Descombes | Analysis of operator splitting in the non-asymptotic regime for nonlinear reaction-diffusion equations. application to the dynamics of premixed flames[END_REF] to separate the resolution of transport in physical space and phase space, for a model involving high order fractional size-moments as well as negative order moments. The way of treating negative order moments properly should be useful for various other problems in the several communities [START_REF] Frenklach | Method of moments with interpolative closure[END_REF]. The main difficulty is to design accurate, robust and realizable numerical schemes, especially for the resolution of the evaporation and transport, where standard solvers fail.

Whereas the basis of the approach has been presented in a simplified modeling framework in order to focus on its applied mathematical aspects and for the sake of clarity, the extension to realistic models and configurations can be achieved. Detailed droplet models can be used in order to cope with predictive simulations, as well as two-way coupling and mesh motion, in internal combustion engines for IF-PEn [START_REF] Kah | High order moment method for polydisperse evaporating spray with mesh movement: application to internal combustion engines[END_REF]. Higher order velocity moment methods, conditioned on droplet size and based on the Levermore hierarchy of hyperbolic system of conservation laws with an entropic structure [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Levermore | The Gaussian moment closure for gas dynamics[END_REF][START_REF] Sabat | Eulerian modeling and numerical methods for the description of turbulent polydisperse sprays[END_REF] allow to cope with turbulence and particle trajectory crossing. Implementation in the massively parallel code able to handle rather complex geometries has been accomplished [START_REF] Drui | Eulerian modeling and simulations of separated and disperse two-phase flows : development of a unified modeling approach and associated numerical methods for highly parallel computations[END_REF] and the use of the present model assessed in a companion paper [START_REF] Essadki | High order moment methods and adaptive mesh refinement for polydisperse evaporating spray simulations[END_REF], which rather focuses on implementation issues, parallelism and results. It shows the potential of the present model for realistic engine simulations.

However, the main challenge still remains the elaboration and design of a unified model for the liquid injection able to predict atomization as well as the generated polydisperse spray and its dynamics and evaporation [START_REF] Essadki | Statistical modeling of the gas-liquid interface using geometrical variables: toward a unified description of the disperse and separated phase flows[END_REF], for which the present contribution offers some key ingredients.

for the existence and uniqueness of the Maximum Entropy reconstruction from any moment vector in the interior of the fractional moment space, as well as an algorithm to compute the parameters of this reconstruction. Finally, using the canonical moments, the extension of the realizable transport scheme developed for integer moments in [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF] is provided.

The existence and uniqueness of the ME solution is first shown, before giving an algorithm to compute the parameters of such a reconstruction. C.1. Existence and uniqueness of the maximum entropy solution. The following Theorem is shown here.

Theorem 6. If the vector m N = (m 0 , m 1/2 , . . . , m (N )/2 ) belongs to the interior of the Nth fractional moment space, then the constrained optimization (63) problem admits a unique solution, which is in the following form:

(64) n M E (S) = exp -λ 0 - N i=1 λ i S i/2 .
We mention that the case of the integer moments has been already treated in [START_REF] Mead | Maximum entropy in the problem of moments[END_REF].

Here, some ideas of this work are used, but the present proof is completely different and simplified. Indeed, Mead and Papanicolaou [START_REF] Mead | Maximum entropy in the problem of moments[END_REF] have used the monotonic properties of the moments, which is a characterization of the integer moment space, to prove existence of the ME solution. In our case, only the definition of the fractional moment space is used. The proof is then done through two Lemmas.

C.1.1. Uniqueness. The first Lemma allows to show the uniqueness and to give the form of the ME reconstruction.

Lemma 7. If the constrained optimization problem (63) admits a solution, then this solution is unique and can be written in the following form:

(65) n M E (S) = exp -λ 0 - N i=1 λ i S i/2 ,
where λ = (λ 0 , .., λ N ) t ∈ R N . Moreover, this problem is equivalent to find an optimum of the potential function:

(66) G(λ 0 , . . . , λ N ) = 1 0 exp(-λ 0 - N i=1 λ i S i/2 )dS + N -1 k=0 λ k m k/2 .
Proof. The Lagrangian function associated to this standard constrained optimization problem is:

(67) L(n, λ) = H[n] -(λ 0 -1) 1 0 n(s)ds -m 0 - N i=1 λ i 1 0 s i/2 n(s)ds -m i/2 ,
where λ = (λ 0 , .., λ N ) is the vector of the Lagrange's multipliers. Let us suppose that, for a given moment vector m N , there exists a density function n M E which is the solution of the ME problem (63). So, there exists a vector λ M E for which the differential of the Lagrange function L(n, λ) at the point (n M E , λ M E ) vanishes:

(68)    DL(n M E , λ) • (h, 0) = 1 0 h(s) -ln(n M E (s)) - N i=0 λ i s i/2 ds = 0, ∂L ∂λi (n M E , λ) = 1 0 s i/2 n M E (s)ds -m i/2 = 0, ,
We choose to present the scheme in a two dimensional space to lighten the notations and we note U = (u, v). For the transport resolution in physical space, we use a dimensional splitting algorithm. In this context we consider a free transport in one direction (we present the x-direction here) of the droplets without the evaporation nor the drag force.

(78)

∂ t m k/2 + ∂ x (m k/2 u)= 0, k = 0, . . . , 3, ∂ t (m 1 U ) + ∂ x (m 1 u U )= 0,
The mathematical structure of the pressureless gas system leads to some singularities (known as δ-shocks). These singularities occur when the monokinetic assumption is violated. This can happen when trajectory crossings take place and lead to particles accumulation in a very small volume. Following the idea of de Chaisemartin [START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF], Kah et al [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF] developed a finite volume kinetic scheme for the EMSM model. We use the same approach to solve numerically the system (78). In the following, we present briefly the main steps to derive the kinetic scheme for the system (78):

A. We write the equivalent kinetic system to the pressureless system (78), as it was proposed in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF]:

(79) ∂ t f + ∂ x (c x f ) = 0 , and f (t, x, c, S) = n(t, x, S)δ(c -U ), B. We use the finite volume discretization of the system (78): ( 80)

m n+1 i = m n i - ∆t ∆x (F i+1/2 -F i-1/2 ), p n+1 i = p n i - ∆t ∆x (G i+1/2 -G i-1/2 ),
where the fluxes are expressed as function of the NDF:

F ± i+1/2 G ± i+1/2 = 1 ∆t tn+1 tn 1 0 R         1 S 1/2 S S 3/2 Su Sv         f (t, x i+1/2 , c, S)dcdSdt,
C. We split the fluxes in two integral parts: F i+1/2 = F + i+1/2 + F - i+1/2 and G i+1/2 = G + i+1/2 + G - i+1/2 . The first (resp. the second) corresponds to the droplet of positive (resp. negative) velocity in x-direction. Then, we use the exact solution of the kinetic system (79), to express the fluxes as function of the NDF at t = t n . D. Finally, the fluxes are expressed as function of a spatial reconstruction of the moments and velocities at t = t n :

(81)

F ± i+1/2 G ± i+1/2 = 1 ∆t x i+1/2 x i-1/2         m 0 (t n , x) m 1/2 (t n , x) m 1 (t n , x) m 3/2 (t n , x) m 1 u(t n , x) m 1 v(t n , x)         1 Σ ± (x)dx,
where Σ ± = x , ±(x i+1/2 -∆tu(t n , x )) < ±x D.2. Second order scheme. Kah et al [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF] developed a realizable second order kinetic scheme. They showed that the canonical moments (in the context of integer moments) are transported variables. Therefore, these quantities satisfy a maximum principle. Since, the canonical moments live in the simple space [0, 1] N (N = 3 in our case). The authors proposed to use linear reconstruction of the canonical moments to design a high order scheme, instead of reconstructing directly the moments, which belong to a complex space. We adopt the same approach with some adaptations for the fractional moments. After reconstructing the variables (velocity and canonical moments) the fluxes are computed by a simple integration. D.2.1. Reconstruction. The reconstructed variables are the moment m 0 , the canonical moments (defined in part B) p 1 , p 2 , p 3 and the velocity. x i-1/2 m 0 (x)p 1 (x)[(1 -p 1 )p 2 + p 1 ](x)u(x)dx.

Compared to the expressions developed in the case of the EMSM model, only the last integral expression is different. In fact, the velocity is weighted with the moment m 1 for both models, but in the new model, m 1 acts as a second order moment. For this reason, the expression of the moment m 1 as function of the canonical moment is different from the one in the case of integer moments.

Kah el al. [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF] show that the bar terms can be written as follows:

(85) p k,i = a k,i + b k,i Dp k,i , u i = a u,i + b u,i Du i , where for each k, a k,i and b k,i are independent of Dp k,i , and a u,i and b u,i are independent of Du i . D.2.2. Slope limitation. In order to satisfy the maximum principle for the transported quantities (the canonical moment and the velocity) and the positivity of the number density m 0 , the slopes should be calculated carefully. Following the development done in [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF], the slopes are calculated as follows: (86) Dm 0,i = φ(m n 0,i-1 , m n 0,i , m n 0,i ) min Using the equations ( 84),( 85) and (86) the slopes and the bar variables can be expressed as functions of the current and neighbor cell variables. However, these algebra relations are quite heavy. Therefore, their calculation is achieved using Maple software. D.2.3. Fluxes Computation. After computing the slopes and the bar variables, the fluxes can be computed as follows: (87) .

|m n 0,i+1 -m n 0,i | ∆x , |m n 0,i -m n 0,i-1 | ∆x , 2m n
F + i+1/2 G + i+1/2 = 1 ∆t x i+1/2 x L i+1/2 m 0         1 p 1 p 1 [(
F - i+1/2 G - i+1/2 = - 1 ∆t x R i+1/2 x i+1/2 m 0         1 p 1 p 1 [(1 -p 1 )p 2 + p 1 ] p 1 (1 -p 1 )( (88)
The expressions inside the integrals are polynomial functions of x of order up to 6, their calculation can be achieved by using four points of the Gauss-Legendre quadrature.

Fig. 1 :

 1 Fig. 1: Initial size distribution (dashed line) and the reconstructed size distribution at t = 0.1: kinetic solution (solid line), original evaporation scheme for fractional moments (cross).

[

  K∆t,1] l where l = 0/2, 1/2 . . . 3/2) and on an infinite set of the moments of support [K∆t, 1] (m [K∆t,1] l

Fig. 2 :

 2 Fig. 2: ME-reconstructed NDF solution of NEMO with n - q = 1 (cross), exact kinetic NDF (solid line) and the initial distribution (dashed line), at time t = 0.1 (left) and t = 0.2 (right).

6. 1 .

 1 Evaporation in 0D simulation.6.1.1. Evaporation with d 2 law for an initial smooth NDF. We consider an initial NDF in the form of the ME-reconstruction NDF, which is the same initial condition as the one used in Figure1:[START_REF] Williams | Spray combustion and atomization[END_REF] n 0 (S) = exp(-16(S 1/2 -1/4) 2 (S 1/2 + 1)).

Fig. 3 :

 3 Fig. 3: Relative moment errors, NEMO (n - q = 1): m 0 (solid line), m 1/2 (dashdotted line), m 1 (cross) and m 3/2 (circle).

Fig. 4 :

 4 Fig. 4: The ME-reconstructed NDF in the solid line and the initial discontinuous NDF in the dashed line.

Fig. 8 :

 8 Fig. 8: Evolution of the moment errors relatively to their initial value, NEMO algorithm: m 0 (solid line), m 1/2 (dashdotted line), m 1 (cross) and m 3/2 (circle).

Fig. 9 :

 9 Fig.9: Spatial profile of the moment m 0 at t = 0.5 with 128 grid cells. First order scheme (cross), second order scheme (circle) and exact solution (solid-line).

Fig. 10 :

 10 Fig. 10: Grid convergence; first order: 64 cells (cross) and 256 cells (plus), second order: 64 cells (circle) and 256 cells (square). Exact solution in blue dashed-line.

Fig. 11 :

 11 Fig. 11: Error curves of m 0 with respect to grid refinement in logarithm scale: first order (cross) and second order (circle) schemes.

Fig. 12 :

 12 Fig. 12: Initial condition for the crossing case. Left: Initial moment fields, the curves represent the moment with decreasing order in terms of value. Right: initial velocity.

  (53) u g (x, y) = sin(2πx) cos(2πy), v g (x, y) = -cos(2πx) sin(2πy), where (x, y) ∈ [0, 1] 2 and with periodic boundary conditions. Initially, the spray is localized in the bottom-left vortex. The initial spatial size-distribution is given by:(54) n(t, x, S) = 1 [a,b] (S)1 {x ,||x -xc||2< √ 2r} (x) exp(-||x -x c || 2 2 /r 2 ),where [a, b] = [0.25, 0.75], x c = (0.15, 0.15) and r = 0.1. The initial Stokes number computed according to the mean size S = m 1 /m 0 is St( S) = 0.05, which is close to the critical Stokes number. The Stokes number decreases over time because of evaporation. The spray evaporation rate is K = 0.5. Figures 14-15 present the computed spatial distribution of the volume fraction at two different times, using the EMSM model (left) and the new fractional moment model (right).

Fig. 14 :

 14 Fig.14:The spatial distribution of the volume fraction for the Taylor-Green simulation at t = 0.5. The computation is carried out in a uniform grid 128 × 128.

Fig. 15 :

 15 Fig. 15: The spatial distribution of the volume fraction for the Taylor-Green simulation at t = 1.0. The computation is carried out in a uniform grid 128 × 128.

Fig. 16 :

 16 Fig. 16: Volume fraction profiles at t = 0.5 with EMSM (dashed red line) and fractional moment model (solid black line) over the segments: x ∈ [0., 0.5] and y = 0.3 (left) and x = 0.3 and y ∈ [0., 0.5].

2 x 2 x i- 1 /2 m 0 p 1 ( 1 -p 1 )( 1 -p 2 )p 2 p 3 + [( 1 -p 1 )p 2 + p 1 ] 2 (

 22111112311212 x) = m 0,i + Dm 0,i (x -x i ), p 1 (x) = p 1,i + Dp 1,i (x -x i ), p 2 (x) = p 2,i + Dp 2,i (x -x i ), p 3 (x) = p 3,i + Dp 3,i (x -x i ), u(x) = u i + Du i (x -x i ), v(x) = v i + Dv i (x -x i ), where x ∈ [x i-1/2 , x i+1/2]. Generally the quantities with the bar are different from the cell averaged quantities p k,i , u k,i and v k,i and they are determined depending on the slopes and the following conservation properties: i-1/2 m 0 (x)p 1 (x)[(1 -p 1 )p 2 + p 1 ](x)dx,

  φ(a, b, c) = 1/2(sgn(b -a) + sgn(c -b)).

  1 -p 1 )p 2 + p 1 ] p 1 (1 -p 1 )(1 -p 2 )p 2 p 3 + [(1 -p 1 )p 2 + p 1 ] 2 p 1 [(1 -p 1 )p 2 + p 1 ]u p 1 [(1 -p 1 )p 2 + p 1 ]v

	
	      	dx,
	and

  1 -p 2 )p 2 p 3 + [(1 -p 1 )p 2 + p 1 ] 2 p 1 [(1 -p 1 )p 2 + p 1 ]u p 1 [(1 -p 1 )p 2 + p 1 ]v

			
			      	dx,
	such that
	(89)	x L i+1/2 = x i+1/2 -∆t x R i+1/2 = x i+1/2 -∆t	(ū i + ∆x 2 Du i ) + 1 + ∆tDu i (ū i+1 -∆x 2 Du i+1 ) + , 1 + ∆tDu i+1

Entropy maximization with fractional moments is used in[37,

[START_REF] Gzyl | Hausdorff moment problem and fractional moments[END_REF][START_REF] Tagliani | Hausdorff moment problem and fractional moments: A simplified procedure[END_REF] as a remedy to illconditioning related to a high number of integer moments[START_REF] Talenti | Recovering a function from a finite number of moments[END_REF]. The considered set of fractional moments is then recovered from the integer ones, and their orders are optimized to minimize the entropy difference with the real function. In our contribution, a known and small set of fractional moments is used, deduced from physical considerations, so that the problematic is different.

We can also use hybrid approach, in this case we take into account the evaporation fluxes coming from the right section: m(tn) -Φ (i) -(tn) + Φ (i+1) -(tn), where i indexes the size-section.

The exact kinetic solution is computed for the initial ME reconstructed size distribution.

EMSM model and present model, with one size-section, are limited in predicting the sizeconditioned dynamics compared with Multi-fluid model with a Stokes number computed from a mean size. The Correlated Size Velocity Moment model[START_REF] Vié | Size-velocity correlations in high order moment methods for polydisperse evaporating sprays: modeling and numerical issues[END_REF] allows to tackle this issue.

[START_REF] De Chaisemartin | Eulerian models for turbulent spray combustion with polydispersity and droplet crossing[END_REF] Switching from 30 variables -multi-fluid

-to 6 variables -EMSM -in 2D resulted in a acceleration factor of more than 4, the maximization of entropy being responsible for not reaching a factor of 5.

Acknowledgments. This research was supported by a PhD grant for M. Essadki from IFPEn, as well as by Ecole Doctorale de Mathématiques Hadamard and EM2C. We want to thank A. Larat and A. Vié for several interesting and helpful discussions.

Appendix A. Introduction. In the main paper, a set of fractional moments are used to describe polydisperse evaporating sprays and correspond to physically identified interface geometrical quantities. In this Appendix, we present some mathematical elements related to the use of these fractional moments. These elements are already known for integer moments and can be extended in our case. First, the fractional moment space is defined and the notion of canonical moments and quadrature is directly extended from the integer moments case. Then, a complete proof is given Appendix B. Fractional moment space, canonical moments and quadrature. The purpose of this part is to extend some definition and relevant properties of the integer moment space to fractional moment space, in order to provide necessary tools to characterize the geometry of this space and thus to design realizable numerical schemes. In the following, we use the normalized fractional moments c k/2 = m k/2 /m 0 .

B.1. Fractional moment space -link with the integer moment space. Definition 5. We define the Nth fractional moment space as follows:

where P([0, 1]) denotes the set of all probability density measures defined on the interval [0, 1] and c k/2 = 1 0 x k/2 dµ(x). The fractional moments can be expressed as integer moments by using the following mapping r 2 = x for x ∈ [0, 1]. Let us then denote µ the measure on the interval [0, 1], which is in a one to one correspondence with the measure µ on [0, 1] through the above transformation, that is: µ([0, r]) = µ([0, r 2 ]). The moments can then be written:

This relation expresses an identification between the fractional moment c k/2 (µ) of the measure µ and the integer moment c k = c k ( µ) of the measure µ. In the following, we use these notations to differentiate between the two natures of the moments, even if they are equal. With this simple identification, we will take benefit from the already existing results of integer moment space to extend them to the case of fractional moment space.

B.2. Canonical moments. For a fractional moment vector c N ∈ M 1/2

N ([0, 1]), we denote by P N (c N ) the set of all measure µ ∈ P([0, 1]), which are the solution of the following finite Hausdorff moment problem:

), we can show that the set P N (c N ) is infinite. Indeed, these results were shown in [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] for integer moments, and its generalization for fractional is straightforward through the identification (55). Furthermore, the set of c (N +1)/2 (µ), where µ ∈ P(c N ), is infinite and the maximum and the minimum of this set are different:

We define the canonical fractional moments as follows:

(58)

.

The canonical moment vector belongs to [0, 1] N , which has a simple geometry compared to the moment space. By using the results obtained on canonical integer moments [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF], we can write the algebraic relation between fractional moments and their corresponding canonical moments, in the case of N = 3, by using the identification (55):

.

belongs to the boundary of the moment space M 1/2 N +1 , and the measure set P N +1 (c ± N +1 ) = {µ ± } has only one element. The measure µ + (resp. µ -) is called the upper (resp. lower) principal representation. In the case of integer moment vector c N , it was shown [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] that the lower and upper principal representations ( µ + and µ -) can be expressed as sum of n s ≤ (N + 1)/2 weighted delta-Dirac functions (we count the abscissas in the interior ]0, 1[ by 1 and the ones in the extremity, 0 or 1, by 1/2):

We recall that the subscript • is used for integer moments and their corresponding measure, which are related to the fractional moments according to the identification (55).

In the case where N is odd (corresponding to an even number of moments), the lower principal representation corresponds to the classical Gauss quadrature. Let us remove the superscripts -to simplify the notations. Then, the weights w i and the abscissas r i (which are in ]0,1[) can be computed from the moments c N , for example with the Product-Difference algorithm [START_REF] Gordon | Error bounds in equilibrium statistical mechanics[END_REF]. In other words, this algorithm (efficiently) solves the following non-linear system:

The identification (55) allows us to relate the quadrature in terms of fractional moments, to the quadrature in terms of integer moments:

where x i = (r i ) 2 and w i = wi . Thus, the Product-Difference algorithm allows to determine the weights and abscissas of the quadrature for fractional moments.

Appendix C. Reconstruction through Entropy Maximization. The Maximum Entropy (ME) problem reads as follows:

where h is a positive distribution. Since the system (68) is valid for all h, it yields:

The problem then consists in finding a vector λ = (λ 0 , .., λ N ) in R N which satisfies the moment equations in the system (69). This problem is equivalent to find an optimum of the potential function G(λ 0 , .., λ N ) defined in (66).

The Hessian matrix H defined by H i,j = ∂ 2 G ∂λi∂λj is a positive definite matrix, which ensures uniqueness of an eventual existing solution.

C.1.2. Existence. The existence of such reconstruction is shown through the following Lemma.

Lemma 8. If the vector m N = (m 0 , m 1/2 , . . . , m (N )/2 ) belongs to the interior of the Nth fractional moment space, then the function G, defined in (66), is a continuous function in R N , and goes to infinity when ||λ|| → +∞.

Proof. Let us suppose that the last assertion is wrong, so there exists a sequence (λ

Hence, there exists A ∈ R such that:

We write for each n ∈ N, λ (n) = λ n (α

N ), such that

Since the sequence (α (n) ) n=0,1.. is a bounded sequence, we can extract a convergent subsequence (α φ(n) ) n , where φ : N → N is an increasing function and:

To simplify the notations, we can directly consider that α

Since the vector m N = (m 0 , m 1/2 , .., m (N )/2 ) is a moment vector, there exists a nonnegative distribution function f such that m k/2 = 1 0 s k/2 f (s)ds for k = 0, . . . N , and

Since the first integral is positive (73)

When n tends to infinity, we get:

(74)

We have Q = 0, f ≥ 0 and f = 0, and since Q is a continuous function, it follows from the inequality (74) that there exists [a, b] ⊂ [0, 1] in which Q(s) ≤ -B and B > 0. Since Q (n) converges uniformly to Q in [0, 1], then, for all s ∈ [a, b] and for the large enough values of n:

Using these results in the inequality (72), as well as α n k m k/2 ≥ -m 0 , we get:

In the limit when n goes to infinity, we get the contradiction +∞ ≤ A, thus concluding the proof.

C.2. Algorithm of the NDF reconstruction through the Entropy Maximization. The reconstruction of the NDF through the maximization of Shannon entropy goes back to finding the Lagrange's multipliers λ 0 . . . λ N such that:

Solving this nonlinear system is equivalent to optimize the convex function G(λ).

We solve the problem by using the Newton iterations as proposed in [START_REF] Mead | Maximum entropy in the problem of moments[END_REF]. The ME reconstruction is then done thanks to the Algorithm 1.

Algorithm 1 ME algorithm

Choose initial guess of the vector λ. The integral computations are done by using Gauss-Legendre quadrature. In [START_REF] Mead | Maximum entropy in the problem of moments[END_REF], it is shown that 24-point quadrature allows to accurately compute the different integral expressions involved in the Algorithm 1. Appendix D. Transport scheme. Here, first and second order schemes are given, which are directly extended from the schemes developed by Kah et al [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse spray[END_REF] for the transport of integer moments. For that, the canonical moments defined for fractional moments in section B are used.