## Abstract

The classical Stefan problem for freezing (or melting) a sphere is usually treated by assuming that the sphere is initially at the fusion temperature, so that heat flows in one phase only. Even in this idealized case there is no (known) exact solution, and the only way to obtain meaningful results is through numerical or approximate means. In this study, the full two-phase problem is considered, and in particular, attention is given to the large Stefan number limit. By applying the method of matched asymptotic expansions, the temperature in both the phases is shown to depend algebraically on the inverse Stefan number on the first time scale, but at later times the two phases essentially decouple, with the inner core contributing only exponentially small terms to the location of the solid–melt interface. This analysis is complemented by applying a small-time perturbation scheme and by presenting numerical results calculated using an enthalpy method. The limits of zero Stefan number and slow diffusion in the inner core are also noted.

## 1. Introduction

One of the most simple moving boundary problems to pose is the classical Stefan problem for the inward solidification of a spherical ball of liquid or, equivalently, the melting of a spherical ice ball. For the solidification problem, consider a sphere of liquid of radius *a*, initially at some constant temperature *V*^{*}, which is greater than or equal to the fusion temperature . Now suppose the outer surface of the sphere is held at a temperature , which is lower than the fusion temperature; the sphere will begin to solidify inwards. The problem is to solve for the temperature field in both the liquid and the solid phases, as well as to track the location of the interface between the two phases. The corresponding melting problem is equivalent (with appropriate trivial adjustments in language), but for definiteness we think of the Stefan problem in terms of solidification, and note that all the results and conclusions in this paper apply equally well for both the cases.

The characterizing feature of this solidification problem, which holds for all Stefan problems, is the moving solid–melt boundary. Latent heat will be liberated at this interface between the two phases, and for the classical Stefan problem, it is assumed to be removed via conduction only. Thus, the model considered here is that the temperature in each phase is governed by the linear heat equation and coupled via a boundary condition (the Stefan condition) that describes the dependence of the latent heat removal on the speed of the interface.

We denote the thermal diffusivity in the solid and liquid phases by and , respectively. By scaling all temperatures, lengths and time with respect to , *a* and , respectively, and by setting the non-dimensional fusion temperature to be zero, the dimensionless solidification problem becomes(1.1)(1.2)with the following fixed boundary conditions:(1.3)(1.4)(1.5)(1.6)the moving boundary condition (the Stefan condition)(1.7)and the initial conditions (for a derivation of (1.7), see Gupta (2003) for example)(1.8)Here, *u*(*r*, *t*) and *v*(*r*, *t*) are the temperature fields in the solid and liquid, respectively, *r* is the radial distance, *t* represents time and *r*=*R*(*t*) describes the location of the solid–melt interface. The three parameters in the problem are the dimensionless initial temperature *V*, the Stefan number and the ratio of thermal diffusivities . Here, *c* is the specific heat of the substance, which is assumed to be the same constant in both the phases, and *L* is the latent heat of fusion. Note that we have neglected any density change upon solidification in our study.

The two-phase problem (1.1)–(1.8) is highly nonlinear with no known exact analytical solution. The overwhelming majority of the literature pertaining to this is concerned with the special case in which the liquid phase is initially at the fusion temperature, meaning *V*=0. In this case, *v*≡0 for all time and the heat flows in the solid phase only. The resulting one-phase problem (1.1), (1.3), (1.4), (1.7) (with *v*≡0) and (1.8)_{2} still does not have an exact solution, although a great deal of analytical progress has been made via large Stefan number asymptotics (Pedroso & Domoto 1973; Riley *et al*. 1974; Stewartson & Waechter 1976; Soward 1980), small-time analysis (Davis & Hill 1982; Hill & Kucera 1983) and near-complete solidification asymptotics (Soward 1980; Herrero & Velázquez 1997). Furthermore, numerous numerical treatments of the one-phase problem have been recorded, such as Tao (1967), Selim & Seagrave (1973) and Liu & McElwain (1997).

