Phase Field Modelling of Microstructure Evolution: Dendrites, Grain Growth, and Precipitates

Phase field modelling is the dominant computational framework for simulating the evolution of microstructure in metallic alloys during solidification, solid-state transformation, and coarsening. By replacing the physically sharp solid-liquid or phase boundary with a diffuse interface described by a continuous order parameter field, the method eliminates the need to explicitly track moving interfaces, enabling simulation of complex morphologies — dendritic networks, polycrystalline grain structures, and multi-particle precipitate distributions — that are inaccessible to sharp-interface or analytical approaches.

Theoretical Foundations of the Phase Field Method

The phase field method originates in the Ginzburg-Landau theory of second-order phase transitions, extended to describe the non-equilibrium evolution of microstructure by Landau and Khalatnikov and later developed for materials applications by Cahn, Hilliard, Allen, Fix, Langer, and Karma among others. The central idea is to represent the state of the material throughout the simulation domain with one or more continuous field variables, called order parameters or phase fields, that take distinct values in each phase and vary smoothly across the interface over a characteristic length scale δ.

The Free Energy Functional

The total free energy of the system is written as a Ginzburg-Landau functional integrating a bulk free energy density plus a gradient energy penalty that penalises spatial variations in the order parameter, giving the interface a finite energy cost per unit area:

F[φ, c] = ∫ { f_bulk(φ, c, T) + (ε²/2)|∇φ|² } dV

where:
  φ    = order parameter (phase field) — 0 in phase A, 1 in phase B
  c    = composition field (mol fraction)
  T    = temperature (K)
  ε    = gradient energy coefficient (J·m⁻¹ or similar)
  f_bulk = local bulk free energy density (J·m⁻³)

Relationship to interface energy σ and width δ:
  σ = ∫ [ (ε²/2)(dφ/dx)² + f_bulk(φ) ] dx = (ε√W)/(3√2)
  δ = ε / √W   (for double-well height W)

The Allen-Cahn Equation (Non-Conserved Order Parameter)

When the order parameter φ is non-conserved — meaning its spatial integral is not required to remain constant — the phase field evolves by the L2 gradient flow (steepest descent) of the free energy functional, giving the Allen-Cahn equation. This applies to grain boundary motion, martensitic transformation, and antiphase domain coarsening.

∂φ/∂t = −L · δF/δφ

Expanded form:
  ∂φ/∂t = L [ ε²∇²φ − ∂f_bulk/∂φ ]

where:
  L   = interface mobility (m·s⁻¹·J⁻¹)
  ∇²φ = Laplacian of the order parameter
  ε²∇²φ = gradient energy term (smoothing)
  ∂f_bulk/∂φ = bulk driving force

Sharp-interface limit: interface velocity v = Mσκ
  where M = interface mobility, σ = interface energy, κ = local mean curvature

The Cahn-Hilliard Equation (Conserved Composition)

When the evolving field variable is a composition that must be globally conserved (no atoms created or destroyed), the appropriate evolution equation is the Cahn-Hilliard equation, written as a continuity equation with a flux driven by gradients in the variational derivative of the free energy (chemical potential). This governs spinodal decomposition, precipitate nucleation, and coarsening.

∂c/∂t = ∇·[ M(c) ∇(δF/δc) ]

Expanded form:
  ∂c/∂t = ∇·[ M(c) ∇( ∂f_bulk/∂c − ε²∇²c ) ]

where:
  c    = local composition (mol fraction)
  M(c) = composition-dependent mobility (m²·s⁻¹·J⁻¹)
  δF/δc = chemical potential μ

Mobility-diffusivity connection:
  M(c) = D(c) / (∂²f_bulk/∂c²)
where D(c) is the interdiffusion coefficient.

Conservation check: ∫ c dV = const (flux divergence form)

CALPHAD Coupling: Replacing Phenomenological Free Energy with Thermodynamic Databases

Early phase field models used polynomial or double-obstacle free energy functions chosen for mathematical convenience rather than physical accuracy. For quantitative simulation of industrial alloys, these must be replaced with thermodynamically assessed free energy functions from CALPHAD databases. The coupling strategy extracts the Gibbs free energy Gα(T, x) for each phase α as a function of temperature and composition from a database (TCFE10 for steels, TTAL8 for aluminium, TCNI12 for nickel superalloys) and computes the driving pressure for transformation and the equilibrium compositions at the interface by interfacing with a Gibbs energy minimisation engine (ThermoCalc, OpenCalphad, Pycalphad).

Driving Force Extraction

ΔG = Gα(T, xα*) − Gβ(T, xβ*) − ∑_i μ_i (xαⅡ* − xβⅡ*)

Partitioning condition at interface (KKS model):
  x = h(φ)·xβ + [1−h(φ)]·xα   (mixture rule)
  μⅡα(T, xα*) = μⅡβ(T, xβ*)        (equal chemical potential)

where:
  xα*, xβ* = compositions in each phase at the interface
  μⅡ        = chemical potential of species i
  h(φ)     = interpolation function, e.g. φ²(3−2φ) (cubic)

The Kim-Steinbach-Beckermann (KKS) model reformulation is the standard approach for CALPHAD-coupled simulations: it decouples the phase field from the composition field, allowing each phase to have its own local composition consistent with the thermodynamic database, connected by the equal chemical potential constraint. This avoids the artificial solute trapping that arises when a single composition field is interpolated across the interface.

Dendritic Solidification: Morphology and Spacing

Dendritic solidification is the most extensively studied application of the phase field method, motivated by its importance in casting, investment casting of superalloys, and additive manufacturing. The characteristic branched snowflake-like morphology of dendrites arises from the Mullins-Sekerka instability: a planar solidification front is unstable to perturbations when thermal or solutal gradients are insufficient to suppress the capillary-driven tip sharpening.

Interface Anisotropy and Crystal Symmetry

Pure phase field models without anisotropy produce isotropic, seaweed-like microstructures. Crystal anisotropy is incorporated through an orientation-dependent gradient energy coefficient:

ε(θ) = ε_0 [ 1 + δ_4 cos(4θ) ]

where:
  θ   = angle between interface normal and crystal axis [100]
  δ_4 = anisotropy strength (typically 0.02–0.05 for metals)
  ε_0 = mean gradient energy coefficient

3D generalisation for cubic anisotropy:
  a_s(n) = 1 + δ_4 [ (n_x⁴ + n_y⁴ + n_z⁴) − 1/3 ]
  where n = (n_x, n_y, n_z) = interface unit normal

Effect: Higher δ_4 produces sharper dendritic tips aligned
with <100> in FCC/BCC metals. At δ_4 > 1/15, the equilibrium
crystal shape (ECS) develops flat facets.

Primary Dendrite Arm Spacing

The primary dendrite arm spacing (PDAS) λ₁ is set by solidification conditions at the growth front. The Hunt-Kurz scaling and the more rigorous Kurz-Fisher model predict:

λ₁ ≈ C₁ (G·V)⁻¹∕⁴

where:
  G   = thermal gradient at the solidification front (K·m⁻¹)
  V   = solidification front velocity (m·s⁻¹)
  C₁ = alloy-dependent constant (K·m·s)¹∕⁴·m

Kurz-Fisher model (more rigorous):
  λ₁ = 4.3 [ ΔT_0 D Γ / k ]¹∕⁴ (G·V)⁻¹∕⁴

where:
  ΔT_0 = solidification range = m_L C_0 (1/k − 1) (K)
  D     = liquid diffusivity (m²·s⁻¹)
  Γ    = Gibbs-Thomson coefficient = σ_SL / ΔS_f (m·K)
  k     = equilibrium partition coefficient
  m_L   = liquidus slope (K·wt%⁻¹)

Phase field simulations quantitatively reproduce this scaling when the interface width is converged (thin-interface limit), confirming that the method correctly captures the underlying physics. For columnar-to-equiaxed transition (CET) problems, multi-seeded simulations with competitive grain elimination reproduce the grain texture selection observed in electron beam welded and laser powder bed fusion parts.

