F Laurent

T T Nguyen

Realizable second-order finite-volume schemes for the advection of moment sets of the particle size distribution

Keywords: population balance equation, moment method, advection, realizable scheme, finite volume, kinetic scheme

The accurate description and robust simulation at relatively low cost of a size polydisperse population of fine particles in a carrier fluid is still a major challenge for many applications. For this purpose, moment methods, derived from a population balance equation, represent a very interesting strategy. However, one of the major issues of such methods is the realizability: the numerical schemes have to ensure that the moment sets stay realizable, i.e. that an underlying distribution exists. This issue is all the more crucial that some moment vectors can be at the boundary of the moment space for practical applications, corresponding to a population of particles with only one or a few sizes. It is then investigated here for the advection operator, for which it is particularly significant. Then second order realizable kinetic finite volume schemes are designed, with two strategies for the fluxes evaluation based on the work of Kah et al.

, which are here completely revisited, extended to take into account the boundary of the moment space and any number of moments, analyzed and compared in a Cartesian mesh context. For a potential easiest generalization to unstructured meshes, simplified but still realizable versions of these schemes are also developed. The high accuracy of all the schemes is then numerically checked on 1D and 2D test cases, with Cartesian meshes, and their robustness is shown, even when some moment vectors are at the boundary of the moment space.

Introduction

Populations of non-inertial particles in a carrier fluid are encountered in several kinds of applications (see [START_REF] Ramkrishna | Population balance modeling: Current status and future prospects[END_REF] and references therein): soot in combustion applications, nanoparticles synthesis, microbubbles in biology processes, aerosol technology, ... Its evolution can be described by a population balance equation (PBE) [START_REF] Ramkrishna | Population balance modeling: Current status and future prospects[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF][START_REF] Rigopoulos | Population balance modelling of polydispersed particles in reactive flows[END_REF][START_REF] Sporleder | On the population balance equation[END_REF], which is a transport equation for the number density function (NDF) of the particles. This NDF depends on time, spatial location and one or several internal coordinates, which can for example describe the size of the particles. The PBE includes usually the spatial transport terms, describing for example advection and diffusion, and some localized source terms describing, at each spatial location, phenomena such as nucleation, aggregation, coagulation, breakup, growth or oxidation/dissolution. It is coupled with the equations, usually Navier-Stokes equations, describing the carrier fluid [START_REF] Zucca | Implementation of the population balance equation in CFD codes for modelling soot formation in turbulent flames[END_REF].

