Phase-Field Modeling of Two-Phase Flows: A Projection-Based Cahn-Hilliard-Navier-Stokes Framework

TL;DR

This paper develops a projection-based phase-field model combining Cahn-Hilliard and Navier-Stokes equations for simulating complex two-phase flows with diffuse interfaces, addressing numerical challenges and validating against benchmark problems.

Contribution

It introduces a decoupled pressure-projection numerical scheme for the coupled CH-NS system with variable properties, enabling stable and accurate multiphase flow simulations.

Findings

Accurately predicts bubble rise velocity and interface shape.

Successfully captures hydrodynamic instabilities like Plateau-Taylor.

Framework is extendable to heat transfer and phase change phenomena.

Abstract

The coupled Cahn-Hilliard and Navier-Stokes (CH-NS) equations provide a powerful framework for modeling multiphase flows with diffuse interfaces, enabling simulations of droplet breakup, bubble dynamics, and hydrodynamic instabilities. These capabilities are vital in boiling heat transfer, microfluidics, coating, additive manufacturing, and oil-water separation, where resolving fluid-fluid interactions is essential. Numerically, the CH-NS system is challenging: the Cahn-Hilliard equation involves higher-order derivatives and nonlinearities, and coupling with Navier-Stokes introduces strong two-way interactions. The velocity field advects the phase field, while the evolving interface alters density and viscosity, feeding back into the flow. Variable-density and variable-viscosity systems further increase complexity, requiring accurate treatment of property contrasts without losing stability or mass conservation. To address this, we employ a decoupled pressure-projection method with finite differences on staggered grids and explicit Euler time stepping. Our formulation extends the CH-NS system to homogeneous and variable-property fluids with consistent hydrodynamic-phase-field coupling. Validation against canonical benchmarks-including bubble rise and Plateau-Taylor instability shows excellent agreement in rise velocity, interface shape, and instability wavelength. This framework establishes a reproducible foundation for multiphysics extensions such as heat transfer, phase change, and electrohydrodynamics in boiling, droplet manipulation, and electronics cooling

Equations

$$\varphi =\frac{m_1-m_2}{m_1+m_2},c=\frac{m_1}{m_1+m_2},$$

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

$$\mathcal{E}(\varphi )={\int }_{\Omega }\left[F(\varphi )+\frac{{ϵ}^2}{2}|\nabla \varphi {|}^2\right]d\mathbf{x},$$

$$\mu =\frac{\delta \mathcal{E}}{\delta \varphi }.$$

$$\mathcal{E}(\varphi )=\int \left[\frac{{ϵ}^2}{2}|\nabla \varphi {|}^2+\frac{1}{4}({\varphi }^2-1{)}^2\right]dV.$$

$$\delta \mathcal{E}=\mathcal{E}(\varphi +\delta \varphi )-\mathcal{E}(\varphi ).$$

$$\mathcal{E}(\phi+\delta\phi)=\int\Bigg{[}\frac{\epsilon^{2}}{2}|\nabla(\phi+\delta\phi)|^{2}\\ +\frac{1}{4}\left((\phi+\delta\phi)^{2}-1\right)^{2}\Bigg{]}\,dV.$$

$$|\nabla (\varphi +\delta \varphi ){|}^2$$

$$=|\nabla \varphi {|}^2+2\nabla \varphi \cdot \nabla \delta \varphi +|\nabla \delta \varphi {|}^2.$$

$$|\nabla (\varphi +\delta \varphi ){|}^2\approx |\nabla \varphi {|}^2+2\nabla \varphi \cdot \nabla \delta \varphi .$$

$$((\varphi +\delta \varphi {)}^2-1{)}^2={\left[{\varphi }^2+2\varphi \delta \varphi +(\delta \varphi {)}^2-1\right]}^2.$$

$$((\varphi +\delta \varphi {)}^2-1{)}^2\approx ({\varphi }^2-1{)}^2+4\varphi ({\varphi }^2-1)\delta \varphi .$$

$$\mathcal{E}(\phi+\delta\phi)\approx\int\Bigg{[}\frac{\epsilon^{2}}{2}\left(|\nabla\phi|^{2}+2\nabla\phi\cdot\nabla\delta\phi\right)\\ +\frac{1}{4}\left((\phi^{2}-1)^{2}+4\phi(\phi^{2}-1)\,\delta\phi\right)\Bigg{]}\,dV.$$

$$\delta \mathcal{E}\approx \int \left[{ϵ}^2\nabla \varphi \cdot \nabla \delta \varphi +\varphi ({\varphi }^2-1)\delta \varphi \right]dV.$$

$$\int {ϵ}^2\nabla \varphi \cdot \nabla \delta \varphi dV=-\int {ϵ}^2({\nabla }^2\varphi )\delta \varphi dV.$$

$${\nabla }^2\varphi =\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}+\frac{{\partial }^2\varphi }{\partial z^2}.$$

$$\delta \mathcal{E}=\int \left[-{ϵ}^2{\nabla }^2\varphi +\varphi ({\varphi }^2-1)\right]\delta \varphi dV.$$

$$\frac{\delta \mathcal{E}}{\delta \varphi }=-{ϵ}^2{\nabla }^2\varphi +\varphi ({\varphi }^2-1).$$

$$\mu =-{ϵ}^2{\nabla }^2\varphi +\varphi ({\varphi }^2-1).$$

$$\rho (\varphi )$$

$$\eta (\varphi )$$

$$\begin{gathered} \rho (\varphi )\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}\right]=-\nabla p+\nabla \cdot \left[\eta (\varphi )D(\mathbf{u})\right] \\ +\mathbf{F}, \end{gathered}$$

$$\begin{gathered} \rho (\varphi )\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}\right]=-\nabla p+\nabla \cdot \left[\eta (\varphi )D(\mathbf{u})\right] \\ +{\mathbf{SF}}_{\text{sing}}+\mathbf{F}, \end{gathered}$$

$${\mathbf{SF}}_{\text{sing}}=\tilde{\sigma }{\epsilon }^{-1}\mu \nabla \varphi ,$$

$$\rho (\varphi )\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}\right]$$

$$+\nabla \cdot \left[\eta (\varphi )D(\mathbf{u})\right]$$

$$+{\mathbf{SF}}_{\text{sing}}+\mathbf{F},$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{u}\cdot \nabla \varphi =\nabla \cdot \mathcal{J},$$

$$\mathcal{J}=-M(\varphi )\nabla \mu .$$

$$\mu =-{ϵ}^2{\nabla }^2\varphi +\varphi ({\varphi }^2-1),$$

$${\mathbf{SF}}_{\text{sing}}=\tilde{\sigma }{ϵ}^{-1}\mu \nabla \varphi ,$$

Full text

††preprint: APS/123-QED

1. I Introduction

2. II Methodology

3. II.1 Phase-Field Representation

4. II.2 Definition of the Phase-Field Variable

5. II.3 Cahn–Hilliard Equation with Advection

6. II.4 Chemical Potential in the Phase-Field Model

7. II.5 Hydrodynamic Model: Navier–Stokes with Surface Tension

8. II.6 Governing Equations in the Phase-Field Framework

9. II.7 Property Variation Across the Interface

10. II.8 Discretization Strategy

11. II.9 Time Discretization of the Momentum Equation

12. II.10 Time Discretization for the Cahn–Hilliard Equation

13. II.11 Space Discretization

