Coupled Cahn-Hilliard Navier-Stokes Model

• The coupled Cahn-Hilliard–Navier-Stokes model is a framework that captures interfacial dynamics, phase transitions, and energy dissipation in multiphase flows.
• It integrates phase-field thermodynamics with incompressible fluid dynamics, modeling applications like droplet coalescence and bubble rise.
• Advanced numerical techniques, such as projection methods and operator-splitting, ensure stability and second-order temporal accuracy in simulations.

The coupled Cahn-Hilliard–Navier-Stokes (CHNS, or NSCH) model is a foundational framework for simulating incompressible multiphase flows with diffuse interfaces, capturing the interfacial dynamics, topological transitions, and two-way phase–hydrodynamics coupling central to multiphysics phenomena such as droplet coalescence, bubble rise, and hydrodynamic instabilities. Below, advanced aspects of the model and its implementation are elaborated, noting both canonical theory and recent methodological advances.

1. Governing Equations and Phase-Field Thermodynamics

The CHNS model describes the time evolution of a conserved order parameter $\phi(\mathbf{x},t)\in[-1,1]$, which demarcates two fluid phases ($\phi=-1$ and $\phi=+1$), coupled to the velocity $\mathbf{u}(\mathbf{x},t)$ and pressure $p(\mathbf{x},t)$ of an incompressible fluid. The total free energy over the domain $\Omega$ is

$$\mathcal{E}[\phi]=\int_{\Omega}F(\phi)+\frac{\epsilon^2}{2}|\nabla\phi|^2\,d\mathbf{x}\,,\quad F(\phi)=\frac{1}{4}(\phi^2-1)^2\,,$$

with $\epsilon>0$ the diffuse-interface thickness parameter. The chemical potential is the variational derivative

$$\mu = \frac{\delta\mathcal{E}}{\delta\phi} = -\epsilon^2\nabla^2\phi + \phi(\phi^2-1)\,.$$

Phase-field evolution obeys the advective Cahn–Hilliard equation,

$$\frac{\partial\phi}{\partial t}+\mathbf{u}\cdot\nabla\phi = \nabla\cdot[M(\phi)\nabla\mu]\,,$$

where $M(\phi)\geq0$ is the mobility, which may be constant or degenerate.

Mass and momentum balance for the incompressible fluid are governed by

$$\rho(\phi)\left[\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] = -\nabla p + \nabla\cdot\left[\eta(\phi)\left(\nabla\mathbf{u}+(\nabla\mathbf{u})^T\right)\right] + \mathbf{f}_\text{body} -\mu\nabla\phi\,,\qquad \nabla\cdot\mathbf{u}=0.$$

Fluid properties are typically interpolated linearly: $$\rho(\phi) = \frac{1+\phi}{2}\rho_1+\frac{1-\phi}{2}\rho_2,\quad \eta(\phi) = \frac{1+\phi}{2}\eta_1+\frac{1-\phi}{2}\eta_2\;.$$ The singular capillary force $-\mu\nabla\phi$ arises naturally from the phase-field free-energy, enforcing surface tension; $\mathbf{f}_\text{body}$ includes external accelerations such as gravity.

2. Mathematical Structure and Coupling Mechanisms

The CHNS model is deeply coupled:

• The fluid flow advects the order parameter ($\mathbf{u}\cdot\nabla\phi$), transporting the interface.
• The phase field determines spatially varying $\rho(\phi)$ and $\eta(\phi)$, which in turn modulate fluid inertia and viscous stresses.
• The chemical potential gradient $-\mu\nabla\phi$ provides the Korteweg (capillary) force in the Navier–Stokes equations.
• For variable-density/variable-viscosity systems, these nonlinearities require careful numerical interpolation at grid faces and momentum-coupling steps.