In this work, only one internal variable is considered, describing the size of the particles, assuming for example that they are spherical. In order to be able to describe the size polydispersion of the particles at a reasonable cost, the use of moment methods seems to be an interesting strategy (see for example [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF][START_REF] Barrett | A comparison of some approximate methods for solving the aerosol general dynamic equation[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF]): only a finite set of moments of the NDF are then transported. It can also be hybridized with a discretization along the internal coordinate [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF][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][START_REF] Nguyen | Solution of population balance equations in applications with fine particles: Mathematical modeling and numerical schemes[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Attarakih | Solution of the population balance equation using the sectional quadrature method of moments (sqmom)[END_REF]. However, two major issues arise for moment methods. The first one is the closure of the moment equations essentially due to the source terms in the PBE. Several strategies were used: some of them provide a functional dependence of the unknown moments using the transported moment set, such as the interpolative closure (MOMIC) [START_REF] Frenklach | Aerosol dynamics modeling using the method of moments[END_REF]. For the other ones, a NDF, or its corresponding measure, is reconstructed from the moment set, allowing evaluation of all the unclosed terms. This reconstruction can be for example the entropy maximization [START_REF] Tagliani | Hausdorff moment problem and maximum entropy: a unified approach[END_REF][START_REF] Mead | Maximum entropy in the problem of moments[END_REF][START_REF] Vié | Size-velocity correlations in hybrid high order moment/multifluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF], a sum of Dirac delta function (quadrature method of moment, QMOM) [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF] or a superposition of kernel density functions (kernel density element method, KDEM [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF] or extended quadrature method of moment EQMOM [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gasparticle flows in a LES framework[END_REF][START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF][START_REF] Madadi-Kandjani | An extended quadrature-based moment method with log-normal kernel density functions[END_REF]). The second major issue of moment methods is the realizability. Indeed, since the set of variables are the moments of a non-negative NDF (or, more rigorously, a positive measure) on R + or a sub-interval of R + , it belongs to a space strictly included in R N + , where N is the number of moments [START_REF] Shohat | The Problem of Moments[END_REF][START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF][START_REF] Gautschi | Orthogonal polynomials: applications and computation[END_REF]. This space is called the moment space. The numerical methods have to ensure that the variables stay in this moment space, i.e. that the moments stay realizable. This issue is not always considered, thus leading to unphysical results (e.g. invalid moment sets). Indeed, the classical schemes for high-order transport in physical space can lead to invalid moment sets [START_REF] Wright | Numerical advection of moments of the particle size distribution in Eulerian models[END_REF][START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF][START_REF] Vikas | Realizable high-order finite-volume schemes for quadraturebased moment methods applied to diffusion population balance equations[END_REF], as well as for the source terms [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Nguyen | Solution of population balance equations in applications with fine particles: Mathematical modeling and numerical schemes[END_REF], even if the closure itself ensures the realizability at the continuous level. This happens all the more easily when some moment sets are at the boundary of the moment space, thus corresponding to a sum of a few weighted Dirac delta functions, as obtained through nucleation. To circumvent this issue, some authors resort to moment correction algorithms [START_REF] Mcgraw | Correcting moment sequences for errors associated with advective transport[END_REF][START_REF] Wright | Numerical advection of moments of the particle size distribution in Eulerian models[END_REF] based on a necessary but eventually not sufficient condition for realizability in order to obtain a valid moment set. The cost of the method then increases and the correction spoils the overall accuracy.

It is then very important to develop realizable schemes, i.e. schemes directly preserving the realizability of the moment set. Moreover, an operator splitting strategy, solving separately the spatial transport of both phases and the source terms was shown to be efficient and well adapted to industrial-oriented codes [START_REF] Doisneau | Numerical strategy for unsteady two-way coupling in polydisperse sprays: application to solid rocket motor instabilities[END_REF]. This allows us to deal separately with the spatial transport and the source terms. Concerning the source terms, realizable schemes were already developed for moment methods where the closure is based on a reconstruction of the NDF [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF][START_REF] Nguyen | Solution of population balance equations in applications with fine particles: Mathematical modeling and numerical schemes[END_REF]. A realizable scheme was also provided for the diffusion operator in the case of QMOM [START_REF] Vikas | Realizable high-order finite-volume schemes for quadraturebased moment methods applied to diffusion population balance equations[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF]. In this work, only advection then is considered. In practice, this operator, at least when using first order explicit finite volume methods, is usually much less costly than the potentially complex source term operator, composed of one ODEs system for each considered moment vector. So, it will be very interesting to minimize the number of degree of freedom by using a high order scheme for advection, as soon as its cost is not prohibitive.

A Lagrangian type of scheme has been developed [START_REF] Attili | Application of a robust and efficient lagrangian particle scheme to soot transport in turbulent flames[END_REF]. The advection of the moments is then obtained through the advection of some numerical particles for which a moment vector is affected. The resulting scheme is then naturally realizable. However, it suffers from the same drawbacks as usual Lagrangian methods: the need of interpolation of the carrier phase properties, the non easy coupling with this phase and the complexity in term of parallelization for high performance computing. Moreover, it could need a large number of numerical particles to converge, for which the eventually costly source terms operator has to be solved. That is why only Eulerian schemes are studied here. On the one hand, a second-order realizable kinetic finite volume method has been developed in a structured mesh context [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF], when the support of the NDF is compact. It was recently applied in a context of a mesh refinement [START_REF] Essadki | Adaptive mesh refinement for polydisperse spray simulation[END_REF]. It is based on a kinetic evaluation of the fluxes thanks to the use of the analytical solution at the kinetic level and on a MUSCL type of reconstruction on the canonical moments, which define a one to one relation between the interior of the moment space and the interior of an hypercube. However, it was only applied for inertial particles and for a four moments method, the algebra being otherwise difficult. On the other hand, a pseudosecond-order realizable finite volume method [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF] has been developed in structured and unstructured mesh contexts. Fluxes computation is then based on a reconstruction of the moments at the cell interfaces; it is obtained thank to the Gauss quadrature of the moments, just reconstructing the weights. However, it was reduced to an even number of moments and suffers from some accuracy reduction when the quadrature points evolve strongly. In this work, the last two schemes are completely revisited, generalizing them to any number of moments and allowing them to deal with the boundary of the moment space without loosing the realizability, which is an hard task. They are also analyzed, especially looking for conditions to obtain the second order of accuracy, and they are compared. This is done for NDF of support included in [0, +∞), the case of a compact support being discussed in the appendix. Kinetic schemes are thus derived in a structured mesh context, first with a reconstruction on variables defining a one to one relation between the interior of the moment space and (0, +∞) N , where N is the number of moments in the set. Algorithms are then adapted to the case of the support included in [0, +∞) and generalized to any number of moments. The weight reconstruction is also considered, for an even or an odd number of moments, in a different way compared to [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF], thus not being dependent of abscissas differences between the cells. Simplified schemes are then derived with the two kinds of reconstructions, the first one being modified for this case. It will be generalizable to unstructured meshes in a cell-centered context.

The paper is then organized as follows. In Section 2, the moment equations for the pure advection case are given, as well as the realizability constraints. Then, in Section 3, realizable finite volume kinetic schemes are provided and their orders of accuracy discussed. Some simplified version of these schemes are also given in Section 4, as well as the new constraint on the CFL to guarantee the realizability. Finally some verifications are given, considering systems with high numbers of moments, first for 1D configurations with steady or unsteady and compressible carrier phase velocity fields in Section 5 and for the 2D configuration of the Taylor-Green vortices in Section 6.

Moment transport equation and realizability

In this section, moment equations for the pure advection case are first recalled. Then, the space in which the moment vector lives is described, as well as the realizability conditions. This can be done directly, using the Hankel determinants. But some interesting tools are also introduced, defining a bijection between the interior of the moment space M N and (0, +∞) N .

Moment equations

Let us consider the NDF, denoted f (t, x, ξ), of some cloud of small particles transported by a carrier fluid. The parameter ξ, which is the size of the particles, lives in the interval [0, +∞). The case of a compact support is discussed in Appendix A. In the case of pure advection, the population balance equation (PBE) then reduces to:

∂ t f + ∂ x (u f ) = 0. ( 1
)
Since the considered particles are non-inertial, u is the velocity of the carrier phase, which is a priori compressible. It is assumed to be a regular function of (t, x) in what follows.

Instead of resolving directly this PBE, one considers a finite set of moments (m k ) k∈{0,1,...,N } of the NDF, the k th order moment m k being defined by:

m k (t, x) = +∞ 0 ξ k f (t, x, ξ)dξ.
(

) 2
These moments are the solution of the following system of equations, denoting m N = (m 0 , . . . , m N ) t :

∂ t m N + ∂ x (u m N ) = 0. ( 3
)
Let us remark that this system is closed, contrary to the often encountered systems in moment methods. However, it is usually only a part of a more complex problem for which some closure is needed and can be provided through the reconstruction of a NDF from the moments (see for example [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF]). Moreover, the moment equations seem here independent from one another. In fact, they are coupled by the fact that the vector m N = (m 0 , m 1 , . . . , m N ) t has to be a moment vector of a positive measure, i.e. has to stay in the moment space. This is the realizability condition, which is detailed in the next section.

Moment space: definition and first characterization

Let us denote P the space of finite positive Borel measures on (0, +∞). And for µ ∈ P, let us denote m N (µ) the vector of moments of µ of order 0 to N , assuming that they are finite:

m N (µ) = (m 0 (µ), . . . , m N (µ)) t , m k (µ) = +∞ 0 x k dµ. ( 4
)
The moment vectors m N = (m 0 , m 1 , . . . , m N ) t lives in the N th-moment space.

Definition 2.1. The N th-moment space M N on the interval (0, ∞) is given by

M N = {m N (µ) | µ ∈ P}.
If a moment vector m N belongs to this space, it is said to be realizable and one then defines

P(m N ) = µ ∈ P | m N = +∞ 0 (1, x, . . . , x N ) t dµ .
If m N belongs to the interior of this space, it is said to be strictly realizable.

This N th-moment space is convex. To characterize it, one can introduce the Hankel determinants, defined by:

H 2n+d = m d . . . m n+d . . . . . . m n+d . . . m 2n+d , (5)
with d = 0, 1; n ≥ 0. Indeed, one has the following theorem, for which a proof can be found in [START_REF] Shohat | The Problem of Moments[END_REF][START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF][START_REF] Gautschi | Orthogonal polynomials: applications and computation[END_REF]:

Theorem 2.1. The vector m N = (m 0 , m 1 , . . . , m N ) t is strictly realizable if and only if

H k > 0, k ∈ {0, 1, . . . , N } (6)
and if it belongs to the boundary of the moment space, then there exists k ≤ N such that

H 0 > 0, . . . , H k-1 > 0, H k = 0, . . . , H N = 0. ( 7
)
In the latter case, k is denoted N (m N ) and P(m N ) is a singleton: a sum of k+1 2 weighted Dirac delta functions 1 .

Let us remark that if N (m N ) is odd, then one of the Dirac delta functions is centered at 0. Moreover, [START_REF] Zucca | Implementation of the population balance equation in CFD codes for modelling soot formation in turbulent flames[END_REF] is not a sufficient condition for a moment to be on the boundary (and not outside) of the moment space: an additional condition is then needed [START_REF] Curto | Recursiveness, positivity, and truncated moment problems[END_REF][START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF].

The moment space has a rather complex geometry, as explained in the next section. Moreover, the Hankel determinants provide algebraic relations to determine if a vector belongs to the moment space but this tool is not easy to use, since we do not want to compute all these determinants. But other quantities can be derived, linked with a one-to-one mapping of the interior of the moment space, thanks to the theory of orthogonal polynomials.

Orthogonal polynomials theory and new characterization of the moment space

Let P be the space of real polynomials. For a positive finite Borel measure µ such that its moments are well defined and for p, q in P, let us define the scalar product:

p, q = R p(x)q(x)dµ. (8)
When considering a measure µ with infinite support, a sequence (π k ) k≥0 of orthogonal polynomials relative to this scalar product ( π k , π p = δ k,p π k , π k ), where π k is of exact degree k, satisfies the following three term recurrence relation, with π 0 = 1, π -1 = 0 and β k > 0 [START_REF] Favard | Sur les polynomes de Tchebicheff[END_REF]:

π k+1 (x) = (x -α k )π k (x) -β k π k-1 (x). (9)
Conversely, if the sequence of polynomials satisfies (9) with β k > 0 for all k ∈ N, then there exists a measure µ on the real line for which the polynomials are orthogonal. Moreover, this measure is supported on [0, +∞) if and only if there exists a sequence (ζ k ) k≥1 of positive numbers such that the coefficients in the recurrence relation [START_REF] Barrett | A comparison of some approximate methods for solving the aerosol general dynamic equation[END_REF] satisfy α 0 = ζ 1 and for all k ≥ 1 [START_REF] Chihara | An introduction to orthogonal polynomials, Gordon and Breach Science Publishers[END_REF]:

β k = ζ 2k-1 ζ 2k , α k = ζ 2k + ζ 2k+1 . (10)
And this measure is supported on [0, 1] if and only if the coefficients ζ k form a chain sequence, i.e. they can be decomposed as

ζ k = p k (1 -p k-1
), with p 0 = 0 and p k ∈ (0, 1) for k ≥ 1 [START_REF] Wall | Analytic Theory of Continued Fractions[END_REF].

When considering a measure with finite support, the sequence of such orthogonal polynomials is finite and only a finite number of the ζ k (or of the p k ) can be defined. These coefficients then allow to characterize the interior of the moment space in the case where the support is included in [0, +∞). Indeed, for a moment vector m N at the boundary of the moment space M N , corresponding to a measure with a finite support, one has ζ N (m N ) = 0 and the ζ k for k > N (m N ) are not defined.

A more geometrical point of view can also be considered [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]. For a vector m k-1 in the interior of the moment space M k-1 , let us then first define

m - k (m k-1 ) = min µ∈P(m k-1 ) m k (µ), m + k (m k-1 ) = max µ∈P(m k-1 ) m k (µ), (11)
as the lower and upper boundary of the admissible interval for the moment m k of order k, the lower order moments m k-1 being known. In what follows, only measures supported on [0, +∞) are considered. In this case,

m + k = +∞ and m - k is finite, with m k -m - k = m 0 k i=1 ζ i . Moreover, m - k
strongly depends on m k-1 , meaning that the interval where the k th order moment lives strongly depends on the value of the lower order moments. Some examples are given in Appendix B. Then, the strict realizability is characterized either by the positivity of the Hankel determinants or by the positivity of the ζ k and induces a link between the moments. Several algorithms allow to compute efficiently the recursion coefficient α k and β k and then the ζ k from the moments: Rutishauser's QD algorithm [START_REF] Rutishauser | Der Quotienten-Differenzen-Algorithmus[END_REF][START_REF] Henrici | The quotient-difference algorithm[END_REF], Gordon's PD algorithm [START_REF] Gordon | Error bounds in equilibrium statistical mechanics[END_REF][START_REF] Gordon | Error bounds in spectroscopy and nonequilibrium statistical mechanics[END_REF] and variation of an algorithm attributed to Chebyshev and given by Wheeler in [START_REF] Wheeler | Modified moments and Gaussian quadrature[END_REF]. Since this last one is said to be slightly more stable in practice [START_REF] Wheeler | Modified moments and Gaussian quadrature[END_REF], let us detail it, as well as the reverse algorithm, giving the moments from the ζ k and m 0 .

Chebyshev algorithm for the computation of the ζ k and reverse algorithm

Consider the matrix Z with elements Z k,p =< π k x p >, which must be zero if k > p and which satisfy Z -1,p = 0, Z 0,p = m p and, thanks to the orthogonal polynomials recursion formula:

Z k+1,p = Z k,p+1 -α k Z k,p -β k Z k-1,p . (12)
Coefficients α k , β k are determined by:

β 0 = m 0 , α 0 = m 1 m 0 , ∀k > 0 β k = Z k,k Z k-1,k-1 , α k = Z k,k+1 Z k,k - Z k-1,k Z k-1,k-1 , (13)
which results from Z k+1,k-1 = 0 = Z k+1,k . Thus, if the Z k,p are known for each k ≤ n and p ≤ k, then one can compute the α n , β n from (13) and then the Z n+1,p for p = 0, . . . , n + 1 from [START_REF] Nguyen | Solution of population balance equations in applications with fine particles: Mathematical modeling and numerical schemes[END_REF].

In what follows, it will be also interesting to compute the moments (m k ) k∈{1,...,N } from the (ζ k ) k∈{1,...,N } and from m 0 . It could be done by reversing the previous algorithm, but Skibinsky [42] showed the following theorem, also proved in [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] by another method. Theorem 2.2. Let S i,j be given by S i,j = 0, 0 ≤ j < i, S 0,j = 1, j ≥ 0 and

S i,j = S i,j-1 + ζ j-i+1 S i-1,j , 1 ≤ i ≤ j.
If m n-1 is in the interior of the moment space, then

m n = m 0 n/2 i=0 S 2 i,n-i n-2i j=1 ζ j . (14)
Let us remark that this theorem was proved for measures on [0, 1]. It is easy to see that it is still available for measures on a compact support [0, ξ max ], thanks to relations given in Appendix A. Moreover, a truncated moment vector on [0, +∞) can always be seen as a moment vector of a measure with a compact support [0, ξ max ] (it can be seen by considering the corresponding quadrature), in such a way that this theorem is still valid in our case. This allows us to prove an important new property of the relation between the ζ k and the moments.

Corollary 2.3. The moment m n can be written

m n = m 0 n k=1 ζ k + P n (ζ 1 , . . . , ζ n-1 ) , ( 15
)
where

P n (ζ 1 , . . . , ζ n-1
) is a polynomial function of degree n.

Proof. It is easy to prove, by recursion, that S i,j , for i ≤ j is a polynomial function of degree i of (ζ k ) k=1,...,j-i+1 , in such a way that S 2 i,n-i is a polynomial function of degree 2i of (ζ k ) k=1,...,n-2i+1 . Then the property (15) follows immediately from [START_REF] Attarakih | Solution of the population balance equation using the sectional quadrature method of moments (sqmom)[END_REF].

Let us remark that the term

m 0 P n (ζ 1 , . . . , ζ n-1 ) corresponds to m - k (m k-1 ), since m n -m - n (m n-1 ) = m 0 n i=1 ζ i ,
as shown in [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]. The polynomials P n are given in Appendix C for n = 1, . . . , 7.

Realizable finite volume kinetic schemes

The major issue of the numerical scheme developed here is to ensure the realizability of the vector m N . Indeed, Wright [START_REF] Wright | Numerical advection of moments of the particle size distribution in Eulerian models[END_REF] showed that independent transport of moments with algorithms of order greater than one in space can result in the generation of invalid moment sets.

Let us develop a realizable numerical scheme for the 1D configuration. A generalization on 2D or 3D problems is straightforward for Cartesian meshes thanks to the method of lines.

General form of the finite volume kinetic scheme

Let us introduce a discretization of the spatial domain into cells [x j-1 2 , x j+ 1 2 ] of center

x j = x j-1 2 +x j+ 1 2 2
and of width ∆x = x j+ 1 2 -x j-1 2 . A time discretization t 0 = 0 < t 1 < • • • < t n < . . . is also used, with the time steps ∆t n = t n+1 -t n . The properties of the averaged value of the moments on a cell is used here to define a numerical scheme.

Equations on the averaged value of the moments over a cell

One defines the characteristics X(t; s, y) as the solution of d t X(t; s, y) = u(t, X(t; s, y)), X(s; s, y) = y.

One also defines J(t; s, y) the derivative of y → X(t; s, y), for fixed values of t and s.

Proposition 3.1. For m N = (m 0 , . . . , m N ) t a moment vector solution of (3), the mean value of each moment of order k at time t n , denoted m n k,j = 1 ∆x

x j+ 1 2 x j-1 2 m k (t n , x)dx, satisfies m n+1 k,j = m n k,j - ∆t n ∆x F n k,j+ 1 2 -F n k,j-1 2 , (16)
with

F n k,j+ 1 2 = 1 ∆t n x j+ 1 2 X(t n ;t n+1 ,x j+ 1 2 ) m k (t n , x)dx. (17)
Proof. Let us first remark that ∂ t J(t; t n+1 , x) = J(t; t n+1 , x)∂ x u(t, X(t; t n+1 , x)). Then, the function t → m k (t, X(t; t n+1 , x))J(t; t n+1 , x) is constant, since its derivative is:

[∂ t m k + u∂ x m k + m k ∂ x u] (t, X(t; t n+1 , x))J(t; t n+1 , x) = 0.
Using the values of this function at time t n and t n+1 , we obtain the k th component of the solution of (3) between t n and t n+1 :

m k (t n+1 , x) = m k (t n , X(t n ; t n+1 , x))J(t n ; t n+1 , x).
Then, its averaged value at time t n+1 can be written, using the change of variables ξ = X(t n ; t n+1 , x):

m n+1 k,j = 1 ∆x x j+ 1 2 x j-1 2 m k (t n , X(t n ; t n+1 , x))J(t n ; t n+1 , x)dx = 1 ∆x X(t n ;t n+1 ,x j+ 1 2 ) X(t n ;t n+1 ,x j-1 2 ) m k (t n , ξ)dξ.
This concludes the proof.

This results allows us to develop a kinetic scheme.

Kinetic scheme

Let us denote m n j = (m n 0,j , . . . , m n N,j ) t an approximation at time t n of the averaged value of the moment vector over the cell j and X j+ 1 2 an approximation of X(t n ; t n+1 , x j+ 1 2 ). In what follows, the CFL like number defined by:

CFL = max n,j u n j,max ∆t n ∆x , u n j,max = max u(t, x), x ∈ [x j-1 2 , x j+ 1 2 ], t ∈ [t n , t n+1 ] ( 18
)
is assumed to be smaller than one in such a way that X(t n ; t n+1 , x j+ 1 2 ) is in the cell j or j + 1. Moreover, the same property is assumed for X j+ 1 2 . First, let us define the approximation of X(t n ; t n+1 , x j+ 1 2 ) that will be used here. For a constant and uniform fluid velocity, the exact value is used

X j+ 1 2 = x j+ 1 2 -u∆t n . (19)
Otherwise, one need a good approximation of the velocity and a resolution of the ODE defining the characteristics with an at least second order of accuracy to obtain third order approximation of X(t n ; t n+1 , x j+ 1 2 ). For a stationary fluid velocity, a linear reconstruction of this velocity is done inside the cell. In what follows, we will assume that the fluid velocity is known at the interfaces and its value inside the cell is given by a linear interpolation. This leads to:

X j+ 1 2 = x j+ 1 2 - u(x j+ 1 2 ) δ 1 -e -δ∆t n , (20)
where, denoting

δ j = u x j+ 1 2 -u x j-1 2 ∆x : δ = δ j if u(x j+ 1 2 ) ≥ 0 and δ = δ j+1 if u(x j+ 1 2 ) < 0. Let us remark that X j+ 1 2 < x j+ 1
2 in the first case and X j+ 1 2 > x j+ 1 2 in the second one. For an unstationary fluid velocity, a linear temporal and spatial interpolation of the fluid velocity can be used. Using an explicit second order Runge et Kutta method for the resolution of the ODEs defining the characteristics, one obtains:

X j+ 1 2 = x j+ 1 2 - ∆t n 2 (1 -δ∆t n ) u t n+1 , x j+ 1 2 + u t n , x j+ 1 2 , (21)
where, denoting

δ j = u t n ,x j+ 1 2 -u t n ,x j-1 2 ∆x : δ = δ j if u(t n+1 , x j+ 1 2 ) ≥ 0 and δ = δ j+1 if u(t n+1 , x j+ 1
2 ) < 0. It leads to a third order approximation of X(t n ; t n+1 , x j+ 1 2 ): Lemma 3.2. Let us assume that the CFL number is smaller than one and that the parameter X j+ 1 2 is defined by ( 19), [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gasparticle flows in a LES framework[END_REF] or [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF], depending of the effective dependence of u(t, x) in t and x. Then X j+ 1 2 is a third order approximation of X(t n ; t n+1 , x j+ 1 2 ):

X j+ 1 2 -X(t n ; t n+1 , x j+ 1 2 ) = O(∆x 3 ). Moreover, X j+ 1 2 -X(t n ; t n+1 , x j+ 1 2 ) -X j-1 2 + X(t n ; t n+1 , x j-1 2 ) = O(∆x 4 ) and x j+ 1 2 -X j+ 1 2 - x j-1 2 + X j-1 2 = O(∆x 2 ). Proof. A Taylor expansion for X(t n ; t n+1 , x j+ 1
2 ) leads to:

X(t n ; t n+1 , x j+ 1 2 ) = x j+ 1 2 + t n t n+1 u(t, X(t; t n+1 , x j+ 1 2 ))dt (22) = x j+ 1 2 -u(t n+1 , x j+ 1 2 )∆t n + (∆t n ) 2 2 ∂ t u(t n+1 , x j+ 1 2 )+u(t n+1 , x j+ 1 2 )∂ x u(t n+1 , x j+ 1 2 ) +O((∆t n ) 3 ).
Formula ( 19) is exact for a constant fluid velocity and, since δ is a first order approximation of ∂ x u(t n , x j+ 1 2 ), the Taylor expansion of [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gasparticle flows in a LES framework[END_REF] in the case of a stationary fluid velocity and of (21) in the general case leads to the same zeroth, first and second order terms as in [START_REF] Madadi-Kandjani | An extended quadrature-based moment method with log-normal kernel density functions[END_REF]. Finally, both last results are shown from the same kind of developments, eventually adding the third order term.

Following the property of the exact solution given in Proposition 3.1, the scheme that we consider here is given by the recursion formula:

m n+1 j = m n j - ∆t n ∆x F n j+ 1 2 -F n j-1 2 , (23)
with

F n j+ 1 2 = 1 ∆t n x j+ 1 2 X j+ 1 2 m n (x)dx, (24)
the function m n (x) = (m n 0 (x), . . . , m n N (x)) t being defined from the moments m n j , in a way described in Section 3.2.

Scheme properties

The properties of the scheme will depends on this reconstruction of the moments over the cells. We impose here the two following properties, assuming that m n j is in the moment space M N :

• [P1] The averaged value of the reconstructed moment vector m n (x) over each cell j is m n j :

1 ∆x

x j+ 1 2 x j-1 2 m n (x)dx = m n j .
• [P2] For each x, the vector m n (x) is in the moment space M N .

For the accuracy, a third property is introduced:

• [P3] From an averaged moment vector m j = 1 ∆x x j+ 1 2 x j-1 2 m(x)dx, the reconstruction m(x) is
second order accurate: there exists a bounded function ϕ m such that

m(x) = m(x) + ∆x 2 ϕ m (x), (25)
with, for all j:

∀δ ∈] -∆x, ∆x[ D j (ϕ m ) def = x j+ 1 2 +δ x j+ 1 2 ϕ m (x)dx - x j-1 2 +δ x j-1 2 ϕ m (x)dx = O(∆x 2 ).
This leads to the following theorem:

Theorem 3.3.
Let us assume that the CFL number is smaller than one. The finite volume scheme defined by [START_REF] Shohat | The Problem of Moments[END_REF][START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] with the properties [P1] and [P2] of the moment reconstruction is realizable.

With the additional property [P3] and using the value of X j+ 1 2 defined by ( 19), ( 20) or ( 21), then a second order is obtained for the consistency error.

Proof. The first property [P1], with the definition (24) of the fluxes implies:

m n+1 j = 1 ∆x X j+ 1 2 X j-1 2 m n (x)dx. ( 26
)
Since the characteristics cannot cross themselves, one have X j-1 2 < X j+ 1 2 . Then [P2] ensures that m n+1 j is in the moment space. Moreover, from the exact solution m = (m k ) k∈{0,1,...,N } of (3), the consistency error is defined by

n j = m n+1 j -m n j ∆t n + F n j+ 1 2 -F n j-1 2 ∆x ,
where

m n+1 j = (m n k,j ) k∈{0,1,...,N } is the exact averaged moment vector, m n k,j = 1 ∆x x j+ 1 2 x j-1 2 m k (t n , x)dx,
and the flux vector

F n j+ 1 2 = 1 ∆t n x j+ 1 2 X j+ 1 2 m n (x)dx is defined with the reconstruction m n (x) obtained
from the m n j . Thanks to Proposition 3.1, this can be written:

n j = 1 ∆t n ∆x [I 1 + I 2 ] ,
with

I 1 = x j+ 1 2 X j+ 1 2 (m n (x) -m(t n , x))dx - x j-1 2 X j-1 2 (m n (x) -m(t n , x))dx and I 2 = X t n ;t n+1 ,x j+ 1 2 X j+ 1 2 m(t n , x)dx - X t n ;t n+1 ,x j-1 2 X j-1 2 m(t n , x)dx.
Thanks to [P3], the first term I 1 can be written:

I 1 = ∆x 2   x j+ 1 2 X j+ 1 2 ϕ m(t n ,.) (x)dx - x j-1 2 X j-1 2 ϕ m(t n ,.) (x)dx   = -∆x 2   x j+ 1 2 +δ x j+ 1 2 ϕ m(t n ,.) (x)dx - x j-1 2 +δ x j-1 2 ϕ m(t n ,.) (x)dx - x j+ 1 2 +X j-1 2 -x j-1 2 X j+ 1 2 ϕ m(t n ,.) (x)dx   , with δ = X j-1 2 -x j-1 2 , which is in ] -∆x, ∆x[.
Then, thanks to the property [P3] and Lemma 3.2, I 1 = O(∆x 4 ). The second term I 2 can be written:

I 2 = m(t n , X j+ 1 2 ) X t n ; t n+1 , x j+ 1 2 -X j+ 1 2 -m(t n , X j-1 2 ) X t n ; t n+1 , x j-1 2 -X j-1 2 +O(∆x 6 ) = m(t n , X j+ 1 2 ) X t n ; t n+1 , x j+ 1 2 -X j+ 1 2 -X t n ; t n+1 , x j-1 2 +X j-1 2 + O(∆x 4 ).
The result of Lemma 3.2 allows to conclude the proof.

Let us remark that a classical MUSCL type reconstruction for each moment separately satisfies the properties [P1] and [P3] except around extrema of the considered moments (see Appendix D), but this would inevitably lead to unrealizable moments m n (x), due to the complex shape of the moment space [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF], in which case the property [P2] would not be satisfied.

Reconstruction of the moments in the cell and flux computation

To complete the scheme [START_REF] Shohat | The Problem of Moments[END_REF][START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF] a spatial reconstruction of the moments from their mean value is provided, verifying at least properties [P1] and [P2]. Different kinds of such reconstructions can be found in the literature. The constant reconstruction is the simplest, but leads to low order schemes. It is given in Section 3.2.1.

Another one is based on a direct reconstruction at the level of the NDF. This was done in the completely different context of linear kinetic equations in slab geometry [START_REF] Garrett | A Comparison of Moment Closures for Linear Kinetic Transport Equations: The Line Source Benchmark[END_REF][START_REF] Schneider | A realizability-preserving high-order kinetic scheme using WENO reconstruction for entropy-based moment closures of linear kinetic equations in slab geometry[END_REF]. In this context, the internal variable had a compact support and the advection velocity depended on the internal variable ξ in such a way that the moment equations had to be closed; this was done through the reconstruction of the NDF from the moments by entropy maximization [START_REF] Mead | Maximum entropy in the problem of moments[END_REF]. Then, for each value of the internal variable, a spatial MUSCL type reconstruction of this NDF was used. And a quadrature was used to compute the fluxes, in such a way that the spatial reconstruction had to be done for a finite number of values of the internal variable. This is applied in our context in a quite different way presented in Appendix F, using EQMOM [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF][START_REF] Nguyen | Solution of population balance equations in applications with fine particles: Mathematical modeling and numerical schemes[END_REF] instead of the entropy maximization, which is not really adapted to NDF defined on [0, +∞). However, it does not allow to deal with the boundary of the moment space, since the NDF has to be a real function, and it induces a high supplementary cost related to the reconstruction (the computation time is then multiplied by at least 65 when using 5 moments, compared to the schemes described in this section). Moreover, it is also shown in Appendix F that the maximum principle on the moments can be lost.

Still in the case of compact support, the moments can be spatially reconstructed through a MUSCL type reconstruction of the canonical moments (see Appendix A for their definition and their link with the ζ k ). This was done in the context of population of inertial particles [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF], their velocity belonging to the internal variables and being then treated differently compare to our application. Moreover, this reconstruction was limited to a model considering only four moments. In Section 3.2.2, since we are interested in moments on [0, +∞), instead of the canonical moments, the variables (ζ k ) k∈{1,...,N } are reconstructed. Moreover, as in [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF], the property [P1] is not directly verified, since the reconstructed moments are no more affine in the cell and the same kind of strategy is used to recover it, with eventually additional corrections here to ensure robustness, even for moment vectors close to the boundary of the moment space. And here, the reconstruction is done for any number of moments, eventually at the boundary of the moment space.

The last type of spatial reconstruction of the moments was done in [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF]. It used a spatial reconstruction of the only weights of the quadratures associated with the moment vectors, and was used in another context, the internal variable being the velocity and evolving in R. However, it was limited to an even number of moments, in the interior of the moment space and led to only a pseudo-second-order of accuracy, due to the fact that the variation of the abscissas was not taken into account in the reconstruction. In Section 3.2.3, a variant of this reconstruction is introduced, adapted to any number of moments, eventually at the boundary of the moment space and allowing a real second-order of accuracy for half the lowest order moments.

The constant reconstruction based kinetic scheme

From the mean value m n k,j of the moment in the cell, the first idea is to consider a constant value of the moments into the cell:

m n (x) = m n j for x ∈]x j-1 2 , x j+ 1 2
[. The flux can then be written:

F n j+ 1 2 = m n j ∆t n x j+ 1 2 -X j+ 1 2 + - m n j+1 ∆t n X j+ 1 2 -x j+ 1 2 + .
where u + = max{0, u}. Its convergence with a first order of accuracy is easily shown with the value of X j+ 1 2 defined by ( 19), [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gasparticle flows in a LES framework[END_REF] or [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF].

The ζ reconstruction based kinetic scheme (ζ kinetic scheme)

Let us denote (ζ n k,j ) k∈{1,N } the ζ k corresponding to the moment vector m n j . For a fixed value of n, instead of reconstructing directly m n (x), the corresponding (ζ k (x)) k∈{1,...,N } as well as m n 0 (x) are reconstructed. For m n 0 , a classical MUSCL type of reconstruction is used:

m n 0 (x) = m n 0,j + D n 0,j (x -x j ) for x ∈]x j-1 2 , x j+ 1 2 [. For each k, a piecewise affine reconstruction ζ n k (x) is defined by: ζ n k (x) = ζn k,j + D n k,j (x -x j ) for x ∈]x j-1 2 , x j+ 1 2 [. Each ζn k,j is computed to ensure [P1]. Thanks to Corollary 2.3, this leads to: ∆x m n k,j = ζn k,j x j+ 1 2 x j-1 2 m n 0 (x) k-1 i=1 ζ n i (x)dx + D n k,j x j+ 1 2 x j-1 2 (x -x j )m n 0 (x) k-1 i=1 ζ n i (x)dx + x j+ 1 2 x j-1 2 m n 0 (x)P k (ζ n 1 (x), . . . , ζ n k-1 (x))dx.
Let us write ζn k,j = a k,j + b k,j D n k,j , where a k,j and b k,j are defined by:

a k,j x j+ 1 2 x j-1 2 m n 0 (x) k-1 i=1 ζ n i (x)dx = ∆x m n k,j - x j+ 1 2 x j-1 2 m n 0 (x)P k (ζ n 1 (x), . . . , ζ n k-1 (x))dx, ( 27
) b k,j x j+ 1 2 x j-1 2 m n 0 (x) k-1 i=1 ζ n i (x)dx = - x j+ 1 2 x j-1 2 (x -x j )m n 0 (x) k-1 i=1 ζ n i (x)dx. ( 28
)
Let us remark that |b k,j | is smaller than ∆x 2 . Numerically, the term

P k (ζ n 1 (x), . . . , ζ n k-1 (x)
) is computed by using the reverse algorithm (see Section 2.4), setting ζ k = 0. Moreover, the integrals, which are integrals of polynomial functions, are computed using a Gauss-Legendre quadrature with N 2 + 1 points, which is enough to obtain their exact value.

For m 0 , a minmod limiter is used:

D n 0,j = minmod m n 0,j+1 -m n 0,j ∆x , m n 0,j -m n 0,j-1 ∆x . (29)
Moreover, as in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF], the slopes

D n k,j are such that ζ n k (x) stays between min{ζ n k,j-1 , ζ n k,j , ζ n k,j+1 } and max{ζ n k,j-1 , ζ n k,j , ζ n k,j+1
} in each cell j, as soon as it is the case for a k,j . One can use [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF]:

D n k,j =        min |ζ n k,j+1 -a k,j | ∆x + 2b k,j , |a k,j -ζ n k,j-1 | ∆x -2b k,j , 2a k,j ∆x -2b k,j if ζ n k,j-1 < ζ n k,j < ζ n k,j+1 , -min |ζ n k,j+1 -a k,j | ∆x + 2b k,j , |a k,j -ζ n k,j-1 | ∆x -2b k,j , 2a k,j ∆x + 2b k,j if ζ n k,j-1 > ζ n k,j > ζ n k,j+1 . (30)
The last condition (in the min) ensures the non-negativity of the reconstruction, as soon as a k,j is non-negative. Unfortunately, a k,j can eventually be negative, meaning that, even with a constant reconstruction for ζ n k (x), one can obtain unrealizable moments. In this case, a correction has to be done on the slopes D n i,j , with i < k: they are reduced successively for i = k 0 , ..., k -1 where k 0 is the maximal value such that a k,j is positive if the D n i,j were equal to zero for i ∈ [k 0 , k -1]. For that, starting with i = k 0 till i = k -1, each D n i,j is multiplied by 0.9 once or several times (at most 5 times here) and if it is not sufficient, it is set to zero. In the worst case, all slopes D n i,j are then set to zero, for i ≥ 1. But when needed, only one limitation is often sufficient. Let us remark that the consequence of such kind of a posteriori corrections, if one or several slopes is set to zero, is the reduction of the polynomial degree of m k (x), as done in other kind of method such as MOOD [START_REF] Clain | A high-order finite volume method for systems of conservation lawsmulti-dimensional optimal order detection (mood)[END_REF], in a quite different way. The algorithm detailing the reconstruction procedure and especially the corrections ensuring the positivity of a k,j is given in Appendix E.

The verification of property [P3] for this reconstruction is complex. Thanks to the polynomial form of the functions, it can be shown for regular initial conditions that, except near the boundary of the moment space (where one of several

ζ n i,j is small), a k,j = ζ n k,j + O(∆x 2 ) and b k,j = O(∆x 2 ), in such a way that ζn k,j = ζ n k,j + O(∆x 2
). This means that, for a small enough value of ∆x and for moments far enough from the boundary of the moment space, the parameter a k,j should be positive in such a way that no correction has to be done. That is why this case was not encountered in [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF] where the boundary of the moment space were not dealt with. In this context, i.e. far from the boundary of the moment space, one can then show that, except near any extrema of the function

ζ k (x), m(x) = m(x) + ∆x 2 ϕ m (x)
, where ϕ m is bounded. However, the additional condition on ϕ m is hard to show, even if we guess that it is valid, except near extrema of any ζ k or near the boundary of the moment space. The order of accuracy of the scheme will be checked numerically on some examples.

The fluxes can then be written, thanks to Corollary 2.3:

F n k,j+ 1 2 = 1 ∆t n x j+ 1 2 X j+ 1 2 m n 0 (x)   k j=1 ζ n j (x) + P k (ζ n 1 (x), . . . , ζ n k-1 (x))   dx, (31)
They are computed using a Gauss-Legendre quadrature with N 2 + 1 points, which gives their exact value since one has to integrate polynomial of degree at most N + 1. The corresponding scheme is called "ζ reconstruction based kinetic scheme", abbreviated to "ζ kinetic scheme". Let us remark that the resolution of the equation on m 0 does not depend on the other moments, even if some corrections are needed for the positivity of the a k,j . Moreover, the corresponding scheme is TVD for the constant velocity case. It is also easy to see that a maximum principle is obtained for

ζ 1 = m 1 /m 0 in this case.

The quadrature weights reconstruction based kinetic scheme (QW kinetic scheme)

Similarly to what was done in [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF], the reconstruction is based here on the quadrature of the moment sets, but in a different way. Let us then define the quadratures weights (w α,j ) α∈{1,...,p} and abscissas (ξ α,j ) α∈{1,...,p} of the moment set (m n k,j ) k∈{0,...,N } , such that:

m n k,j = p α=1 w α,j ξ k α,j , k ∈ {0, . . . , N }, with p = N 2 + 1.
If the number of moments is even, it is the classical Gauss quadrature with N = 2p -1. Then, abscissas and weights are deduced from the eigenvalues and eigenvectors of the Jacobi matrix with coefficients α k and √ β k [START_REF] Press | Numerical recipes[END_REF], where α k and β k are the coefficient of the three term recurrence relation [START_REF] Barrett | A comparison of some approximate methods for solving the aerosol general dynamic equation[END_REF]. Otherwise, N = 2(p -1) and a Gauss-Radau quadrature is used, meaning that an abscissa is set to zero: ξ j,1 = 0. They are computed in the same way as for the classical Gauss quadrature, except that α p-1 is set to -β p-1 Πp-1(0) Πp(0) [START_REF] Press | Numerical recipes[END_REF]. Like in [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF], the reconstruction of the moments inside the cell j is done through a reconstruction of the weights:

m n k (x) = p α=1 w α (x)ξ k α,j , w α (x) = w α,j + D n α,j (x -x j ), x ∈ [x j-1 2 , x j+ 1 2 ].
The averaged value of m n k (x) on the j th cell is then automatically m n k,j . But to define the slopes, unlike in [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF], one uses here a quadrature of the neighbor cells with the same abscissas:

m n k,j-1 = p α=1 w - α,j ξ k α,j , m n k,j+1 = p α=1 w + α,j ξ k α,j , k ∈ {0, . . . , p -1}. ( 32
)
This allows also to consider moments at the boundary of the moment space, where p is then reduced, which was not possible with the method used in [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF], where the case of the neighboring cells with a representation by a quadrature with a different number of abscissas could not be taken into account easily. Here, the weights (w - α,j ) α∈{1,...,p} and (w + α,j ) α∈{1,...,p} are well defined since they are the solution of a linear system with a Vandermonde matrix of coefficients ξ α,j , which are distinct. If a regular initial condition is considered for (3) and if (m n k,j ) k∈{0,...,N } is far from the boundary of the moment space, then these weights (w - α,j ) α∈{1,...,p} and (w + α,j ) α∈{1,...,p} are usually non-negative since the moment sets for j -1, j and j + 1 are close to each other and the weights w α,j are far from zero. But if it is not the case, the most negative w + α,j or w - α,j is set to zero and the corresponding linear system in [START_REF] Curto | Recursiveness, positivity, and truncated moment problems[END_REF] is solved for k ∈ {0, . . . , p -2}, thus eliminating the corresponding abscissa ξ α,j . This operation is reproduced till all the w ± α,j are non-negative. The corresponding algorithm is detailed in Appendix E. At worst, only one abscissa will stay in [START_REF] Curto | Recursiveness, positivity, and truncated moment problems[END_REF]. But for the distributions used in this paper, when needed, only one or two abscissas has to be eliminated when using 10 moments. One then defines the slopes by using a minmod limiter:

D n α,j = 1 2 sgn(w + α,j -w α,j ) + sgn(w α,j -w - α,j ) min |w + α,j -w α,j | ∆x , |w α,j -w - α,j | ∆x .
In order to verify the property [P3] let us rewrite this slope in the following way, in the same way as in Appendix D:

D n α,j = w + α,j -w α,j ∆x Φ(θ αj ), θ α,j = w α,j -w - α,j w + α,j -w α,j
.

The function Φ then corresponds to the classical minmod flux limiter [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]: Φ(θ) = minmod(1, θ).

For the k th order moment, this corresponds to the following reconstruction, for x ∈ (x j-1 2 , x j+ 1 2 ):

m n k (x) = m n k,j + (x -x j )D n k,j , D n k,j = p α=1 D n α,j ξ k α,j = p α=1 w + α,j -w α,j ∆x ξ k α,j Φ(θ α,j ).
Thanks to the definition of the weights, the corresponding slope can be rewritten, if k ≤ p -1:

D n k,j = m n k,j+1 -m n k,j ∆x + p α=1 w + α,j -w α,j ∆x ξ k α,j [Φ(θ α,j ) -1] .
With only the first term of the slope, the reconstruction of the moment would correspond to a second order approximation without any slope limiter. Moreover, let us denote w α,j (x) the weights of a quadrature corresponding to the analytical solution with imposed abscissas (ξ α,j ) k∈{1,...,p} : m k (x) = p α=1 w α,j (x)ξ k α,j for k = 0, . . . , p -1. Except near the extrema of this weights, the term Φ(θ α,j )-1 is first order accurate. Then, in this case, the reconstruction m n k (x) from the analytical averaged moments is of the form given in

[P3], m k (x) = m k (x)+∆x 2 ϕ m k (x). Moreover, for example if δ ∈]0, ∆x[: D j (ϕ m k ) = O(∆x 2 )+ δ(∆x-δ) 2∆x 2 p α=1 w + α,j+1 -w α,j+1 ∆x ξ k α,j+1 [Φ(θ α,j+1 )-1]- p α=1 w + α,j -w α,j ∆x ξ k α,j [Φ(θ α,j )-1] .
Thanks to the use of the w α,j (x), the previous expression can be written:

D j (ϕ m k ) =O(∆x 2 ) + δ(∆x-δ) 2∆x 2 p α=1 w α,j+1 (x j+ 1 2 ) -w α,j (x j+ 1 2 ) ξ k α,j+1 [Φ(θ α,j+1 )-1] + p α=1 w α,j (x j+ 1 2 ) ξ k α,j+1 -ξ k α,j [Φ(θ α,j+1 )-1] + p α=1 w α,j (x j+ 1 2 )ξ k α,j+1 [Φ(θ α,j+1 )-Φ(θ α,j )] . .
Except near the extrema of the weight w α,j (x), the term Φ(θ α,j+1 )-Φ(θ α,j ) is second order accurate. Moreover, except near the boundary of the moment space, the quadrature reconstruction is stable, in such a way that ξ k α,j+1 -ξ k α,j and w α,j+1 (x j+ 1 2 ) -w α,j (x j+ 1 2 ) are then also first order terms. Thus, except near the boundary of the moment space and near the extrema of the quadrature weights, the property [P3] is checked for the moments of order 0 to p -1.

Finally, the fluxes are written:

F n j+ 1 2 = p α=1 w α,j + 1 2 D n α,j (X j+ 1 2 -x j-1 2 ) Ξ α,j ∆t n x j+ 1 2 -X j+ 1 2 + - p α=1 w α,j+1 -1 2 D n α,j+1 (x j+ 3 2 -X j+ 1 2 ) Ξ α,j+1 ∆t n X j+ 1 2 -x j+ 1 2 + .
with Ξ α,j = (1, ξ α,j , . . . , ξ N α,j ) T . This scheme is called "quadrature weights (QW) reconstruction based kinetic scheme", abbreviated to "QW kinetic scheme". Let us remark that, since the limitation is done separately for each weight, the global reconstruction can not guarantee a maximum principle on any moments. One can however remark that the ξ α,j are bounded by the maximal and minimal values of the abscissas corresponding to the initial moments. Indeed, the new moments m n+1 j correspond to a sum of weighted Dirac delta functions at ξ α,j and ξ j±1,α with α = 1, . . . , p. Then, the new abscissas are necessarily inbetween these values.

Dealing with the boundary of the moment space

It can happen that, at some point, the moment vector is at the boundary of the moment space, with eventually for one of the neighbor cells, a moment vector at another boundary (different value of N (m N )) or in the interior of the moment space. In this case, the developed schemes can locally loose their second order of accuracy. But the key point is their realizability, which can be very difficult to preserve in this case, which was not taken into account in [START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF] or [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF]. However, here, the developed schemes are still realizable if some precautions are taken. For the ζ reconstruction, the ζ k for k > N (m N ) are not defined but are set to zero and eventually some corrections are done on the slopes to guaranty the positivity of a k,j , as explained in Section 3.2.2. For the QW reconstruction, the value of p has just to be reduced to N (m N )+1 2 . Let us remark that for the quadrature, such kind of adaptation of the number of weights were introduced [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF], but here it is directly based on the evaluation of N (m n N ), through the computation of the ζ k , as for the ζ reconstruction. Since they are deduced from the coefficients of the three term recurrence relation used for the computation of the quadrature thanks to [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF], they are computed in any cases, with at worst a marginal increase of the cost.

Numerically, an additional difficulty comes from the detection of the boundary of the moment space, i.e. the evaluation of N (m N ), the computation of the ζ k being ill-conditioned near this boundary. Indeed, when a ζ k is positive but very small, the computation of the next ζ i can give anything. Then, a small parameter is introduced here, and N (m N ) is replaced by N (m N ) defined as the minimum value of k such that ζ k < (set to N + 1 if all the ζ k are greater than ). Moreover, for i ≥ N (m N ), ζ i is then set to zero and the corresponding vector is then projected: the moments are computed from the new set of (ζ k ) k∈{1,...,N } . In practice, we used in this work = 10 -7 .

A realizable simplified finite volume scheme

The accuracy of the kinetic schemes comes from the accuracy of the spatial reconstruction as well as the correct evaluation of the characteristics, to characterize the part of the cell containing the particles that will be transferred to another one during the time step. For this last point, a generalization to multi-dimensional unstructured meshes is not evident. In this section, we then consider semi-discretized finite volume schemes of the form:

d t m j (t) = - 1 ∆x F j+ 1 2 (t) -F j-1 2 (t) , F j+ 1 2 (t) = u(t, x j+ 1 2 )m j+ 1 2 (t), (33)
where only the reconstructed value m j+ 1 2 (t) of the moment at the interface x j+ 1 2 is considered. The system [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] will then be solved thanks to a strong stability preserving (SSP) explicit Runge-Kutta method [START_REF] Gottlieb | Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations[END_REF], a second order one here. It is a convex combination of Euler explicit time steps. The realizability then only have to be shown for the explicit Euler method:

m n+1 j = m n j - ∆t n ∆x F n j+ 1 2 -F n j-1 2 , (34)
with

F n j+ 1 2 = u(t n , x j+ 1 2 )m n j+ 1 2
. Moreover, to deal with the boundary of the moment space, a numerical projection is also done after each of these Euler explicit time steps, as described in Section 3.2.4.

Two types of reconstruction are done to define the moment vector m j+ 1 2 at the interface, as for the kinetic scheme: a first one using the ζ k variables and a second one using the weights of the quadrature. The realizability of the scheme is then shown in each case.

The ζ reconstruction based simplified scheme (ζ simplified scheme)

Here, the reconstruction is quite different from the one of Section 3.2.2 in order to be able to ensure the realizability with the explicit Euler method. Indeed, from the idea of Berthon [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF], the cell j is split in three parts of size ∆x/3 here, on which the reconstruction of the moment vector is constant. It is denoted m - j for the left part, m * j for the middle part and m + j for the right part, in such a way that:

m n j = 1 3 m - j + m * j + m + j .
The values of the ζ k corresponding to the moment vectors m n j , m - j , m * j and m + j are denoted respectively ζ n k,j , ζ - k,j , ζ * k,j and ζ + k,j . Let us remark that Berthon [START_REF] Berthon | Stability of the MUSCL schemes for the Euler equations[END_REF] used this kind of reconstruction in a completely different context of the Euler equations (i.e. considering moments of order 0, 1 and 2 in velocity, on R). This has then to be adapted in our context of high order moments in size. Moreover, here, the role of m * j is just to guaranty the realizability of the scheme by doing a reconstruction of realizable moments in the entire cell (property [P2]). But this kind of reconstruction is not used in the context the kinetic scheme of Section 3, since it would necessarily induce a loss of accuracy, the reconstruction being only accurate at the bounds of the cell.

The reconstruction is then done in the following way. For the left and right values, a classical MUSCL type of reconstruction is used:

m ± 0 = m n 0,j ± D n 0,j ∆x 2 , ζ ± k,j = ζ n k,j ± D n k,j ∆x 2 , k ∈ {1, . . . , N }, (35)
with the minmod limiters:

D n 0,j = minmod m n 0,j+1 -m n 0,j ∆x , m n 0,j -m n 0,j-1 ∆x , (36)
D n k,j = minmod ζ n k,j+1 -ζ n k,j ∆x , ζ n k,j -ζ n k,j-1 ∆x . (37)
Then, the values of m + j and m - j can be computed. However, the vector m * j , given by 3m n jm + j -m - j , is not necessarily realizable. In this case, some values of the D n k,j , for k > 0, have to be reduced. To ensure the realizability of (m * 0,j , m * 1,j ), the following limitation on D n 1,j has to be added:

D n 0,j D n 1,j < 2 ∆x 2 m n 1,j . (38)
For the higher order moments, it is harder to derive analytical formulas. Instead, if m * j is not in the moment space, a correction has to be done on the slopes D n i,j , with i < k: they are reduced successively for i = k 0 , ..., k-1 where k 0 is the maximal value such that (m * 0,j , . . . , m * k,j ) is realizable if the D n i,j were equal to zero for i ∈ [k 0 , k -1]. For that, starting with i = k 0 till i = k -1, each D n i,j is divided by 2 or set to zero if it is not sufficient. The corresponding algorithm is detailed in Appendix E.

Once the reconstruction is done, the value of the fluxes is obtained by

F n j+ 1 2 = max{u(t n , x j+ 1 2 ), 0}m + j + min{u(t n , x j+ 1 2 ), 0}m - j+1 . (39)
The corresponding scheme is called "ζ reconstruction based simplified scheme", abbreviated to "ζ simplified scheme". Under some limitation on the CFL number, it is realizable, as shown in the following theorem.

Theorem 4.1. From realizable moment vectors (m n j ) j , the equations (34,39) define a realizable moment vector m n+1 j if the CFL number is smaller than 1/3.

Proof. The vector m n+1 j is given by:

m n+1 j = m n j -λ j m + j -μj-1 m - j + µ j m - j+1 + λj-1 m + j-1 ,
where λ j = ∆t n ∆x max{u(t n , x j+ 1 2 ), 0}, µ j = -∆t n ∆x min{u(t n , x j+ 1 2 ), 0}, λj-1 = ∆t n ∆x max{u(t n , x j-1 2 ), 0}, μj-1 = -∆t n ∆x min{u(t n , x j-1 2 ), 0} are coefficients between 0 and 1/3. The last two terms define some moment vectors and the rest can be written 1 3 -λ j m + j + 1 3 -μj-1 m - j + 1 3 m * j and is also realizable.

The quadrature weights reconstruction based simplified scheme (QW simplified scheme)

The same reconstruction of the weights is done here as for the kinetic scheme in Section 3.2.3. Using the same notation, the flux is then given by:

F j+ 1 2 (t) = max{u(t, x j+ 1 2 ), 0} p α=1 w α,j + ∆x 2 D n α,j Ξ α,j (40)
+ min{u(t, x j+ 1 2 ), 0}

p α=1 w α,j+1 - ∆x 2 D n α,j+1 Ξ α,j+1 .
The corresponding scheme is called "quadrature weights (QW) reconstruction based simplified scheme", abbreviated to "QW simplified scheme" and is then realizable under some less restrictive limitation on the CFL number than for the ζ-simplified scheme.

Theorem 4.2. From realizable moment vectors (m n j ) j , the equations (34,40) define a realizable moment vector m n+1 j if the CFL number is smaller than 1/2.

Proof. The vector m n+1 j can be written:

m n+1 j = p α=1 w α,j -λ j w α,j + ∆x 2 D n α,j -μj-1 w α,j - ∆x 2 D n α,j Ξ α,j + µ j p α=1 w α,j+1 - ∆x 2 D n α,j+1 Ξ α,j+1 + λj-1 p α=1 w α,j-1 + ∆x 2 D n α,j-1 Ξ α,j-1 ,
with the same definition of λ j , µ j , λj-1 and μj-1 as in the proof of Theorem 4.1 but now belonging to [0, 1/2]. The only thing to check is the non-negativity of the weights in the first summation. If w α,j is zero, then D n α,j = 0 and this term is zero. Otherwise, w α,j is positive and this term can be written:

w α,j 1 -1 + ∆x 2w α,j D n α,j λ j -1 - ∆x 2w α,j D n α,j μj-1 , which is non-negative, since necessarily |D n α,j | ≤ 2wα,j
∆x and λ j and μj-1 are smaller than 1/2.

Verification: comparisons with analytical solution in 1D configurations

Two configurations are studied here for which analytical solutions are available. For the first one, moments are transported by a constant fluid velocity in a periodic domain, whereas for the second one, moments are transported by an unsteady and compressible fluid velocity. In both cases, three initial conditions can be considered, defined either from the NDF or from the ζ k : a regular one with moments in the interior of the moment space, an oscillating one where the ζ k oscillate at different frequencies to test the robustness of the methods and a multi-modal one, reaching the boundary of the moment space and typical of what can be obtained in physical systems were nucleation and aggregation occur. These initial conditions are first described before comparing the results of our schemes to the analytical solution in both test cases. Then, the schemes are tested in 1D configurations with a constant or a compressible fluid velocity, considering a high number of moments, equal to 10 unless mentioned otherwise.

Initial moments 5.1.1. Regular initial NDF

A regular initial condition is defined through the following beta-NDF:

f 0 (ξ, x) = 16x 2 (1 -x) 2 ξ λ(x) (1 -ξ) µ(x) β(λ(x), µ(x)) , with λ(x) = 7 2 + 3 2 sin(2πx), µ(x) = 7 2 - 3 2 cos(2πx)
and the initial moments are defined by m 0 k (x) =

Oscillating initial ζ k

An initial solution can be also defined from a choice of the ζ k and of m 0 . Here, oscillating functions are used for ζ k , with a frequency increasing with k, whereas m 0 is a polynomial function:

ζ 0 k (x) = x 2 1.01 + cos πk 2 x , m 0 0 (x) = 16x 2 (1 -x) 2 . ( 41
)
This allows to recover all the initial moments, which are in the interior of the moment space, but quite close to the boundary at some points. To test the ability of the methods to deal with the boundary of the moment space and the transition with the interior of this space, another initial distribution is introduced. It is detailed in Appendix G and the corresponding moments are regular (C 2 ) functions of x. Then, the distribution is only one Dirac delta function at ξ = ξ 1 = 0.02 for x ∈ 0, 1 4 , in such a way that only ζ 1 is not zero, as seen in Fig. 2(left). It represents a distribution obtained through nucleation. For x ∈ 1 4 , 1 3 , it is a sum of two Dirac delta functions at ξ = ξ 1 and ξ = 2ξ 1 and ζ 2 and ζ 3 are positive, with a discontinuity at x = 1 3 for ζ 2 . The second peak represents particles obtained by aggregation of two initial nuclei. A continuous Rosin-Rammler distribution is added for the rest of the domain, in such a way that the moments are then in the interior of the moment space. The obtained distribution is plotted in Fig. 2(right) for three different positions. Moreover, the corresponding moments are plotted in Fig. 7.

Results with a constant fluid velocity

To test the different schemes, let us first consider a case with a constant velocity. With no restriction, one assumes u = 1. Moreover, we consider the spatial domain [0, 1] with periodic boundary conditions. The analytical solution at time t = 2 or t = 5 is then equal to the initial condition.

Numerical accuracy in a regular case

When considering the regular initial NDF, with moments far from the boundary of the moment space, simulations are done with the three kinetic schemes on uniform meshes with a CFL number equal to 0.8 and are compared with the analytical solution at the final time t = 2. First, using a 100 points spatial discretization, the moments of order 0 to 6 obtained at the final time with the first order scheme are plotted in Fig. 1(left), showing its numerical diffusion whereas the other kinetic schemes lead to very precise results (not shown here).

-

logarithm of the error logarithm of the cell width -

logarithm of the error logarithm of the cell width To evaluate more precisely the accuracy of the methods, the L 1 norm of the errors, divided by the L 1 norm of the corresponding moments, are plotted as a function of the cell width in Fig. 3, using from 16 to 4096 cells. The L 1 norm is used since we saw that a local loss of accuracy can happen near extrema of the reconstructed variables or near zones where the moments are at the boundary of the moment space. Both the methods are numerically almost second order accurate: the order is about 1.93 for all moments with the ζ kinetic scheme and 1.91 for the moments of order 0 to 5 with the QW kinetic scheme. Moreover, the ζ kinetic scheme does not introduce any dependance on the number of considered moments in this case where no corrections are needed on the slope to ensure the positivity of the a k,j . On the contrary, for the QW kinetic scheme, the number of considered moments has an influence on the accuracy on all moments. It can be seen in Fig. 4, where the same normalized L 1 norm of the errors are plotted when using 6 and 5

-6 -5 -4 -3 -2 -1 -4 -3.5 -3 -2.5 -2 -1.5 -1
logarithm of the error logarithm of the cell width - moments with this scheme. Moreover, as explained in Section 3.2.3, the second order of accuracy is only guarantied for half the first order moments and indeed, here, one can see that for half of the highest order moments, the order of accuracy degenerates to 1 when a fine enough mesh is used (this degenerescence did not yet really appear for the finest mesh used in the case with 10 moments). Let us also remark that, for this case, no correction is needed to ensure the positivity of the weights of the reconstruction (i.e. no elimination of weights for the neighbor cells), as expected since the moments are far from the boundary of the moment space.

When considering the simplified schemes, the CFL number is set to 0.3. The behavior of the schemes are similar (and not represented here): both the methods are also numerically almost second order accurate: the order is about 1.93 for all moments with the ζ simplified scheme and about 1.9 for the moments of order 0 to 6 with the QW simplified scheme. They are also compared to the kinetic schemes, using the same CFL equal to 0.3. The errors on the moments of order 0 and 1 are then plotted in Fig. 5. One can first remark that the accuracy of the kinetic schemes is smaller than for a CFL equal to 0.8, with more than half an order of magnitude for the difference.

It can be checked in this case that the results for the simplified schemes are very little sensitive to the CFL number. Moreover, the accuracy obtained with the non-constant reconstructions of the moments are much higher than with the constant reconstruction, allowing to reduce drastically the number of cells needed for a given accuracy (about a factor 10 on the cell number for an error of 1%). This is an important point for more complex problems where a quite costly ODE system has to be solved on each cell for the resolution of the source terms. Moreover, the ζ reconstruction gives the best results for most of the discretizations. Let us consider more challenging test cases. The final time, t = 5, is larger here to amplify the phenomena. The oscillating initial ζ k is first considered. One can see that the ζ schemes, using 100 points and a CFL number equal to 0.3, give accurate results on the moments of order 0 to 4 (see Fig. 6 middle). However, for the highest order moments, the maximum principle is no more respected and the error is larger than for the first order moments but still smaller than with the other reconstructions, at least for the kinetic scheme. With the QW schemes using 10 moments and the same discretization (Fig. 6 bottom), the results are a little bit less good: the maximum principle is no more respected for the first order moments with some kind of oscillations around the analytical solution. However, these oscillations do not blow up and the accuracy is still better than with the first order scheme, for which the numerical diffusion flatten a lot the results (Fig. 6 top). Moreover, the accuracy of the QW schemes for the first order moments is improved by the use of a smaller number of moments. Computations were done also with the multi-modal initial NDF, representing a more realistic NDF. The analytical solution as well as the solution computed with the ζ and QW schemes with 100 points and a CFL number equal to 0.3 are represented in Fig. 7. Both schemes do not encounter any problem simulating this test case, even if the moment vector is at the boundary of the moment space in some part of the domain. Moreover, in this more realistic case, the maximum principles are respected and the ζ reconstruction gives the most accurate results when using 10 moments, but the use of only 4 moments allow to improve the accuracy for the QW schemes. As it will be shown for the unsteady case, the difference of accuracy on the zeroth and first order moments is higher between the ζ and QW reconstructions when a finer discretization is used for simulations with 10 moments and a second order of accuracy is still recovered. Moreover, all methods are still always much more accurate than the schemes using a constant reconstruction.

Results with an unsteady and compressible fluid velocity

An unsteady fluid velocity u and the corresponding characteristics are given, for x ∈ [0, 1] by:

u(t, x) = 1 -x 1 + t , X(t; s, x) = 1 + x -1 1 + t -s .
The analytical solution is given by:

m N (t, x) = (1 + t)m 0 N (1 + (x -1)(1 + t)) .
For the simulations, we use N = 9 and the final time is t = 1. The same kind of simulations are done in the unsteady case as with the stationary one. For the ζ and QW kinetic schemes and the regular initial NDF, the errors on the moments are shown in Fig. 8. The numerical order of accuracy is about 1.93 when using the ζ reconstruction and about 1.91 for the moments of order 0 to 5, when using the QW reconstruction. Moreover, in this case, the influence of the CFL is much smaller than with the constant fluid velocity, probably due to the interpolation error of the fluid velocity in the case of the kinetic schemes. All methods then give similar results compared to the constant fluid velocity test case, the ζ reconstruction based schemes being still a little bit more accurate.

Numerical accuracy

The difference between the two kinds of reconstructions is higher when considering the multimodal initial NDF (see Fig. 9 left), this difference being reduced when considering only 4 moments with the QW reconstruction based schemes. In any cases, the accuracy stays much higher than with the constant reconstruction. In Fig. 9(right), it can be seen that the ζ and QW kinetic schemes well capture the mean value of the distribution (corresponding to ζ 1 ), as well as its variance (corresponding to ζ 1 ζ 2 ) with only 100 cells. The corresponding plot for the simplified scheme is very similar and not shown here. Moreover, one can see in Fig. 10 that the second order of accuracy is still obtained in this case, for both the ζ and the QW kinetic schemes. for the exact solution (solid black lines) and for the simulation with the kinetic scheme and a CFL equal to 0.8, using the constant reconstruction (dashed red lines), the ζ reconstruction (blue dashed line) and the QW reconstruction (dotted green line), with 100 points.

-

Computational time

In the case of the multi-modal initial NDF, the computational time of the different schemes with the use of several numbers of moments are given in Table 1, when using a CFL number equal to 0.3 for all methods. This case is chosen since it is representative of physical cases and induces a fair comparison of the methods: some corrections of the slopes are needed for the ζ schemes as well as some abscissa eliminations for a few points for the QW schemes.

It shows that the computational time for the high order schemes is about 2.5 to 3 times higher than for the corresponding first order scheme, for the same discretization, except for the simplified QW scheme and for the ζ kinetic scheme with N ≥ 8 where it is about 4 to 5 times higher. Indeed, when considering the ζ kinetic scheme, Gauss-Legendre quadrature with N 2 + 1 points are used to compute all the integrals, thus still increasing the cost when the number of moments increase. Moreover, for the simplified schemes, the computational times are larger due to the use of the second order SSP Runge-Kutta method, thus needing two reconstructions per time step. However, N the ζ simplified scheme is still competitive, at least for the same CFL, due to the very simple reconstruction. But the CFL cannot be increased beyond 1/3, whereas the global computational time of the other methods can be decreased by an increasing of the CFL number. When considering a global problem, including source terms, the time step can however have to be limited by the coupling characteristic time between the operators [START_REF] Descombes | Adaptive time splitting method for multi-scale evolutionary partial differential equations[END_REF][START_REF] Doisneau | Numerical strategy for unsteady two-way coupling in polydisperse sprays: application to solid rocket motor instabilities[END_REF]. Finally, due to the possible reduction of the number of cells, the cost of the schemes developed here is lower than the one of the first order scheme for the same accuracy. Moreover, the number of degrees of freedom being then reduced, it implies a reduction of the global cost of the complete problem with source terms.

Results with 2D, steady and incompressible fluid velocity

Let us consider a 2D configuration of particles in a vortex. The unsteady fluid velocity is then defined for x ∈ [0, 1/2] and y ∈ [0, 1/2] by: u x (t, x, y) = sin(2πx) cos(2πy), u y (t, x, y) = -cos(2πx) sin(2πy).

Let us introduce the distance r to the point (1/8, 1/8). A population of particles is initially present in the vortex, in a disk defined by r < 1/8. Its distribution is given by f 0 (ξ, 1 -8r), where f 0 corresponds to the multi-modal NDF defined by (G.1), with a weight w 3 (x) now being equal to 1 for x ∈ [2/3, 1]. This population is then transported by the fluid till the time t = 0.8. Due to the incompressibility of the fluid phase, the value of the moments are conserved along the characteristics. A reference solution is then computed by solving reverse characteristics from each position, using a high order ODE solver with time step adaptation. For example, the zeroth order moments obtained with this method is given in Fig. 11(top left). Simulations are done for moments of order 0 to 9, using a 200 × 200 uniform discretization and a CFL number equal to 0.3. To solve the 2D problem, a dimensional splitting is used: a 1D scheme is used alternatively on one-dimensional problems in the x and y directions. To obtain a second order of accuracy with this splitting method but using the same CFL number for each operator, steps of length ∆t are used on each problem, but alternating the order of these steps in alternate time steps [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF].

Since the kinetic and the simplified schemes give very similar results, only the results with the kinetic scheme are presented here. The zeroth order moment is plotted in Fig. 11. One can see the good level of accuracy of the ζ and QW kinetic schemes, especially compared to the first order scheme. For going further in the comparisons, the mean value of the distribution, i.e. m 1 /m 0 and its variance, i.e. ζ 1 ζ 2 are plotted in Fig. 12 as functions of x at y = 0.4. The high numerical diffusion of the first order scheme is still remarkable, whereas the ζ and QW kinetic schemes allow to capture accurately these quantities.

Finally, for the kinetic schemes, the L 1 norm of the errors, divided by the L 1 norm of the corresponding moments, are plotted as a function of the cell width in Fig. 13, using from 100 to 1600 cells per direction. This shows the convergence of the methods and their high accuracy compared to the first order scheme, even if the second order of accuracy is not attained for the considered meshes, the solution having high gradients.

Conclusion

In this paper, we have provided four realizable accurate finite volume schemes for the transport of moments by a given velocity field, in Cartesian mesh context: for the ones, the flux computation is based on a follow-up of the characteristics (kinetic schemes) and for the other, it is based on the value of the moments at the interface (simplified schemes). In any case, a spatial reconstruction of the moments is needed and done by reconstructing variables that only have to be non-negative: either the corresponding ζ k variables or the weights of the corresponding quadrature. Unlike previous developed realizable Eulerian schemes, they are able to deal with moment vectors of all sizes, at least till 10 moments, possibly at the boundary of the moment space, which can occur in practical applications. The verifications have been done on various test cases, 1D and 2D, with steady or unsteady and compressible or incompressible fluid velocities. The developed schemes then showed their high accuracy, compared to the first order scheme, and their second order of accuracy for all moments in the case of the ζ reconstruction and for half the lowest order ones for the QW reconstruction. The kinetic schemes are more accurate than the simplified schemes, especially due to the fact that a higher CFL number can be used. But their generalization to unstructured meshes seems less easy than for the simplified schemes. Moreover, the ζ reconstruction leads to slightly better accuracy than the QW one in most cases. And the accuracy on the QW schemes is much more sensitive to the number of considered moments and can eventually decrease when this number increases. Finally, for a given accuracy, the cost of all the developed schemes is usually lower than the one of the first order scheme. And since they all allow a high reduction of the number of degrees of freedom, the global cost of the complete problem with source terms can be drastically reduced by using these schemes, compared to the first order ones, thus showing the great interest of such kind of schemes. Moreover, we are studying an implementation of a simplified scheme in the open-source computational fluid dynamics toolbox OpenFOAM as part of the Open-QBMM project [START_REF] Openqbmm | An open-source implementation of quadrature-based moment methods[END_REF][START_REF] Passalacqua | An open-source quadrature-based population balance solver for OpenFOAM[END_REF].

In this case of compact support, the QW schemes can be applied without any change. However, the ζ-schemes have to be adapted: the reconstruction has then to be done on the canonical moments to ensure the realizability for the corresponding support. In the case of the reconstruction for the kinetic scheme, the reconstructed ζ k used for the flux computation in [START_REF] Essadki | Adaptive mesh refinement for polydisperse spray simulation[END_REF] are then ζ n k (x) = ξ max p n k (x)(1 -p n k-1 (x)) where p n k (x) = pn k,j + D n k,j (x -x j ) on the j th cell and pn k,j = a k,j + b k,j D n k,j . The coefficients a k,j and b k,j are now given by: a k,j 2 , the slope D n k,j being obtained through a minmod limitation from the p n k,i . The values of the fluxes are then given by [START_REF] Gordon | Error bounds in equilibrium statistical mechanics[END_REF]. Let us remark that the corresponding schemes can be used in the more general case of a support in [0, +∞) by using for ξ max the largest value of the abscissas obtained by a quadrature on each moment set at time t n . Moreover, this value can only decrease with the transport schemes and the moments are then always bounded. with m j = m(x j ) and m j = m (x j ), it is easy to see that except near extreme points of m: θ j = 1 + O(∆x) and θ j+1 -θ j = O(∆x 2 ) and m(x) -m(x) = O(∆x 2 ) for x ∈ (x j-1 2 , x j+ 1 2 ). Without loss of generality, let us assume that 0 ≤ δ < ∆x and let us denote We then have to see if I j-1 2 -I j-1 2 = O(∆x 2 ). First, using Taylor expansion: Then, except near extreme points of m, the property [P3] is verified with this kind of reconstruction. Moreover, let us remark that the use of the minmod limiter allows the reconstruction to have a smaller total variation than m. if a k,j < 0 then k 0 ← max i, Ψ m n 0,j , D n 0,j , . . . , a i-1,j , b i-1,j , D n i-1,j , a i,j , b i,j , 0, . . . , a k-1,j , b k-1,j , 0 > 0 ; n c ← 0; while a k,j < 0 do if n c < 5 then n c ← n c + 1; D n k0,j ← 0.9 D n k0,j ; ζn k0,j ← a k0,j + b k0,j D n k0,j ; Compute a p,j , b p,j and D n p,j from (27,28) and ( 30) for p = k 0 + 1, . . . , k; else n c ← 0; D n k0,j ← 0; ζn k0,j ← a k0,j ; k 0 ← k 0 + 1; ζn k,j ← a k,j + b k,j D n k,j ;

I j-
The reconstruction used for the ζ simplified scheme is detailed in Algorithm 2, especially the corrections done in the case where m * j is not in the interior of the same moment space as m n j , i.e. where N (m * j ) < N (m n j ).

with ξ 1 = 0.02 and where R is the Rosin-Rammler pdf:

1 0Figure 1 :

 11 Figure 1: Constant fluid velocity test case with the regular initial solution; Right: space evolution of the moments of order 0 to 6 (top to bottom lines) for the exact initial solution (solid black lines) and for the simulation with the first order kinetic scheme at time t = 2 (red dashed lines). Left: initial values of the ζ k , k ∈ {1, . . . , 6}.

Figure 2 :

 2 Figure 2: Left: ζ k for k=1,2,3,4 as functions of the spatial location x. Right: initial distribution f 0 (ξ, x) as function of ξ for x = 0.5, x = 0.6 and x = 0.7 (for the Dirac delta functions, the absolute value of the heights of the scaled weights is arbitrary).

Figure 3 :

 3 Figure 3: Constant fluid velocity test case with the regular initial NDF: L 1 norm of the error on the moments of order 0 to 9 (solid red lines, bottom to top) with the ζ kinetic scheme (left) and the QW kinetic scheme (right) for a CFL equal to 0.8. The same line of slope 2 is represented by black dashed lines on the two figures.

Figure 4 :

 4 Figure 4: Constant fluid velocity test case with the regular initial solution: L 1 norm of the error on the moments of order 0 to N (solid red lines, bottom to top) with the QW kinetic scheme for a CFL equal to 0.8 and for N = 5 (left) and N = 4 (right). The same line of slope 2 is represented by black dashed lines on the two figures.

Figure 5 :

 5 Figure 5: Constant fluid velocity test case with the regular initial NDF: L 1 norm of the error on the moments of order 0 (left) and 1 (right) with the kinetic (solid lines) of the simplified (dashed lines) schemes for a CFL equal to 0.3, using the constant (black +), the ζ (red) or the QW (blue ) reconstruction.

Figure 6 :

 6 Figure6: Constant fluid velocity test case with the oscillating initial ζ k : space evolution of the moments of order 0 to 4 (left, top to bottom lines) and of order 5 to 8 (right, top to bottom lines for x < 0.7) at time t = 5 for the exact solution (solid black lines) and for the simulation with the constant reconstruction (top), with the ζ reconstruction (middle) and the QW reconstruction (bottom) for the kinetic scheme (red dashed line or green dots for QW with 4 moments) and the simplified scheme (blue dotted lines) with 100 points and a CFL equal to 0.3.

Figure 7 :

 7 Figure7: Constant fluid velocity test case with the multi-modal initial NDF: space evolution of the moments of order 0 to 4 (left, top to bottom lines) and of order 5 to 8 (right, top to bottom lines) at time t = 5 for the exact solution (solid black lines) and for the simulation with the ζ reconstruction (top) and the QW reconstruction (bottom) for the kinetic scheme (red dashed line or green dots for QW with 4 moments) and the simplified scheme (blue dotted lines) with 100 points and a CFL equal to 0.3.

Figure 8 :

 8 Figure 8: Unsteady fluid velocity test case with the regular initial NDF: L 1 norm of the error on the moments of order 0 to 9 (solid red lines, bottom to top) with the ζ kinetic scheme (left) and the QW kinetic scheme (right) for a CFL equal to 0.8. The same line of slope 2 is represented by black dashed lines on the two figures.

Figure 9 :

 9 Figure9: Unsteady fluid velocity test case with the multi-modal initial NDF. Left: L 1 norm of the error on the moments of order 0 with the kinetic (solid lines) of the simplified (dashed lines) schemes for a CFL equal to 0.3, using the constant (black +), the ζ (red) or the QW (blue when N = 9 and green • when N = 3) reconstruction. Right: space evolution of the mean value ζ 1 (top curves) and the variance ζ 1 ζ 2 (bottom curves) of the NDF at time t = 1 for the exact solution (solid black lines) and for the simulation with the kinetic scheme and a CFL equal to 0.8, using the constant reconstruction (dashed red lines), the ζ reconstruction (blue dashed line) and the QW reconstruction (dotted green line), with 100 points.

Figure 10 :

 10 Figure 10: Unsteady fluid velocity test case with the multi-modal initial NDF: L 1 norm of the error on the moments of order 0 to 9 (solid red lines, bottom to top) with the ζ kinetic scheme (left) and the QW kinetic scheme (right) and with the first order scheme (both figures, dotted green lines) for a CFL equal to 0.8. The same lines of slope 1 and 2 are represented by black dashed lines on the two figures.

Figure 11 :

 11 Figure 11: 2D test case: moment m 0 for the reference solution (top left), for the simulations on a 200 × 200 uniform mesh with a CFL number equal to 0.3, using the first order kinetic scheme (top right), the ζ kinetic scheme (bottom left) and the QW kinetic scheme (bottom right).

Figure 12 :

 12 Figure 12: 2D test case: mean value (left) and variance (right) of the distribution as a function of x at y = 0.4 for the reference solution (black solid line), for the simulations on a 200 × 200 uniform mesh and a CFL number equal to 0.3, with the first order kinetic scheme (green dots), with the ζ kinetic scheme (red dashed line) and with the QW kinetic scheme (blue dotted line)

Figure 13 :

 13 Figure13: 2D test case: L 1 norm of the error on the moments of order 0 to 9 (solid red lines, bottom to top) with the ζ kinetic scheme (left) and the QW kinetic scheme (right) and with the first order scheme (both figures, dashed green lines) for a CFL equal to 0.8.

  n i (x))dx.The same kind of slope limiters and corrections are used as fo the ζ kinetic scheme. For the simplified scheme, the reconstructedζ k are ζ ± k,j = ξ max p ± k,j (1 -p ± k-1,j ) with p ± k,j = p n k,j ± D n k,j

Figure A. 14 :

 14 Figure A.14: Constant fluid velocity test case with the multi-modal initial NDF: space evolution of the moments of order 0 to 4 (left, top to bottom lines) and of order 5 to 8 (right, top to bottom lines) at time t = 5 for the exact solution (solid black lines) and for the simulation with the canonical moments reconstruction for the kinetic scheme (red dashed lines) and the simplified scheme (blue dotted lines) with 100 points and a CFL equal to 0.3.These schemes are used in the constant fluid velocity test case with the multi-modal initial NDF. The 10 moments obtained by the simulation with 100 cells are plotted in Fig. A.14. It shows a similar accuracy compared to the ζ schemes.

Algorithm 1 :

 1 Reconstruction in the cell j for the ζ kinetic scheme Data: m n 0,j , m n 0,j-1 , m n 0,j+1 , ζ n k,j , ζ n k,j-1 , ζ n k,j+1 k∈{1,...,N } Result: Reconstruction parameters: ζn k,j k∈{1,...,N } and D n k,j k∈{0,...,N } Compute D 0,j from (29); for k ← 1 to N do Compute a k,j , b k,j from (27,28); Compute D k,j from (30) ;

.Figure G. 16 :

 16 Figure G.16: Left: weights w 1 , w 2 and w 3 as functions of the spatial location x. Right: functions λ(x) (black solid line) and k(x) (red dashed line).

Table 1 :

 1 Unsteady fluid velocity test case with the multi-modal initial NDF: normalized computational time for a CFL number equal to 0.3.

	4	5	6	7	8	9

r denotes here the largest integer less than or equal to the real number r

Acknowledgments

This work was supported by a grant from French National Research Agency: ASMAPE project (ANR-13-TDMO-02 ASMAPE, PI IFPen: O. Colin; PI for EM2C Lab.: F. Laurent).

The authors would like also to thank Professor R. O. Fox for several helpful discussions and for the precious feedback he provided about this paper.

, with d = 0, 1; n ≥ 0. Similarly to the case with support in [0, +∞), one has the following characterization of the moment space [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF]: the vector m N = (m 1 , . . . , m N ) t is realizable (i.e. in the moment space) if and only if H k ≥ 0 and H k ≥ 0, k ∈ {0, 1, . . . , N }.

Moreover, it is strictly realizable (i.e. in the interior of the moment space) if and only if these Hankel determinants are positive.

From the ζ k corresponding to the moments m k (then the ζ k /ξ max correspond to the mk ), one can defined the canonical moments

) for µ ∈ P(m k-1 ), (i.e. p k is 0 or 1), the measure µ is a sum of weighted Dirac distributions and m k belongs to the boundary of the moment space.

Appendix B. Examples of realizability constraints for moments on [0, +∞)

The first constraints [START_REF] Sporleder | On the population balance equation[END_REF] for the strict realizability of a moment vector (m 0 , m 1 , . . . ) t can be written, making appear the values of m - k (m k-1 ): m 0 > 0, m 1 > 0,

Appendix C. Moments as functions of the ζ k

The moments can be written from the

as shown in Corollary 2.3. The first polynomial functions P k are given by:

Appendix D. Property [P3] for classical MUSCL reconstruction

Let us consider just one moment, denoted m, without any indice here. Let us denote m(x) the exact solution of the transport equation at a given time and m j = 1 ∆x x j+ 1 2

x j-1 2 m(x)dx its averaged value in the cell j, for any j. Let us then consider the affine reconstruction m(x) from the m j , corresponding to a MUSCL scheme:

). The slope D j is then an approximation of the derivative of m, using a limiter [START_REF] Leveque | Numerical methods for conservation laws[END_REF]:

where the function Φ defines the flux limiter, such that Φ(1) = 1 and Φ is bounded and Lipschitz continuous at θ = 1 . Using Taylor expansions of m around x j , for example

The reconstruction used for the ζ kinetic scheme is detailed in Algorithm 1, especially the corrections done in the case a k,j < 0. More precisely, this algorithm computes the parameters of the reconstruction in the cell j ( ζn k,j for k ∈ {1, . . . , N } and the slopes D n k,j for k ∈ {0, . . . , N }), providing the zeroth order moments and the values of the ζ k corresponding to m n j , m n j-1 and m n j+1 . For that, the formula [START_REF] Vikas | Realizable high-order finite-volume schemes for quadraturebased moment methods applied to diffusion population balance equations[END_REF] leading to the computation of a k,j from the reconstruction of m 0 and the ζ i for i < k is denoted: Compute the moments m ± j from m ± 0,j and (ζ ± k,j ) k∈{1,...,N } (reverse algorithm);

The reconstruction used for the QW schemes is detailed in Algorithm 3. For that, let us denote V (ξ α ) α∈{1,...,n} the Vandermonde matrix corresponding to the coefficients (ξ 1 , . . . , ξ n ), i.e. the matrix of coefficients ξ i-1 j on the row i ∈ {1, . . . , n} and column j ∈ {1, . . . , n}.

Algorithm 3: Reconstruction in the cell j for the QW schemes Data: m n j , m n j-1 , m n j-1

Result: w α,j , ξ α,j , D n α,j α∈{1,...,

Compute (w α,j , ξ α,j ) α∈{1,...,p} from (m n k,j ) k∈{0,...,N (m n j )-1} (Chebychev algorithm); C = {1, . . . , p};

while any w ± α,j α∈C < 0 do p ← p -1; Determine α 0 such that w ± α0,j = min{w ± α,j , α ∈ C};

Appendix F. Reconstruction at the level of the NDF Another way to reconstruct the moments m n (x) is to consider a spatial reconstruction at the NDF level [START_REF] Garrett | A Comparison of Moment Closures for Linear Kinetic Transport Equations: The Line Source Benchmark[END_REF][START_REF] Schneider | A realizability-preserving high-order kinetic scheme using WENO reconstruction for entropy-based moment closures of linear kinetic equations in slab geometry[END_REF]. Let us then reconstruct a NDF f n j from the moment set (m n k,j ) k∈{0,...,N } , using EQMOM or entropy maximization and let us define

]. The slope is defined thanks to a minmod limiter:

and the reconstructed moments for the kinetic scheme are then m n k (x) = +∞ 0 ξ k f n (x, ξ)dξ . Practically, a quadrature on the ξ variable is used to compute the fluxes. But here, to ensure the [P1] property, this quadrature has to correspond to the measure f j (ξ)dξ, i.e.

where we define w α,j = wα,j /f n j (ξ α,j ), since f n j (ξ α,j ) cannot be zero. In the case of EQMOM reconstruction, the weights wα,j and abscissas ξ α,j can correspond to the secondary quadrature [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF][START_REF] Nguyen | Solution of population balance equations in applications with fine particles: Mathematical modeling and numerical schemes[END_REF]. This then leads to:

The corresponding kinetic scheme is then realizable and the flux are written: Figure F.15: Constant fluid velocity test case with a multi-modal initial NDF: space evolution of the moments of order 0 to 6 (top to bottom) at time t = 1 for the exact solution (solid black lines) and for the simulation with the kinetic scheme using a reconstruction of the NDF with the gamma-EQMOM (left, red dashed lines) and the log-normal-EQMOM (right, red dashed lines) with 100 cells and a CFL equal to 0.8.

However, there are three drawbacks for this kind of scheme. First, it depends on the choice of the NDF reconstruction, since several reconstructions are possible, and it also adds the cost of the reconstruction itself, which is not negligible. However, such kind of reconstruction is usually needed for the other operators present in the complete physical problem. The second and more problematic drawback comes from the impossibility to deal with the boundary of the moment space with this kind of method. Finally, the maximum principle on the moments may not be preserved. Indeed, for example, f n (x, ξ) is smaller than max{f n j (ξ), f n j-1 (ξ), f n j+1 (ξ)} for any value of ξ. But the moments of the distribution ξ → max{f n j (ξ), f n j-1 (ξ), f n j+1 (ξ)} has no reason to be limited by max{m n k,j , m n k,j+1 , m n k,j-1 }. This can be illustrated on a multi-modal but regular initial NDF: f 0 (ξ, x) = 9x 2 (2 -3x) 2 1 [0,2/3] (x)R(ξ, λ 1 , k 1 ) + 9(3x -1) 2 (1 -x) 2 1 [1/3,1] (x)R(ξ, λ(x), k(x)), where R is the Rosin-Rammler pdf and with λ 1 = 0.03, k 1 = 2. The parameters λ(x) and k(x) are plotted in Fig. G.16(right). Here, both the gamma and the log-normal EQMOM reconstructions [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF][START_REF] Madadi-Kandjani | An extended quadrature-based moment method with log-normal kernel density functions[END_REF][START_REF] Nguyen | Solution of population balance equations in applications with fine particles: Mathematical modeling and numerical schemes[END_REF] are used to computed fluxes defined by (F.1). The result of the simulations for 5 moments with the corresponding kinetic scheme with 100 cells and a CFL number equal to 0.8 are plotted in Fig. F.15. The maximal principle is not respected for the zeroth order moment and the results depend on the reconstruction, with a slightly better behavior here when using the gamma-EQMOM.

Appendix G. Multi-modal initial NDF

The multi-modal NDF is defined by: f 0 (ξ, x) = w 1 (x)δ ξ1 (ξ) + w 2 (x)δ 2ξ1 (ξ) + w 3 (x)R(x, λ(x), k(x)), (G.1)