The Thin-Interface (Karma-Rappel) Formulation

A fundamental numerical challenge is that physical solid-liquid interface widths in metals are on the order of 1–10 nm, while dendrite arm spacings are 1–100 μm. Resolving both simultaneously would require impractically fine grids. Karma and Rappel (1996–1998) showed that the interface can be artificially widened without loss of quantitative accuracy by introducing an anti-trapping current →j_at in the composition evolution equation that exactly cancels the spurious solute trapping that the widened interface would otherwise introduce:

Modified Cahn-Hilliard (with anti-trapping):
  ∂c/∂t = ∇·[ D(φ)∇c + W·a₁(cβ−cα)(∂φ/∂t)(∇φ/|∇φ|) ]

The anti-trapping term:
  →j_at = −W·a₁(cβ−cα)(∂φ/∂t) n̂

where:
  W     = phase field interface width (numerical parameter)
  a₁   = numerical coefficient = 1/(2√2)
  cβ,cα = equilibrium compositions in each phase
  n̂    = interface unit normal = ∇φ/|∇φ|

Effect: Eliminates O(W) spurious solute trapping; allows
use of W up to 20× the capillary length without accuracy loss.

Grain Growth and Recrystallisation

Grain growth in polycrystalline metals is driven by the reduction in total grain boundary area and hence total grain boundary energy. Normal grain growth follows parabolic kinetics: mean grain diameter d scales as d² − d₀² = Kt, where K is a temperature-dependent rate constant controlled by boundary mobility M and specific boundary energy γ. The phase field representation of this process uses a multi-order-parameter approach in which each grain is assigned its own scalar field φ_i.

Multi-Order-Parameter Grain Growth Model

F = ∫ { ∑_i [ -α/2 φ_i² + β/4 φ_i⁴ ] + γ ∑_{i≠j} φ_i²φ_j² + ∑_i (ε²/2)|∇φ_i|² } dV

Evolution (Allen-Cahn for each grain field):
  ∂φ_i/∂t = −L [ −αφ_i + βφ_i³ + 2γφ_i ∑_{j≠i} φ_j² − ε²∇²φ_i ]

Parabolic grain growth law:
  <d>² − <d>_0² = K(T)·t

Arrhenius rate constant:
  K(T) = K_0 exp(−Q_gb / RT)

Typical Q_gb values:
  Pure Fe:      150–170 kJ·mol⁻¹
  Austenitic SS: 280–320 kJ·mol⁻¹
  Cu:           100–120 kJ·mol⁻¹

Zener Pinning by Second-Phase Particles

Second-phase particles (oxides, nitrides, carbides) exert a pinning pressure on moving grain boundaries that opposes the capillary driving force. The Zener pinning pressure per unit boundary area is:

P_Z = (3γ_gb f_v) / (2r_p)

Limiting (pinned) grain size — Zener relation:
  d_lim = (4r_p) / (3f_v)

where:
  γ_gb = specific grain boundary energy (J·m⁻²)
  f_v   = volume fraction of pinning particles
  r_p   = mean particle radius (m)

In phase field: pinning particles are represented as
fixed φ = 1 regions that cannot evolve; their curvature
interaction with the grain boundary interface naturally
reproduces Zener pressure without explicit force terms.

Aluminium alloys processed by annealing rely on dispersoid control (Al₃Zr, Al₃Sc particles) to pin subgrain boundaries during recovery and restrict recrystallisation. Phase field simulations of these systems correctly predict the strong particle radius dependence of the pinned grain size.

Precipitate Nucleation, Growth, and Coarsening

The evolution of second-phase precipitates in metallic alloys — strengthening phases, deleterious intermetallics, carbides, nitrides — passes through three stages: nucleation, diffusion-controlled growth, and Ostwald ripening (coarsening). Phase field modelling can simulate the growth and coarsening stages quantitatively; nucleation is typically handled by augmenting the deterministic phase field equations with stochastic Langevin noise terms or a classical nucleation theory (CNT) insertion algorithm.

