A simple, fully-discrete, unconditionally energy-stable method for the two-phase Navier-Stokes Cahn-Hilliard model with arbitrary density ratios

The two-phase Navier-Stokes Cahn-Hilliard (NSCH) mixture model is a key framework for simulating multiphase flows with non-matching densities. Developing fully discrete, energy-stable schemes for this model remains challenging, due to the possible presence of negative densities. While various methods have been proposed, ensuring provable energy stability under phase-field modifications, like positive extensions of the density, remains an open problem. We propose a simple, fully discrete, energy-stable method for the NSCH mixture model that ensures stability with respect to the energy functional, where the density in the kinetic energy is positively extended. The method is based on an alternative but equivalent formulation using mass-averaged velocity and volume-fraction-based order parameters, simplifying implementation while preserving theoretical consistency. Numerical results demonstrate that the proposed scheme is robust, accurate, and stable for large density ratios, addressing key challenges in the discretization of NSCH models.

💡 Research Summary

The paper addresses the longstanding difficulty of constructing fully discrete, unconditionally energy‑stable numerical schemes for the two‑phase Navier‑Stokes‑Cahn‑Hilliard (NSCH) model when the two fluids have arbitrary, possibly very different, densities. While many energy‑stable methods exist for the equal‑density case, extending them to non‑matching densities is problematic because the density may become negative under phase‑field modifications, and the standard formulations involve a non‑divergence‑free velocity field that complicates the analysis.

To overcome these obstacles, the authors first reformulate the NSCH system in terms of the mass‑averaged velocity and a volume‑fraction‑based order parameter (φ). By linearly combining the mass‑balance equation with the phase‑field evolution equation they obtain an explicit expression for div v that includes a corrective term involving the mobility tensor. Substituting this expression back into the momentum equation yields an alternative but mathematically equivalent strong form (equations (10a)–(10d)). This reformulation has two crucial advantages: (i) the weighting functions used in the energy‑stability proof now belong to standard Sobolev spaces, and (ii) the kinetic‑energy term can be written with a positively extended density ρ⁺(φ), guaranteeing that the density entering the energy functional is always non‑negative.

The authors then derive a variational (weak) formulation (Lemma 5) that respects the underlying thermodynamic structure: mass is conserved, the total density is conserved, and the total free energy E(φ,v)=∫Ω