Moving boundary problems for quasi-steady conduction limited melting

TL;DR

This paper models the melting of dendritic crystals as a quasi-steady Stefan problem, deriving asymptotic results and analyzing stability, supported by numerical simulations and motivated by microgravity experiments.

Contribution

It extends previous asymptotic results for crystal melting, incorporates surface tension effects, and models non-monotonic aspect ratio behavior observed in experiments.

Findings

Asymptotic results for complete melting derived using Baiocchi transform.

Numerical level set simulations support theoretical predictions.

Model reproduces experimentally observed non-monotonic aspect ratio changes.

Abstract

The problem of melting a crystal dendrite is modelled as a quasi-steady Stefan problem. By employing the Baiocchi transform, asymptotic results are derived in the limit that the crystal melts completely, extending previous results that hold for a special class of initial and boundary conditions. These new results, together with predictions for whether the crystal pinches off and breaks into two, are supported by numerical calculations using the level set method. The effects of surface tension are subsequently considered, leading to a canonical problem for near-complete-melting which is studied in linear stability terms and then solved numerically. Our study is motivated in part by experiments undertaken as part of the Isothermal Dendritic Growth Experiment, in which dendritic crystals of pivalic acid were melted in a microgravity environment: these crystals were found to be prolate spheroidal in shape, with an aspect ratio initially increasing with time then rather abruptly decreasing to unity. By including a kinetic undercooling-type boundary condition in addition to surface tension, our model suggests the aspect ratio of a melting crystal can reproduce the same non-monotonic behaviour as that which was observed experimentally.

Equations

$$\beta =\frac{L}{c(u_{\infty }^{\ast }-u_{\mathrm{m}}^{\ast })},$$

$$u^{\ast }=u_{\mathrm{m}}^{\ast }(1-\gamma {\kappa }^{\ast })\text{mbox}on\partial {\Omega }^{\ast },$$

$$t=\frac{k}{{\ell }^2\beta }{t}^{\ast },\mathbf{x}=\frac{1}{\ell }{\mathbf{x}}^{\ast },u=\frac{u^{\ast }-u_{\mathrm{m}}^{\ast }}{u_{\infty }^{\ast }-u_{\mathrm{m}}^{\ast }},$$

$$\text{mbox}in{\mathbb{R}}^3\setminus \Omega (t):$$

$$\text{mbox}on\partial \Omega :$$

$$\text{mbox}asr\to \infty :$$

$$\text{mbox}in{\mathbb{R}}^3\setminus \Omega (t):{\nabla }^{2}u=0,$$

$$\text{mbox}in{\mathbb{R}}^3\setminus \Omega (0):$$

$$\text{mbox}in\Omega (0)\setminus \Omega (t):$$

$$\text{mbox}in{\mathbb{R}}^3\setminus \Omega (0):$$

$$\text{mbox}in\Omega (0)\setminus \Omega (t):$$

$$\text{mbox}on\partial \Omega :$$

$$\text{mbox}asr\to \infty :$$

$$x^2+y^2+\frac{z^2}{z_0(0{)}^2}=1,$$

$$\frac{x^2+y^2}{{\rho }_0(t{)}^2}+\frac{z^2}{z_0(t{)}^2}=1,$$

$$\frac{z_0(t)}{z_0(0)}={\rho }_0(t)=\sqrt{1-\frac{t}{t_{\mathrm{e}}}},$$

$$t_{\mathrm{e}}=\frac{z_0(0)}{4\sqrt{z_0(0{)}^2-1}}\ln \left(\frac{z_0(0)+\sqrt{z_0(0{)}^2-1}}{z_0(0)-\sqrt{z_0(0{)}^2-1}}\right).$$

$$\text{mbox}in{\mathbb{R}}^3\setminus \Omega (0):$$

$$\text{mbox}in\Omega (0):$$

$$\text{mbox}asr\to \infty :$$

$$W=-\frac{1}{4\pi }\int \int {\int }_{\Omega (0)}\frac{1}{|\mathbf{x}-{\mathbf{x}}^{\prime }|}\mathrm{d}{V}^{\prime },$$

$$w_{\mathrm{e}}\sim a(x^2+y^2)+\left(\frac{1}{2}-2a\right)z^2\text{mbox}asr\to 0,$$

$$\text{mbox}in{\mathbb{R}}^3\setminus {\Omega }_0(0):$$

$$\text{mbox}on\partial {\Omega }_0:$$

$$\Phi \sim a(X^2+Y^2)+\left(\frac{1}{2}-2a\right)Z^2-d+\frac{1}{3R},$$

$$w\sim a(x^2+y^2)+\left(\frac{1}{2}-2a\right)z^2-d{T}^2+\frac{T^3}{3r}\text{mbox}asr\to 0.$$

$$\frac{X^2+Y^2}{q_0^2-1}+\frac{Z^2}{q_0^2}=\frac{1}{q_0^{2/3}(q_0^2-1{)}^{2/3}},$$

$$a=\frac{1}{4}{q}_0^2-\frac{1}{8}{q}_0(q_0^2-1)\ln \left(\frac{q_0+1}{q_0-1}\right).$$

$$d=\frac{1}{4}{q}_0^{1/3}(q_0^2-1{)}^{1/3}\ln \left(\frac{q_0+1}{q_0-1}\right).$$

$$w\sim w_{\mathrm{e}}-(t-t_{\mathrm{e}})+\frac{T^2}{3r}\text{mbox}ast\to t_{\mathrm{e}}^{-}.$$

$$t=t_{\mathrm{e}}-d{T}^2+\mathcal{O}(T^5)\text{mbox}asT\to 0,$$

$$T\sim \frac{1}{\sqrt{d}}(t_{\mathrm{e}}-t{)}^{1/2}\text{mbox}ast\to t_{\mathrm{e}}^{-}.$$

$$\mathcal{A}(t_{\mathrm{e}})=\frac{q_0}{\sqrt{q_0^2-1}},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\int \int {\int }_{{\mathbb{R}}^3\setminus \Omega (t)}\Phi (\mathbf{x})\mathrm{d}V=0,$$

Full text

{AMS}

35R37, 80A22, 65M99

1 Introduction

While there is a variety of simple models to approximate the shape of a melting particle [33, 38], the traditional approach from a mathematical perspective is to employ a Stefan problem, which involves the linear heat equation subject to appropriate boundary conditions on the solid-melt interface. These moving boundary problems are well studied via rigorous analysis, asymptotic techniques, some exact solutions and numerical computation. Almost all of the analytical progress has been made for one-dimensional problems or those with radial symmetry [23, 39, 47, 48, 55], although there have been successful studies in which the symmetry is broken [37, 44, 46, 56]. We continue this direction in the present study, focusing on the melting of an axially symmetric dendritic crystal. We employ both analytical and numerical techniques to study the shape of the evolving crystal, focussing on the very final stages of melting.

A key aspect of a traditional Stefan problem is that the effects of convection are ignored. An excellent example of a relevant physical application involves certain experiments undertaken on the space shuttle Columbia, as part of the so-called Isothermal Dendritic Growth Experiment (IDGE) [21, 22, 43], in which convection is not an issue. The conduction-limited melting that was studied in those experiments provides a physical motivation for the kind of theoretical Stefan problems considered here. A brief summary of these experiments is as follows. A pure liquid melt, pivalic acid, is held at a temperature $u^{*}>u_{\mathrm{m}}^{*}$, where $u_{\mathrm{m}}^{*}\approx 35.9$ ∘C is the equilibrium melting temperature. The temperature is then reduced to slightly supercool the melt so that $u\lesssim u_{\mathrm{m}}$ throughout. The growth of dendrites is initiated by activating a thermoelectric cooler to chill a small isolated volume of the melt, leading to a dendritic mushy zone. Finally, the temperature is raised to remelt the crystals, returning the system to a stable melt phase.

We are particularly interested in the final component of the IDGE. After sufficient melting of the mushy zone had occurred, the remaining fragments consisted of isolated crystallites that resulted from partially melting dendritic side branches. Typically these were roughly prolate spherical in shape (see Figure 1). For the final minute of melting of a particular crystal, video data (filmed at 30 frames per second) was analysed to determine the aspect ratio at each time. For the examples presented by Glicksman and co-workers [22, 43], the aspect ratio of the needle-shaped crystals increased with time from about 7 at $t_{\mathrm{e}}-t=60$ s to 17 at $t_{\mathrm{e}}-t=10$ s, where $t_{\mathrm{e}}$ is the final melting time (also referred to as the extinction time). After about $t_{\mathrm{e}}-t=10$ s, the aspect ratio rapidly decreased, and appeared to approach unity as $t\rightarrow t_{\mathrm{e}}^{-}$, meaning that the crystals were spherical just before extinction.

In order to make analytical progress, Glicksman et al. [21] model the process with a one-phase quasi-steady problem, which results by ignoring heat conduction within the crystals and assuming an infinite Stefan number. Here, the Stefan number is defined by

[TABLE]

where $c$ is the specific heat, $L$ the latent heat of fusion per mass and $u_{\infty}^{*}-u_{\mathrm{m}}^{*}$ is the temperature difference between the melt away from the crystal and the melting temperature. In reality, for this particular experiment the parameter values were $L/c\approx 10.99$ K, $u_{\infty}^{*}-u_{\mathrm{m}}^{*}\approx 1.8$ K, so $\beta\approx 6.1$, which is not reasonably large. Glicksman et al. [21] derive an exact solution to the infinite-Stefan-number problem in an infinite domain in prolate spheroidal coordinates, which applies under the further assumption that the aspect ratio of the dendrite remains constant. This solution is a special case of that presented earlier by Ham [25] and Howison [28] (which holds for the more general shape of an ellipsoid with constant aspect ratios), for example, and that derived using the Baiocchi transform by McCue et al. [45] (again, for an ellipsoid). The solution was used by Glicksman et al. [21] to approximate the time-dependence of the melting process, with quite good agreement with experimental results.

Glicksman and co-workers [21, 22, 43] did not provide an explanation for the observed increase in aspect ratio during the first 50 s of melting; however, the subsequent decrease in aspect ratio (during the final 10 s of melting) was accounted for by noting that by this stage of the melting process the crystals had become small enough for surface tension effects to begin to dominate [22, 43]. As a consequence, the needle tips with high curvature melted more quickly than the remainder of the crystals, in accordance with the Gibbs-Thomson law

[TABLE]

which states that the actual melting temperature on a curved surface is not constant, but instead depends weakly on the mean curvature $\kappa^{*}$ (defined to be positive for a sphere) via the surface tension coefficient $\gamma$ (defined to be $\gamma=2\sigma^{*}/\rho_{s}L$, where $\sigma^{*}$ measures surface energy effects with dimensions Nm*-1* or Jm*-2* and $\rho_{s}$ is the density of the solid phase) [3]. Here $\partial\Omega^{*}$ denotes the solid-melt interface. For the IDGE experiments, the surface tension coefficient is roughly $\gamma\sim 10^{-10}\,$ m.

In this article, we are motivated by these issues to undertake a theoretical study of the one-phase quasi-steady Stefan problem. The mathematical problem is re-formulated in Section 2 with a Baiocchi transform for the special zero-surface-tension case. In Section 3, we go on to provide a near extinction analysis for a general shaped initial crystal, including numerical results for cases in which crystals ultimately melt to a single point or pinch off and break into two separate pieces. The role of surface tension is then explored in Section 4, while in Section 5 we consider an additional effect on the moving boundary, kinetic undercooling. We show that kinetic undercooling acts as a de-stabilising term, and is effectively in competition with surface tension. When these two terms are considered simultaneously, we find that the aspect ratio of a prolate spheroid can initially increase before decreasing suddenly to unity in the extinction limit, which is the same behaviour as observed in the IDGE. We close in Section 6 with a summary of the key results and a brief discussion of how our work relates to the experiments described by Glicksman and co-workers [21, 22, 43]. An important point to note is that the quasi-steady assumption used in this article leads to a moving boundary problem that also describes bubble contraction in a porous medium [12, 28, 45]. Thus our study also describes the effect that surface tension has on the shape of a bubble in the limit that it contracts to a point. This connection is revisited in Section 6.

2 Quasi-steady formulation with zero surface tension

2.1 Governing equations

Consider a solid substance (the crystal dendrite), initially at melting temperature $u_{\mathrm{m}}^{*}$ occupying the region $\Omega^{*}(0)$, surrounded by the same substance in liquid form in $\mathbb{R}^{3}\setminus\Omega^{*}$. In the far field, a higher temperature $u_{\infty}^{*}$ is applied, and thus melting proceeds until the crystal melts completely at the extinction time $t_{\mathrm{e}}^{*}$.

Setting $k$ to be the thermal diffusivity, we scale variables using

[TABLE]

where $\ell$ is a characteristic length scale of the initial crystal shape, and $\beta$ is the Stefan number (1). The resulting one-phase Stefan problem for melting the crystal is

