A robust computational framework for the mixture-energy-consistent six-equation two-phase model with instantaneous mechanical relaxation terms

We present a robust computational framework for the numerical solution of a hyperbolic 6-equation single-velocity two-phase system. The system’s main interest is that, when combined with instantaneous mechanical relaxation, it recovers the solution of the 5-equation model of Kapila. Several numerical methods based on this strategy have been developed over the years. However, neither the 5- nor 6-equation model admits a complete set of jump conditions because they involve non-conservative products. Different discretizations of these terms in the 6-equation model exist. The precise impact of these discretizations on the numerical solutions of the 5-equation model, in particular for shocks, is still an open question to which this work provides new insights. We consider the phasic total energies as prognostic variables to naturally enforce discrete conservation of total energy and compare the accuracy and robustness of different discretizations for the hyperbolic operator. Namely, we discuss the construction of an HLLC approximate Riemann solver in relation to jump conditions. We then compare an HLLC wave-propagation scheme which includes the non-conservative terms, with Rusanov and HLLC solvers for the conservative part in combination with suitable approaches for the non-conservative terms. We show that some approaches for the discretization of non-conservative terms fit within the framework of path-conservative schemes for hyperbolic problems. We then analyze the use of various numerical strategies on several relevant test cases, showing both the impact of the theoretical shortcomings of the models as well as the importance of the choice of a robust framework for the global numerical strategy.

💡 Research Summary

The paper presents a comprehensive computational framework for solving the hyperbolic six‑equation single‑velocity two‑phase flow model, particularly when instantaneous mechanical relaxation is applied so that the model reduces to the classical five‑equation Kapila formulation. The authors begin by highlighting a fundamental difficulty: both the six‑equation and five‑equation models contain non‑conservative products, which prevents a complete set of jump conditions from being uniquely defined. Consequently, the design of robust Riemann solvers and the choice of discretization for the non‑conservative terms become critical, especially in the presence of shocks and strong rarefactions.

To address this, the paper investigates several discretization strategies for the hyperbolic operator. The first strategy couples a conventional approximate Riemann solver (either Rusanov or HLLC) for the conservative part of the system with dedicated treatments for the non‑conservative terms, following approaches previously used for diffusion operators. The second strategy adopts a wave‑propagation method in which the HLLC approximate Riemann solver is extended to include the non‑conservative contributions directly, as originally proposed for the six‑equation model with instantaneous relaxation. In both cases, the phasic total energies are taken as prognostic variables, guaranteeing discrete conservation of the mixture total energy without the need for additional correction steps that are required when phasic internal energies are used.

A key theoretical contribution is the reinterpretation of several non‑conservative discretizations as path‑conservative schemes. By expressing the non‑conservative products as integrals along prescribed paths in state space, the authors show that the resulting schemes satisfy a generalized form of the Rankine‑Hugoniot condition and inherit the stability properties of well‑known path‑conservative methods for hyperbolic systems. This insight explains why some discretizations are markedly more robust than others.

The authors then conduct an extensive numerical assessment on a suite of benchmark problems: (i) one‑dimensional shock–tube tests with large pressure jumps, (ii) two‑dimensional cylindrical explosion, (iii) problems involving pressure non‑equilibrium and finite‑rate mechanical relaxation, and (iv) cases with extreme rarefaction waves. Across all tests, the wave‑propagation HLLC scheme with path‑conservative treatment of the non‑conservative terms consistently preserves the physical bounds of the volume fraction (0 ≤ α ≤ 1), yields shock profiles that match the reference five‑equation solution, and avoids spurious oscillations. The Rusanov/HLLC combination also produces accurate results for many cases, but in highly non‑equilibrium situations it can generate slight overshoots or undershoots in the volume fraction and in the distribution of phasic energies, leading to noticeable differences in the post‑shock state.

Importantly, the paper demonstrates that the six‑equation model can be used as an independent solver rather than merely an auxiliary tool for the five‑equation model. By applying the instantaneous mechanical relaxation step after each hyperbolic update, the framework remains applicable when the relaxation rate μ is finite, opening the way to simulations of finite‑rate relaxation phenomena. The authors also integrate their methods into the SAMURAI finite‑volume library, showing that the approach is compatible with existing high‑performance computing infrastructures.

In conclusion, the study provides a rigorous analysis of the impact of non‑conservative term discretization on the solution of the six‑equation two‑phase model, proposes a robust path‑conservative HLLC wave‑propagation scheme, and validates its superiority through comprehensive testing. The resulting framework offers improved accuracy and stability for simulations involving shocks, phase transitions, and mechanical relaxation, and it lays a solid foundation for future extensions to multi‑dimensional, multi‑physics two‑phase flow problems.