A considerably lesser time has been spent dealing with the full two-phase problem (1.1)–(1.8). The upper and lower bounds on the temperature and the location of the solid–melt interface are derived in Dewynne & Hill (1986) with the use of an integral formulation, and these are compared with the numerical results obtained by applying an enthalpy method (Voller & Cross 1981). Further analytical treatment is given in Kucera & Hill (1986), where an attempt is made to generalize the small-time analysis of Davis & Hill (1982; which is for the single-phase case only). For the corresponding problem in a cylindrical domain, a large Stefan number analysis is detailed in Jiji & Weinbaum (1978). We note that variations of (1.1)–(1.8) have also been considered in the literature. For example, with the Dirichlet condition (1.3) replaced by a Neumann condition or a radiation–convection-type condition, analytical results valid for small time are compared with numerical solutions in Gupta (1987) and Gupta & Arora (1992).

In §2, we extend the previous analyses dedicated to large Stefan number asymptotics (*β*≫1), which (with the exception of Jiji & Weinbaum 1978) were all for the idealized one-phase problem, to allow for the second phase. We find that, in general, the solid–melt interface *r*=*R*(*t*) moves slowly (when compared with the rate at which heat diffuses), with the process in the solid phase being roughly quasi-steady. This is because, since *β*≫1, a large amount of latent heat is being produced at the interface, and this heat needs to be diffused away for the evolution of interface. On the other hand, the temperature in the liquid core is found to decay relatively quickly to zero (the dimensionless fusion temperature), with the magnitude of the temperature gradient ∂*v*/∂*r* on the interface being much less than the magnitude of ∂*u*/∂*r*. Thus, at least for *β*≫1, the two phases are only weakly coupled via (1.7); the degree of coupling is quantified through matched asymptotic expansions in §2. Our approach is compared with that given in Jiji & Weinbaum (1978; for the corresponding cylindrical problem), which does not describe the manner in which the temperature decays in the liquid phase.

A simple exact solution to (1.1)–(1.8) for the zero latent heat limit *β*=0 is noted in §3, while the singular regime *κ*≪0 is considered briefly in §4. The latter limit corresponds to the heat diffusing very slowly in the liquid core, leading to an interior layer near the moving interface, and a classical one-phase problem in the solid phase. These ideas are discussed in Evans & King (2000), King & Evans (2005) and Struckmeier & Unterreiter (2001).

Care needs to be taken during the initial stages of the solidification process, since d*R*/d*t*→−∞ as *t*→0^{+}. In §5, we seek to understand the small-time behaviour by refining the analysis of Kucera & Hill (1986), which is not accurate for larger values of the Stefan number. All of the analytical results are compared with numerical solutions in §6, and finally a brief discussion is included in §7.

## 2. Large Stefan number limit *β*≫1

In this section, we generalize the large Stefan number analysis of Pedroso & Domoto (1973), Riley *et al*. (1974), Stewartson & Waechter (1976) and Soward (1980) to hold for the full two-phase problem (1.1)–(1.8). Our main goal is to determine the extent to which the inner liquid phase affects the outer solid phase and the evolution of the solid–melt interface. In doing so, we wish to establish under what circumstances the second phase (the liquid) can be ignored.

For *β*≫1, a large amount of latent heat is being produced at the solid–melt interface *r*=*R*(*t*), the result being that generally the interface moves very slowly. The two exceptions are for 1−*R*≪1 or *R*≪1; in these two regimes, the interface speeds up, with d*R*/d*t*→−∞ as *t*→0^{+} and , where *t*_{f} is the time it takes for the sphere to completely solidify.

In the large Stefan number limit, the two-phase problem breaks up into a number of different time scales, which are described below.

### (a) Time scale 1, *t*=*O*(1)

