Soot growth by monodisperse particle dynamics coupled with computational fluid dynamics

A multiscale modeling framework, integrating MD (molecular dynamics)-derived soot nucleation and surface growth rates into a coupled computational fluid dynamics (CFD)–monodisperse particle dynamics (PD) model, is implemented and benchmarked for a premixed ethylene burner-stabilized stagnation (BSS) flame. The proposed coupled model is validated by comparing the soot number density, volume fraction, and particle size with measurements across the BSS flame, as well as with the results obtained by CFD-PD using a semi-empirical nucleation rate by Moss–Brookes and the hydrogen abstraction-carbon addition (HACA) surface growth rate. Incorporation of the MD-derived nucleation rate is in excellent agreement with both experimental data and a detailed sectional model from the literature, especially in the post-flame region. The proposed MD-informed CFD-PD model is computationally efficient compared to detailed population balance equation models as it does not rely on reaction kinetic modeling. To the best of our knowledge, this work is the first systematic assessment of the potential of MD-derived nucleation and surface growth rate constants for soot modeling. In addition, the proposed CFD-coupled monodisperse PD framework can serve as a predictive tool for soot modeling and design-oriented simulations of practical combustion and aerosol systems.

1 Introduction

Soot formation in combustors and flames is inherently multiscale, involving coupled chemical and physical processes (Michelsen, 2017) that span several orders of magnitude in length and time. At the molecular zone, soot formation begins with reactions of small radicals and hydrocarbon fragments, only a few Angstroms in size, that form polycyclic aromatic hydrocarbons (PAHs), the main precursors of soot (Frenklach, 2002) as shown by atomic force microscopy (Schulz et al., 2019). These PAHs evolve into sub-2 nm soot nuclei through physical or chemical nucleation, which grow by surface reactions into primary particles (Mitchell and Frenklach, 2003). The primary particles coagulate by aggregation and/or agglomeration to form large (>100 nm) aggregates within microseconds to milliseconds (Michelsen, 2017). Reactive flow in combustors and flames extends over dimensions on the order of meters and particle residence times on the order of seconds, which strongly influence soot evolution (Trivanovic and Pratsinis, 2024; Kelesidis et al., 2017).

Hierarchical modeling provides a pathway to address the multiscale nature of soot formation and growth by linking information across different length and timescales (Goudeli, 2019). Molecular dynamics (MD) simulations can capture the onset and early stages of soot nucleation (Mao et al., 2017), as well as surface growth (Chen and Luo, 2020) and coalescence (Sharma et al., 2021), but remain limited to a few nanometers and hundreds of nanoseconds. Traditionally, soot nucleation has been modeled using physical PAH–PAH interactions, such as with the POLIMI soot nucleation sub-model (Saggese et al., 2015), whereas more recent approaches adopt hybrid (physical–chemical) reactive PAH dimerization models (Kholghy et al., 2018). More recently, reactive molecular dynamics (MD) has been employed to determine the critical nucleus size and nucleation rate of carbonaceous nanoparticles (Fakharnezhad et al., 2025). Particle-scale models, such as event-driven simulations (Patterson and Kraft, 2007; Goudeli et al., 2015; Kelesidis et al., 2017), can resolve particle interactions, while continuum models, such as population balance equations (PBEs) and computational fluid dynamics (CFD) can track particle growth and transport (Saggese et al., 2016), respectively, reaching length and timescales of micrometers and seconds. For example, event-driven mesoscale simulations of soot nanoparticles have been used to determine the coagulation rate in the transition regime, yielding an enhancement of 82 % over the harmonic mean interpolation between the theoretical kernels of the free molecular and continuum regime (Kelesidis and Goudeli, 2021). Such corrections have been implemented into monodisperse particle dynamics simulations (Fakharnezhad et al., 2026) to account for the primary particle and aggregate polydispersity. A similar hierarchical modeling framework has also been applied for TiO₂ nanoparticle synthesis, where fractal dimension and coagulation rate equations derived from Langevin dynamics (Chen et al., 2023) were incorporated into a CFD-PBE model (Chen et al., 2023), along with sintering rates obtained by MD simulations (Buesser et al., 2011), linking molecular-, particle- and continuum-scale simulations.

The evolution of soot particle population can be described by the general dynamic equation, a partial integro-differential equation, requiring the implementation of numerical methods. One highly accurate numerical approach relies on the sectional method, where the particle population is divided into size bins and transport equations are solved for each class (Saggese et al., 2015), but its computational cost increases with the number of sections, making it impractical for multidimensional CFD simulations (Saggese et al., 2016; Kholghy and Kelesidis, 2021). A computationally less demanding alternative is the method of moments (MOM) (Hulburt and Katz, 1964), which tracks only a few statistical moments of the distribution (e.g., number, mass, surface area). However, MOM relies on closure assumptions (Kholghy and Kelesidis, 2021), limiting accuracy when soot distributions deviate from the assumed shape.

