A HHO formulation for variable density incompressible flows where the density is purely advected

Imagine you are trying to simulate a swirling storm of two different liquids, like oil and water, or hot air and cold air, inside a computer. These liquids don't mix, but they push against each other, and their densities (how heavy they are for their size) change depending on where they are.

This paper presents a new, super-smart mathematical recipe (called HHO-ESDIRK) for a computer to solve this puzzle. The authors, Lorenzo Botti and Francesco Carlo Massa, have built a tool that is faster, more accurate, and more stable than previous methods.

Here is the breakdown using simple analogies:

1. The Problem: The "Heavy vs. Light" Dance

Think of the Rayleigh-Taylor Instability (the main test case in the paper) as a heavy blanket (dense fluid) sitting on top of a light pillow (light fluid). Gravity wants to pull the heavy blanket down, but the light pillow wants to float up. They start to mix in a chaotic, swirling dance.

The difficulty for computers is that:

The liquids must stay "incompressible" (you can't squeeze them into a smaller box).
The density of the liquid moves only with the flow (like a leaf floating down a river; it doesn't magically appear or disappear).
If the computer makes a tiny mistake, the simulation can explode or give nonsense results.

2. The Solution: The "Hybrid High-Order" (HHO) Method

The authors use a method called HHO. Imagine you are trying to paint a picture of a complex landscape.

Old methods were like painting with a single giant brushstroke for the whole picture, or using tiny, disconnected tiles that didn't quite line up.
The HHO method is like having a team of artists. They paint the inside of each small tile (the cell) and the edges of the tiles (the skeleton) separately, but they constantly talk to each other to make sure the picture looks seamless.

Why is this special?

Volume Conservation: The method ensures that the total amount of "stuff" in the box never changes. It's like a leak-proof bucket; no water is lost or gained by accident.
Pure Advection: Because the bucket is leak-proof, the density (the "color" of the fluid) just gets carried along by the current without getting smeared or distorted. It's like a perfect conveyor belt.
Pressure Robustness: In fluid physics, pressure is often the "boss" that causes errors. This method makes the velocity (speed) and density calculations immune to pressure mistakes. It's like having a car where the engine performance doesn't drop just because the speedometer is slightly off.

3. The Time Machine: ESDIRK

Simulating fluids happens over time. You can't just jump from "now" to "later"; you have to take steps.

The authors use a specific type of time-stepping called ESDIRK. Think of this as a very careful hiker climbing a steep mountain. Instead of taking giant, risky leaps, the hiker takes small, calculated steps, checking their footing at every stage before moving forward.
This allows the computer to take larger steps in time without losing accuracy, making the simulation run much faster.

4. The "Magic Trick": Static Condensation

One of the biggest headaches in these simulations is memory. The computer needs to remember millions of numbers to solve the equations.

The Trick: The HHO method uses a technique called Static Condensation. Imagine you are solving a giant jigsaw puzzle. Instead of trying to hold every single piece in your hands at once, you solve the pieces inside each small section of the puzzle first, and only keep the "edge pieces" to connect the sections together.
The Result: This drastically reduces the amount of memory the computer needs. It's like shrinking a massive library down to a few key index cards without losing the story.

5. The Results: What Did They Prove?

The authors tested their new recipe on two main scenarios:

The Smooth Test (Kovasznay Flow): They checked if the math worked perfectly on a known, smooth flow. It did. The computer got the answer right down to the last decimal place, and the error dropped rapidly as they made the grid finer.
The Chaotic Test (Rayleigh-Taylor Instability): They simulated the heavy blanket vs. light pillow scenario.
- Low Density Difference: They showed that using a "high-order" method (very detailed, complex math) on a coarse (low-resolution) grid gave better results than a "low-order" method on a super-fine grid. This is huge because it saves massive amounts of computing power.
- High Density Difference: When the fluids were very different (heavy vs. very light), the simulation can become unstable. They proved that their method, even at the simplest level, keeps the density values realistic (it doesn't accidentally calculate negative density, which is physically impossible).

The Bottom Line

This paper introduces a new, highly efficient way to simulate mixing fluids. It combines a clever spatial grid (HHO) with a smart time-stepping strategy (ESDIRK) and a memory-saving trick (Static Condensation).

In everyday terms: It's like upgrading from a bicycle with a flat tire to a high-speed electric scooter that can navigate rough terrain without losing its balance, all while using less battery power. This allows scientists to simulate complex fluid mixtures (like oil spills, weather patterns, or industrial mixing) faster and more accurately than ever before.

Here is a detailed technical summary of the paper "A HHO formulation for variable density incompressible flows where the density is purely advected" by Botti and Massa.

1. Problem Statement

The paper addresses the numerical simulation of variable-density incompressible flows, specifically focusing on mixtures of immiscible fluids. The governing equations are the incompressible Navier–Stokes equations coupled with a mass conservation equation where density is treated as a purely advected scalar.

Key challenges in this regime include:

Exact Incompressibility: Ensuring the velocity field is divergence-free to satisfy the volume conservation constraint.
Pressure-Robustness: Ensuring that errors in the velocity and density fields are not polluted by the accuracy of the pressure approximation, particularly when irrotational body forces are present.
Convection Dominance: Maintaining stability and preventing spurious oscillations in high Reynolds number flows or at sharp interfaces (e.g., Rayleigh-Taylor instability).
High-Order Accuracy: Achieving high-order convergence in both space and time while managing computational costs.

2. Methodology

The authors propose a Hybrid High-Order (HHO) formulation coupled with Explicit Singly Diagonal Implicit Runge-Kutta (ESDIRK) time integration.

Spatial Discretization (HHO)

Hybrid Spaces: The method utilizes hybrid unknowns defined on both mesh cells ( $T$ $T$ ) and mesh faces ( $F$ $F$ ).
- Velocity: Cell unknowns are polynomials of degree $k+1$ , and face unknowns are degree $k$ (or $k+1$ on Neumann boundaries).
- Density: Cell and face unknowns are polynomials of degree $k$ .
- Pressure: Cell unknowns are degree $k$ , and face unknowns are degree $k+1$ .
H(div)-Conformity: A key feature is the reconstruction of an $H(\text{div})$ -conforming velocity field from non-conforming polynomial spaces. This ensures that the discrete velocity is pointwise divergence-free, allowing the mass conservation equation to be formulated as a pure advection of density ( $\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0$ ).
Stabilization:
- Momentum: Uses skew-symmetric formulations for convection terms with upwind stabilization. Two variants are proposed: HHO- $\pi^k$ (using degree $k$ projections for advected velocity, provably stable) and HHO- $\pi^{k+1}$ (using degree $k+1$ projections, offering higher accuracy in convection-dominated regimes).
- Mass Transport: Employs upwind fluxes and a specific time-derivative stabilization term to ensure robustness when normal velocity vanishes on faces.
Boundary Conditions: Weakly imposed via Nitsche's method and appropriate flux formulations.

Temporal Discretization (ESDIRK)

The system is treated as a set of Differential-Algebraic Equations (DAEs) of index 2 due to the incompressibility constraint.
ESDIRK Schemes: High-order accurate, implicit time-stepping schemes are used. These are specifically designed to handle the algebraic constraints of incompressible flows efficiently.
Algebraic Structure: The method employs static condensation to eliminate cell unknowns locally, reducing the global system size to only face unknowns. This significantly reduces memory footprint and improves the sparsity of the Jacobian matrix.
Solver Strategy: Iterative solvers (FGMRES) are accelerated using $p$ -multigrid preconditioners, leveraging the hierarchical structure of the polynomial spaces.

3. Key Contributions

Exact Volume Conservation & Pure Advection: The formulation enforces exact volume conservation cell-by-cell up to machine precision. Consequently, the density is purely advected without artificial diffusion or compression artifacts inherent in some other methods.
Pressure-Robustness: The method is proven to be pressure-robust. Irrotational body forces affect only the pressure field, leaving velocity and density errors insensitive to pressure approximation errors.
Stability Analysis: The authors provide a discrete stability proof showing that the scheme is dissipative for homogeneous Dirichlet boundary conditions, utilizing skew-symmetric forms and specific stabilization terms.
Bounds Preservation at Low Order: The $k=0$ variant is shown to be bounds-preserving, preventing non-physical negative densities in high-density-ratio scenarios.
Efficiency: Through static condensation and $p$ -multigrid, the method achieves high-order accuracy with a reduced memory footprint compared to standard DG methods.

4. Numerical Results

The method was validated through three distinct test cases:

Kovasznay Flow (Steady, Constant Density):
- Demonstrated optimal spatial convergence rates ( $k+1$ for pressure, $k+2$ for velocity in $L^2$ , $k+1$ for gradients).
- Confirmed that the density remains constant (pure advection) regardless of discretization parameters.
- Showed that HHO- $\pi^{k+1}$ offers superior accuracy for low-order ( $k=0$ ) approximations compared to HHO- $\pi^k$ .
Guermond and Quartapelle Manufactured Solution (Unsteady, Variable Density):
- Validated both spatial and temporal convergence rates.
- Confirmed pressure-robustness: Velocity and density errors remained independent of pressure errors.
- Demonstrated that high-order ESDIRK schemes (3-stage and 6-stage) achieve their theoretical temporal convergence rates (2nd and 4th order, respectively).
Rayleigh-Taylor Instability (RTI):
- Low Atwood Number ( $At=0.5$ ): Simulated at $Re=1,000$ and $Re=5,000$ . High-order ( $k=6$ ) solutions on coarse meshes matched low-order ( $k=1$ ) solutions on fine meshes, demonstrating significant computational efficiency. However, at high $Re$ , spurious oscillations appeared without oscillation control strategies.
- High Atwood Number ( $At=0.75$ ): Simulated at $Re=1,000$ with a density ratio of 7:1. The $k=0$ formulation successfully preserved density bounds (no negative densities), proving stability for high-density ratios. While spatial accuracy was limited at $k=0$ , the main plume structure was correctly captured.

5. Significance

This work represents a significant advancement in the simulation of variable-density incompressible flows. By combining the HHO framework with ESDIRK time integration, the authors have created a method that:

Unifies Accuracy and Stability: It achieves high-order accuracy while maintaining the stability required for convection-dominated and high-density-ratio flows.
Enables Efficient High-Order Simulation: The use of static condensation and $p$ -multigrid makes high-order simulations computationally feasible, offering a "best of both worlds" scenario where coarse high-order meshes can outperform fine low-order meshes.
Provides a Robust Foundation for Multiphase Flow: The ability to handle pure advection of density with exact incompressibility makes this formulation a strong candidate for future extensions to complex multiphase systems, such as the Cahn-Hilliard-Navier-Stokes (CHNS) system.

The paper concludes that the proposed ESDIRK-HHO formulation is a powerful tool for simulating complex fluid mixtures, offering a unique balance of physical fidelity (pressure-robustness, exact volume conservation) and computational efficiency.