Modeling of silica synthesis in a laminar flame by coupling an extended population balance model with computational fluid dynamics

Abstract

In the present study, we propose a novel extended population balance equation (PBE) model for aggregation and sintering and couple it with computational fluid dynamics (CFD) to investigate synthesis of silica nanoparticles in a laminar diffusion flame. The extended PBE includes finite-rate sintering of primary particles by solving the PBE together with a transport equation for the number concentration of primary particles. In the process simulated, the particles are formed via the oxidation of a vapor precursor, hexamethyldisiloxane (HMDSO), and the aerosol processes include nucleation, condensation, aggregation and sintering. The model is validated with detailed experimental in-situ SAXS data found in the literature and is also compared with a monodisperse and a two-PBE approach. Good agreement is found between the extended one-PBE and two-PBE models, while both of them provide a substantial improvement over the monodisperse one. Furthermore, the coupled CFD-PBE simulation with the extended one-PBE model reduces substantially the computational time as compared with the two-PBE model and requires less than twice the time needed for the monodisperse model. Excellent agreement is found between numerical predictions and experimental data for temperature along the centerline and reasonably good agreement is found between numerical predictions and SAXS data for primary particle diameters. While results for the particle number concentration are in qualitative agreement with the experimental data, the particle formation rate is overpredicted, leading to an overestimation of the number concentration of the primary particles. This is attributed to uncertainties in the experimental data and precursor decomposition kinetics.

EDITOR:

• Mark Swihart

1. Introduction

Aerosol flame synthesis is one of the main methods for the industrial production of nanomaterials and has received increased attention over the last few decades. Flame reactors are simple, scalable and most favorable in terms of energy use (Buesser and Pratsinis 2012). Today, 80–90% of the commercially available flame-made nanomaterials are manufactured in diffusion or spray flame reactors (Wegner and Pratsinis 2003b). The nanoparticles produced find a wide range of applications in gas sensors, energy storage materials, nanotoxicity studies, catalysts, biomaterials, and electroceramics (Meierhofer and Fritsching 2021), among others. It is thus desirable to develop models that can aid the production of tailor-made nanoparticles by predicting the particle size and morphology and linking them to process design and operating conditions.

A number of complex transport and kinetic processes take place in flame synthesis of nanoparticles. Fluid mechanics ‘set the stage’ by establishing the temperature and reactant distribution in the domain. Gas-phase chemical kinetics govern the fuel combustion and the precursor decomposition and oxidation. Aerosol processes include nucleation, condensation, aggregation, and sintering. In spite of extensive experimental and theoretical research over the last few decades, modeling, and simulation of the overall aerosol flame synthesis process remains a challenge (Buesser and Gröhn 2012). Prediction of particle morphology in different flow regimes, such as laminar or turbulent, is a particularly difficult task (Camenzind et al. 2008).

Flame-made nanoparticles are polydisperse in terms of particle sizes, and their polydispersity is characterized by the particle size distribution (PSD). Aerosol dynamics can be described by the population balance equation (PBE), also known as the general dynamic equation (GDE) for aerosols, which describes the spatial and temporal evolution of the PSD due to all the aforementioned physicochemical phenomena encountered in nanoparticle flame synthesis. Furthermore, spatial phenomena such as particle convection, diffusion, and thermophoresis cause particles to follow different trajectories with different temperature and species-concentration histories, leading to a product with a wide variety in morphology. Therefore, the comprehensive modeling of flame synthesis requires coupling of PBE with fluid dynamics, species transport and chemical kinetics to predict nanoparticle properties in laminar and turbulent flows. The role of CFD-PBE simulations is twofold: first, they allow comparison with experimental data to reveal sources of uncertainties and, thus, guide theoretical and experimental research, and second, they can predict the particle properties (primary particle and aggregate sizes, number concentrations, etc.) in the whole reactor domain and thus facilitate the design and scale-up of reactors.

A number of CFD-based studies of flame synthesis of nanoparticles have been reported in the literature. In the present work, the proposed methodology is applied to silica synthesis; hence we survey the relevant studies in chronological order, focusing on the modeling approach adopted in each case. For studies on other inorganic nanoparticles, we refer to the review articles of Buesser and Pratsinis (2012); Buesser and Gröhn (2012), Meierhofer and Fritsching (2021) and Raman and Fox (2016). Furthermore, a review of experimental studies including in-situ measurements of silica nanoparticle synthesis in flames has been conducted, and the studies are summarized in Table 1.

Table 1. Comparison of reported flame-made silica characteristics found in selected studies in the literature measured with in-situ techniques.

Reference	Experimental technique	Synthesis method	Precursor	$d_p$ [nm]	$R_g$ [nm]	$D_f$
Scheckman, McMurry, and Pratsinis (2009)	Mass-Mobility	DF	HMDSO	$12-93$	n/a	$1.7-2.4$
Camenzind et al. (2008) ‘S-2’	SAXS + TS	LF DF	HMDSO	$5-90$	$52-121$	$1.7-2.9$
Camenzind et al. (2008) ‘S-10’	SAXS + TS	TF DF	HMDSO	$9-12$	$45-93$	$1.75-2.13$
Kammler et al. (2005) ‘Cold’	USAXS	TF PM	HMDSO	$5-11$	$65-135$	$2.1-2.35$
Kammler et al. (2005) ‘Hot’	USAXS	TF PM	HMDSO	$8-27$	$40-170$	$1.95-2.58$
Kim and Choi (2003)	plane laser LS	flat PM	SiCl4	n/a	$50-220$	$1.3-2$
Cho and Choi (2000)	LS + TS	$H_2$ DF	SiCl4	$15-50$	$25-140$	$1.62-1.76$
Choi et al. (1999)	LS + TS	$H_2$ DF	SiCl4	$19-49$	$55-110$	$1.35-1.85$
Chang and Biswas (1992)	LS	DF	SiCl4	n/a	$292-303$	$2.17-2.8$
Chung, Sheu, and Tsai (1992)	laser LS	CF DF	SiCl4	$20-30$	$50-200$	n/a
Zachariah et al. (1989)	dynamic LS	CF DF	SiH4	$20-30$	n/a	n/a
Ulrich and Rieh (1982)	laser LS	PM	SiCl4	$7.6-18$	n/a	n/a

The data presented in the table correspond to: experimental technique (SAXS: small-angle X-ray scattering, USAXS: untra small-angle X-ray scattering, TS: thermophoretic sampling, LS: light scattering), method of synthesis (LF: laminar flow, TF: turbulent flow, CF: counterflow, DF: diffusion flame, PM: premixed flame), material used as a precursor, primary particle $d_p$ and aggregate $R_g$ sizes, fractal dimension $D_f.$

The coupling of aerosol and fluid dynamics for modeling silica nanoparticle formation in laminar flows was explored in the studies of Kim and Pratsinis (1988), Pratsinis and Kim (1989) and Kim and Pratsinis (1989). A method of moments, where the PSD was assumed to adopt a lognormal shape, was employed to describe aerosol dynamics, while only spherical particles were considered.

Lee, Oh, and Choi (2001) investigated flame synthesis of silica non-spherical nanoparticles in a laminar flat flame. A one-dimensional flame model was employed, together with the sectional two-dimensional PBE model of Xiong and Pratsinis (1993), which accounts for the evolution of aggregate structure. Different expressions for the sintering characteristic time were tested, including a hybrid model between viscous flow and atomistic diffusion sintering models (Ehrman 1999). Numerical predictions for number concentrations and particle volume fractions were compared with the in-situ experimental data of Chang and Biswas (1992). Later, Kim et al. (2003) performed two-dimensional CFD simulations to investigate flame synthesis of silica nanoparticles by oxidation of silicon tetrachloride (SiCl₄) in a co-flow diffusion laminar flame. A two-PBE sectional model was coupled with CFD and the hybrid sintering model was employed. Their model slightly underestimated the primary particle sizes.

Subsequently, Ji et al. (2007) investigated flame synthesis of silica nanoparticles via CFD. They implemented a moment model and coupled it with commercial CFD software to simulate synthesis of silica nanoparticles through oxidation of tetraethylorthosilicate (TEOS) in a hydrogen diffusion flame reactor. Good agreement with the experimental ex-situ data of Jang (2001) was achieved by performing a parametric analysis and adjusting the particle formation and growth rate constants. However, their model did not distinguish between primary particles and aggregates. Buddhiraju and Runkana (2012) advanced the analysis of Ji et al. (2007) by using a monodisperse population balance model that accounts for fractal aggregates. Results were compared again with the ex-situ experimental data of Jang (2001) and were in close agreement.

Gröhn et al. (2011) incorporated a monodisperse moment model that accounts for fractal aggregates within CFD. They used the Reynolds-averaged Navier-Stokes equations (RANS) approach with a realizable k-$ϵ$ turbulence model to model flame synthesis of silica nanoparticles through oxidation of hexamethyldisiloxane (HMDSO) vapor in turbulent methane/oxygen diffusion flames. Furthermore, they investigated different sintering-time expressions, and the results were compared with experimental ex-situ particle size data. Their model underpredicted the primary particle sizes at high oxygen flow rates and at high silica production rates. Olivas-Martinez et al. (2015) used a commercial CFD software to model flame spray pyrolysis (FSP) for synthesis of silica nanopowder and demonstrated the capability of CFD models to be used in the design process of FSP reactors. The standard k-$ϵ$ turbulence model was employed to describe fluid dynamics, spray droplet evolution was described from the Lagrangian point of view and the quadrature method of moments (QMOM) was used to solve the PBE. Their model incorporated nucleation and coagulation of nanoparticles. However, finite-rate sintering was not considered and, therefore, primary particle and aggregate sizes were not distinguished. Instead, the particles were assumed to be spherical, and a finite value of a sticking coefficient was incorporated in the coagulation kernel.

Vo et al. (2017) developed a sparse-Lagrangian multiple mapping conditioning (MMC) model in the context of large eddy simulation (LES), thus accounting for turbulence-chemistry interaction. This model was applied to describe silica nucleation in a temporal mixing layer and the results were validated against data obtained via direct numerical simulations (DNS). Subsequently, this approach was employed by Neuber et al. (2019) who presented a joint experimental and numerical study of flame synthesis of silica nanoparticles. They showed the modeling capabilities of the PBE-MMC-LES approach by simulating flame synthesis of silica nanoparticles in a turbulent reacting jet. This approach was based on a sectional model for aerosol dynamics and also considered the interaction between turbulence, chemistry and particle dynamics. Aggregation of non-spherical particles was considered by assuming a constant primary particle diameter of 0.98 nm. To an extent, a fair agreement was found with the experimental data, while discrepancies were attributed to uncertainties in the precursor-decomposition reaction mechanism.

