Three-Phase Cahn-Hilliard/Allen-Cahn Model

• The three-phase Cahn-Hilliard/Allen-Cahn system is a multiphysics framework that combines mass-conserving diffusive and order parameter dynamics via fourth- and second-order PDEs to model phase separation and coarsening phenomena.
• It employs a variational free-energy formulation with pairwise surface tensions and SPD constraints, ensuring physically consistent predictions and well-posedness of the governing equations.
• Recent advances using physics-informed neural operators (PINOs) demonstrate enhanced computational accuracy and efficiency in simulating complex interfacial dynamics in three-phase mixtures.

The three-phase Cahn–Hilliard/Allen–Cahn system is a multiphysics continuum framework for simulating microstructural evolution in mixtures and alloys containing three distinct phases. It couples mass-conserving and @@@@2@@@@ via fourth- and second-order partial differential equations (PDEs), capturing both diffusive and interfacial processes relevant to mesoscale phase separation and coarsening phenomena.

1. Mathematical Formulation

The core model considers either three incompressible fluids (volume fractions $\mathbf{c}=(c_1,c_2,c_3)^T$, $c_i\ge0$, $\sum_{i=1}^3c_i=1$) or, in metallurgical applications, a composition field $c(\bm x, t)$ and two long-range order parameters $\eta_1(\bm x, t)$, $\eta_2(\bm x, t)$. The coupling arises from a free-energy functional embodying bulk and interfacial thermodynamics as well as elastic effects: $$\mathcal{F}[c,\eta_1,\eta_2] = \int_{\Omega} \big[ f(c,\eta_1,\eta_2) + \kappa_c |\nabla c|^2 + \sum_{i=1}^2 \kappa_{\eta_i} |\nabla \eta_i|^2 + w_{el}(c,\eta_1,\eta_2) \big]\,d\bm x$$ with $f$ typically being a sixth-order polynomial, and $w_{el}$ the elastic energy density.

The gradient-flow dynamics are given by:

• Cahn–Hilliard equation for $c$ (fourth order in space),

$$\frac{\partial c}{\partial t} = M \nabla^2\left(\frac{\partial f}{\partial c}\right) - 2 \kappa_c M \nabla^4 c$$

• Allen–Cahn equations for $\eta_i$ ($i=1,2$, second order in space),

$$\frac{\partial \eta_i}{\partial t} = -L\left[\frac{\partial f}{\partial \eta_i} - 2\kappa_{\eta_i} \nabla^2 \eta_i + \frac{\delta F_{el}}{\delta \eta_i}\right]$$

Here, $M$ denotes mobility, $L$ is the interface kinetics, $\kappa_c$, $\kappa_{\eta_i}$ are gradient-energy coefficients, and $\frac{\delta F_{el}}{\delta \eta_i}$ is computed efficiently in Fourier space with pseudo-spectral methods (Gangmei et al., 24 Jul 2025).

2. Free-Energy Structure and Pairwise Surface Tensions

In the generic three-phase setting, the local free-energy density is constructed using pairwise surface tension parameters $\sigma_{ij}$ and a symmetric matrix $\Lambda$: $$W(\boldsymbol\phi, \nabla\boldsymbol\phi) = \frac{\eta}{2} \nabla\boldsymbol\phi : \Lambda \nabla\boldsymbol\phi + \frac{1}{\eta} F(\mathbf{c})$$ Bulk potential $F(\mathbf{c})$ is formulated so that for two-phase configurations, the correct interfacial energy scaling with $\sigma_{ij}$ is preserved.

The coefficient matrix $\Lambda$ is determined by the system of equations

$$\mathbf{L}_{kl} \mathbf{L}_{kl}^T : \Lambda = \frac{9}{2} \sigma_{kl}$$

for vectors $\mathbf{L}_{kl} = \mathbf{a}_k - \mathbf{a}_l$, encapsulating the geometric effect of surface tensions among phases (Wu et al., 2016).

A crucial property is the symmetric positive-definite (SPD) requirement of $\tilde\Lambda$ (the projection of $\Lambda$ onto the tangent space of concentration variables), which guarantees well-posedness of the PDEs. This SPD property is equivalent to classical triangle inequalities for the $\{\sqrt{\sigma_{ij}}\}$, geometrically ensuring the physicality of the system.

3. Variational Gradient-Flow Equations

