Large Eddy Simulation of Turbulent Flame Synthesis of Silica Nanoparticles with an Extended Population Balance Model

In the present study, a recently proposed extended population balance equation (PBE) model for aggregation and sintering is incorporated into a large eddy simulation-probability density function (LES-PDF) modelling framework to investigate synthesis of silica nanoparticles in a turbulent diffusion flame. The stochastic field method is employed to solve the LES-PBE-PDF equations, characterising the influence of the unresolved sub-grid scale motions and accounting for the interactions between turbulence, chemistry and particle dynamics. The models for gas-phase chemistry and aerosol dynamics are the same as those recently used by the authors to simulate silica synthesis in a laminar flame (Tsagkaridis et al. in Aerosol Sci Technol 57(4):296–317, 2023). Thus, by retaining the same kinetics without any adjustments in parameters, we focus on the modelling issues arising in silica flame synthesis. The LES results are compared with experimental in-situ small-angle X-ray scattering (SAXS) data from the literature. Good agreement is found between numerical predictions and experimental data for temperature. However, the LES model underestimates the SAXS data for the primary particle diameter by a factor of two. Possible reasons for this discrepancy are discussed in view of the previous laminar flame simulations.


Introduction
Aerosol flame synthesis is one of the most widely used methods for the manufacturing of nanoparticles such as carbon blacks, fumed silica ( SiO 2 ), pigmentary titania ( TiO 2 ), zirconia ( ZnO 2 ) and alumina ( Al 2 O 3 ) (Rosner 2005;Buesser and Pratsinis 2012;Park and Park 2015).Such nanoparticles find a wide range of applications in gas sensors, energy storage materials, nanotoxicity studies, catalysts, biomaterials, and electroceramics (Meierhofer and Fritsching 2021), amongst others.Accurate modelling of all the processes involved in flame synthesis of nanoparticles is of paramount importance and can lead to better control over particle properties.This, in turn, allows synthesizing nanoparticles with specific characteristics tailored to particular applications.
Several processes take place during the formation of nanoparticles in flames, including fuel combustion, precursor decomposition and oxidation, nucleation, condensation, aggregation and sintering.In most industrial cases, aerosol flame synthesis occurs in high-Reynolds-number turbulent flows with typical Re numbers of the order of 10 6 (Buesser and Pratsinis 2012).Product characteristics are highly affected by flow configuration, turbulent mixing of reactants, flame temperature, precursor loading, etc (Raman and Fox 2016).Prediction of particle morphology in different flow regimes, such as laminar or turbulent, remains a challenge (Camenzind et al. 2008).
Flame-made aerosols are polydisperse, i.e. they are composed of particles of different sizes, and the particle size distribution (PSD) is an essential characteristic of nanoparticle products.The spatial and temporal evolution of the PSD, which is determined by aerosol and fluid dynamics, can be described by the population balance equation (PBE), also known as the general dynamic equation (GDE) for aerosols.Numerical methods for the solution of the PBE include moment, Monte Carlo and discretisation (also called sectional) methods.A review of these methods can be found in Rigopoulos (2010).
Simulation of flame synthesis of nanoparticles in turbulent flows should incorporate several models and solution methods.More specifically, the individual model components for the entire synthesis process include the description of turbulent flow, gas-phase chemical kinetics, precursor decomposition and oxidation kinetics, PSD evolution, description of physicochemical aerosol processes, and interactions between all these elements.The accuracy of the simulation is subject to the uncertainties and assumptions involved in all of these elements.The use of simplified models for the reaction kinetics, aerosol dynamics, and their interaction due to turbulent flow may lead to compensating errors and complicate the validation efforts.
In the following, we review existing strategies for CFD modelling of nanoparticle flame synthesis in turbulent flows, focusing on the modelling aspects adopted in each case.A more extensive discussion of CFD-based studies dealing with silica synthesis have been presented in Tsagkaridis et al. (2023) and is summarised here in Table 1.We also refer the reader to the review papers of Buesser and Pratsinis (2012), Buesser and Gröhn (2012), Meierhofer and Fritsching (2021) and Raman and Fox (2016).
Early CFD-based studies of flame synthesis of nanoparticles were based on the Reynolds-averaged Navier-Stokes (RANS) approach.Equations describing particle dynamics (Kruis et al. 1993) were integrated along characteristic trajectories through the flame reactors.Researchers employed this methodology to investigate the flame synthesis of palladium (Tsantilis et al. 1999), alumina (Johannessen et al. 2000) and titania (Johannessen et al. 2001) nanoparticles.Silica flame synthesis was investigated by Lee et al. (2001); in their model, a one-dimensional laminar flame was considered while particle

