 |
SootLib
|
|
SootLib
|
Soot nucleation arises from gas phase species reacting to form the smallest soot particles. The rate ( \(\#/\)m \(^3\)kg s) is written as $$R_n(m) = J_n\delta(m-m_n),$$ where \(J_n\,(=)\,\#/\)m \(^3\)s and \(m_n\) is mass of the soot nucleate.
For moment methods (e.g., MONO, QMOM, LOGN, MOMIC), the moment \(r\) source term is $$R_{r,n} = \int m^kR_n(m)dm = \int m^kJ_n\delta(m-m_n)dm,$$ $$ \color{#1779c4} { R_{r,n} = m_n^kJ_n. } $$
\(J_n\) is written in terms of the gas state \(T\), \(P\), and \(y_i\) and various models are available.
SootLib's most basic chemistry model is the simplified kinetic mechanism presented by Leung and Lindstedt (LL), which consists of four Arrhenius-style rate expressions: one each for soot nucleation, surface growth, oxidation, and coagulation [10]. The LL nucleation rate depends only on the concentration of gaseous acetylene, the gas temperature, and an empirically-determined rate constant with units of (#/m \(^3\)s): $$R_{nuc} = 0.1\times 10^5 e^{-21100/T} [\text{C}_2\text{H}_2]\cdot 2N_a/C_{min}.$$
Lindstedt later proposed an alteration to the LL nucleation step's pre-exponential factor to increase accuracy without changing the expression's form (LIN) [12] , giving the following expression for the rate in units of (#/m \(^3\)s): $$R_{nuc} = 0.64\times 10^{4} e^{-21100/T} [\text{C}_2\text{H}_2]\cdot 2N_a/C_{min}.$$
Lindstedt proposed a nucleation rate for soot formation from benzene [13] with the following rate expression with units of units of (#/m \(^3\)s): $$R_{nuc} = 0.75\times 10^5 e^{-21000/T} [\text{C}_6\text{H}_6]\cdot 6N_a/C_{min}.$$
SootLib implements the PAH nucleation model presented by Blanquart and Pitsch, in which nascent soot particles are created by the collision of two PAH dimers, themselves created by the collision of two PAH molecules [2]. Assuming free-molecular coagulation for the self-collision of PAH molecules, the effective rate of dimerization is given by $$\doot{n}_D = \sum_i \gamma_i \frac{1}{2} \beta_i n_i^2,$$ where $$\beta_i = 4\left( \frac{\pi k_B T}{m_i} \right)^{1/2} \left( \frac{6m_i}{\pi \rho_s} \right)^{2/3}$$ and \(k_B\) is the Boltzmann constant, \(T\) is gas temperature, \(m_i\) is the mass of PAH species \(i\), \(n_i\) is the number density of PAH species \(i\), and \(\rho_s\) is the solid soot density. The sticking coefficient \(\gamma_i\), which scales with \(m_i^4\), accounts for the difference between the observed self-collision rate and the rate predicted by kinetic molecular theory. The following table lists the PAH species considered by this PAH nucleation model, their masses \(m_i\), and their sticking coefficients \(\gamma_i\):
| PAH species | Formula | \(m_i\) (amu) | \(\gamma_i\) |
|---|---|---|---|
| Naphthalene | \(\text{C}_{10}\text{H}_8\) | 128 | 0.0010 |
| Acenaphthylene | \(\text{C}_{12}\text{H}_8\) | 152 | 0.0030 |
| Biphenyl | \(\text{C}_{12}\text{H}_{10}\) | 154 | 0.0085 |
| Phenathrene | \(\text{C}_{14}\text{H}_{10}\) | 178 | 0.0150 |
| Acephenanthrylene | \(\text{C}_{16}\text{H}_{10}\) | 202 | 0.0250 |
| Pyrene | \(\text{C}_{16}\text{H}_{10}\) | 202 | 0.0250 |
| Fluoranthene | \(\text{C}_{16}\text{H}_{10}\) | 202 | 0.0250 |
| Cyclo[cd]pyrene | \(\text{C}_{18}\text{H}_{10}\) | 226 | 0.0390 |
Rather than computing the properties of every possible dimer combination, the model evaluates only the total rate of dimer formation \(\doot{n}_D\) and the average carbon content of each dimer. Because PAH and dimer collision rates are high, we compute dimer concentration assuming a quasi-steady state. The PAH dimer concentration can be found by solving $$\doot{n}_D = \beta_{D,D} n_D^2 + \left( \sum_{j} \beta_{j,D} n_j \right) n_D $$ for the steady state dimer concentration \(n_D\). The second term on the right hand side of the equation represents the condensation of PAH dimers onto the surface of soot particles, where \(n_j\) is the number density of soot particles of size \(j\) and \(\beta_{j,D}\) is the collision coefficient as applied to the interaction of soot particles and PAH dimers. Once the steady state dimer concentration has been computed, the nucleation rate is given by $$R_{nuc} = \frac{1}{2}\beta_{D,D} n_D^2.$$
This PAH nucleation model also accounts for condensation of PAH dimers onto existing soot particles, which contributes to chemical surface growth. The PAH condensation rate is computed similarly to the nucleation rate, but \(\beta\) is computed as the collision rate between PAH dimers and soot particles rather than the self-collision rate of PAH dimers, and the mass concentration of both the PAH dimers and soot particles is taken into account.
For a continuous particle size distribution (assumed here), the growth rate ( \(\#/\)m \(^3\)kg s) is given by $$R_g(m) = -\frac{\partial}{\partial m}(v_gn(m)),$$ where \(v_g\) is a velocity in the mass coordinate (the so-called internal coordinate) of the PSD and has units of kg/s. \(v_g\) is given by $$v_g = k_sS(m),$$ where \(k_s\) is the mass growth rate per surface area \(k_s\,(=)\,\)kg/m \(^2\)s, and \(S(m)\) is the soot surface area per particle. The expression for \(R_g(m)\) follows from $$R_g(m) = \lim_{\Delta m\rightarrow 0}\frac{k_s}{\Delta m}(N_{i-1}S_{i-1}-N_iS_i),$$ where \(N_i\) is a discrete particle number density \(\#/\)m \(^3\), and \(S_i=S(m_i)\).
The moment source term is $$R_{r,g} = \int m^kR_n(R)dm = -\int m^k\frac{\partial}{\partial m}(v_gn(m))dm.$$ Integrating by parts and noting that \(n(m)=0\) at the bounds of the integration gives $$R_{r,g} = k\int_0^\infty v_gm^{k-1}n(m)dm.$$ Now, use \(v_g = k_sS(m)\) with \(S = \pi d^2\), \(m=\rho_s\pi d^3/6\), or $$S(m) = \pi\left(\frac{6}{\pi\rho_s}\right)^{2/3}m^{2/3}.$$ Then we have $$ \color{#1779c4} { R_{r,g} = k_s\pi\left(\frac{6}{\pi\rho_s}\right)^{2/3}kM_{k-1/3}. } $$ Here, the \(M_{k-1/3}\) makes it clear that we will need expressions for fractional moments, and evaluation of such moments requires moment closure.
Most soot surface growth models rely on acetylene ( \(\text{C}_2\text{H}_2\)) as the primary source of gaseous carbon, though other species may also contribute. Additionally, surface growth models tend to include some dependence on the soot particle's surface area since the availability of sites for the addition of carbon atoms to existing particles tends to be a limiting factor in rate calculations [19].
The Leung and Lindstedt mechanism for soot surface growth (LL) [10] depends on both the gaseous concentration of acetylene and the particle surface area available for surface reactions. The overall rate of soot surface growth is computed as $$R_{grw} = 0.6\times 10^4 e^{-21100/T} f(A_s) [\text{C}_2\text{H}_2].$$ \(f(A_s)\) is a function of the available particle surface area given by $$f(A_s) = \sqrt{\pi \left( \frac{6MW_C}{\pi \rho_{s}} \right) ^{2/3}} \left[ \frac{\rho Y_{s}}{MW_C} \right]^{1/3} [\rho N]^{1/6},$$ where \(MW_{\text{C}}\) is the molar mass of carbon, \(\rho\) is the gas density, \(\rho_{s}\) is the solid density of soot, \(Y_{s}\) is the mass fraction of soot, and \(N\) is the soot particle number density. Results showed that the normal square dependence of the rate on available surface sites does not accurately predict the soot number density throughout the flame [10], and later models explored alternative methods.
Lindstedt proposed and tested several alternative models for the surface growth rate, concluding that the model in which the growth rate is dependent on soot number density but independent of soot surface area produced better results [13]. This modified growth rate (LIN) is given by $$R_{grw} = 750 e^{-12100/T} [\text{C}_2\text{H}_2] 2N_s MW_C,$$ where \(N_s\) is the number density of soot particles and \(MW_{\text{C}}\) is the molar mass of carbon.
The HACA growth mechanism consists of a repeating sequence of elementary reaction steps with individual Arrhenius-style rate expressions, listed in the following table.
| Reaction | Rate constant | |
|---|---|---|
| (1) | C(s)-H + H \(\leftrightarrow\) C(s) \(^.\) + H \(_2\) | \( k_{1} = 4.2\times 10^{13} e^{-13/RT}\) |
| \( k_{r1} = 3.9\times 10^{12} e^{-11/RT} \) | ||
| (2) | C(s)-H + OH \(\leftrightarrow\) C(s) \(^.\) + H \(_2\)O | \( k_{2} = 1.0\times 10^{10} T^{0.734} e^{-1.43/RT}\) |
| \( k_{r2} = 3.68\times 10^{8} T^{1.139} e^{-17.1/RT} \) | ||
| (3) | C(s) \(^.\) + H \(\rightarrow\) C(s)-H | \( k_{3} = 2.0\times 10^{13}\) |
| (4) | C(s) \(^.\) + C \(_2\)H \(_2\) \(\rightarrow\) C(s)-H + H | \( k_{4} = 8.0\times 10^{7} T^{1.56} e^{-3.8/RT} \) |
| (5) | C(s) \(^.\) + O \(_2\) \(\rightarrow\) 2CO + products | \( k_{5} = 2.2\times 10^{12} e^{-7.5/RT}\) |
| (6) | C(s)-H + OH \(\rightarrow\) CO + products | \( k_{6} = 0.13\cdot 1290 T^{-1/2} \) |
Reactions~(1)–(4) represent surface growth, while reactions (5) and (6) represent oxidation. Reactions~(1)–(5) are required for calculation of the steady state radical site number density calculation, as follows.
The HACA model assumes that the reactive \(\text{C-H}\) surface sites on a soot molecule are located at aromatic bays, where steric hindrance and molecular deformation weaken the \(\text{C-H}\) bonds at the bay sites [1]. The number density of \(\text{C-H}\) surface sites is estimated to be \(\chi_{C-H}=2.3\times 10^{15}\) sites\m \(^2\), and the the number density of radical surface sites \(\chi_{C^.}\) is calculated assuming steady state, which leads to $$\chi_{C^.} = \left[ \frac{k_1[\text{H}] + k_2[\text{OH}]}{k_{r1}[\text{H}_2] + k_{r2}[\text{H}_2\text{O}] + k_3[\text{H}] + k_4[\text{C}_2\text{H}_2] + k_5[\text{O}_2]} \right] \chi_{C-H},$$ where \(k_i\) represents the rate constant of reaction \(i\) in the sequence above and subscripts prefixed with \(r\) indicate reverse reactions [3] [1]. Once \(\chi_{C^.}\) has been computed, the soot surface growth rate by the HACA mechanism is given by the forward rate of Reaction (4) in the table above: $$R_{grw} = k_4 \chi_{C^.} [\text{C}_2\text{H}_2].$$
Oxidation follows the same form as growth but the positive \(k_s\) is replaced with the negatively signed \(k_o\).
Similar to surface growth, oxidation mechanisms often depend on the available surface area of oxidizing soot particles, which may or may not be a limiting factor depending on the thermodynamic state, the composition of the gas mixture, and the amount of soot present. Early soot oxidation models tend to rely on \(\text{O}_2\) as the principal oxidant, though experimental studies show that OH and sometimes O can also contribute significantly to soot oxidation rates.
The global rate expression presented by Lee et al. [8] is $$R_{oxi} = 1.085\times 10^4 P_{\text{O}_2} T^{-1/2} e^{-1.977\times 10^4/T},$$ where \(P_{\text{O}_2}\) is the partial pressure of \(\text{O}_2\).
In the model presented by Nagle and Strickland-Constable, the oxidation rate depends on a nonlinear combination of Arrhenius-style rate constants [Nagle_1962:] $$R_{oxi} = \rho_{s} \left[ k_A P_{\text{O}_2} \left( \frac{x}{1+k_Z P_{\text{O}_2}}\right) + k_B P_{\text{O}_2} (1-x) \right]$$ where $$x=\frac{1}{1+\frac{k_T}{k_B P_{\text{O}_2}}}$$ and $$k_A = 20e^{-15098/T},$$ $$k_B = 4.46\times 10^{-3} e^{-7650/T},$$ $$k_T = 1.51\times 10^5 e^{-48817/T},$$ $$k_Z = 21.3e^{2063/T}.$$
The Leung and Lindstedt oxidation expression is based on the earlier Lee et al. model, but uses a different temperature dependence and modified empirical factors [Leung_1991:] $$R_{oxi} = 0.1\times 10^5 T^{1/2} e^{-19680/T} [\text{O}_2].$$
These early soot oxidation models either do not account for the influence of oxidation by OH or account for it empirically as part of the rate constant. Neoh et al. use the expression $$R_{oxi}=0.13\cdot 1290 P_{\text{OH}} T^{-1/2}$$ to represent soot oxidation by OH [16] [15].
To account for the lack of consideration for oxidation by OH in the previous models, SootLib adds the preceding expression to those presented by Lee et al. and Nagle and Strickland-Constable, resulting in the combined LEE_NEOH and NSC_NEOH options for soot oxidation in SootLib. Leung and Lindstedt explicitly acknowledge the lack of oxidation by OH in their model and consider their one-step oxidation mechanism sufficient for their purposes; in light of their reasoning, SootLib does not alter or add to the existing Leung and Lindstedt oxidation step.
Guo, Anderson, and Sunderland [6] examined 12 experimental studies of soot oxidation and found optmized rates for oxidation from OH, and O2. The OH rate is the same as that of Neoh, but with a slightly different reaction efficiency of 0.1, $$R_{oxi,OH}=0.1\cdot 1290 P_{\text{OH}} T^{-1/2}.$$ The oxidation rate by O2 is given by $$R_{oxi,O2} = 15.8P_{O2}T^{1/2}e^{-23450/T}/1000,$$ with units of (kg/m \(^2\)s), where the oxygen pressure is in Pascals.
Josephson et al. [7] examined 13 experimental studies of soot oxidation and found optmized rates for oxidation from OH, and O2 using Bayesian methods. The OH rate is the same as that of Neoh, but with a different reaction efficiency of 0.15, $$R_{oxi,OH}=0.15\cdot 1290 P_{\text{OH}} T^{-1/2}.$$ The oxidation rate by O2 is given by $$R_{oxi,O2} = 0.798P_{O2}T^{1/2}e^{-21294/T}/1000,$$ with units of (kg/m \(^2\)s), where the oxygen pressure is in Pascals.
SootLib also includes the HACA mechanism for soot oxidation, represented by reactions (5) and (6) in the HACA table of reactions in the surface growth chemistry section. The HACA mechanism accounts for oxidation by both \(\text{O}_2\) and OH and includes the previously-discussed dependence on available particle surface sites. The HACA rates of oxidation by \(\text{O}_2\) and OH are given by [1] $$R_{oxi} = k_5 \chi_{C^.} [\text{O}_2]$$ $$R_{oxi} = k_6 \chi_{C-H} P_{OH}.$$
Coagulation here is based on the presentation in [18]. See also [5].
Particle coagulation between two particles of masses \(m\) and \(\mu\) results in the rate ( \(\#/\)m \(^3\)kg s) $$R_c(m) = \frac{1}{2}\int_0^m\beta(\mu, m-\mu)n(\mu)n(m)d\mu - \int_0^\infty\beta(m,\mu)n(\mu)n(m)d\mu.$$ The first term represents creation of particles of size \(m\) by particles smaller than size \(m\), and the second term represents destruction of particles of size \(m\) by collision of size \(m\) particles with all particles \(\mu\). \(\beta\) is the coagulation coefficient.
The coagulation rate for a moment \(M_r\) is $$ R_{r,c} = \frac{1}{2}\iint\beta(m,\mu)n(m)n(\mu)\left[-2m^k+\sum_{j=0}^r\left(r\atop j\right)m^j\mu^{r-j}\right]dmd\mu. $$ This simplifies to $$ \color{#1779c4} { \begin{align} &R_{0,c} = -\frac{1}{2}\iint\beta(m,\mu)n(m)n(\mu)dmd\mu, \\ &R_{1,c} = 0, \\ &R_{r\ge 2,c} = \frac{1}{2}\iint\beta(m,\mu)n(m)n(\mu)\sum_{j=1}^{r-1}\left(r\atop j\right)m^j\mu^{r-j}dmd\mu. \end{align} } $$
The coagulation coefficient \(\beta\) is generally written for two regimes: free molecular, and continuum, depending on ratio of the mean free path of the particles to the particle diameter (particle Knudsen number).
The coagulation coefficient in the free molecular regime (large Knudsen number) is given by $$ \color{#1779c4} { \beta^{FM}(m,\mu) = \epsilon_c\left(\frac{\pi k_BT}{2}\right)\left(\frac{6}{\pi\rho_s}\right)^{2/3}\left(\frac{1}{m}+\frac{1}{\mu}\right)^{1/2}(m^{1/3}+\mu^{1/3})^2. } $$ Here, \(k_B\) is the Boltzmann constant, \(\rho_s\) is the soot particle density, and \(\epsilon_\text{C}\) is a van der Waals enhancement factor ( \(\epsilon_c = 2.2\)). The term \((1/m+1/\mu)^{1/2}\) is the reduced mass \(m_1m_2/(m_1+m_2)\), and \((6/\pi\rho_s)^{2/3}(m^{1/3}+\mu^{1/3})^2 = (d_m+d_\mu)^2\), where \(d_m\) and \(d_\mu\) are the diameters of particles with masses \(m\) and \(\mu\), respectively.
The coagulation coefficient in the continuum regime (small Knudsen number) is given by $$ \color{#1779c4} { \beta^C(m,\mu) = \frac{2k_BT}{3\mu_v}\left(\frac{C_m}{m^{1/3}}\frac{C_\mu}{\mu^{1/3}}\right)(m^{1/3}+\mu^{1/3}), } $$ where \(\mu_v\) is the gas dynamic viscosity, and \(c_m\) (and similarly for \(c_\mu\)) is $$ \color{#1779c4} { C_m = 1+1.657Kn_m, } $$ $$ \color{#1779c4} { Kn_m = 2\lambda_g/d_m, } $$ $$ \color{#1779c4} { \lambda_g = \nu\left(\frac{\pi M}{2RT}\right)^{1/2}, } $$ where \(C_m\) is the Cunningham slip correction factor, \(\lambda_g\) is the gas mean free path, \(\nu\) is the kinematic viscosity, \(M\) is the gas mean molecular weight, and \(R\) is the gas constant.
A harmonic mean can be used to transition between the FM and C regimes: $$ \color{#1779c4} { \beta^{HM} = \frac{\beta^{FM}\beta^C}{\beta^{FM}+\beta^C}. } $$
Fuchs uses the following form that transitions between the FM and C regime. $$ \color{#1779c4} { \beta^{F}(m,\mu) = 2\pi(D_m+D_\mu)(d_m+d_\mu)\left( \frac{d_m+d_\mu}{d_m+d_\mu+2(g_m^2+g_\mu^2)^{1/2}} + \frac{8(D_m+D_\mu)/\epsilon_c}{(\bar{c}_m^2+\bar{c}_\mu^2)^{1/2}(d_m+d_\mu)} \right)^{-1}, } $$ $$ \color{#1779c4} { \bar{c}_m = \left(\frac{8k_BT}{\pi m}\right)^{1/2}, } $$ $$ \color{#1779c4}{ g_m = \frac{\sqrt{2}}{3d_ml}\left[(d_m+l_m)^3-(d_m^2+l_m^2)^{3/2}\right] - \sqrt{2}d_m, } $$ $$ \color{#1779c4} { l_m = \frac{8D_m}{\pi\bar{c}_m}, } $$ $$ \color{#1779c4} { D_m = \frac{k_BTC_c}{3\pi\mu_vd_m}, } $$ where \(D_m\) is particle diffusivity, \(l_m\) is a particle mean free path, \(\bar{c}_m\) is the mean particle velocity, and \(g_m\) is defined in [5] page 291. (Note that Seinfeld erroneously leaves off the \(\sqrt{2}\) in the second term of \(g_m\).)
Here is a link to a Jupyter notebook with coagulation kernels plotted.
1.9.5