A discrete soot particle size distribution (PSD) is denoted \(n_i\equiv n(m_i)\), where \(n_i\) has units of \(\#/\text{m}^3\). \(m_i\) is the mass of size \(i\). The particle size distribution evolves as a rate equation with source terms for
Moments are defined as \[M_r = \sum_{i=1}^\infty m_i^rn_i. \tag{1}\]
So called size moments are \(\mu_r = M_r/M_0\).
\(M_r\) has units of \(\text{kg}^r/\text{m}^3\).
\(\mu_r\) has units of \(\text{kg}^r/\text{particle}\).
The rate expression for the moments is \[\frac{dM_r}{dt} = \hat{N}_r + \hat{G}_r + \hat{C}_r,\] where \[\hat{N}_r = \sum_im_i^r\dot{N}_i,\]\[\hat{G}_r = \sum_im_i^r\dot{G}_i,\]\[\hat{C}_r = \sum_im_i^r\dot{C}_i.\]
Denote \(m_0\) as the mass of the smallest soot particle, assumed spherical, with surface area \(S_0\), and volume \(V_0\).
This is not ideal notation since \(m_0\ne m_{i=0}\), where \(m_{i=0}\) would correspond to a particle of size zero with \(m_{i=0}=0\). But note that, in defining \(M_r\), the sum doesn’t include \(i=0\).
Particle sizes are spaced by \(\dm\), with \(m_1=\dm\).
Here, \(\d_i\) is a shape factor (not diameter), with \(d_i\in[2/3, 1]\). When \(\d_i=2/3\), we have spherical particles, and when \(\d_i=1\), we have \(S_i/S_0=V_i/V_0= n_p\), where \(n_p\) is the number of primary particles, of size \(m_0\) in this case. That is, \(\d_i=1\) corresponds to an aggregate of nonoverlapping uniform spheres.
Balthasar and Frenklach assume \(\d_i=\da\) at any given spatial location and time. That is, at a given location and time \(\d_i\) is replaced by its ensemble average over sizes \(i\).
For the d-moment, nucleated particles are assumed spherical,
\[\hat{M}_\da = m_0^{2/3}J_{nuc}.\]
Growth
Consider particle size bins:
| i-1 | i | i+1 |
Particles in bin \(i-1\) grow into bin \(i\) as a source term, and particles in bin \(i\) grow out of \(i\) into \(i+1\) as a sink term. This is written generically as
The offset indices \(i\) and \(i-1\) in the first summation are inconvenient. For \(i=1\) the first term is in the summation is zero since \(n_0=0\). So, we can write the summation as \(\sum_{i=2}^\infty m_i^rR_{i-1}n_{i-1} = \sum_{i=1}^\infty m_{i+1}^rR_in_i\), and use \(m_{i+1}=m_i+\dm\) to give
Moment methods often include the first six moments for \(r\in[0,\,5]\). Comparing terms in the sum over \(k\) for \(i=16\) corresponding to pyrene gives relative errors ranging from 0 for \(r=1\) to 12% for \(r=5\). For \(i=32\), the relative error for \(r=5\) is 6%. For \(i=64\), corresponding to the soot nucleation size of two pyrene dimers, the relative error for \(r=5\) is 3%.
import numpy as npfrom scipy.special import binomi =64# number of carbon atoms in particle (16 is pyrene, 64 is two pyrene dimers)print(f"i = {i}\n")print(f"r = {0} is not applicable")for r inrange(1,6):print(f"r = {r}")sum=0.0for k inrange(r): term = binom(r,k)*i**kprint(f"\tk = {k}: term = {term:.1f}")sum+= term rel = np.abs((sum-term)/sum)print(f"\trelative error keeping only the last term = {rel*100:.1f}%")
i = 64
r = 0 is not applicable
r = 1
k = 0: term = 1.0
relative error keeping only the last term = 0.0%
r = 2
k = 0: term = 1.0
k = 1: term = 128.0
relative error keeping only the last term = 0.8%
r = 3
k = 0: term = 1.0
k = 1: term = 192.0
k = 2: term = 12288.0
relative error keeping only the last term = 1.5%
r = 4
k = 0: term = 1.0
k = 1: term = 256.0
k = 2: term = 24576.0
k = 3: term = 1048576.0
relative error keeping only the last term = 2.3%
r = 5
k = 0: term = 1.0
k = 1: term = 320.0
k = 2: term = 40960.0
k = 3: term = 2621440.0
k = 4: term = 83886080.0
relative error keeping only the last term = 3.1%
For oxidation, particles in bin \(i+1\) decay into bin \(i\) as a source term, and particles in bin \(i\) decay out of \(i\) into \(i-1\) as a sink term. This is written generically as
We treat the offset indices in the first summation. Since \(m_0=0\), we can start the sum at \(i=0\), and write \(\sum_{i=0}^\infty m_i^rR_{i+1}n_{i+1} = \sum_{i=1}^\infty m_{i-1}^rR_iS_i\). Use \(m_{i-1} = m_i-\dm\) to give
Use \(m_i = i\dm\), then apply the Binomial Theorem to the resulting \((i-1)^r\) term. The last term in the binomial summation will cancel with the \(i^r\) term. We then replace \(\dm^r\) with \(\dm^{r-k}dm^k\), and combine the \(dm^ki^k\) to \(m_i^k\). The result is
which only differs from Equation 7 for \(\hat{G}_r\) in the sign.
d-moment
For the \(\hat{O}_\da\), we start with Equation 11, with \(r=\da\), insert Equation 4 for \(R_i\), and use the moment definition for the resulting \(m_i^{2\da}\) term to get
This differs from Equation 8 in the \((m_i+\dm)\) versus \((m_i-\dm)\) term in the summation.
OH oxidation
For oxidation by OH, the reaction rate is based on the collision rate between OH molecules and a soot particle, with a constant efficiency factor \(\gamma_{OH}\) used. (Similarly for oxidation by O, the oxygen radical.) \(R_i\) is given by
\[R_i = \gamma_{OH}\beta_{OH}n_{OH}, \tag{13}\]
where \(\gamma_{OH}\) is the reaction efficiency, \(n_{OH}\) is the number density of OH molecules, and \(\beta_{OH,i}\) is the collision frequency in the free-molecular regime, given by
As for O\(_2\) oxidation, start with Equation 11, with \(r=\da\), insert Equation 13 for \(R_i\) using Equation 14, and use the moment definition for the resulting \(m_i^{2/3+\da}\) term to get
where \(D_A=1.395\sqrt{3}\) Å, and \(d_A\) is the size of a single aromatic ring, and 1.395 Å is the length of an aromatic C-C bond. See Frenklach and Wang in Bockhorn 1994, page 176 (note the typo in equation 10.19).
Integer moments
Insert the expression for \(\beta_{PAH,i}\) in the expression for \(\hat{P}_r\) and use the moment definition to give
The first source term creates size \(i\) by combination of smaller sizes and the second sink term removes size \(i\) by it’s collision with any other size \(j\), including \(i=j\). The factor \(1/2\) on the source term corrects for double counting on the sum.
\(\beta\) is a collision frequency and depends on whether particles are in a free-molecular \(\beta^{FM}\) or continuum \(\beta^C\) regime. The free-molecular regime corresponds to particles whose diameter is smaller than their mean free path. The continuum regime corresponds to particles whose diameter is larger than the particle mean free path.
Continuum regime
The collision frequency in the continuum regime is given by
where \(D_c\) is the collision diameter, \(\mu\) is viscosity, and \(C_c\) is the Cunninham slip correction factor. See Kazakov 1998 page 483, and Kruis 1993 for discssion of using \(D_c\) for all diameters in the equation. In writing the equation in terms of \(m_i^{1/3}\) instead of \(D_{c,i}\), the proportionality constants cancel. \(C_c\) is given by
\[C_{c,i} = 1+1.257\frac{2\lambda_f}{D_{c,i}},\]
where \(2\lambda_f/D_{c,i}\) is the Knudsen number, and \(\lambda_f\) is the gas mean free path,
where \(M\) is mean molecula weight, \(\nu\) is kinematic viscosity, and \(R_g\) is the gas constant. (See Atmospheric Chemistry and Physics 3\(^\text{rd}\) edition, by Seinfeld and Pandis, pages 365, 372.) Using \(D_{c,i} = C_aC_s m_i^{1/3}\) gives
Here, \(k_B\) is the Boltzmann constant, and \(\epsilon_c\) is a van der Waals enhancement factor, with \(\epsilon_c = 2.2\) taken from Harris and Kennedy 1988.
Integer moments
The moment source terms for coagulation are given by
In the first term, we can reorder the sums as \(\sum_{i=1}^\infty\sum_{j=1}^{i-1} = \sum_{j=1}^\infty\sum_{i=j-1}^\infty\). Then introduce \(k=i-j\) which gives \(i=k+j\), and use \(m_i=m_{k+j}=m_k+m_j,\)
Now, on the first term, change notational indices \(j\rightarrow i\), \(k\rightarrow i\), and use \((m_i+m_j)^r = \sum_{k=0}^r\binom{r}{k}m_i^km_j^{r-k}\),
These rates cannot be expanded further in terms of moments because of the fractional power on the last term in parentheses. In MOMIC, grid functions are used and interpolation is done on integer powers. For monodispersed, QMOM, and DQMOM, the terms can be evaluated directly.
d-Moment
The coagulation rate for the d-moment is neglected since coagulation is assumed to involve only point contacts.