Fast integration method for averaging polydisperse bubble population dynamics

TL;DR

This work addresses the computational bottleneck of evaluating ensemble-averaged moments in polydisperse bubbly flows, where highly oscillatory integrands arise from a distribution of bubble sizes. It introduces Levin's collocation method, which separates the oscillatory bubble-dynamics component from a smooth amplitude by solving a non-oscillatory differential system using spectral (Chebyshev) collocation; the resulting moments are obtained from boundary evaluations, avoiding the need to solve advection PDEs at many quadrature points. The method demonstrates substantial efficiency gains, achieving comparable or superior accuracy with orders of magnitude fewer evaluation points (e.g., 10^2–10^4× fewer nodes) and maintaining constant cost as oscillations grow over time, thereby enabling long-time, large-scale bubbly-flow simulations. While the approach relies on linearized dynamics and neglects several complex phenomena (translation, coalescence, breakup), it offers a robust route to fast moment computation and can be extended to include more complex dynamics and time-dependent distributions, with broad implications for Euler–Euler bubbly flow solvers.

Abstract

Ensemble-averaged polydisperse bubbly flow models require statistical moments of the evolving bubble size distribution. Under step forcing, these moments reach statistical equilibrium in finite time. However, the transitional phase before equilibrium and cases with time-dependent forcing are required to predict flow in engineering applications. Computing these moments is expensive because the integrands are highly oscillatory, even when the bubble dynamics are linear. Ensemble-averaged models compute these moments at each grid point and time step, making cost reduction important for large-scale bubbly flow simulations. Traditional methods evaluate the integrals via traditional quadrature rules. This approach requires a large number of quadrature nodes in the equilibrium bubble size, each equipped with its own advection partial differential equation (PDE), resulting in significant computational expense. We formulate a Levin collocation method to reduce this cost. Given the differential equation associated with the integrand, or moment, the method approximates it by evaluating its derivative via polynomial collocation. The differential matrix and amplitude function are well-suited to numerical differentiation via collocation, and so the computation is comparatively cheap. For an example excited polydisperse bubble population, the first moment is computed with the presented method at $10^{-3}$ relative error with 100 times fewer quadrature nodes than the trapezoidal rule. The gap increases for smaller target relative errors: the Levin method requires $10^4$ times fewer points for a relative error of $10^{-8}$. The formulated method maintains constant cost as the integrands become more oscillatory with time, making it particularly attractive for long-time simulations.

Figures (6)

• Figure 1: Schematic of ensemble-averaging approach for polydisperse bubble populations. Individual bubble dynamics are replaced by statistical moments of the bubble size distribution, which couple to the continuous liquid phase through mixture-averaged properties.
• Figure 2: Example integrands $f(R_o) R^i(t,R_o)$ at $t=200$ for different moment orders $i$, showing the oscillatory behavior that emerges from polydisperse bubble populations. The complex oscillatory structure makes traditional quadrature methods computationally expensive.
• Figure 3: Integrands for the Bessel oscillator verification case showing $J_k(rx)/(x^2 + 1)$ for different oscillation parameters $r$. As $r$ increases, the integrand becomes more oscillatory, demonstrating the challenge for traditional quadrature methods.
• Figure 4: Convergence comparison between the Levin collocation method and trapezoidal rule for the Bessel oscillator integral. The Levin method demonstrates exponential convergence and requires orders of magnitude fewer function evaluations to achieve comparable accuracy, especially for highly oscillatory integrands (larger $r$).
• Figure 5: Structure of differential matrices $\bm{A}$ for different moments $i$. As required for the efficient use of the Levin method, the matrices exhibit smooth, non-oscillatory behavior, with the oscillatory parts isolated into the bubble dynamics terms. Results correspond to the parameter values mentioned above.