[TABLE]

For what follows we shall take the quasi-steady limit $\beta=\infty$, which is an appropriate approximation for experiments in which the latent heat is large or the specific heat is small. As a result, the parabolic equation (4a) becomes Laplace’s equation

[TABLE]

and thus we do not require an initial condition for $u$.

As mentioned in the Introduction, the governing equations (4e) with (4b)-(4d) are also relevant for the problem of a bubble that is forced to contract in a saturated medium, where the fluid flow is governed by Darcy’s law [12, 28, 45], as well as the two-dimensional analogue for Hele-Shaw flow [15, 14, 42]. These equations also arise in other moving boundary problems, for example the small Péclet number limit of advection-diffusion-limited dissolution/melting models [6, 27, 32, 53, 57], for which it is also of interest to track the moving boundary and predict its shape and location (the collapse point [53]) close to the extinction time; other closely related advection-diffusion-like moving boundary problems in potential flow have similar governing equations in the small Péclet number limit [4, 7].

2.2 Baiocchi transform

We use the Baiocchi transform defined by

[TABLE]

where we are using the notation $t=\omega({\bf x})$ to denote the solid-melt interface $\partial\Omega$. The Baiocchi transform is widely used in the analysis of moving boundary problems with boundary conditions of the form (4b)-(4c), for example [8, 13, 31, 36, 40, 45]. Note that while here we restrict ourselves to (4e), the approach is also applicable to (4a) [44, 46].

Transforming the governing equations (4e) with (4b)-(4d), we derive the nonlinear moving boundary problem for $w$:

[TABLE]

Once a solution for the Baiocchi variable $w$ is determined, the temperature $u$ can be recovered via $u=\partial w/\partial t$. We note that an advantage of the Baiocchi transform is that it transforms the inhomogeneous boundary condition (4c) into a homogeneous boundary condition. Another is that time appears as a parameter in (6a)-(6e), so that the problem can be solved at any time without knowledge of the solution at previous times.

2.3 Exact solution for prolate spheroid

For the case in which the initial crystal shape $\partial\Omega(0)$ is an ellipsoid, (6a)-(6e) can be solved in ellipsoidal coordinates exactly, as done as part of the analysis by McCue et al. [45]. The solution for the interface $\partial\Omega(t)$ remains ellipsoidal with constant aspect ratios for all time until extinction. An equivalent solution without the Baiocchi transform is provided in Howison [28].

We present here a summary of this exact solution in the special case for which the initial crystal shape $\partial\Omega(0)$ is the prolate spheroid

[TABLE]

with initial aspect ratio $\mathcal{A}(0)=z_{0}(0)$. (This special case, together with the case in which the crystal is initially an oblate spheroid, is also recorded by McCue et al. [45].) The exact solution is that $\partial\Omega(t)$ retains its prolate spheroidal shape as

[TABLE]

where $z_{0}(t)>0$ and $\rho_{0}(t)>0$ measure the major and minor axes of the dendrite, respectively, with constant aspect ratio $\mathcal{A}(t)=z_{0}(t)/\rho_{0}(t)=z_{0}(0)$ (here the length scale $\ell$ is chosen so that $\rho_{0}(0)=1$). The full solution has the time-dependence

[TABLE]

where

[TABLE]

The result (10) is also derived in Glicksman et al. [21]. Although, as mentioned above, the aspect ratio of the melting crystals in the Isothermal Dendritic Growth Experiment was not constant, these authors make a rough guess for the average value of the aspect ratio over the first 50 seconds of melting, and then compare (9) with experimental results. Their agreement is quite good, reflecting the square root of time dependence near extinction.

3 Analysis of zero-surface-tension problem

McCue et al. [45] were concerned primarily with analysing the near extinction behaviour for a variation of (6a)-(6e) in which $\Omega(0)$ coincides with an outer boundary (i.e., a finite-domain problem in which the crystal initially occupies the entire domain). Here we provide equivalent results for the full infinite-domain problem (6a)-(6e) and apply the level set method to support these findings.

3.1 Extinction time and extinction points

For a given initial crystal shape $\Omega(0)$, we wish to determine how long it takes to melt (the extinction time $t_{\mathrm{e}}$) and the point at which the crystal eventually vanishes as $t\rightarrow t_{\mathrm{e}}^{-}$ (the extinction point ${\bf x}_{\mathrm{e}}$). The convenient framework for this analysis is via the Baiocchi transform. As mentioned above, time appears as a parameter in (6a)-(6e), meaning we can skip to the extinction time to compute $w_{\mathrm{e}}({\bf x})=w({\bf x},t_{\mathrm{e}})$. It is convenient to set $w_{\mathrm{e}}=W({\bf x})+t_{\mathrm{e}}$, so $W$ satisfies the linear problem

[TABLE]

The extinction point ${\bf x}_{\mathrm{e}}$ is then the local maximum of $W$, and the extinction time is recovered via $t_{\mathrm{e}}=-W({\bf x}_{\mathrm{e}})$. As noted by Entov & Etingof [15], (11a)-(11c) defines the dimensionless gravity potential of $\Omega(0)$, thus

[TABLE]

which provides an interesting connection between our problem and gravity potential generated by a uniform body.

Whilst in practice it is not feasible to compute $W$ analytically for a general initial crystal shape $\Omega(0)$, such a calculation can be performed numerically. Indeed, we provide a number of simple examples in Section 3.4 in which we compute $W$ for both convex and non-convex initial shapes. We include in those examples cases for which $W$ has two local maxima. In such instances, if the two local maxima are equal, then the crystal must pinch off into two, with the local maxima corresponding to the extinction points for each of the two satellite crystals. We also provide an example of the more complicated case in which there are two local maxima that are not equal; here, the use of $W$ can only predict the final extinction for the largest of the two satellite crystals.

3.2 Near-extinction analysis

For the case of an axially symmetric initial crystal with the $z$ axis pointing down the centreline, we can translate the coordinate system so that the extinction point ${\bf x}_{\mathrm{e}}$ lies on the origin. Since $w_{\mathrm{e}}=0$ at ${\bf x}={\bf x}_{\mathrm{e}}$ and ${\bf x}_{\mathrm{e}}$ is a local maximum of $w_{\mathrm{e}}$, a simple Taylor series for this axially symmetric geometry implies that $w_{\mathrm{e}}\sim a(x^{2}+y^{2})+bz^{2}$ as $r\rightarrow 0$. Further, as a consequence of (11b), we then have

[TABLE]

