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.
Definition 1.1 (Classical Two-Phase Stefan Problem). Find the temperature \(T(x,t)\) and the free boundary \(\Gamma(t)\) such that:
\[
\rho_i c_i \frac{\partial T}{\partial t} = \nabla \cdot (k_i \nabla T), \quad x \in \Omega_i(t), \quad i \in \{s, l\},
\]
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.
Definition 1.2 (One-Phase Stefan Problem). When the temperature in one phase is uniformly at the melting point \(T_m\), only one heat equation needs to be solved. For solidification into a supercooled liquid at temperature \(T_\infty < T_m\), we solve:
\[
\frac{\partial T}{\partial t} = \kappa \nabla^2 T, \quad x \in \Omega_l(t),
\]
\[
T = T_m \quad \text{on } \Gamma(t), \qquad L\rho\, v_n = -k \frac{\partial T}{\partial n}\bigg|_l \quad \text{on } \Gamma(t).
\]
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
\[
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.
Theorem 1.1 (Stefan Condition from Energy Conservation). Let \(\Gamma(t)\) be a smooth interface separating solid and liquid phases, with unit normal \(\hat{n}\) pointing into the liquid. If \(T\) satisfies the heat equation in each phase and is continuous across \(\Gamma(t)\), then energy conservation at the interface implies
\[
L\rho\, v_n = [k \nabla T \cdot \hat{n}]_s^l = k_s \nabla T_s \cdot \hat{n} - k_l \nabla T_l \cdot \hat{n},
\]
where \(v_n\) is the normal velocity of \(\Gamma(t)\) in the direction of \(\hat{n}\).
Proof. Apply the Reynolds transport theorem to the energy integral over a control volume \(V(t)\) containing a segment of the interface. The internal energy density has a jump \([\rho e] = \rho L\) across \(\Gamma(t)\). Writing the energy balance:
\[
\frac{d}{dt} \int_{V(t)} \rho\, e\, dV = -\oint_{\partial V(t)} \mathbf{q} \cdot \hat{n}\, dS,
\]
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.
Example 1.1 (Neumann Solution for Solidification in a Half-Space). Seek a similarity solution with the interface at \(x = s(t) = 2\lambda\sqrt{\kappa_s t}\), where \(\lambda\) is a dimensionless constant. Introducing the similarity variable \(\eta = x/(2\sqrt{\kappa_i t})\) in each phase:
\[
T_s(\eta) = T_w + (T_m - T_w)\frac{\operatorname{erf}(\eta)}{\operatorname{erf}(\lambda)}, \quad 0 < \eta < \lambda,
\]
\[
T_l(\eta) = T_0 - (T_0 - T_m)\frac{\operatorname{erfc}(\eta\sqrt{\kappa_s/\kappa_l})}{\operatorname{erfc}(\lambda\sqrt{\kappa_s/\kappa_l})}, \quad \eta > \lambda.
\]
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).
Theorem 1.2 (Existence and Uniqueness — Friedman). For the one-phase Stefan problem in one space dimension with smooth initial data and compatible boundary conditions, there exists a unique classical solution \((T, s) \in C^{2,1}(\Omega) \times C^1([0,T_0])\) for some \(T_0 > 0\). The free boundary \(s(t)\) is a \(C^\infty\) function of time for \(t > 0\).
In higher dimensions, the classical theory becomes more delicate. Interfaces can develop corners, cusps, and topological changes. This motivates the weak (enthalpy) 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.
Definition 1.3 (Enthalpy Formulation). Define the enthalpy as:
\[
H(T) = \begin{cases} \rho c_s T & \text{if } T < T_m, \\ \rho c_s T_m + \rho L \xi & \text{if } T = T_m, \; \xi \in [0,1], \\ \rho c_l T + \rho L & \text{if } T > T_m, \end{cases}
\]
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.
Theorem 1.3 (Existence of Weak Solutions — Kamenomostskaya/Oleinik). For the Stefan problem in \(\mathbb{R}^n\) with \(L^2\) initial data and appropriate boundary conditions, there exists a unique weak solution in the enthalpy sense. The temperature \(T\) is continuous, and the enthalpy \(H\) belongs to \(L^\infty(0,T; L^2(\Omega))\).
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.
Example 1.4 (Enthalpy Method on a Uniform Grid). Consider the one-dimensional Stefan problem on \([0, 1]\) with \(T(0,t) = -1\), \(T(1,t) = 1\), and initial temperature \(T(x,0) = 2x - 1\), so the initial interface is at \(x = 1/2\). Let \(T_m = 0\), \(\kappa = 1\), and \(\rho L = 1\). Discretize with \(N\) cells of width \(\Delta x = 1/N\) and time step \(\Delta t\). The implicit enthalpy scheme reads:
\[
\frac{H_j^{n+1} - H_j^n}{\Delta t} = \frac{T_{j+1}^{n+1} - 2T_j^{n+1} + T_{j-1}^{n+1}}{\Delta x^2},
\]
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.
Definition 1.4 (Stefan Number). The Stefan number is defined as:
\[
\mathrm{St} = \frac{c\,\Delta T}{L},
\]
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.
Example 1.2 (Quasi-Steady Approximation). For \(\mathrm{St} \ll 1\) in one dimension with a cold wall at \(x = 0\) and melting temperature at the interface \(x = s(t)\), the quasi-steady temperature profile is linear:
\[
T(x,t) \approx T_w + (T_m - T_w)\frac{x}{s(t)}, \quad 0 < x < s(t).
\]
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}\).
Example 1.5 (Perturbation Expansion in the Stefan Number). For the one-phase problem in a half-space with \(\mathrm{St} \ll 1\), write the interface position as \(s(t) = s_0(t) + \mathrm{St}\, s_1(t) + O(\mathrm{St}^2)\). The leading-order term is the quasi-steady result \(s_0 = \sqrt{2\,\mathrm{St}\,\kappa\, t}\). At next order, the sensible heat stored in the solid (which the quasi-steady approximation neglects) contributes a correction. Expanding the exact Neumann transcendental equation for small \(\lambda\):
\[
\lambda = \sqrt{\mathrm{St}/\pi}\left(1 + \frac{\mathrm{St}}{2\pi} + \frac{\mathrm{St}^2}{3\pi^2} + \cdots\right),
\]
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.
Theorem 1.4 (Comparison Principle for the Stefan Problem). Let \((T_1, \Gamma_1)\) and \((T_2, \Gamma_2)\) be two solutions of the one-phase Stefan problem with initial interfaces \(\Gamma_1(0) \subset \Gamma_2(0)\) (i.e., the solid region for solution 1 is contained in the solid region for solution 2) and with \(T_1 \leq T_2\) initially. Then \(\Gamma_1(t) \subset \Gamma_2(t)\) for all \(t > 0\): the ordering of free boundaries is preserved in time.
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.
Example 1.3 (Bounding the Freezing Time of a Slab). Consider a slab of water of thickness \(2a\) initially at the melting temperature \(T_m = 0^\circ\)C, with both faces held at temperature \(T_w < 0\). The two solidification fronts advance inward from \(x = 0\) and \(x = 2a\). By the comparison principle, each front is bounded above by the Neumann solution for a half-space (which ignores the finite slab width) and bounded below by the quasi-steady approximation. The slab is completely frozen at time \(t_f\) satisfying:
\[
\frac{a^2}{2\,\mathrm{St}\,\kappa} \leq t_f \leq \frac{a^2}{4\lambda^2\kappa},
\]
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.
Definition 2.1 (Hele-Shaw Flow). In a Hele-Shaw cell with gap width \(b\), the gap-averaged velocity \(\mathbf{u}\) and pressure \(p\) satisfy:
\[
\mathbf{u} = -\frac{b^2}{12\mu}\nabla p, \qquad \nabla \cdot \mathbf{u} = 0,
\]
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.
Theorem 2.1 (Saffman-Taylor Instability). Consider a planar interface \(y = Y(t)\) in a Hele-Shaw cell of width \(W\), with the less viscous fluid above and the more viscous fluid below. A sinusoidal perturbation of wavenumber \(k\):
\[
y = Y(t) + \hat{\epsilon}\, e^{\sigma t}\cos(kx)
\]
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}\).
Proof. The base state is a planar interface advancing at speed \(V_0\), with pressure \(p_0 = -12\mu V_0 y/b^2\) in the viscous fluid. Perturbing the interface to \(y = Y_0(t) + \hat{\epsilon}\,e^{ikx}\) and linearizing Laplace's equation for the pressure perturbation \(\hat{p}\) gives \(\hat{p} = A\,e^{-|k|(y - Y_0)}\). The kinematic condition and the pressure boundary condition (with curvature) yield
\[
\sigma = V_0|k| - \frac{b^2\sigma_{\text{st}}}{12\mu}|k|^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.
Example 2.3 (Viscous Fingering in Oil Recovery). In secondary oil recovery, water is injected into a porous reservoir to displace the more viscous oil toward production wells. The displacement is governed by Darcy's law, making it mathematically identical to a Hele-Shaw flow. When the mobility ratio \(M = (\mu_{\text{oil}}/k_{\text{oil}})/(\mu_{\text{water}}/k_{\text{water}}) > 1\), the displacement front is Saffman-Taylor unstable, and viscous fingers of water bypass large regions of oil, dramatically reducing the sweep efficiency. For a typical waterflood with \(M \approx 10\), the linear theory predicts instability at all wavelengths larger than a capillary length of order millimeters, and nonlinear finger competition leads to a few dominant fingers that control the breakthrough time. Polymer flooding (adding polymer to the injection water to increase its viscosity and reduce \(M\)) is one of the primary enhanced oil recovery strategies, motivated directly by the need to suppress this instability.
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.
Theorem 2.2 (Richardson's Moment Theorem). For zero-surface-tension Hele-Shaw flow driven by injection at the origin with flux \(Q\), the moments
\[
M_k(t) = \iint_{\Omega(t)} z^k\, dA, \quad k = 0, 1, 2, \ldots,
\]
(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.
Definition 2.2 (Polubarinova-Galin Equation). Let \(z = f(\zeta, t)\) be the conformal map from the unit disk \(|\zeta| < 1\) to the fluid domain \(\Omega(t)\), normalized so that \(f(0,t) = 0\) and \(f'(0,t) > 0\). Then
\[
\operatorname{Re}\left[\dot{f}(\zeta,t)\,\overline{\zeta f'(\zeta,t)}\right] = \frac{Q}{2\pi}, \quad |\zeta| = 1.
\]
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)\).
Example 2.1 (Exact Solution: Expanding Circle). For injection at the origin with flux \(Q\) into an initially circular domain of radius \(R_0\), the conformal map is \(f(\zeta, t) = R(t)\zeta\) and the Polubarinova-Galin equation reduces to:
\[
R\dot{R} = \frac{Q}{2\pi}, \qquad \Rightarrow \qquad R(t) = \sqrt{R_0^2 + \frac{Qt}{\pi}}.
\]
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.
Example 2.2 (Polynomial Conformal Maps and Finite-Time Singularities). Consider a conformal map of the form \(f(\zeta, t) = a_1(t)\zeta + a_N(t)\zeta^N\). The Polubarinova-Galin equation reduces to a system of ODEs for the coefficients \(a_1(t)\) and \(a_N(t)\). For certain initial conditions, the map ceases to be univalent (one-to-one) in finite time, corresponding to cusp formation on the free boundary. This finite-time blow-up is a manifestation of the ill-posedness of the zero-surface-tension Hele-Shaw problem and demonstrates the need for regularization (surface tension, kinetic undercooling, or a mushy zone).
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
\[
\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.
Theorem 3.1 (Mullins-Sekerka Instability). A sinusoidal perturbation \(\hat{\epsilon}\,e^{\sigma t + ikx}\) of the planar front has growth rate
\[
\sigma(k) = V|k|\left(\frac{G_l}{V|k|/\kappa + V/\kappa} - \frac{\Gamma_{\text{GT}} k^2}{V|k|/\kappa + V/\kappa}\right) - V|k|\frac{G_s}{V|k|/\kappa + V/\kappa},
\]
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.
Example 3.2 (Critical Wavelength for Ice Solidification). For ice solidifying from pure water at modest undercoolings (\(\Delta T \approx 1\) K), the capillary length is \(\Gamma_{\text{GT}} \approx 2.7 \times 10^{-8}\) m, the thermal diffusivity is \(\kappa \approx 1.3 \times 10^{-7}\) m\(^2\)/s, and a typical growth velocity is \(V \approx 10^{-5}\) m/s. The critical wavelength is:
\[
\lambda_c = \frac{2\pi}{k_c} = 2\pi\sqrt{\frac{2\kappa\Gamma_{\text{GT}}}{V}} \approx 2\pi\sqrt{\frac{2 \times 1.3 \times 10^{-7} \times 2.7 \times 10^{-8}}{10^{-5}}} \approx 5 \times 10^{-5} \text{ m} = 50 \;\mu\text{m}.
\]
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.
Definition 3.1 (Gibbs-Thomson Condition). The equilibrium temperature of a curved solid-liquid interface with mean curvature \(\kappa\) is:
\[
T_{\Gamma} = T_m - \Gamma_{\text{GT}} \kappa - \frac{V}{\mu_k},
\]
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.
Proof. (Derivation of the Gibbs-Thomson Condition from Thermodynamics.) Consider a small patch of the solid-liquid interface with area \(A\) and mean curvature \(\kappa\). The Gibbs free energy of the system includes a bulk term and a surface term:
\[
G = G_{\text{bulk}} + \sigma_{sl} A.
\]
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.
Example 3.1 (Ivantsov Paraboloid Solution). Ivantsov (1947) found an exact steady-state solution for a paraboloidal dendrite tip growing at velocity \(V\) into a supercooled melt. In parabolic coordinates centered on the tip, the temperature field satisfies:
\[
\operatorname{Iv}(\operatorname{Pe}) = \operatorname{Pe}\, e^{\operatorname{Pe}} E_1(\operatorname{Pe}) = \frac{\Delta T c_l}{L} \equiv \mathrm{St},
\]
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).
Theorem 3.2 (Solvability Theory). For an anisotropic Gibbs-Thomson condition of the form \(\Gamma_{\text{GT}}(\theta) = \Gamma_0(1 - \epsilon_4 \cos 4\theta)\) with fourfold crystalline anisotropy \(\epsilon_4 > 0\), there exists a discrete set of steady-state needle crystal solutions. Only the fastest-growing solution (with the largest tip velocity) is linearly stable. The selected operating point satisfies
\[
\sigma^* = \frac{2\kappa d_0}{V R^2} = \text{const}(\epsilon_4),
\]
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.
Example 3.3 (Tip Selection in Succinonitrile). Succinonitrile (SCN) is the preferred model material for dendritic growth experiments because it is transparent, has a low melting point (\(T_m = 58.08^\circ\)C), and solidifies with cubic symmetry like metals. The IDGE (Isothermal Dendritic Growth Experiment) conducted aboard the Space Shuttle by Glicksman, Koss, and Winsa (1994) measured dendrite tip velocities and radii in microgravity (eliminating convective transport) at various undercoolings. For an undercooling of \(\Delta T = 1\) K, the measured tip velocity was \(V \approx 6.5\;\mu\)m/s with tip radius \(R \approx 21\;\mu\)m, giving a stability parameter \(\sigma^* = 2\kappa d_0/(VR^2) \approx 0.020\). This agreed with the solvability theory prediction to within experimental uncertainty, providing definitive validation of microscopic solvability as the selection mechanism for dendritic growth.
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.
Proposition 3.1. In the WKB approximation, a perturbation created at the tip with wavenumber \(k\) and amplitude \(A_0\) is amplified as it travels along the dendrite surface to position \(s\) (arc length from the tip) as:
\[
A(s) = A_0 \exp\left(\int_0^s \sigma(k(s'), s')\, \frac{ds'}{V_{\text{tip}}}\right),
\]
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.
Definition 4.1 (Liquidus and Solidus). For a binary alloy with composition \(C\) (mass fraction of solute B):
\[
T_L(C) = T_m - m_L C \quad \text{(liquidus)}, \qquad T_S(C) = T_m - \frac{m_L}{k_p} C \quad \text{(solidus)},
\]
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.
Example 4.2 (Pb-Sn Solder Alloy). The Pb-Sn system is the classic example used in introductory metallurgy courses. With \(T_m(\text{Pb}) = 327^\circ\)C, \(m_L \approx 2.3\) K/wt%, and \(k_p \approx 0.31\) for dilute Sn in Pb, an alloy with 10 wt% Sn begins to solidify at \(T_L = 327 - 2.3 \times 10 = 304^\circ\)C. The first solid contains only \(k_p \times 10 = 3.1\) wt% Sn, so the remaining liquid is enriched in tin. As solidification proceeds, the liquid composition and the solid composition both increase, tracing out the liquidus and solidus curves respectively. The resulting microsegregation — the variation of composition within a single dendrite arm — is described by the Scheil equation (assuming no diffusion in the solid):
\[
C_s(f_s) = k_p C_0 (1 - f_s)^{k_p - 1},
\]
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.
Definition 4.2 (Worster's Mushy-Layer Equations). Let \(\phi(x,t)\) be the local solid fraction, \(T(x,t)\) the temperature, and \(C(x,t)\) the interstitial liquid composition. In the mushy zone, local thermodynamic equilibrium gives:
\[
T = T_L(C) = T_m - m_L C.
\]
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.
Theorem 4.1 (Linear Stability of the Mushy Layer — Worster). The quiescent base state of unidirectional solidification of a binary alloy is governed by a Rayleigh number:
\[
\mathrm{Ra}_m = \frac{\rho_l \beta_C \Delta C\, g\, \Pi_0 h}{\mu \kappa},
\]
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.
Proposition 4.1. The effective thermal conductivity \(k_{\text{eff}}\) of a mushy zone with solid fraction \(\phi\) satisfies the Hashin-Shtrikman bounds:
\[
k_s + \frac{(1-\phi)}{\frac{1}{k_l - k_s} + \frac{\phi}{n k_s}} \leq k_{\text{eff}} \leq k_l + \frac{\phi}{\frac{1}{k_s - k_l} + \frac{(1-\phi)}{n k_l}},
\]
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.
Example 4.3 (Rayleigh Number for NH\(_4\)Cl-H\(_2\)O). The aqueous ammonium chloride system (NH\(_4\)Cl-H\(_2\)O) is the canonical laboratory model for mushy-layer convection, since it is transparent and solidifies at convenient temperatures. For a 26 wt% NH\(_4\)Cl solution cooled from below at rate \(\Delta T/h \approx 1\) K/cm with mushy-layer thickness \(h \approx 2\) cm, permeability \(\Pi_0 \approx 10^{-9}\) m\(^2\), and dynamic viscosity \(\mu \approx 10^{-3}\) Pa\(\cdot\)s, the mushy-layer Rayleigh number is
\[
\mathrm{Ra}_m = \frac{\rho g \beta_C \Delta C\, \Pi_0 h}{\mu \kappa} \approx \frac{10^3 \times 10 \times 3 \times 10^{-2} \times 10^{-9} \times 2 \times 10^{-2}}{10^{-3} \times 10^{-7}} \approx 60.
\]
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.
Definition 4.3 (Eutectic Solidification). At the eutectic composition \(C_E\) and temperature \(T_E\), the liquid transforms simultaneously into two distinct solid phases \(\alpha\) and \(\beta\):
\[
\text{Liquid}(C_E) \xrightarrow{T_E} \alpha(C_\alpha) + \beta(C_\beta).
\]
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.
Example 4.1 (Al-Si Eutectic). The Al-Si system has a eutectic at \(C_E \approx 12.6\) wt% Si and \(T_E = 577^\circ\)C. During directional solidification at typical velocities of \(V \approx 10\;\mu\)m/s, the Jackson-Hunt theory predicts lamellar spacings of \(\lambda \approx 5{-}10\;\mu\)m, in good agreement with experiments. The Al-Si eutectic is widely used in automotive castings (engine blocks, pistons) where the fine eutectic microstructure provides good wear resistance.
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.
Definition 5.1 (Baiocchi Transform). Define the transformed variable:
\[
w(x,t) = \int_0^t T(x,\tau)\, d\tau.
\]
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.
Definition 5.2 (Classical Obstacle Problem). Given a domain \(\Omega \subset \mathbb{R}^n\), an obstacle \(\psi \in C^2(\Omega)\) with \(\psi < 0\) on \(\partial\Omega\), and boundary data \(g \geq \psi\) on \(\partial\Omega\), find \(u \in H^1(\Omega)\) minimizing:
\[
J[u] = \frac{1}{2}\int_\Omega |\nabla u|^2\, dx
\]
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\}\).
\[
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\).
\[
\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.
Example 5.1 (Elastic Membrane over an Obstacle). A stretched elastic membrane is pressed against a rigid obstacle from above. The displacement \(u(x)\) of the membrane minimizes its elastic energy (proportional to \(\int|\nabla u|^2\)) subject to staying above the obstacle \(u \geq \psi\). The contact region is the coincidence set, and the free boundary is the edge of contact. In one dimension with obstacle \(\psi(x) = 1 - x^2\) on \((-2,2)\) and \(u = 0\) on the boundary, the solution is:
\[
u(x) = \begin{cases} 1 - x^2 & |x| \leq 1, \\ 2 - 2|x| + x^2 - 1 & 1 < |x| < 2, \end{cases}
\]
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\}\)?
Theorem 5.1 (Optimal Regularity — Caffarelli). Let \(u\) solve the obstacle problem with \(\psi \in C^2(\Omega)\). Then:
\[
u \in C^{1,1}(\Omega),
\]
that is, \(u\) has Lipschitz continuous first derivatives. This regularity is optimal: \(u\) is generally not \(C^2\) across the free boundary.
Proof. (Sketch.) The key tool is the monotonicity formula. Define, for a free boundary point \(x_0 \in \Gamma\),
\[
\Phi(r) = \frac{1}{r^4}\int_{B_r(x_0)} |\nabla u|^2\, dx.
\]
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.
Theorem 5.2 (Free Boundary Regularity — Caffarelli). Let \(u\) solve the obstacle problem in \(\mathbb{R}^n\). Then the free boundary \(\Gamma = \partial\{u > \psi\}\) can be decomposed as \(\Gamma = \Gamma_{\text{reg}} \cup \Gamma_{\text{sing}}\), where:
(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.
Example 5.2 (Blow-Up Analysis in Two Dimensions). Consider the obstacle problem in \(\mathbb{R}^2\) with obstacle \(\psi \equiv 0\). At a regular free boundary point \(x_0 \in \Gamma\), define the rescaled functions \(u_r(x) = u(x_0 + rx)/r^2\). As \(r \to 0\), these converge (up to rotation) to the half-plane solution:
\[
u_0(x_1, x_2) = \frac{1}{2}\max(x_1, 0)^2.
\]
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.
Definition 5.3 (Thin Obstacle Problem). Let \(\Omega \subset \mathbb{R}^n\) be a domain bisected by a hyperplane \(\{x_n = 0\}\). Given a thin obstacle \(\psi\) on \(\{x_n = 0\}\), find \(u\) minimizing \(\int_\Omega |\nabla u|^2\, dx\) subject to \(u(x', 0) \geq \psi(x')\) on \(\{x_n = 0\} \cap \Omega\).
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.
Example 5.3 (The Signorini Problem in Elasticity). The original motivation for the thin obstacle problem comes from contact mechanics. Consider an elastic body \(\Omega\) resting on a rigid surface \(\{x_n = 0\}\) under the action of gravity. The displacement \(u\) minimizes the elastic energy subject to the non-penetration constraint \(u(x', 0) \geq 0\). The contact region \(\{u(x', 0) = 0\}\) and its boundary (the contact line) are unknown a priori. In the scalar simplification, this is exactly Definition 5.3 with \(\psi = 0\). The contact pressure \(-\partial u/\partial x_n|_{x_n = 0^+}\) is nonnegative in the contact region and vanishes outside it, giving the complementarity conditions:
\[
u \geq 0, \quad -\frac{\partial u}{\partial x_n} \geq 0, \quad u \cdot \frac{\partial u}{\partial x_n} = 0 \quad \text{on } \{x_n = 0\}.
\]
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.
Definition 6.1 (Ginzburg-Landau Free Energy). The Ginzburg-Landau (or Van der Waals-Cahn-Hilliard) free energy is:
\[
\mathcal{F}[\phi] = \int_\Omega \left(\frac{\epsilon^2}{2}|\nabla\phi|^2 + W(\phi)\right)\, dx,
\]
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\).
Definition 6.2 (Allen-Cahn Equation). The Allen-Cahn equation is the \(L^2\) gradient flow of \(\mathcal{F}\):
\[
\tau \frac{\partial\phi}{\partial t} = \epsilon^2 \nabla^2\phi - W'(\phi) = \epsilon^2 \nabla^2\phi + \phi - \phi^3,
\]
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.
Definition 6.3 (Cahn-Hilliard Equation). The Cahn-Hilliard equation for spinodal decomposition is:
\[
\frac{\partial\phi}{\partial t} = \nabla \cdot \left(M \nabla \mu\right), \qquad \mu = \frac{\delta\mathcal{F}}{\delta\phi} = -\epsilon^2\nabla^2\phi + W'(\phi),
\]
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}\).
Example 6.1 (Spinodal Decomposition in a Binary Alloy). Consider a homogeneous binary alloy with composition \(C_0 = 0.5\) (symmetric quench) inside the spinodal curve, where \(W''(C_0) < 0\). Linearizing the Cahn-Hilliard equation about the uniform state \(\phi = C_0\), a Fourier mode \(\hat{\phi}_k e^{ikx}\) grows at rate
\[
\sigma(k) = -M k^2\left(W''(C_0) + \epsilon^2 k^2\right).
\]
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.
Theorem 6.1 (Sharp-Interface Limit of Allen-Cahn). Consider the Allen-Cahn equation with a coupling to the temperature field:
\[
\tau\epsilon \frac{\partial\phi}{\partial t} = \epsilon^2\nabla^2\phi - W'(\phi) + \epsilon\, \lambda\, u,
\]
\[
\frac{\partial u}{\partial t} + \frac{L}{2c}\frac{\partial\phi}{\partial t} = \kappa\nabla^2 u,
\]
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.
Proof. (Sketch via matched asymptotics.) Divide the domain into an outer region (away from the interface, where \(\phi \approx \pm 1\)) and an inner region (within \(O(\epsilon)\) of the interface). In the outer region, expand \(\phi = \pm 1 + \epsilon\phi_1 + \cdots\) and \(u = u_0 + \epsilon u_1 + \cdots\). At leading order, \(W'(\pm 1)\phi_1 = \lambda u_0\), and \(u_0\) satisfies the heat equation. In the inner region, introduce the stretched coordinate \(\xi = d(x,t)/\epsilon\), where \(d\) is the signed distance to the interface. At leading order, \(\phi_0(\xi) = \tanh(\xi/\sqrt{2})\), the equilibrium profile. At next order, a solvability condition (Fredholm alternative applied to the linear operator \(\mathcal{L} = \partial_\xi^2 - W''(\phi_0)\)) yields the Gibbs-Thomson condition relating the interface temperature to the curvature and velocity. Matching the inner and outer temperature expansions produces the Stefan condition. \(\square\)
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).
Example 6.2 (Matched Asymptotics: Inner Region Solvability). The solvability condition in the matched-asymptotic derivation deserves a closer look, as it is the mechanism by which the Gibbs-Thomson condition emerges. In the inner region, at first order the perturbation \(\phi_1(\xi)\) satisfies the linear equation
\[
\mathcal{L}\phi_1 = \left(\frac{d^2}{d\xi^2} - W''(\phi_0(\xi))\right)\phi_1 = g(\xi),
\]
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.
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.
Definition 6.4 (Multi-Phase Field Model). For \(N\) phases with order parameters \(\phi_i(x,t)\), \(i = 1, \ldots, N\), subject to the constraint \(\sum_{i=1}^N \phi_i = 1\), the free energy is:
\[
\mathcal{F}[\{\phi_i\}] = \int_\Omega \sum_{i < j} \left(\frac{\epsilon_{ij}^2}{2}|\phi_i \nabla\phi_j - \phi_j \nabla\phi_i|^2 + W_{ij}(\phi_i, \phi_j)\right)\, dx,
\]
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
\[
\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.
Definition 7.1 (Level Set Equation). If the interface \(\Gamma(t)\) moves with normal velocity \(v_n\), the level set function evolves according to:
\[
\frac{\partial\phi}{\partial t} + v_n|\nabla\phi| = 0,
\]
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.
Theorem 7.1 (Viscosity Solution Framework — Crandall-Lions). The level set equation with velocity \(v_n = v_n(\kappa, x, t)\) admits a unique viscosity solution \(\phi \in C([0,T]; \text{Lip}(\mathbb{R}^n))\) provided \(v_n\) satisfies appropriate regularity and monotonicity conditions. The zero level set of this viscosity solution defines the generalized evolution of the interface, which agrees with the classical motion when the latter exists and provides a weak notion of motion past singularities.
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.
Example 7.2 (Narrow Band Level Set for Stefan Problem). For a two-dimensional Stefan problem on a \(256 \times 256\) grid, the full level set method requires solving the level set equation on all \(256^2 = 65{,}536\) grid points, even though only points near the interface contribute to the free boundary dynamics. The narrow band approach (Adalsteinsson and Sethian, 1995) restricts the level set computation to a band of width \(\sim 6\Delta x\) around the zero level set, reducing the cost from \(O(N^2)\) to \(O(N^{3/2})\) in two dimensions (or \(O(N^{5/3})\) vs. \(O(N^2)\) in three dimensions, where \(N\) is the number of grid points per dimension). The band is reconstructed periodically as the interface moves. Combined with the fast marching method (Sethian, 1996) for reinitialization, which solves the Eikonal equation \(|\nabla\phi| = 1\) in \(O(N^n \log N)\) time using an upwind Dijkstra-like algorithm, this makes level set methods competitive with front-tracking approaches for moderately complex interface geometries.
Example 7.1 (Motion by Mean Curvature). For an interface moving with normal velocity proportional to its curvature, \(v_n = -\kappa\), the level set equation becomes:
\[
\frac{\partial\phi}{\partial t} = |\nabla\phi|\nabla \cdot \left(\frac{\nabla\phi}{|\nabla\phi|}\right).
\]
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\).
Definition 7.2 (Reinitialization Equation). To restore the signed distance property, solve to steady state in fictitious time \(\tau\):
\[
\frac{\partial\phi}{\partial\tau} + \operatorname{sgn}(\phi_0)(|\nabla\phi| - 1) = 0,
\]
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.
Definition 7.3 (Volume-of-Fluid Method). The VOF method (Hirt and Nichols, 1981) represents the interface through the volume fraction \(f(x,t) \in [0,1]\), defined in each computational cell as the fraction occupied by one of the two fluids. The volume fraction satisfies:
\[
\frac{\partial f}{\partial t} + \nabla \cdot (f\mathbf{u}) = 0.
\]
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.
Example 7.3 (Mass Loss in Level Set vs. VOF for a Shrinking Bubble). Consider a circular gas bubble of initial radius \(R_0\) in a viscous liquid, shrinking under surface tension according to \(R(t) = R_0 - \sigma t/\mu\). On a \(128 \times 128\) grid with \(R_0 = 0.25\) (in a unit domain), a standard level set simulation accumulates a mass error of approximately 5--10% by the time the bubble has shrunk to half its initial radius. The VOF method preserves mass to machine precision by construction, since the volume fraction update is a conservative finite-volume scheme. The coupled level set/VOF (CLSVOF) method of Sussman and Puckett (2000) combines the mass conservation of VOF with the geometric accuracy of the level set (for curvature and normal computation), achieving mass errors below 0.1% while maintaining second-order accurate curvature estimates.
Definition 7.4 (Front-Tracking Method). In the front-tracking method (Tryggvason, Bunner, Esmaeli, et al.), the interface \(\Gamma(t)\) is represented explicitly by a set of connected marker points (in 2D) or triangulated surface elements (in 3D). The markers are advected with the local fluid velocity, and the interface forces (surface tension) are distributed to the underlying Eulerian grid using a regularized delta function.
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.
Definition 7.5 (Ghost Fluid Method). In the GFM for the Stefan problem, the temperature field is extended across the interface using "ghost" values that enforce the jump conditions. At grid points adjacent to the interface (identified by sign changes in the level set function), the finite difference stencil is modified so that the Stefan condition
\[
L\rho\,v_n = [k\nabla T \cdot \hat{n}]
\]
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.
Example 8.1 (Czochralski Crystal Growth). In the Czochralski process, a seed crystal is dipped into a crucible of molten silicon and slowly withdrawn while rotating. The solidification front is a free boundary whose shape is determined by the balance of heat transport (conduction in the crystal, conduction and convection in the melt, radiation from free surfaces), surface energy, and the imposed pulling rate. The mathematical model involves coupled free boundary problems for:
(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.
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.
Example 8.2 (Snowflake as a Free Boundary Problem). The growth of an ice crystal from vapor is governed by:
\[
\nabla^2 c = 0 \quad \text{in the vapor (quasi-steady)},
\]
\[
c = c_{\text{eq}}(\kappa, T) \quad \text{on the crystal surface } \Gamma(t),
\]
\[
v_n = D\frac{\partial c}{\partial n}\bigg|_\Gamma \quad \text{(mass conservation)},
\]
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).
Example 8.3 (Sea Ice as a Mushy Layer). First-year sea ice grows from the ocean surface when the atmospheric temperature drops below the freezing point of seawater (\(\approx -1.8^\circ\)C for 35 psu salinity). The mathematical model includes:
(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.
Example 8.4 (Continuous Casting of Steel). In a continuous caster operating at withdrawal speed \(V\), the shell thickness \(s(z)\) at distance \(z\) below the meniscus satisfies approximately
\[
s(z) \approx K\sqrt{z/V},
\]
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.
Example 8.6 (Melt Pool Geometry in Selective Laser Melting). For a Ti-6Al-4V powder bed irradiated by a laser of power \(P = 200\) W and scan speed \(V_s = 1\) m/s with spot radius \(r_0 = 50\;\mu\)m, the steady-state melt pool can be approximated by the Rosenthal solution — the three-dimensional analogue of the Neumann similarity solution for a moving point source. The temperature field in a frame moving with the laser is:
\[
T(x,y,z) = T_0 + \frac{P}{2\pi k R}\exp\left(-\frac{V_s(R + x)}{2\kappa}\right), \quad R = \sqrt{x^2 + y^2 + z^2},
\]
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).
Example 8.5 (Tumor Growth Model). A simplified model of avascular tumor growth treats the tumor as a region \(\Omega(t) \subset \mathbb{R}^3\) with a free boundary \(\Gamma(t)\). The nutrient concentration \(c(x,t)\) (e.g., oxygen) satisfies:
\[
\nabla^2 c = \lambda c \quad \text{in } \Omega(t) \quad \text{(consumption)},
\]
\[
c = c_\infty \quad \text{on } \Gamma(t) \quad \text{(supply from vasculature)},
\]
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.
Example 8.7 (Lava Flow as a Free Boundary Problem). When lava flows down a slope, it forms a crust by cooling from its exposed surface. The crust is a solid shell surrounding a liquid interior, and the growth of the crust is a Stefan problem coupled to the fluid dynamics of the underlying flow. For a basaltic lava with thermal diffusivity \(\kappa \approx 5 \times 10^{-7}\) m\(^2\)/s, latent heat \(L \approx 4 \times 10^5\) J/kg, and surface cooling rate sufficient to give an effective Stefan number \(\mathrm{St} \approx 0.5\), the crust thickness grows as \(h(t) \approx 2\lambda\sqrt{\kappa t}\), with \(\lambda \approx 0.5\). After 1 hour, this gives \(h \approx 0.04\) m, and after 1 day, \(h \approx 0.2\) m. The coupling between crust growth and lava rheology (the lava thickens as it cools, eventually stopping when the crust becomes thick enough to resist the driving gravitational stress) makes this a rich nonlinear free boundary problem with important implications for volcanic hazard assessment.
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.
