GitHub

Cahn-Hilliard Model

Overview

The Cahn-Hilliard equation describes the evolution of a conserved order parameter. It is used for modeling spinodal decomposition, phase separation, and coarsening.

The concentration field c represents the local composition (mole fraction).

Mathematical Formulation

The Cahn-Hilliard equation:

\[\frac{\partial c}{\partial t} = M \nabla^2 \mu\]

where the chemical potential μ is:

\[\mu = \frac{df}{dc} - \kappa \nabla^2 c\]

Parameters:

  • M: Mobility coefficient
  • κ: Gradient energy coefficient
  • f(c): Bulk free energy density (e.g., double-well)

Combining these gives the 4th-order equation:

\[\frac{\partial c}{\partial t} = M \nabla^2 \left( \frac{df}{dc} \right) - M \kappa \nabla^4 c\]

Quick Example

Basic Usage

using PhaseFields

# Create model and free energy
model = CahnHilliardModel(M=1.0, κ=1.0)
f = DoubleWellFreeEnergy(A=1.0, c_eq_α=0.2, c_eq_β=0.8)

# Compute chemical potential
c = 0.5
∇²c = -0.01
μ = cahn_hilliard_chemical_potential(f, model, c, ∇²c)

# Compute concentration change rate
∇²μ = 0.001
dcdt = cahn_hilliard_rhs(model, ∇²μ)

Free Energy Functions

# Double-well free energy: f(c) = A(c - c_α)²(c - c_β)²
f = DoubleWellFreeEnergy(A=1.0, c_eq_α=0.2, c_eq_β=0.8)

# Evaluate
energy = free_energy_density(f, 0.5)
μ_bulk = chemical_potential_bulk(f, 0.5)

Interface Properties

# Estimate interface width
W = cahn_hilliard_interface_width(model, f)

# Stability constraint for explicit time stepping
dt_max = cahn_hilliard_stability_dt(model, f, dx)

API Reference

PhaseFields.CahnHilliardModel — Type
CahnHilliardModel

Cahn-Hilliard model for conserved order parameter (concentration) evolution.

The Cahn-Hilliard equation: ∂c/∂t = ∇·(M∇μ) μ = df/dc - κ∇²c

where:

  • c: Concentration field (conserved quantity)
  • M: Mobility [m²/(J·s)]
  • κ: Gradient energy coefficient [J·m²/mol] or [J/m] depending on normalization
  • μ: Chemical potential [J/mol]
  • f(c): Bulk free energy density [J/mol]

Fields

  • M::Float64: Mobility
  • κ::Float64: Gradient energy coefficient

Example

model = CahnHilliardModel(M=5.0, κ=2.0)
source
PhaseFields.DoubleWellFreeEnergy — Type
DoubleWellFreeEnergy

Symmetric double-well free energy density for binary phase separation.

f(c) = ρs · (c - cα)² · (cβ - c)²

Minima at c = cα (α-phase) and c = cβ (β-phase).

Fields

  • ρs::Float64: Barrier height [J/mol]
  • cα::Float64: α-phase equilibrium concentration
  • cβ::Float64: β-phase equilibrium concentration

Example

# PFHub BM1 parameters
f = DoubleWellFreeEnergy(ρs=5.0, cα=0.3, cβ=0.7)
source
PhaseFields.free_energy_density — Function
free_energy_density(f::DoubleWellFreeEnergy, c::T) where T<:Real

Compute bulk free energy density f(c).

f(c) = ρs · (c - cα)² · (cβ - c)²

This function is AD-compatible.

Arguments

  • f: DoubleWellFreeEnergy parameters
  • c: Concentration

Returns

  • Free energy density [J/mol]

Example

f = DoubleWellFreeEnergy(ρs=5.0, cα=0.3, cβ=0.7)
energy = free_energy_density(f, 0.5)  # At spinodal point
source
PhaseFields.chemical_potential_bulk — Function
chemical_potential_bulk(f::DoubleWellFreeEnergy, c::T) where T<:Real

Compute bulk chemical potential df/dc.

df/dc = 2ρs · (c - cα) · (cβ - c) · (cα + cβ - 2c)

This function is AD-compatible.

Arguments

  • f: DoubleWellFreeEnergy parameters
  • c: Concentration

Returns

  • Bulk chemical potential [J/mol]

Example

f = DoubleWellFreeEnergy(ρs=5.0, cα=0.3, cβ=0.7)
μ_bulk = chemical_potential_bulk(f, 0.5)  # = 0 at c = (cα+cβ)/2
source
PhaseFields.cahn_hilliard_chemical_potential — Function
cahn_hilliard_chemical_potential(model::CahnHilliardModel, f::DoubleWellFreeEnergy,
                                 c::T, ∇²c::T) where T<:Real

Compute total chemical potential including gradient term.

μ = df/dc - κ∇²c

Arguments

  • model: CahnHilliardModel parameters
  • f: DoubleWellFreeEnergy for bulk contribution
  • c: Concentration at current point
  • ∇²c: Laplacian of concentration

Returns

  • Total chemical potential [J/mol]

Example

model = CahnHilliardModel(M=5.0, κ=2.0)
f = DoubleWellFreeEnergy(ρs=5.0, cα=0.3, cβ=0.7)
μ = cahn_hilliard_chemical_potential(model, f, 0.5, -0.01)
source
PhaseFields.cahn_hilliard_rhs — Function
cahn_hilliard_rhs(model::CahnHilliardModel, ∇²μ::T) where T<:Real

Compute the right-hand side of the Cahn-Hilliard equation.

∂c/∂t = M∇²μ

Note: This assumes ∇·(M∇μ) ≈ M∇²μ for constant mobility.

Arguments

  • model: CahnHilliardModel parameters
  • ∇²μ: Laplacian of chemical potential

Returns

  • Time derivative ∂c/∂t

Example

model = CahnHilliardModel(M=5.0, κ=2.0)
dcdt = cahn_hilliard_rhs(model, 100.0)
source
PhaseFields.cahn_hilliard_interface_width — Function
cahn_hilliard_interface_width(model::CahnHilliardModel, f::DoubleWellFreeEnergy)

Estimate equilibrium interface width.

W ≈ √(κ / (ρs · (cβ - cα)²))

Arguments

  • model: CahnHilliardModel parameters
  • f: DoubleWellFreeEnergy parameters

Returns

  • Approximate interface width
source
PhaseFields.cahn_hilliard_stability_dt — Function
cahn_hilliard_stability_dt(model::CahnHilliardModel, dx::Float64)

Estimate maximum stable time step for explicit Euler scheme.

For the 4th-order Cahn-Hilliard equation, stability requires: dt < dx⁴ / (16 · M · κ)

Arguments

  • model: CahnHilliardModel parameters
  • dx: Grid spacing

Returns

  • Maximum stable time step

Example

model = CahnHilliardModel(M=5.0, κ=2.0)
dt_max = cahn_hilliard_stability_dt(model, 1.0)  # ≈ 0.00625
source

See Also