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

TL;DR

This work extends the Stefan problem to a three-phase spherical geometry with finite-sized particles, incorporating density and kinetic-energy jumps across moving solid–liquid and liquid–vapor interfaces and curvature effects via Gibbs–Thomson relations. A fixed-grid sharp-interface method in spherical coordinates is developed and shown to achieve second-order accuracy, with small-time analytical solutions used to initialize simulations in both two- and three-phase settings. Validations against established two-phase results and detailed analyses of kinetic-energy contributions reveal that kinetic terms delay melting for nano- to sub-micron particles, while their impact diminishes for micron-scale particles, highlighting scale-dependent physics relevant to metal additive manufacturing. The approach provides a robust, open-source framework for fully-resolved CFD studies of simultaneous melting, boiling, and evaporation in spherical phase-change systems.

Abstract

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).

Figures (8)

• Figure 1: Schematic of the \ref{['fig_PCMGraphic_A']} two-phase and \ref{['fig_PCMGraphic_B']} three-phase spherical Stefan problem in which an initially solid PCM melts and boils, respectively due to an imposed temperature condition at the outer surface ($r = R_b(t)$).
• Figure 2: Schematic of the 1D grid used to discretize the spherical heat equations using the sharp interface technique. Irregular cells $j$ and $j+1$, $p$ and $p+1$, and $N$ abut interfaces separating phases $A$ and $B$, and phases $B$ and $C$, and outer free boundary, respectively. The cell centers are marked with $(\bullet)$, and the evolving interfaces and free boundary are marked with dashed lines (---).
• Figure 3: Melting of a gold nanoparticle of initial radius $R_0 = 100$ nm for two Stefan numbers: $\beta_m = 100$ and $\beta_m = 10$. $\rho^{\rm SL} = 1$ implies gold solid and liquid phases have the same density, and $\rho^{\rm SL} \ne1$ implies different densities. \ref{['fig_CompBeta100']} Time evolution of the melt front for Stefan number $\beta_m = 100$ and \ref{['fig_CompBeta10']} for Stefan number $\beta_m = 10$. \ref{['fig_CompTemp']} Dimensional temperature variation within the nanoparticle at two different physical times for Stefan number $\beta_m = 100$. Results are compared against those reported in Font et al. font2015nanoparticle. The dotted line (..) shows the melt temperature variation.
• Figure 4: Convergence rate of the numerical error for \ref{['fig_IntfConv3P']} interface position $\widehat{R}_1$, \ref{['fig_IntfR2Conv3P']} interface position $\widehat{R}_2$ and \ref{['fig_TempConv3P']} temperature distribution as a function of grid size $N$. For interface positions, errors are computed over the time interval $\widehat{t} = 0.1911$ to $\widehat{t} = 0.1930$. For temperature, the error is computed at $\widehat{t} = 0.19$.
• Figure 5: Time evolution of interface positions and temperature distribution within the domain for a 100 nm Al-like particle, with and without considering the kinetic energy terms in the Stefan conditions. The vapor phase density is taken to be $\rho^{\rm V} = 500~\mathrm{kg/m^3}$. Time evolution of \ref{['fig_R1']} the melt front, \ref{['fig_R2']} the boiling front, and the \ref{['fig_Rb']} outer free boundary. \ref{['fig_TempComparison3P']} Temperature distribution within the particle at $\widehat{t} = 0.1$.