14. i Scalar Variables on Cell Centers

15. ii Velocity Variables on a Staggered Grid

16. II.12 Constant Density and Viscosity

17. i Space Discretization of the Momentum Equation

18. ii Space Discretization of the Pressure Computation Equation

19. iii Space Discretization of the Chemical Potential

20. iv Space Discretization of the Phase Field

21. II.13 Spatial Discretization: Non-Matched Properties

22. II.14 Pressure Correction for Variable Density and Viscosity

23. i Expansion of Term (i)

24. ii Expansion of Term (ii)

25. II.15 Inhomogeneous Density Correction in Divergence-Free Velocity Evaluation

26. III Results & Discussion

27. III.1 Rising Bubble with Density and Viscosity Contrast

28. III.2 Rayleigh–Taylor Instability

29. IV Conclusion

30. V Appendix

31. V.1 Expression for Free Energy Functional

32. i Derivation of $\kappa_{ij}^{(1)}$

33. ii Derivation of $\kappa_{ij}^{(2)}$

34. V.2 Advective Cahn–Hilliard Formulation

35. V.3 Derivation of Linear Interpolation Formulae

36. i Viscosity interpolation

37. ii Density interpolation

I Introduction

Multiphase flows involving two or more immiscible fluids are central to a wide range of natural processes and engineering applications, including boiling heat transfer Li et al. (2025), microfluidics Raynaldo et al. (2024), inkjet printing Wang et al. (2025), coating flows Li et al. (2024), additive manufacturing Liu et al. (2025), and oil–water separation Wu et al. (2001). The dynamics of these systems are governed by the interplay between fluid motion, interfacial tension, and property variations across the interface. Accurate prediction of interfacial evolution is therefore critical for understanding and optimizing such processes.

A variety of numerical strategies have been developed to model multiphase flows. In sharp-interface methods such as the level set (LS) Gibou et al. (2018) and volume of fluid (VOF) Gu et al. (2019), the interface is assumed to have zero thickness. In contrast, diffuse-interface (phase-field) methods represent the interface as a smooth transition region of finite thickness. This concept, dating back to Rayleigh Rayleigh (1892) and van der Waals Van der Waals (1979), has since become a powerful and widely used tool for studying interfacial dynamics Steinbach (2009); Kim (2012); Lamorgese et al. (2011). A key advantage of the diffuse-interface approach is that the governing equations can be derived from an energy-based variational framework consistent with hydrodynamic principles. However, ensuring that numerical schemes respect the discrete energy law is critical: violations lead to artificial dissipation or dispersion near the interface and can generate large numerical errors Wang and Wise (2011). Thus, despite its success, the development of accurate, energy-stable, and general-purpose multiphase flow solvers remains an active research area.

In the diffuse-interface framework, the Cahn–Hilliard equation Cahn (1961) governs the evolution of the order parameter, naturally capturing interface motion, deformation, and topological transitions. When coupled with the Navier–Stokes equations, it yields the Cahn–Hilliard–Navier–Stokes (CH–NS) system, a tightly coupled nonlinear model in which the velocity field advects the phase field, while the evolving interface modifies density and viscosity distributions that feed back into the flow. This bidirectional coupling introduces significant challenges, particularly in variable-density and variable-viscosity systems.

Over the years, a variety of numerical methods have been proposed to address these challenges. Convex splitting schemes Han and Wang (2015) are unconditionally stable and widely used for the nonlinear chemical potential, though they can be computationally expensive. Shen and Yang Shen and Yang (2010) developed energy-stable schemes for both constant- and variable-density CH–NS flows, while Liu et al. Liu et al. (2022) proposed a fully decoupled and unconditionally stable scheme for the Allen–Cahn–Navier–Stokes system. The scalar auxiliary variable (SAV) approach Shen et al. (2018) introduces an auxiliary variable to reformulate the free energy, simplifying the discretization, but at the cost of modifying the original energy structure. Similarly, invariant energy quadratization (IEQ) Yang and Yu (2018) offers second-order accuracy for complex nonlinearities, though its link to the original energy dissipation law holds only under certain conditions Zhang et al. (2021). Chen and Zhao Chen and Zhao (2020) introduced linear, second-order approaches with modified leap-frog time marching. Adaptive methods, such as those of Chen and Shen Chen and Shen (2016), refine both space and time dynamically to better capture localized interfacial dynamics while maintaining stability.

For the Navier–Stokes component, a wide family of splitting and projection methods has been developed. Guermond and Salgado Guermond and Salgado (2009) proposed a penalty-based splitting method requiring only a constant-coefficient Poisson solve per time step, while Guermond and Minev Guermond and Minev (2019) introduced adaptive third-order time stepping for incompressible flows. Guermond et al. Guermond et al. (2006) summarized major classes of projection schemes, including non-incremental pressure correction, standard incremental pressure correction, velocity-correction, and consistent splitting methods. Fractional-step approaches have also been extended to moving contact line problems with variable density and viscosity Gao and Wang (2014), where convex splitting of the Cahn–Hilliard equation ensures stability under mild assumptions.

The CH–NS framework has been employed to study a wide range of canonical multiphase flow problems. Shen and Yang Shen and Yang (2010) examined the dynamics of single rising bubbles in heavier media under both constant- and variable-density formulations, while Shen et al. Shin et al. (2020) investigated bubble deformation in shear flows using high-order polynomial free energies. Yang and Kim Yang and Kim (2020) analyzed bubble dynamics with variable mobility, highlighting the role of transport properties in interfacial evolution. In addition, adaptive phase-field solvers Chen and Shen (2016) and spectral-element formulations Xiao et al. (2022) have been developed to improve computational efficiency and accuracy. Benchmark studies of rising bubbles have also been carried out across multiple numerical approaches Hysing et al. (2009), including finite-element implementations Aland and Voigt (2012), providing valuable reference data for validation and comparison.

In this work, we consolidate and present a self-contained derivation of the CH–NS system for both homogeneous and variable density and viscosity two-fluid systems, together with a comprehensive description of our discretization strategies. Our implementation employs a finite-difference framework with staggered grids for spatial discretization and explicit Euler integration in time. A decoupled pressure-projection method is used to stabilize the hydrodynamic coupling. Emphasis is placed on mass conservation, numerical stability, and computational efficiency. To validate the solver, we present results for two benchmark problems: (i) the rising bubble, and (ii) the Plateau–Rayleigh instability. Quantitative comparisons with reference data demonstrate excellent agreement in rise velocity, interface deformation, and instability wavelength. This reproducible framework provides a solid foundation for future extensions to multiphysics problems such as heat transfer, phase change, and electrohydrodynamic effects Szczukiewicz et al. (2014); Rothman and Zaleski (1994); Azizian et al. (2019).

The remainder of this article is organized as follows. Sections II–II.6 introduce the governing equations of the CH–NS system for both constant- and variable-density/viscosity formulations, together with the relevant dimensionless parameters. Section II.8 describes the spatial and temporal discretization strategies, including the treatment of coupling terms and property variations. Section III presents the benchmark problems, numerical setup, and corresponding results, with comparisons to reference data. Finally, Section IV summarizes the main findings and outlines potential directions for future work.

II Methodology

II.1 Phase-Field Representation

