T T Nguyen

F Laurent

R O Fox
email: rofox@iastate.edu

M Massot
email: marc.massot@centralesupelec.fr

Solution of population balance equations in applications with fine particles: mathematical modeling and numerical schemes

Keywords: aerosol, population balance equation, quadrature-based moment method, sectional method, hybrid method

The accurate description and robust simulation, at relatively low cost, of global quantities (e.g. number density or volume fraction) as well as the size distribution of a population of fine particles in a carrier fluid is still a major challenge for many applications. For this purpose, two types of methods are investigated for solving the population balance equation with aggregation, continuous particle size change (growth and size reduction), and nucleation: the extended quadrature method of moments (EQMOM) based on the work of Yuan et al. (J. Aerosol Sci., 51:1-23, 2012) and a hybrid method (TSM) between the sectional and moment methods, considering two moments per section based on the work of Laurent et al. (Commun. Comput. Phys., accepted, 2016). For both methods, the closure employs a continuous reconstruction of the number density function of the particles from its moments, thus allowing evaluation of all the unclosed terms in the moment equations, including the negative flux due to the disappearance of particles. Here, new robust and efficient algorithms are developed for this reconstruction step and two kinds of reconstruction are tested for each method. Moreover, robust and accurate numerical methods are developed, ensuring the realizability of the moments. The robustness is ensured with efficient and tractable algorithms despite the numerous couplings and various algebraic constraints thanks to a tailored overall strategy. EQMOM and TSM are compared to a sectional method for various simple but relevant test cases, showing their ability to describe accurately the fine-particle population with a much lower number of variables. These results demonstrate the efficiency of the modeling and numerical choices, and their potential for the simulation of real-world applications.

Introduction

The evolution of a population of fine, that is non-inertial, particles in a carrier fluid can be described by a population balance equation (PBE) [START_REF] Marchisio | Quadrature method of moments for population-balance equations[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF][START_REF] Piton | Simulation of fine particle formation by precipitation using computational fluid dynamics[END_REF][START_REF] Raman | Modeling of fine-particle formation in turbulent flames[END_REF][START_REF] Ramkrishna | Population balance modeling: Current status and future prospects[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][START_REF] Sung | Large-eddy simulation modeling of turbulent flame synthesis of titania nanoparticles using a bivariate particle description[END_REF][START_REF] Wang | Application of in situ adaptive tabulation to CFD simulation of nanoparticle formation by reactive precipitation[END_REF][START_REF] Zucca | Implementation of the population balance equation in CFD codes for modelling soot formation in turbulent flames[END_REF]. There are many potential applications such as soot modeling, aerosol technology, nanoparticle synthesis, microbubbles, reactive precipitation, and coal combustion (see [START_REF] Ramkrishna | Population balance modeling: Current status and future prospects[END_REF] and references therein). The PBE is a transport equation for the number density function (NDF) of the particles. The NDF depends on time, spatial location and the internal coordinates, which can include, for example, volume, surface area or chemical composition. The mathematical form of a typical PBE includes spatial transport (e.g. convection and diffusion), derivative source terms for continuous particle size change (e.g. oxidation/dissolution and surface growth), integral terms (e.g. aggregation and breakage), and Dirac-delta-function source terms describing the formation of particles (e.g. nucleation). Moreover, it is usually not only important to predict the evolution of global quantities of the particle population, but also to have some information on the NDF. For example, when considering soot, the total produced mass or volume fraction, as well as the size-dependent NDF represent essential elements in present and future emission regulations. But since the evolution of the NDF is usually coupled with the resolution of the Navier-Stokes equation for the carrier fluid [START_REF] Zucca | Implementation of the population balance equation in CFD codes for modelling soot formation in turbulent flames[END_REF], the cost of its resolution has to be reasonable so that the global simulation will be affordable. We therefore seek a robust method able to describe accurately some global quantities of the particle population, but also able to give a good idea of the shape of the NDF, at a reasonable cost.

In this work, only size is considered as the internal variable. Moreover, even if the geometry of fine particles can be complex, we assume that only one variable, e.g. volume v, is needed to describe it, eventually taking into account the more complex shape of large particles compared to smaller ones through a fractal dimension depending on size. This work thus represents a first step before considering more complex models that add another internal coordinate variable. Different methods are available in the literature to solve the PBE. The Monte-Carlo method [START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF] is usually too costly to be coupled with a flow solver, especially when considering particle interactions like aggregation. We therefore focus on deterministic methods. With one internal variable, deterministic methods can be based on a discretization along the size variable. Equations are written for the total number density or the total mass density of the particles inside each interval of the size-discretization. These intervals are called sections in what follows, in reference to the sectional methods, which fall in this category. A large literature is devoted to this type of method, especially for the resolution of the aggregation and/or breakage PBE (see e.g. [START_REF] Kumar | On the solution of population balance equations by discretization-I. A fixed pivot technique[END_REF][START_REF] Vanni | Approximate population balance equations for aggregation-breakage processes[END_REF][START_REF] Kumar | Convergence analysis of a finite volume scheme for solving non-linear aggregation-breakage population balance equations[END_REF] and references therein). Among them are the fixed-pivot [START_REF] Kumar | On the solution of population balance equations by discretization-I. A fixed pivot technique[END_REF] and the cell-average [START_REF] Kumar | Improved accuracy and convergence of discretized population balance for aggregation: The cell average technique[END_REF] techniques, which consider that the particle population of one section is represented by only one size (pivot size), the new particles after collision or breakup being distributed over the sections in such a way that the discrete equations are consistent with the global number and mass (they are said to be "moment preserving"), or also with higher order moments [START_REF] Alopaeus | Solution of population balances with breakage and agglomeration by high-order moment-conserving method of classes[END_REF]. These methods have been generalized to growth, nucleation and aggregation [START_REF] Kumar | On the solution of population balance equations by discretization-III. Nucleation, growth and aggregation of particles[END_REF][START_REF] Alopaeus | Solution of population balances with growth and nucleation by high order moment-conserving method of classes[END_REF]. Some other methods are based on a "conservative form" of the PBE for the mass density function [START_REF] Filbet | Numerical simulation of the Smoluchowski coagulation equation[END_REF][START_REF] Bourgade | Convergence of a finite volume scheme for coagulation-fragmentation equations[END_REF]: a conservative finite-volume method developed for aggregation and breakup [START_REF] Filbet | Numerical simulation of the Smoluchowski coagulation equation[END_REF][START_REF] Bourgade | Convergence of a finite volume scheme for coagulation-fragmentation equations[END_REF] and extended to growth and nucleation [START_REF] Qamar | On the solution of population balances for nucleation, growth, aggregation and breakage processes[END_REF],

or some moment-preserving methods [START_REF] Kumar | Moment preserving finite volume schemes for solving population balance equations incorporating aggregation, breakage, growth and source terms[END_REF]. These methods have been shown to be convergent when considering only aggregation or breakup, with first-order accuracy for the finite-volume methods [START_REF] Filbet | Numerical simulation of the Smoluchowski coagulation equation[END_REF][START_REF] Bourgade | Convergence of a finite volume scheme for coagulation-fragmentation equations[END_REF] and second-order accuracy for the fixed-pivot and cell-average techniques [START_REF] Kumar | Convergence analysis of sectional methods for solving breakage population balance equations. I. The fixed pivot technique[END_REF][START_REF] Kumar | Convergence analysis of sectional methods for solving breakage population balance equations. II. The cell average technique[END_REF][START_REF] Giri | Convergence analysis of sectional methods for solving aggregation population balance equations: the fixed pivot technique[END_REF]. But a large number of sections is always used and the ability of these methods to describe adequately the NDF with a small number of sections has not been fully explored, especially for the complete problem with nucleation and growth. Moreover, to our knowledge, such methods have never been reported for cases where the particle size is decreasing through a continuous process.

A different kind of method, the only one that will be called sectional here (even if some of the previous ones are also called sectional in the literature), is based on a closure through a continuous reconstruction of the NDF inside each section. This reconstruction can be constant [START_REF] Gelbard | Sectional representations for simulating aerosol dynamics[END_REF][START_REF] Laurent | Multi-fluid modelling of laminar polydisperse spray flames: origin, assumptions and comparison of sectional and sampling methods[END_REF][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF] or affine [START_REF] Netzell | Calculating the soot particle size distribution function in turbulent diffusion flames using a sectional method[END_REF]. When considering sprays, sectional methods are also called "Eulerian multi-fluid methods" [START_REF] Laurent | Multi-fluid modelling of laminar polydisperse spray flames: origin, assumptions and comparison of sectional and sampling methods[END_REF] or "one-size moment" method (OSM), and they are developed after reduction of the internal variables to only size thanks to velocity moments and a mono-kinetic closure. The corresponding model is a finite-volume method. It was shown to be first-order accurate in the pureevaporation case [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF] and exhibited first-order numerical accuracy for the investigated cases, taking into account the collisions (coalescence) [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF]. Moreover, in this approach, an affine reconstruction of MUSCL type was tested in the pure-evaporation case. However, if its order of accuracy is higher, the effective accuracy is not much improved compared to the first-order method, except with a large number of sections. Then, as with other discretized methods, sectional methods lead to an accurate prediction of the NDF with a large enough number of sections. However, for many applications, physical transport must also be considered, as well as coupling with the carrier fluid. This can be done thanks to the use of an operator-splitting method (see e.g. [START_REF] Doisneau | Numerical strategy for unsteady two-way coupling in polydisperse sprays: application to solid rocket motor instabilities[END_REF]), but if a large number of sections has to be considered, the computational cost can be prohibitive, since one has to transport at least one variable per section.

In contrast, moment methods do not use a direct resolution of the NDF, but rather transport a finite set of its moments, usually the first few integer order ones. Since they are the moments of a non-negative NDF (or, more rigorously, a positive measure), this moment set belongs to a space strictly included in R N + , where N is the number of moments [START_REF] Shohat | The Problem of Moments, 4th Edition[END_REF][START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF]34]. This space is called the moment space. The NDF cannot be recovered from this finite moment set: there is an infinite number of possibilities in non-degenerate cases, i.e. when the moment set is in the interior of moment space, whereas a unique sum of weighted Dirac delta functions is possible for the degenerate cases, i.e. for the boundary of moment space. One can remark that the degenerate case can appear in problems of interest due to nucleation of fine particles just as they begin to aggregate. Most importantly, moment methods give access to some important properties of the NDF.