The model exhibits a formal energy-dissipation law: $$\frac{d}{dt}\left(\frac{1}{2}\int_\Omega\rho|\mathbf{u}|^2\,d\mathbf{x}+\mathcal{E}[\phi]\right) = -\int_\Omega[\eta|\nabla\mathbf{u}|^2+M|\nabla\mu|^2]\,d\mathbf{x}\leq0,$$ guaranteeing thermodynamic consistency at the PDE level (Manna et al., 27 Aug 2025, Zhao, 2021).

3. Numerical Schemes: Projection Methods and Decoupling

A widely adopted simulation strategy is the pressure-velocity projection splitting (Manna et al., 27 Aug 2025):

1. Velocity Prediction: Advance $\mathbf{u}^n\rightarrow\tilde{\mathbf{u}}^{n+1}$ without pressure, treating advection, viscosity, external force, and capillarity with explicit Euler or higher-order (semi-)implicit schemes.
2. Pressure Correction: Solve a variable-coefficient elliptic (Poisson) problem to determine $p^{n+1}$ ensuring $\nabla\cdot\mathbf{u}^{n+1}=0$.
3. Velocity Update: Project the predictor velocity onto the divergence-free space.

Spatial discretization commonly uses staggered (MAC) finite-difference grids for mass, phase, and momentum conservation. Scalar quantities $(\phi,\mu,p)$ are cell-centered, velocities are face-centered. Convective terms are discretized in a skew-symmetric or central-difference form to control numerical dissipation and conserve kinetic energy. Laplacian terms leverage standard finite-difference stencils, and variable properties are interpolated arithmetically at cell faces to ensure stability and consistency.

Explicit Euler is widely used for temporal integration of both momentum and phase-field equations; the time step is constrained by both fluid CFL and fourth-order Cahn–Hilliard diffusion: $$\Delta t \leq C_N h/|\mathbf{u}|_{\max},\qquad \Delta t \leq C_\text{CH} h^4/(M\epsilon^2).$$ Recent advances provide unconditionally energy-stable, second-order-accurate decoupled operator-splitting schemes utilizing Strang–Marchuk splitting, enabling efficient wetted interface dynamics, second-order temporal accuracy, and rigorous discrete energy dissipation (Zhao, 2021).

4. Variable Material Properties and Stabilization

High density or viscosity contrasts present stiff numerical challenges due to sharp interfacial gradients. These are addressed by:

• Averaging at Faces: For pressure-correction and viscous terms, densities and viscosities at faces are constructed via arithmetic averaging.
• Smoothing: Slightly smoothing or clipping the phase-field $\phi$ prior to property interpolation mitigates convergence degradation in the pressure Poisson solve.
• Preconditioning: Multigrid preconditioners tailored for variable-coefficient problems enhance linear solver performance.
• Interface Thickness: Numerically thickened interfaces (large $\epsilon$ relative to physical) help buffer extreme property jumps, alleviating fourth-order time-step constraints.

Unconditional energy-stability, mass conservation, and global free-energy monotonicity can be numerically demonstrated for sufficiently small time steps (Manna et al., 27 Aug 2025, Zhao, 2021).

5. Physical Benchmarking and Application Domains

Projection-based NSCH frameworks have been validated on canonical multiphase benchmarks:

• Rising Bubble: Simulations of gas bubbles in liquid channels (e.g., $R=0.25$, $\rho_{liq}=1000$, $\rho_{gas}=100$, $\eta_{liq}=10$, $\eta_{gas}=1$, $\sigma=24.5$) reproduce dome curvature and terminal velocities in microchannels to within $2\%$ of experimental values (Bhaga & Weber 1981). Wake vorticity and vortex-pair structure agree with experimental and numerical references.
• Hydrodynamic Instabilities: Plateau–Rayleigh (Rayleigh–Taylor) instabilities are accurately reproduced. Growth rates, instability wavelengths, and nonlinear interface evolution (spike/bubble shapes, vortex roll-up) are consistent with linear stability theory and high-resolution DNS reference simulations (Manna et al., 27 Aug 2025).
• Energy and Mass Tracking: Staggered-grid projection schemes maintain mass conservation and can demonstrate near-monotonic decay of global energy, provided central-difference fluxes and adequately small time steps are used.