We consider a two-dimensional microchannel domain $\Omega$ containing two immiscible, incompressible, viscous fluids with constant densities $\rho_{1}$ and $\rho_{2}$ and dynamic viscosities $\eta_{1}$ and $\eta_{2}$, respectively. The spatial distribution of the two fluids is described by a phase-field variable $\phi(\mathbf{x},t)$, where $\mathbf{x}$ denotes the spatial coordinate and $t$ the time. The order parameter $\phi$ takes the bulk value $\phi=-1$ in Fluid 1 and $\phi=+1$ in Fluid 2, with a diffuse transition layer of finite thickness $\epsilon$ representing the fluid–fluid interface.

Figure 1 illustrates this diffuse-interface description. Panel (a) shows the microchannel geometry, where a gas bubble (Fluid 1) is surrounded by a liquid (Fluid 2). Panel (b) provides an enlarged view of the bubble–liquid system, highlighting the subdomains $\Omega_{1}$ and $\Omega_{2}$ occupied by Fluid 1 and Fluid 2, respectively. The interface $\Gamma$ separates the two fluids, with the unit normal vector $\mathbf{n}$ defined as pointing from the liquid (Fluid 2) toward the gas bubble (Fluid 1). Panel (c) depicts the smooth variation of the order parameter $\phi$ across the interfacial region (blue shaded layer), which enables the natural treatment of interfacial dynamics and material property variations without explicit interface tracking.

II.2 Definition of the Phase-Field Variable

In the phase-field method, the interface emerges naturally from the continuous variation of an order parameter between two bulk values. For a binary incompressible mixture of a liquid (Fluid 1, also denoted A) and a gas (Fluid 2, also denoted B), each with distinct densities $\rho_{1},\rho_{2}$ and viscosities $\eta_{1},\eta_{2}$, the phase-field variable Liu and Shen (2003); Kim (2005) can be defined in terms of the local masses $m_{1}$ and $m_{2}$ of the two fluids as

[TABLE]

where $c$ denotes the volume fraction of Fluid 1 (liquid A). By construction, $\phi\in[-1,1]$ and $c\in[0,1]$. In this work, $\phi$ serves as the primary order parameter, while $c$ is used when convenient for interpretation.

II.3 Cahn–Hilliard Equation with Advection

The temporal evolution of the phase field is governed by the Cahn–Hilliard (CH) equation Cahn (1961), which models the transport of a conserved scalar field subject to advection by the bulk velocity $\mathbf{u}$ and non-Fickian diffusion driven by gradients in chemical potential $\mu$:

[TABLE]

where $M(\phi)\geq 0$ is the mobility. The choice of $M(\phi)$—whether taken as a constant or made phase-dependent—determines whether bulk or interfacial diffusion dominates.

The chemical potential $\mu$ is obtained as the variational derivative of the total free energy functional (see Appendix V.1 for the complete derivation).

[TABLE]

where $F(\phi)=\frac{1}{4}(\phi^{2}-1)^{2}$ is the double-well Helmholtz free energy density and $\epsilon$ is related to the interfacial thickness. The first term in $\mathcal{E}$ represents the bulk free energy, while the second term accounts for interfacial energy by penalizing sharp gradients in $\phi$.

For the concentration variable $c$, an equivalent form $F(c)=\frac{1}{4}c^{2}(c-1)^{2}$ may be used. Both formulations capture phase separation and interfacial tension effects within the same energetic framework. The detailed derivation of Eq. (2) are shown in Appendix V.2.

II.4 Chemical Potential in the Phase-Field Model

Given the free-energy functional $\mathcal{E}(\phi)$ (Section V.1), the chemical potential $\mu$ follows from its variational derivative with respect to the order parameter $\phi$.

[TABLE]

The free energy functional is given by

[TABLE]

Let $\delta\phi$ be a small perturbation in $\phi$. The change in $\mathcal{E}(\phi)$ is

[TABLE]

Substituting $\phi+\delta\phi$ into Eq. (5):

[TABLE]

For the gradient term:

[TABLE]

Neglecting the higher-order term $|\nabla\delta\phi|^{2}$:

[TABLE]

For the polynomial term:

[TABLE]

Neglecting $(\delta\phi)^{2}$:

[TABLE]

Substituting these into $\mathcal{E}(\phi+\delta\phi)$:

[TABLE]

Subtracting $\mathcal{E}(\phi)$ from both sides:

[TABLE]

Applying integration by parts to the gradient term and assuming $\delta\phi=0$ on $\partial\Omega$:

[TABLE]

Here, $\nabla^{2}\phi$ is the Laplacian:

[TABLE]

Combining results:

[TABLE]

From the definition of the functional derivative:

[TABLE]

Thus, the chemical potential is:

[TABLE]

II.5 Hydrodynamic Model: Navier–Stokes with Surface Tension

Following the phase-field formulation and chemical potential introduced in Section II.4, we couple the advective Cahn–Hilliard equation to the incompressible Navier–Stokes equations. Such coupling is a standard feature of diffuse-interface methods Kim (2012); Anderson et al. (1998); He and Kasagi (2008), allowing the accurate treatment of interfacial dynamics in multiphase flows.

Let $\rho_{1}$ and $\rho_{2}$ denote the densities of the two fluids (bubble and liquid). In the phase-field framework, density and viscosity are expressed as linear functionsLiu and Shen (2003) of the order parameter $\phi$:

[TABLE]

With these definitions, the governing equations for an unsteady, viscous, incompressible, and immiscible two-fluid system take the form:

[TABLE]

where $\partial_{t}\mathbf{u}$ is the local acceleration, $(\mathbf{u}\cdot\nabla)\mathbf{u}$ is the convective acceleration, $-\nabla p$ is the pressure gradient force, $\nabla\cdot[\eta(\phi)D(\mathbf{u})]$ represents the viscous force per unit volume, and $\mathbf{F}$ accounts for external body forces (e.g., gravity).

Here $D(\mathbf{u})=\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}$ is the rate-of-deformation tensor, i.e., the symmetric part of the velocity gradient.

At the interface $\Gamma$ between the two immiscible fluids (see Fig. 1), discontinuities in pressure and viscous stresses arise due to surface tension and property contrasts Kim (2012). Surface tension balances the jump in normal stresses, ensuring mechanical equilibrium Kim (2012); Uzgoren et al. (2009). Incorporating this effect, the Navier–Stokes equations become:

[TABLE]

where $\mathbf{SF}_{\text{sing}}$ denotes the singular surface tension force.

In our phase-field formulation, we adopt the surface tension model Kim (2005):

[TABLE]

where $\tilde{\sigma}$ is a scaled surface tension parameter (proportional to the physical surface tension $\sigma$), $\mu$ is the chemical potential from Eq. (18), and $\varepsilon$ is the characteristic length scale over which the interface is diffused.

II.6 Governing Equations in the Phase-Field Framework

Extending the chemical potential formulation from the previous section, we couple the phase-field model to the Navier–Stokes equationsJacqmin (1999) for two immiscible, viscous fluids. The resulting system enforces momentum conservation, incompressibility, and phase-field evolution, enabling a consistent description of bulk flow and interfacial physics.

Momentum conservation.

The momentum equation is derived from the incompressible Navier–Stokes equations. It accounts for density and viscosity variations through the phase field $\phi$, and includes contributions from surface tension and external body forces:

[TABLE]

where $\rho(\phi)$ is the local density, $\mathbf{u}$ is the velocity field, $p$ is the pressure, $\eta(\phi)$ is the viscosity, $D(\mathbf{u})$ is the rate-of-deformation tensor, $\mathbf{SF}_{\text{sing}}$ is the singular surface tension force, and $\mathbf{F}$ represents external body forces such as gravity.