Gibbs-Thomson Effect and Precipitate Coarsening

The Gibbs-Thomson effect describes the elevation of solute solubility around a curved precipitate-matrix interface relative to a flat interface. Smaller precipitates have higher curvature and therefore higher surrounding matrix solute concentration, creating a concentration gradient that drives diffusion from small to large precipitates — the Ostwald ripening mechanism.

Gibbs-Thomson solubility at precipitate radius r:
  c_r = c_∞ exp(2Γ / r) ≈ c_∞ (1 + 2Γ/r)   [for 2Γ/r << 1]

where:
  c_∞ = equilibrium solubility at flat interface
  Γ   = Gibbs-Thomson coefficient = 2V_mσ / RT (m)
  σ   = precipitate-matrix interface energy (J·m⁻²)
  V_m  = molar volume of precipitate (m³·mol⁻¹)

Lifshitz-Slyozov-Wagner (LSW) coarsening law:
  <r>³ − <r>_0³ = K_LSW · t

  K_LSW = (8 D Γ c_∞) / (9 RT)

  where D = diffusivity of rate-controlling species (m²·s⁻¹)

Phase field reproduction: Γ enters via gradient energy term
  σ = ε√W/(3√2)  →  Γ = 2V_mσ / RT

The t^1/3 coarsening law is automatically reproduced by phase field simulations without any explicit coarsening model — it emerges from the combination of the Gibbs-Thomson capillary effect (embedded in the gradient energy term) and diffusional transport (governed by the Cahn-Hilliard equation). This is a key validation test for any phase field implementation.

Precipitate Morphology and Elastic Strain Effects

When a precipitate has a crystal structure or lattice parameter that differs from the matrix, coherency strains arise at the interface. These elastic strain energies compete with interface energy and bulk chemical driving force to determine the equilibrium precipitate morphology:

F_total = F_chemical + F_interface + F_elastic

F_elastic = ½ ∫ C_ijkl ε_ijᵉ ε_klᵉ dV

Coherency strain:
  ε_ijᵉ = ε_ij − ε_ij⁰ = ε_ij − δ_ij · Δa/a

where:
  C_ijkl = elastic stiffness tensor
  ε_ij⁰  = stress-free transformation strain (misfit)
  Δa/a  = lattice misfit (e.g. 0.003 for γ’ in Ni superalloys)

Morphology map:
  Low misfit + low anisotropy: spherical precipitates
  Moderate misfit: cuboidal (γ’ in Ni-Al, Ni-Al-Ti)
  High misfit: plate/disc (Θ phase in Al-Cu, M₂₃C₆ in steels)

The characteristic cuboidal morphology of γ’ (Ni₃Al) precipitates in nickel superalloys, and their alignment along the <100> elastically soft directions, is reproduced by phase field simulations that include the elastic energy functional with an anisotropic stiffness tensor. See the martensite formation article for discussion of transformation strains in steels.

Multi-Component, Multi-Phase Systems

Real engineering alloys are multicomponent (Fe-Mn-Si-C steels, Al-Cu-Mg-Zn alloys, Ni-Al-Cr-Co-Ti superalloys) and may contain three or more thermodynamically stable phases simultaneously. The Steinbach multi-phase field model extends the two-phase formulation to N phases, each described by its own phase field φ_α, with the constraint ∑φ_α = 1 enforced at every point in the domain.

Sum constraint:  ∑_α φ_α = 1   (at all x, t)

Multi-phase free energy:
  F = ∫ { ∑_{α<β} [ σ_αβ ( -(ε²/2)(∇φ_α·∇φ_β) + f_dw(φ_α,φ_β) ) ] + ∑_α h(φ_α)fα(c,T) } dV

Phase field evolution:
  ∂φ_α/∂t = − (M_φ/N) ∑_β [δF/δφ_α − δF/δφ_β]

