Sensitivity Analysis of Key Parameters for Population Balance Based Soot Model for Low-Speed Diffusion Flames

Abstract

In this article, the evolution of in-flame soot species in a slow speed, buoyancy-driven diffusion flame is thoroughly studied with the implementation of the population balance approach in association with computational fluid dynamics (CFD) techniques. This model incorporates interactive fire phenomena, including combustion, radiation, turbulent mixing, and all key chemical and physical formation and destruction processes, such as particle inception, surface growth, oxidation, and aggregation. The in-house length-based Direct Quadrature Method of Moments (DQMOM) soot model is fully coupled with all essential fire sub-modelling components and it is specifically constructed for low-speed flames. Additionally, to better describe the combustion process of the parental fuel, ethylene, the strained laminar flamelet model, which considers detailed chemical reaction mechanisms, is adopted. Numerical simulation is validated against a self-conducted co-flow slot burner experimental measurement. A comprehensive assessment of the effect of adopting different nucleation laws, oxidation laws, and various fractal dimension and diffusivity values is performed. The results suggest the model employing Moss law of nucleation, modified NSC law of oxidation, and adopting a fractal dimension value of 2.0 and Schmidt number of 0.9 yields the simulation result that best agreed with experimental data.

1. Introduction

Soot particles are fine carbonaceous particulates that are formulated through multiple physical binning and chemical reaction processes of intermediate combustion products. They are invariably incepted and evolved in the flame region, and subsequently dispersed in the smoke layer. The soot formation mechanisms begin with the soot nuclei inception from precursor species and it follows by the formation of lager fractal object via particle-particle collision events, including coalescence and aggregation, as well as via surface growth reaction of the soot particulates. It ends with the breakage of the clusters into fragments that are eventually consumed via oxidation processes [1,2]. Soot particulates can be carcinogenic and they endanger our lives through toxication, suffocation, or causing permanent damage to our respiratory systems, such as trachea, blood vessels, and lungs [3]. As a strong emitter and absorber of radiative heat, soot particulates significantly affect the global heat transfer process, both in the inflame and surrounding regions [2,4]. This is especially hazardous in compartment fire scenarios, where the inhalation of the toxic compounds, including flame generated gas mixture and soot particulates, combined its ability to obscure visibility through a combination of absorption and light scattering, can significantly interfere with the evacuation process [5]. The recent development of bio-based fire retardants material could significantly reduce the generation of fire and smoke [6,7], however, there are still a considerable amount of other combustible materials, such as wood and plastic consumables, which could potentially pose a threat to building fire safety. Hence, the evolution of soot particles is of great interest to both research and industrial communities. Is useful in the design of fire protection engineering systems, such as smoke extraction, smoke curtain, pressurised zone, and sprinkler systems. Accordingly, this work aims to investigate the formation and destruction processes of soot particulates in a non-premixed buoyant flame that is commonly encountered in compartment or wildland fires.

Flame generated soot particles have been commonly studied via experimental approaches. For instance, soot particulates can be locally collected using thermophoresis sampling technique, and then be visualized with transmission electron microscope (TEM) imaging [8]. Soot morphological information, e.g., projected area, radius of gyration, fractal dimension of the aggregates, primary particle size etc., can be collectively characterised. The second type of diagnostic approach is the optical based non-intrusive measuring technique, e.g., laser-induced incandescence (LII), elastic light scattering (ELS), small-angle neutron scattering (SANS), and small-angle X-ray scattering (SAXS), which provides insight into a wide variety of soot related properties, including size, structure, and composition of particulates [9,10,11,12]. However, spatial resolution of the measurement and the potential perturbation of the very combustion process under investigation when intrusive approaches are used limited the experimental approaches. Most of the optical and laser-based diagnostics are dependent on optical properties of the materials, which are very difficult characterise. Numerical modelling of soot particles using the population balance approach (PBA), on the other hand, can deliver informative data that are not assessable by experimental measurements alone, such as displaying the full spectrum of particle size and number density distribution for any location and time. Therefore, the use of numerical approaches can provide more theoretical insight that can significantly enhance our understanding of the soot formation mechanisms and the dispersion of soot particles for the research communities.

Numerical simulation of soot involves the appropriate description of the evolution of the fine hydrocarbon particulates as a result of incomplete combustion, which were subsequently immersed or recirculated back into the flammable gas mixture. This phenomenon could be resolved on various scales depending on the aim and focus of the investigation. For example, one conventional method that is commonly used in bushfire or compartment fire modelling is the macroscale level approach, in which only the hydrodynamics behaviours of the phase are resolved, and no more than a soot particle concentration (i.e., mass/volume fraction) term is resolved. However, such macroscale models were initially designed to assess industrial-scale fires and are usually limited to only simulate smoke movement, flow velocities, and temperatures changes due to soot radiation. Furthermore, the deviation of the predicted quantities is often beyond 20% when comparing with experimental measurements due to the lack of detail description of the actual physical and chemical processes that are involved [13]. The microscale Eulerian–Lagrangian approach, on the other hand, tracks the motion of discrete elements by considering the interaction of discrete particles within the computational domain. It provides more fundamental insights into the detailed description of the soot species evolution, thus enabling the deliver simulation with improved prediction accuracy [14]. Nevertheless, it is limited to microscopic level due to its immense requirement for small mesh sizes and it is not practical for medium-to large-scale fire modelling for the enormous computational cost. In practical engineering applications, most of the soot models used at present in commercial computational fluid dynamic (CFD) packages resolve the issue from an intermediate level, i.e., mesoscale, which focus on the local structural development and are aimed to investigate the evolvement of cluster groups e.g., droplets, bubbles, and particles [15,16,17,18,19]. The mesoscale model categorised all sooting processes into the consideration of critical mechanisms and it limits its scope to predict key soot-related properties with reasonable accuracy, whilst notably reducing the computational effort required [20,21,22]. The Direct Quadrature Method of Moments (DQMOM) soot model presented in this paper, therefore, is targeting to resolve soot evolution on a mesoscale level by taking into account all necessary soot processes that are involved the formation pathway, namely particle inception, aggregations, coagulation, surface growth, and oxidation.

One of the major drawbacks of common mesoscale soot models is the difficulty in quantifying case dependent properties. In practice modelling, such properties are often represented by an empirical value or are approximated by empirical formulations. This assumption or simplification, if not carefully justified, can generate large numerical errors. In fire modelling, the particulate size of soot species is one of the properties that has been oversimplified by most current mesoscale models. It is generally agreed that the soot particulate size can vary between a broad range, depending on various aspects, such as parent fuel composition, burner configuration, ambient condition, etc. The two-equations based semi-empirical model, e.g., Moss–Brookes soot model used in Fluent, however, considers a uniform soot particle size distribution (PSD) that is set based on a predefined empirical data, i.e., 35 nm [20]. The population balance based Standard Method of Moments (SMM) soot model evaluates particle size distribution with the use of a presumed log-normal function that is based on experimental data collected under certain flaming conditions [23]. The quadrature method of moment (QMOM) is an attractive alternative, which effectively resolves the enclosure by adopting a moment construction and inversion algorithm. This mathematical approach closed the population balance equation with limited computational burden and it proved to be numerically stable in most of its applications. Nevertheless, it is formulated based on a presumed particle size distribution method that is similar to the concept of SMM. The abovementioned numerical soot models could not directly evaluate the distribution of particle size that varies according to the combustion conditions, which could potentially lead to inaccuracies in determining the other soot properties that are particle size dependent.