Feroughi et al. (2017) studied experimentally and numerically a premixed laminar low-pressure flame seeded with low concentrations of HMDSO, and they proposed a reaction scheme for the decomposition and combustion of HMDSO. This scheme was incorporated, along with a monodisperse moment model, into CFD to simulate flame synthesis of nanoparticles. A fair agreement between numerical predictions and experimental particle sizes was found, while the model slightly underpredicted the primary particle sizes at high precursor concentrations. The proposed reaction mechanism was utilized by (Rittler et al. 2017) where they developed and presented a model within the LES framework for describing flame spray pyrolysis for synthesis of nanoparticles; however, a comparison with experimental data was not attempted. Recently, Dasgupta et al. (2022) used a commercial CFD software to perform unsteady RANS and study silica nanoparticle synthesis from flame spray pyrolysis. Nanoparticle dynamics was described with the hybrid method of moments while the model of Tsantilis, Briesen, and Pratsinis (2001) for sintering time was employed. Different operating conditions of a complex geometry reactor were investigated, and mean particle diameters were compared with ex-situ experimental data.

Several points can be concluded from the preceding literature review. First, there is a need to develop models capable of predicting the change in particle morphology at different flow regimes. Second, most of the reported CFD studies on turbulent flames were within the context of RANS, which adds extra uncertainty due to the need for modeling turbulent flow, turbulence-chemistry interaction and turbulence-particle interaction where extra terms arise from the Reynolds averaging of the PBE (Rigopoulos 2007; Tsagkaridis, Rigopoulos, and Papadakis 2022). The modeling challenges arising in turbulent-flow simulations are discussed in more detail in the review articles of Rigopoulos (2010) and Raman and Fox (2016). The extra uncertainty introduced by turbulence closure models could lead to compensation of modeling errors and this complicates the efforts in identifying the precise sources of discrepancies when comparing with experiments. Furthermore, in most of the studies, a monodisperse or moment method was employed for the solution of PBE, which involves further modeling assumptions and does not predict the evolution of the PSD directly. Finally, in most cases, the validation of previous CFD models was done by comparing final primary particle sizes with ex-situ experimental data for particles, while comparison with in-situ data focusing on the evolution of particle growth has been overlooked.

In light of the aforementioned points, the purpose of the present article is twofold. First, we propose a novel extended population balance model for simulating aerosol flame synthesis with the following key characteristics: it distinguishes between primary particles and aggregates, accounts for finite-rate sintering and is computationally efficient. To accomplish this, a sectional method is employed to solve the PBE, which is solved together with a transport equation for the number concentration of primary particles. This approach is more computationally efficient (especially with respect to coupling with CFD) than a two-dimensional or a two-equation PBE while still being able to predict the evolution of primary particle sizes. Second, the aforementioned PBE is coupled with a detailed CFD model and employed to predict the detailed experimental in-situ small-angle X-ray scattering (SAXS) measurements of Camenzind et al. (2008) in a laminar silica synthesis diffusion flame. The quality of CFD simulations is subject to the uncertainties and assumptions involved in the mathematical model. Laminar flames do not pose the complications of turbulence and turbulence-kinetics interaction, and furthermore, the sectional approach employed is free of the closure models required by monodisperse and moment methods. By employing an approach with minimum assumptions with respect to fluid dynamics and population balance modeling, we aim to limit the modeling uncertainties to the selection of aerosol kinetic models and, therefore, shed light on the nature of discrepancies between predictions and experimental results.

The rest of the article is organized as follows: First, the governing equations and population balance models employed in this study are introduced. The subsequent sections describe the numerical solution method, the kinetic models employed and the flame configuration and computational setup. Results are then presented and discussed, and the article concludes with a summary of the main findings.

2. Mathematical formulation

2.1. Reacting flow

The variable density and low Mach number formulation of the conservation equations for mass, momentum, species and energy is employed and summarized below. The equations are presented with Cartesian tensor notation. The continuity equation is: (1) $$\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u_i)}{\partial x_i}=0\hfill \end{array}$$(1)

The momentum equation is (2) $$\begin{array}{c}\frac{\partial (\rho u_i)}{\partial t}+\frac{\partial (\rho u_j{u}_i)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial {\sigma }_{ij}}{\partial x_j}+\rho g_i\hfill \end{array}$$(2) where $u_i$ is the velocity component in the i-th direction, p is the pressure, $\rho$ is the density and $g_i$ is the gravitational acceleration. For a Newtonian fluid, the stress tensor is given by ${\sigma }_{ij}=2\mu \left(S_{ij}-\frac{1}{3}{S}_{kk} {\delta }_{ij}\right),$ where $\mu$ is the dynamic viscosity of the gas mixture, ${\delta }_{ij}$ is the Kronecker delta, and the strain rate is defined by $S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$ The species transport equation is (3) $$\begin{array}{c}\frac{\partial (\rho Y_{\alpha })}{\partial t}+\frac{\partial (\rho u_j{Y}_{\alpha })}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\rho D_{\alpha }\frac{\partial Y_{\alpha }}{\partial x_j}\right)+\rho \hspace{0.17em}{\dot{\omega }}_{\alpha }\hfill \end{array}$$(3) where $Y_{\alpha },$ $D_{\alpha }$ are the mass fraction and diffusion coefficient of a chemical species into the gas mixture, respectively,¹ while the mass source due to chemical reaction is denoted ${\dot{\omega }}_{\alpha }.$ Fick’s diffusion law and the Hirschfelder and Curtiss approximation (Hirschfelder, Curtiss, and Bird 1955; Poinsot and Veynante 2005) have been employed to calculate $D_{\alpha }.$ A correction velocity is employed to ensure mass conservation (Poinsot and Veynante 2005). The species transport coefficients are estimated with relationships derived from the kinetic theory of gases. The dynamic viscosity of the mixture is obtained by Wilke’s mixing rule (Wilke 1950). Finally, the energy equation is written in terms of specific enthalpy that includes the thermal and chemical enthalpies, as follows: (4) $$\begin{array}{c}\frac{\partial (\rho h)}{\partial t}+\frac{\partial (\rho u_{j}h)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\frac{k}{c_p}\frac{\partial h}{\partial x_j}\right)-\sum _{\alpha =1}^M{h}_{\alpha }\frac{\partial }{\partial x_j}\left[\rho \left(\frac{k}{\rho c_p}-D_{\alpha }\right)\frac{\partial Y_{\alpha }}{\partial x_j}\right]+\rho \hspace{0.17em}{\dot{q}}_{\mathit{rad}}+\rho \hspace{0.17em}{\dot{q}}_{\mathit{rad},p}\hfill \end{array}$$(4) where ${\dot{q}}_{\mathit{rad}}$ and ${\dot{q}}_{\mathit{rad},p}$ denote enthalpy losses due to the thermal radiation from gas species and particles, respectively, and M is the number of chemical species. Since both chemical and thermal enthalpies are included in the transported quantity, the interchange between them due to reaction does not appear explicitly. The thermal conductivity, k, is obtained by the averaging formula of Mathur and Saxena (1967). Viscous heating, Soret and Dufour effects have been neglected in the overall formulation.

2.2. Population balance

A population balance equation (PBE) is employed to describe aerosol dynamics. Particle morphology is an important factor and its prediction is one of the objectives of the present article. A two-dimensional formulation (Koch and Friedlander 1990; Xiong and Pratsinis 1993) introduces an additional dimension such as particle surface area, but its Discretized form is too expensive for coupling with CFD. A two-PBE formulation employs two sets of PBEs, one for the number density of particles and another one for number of primary particles (Rogak 1997) or particle surface area (Tsantilis and Pratsinis 2000); such an approach is computationally feasible and able to incorporate morphology effects, but still requires twice the number of equations as compared with the one-PBE approach. A monodisperse model (Kruis et al. 1993) has also been proposed, although this entails several assumptions. In the present section, we first describe the monodisperse and two-PBE approaches employed in our study. Subsequently, a new formulation is proposed that, while employing one PBE, incorporates elements of the two-PBE approach with respect to the description of morphology.

Before proceeding, we first summarize the description of particle morphology adopted in the present study. Coalescence and sintering give rise to fractal aggregates, which are clusters of primary particles. The aggregation kernel must, therefore, be formulated in terms of the aggregate collision diameter, $d_g.$ Fractal-like aggregates follow a power-law relationship that relates their collision diameter with the number of primary particles per aggregate, $n_p,$ and their average diameter, $d_p,$ and assumes the following form (Friedlander 2000; Buesser and Pratsinis 2012): (5) $$\begin{array}{c}{d}_g=d_p{\left(\frac{n_p}{k_f}\right)}^{1/D_f}\hfill \end{array}$$(5) where $k_f$ is a prefactor associated with the aggregate anisotropy (Heinson, Sorensen, and Chakrabarti 2010) and $D_f$ is the fractal dimension, which has values close to 3 for compacted aggregates and close to 1 for chain-like structures (Friedlander 2000). The evolution of $D_f$ has little effect on final particle characteristics (Goudeli, Eggersdorfer, and Pratsinis 2015). In the present simulations, the values $k_f=1.4$ and $D_f=1.91$ were employed corresponding to ballistic cluster–cluster aggregation in the free molecule regime (Ball and Jullien 1984; Friedlander 2000).

2.2.1. The monodisperse model

