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

Abstract: This paper presents a high-order bound-preserving oscillation-eliminating discontinuous Galerkin (BP-OEDG) scheme for simulating gas-gas and gas-liquid two-phase flows governed by the Kapila five-equation model with the Tammann equation of state (EOS). The primary computational bottleneck arises from the severe CFL restriction imposed by the stiff $κ$-source term in the volume fraction equation. To circumvent this, we propose a novel operator-splitting strategy that decouples the system into a transport model and a stiff $κ$-source term. The former is discretized via a quasi-conservative DG method \cite{cheng2020quasi}, while the latter is resolved by an adaptive implicit strategy hybridizing the backward Euler and SDIRK2 methods. We rigorously prove that this implicit treatment is unconditionally BP, effectively removing the stiffness-induced stability constraints inherent in traditional explicit schemes. To further enhance precision, a velocity divergence reconstruction inspired by the Local Discontinuous Galerkin (LDG) method is integrated into the implicit solver. Furthermore, an OE limiter is employed to suppress spurious oscillations without characteristic decomposition, complemented by a BP limiter to ensure the BP property of partial densities, pressure, and volume fraction. Crucially, we prove that the proposed BP-OEDG scheme, integrated with the splitting strategy, strictly satisfies the Abgrall condition. Extensive numerical experiments, including challenging water-air shock-bubble interactions, demonstrate the superior robustness and efficiency of the method.

• The paper develops a high-order bound-preserving, oscillation-eliminating DG method for the Kapila five-equation model, using implicit treatment of the stiff volume-fraction source to avoid severe CFL restrictions.
• The scheme proves convexity of a pressure-positive admissible state set, preserves the Abgrall pressure–velocity equilibrium condition, and achieves second- and third-order convergence with DG(P1) and DG(P2), respectively.
• Numerical tests show stable results for extreme gas–liquid problems, including shock–bubble interactions, while Strang splitting remains limited to second-order interface accuracy and motivates future fully coupled IMEX methods.

This paper develops a high-order bound-preserving oscillation-eliminating discontinuous Galerkin (BP-OEDG) scheme for the Kapila five-equation model of compressible gas-gas and gas-liquid two-phase flows, closed by the Tammann equation of state (EOS). The central obstacle addressed is the severe CFL restriction induced by the stiff $\kappa$-source term $\kappa\nabla\cdot\mathbf{u}$ in the volume fraction equation, where $\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$ measures the compressibility contrast between the phases. The authors combine an operator-splitting strategy with an adaptive implicit solver for the stiff term, an oscillation-eliminating (OE) procedure, and a bound-preserving (BP) limiter, and they prove that the resulting fully discrete scheme satisfies the Abgrall condition for the Tammann EOS.

Admissible state space

The paper departs from prior DG work for two-medium flows by replacing the conventional admissible set $$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$ with the stricter set $\mathcal{G}^p$ that additionally requires positive mixture pressure $p>0$, appropriate for gas-gas and gas-liquid flows. A concrete counterexample in the appendix shows $\mathcal{G}$ is non-convex under the Tammann EOS: two individually admissible states (one with negative pressure but positive sound speed) can have an inadmissible convex combination with $c^2<0$. This is a notable claim, since convexity of the admissible set is typically the foundation of Zhang–Shu-type BP analysis. Theorem 1 proves $\mathcal{G}^p$ is convex, using the Cauchy–Schwarz inequality to establish concavity of the negative kinetic energy. This convexity underwrites the BP limiter design and is an independent structural contribution beyond the splitting strategy.

Operator splitting and the implicit $\kappa$-solver

The Kapila model is split via second-order Strang splitting into (i) the five-equation transport model, solved explicitly by the quasi-conservative DG method with Lax–Friedrichs fluxes, and (ii) the stiff ODE $\kappa\nabla\cdot\mathbf{u}$0, in which only the volume fraction evolves. The implicit solver is a backward Euler discretization in which $\kappa\nabla\cdot\mathbf{u}$1, the mixture parameters, and the pressure are all evaluated at the new time level; Theorem 2 proves existence of a bounded solution $\kappa\nabla\cdot\mathbf{u}$2 for any $\kappa\nabla\cdot\mathbf{u}$3 via a sign argument on an auxiliary function, with the root located by bisection at each quadrature point. The practical consequence is that the BP CFL constraint imposed by the volume fraction equation—previously shrinking the time step by orders of magnitude—is removed entirely, so the scheme runs at a standard CFL number.