The Direct Quadrature Method of Moments (DQMOM) that was developed by Marchisio and Fox [24] directly solves the nodes and weights of the quadrature approximation of the number density function for an internal coordinate, e.g., particle sizes. This novel approach enables the prediction of the real-time and spatial size distribution and, hence, further improves the accuracy of soot mass fraction predictions. Previous numerical studies that implemented the DQMOM approach in soot modelling have showcased its competency to provide promising results with moderate computational burden [25,26,27,28]. However, those studies are mainly focused on the investigation of inertia-driven highly turbulent flame with relatively high fuel/oxidiser injection rate, for example, with 50 ms⁻¹ fuel injection rate. A DQMOM-based soot model targeting low-speed buoyancy-driven diffusion flame, which corresponds better to natural fire occurrences in open configuration, e.g., bushfire and compartment fire, has not been formulated and tested. Moreover, it has been reported that DQMOM algorithm can be numerically unstable and difficult to implement, i.e., generating unrealizable (negative) weights and nodes, if the system is not cautiously designed. Alternative approaches that were proposed by Mueller et al. and Chittipotula et al. [28,29] could potentially resolve the numerical issue. However, the alternative approaches often require solving additional sub-modules that are integrated to the framework or correction algorithms, which could significantly increase the computational cost or need the cumbersome coupling of the CFD code with an external mathematical package. This challenge could become more predominant when the investigated case has relatively small spatial-diffusion terms [29], like the flame of the interest in this work. Therefore, it is critical to systematically investigate the effect of all modelling aspects to ensure the numerical robustness and simulation accuracy.

The information about soot particle size distribution could also aid the design of fire protection systems in practical engineering. Among previous works on smoke control systems [30], it is discovered that the smoke detection system is mainly categorised into three types: (i) ionisation smoke alarms, (ii) light obscuration photoelectric alarms, and (iii) light scattering photoelectric alarms. When the soot particles reach a particular size, it will trigger the alarms by either (i) providing a voltage signal for ionisation alarms, (ii) blocking the light signals to the receiver for obscuration type photoelectric alarms, or (iii) changing the light scattering for scattering type photoelectric alarms. Nevertheless, the current drawbacks of the use of numerical models to predict the detection of smoke particulates are lacking in consideration of the particle size distribution. In essence, the alarms in numerical codes such as FDS is trigger by either gas temperature or smoke volume fraction, which often lead to inaccurate response time. With the application of the DQMOM approach, providing an enhanced understanding of soot formation and destruction processes, especially for the particle size distribution, will contribute to the prediction of fire protection systems.

In light of the above gaps in knowledge, this work will present, for the very first time, a detailed parametric study on the effects of key modelling aspects of an in-house soot model that is built based upon the DQMOM algorithm that enables the real-time tracking of the evolution of PSD. The proposed soot model is fully coupled with a gas-phase combustion model incorporating detailed chemistry, and it is specifically formulated for a for a low-speed buoyant diffusion flame configuration with optimised parameters. The results of the present model are validated against measurements from a self-constructed experiment in the Combustion Laboratory, School of Mechanical and Manufacturing Engineering, the University of New South Wales. Comprehensive assessments will be performed on a collection of key parameters within the DQMOM framework towards the performance of soot modelling for turbulent diffusion flames, and the objectives can be summarised, as follows:

(i)

comparison of three nucleation laws, including Moss, Leung, and Fairweather, and discuss their effectiveness in predicting the changes in particle size due to the particle generations;

(ii)

study of two different oxidation law including modified NSC and Said to investigate the appropriate soot particle reduction or disruption mechanisms for DQMOM;

(iii)

investigating the influences of fractal dimensions towards the aggregation and surface growth mechanisms, as well as the effect on particle size distribution;

(iv)

study the effect of diffusion coefficient for soot quantities and number density towards the dispersion of soot particles within the computational field; and,

(v)

provide a deeper understanding of the soot formation process, in particular, within the flaming/soot nucleating region.

2. Mathematical Model

Numerical simulations have been carried out using an in-house population balance approach (PBA) DQMOM code coupled with Fluent solver (i.e., version 19.1). This computational fluid dynamics model incorporates sub-modelling components, including soot formation, turbulence, combustion, and radiation heat exchange to describe the fire phenomena. To take into account the highly irregular nature of the combustion process, several assumptions are made, including (i) the low Mach number flow equations are considered, (ii) the thermal-physical properties are constant, (iii) the ratio between mass and thermal diffusivity (Lewis number) is unity, (iv) no external body force and heat source are applied towards the flow, and (vi) non-adiabatic non-premixed combustion of the energy equation is adopted.

2.1. Governing Equations

The fully coupled, non-linear, and interactive physical and chemical processes that are involved in the turbulent reacting gas mixture caused by the pulsating diffusion flame is modelled using computational fluid dynamics techniques. With the application of the Favre averaging approach, the following governing equations, i.e., the continuity, Navier–Stokes, and scalar transport equations are utilised in this simulation study are given as:

$$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial }{\partial x_i}\overline{\rho }{\tilde{U}}_i=0$$

(1)

$$\frac{\partial }{\partial }\overline{\rho }\widetilde{U_i}+\frac{\partial }{\partial }{\tilde{U}}_i{\tilde{U}}_j+\frac{\partial }{\partial }\overline{\rho }\widetilde{u_i{u}_j}=-\frac{\partial \overline{p}}{\partial x_i}+\left(\nu +{\nu }_t\right)\frac{{\partial }^2\overline{\rho }{\tilde{U}}_i}{\partial x_j\partial x_j}$$

(2)

$$\frac{\partial }{\partial }\overline{\rho }{\tilde{\Phi }}_{\alpha }+\frac{\partial }{\partial }\overline{\rho }{\tilde{\Phi }}_{\alpha }{\tilde{U}}_i+\frac{\partial }{\partial }\overline{\rho }\widetilde{{\varphi }_{\alpha }{u}_i}=\left(\Gamma +{\Gamma }_t\right)\frac{{\partial }^2\overline{\rho }{\tilde{\Phi }}_{\alpha }}{\partial x_i\partial x_i}+\widetilde{S\left(\Phi \right)}$$

(3)

where $\overline{\rho }$ and $\overline{p}$ are the mean density and pressure, $\nu$ and ${\nu }_t$ is the kinematic and turbulent viscosity, ${\tilde{U}}_i$ and $u_i$ are the Favre-averaged and fluctuation value of the mean fluid velocity, ${\tilde{\Phi }}_{\alpha }$ and ${\varphi }_{\alpha }$ are the mean and fluctuation value of the $\mathsf{\alpha }$th scalar, $\Gamma$ and ${\Gamma }_t$ are molecular and turbulent diffusivity, and $\widetilde{S\left(\Phi \right)}$ is the Favre-averaged chemical reaction source term.

2.2. Turbulence and Combustion Modelling

Literature works have suggested that the standard $k$-$\epsilon$ model is usually suitable for modelling non-premixed diffusion and non-swirling flames, like the one under investigation in this work [25,26,27,28]. The chemical reaction source term in the transport equations of the involved reacting scalars is determined by the strained laminar flamelet approach, in which the combustion chemistry is a pre-assumed probability density function (pre-PDF) of the mixture fraction ($f$) and the scalar dissipation ($\mathsf{\chi }$). In essence, the mixture fraction governs the amount of the fuel mixture in each control volume element in the simulation domain. The scalar dissipation is a term that is introduced to describe the strain and extinction of the flame, in which the magnitude of this quantity depicts the departure of the combustion process away from its chemical equilibrium [31]. It should be noted that, in the present work, the GRI-MECH 3.0 detailed chemical reaction mechanisms, which includes 325 reaction steps and 53 chemical species [32], was implemented to formulate the flamelet library for the strained laminar flamelet model, with ethylene (C₂H₄) being selected as the parental fuel. Such an approach of resolving turbulence-chemistry interaction has been demonstrated in previous studies to provide a reasonable result with moderate computational burden [18,33].