Table 1
Selected CFD-based studies of flame synthesis of silica nanoparticles and summary of the main features of the modelling approach adopted in each case The information presented in the table corresponds to: reactor type (DF: diffusion flame, PM: premixed flame, TML: temporal mixing layer, FSP: flame spray pyrolysis), the precursor, the CFD approach for the flow simulation (Lam: laminar flow), combustion model (FRC: finite-rate-chemistry, ED: eddy dissipation model, SE: species equilibrium approach), models for describing particle dynamics (whether turbulence-chemistry and turbulence-particle dynamics interaction were considered in the model or not, whether a comparison with experimental data was attempted or not), and a comment dynamics were described by a sectional two-dimensional PBE model due to Xiong and Pratsinis (1993).
A fully coupled fluid and particle dynamics study of flame synthesis via CFD was conducted by Schild et al. (1999).A simple monodisperse model due to Kruis et al. (1993) was coupled with an available commercial CFD solver to perform laminar-flow simulations of TiO 2 nanoparticle formation in a hot-wall reactor.Afterwards, Mühlenweg et al. (2002) extended this methodology, and CFD was coupled with a sectional PBE model showing the feasibility of this approach.
Synthesis of TiO 2 nanoparticles in a turbulent diffusion flame was simulated by Yu et al. (2008aYu et al. ( , 2008b) by employing a k-turbulence model together with the eddy dissipation concept for turbulent combustion modelling, and good agreement was found between numerical predictions for temperature and experimental data.Particle dynamics was described with an efficient quadrature method of moments (QMOM).Flame synthesis of SiO 2 nanoparticles in turbulent diffusion flames was simulated by Gröhn et al. (2011).
A realisable k-model together with a chemical equilibrium assumption were employed to account for flame dynamics, while a monodisperse model described particle dynamics.Later, Gröhn et al. (2012) extended the methodology to model flame spray pyrolysis (FSP) of ZnO 2 nanoparticles.
A methodology for describing flame synthesis of nanoparticles that also accounts for the turbulence-chemistry and turbulence-nanoparticles interaction (in the context of RANS) was developed by Akroyd et al. (2011).They employed a projected fields method to solve the joint composition probability density function (PDF) transport equation in conjunction with a method of moments with interpolative closure (MOMIC) for aerosol dynamics.
An approach based on large eddy simulation (LES) for TiO 2 nanoparticle synthesis in turbulent flames was proposed by Sung et al. (2011).Detailed chemical reaction mechanisms for both fuel combustion and precursor oxidation were employed, while the turbulence-chemistry interaction was modelled using on a flamelet approach.Aerosol dynamics were described by a univariate QMOM method, while sintering was not included in the model.Subsequently, Sung et al. (2014) improved on this approach by using a bivariate conditional quadrature method of moments (CQMOM) that can incorporate the sintering process into the model and thus describe aggregate morphology.The prediction of particle surface area was significantly improved compared to the results with the univariate model, while the primary particle size was under-predicted.The discrepancies were attributed to uncertainties involved in the turbulent inflow conditions.In both studies, the terms associated with the unresolved scales arising in the LES-filtering of the PBE equation were neglected.
An LES methodology for describing flame spray pyrolysis for synthesis of SiO 2 nano- particles was presented by Rittler et al. (2017).A model based on a premixed flamelet generated manifold approach (PFGM) along with an artificial flame thickening method was proposed for the system, and particle-dynamics evolution was modelled with a monodisperse model.Later, this methodology was extended by Wollny et al. (2020), who presented a multiscale model within the LES framework that included a sectional method for solving aerosol dynamics.However, the effect of the unresolved turbulent scales on aerosol dynamics was neglected.Vo et al. (2017) presented an LES methodology based on a sparse-Lagrangian multiple mapping conditioning (MMC) model for describing nanoparticle formation in turbulent flows.Later, Neuber et al. (2019) extended this methodology and coupled it with a PBE model solved via a discretisation method.Their PBE-MMC-LES approach, called accounts for the interaction between turbulence, chemistry and particle dynamics by solving the joint scalar PDF of gas species and discretised PBE.This methodology was employed to simulate flame synthesis of SiO 2 nanoparticles in a jet by oxidation of SiH 4 .However, finite-rate sintering was not included in the model and a constant primary particle diameter of 0.98 nm was considered.Results were compared with experimental data and the discrepancies were attributed to uncertainties involved in the precursor-decomposition reaction mechanism.
Direct Numerical Simulation (DNS) studies of flame synthesis of nanoparticles have also started to emerge in the literature.Abdelsamie et al. (2020) performed three-dimensional simulations to describe ( TiO 2 ) flame synthesis in conditions mimicking the Spray- Syn burner (Schneider et al. 2019), where the particle evolution was described with the monodisperse model of Kruis et al. (1993).Tsagkaridis et al. (2022) investigated turbulence-coagulation interaction via DNS and discretised PBE in the context of a flow laden with aerosol particles (not including particle formation processes).Flame synthesis of iron oxide nanoparticles was investigated via two-dimensional DNS and a sectional model by Cifuentes et al. (2020), where particular emphasis was given to flame-vortex interactions and the effect of high particle Schmidt numbers.
From the preceding review, it can be seen that, in most of the reported CFD studies on turbulent flame synthesis, the complex interaction between turbulence, chemistry and particle dynamics was not fully accounted for.Furthermore, simplified methods were often employed with regards to particle morphology.The main objective of the present study is to investigate numerically the synthesis of silica nanoparticles in a turbulent flame via a comprehensive LES methodology that also accounts for the turbulence-chemistry and turbulence-particle dynamics interaction.Particle morphology will be described by an extended population balance model recently proposed by Tsagkaridis et al. (2023).The extended PBE comprises a sectional method where the PBE is solved together with a transport equation for the number concentration of primary particles.In the present work, this model is introduced into the filtered probability density function/stochastic field modelling framework (denoted as LES-PBE-PDF) to predict the detailed experimental in-situ smallangle X-ray scattering (SAXS) measurements of Camenzind et al. (2008) for the "S-10"1 turbulent diffusion flame.The LES methodology utilises the same set of models for gasphase chemistry and aerosol dynamics that was recently employed by the authors (Tsagkaridis et al. 2023) to simulate the S-2 laminar diffusion flame of Camenzind et al. (2008).Results and possible sources of uncertainties in the model are discussed in light of the aforementioned laminar flame simulation results (Tsagkaridis et al. 2023).
The rest of the paper is structured as follows.The LES-PBE-PDF formulation is introduced in Sect. 2 together with a summary of the models employed for chemical kinetics and aerosol dynamics.This is followed by a description of the flame configuration and computational setup in Sect.3. Subsequently, the LES results are presented in Sect.4, where the discussion focusses on the discrepancies with the experimental SAXS data.An investigation of the effect of the particle sub-grid scale mixing is also shown, before the work concludes with a summary of the main findings in Sect. 5.

Governing Equations
The equations are written in a Cartesian frame of reference using tensor notation, where summation over repeated indices is implied.The conservation equations for mass and momentum for low-Mach-number flows are: where is the dynamic viscosity of the gas mixture.The isotropic part of the viscous stress has been adsorbed into the pressure.
The mass conservation equation for the chemical species is2 : where Y is the mass fraction of species while ω is the chemical reaction rate.Equal dif- fusivity, D, is assumed for all species and related to the dynamic viscosity of the mixture through the Schmidt number, Sc = ∕( D) , while the value Sc = 0.7 is employed for all gas species.The energy equation is expressed in terms of specific enthalpy that includes the thermal and chemical enthalpies: where ̇qrad and ̇qrad,p are the sink terms for enthalpy losses due to the thermal radiation from gas species and particles, respectively.A unit Lewis number has been assumed (Peters 2000) meaning Sc h = Sc k , while acoustic interactions and viscous heating have been neglected in the overall formulation.In the present work, multicomponent ideal gas flow was considered and the mixture density is a function of the local values of the reactive scalars and temperature.
The evolution of the particle size distribution is described by the PBE (Friedlander 2000).Anticipating the application of the PBE to variable density flows (Sewerin and Rigopoulos 2018), we consider the number density of particles per unit of mixture mass and the PBE assumes the following form: (1) t + ( u j ) x j = 0 (2) where v and w denote volume of aggregates, n = n(x, t, v) is the number density function per unit of mixture mass (denoting number of aggregates per unit of mixture mass and unit of particle volume).The explicit dependence of n on x and t has been omitted.B and C represent the nucleation and condensation rates, respectively, while (v − v nuc ) is the Dirac delta function and (v, w) is the aggregation kernel.The detailed expressions of the kinetic models employed will be discussed in Sect.2.2.The velocity, U j , is equal to the sum of the flow velocity u j and the thermophoretic velocity, and the latter is given by (Friedlander 2000): Equal diffusivity coefficients have been assumed for all particle sizes, while a constant value Sc p = 1000 is assigned for the number density to account for the low diffusivity of particles.
The evolution of particle-morphology can be described (Tsagkaridis et al. 2023) by solving the PBE together with a transport equation for the total number concentration of primary particles per unit of mixture mass, N p : where J is the nucleation rate, v p,av is the average primary particle volume, V the particle volume concentration per unit of mixture mass, M 23 = ∫ ∞ 0 v 2∕3 n(v) dv is the 2∕3 fractional moment of the PSD and s is a sintering characteristic time.Description of particle morphology can then be obtained from knowledge of n(v) and N p .In particular, v p,av is approxi- mated as: where a monodisperse size distribution has been assumed for the primary particles.The primary particle diameter d p,av is then obtained from v p,av by assuming spherical particles.Furthermore, fractal-like aggregates follow a power-law relationship (Friedlander 2000;Buesser and Pratsinis 2012).The aggregate collision diameter, d g , of aggregates of volume v is given by: where k f is a prefactor associated with the aggregate anisotropy (Heinson et al. 2010) and D f is the fractal dimension.The values k f = 1.4 and D f = 1.91 were assigned to the aggre- gate parameters corresponding to ballistic cluster-cluster aggregation in the free molecule (5) (Ball and Jullien 1984;Friedlander 2000).The average radius of gyration of aggregates is evaluated as: where dv is the number concentration of aggregates and v d is the volume of the smallest detectable particles (1 nm) by SAXS measurements (Camenzind et al. 2008).
As discussed later, a sectional model is employed for the numerical solution of the PBE and therefore, for the purposes of the mathematical formulation, a discretised form of the number density is introduced.The fundamental idea of sectional methods is to discretise the particle size dimension (in our case the particle volume) into discrete sections and evaluate the particle number densities at these sections.In this manner, the number density function n(x, t, v) can be represented by an array of m scalars n i (x, t), i = 1, m , each expressing the number density of particles within a section i with volume in the range v i−1 , v i .Based on the above, equations ( 3) to ( 7) can be expressed as generic reactive scalar transport equations in terms of a scalar , = 1, N s where N s is the number of scalars required to describe the system, i.e. the number of species plus enthalpy plus m + 1 scalars required after the discretisation in v of the PBE model: where J ,j is the scalar diffusion flux given by: where Sc is the Schmidt/Prandtl number, while ωa ( ) is the source term for Y , h, n i or N p .

Models for Chemical Kinetics and Aerosol Dynamics
Tsagkaridis et al. ( 2023) simulated synthesis of silica nanoparticles by oxidation of hexamethyldisiloxane (HMDSO) precursor in an oxygen/methane laminar diffusion flame.The objective of the present study is to simulate synthesis of silica nanoparticles in a turbulent flame and employ the same set of models for the description of chemistry and aerosol dynamics as in the laminar flame simulation (Tsagkaridis et al. 2023).For this reason, the details of the models employed are discussed only briefly here.
The Gas Research Institute (GRI) 1.2 reaction mechanism (Frenklach et al. 1995b, a) containing 175 steps and 31 species is used to describe methane combustion.The global two-step mechanism proposed by Feroughi et al. (2017) and shown in Table 2 is employed to describe HMDSO combustion.The two-step HMDSO mechanism was also employed in our laminar flame simulation (Tsagkaridis et al. 2023) and has been used previously in the literature to study the synthesis of silica nanoparticles via flame spray pyrolysis (Rittler et al. 2017).Thermal radiation from gas species is considered in the simulation by employing the optically thin radiation model of Barlow et al. (2001).
Following an approach similar to that of Shekar et al. (2012), the nucleation process involves the formation of dimers and the rate is given by: where N A is Avogadro's number, C SiO 2 is the gas-phase concentration of SiO 2 monomers, K fm is the free molecule coagulation kernel, E F = 1.88 is the Van der Waals enhancement factor, k b is the Boltzmann constant, T is the local gas-mixture temperature, is the local gas-mixture density, and m SiO 2 and d o are the mass and diameter of a SiO 2 (g) molecule, respectively.It is assumed that SiO 2 (g) molecules have diameter d o = 0.44 nm and that the density of silica particles is p = 2196 kg∕m 3 .After discretisation of the PBE in the parti- cle volume domain, the nucleation source term in Eq. ( 5) takes the form: where v nuc is the volume of dimers (volume of nuclei), v m,nuc is the representative volume (midpoint) in the section that covers the volume of nuclei and dv nuc is the volume range of that section.
Condensation is considered as a collision process in the free-molecular regime between silica monomers and particles.The focus here is on the physical description of gas-to-particle mass conversion due to condensation.We refer to the papers of Sun et al. (2021) and Tsagkaridis et al. (2023) for the complete details on the implementation of the condensation model into the PBE.The rate of change of the particle volume concentration per unit of mixture mass due to condensation can be evaluated as: Table 2 Global two-step reaction mechanism for HMDSO combustion (Feroughi et al. 2017) The reaction rate constant is given in the form: where T stands for temperature, R for the gas constant, and E a for the activation energy.The units of the parameters A are expressed in terms of cm, s, and mol while E a is given in cal/mol Step where cond (v 0 , v) is the collision frequency of a monomer of volume v 0 with a particle of volume v and d g is the collision diameter of particles of volume v evaluated according to Eq. ( 9).The amplification factor of 1.28 is due to Van der Waals interactions.The expression for the sintering characteristic time of silica nanoparticles given by Tsantilis et al. (2001) is used in the present study: where d p is the primary particle diameter in meters, T the gas temperature in Kelvin and d p,min = 1 ⋅ 10 −9 m is the cut-off size below which particles are assumed to coalesce instantaneously.
Aggregation in the free molecule, transition and continuum regime is considered and appropriate expressions for the aggregation kernel, (v, w) , are employed depending on the values of the Knudsen number, Kn = 2 f ∕d g , where f is the mean free path.The expres- sions for the free molecule particle size regime (Kn < 0.1) and the continuum regime (Kn > 10) are (Friedlander 2000): where C c (v) = 1 + 1.257Kn is the Cunningham slip correction factor and has been included in the expressions since the number density function has been defined per unit of mixture mass.A harmonic mean of the two expressions (Pratsinis 1988) is employed for the transition regime (0.1 < Kn < 10).
Thermal losses due to the radiation from nanoparticles are included in the model and evaluated as:

Large Eddy Simulation
In LES, the larger energy-containing turbulent motions are resolved directly, whereas the effect of the unresolved sub-grid scale (SGS) motions needs to be modelled.The separation of the scales is achieved through a spatial low-pass filtering operation, which for a function Q = Q(x, t) is defined as its convolution with a filter function G according to: where integration is defined over the entire flow domain Ω .In the present study, a box or "top hat" filter is employed and Δ is a characteristic filter width which may vary with position.
In LES of combusting flows, the use of density-weighted (Favre) filtering, defined by Q = Q∕ , is usually employed to account for density variations in the unresolved scales.
Applying the density-weighted filtering operation to the governing transport equations results in the following filtered equations (Jones and Navarro-Martinez 2007;Jones and Prasad 2010): is the resolved rate of strain tensor and r ij is the anisotropic part of the sub-grid scale stress tensor sgs ij = (ũ i u j − ũi ũj ) .The isotropic part of the sub-grid stress has been absorbed into the pressure.Following standard methods, r ij is determined via a dynamic eddy-viscosity model: where and ‖ Sij ‖ = � 2 Sij Sij is the Frobenius norm of the resolved rate of strain tensor.The filter width Δ is defined as the cubic root of the local grid cell volume and the parameter C S is determined through the dynamic procedure of Piomelli and Liu (1995).In Eq. ( 25), the terms J sgs ,j = ũj − ũj ̃ and  ω ( ) (filtered net scalar formation/consumption rate) remain unclosed, as for instance  ω ( ) ≠  ω ( ) , and require modelling.The former term represents transport due to sub-grid fluctuations, while the latter is related to turbulence-chemistry and turbulence-particle dynamics interaction.

Joint Sub-grid Probability Density Function
An approach based on the joint sub-grid probability density function (PDF) method is employed in the present study to take into account the effect of the unresolved scales on the evolution of the scalar quantities.The interaction of turbulence and particle formation processes was analysed by Rigopoulos (2007), who developed a PBE-PDF equation for taking into account the unclosed correlations that result from the non-linear terms in the PBE, initially in a RANS context.The LES-PDF method for reacting flows was pioneered by Gao and O'Brien (1992) and Gao and O'Brien (1994), while the PBE was incorporated into the LES-PDF framework by Sewerin and Rigopoulos (2017) for constant density flows and by Sewerin and Rigopoulos (2018) for variable density flows.The formulation adopted in the present study is based on that given by Jones and Navarro-Martinez (2007) and Jones and Prasad (2010) and incorporates elements given by Sewerin and Rigopoulos (2018).
Following Gao and O'Brien (1992) and using the filtering operation defined in Eq. ( 22), a joint sub-grid (or filtered) PDF for the N s scalar quantities is defined as: where g( ; 'Brien 1980; Pope  1985), is the Dirac delta function, and denotes the phase or sample space of the scalar quantities .Based on the definition of f ( ; x, t) , the probability of observing values of the scalar variable in the range   <   <   + d  within the filter volume is given by f ( ; x, t) d .For variable density flows, a density-weighted PDF f ( ; x, t) is usually introduced, and a transport equation for this quantity can be derived (Gao and O'Brien 1992;Colucci et al. 1998;Jaberi et al. 1999;Sewerin and Rigopoulos 2018) from the appropriate conservation equations using standard methods: where the arguments of f ( ; x, t) have been omitted for compactness.The turbulent transport term has been approximated by a simple gradient closure directly analogous to the Smagorinsky model, and a constant turbulent Schmidt/Prandtl number Sc sgs = 0.7 is adopted in the present study (Jones and Prasad 2010).Note that in Eq. ( 29) the scalar formation/consumption term appears in closed form while the last term on the right-hand side in Eq. ( 29) represents sub-grid scale mixing and requires modelling.There are several modelling proposals for closing the micromixing term (Fox 2003).Following Jones and Prasad (2010), the linear mean square estimation (LMSE) micromixing model (Dopazo and O'Brien 1974;Dopazo 1975Dopazo , 1979) is adopted in the present study: where the sub-grid mixing (or micromixing) times scale is: The value C d = 2 is usually assigned to the micromixing constant for the chemical species and enthalpy (Jones and Prasad 2010).However, the diffusivity of nanoparticles is much lower than that of gas species; the diffusion coefficient of particles can be estimated by the Stokes-Einstein equation (Friedlander 2000).Therefore to account for this effect, a simple approach is adopted here and a different sub-grid mixing time scale is applied for the discretised number densities, n i (x, t) , which is assumed to scale with the particle Schmidt number Sc p = 1000 .A similar approach was employed by Vo et al. (2017) for the definition of a Lagrangian mixing time for the particles.In the present model, the different micromixing time of the particles can be expressed via the micromixing constant C d .The value C d,p = 0.002 is employed for the discretised number densities, while C d,s = 2 for the rest of the scalars.The possibility of using the same value of C d for all scalars is also explored in order to investigate the effect of the micromixing constant in Sect.4.3.

Stochastic Field Method
The stochastic field method, proposed by Valiño (1998) and Sabel'nikov and Soulard (2005), is employed to solve Eq. ( 29).The LES formulation of this method used here follows Jones and Navarro-Martinez (2007).The underlying idea of this method is the derivation of a system of stochastic differential equations (SDE) equivalent to the closed form of the filtered PDF transport equation.f ( ; x, t) is then represented by an ensemble of N f stochastic fields (x, t), = 1, N f for each of the reactive scalar = 1, N s .The evolution of the stochastic fields is described by the following equation: where dW j (t) is the increment of a three-dimensional Wiener process.dW j (t) is different for each field but not a function of the spatial location x and can be approximated by time-step increments in the j th direction as j √ dt where j is a {−1, 1} dichotomic random vector.The Favre filtered value of each reactive scalar , represented by the first moment of the scalar fields, can then be computed by averaging over the stochastic field values:

Solution of the PBE
A conservative finite-volume sectional method developed by Liu and Rigopoulos (2019); O'Sullivan and Rigopoulos (2022) was employed for the numerical solution of the PBE.The advantage of the method is that it provides an accurate prediction of the particle size distribution and conserves the first moment (with respect to volume) during the aggregation process, an important feature for mass balance.The method has been shown to be robust and efficient when coupled with CFD (Sun et al. 2021;Sun and Rigopoulos 2022;Tsagkaridis et al. 2023).

