Enthalpy Method for Phase Change Analysis

Standard heat conduction analysis solves $\rho c_p \partial T/\partial t = \nabla \cdot (k\nabla T) + Q$. This works fine for solid or liquid phases separately, but it fundamentally cannot represent the solid-liquid phase transition. When metal solidifies, latent heat is released at constant temperature — that's heat flowing without any temperature change. The standard equation has no mechanism for this. The enthalpy method handles it by reformulating with total enthalpy $H$ as the primary variable, which naturally absorbs both sensible heat and latent heat. It's essential for any simulation involving melting or solidification: casting, welding, PCM thermal storage, ice formation, and additive manufacturing melt pools.

Classic example: ice melting at 0°C. You put ice at exactly 0°C on a warm surface. Heat flows from the surface into the ice. But the ice stays at 0°C — all the incoming energy goes into breaking the crystal lattice bonds (latent heat of fusion = 334 kJ/kg for water). Only after all the ice has melted does the temperature start rising. If you were running a standard heat conduction code, it would predict the temperature jumping immediately since no $c_p$ spike at exactly 0°C is defined. The enthalpy formulation handles this correctly by tracking $H$ — when $H$ is in the phase change range, $T$ stays constant while $H$ increases.

The enthalpy transport equation using total enthalpy $H$ as the primary variable:

Total enthalpy $H$ separates into sensible and latent contributions:

where $f_l \in [0, 1]$ is the liquid fraction and $L_f$ is the latent heat of fusion (J/kg). In the liquid phase $f_l = 1$, solid phase $f_l = 0$, mushy zone $0 < f_l < 1$. The temperature is recovered from $H$ through the inverse relation $T(H)$, which requires nonlinear iteration when using $H$ as the marching variable.

An alternative formulation absorbs latent heat into an effective specific heat:

The second term is the "heat capacity spike" at the phase transition. Both formulations are equivalent mathematically; they differ in implementation and numerical behavior.

Pure substances (pure aluminum, pure water) melt at a single temperature — a sharp phase transition. Alloys and polymers have a two-phase region between solidus temperature $T_s$ (where solidification begins on cooling) and liquidus temperature $T_l$ (where it is complete). This is the mushy zone, so named because the material has a semi-solid, paste-like consistency in this range. The liquid fraction varies with temperature:

This is the linear lever rule approximation — more accurate Scheil or Clyne-Kurz models account for the nonlinear partitioning of solute between phases during solidification. The wider the mushy zone temperature range, the more numerically tractable the phase change is. For alloys with $T_l - T_s > 20$°C (common in aluminum casting alloys), the enthalpy method converges readily without special treatment.

The classical Stefan problem is a special case: a pure substance with a sharp melting point, where the solid-liquid interface is a mathematically sharp boundary that moves over time. The interface position $s(t)$ must be tracked explicitly, satisfying the Stefan condition at the interface:

This says: the difference in heat flux across the interface equals the latent heat released per unit area per unit time by the moving boundary. The Stefan problem is difficult to solve numerically because you must track the moving interface. The enthalpy method's breakthrough was to eliminate explicit interface tracking entirely — the interface location is implicit in the $f_l$ field and emerges naturally from the solution without any special treatment. For a pure substance, the mushy zone collapses to a sharp interface in the limit, and the enthalpy method recovers the correct Stefan condition automatically.

Here's a practical comparison:

In practice, the latent heat source term method is the most widely used in commercial codes (Fluent, OpenFOAM) and gives the best combination of accuracy and robustness. The $H$-based formulation is preferred in dedicated casting software (ProCAST) where it is tightly integrated with phase diagram data.

The key stabilization strategies from practice:

• Temperature smearing (mushy zone widening): Widen the mushy zone by ±1–3°C beyond the measured solidus-liquidus range. This reduces the height of the $\partial f_l/\partial T$ spike and significantly improves convergence for alloys with narrow mushy zones.
• Under-relaxation on liquid fraction: Apply a relaxation factor $\omega = 0.7$–0.9 to $f_l$ updates: $f_l^{new} = f_l^{old} + \omega(f_l^{computed} - f_l^{old})$. Prevents the overshooting that causes oscillations.
• Adaptive time stepping: Automatically reduce $\Delta t$ in cells where $|\partial f_l/\partial t|$ is large (i.e., cells currently undergoing phase change). Most commercial codes support this automatically.
• Isothermal banding: For pure metals (zero mushy zone width), artificially define a narrow band of 0.5–1°C to prevent the delta-function singularity. Values smaller than this are not physically meaningful given measurement uncertainties.
• Iterative $f_l$ update (fixed-point iteration): At each time step, iterate the $T$–$f_l$ coupling until $|f_l^{n+1} - f_l^n| < 10^{-4}$ before advancing to the next time step. This is the standard approach in source-based methods.

The mold-casting interface heat transfer coefficient (HTC) is the single most important and uncertain parameter in casting solidification analysis. Physical stages:

1. Liquid metal contact (0–few seconds): HTC is highest — 2,000–10,000 W/(m²K) for permanent steel molds, reflecting direct liquid metal contact.
2. Solidification and air gap formation: As the casting solidifies and contracts, a microscopic air gap forms between the casting surface and the mold wall. HTC drops dramatically to 200–500 W/(m²K) and continues decreasing as the gap grows.
3. Full solidification: Radiative heat transfer across the air gap dominates. HTC continues to decrease.

In simulations, implement a time-varying or temperature-varying HTC via a user-defined function. Using a constant HTC that ignores air gap formation will produce 20–40% errors in solidification front position and defect prediction. For sand casting, HTC is 250–750 W/(m²K) and more uniform throughout.

Here's the complete step-by-step workflow for aluminum casting solidification:

1. Material data collection: For the alloy (e.g., A356 aluminum), you need: $T_s = 555$°C, $T_l = 615$°C, $L_f = 397$ kJ/kg, temperature-dependent $k(T)$ and $c_p(T)$ for both solid and liquid phases. For the mold (steel or sand), standard thermal properties.
2. Geometry and mesh: Include both casting and mold in a single coupled thermal model. Cell size in the mushy zone region: 1–3 mm for typical die-cast components. Use finer cells near the gate and runner where solidification begins and ends last.
3. Initial conditions: Casting initialized at the pouring temperature (e.g., 700°C). Mold initialized at the preheated mold temperature (150–250°C for die casting).
4. Boundary conditions: Convective cooling on the external mold surface (h = 20–50 W/(m²K) for air cooling; 500–5000 W/(m²K) for water-cooled inserts). Time-varying HTC at the mold-casting interface.
5. Time stepping: Start with $\Delta t = 0.05$–0.1 s. Enable adaptive time stepping — the solver will automatically reduce $\Delta t$ near the solidification front.
6. Post-processing: Generate solidification sequence maps (time for each cell to solidify). Identify hot spots — regions that solidify last are prone to shrinkage porosity. Feeding paths should be maintained from the riser to all hot spots.
7. Validation: Compare predicted total solidification time against measurement. Compare thermal history at thermocouple positions if data is available. Shrinkage porosity predictions should correlate with X-ray or sectioning inspection of trial castings.

Almost certainly numerical, not physical. The real melting front is smooth. What you're seeing is the effective specific heat spike causing overshoot-undershoot oscillations in the temperature field. Definitive diagnosis: if the oscillations are periodic with a spatial wavelength of 2–3 cells and appear only near the melting temperature, they're numerical. Fix sequence to try in order: (1) Widen the mushy zone artificially (add ±2°C to your $T_s$ and $T_l$); (2) Reduce the under-relaxation factor for $f_l$ to 0.6–0.7; (3) Reduce the time step by 2×; (4) Switch from effective specific heat to the latent heat source formulation. Most cases resolve with steps (1) and (2) alone.

Here's a comparison of the main tools:

Several active frontiers are transforming how phase change is simulated:

• Additive Manufacturing Melt Pool Simulation: Laser Powder Bed Fusion (LPBF) and Directed Energy Deposition (DED) processes involve rapid, repeated melting and solidification at the microscale. Multi-physics simulations coupling enthalpy method heat transfer with fluid flow (Marangoni convection from surface tension gradients), vaporization, and powder particle dynamics are the state of the art. Enthalpy-LBM (Lattice Boltzmann) methods are particularly active for melt pool microstructure prediction.
• Phase-Field Method Integration: Phase-field models represent the solid-liquid interface as a diffuse layer with a scalar order parameter. They naturally handle dendrite growth, eutectic solidification, and grain structure evolution. Coupling phase-field with enthalpy method provides both macroscale temperature accuracy and microscale microstructure information — critical for predicting mechanical properties of castings.
• CALPHAD-Coupled Enthalpy: CALPHAD (CALculation of PHAse Diagrams) thermodynamic databases provide $T_s$, $T_l$, $L_f$, and $f_l(T)$ directly from alloy composition for multi-component systems. ProCAST and JMatPro integrate CALPHAD with the enthalpy method for accurate multi-component alloy solidification.
• PCM Optimization for Thermal Storage: Phase Change Materials (PCMs) in building thermal mass, battery thermal management, and concentrated solar power require enthalpy method analysis at the system scale. ML-assisted PCM formulation optimization (matching melting point and latent heat to application) is very active.

Here are the most common problems and their remedies: