Frédérique Laurent
email: frederique.laurent@ecp.fr

Alaric Sibra
email: alaric.sibra@gmail.com

François Doisneau
email: francois.doisneau@centraliens.net

Two-size moment multi

Keywords: AMS subject classifications: 35L03, 65M08, 65M12, 76T10 Polydisperse sprays, Moment methods, Multi-Fluid, Evaporation, Coalescence

come

Introduction

Two-phase flows constituted of a gaseous phase carrying a disperse condensed phase play a key role in many industrial and scientific applications e.g. spray combustion in Diesel engines or aeronautical combustion chambers, heterogeneous energetic materials in solid rocket motors, fluidized beds, etc. In all these applications the disperse phase is composed of particles/droplets of various sizes that can possibly evaporate, breakup, coalesce or aggregate, and have their own inertia and size-conditioned dynamics and heating. The needs of research and industrial applications and the availability of powerful computational means appeal for high-fidelity, robust and efficient descriptions of the disperse phase in such flows.

We choose to describe the dynamics of the disperse phase in a statistical sense using a kinetic approach because of its accuracy and flexibility. The disperse phase information is completely contained in the so-called Number Density Function (NDF). The NDF measures an ensemble average (over a given set of initial conditions) number of particles at a specific location in the phase space at a given time. The phase space is determined by the number of internal coordinates that describe the particle state: position, velocity, size, temperature, etc. These variables evolve due to physical phenomena: transport, drag force, evaporation, heating, collisions etc. which are accounted for through the Williams equation [START_REF] Williams | Spray combustion and atomization[END_REF]. Moreover, in many applications, the coordinate of the phase space that is the most essential to deal with is size, because the other internal coordinates, such as the velocity or the temperature, are strongly conditioned on it.

