3. Soot Models

3.1. Soot Source Terms

Spitfire supports a few methods for computing coefficients of soot growth, oxidation, nucleation, etc. in a flamelet library. These coefficients can be used in the source terms of transport equations for the number density and mass fraction of soot particles. This two-equation approach arises from representing soot as a monodisperse population of spherical particles. The source term for the soot number density (\(N_s\), in kmol soot / kg air) involves nucleation and coagulation while the source term for the soot mass fraction (\(Y_s\)) involves nucleation, condensation, surface growth, and oxidation. These can be generally written for laminar flow as

(3.1)\[\begin{split}\dot{S}_{N_s} &= (\mathrm{molar\ nucleation}) - (\mathrm{coagulation}) N_s^{11/6}Y_s^{1/6} \\ \dot{S}_{Y_s} &= (\mathrm{mass\ nucleation}) + (\mathrm{mass\ condensation})N_s + ((\mathrm{surface\ growth}) - (\mathrm{oxidation}) )N_s^{1/3}Y_s^{2/3}\end{split}\]

When we compute filtered source terms for turbulent flow, we convolve each tabulated coefficient separately, then substitute these into (3.1) along with a filtered \(\tilde{Y}_s\) and \(\tilde{N}_s\). This can be thought of as assuming delta PDFs for \(Y_s\) and \(N_s\) and including them as dimensions in the convolutions with statistical independence from the other table dimensions (mixture fraction, stoichiometric dissipation rate, etc.).

Spitfire supports three models for computing these source term coefficients:

  • Aksit_Moss

  • Aksit_Moss_Benzene

  • PAH_alpha

The coefficients are added to a library, given a chemical mechanism and soot model name, through the function tabulate_soot_source_coefficients. The corresponding source term function can then be grabbed using get_soot_source_term_evaluator and, likewise, a source term Jacobian evaluator can be obtained using get_soot_source_term_jacobian_evaluator.

3.1.1. Aksit-Moss models

The Aksit-Moss models [AksitMoss] use acetylene-based nucleation and surface growth rates with empirically-determined parameters. They can also incorporate aromatic contributions based on benzene concentrations.

[AksitMoss]

J.B. Moss and I.M. Aksit, Modelling soot formation in a laminar diffusion flame burning a surrogate kerosene fuel, Proceedings of the Combustion Institute, Volume 31, 2007

In both Aksit-Moss models, the coagulation coefficient (coag) is given by

(3.2)\[\mathrm{coag} = \sqrt{\frac{24 V T R_u N_a}{\rho_s}}\rho^2\]

where \(\rho_s\) is the density of soot (default \(1800 \ \mathrm{kg}/\mathrm{m}^3\)) and \(V\) is the volume factor

(3.3)\[V = \left(\frac{6}{\pi \rho_s N_a}\right)^{1/3}\]

Thermodynamic instability of PAH, along with oxidation below temperatures of about 2000 K, leads to the phenomenon referred to as the soot bell curve. We use a correction factor \(B(T)\) applied multiplicatively to the nucleation and surface growth rates to account for this. This bell curve multiplier is written as

(3.4)\[\begin{split}B(T) &= \mathrm{max}\left(1 - \exp{\left( \frac{\Delta T_\mathrm{act}}{T_c} - \frac{\Delta T_\mathrm{act}}{T}\right)}, 0\right) \\ \Delta T_\mathrm{act} &= 16000\ \mathrm{K} \\ T_c &= 2050\ \mathrm{K}\end{split}\]

For the Aksit_Moss model without benzene, molar nucleation (nucl) and surface growth (surf) depend solely on the concentration of acetylene (\(\mathrm{C}_2\mathrm{H}_2\)) with a linear relationship between the number of active sites and the soot surface area:

(3.5)\[\begin{split}\mathrm{nucl_{\mathrm{C}_2\mathrm{H}_2}} &= B(T) k_\mathrm{nucl}[\mathrm{C}_2\mathrm{H}_2] \\ k_\mathrm{nucl} &= 54 \exp{\left(\frac{-21100}{T}\right)}\end{split}\]
(3.6)\[\begin{split}\mathrm{surf_{\mathrm{C}_2\mathrm{H}_2}} &= \rho A B(T) k_\mathrm{\mathrm{C}_2\mathrm{H}_2}[\mathrm{C}_2\mathrm{H}_2] \\ k_\mathrm{\mathrm{C}_2\mathrm{H}_2} &= 11700 \exp{\left(\frac{-12000}{T}\right)}\end{split}\]

