Local Fixed Pivot Quadrature Method of Moments for Solution of Population Balance Equation

Abstract

A local fixed pivot quadrature method of moments (LFPQMOM) is proposed for the solution of the population balance equation (PBE) for the aggregation and breakage process. First, the sectional representation for aggregation and breakage is presented. The continuous summation of the Dirac Delta function is adopted as the discrete form of the continuous particle size distribution in the local section as performed in short time Fourier transformation (STFT) and the moments in local sections are tracked successfully. Numerical simulation of benchmark test cases including aggregation, breakage, and aggregation breakage combined processes demonstrate that the new method could make good predictions for the moments along with particle size distribution without further assumption. The accuracy in the numerical results of the moments is comparable to or higher than the quadrature method of moment (QMOM) in most of the test cases. In theory, any number of moments can be tracked with the new method, but the computational expense can be relatively large due to many scalar equations that may be included.

1. Introduction

Determination of the interfacial area with high accuracy between different phases in a dispersed system is critical to the prediction of the flow behaviors and mass transfer. Population balance equation (PBE)—as an essential tool to describe the multiphase system and capable of predicting the interfacial area by tracking the particle size distribution and describing the micro-behaviors that influence the particle size distribution of the disperse phase—has been widely used in scientific and engineering fields. Unfortunately, such an equation is rather complex, and the numerical method is the only choice in most cases. Several predominate numerical methods have been proposed previously, which can be divided broadly into three categories: discretization methods (DM) [1,2,3,4,5,6], methods of moment (MOM) [7,8,9,10,11,12], and stochastic methods (SM) [13,14,15]. However, it is difficult for these methods to track both particle size distribution and its moments with high accuracy. In a recent paper, a promising numerical scheme—the sectional quadrature method of moment (SQMOM)—was proposed by Attarakih et al. for solving the PBEs for aggregation and breakage processes [16]. With the new scheme, both the particle size distribution and its moments could be tracked. However, several problems accompany these new methods. (1) In SQMOM, the semi-infinite solution domain [0, ∞) is truncated into a finite domain [V_min, V_max]. Using a finite domain to represent an infinite domain will inevitably introduce errors, which are referred to as finite domain errors. As a result, the moments cannot be tracked exactly. This is inherent in almost all DMs. In this sense, SQMOM is a kind of DM more than MOM. (2) If a new born particle with a volume smaller than V_min or larger than V_max is introduced by a certain event (aggregation or breakage), the event will be omitted. Therefore, the new method may influence the physical models in some extreme situations. (3) In SQMOM, the product-different (PD) algorithm [17] must be employed to obtain the abscissas and the weights in each section when a relatively large number of moments (≥6) need to be tracked. However, in PD, the values of abscissas are not restricted to the section they should belong to in mathematics. Hence, with the PD algorithm, the abscissas obtained might be outside the section they originally belong to. Until now, there is no effective solution to such a problem. Thus, it is almost impossible to track many moments with SQMOM.

In this work, the basic idea of STFT is applied to the aggregation and breakage processes, and a novel method, namely, the local fixed pivot quadrature method of moments (LFPQMOM) is proposed, with which moments can be tracked exactly. Meanwhile, the particle size distribution is constructed with high accuracy without introducing further assumption. This work is arranged as follows: a general form for population balance equation is presented, followed by the sectional representation for PBE including aggregation and breakage. Then, the local fixed pivot quadrature method of moments along with the treatment of the aggregation and breakage process is given. Finally, several benchmark test cases including pure aggregation, pure breakage, and aggregation and breakage combined processes are presented to validate the new method.

2. Population Balance Equation

Population balance equation is the continuous form of the number density function. In any disperse system, volume based PBE can be written as:

$$\frac{\partial f(V;\mathbf{x},t)}{\partial t}+\frac{\partial <u_i{>}_{V}f(V;\mathbf{x},t)}{\partial x_i}=S(V;t),$$

(1)

where f(V; x, t) is the number density function, describing the particle number distribution in property space referred to as internal coordinates, location, and time space referred as external coordinates, V is the particle volume, x is space location, t is time, and <u_i>_V is the particle velocity in the ith direction conditioned on volume V. S(V; t) describes the micro-behaviors of the particle swarm. In this work, we focus our attention on the solution method for the population balance equation for tracking particle size distribution and its moment due to aggregation and breakage at a given location x. To shorten the symbolic expression, space coordinate x is neglected in the discussion below.

A general description of aggregation in any particle system can be expressed as Equation (2), but with a different aggregation kernel.

$$S(V)=\frac{1}{2}{\displaystyle {\int }_0^V\beta (V-\upsilon ,\upsilon )f(V-\upsilon ;t)f(\upsilon ;t)\mathrm{d}\upsilon }-f(V;t){\displaystyle {\int }_0^{\infty }\beta (V,\upsilon )f(\upsilon ;t)\mathrm{d}\upsilon },$$

(2)

where β(V, υ) is the aggregation kernel, the first term on the right hand side is the rate of birth of the particles with volume V due to aggregation of smaller particles, whereas the second is the rate of death of the particles with volume V due to their aggregation with other particles. The general description for breakage in any particle system can be expressed as follows, but with different breakage kernel or daughter distribution.

$$S(V)={\displaystyle {\int }_V^{\infty }a(\upsilon )b(V|\upsilon )f(\upsilon ;t)\mathrm{d}\upsilon }-a(V)f(V;t),$$

(3)

where a(V) is the breakage kernel and b(V| υ) is the daughter distribution. The first term on the right-hand side is the rate of birth of particles with volume V due to the breakup of larger particles, whereas the second is the rate of death of particles with volume V due to their breakup.

2.1. Sectional Representation for Population Balance Equation

Multiplying Equation (1) by a property function ψ(V) and integrating it over arbitrary domain [a, b] yields

$${\displaystyle {\int }_a^b\psi (V)\left[\frac{\partial f(V;t)}{\partial t}+\frac{\partial <u_i{>}_{V}f(V;t)}{\partial x_i}\right]\mathrm{d}V}={\displaystyle {\int }_a^b\psi (V)S(V;t)\mathrm{d}V},$$

(4)