There are several strategies to solve this kinetic equation. Lagrangian-Monte-Carlo approach [START_REF] Amsden | Kiva II, a computer program for chemically reactive flows with sprays[END_REF][START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF] allows to approximate the NDF by a sample of discrete numerical parcels describing particles of various internal coordinates. It is called Direct Simulation Monte-Carlo method (DSMC) in [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF] and is generally considered to be the most accurate method for solving this type of equation; it is specially suited for Direct Numerical Simulations (DNS) on canonical configurations since it does not introduce any numerical diffusion. However, the number of parcels required to achieve a satisfactory statistical convergence comes to be high in 3D cases, especially for unsteady configurations, when the size distribution has to be well approximated in addition to the spatial repartition of the spray. To overcome this limitation, Eulerian methods offer a promising alternative. The main objective is then to describe both the size distribution and the velocity (and eventually the temperature) conditioned on size. The size distribution, for the spray as well as for non-inertial particles (aerosol, soots) can be modeled thanks to three types of methods: 1) the sectional methods [START_REF] Gelbard | Sectional representations for simulating aerosol dynamics[END_REF][START_REF] Tambour | A Lagrangian sectionnal approach for simulating droplet size distribution of vaporizing fuel in a turbulent jet[END_REF] introducing a discretization of the size variable into intervals called sections, 2) the moment methods (see for example [START_REF] Kah | High order moment method for polydisperse evaporating sprays with mesh movement: application to internal combustion engines[END_REF][START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF][START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF][START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF][START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF]), which consist in writing equations on some moments of the NDF, 3) the class methods, sampling the size variable (see citations of [START_REF] Doisneau | Eulerian modeling and simulation of polydisperse moderately dense coalescing spray flows with nanometric-to-inertial droplets : application to Solid Rocket Motors[END_REF] and [START_REF] Sibra | Modélisation et étude de l'évaporation et de la combustion de gouttes dans les moteurs à propergol solide par une approche eulérienne Multi-Fluide[END_REF]). The velocity and temperature can be easily considered by class methods. However, these methods suffer from their inability to tackle integral terms except at the cost of strong modal simplifications. For the other methods, moments are used to capture the velocity (or temperature) conditioned on size. For the spray, a sectional type method, called multi-fluid, was then rigorously derived from the kinetic model [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], considering only zero and first order moments in velocity. It leads to a good description of the size polydispersity and size-velocity correlations of evaporating sprays with small enough droplets, at a reduced cost compared to DSMC [START_REF] De Chaisemartin | Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporating sprays[END_REF], since deterministic methods of discretization in the spatial dimension are affordable and efficient. For bigger droplets, higher order moments have to be considered [START_REF] Kah | Eulerian quadrature-based moment models for polydisperse evaporating sprays[END_REF], such as the velocity dispersion, for example using a Gaussian closure [START_REF] Dupif | Modélisation et simulation de l'écoulement diphasique dans les moteurs à propergol solide par une approche eulérienne polydispersée en taille et en vitesse[END_REF][START_REF] Sabat | Modèles Eulériens et méthodes numériques pour la description des sprays polydispersés turbulents[END_REF][START_REF] Sabat | On the development of high order realizable schemes for the Eulerian simulation of disperse phase flows: a convex-state preserving Discontinuous Galerkin method[END_REF][START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF]. This multi-fluid methods can reproduce accurately the size distribution, as well as the size-conditioned velocity if a large enough number of sections are used, but they remain usually too costly for practical applications. The moments methods allows the use of a smaller number of variables, a closure being often done by reconstructing a NDF from these moments, but it does usually not allow to accurately solve the NDF. Moreover, this induces complications when there are several internal variables and realizable numerical schemes are difficult to design, i.e. schemes ensuring that the variables are always moments of a positive measure. Still some satisfactory solutions are given for particular cases in [START_REF] Emre | Eulerian modeling of a polydisperse evaporating spray under realistic internal-combustion-engine conditions[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF].

To improve multi-fluid methods, several size moments can be considered in each section in order to be able to decrease the number of sections and then the overall cost of the method, for the same accuracy. It leads to some hybrid method between sectional and moment methods and the NDF is reconstructed inside the section. This reconstruction has to be efficient and accurate since it has to be done in each cell at each time step to close the equations. Such kind of hybrid method was developed in [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] for the description of the evaporation, using a polynomial reconstruction of the size NDF to close the equations, or a bi-affine one for the case of the two moment method, in order to ensure the positivity of this NDF. The interest of using two size moments inside the section, instead of an affine NDF reconstruction of MUSCL type from one moment per section, was also clearly shown on the accuracy of the method. However, only integer moments in the droplet surface variable were considered, in such a way that the mass density were only recovered thanks to the reconstruction. Moreover, the numerical analysis of the method were done only for the case the polynomial reconstruction in the section, not for the biaffine one, which is always positive. Another two size moment (TSM) method [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF][START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF] uses two important variables in each section: the number density and the mass density. To automatically ensure its positivity, an exponential function is used as the reconstruction in the sections, similarly to what is obtained with the entropy maximization strategy [START_REF] Mead | Maximum entropy in the problem of moments[END_REF]. It was used for for the simulation of dynamics and coalescence of particles in solid propellant combustion [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] but two limitations were exposed in that work: first the complexity and cost of the inversion algorithm and second its lack of accuracy -especially on steep and strongly varying size distributions-as a result of a numerical limitation to prevent overflows during inversion. More than two moments were also used for the description of evaporating dilute sprays, with the entropy maximization strategy for the reconstruction [START_REF] Kah | High order moment method for polydisperse evaporating sprays with mesh movement: application to internal combustion engines[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF]. This method showed its potential for the simulation of the injection of a liquid disperse phase into a gas for automotive engine applications. However, the more moments are used, the more realizability conditions has to be respected, demanding the development of realizable numerical methods and algorithms, which do not degrade the computational efficiency. Moreover, all the hybrid methods reviewed above use a constant reconstruction of the velocity in the section, thus limiting the description of size-velocity correlations inside the section, except [START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF], where a size-velocity moment is added, which brings in additional difficulties.

In this paper, the objective is to write a model on physically important variables in the sections, with a limited number of constraints on it (the realizability conditions) and with a closure thanks to a very simple, robust and accurate reconstruction. To limit the constraints, only the TSM methods are then considered. Moreover, a piecewise affine reconstruction for the size distribution, from the number and mass density, is considered, in the simplest way allowing to ensure its positivity. Since it is a keypoint both in terms of cost and accuracy of the method, this reconstruction is analyzed and compared to the exponential one. The velocity description inside the section is also improved, without considering any additional moment. Accurate and realizable schemes are then developed and validated for two categories of physical phenomena that are evaporation -a continuous drop size variation (derivative operator)-and coalescence -a discontinuous, non-linear, drop size variation (integral source term), in conjunction with the treatment of drag. This is done for each of such operator separately, in the perspective of the use of operator splitting for the complete problem and it comforts us on the fact that other phenomena can be treated as well (condensation, break-up, heating) but this is not discussed further for the sake of clarity and size of the paper. Some realizable transport schemes being given in [START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF] for structured meshes and in [START_REF] Le Touze | Multislope {MUSCL} method for general unstructured meshes[END_REF] for TSM on unstructured ones, this issue is no more discussed here. Finally, a comprehensive test case, exhibiting the size/velocity couplings [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] through drag, evaporation and coalescence, is performed with the chosen approach, allowing to assess its potential on realistic problems.

The paper is organized in the following way. In section 2 we introduce the multi-fluid model for a disperse phase. In section 3 we detail reconstruction techniques for the size and size-conditioned variables with a particular focus on numerical properties of accuracy, stability and cost. In section 4 we present and we systematically assess numerical schemes for the resolution of the evaporation part of the multi-fluid system. In section 5 a numerical scheme for the coalescence part of the multi-fluid system is devised and assessed. Finally in section 6 the full strategy is tested on a 2D unsteady case, which proves the efficiency of the modeling and numerical choices.

Eulerian spray modeling

Eulerian multi-fluid models can be derived from a kinetic model describing transport, drag, evaporation, heating, collisions and breakup [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF][START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[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][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. Here, it is done from a simplified kinetic model, extracting the main difficulties by considering the different categories of physical phenomena in the simplest way: transport, drag force inducing a size-velocity coupling, evaporation, which is an advection in the size phase space and coalescence, which is an integral source term. Moreover a particular attention will be devoted in this paper to ensure that the methods developed are applicable to more complex laws and to other phenomena, this application effort being done in another piece of work for complex drag, evaporation and heating laws [START_REF] Sibra | Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using eulerian multi-fluid methods[END_REF].

A kinetic description: the Williams equation

A dilute or moderately dense spray can be described by its number density function (NDF) f φ , where f φ (t,x,φ,u)dxdudφ denotes the average number of droplets (in a statistical sense) at time t, in a volume of size dx around a space location x, with a velocity in a du-neighborhood of u and with a size in a dφ-neighborhood of φ. The droplets are assumed to be spherical and characterized by one geometry parameter φ, which can be the radius R, the surface S of the volume V. All variables being considered as nondimensional, the relation between these size parameters are given by S = R 2 and V = R 3 . The NDF being the same in all three representations, we notice that f R dR= f S dS= f V dV.

For the sake of simplicity and for the purpose of this paper, we consider that the evaporation process is described by a d 2 law without convective corrections, that the drag force is given by a Stokes law, and finally that the unsteady heating of the droplets does not need to be modeled so that the droplet temperature was not considered here. We refer to [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF][START_REF] Sibra | Modélisation et étude de l'évaporation et de la combustion de gouttes dans les moteurs à propergol solide par une approche eulérienne Multi-Fluide[END_REF] for a consideration of heating.

The evolution of the NDF is described by the following Boltzmann-type equation [START_REF] Williams | Spray combustion and atomization[END_REF]:

∂ t f φ +∂ x •(u f φ )+∂ u •(F f φ )-∂ φ (K φ f φ ) = 1 Kn Γ( f φ , f φ ), (2.1)
where Kn is the Knudsen number. The Stokes drag force per unit mass is given by F = (u g -u)/St, where u g is the gas velocity and the Stokes number St is linear as a function of the droplet surface S(φ): St=A.S(φ). The evaporation rate K φ =K S /S (φ) is such that K S is a constant, denoted K in the following. The kinetic modeling for the collision operator is taken from [START_REF] Hylkema | A random particle method to simulate coalescence phenomena in dense liquid sprays[END_REF][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. Then, in our context of dilute to moderately dense sprays, only binary collisions are taken into account and the mean collision time is assumed very small compared to the inter-collision time. Moreover, the mass and the momentum of droplets are supposed to be preserved during the collision. Finally, we assume that every collision leads to coalescence. This last assumption is discussed in [START_REF] Achim | Simulation de collisions, coalescence et rupture de gouttes par une approche lagrangienne: application aux moteurs à propergol solide[END_REF][START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF] since droplet rebound or separation can happen if the remaining internal kinetic energy of the new droplet is too high. It brings in no limitation of this proposed framework. The collision operator resulting from these hypotheses is developed as Γ( f V , f V ) = Q + +Q -where Q + and Q - respectively correspond to the quadratic integral operators associated with creation and destruction of droplets due to coalescence [START_REF] Hylkema | Modélisation cinétique et simulation numérique d'un brouillard dense de gouttelettes[END_REF]:

Q -= - V * u * f V (t,x,V,u) f V (t,x,V * ,u * )β (V) 1 3 ,(V * ) 1 3 |u-u * |dV * du * Q + = 1 2 V * ∈[0,V] u * f V (t,x,V ,u ) f V (t,x,V * ,u * )β (V * ) 1 3 ,(V ) 1 3 |u -u * |JdV * du *
where, considering two precursor droplets of volumes V * and V colliding to form a new droplet of volume V, the mass and momentum conservation are written: V = V -V * and V u = Vu-V * u * . Moreover, J is the Jacobian of the mapping (V,u) → (V ,u ) and β(r * ,r ) = π (r * +r ) 2 is the impact parameter. A collision efficiency can be considered in the coalescence kernel [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] but is equal to one here, for the sake of simplicity.

This equation is not solved directly, but a system of equations on moments of f φ is derived from it. The formalism and the associated assumptions needed to derive the Eulerian multi-fluid models are introduced 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 in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] when including the coalescence. We recall the two main steps, which are the semi-kinetic and the sectional integrations in order to precisely introduce the new models considered in the present contribution.

Semi-kinetic model

In a first step, we reduce the phase space to the only droplet size variable. We then only consider the first two moments in the velocity variable:

n φ (t,x,φ) = u∈R f φ (t,x,u,φ)du, n φ (t,x,φ)u d (t,x,φ) = u∈R u f φ (t,x,u,φ)du. (2.
2)

The mono-kinetic assumption [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF] is used to close the equations, meaning that, at a given time and space location, the droplets having the same sizes have the same velocities, with no dispersion:

f φ (t,x,u,φ) = n φ (t,x,φ)δ(u-u d (t,x,φ)). (2.3)
This assumption is reasonable for dispersed sprays when the Stokes number St is small enough [START_REF] De Chaisemartin | Eulerian models and numerical simulation of turbulent dispersion for polydisperse evaporating sprays[END_REF]. For more inertial droplets, higher order velocity moments have to be considered [START_REF] Dupif | Modélisation et simulation de l'écoulement diphasique dans les moteurs à propergol solide par une approche eulérienne polydispersée en taille et en vitesse[END_REF][START_REF] Kah | Eulerian quadrature-based moment models for polydisperse evaporating sprays[END_REF][START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF], but this is not in the scope of this paper. For denser sprays, for which collisions can occur, this monokinetic assumption has been discussed in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. Indeed, the coalescence source term induces a velocity dispersion, since there is no reason for a droplet created by the coalescence of two droplets of various sizes, which is deduced from momentum conservation, to exactly match the velocity corresponding to its new size. In [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF], the Eulerian multi-fluid model was considered as a projection of the original distribution function at the kinetic level onto a 1D sub-manifold. The semi-kinetic model, giving equation on n φ and n φ u d is then obtained by assuming a Gaussian velocity distribution and taking the limit when the dispersion of this distribution tends to zero. This leads to the following system, written for φ = V for the sake of simplicity, dropping the V superscript in n V to simplify the notations:

∂ t n+∂ x (nu d )-∂ V (K V n) = 1 Kn -n(V) ∞ 0 n(V * )β(V,V * )|u d (V)-u d (V * )|dV * + 1 2 V 0 n(V-V * )n(V * )β(V-V * ,V * )|u d (V * )-u d (V-V * )|dV * , ∂ t (nu d )+∂ x (nu d ⊗u d )-∂ V (K V nu d ) = n(V) u g -u d (V) St(V) + 1 Kn -n(V)u d (V) ∞ 0 n(V * )β(V,V * )|u d (V)-u d (V * )|dV * + 1 2 V 0 n(V-V * )n(V * )β(V-V * ,V * ) (V-V * )u d (V-V * )+V * u d (V * ) V |u d (V * )-u d (V-V * )|dV * .

Multi-fluid models

The second modeling step leading to the Eulerian multi-fluid model is described in [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[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][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. It relies on the choice of a discretization 0 = φ 0 < φ 1 < ••• < φ N = φ max for the droplet size phase space, where φ max can eventually be equal to infinity. A system of conservation equations is then written for each fixed interval [φ k-1 ,φ k ] called section. The set of droplets in one section can then be seen as a "fluid" exchanging mass and momentum. In the original Eulerian multi-fluid model [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], only the mass density and the momentum corresponding to each section are considered, dimensionless here:

m k = φ k φ k-1 V(φ)n φ (φ)dφ, m k u k = φ k φ k-1 V(φ)u d (φ)n φ (φ)dφ, (2.4)
where V(φ) is the volume corresponding to the size variable:

V(R) = R 3 , V(S) = S 3/2 .
Only one size moment is then considered: the moment of order 3/2 if the size variable is the surface, as advised in [START_REF] Laurent | Analyse numérique d'une méthode multi-fluide Eulérienne pour la description de sprays qui s'évaporent[END_REF]. This method is referred to as the One Size Moment (OSM) method in the following. For a better description of the size distribution inside each section, some extension of the model using two moments were developed in [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF], adding the number density of the section, which is the moment of order 0:

n k = φ k φ k-1 n φ (φ)dφ. (2.5)
This model is referred to as Two Size Moment (TSM) model and it is the subject of the present paper. Since two moments are involved, a realizability constraint appears, so that they remain moments of a positive NDF on the considered interval:

0 < V(φ k-1 )n k < m k < V(φ k )n k , or n k = m k = 0. (2.6)
To close the equations for the OSM or TSM method, the size distribution and the velocity are reconstructed from the moments in the section:

n φ (φ) = N ∑ k=1 κ k (φ)1 [φ k-1 ,φ k ] (φ), u d (φ) = N ∑ k=1 χ k (φ)1 [φ k-1 ,φ k ] (φ). (2.7)
The resulting equations are then:

∂ t n k +∂ x (n k u k ) = N (k+ 1 2 ) -N (k-1 2 ) + 1 Kn Q n+ k -Q n- k , (2.8
)

∂ t m k +∂ x (m k u k ) = E (k+ 1 2 ) -E (k-1 2 ) -M (k) + 1 Kn Q m+ k -Q m- k , (2.9) ∂ t (m k u k )+∂ x (m k U (k) 2 ) = E (k+ 1 2 ) u -E (k-1 2 ) u -M (k) u +m k F (k) + 1 Kn Q mu+ k -Q mu- k , ( 2.10)
where the fluxes for the evaporation, with an upwind formulation, are written:

N (k-1 2 ) = Kκ k (φ k-1 ), E (k-1 2 ) = KV(φ k-1 )κ k (φ k-1 ), M (k) = φ k φ k-1 V (φ)Kκ k (φ)dφ, E (k-1 2 ) u = KV(φ k-1 )κ k (φ k-1 )χ k (φ k-1 ), M (k) u = φ k φ k-1 V (φ)Kκ k (φ)χ k (φ)dφ.
The transport term and the drag term are given by

m k U (k) 2 = φ k φ k-1 κ k (φ)χ k (φ)⊗χ k (φ)dφ, m k F (k) = φ k φ k-1 V(φ) u g -χ k (φ) St(φ) κ k (φ)dφ.
and the coalescence source terms are:

Q n+ k = 1 2 N ∑ i=1 N ∑ j=1 Q 1 ijk , Q n- k = N ∑ i=1 N ∑ j=1 Q 1 kij , Q m+ k = N ∑ i=1 N ∑ j=1 Q 2 * ijk , Q m- k = N ∑ i=1 N ∑ j=1 Q 2 * kij , Q mu+ k = N ∑ i=1 N ∑ j=1 Q 3 * ijk , Q mu- k = N ∑ i=1 N ∑ j=1 Q 3 * kij , with    Q 1 ijk Q 2 * ijk Q 3 * ijk    = D ijk χ j r 2 -χ j r 2   1 r 3 r 3 χ j r 2   κ i (r)κ j (r )β(r,r )drdr . ( 2.11)
where

D ijk = {(r * ,r ) ∈ [r i-1 ,r i ]×[r j-1 ,r j ], (r * ) 3 +(r ) 3 ∈ [r 3 k-1 ,r 3 k ]} [9]
. One obtains a system on m k and m k u k for OSM, with the constraint m k ≥ 0 and on n k , m k and m k u k for TSM under the realizability constraint (2.6). This system is closed as soon as the size and velocity reconstructions are given. It means that one has to define some κ k ≥ 0 and χ k from the considered variables with the constraints (2.4) and (2.5) for TSM. For OSM [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], constant reconstructions in the sections were used. The TSM leads to a more accurate description of the size distribution inside the section, depending on how the reconstruction is done. This key point is thoroughly analyzed in the next section and will condition the good behavior and the accuracy of the model. Moreover, a key improvement of the velocity reconstruction is proposed, without adding any moments, which will be essential for an accurate and robust description of the spray dynamics.

Size and velocity reconstruction from the moments

We first deal with the size reconstruction, which is the most sensitive part in term of robustness and accuracy. Both the exponential reconstruction [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] and a new affine reconstruction are considered. The analysis of the latter is then given in terms of stability and accuracy. The two reconstructions are then compared on specific test cases. Finally, two types of velocity reconstructions are described: constant or affine in the section.

In the following, the size variable is the droplet surface S since it was shown to be the best choice for purely evaporating cases (see [START_REF] Laurent | Analyse numérique d'une méthode multi-fluide Eulérienne pour la description de sprays qui s'évaporent[END_REF] and Appendix A). Indeed, the influence of the behavior of the NDF at zero size induced by this choice is important when considering evaporation, whereas it has no real impact on phenomena like coalescence.

Size reconstruction

We consider the discretization 0= S 0 < ••• < S N = S max of the size interval [0,S max [ and one denotes ∆S k = S k -S k-1 for each k ∈ {1,...,N} and ∆S = min k∈{1,...,N} ∆S k . The distribution κ k is here reconstructed from the moments of order 0 and 3/2, n k and m k , such that:

n k = S k S k-1 κ k (S)dS, m k = S k S k-1 S 3/2 κ k (S)dS.
The realizability condition (2.6) is then written, here:

0<S 3/2 k-1 n k <m k <S 3/2 k n k or n k =m k =0.
An exponential form of the reconstruction was first used [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF][START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF], automatically ensuring non-negativity for the reconstruction. But due to its computational drawbacks, we decide to introduce here a simple affine reconstruction in a way also ensuring non-negativity.

Exponential reconstruction

In [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF][START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF], a reconstruction of the distribution is chosen as an exponential in each section:

f a (S) = a k exp(-b k (S-S k-1 )), S ∈ [S k-1 ,S k ]. (3.1)
The values of the two parameters a k and b k are computed so that the moments of order 0 and 3/2 of f a on the section [S k-1 ,S k ] are exactly n k and m k . It can then be shown [START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF] that b k is the unique zero of the function b → g k (b)-

m k n k
, where g k is defined for b ∈ R by:

g k (b) = S k S k-1 S 3/2 e -bS dS S k S k-1 e -bS dS (3.2) = 1 1-e -b∆S k S 3/2 k-1 -S 3/2 k e -b∆S k + 3 2b S k-1 -S k e -b∆S k + 3e bS k-1 2b √ S k √ S k-1 e -br 2 dr
and a k is given by:

a k = b k 1-exp(-b k ∆S k ) n k . (3.3)
Numerically, Ridders' method [START_REF] Ridders | A new algorithm for computing a single root of a real continuous function[END_REF] is used to compute b k by solving the non-linear system

g k (b) = m k n k
. This function g k (b) is computed using the Dawson's integral for b < 0, the exponentially-scaled complementary error function for b > 0 and a Taylor expansion is used to spare cost for small values of |b|, a particular attention being devoted to the junctions. The function g k is then quite costly to compute and the resolution of the nonlinear problem induces usually about 7 to 17 † evaluations of g k to converge with an error smaller than 10 -4 . Moreover, the reconstruction's accuracy and stability are threatened by overflow for large values of |b| so that b k is searched in the interval -

X ∆S k , X ∆S k ,
where e X is the overflow constant. Moreover, errors are still encountered at the inversion step for large values of |b|. It corresponds to sections where the moment vector (n k ,m k ) is closed to the boundary of the moment space [START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF], i.e. where the mean volume

m k n k is closed to S 3/2 k-1 or S 3/2 k .
Indeed, the derivative of g k (b) is then very low, behaving as b -2 or even b -2.5 in the first section. This means that a small error on the ratio m k n k induces quite a large error on b k . And the error on the parameter b k induces an error of the same order times ∆S k , for the value of the reconstructed function at one bound of the section, since n k = O(∆S k ) and:

f a (S k-1 ) = b k 1-exp(-b k ∆S k ) n k , f a (S k ) = b k exp(b k ∆S k )-1 n k .
All that limits the accuracy of the method if steep distributions are handled. Concerning the cost, the reconstruction itself can be costly due to the large number of evaluation of a costly function. Moreover, computing source terms such as (2.11) can also be costly.

Since no analytical formula can be written, computation based on a quadrature of the measure κ k (S)dS was proposed in [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF], allowing a good accuracy, even when large |b k | are involved, but at the price of a computation of quadrature points at each time step and for each section. Then, because of its possible lack of accuracy and high cost, which will be experienced in section 3.1.5, (and also in sections 4.1.4 and 5.2 for the cost for the evaporation and coalescence operators), an alternative reconstruction is proposed.

Affine reconstruction

There is a great interest in using as simple functions as polynomial ones for the reconstruction, hopping then an easy computation of the source terms and also a simplification of the computation of this reconstruction. However, their positivity is no more ensured.

Here, for TSM, a simplified version of the positive bi-affine reconstruction of [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] is proposed, from the moment of order 0 and 3/2 whereas moments of order 0 and 1 were considered in [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF]. So we will first show how this kind of reconstruction can always be done. Then, some nice stability and accuracy properties are shown for this new piecewise affine reconstruction, in the context of the approximation of regular functions, which are essential for a good accuracy of our spray simulations. The reconstructed function f a is then taken continuous inside each section ]S k-1 ,S k [ and affine on the whole section or on a part of it, the function being zero elsewhere. This function can then be written:

f a (S) 6 Z Z Z Z Z Z α k S k-1 =S (k) a β k = 0 S (k) b S k f a (S) 6 α k β k S k-1 =S (k) a S k =S (k) b f a (S) 6 β k S k-1 α k = 0 S (k) a S k =S (k) b
f a (S) =        0 if S k-1 ≤ S ≤ S (k) a α k +(β k -α k ) S-S (k) a S (k) b -S (k) a if S (k) a ≤ S ≤ S (k) b 0 if S (k) b ≤ S ≤ S k i ∈ {1,...,N}, (3.4)
with only three considered cases illustrated in Figure 3:

[C.1] S (k) b ∈]S k-1 ,S k [, S (k)
a = S k-1 and β k = 0. [C.2] S (k) a = S k-1 and S (k) b = S k . [C.3] S (k) a ∈]S k-1 ,S k [, S (k)
b = S k and α k = 0. Let us denote R k-1 = S k-1 , R b = S (k) b and µ (k) in f = 2 ∆S 2 k S k S k-1 (S k -S)S 3/2 dS, µ (k) sup = 2 ∆S 2 k S k S k-1 (S-S k-1 )S 3/2 dS.
These two last scalars belong to ]S 3/2 k-1 ,S 3/2 k [. The values of the parameters of the reconstruction (3.4) are then given in the following proposition.

Proposition 3.1 (Piecewise affine reconstruction

). Let n k > 0 and m k ∈]S 3/2 k-1 n k ,S 3/2 k n k [, for some k ∈ {1,2,...,N}. Then, there exists on ]S k-1 ,S k [ a unique non-negative continuous function f a given by (3.4) with [C.1], [C.2] or [C.3] and such that S k S k-1 f a (S)dS = n k , S k S k-1 S 3/2 f a (S)dS = m k . (3.5)
The function's parameters are given by:

• If S 3/2 k-1 < m k n k < µ (k) in f , then S (k) a = S k-1 , β k = 0, α k = 2n k S (k) b -S k-1 and S (k) b is the root in ] S k-1 , √ S k [ of P k (X) =2X 5 +4R k-1 X 4 +6R 2 k-1 X 3 + 8R 3 k-1 - 35m k 4n k X 2 + 10R 3 k-1 - 35m k 2n k R k-1 X + 5R 3 k-1 - 35m k 4n k R 2 k-1 . (3.6) • If µ (k) in f ≤ m k n k ≤ µ (k) sup , then S (k) a = S k-1 , S (k)
b = S k and

α k = 2 µ (k) sup n k -m k ∆S k µ (k) sup -µ (k) in f , β k = 2 m k -µ (k) in f n k ∆S k µ (k) sup -µ (k) in f .
(3.7)

• If µ (k) sup < m k n k < S 3/2 k , then S (k) b = S k , α k = 0, β k = 2n k S k -S (k) a and S (k) a is the root in ] S k-1 , √ S k [ of Q k (X) =2X 5 +4R k X 4 +6R 2 k X 3 + 8R 3 k - 35m k 4n k X 2 + 10R 3 k - 35m k 2n k R k X + 5R 3 k - 35m k 4n k R 2 k . (3.8)
Proof. If such a function f a exists, then the constraints (3.5) can be written

n k = α k + β k 2 (S (k) b -S (k) a ), m k = α k S (k) b S (k) a S (k) b -S S (k) b -S (k) a S 3/2 dS+ β k S (k) b S (k) a S-S (k) a S (k) b -S (k) a S 3/2 dS.
In the case [C.2], the resolution of this system on α k and β k gives (3.7) and they are nonnegative if and only if µ

(k) in f ≤ m k n k ≤ µ (k) sup . In the case [C.1], S (k)
b has to be such that

m k n k = g(S (k) b ), g(x) = x S k-1 (x-S)S 3/2 dS x S k-1 (x-S)dS = 2 (x-S k-1 ) 2 x S k-1 (x-S)S 3/2 dS and the ratio m k n k is smaller than µ (k) in f = g(S k ) since the function g is an increasing function on ]S k-1 ,+∞[. Let us denote R b = S (k) b , R k-1 = S k-1 , R k = √ S k . Dividing by (R b - R k-1 ) 2 , it leads to: m k n k = 4 35(R b +R k-1 ) 2 2R b 5 +4R k-1 R b 4 +6R 2 k-1 R b 3 +8R 3 k-1 R b 2 +10R 4 k-1 R b +5R 5 k-1 .
It is then easy to show that R b should be a root of the polynomial function defined by (3.6). Moreover, let us remark that:

P k (R k-1 ) = 35R 2 k-1 R 3 k-1 - m k n k < 0, P k (R k ) = 35 4(R k -R k-1 ) 2 ∆S 2 k µ (k) in f - m k n k > 0.
Then P k has at least one root on ]S k-1 ,S k [. Moreover, using Descartes' rule of signs, we show that P k has at most one positive root since the coefficients are either positive or of the same sign as 32 35

R 3 k-1 -m k n k or to 20 35 R 3 k-1 -m k n k
for the two last ones. The last case is similar to the previous one. Indeed, for [C.3], one necessarily has

m k n k > µ (k) sup > 2 7 R 3 k and R a = S (k) a is a root of Q k . Moreover, Q k has only one root in ]R k-1 ,R k [ since Q k (R k-1 ) < 0, Q k (R k ) > 0 and thanks to Descartes' rule of signs. Let us denotes: S m,k = 1 2 35 m k n k 2/3 . One can remark that, in the case [C.1], P( S m,k ) is positive, meaning that S (k) b < S m,k
. This allows to treat easily the last section in the case S max = +∞, which is a great advantage over the OSM method.

Moreover, a root of a polynomial has to be computed in cases [C.1] and [C.3]. Numerically, the second order modified Mueller method [START_REF] Costabile | A modification of Muller's method[END_REF] is used to compute it. It usually requires between 5 and 9 evaluations of the polynomial to obtain an error smaller than

10 -4 on S (k) a /S max or S (k)
b /S max . It is the heaviest part of the inversion but it occurs generally for a small number of sections as experimented and seen through the following result for narrow enough sections. Proposition 3.2. Let f be a C 2 function on [S k-1 ,S k ] and n k and m k be its moments of order 0 and 3/2. Then µ

(k) in f ≤ m k n k ≤ µ (k) sup as soon as ∆S k f ∞ + 7 12 ∆S 2 k f ∞ ≤ f (S k-1 ), where . ∞ is the norm of L ∞ (S k-1 ,S k ).
Proof. Let us compute the following ratio:

R = 2 µ (k) sup -µ (k) in f   m k n k - µ (k) in f +µ (k) sup 2   = ∆S k S k S k-1 f (S) S 3 2 ∆S k - S k S k-1 σ 3 2 dσ dS S k S k-1 f (S)dS S k S k-1 (2σ-S k -S k-1 )σ 3 2 dσ . Then, denoting d 2 (S) = f (S)-f (S k-1 )-(S-S k-1 ) f (S k-1 ) : |R| = S k S k-1 d 2 (S) S 3 2 ∆S k - S k S k-1 σ 3 2 dσ dS+ S k S k-1 (S-S k-1 ) f (S k-1 ) S 3 2 ∆S k - S k S k-1 σ 3 2 dσ dS 1 ∆S k ∆S k f (S k-1 )+ S k S k-1 ( f (S)-f (S k-1 ))dS S k S k-1 (2σ-S k -S k-1 )σ 3 2 dσ ≤ f ∞ δ 2 + f ∞ |δ 1 | f (S k-1 )-∆S k 2 f ∞ , as soon as f (S k-1 ) > ∆S k 2
f ∞ and with

δ 1 = S k S k-1 (S-S k-1 ) S 3 2 ∆S k - S k S k-1 σ 3 2 dσ dS S k S k-1 (2σ-S k -S k-1 )σ 3 2 dσ = ∆S k 2
and

δ 2 = S k S k-1 (S-S k-1 ) 2 2 S 3 2 ∆S k - S k S k-1 σ 3 2 dσ dS S k S k-1 (2σ-S k -S k-1 )σ 3 2 dσ ≤ ∆S 3 k 6 -S 3 2 k-1 ∆S k + S k S k-1 σ 3 2 dσ S k S k-1 (2σ-S k -S k-1 )σ 3 2 dσ Introducing R k-1 = √ S k-1 and R k = √
S k and after some algebra:

δ 2 ≤ 7 12 ∆S 2 k 2R 4 k +6R 3 k R k-1 +10R 2 k R 2 k-1 +9R k R 3 k-1 +3R 4 k-1 3R 4 k +9R 3 k R k-1 +11R 2 k R 2 k-1 +9R k R 3 k-1 +3R 4 k-1 ≤ 7 12 ∆S 2 k . If ∆S k f ∞ + 7 12 ∆S 2 k f ∞ ≤ f (S k-1 ), then |R| ≤ 1 and consequently µ (k) in f ≤ m k n k ≤ µ (k) sup .
This result induces that the case [C.2] will concern an increasing fraction of the sections, when their width tends to zero. Moreover, in this case, formulas for the reconstruction are analytical, whereas in the other cases, the root of a fifth order polynomial has to be computed. This could induce a lower cost of the reconstruction itself, compared to the exponential reconstruction, which will be checked in Section 3.1.5 with some examples. Let us now analyze this method as regards stability and accuracy.

Stability of the affine reconstruction

Since numerical methods will induce errors on the moments, let us analyzed the effect of such errors on the reconstruction. But instead of the moments itself, one considers the following independent variables:

(n k ,µ k ) k∈{1,...,N} |∀k, n k ≥0, µ k ∈]S 3/2 k-1 ,S 3/2
k [} to the set of piecewise C 1 functions equipped with the L ∞ norm. Moreover, its partial derivatives are given by:

∂ µ k f a (S) = n k 2S-S (k) a -S (k) b S (k) b S (k) a σ 3 2 (2σ-S (k) a -S (k) b )dσ 1 [S (k) a ,S (k)
b ] (S),

∂ n k f a (S)=                            2 S (k) b -S S (k) b -S k-1 2 1 [S k-1 ,S (k) b ] (S) if µ k ∈]S 3 2 k-1 ,µ (k) in f [,   2 S k -S ∆S 2 k µ (k) sup -µ k µ (k) sup -µ (k) in f +2 S-S k-1 ∆S 2 k µ k -µ (k) in f µ (k) sup -µ (k) in f   1 [S k-1 ,S k ] (S) if µ k ∈ [µ (k) in f ,µ (k) sup ], 2 S-S (k) a S k -S (k) a 2 1 [S (k) a ,S k ] (S) if µ k ∈]µ (k) sup ,S 3 2 k [. Any moment of ∂ n k f a (S) or ∂ µ k f a (S)
is bounded, independently from ∆S (one can even remark that the moments of order 0, 1/2 or 3/2 of ∂ n k f a (S) are smaller than a constant times n k ). This ensures the stability of the reconstruction.

Accuracy of the affine reconstruction

Since we use the reconstruction to close our equations, we can state that the closer the reconstruction is to the exact distribution, the more accurate our model will be. Here, the accuracy of the reconstruction is analyzed assuming that the values of the moments are exact. The effect of the error on the moments was analyzed in the previous section, devoted to stability. In the asymptotic limit of a large number of sections, one has the following result. Theorem 3.1 (Accuracy of the piecewise affine reconstruction). Let f be a C 2 fonction on [0,S max ] with S max < ∞ with a finite number of zeros. Let n k and m k be its moments of order 0 and 3/2 on ]S k-1 ,S k [, for k ∈{1,2,...,N} and f a the corresponding piecewise affine reconstruction given by the proposition 3.1. Then, f a is a second order approximation of f when the maximal section width ∆S tends to zero: ff a L 1 (0,S max ) = O(∆S 2 ). 25 12 ∆S 2 f L ∞ (0,S max ) } in such a way that ∪ k ∈I 2 [S k-1 ,S k ]⊂I ∆S , thanks to the proposition 3.2. For a small enough value of ∆S, I ∆S is a union of intervals [a,b] around each zero λ of f , and

I ∆S = {S ∈ [0,S max ] | f (S)≤2∆S| f (S)|+
b a | f (x)|dx = λ a f (x)dx + b λ f (x)dx = | f (a)|+| f (b)| = O(∆S).

It means that

I ∆S | f (S)|dS = O(∆S) and ∑ k ∈I 2 S k S k-1 f (S)dS ≤ I ∆S f (S)dS ≤ I ∆S 2∆S| f (S)|+ 25 12 ∆S 2 f L ∞ (0,S max ) dS = O(∆S 2 ).
Moreover, since f and f a have the same integrals on each section:

∑ k ∈I 2 S k S k-1 | f (S)-f a (S)|dS ≤ 2 ∑ k ∈I 2 S k S k-1 f (S)dS = O(∆S 2 ).
For k ∈ I 2 , the reconstruction corresponds to the case [C.2]. Let then us define the function u, C 2 on [0,1], by: u

(x) = ( f -f a )(S k-1 +x∆S k ) for x ∈ [0,1]. One have u = ∆S 2
k f since f a is affine on the entire section. Moreover, since f and f a have the same moments of order 0 and 3/2 on ]S k-1 ,S k [, the function u is such that

1 0 u(x)dx = 0, 1 0 (S k-1 +x∆S k ) 3/2 u(x)dx = 0.
(3.9)

One introduces a Taylor expansion on it, between the values y ∈ [0,1] and x ∈ [0,1]:

u(x) = u(y)+(x-y)u (y)+ (x-y) 2 2 u (v x,y ),
with v x,y a value between x and y. This expansion is introduced in (3.9) and then

A u(y) u (y) = - 1 2 1 0 (x-y) 2 u (v x,y )dx 1 0 (x-y) 2 (S k-1 +x∆S k ) 3/2 -(S k-1 +y∆S k ) 3/2 u (v x,y )dx , with A = 1 1 2 -y α β , α β = 1 0 1 x-y (S k-1 +x∆S k ) 3/2 -(S k-1 +y∆S k ) 3/2 dx.
One can then solve the system and write u(y):

u(y)=- β 1 0 (x-y) 2 u (v x,y )dx+(y-1/2) 1 0 (x-y) 2 (S k-1 +x∆S k ) 3/2 -(S k-1 +y∆S k ) 3/2 u (v x,y )dx 2 1 0 (x-1/2)(S k-1 +x∆S k ) 3/2 dx . A majoration is then done, introducing γ = S k-1 /dS: |u(y)| ≤ u L ∞ (0,1) F(γ) with F(γ) = 1 3 1 0 (x-y)[(γ+x) 3/2 -(γ+y) 3/2 ]dx+ 1 2 1 0 (x-y) 2 (γ+x) 3/2 -(γ+y) 3/2 dx 1 0 (2x-1)(γ+x) 3/2 dx .
If γ is not zero, then it can be written:

F(γ)= 2 1 0 (x-y)[(1+x/γ) 3/2 -(1+y/γ) 3/2 ]dx+3 1 0 (x-y) 2 (1+x/γ) 3/2 -(1+y/γ) 3/2 dx 6 1 0 (2x-1)[(1+x/γ) 3/2 -1]dx .
Then, F is a continuous function of γ defined on R + and bounded since it tends to the following constant, when γ tends to infinity:

F(γ) → 2 1 0 (x-y) 2 dx+3 1 0 |x-y| 3 dx ≤ 2 3 + 3 4 = 17 12 . Then f -f a L ∞ (S k-1 ,S k ) ≤ C∆S 2 k f L ∞ (S k-1 ,S k ) , with a constant C independent from ∆S k . And finally f -f a L 1 (0,S max ) ≤ ∑ k ∈I 2 S k S k-1 | f (S)-f a (S)|dS+ ∑ k∈I 2 f -f a L ∞ (S k-1 ,S k ) ∆S k = O(∆S 2 ).
One can remark that the second order accuracy is proven in L ∞ norm in each section corresponding to case [C.2]. Moreover, the previous result is easily generalizable to the case where f is equal to zero in a finite number of sub-intervals of [0,S max ]. The aim of this section is to compare the accuracy of the two types of TSM reconstructions and of the OSM reconstruction on various probability density functions (PDFs): their numerical order of accuracy as well as their accuracy for a small number of sections. This allows to test the quality of the closure given by the models.

Let us then introduce the following PDFs f 0 , defined on [0,1], and plotted in Figure 2:

• a regular PDF: 2 with I 0 = 0.9963

f 0 (S) = (1+8S)(1-S) 2 I 0 exp 0.001 1- 1 (1-S)
• a bimodal PDF, defined as a sum of two beta PDFs:

f 0 (S)=10 2S(1-S) 4 +S 4 (1-S)
• a beta PDF:

f 0 (S) = 105 S 4 (1-S) 2
• a gamma PDF: f 0 (S) = One computes moments of order 0 and 3/2 of f 0 on each interval of a regular discretization of [0,1] in N sections. The section width is then ∆S = 1/N. From these moments, the affine reconstruction and the exponential one are computed and compared with the original function f 0 . An example of such kind of comparison is given in Figure 3, for the bimodal PDF, with 4 and 8 sections. The corresponding L ∞ norm of the difference between f 0 and its reconstruction, depending of the section width, is then plotted in Figure 4 for the two types of reconstruction. The second order of accuracy of the affine reconstruction is observed. It is a stronger result than the one shown in theorem 3.1, since L 1 norm was considered in this theoretical result, or L ∞ one but for the sections where an affine reconstruction is possible in the "complete" section (case [C.2]). It is interesting since the evaporation induces the use of pointwise values of the reconstruction. Moreover, the case [C.2] is the most usual: the cases [C.1] and [C.3] only appear for at most 3 sections for the bimodal, beta and gamma PDFs and on a decreasing proportion of the sections for the lognormal and regular PDFs, up to 0.2% or 0.5% of the total number of sections for the most discretized case. For the exponential reconstruction, there is a reduction of the order of accuracy for the bimodal and the beta PDFs. This is due to the loss of accuracy of the computation of the parameter b when |b| is large enough and to the limitation of this parameter for the computation. Indeed, the limitation occurs in a small number of sections (about 0.1%) but |b| is larger than 200 in about 1% of the sections and larger than 100 in about 2%. Moreover, one can also remark that the accuracy of the affine reconstruction is higher than the exponential one in all tested cases and is high even with a small number of sections. For example, the error is about 1% for 16 sections with the regular PDF. Concerning the computation cost, the affine reconstruction is usually more than twice as fast as the exponential one.

Velocity reconstruction

In the previous works [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[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][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF], a constant velocity was always considered in the section, except in [START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF]: χ k (S) = u k . This leads to:

U (k) 2 = u k ⊗u k , E (k-1 2 ) u = E (k-1 2 ) u k , M (k)
u = M (k) u k , F (k) = u g -u k St , 1 St = φ k φ k-1 V(φ) St(φ) κ k (φ)dφ φ k φ k-1 V(φ)κ k (φ)dφ and Q 3 * ijk = Q 2 * ijk u i , Q 1 ijk Q 2 * ijk = |u i -u j | D ijk 1 (r * ) 3 κ i (r * )κ j (r )β(r * ,r )dr * dr .
(3.10)

To increase the accuracy on the description of velocities as conditioned by size, a velocity distribution that is polynomial in size was considered in the sections in [START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF], at the cost of transporting extra size-velocity crossed-moments. Here, we test a MUSCL type reconstructions to increase the accuracy with a smaller increase of cost. The size conditioned velocity is written as a function of surface and reads:

u d (S) = N ∑ k=1 χ k (S)1 [S k-1 ,S k ] (S), χ k (S) = u k +γ (k) (S-S (k) avg ). (3.11)
The parameter S (k)

avg is chosen such that m k u k = S k S k-1 S 3/2 χ k (S)κ k (S)
dS and the slope γ (k) is computed thanks to a minmod limiter:

S (k) avg = S k S k-1 S 5/2 κ k (S)dS S k S k-1 S 3/2 κ k (S)dS ∈ ]S k-1 ,S k [, γ (k) = minmod u k+1 -u k S (k+1) avg -S (k) avg , u k -u k-1 S (k) avg -S (k-1) avg .
Let us remark that for the affine TSM, S

avg can be written as a rational function of S (k) a and S (k) b , whereas its analytical expression is more complex for the exponential TSM. In all cases, this velocity reconstruction should give us a second order reconstruction of the velocity, as a function of the droplet surface.

Conclusion on the reconstructions

At this point, the system (2.8-2.10) is closed by a reconstruction of the size distribution κ k and of the velocity χ k from the moments n k , m k and m k u k . Considering the size distribution, the affine reconstruction introduced in this paper is shown to be stable and second order accurate in all cases and not only in the case [C.2] where the reconstruction linearly depends on the moments. Moreover, its cost is twice as low as the exponential reconstruction and its accuracy is better on tested cases. Finally, its simple shape makes the source term computation easier. Concerning the velocity, an affine reconstruction is proposed, instead of the usually used constant one, in order to improve the description of this quantity inside the section.

The system (2.8-2.10) has now to be solved numerically under the constraint of preserving the realizability of the moments in each section. To do so a splitting strategy is used. Then dedicated schemes are developed separately for each operator. Since transport schemes were already developed in [START_REF] Dufour | Modélisation multi-fluide eulérienne pour les écoulements diphasiques à inclusions dispersées[END_REF][START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF], one focus on evaporation and drag and on coalescence.

Scheme for evaporation and drag

The system corresponding to evaporation and drag is the following ODE system:

d t n k =Kκ k+1 (t,S k )-Kκ k (t,S k-1 ), (4.1)
d t m k =KS 3/2 k κ k+1 (t,S k )-KS 3/2 k-1 κ k (t,S k-1 )- 3 2 S k S k-1 √ SKκ k (t,S)dS, (4.2)
d t (m k u k ) =KS 3/2 k κ k+1 (t,S k )χ k+1 (t,S k )-KS 3/2 k-1 κ k (t,S k-1 )χ k (t,S k-1 ) (4.3) - 3 2 S k S k-1 √ SKκ k (t,S)χ k (t,S)dS+ S k S k-1 √ S u g -χ k (S) A dS,
where K = K S is the constant evaporation rate and A = St(S)/S is a constant for the considered Stokes drag force. In fact, there is one ODE system for each point x (or each cell of the mesh if a spatial discretization is used). But, since these systems are independent, the spatial dependence is skipped for the sake of simplicity and without loss of generality. We will first deal with the two first equations, corresponding to the pure evaporating case, since they are independent of the velocity, here. Let us first remark that directly using on this system an explicit Euler scheme, or any explicit SSP Runge-Kutta method, may return non-realizable moments and will not be robust, as shown in Appendix B. Various realizable schemes can be used. When considering integer moments, a kinetic scheme was introduced in [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] for the constant evaporation law. Another robust scheme is developed in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] for high order moment methods for eventually non-constant evaporation rate. In the following, this last scheme is referred to as the "quadrature kinetic scheme" since we show that it can be interpreted as a quadrature of the kinetic scheme's formulae. Here, in our context of fractional moment method, the kinetic scheme is first adapted. Then, the quadrature kinetic scheme and a simplified version of it are developed and numerically tested. Finally, the quadrature kinetic scheme is extended to the additional consideration of the drag force.

Case with evaporation alone 4.1.1 Kinetic scheme (KS) for the evaporating case

The system (4.1-4.2) is related to the following kinetic equation on the droplet distribution n(t,S) = f (t,u,S)du:

∂ t n-∂ S (Kn) = 0. (4.4)
This kinetic equation has the following analytical solution: n(t,S) = n(0,S+Kt). Let us denotes n n k and m n k some approximation of the moments of order 0 and 3/2 of n(t n ,.) on the section [S k-1 ,S k ]. Then, a kinetic scheme (KS) is defined which consists in solving exactly (4.4) on a time step, from a reconstruction This can be rewritten, for k > 1, in such a way that some flux terms appear:

∑ N k=1 κ n k (S)1 [S k-1 ,S k ] (S)
n n+1 k =n n k - S k-1 +K∆t S k-1 κ n k (S)dS+ S k +K∆t S k κ n k+1 (S)dS (4.7) m n+1 k =m n k + S k S k-1 (S-K∆t) 3/2 -S 3/2 κ n k (S)dS (4.8) - S k-1 +K∆t S k-1 (S-K∆t) 3/2 κ n k (S)dS+ S k +K∆t S k (S-K∆t) 3/2 κ n k+1 (

S)dS

This scheme, applied to the use of two integer moments, was analyzed in [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF]: with the use of an affine reconstruction on the complete section (corresponding to case [C.2] here), the scheme was shown to be second order accurate. Here, thanks to Theorem 3.1, it can be shown easily, and in all cases, that the scheme is second order consistent in L 1 norm (and in L ∞ norm for sections in case [C.2]). Proposition 3.3 induces its stability. This kinetic scheme then highly accurate, due to the fact that the kinetic equation is solved exactly during the time step. However, it can be not easily generalized to non constant evaporation laws. It is the reason of the development of simplified schemes.

Quadrature kinetic scheme (QKS) for the evaporating case

The scheme presented in [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] is derived for a system of moments M i,k = S k S k-1 S i f (S)dS, with i = 0,1,...,2N m -1 on each section k. It is seen as a 3-step method:

1. A distribution f n is reconstructed from the moments (M n k,i ) i=0,...,2N m -1,k=1,...,N at t = t n in order to compute the "fluxes",

Φ n k,i = S k-1 +K∆t S k-1 S i f n (S)dS 2.
For each k = 1,..., N, the N m weights w k,α and abscissas σ k,α , corresponding to the quadrature of the moment vector (M k,i -Φ n k,i +Φ n k+1,i ) i=0,...,2N m -1 and then to the measure f n (S)1 ]S k-1 +K∆t,S k +K∆t[ (S)dS, are computed.

The moments (M n+1

k,i ) i=0,...,2N m -1 corresponding to the weights w k,α and the abscissas σ k,α -K∆t are computed:

M n+1 k,i = ∑ N m α=1 w k,α (σ k,α -K∆t) i .
In fact, this can be seen as a quadrature formula of the kinetic scheme given by:

M n+1 k,i = S k S k-1 S i f n (S+K∆t)dS = S k +K∆t S k-1 +K∆t (S-K∆t) i f n (S)dS.
Let us remark that since only integer moments are considered, the scheme is identical to the kinetic one, no erreor being introduced by the quadrature. Moreover, this interpretation of the scheme as a quadrature formula on the kinetic scheme is still valid for the case of size dependant evaporation coefficient K, but the use of the quadrature then induces an additional error compared to the kinetic scheme.

Here, the moments of order 0 and 3/2 are used. Quadrature formulas are used on the kinetic scheme (4.5-4.6). But, here, the interval ]S k-1 +K∆t,S k +K∆t[ is divided into the two following intervals ]S k-1 +K∆t,S k [ and ]S k ,S k +K∆t[ and the quadrature is done on each one. This sub-division is relevant for the global problem, including the drag force or a more complex evaporation rate. Indeed, since each section has its own velocity and could have its own evaporation rate as in [START_REF] Sibra | Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using eulerian multi-fluid methods[END_REF], it will be interesting to consider separately the droplets coming from their respective sections. Moreover, to obtain a good approximation, while still conserving the moment of order 3/2, we use two quadrature points and moments of order 0, 1/2, 1, 3/2, these moments being computed analytically. It can be seen as a classical two-point quadrature on the measure obtained by the change of variable R = √ S. An analytical formula can then be used [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF]. Let us then denote (w 1 ,w 2 ) and (

√ σ 1 , √ σ 2 ) the quadrature points of 2Rκ n k (R 2 )1 ]S k-1 +K∆t,S k [ (R 2
)dR and ( w1 , w2 ) and

( √ σ1 , √ σ2 ) the quadrature points of 2Rκ n k+1 (R 2 )1 ]S k ,S k +K∆t[ (R 2 )dR.
They are then such that:

S k S k-1 +K∆t S j/2 κ n k (S)dS = 2 ∑ i=1 w i (σ 0 i ) j/2 S k +K∆t S k S j/2 κ n k (S)dS = 2 ∑ i=1
wi ( σ0 i ) j/2 j ∈ {0,1,2,3}. (4.9) Then, the update values of the moments are:

n n+1 k m n+1 k = 2 ∑ i=1 w i 1 (σ i -K∆t) 3/2 + 2 ∑ i=1 wi 1 ( σi -K∆t) 3/2
Compared to the kinetic scheme, an additional error comes from the use of the quadrature for integration of some function g on the measure described above. Since the support widths of the measures are smaller than ∆S k /(R k +R k-1 ), this error will behave at most as O(∆S 4 k /(R k +R k-1 ) 4 ) as soon as the function g is C 4 , its fourth order derivative appearing in the bound. However, for the mass density, g(x) = (x 2 -K∆t) 3/2 and g (4) (x) is singular at x = √ K∆t. This means that some accuracy will be lost in the first section.

Simplified quadrature kinetic scheme (SQKS) for the evaporating case

For practical applications, a still simpler algorithm can be needed to reduce the cost, using a "quadrature" with only one point, i.e. a mean surface. This can be done directly on the moments of order 0 and 3/2. Then, one compute the values of the moments of order 0 and 3/2 on each of the separated intervals ]S k-1 +K∆t,S k [ and ]S k ,S k +K∆t[ and we deduce the value of a mean surfaces σ and σ:

n m = S k S k-1 +K∆t κ n k (S)dS, n m σ 3/2 = S k S k-1 +K∆t S 3/2 κ n k (S)dS, ñm = S k +K∆t S k κ n k+1 (S)dS, ñm σ3/2 = S k +K∆t S k S 3/2 κ n k+1 (S)dS.
The update values of the moments are then:

n n+1 k m n+1 k = n m 1 (σ-K∆t) 3/2 + ñm 1 ( σ-K∆t) 3/2
Here, compared to the kinetic scheme, the additional error coming from the quadrature is now dominant. Indeed, it leads to an error behaving as O(∆S 2 ) and then to a global first order consistency error.

Results on the purely evaporating case

To compare the schemes, the purely evaporating problem with a dimensionless evaporation rate K = 1 is resolved using the initial dimensionless distributions defined in section 3.1.5. More precisely, the regular and the bimodal distributions are used since they offer a different behavior in terms of accuracy of the affine and exponential reconstructions, as shown above. A uniform discretization in S on [0,1] is used and the CFL-like condition is fixed such that K∆t/∆S = 0.8. Computations are done from t = 0 to the end of the evaporation, t = 1, with the different kinds of schemes (KS, QKS and SQKS) for TSM with the affine reconstruction (aff TSM) and only KS for TSM with the exponential reconstruction (exp TSM) and for OSM, in order to assess the effect of the reconstruction. One then compares the reconstructed distribution f n r (S) at each time t n from the computed moments (n n k ,m n k ) k to the analytical solution f (t n ,S)= f 0 (S+Kt n ). Let us define the error by max n f n rf (t n ,.) L 1 (0,1) . This error is then presented in Figure 5, as a function of the section width ∆S = 1/N for a number N of sections from 2 to 1024 and for the two initial distributions. As expected, the KS is second order accurate with the piecewise affine reconstruction. In fact, one would have noticed the same order of accuracy by using the L ∞ norm instead of the L 1 norm. With the exponential one, it is also second order accurate when considering the regular initial distribution but there is a loss of accuracy for the bimodal initial distribution, resulting from the loss of accuracy of the reconstruction in this case. The KS for OSM leads to a first order of accuracy corresponding to the first order of accuracy of the reconstruction. Concerning the other schemes, one can see that the accuracy of the QKS is really close to the one of the KS. In fact, the L ∞ norm of the error on the distribution saturates with the QKS, due only to the first section: when eliminating this section in the evaluation of the error, one recovers a second order of accuracy with this norm. Finally, the SQKS is first order accurate as expected, but with an accuracy still at least ten times better than OSM with the KS. So, using TSM, even with this simplified scheme, still induces a much better accuracy for the consideration of the evaporation than OSM with the KS and twice as much sections, i.e. with the same number of size moments (but with a larger number of total variables for OSM when considering also the velocity).

Two global quantities are also considered since they are important for practical appli- cations: the total number density ∑ N k=1 n k and the total mass density ∑ N k=1 m k . The errors on these two variables, divided by their initial values, are plotted in Figure 6 for the regular initial distribution. First, one observes that KS is at least second order accurate with TSM and first order accurate with OSM. For the total mass density, it is even third order accurate. This is due to a compensation of fluxes between the sections, since their is no mass flux corresponding to disappearing droplets (i.e. no mass flux at zero size). It means that the mass density transferred to the gaseous phase is accurately described. Using TSM with the piecewise affine reconstruction, the QKS induces a loss of accuracy for the first section only visible on the total number density when more than 256 sections are used. The SQKS is second order accurate for the total mass density but only first order accurate for the total number density. However, the error induced by this scheme is smaller than 1% as soon as at least 3 sections are used.

Finally, with the QKS, the affine TSM is more than two times faster than the exponentiel TSM with a better accuracy (results of QKS with the exponentiel TSM are not represented here, for the sake of clarity).

Case with evaporation and drag

Quadrature kinetic scheme (QKS) for evaporation and drag

In the following, we use the QKS with the affine reconstruction for their accuracy, their simplicity, and their reasonable cost. The drag force is solved together with the evaporation in the following way . First, the quadrature points (w i ) i=1,2 , (σ 0 i ) i=1,2 and ( wi ) i=1,2 , ( σ0 i ) i=1,2 are defined as in section 4.1.2 by (4.9). Then, the evolution of the abscissae dur-ing the time step is given by solving the systems:

           d t σ i (t) = -K, σ i (0) = σ 0 i , d t v i (t) = u g -v i (t) Aσ i (t) , v i (0) = χ k (σ 0 i ),            d t σi (t) = -K, σi (0) = σ0 i , d t ṽi (t) = u g -ṽi (t) A σi (t) , ṽi (0) = χ k+1 ( σ0 i ).
Let us also notice that we can provide an analytical solution if the gas velocity u g is constant.

As an example for the first one, the analytical solution reads:

σ i (t) = σ 0 i -Kt v i (t) = u g +(u k -u g ) 1- Kt σ 0 i 1 KA
.

The updated moment values are then:

  n n+1 k m n+1 k m k u n+1 k+1   = 2 ∑ i=1 w i   1 (σ i (∆t)) 3/2 (σ i (∆t)) 3/2 v i (∆t)   + 2 ∑ i=1 wi   1 ( σi (∆t)) 3/2 wi ( σi (∆t)) 3/2 ṽi (∆t)  
Let us remark that this scheme still corresponds to a quadrature approximation of a kinetic scheme, using the analytical solution of the kinetic equation for evaporation and drag. Moreover, the same strategy can be used to solve the drag force together with the SQKS for the evaporation.

Result for evaporation and drag

Let us consider a test case with evaporation and drag, with 1 KA =9.89. As in section 4.1.4, we use the dimensionless evaporation rate K = 1 and the regular PDF f 0 , defined in section 4.1.4, as initial dimensionless size distributions. The mono-kinetic assumption is used for the initial condition, so that the initial NDF is:

f (0,S,u) = f 0 (S)δ(u-u 0 (S)), u 0 (S) = u g 1+2S 2 - 4 3 S 3 + S 4 4 , (4.10)
where u g = 1 is the gas velocity, assumed to be constant. The velocity profile u 0 , represented in Figure 8-right (solid line with circles), is typical of the velocity of a spray having evolved in a decelerating gas velocity field. The analytical solution is then:

f (t,S,u) = f 0 (S-Kt)δ(u-u d (t,S)), u d (t,S) = u g + u 0 (S+Kt)-u g 1+ Kt S -1 KA .
Resolution of system (4.1-4.3) is conducted using the QKS described above for OSM and for the affine TSM, with constant or affine reconstruction of the velocity. The same discretizations as in section 4.1.4 are used and only the accuracy on momentum is studied here, since the other variables were already considered in section 4.1.4.

The error on the mean velocity of the spray (∑m k u k /∑m k ) and on the momentum, divided by the initial value of these quantities, are given in Figure 7. All methods involving a constant velocity reconstruction in each section are first order accurate. However, the error of TSM is much lower than the one of OSM for the momentum but similar for the mean velocity due to error compensations in this case. The affine reconstruction of the velocity does not improve the error of OSM: it is the same for the momentum and it is even worse for the mean velocity. For the affine TSM, it increases the accuracy: the method is second order accurate, with a much smaller error for all the discretizations.

Scheme for coalescence

The model for the purely coalescing problem is written:

d t   n k m k m k u k   = - 1 Kn N ∑ i=1 N ∑ j=1    Q 1 kij Q 2 * kij Q 3 * kij   + 1 Kn N ∑ i=1 N ∑ j=1    1 2 Q 1 ijk Q 2 * ijk Q 3 * ijk   .
(

Both cases of a constant and a linear reconstruction of the velocity in each section are considered here, in such a way that the terms Q 1 ijk , Q 2 * ijk and Q 3 * ijk are given by (3.10) and

Q 3 * ijk = Q 2 *
ijk u i or by (2.11). The complete convergence analysis of this problem, when the size width tends to zero, is out of the scope of this paper. But it will be numerically tested, after giving a numerical ODE solver able to ensure the realizability.

Numerical scheme

Let us first consider an explicit Euler method for this ODE system. Proposition 5.1 (Realizability for the explicit Euler method on the coalescing problem). From a set (n n k ,m n k ) k=1,...,N of realizable moments and from (u n k ) k=1,...,N , one step of the explicit Euler scheme on ODE system (5.1) leads to realizable moments

(n n+1 k ,m n+1 k ) k=1,...,N if ∆t ≤ 1 max j,k |u n k -u n j |max k r * >0 β(r k ,r * ) f a (r * 2 )2r * dr * (5.2)
where

f a (S) = ∑ N k=1 κ n k (S)1 [S k-1 ,S k ] (S) is the reconstructed distribution from the moments (n n k ,m n k ) k∈{1,...,N} . Proof. Let us define a n k = m n k -V k-1 n n k and b n k = V k n n k -m n k .
These variables are positive. Moreover, with the change of variable V = S 3/2 and f v (V)dV = f a (S)dS:

a n+1 k = a n k + ∆t 2 V k V k-1 (v-V k-1 ) v * ∈[0,v] f v (v-v * ) f v (v * )β(v-v * ,v * )|u d (v * )-u d (v-v * )|dv * dv -∆t V k V k-1 (v-V k-1 ) f v (v) v * ≥0 f v (v * )β(v,v * )|u d (v * )-u d (v)|dv * dv.
Since the velocity reconstruction u d does not introduce values bigger of smaller than the (u n k ) k∈{1,...,N} , one can write:

a n+1 k ≥ a n k -∆t max j,k |u n k -u n j | V k V k-1 (v-V k-1 ) f v (v) v * ≥0 f v (v * )β(v,v * )dv * dv ≥ a n k 1-∆t max j,k |u n k -u n j | v * ≥0 f v (v * )β(V k ,v * )dv *
This quantity is then positive if ∆t is limited as in (5.2). In the same way, we prove that b n+1 k is positive in this case, which concludes the proof.

Then any SSP method can be applied to solve the ODE system (5.1), as soon as the time step is limited as in (5.2). We use here the third order three stage Runge-Kutta SSP method, with the same limitation (5.2) on the time step as for the explicit Euler method (see for example [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF]). Let us notice that a sufficient limitation is given by

∆t ≤ 1 max j,k |u n k -u n j |π V 2/3 k n n tot +2V 1/3 k (n n tot ) 2/3 (m n tot ) 1/3 +(n n tot ) 1/3 (m n tot ) 2/3
where n n tot = ∑ N k=1 n n k (which decreases with time) and m n tot = ∑ N k=1 m n k (which is constant). This limitation is used for our computations.

Moreover, one has also to compute the terms Q 1 ijk , Q 2 * ijk and Q 3 * ijk . To compute them analytically, an algorithm was given for OSM in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF], allowing a pre-calculation of part of them. It could be generalized for the affine TSM, when the velocity is considered constant in the sections, and the pre-calculation could also be done for the sections in case [C.2], i.e. for most of the sections. But when an affine reconstruction of the velocity is used in the section, this cannot be done so easily. Moreover, analytical formulas for the integrals are no more valid when realistic collision efficiencies are considered. So, here, these integrals are numerically approximated by the use of quadrature. The previous proof can still be done with these approximations, as soon as the chosen quadrature gives the exact value of integrals of f a (S)dS and S 3/2 f a (S)dS on the sections. One can remark that this condition is not rigorously satisfied by the quadrature used till now for the exponential TSM: neither Gauss-Legendre quadrature nor adapted quadrature (corresponding to the measure κ k (S)dS) [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. Realizability at the computational level is then still an open question in this case. For the affine TSM, a Gauss-Legendre quadrature in the R variable is used, with at least 4 quadrature points, ensuring then the realizability of the scheme.

Results on the purely coalescing case

In order to test the accuracy of the method for the description of coalescence, let us consider a fictitious homogeneous configuration with a spray constituted of droplets initially at a velocity depending on their size, in such a way that coalescence will occur. More precisely, the Knudsen number is Kn = 0.0273 and the initial distribution is given by (4.10), similarly to the test case on evaporation and drag (section 4.2.2), where f 0 is the lognormal PDF defined in section 3.1.5. Computations are done with the OSM and the affine TSM models, with a constant and an affine reconstruction of the velocity inside the sections. For OSM with a constant reconstruction of the velocity, the coalescence source terms are computed thanks to the analytical formulae given in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF], whereas for the other methods, a 5 or 10 point Gauss-Legendre quadrature is used to compute the terms Q 1 ijk , Q 2 * ijk and Q 3 * ijk . The sections are equidistributed for S ∈ [0,S max ] with S max =4 and the last section is [S max ,+∞[. This value of S max is large enough so that the mass of the last section will be small during all the computation (0.036% of the mass corresponds to droplets of surface bigger than S max at the end of the computation). The sections width can vary but the time step is constant equal to 10 -2 s, for a time varying between 0 and 5. Indeed, the objective here is not to evaluate the well known ODE solver but to evaluate the model: the error induced by this choice is less than 3.10 -6 for the reference solution on all the variables we will look at, as checked thanks to a computation with a ten times smaller time step. Moreover, the time step limitation (5.2) does not directly depend on the section width and we check during the computation that this limitation is always enforced for the chosen time step. Since there is no analytical solution, the comparison will be done with a reference computation. This reference is given by a computation using OSM with a constant velocity reconstruction in the section. Moreover, the discretization involves 5121 sections for S max = 5, i.e. with an 8 times smaller width of the sections than the most refined of the other computations. The evolution of the size distribution as well as of the velocity profile is represented in Figure 8, using this solution.

-2,5 -2 -1,5 -1 -0,5 0 logarithm of the section width The global mass and the momentum are conserved by construction of the scheme. To evaluate the accuracy of the different methods, we look first at the evolution of some other global variables: the total number density and the Sauter mean radius, defined by

+∞ 0 S 3/2 f (S)dS +∞ 0 S f (S)dS
. The L 1 norm of the error on these quantities (normalized by their initial value) as a function of the section width divided by S 0 is plotted in Figure 9. A 5 point Gauss-Legendre quadrature is used to compute the coalescence terms, here. But when using 10 point Gauss-Legendre quadratures, almost the same results are obtained. One can see that the accuracy of TSM with an affine velocity reconstruction in the section is much higher than for all the other methods, with a numerical order of accuracy almost equal to 3 and an accuracy better than the other methods with a two times smaller section width, for most of the discretizations. With this method, 17 sections are enough so that the error on the total number is less than 1%. One can also remark that the TSM method with a constant reconstruction of the velocity in the section is less accurate for the evaluation of the mean Sauter radius than OSM. This is due to the fact that, in a section (for an increasing velocity profile, like here), the velocity is underestimated for the largest droplets whereas it is overestimated for the smallest ones. If the coalescence is dominant with the biggest droplets, it means that the coalescence is underestimated for the largest droplets whereas it is overestimated for the smallest ones and this effect is reinforced by the affine representation of the distribution in the section. We also look at the repartition of mass and momentum in the section by evaluating the L 1 norm of the errors on the number, mass and momentum:

1 N t N t ∑ n=1 N ∑ k=1 |n n k -nn k | 1 N t N t ∑ n=1 N ∑ k=1 |m n k -mn k | 1 N t N t ∑ n=1 N ∑ k=1 |m n k u n k -mn k ûn k |
where nn k , mn k , mn k ūn k are the number, mass and momentum given by the reference solution and corresponding to section k of the computation. Here, the results are shown for computations using a 10 point Gauss-Legendre quadrature for the coalescence terms, except for TSM with the affine velocity reconstruction were both kinds of quadratures are used to see their effects. One still observe the better accuracy of TSM with the affine velocity reconstruction compared to the other methods. A saturation of the error is observed: it is due to the use of the quadrature for the computation of the coalescence source terms. This saturation occurs at a level that depends indeed on the number of quadrature points: the error on the computation of the coalescence source terms becomes then significant compared to the error due to the size discretization.

Finally, the results of the exponential TSM are not presented here, but it was tested with the adapted quadrature [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] using two or three nodes. This quadrature points needs to be computed for each reconstruction and for each domain D ijk . Analytical formula can be used with only two nodes but the error of this approximation becomes significant as soon as more than 30 sections are used. Then, for similar accuracy, the cost of this model is three times till twenty times larger than the affine TSM using 5 quadrature nodes.

Validation on an academic nozzle configuration featuring sizevelocity coupling

The models and the corresponding schemes have been evaluated separately for each of the evaporation, drag and coalescence processes. Let us then compare them in a configuration coupling all these phenomena. In non-uniform flows the evolution of velocities depends on sizes through the drag process while the evolution of sizes depends on the size-conditioned velocity differences through evaporation and coalescence. A decelerating nozzle is considered as a baseline non-uniform flow since it exhibits the above coupling. By assessing how efficiently this case is rendered with the various strategies, we prove the ability of our moment method to tackle the difficulties that occur in industrial cases such as Diesel sprays [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF] and particle phases in solid rocket motors [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF].

Description of the configuration

Let us consider the decelerating nozzle test case with an unsteady injection of droplets defined in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. The OSM and affine TSM methods are compared on this test case where evaporation and coalescence significantly occur. Only essential characteristics of the configuration are given here, in a dimensionless way ‡ . A spray, transported by an incompressible gas flow, is injected in a 2D axi-symmetrical diverging conical nozzle. The streamlines of the gas as well as the trajectories of the droplets are assumed to be straight lines, so that only the dependance in the axial coordinate z has to be considered. Moreover the effect of the droplets on the gaseous phase is not taken into account (see [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] for a discussion of the validity of the assumptions). Then, the gaseous velocity is characterized by its axial component u z (z) = z 0 z 2 , where z 0 = 3.1 is the coordinate of the nozzle inlet. The droplets slow down because of the deceleration of the gas flow, at a rate depending on their size and inertia, inducing velocity differences and therefore coalescence.

The NDF of the spray at the inlet is given by:

f (t,z 0 ,S,u) = f 0 (S)[1+0.9sin(2πt)]δ(u-u z (z 0 )).
The injected droplet distribution f 0 was constant as a function of radius on the interval R ∈ [1/2,1] and zero elsewhere in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. But in this case, the accuracy is limited by the ‡ Compared to the case described in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF], dimensionless variables are defined thanks to the use of typical values for the time (T ≡ 0.01611s, the period of the oscillations at the injection), for the droplet surface (S 0 = 4πR 2 0 with R 0 = 30µm), for the velocity (u z (z 0 ) = 4ms -1 , the velocity at the nozzle inlet) and for the droplet number density (n 0 = 2.9 10 10 m -3 ). The Knudsen number is then defined by Kn = 1/(T u z (z 0 )S 0 n 0 ). singularity of the distribution. To avoid that, a quite similar but regular distribution is used: the beta PDF defined in section 3.1.5. The Knudsen number is Kn = 0.0471, the evaporation rate is K = 0.283 and the Stokes number is characterized by A = St(S) S = 0.357. The spatial domain is [z 0 ,2z 0 ] so that only a few droplets go out of this domain, most of them being completely evaporated before they could go through the domain. The computations are done for a time corresponding to 9 periods of the injected mass oscillations, since after 8 periods, the solution becomes periodic in the whole domain.

Numerical strategy and results

Concerning the numerical strategy used for the computations with OSM or TSM, a Strang splitting is used between the transport in the physical space and the transport in the phase space (source terms of evaporation, drag and coalescence). All numerical schemes have to preserve the non-negativity of the mass densities of the sections and also the realizability conditions for TSM. For the transport in the physical space, a second order kinetic scheme is used, based on the one given in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] but specifically adapted to the nozzle configuration. This scheme is described in Appendix C. Let us highlight in the particular case of the model considering an affine reconstruction of the velocity in the section, that the velocity is assumed to be constant in the sections for this transport operator, since it is not the aim of this paper to develop a transport scheme for this model in the nozzle configuration. Moreover, a scheme adapted to this model can be found in [START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF] for more classical configurations, with cartesian meshes. For the transport in phase space, a splitting is done between drag/evaporation and coalescence. The scheme for evaporation and drag is given in section 4.2.1 and the scheme for coalescence is given in section 5.

A 600-point uniform discretization of the spatial domain is used. The splitting time step is equal to 0.0025 in such a way that the convective CFL is about 0.48. This choice leads to small errors due to the time and space discretizations compared to the error due to the size discretization. A reference solution is obtained by using OSM with 1024 sections equidistributed in surface between 0 and S max = 7.111, plus the last section. The choice of S max allows to obtain a completely negligible mass in the last section, this section being badly represented by OSM. For the other computations, sections with a constant width ∆S = 0.25 are considered in [0,S max ], plus a last section [S max ,+∞[, with a maximal value of the surface then depending on the number N of sections: S max = (N -1)∆S. Considering OSM, a total of 20 sections are used, whereas only 12 sections are used for TSM, with a constant or an affine reconstruction of the velocity. Indeed, adding more sections only leads to a slight accuracy improvement for this method, which is able to describe correctly the last section. The profiles at t = 8.44 (in the ninth period) of the total mass density, the mass averaged velocity and the Sauter mean radius are plotted in Figure 11. One can see the accuracy of TSM compared to OSM even with a smaller number of sections. Moreover, with the affine reconstruction of the velocity in the section, TSM leads to even better results, very close to the reference solution.

Conclusion

In this paper, we have provided a very robust affine closure for the Two Size Moment method for polydisperse spray simulations, which allows to reduce the overall cost and increase accuracy and stability compared to previous approaches. In the pursuance of high order coupled methods, the velocity reconstruction has also been improved. The new method allows an easy computation of the source terms. Moreover, robust and accurate schemes were developed for simplified source terms, ensuring the realizability of the moments and easily extendable to more complex physical models, thanks to the flexibility of the present strategy. The verification of the method has been done on various cases, with evaporation, drag and coalescence, and especially on an unsteady 2D case. The high order of the reconstruction and of the schemes have then been shown in the asymptotics of a large number of sections. Also the accuracy of the method was exhibited when a reduced number of sections is used. So the description of the spray, of its size polydispersity, of its size-conditioned velocities, and of their evolution through drag, evaporation, coalescence etc., can be suitably achieved with a Two Size Moment method based on the present affine reconstruction with a low number of sections (below 10). It is important for industrial applications such as spray combustion, spray in solid rocket motors and many other industrial needs, since the numerical cost is a crucial issue. In our point of view, the present strategy is a good trade-off between complexity/accuracy and flexibility/cost without renouncing to robustness, which is a key feature of numerical simulation. In addition, this method is a good candidate for being coupled to the Anisotropic Gaussian velocity closure introduced to treat particle trajectory crossings [START_REF] Vié | On the Anisotropic Gaussian closure for the prediction of inertial-particle laden flows[END_REF]. This combination is expected to describe both the size polydispersion and also the velocity dispersion induced by some droplet trajectory crossing especially occurring in turbulent flows [START_REF] Dupif | Modélisation et simulation de l'écoulement diphasique dans les moteurs à propergol solide par une approche eulérienne polydispersée en taille et en vitesse[END_REF][START_REF] Dupif | Two-size moment eulerian multi-fluid method describing the statistical trajectory crossing: modeling and numerical scheme[END_REF].

A Affine reconstruction in volume variable

When considering the volume variables, mass and number densities in the section are then integer moments of the NDF:

n k = V k V k-1 n V (V)dV and m k = V k V k-1
Vn V (V)dV. The same kind of reconstruction as in section 3.1.2 can then be considered:

f v a (V) =        0 si V k-1 ≤ V ≤ V (k) a α k +(β k -α k ) V-V (k) a V (k) b -V (k) a si V (k) a ≤ V ≤ V (k) b 0 si V (k) b ≤ V ≤ V k (A.1)
with the same three cases as for the affine in S reconstruction, meaning that f v a is continuous inside the section. One can then easily prove the following proposition: Proposition A.1 (Affine reconstruction from integer moments). Let n k >0 and m k ∈]V k-1 n k ,V k n k [, for some k∈{1,2,..., N}. Then, their exists on ]V k-1 ,V k [ a unique non-negative continuous function f v a given by (3.4) and such that

V k V k-1 f v a (V)dV = n k , V k V k-1 V f v a (V)dV = m k . (A.2)
Its parameters are given by:

• If V k-1 < m k n k < V k +2V k-1 3 , then V (k) a = V k-1 , β k = 0 and V (k) b = 3 m k n k -2V k-1 , α k = 2n 2 k 3(m k -V k-1 n k ) . (A.3) • If V k +2V k-1 3 ≤ m k n k ≤ 2V k +V k-1 3 , then V (k) a = V k-1 , V (k) b = V k and α k = 2 (∆V k ) 2 [(2V k +V k-1 )n k -3m k ], β k = 2 (∆V k ) 2 [3m k -(V k +2V k-1 )n k ].
(A.4)

• If 2V k +V k-1 3 < m k n k < V k , then V (k) b = V k , α k = 0 and V (k) a = 3 m k n k -2V k , β k = 2n 2 k 3(V k n k -m k ) . (A.5)
The same kind of analysis as in sections 3.1.3 and 3.1.4 could be done for this reconstruction. For example, one can show that, in case 2:

f -f v a L ∞ (V k-1 ,V k ) ≤ 17 12
∆V 2 f L ∞ (0,1) .

But the evaporation with a constant rate K = K S generates NDFs f S that are nonnecessarily zero at S =0, since the solution is shift of the initial NDF. Let us then consider such kind of function (in fact, the gamma and regular functions defined in section 3.1.5) and let us do its piecewise affine reconstruction f v a from its moments. The corresponding function f S (S)= 3 2

√

S f v a S 3/2 is plotted in Figure 12 for 4 or 16 sections. When the initial distribution if zero at S = 0, as the gamma one, the reconstruction is good. But when it is not the case, the accuracy is affected, especially close to zero and the evaporation fluxes will not be accurately evaluated. This shows that this kind of reconstruction is not adapted to an accurate resolution of the evaporation.

B Backward Euler scheme for evaporation and realizability

Let us evidence the realizability problems with the use of an explicit Euler scheme directly on the purely evaporating system (4.1,4. , which tends to zero. So, in this case of an emptying section, the time step has eventually to become very small to ensure the realizability condition (2.6). This shows that this ODE solver, as well as any explicit SSP Runge-Kutta method, may introduce some realizability failures and will then not be robust.

C Transport scheme for the nozzle test case

In order to discretize the transport part of the multi-fluid system (2.8-2.10), that is, the left-hand side of the system of conservation equations for unsteady configurations, we consider a second-order finite volume kinetic scheme similar to the one introduced by Bouchut et al. for pressureless gas dynamics equations [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] and already adapted for OSM in the nozzle configuration in [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. For the sake of simplicity, the variables of each section k are denoted ρ for the mass density m k , n for the number density n k and u for the axial velocity u k . Moreover, the momentum is denoted q = ρu. In the axisymetrical configuration, there are no azimuthal velocity and thanks to the similarity assumption the axial velocity only depends on the axial coordinates z and the reduced radial velocity u r /r is u/z. The equations then reduce to the following 1D unsteady system:

∂ t   n ρ q   +2 u z   n ρ q   +∂ z (u   n ρ q   ) = 0.
As in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF], a correspondance is done with a kinetic equation, adapted here to our configuration. Indeed, starting at time t n with the functions ρ n (z), n n (z)=ρ n (z)µ n (z) and u n (z), the functions at time t n+1 can be obtained by solving the transport equation derived from the similarity assumption:

∂ t f +2 ξ z f +ξ∂ z f = 0, (C.
f (z,ξ,ζ;t) = z-ξ(t-t n ) z 2 f n (z-ξ(t-t n ),ξ,ζ) (C.2)
The new moments are then obtained by projecting the exact solution:

  ρ n+1 (z) n n+1 (z) q n+1 (z)   = ρ n+1 (z)   1 µ n+1 (z) u n+1 (z)   = R 2   1 ζ ξ   f (z,ξ,ζ;t n+1 )dξdζ
Let us notice here that the variable ζ is just added here to be able to write the equation on the number n. It represents the ratio n/ρ, ie n k /m k and it then belongs to the interval

1 S 3/2 k , 1 S 3/2 k-1 .
Here, a finite volume is used with a uniform spatial discretization, the size of the cells being denoted ∆z and the averaged values of ρ n (z), n n (z) and q n (z) on the cell j defined by [z j-1 2 ,z j+ 1 2 ] and denoted ρ n j , n n j = ρ n j µ n j and q n j = ρ n j u n j . The scheme is then: The obtained fluxes rely, through (C.2), on the knowledge of functions ρ n (z), µ n (z) and u n (z), which have to be determined from the discrete values ρ n j , n n j and q n j . We use piecewise constant over half-cells, as for the simplified second-order kinetic scheme defined in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF]. If we do not write, on purpose for the sake of legibility, the n exponent, the moments take the form

ρ(z) = ρ R j-1 2 µ(z) = µ R j-1 2 u(z) = u R j-1 2 z j-1 2 < z < z j , ρ(z) = ρ L j+ 1 2 µ(z) = µ L j+ 1 2 u(z) = u L j+ 1 2 z j < z < z j+ 1 2 ,      log z j z j-1 2 +log   1+ β R j-1 2 1+α R j-1 2   - α R j-1 2 2(1+α R j-1 2 ) + β R j-1 2 2(1+ β R j-1 2 )   +ρ L j+ 1 2 u L j+ 1 2    1 µ L j+ 1 2 u L j+ 1 2      log z j+ 1 2 z j -log 1+ β L j+ 1 2 +α L j+ 1 2   1- α L j+ 1 2 6   - β L j+ 1 2 2(1+ β L j+ 1 2 )   ,

Figure 1 :

 1 Figure 1: The three types of affine reconstructions.

Proposition 3 . 3 .

 33 n k and µ k = m k n k for each k. This change allows to separate the effect of a change of level of the function from the effect of a change of the mean volume. Computation in each of the cases [C.1], [C.2] and [C.3] allows to show the following property. The piecewise affine reconstruction defined by proposition 3.1 is a C 1 application from the space {

  Proof. In order to treat separately the case [C.2], let us denote I 2 the set of all k∈{1,2,..., N} such that the reconstruction corresponds to the case [C.2]. One also denotes

3. 1 . 5 Figure 2 :

 152 Figure 2: Model size distributions.

Figure 3 :

 3 Figure 3: Bimodal distribution (black line) and its exponential (left) and affine (right) reconstructions with 4 (red dashed line) and 8 (blue solid line) sections.

Figure 4 :

 4 Figure 4: L ∞ norm of the error for the affine (left) and the exponential (right) reconstructions.

  of the NDF corresponding to the moments n n k and m n k and in taking its moments. With the CFL-like condition K∆t ≤ ∆S k and with the convention κ n N+1 = 0, we get:

Figure 5 :

 5 Figure 5: Purely evaporating case with the regular initial distribution (left) and the bimodal one (right): maximal value in time of the L 1 norm of the error on the NDF depending on the section width. Lines of slope 1 and 2 are represented by black dashed lines.

Figure 6 :

 6 Figure 6: Purely evaporating case with the regular initial distribution: maximal value of the error, depending on the section width, on the total number (left) and on total mass (right) divided by their initial values. Lines of slope 1, 2 (left and right) and 3 (right) are represented by black dashed lines.

Figure 7 :

 7 Figure 7: Case with evaporation and drag for the regular initial distribution: maximal value of the error, depending on the section width, on the mean velocity (left) and on the momentum (right), divided by their initial values. Lines of slope 1 and 2 are represented by black dashed lines.

Figure 8 :

 8 Figure8: Size distribution (left) and velocity profile (right) for the pure coalescence test case with the lognormal initial distribution at t=0 (solid line with circles), at t = 0.5 (dashed line) and at t =1, 2, 3, 4, 5 (solid lines).

Figure 9 :

 9 Figure 9: Coalescence test case with the lognormal initial distribution: L 1 norm of the error, depending on he section width, on the total number (left) and Sauter mean diameter (right) divided by their initial value. The straight lines of slope 2 and 3 are represented by black dashed lines.

  in u TSM -affine in u -5 quad points

Figure 10 :

 10 Figure 10: Coalescence test case with the lognormal initial distribution: L 1 norm of the error, depending on the section width, on the mass (left) and on the momentum (right). The straight line of slope two is represented in black dashed line.

Figure 11 :

 11 Figure 11: Results for the nozzle test case with OSM and affine TSM with a constant or an affine reconstruction of the velocity in the section: total mass density (top left), mean velocity minus gas velocity (top right), Sauter mean radius (bottom).

Figure 12 :

 12 Figure 12: Distributions as a function of S corresponding to a linear reconstruction in V with 4 (red curve) and 16 (blue curve) sections of the gamma function (left) and of the regular one (right).

1 S 1 / 2 1 S 1 / 2 1 S 1 / 2 S k- 1 S 3 / 2 k 1 S 1 /

 11211211213211 k +∆t[Kκ k+1 (t,S k )-Kκ k k-1 )], k +∆t KS 3/2 k κ k+1 (t,S k )-KS 3/2 k-1 κ k (t,S k-1 )-Kκ k (t,S)dS . [κ k (S)-κ k+1 (S k )]dS ≤ [κ k (S k-1 )-κ k (S)]dS ≤ S k -S 3/2 κ k (S)dSLet us consider, for example κ k+1 (S k ) = 0, the first constraint is written 2 κ k (S)dS and, if m n k /n n k tends to zero (which corresponds to case 1 of the reconstruction, where S (k) b tends to S k-1 ), then the right-hand side term is equivalent to

1 )f

 1 (z,ξ,ζ;t n ) = f n (z,ξ,ζ) = M(ρ n (z),µ n (z),u n (z),ξ,ζ),where the Maxwellian M is defined for any ρ > 0, (µ,u,ξ,ζ) ∈ R4 by M(ρ,µ,u,ξ,ζ) = ρδ((ξ,ζ)-(u,µ)) and by projecting the obtained distributions in order to get their new moments. The analytical solution of the problem at the kinetic level is:

  z,ξ,ζ;t)dξdζdtdz

† For the result fo section 3.1.5, the accuracy is obtained thanks to a smaller tolerance(10 -8 ) on the evaluation of the b k , then demanding 9 to 25 evaluations of g k in most cases.

Acknowledgements

The present research was done thanks to a Ph.D. Grant from DGA, Ministry of Defense (M.S. Amiet, Technical Monitor) and a Ph. D. Grant from ONERA. The authors would also like to thank Prof. Marc Massot for several helpful discussions and for the precious feedback he provided about this manuscript.

the values in each half cells being computed in the same way as in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF]. Let us denote

Under the CFL-like condition ∆t ∆z sup z |u n (z)|≤ 1 2 , and assuming that all velocities are positive, the fluxes can be written: