Alexandre Bouaniche

Jerome Yon

Pascale Domingo

Luc Vervisch
email: vervisch@coria.fr

Analysis of the soot particle size distribution in a laminar premixed flame: A hybrid stochastic/fixed-sectional approach

Keywords: Soot, Stochastic method, Sectional method, Probability density function, Particle size distribution, Shape modeling

The numerical simulation of the soot particle size distribution in flames is addressed by solving the balance equations for total number density and the probability density function (PDF) of particle sizes, leading to a hybrid stochastic/fixed-sectional method for solving soot population balance equation. Well established models are introduced in these equations for the chemistry of PAH in ethylene combustion, for particle nucleation, growth, agglomeration and oxidation. These closures are combined with a stochastic approach, which drives the evolution of a fixed number of computational particles used to solve for the particle size distribution with a control of agglomeration and numerical roundoff error through a fixed sectional discretisation. A laminar sooting flame is simulated to compare the results against measurements and previous numerical simulations, confirming the validity of the novel approach in terms of accuracy and CPU efficiency. The relation between the mobility diameter, measured in the experiments, and the equivalent sphere diameter, introduced in the modeling, is discussed under this novel numerical framework. The influence of the fractal particle shape on the simulated particle size distribution is explored. Finally, Particle Size Distributions obtained from the hybrid method are compared to the ones obtained with a representative fixed-sectional method.

Introduction

The prediction of soot particles emission motivates intense research, with direct application for human health and environment [START_REF] Donaldson | Combustion-derived nanoparticles: A review of their toxicology following inhalation exposure[END_REF][START_REF] Jensen | The potential impact of soot particles from aircraft exhaust on cirrus clouds[END_REF]. The development of novel techniques of synthesis of nano-structured materials also benefits from these works [START_REF] Kelesidis | Flame synthesis of functional nanostructured materials and devices: Surface growth and aggregation[END_REF]. In terms of numerical modeling, this is a very challenging field, in which complex hydrocarbon chemistry is combined with the description of a reacting solid phase. Usual flow transport phenomena need to be coupled with mechanisms of nucleation, growth, agglomeration and oxidation of solid particles, which are driven by the intricate coupling existing between particle dynamics, gaseous and surface chemistry.

Models based on mono-disperse assumption have been widely used for their relative simplicity. They are usually cast in two-equation approaches, solving for total number density of particles and soot volume or mass fraction [START_REF] Ma | Optimization of soot modeling in turbulent nonpremixed ethylene/air jet flames[END_REF]. More recently, advanced methods based on moments of the particle size distribution have been discussed in the literature with success [START_REF] Frenklach | Method of moments with interpolative closure[END_REF][START_REF] Lindstedt | Joint-scalar transported PDF modeling of soot formation and oxidation[END_REF][START_REF] Mueller | Hybrid method of moments for modeling soot formation and growth[END_REF][START_REF] Salenbauch | A numerically robust method of moments with number density function reconstruction and its application to soot formation, growth and oxidation[END_REF][START_REF] Wick | Modeling soot oxidation with the extended quadrature method of moments[END_REF][START_REF] Zucca | Implementation of the population balance equation in cfd codes for modelling soot formation in turbulent flames[END_REF]. In some applications, the knowledge of the particle size distribution (PSD) is mandatory and this information cannot always be reconstructed accurately from only a few moments. Along these lines, sectional methods provide a direct access to detailed information on the PSD, through a discretisation of the range of particle sizes. These methods have been widely used in the field of soot modeling [START_REF] Aubagnac-Karkar | Sectional soot model coupled to tabulated chemistry for Diesel RANS simulations[END_REF][START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF]. Fixed-sectional methods featuring a fixed discretisation of particle size were introduced [START_REF] Gelbard | Numerical solution of the dynamic equation for particulate systems[END_REF][START_REF] Hounslow | A discretized population balance for nucleation, growth, and aggregation[END_REF]. Subsequently, significant effort was dedicated to reducing the impact of numerical diffusion linked to growth in the space of particle size. Different strategies were proposed like the use of higher order schemes [START_REF] Gunawan | High resolution algorithms for multidimensional population balance equations[END_REF][START_REF] Ma | High-resolution simulation of multidimensional crystal growth[END_REF][START_REF] Park | A novel fixed-sectional model for the formation and growth of aerosol agglomerates[END_REF] or moving/adaptive grids [START_REF] Kumar | On the solution of population balance equations by discretization nucleation, growth and aggregation of particles[END_REF][START_REF] Sewerin | An explicit adaptive grid approach for the numerical solution of the population balance equation[END_REF].

Other approaches have been based on direct simulation Monte Carlo (DSMC) techniques, in which the physical population of gas particles is represented by a collection of computational particles. The particle dynamics (nucleation, growth, coagulation, flow transport) are represented at the level of these computational particles in a systematic way, allowing for the number of computational particles to increase or decrease according to soot evolution. DSMC was successfully applied to particle nucleation, agglomeration, and surface growth in [START_REF] Kruis | Direct simulation monte carlo method for particle coagulation and aggregation[END_REF][START_REF] Liffman | A direct simulation monte-carlo method for cluster coagulation[END_REF][START_REF] Maisels | Direct simulation monte carlo for simultaneous nucleation, coagulation, and surface growth in dispersed systems[END_REF] and applied to soot modeling in laboratory flames in [START_REF] Balthasar | A stochastic approach to calculate the particle size distribution function of soot particles in laminar premixed flames[END_REF]. Several techniques have been proposed to enhance Monte Carlo methods for predicting the Particle Size Distribution and lower associated computational cost. Constant number Monte Carlo [START_REF] Khalili | Constant number monte carlo simulation of population balances with multiple growth mechanisms[END_REF][START_REF] Lin | Solution of the population balance equation using constant-number monte carlo[END_REF][START_REF] Smith | Constant-number monte carlo simulation of population balances[END_REF], Event-Driven Constant-Volume Monte Carlo [START_REF] Zhao | A new event-driven constant-volume method for solution of the time evolution of particle size distribution[END_REF], majorants and binary trees [START_REF] Patterson | Extending stochastic soot simulation to higher pressures[END_REF] and stochastic weighted methods [START_REF] Deville | Weighted flow algorithms (wfa) for stochastic particle coagulation[END_REF][START_REF] Eibeck | Stochastic particle approximations for Smoluchoski coagulation equation[END_REF][START_REF] Haibo | Multi-monte carlo approach for general dynamic equation considering simultaneous particle coagulation and breakage[END_REF][START_REF] Patterson | Stochastic weighted particle methods for population balance equations[END_REF].

Direct Monte Carlo methods are accurate, but still require quite large computing resources. Even though computational cost was reduced, it might not be in line with unsteady three-dimensional simulations of real combustion chambers.

An attempt is made in this work to tackle the problem through the soot size probability density function (PDF). This approach features a fixed number of numerical particles to build the PDF at the location (x, t) and it is organised to facilitate its coupling with large eddy simulation (LES). This novel approach, combining stochastic Monte Carlo (numerical particles) and fixed-sectional methods, has been found to minimise the discretisation errors when solving the surface growth/loss term of the population balance equation for reference canonical problems [START_REF] Bouaniche | A hybrid stochastic/fixed-sectional method for solving the population balance equation[END_REF]. The solution method relies on a fixed number of stochastic particles and sections, with a numerical algorithm organised to minimise errors even for a moderate number of stochastic particles and sections. Setting a fixed number of stochastic particles per grid cell brings many advantages on the computing side, but also poses specific challenges to simulate nucleation and agglomeration. These points are discussed and the algorithm proposed is applied to the prediction of the soot size distribution in a one-dimensional freely propagating premixed flame, previously investigated both experimentally and numerically by Zhao et al. [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF]. It is important to note that the focus of the present paper is not to discuss new physical models for soot particle nucleation, agglomeration or surface growth or loss. The main objective is to propose an efficient and accurate numerical scheme for solving the Population Balance Equation in the context of sooting flames. The performance of this method has been previously assessed against a series of analytical test cases in [START_REF] Bouaniche | A hybrid stochastic/fixed-sectional method for solving the population balance equation[END_REF]. In the present work, the applicability of the hybrid method to sooting flames is demonstrated. It is first validated using the physical models for soot source terms as in the Monte Carlo simulation of [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF]. Then, the presented hybrid method is shown to be more efficient than a reference fixed-sectional method. Finally, the influence of aggregate shape modelling on the Particle Size Distribution is studied.

2 The hybrid stochastic/fixed-sectional method for soot population balance equations

Statistical description

Let us define the particle size distribution (PSD) n(v; x, t), number of particles of characteristic size v (in terms of volume or mass), per unit of flow volume and per unit of characteristic size; N i (x, t), the number of soot particles of characteristic size v i per unit of flow volume and N T (x, t), the total number of soot particles per unit of flow volume. These quantities are related according to

N i (x, t) = Iv i n(v; x, t)dv , (1)
N T (x, t) = ∞ vo n(v; x, t)dv = M -1 i=0 N i (x, t) , (2)
where the interval

I vi ≡ [v inf i , v sup i ] defines the i-th fixed-section of size, for i = 0, • • • , M -1, with M the number of fixed sections and v inf 0 is denoted v o for brevity.
Let us consider N T (x, t), the total number of particles per unit volume, and P (v * ; x, t), the probability density function (PDF) of the particles characteristic size, where v * ∈ [v o , ∞] denotes the sample space variable associated to v, which is seen as a random variable. The motivation for using this PDF is to be able to describe the surface growth or loss term as a linear term instead of directly solving for the Population Balance Equation, in which the growth/loss term appears in a strongly non-linear form.

The relation between n(v; x, t), the particle number density per unit size, N i (x, t), the number density of particles whose size is in the section I vi (v ∈ I vi ) at the flow position 'x' at time 't' (Eq. ( 1)), and P (v * ; x, t), the PDF of the particles sizes reads:

Iv i n(v * ; x, t)dv * = N i (x, t) = N T (x, t) Iv i P (v * ; x, t)dv * , (3)
where

Iv i P (v * ; x, t)dv * (4)
is the probability to find particles of sizes v ∈ I vi . Because (3) should be valid whatever I vi , in particular in the limit where

dv i = v sup i -v inf i
goes to zero:

n(v * ; x, t) = N T (x, t)P (v * ; x, t) . (5)
The function

δ(v -v * ) = lim dv→0 1/dv if v ∈ [v * -dv/2, v * + dv/2] (6)
= 0 otherwise ,

is introduced and

P (v * ; x, t) = δ(v(x, t) -v * )
, where • denotes a statistical average [START_REF] Dopazo | Relaxation of initial probability density functions in the turbulent convection of scalar fields[END_REF][START_REF] Dopazo | Statistical description of the turbulent mixing of scalar fields[END_REF][START_REF] Kollmann | The pdf approach to turbulent flow[END_REF][START_REF] Lundgren | Distribution function in the statistical theory of turbulence[END_REF]. The PDF P (v * ; x, t) of soot sizes is discretised over a set of N P stochastic particles, each carrying information on the particle size, i.e., v = v k for k = 1, • • • , N P and P (v * ; x, t) = (1/N P )

N P k=1 δ(v k (x, t) -v * ).
The representative size v k of a stochastic particle can take any value between the considered size bounds [v o , v M ], whatever the fixed sectional mesh, and N P the total number of stochastic particles is fixed. The mesh is divided into M sections, each section being defined by the interval I vi of length ∆v i = v sup i -v inf i . Within the set of N P particles, an integer number n Pi (x, t) of stochastic particles have sizes so that

v k ∈ I vi ≡ [v inf i , v sup i
]. This number of stochastic particles relates to the PDF,

Iv i P (v * ; x, t)dv * = n Pi (x, t) N P Stochastics = Physical particles N i (x, t) N T (x, t) . (8)
Therefore, simulating the evolution of n Pi (x, t) and N T (x, t) allows for fully characterising the Particle Size Distribution.

Balance equations

The population balance equation followed by n(v * ; x, t) may be written in the limit of dv * → 0 [START_REF] Bouaniche | A hybrid stochastic/fixed-sectional method for solving the population balance equation[END_REF][START_REF] Ramkrishna | The status of population balances[END_REF][START_REF] Ramkrishna | Population balances. Theory and applications to particulate systems in engineering[END_REF][START_REF] Solsvik | The foundation of the population balance equation: A review[END_REF],

∂n(v * ; x, t) ∂t + u • ∇n(v * ; x, t) + ∂ ∂v [G(v * )n(v * ; x, t)] (9) = Ḣo (x, t)δ(v o -v * ) + Ȧi (x, t)δ(v i -v * ) ,
where G is the soot growth/oxidation rate. The nucleation source per unit of flow volume Ḣo (x, t) and the Smoluchowski agglomeration source/sink for the i-th section Ȧi (x, t) are defined as follows from the nucleation ḣ(v; x, t) and agglomeration ȧ(v; x, t) sources per unit of flow volume and per unit of characteristic size:

Ḣo (x, t) = Iv o ḣ(v; x, t)dv . ( 10
) Ȧi (x, t) = Iv i ȧ(v; x, t)dv , (11) with ȧ
(v; x, t) (12) = 1 2 v 0 β(v -v, v)n(v -v; x, t)n(v; x, t)dv -n(v; x, t) ∞ 0 β(v, v)n(v; x, t)dv ,
and A T is the total sink due to agglomeration over all particles, thus the sum of A i (x, t) over all sections

ȦT (x, t) = ∞ vo ȧ(v; x, t)dv = M -1 i=0 Ȧi (x, t) . (13)
The transport equation for the total number particle density N T reads

∂N T (x, t) ∂t + u(x, t) • ∇N T (x, t) = Ḣo (x, t) + ȦT (x, t) , (14)
Combining the relation ( 8) with ( 9) and ( 14) provides the balance equation for the PDF [START_REF] Bouaniche | A hybrid stochastic/fixed-sectional method for solving the population balance equation[END_REF],

∂P (v * ; x, t) ∂t + u(x, t) • ∇P (v * ; x, t) = (i) - ∂ ∂v * G(v * )P (v * ; x, t) + Ḣo (x, t) N T (x, t) δ(v o -v * ) -P (v * ; x, t) (ii) + 1 N T (x, t) Ȧi (x, t)δ(v i -v * ) -ȦT (x, t)P (v * ; x, t) (iii) , (15)
with (i) the surface growth/loss term transporting the probability in size space.

(ii) the nucleation rate which increases the probability to find particles at the smallest size considered and thereby decreases the probability density of other sizes and (iii) the agglomeration also redistributing the probability according to the agglomeration rates.

Stochastic/fixed-sectional algorithm

The stochastic particles are advanced in time following the PDF equation [START_REF] Gelbard | Numerical solution of the dynamic equation for particulate systems[END_REF]. Surface growth or loss is applied to get v k (t n+ 1 2 ),

dv k (t) dt = G(v k (t)) , k = 1, • • • , N P . (16)
The number density evolves to N T (x, t n+1 ) with

dN T (t) dt = Ḣo (t) + ȦT (t) . (17)
The PDF equation ( 15) is also advanced in time with nucleation and agglomeration:

P (v * ; t n+1 ) = α Ho δ(v o -v * ) + α Ai δ(v -v * ) + (1 -α Ho -α A T )P (v * ; t n+ 1 2 ) . (18)
From Eq. ( 15):

α Ho = Ḣo (t n+ 1 2 ) N T (t n+1 ) • δt , (19)
α Ai = Ȧi (t n+ 1 2 ) N T (t n+1 ) • δt , (20)
α A T = ȦT (t n+ 1 2 ) N T (t n+1 ) • δt . (21)
Multiplying the PDF evolution by N P × ∆v i returns the variation of the number of stochastic particles per section (see Eq. ( 8)),

n Pi (t n+1 ) = n Pi (t n+ 1 2 ) + ∆n Pi (t n+ 1 2 ) , (22)
n P is an integer and ∆n P a real number with

∆n Po (t n+ 1 2 ) = (α Ho + α Ao ) N P -(α Ho + α A T )n Po (t n+ 1 2 ) + α Ro N P , (23)
∆n Pi (t

n+ 1 2 ) = α Ai N P -(α Ho + α A T )n Pi (t n+ 1 2 ) + α Ri N P for i = o . ( 24
)
α Ri represents the accumulation of the round-off error, vanishing for N P → ∞. The fractional part {∆n Pi (t n+ 1 2 )} is defined from the nearest integer ∆n Pi (t n+ 1 2 ) ,

{∆n Pi (t n+ 1 2 )} = ∆n Pi (t n+ 1 2 ) -∆n Pi (t n+ 1 2 ) . (25)
The following Monte Carlo algorithm is applied:

-For ∆n Pi (t n+ 1 2 ) < 0, a number -∆n Pi (t n+ 1 2 ) of stochastic particles is randomly chosen among the n Pi (t n+ 1 2 ) ∈ I vi . -The selected particles from all I vi intervals (i = 0, • • • , M -1) form an ensemble P(t n+ 1 2 ) of particles whose characteristic size must change. -For ∆n Pi (t n+ 1 2 ) > 0, a number ∆n Pi (t n+ 1
2 ) of particles are randomly taken from P(t n+ 1

2 ) and allocated to I vi around the representative size v i (t n+ 1 2 ), defined to conserve mass, as discussed thereafter.

The residual decimal part {∆n Pi (t n+ 1 2 )} defines N R i (t n ), a residual number density of physical particles in the section I vi , which is computed at time t n according to [START_REF] Donaldson | Combustion-derived nanoparticles: A review of their toxicology following inhalation exposure[END_REF],

N R i (t n ) = {∆n Pi (t n-1 2 )} N P N T (t n ) , (26)
∆n Pi (t n-1 2 ) denotes ∆n P of the previous iteration in time. The growth/loss of the physical particles represented by this number density residual N R i (t n ) is solved separately between t n and t n+ 1 2 . This is done with a sectional method based on the 3-point discretisation for particle growth/loss [START_REF] Park | A novel fixed-sectional model for the formation and growth of aerosol agglomerates[END_REF]. α Ri is then obtained from

α Ri = N R i (t n+ 1 2 ) N T (t n+1 ) , (27)
to compute ∆n Pi (t n+ 1 2 ) from Eqs. ( 23) and ( 24). The agglomeration source Ȧi (t) of Eq. ( 20) is computed from the method proposed in [START_REF] Kumar | On the solution of population balance equations by discretization -I: A fixed pivot technique[END_REF]. For any colliding particles of volume v in section i and v in section j, the collision kernel

β(v, v) is assumed fixed to β(v i , v j ) = β i,j .
Particles formed by agglomeration are distributed in the sections so that the zeroth and first moments (number and mass) of the PSD are conserved. This method is computationally efficient by evaluating the double integrals of the collision kernel at the aggregated level of M sections. This results in a M × M calculation for the agglomeration source terms in equation [START_REF] Kruis | Direct simulation monte carlo method for particle coagulation and aggregation[END_REF]. It is to be noted that if every stochastic particles pairs were considered this would result in a more expensive N P × N P calculation. Methods as majorants and binary trees [START_REF] Patterson | Extending stochastic soot simulation to higher pressures[END_REF] have been developed for Monte Carlo methods to avoid the cost of a full N P × N P calculation. Our strategy in the proposed method is different as we calculate the agglomeration source terms just as in a sectional method [START_REF] Kumar | On the solution of population balance equations by discretization -I: A fixed pivot technique[END_REF]. The agglomeration source used in [START_REF] Hounslow | A discretized population balance for nucleation, growth, and aggregation[END_REF] reads

Ȧi (t) = k≤j≤i j,k vi-1≤vj +v k ≤vi+1 1 - δ j,k 2 ηβ j,k N j (t)N k (t) -N i (t) M -1 k=0 β i,k N k (t) , (28)
with

η =        v i+1 -(v j + v k ) v i+1 -v i if v i ≤ v j + v k ≤ v i+1 , v i-1 -(v j + v k ) v i-1 -v i if v i-1 ≤ v j + v k ≤ v i , (29)
The characteristic volume v i must be representative of the average mass contained in the i-th section. v i is computed from the volumes of the stochastic particles contained in both the section and the residual terms resulting from the roundoff,

v i (t n+ 1 2 ) = (N T (t n )/N P ) n P i (t) k=1 v k i (t n+ 1 2 ) + N R i (t n+ 1 2 )v i (t n ) (N T (t n )/N P )n Pi (t n+ 1 2 ) + N R i (t n+ 1 2 ) , (30)
where

v k i = v k if v k ∈ I vi and v k i = 0 otherwise, n Pi (t n+ 1 2
) is the number of stochastic particles in the i-th section (Eq. ( 8)) and N R i (t n+ 1 2 ) is the residual number density of the particles in the section after applying surface gross or loss. v i needs to be updated again after reallocation of the particles due to agglomeration, to provide v i (t n+1 ) from [START_REF] Kumar | On the solution of population balance equations by discretization nucleation, growth and aggregation of particles[END_REF] with

N T (t n+1 ), v k i (t n+1 ), N R i (t n+1 ), v i (t n+ 1 2 ), n Pi (t n+1
). Once the agglomeration source is known from Eq. ( 28), in the sections for which ∆n Pi (t n+ 1 2 ) > 0, the stochastic particles to be reassigned are chosen according to the Monte Carlo procedure discussed above. In practice, this is done in two steps. First, the ∆n Pi (t n+ 1 2 ) particles are allocated randomly within the section at sizes v k (t n+ 3 4 ), which are samples of a random variable v following a target piecewise linear distribution [21] defined by the probability density function,

p(v | v inf i , v sup i , w i , w i+1 ) = 2 w i (v sup i -v) + w i+1 (v -v inf i ) (w i + w i+1 )∆v 2 i . (31)
In this distribution, the weights, w i , are calculated from the variations of the number densities at v i ,

w i = ∆n(v i-1 ; t n+ 1 2 ) + ∆n(v i ; t n+ 1 2 ) -∆n(v i-1 ; t n+ 1 2 ) v i -v i-1 (v inf i -v i-1 ) , w i+1 = ∆n(v i ; t n+ 1 2 ) + ∆n(v i+1 ; t n+ 1 2 ) -∆n(v i ; t n+ 1 2 ) v i+1 -v i (v sup i -v i ) ,
with ∆n i (t) = ∆n Pi (t)N T (t)/(N P ∆v i ) (Eq. ( 8)). In the process of applying the sources to the PDF evolution, such random piecewise linear distribution secures a continuous distribution of the stochastic particles. However, it does not guarantee strict volume/mass conservation by itself. This is achieved in a second step by calculating a corrective factor K i

K i = v i (t n+ 1 2 ) (1/n Pi (t n+ 3 4 )
)

n P i (t n+1 ) k=1 v k i (t n+ 3 4 ) , (32)
then,

v k i (t n+1 ) = K i v k i (t n+ 3 4 ) , ( 33
)
and mass is conserved through the reallocation process. Nucleation size is set as the lower boundary of the smallest size section v o . As the numerical steps corresponding to nucleation/agglomeration and growth are sequential in the present model, it is necessary to account for a dispersion of effective nucleation sizes due to particle growth during the nucleation/agglomeration time step. For ∆n Po (t n+ 1 2 ) > 0, the ∆n Po (t n+ 1 2 ) particles are therefore allocated randomly following a target uniform distribution between v o and v o + G(v o )δt.

A systematic analysis of this hybrid stochastic/fixed-sectional method for solving the population balance equation, along with comparisons against analytical solutions of canonical PSD evolutions may be found in [START_REF] Bouaniche | A hybrid stochastic/fixed-sectional method for solving the population balance equation[END_REF].

Modeling of nucleation, growth and agglomeration source terms

Each stochastic particle is associated with a volume v k . Soot particles are considered spherical (usual assumption in the range of sizes studied [START_REF] Aubagnac-Karkar | Sectional soot model coupled to tabulated chemistry for Diesel RANS simulations[END_REF][START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF]), except for the results presented in section 5.

Well established nucleation modeling by the collision of two pyrene (C 16

H 10 ) is used [1], Ḣ(x, t) = 0.5β py N 2 py , (34)
with

N py = [C 16 H 10 ]N A , (35)
the volume number of pyrene molecules, N A is the Avogadro constant.

The surface growth rate G is controlled by (i) condensation of pyrene molecules on soot particles (G Cond ), (ii) C 2 H 2 (acetylene) addition by the HACA mechanism and surface oxidation by O 2 and OH (G HACA,Oxi ) [START_REF] Appel | Kinetic modeling of soot formation with detailed chemistry and physics: laminar premixed flames of C2 hydrocarbons[END_REF]. [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF]. The condensation source term is:

G(v k ; x, t) = G Cond (v k ; x, t) + G HACA,Oxi (v k ; x, t) is applied to each particle v k through equation
G Cond (v k ; x, t) = m py ḢCond (v k ; x, t)/N i (x, t) , (36)
(for v k ∈ I vi ), with m py the mass of one pyrene molecule and

ḢCond (v k ; x, t) = β v k ,py N i (x, t)N py . (37)
The surface reaction source term is

G HACA,Oxi (v k ; x, t) = ( ωC2H2 + ωO2 + ωOH )/N i (x, t) , (38)
with the chemical sources,

ωC2H2 = 2W C k C2H2 [R v k ][C 2 H 2 ] , (39) ωO2
= -2W C k O2 [R v k ][O 2 ] , (40) ωOH
= -W C k OH [S v k ][OH] , (41)
where W C is the molar mass of Carbon, and k C2H2 , k O2 , k OH are calculated using the kinetic parameters given in [START_REF] Appel | Kinetic modeling of soot formation with detailed chemistry and physics: laminar premixed flames of C2 hydrocarbons[END_REF].

[S v k ] and [R v k ], are obtained as in [1, 2] through [S v k ] + [R v k ] = s k χ α HACA N i (x, t)/N A , (42)
s k is the particle surface and χ is the number of sites per unit surface of soot.

[S v k ] is representative of stable sites at soot particles surface and [R v k ] of radical sites. [R v k ] is calculated assuming that radical sites are in quasi-steady state. α HACA (v k ) is adapted from [START_REF] Appel | Kinetic modeling of soot formation with detailed chemistry and physics: laminar premixed flames of C2 hydrocarbons[END_REF], calculating its value for each stochastic particle instead of using the first size moment to calculate a unique value of α HACA . In practice

α HACA (v k ) = tanh a log(ρ soot × v k /(W C /N A ) + b , (43)
where a = 12.65 -0.00563T and b = -1.38 + 0.00068T . This dependence of α HACA (v k ) on particle size expresses the fact that mature particles feature a lower proportion of active sites per unit of surface. Agglomeration is described by the standard Smoluchowski equation. The collision rates entering the Smoluchowski equation expressing Ȧi (x, t), are calculated depending on the Knudsen number as in [START_REF] Blanquart | A joint volume-surface-hydrogen multi-variate model for soot formation[END_REF][START_REF] Kazakov | Dynamic modeling of soot particle coagulation and aggregation: Implementation with the method of moments and application to high-pressure laminar premixed flames[END_REF], with β f,m vi,vj in the free molecular regime and β cont vi,vj in the continuum regime. The harmonic mean of the asymptotic values is taken in the transition regime β f,m vi,vj β cont vi,vj /(β f,m vi,vj + β cont vi,vj ). Collision frequencies are calculated in the same manner for collisions between pyrene molecules (nucleation) or between pyrene molecules and soot particles (condensation). The Smoluchowki equation is discretised as in [START_REF] Aubagnac-Karkar | Sectional soot model coupled to tabulated chemistry for Diesel RANS simulations[END_REF][START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF] to account for the source term Ȧi (x, t) in each I vi interval (or section). 4 Simulation of a one-dimensional freely propagating laminar premixed flame

Flow configuration and model parameters

A one-dimensional fuel-rich (φ = 2.07) laminar premixed ethylene-argonoxygen flame [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF] is simulated at ambient pressure. The velocity, temperature and mole fractions in the fresh gases are 8.26 cm/s, 300 K, X C2H4 = 0.133, X O2 = 0.193 and X Ar = 0.674, respectively. Soot particles whose mobility diameter is superior to 3 nm were detected through SMPS, to provide PSD at several heights above burner. Both measurements and simulations reported in [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF] will serve as reference to evaluate the proposed soot modeling. The gaseous one-dimensional flame is first computed with complex molecular transport and the detailed chemical scheme by Appel et al. [START_REF] Appel | Kinetic modeling of soot formation with detailed chemistry and physics: laminar premixed flames of C2 hydrocarbons[END_REF] (101 species and 544 elementary reactions) using CANTERA [START_REF] Goodwin | Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF]. To limit potential compensation of errors between heat transfer and chemistry, as in previous works [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF], the measured temperature profile is imposed.

The gas phase species mass fractions (Figures 1(a) and 1(b)) agree with previous calculations, in which soot chemistry was fully coupled. In our case, only pyrene consumption is considered when solving Eqs. ( 14)-( 15) for soot, without much impact on most influential species. As in [START_REF] Balthasar | A stochastic approach to calculate the particle size distribution function of soot particles in laminar premixed flames[END_REF][START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF], soot dynamics are simulated in a moving reference frame, following the fluid particles from fresh to burnt gases. For the sake of comparison, as in [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF], a shift by +0.25 cm of the computational distance above the burner is added, to account for the probe cooling effect.

To compare against experiments, soot number density is expressed as n Particle volume freely evolves according to Eq. ( 16) during growth and oxidation. According to the algorithm discussed above, for nucleation and agglomeration, intervals of sizes I vi must be defined. The nucleation mass m o = 2m py , corresponds to the mass of nascent soot particles from the collision of two pyrene molecules. The left boundary of

I vo is v o = m o × ρ soot .
This corresponds to a diameter of 0.88 nm. A geometric grid is used, following:

v inf i = v o F i s , (44)
Here the value of F s = 1.5 is used. The upper diameter of the grid is 51 nm. Another size range and other values of F s are also tested in section 6. The stochastic particles are initialised at the lowest section and N T (t = 0) = 1.0 cm -3 . This initial distribution is only necessary to start the calculation but it is largely negligible as compared to N T levels reached during the calculation (see Figure 2(a)). This approach has been validated elsewhere against analytical test cases featuring simultaneous nucleation and growth [START_REF] Bouaniche | A hybrid stochastic/fixed-sectional method for solving the population balance equation[END_REF].

Results

Figure 2(a) shows that the simulation approximates well N T , the total number density. The satisfactory prediction of pyrene mass fraction, aside from a small over-prediction for H > 1 cm (Fig. 1(b)), brings confidence on the calculation of the nucleation term. Therefore, observing good results on total number density suggests that agglomeration is well resolved. The over-prediction of N T against the reference calculation [58] postprocessing all diameters seen for higher values of H, the distance to the burner, may be explained by the stronger nucleation induced by the small over-prediction of pyrene. Soot volume fraction also agrees with the previously published simulations (Fig. 2(b)). This was expected as the same physical models for the nucleation, agglomeration and surface growth source terms were used. However, compared to measurements, soot volume fraction is significantly underestimated in both simulations above H = 0.75 cm, which may indicate an underestimation of surface growth. A possible cause of surface growth underestimation could be the soot formation and growth models. Another possible cause could be the models used for aggregate geometry through the calculated collision diameters and reactive surface. In the present paper, new source terms models are out of scope. Only the α parameter is varied and the influence of aggregate shape modelling on the Particle Size Distribution is looked at (section 5). Indeed, it is not obvious that soot particles can be modelled as perfect spheres in premixed ethylene flames as reported in [START_REF] Schenk | Morphology of nascent soot in ethylene flames[END_REF]. Figure 3 shows PSD and confirms the capability of the model to capture the soot size distribution. A unique nucleation mode is seen in both reference and present simulation for the lower streamwise position H. The transition to the bimodal distribution, with an additional accumulation mode, takes place approximately at the same height above burner. Finally, at H = 1.2 cm both PSD feature similar shapes. The present calculation shows higher number densities for small particles, which is in line with the previous observation that nucleation and total number density are slightly enhanced in the current simulation for higher values of H. In this first attempt, the numerical noise was minimized by performing the simulation with a set of N P = 10000 stochastic particles. The crosses in Fig. 3 show results obtained with N P = 1000, without much differences. As further assessed in section 6, the moderate CPU cost of the proposed method allows to envision the use of the hybrid strategy in LES within either Eulerian or Langrangian context for transporting the stochastic particles, defining N P , the total number of stochastic particles, from those present within a mesh cell [START_REF] Haworth | Progress in probability density function methods for turbulent reacting flows[END_REF].

Figure 4 presents PSD with number densities normalized by total number density. The computations (current and reference one) predict the right trend of the evolution, from a pure nucleation mode toward a bimodal distribution. The localization of the number density dip between both modes at about 4 nm compares well to experiment. Nevertheless, both computations fail to reproduce the number densities in the accumulation mode. Measured diameters are significantly larger than simulated ones, which goes with the underestimation of soot volume fraction in both simulations.

As discussed in [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF], this underestimation could result from the disturbance by the probe, responsible of flow stagnation in front of the sampling orifice, increasing particles residence time and therefore their sizes. A second candidate for this difference is the mobility diameters measured by SMPS, which might not be directly comparable to diameters simulated assuming perfect spheres and constant density of the solid. A point which is now further examined.

Parametric analysis of aggregate shape modeling

The influence of particle shape modeling on computed PSD and on integral values, as soot volume fraction, is now explored. The measured mobility diameters, which were directly used for comparison in Fig. 4, are converted into equivalent sphere diameters assuming specific fractal shapes, thus enabling more meaningful comparisons between experiment and simulations. Then, the agglomeration modeling is improved with a particle shape based on the fractal law, to better estimate, in additional simulations, the collision diameters in the calculation of the β vi,vj , the frequency of collision of the particles and the aggregate surfaces. Using the semi-empirical closure for the effective density measurements of the fractal aggregates proposed by Yon et al. [START_REF] Yon | A simple semi-empirical model for effective density measurements of fractal aggregates[END_REF], an equivalent sphere did ev can be inferred from d m , the measured mobility diameters. The calculation is made under the pressure and temperature conditions of the experiment, with a primary particle density ρ soot = 1800 kg.m -3 and a fractal dimension D f = 1.73, i.e. corresponding to the soot properties determined in the conversion model used in the present paper (Table 1 of [START_REF] Yon | A simple semi-empirical model for effective density measurements of fractal aggregates[END_REF] for 'CAST' aggregates). Consequently, the corresponding parameters have been determined in a propane diffusion flame. Nevertheless, the proposed range is in agreement with most of the values reported in the literature. Furthermore, by keeping in mind the large experimental uncertainties related to the determination of these parameters, there is, up to now, no evidences of a possible variability of these parameters to the fuel and the nature of the flame. For example, the authors of [START_REF] Schenk | Morphology of nascent soot in ethylene flames[END_REF] measured D f ranging from 1.2 to 1.8 in a series of premixed ethylene flames. Similar values have been used for simulating ethylene premixed flames, for example in [START_REF] Saggese | Kinetic modeling of particle size distribution of soot in a premixed burner-stabilized stagnation ethylene flame Combust[END_REF] the authors used a constant value of D f = 1.8. Concern-ing the density of the primary particles, ρ soot = 1800 kg.m -3 is the most encountered value related to black carbon and mature soot. The existence of a strong dependence of this parameter to the content of organic carbon has been recently shown [START_REF] Ouf | True density of combustion emitted particles: A comparison of results highlighting the influence of the organic contents[END_REF]. This reference reports a density up to 1834 kg.m -3 for low OC content (below 5 percent) and smaller values for organic particles. Nevertheless, up to now, it is difficult to assess the amount of organic carbon of soot in flames. In consequence, the hypothesis of mature and inorganic soot is done. This hypothesis is done by experimentalists, for example, the value of 1860 kg.m -3 was used for the particle size determination based on LII models [START_REF] Bladh | Probing the smallest soot particles in low-sooting premixed flames using laser-induced incandescence[END_REF]. It is also a common value for numerical studies on ethylene premixed flames in the literature. Similarly, the values of 1800 kg.m -3 and 1860 kg.m -3 were used in [START_REF] Eberle | Large Eddy Simulations of a Sooting Lifted Turbulent Jet-Flame[END_REF] and [START_REF] Aubagnac-Karkar | Sectional soot model coupled to tabulated chemistry for Diesel RANS simulations[END_REF] respectively. The number density distribution versus d ev , calculated from the measurements of d m , are shown in Fig. 5 for various d p , the diameter of the primary particles. For the higher values of d p , the particles are mostly spherical and d ev does not significantly differ from the measured mobility diameter. The departure from the spherical shape becomes significant as the characteristic size d p of the primary particles decreases, leading to aggregates composed of many primary particles.

Additional simulations are performed in which the mobility d m and collision d c diameters, which enter the calculation of the collision frequencies β vi,vj , are estimated from d g , the gyration diameter. The constant of proportionality is evaluated by considering the limit of a single sphere, as in [START_REF] Blanquart | A joint volume-surface-hydrogen multi-variate model for soot formation[END_REF][START_REF] Kazakov | Dynamic modeling of soot particle coagulation and aggregation: Implementation with the method of moments and application to high-pressure laminar premixed flames[END_REF][START_REF] Kruis | A simple model for the evolution of the characteristics of aggregate particles undergoing coagulation and sintering[END_REF]. d g is estimated according to the fractal law

n p = k f (d g /d p ) D f . ( 45
)
n p is the number of primary particles in the aggregates, D f is the fractal dimension, and k f is the fractal pre-factor. D f = 1.73 and k f = 1.94 are assumed constant for aggregates verifying n p > 10 [START_REF] Lattuada | A simple model for the structure of fractal aggregates[END_REF][START_REF] Yon | A simple semi-empirical model for effective density measurements of fractal aggregates[END_REF]. n p = v k /(πd 3 p /6) can be inferred from v k , the aggregate volume solved with the PDF of sizes. The only parameter left to determine d g from ( 45) is then d p . An attempt to estimate d p is done by performing a preliminary simulation accounting only for nucleation and surface growth, thus without agglomeration, leading to d p of the order of 5 nm at the streamwise position H = 1.2 cm.

To analyze the sensitivity of the modeling to these various particle parameters, simulations including the fractal law are performed first fixing d p = 5 nm and varying α HACA , calibrating in the modeling the surface reactivity of the particles (Fig. 6). In a second series of simulations (Fig. 7), α HACA = 1 and d p is varied between 5 nm up to an upper value for which the fractal law provides results similar to those of spherical particles (denoted 'spheres', dotted line in Fig. 7). It is important to note that we considered the particles to be spherical when the corresponding volume is lower or equal to the one of a primary spherule i.e d ev ≤ d p .

Experimental results are fairly well reproduced with α HACA = 1.0 and d p = 5 nm (Figs. 6 and7). For the lower values of d p , the reactive surface available for C 2 H 2 addition is higher, which explains higher values of the volume fraction as observed in Figure 8. Moreover, for small values of d p , aggregates are more open and d c , the collision diameter, increases, which in turn increases collision rates (in particular in the free molecular regime) and the distribution is shifted towards larger particles (Fig. 7). This exercice illustrates the strong impact of particle shape modeling when comparing against experiments. In particular, spherical particles might not always be the best choice. Also, bi-variate [START_REF] Patterson | Models for the aggregate structure of soot particles[END_REF] or even tri-variate modeling [START_REF] Blanquart | A joint volume-surface-hydrogen multi-variate model for soot formation[END_REF] could constitute an interesting alternative to calculate d p .

Comparison between hybrid and sectional methods

In this section, the present hybrid method is compared to a representative sectional method. The source terms for nucleation, growth and agglomeration are calculated in the same way for both methods. While growth is solved directly through stochastic particles in the hybrid method, a discretization scheme is used for the sectional method as in [START_REF] Park | A novel fixed-sectional model for the formation and growth of aerosol agglomerates[END_REF]:

∂N i (x, t) ∂t = G(v i-1 )N i-1 (x, t) v i -v i-1 - G(v i )N i (x, t) v i+1 -v i . (46)
The impact of mesh resolution is studied hereafter. Again, a geometric grid is used as defined by equation [START_REF] Patterson | Stochastic weighted particle methods for population balance equations[END_REF]. The centered section sizes range from 0.88 to 225 nm. This corresponds to a base case grid with F s = 2 and 25 sections which are typical values in the literature as for example in [START_REF] Eberle | A sectional PAH model with reversible PAH chemistry for CFD soot simulations[END_REF]. As the number of sections M is varied, the size range is kept constant and F s is varied accordingly. The results for height above the burner of 1.2 cm can be seen in Figures 9 and10 respectively in logarithmic and in linear scales. For sake of clarity, the figures do not include the results for all grids for which calculations were run.

Considering the 200 sections simulation as a reference, the hybrid method yields better results than the sectional one using the same number of sections. For more quantitative comparison, the error against the 200 sections simulation can be calculated for each mesh. The Earth Mover's Distance (EMD) [START_REF] Rubner | A metric for distributions with applications to image databases[END_REF], also known as Wasserstein-1 distance, is used as metric of error. It is calculated using the 'emdist' package of the CRAN project [START_REF] Urbanek | Package emdist[END_REF]. The EMD is then normalised by the standard deviation of the reference distribution. The evolution of this normalised EMD is given as a function of the number of sections in Fig- ure 11. The hybrid method yields lower error (with the same mesh refinement) for the whole range of tested grids. On the other hand, it presents a moderate CPU cost overhead as compared to the sectional method, as can be seen in Figure 12. The normalized CPU cost presented here is the CPU cost of the simulation divided by the CPU cost of the base case simulation (sectional, 25 sections). In order to assess the tradeoff between accuracy and cost the EMD error metric is plotted as a function of normalized CPU cost in Figure 13. For better readability, only the calculations up to 45 sections are shown on this last plot. Better accuracy is reached for lower CPU cost in the range of tested meshes. The hybrid method seems to offer an advantageous tradeoff in terms of CPU cost and accuracy, at least on this test case. Generalization of these conclusions will depend on the relative importance of agglomeration and growth source terms as the advantage of the present method over fixed sectional methods lays in accurate surface growth solving. A hybrid stochastic/fixed-sectional approach has been proposed to simulate the dynamics of soot particles in flames. The method relies on a fixed number of stochastic particles, whose evolution allows for solving the equation for the soot size probability density function, which is jointly considered with the equation for total number density. This hybrid approach is coupled to PAH chemistry and well-established modeling for nucleation, growth, agglomeration and soot oxidation and applied to the simulation of a premixed sooting flame of the literature [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF]. The comparison against previous numerical simulation and experiments, confirms the potential of the proposed strategy. In a second part, the simulations are used to investigate the influence of the modeling of the aggregate shape. It is in particular illustrated how this impacts on particle size distribution. Finally, the interest of the novel approach in term of tradeoff between cpu cost and accuracy has been quantified. Particle size distributions obtained from the hybrid method are comparable to those measured and calculated in [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF]. In the present work, the use of the hybrid method enabled a drastic reduction of the number of stochastic particles used (10 000 or 1000 depending on the runs as opposed to 200 000 in [START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF]). The extension to three dimensional flows of the soot population balance equation solution method would require the transport of stochastic particles in physical space. This can easily be done using well-established approaches, either in Lagrangian [START_REF] Haworth | Progress in probability density function methods for turbulent reacting flows[END_REF] or Eulerian context [START_REF] Pope | Monte Carlo method for the PDF equations of turbulent reacting flow[END_REF].

Fig. 1

 1 Fig.1Species mass fractions. Lines: present simulation. Symbols: reference[START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF].

Fig. 2

 2 Fig. 2 Soot number density and volume fraction vs height above burner. Symbols: experiments. (a) and (b) (diameters larger than 3 nm) Solid line: present simulation. Dashed line: reference simulation [58]. (a) (for all diameters) Dash-dot line: present simulation. Dotted line: reference [58].

Fig. 3

 3 Fig. 3 Particle Size Distribution for several heights above the burner. Solid line: present simulation with 10000 stochastic particles. Cross: 1000 stochastic particles. Dashed line: reference simulation [58].

Fig. 4

 4 Fig. 4 Normalized Particle Size Distribution for several heights above the burner. Solid line: present simulation. Dashed line: reference simulation [58]. Symbols: experiments.

Fig. 5

 5 Fig. 5 Normalized Number Density at H = 1.2 cm vs the sphere equivalent diameter Circles: dm measured mobility diameters. Crosses: dev with dp = 30 nm. Triangles: dev with dp = 15 nm. Plus symbols: dev with dp = 5 nm.

Fig. 6

 6 Fig. 6 Normalized PSD at H = 1.2 cm. dp = 5 nm. Solid line: α HACA = 1.0. Dashed line: α HACA = 0.9. Dash-dot line: α HACA = 0.8. Dotted line: α HACA = 0.7. Crosses: experiments dev with dp = 30 nm. Plus symbols: experiments dev with dp = 5 nm.

Fig. 7

 7 Fig. 7 Normalized PSD at H = 1.2 cm. α HACA = 1.0. Solid line: dp = 5 nm. Dashed line: dp = 7 nm. Dash-dot line: dp = 10 nm. Dotted line: sphere model. Crosses: experiments dev with dp = 30 nm. Plus symbols: experiments dev with dp = 5 nm.

Fig. 8

 8 Fig.8Soot volume fraction vs height above burner. Solid line: dp = 5 nm. α HACA = 1.0. Dash-dot line: spheres. α HACA = 1.0. Dotted line: spheres. α HACA from[START_REF] Appel | Kinetic modeling of soot formation with detailed chemistry and physics: laminar premixed flames of C2 hydrocarbons[END_REF] as in[START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF]. Dashed line: reference[START_REF] Zhao | Measurement and numerical simulation of soot particle size distribution functions in a laminar premixed ethylene-oxygen-argon flame[END_REF]. Symbols: experiment.

3 ]Fig. 9

 39 Fig. 9 PSD at H = 1.2 cm. Logarithmic scale. Dotted line: Sectional M = 25. Dashed line: Sectional M = 45. Solid line: Sectional M = 200. Line with empty circles: Hybrid M = 25. Line with full circles: Hybrid M = 45.

Fig. 10

 10 Fig. 10 PSD at H = 1.2 cm. Linear scale. Dotted line: Sectional M = 25. Dashed line: Sectional M = 45. Solid line: Sectional M = 200. Line with empty circles: Hybrid M = 25. Line with full circles: Hybrid M = 45.

Fig. 11

 11 Fig. 11 Normalised EMD as a function of the number of sections M (reference: 200 sections, solution at x = 1.2 cm). Dashed line: Sectional method. Solid line: Hybrid method.

Fig. 12

 12 Fig. 12 Normalised CPU cost as a function of the number of sections M . Dashed line: Sectional method. Solid line: Hybrid method.

Fig. 13

 13 Fig. 13 Normalised EMD as a function of normalised CPU cost. Dashed line: Sectional method. Solid line: Hybrid method. Labels: Number of sections

Acknowledgments

This work is supported by the European Union under the project SOPRANO, Horizon 2020 Grant Agreement No. 690724

Compliance with Ethical Standards

Conflict of Interest The authors declare that they have no conflict of interest.