Equation (4) is similar to, but not a weak form of, the solution of a general differential equation. A weak form is usually used to construct the algebraic relation of the original variable (particle size distribution f(V;t) in this work) through the integration relation of a differential equation over a cell (a section in this work). However, Equation (4) is used to construct the equations for the moments, which generally represent a physical quantity (total particle number, total particle volume, for instance)

Let the following relation be:

$$m(t)={\displaystyle {\int }_a^b\psi (V)f(V;t)\mathrm{d}V},$$

(5)

where m(t) can be regarded as the local moments of the number density function over the domain of [a, b]. ψ(V) is the property function, and, in the traditional moment methods, monomials are usually adopted. After some simple manipulations, Equation (4) can be written as:

$$\frac{\partial m(t)}{\partial t}+\frac{\partial \left(<u_i>m(t)\right)}{\partial x_i}={\displaystyle {\int }_a^b\psi (V)S(V;t)\mathrm{d}V},$$

(6)

where <u_i> is the average velocity. The term on the right-hand side is the moment transformation for the micro-behaviors, aggregation, and breakage, for instance, influencing the evolution of m(x, t) in state space. Before deriving the actual expression of the source term, we introduce a switch function with a Boolean parameter.

$$\theta (s)=\{\begin{array}{c}1\hspace{1em}if\hspace{0.17em}s\hspace{0.17em}is\hspace{0.17em}true\\ 0\hspace{1em}if\hspace{0.17em}s\hspace{0.17em}is\hspace{0.17em}false\end{array},$$

(7)

2.2. Aggregation

Here, we want to derive a general conservation equation for m(t). The contribution to the variation of the property m(t) in domain [a, b], due to aggregation, can be divided into six parts, as illustrated in Figure 1. In the figure, three kinds of particles can be made: A, particles with volume smaller than a; B, particles with volume between a and b; and C, particles with volume larger than b.

1. A + A → B: Particles with volumes lower than a aggregate to form particles with volume in [a, b). The net rate for the variation in m(t) due to this contribution is given by:

$$S_1^a(t)=\frac{1}{2}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^a\psi (V+\upsilon )\theta (a\le V+\upsilon <b)\beta (V;\upsilon })f(V;t)f(\upsilon ;t)\mathrm{d}V\mathrm{d}\upsilon },$$

(8)

2. A + B → B: Particles with volume lower than a aggregate with particles with volume in [a, b) to form particles smaller than b. This contribution is given by:

$$S_2^a(t)={\displaystyle {\int }_0^a{\displaystyle {\int }_a^b\left(\psi (V+\upsilon )-\psi (V)\right)\theta (V+\upsilon <b)\beta (V;\upsilon })f(V;t)f(\upsilon ;t)\mathrm{d}V\mathrm{d}\upsilon },$$

(9)

3. A + B → C: Particles with volume lower than a aggregate with particles with volume in [a, b) to form particles larger than b. This will cause a net loss for m(t) which can be written:

$$S_3^a(t)=-{\displaystyle {\int }_0^a{\displaystyle {\int }_a^b\psi (V)\theta (V+\upsilon \ge b)\beta (V;\upsilon })f(V;t)f(\upsilon ;t)\mathrm{d}V\mathrm{d}\upsilon },$$

(10)

4. B + B → B: Particles with volume in [a, b) aggregate to form particles in the same range. This contribution can be written as:

$$S_4^a(t)=\frac{1}{2}{\displaystyle {\int }_a^b{\displaystyle {\int }_a^b\left(\psi (V+\upsilon )-\psi (V)-\psi (\upsilon )\right)\theta (V+\upsilon <b)\beta (V;\upsilon })f(V;t)f(\upsilon ;t)\mathrm{d}V\mathrm{d}\upsilon },$$

(11)

5. B + B → C: Particles with volume in [a, b) aggregate to form particles larger than b and this contribution is given by:

$$S_5^a(t)=-\frac{1}{2}{\displaystyle {\int }_a^b{\displaystyle {\int }_a^b\left(\psi (V)+\psi (\upsilon )\right)\theta (V+\upsilon \ge b)\beta (V;\upsilon })f(V;t)f(\upsilon ;t)\mathrm{d}V\mathrm{d}\upsilon },$$

(12)

6. B + C → C: Particles from the domain [a, b) aggregate with particles larger than b to form particles larger than b, which causes a net loss for m(t) in [a, b), and can be written as:

$$S_6^a(t)=-{\displaystyle {\int }_a^b{\displaystyle {\int }_b^{\infty }\psi (V)\beta (V,\upsilon })f(V;t)f(\upsilon ;t)\mathrm{d}V\mathrm{d}\upsilon },$$

(13)

The net rate for the variation in m(t) over the domain [a, b) due to aggregation is the summation of Equations (8)–(13).

2.3. Breakage

For the breakage process, the contribution to the change in m(t) over [a, b) can be divided into two parts, as illustrated in Figure 2. The particle classification is the same as that for aggregation process.

1. C → A + B + C: Particles with volume larger than b break up to form the particles in the domain [a, b). This contribution is given by:

$$S_1^b(t)={\displaystyle {\int }_b^{\infty }b(\upsilon )f(\upsilon ;t){\displaystyle {\int }_a^b\psi (V)b(V|\upsilon )}}\mathrm{d}V\mathrm{d}\upsilon ,$$

(14)

2. B → B + C: Particles in the domain [a, b) break up to form the particles in the same range. This part can be written as:

$$S_2^b(t)=-{\displaystyle {\int }_a^b\psi (V)b(V)f(}V;t)\mathrm{d}V+{\displaystyle {\int }_a^{b}b(\upsilon )f(\upsilon ;t){\displaystyle {\int }_a^{\upsilon }\psi (V)p(V|\upsilon )\mathrm{d}V}\mathrm{d}\upsilon },$$

(15)

The first term on the right-hand side is the removal of particles due to the particle breakage in [a, b) and the second is the addition of particle properties due to the particle birth.

The total contribution to the variation of m(t) due to the particle breakage is the summation of Equations (14) and (15).

3. Local Fixed Pivot Quadrature Method of Moments

To track the moments exactly in this work, the particle size distribution was assumed to have the form.

