Source-linked AI summary
A family of second order, linear, unconditionally stable methods for the Cahn-Hilliard-Navier-Stokes equations
Daozhi Han, Nan Jiang, Jonah H. Nissan, Sayantan Sarkar
TL;DR
The paper tackles restrictive time stepping and solver costs in CHNS simulations while preserving discrete energy dissipation. It proposes a linear second-order scheme with extrapolation, auxiliary variables, and curvature regularization, proving unconditional long-time stability and reporting approximately second-order convergence with conservation and dissipation in benchmark flows. The reported scope also includes robustness behavior and representative interfacial-flow simulations.
Problem
The stiff fourth-order Cahn-Hilliard operator and nonlinear fluid coupling can make explicit time stepping severely restrictive, while energy-stable alternatives may require nonlinear solves.
Method
The method combines extrapolated nonlinear terms, an auxiliary-variable free-energy reformulation, and temporal-curvature regularization in a fully linear second-order semi-discrete scheme.
Results
The scheme satisfies a discrete energy estimate guaranteeing unconditional long-time stability and requires only linear solves at each time step.
Takeaways & Limitations
Numerical experiments report approximately second-order temporal convergence, accurate mass conservation, energy dissipation, and simulations of representative interfacial flows.
Takeaways & Limitations
The interpretation of earlier loss of boundedness at smaller time steps is based on observed numerical behavior, not a rigorous instability mechanism established by the analysis.
Abstract
from arXiv · showhide
We present a family of second-order, linear, unconditionally stable implicit-explicit (IMEX) methods for the Cahn-Hilliard-Navier-Stokes (CHNS) equations modeling matched-density two-phase flows. The proposed semi-discrete scheme combines extrapolation of the nonlinear terms with an auxiliary-variable formulation of the nonlinear free-energy term and a temporal-curvature regularization controlled by a parameter $ε$. We establish a discrete energy estimate showing unconditional long-time stability of the method for $θ\in(1/2,1]$ and $ε\geq0$. The resulting scheme requires only linear solves at each time step. Numerical experiments demonstrate approximately second-order temporal convergence and examine mass conservation, energy dissipation, numerical robustness, and several representative interfacial-flow problems, including spinodal decomposition, droplet shape relaxation, two-phase lid-driven cavity flow, and Rayleigh-Taylor instability.
1. Introduction.
The paper addresses the difficulty of simulating strongly coupled CHNS flows by proposing a unified family of second-order, linear, unconditionally stable methods. The approach combines extrapolated nonlinear terms, auxiliary-variable free-energy treatment, and tunable temporal-curvature regularization.
- Diffuse-interface models represent material interfaces through smooth transitions and naturally accommodate droplet coalescence and pinch-off.
- Explicit treatment of the fourth-order Cahn-Hilliard operator can impose severe restrictions, as strong as Δt ∼ O(Δx4).
- Existing alternatives trade computational structure against stability or solver cost: fully implicit and convex-splitting schemes may require nonlinear systems or iterations.
- The proposed family targets second-order temporal accuracy, unconditional stability, linear nonlinear couplings, and tunable regularization in one framework.
- The fully linear algorithm uses extrapolated stiff couplings and temporal-curvature regularization, proves a discrete energy-dissipation law without a time-step restriction, and requires only linear solves.
- Numerical experiments examine temporal convergence, energy behavior, robustness, conservation, and representative multiphase-flow benchmarks.
2. The model and the numerical scheme.
The paper formulates matched-density CHNS flow for two immiscible incompressible fluids and introduces a second-order semi-discrete family using an auxiliary free-energy variable and curvature regularization.
- The model considers two immiscible incompressible fluids in a bounded Lipschitz domain with matched density set to unity.
- The phase field uses a Ginzburg-Landau free energy balancing mixing interactions against double-well separation energy.
- Competition between these interactions produces a diffusive interface whose thickness is proportional to η.
- The formulation imposes no-slip velocity and homogeneous normal-derivative conditions for the phase field and chemical potential.
- The continuous CHNS system satisfies an energy law under gravity as the only external forcing.
- The scheme introduces q = (ϕ2 − 1)/η2 so that the nonlinear free-energy derivative is represented as f(ϕ) = ϕq, then proposes a second-order linear family with ε ≥ 0.
3. Unconditional Stability.
The stability analysis establishes unconditional long-time stability for the proposed semi-discrete method through discrete norms, energy estimates, and cancellation of incompressibility-related terms.
- Theorem 3.2 states that the method is unconditionally long-time stable.
- The proof combines inner-product identities, summation over time levels, and lower bounds from Lemma 3.1 to derive the discrete energy estimate.
- Convection and pressure terms vanish through incompressibility and the homogeneous velocity boundary condition.
- The estimate includes weighted norms of q and u from the underlying three-level formulation.
- The numerical implementation uses Taylor–Hood P2−P1 elements for velocity and pressure and P2−P2 elements for phase and chemical-potential variables.
4. Numerical experiments.
Temporal convergence is examined with a Method of Manufactured Solutions on the unit square, using artificial forcing to make a prescribed analytical solution exact.
- The convergence test applies the Method of Manufactured Solutions on Ω = [0, 1]2.
- Artificial forcing terms are added to the continuous CHNS equations so prescribed functions form an exact analytical solution.
4.1. Convergence Analysis: Temporal Accuracy.
The temporal-convergence tests use refined spatial meshes, two mixed finite-element configurations, and successively halved time steps. Across both discretizations, the observed rates are close to second order for velocity, pressure, and phase field.
- Test setup: The simulations advance to T = 1.0 with exact-profile initial and boundary data and fixed parameters including θ = 0.8 and ϵ = 10−5.The exact velocity and phase-field profiles are prescribed for the convergence tests.
- Test setup: Two mixed finite-element configurations are tested: P2 −P1 −P2 −P2 with degree k = 2 and P3 −P2 −P3 −P3 with degree k = 3.
- Temporal refinement: The spatial mesh is fixed at Nx = 128 while the time step is halved from ∆t = 0.1 to ∆t = 0.0125.Discrete L2 errors are reported for velocity, pressure, and phase field.
- Results: Observed temporal-convergence rates are close to second order for velocity, pressure, and phase field in the P2-based tests.The corresponding P2 temporal-error decay is consistent with approximately second-order accuracy over the tested range.
- Results: The P3-based temporal-error decay is also consistent with approximately second-order temporal accuracy over the tested range.
4.2. Robustness of the Scheme.
The proposed scheme retains approximately second-order temporal behavior while improving robustness through curvature regularization and the dissipative parameter θ. In the lid-driven cavity tests, unregularized weakly dissipative settings can lose boundedness, whereas moderate regularization keeps energy and enstrophy bounded.
- Temporal Cauchy Convergence Test: The measured interior Cauchy rates are close to second order for most configurations and for the velocity, phase, and pressure variables.Moderate curvature regularization and extrapolated nonlinear terms show no evident first-order temporal contamination over the tested range.
- Temporal Cauchy Convergence Test: Very large regularization parameters, including ϵ = 5ν and ϵ = 10ν, produce irregular apparent Cauchy rates that should not be interpreted as genuine high-order convergence.Strong regularization can make successive differences extremely small, increasing sensitivity to cancellation, round-off, and pre-asymptotic effects.
- Dependence of Performance on ϵ and θ: The θ = 0.51, ϵ = 0 configuration loses boundedness for all three displayed time steps, with instability occurring earlier in physical time as ∆t is reduced.Thus, time-step refinement alone does not restore robustness in this weakly dissipative regime.
- Dependence of Performance on ϵ and θ: Positive curvature regularization keeps energy and enstrophy bounded over the tested interval, including for relatively large time steps.The results suggest that regularization damps temporal oscillatory components associated with loss of robustness in the unregularized case.
- Dependence of Performance on ϵ and θ: Increasing θ increases numerical dissipation, while ϵ provides a separate curvature-regularization mechanism that improves robustness near the weakly dissipative Crank–Nicolson limit.The θ = 1 BDF2 limit is substantially more robust even when ϵ = 0.
- Spinodal Decomposition: In spinodal decomposition, phase mass remains constant to numerical precision and computed energy decreases monotonically for each tested time step.The benchmark examines long-time behavior, energy dissipation, and phase-mass conservation.
4.3. Spinodal Decomposition, Energy Dissipation, and Mass Conservation.
The spinodal-decomposition benchmark tests long-time dissipation and phase-mass conservation. Energy decreases while phase mass remains constant to numerical precision, with fine-step energy evolution effectively converged.
- The simulation uses a unit-square domain with random-noise initialization, zero velocity and pressure, no-slip walls, and parameters ν = 0.1, η = 0.02, λ = 0.001, M = 0.01, ϵ = 10−5.
- For ∆t = 0.05, 0.005, and 0.001, total phase mass remains constant at approximately 0.20 while computed energy decreases monotonically.
- The energy curves for ∆t = 0.005 and ∆t = 0.001 are visually almost indistinguishable, indicating time-step convergence at the plotted resolution.
- During 0 ≤t ≤20, the initially mixed state rapidly separates into an interconnected network of phase domains.
- During 40 ≤t ≤100, coarsening continues as curved structures relax, thin connections break, and smaller domains merge into larger ones while energy decreases and phase mass remains conserved.
4.4. Shape Relaxation of a Square Droplet.
The square-droplet benchmark examines surface-tension-driven relaxation. High-curvature corners smooth rapidly, the droplet approaches a circular equilibrium, and enclosed phase mass is preserved.
- The computation uses a 256 × 256 uniform triangular mesh, zero initial velocity and pressure, no-slip walls, and ∆t = 0.005.
- The square corners rapidly smooth, rounded-square configurations persist for approximately 0.035 ≤t ≤0.14, and an approximately circular equilibrium is approached by t = 1.0.
- The computation remains stable throughout the simulated interval and preserves enclosed phase mass to numerical precision.
4.5. Two-Phase Lid-Driven Cavity Flow.
The lid-driven cavity test probes strong shear and substantial interfacial deformation. The moving lid creates a clockwise vortex that stretches, rolls, folds, and transports the interface while the computation remains stable.
- A horizontal diffuse interface is initialized at y = 0.5, with no-slip side and bottom walls and a smoothly vanishing driven velocity at the upper corners.
- The simulation uses ν = 0.002, η = 0.01, λ = 2 × 10−6, M = 0.005, ϵ = 10−5, and P2 −P1 −P2 −P2 elements.
- The moving lid generates a clockwise primary vortex that deforms the initially horizontal interface into a pronounced interfacial wave by t = 5.
- From approximately 7.5 ≤t ≤10, the interface is stretched, rolled, and transported upward before forming a thin spiraling structure through substantial folding.
- The computation remains stable over the full simulated interval and resolves the strongly deformed diffuse interface without visible spurious oscillations.
4.6. Rayleigh–Taylor Instability.
The Rayleigh–Taylor benchmark tests density-difference-driven interfacial dynamics across two viscosity regimes. Higher viscosity yields smoother evolution, whereas lower viscosity produces secondary roll-up structures, with stability maintained in both cases.
- The benchmark places a denser fluid above a lighter fluid in a gravitational field and uses viscosity values ν = 0.01 and ν = 0.001.
- For ν = 0.01, viscous damping suppresses small-scale roll-up, producing a relatively smooth symmetric spike and broad bubble-like structures.
- For ν = 0.001, stronger inertial effects intensify side shear layers and produce secondary roll-up structures by approximately t = 0.72.
- The low-viscosity structures continue deforming and entraining surrounding fluid as the dense spike descends.
- The discretization remains numerically stable for both viscosity values and captures the qualitative transition from strongly damped to more intricate interfacial flow.