QBMMlib: A library of quadrature-based moment methods

QBMMlib is an open source Mathematica package of quadrature-based moment methods and their algorithms. Such methods are commonly used to solve fully-coupled disperse flow and combustion problems, though formulating and closing the corresponding governing equations can be complex. QBMMlib aims to make analyzing these techniques simple and more accessible. Its routines use symbolic manipulation to formulate the moment transport equations for a population balance equation and a prescribed dynamical system. However, the resulting moment transport equations are unclosed. QBMMlib trades the moments for a set of quadrature points and weights via an inversion algorithm, of which several are available. Quadratures then closes the moment transport equations. Embedded code snippets show how to use QBMMlib, with the algorithm initialization and solution spanning just 13 total lines of code. Examples are shown and analyzed for linear harmonic oscillator and bubble dynamics problems.

PBEs can be solved by the method of classes [22,23] or the method of moments (MOM) [24,25]. QBMMlib employs the MOM because it can more naturally handle problems with multiple internal coordinates (e.g. velocities). Figure 1 shows a typical QBMM-based solution procedure. The MOM represents the NDF via a set of statistical moments and the transport equations for them follow from Email address: spencer@caltech.edu (Spencer H. Bryngelson) the PBE. Inverting the moments to a set of weights and abscissas provides a basis for approximating the unclosed transport equations via quadrature (QMOM) [26].
There is one other actively developed open source QBMM solver: OpenQBMM [38,39]. It is a library for Open-FOAM [40] and implements CQMOM and (3-node) CHyQ-MOM. MFiX [32,41] and Fluidity [42] use DQMOM, though modern conditional methods (e.g. CQMOM and CHyQMOM) generally outperform it [7]. Note that these are fully-coupled flow solvers. QBMMlib instead decouples these problems and solves the moment transport equations directly for an input dynamical system. This makes it preferable for prototyping and testing on novel physical problems. In pursuit of this, QBMMlib places emphasis  on expressive programming, simple interfaces, and symbolic computation where possible. As a result, it can solve PBE-based problems with just a few lines of code. Section 2 describes QBMMlib's implementation of the PBE and QBMMs. Section 3 verifies its methods and demonstrates its capabilities for three example problems. Section 4 discusses the utility and novelty of QBMMlib, which concludes the paper.

Software description
QBMMlib is a collection of Mathematica functions for solving PBEs via QBMMs. Table 2  MomentInvert inverts these moments to weights and abscissas that close the moment set via quadrature (ComputeRHS) and project it onto a realizable moment space (Project). The next sections describe the details of these routines.

Population balance equations (TransportTerms)
QBMMlib can solve one-and two-dimensional populations balance equations. A third direction can be added if its NDF is stationary. The two-dimensional case is detailed here without loss of generality. For illustration, consider where the dots indicate partial time derivatives, g is a function, and ξ = {x,ẋ} are the internal coordinates. The number density function f describes the state and statistics of this system in the ξ-space. A population balance equation (PBE) governs f as where the zero right-hand-side indicates conservation of f , though sinks and sources can model aggregation and breakup [35]. Quadrature-based methods to solving ( where p j (for j = 1, . . . , N ξ ) are the moment indices, N ξ is the number of internal coordinates (N ξ = 2 in (1)), and Ω MomentInvert Input: Moment set (moments, M ) and its indices Output: Optimal set of abscissas (xi, ξ) and weights (w, w) Options: Method, Permutation (default: 12 (ξ 1 |ξ 2 ), N ξ > 1 only)

OutAbscissa
Input: Abscissas Output: Threaded abscissas Table 2: Example public-facing routines. Parenthetical variables correspond to the code snippets and notation of section 2.
is the domain of f . These moments evolve as where, for (1), This forcing follows from the PBE via integration-by-parts [43].
For the prescribed dynamicsẍ (as in (1)), the integral term of (5) is equivalent to a sum of moments. For example, if The routine TransportTerms manipulates the PBE (as in 6) to determine the coefficients and moment indices that constitute F . The code snippet below demonstrates this functionality for the dynamics of 1. Here, the unassigned coefficients c [1] and c [2] correspond to the moment indices l and m of 5.

Moment inversion and quadrature weights (MomentIndex, MomentInvert)
MomentInvert inverts the set of raw moments M into a set of quadrature weights w and abscissas (nodes) ξ: Many algorithms can perform this procedure, each with its own relative merits, as discussed in section 1. Common approaches for one-dimensional moment sets (N ξ = 1) are QMOM [24,26] and hyperbolic QMOM (HyQMOM) [31]. For higher-dimensional moment sets (N ξ > 1) conditioned moment methods are cheaper and more stable than performing QMOM on each coordinate direction individually [36]. Such conditioned methods perform 1D moment inversion in one coordinate direction, then condition the next directions on the previous ones [36]. Examples of these are conditional-QMOM (CQMOM) and conditional-HyQMOM (CHyQMOM) [31,37]. The order that this conditioning is done is called the permutation. For 2D problems their are two permutations, ξ 1 |ξ 2 (coordinate direction ξ 2 conditioned on ξ 1 ) and the reverse, ξ 2 |ξ 1 .
The indices that makeup a so-called optimal set M depend on the moment inversion method and the number of quadrature nodes used in each internal coordinate direction N ξ . Here, "optimal" constrains the number of moments (and their order) required to yield a full-rank and square coefficient matrix [44]. Optimal moment sets are more stable and smaller (and so cheaper) than nonoptimal ones. For  Table 3 shows the optimal moment sets QBMMlib uses for N ξ = 2 two-dimensional problems [31,36]. MomentIndex computes these moment indices given the inversion algorithm (method) and number of quadrature points (n) in each internal coordinate direction (N ξ ). The Table 3: Example optimal N ξ = 2 moment sets for permutations as labeled. MomentInvert then inverts the moment set M to a set of weights and abscissas (as in (7)). In the following Mathematica code snippet, the moment set M (moments) is initialized via a two-dimensional Gaussian distribution (BinormalDistribution), though in principle M can be any realizable moment set. The method inversion algorithm then converts it to quadrature weights (w) and abscissas (xi).  Figure 2 shows the abscissas for different moment inversion algorithms. Their locations in the internal coordinate space ξ are different, though their weights are such that each quadrature reproduces the exact (up to) second-order moments. We verify that QBMMlib has this property to the carried precision. We discuss the QBMMlib quadrature routines used for this next.

Moment system closure via quadrature (ComputeRHS)
Quadrature approximates the raw moments defined in (3) and required by (5) as where ξ j,i (for i = 1, . . . , N ξj ) are the abscissa for internal coordinate direction ξ j (for j = 1, . . . , N ξ ). These quadrature approximations build F of (5) (F) via the QBMM function ComputeRHS. ComputeRHS approximates (via quadrature) and sums the required moments (exps) and their coefficients (coefs) for each moment index (momidx). The code snippet below shows this implementation in QBMMlib.

Realizable time integration (Project)
Stable and realizable time integration of (4) requires recasting the moment set M from its weights and abscissas [13]. QBMMlib function Project performs this projection. A time integrator (e.g. Euler's method) then computes the next iteration of the moment set (moments). The code snippet below shows an example QBMMlib projection and Euler time step (with time step size dt). QBMMlib also includes an adaptive strong-stability-preserving (SSP) third-order-accurate Runge-Kutta (RK) time integrator, RK23 [45]. The difference between the SSP-RK3 solution an embedded second-order-accurate SSP-RK solution provides a first-order approximation of the time step error. The time step size adjustment is then proportional to this error. The illustrative examples of the next section use this adaptive time stepping procedure.

Linear harmonic oscillator
The example case of section 2, including 6, is a linear harmonic oscillator. This moment system is linear and thus closed, so it can also verify the solution methods of QBMMlib via comparison to Monte Carlo simulations. Figure 3 shows the evolution of CQMOM and CHyQ-MOM moment sets and compares them to Monte Carlo surrogate truth solutions. The behavior of the first-order moments (M 10 and M 01 ) match the positions and velocities expected of a linear oscillator. Further, the QBMMlib solutions match the moments of the Monte Carlo simulations to plotting accuracy. The L 2 norm · 2 of the error quantifies this matching as where superscripts (MC) and (QMOM) are correspond to Monte Carlo and QBMMlib simulations, respectively. Figure 4 shows ε MC 2 for varying Monte Carlo ensemble size N MC and QBMMlib method CHyQMOM, though the Monte Carlo moment errors dominate the QBMM ones and so CQMOM has the same results. Indeed, the error converges at the expected rate.

Bubble cavitation
The dynamics of a cavitating gas bubble dispersion serves as a two-internal-coordinate nonlinear example problem. The Rayleigh-Plesset equation models the bubble dynamics [46]: where R is the bubble radius, Re is the Reynolds number (dimensionless ratio of inertial to viscous effects), and C p is the dimensionless pressure ratio between the suspending fluid and bubbles. Thus, R andṘ are the two internal coordinates (ξ). For our purposes it suffices to ignore surface tension effects (following (10)) and use C p = 1/0.3 to represent a relatively large pressure ratio. This formulation is non-dimensionalized by the (monodisperse) equilibrium bubble radius and suspending fluid density and pressure. The initial NDF is a log-normal distribution in the R-coordinate (shape parameters µ R = 1, σ R = 0.2) and a normal distribution in theṘ coordinate (µṘ = 0, σṘ = 0.1). The NDF is initially uncorrelated. Figure 5 shows the moment dynamics for two bubble dispersions problems: (a) viscous Re = 10 and (b) inviscid Re → ∞. Here, Re = 10 is the Reynolds number that corresponds to 1 µm bubbles in water and Re → ∞ represents ignoring viscous effects. Invoking Re → ∞ is not appropriate for most cavitation problems of physical relevance, though it provides a useful reference. In both cases the mean bubble radius M 10 oscillates and damps. This damping is more significant in (a) than (b) due to viscous effects, as expected. In the Re = 10 case, this is sufficient for the QBMMlib-predicted moments to match the Monte Carlo results. However, for Re → ∞, the evolving moment set is unable to faithfully represent the bubble oscillations, particularly at long times. Indeed, a mismatch between the Monte Carlo and QBMMlib results is clear for M 02 and M 11 . These differences are qualitatively similar for both the CQMOM and CHyQMOM algorithms. This is because closing the moment system requires extrapolating out of the represented moment space, which is of similar fidelity for both algorithms.

Impact and conclusions
This paper introduced QBMMlib, a library for solving PBEs using quadrature-based moment methods. It is a Wolfram Language package, which is useful for automating the procedure of using QBMMs for simulating phenomena like bubble and particle dynamics. This includes constructing a moment set for a given QBMM, determining the righthand-side functions corresponding to a governing equation automatically, and inverting the moment set for quadrature points to close the system. These routines leverage Mathematica's symbolic algebra features and include modern QMOM and conditional-QMOM methods. Having these features available in a unified framework is helpful, particularly when it is unclear what QBMM will be appropriate (or stable) for the model dynamics. Our searches suggest that QBMMlib is the only library, open source or otherwise, that provides such capabilities. Given this, QBMMlib should help researchers prototyping QBMMs for their physical problems (or developing new QBMMs entirely). Indeed, the authors used QBMMlib to guide the implementation of CHyQMOM for phase-averaged bubble cavitation into MFC, the first flow solver with this capability [47,48].

Conflict of Interest
We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.