🔬 physics

Self-similar structure of non-isothermal variable-density mushy Stefan problems and an improved low-Mach enthalpy method

This paper analyzes the self-similar structure of finite-width mushy Stefan problems to improve the accuracy and stability of the low-Mach enthalpy method for simulating non-isothermal phase-change processes with variable densities.

Original authors: Yavkreet Swami, Amneet Pal Singh Bhalla

Published 2026-08-04

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Yavkreet Swami, Amneet Pal Singh Bhalla

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the world of materials science as a giant, bustling kitchen where scientists are trying to master the art of cooking with extreme temperatures. Sometimes, they need to melt metal to cast a car part, freeze water to make ice cream, or even boil liquids to create steam. The tricky part isn't just heating things up; it's understanding what happens right at the edge where a solid turns into a liquid, or vice versa. This boundary is called a "phase-change front." In the real world, this isn't a razor-sharp line like a knife cutting paper. Instead, it's more like a fuzzy, melting zone where solid and liquid mix together, kind of like a slushy drink that is half ice and half water.

For decades, scientists have used a clever mathematical trick called the "enthalpy method" to simulate these melting and freezing processes on computers. Think of this method as a way to track the "heat energy" of a material without having to constantly redraw the map of where the ice ends and the water begins. It's a favorite tool for engineers designing everything from airplane engines to 3D printers. However, there's a catch: most of these computer models assume that the solid and the liquid weigh the same amount. But in reality, when things melt or freeze, they often shrink or expand, changing their density. This is like a crowd of people suddenly deciding to hold hands and squeeze together, or spread out and dance apart. When this happens, the material starts to flow, creating currents that can mess up the melting process. Until now, it wasn't clear if the old, simple computer models could handle these "heavy" changes accurately, especially when the density difference is huge, like between a liquid and a gas.

This paper dives deep into that fuzzy, melting zone to see if the old rules still hold up when things get heavy and flowy. The authors, Yavkreet Swami and Amneet Pal Singh Bhalla, decided to treat the melting zone not as a mistake to be ignored, but as a real, physical place with its own rules. They discovered that this "mushy" region has a hidden, self-similar structure—meaning if you zoom in or out, the pattern of how it grows looks the same, just scaled up or down, like a fractal snowflake. By using this pattern, they were able to solve the math for how the temperature and the flow of the material behave when the solid and liquid have very different weights.

Their big finding is that the old way of thinking about these problems was slightly off. They showed that the computer models (enthalpy methods) are actually solving a problem with a "fuzzy" boundary, not the sharp, thin line that classical physics often assumes. For materials that don't change weight much (like most metals), the fuzzy models work great and match the sharp-line predictions. But for materials that change weight a lot—like when a liquid turns into a gas, or a solid turns into a liquid with a massive density shift—the old models start to wobble and become unstable. The authors proved that to get the right answer for these tricky cases, you have to make the "fuzzy" zone in the computer model very, very thin, or else the results will drift away from reality.

To fix this, they tweaked their computer code, creating an "improved low-Mach enthalpy method." Think of this as upgrading a recipe to handle a new, difficult ingredient. They changed how they calculate the energy inside that mushy zone, making sure the math respects the fact that the material is flowing because it's changing weight. When they tested this new method against their own new mathematical "gold standard" (the similarity solution they derived), it worked perfectly. It stayed stable and accurate even when the density difference was huge, reaching ratios of 540 to 1 (like comparing the weight of a feather to a heavy rock).

The paper also ran a series of tests to see how fast and accurate the new method is. They found that as they made their computer grid finer (like using a higher resolution camera), the results got better and better, converging at a rate between first and second order. This means the new method is reliable and precise. The authors are confident that this work opens the door for using these powerful simulation tools on much more extreme problems, like boiling and condensation, where the density changes are massive. They didn't just find a bug; they rewrote the rulebook for how to simulate melting and freezing when the material decides to change its weight and start flowing.

Technical Summary: Self-similar structure of non-isothermal variable-density mushy Stefan problems and an improved low-Mach enthalpy method

Problem Statement
The enthalpy method is a widely adopted approach for simulating phase-change problems (melting and solidification) on fixed grids without explicit interface tracking. It typically models the phase transition over a finite temperature interval, $\Delta T = T_{liq} - T_{sol}$ , creating a "mushy" region. While the method is generally assumed to converge to the classical sharp-interface Stefan problem as $\Delta T \to 0$ , this assumption has largely been validated under the simplifying condition of equal solid and liquid densities ( $\rho_S = \rho_L$ ).