Incompressibility condition.

For incompressible flows, the velocity field must be divergence-free:

[TABLE]

ensuring conservation of mass in the fluid domain.

Phase-field evolution.

The interface motion is governed by the advective Cahn–Hilliard equation:

[TABLE]

where $\mathcal{J}$ is the phase-field flux, given by

[TABLE]

Chemical potential.

Using the variational derivative of the free energy functional $\mathcal{E}(\phi)$ (see Section II.4), the chemical potential is

[TABLE]

where $\epsilon$ is the interface thickness parameter.

Surface tension force.

The singular surface tension force is modeled as

[TABLE]

where $\tilde{\sigma}$ is a scaled surface tension parameter proportional to the physical surface tension $\sigma$.

Equations (24)–(29) form the coupled Navier–Stokes–Cahn–Hilliard system used to simulate bubble dynamics in the present study.

II.7 Property Variation Across the Interface

As described in Section II.6, the order parameter $\phi$ distinguishes the two phases:

[TABLE]

The densities $(\rho_{1},\rho_{2})$ and viscosities $(\eta_{1},\eta_{2})$ of these phases vary smoothly across the interface, and are obtained via linear interpolation (full derivation in Appendix V.3).

Truncated phase field.

Solutions of the advective Cahn–Hilliard equation may yield $\phi$ values outside the physical range $[-1,1]$. This occurs because the equation does not strictly satisfy the maximum principle Shen and Yang (2015). To prevent unphysical overshoots, we define a truncated phase field:

[TABLE]

where $\mathrm{sign}(\phi)$ returns $+1$ or $-1$.

Interpolated properties.

Replacing $\phi$ by $\hat{\phi}$ in the interpolation formulas gives:

[TABLE]

These expressions ensure smooth property variation across the diffuse interface while keeping $\phi$ within physical bounds.

II.8 Discretization Strategy

The coupled Cahn–Hilliard–Navier–Stokes (CH–NS) system governs both interfacial dynamics and bulk fluid motion in two-phase flow. Solving the fully coupled equations in a monolithic fashion is computationally demanding and can suffer from stability issues, particularly when large density or viscosity contrasts are present Guermond et al. (2006). To improve numerical stability, reduce memory requirements, and efficiently leverage existing solvers, we employ a staggered time-stepping (operator-splitting) scheme.

At each time step $t^{n}\to t^{n+1}$, the solution procedure is divided into two stages:

1. Momentum and pressure update. The Navier–Stokes equations (24) are solved using the phase field $\phi^{n}$ and chemical potential $\mu^{n}$ from the previous step, yielding the updated velocity $\mathbf{u}^{n+1}$ and pressure $p^{n+1}$.

2. Phase-field update. With $\mathbf{u}^{n+1}$, the Cahn–Hilliard equation (26) is advanced to compute $\phi^{n+1}$ and $\mu^{n+1}$, ensuring accurate interface transport and diffusion while preserving incompressibility.

This fractional-step strategy decouples the CH and NS systems, allowing each subproblem to be treated with methods best suited to its characteristics. In practice, this approach maintains second-order temporal accuracy, good mass conservation, and robustness across a range of multiphase flow conditions Shen and Yang (2010).

II.9 Time Discretization of the Momentum Equation

In the first stage of the staggered update, we advance the velocity and pressure fields. The momentum conservation equation is

[TABLE]

where $\rho(\phi)$ and $\eta(\phi)$ are the phase-dependent density and viscosity, $\mathbf{u}$ is the velocity field, $p$ is the pressure, $D(\mathbf{u})=\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}$ is the rate-of-deformation tensor, $\mathbf{F}$ denotes external body forces (e.g., gravity), $\tilde{\sigma}$ is the scaled surface tension coefficient, $\varepsilon$ is the interface thickness parameter, $\mu$ is the chemical potential, and $\phi$ is the phase-field variable.

We approximate the time derivative as

[TABLE]

where $\mathbf{u}^{n}$ is the velocity at the current time step, $\tilde{\mathbf{u}}^{n+1}$ is an intermediate velocity field, and $\Delta t$ is the time step size.

Substituting into the momentum equation gives

[TABLE]

Here, Term I: pressure gradient contribution.

Term II: viscous stress, external body forces, and surface tension effects.

We compute the intermediate velocity in four steps:

Step 1: Intermediate velocity without pressure Omitting Term I:

[TABLE]

Rewriting:

[TABLE]

Multiplying by $\Delta t$:

[TABLE]

Thus:

[TABLE]

This intermediate $\tilde{\mathbf{u}}^{n+1}$ is not divergence-free.

Step 2: Pressure correction

[TABLE]

Rearranging:

[TABLE]

Step 3: Enforcing incompressibility Taking divergence and setting $\nabla\cdot\mathbf{u}^{n+1}=0$:

[TABLE]

with Neumann boundary condition:

[TABLE]

Step 4: Time-dependent pressure equation We can write:

[TABLE]

Using explicit Euler with $\Delta\tilde{t}=\Delta t/K$:

[TABLE]

The resulting divergence-free $\mathbf{u}^{n+1}$ is then used in the second stage to advance the Cahn–Hilliard equation.

II.10 Time Discretization for the Cahn–Hilliard Equation

In the second stage of the staggered scheme, we update the phase field $\phi$ using the velocity field $\mathbf{u}^{n+1}$ computed from the Navier–Stokes step. The chemical potential at time step $n$ is given by:

[TABLE]

Step 1: Explicit Euler update of the CH equation We discretize the advective Cahn–Hilliard equation:

[TABLE]

Multiplying through by $\Delta t$:

[TABLE]

Step 2: Isolating $\phi^{n+1}$ Rearranging:

[TABLE]

Here:

• The first term $\phi^{n}$ represents the phase field from the previous time step.
• The second term accounts for advective transport by $\mathbf{u}^{n+1}$.
• The third term describes diffusive relaxation driven by the chemical potential gradient.

Step 3: Time-step restriction Due to the explicit treatment of the fourth-order spatial derivative in the Cahn–Hilliard equation, the time step must satisfy Eyre (1998):

[TABLE]

where $\Delta t$ is the time-step size, $h$ is the spatial grid spacing, $m$ is the mobility coefficient, $\varepsilon$ is the diffuse-interface thickness parameter, and $C$ is a stability constant of order unity.

In practice, this stability restriction is often less severe than it appears. The interface thickness $\varepsilon$ is very small: for realistic fluid–fluid interfaces, the physical thickness is on the order of nanometers (except near the critical point). Directly resolving such thin interfaces in numerical simulations becomes infeasible once the characteristic domain length exceeds the micrometer scale. To overcome this, many phase-field simulations employ a numerically broadened interface Jacqmin (1999), which retains the correct macroscopic dynamics while permitting feasible grid sizes and time steps.

II.11 Space Discretization

The time-discrete scheme presented in the previous section must be further discretized in space to enable numerical implementation. We employ a finite difference method (FDM) on an equidistant rectangular grid with uniform spacing $h$ in both $x$- and $y$-directions.

Although the discussion below is restricted to two dimensions for clarity, the method extends naturally to three dimensions without loss of generality.

The computational domain is defined as:

[TABLE]

which is subdivided into $N_{x}\times N_{y}$ square cells:

[TABLE]

This structured arrangement provides a simple, memory-efficient, and highly parallelizable framework for solving the coupled Cahn–Hilliard–Navier–Stokes system.

i Scalar Variables on Cell Centers

Scalar variables — such as the phase field $\phi$, chemical potential $\mu$, and pressure $p$ — are stored at the centers of the computational cells. This collocation allows straightforward computation of Laplacians and other scalar operators. To handle boundary conditions consistently, we extend the domain with one layer of ghost points outside each physical boundary.

The discrete representation of a scalar field $\phi(x,y)$ is given by:

[TABLE]

where $\phi$ denotes a generic scalar quantity and the index ranges include the ghost points.

ii Velocity Variables on a Staggered Grid

When solving incompressible flow problems, a direct collocation of pressure and velocity can lead to odd–even decoupling and spurious checkerboard patterns in the pressure field. To avoid these numerical artifacts, we adopt a staggered Cartesian grid, also known as the Marker-and-Cell (MAC) grid Ferziger et al. (2019); Harlow et al. (1965).

In this arrangement:

• •

The $u_{0}$ (or $u$) velocity component — aligned with the $x$-direction — is stored at the midpoints of vertical cell faces.

• •

The $u_{1}$ (or $v$) velocity component — aligned with the $y$-direction — is stored at the midpoints of horizontal cell faces.

This placement allows for natural, centered differences when computing the divergence of velocity or the gradient of pressure, leading to improved accuracy and stability.

Formally, the discrete velocity locations are:

[TABLE]

The index ranges include ghost points for the velocity variables, which simplify the treatment of velocity boundary conditions. Ghost points are artificial nodes introduced just outside the physical domain that allow uniform indexing of interior and boundary nodes. They enable the straightforward imposition of Dirichlet, Neumann, or periodic boundary conditions without special-case handling, and they eliminate conditional logic in stencil operations. The additional memory required is negligible compared with the benefits in implementation efficiency and clarity.

In this work, we also employ a staggered, or Marker-and-Cell (MAC) as shown in Figure 4, grid arrangement Harlow et al. (1965), where scalar quantities ($\phi,\mu,p$) are stored at cell centers, while velocity components are defined on cell faces ($u_{0}$ at vertical faces and $u_{1}$ at horizontal faces). This configuration offers several advantages: it improves stability by suppressing spurious pressure oscillations, ensures a more accurate coupling between divergence and gradient operators for enhanced mass conservation, and yields better-conditioned linear systems for pressure correction Guermond et al. (2006). These properties are particularly important in the coupled Cahn–Hilliard–Navier–Stokes framework, where strong interaction between velocity, pressure, and interface dynamics is essential for both physical fidelity and numerical robustness. This spatial discretization framework forms the foundation for our numerical method: scalar variables at cell centers ensure accurate phase field and pressure calculations, while staggered velocities at cell faces enable stable and physically consistent fluid–structure interaction.

In the following section, we describe the finite difference stencils used to approximate spatial derivatives on this grid.

II.12 Constant Density and Viscosity

i Space Discretization of the Momentum Equation

We begin with the governing equation for the intermediate velocity field:

[TABLE]

Under the assumption of constant density and viscosity, $\eta(\phi)=\eta$ and $\rho(\phi)=\rho$, the equation becomes:

[TABLE]

Here, $\mathbf{u}=[u_{0},u_{1}]$ is the velocity vector in 2D. We next simplify Term (ii).

Rate-of-deformation tensor.

The rate-of-deformation tensor $D(\mathbf{u})$ is defined as:

[TABLE]

In Cartesian coordinates:

[TABLE]

The components may also be written as:

[TABLE]

Divergence of $D(\mathbf{u})$.

Taking the divergence gives:

[TABLE]

Equivalently:

[TABLE]

For incompressible flow, $\nabla\cdot\mathbf{u}=0$.

[TABLE]

This implies

[TABLE]

Substituting, the $x$-component becomes:

[TABLE]

and the $y$-component:

[TABLE]

Thus:

[TABLE]

Relation to the Laplacian.

For an incompressible velocity field $\mathbf{u}=(u_{0},u_{1})^{T}$, the vector Laplacian is

[TABLE]

Using the incompressibility condition, it follows that

[TABLE]

Therefore, when the viscosity is constant,

[TABLE]

Discrete Laplacian for $u_{0}$.

[TABLE]

Discrete Laplacian for $u_{1}$.

[TABLE]

Term (iii) – Chemical potential term.

In $x$-direction:

[TABLE]

In $y$-direction:

[TABLE]

Term (i) – Convective term.

The continuous form:

[TABLE]

[TABLE]

Final assembly

For each component $\alpha\in\{0,1\}$, the update is:

[TABLE]

ii Space Discretization of the Pressure Computation Equation

We begin with the semi-implicit pressure update equation:

[TABLE]

where $k=0,\dots,K-1$ denotes the sub-iterations within a single time step.

For a homogeneous system, $\rho(\phi)=\rho$, $\Delta t=dt$, and $\Delta\tilde{t}=\frac{\rho h^{2}}{8}$, which yields:

[TABLE]

Since $\rho$ is constant, the diffusion term simplifies as:

[TABLE]

Substituting, we obtain:

[TABLE]

Spatial discretization.

The Laplacian and divergence operators are expanded in two dimensions as:

[TABLE]

Including the $\rho$ and $dt$ scaling for the divergence term:

[TABLE]

Finite difference form.

The second-order central difference approximation for the Laplacian is:

[TABLE]

The forward-difference approximation for the divergence is:

[TABLE]

Hence:

[TABLE]

Fully discretized pressure update.

Combining all terms, we have:

[TABLE]

Velocity correction step.

After the pressure update, the divergence-free velocity field is computed as:

[TABLE]

With:

[TABLE]

The corrected velocity at each grid point is:

[TABLE]

In compact vector form:

[TABLE]

iii Space Discretization of the Chemical Potential

The chemical potential is expressed as:

[TABLE]

where $\phi$ is the phase-field order parameter and $\epsilon$ is a characteristic length scale.

Laplacian operator in two dimensions.

On a uniform grid with spacing $h$ in both $x$ and $y$ directions, the Laplacian can be written as:

[TABLE]

Using the central difference approximation:

[TABLE]

Thus, the discrete Laplacian is:

[TABLE]

Discretized chemical potential.

Substituting the above discrete Laplacian into the expression for $\mu$ gives:

[TABLE]

This final expression represents the fully discretized chemical potential at grid point $(i,j)$, incorporating both the diffusive term (scaled by $\epsilon^{2}$) and the local nonlinear term $\phi(\phi^{2}-1)$.

iv Space Discretization of the Phase Field

The evolution of the phase field (order parameter) $\phi$ is governed by:

[TABLE]

where the flux $\mathcal{J}$ is given by:

[TABLE]

with $M$ denoting the mobility and $\mu$ the chemical potential.

Expanded form.

Substituting the flux expression into Eq. (176) yields:

[TABLE]

which, for constant $M$, simplifies to:

[TABLE]

(a) Time derivative term.

Using a forward difference in time:

[TABLE]

(b) Advection terms.

Using central differences in space:

[TABLE]

(c) Diffusion term.

[TABLE]

Combined discretization.

[TABLE]

Isolating $\phi_{i,j}^{n+1}$.

[TABLE]

Compact final form.

[TABLE]

II.13 Spatial Discretization: Non-Matched Properties