where $1/6<a<1/4$. As we shall see, the parameter $a$ is effectively all the melting crystal “remembers” from its initial condition; it is this single parameter that controls the aspect ratio of the crystal at extinction. Note that the higher order terms in (13) are not required in the following analysis (they would be for the special case $a=1/4$, which represents the borderline between the type of extinction considered in this section and when a bubble breaks up into two, as treated in Section 3.4).

In the limit $t\rightarrow t_{\mathrm{e}}^{-}$, the inner region is for $r=\mathcal{O}(T)$, where $T(t)$ is a length scale defined so that the volume of the melting crystal is fixed to be $4\pi T^{3}/3$. We write $w\sim T^{2}\Phi({\bf X})$ as $t\rightarrow t_{\mathrm{e}}^{-}$, where ${\bf X}={\bf x}/T$, so that

[TABLE]

as $R\rightarrow\infty$, where $d$ is a constant found as part of the solution to (14a)-(14c). We see from (14c) that a matching condition for the outer region is

[TABLE]

The solution to (14a)-(14c) in prolate spheroidal coordinates is provided in Appendix A. According to this solution the dendrite boundary $\partial\Omega_{0}$ is described by

[TABLE]

where $q_{0}$ is a parameter that is related to the special constant $a$ by

[TABLE]

Further, the constant $d$ in (14c) is related implicitly to $a$ by

[TABLE]

Note that the prolate spheroid approaches a perfect sphere in the limit $a\rightarrow 1/6^{+}$, in which case $d\rightarrow 1/2^{+}$.

The outer region is for $r=\mathcal{O}(1)$, for which

[TABLE]

Matching with the inner gives the time-dependence

[TABLE]

or, in other words,

[TABLE]

Thus we see that, regardless of the shape of the initial crystal, the square root of time scaling determined experimentally in Glicksman et al. [21] is as expected.

In summary, the zero-surface-tension model predicts that, provided there is no pinch-off, an axially symmetric dendrite will melt to a spheroid in the extinction limit. While this spheroid could be prolate or oblate, we concentrate here on the prolate case, as this is the one observed in the IDGE [21, 22, 43]. The aspect ratio of the prolate spheroid at extinction is given by

[TABLE]

which provides an implicit dependence of $\mathcal{A}$ on the constant $a$ via (17). Here $a$ is the only parameter that is required to characterise the initial dendrite shape (it is found by solving (12) and expanding $w_{\mathrm{e}}$ about ${\bf x}_{\mathrm{e}}$). The time-dependence of the melting is given by (21), where the volume of the dendrite shrinks like $4\pi T^{3}/3$ (in other words, $T$ provides a natural length scale for the melting dendrite). Again, this time-dependence is related to the initial dendrite shape via the parameter $a$ (since $d$ is given by $a$ through (18) and (17)).

In the special case in which the dendrite is initially the prolate spheroid (7), then it retains its aspect ratio. This is, of course, the exact solution listed in Section 2.3.

Finally, for sufficiently symmetric crystals we have $a=1/6$ which gives $d=1/2$. Here $\Phi=R^{2}/2-1/2+1/3R$ and the dendrite becomes spherical in the limit with $T\sim\sqrt{2}(t_{\mathrm{e}}-t)^{1/2}$. The special case of an initially spherical dendrite remains spherical.

At this point it is worth mentioning that for large Stefan numbers, $\beta\gg 1$, the scaling (21) eventually ceases to hold for the full classical Stefan problem with (4a) instead of (4e) [46]. However, this discrepancy would not be observed on the scale of the IDGE experiments.

3.3 Null quadrature domains

It is worth relating some of the above arguments to well-known and long-established results [12, 18, 28]. First, by applying Green’s theorem it can be shown that

[TABLE]

where $\Phi$ is a suitable harmonic function and $\Omega(t)$ is the shape of a melting crystal from the infinite-domain problem (6a)-(6e) (Howison [28]). Noting that these quasi-steady problems with zero surface tension are time-reversible, we can seek so-called ‘ancient’ solutions for which the entire domain $\mathbb{R}^{3}\setminus\Omega(t)$ vanishes in the limit $t\rightarrow-\infty$. From (23) it follows that for these ancient solutions $\mathbb{R}^{3}\setminus\Omega(t)$ must be a null quadrature domain. The only suitable such domain is the exterior of an ellipsoid (see Karp [34] for a discussion on null quadrature domains). For any other initial crystal shape, the backwards problem with $t$ decreasing leads to some kind of finite-time blow-up or perhaps a scenario in which part of the crystal boundary expands infinitely leaving behind ‘fjords’ or ‘tongues’ (these scenarios are much better understood in the two-dimensional Hele-Shaw problem; see also Howison [29, 30] for explicit examples of each case).

As discussed in Section 3, for a melting crystal (of general initial shape) the generic limiting behaviour is that it becomes ellipsoidal in shape as $t\rightarrow t_{\mathrm{e}}^{-}$. This result can also be derived using an alternative approach, as suggested more recently by King & McCue [36], who treated the two-dimensional Hele-Shaw case. First, we see that for the integral in (23) to converge we could choose $\Phi=r^{\ell}Y_{\ell}^{m}$, where $Y_{\ell}^{m}$ are spherical harmonics and $\ell$ is an integer such that $\ell\leq-4$. Rescaling lengths such that $\bar{r}=r/T$, we have from (23) that

[TABLE]

Thus, the left-hand side vanishes as $T\rightarrow 0$, or $t\rightarrow t_{\mathrm{e}}^{-}$, meaning that the exterior of the crystal approaches a null quadrature domain in the limit, and thus the crystal itself approaches an ellipsoid in shape.

3.4 Numerical examples

We present some numerical examples that demonstrate the key features discussed above. To solve (11a)-(11c) numerically, we formulate a level set function, $\phi(\mathbf{x})$, such that $\phi>0$ for $\mathbf{x}\in\Omega(0)$ and $\phi<0$ for $\mathbf{x}\in\mathbb{R}^{3}\backslash\Omega(0)$. Thus we can reformulate (11a) and (11b) as

[TABLE]

where $H$ is the Heaviside function. We note that $H(\phi)$ is discontinuous at $\mathbf{x}\in\partial\Omega(0)$, so for numerical purposes we implement a smoothed Heaviside function

[TABLE]

where $\delta=1.5\Delta x$. For this purpose, it is convenient to work in spherical polar coordinates $(r,\theta,\varphi)$ and represent the axially symmetric moving boundary $\partial\Omega$ by $r=s(\theta,t)$. Thus, (25) becomes