Coupling of PBE with CFD
The PBE model is incorporated within an in-house block-structured, boundary conforming coordinate LES code, BOFFIN-LES (Jones et al. 2002).The LES code has been used to simulate a wide range of combustion flows (Jones et al. 2015) second-order-accurate finite volume method based on a low-Mach number pressure-based formulation using a collocated variable arrangement.A second-order central discretisation scheme is employed for the spatial derivatives, while a second-order accurate Crank-Nicolson method is employed for time discretisation.A total variation diminishing (TVD) scheme based on van Leer's limiter is used for the scalar equation convection terms to avoid unphysical values.The convection and viscous terms are treated implicitly.A twostep approximate factorization (SIMPLE type) method for the pressure-velocity coupling is employed to ensure mass conservation.The stochastic field equations are solved using a weak first-order temporal Euler-Maruyama scheme (Kloeden and Platen 1992).A stochastic noise reduction procedure given by Prasad (2011) was employed in the present study.An approximate fractional step method is applied for the solution of the transport equations of the reactive scalars and discretised number densities.More specifically, separate fractional steps are employed for convection-diffusion, chemical reaction and PBE integration.The chemical reaction rate equations are integrated using an implicit backward Euler scheme, while the PBE is integrated with an explicit Euler scheme since the discretised PBE is not stiff.The reaction mechanism is hard-coded in order to accelerate the computations.The code is fully parallelized, utilising domain splitting and Message Passing Interface (MPI) routines.The chemical reaction fractional step is additionally parallelized to distribute reaction source computations evenly across all processors.

Flame Configuration
LES calculations were performed for the S-10 turbulent flame of Camenzind et al. (2008).To the best of the authors' knowledge, the present work is the first simulation of this experiment.Silica nanoparticles were synthesised by oxidation of HMDSO in oxygen/methane diffusion flames.The same co-flow diffusion flame reactor was used for both the laminar (S-2) and turbulent (S-10) cases.The reactor comprised three concentric stainless-steel tubes.The inner tube diameters were D = 1.8 mm, D 1 = 3.5 mm and D 2 = 4.8 mm, while the tube wall thickness was 0.3 mm (Wegner and Pratsinis 2003).To bring HMDSO to vapour phase, 0.3 l/min of N 2 passed through a reservoir filled with liquid HMDSO.Sub- sequently, the HMDSO-laden N 2 gas flow was mixed with 0.5 l/min of methane CH 4 and introduced into the reactor through the centre tube.The consumption rate of HMDSO was 6.5 g/h.The 1 st annulus was fed with 0.5 l/min of N 2 to slightly lift the flame and prevent particle deposition on the nozzle tips.The 2 nd annulus was fed with 10 l/min of O 2 .The temperature of the tubes was maintained at 75 • C to prevent precursor condensation on the tube walls.The Reynolds number, based on the bulk velocity and the outer diameter of the burner, was Re ≈ 4280 and the experiment was conducted under atmospheric pressure.It should be stressed that the only difference between the S-2 and S-10 flames is the O 2 flow rate through the 2 nd annulus.In the S-2 case, the O 2 flow rate was 2 l/min, while in the S-10 case, the O 2 flow rate was five times larger, resulting in a turbulent flame.

Computational Setup
The solution domain is cylindrical with a diameter 40D, where D = 1.8 mm is the inner tube diameter, and extends 60D in the streamwise direction.The simulation was performed on a non-uniform computational grid with approximately 4.65 million cells.The domain was split into 528 blocks (sub-domains) with 8800 cells per block, and a central processing unit (CPU) was assigned to each block.In total, 264 cells were employed for the discretisation of the domain in the streamwise direction and were distributed based on a geometric expansion formula, resulting in a minimum cell width of 8 × 10 −5 m in the vicinity of the nozzle.The cell width was 7.25 × 10 −5 m in the other two directions in the central region of the jet (r < 1.33D) .Grid refinement was employed at the height of the second annulus in the radial direction to ensure that at least eight cells were used to resolve the strong gradients that arise near the inlet.The cell width expanded smoothly in the radial direction for r > 1.5D.
A free slip boundary condition was used for the lateral surface of the cylindrical domain, and a convective outflow boundary condition was applied at the exit plane.In the simulation, the burner tubes were excluded from the solution domain, while the tube thickness (0.3 mm) was considered by applying a no-slip boundary condition at these radii at the inlet plane.At the inflow, laminar flow velocity profiles were employed for all the jets issuing into the reactor.A laminar flow profile was also employed for the jet coming from the second annulus; this choice was made in the absence of further information and on the basis that the estimated Reynolds number, Re, of the flow within the annulus was smaller than the critical Re c number for transition to turbulence in annular flow (Dou et al. 2010).Therefore, the turbulent flow downstream was generated due to shear (velocity gradients) and not by applying artificial turbulence at the inlet.In the absence of any information for the entrainment velocity, a uniform velocity U air = 0.8 m/s was used for the air co-flow stream at the inlet.Preliminary test cases showed that U air has a minor effect on the results.The inlet boundary conditions for the bulk streamwise velocity, temperature and mole fractions for all jet streams are shown in Table 3.
Following Mustata et al. (2006) and Jones and Navarro-Martinez (2007), eight stochastic fields ( N f = 8 ) were employed to approximate the sub-grid PDF introduced in Sect.2.4.A comparison of results obtained with eight and sixteen stochastic fields will be shown in Sect.4.4.The time step was fixed at Δt = 5 × 10 −7 s, resulting in a maximum Courant- Friedrichs-Levy (CFL) number of 0.4.The initial flow was allowed to reach a statistically steady state.Following this, the simulation was run for a period of at least three flowthrough times, based on the bulk velocity, for the collection of the statistics.Azimuthal averaging of statistical quantities was also performed.Each simulation required approximately 200,000 CPU hours (AMD EPYC 7742, 2.25 GHz, 64-core processor).
Regarding the PBE numerical parameters, an exponential grid with 60 intervals was employed for the discretisation of the particle volume domain and covered a particle size range from 0.554 nm to 1.4 m in diameter.Convergence studies with a perfectly stirred reactor showed that this grid was sufficient for an accurate description of the PSD evolution.The discretisation of the particle volume domain allows for the tabulation of the aggregation kernel.Specifically, the aggregation kernel for several particle-volume pairs and primary particle diameters was pre-computed and stored, leading to a substantial reduction in computational cost (Tsagkaridis et al. 2023).

Results and Discussion
In the present section, the LES results are presented and discussed.A mass-based definition of the number density function was presented in Sect. 2. However, particle concentrations are expressed per unit of mixture volume to make the results consistent with the experimental data and also allow for a direct comparison with the simulation data presented in our laminar flame simulation (Tsagkaridis et al. 2023).Comparisons of numerical predictions with the experimental SAXS data of Camenzind et al. (2008) are also presented.Time-averaged quantities are presented in all figures unless specified otherwise.