In the general case of non-matched density and viscosity, the discretization of the velocity update equation becomes more elaborate. We define:

[TABLE]

Starting from the velocity predictor:

[TABLE]

Viscous term for non-homogeneous viscosity.

The rate-of-deformation tensor is:

[TABLE]

In Cartesian coordinates:

[TABLE]

Multiplying by $2\eta$, we obtain:

[TABLE]

Divergence of the viscous stress.

The $x$-component is:

[TABLE]

The $y$-component is:

[TABLE]

Full tensor form.

[TABLE]

Spatial discretization — $x$-component.

For the term $2(\partial_{x}\eta_{ij}\,\partial_{x}u_{0}+\eta_{ij}\,\partial^{2}_{x}u_{0})$:

[TABLE]

Analogous formulas hold for the $y$-derivatives and cross-derivatives.

The discretized $x$-component becomes:

[TABLE]

Spatial discretization — $y$-component.

Similarly:

[TABLE]

II.14 Pressure Correction for Variable Density and Viscosity

In the case of non-homogeneous density $\rho(\mathbf{x},t)$ and viscosity, the pressure equation must be modified to account for spatial variations in these quantities. The iterative pressure update at sub-step $k$ of the fractional-step method is written as:

[TABLE]

i Expansion of Term (i)

Multiplying Term (i) by $\rho^{n}$, we write:

[TABLE]

Expanding the divergence operator:

[TABLE]

Using the product rule, we obtain:

[TABLE]

Discrete Laplacian $\Delta p^{n+\frac{k}{K}}$:

[TABLE]

Gradient term $\rho^{n}\nabla(1/\rho^{n})\cdot\nabla p^{n+\frac{k}{K}}$:

Central differences for $\nabla(1/\rho^{n})$:

[TABLE]

Central differences for $\nabla p$:

[TABLE]

Combining these:

[TABLE]

Complete discrete form of Term (i):

[TABLE]

ii Expansion of Term (ii)

We have:

[TABLE]

Discrete form:

[TABLE]

Thus:

[TABLE]

The contribution to the pressure equation is then:

[TABLE]

II.15 Inhomogeneous Density Correction in Divergence-Free Velocity Evaluation

When density varies spatially, the velocity correction step becomes:

[TABLE]

In discrete form:

[TABLE]

Replacing $\rho$ by local averages:

[TABLE]

Where:

[TABLE]

Explicitly:

[TABLE]

III Results & Discussion

Using the numerical framework outlined in Section II, we investigate two benchmark problems: (a) the rise of a single bubble in a non-homogeneous liquid, and (b) the Plateau–Rayleigh instability. These cases were selected because they test different aspects of the phase-field formulation: interfacial deformation under buoyancy in confined geometries, and capillary-driven instability in slender liquid structures.

III.1 Rising Bubble with Density and Viscosity Contrast

The rise of bubbles in immiscible fluids is a canonical problem in multiphase flow, with relevance to microfluidic devices, chemical reactors, and natural systems such as gas release from sediments. In confined microchannels, bubble motion is governed by three competing effects: (a) Buoyancy, which drives upward motion, (b) Surface tension, which resists deformation and promotes minimal interface curvature, and (c) Viscous forces, which damp motion and can induce asymmetric deformation near channel walls.

In this work, we investigate these interactions using our coupled Cahn–Hilliard–Navier–Stokes solver, validated against benchmark rising-bubble tests, to capture the balance of buoyancy, surface tension, and viscous stresses in confined geometries.

Physical setup and initial conditions.

The computational domain is a rectangular region $\Omega=[1,0]\times[2,0]$ with uniform grid spacing $h$ in both $x$ and $y$ directions. The domain is initially filled with liquid ($\phi\approx-1$), and a bubble of radius $0.25$ is centered at $(0.5,\,0.5)$ with $\phi\approx 1$. The initial phase-field distribution is given by:

[TABLE]

ensuring a smooth diffuse interface of width proportional to $\epsilon$.

Material properties as shown in Figure 3 are defined as:

[TABLE]

where the subscripts $-1$ and $+1$ refer to the outer liquid and the inner vapor phase, respectively. The gravitational force is specified as:

[TABLE]

and the surface tension coefficient is $\sigma=24.5$. Boundary conditions are no-slip at the top and bottom walls and free-slip at the lateral boundaries, mimicking a vertically confined microchannel.

Simulation results and flow features.

Figure 5 (top row, i–iv) shows the evolution of the bubble interface as $\phi$ transitions smoothly between $\phi=1$ (vapor) and $\phi=-1$ (liquid). At $t=0$, the bubble is perfectly circular. As buoyancy accelerates the lighter vapor phase upward, deformation begins:

• Step 250: Slight elongation along the vertical axis, with the lower interface flattening due to upward acceleration.
• Step 750: The bubble exhibits a pronounced dome-like shape, with increased curvature near the top and compression at the bottom.
• Step 1000: The deformation becomes more substantial, reflecting the balance between buoyancy-driven stretching and surface tension’s restorative effect.

The bottom row (v–viii) in Figure 5 presents velocity vectors overlaid on vorticity contours. Two symmetric counter-rotating vortices form in the bubble’s wake soon after motion begins. These vortices grow in size and intensity before stabilizing, generating a low-pressure zone behind the bubble. This vortex-induced pressure drop modifies the drag force, affecting the rise velocity and shape evolution. The gradual change in vortex strength over time illustrates the interplay between inertia and viscous dissipation in confined two-phase flows.

Validation with experimental data.

Figure 6 compares the simulated bubble shape with experimental microchannel images Bhaga and Weber (1981). Both show a dome-shaped bubble with similar curvature, aspect ratio, and wall proximity effects. The agreement demonstrates that:

1. The diffuse-interface formulation accurately captures interfacial deformation and capillary effects in confined geometries.
2. The coupled hydrodynamic–phase-field approach reproduces wake-induced flow patterns consistent with experimental visualization.
3. The simulation maintains numerical stability and mass conservation over long timescales, essential for predictive multiphase modeling.

Overall, this test confirms the model’s robustness for simulating buoyancy-driven interfacial dynamics in confined microscale channels.

III.2 Rayleigh–Taylor Instability

We now apply the phase-field formulation to the Rayleigh–Taylor instability (RTI), a canonical benchmark for studying buoyancy-driven interfacial instabilities in multiphase flows. RTI occurs when a denser fluid is placed above a lighter fluid in a gravitational field, making the configuration unstable. Small perturbations at the interface grow over time as the heavier fluid descends into the lighter one (spikes) and the lighter fluid rises into the heavier one (bubbles). This instability is relevant to a wide range of natural and engineering processes, including oceanic mixing, astrophysical flows (e.g., supernova remnants), and inertial confinement fusion.

Physical mechanism.

The instability is driven by the gravitational potential energy difference between the two fluid layers. Any small perturbation to the initially flat interface is amplified because gravity accelerates the heavier fluid downward and the lighter fluid upward. Surface tension can act to stabilize short-wavelength disturbances, while viscosity modifies the growth rate and morphology.

Numerical setup.

The computational domain is a rectangular box of size $\Omega=[0,L_{x}]\times[0,L_{y}]$, initially filled with two immiscible fluids. The initial phase-field distribution is given by:

[TABLE]