$$f(V;t)={\displaystyle \sum _{i=0}^{n-1}{\omega }_i(t)\delta [V-V_i]},$$

(16)

where ω_i(t) are the weights; V_i are the quadrature nodes (pivots); n is the number of the quadrature nodes; and δ(V) is the Dirac Delta function.

It should be noted that Fourier transform of the Dirac Delta function is a constant of unity. The Dirac Delta functions in Equation (16) overlap one another in frequency space. Any changes in particle size distribution would be smoothed out; hence, it would not be possible to obtain the particle size distribution if Equation (16) is applied to the whole solution domain. As a result, we could not get the particle size distribution with traditional moment methods such as QMOM or FPQMOM. To retain information of particle size distribution, Equation (16) should be employed in the local domain as performed in short-time Fourier transformation (STFT).

To apply the idea of STFT to the new method, the volume space was divided into [0, υ₁), [υ₁, υ₂), …, [υ_i, υ_i+1), …, [υ_N−1, ∞), N sections. Note that the section in the present work was equivalent to the cell (or bin) in the discretization methods and resembled the window in STFT. υ_i are the abscissas of the section ends, i.e., nodes. In this work, monomials were adopted as property function, i.e., ψ(V) = V^k (k = 0, 1, 2, …, n − 1). In the ith section, the moment equation (Equation (5)) becomes

$$m_k^{(i)}(t)={\displaystyle {\int }_{{\upsilon }_i}^{{\upsilon }_{i+1}}f(V;t)V^k\mathrm{d}V},$$

(17)

Introducing Equation (16) into the ith section yields

$$m_k^{(i)}(t)={\displaystyle \sum _{j=0}^{n_i-1}{\omega }_j^{(i)}(t)}{V}_j^{(i)k},$$

(18)

where $m_k^{(i)}$ is the kth moment of section i; ${\omega }_j^{(i)}$ are the weights of the pivots in section i; $V_j^{(i)}$ is the jth pivot in section i and has to be specified beforehand; and n_i is the number of the pivots in section i.

Equation (18) could be written in matrix form if k = 0, …, n_i − 1 are adopted

$$\left[\begin{array}{cccc}1& 1& \cdots & 1\\ V_0^{(i)}& V_1^{(i)}& \cdots & V_{n_i-1}^{(i)}\\ ⋮& ⋮& \ddots & ⋮\\ V_0^{(i)n_i-1}& V_1^{(i)n_i-1}& \cdots & V_{n_i-1}^{(i)n_i-1}\end{array}\right]\left[\begin{array}{c}{\omega }_0^{(i)}\\ {\omega }_1^{(i)}\\ ⋮\\ {\omega }_{n_i-1}^{(i)}\end{array}\right]=\left[\begin{array}{c}{m}_0^{(i)}\\ m_1^{(i)}\\ ⋮\\ m_{n_i-1}^{(i)}\end{array}\right],$$

(19)

Equation (19) is a vandermonde linear system, the coefficient matrix of which is rather ill-conditioned. Direct inversion of such a coefficient matrix or solving it using traditional numerical methods would inevitably introduce errors. We must resort to the special algorithm for such an equation set, and one such algorithm was recommended by Golub [18].

In theory, the locations of pivots in section i can be specified arbitrarily within [υ_i, υ_i+1) as long as each location has a unique value. It is worth pointing out that the main concern here is the numerical accuracy of the quadrature over [υ_i, υ_i+1) or [υ_n−1, ∞). It is well known that numerical integration of arbitrary smooth functions over a bounded interval is most accurate when carried out with quadrature points chosen as zeros of orthogonal polynomials. Several classical orthogonal polynomials are available [19] such as Gauss-Legendre over (−1, 1) with weighted function W(x) = 1, Gauss-Chebyshev over (−1, 1) with weighted function W(x) = (1 − x²)^−1/2, Gauss-Laguerre over (0, ∞) with weighted function W(x) = x^α e^−x (α > −1), and Gauss-Jacobi over (−1, 1) with weighted function (1 − x)^α(1 + x)^β. The parameters, α in Gauss-Laguerre and α and β in Gauss-Jacobi, controlling the actual form of the polynomials can be used to optimize the distribution of the particle ensembles according to the actual distribution. Note that the orthogonal domains of Gauss-Legendre, Gauss-Chebyshev, and Gauss-Jocobi are not [υ_i, υ_i+1). Thus, a transformation is required before these polynomials can be used in the simulation. The following relation can serve this purpose if υ_i+1 ≠ ∞

$$s=\frac{{\upsilon }_{i+1}-{\upsilon }_i}{2}t+\frac{{\upsilon }_{i+1}+{\upsilon }_i}{2},$$

(20)

Equation (20) is used to transform the pivots t in [−1, 1] to s in [υ_i, υ_i+1). The zeros of the orthogonal polynomial of Gauss-Laguerre with α = 0, are taken as the locations of pivots after transforming t in [0, ∞) to s in [υ_N−1, ∞) using a simple linear relation s = t + υ_N−1.

After transformation using this method, the Gauss-Chebyshev orthogonal polynomial is applied to the first n − 1 sections and Gauss-Laguerre to the last section. An exponential form is adopted when discretizing the property domain, i.e.,

$${\upsilon }_i=s^{i-1}{\upsilon }_1,$$

(21)

with υ₀ = 0 and υ_N = ∞. s is the scale factor.

3.1. Aggregation

For pure aggregation process, there are six terms (Equations (8)–(13)) to contribute the variation for the moments m_k in each section. In this work, monomials were adopted as the property functions, i.e., ψ(V) = V^k (k = 0, 1, 2, …, n − 1). In the lth (l = 0, 1, …, N − 1) section, after introducing the approximation of Equation (16) and some manipulations, Equation (8) yields

$$S_1^a(t)=\frac{1}{2}{\displaystyle \sum _{i=0}^{l-1}{\displaystyle \sum _{j=0}^{l-1}{\displaystyle \sum _{q=0}^{n_j}{\displaystyle \sum _{p=0}^{n_i}{\omega }_{jq}{\omega }_{ip}}\theta ({\upsilon }_l<u_{jq}+u_{ip}<{\upsilon }_{l+1}){(u_{jq}+u_{ip})}^k\beta (u_{jq},u_{ip})}}},$$

(22)

Equation (9) yields

$$S_2^a(t)={\displaystyle \sum _{i=0}^{l-1}{\displaystyle \sum _{q=0}^{n_l}{\displaystyle \sum _{p=0}^{n_i}{\omega }_{lq}{\omega }_{ip}}\theta (u_{lq}+u_{ip}<{\upsilon }_{l+1})[{(u_{lq}+u_{ip})}^k-u_{lq}^k]\beta (u_{lq},u_{ip})}},$$

(23)

Equation (10) yields

$$S_3^a(t)=-{\displaystyle \sum _{i=0}^{l-1}{\displaystyle \sum _{q=0}^{n_l}{\displaystyle \sum _{p=0}^{n_i}{\omega }_{lq}{\omega }_{ip}}\theta (u_{lq}+u_{ip}>{\upsilon }_{l+1})u_{lq}^k\beta (u_{lq},u_{ip})}},$$

(24)

Equation (11) yields

$$S_4^a(t)=\frac{1}{2}{\displaystyle \sum _{p=0}^{n_l}{\omega }_{lp}}{\displaystyle \sum _{q=0}^{n_l}{\omega }_{lq}}\beta (u_{lp},u_{lq})\theta (u_{lp}+u_{lq}<{\upsilon }_{l+1})[{(u_{lp}+u_{lq})}^k-u_{lp}^k-u_{lq}^k],$$

(25)

Equation (12) yields

$$S_5^a(t)=-\frac{1}{2}{\displaystyle \sum _{p=0}^{n_l}{\omega }_{lp}}{\displaystyle \sum _{q=0}^{n_l}{\omega }_{lq}}\beta (u_{lp},u_{lq})\theta (u_{lp}+u_{lq}>{\upsilon }_{l+1})[u_{lp}^k+u_{lq}^k],$$

(26)

and Equation (13) yields

$$S_6^a(t)=-{\displaystyle \sum _{i=l+1}^{N-1}{\displaystyle \sum _{q=0}^{n_l}{\omega }_{lq}}{\displaystyle \sum _{p=0}^{n_i}{\omega }_{ip}}{u}_{lq}^k\beta (u_{lq},u_{ip})},$$

(27)

where N is the total number of the sections; n_i is the count of the pivots in section i; u_ij are the abscissas of jth pivot in section i; and ω_ij is the weight of jth pivot in section i.

The summation of Equations (22)–(27) can be abbreviated

$$\begin{array}{l}{S}^a(t)=\frac{1}{2}{\displaystyle \sum _{i=0}^l{\displaystyle \sum _{j=0}^l{\displaystyle \sum _{q=0}^{n_j}{\displaystyle \sum _{p=0}^{n_i}{\omega }_{jq}{\omega }_{ip}}\theta ({\upsilon }_l\le u_{jq}+u_{ip}<{\upsilon }_{l+1}){(u_{jq}+u_{ip})}^k\beta (u_{jq},u_{ip})}}}\\ \hspace{1em}\hspace{1em}\hspace{1em}-{\displaystyle \sum _{i=0}^{N-1}{\displaystyle \sum _{q=0}^{n_l}{\displaystyle \sum _{p=0}^{n_i}{\omega }_{lq}{\omega }_{ip}}\beta (u_{lq},u_{ip})u_{lq}^k}}\end{array}$$

(28)

3.2. Breakage

Two terms contribute the variation of the moments m_k(t) in the pure breakage process. In the lth (l = 0, 1, …, N − 1) section, after introducing the approximation of Equation (16) and some manipulations, Equation (14) yields

$$S_1^b(t)={\displaystyle \sum _{i=l+1}^{N-1}{\displaystyle \sum _{p=0}^{n_i}{\omega }_{ip}b(u_{ip})}{\displaystyle {\int }_{{\upsilon }_l}^{{\upsilon }_{l+1}}{\upsilon }^{k}p(\upsilon |u_{ip})\mathrm{d}\upsilon }},$$

(29)

and Equation (15) yields

$$S_2^b(t)={\displaystyle \sum _{p=0}^{n_l}{\omega }_{lp}}b(u_{lp}){\displaystyle {\int }_{{\upsilon }_l}^{u_{lp}}{\upsilon }^{k}p(\upsilon |u_{lp})\mathrm{d}\upsilon }-{\displaystyle \sum _{p=0}^{n_l}{\omega }_{lp}}b(u_{lp})u_{lp}^k,$$

(30)

The total contribution to the variation of m_k(t) due to the breakage of particles is the summation of Equations (29) and (30).

$$\begin{array}{l}{S}^b(t)={\displaystyle \sum _{i=l+1}^{N-1}{\displaystyle \sum _{p=0}^{n_i}{\omega }_{ip}b(u_{ip})}{\displaystyle {\int }_{{\upsilon }_l}^{{\upsilon }_{l+1}}{\upsilon }^{k}p(\upsilon |u_{ip})\mathrm{d}\upsilon }}-{\displaystyle \sum _{p=0}^{n_l}{\omega }_{lp}}b(u_{lp})u_{lp}^k\\ \hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{p=0}^{n_l}{\omega }_{lp}}b(u_{lp}){\displaystyle {\int }_{{\upsilon }_l}^{u_{lp}}{\upsilon }^{k}p(\upsilon |u_{lp})\mathrm{d}\upsilon }\end{array}$$

(31)

3.3. Reconstruction of the Moments

In the new method, the moments in the local domain (section) are tracked. The moments of the particle size distribution are the summation of the moments in all the sections, i.e.,

$$m_k(t)={\displaystyle \sum _{i=0}^{N-1}{m}_k^{(i)}(t)},$$

(32)

3.4. Reconstruction of the Particle Size Distribution

The average particle number density in the ith section can be written as

$$\overline{f}{(V;t)}_i=\frac{m_0^{(i)}(t)}{{\upsilon }_{i+1}-{\upsilon }_i},$$

(33)

The average particle size is defined as the weighted average of particles in section i

$$\overline{V_i}=\frac{{\displaystyle \sum _{j=0}^{n_i-1}{\omega }_j^{(i)}{V}_j^{(i)}}}{{\displaystyle \sum _{j=0}^{n_i-1}{\omega }_j^{(i)}}}=\frac{m_1^{(i)}(t)}{m_0^{(i)}(t)},$$

(34)

It is worth pointing out that, in traditional methods of moment (MOMs) such as quadrature method of moment (QMOM) [7] or fixed pivot quadrature method of moment (FPQMOM) [12], the assumption of continuous summation of the Dirac Delta function is adopted in the overall domain, thus the particle size distribution is lost in the simulation. With traditional MOMs, particle size distribution cannot be obtained. However, with the new method detailed here, the assumption of the continuous summation of the Dirac Delta function is adopted in the subsection, detailed information in the local domain is retained in the simulation, and, with this information, the particle size distribution can be reconstructed without any additional assumption. If only one section is adopted in the simulation, the method here is the same as FPQMOM.

4. Test Cases

In all test cases presented here, the volume space was divided into 40 sections using Equation (21) with υ₁ = 0.01 and s = 1.2. In each section the six pivots were specified (i.e., six moments were tracked). A time step of 0.01 s was adopted. All simulations in this work were carried out on a desktop computer with Intel Core i5 CPU (3.1 GHz) and 8 G memory.

4.1. Pure Aggregation

Five test cases were carried out to validate the new method for the pure aggregation process, in the first of which a constant aggregation kernel was adopted, having the form

$$\beta (V,\upsilon )=C_0,$$

(35)

with C₀ = 1, together with exponential initial particle size distribution

$$f(V)=\frac{N_0}{V_0}\exp (-\frac{V}{V_0}),$$

(36)

with N₀ = 1, V₀ = 1. In this test case, an analytical solution was available [20]

$$f(V,T)=\frac{{(1-T)}^2{N}_0}{V_0}\exp (-(1-T)\frac{V}{V_0}),$$

(37)

with T = C₀N₀t/(2 + C₀N₀t). The moments could be evaluated analytically. Figure 3 depicts the evolution of the percentage errors for the first six moments with: (A) LFPQMOM; (B) QMOM; and (C) FPQMOM. In the figure, the evolution for the relative errors for all the moments were similar in both methods, but a little larger than FPQMOM. Errors increased rapidly initially and leveled off after 2 s. The maximum percentage error with QMOM was 1.12 × 10⁻⁹%, whereas with the new method, it was 1.10 × 10⁻⁹%, a small decrease. FPQMOM had the highest accuracy. Higher rank moments had larger percentage errors with both methods. It is worth pointing out that, with LFPQMOM, the relative errors for m₁ were always zero within the simulation time, but with QMOM, they were lower than 10⁻¹²%, but larger than zero. With the new method, m₁ could be predicted without any error. In the figure, percentage errors for m₁ were not included due to their low values.

Figure 4 gives the particle size distribution predicted by LFPQMOM at 0, 2, 5, and 10 s. In the figure, the numerical PSDs agreed exactly with the analytical PSDs. The number density approached zero when the dimensionless volume was larger than 10.

The computational time was about 28,360, 1795, and 625 ms for LFPQMOM, QMOM, and FPQMOM, respectively. The computational load of LFPQMOM was much larger than QMOM and FPQMOM, since, in LFPQMOM, 40 sections were adopted in this work and, in each section, six moments were tracked. The total number of variables included in the simulation with LFPQMOM was 240, whereas, in QMOM and FPQMOM, only six variables were included. The computational time for each variable on average was 118, 299, and 104 ms for LFPQMOM, QMOM, and FPQMOM, respectively. From the view of a one-variable computational load, LFPQMOM and FPQMOM were more efficient than QMOM because of the vandermode linear system solver used.

In the second test case, a constant aggregation kernel was used, i.e., Equation (35) with C₀ = 1, but with Gaussian initial particle size distribution, having the form

$$f(V)=\frac{N_{0}V}{V_0^2}\exp (-\frac{V}{V_0}),$$

(38)

with N₀ = 1 and V₀ = 1. In this test case, analytical particle size distribution is given by [20]

$$f(V,T)=\frac{{(1-T)}^2{N}_0}{\sqrt{T}{V}_0}\exp (-\frac{1}{V_0}V)sinh(\frac{\sqrt{T}}{V_0}V),$$

(39)

with T = C₀N₀t/(2 + C₀N₀t). The moments can be evaluated analytically.

Figure 5 illustrates the relative errors for the first six moments with: (A) LFPQMOM; (B) QMOM; and (C) FPQMOM. With LFPQMOM, the relative errors rose rapidly initially, and became flat after 2 s. The relative errors for m₂ after 6 s and m₄ and m₅ after 9 s began to rise, but for m₀ and m₃ they began to fluctuate after 7 s. With QMOM, the relative errors for all the moments rose sharply initially, and leveled off after 2 s, the evolution of which was very similar to that in the previous case. With the new method, the prediction for the moments of m₂, m₄, and m₅ was poor relative to QMOM, but was lower 10⁻⁶% within the simulation time, which could still satisfy industrial needs. With FPQMOM, the numerical prediction was the best. In this test case, m₃ could be predicted exactly without any error within the simulation time with the new method. The percentage errors for m₃ were not included in the figure due to their low values with both methods. Figure 6 gives the particle size distribution predicted by LFPQMOM at 0, 2, 5, and 10 s, along with analytical solutions for comparison. In the figure, the numerical PSDs with our new method agreed well with the analytical PSDs. Again, the new method made an excellent prediction for both moments and PSD. The computational time was about 71,766, 5024, and 1709 for LFPQMOM, QMOM, and FPQMOM, respectively, in this test.