with an area factor

(3.7)\[A = V^2 \pi N_a\]

The mass of an incipient particle is given directly as \(M_p =144\) kg/kmol, which multiplies the molar \(\mathrm{nucl_{\mathrm{C}_2\mathrm{H}_2}}\) to get the mass nucleation rate.

The Aksit_Moss_Benzene model uses a nucleation coefficient that is twice that of Aksit_Moss in order to account for the contribution of aromatic components. The surface growth in Aksit_Moss_Benzene also incorporates aromatic contributions by depending on both acetylene and benzene (A1). The coefficient for acetylene growth is identical to Aksit_Moss while the coefficient for benzene growth is given by

(3.8)\[\mathrm{surf}_{\mathrm{A}1} = 20 \rho A B(T) k_\mathrm{\mathrm{C}_2\mathrm{H}_2}[\mathrm{A}1]\]

In both Aksit-Moss models, oxidation depends on the concentrations of both oxygen (\(\mathrm{O}_2\)) and the hydroxyl radical (\(\mathrm{OH}\)):

(3.9)\[\begin{split}\mathrm{surf}_{\mathrm{O}_2} &= \rho A k_{\mathrm{O}_2} [\mathrm{O}_2] \\ k_{\mathrm{O}_2} &= 500 \sqrt{T} \exp{\left(\frac{-20000}{T}\right)} \\ \mathrm{surf}_{\mathrm{OH}} &= \rho A k_{\mathrm{OH}} [\mathrm{OH}] \\ k_{\mathrm{OH}} &= 4.2325 \sqrt{T}\end{split}\]

The get_soot_source_term_evaluator function will return the soot_src_aksit_moss function for library coefficients generated with the Aksit_Moss model and will return the soot_src_aksit_moss_benzene function for library coefficients generated with the Aksit_Moss_Benzene model.

3.1.2. Newer PAH-based models

The PAH_alpha model involves the effects of larger aromatics on the production of soot.

Gas phase molecular collisions are modeled for the free molecular regime, with self-collisions between molecules of species i having frequency

(3.10)\[\beta_{i-i} = 2.2 \sqrt{ \frac{4 \pi R_u T }{ M_C N_{C,i}}} d_p^2\]

where \(M_C\) is the molecular weight of carbon, \(N_{C,i}\) is the number of carbon atoms in PAH species i, and \(d_p\) is the average diameter of the PAH particles:

(3.11)\[d_p = \left(6 \frac{2 \overline{MW_{\mathrm{PAH}}}}{N_a \pi \rho_s} \right)^{1/3}\]

using the density of soot \(\rho_s\) and average molecular weight (MW) of the PAH species.

Mass nucleation is given by the sum over PAH species \(i\):

(3.12)\[\mathrm{dimer}\_\mathrm{src} = \sum_i 2 MW_i (\gamma_i \beta_{i-i} [\mathrm{PAH}_i]^2 N_a)\]

for the following species with their sticking coefficients \(\gamma_i\) in parenthesis: naphthalene (0.001), acenaphthalene (0.003), biphenyl (0.0085), phenanthrene (0.015), acephenanthrylene (0.025), pyrene (0.025), flouranthene (0.025), and cyclopentapyrene (0.039). Molar nucleation is then given by \(\frac{0.5}{361.4} \mathrm{dimer}\_\mathrm{src}\). We note that these nucleation rates tend to dominate and over-predict soot production. Therefore, when evaluating the source terms with soot_src_pah_alpha a factor of \(10^{-6}\) is applied to both nucleation rates by default through the argument dimer_fac. This seems like the result of a fundamental inconsistency and is why this model is not supported beyond alpha-stage readiness. Assistance from soot experts would be terrific in productionizing PAH-based models such as this.

Coagulation is the same as in the Aksit-Moss models.

The surface area factor for the surface growth and oxidation rates is

(3.13)\[G_A = \rho N_A V^2\]

