Method of Moments

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

\[\frac{dn_i}{dt} = \dot{N}_i + \dot{G}_i + \dot{C}_i.\]

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\).

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\).

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 \[C_s = \left(\frac{6}{\pi\rho_s}\right)^{1/3}.\]

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\).

\[S_i = S_0\left(\frac{m_i}{m_0}\right)^\da = \pi C_s^2m_0^{\frac{2}{3}-\da}m_i^\da. \tag{2}\]

Balthasar and Frenklach further assume \[\mu_\da = \mu_1^\da,\] where \[\mu_\da = M_\da/M_0.\] This gives

\[\da = \log\mu_\da/\log\mu_1.\]

Collision diameter

Particle collision rates depend on the particle size. For non-spherical particles, this is given as

\[\begin{align} &D_{c,d} = C_aD_{c,d23}, \\ &C_a = (3-3\da) + C_d(3\da-2), \\ &C_d = 1.9124, \end{align}\]

where \(D_{c,d23}\) is the radius of a volume-equivalent sphere, that is, \[D_{c,d23,i} = C_sm_i^{1/3}.\]

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

\[\dot{N}_i = J_{nuc}\delta(m_i-m_0),\]

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

\[\hat{M}_r = \sum_im_i^r \dot{N}_i = m_0^rJ_{nuc}.\]

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

\[\dot{G}_i^+ = R_{i-1}n_{i-1} - R_in_i. \tag{3}\]

\(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%.

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%

Integer moments

Using Equation 4 in Equation 5 along with the moment definition gives

\[\hat{G}_r = k_s\pi C_s^2m_0^{\frac{2}{3}-\da} rM_{r-1+\da}. \tag{7}\]

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

\[\hat{G}_\da = \frac{k_s}{\dm}\pi C_s^2m_0^{\frac{2}{3}-\da}\left[\left(\sum_{i=1}^\infty (m_i+\dm)^\da m_i^\da n_i\right) - M_{2\da}\right]. \tag{8}\]

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

\[\dot{O}_i = R_{i+1}n_{i+1} - R_in_i. \tag{9}\]

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

\[\hat{O}_r^{O_2} = -k_s\pi C_s^2m_0^{\frac{2}{3}-\da} rM_{r-1+\da}, \tag{12}\]

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

\[\hat{O}_\da^{O_2} = \frac{k_s}{\dm}\pi C_s^2m_0^{\frac{2}{3}-\da}\left[\left(\sum_{i=1}^\infty (m_i-\dm)^\da m_i^\da n_i\right) - M_{2\da}\right].\]

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:

\[\hat{O}_r^{OH} = -\gamma_{OH}\sqrt{\frac{\pi k_BT}{2m_{OH}}}n_{OH}r\dm C_s^2C_a^2M_{r-\frac{1}{3}}.\]

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

\[\hat{O}_\da^{OH} = \gamma_{OH}\sqrt{\frac{\pi k_BT}{2m_{OH}}}n_{OH}C_s^2C_a^2\left[\left(\sum_{i=1}^\infty(m_i-\dm)^\da m_i^{\frac{2}{3}}n_i\right) - M_{\frac{2}{3}+\da}\right].\]

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

\[\hat{P}_r = \epsilon_c\sqrt{\frac{\pi k_BT}{2m_{PAH}}}r\dm n_{PAH}\left(D_{PAH}^2M_{r-1} + 2D_{PAH}C_sC_aM_{r-\frac{2}{3}} + C_s^2C_a^2M_{r-\frac{1}{3}}\right).\]

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

\[\begin{align} \hat{P}_\da = \epsilon_c\sqrt{\frac{\pi k_bT}{2m_{PAH}}}n_{PAH} & \left[\left( \sum_{i=1}^\infty n_i(m_i+\dm)^\da\left(D_{PAH}^2 + 2D_{PAH}C_sC_am_i^{\frac{1}{3}}+C_s^2C_a^2m_i^{\frac{2}{3}}\right) \right)\right. - \\ & \left. \phantom{x}\left(D_{PAH}^2M_\da + 2D_{PAH}C_sC_aM_{\frac{1}{3}+\da} + C_s^2C_a^2M_{\frac{2}{3}+\da}\right) \right]. \end{align}\]

Coagulation

The coagulation rate is given by

\[\dot{C}_i = \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.\]

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

\[\beta_{i,j}^C = K_C\left[ 2+\frac{m_i^{1/3}}{m_j^{1/3}}+\frac{m_j^{1/3}}{m_i^{1/3}} + K_C^\prime\left(\frac{1}{m_i^{1/3}} + \frac{1}{m_j^{1/3}} + \frac{m_i^{1/3}}{m_j^{2/3}} + \frac{m_j^{1/3}}{m_i^{2/3}}\right) \right], \tag{15}\]

where

\[K_C = \frac{2k_BT}{3\mu},\] \[K_C^\prime = \frac{2.514\lambda_f}{C_aC_s}.\]

Free molecular regime

The collision frequency in the free-molecular regime is given by

\[\beta^{FM}_{i,j} = K_{FM}\left(\frac{1}{m_i}+\frac{1}{m_j}\right)^{1/2}(m_i^{1/3} + m_j^{1/3})^2, \tag{16}\]

where

\[K_{FM} = \epsilon_c\sqrt{\frac{\pi k_BT}{2}}C_s^2C_a^2.\]

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

\[\begin{align} \hat{C}_{r>1}^C = \frac{K_C}{2}\sum_{k=1}^{r-1}\binom{r}{k}&\left[2M_kM_{r-k} + M_{k+\frac{1}{3}}M_{r-k-\frac{1}{3}} + M_{k-\frac{1}{3}}M_{r-k+\frac{1}{3}} + \right. \\ &\left. K_C^\prime\left( M_{k-\frac{1}{3}}M_{r-k} + M_kM_{r-k-\frac{1}{3}} + M_{k+\frac{1}{3}}M_{r-k-\frac{2}{3}} + M_{k-\frac{2}{3}}M_{r-k+\frac{1}{3}} \right)\right]. \end{align}\]

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

\[\hat{C}_0^C = -K_C\left[ M_0^2 + M_\frac{1}{3}M_{-\frac{1}{3}} + K_C^\prime(M_0M_{-\frac{1}{3}} + M_\frac{1}{3}M_{-\frac{2}{3}}) \right].\]

Free molecular regime

Inserting Equation 16 into Equation 17 and using the moment definition gives

\[\hat{C}_{r>1}^{FM} = \frac{K_{FM}}{2}\sum_{i=0}^\infty\sum_{j=0}^\infty\sum_{k=1}^{r-1}\binom{r}{k}m_i^km_j^{r-k}\left(m_i^{1/3}+m_j^{1/3}\right)^2\left(\frac{1}{m_i}+\frac{1}{m_j}\right)^{1/2}n_in_j.\]

For \(r=0\), we insert Equation 16 into Equation 18, giving

\[\hat{C}_0^{FM} = -\frac{K_{FM}}{2}\sum_{i=0}^\infty\sum_{j=0}^\infty\left(m_i^{1/3}+m_j^{1/3}\right)^2\left(\frac{1}{m_i}+\frac{1}{m_j}\right)^{1/2}n_in_j.\]

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.