The particle dynamics equation can be greatly simplified when assuming monodisperse particles or aggregates at a given location (Kruis et al., 1993). So, the soot population can be represented by a single average size, requiring tracking of the evolution of only the number density, surface area, and volume (Kholghy and Kelesidis, 2021). Even though PDs neglect the spread of soot size distribution, this simplification typically results in deviations of up to 30 % (Kholghy and Kelesidis, 2021) in the predicted evolution of soot volume fraction, number density, and average particle diameters compared to polydisperse PBEs, owing to the rapid attainment of the self-preserving size distribution (Kholghy and Kelesidis, 2021) and an asymptotic collision frequency function therein, regardless of the primary particle polydispersity (Goudeli et al., 2016). Therefore, monodisperse PD models can provide average aggregate characteristics reasonably well at a much lower computational cost than sectional methods or MOM, making them attractive for coupling with multidimensional CFD simulations in realistic geometries.

Another challenge in soot kinetic models, such as POLIMI (Saggese et al., 2015), is that nucleation and surface growth rates are derived directly from detailed gas-phase chemistry. For example, nucleation is typically described through PAH–PAH interactions, which requires prior calculation of PAH concentrations using large kinetic mechanisms that involve many species and reactions (Ranzi et al., 2012). Similarly, surface growth is governed by detailed chemistry, proceeding either through radical-driven pathways or via condensation of PAHs and their fragments onto soot particles (Saggese et al., 2015). As such, both pathways are tied to extensive gas-phase mechanisms, which can be computationally prohibitive for multidimensional CFD simulations in realistic geometries. Traditionally, however, soot formation was described using semi-empirical models, in which nucleation and surface growth were correlated directly to acetylene rather than PAH concentrations (Brookes and Moss, 1999). While far more efficient and straightforward to implement, these correlations were developed for specific flame conditions and lack generalizability across different fuels and operating conditions (Johnson et al., 2021; Fakharnezhad et al., 2026). Thus, both detailed PAH-based and semi-empirical acetylene-based models face critical limitations when applied to practical design applications.

Here, lumped, acetylene-concentration-dependent nucleation (Fakharnezhad et al., 2026) and surface growth (Goudeli et al., 2026) rates derived from reactive MD are implemented into an easy-to-use monodisperse PD model (Fakharnezhad et al., 2026), which accounts for simultaneous nucleation, surface growth, and coagulation. The PD model is coupled with two-dimensional CFD to account for the flame geometry and fluid flow. These MD-derived rates, obtained without assumptions about the contributing reaction pathways or precursor species (Fakharnezhad et al., 2025), do not require simultaneous implementation or use of detailed gas-phase reaction mechanisms. The performance of these MD-derived rates is systematically assessed against established nucleation and surface growth rates from the literature implemented in the developed CFD–monodisperse framework, and the results are benchmarked against a detailed sectional model (Saggese et al., 2015, 2016) and experiments from a burner-stabilized stagnation (BSS) flame (Camacho et al., 2015). The results are also compared to those obtained from the standard semi-empirical Moss–Brookes model for spherical soot nanoparticles implemented in Fluent, which is also based on a monodisperse particle dynamics formulation but assumes coagulation only by full coalescence.

2 Theory

2.1 CFD model

ANSYS Fluent (V 2023 R1) was applied to model soot formation in a BSS laminar ethylene flame, with an unburned mixture consisting of 16.3 mol % ethylene and 23.7 mol % oxygen diluted in argon (Camacho et al., 2015).

2.1.1 Geometry and mesh

This flame was chosen for its well-defined flame geometry, consisting of a porous plug burner and featuring nitrogen shielding for isolation from ambient air and an adjustable cooled stagnation plate with an embedded orifice, facing the incoming flow for sampling (Camacho et al., 2015), with dimensions provided in Fig. 1. The flow geometry is axisymmetric, and so the velocity, temperature, and species composition depend solely on the axial and radial directions. Consequently, the computational domain is reduced to a two-dimensional rectangular section, representing a radial–axial slice (Fig. 1: side 6), which is revolved around the axis of symmetry to model the full cylindrical flame geometry. This 2D domain was discretized into 46 400 quadrilateral cells using a non-uniform mesh. The cell size ranges from ∼0.05 mm in the axial direction near the reactant inlet and stagnation plate up to 0.2 mm away from them.

Figure 1 (a) Schematic of the computational domain with applied boundary conditions and mesh configuration used in the CFD simulation, along with (b) a representative contour plot of the resulting temperature field.

2.1.2 Governing equations

The steady-state Navier–Stokes equations were solved using a pressure-based solver that ensures mass conservation through a pressure correction equation. The equation for conservation of momentum (Batchelor, 2000) is as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }\mathit{v}\right)+\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{v}\mathit{v}\right)=-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\mathit{\tau }}}\right)+\mathit{\rho }\mathit{g}+\mathit{F},\end{array}$$

where ρ is the gas density, v is the velocity vector, p is the static pressure, and $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\mathit{\tau }}}$ is the viscous stress tensor. The terms ρg and F represent the gravitational and external body forces, respectively. The energy equation was also included to capture thermal variations:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}& \left(\mathit{\rho }\left(e+{\displaystyle \frac{v^{\mathrm{2}}}{\mathrm{2}}}\right)\right)+\mathrm{\nabla }\cdot \left(\mathit{\rho }v\left(h+{\displaystyle \frac{V^{\mathrm{2}}}{\mathrm{2}}}\right)\right)=\mathrm{\nabla }\\ & \cdot \left(k_{\mathrm{eff}}\mathrm{\nabla }T-\sum _j{h}_j{\mathit{J}}_j+{\stackrel{\mathrm{‾}}{\mathit{\tau }}}_{\mathrm{eff}}\cdot \mathit{v}\right)+S_{\mathrm{h}}.\end{array}\end{array}$$

The e and h on the left-hand side of Eq. (2) denote the internal energy and the enthalpy, respectively. On the right-hand side, k_eff is the effective thermal conductivity, J_j is the diffusion flux of species j, and S_h accounts for volumetric heat sources including the heat released by chemical reactions.

Species transport was modeled using the volumetric reaction approach (Kee et al., 2003) with the DRM22 reaction mechanism (Kazakov and Frenklach, 1995), including 22 species and 104 reactions. The DRM22 was developed as a reduced form of GRI-Mech 1.2 and has been validated against detailed mechanism predictions of ignition delays and laminar premixed methane-air flame properties over a wide range of conditions, showing deviations typically below 5 % (Kazakov and Frenklach, 1995). This mechanism is further validated here, as discussed in Sect. 3.1.1. The CHEMKIN-CFD chemistry solver was employed, including thermal diffusion effects.

2.1.3 Boundary conditions

At the premixed inlet (Fig. 1: boundary 6), the velocity components, temperature, and species mass fluxes are prescribed. The axial velocity is set to v=0.14 m s⁻¹, and the temperature is fixed at T=473 K (Saggese et al., 2016). The species fluxes account for both convective and diffusive contributions, ensuring mass conservation at the boundary:

$$\begin{array}{}\text{(3)}& \mathit{\rho }v{Y}_k+\mathit{\rho }{Y}_k{V}_k=(\mathit{\rho }v{Y}_k{)}_{z=\mathrm{0}},\end{array}$$

where ρ is the gas density, Y_k is the mass fraction of species k, and V_k is the diffusion velocity. The species balance fluxes that account for both convective and diffusive contributions were obtained from Chemkin-Pro under identical conditions and were imposed as the inlet boundary condition in the Fluent solver. A pure nitrogen inlet (Fig. 1: boundary 3) is defined with an inlet velocity of 0.4 m s⁻¹ and temperature of 473 K (Saggese et al., 2016). The stagnation surface (Fig. 1: boundary 5) is treated as a no-slip wall with a fixed temperature of 500 K (Saggese et al., 2016). At the outlet (Fig. 1: boundary 4), an atmospheric pressure condition is applied, allowing the flow to exit the domain naturally.

2.1.4 Initial conditions

All simulations were first initialized using the Hybrid Initialization method used to generate a physically representative initial flow field through boundary-based interpolation. Then, Fluent simulations were performed under cold-flow conditions (no combustion chemistry considered) until the flow reached steady state. Upon convergence of the cold-flow solution, the volumetric reactions were activated, and a temperature of 2000 K was imposed in the whole region above the burner surface. This temperature patch acts as an ignition source to initiate flame formation in the combustion zone. This two-step procedure ensured numerical stability and provided a physically realistic initial condition for the reacting flow simulations.

2.2 Soot dynamics

The evolution of soot was described using a monodisperse PD framework coupled to the CFD solver in ANSYS Fluent through user-defined functions. Specifically, the evolution of total particle number density N and the total carbon molar concentration C of carbon contained in soot particles (obtained by the soot mass concentration using the molar mass of carbon, MW_c=12 g mol⁻¹) was modeled using CFD transport equations, with nucleation, surface growth, and coagulation rates incorporated as custom source terms through user-defined functions in Fluent. In this implementation, the transport of total particle number density N is described by

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial \left(\mathit{\rho }{u}_{i}N\right)}{\partial x_i}}-{\displaystyle \frac{\partial }{\partial x_i}}\left({\mathrm{\Gamma }}_{k,N}{\displaystyle \frac{\partial N}{\partial x_i}}\right)=S_{N,\mathrm{Nucl}.}+S_{N,\mathrm{Coag}.},\end{array}$$

and the transport of carbon molar concentration C is defined by

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial \left(\mathit{\rho }{u}_{i}C\right)}{\partial x_i}}-{\displaystyle \frac{\partial }{\partial x_i}}\left({\mathrm{\Gamma }}_{k,C}{\displaystyle \frac{\partial C}{\partial x_i}}\right)=S_{C,\mathrm{Nucl}.}+S_{C,\mathrm{SG}}.\end{array}$$

The terms on the left-hand side of Eqs. (3) and (4) represent convection and diffusion in the flow, where u_i denotes the velocity, and Γ_k,N and Γ_k,C are the effective diffusion coefficients for number density and carbon concentration, respectively. The right-hand side of Eq. (4) includes the nucleation term, $S_{N,\mathrm{Nucl}.}$, which describes the formation of new particles, and the coagulation term, $S_{N,\mathrm{Coag}.}$, which describes the reduction in the number of particles due to agglomeration. In Eq. (3), the source terms correspond to carbon addition due to particles newly formed by nucleation, $S_{C,\mathrm{Nucl}.}$, which can be obtained by

$$\begin{array}{}\text{(6)}& S_{C,\mathrm{Nucl}.}=S_{N,\mathrm{Nucl}.}{\displaystyle \frac{n^{\ast }}{N_{\mathrm{Av}}}},\end{array}$$

where n^∗ is the number of carbon atoms within the soot critical nucleus, and N_Av is Avogadro's number. S_C, SG is the accumulation of carbon on existing particles due to surface growth.

The total soot surface area concentration, A, is not explicitly solved in this implementation. Instead, it is reconstructed from N and C using a morphological relationship (Xiong and Pratsinis, 1993):

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{A/N}{a_{\mathrm{0}}}}={\left({\displaystyle \frac{V/N}{v_{\mathrm{0}}}}\right)}^{D_{\mathrm{s}}/\mathrm{3}},\end{array}$$

where a₀ and v₀ denote the surface area and volume of a primary soot particle, $V=\frac{C{\mathrm{MW}}_{\mathrm{c}}}{{\mathit{\rho }}_{\mathrm{soot}}}$ is the total soot volume concentration of the whole particle population, MW_c=12 g mol⁻¹, and D_s is the surface fractal dimension. For spherical particles, D_s=2, while for agglomerates and aggregates, D_s is between 2 and 3. The D_s correlation and the coagulation kernel employed in this study, are detailed in Fakharnezhad et al. (2026).

The source terms for nucleation and surface growth follow established literature models in the baseline framework (as listed in Table 1), where the Moss–Brookes acetylene-based nucleation rate (Brookes and Moss, 1999) is combined with the hydrogen abstraction-carbon addition (HACA) surface growth mechanism (Appel et al., 2000) (Model I). To assess the impact of molecular-level rate expressions, three additional cases were considered: in Model II, the Moss–Brookes nucleation rate was replaced with the MD-derived nucleation rate (Fakharnezhad et al., 2026) while keeping HACA unchanged; in Model III, nucleation is described by the Moss–Brookes equation, while surface growth was substituted with an MD-derived rate (Goudeli et al., 2026); and, in Model IV, both nucleation and surface growth were described using MD-derived expressions. The MD-derived nucleation rate (Model II) was developed from MD simulations of isothermal acetylene pyrolysis at temperatures of 1200–1800 K and high pressures of 189–568 atm in Fakharnezhad et al. (2026). The nucleation rate was determined by monitoring the temporal evolution of the number of incipient soot clusters for various initial acetylene concentrations using a single lumped-rate expression (Table 1: $S_{N,\mathrm{Nucl}.}$ of Model II). This MD-derived rate equation corresponds to the overall nucleation rate and does not require the identification of specific intermediates or tracking their concentration. As such, in contrast to traditional approaches where existing kinetic models are interfaced with CFD, here, a detailed PAH-based mechanism does not need to be employed.

Similarly, the MD-derived surface growth rate, k_MD, SG (Goudeli et al., 2026), in Model III was obtained during isothermal surface growth of incipient soot nanoparticles by acetylene pyrolysis at a comparable temperature (1350–1800 K) and pressure range (106–568 atm). The steady-state particle mass growth rate, ${\left(\frac{\mathrm{d}m}{\mathrm{d}t}\right)}^{\mathrm{SG}}$, was obtained by distinct simulations with different initial acetylene concentrations [C₂H₂] and particle sizes (or surface areas, a), and the rate constant k_MD, SG was determined by the following (Goudeli et al., 2026):

$$\begin{array}{}\text{(8)}& k_{\mathrm{MD},\mathrm{SG}}={\displaystyle \frac{{\left(\frac{\mathrm{d}m}{\mathrm{d}t}\right)}^{\mathrm{SG}}}{{\mathrm{MW}}_{\mathrm{c}}\cdot \left[{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{2}}\right]\cdot a}},\end{array}$$

resulting in an Arrhenius dependence with temperature (Goudeli et al., 2026),

$$\begin{array}{}\text{(9)}& \ln \left(k_{\mathrm{MD},\mathrm{SG}}\right)=\mathrm{3.946}-{\displaystyle \frac{\mathrm{2017.3}}{T}}.\end{array}$$

Equation (8) originates from the HACA equation (Frenklach and Wang, 1991), where $k_{\mathrm{s}}\cdot {\mathit{\chi }}_{\mathrm{soot}}/N_{\mathrm{Av}}$ has been replaced by the MD-derived lumped-rate constant, k_MD, SG (Goudeli et al., 2026).

For comparison, the Moss–Brookes soot model (Brookes and Moss, 1999), an early model developed in 1999, implemented in Fluent (Model V), is also considered as it is widely used in design applications (Wang, 2020). The latter is a semi-empirical, monodisperse formulation coupled with CFD, assuming full coalescence of spherical particles during coagulation in the free-molecular regime. Additionally, a detailed sectional CFD-PBE model from the literature (Saggese et al., 2015, 2016) (Model VI) is included, representing 1D and 2D simulations with full gas-phase chemistry and soot dynamics. The results of this model reported in the literature (Saggese et al., 2015; Saggese et al., 2016) are used to benchmark the proposed reduced framework against detailed sectional simulations. Overall, these six models (see Table 1) encompass a range of physical assumptions. For example, the MD-derived nucleation (Fakharnezhad et al., 2026) and surface growth rates (Goudeli et al., 2026) have been obtained at high pressures and within a narrow temperature range of approximately 1200–1800 K. All are systematically benchmarked against experiments in Sect. 3.

3 Results and discussion

3.1 Validation

3.1.1 Reaction mechanism

The reduced DRM22 mechanism (Kazakov and Frenklach, 1995) was first benchmarked in one-dimensional Chemkin-Pro simulations against a well-characterized premixed ethylene flame (Φ≈1.7) (Gerasimov et al., 2012), where detailed experimental measurements are available, providing a suitable reference for validation under conditions comparable to the flame under investigation. Figure 2 shows the mole fraction profiles of major and intermediate species as a function of the height above the burner (HAB), including (a) H₂O, H₂, and O₂; (b) CO, CO₂, and C₂H₄; and (c) C₂H₂ and CH₄. The predictions are in excellent agreement with measurements (Gerasimov et al., 2012) of O₂ and C₂H₄ consumption and H₂O, CO, H₂, and CO₂ formation (Fig. 2a and b), validating the employed chemical reaction scheme for oxidation. The mechanism also captures the trend and order of magnitude of the mole fraction of C₂H₂ and CH₄ intermediates (Fig. 2c), validating further the DRM22 mechanism for the chemistry of the present premixed ethylene flame.

Table 1 Summary of soot nucleation and surface growth models employed in this study.

Figure 2 Species profiles as a function of the height above the burner, HAB, in the premixed ${\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{4}}/{\mathrm{O}}_{\mathrm{2}}/\mathrm{Ar}$ sooting flame (equivalence ratio, Φ=1.7) obtained by experiments (Gerasimov et al., 2012: symbols) and kinetic modeling (lines).

3.1.2 Computational fluid dynamics (CFD) model

Figure 3 shows the CFD-predicted centerline temperature profile (line) as a function of HAB, showing excellent agreement with temperature measurements in BSS flames (Camacho et al., 2015: symbols). The temperature increases abruptly with HAB up to 1 mm and reaches a maximum of approximately 1700 K between HAB = 1–4 mm, before it drops downstream in the flame (HAB > 4 mm) as it approaches the stagnation plate at HAB = 8 mm. The CFD-obtained temperature contour (Fig. 1b) also nicely reproduces the overall shape of the real flame by experiments (Fig. 1a), confirming that the thermal field is described reliably.

3.2 Soot growth by coupled CFD-PD simulations

The coupled CFD-PD model (Eqs. 4–7) tracks the soot number density, carbon concentration, and surface area evolution, accounting for nucleation, surface growth, and coagulation, as described in Sect. 2.3. Figure 4 shows contours of the soot number density, N (Fig. 4a), and volume fraction, f_v (Fig. 4b), in the sooting flame for burner-to-stagnation plate separation, H_p, of 8 mm, obtained using the Moss–Brookes nucleation rate and HACA surface growth rate (Model I: baseline model). Although a 2D axisymmetric domain was used, the analysis focuses on the flame centerline, consistently with the 1D approximation of BSS flames (Saggese et al., 2016). Near the burner exit (HAB = 0–1 mm), where the temperature of fresh premixed gas starts to increase as it combusts and intermediate species (such as C₂H₂) form, limited soot nucleation may occur. However, the soot precursors still form, and nucleation is not yet dominant. Downstream in the flame, at HAB = 1–5 mm, the temperature reaches its maximum, and ethylene decomposes to acetylene and other soot precursors that result in pronounced nucleation and subsequent particle growth.

Figure 3 Centerline temperature profile as a function of HAB obtained by CFD (line) and experiments (Camacho et al., 2015).

Figure 4 Contours of soot (a) number density, N, and (b) volume fraction, f_v, in the sooting flame with burner-to-stagnation plate separation of H_p=8 mm (predicted using Model I: Moss–Brookes nucleation and HACA surface growth rates).

The soot number density increases rapidly after the onset of nucleation, reaching a maximum of $N=\mathrm{9}\times {\mathrm{10}}^{\mathrm{16}}$ # m⁻³ at HAB ≈3.6 mm, while, at HAB > 3.6 mm, it decreases as coagulation takes over (Fig. 4a). In contrast, the soot volume fraction increases continuously downstream in the flame as particle mass grows through reactions of gaseous reactive species with the soot particle surface, even though their number density decreases by coagulation. The highest soot loading is observed close to the stagnation plate (Fig. 4b). Since the probe was positioned along the centerline but not explicitly included in the 2D simulation, its effect was accounted for by shifting the sampling position upstream along the axial direction following the correction of Saggese et al. (2016). The required spatial shift ranges from 1 to 1.6 mm, depending on the burner-to-plate separation, as obtained from the probe effect analysis (Saggese et al., 2016). For this case, the highest soot loading is ∼100 ppm after applying a 1.5 mm shift. This correction was applied to all results to ensure consistency with experiments.

