🔢 mathematics

Finite element approximation of the enthalpy formulation for Stefan problems on evolving surfaces

This paper proposes and analyzes a spatially discrete evolving surface finite element method for the enthalpy formulation of two-phase Stefan problems on evolving surfaces, establishing $\mathcal{O}(\sqrt{h})$ error bounds without mass-lumping and introducing a novel implementation that avoids numerical quadrature errors.

Original authors: Philip J. Herbert, Thomas Sales, Chandrasekhar Venkataraman

Published 2026-07-07

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Philip J. Herbert, Thomas Sales, Chandrasekhar Venkataraman

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 you have a soap bubble floating in the air. Now, imagine that bubble is freezing. As it freezes, the ice forms a boundary between the liquid part and the solid part. This boundary isn't static; it moves, shifts, and changes shape as the bubble itself might be expanding, shrinking, or wobbling in the wind.

This paper is about creating a computer simulation to predict exactly how that freezing process happens on a surface that is constantly changing shape.

Here is a breakdown of what the authors did, using simple analogies:

1. The Problem: The "Melting Ice Cube" on a Wobbly Trampoline

The scientific problem they are solving is called the Stefan problem. In everyday terms, think of an ice cube melting. You have two phases: solid ice and liquid water. The tricky part is the "free boundary"—the exact line where the ice turns to water.

Usually, scientists simulate this on a flat, stationary table. But this paper tackles a much harder version: What if the table itself is a wobbly, stretching trampoline?

The surface (the bubble or the trampoline) is moving.
The material on it (the heat/enthalpy) is changing.
The boundary between the "frozen" and "thawed" parts is moving along with the surface.

2. The Tool: A Digital Mesh That Moves

To solve this, the authors used a method called Finite Element Analysis.

The Analogy: Imagine trying to draw a map of a moving, stretching rubber sheet. You can't just use a fixed grid of graph paper because the paper would tear or crumple.
The Solution: Instead, you place a net of tiny triangles (a mesh) directly onto the rubber sheet. As the sheet stretches and moves, the knots in your net move with it. The computer calculates the physics on these moving triangles. This is called an Evolving Surface Finite Element Method (ESFEM).

3. The Innovation: A New "Translator" for the Math

The authors faced a specific mathematical headache. When the surface moves, the math gets very messy, especially because the "freezing" part of the equation is "degenerate" (it behaves strangely when things are exactly at the freezing point).

The Old Way: Previous methods often relied on "mass-lumping," which is like simplifying a complex recipe by ignoring some ingredients to make the math easier. But on a moving surface, this simplification can break the rules of physics (like making heat flow backward).
The New Way: The authors invented a new mathematical "translator" (a projection-type operator).
- The Analogy: Imagine you have a message written on a stretching rubber band (the computer mesh). You need to read it on a rigid wall (the real surface). Usually, stretching the rubber band distorts the message. The authors built a special lens that looks at the message on the rubber band and translates it perfectly onto the wall, correcting for the stretch without losing any information.
- Why it matters: This allowed them to prove their method is stable and accurate without needing to simplify the physics or assume the solution is perfectly smooth (which it rarely is in real life).

4. The "Exact" Trick: No More Guessing

When computers do math, they often have to guess the area under a curve (a process called numerical quadrature). It's like estimating the amount of paint needed for a wall by counting tiles rather than measuring the wall.

The Innovation: The authors proposed a new way to implement the code that avoids these "guesses" entirely for certain types of problems.
The Analogy: Instead of estimating how much paint you need by looking at the tiles, they found a way to calculate the exact amount of paint needed for every single shape on the wall, even if the wall is moving. They call this an "exact discretisation." This removes a whole layer of potential computer error.

5. The Surprising Discoveries

The paper includes computer experiments that showed some surprising things that only happen on moving surfaces:

Spontaneous Freezing/Thawing: On a stationary table, if you start with a specific temperature, the laws of physics say the temperature can't suddenly jump to a new state without an outside heat source. However, on a moving surface, the motion itself acts like a heat source.
- The Analogy: Imagine a crowd of people on a dance floor. If the floor suddenly shrinks, the people get squeezed together (density increases). If the floor expands, they spread out. The authors found that just by shrinking or expanding the surface, the "temperature" of the material could spontaneously change, causing new phases (like ice forming in liquid water) to appear out of nowhere, even if no external heat was added.
Mushy Regions: They also observed "mushy regions" (areas that are half-frozen, half-liquid) forming spontaneously on moving surfaces, something that usually requires external heat sources to happen on stationary ones.

6. The Result: How Accurate is it?

The authors proved mathematically that their method is accurate.

They showed that as they made the "triangles" in their net smaller (increasing the resolution), the error in their simulation decreased.
They proved the error shrinks at a rate of roughly the square root of the mesh size ( $O(\sqrt{h})$ ).
The Happy Surprise: In their computer tests, the method actually performed better than their math predicted, showing a much faster rate of accuracy ( $O(h^2)$ ) for the specific smooth examples they tested.

Summary

In short, this paper builds a robust, high-precision computer engine to simulate phase changes (like freezing and melting) on surfaces that are moving and changing shape. They did this by inventing a new mathematical "lens" to handle the moving geometry and a new coding trick to eliminate calculation errors. They proved the engine works and discovered that moving surfaces can cause freezing and thawing to happen in ways that stationary surfaces never could.

