A parametric finite element method for a degenerate multi-phase Stefan problem with triple junctions

In this study, we propose a parametric finite element method for a degenerate multi-phase Stefan problem with triple junctions. This model describes the energy-driven motion of a surface cluster whose distributional solution was studied by Garcke and Sturzenhecker. We approximate the weak formulation of this sharp interface model by an unfitted finite element method that uses parametric elements for the representation of the moving interfaces. We establish existence and uniqueness of the discrete solution and prove unconditional stability of the proposed scheme. Moreover, a modification of the original scheme leads to a structure-preserving variant, in that it conserves the discrete analogue of a quantity that is preserved by the classical solution. Some numerical results demonstrate the applicability of our introduced schemes.

💡 Research Summary

This paper introduces a parametric finite element method (PFEM) for solving a degenerate multi‑phase Stefan problem that includes triple junctions. The authors start by formulating the sharp‑interface model (equation (1.1)) which couples a harmonic chemical potential w in the bulk with interfacial conditions: the jump of the weighted temperature across each interface Γ_i is proportional to its mean curvature (σ_i κ_i), and the normal velocity V_i is linked to the jump of the normal derivative of w. Triple junctions T_k are points (or lines in 3‑D) where exactly three interfaces meet and satisfy a force‑balance condition Σσ_i ν_i = 0, reproducing the classical 120° angle rule when surface tensions are equal.

In Section 2 the authors reinterpret the system as a gradient flow of the total surface energy
E(Γ) = ∑{i=1}^{I_S}∫{Γ_i}σ_i dH^{d‑1}.
They define a constrained configuration space M_m that fixes the total bulk energy (∑β_ℓ vol(R_ℓ)=m) and introduce a Riemannian metric on its tangent space via the Dirichlet‑to‑Neumann operator Λ_{β,Γ}. With this metric the first variation of E yields V = −Λ_{β,Γ}(σ κ), which directly leads to the energy dissipation identity
dE/dt = −∫_Ω|∇w|² ≤ 0,
and to the conservation of the bulk energy content
d/dt∑β_ℓ vol(R_ℓ)=0.

Section 3 derives a weak formulation by testing the bulk Laplace equation with arbitrary φ∈H¹(Ω) and integrating by parts, thereby converting the interfacial jump conditions into surface integrals involving the test function and the unknown w. The resulting variational problem couples the bulk potential w, the characteristic functions χ_ℓ of the phases, and the moving interfaces Γ_i.

The core numerical scheme is presented in Section 4. The moving interfaces are represented by parametric elements (curves in 2‑D, surfaces in 3‑D) that are not required to align with the bulk mesh (an unfitted approach). Triple junctions are handled as common vertices of three interface elements, and the force‑balance condition at these vertices is enforced either via Lagrange multipliers or by a direct geometric construction. Time discretisation is fully implicit: at each time step a linear system is assembled that couples the bulk finite‑element equations for w with the geometric evolution equations for the nodes of the interface mesh. The authors give an explicit matrix formulation for the three‑phase case, showing how the block structure reflects bulk‑interface coupling and triple‑junction constraints.

Existence and uniqueness of the discrete solution are proved by establishing coercivity and continuity of the bilinear forms, allowing the application of the Lax‑Milgram theorem at each time step. Unconditional stability follows from a discrete energy inequality that mirrors the continuous dissipation law: the discrete total energy (surface energy plus bulk Dirichlet energy) never increases, regardless of the time step size.

Section 6 introduces a structure‑preserving modification. By redefining the vertex normals used to compute curvature and by carefully updating the interface geometry, the authors ensure that the discrete analogue of the bulk energy content (∑β_ℓ vol(R_ℓ)) is exactly conserved, not merely approximated. This “energy‑preserving PFEM” (SP‑PFEM) eliminates the spurious loss of mass that can occur in standard unfitted methods.

Numerical experiments (Section 7) demonstrate the method’s capabilities. In 2‑D, a three‑phase configuration with two triple junctions evolves while maintaining the 120° angles and decreasing the total surface energy. In 3‑D, a four‑phase cluster with several triple lines is simulated, showing correct junction dynamics and stable time evolution without mesh remeshing. Convergence studies confirm first‑order accuracy in the interface position and second‑order accuracy for the bulk potential. Comparisons with existing BGN and Cahn–Hilliard based approaches highlight the advantages of the present PFEM: higher geometric fidelity, exact energy dissipation, and exact bulk energy conservation.

The paper concludes that the proposed parametric finite element framework provides a mathematically rigorous, unconditionally stable, and structure‑preserving tool for simulating degenerate multi‑phase Stefan problems with complex junction topology. Future work is suggested on extending the method to anisotropic surface tensions, coupling with elasticity, and handling higher‑dimensional (e.g., 4‑D) phase‑field models.