Figure 5 shows the soot number density and soot volume fraction predicted by the CFD-coupled monodisperse PD model using different nucleation and surface growth rates (Table 1: Models I–IV), as well as by the Moss–Brookes soot model implemented in Fluent (Model V) at probe locations, H_p, of 8 and 4 mm above the burner. The results are compared to a 2D CFD–sectional PBE model by Saggese et al. (2016; Table 1: Model VI), who reported N and f_v in these two probe locations, and to experiments (Camacho et al., 2015) at the same H_p of the BSS ethylene flame.

Figure 5 (a, c) Soot number density (N) and (b, d) volume fraction (f_v) at (a, b) H_p=8 and (c, d) 4 mm obtained by CFD-coupled monodisperse PD with various nucleation and surface growth rates (Models I to IV) and by the Moss–Brookes soot model implemented in Fluent (Model V). The results are compared to a 2D sectional model (Saggese et al., 2016): Model VI) and experiments (Camacho et al., 2015) in the same BSS ethylene flame.

At H_p=8 mm (Fig. 5a), the CFD-coupled monodisperse PD model with the Moss–Brookes nucleation rate and HACA surface growth rate (Model I, baseline) yields $N=\mathrm{3.66}\times {\mathrm{10}}^{\mathrm{16}}$ # m⁻³, overpredicting the experimental value of 1.66×10¹⁶ # m⁻³. However, the prediction is comparable (within ∼6 %) to that reported by Saggese et al. (2016) (Model VI; $N=\mathrm{3.44}\times {\mathrm{10}}^{\mathrm{16}}$ # m⁻³). Replacing the nucleation and/or surface growth rates with the MD-derived expressions (Models II–IV) results in higher N values (up to about a factor of 3) relative to Model I, while the differences among all models, including Model V, remain within the experimental uncertainty at this H_p.

The baseline model (Model I) predicts f_v=54 ppb (Fig. 5b), which is about twice that obtained from the sectional model (Model VI; f_v=26.1 ppb; Saggese et al., 2016) but within 1 order of magnitude of the measured f_v of 7 ppb. This indicates that the present framework reproduces soot mass loading with an accuracy comparable to that of detailed sectional modeling. When nucleation is described by the MD-derived rate (Model II), f_v decreases down to 28.4 ppb, matching the sectional model prediction within ∼4 % and improving agreement with the measured f_v by ∼57 % compared to Model I. Accounting for the MD-derived surface growth rate (Models III and IV) decreases the f_v down to 5.7 (Model III) and 3.8 ppb (Model IV), below the measured soot mass loading. The underprediction of f_v in Models III and IV, which incorporate the MD-derived surface growth rate, can be attributed to their substantially lower surface growth contribution compared to the baseline model using the HACA mechanism. The carbon addition rate from surface growth, ${\left(\frac{\mathrm{d}C}{\mathrm{d}t}\right)}^{\mathrm{SG}}$, in these models is almost 2 orders of magnitude lower than in Model I (Fig. S1 in the Supplement), especially across the main soot-forming region (1–7 mm). This reduction in carbon addition explains the lower soot mass predicted by the MD-derived surface growth models.

The Moss–Brookes soot model (Model V) yields $N=\mathrm{5.11}\times {\mathrm{10}}^{\mathrm{16}}$ # m⁻³ (Fig. 5a), which is comparable to all of the other models (Models I–IV) considered, yet the corresponding f_v (4572 ppb) is nearly 3 orders of magnitude higher than the measured f_v and the highest among all models. This overprediction stems from the excessive surface growth contribution (Fig. S1), which is nearly 3 orders of magnitude larger than in the baseline case. These results confirm that surface growth dominates soot mass formation and strongly influences f_v predictions, consistently with the literature (Harris and Weiner, 1985).

At H_p=4 mm (Fig. 5c, d), the baseline model yields $N=\mathrm{3.44}\times {\mathrm{10}}^{\mathrm{16}}$ # m⁻³ and f_v=3.43 ppb, about 1 order of magnitude higher than those obtained by the sectional model ($N=\mathrm{7.82}\times {\mathrm{10}}^{\mathrm{15}}$ # m⁻³ and f_v=0.24 ppb) and 2 and 3 orders of magnitude higher than the experimentally measured N of 1.22×10¹⁴ # m⁻³ and f_v of 0.002 ppb. Models II to V show comparable trends for both soot number density and volume fraction, as also observed at H_p=8 mm (Fig. 5a, b). Overall, all models substantially overpredict the experimental data at H_p=4 mm, with deviations of approximately 2 to 3 orders of magnitude for both N and f_v.