[TABLE]

The spatial derivatives in (3.4) are approximated using central finite differencing, with homogeneous Neumann boundary conditions applied at $r=0$, $\theta=0$, and $\theta=\pi$. The far-field boundary condition (11c) is incorporated using a Dirichlet-to-Neumann map described in Section B.2.2.

3.4.1 Symmetric initial condition

We consider a selection of initial conditions to illustrate a few different qualitative behaviours. Again, using spherical polar coordinates $(r,\theta,\varphi)$ with $\partial\Omega$ denoted by $r=s(\theta,t)$, the first is the prolate spheroid

[TABLE]

where $r_{0}$ describes the initial aspect ratio. The second initial condition is a peanut-shaped interface described by

[TABLE]

where $r_{0}$ can be interpreted as a measure of the depth of the pinch in the middle of the peanut. Following Garzon et al. [19], the third initial condition considered is a dumbbell shaped interface of the form $s(\theta,0)=(\rho^{*}(\theta)^{2}+z^{*}(\theta)^{2})^{1/2}$, where

[TABLE]

for $0\leq\theta\leq\pi/2$; for $\pi<\theta\leq 2\pi$ this initial condition is made symmetric by reflecting about $\theta=\pi/2$.

In Figure 2, we illustrate some numerical results by choosing parameter values from these three initial conditions. For the prolate spheroid (27) we provide results for $r_{0}=0.8$, noting that this initial condition is obviously convex. For the peanut shaped surface (28), we choose $r_{0}=0.5$, which is not convex but is instead mean convex. Finally, for the dumbbell shape (29a)-(29d), we choose $r_{0}=4.75$, which again corresponds to a nonconvex shape which is still mean convex, but this time with a particularly thin neck region. In all of these case, we show in Figure 2 the initial shape, the numerical solution to (4b)-(4e) shortly before the extinction time, and the corresponding solution to the Baiocchi transform problem (11c) and (3.4).

For both of the first two examples in Figure 2, namely (27) with $r_{0}=0.8$ and (28) with $r_{0}=0.5$, the solution to (4b)-(4e) contracts to a single point at extinction. By observing the third column of Figure 2, we see this is consistent with the solution of (11c) and (3.4), which shows $|W|$ having one local maximum at the origin, predicting one point at extinction. This comparison highlights that convex shapes and some nonconvex shapes will contract to a single point. The extinction time predicted by the Baiocchi transform is computed by evaluating $|W|$ at ${\bf x}_{\mathrm{e}}$ (which, for this problem is the origin) giving the values $t_{\mathrm{e}}=0.370$ and $t_{\mathrm{e}}=0.233$ for (27) with $r_{0}=0.8$ and (28) with $r_{0}=0.5$, respectively. Comparing this to the extinction times computed from the numerical solution to (4b)-(4e), we find there is less than $0.1\%$ relative difference, suggesting excellent agreement.

The equation (28) with $r_{0}=0.5$ provides a good test for the prediction (22). For this purpose we take the solution to the Baiocchi transform problem (11), which in this case predicts that $q_{0}=1.100$ and $a=0.215$. As such, our prediction for the aspect ratio at extinction is $\mathcal{A}=2.395$. The time-dependent behaviour of the aspect ratio for our numerical solution to the full problem (using the level set method) is presented in Figure 3. This figure demonstrates how well these two results agree with other.

For initial condition (29a)-(29d) with $r_{0}=4.75$, Figure 2 shows different qualitative behaviour. Here, we see that solutions to (4b)-(4e) will undergo pinch-off and ultimately the two satellite crystals will contract to separate points of extinction. Again, this is consistent with the solution to (11c) and (3.4) as the third column of Figure 2 indicates that $|W|$ has two local maxima. By approximating the locations of these maxima and the values of $|W|$ at these points, we find the Baiocchi transforms predicts that the interface will contract to extinction points at $z_{\mathrm{e}}=\pm 0.577$ at time $t=0.100$. Comparing these results with the extinction locations and times approximated from the numerical solution to (4b)-(4e), we find a relative difference less than $0.2\%$. This example shows, for symmetric initial conditions, how well the Baiocchi transform approach can be used to predict whether pinch-off will occur, as well as the extinction points and time.

In summary, these numerical results indicate that for a given initial interface, $\partial\Omega(0)$, each of the aspect ratio at extinction, the extinction time and location of the extinction point for an interface evolving according to (4b)-(4e) can be predicted from the solution to (11a)-(11c). Further, the indication is that this is true both for interfaces that contract to a single point of extinction, or undergo pinch-off and contract to multiple points of extinction, at least for symmetric initial conditions. Finally, these results illustrate the capacity of the level set based numerical scheme, presented in Appendix B, to accurately describe the dynamics of the interface once a change in topology has occurred.

3.4.2 Asymmetric initial condition

The numerical solutions of (4b)-(4e) presented in Section 3.4.1 indicate that when $\partial\Omega(t)$ is sufficiently non-convex then the interface will undergo a change in topology. As initial conditions considered in Section 3.4.1 are symmetric along the major axis (about $\theta=\pi/2$), the two interfaces which form after pinch-off will have the same extinction time. We now investigate a class of asymmetric initial conditions that undergo pinch-off into two surfaces of differing volumes. We expect the smaller of the two volumes to contract to a point first, followed by the larger, thus giving two distinct extinction times.

We again consider an initial condition of the form of (29a)-(29d), but this time for $0\leq\theta\leq\pi$. In Figure 4, we plot the time evolution of the numerical solution to (4b)-(4e) and the corresponding numerical solution to (11a)-(11c) for the representative case $r_{0}=5.1$. We observe that the full time-dependent solution to (4b)-(4e) undergoes a change in topology at approximately $t=0.076$, with crystal domain $\Omega(t)$ pinching off into two. The smaller of the two satellite crystals contracts to a point at $z_{\mathrm{e}}=0.564$ when $t=0.086$, followed by the remaining larger satellite crystal which contracts to a point at $z_{\mathrm{e}}=-0.773$ when $t=0.127$. The corresponding numerical solution to the Baiocchi transform problem (3.4), Figure 4 indicates that $|W|$ has two local maxima, located at $z_{\mathrm{e}}=0.506$ and $z_{\mathrm{e}}=-0.767$, with $|W|$ equal to 0.092 and 0.127 at these points, respectively. Thus we see that the predicted values of the extinction points and times agree well for the larger of the two satellite crystals (as it should) but not at all for the smaller crystal. redThat our approach can only provide information about the extinction time and point for the largest satellite crystal is a minor limitation to the Baiocchi transform framework.

