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%
Method of Moments
\[ \newcommand\tprime{\prime\prime\prime} \newcommand{\prtl}[2]{\frac{\partial #1}{\partial #2}} \newcommand\d{\mathbf{d}} \newcommand\da{{\langle\mathbf{d}\rangle}} %\newcommand\dm{{\Delta m}} \newcommand\dm{{\Delta_m}} \]
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
- Nucleation, \(\dot{N}_i\)
- Growth, \(\dot{G}_i\)
- Coagulation, \(\dot{C}_i\)
Moments are defined as
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
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\).
Mass, volume, surface area
Spherical particles
\[m_i = \rho_sV_i,\] \[V_i = \frac{\pi}{6}d_i^3,\] \[S_i = \pi d_i^2,\] \[d_i = \left(\frac{6m_i}{\pi\rho_s}\right)^{1/3},\] \[S_i = \pi\left(\frac{6}{\pi\rho_s}\right)^{2/3}m_i^{2/3}.\]
Define
Non-spherical particle aggregates
Following Balthasar and Frenklach 2005,
\[S_i = S_0\left(\frac{V_i}{V_0}\right)^{\d_i} = S_0\left(\frac{m_i}{m_0}\right)^{\d_i}.\]
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\).
Balthasar and Frenklach further assume \[\mu_\da = \mu_1^\da,\] where \[\mu_\da = M_\da/M_0.\] This gives
Collision diameter
Particle collision rates depend on the particle size. For non-spherical particles, this is given as
where \(D_{c,d23}\) is the radius of a volume-equivalent sphere, that is,
In the following, in all cases, setting \(\da=2/3\) results in a spherical particle assumption.
Nucleation
Nucleation to the smallest size is given by
where \(\delta\) is the Dirac function. \(J_{nuc}\) has units of \(\text{particles}/\text{m}^3\text{s}\).
The moment source term is given by
For the d-moment, nucleated particles are assumed spherical,
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
\(R_i\) is a rate of formation of particles of size \(i\), per particle. For reaction-rate-limited growth, e.g., by HACA, \(R_i\) is given by
\[R_i = \frac{k_s}{\dm}S_i = \frac{k_s}{\dm}\pi C_s^2m_0^{\frac{2}{3}-\da} m_i^\da. \tag{4}\]
The moment rate \(\hat{G}_r\) is given by
\[\hat{G}_r = r\dm\sum_{i=1}^\infty m_i^{r-1}R_in_i. \tag{5}\]
Derivation
This is derived as follows. Apply the moment definition to Equation 3:
\[\hat{G}_r = \sum_{i=1}^\infty m_i^rR_{i-1}n_{i-1} - \sum_{i=1}^\infty m_i^rR_in_i.\]
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
\[\hat{G}_r = \sum_{i=1}^\infty\left[(m_i+\dm)^r-m_i^r\right]R_in_i. \tag{6}\]
Also, write \(m_i^r=\dm^ri^r\), and \(m_{i+1}^r = \dm^r(i+1)^r\) to get
\[\hat{G}_r = \dm^r\sum_{i=1}^\infty\left[(i+1)^r-i^r\right]R_in_i.\]
Use the Binomial Theorem to write
\[\hat{G}_r = \dm^r\sum_{i=1}^\infty\left[\sum_{k=0}^r\binom{r}{k}i^k-i^r\right]R_in_i.\]
The last term of the inner sum cancels with the \(i^r\) term. Replace \(\dm^r\) with \(\dm^{r-k}\dm^k\), and \(\dm^ki^k=m_i^k\),
\[\hat{G}_r = \sum_{i=1}^\infty\sum_{k=0}^{r-1}\binom{r}{k}\dm^{r-k}m_i^kR_in_i.\]
Only the last term in the sum over \(k\) is significant for small \(\Delta m\), which gives Equation 5.
The assumption retaining the last term in the sum is demonstrated below for integer moments using Equation 4. In that case, the growth rate is
\[\hat{G}_r \sim \sum_{i=0}^\infty\sum_{k=0}^{r-1}\binom{r}{k} \frac{m_i^{k+\da}n_i}{\dm^k} \sim \sum_{i=0}^\infty i^\da n_i \sum_{k=0}^{r-1}\binom{r}{k}i^k.\]
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%.
Integer moments
Using Equation 4 in Equation 5 along with the moment definition gives
d-moment
For \(\hat{G}_\da\) a different expression than Equation 7 is needed since \(\da\) is not an integer, unlike \(r\).
Start with Equation 6 with \(r=\da\), insert Equation 4 for \(R_i\), and use the moment definition for the resulting \(m_i^{2\da}\) term to get
Oxidation
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
The moment rate \(\hat{O}\) is given by
\[\hat{O} = -r\dm\sum_{i=1}^\infty m_i^{r-1}R_in_i. \tag{10}\]
This differs from Equation 5 for \(\hat{G}_r\) only in the sign, as expected.
Derivation
The derivation is similar to that of the growth rate; an outline is given below.
\[\hat{O} = \sum_{i=1}^\infty m_i^rR_{i+1}n_{i+1} - \sum_{i=1}^\infty m_i^rR_in_i.\]
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
\[\hat{O} = \sum_{i=1}^\infty\left[(m_i-\dm)^r - m_i^r\right]R_in_i. \tag{11}\]
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
\[\hat{O} = \sum_{i=1}^\infty\sum_{k=0}^{r-1}\binom{r}{k}\dm^{r-k}m_i^k(-1)^{r-k}R_in_i.\]
Note the alternating signs of the terms in the summation over \(k\). Only the last term of the sum over \(k\) is significant, giving Equation 10.
Oxidation can be reaction-rate-limited, e.g., by O\(_2\), or collisional, e.g., by OH or O.
O\(_2\) oxidation
Integer moments
For O\(_2\) oxidation, we can use Equation 4 for \(R_i\) in Equation 10, along with the moment definition to give
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
\[\beta_{OH,i} = \sqrt\frac{\pi k_BT}{2m_{OH}}C_s^2C_a^2m_i^{2/3}, \tag{14}\]
where \(k_B\) is the Boltzmann constant. This collision rate is derived from the more general expression assuming \(m_{OH}\ll m_i\).
Integer moments
Use the above relations for \(R_i\) and \(\beta_{OH,i}\) in Equation 10:
d-moment
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
PAH Condensation
PAH condensation is assumed to occur in the free-molecular regime. We use Equation 5 with \(\hat{P}_r=\hat{G}_r\) and \(R_i=\beta_{PAH,i}n_{PAH}\):
\[\hat{P}_r = r\dm\sum_im_i^{r-1}n_i\beta_{PAH,i}n_{PAH}.\]
The free-molecular coagulation rate constant is
\[\beta_{PAH,i} = \epsilon_c\sqrt{\frac{\pi k_BT}{2}}\left(\frac{1}{m_{PAH}}+\frac{1}{m_i}\right)^{1/2}(D_{PAH} + D_i)^2.\]
We assume \(m_{PAH}\ll m_i\) in the reduced mass term, the first term in parentheses, and use the collision diameter for \(D_i\) to give
\[\beta_{PAH,i} = \epsilon_c\sqrt{\frac{\pi k_BT}{2m_{PAH}}}(D_{PAH} + C_sC_am_i^{\frac{1}{3}})^2.\]
Frenklach models the collision diameter of PAH as
\[D_{PAH} = D_A\sqrt{\frac{2m_{PAH}}{3m_C}},\]
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
d-moment
Start with Equation 6 with \(\hat{P}_r=\hat{G}_r\), \(r=\da\), \(R_i=\beta_{PAH,i}n_{PAH}\), and using the moment definition to obtain
Coagulation
The coagulation rate is given by
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
\[\begin{align} \beta_{i,j}^C &= \frac{2k_BT}{3\mu}(D_{c,i}+D_{c,j})\left(\frac{C_{c,i}}{D_{c,i}} + \frac{C_{c,j}}{D_{c,j}}\right), \\ &= \frac{2k_BT}{3\mu}(m_i^{1/3}+m_j^{1/3})\left(\frac{C_{c,i}}{m_i^{1/3}} + \frac{C_{c,j}}{m_j^{1/3}}\right), \end{align}\]
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,
\[\lambda_f = \nu\left(\frac{\pi M}{2R_gT}\right)^{1/2},\]
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
where
Free molecular regime
The collision frequency in the free-molecular regime is given by
where
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
\[\hat{C}_r = \frac{1}{2}\sum_{i=1}^\infty\sum_{j=1}^\infty\sum_{k=1}^{r-1}\binom{r}{k}m_i^km_j^{r-k}\beta_{i,j}n_in_j. \tag{17}\]
Derivation
\[\begin{align} \hat{C}_r = \sum_{i=1}^\infty m_i^r\dot{C}_i &= \sum_{i=1}^\infty m_i^r\left(\frac{1}{2}\sum_{j=1}^{i-1}\beta_{j,i-j}n_jn_{i-j} - \sum_{j=1}^\infty\beta_{i,j}n_in_j\right), \\ &= \frac{1}{2}\sum_{i=1}^\infty\sum_{j=1}^{i-1}m_i^r\beta_{j,i-j}n_jn_{i-j} - \sum_{i=1}^\infty\sum_{j=1}^\infty m_i^r\beta_{i,j}n_in_j. \end{align}\]
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,\)
\[\hat{C}_r = \frac{1}{2}\sum_{j=1}^\infty\sum_{k=1}^\infty (m_k+m_j)^r\beta_{j,k}n_jn_k - \sum_{i=1}^\infty\sum_{j=1}^\infty m_i^r\beta_{i,j}n_in_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}\),
\[\hat{C}_r = \frac{1}{2}\sum_{i=1}^\infty\sum_{j=1}^\infty\sum_{k=0}^r\binom{r}{k}m_i^km_j^{r-k}\beta_{i,j}n_in_j - \sum_{i=1}^\infty\sum_{j=1}^\infty m_i^r\beta_{i,j}n_in_j. \tag{18}\]
The \(k=0\) and \(k=r\) components of the first term (combined), cancel with the second term, giving Equation 17.
Continuum regime
Inserting Equation 15 into Equation 17 and using the moment definition, gives
For \(r=1\) we have \(\hat{C}_1=0\), since mass is conserved.
When \(r=0\), we use Equation 18, where the sum is one term with \(k=0\), and \(\binom{0}{0}=1\). This results in
Free molecular regime
Inserting Equation 16 into Equation 17 and using the moment definition gives
For \(r=0\), we insert Equation 16 into Equation 18, giving
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.