A bound-preserving oscillation-eliminating discontinuous Galerkin scheme for compressible two-phase flow

This paper proposes a high-order, bound-preserving oscillation-eliminating discontinuous Galerkin scheme for compressible two-phase flows that overcomes severe stiffness-induced CFL restrictions through a novel operator-splitting strategy with an adaptive implicit solver, while rigorously ensuring stability, accuracy, and adherence to the Abgrall condition.

The Gist

Imagine you are trying to simulate a high-speed collision between two very different fluids, like a shockwave in water smashing into a bubble of air. In the world of computer simulations, this is a nightmare. The fluids behave differently, they squash and stretch at different rates, and the math governing their interaction is incredibly "stiff."

Think of "stiffness" here like trying to drive a car where the brakes are stuck on the floor. If you try to move forward even a tiny bit (a small time step in the simulation), the brakes fight back so hard that the car might flip over or the engine might explode. In computer terms, this forces the simulation to take steps so incredibly small that it would take years to simulate a split second of real time.

This paper introduces a new, smarter way to drive that car. Here is the breakdown of their solution using simple analogies:

1. The Problem: The "Stiff" Brake

The authors are working with a specific set of rules (the Kapila five-equation model) that describes how two fluids mix and move. The trouble comes from one specific rule (the $\kappa$-source term) that handles how the fluids compress. When a shockwave hits the boundary between water and air, this rule goes into overdrive.

If the computer tries to solve everything all at once (the traditional way), it gets stuck. To keep the math from breaking, it has to slow down the simulation time so drastically that the calculation becomes impossible.

2. The Solution: The "Split-Second" Strategy

The authors propose a clever trick called Operator Splitting. Imagine you are trying to bake a cake while simultaneously fixing a leaking pipe. Doing both at the exact same moment is chaotic and likely to fail. Instead, you do them in separate, focused steps:

• Step A: Fix the pipe (solve the "stiff" compression part).
• Step B: Bake the cake (solve the movement and flow part).

By separating these two tasks, the computer can handle the "leaky pipe" (the stiff math) using a special, slow-and-steady implicit method that never breaks, and then handle the "baking" (the flow) using a fast, high-precision method.

3. The "Bound-Preserving" Safety Net

In these simulations, numbers represent physical things like density and pressure. If the math goes wrong, the computer might calculate that air has a negative density or that a bubble has 150% of its volume (which is impossible). This causes the simulation to crash.

The authors built a Bound-Preserving (BP) limiter. Think of this as a bouncer at a club. If a number tries to leave the "safe zone" (e.g., a volume fraction trying to go above 100% or below 0%), the bouncer immediately kicks it back inside the safe zone. This ensures the simulation never produces "nonsense" physics, even when things get chaotic.

4. The "Oscillation-Eliminating" Shock Absorber

When a shockwave hits a bubble, it creates sharp edges and ripples. Standard math often creates fake, jagged "ghost waves" (oscillations) around these sharp edges, making the picture look noisy and wrong.

The authors use an Oscillation-Eliminating (OE) technique. Imagine driving over a bumpy road. A standard car might bounce wildly. This new method acts like a high-tech suspension system that smooths out the ride without losing the detail of the bumps. It removes the fake noise while keeping the real physics sharp, and it does this without needing to do complex, slow calculations to figure out the direction of the waves.

5. The Result: A Smooth, Fast Ride

The authors tested their new method on some very difficult scenarios:

• Shock hitting a helium bubble: Like a sonic boom hitting a soap bubble.
• Water shock hitting an air bubble: A massive underwater explosion hitting a pocket of air.

In these tests, their method was able to run fast (using standard time steps) without crashing, while keeping all the numbers physically realistic. It captured the complex shapes of the bubbles and the shockwaves with high precision, proving that you can simulate these extreme events without the computer getting stuck in "slow motion."

In summary: The paper presents a new mathematical engine that splits a difficult problem into manageable chunks, uses a safety net to keep numbers realistic, and smooths out the noise. This allows computers to simulate violent collisions between different fluids quickly and accurately, solving a problem that previously required impossible amounts of computing power.

Technical Summary

1. Problem Statement

The paper addresses the numerical simulation of compressible gas-gas and gas-liquid two-phase flows governed by the Kapila five-equation model. This model is a reduced form of the Baer-Nunziato model, assuming mechanical and thermal equilibrium (equal velocities and pressures) between phases.

Key Challenges:

• Stiffness of the $\kappa$-source term: The volume fraction equation contains a non-conservative source term ($\kappa \nabla \cdot \mathbf{u}$) that accounts for the compressibility differences between fluids. When shock or rarefaction waves interact with material interfaces, this term becomes extremely stiff.
• CFL Restriction: Traditional explicit time-integration schemes require time steps restricted by the stiffness of the $\kappa$-term, often reducing the allowable time step by several orders of magnitude compared to the standard CFL condition. This makes simulating extreme cases (e.g., water-air shock-bubble interactions) computationally prohibitive.
• Physical Bounds and Oscillations: High-order Discontinuous Galerkin (DG) methods are prone to generating spurious oscillations near discontinuities and can violate physical bounds (e.g., negative pressure, volume fractions outside $[0,1]$), leading to numerical breakdown.
• Abgrall Condition: Maintaining pressure and velocity equilibrium at isolated material interfaces is critical to avoid non-physical artifacts, a property known as the Abgrall condition.

2. Methodology