In practical applications involving significant volume changes (e.g., boiling, condensation, or specific metal solidification), the density ratio $\rho_L/\rho_S$ can deviate substantially from unity. In such cases, phase change induces fluid motion due to density variations, coupling the thermal problem with mass and momentum conservation. Existing literature lacks a systematic investigation into the accuracy and convergence properties of enthalpy methods for high-density-ratio problems with finite $\Delta T$ . Furthermore, previous validations against the sharp-interface limit may be insufficient for assessing methods designed to handle density-change-induced flow, as the mushy Stefan problem (finite $\Delta T$ ) is fundamentally distinct from the classical Stefan problem (infinitesimal $\Delta T$ ).

Methodology
The authors address these gaps through a combination of theoretical analysis and numerical development:

Theoretical Formulation of the Mushy Stefan Problem:
The authors formulate a non-isothermal, variable-density mushy Stefan problem. They demonstrate that the governing equations possess a self-similar structure, reducing the partial differential equations to a boundary-value problem involving two unknown similarity parameters ( $\lambda_s$ and $\lambda_l$ ) that define the solidus and liquidus front positions.
- Equal-Density Case: A closed-form analytical similarity solution is derived.
- Unequal-Density Case: The problem reduces to a coupled nonlinear system governing temperature and velocity fields in the mushy and liquid regions. This system is solved numerically using a boundary-value solver (MATLAB's bvp5c) with a continuation strategy to handle extreme density ratios (up to $\rho_L/\rho_S = 2700$ ) and small $\Delta T$ .
Analysis of Enthalpy Definitions:
The study examines two definitions for specific enthalpy in the mushy region. The "original" definition leads to a nonlinear relationship between liquid fraction ( $\phi$ ) and temperature ( $T$ ) when densities differ, while the "improved" definition yields a linear $\phi$ - $T$ relation. The authors show that the original formulation exhibits spurious instabilities at high density ratios, whereas the improved formulation remains stable.
Improved Low-Mach Enthalpy Method:
Based on the theoretical insights, the authors propose an improved low-Mach enthalpy method. Key features include:
- A revised specific enthalpy definition that ensures a linear liquid-fraction–temperature relation.
- A kinematic velocity-divergence constraint derived directly from mass conservation to account for density-change-induced flow.
- Implementation within the IBAMR library using a staggered Cartesian grid, Newton iterations for the nonlinear enthalpy equation, and a projection method with a preconditioner to handle the stiff Carman–Kozeny drag force in the solid region.

Key Results

Self-Similarity and Convergence: The analysis confirms that the enthalpy method solves the mushy Stefan problem, not the classical Stefan problem. The width of the mushy region grows proportionally to $\sqrt{t}$ , whereas the classical Stefan front remains infinitesimally thin.
Density Ratio Sensitivity:
- For equal and low-density ratios ( $\rho_L/\rho_S \approx 1$ to $2$), the mushy Stefan solution closely approximates the sharp-interface solution over a wide range of $\Delta T$ values. Consequently, equal-density benchmarks are not stringent enough to validate methods for density-change-induced flow.
- For large-density ratios ( $\rho_L/\rho_S \gg 1$ ), the phase-front location becomes highly sensitive to $\Delta T$ . The mushy solution only approximates the sharp-interface limit when $\Delta T$ is sufficiently small.
Method Performance:
- The improved low-Mach enthalpy method demonstrates excellent agreement with the derived mushy Stefan similarity solutions across a wide range of $\Delta T$ values and density ratios (up to $\rho_L/\rho_S \approx 540$ in verification cases).
- In contrast, the original formulation produces interfacial instabilities and unphysical vortical structures in the liquid region for large density ratios.
- Convergence Rates: The improved method achieves spatial-temporal convergence rates between first and second order for both phase-front location and temperature fields, even at density ratios of $O(3)$ .

Significance and Claims
The paper claims to fill a longstanding gap in the enthalpy method literature by providing a rigorous theoretical benchmark (the similarity solution of the non-isothermal variable-density mushy Stefan problem) for validating numerical methods at finite $\Delta T$ .

The authors assert that:

Benchmark Validity: Large-density-ratio cases provide a significantly more stringent validation benchmark for modern enthalpy methods than equal-density cases, as the latter fail to capture the sensitivity of the solution to the phase-change interval.
Methodological Improvement: The proposed improvements to the low-Mach enthalpy method (specifically the revised enthalpy definition and linear $\phi$ - $T$ relation) are essential for stability and accuracy when simulating phase change with substantial density contrasts.
Extended Applicability: These theoretical and numerical developments extend the applicability of enthalpy methods beyond traditional metal solidification (where density ratios are near unity) to processes involving much larger density ratios, such as boiling and condensation.

The work does not propose new experimental applications but rather establishes a robust numerical framework and theoretical foundation for extending enthalpy-based simulations to a broader class of variable-density phase-change problems.

More like this