The first time scale is for *t*=*O*(1). On this short time scale (when compared with *t*_{f}, which turns out to be *O*(*β*)), the temperature in the liquid phase is approximately equal to *V* for most values of *r*, but rapidly decreases to *v*=0 near the free boundary *r*=*R*(*t*). Thus, there are two length scales to consider: one away from the free boundary and the other near the free boundary.

Note that the analysis on this time scale is not necessary in order to examine the next time scale, which is for *t*=*O*(*β*), since it happens that on the next time scale we are able to apply the initial conditions directly (without matching). However, the results for this time scale do shed light on the coupling between the two phases for small time, which is something we are most interested in.

#### (i) Inner region, 1−*r*=*O*(*β*^{−1/2})

For the inner region, we scale the spatial variables as and , and writeThe leading order partial differential equations are(2.1)with boundary conditions(2.2)(2.3)(2.4)The function in the last condition comes from matching onto the outer region, as subsequently described.

The next order problems are(2.5)with boundary conditions(2.6)(2.7)(2.8)Again, the function will be determined by matching with the outer region.

In terms of and the solutions to (2.1)–(2.4) and (2.5)–(2.8) arewith and satisfying the differential equations(2.9)and initial conditions , .

#### (ii) Outer region, 1−*r*=*O*(1)

The outer region is for 1−*r*=*O*(1). Here, we write , whereThe solution for is(2.10)

#### (iii) Matching between regions

By rewriting (2.10) in inner variables and expanding as *β*→∞, we findThus, matching between the two regions givesWe may now solve (2.9) for the moving boundary location, yielding

#### (iv) Small-time limit of time scale 1 solution

Although the solutions given in §2*a*(iii) involve infinite series, in practice for *t*=*O*(1) only a few terms are necessary to give an accurate solution. However, as *t* decreases, the terms begin to decay much more slowly, so that for *t*≪1 many terms are required to give an accurate solution.

This problem with convergence for *t*≪1 is overcome by noting that the function(2.11)has the property (Whittaker & Watson 1952, p. 273)(2.12)Thus(2.13)It follows that in the inner region we havewhere *E*_{1}(*z*) is the exponential integral defined byThe infinite sums in these expressions contain terms that decay very quickly for *t*≪1, and in particular, we now see thatas *t*→0^{+}, *β*→∞. By solving this equation asymptotically for *t*, we find(2.14)

#### (v) Summary of time scale *t*=*O*(1)

On the time scale *t*=*O*(1), we see that near the moving boundary the temperature in both the phases has an algebraic dependence on the small parameter *β*^{−1/2}. Furthermore, the two phases are coupled, although to a leading order both the temperature in the solid and the location of the free boundary are independent of the liquid phase.

### (b) Time scale 2, *t*=*O*(*β*)

The second time scale is for *t*=*O*(*β*); therefore, we rescale time as , where . As mentioned earlier, we can describe this time scale without reference to the previous time scale *t*=*O*(1), since we are able to apply the initial conditions (1.8) directly, and thus it is not necessary to match back onto *t*=*O*(1).

#### (i) Liquid phase

In the liquid phase, we have the heat conduction problemwith boundary conditionsand initial conditionswhere we may think of *R* as being a given function of .