4 Effects of surface tension

An inevitable consequence of melting a small crystal is that eventually the curvature will become large enough so that surface tension effects become important. For what follows, instead of (4b) we use the dimensionless version of the Gibbs-Thomson law (2), which is

[TABLE]

where $\sigma=\gamma u_{m}^{*}/\ell(u^{*}_{\infty}-u^{*}_{m})$ is the dimensionless surface tension coefficient, and $\kappa$ is the dimensionless signed mean curvature.

4.1 Linear stability analysis for near spherical crystal

It proves useful to outline the linear stability analysis for interfaces evolving according to (4c)-(4e) and (30) with a near-spherical initial condition. In spherical polar coordinates $(r,\theta,\varphi)$, we represent the axially symmetric moving boundary $\partial\Omega$ by $r=s(\theta,t)$, so that our problem is

[TABLE]

We seek a perturbed spherical solution to (31a)-(31d) of the form

[TABLE]

where $\varepsilon\ll 1$. The leading order solution is

[TABLE]

where $s_{0}(0)=r_{0}$. For the $\mathcal{O}(\varepsilon)$ system,

[TABLE]

the solutions are of the form

[TABLE]

where $A_{n}$ is a sequence of unknown coefficients, $P_{n}$ is the $n$th Legendre polynomial, and $\gamma_{n}$ is the $n$th mode of perturbation to the sphere. We are able to eliminate $A_{n}$ to obtain

[TABLE]

Since $(1/\gamma_{n})\mathrm{d}\gamma_{n}/\mathrm{d}s_{0}\to 0$ in the limit that $s_{0}\to 0$ for $n\geq 2$, we see that each mode of perturbation is stable, and a perturbed sphere will evolve to a sphere in the extinction limit, as expected.

The special case in which the initial condition is the prolate spheroid with major and minor axes $r_{0}+\varepsilon$ and $r_{0}$, respectively, then

[TABLE]

That is, $\gamma_{2}(0)=2/3$ and $\gamma_{n}(0)=0$ for $n\geq 3$. This initial condition has an aspect ratio of $1+\varepsilon/r_{0}+\mathcal{O}(\varepsilon^{2})$. The exact solution for $\gamma_{2}$ is

[TABLE]

and the aspect ratio for this particular initial condition therefore becomes

[TABLE]

Note that when $\sigma=0$, then $3\gamma_{2}/2s_{0}=1/r_{0}$, resulting in the aspect ratio remaining constant, which is consistent with the known exact solution of Section 2.3. Otherwise, for $\sigma>0$, the aspect ratio decreases monotonically to unity, as shown later in Figure 5.

4.2 Long thin needle problem

We consider here the limit of a long thin melting dendrite. Suppose the axially-symmetric shape of the dendrite is given by $\rho=S(z,t)$ where $\rho^{2}=x^{2}+y^{2}$. Suppose also that $S_{0}(z)=S(z,0)$, $\rho_{0}(t)=S(0,t)$, $S(z_{0}(t),t)=0$, where $\alpha=z_{0}(0)/\rho_{0}(0)\ll 1$ such that the initial aspect ratio, $\mathcal{A}(0)=1/\alpha$, is large.

The inner region is for $r=\mathcal{O}(\rho_{0}(t))$. Here the melting is almost two-dimensional with $\partial u/\partial z\ll 1$ and $\partial S/\partial z\ll 1$ so that, to leading order,

[TABLE]

The solution to (40a)-(40c) is

[TABLE]

where the form for $S$ is determined by the missing far-field condition, which is found by considering the outer region.

In this outer region, which is for $r=\mathcal{O}(z_{0}(t))$, the dendrite appears as a slit. We scale $\tilde{\rho}=\rho/(\alpha\rho_{0}(t))$, $\tilde{t}/\ln\alpha$ and rewrite the inner solution (41) to be

[TABLE]

The leading order solution in the outer region is $u=1$, thus, after matching with the leading order term in (42) as $\alpha\to\infty$, we find

[TABLE]

For the zero surface tension case $\sigma=0$, we can solve (43) explicitly to give

[TABLE]

again providing square root time dependence.

Of particular interest is the special case in which the initial dendrite is the prolate spheroid (7). Here $\rho_{0}=\alpha$ and $z_{0}(0)=1$, so initially the dendrite has the aspect ratio $\mathcal{A}(0)=1/\alpha$. From (43) we find the interface is given implicitly by

[TABLE]

Note that the small parameter in this limit is $1/\ln\alpha$, which suggests the analysis here is valid only for extremely large aspect ratios.

4.3 Numerical results for canonical problem

For the melting prolate spheroidal crystal considered in Section 2.2, whose surface is (8), we find the mean curvature is largest near the tip, given by

[TABLE]

Thus the right hand side of (40b) becomes $\mathcal{O}(1)$ when $t_{\mathrm{e}}-t=\mathcal{O}(\sigma^{2})$, suggesting we rescale according to

[TABLE]

and treat the following problem

[TABLE]

when $\hat{t}=\mathcal{O}(1)$, $|\mathbf{\hat{x}}|=\mathcal{O}(1)$, where hats denote scaled quantities. For the case in which the initial crystal, $\hat{\Omega}(0)$, is a prolate spheroidal in shape, this is a canonical problem for melting a solid. This one parameter in the problem is the initial aspect ratio.

Using the numerical scheme described in Appendix B, we solve (48a)-(48d) for $\hat{u}$ and $\hat{\Omega}$. We first consider a near spherical prolate spheroid initial condition such that the initial aspect ratio is close to unity. Figure 5 compares the aspect ratio of the numerical solution to (48a)-(48d) with $\alpha=0.85$ with the aspect ratio as predicted by linear stability analysis given by (39). This figure shows excellent agreement between the numerical solution and linear stability analysis, confirming that the numerical scheme presented in Appendix B is able to describe the behaviour of the interface as the aspect ratio decreases to unity. Further, we numerically solve (48a)-(48d) with $\alpha=1/6$, and plot the time evolution of the solution and corresponding aspect ratio in Figure 6. As expected, this figure shows that the aspect ratio decays to unity in the limit that $t\to t_{\mathrm{e}}^{-}$.

5 Kinetic undercooling

In this section, we very briefly consider the effects of extending the dynamic boundary condition (30) to include a kinetic undercooling-type term:

[TABLE]

