MOVING BOUNDARY PROBLEMS FOR QUASI-STEADY 1 CONDUCTION LIMITED MELTING

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


Introduction.
While there is a variety of simple models to approximate the 21 shape of a melting particle [33,38], the traditional approach from a mathematical 22 perspective is to employ a Stefan problem, which involves the linear heat equation  44,46,56]. We continue this direction in the present study, focusing on the 29 melting of an axially symmetric dendritic crystal. We employ both analytical and 30 numerical techniques to study the shape of the evolving crystal, focussing on the very 31 final stages of melting.

32
A key aspect of a traditional Stefan problem is that the effects of convection zone. Finally, the temperature is raised to remelt the crystals, returning the system 44 to a stable melt phase.
t → t − e , meaning that the crystals were spherical just before extinction. 56 In order to make analytical progress, Glicksman et al. [21] model the process 57 with a one-phase quasi-steady problem, which results by ignoring heat conduction 58 within the crystals and assuming an infinite Stefan number. Here, the Stefan number 59 is defined by where c is the specific heat, L the latent heat of fusion per mass and u * ∞ − u * m is 62 the temperature difference between the melt away from the crystal and the melting 63 temperature. In reality, for this particular experiment the parameter values were 64 L/c ≈ 10.99 K, u * ∞ −u * m ≈ 1.8 K, so β ≈ 6.1, which is not reasonably large. Glicksman Glicksman and co-workers [21, 22, 43] did not provide an explanation for the ob-74 served increase in aspect ratio during the first 50 s of melting; however, the subsequent 75 decrease in aspect ratio (during the final 10 s of melting) was accounted for by noting 76 that by this stage of the melting process the crystals had become small enough for 77 surface tension effects to begin to dominate [22,43]. As a consequence, the needle 78 tips with high curvature melted more quickly than the remainder of the crystals, in In this article, we are motivated by these issues to undertake a theoretical study 88 of the one-phase quasi-steady Stefan problem. The mathematical problem is re-89 formulated in Section 2 with a Baiocchi transform for the special zero-surface-tension 90 case. In Section 3, we go on to provide a near extinction analysis for a general shaped 91 initial crystal, including numerical results for cases in which crystals ultimately melt 92 to a single point or pinch off and break into two separate pieces. initially at melting temperature u * m occupying the region Ω * (0), surrounded by the 109 same substance in liquid form in R 3 \ Ω * . In the far field, a higher temperature u * ∞ is 110 applied, and thus melting proceeds until the crystal melts completely at the extinction 111 time t * e .

112
Setting k to be the thermal diffusivity, we scale variables using where is a characteristic length scale of the initial crystal shape, and β is the Stefan

115
This manuscript is for review purposes only. For what follows we shall take the quasi-steady limit β = ∞, which is an appro-124 priate approximation for experiments in which the latent heat is large or the specific 125 heat is small. As a result, the parabolic equation (2.2a) becomes Laplace's equation and thus we do not require an initial condition for u.

129
As mentioned in the Introduction, the governing equations (2.2e) with (2.2b)-

130
(2.2d) are also relevant for the problem of a bubble that is forced to contract in a 131 saturated medium, where the fluid flow is governed by Darcy's law [12,28,45], as 132 well as the two-dimensional analogue for Hele-Shaw flow [15,14,42]. These equations 133 also arise in other moving boundary problems, for example the small Péclet num-134 ber limit of advection-diffusion-limited dissolution/melting models [6, 27, 32, 53, 57], 135 for which it is also of interest to track the moving boundary and predict its shape where we are using the notation t = ω(x) to denote the solid-melt interface ∂Ω. The

148
Transforming the governing equations (2.2e) with (2.2b)-(2.2d), we derive the 149 nonlinear moving boundary problem for w: in Ω(0) \ Ω(t) : Once a solution for the Baiocchi variable w is determined, the temperature u can be 157 recovered via u = ∂w/∂t. We note that an advantage of the Baiocchi transform is 158 that it transforms the inhomogeneous boundary condition (2.2c) into a homogeneous 159 boundary condition. Another is that time appears as a parameter in (2.4a)-(2.4e), 160 so that the problem can be solved at any time without knowledge of the solution at 161 previous times.  We present here a summary of this exact solution in the special case for which 169 the initial crystal shape ∂Ω(0) is the prolate spheroid 170 (2.5) x 2 + y 2 + z 2 z 0 (0) 2 = 1, The exact solution is that ∂Ω(t) retains its prolate spheroidal shape as where z 0 (t) > 0 and ρ 0 (t) > 0 measure the major and minor axes of the dendrite, 176 respectively, with constant aspect ratio A(t) = z 0 (t)/ρ 0 (t) = z 0 (0) (here the length 177 scale is chosen so that ρ 0 (0) = 1). The full solution has the time-dependence which provides an interesting connection between our problem and gravity potential 209 generated by a uniform body.

