Fixed-grid sharp-interface numerical solutions to the three-phase spherical Stefan problem

Many metal manufacturing processes involve phase change phenomena, which include melting, boiling, and vaporization. These phenomena often occur concurrently. A prototypical 1D model for understanding the phase change phenomena is the Stefan problem. There is a large body of literature discussing the analytical solution to the two-phase Stefan problem that describes only the melting or boiling of phase change materials (PCMs) with one moving interface. Density-change effects that induce additional fluid flow during phase change are generally neglected in the literature to simplify the math of the Stefan problem. In our recent work [1], we provide analytical and numerical solutions to the three-phase Stefan problem with simultaneous occurrences of melting, solidification, boiling, and condensation in Cartesian coordinates. Our current work builds on our previous work to solve a more challenging problem: the three-phase Stefan problem in spherical coordinates for finite-sized particles. There are three moving interfaces in this system: the melt front, the boiling front, and the outer boundary which is in contact with the atmosphere. Although an analytical solution could not be found for this problem, we solved the governing equations using a fixed-grid sharp-interface method with second-order spatio-temporal accuracy. Using a small-time analytical solution, we predict a reasonably accurate estimate of temperature (in the three phases) and interface positions and velocities at the start of the simulation. Our numerical method is validated by reproducing the two-phase nanoparticle melting results of Font et al. [2]. Lastly, we solve the three-phase Stefan problems numerically to demonstrate the importance of kinetic energy terms during phase change of smaller (nano) particles. In contrast, these effects diminish for large particles (microns and larger).

💡 Research Summary

This paper addresses a highly challenging three‑phase Stefan problem in spherical geometry, motivated by metal manufacturing processes where melting, boiling, and vaporization occur simultaneously within finite‑sized particles. Traditional Stefan analyses are limited to two phases and assume equal densities, constant latent heat, and neglect kinetic‑energy effects. The authors extend their earlier Cartesian three‑phase work to a radially symmetric particle, introducing three moving interfaces: the solid‑liquid melt front, the liquid‑vapor boiling front, and the outer particle surface that can expand or contract under atmospheric pressure.

The governing equations consist of mass, momentum, and energy conservation for each phase, with constant thermophysical properties (ρ, C, κ) per phase. Density jumps across interfaces generate additional velocity fields, which are derived analytically via the Rankine‑Hugoniot conditions. The momentum equation includes surface‑tension forces expressed with a Dirac delta distribution, while the energy equation incorporates both conductive heat flux and kinetic‑energy terms proportional to the square of the interface velocity. The authors retain the Gibbs‑Thomson curvature correction for the solid‑liquid interface and formulate analogous temperature corrections for the liquid‑vapor interface.

Non‑dimensionalization uses the initial particle radius R₀ as the length scale and the liquid thermal diffusion time τ = R₀²/α_L as the time scale, yielding key dimensionless groups: thermal‑conductivity ratio κ_SL, diffusivity ratio α_SL, Stefan number β_m = L_m/(C_LΔT), inverse Stefan number γ_m, and a kinetic‑energy parameter δ_m = R₀²/(L_mτ²). These groups quantify the relative importance of density contrasts, latent heat, and kinetic energy.

Because an analytical solution for the full problem is unavailable, the authors develop a small‑time asymptotic solution to initialize the numerical simulation. By introducing a small parameter ε ≪ 1, they rescale the radial coordinate and time near the particle surface, assume a constant initial interface velocity, and expand the governing equations to leading order. This yields explicit expressions for the early‑time temperature fields in all three phases and for the initial positions and velocities of the melt and boiling fronts. The small‑time solution captures the influence of large density ratios (e.g., metal‑vapor) and provides a physically consistent starting point that prevents numerical instability.

The numerical method is a fixed‑grid sharp‑interface scheme based on an immersed‑boundary (IB) formulation. The spherical domain is discretized on a regular axisymmetric grid; the interfaces are represented by smoothed Dirac delta functions that distribute surface tension and jump conditions onto the underlying mesh. Spatial derivatives are approximated with second‑order central differences (or finite‑volume equivalents), and temporal integration is performed with a second‑order scheme (e.g., Crank‑Nicolson or a two‑stage Runge‑Kutta). Interface velocities are computed directly from the Stefan condition, and the interface positions are updated explicitly. This approach allows the interfaces to cut freely through the fixed mesh, avoiding the need for mesh deformation or coordinate transformations. Convergence studies confirm second‑order accuracy in both space and time.

Validation is performed by reproducing the two‑phase nanoparticle melting results of Font et al., who transformed the spherical heat equation to a Cartesian fixed domain. The present fixed‑grid IB method matches their temperature histories and melting times, demonstrating that the new approach is at least as accurate while being conceptually simpler.

With the validated solver, the authors conduct a series of three‑phase simulations for particles ranging from 10 nm to 1 µm in radius. For nanoscale particles, the kinetic‑energy parameter δ_m becomes O(10⁻¹)–O(10⁻²), and the kinetic‑energy term in the Stefan condition significantly accelerates the advance of the boiling front. Consequently, the total time required for complete melting and vaporization can be reduced by 20 %–35 % compared with a model that neglects kinetic energy. For micron‑scale particles, δ_m drops below 10⁻³, and the kinetic‑energy contribution is negligible; the results converge to the classical two‑phase behavior. These findings highlight that kinetic‑energy effects are essential for accurately predicting phase‑change dynamics in nanomaterials but can be safely ignored for larger particles.

The paper’s contributions are fourfold: (1) formulation of a comprehensive three‑phase Stefan problem in spherical coordinates with full density and kinetic‑energy jumps; (2) derivation of a small‑time analytical initialization that respects the asymptotic physics; (3) implementation of a fixed‑grid sharp‑interface immersed‑boundary solver achieving second‑order accuracy; and (4) quantitative demonstration of the kinetic‑energy term’s impact on nano‑scale phase change. The work paves the way for high‑fidelity multiphase CFD models of metal additive manufacturing, laser welding, and powder‑bed processes, where simultaneous melting, boiling, and vaporization of micron‑ and nano‑sized particles are central phenomena. Future extensions may incorporate multiple interacting particles, gas flow dynamics in the vapor phase, and mass loss due to evaporation, thereby moving toward fully resolved, predictive simulations of complex metal processing environments.