GitHub

Thermal Phase Field Model

Overview

The thermal phase field model couples solidification with heat transfer, capturing latent heat release and thermal gradients. It is essential for modeling realistic solidification where temperature evolution affects interface kinetics.

Mathematical Formulation

Phase Field Equation

\[\tau \frac{\partial \phi}{\partial t} = W^2 \nabla^2 \phi - g'(\phi) + \lambda h'(\phi) u\]

Heat Equation

\[\frac{\partial u}{\partial t} = \alpha \nabla^2 u + \frac{1}{2} \frac{\partial \phi}{\partial t}\]

Dimensionless Temperature

\[u = \frac{T - T_m}{L / C_p}\]

where:

  • T: Physical temperature [K]
  • T_m: Melting temperature [K]
  • L: Latent heat [J/m³]
  • C_p: Heat capacity [J/(m³·K)]
  • α: Thermal diffusivity [m²/s]

Stefan Number

\[St = \frac{C_p \Delta T}{L}\]

Quick Example

Basic Usage

using PhaseFields

# Create model (nickel-like parameters)
model = ThermalPhaseFieldModel(
    τ = 1e-6,       # Relaxation time [s]
    W = 1e-6,       # Interface width [m]
    λ = 2.0,        # Coupling strength
    α = 1e-5,       # Thermal diffusivity [m²/s]
    L = 2.35e9,     # Latent heat [J/m³]
    Cp = 5.42e6,    # Heat capacity [J/(m³·K)]
    Tm = 1728.0     # Melting point [K]
)

# Convert temperature
T = 1700.0  # Physical temperature [K]
u = dimensionless_temperature(model, T)

# Convert back
T_back = physical_temperature(model, u)

With DifferentialEquations.jl

using PhaseFields
using OrdinaryDiffEq

model = ThermalPhaseFieldModel(
    τ=1e-6, W=1e-6, λ=2.0,
    α=1e-5, L=2.35e9, Cp=5.42e6, Tm=1728.0
)

grid = UniformGrid1D(N=101, L=1e-4)
bc_φ = NeumannBC()
bc_u = NeumannBC()

# Initial conditions
φ0 = [x < grid.L/2 ? 1.0 : 0.0 for x in grid.x]  # Solid seed
u0 = fill(-0.05, grid.N)  # Slight undercooling

# Create and solve
prob = create_thermal_problem(model, grid, bc_φ, bc_u, φ0, u0, (0.0, 1e-3))
sol = solve(prob, QNDF(autodiff=false))

# Extract solution
φ_hist, u_hist = extract_thermal_solution(sol, grid.N)

Stefan Problem

The model can simulate the classical Stefan problem where a planar interface advances into an undercooled melt, releasing latent heat.

# Undercooling determines interface velocity
ΔT = 10.0  # [K]
u0 = dimensionless_temperature(model, model.Tm - ΔT)

API Reference

PhaseFields.ThermalPhaseFieldModel — Type
ThermalPhaseFieldModel

Parameters for thermal phase field model coupling phase evolution with heat transfer.

Fields

  • τ::Float64: Relaxation time [s]
  • W::Float64: Interface width parameter [m]
  • λ::Float64: Thermal coupling strength (dimensionless)
  • α::Float64: Thermal diffusivity [m²/s]
  • L::Float64: Latent heat [J/m³]
  • Cp::Float64: Volumetric heat capacity [J/(m³·K)]
  • Tm::Float64: Melting temperature [K]
  • K::Float64: Latent heat coefficient L/Cp [K] (computed)

Example

model = ThermalPhaseFieldModel(
    τ = 1.0,
    W = 1.0,
    λ = 1.0,
    α = 1e-5,
    L = 2.35e9,
    Cp = 5.42e6,
    Tm = 1728.0
)
source
PhaseFields.dimensionless_temperature — Function
dimensionless_temperature(T, Tm, L, Cp)

Convert physical temperature to dimensionless temperature.

Arguments

  • T: Physical temperature [K]
  • Tm: Melting temperature [K]
  • L: Latent heat [J/m³]
  • Cp: Volumetric heat capacity [J/(m³·K)]

Returns

  • u: Dimensionless temperature u = Cp(T - Tm)/L

Notes

  • u = 0 at melting point
  • u < 0 for undercooling (T < Tm)
  • u > 0 for superheating (T > Tm)
source
PhaseFields.physical_temperature — Function
physical_temperature(u, Tm, L, Cp)

Convert dimensionless temperature back to physical temperature.

Arguments

  • u: Dimensionless temperature
  • Tm: Melting temperature [K]
  • L: Latent heat [J/m³]
  • Cp: Volumetric heat capacity [J/(m³·K)]

Returns

  • T: Physical temperature [K]
source
PhaseFields.stefan_number — Function
stefan_number(ΔT, L, Cp)

Calculate Stefan number.

Arguments

  • ΔT: Temperature difference (undercooling or superheating) [K]
  • L: Latent heat [J/m³]
  • Cp: Volumetric heat capacity [J/(m³·K)]

