Quadrature method of moments (QMOM)

Problem overview

Consider a continuous size distribution $n(m)$ with units of number per volume per mass. The size at any mass evolves as

$$ \frac{dn(m)}{dt} = S(n(m),m), $$

where $S(n(m),m)$ is the source term that includes nucleation, coagulation, growth, and oxidation, and depends on the size distribution $n(m)$, the particle size $m$, and the gas state (which is not explicitly shown in the function arguments here).

Soot moments of order $k$ are defined as

$$ M_k = \int n(m)m^kdm. $$

Here, $M_0$ is just the total number density of particles, #/m$^3$, and $M_1$ is the soot mass density, kg soot/m$^3$.

The evolution equation for $M_k$ is found by multiplying $dn/dt=S$ by $m^k$ and integrating over all soot sizes. The integral commutes with the derivative to give

$$ \frac{dM_k}{dt} = \int m^kS(n(m),m)dm. $$

To solve for the moments, we have to compute the integral of the source term, but the source term depends on $n(m)$. But we don’t know $n(m)$ when using the method of moments, we only know some number of moments $M_k$. That is, we need to solve equations of the form $dM_k/dt=S_M(M_k)$. The moment equation is unclosed. We use QMOM to close the equations by writing the integral of the source term in terms of $M_k$.

QMOM formulation

In general, when integrating a function $F(m)$ using quadrature, we assume the function can be split into a product $F(m)=W(m)f(m)$, where $f(m)$ is well-represented by a polynomial over the integration domain. The quadrature definition is then

$$ \int F(m) = \int W(m)f(m)dm \approx \sum_{i=0}^{N_Q-1}w_if(m_i), $$

where $N_Q$ is the number of quadrature points, $m_i$ are abscissas, or function evaluation points, and $w_i$ are corresponding weights. The problem is then to evaluate weights $w_i$ and abscissas $m_i$, after which the integral can be computed by simply evaluating the sum of $w_i$ times the function $f$ evaluated at the quadrature points $m_i$.

In application to soot formation, the weight function $W(m)$ is just the particle size distribution $n(m)$.

To evaluate $m_i$ and $w_i$, use the moments of $W(m)$. Consider the first four moments:

$$ \begin{align} M_0 &= \int n(m)m^0dm = w_0m_0^0 + w_1m_1^0, \\\\ M_1 &= \int n(m)m^1dm = w_0m_0^1 + w_1m_1^1, \\\\ M_2 &= \int n(m)m^2dm = w_0m_0^2 + w_1m_1^2, \\\\ M_3 &= \int n(m)m^3dm = w_0m_0^3 + w_1m_1^3. \end{align} $$

If the moments $M_0$ to $M_3$ are known, this gives a system of four equations in four unknowns $w_0,$ $w_1,$ $m_0,$ $m_1$ at two quadrature points.
Note that we don’t need to know $n(m),$ only the $M_k$.
The closure problem is then solved: given the moments, solve the system of equations for the $m_i$ and $w_i$, then use these to evaluate integrals for the moments.
In general, the number of quadrature points $N_Q$ is half the number of moments considered, and the number of moments is even.
The system of equations given is ill-conditioned, and specialized methods, such as the Wheeler algorithm, are used to solve it. See Marchisio and Fox Appendix A for details.

Relation to the size distribution

Note that the QMOM formulation is equivalent to assuming that the particle size distribution is given by

$$n(m)\approx\sum_i^{N_Q-1}w_i\delta(m-m_i),$$

where $\delta(m-m_i)$ is the Dirac delta function, which is a generalized function whose value is zero everwhere except where $m=m_i$, and whose integral is unity, so that $\int w_i\delta(m-m_i)dm=w_i,$ implying that $\sum_iw_i=1.$ Furthermore, $\int w_if(m)\delta(m-m_i)dm = w_if(m_i).$

In practice, this formuation means that the equation for the moment is given by

$$ \frac{dM_k}{dt} = \int m^kS(n(m),m)dm \approx\sum_{i=0}^{N_Q-1}m_i^kS(w_i,m_i). $$

This greatly simplifies the evaluation of the moment source terms, since no integration is needed and the source terms can be evaluated directly at the quadrature points as if using a sectional model with $n(m_i)$ given. That is, we replace the integral with a sum, and

wherever $n(m)$ appears, it is replaced with $w_i$
wherever $m$ appears, it is replaced with $m_i$.