210
Whilst in practice it is not feasible to compute W analytically for a general initial  and x e is a local minimum of w e , a simple Taylor series for this axially symmetric 224 geometry implies that w e ∼ a(x 2 + y 2 ) + bz 2 as r → 0. Further, as a consequence of 225 (3.1b), we then have where 1/6 < a < 1/4. As we shall see, the parameter a is effectively all the melting 228 crystal "remembers" from its initial condition; it is this single parameter that controls 229 the aspect ratio of the crystal at extinction. Note that the higher order terms in (3.3) 230 are not required in the following analysis (they would be for the special case a = 1/4, 231 which represents the borderline between the type of extinction considered in this 232 section and when a bubble breaks up into two, as treated in Subsection 3.4).

233
In the limit t → t − e , the inner region is for r = O(T ), where T (t) is a length scale 234 defined so that the volume of the melting crystal is fixed to be 4πT 3 /3. We write where Ω 0 denotes the crystal which has volume 4π/3 in these self-similar coordinates, 240 and N denotes a normal direction. In order to match with (3.3) we require that as R → ∞, where d is a constant found as part of the solution to (3.4a)-(3.4c). We 243 see from (3.4c) that a matching condition for the outer region is The solution to (3.4a)-(3.4c) in prolate spheroidal coordinates is provided in Ap-246 pendix A. According to this solution the dendrite boundary ∂Ω 0 is described by where q 0 is a parameter that is related to the special constant a by Note that the prolate spheroid approaches a perfect sphere in the limit a → 1/6 + , in The outer region is for r = O(1), for which Matching with the inner gives the time-dependence Thus we see that, regardless of the shape of the initial crystal, the square root of time 261 scaling determined experimentally in Glicksman et al.

262
In summary, the zero-surface-tension model predicts that, provided there is no 263 pinch-off, an axially symmetric dendrite will melt to a spheroid in the extinction limit.

264
While this spheroid could be prolate or oblate, we concentrate here on the prolate 265 case, as this is the one observed in the IDGE [21, 22, 43]. The aspect ratio of the 266 prolate spheroid at extinction is given by which provides an implicit dependence of A on the constant a via (3.7). Here a is 269 the only parameter that is required to characterise the initial dendrite shape (it is  .7)).

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

278
Here Φ = R 2 /2 − 1/2 + 1/3R and the dendrite becomes spherical in the limit with (3.14) Thus, the left-hand side vanishes as T → 0, or t → t − e , meaning that the exterior of 309 the crystal approaches a null quadrature domain in the limit, and thus the crystal 310 itself approaches an ellipsoid in shape. late a level set function, φ(x), such that φ > 0 for x ∈ Ω(0) and φ < 0 for x ∈ R 3 \Ω(0).

314
Thus we can reformulate (3.1a) and (3.1b) as where H is the Heaviside function. We note that H(φ) is discontinuous at x ∈ ∂Ω(0), 317 so for numerical purposes we implement a smoothed Heaviside function where δ = 1.5∆x. For this purpose, it is convenient to work in spherical polar 320 coordinates (r, θ, ϕ) and represent the axially symmetric moving boundary ∂Ω by 321 r = s(θ, t). Thus, (3.15) becomes The spatial derivatives in (3.17) are approximated using central finite differencing, 325 with homogeneous Neumann boundary conditions applied at r = 0, θ = 0, and θ = π.

326
The far-field boundary condition (3.1c) is incorporated using a Dirichlet-to-Neumann 327 map described in Appendix B.2.2. 328 3.4.1. Symmetric initial condition. We consider a selection of initial condi-329 tions to illustrate a few different qualitative behaviours. Again, using spherical polar 330 coordinates (r, θ, ϕ) with ∂Ω denoted by r = s(θ, t), the first is the prolate spheroid where r 0 describes the initial aspect ratio. The second initial condition is a peanut- for 0 ≤ θ ≤ π/2; for π < θ ≤ 2π this initial condition is made symmetric by reflecting 346 about θ = π/2.