The simplest and probably the most widely used model for describing aerosol dynamics in flame-synthesis studies is the monodisperse model of Kruis et al. (1993). In this model, only transport equations for the moments of the PSD are solved, while both aggregate and primary particle size distributions are assumed to be monodisperse. Its success in describing relatively well average particle properties can be attributed to some unique features of nanoparticle flame synthesis, in particular the rapid attainment of the self-preserving aerosol size distribution (SPSD) by coagulation, which simplifies the modeling strategy Buesser and Pratsinis (2012). While the monodispersed model was initially formulated in terms of surface area, in the present work we reformulate it in terms of the number concentration of primary particles for consistency with the detailed PBE models to be discussed later. Specifically, transport equations for the number concentration of aggregates, N and particle volume fraction, and V, are solved along with a transport equation for the number concentration of primary particles, $N_p.$ This formulation is equivalent to that of (Kruis et al. 1993; Gröhn et al. 2011) and the model implementation is validated with the data of Goudeli, Eggersdorfer, and Pratsinis (2015) (Figure S1 of the supplementary information). The equations to be solved are: (6) $$\begin{array}{c}\frac{\partial N}{\partial t}+\frac{\partial (U_{j}N)}{\partial x_j}-\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial N}{\partial x_j}\right)=J(\mathbf{Y},T)-\frac{1}{2}\beta (v_m,v_m)N^2\hfill \end{array}$$(6) (7) $$\begin{array}{c}\frac{\partial V}{\partial t}+\frac{\partial (U_{j}V)}{\partial x_j}-\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial V}{\partial x_j}\right)=J(\mathbf{Y},T)v_{\mathit{nuc}}+\beta (v_m,v_o)N\frac{p_{{\text{SiO}}_2}}{k_{B}T}{v}_o\hfill \end{array}$$(7) (8) $$\begin{array}{c}\frac{\partial N_p}{\partial t}+\frac{\partial (U_j{N}_p)}{\partial x_j}-\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial N_p}{\partial x_j}\right)=J(\mathbf{Y},T)-\frac{3}{{\tau }_s}\left(\frac{V}{v_{p,av}}-\frac{M_{23}}{v_{p,av}{}^{2/3}}\right)\hfill \end{array}$$(8) where $v_o$ is the volume of a SiO₂ monomer, the particle mean size is given by $v_m=V/N,$ $v_{\mathit{nuc}}$ is the volume of a nucleus particle, $D_p$ is the diffusion coefficient of particles of volume v estimated by the Stokes-Einstein equation (Friedlander 2000), J is the nucleation rate, $\beta$ is the aggregation kernel and ${\tau }_s$ is the sintering characteristic time. The kinetic expressions employed will be shown in Section 4. The last term in Equation (7)(7) $$\begin{array}{c}\frac{\partial V}{\partial t}+\frac{\partial (U_{j}V)}{\partial x_j}-\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial V}{\partial x_j}\right)=J(\mathbf{Y},T)v_{\mathit{nuc}}+\beta (v_m,v_o)N\frac{p_{{\text{SiO}}_2}}{k_{B}T}{v}_o\hfill \end{array}$$(7) represents the condensation process, where $p_{{\text{SiO}}_2}$ is the partial pressure of the SiO₂ gas species and $k_b$ is the Boltzmann constant. In the monodisperse model, the $2/3$ fractional moment of the PSD is given by $M_{23}=v_m^{2/3}N=V^{2/3}{N}^{1/3}.$ The velocity, $U_j,$ is equal to the sum of the flow velocity, $u_j,$ and the thermophoretic velocity, $u_{j,th},$ and the latter can be estimated from the following relationship (Friedlander 2000): (9) $$\begin{array}{c}{u}_{j,th}=-0.55\frac{\mu }{\rho \hspace{0.17em}T}\frac{\partial T}{\partial x_j}\hfill \end{array}$$(9) where T is the local temperature of the gas mixture. The average primary particle size, $v_{p,av},$ and the average number of primary particles per aggregate, $n_{p,av},$ can be calculated as follows: (10) $$\begin{array}{c}{v}_{p,av}=\frac{V}{N_p}\hfill \end{array}$$(10) (11) $$\begin{array}{c}{n}_{p,av}=\frac{N_p}{N}\hfill \end{array}$$(11)

Figure 1. Experimental setup for the synthesis of the SiO₂ nanoparticles with a co-flow diffusion flame burner. Reprinted from ‘Nanostructure evolution: from aggregated to spherical SiO₂ particles made in diffusion flames’, by A. Camenzind, H. Schulz, A. Teleki, G. Beaucage, T. Narayanan and S.E. Pratsinis, European Journal of Inorganic Chemistry, 2008, vol. 6, p. 911–918. Copyright Wiley-VCH GmbH. Reproduced with permission.

2.2.2. The two-PBE model for aggregation and sintering

A more comprehensive approach to the prediction of the PSD constitutes the direct use of the PBE and its solution by Discretization (the so-called sectional approach). Three main augmentations to the basic PBE have been proposed to account for particle morphology: the two-dimensional PBE, the modified one-dimensional PBE and the coupled two-PBE approach. The two-dimensional approach is too expensive for coupling with a CFD simulation; hence CFD-based efforts have employed the other two methods, a comparison between which can be found in Sun, Rigopoulos, and Liu (2021) in the context of soot formation. In the present work, the two-PBE approach will be used, and its equations are briefly summarized below. The PBE for the aggregates is (Friedlander 2000): (12) $$\begin{array}{c}\frac{\partial n}{\partial t}+\frac{\partial (U_{j}n)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C(\mathbf{Y},v)\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)dw\hfill \end{array}$$(12) where v and w express volume of aggregates, while $n=n(\mathbf{x},t,v)$ is the particle number density concentration (denoting number of particles per unit of particle volume and mixture volume). B and C denote the nucleation² and condensation rates, respectively, while $\delta (v-v_{\mathit{nuc}})$ is the Dirac delta function. In the two-PBE approach, an additional PBE must be solved to account for particle morphology. In the present formulation, this is an equation for the number density of primary particles and takes the following form (Sun, Rigopoulos, and Liu 2021): (13) $$\begin{array}{c}\frac{\partial n_p}{\partial t}+\frac{\partial (U_j{n}_p)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(n_{av}{D}_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C_p(\mathbf{Y},v)-\frac{3}{{\tau }_s}\left[\frac{n_p}{n}-{\left(\frac{n_p}{n}\right)}^{2/3}\right]\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)(n_{av}(w)+n_{av}(v-w))dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)n_{av}(v)dw\hfill \end{array}$$(13) where $n_p=n_p(\mathbf{x},t,v)$ is the number density function of primary particles that form an aggregate of volume v per unit volume, $n_{av}(\mathbf{x},t,v)$ is the average number of primary particles per aggregate of volume v and $C_p$ is the condensation rate on primary particles. By definition, the $n_{av}$ is given by: (14) $$\begin{array}{c}{n}_{av}(v,t)=\frac{n_p(v,t)}{n(v,t)}\hfill \end{array}$$(14)

Note that, unlike n and $n_p,$ $n_{av}$ is not a distribution. It is also worth mentioning that the integration of Equations ((12)(12) $$\begin{array}{c}\frac{\partial n}{\partial t}+\frac{\partial (U_{j}n)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C(\mathbf{Y},v)\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)dw\hfill \end{array}$$(12) and (13))(13) $$\begin{array}{c}\frac{\partial n_p}{\partial t}+\frac{\partial (U_j{n}_p)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(n_{av}{D}_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C_p(\mathbf{Y},v)-\frac{3}{{\tau }_s}\left[\frac{n_p}{n}-{\left(\frac{n_p}{n}\right)}^{2/3}\right]\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)(n_{av}(w)+n_{av}(v-w))dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)n_{av}(v)dw\hfill \end{array}$$(13) over v results in transport equations for the number concentration of aggregates, N, and the number concentration of primary particles, $N_p,$ which are the variables solved for in the monodisperse model discussed earlier.

Sintering of incipient particles smaller than a minimum particle diameter can be considered instantaneous (Tsantilis, Briesen, and Pratsinis 2001). A critical diameter $d_{p,\min }$ that distinguishes particles that coalesce instantaneously from those that fuze with a finite rate can be introduced in the two-PBE model. This way, particles smaller than the critical diameter are assumed to be spherical, while the number of equations to be solved in Discretization methods is reduced since $n_p(v<v_{p,\min })\equiv n(v<v_{p,\min }),$ where $v_{p,\min }=\pi d_{p,\min }^3/6.$ This value was selected to be 1 nm in the present study.

From knowledge of $n(v)$ and $n_p(v),$ we estimate the number concentration of primary particles, $N_p,$ the average primary particle diameter, $d_{p,av},$ the average number of primary particles per aggregate, $n_{p,av},$ and the average radius of gyration of aggregates as follows: (15) $$\begin{array}{cc}{N}_p\hfill & =\underset{v_d}{\overset{\infty }{\int }}{n}_p(v)dv\hfill \end{array}$$(15) (16) $$\begin{array}{c}{d}_{p,av}=\frac{\underset{v_d}{\overset{\infty }{\int }}{d}_p(v)n_p(v)dv}{N_p}\hfill \end{array}$$(16) (17) $$\begin{array}{cc}{n}_{p,av}\hfill & =\frac{N_p}{N}\hfill \end{array}$$(17) (18) $$\begin{array}{c}{R}_{g,av}=\frac{1}{2}\frac{\underset{v_d}{\overset{\infty }{\int }}{d}_g(v)n(v)dv}{N}\hfill \end{array}$$(18) where $d_p(v)$ is the diameter of the primary particles in an aggregate of volume v, $d_g(v)$ is the aggregate diameter given by Equation (5)(5) $$\begin{array}{c}{d}_g=d_p{\left(\frac{n_p}{k_f}\right)}^{1/D_f}\hfill \end{array}$$(5) , and $v_d$ is the volume of the smallest detectable particles (1 nm) by SAXS measurements (Camenzind et al. 2008).

The assumption inherent in the two-PBE model is that all aggregates of particle volume v consist of equally-sized spherical primary particles. While this is not strictly true, similar assumptions have been used in experimental techniques (Iyer et al. 2007), while transmission electron microscope (TEM) images of aggregates have repeatedly shown that the distribution of primary particles within an aggregate is relatively uniform (Kammler et al. 2005), supporting the assumptions adopted here. The two-PBE model is a comprehensive approach that yields information for both aggregate and primary particle size distributions and their evolution at a reduced cost as compared with a two-dimensional PBE, but the number of equations to be solved is still twice that of a one-PBE approach, which is significant when it comes to coupling with CFD.

2.2.3. A Novel extended one-PBE model for aggregation and sintering

In the present section, we propose an extended one-PBE model that captures aggregation dynamics and can describe particle-morphology evolution while maintaining computational cost to a level commensurable with other one-PBE approaches. One-PBE models have been employed in the modeling of soot formation in laminar and turbulent flows by simply introducing a predefined cutoff size (Netzell, Lehtiniemi, and Mauss 2007; Rodrigues et al. 2018) separating primary particles and aggregates that corresponds to the average primary particle size observed experimentally. This approach does not include a finite-rate sintering model, and a detailed discussion of its shortcomings can be found in Sun, Rigopoulos, and Liu (2021). Here, we overcome this problem by solving the PBE together with a transport equation for the number concentration of primary particles, $N_p.$ In this way, the sintering process is incorporated into the model and the one-PBE model can predict the evolution of primary particle size.