Technical Summary: Finite Element Approximation of the Enthalpy Formulation for Stefan Problems on Evolving Surfaces

Problem Statement
This paper addresses the numerical approximation of the two-phase Stefan problem posed on an evolving surface, $\Omega(t) \subset \mathbb{R}^3$ . The problem models phase transitions (e.g., freezing soap bubbles, aircraft icing, cell polarization) where the surface moves with a prescribed velocity field $V$ . The authors focus on the enthalpy formulation, which avoids explicit tracking of the free boundary $\Gamma(t)$ by treating the enthalpy $e$ and temperature $u$ as a coupled system:
$\partial_\bullet e + e(\nabla_\Omega \cdot V) - \Delta_\Omega u = f, \quad e \in \beta(u),$
where $\partial_\bullet$ is the material derivative, $\nabla_\Omega$ and $\Delta_\Omega$ are surface operators, and $\beta$ is a maximal monotone set-valued map representing the phase change (e.g., latent heat). The primary challenge lies in the degeneracy of the resulting parabolic equation and the lack of regularity in solutions, particularly when "mushy regions" (where the free boundary has an interior) form. To date, no literature has provided a numerical approximation for Stefan-type problems on evolving surfaces.

Methodology
The authors propose a spatially discrete Evolving Surface Finite Element Method (ESFEM) based on the work of Dziuk and Elliott. The methodology involves:

Spatial Discretization: The evolving surface is approximated by a triangulated surface $\Omega_h(t)$ with piecewise linear finite element spaces $S_h(t)$ . The method utilizes "lifts" to map functions between the discrete and continuous surfaces to handle the variational crime of non-conformity.
Enthalpy Formulation: The authors avoid mass-lumping techniques and discrete maximum principles, which are often required in stationary settings but are problematic on evolving surfaces due to mesh deformation. Instead, they rely on the monotonicity of the enthalpy function $\beta$ .
New Projection Operator: A central methodological innovation is the introduction of a new projection-type operator, $P_h$ , mapping discrete functions to the continuous space. This operator is designed to preserve mean values and commute effectively with the material derivative, allowing the authors to bypass the strong regularity assumptions (e.g., $L^2H^2$ ) typically required in ESFEM error analysis.
Implementation Strategies: The paper discusses two fully discrete schemes (using backward Euler in time):
- A standard approach using numerical quadrature for the nonlinear term.
- A novel "exact discretization" approach for piecewise polynomial enthalpy functions. This method avoids quadrature errors by analytically integrating the nonlinear term over sub-regions of the mesh elements, a technique not previously considered even in stationary settings.

Key Contributions

First Error Bounds on Evolving Surfaces: The paper provides the first error analysis for a free boundary problem posed on an evolving surface. The main result (Theorem 5.18) establishes an $O(\sqrt{h})$ error bound for the temperature in the $L^2L^2$ norm under minimal regularity assumptions.
Novel Projection Operator: The introduction of the operator $P_h$ (Definition 5.11) enables the analysis of degenerate problems without requiring high regularity of the time derivative of the solution, a significant departure from standard ESFEM analysis.
Quadrature-Free Implementation: The authors propose and analyze a numerical implementation that avoids errors due to numerical quadrature for piecewise polynomial enthalpy functions. This method is applicable to both evolving and stationary flat domains.
Stability and Continuous Dependence: The paper proves the stability of the semi-discrete scheme and establishes continuous dependence on the data (Lemma 5.5), ensuring the numerical solution behaves consistently with perturbations in initial data and source terms.

Results and Numerical Experiments

Theoretical Bounds: The derived error bound is $O(\sqrt{h})$ in the $L^2L^2$ norm for temperature, matching the order of convergence found in stationary, flat settings (e.g., [24, 52]).
Experimental Order of Convergence (EOC): Numerical experiments on expanding and rotating spheres demonstrate that the method achieves an observed convergence rate of $O(h^2)$ for the temperature in the $L^2L^2$ norm. The authors attribute this super-convergence to the smoothness of the specific exact solutions constructed for the tests, noting that constructing rougher solutions for surface Stefan problems is difficult.
Phenomena on Evolving Domains: The experiments illustrate unique behaviors arising only on evolving domains:
- Nucleation without Heat Sources: Surface evolution (expansion or contraction) can force the enthalpy to cross phase thresholds, causing new phases to nucleate even in the absence of external heat sources ( $f=0$ ).
- Mushy Region Formation: Contrary to stationary settings where mushy regions typically require heat sources, the numerical results show that mushy regions can spontaneously form on evolving surfaces even with $f=0$ , driven by the geometric evolution of the domain.

Significance and Claims
The authors position this work as an important step toward the numerical analysis of degenerate parabolic PDEs on evolving surfaces. They claim that their approach, which avoids mass-lumping and relies on a new projection operator, is robust and extends the applicability of ESFEM to singular and degenerate problems (such as the porous medium equation or Cahn-Hilliard equation) where regularity is limited. While the theoretical error bound is $O(\sqrt{h})$ , the paper notes that the method performs comparably to quadrature-based methods and may exhibit higher convergence rates in practice for smooth solutions. The work is presented as a foundational contribution, with the authors explicitly stating that the specific case of the Stefan problem is their primary focus, but the analytical tools developed are expected to have significant impact beyond this specific application.