347
In Figure 2, we illustrate some numerical results by choosing parameter values 348 from these three initial conditions. For the prolate spheroid (3.18) we provide results 349 for r 0 = 0.8, noting that this initial condition is obviously convex. For the peanut 350 shaped surface (3.19), we choose r 0 = 0.5, which is not convex but is instead mean  For initial condition (3.20a)-(3.20d) with r 0 = 4.75, Figure 2 shows different qual-375 itative behaviour. Here, we see that solutions to (2.2b)-(2.2e) will undergo pinch-off 376 and ultimately the two satellite crystals will contract to separate points of extinction. Baiocchi transform approach can be used to predict whether pinch-off will occur, as 385 well as the extinction points and time.

386
In summary, these numerical results indicate that for a given initial interface,  interfaces which form after pinch-off will have the same extinction time. We now 400 investigate a class of asymmetric initial conditions that undergo pinch-off into two 401 surfaces of differing volumes. We expect the smaller of the two volumes to contract 402 to a point first, followed by the larger, thus giving two distinct extinction times. 403 We again consider an initial condition of the form of (3.20a)-(3.20d), but this 404 time for 0 ≤ θ ≤ π. In Figure 4, we plot the time evolution of the numerical so-405 Fig. 3: The evolution of the aspect ratio for the example initial condition (3.19) with r 0 = 0.5 is presented as a solid (blue) curve. The (red) dashed curve is the predicted aspect ratio at extinction, given by (3.12). The agreement is quite good.   We seek a perturbed spherical solution to (4.2a)-(4.2d) of the form 436 u(r, θ, ϕ, t) = u 0 (r, t) + εu 1 (r, θ, t) + O(ε 2 ), (4.3a) A n r −n P n (cos θ), where A n is a sequence of unknown coefficients, P n is the nth Legendre polynomial, 450 and γ n is the nth mode of perturbation to the sphere. We are able to eliminate A n 451 to obtain 452 (4.7) 1 γ n dγ n ds 0 = (n − 1)((n 2 + 3n + 4)σ + s 0 ) s 0 (s 0 + 2σ) .

453
Since (1/γ n )dγ n /ds 0 → 0 in the limit that s 0 → 0 for n ≥ 2, we see that each mode of 454 perturbation is stable, and a perturbed sphere will evolve to a sphere in the extinction 455 limit, as expected.

467
Note that when σ = 0, then 3γ 2 /2s 0 = 1/r 0 , resulting in the aspect ratio remaining   The inner region is for r = O(ρ 0 (t)). Here the melting is almost two-dimensional 477 with ∂u/∂z 1 and ∂S/∂z 1 so that, to leading order, 478 in ρ > S(z, t) : The solution to (4.11a)-(4.11c) is where the form for S is determined by the missing far-field condition, which is found 485 by considering the outer region.

486
In this outer region, which is for r = O(z 0 (t)), the dendrite appears as a slit. We 487 scaleρ = ρ/(αρ 0 (t)),t/ ln α and rewrite the inner solution (4.12) to be 488 (4.13) The leading order solution in the outer region is u = 1, thus, after matching with the 490 leading order term in (4.13) as α → ∞, we find 491 (4.14) For the zero surface tension case σ = 0, we can solve (4.14) explicitly to give again providing square root time dependence.

495
Of particular interest is the special case in which the initial dendrite is the prolate 496 spheroid (2.5). Here ρ 0 = α and z 0 (0) = 1, so initially the dendrite has the aspect 497 ratio A(0) = 1/α. From (4.14) we find the interface is given implicitly by Note that the small parameter in this limit is 1/ ln α, which suggests the analysis here 500 is valid only for extremely large aspect ratios. spheroidal crystal considered in Subsection 2.2, whose surface is (2.6), we find the 503 mean curvature is largest near the tip, given by the initial crystal,Ω(0), is a prolate spheroidal in shape, this is a canonical problem 516 for melting a solid. This one parameter in the problem is the initial aspect ratio.

517
Using the numerical scheme described in Appendix B, we solve (4.19a)-(4.19d) 518 forû andΩ. We first consider a near spherical prolate spheroid initial condition 519 such that the initial aspect ratio is close to unity. Figure 5 compares the aspect ra-  suggesting that an initially prolate spheroidal crystal will tend to a sphere in the 547 extinction limit. This conclusion is that same as before in Subsection 4.1 when c = 0.