These techniques apply broadly to multiphysics scenarios including boiling, droplet manipulation in microfluidics, additive manufacturing, and multiphase separation processes.

6. Extensions, Open Problems, and Analytical Developments

Several directions have been established or are active in the literature:

• Thermal Coupling: Extensions to non-isothermal models incorporate additional equations for (inverse) temperature and internal energy, maintaining thermodynamic consistency through appropriate coupling of temperature, chemical potential, and velocity (Berti et al., 2011, Brunk et al., 2024).
• Nonlocality and Boundary Coupling: Models with nonlocal free energy kernels and/or bulk-surface coupling accommodate more complex phase interactions and dynamic boundary effects (Frigeri et al., 2013, Stange, 10 Nov 2025).
• Thermodynamic Consistency: All models preserve a global energy dissipation structure and, under no-flux and no-slip boundary conditions, satisfy discrete (or continuous) analogs of the Clausius–Duhem inequality.
• Regularity Theory: Recent deterministic PDE analysis has established global well-posedness, strong solution uniqueness, and attractor structures in two spatial dimensions. For example, well-posedness with non-degenerate mobilities and singular potentials, and eventual regularization of weak solutions to strong, is now available for nonlocal CHNS systems (Frigeri et al., 2013, Stange, 10 Nov 2025).
• Optimal Control and Reduced-Order Modeling: Formulations for boundary or distributed control, as well as model-order reduction and optimality systems for coupled CHNS, are advancing toward computationally efficient simulation and control (Signori et al., 26 Sep 2025, Bag et al., 2023, Gräßle et al., 2019).

Open technical challenges persist for:

• Global regularity and long-time dynamics in three-dimensional, nonlocal, and variable-property regimes.
• Handling extreme property ratios, singular potentials, and degenerate mobility cases with guaranteed unconditional stability.
• Coupling with additional physics (contact lines, electrowetting, reactive flows, etc.) in both modeling and numerics.

7. Summary Table: Representative Model Components

Aspect	Mathematical Structure	Reference Equation(s)
Free Energy Functional	$$\mathcal{E}[\phi]=\int F(\phi)+\frac{\epsilon^2}{2}|\nabla\phi|^2$$	(1) (Manna et al., 27 Aug 2025)
Chemical Potential	$\mu=-\epsilon^2\nabla^2\phi+\phi(\phi^2-1)$	(2) (Manna et al., 27 Aug 2025)
Phase-Field Evolution	$$\partial_t\phi+\mathbf{u}\cdot\nabla\phi=\nabla\cdot[M(\phi)\nabla\mu]$$	(3) (Manna et al., 27 Aug 2025)
Navier–Stokes–Cahn–Hilliard	$$\partial_t(\rho\mathbf{u})+\cdots=-\nabla p+\nabla\cdot(\eta D(\mathbf{u}))-\mu\nabla\phi$$	(5) (Manna et al., 27 Aug 2025)
Property Interpolation	$\rho(\phi),\ \eta(\phi)$ via linear maps	(4) (Manna et al., 27 Aug 2025)
Pressure Correction	$$\nabla\cdot\left(\frac{1}{\rho^n}\nabla p^{n+1}\right) = \frac{1}{\Delta t}\nabla\cdot\tilde{\mathbf{u}}^{n+1}$$	(9) (Manna et al., 27 Aug 2025)
Energy Dissipation Law	$$\frac{d}{dt}\left(E_{\text{tot}}\right) = -\int(\cdots)\leq0$$	(Zhao, 2021, Manna et al., 27 Aug 2025)

For the detailed validation, discretization, and coupling mechanisms, see (Manna et al., 27 Aug 2025) and (Zhao, 2021). These provide the canonical, reproducible foundation for multiphase flow simulation with diffuse interfaces, and underpin numerous recent analytical and computational advances in phase-field hydrodynamics.