AMATH 858: Free Boundary Problems and Phase Transitions
Estimated study time: 2 hr 25 min
Table of contents
These notes synthesize material from A. Friedman’s Variational Principles and Free-Boundary Problems, S.H. Davis’s Theory of Solidification, J. Crank’s Free and Moving Boundary Problems, A. Visintin’s Models of Phase Transitions, G.P. Galdi and R. Rannacher’s Fundamental Trends in Fluid-Structure Interaction, and supplementary material from Cambridge DAMTP Solidification and Phase Change notes (Worster), NYU Courant free boundary problem lecture notes, Oxford applied mathematics materials, and MIT OCW materials on phase transitions.
Chapter 1: The Stefan Problem
The classical Stefan problem is the prototypical free boundary problem: find the temperature field in a medium undergoing a phase change, where the boundary between liquid and solid is itself unknown and must be determined as part of the solution. Named after Josef Stefan, who in 1889 studied the melting of the polar ice cap, the problem in fact dates back to Lamé and Clapeyron (1831), who considered the solidification of a liquid in a half-space. The Stefan problem illustrates why free boundary problems are fundamentally harder than standard boundary value problems — the domain of the PDE is coupled to its solution, creating an inherent nonlinearity even when the governing equation in each phase is linear.
1.1 Solidification and Melting as a Free Boundary Problem
Consider a pure substance occupying a region \(\Omega \subset \mathbb{R}^n\) that can exist in two phases: solid and liquid. Let \(\Omega_s(t)\) and \(\Omega_l(t)\) denote the solid and liquid regions at time \(t\), separated by an interface \(\Gamma(t)\). In each phase the temperature \(T(x,t)\) satisfies the heat equation, but the interface position is unknown.
with the conditions on the interface \(\Gamma(t)\):
\[ T = T_m \quad \text{on } \Gamma(t), \]\[ L \rho v_n = k_s \frac{\partial T}{\partial n}\bigg|_s - k_l \frac{\partial T}{\partial n}\bigg|_l \quad \text{on } \Gamma(t), \]where \(\rho_i, c_i, k_i\) are the density, specific heat, and thermal conductivity in phase \(i\), \(T_m\) is the melting temperature, \(L\) is the latent heat per unit mass, and \(v_n\) is the normal velocity of the interface.
The first interface condition states that the temperature equals the equilibrium melting temperature on the phase boundary. The second condition — the Stefan condition — expresses conservation of energy at the interface: the jump in heat flux across the boundary must account for the latent heat released or absorbed as the interface moves.
In the simplest setting, we assume equal and constant material properties in both phases (\(\rho_s = \rho_l = \rho\), \(c_s = c_l = c\), \(k_s = k_l = k\)), so the thermal diffusivity \(\kappa = k/(\rho c)\) is uniform. The problem then reduces to the heat equation in each phase with the Stefan condition coupling them through the interface.
The one-phase Stefan problem arises naturally in many applications: the freezing of a lake (where the ice is essentially at \(0^\circ\)C throughout), or the solidification of a casting where the solid is maintained at the melting temperature by the latent heat release. The mathematical distinction between the one-phase and two-phase problems is significant: in the one-phase problem, the free boundary condition involves the gradient from only one side, simplifying the analysis considerably. The one-phase problem also arises in oxygen diffusion in biological tissue (the “oxygen consumption” problem), where oxygen diffuses through tissue and is consumed at a rate that depends on its concentration, with a free boundary separating oxygenated from anoxic regions.
1.2 The Stefan Condition and Energy Balance
The Stefan condition deserves careful derivation from first principles. Consider a pillbox control volume straddling the interface \(\Gamma(t)\), of thickness \(2\varepsilon\) and cross-sectional area \(A\). The total energy in the pillbox is
\[ E = \int_V \rho\, e(T)\, dV, \]where \(e(T)\) is the internal energy per unit mass. As the interface sweeps through the pillbox, the discontinuity in internal energy (due to the latent heat \(L\)) must be balanced by the net heat flux through the pillbox faces.
where \(v_n\) is the normal velocity of \(\Gamma(t)\) in the direction of \(\hat{n}\).
where \(\mathbf{q} = -k\nabla T\) is the heat flux. Using the generalized divergence theorem for a domain with a moving discontinuity, the left-hand side produces a surface integral over \(\Gamma(t)\) involving the jump in \(\rho e\) multiplied by \(v_n\). Shrinking the pillbox (\(\varepsilon \to 0\)) eliminates volume contributions, leaving
\[ \rho L\, v_n = k_s \frac{\partial T_s}{\partial n} - k_l \frac{\partial T_l}{\partial n} \]on \(\Gamma(t)\). \(\square\)
This derivation highlights a key structural feature: the Stefan condition is not an additional assumption but a consequence of energy conservation applied to a system with a latent heat discontinuity. The same reasoning applies to any conservation law with a moving discontinuity — the Rankine-Hugoniot conditions in gas dynamics are the exact analogue.
1.3 Exact Solutions: The Neumann Solution
The most celebrated exact solution of the Stefan problem is the Neumann solution (1860) for a half-space. Consider a semi-infinite liquid at uniform temperature \(T_0 > T_m\), with the boundary \(x = 0\) suddenly cooled to \(T_w < T_m\) at \(t = 0\). A solid phase grows from \(x = 0\) into the liquid.
The Stefan condition yields the transcendental equation for \(\lambda\):
\[ \frac{k_s(T_m - T_w)}{\sqrt{\pi\kappa_s}\,\operatorname{erf}(\lambda)} e^{-\lambda^2} - \frac{k_l(T_0 - T_m)}{\sqrt{\pi\kappa_l}\,\operatorname{erfc}(\lambda\sqrt{\kappa_s/\kappa_l})} e^{-\lambda^2 \kappa_s/\kappa_l} = \rho L \lambda \sqrt{\kappa_s}. \]This equation has a unique positive root \(\lambda\) for any physically meaningful parameter values.
The key insight from the Neumann solution is the \(\sqrt{t}\) growth law: the solidification front advances as \(s(t) \sim \sqrt{t}\). This parabolic scaling is universal for diffusion-controlled free boundaries and persists even when the exact solution is unavailable. Experimentally, measuring the departure from \(\sqrt{t}\) scaling is a diagnostic for convective effects or other non-diffusive processes.
1.4 Existence, Uniqueness, and Regularity
The mathematical theory of the Stefan problem was placed on rigorous footing by Friedman, Kamenomostskaya (Kamin), and others in the 1960s. The classical formulation requires smooth interfaces and classical solutions of the heat equation in each phase, but this breaks down when interfaces develop singularities (e.g., topological changes such as merging of two solidification fronts).
In higher dimensions, the classical theory becomes more delicate. Interfaces can develop corners, cusps, and topological changes. This motivates the weak (enthalpy) formulation.
1.5 The Enthalpy Method and Weak Formulation
To handle irregular interfaces and topological changes, we reformulate the Stefan problem in terms of the enthalpy \(H(T)\), which incorporates the latent heat as a discontinuity in the energy-temperature relationship.
where \(\xi\) represents the liquid fraction. The Stefan problem is equivalent to finding \(H\) and \(T\) satisfying:
\[ \frac{\partial H}{\partial t} = \nabla \cdot (k(T)\nabla T) \quad \text{in } \Omega \times (0,\infty) \]in the distributional sense, with \(T = T(H)\) given by the inverse of the enthalpy relation.
The enthalpy formulation is a single equation valid throughout \(\Omega\), with no explicit reference to the interface. The free boundary is recovered as the level set \(\{x : H(x,t) \in [\rho c_s T_m, \rho c_s T_m + \rho L]\}\). This weak formulation permits solutions with irregular or mushy interfaces and is the basis for most numerical methods.
The enthalpy method is also the starting point for numerical approximation. Discretizing the enthalpy equation directly avoids the need to track the interface explicitly, and the latent heat appears as a source term concentrated near \(T = T_m\). This approach, pioneered by Voller and others in the 1980s, remains the workhorse of computational phase-change modeling.
where \(T_j = T(H_j)\) is obtained by inverting the enthalpy relation. Since \(H(T)\) is multivalued at \(T = T_m\), this nonlinear system is solved by iteration (e.g., the Voller-Swaminathan source-based method, which linearizes the latent heat source). For \(N = 100\) and \(\Delta t = 10^{-4}\), the computed interface position agrees with the exact (Neumann) solution to within \(O(\Delta x)\), demonstrating the first-order accuracy of the enthalpy method at the interface — a consequence of smearing the latent heat over one cell width.
1.6 The Stefan Number and Dimensional Analysis
The key dimensionless parameter governing the Stefan problem is the Stefan number, which measures the ratio of sensible heat to latent heat.
where \(\Delta T\) is the characteristic temperature difference (e.g., \(T_m - T_w\) for solidification against a cold wall) and \(L\) is the latent heat per unit mass. When \(\mathrm{St} \ll 1\), latent heat dominates and the interface moves slowly; when \(\mathrm{St} \gg 1\), sensible heat dominates and the problem approaches the pure conduction limit without phase change.
For water at atmospheric pressure, \(L \approx 334\) kJ/kg and \(c \approx 4.2\) kJ/(kg K), so a 10 K undercooling gives \(\mathrm{St} \approx 0.13\). For metals, \(\mathrm{St}\) is typically in the range 0.1–1. The small Stefan number limit (\(\mathrm{St} \to 0\)) is a singular perturbation in which the temperature field adjusts quasi-statically to the slowly moving boundary — this is the quasi-steady approximation widely used in engineering.
The Stefan condition then gives:
\[ \rho L \dot{s} = k\frac{T_m - T_w}{s}, \qquad \Rightarrow \qquad s(t) = \sqrt{\frac{2k(T_m - T_w)t}{\rho L}} = \sqrt{2\,\mathrm{St}\,\kappa\, t}. \]This recovers the \(\sqrt{t}\) growth law with an explicit prefactor proportional to \(\sqrt{\mathrm{St}}\).
The quasi-steady approximation is the foundation of many engineering correlations for solidification and melting. Its accuracy improves as \(\mathrm{St} \to 0\), and correction terms can be computed systematically via a perturbation expansion in \(\mathrm{St}\).
so that \(s(t) = \sqrt{2\,\mathrm{St}\,\kappa\, t}\left(1 + \mathrm{St}/(2\pi) + \cdots\right)\). The first correction increases the interface speed by a factor \(1 + \mathrm{St}/(2\pi) \approx 1.02\) for \(\mathrm{St} = 0.13\) (water), confirming that the quasi-steady approximation is excellent for small Stefan numbers.
1.7 The Multi-Dimensional Stefan Problem and Comparison Principles
In multiple space dimensions, the Stefan problem has a much richer structure. The interface is now a surface in \(\mathbb{R}^n\) that can develop complex geometries, and the interplay between curvature and the thermal field introduces new phenomena.
The comparison principle is a fundamental tool for obtaining a priori bounds on the free boundary location. It implies, for example, that the solidification front in a bounded domain always lies between the fronts computed with maximum and minimum boundary temperatures, providing rigorous upper and lower bounds on the interface position.
where \(\lambda\) is the Neumann root and \(\mathrm{St} = c|T_w|/L\). For a 10 cm slab of water with \(T_w = -10^\circ\)C (\(\mathrm{St} \approx 0.13\)), these bounds give \(t_f\) between roughly 27 and 36 hours, consistent with the everyday experience of freezing a thick layer of water.
Chapter 2: Hele-Shaw Flows
The Hele-Shaw cell — two closely spaced parallel plates with a viscous fluid between them — provides one of the most elegant physical realizations of a free boundary problem. The mathematical structure is remarkably clean: the pressure satisfies Laplace’s equation in the fluid domain, and the interface moves with velocity proportional to the pressure gradient. This Laplacian growth model connects solidification theory to potential theory, conformal mapping, and integrable systems.
2.1 Darcy’s Law and the Hele-Shaw Cell
Henry Selby Hele-Shaw introduced his cell in 1898 as an analogue device for visualizing potential flow around obstacles. The gap-averaged velocity in a Hele-Shaw cell satisfies Darcy’s law, the same equation governing flow through a porous medium.
so that \(p\) is harmonic in the fluid domain:
\[ \nabla^2 p = 0 \quad \text{in } \Omega(t). \]On the free boundary \(\Gamma(t)\) separating fluid from air (or a less viscous fluid):
\[ p = -\sigma \kappa \quad \text{on } \Gamma(t), \qquad v_n = -\frac{b^2}{12\mu}\frac{\partial p}{\partial n} \quad \text{on } \Gamma(t), \]where \(\sigma\) is the surface tension, \(\kappa\) is the interface curvature, and \(v_n\) is the normal velocity of \(\Gamma(t)\).
When surface tension is neglected (\(\sigma = 0\)), the pressure is simply \(p = 0\) on the free boundary, and the problem reduces to finding a domain \(\Omega(t)\) in which a harmonic function vanishes on the boundary and has prescribed flux there. This is the zero-surface-tension Hele-Shaw problem, which is mathematically identical to the one-phase Stefan problem with \(T_m = 0\) after the identification \(p \leftrightarrow T\), \(v_n \leftrightarrow v_n\).
This correspondence between Hele-Shaw flow and the Stefan problem is more than a formal analogy. It means that any analytical or numerical technique developed for one problem can be transferred to the other. For example, the comparison principles and weak formulations of the Stefan problem (Chapter 1) apply directly to Hele-Shaw flow, while the conformal mapping techniques and conserved moments that are natural for Hele-Shaw flow can be used to construct exact solutions of the Stefan problem. The correspondence also extends to the instability theory: the Saffman-Taylor instability (Section 2.2) and the Mullins-Sekerka instability (Chapter 3) are the same mathematical phenomenon expressed in different physical variables.
2.2 Saffman-Taylor Instability
In 1958, Saffman and Taylor demonstrated experimentally and theoretically that when a less viscous fluid (e.g., air) displaces a more viscous fluid in a Hele-Shaw cell, the planar interface is unstable. This is the viscous fingering instability, one of the most visually dramatic pattern-forming instabilities in fluid mechanics.
grows at rate
\[ \sigma = V_0 k - \frac{b^2 \sigma_{\text{st}}}{12\mu} k^3, \]where \(V_0\) is the unperturbed interface speed and \(\sigma_{\text{st}}\) is the surface tension. All wavelengths with \(k < k_c = \sqrt{12\mu V_0/(b^2 \sigma_{\text{st}})}\) are unstable, with the most dangerous wavenumber at \(k_c/\sqrt{3}\).
Since \(\sigma > 0\) for small \(|k|\) and \(\sigma < 0\) for large \(|k|\), there is a band of unstable wavenumbers. \(\square\)
The Saffman-Taylor instability is the Hele-Shaw analogue of the Mullins-Sekerka instability in solidification (Chapter 3). Both arise from the same mathematical mechanism: a diffusive field (temperature or pressure) satisfying Laplace’s equation, with growth favoring protrusions because they experience enhanced gradients, stabilized only by surface energy at short wavelengths.
where \(\eta\) parameterizes the finger boundary. Experiments, however, robustly select \(\lambda \approx 1/2\) (the finger occupies half the channel width). This selection puzzle, resolved by the inclusion of surface tension as a singular perturbation (McLean and Saffman, 1981; Combescot et al., 1986; Shraiman, 1986), is the Hele-Shaw analogue of the dendrite tip selection problem (Chapter 3). In both cases, surface tension — no matter how small — is a singular perturbation that selects a unique solution from a continuum.
2.3 Richardson’s Moments and Exact Solutions
A remarkable feature of the zero-surface-tension Hele-Shaw problem is the existence of an infinite number of conserved quantities, discovered by Richardson (1972). These moments constrain the evolution of the fluid domain and connect Hele-Shaw flow to integrable systems.
(where \(z = x + iy\) is the complex coordinate) satisfy
\[ \dot{M}_0 = Q, \qquad \dot{M}_k = 0, \quad k \geq 1. \]That is, all moments except the area are conserved.
Richardson’s moments provide a powerful tool for constructing exact solutions via conformal mapping. If the initial domain is described by a conformal map \(z = f(\zeta, t)\) from the unit disk, the Polubarinova-Galin equation governs the evolution of this map. The conservation of the moments \(M_k\) for \(k \geq 1\) severely constrains the possible domain evolutions. For example, the centroid \(M_1/M_0\) and the moments of inertia \(M_2/M_0\) are determined entirely by their initial values and the area growth rate. This means that an initially circular domain driven by off-center injection will not remain circular, but its centroid will remain fixed and its aspect ratio will evolve in a prescribed manner.
The Schwarz function approach, developed by Davis, Crowdy, and others, provides an alternative and often more powerful framework. The Schwarz function \(S(z)\) of a curve \(\Gamma\) satisfies \(\bar{z} = S(z)\) on \(\Gamma\), and the Hele-Shaw evolution translates into an elegant evolution equation for \(S(z)\).
where \(z_0\) is the injection point and \(Q\) the injection rate. The singularities of \(S\) (which lie inside the fluid domain) move according to a Hamiltonian dynamics that is equivalent to the dispersionless limit of the Toda lattice hierarchy. This connection to integrable systems, discovered by Mineev-Weinstein, Wiegmann, and Zabrodin (2000), reveals that Hele-Shaw flow possesses infinitely many conserved quantities (Richardson’s moments) and a tau-function structure familiar from soliton theory.
The domain remains circular for all time, and the pressure is:
\[ p(r,t) = \frac{Q\mu}{2\pi b^2/12}\ln\frac{r}{R(t)}, \quad r < R(t). \]This exact solution provides a benchmark for numerical methods and a base state for stability analysis.
2.4 Connection to DLA and Laplacian Growth
In the zero-surface-tension limit, Hele-Shaw flow is a deterministic version of diffusion-limited aggregation (DLA), the stochastic growth model introduced by Witten and Sander (1981). In DLA, random walkers are released from infinity and stick upon contact with a growing cluster. The probability of attachment at a point on the boundary is proportional to the harmonic measure — the same quantity that controls the growth rate in Hele-Shaw flow. Both models produce ramified, fractal-like structures. The connection between Hele-Shaw flow, DLA, and the broader class of Laplacian growth models remains an active area of research, with deep connections to random matrix theory, integrable systems, and conformal field theory.
Chapter 3: Morphological Instability
When a solidification front is not perfectly planar, small perturbations can grow or decay depending on the competition between the destabilizing effect of the thermal (or solutal) gradient and the stabilizing effect of surface energy. The linear stability analysis of a planar front, carried out by Mullins and Sekerka in their landmark 1964 paper, provides the theoretical foundation for understanding dendritic growth, cellular solidification, and the rich morphological zoo observed in experiments.
3.1 Linear Stability of a Planar Front
Consider a planar solidification front advancing at steady velocity \(V\) into a supercooled melt. In a frame moving with the front, the temperature in the liquid satisfies:
\[ \kappa \nabla^2 T + V \frac{\partial T}{\partial z} = 0, \]where \(z\) is the coordinate normal to the front. The base state solution is \(T_0(z) = T_\infty + \Delta T\, e^{-Vz/\kappa}\) for \(z > 0\) (liquid), where \(\Delta T = T_m - T_\infty\) is the undercooling.
where, in the symmetric model (\(k_s = k_l\), \(\kappa_s = \kappa_l\)), this simplifies to
\[ \sigma(k) = \frac{V|k|}{2}\left(1 - \frac{2\kappa \Gamma_{\text{GT}}}{V} k^2 - \frac{2\kappa}{V^2}|k|G\right), \]with \(G\) the (stabilizing) temperature gradient in the liquid ahead of the front and \(\Gamma_{\text{GT}} = \sigma_{sl} T_m / (\rho L)\) the Gibbs-Thomson capillary length. There exists a critical wavenumber
\[ k_c = \sqrt{\frac{V}{2\kappa \Gamma_{\text{GT}}}} \]above which perturbations are stable, and a most dangerous wavenumber \(k_{\max} \approx k_c/\sqrt{3}\).
The physical mechanism is clear: a bump on the solidification front protrudes into the supercooled liquid, where it experiences a larger temperature gradient and therefore solidifies faster, amplifying the perturbation. Surface energy opposes this by penalizing high curvature (the Gibbs-Thomson effect), creating a short-wavelength cutoff. The result is a band of unstable wavelengths, with the most dangerous wavelength setting the initial scale of the emerging pattern.
where \(G_C\) is the solute concentration gradient at the interface. When the temperature gradient is insufficient to keep the liquid ahead of the interface above the liquidus temperature, constitutional supercooling occurs and the planar front is unstable.
This is consistent with the observed scale of initial dendritic instabilities in ice.
3.2 The Gibbs-Thomson Condition
The key stabilizing mechanism is the dependence of the equilibrium melting temperature on interface curvature, known as the Gibbs-Thomson condition. This thermodynamic relation accounts for the energetic cost of creating curved interfaces.
where \(\Gamma_{\text{GT}} = \sigma_{sl} T_m / (\rho L)\) is the capillary length (typically \(\sim 10^{-7}\) m for metals), \(\sigma_{sl}\) is the solid-liquid interfacial energy, and \(\mu_k\) is the interface kinetic coefficient. For slowly moving interfaces, the kinetic term is often neglected.
The capillary length \(\Gamma_{\text{GT}}\) is tiny — on the order of angstroms for most metals — yet it plays a crucial role in pattern selection. Without surface energy, the Mullins-Sekerka instability would predict unbounded growth rates at short wavelengths, leading to an ill-posed problem. The Gibbs-Thomson condition regularizes the problem and selects the characteristic scale of dendritic sidebranches.
If the interface advances by a distance \(\delta n\) normal to itself, the volume of solid increases by \(\delta V = A\,\delta n\) and the area changes by \(\delta A = -\kappa A\,\delta n\) (by the first variation of area formula). At temperature \(T\) near \(T_m\), the bulk free energy change per unit volume upon solidification is \(\Delta g = \rho L(T_m - T)/T_m\). At equilibrium (\(\delta G = 0\)):
\[ 0 = \Delta g\, A\,\delta n + \sigma_{sl}(-\kappa A\,\delta n) = \left(\rho L\frac{T_m - T}{T_m} - \sigma_{sl}\kappa\right)A\,\delta n. \]Solving for the equilibrium temperature gives \(T = T_m - \sigma_{sl}T_m\kappa/(\rho L) = T_m - \Gamma_{\text{GT}}\kappa\), which is the Gibbs-Thomson condition without kinetic effects. The kinetic undercooling \(V/\mu_k\) arises from the non-equilibrium thermodynamics of the attachment process at the interface, modeled as a thermally activated rate process. \(\square\)
3.3 Dendritic Growth and the Ivantsov Solution
Dendrites — the tree-like crystal structures that form during solidification — are perhaps the most iconic morphology in materials science. The mathematical theory of dendrite tip selection is one of the great achievements of twentieth-century applied mathematics.
where \(\operatorname{Pe} = VR/(2\kappa)\) is the Péclet number based on the tip radius \(R\), \(E_1\) is the exponential integral, and \(\mathrm{St}\) is the Stefan number. The Ivantsov relation determines only the product \(VR\), not \(V\) and \(R\) separately.
The Ivantsov solution reveals a fundamental degeneracy: for a given undercooling, there is a one-parameter family of paraboloidal solutions parameterized by \(VR = \text{const}\). Fast, sharp dendrites and slow, blunt dendrites are equally valid solutions of the transport problem alone. An additional selection mechanism is needed, and it is provided by crystalline anisotropy.
3.4 Tip Selection and Microscopic Solvability
The tip selection problem — determining which member of the Ivantsov family is actually observed — was resolved in the 1980s by the theory of microscopic solvability (Kessler, Koplik, Levine; Barbieri, Hong, Langer; Ben Amar, Pomeau).
where \(d_0 = \Gamma_{\text{GT}}c_l/L\) is the chemical capillary length and \(\sigma^*\) depends on the anisotropy strength.
The essential insight is that anisotropy, no matter how small, breaks the degeneracy of the Ivantsov family. Without anisotropy (\(\epsilon_4 = 0\)), there is no steady needle crystal solution — the isotropic problem is structurally unstable. This is a singular perturbation: the anisotropy enters at the tip, where the curvature is largest, and selects a unique operating point. The theory has been confirmed by detailed numerical simulations and by the beautiful microgravity experiments of Glicksman and collaborators.
3.5 Dendritic Sidebranching and Pattern Formation
Beyond the steady-state tip, real dendrites exhibit sidebranches — secondary arms that grow perpendicular to the primary trunk. Sidebranches are generated by the amplification of noise (thermal fluctuations) at the tip through the Mullins-Sekerka mechanism as perturbations are advected away from the tip along the dendrite sides.
where \(\sigma(k, s')\) is the local Mullins-Sekerka growth rate evaluated with the local curvature and velocity of the dendrite surface. The most amplified perturbations have wavelengths comparable to the tip radius, consistent with experimental observations.
The competition between deterministic tip selection and stochastic sidebranching produces the characteristic “Christmas tree” morphology of dendrites. At large distances from the tip, sidebranches undergo competitive growth (coarsening), with some arms growing at the expense of their neighbors. The resulting microstructure — the spacing, orientation, and distribution of dendrite arms — controls the mechanical and thermal properties of the solidified material.
The sidebranching phenomenon also provides a natural amplifier for noise. Brener and Temkin (1995) showed that the selective amplification mechanism means that sidebranches originating from thermal noise have a universal spacing (independent of the noise amplitude, depending only logarithmically on it), explaining why the sidebranch spacing in experiments is remarkably reproducible even though the underlying fluctuations are stochastic. This insensitivity to the noise amplitude is a hallmark of the convective amplification mechanism: perturbations grow exponentially as they travel down the dendrite, so the final amplitude is dominated by the growth factor rather than the initial seed.
Chapter 4: Mushy Zones and Alloy Solidification
Real materials are rarely pure substances. Most engineering alloys, geological melts, and biological fluids are multi-component systems in which solidification produces a two-phase mushy zone — a porous matrix of solid dendrites bathed in interdendritic liquid. The mathematical description of mushy zones, pioneered by Worster in the late 1980s and 1990s at Cambridge, combines ideas from phase-change thermodynamics, porous media flow, and convective instability.
4.1 Binary Alloy Phase Diagrams
The solidification of a binary alloy (two components, A and B) is governed by the equilibrium phase diagram, which specifies the phases present at each temperature and composition.
where \(m_L > 0\) is the (absolute) liquidus slope, \(k_p < 1\) is the equilibrium partition coefficient (the ratio of solid to liquid compositions at equilibrium), and \(T_m\) is the melting temperature of pure A. Between the liquidus and solidus, solid and liquid coexist in the mushy zone.
When a binary alloy at initial composition \(C_0\) is cooled from above the liquidus, solidification begins at the liquidus temperature \(T_L(C_0)\). The first solid to form has composition \(k_p C_0 < C_0\), rejecting solute into the liquid. This solute rejection drives constitutional supercooling, which is the dominant mechanism for morphological instability in alloy solidification.
where \(f_s\) is the fraction solidified. The Scheil equation predicts that the last liquid to solidify has very high tin content, approaching the eutectic composition, which explains the observation of eutectic pockets between primary dendrite arms.
4.2 The Mushy Layer as a Porous Medium
When constitutional supercooling is sufficiently strong, the planar solidification front breaks down and a mushy layer forms. This layer consists of a fine dendritic network with interstitial liquid, and on length scales much larger than the dendrite spacing, it can be treated as a porous medium.
Conservation of heat and solute yield:
\[ \rho c \frac{\partial T}{\partial t} + \rho L \frac{\partial \phi}{\partial t} = \nabla \cdot (k_{\text{eff}} \nabla T), \]\[ (1 - \phi)\frac{\partial C}{\partial t} + C(1 - k_p)\frac{\partial \phi}{\partial t} = \nabla \cdot (D_{\text{eff}}(1-\phi)\nabla C) + (1-\phi)\mathbf{u} \cdot \nabla C, \]where the interstitial liquid velocity \(\mathbf{u}\) satisfies Darcy’s law:
\[ \mathbf{u} = -\frac{\Pi(\phi)}{\mu}(\nabla p + \rho_l \beta_C C\, \mathbf{g}), \]with permeability \(\Pi(\phi)\) (e.g., the Kozeny-Carman relation \(\Pi = \Pi_0 (1-\phi)^3/\phi^2\)).
The mushy-layer equations couple heat and solute transport to the solid fraction through the liquidus constraint. The resulting system is a degenerate parabolic system — degenerate because the solid fraction \(\phi\) appears as a coefficient that can vanish (at the mush-liquid boundary) or reach unity (at the mush-solid boundary).
The local thermodynamic equilibrium assumption \(T = T_L(C)\) is the key simplification that makes the mushy-layer equations tractable. It states that at every point within the mush, the temperature and interstitial liquid composition lie on the liquidus curve — the solid and liquid phases are in local chemical equilibrium. This assumption is valid when the diffusion time across a dendrite arm spacing \(d_2\) is much shorter than the macroscopic evolution time: \(d_2^2/D_l \ll h/V\), where \(h\) is the mush thickness and \(V\) the growth velocity. For metallic alloys with \(D_l \sim 10^{-9}\) m\(^2\)/s, \(d_2 \sim 10^{-4}\) m, \(h \sim 10^{-2}\) m, and \(V \sim 10^{-4}\) m/s, the ratio is \(d_2^2 V/(D_l h) \sim 10^{-3} \ll 1\), confirming the validity of the equilibrium assumption. For rapid solidification or very coarse microstructures, departures from local equilibrium must be modeled explicitly.
4.3 Chimney Convection
One of the most striking phenomena in mushy-layer solidification is the formation of chimneys — narrow, liquid-filled channels that penetrate the mushy layer and transport solute-rich liquid upward (or downward, depending on the buoyancy contrast). Chimneys are responsible for freckle defects in nickel-base superalloy castings and for brine channels in sea ice.
where \(\Delta C\) is the characteristic composition difference across the mushy layer, \(h\) is its thickness, and \(\Pi_0\) is a reference permeability. Convective instability sets in when \(\mathrm{Ra}_m\) exceeds a critical value \(\mathrm{Ra}_{m,c}\) that depends on the concentration ratio \(\mathcal{C} = C_0(1 - k_p)/(C_0 - C_e)\), the Stefan number, and the Darcy number.
4.4 Homogenization and Effective Properties
The derivation of the mushy-layer equations from the microscopic Stefan problem is a homogenization problem: averaging over the dendritic microstructure to obtain effective equations on the macroscopic scale.
where \(n\) is the spatial dimension. A common approximation is the volumetric average \(k_{\text{eff}} = \phi k_s + (1-\phi)k_l\).
Similar bounds exist for the solutal diffusivity and permeability, though the permeability is typically far more sensitive to the microstructural geometry. The Kozeny-Carman relation provides a useful empirical approximation, but more sophisticated models based on homogenization theory (e.g., using periodic cell problems) are needed for quantitative predictions.
This exceeds the critical value \(\mathrm{Ra}_{m,c} \approx 25\) found by Worster (1992), confirming the experimental observation of convective chimney plumes. Increasing the initial salt concentration or thickening the mush layer drives the system further into the convective regime, producing more vigorous plumes and more pronounced chimney channels.
4.5 Eutectic and Peritectic Solidification
Beyond simple binary alloys with a single mushy zone, more complex phase diagrams exhibit eutectic and peritectic reactions that produce qualitatively different solidification behavior.
The resulting microstructure consists of alternating lamellae (or rods) of \(\alpha\) and \(\beta\) with a characteristic spacing \(\lambda\) determined by the competition between interfacial energy and diffusion. The Jackson-Hunt theory (1966) predicts:
\[ \lambda^2 V = \text{const}, \]relating the lamellar spacing to the growth velocity.
Eutectic solidification involves two coupled free boundaries (the \(\alpha\)-liquid and \(\beta\)-liquid interfaces) that advance cooperatively. The solute rejected by one phase is consumed by the adjacent phase, creating a lateral diffusion coupling that stabilizes the lamellar pattern. The mathematical analysis involves solving a diffusion problem in the liquid with periodic boundary conditions matching the lamellar structure.
where the first term represents the solutal undercooling (proportional to the diffusion distance \(\lambda\)) and the second represents the capillary undercooling (inversely proportional to the radius of curvature, which scales with \(\lambda\)). The minimum undercooling occurs at \(\lambda_{\min} = \sqrt{K_2/(K_1 V)}\), giving the \(\lambda^2 V = \text{const}\) scaling. Experimentally, eutectics are observed to grow near (but not exactly at) the minimum-undercooling spacing, suggesting an extremum principle as the selection mechanism — an open theoretical question that has been addressed by stability analyses (Datye and Langer, 1981; Karma and Sarkissian, 1996).
Chapter 5: Variational Inequalities and Obstacle Problems
The Stefan problem admits an elegant reformulation as a variational inequality, connecting free boundary problems to the classical obstacle problem of potential theory. This perspective, developed extensively by Duvaut, Lions, Brezis, Caffarelli, and Kinderlehrer, provides the most powerful regularity theory for the free boundary and reveals deep connections to optimization, optimal control, and mathematical finance.
5.1 The Stefan Problem as a Variational Inequality
The key observation is that after a Baiocchi-type transformation, the Stefan problem becomes equivalent to an obstacle problem. Consider the one-phase Stefan problem in the enthalpy formulation.
If \(T\) solves the one-phase Stefan problem, then \(w\) satisfies:
\[ w \geq 0, \qquad -\Delta w + \frac{\partial w}{\partial t} \leq f, \qquad w\left(-\Delta w + \frac{\partial w}{\partial t} - f\right) = 0, \]where \(f\) encodes the initial and boundary data. The free boundary is \(\Gamma(t) = \partial\{w > 0\}\), the boundary of the positivity set of \(w\).
This reformulation recasts the free boundary as the boundary of a “contact set” — the region where the solution touches the obstacle \(w = 0\). The complementarity condition \(w(\mathcal{L}w - f) = 0\) states that either the solution is strictly positive (and the PDE holds with equality) or the solution is zero (and the PDE holds as an inequality). This structure is shared by many problems in mechanics (contact problems), finance (American options), and optimal control.
5.2 The Obstacle Problem
The obstacle problem is the simplest and most studied variational inequality. It serves as a model for all the analytical techniques used in free boundary regularity theory.
subject to \(u \geq \psi\) in \(\Omega\) and \(u = g\) on \(\partial\Omega\). The free boundary is \(\Gamma = \partial\{u > \psi\} \cap \Omega\), the boundary of the coincidence set \(\{u = \psi\}\).
The Euler-Lagrange conditions for the obstacle problem are:
\[ u \geq \psi, \qquad -\Delta u \geq 0, \qquad (u - \psi)(-\Delta u) = 0 \quad \text{in } \Omega. \]In the region where \(u > \psi\), the solution is harmonic. In the coincidence set \(\{u = \psi\}\), we have \(u = \psi\), and \(-\Delta u = -\Delta\psi \geq 0\). The transition between these two regimes occurs at the free boundary \(\Gamma\).
The obstacle problem admits an equivalent formulation as a variational inequality: find \(u \in K = \{v \in H^1_0(\Omega) : v \geq \psi\}\) such that
\[ \int_\Omega \nabla u \cdot \nabla(v - u)\, dx \geq 0 \quad \text{for all } v \in K. \]The existence and uniqueness of the solution follow from the Lions-Stampacchia theorem (1967), which guarantees that a coercive, continuous bilinear form on a closed convex subset of a Hilbert space has a unique minimizer. This abstract framework, developed by Lions and collaborators in the 1960s, unified a large class of problems in mechanics, physics, and optimization under the umbrella of variational inequalities.
with free boundary points at \(x = \pm 1\), where \(u\), \(u'\) are continuous but \(u''\) has a jump.
5.3 Regularity of the Free Boundary: Caffarelli’s Theory
The regularity theory for the free boundary in obstacle problems, developed by Luis Caffarelli beginning in the late 1970s, is one of the deepest achievements in the analysis of PDEs. The central question is: how smooth is the free boundary \(\Gamma = \partial\{u > \psi\}\)?
that is, \(u\) has Lipschitz continuous first derivatives. This regularity is optimal: \(u\) is generally not \(C^2\) across the free boundary.
Caffarelli showed that \(\Phi(r)\) is monotone nondecreasing in \(r\). Since \(u\) vanishes on the obstacle to first order (i.e., \(u(x) - \psi(x) = O(|x-x_0|^2)\) near \(x_0 \in \Gamma\)), the monotonicity formula implies \(|\nabla^2 u| \leq C\), giving \(C^{1,1}\) regularity. That this is sharp follows from the explicit one-dimensional solution, where \(u'' = -\Delta\psi\, \chi_{\{u = \psi\}}\) has a jump across the free boundary. \(\square\)
Beyond the regularity of the solution itself, Caffarelli classified the free boundary into regular and singular parts.
(i) The regular set \(\Gamma_{\text{reg}}\) is locally a \(C^{1,\alpha}\) surface (and \(C^\infty\) or even analytic if the data are smooth enough).
(ii) The singular set \(\Gamma_{\text{sing}}\) (where the coincidence set has zero density) is contained in a \(C^1\) manifold of dimension \(n - 1\), and is small in both topological and measure-theoretic senses.
The classification proceeds via blow-up analysis at free boundary points. At a regular point, the blow-up limit (rescaled solution as one zooms in) is a half-space solution \(u_0(x) = \frac{1}{2}\max(x \cdot e, 0)^2\), and the implicit function theorem gives smoothness. At a singular point, the blow-up limit is a homogeneous quadratic polynomial that vanishes on a lower-dimensional set.
The coincidence set \(\{u = 0\}\) converges to the half-space \(\{x_1 \leq 0\}\), confirming that the free boundary is locally a smooth curve near \(x_0\). At a singular point, the blow-up limit takes the form \(u_0 = \frac{1}{2}(a x_1^2 + b x_2^2)\) with \(a + b = \lambda\) (the Laplacian of the right-hand side) and \(a, b \geq 0\). When \(a > 0\) and \(b > 0\), the coincidence set degenerates to the single point \(\{0\}\), and the free boundary has a cusp. These singular points are isolated in two dimensions (Monneau, 2003), so the free boundary is a smooth curve except at finitely many cusp points.
5.4 The Thin Obstacle Problem
A natural generalization is the thin obstacle problem (also called the Signorini problem), where the obstacle constraint is imposed on a codimension-one surface rather than throughout the domain.
The thin obstacle problem arises in the study of semipermeable membranes, in the pricing of American options on assets with jump processes, and in the regularity theory for the fractional Laplacian. The optimal regularity is \(C^{1,1/2}\) (Athanasopoulos and Caffarelli), and the free boundary regularity theory is considerably more subtle than the classical obstacle problem due to the lower regularity and the possibility of polynomial degeneracies.
This problem was formulated by Antonio Signorini in 1933 and solved by Gaetano Fichera in 1963, making it one of the earliest variational inequalities to be studied rigorously. The Signorini problem has experienced a resurgence of interest due to its connection to the fractional obstacle problem: the extension technique of Caffarelli and Silvestre (2007) shows that the obstacle problem for the fractional Laplacian \((-\Delta)^s\) in \(\mathbb{R}^n\) is equivalent to a thin obstacle problem for a degenerate elliptic operator in \(\mathbb{R}^{n+1}\). This reformulation has enabled the application of Caffarelli’s regularity machinery to nonlocal free boundary problems, opening a rapidly growing area of analysis.
Chapter 6: Phase-Field Models
Phase-field models replace the sharp interface \(\Gamma(t)\) with a diffuse interfacial region of thickness \(\epsilon\), described by an order parameter \(\phi(x,t)\) that varies smoothly from one phase to another. This diffuse-interface approach eliminates the need to track the interface explicitly, making it particularly attractive for numerical simulation of complex morphologies (dendrites, grain boundaries, multi-phase systems). The mathematical challenge is to demonstrate that in the sharp-interface limit (\(\epsilon \to 0\)), the phase-field model recovers the original free boundary problem.
6.1 The Allen-Cahn Equation
The simplest phase-field model is the Allen-Cahn equation, introduced by Allen and Cahn (1979) to model the motion of antiphase boundaries in crystalline solids. It is the \(L^2\) gradient flow of the Ginzburg-Landau free energy.
where \(\phi(x) \in [-1, 1]\) is the order parameter (\(\phi = +1\) in one phase, \(\phi = -1\) in the other), \(\epsilon > 0\) controls the interface thickness, and \(W(\phi)\) is a double-well potential with minima at \(\phi = \pm 1\). The standard choice is \(W(\phi) = \frac{1}{4}(1 - \phi^2)^2\).
where \(\tau > 0\) is a relaxation time. The equilibrium interface profile is the heteroclinic connection:
\[ \phi_0(z) = \tanh\left(\frac{z}{\sqrt{2}\,\epsilon}\right), \]where \(z\) is the signed distance to the interface. The interface width scales as \(\epsilon\), and the interfacial energy is
\[ \sigma = \int_{-\infty}^{\infty} \epsilon|\phi_0'|^2 + \frac{1}{\epsilon}W(\phi_0)\, dz = \frac{2\sqrt{2}}{3}\epsilon. \]The Allen-Cahn equation does not conserve the total amount of each phase — the volume of one phase can shrink at the expense of the other. This is appropriate for order-disorder transitions but not for systems where mass is conserved (e.g., spinodal decomposition in a binary mixture).
6.2 The Cahn-Hilliard Equation
For conserved order parameters, the appropriate gradient flow is the \(H^{-1}\) gradient flow of \(\mathcal{F}\), which yields the Cahn-Hilliard equation.
where \(\mu\) is the chemical potential and \(M > 0\) is a mobility coefficient. The total order parameter \(\int_\Omega \phi\, dx\) is conserved. The equation is fourth-order in space, requiring two boundary conditions on each boundary (typically \(\partial\phi/\partial n = \partial\mu/\partial n = 0\) for no-flux conditions).
The Cahn-Hilliard equation was derived by Cahn and Hilliard (1958) to model spinodal decomposition: the spontaneous phase separation of a binary mixture quenched below its critical temperature. The equation predicts the formation of a fine-scale lamellar or droplet pattern that subsequently coarsens — a process beautifully captured by the Lifshitz-Slyozov-Wagner (LSW) theory of Ostwald ripening, which predicts the characteristic length scale grows as \(t^{1/3}\).
Since \(W''(C_0) < 0\), modes with \(k < k_c = \sqrt{|W''(C_0)|}/\epsilon\) are unstable, and the most dangerous wavenumber is \(k_{\max} = k_c/\sqrt{2}\), giving the initial characteristic length scale \(\lambda_{\max} = 2\pi\sqrt{2}\,\epsilon/\sqrt{|W''(C_0)|}\). For a Fe-Cr alloy at 475\(^\circ\)C (the “475\(^\circ\)C embrittlement” regime), \(\epsilon \sim 1\) nm and \(|W''| \sim 10^8\) J/m\(^3\), giving \(\lambda_{\max} \sim 10\) nm, consistent with atom probe tomography measurements of the early stages of spinodal decomposition.
6.3 Sharp-Interface Limit via Matched Asymptotics
The central theoretical result connecting phase-field models to the classical free boundary problem is the sharp-interface limit, obtained by matched asymptotic analysis as \(\epsilon \to 0\). This was first carried out rigorously by Caginalp (1989) and subsequently refined by many authors.
where \(u = (T - T_m)/(L/c)\) is the dimensionless temperature. In the limit \(\epsilon \to 0\) with \(\tau, \lambda\) appropriately scaled, the phase-field model reduces to the two-phase Stefan problem with Gibbs-Thomson condition:
\[ \kappa\nabla^2 u = \frac{\partial u}{\partial t} \quad \text{in each phase}, \]\[ u = -\sigma_{sl}\kappa/(\rho L) - \beta\, v_n \quad \text{on } \Gamma(t), \]\[ L\rho\, v_n = [k\nabla u \cdot \hat{n}] \quad \text{on } \Gamma(t), \]where \(\kappa\) is the interface curvature, \(v_n\) the normal velocity, and \(\beta\) a kinetic coefficient.
The matched-asymptotic approach is powerful but formal. Rigorous justification of the sharp-interface limit has been achieved in various settings using techniques from geometric measure theory and viscosity solutions (Chen, 1992; de Mottoni and Schatzman, 1995; Ilmanen, 1993).
where \(g\) contains terms involving the curvature \(\kappa\), the normal velocity \(v_n\), and the outer temperature \(u_0\). The operator \(\mathcal{L}\) is self-adjoint with a one-dimensional null space spanned by \(\phi_0'(\xi) = \operatorname{sech}^2(\xi/\sqrt{2})/\sqrt{2}\) (the translation mode, arising from the invariance of the equilibrium profile under shifts). By the Fredholm alternative, a bounded solution \(\phi_1\) exists if and only if \(g\) is orthogonal to the null space:
\[ \int_{-\infty}^{\infty} g(\xi)\,\phi_0'(\xi)\, d\xi = 0. \]Evaluating this integral yields \(u_0|_{\Gamma} = -\sigma\kappa/(\rho L) - \beta v_n\), which is precisely the Gibbs-Thomson condition with kinetic undercooling. The interfacial energy \(\sigma\) and kinetic coefficient \(\beta\) are expressed as integrals of the equilibrium profile, providing an explicit connection between the phase-field parameters and the sharp-interface quantities.
6.4 Quantitative Phase-Field Models
A key practical issue is that standard phase-field models introduce spurious interface kinetics (the \(\beta\, v_n\) term) that are absent in the physical problem. The quantitative phase-field model of Karma and Rappel (1996, 1998) eliminates this artifact by a clever choice of the coupling function and time scale.
is added to the solute (or temperature) equation to cancel the anomalous interface effects that arise at finite \(\epsilon\). This allows computations with interface widths much larger than the physical capillary length, dramatically reducing computational cost while maintaining quantitative accuracy.
6.5 Multi-Phase and Multi-Component Extensions
Phase-field models extend naturally to multi-phase systems (e.g., grain growth, where many grains with different crystallographic orientations coexist) and multi-component alloys.
with evolution equations respecting the constraint. This formulation, due to Steinbach, Nestler, and others, allows independent specification of interfacial energies \(\sigma_{ij}\) for each pair of phases.
Multi-phase field models are now the standard computational tool for simulating microstructure evolution in multi-component alloys, eutectic solidification, and grain growth. The mathematical analysis of these systems — existence, regularity, and sharp-interface limits — remains an active area of research.
Chapter 7: Level Set and Front-Tracking Methods
Numerical solution of free boundary problems requires a method for representing and evolving the interface. Broadly, there are two approaches: Lagrangian methods, which track the interface explicitly using markers or meshes, and Eulerian methods, which represent the interface implicitly as a level set of a function defined on a fixed grid. Each approach has strengths and limitations, and the choice depends on the problem at hand.
7.1 The Level Set Method
The level set method, introduced by Osher and Sethian (1988), represents the interface \(\Gamma(t)\) as the zero level set of a function \(\phi(x,t)\):
\[ \Gamma(t) = \{x : \phi(x,t) = 0\}. \]The function \(\phi\) is typically initialized as the signed distance to the interface: \(\phi > 0\) in one phase, \(\phi < 0\) in the other.
or equivalently, if the velocity field \(\mathbf{v}\) is extended to a neighborhood of the interface:
\[ \frac{\partial\phi}{\partial t} + \mathbf{v} \cdot \nabla\phi = 0. \]Geometric quantities are recovered from \(\phi\): the unit normal is \(\hat{n} = \nabla\phi/|\nabla\phi|\), and the mean curvature is \(\kappa = \nabla \cdot (\nabla\phi/|\nabla\phi|)\).
The great advantage of the level set method is that topological changes (merging and splitting of interfaces) are handled automatically, with no special logic required. The interface is never explicitly parameterized, so there is no mesh tangling or marker redistribution.
7.2 Hamilton-Jacobi Equations and Numerical Schemes
The level set equation is a Hamilton-Jacobi equation, and its numerical treatment requires special care to handle the development of kinks (discontinuities in \(\nabla\phi\)) in finite time, even from smooth initial data.
The standard numerical discretization uses upwind schemes from hyperbolic conservation laws. The first-order scheme uses the Godunov Hamiltonian; higher-order accuracy is achieved via ENO (essentially non-oscillatory) or WENO (weighted essentially non-oscillatory) spatial discretizations combined with TVD (total variation diminishing) Runge-Kutta time integration.
A circle of initial radius \(R_0\) shrinks as \(R(t) = \sqrt{R_0^2 - 2(n-1)t}\) and collapses to a point at time \(t^* = R_0^2/(2(n-1))\) in \(\mathbb{R}^n\). This provides a benchmark for level set codes. Note the connection to the Allen-Cahn equation: in the sharp-interface limit \(\epsilon \to 0\), the Allen-Cahn equation produces motion by mean curvature.
7.3 Reinitialization
A practical difficulty with the level set method is that the function \(\phi\) ceases to be a signed distance function as it evolves, even if initialized as one. Steep and flat gradients develop, degrading the accuracy of curvature computations and other geometric quantities that depend on \(|\nabla\phi| = 1\).
where \(\phi_0 = \phi(x, t)\) is the level set function at the current physical time. This is an Eikonal-type equation that redistributes the level sets of \(\phi\) to be equally spaced (in the normal direction) without moving the zero level set.
7.4 Volume-of-Fluid and Front-Tracking Methods
The level set method is not the only approach to interface tracking. Two important alternatives are the volume-of-fluid (VOF) method and front-tracking methods.
The interface is reconstructed in each cell using a piecewise-linear interface calculation (PLIC).
The VOF method has the important property of exact volume conservation (by construction), which the level set method does not guarantee without additional corrections. This makes VOF the method of choice for applications where mass conservation is critical, such as multiphase flow with surface tension.
7.5 Ghost Fluid and Immersed Interface Methods
An important class of methods that bridges the gap between sharp and diffuse interface approaches is the ghost fluid method (GFM) of Fedkiw, Aslam, Merriman, and Osher (1999), designed to impose jump conditions sharply on a fixed Cartesian grid.
and the Gibbs-Thomson condition \(T = T_m - \Gamma_{\text{GT}}\kappa\) are imposed to first- or second-order accuracy without smearing. The interface location within each cell is determined from the level set function, and the jump in the normal derivative of \(T\) is computed from the Stefan condition.
The GFM achieves sharp resolution of the interface conditions on a Cartesian grid, avoiding the diffusion of the interface inherent in enthalpy methods. The closely related immersed interface method (IIM) of LeVeque and Li (1994) takes a more systematic approach, modifying the finite difference coefficients near the interface to maintain uniform accuracy throughout the domain, including at irregular grid points adjacent to the interface.
7.6 Comparison of Methods
Each interface-tracking method has distinct advantages and disadvantages:
The level set method handles topological changes naturally and provides easy access to geometric quantities (normal, curvature) but does not conserve volume exactly and requires reinitialization. The VOF method conserves volume exactly but provides poor accuracy for curvature computation and requires complex geometric reconstruction algorithms. Front-tracking methods provide high accuracy for smooth interfaces and exact representation of surface tension but struggle with topological changes (merging and breakup require explicit surgical procedures). Hybrid methods, such as the coupled level set/VOF (CLSVOF) method, attempt to combine the advantages of both approaches.
For phase-change problems specifically, the choice of interface method is often secondary to the treatment of the Stefan condition. The enthalpy method (Chapter 1) avoids explicit interface tracking entirely by smearing the latent heat release over a finite temperature interval, while the phase-field method (Chapter 6) replaces the sharp interface with a diffuse one. These approaches are generally preferred for solidification problems with complex morphologies.
Chapter 8: Applications: Crystal Growth, Ice, and Industrial Processes
Free boundary problems are not merely mathematical abstractions — they arise in a vast range of physical settings, from semiconductor manufacturing to climate science. In this final chapter, we survey several important applications, emphasizing how the mathematical tools developed in the preceding chapters inform the understanding of each physical system.
8.1 Crystal Growth from the Melt
The growth of single crystals from the melt is a cornerstone of the semiconductor industry. The Czochralski process, invented in 1916, remains the dominant method for producing silicon wafers.
(i) The crystal-melt interface (a Stefan problem with convection);
(ii) The melt-gas and crystal-gas free surfaces (capillary surfaces determined by the Young-Laplace equation);
(iii) The thermal field in the melt (Navier-Stokes with buoyancy, Marangoni, and forced convection).
The shape of the crystal-melt interface in Czochralski growth is critical for crystal quality. A convex interface (bowing into the melt) promotes the incorporation of dislocations and thermal stresses, while a flat or slightly concave interface is preferred. Controlling the interface shape through manipulation of the thermal environment is a central goal of crystal growth engineering.
8.2 Snowflake Formation
The formation of snowflakes (ice crystals growing from water vapor) is one of the most beautiful examples of dendritic growth in nature. The morphology of a snowflake — whether it forms plates, columns, needles, or dendrites — depends sensitively on the temperature and supersaturation of the ambient vapor, as summarized in the Nakaya diagram (1954). Ukichiro Nakaya, working at Hokkaido University in the 1930s, was the first to grow artificial snow crystals under controlled conditions, systematically mapping the crystal habit as a function of temperature and supersaturation. His observation that “snow crystals are letters sent from heaven” captured the public imagination and remains a motivating image for the mathematical study of pattern formation. The physical origin of the dramatic morphological transitions in the Nakaya diagram (e.g., plates at \(-2^\circ\)C, columns at \(-5^\circ\)C, dendrites at \(-15^\circ\)C) is the strong temperature dependence of the ice crystal’s surface energy anisotropy and attachment kinetics, which change the relative growth rates of the basal and prism faces.
where \(c\) is the water vapor concentration, \(c_{\text{eq}}\) incorporates the Gibbs-Thomson effect and kinetic anisotropy, and \(D\) is the vapor diffusivity. The sixfold symmetry of ice (\(I_h\) lattice) enters through the anisotropy of the surface energy and kinetic coefficient, producing the characteristic hexagonal symmetry of snowflakes. Phase-field simulations by Gravner and Griffeath (2009) and others have reproduced the full diversity of the Nakaya diagram.
8.3 Sea Ice Growth
Sea ice covers approximately 7% of the Earth’s surface at its maximum extent and plays a crucial role in the global climate system through the ice-albedo feedback, modification of ocean-atmosphere heat exchange, and brine rejection into the ocean. The growth of sea ice is a mushy-layer problem (Chapter 4), with seawater as the binary alloy (water + NaCl).
(i) Heat conduction through the ice, with Stefan condition at the ice-ocean interface;
(ii) Brine channel formation and drainage (chimney convection in the mushy layer);
(iii) Snow cover as an insulating layer with variable thermal conductivity;
(iv) Radiative forcing at the upper surface. The classic Stefan model (with appropriate effective thermal properties) predicts ice thickness growing as \[ h(t) \approx \sqrt{\frac{2k_i(T_m - T_a)t}{\rho L}}, \]
where \(T_a\) is the atmospheric temperature. This \(\sqrt{t}\) law (cf. Chapter 1) was verified by early Arctic observations and remains the basis of simple sea ice models.
8.4 Continuous Casting
Continuous casting is the process by which most of the world’s steel is produced — over 95% of the world’s 1.9 billion tonnes of annual steel output passes through a continuous caster. Molten steel is poured into a water-cooled copper mold, and the partially solidified strand is continuously withdrawn at speeds of 1–2 m/min. The solidification front (mush-liquid boundary) is a free boundary whose position determines the shell thickness and hence the mechanical integrity of the strand. If the shell is too thin when the strand exits the mold, a breakout occurs — the liquid steel bursts through the shell, causing equipment damage and production shutdowns costing millions of dollars. The mathematical prediction of shell thickness using the Stefan model is therefore of direct industrial importance.
where \(K\) is a solidification constant depending on the steel grade, mold heat transfer, and superheat. This is the Neumann \(\sqrt{t}\) law (Chapter 1) with \(t = z/V\), reflecting the dominant role of diffusion. Deviations from \(\sqrt{z}\) growth indicate the onset of convection in the liquid pool, non-uniform mold heat transfer, or the formation of an air gap between the shell and the mold.
8.5 Injection Molding and Additive Manufacturing
Free boundary problems also arise in polymer processing (injection molding, where a viscous polymer melt fills a mold cavity) and additive manufacturing (3D printing, where a heat source selectively melts and resolidifies material).
In selective laser melting (SLM), a laser beam scans a powder bed, creating a small melt pool that solidifies rapidly behind the beam. The melt pool boundary is a free boundary governed by the Stefan problem with a moving heat source (the laser), Marangoni convection in the melt pool, phase change (melting and resolidification), and possible evaporation at high laser powers. The mathematical formulation couples the heat equation with a moving source term \(Q(x - V_s t)\) (where \(V_s\) is the scan velocity and \(Q\) represents the laser energy deposition profile, typically Gaussian) to the Navier-Stokes equations in the melt pool, with the melt pool boundary itself being a free boundary determined by the solidus isotherm.
The rapid solidification rates (\(10^5\)-\(10^7\) K/s) place the process far from equilibrium, and the Gibbs-Thomson condition must be supplemented with kinetic undercooling and solute trapping effects. Understanding the melt pool geometry and the resulting microstructure is essential for controlling the quality of additively manufactured parts.
where \(k \approx 7\) W/(m\(\cdot\)K) and \(\kappa \approx 3 \times 10^{-6}\) m\(^2\)/s for Ti-6Al-4V. The melt pool boundary is the isotherm \(T = T_m = 1660^\circ\)C. For the parameters above, the melt pool is approximately 150 \(\mu\)m long, 80 \(\mu\)m wide, and 40 \(\mu\)m deep — an elongated teardrop shape. The solidification front at the trailing edge of the melt pool experiences cooling rates of order \(V_s \times G \sim 10^6\) K/s, producing an extremely fine microstructure with columnar \(\beta\)-grains oriented along the steepest thermal gradient.
8.6 Tumor Growth as a Free Boundary Problem
In a striking application outside traditional materials science, the growth of solid tumors has been modeled as a free boundary problem since the work of Greenspan (1972) and Byrne and Chaplain (1996).
and the tumor boundary moves according to:
\[ v_n = -\frac{\partial p}{\partial n}\bigg|_\Gamma, \qquad -\nabla^2 p = S(c) \quad \text{in } \Omega(t), \]where \(p\) is the internal pressure and \(S(c)\) is the net proliferation rate (cell birth minus death). When \(S(c) = \mu(c - \tilde{c})\), with \(\tilde{c}\) a threshold, the model exhibits a steady spherical solution whose stability can be analyzed by the Mullins-Sekerka-type analysis of Chapter 3.
The tumor growth model is structurally similar to the Hele-Shaw problem (Chapter 2): the pressure is determined by a Poisson equation inside the domain, and the boundary moves with velocity proportional to the pressure gradient. The morphological instability of the tumor boundary (analogous to the Saffman-Taylor instability) may be relevant to the development of invasive morphologies. Recent work has incorporated mechanics, angiogenesis, and the immune response, leading to increasingly realistic free boundary models of cancer growth.
8.7 Geophysical Free Boundary Problems
Beyond sea ice, free boundary problems arise throughout the geosciences. The solidification of the Earth’s inner core is a planetary-scale Stefan problem: the inner core grows by freezing from the iron-rich liquid outer core, releasing latent heat and light elements (silicon, oxygen, sulfur) that drive the compositional convection powering the geodynamo. The growth rate of the inner core — estimated at roughly 1 mm per year — is controlled by the balance between the secular cooling of the core and the release of latent heat, precisely the Stefan condition applied at the inner core boundary.
8.8 Outlook
Free boundary problems remain at the frontier of applied mathematics, with active research in:
- Regularity theory for the Stefan problem in higher dimensions and for systems (Caffarelli, Salsa, Shahgholian, and their schools);
- Quantitative phase-field models for multi-component, multi-phase systems (Karma, Plapp, Steinbach);
- Data-driven and machine-learning approaches to free boundary problems (physics-informed neural networks for Stefan problems, neural operators for parametric free boundary problems);
- Stochastic free boundary problems (random interfaces, fluctuating hydrodynamics at phase boundaries);
- Free boundary problems in biology (wound healing, tissue growth, biofilm formation);
- Shape optimization and inverse free boundary problems (determining material properties or boundary conditions from observed interface motion);
- Free boundary problems on evolving surfaces and in complex geometries (phase transitions on curved membranes, solidification in porous media).
The mathematical framework developed in these notes — the interplay of PDE theory, variational principles, asymptotic analysis, and numerical methods — provides the essential toolkit for all of these directions. The field’s vitality is reflected in the continuing stream of Fields Medal and Abel Prize citations that reference free boundary problems: Caffarelli’s 2023 Abel Prize for regularity theory, Figalli’s 2018 Fields Medal for work on optimal transport (closely connected to obstacle problems), and Villani’s 2010 Fields Medal for contributions to kinetic theory and entropy methods that apply to phase-field models. Free boundary problems remain a meeting ground where pure analysis, physical modeling, and computational science inform and enrich one another.
