A robust high-resolution algorithm for quadrature-based moment methods applied to high-speed polydisperse multiphase flows

This paper presents a robust, high-resolution Eulerian algorithm for simulating high-speed polydisperse granular multiphase flows by coupling compressible gas dynamics with mass-based moment equations closed via the generalized quadrature method of moments, demonstrating its effectiveness through various complex shock-driven numerical experiments.

The Gist

Imagine you are trying to predict how a cloud of dust behaves when a shockwave hits it. Now, imagine that this dust isn't just made of identical grains of sand. Some grains are tiny like flour, some are medium like sugar, and some are huge like pebbles.

In the real world, when a shockwave hits this mix, the tiny grains zip away instantly, while the heavy pebbles lag behind. They separate, swirl, and interact in complex ways.

The Problem:
For a long time, computer models used to simulate these flows had to cheat. They assumed every particle was the exact same size (like a bag of identical marbles). This made the math easy, but it missed the most interesting part: size segregation. If you want to predict a volcanic eruption, a dust explosion, or a rocket engine's performance, you need to know how different-sized particles behave differently.

The Solution:
The authors of this paper built a new, super-smart computer algorithm. Think of it as a high-resolution camera for fluid dynamics. Instead of taking a blurry photo where everything looks the same, this camera can see the fine details of how tiny particles and big particles move separately, even when they are mixed together in a chaotic storm.

Here is how they did it, using some everyday analogies:

1. The "Magic List" (Quadrature Method of Moments)

Usually, to track every single particle in a cloud, you'd need a computer the size of a planet. That's impossible.
Instead, this algorithm uses a "Magic List." Imagine you have a bag of mixed nuts. Instead of counting every single almond and cashew, you just track a few key statistics: the average size, how spread out the sizes are, and how "lumpy" the distribution is.
The paper uses a mathematical trick called Quadrature Method of Moments (QMOM). It takes these statistics and turns them back into a "virtual list" of a few representative particle sizes (nodes). It's like saying, "We don't need to track 1 million particles; we just need to track 3 'super-particles' that represent the small, medium, and large groups perfectly."

2. The "Traffic Cop" (Riemann Solvers)

When a shockwave hits, things get chaotic. Information travels at different speeds.
The algorithm acts like a super-efficient Traffic Cop at a busy intersection. It solves a "Riemann Problem" (a fancy math term for "what happens when two different states crash into each other?") separately for the gas and for the particles.

• The Gas: It uses a high-speed solver (HLLC) to handle the air.
• The Particles: It treats each "super-particle" on its Magic List as its own little traffic system. It calculates how the small particles react to the shock, then how the medium ones react, and so on.
  By solving these separately and then stitching the answers together, the computer avoids getting confused or crashing (which is called "numerical instability").

3. The "Smart Zoom" (High-Order Reconstruction)

Older methods were like looking at a pixelated image. If a wall of dust moved, the computer would smear it out, making the edge look fuzzy.
This new method uses High-Order Reconstruction. Imagine drawing a line on a piece of paper. A low-order method draws a jagged, blocky line. A high-order method draws a perfectly smooth, sharp line.
The authors added a safety net: if the computer gets confused (like when a patch of dust is surrounded by empty air, creating an "island"), it automatically slows down to a simpler, safer method to prevent errors, then speeds back up when things get stable.

4. The "Real-World Tests"

The team didn't just write the code; they put it through the wringer with extreme scenarios:

• The Dust Curtain: A shockwave hits a wall of dust. The model correctly showed the small dust flying ahead while the big dust lagged behind.
• The Dust Explosion: They simulated a shockwave hitting a pile of dust. The model showed how the dust "rolls up" into vortices (swirls), just like real experiments.
• The Fireball: They simulated a high-pressure gas exploding inside a shell of particles. The model showed how the explosion pushed the particles out, creating "fingers" of dust and separating sizes as the fireball expanded and then collapsed.

Why Does This Matter?

This isn't just about dusty clouds. This technology helps us:

• Design safer rockets: Understanding how fuel particles burn and mix.
• Predict volcanic eruptions: Knowing how ash clouds travel and separate by size helps predict which areas will be covered in heavy rocks vs. fine ash.
• Prevent industrial accidents: Simulating dust explosions in factories to design better safety systems.

In a nutshell:
This paper gives us a high-definition, multi-lens camera for simulating how mixtures of different-sized particles behave in extreme, high-speed environments. It replaces the old "blurry" models with a sharp, robust tool that can finally see the difference between a grain of sand and a pebble when the world is shaking.

Technical Summary

1. Problem Statement

High-speed polydisperse granular multiphase flows are critical in natural phenomena (volcanic plumes, protoplanetary disks) and industrial/defense applications (propellants, dust explosions). While Lagrangian methods can handle polydispersity by tracking individual particles, they become computationally prohibitive for systems with vast particle numbers. Eulerian approaches are more efficient but traditionally rely on monodisperse assumptions (treating all particles as having the same size), which fails to capture critical physical phenomena driven by size distributions, such as:

• Size segregation: Different particle sizes accelerating at different rates under drag forces.
• Triboelectric charging and lightning generation.
• Complex flow structures: Vortex formation and shock interactions that vary by particle size.

Existing Quadrature-Based Moment Methods (QBMM) allow for continuous size distributions by transporting statistical moments (variance, skewness, etc.) rather than discrete bins. However, previous implementations for high-speed compressible flows relied on first-order numerical methods, which introduced excessive numerical dissipation, smearing out sharp interfaces, vortices, and contact surfaces essential for accurate high-speed flow simulation.

The Core Challenge: Developing a robust, high-order Eulerian numerical algorithm capable of simulating compressible, polydisperse multiphase flows with minimal numerical dissipation while maintaining stability in extreme conditions (e.g., near-packing limits, strong shocks).

2. Methodology

The authors extended their previous monodisperse high-order framework to a fully polydisperse model using the Generalized Quadrature Method of Moments (GQMOM).

A. Governing Equations

The model couples compressible gas equations (mass, momentum, energy, species) with a Generalized Population Balance Equation (GPBE) for the particle phase.

• Particle Phase: The GPBE tracks the evolution of the Particle Number Density Function (NDF) in phase space (position, velocity, mass, internal energy).
• Closure: The NDF is closed using Gaussian Quadrature, approximating the distribution with $N$ nodes (abscissae $m_\alpha$ and weights $f_\alpha$).
• Physics Included: The model accounts for drag, lift, convective heat transfer, particle-particle collisions (inelastic), frictional collisions (at high volume fractions), particle-fluid-particle pressure, and finite-size particle forces.
• Equation of State: Includes kinetic, collisional, and frictional pressure components. A specific impact-velocity-dependent coefficient of restitution is used to model particle deformation.

B. Numerical Scheme

The governing equations are solved using a Method of Lines approach with Strang operator splitting to separate hyperbolic (advection/wave) processes from source terms (drag, collisions, heat transfer).

1. Hyperbolic Solver (Advection):

   • Gas Phase: Solved using the HLLC approximate Riemann solver with 5th-order MUSCL reconstruction and a low-Mach correction.
   • Granular Phase: Since an analytical solution for the polydisperse Riemann problem is impossible, the authors transform the moment equations into a set of decoupled monodisperse granular equations at each quadrature node.
   • Granular Riemann Solver: A modified AUSM+-up solver is applied at each quadrature node to compute fluxes.
   • Reconstruction Strategy:
     • Weights and Abscissae: Reconstructed using 1st-order upwinding to strictly maintain realizability (preventing negative weights or masses).
     • Conditional Variables (Velocity, Temperature): Reconstructed using 5th-order WENO (Weighted Essentially Non-Oscillatory) schemes.
     • Robustness Checks: The algorithm dynamically degrades the order of reconstruction (from 5th to 3rd to 1st order) near "particle islands/lakes" (vacuum interfaces) or where particle diameter varies sharply within the stencil to prevent instability. A Rusanov flux is used as a fallback for degraded edges.

2. Source Term Solver:

   • Drag, heat transfer, and collision terms are integrated using analytical solutions (matrix exponentials) or stiff ODE solvers (VODE) where analytical solutions are not feasible.
   • The integration is split into sub-steps (Drag $\to$ Heat Transfer $\to$ Collisions $\to$ Friction) to improve stability.

3. Moment Inversion:

   • The Beta-GQMOM algorithm is employed to reconstruct the quadrature weights and abscissae from the transported moments. This reduces the number of required transport equations compared to standard QMOM and allows for accurate representation of continuous distributions.

3. Key Contributions

• High-Order Polydisperse Algorithm: First extension of QBMM to high-order (5th order) accuracy for high-speed compressible multiphase flows, significantly reducing numerical dissipation compared to previous first-order methods.
• Decoupled Riemann Solver Strategy: A novel approach to solving the polydisperse granular Riemann problem by treating each quadrature node as an independent monodisperse system, enabling the use of robust monodisperse solvers (AUSM) within a polydisperse framework.
• Dynamic Reconstruction Degradation: A robust strategy to handle "vacuum" regions and sharp gradients by automatically lowering the interpolation order near particle islands, lakes, or large diameter variations, preventing numerical blow-up.
• Comprehensive Physics Integration: Full coupling of drag, lift, heat transfer, inelastic collisions, friction, and finite-size effects within a high-order Eulerian framework.

4. Results and Validation

The algorithm was validated against a series of 1D and 2D test cases involving aluminum particles (H-2 and H-10 size distributions):

1. Linear Advection of Particle Curtain: Demonstrated the method's ability to preserve sharp phase discontinuities without diffusing the interface or disturbing constant pressure/temperature fields (error $\sim 10^{-9}$).
2. Dilute Polydisperse Shock Tube: Captured size segregation effects where smaller particles accelerated faster than larger ones behind the shock, creating a gradient in mean particle diameter. The results matched theoretical expectations and showed superior resolution compared to monodisperse models.
3. Shock-Particle Curtain Interaction: Simulated a Mach 2 shock hitting a dense curtain ($\alpha_p = 25\%$). Successfully captured the widening of the curtain and the smooth gradient of particle size due to differential acceleration.
4. Shock-Dilute Dust Layer: Modeled a Mach 1.6 shock interacting with a dust layer. The polydisperse model produced a single, diffused vortex (unlike the two distinct vortices in monodisperse models) and correctly predicted size segregation at the shock front.
5. Shock-Dense Dust Pile: Tested a Mach 3 shock on a dense pile ($\alpha_p = 47\%$). The model handled the "locking" of particles near the packing limit and the formation of granular clusters (streamers) due to non-linear drag and inelastic collisions.
6. Shock-Induced Dust Lifting: Validated against experimental data for a Mach 1.53 shock lifting a dense dust layer. The simulation accurately predicted the dust rise height and time history, capturing the size segregation where small particles remained in the boundary layer while larger particles traveled further.
7. Spherical Dispersal of Dense Particle Bed: An extreme test involving a high-pressure gas explosion dispersing a dense particle shell. The model captured granular shocks, radial expansion, and complex size segregation where smaller particles were trapped in the fireball while larger particles escaped.

5. Significance

This work represents a major advancement in computational multiphase flow. By enabling high-resolution, low-dissipation simulations of polydisperse flows, the method allows researchers to:

• Accurately predict size segregation, a phenomenon completely missed by monodisperse models but critical in volcanic flows, planetary formation, and industrial processing.
• Simulate extreme conditions (high pressures, dense packing, strong shocks) with robustness, making it suitable for defense applications (propellants, explosives) and safety analysis (dust explosions).
• Provide a computationally efficient alternative to Lagrangian tracking for systems with billions of particles, bridging the gap between detailed physics and practical engineering simulation.

The developed framework sets a new standard for Eulerian multiphase flow solvers, offering a pathway to include complex physical phenomena (combustion, breakage, phase change) in high-fidelity simulations of real-world polydisperse systems.