One should note that, with LFPQMOM, we did not track the moments of the PSD over the whole domain, but the moments in each section. As a result, to compute relative error with LFPQMOM, we had to resort to Equation (32) to first compute the moments of the PSD. Thus, the numerical errors of each moment presented in the figure for LFPQMOM are the compound errors, but not for one variable. For instance, to calculate errors of m₃, we calculated the m₃ with m_3i of the whole PSD using Equation (32) first, and then compared it with its analytical counterpart where m_3i is the tracked variable not m₃. In each section, a fixed pivot quadrature method of moment (FPQMOM) was used. FPQMOM was more accurate than QMOM, as can be seen in the figure. As a result, the relative errors of the moments with FPQMOM do not deviate from QMOM or LFPQMOM too much. If the relative error for a tracked variable is too small, the values of the relative errors can be random. The compound relative errors can also be random. The fluctuations in the errors of m₀ in Figure 5A may be related to this.

In the third test case, a sum aggregation kernel was used

$$\beta (V,\upsilon )=C_0(V+\upsilon ),$$

(40)

with C₀ = 0.01, together with the exponential initial distribution, i.e., Equation (36) with N₀ = 1 and V₀ = 1. Analytical solution was given by [20]

$$f(x,\tau )=\frac{(1-\tau )\exp (-x(1+\tau ))}{x{\tau }^{1/2}}{I}_1(2x{\tau }^{1/2}),$$

(41)

with x = V/V₀ and τ = 1 − exp(−C₀N₀V₀t). I₁(x) is the modified Bessel function of the first kind. Moments can be evaluated using the numerical integration method.

Figure 7 shows the evolution of percentage errors for the first six moments with: (A) LFPQMOM; (B) QMOM; and (C) FPQMOM. In the figure, within the simulation time, the relative errors for all the moments were lower than 10⁻⁵%. It should be noted that the relative errors for m₁ were not lower than 10⁻¹²%, even though the conservation of the mass was considered in the new method. This is because no analytical solutions for the moments were available. The numerical errors for m₁ were not caused by the new method or QMOM, but came from the evaluation of the moments through numerical integration. The sharp drops in the relative errors were related with the numerical errors of the evaluation of the moments. Figure 8 illustrates the particle size distribution at 0, 2, 10, and 30 s. It can be observed that the numerical predictions with the new method were in excellent agreement with the analytical counterparts. Number density with dimensionless volume larger than 10 approached zero. The computational time was about 86,314, 5570, and 1948 ms for LFPQMOM, QMOM, and FPQMOM, respectively.

In the fourth test case, a product aggregation kernel was adopted, having the form

$$\beta (V,\upsilon )=C_{0}V\upsilon ,$$

(42)

with C₀ = 0.01, together with the exponential initial PSD, i.e., Equation (36) with N₀ = 1 and V₀ = 1. Analytical solution for this test case was given by [20]

$$f(x,T)=\frac{N_0}{V_0}\exp (-x(1+T)){\displaystyle \sum _{k=0}^{\infty }\frac{T^k{x}^{3k}}{(k+1)!\mathsf{\Gamma }(2k+2)}},$$

(43)

where $T=C_0{N}_0{V}_0^{2}t$, x = V/V₀ and Γ (x) is the gamma function.

Note that the product kernel is a gelling kernel, and the second and higher moments will become infinite within finite time [21]. Most existing moment methods cannot predict the moments exactly, especially near the gelling point. The analytical solutions for moments m₀, m₁, and m₂ are given in

Table 1. Analytical solution for the first three moments for the aggregation process with product kernel.

Moments	Analytical Solution
m0(t)	$m_0(t)=m_0(0)-\frac{1}{2}{C}_0{m}_1^2(0)t,\hspace{0.17em}\hspace{0.17em}0<t<t_{gel}$
m1(t)	$m_1(t)=m_1(0),\hspace{0.17em}\hspace{0.17em}0<t<t_{gel}$
m2(t)	$m_2(t)=\frac{m_2(0)}{1-C_0{m}_2(0)t},\hspace{0.17em}\hspace{0.17em}0<t<t_{gel}$

. Based on the solution for m₂, the actual value for t_gel can be calculated

$$t_{gel}=\frac{1}{C_0{m}_2(0)},$$

(44)

The gelling time t_gel is 50 s in this test case. To predict the moments as exactly as possible, especially near the gelling point, two nodes were adopted for QMOM and three pivots in each section for LFPQMOM. The percentage errors of the first three moments are shown in Figure 9A–C for LFPQMOM, QMOM, and FPQMOM, respectively. An interesting feature with the new method is that the simulation can be carried out for as long as possible before the gelling point. The maximum time that can be used as the simulation time is t_gel − Δt. Δt is the time step. In this test case, a time step of 0.01 s was adopted, thus any time of t ≤ 49.99 s could be adopted as the simulation time. With QMOM, only an even number of moments could be tracked, four moments in this test case. We could only carry out the simulation before the time of 49.6 s when m₃ diverged. After that time, QMOM was not appropriate, and the new method could also make a good prediction, with a relative error smaller than 10⁻⁵%, as can be seen in Figure 9A. Relative errors with FPQMOM were similar to those with LFPQMOM. The relative errors for m₁ with all methods here were lower than 10⁻¹³%, hence were not included in the figures. Figure 10 shows a good agreement of numerical PSD predicted by the new method at 0, 10, 20, and 50 s with the analytical solution. Number density approached zero with a dimensionless volume larger than 10. The computational time was about 143,962, 8957, and 3199 ms for LFPQMOM, QMOM, and FPQMOM, respectively, in this test.

In the final test case for pure aggregation, a product kernel was used, i.e., Equation (42) with C₀ = 0.01, but with a Gaussian initial PSD, i.e., Equation (38) with N₀ = 1 and V₀ = 1. Analytical solution for PSD is given by [21]

$$f(x,T)=\frac{N_0}{V_0}\exp (-x(2+T)){\displaystyle \sum _{k=0}^{\infty }\frac{T^k{x}^{4k+1}{2}^{3k+3}}{(k+1)!\mathsf{\Gamma }(3k+3)}},$$

(45)

