By Evalf, and other Nutils contributors

officialCahn-Hilliard

Solves the Cahn-Hilliard equation, which models the unmixing of two phases (φ=+1 and φ=-1) under the influence of surface tension. It is a mixed equation of two scalar equations for phase φ and chemical potential η:

dφ/dt = -div(J(η))
η = ψ'(φ) σ / ε - σ ε Δ(φ)

along with constitutive relations for the flux vector J = -M ∇η and the double well potential ψ = .25 (φ² - 1)², and subject to boundary conditions ∂ₙφ = -σd / σ ε and ∂ₙη = 0. Parameters are the interface thickness ε, fluid surface tension σ, differential wall surface tension σd, and mobility M.

Cahn-Hilliard is a diffuse interface model, which means that phases do not separate sharply, but instead form a transition zone between the phases. The transition zone has a thickness proportional to ε, as is readily confirmed in one dimension, where a steady state solution on an infinite domain is formed by η(x) = 0, φ(x) = tanh(x / √2 ε).

The main driver of the unmixing process is the double well potential ψ that is proportional to the mixing energy, taking its minima at the pure phases φ=+1 and φ=-1. The interface between the phases adds an energy contribution proportional to its length. At the wall we have a phase-dependent fluid-solid energy. Over time, the system minimizes the total energy:

E(φ) := ∫_Ω ψ(φ) σ / ε + ∫_Ω .5 σ ε ‖∇φ‖² + ∫_Γ (σm + φ σd)
        \                \                  \
         mixing energy    interface energy   wall energy

Proof: the time derivative of E followed by substitution of the strong form and boundary conditions yields dE/dt = ∫_Ω η dφ/dt = -∫_Ω M ‖∇η‖² ≤ 0. □

Switching to discrete time we set dφ/dt = (φ - φ0) / dt and add a stabilizing perturbation term δψ(φ, φ0) to the double well potential for reasons outlined below. This yields the following time discrete system:

φ = φ0 - dt div(J(η))
η = (ψ'(φ) + δψ'(φ, φ0)) σ / ε - σ ε Δ(φ)

For stability we wish for the perturbation δψ to be such that the time discrete system preserves the energy dissipation property E(φ) ≤ E(φ0) for any timestep dt. To derive suitable perturbation terms to this effect, we define without loss of generality δψ'(φ, φ0) = .5 (φ - φ0) f(φ, φ0) and derive the following condition for unconditional stability:

E(φ) - E(φ0) = ∫_Ω .5 (1 - φ² - .5 (φ + φ0)² - f(φ, φ0)) (φ - φ0)² σ / ε
             - ∫_Ω (.5 σ ε ‖∇φ - ∇φ0‖² + dt M ‖∇η‖²) ≤ 0

The inequality holds true if the perturbation f is bounded from below such that f(φ, φ0) ≥ 1 - φ² - .5 (φ + φ0)². To keep the energy minima at the pure phases we additionally impose that f(±1, φ0) = 0, and select 1 - φ² as a suitable upper bound which we will call "nonlinear".

We next observe that η is linear in φ if f(φ, φ0) = g(φ0) - φ² - (φ + φ0)² for any function g, which dominates if g(φ0) ≥ 1 + .5 (φ + φ0)². While this cannot be made to hold true for all φ, we make it hold for -√2 ≤ φ, φ0 ≤ √2 by defining g(φ0) = 2 + 2 |φ0| + φ0², which we will call "linear". This scheme further satisfies a weak minima preservation f(±1, ±|φ0|) = 0.

We have thus arrived at the three stabilization schemes implemented here:

• nonlinear: f(φ, φ0) = 1 - φ²
• linear: f(φ, φ0) = 2 + 2 |φ0| - 2 φ (φ + φ0)
• none: f(φ, φ0) = 0 (not unconditionally stable)

Finally, we observe that the weak formulation:

∀ δη: ∫_Ω [ dt J(η)·∇δη - δη (φ - φ0) ] = 0
∀ δφ: ∫_Ω [ δφ (ψ'(φ) + δψ'(φ, φ0)) σ / ε + σ ε ∇(δφ)·∇(φ) ] = ∫_Γ -δφ σd

is equivalent to the optimization problem ∂F/∂φ = ∂F/∂η = 0, where

F(φ, φ0, η) := E(φ) + ∫_Ω [ .5 dt J(η)·∇η + δψ(φ, φ0) σ / ε - η (φ - φ0) ]

For this reason, the stab enum in this script defines the stabilizing term δψ, rather than f, allowing Nutils to construct the residuals through automatic differentiation using the optimize method.