where:
  σ_αβ = interface energy between phases α and β
  N     = number of phases in contact at the point

Phase Field Software for Industrial Metallurgy

Industrial Applications

Phase field modelling has transitioned from purely academic research to industrial process design in several sectors. The following applications illustrate the current state of industrial deployment:

Casting and Solidification

Automotive aluminium cylinder block castings, aerospace turbine blade castings, and large steel ingots all suffer from macro- and microsegregation arising from non-equilibrium dendritic solidification. MICRESS-based simulations coupled to casting process codes (MAGMA, ProCAST) predict local dendrite arm spacing, microsegregation profiles, and freckle formation risk in directionally solidified components. For Ni superalloys, CALPHAD-coupled phase field predictions of the mushy zone extent guide the design of withdrawal rates in Bridgman furnaces. See the investment casting of superalloys article for process context.

Steel Heat Treatment

Phase field models of austenite decomposition during continuous cooling predict the ferrite start temperature, transformation kinetics, and the resulting volume fractions of ferrite, pearlite, and bainite as a function of cooling rate and alloy composition. This is particularly valuable for linepipe steels and high-strength structural steels where toughness requirements impose tight constraints on phase balance. Coupling with the Fe-C phase diagram and multicomponent CALPHAD databases enables prediction of the effect of microalloying additions (Nb, V, Ti) on transformation kinetics.

Additive Manufacturing

Laser powder bed fusion (L-PBF) and directed energy deposition (DED) involve rapid solidification at cooling rates of 10⁵–10⁷ K/s, far from equilibrium, generating columnar grain textures, cellular solidification microstructures, and solute segregation patterns that strongly influence mechanical properties. Phase field simulations — particularly those using the Karma-Rappel thin-interface formulation — predict the grain morphology, cell spacing, and microsegregation as functions of laser parameters (power, scan speed, spot size), enabling computationally guided process parameter selection without exhaustive experimental trials.

Welding Metallurgy and HAZ

The heat-affected zone (HAZ) in steel welds undergoes rapid thermal cycling that drives austenite grain growth, dissolution of carbides and nitrides, and subsequent decomposition on cooling. Multi-order-parameter phase field models of grain growth in the coarse-grained HAZ, incorporating Zener pinning by undissolved TiN particles, predict the grain size distribution as a function of peak temperature and holding time — an input to HAZ toughness predictions via the Hall-Petch relationship. For hydrogen-induced cracking susceptibility assessment, phase field models of martensite formation, combined with CALPHAD-predicted martensite start temperature, provide the local hardness distribution that governs cracking risk.

Creep and High-Temperature Service

In service at elevated temperatures, strengthening precipitates coarsen, rafting of γ’ occurs in single-crystal nickel superalloys under creep stress, and carbide coarsening reduces the creep resistance of ferritic steels. Phase field models that couple elasticity (for directional rafting) with diffusion and thermodynamics (for coarsening kinetics) quantitatively predict these degradation processes, providing the basis for remaining life assessment of turbine blades and high-temperature pressure vessels.

Numerical Implementation and Best Practices

Spatial Discretisation

Phase field equations are most commonly solved on regular Cartesian grids using finite differences, which offer straightforward implementation and high computational efficiency. The gradient energy term requires at least 4–5 grid points across the diffuse interface for numerical accuracy. Finite element methods (MOOSE, PRISMS-PF) enable adaptive mesh refinement (AMR) that concentrates computational resources at the diffuse interface, reducing total computational cost by 10–100x for systems with large matrix regions.

Time Integration

Explicit (forward Euler) time integration is simplest but requires a time step constrained by both the CFL stability condition and the interface relaxation time. Semi-implicit or fully implicit schemes (Crank-Nicolson, backward Euler with iterative solvers) allow time steps 10–100x larger at the cost of solving a linear or nonlinear system at each step. For stiff problems involving very different time scales (fast interface evolution, slow long-range diffusion), operator-splitting approaches separate the fast and slow components.

Parameter Calibration