where $T=C_0{N}_0{V}_0^{2}t$, x = V/V₀ and Γ (x) is the gamma function. The analytical moments are listed in Table 1. The gelling time t_gel was 16.67 s. Similar to the previous case, two nodes were adopted with QMOM and three pivots in each section with LFPQMOM to track the moments exactly. Figure 11 depicts the evolution of the relative errors for the first six moments predicted by: (A) LFPQMOM; (B) QMOM; and (C) FPQMOM. The gelling time was 16.67 and a simulation time of 16.67 − Δt (=16.66) could be carried out with the new method, when m₃ diverged. A maximum simulation time of 16.65 s could only be adopted with QMOM. Within a simulation of 16.65 s, the percentage errors with both methods were comparable, which were lower than 10⁻⁶%, but, after that time, QMOM was not appropriate for this test case. Similar to the previous case, the relative errors for m₂ near the gelling point rose sharply. The relative errors for m₁ were not included in the figure for their low values. Again, numerical errors with FPQMOM were similar to LFPQMOM. Figure 12 gives the particle size distributions at 0, 2, 10, and 16.66 s predicted by LFPQMOM along with their analytical counterparts for comparison. The numerical PSDs were in excellent agreement with the analytical PSDs on all the time points sampled. The computational time was about 46,042, 2818, and 996 ms for LFPQMOM, QMOM, and FPQMOM, respectively, in this test.

4.2. Pure Breakage

Three test cases were conducted for the pure breakage process. In the first test case, constant breakage kernel was adopted, i.e.,

$$a(V)=C_0,$$

(46)

with C₀ = 0.1. Uniform distribution was used as the daughter distribution

$$b(V|\upsilon )=2/V,$$

(47)

Using Equation (36) with N₀ = 1 and V₀ = 1 as the initial distribution, the analytical solution for the moments was given by [22]

$$m_k=m_k(0)\exp [\frac{1-k}{1+k}{C}_{0}t],$$

(48)

Unfortunately, analytical PSD was not available in this test case. Figure 13 shows the evolution of the relative errors for the first six moments predicted by: (A) LFPQMOM; (B) QMOM; and (C) FPQMOM. It was observed that the evolutions of the percentage errors of LFPQMOM and QMOM were very similar. The relative errors rose sharply initially and leveled off after 15 s. Lower rank moments had bigger numerical errors and within the simulation time, the relative errors for all the moments were lower than 2 × 10⁻⁹%. The percentage errors for m₃ were not included in the figure due to their low values. Similar to the previous cases, m₃ could be predicted without any error with the new method. Numerical predictions with FPQMOM were the best in this test. Figure 14 depicts the particle size distribution at 0, 2, 10, and 50 s. Due to lack of analytical PSD, a comparison of the numerical PSDs with their analytical counterparts was impossible, and thus, in this test case, only numerical PSDs are given at 2, 10, and 50 s. The computational time was about 14,850, 1037, and 345 ms for LFPQMOM, QMOM, and FPQMOM, respectively, in this test.

In the second test case for breakage process, linear breakage kernel was used, i.e.,

$$a(V)=a_{0}V,$$

(49)

with a₀ = 0.1, together with a uniform daughter distribution, i.e., Equation (47), exponential initial distribution, i.e., Equation (36) with N₀ = 1 and V₀ = 1. Analytical solution was given by [22]

$$f(V;t)=\frac{N_0}{V_0}{(1+a_0{V}_{0}t)}^2\exp \left(-\frac{V}{V_0}(1+a_0{V}_{0}t)\right),$$

(50)

Moments can be evaluated by numerical integration or direct integration analytically. Percentage errors for the first six moments—predicted by: (A) LFPQMOM; (B) QMOM; and (C) FPQMOM—are shown in Figure 15. In Figure 15B, the relative errors for the moments predicted by QMOM rose sharply at incipient time, especially for larger rank moments, m₄ and m₅ for instance. However, with the new method, the relative errors rose relatively slowly at incipient time, as can be seen in Figure 15A. With the simulation time, the maximum error for all the moments with FPQMOM was 1.88% whereas with QMOM, the maximum was 3.66%. With the new method, the numerical accuracy was improved relative to QMOM in this test case. The relative errors for m₀ were comparable. Again, with the new method, m₁ was predicted without any error, but not with QMOM. Relative errors for m₁ were not included in the figure. FPQMOM was again the best in numerical prediction, and the relative errors for all the moments were less than 0.1% within the simulation time. Figure 16 gives the particle size distribution at 0, 2, 5, and 10 s, from which we observed that the numerical predictions with the new method agreed very well with the analytical PSDs at the time sampled. The computational time was about 3020, 190, and 66 ms for LFPQMOM, QMOM, and FPQMOM, respectively, in this test.

In the third test case for pure breakage, the power kernel was used as the breakage kernel having the form

$$a(V)=C_0{V}^2,$$

(51)

with C₀ = 0.01. Uniform distribution was used for the daughter distribution, i.e., Equation (47), together with the exponential initial distribution, i.e., Equation (36) with N₀ = 1 and V₀ = 1. Analytical solution was given by [22]

$$f(V;t)=\frac{N_0}{V_0}[1+2C_{0}t(V_0^2+V_{0}V)]\exp (-\frac{V}{V_0}-C_{0}t{V}^2),$$

(52)

The numerical method was used to evaluate the moments. Evolutions of the relative errors for the first six moments are shown in Figure 17A,B for LFPQMOM and QMOM, respectively. In the figure, it was observed that, with the new method, all the relative errors were lower than 0.3% before the time of 20 s, and then fluctuated for m₄ and m₅, and rose for m₀–m₃. Within the simulation time, all the relative errors were lower than 2%. However, with QMOM, all the relative errors for the moments except m₀ and m₁ rose up to 0.1 within 5 s initially. Within the simulation time, the maximum percentage error for all the moments was 13.8%. With the new method, the numerical accuracy was improved in this test case. Again, LFPQMOM made an exact prediction for m₁, similar to the previous cases. The relative errors for m₁ were not included in the figure. Again, LFPQMOM had the smallest errors. Similar to the third case for the aggregation process, the sharp drops in the relative errors were related to the numerical evaluation of the moments. Figure 18 depicts the particle size distribution predicted by LFPQMOM at 0, 10, 20, and 50 s along with analytical counterparts for comparison. The numerical predictions exactly agreed with the analytical PSDs at different time points. The computational time was about 15,050, 876, and 338 ms for LFPQMOM, QMOM, and FPQMOM, respectively, in this test.

4.3. Aggregation and Breakage