The evolution of the PSD of aggregates is described by the same equation (Equation (12)(12) $$\begin{array}{c}\frac{\partial n}{\partial t}+\frac{\partial (U_{j}n)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C(\mathbf{Y},v)\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)dw\hfill \end{array}$$(12) ) as in the two-PBE model. To obtain the additional equation, we assume that the primary particle size is monodisperse (i.e., it is the same for aggregates of all sizes, unlike in the case of the two-PBE model). Based on this assumption, $n_{av}$ can be expressed as: (19) $$\begin{array}{c}{n}_{av}(v,t)=\frac{v}{v_{p,av}(t)}\hfill \end{array}$$(19) where v is the volume of an aggregate and $v_{p,av}=\pi d_{p,av}^3/6$ is the average primary particle volume. By substituting Equation (19)(19) $$\begin{array}{c}{n}_{av}(v,t)=\frac{v}{v_{p,av}(t)}\hfill \end{array}$$(19) into Equation (13)(13) $$\begin{array}{c}\frac{\partial n_p}{\partial t}+\frac{\partial (U_j{n}_p)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(n_{av}{D}_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C_p(\mathbf{Y},v)-\frac{3}{{\tau }_s}\left[\frac{n_p}{n}-{\left(\frac{n_p}{n}\right)}^{2/3}\right]\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)(n_{av}(w)+n_{av}(v-w))dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)n_{av}(v)dw\hfill \end{array}$$(13) and integrating both sides of (13)(13) $$\begin{array}{c}\frac{\partial n_p}{\partial t}+\frac{\partial (U_j{n}_p)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(n_{av}{D}_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C_p(\mathbf{Y},v)-\frac{3}{{\tau }_s}\left[\frac{n_p}{n}-{\left(\frac{n_p}{n}\right)}^{2/3}\right]\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)(n_{av}(w)+n_{av}(v-w))dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)n_{av}(v)dw\hfill \end{array}$$(13) over v, we derive a transport equation for the number concentration of primary particles: (20) $$\begin{array}{c}\frac{\partial N_p}{\partial t}+\frac{\partial (U_j{N}_p)}{\partial x_j}-\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial N_p}{\partial x_j}\right)=J(\mathbf{Y},T)-\frac{3}{{\tau }_s}\left(\frac{V}{v_{p,av}}-\frac{M_{23}}{v_{p,av}{}^{2/3}}\right)\hfill \end{array}$$(20) where V the particle volume fraction and $M_{23}$ is the $2/3$ fractional moment of the PSD. The condensation term is absent since it does not alter the number concentration of primary particles, while the aggregation terms cancel out after the integration of (13)(13) $$\begin{array}{c}\frac{\partial n_p}{\partial t}+\frac{\partial (U_j{n}_p)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(n_{av}{D}_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C_p(\mathbf{Y},v)-\frac{3}{{\tau }_s}\left[\frac{n_p}{n}-{\left(\frac{n_p}{n}\right)}^{2/3}\right]\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)(n_{av}(w)+n_{av}(v-w))dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)n_{av}(v)dw\hfill \end{array}$$(13) over v. Equation (20)(20) $$\begin{array}{c}\frac{\partial N_p}{\partial t}+\frac{\partial (U_j{N}_p)}{\partial x_j}-\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial N_p}{\partial x_j}\right)=J(\mathbf{Y},T)-\frac{3}{{\tau }_s}\left(\frac{V}{v_{p,av}}-\frac{M_{23}}{v_{p,av}{}^{2/3}}\right)\hfill \end{array}$$(20) is solved along with (12)(12) $$\begin{array}{c}\frac{\partial n}{\partial t}+\frac{\partial (U_{j}n)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C(\mathbf{Y},v)\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)dw\hfill \end{array}$$(12) but, unlike the second equation in the two-PBE model, does not require Discretization in the particle volume domain. The $v_{p,av}$ and $n_{p,av}$ are then estimated by Equations ((10)(10) $$\begin{array}{c}{v}_{p,av}=\frac{V}{N_p}\hfill \end{array}$$(10) and (11))(11) $$\begin{array}{c}{n}_{p,av}=\frac{N_p}{N}\hfill \end{array}$$(11) . The major advantage of Equation (20)(20) $$\begin{array}{c}\frac{\partial N_p}{\partial t}+\frac{\partial (U_j{N}_p)}{\partial x_j}-\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial N_p}{\partial x_j}\right)=J(\mathbf{Y},T)-\frac{3}{{\tau }_s}\left(\frac{V}{v_{p,av}}-\frac{M_{23}}{v_{p,av}{}^{2/3}}\right)\hfill \end{array}$$(20) over the conventional one-PBE model is that it brings a finite-rate sintering model into the PBE. Compared with the two-PBE model, the drawback is the assumption of monodisperse primary particle size. The importance of this assumption will be discussed in Section 6, where results from both models will be compared with experimental data.

3. Numerical methods

A sectional method is employed for the numerical solution of the PBE (Liu and Rigopoulos 2019). It is a conservative finite volume method that was recently extended to the two-PBE formulation by Sun, Rigopoulos, and Liu (2021). The method provides an accurate prediction of the distribution with a small number of sections while conserving the first moment (with respect to volume) during the aggregation process, which is important for satisfying the mass balance. The method can employ an arbitrary non-uniform grid for Discretization along the particle volume domain and has also been shown to be robust and efficient when incorporated within a CFD framework (Liu and Rigopoulos 2019; Sun, Rigopoulos, and Liu 2021; Sun and Rigopoulos 2022). This approach is used for both the one- and two-PBE formulations.

The PBE models have been implemented in our in-house CFD code BOFFIN, which is based on the variable-density low Mach number formulation of the Navier-Stokes equations and has been used to simulate laminar and turbulent flames with soot formation that have similar features to the problem considered here (Sun, Rigopoulos, and Liu 2021; Sun and Rigopoulos 2022). The Navier-Stokes equations are discretized by applying the finite volume method on a staggered grid. A second-order central Discretization scheme is employed for the convection and viscous terms, while a second-order accurate Crank-Nicolson method is employed for the temporal integration. The convection and viscous terms are treated implicitly. The pressure-velocity coupling is dealt with an iterative SIMPLE-type scheme. The transport equations for the reactive scalars and Discretized PBE are solved with a fractional step method, where separate fractional steps were employed for convection-diffusion, chemical reaction, and PBE integration. A backward Euler scheme is used for the integration of the chemical reaction step, while the explicit Euler scheme was selected for the PBE step since the Discretized PBE is not stiff. The convection-diffusion term is solved as described above with the addition that van Leer’s TVD scheme was employed for the convection terms to ensure boundedness for the reactive scalars.

4. Gas-phase and aerosol kinetics

4.1. Gas-phase kinetics and radiation

The reaction mechanism for methane combustion used in the present study is the Gas Research Institute (GRI) 1.2 mechanism (Frenklach et al. 1995), containing 175 steps and 31 species. The reduced two-step mechanism proposed by Feroughi et al. (2017) is employed for the HMDSO (precursor) decomposition and oxidation and is shown in Table 2. This global two-step HMDSO mechanism has been used previously in the literature to study the synthesis of silica nanoparticles via flame spray pyrolysis (Rittler et al. 2017).

Table 2. Global two-step reaction mechanism for HMDSO combustion (Feroughi et al. 2017).

Step	Reaction	A	n	$E_a$
$G_1$	HMDSO + OH → 2SiO + 6CH3 + H	$6\times {10}^{12}$	0.46	15000
$G_2$	SiO + H2O → SiO2(g) + H2	$8.5\times {10}^{10}$	0	5650

The reaction rate constant is given in the form $k=A{T}^n\hspace{0.17em}\exp [-E_a/(RT)],$ where $E_a$ stands for the activation energy. The units of the parameter A are expressed in terms of cm, s, and mol, while $E_a$ is given in cal/mol.

Following Barlow et al. (2001), the optically thin hypothesis is invoked in the present study to model thermal radiation from gas species. The radiative loss term in the enthalpy h equation due to the radiation of the most abundant species H₂O, CO₂, CH₄, and CO can be computed according to: (21) $$\begin{array}{c}{\dot{q}}_{\mathit{rad}}=-\frac{4\sigma }{\rho }(T^4-T_b{}^4)\sum _i{a}_{p,i}{p}_i\hfill \end{array}$$(21) where $\sigma =5.669·{10}^{-8}$ W/${\mathrm{m}}^2/{\mathrm{K}}^4$ is the Stefan-Boltzmann constant, $T_b=295$ K is the ambient background temperature, i-index indicates summation over the four major radiating species, $p_i$ is the partial pressure of species i and $a_{p,i}$ is the Planck mean absorption coefficients expressed as polynomials in temperature (Barlow et al. 2001; Grosshandler, 1993).

4.2. Nucleation

We follow an approach similar to that described in Shekar et al. (2012) to describe the nucleation process. We assume that two molecules (monomers of SiO₂) in the gas phase collide to form a dimer. The dimers are considered to be the first primary particles in the simulation. The free molecule regime kernel is used to describe the self-collision frequency of the monomers. Based on these considerations, the nucleation rate is given by: (22) $$\begin{array}{c}J=\frac{1}{2}{K}_{fm}{N}_A^2{C}_{{\text{SiO}}_2}{}^2\hfill \end{array}$$(22) where $N_A$ is Avogadro’s number, $C_{{\text{SiO}}_2}$ is the gas-phase concentration of SiO₂ monomers and $K_{fm}$ is the free molecule coagulation kernel given by: (23) $$\begin{array}{c}{K}_{fm}=E_F·4\sqrt{\frac{\pi k_{b}T}{m_{{\text{SiO}}_2}}}{d}_o{}^2\hfill \end{array}$$(23) where $E_F$ is the Van der Waals enhancement factor whose value is approximated to 1.88 based on the observations of Kennedy and Harris (1990), $k_b$ is the Boltzmann constant, and $m_{{\text{SiO}}_2}$ and $d_o$ are the mass and diameter of a SiO₂(g) molecule, respectively. It is assumed that SiO₂(g) molecules have diameter $d_o=0.44$ nm and that the density of silica particles is ${\rho }_p=2196$ kg/${\mathrm{m}}^3.$ After Discretization of the PBE, the nucleation source term in Equation (13)(13) $$\begin{array}{c}\frac{\partial n_p}{\partial t}+\frac{\partial (U_j{n}_p)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(n_{av}{D}_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C_p(\mathbf{Y},v)-\frac{3}{{\tau }_s}\left[\frac{n_p}{n}-{\left(\frac{n_p}{n}\right)}^{2/3}\right]\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)(n_{av}(w)+n_{av}(v-w))dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)n_{av}(v)dw\hfill \end{array}$$(13) takes the form: (24) $$\begin{array}{c}B=J \frac{v_{\mathit{nuc}}}{d{v}_{\mathit{nuc}} v_{m,\mathit{nuc}}}\hfill \end{array}$$(24) where $v_{\mathit{nuc}}$ is the volume of dimers (volume of nuclei), $v_{m,\mathit{nuc}}$ is the representative volume (midpoint) in the section that covers the volume of nuclei and $d{v}_{\mathit{nuc}}$ is the interval of that section, respectively.

The approach described above will be used for the results in Section 6, except for Section 6.4.2 where some alternative nucleation models will be explored in order to investigate the sensitivity of the simulation to the choice of nucleation model.

4.3. Condensation

