1 Introduction↩︎

We are interested in the analysis of evolving surface finite element methods (ESFEM) for two-phase Stefan problems posed on an evolving surface, \(\Omega(t) \subset \mathbb{R}^3\), moving with a prescribed velocity, \(\mathbf{V}\). In particular, we consider the so-called enthalpy formulation of the Stefan problem, \[\tag{1} \begin{gather} \partial^\bullet e + e (\nabla_{\Omega}\cdot \mathbf{V}) - \Delta_{\Omega}u = f, \text{ on } \Omega(t), \tag{2}\\ e \in \beta(u), \tag{3} \end{gather}\] with initial condition \(e(0) = e_0 \in L^2(\Omega(0))\) and where the enthalpy \(\beta\) is a set-valued map, as discussed in Section 2. Here the function \(e\) is the enthalpy function (i.e. the heat content) and \(u\) is the temperature distribution. We shall assume that \(\Omega(t)\) is a sufficiently smooth surface without boundary, and hence there are no boundary conditions associated with 1 . We defer discussion of the differential operators used in 1 until Section 2.

Problems of the form 1 were proposed and analysed in [1], wherein the authors study the well-posedness for suitably defined solutions. Stefan-type problems posed on a surface appear in nature as models of phase transitions on surfaces, for example when a soap bubble freezes [2], [3], in free boundary limits of bulk-surface models of ligand-receptor dynamics [4], [5], and in models of cell-polarisation [6], [7]. There is also interest in industrial applications of Stefan problems on surfaces, such as welding [8]–[10] and aircraft icing [11], [12], which can be modelled as Stefan problems on surfaces (which may also be deforming). This enthalpy formulation provides a generalised notion of solution for the Stefan problem and has been used extensively since initial work in the 1960s by Oleı̆nik [13], Kamenomostskaja [14], and Friedman [15]. This weak notion of solution allows one to implicitly capture the evolution of the free boundary, \(\Gamma(t)\) as illustrated in Figure 2, through the definition of the graph \(\beta\) (sometimes referred to as the generalised enthalpy [16]), instead of explicitly capturing the free boundary, as in front-tracking methods [17]. We note that sufficiently regular enthalpy solutions do indeed solve the more standard strong formulation of the Stefan problem 7 , provided the free boundary, \(\Gamma(t)\), does not develop an interior (a so-called mushy region), cf. [1]. We refer the reader to [16], [18]–[21] for further details on the enthalpy formulation of the Stefan problem, as well as [22]–[28] for results on finite element approximations of the two-phase Stefan problem. To our knowledge, there is no literature concerning the numerical approximation of Stefan-type problems posed on an evolving surface. We do however, refer the reader to recent work by Garcke and Nürnberg [3] where they study a free boundary problem governing anisotropic crystal growth on stationary surfaces. The corresponding system is a surface Stefan problem with surface tension and kinetic undercooling, and the authors propose a finite element approximation in the style of Barrett–Garcke–Nürnberg type discretisations, cf. [29].

The ESFEM developed by Dziuk and Elliott [30] is a robust and efficient numerical method for the solution of PDEs posed on evolving surfaces. The method has been applied to a number of different equations and its analysis is a burgeoning area of current numerical analysis research, see e.g., [31] for a review. Despite their relevance in analysis and applications [32] the numerical approximation of degenerate parabolic PDEs on evolving surfaces with the ESFEM or other methods has, to our best knowledge, not been considered previously. This work is therefore, an important step in this direction and is expected to have significant impact beyond the specific case of the Stefan problem on an evolving surface which is our primary focus, cf. Remark 8.
Main contributions:

• Our main contribution (Theorem 10) is an error bound for a semi-discrete ESFEM approximation of 1 . Our approach does not rely on mass-lumping techniques, and obtains the same order error as in the stationary, flat setting [23], [28].

• We introduce a new projection-type operator (Definition 5), which enables us to avoid the strong (and in this case infeasible) regularity assumptions often required in ESFEM error analysis [33]. We believe this approach will be useful in the analysis of ESFEM for other degenerate problems, cf. Remark 8.

• We introduce a new numerical method for solving the fully discrete problem by using an “exact discretisation” that avoids errors due to numerical quadrature under the assumption that \(\beta\) is piecewise polynomial. This numerical method is also applicable and novel in the stationary, flat setting.

Outline: The outline of this paper is as follows. In Section 2 we discuss the two-phase Stefan problem on an evolving surface and state our assumptions. Section 3 recalls relevant material on evolving surfaces, and the evolving function space theory developed in [34], [35]. Section 4 recalls aspects of the evolving surface finite element method (ESFEM) of Dziuk and Elliott [31], and states some previous results to be used in our subsequent analysis. In Section 5 we propose a spatially discrete numerical method for the approximation of 1 for which we prove stability (Lemma 8), continuous dependence on the data (Lemma 10), and error bounds (Theorem 10) under minimal regularity assumptions by introducing a new projection-type operator (Definition 5). To the authors’ knowledge, this is the first error bound for an approximation of a free boundary problem posed on an evolving surface. In Section 6 we propose a fully discrete numerical method and describe a novel implementation that avoids errors due to numerical quadrature of the method. This implementation has not previously been considered in the literature, even in the flat, stationary setting. Our numerical experiments indicate this method performs similarly to quadrature-based methods. This section is closed with numerical experiments illustrating solution behaviour that can only arise on an evolving domain, and tests that demonstrate the experimental order of convergence of the method.

2 The Stefan problem on an evolving surface↩︎

We shall assume throughout that \(\Omega(t)\) is a \(C^2\) evolving surface with known surface evolution, i.e. the surface evolution is independent of the solution to 1 . By \(\partial^\bullet\), we denote the material time derivative following the flow of \(\mathbf{V}\), \(\nabla_{\Omega}\) denotes the tangential gradient on \(\Omega(t)\), and \(\Delta_{\Omega}\) denotes the Laplace–Beltrami operator on \(\Omega(t)\). Note that due to the surface evolution the operators \(\nabla_{\Omega}\) and \(\Delta_{\Omega}\) vary in time, but we shall ignore this in our notation. We refer the reader to [31] for further details on these differential operators. We employ notion of the form \(L^p_{X}\), and \(H^1_{X}\), for evolving Bochner spaces, and evolving Sobolev–Bochner spaces respectively, where \(\{X(t)\}_{t \in [0,T]}\) is a family of time-dependent Banach spaces. This is defined in greater detail in Section 3.

2.1 Weak formulation↩︎

Given data \(e_0 \in L^\infty(\Omega(0))\) and \(f \in L^\infty_{L^\infty}\), one seeks a weak solution pair \((e, u) \in L^{\infty}_{L^\infty} \times L^2_{H^1}\) solving 2 , subject to the inclusion 3 , in the sense that \[\begin{align} -\int_0^T \int_{\Omega(t)} e \partial^\bullet\eta + \int_0^T \int_{\Omega(t)} \nabla_{\Omega}u \cdot \nabla_{\Omega}\eta = \int_0^T \int_{\Omega(t)} f \eta + \int_{\Omega(0)} e_0 \eta(0), \quad \text{such that} \quad e \in \beta(u), \end{align}\] for all \(\eta \in L^2_{H^1} \cap H^1_{L^2}\) with \(\eta(T) = 0\). We defer discussion on the function spaces used here to Section 3. Here \(\beta: \mathbb{R} \rightrightarrows \mathbb{R}\) is a set-valued map which satisfies the following assumptions.