where $v_{n}$ is the normal velocity of $\partial\Omega$ and $c$ is the kinetic coefficient. An argument for this extended boundary condition is that (30) can be derived under equilibrium conditions, while (49) is a corrected version that takes into account nonequilibrium kinetic effects [24, 41]. Physically, a nonzero kinetic coefficient $c>0$ penalises high interface speeds, which is important near extinction since our interface speed scales like $(t_{\mathrm{e}}-t)^{-1/2}$. A wide variety of studies of Stefan problems have considered kinetic undercooling [2, 3, 10, 11, 16, 17, 35]. The other important previous study is Dallaston & McCue [9], where the two-dimensional analogue of the quasi-steady problem (4e), (49), (4c)-(4d) is treated in some detail.

Following the linear stability analysis outlined in Section 4.1 using (49) with $c>0$, we find the second mode of perturbation satisfies

[TABLE]

from which we see that

[TABLE]

suggesting that an initially prolate spheroidal crystal will tend to a sphere in the extinction limit. This conclusion is that same as before in Section 4.1 when $c=0$. On the other hand, a significant difference in qualitative behaviour is that the aspect ratio with $c>0$ may first increase and then decrease (to unity), which is a feature not observed when $c=0$. The turning point can be calculated via

[TABLE]

Given $s_{0}$ is defined on the domain $0\leq s_{0}\leq r_{0}$, the aspect ratio will monotonically decrease to unity if

[TABLE]

otherwise, the aspect ratio will be non-monotone.

Our work is motivated in part by a series of experiments performed as part of the IDGE [21, 22, 43]. In these experiments, it was observed that the aspect ratio of melting crystals increased for a period of time before decreasing to unity at extinction. In the context of the results presented in this section, Figure 7 illustrates the aspect ratio of a (near-spherical) prolate spheroid predicted by linear stability analysis and the aspect ratio of the melting PVA crystals [22]. This figure shows that when both the effects of surface tension and kinetic undercooling are considered, the solution to (4c)-(4e) and (49) is qualitatively similar to the experimental results (while of course the scale is different).

6 Discussion

In this paper, we have studied a quasi-steady one phase Stefan problem for melting an axially symmetric crystal. In Section 3 we treat a zero-surface-tension model and use analytical tools to show that axially symmetric crystals will tend to prolate spheroids in the limit that they melt completely, namely $t\to t_{\mathrm{e}}^{-}$, with an aspect ratio that depends on the initial condition. The point to which the crystals ultimately shrink, together with the melting time, is predicted by this analysis and confirmed using a novel numerical scheme based on the level set method (presented in Appendix B). An advantage of this scheme is that we are also able to present numerical results for crystals that undergo pinch-off and contract to multiple points of extinction.

We consider the effects of surface tension by the Gibbs-Thomson law (2) in Section 4. By performing linear stability analysis on the spherical solution, we show that surface tension acts to smooth out perturbations to the interface, suggesting it becomes spherical in the extinction limit. A numerical study of canonical problem confirms this prediction. These results are as expected and also indicated by the experimental results summarised by Glicksman and co-workers [21, 22, 43]. However, the one feature of the IDGE not described by the model with surface tension is the non-monotonic behaviour of the aspect ratio, where the aspect ratio first increases as the crystal becomes very long and thin, and then very quickly decreases to unity as surface tension ultimately acts to produce a perfect sphere in the extinction limit. In order to mimic this non-monotonic behaviour, we have included the effects of kinetic undercooling in the model in Section 5, which shows that the competition between kinetic effects and surface tension does indeed produce the qualitative behaviour observed.

A key assumption in our paper is that the Stefan number in (4a) is taken to be large, namely $\beta\gg 1$, so that (4a) reduces to (4e) and our moving boundary is therefore quasi-steady. There are two issues related to this assumption that we wish to mention. First, our problem for melting a crystal is the same as that for a bubble contracting in a porous medium where the flow is governed by Darcy’s law [12, 28, 45], although in that context the far-field (Dirichlet-type) boundary condition (4d) should probably be replaced with a flux condition that dictates how quickly the bubble volume is decreasing (in two dimensions the equations describe bubble contraction in a Hele-Shaw cell [15, 14, 42]). For the case in which a bubble pinches off to produce two shrinking bubbles, the problem formulation would also need to consider two points of extraction that coincide with the eventual extinction points. The second issue is that, strictly speaking, for the extremely late stages of melting, our quasi-steady model with (4e) is no longer applicable in the large Stefan number limit, and instead (4a) must be retained. The mathematical details of such an exponentially short final-melting stage have been recorded in a number of previous studies [1, 26, 44, 46, 55].

Acknowledgments

SWM and LCM acknowledge the support of the Australian Research Council Discovery Project DP140100933. We are grateful to the anonymous referees for their helpful feedback.

Appendix A Prolate spheroids with constant aspect ratio

To solve the inner problem (14a)-(14c) we employ prolate spheroidal coordinates $(\xi,\eta,\phi)$ defined by

[TABLE]

where $\xi\geq 0$, $0\leq\eta\leq\pi$, $0\leq\phi<2\pi$, and $k$ is a constant to be determined below. The crystal boundary $\partial\Omega_{0}$ is described by $\xi=\xi_{0}$ or, equivalently,

[TABLE]

Motivated by the relationship

[TABLE]

we look for a solution of the form

[TABLE]

where $q=\cosh\xi$ and $q_{0}=\cosh\xi_{0}$ and obtain a coupled system of two second order (Legendre-type) differential equations for $f_{1}$ and $f_{2}$. These (and the constant $k$) are solved subject to the four conditions $f_{1}=f^{\prime}_{1}=f_{2}=f^{\prime}_{2}=0$ on $q=q_{0}$, and the far-field condition (14c) to give

[TABLE]

where

[TABLE]

and $d$ is given by (18). The important relationship between $q_{0}$ and the special constant $a$ is given by (17).

Appendix B Numerical solution - A level set approach

To find the numerical solution of (4b)-(4e), we implement a level set based approach. The level set method (LSM), first proposed by Osher and Sethian [52], is a tool used to study a wide range of moving boundary problems. We refer the reader to Osher & Fedkiw [51] and Sethian [54] for comprehensive overviews of implementation strategies and applications. The LSM utilises an Eulerian approach by representing an $n$-dimensional interface, $\partial\Omega(t)$, as the zero level set of a $n+1$-dimensional surface, $\phi(\mathbf{x},t)$, such that

[TABLE]

By representing the interface implicitly, the LSM can be used to describe complex behaviour such as the changes in topology observed in Figure 2, while operating on a simple regular two-dimensional grid.

The evolution of the level set function $\phi$ is described by the level set equation

[TABLE]

where $F$ is a continuous function defined on all of the computational domain, satisfying $F=V_{n}$ on $\mathbf{x}=\partial\Omega(t)$. In the context of (4b)-(4e), by noting that the outward normal of $\phi$ is $\mathbf{n}=\nabla\phi/|\nabla\phi|$, a suitable expression for $F$ on and outside the interface is

[TABLE]

This leaves the matter of defining a suitable extension of $F$ to inside the interface. Among several possibilities in the level set literature, we opt for a biharmonic extension as proposed by Moroney et al. [49], and compute $F$ inside the interface to satisfy

[TABLE]

together with the boundary conditions that $F$ and $\partial F/\partial n$ are continuous across $\partial\Omega(t)$. This method of extension shares the main property of the LSM itself, in not requiring the location of the interface to be calculated explicitly. To solve (64), we formulate the biharmonic stencil over the entire domain, which is then modified so that values of $F$ outside the interface, whose location is determined from the sign of $\phi$, are not overwritten. The resulting linear system is solved using LU decomposition. This extension is a variant of a two-dimensional thin plate spline interpolant defined on the level set grid.

B.1 General algorithm

The algorithm used to solve (4b)-(4e) numerically is outlined as follows:

• Step 1

For a given initial condition $s(\theta,0)$, construct a level set function $\phi(r,\theta,0)$ such that $\phi<0$ inside the interface and $\phi>0$ outside the interface. This function is then converted to a signed distance function using the method of crossing times as described by Osher & Fedkiw [51].

• Step 2

Compute the temperature, $u$, on the domain $r\geq s(\theta,t)$ using the procedure described in Section B.2.

• Step 3

Compute $F$ according to (63), where the derivatives are evaluated using central finite differences. $F$ is extended over the entire computational domain by solving (64) at nodes where $\phi<0$, with boundary data from step 3.

• Step 4

Update $\phi$ by advancing the level set equation given by (62), where the time step is $\Delta t=0.25\times\Delta x/\max|F|$. We discretise the spatial derivatives in (62) using a ENO2 scheme for the spatial derivatives and integrate in time using second order Runge-Kutta where $\Delta t=0.25\times\Delta r/\max|F|$.

• Step 5

Reinitialise $\phi$ every 5 time-steps to a signed distance function by solving the reinitialisation equation

[TABLE]

where

[TABLE]

We use 5 pseudo-timesteps with $\Delta\tau=0.2\Delta r$.

• Step 6

Repeat steps 2-5 until the desired simulation time is attained.

B.2 Solving for temperature

Evaluating the speed function $F$ in the level set equation (62) requires first calculating the temperature $u$. This is achieved by using a modified finite difference stencil for Laplace’s equation in the region outside the interface. For nodes away from the interface, a standard 5-point stencil is used such that the discrete equation is

[TABLE]

For the singularity at $\theta=0$, noting that $\partial u/\partial\theta=0$ and using L’Hoptial’s rule then

[TABLE]

The same procedure is applied at $\theta=\pi$. Difficulties arise when attempting to incorporate the dynamic condition (4b) on the interface and the far-field boundary condition (4d). We detail the methodology used to overcome each of these difficulties in Sections B.2.1 and B.2.2, respectively. A schematic of the problem is given in Figure 8, which illustrates the different equations to be solved in each part of the computational domain.

We note that since the governing equation for temperature satisfies Laplace’s equation, an alternative approach for computing the temperature $u$ is the boundary integral method, which can be coupled with the level set method to solve problems where changes in topology occur [19]. However, an advantage of using a finite difference stencil is that it can easily be adapted to problems where the boundary integral method is not applicable. For example, we have used a similar method to the one presented in this section to study non-standard Hele-Shaw flow where pressure is not harmonic and for which the boundary integral method is much less suitable [50].

B.2.1 Incorporating the dynamic boundary condition

Special consideration must be taken when solving for nodes adjacent to the interface as we can no longer use the second order central differencing scheme (67). Instead we follow the work of Chen et al. [5] and approximate the spatial derivatives by fitting a quadratic polynomial from values on and near the interface and differentiating this polynomial twice. Supposing the interface is located between two nodes $(i-1,j)$ and $(i,j)$, the quadratic is fitted using the three points $(r_{b},u_{b})$, $(r_{i,j},u_{i,j})$, and $(r_{i+1,j},u_{i+1,j})$. Here $r_{b}$ denotes the location of the interface and $u_{b}$ is the temperature at the interface. The value of $r_{b}$ is found by noting that $\phi$ is a signed distance function and so the distance between $r_{b}$ and $r_{i,j}$, denoted $h$, can be calculated by

[TABLE]

Thus (67) becomes

[TABLE]

The same procedure is applied if the interface is between $r_{i}$ and $r_{i+1}$, or in the azimuthal direction.

The value of $u_{b}$ is determined by the dynamic condition (49), where in the case of surface tension the mean curvature term

[TABLE]

is approximated using central finite differences, while the normal velocity from the previous time step is used for the kinetic undercooling term.

B.2.2 The far-field condition

Special consideration must also be given when considering the boundary condition at $r\to\infty$. One method for simulating this far-field condition is to make the computational domain much larger than the radius of the interface and then impose $u=1$ on the outer boundary. However, this is computationally expensive as very large domains must be used to form an accurate solution. Instead, we simulate the far-field condition using the Dirichlet-to-Neumann (DtN) method [20]. This method is implemented by introducing a spherical artificial boundary, $R$, which is larger than the radius of the interface, i.e. $R>s(\theta,t)$. Outside of this boundary

[TABLE]

holds, where $f(\theta)$ is an unknown function. This problem can be solved exactly via separation of variables giving

[TABLE]

where

[TABLE]

and $P_{n}$ denotes the $n$th Legendre polynomial. Matching this outer solution with the inner numerical solution on the artificial boundary $R$ provides the necessary Neumann boundary conditions for the numerical scheme. By taking the derivative of (72) with respect to $r$ at $r=R$ and evaluating (73) using the trapezoidal rule, the finite difference stencil for the radial derivatives is updated with

[TABLE]

where

[TABLE]

and

[TABLE]

From a practical perspective, we cannot, of course, evaluate the series in (76) using an infinite number of terms, but have found that using 10 terms gives sufficient accuracy. Furthermore, it is a straightforward exercise to use the DtN method for other types of far-field boundary conditions such as flux condition for fluid flow whereby $\partial u/\partial r\sim 1/r^{2}$ as $r\to\infty$.