Flow Field and Temperature
The instantaneous and time-averaged fields for streamwise velocity and temperature are presented in Fig. 1a and b, respectively.It is observed that the flame is slightly lifted and not attached to the burner nozzle.In the experiment, N 2 was delivered into the reac- tor through the 1 st annulus in order to slightly lift the flame and prevent the deposition of particles on the nozzle tips, and this is captured well by the simulation.The flame is stable at low heights above burner (HAB), i.e. less than 0.4 cm; however, instabilities develop early on, leading to vortical structures that stretch the flame front and eventually break the jet.Intense turbulent mixing and a wide range of turbulent length scales are evident at downstream locations.Furthermore, the velocity field exhibits a recirculation The instantaneous and time-averaged fields for particle volume concentration (or silica volume fraction), V, are presented in Fig. 1c.Particles are formed by nucleation quite early in the flame.As will be shown later, most of the gas-to-particle conversion occurs in the flame reaction zone at low HABs (HAB < 2 cm).The maximum time-aver- aged silica volume fraction is found in the flame reaction zone at HAB < 0.4 cm, where the flame is relatively stable.At downstream locations, the jet stream becomes unstable and breaks into smaller vortices.This causes the appearance of fluid patches with large values of V, larger than the time-averaged value, at downstream locations.
A comparison of numerical predictions for temperature with the experimental Fourier transform infrared (FTIR) spectroscopy data of Camenzind et al. (2008) is shown in Fig. 2. As discussed in our earlier work (Tsagkaridis et al. 2023), the radial maximum temperatures are more representative of the experimental FTIR data; thus, numerical results for the radial maximum value of time-averaged temperatures are plotted in Fig. 2 against HAB.Error bars of ±60 K (Kammler et al. 2002) have been added in the figure for the experimental FTIR data.Overall, the numerical results for flame temperature are in good agreement with the experimental data.Focussing on the SAXS data at low HABs, it is observed that the radial maximum temperature increases very rapidly at HAB = 0.8 cm and decreases afterwards; this behaviour is reproduced well by the simulation.The simulation slightly underpredicts the experimental data at downstream locations (3 < HAB < 8 cm).This is likely due to uncertainties in the inflow bound- ary conditions and the micromixing model (Eq.( 30)).As mentioned earlier, no turbulent fluctuations were superimposed at the inflow velocity profiles.Nevertheless, the flame temperature is predicted accurately at low HABs in the region of interest where  (2008) for flame temperature along the centreline.An average experimental error of ±60 K given by Kammler et al. (2002) was used in the figure for the FTIR data nucleation, condensation, sintering and aggregation occur simultaneously.At downstream locations, the effect of all aerosol processes except aggregation is diminished, and furthermore SAXS data are not available for these locations.Based on these considerations, it is expected that the underestimation of temperature at downstream locations will not affect the further discussion of results in the present study.

Comparison with Experimental SAXS Data
In the present section, numerical predictions for particle quantities are compared with the experimental SAXS data of Camenzind et al. (2008) and the nature of the discrepancies is discussed.Following Tsagkaridis et al. (2023), we present simulation results along the pathline of silica volume fraction, where intense particle growth and high temperatures are found.It should be mentioned that the maximum silica volume fractions are found on the centreline at HAB > 1 cm due to the effect of turbulent mixing; therefore, the simulation results presented for HAB > 1 cm correspond to centreline results.
Numerical predictions for the primary particle diameters, d p,av , along the pathline of maximum silica volume fraction are shown in Fig. 3a, while numerical predictions for the average radius of gyration of aggregates, R g,av , are shown in Fig. 3b.R g,av is estimated by Eq. ( 10).The simulation predicts fully fused particles (spherical structures) at low HABs, while aggregate structures are predicted at downstream locations (HAB > 1.4 cm).This can also be seen from the small values of R g,av at low HABs and it is the aftermath of the high sintering rates predicted by the simulation at low HABs (see also Fig. 7 to be discussed later).On the contrary, the SAXS data suggest that fractal aggregates consisting of a large number of primary particles O(100) are formed very early in the flame (Camenzind et al. 2008).Aside from the experimental uncertainty of the SAXS data for d p,av , it can be seen from Fig. 3a that the LES simulation underpredicted particle sizes by a factor of 2. It is worth mentioning that the relative difference of d p,av with the SAXS data in our laminar flame simulation with the two-PBE model was about 15 % (see Fig. 6a of Tsagkaridis et al.  (2023)).The underprediction of R g,av observed in Fig. 3b is a consequence of the under- prediction of d p,av ; larger d p,av would result in smaller sintering rates, leading to aggregate structures.At present, it is not possible to determine the reasons for the discrepancy observed in Fig. 3a in a conclusive manner.However, some light can be shed by considering the extent to which the elements of the methodology have been validated in previous works.Primary particles can grow in size by sintering, the combined effect of aggregation and sintering (coagulation), and condensation.For clarification, coagulation is defined here as the particle collision process that results in fully fused particles, while aggregation is defined as the collision process that results in aggregated particles.The uncertainty involved in the models employed in the present study (presented in Sect.2.2) for the description of these aerosol processes will be discussed below.
Focussing firstly on the sintering process, it is anticipated that the selection of a different sintering model, which would give larger sintering rates, would not result in higher estimates for d p,av , since the simulation already predicts fully fused particles at HAB < 1.4 cm.In our laminar flame simulation, results for primary particles were in good agreement with the SAXS data (see Fig. 6a of Tsagkaridis et al. (2023)), giving evidence that the sintering process is represented well by the model.Good agreement between simulation and experimental data was also found for the geometric standard deviation (see Fig. 7 of Tsagkaridis et al. (2023)), suggesting that aggregation dynamics is described well by the extended one-PBE model.The condensation model described in Sect.2.2 is based on the assumption that collisions between particles and product molecules occur in the free molecule regime, an assumption readily satisfied at these elevated flame temperatures found at low HABs.From this discussion, it can be inferred that the models for sintering, aggregation, and condensation cannot explain the discrepancies observed in Fig. 3.
The evolution of primary particle sizes is also affected by the nucleation model; particles grow due to coagulation and coagulation/aggregation rates qualitatively follow nucleation rates at the early stages of nanoparticle synthesis (see, for example, Fig. 7b to be discussed later and Fig. 8a of Tsagkaridis et al. (2023)).In the present study, nucleation is treated as a collision process of gas molecules (monomers), resulting in the formation of dimers that serve as the first particles.Similar approaches are found in the literature (Shekar et al. 2012).The dimers-based model described adequately particle-morphology evolution in our laminar flame simulation and was in better agreement with the experimental data compared to the instantaneous-nucleation-assumption and classical-nucleationtheory models (Tsagkaridis et al. 2023).However, it should be acknowledged that accurate prediction of particle nucleation rates remains an unresolved problem in the field of nanoparticle flame synthesis (Rosner 2005).
Another source of uncertainty is associated with the precursor decomposition/oxidation reaction mechanism.The two-step chemical reaction mechanism of Feroughi et al. (2017) for HMDSO combustion (presented here in Table 2) was developed by optimising model parameters to fit experimental data obtained with low-pressure flames and low HMDSO concentrations.The experiment under investigation was conducted at atmospheric pressure and with relatively higher precursor concentrations.
There are also uncertainties in the closure models for describing turbulence-chemistry and turbulence-particle dynamics interaction.In the present analysis, the PDF-stochastic field method was employed to describe the influence of the unresolved scales and modelling was required for the micromixing term (Eq.( 30)).While uncertainty regarding the micromixing model for nanoparticles remains, it is likely that its effect is less important than the aspects of kinetics discussed above.This is consistent with the fact that the primary particles are underpredicted even at HAB <0.6 cm where the flame is still relatively stable and turbulence is not yet developed (see Fig. 1a and b).Furthermore, two different values of the micromixing constant for number densities were tested and results showed that d p,av is relatively insensitive to the micromixing constant C d,p , at least for the present flame configuration (see Fig. 5b to be discussed in Sect.4.3).Another parameter in the PDF-stochastic field method is the number of stochastic fields.Results obtained with eight and sixteen stochastic fields (to be discussed in Sect.4.4) were almost identical, indicating that this model parameter is not associated with the discrepancies observed in Fig. 3.
In light of the aforementioned arguments, the discrepancy between the simulation data for d p,av and the experimental SAXS data is attributed mainly to uncertainties involved in the nucleation model and the precursor decomposition kinetics.
Numerical predictions for the number concentration of primary particles, N p , along the pathline of maximum silica volume fraction are compared with the experimental data of Camenzind et al. (2008) in Fig. 4a.The model overpredicts the experimental data, which is consistent with our laminar flame study (see Fig. 6b of Tsagkaridis et al. (2023)).This discrepancy is most likely related to uncertainties involved in the nucleation model and the SAXS measurements.
Numerical predictions for the silica volume fraction, V, are displayed in Fig. 4b together with experimental data given by Camenzind et al. (2008).It should be emphasised that SAXS data analysis gives only an indication of the total volume fraction of silica species in the aerosol (Kammler et al. 2005); therefore, quantitative comparisons are not expected at this point.It is, however, interesting to examine the data qualitatively.In a series of preliminary simulations for both turbulent and laminar cases (Tsagkaridis et al. 2023), with different nucleation models, it was observed that the total gas-to-particle mass conversion predicted by the model was relatively insensitive to the nucleation model.It was the contribution of condensation versus nucleation that was altered but not the total gas-to-particle conversion.This implies that the qualitative discrepancy observed in Fig. 3b is most likely related to the precursor decomposition/oxidation reaction mechanism (presented in Table 2).Specifically, the experimental data suggest that V increases gradually, while the LES model predicts that V reaches a maximum value at HAB = 0.3 cm and decreases afterwards.This suggests that the two-step HMDSO reaction mechanism of Feroughi et al. (2017) is rather fast, at least for the flame conditions examined here.