Figure 1 presents the post-processed flamelet profile results for major and minor species, respectively, at near-equilibrium and near-extinction scalar dissipation rates, being generated for 20 different scalar dissipation rates that range from 10⁻⁶ to 180. As can be seen, the key significant combustion products are oxygen, carbon monoxide, carbon dioxide, water vapour, hydrogen, and acetylene, while minor species include hydrocarbon compounds, oxygen, and hydrogen molecules.

2.3. Direct Quadrature Method of Moments (DQMOM) Model

The DQMOM model is written in numerical code and then adapted into the CFD framework. It fully resolves the evolution of the soot particulates in turbulent flames and its size distribution to give a better representation of the soot formation and disruption mechanisms by the DQMOM model. The population balance equation for the Favre-averaged number density $\tilde{n}$ is firstly described, similar to the transport equation for moments, as:

$$\frac{\partial }{\partial t}\left(\overline{\rho }\tilde{n}\right)+\frac{\partial }{\partial x_j}\left(\overline{\rho }\widetilde{u_j}\tilde{n}\right)=\frac{\partial }{\partial x_j}\left[\left({\Gamma }_T+\Gamma \right)\frac{\partial \tilde{n}}{\partial x_j}\right]+{\overline{S}}_{\tilde{n}}$$

(4)

where $\tilde{n}$ is the Favre-averaged number density function (NDF) of the internal coordinate vector $\xi$ and ${\overline{S}}_{\tilde{n}}$ is the source terms taken into the consideration of all chemical and physical mechanisms that are involved in soot evolution.

The main idea behind this methodology is to resolve the closure by introducing a quadrature approximation to evaluate the development of the moments of NDF. For a length-based system with mono-variable $\left(\xi =L\right)$, the NDF can be expressed as:

$$\tilde{n}\approx {\displaystyle \sum }_{\alpha =1}^N{\omega }_{\alpha }\delta \left[L-L_{\alpha }\right]$$

(5)

where ${\omega }_{\alpha }$ and $L_{\alpha }$ are the weight and abscissas of node $\alpha$ of the quadrature approximation, $N$ is the number of node points defined in the distribution of the internal coordinate, and $\delta$ indicates the Dirac delta function. Various soot-related properties can be subsequently derived from the $k$th moment of NDF. For example, with a length-based NDF, the particle number density, mean particle size, and volume fraction corresponding to the zeroth (${\tilde{m}}_0$), first (${\tilde{m}}_1$), and third moment (${\tilde{m}}_3$), and can be resolved as:

$${\tilde{m}}_k={{\displaystyle \int }}_0^{+\infty }\tilde{n}{L}^{k}dL\approx {\displaystyle \sum }_{\alpha =1}^N{\omega }_{\alpha }{L}_{\alpha }^k$$

(6)

Meanwhile, the weight and abscissas are solved through the DQMOM transport equations:

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\omega }_{\alpha }\right)+\frac{\partial }{\partial x_j}\left(\overline{\rho }\widetilde{u_j}{\omega }_{\alpha }\right)=\frac{\partial }{\partial x_j}\left[{\Gamma }_T\frac{\partial {\omega }_{\alpha }}{\partial x_j}\right]+{\alpha }_{\alpha }$$

(7)

$$\frac{\partial }{\partial t}\left(\overline{\rho }{ℒ}_{\alpha }\right)+\frac{\partial }{\partial x_j}\left(\overline{\rho }\widetilde{u_j}{ℒ}_{\alpha }\right)=\frac{\partial }{\partial x_j}\left[{\Gamma }_T\frac{\partial {ℒ}_{\alpha }}{\partial x_j}\right]+b_{\alpha }$$

(8)

where ${ℒ}_{\alpha }={\omega }_{\alpha }\times L_{\alpha }$ is defined as the $\alpha$th weighted abscissa, and ${\alpha }_{\alpha }$ and $b_{\alpha }$ are the source terms for weight and weighted abscissa, respectively. For a numerical system with $N=2$, the source terms of ${\alpha }_{\alpha }$ and $b_{\alpha }$ can be determined via solving the following rank four linear system of the first four moments:

$${\alpha }_1+{\alpha }_2={\overline{S}}_0$$

(9)

$$b_1+b_2={\overline{S}}_1$$

(10)

$$-L_1^2{\alpha }_1-L_2^2{\alpha }_2+2L_1{b}_1+2L_2{b}_2={\overline{S}}_2+{\overline{C}}_2$$

(11)

$$-2L_1^3{\alpha }_1-2L_2^3{\alpha }^2+3L_1^2{b}_1+3L_2^2{b}_2={\overline{S}}_3+{\overline{C}}_3$$

(12)

where ${\overline{C}}_k$ is a correction term of diffusion in real space expressed as:

$${\overline{C}}_k=k\left(k-1\right){\displaystyle \sum }_{\alpha =1}^N{L}_{\alpha }^{k-2}{\omega }_{\alpha }{\Gamma }_T\frac{\partial D_{\alpha }}{\partial x_j}\frac{\partial D_{\alpha }}{\partial x_j}$$

(13)

and ${\overline{S}}_k$ is the integrated source term of the $k$th moments as:

$${\overline{S}}_k={{\displaystyle \int }}_{-\infty }^{+\infty }{L}^k{\overline{S}}_{\tilde{n}}dL$$

(14)

In this particular soot modelling application, the integrated source term is evaluated as a summation of several contributions that account for critical physical and chemical mechanisms that are involved in soot evolution, e.g., nucleation, coagulation, surface growth, and oxidation. Each contribution can be expressed as a function of weight, weighted abscissas and the reaction rate of the mechanism. The formulation of the individual terms in the integrated source term is further described in detail in the following section.

2.4. Soot Formation Kinetics

In practical soot modelling, the following mechanisms are often defined as the key sooting process that is involved in the soot formation pathway, namely nucleation, aggregation, surface growth, and oxidation. Therefore, those mechanisms are considered in the formulation of the integrated source term, ${\overline{S}}_k$.

2.4.1. Nucleation

The in-flame soot formation commences with the inception process that solid-phase soot nuclei formed from the basic structural units (BSUs) consist of fuel fragments as a result of pyrolysis of the parent fuel [34]. The nucleation is, so far, the most complicated and least-understood process in soot formation. A generally accepted description of the phenomena is that the BSUs would polymerise and grow in size to formulate liquid-like polycyclic aromatic hydrocarbons (PAHs) precursor species with no internal structure, and subsequently transformed into solid-phase soot nuclei [2,35,36,37,38,39], under certain flaming condition.

Numerically, the rate of nucleation is often expressed as a function of the concentration of acetylene, the main BSUs of hydrocarbon flame [40]. The following nucleation kinetics, while employing Moss [20], Leung [40], and Fairweather [27] approach to yield the reaction rates, refer to Equations (15)–(17), respectively, are considered in this study:

$$J_{Moss}=6\times {10}^6{\overline{\rho }}^2{N}_A\sqrt{T}\exp \left(\frac{-46100}{T}\right)X_{C_2{H}_2}$$

(15)

$$J_{Leung}=\frac{2}{C_m}{N}_{A}1.7\exp \left(\frac{-7548}{T}\right)C_{C_2{H}_2}$$

(16)

$$J_{Fairweather}=\frac{2.7\times {10}^6}{{\rho }_s{C}_m}{N}_A\exp \left(\frac{-46100}{RT}\right)C_{C_2{H}_2}$$

(17)

where $\overline{\rho }$ and ${\rho }_s$ are local mixture density and the soot density ($1800{\mathrm{kg}\text{·}\mathrm{m}}^{-3}$), $N_A$ is the Avogadro number, $R$ is the gas constant, $T$ is the local temperature of the mixture, and $C_{C_2{H}_2}$ and $X_{C_2{H}_2}$ are mole concentration and mole fraction of acetylene, respectively. The $C_m$ term is defined as the minimum number of carbon atoms in the soot nucleus, which is proportional to the pre-defined nucleus size. As suggested in the literature, the value of $C_m$ is defined as 700 in the Leung approach representing soot nucleus with size of about 2.4 nm, and $2.9\times {10}^5$ in the Fairweather approach corresponding to soot nucleus with size of about 18 nm [27,40].

The source term of $k$th moments due to particle inception can subsequently be evaluated as:

$${\overline{S}}_k^{pi}\approx \frac{L_{\epsilon }^k}{k+1}\overline{J}$$

(18)

where $L_{\epsilon }$ is the size of nuclei that is set based on the definition of $C_m$ value.

2.4.2. Aggregation

As nascent soot particulates proceed to flow downstream, they evolve to form larger soot cluster with fractal structure, due to particle-particle collision that has insufficient time to complete particle merging (i.e., coalescence), and such a mechanism is referred to as aggregation [1]. The rate of aggregation or collision efficiency of two soot particulates with collision radius R_c1 and R_c2 can be expressed with a single expression that interpolates the coagulation process between free molecular, continuum, as well as the transition regimes, as [41]:

$$\begin{array}{ll}\beta \left(R_{c1},R_{c2}\right)& =4\pi \left(D_1+D_2\right)\\ & \times \left(R_{c1}+R_{c2}\right){\left[\frac{\left(R_{c1}+R_{c2}\right)}{R_{c1}+R_{c2}+\sqrt{g_1^2+g_2^2}}+\frac{4\left(D_1+D_2\right)}{\sqrt{c_1^2+c_2^2}+R_{c1}+R_{c2}}\right]}^{-1}\end{array}$$

(19)

where

$$c_i=\sqrt{\frac{8k_b\tilde{T}}{\pi m_i}}$$

(20)

$$D_i=\frac{k_b\tilde{T}}{6\pi \mu R_{ci}}\left(\frac{5+4K{n}_i+6K{n}_i^2+18K{n}_i^3}{5-k{n}_i+\left(8+\pi \right)K{n}_i^2}\right)$$

(21)

$$l_i=\frac{8D_i}{\pi c_i}$$

(22)

$$g_i=\frac{{\left(2R_{ci}+l_i\right)}^{\frac{1}{3}}-{\left(4R_{ci}^2+l_i^2\right)}^{\frac{3}{2}}}{6R_{ci}{l}_i}-2R_{ci}$$

(23)

With $k_b$ the Boltzmann constant, $Kn$ the Knudsen number defined as the ratio between the molecular mean free path of gas mixture molecules and mean particle radius, $m$ the soot particle mass, and $R_c$ the collision radius. It is worth noting that, rather than adopting the diameter of soot particles, the collision radius $R_c$ is a characteristic value that resembles the radius of gyration, and it is used to represent the characteristic dimension used to determine collision efficiency. The collision radius $R_c$ is evaluated based on [26]:

$$R_c=\frac{D_0}{2}{\left(\frac{V}{V_0}\right)}^{\frac{1}{D_f}}$$

(24)

where $D_0$ and $V_0$ are diameter and volume of a representative soot primary particle, $V$ is the volume of the soot aggregate, and $D_f$ is the fractal dimension, a morphological property that describes the compactness of a fractal object such as soot particulates.

The source term of $k$th moments due to aggregation can subsequently be formulated as:

$${\overline{S}}_k^{Ag}\approx \frac{1}{2}{\displaystyle \sum }_{\alpha =1}^N{\displaystyle \sum }_{\gamma =1}^N{\left(L_{\alpha }^3+L_{\gamma }^3\right)}^{\frac{k}{3}}{\beta }_{\alpha \gamma }{\omega }_{\alpha }{\omega }_{\gamma }-{\displaystyle \sum }_{\alpha =1}^N{\displaystyle \sum }_{\gamma =1}^N{L}_{\alpha }^k{\beta }_{\alpha \gamma }{\omega }_{\alpha }{\omega }_{\gamma }$$

(25)

2.4.3. Molecular Growth

The mass and size of soot particulates can be increased via chemical and physical surface reactions and processes, for example, the hydrogen abstraction carbon addition based surface growth reaction and the condensation of precursor species onto the particle surfaces. Additionally, they could be decreased as the soot is consumed by the oxidation reactions with ${\mathrm{O}}_2$ molecules and $\mathrm{OH}$ radicals attacking soot clusters as they enter the soot-oxidation zone [2].

In practical soot modelling, as an alternative to the abovementioned complex surface reactions, a simplified reaction model, based on the use of BSUs, i.e., acetylene, to link the gas-phase chemistry to the continuous size change [42], is used to formulate the reaction rate of surface growth $G^{sg}$, as:

$$G^{sg}=\frac{6}{D_f{\rho }_s}{\left(\frac{R_c}{R_{c0}}\right)}^{\frac{3-D_f}{3}}2M_s\times 6\exp \left(\frac{-6038}{T}\right)C_{C_2{H}_2}$$

(26)

where $R_c$ is the collision radius, the subscript “$0$” indicates the primary particle, thus the collision radius of a monomer will be equivalent to its radius, and $M_s$ is the soot molecular weight (12 g·mol⁻¹). It should be noted that the term $\frac{6}{D_f{\rho }_s}{\left(\frac{R_c}{R_{c0}}\right)}^{\frac{3-D_f}{3}}$ is referred to as the soot surface area coefficient that is dependent on the size, as well as the fractal property of the soot aggregates [26,40].

The oxidation process could also be simulated based on the concentration of reacting species and ambient conditions. Two oxidation approaches that have been widely accepted in a commercial package and research code, namely the Said approach [43] and modified NSC approach [20,27], have been investigated in this work. The reaction rate determined based on the two approaches can be expressed as:

$$G_{Said}^{ox}=-\frac{p}{D_f{\rho }_s}{T}^{-\frac{1}{2}}6.5\exp \left(\frac{-26500}{T}\right)Y_{O_2}$$

(27)

and

$$G_{NSC}^{ox}=\left(r{M}_s\frac{6}{D_f{\rho }_s}{\left(\frac{R_c}{R_{c0}}\right)}^{\frac{3-D_f}{3}}\right)$$

(28)

$$r=120\left[\frac{k_A{X}_{O_2}\chi }{1+k_Z{X}_{O_2}}+k_B{X}_{O_2}\left(1-\chi \right)\right]f_{O_2}$$

(29)

$$\chi ={\left(1+\frac{k_T}{k_B{X}_{O_2}}\right)}^{-1}$$

(30)

where $Y_{O_2}$ and $X_{O_2}$ are the mass and mole fraction of oxygen and $k_A, k_B, k_T$ and $k_Z$ are rate constants adopted from semi-empirical investigation [20,40]. Herein, a temperature dependent correction factor $f_{O_2}$, introduced by Liu et al., which tuned to best match the experimental observation of soot formation in smoking flame [40,44,45], is defined as:

$$f_{O_2}={\left(1+\exp \left(-\frac{T-1650}{80}\right)\right)}^{-1}$$

(31)

It should be noted that Liu’s correction factor becomes negligible as the local temperature is below $1300\mathrm{K}$ and it approaches unity when the temperature increases to $2000\mathrm{K}$. Therefore, the correction factor also serves as a cut-off function that limits the oxidation reaction to regions with relatively high temperature.

The source terms of $k$th moments, due to the combined effects of surface growth and oxidation, applying the quadrature approximation, are:

$${\overline{S}}_k^{sg+ox}\approx k{\displaystyle \sum }_{\alpha =1}^N{\omega }_{\alpha }{L}_{\alpha }^{k-1}\left(G^{sg}+G^{ox}\right)$$

(32)

2.5. Numerical and Case Configurations

The geometry of the model that is used in the current study replicates a laboratory scaled co-flow diffusion flame generated from an in-house built co-flow slot burner with a centre slot for fuel injection and two side slots for oxidiser supply. Both ethylene and air are injected at $0.07{\mathrm{ms}}^{-1}$, normal to the burner front through the fuel and oxidiser slots, respectively, and the combustion occurs at atmospheric pressure. It should be noted that, the injection rate of fuel applied in this work is significantly lower than the $50{\mathrm{ms}}^{-1}$ that is used in other numerical studies adopting the DQMOM algorithm [25,26,28]. A computational domain symmetrical with respect to burner centreline, with dimensions of $0.8\mathrm{m}$ and $1.2\mathrm{m}$, is defined and applied with a non-uniform mesh of about $800,000$ elements. Higher mesh resolution is applied at regions near the axis and close to burner front where both the flow development and intense chemical reaction are expected to take place. The dimension of the domain is set based on the observation of flame shape and the mesh configuration has been verified to ensure that the yield simulation result is independent of grid resolution.

One of the potential challenges of the DQMOM approach is that the system could be ill-conditioned when the weight and abscissas of the defined internal coordinate become null, i.e., raised when considering the regions with no soot particles or particulates with infinitely small particle size, which leads to a singular or near-singular condition when solving the linear system (Equations (9)–(12)). This ill-conditioned system could result in either the linear system being unable solve or solved with large numerical error. The issue can be resolved by seeding the domain with particles that have a size distribution that is predefined based upon the PD algorithm that was proposed by Gordon et al. [46]. For mono-variate case with $N=2$, the pre-defined particle size $L_{\alpha }$ can be set as $L_1=0.2113\times L_{\epsilon }$ and $L_2=0.7887\times L_{\epsilon }$ and then applied to regions where the weights are null or abscissas are not defined.

The DQMOM coefficient matrix could also become near singular when the values of any pair of abscissas are equal, or getting too close to each other, resulting in numerical instability and yielding significant numerical error. To tackle the issue, a small perturbation was introduced to the abscissas to ensure non-distinct values of $L_{\alpha }$ in the DQMOM coefficient matrix, as suggested by Marchisio and Fox [24].

Soot radiative heat transfer is taken into consideration by introducing an additional source term into the energy transport equation, which is defined as:

$$Q_{rad}=-\sigma a_s\left(T^4-T_{ref}^4\right)$$

(33)

where $\sigma$ is the Stefan–Boltzmann constant, $T_{\mathit{ref}}$ is the ambient reference temperature that is defined as $300\mathrm{K}$, and the $a_s$ is the absorption coefficient expressed as:

$$a_s=2370T{f}_v$$

(34)

where $f_v$ is the soot volume fraction, as evaluated from the third moment of the NDF.

2.6. Sensitivity Analysis

A comprehensive assessment is carried out to investigate the effect of the key components and submodules that are integrated within the framework on both the simulation result as well as the numerical stability of the system. The components and submodules of interest in this work are defined as: nucleation law, oxidation law, value of fractal dimension, and diffusivity. Nucleation and oxidation are fundamental soot mechanisms through the soot evolution. Fractal dimension is directly correlated to the evaluation of rate of soot mechanisms involving surface reactions. Diffusivity governs the spatial transportation and distribution of the soot species.

A summary of all simulation cases under investigation is tabulated in

Table 1. Summary of cases in the parametric study performed in this work.

Case	Nucleation	Oxidation	Fractal dimension	Schmidt Number
A	Moss	modified NSC	2.0	0.7
B	Leung	modified NSC	2.0	0.7
C	Fairweather	modified NSC	2.0	0.7
D	Moss	Said	2.0	0.7
E	Moss	modified NSC	1.8	0.7
F	Moss	modified NSC	2.2	0.7
G	Moss	modified NSC	2.0	0.5
H	Moss	modified NSC	2.0	0.9

. Please note that, the simulation result generated from the baseline case (Case A) is compared with that evaluated from other cases (Case B to Case H), which have only one change in the modelling aspect made based on Case A, as well as the experimental data obtained from an in-house measurement. As mentioned earlier, the experiment measures the flame temperature and examines the soot particles using thermophoresis soot sampling and TEM imaging, at three selective height above burner (HAB). The selected HABs represent different typical stages of soot evolution within the flame, namely soot inception region ($20\mathrm{mm}$ HAB), soot growth region ($30\mathrm{mm}$ HAB), and soot mature region ($40\mathrm{mm}$ HAB), and they are referred to as ‘exp_L’, ‘exp_M’, and ‘exp_H’ in the following sections, respectively. Morphological information of the soot sample, such as size of primary particles, projected area, radius of gyration, and fractal dimension of the agglomerates, as well as soot volume fraction, are subsequently derived via either direct characterisation of the TEM images or an automatic image processing algorithm. Details of the experimental setup and post image processing are provided in our previous work [47,48,49].

3. Results

In this section, the details of the flame that is generated via flamelet approach bundled with detailed chemistry is firstly presented and validated with experimental measurement. The DQMOM soot models, i.e., Case A to Case H are subsequently coupled with the combustion modules. The result of soot related properties that are generated by the proposed models are presented and assessed to evaluate the effect of the various modelling aspects on the prediction of in-flame soot evolution.

3.1. Flamelet and Combustion Process

Figure 2 and Figure 3 report typical examples of contour plots of flame temperature and reacting species concentration. These results are generated from the baseline case (Case A), with the standard $k$-$\epsilon$ turbulence model, flamelet approach for combustion, Moss law and modified NSC law for nucleation and oxidation rate evaluation, fractal dimension of $2.0$, and turbulent Schmidt number of $0.7$.

As mentioned in the previous section, in this co-flow diffusion flame configuration, the parent fuel (ethylene) is injected from the centre inlet into the computational domain and it reacted with the oxidiser (air) entrained from co-flow inlet and the surroundings. Such a rapid combustion process generates a typical diffusion combustion temperature profile, as presented in Figure 2, which is in good agreement with the thermocouple measurements of flame temperatures at several lateral positions and repeats at three heights above the burner. For instance, the peak flame temperature at various HABs, which were obtained numerically and experimentally, both occur at the intersection where fuel and oxidiser encounter and react. The discrepancy in the comparison may be attributed to the inevitable experimental uncertainties, e.g., the spatial resolution of the measurement, averaging of temperature fluctuation that is caused by flame flickering, etc.

The concentration of reacting species, as generated via detailed chemistry coupled with the flamelet approach, is presented in Figure 3. The parent fuel ethylene (C₂H₄) generally resides at low HAB i.e., regions close to burner front. It subsequently de-composed and reformed into other intermediate species via the pyrolysis process as it flows to the downstream. Acetylene (C₂H₂), the BSUs for soot particle nucleation and surface growth reaction, is capsulated within the flame sheet and its concentration increases with the heights at the low-intermediate flame region before approaching the high temperature oxidation zone. The concentration of oxidant (O₂) peaks at regions that are away from the flame reaction zone and gradually decreases as the lateral position approaches the flame centreline. Unfortunately, due to the lack of available experimental resources, the validation of the concentration of reacting species is not carried out. However, the result presented above is in line with the general expectation of the spatial distribution of flammable gas mixture, as well as oxidant species of a co-flow diffusion flame configuration.

The DQMOM based soot model was subsequently coupled with the established flame profile to characterise the evolution of soot particulates. An initial assessment on the proposed model has demonstrated significant improvement in predicting soot related properties, when compared with that of a mono-dispersed PSD based Moss–Brookes soot model. A comprehensive assessment on the selective modelling aspects, i.e., nucleation law, oxidation law, value of fractal dimension, and diffusivity, is then performed to further enhance the numerical robustness and simulation accuracy of the proposed model.

3.2. Sensitivity Analysis of DQMOM Key Modelling Parameters

Herein, four key modelling aspects were integrated into the DQMOM framework, namely nucleation law, oxidation law, value of fractal dimension, and diffusivity, are taken into consideration in this sensitivity analysis.

3.2.1. Nucleation

Nucleation involves the transformation of soot particles from gas-phase BSUs to amorphous intermediate species (e.g., PAH) and finally to solid phase nuclei, and, to date, it is the least understood mechanism in the soot process. From the numerical modelling aspect, the rate of nucleation dominates the overall quantity of soot been formulated, such as number density and concentration. In the framework of the DQMOM model, the rate of nucleation governs the soot number density (i.e., weight). This weight of the internal coordinate is a key parameter in the formulation of integrated source term of $k$th moments describing all major soot mechanisms, such as nucleation, coagulation, surface growth, and oxidation (as stated in Equations (18), (25) and (32)). Therefore, the nucleation law that was employed in DQMOM framework could directly affect the evaluation of all weights and weighted abscissas, as well as the internal coordinate (i.e., particle size). It is the most significant aspect and it has been investigated firstly in this sensitivity analysis. Three nucleation laws namely Moss law (Case A), Leung law with nuclei size of $2.4\mathrm{nm}$ (Case B), and Fairweather law with nuclei size of $18\mathrm{nm}$ (Case C), have been adopted for assessment. The results of soot-related properties are presented in Figure 4, Figure 5, Figure 6 and Figure 7.

Figure 4 presents the rate of nucleation that is derived from the investigated nucleation laws. When compared with the Leung law and Fairweather law, the region nucleation dominants given by Moss law is more confined and it is located further downstream where the combustion process is well developed with higher flame temperature. Studies on the activation temperature of the nucleation process suggested that intensive soot formation initiates when the local temperature roughly exceeds $1250\mathrm{K}$ and terminates where the concentration level of precursor becomes low [50]. The nucleation rate profile that was generated by Leung and Fairweather appears to be less realistic when compared with the Moss law, in the implementation of the current burner configuration. Both of the approaches tend to overpredict the nucleation process at regions with high concentration of acetylene but they have relatively lower flame temperature, i.e., closer to burner front or away from reaction zone on the side towards flame centreline. The difference in spatial distribution of the nucleation dominant regions may be attributed to the deviation of the empirical constant that is used in the rate formulations, such as frequency factor and activation energy.

The soot particle number density that is reported in Figure 5 presents a profile similar to that of the nucleation rate and confirms the speculation. It is noticeable that the profile of number density of model adopting Moss law is constrained to the high temperature region. On the other hand, the profile that was predicted by the Leung Law and Fairweather law shows that a large number of soot particulates accumulate at low-intermediate flame height as well as at regions closer to the flame centreline, where local temperature is relatively low, which differs from the findings in the experimental study of similar flame configuration [8,51,52].

In the framework of the DQMOM soot model, the correlation between the nucleation rate and particle size is intricate. The size of soot particulate, calculated as a ratio between the abscissas and weight (i.e., number density), is inversely proportional to the number of particulates in the control volume. On the other hand, the integrated source term of surface reaction, such as surface growth and oxidation, is expressed as a function of number density (as stated in Equation (32)). It implies that the occurrence frequency of collision events i.e., BSU’s adhere to or oxidation radicals attack the agglomerates, should be proportionate to the concentration of the soot particulates within the control volume. Therefore, the effect of the nucleation rate on particle size is numerically ambiguous.

Nonetheless, the particle size predicted by models utilizing the three nucleation law is presented in Figure 6. From the contour plot, it can be seen that the soot particle size determined by the Fairweather law is relatively lower, whilst that provided by the Moss law and Leung law are more comparable, when compared with the experimental data at the sampling points. However, at the flame centreline, the particle size that is predicted by Moss law approach at three HABs is more scattered when compared with that estimated using the Leung law and Fairweather law. It captures the feature that the size of soot particulates increases with flame heights before approaching flame tip. Such a tendency could also be found in previous studies on in-flame soot formation [8,34].

The simulation result of soot volume fraction, when compared with experimental data, is presented in Figure 7. It is evident that the Fairweather Law severely underestimates the soot volume fraction due to the combined effects of underpredicting soot number density as well as soot particle size. The value of soot volume fraction from Moss law and Leung law is again comparable with experimental data at the sampling locations. However, a more distinguished difference in soot volume fraction between the flame centreline and the expected soot zone can be observed from the Moss law when compared with that predicted by the Leung law. Such a profile indicates that most soot species are confined in a narrow sooting region and its concentration sharply decreases as the lateral position moves away in both directions: towards flame centreline and surrounding. Again, this observation agrees well with the experimental data [8,34,53] as well as the general expectation. Therefore, based on the abovementioned analysis and comparison, the Moss law is selected as the preferred mechanism of nucleation and it is adopted in the framework of DQMOM soot model.

3.2.2. Oxidation

The effect of implementing a different law of oxidation, which functioned as the destructive mechanism of soot particle size, is assessed in this study. The analysis herein consists of the comparison of models employing a modified NSC law (Case A) and Said law (Case D), and the results are reported in from Figure 8, Figure 9, Figure 10 and Figure 11.

The oxidation rate estimated from both cases presented in Figure 8, shows a notable difference regarding the oxidation dominant region. It is generally accepted that the rate of oxidation in diffusion flame is expected to be predominant at regions where flame temperature and oxygen concentration become high i.e., alongside with the reaction zone. Despite the difference in the magnitude of oxidation rate, the oxidation profile that is provided by the Said law predominantly occurs at the downstream flame region and close to flame front where a peak temperature and high concentration of oxygen persist. This could be justified as the oxidation rate in the Said law is formulated based on the Arrhenius equation as a function of local temperature and concentration of oxidant species, as referred to in Equation (27). Alternatively, the oxidation rate evaluated from modified NSC law peak at low-intermediate flame height in the vertical direction and constrained in a narrow band close to soot layer in the lateral direction. This could be well explained as the effect of the additional coefficient of surface area that is introduced in the modified NSC law that limits the reaction region of oxidation to where larger soot particulates reside.

As mentioned previously, as a surface reaction, oxidation occurs via oxidant species attack the soot agglomerate on its active surface site. Apart from the common aspects, such as ambient condition and activation energy, the oxidation rate should be proportionate to the surface area of the soot agglomerate where the reaction is taking place. To be more specific, under the same ambient environment, e.g., with the same oxidant species concentration and local temperature, the control volume with the presence of larger particulates is expected to have a higher oxidation rate, when compared with that of smaller particulates. The Said law neglects such a correlation and hence its evaluated oxidation rate peaks near the flame tip where the soot particulates could merely be observed due to the intensive oxidation reaction. However, it underestimates the oxidation rate at a high temperature reaction region at low-intermediate HABs where larger soot particulates are formed. On the other hand, modified NSC law takes the factor of the active surface site into consideration and enhances the oxidation intensity where larger soot particulates persist, i.e., in the vicinity of soot layer. This deduction is verified by comparing the oxidation rate at low-intermediate flame height, as shown in Figure 9. It can be observed that, when compared with the modified NSC law, the Said law overlooks the effect of oxidation from flame centreline to soot and reaction zone and beyond in the lateral direction, at all three heights of interest.

A similar conclusion could be drawn when comparing the soot particle sizes that were obtained from the two oxidation kinetics, as presented in Figure 10. The particle size evaluated from the two approaches have a nearly identical profile in spatial distribution and comparable magnitude at the downstream flame regions. However, a larger size particulate has been evaluated by Said law at upstream flame heights, i.e., about $20\mathrm{mm}$ HAB in the vertical direction and from flame centreline to and beyond the reaction zone in the lateral direction, when compared with that of modified NSC law. This inaccurate description of oxidation phenomenon leads to a distinct thin layer of un-oxidized soot particles at the lateral location near and beyond the flame reaction zone in the profile of soot volume fraction generated by Said law, as presented in Figure 11. Such an un-oxidated soot layer is contrary to the findings of both the numerical and experimental studies of similar flame configuration [8,25,26,27,34,40]. Based on the abovementioned comparison in reaction rate, soot particle size, as well as the soot volume fraction, it is believed that the modified NSC law is more applicable for the flame of interest, and therefore it is selected as the preferred oxidation law to be implemented in the DQMOM model.

3.2.3. Fractal Dimension

The fractal dimension, $D_f$, gives insights into the geometrical characteristics of fractal object, e.g., soot clusters. Soot aggregates with higher degree of compactness are generally assumed to have larger value of $D_f$, whereas aggregates with more linear and open branched structures are characterized by a smaller $D_f$ value [54]. In the numerical simulation of soot formation, the fractal dimension is closely related to the prediction of rate of processes that involve radical species attack on reactive site and event of collision, for example, surface reactions and coagulation. Fractal dimension is a key factor in quantifying size-related properties, e.g., radius of gyration, collision radius, effective diameter, etc. For computational simplicity, a unique value of $D_f$ of soot aggregates is generally accepted and assumed in mono-variate population balance model, with a typical value that ranges between $1.8$ and $3.0$. Data from previous studies on a non-premixed diffusion flame have demonstrated that the selection of fractal dimension could have a marked impact on the final prediction [26,28]. Therefore, it is important to investigate $D_f$ and evaluate its influence on the current burner configuration.

The fractal dimension of a fractal object could be quantified based on the morphological information that was obtained experimentally by examining the TEM image and is derived by least-square regression fittings to the power law relationship [55,56]:

$$n=k_f{\left(\frac{r_g}{d_p}\right)}^{D_f}$$

(35)

whereby $k_f$ is the empirical fractal prefactor and $d_p$ is mean diameter of primary particle within the aggregate. $n$ is the number of primary particles per aggregate, determined based on the ratio of the projected area of the aggregate $A_a$ and average projection area of the primary particles $A_p$ observed within the aggregate, as well as the empirical constant $k_a$ and exponent $\alpha$ set as $1$ and $1.09$, respectively:

$$n=k_a{\left(\frac{A_a}{A_p}\right)}^{\alpha }$$

(36)

$r_g$ refers to the radius of gyration of the soot aggregates, and it is determined using:

$$r_g^2=\frac{1}{m}{{\displaystyle \sum }}^{\text{}}{r}_i^2$$

(37)

where $r_i$ is the distance from the centroid of the aggregate to individual pixels and $m$ is the number of pixels within the projection image of an aggregate. The slope of the least-square linear fits of $\log \left(n\right)$ and $\log \left(\frac{r_g}{d_p}\right)$ were subsequently performed to yield the value of fractal dimension.

Examining of the TEM images of soot samples that were collected at three HABs suggests that the fractal dimension of the soot clusters that were generated from the flame of interest is within the range of $1.8$ and $2.0$, which agrees well with the results that were reported in the literature [57]. Based on this observation, three values of $D_f$, i.e., $1.8$ (Case E), $2.0$ (Case A), and $2.2$ (Case F), have been adopted in the DQMOM soot models for assessment, and the result is reported in Figure 12, Figure 13, Figure 14 and Figure 15.

As mentioned previously, within the DQMOM soot model framework, the fractal dimension primarily affects the determination of the rate of surface reaction mechanisms. For fractal objects with equal size, i.e., the longest distance could be made from the TEM image between any two pixels within a soot aggregate, the particulates with larger $D_f$ tend to be more morphologically compact. Its effective diameter, quantified as collision radius (refers to Equation (24)), therefore will be smaller when compared with that of a particle with smaller value of $D_f$. This implies that aggregate with larger value of $D_f$ could have limited active site for surface reaction to take place, and as a result restrain the rate of surface reaction, such as oxidation and surface growth, when compared with that of smaller value of $D_f$. The rate of oxidation and surface growth, therefore, are compared and presented in Figure 12 and Figure 13, respectively.

The oxidation rates presented in Figure 12 generated by models with three $D_f$ are almost identical, regarding the spatial distribution of reaction region as well as the magnitude of the reaction rate. This could be attributed to the temperature dependent correction factor $f_{O_2}$ in the modified NSC law. This correction factor restrains the oxidation to occur only at where the local temperature exceeds $1300\mathrm{K}$ and cut off the oxidation rate to null when the local temperature is below the threshold. This cut-off correction regulates the profile and the intensity of oxidation based on the local temperature and concentration of oxidant species, hence it dampens the effect of the fractal dimension to the evaluation of oxidation rate.

On the other hand, the rate of surface growth of the model with smaller fractal dimension ($D_f=1.8$) is about nine-fold larger than that yield from the model with larger $D_f$ value ($D_f=2.2$) as shown in Figure 13. The prominent increase in surface growth rate with the decrease of $D_f$ value could be explained as the surface growth rate is formulated based on the local temperature, concentration of BSUs, as well as the coefficient of surface area, as stated in Equation (26). With the previously mentioned hypothesis, the rate of surface growth is proportionate to the surface area coefficient, $\frac{6}{D_f{\rho }_s}{\left(\frac{R_c}{R_{c0}}\right)}^{\frac{3-D_f}{3}}$. The decrease of $D_f$ value will exponentially increase the surface area coefficient, and hence intensify the ensuing rate of surface growth. Similarly, the collision radius, $R_c$, is proportional to the number of soot primary particles within the aggregate to the power of $\frac{1}{D_f}$ (refer to Equation (24)). With the implementation of a smaller $D_f$ value, the collision radius will also exponentially increase and subsequently result in an exaggerated prediction of the surface growth rate.

The particle size contour and comparison between experimental data and simulation results at three measuring points are presented in Figure 14. The model adopting fractal dimension of 1.8 evaluates the particle with a size significantly larger due to the excessively intense surface growth reaction, i.e., increased about three times when compared with that of fractal dimension of 2.2. This observation is consistent with the sensitivity analysis reported in the literature that compares the particle size at flame centreline at various HAB, generated from models with a fractal dimension ranging from 1.8 to 3.0 [26]. The ensuing estimation of the soot volume fraction based on the number density and particulate size, as presented in Figure 15, again agrees well with the abovementioned speculation. The soot volume fraction evaluated with the adoption of a fractal dimension of 1.8 is significantly larger, i.e., by four orders of magnitude at flame centreline, when compared with that of model implementing the value of 2.2.

From the comparison of soot particle size and soot volume fraction with the experiment data, it could be concluded that the $D_f$ value of 2.0 provides the best agreement, whilst the $D_f$ value of 1.8 and 2.2 either overestimates and underestimates the two soot related properties. This finding is supported by the literature [26,28], and therefore value of 2.0 is selected as the preferred settings for fractal dimension and it is integrated into the proposed DQMOM soot model framework.

3.2.4. Diffusivity

The diffusion term, accounts for the transport of scalar due to concentration gradients, is also investigated in this study. In dilute species transport, the flux due to diffusion is given by Fick’s first law. Its intensity is dependent upon the steepness of the concentration gradient and it is numerically proportional to a coefficient that is referred to as diffusivity. This diffusivity could be evaluated as the summation of molecular diffusivity, $\Gamma$, and the turbulent diffusivity, ${\Gamma }_t$. For the current study of soot formation in the co-flow diffusion flame, the molecular diffusion can be safely neglected as the turbulent diffusivity is dominant when compared with molecular one [25,26]. The turbulent diffusivity is calculated as the ratio between eddy viscosity, ${\nu }_t$, and turbulent Schmidt number, $S{c}_t$, which describes the relative diffusion of momentum and mass due to turbulence. Therefore, the rate of mass diffusion is inversely proportional to the turbulent Schmidt number. The turbulent Schmidt number is considered as an empirical constant in the most practical fire simulation. For computational simplicity, a global value of Schmidt number i.e., $S{c}_t=0.7$ is generally prescribed as a widely accepted value in the modelling of jet flows and gas turbine combustor modelling [58]. However, other values of turbulent Schmidt number, ranging from $0.2$ to $1.0$, have been used in a number of applications, e.g., the modelling of gas turbine combustor, combustor with rich burn, quick quench, and lean burn, etc. [21].

The soot evolution and the spatial transport concurrent with the flammable gases in the current study have been resolved via a mixture model approach. To be more specific, the particulates that are immersed in the gas-phased mixture are assumed to be infinitely small; hence they have negligible effect on the global flow field. A typical value of $S{c}_t=0.7$ was initially implemented when coding the transport equation for both soot species and gas mixtures. However, soot particulates exhibited a complexed phase evolution that depends on the stage in the formation pathway. i.e., from liquid-alike nascent soot precursor or PAH with amorphous internal structure, to graphitized mature soot particulate. The tendency in spatial transport due to diffusion of the soot particulate should be different when compared with gas-phased mixture. The assessment of the effect of the turbulent Schmidt number on the prediction of soot properties of the current flame configuration is therefore carried out.

Three values of $S{c}_t$, $0.5$ (Case G), $0.7$ (Case A), and $0.9$ (Case H) have been selected in this assessment. The contour plot of soot number density, as evaluated by the three models, is firstly presented in Figure 16, which are almost identical in terms of both spatial distribution as well as peak value. It is generally agreed that the effect of diffusion term on the prediction result differs from the previously exanimated soot kinetics. For example, the variation of nucleation law, oxidation law, and the value of fractal dimension vary the formulation of integrated source term and affect the ensuing prediction of the sooting processes. However, the diffusion term in transport equation does not have a direct impact on the rate of birth and death process of soot evaluation. As a result, non-distinct discrepancy in terms of peak value of number density of the three cases can be observed.

Nevertheless, as depicted in Figure 17, the quantity of soot number density evaluated by the three cases significantly varies at regions away from the soot layer, and such a trend can be observed in all selected HABs”. The result clearly demonstrates that the dispersion of soot particulates is notably affected by the diffusivity, as it governs the spatial transport of soot species from their inception region to the surrounding. For instance, with the smaller Schmidt number, i.e., $S{c}_t=0.5$, there are more soot particles being spread away from the nucleating region towards non-active regions, such as flame centreline as well as regions away from flamelet. Whilst with a lager Schmidt number, i.e., $S{c}_t=0.9$, most of the soot particulates are confined in a narrow band with a relatively low quantity of soot particles being observed elsewhere. The result implies that denser solid soot particles with higher $S{c}_t$ values are inherently more resistant to the spatial movement due to diffusion, when compared with that of a lower $S{c}_t$ value.

The soot particle size and the soot volume fraction are subsequently derived as the first and third moment of the number density function, and as a result, a similar trend of the two quantities is observed. For instance, the contour plot of soot particle size and soot volume fraction shown in Figure 18 and Figure 19 depict a relatively similar profile in terms of the quantity as well as the general spatial distribution. However, a close observation clearly indicates that, with the increase of $S{c}_t$ number from $0.5$ to $0.9$, the profile of both soot particle size and soot volume fraction at low-intermediate HAB transformed from an oval-shaped distribution that peaks at the centreline to an A-shaped profile that confines the two quantities to regions close to the expected soot band. Similarly, the X-Y plot of soot particle size and soot volume fraction at the selected HABs again indicates that both quantities at all the heights of interest shift from a typical decay curve peaks at the flame centreline, to a positively skewed bell-shaped curve, with the increase of the $S{c}_t$ value.

From the abovementioned comparison, the conclusion could be made that, with the implementation of $S{c}_t$ of $0.9$, the simulation results of both soot particle size as well as soot volume fraction agree best with in-house experimental data at all sampling heights. The result evaluated with larger $S{c}_t$ also resembles the planar measurement from the LII techniques. This observation is again aligned with the expectation that denser soot particulate should be inherently more resistant toward spatial transport due to diffusion, as compared with other gas-phase mixtures. Therefore, the $S{c}_t$ of $0.9$ is selected as the preferred value for diffusivity and it is integrated into the proposed DQMOM soot model framework.

4. Conclusions

In this article, a fire field model that is associated with an in-house Direct Quadrature Method of Moments (DQMOM) code was specifically developed for buoyancy-driven low-speed turbulent diffusion flames. The DQMOM model embraced all major mechanisms for soot particle evolution/disruption, including nucleation, surface growth, and oxidation and coagulation. Furthermore, it was fully coupled with detailed chemistry gas-phase combustion model, which allows for predictions of intermediate chemical species to enhance the real-time tracking of the particle size distribution. The proposed model was not only capable of providing additional information for the soot particle distribution in the computational field, but it also yielded other important secondary properties, i.e., particle number density and local volume fraction with good accuracy and effective computational cost.

Using the proposed framework, a collection of soot nucleation and oxidation kinetics were numerically investigated and then compared against a self-conducted experiment while using ethylene/air burner with soot diagnostics imaging measurements. When compared with the Leung law and Fairweather law, it was found that the Moss nucleation law provided more accurate prediction of nucleation rate and the subsequent particle number density distribution. The Moss law was capable of limiting the soot inception to high flame temperature region. With regards to oxidation kinetics, the modified NSC law demonstrated better results when compared to the Said law. The modified NSC law constrained the oxidation rate at the surrounding region with the highest concentration of soot particulates, while the Said law under-predicted the oxidation rate at low-intermediate flame heights. In addition, the fractal dimension and Schmidt number were also examined in this study, and it was discovered that the value set of 2.0 and 0.9 yields the best result when compared with the experimental data.

In summary, combining the advantage of the improved DQMOM approach with the detailed chemistry combustion model, the proposed in-house CFD model significantly enhanced the temperature field, intermediate chemical species, and soot yield predictions, as well as providing a more comprehensive physical description of the evolution and disruption processes for soot particles. This enabled a better understanding of the nucleating region for soot particles and it agreed well against the experimental measurements. The high-fidelity prediction of soot particles will be a major advancement that is applicable to the design of fire protection engineering systems, including smoke extraction, curtains, and alarms.