For moment methods, two major issues arise. The first is closure of the moment equations due to the nonlinear source terms in the PBE. This includes the negative flux due to the disappearance of particles when continuous size reduction is considered (e.g. oxidation or evaporation), which requires a point-wise evaluation of the NDF [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]. Two kinds of closures are used in the literature: (i) a functional dependence of the unclosed terms (usually expressed through some fractional moments) is provided using the moment set; (ii) a NDF, or its corresponding measure, is reconstructed from the moment set, allowing evaluation of all the unclosed terms. In the first category, one finds the interpolative closure (MOMIC) [START_REF] Frenklach | Aerosol dynamics modeling using the method of moments[END_REF][START_REF] Frenklach | Method of moments with interpolative closure[END_REF] widely used in the soot community and extended to the bivariate case [START_REF] Mueller | A joint volume-surface model of soot aggregation with the method of moments[END_REF]. MOMIC is based on an interpolation along the order of the moments. However, this kind of method does not allow one to deal with the disappearance fluxes, except for the hybrid method of moments (HMOM) [START_REF] Mueller | Modeling the oxidation-induced fragmentation of soot aggregates in laminar flames[END_REF], which is a combination of MOMIC and DQMOM described below and was developed for a bivariate case. Moreover, these methods do not guarantee that the closure can correspond to any NDF. This is why a second way to close the moment equations has been developed using quadrature-based moment methods (QBMM) [START_REF] Marchisio | Quadrature method of moments for aggregationbreakage processes[END_REF][START_REF] Marchisio | Quadrature method of moments for population-balance equations[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF]. Because the internal variable (size) is assumed to live in the space [0, ∞), the problem of NDF reconstruction from a finite set of moments is known as the truncated Stieljes moment problem [START_REF] Shohat | The Problem of Moments, 4th Edition[END_REF]34] 1 . Among these reconstructions, a sum of weighted Dirac delta functions can be used, leading to the widely employed quadrature method of moments (QMOM) [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF]. This reconstruction for an even set of integer moments (m 0 , m 1 , . . . , m 2N -1 ) is the lower principal representation, i.e. the corresponding moment m 2N is minimal. In its variant, the direct quadrature method of moments (DQMOM) [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF], the equations are directly written for the weights and abscissas of the reconstruction, leading however to some shortcomings related to the conservation of moments. Although these methods have the advantage of being applicable to multivariate cases (directly for DQMOM or using the conditional quadrature method of moments (CQMOM) [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF]), they still cannot deal with disappearance fluxes [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: Direct quadrature method of moments and multi-fluid method[END_REF]. A continuous NDF reconstruction must be considered for this purpose.

One such reconstruction is entropy maximization, which is well defined for the truncated Hausdorff moment problem [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], and an algorithm is available to reproduce any moment set with reasonable accuracy [START_REF] Vié | Size-velocity correlations in hybrid high order moment/multifluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF]. Entropy maximization has been used for sprays [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[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][START_REF] Essadki | Adaptive mesh refinement for polydisperse spray simulation[END_REF][START_REF] Essadki | High order moment model for polydisperse evaporating sprays towards interfacial geometry description[END_REF]. However, when considering the truncated Stieljes moment problem, it fails to reproduce some moment sets [START_REF] Tagliani | Maximum entropy solutions and moment problem in unbounded domains[END_REF]. Some other types of reconstructions have been proposed (see [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF] and references therein), but either with a too low number of moments in such a way that multi-nodal distributions, which often appear in the problems of interest here, cannot be described, or with almost always some negative values for the reconstructed NDF, thus leading to potential instability issues for their numerical resolution. For the truncated Hausdorff moment problem, a nonnegative reconstruction was developed using a superposition of kernel density functions (KDF) (kernel density element method, KDEM) [START_REF] Athanassoulis | The truncated Hausdorff moment problem solved by using kernel density functions[END_REF]. However, KDEM only guarantees that some of the low-order moments are exactly preserved [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF]. First for the truncated Hamburger moment problem [START_REF] Chalons | A multi-Gaussian quadrature method of moments for gasparticle flows in a LES framework[END_REF], and then for the truncated Hausdorff and Stieljes moment problems [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], a nonnegative reconstruction was developed using QBMM, allowing to exactly preserve all the moments except sometimes the last one. The reconstructed NDF is then a sum of nonnegative weighted KDF, able to converge to the Dirac delta function, even if this transition is not yet numerically effective in the proposed algorithms.

The QBMM corresponding to this closure is called the extended quadrature method of moments (EQMOM) and was, for example, applied to soots in [START_REF] Salenbauch | Modeling soot formation in premixed flames using an extended conditional quadrature method of moments[END_REF].

The second major issue associated with moment methods is realizability: the moments must remain in moment space, which is a convex space [START_REF] Shohat | The Problem of Moments, 4th Edition[END_REF][START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF]34]. This issue is not always considered, especially when the first type of closure is used, thus leading to unphysical results (e.g. invalid moment sets). Indeed, even if the closure itself ensures the realizability at the continuous level, the classical schemes for high-order transport in physical space can lead to an invalid moment set [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 transport in phase space, especially when considering continuous particle size reduction [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF] and/or the transition between a Dirac distribution (due to the nucleation) and a smooth distribution (due to the aggregation/coagulation). 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 not sufficient condition for realizability [START_REF] Shohat | The Problem of Moments, 4th Edition[END_REF] in order to obtain a valid moment set. The cost of the method then increases and the correction spoils the overall accuracy. This is why, in this paper, the developed schemes directly preserve the realizability of the moment set [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Essadki | High order moment model for polydisperse evaporating sprays towards interfacial geometry description[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF].

A third type of method, which is a hybrid method between the sectional and moment methods, has also been developed. It consists in using more than one moment per section. The idea is to have a better representation of the NDF in the sections, allowing for a smaller number of sections compared to the sectional methods. The moving pivot technique [START_REF] Kumar | On the solution of population balance equations by discretization-II. A moving pivot technique[END_REF] can be seen as belonging to this category, considering two moments per section (but equations on the zeroth order moment and the pivot, which is a mean size), as well as its generalizations to any number of moments per section: the sectional QMOM (SQMOM) [START_REF] Attarakih | Solution of the population balance equation using the sectional quadrature method of moments (sqmom)[END_REF] and the sectional DQMOM (SDQMOM) [START_REF] Bruns | Development of a hybrid sectional quadrature-based moment method for solving population balance equations[END_REF].

The NDF is then represented by one or a sum of weighted Dirac delta functions. However, these methods were essentially developed for aggregation and breakup. Their adaptation to growth could be done through a method of characteristics, thus inducing a movement of the section bounds and the eventual need of a supplementary section if the nucleation is also considered [START_REF] Kumar | On the solution of population balance equations by discretization-III. Nucleation, growth and aggregation of particles[END_REF]. But, this can usually not be done for the consideration of size decreasing without the disappearance of some pivots or some sections, thus inducing some difficulties, especially with the physical transport. Otherwise, for fixed sections, these kinds of hybrid methods suffer from the difficulty of quadrature moment methods to deal with the fluxes between sections due to growth or size reduction phenomena, as for the disappearance flux at zero size in case of size decreasing [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]. More recently, other kinds of hybrid methods has been developed in the context of sprays [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] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[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], using a continuous NDF reconstruction on each section, thus allowing to evaluate the fluxes between the sections. Moreover, if the realizability issue has also to be considered, then contrary to pure moment methods, its complexity is often low since only a few moments are considered (typically two or eventually four). Among these hybrid methods, in spray modeling context, the two-size moment (TSM) method has been shown to be very accurate for evaporation (i.e. size reduction), as well as for coalescence (i.e. collisions) thanks to a second-order accurate reconstruction of the NDF from the moments in the section [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF], leading also to a good representation of the NDF with a small number of sections. Thus, TSM is a useful method for the fine-particle applications of interest here.

In this paper, we focus on both EQMOM and TSM in the context where the following physical phenomena are considered (i.e. all types of processes except spatial transport): nucleation, aggregation, continuous growth and size reduction. New robust reconstruction algorithms are provided for each method, able to deal with the boundary of the moment space in the case of EQMOM. Moreover, since operator time-splitting techniques will be used for the complete problem, realizable schemes are developed for each operator separately, especially a new one for continuous particlesize change. They are also tested on simple but relevant test cases, isolating the most challenging aspects when using moment methods. The remainder of the paper is organized as follows. Section 2 is dedicated to the description of EQMOM and TSM: the closures are given, as well as efficient algorithms to compute them. Then, realizable numerical schemes are provided and the methods are compared to a sectional method, using a constant reconstruction in each section, the OSM method, when considering particle size reduction and growth (Sec. 3), aggregation (Sec. 4), and a combination of nucleation, aggregation and size reduction (Sec. 5). Conclusions are drawn in Sec. 6.

Mathematical Model

The hybrid and quadrature-based moment methods that will be used in this work are derived from the spatially homogeneous and mono-variate PBE, considering nucleation, aggregation and continuous particle size change. For the hybrid method, TSM is used, with a reconstruction in the volume or radius variable. For the QBMM, EQMOM is used with a gamma or a log-normal KDF. The PBE is first recalled and the equations for the moments of the NDF on an arbitrary interval are given. To close the equations, a non-negative NDF is reconstructed and the details of this reconstruction are provided for each method. Moreover, the moment space is described in each case.

Population balance equation (PBE)

Let us consider the volume v as the only internal variable. In the case of spatial homogeneity, the NDF f (t, v) is only a function of time t and v and the PBE reads

∂ t f (t, v) = S nuc + S agg + S gro + S red (1)
with t ≥ 0 and v ≥ 0. The source terms corresponding to nucleation S nuc , aggregation S agg , surface growth S gro , and continuous size reduction S red are given, respectively, by

S nuc = j(t) δ(v -V nuc ) (2a)
S agg = 1 2 v 0 β(t, v -v , v ) f (t, v -v ) f (t, v ) dv -f (t, v) ∞ 0 β(t, v, v ) f (t, v ) dv (2b) S gro = ∂ v [R gro (t, v) f (t, v)] S red = ∂ v [R red (t, v) f (t, v)] (2c)
where the nuclei volume V nuc is constant and j(t) is the nucleation rate [START_REF] Blanquart | Analyzing the effects of temperature on soot formation with a joint volume-surface-hydrogen model[END_REF][START_REF] Mueller | A joint volume-surface model of soot aggregation with the method of moments[END_REF][START_REF] Mueller | Modeling the oxidation-induced fragmentation of soot aggregates in laminar flames[END_REF]. Moreover,

β(t, v, v
) is the aggregation kernel, R gro (t, v) ≥ 0 the surface growth rate, and R red (t, v) ≤ 0 the rate of size reduction. From the literature there are several types of aggregation kernels such as sum, product or Brownian [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF]. Only the sum and Brownian kernels will be considered in what follows, the first for verification purposes since some analytical solutions are available, and the second for a more physical dependence on the size. Time-dependent kernels will not be considered.

The drift terms (2c) depend on the particle geometry as well as on the size change process. They represent the rate dv dt of change of the particle volume. In the case of particle size reduction, R red (t, v) is usually proportional to the particle surface. For a spherical particle, this leads to

R red (t, v) = -c red (t) v 2/3 1 R+ (v) ≤ 0 (3)
where c red (t) depends only on time. While the small particles are nearly spherical, this is usually not the case of the larger ones. A fractal dimension depending on the volume can be used in order to express the surface area as a function of the volume. However, this is not done here: the methods will be evaluated with simple models in the framework of this paper. For surface growth, we consider diffusion-controlled growth [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF], meaning that the radius r of the particle increases proportionally to 1/r. For spherical particles, this leads to the following volume change rate (see Appendix A):

R gro (t, v) = c gro (t) v 1/3 ≥ 0 (4)
where c gro (t) is independent of the volume. In what follows, the variable v 1/3 , which is proportional to the radius when considering spherical particles, will be denoted r and used as the variable of interest in our example results.

A typical solution of the PBE in fine-particle applications starts with a Dirac delta function, due to nucleation [START_REF] Raman | Modeling of fine-particle formation in turbulent flames[END_REF]. Then, the aggregation causes particles of larger sizes to appear, whereas the drift terms make the sizes of all particles evolve continuously. So, from monodisperse, the distribution becomes polydisperse, with eventually a short period of time where there are only a few sizes. One then has to deal with all these cases with our approximate models based on QBMM and hybrid methods.

Moment equations and related issues

Let us consider the interval (V min , V max ), which will be the space [0, ∞) for QBMM or a section for the hybrid method. Consider the moment of order k of the NDF on this interval: m k (t) = Vmax Vmin v k f (t, v) dv, with k ∈ {0, 1, . . . , N }. Multiplying the PBE (1) by v k and integrating over the support, we obtain the following ordinary differential equations (ODE) for m k :

d t m k = S nuc , v k + S agg , v k + S gro , v k + S red , v k (5)
where

S nuc , v k = V k nuc j(t) 1 (Vmin,Vmax) (V nuc ) (6a) S agg , v k = 1 2 Ω (v + v ) k β(t, v, v )f (t, v) f (t, v )dvdv - Vmax Vmin v k f (t, v) ∞ 0 β(t, v, v ) f (t, v )dv dv (6b) S gro , v k = V k max R gro (t, V max )f (t, V max )-V k min R gro (t, V min )f (t, V min )-k Vmax Vmin v k-1 R gro (t, v)f (t, v)dv (6c) S red , v k = V k max R red (t, V max )f (t, V max )-V k min R red (t, V min )f (t, V min )-k Vmax Vmin v k-1 R red (t, v)f (t, v)dv (6d)
and

Ω = {(v, v ) > 0 V min ≤ v + v ≤ V max }. This set of ODEs for k ∈ {0, 1, . . . , N } is unclosed
since the NDF is unknown. Here, in order to be sure that the source terms are physical, the NDF will be reconstructed from its moments. This reconstruction will have to be well defined for all the physically possible cases, including monodisperse or discrete polydisperse cases, and will have to give a representative value of the NDF at each bound of the interval in cases where this bound is not infinite. Finally, the set of moments lives in a convex space, which is described in the next subsection, and the numerical scheme will have to guarantee that the computed moment set will stay in this space.

Let us remark that the volume variable is used for the NDF as well as for the moment definition. However, for size reduction, the variable r = v 1/3 is also interesting since it decreases in an affine way. This is why integer moments for this variable will also be considered, leading to fractional moments in v and/or a reconstruction of the NDF as a function of r (see the corresponding change of variable in Appendix A).

Moment space and realizability

In order to give a clear picture of the moment space and the realizability problem, we recall some knowledge from the theory of moments. For simplicity, we drop time t hereinafter in this section.

For hybrid methods, only two moments m 0 and m 1 are considered and the space in which they live is clearly

M 1 (V min , V max ){(m 0 , m 1 ) ∈ R 2 + V min m 0 ≤ m 1 ≤ V max m 0 }. If m 0 is positive, then the equalities V min m 0 = m 1 and V max m 0 = m 1 correspond
to the degenerate cases where the NDF is a Dirac delta function at V min or V max , respectively. There is no reason why such cases should appear if V min and V max do not coincide with the nuclei size V nuc . Only the degenerate case (m 0 , m 1 ) = (0, 0) will then be considered, leading to a zero NDF.

For QBMM, the space (0, ∞) is considered for the interval (V min , V max ), as well as moments of order 0 to N on this interval in the variable ξ, which is the volume v or the variable r = v 1/3 . For any positive measure dµ(ξ), i.e. induced by a nondecreasing function µ(ξ) on [0, ∞), let us denote m N (dµ) the vector of moments defined by

m N (dµ) = (m 0 (dµ), m 1 (dµ), . . . , m N (dµ)) T , m k (dµ) = ∞ 0 ξ k dµ(ξ), k ≥ 0, (7)
assuming that such moments are finite. One then has the following definition of the moment space:

Definition 2.1 (Moment space). The moment space M N (0, ∞) is defined as the set of all moment vectors m N (dµ), where dµ(ξ) is a positive measure having finite moments of order 0 to N .

To simplify the notation, M N (0, ∞) is also denoted M N . The moment space can also be defined by the moments of the probability measure [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF]34]. Let us denote it by MN in this case. There is a one-to-one relation between M N -(0, . . . , 0) T and MN through the division by the zero-order moment. Here, the zero moment vector can be used and this is why we employ the definition for positive measures. Moment space is convex [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF] and can be characterized by the Hankel determinants defined by Since the boundary of moment space can be attained (monodisperse or discrete polydisperse cases with a small number of sizes), it is important to be able to determine N (m N ) from a numerical point of view. In practical applications, such boundary detection with Hankel determinants is costly and can be inaccurate when the moments are close to the boundary, due to numerical errors. A more efficient algorithm can be extracted from the theory of orthogonal polynomials [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF][START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF][START_REF] Chihara | An Introduction to Orthogonal Polynomials[END_REF]. A sequence of polynomials {P k } i∈N , where P k is of exact degree k and orthogonal with respect to some positive measure dµ(ξ) on the support [0, ∞) if and only if ∞ 0 P k (x) P l (x) dµ(x) = 0 for any k = l. It is then well known that this sequence satisfies a three-term recurrence relation of the form

H 2n+d = m d . . . m n+d . . . . . . . . . m n+d . . . m 2n+d (8)
P k+1 (x) = (x -a k+1 )P k (x) -b k+1 P k-1 (x), k ∈ N, (9)
with b k+1 > 0 and P -1 (x) = 0, P 0 (x) = 1. Conversely, if the sequence of polynomials satisfies ( 9) with b k+1 > 0 for all k ∈ N, then there exists a measure on the real line for which the polynomials are orthogonal. It is also well known [START_REF] Chihara | An Introduction to Orthogonal Polynomials[END_REF] that the measure µ is supported on [0, ∞) if and only if there exists a sequence of positive numbers {ζ n } n∈N * such that the coefficients in the recurrence relation ( 9) satisfy for all k ≥ 1: Moreover, the link between {ζ n } n∈N and Hankel determinants [START_REF] Dette | The Theory of Canonical Moments with Applications in Statistics, Probability, and Analysis[END_REF] is given by

b k = ζ 2k-1 ζ 2k , a k = ζ 2k + ζ 2k+1 . (10)
ζ k = H k H k-3 H k-1 H k-2 (11)
where we use

H k = 1 if k ≤ 0.
Then, one has the following properties:

Proposition 2.2 (Realizability). Let us consider a moment vector m N and the corresponding

(ζ k ) k=1,...,N . Then m N is in the interior Int(M N ) of moment space if and only if ζ k > 0 for k = 1, . . . , N . Moreover if m N is on the boundary ∂M N of moment space, with N (m N ) = n then ζ k > 0 for k = 1, . . . , n -1 and ζ n = 0.
Several algorithms can be employed to compute efficiently the recurrence coefficients a k and b 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 quadratures[END_REF]. Since it is found to be slightly more stable in practice [START_REF] Wheeler | Modified moments and Gaussian quadratures[END_REF], the last one is used here and referred to as the Chebyshev algorithm (see Appendix C for the description of the ζ-Chebyshev algorithm, coupling this Chebyshev algorithm to the computation of the ζ k ). Moreover, one can remark that, given a moment vector m 2n-1 ∈ Int(M 2n-1 ), the quadrature can be evaluated based on the tridiagonal Jacobi matrix [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF] formed by the recurrence coefficients a k and b k . Thus, we obtain a consistent and robust way to compute the quadrature and verify the realizability of an arbitrary moment set, including boundary detection.

In the rest of the paper, the less rigorous notation f (ξ)dξ is used for the measure, instead of dµ(ξ), making the NDF apparent.

EQMOM

When considering QBMM, the EQMOM reconstruction is able to close the moment transport equations [START_REF] Rigopoulos | Population balance modelling of polydispersed particles in reactive flows[END_REF], since it gives a value of the NDF at zero size. Moreover, it is able to degenerate to a sum of weighted Dirac delta functions in the case where the moment vector is on the boundary of moment space. In this section, the principle of this method is recalled. Moreover, a more robust and efficient EQMOM moment-inversion algorithm than the one given in [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF] is developed here, able to deal with the boundary of moment space.

In EQMOM, the NDF is represented by weighted sum of KDF [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF]:

f (ξ) = N α=1 w α δ σ (ξ, ξ α ) ( 12
)
where σ is a unique nonnegative parameter shared by all KDF. This representation also captures the tail of the distribution corresponding to large particles. To deal with the boundary of moment space, the following informal convergence is imposed through the choice of the KDF: lim σ→0 δ σ (ξ, ξ α ) = δ(ξξ α ). Let us define ξ k,α , m n = (m 0 , . . . , m n ) T and m * n = (m * 0 , . . . , m * n ) T , respectively, using

ξ k,α = ∞ 0 ξ k δ σ (ξ, ξ α ) dv, m k = N α=1 w α ξ k,α , m * k = N α=1 w α ξ k α . (13)
A second constraint on the KDF is that for any k ≥ 1, there exists an invertible matrix A k (σ), independent of the weights w α and abscissas ξ α such that m k = A k (σ)m * k . The purpose of this constraint is to allow for the use of the quadrature based on the Chebyshev algorithm. Indeed, for

any value of σ, one can compute m * 2N -1 (σ) = A 2N -1 (σ) -1 m 2N -1
, and then use the quadrature algorithm to compute the weights (w α (σ)) N α=1 and abscissas (ξ

α (σ)) N α=1 if m * 2N -1 (σ)
is in the interior of moment space. The moments of orders 0 to 2N -1 of the EQMOM reconstruction corresponding to these parameters for a given σ, (w α (σ)) N α=1 and (ξ α (σ)) N α=1 are m 0 to m 2N -1 , and the value of σ has to be adapted in order that its 2N th-order moment is m 2N (which is not always possible).

The original EQMOM moment-inversion algorithm is given in [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF]. There are two open issues with this algorithm. First, the boundary of moment space was not really dealt with, as well as the transition with the interior of moment space since only moment vectors far from this boundary were considered. This point is however essential here for the robustness and accuracy of our computations, considering the typical solution of the PBE described in Sec. 2.1. Second, the iterative algorithm for the computation of σ was not optimal. This is why we propose an improved version of the algorithm, summarized in Fig. 1.

Here, if m 2N is in the interior of moment space, one defines the function m2N of σ in the following way: if If m * 2N -1 (σ) is not in the interior of moment space, then the value of σ is invalid. The function m2N (σ) is then set to a very high value, e.g. 10 100 , so that this case will be automatically eliminated.

m * 2N -1 (σ) = A 2N -1 (σ) -1 m 2N is in the interior Int(M 2N -1 )
Hence, the value of σ is obtained by solving the following scalar nonlinear problem:

D 2N (σ) = 0 where D 2N (σ) = m 2N -m2N (σ). ( 14
)
Let us remark that, in the case N = 2, an analytical condition on σ of the form σ ≤ σ

max for the realizability of m * 2N -1 (σ) can usually be obtained, ensuring the non-negativity of the corresponding Hankel determinants H * 2 and H * 3 . Moreover, D 2N (0) is positive if m 2N is in the interior of moment space and the problem ( 14) either has a solution or a nonlinear solver will give the value corresponding to an upper limit σ max for σ, when a singularity appears for the function D 2N (σ). In this last case, one just has to minimize the error on the last moment by minimizing D 2N (σ) 2 .

The global algorithm of EQMOM reconstruction from a moment vector m 2N , able to deal with the boundary of moment space, is then:

1. Determine N (m 2N ) by computing the (ζ k ) k=1,...,2N with the ζ-Chebyshev algorithm. 1.1. If N (m 2N
) is an odd number 2n -1, then σ = 0 and the quadrature points are obtained from three-term recurrence coefficients given by the (ζ k ) k=1,...,2n-1 . 1.2. Otherwise, let N (m 2N ) = 2n where n ≤ N and move to step 2.

f (ξ) 6 Z Z Z Z Z Z α k Ξ k-1 = Ξ (k) a β k = 0 Ξ (k) b Ξ k f (ξ) 6 α k β k Ξ k-1 = Ξ (k) a Ξ k = Ξ (k) b f (ξ) 6 β k Ξ k-1 α k = 0 Ξ (k) a Ξ k = Ξ (k) b
2. Solve the scalar nonlinear problem [START_REF] Kumar | Improved accuracy and convergence of discretized population balance for aggregation: The cell average technique[END_REF] in the interval (0, σ

max ) using Ridders' method [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF]. 2.1. If D 2n (σ) = 0 then go to step 4.

2.2. Otherwise, in case of not being able to obtain σ, set σ max = σ and go to step 3. 3. Use Brent's method [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF] on [0, σ max ] to minimize D 2n (σ) 2 . 4. Once the parameter σ is obtained, the weights and abscissas are computed with the quadrature algorithm, from m * 2n-1 = A 2n-1 (σ) -1 m 2n-1 . Two types of KDF are considered for the Stieljes problem: gamma [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF] and log-normal [START_REF] Madadi-Kandjani | An extended quadrature-based moment method with log-normal kernel density functions[END_REF].

The matrix A k (σ) can then be given explicitly, as well as the interval on which σ is located. This is done in Appendix D.

TSM method

For TSM, a discretization 0 = V 0 < V 1 < . . . < V Ns = +∞ is introduced. On each section [V k-1 , V k ],
the moments of orders 0 and 1 in volume are considered:

N k = V k V k-1 f (v) dv, M k = V k V k-1 v f (v) dv. ( 15
)
Physically, they are usually the main variables of interest since the global moment of order 0 is conserved by the growth term and the global moment of order 1 is conserved by aggregation, for example [START_REF] Kumar | On the solution of population balance equations by discretization-I. A fixed pivot technique[END_REF]. This is one of the advantages of the hybrid TSM method over the sectional methods, where a single variable has to be chosen, thus usually loosing some local conservation properties. Moreover, the last section here is unbounded: a maximal size does not have to be determined as for the standard sectional methods.

In practice, the change of variable v = r 3 can be interesting to consider when dealing with size reduction. Indeed, the variable r = v 1/3 then decreases at a constant rate, making it the best variable for describing such phenomena [START_REF] Laurent | Numerical analysis of Eulerian multi-fluid models in the context of kinetic formulations for dilute evaporating sprays[END_REF][START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF]. The considered variables are then moments of orders 0 and 3 of the corresponding NDF: f r (r) (with f r (r)dr = f (v)dv):

N k = R k R k-1 f r (r) dr, M k = R k R k-1 r 3 f r (r) dr, ( 16
)
with R 3 k = V k for k ≥ 0. The two types of reconstruction are considered in this work. Let us then denote ξ the considered variable (ξ = v or ξ = r = v 1/3 ) and Ξ k the corresponding bounds of the sections (

Ξ k = V k or Ξ k = R k = V 1/3 k
). The principle of the reconstruction is given in [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF], with a different choice for the variable (r 2 in this article). The NDF, as a function of ξ, is then approximated by an affine by part function on each section:

f ξ (ξ)| [Ξ k-1 ,Ξ k ] = α k + (β k -α k ) ξ -Ξ (k) a Ξ (k) b -Ξ (k) a 1 [Ξ (k) a ,Ξ (k)
b ] (ξ).

(

) 17
Among the four parameters α k , β k , Ξ

a and Ξ (k) b used for the reconstruction, two are always fixed, depending on the three possible cases, which are illustrated in Fig. 2:

1. Ξ (k) a = ξ k-1 , Ξ (k) b < ξ k and β k = 0. 2. Ξ (k) a = ξ k-1 and Ξ (k) b = ξ k . 3. Ξ (k) a > ξ k-1 , Ξ (k) b = ξ k and α k = 0.
The two other ones are determined in such a way that the moments of orders 0 and

1 if ξ = v or 3 if ξ = r on [Ξ k-1 , Ξ k ] of f ξ are N k
and M k , respectively. They can be computed analytically and the formula are given in Appendix E.

Continuous particle-size changes

We now demonstrate accurate and realizable schemes for the description of continuous particlesize changes. They are based on the scheme developed for spray evaporation with moment methods by Massot et al. [START_REF] Massot | A robust moment method for evaluation of the disappearance rate of evaporating sprays[END_REF], and adapted to TSM by Laurent et al. [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF]. However, for moment methods, this scheme was developed for an even number of moments. An adaptation is done here, based on an interpretation of the scheme as an evaluation by a quadrature method of the integrals deduced from a kinetic scheme, similarly as the one given in [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF] in the case of a constant size decreasing rate. Therefore, the kinetic scheme and the analytical solution of the PBE on which it is based are first recalled. Then, the quadrature-kinetic scheme (QKS) is given for both EQMOM and TSM. Finally, the ability of the two methods to describe continuous particle-size changes is then evaluated, separately for size reduction and growth.

Analytical solution and kinetic scheme

The rate is denoted by R in general for both growth and size reduction. The PBE then reads

∂ t f (t, ξ) + ∂ ξ R(t, ξ) f (t, ξ) = 0 ( 18
)
where ξ can be volume v or the variable r = v 1/3 . If f (0, ξ) = f 0 (ξ) is the initial NDF, then the analytical solution of this equation is given by

f (t, ξ) = f 0 (Ξ(0; t, ξ))J(0; t, ξ) (19)
where the characteristics Ξ(t; s, ξ) of ( 18) are defined by the evolution of particle size:

dΞ(t; s, ξ) dt = R(Ξ(t; s, ξ)) with Ξ(s; s, ξ) = ξ, (20)
and J(t; s, ξ) is the Jacobian of the transformation ξ → Ξ(t; s, ξ). For the considered growth and size reduction models, these analytical solutions are provided in Appendix A.

The principle of the kinetic scheme is to use the exact solution of the PBE inside a time step. Then, from the moments at time t n , a NDF f n (ξ) is reconstructed for ξ ∈ R + , using the algorithms given in Sec. 2.4 for EQMOM and in Sec. 2.5 for TSM. The exact solution of the PBE (18) between t n and t n+1 starting from this NDF is f (t, ξ) = f n (Ξ(t n ; t, ξ))J(t n ; t, ξ). The moments at time t n+1 are deduced from the NDF f (t n+1 , ξ). Since they are moments on some interval (ξ min , ξ max ) of the NDF, they are written as

ξmax ξmin ξ k f n (Ξ(t n ; t n+1 , ξ))J(t n ; t n+1 , ξ)dξ = Ξ(tn;tn+1,ξmax) Ξ(tn;tn+1,ξmin) (Ξ(t n+1 ; t n , ζ)) k f n (ζ)dζ ( 21
)
where the second expression is obtained thanks to the change of variable ζ = Ξ(t n ; t n+1 , ξ), equivalent to ξ = Ξ(t n+1 ; t n , ζ). Depending of the complexity of R(t, ξ), the solution and its moments are not always easy to compute analytically. This is why the QKS scheme is used.

QKS

The principle of QKS is the same as the kinetic scheme, except that the integrals ( 21) are now evaluated thanks to a Gauss quadrature of the measure f n (ζ)1 [Ξ(tn;tn+1,ξmin),Ξ(tn;tn+1,ξmax)] (ζ)dζ. Indeed, from the moments of orders 0 to 2N +1 of f n on the interval [Ξ(t n ; t n+1 , ξ min ), Ξ(t n ; t n+1 , ξ max )], quadrature points {w α , ξ α } N +1 α=1, are determined and the previous integral is approximated by

N +1 α=1 w α Ξ(t n+1 ; t n , ξ α ) k . (22)
The number N is equal to 1 for TSM, and corresponds to the number of KDF used for EQ-MOM. Using at least N + 1 quadrature points is essential so that the moment vector at time t n+1 is in the interior of the corresponding moment space. Moreover, if the abscissas ξ α belong to (Ξ(t n ; t n+1 , ξ min ), Ξ(t n ; t n+1 , ξ max )), then Ξ(t n+1 ; t n , ξ α ) belongs to (ξ min , ξ max ) and the scheme is realizable. Then, the ODE [START_REF] Qamar | On the solution of population balances for nucleation, growth, aggregation and breakage processes[END_REF] defining the characteristics has to be solved once in reverse time to compute each bound Ξ(t n ; t n+1 , ξ min ), when the result is nontrivial, and N + 1 times to determine the evolution of the quadrature abscissas. In a one-way coupling context, this can be done easily, even if the drift rate is time dependent due to a dependence on the gas variables. In a two-way coupling context with a time-dependent drift rate, either the evolution of the drift term during the time step can be estimated or it is set constant.

QBMM

For QBMM, the evolution of the moments due to continuous size reduction is as follows:

1. Given the moment set m n k 2N k=0 at time t n , reconstruct a NDF f n (ξ).

Compute the moments mk

2N +1 k=0 of f n on [Ξ(t n ; t n+1 , 0), ∞): mk = ∞ Ξ(tn;tn+1,0) ξ k f n (ξ)dξ.
m n+1 k = N +1 α=1 w α Ξ(t n+1 ; t n , ξ α ) k , k = 0, . . . , 2N. (23)
Let us remark that the moments mk are computed analytically. The moments of the NDF between 0 and Ξ(t n ; t n+1 , 0) correspond to the disappearance fluxes. When considering surface growth, these fluxes are zero and mk = m k . The same algorithm can be used, then computing one more moment m 2N +1 and the N corresponding quadrature points. However, instead of these quadrature points, the secondary quadrature defined in Sec. 4 can be used, as soon as the number of secondary quadrature points N α is greater than N + 1. This was done in [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF], which required solving the ODE N × N α times. The use of secondary quadrature points with a suitability chosen N α is especially recommended for cases where the function R(t, ξ) takes on both positive (growth) and negative (size reduction) values for varying ξ, as soon as there is no disappearance fluxes.

TSM

For TSM, the scheme is similar, but dividing the interval [Ξ(t n ; t n+1 , ξ min ), Ξ(t n ; t n+1 , ξ max )] into two sub-intervals on which the reconstructed NDF is smooth. This allows for a better accuracy of the scheme [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF]. Let us denote by V (ξ) the volume corresponding to the variable ξ (i.e. V (v) = v and V (r) = r 3 ) and Ξ k the bounds of the sections in terms of the variable ξ, in such a way that V (Ξ k ) = V k . Then, the evolution of the moments through continuous size reduction is as follows:

1. Given the moments {N n k , M n k } Ns k=1 at time t n , reconstruct a NDF f n (ξ). 2. Compute the moments of orders 0 to 3 of f n on each interval [Ξ(t n ; t n+1 , Ξ k-1 ), Ξ k ] and on [Ξ k , Ξ(t n ; t n+1 , Ξ k )] and the corresponding quadrature points {w

(k,I) α , ξ (k,I) α } 2 α=1 and {w (k,II) α , ξ (k,II) α } 2 α=1 .
3. Update the moments on section k:

N n+1 k = 2 α=1 w (k,I) α + 2 α=1 w (k,II) α and M n+1 k = 2 α=1 w (k,I) α V Ξ(t n+1 ; t n , ξ (k,I) α ) + 2 α=1 w (k,II) α V Ξ(t n+1 ; t n , ξ (k,II) α ) .
A similar scheme is employed for growth, but using the two intervals [Ξ(t n ;

t n+1 , Ξ k-1 ), Ξ k-1 ] and [Ξ k-1 , Ξ(t n ; t n+1 , Ξ k )] in step 2.
Let us remark that the time step ∆t n = t n+1t n has to be limited in order for the obtained moment set to stay in moment space. Indeed, the value Ξ(t n+1 ; t n , ξ) has to stay in

[Ξ k-1 , Ξ k ] if the variable ξ is in [Ξ(t n ; t n+1 , Ξ k-1 ), Ξ(t n ; t n+1 , Ξ k )].
In the case of a size-reduction rate given by [START_REF] Piton | Simulation of fine particle formation by precipitation using computational fluid dynamics[END_REF] or for the growth rate given by (4), this is realized, respectively, through the following conditions: max k∈{1,...,Ns}

c red ∆t n 3(V 1/3 k -V 1/3 k-1 ) ≤ 1 max k∈{1,...,Ns} 2c gro ∆t n 3(V 2/3 k -V 2/3 k-1 ) ≤ 1. ( 24
)
The coefficient on the left-hand side of each inequality is called here the CFL number.

Numerical results for size reduction

We discuss in this section numerical results for continuous size reduction using EQMOM and TSM, with an initial NDF typical of the ones induced by nucleation and aggregation. The test case is described in Appendix B.1. The simulations are done for time t between 0 and T = 0.1 and compared to the analytical solution. Some global quantities are considered: the global v-moments m 0 , m 1 , m 2 , the average volume m 1 /m 0 , and the variance

σ 2 v = m2m0-m 2 1 m 2 0
. The maximum values in time of the errors between numerical q nu and exact q ex quantities, normalized by the maximal value of the exact solution are computed: Moreover, the reconstructed NDF, f nu (t, v), as well as the corresponding volume density function (VDF) vf nu (t, v) are also compared with the exact solution, f ex (t, v) and vf ex (t, v), respectively. The maximal value of the L 1 norm of the difference between the computed and exact functions, normalized by the maximal value of the L 1 norm of the exact solution, is then computed:

err q = max t∈[0,T ] |q nu (t) -q ex (t)| max t∈[0,T ] q ex (t) . ( 25
err N DF = max t∈[0,T ] ∞ 0 |f nu (t, v) -f ex (t, v)| dv max t∈[0,T ] ∞ 0 f ex (t, v) dv (26) err V DF = max t∈[0,T ] ∞ 0 v |f nu (t, v) -f ex (t, v)| dv max t∈[0,T ] ∞ 0 vf ex (t, v) dv (27)

Use of v-reconstruction

Computations are done with gamma-and Ln-EQMOM with a reconstruction in v using N = 4 KDF. The considered moments in v are (m k ) 8 k=0 . TSM is also used with a uniform discretization in the volume between V 0 = 0 and V Ns-1 = V max , the last section being [V max , ∞). The number of sections is N s = 41 and V max = 50. QKS is used with a constant time step given by ∆t = 0.001, corresponding to a CFL of 0.05 for TSM. A second time step ∆t = 0.018 is also used for TSM, corresponding to a larger CFL, equal to 0.9.

In all the cases, the global moments of orders 1, 2 and 3 (not shown) are well reproduced by all methods. However, the errors on the 0 th -order moment, plotted in Fig. 3 to t = 0.033), one can see that the shape of this function is well captured by TSM, except an accumulation at zero size in the first section, whereas a shift is observed for EQMOM. This is due to the singularity in the exact solution (B.2) at v = 0. Indeed, the reconstructed NDF is regular at v = 0 for TSM, equal to zero at v = 0 for Ln-EQMOM and either equal to zero or singular at v = 0 for gamma-EQMOM. In all cases, the behavior at v = 0 is very different compared to the exact solution, which behaves like v -2/3 . This means that even if the considered moments had good values, the reconstruction close to zero would be far from the real one and the evolution of m 0 would then be badly reproduced. This also explains the accumulation of particles of small size until some time-step dependent time when they suddenly disappear. For gamma-EQMOM, the effect is smoother due to the possible singularity at zero. For TSM, it concerns only the first section, which is why the NDF is well reproduced in this case, except close to zero, even with a wrong value for the global 0 th -order moment. In order to overcome this issue, we use the r = v 1/3 variable so that the exact NDF turns out to be advection of the initial NDF in phase space.

Numerical results with r-reconstruction

The computations are repeated with gamma-and Ln-EQMOM using the variable r and N = 2, 3 and 4 KDFs. The considered moments in v are (m k/3 ) 2N k=0 and f r is now reconstructed. TSM is also used with the same moments as in the previous section, but with a reconstruction in r and several uniform discretizations in r between R 0 = 0 and R Ns-1 = R max , the last section being [R max , ∞), where R max = V 1/3 max and V max = 50. The OSM method, is also used, considering the same kind of discretization as TSM (without the last section, which is almost empty), using the variables M k and a constant reconstruction in the variable r. For all methods, QKS is employed with a constant time step, ∆t = 0.001, small enough so that the simulations are converged in time.

The errors induced by the EQMOM simulations are gathered on Table 1, showing first that these methods reproduce accurately the moments of orders 1 and 2 in v. The error on the 0 th -order moment is smaller than with the v-reconstruction, but it is still between 14 and 23 % for Ln-EQMOM, and 10 and 16 % for gamma-EQMOM. When looking at the evolution of m 0 for N = 4 (see Fig. 5), one observes a discontinuous behavior for Ln-EQMOM, with four discontinuities, corresponding to some discontinuities of the abscissas (see Fig. 6-left). Indeed, in this case, the reconstructed NDF is always equal to zero at zero size. The corresponding fluxes are then very small, especially if the time step is small. But size reduction induces in particular a decrease of the first abscissa and then a concentration of the first KDF close to zero, making it suddenly disappear when it is too small (and the moment set is then briefly close to the boundary of the moment space, also showing the robustness of the reconstruction algorithm). With gamma-EQMOM, the evolution of m 0 as well as the abscissas is smoother (see Fig. 6-right) thanks to the possibility of the reconstruction to not be zero, even if it is singular in this case (but the flux itself is not singular thanks to the integration over a small interval).

Concerning TSM, errors for three discretizations (N s = 5, 9, 41) are given in Table 2. When compared to gamma-EQMOM with N = 4 (using 9 moments), which gives the best results among the EQMOM simulations, TSM with only 5 sections (10 variables) is less accurate for m 1 and m 2 , but gives slightly better accuracy on m 0 and on the mean volume m 1 /m 0 , with errors close to 8 % and 4 % respectively. The accuracy, however, rapidly increases with the number of sections, as shown in Fig. 7, where the errors on m 0 , m 1 and m 2 are plotted as functions of the section size, using 3 to 65 sections. One then sees at least second-order accuracy for TSM. Moreover, in Table 2, errors induced by OSM for N s = 10, 18, 82 are shown, thus using double the number of sections and then the same number of variables as TSM. The errors on all quantities are higher than with TSM and the convergence, when increasing the number of sections, is slow (first-order accuracy for m 1 and m 0 and smaller order for m 2 , as shown on Fig. 7).

The NDF f (t, v) at t = 0.033 is plotted on Fig. 8, from the reconstruction f r (t, r) obtained from EQMOM and TSM with f (t, v) = f r (t, v 1/3 )/(3v 2/3 ). For N = 4, the reconstruction obtained with gamma-EQMOM is much better than when considering v-reconstruction and quite similar to the reconstruction with Ln-EQMOM, except near zero. Both are close to the exact solution, with some shift for the second mode of this bimodal function. The error on this distribution and on the VDF (eliminating the singularity at zero) is indeed quite good for gamma-EQMOM (see Table 1) and also for Ln-EQMOM. When considering TSM, this error is higher than with EQMOM using N = 4 KDFs, except with more than 21 section and the NDF is very well reproduced with 41 sections. However, TSM is more accurate than OSM, which is not able to reproduce correctly the NDF with 82 sections, the error on the NDF and VDF being then larger than the one with gamma-EQMOM.

In summary, when considering size reduction, one can see from the examples presented in this section that gamma-EQMOM is always more accurate than Ln-EQMOM. Moreover, gamma-EQMOM gives a quite good estimate of the NDF. The convergence of EQMOM in terms of the number of moments is, however, quite slow and induces an increasing complexity of the moment space, so that considering more than nine moments (i.e. N = 4) may not be worthwhile for most applications. For better accuracy, TSM can be used, with a larger number of variables, which would have to be transported in space in many applications. However, the good convergence of TSM allows the use of a much smaller number of sections (and also of variables) than OSM for the same accuracy.

Numerical results for diffusion-controlled growth

The test case described by McGraw [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF], and described in Appendix B.2 is used here to test the schemes for diffusion-controlled growth. Simulations are done between time 0 and T = 20s, with EQMOM and TSM using r-reconstruction, and OSM using a constant reconstruction in v (the method is unstable with the constant reconstruction in r). For OSM and TSM, a uniform discretization is used in the r variable, between 0 and R max = 60 µm. Moreover, for all methods, QKS is used with a constant time step, ∆t = 0.01, small enough so that the simulations are converged in time. The same kinds of post-treatments are done as with the size-reduction case, except that the NDF in r in now considered.

The normalized errors for all methods are given in Tables 3 and4. For Ln-and gamma-EQMOM, the moment errors are very small, even with N = 2. One can remark that the growth model conserves the zeroth-order moment and yields closed equations for even-order r-moments. For EQMOM, the even-order r-moments are thus very accurately reproduced, e.g. m 2 when N = 4, except eventually for the last moment since it is not always well reproduced by the reconstruction. For TSM, the zeroth-order moment is also conserved, such that its error is very small and the accuracy of the moments of orders 1 and 2 is high, even with only 15 sections. For OSM, a large number of sections is needed to reproduce the moments (even at zero order, which is not conserved) with an error of a few percent. The error curves for the moments found with OSM and shown in Fig. 9 as functions of the section size, using 3 to 129 sections. It can be clearly observed that the convergence rates are at least first and second order, respectively, for OSM and TSM. Results for the reconstructed NDFs at time T are shown in Fig. 10. As can be observed, the EQMOM results are very similar for Ln-and gamma-EQMOM, and in relatively good agreement with the exact solution. The errors on the NDF and VDF (see Table 3 ) do not really depend on the number of KDFs and are about 30 or 40 %. On the other hand, the OSM result does not capture correctly the sharp jump in NDF at r ≈ 5.5, except with a very large number of sections (e.g. 500). In contrast, TSM with 45 sections does a good job of reproducing the exact NDF.

In summary, when considering surface growth, one can see from the examples presented in this section that gamma-and Ln-EQMOM give equivalent results, which are very good for the moments and relatively good for the NDF, and not significantly improved using a larger number of moments. Discretized methods, especially OSM, are disfavored by the fact that the NDF is shifted toward the larger sizes and by the sharp jump of the exact NDF. However, TSM is able to reproduce well moments and the sharp NDF with a quite small number of sections, much smaller than for OSM.

Aggregation

In this section we apply EQMOM and the discretized methods to solve a PBE with only aggregation. Such applications are known to be particularly challenging for discretized methods because the mean particle size and variance grow continuously with time. From a numerical perspective, the aggregation operator has an integro-differential form that strongly couples all points in size phase space. In real applications, numerical simulation of the aggregation term is often the most computationally expensive operation.

EQMOM

As first shown by Marchisio et al. [START_REF] Marchisio | Quadrature method of moments for aggregationbreakage processes[END_REF], QBMM are particularly well suited for aggregation because the abscissas evolve in phase space to adapt to the changing shape of the NDF. After using the EQMOM moment-inversion algorithm, we obtain a closure for ( 5) and [START_REF] Rigopoulos | Population balance modelling of polydispersed particles in reactive flows[END_REF]. Solving these 565 moment equations requires the numerical approximation of the integrals in (6) since the KDF δ σ (ξ, ξ α ) is a smooth function. This can be done by using the quadrature specific to each KDF:

∞ 0 g(v) δ σ (ξ, ξ α ) dv ≈ Nα β=1 ω αβ g(v αβ ) (28)
where N α is the number of secondary quadrature points. N α must be larger than N (in practice N α = N + 1), in such a way that ( 28) is exact when g is a polynomial of degree less than or equal to 2N . These weights ω αβ and abscissas v αβ can be easily computed from the known recurrence coefficients of the orthogonal polynomials corresponding to the KDF: the generalized Laguerre [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF] and Stieltjes-Wigert [START_REF] Wilck | A general approximation method for solving integrals containing a lognormal weighting function[END_REF] for gamma-and Ln-EQMOM, respectively. The approximation of the integral for any arbitrary function with respect to the reconstructed NDF is then

∞ 0 g(v) f (v) dv = N α=1 Nα β=1 w αβ g(v αβ ) (29)
where w αβ = w α ω αβ . Finally, we obtain the system of moment equations for aggregation, for k ∈ N:

dm k dt = 1 2 N α1=1 Nα 1 β1=1 N α2=1 Nα 2 β2=1 w α1β1 w α2β2 (v α1β1 + v α2β2 ) k -v k α1β1 -v k α2β2 β(v α1β1 , v α2β2 ) (30)
A realizable ODE solver must be used to solve this system.

TSM

The equations for the moments in one section [V min , V max ] = [V k-1 , V k ] are given by ( 5) and ( 6). For aggregation, the integration domain

Ω k = {(v * , v ) > 0 V k-1 < v * + v < V k } then appears.
Since the closure consists in a reconstruction of the NDF on each section, this domain is divided into some elementary sub-domains: [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF]. One then defines the following elementary integrals:

D ijk = Ω k ∩ ([V i-1 , V i ] × [V j-1 , V j ])
Q n ijk Q * ijk = D ijk 1 v * β(v * , v ) f (v * ) f (v ) dv * dv , (31)
The equations for the moments of orders zero and one in v in section k are then

∂ t N k = 1 2 Ns i=1 Ns j=1 Q n ijk - Ns i=1 Ns j=1 Q n kij (32a) ∂ t M k = Ns i=1 Ns j=1 Q * ijk - Ns i=1 Ns j=1 Q * kij (32b)
The elementary integrals are computed using a 5-point Gauss-Legendre quadrature in the variable r, in each direction. Let us also remark that in the case of a time-independent collision kernels, a numerically efficient way to compute the source terms can be developed [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF], as for the sectional method in the spray context [START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF], using pre-computation of some terms and a compact storage of the variables. Finally, as for QBMM, the realizability constraint has to be respected in the resolution of this ODE system.

Realizable, adaptive time step

Due to the complexity of aggregation kernels in practical problems [START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF], the ODE system can be stiff and we then want to use an adaptive time step to limit the cost of the computations. Concerning realizability, since the aggregation operator does not move the moment vector to the boundary of moment space, explicit Euler methods with small enough time steps can be used as well as any convex combination of Euler explicit time steps, such as the strong-stability-preserving (SSP), explicit Runge-Kutta methods [START_REF] Gottlieb | Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations[END_REF]. We therefore use an adaptive time-step algorithm, based on embedded SSP explicit Runge-Kutta methods, with a time-step selection designed both to control the error and to ensure the realizability of the moment set (see Appendix F for details).

Ln

Numerical results

We present results in this section for two aggregation test cases from [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF]: one, described in Appendix B.3, using a sum kernel (case 3) and the other one, described in Appendix B.4, using the Brownian kernel (case 4). For TSM, we consider the following geometric discretization:

V 0 = 0, V k = V max (RN s ) Ns-k+1 Ns-1 , k ∈ {1, . . . , N s -1}, V Ns = ∞ ( 33
)
where R is a user-defined parameter, as well as V max : here, V max = 25000 and R = 5000. The same kind of discretization is used for OSM, without the last section and then with a larger value for V max (V max = 30000, R = 10 4 ). A dependence on the number of sections is then introduced in the ratio between two successive section sizes, (RN s )

1

Ns-1 , in such a way that the discretization tends to a uniform discretization when N s tends to infinity. Aggregation with the sum kernel yields quantitatively similar results for r-and v-moments. In general, the results found with v-reconstruction are slightly more accurate, but since the r variable gives better accuracy for other phenomena (e.g. size reduction), all results presented here are obtained using r-reconstruction. The normalized errors are given in Tables 5 and6. For the moments, EQMOM yields very small errors for both KDFs, even with N = 2. These results are consistent with previous work using QBMM [START_REF] Marchisio | Quadrature method of moments for aggregationbreakage processes[END_REF][START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF] and are due to the fact that for the sum kernel, the equations on the moments are closed, except for the last one. In comparison, for OSM and TSM the errors on the moments are significantly larger than with EQMOM. This difference is mainly due to the nature of aggregation where the NDF moves to larger and larger r as time progresses (if it is not in competition with size-reduction phenomena, as it might be in some applications). With QBMM, the abscissas follow the movement of the NDF in phase space, while for discretized methods the discretization of phase space is fixed and the value of V max has to be adapted to the final time. Nevertheless, for finite times, TSM with as few as 15 sections provides a much more accurate reconstruction of the NDF as compared to EQMOM.

For Brownian aggregation, the time evolutions of selected moments are shown in Fig. 11, and the NDFs at t = 10 are given in Fig. 12. For the moments, EQMOM with either KDF yields very accurate predictions (not shown), confirming previous observations [START_REF] Marchisio | Quadrature method of moments for aggregationbreakage processes[END_REF]. OSM converges to the reference solution after approximately 90 sections, while TSM gives accurate predictions with as few as 15 sections. With a small number of sections, OSM tends to over-predict the variance σ 2 v by a large amount. This is due to the behavior of the reconstructed NDF in the largest sections, which for OSM decays too slowly with increasing v. For this reason, the errors in moments of orders greater than two are significant, making OSM unsuitable for aggregation processes. The reconstructed NDFs in Fig. 12 confirm the general trend that gamma-EQMOM yields better results than Ln-EQMOM. The NDF found with gamma-EQMOM is quite accurate, as previously reported [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF]. In order to obtain a comparable NDF with OSM at least 30 sections are required, while for TSM 15 sections yield an accurate NDF. . OSM (bottom left) with 10 (dashed magenta line with ×), 30 (dashed green line with ) and 90 (dashed blue line with •) sections. TSM (bottom right) with 5 (dashed green line with ), and 15 (dashed blue line with •) sections.

5.

Operator time-splitting method for general problems

Strang time-splitting method

We consider in this section an algorithm to solve moment equations with all the types of source terms appearing in fine-particle applications. The idea is to use an operator time-splitting method that solves successively in one time step the nucleation-aggregation and particle-size-change operators. In this work, we use the Strang time-splitting method [START_REF] Strang | On the construction and comparison of difference schemes[END_REF], which is second order in time: we solve successively, in one time step ∆t = t n+1t n , 1. the nucleation-aggregation operator on t n , t n + ∆t 2 2. the particle-size-change operator on [t n , t n+1 ] 3. the nucleation-aggregation operator on t n + ∆t 2 , t n+1 Since the nucleation-aggregation operator is stiffer than the particle-size-change operator, in general, the splitting time step ∆t is chosen to satisfy the stability of the latter. Then, we apply an adaptive time step for the nucleation-aggregation operator. For each operator, we use the realizable and robust schemes given in Secs. 3 and 4 for particle-size-change and nucleation-aggregation processes, respectively.

Numerical results

We investigate here a test case with nucleation, aggregation and size reduction described in Appendix B.5 For this test case, a large number of nuclei are produced at the start, leading to significant particle growth due to aggregation. Eventually, growth is moderated by size reduction, resulting in a maximum in time of the mean particle diameter m 1/3 /m 0 . For sufficiently long times, size reduction would consume all particles. Here, we show results only up to the time where size reduction begins to overtake aggregation (t ≈ 6).

For TSM, we use the following geometric discretization, similar to the one classically used for soot applications with the sectional method:

V k = V V max V k-1 Ns-2 , k ∈ {1, . . . , N s -1}, (34)
with V 0 = 0 and V Ns = ∞. The volume V is slightly larger than the nuclei size ( V = 1.001) and V max = 10 4 . For OSM, the last section is skipped, V max = 2 × 10 4 and the first section is divided in two sections (N s = 10) or three (N s = 30) or seven (N s = 90) where the last one ([0.95 1.05], V = 1.05) contains the nuclei size. All simulations are done with a constant splitting time step ∆t = 0.003. A reference solution is also computed, using TSM with 101 sections and a ten time smaller splitting time step. Except for OSM, all methods yield essentially the same results for m 0 and the mean particle diameter as seen in Fig. 13. With OSM, the mean particle diameter during the nucleation period (t < 0.5) is under-predicted and m 0 is over-predicted even when a relatively large number of sections is used. For the variance σ 2 v , the EQMOM results for N = 3 and 4 are essentially identical, and equal to the TSM result found with 101 sections (see Fig. 14). In contrast, OSM converges slowly towards the reference solution with increasing numbers of sections, while TSM requires approximately 45 sections to attain convergence. As seen in Fig. 15, the reconstructed NDFs differ significantly between the methods. In general, Ln-EQMOM has more obvious modal peaks as compared to gamma-EQMOM, which are associated with the behavior of the KDFs at v = 0. Increasing N from 3 to 4 KDFs shows some improvement (especially for gamma-EQMOM), but the convergence is slow. In contrast, TSM yields a better NDF reconstruction as compared to EQMOM, even with as few as 5 sections, while OSM requires a much larger number of sections to obtain reasonably accurate results. Moreover, with our not yet optimized codes, CPU time for TSM with 15 sections, which also accurately describe the evolution of the moments of order 0 to 2, is of the same order than for EQMOM with N = 4, but using 30 variables instead of 9, which will have to be transported in heterogeneous cases. TSM with 5 sections is a little bit less accurate on these moments but is at least 5 times faster and use 10 variables. Finally, OSM with 30 sections is 4 times slower than TSM with 15 sections for the same number of variables and a much lower accuracy.

Conclusions

In this work, we have developed quadrature-based moment methods (EQMOM with gamma and log-normal KDF) and a hybrid method (TSM) for modeling fine particles undergoing nucleation, aggregation and continuous changes (growth and size reduction). For both classes of methods, efficient algorithms are given to reconstruct the NDF, using analytical formula for TSM and an improved algorithm compared to [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF] for EQMOM, which are able to deal with the boundary of moment space. In addition, realizable numerical schemes are presented for each individual term in the PBE, as well as an operator-splitting approach to resolve the complete problem. This represents a key point for their use in real applications.

In the considered test cases, for EQMOM the gamma KDF leads to more accurate results than the log-normal KDF. Moreover, better results are obtained using fractional (r) moments. For example, more accurate results are found for cases with size reduction and the results for aggregation are not significantly different with fractional moments. Consistent with previous results in the literature [START_REF] Marchisio | Quadrature method of moments for population-balance equations[END_REF][START_REF] Marchisio | Quadrature method of moments for aggregationbreakage processes[END_REF], we find that EQMOM is well suited for accurately predicting the evolution of the moments for general problems, especially when aggregation processes are involved. Furthermore, EQMOM with N = 4 predicts the NDF with a reasonably good accuracy. We also showed that the hybrid TSM method is much more accurate than OSM, a widely used sectional method. Indeed, TSM conserves both the global number and mass of the particles, and can deal with a non compact support of the NDF, even if the size-space discretization must be adapted to the expected range of particle sizes for good accuracy. Because the NDF reconstruction found with TSM is much better in each section, the number of sections needed to obtain good accuracy for the first few moments and the NDF is much smaller as compared to OSM. In general, discretized methods are usually not as efficient as EQMOM for predicting the moments, but TSM is able to accurately predict the NDF with a reasonably small number of sections (from 5 to 15). In contrast, EQMOM converges slowly to the exact NDF as the number of moments increases and can not deal with too large a number of moments due to the increasing complexity of moment space and the numerical accuracy of the Chebyshev algorithm.

Further studies are needed to extend the accurate and realizable numerical schemes presented here to multivariate PBE models for fine-particle formation [START_REF] Raman | Modeling of fine-particle formation in turbulent flames[END_REF][START_REF] Sung | Large-eddy simulation modeling of turbulent flame synthesis of titania nanoparticles using a bivariate particle description[END_REF][START_REF] Mehta | On the role of gas-phase and surface chemistry in the production of titania nanoparticles in turbulent flames[END_REF]. This can be accomplished by since the characteristics are given by Ξ(t; s, v) = v 1/3 -c red 3 (ts)

3 .

(B.3)

Let us remark that, when considering r = v 1/3 instead of v and starting from a r-NDF f r 0 , the solution is f r (t, r) = f r 0 r + c red 3 t (B.4)

M k N k ∈ V k , V (k) min V (k) min , V (k) max V (k) max , V k V (k) a V k-1 V k-1 3 M k N k -2V k V (k) b 3 M k N k -2V k-1 V k V k α k 2N k V (k) b -V k-1 6 (V k -V k-1 ) 2 V (k) max N k -M k 0 β k 0 6 (V k -V k-1 ) 2 -V (k) min N k + M k 2N k V k -V (k) a
Table E.7: Coefficients of the TSM reconstruction in v, depending of the ratio M k /N k .

M k N k ∈ V k , V (k) min V (k) min , V (k) max V (k) max , V k R (k) a R k-1 R k-1 R (1) k R (k) b R (2) k R k R k α k 2N k R (k) b -R k-1 2 (R k -R k-1 ) V (k) max -V (k) min V (k) max N k -M k 0 β k 0 2 (R k -R k-1 ) V (k) max -V (k) min -V (k) min N k + M k 2N k R k -R (
M k N k + -7R 3 k-1 + 27 M k N k 2 + 5R 6 k-1     1/3 +   5 27   -7R 3 k-1 + 27 M k N k - -7R 3 k-1 + 27 M k N k 2 + 5R 6 k-1     1/3 - 2R k-1 3 (E.4b)
Appendix F. Realizable, adaptive time-step algorithm

Let us consider the system of ODEs arising from solving the moment equations:

dy(t) dt = f (t, y(t)) (F.1)
where t ∈ R + , y, f ∈ R d . Explicit Runge-Kutta methods (ERK) [START_REF] Gottlieb | Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations[END_REF][START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems, 2nd Edition[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems[END_REF] are used to solve this equation: from the approximate solution y (n) at time t n , the solution at time t n+1 = t n + h is approximated by

y (n+1) = y (n) + h s i=1 b i f (τ i , η i ) (F.2)
where s is the number of stages of the ERK method, τ i and η i are given by τ i = t n + c i h (F.3a)

η i = y (n) + h i-1 j=1 a ij f (τ j , η j ) (F.3b)
and the coefficients b i , a ij , c i are given by Butcher tables [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems, 2nd Edition[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems[END_REF]].

An adaptive time-step method can be obtained with an embedded-explicit Runge-Kutta method (EERK) [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations I. Nonstiff Problems, 2nd Edition[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems[END_REF]. Then, two ERK methods are considered, using the same auxiliary values η i , but with a different linear combination of these values so that one is p th -order accurate and the other (p + 1) th -order accurate. This provides an error estimate. The method can be written in the form of a combined Butcher table: where the two last lines stand for ERK of order p and p + 1, respectively, and the solutions read

y (n+1) = y (n) + h s i=1 b i f (τ i , η i ) (F.4a) ŷ(n+1) = y (n) + h s i=1 b i f (τ i , η i ) (F.4b)
where s = s or s = s -1. In this work, we use strong stability preserving (SSP) Runge-Kutta methods [START_REF] Gottlieb | Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations[END_REF]: one of order-two embedded in the optimal third-order SSP Runge-Kutta method. The coefficients are given by where Atol i and Rtol i are, respectively, user-defined absolute and relative tolerances. As a measure of the error we take

err = 1 n d i=1 y (n) i - ŷ(n) i sc i 2 (F.7)
Then err is compared to 1 to find an "optimal" step size:

h new =

  with d ∈ {0, 1} and n ∈ N. Indeed, m N = (m 0 , m 1 , . . . , m N ) T is in the interior Int(M N ) of the moment space M N if and only if the Hankel determinants H k are positive for k ∈ {0, 1, . . . , N } [32]. Moreover, if the moment vector is on the boundary of moment space ∂M N , then the Hankel determinants H k are zero for k ≥ n. The integer n is then denoted N (m N ) and H k is positive for k < N (m N ), whereas H N (m N ) = • • • = H N = 0. In this case, the only measure corresponding to this moment vector is a sum of k weighted Dirac delta functions, with N (m N ) = 2k -1 if it is an odd number. If the considered Hankel determinants are all positive, i.e. the moment vector m N is in the interior of the moment space M N then, by convention, N (m N ) = N + 1.

Figure 1 :

 1 Figure 1: EQMOM moment-inversion algorithm.

  of moment space, then one computes the corresponding quadrature weights (w α (σ)) N α=1 and abscissas (ξ α (σ)) N α=1 and one defines m * 2N (σ) = 2N α=1 w α (σ)ξ α (σ) 2N . Let us denote m * 2N (σ) = (m * 0 , . . . , m * 2N ) T . The value of m2N (σ) is then deduced from this vector thanks to the matrix A 2N (σ) (just its last line in fact): m2N (σ) = (0, . . . , 0, 1)A 2N (σ)m * 2N (σ).

Figure 2 :

 2 Figure 2: The three possible TSM reconstructions.

3 . 4 .

 34 Compute quadrature points to obtain {w α , ξ α } N +1 α=1 from the new moment set mk 2N +1 k=0 . Update the moment set thanks to the characteristics

Figure 3 :

 3 Figure 3: Case 1. Evolution of m 0 for size reduction using v-reconstruction: exact (solid red line) and numerical (dashed blue line) solutions. Ln-EQMOM N = 4 (top left). Gamma-EQMOM N = 4 (top right). TSM with 41 sections: CFL=0.9 (bottom left) and CFL=0.05 (bottom right).

Figure 4 :

 4 Figure 4: Case 1. NDF at t = 0.033 for size reduction using v-reconstruction: exact (solid red line) and numerical (dashed blue line) solutions. Ln-EQMOM N = 4 (top left). Gamma-EQMOM N = 4 (top right). TSM with 41 sections: CFL=0.9 (bottom left) and CFL=0.05 (bottom right).

Figure 5 :

 5 Figure 5: Case 1. Evolution of m 0 for size reduction using r-reconstruction: exact (solid red line) and numerical (dashed blue line) solutions. Ln-EQMOM N = 4 (top left). Gamma-EQMOM N = 4 (top right). TSM with 9 (bottom left) and 41 sections (bottom right).

Figure 6 : 2 Figure 7 :

 627 Figure 6: Case 1. Size reduction using r-reconstruction. Abscissas of Ln-EQMOM (left) and gamma-EQMOM (right) with N = 4.

Figure 8 :

 8 Figure 8: Case 1. NDF at t = 0.033 for size reduction using r-reconstruction: exact (solid red line) and numerical (dashed blue line) solution. Ln-EQMOM N = 4 (top left). Gamma-EQMOM N = 4 (top right). TSM with 9 (bottom left) or 41 (bottom right) sections.

2 Figure 9 :

 29 Figure 9: Case 2. Error curves of moments for diffusion-controlled growth using OSM (solid red line with ×) and TSM (solid blue line with ) reconstructions. Dashed black line: order 1 and 2.

Figure 10 :

 10 Figure 10: Case 2. NDF for diffusion-controlled growth at t = 20: numerical (dashed blue line) vs. exact (solid red line). EQMOM reconstruction with r-moments, N = 4: Ln (top left), gamma (top right). OSM reconstruction (bottom left) with 90 (dot green line) and 500 (dashed blue line) sections. TSM reconstruction (bottom right) with 15 (dot green line) and 45 (dashed blue line) sections.

Figure 11 :

 11 Figure11: Case 4. Brownian aggregation with self-similar NDF. Evolution of m 0 (top left), mean particle diameter (top right) and variance (bottom) versus reference solution (red solid line). OSM (top, bottom left) with 10 (dashed magenta line with ×), 30 (dashed green line with ) and 90 (dashed blue line with •) sections. TSM (bottom right) with 5 (dashed green line with ), and 15 (dashed blue line with •) sections.

Figure 12 :

 12 Figure12: Case 4. Brownian aggregation with self-similar NDF (solid red line). Ln-EQMOM (top left): N = 4 (dashed blue line with ×). Gamma-EQMOM (top right): N = 4 (dashed blue line with ×). OSM (bottom left) with 10 (dashed magenta line with ×), 30 (dashed green line with ) and 90 (dashed blue line with •) sections. TSM (bottom right) with 5 (dashed green line with ), and 15 (dashed blue line with •) sections.

Figure 13 :

 13 Figure13: Case 5. Time evolution of m 0 (left) and mean particle diameter (right). Ln-EQMOM (top row): N = 3 (dashed green line with ), N = 4 (dashed blue line with ×). OSM (bottom row): 10 (dashed magenta line with •), 30 (dashed green line with ), 90 (dashed blue line with ×) sections. TSM, with 101 sections: (red solid line).

Figure 14 :

 14 Figure 14: Case 5. Time evolution of σ 2v . Gamma-EQMOM (top): N = 3 (dashed green line with ), N = 4 (dashed blue line with ×). OSM (bottom left): 10 (dashed magenta line with •), 30 (dashed green line with ), 90 (dashed blue line with ×) sections. TSM (bottom right): 5 (dashed green line with ), 15 (dashed magenta line with •), 45 (dashed blue line with ×), 101 (solid red line).

Figure 15 :

 15 Figure15: Case 5. Reconstructed NDF at t = 6. Ln-EQMOM (left): N = 3 (dot green line) and N = 4 (dashed blue line). Gamma-EQMOM (right): N = 3 (dot green line) and N = 4 (dashed blue line). OSM (bottom left): 10 (dashed magenta line with •), 30 (dashed green line with ), 90 (dashed blue line with ×) sections. TSM (bottom right): 5 (dashed green line with ), 15 (dashed blue line with ×), 101 (solid red line).

k) a Table E. 8 :

 8 Coefficients of the TSM reconstruction in r, depending of the ratio M k /N k .For any N k and Mk such that 0 < V k-1 N k < M k < V k N k ,the value of the coefficients are given in Table E.8, depending of the ratio M k /N k , with R

0 c 2 a 21 c 3 a 31 a 32 c

 32 s a s1 a s2 . . . a s,s

  fact that SSP Runge-Kutta yields a convex combination of some explicit Euler steps, we are able to ensure the realizability of the final moment set by checking each of the intermediate sets and adjusting the time step such that they are realizable. The parameter ierr (see Fig. F.16) is set to ierr = 0 when the corresponding intermediate moment sets are realizable.The local error estimate reads= y (n)ŷ(n) = h 3 f (τ 1 , η 1 ) + f (τ 2 , η 2 ) -2f (τ 3 , η 3 ) (F.5)We want this error to satisfy component-wise|y (n) i ŷ(n) i | ≤ sc i , sc i = Atol i + max |y

Figure F. 16 :

 16 Figure F.16: Realizable, adaptive time-step algorithm.

  )

	1								1
	0.9								0.9
	0.8								0.8
	0.7								0.7
	0.6								0.6
	0.5	0	0.02	0.04	t	0.06	0.08	0.1	0.5	0	0.02	0.04	t	0.06	0.08	0.1
	1								1
	0.9								0.9
	0.8								0.8
	0.7								0.7
	0.6								0.6
	0.5	0	0.02	0.04	t	0.06	0.08	0.1	0.5	0	0.02	0.04	t	0.06	0.08	0.1

Table 1 :

 1 Case 1. Normalized L ∞ norm in time of errors for size reduction with EQMOM using r-reconstruction.

				Ln-EQMOM		gamma-EQMOM
	N		2	3	4		2	3	4
	m 0	0.23	0.15	0.14	0.16	0.11	0.10
	m 1 m 2 m1 m0 σ 2 v NDF	1.2 × 10 -3 1.5 × 10 -4 -1.9 × 10 -5 5.8 × 10 -6 10 -4 0.13 6.8 × 10 -2 7.5 × 10 -2 8.6 × 10 -2 5.8 × 10 -2 5.1 × 10 -2 1.2 × 10 -3 1.5 × 10 -4 9.3 × 10 -5 -1.8 × 10 -5 5.5 × 10 -6 -1.3 × 10 -2 1.6 × 10 -2 -1.1 × 10 -2 9.6 × 10 -3 1.6 1.1 1.0 1.2 0.82 0.69
	VDF	0.97	0.33	0.28	0.95	0.31	0.27
				OSM				TSM
	N s		10	18	82	5	9	21	41
	m 0		0.12	0.11	2.7×10 -2 8.2×10 -2 2.8×10 -2 3.8×10 -3 5.5×10 -4
	m 1	2.4×10 -2 1.5×10 -5 2.2×10 -3 7.5×10 -3 2.8×10 -3 4.6×10 -4 4.4×10 -5
	m 2		0.40	0.38	0.34	4.4×10 -2 1.4×10 -2 1.5×10 -3 8.1×10 -5
	m1 m0	7.8×10 -2 6.5×10 -2 1.3×10 -2 4.4×10 -2 1.4×10 -2 1.9×10 -3 2.5×10 -4
	σ 2 v		1.7	1.3	1.3	0.10	3.1×10 -2 3.5×10 -3 2.2×10 -4
	NDF		1.3	1.4	0.85	1.4	1.0	0.42	0.18
	VDF		2.0	1.7	0.89	1.5	1.1	0.44	0.19

Table 2 :

 2 Case 1. Normalized L ∞ norm in time of errors for size reduction with TSM using r-reconstruction.

  × 10 -14 1.3 × 10 -13 7.0 × 10 -14 1.3 × 10 -13 1.6 × 10 -13 8.0 × 10 -14

			Ln-EQMOM			gamma-EQMOM
	N	2	3	4	2	3	4
	m 0 6.2 m 1 6.3 × 10 -4 m 2 -m1 m0 6.3 × 10 -4 σ 2 v -NDF 0.36	1.3 × 10 -4 1.4 × 10 -3 4.3 × 10 -13 3.5 × 10 -5 1.3 × 10 -4 3.5 × 10 -5 1.7 × 10 -3 2.1 × 10 -5 0.34 0.34	4.1 × 10 -4 -4.1 × 10 -4 -0.42	1.2 × 10 -4 7.6 × 10 -11 6.4 × 10 -13 4.1 × 10 -5 1.2 × 10 -4 4.1 × 10 -5 9.0 × 10 -5 2.6 × 10 -5 0.42 0.43
	VDF	0.25	0.26	0.30	0.38	0.41	0.47
							TSM are

Table 3 :

 3 Case 2. Normalized L ∞ norm in time of errors for diffusion-controlled growth with EQMOM using rreconstructions.

			OSM			TSM
	N s	10	30	90	5	15	45
	m 0 m 1 m 2 m1 m0 σ 2 v NDF 0.95 0.53 0.22 8.4 × 10 -2 2.1 × 10 -2 0.19 3.0 × 10 -2 3.6 × 10 -14 4.6 × 10 -13 1.9 × 10 -13 0.11 4.7 × 10 -3 4.7 × 10 -4 1.1 0.22 6.0 × 10 -2 0.23 1.1 × 10 -2 3.4 × 10 -4 0.87 0.17 2.6 × 10 -2 0.11 4.7 × 10 -3 4.7 × 10 -4 4.7 0.68 0.14 0.54 1.2 × 10 -2 1.4 × 10 -4 0.75 0.42 1.2 0.54 0.19
	VDF 1.1	0.53	0.19	0.67	0.29	6.9 × 10 -2

Table 4 :

 4 Case 2. Normalized L ∞ norm in time of errors for diffusion-controlled growth with TSM using rreconstructions.

  × 10 -7 1.1 × 10 -7 3.3 × 10 -8 4.7 × 10 -7 1.2 × 10 -7 3.3 × 10 -8

			-EQMOM			gamma-EQMOM
	N	2	3	4	2	3	4
	m 0 4.5 m 2 m1 m0 4.5 × 10 -7 1.1 × 10 -7 3.3 × 10 -8 4.7 × 10 -7 1.2 × 10 -7 3.3 × 10 -8 -1.8 × 10 -6 5.1 × 10 -7 -1.9 × 10 -6 5.2 × 10 -7 σ 2 v -1.7 × 10 -6 5.0 × 10 -7 -1.8 × 10 -6 5.0 × 10 -7 NDF 0.22 0.36 0.66 0.25 0.18 0.12
	VDF	0.17	0.20	0.27	0.36	0.27	0.20

Table 5 :

 5 Case 3. Aggregation with the sum kernel: L ∞ norm in time of relative errors with EQMOM using rreconstruction.

			OSM			TSM
	N s	10	30	90	5	15	45
	m 0 m 2 m1 m0 σ 2 v NDF VDF	0.25 1.5 × 10 3 0.20 1.4 × 10 3 0.55 0.94	0.13 5.3 × 10 -2 1.7 × 10 -3 4.5 × 10 -3 6.3 × 10 -3 31 0.26 0.73 0.13 2.4 × 10 -3 0.12 5.0 × 10 -2 1.7 × 10 -3 4.5 × 10 -3 6.3 × 10 -3 29 0.22 0.75 0.13 3.7 × 10 -3 0.22 8.4 × 10 -2 0.32 2.7 × 10 -2 1.5 × 10 -2 0.37 8.3 × 10 -2 0.79 9.1 × 10 -2 1.1 × 10 -2

Table 6 :

 6 Case 3. Aggregation with the sum kernel: L ∞ norm in time of relative errors with OSM and TSM using r-reconstruction.

It is the truncated Hausdorff moment problem if the internal variable lives in a compact support, and the Hamburger problem if it lives on the real line.

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).

conditioning the other internal variable by the volume, using the mono-kinetic assumption, along the lines of what has been done with the size-velocity internal variables for sprays [START_REF] Laurent | Multi-fluid modelling of laminar polydisperse spray flames: origin, assumptions and comparison of sectional and sampling methods[END_REF][START_REF] Laurent | Eulerian multi-fluid modeling for the numerical simulation of coalescence in polydisperse dense liquid sprays[END_REF][START_REF] Vié | Size-velocity correlations in hybrid high order moment/multifluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF][START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF] and bubbly flows [START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF][START_REF] Buffo | Multivariate quadrature-based moment methods for turbulent polydisperse gas-liquid systems[END_REF][START_REF] Yuan | An extended quadrature-based mass-velocity moment model for polydisperse bubbly flows[END_REF]. Another very important topic for fine-particle applications involving spatial transport of the moments is the development of realizable, high-order advection [START_REF] Vikas | Realizable high-order finite-volume schemes for quadrature-based moment methods[END_REF][START_REF] Kah | A high order moment method simulating evaporation and advection of a polydisperse liquid spray[END_REF][START_REF] Laurent | Realizable high-order finite-volume schemes for the advection of moment sets of the particle size distribution[END_REF][START_REF] Essadki | Adaptive mesh refinement for polydisperse spray simulation[END_REF] and diffusion [START_REF] Vikas | Realizable high-order finite-volume schemes for quadraturebased moment methods applied to diffusion population balance equations[END_REF] schemes, which can be coupled with a Navier-Stokes solver for the carrier fluid. Ideally, such schemes should accommodate unstructured computational grids and both explicit and implicit time stepping.

Appendix A. Change of variable in the PBE

In this work, the chosen size variable is volume v, which is then the internal variable of the NDF: f (t, v). However, for spherical particles, the radius r = 3v 4π 1/3 can be equivalently used. Since the corresponding NDF f r represents the same particle population, one has f (t, v)dv = f r (t, r)dr, and thus

While the change of variable is easy for most of the terms in the PBE, the drift terms are written Only continuous size reduction is considered here, with a rate given by R red (v) = -c red v 2/3 1 R+ (v), with c red = 50 The initial NDF is the sum of two log-normal distributions:

where

The initial NDF is bimodal, as is typical of the ones induced by nucleation and aggregation. Since c red is constant, the analytical solution, using the volume as the internal variable, is given by

. Case 2: diffusion-controlled growth test case

The test case of pure diffusion-controlled growth described by McGraw [START_REF] Mcgraw | Description of aerosol dynamics by the quadrature method of moments[END_REF] is considered here. The initial distribution is f r (0, r) = ar 2 exp(-b r) where r is the particle radius (µm), a = 0.108 µm -3 cm -3 and b = 0.6 µm -1 . The growth rate is given by R gro (t, v) = c gro v 1/3 , with c gro = 2.34 µm 2 /s (considering v = r 3 ).

Since the coefficient c gro is constant, the analytical solution, using the volume as the internal variable is given by:

since the characteristics are given by

Using the variable r = v 1/3 for ξ, we find

Appendix B.3. Case 3: aggregation test case with the sum kernel Pure aggregation test case is considered here, using the sum kernel β(v, v ) = v + v and with initial condition f (0, v) = exp(-v). An analytical NDF for this problem is given in [START_REF] Gelbard | Numerical solution of the dynamic equation for particulate systems[END_REF].

Appendix B.4. Case 4: aggregation test case with the Brownian kernel

Pure aggregation test case is considered here, using the more complex Brownian kernel: v) is used as for case 3. In this case, a self-similar NDF is found for large enough time [START_REF] Vemury | Self-preserving size distributions of agglomerates[END_REF].

Appendix B.5. Case 5: test case with nucleation, aggregation and size reduction

We investigate here a test case with nucleation, aggregation and size reduction where the nuclei size is V nuc = 1 and the nucleation rate is j(t) = 10 3 1 [0,0.5] (t). The Brownian kernel is used for aggregation, and the size-reduction rate is given by R red (v) = -c red v 2/3 1 R+ (v) with c red = 3.

Appendix C. ζ-Chebyshev algorithm

In order to compute the coefficients {ζ k } k∈N * from the moments {m k } k∈N , the Chebyshev algorithm is used, combined to formula [START_REF] Zucca | Implementation of the population balance equation in CFD codes for modelling soot formation in turbulent flames[END_REF]. Then, setting σ 0,k = m k for k ≥ 0, a 0 = m 1 /m 0 and ζ 1 = a 0 , one computes, for l ≥ 1: σ 1,l = σ 0,l+1a 0 σ 0,l .

For k ≥ 1, one then computes successively:

Appendix D. KDF for EQMOM

Two types of KDF are used for EQMOM reconstruction, when considering the Stieljes problem: gamma [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF] and log-normal [START_REF] Madadi-Kandjani | An extended quadrature-based moment method with log-normal kernel density functions[END_REF]. Let us recall here their corresponding matrix A k (σ), as well as the interval on which σ is located.

Appendix D.1. Gamma-EQMOM

The gamma KDF reads

where Γ(•) is the gamma function, σ > 0, k α > 0 and k α = ξ α /σ. The moment of order k of the gamma KDF is given by

where a k,k = 1 and a ki is given by the following recurrence formula:

Thus, A k (σ) is a triangular matrix of coefficients [A k (σ)] i,j = a ij σ i-j for 0 ≤ j ≤ i, which can be easily inverted. Moreover, the coefficients a ij can be precomputed and stored. Finally, the maximum value for σ, in the case N = 2 is given by

. Ln-EQMOM

The log-normal (Ln) KDF reads

where σ and ξ α are positive. The moment of order k of the Ln KDF is ξ k,α = ξ k α e k 2 σ 2 /2 . The matrix A k (σ) is a diagonal matrix of coefficients [A k (σ)] i,i = e i 2 σ 2 /2 so that the relation between m k and m * k is easy to compute. Moreover, a change of variable is used in [START_REF] Kumar | Improved accuracy and convergence of discretized population balance for aggregation: The cell average technique[END_REF]. The nonlinear problem is then solved for b = e σ 2 /2 , in the interval (1, b

Appendix E. TSM method

Let us give here the formula of the TSM reconstruction in the two considered cases: one using the volume variable v and one using the variable r = v 1/3 .

Appendix E.1. v-reconstruction

The NDF is approximated by an affine function on each section

] (v). (E.1)

Let us define the two following mean volumes:

For any N k and M k such that 0 < V k-1 N k < M k < V k N k , the value of the coefficients are given in Table E.7, depending of the ratio M k /N k [START_REF] Laurent | Two-size moment multi-fluid model: a robust and highfidelity description of polydisperse moderately dense evaporating sprays[END_REF].

Appendix E.2. r-reconstruction When considering the r variable, the NDF is approximated by the following function on each section I k :

b ] (r)

Let us define the following two mean volumes:

10 (E.3b)