A numerical comparison of the method of moments for the population balance equation

TL;DR

This paper compares various moment closure methods for solving the population balance equation, analyzing their accuracy, realizability, and computational efficiency through diverse numerical examples.

Contribution

It provides a comprehensive numerical comparison of polynomial, maximum entropy, and quadrature moment methods for population balance equations, including implementation insights.

Findings

Quadrature method of moments (QMOM) shows favorable accuracy and efficiency.

Maximum entropy closures offer good realizability properties.

Method of moments accuracy depends on the closure type and problem complexity.

Abstract

We investigate the application of the method of moments approach for the one-dimensional population balance equation. We consider different types of moment closures, namely polynomial (P_N) closures, maximum entropy (M_N) closures and the quadrature method of moments QMOM_N. Realizability issues and implementation details are discussed. The numerical examples range from spatially homogeneous cases to a population balance equation coupled with fluid dynamic equations for a lid-driven cavity test case. A detailed numerical discussion of accuracy, order of the moment method and computational time is given.

Tables

Table 1. Table 1: Parameter choice for the Newton algorithm in the case of the maximum entropy closure.

Parameter	$k_{max}$	$ϵ$	$\chi$	$\tau$	$\left\{r_l\right\}$
Value	$400$	$2^{-52}$	$3/5$	${10}^{-9}$	$\{0,{10}^{-10},{10}^{-8},{10}^{-6},{10}^{-4},{10}^{-2},0.1,0.5,1\}$

Table 2. Table 2: Parameter choice for the volume domain and the initial condition.

Parameter	$v_{min}$	$v_{max}$	${\sigma }_{init}$	${\alpha }_0$	$d_0$	$N_0$	$v_0$
Value	$\frac{1}{6}\pi {0.001}^3$	$\frac{1}{6}\pi {0.009}^3$	$0.1\left(v_{max}-v_{min}\right)$	$0.01$	${\left(\frac{3}{\pi }\left(v_{min}+v_{max}\right)\right)}^{\left(1/3\right)}$	$\frac{6{\alpha }_0}{\pi d_0^3}$	$\frac{\pi }{6}{d}_0^3$

Table 3. Table 3: Parameter choice for the breakage frequency.

Parameter	$C_1$	$C_2$	${\alpha }_d$	$ϵ$	${\rho }_d$	$\sigma$
Value	$0.12$	$0.078$	${\alpha }_0$	$0.004$	$865.6$	$0.0361$

Table 4. Table 4: Parameter choice for the aggregation kernel.

Parameter	$C_{\Omega }$	$k_{{\omega }^{\prime }}$	${\alpha }_d$	$ϵ$	${\rho }_c$	$\sigma$	${\eta }_c$
Value	$41.2$	$1.33\cdot {10}^{10}$	${\alpha }_0$	$0.004$	$1000$	$0.0361$	$0.001$

Equations

$${\partial }_{t}f(t,x$$

$$=$$

$$+\underset{\text{=:B+[f](t,x,v)}}{\underbrace{{\int }_v^{v_{max}}\Gamma (t,x,v^{\prime })\beta (v,v^{\prime })f(t,x,v^{\prime })d{v}^{\prime }}}-\underset{\text{=:B−[f](t,x,v)}}{\underbrace{\Gamma (t,x,v)f(t,x,v)}}.$$

$$\beta (v,v^{\prime })$$

$${\int }_{v_{min}}^{v^{\prime }}\beta (v,v^{\prime })dv$$

$${\int }_{v_{min}}^{v^{\prime }}v\beta (v,v^{\prime })dv$$

$$\omega (t,x,v,v^{\prime })$$

$$\gamma (t,x):=⟨mf(t,x,\cdot )⟩:={\int }_{v_{min}}^{v_{max}}m(v)f(t,x,v)dv,$$

$${\partial }_t\gamma (t,$$

$$=$$

$$⟨m_i(v)f(v)⟩={\gamma }_i,0\le i\le N.$$

$${\text{R}}_N=\{\gamma \in {\mathbb{R}}^{N+1}|\exists f\ge 0,⟨mf⟩=\gamma \}$$

$$f_{P_N}(t,x,v)=\sum _{i=0}^N{a}_i(t,x)m_i(v).$$

$${\gamma }_i(t,x)=\sum _{j=0}^N{a}_j(t,x)⟨m_i(v)m_j(v)⟩,i=0,...,N.$$

$$m_0(v)=1,m_1(v)=\frac{2}{v_{max}-v_{min}}$$

$$m_i(v)=\frac{2i-1}{i}\left(\frac{2}{v_{max}-v_{min}}(v-v_{min})-1\right)$$

$${\gamma }_i(t,x)$$

$$f_{P_N}(t,x,v)=\sum _{i=0}^N{\gamma }_i(t,x)\frac{2i+1}{v_{max}-v_{min}}{m}_i(v).$$

$$\underset{g}{\min }$$

$$⟨mg⟩=\gamma .$$

$$\underset{\alpha \in {\mathbb{R}}^{N+1}}{\min }\left(⟨\eta (g)+{\alpha }^{T}mg⟩-\alpha \gamma \right).$$

$$\hat{\alpha }(\gamma )=\underset{\alpha \in {\mathbb{R}}^{N+1}}{\mathrm{argmin}}\left(⟨{\eta }_{\ast }({\alpha }^{T}m)⟩-\alpha \gamma \right).$$

$$f_{M_N}=G_{\hat{\alpha }}:={\eta }_{\ast }^{\prime }({\hat{\alpha }}^{T}m).$$

$$G_{\alpha }=\exp ({\alpha }^{T}m).$$

$$h(\alpha )=⟨G_{\alpha }⟩-{\alpha }^T\gamma .$$

$$f_{QMO{M}_N}=\sum _{j=1}^n{w}_j\delta (v-v_j),$$

$${\gamma }_i=\sum _{j=1}^n{w}_j{v}_j^i,i=0,...,N.$$

$${\partial }_t(vf(t,x,v))+{\nabla }_x\cdot (v_d(t,x)f(t,x,v))-{\nabla }_x\cdot (D(t,x){\nabla }_{x}f(t,x,v))+{\partial }_v\left(F^{agg}(t,x,v)+F^{brk}(t,x,v)\right)=0,$$

$$F^{agg}(t,x,v)={\int }_{v_{min}}^v{\int }_{v-u}^{v_{max}}u\omega (t,x,u,w)f(t,x,u)f(t,x,w)dwdu,$$

$$F^{brk}(t,x,v)=-{\int }_v^{v_{max}}{\int }_{v_{min}}^{v}u\beta (u,w)\Gamma (t,x,w)f(t,x,w)dudw.$$

$${\Lambda }_i:=[v_{i-\frac{1}{2}},v_{i+\frac{1}{2}}[,i=1,...,n_v\text{ with }{v}_{i+\frac{1}{2}}=v_{min}+i\Delta v.$$

$$f_i(t,x)=\frac{1}{\Delta v}{\int }_{v_{i-\frac{1}{2}}}^{v_{i+\frac{1}{2}}}f(t,x,v)dv.$$

Full text

††journal: arXiv.org

1 Introduction

Population balance equations are widely used in engineering applications, including aerosol physics, high shear granulation, pharmaceutical industries, polymerization and emulsion processes, evaporation and condensation processes in bubble column reactors, bioreactors, turbulent flame reactors and many others, see [44, 43, 45, 13, 19, 24, 24] and references therein. These polydisperse processes are characterized by two phases: one of them is the continuous and the second is a dispersed phase consisting of particles. The particles can take different forms like crystals, drops or bubbles with several possible properties such as volume, chemical composition, porosity and enthalpy. In this work only the volume is considered.

The dynamic evolution of the particle number distribution which is described by the population balance equation (PBE) of the dispersed phase depends not only on the particle-particle interactions, but also on the continuous phase, due to interaction of these particles with the continuous flow field in which they are dispersed. These interactions usually result in the common mechanisms aggregation, breakage, condensation, growth and nucleation. In our work we concentrate on binary aggregation, namely the Smoluchowski aggregation [42] and multiple breakages, since binary breakage is not sufficient for some of these applications. Binary aggregation is the process of merging two particles to a larger particle, whereas in a breakage process, a particle breaks into several smaller fragments. Often included in a typical PBE are spatial transport terms, i.e. advection and diffusion terms. We use the simplest case and assume that the mean particle velocity is the same as those of the fluid. The resulting population balance equations range from integro-differential to partial integro-differential equations of hyperbolic or parabolic type in phase space, in addition to the differential terms for advection and diffusion in the physical space.

The complex structure of the PBE allows for an analytical solution only for some simple breakage and aggregation kernels, see [20, 60] and references therein. Thus, numerical solution methods have been extensively studied. In literature, several different classes of methods are discussed like Monte-Carlo methods, see, for example [35] or finite element techniques, see for example [38] and references therein. Other numerical techniques are the method of successive approximations [43], the methods of classes [43, 30, 51], and the method of weighted residuals [43, 57]. Another approach is to apply finite volume schemes (FVS), which are frequently used for solving conservation laws. In [22, 21], this approach is used for solving aggregation PBEs using a conservative formulation of the equation. Later, the scheme has been applied to solve the combined breakage and aggregation PBE in [31]. In [33], sectional methods have been developed, among them the cell average technique for both uniform and non-uniform grids.

Another class of methods for solving the complex PBE is the method of moments. Instead of solving the PBE in full phase space, moment models reduce the high dimensionality of the original equation by using a finite set of moment equations, which are obtained by taking volume averages with respect to some basis of the volume space. The resulting moment equations require additional information about the unknown particle number distribution that must be approximated via a closure. By the use of such an ansatz function the distribution function can be reconstructed from the available moments. One big issue regarding the method of moments is realizability. Roughly speaking, realizable moments are those moments associated with non-negative distribution functions. The set of realizable moments forms a convex cone in the set of all moment vectors. Several types of closures are available. One that is often used for population balance equations is the quadrature method of moments ($QMOM$), first introduced in the PBE context in [40]. It uses a non-linear Gaussian-like quadrature rule to compute weights and abscissas from the available moments, assuming that the underlying distribution function is a weighted sum of Dirac $\delta-$functions. The computation of the quadrature rule involves solving a non-linear system. There also exist several extensions of $QMOM$, such as $SQMOM$ [8, 12, 9] and $EQMOM$ [59].

The standard moment method in linear radiative transport is the polynomial closure $P_{N}$ [37], where $N$ is the order of the highest-order moment. For multi-variable cell population balance equations, $P_{N}$ closures of different types where discussed in [39]. Using the polynomial closure with respect to an orthogonal basis ensures a non-expensive moment method, at the price of a potentially negative, and thus physically meaningless, number density. Another moment approach is the use of entropy-based closures $M_{N}$ which approximate the full distribution by an ansatz that solves a constrained, convex optimization problem. Using physically relevant entropies, the big advantage of $M_{N}$ is the preservation of many important properties like positivity, entropy dissipation, and hyperbolicity [36]. The disadvantage is the computational cost, especially compared to the polynomial closure, since an optimization problem has to be solved numerically in every point of the space-time grid. Realizability plays a big role for $M_{N}$, since the solving the optimization problem near the realizability boundary is especially expensive or even impossible.

In the present paper we aim at a numerical comparison of the above described closure approaches. The outline of the paper is the following. First, in section 2 we introduce the population balance equation used in this paper. Then we shortly state the finite volume scheme, based on the work in [34]. After that in section 3 we explain the method of moments and the concept of realizability. Also, contained in this section are the different types of moment closures for which we compute the PBE, namely the polynomial closure $P_{N}$, the maximum entropy closure $M_{N}$ and the quadrature method of moments $QMOM_{N}$. This is followed by a description of the numerical methods, supplemented with further implementation details in section 4. The numerical results for the homogeneous case neglecting diffusion and advection are presented in section 5. These results are followed by two-dimensional results on a PBE coupled to fluid dynamic equations for a lid-driven cavity test. In section 6 we summarize our conclusions.

2 Population Balance Equation

We consider a one-dimensional population balance equation (PBE). It describes the time evolution of the particle number distribution function $f(t,x,v)$ under the simultaneous effect of binary aggregation, multiple breakage, diffusion and advection processes [34]. The PBE reads

[TABLE]

The particle number density $f(t,x,v)\geq 0$ depends on the time $t\in\mathbb{R}^{+}$, the space $x\in\mathbb{R}^{n_{d}},~{}n_{d}\geq 1$ and the volume $v\in[v_{min},v_{max}]\subseteq\mathbb{R}^{+}$ of the particles. The particle volume is limited through the minimal volume $v_{min}$ and the maximal volume $v_{max}$, which are physical properties determined by experiments. In theoretical works, $v_{min}=0$ and $v_{max}=\infty$ is widely used. The left-hand side of the equation is devoted to the space and time evolution of the density: it includes an advection term with the droplet velocity $v_{d}(t,x)$ and a diffusion term including a diffusion rate $D(t,x)$. On the right hand side $A^{+}$ and $A^{-}$ characterize the binary aggregation, where the birth term $A^{+}$ determines the creation of particles of volume $v$ and the death term $A^{-}$ describes disappearance of particles of volume $v$ in the population balance. $B^{+}$ and $B^{-}$ are the breakage terms [34]. The aggregation and breakage processes are modeled via the aggregation kernel $\omega(t,x,v,v^{\prime})$, the breakage frequency $\Gamma(t,x,v)$ and the daughter droplet distribution function $\beta(v,v^{\prime}).$ The breakage frequency $\Gamma(t,x,v)\geq 0$ represents the fractional number of breakage events per unit time of droplets of size $v$ while the daughter droplet distribution function $\beta(v,v^{\prime})\geq 0$ is the probability function for the creation of particles of size $v$ from particles of size $v^{\prime}$ [29].

We assume the following properties for the daughter droplet distribution function:

[TABLE]

The first property (2.2a) comes from the fact that no particle can split into larger particles. The function $N(v^{\prime})\in\mathbb{R}^{+}$ in (2.2b) represents the number of fragments obtained from the splitting of a particle of size $v^{\prime}$ [29]. Property (2.2c) takes into account that over time the total mass should be conserved by breakage events, while the total number of particles increases.

We furthermore set $\beta(v,v^{\prime})\equiv 0,~{}\Gamma(t,x,v)\equiv 0$ for all $v,v^{\prime}\notin[v_{min},v_{max}]$.

Binary aggregation appears if two particles of size $v_{1}$ and $v_{2}$ collide and merge to a particle of size $v_{1}+v_{2}$. The aggregation kernel $\omega$ provides the probability that a collision of two particles results in a successful merge. We again assume the following conditions:

[TABLE]

The symmetry condition (2.3b) comes from the fact that the two merging events $(v_{1},v_{2})\to v_{1}+v_{2}$ and $(v_{2},v_{1})\to v_{1}+v_{2}$ are equivalent (symmetry of binary aggregation events).

Formulas for our choice of kernels can be found in section 5.

3 Moment Models and Realizability

The method of moments is commonly used to derive reduced models for kinetic transport equations [47, 58, 23, 25, 53]. To solve such a reduced model is often less expensive, regarding computational time, than to solve the complete model. Equation (2.1) depends on time, space and the volume. Here, the method of moments consists in a (nonlinear) projection of the solution $f$ with respect to a polynomial basis in the volume variable $v$.

We start with some definitions that closely follow [36, 5, 4].

We denote by $\gamma=\left(\gamma_{0},\gamma_{1},...,\gamma_{N}\right)^{T}:\mathbb{R}^{+}\times\mathbb{R}^{n_{d}}\to\mathbb{R}^{N+1}$ the vector containing entries of the first $N+1$ moments of $f$ with respect to some polynomial basis $m=m(v)=\left(m_{0}(v),...,m_{N}(v)\right)^{T}$ of $\mathbb{P}_{N}([v_{min},v_{max}])$111Other bases are also possible, see, e.g., [47]., that means

[TABLE]

where the integration is applied component-wise. Then the system of $N+1$ moment equations for the population balance equation is derived by multiplying (2.1) by the basis $m$ and integrating over the volume domain. This results in

[TABLE]

The moments $\gamma_{i}$ depend on the distribution function that appears on the right hand side of (3.2) in the aggregation and breakage terms. Here, an ansatz for the distribution function, which depends on the known moments, is needed to close the equations.

3.1 Realizability

When dealing with closures, the question arises which moments can actually be realized by a physical (non-negative) distribution function, compare [18, 3]. In other words, for given moments $\gamma_{i}$ one has to find a non-negative distribution function $f(v)\geq 0$ such that

[TABLE]

The realizability domain

[TABLE]

is defined as the set of vectors $(\gamma_{0},...,\gamma_{N})^{T}\in\mathbb{R}^{N+1}$ such that this problem has a solution. We refer to [18, 3] for further reading and explicit characterizations of $\textbf{{R}}_{N}$ in terms of the moments $\gamma$. Note that the realizable set changes when $\langle\cdot\rangle$ is replaced by a numerical quadrature, see section 4.3 and [4, 6]

3.2 Moment Closures

Assuming a specific form of the distribution function $f$ in the integrals appearing in (3.2), which depends on the available moments, one can close the equations. There exist many types of closures in literature, each with their own advantages and disadvantages. In this section we will explain the polynomial closure $P_{N}$, the maximum entropy closure $M_{N}$ as well as the quadrature method of moments $QMOM$.

3.2.1 Polynomial Closure

The $P_{N}$ equations can be most easily understood as a Galerkin semi-discretization in $v$. The idea is to expand the distribution function in terms of a truncated series expansion

[TABLE]

Inserting $f_{P_{N}}$ into the definition of moments (3.1) leads to

[TABLE]

The basis $m$ can be chosen to be any kind of polynomial basis, for example the monomials $v^{i}$. But to simplify the solution of the linear system of equations (3.3) for the expansion coefficients $a_{i}$, we use an orthogonal basis. Choosing the Legendre polynomials and shifting them to the interval $[v_{min},v_{max}]$ gives us the following recursion formula

[TABLE]

Therefore, the moments $\gamma_{i}$ for $i=0,...,N$ can be expressed as

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. Solving for the coefficients $a_{i}$, the approximated distribution function for the $P_{N}$ closure looks like

[TABLE]

The polynomial closure does not necessarily provide a non-negative distribution function, thereby the computed moments could get physically meaningless, like e.g. a negative particle number.

3.2.2 Maximum Entropy

Another kind of closure for the approximation of the distribution $f$ are entropy-based methods. Here, $f_{M_{N}}$ is the solution of a constrained, convex optimization problem [25]. Namely, for strictly convex $\eta:\mathbb{R}\to\mathbb{R}$ the optimization problem looks like

[TABLE]

Transforming (3.6) into its unconstrained, strictly convex dual problem we end at searching for the Lagrange multipliers $\hat{\alpha}(\gamma)\in\mathbb{R}^{N+1}$ that solve

[TABLE]

Using the Legendre dual $\eta_{*}:\mathbb{R}\to\mathbb{R}$ of $\eta$, this is equivalent to

[TABLE]

So, if a solution of problem (3.7) exists, it takes the form [36]

[TABLE]

Choosing the Maxwell-Boltzmann entropy $\eta(g)=g\log(g)-g$, the Legendre dual and its derivative are of the form $\eta_{*}(y)=\eta_{*}^{\prime}(y)=e^{y}$ and therefore the solution of (3.7) can be written as

[TABLE]

We define the dual objective function $h:\mathbb{R}^{N+1}\to\mathbb{R}$ (the function which gets minimized) by

[TABLE]

Through the properties of the exponential function it is ensured that $f_{M_{N}}$ is always a positive function. Because of the strict convexity, problem (3.6) does have a unique solution whenever $\gamma$ is realizable [25].

3.2.3 QMOM

In $QMOM$ the general idea is to use an n-atomic discrete measure as ansatz. This means that the approximated distribution function consists only of a linear combination of Dirac $\delta$-functions

[TABLE]

with $w_{j}>0$ and $v_{j}\in[v_{min},v_{max}]$ for all $j=1,...,n.$ Plugging $f_{QMOM_{N}}$ into (3.1) and choosing $m_{i}(v)=v^{i}$ leads to

[TABLE]

The above nonlinear system can be solved using the Wheeler-algorithm [56, 59, 18, 3] which diagonalizes a tridiagonal matrix to find the weights $w_{i}$ and abscissas $v_{i}$. This results in a robust and efficient algorithm for the inversion of the moment problem, which will only succeed if $\gamma$ is realizable. It can be shown that a moment vector on the realizability boundary can be uniquely represented by an atomic distribution function. Thus, $QMOM$ is able to exactly reproduce this behavior in such a case [53, 28].

4 Numerical Implementation

4.1 Finite Volume Scheme for the PBE

We will shortly recall the finite volume scheme (FVS) for the one-dimensional PBE with binary aggregation and multiple breakages on uniform meshes, as introduced in [21, 10]. Its solution will be used as a reference for our numerical comparison of the moment closure schemes.

Extensions of this scheme to different types of uniform and non-uniform grids and analytical results can be found in [32, 34, 30].

At first, equation (2.1) is rewritten in conservative form [34] for the mass density $vf(t,v)$,

[TABLE]

with the continuous aggregation flux

[TABLE]

and the continuous breakage flux

[TABLE]

The volume domain $[v_{min},v_{max}]$ is discretized into equidistant cells

[TABLE]

The grid size is denoted by $\Delta v=(v_{max}-v_{min})/n_{v}$ and the midpoints $v_{i}$ are computed by $v_{i}=\left(v_{i+\frac{1}{2}}+v_{i-\frac{1}{2}}\right)/2$. Let $f_{i}(t)$ be the semi-discrete approximation of the cell average of the solution on cell $i$ [32]

[TABLE]

The FVS is derived by integrating the conservation law over every cell $\Lambda_{i},$

[TABLE]

The numerical fluxes $J_{i+\frac{1}{2}}$ are chosen as appropriate approximations of the continuous flux functions [32, 33, 34, 29, 21, 21]. The numerical breakage flux takes the form [34]

[TABLE]

Moreover, $J_{1/2}^{brk}(t,x)=J_{n_{v}+1/2}^{brk}(t,x)=0$ since at the boundary we have $F^{brk}(t,x,v_{min})=F^{brk}(t,x,v_{max})=0.$

The numerical aggregation flux is given by [21, 34]

[TABLE]

with $\alpha_{i,l}=i-l-\zeta+2$. The index $\zeta$ is determined by $v_{i+\frac{1}{2}}-v_{l}=v_{i-l+1}-v_{min}\in\Lambda_{i-l+1-\zeta}$ and

[TABLE]

Again, $J_{1/2}^{agg}(t,x)=J^{agg}_{n_{v}+1/2}(t,x)=0$, since at the boundaries we have $F^{agg}(t,x,v_{min})=F^{agg}(t,x,v_{max})=0$. We remark that restrictions of the size of the time-step to ensure positivity of the distribution function have been investigated in [21, 29].

4.2 Numerical approximation of the moment system

We reconsider the moment equations (3.2). For a better numerical approximation we transform the second term on the right hand side of (3.2) towards

[TABLE]

with $g_{A^{-}}(v,x)=(v_{max}-v)(x-v_{min})/(v_{max}-v_{min})+v_{min}$. In general, one is not able to compute the above integrals exactly. They have to be evaluated numerically. As quadrature rule for the computation of the terms on the right hand side of (3.2) we use a Gauß Lobatto formula, which integrates polynomials up to degree $2n_{Q}-3$ exactly for $n_{Q}$ quadrature points.

For the $QMOM_{N}$, whose closure is based on Dirac $\delta$ functions, we approximate the integrals by

[TABLE]

Using property (2.2c) of the $\beta-$kernel, we can easily see that mass conservation is exactly fulfilled in case of the quadrature moment method, when $m$ is the monomial basis:

[TABLE]

4.3 Numerical Realizability

Since for maximum-entropy models almost all of the integrals $\langle\cdot\rangle$ have to be evaluated numerically, we need to adapt the definition of the realizable set.

For a function $g:[v_{min},v_{max}]\to\mathbb{R}$ the nodes $\left\{v_{i}\right\}_{i=1}^{n_{Q}}$ and weights $\left\{w_{i}\right\}_{i=1}^{n_{Q}}$ of a quadrature rule Q are chosen such that $\langle g\rangle$ is approximated by

[TABLE]

By abuse of notation, $\langle\cdot\rangle$ should be understood as $\textbf{{Q}}(\cdot)$ whenever needed.

For an arbitrary quadrature rule Q the Q-realizable set is defined as [4, 6]

[TABLE]

It is worth mentioning that $\textbf{{R}}_{N}^{\textbf{{Q}}}$ is a strict (polytopic) subset of $\textbf{{R}}_{N}$ and like $\textbf{{R}}_{N}$, it is an open convex cone [4]. $\textbf{{R}}_{N}^{\textbf{{Q}}}$ depends on the choice of the quadrature nodes and, in particular, on the number of points $n_{Q}$.

4.4 Realizability-preserving property of the schemes

As we have noted above, the maximum-entropy moment problem (3.1) has a solution if and only if the moments $\gamma$ are realizable (in $\textbf{{R}}_{N}$ or $\textbf{{R}}_{N}^{\textbf{{Q}}}$, respectively). Since we need to be able to solve the moment problem to evaluate the right hand side of (3.2), it is crucial to maintain realizability throughout the computation. Similar problems have been treated in the context of radiative transfer, e.g., in [6, 46, 48, 16, 41].

Theorem 4.1**.**

Under the non-negativity of the kernels $\Gamma$ and $\omega$ for all $v\in[v_{min},v_{max}]$ and neglected space dependency, the explicit Euler scheme discretization of the moment equations $\eqref{ME}$, combined with $M_{N}$ or $QMOM_{n}$, preserves realizability under the CFL-like condition

[TABLE]

that is, given initial data in $\textbf{{R}}_{N}^{\textbf{{Q}}}$ ($\textbf{{R}}_{N}$, respectively), one step of the forward Euler scheme provides updated values that are also in $\textbf{{R}}_{N}^{\textbf{{Q}}}$ ($\textbf{{R}}_{N}$, respectively).

The proof of this theorem can be obtained similarly to the work [21, 30] and [5, 6] and is therefore omitted.

4.5 Implementation details for maximum entropy models

For the numerical implementation of the maximum entropy model two big questions arise: how to solve the optimization problem and what to do if the algorithm fails to converge (in a reasonable time). The optimization problem is commonly solved by variations of the Newton algorithm [4, 5, 1, 52]. Problems occur for moment vectors which are near the realizability boundary. In this case, the algorithm may require a large number of Newton iterations to converge or does not converge at all. Also, the solution will be sensitive to small changes in the moments, caused by an ill-conditioned Hessian matrix of the dual objective function. The problem is impaired by the use of a quadrature rule and the finite-precision arithmetic to evaluate the integrals [4]. There are several variations of Newton’s method available that deal with the singularity of the Hessian matrix [50, 2, 5]. We follow the ideas in [4]. There, the authors use an adaptive change of basis together with a regularization method.

4.5.1 Newton’s Method with an Adaptive Change of Basis

We slightly modify the ideas and algorithms proposed in [4], where more details and results can be found. The optimization algorithm we use is based on Newton’s method stabilized by an Armijo backtracking line search [7] and an adaptive change of basis to improve the condition number of the Hessian. It computes an approximation $\alpha$ of the true solution $\hat{\alpha}.$ If the algorithm fails to converge we use some regularization technique.

Let us recall the dual objective function

[TABLE]

Its gradient $g:\mathbb{R}^{N+1}\to\mathbb{R}^{N+1},$ Hessian $H:\mathbb{R}^{N+1}\to\mathbb{R}^{(N+1)\times(N+1)}$ and its Newton direction $d:\mathbb{R}^{N+1}\to\mathbb{R}^{N+1}$ are given by

[TABLE]

4.5.2 Initial Guess of the Optimization Algorithm and Regularization

First, we define two bases $m_{mono}=(1,v,v^{2},...,v^{N})^{T}$ and $m_{o}$, whereas $m_{mono}$ is used to evaluate the integral terms of (3.2) and to compute the moments. The basis $m_{o}$ is needed in the optimization algorithm as the initial orthogonal triangular basis to start the algorithm. Then the iterative basis remains orthogonal and triangular. We choose as basis the Legendre basis (3.4). For the Hessian matrix to have full rank it is necessary that $n_{Q}\geq N+1$ [4]. In [4] the authors proposed to initially start with the Legendre basis and then save the new adapted basis computed in the Newton algorithm to use it in the next time step as initial basis. Similarly, they handled the approximated Lagrange multipliers $\alpha$. We proceed in a slightly different way. We use in every Newton algorithm the Legendre basis $m_{o}$ as initial basis no matter what time step we have. We then compute an initial guess $\alpha_{M_{1}}$ for the Lagrange multipliers involving only the first two moments $\gamma_{0}$, $\gamma_{1}$ with the help of the Newton algorithm. Here we use as initial guess the Legendre basis and $\alpha_{0,M_{1}}=\left(\ln(\gamma_{0}),0\right)^{T}.$ Our initial guess for the total set of moments is $\alpha_{0}=\left(\alpha_{M_{1}},0,...,0\right)^{T}\in\mathbb{R}^{N+1}.$ Of course for the actual optimization algorithm we need those multipliers with respect to the orthogonal basis $\beta_{0}=(m_{o}^{-T}\alpha_{0}).$ If the moments are too close to the realizability boundary, the algorithm could fail, therefore in this case a regularization is used. We define a realizable moment $Q_{mono}$ such that we can derive a new realizable moment which is farther away from the realizability boundary through a convex combination of the original moment $\gamma$ and $Q_{mono}.$ Note that since $\textbf{{R}}_{N}$, $\textbf{{R}}_{N}^{\textbf{{Q}}}$ both are convex cones, convex combination of realizable moments are again realizable. Therefore, we again take advantage of $\alpha_{0}$ and compute the corresponding moments both in the monomial basis and in the Legendre basis

[TABLE]

So starting the Newton algorithm with some moments in the monomial basis $\gamma$ and their corresponding Legendre moments $\tilde{\gamma},$ the regularization technique we use for $r\in[0,1]$ has the following form

[TABLE]

It is implemented in this way because for the reduced moment equations (3.2), mass conservation should be guaranteed, which means $\gamma_{1}$ remains constant. The Lagrange multiplier $\alpha_{0}$ exactly reproduces the zeroth, and first moment. So despite the regularization technique they do not change.

4.5.3 Stopping Criterion

For the stopping criterion, we compute the gradient in the $k-$th iteration by

[TABLE]

We stop the algorithm if

[TABLE]

The last two terms arise from the regularization technique. If the algorithm stops and regularization was needed we replace our original moments $\gamma$ by $(1-r)\gamma+rQ_{mono}.$ The rest of the algorithm follows the one shown in [4].

4.6 Numerical Treatment of the population balance equation and its moment equations

We consider the two-dimensional case. Let $\Omega=[x_{min},x_{max}]\times[y_{min},y_{max}]$ be our rectangular spatial domain. We discretize it by a Cartesian equidistant grid with cells $C_{i,j}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}],$ where $x_{i+1\frac{1}{2}}=x_{min}+i\Delta x,$ $i=0,...,N_{x}$ and $y_{j+\frac{1}{2}}=y_{min}+j\Delta y,~{}j=0,...,N_{y}$ with the step sizes $\Delta x=(x_{max}-x_{min})/N_{x}$ and $\Delta y$ correspondingly. The cell centers of such a cell $C_{i,j}$ are computed by $(x_{i},y_{j})=(\frac{1}{2}(x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}}),\frac{1}{2}(y_{j+\frac{1}{2}}+y_{j+\frac{1}{2}}))$. The numerical approximation of the average of $f$ over a cell $C_{i,j}$ is defined by

[TABLE]

$T$ is the final time until which the numerical evaluation of the solution is calculated. The time domain $[0,T]$ is discretized by the points of time $t^{k}=k\Delta t$ with the time step size $\Delta t.$

4.6.1 Positivity-Preserving Discretization of the Population Balance Equation

In this section we derive a positivity-preserving numerical scheme for the spatially inhomogeneous population balance equations neglecting here for simplicity the aggregation and breakage terms, which have been treated in a previous section.

[TABLE]

The droplet velocity is defined by $v_{d}(t,x)=(u(t,x),z(t,x))^{T}$ and the diffusion rate by $D=D(t,x)>0.$ Integrating equation (4.2) over the cells $[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]\times[t_{k},t_{k+1}]$ and the use of the mid point rule gives the finite volume formulation

[TABLE]

We exemplary show for $F_{i+\frac{1}{2},j}^{k}$ and $H_{i+\frac{1}{2},j}^{k}$ how to handle these terms. $G$ and $M$ work similarly. Let us first address the advection flux in $x$-direction $F_{i+\frac{1}{2},j}^{k}$. We split droplets into left- and right-moving particles and get for the numerical flux

[TABLE]

where $u_{i+\frac{1}{2},j}^{k}=1/2\left(u_{i,j}^{k}+u_{i+1,j}^{k}\right).$ The values $f^{+,k}_{i+\frac{1}{2},j}$, $f^{-,k}_{i+\frac{1}{2},j}$ on the right and left sides of the cell edge at $(x_{i+\frac{1}{2}},y_{j})$ [25, 23], are approximated by

[TABLE]

with $\sigma_{i,j}^{x,k}$ the approximation of the slope in the $x-$direction in cell $C_{i,j}$ at time $t_{k},$

[TABLE]

The minmod limiter does not only guarantee non-negativity of $f^{+,k}_{i+\frac{1}{2},j}$, $f^{-,k}_{i+\frac{1}{2},j}$ but also suppresses spurious oscillations [49]. The diffusion term is treated in a similar way as in [15],

[TABLE]

4.6.2 Realizability Preserving Discretization of the Moment Equations

For the maximum entropy closure or the quadrature method of moments we did discuss the realizability preservation of the spatially homogeneous case in subsection 4.3. Here, a realizability preserving scheme for the advection-diffusion equation neglecting the breakage and aggregation terms is given. Consider

[TABLE]

Assuming that $\gamma_{i,j}^{k}$ is realizable and that $\tilde{f}$ is the corresponding distribution function approximated by the maximum entropy model. Then to derive a realizable scheme we multiply equation (4.3) with the basis $m$ and integrate over the whole volume domain [25, 48], which gives

[TABLE]

the terms $\tilde{F}_{i\pm\frac{1}{2},j}^{k},~{}\tilde{G}_{i,j\pm\frac{1}{2}}^{k},~{}\tilde{H}_{i\pm\frac{1}{2},j}^{k},~{}\tilde{M}_{i,j\pm\frac{1}{2}}^{k}$ are in the same way defined as $F_{i\pm\frac{1}{2},j}^{k},~{}G_{i,j\pm\frac{1}{2}}^{k},~{}H_{i\pm\frac{1}{2},j}^{k},~{}M_{i,j\pm\frac{1}{2}}^{k}$ but with respect to the approximate distribution function $\tilde{f}$. For the polynomial closure $P_{N}$ we can also use this scheme, but there no realizability can be guaranteed, since the distribution function can get negative.

A similar scheme as in (4.5) applied to the moment equations (4.4) with the $QMOM$ closure without aggregation and breakage terms can be shown to be realizability preserving in every time step [54, 55]. Recall that the $QMOM-$distribution function is of the form $\tilde{f}(t,x,v)=\sum_{\alpha}w_{\alpha}(t,x)\delta\left(v-v_{\alpha}(t,x)\right).$ The discretization is now similar to the discretization method (4.5). One only has to change the definition of the terms $\tilde{f}_{i\pm\frac{1}{2},j}^{+/-,k},~{}\tilde{f}_{i,j\pm\frac{1}{2}}^{+/-,k}$. They are determined by

[TABLE]

The interface weights are defined as

[TABLE]

Again we use the minmod limiter to reconstruct the slopes as

[TABLE]

The discretization in the $y-$direction works analogously.

Theorem 4.2**.**

To ensure either a non-negative distribution function for the advection-diffusion equation (4.2) discretized by the finite volume scheme (4.3) or a realizability preservation of the schemes for the moment equations with $QMOM$, or the maximum entropy closures of the advection diffusion problem (4.4), the following time step restriction has to be fulfilled

[TABLE]

with

[TABLE]

and

[TABLE]

5 Numerical Results

We compare the zeroth moment of the different closures $QMOM$ (3.9), maximum entropy (3.8) and $P_{N}$ (3.5) with each other. As reference solution we use the finite volume scheme explained in section 4.1. The zeroth moment for the finite volume scheme is computed by

[TABLE]

We will compare them regarding the $L_{2}$-error as well as computation time. The relative error in the $L^{2}$ sense for two number densities $\gamma_{0,\text{ref}}$, $\gamma_{0,\text{clo}}$ is evaluated by the formula

[TABLE]

where $\Omega$ is the spatial domain and $T$ the final time until which we compute the numerical solution. For the Newton-algorithm in the maximum entropy case we need to choose several parameters which can be found in table 1. For the explanation of the parameters see [4].

For all examples we use the a daughter droplet distribution which is based on the purely statistical daughter droplet distribution function of Hill and Ng [26]. Recall that $N(v^{\prime})$ is the average number a particle of size $v^{\prime}$ splits. Our daughter droplet distribution is a weighted sum of daughter droplet functions

[TABLE]

where the $i-$breakup process (a mother droplet splits in average into $i\in\mathbb{N}$ smaller droplets) daughter distribution functions for $v_{max}\geq v^{\prime}>pv_{min}$ are defined as

[TABLE]

[TABLE]

where $m\in\mathbb{N}$ determines the shape of $\tilde{\beta}_{i}.$ The no-breakup daughter distribution function is defined for $v_{max}\geq v^{\prime}\geq v_{min}$ as

[TABLE]

For all $v^{\prime}\in[v_{min},v_{max}]$ the weight functions $g_{i}(v^{\prime})\in[0,1]$ have to fulfill

[TABLE]

and $c(v^{\prime})\in[2,\infty)$. We choose $c(v^{\prime})$ as $c(v^{\prime})=2N(v^{\prime})-1.$ For the computation of the weight functions let $\textbf{1}_{U}(x)$ be the indicator function on $U$. We choose $g_{N(v^{\prime})-j}=g_{N(v^{\prime})+j}$ for all $j=1,...,N(v^{\prime})-1$ and for all $v^{\prime}\in[v_{min},v_{max}]$ and assume that the weight functions are of the form

[TABLE]

with $k(v^{\prime},N(v^{\prime}))\in\mathbb{N}$ determining that the first $k(v^{\prime},N(v^{\prime}))-1$ and last $k(v^{\prime},N(v^{\prime}))-1$ weight functions of the series $g_{1}(v^{\prime}),..,g_{2N(v^{\prime})-1}(v^{\prime})$ are zero for $v^{\prime}$. With the first property (5.2) of the weight functions, the variable $j$ is determined by

[TABLE]

Let us define the intervals

[TABLE]

Then, if $v^{\prime}\in I_{l}$ for some $l\in\mathbb{N}$ we define

[TABLE]

This definition fulfills all the properties mentioned in section 2 .

Although the choice of the weight functions seems complicated there are just chosen such that they are symmetric about $N(v^{\prime})$ and the largest weight corresponds to $N(v^{\prime})$ for all $v^{\prime}$ and such that all properties mentioned in section 2 are fulfilled. The choice also guarantees that droplets cannot break into droplets whose volume is not contained in $\left[v_{min},v_{max}\right]$. Figure 1 demonstrates the behavior of the weight functions $g_{i}(v^{\prime})$ on the example of quasi-ternary breakage.

5.1 Pure Breakage

First we only want to consider a pure breakage problem without spatial dependency. This means that we set $\omega\equiv 0$ in (3.2) and (2.1). Our initial condition is a normal distribution

[TABLE]

For the daughter droplet distribution function we choose binary breakage $p=2$ and $m=2$. Table 2 contains some of the parameters chosen for this example.

The time step size is for all methods $\Delta t=0.01$. We compute the solutions until $T=4$ is reached. The number of points of the quadrature rule is $n_{Q}=100$ and for the reference solution computed by the finite volume scheme we use $n_{v}=5000$ grid points. As breakage frequency we choose the one from Coulaloglou and Tavlarides [17]

[TABLE]

Table 3 shows the values corresponding to the parameters for the breakage frequency.

In Figure 2 the results are shown for the comparison of $P_{N}$ (magenta line), $M_{N}$ (green line) and $QMOM_{N}$ (black line). We computed the density for all three models with different values for the order, namely up to order thirteenth. The figure is divided into three subfigures: (a) showing the order $N$ plotted against the computation time, (b) illustrating the relation between the order and the relative $L_{2}$ error (5.1) and (c) showing the computation time plotted against the relative $L_{2}$ error (5.1). It is demonstrated that with increasing order the computation time also increases. For this situation $M_{N}$ and $QMOM_{N}$ have much higher computation times than $P_{N}$. We note that the optimization algorithm for $M_{N}$ is the dominating part of the computation. In case of $QMOM_{N}$ the solution of the nonlinear system is the costly part. This results in the higher computation times for $M_{N}$ and $QMOM_{N}$. On the other hand, as mentioned above, the $P_{N}$ approach might lead to unphysical moments.

5.2 Pure Aggregation

For the pure aggregation problem we choose $\Gamma\equiv 0,~{}\beta\equiv 0$ in (3.2) and (2.1) and neglect spatial dependency. We use an aggregation kernel proposed by Coulaloglou and Tavlarides [17]

[TABLE]

The parameters remain the same as in the pure breakage case see table 2. Table 4 shows the parameters for the aggregation kernel.

As before, we choose $n_{Q}=100$, $\Delta t=0.01$, $n_{v}=5000$. The comparison of $P_{N}$ (magenta line), $M_{N}$ (green line) and $QMOM_{N}$ (black line) is illustrated in figure 3. Again the densities for all three models are computed with changing values for the order, namely up to the ninth order. The figure is subdivided in the same way as figure 2. It is demonstrated that with increasing order the computation time also increases. In this case again the computation times of $M_{N}$ and $QMOM_{N}$ are much larger than those of $P_{N}$. As before the optimization algorithm is the dominating part of $M_{N}$, resulting in the higher computation time. If we compare figures 2 and 3, one observes that the relative error in the $QMOM_{N}$ case is much higher. This is due to the fact, that the breakage integrals are less complex then the aggregation integrals, where the distribution function has to be evaluated twice.

5.3 Coupled Diffusion-Advection-Breakage-Aggregation

For this example we consider equations (2.1) and (3.2), we choose the diffusion rate to be constant $D(t,x)=D=0.001.$ We assume that the droplet velocity $v_{d}(t,x)$ is the same as the velocity of the fluid computed by the lid driven cavity problem [11] with a Reynold’s number of $Re=5.$ The spatial domain is a square $\Omega=[0,1]\times[0,1]$. We discretize it with $50$ gridpoints in both directions. We choose as breakage frequency (5.4) with the same parameters as in table 3. For the daughter droplet distribution function we choose again $p=2$ and $m=2$. The aggregation kernel is the one from (5.5) of the pure aggregation problem with the parameters from table 4. Again, $n_{Q}=100.$ The reference solution is computed via the finite volume scheme from section 4.1 with $n_{v}=2000$.

We compare the densities of the finite volume schemes with those of the moment methods, namely $QMOM$ (3.9), maximum entropy (3.8) and $P_{N}$ (3.5). The numerical simulation for all methods runs until the final time $T=5$ is reached. As initial distribution function we choose a Gaussian distribution,

[TABLE]

the parameters for $f$ are chosen as in table 2.

In figure 4 and 5 the time evolution of the densities for the different methods, namely FVS (a), $M_{1}$ (b), $P_{1}$ (c), $QMOM_{2}$ (d) are exemplary shown for two different points of time.

With increasing time the densities move from the lower left corner to the upper left corner, following the surrounding fluid. We can see that $P_{1}$ and $M_{1}$ look similar to the reference solution computed by the FVS. For $QMOM_{2}$ the solution shows bigger differences compared to the reference solution.

The comparison of $P_{N}$ (magenta line), $M_{N}$ (green line) and $QMOM_{N}$ (black line) is illustrated in figure 6. Again, the densities for all three models are computed with changing values for the order. The numerical solution is computed up to ninth order. The figure is subdivided in the same way as figure 2. It is demonstrated that with increasing order the computation time also increases and the relative $L_{2}$ error decreases. Again the Newton algorithm is responsible for the big difference in computation time between $M_{N}$ and $P_{N}.$

Remark 5.1**.**

The $P_{N}$ method is equivalent to the well-known discrete ordinates in case of radiative transfer equations [14], implying that meaningful initial conditions (i.e. those where the underlying $P_{N}$ distribution is positive for all $v$) will stay meaningful. While this is the case in the examples presented, this cannot be guaranteed in all cases. This requires to carefully choose the situations where the $P_{N}$ closure can be applied, if positivity of the kinetic solution is required. If the initial conditions are such that the underlying $P_{N}$ distribution is not positive, one has to choose $M_{N}$ or the $QMOM$ model which do not suffer from the above mentioned drawback [27].

6 Conclusions

We have compared a linear moment closure, a nonlinear maximum entropy moment closure and the quadrature method of moments numerically in terms of accuracy, order of the moments and computation times for a population balance equation. In the cases where breakage is involved the computation time of the quadrature method of moments and the maximum entropy closure are alike, in contrast to the pure aggregation example. In all cases, the linear closure is the fastest of the three methods. On the other hand, for the linear moment closure the underlying distribution function is not necessarily positive, which can lead to meaningless solutions in some applications. The error produced by the linear closure and the maximum entropy closure are smaller than the ones produced by the quadrature method of moments.

Future research should include other versions of $QMOM$, like $EQMOM$ or $SQMOM$ [8, 59], investigating if the classical moment models maintain their superiority. Furthermore, low-order ($N\leq 2$) partial moment methods on a partition of the velocity space should be taken into account, simplifying the inversion of the moment problem. These models should be compared with similar low-order $QMOM$ models regarding error and run-time. Additionally, multi-dimensional PBEs and their moment approximations should be investigated. Finally, the comparison of our simulations with real experiments would be interesting.

Acknowledgements

Funding by the Deutsche Forschungsgemeinschaft (DFG) within the RTG GrK 1932 ”Stochastic Models for Innovations in the Engineering Sciences” and by BMBF Verbundprojekt 05M2016-proMT is gratefully acknowledged.