The performance of the reduced models is assessed across a broader range of burner-to-stagnation plate separations investigated experimentally, as shown in Fig. 6. The results are now compared against the 1D sectional model of Saggese et al. (2015), which employed the same soot sub-model as in 2D simulations (Fig. 5) but covers a wider range of H_p values from 4 to 12 mm. Figure 6 shows (a) the soot number density and (b) soot volume fraction as a function of H_p, as predicted by Models I–V, the 1D sectional model (Model VI: Saggese et al., 2015), and experiments (Camacho et al., 2015). The sectional model (Model V) shows the best agreement with number density and soot volume fraction measurements, outperforming all reduced models, as expected. The soot number density predicted by Models II–V is practically constant across all H_p values. Model I suggests a slight decrease in the soot number concentration with H_p, which is qualitatively consistent with the detail sectional model and in reasonable agreement with measurements (circles) for H_p≥6 mm. At H_p=4 mm, however, all models overpredict N by 2–3 orders of magnitude. This discrepancy has been partially attributed to simplifications made in employed nucleation mechanisms (Hou et al., 2019), even though it may also stem from the presence of sub-2 nm soot nuclei (Fakharnezhad et al., 2026) undetected during measurements by standard butanol- or water-based condensation particle counters (Camacho et al., 2015).

Figure 6 (a) Soot number density, N, and (b) volume fraction, f_v, as a function of H_p as obtained by CFD-coupled monodisperse PD with various nucleation and surface growth rates (Models I to IV) and by the Moss–Brookes soot model implemented in Fluent (Model V) (filled symbols). The results are compared to 1D sectional model (Model VI: orange line; Saggese et al., 2015) and soot measurements (circles) in ethylene premixed flames (Camacho et al., 2015). The error bars for the measured f_v (not shown) are within the symbol size.

Even though all models suggest soot number concentrations comparable to each other (within 1 order of magnitude), they result in distinct soot volume fractions that span over up to 5 orders of magnitude (Fig. 6b). Both Models I and II yield comparable f_v predictions, in good agreement with the sectional model (Saggese et al., 2015) and measurements (Camacho et al., 2015) within up to 1 order of magnitude for H_p≥6 mm. In contrast, the Moss–Brookes soot model (Model V) overpredicts the soot volume fraction by nearly 3 orders of magnitude, whereas Models III and IV underpredict the measurements by about 2 orders of magnitude. As previously shown at H_p=8 mm (Fig. S1), this underprediction is caused by the lower carbon addition rate from surface growth, which primarily arises from the fact that the MD-derived effective surface growth rate constant (Goudeli et al., 2026) is approximately 3 times lower than that predicted by the HACA-based correlation of Appel et al. (2000). This slower surface growth rate limits carbon accumulation in the soot phase, resulting in lower soot volume fractions. At H_p=4 mm, all reduced models overpredict the measurements by 2–5 orders of magnitude. Such a significant overprediction is not unique to the low-order modeling framework (Models I–V) as it is observed even with detailed sectional models (Model VI), which has been attributed to incomplete description of soot nucleation PAH chemistry or the assumption of a non-variable density for nascent soot (Hou et al., 2019).

Figure 7 shows the soot mobility diameter, d_m, predicted by Models I–IV using the CFD-coupled monodisperse PD framework developed in this study. The volume-equivalent diameter, d_v, is also shown as predicted by the Moss–Brookes soot model (Model V), which assumes coagulation by full coalescence, neglecting the formation of particles with fractal-like structure (i.e., $d_{\mathrm{p}}=d_{\mathrm{m}}=d_{\mathrm{v}}$). The model predictions are compared with the average d_m from the detailed 1D sectional model (Model VI: orange line) derived from the reported mobility size distributions in Saggese et al. (2015: Fig. 11), as well as with the average d_m obtained by scanning mobility particle sizer (SMPS) measurements in premixed ethylene flames (Camacho et al., 2015: symbols). Even though the present PD model does not account for the agglomerate and/or aggregate polydispersity, the mobility-based self-preserving size distribution can be derived based on the average d_m (Fig. 7) and the mobility-based self-preserving geometric standard deviation of soot (${\mathit{\sigma }}_{\mathrm{g},\mathrm{m}}=\mathrm{1.48}\pm \mathrm{0.03}$), obtained by discrete element method simulations (Kelesidis and Goudeli, 2021: Fig. 2). All models capture the expected increase in particle size with increasing burner-to-stagnation plate separation. Similarly to soot volume fraction and number concentration (Fig. 6), the 1D sectional model predicts most accurately the mobility diameter measurements. Among the reduced models, the MD-informed nucleation model (Model II) shows the best agreement with the measured mobility diameters and with the 1D sectional model (Model VI) predictions, especially for H_p≥6 mm, followed by the baseline Model I, which predicts slightly larger diameters at all H_p values. Both Models I and II predict the formation of fractal-like agglomerates (as shown in Fig. S2 for Model II) with $\mathrm{2}\le D_{\mathrm{s}}\le \mathrm{3}$, capturing the fractal-like morphology of soot aggregates produced in flames. Models III and IV, which employ the MD-derived surface growth rate, suggest limited particle growth (<3 nm) and nearly spherical particle formation (D_s=2, based on the D_s correlation used in Fakharnezhad et al., 2026) even at large separation distances (H_p≥8 mm). In contrast, Model V overestimates particle size substantially due to the much higher carbon addition during surface growth compared to HACA-based or MD-derived correlations. This model inherently predicts spherical particles as it explicitly assumes coagulation by full coalescence (Brookes and Moss, 1999).