Effect of Particle Sub-grid Mixing
It is of interest at this point to investigate the influence of the sub-grid contribution and examine the effect of the associated model parameters.The model parameter that directly affects the sub-grid contribution is the micromixing constant C d (Eq. ( 31)).Jones and  Prasad (2010) employed the PDF-stochastic field method to simulate a series of particlefree turbulent diffusion flames and examined the effect of the micromixing constant C d on the LES results.Specifically, Jones and Prasad (2010) demonstrated that simulations with different values for C d , i.e two and three, resulted in negligible differences in the pre- dicted flame temperatures.This suggests that the significance of the micromixing model is severely diminished in LES as compared to RANS (Jones and Prasad 2010).However, the impact of the micromixing model on particulate flames is more questionable due to the very low diffusivity of nanoparticles.
The diffusivity of nanoparticles is much lower than that of gas-phase species.To account for this effect in the present study, the value C d,p = 0.002 was assigned to the dis- cretised number densities, while the value C d,s = 2 was used for the chemical gas-species and enthalpy (see Sect. 2.4).To further investigate the effect of the micromixing constant, a simulation with C d,p = 2 (and C d,s = 2 ) was performed and comparisons of the two cases are presented in Fig. 5.
Numerical predictions for the flame temperature obtained with the two different values of the micromixing constant for number densities are presented in Fig. 5a.The flame temperature is relatively insensitive to the value of C d,p .This can be partly explained from that the fuel-combustion reactions are only weakly coupled to the precursor decomposition chemical reactions (see Table 2).Furthermore, the effect of thermal losses due to the radiation from nanoparticles, which could affect the flame temperature profile, seems to be less significant at downstream locations where lower temperatures are encountered; q rad,p scales with the fifth power of temperature (see Eq. ( 21)).Results for the primary particle diameter, d p,av , for the two cases, are shown in Fig. 5b.The two cases with the different values of C d,p resulted in identical results for d p,av .This is to be expected since sintering becomes negligible at roughly HAB < 2.1 cm as will be discussed in Sect. 4.5 and shown in Fig. 7b.
Numerical predictions for N p and V are presented in Fig. 5c and d, respectively.Small differences are found between the two cases at low HABs in the core region of the flame ( 0.8 <HAB < 2 cm) where nucleation, condensation, sintering, and aggregation occur simultaneously.Deviations between the two cases emerge at HAB > 2 cm, where intense turbulent mixing occurs.The case with the lower value of the micromixing constant ( C d,p = 0.002 ) resulted in larger values for N p and V at downstream locations.However, experimental data for particle properties are not available for most of this region.

Effect of the Number of Stochastic Fields
The results presented so far were obtained using N f = 8 stochastic fields for the solution of Eq. ( 29) through Eq. ( 32).The choice of the number of stochastic fields, N f , has been investigated by Mustata et al. (2006).Specifically, Mustata et al. (2006) used the PDFstochastic field method to describe a turbulent diffusion flame.By performing simulations with eight and sixteen stochastic fields, they showed that the two cases produced results with negligible differences.Similar findings were reported by Jones and Navarro-Martinez (2007).However, these studies dealt with particle-free flames.A similar analysis is performed in the present study but in the context of the LES-PBE-PDF method.Comparisons of LES results with eight and sixteen stochastic fields (for some representative variables) are presented in Fig. 6.Both simulation cases produced almost identical results, suggesting that eight stochastic fields are sufficient to obtain convergence and represent the sub-grid scale effects for both species and number densities.