This problem is similar to the one treated in Nayfeh (1973, p. 150). We set *ρ*=*r*/*R* and look for a solution of the formThus, to a leading order, we obtain the eigenvalue problemwhere the prime symbol denotes a derivative with respect to . The eigenvalues must therefore be of the form , where *n* is an integer, and therefore, after satisfying the initial condition *v*=*V* at , we find(2.15)

It is worth comparing our approach with that described by Jiji & Weinbaum (1978), who consider the corresponding problem in a cylindrical domain. Jiji & Weinbaum write down a perturbation series for the liquid phase in powers of *β*^{−1/2}, which by its very nature, can never capture the exponentially decaying behaviour given in (2.15). As a result, Jiji & Weinbaum find that the temperature at each order in the liquid phase is identically zero.

#### (ii) Solid phase and free boundary location

In the solid phase, we havewith boundary conditions(2.16)From the analysis in the liquid phase, we see that(2.17)thus, by noting the form of the Stefan condition (2.16), we expect the liquid phase to contribute exponentially small terms to the location of the free boundary.

With this in mind, and noting the large Stefan number analysis for the one-phase problem given in Pedroso & Domoto (1973) and Riley *et al*. (1974), we treat *r* and *R* as being the two independent variables, and write(2.18)The ellipses in (2.18) denote terms that are *O*(*β*^{−3}) and independent of *V*, while the term is exponentially small in *β* (and will depend on *V*).

We may read off the solutions(2.19)from Pedroso & Domoto (1973) and Riley *et al*. (1974) and, given (2.17) and the initial condition *R*=1 at *t*=0, we find from (2.16)_{2} that(2.20)But, given (2.18), we may rewrite the integral in the curly brackets as(2.21)The temperature in the inner core is now found by combining (2.15) and (2.21).

#### (iii) Small-time limit of time scale 2 solution

As with many of the expressions given in §2*a*(iii), parts of the solutions in §2*b*(ii) involve infinite series. For *t*=*O*(*β*), these series contain terms that decay extremely quickly, with the first term being sufficient for practical purposes. However, as we have *R*→1^{−}, andas *R*→1 and *β*→∞. Thus, care must be taken in this limit, as an increasing number of terms are required to obtain an accurate solution.

To obtain the limiting behaviour as *R*→1^{−} of our large Stefan number solution, we again use the property (2.12) of the function defined in (2.11). By expanding and then integrating term by term, we findas *R*→1^{−}. Thus, from (2.18) and (2.19) we can derive the expression (2.14).

#### (iv) Summary of time scale *t*=*O*(*β*)

On the time scale *t*=*O*(*β*), we see that the temperature in the solid phase has algebraic dependence on the small parameter *β*^{−1}, while the temperature in the liquid phase is exponentially small in *β*. Thus, the solid and liquid phases essentially decouple, with the liquid phase making only exponentially small contributions to the location of the solid–melt interface (the result is that for *t*=*O*(*β*), the solid phase and the free boundary location look almost identical to those found for the well-studied one-phase problem). In the one-phase limit *V*=0, our results reduce to those given previously in the literature (Pedroso & Domoto 1973; Riley *et al*. 1974), with the agreement being to the same order as that detailed in those studies.

### (c) Time scale 3, *t*−*t*_{f}=*O*(1)

Recall that *t*_{f} is the time to complete solidification and *t*_{f}=*O*(*β*). On the time scale *t*−*t*_{f}=*O*(1) (i.e. at times close to complete solidification), the interface no longer moves slowly and so the leading order behaviour in the solid away from the interface is no longer quasi-steady. We expect that at this stage the two phases completely decouple, and that the appropriate analysis coincides with the one-phase problem. The details are given in Riley *et al*. (1974), Stewartson & Waechter (1976) and Soward (1980).

### (d) Time scale 4,

As with the one-phase problem, further analysis is required to remove the non-uniformities that arise on the time scale *t*−*t*_{f}=*O*(1). This leads to an exponentially short final time scale, which is for . The details are presented by Stewartson & Waechter (1976), Soward (1980) and Herrero & Velázquez (1997), and are not repeated here.

### (e) Time to complete solidification

It is of interest to obtain an approximate value for *t*_{f}, the time to complete solidification. Unfortunately, the expansion (2.18) breaks down before *R*=0 and thus further analysis is required on the third time scale *t*−*t*_{f}=*O*(1), as mentioned above. However, we can read off the value of *t*_{f} for the one-phase problem (*V*=0) from Riley *et al*. (1974), Stewartson & Waechter (1976) and Soward (1980), and simply add the appropriate correction for the second phase (the liquid).

The leading order behaviour of the additional term is found by considering as *β*→∞. This approach works because the second phase contributes algebraic corrections to the solidification time only during the first time scale *t*=*O*(1) (the contributions from the later time scales are exponentially small, as mentioned above). The result is that(2.22)as *β*→∞, and thuswhere *ζ*(*z*) is the Riemann zeta function. The above result for *t*_{f} with *V*=0 is given by Riley *et al*. (1974), Stewartson & Waechter (1976) and Soward (1980). The derivation of (2.22) is detailed in appendix A.

### (f) Summary

In summary, for Stefan number *β*≫1 the extinction time is *O*(*β*). On the first time scale *t*=*O*(1), the leading order behaviour near the interface coincides with the one-phase problem, with the effects of the liquid phase only coming in at the first correction term. By the second time scale *t*=*O*(*β*), the liquid phase contributes only exponentially small terms to the location of the solid–melt interface, with the analysis in the solid phase coinciding with the one-phase problem to all algebraic orders. At times close to complete solidification, the temperature in the liquid essentially vanishes and the analysis follows the one-phase problem.

We note that this sort of qualitative behaviour is essentially implied by Jiji & Weinbaum (1978), who consider the equivalent problem for a cylinder. However, with their approach, Jiji & Weinbaum are neither able to describe the temperature field in the liquid for *t*=*O*(*β*) nor obtain an expression for the time to complete solidification.

## 3. Zero Stefan number solution *β*=0

When *κ*=1, there is an exact solution for the case *β*=0, i.e.(3.1)with the location of the moving boundary *r*=*R*(*t*) given implicitly byProvided *V*=*O*(1), the first term in the infinite sums above is usually sufficient, at least for practical purposes.

Using this exact solution, we find the complete solidification time *t*_{f} to be given implicitly byAs *V*→0, we have *t*_{f}→0 (very slowly), and in factalthough this leading order approximation is only accurate for extremely small values of *V*. Thus, the sphere freezes instantly in the zero Stefan number limit of the one-phase problem. This is to be expected, since with *β*=0 there is no latent heat liberated by the interface, and with *V*=0 no heat must diffuse out in order to lower the temperature in the core to the fusion temperature.

On the other hand, using a formula similar to (2.13) (see Whittaker & Watson (1952, p. 124), with *a*=1/2), we findwhich is of course consistent with the notion that an increase in the initial temperature should increase the freezing time.

## 4. Slow diffusion limit *κ*≪1

In the limit *κ*→0, which corresponds to slow diffusion in the liquid phase, the two-phase problem reduces to a one-phase problem as described by Evans & King (2000), King & Evans (2005) and Struckmeier & Unterreiter (2001).

The important point is that in order to derive a leading order model for *V*≠0, *κ*≪1, one cannot simply set *κ*=0 in the Stefan condition (1.7), since putting *κ*=0 into (1.2) suggests that *v*≡const., which is incorrect. The problem is thus singularly perturbed, with an interior layer developing near the solid–melt interface.

The interior layer is for , where . By writing in this region, we havewhich, provided that , implies that where we have matched with the region away from the sold–melt interface. It follows that for *κ*≪1, the approximate one-phase problem is:(4.1)(4.2)(4.3)(4.4)The main point here is that this free boundary problem is the same as the well-studied classical one-phase problem (relevant for *κ*=*O*(1) and *V*=0) except that the Stefan number *β* is replaced by *β*+*V*. It is now the size of *β*+*V* which determines the speed of the solid–melt interface.

## 5. Small-time perturbation *t*≪1

In this section we derive an approximate solution to (1.1)–(1.8), which is valid for small time. The approach used is similar to the one outlined in Kucera & Hill (1986); however, a number of key refinements have been made, which for certain parameter values lead to a much more accurate approximation.

### (a) Review of Kucera & Hill (1986)

In Kucera & Hill (1986), a small-time solution to (1.1)–(1.8) is sought by first making the domain fixing transformation , , and then applying the ansatz(5.1)as *τ*→0^{+}. With this change of variables, the liquid phase is transformed to 1≤*ξ*<∞, the solid phase is for 0≤*ξ*≤1 and the time domain is 0<*τ*<∞, with *τ*→0^{+} as *t*→0^{+}.

There are two main concerns with this approach. The first is that the ansatz (5.1) assumes that to leading order the solution in both the solid and the liquid phases is almost self-similar; however, this is only true for the outer solid phase. As a consequence of their approach, Kucera & Hill (1986) were unable to satisfy the initial condition (1.8)_{1}.

The second concern with using (5.1) is that, as a part of the analysis, a Taylor expansion is employed under the assumption that *ξτ*≪1 for *τ*≪1. However, the quantity is independent of *τ*, and in fact becomes very large close to the centre of the sphere, even for a very small time. This means that the scheme is unlikely to converge in the appropriate regime, regardless of how many extra terms in (5.1) are included.

In §6 of Kucera & Hill (1986) a second method is proposed, which involves small-time perturbation expansions of the formas *Y*→0^{+}, where(5.2)With this new transformation, the liquid phase is now for 1≤*X*<(1−*R*)^{−1} (note that it is not fixed), the solid phase is for 0≤*X*≤1 and the time domain is for 0<*Y*<1, with *Y*→0^{+} as *t*→0^{+}. In this case, Kucera & Hill (1986) were able to satisfy the initial condition (1.8)_{1} but not the no flux condition (1.5).

### (b) An alternative method

Since the small-time approximation derived by Kucera & Hill (1986) does not agree well with the numerical results for moderate to large values of the Stefan number *β* (figure 1), we endeavour to present an alternative approximation in this section.

#### (i) A summary of the details

As mentioned above, any small-time expansion which assumes that the leading order solution in the liquid phase is self-similar can never satisfy every one of the boundary conditions (1.3)–(1.7) and the initial conditions (1.8). Thus, this sort of approach can only be approximate in nature, as one boundary condition must be sacrificed, or satisfied only in the limit *t*→0^{+}.

The starting point for our alternative scheme is to look for solutions of the form(5.3)(5.4)as *Y*→0^{+}, where *X* and *Y* are defined in (5.2). We arrive at this ansatz by taking note of the solutions presented in §§2 and 3, as well as the small-time perturbation series developed by Davis & Hill (1982) for the one-phase problem. After substituting (5.3) and (5.4) in (1.1)–(1.4) and (1.6)–(1.8), we find that *A*_{0}, *A*_{1}, *B*_{0} and *B*_{1} satisfy a series of coupled ordinary differential equations with boundary conditionsThe leading order solutions are(5.5)where *γ* is the solution to the transcendental equation(5.6)with *L*_{s} and *L*_{ℓ} defined byEquation (5.6) also arises in the well-known Neumann solution (Carslaw & Jaeger, 1973, p. 285), which means that for small time the solidification process in the neighbourhood of *r*=1 behaves in the same way as the one-dimensional process in a semi-infinite domain.

The next order solutions are(5.7)whereNote that this approximate solution does not satisfy the boundary condition (1.5), and hence for each value of time *t*≪1, there will be a small region near *r*=0 for which the solution is not appropriate. We revisit this point in §5*c*.

An approximate location of the moving boundary is given by(5.8)which for small time gives(5.9)

We remark that for the one-phase problem with *V*=0, Hill & Kucera (1983) consider two more terms in (5.3), to include . Owing to the constraints in algebraic manipulations, such higher order terms have not been attempted here for the full two-phase problem.

#### (ii) Large Stefan number limit of small-time perturbation

In the limit *β*→∞, we have from (5.6) that *γ*→0, and in factSubstituting this expression into (5.9) givesThis statement extends the result (2.14).

Some lengthy calculations confirm that the large Stefan number limit of the present small-time perturbation for the solid phase agree with the small-time limit of the large Stefan number solutions given in §2. However, there is not the same level of agreement in the liquid phase, which is not surprising, as the solution in the liquid phase is not self-similar and the small-time perturbation is not uniformly valid there.

#### (iii) Zero Stefan number limit of small-time perturbation

Substituting *β*=0 and *κ*=1 into (5.6) givesThus, we find that, as *t*→0^{+}These solutions are in close agreement with the exact solutions (3.1), provided that *t*≪1.

#### (iv) Slow diffusion limit of small-time perturbation

By considering the limiting behaviour of the error function (Abramowitz & Stegun 1970, p. 298), we find thatwhich implies that for *κ*≪1, (5.6) is given approximately byThis well-known transcendental equation arises in the study of one-phase problems (Carslaw & Jaeger, 1973, p. 286), except that the usual Stefan number has been replaced by *β*+*V*. The conclusion is that for *κ*≪1 and *t*≪1, the temperature in the solid phase is essentially self-similar and governed by the usual one-phase equations with an adjusted Stefan number *β*+*V*. These one-phase equations are given by (4.1)–(4.4), and thus the slow diffusion limit of the small-time perturbation solutions are totally consistent with the analysis given in §4.

### (c) Further approximation

The temperature profiles generated by the small-time approximation of §5*b* do not satisfy the no flux condition (1.5), and in fact, this approach has the temperature in the liquid phase *v* blowing up in the limit *r*→0. In order to satisfy (1.5), we adjust (5.4) by writingwhere *B*_{0} and *B*_{1} are given by (5.5) and (5.7), respectively. Given that , this alteration can be thought of as adding an image solution in the domain −1≤*r*<0. The effect of this change is that the new solution now satisfies (1.1)–(1.5) and (1.8) exactly, but satisfies (1.6) and (1.7) only approximately; however, the errors from (1.6) and (1.7) are exponentially small in *Y* as *Y*→0^{+}. The result is that we are able to derive an excellent approximate solution to the two-phase Stefan problem, valid for *t*≪1.

## 6. Numerical results

The two-phase problem (1.1)–(1.8) is solved numerically using an enthalpy method, which we summarize here. The fundamental idea is to recast (1.1) and (1.2) in the form(6.1)where the enthalpy *h* is related to the temperature *T* by(6.2)and the fixed boundary and initial conditions become , and . A finite-difference scheme is applied to (6.1) to give(6.3)where the standard notation is used. After each time step the temperature is updated via (6.2). The location of the solid–melt interface can then be determined, again with the use of (6.2), as described in Voller & Cross (1981).

Figure 1 provides a number of temperature profiles for *β*=10, *κ*=1 and *V*=1, which is representative of a large Stefan number solution. Included in the figure are numerical results from the enthalpy method, the asymptotic results from §2*b*, the small-time approximation from Kucera & Hill (1986) and the small-time approximation described in §5*c*. In figure 1*a*, three profiles are drawn, for *R*=0.95, 0.875 and 0.8, respectively, to illustrate the small-time behaviour, while in figure 1*b*, profiles are drawn for later times. In figure 1*b*, only the outer solid phase is included since the temperature in the liquid phase for these times is essentially zero for all *r* (see (2.15) and (2.21)).

We observe in figure 1 that the asymptotic results derived in §2 agree extremely well with those from the numerical scheme, especially for smaller values of time. As the solid–melt interface approaches the centre of the sphere, we expect this approach to break down, as described in §2*c*, and this tendency can be observed from the figure. Further, we see that at least for these parameter values, the small-time approximation described in §5*c* provides a much closer agreement with the numerical solution than the one given by Kucera & Hill (1986), especially in the liquid phase. In the outer solid phase, there is even a reasonable agreement between our small-time solution and the numerical results for larger values of time.

In figure 2, the dependence of the solid–melt location on time is shown for the same parameter values as those in figure 1. The agreement between the numerical and analytical results is good for small time, while the large Stefan number solution (2.18) also works well for larger times. Note that only the first two terms in (2.18) are included in this figure.

To illustrate the behaviour for small Stefan numbers, temperature profiles are shown in figure 3 for the case *β*=0.1, *κ*=1 and *V*=1. In addition to the numerical solution, the exact zero Stefan number solution is drawn, as well as the two small-time approximations (from Kucera & Hill (1986) and §5*c*). Again, the agreement between each approach is good for small times, but it is worth noting that for small values of the Stefan number *β*, the small-time solutions do not work well for larger times (unlike the case with larger Stefan numbers), which is of course not unexpected.

Finally, to demonstrate the singular behaviour for slow conduction, temperature profiles are shown in figure 4 for *β*=1, *κ*=0.01 and *V*=1. In addition to the numerical solution and the two small-time approximations, results are shown from §4. That is, in the solid phase, the plot shown is the numerical solution to the one-phase problem with *β* replaced by *β*+*V*, while in the solid phase, the plot shown is for the interior layerwhere is obtained by solving the one-phase problem numerically. We see that (in the case of constant *V*) the leading order asymptotic solution in the interior layer provides an excellent approximation for the entire liquid phase (including the outer region). Note that the agreement between the small-time results from §5*c* and the numerical solution is very good for early times; however, for later times an artefact of the scheme presented in §5 is that it produces some slightly spurious behaviour in the liquid phase.

## 7. Discussion

It is somewhat surprising that the two-phase Stefan problem for spheres has received such little attention, given the significant amount of literature devoted to the idealized one-phase case. For large values of the Stefan number *β*, it can be loosely argued that since the solid–melt interface moves slowly relative to the speed at which heat conducts in the solid phase, the temperature in the inner liquid phase decays very quickly to the fusion temperature. Thus, except for small values of time, the liquid phase does not affect the solidification process and can be essentially ignored.

Using the method of matched asymptotic expansions, we have explored this structure in some detail. We have been able to confirm that on the time scale at which heat diffuses to the centre of the sphere, the leading order description of both the temperature in the solid phase and the location of the solid–melt interface is independent of the inner liquid phase, although there is a coupling between the phases at higher orders. However, on the much longer time scale at which solidification occurs, the solid and liquid phases essentially decouple, with the temperature in the liquid phase being exponentially small. This work extends the existing analysis for the one-phase problem.

While the present study is limited to the simple spherical geometry, the scalings for the large Stefan number limit would be expected to carry through to the more general problem of the inward solidification of a truly three-dimensional region of liquid. Thus, we might expect the near-complete solidification results of McCue *et al*. (2005) for the one-phase problem with arbitrary geometry to also hold in the two-phase case.

Consideration has also been given to deriving an approximate solution valid for small time. In the one-phase problem, this can be done by employing (5.2) and then writing out the temperature *u* in the form (5.3), as in Davis & Hill (1982). This small-time solution has been extended to include four terms *A*_{i}, *i*=1, 2, 3, 4 by Hill & Kucera (1983). Adapting this approach to the two-phase problem is not straightforward since the leading order behaviour for the two-phase problem is not self-similar. However, with the use of an image solution, we have been able to derive an approximate solution that works very well for small values of time.

We conclude by noting that an analogous study to the one presented in this paper could be undertaken for cylindrical geometry. In particular, the large Stefan number analysis of Jiji & Weinbaum (1978) is incomplete, and for example does not describe the manner in which the temperature of the liquid phase decays. We would expect many of the scalings identified in §2 to also hold for the cylindrical case, the most notable exceptions being for the later time scales mentioned in §2*c*,*d* (e.g. Soward 1980). For a more general two-dimensional geometry, the corresponding large Stefan number behaviour at times just before complete solidification is described in McCue *et al*. (2003).

## Acknowledgments

The authors are very grateful to Pei Tillman for reading through an earlier version of the manuscript, and for many fruitful discussions. B.W. and J.M.H. acknowledge support from the Australian Research Council. J.M.H. is grateful for the provision of an Australian Professorial Fellowship.

## Footnotes

- Received November 10, 2007.
- Accepted March 13, 2008.

- © 2008 The Royal Society