Figure 7 Soot mobility diameter (d_m) predicted by the CFD-coupled monodisperse PD framework (Models I–IV) and volume-equivalent diameter (d_v) predicted by the Moss–Brookes soot model (Model V; full coalescence) (filled symbols). The model predictions are compared to the 1D sectional model (Model VI, Saggese et al., 2015: orange line) and to the average d_m from SMPS measurements in premixed ethylene flames (Camacho et al., 2015).

Among the reduced models (Models I–IV), the baseline framework (Model I), accounting for the semi-empirical nucleation rate and the HACA surface growth rate, captures the overall trend of the soot number density (Fig. 6a), volume fraction (Fig. 6b), and particle size (Fig. 7). Accounting for the MD-derived nucleation rate, Model II improves the d_m predictions and nicely matches the soot volume fraction and SMPS measurements in this flame. Model II predicts comparable N (within 1 order of magnitude) with the detailed sectional model across all H_p but is outperformed by Models I (Moss–Brookes nucleation rate and HACA-predicted surface growth rate) and V (Moss–Brookes nucleation and surface growth rates). It should be noted that, even though the Moss–Brookes nucleation rate in Model I and the lumped MD-derived nucleation rate employed in Model II cannot capture the level of chemical detail accounted for in sectional models interfaced with chemical kinetic reaction mechanisms, they do capture f_v and d_m measurements with reasonable accuracy, especially at HAB > 5 mm, at significantly lower computational cost.

4 Conclusions

A multiscale modeling framework was developed that couples a monodisperse particle dynamics (PD) model with computational fluid dynamics (CFD), describing soot formation and growth in a premixed ethylene burner-stabilized stagnation (BSS) flame. The coupled PD-CFD model, accounting for simultaneous nucleation, surface growth, and coagulation, is used to evaluate the effect of nucleation and surface growth rate equations developed semi-empirically or by molecular dynamics (MD) on the soot number density, volume fraction, and soot mobility diameter. Accounting for the semi-empirical Moss–Brookes nucleation rate and the hydrogen abstraction-carbon addition (HACA) surface growth rate, Model I (baseline) reproduces the soot number density, volume fraction, and particle size measurements reasonably well across the BSS flame. Incorporating the MD-derived nucleation rate and the HACA surface growth model (Model II) further improved soot volume fraction and mobility diameter predictions compared to those obtained by semi-empirical rates, nicely matching measurements and a sectional model with the detailed chemical kinetic mechanism (Saggese et al., 2015) in the post-flame region (H_p≥6 mm). Model II, encompassing the MD-informed nucleation rate, reproduced mobility diameter measurements most accurately among all reduced models and yielded soot volume fractions within 1 order of magnitude from experimental data, outperforming reduced models employing the MD-derived or Moss–Brookes surface growth rates. Even though the developed MD-informed CFD-PD framework does not include the detail accounted for in chemical kinetic mechanisms, it captures the measured f_v and d_m reasonably well, resulting in even comparable predictions to those by a detailed sectional model in the post-nucleation flame region at significantly lower computational cost. Even though this framework is benchmarked in a premixed ethylene BSS flame where the centerline temperature is mostly within the range in which the MD rate constants were derived, it provides a promising foundation for predictive, design-oriented simulations of combustion and aerosol synthesis systems. Additional benchmarking would be necessary before application to different flame configurations.

Appendix A: Nomenclature

A	Total soot surface area concentration
Ap	Pre-exponential factor
a	Surface area of a soot particle
a0	Surface area of a primary soot particle
C	Total carbon molar concentration
Cα	Model constant for soot inception rate
Cγ	Model constant for surface growth rate
dm	Soot mobility diameter
dp	Soot primary particle diameter
dv	Soot volume-equivalent diameter
Ds	Surface fractal dimension
e	Internal energy
E	Activation energy
F	External body force
fv	Soot volume fraction
g	Gravitational acceleration vector
h	Enthalpy
HAB	Height above burner
Hp	Burner-to-stagnation plate separation
Jj	Diffusion flux of species j
keff	Effective thermal conductivity
kf	Forward reaction rates of soot with H
kb	Backward reaction rates of soot with H
$k_{\mathrm{MD},\mathrm{Nucl}.}$	MD-derived nucleation rate constant
kMD, SG	MD-derived SG rate constant
ks	Surface reaction rate coefficient
M	Soot mass fraction
MWc	Molar mass of carbon
n	Temperature exponent
n∗	Number of carbon atoms within
	critical nucleus of soot
N	Total particle number density
NAv	Avogadro number
p	Static pressure
R	Universal gas constant
Sh	Volumetric heat sources
T	Temperature
Tα	Activation temperature of soot inception
Tγ	Activation temperature of soot surface growth rate
v	Axial velocity
v	Velocity vector
v0	Volume of a primary soot particle
Vk	Diffusion velocity
Yk	Mass fraction of species k
Greek Letter
α	Fraction of available sites
Γk	Effective diffusion coefficient
χsoot−H	Total number of available surface sites
ρ	Gas density
ρsoot	Soot density
$\stackrel{\mathrm{‾}}{\mathit{\tau }}$	Viscous stress tensor
σg,m	Mobility-diameter-based geometric
	standard deviation