Returns

  • St: Stefan number St = Cp·ΔT/L (sensible heat / latent heat ratio)

Notes

  • St << 1: Latent heat dominates, slow interface motion
  • St ~ 1: Comparable sensible and latent heat
  • St >> 1: Sensible heat dominates, fast interface motion
source
PhaseFields.thermal_phase_rhs — Function
thermal_phase_rhs(model, φ, ∇²φ, u)

Compute right-hand side of phase field evolution equation with thermal driving.

Arguments

  • model::ThermalPhaseFieldModel: Model parameters
  • φ: Phase field value (0=liquid, 1=solid)
  • ∇²φ: Laplacian of phase field
  • u: Dimensionless temperature

Returns

  • ∂φ/∂t: Time derivative of phase field (divide by τ is already done)

Equation

τ ∂φ/∂t = W²∇²φ - g'(φ) + λu h'(φ)

Returns the RHS divided by τ:

∂φ/∂t = (W²∇²φ - g'(φ) + λu h'(φ)) / τ

Notes

  • g(φ) = φ²(1-φ)² double-well potential, g'(φ) = 2φ(1-φ)(1-2φ)
  • h(φ) = φ²(3-2φ) interpolation function, h'(φ) = 6φ(1-φ)
  • Undercooling (u < 0) drives solidification (increases φ)
  • Superheating (u > 0) drives melting (decreases φ)
source
PhaseFields.thermal_heat_rhs — Function
thermal_heat_rhs(model, u, ∇²u, dφdt)

Compute right-hand side of heat equation with latent heat source.

Arguments

  • model::ThermalPhaseFieldModel: Model parameters
  • u: Dimensionless temperature
  • ∇²u: Laplacian of dimensionless temperature
  • dφdt: Time derivative of phase field

Returns

  • ∂u/∂t: Time derivative of dimensionless temperature

Equation

∂u/∂t = α∇²u + (1/2) ∂φ/∂t

Notes

  • Solidification (dφdt > 0) releases latent heat (u increases)
  • Melting (dφdt < 0) absorbs latent heat (u decreases)
  • Factor 1/2 comes from h(φ) normalization where solid fraction = h(φ)
source
PhaseFields.thermal_stability_dt — Function
thermal_stability_dt(model, dx)

Compute stable time step for coupled thermal phase field evolution.

Arguments

  • model::ThermalPhaseFieldModel: Model parameters
  • dx: Grid spacing [m]

Returns

  • dt: Stable time step [s]

Notes

Combines stability conditions:

  • Thermal diffusion: dt < dx²/(2α) in 1D, dx²/(4α) in 2D
  • Phase field: dt < τ·dx²/(2W²)

Returns conservative (1D) limit with safety factor.

source
PhaseFields.ThermalODEParams — Type
ThermalODEParams{M, G, BC1, BC2, A}

Parameters for coupled thermal phase field ODE integration.

Type Parameters

  • M: Model type
  • G: Grid type (UniformGrid1D or UniformGrid2D)
  • BC1: Phase field boundary condition type
  • BC2: Temperature boundary condition type
  • A: Laplacian workspace array type

Fields

  • model::M: Thermal phase field model
  • grid::G: Spatial discretization (1D or 2D)
  • bc_φ::BC1: Phase field boundary conditions
  • bc_u::BC2: Temperature boundary conditions
  • ∇²φ::A: Pre-allocated Laplacian workspace for φ
  • ∇²u::A: Pre-allocated Laplacian workspace for u
  • dφdt::A: Pre-allocated dφ/dt workspace
source
PhaseFields.thermal_ode! — Function
thermal_ode!(dy, y, p::ThermalODEParams, t)

In-place RHS function for coupled thermal phase field equations.

State vector y = [φ; u] where:

  • φ: Phase field (indices 1:N for 1D, 1:Nx*Ny for 2D)
  • u: Dimensionless temperature (indices N+1:2N for 1D)
source
PhaseFields.create_thermal_problem — Function
create_thermal_problem(model, grid, bc_φ, bc_u, φ0, u0, tspan)

Create ODEProblem for coupled thermal phase field equations.

Arguments

  • model::ThermalPhaseFieldModel: Thermal phase field model
  • grid::UniformGrid1D: Spatial grid
  • bc_φ::BoundaryCondition: Phase field boundary conditions
  • bc_u::BoundaryCondition: Temperature boundary conditions
  • φ0: Initial phase field
  • u0: Initial dimensionless temperature
  • tspan: Time span (tstart, tend)

Returns

  • ODEProblem with state vector [φ; u]
source
PhaseFields.extract_thermal_solution — Function
extract_thermal_solution(sol, N)

Extract φ and u from thermal ODE solution.

Arguments

  • sol: ODESolution from solve()
  • N: Number of grid points

Returns

  • (φ_history, u_history): Arrays of shape (N, n_times)
source

See Also