A more fundamental approach for modeling surface growth, compared to the Aksit-Moss implementation, employs elementary reactions to describe the surface chemistry, along with a quasi-steady concentration of radical surface sites. Following the formulation from [pahsg], the surface growth rate for our PAH model is

(3.14)\[\begin{split}\mathrm{sg} &= 1.7 \cdot 10^{19} k_4 f^* [\mathrm{C}_2\mathrm{H}_2] \frac{2 M_C}{N_A} G_A \\ f^* &= \frac{ 1 }{ 1 + \frac{k_{1b} [\mathrm{H}_2\mathrm{O}] + k_{2b} [\mathrm{H}_2] + k_{3b} [\mathrm{H}] + k_4 [\mathrm{C}_2\mathrm{H}_2]}{k_{1f} [\mathrm{OH}] + k_{2f} [\mathrm{H}] + k_{3f}} } \\ k_{1f} &= 6.72 \cdot 10^{-2} T^{3.33} \exp\left( \frac{-6.09 \cdot 10^6}{R_u T} \right) \\ k_{1b} &= 6.44 \cdot 10^{-4} T^{3.79} \exp\left( \frac{-27.96 \cdot 10^6}{R_u T} \right) \\ k_{2f} &= 1.08 \cdot 10^{5} T^{1.8} \exp\left( \frac{-68.42 \cdot 10^6}{R_u T} \right) \\ k_{2b} &= 8.68 \cdot 10^{1} T^{2.36} \exp\left( \frac{-25.46 \cdot 10^6}{R_u T} \right) \\ k_{3f} &= 1.13 \cdot 10^{16} T^{-0.06} \exp\left( \frac{-476.05 \cdot 10^6}{R_u T} \right) \\ k_{3b} &= 4.17 \cdot 10^{10} T^{0.15} \\ k^4 &= 2.52 \cdot 10^6 T^{1.1} \exp\left( \frac{-17.13 \cdot 10^6}{R_u T} \right)\end{split}\]

where each of the rates \(k_i\) are in Hz.

[pahsg]

M. Mueller, G. Blanquart, and H. Pitsch, Hybrid method of moments for modeling soot formation and growth, Combustion and Flame, Volume 156, 2009

The oxidation rates depend on the concentrations of O, OH, and \(\mathrm{O}_2\) with linear dependence on soot surface area. The rate for oxygen has a modified Arrhenius form while the rates for OH and O are based on Brownian collision of a radical species with a large soot particle. In the equations below, \(M_{\mathrm{O}_2}\) is the molecular weight of oxygen, \(M_{\mathrm{OH}}\) is the molecular weight of OH, and \(M_\mathrm{O}\) is the molecular weight of O.

(3.15)\[\begin{split}\mathrm{ox}_{\mathrm{O}_2} &= 1.7 \cdot 10^{19} k_{\mathrm{O}_2} \frac{2 M_C}{N_A} [\mathrm{O}_2] G_A \\ k_{\mathrm{O}_2} &= 2.2 \cdot 10^6 \exp \left( \frac{-31.38 \cdot 10^6}{R_u T} \right)\end{split}\]
(3.16)\[\mathrm{ox}_{\mathrm{OH}} = 0.04 M_C \sqrt{\left( \frac{R_u T}{2 \pi M_{OH}} \right)} [\mathrm{OH}] G_A\]
(3.17)\[\mathrm{ox}_{\mathrm{O}} = 0.2 M_C \sqrt{\left( \frac{R_u T}{2 \pi M_O} \right)} [\mathrm{O}] G_A\]

The get_soot_source_term_evaluator function will return the soot_src_pah_alpha function for library coefficients generated with the PAH_alpha model.

3.1.3. Correcting for soot leakage

Soot particles are formed at fuel-rich mixtures and then oxidized as they move toward fuel-lean mixtures. In many scenarios, such as smoking flames, this oxidation happens at a finite rate leading to unoxidized soot leaking into lean mixtures. One manner in which this can happen is when radiative heat losses are large enough to substantially lower oxidation rates through the lower temperatures, such as happens in flames with large residence times and Damköhler numbers. Local extinction events, such as through fast turbulent transport, can also lead to suppressed oxidation. Finite-rate oxidation models are important for modeling these kinds of phenomena.

Recent developments in soot modeling for turbulent flames have introduced the concept of a soot subfilter PDF. This can account for the effect of soot leaking into the oxidizer side of the flame on both the soot growth and oxidation rates, which in turn reduces the tendency of flamelet models to over-predict soot oxidation [Col2023]. The approach starts by finding the mixture fraction value, termed \(Z_\mathrm{soot}\), at which the soot surface growth and soot oxidation rates are equivalent. The subfilter is then defined as a hyperbolic tangent function of mixture fraction, parameterized by \(Z_\mathrm{soot}\) and a soot oxidation Damköhler number , \(\zeta\), which determines the width of the transition

(3.18)\[\zeta = \frac{\chi (Z_\mathrm{soot})}{2 k_\mathrm{ox}(Z_\mathrm{soot})}.\]

In (3.18), \(k_\mathrm{ox}\) is the sum of oxidation coefficients at a given mixture fraction, and \(\chi\) is the scalar dissipation rate function of mixture fraction. In practice, we interpolate \(k_\mathrm{ox}\) over mixture fraction for a given Spitfire library in order to evaluate it at \(Z_\mathrm{soot}\) while the algebraic form of \(\chi(Z)\) is known.

[Col2023] (1,2)

Hernando Maldonado Colmán, Antonio Attili, and Michael E. Mueller, Large Eddy Simulation of turbulent nonpremixed sooting flames: Presumed subfilter PDF model for finite-rate oxidation of soot, Combustion and Flame, Volume 258, 2023

Spitfire builds on the work of Colmán et. al. [Col2023], by including a soot_leakage_fraction as an additional parameter to the soot subfilter \(C_s(Z)\)

(3.19)\[C_s(Z) = 0.5(1+\mathrm{soot\_leakage\_fraction}) + 0.5(1-\mathrm{soot\_leakage\_fraction})\tanh((Z - Z_\mathrm{soot} + \zeta)/\zeta).\]

Running tabulate_soot_source_coefficients with \(0\leq\mathrm{soot}\_\mathrm{leakage}\_\mathrm{fraction}<1\) creates a subfilter that multiplies all surface growth and oxidation coefficients in the library and gets added itself to the library as the property ‘soot_subfilter’. A soot leakage fraction of 1 is equivalent to not applying a soot subfilter, in which case oxidation and growth rates are not modified compared to the original soot model. Decreasing this fraction corresponds to increasing oxidation rates, which further limits the amount of soot that would be available for oxidation at \(Z<Z_\mathrm{soot}\). This leads to smaller tabulated oxidation coefficients for lean mixtures. A fraction of 0 represents infinitely fast oxidation, where no soot is present for \(Z<Z_\mathrm{soot}\) and the subfilter resembles a Heaviside. When building a turbulent library, the corrected surface growth and oxidation coefficients (\(\mathrm{surf}_i C_s\)) are convolved with the mixture fraction PDF and normalized by the convolution of the soot subfilter itself with that same PDF:

(3.20)\[\widetilde{\mathrm{surf}_i} = \frac{\int \mathrm{surf}_i C_s P_z dZ}{\int C_s P_z dZ}\]

3.2. Gas Radiation Models

Absorption and radiation gas coefficients can be added onto a library through the function tabulate_gas_radiation_model. Spitfire supports a ‘Barlow_GreyGas’ model for gas radiation and absorption.

The absorption coefficient is computed as

(3.21)\[\mathrm{gas\_abs\_coeff} = \frac{R_u}{101325.} \rho T \sum_i (\mathrm{abs}_i(T) y_i/MW_i) \ \ (\mathrm{m}^{-1})\]

for the sum over the species \(\mathrm{CH}_4\), \(\mathrm{H}_2\mathrm{O}\), \(\mathrm{CO}\), \(\mathrm{CO}_2\) absorption fits (if present in the chemical mechanism) given by