The authors propose a Bound-Preserving Oscillation-Eliminating Discontinuous Galerkin (BP-OEDG) scheme integrated with a novel operator-splitting strategy.

A. Operator-Splitting Strategy

To circumvent the stiffness-induced stability constraints, the Kapila five-equation model is decoupled into two sub-problems:

1. Transport Model (Five-equation transport): Contains the conservative laws for mass, momentum, and energy, plus the advection of volume fraction (without the stiff source).
   • Discretization: Solved using a quasi-conservative DG method (based on Zhang et al. [8]) that satisfies the Abgrall condition.
2. Stiff Source Term ($\kappa$-term): Contains only the stiff source term for the volume fraction equation.
   • Discretization: Solved using an adaptive implicit strategy.

B. Implicit Solver for the Stiff Term

• Hybrid Implicit Scheme: The authors propose a hybrid of the Backward Euler (first-order, unconditionally stable) and SDIRK2 (Second-order Diagonally Implicit Runge-Kutta) methods.
• Adaptive Mechanism: The scheme initiates with a Backward Euler step. If the predicted solution falls within the physical bounds $[0, 1]$, it proceeds to the second stage of SDIRK2 to maintain second-order temporal accuracy. If the prediction violates bounds, it degenerates locally to a robust Backward Euler step.
• Unconditional BP: The authors rigorously prove that this implicit treatment is unconditionally bound-preserving (BP), effectively removing the stiffness-induced time-step restriction.
• Velocity Divergence Reconstruction: To restore optimal high-order spatial accuracy (which is lost when directly differentiating DG polynomials), a Local Discontinuous Galerkin (LDG) inspired reconstruction is used to compute $\nabla \cdot \mathbf{u}$ within the implicit solver.

C. Limiters

• Oscillation-Eliminating (OE) Limiter: An OE procedure is applied after each Runge-Kutta stage. Unlike traditional limiters (e.g., WENO, TVD) that require characteristic decomposition, the OE limiter suppresses spurious oscillations without decomposing the system, significantly reducing computational complexity for multiphase flows.
• Bound-Preserving (BP) Limiter: A scaling limiter is applied to ensure partial densities, pressure, and volume fractions remain within the strict physical admissible set $G_p = \{U : z_1\rho_1 > 0, z_2\rho_2 > 0, p(U) > 0, z_1 \in [0, 1]\}$.

D. Theoretical Proofs

• Convexity of Admissible Set: The paper proves that the set $G_p$ is convex under the Tammann Equation of State (EOS), whereas the looser set $G$ (based on sound speed squared) is non-convex.
• Abgrall Condition Preservation: The authors prove that the entire framework (splitting, implicit solver, OE limiter, and BP limiter) strictly satisfies the Abgrall condition, ensuring that constant pressure and velocity are preserved at isolated material interfaces.

3. Key Contributions

1. Novel Operator-Splitting DG Framework: First application of a high-order DG method to the gas-liquid Kapila model that successfully circumvents the extreme CFL restrictions caused by the stiff $\kappa$-source term via operator splitting.
2. Unconditionally BP Implicit Solver: Development of an adaptive implicit solver (hybrid Backward Euler/SDIRK2) that guarantees the volume fraction remains in $[0,1]$ regardless of the time step size.
3. Oscillation Elimination without Characteristic Decomposition: Implementation of an OE limiter that effectively suppresses oscillations in complex multiphase systems without the expensive and complex characteristic decomposition required by traditional limiters.
4. Rigorous Theoretical Guarantees: Proof of the convexity of the physical admissible set $G_p$ and the preservation of the Abgrall condition for the fully discrete scheme.
5. High-Order Accuracy: Restoration of optimal $(K+1)$-th order accuracy for the divergence term via LDG-inspired reconstruction within the implicit solver.

4. Numerical Results

The method was validated through extensive 1D and 2D benchmark cases:

• Smooth Gas-Gas Flow: Demonstrated second-order (P1) and third-order (P2) convergence rates, matching reference solutions.
• Isolated Interface: Confirmed the non-disturbing property of pressure and velocity fields at a gas-liquid interface.
• Double Rarefaction & Riemann Problems: Showed the scheme's ability to handle strong rarefaction waves and high-pressure ratios while maintaining physical bounds (no negative pressure).
• Shock Tube: Successfully captured shocks, contact discontinuities, and rarefaction waves in a high-pressure gas-liquid shock tube.
• Shock-Bubble Interactions:
  • Air shock hitting Helium bubble: Results matched experimental data (Hass & Sturtevant) with high fidelity and no spurious oscillations.
  • Water shock hitting Air bubble: Successfully simulated the extreme stiffness of water-air interactions, a case often computationally prohibitive for explicit DG methods, demonstrating superior robustness.

5. Significance

This work represents a significant advancement in computational fluid dynamics for multiphase flows. By combining operator splitting with high-order DG methods, it resolves the long-standing trade-off between numerical stability (handling stiffness) and computational efficiency (allowing larger time steps). The ability to simulate extreme gas-liquid interactions (like shock-bubble collisions) with high-order accuracy and strict physical bounds makes this scheme a powerful tool for applications in aerospace, energy, and environmental sciences where such complex multiphase dynamics occur.

The paper acknowledges that while the splitting method introduces local order degradation at interfaces, it is sufficient for the targeted extreme compression dynamics. Future work aims to develop a fully coupled IMEX-OEDG framework to eliminate splitting errors entirely.