548
On the other hand, a significant difference in qualitative behaviour is that the aspect 549 ratio with c > 0 may first increase and then decrease (to unity), which is a feature Given s 0 is defined on the domain 0 ≤ s 0 ≤ r 0 , the aspect ratio will monotonically 553 decrease to unity if otherwise, the aspect ratio will be non-monotone.

556
Our work is motivated in part by a series of experiments performed as part of 557 the IDGE [21, 22, 43]. In these experiments, it was observed that the aspect ratio of 558 melting crystals increased for a period of time before decreasing to unity at extinction.

559
In the context of the results presented in this section, Figure 7 illustrates the aspect  6. Discussion. In this paper, we have studied a quasi-steady one phase Stefan 566 problem for melting an axially symmetric crystal. In Section 3 we treat a zero-surface-567 tension model and use analytical tools to show that axially symmetric crystals will 568 tend to prolate spheroids in the limit that they melt completely, namely t → t − e , with 569 an aspect ratio that depends on the initial condition. The point to which the crystals 570 ultimately shrink, together with the melting time, is predicted by this analysis and 571 confirmed using a novel numerical scheme based on the level set method (presented 572 in Appendix B). An advantage of this scheme is that we are also able to present 573 numerical results for crystals that undergo pinch-off and contract to multiple points 574 of extinction.

589
A key assumption in our paper is that the Stefan number in ( contraction in a Hele-Shaw cell [15,14,42]). For the case in which a bubble pinches 598 off to produce two shrinking bubbles, the problem formulation would also need to 599 consider two points of extraction that coincide with the eventual extinction points.

600
The second issue is that, strictly speaking, for the extremely late stages of melting, where ξ ≥ 0, 0 ≤ η ≤ π, 0 ≤ φ < 2π, and k is a constant to be determined below.

616
The crystal boundary ∂Ω 0 is described by ξ = ξ 0 or, equivalently, Motivated by the relationship   cally is outlined as follows:

667
Step 1 For a given initial condition s(θ, 0), construct a level set function φ(r, θ, 0) 668 such that φ < 0 inside the interface and φ > 0 outside the interface. This 669 function is then converted to a signed distance function using the method of 670 crossing times as described by Osher & Fedkiw [51].

671
Step 2 Compute the temperature, u, on the domain r ≥ s(θ, t) using the procedure 672 described in Appendix B.2.

673
Step 3  Step 5 Reinitialise φ every 5 time-steps to a signed distance function by solving the

686
Step 6 Repeat steps 2-5 until the desired simulation time is attained.
For the singularity at θ = 0, noting that ∂u/∂θ = 0 and using L'Hôpital's rule then  Fig. 8: Schematic of how the speed function, F , is computed for each time step. Blue region denotes where temperature, u, is solved for using finite differences. This finite difference stencil must be adjusted to incorporate the dynamic boundary condition (Appendix B.2.1). To incorporate the far-field boundary condition, we impose an artificial boundary at r = R and implement a Dirichlet to Neumann mapping (Appendix B.2.2). F is computed outside the interface using (B.3), and is extended to be defined over the entire computational domain by solving the biharmonic equation.
The same procedure is applied at θ = π. Difficulties arise when attempting to in-697 corporate the dynamic condition (2.2b) on the interface and the far-field boundary 698 condition (2.2d). We detail the methodology used to overcome each of these difficul- where changes in topology occur [19]. However, an advantage of using a finite differ-706 ence stencil is that it can easily be adapted to problems where the boundary integral 707 method is not applicable. For example, we have used a similar method to the one 708 presented in this section to study non-standard Hele-Shaw flow where pressure is not 709 harmonic and for which the boundary integral method is much less suitable [50]. twice. Supposing the interface is located between two nodes (i − 1, j) and (i, j), the 716 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 717 r b denotes the location of the interface and u b is the temperature at the interface.

718
The value of r b is found by noting that φ is a signed distance function and so the 719 distance between r b and r i,j , denoted h, can be calculated by 720 (B.9) h = ∆r φ i,j φ i,j − φ i−1,j . The same procedure is applied if the interface is between r i and r i+1 , or in the 724 azimuthal direction. (n + 1)P n (cos θ j )P n (cos θ k ) sin θ k .

758
From a practical perspective, we cannot, of course, evaluate the series in (B.16) using 759 an infinite number of terms, but have found that using 10 terms gives sufficient accu-