(3.22)\[\begin{split}\mathrm{abs}_{\mathrm{CH}_4}(T) &= c_0 + c_1 T + c_2 T^2 + c_3 T^3 + c_4 T^4 \\ c_0 &= 6.6334 \\ c_1 &= -0.0035686\\ c_2 &= 1.6682\cdot 10^{-8}\\ c_3 &= 2.5611 \cdot 10^{-10}\\ c_4 &= -2.6558 \cdot 10^{-14}\end{split}\]
(3.23)\[\begin{split}\mathrm{abs}_{\mathrm{H}_2\mathrm{O}}(T) &= c_0 + \frac{1000}{T} (c_1 + \frac{1000}{T} (c_2 + \frac{1000}{T} (c_3 + \frac{1000}{T} (c_4 + \frac{1000}{T} + c_5)))) \\ c_0 &= -0.23093\\ c_1 &= -1.12390\\ c_2 &= 9.41530\\ c_3 &= -2.99880\\ c_4 &= 0.51382\\ c_5 &= -1.86840 \cdot 10^{-5}\end{split}\]
(3.24)\[\begin{split}\mathrm{abs}_{\mathrm{CO}_2}(T) &= c_0 + \frac{1000}{T} (c_1 + \frac{1000}{T} (c_2 + \frac{1000}{T} (c_3 + \frac{1000}{T} (c_4 + \frac{1000}{T} + c_5)))) \\ c_0 &= 18.741\\ c_1 &= -121.310\\ c_2 &= 273.500\\ c_3 &= -194.050\\ c_4 &= 56.310\\ c_5 &= -5.8169\end{split}\]
(3.25)\[\begin{split}\mathrm{abs}_{\mathrm{CO}} &= c_0 + T (c_1 + T (c_2 + T (c_3 + T c_4))) \\ c_0 &= \begin{cases} 10.09 \ \ \mathrm{if} \ \ T \geq 750 \\ 4.7869 \ \ \mathrm{otherwise} \end{cases} \\ c_1 &= \begin{cases} -0.01183 \ \ \mathrm{if} \ \ T \geq 750 \\ -0.06953 \ \ \mathrm{otherwise} \end{cases} \\ c_2 &= \begin{cases} 4.7753 \cdot 10^{-6} \ \ \mathrm{if} \ \ T \geq 750 \\ 2.95775\cdot 10^{-4} \ \ \mathrm{otherwise} \end{cases} \\ c_3 &= \begin{cases} -5.87209 \cdot 10^{-10} \ \ \mathrm{if} \ \ T \geq 750 \\ -4.25732\cdot 10^{-7} \ \ \mathrm{otherwise} \end{cases} \\ c_4 &= \begin{cases} -2.5334 \cdot 10^{-14} \ \ \mathrm{if} \ \ T \geq 750 \\ 2.02894\cdot 10^{-10} \ \ \mathrm{otherwise} \end{cases}\end{split}\]

The radiation source is then related to the absorption as

(3.26)\[\mathrm{gas\_rad\_source} = \mathrm{gas\_abs\_coeff} \sigma_B T^4 / \pi \ \ (\mathrm{MW}/\mathrm{m}^3)\]

3.3. Soot Radiation Models

Absorption and radiation soot coefficients can be added onto a library through the function tabulate_soot_radiation_model. Spitfire currently supports a ‘linear_absorption_coefficient’ model for soot radiation and absorption.

The absorption coefficient is computed as

(3.27)\[\mathrm{soot\_abs\_coeff} = \mathrm{max}(\mathrm{intercept} + \mathrm{slope} * T, \mathrm{min\_value}) \frac{\rho}{\rho_s} \ \ (\mathrm{m}^{-1})\]

where the parameters ‘intercept’, ‘slope’ and ‘min_value’ are provided by the user as desired. Default values for these parameters are ‘intercept’ = -3.75e5 1/m, ‘slope’ = 1735 1/(m-K) and ‘min_value’ = 0 1/m.

The radiation source is then related to the absorption as

(3.28)\[\mathrm{soot\_rad\_source} = \mathrm{soot\_abs\_coeff} \sigma_B T^4 / \pi \ \ (\mathrm{MW}/\mathrm{m}^3)\]

In a CFD simulation, the effects of soot on the total absorption coefficient can be represented with

(3.29)\[\mathrm{absorption\_coeff} = \mathrm{gas\_abs\_coeff} + Y_s \mathrm{soot\_abs\_coeff}\]

and the effects of soot on the total radiation source can be represented with

(3.30)\[\mathrm{radiation\_source} = \mathrm{gas\_rad\_source} + Y_s \mathrm{soot\_rad\_source}\]

3.4. Tutorials