Investigation of Gas-to-Particle Conversion
In the present section, the silica formation processes in the turbulent flame, as predicted by the LES method, are further analysed to provide further insight into the results demonstrated so far.We start with describing the SiO 2 (g) formation rate as it is driving many of the subsequent processes.The formation rate of SiO 2 (g) species, R G2 , (see G2 reaction step in Table 2) along the pathline of maximum silica volume fraction is shown in Fig. 7a.It is observed that the peak values are found at low HABs (HAB < 0.6 cm) where the flame is still relatively stable, while most of the precursor is converted to SiO 2 (g) at HAB < 1.65 cm.
The rate of change of particle number concentration due to aerosol processes along the pathline of maximum silica volume fraction is given in Fig. 7b.At HAB < 1.4 cm, the sintering rate is many orders of magnitude higher than the rates of the other processes.At these HABs, the simulation predicts fully fused particles, as already discussed in Sect.4.2a.Subsequently, the sintering rate decreases rapidly by several orders of magnitude due to the drop in flame temperature, and aggregated particles start to form at HAB > 1.4 cm.The effect of sintering becomes negligible at downstream locations (see intersection of curves at HAB = 2.1 cm), which explains why the primary particle sizes remain constant at these locations (see Fig. 3a).The nucleation rate reaches maximum values at HAB < 1.2 cm in the flame reaction zone where new SiO 2 (g) species are constantly formed, while it decreases by several orders of magnitude downstream where precursor conversion is almost complete.The predicted d p,av slightly decreases at 1.8 <HAB < 3 cm partly because of the effect of nucleation.Interestingly, nucleation of particles continues to take place until HAB ≈ 7 cm; however, its effect on d p,av and other particle quantities at 3 <HAB < 7 cm is less profound.The aggregation rate follows the nucleation rate qualitatively, while aggregation becomes the dominant mechanism at HAB > 7 cm.The PSD reaches a self-preserving state (Vemury  and Pratsinis 1995) at HAB > 7.3 cm (not shown here).
The rate of change of the silica volume fraction due to aerosol processes is portrayed in Fig. 8a, while results along the pathline of maximum silica volume fraction are given in Fig. 8b.There are two features worth noting.The first one is that most of the gas-to-particle conversion occurs at low HABs (HAB < 2 cm).Intense particle formation and growth take place in the flame reaction zone at HAB < 0.5 cm where the flame is relatively stable; the maximum silica volume fraction was found at HAB = 0.22 cm (see Fig. 1c and b).The sec- ond one is that the condensation rate is almost five times larger than the nucleation rate.Thus, condensation is responsible for the largest proportion of the mass addition into the particle phase.The same was also observed in our laminar flame simulation.Finally, a spatial correlation between nucleation and condensation is also observed in Fig. 8a, as expected.

Conclusions
Synthesis of silica nanoparticles in a turbulent diffusion flame has been investigated numerically.The main objectives of the present study was to apply a comprehensive LES methodology using the same set of models and parameters for gas-phase chemistry and aerosol dynamics as those employed recently by the authors (Tsagkaridis et al. 2023) in a laminar flame and finally compare the LES results with experimental data found in the literature.Particular emphasis was given to the nature of the discrepancies between predictions and experimental data, while possible sources of uncertainties in the model were discussed in light of our aforementioned laminar flame simulation results.
The LES-PBE-PDF approach was employed for the simulations.The model components encompass a comprehensive description of gas-phase chemistry and aerosol dynamics.The evolution of particle morphology was described with a recently proposed extended population balance model, whose solution is accomplished via the conservative finite-volume sectional method.Furthermore, the model predicts the evolution of primary particle sizes by solving the PBE together with a transport equation for the number concentration of primary particles.The extended one-PBE model was incorporated into the LES-PDF approach, while the stochastic field method was employed for solving the LES-PBE-PDF equations.This approach accounts for the interaction between turbulence, chemistry and particle dynamics.
The LES-PBE-PDF results were compared with the experimental SAXS data of Camenzind et al. (2008) for the S-10 turbulent flame.To the best of the authors' knowledge, this was the first simulation of this experiment.Overall, good agreement was found between numerical predictions and experimental data for time-averaged temperature along the centreline.However, the LES-PBE-PDF model underestimated primary particle sizes by a factor of two.By examining the components of the methodology and leveraging knowledge that has been gained from our laminar flame simulations, the discrepancy was attributed to uncertainties associated with the nucleation model and precursor decomposition kinetics.Furthermore, the number concentration of primary particles and the silica volume fraction were overestimated by the model, indicating that the two-step HMDSO reaction mechanism employed is rather fast, at least for the flame conditions examined here.
The effect of the particle sub-grid scale mixing was investigated by comparing results obtained with two different values of the micromixing constant for number densities.A lower value of the micromixing constant, which reflects the low diffusivity of nanoparticles, resulted in larger values for the silica volume fraction at downstream locations, but experimental data are not available for most of that region.Finally, results obtained with eight stochastic fields were almost identical to those obtained with sixteen stochastic fields, suggesting that eight stochastic fields are sufficient to represent the sub-grid scale effects for both species and number densities.
An investigation of the gas-to-particle conversion processes was also conducted.The LES-PBE-PDF predicted that most of the precursor was fully converted to product species in the region above the burner.Furthermore, the results showed that the contribution of condensation to the total gas-to-particle conversion was more significant than that of nucleation.
Overall, the LES-PBE-PDF method has been shown to be a promising tool for simulating nanoparticle synthesis in turbulent flames.It is, however, clear that further research is required to reach quantitative agreement, particularly with respect to kinetic models.

Fig. 1
Fig. 1 Contour plots of the instantaneous (left) and time-averaged (right) fields for a streamwise velocity, b temperature and c silica volume fraction.HAB stands for height above burner

Fig. 2
Fig. 2 Comparison of numerical predictions with the experimental FTIR data of Camenzind et al. (2008) for flame temperature along the centreline.An average experimental error of ±60 K given by Kammler et al. (2002) was used in the figure for the FTIR data

Fig. 3
Fig. 3 Comparison of the experimental SAXS data of Camenzind et al. (2008) with numerical predictions for a primary particle diameter and b radius of gyration of aggregates along the pathline of maximum silica volume fraction

Fig. 4
Fig. 4 Comparison of the experimental data of Camenzind et al. (2008) with numerical predictions for a primary particle number concentration and b silica volume fraction along the pathline of maximum silica volume fraction.A logarithmic scale was used for the vertical axes

Fig. 5
Fig.5Comparison of numerical predictions obtained with two different values of the micromixing constant for number densities, i.e.C d,p = 0.002 and C d,p = 2 : a flame temperature; b primary particle diameter; c primary particle number concentration; d particle volume fraction.In both cases, the micromixing constant for the gas-phase species and enthalpy was C d,s = 2 .Comparison is also made with the experimental data ofCamenzind et al. (2008)

Fig. 6 Fig. 7 a
Fig.6Comparison of numerical predictions obtained with eight and sixteen stochastic fields: a flame temperature; b primary particle diameter; c primary particle number concentration.For both cases, the micromixing constants were C d,s = 2 and C d,p = 0.002 .Comparison is also made with the experimental data ofCamenzind et al. (2008)

Fig. 8
Fig. 8 Rate of change of silica volume fraction due to nucleation and condensation.a Contour plots; b results along the pathline of maximum silica volume fraction

Table 3
Inlet boundary conditionsu inlet denotes the stream bulk inlet velocity in the streamwise direction