RasmahRasmah

Cahn–Hilliard Phase Separation

Overview

Many alloys, polymer blends and emulsions separate into two phases — oil and water, or the two species of a binary mixture — without any external forcing. The Cahn–Hilliard equation (Cahn & Hilliard 1958) is the canonical model of this spinodal decomposition: the concentration c obeys a conserved gradient flow of a free energy that penalises both the bulk energy f(c) and the interface energy (ε²/2)|∇c|²:

\[\partial_t c = \nabla \cdot (M \nabla \mu), \qquad \mu = f'(c) - \varepsilon^2 \nabla^2 c, \qquad f(c) = \frac{(c^2 - 1)^2}{4}\]

where M is the mobility, ε the interface width, and μ the chemical potential. Two properties define a valid solution: the total mass ∫c dx is conserved exactly (the flux M∇μ is a divergence), and the free energy decreases monotonically — the flow is a descent in a gradient metric.

Rasmah solves the mixed (c, μ) form with the energy-stable convex-splitting scheme of Eyre (1998): the convex part c⁴/4 is treated implicitly, the concave −c²/2 explicitly.

using Rasmah

Build a square mesh

Start from a structured grid of points on [0,1]² and triangulate it with a Delaunay mesh. This is the domain the two phases will separate on.

n = 30
grid = [[i / n, j / n] for j in 0:n for i in 0:n]
mesh = delaunay_triangulation(reduce(hcat, grid))
println("nodes = ", size(mesh.vertices, 2), ", triangles = ", size(mesh.faces, 2))
nodes = 961, triangles = 1800

Initial perturbation

A perfectly uniform field is a fixed point (the free energy is at a maximum but has zero gradient). To trigger separation we seed a small, zero-mean cosine perturbation. solve(CahnHilliard(), …) advances the conserved flow and returns the concentration field at every step.

c0 = [0.05 * sin(3π * x) * sin(3π * y) for (x, y) in eachcol(mesh.vertices)]

res = solve(CahnHilliard(), mesh, c0; nsteps=150, dt=2e-4, mobility=10.0, ε=0.03)
c_end = res.cs[end]
println("c range: ", round(minimum(c_end); digits=3), " … ", round(maximum(c_end); digits=3))
c range: -1.053 … 0.997

The two phases separate to the equilibrium values c = ±1:

vtk(mesh; field = c_end, fieldname = "concentration c")

Verification — mass conservation and energy decrease

The defining checks are the conserved total mass and the monotonically decreasing free energy F(c) = ∫ [f(c) + (ε²/2)|∇c|²] dx:

sp = FESpace(mesh, ReferenceFE(2, 1))
mass0 = sum(c0) / length(c0)
mass1 = sum(c_end) / length(c_end)
F0 = cahnhilliard_energy(sp, c0, 0.03)
F1 = cahnhilliard_energy(sp, c_end, 0.03)
println("mean c:   ", round(mass0; digits = 6), " → ", round(mass1; digits = 6))
println("free energy: ", round(F0; digits = 4), " → ", round(F1; digits = 4), "  (decreases: ", F1 < F0, ")")
mean c:   0.002074 → 0.004572
free energy: 0.2497 → 0.0974  (decreases: true)

Next steps