where $y_{0}$ denotes the initial interface height, and the sinusoidal perturbation of amplitude $0.1$ seeds the instability. The density and viscosity are defined as $\rho_{+1}$, $\eta_{+1}$ for the heavy fluid ($\phi\approx 1$) and $\rho_{-1}$, $\eta_{-1}$ for the light fluid ($\phi\approx-1$). The gravitational force is applied in the negative $y$-direction, $\mathbf{f}=(0,-g(\rho(\phi)-\rho_{-1}))$, with $g$ the gravitational acceleration. No-slip boundary conditions are imposed at the top and bottom, and free-slip boundaries at the sides. The coupled Navier–Stokes–Cahn–Hilliard system (Section II) governs the evolution.

Simulation results.

Figure 7 presents the time evolution of the instability in six snapshots (a)–(f). Initially, the interface is nearly flat, with only the imposed sinusoidal perturbation visible. As time progresses, the heavier fluid penetrates downward in spike-like structures, while the lighter fluid rises upward in bubble-like protrusions. At early times, the perturbation amplitude grows exponentially, consistent with linear Rayleigh–Taylor theory Rayleigh (1882); Taylor (1950). As the system enters the nonlinear regime, spikes sharpen and accelerate under gravity, while bubbles broaden as they rise Youngs (1984). The shear that develops along the flanks of spikes and bubbles generates counter-rotating vortices, which in turn intensify and promote the formation of a mixing zone where small-scale interpenetration of the fluids occurs Youngs (1984).

Discussion.

The phase-field method captures both the early linear growth and the highly nonlinear deformation without requiring explicit interface tracking. The diffuse-interface representation ensures smooth curvature evaluation and accurate surface tension forces. The simulated patterns, including spike and bubble morphology, agree qualitatively with experimental and theoretical predictions of RTI evolution. This validates the capability of the present formulation to model buoyancy-driven instabilities with strong density and viscosity contrasts.

IV Conclusion

In this study, we employed the phase-field method within the framework of the coupled Navier–Stokes–Cahn–Hilliard (NSCH) systemJacqmin (1999) to investigate the dynamics of immiscible two-phase flows. Two representative benchmark problems were considered: (i) the rise of a bubble in a non-homogeneous liquid, and (ii) the Rayleigh–Taylor instability between two fluids of different densities.

For the rising-bubble case, the simulations accurately reproduced the full sequence of dynamic behavior, from the initial spherical shape to progressive deformation under the competing influences of buoyancy, surface tension, and viscous drag. The formation and evolution of counter-rotating vortices in the wake, as well as the gradual shape transition to a dome-like interface, were captured with high fidelity, underscoring the method’s ability to resolve subtle interfacial curvature effects in confined geometries.

For the Rayleigh–Taylor instability, the model reproduced both the qualitative morphology and the quantitative growth rates predicted by linear stability theory Youngs (1984) during the early stages of instability development. The dimensionless growth rates, extracted from interface amplitude evolution, were consistent with theoretical predictions and benchmark numerical studies Youngs (1984). In the nonlinear regime, the simulation captured the emergence of roll-up and finger-like structures, in agreement with high-resolution studies of RTI dynamics, thereby demonstrating the method’s robustness in handling topological changes and interface reconnection without explicit interface tracking.

The NSCH framework proved to be a versatile and robust approach for modeling immiscible two-phase flows with complex interface dynamics, large deformations, and topological changes. Its diffuse-interface formulation naturally accommodates interfacial curvature effects, breakup, and coalescence, making it well-suited for multiphysics applications where phase separation couples strongly with fluid motion. While this work focused on isothermal conditions, the methodology can be extended to incorporate thermal effects, enabling the investigation of thermally driven phase transformations and boiling phenomena. Such extensions are particularly relevant to microelectronics cooling Lorenzini and Joshi (2019, 2024), where localized heating induces bubble nucleation and departure, significantly affecting heat transfer efficiency. The present formulation, with appropriate coupling to an energy equation, offers a direct path toward predictive simulations of such thermofluidic processes.

In summary, this study demonstrates that the phase-field method, implemented within the NSCH system, is capable of delivering both qualitative and quantitative agreement with established theory and experiments for immiscible two-phase flows. The findings not only validate the computational framework but also position it as a powerful tool for future research into multiphysics problems involving interfacial dynamics, thermal effects, and complex geometries relevant to industrial and microscale applications.

V Appendix

Taylor Series Expansion

The Taylor series expansion of a scalar function $f$ about a point $x_{0}$ is given by

$\displaystyle f(x)$ \displaystyle=f(x_{0})+\frac{\partial f}{\partial x}\Big{|}_{x=x_{0}}(x-x_{0})

\displaystyle\quad+\frac{1}{2!}\frac{\partial^{2}f}{\partial x^{2}}\Big{|}_{x=x_{0}}(x-x_{0})^{2}+\cdots

(137)

For a multivariable function expanded about the point $\mathbf{x_{0}}=(x_{1,0},x_{2,0},\ldots,x_{n,0})$, we have

$\displaystyle f(\mathbf{x})$ \displaystyle=f(\mathbf{x_{0}})+\sum_{i}\frac{\partial f}{\partial x_{i}}\Big{|}_{\mathbf{x}=\mathbf{x_{0}}}(x_{i}-x_{i,0})

\displaystyle\quad+\frac{1}{2!}\sum_{i,j}\frac{\partial^{2}f}{\partial x_{i}\partial x_{j}}\Big{|}_{\mathbf{x}=\mathbf{x_{0}}}(x_{i}-x_{i,0})(x_{j}-x_{j,0})

$\displaystyle\quad+\cdots$

(138)

Here, $\mathbf{x}=(x_{1},x_{2},\ldots,x_{n})$ denotes the vector of variables and $\mathbf{x_{0}}=(x_{1,0},x_{2,0},\ldots,x_{n,0})$ is the expansion point.

V.1 Expression for Free Energy Functional

This section outlines the development of a general equation for free energy in nonuniform systems, crucial for modeling complex interfacial phenomena. Such systems exhibit spatial variations in the phase field, and the local free energy density $f$ depends not only on the local phase field variable $\phi$ but also on its gradients. To capture this, we decompose $f$ into two contributions: (1) a function of the local phase field variable, and (2) functions of its derivatives. This approach ensures both compositional and gradient effects are represented. To formalize this idea, we expand $f$ in a Taylor series around the reference state $f_{0}(\phi)$, corresponding to a uniform phase field ($\nabla\phi=0$) Cahn (1961).

Taylor Expansion of Free Energy

[TABLE]

Coefficient Definitions

To simplify notation, we define:

[TABLE]

Final Expression

Using these definitions, the expansion of $f$ is:

[TABLE]

Symmetry Constraints

In general, $\kappa_{ij}^{(1)}$ and $\kappa_{ij}^{(2)}$ are tensors reflecting the symmetry of the fluid mixture, while $L_{i}$ represents gradient-dependent contributions. For an isotropic and homogeneous binary mixture, the system must remain invariant under spatial reflections and rotations.

Consider a reflection about the $x_{i}$-axis:

[TABLE]

The first-order term transforms as:

[TABLE]

For invariance, this requires

[TABLE]

Thus, gradient-linear contributions vanish by symmetry, and the leading correction to $f$ arises from quadratic gradient terms.

i Derivation of $\kappa_{ij}^{(1)}$

Symmetry Conditions on $\kappa_{ij}^{(1)}$