By using the inverse function \(\mathcal{U}\) one finds that \(u = \mathcal{U}(e)\), and hence one may equivalently write 1 as a single equation \[\begin{align} \partial^\bullet e + e (\nabla_{\Omega}\cdot \mathbf{V}) - \Delta_{\Omega}\mathcal{U}(e) = f, \text{ on } \Omega(t). \label{eqn:32stefan32problem2} \end{align}\tag{5}\] In this setting one readily observes that since \(\mathcal{U}'(\cdot) \geq 0\) this equation is a degenerate, quasilinear, parabolic equation for \(e\). Using 5 one can now define a weak solution to 5 to be a function \(e \in L^2_{H^{-1}} \cap L^\infty_{L^2}\) such that \(\mathcal{U}(e) \in L^2_{H^1}\) and \[\begin{align} \left \langle \partial^\bullet e, \phi \right \rangle_{H^{-1}(\Omega(t)) \times H^{1}(\Omega(t))} + \int_{\Omega(t)} e \phi (\nabla_{\Omega}\cdot \mathbf{V}) + \int_{\Omega(t)} \nabla_{\Omega}\mathcal{U}(e) \cdot \nabla_{\Omega}\phi = \int_{\Omega(t)} f \phi, \label{eqn:32weak32stefan} \end{align}\tag{6}\] for all \(\phi \in H^1(\Omega(t))\) and almost all \(t \in [0,T]\) such that \(e(0) = e_0\), for some \(e_0 \in L^2(\Omega(0))\). Here angled brackets denote the duality pairing between \(H^{-1}(\Omega(t))\) and \(H^1(\Omega(t))\). This formulation will be the basis of our analysis. We note that in the case of weaker data, namely \(e_0 \in L^1(\Omega(0))\) and \(f \in L^1_{L^1}\), one may define a weaker notion of solution [1] for which our analysis is not applicable.

2.2 Strong formulation of a two-phase Stefan problem↩︎

As in [1] when \(\beta\) is given by ?? then a sufficiently smooth solution, \(u \in H^1_{L^2} \cap L^2_{H^2}\), can be shown to solve the classical strong formulation of the Stefan problem \[\tag{7} \begin{align} \partial^\bullet u + (u+1)(\nabla_{\Omega}\cdot \mathbf{V}) - \Delta_{\Omega}u = f, &\qquad \text{on } \bigcup_{t \in [0,T]} \Omega_+(t) \times \{t\}, \tag{8}\\ \partial^\bullet u + u(\nabla_{\Omega}\cdot \mathbf{V}) - \Delta_{\Omega}u = f, &\qquad \text{on } \bigcup_{t \in [0,T]} \Omega_-(t) \times \{t\}, \tag{9}\\ -(\nabla_{\Omega}u_+ - \nabla_{\Omega}u_-) \cdot \boldsymbol{\mu} = V_\mu, &\qquad \text{on } \bigcup_{t \in [0,T]} \Gamma(t) \times \{t\}, \tag{10}\\ u = 0, &\qquad \text{on } \bigcup_{t \in [0,T]} \Gamma(t) \times \{t\}, \tag{11} \end{align}\] pointwise almost everywhere¹, where \(u_\pm = u|_{\overline{\Omega_{\pm}}}\). Here \(\boldsymbol{\mu}\) denotes the unit conormal vector to \(\Gamma(t)\) (i.e. the vector which is tangential to \(\Omega(t)\), normal to \(\Gamma(t)\), and pointing into \(\Omega_+(t)\)), and \(V_\mu\) denotes the conormal velocity of \(\Gamma(t)\). We illustrate the relations between \(\Omega_{\pm}(t)\) and \(\Gamma(t)\) appearing in this strong formulation of the Stefan problem in Figure 2. This holds under the assumption that there are no so-called “mushy regions”, i.e. regions where the interior of \(\Gamma(t)\) is non-empty. In the presence of heat sources it is known that this assumption may not hold, cf. [16], [36], and indeed we shall observe numerically that mushy regions may form on an evolving surface even in the absence of a heat source (cf. Figure 4). In this setting it is clear that the enthalpy solutions we are approximating are generalised solutions of the strong formulation of Stefan problem, which does not make sense in the presence of mushy regions. Moreover, when a solution to the strong formulation 7 exists it is well known (see for instance [16], [19]) that this is also an enthalpy solution.

3 Evolving surfaces and function spaces↩︎

We assume that the evolving surface, \(\Omega(t)\), is \(C^2\) for all \(t \in [0,T]\), with a velocity field \(\mathbf{V} \in C^1([0,T];\mathbf{C}^2(\mathbb{R}^3;\mathbb{R}^3))\). It particular one finds that \[\sup_{t \in [0,T]} \|\mathbf{V}\|_{\mathbf{C}^2(\Omega(t);\mathbb{R}^3)} \leq C < \infty.\] We emphasise that this is the material velocity of the surface, \(\Omega(t)\), and is known a priori and hence independent of the free boundary.

One defines the parametrisation following the flow of \(\mathbf{V}\) as the unique solution, \(\Phi\colon [0,T] \times \Omega(0) \rightarrow \mathbb{R}^3\), of \[\frac{\mathrm{d}}{\mathrm{d}t}\Phi(t;\mathbf{x}) = \mathbf{V}(t; \Phi(t;\mathbf{x})) \quad \forall (t,\mathbf{x}) \in [0,T]\times\Omega(0), \quad \Phi(0;\mathbf{x}) = \mathbf{x} \quad \forall \mathbf{x} \in \Omega(0).\] By construction this is such that \(\Phi(t;\Omega(0)) = \Omega(t)\), and hence one has a parametrisation of the evolving surface. Using this parametrisation we may pushforward functions on \(\Omega(0)\) by \[\Phi_t \psi \mathrel{\vcenter{:}}= \psi \circ \Phi(t)^{-1} \quad \forall \psi: \Omega(0) \rightarrow \mathbb{R},\] and pullback functions on \(\Omega(t)\) by \[\Phi_{-t} \chi \mathrel{\vcenter{:}}= \chi \circ \Phi(t) \quad \forall \chi: \Omega(t) \rightarrow \mathbb{R}.\] Using these pushforwards/pullbacks one can define evolving Bochner spaces, which we shall denote as \(L^p_X\), as in [34], [35] as follows. Given a family of Banach spaces, \(\{X(t)\}_{t \in [0,T]}\), consisting of functions \[\chi: \bigcup_{t \in [0,T]} \Omega(t) \times \{ t \} \rightarrow \mathbb{R},\] one can define \(L^p_X\), for \(p\in [1,\infty]\), as the set of functions such that \(\Phi_{-t} \chi \in L^p([0,T];X(0))\), and denote the corresponding norm as \(\|\cdot\|_{L^p_X}\). We refer the reader to [34], [35] for further details on evolving Bocher spaces. We say that the function space-parametrisation pairs, \(\{(X(t), \Phi_t)\}_{t \in [0,T]}\), are compatible if we have also that:

1. there exists a constant \(C_X\) independent of \(t\) such that \[\begin{align} \|\Phi_t \psi\|_{X(t)} \leq C_X \|\psi\|_{X(0)} &\quad \forall \psi \in X(0),\\ \|\Phi_{-t} \chi\|_{X(0)} \leq C_X \|\chi\|_{X(t)} &\quad \forall \chi \in X(t), \end{align}\]

2. for all \(\chi \in X(0)\) the map \(t \mapsto \|\Phi_t \chi\|_{X(t)}\) is measurable.

We refer the reader to [34], [35] for further details on evolving function spaces. For compatible pairs, \(\{(X(t), \Phi_t)\}_{t \in [0,T]}\), \(L^p_X\) is a Banach space when equipped with the \(\|\cdot\|_{L^p_X}\) norm defined by \[\|\chi\|_{L^p_X} \mathrel{\vcenter{:}}= \begin{cases} \left( \int_0^T \|\chi(t)\|_{X(t)}^p \, \mathrm{d}t\right)^{\frac{1}{p}}, & p \in [1,\infty),\\ \underset{t \in [0,T]}{\operatorname{ess\,sup}}\, \|\chi(t)\|_{X(t)}, & p = \infty. \end{cases}\] Moreover, if \(p = 2\) and \(\{X(t)\}_{t \in [0,T]}\) is a family of Hilbert spaces then \(L^2_X\) is a Hilbert space, where one obtains the corresponding inner product by polarisation. Under the assumption of compatibility the function spaces \(L^p_X\) inherit many of the nice properties of the usual Bochner spaces — we refer the reader to [34] for details. We will be interested in the case where \(X(t)\) is a Sobolev space defined over \(\Omega(t)\), which we shall denote as \(W^{k,q}(\Omega(t))\) for \(k \in \mathbb{N}\) and \(q \in [1,\infty]\). It is well known [34], [35], [37] that our assumption on \(\mathbf{V}\) means that the pairs \(\{(W^{k,q}(\Omega(t)),\Phi_t)\}_{t\in[0,T]}\) are compatible in the above sense for \(k \in \{0,1\}\) and \(q \in [1,\infty]\). In the case \(q = 2\) we shall write \(H^k(\Omega(t)) \mathrel{\vcenter{:}}= W^{k,2}(\Omega(t))\). We refer the reader to [38] for further details on Sobolev spaces defined on Riemannian manifolds, and [37] for further details on Sobolev spaces on evolving hypersurfaces.

The natural notion of a time derivative on the evolving surface must take into account both the variation of the function in time, but also effects due to evolution of the surface. As such one works with the material derivative, defined as \[\partial^\bullet\chi \mathrel{\vcenter{:}}= \Phi_t \left( \frac{\mathrm{d}}{\mathrm{d}t} \Phi_{-t} \chi \right),\] for all sufficiently smooth functions \(\chi\). This notion can be generalised to a weak material derivative analogously to the usual time derivative, cf. [34]. We are particularly interested in the Gelfand triple setting, \(H^1(\Omega(t)) \subset L^2(\Omega(t)) \equiv (L^2(\Omega(t)))^* \subset H^{-1}(\Omega(t))\), and a weak time derivative taking values in \(H^{-1}(\Omega(t))\), the dual space to \(H^1(\Omega(t))\). We shall denote by \(H^1_{H^{-1}}\) the evolving Sobolev–Bochner space consisting of functions \[H^1_{H^{-1}} \mathrel{\vcenter{:}}= \{ \chi \in L^2_{L^2} \mid \partial^\bullet\chi \in L^2_{H^{-1}} \},\] and when the material derivative has further regularity \(\partial^\bullet\chi \in L^2_{L^2}\) we shall write \(\chi \in H^1_{L^2}\).

We end this section by recalling the transport theorem on an evolving surface, which we will use in our later analysis. For this we firstly introduce the following notation. \[\begin{align} m(t; \phi, \psi) &\mathrel{\vcenter{:}}= \int_{\Omega(t)} \phi \psi,\\ m_*(t; \mathcal{L}, \psi) &\mathrel{\vcenter{:}}= \left \langle \mathcal{L}, \psi \right \rangle_{H^{-1}(\Omega(t)) \times H^1(\Omega(t))}\\ g(t; \phi, \psi) &\mathrel{\vcenter{:}}= \int_{\Omega(t)} \phi \psi (\nabla_{\Omega}\cdot \mathbf{V}),\\ a(t; \phi, \psi) &\mathrel{\vcenter{:}}= \int_{\Omega(t)} \nabla_{\Omega}\phi \cdot \nabla_{\Omega}\psi,\\ b(t; \phi, \psi) &\mathrel{\vcenter{:}}= \int_{\Omega(t)} \left((\nabla_{\Omega}\cdot \mathbf{V}) \mathbb{I} - \nabla_{\Omega}\mathbf{V} - (\nabla_{\Omega}\mathbf{V})^T\right) \nabla_{\Omega}\phi \cdot \nabla_{\Omega}\psi, \end{align}\] for all sufficiently smooth functions \(\phi, \psi\), linear functionals \(\mathcal{L} \in H^{-1}(\Omega(t))\), and where \(\mathbb{I}\) denotes the identity matrix. We will omit the argument \(t\) throughout, as it will be clear from context.

4 The evolving surface finite element method↩︎

4.1 ESFEM and geometric perturbation estimates↩︎

Let us now briefly recap some of the details of the evolving surface finite element method. We refer the reader to [31], [33] for further details. Given a sufficiently smooth surface, \(\Omega(0)\), and a set of vertices \(\{\mathbf{x}_i(0)\}_{i \in \{1,\ldots, N_h\}} \subset \Omega(0)\), we may construct a triangulated domain by appropriately connecting these vertices by edges. Here \(N_h\) denotes the number of degrees of freedom. We denote the corresponding triangulation as \(\mathcal{T}_h(0)\), which in turn defines a triangulated surface, \(\Omega_h(0)\), via \[\Omega_h(0) \mathrel{\vcenter{:}}= \bigcup_{K \in \mathcal{T}_h(0)} K.\] One then may evolve the vertices by using the parametrisation, \(\Phi\), to obtain vertices \(\{\mathbf{x}_i(t)\}_{i \in \{1, \ldots, N_h\}}\) and a corresponding triangulation \(\mathcal{T}_h(t)\). One then defines the triangulated surface \[\Omega_h(t) \mathrel{\vcenter{:}}= \bigcup_{K(t) \in \mathcal{T}_h(t)} K(t).\] Notice that this construction yields a discrete velocity field, \(\mathbf{V}_h\), where is it straightforward, cf. [33], to see that this is the Lagrange interpolant of \(\mathbf{V}\). We shall denote the mesh size of our family of triangulated surfaces as \[h \mathrel{\vcenter{:}}= \sup_{t \in [0,T]} \max_{K(t) \in \mathcal{T}_h(t)} \operatorname{diam}(K(t)).\] Throughout our analysis we will assume that our triangulations of the evolving surface are uniformly quasi-uniform in the sense of [33].

We now denote our (continuous Lagrange) finite element spaces as \[S_h(t) \mathrel{\vcenter{:}}= \{ \phi: \Omega_h(t) \rightarrow \mathbb{R} \mid \phi|_{K(t)} \text{ is affine linear } \forall K(t) \in \mathcal{T}_h(t) \}.\]

4.1.1 Lifts and geometric perturbation estimates↩︎

Since our triangulated surface is not the true surface our associated finite element families are non-conforming. As such we are committing a so-called variational crime (cf. [39]) which we will mitigate by the use of lifts. For a surface, \(\Omega(t)\), there exists a neighbourhood of \(\Omega(t)\), denoted \(\mathcal{N}(\Omega(t)) \subset \mathbb{R}^3\) such that each \(\mathbf{x} \in \mathcal{N}(\Omega(t))\) may be uniquely expressed in Fermi coordinates \[\mathbf{x} = \mathbf{p}(t;\mathbf{x}) + d(t,\mathbf{x}) \boldsymbol{\nu}(t; \mathbf{p}(t; \mathbf{x})).\] Here \(\mathbf{p}(t;\cdot)\) denotes the closest point projection onto \(\Omega(t)\), \(d(t;\cdot)\) denotes the signed distance function of \(\Omega(t)\), and \(\boldsymbol{\nu}(t;\cdot)\) denotes the outward unit normal vector on \(\Omega(t)\). We now use this to “lift” a function from \(\Omega_h(t)\) onto \(\Omega(t)\) implicitly via \[\eta_h^\ell(\mathbf{p}(t;\mathbf{x})) = \eta_h(\mathbf{x}) \quad \forall \mathbf{x} \in \Omega_h(t), \eta_h: \Omega_h(t) \rightarrow \mathbb{R}.\] One can similarly define an inverse lift from \(\Omega(t)\) onto \(\Omega_h(t)\) by \[\eta^{-\ell}(\mathbf{x}) = \eta(\mathbf{p}(t; \mathbf{x})) \quad \forall \mathbf{x} \in \Omega_h(t), \eta : \Omega(t) \rightarrow \mathbb{R}.\] These operations are stable in \(W^{k,q}(\Omega(t))\) for \(k \in \{0,1\}\) and \(q \in [1,\infty]\) in the sense that there exist constants \(C_{k,q} > 0\), independent of \(h\), such that \[\begin{align} \frac{1}{C_{k,q}}\|\eta^{-\ell}\|_{W^{k,q}(\Omega_h(t))} \leq \|\eta\|_{W^{k,q}(\Omega(t))} \leq {C_{k,q}}\|\eta^{-\ell}\|_{W^{k,q}(\Omega_h(t))} \quad \forall \eta \in W^{k,q}(\Omega(t)). \label{eqn:32lift32stability} \end{align}\tag{12}\] We refer the reader to [33] for further details.

The evolution of the nodes of \(\Omega_h(t)\) by \(\mathbf{V}_h\) induces a parametrisation, \(\Phi^h(t): \Omega_h(0) \rightarrow \Omega^h(t)\), analogously to the definition of \(\Phi(t)\). This in turn defines a (strong) discrete material derivative by \[\partial^\bullet_h \chi_h \mathrel{\vcenter{:}}= \Phi_t^h \left( \frac{\mathrm{d}}{\mathrm{d}t} \Phi_{-t}^{h} \chi_h \right),\] for sufficiently smooth functions, \(\chi_h\), defined on \(\Omega_h(t)\). and one may define a weak discrete material derivative accordingly. Using this, one may state a discrete transport theorem, analogous to Lemma 1. For this we introduce the following notation for bilinear forms. \[\begin{align} m_h(t; \phi_h, \psi_h) &\mathrel{\vcenter{:}}= \int_{\Omega_h(t)} \phi \psi,\\ g_h(t; \phi_h, \psi_h) &\mathrel{\vcenter{:}}= \int_{\Omega_h(t)} \phi_h \psi_h (\nabla_{\Omega_h}\cdot \mathbf{V}_h),\\ a_h(t; \phi_h, \psi_h) &\mathrel{\vcenter{:}}= \int_{\Omega_h(t)} \nabla_{\Omega_h}\phi_h \cdot \nabla_{\Omega_h}\psi_h,\\ b_h(t; \phi_h, \psi_h) &\mathrel{\vcenter{:}}= \int_{\Omega_h(t)} \left((\nabla_{\Omega_h}\cdot \mathbf{V}_h) \mathbb{I} - \nabla_{\Omega_h}\mathbf{V}_h - (\nabla_{\Omega_h}\mathbf{V}_h)^T\right) \nabla_{\Omega_h}\phi_h \cdot \nabla_{\Omega_h}\psi_h. \end{align}\]

By combining the parametrisation \(\Phi^h(t)\), and the lifts one may define a lifted material derivative as \[\partial^\bullet_\ell \chi = (\partial^\bullet_h \chi^{-\ell})^\ell,\] for all sufficiently smooth functions, \(\chi\), defined on \(\Omega(t)\). We refer the reader to [33] for further details, and the proof of the following result relating \(\partial^\bullet\) and \(\partial^\bullet_\ell\).

One may equivalently formulate Lemma 1 using the lifted material derivative as follows.

Finally, to quantify the error induced by lifting functions we have the following geometric perturbation results.

4.2 Inverse Laplacians↩︎

In our later analysis we will make extensive use of inverse Laplacian operators. We now recall the definition of these operators.

We also have the following result concerning time-differentiability.

We will also require notions of a finite element inverse Laplacian defined on \(\Omega_h(t)\) as in the following definitions.

Next we recall that for \(\Sigma(t)\) a \(\mathcal{H}^2\)-measurable set we shall denote its \(\mathcal{H}^2\)-measure by \(|\Sigma(t)|\). For such a region, and a function \(z \in L^1(\Sigma(t))\), we define the mean value of \(z\) over \(\Sigma(t)\) as \[\mathchoice {{\setbox 0=\displaystyle{\textstyle-}{\int}\vcenter{\textstyle- }\kern-.6\wd 0}} {{\setbox 0=\textstyle{\scriptstyle-}{\int}\vcenter{\scriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptscriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} \!\int_{\Sigma(t)} z \mathrel{\vcenter{:}}= \frac{1}{|\Sigma(t)|}\int_{\Sigma(t)} z.\]

We can relate \(\mathcal{G}\) and \(\mathcal{G}_h\) through the following error estimate, in analogy to the usual inverse Laplacians on Euclidean domains.

One can similarly relate \(\mathcal{G}\) and \(\mathcal{G}_h\) to \(\mathcal{G}_{S_h}\), but this will not be useful in our subsequent analysis.

5 A semi-discrete evolving surface finite element method↩︎

We now introduce our spatially-discrete finite element method to be analysed. Given initial data, \(e_h^0 \in S_h(0)\), one finds \(e_h \in \mathcal{S}_h\), where \[\mathcal{S}_h \mathrel{\vcenter{:}}= \{ \chi_h \mid \Phi_{-t}^h \chi_h \in C^1([0,T];S_h(0)) \},\] such that \[\begin{align} m_h(\partial^\bullet_h e_h, \phi_h) + g_h(e_h, \phi_h) + a_h(\mathcal{U}(e_h), \phi_h) = m_h(f_h, \phi_h) \quad \forall \phi_h \in S_h(t), \label{eqn:32semidiscrete32stefan} \end{align}\tag{13}\] for almost all \(t \in [0,T]\), and such that \(e_h(0) = e_h^0\).

5.1 Stability↩︎

Next we will show a result concerning continuous dependence on the data. For this we require the following Ritz projection and error bound.

5.1.1 A method with quadrature↩︎

We now end this section with a result comparing our numerical method to a numerical method using a quadrature rule on the nonlinear term. Notice that we do note consider the effects of mass-lumping since our results will not use discrete maximum principle arguments. This is in contrast to previous numerical analysis of the Stefan problem [23]. For this we firstly recall the following result.

5.2 Error bounds↩︎

We now prove our main result, Theorem 10, for which we shall make the following assumption on the data.

Examples of data satisfying this are \(e_h^0 = \Lambda_h(0) e_0\) and \(f_h(t) = \Lambda_h(t) f(t)\), where \(\Lambda_h(t) : L^2(\Omega(t)) \rightarrow S_h(t)\) denotes the \(L^2\) projection onto \(S_h(t)\), or suitably defined² Ritz projections.

5.2.1 Preliminaries↩︎

We now introduce some preliminary results to be used in our subsequent error analysis. Firstly we define an \(L^2\) projection-type operator from \(\Omega_h(t)\) onto \(\Omega(t)\), which to our knowledge has not previously been used in the analysis of surface finite element methods.

Clearly such a function exists by the Riesz representation theorem. We note that this operator is not truly a projection since \(L^2(\Omega(t))\) is not a subset of \(L^2(\Omega_h(t))\) — however, it is “almost a projection” as we shall see in Lemma 13. We now state some of the basic properties of this operator.

The following technical result we prove concerns the differentiability in time of \(\mathscr{P}_hz_h\) for sufficiently smooth (in time) finite element functions \(z_h\).

5.2.2 Error analysis↩︎

To now derive suitable error equations let us introduce shorthand notation \[\mathscr{E}_h\mathrel{\vcenter{:}}= e - \mathscr{P}_he_h.\] We note that we may test 5 with \(\phi = \mathcal{G}\mathscr{E}_h\), which we observe is well-defined since \[\int_{\Omega(t)} e(t) = \int_{\Omega(0)} e_0 = \int_{\Omega_h(0)} e_h^0 = \int_{\Omega_h(t)} e_h(t) = \int_{\Omega(t)} \mathscr{P}_he_h(t).\] Likewise we may test we test 13 with \(\phi_h = \mathcal{G}_h \mathscr{E}_h\) which is well-defined for the same reason. Doing this one obtains \[\begin{gather} m_*(\partial^\bullet e, \mathcal{G}\mathscr{E}_h) + g(e, \mathcal{G}\mathscr{E}_h) + m(\mathcal{U}(e), \mathscr{E}_h) = m(f, \mathcal{G}\mathscr{E}_h), \tag{22}\\ m_h(\partial^\bullet_h e_h, \mathcal{G}_h\mathscr{E}_h) + g_h(e_h, \mathcal{G}_h\mathscr{E}_h) + m_h\left(R_h\mathcal{U}(e_h), \mathscr{E}_h^{-\ell} - \mathchoice {{\setbox 0=\displaystyle{\textstyle-}{\int}\vcenter{\textstyle- }\kern-.6\wd 0}} {{\setbox 0=\textstyle{\scriptstyle-}{\int}\vcenter{\scriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptscriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} \!\int_{\Omega_h(t)} \mathscr{E}_h^{-\ell}\right) = m_h(f_h, \mathcal{G}_h\mathscr{E}_h), \tag{23} \end{gather}\] where we have used the Definition 1, Definition 3 and Definition 4. We now rewrite 23 in terms of \(\mathscr{P}_he_h\) as \[\begin{align} m_*(\partial^\bullet\mathscr{P}_he_h, \mathcal{G}_h^\ell \mathscr{E}_h) + g(\mathscr{P}_he_h, \mathcal{G}_h^\ell\mathscr{E}_h) + m(\mathcal{U}(\mathscr{P}_he_h),\mathscr{E}_h) = m_h(f_h, \mathcal{G}_h \mathscr{E}_h) + \sum_{n=1}^5 I_n \label{eqn:32error32pf3} \end{align}\tag{24}\] where we have defined consistency errors \[\begin{gather} I_1 \mathrel{\vcenter{:}}= m_*(\partial^\bullet\mathscr{P}_he_h, \mathcal{G}_h^\ell \mathscr{E}_h) - m_h(\partial^\bullet_h e_h, \mathcal{G}_h\mathscr{E}_h),\\ I_2 \mathrel{\vcenter{:}}= g(\mathscr{P}_he_h, \mathcal{G}_h^\ell \mathscr{E}_h) - g_h(e_h, \mathcal{G}_h\mathscr{E}_h),\\ I_3 \mathrel{\vcenter{:}}= m(\mathcal{U}(\mathscr{P}_he_h), \mathscr{E}_h) - m_h(\mathcal{U}(e_h), \mathscr{E}_h^{-\ell}),\\ I_4 \mathrel{\vcenter{:}}= m_h(\mathcal{U}(e_h) - R_h \mathcal{U}(e_h), \mathscr{E}_h^{-\ell}),\\ I_5 \mathrel{\vcenter{:}}= \left( \int_{\Omega_h(t)} \mathcal{U}(e_h) \right)\left( \mathchoice {{\setbox 0=\displaystyle{\textstyle-}{\int}\vcenter{\textstyle- }\kern-.6\wd 0}} {{\setbox 0=\textstyle{\scriptstyle-}{\int}\vcenter{\scriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptscriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} \!\int_{\Omega_h(t)} \mathscr{E}_h^{-\ell} \right), \end{gather}\] so that when we subtract 24 from 22 we obtain \[\begin{align} m_*(\partial^\bullet\mathscr{E}_h, \mathcal{G}\mathscr{E}_h) + g(\mathscr{E}_h, \mathcal{G}\mathscr{E}_h) + m(\mathcal{U}(e) - \mathcal{U}(\mathscr{P}_he_h), \mathscr{E}_h) = \sum_{n = 1}^8 I_n, \label{eqn:32error32pf4} \end{align}\tag{25}\] where we have introduced further consistency errors defined by \[\begin{gather} I_6 \mathrel{\vcenter{:}}= m_*(\partial^\bullet\mathscr{P}_he_h, \mathcal{G} \mathscr{E}_h- \mathcal{G}_h^\ell \mathscr{E}_h),\\ I_7 \mathrel{\vcenter{:}}= g(\mathscr{P}_he_h, \mathcal{G} \mathscr{E}_h- \mathcal{G}_h^\ell \mathscr{E}_h),\\ I_8 \mathrel{\vcenter{:}}= m(f , \mathcal{G}\mathscr{E}_h) - m_h(f_h, \mathcal{G}_h \mathscr{E}_h). \end{gather}\] The main component of our error analysis is the following lemma concerning these consistency errors.

We now state, and prove, our main result.

6 Numerical examples↩︎

6.1 Numerical methods↩︎

In this section we consider the discretisation in time of 13 by a backward Euler time discretisation, similar to that of [49]. In the following we consider a uniform timestep size \(\tau = \frac{T}{N_T}\) for some \(N_T \in \mathbb{N}\). We also introduce some shorthand notation for the current time, \(t_n \mathrel{\vcenter{:}}= n \tau\), and the space of finite element functions at time \(t_n\), \(S_h^n \mathrel{\vcenter{:}}= S_h(t_n)\). The fully discrete problem is now as follows: Given data \(e_h^{n-1} \in S_h^{n-1}\) and \(f_h^n \in S_h^n\), find \(e_h^n \in S_h^n\) such that \[\begin{align} \frac{1}{\tau} \left(m_h(t_n; e_h^n, \phi_h^n) - m_h(t_{n-1}; e_h^{n-1}, \underline{\phi_h^n})\right) + a_h(t_n;\mathcal{U}(e_h^n), \phi_h^n) = m_h(t_n; f_h^n, \phi_h^n), \label{eqn:32fullydiscrete32stefan} \end{align}\tag{31}\] for all \(\phi_h^n \in S_h^n\). Here \(\underline{\phi_h^n} \in S_h^{n-1}\) denotes the vector with the same nodal values as \(\phi_h^n \in S_h^n\) but defined over the previous surface. We do not analyse this fully discrete numerical method, but we expect that the analysis will follow by combining techniques introduced in this present work with those developed for fully discrete ESFEM, cf. [44], [49]–[52]. We now introduce two numerical methods to approximate 31 .

6.1.1 A method using quadrature↩︎

For our first numerical method, we introduce a quadrature rule for the nonlinear term, as we did in Lemma 12. In matrix-vector form, the fully discrete scheme 31 with a quadrature rule may be written as \[\begin{align} M^n \boldsymbol{\mathsf{e}}^n + \tau A^n \mathcal{U}(\boldsymbol{\mathsf{e}}^n) = M^{n-1} \boldsymbol{\mathsf{e}}^{n-1} + \tau M^n \boldsymbol{\mathsf{f}}^n. \label{eqn:32matrix32vector32form32quadrature} \end{align}\tag{32}\] where \(\boldsymbol{\mathsf{e}}^n\) and \(\boldsymbol{\mathsf{f}}^n\) denote the vector of nodal values for \(e_h^n \in S_h^n\) and \(f_h^n \in S_h^n\) respectively, and the mass and stiffness matrices are given by components \[M^n_{ij} = m_h(t_n; \phi_i^n, \phi_j^n), \quad A^n_{ij} = a_h(t_n; \phi_i^n, \phi_j^n),\] for \(\phi_i^n\) the ’\(i\)’th basis function of \(S_h^n\). We implement this method in DUNE [53], solving the nonlinear problem by using a non-smooth Newton method [54] where the corresponding linear problems are solved with an exact solver and our tolerance for the Newton solver is chosen as \(10^{-7}\).

6.1.2 A method using an exact discretisation↩︎

We now consider a new approach to the discretisation of this problem, based on writing 5 as \[\begin{align} \partial^\bullet e + e (\nabla_{\Omega}\cdot \mathbf{V}) - \nabla_{\Omega}\cdot( \mathcal{U}'(e) \nabla_{\Omega}e) = f. \end{align}\]

In matrix-vector form this yields a system of the form \[\begin{align} M^n \boldsymbol{\mathsf{e}}^n + \tau \mathcal{A}(t_n;\boldsymbol{\mathsf{e}}^n)\boldsymbol{\mathsf{e}}^n = M^{n-1} \boldsymbol{\mathsf{e}}^{n-1} + \tau M^n \boldsymbol{\mathsf{f}}^n, \label{eqn:32matrix32vector32form32exact} \end{align}\tag{33}\] where we have defined the solution-dependent matrix \[\mathcal{A}(t_n;\boldsymbol{\mathsf{e}}^n)_{ij} = \int_{\Omega_h(t_n)} \mathcal{U}'(e_h^n) \nabla_{\Omega_h}\phi_i \cdot \nabla_{\Omega_h}\phi_j.\] To construct this matrix for general functions \(e_h^n\), or \(\mathcal{U}\) may be challenging. In the case that \(e_h^n\) is piecewise linear, one may divide each of the simplices into smaller regions in which the integration of \(\mathcal{U}'(e_h^n)\) follows via simple quadrature. For example, when considering \(\mathcal{U}\) as given in ?? , a piecewise polynomial function, one divides each of the simplices into the polygons/polyhedra in which \(e_h^n<0\) and \(e_h^n\geq1\), in these regions, exact integration is the result of a classical quadrature rule. This strategy extends to a wide class of functions \(\mathcal{U}\) which are piecewise smooth and exact quadrature rules available on each smooth region. We note that in applications it is typically the case that \(\mathcal{U}\) is piecewise linear, cf. [19].

With this matrix assembled, the problem is solved by a fixed point iteration where we solve sequences of linear problems \[M^n \boldsymbol{\mathsf{e}}^{n,k} + \tau \mathcal{A}(t_n;\boldsymbol{\mathsf{e}}^{n,k-1})\boldsymbol{\mathsf{e}}^{n,k} = M^{n-1} \boldsymbol{\mathsf{e}}^{n-1} + \tau M^n \boldsymbol{\mathsf{f}}^n,\] with an exact solver, provided \(\tau\) is sufficiently small. We choose our initial guess as \(\boldsymbol{\mathsf{e}}^{n,0} = \boldsymbol{\mathsf{e}}^{n-1}\), and choose our stopping criteria to be \[|M^n \boldsymbol{\mathsf{e}}^{n,k} + \tau \mathcal{A}(t_n;\boldsymbol{\mathsf{e}}^{n,k})\boldsymbol{\mathsf{e}}^{n,k} - M^{n-1} \boldsymbol{\mathsf{e}}^{n-1} - \tau M^n \boldsymbol{\mathsf{f}}^n| < \mathsf{tol},\] where we set the tolerance for our fixed point solver to be \(\mathsf{tol} = 10^{-7}\).

This method is somewhat similar to a recent method designed in [55] where the authors solve a linear system determined by a “flag” indicating which phase the mesh point is in. This approach also relies on the piecewise linear nature of \(\mathcal{U}\).

Due to the fixed point nature of this algorithm, we found that this approach was slower than the approach using quadrature. Our implementation of this fixed point method found that, for sufficiently small \(\tau\), this method converged within a few iterations. We leave the topic of more efficient solvers based on this approach for future work.

6.2 Examples on an evolving surface↩︎

In this subsection we demonstrate some phenomena exhibited by the evolving surface Stefan problem which arise only on evolving domains.

6.2.1 Nucleation in the absence of external heat sources↩︎

We now demonstrate a phenomenon which cannot happen on a stationary domain. In the following we shall assume \(f \equiv 0\). Consider the graph \(\beta\) given by \[\begin{align} \beta(r) \mathrel{\vcenter{:}}= \begin{cases} \{r\}, & r < 1,\\ [1,2], & r = 1,\\ \{ r + 1\}, & r > 1, \end{cases} \label{eqn:32nucleation95graph} \end{align}\tag{34}\] and initial data \(e_0 \in L^\infty(\Omega(0))\) which takes values in \([1,1+\delta]\) for some \(\delta \in (0,1)\). On a stationary domain the maximum principle (see for instance [56]) implies that the solution to 5 must be valued in \([1,1+\delta]\), and in particular the temperature is \(u \equiv 1\) for all time. However, on an evolving domain the usual maximum principle does not apply and indeed the surface evolution can now force the enthalpy to take arbitrary values in \(\mathbb{R}^+\). To illustrate this if one has initial data as above, then testing 6 with \(\phi \equiv 1\) and using Lemma 1 one finds \[0 = m_*(\partial^\bullet e ,1) + g(e,1) = \frac{\mathrm{d}}{\mathrm{d}t} \int_{\Omega(t)} e(t).\] Hence one has that \[\int_{\Omega(t)} e(t) = \int_{\Omega(0)} e_0 \quad \text{for a.e. } t \in [0,T],\] which one may rewrite as \[\mathchoice {{\setbox 0=\displaystyle{\textstyle-}{\int}\vcenter{\textstyle- }\kern-.6\wd 0}} {{\setbox 0=\textstyle{\scriptstyle-}{\int}\vcenter{\scriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptscriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} \!\int_{\Omega(t)} e(t) = \frac{|\Omega(0)|}{|\Omega(t)|} \mathchoice {{\setbox 0=\displaystyle{\textstyle-}{\int}\vcenter{\textstyle- }\kern-.6\wd 0}} {{\setbox 0=\textstyle{\scriptstyle-}{\int}\vcenter{\scriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptscriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} \!\int_{\Omega(0)} e_0.\] In particular if the surface area, \(|\Omega(t)|\), decreases enough so that there exists \(t^* \in [0,T]\) such that \[|\Omega(t^*)| < \frac{1}{2}\int_{\Omega(0)} e_0,\] then one finds \[\mathchoice {{\setbox 0=\displaystyle{\textstyle-}{\int}\vcenter{\textstyle- }\kern-.6\wd 0}} {{\setbox 0=\textstyle{\scriptstyle-}{\int}\vcenter{\scriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptscriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} \!\int_{\Omega(t^*)} e(t^*) = \frac{|\Omega(0)|}{|\Omega(t^*)|} \mathchoice {{\setbox 0=\displaystyle{\textstyle-}{\int}\vcenter{\textstyle- }\kern-.6\wd 0}} {{\setbox 0=\textstyle{\scriptstyle-}{\int}\vcenter{\scriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptscriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} \!\int_{\Omega(0)} e_0 > 2.\] Hence there must exist a region of positive measure on \(\Omega(t^*)\) such that \(e(t^*) > 2\), meaning that the surface evolution has caused the \(\{u > 1\}\) phase to nucleate. Conversely, if we consider an expanding surface, where there exists some \(t^* \in [0,T]\) such that \[|\Omega(t^*)| > \int_{\Omega(0)} e_0,\] then one finds that \[\mathchoice {{\setbox 0=\displaystyle{\textstyle-}{\int}\vcenter{\textstyle- }\kern-.6\wd 0}} {{\setbox 0=\textstyle{\scriptstyle-}{\int}\vcenter{\scriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptscriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} \!\int_{\Omega(t^*)} e(t^*) = \frac{|\Omega(0)|}{|\Omega(t^*)|} \mathchoice {{\setbox 0=\displaystyle{\textstyle-}{\int}\vcenter{\textstyle- }\kern-.6\wd 0}} {{\setbox 0=\textstyle{\scriptstyle-}{\int}\vcenter{\scriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} {{\setbox 0=\scriptscriptstyle{\scriptscriptstyle-}{\int}\vcenter{\scriptscriptstyle- }\kern-.6\wd 0}} \!\int_{\Omega(0)} e_0 < 1,\] and hence the \(\{u < 1\}\) phase has nucleated.

We demonstrate this phenomena in Figure 3 in the following numerical example, solving 32 , on an ellipsoid given by the \(0\)-level set of \[\phi(x,y,z;t) = \frac{2x^2}{1 + \exp\left(\frac{3t}{2}\right)} + \frac{y^2}{(2 - \tanh(2t))^2} + z^2 - 1,\] over a time interval \(t \in [0,1]\), with initial data \(e_0(x,y,z) \equiv 1.5\). Here we consider a mesh with \(h \approx 0.3878\), and a timestep size \(\tau = 2\cdot 10^{-4}\). The surface evolution here has no tangential component, and hence the nucleation of phases is due to geometric motion of the domain rather than advective effects.

6.2.2 Formation of mushy regions in the absence of external heat sources↩︎

In the stationary, Euclidean setting it is known [16] that one may develop a mushy region in the presence of an external heat source. Conversely, it is known [20], [57], that in the absence of external heat sources, and sufficiently smooth data, that mushy regions do not spontaneously develop. However, on an evolving surface it appears to be the case that mushy regions can spontaneously develop even when one does not have an external heat source. We demonstrate this in Figure 4, where we solve 5 with \(f \equiv 0\), and choose \(\beta\) given by ?? . Heuristically this can be explained by thinking of the surface evolution of \(\Omega(t)\) acting as an external heat source in some sense. This implies that the existence of strong solutions, for which there cannot be a mushy region, to the Stefan problem on an evolving surface is quite delicate, and to our knowledge there are indeed no results of this kind.

In Figure 4 we demonstrate the formation of such a mushy region in the absence of external heat sources. Here we consider an evolving torus, with major radius \(R(t) = 0.75 + 0.75t\) and minor radius \(r(t) \equiv 0.25\), and initial data \[e_0(x,y,z) = \begin{cases} \tanh(10(x-0.4)), & x \leq 0.4,\\ \tanh(10(x-0.4)) + 1, & x > 0.4. \end{cases}\] We consider a mesh with \(h \approx 0.0772\) and a timestep size \(\tau = 10^{-4}\).

6.3 Experimental order of convergence↩︎

In the following we choose our timestep size to be \(\tau = \mathcal{O}(h^2)\) so that the error we observe is dominated to the spatial discretisation, as considered in our analysis. To approximate the \(L^2_{L^2}\) norm of a function \(z_h\) (with nodal vector \(\boldsymbol{\mathsf{z}}^n\) at \(t = t_n\)) we shall use the approximation \[\|z_h\|_{L^2_{L^2}} \approx \left( \sum_{n=1}^{N_T} \tau \|z_h^n\|_{L^2(\Omega_h(t))}^2 \right)^{\frac{1}{2}} = \left( \sum_{n=1}^{N_T} \tau \boldsymbol{\mathsf{z}}^n \cdot M^n \boldsymbol{\mathsf{z}}^n \right)^{\frac{1}{2}}.\] Similarly we approximate the \(L^{\infty}_{H^{-1}}\) norm via \[\|z_h\|_{L^{\infty}_{H^{-1}}}^2 \approx \max_{n=1,\ldots,N_T} \left(\boldsymbol{\mathsf{z}}^n \cdot M^n (M^n + A^n)^{-1}M^n \boldsymbol{\mathsf{z}}^n \right).\] This right-hand expression can be viewed as a finite element approximation \(m(t_n;z,\widetilde{\mathcal{G}}z) \approx \|z\|_{H^1(t_n)}^2\) for a suitably defined inverse Laplacian operator \(\widetilde{\mathcal{G}}\). One can then justify the use of this approximation by appealing to well-known error bounds for surface finite elements [31] — we omit such justification.

Our experimental order of convergence (EOC) plots (Figure [fig:32EOC95expanding95sphere] and Figure [fig:32EOC95rotating95sphere]) indicate a higher order of convergence than that predicted by Theorem 10. This has previously been observed in [57], where the authors observe \(\mathcal{O}(h)\) convergence for the temperature in the \(L^2_{L^2}\) norm. For comparison, in our experiments we observe \(\mathcal{O}(h^2)\) convergence for the temperature in the \(L^2_{L^2}\) norm — we expect this is due to our exact solution being quite smooth, as it is difficult to construct rougher closed-form solutions to a surface Stefan problem.

6.3.1 Stationary interface on a expanding sphere↩︎

Here we consider our spatial domain to be an expanding sphere, given by the zero level set of \[\phi(x,y,z;t) = x^2 + y^2 + z^2 - e^{4t},\] and choose our initial data to be given by the Lagrange interpolant of \(e_0\), where \[e_0(x,y,z) = \begin{cases} P_1(y) + \frac{1}{2}P_3(y), & y \leq 0,\\ 1 + P_1(y) + \frac{1}{2}P_3(y), & y > 0, \end{cases}\] where \(P_n(\cdot)\) denotes the ’\(n\)’th Legendre polynomial. For this given choice of initial data, and domain evolution, an exact solution of the Stefan problem is known and constructed in §7.1.

6.3.2 Rotating sphere↩︎

Here we demonstrate experimental order of convergence results using the rotating sphere solution constructed in §7.2. We consider the evolution of the unit sphere under the map \(\Phi\) given by \[\Phi(\mathbf{x};t) = \begin{pmatrix} 1 & 0 & 0\\ 0 & \cos(t) & \sin(t)\\ 0 & -\sin(t) & \cos(t) \end{pmatrix} \mathbf{x}.\] Here we choose our initial data to be given by the Lagrange interpolant of \(e_0\), where \[e_0(x,y,z) = \begin{cases} P_1(y) + \frac{1}{2}P_3(y), & y \leq 0,\\ 1 + P_1(y) + \frac{1}{2}P_3(y), & y > 0, \end{cases}\] where \(P_n(\cdot)\) denotes the ’\(n\)’th Legendre polynomial.

Acknowledgments↩︎

The authors would like to thank Vanessa Styles and James Van Yperen for their comments on an early draft of this manuscript. TS is supported by the UK Engineering and Physical Sciences Research Council (Grant number: EP/Z535138/1). PJH, TS, and CV are all supported by a UK Engineering and Physical Sciences Research Council Mathematical Sciences Small Grant.

Conflict of interest↩︎

There are no conflicts of interest to declare.

Data availability↩︎

The code used for the numerical examples in this paper is available upon reasonable request to the authors.

7 Exact solutions for the Stefan problem on an evolving sphere↩︎

In this appendix we construct some exact solutions for the Stefan problem on an evolving sphere. We consider two kinds of evolution: uniform expansion/dilation in the normal direction, and tangential motion given by a uniform rotation. Note that this latter case corresponds to a Stefan problem with advection posed on a stationary surface. We also refer the reader to [3], [58] wherein exact solutions to surface free boundary problems (namely the Mullins–Sekerka problem) are constructed for use as a benchmark.

7.1 A sphere with varying radius↩︎

Here we manufacture a solution for the the strong formulation of the Stefan problem, 7 , on a shrinking sphere, with a fixed interface at \(z=0\), to be used in the experimental order of convergence calculations in Section 6. Notice that although the calculations in [1] do not consider an external heat source, given that our manufactured solution will be chosen such that there are no mushy regions we will be able to verify that it is also a solution of the enthalpy formulation. We shall consider a ball of radius \(\rho(t)\), centred at the origin, which we shall denote as \(S^2(\rho(t))\) To construct our solution we will firstly solve the heat equation 8 on the upper hemisphere \(S^2(\rho(t)) \cap \{z \geq 0\}\) subject to homogeneous Dirichlet boundary conditions on \(S^2(\rho(t)) \cap \{z=0\}\). For this we pullback 8 onto the unit sphere \(S^2(1)\), noting that as we assume evolution is exclusively in the normal direction one finds \[\nabla_{\Omega}\cdot \mathbf{V} = HV_N = \frac{2\rho'(t)}{\rho(t)},\] where \(V_N = \rho'(t)\) is the normal velocity of the sphere and \(H = \frac{2}{\rho(t)}\) is (twice) the mean curvature. By using the definition of the material derivative, and the pullback equation for the Laplace–Beltrami operator [59], one can readily observe that the pullback \(\widetilde{u}\) solves \[\frac{\partial \widetilde{u}}{\partial t} + 2(\widetilde{u} + 1)\frac{\rho'(t)}{\rho(t)} - \frac{1}{\rho(t)^2} \Delta_{S^2} \widetilde{u} = \widetilde{f},\] where \(\widetilde{f}\) is the pullback of \(f\) onto \(S^2(1)\), and \(\Delta_{S^2}\) is the Laplace–Beltrami operator on \(S^2(1)\).

We now choose \(\widetilde{u}\) as \[\widetilde{u}(x,y,z;t) = e^{-2t}P_{1}(z) + \frac{1}{2}e^{-12t}P_{3}(z),\] where \(P_n(\cdot)\) denotes the ’\(n\)’th Legendre polynomial. It is well-known (see for instance [60]) that the eigenfunctions of the Laplace–Beltrami operator on \(S^2(1)\) are the spherical harmonics, which includes the family \(\{P_n(z)\}_{n \in \mathbb{N} \cup \{0\}}\), where one finds that \[\Delta_{S^2} P_n(z) = -n(n+1) P_n(z).\] Moreover, by choosing only odd degree spherical harmonics we also satisfy the Dirichlet boundary condition, owing to the property that \(P_{2n+1}(0) = 0\) for all \(n \in \mathbb{N}\). As such, we now take \(\widetilde{f}\) on the upper hemisphere \(S^2(1) \cap \{ z \geq 0\}\) to be given by \[\begin{align} \widetilde{f}(x,y,z;t) &= (-2e^{-2t}P_{1}(z)-6 e^{-12t}P_{3}(z))\left(1-\frac{1}{\rho(t)^2}\right)+ (2+ 2e^{-2t}P_{1}(z) + e^{-12t}P_{3}(z))\frac{\rho'(t)}{\rho(t)}\\ &=: \widetilde{f}_+(z;t). \end{align}\]

Repeating essentially the same calculations for 9 we find that \[\begin{align} \widetilde{f}(x,y,z;t) &= (-2e^{-2t}P_{1}(z)-6 e^{-12t}P_{3}(z))\left(1-\frac{1}{\rho(t)^2}\right)+ (2e^{-2t}P_{1}(z) + e^{-12t}P_{3}(z))\frac{\rho'(t)}{\rho(t)}\\ & =: \widetilde{f}_-(z;t), \end{align}\] on the lower hemisphere \(S^2(1) \cap \{ z < 0\}\). All that remains is to push-forward our functions onto \(S^2(\rho(t))\), for which we find that \[u(x,y,z;t) = e^{-2t}P_{1}\left(\frac{z}{\rho(t)}\right) + \frac{1}{2}e^{-12t}P_{3}\left(\frac{z}{\rho(t)}\right),\] and \[\begin{align} f(x,y,z;t) = \begin{cases} \widetilde{f}_+\left( \frac{z}{\rho(t)}; t \right), & z \geq 0,\\ \widetilde{f}_-\left( \frac{z}{\rho(t)}; t \right), & z < 0. \end{cases} \label{eqn:32manufactured32solution32rhs} \end{align}\tag{35}\] We note that the values \(f\) takes on the interface \(S^2(\rho(t)) \cap \{z = 0 \}\) are seemingly irrelevant as the PDEs are defined away from this interface. Moreover, since the Laplace–Beltrami operator is invariant under isometries, we may choose the interface to be any great circle rather than specifically \(\{z = 0\}\). We shall do this in our examples, where we instead choose the interface to be \(\{y = 0\}\), so that the interface is not fitted to our mesh. In our numerical experiments in Section 6 we shall use this manufactured solution with \(\rho(t) = e^{2t}\).

7.2 A rotating sphere↩︎

As a second benchmark in Section 6 we shall consider \(\Omega(t)\) to be a rotating unit sphere. In this case the evolution of the sphere is given by \[\mathbf{V}(\mathbf{x};t) = A \mathbf{x}, \quad \forall \mathbf{x} \in S^2(1),\] where \(\mathbf{A} \in \mathbb{R}^{3 \times 3}\) is some given antisymmetric matrix. It is a straightforward calculation to verify that \(\nabla_{\Omega}\cdot \mathbf{V} = 0\), and hence the temperature, \(u\), solves the same PDE on both sides of the interface. We will solve the PDE by pulling back onto \(S^2(1)\), as we did above, where we note that \[\Phi_{-t}(\Delta_{\Omega}u) = \Delta_{S^2} (\Phi_{-t}u),\] since \(\Phi(t) : S^2(1) \rightarrow \Omega(t)\) is an isometry. Thus, the pullback of 8 is now \[\frac{\partial \widetilde{u}}{\partial t} - \Delta_{S^2} \widetilde{u} = 0,\] where \(\widetilde{u} = \Phi_{-t}u\) is the pullback of \(u\) onto \(S^2(1)\). Again by noting that the functions \(P_n(z)\) are eigenfunctions of the Laplace–Beltrami operator on \(S^2(1)\) we find that the above PDE is solved by \[\widetilde{u}(x,y,z;t) = e^{-2t} P_1(z) + \frac{1}{2} e^{-12t} P_3(z),\] which also satisfies the homogeneous Dirichlet condition on \(z = 0\). Moreover, this is a solution on both sides of the interface, i.e. \(S^2(1) \cap \{ \widetilde{u} < 0\}\) and \(S^2(1) \cap \{ \widetilde{u} > 0\}\). Thus a strong solution to the Stefan problem on \(\Omega(t)\) is given by \[u(x,y,z;t) = e^{-2t} P_1(\Phi_t z) + \frac{1}{2} e^{-12t} P_3(\Phi_t z),\] where the free boundary is the curve \(\Gamma(t) = S^2(1) \cap \{\Phi_t z = 0\}\). As above, we may we may choose the initial interface to be any great circle rather than \(\{z = 0\}\).

References↩︎

---

1. Here we understand almost everywhere to mean up to a set of zero \(\mathcal{H}^2\) measure for 8 , 9 , and up to a set of zero \(\mathcal{H}^1\) measure for 10 , 11 .↩︎

2. For example the Ritz projection of [33], rather than that in our Definition 4.↩︎