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

Frédérique Laurent

Damien Kah

Stephane De Chaisemartin

S T Éphane De Chaisemartin

A robust moment method for evaluation

Keywords: Polydisperse sprays, aerosols, evaporation, Quadrature Method of Moments, Eulerian multi-fluid model, moment space, population balance equation, Maximum Entropy reconstruction, kinetic schemes AMS subject classifications. 35Q35, 65M12, 65M99, 76T10

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction and mathematical framework. The modeling and numerical simulation of polydisperse multiphase systems have become a crucial issue in a number of applications such as combustion, chemical engineering science (precipitation and crystallization, particle formation in aerosols, fluidization...) but also in biology and biochemistry (modeling of sprays in the lung...). In many industrial combustion applications such as Diesel engines, fuel is stocked in condensed form and burned as a dispersed liquid phase carried by a gaseous flow. Two phase effects as well as the polydisperse character of the droplet size distribution (since the droplets dynamics depend on their inertia and are conditioned by size) can significantly influence flame structure [START_REF] Massot | Counterflow spray diffusion flames of heptane: computations and experiments[END_REF]. The cloud of particles can also be submicronic particles such as soots produced in a flame which require some fine description of size distributions [START_REF] Zucca | Implementation of the population balance equation in cfd codes for modelling soot formation in turbulent flames[END_REF][START_REF] Hall | Predictions of soot dynamics in opposed jet diffusion flames, a tribute to irvin glassman[END_REF]. Size distribution effects are also encountered in a crucial way in solid propellant rocket boosters, where the cloud of alumina particles experiences coalescence and become polydisperse in size, thus determining the global dynamical behavior [START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF]. Besides there exists a large literature devoted to aggregation-breakage-growth phenomena in aerosol science and we refer to [START_REF] Ramkrishna | Population Balances: Theory and Applications to Particulate Systems in Engineering[END_REF] and [START_REF] Fox | Multiphase Reacting Flows: Modelling and Simulation[END_REF] for details and references. In all these cases, the fact that the cloud is polydisperse has to be modeled since it strongly influences the physics of the problem. Consequently, it is important to have reliable models and numerical methods in order to be able to describe precisely the physics of such flows.

The context is thus the following : a dispersed phase which can be liquid (cloud of droplets called spray or aerosol depending on their inertia) or solid (cloud of particles such as soots) interacts with a carrier phase. This dispersed phase is polydisperse, i.e. the particle size distribution (PSD) covers a large size spectrum. The modeling of such polydisperse multiphase systems relies most of the time on a mesoscopic description by a population balance equation (PBE) in chemical engineering science [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF][START_REF] Ramkrishna | Population Balances: Theory and Applications to Particulate Systems in Engineering[END_REF], also called a Williams equation in the spray modeling community [START_REF] Williams | Spray combustion and atomization[END_REF][START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], on the number density function (NDF). This approach describes the droplets/particles as a cloud of point particles which experience [START_REF] Abramov | A practical computational framework for the multidimensional momentconstrained maximum entropy principle[END_REF] transport in real space, (2) droplet evaporation/particle sintering or oxidation, molecular growth (3) acceleration of droplets due to drag when they have their own inertia conditioned by size, and (4) coalescencebreakup of droplets or aggregation-breakage of particles influencing polydispersity. The details of the interface behavior, angular momentum of droplets, detailed internal temperature distribution inside the droplet/particles, etc., are not predicted. Instead, a finite set of global properties such as for example size, velocity, temperature are modeled and called internal variables. This level of description is relevant since it is the only one for which numerical simulations at the scale of a combustion chamber or in a free jet can be conducted.

In this mesoscopic framework, there exists considerable interest in the development of numerical methods for simulating multiphase flows. In the literature, the Lagrangian Monte-Carlo approach [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF], called Direct Simulation Monte-Carlo method (DSMC) in [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF], is generally considered to be the most accurate for solving the PBE. However, its computational cost is high, especially in unsteady, spatially extended configurations coupled to a CFD code for the carrier phase. Besides the method is difficult to couple accurately with a Eulerian description of the carrier phase. There is thus considerable impetus to develop Eulerian methods based upon moment closure.