The coefficients $\kappa_{ij}^{(1)}$ and $\kappa_{ij}^{(2)}$ must respect the symmetry and isotropy of a binary fluid mixture. These requirements ensure that the free energy functional remains invariant under spatial transformations such as rotations and reflections.

The term involving $\kappa_{ij}^{(1)}$ in the Taylor series expansion is given by:

[TABLE]

For an isotropic and homogeneous binary fluid mixture, invariance under rotations and reflections implies that:

[TABLE]

Definition of $\kappa_{1}$

The isotropic coefficient $\kappa_{1}$ is defined as the derivative of the free energy density with respect to the Laplacian of $\phi$, evaluated at the reference state:

[TABLE]

Therefore, the conditions for $\kappa_{ij}^{(1)}$ can be summarized compactly as:

[TABLE]

ii Derivation of $\kappa_{ij}^{(2)}$

The Term Involving $\kappa_{ij}^{(2)}$

The term involving $\kappa_{ij}^{(2)}$ in the Taylor series expansion is:

[TABLE]

Similar to $\kappa_{ij}^{(1)}$, the coefficients $\kappa_{ij}^{(2)}$ must be identical for all diagonal elements and vanish for off-diagonal elements in order to maintain invariance under rotation and reflection.

Definition of $\kappa_{2}$

The coefficient $\kappa_{2}$ is defined as the second derivative of the free energy density with respect to $(\partial\phi)^{2}$, evaluated at the reference state:

[TABLE]

Thus, the conditions for $\kappa_{ij}^{(2)}$ are:

[TABLE]

Free Energy Expression

Integrating over a volume $V$ of the solution, we obtain the total free energy:

[TABLE]

Hence, for a binary fluid mixture, the free energy density reduces to:

[TABLE]

Simplifying the Integral Term

Consider the contribution:

[TABLE]

Since $\kappa_{1}$ is a function of $\phi$, we apply the product rule:

[TABLE]

Substituting into Eq. (157), we obtain:

[TABLE]

Divergence Theorem

The divergence theorem (Gauss’s theorem) relates the flux of a vector field through a closed surface to the divergence of the field within the enclosed volume:

$\int_{V}(\nabla\cdot\mathbf{A})\,dV=\int_{S}(\mathbf{A}\cdot\mathbf{n})\,dS,$

where $V$ is the volume, $S$ is the enclosing surface, $\mathbf{A}$ is a vector field, and $\mathbf{n}$ is the outward normal unit vector on $S$.

Applying this to the first term yields:

[TABLE]

Thus,

[TABLE]

Simplifying the Second Term

Since $\kappa_{1}=\kappa_{1}(\phi)$:

[TABLE]

which gives:

[TABLE]

Substituting back, we obtain:

[TABLE]

Finally, neglecting boundary effects by choosing $\nabla\phi\cdot\mathbf{n}=0$ at the boundary, the surface integral vanishes, yielding:

[TABLE]

Final Expression for the Free Energy Functional

Substituting this back into Eq. (155), we obtain:

[TABLE]

where

[TABLE]

Equation (166) is pivotal: it shows that, to first approximation, the free energy of a small volume of a nonuniform solution can be expressed as the sum of two contributions—the homogeneous free energy density $f_{0}$ and a “gradient energy” term that depends on the local order parameter field $\phi$.

V.2 Advective Cahn–Hilliard Formulation

In this section, we derive the Cahn–Hilliard equation in the presence of fluid motion by incorporating an advective term. This extended form, denoted as Eq. 176 in Section II.3, is essential when modeling mixtures of incompressible and immiscible fluids where transport occurs not only by diffusion but also by bulk advection. To illustrate, we consider two fluids, A and B, with distinct densities and viscosities. Their local composition within a control volume can be represented through a volume fraction approach, akin to the Volume of Fluid (VOF) method. Specifically, we use the volume fraction of fluid A, denoted by $C$ ($0\leq C\leq 1$), as the primary descriptor of the mixture.

Density definitions.

Given $C$, the local densities of fluids A and B in a volume element are:

[TABLE]

where $\rho_{A}$ and $\rho_{B}$ are the bulk densities. The average density is thus:

[TABLE]

Mass conservation.

To establish the continuity equation, consider the mass of species A within a control volume $V$:

[TABLE]

Its time rate of change is

[TABLE]

Applying the divergence theorem, the rate of change must balance the net flux across the surface $S$ bounding $V$:

[TABLE]

where $\mathbf{n}_{A}$ is the flux of species A and $\mathbf{n}$ the outward normal. Since $V$ is arbitrary:

[TABLE]

Flux decomposition.

The total flux of species A comprises advective and diffusive contributions:

[TABLE]

where $\mathbf{u}$ is the fluid velocity and $\mathbf{j}_{A}$ is the diffusive flux. Neglecting diffusion ($\mathbf{j}_{A}=0$) reduces $\mathbf{n}_{A}$ to its advective form.

Continuity equation for volume fraction.

Substituting $\tilde{\rho}_{A}=C\rho_{A}$ into Eq. (170) gives

[TABLE]

Normalizing by $\rho_{A}$ yields

[TABLE]

for component A.

Similarly, for fluid B we obtain

[TABLE]

for component B.

Divergence-free velocity condition.

Adding Eqs. (172) and (173) gives

[TABLE]

Since interfacial diffusion ensures

[TABLE]

we obtain

[TABLE]

which enforces incompressibility of the velocity field.

Final advective Cahn–Hilliard equation.

The governing evolution equation for the volume fraction $C$ becomes

[TABLE]

where $\mathcal{J}=\mathbf{j}_{A}=-\mathbf{j}_{B}$ is the diffusive flux.

By replacing the volume fraction $C$ with the order parameter $\phi$, which ranges from $-1$ (pure B) to $+1$ (pure A), Eq. (175) transforms into the convective Cahn–Hilliard equation:

[TABLE]

This formulation is general and remains valid even for fluids with different densities and viscosities, since these contrasts do not alter the fundamental structure of the Cahn–Hilliard equation. The order parameter thus serves as a natural descriptor of the fluid mixture, linking advective transport with interfacial diffusion in a unified framework.

V.3 Derivation of Linear Interpolation Formulae

i Viscosity interpolation

For viscosity, set $(x_{0},y_{0})=(-1,\eta_{1})$ and $(x_{1},y_{1})=(1,\eta_{2})$:

[TABLE]

Replacing $\phi$ with $\hat{\phi}$ gives Eq. (31)

Linear Interpolation

Linear interpolation estimates a value between two known points $(x_{0},y_{0})$ and $(x_{1},y_{1})$ using:

$y=y_{0}+\frac{x-x_{0}}{x_{1}-x_{0}}\,(y_{1}-y_{0}).$

ii Density interpolation

Similarly, for density $(\rho_{1},\rho_{2})$:

[TABLE]

which matches Eq. (32) after truncation.

Acknowledgements

This research was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, through the Data, Artificial Intelligence, and Machine Learning at DOE Scientific User Facilities program under Award No. 34532 (Digital Twins). Portions of this work were carried out at the Center for Nanoscale Materials and the Advanced Photon Source, both DOE Office of Science User Facilities supported by the Office of Basic Energy Sciences, under Contracts DE-AC02-06CH11357 and DE-AC05-00OR22725, respectively. Computational resources were provided by the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported under Contract DE-AC02-05CH11231, as well as by the Laboratory Computing Resource Center (LCRC) at Argonne National Laboratory.