A family of second order, linear, unconditionally stable methods for the Cahn-Hilliard-Navier-Stokes equations

Technical Summary: A Family of Second-Order, Linear, Unconditionally Stable Methods for the Cahn-Hilliard-Navier-Stokes Equations

Problem Statement
The accurate numerical simulation of complex multiphase fluid dynamics, specifically those modeled by the Cahn-Hilliard-Navier-Stokes (CHNS) equations, presents significant challenges due to the mathematical stiffness of the fourth-order Cahn-Hilliard operator and its nonlinear coupling with the incompressible Navier-Stokes equations. Standard explicit time-stepping schemes face severe time-step restrictions (e.g., $\Delta t \sim O(\Delta x^4)$) to maintain stability. While fully implicit or convex-splitting methods offer unconditional energy stability, they often require solving large, coupled nonlinear systems at each time step, increasing computational cost. There is a need for algorithms that combine second-order temporal accuracy, unconditional stability, and linear solvability without sacrificing the physical conservation properties of the system.

Methodology
The authors propose a family of second-order, linear, unconditionally stable implicit-explicit (IMEX) methods for matched-density two-phase flows. The semi-discrete scheme integrates three primary components:

1. Extrapolation of Nonlinear Terms: Nonlinear convective and coupling terms are treated explicitly using extrapolated quantities (denoted as $H_{n+\theta}$), avoiding the need for nonlinear iterations.
2. Auxiliary Variable Formulation: The nonlinear free-energy term is handled via an auxiliary variable $q$ (derived from the double-well potential $F(\phi)$), reformulating the problem to maintain linearity.
3. Temporal-Curvature Regularization: A stabilization mechanism controlled by a parameter $\epsilon \geq 0$ is introduced. This involves a specific interpolation operator $J^\epsilon_{n+\theta}$ that incorporates the discrete temporal curvature of the solution, a concept adapted from previous work on Navier-Stokes regularization.

The resulting algorithm (Algorithm 2.1) requires solving only linear systems with constant coefficients at each time step. The scheme is parameterized by $\theta \in (1/2, 1]$ and $\epsilon \geq 0$.

Key Contributions and Theoretical Results

• Unconditional Stability: The paper establishes a rigorous discrete energy estimate proving that the proposed scheme is unconditionally long-time stable for $\theta \in (1/2, 1]$ and $\epsilon \geq 0$. The proof utilizes symmetric positive definite matrices to define discrete norms and demonstrates that the scheme satisfies a discrete energy-dissipation law without time-step restrictions.
• Linear Solvability: Unlike fully implicit or convex-splitting approaches that may require nonlinear solvers, this method yields linear systems at every time step, significantly reducing computational complexity.
• Second-Order Accuracy: Theoretical analysis and numerical experiments confirm that the method retains second-order temporal accuracy.

Numerical Results
The authors validate the method through a series of benchmark computations:

• Convergence Analysis: Using the Method of Manufactured Solutions (MMS), the scheme demonstrates approximately second-order temporal convergence for velocity, pressure, and phase field variables across different finite element configurations ($P2-P1-P2-P2$ and $P3-P2-P3-P3$).
• Robustness and Regularization: In the two-phase lid-driven cavity problem, the study examines the interplay between $\theta$ and $\epsilon$. It is observed that for $\theta$ values close to the Crank-Nicolson limit (e.g., $\theta=0.51$) with $\epsilon=0$, the scheme can become unstable due to weakly damped temporal oscillations. The introduction of positive curvature regularization ($\epsilon > 0$) effectively damps these oscillations, restoring boundedness and robustness even for larger time steps.
• Physical Benchmarks:
  • Spinodal Decomposition: The method accurately captures the energy dissipation and mass conservation over long-time simulations ($T=100$), showing effective phase separation and coarsening.
  • Droplet Relaxation: The algorithm successfully resolves surface-tension-driven interface motion, transitioning a square droplet to a circular equilibrium while preserving mass.
  • Lid-Driven Cavity & Rayleigh-Taylor Instability: The scheme handles strong shear, interfacial deformation, and density-driven instabilities (including Kelvin-Helmholtz roll-up in low-viscosity regimes) without visible spurious oscillations.

Significance
The paper claims that the proposed family of methods offers a practical and efficient alternative to fully nonlinear discretizations for CHNS systems. By decoupling the nonlinearities through extrapolation and auxiliary variables while stabilizing the temporal discretization via curvature regularization, the method achieves a balance of second-order accuracy, unconditional stability, and linear computational cost. The authors emphasize that the scheme preserves the principal physical structures of the model, including energy dissipation and mass conservation, making it suitable for simulating complex interfacial flows in various regimes.