Two test cases were conducted to validate the new method for the aggregation and breakage combined processes. The resulting source for Equation (6) is the summation of Equations (28) and (31). In the first test case, a constant kernel for aggregation, i.e., Equation (35) with C₀ = 1, a linear kernel for breakage, i.e., Equation (49) with a₀ = 0.1, and uniform daughter distribution, i.e., Equation (47), along with the following initial particle size distribution were adopted

$$f(L)=\frac{N_0^2}{V_0}\exp (-\frac{N_0}{V_0}V),$$

(53)

with N₀ = 1, V₀ = 1. In this test case, analytical PSD was available [23]

$$f(V;t)=\frac{N_0^2}{V_0}{[\mathsf{\Phi }(\tau )]}^2\exp [-\frac{N_0}{V_0}V\mathsf{\Phi }(\tau )],$$

(54)

with Φ(τ) = Φ(∞)[1 + Φ(∞)tanh(Φ(∞)τ/2)]/[Φ(∞) + tanh(Φ(∞)τ/2)], Φ(∞) = (2a₀V₀/C₀)^1/2/N₀, and τ = C₀N₀t. Moments could be evaluated analytically. Evolutions of percentage errors for the moments are shown in Figure 19 with: (A) LFPQMOM; (B) QMOM; and (C) FPQMOM. In the figure, the evolutions for the relative errors were very similar with the three methods. Initially, all the relative errors rose sharply and leveled off after the time of 20 s, except the errors for m₃, which had a continuous drop as the time advanced. With the new method and FPQMOM, there were two fluctuations in the drop of relative errors for m₃, but not with QMOM. Within the simulation time, the maximum relative error with both methods were comparable. In the figure, the relative errors for m₀ and m₁ were not included due to their low values. Figure 20 gives the particle size distributions at 0, 2, 10, and 50 s. It was observed that the numerical PSDs agreed very well with the analytical PSDs. Note that the PSD at 10 s coincided with the PSD at 50 s, proving that a steady state had been reached by the time of 10 s. Note that, with different physical parameters, the steady state can be different. The computational time was about 497,950, 32,275, and 10,374 ms for LFPQMOM, QMOM, and FPQMOM, respectively, in this test.

In the second test case for the aggregation and breakage combined process, a constant aggregation kernel was used, i.e., Equation (35) with C₀ = 1, a linear kernel for breakage, i.e., Equation (49) with a₀ = 0.25, and uniform daughter distribution, i.e., Equation (47), along with the following initial particle size distribution were adopted

$$f(V)=N_0{\left[2\frac{N_0}{V_0}\right]}^{2}V\exp \left(-2\frac{N_0}{V_0}V\right),$$

(55)

with N₀ = 1 and V₀ = 2. When the following relation was satisfied, the analytical solution was available [23]

$$2a_0{V}_0/C_0=N_0^2,$$

(56)

The analytical solution for particle size distribution was given by [23]

$$f(V;t)=\frac{N_0^2}{V_0}{\displaystyle \sum _{i=1}^2\frac{K_1(\tau )+p_i{K}_2(\tau )}{L_2(\tau )+4p_i}\exp \left(p_i\frac{N_0}{V_0}V\right)},$$

(57)

with K₁(τ) = 7 + τ + exp(−τ), K₂(τ) = 2 − 2exp(−τ), L₂(τ) = 9 + τ − exp(−τ), p_1,2 = 0.25(e^−t−τ − 9 ± (d(τ))^0.5), d(τ) = τ² + (10τ − 2τe^−τ) + 25 − 26e^−τ + e^−2τ, τ = C₀N₀t. The moments can be evaluated analytically.

Figure 21 depicts the evolutions of percentage errors for m₂–m₅ with: (A) LFPQMOM; (B) QMOM; and (C) FPQMOM. The figure is very similar to the previous case. Initially, the errors rose sharply and leveled off for m₂, m₄, and m₅, but decreased for m₃ continuously. There was one difference from the previous case: the maximum percentage error was 0.42% with the new method, and 1.00% with QMOM. The numerical accuracy was improved with the new method relative to QMOM, which was similar to FPQMOM. Relative errors for m₀ and m₁ were not included in the figure given their low values. The numerical and the analytical particle size distributions at 0, 0.5, 2, 10, and 20 s are shown in Figure 22. The numerical PSDs agreed well with the analytical PSDs. The PSDs at 10 and 20 s coincided with each other, demonstrating that, by the time of 10 s, a steady state had been reached. Initially, the PSD was a Gaussian distribution, but, as the time advanced, the distributions became exponential, as can be seen in the figure. With the other parameters, the particle size distributions and the steady state could be different. The computational time was about 198,260, 12,048, and 4198 ms for LFPQMOM, QMOM, and FPQMOM, respectively, in this test.

5. Conclusions

In this work, a new numerical method, namely, the local fixed pivot quadrature method of moment (LFPQMOM) was proposed for the aggregation and breakage process. In the new method, particle size distribution in the local section was assumed to be a continuous summation of the Dirac Delta function with fixed pivots as performed in short time Fourier transformation (STFT). This assumption was applied to the sectional representation of aggregation or breakage, and the moments of the particle size distribution in the local domain were tracked successfully. The moments in the whole domain and the particle size distribution were constructed based on the moments in the local section without introducing further assumption.

Simulation of benchmark test cases including pure aggregation, pure breakage, and aggregation breakage combined process demonstrated that: (1) for the pure aggregation process, constant breakage process, and the aggregation and breakage combined process with exponential initial PSD, the accuracy of the prediction for the moments with the new method was comparable to quadrature method of moment (QMOM); (2) for the linear or square breakage process, and aggregation and breakage combined process with Gaussian initial PSD, the numerical accuracy for the moments was improved further; and (3) the new method could make good predictions for PSD in all of the test cases in this work.

A special algorithm was employed to solve the vandermode linear system, with which the coefficient matrix did not need to be constructed and, therefore, the influence of the ill-conditional feature of such a matrix on the numerical accuracy could be avoided. In theory, any number of moments can be tracked. Due to the large number of scalar equations included, the computational expense can be relatively large. Extension of the new algorithm to the solution of the multi-variable population balance equation will be investigated in depth in our future work.