For the treatment of the polydispersity and the precise description of the PSD, there exists basically two options : 1-to discretize the size phase space into size intervals called "sections" since the work of Tambour and coworkers (see for example [START_REF] Greenberg | On the origin of spray sectional conservation equations[END_REF] and references therein). There exists a large variety of applications and methods such as multi-fluid models for sprays (see [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF][START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] extended from sectional models of Tambour, Greenberg et al. and [34,[START_REF] Massot | Eulerian multi-fluid models : modeling and numerical methods[END_REF][START_REF] De Chaisemartin | Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporation sprays[END_REF] for a review of recent results) or sectional approach for aerosols (see [START_REF] Vanni | Approximate population balance equations for aggregation-breakage processes[END_REF][START_REF] Kostoglou | Extended cell average technique for the solution of coagulation equation[END_REF] and references therein), also called class methods in [START_REF] Bove | A novel algorithm for solving population balance equations: the parallel parent and daughter classes. derivation, analysis and testing[END_REF]. 2-to adopt a moment point of view and to close the set of evolution equations through quadrature methods such as QMOM, quadrature method of moment initiated in [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF][START_REF] Wright | Bivariate extension of the quadrature method of moments for modeling simultaneous coagulation and sintering of particle populations[END_REF] or DQMOM, direct quadrature method of moment introduced by Marchisio and Fox (see [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF] and the recent review in [START_REF] Fox | Multiphase Reacting Flows: Modelling and Simulation[END_REF]). In a recent paper, Fox, Laurent and Massot have presented a comprehensive comparison between DQMOM and the Eulerian multi-fluid models for polydisperse evaporating sprays with coalescence [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF].

On the one hand, [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] had demonstrated the capability of an Eulerian multi-fluid model to capture the physics of polydisperse dense evaporating sprays. This approach relies on the derivation of a semi-kinetic model from the Williams equation using a moment method for velocity conditioned by size, but keeping the size as a variable [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. This function is then discretized using a "finite-volume" approach that yields conservation equations for mass, momentum (and eventually other properties such as temperature) of droplets in fixed size intervals called "sections" extending the original work of Tambour, Greenberg and collaborators [START_REF] Greenberg | An opposed jet quasi-monodisperse spray diffusion flame[END_REF][START_REF] Greenberg | On the origin of spray sectional conservation equations[END_REF] for evaporating sprays. Even though this approach has recently been extended to second order by Laurent [START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] and Dufour [START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF][START_REF] Dufour | A second-order multi-fluid model for evaporating sprays[END_REF] for evaporating sprays, the necessity to discretize the size phase space can be a stumbling block for applications implying a drift velocity in the size phase space such as evaporation since it introduces intrinsic numerical diffusion and then requires fine grids in size phase space. Moment methods such as DQMOM introduced by [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF], on the other hand, do not encounter this limitation and have proved in [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF] to perform quite well on a number of configurations.

Nevertheless, there exists one configuration for which the moment methods experience a difficulty as pointed out in [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF]. In some situations, a drift velocity, that is, the rate of change due to continuous phenomena of the internal coordinate such as the size of the particles, has sometimes to be taken into account; it can be either positive as for molecular growth, or negative such as for evaporation of droplets or oxidation of soots. When negative, it leads to the disappearance of droplets/particles thus creating a negative flux at zero size. Its closure requires an evaluation of the reconstructed NDF at zero size from the knowledge of some moments. Except in the particular case where the NDF is a sum of Dirac delta functions and where the flux at zero size is very singular, the nature of this information, i.e. pointwise in internal coordinate, and its influence on moment dynamics results in a difficulty from both a modeling and a numerical point of view. Such a problem is not encountered with multi-fluid modeling [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF][START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF] because of the natural averaging procedure in one section. Some preliminary solution was proposed in [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF], which offered a stable numerical scheme but not entirely satisfactory since it was based on a number of constraints and did not provide a good enough level of accuracy for continuous NDF in the size phase space. In the present contribution we obtain a comprehensive solution to this important issue.

We make use of a diversity of ingredients. First, the evaluation of pointwise information at a given size from the knowledge of a set of moments is conducted through Maximum Entropy reconstructions [START_REF] Mead | Maximum entropy in the problem of moments[END_REF] and associated algorithms. The algorithm relies on convex optimization under constraints and is very efficient when the NDF is far enough from the frontier of the moment space. In fact, we characterize this notion of distance from the frontier in order to make a clear separation between cases where the modeling of the problem relies on a sum of Dirac delta functions and is discrete in the size phase space, from the case where a continuous description is made possible even if the distribution reaches regions in the moment space very close to its frontier. In most cases, the discrete case can be seen as an idealization of a continuous distribution with little variance around one or several peaks. Second, once the flux at zero size can be evaluated, we obtain a closed system of ordinary differential equations (ODEs) and we have to design an algorithm which guarantees the preservation of the moment space, the geometry of which is far from obvious and based on algebra properties; it can not be achieved by using standard ODEs integrators. We couple the resolution along the characteristics (already used in [START_REF] Ramkrishna | Population Balances: Theory and Applications to Particulate Systems in Engineering[END_REF]) to the ideas of the kinetic schemes for hyperbolic equations [START_REF] Perthame | Kinetic formulation of conservation laws[END_REF] and make the link with QMOM [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF] and DQMOM [START_REF] Marchisio | Quadrature method of moments for aggregationbreakage processes[END_REF] ap-proaches. The obtained stability so far strongly relies on a invariance property based on a special constant surface regression rate. The proposed method is then shown to be well-suited for general evaporation laws/negative drift velocities and its efficency is then exhibited through several numerical results for various size distributions.

Since we introduce some new tools in order to resolve the flux evaluation, we are also able to introduce a new Eulerian type of description which will combine both the flexibility of Eulerian multi-fluid models and the efficiency of DQMOM. Consequently, we still divide the size phase space into "sections", i.e. fixed size intervals as in the multi-fluid approach [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], but we need a very precise description of the NDF inside this section through a set of moments as in the QMOM or DQMOM. Such a perspective is very attractive, however, it also yields the difficult task of evaluation fluxes at the section boundaries from the knowledge of internal moments; such a task is equivalent as the evaluation of the flux of droplets at zero size and the same tools can be used. It yields a very precise description of evaporation with a minimal number of variables, since only one, two or three sections are needed, as well as a low computational cost. It offers a very stable numerical scheme as compared to other solutions such as for example presumed NDF [START_REF] Mossa | Extension polydisperse pour la description euler-euler des écoulements diphasiques réactifs[END_REF]. To our knowledge, the results presented in this paper allow a real breakthrough in the field by originally combining a variety of mathematical tools and set the basis for more complicated problems in spatially extended configurations and multi-dimensional problems. A special attention has to be devoted to transport in physical space; a first study is provided in [START_REF]Modeling of polydisperse sprays using a high order size moment method for the numerical simulation of advection and evaporation[END_REF][START_REF] Massot | Eulerian multi-fluid models : modeling and numerical methods[END_REF] and it will be the subject of a companion paper [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF].

The paper is organized as follows. In the remaining of this first Section, we will consider the PBE of particles or the Williams equation for sprays and obtain a model equation which is characteristic of the difficulties we will have to face with a negative drift velocity such as evaporation of droplets or oxidation of soot particles. Considering the dynamics of moments for such an equation, we identify the critical issue of obtaining the flux of disappearing particles or droplets from the moments and show that this flux, if not zero, will define the moment dynamics. The question of the flux at zero size is similar to a more general problem which is presented in § 1.3 and relates to a Eulerian multi-fluid method. Section 2 tackles the issue of the Hausdorff finite moment problem as well as the geometry of the finite moment space, lower principal representation and NDF reconstruction through the Maximum Entropy formalism; the limitation of such an approach near the frontier of the moment space is precisely described. With these ingredients, Section 3 is devoted to the detailed presentation of the new method allowing an accurate and stable treatment of the flux problem in the case of a constant drift velocity based on a special invariant property of the moment space during transport in size phase space at contant velocity. The link with a DQMOM type of approach is provided in § 4 and allows to highlight the novelty of the proposed approach. Two extensions of the proposed approach are provided in § 5 for piecewise constant drift velocities based on pressureless gas dynamics, and 6 for treating arbitrary drift velocities. The results Section 7 is devoted to the numerical illustration of the proposed algorithm for a series of representative test-cases. Two Appendices allow to explain why classical ODEs integrators are helpless for the proposed problem and provide illustrations of the ME reconstruction behavior close to the frontier of the moment space.

1.1. Droplets, aerosols or particles : the same model problem. Let us start with the Williams transport equation [START_REF] Williams | Spray combustion and atomization[END_REF] for the joint size, velocity number density function f (t, x; S, u) is

∂ t f + u • ∂ x f -∂ S (Rf ) + ∂ u • (F f ) = Γ, (1.1)
where R is the evaporation rate, F is the drag force acting on the droplet, and Γ is the coalescence term such as in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF][START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF].

Since we focus in this paper on the difficulties associated with evaporation process, we will reduce this equation to the "model" partial differential equation (PDE) on a homogeneous number distribution function n(t; S) which satisfies :

∂ t f -∂ S (R f ) = 0, (1.2)
where the evaporation coefficient R(S) is a positive function of S. Such a PDE can be recovered from equation (1.1) by assuming f (t, x; u, S) = n(t; S)δ(u), F = 0 and Γ = 0. Various levels of modeling can be considered; the first one is the classical

d 2 evaporation law R(S) = R 0 1 I [0,+∞[ (S)
, where R 0 is a constant, which does not depend on droplet size; it is a simplified law for droplet evaporation. This model will be investigated in a first part of the paper. Then, in the context of multi-fluid models, a size discretization is introduced and the evaporation law thus become constant in each size interval [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] in such a way that, even considering the d 2 law in each size section results in a discontinuous evaporation law at a global level. Finally, the third one is an arbitrary evaporation law.

In the case of aerosols or soots particles, following [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF], the PBE on the NDF f (t, x; ξ), where ξ denote the internal variables, reads :

∂ t f + ∂ x • [< u|ξ > f ] -∂ x • (D x ∂ x f ) = S ξ , (1.3)
where < u|ξ > denotes the mean velocity conditioned on the internal coordinate ξ. The diffusion coefficient reads D x and the source term S ξ (t, x; ξ) is due to aggregation, breakage and drift. The term we are interested in is the drift velocity in the source term S ξ . Assuming the system homogeneous, < u|ξ >= 0 and assuming no aggregation-breakage of the particles, the variable describing the internal coordinate can be considered to be merely the size of the particles. We consider a negative drift velocity, without aggregation nor breakage, S ξ (t; ξ) = ∂ ξ (R(ξ) f (t; ξ)). If the size is described by the particle surface, for soot particles for example, we end-up with the same model as the d 2 law since oxidation is proportional to the particle surface.

Finally, we could also consider the general dynamic equations for the continuous distribution function proposed in Friedlander [START_REF] Friedlander | Smoke, Dust, and Haze, Fundamental of Aerosol Dynamics[END_REF] for aerosols where the drift velocity is denoted "the migration velocity through the size space" (Chapter 11) and neglect some phenomena in order to end-up with the proposed the identified PDE model, thus showing the common background of the proposed study.

Note: Let us notice that it is reasonable to start working with such a simplified configuration since we aim at isolating the difficulty of evaluating the flux associated to the drift velocity R. This model will be relevant for droplets when size-conditioned dynamics and coalescence/break-up are added such as in [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF] or for soots as in [START_REF] Zucca | Implementation of the population balance equation in cfd codes for modelling soot formation in turbulent flames[END_REF]. However, using either multi-fluid models or DQMOM approaches, we know how to treat transport and drag as well as droplet/droplet interactions [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF][START_REF] De Chaisemartin | Evaluation of Eulerian Multi-Fluid versus Lagrangian methods for the ejection of polydisperse evaporating sprays by vortices[END_REF][START_REF] Massot | Eulerian Multi-fluid method for the numerical simulation of evaporating polydisperse sprays : modelling and numerical issues in multi-dimensional configurations[END_REF][START_REF] Fréret | Pulsated free jets with polydisperse spray injection : Experiments and numerical simulations[END_REF] (or aggregation, breakage and diffusion) so that solving the flux problem is a necessary starting point. We consider a monovariate NDF and assume that the internal coordinate can only be positive (the moments will not be taken over the whole real line such as in [START_REF] Zucca | Implementation of the population balance equation in cfd codes for modelling soot formation in turbulent flames[END_REF]); the extension to bi-variate distributions, such as in [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF] for soot geometry modeling, as well as the coupling with the other physical processes, will be discussed in the conclusion.

1.2. Flux at zero size for moment dynamics : a model equation. From the previous subsection, we will focus on the simple model PDE on the number distribution function f (t, S) which satisfies :

∂ t f -∂ S (R f ) = 0, (1.4)
where the drift velocity R reads, in the first part of this paper, R = R 0 1 I [0,+∞[ (S), with R 0 = 1, that is the drift velocity is negative and has been non-dimensionalized to -1. The method will be first derived for this drift velocity before being extended to more general one.

It can be noted that this is a transport equation with discontinuous coefficients and that it has well-defined measure solutions [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF]. Roughly speaking, the jump in velocity from zero for negative sizes to one for positive sizes creates the crossing of characteristics and results in a measure concentration at zero size. Consequently, it can be shown that the exact solution of equation (1.4) with a positive initial data f 0 (S) on [0, +∞) (eventually a positive measure) reads :

f (t, S) = n 0 (t) δ(S) + f 0 (S + t), S ≥ 0, (1.5)
where n 0 (t) = t 0 f 0 (S) dS is exactly the number of droplets/particles which have disappeared at time t. Thus this number density concentration is well defined mathematically in the sense of duality solution [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF].

The principle of moment methods is to be able to write, eventually using quadrature formulas, a set of dynamical equations on a set of moments of the NDF. We illustrate this point of view in the context of the model PDE. For N a given integer, let M denote the vector of surface moments :

M =    m 0 . . . m N    , m k (t) = Smax 0 S k f (t, S) dS. (1.6)
where S max is a given maximal size. Using a natural integration by parts, we get the following set of ordinary differential equations for the evolution of vector of moments M :

d t M = -A M -φ -, (1.7)
where besides the flux at zero size φ -, we introduce the translation nilpotent matrix A :

φ -= f (t, 0)      1 0 . . . 0      , A =         0 0 1 0 2 . . . . . . . . . 0 N 0         .
It is important to note that this form of the system is not closed since the value of the flux can not be found directly from the values of M. The knowledge of the value of f (t, S) at the zero boundary for all t is equivalent to fully solving the original "kinetic" problem. These values then have to be provided as functions of M(t), which constitutes the challenge. The problem we just raised is similar if one wants to use a DQMOM method, since the source terms for the evolution of abscissas and weights are found from the equations of the moments and the DQMOM dynamics is then completely related to a choice for the form of the flux f (t, 0) in terms of moments or abscissas and weights [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF].

The measure concentration at zero size has no impact on higher moments than the total number density and does not result in any complex behavior, except that the flux of disappearing droplet which concentrate in a Dirac delta function at zero size results in a modification of the whole set of moments through the coupling matrix A.

1.3. Similar problematic : combined multi-fluid and moment dynamics. In fact, the question of evaluating the flux at a given size can be considered in a more general framework. We consider, as in the multi-fluid model, "sections", a section being defined by a fixed size interval [S min , S max ]. For N a given integer, let M denote the vector of droplet surface moments which is the same as in (1.6) but with m k (t) = Smax Smin S k f (t, S) dS. The size phase space can be eventually discretized into several sections [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF][START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF]. The moment dynamics then read :

d t M = -A M -φ -+ φ + , (1.8) φ -= f (t, S min )      1 S min . . . S N min      , φ + = f (t, S max )      1 S max . . . S N max      ,
Solutions have already been proposed in order to define reasonable fluxes up to order two by taking moment vectors of dimension one [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF][START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] or two [START_REF] Dufour | A second-order multi-fluid model for evaporating sprays[END_REF][START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF]. However, in the present paper we are not concerned with the perspective of numerical analysis which characterizes the behavior of a numerical scheme as the section width is approaching zero, but with the design of a new scheme for which the use of only up to 3 sections is enough for a very good accuracy and for which we can use some ingredients from QMOM or DQMOM. Before presenting our scheme, let us recall briefly the framework of finite moment spaces and NDF reconstruction from its moments.

2. Finite moment spaces and Hausdorff finite moment problem. In order to solve the problem (1.7) or a similar problem, a reconstruction of the number density function f can be used to have an expression of the flux f (0, t). Such a reconstruction is possible if and only if the numerical method allow to stay in the moment space. We then first give a brief overview of the definition and properties of this space, from [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]. Then, we will show how a positive function can be reconstructed from a finite number of its moments, which is the Hausdorff finite moment problem.

2.1. Moment space. For the sake of simplicity, we consider the non-dimensional size interval [0, 1]. If P denotes the set of all probability measures on the Borel sets of the interval [0, 1], then the N th-moment space M N on the interval [0, 1] denotes the set of moment vector of dimension N , M N ⊂ [0, 1] N :

M N = {c N (µ)|µ ∈ P}, c N (µ) = (c 1 (µ), . . . , c N (µ)) t , c k (µ) = 1 0 x k dµ(x).
Let us recall that, since we consider probability measures, we always have c 0 = 1. This set is convex but has a complex geometry.

For c N ∈ M N , let P(c N ) denote the set of all measures µ ∈ P with moments up to the order N equal to c N = (c 1 , . . . , c N ) t and

c + N +1 (c N ) = max µ∈P(cN ) c N +1 (µ), c - N +1 (c N ) = min µ∈P(cN ) c N +1 (µ),
If c N is in the interior of M N then c - N +1 < c + N +1 and P(c N ) is infinite [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]. A canonical moment sequence (p k ) k≤N +1 can then be defined by :

p k = c k -c - k (c k-1 ) c + k (c k-1 ) -c - k (c k-1 ) . (2.1)
Note that the canonical moments vary in the interval [0, 1] and they remain invariant under a linear transformation of the measure. Moreover, p N +1 is 0 or 1 if and only if c N +1 is at the boundary of M N +1 and P(c N +1 ) is then a singleton. An algorithm is given in [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] in order to go from moments to canonical moments and conversely.

In our case, since we are dealing with number density functions over an interval (S min , S max ). Let us first denote M N (S min , S max ) the N th-moment space corresponding to this interval. For any non zero moment vector M = (m 0 , m 1 , . . . , m N ) t ∈ R N +1 , we define the vector of normalized moments M ∈ M N (S min , S max ) ⊂ R N by M = (m 1 /m 0 , . . . , m N /m 0 ) t . Let us then denotes M N (S min , S max ) the N th-moment space of such moment vectors M. Moreover, the first canonical moments can be written :

p 1 = m 1 -S min m 0 (S max -S min )m 0 , p 2 = m 0 m 2 -m 2 1 (m 1 -S min m 0 )(S max m 0 -m 1 ) , p 3 = (S max m 0 -m 1 )[(m 1 -S min m 0 )m 3 -(m 2 -S min m 1 -S 2 min m 0 )m 2 -S 2 min m 2 1 ] (S max -S min )(m 0 m 2 -m 2 1 )(-m 2 + (S max + S min )m 1 -S max S min m 0 )
.

The geometry of the space of the canonical moments vectors is much more simpler than the one of the moments since it is then the cube [0, 1] N . The canonical moments are then very useful in order to check the belonging to the moment space.

2.2. Some solutions to the Hausdorff finite moment problem. The Hausdorff finite moment problem for the moments M is : finding a positive real valued function f M defined on [S min , S max ] such that

M = Smax Smin f M (x)      1 x . . . x N      dx (2.2)
As seen in the previous Section, if M belongs to the interior of M N (S min , S max ), there is an infinity of solutions.

Similar problems were tackled in different ways in the context of solving the problem (1.4). With the use of only one moment, a fixed profile can be used [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], like a constant, leading to an at most first order method in term of the size discretization [START_REF] Laurent | Analyse numérique d'une méthode multi-fluide Eulérienne pour la description de sprays qui s'évaporent[END_REF][START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF]. Note that the moment of interest is then the moment of order 3/2 and the positivity of such moment is sufficient in order to stay in the moment space. Two kinds of extension of such methods were developed. First, in [START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF], N moments were used, with some polynomial functions for f M , in order to obtain a N th order method. But the moment space of positive polynomial functions is only a sub-space of the moment space and a modified method is introduced in this article in the case of two moments, using a reconstruction with a bi-affine function instead of an affine one. An other extension were developed in [START_REF] Dufour | A second-order multi-fluid model for evaporating sprays[END_REF] with the use of moments of order 0 and 3/2. It is then show that a function of the type exp(ax + b) can be used and a preservation of the moment space was also taken for.

Here, some reconstructions are introduced for an arbitrary number of moments in link with the structure of the moment space.

Lower principal representation.

For any point in the interior of the N th-moment space M N (S min , S max ), with N = 2n -1, it can be shown that there exists a unique lower principal representation (i.e. for which the moment m N +1 is minimal; (m 0 , . . . , m N +1 ) then belongs to the boundary of M N +1 (S min , S max )) under the form of a sum of n weighted Dirac delta function, the support or roots of which are in (S min , S max ). This support can be shown to be the root of some orthogonal polynomials associated to Hankel determinants [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]. It is this lower principal representation that is used in quadrature method of moments (QMOM) introduced in [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF] and further used in the DQMOM [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF].

Then, for any M ∈ Int(M N (S min , S max )), there exists one unique representation of the vector of moments by using weights (ω i ) i∈ [1,n] and abscissas (S i ) i∈ [1,n] such that :

M =      m 0 m 1 . . . m N      =      n i=1 ω i n i=1 ω i S i . . . n i=1 ω i S N i      , (2.3)
where we recall that N + 1 = 2n. The corresponding number density function is then :

f - M (S) = n i=1 ω i δ(S -S i ). (2.4)
Abscissas S i are the roots of the polynomial

P = X n + n-1 k=0 σ k X k with the σ k given by :       m 0 m 1 . . . m n-1 m 1 . . . . . . . . . m n-1 . . . . . . m 2n-1            σ 0 σ 1 . . . σ n-1      = -      m n m n+1 . . . m 2n      .
The weights are then solutions of a linear system. A QD algorithm is given in [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] to compute these weights and abscissas. A such reconstruction is used for quadrature formulas [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF] but cannot be used in order to find the value of the flux, which is a punctual value of the NDF.

NDF reconstruction through the Maximum Entropy formalism.

A smooth approximated NDF can be obtained by the Maximum Entropy (ME) reconstruction [START_REF] Mead | Maximum entropy in the problem of moments[END_REF]. The ME method yields a smooth distribution which maximizes the following Shannon entropy from information theory :

H[f ] = - Smax Smin f (x) ln f (x)dx.
The existence of such ME distribution is shown in [START_REF] Tagliani | Hausdorff moment problem and maximum entropy: a unified approach[END_REF] as soon as the vector of moments M belongs to the interior of the moment space M N (S min , S max ). This is a standard constrained optimization problem, leading to the following explicit representation of the ME approximate :

f ME M (x) = exp   - N j=0 ξ j x j   . (2.5)
to be supplemented by the condition (2.2). It can be seen as a generalization of the reconstruction used in [START_REF] Dufour | A second-order multi-fluid model for evaporating sprays[END_REF] except that the considered moments were moments of order 0 and 3/2 instead of integer moments here. The coefficients ξ 0 , . . . , ξ N are the Lagrange's multipliers. We then just need to minimize the following convex potential :

∆ = Smax Smin   exp   - N j=0 ξ j x j   -1   dx + N j=0 ξ j m j .
Indeed, its stationary points are given by

∂∆ ∂ξ i = 0 ⇒ Smax Smin x i exp   - N j=0 ξ j x j   dx = m i .
Numerically, a Newton method is used, as in [START_REF] Mead | Maximum entropy in the problem of moments[END_REF] : starting from initial choices ξ = (ξ 0 , . . . , ξ N ) ′ , updated ξ's are defined from

ξ + = ξ -H -1 (M 0 -< X > ξ ),
where < X > ξ = (< x 0 > ξ , . . . , < x N > ξ ) ′ is the vector of approximated moments, with

< x k > ξ = Smax Smin x i exp   - N j=0 ξ j x j   dx (2.6)
and H is the Hessian matrix defined by H i,j = ∂∆ ∂ξi∂ξj =< x i+j > ξ for i, j = 0, . . . , N . The used numerical procedure is the same as in [START_REF] Mead | Maximum entropy in the problem of moments[END_REF]. In this paper, they found that a double-precision 24-point Gaussian quadrature method very efficiently produces the demanded accuracy for < x k > ξ . The quadrature points are computed once for the interval [0, 1] and a change of variable is used for all the integrals in order to come back to this interval. Because H is symmetric positive definite, a classical Choleski decomposition is used in order solve the linear system [START_REF] Press | Numerical Recipes in Fortran 77[END_REF]. A number of iterations from 4 to 15 is needed for the Newton method to converge, with an accuracy of 10 -6 on the normalized moments, from the initial value ξ = (-ln(m 0 )/(S max -S min ), 0, . . . , 0) ′ corresponding to a constant reconstruction preserving the moment of order 0. This number of iterations depends on the distance to the boundary of the moment space : the closest the moment vector is to this boundary, the largest this number of iterations is. Since this issue is not treated in [START_REF] Mead | Maximum entropy in the problem of moments[END_REF], we propose to treat it in the next paragraph in details. Finally such a reconstruction provides a value for the flux at zero size and thus a closure for problem (1.7) in terms of the moment vector.

An example is given in Fig. 2.1-right in case N = 3 and we propose a comparison in Fig. 2.1-left with a one moment (of order 3/2) constant recontruction in 12 sections, like in [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] and a reconstruction in Fig. 2.1-top-right using two moments and a biaffine shape on each 4 sections, like in [START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF]. Let us note that the numerical method using the first reconstruction needs 12 moments, while with the second one we need 8 moments and with the ME reconstruction we need 4 moments. The accuracy of the resulting approaches will be analysed in § 7 where all the results and discussion are gathered.

2.2.3.

Behavior at the frontier of the moment space. Let us note that such an algorithm can lead to numerical difficulties due to ill-conditioned H matrices appearing in two cases. Such a bad condition number is encountered for moment vectors too close to the boundary of the moment space and results in highly oscillating ξ coefficients in the ME approximation for a reasonable number of moments. It can also appear in the case of a large vector of moments as a consequence of the peculiar high-dimensional geometry of the moment space. However, in the present contribution, we will remain low-dimensional, typically considering up to 6 moments, for which the algorithm will prove to be very efficient as long as we stay away from the frontier of the moment space. Since distributions at the frontier of the moment space consist in a sum of Dirac delta functions for which the ME reconstruction is not possible, we have to clarify the distance between such measure distributions and the closed subset of the interior of the moment space for which our formalism is operational, and eventually to make the link between the frontier and such a subset.

A first piece of answer is provided with Figure 2.2. Indeed, for a reconstruction with 4 moments, the ME reconstruction can be done, even for moment vector quite close to the boundary of the moment space. For three values of the first normalized moment, 0.1, 0.5 and 0.9, we compare a 2D-cut of the moment space in the (m 2 , m 3 ) plane (the three traces of these cuts in the (m 1 , m 2 ) plane are presented in Figure 2.2 bottom-right) with the set of moments in this place which can be attained with the ME formalism within an error for the moment vector evaluation below 1% with the presented Gaussian quadrature approach. In fact the quasi-totality of the the moment space can be covered using such an approach.

However, there is a second numerical difficulty when approaching the frontier of the moment space and it comes from the lack of accuracy of the Gaussian quadrature for the computation of the integrals used in the algorithm through equation (2.6). It is quite obvious when seeing the type of reconstruction obtained in such cases (see subsection 7.1 and Appendix B). Moreover, using a more precise integration method for the evaluation of the Hessian matrix allows ME reconstructions for moment vectors close at least till 1% of the boundary (that is to say, at least when all canonical moment are in the interval [0.01, 0.99]). Thus, it becomes clear that we are able to reach almost the entire moment space with the proposed formalism and we even can quantify the small subset for which a ME formalism is out of reach. Within this context, the last question to be answered is how to use the algorithm proposed in this paper for this small subset. We will provide a final answer to this question in Appendix B, where such distributions are also illustrated for the sake of completeness.

Eventually, let us remind the reader that the whole purpose of the paper is to provide a comprehensive answer for the flux evaluation and associated numerical method; when we are at the boundary of the moment space, the flux evaluation is not a problem any more; consequently, in the following we will stick with smooth enough size distribution functions with moment vector within the limits previously indicated and this will not be restrictive.

3.

A new numerical scheme with flux evaluation for stable and accurate moment dynamics in the case of constant evaporation law. The case of a constant evaporation law is first studied in order to present the principle of the new numerical schemes. In the following, we will see how this schemes can be extended to more general evaporation laws.

As soon as a reconstruction is done for the NDF from its moments, problem (1.8) could be solved directly with the pointwise value of this reconstruction at S = S min , using a dedicated stiff ODE-solver. It is done in Appendix A and it is then shown that the stability is not warranted.

Before introducing the first key ingredient of the numerical scheme which guarantee this stability, let us underline an invariance property which is strongly related to the linear transport equation (1.4).

3.1. Invariance property of moment space through evaporation with zero fluxes. Lemma 3.1. As long as the fluxes φ -(t) and φ + (t) remain zero, the evolution of the NDF at the kinetic level is a pure translation in the size phase space and the dynamical evolution of the moments is in closed form. Moreover, we have :

M(t) = exp(-t A) M(0), exp(t A) =         1 0 t 1 t 2 2 t . . . . . . C j-1 i-1 t i-j . . . . . . t N . . . . . . N t 1         .
Let us note that the assumption φ -(t) remain zero implies that f (0, η) = 0 for η ∈ [S min , S min + t]. This lemma has a strong implication : when the dynamics of evaporation is described by a translation in the size phase space (i.e. no fluxes of droplets appearing or disappearing) and when we choose a set of incremental moments in this "natural" evaporation variable, the moment space is invariant by this dynamics. Thus, it is not necessary to invoque quadrature methods since the dynamics is closed. It is then interesting to characterize the related evolution of the abcissas and weights since we are sure to remain in the interior of the moment space of order N .

Lemma 3.2. The evolution equation on the moments can be written in terms of the abcissas and weights in the following form :

     n i=1 ω i (t) n i=1 ω i (t) S * i (t) . . . n i=1 ω i (t) S * N i (t)      =      n i=1 ω i (0) n i=1 ω i (0) S i (0) . . . n i=1 ω i (0) S N i (0)      (3.1)
where the transport term associated to the exponential has been directly included in the new abcissas S * i (t) = S i (t) + t. Proof. The proof is straightforward once it has been noticed that we have the following equality :

exp(t A)      n i=1 ω i (t) n i=1 ω i (t) S i (t) . . . n i=1 ω i (t) S N i (t)      =      n i=1 ω i (t) n i=1 ω i (t)(S i (t) + t) . . . n i=1 ω i (t)(S i (t) + t) N      . (3.2)
3.2. Integrated version of the dynamical system. Since the pointwise definition of the boundary fluxes is difficult to evaluate, we change the point of view and try to construct, like in the kinetic schemes for hyperbolic equations, an integral formulation making use of the underlying kinetic equation.

The integral formulation of the preceding system (1.7) of ODE's is the following one :

exp(t A) M(t) = M(0) + t 0 exp(η A) (-φ -(η) + φ + (η)) dη. (3.3)
We then want to use the peculiar structure of the integral form of the solution in order to obtain a more convenient form of the fluxes through two arguments : the characteristic solution of the original conservation equation (such as in the kinetic schemes [START_REF] Perthame | Kinetic formulation of conservation laws[END_REF]) as well as a moment formulation of the fluxes. This is the purpose of the following lemma.

Lemma 3.3. The integral form of the solution can be rewritten :

exp(t A) M(t) = M(0) -Ψ -(t) + Ψ + (t) (3.4)
where the fluxes read :

Ψ -(t) = Smin+t Smin f (0, β)      1 β . . . β N      dβ, Ψ + (t) = Smax+t Smax f (0, β)      1 β . . . β N      dβ.
Proof. Let us first focus on the flux of disappearing droplets :

t 0 exp(η A) φ -(η) dη = t 0 f (η, S min )         1 η 1 0 η 2 2 η . . . . . . C j-1 i-1 η i-j . . . . . . η N . . . . . . N η 1              1 S min . . . S N min      dη
Using on the one side :

        1 η 1 0 η 2 2 η . . . . . . C j-1 i-1 η i-j . . . . . . η N . . . . . . N η 1              1 S min . . . S N min      =      1 η + S min . . . (η + S min ) N      ,
and the fact that along the characteristics curves, the solution is constant so that we have f (η, S min ) = f (0, η + S min ); we then get :

t 0 exp(η A) φ -(η) dη = t 0 f (0, η + S min )      1 η + S min . . . (η + S min ) N      dη.
A simple change of variable β = η + S min allows to conclude on the first equality. The treatment of the gain term is readily the same. There are several important points that have to be underlined. First, it can be easily shown that in the configuration where the fluxes are both set to zero, the right variable is then exp(t A) M(t), which is invariant, the evolution of the moments being then only modified by transport through the multiplication by exp(-t A). The integral form of the ODE system (3.4) then decouples the transport per se from the evolution of the moments through the fluxes.

The second point relates to the numerical scheme to be used in order to resolve the system (1.7). It should be noted that in the situation where f (t, S max ) = 0, i.e. when there is no flux coming from the upper sections and when S min = 0, a flux based on pointwise value of f only participates in the evolution of the zeroth order moment, since S N min = 0 for N greater than zero. With this new integral form, we then get another way of integrating the system of ODE's without the need for a pointwise reconstruction for which the stability is not warranted as shown in Appendix A : it is seen that the resolution of the ODE system with classical methods like Runge Kutta or Euler leads to an exit of the moment space and a divergence of the method. The integral form allows to remain in the moment space M N (S min , S max ) of order N as presented in the following proposition.

Proposition 3.4. For a given set of moments M(0) = M 0 at time t = 0 such that m 0 > 0, we assume that the vector of normalized moments M 0 belongs to the interior of the moment space of order N associated to the interval [S min , S max ], that is M 0 ∈ Int(M N (S min , S max )). Then, for any measurable positive function f M 0 defined on the interval [S min , S max + T ]and solution of the Hausdorff finite moment problem (2.2) in the interval [S min , S max ], we have

M(t) ∈ Int(M N (S min , S max )), (3.5) as long as t ∈ [0, T ].
Proof. We use the integral form (3.4) and notice that

M 0 = Smax Smin f M 0 (η)      1 η . . . η N      dη.
Consequently, we obtain that :

exp(t A) M(t) = Smax+t Smin+t f M 0 (η)      1 η . . . η N      dη.
This last equality implies that exp(t A) M(t) belongs to the interior of the moment space M N (S min + t, S max + t) so that naturally, M(t) belongs the moment space M N (S min , S max ) and the proof is complete.

The choice of such function f M 0 was discussed in previous Section and we use then the ME reconstruction in such a way that the flux can be given by :

Ψ -(t) = Smin+t Smin f ME M 0 (β)      1 β . . . β N      dβ. (3.6)
Let us remark that, as explained in § 2.2.2, the ME reconstruction is not easy to obtained for moment vectors closed to the boundary of the moment space. But, as presented in Appendix B, when the moment vector become closer to a point of this boundary, the ME reconstruction becomes closer and closer of the lower principal representation and of the sole solution of the Hausdroff moment problem corresponding this point, that is to say a sum of Dirac delta functions. Let us denote i ω i δ Si one of these functions. Then, the flux can be compute from this reconstruction :

Ψ -(t) = ω i 1, S i , . . . , S N i t , (3.7)
where the sum is taken for all i such that S i ∈ [S min , S min + t]. It can be seen that this type of flux will allow the preservation of the moment space as it will be shown in the following for the flux (3.6).

3.3. Link with abcissas and weights of the lower principal representation. We know that for any M 0 ∈ Int(M N ), there exists one unique lower principal representation of the vector of moments by using weights (ω i ) i∈ [1,n] and abcissas (S i ) i∈ [1,n] . Since we know that the evolution of the vector of moments takes place in the interior of the moment space Int( M N ), we can write :

exp(t A)      n i=1 ω i (t) n i=1 ω i (t) S i (t) . . . n i=1 ω i (t) S N i (t)      =      n i=1 ω i (0) n i=1 ω i (0) S i (0) . . . n i=1 ω i (0) S N i (0)      -Ψ -(t) + Ψ + (t).
This formulation can be considered as very interesting since it can be readily seen from it that in the presence of zero fluxes, the weights remain unchanged and the abcissas are only translated at the right velocity. The key variable is not S i (t) but in fact S i (t) + t. This will be confirmed in the next Section where we make the link between the present approach and the DQMOM approach.

Moreover, this formulation will allow us to emphasize an algorithm for the scheme. It needs the calculation of abscissas and weights from moments. It can be noted that the computation of M(t) from M(0) does not need the separation between flux and transport and could be done, like in [START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] or [START_REF] Dufour | A second-order multi-fluid model for evaporating sprays[END_REF], using directly the following formula :

M(t) = exp(-t A) M(0) -exp(-t A) Ψ -(t) + exp(-t A) Ψ + (t).
However, the calculation of abscissas and weights is not an expensive operation if we only have n = 2 or n = 3 like for the results presented in this paper. Moreover, they can be very useful for more general problems where phenomena like collision have to be taken into account, which can be done through a quadrature method using these abscissas and weights.

3.4. Algorithm of the new scheme. The new scheme introduces a discretization 0 = S 0 < S 1 < . . . < S K = S max of the support [0, S max ] of the distribution function. The initial moments M k (0) on each interval [S k-1 , S k ], called section as for the classical multi-fluid method, are known. The resolution algorithm corresponding to the integral formula (3.4) is used in order to find the moments in each section after a time t. This algorithm is :

1. In each section [S k-1 , S k ], a reconstruction f k is done of the distribution from its moments M k (0) by the ME method and the flux Ψ k -(t) with the section k -1 is computed by (3.6)

with f ME M k (0) = f k . 2.
The weights ω k i and the abscissas S k i corresponding to the moments M k (0) -Ψ k -(t) + Ψ k+1 -(t) are computed using the QD algorithm. Since these moments are in fact the moments of

i f i (S)1 [Si-1,Si[ (S) on [S k-1 + t, S k + t], the abscissas are in ]S k-1 + t, S k + t[. 3.
The moments M k (t) corresponding to the weights ω k i and the abscissas S k it are computed. It gives an approximation of the moments of f (t, .) on ]S k-1 , S k [. In practice, a time discretization associated to the CFL like condition ∆t ≤ ∆S is introduced and the previous algorithm governs the time evolution of both the moments in the various sections as well as the associated abcissas and weights.

The cost of this algorithm is essentially the cost of the ME reconstruction in each section by the Newton method and the cost of the QD algorithm. Since the number of moment we will use is small, and since the Newton method generally only requires around five iterations in order for the moments of the approximated number density function to be at a maximum relative error up to 10 -6 from the original moments, the cost of the method is very moderate.

Let us underline once more that such a method has the major advantage, with a very moderate cost, to preserve the moment vector in the singular transport equation with a very stable model for the flux of disappearing droplets/particles. Moreover, after making the link between such a scheme and the DQMOM approach for a better in understanding the novelty of the present contribution, we will see how this algorithm can be easily extended to more general evaporation laws, with the same type of steps.

4. Link with a DQMOM approach. For the sake of completeness, let us first recall the DQMOM approach in order to resolve the evaporation dynamics in the size phase space and then adapt it to the particular structure of our problem.

4.1. DQMOM and new DQMOM formalism of the evaporation process. In order to solve the transport equation (1.4) using DQMOM, we consider the following approximation of f :

f (t, S) = n i=1 ω i (t) δ(S -S i (t)), (4.1)
from which we have :

∂ t f = n i=1 a i (t) δ(S -S i (t)) - n i=1 (b i (t) -S i (t)a i (t)) δ ′ (S -S i (t)), (4.2)
where we have defined for commodity

a i (t) = d t ω i (t) b i (t) = d t (ω i (t) S(t)) (4.3)
The evolution of the approximated moments mk depends linearly on the a i and b i :

d t mk = n i=1 (1 -k) a i (t) S k i (t) + n i=1 k b i (t) S k-1 i (t) (4.4)
However, the variables for which we have the simplest evolution equation are exp(t A)M(t), so that it is interesting to look at the system of ODE's satisfied by exp(t A) M(t).

Lemma 4.1. The evolution equation of the translated approximated moments reads :

d t exp(t A) M = n i=1 (1 -k) a i (S i + t) k + n i=1 k b i + ω i + t a i (S i + t) k-1 (4.5)
Considering the new variables :

b * i = b i + ω i + t a i = d t [ω i (S i + t)], S * i = S i + t, (4.6)
the new constraints of DQMOM defining the source terms are given by :

n i=1 (1 -k) a i (S * i ) k + n i=1 k b * i (S * i ) k-1 = -f (t, S min ) S k min + f (t, S max ) S k max . (4.7)
Proof. Let us first notice that :

d t [exp(t A) M] = exp(t A)[A M + d t M],
and

(A M + d t M) k = n i=1 (1 -k) a i S k i + k b i S k-1 i + k ω i S k-1 i ,
Consequently, we obtain :

(exp(t A)[A M + d t M]) k = n i=1 (1 -k) a i (S i + t) k + k (b i + ω i + t a i ) (S i + t) k-1 ,
and the proof is complete.

Let us emphasize that we obtain exactly the same form of linear system as the one obtained in the classical formulation. However, two advantages are to be underlined. First, we have isolated the transport part from the flux part as already mentioned. Second, and it is clearly related, the abscissas appearing in the matrix to be solved do not evolve with transport in time but only change due to the fluxes so that the conditioning of the matrix, for scaled abscissas, can be controlled easily [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF].

Finally, it is important to realize that we still have a differential form of the DQMOM system of equations, a remark we already mentioned earlier as far as the moment equations were concerned. In the situation of a zero upper flux and with S min = 0, we see that the flux will only influence the zeroth order moment. The difference between the two approaches can be seen to be the same one as the difference between the kinetic schemes (for which we work with a time integrated version of the dynamical system and invoke the kinetic evolution in order to evaluate the fluxes) and the method of lines (for which the integration in only performed in the size variable and where the time is kept a continuous variable, thus leading to a system of ordinary differential equation to be resolved by a standard ODE integrator) in the framework of hyperbolic system of conservation laws [START_REF] Leveque | Numerical methods for conservation laws[END_REF][START_REF]Finite volume methods for hyperbolic problems[END_REF].

Relation between the two formulations.

Instead of working with a differential form of the system of abcissas and weights and plugging it in a LSODE like solver for ordinary differential equations, and encountering difficulties as far as stability is concerned, we make the link between the previous DQMOM formulation and the original integral form of the system developed in the previous Section in the following Lemma.

Lemma 4.2. Remembering the DQMOM new formulation in differential form (4.7), one can find a first integral for each k :

t 0 n i=1 (1 -k) a i (t ′ )S * k i (t ′ ) + n i=1 k b * i (t ′ ) (S * k-1 i (t ′ )) dt ′ = t 0 d t ′ n i=1 ω i (t ′ ) (S * i (t ′ )) k dt ′ Proof. Let us define B = k (S * i ) k-1 d t ′ (ω i S * i ) + (1 -k) d t ′ [ω i (S * i ) k ].
It is sufficient to note that :

B = k (S * i ) k-1 [d t ′ (ω i S * i ) -S * i d t ′ ω i ] + (S * i ) k d t ′ ω i = k (S * i ) k-1 ω i d t ′ S * i + (S * i ) k d t ′ ω i = ω i d t ′ (S * i ) k + (S * i ) k d t ′ ω i
which completes the proof. Thus, we have in fact an exact solution in time of the modified DQMOM differential system (4.5) through equation (3.1) once the fluxes Ψ ± have been evaluated. The link with a DQMOM formulation is then complete and we are also able to fully compare the two approaches and appreciate the novelty of the proposed numerical scheme.

5. Evaporation rate depending on size, using a piecewise constant approximation. In this Section, our aim is to treat a general evaporation law of the form

∂ t f -∂ S (R f ) = 0, (5.1)
with a piecewise constant coefficient R(S) which can approximate a continuous but non constant law. This appears naturally when considering more realistic evaporation law like in [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], which continuously depends on droplet temperature; this temperature, continuously depending of the droplet size because of the thermal inertia, is approximated by a piecewise constant function in multi-fluid method developed there, in such a way that the evaporation rate is then also a piecewise constant function. Our goal is to show how this phenomena can be simulated by the algorithm developed in this article.

In order to get a numerical solution based on the evolution of the moments, without lost of completeness, we will consider two size sections [0, 0.5] and [0.5, 1], in which the evaporation is taken constant. So far, the method based on moment dynamics has been displayed for a case where the evaporation rate was constant all over the size interval. The algorithm presented in § 3.4 must be modified in order to preserve the moment space. Indeed, in step 3, the abscissas, intially in [S k-1 +t, S k +t], are shifted by -t so that they lie in [S k-1 , S k ] at the end of the algorithm. Let us take in the first section R = -1 and in the second section R = -2. With the previous algorithm, the size abscissas of the first section would then be in [S 0 + t, S 1 + 2t], so that at then end they lie in [S 0 , S 1 + t]. There is no guarantee that the corresponding vector of moment in [S 0 , S 1 + t] is a vector of moment in [S 0 , S 1 ].

In order to devise the algorithm to solve that problem, we do an analogy with the pressureless gaz. Let consider a NDF f (t, S, v) depending of the evaporation velocity v. It is solution of the following equation :

∂ t f + v∂ S f = 0.
We then state that f (t, S, v) = n(t, S)δ(R(t, S)-v), in such a way that R is a velocity in the size phase space. If we take the 0 th and 1 st order moment with respect to R of the kinetic transport equation, we get

∂ t n(t, S) + ∂ S n(t, S)R = 0 ∂ t n(t, S)R + ∂ S n(t, S)R 2 = 0 (5.2)
With the intial condition that R is discontinuous, this appears to be a Riemann problem for the density n and the evaporation rate R (which is also a velocity). The scheme we use consists in solving a Riemann problem between each size interval and for each time step. Let us notice that the second equation of the system (5.2) is not an equation of the initial problem. The evaporation rate just depends on the size and is constant over time. The reason why we artificially transport R during a time step is because the piecewiese approximation on R generate shocks in the scheme. If R is not transported, then a shock occurs at the limit between the sections so that it is impossible to know in which section it belongs to. Moreover, it is compatible with a Lagrangian view of the droplets which goes from a section to another. Figure 5.1 displays the solution in the sense of the characteristics that is computed. The first step of the algorithm is exactly the same as for the one described in § 3.4. But contrary to the previous algorithm, we don't have to calculate two fluxes any more (Ψ k (t) and Ψ k+1 (t)), but five:

1. the evaporative flux going out of section 1 ( corresponding to Ψ k (t)) 2. the flux of the droplets from section 1 whose characteristics go into the δ-shock 3. the flux resulting from the droplets from section 1 which are in the δ-shock 4. the flux of the droplets from section 1 whose characteristics go into the δ-shock 5. the evaporative flux entering section 1 but not going into the δ-shock ( corresponding to Ψ k+1 (t)) By adding theses fluxes computing the corresponding weights and abscissas and shifting them at the correct evaporation rate, we are able to construct the solution of the Riemann problem. The velocity of the δ-shock is given in [START_REF] Bouchut | On zero pressure gas dynamics[END_REF]

u =        R + (0.5)( √ f + (0.5)+R -(0.5) √ f -(0.5) √ f + (0.5)+ √ f -(0.5) if f -(0.5) = f + (0.5) R + (0.5)+R -(0.5) 2 if f -(0.5) = f + (0.5) (5.3)
6. Treatment of arbitrary evaporation laws. Here, our aim is to treat a general evaporation law of the form

∂ t f -∂ ξ (R f ) = 0, (6.1)
for which the form of the evaporation coefficient is given as a function of ξ in general form R = R(ξ). Our aim in this Section is to directly generalize the algorithm defined in § 3.4 to this kind of arbitrary but smooth evaporation laws in a sense that will be defined in the following.

First, we will evaluate the characteristics and show that some restrictions have to be imposed on the profile of function R from both the modeling and the mathematical point of view. The key issue is the behavior near zero size. Let us make the following set of assumptions on the evaporation law. We assume that we deal with real evaporation, that is, R(ξ) > 0 for ξ > 0 and since droplets/particles of negative size do not evaporate, we assume R(ξ) = 0, ξ < 0. Consequently, there are two cases. The first one is the case of a discontinuous evaporation coefficient, where we assume that R(ξ) is globally Lipschitz on the interval [0, S max ] and admits a limit at zero size which is positive, thus yielding measure solution with a measure concentration of droplet number density at zero size. The second case, where R(ξ) is globally Lipschitz on every compact subset of the interval (0, ξ max ], includes the continuous one, where R(0) = 0, as well as cases with infinite values at zero size (the typical example is the d 2 -law with (6.1) written for the radius of the droplets ξ, in such a way that R(ξ) = cst 8πξ ). 6.1. Change of variable -characteristics. In this context, let us consider the characteristic curves defined by

dΞ dt = R(Ξ). (6.2)
For any ξ 0 > 0, with the considered regularity of R, they are defined for all t as long as the curve does not crosses ξ 0 . We then define γ(ξ) a primitive of the function 1/R(ξ); it is a bijection of [ξ 0 , ξ max ] onto its range. We also define Ψ(γ), a function such that f (t, ξ)dξ and Ψ(t, γ)dγ define the same measure on their respective intervals; this implies that Ψ(t, γ) = R(ξ) f (t, ξ) and

∂ t Ψ -∂ γ Ψ = 0. (6.3)
We then obtain the usual transport equation (1.4). However, two key questions arise, the first one concerns the behavior of the evaporation law around zero size and the second, related to the first, do the droplet disappear in finite time due to evaporation, or do their size approach zero asymptotically?

It is clear that the disappearance in finite time of the droplet occurs if the limit of R as ξ approaches zero is positive and even infinite. Assuming that the evaporation law behaves like a ξ α in the neighborhood of zero, the limiting value of α is 1. For α = 1, the size of one droplet approaches zero asymptotically and exponentially. For any positive ε, α = 1+ε leads also to an asymptotic decrease towards zero, nevertheless algebraically, whereas α = 1ε leads to a disappearance of droplets in finite time. This kind of behavior was observed already in [START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF]. However, the problem of the flux of droplet number at zero size has only to be treated when α = 1ε. Consequently we restrict our study to this case. In such a case we can take γ as the primitive of 1/R which is zero at ξ = 0; this yields a one to one correspondence between ξ and the characteristic variable γ, for ξ ∈ [0, ξ max ]. Finally, the theory presented before can easily be extended to the characteristic system in the Ψ variable instead of f . Examples : R(ξ) = c/ξ, γ = ξ 2 /2, this is the case of the droplet radius with c = 1/8π and γ is the droplet surface. ε = 1/2, a = 1, Ξ(t) = (Ξ(0) 1/2t/2) 2 which leads to an evaporation time of 2 if the initial size is 1, γ = 2ξ 1/2 . 6.2. Algorithm of the new scheme. The algorithm presented in § 3.4 can be generalized to arbitrary evaporations laws, still assuring the preservation of the moment space, in the following way : a discretization 0 = ξ 0 < ξ 1 < . . . < ξ N = ξ max of the support [0, ξ max ] is introduced. Let us assume that the moments M k (t) on each interval [ξ k-1 , ξ k [ and at time t are known. The moments at time t + ∆t is computed by : 1. In each section [ξ k-1 , ξ k ], a reconstruction f k is done of the distribution from its moments M k (t) by the ME method. Let then define the characteristics Ξ k by :

dΞ k dt = R(Ξ k ), Ξ k (t) = ξ k . (6.4)
A ODE solver can be used to obtain Ξ k (t + ∆t). The flux Ψ k -(t + ∆t) with the section k -1 is computed by :

Ψ -(t + ∆t) = Ξ k-1 (t+∆t) ξ k-1 f k (β)      1 β . . . β N     
dβ. (6.5) 2. The weights ω i,k and the abscissas S i,k corresponding to the moments M k (t)-Ψ k -(t + ∆t) + Ψ k+1 -(t + ∆t) are computed using the QD algorithm. Since these moments are in fact the moments of

i f i (ξ)1 [ξi-1,ξi[ (S) on [Ξ k-1 (t + ∆t), Ξ k (t + ∆t)], the abscissas are in ]Ξ k-1 (t + ∆t), Ξ k (t + ∆t)[. 3. Let then define the characteristics Y i,k by : dY i,k dt = -R(Y i,k ), Y i,k (0) = S k i . (6.6)
The moments M k (t) corresponding to the weights ω i,k and the abscissas Y i,k (∆t) are computed. The abscissas Y i,k (∆t) necessarily belong to the interval ]ξ k-1 , ξ k [, so that M k (t) belongs to the moment space corresponding to this interval. It gives then an approximation of the moments of f (t + ∆t, .) on ]ξ k-1 , ξ k [. This algorithm is valid for a CFL-like condition max ξ (R(ξ))∆t < ∆ξ. In fact, it is still valid for a less restrictive condition :

Ξ k (τ ) ∈ [ξ k , ξ k+1 ] for τ ∈ [t, t + ∆t].

Results.

7.1. Constant evaporation law. The method described in § 3 is tested on three kinds of initial distributions : a smooth one, presented in Fig. 7.1, a discontinuous one, presented in Fig. 7.4 and a distribution corresponding to moments close to the boundary of the moment space M 2 in Fig. 7.7 (it is the ME-reconstruction corresponding to canonical moments p 1 = 0.8, p 2 = 0.01 and p 3 = 0.5, see Appendix B).

A discretization into sections is introduced and the moments inside each section are computed for the initial distribution in order to initialize the simulation with our moment method. The evolution of the obtained moments can be compared with the moments of the exact solution (1.5). Moreover, since, at each time step a ME reconstruction f ME (S) of the distribution is performed, we can easily compare such an approximated distribution from the data of moments to the exact distribution, even if the computed variables are the sole moments. Finally, the variable of interest for the spray is its mass density, which corresponds to a moments of order 3/2. This fractional moment is not resolved but can be calculated through the ME reconstruction of the distribution by Smax 0 S 3/2 f ME (S)dS. This variable will be compared with the moment of order 3/2 of the exact solution but also with the one obtained with another numerical method : a multi-fluid method of order two derived in [START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF]. Smooth NDF (solid), ME reconstructed NDF obtained with the moments computed with the moment method with 1 section and 4 moments (dashed line) and corresponding peaks, at times t = 0, t = 0.25, t = 0.5 and t = 0.75. error on the mass (in percent) moment method, 1 section moment method, 2 sections order 2 method, 4 sections order 2 method, 8 sections Fig. 7.3. Evolution of the error on the mass density relatively to its initial value, calculated with the moment method with one or two sections and the multi-fluid method of order 2 with 4 or 8 sections.

In Fig. 7.1, one can see the evolution of the smooth distribution. A simulation is also done with the scheme previously described with only 1 section and 4 moments. The time step is here dt = 0.025 but it can be seen that the accuracy results do not depend of the time discretization. The ME reconstruction as well as the peaks showing the abscissas and the weights corresponding to the moments are drawn. One can see that the ME reconstruction follows very well the exact distribution, with a total evaporation at approximatively the same time (see the scales of the figures). Moreover, the error on the moments is smaller than one percent of each initial moment as shown in Fig. 7.2-left. This error decreases under 0.2% of each initial moment if two sections [0, 1/2] and [1/2, 1] are used (see Fig. 7.2-right). With these two sections, the ME reconstruction can almost not be distinguished from the exact solution. Concerning the mass density, the error is smaller than 0.25% for a simulation with the moment method with 1 section and 4 moments and still ten times smaller with 2 sections, as seen in Fig. 7.3. In order to have the same level of accuracy, four sections are needed with the multi-fluid method of order two, that is to say twice more moments. Moreover, the multi-fluid method giving directly equations on moments of order 3/2 on each section, developed in [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF], were also tested, using a constant reconstruction on each section. It can be seen that our new scheme is much accurate than this method, using a number of moments three times smaller (like in Fig. 2.1-left).

In Fig. 7.4, one can see the evolution of the discontinuous distribution. A simulation is also done with the moment method with only 2 sections and 4 moments per section. One can see that the method is still able to reproduce a good estimation of this kind of distribution but with a lower level of accuracy. However, the accuracy of the moment is good as shown in Fig. 7.5 for calculations with one or two sections. In fact, except for the moment of order 0, the accuracy is even very good. But the difference between the two discretizations is not as large as for the regular case. Concerning the mass density, the error is smaller than 0.3% for a simulation with the moment method with 1 section and 4 moments and two times smaller with 2 sections, as seen in Fig. 7.3. In order to have the same level of accuracy, six sections are needed with the multi-fluid method of order two in the first case and twelve in the second case, that is to say a number of moments three of two times larger.

In Fig. 7.7, one can see the evolution of the distribution corresponding to moments close to the boundary of the moment space M 2 . The initial distribution is chosen as the ME-reconstruction corresponding to canonical moments p 1 = 0.8, p 2 = 0.01 and error on the moments (in percent) moment of order 0 moment of order 1 moment of order 2 moment of order 3 Fig. 7.5. Evolution of the error on the moments calculated with the moment method with one section (left) or two sections (right) relatively to their initial value. p 3 = 0.5, as done in Appendix B. In this case, the computation of integrals < x k > ξ , Eq. (2.6), as well as the flux (3.6) are no more done with a Gaussian quadrature, which is not precise enough. Instead, Simpson formula are used here on 1000 sub-intervals of the integration interval; it is too much, but the aim was just to test the method.

The results are very good, since the analytical evolution of the initial distribution and the ME-reconstruction of the moments obtained by the computation can not be 7.2. Discontinuous evaporation law. We assume here that the evaporation rate has a linear profile and writes R = R 0 (1 + αS), with R 0 = -0.5 and α = 2. We show now a comparison between the analytical solution and the numerical solution of Eq. (6.1). Using the change of variable proposed in § 6, the solution writes :

f (t, S) = R(γ -1 (γ(S) + t)) R(S) f 0 (t, γ 1 (γ(S) + t) (7.1)
where

γ(S) = S 1 dS ′ R(S ′ ) = 1 αR 0 log 1 + αS 1 + α (7.2)
as we consider the size interval [0, 1].

We approximate R by its averaged values on each interval : -0.75 on the first interval and -1.25 on the second one. Simulation is done with a time step ∆t = 0.025. Figure 7.9 displays the analytical and numerical number density at different times. The ME reconstruction as well as the peaks showing the abscissas and the weights are drawn. One can see that the ME reconstruction describes well the exact distribution. Figure 7.10 displays the errors on the moments, compared to the analytical solution with the discontinuous evaporation rate or to the one with the affine evaporation rate. The errors compared to the affine evaporation rate are smaller than 8% of each initial moment; moreover, this error is smaller than 1% for moments of order greater than one. The accuracy is lower than in the case of a uniform evaporation rate, but this is a very difficult case because the rate R has a large variation and is approximated by piecewise constant function only on 2 sections. This result is quite encouraging and validates our method. 7.3. Arbitrary evaporation laws. Algorithm of Section 6.2 is tested here for different kind of evaporation laws, with always a smooth initial NDF, presented in Fig. 7.1.

7.

3.1. Smooth evaporation law. The algorithm described in the previous Section is tested on two kinds of coefficient of evaporation : an affine one R(ξ) = 0.5 + ξ and a "square root" one R(ξ) = √ 0.5 + ξ. Let us note f 0 the initial NDF. Using the characteristics, it can be shown that the solution of Eq. (6.1) is then given by : Fig. 7.9. Smooth NDF (solid), ME reconstructed NDF obtained with the moments computed with the moment method with discontinuous evaporation rates with 2 sections and 4 moments (dashed line), and corresponding peaks, at times t = 0, t = 0.17, t = 0.33, t = 0.5, t = 0.67 and t = 0.83 where Ξ(t; 0, ξ) is the solution at time t of Eq. (6.2) with Ξ(0) = ξ. For the affine coefficient of evaporation, it is given by : Ξ(t; 0, ξ) = (0.5 + ξ) exp(t) -0.5 (7.3) and for the square The initial NDF is the smooth one, presented in Fig. 7.1.

f (t, ξ) = f 0 (Ξ(t; 0, ξ) R(Ξ(t; 0, ξ) R(ξ) ,
The moment method is then applied with only one section and four moments. In practice, a second order explicit Runge-Kutta method is used to solve the ODEs (6.4) and (6.6), the explicit Euler method inducing a double value of the error. Moreover, an analytic solution of such equations, given here by Eq. ( 7.3) or Eq. (7.4) does not bring more accuracy here. The important property to be verified by the ODE solver is that it is accurate enough to avoid characteristic crossing (which will by always the case with a monotone function R).

The evolution of the error on the 4 considered moments are plotted in Fig. 7.11-left for the affine coefficient of evaporation and in Fig. 7.11-right for the square root one. It shows a level of accuracy similar to the one obtained with the constant evaporation law in § 7.1.

Note : An approximation of the integral of Eq. (6.5) by

(Ξ k-1 (t + ∆t) -ξ k-1 )f k (β)      1 β . . . β N      , β = Ξ k-1 (t + ∆t) + ξ k-1 2
could be envisioned since it is less costly that the computation of this integral by the double-precision 24-point Gaussian quadrature method. But it is not accurate enough to allow the conservation of the moment space, since it is like approximating f k by an affine function, which can be negative.

7.3.2. Discontinuous evaporation law. Like in §7.2, a piecewise constant evaporation law can be used as an approximation of the affine coefficient of evaporation from the previous Section. The algorithm of Section 6.2 can then be applied directly. The precision of such computation, with 2 sections and 4 moments in each section, is presented in Fig. 7.12-left. The error is smaller of the percent. However, as seen in Fig. 7.12-right, the accuracy of this kind of approximation of the affine coefficient of evaporation by a piecewise constant one induces an error of 10% on the number and less than 2% for the other moments. It induces the same level of accuracy between the simulation and the analytical result with the affine coefficient. error on the moments (in percent) moment of order 0 moment of order 1 moment of order 2 moment of order 3 Fig. 7.12. Evolution of the error on the moments calculated with the moment method with one section for the discontinuous coefficient of evaporation relatively to their initial value (left) and absolute value of the difference between the moments obtained analytically with the affine coefficient of evaporation and the moments obtained analytically with the discontinuous one, relatively to their initial value (right). 7.3.3. Regularity and "natural" variables. A more general formulation of the kinetic equation can be written using a variable φ for the droplets size. The NDF is then f φ (t, x, φ, u). The parameter φ denotes for example the volume V as in [START_REF] Greenberg | On the origin of spray sectional conservation equations[END_REF], the surface S as in [START_REF] Domelevo | The kinetic sectional approach for noncolliding evaporating sprays[END_REF] or the radius R (f R dR = f S dS = f V dV ). Let us assume that φ = S β . In the case of a d 2 evaporation law for a purely evaporative case, the simple solution f S (t, S) = f S (0, S + t) + n 0 (t)δ(S) becomes for the variable φ :

f φ (t, φ) = (φ 1/β + t) β-1 φ 1-1/β f φ (0, (φ 1/β + t) β ) + n 0 (t)δ(φ).
This leads to a singularity as soon as β > 1, like for the choice of the volume as the size variable.

For the choice of the radius as the size variable, the generalized algorithm described in § 6.2 can be directly applied. The coefficient of evaporation is then R(R) = 1/(8πR) and the CFL-like condition can be written ∆t < 4π(∆R) 2 . The evolution of the corresponding moments are plotted in Fig. 7.13 (from the biggest to the lowest, there are the moments of order 0, 1, 2 and 3). It can be seen that computation with one section (but with the exact resolution of the ODEs for the characteristics) is not very accurate. In fact, as seen on Fig. 7.14 it is necessary to use at least 3 sections to have the same level of accuracy on the number and the mass (moments of order 0 and 3, here) as for the resolution on the equation on the surface, § 7.1. In Fig. 7.15, the evolutions of corresponding analytical NDF and reconstructed NDF are plotted, illustrating the loss of accuracy of the computation with one section, when time increases. The conclusion of this test case is that it is more efficient to use the "natural" variable S than other variables like the radius.

8.

Conclusions. An extension of Eulerian multi-fluid method is given here in a purely evaporating case, using ideas from DQMOM Methods and considering at least 4 moments of the NDF in each section instead of one in [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] or two in [START_REF]Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF][START_REF] Dufour | A second-order multi-fluid model for evaporating sprays[END_REF]. Two difficulties have been resolved. The first one is the evaluation of disappearance rate of evaporating spray which is represented by a pointwise value of the NDF whereas the considered variables are its moments. The second one is the derivation of a scheme insuring the preserving of the moment space. This is done through a kinetic scheme whereas classical ODE solvers are unstable. The resulting scheme is shown to be very precise, with a total number of variables smaller than for the previous methods.

This work is inserted in the development of general scheme for the complete NDF equation, physical transport and transport in space phase being treated separately through splitting method [START_REF] Massot | Eulerian Multi-fluid method for the numerical simulation of evaporating polydisperse sprays : modelling and numerical issues in multi-dimensional configurations[END_REF][START_REF] De Chaisemartin | Evaluation of eulerian multi-fluid versus lagrangian methods for ejection of polydisperse evaporating sprays by vortices[END_REF][START_REF] Massot | Eulerian multi-fluid models : modeling and numerical methods[END_REF]. The total number of variable is then important for the total cost of the method since these variables have to be transported in the physical space [START_REF]Modeling of polydisperse sprays using a high order size moment method for the numerical simulation of advection and evaporation[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Massot | Eulerian multi-fluid models : modeling and numerical methods[END_REF].

Moreover, some extensions of our new scheme are straightforward : the extension to arbitrary evaporation laws (see § 6) or an evaporation law depending on droplet size (see § 5) represented by a discretization of this law inside sections. The ideas developed in this paper can easily be extended to two internal coordinates for example for soots, the only limiting point being the ability of the ME reconstruction to work in such configuration; however, this has already been treated in the literature and does not provide additional difficulties [START_REF] Abramov | A practical computational framework for the multidimensional momentconstrained maximum entropy principle[END_REF] except that the study at the boundary of the moment space should be updated to this new configuration and in the context of Analytical solution (solid lines), its reconstruction from its moments by maximization of entropy (dotted lines) and solution of (A.1) computed with the Runge-Kutta method for the singular initial distribution, at t = 0.6.

A.2. Explicit Euler method. Let us use the explicit Euler method with a fixed time step. We show that, as for the Runge-Kutta solver, the resolution breaks down, and we study now the cause of this collapse. Indeed, we see what happens for the canonical moments. As seen in § 2, the canonical moments indicate the position of the moment vector relative to the moment space boundaries. If the vector goes out of the moment space, that means that there is no function which successive moments are equal to the components of the vector, and thus no more ME reconstruction is available. The canonical moments must stay between 0 and 1.

In order to show that the resolution does lead to an improper solution regardless of the time step, we present two results obtained with three different time steps : dt = 0.1 and dt = 0.001. The evolution of the canonical moments is represented to show that the collapse is due to the (negative) value of one of them. For the smooth distribution, it can be seen in Fig. A.3 the canonical moments of the analytical solution. The moment space becomes narrower and narrower as time goes on. This is due to the fact that support of the distribution decreases as time goes on. For dt = 0.1, we can notice that the resolution collapses when the third canonical moment becomes negative (Fig. A.3-left). The error made by considering the evaporative flux constant during a whole time step is too big to preserve the moment space. When we consider a time step equal to 0.001, the stop time is bigger than 1 (Fig. A.3-right), because contrary to the analytical solution where all the droplets have been evaporated, there are still droplets left in the Explicit Euler solution because of the numerical diffusion. The Explicit Euler solution sticks very well to the analytical solution. That is not amazing because the smaller the time step is, the more accurate the Explicit Euler solution is, because the error made on the evaporative flux is smaller. It must be noticed that this result is possible because the initial distribution we consider has its canonical moments quite far from the moment space boundaries. So, even if there is a small error, the moment space will be preserved for a longer time. But problems appear when there is almost no droplets left. Indeed, the moment space becomes always narrower, and there is a time when the moment space is not preserved, which is revealed by the fact that at least one canonical moment becomes negative (Fig. A.3-right). Nevertheless, one can argue that this method gives precise results, before the moment space is no more preserved. However, its cost becomes prohibitive without insuring a global stability property. Thus the scheme breeds an intrinsic error which inevitably leads to the collapse of the resolution. We present in Fig. A.4 the canonical moments for the exact solution with the singular distribution. The conclusions are the same as the one for the smooth distribution. The smaller the time step is, the longer the space moment is preserved, but there is always a time when the errors make the moment space not preserved any more, as we can see for the time steps dt = 0.05 and dt = 0.001.

The direct resolution of system (A.1) with classical ODE solvers leads then to stability problems due to the fact that the moment space is not conserved. That is why kinetic scheme is developed in § 3. the moment of order 0 over an interval [0, S b ] with S b bigger than 0.02 of the two reconstructions only differ of less than 0.3% and are equal to 1 (the value for the corresponding point of the boundary) with an accuracy greater than 1%. Moments of order bigger than 1 over [0, S b ] are very small in two cases. In the context of a kinetic scheme for the resolution of the evaporation, like in § 3.4 or § 6.2, it is quite easy to see that, for moment vector close to the boundary (even close enough so that the Newton algorithm of the ME reconstruction does not converge) the flux can be computed with the lower principal representation, or eventually with the NDF corresponding to the point of the boundary. for two moment vectors close to these boundaries : they correspond to the canonical moments p 1 = 0.8, p 2 = 0.01 or p 2 = 0.99 and p 3 = 0.5. Like previously, the ME reconstructions and the lower principal representations give NDF close to the expected Dirac delta functions. Moreover, for the first case, p 2 = 0.01, the moments of the two reconstructions over an interval [0, S b ], with S b small enough are close to 0 with an accuracy greater than 2%. For the second case, p 2 = 0.99, the moments of order 0 of the two reconstructions over [0, S b ], with S b bigger than 0.03, only differ of less than 0.04% and differ of the value 1p 1 = 0.2 (the moment of order 0 of the NDF corresponding to the point af the boundary of M 2 ) of less than 0.04%. The other moments are close to zero with the same level of accuracy. Finally, in both cases, the flux can still be computed from the lower principal representation or eventually with the NDF corresponding to the moment vector at the boundary, with a good accuracy. And we can expect than the accuracy will increase when the moment vector will be closer of the boundary.

Fig. 2 . 1 .

 21 Fig. 2.1. Smooth NDF (solid line) and the constant reconstruction using its moments of order 3/2 on 12 sections (top-left figure, dashed line) or the bi-affine reconstruction using its 2 first moments on 4 sections (top-right figure, dashed line) or ME reconstruction using its 4 first moments (bottom figure, dashed line).

2 Fig. 2 . 2 .

 222 Fig. 2.2. Cut of the moment space in the m 2 -m 3 plan (bounded by solid lines) and cut of the moment space numerically attained with an error smaller than 1%, with the ME reconstruction using the Gaussian quadrature for the computation of the integrals (bounded by dashed lines) for m 1 = 0.1 (top-left figure), m 1 = 0.5 (top-right figure) and m 1 = 0.9 (bottom-right figure); (bottom-left figure) projection of the moment space boundary in the m 1 -m 2 plan (solid lines) and representation in this plan of the cuts (dashed lines).

Fig. 5 . 1 .

 51 Fig. 5.1. Numerical solution in the sense of characteristics. The characteristics of section 1 are represented by dashed lines and the ones of section 2 by dotted lines. The solid line stands for the δ-shock. The characteristics close to the border enter into the shock.

  Fig.7.1. Smooth NDF (solid), ME reconstructed NDF obtained with the moments computed with the moment method with 1 section and 4 moments (dashed line) and corresponding peaks, at times t = 0, t = 0.25, t = 0.5 and t = 0.75.

Fig. 7 . 2 .

 72 Fig. 7.2. Evolution of the error on the moments calculated with the moment method with one section (left) or two sections (right) relatively to their initial value.

Fig. 7 . 4 .

 74 Fig.7.4. Discontinuous NDF (solid), ME reconstructed NDF obtained with the moments computed with the moment method with 2 sections and 4 moments (dashed line) and corresponding peaks, at times t = 0, t = 0.25, t = 0.5 and t = 0.75.

Fig. 7 . 6 .Fig. 7 . 7 .Fig. 7 . 8 .

 767778 Fig. 7.6. Evolution of the error on the mass density relatively to its initial value, calculated with the moment method with one or two sections and the multi-fluid method of order 2 with 6 or 8 sections.

Fig. 7 . 10 .

 710 Fig. 7.10. Evolution of the errors on the moments calculated with the moment method with discontinuous evaporation rates, with 2 sections and compared to the analytical solution with the discontinuous evaporation law in absolute value (left) or to the analytical solution with the affine evaporation law (right).

Fig. 7 . 11 .

 711 Fig. 7.11. Evolution of the error on the moments calculated with the moment method with one section for the affine coefficient of evaporation (left) or for the square root coefficient of evaporation (right) relatively to their initial value.

Fig. 7 .Fig. 7 . 14 .

 7714 Fig. 7.13. (left figure) Evolution of the exact moments for the equation on the radius (solid lines) and the moments calculated with the moment method with one section (dashed lines) or two sections (dotted lines); (right figure) Evolution of the error on the moments calculated with the moment method for the equation on the radius, with one section relatively to their initial value.

Fig. 7 . 15 .

 715 Fig.7.[START_REF] Fox | Numerical simulation of spray coalescence in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF]. Smooth NDF in radius (solid), ME reconstructed NDF obtained with the moments computed with the moment method for 1 section and 4 moments (dashed line) or for 2 sections and 4 moments (dotted lines) , at times t = 0, t = 0.25, t = 0.5 and t = 0.75.

  Fig. A.2. Analytical solution (solid lines), its reconstruction from its moments by maximization of entropy (dotted lines) and solution of (A.1) computed with the Runge-Kutta method for the singular initial distribution, at t = 0.6.

Fig. A. 3 .

 3 Fig. A.3. canonical moments of the Euler Explicit solution for dt = 0.1 and dt = 0.001 for the smooth initial distribution.

Fig. A. 4 .

 4 Fig. A.4. canonical moments of the Euler Explicit solution for dt = 0.05 and dt = 0.001 for the singular initial distribution.

Fig. B. 1 .

 1 Fig. B.1. ME reconstruction with 4 moments, for p 1 = 0.01 (left) or p 1 = 0.99 (right), p 2 = 0.5 and p 3 = 0.5 and the Dirac delta functions corresponding to the lower principal representation.

Fig. B. 2 .

 2 Fig. B.2. ME reconstruction with 4 moments, for p 1 = 0.2, p 2 = 0.01 (left) or p 2 = 0.99 (right) and p 3 = 0.5 and the Dirac delta functions corresponding to the lower principal representation.

Fig. B. 3 . 1 .Fig. B. 4 .

 314 Fig. B.3. ME reconstruction with 4 moments, for p 1 = 0.2, p 2 = 0.5 and p 3 = 0.01 (left) or p 3 = 0.99 (right) and the Dirac delta functions corresponding to the lower principal representation.

Acknowledgments. We would like to thank Professor R.O. Fox for several helpful discussions and also the referee 1 for his/her useful comments allowing a substantial improvement of the paper.

* This work was supported by a Young Investigator Award for M. Massot (ANR-05-JCJC-0013 -jéDYS -2005-2009) from ANR in France (National Research Foundation), by the PPF TAMARIS-HP from French Ministry of Research (2006-2010, coordination F. de Vuyst), by a Young Investigator Award (ACI Nouvelles Interfaces des Mathématiques 2003-2006) to M. Massot from the French Ministry of Research and by a CNRS financial support (PEPS "Projet Exploratoire Pluridisciplinaire" 2007-2008, from the ST2I and MPPU Departments of CNRS, coordination: A. Bourdon and F. Laurent) .

Appendix A. Resolution of the ODE system with classical methods. The purpose of this appendix is to show the instabilities of a resolution of ODE systems of type (1.7) with classical schemes and a pointwise reconstruction. Let then use the ME reconstruction with N = 3 and S min = 0, S max = 1. If M(t) = (m 0 (t), m 1 (t), m 2 (t), m 3 (t)), the system which must be solved is the following one :

We will thus compare the results given by some classical ODE solvers with the analytical solution of the equation (1.4). The functions that we choose as references are a smooth and a singular one given in Fig. 7.1 and Fig. 7.4 respectively.

A.1. Explicit Runge Kutta scheme. We use a solver based on an explicit Runge-Kutta (4,5) formula with an adaptative time step [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF]. We present the results of the resolution scheme at two different times in order to illustrate the amplification of the errors generated.

For the smooth distribution, we can see in ). At each time step, an error is made for the evaporative flux. This error amplifies as it can be seen at the two times t = 0.15 and t = 0.2, and finally the moment vector doesn't belong to the moment space any more. In other words, there is no function whose successive moments are equal to the components of the moment vector and it exemplifies the stability problem encountered in the moment method. The problem of the solver is not a problem of stiffness of the system, but a problem of preservation of the moment space, as it is explained in the next subsection ( § A.2). The same conclusions hold for the singular distribution. We present the situation in Fig. A.2 at t = 0.6. Indeed when we compute the solution for a longer time, the resolution collapses.

Appendix B. ME reconstruction and boundaries of the moment space.

The aim of this part is to show the behavior of the ME reconstruction close to the boundary of the moment space. In order to evaluate the distance to this boundary, we use the canonical moments presented in § 2.1

At the boundary of M 1 , the first canonical moment is equal to 0 or 1. The solution of the Hausdorff problem is then a Dirac delta function δ 0 if p 1 = 0 or δ 1 if p 1 = 1, for a moment of order 0 m 0 equal to 1. The ME reconstructions with 4 moments for two moment vectors close to these boundaries are plotted in Fig. B.1; they correspond to the canonical moments p 1 = 0.01 or p 1 = 0.99, p 2 = 0.5 and p 3 = 0.5. The lower principal representations are also plotted (with an arbitrary scale). We can see that the ME reconstructions, as well as the lower principal representations, give NDF close to the expected Dirac delta functions. It can also been shown that