The Allen–Cahn and Cahn–Hilliard dynamics arise as $L^2$ and $H^{-1}$ gradient flows of the total free energy, respectively, subjected to mass-conservation and capillarity constraints. The equations for concentrations $\mathbf{c}$ in strong form are: $$\gamma P_c \tilde\Lambda_c \partial_t \mathbf{c} - \nabla \cdot (\eta \tilde\Lambda_c \nabla \mathbf{c}) + \frac{1}{\eta} P_c F_{\mathbf{c}}(\mathbf{c}) = 0$$ and, for the Cahn–Hilliard sector,

$$P_c \partial_t \mathbf{c} = -\nabla \cdot (M_0 \tilde\Lambda_c^\dagger \nabla \mathbf{w}), \qquad P_c \mathbf{w} = -\nabla\cdot(\eta \tilde\Lambda_c\nabla\mathbf{c}) + \frac{1}{\eta} P_c F_{\mathbf{c}}(\mathbf{c})$$

where $P_c$ is the projector to concentrations with vanishing sum and $\tilde\Lambda_c$ is the projected surface tension matrix (Wu et al., 2016).

A fundamental result is invariance under phase-variable change: all physically measurable dynamics (concentration evolution, interface positions) only depend on $\tilde\Lambda_c$, not on the specific choice of auxiliary variables, ensuring coordinate-free physical predictions (Wu et al., 2016).

4. Numerical Methods and Stability Properties

Finite-element discretizations employ piecewise linear subspaces conforming to the appropriate tangent spaces (e.g., $H^1(\Omega; T\Sigma)$), with time-stepping schemes for both Allen–Cahn and Cahn–Hilliard flows. Three classes of time integration are rigorously analyzed:

• Semi-implicit (first order)
• Fully implicit
• Modified Crank–Nicolson

Discrete energy laws are proved under natural CFL-type step size restrictions, with unconditional stability for Crank–Nicolson schemes. All variants uphold monotonic decrease of the discrete free energy, with exact telescoping in the latter case.

The general structure is

$$E(\boldsymbol\phi_h^n) - E(\boldsymbol\phi_h^{n-1}) \le -C \|\mathbf{c}_h^n - \mathbf{c}_h^{n-1}\|^2$$

for appropriate constants, guaranteeing dissipativity (Wu et al., 2016).

5. Physics-Informed Neural Operator Approach

Recent advances demonstrate that Physics-informed Neural Operators (PINOs) can learn and solve three coupled phase-field equations simultaneously. The PINO architecture projects initial $128\times128\times3$ field data into feature space via a neural lift, processes with a stack of Fourier layers (each layer truncates the FFT to a fixed number of modes to compute spatial derivatives), and projects back to physical concentration and order parameter fields.

Fourier-domain differentiation enables rapid and accurate computation of high-order spatial derivatives: $$\partial_x^n u(\bm x) = \mathcal{F}^{-1}\big[(ik_x)^n \widehat{u}(\bm k)\big]$$

A key finding is that pseudo-spectral Fourier differentiation for the Cahn–Hilliard term reduces PDE loss by approximately $12$ orders of magnitude compared to finite-difference schemes, with no additional boundary-condition term required: periodicity is built into the architecture. PINO generalizes well when trained on a modest number of parameter instances, yielding $<6$\% relative $L^2$ error on composition prediction for unseen scenarios (Gangmei et al., 24 Jul 2025).

6. Representative Simulations and Physical Implications

Representative numerical experiments validate both the continuum and PINO frameworks:

• Spinodal decomposition in three-phase mixtures ($160 \times 160 \times 2$ mesh, random initial perturbations) yields symmetric or biased coarsening patterns, manifesting surface tension effects and mass conservation to machine precision (Wu et al., 2016).
• Triple junction evolution in quaternary systems illustrates the geometric role of pairwise tensions, with junction angles and interface topology responding to inhomogeneous parameter choices as predicted by Young’s law.

In PINO-based predictions for $\theta'$ precipitate growth in Al–Cu alloys, Fourier derivatives naturally enforce periodic BCs, maintain stability for high-order Cahn–Hilliard terms, and deliver accurate phase evolution trajectories across a range of supersaturations and initial field seeds (Gangmei et al., 24 Jul 2025).

7. Limitations, Open Questions, and Future Directions

Identified limitations include the lack of observed computational speed-up at moderate ($128\times128$) 2D resolutions for PINOs; anticipated advantages may manifest in larger or 3D domains. Temporal extrapolation and inference at finer resolutions ("zero-shot resolution transfer") remain unvalidated. Extension to complex microstructures with anisotropic tensions is ongoing work.

This suggests that scaling studies and algorithmic innovations are vital for harnessing PINO flexibility at scale. The full mathematical classification of admissible surface tension matrices, as well as phase-space topology transitions for more than three components, continues to be an active area of investigation (Wu et al., 2016).