Condensation is considered as a collision process between silica monomers and particles. Let section i denote the range $(v_{i-1},v_i).$ In the Discretized PBE, three different condensation scenarios are considered that can change the number density, $n_i,$ within section i: (1)(1) $$\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u_i)}{\partial x_i}=0\hfill \end{array}$$(1) condensation of monomers on particles of section i, ${\omega }_{\mathit{cond},\mathrm{ }0\mathrm{ }\to \mathrm{ }i},$ (2)(2) $$\begin{array}{c}\frac{\partial (\rho u_i)}{\partial t}+\frac{\partial (\rho u_j{u}_i)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial {\sigma }_{ij}}{\partial x_j}+\rho g_i\hfill \end{array}$$(2) condensation of monomers on particles of section $i-1$ that will enter into section i, ${\omega }_{\mathit{cond},\mathrm{ }i-1\mathrm{ }\to \mathrm{ }i},$ and (3)(3) $$\begin{array}{c}\frac{\partial (\rho Y_{\alpha })}{\partial t}+\frac{\partial (\rho u_j{Y}_{\alpha })}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\rho D_{\alpha }\frac{\partial Y_{\alpha }}{\partial x_j}\right)+\rho \hspace{0.17em}{\dot{\omega }}_{\alpha }\hfill \end{array}$$(3) condensation of monomers on particles of section i that will go to section $i+1,$ ${\omega }_{\mathit{cond},\mathrm{ }i\mathrm{ }\to \mathrm{ }i+1}.$ The condensation rate for each of the scenarios mentioned above is evaluated according to (Sun, Rigopoulos, and Liu 2021; Rodrigues et al. 2018): (25) $$\begin{array}{c}{\omega }_{\mathit{cond},\mathrm{ }0\mathrm{ }\to \mathrm{ }i}=C_{{\text{SiO}}_2} N_A\underset{v_{i-1}}{\overset{v_i-v_o}{\int }}{\beta }_{0,i} n(w)dw\hfill \end{array}$$(25) (26) $$\begin{array}{c}{\omega }_{\mathit{cond},\mathrm{ }i-1\mathrm{ }\to \mathrm{ }i}=C_{{\text{SiO}}_2} N_A\underset{v_{i-1}-v_o}{\overset{v_{i-1}}{\int }}{\beta }_{0,i-1} n(w)dw\hfill \end{array}$$(26) (27) $$\begin{array}{c}{\omega }_{\mathit{cond},\mathrm{ }i\mathrm{ }\to \mathrm{ }i+1}=C_{{\text{SiO}}_2} N_A\underset{v_i-v_o}{\overset{v_i}{\int }}{\beta }_{0,i} n(w)dw\hfill \end{array}$$(27) where the collision frequency of a monomer and a particle in section i is evaluated as follows: (28) $$\begin{array}{c}{\beta }_{0,i}=1.28\sqrt{\frac{\pi k_{B}T}{2{\rho }_p}} {(d_o+d_{g,i})}^2 \sqrt{\frac{1}{v_o}+\frac{1}{v_{m,i}}}\hfill \end{array}$$(28) where $v_{m,i}$ is the volume at the middle of section i and $d_{g,i}$ is the collision diameter of the particles at section i evaluated according to Equation (5)(5) $$\begin{array}{c}{d}_g=d_p{\left(\frac{n_p}{k_f}\right)}^{1/D_f}\hfill \end{array}$$(5) . The amplification factor of 1.28 is due to Van der Waals interactions and its value was approximated based on the observations of Kennedy and Harris (1990). Following Sun, Rigopoulos, and Liu (2021), in terms of Discretized PBE, the condensation source terms in Equations (12)(12) $$\begin{array}{c}\frac{\partial n}{\partial t}+\frac{\partial (U_{j}n)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(D_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C(\mathbf{Y},v)\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)dw\hfill \end{array}$$(12) and (13)(13) $$\begin{array}{c}\frac{\partial n_p}{\partial t}+\frac{\partial (U_j{n}_p)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(n_{av}{D}_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C_p(\mathbf{Y},v)-\frac{3}{{\tau }_s}\left[\frac{n_p}{n}-{\left(\frac{n_p}{n}\right)}^{2/3}\right]\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)(n_{av}(w)+n_{av}(v-w))dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)n_{av}(v)dw\hfill \end{array}$$(13) in section i can be evaluated as: (29) $$\begin{array}{c}{C}_i(\mathbf{Y},v)=\left({\omega }_{\mathit{cond},\mathrm{ }0\mathrm{ }\to \mathrm{ }i}+{\omega }_{\mathit{cond},\mathrm{ }i-1\mathrm{ }\to \mathrm{ }i}-{\omega }_{\mathit{cond},\mathrm{ }i\mathrm{ }\to \mathrm{ }i+1}\right)\frac{v_o}{v_{m,i} d{v}_i}\hfill \end{array}$$(29) (30) $$\begin{array}{c}{C}_{p,i}(\mathbf{Y},v)=({\omega }_{\mathit{cond},\mathrm{ }0\mathrm{ }\to \mathrm{ }i} n_{av,i}+{\omega }_{\mathit{cond},\mathrm{ }i-1\mathrm{ }\to \mathrm{ }i} n_{av,i-1}-{\omega }_{\mathit{cond},\mathrm{ }i\mathrm{ }\to \mathrm{ }i+1} n_{av,i})\frac{v_o}{v_{m,i} d{v}_i}\hfill \end{array}$$(30) where $n_{av,i}$ is the average number of primary particles per aggregate for particles in section i and $d{v}_i$ is the particle volume range for that section.

4.4. Sintering

Sintering of two particles in contact is driven by the tendency of the system to reduce its surface free energy (Koch and Friedlander 1990). Thus, the sintering term in Equation (13)(13) $$\begin{array}{c}\frac{\partial n_p}{\partial t}+\frac{\partial (U_j{n}_p)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(n_{av}{D}_p(v,T)\frac{\partial n}{\partial x_j}\right)+B(\mathbf{Y},T)\delta (v-v_{\mathit{nuc}})+C_p(\mathbf{Y},v)-\frac{3}{{\tau }_s}\left[\frac{n_p}{n}-{\left(\frac{n_p}{n}\right)}^{2/3}\right]\hfill \\ +\frac{1}{2}\underset{0}{\overset{v}{\int }}\beta (w,v-w)n(w)n(v-w)(n_{av}(w)+n_{av}(v-w))dw-\underset{0}{\overset{\infty }{\int }}\beta (v,w)n(v)n(w)n_{av}(v)dw\hfill \end{array}$$(13) reduces the number of primary particles. Silica nanoparticles are amorphous materials with low viscosity, and viscous flow (Frenkel 1945) is considered the dominant mechanism when sintering takes place (Friedlander 2000). Many different expressions for the sintering characteristic time, ${\tau }_s,$ of silica nanoparticles appear in the literature without a consensus on which ${\tau }_s$ is more suitable for flame synthesis modeling (Park and Park 2015). Goudeli, Eggersdorfer, and Pratsinis (2015) discussed different sintering times of silica found in the literature and tested them in a perfectly stirred reactor (PSR).

In the present study, all four expressions shown in Table 3 were tested. The Tsantilis et al. expression provided the best results and was thus employed for all the results shown in Section 6, except for Section 6.4.1 where results with all four expressions are compared.

Table 3. Characteristic sintering times for silica nanoparticles tested in the current study.

Reference	Characteristic sintering time, ${\tau }_s$
Tsantilis, Briesen, and Pratsinis (2001)	$6.5·{10}^{-13}{d}_p\hspace{0.17em}\exp \hspace{0.17em}\left[\frac{8.3·{10}^4}{T}\left(1-\frac{d_{p,0}}{d_p}\right)\right]$
Ehrman, Friedlander, and Zachariah (1998)	$6.3·{10}^{-8}{d}_p\hspace{0.17em}\exp \hspace{0.17em}\left[\frac{6.1·{10}^4}{T}\right]$
Kirchhof et al. (2012)	$\frac{(1-1/e)r_s^2}{2^{2/3}(3.5·{10}^{-4})}\exp \hspace{0.17em}\left[\frac{4.4·{10}^4}{T}\right]$
Shekar et al. (2012)	$3.49·{10}^{-14}{d}_p\hspace{0.17em}\exp \hspace{0.17em}\left[\frac{11.44·{10}^4}{T}\left(1-\frac{5.76·{10}^{-9}}{d_p}\right)\right]$

The critical particle diameter in Tsantilis, Briesen, and Pratsinis (2001) expression is $d_{p,0}=1·{10}^{-9}$ m; $r_s$ in Kirchhof et al. (2012) expression denotes the radius of a completely coalesced spherical particle. The units are: m for $d_p,$ $d_{p,0},$ and $r_s,$ K for T and s for ${\tau }_s.$

4.5. Aggregation

Expressions for the aggregation kernel, $\beta (v,w),$ depend on the Knudsen number, $Kn=2{\lambda }_f/d_g,$ where ${\lambda }_f$ is the mean free path of the gas estimated by Willeke’s formula (Willeke 1976). The expressions for the free molecule particle size regime $(Kn<0.1)$ and the continuum regime $(Kn>10)$ are (Friedlander 2000): (31) $$\begin{array}{c}{\beta }_{fm}(v,w)=E_f\sqrt{\frac{\pi k_{B}T}{2{\rho }_s}}{\left(\frac{1}{v}+\frac{1}{w}\right)}^{1/2}{\left(d_{g,v}+d_{g,w}\right)}^2\hfill \end{array}$$(31) (32) $$\begin{array}{c}{\beta }_c(v,w)=\frac{2k_{B}T}{3\mu }\left(\frac{C_c(v)}{d_{g,v}}+\frac{C_c(w)}{d_{g,w}}\right)\left(d_{g,v}+d_{g,w}\right)\hfill \end{array}$$(32) where $C_c(v)=1+1.257Kn$ is the Cunningham slip correction factor. A harmonic mean of the two expressions (Pratsinis 1988) is employed for the transition regime $(0.1<Kn<10):$ (33) $$\begin{array}{c}{\beta }_{tr}(v,w)=\frac{{\beta }_{fm}(v,w){\beta }_c(v,w)}{{\beta }_{fm}(v,w)+{\beta }_c(v,w)}\hfill \end{array}$$(33)

4.6. Thermal radiation from particles

In heavily particle-laden flames, thermal losses due to the radiation from nanoparticles need to be accounted for in the model (Modest 2003). The radiative loss term due to the radiation from the nanoparticles is expressed as: (34) $$\begin{array}{c}{\dot{q}}_{\mathit{rad},p}=-\frac{4\sigma }{\rho } C_s f_v (T^5-T_{\mathit{amb}}^5)\hfill \end{array}$$(34) where $f_v$ is the silica volume fraction and $C_s$ is a proportionality constant. It can be shown that, in the optically thin limit, the heat loss due to emission of thermal radiation from a spherical object surrounded by a vacuum is proportional to the object’s volume (Modest 2003). Therefore, it is reasonable that ${\dot{q}}_{\mathit{rad},p}$ is proportional to $f_v$ in Equation (34)(34) $$\begin{array}{c}{\dot{q}}_{\mathit{rad},p}=-\frac{4\sigma }{\rho } C_s f_v (T^5-T_{\mathit{amb}}^5)\hfill \end{array}$$(34) . Furthermore, Lee (2011) studied the emission properties of silica aggregates at high temperatures and employed the below expression for the spectral radiant intensity $E_{\lambda }$ from the particles in a flame: (35) $$\begin{array}{c}{E}_{\lambda }=\frac{A·2\pi h{c}^2ϵ(\lambda ,T)f_v}{{\lambda }^5\left[\exp \hspace{0.17em}\left(\frac{hc}{\lambda k_{b}T}\right)-1\right]}\hfill \end{array}$$(35) where h is Planck’s constant, c the speed of light in the medium, $\lambda$ the wavelength, $ϵ(\lambda ,T)$ the emissivity of silica and A is an arbitrary constant. Lee (2011) observed that $ϵ$ is a linear increasing function of temperature, while $ϵ(\lambda ,T)$ is expected to be a weak function of $\lambda$ at the wavelengths of interest; therefore, the emissivity of silica is assumed to follow the relationship $ϵ\propto (C_s/A)T.$ Based on these considerations, integration of Equation (35)(35) $$\begin{array}{c}{E}_{\lambda }=\frac{A·2\pi h{c}^2ϵ(\lambda ,T)f_v}{{\lambda }^5\left[\exp \hspace{0.17em}\left(\frac{hc}{\lambda k_{b}T}\right)-1\right]}\hfill \end{array}$$(35) over $\lambda$ and all possible directions results in an expression for the radiative heat flux similar to Equation (34)(34) $$\begin{array}{c}{\dot{q}}_{\mathit{rad},p}=-\frac{4\sigma }{\rho } C_s f_v (T^5-T_{\mathit{amb}}^5)\hfill \end{array}$$(34) , giving support to the radiation model adopted here. In the absence of any information for the value of A, the value $C_s=1100$ m${}^{-1}$ K${}^{-1}$ was selected in the present study, which is lower than the value of 1307 m${}^{-1}$ K${}^{-1}$ used in sooting flame simulations and gave good agreement with the experimental results, as will be shown in Section 6.1.

5. Flame configuration and computational setup

5.1. Flame configuration

The flame under investigation is the ‘S-2’³ laminar case of Camenzind et al. (2008) and, to the best of our knowledge, the present work is the first simulation of this experiment. The synthesis of silica nanoparticles takes place by oxidation of HMDSO in oxygen/methane diffusion flames. This flame was selected because it is one of the few laminar flames for which in-situ experimental data for the particles are available in the literature (Table 1); Camenzind et al. (2008) provided a plethora of particle data for several heights above the burner (HAB) obtained via in-situ SAXS measurements.

The experimental setup is shown in Figure 1. The co-flow diffusion flame reactor used in the experiment 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 2003a). Mass flow controllers were employed to maintain constant levels of all gas flow rates, which are given below at standard temperature and pressure (STP) conditions. In particular, 0.3 l/min of N₂ passed through a reservoir filled with liquid HMDSO to deliver HMDSO vapor. The HMDSO consumption rate was 6.5 g/h, which corresponds to a SiO₂ production rate of 4.8 g/h. Subsequently, the HMDSO-laden N₂ gas flow was mixed with 0.5 l/min of methane CH₄ and delivered in the reactor through the center tube. The first annulus was fed with 0.5 l/min of N₂ to prevent particle deposition on the reactor tips. The second annulus was fed with 2 l/min of O₂. The temperature of the reactor and 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 reactor, was $Re\approx 1414.$

5.2. Computational setup

Two-dimensional axisymmetric simulations were performed, with an axisymmetric boundary condition being applied along the centerline and a symmetry boundary condition being applied at the opposite side. A zero-gradient boundary condition was applied at the outlet boundary in the streamwise direction. In the plots, the r coordinate is defined to denote the cross-stream direction. The stream bulk velocities $u_{\mathit{inlet}}$ of the jets coming from the inner tube, the first and second annulus were $U=6.715$ m/s, $U_1=2.057$ m/s, and $U_2=8.57$ m/s, respectively. The inlet velocities were estimated based on the given mass flow rates reported earlier and by assuming that the temperature of the streams was 75${}^{°}$ C. At the inflow, laminar flow velocity profiles were employed for the inner jet as well as through the annular sections for the other two jet streams. The inlet velocity profiles are shown in Figure S2 of the supplementary information. In the absence of any information for the entrainment velocity, a uniform velocity $U_{\mathit{air}}=0.4$ m/s was used for the air co-flow stream at the inlet and the ambient temperature was set to $T_{\mathit{amb}}=20$${}^{°}$C. Test cases showed that $U_{\mathit{air}}$ has a minor effect on the results (Figure S3 of the supplementary information). The mole fractions of the species for the inner stream were set to $X_{{\text{CH}}_4}=61.338\%,$ $X_{{\mathrm{N}}_2}=36.803\%$ and $X_{\mathit{HMDSO}}=1.859\%.$ The thickness of the burner tubes was 0.3 mm, which was taken into account in the simulation by applying a no-slip boundary condition at these radii. The inlet boundary conditions employed in this study are summarized in Table 4.

Figure 2. Contour plots of (a) velocity and (b) temperature fields.

Figure 3. Contour plots of (a) silica volume fraction and (b) primary particle diameter fields.

Table 4. Inlet boundary conditions, where $u_{\mathit{inlet}}$ denotes the stream bulk inlet velocity in the streamwise direction.

Stream	Profile	$u_{\mathit{inlet}}$ [m/s]	T [${}^{°}$C]	Mole fractions [%]
Inner	Laminar	$U=6.715$	75	$X_{C{H}_4}=61.338;$ $X_{N_2}=36.803;$ $X_{\mathit{HMDSO}}=1.859$
$1^{st}$ annulus	Laminar	$U_1=2.057$	75	$X_{N_2}=100$
$2^{nd}$ annulus	Laminar	$U_2=8.57$	75	$X_{O_2}=100$
Co-flow air	Uniform	$U_{\mathit{air}}=0.4$	20	$X_{N_2}=76.83;$ $X_{O_2}=23.17$

The inlet velocity profile at the inflow is presented in Figure S2 of the supplementary information.

The size of the computational domain is $60D\times 15D,$ where $D=1.8$ mm is the inner tube diameter. A computational grid with $264\times 162$ cells with spacing increasing exponentially by factors 1.0132 and 1.03 in the axial and radial directions, respectively, was employed. Images of the computational domain and grid are shown in Figure S4 of the supplementary information. Grid refinement was employed at the height of the second annulus to ensure that at least ten cells were used to capture the inlet velocity profile. Tests with a finer grid comprising about 4.5 times more cells showed that grid independence had been achieved, and the results are shown in Figure S5 of the supplementary information. The time step was $\Delta t={10}^{-7}$ s, resulting in a maximum Courant-Friedrichs-Levy (CFL) number of 0.03. Each simulation was run until the solution reached a steady state independent of the initial conditions. Following this, the simulations were run for more than 3 million time steps, corresponding roughly to 0.4 s of the simulated process.

Figure 4. Comparison of numerical predictions (two-PBE model) with the experimental FTIR data of Camenzind et al. (2008) for flame temperature along the centerline and effect of radial averaging. An average experimental error of $\pm$60 K given by Kammler et al. (2002) was used in the figure for the FTIR data.

Figure 5. Comparison of numerical predictions (two-PBE model) with experimental data for flame temperature and effect of emission of thermal radiation from gas species and from particles. Results of the radial maximum temperature for each case are compared with the experimental FTIR data of Camenzind et al. (2008). An average experimental error of $\pm$60 K given by Kammler et al. (2002) was used in the figure for the FTIR data.

Regarding the Discretization of the PBE, the particle volume domain covered particle sizes from 0.44 nm to 1 $\mu$ m in diameter and was discretized with an exponential grid comprising 60 intervals. Convergence studies with a perfectly stirred reactor showed that this grid sufficed to describe the PSD evolution accurately.

Finally, to ensure that conservation of mass was achieved in the simulation, the silica mass flow rate at the exit of the computational domain was calculated as ${\int }_{A_{\mathit{out}}} {\rho }_{\mathit{SiO}2} f_v \mathbf{u}·\mathbf{dA}=4.755$ g/h. The relative difference with the reported (Camenzind et al. 2008) silica production rate (4.8 g/h) is less than 1%.

6. Results and discussion

In the present section, the results are presented and discussed. Particular emphasis is given to the validation of the modeling strategy adopted in this study by comparing numerical predictions of particle field quantities with reference in-situ experimental data. Numerical simulations performed with the three population balance models described in Section 2 are referred to as ‘2PBE’, ‘1PBE’, and ‘Mono’, while ‘Exp’ refers to experimental results. An analysis of the calculated timescales and rates is presented to draw insight into the model predictions. Results obtained from alternative sintering and nucleation models are also presented in order to demonstrate the sensitivity of the predictions to modeling choices. Finally, the computational performance of the PBE models is analyzed.

6.1. Flow field and temperature

Contour plots of the velocity and temperature fields are shown in Figure 2. The numerical results from the three models tested in this study were in good agreement with each other in terms of temperature and velocity fields; therefore, contour plots only for the two-PBE model are shown. The simulation was allowed to reach a steady state before the results were extracted. The volume fraction throughout the domain is shown in Figure 3a. Particles are formed within the flame reaction zone, which is located approximately 1.4 mm away from the flame centerline. It can be seen that particles stay within this thin layer at all downstream locations. This is reasonable considering the very low diffusivity of particles (estimated from the Stokes-Einstein relation (Friedlander 2000)) and the laminar flow in the experiment.

Furthermore, a contour plot of the primary particle diameter is shown in Figure 3b. Primary particles have diameters less than 15 nm at low HABs, while they grow in size rapidly at $3.5<$ HAB $<5$ cm at the flame front and eventually reach constant values downstream. It is observed from Figure 3a,b that the maximum primary particle sizes appear within the thin layer where the maximum silica volume fraction is located. As will be explained later, in the vicinity above the burner, the particle formation rate is comparable to particle aggregation rates, resulting in small average primary particle sizes. After the flame front at $3.5<$ HAB $<5$ cm, conversion of the precursor is almost complete, the effect of particle nucleation is less significant and sintering becomes the dominant mechanism, leading to rapid primary-particle growth. At HAB $>5$ cm, sintering ceases, and aggregation becomes the dominant mechanism. This explains why the primary particle size remains almost constant at these locations.

A comparison of numerical results for temperature with the experimental data of Camenzind et al. (2008) is shown in Figure 4. As before, only results with the two-PBE model are shown in Figures 4 and 5 since all of the three PBE models resulted in similar temperature distributions. Accurate prediction of temperature is the cornerstone of flame synthesis modeling, as it is driving many of the subsequent processes. In the present study, numerical results for flame temperature are in good agreement with the experimental data, largely due to the direct coupling of detailed chemistry with flow, the absence of modeling uncertainties involved in turbulent flame simulations (turbulence modeling and turbulence-chemistry interaction), and the inclusion of the model for the thermal radiation from particles.

The experimental temperature data were obtained using Fourier transform infrared (FTIR) spectroscopy (Griffiths and De Haseth 2007; Morrison et al., 1997) which is a line-of-sight technique over the total flame width. It has been documented that the FTIR data represent a radial-average temperature of the diffusion flame in the region close to the burner tip, where reactants are not mixed (Camenzind et al. 2008) (in fact, in this region, the temperature on the flame centerline is well below the reported values, which can also be observed in Figure 2b), while further downstream the radial average does not differ significantly from the centerline value. For this reason, numerical results for the radial maximum value of temperature at each height above the burner are plotted (solid black line) in Figure 4, instead of the centerline values. Radial maximum temperatures are more representative of the experimental data. Emitted radiation (relevant to FTIR measurements) scales with the fourth power of temperature; thus, it can be expected that regions with the highest temperature would dominate in the FTIR data analysis. Nevertheless, radial averaging of values of cells in the flame reaction zone was performed, and results are also presented in Figure 4 (solid magenta line). The flame reaction zone was approximated with the cells whose value of CO₂ concentration is greater than 80% of the radial maximum value at each height above the burner. However, both lines in Figure 4 lie within an average experimental error of $\pm$60 K (Kammler et al. 2002) for most HABs, indicating that the flame temperature was predicted accurately by the simulation.

In order to quantify the effect of thermal radiation from gas species and from particles, two alternative simulations were performed and are compared with the one including radiation in Figure 5. In one case (red line), all radiation thermal losses were neglected, while in another (green line) radiation from gas species only was considered. When only radiation from gas species is considered, then the CFD model considerably overestimates the experimental data at downstream locations. A better agreement with the experimental data is achieved when the emission of radiation from particles is included. These results suggest that thermal losses due to radiation from particles have a non-trivial effect on the evolution of the temperature field.

6.2. Comparison with experimental SAXS data

Comparisons of numerical predictions with the experimental SAXS data of Camenzind et al. (2008) will now be shown. In SAXS measurements, a monochromatic X-ray beam penetrates the flame and is collected on a SAXS detector located a few meters away from the flame. However, particle quantities cannot be uniform in the radial direction in the vicinity of the burner tip where reactants are not mixed. As discussed earlier, the flame reaction zone, where particles are formed, is located slightly off the centerline. The experimental particle data refer to the thin flame zone where the scattered intensity (relevant to SAXS measurements) should be maximum. This was also discussed in Sztucki, Narayanan, and Beaucage (2007) where it was found that the maximum scattered intensity was achieved when the X-ray beam penetrated the flame in the soot growth zone located slightly off the centerline of the flame. For this reason, we present data along the pathline of maximum silica volume faction, accounting also for the radius of the beam 0.1 mm. It should also be mentioned that the smallest detectable particle size by the SAXS measurements was 1 nm, and this was taken into account in the post-processing of the simulation results by not accounting for particles smaller than that when averaging the PSD.

The average primary particle diameter, $d_{p,av},$ the average aggregate radius of gyration, $R_{g,av},$ the number concentration of primary particles, $N_p,$ and the average number of primary particles per aggregate, $n_{p,av},$ along the pathline of maximum silica volume fraction are shown in Figure 6a–d, respectively. The plots show the results from the three PBE models and the measurements. Focusing initially on the results obtained with the more comprehensive two-PBE model, we can see that good agreement with experiments is found for primary particle diameters. In particular, the model predicts relatively well $d_{p,av}$ at low HABs while overestimating by about 15% the experimental data at downstream locations (HAB $>5$ cm). This indicates that aggregation and sintering are represented well by the model. Furthermore, fair agreement can be observed for $R_g$ and $n_{p,av}.$ In particular, the model overestimates $R_g$ (by about 43% at HAB $=3$ cm), most likely due to uncertainties involved in the precursor decomposition kinetics/nucleation model and the questionable validity of Equation (5)(5) $$\begin{array}{c}{d}_g=d_p{\left(\frac{n_p}{k_f}\right)}^{1/D_f}\hfill \end{array}$$(5) at the first stages of particle synthesis. Furthermore, $n_{p,av}$ is overestimated, and this is consistent with the overprediction of $N_p$ discussed below.

Figure 6. Comparison of the experimental SAXS data of Camenzind et al. (2008) with numerical predictions for (a) primary particle diameter, (b) radius of gyration of aggregates, (c) primary particle number concentration and (d) average number of primary particles per aggregate along the pathline of maximum silica volume fraction. Results were obtained with the two-PBE, the extended one-PBE, and the monodisperse model (described in Section 2).

Numerical predictions for $N_p$ are in qualitative agreement with the SAXS data. However, the particle formation rate is overpredicted, leading to an overestimation of $N_p$ by approximately 1.5 orders of magnitude (Figure 6c). This discrepancy is most likely related to uncertainties involved in the SAXS measurements, precursor decomposition kinetics and the nucleation model. Regarding the precursor decomposition/oxidation reaction mechanism in particular, data obtained with low-pressure flames and low HMDSO concentrations were employed for the development of the two-step chemical reaction mechanism of Feroughi et al. (2017), while the flame examined here was at atmospheric pressure and the precursor concentration was relatively high. Nevertheless, it should be mentioned that the good agreement with primary particle diameter is more important, as this quantity is of more practical interest.

Aggregates are formed at low HABs where $n_{p,av}$ takes values close to 200, while compact - almost spherical - particles are formed at downstream locations where $n_{p,av}$ declines quickly and reaches values close to 3 (Figure 6d). The reduction of $n_{p,av}$ is also observed in the SAXS data at HAB $=4$ cm. Experimental data for downstream locations are unavailable, but the model suggests that the residence time of particles at high temperatures (i.e., at $3.5<$ HAB $<5.5$ cm) is large enough for almost complete coalescence of aggregates, leading to compact structures at HAB $>6$ cm. This can also be seen in Figures 6a,b by comparing the values of $d_{p,av}$ and $R_{g,av}$ in the two-PBE model predictions at downstream locations; twice $R_{g,av}$ is only slightly greater than $d_{p,av}.$

Results obtained with the extended one-PBE and the monodisperse model is also presented in Figure 6. Overall, good agreement was found between the two-PBE and the extended one-PBE models. Specifically, the extended one-PBE model somewhat overestimated primary particle sizes at downstream locations (by about 35% relative to the SAXS data and by about 10% relative to the two-PBE results at HAB $=10$ cm), while excellent agreement was found between the two models for the rest of the particle quantities. On the other hand, the monodisperse model underestimates $d_{p,av}$ (Figure 6a) (by about 26% relative to the SAXS data and by about 39% relative to the two-PBE results at HAB $=10$ cm). Furthermore, the monodisperse model underestimates $R_{g,av}$ (Figure 6b) and $n_{p,av}$ (Figure 6c), at the early stages of particle synthesis, i.e., at HAB $<5$ cm, while the two-PBE and the extended one-PBE models describe the evolution of particle morphology in a more accurate manner (both quantitatively and qualitatively).

Measurements are also available for the geometric standard deviation, ${\sigma }_g,$ of the primary particle size distribution. This quantity can be predicted only by the two-PBE model. However, as explained earlier, aggregates at downstream locations are quite compact and comprise only a few primary particles, most likely connected by sinter bonds; thus, they would be difficult to distinguish by SAXS measurements. It is therefore worth comparing with predictions of ${\sigma }_g$ for both primary particles and aggregates. Figure 7 shows the measurements and the predictions by the two-PBE (for both primary particles and aggregates) and the one-PBE (for aggregates only) models; results with the monodisperse model are not included since ${\sigma }_g$ in that case is unity by definition. It can be seen that better agreement is found when comparing with the simulation results for aggregates than for primary particles. In particular, the figure shows good agreement between the numerical predictions of ${\sigma }_g$ for aggregates with the SAXS data (approximately 1.4% relative difference at HAB $=10$ cm) at downstream locations where aggregation is the dominant mechanism and the PSD attains the self-preserving distribution (the ${\sigma }_g$ for aggregates approaches an asymptotic value of approximately 1.38). The two-PBE and one-PBE models are in good agreement with each other except at low HABs, where both overpredict the SAXS data. At that stage, nucleation, condensation and aggregation take place simultaneously and the PSD is far from the self-preserving state. The disagreement can be partly explained by the fact that ${\sigma }_g$ was estimated from SAXS measurements based on the measured polydispersity index and assuming a lognormal particle size distribution of spherical primary particles (Kammler et al. 2005), while at these HABs the PSD was far from lognormal.

Figure 7. Comparison of the experimental SAXS data of Camenzind et al. (2008) with numerical predictions for the geometric standard deviation of aggregate (agg) and primary particle (pp) size distributions along the pathline of maximum silica volume fraction.

6.3. Timescale analysis

The coupling of CFD with population balance allows the calculation of the timescales of different processes at all points in the domain, which can provide further insight on the results presented so far and furthermore help identify the dominant mechanisms. Timescales are defined for nucleation, aggregation and sintering in a consistent way and summarized in Table 5; in all cases, the total number of particles (primary or aggregates, depending on the process) is divided by the rate of the process (or its absolute value when needed). We also define a flow timescale as the residence time, ${\tau }_{\mathit{res}},$ of particles along the pathline of maximum silica volume fraction.

Table 5. Definition of timescales for processes encountered in silica flame synthesis.

Process	Nucleation	Aggregation	Sintering	Flow
Timescale, $\tau$	$\frac{{\int }_0^{\infty }n(v)dv}{J}$	$\frac{{\int }_0^{\infty }n(v)dv}{\left|-\frac{1}{2}{\int }_0^{\infty }{\int }_0^{\infty }n(w)n(v)\mathit{dwdv}\right|}$	$\frac{{\int }_0^{\infty }{n}_p(v)dv-{\int }_0^{\infty }n(v)dv}{\left|-\frac{3}{{\tau }_s}\left(\frac{V}{v_{p,av}}-\frac{M_{23}}{v_{p,av}{}^{2/3}}\right)\right|}$	${\int }_0^x\frac{1}{u_x(x^{\prime })}d{x}^{\prime }$

Results for the nucleation, aggregation and sintering rates along the pathline of maximum silica volume fraction are presented in Figure 8a, while results for the timescales are depicted in Figure 8b. The value of $N_p$ is quite constant at HAB $<3.5$ cm, while it starts to drop at downstream locations (Figure 6c). At low HABs, the flame front is positioned slightly off the centerline and particles are constantly formed within the flame reaction zone. Nucleation and aggregation rates reach comparable values (Figure 8a) within this zone, resulting in constant primary particle number concentrations in this region (Figure 6c) and relatively small average primary particle sizes (Figure 6a). On the other hand, after the flame front on the centerline (HAB $\approx 3.5$ cm), precursor conversion is almost complete; the nucleation rate decreases by several orders of magnitude and sintering of particles becomes the dominant mechanism. This is demonstrated by Figure 8a, where it can be seen that the sintering rate is higher than the nucleation and aggregation rates at $4<$ HAB $<5.6$ cm. In terms of timescales, the sintering timescale is smaller than (or comparable to) the nucleation and aggregation timescales at $4<$ HAB $<6$ cm (Figure 8b), indicating that the effect of sintering is more profound in this region. This explains the reduction of $N_p$ and the apparent growth of the average primary particle sizes (Figure 6a) at these HABs. At HAB $>5$ cm, the temperature decreases, attenuating the effect of sintering and resulting in constant $d_{p,av}$ values. Aggregation becomes the dominant mechanism at HAB $\approx 6$ cm (see the intersection of the lines in Figure 8b). The sintering timescale becomes greater than the residence time of particles in the flow at HAB $>7$ cm, indicating that the effect of sintering is negligible in that region.

Figure 8. Rates (absolute values) of kinetic processes (a) and timescales (b) along the pathline of maximum silica volume fraction. Results were obtained with the two-PBE model. Note that, for the sintering process in (a), the rate refers to the number concentration of primary particles, $N_p.$

6.4. Investigation of alternative model choices

6.4.1. Models for the characteristic sintering time

As explained in Section 4.4, the results shown so far were obtained with the Tsantilis, Briesen, and Pratsinis (2001) expression for the characteristic sintering time, ${\tau }_s.$ Simulations were also performed with the other models shown in Table 3, and the results for $d_{p,av}$ and $R_{g,av}$ are shown in Figure 9. The model of Ehrman, Friedlander, and Zachariah (1998) predicts greater sintering times (smaller coalescence rates) compared to that of Tsantilis, Briesen, and Pratsinis (2001), leading to an underestimation of $d_{p,av}$ for all HABs. The experimentally derived expression of Kirchhof et al. (2012) fails to capture the coalescence of the nucleus-size nanoparticles (it must be noted that particle diameters larger than 10 nm were considered by Kirchhof et al., 2012) and thus does not predict well the particle morphology. Shekar et al. (2012) proposed an expression where the sintering parameters were acquired by fitting their model to the experimental data of Seto et al. (1997). The resulting sintering time leads to an overestimation of $d_{p,av}$ at low HABs and an underestimation of $d_{p,av}$ at large HABs. The use of the expression of Tsantilis, Briesen, and Pratsinis (2001) slightly overpredicts $d_{p,av}$ at the downstream locations. The model of Tsantilis, Briesen, and Pratsinis (2001) accounts for the particle-size dependence of the melting point of silica and the parameter $d_{p,0}=1·{10}^{-9}$ m (Table 3) was proposed as a first approximation. Overall, results obtained with the ${\tau }_s$ of Tsantilis, Briesen, and Pratsinis (2001) are in better agreement with the experimental data compared to the other ${\tau }_s$ cases, at least for the flame conditions examined in the present study. Regarding $R_{g,av},$ the sintering models predict similar trends and slightly overestimate the SAXS data, with the exception of the Kirchhof et al. (2012) model, which yields a greater overprediction.

Figure 9. Effect of sintering time on (a) primary particle diameter, (b) aggregate radius of gyration. Results were obtained with the two-PBE model. Comparison is also made with the experimental data of Camenzind et al. (2008).

6.4.2. Nucleation models

The reliable prediction of particle nucleation rates remains an unresolved problem in the field of flame synthesis of nanomaterials (Rosner 2005). In the present section, we investigate some alternatives to the approach described in Section 4.2 and used to obtain the results presented so far.

It is customary in flame-synthesis studies to assume that fast precursor oxidation leads to silica monomers (molecules) that serve as the first primary particles (Ulrich 1971). This approach will be referred to here as the instantaneous nucleation assumption. The reasoning behind this approach is that monomer vapor is usually in a supersaturated state, the nucleation barriers can be neglected (Xiong and Pratsinis 1991) and, even if the critical clusters (nuclei) comprise few molecules, coagulation quickly extinguishes initial differences (Gröhn et al. 2011; Tsantilis, Kammler, and Pratsinis 2002). Theoretical tools to describe nucleation rate have also been developed, with the most celebrated one being the classical nucleation theory, a discussion of which can be found, for example, in Ford (2004). Here, we compare instantaneous nucleation, the classical nucleation theory as formulated by Girshick and Chiu (1989) and the dimers-based approach that was discussed in Section 4 and used for the results presented in Section 6. Results obtained with the three nucleation models for $d_{p,av}$ and $R_{g,av}$ along the pathline of maximum silica volume fraction are presented in Figure 10a,b, respectively. The classical theory overestimates $d_{p,av}$ at low HABs, while the instantaneous model overestimates $R_{g,av}$ at low HABs. The dimers-based approach seems to describe better the evolution of the particle morphology; however, all models predict similar particle sizes at downstream locations.

Figure 10. Effect of nucleation model on (a) primary particle diameter and (b) aggregate radius of gyration. Three models are compared: instantaneous nucleation (InstNuc), classical nucleation theory (CNT) and the dimers-based approach that was used for the rest of the results. Results were obtained with the two-PBE model. Comparison is also made with the experimental data of Camenzind et al. (2008).

6.5. Analysis of computational performance

The total CPU time per timestep for simulating all processes (including flow field, scalar convection/diffusion, gas phase reaction, and PBE solution) for each of the three models is presented in Table 6 and it can be seen that the simulation with the one-PBE model requires about a third of the CPU time spent on the two-PBE model. Remarkably, it is even comparable with the time required by the monodisperse model (less than twice). This difference is much greater than suggested by the number of equations in each model (116 for the two-PBE compared with 61 equations for the extended one-PBE) and is explained by the fact that the Discretization of the particle volume domain in the one-PBE model allows for the tabulation of the aggregation kernel. On the other hand, tabulation of the aggregation kernel for the two-PBE was not employed because it would be very memory-intensive; in particular, tabulation of the kernel in the one-PBE requires storing a three-dimensional matrix (for two volumes and $d_p$) since, at one time instant, $d_p$ is the same for all particle volume pairs, while the two-PBE requires a four-dimensional matrix since every section has a different $d_p.$ This discussion demonstrates that coupling CFD with population balance is computationally feasible, while the extended one-PBE model is a major step forward in terms of balancing accuracy and computational speed. It should be emphasized that, while this reduction in CPU time is not essential for the simulation of the present laminar flame as the total CPU time is manageable with all approaches, the advantages offered by the extended one-PBE model will be crucial when it comes to the coupling of detailed population balance modeling with LES for the simulation of turbulent flame synthesis.

Table 6. Comparison of PBE models for describing aerosol dynamics and particle morphology.

Model	two-PBE	Extended one-PBE	Mono
No. of equations	116	61	3
PSD of aggregates	polydisperse	polydisperse	monodisperse
PSD of primary particles	polydisperse	monodisperse	monodisperse
CPU time (s) / timestep	1.14	0.32	0.187

The total CPU time per timestep includes all processes (flow field, scalar convection/diffusion, gas phase reaction, and PBE solution). Tabulation of the aggregation kernel was employed for the one-PBE model.

7. Conclusions

In the present work, the synthesis of silica nanoparticles in a laminar diffusion flame was studied via population balance modeling coupled with CFD. Three population balance models were investigated: a monodisperse model, a detailed two-PBE model and an extended one-PBE model. The last of these constitutes a new development and is an augmentation of the basic PBE with a finite-rate sintering expression, while an additional transport equation is solved for the total number of primary particles. The extended one-PBE model can predict the evolution of particle morphology while reducing substantially the computational cost as compared with the two-PBE model.

The models were coupled with an in-house CFD code in order to conduct two-dimensional simulations of silica synthesis via oxidation of HMDSO in a laminar diffusion jet flame. The methodology was validated against the experimental in-situ SAXS data of Camenzind et al. (2008). Results with all models were in excellent agreement with the experimental data for the temperature. Thermal losses due to radiation from particles were found to have a non-trivial effect on the evolution of the temperature field. The predictions by the two-PBE and extended one-PBE model for primary particle sizes were in good agreement with the SAXS data indicating that aggregation and sintering were described relatively well, while the models slightly overestimated the aggregate sizes and the number of primary particles per aggregate. By contrast, the monodisperse model did not capture particle morphology at the early stages of the synthesis. The particle formation rate was overpredicted by all models, leading to an overestimation of the number concentration of primary particles. This discrepancy is most likely due to uncertainties in the experimental data, precursor decomposition kinetics and the nucleation model. The coupling of CFD and PBE also allowed the calculation of the rates and timescales of all processes along the domain. This analysis shed light on the mechanisms underlying the evolution of particle properties and the dominant mechanisms at different spatial locations. An investigation of alternative nucleation and sintering models was also conducted and the analysis showed that, overall, results obtained with the dimers-based nucleation model and the sintering characteristic time of Tsantilis, Briesen, and Pratsinis (2001) were in better agreement with the experimental data compared to those obtained with the other models.

The one-PBE model yielded results close to those of the two-PBE model at a substantially reduced computational cost, indicating that it provides a good balance of physical detail and speed. Most importantly, the time taken for the CFD-PBE simulation with the one-PBE model was less than twice the time needed for the CFD-monodisperse model, while the former provided a much better description of the evolution of particle morphology. The superior speed of the one-PBE model was largely owing to the tabulation of the aggregation kernel.

Overall, the results showed that the coupling of CFD and detailed PBE is capable of describing the evolution of particle properties during silica flame synthesis. The computational advantage offered by the extended one-PBE model paves the way for its coupling with LES and the comprehensive simulation of aerosol synthesis in turbulent flames, which will be the objective of future work.

Acknowledgments

The authors gratefully acknowledge Prof. Sotiris E. Pratsinis for his help on the interpretation of the experimental results of Camenzind et al. (2008) and for insightful discussions on silica synthesis, as well as Dr Eirini Goudeli for her help on the validation of the monodisperse model.

Disclosure statement

The authors report no conflict of interest.

Additional information

Funding

The authors are grateful to the Leverhulme Trust for financial support (grant reference RPG-2018-101), as well as to the UK Materials and Molecular Modelling Hub for computational resources on YOUNG, which is partially funded by EPSRC (EP/P020194/1 and EP/T022213/1) and also supported by Engineering and Physical Sciences Research Council.

Notes

1 Chemical species are denoted with the subscript $\alpha$ and no summation is implied by repeated Greek subscripts.

2 Unlike J in Equations (6)–(8), B is a source of particle number density rather than particle number.

3 Camenzind et al. (2008) used the notation “S-y” to refer to a series of SiO₂-producing flames with different O₂ flow rates, y [l/min].