To recover second-order temporal accuracy, an adaptive strategy hybridizes SDIRK2 with backward Euler: the first stage uses backward Euler over $\kappa\nabla\cdot\mathbf{u}$4 with $\kappa\nabla\cdot\mathbf{u}$5; if the explicitly extrapolated predictor remains in $\kappa\nabla\cdot\mathbf{u}$6, the second SDIRK2 stage is taken, otherwise the step locally degenerates to backward Euler. This preserves unconditional BP while retaining second-order accuracy wherever the solution is far from the bounds. The authors note that the second-order limit of Strang splitting is localized to the two-phase interface, since single-phase regions ($\kappa\nabla\cdot\mathbf{u}$7) bypass the implicit steps.

To restore optimal $\kappa\nabla\cdot\mathbf{u}$8-order spatial accuracy for the divergence term—direct differentiation of a DG($\kappa\nabla\cdot\mathbf{u}$9) velocity yields only $\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$0—the paper reconstructs $\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$1 through an LDG-style weak formulation with a purely upwind flux, and substitutes this $\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$2 into the implicit source term.

Limiters and preservation of the Abgrall condition

After each RK stage, an OE procedure exponentially damps high-order modal coefficients based on jump-based oscillation indicators, eliminating spurious oscillations without characteristic decomposition (a simplification relative to WENO or TVB limiters). A two-step BP limiter then enforces volume fraction and partial density bounds, followed by pressure positivity, via convex scaling toward the cell average; by Theorem 1 this keeps the solution in $\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$3 while preserving cell means.

The key theoretical result (Theorem 4) is that the fully discrete 2D scheme—quasi-conservative DG, implicit $\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$4-solver, SSP-RK time integration, OE operator, and BP limiter—satisfies the Abgrall condition (pressure and velocity equilibrium preserved at an isolated material interface). The proof exploits the linearity-invariance of both the OE procedure and the BP limiter: constant-pressure/equilibrium states are expressible as a linear identity $\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$5, which both operators preserve. The same argument applies to the implicit solver, which trivially satisfies the condition since $\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$6 freezes the state.

Numerical results

Convergence is verified on a 1D smooth gas-gas case with periodic boundaries. With $\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$7 at CFL 0.1, the DG($\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$8) scheme achieves clean second-order convergence and DG($\kappa = z_1 z_2(\nu_1-\nu_2)/\nu$9) achieves third order (e.g., $$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$0 error $$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$1 at 320 cells, order 2.99), confirming that the second-order splitting error remains subdominant to the spatial error.

Case	Purpose	Outcome
Gas-liquid isolated interface	Non-disturbing property	Pressure/velocity undisturbed at $$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$2 volume-fraction contrast
Double rarefaction	BP property	No negative pressure; Kapila and 5-eq models differ in $$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$3 due to compressibility
Gas-liquid Riemann problem	Robustness	Good agreement with fine-mesh reference
Extreme shock tube ($$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$4 ratio $$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$5, $$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$6 ratio 200)	Stiffness	Shock, contact, rarefaction captured stably
Air shock on helium bubble (Haas–Sturtevant)	Experimental validation	Bubble shapes agree at 62–983 $$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$7s
Water shock on air bubble ($$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$8)	Extreme gas-liquid case	Major features captured without spurious oscillations

The water-shock/air-bubble case is significant because such extreme gas-liquid computations are described as absent from the existing DG literature due to stiffness-induced CFL collapse; the splitting scheme runs them at standard CFL. The authors state this directly as a capability claim: the method "removes the stiffness-induced stability constraints inherent in traditional explicit schemes."

Limitations and open questions

The paper is candid that Strang splitting caps temporal accuracy at second order at the two-phase interface, and, more seriously, that in cases of extreme $$\mathcal{G}=\{z_1\rho_1>0,\ z_2\rho_2>0,\ c^2>0,\ z_1\in[0,1]\}$$9-stiffness the splitting may suffer local order degradation or fail to recover the correct physical solution, citing classical splitting analysis. The BP property is proven for the implicit substep, but the combined explicit transport step's BP property relies on the limiter rather than on a fully discrete CFL-type theorem for the coupled system. The authors identify development of a fully coupled IMEX-OEDG framework—avoiding splitting altogether—as the outstanding challenge and the focus of their ongoing work. Whether the adaptive SDIRK2/backward-Euler hybrid retains second-order accuracy uniformly in the stiffness parameter is not analyzed.

Conclusion

The paper delivers a rigorously supported high-order DG framework for the Kapila five-equation model that removes the stiffness-induced time-step restriction through an unconditionally BP implicit treatment of the $\mathcal{G}^p$0-source term, while provably maintaining the Abgrall condition and physical bounds under the Tammann EOS. Its identification of the non-convexity of the standard admissible set, and the convex alternative $\mathcal{G}^p$1, is a substantive observation for BP analysis of two-medium flows. The remaining open question is whether a fully coupled IMEX formulation can achieve the same robustness without the second-order splitting constraint at material interfaces.