Thermal Effect of Surface Tension on the Inward Solidification of Spheres

(1)

Thermal eﬀect of surface tension on the inward

solidiﬁcation of spheres

T. Wu

*

_{, H.-C. Liaw, Y.-Z. Chen}

Department of Mechanical Engineering, National Taiwan University, Taipei 106, Taiwan Received 13 April 2001; received in revised form 21 September 2001

Abstract

Effect of surface tension (Gibbs–Thomson effect) on the inward solidification of a liquid in a spherical container is investigated analytically by solving the unsteady heat equation via a small-time series expansion technique. A nonlinear Shanks transformation is adopted to improve the convergence property of the series solution at large time. The results show that at fixed Stefan number, the effect of surface tension is to increase the growth rate of the freezing front. A local minimum in the freezing rate is found to develop for all surface tension parameter values considered in this study. Also, analytic expressions for the relations between the growth rate of the freezing front, Stefan number and surface tension parameter are derived under the asymptotic condition of small Stefan number. Ó 2002 Elsevier Science Ltd. All rights reserved.

1. Introduction

Solidification is a common process in many practical engineering applications such as refrigeration, ice for-mation, casting of alloys, and crystallization of pure substances. Generally speaking, such a process is fea-tured by the existence of a moving interface between different phases at which the thermal energy in the form of latent heat is liberated. Owing to the unknown loca-tion of this solid/liquid interface and the nonlinear form of the thermal energy balance condition at the interface, analytical solutions are difficult to obtain except for few simple configurations. One of the few earliest exact so-lutions is that by Neumann and separately, by Stefan (see [1]) who studied the cooling of fluid from one end of a one-dimensional semi-infinite region. They obtained a closed solution for the temperature distribution in a self-similar form and deduced that the location of the in-terface advances with the square root of time. The specific solution has since been the starting point for the subsequent numerous approximate solutions to other types of boundary conditions and geometries.

Hence-forth, the problems involving phase changes and moving interfaces were traditionally categorized as ‘Stefan problems’.

References relevant to the one-dimensional Stefan problem were well documented in the book by Hill [2]. Specifically, various techniques and approximations have been attempted to study the problem of an inward solidification with spherical symmetry. These include the power-series approximation method of Kreith and Ro-mie [3]; the closed-form solution by Langford [4] under the assumption of constant solidification rate; the ap-proximate integral method of Poots [5]; the numerical solution by Tao [6], Li [7], and Caldwell and Chan [8]; the perturbation method by Huang and Shih [9]; the small parameter expansion solution by Pedroso and Domoto [10]; the strained coordinates method by Pedroso and Domoto [11]; the matched asymptotic ex-pansions solution by Riley et al. [12]; the asymptotic solution (near the final stage of complete solidification) of Stewartson and Waechter [13]; the small-time ex-pansion solutions by Davis and Hill [14], and Hill and Kucera [15].

However, all the above results did not take into account the thermal eﬀect of the surface tension at the interface. Surface tension aﬀects the equilibrium temp-erature (freezing or fusion temptemp-erature) at the interface of the two phases via the so-called Gibbs–Thomson law. If

www.elsevier.com/locate/ijhmt

*_{Corresponding author. Tel.: +886-2-236-302-31x2413; fax:}

+886-2-236-317-55.

E-mail address: tywu@ccms.ntu.edu.tw (T. Wu).

(2)

we let c denote the surface tension coeﬃcient, L the latent heat per unit volume of the solid phase, R_{the radius of}

the spherical interface, and T

f the freezing/fusion

tem-perature when the interface is ﬂat; then the equilibrium temperature on the spherical freezing front is given by the following Gibbs–Thomson condition ([16]):

T¼ T f 1 þ 2c LR : ð1Þ

That is, the equilibrium temperature on a curved inter-face (with center of radius lies in the liquid phase) is higher than that on a planar surface, and is increasing with decreasing radius of curvature of the interface. It should be noted that in deriving Eq. (1), a small value in ðT_T

fÞ=Tf has been assumed. Thus expression (1) is

not valid when the second term 2c=LR _{in the}

paren-thesis is of Oð1Þ, i.e., not suitable for freezing front with very small radius of curvature. In the present study, we shall investigate the inward freezing of liquid in a spherical domain (such as liquid in a spherical container, or a liquid droplet) under the modiﬁed equilibrium condition (1) as imposed by the surface tension at the interface.

As in all the previous works mentioned above, for-mulation of the present study is based on the classical one-phase model ([2]) that has been traditionally applied to most of the Stefan problems. In this one-sided

ap-proach, the latent heat liberated from the freezing front is assumed to be carried away by the solidified phase only. This simplification leads to the familiar unsteady heat-diffusion equation for the solid phase with moving interfacial boundary. The equation describing the pro-gression of the interface is derived from equating the amount of heat released per increase of the freezing layer to that diffused into the solidified field (Stefan condi-tion). Since the Gibbs–Thomson condition is not valid when the freezing front comes near to the point of complete solidification, the set of equations is solved by means of a small-time series expansion technique. A nonlinear Shanks transformation ([17]) is applied to the series in order to improve the convergence of the solu-tion at large time. The effect of surface tension on the temperature distribution in the solid field is discussed. Also, the asymptotic behavior of the freezing rate in terms of relevant physical parameters is determined to the leading-order approximation of small cooling pa-rameter (Stefan number) value.

2. Formulation of the problem

The governing equation is that of the one-dimen-sional unsteady heat conduction equation written in the spherical coordinate r

Nomenclature

An; Bn integration constants in nth-order solution,

deﬁned in Eqs. (23) and (30)

c heat capacity (per unit volume) of solid phase eðkÞ

m sequence resulted from the kth-order Shanks

transformation

L latent heat (per unit volume) of freezing M conﬂuent hypergeometric function deﬁned in

Eq. (28)

p; q orders of the conﬂuent hypergeometric function

R nondimensional radial position of the freezing front

Rs nondimensional radius of the sphere

R0

m sequence of temperature gradient deﬁned in

Eq. (36)

r nondimensional spherical coordinate

Sm m-term partial sum of the series solution (18)

S0

m m-term partial sum of the temperature

gradient kn

T nondimensional temperature

Tm sequence of temperature solution deﬁned in

Eq. (32)

t nondimensional time

Greek symbols

a thermal diﬀusivity of solid phase

b nondimensional surface tension parameter D Stefan number

e small quantity in the Gibbs–Thomson law, deﬁned in Eq. (A.1)

c surface tension at the solid/liquid interface

g transformed radial coordinate deﬁned in Eq. (11)

kn temperature gradient at interface evaluated

from the nth-order solution

H transformed temperature variable deﬁned in Eq. (13)

Hn transformed nth-order temperature solution

r parameter combining the eﬀects of surface tension and Stefan number

s time-like coordinate deﬁned in Eq. (11) sc critical time at which the freezing rate is

minimum

n transformed variable used in the conﬂuent hypergeometric function (28)

Superscripts

(3)

oT ot ¼ a 2 r oT or þo 2_T or2 ; R6r6R_s; ð2Þ where T _{represents the temperature, a is the thermal}

diﬀusion coeﬃcient of the solid, R

s denotes the radius of

the sphere, and R_ðt_{Þ the instantaneous location of the}

solid/liquid interface (see Fig. 1). Suppose that the temperature at the outer surface of the sphere is sud-denly dropped to T

s below the freezing temperature at

t_{¼ 0}þ _{and kept constant afterwards, then the proper}

boundary condition is T_{¼ T}

s at r¼ Rs: ð3Þ

The equilibrium temperature at the interface is that due to the Gibbs–Thomson condition:

T_{¼ T} f 1 þc L 2 R_ðt_Þ at r_{¼ R}_ðt_Þ: _ð4Þ

Balancing the heat releasing rate and the heat conduc-tion rate at the freezing front gives the equaconduc-tion for the motion of the interface:

LdR dt ¼ ac oT or R ; ð5Þ

where c denotes the heat capacity (per unit volume) of the solid phase.

Deﬁning relevant dimensionless variables as

T¼T _T f T s Tf ; r¼ r R s ; R¼R R s ; t¼at R2 s ; ð6Þ the nondimensional form of the Eq. (2) becomes oT ot ¼ 2 r oT orþ o2_T or2; R 6 r 61: ð7Þ

The associated boundary conditions are

T¼ 1 at r¼ 1; ð8Þ T ¼ b D 1 RðtÞ¼ r RðtÞ at r¼ RðtÞ; ð9Þ and the nondimensional growth rate of the freezing front is given by dR dt ¼ D oT or R : ð10Þ

In the above expressions, the dimensionless parame-ter D¼ cðT

f TsÞ=L is called the Stefan number, which

measures the degree of cooling. The parameter b¼ 2ccT

f=ðRsL

2_{Þ represents the surface tension eﬀect}

(Gibbs–Thomson effect), and r¼ b=D is merely a con-venient dimensionless parameter, which can also be taken as a measure of the surface tension when the Stefan number D is held fixed. Surface tension force is significant only when the radius of curvature is small. Water has a small surface tension coefficient (c¼ 0:077 N/m). Considering a liquid droplet of size 103_{m, the}

parameter r for water–ice system is then about 3 104_.

Melted metals have large surface tension (e.g., c¼ 1:9 N=m for iron, 2.0 N/m for nickel), typical values of r for common metals melts are around 102_{. In view}

of Eq. (9), surface tension has therefore a nonnegligible eﬀect in solidiﬁcation, particularly, when the radius of the freezing front RðtÞ is small.

3. Solution to the problem

The problem involves a moving boundary (interface) whose position is part of the solution and hence is not known a priori. A boundary-ﬁxed transformation [2]

g¼ 1 r 1 RðtÞ¼

1 r

sðtÞ ; ð11Þ

is adopted ﬁrst to map the unknown position of the freezing front to a ﬁxed value in g, g¼ 1. In Eq. (11), s¼ 1 R is another variable describing the instanta-neous position of the freezing front, which can also be treated as a time-like coordinateð0 6 s 6 1Þ. By using the chain rule, the original Eq. (7) can be rewritten in terms of the new coordinatesðg; sÞ

dR dt oT osþ g s dR dt oT og¼ 2 ð1 gsÞs oT ogþ 1 s2 o2_T og2: ð12Þ

A further transformation of temperature T will put the equation in a form amenable to the series-expansion technique,

T ¼H r ¼

H

1 gs: ð13Þ

In terms of the new temperature variable H, Eq. (12) becomes

(4)

s2dR dt oH osþ gs dR dt oH og ¼ o2_H og2 : ð14Þ

The corresponding boundary and interface conditions are

H¼ 0 at outer surface g¼ 0; ð15Þ H¼ r at interface g¼ 1; ð16Þ and the growth rate of the interface is governed by

dR dt ¼ D s 1 1 s oH og g¼1 þ s ð1 sÞ2H g¼1 ! : ð17Þ The forms of Eqs. (14) and (17) suggest that we can expand the solution into a power series of the time-like variable s, i.e.,

Hðg; sÞ ¼X

n¼0

sn_H

nðgÞ: ð18Þ

Substituting (18) into Eqs. (14) and (17) and collecting the like-power terms for each diﬀerent order of s, we obtain the following system of equations:

H00_nþ Dk0½gH0nþ nHn ¼ H00 n1þ D Xn1 k¼0 gðkkþ1 ( rÞH0 nk1 X n2 k¼0 ðn k 1Þðkkþ1 rÞHnk1 ) ; n2 I; nP0; ð19Þ where Pm_k¼0ðÞ ¼ 0 is deﬁned whenever m < k, and kn¼ ðdHn=dgÞj_g¼1 denotes the temperature gradient

evaluated at the interface at each diﬀerent order n. In Eq. (19), the prime symbols are shorthand notations for diﬀerentiations with respect to the coordinate g. The corresponding boundary conditions at each order n are

H0ð0Þ ¼ 1; Hnð0Þ ¼ 0 for n P 1; ð20Þ

H0ð1Þ ¼ r; Hnð1Þ ¼ 0 for n P 1: ð21Þ

The solution to the zeroth-order equation is easily obtained as H0¼ A0 Z g 0 eðDk0=2Þg2_dg_{þ B} 0; ð22Þ where A0¼ r 1 R1 0eðDk0=2Þg 2 dg; B0¼ 1; ð23Þ are the two integration constants determined from the boundary conditions (20) and (21). The temperature gradient at the interface is then obtained through dif-ferentiating Eq. (22): k0 dH0 dg _g¼1¼ ð1 þ rÞ R1 0 eðDk0=2Þg 2 dge ðDk0=2Þ_: _ð24Þ

Notice that the temperature gradient k0appears at both

sides of Eq. (24). An iterative procedure is then required to calculate the numerical value of k0.

The homogeneous parts of all subsequent higher-order equations (n P 1) can be recast into a normal form via the coordinate transformation n¼ Dk0g2=2

nH00_nþ 1 2 n H0_nþn 2Hn¼ 0: ð25Þ Equation (25) falls into the class of the conﬂuent hy-pergeometric equation ([18]), whose general form is

nH00þ ðq nÞH0_{pH ¼ 0:} _ð26Þ

The two independent solutions to the above equation are

Mðp; q; nÞ and n1qMðp q þ 1; 2 q; nÞ ð27Þ in which Mðp; q; nÞ is the conﬂuent hypergeometric function of orderðp; qÞ deﬁned as

Mðp; q; nÞ ¼X 1 k¼0 ðpÞ_k ðqÞk nk k! ð28Þ with ðpÞk¼ p ðp þ 1Þ ðp þ 2Þ ðp þ k 1Þ for k P 1; ðqÞ_k¼ q ðq þ 1Þ ðq þ 2Þ ðq þ k 1Þ for k P 1; ðpÞ₀¼ 1; ðqÞ₀¼ 1:

Note that the power series in (28) terminates at ﬁnite terms for integer value of p. The solution to our present Eq. (25) corresponds to the case ðp; qÞ ¼ ðn=2; 1=2Þ. For certain special values ofðp; qÞ, the hypergeometric function (28) can be related to the more familiar expo-nential and error functions [18]; e.g.

Mðp; p; nÞ ¼ en _{for any number p;}

M 1 2; 3 2;n ¼ 1ffiffiffiffiffiffiffi n p Z pffiffiffiffin 0 ez2 dz¼ ffiffiffi p p 2pffiffiffiffiffiffiffinerfð ffiffiffiffiffiffiffi n p Þ: Thus by using proper recursive formulae ([18]) between different orders of the hypergeometric function Mðp; q; nÞ, solution Hnto the homogeneous Eq. (25) can

be rewritten into polynomial and error functions in g. Having obtained the solution to the homogeneous part of Eq. (19) for all n, the particular solution at each diﬀerent order n can be found by using the method of variation of parameters. The complete solution at each order n is then the sum of the two.

Since the expression of the solution becomes rather lengthy and cumbersome very quickly with increasing

(5)

order n, only the solution corresponding to the ﬁrst or-der in s (i.e. OðsÞ) is given here:

H1¼ A1g B1 eðDk0=2Þg 2 Dk0g Z g 0 eðDk0=2Þg2_dg þ A0 3k0 k0 ð þ k1 rÞgeðDk0=2Þg 2 ð29Þ with integration constants

A1¼

A0

3k0

ðk0þ k1 rÞeðDk0=2Þ; B1¼ 0: ð30Þ

At this order, the temperature gradient k1 at the

inter-face is given by k1¼ dH1 dg _g¼1¼ ðA0D=3Þðk0 rÞeðDk0=2Þ 1 ðA0D=3ÞeðDk0=2Þ : ð31Þ

In Eq. (31), k1 is an explicit function of k0 and A0, so is

the integration constant A1. Thus having determined the

values of k0 and A0 from the lowest-order solution, all

temperature gradients knand integration constants Anof

the subsequent higher-order solutions can be obtained explicitly in terms of theses fundamental values.

In the present study, the solution has been deter-mined up to Oðs4_{Þ. We shall omit their expressions for}

brevity.

4. Results and discussion

4.1. Temperature proﬁles

We ﬁrst deﬁne temperature solutions of increasingly higher orders by a sequence Tm

Tm¼

Sm

1 gs; ð32Þ

where Sm¼P m

n¼0snHn denotes the ﬁrst m-term partial

sum of the series solution (18). Successive Tm indicate

that terms of next order of s are added to the solution series. In our present study, the value of m is from 0 to 4. Fixing the parameter values at D¼ 0:1, r ¼ 0:5, Figs. 2 and 3 plot the sequences of temperature proﬁles Tm at

time s¼ 0:1 and 0.75 (which correspond to 10% and 75% of solidiﬁcation of the sphere), respectively. It can be seen that the series converges very fast at small time. The convergence rate of the series solution deteriorates gradually for large s, as shown in the insert of Fig. 3, which is a local magniﬁcation of the temperature dis-tribution. Owing to the slow convergence, a large number of terms are needed in the expansion (18) to acquire a solution with satisfactory accuracy.

In order to improve the convergence rate of the series solution so that it can be applied to a much broader range of s, a nonlinear transformation proposed by

Shanks [17] is adopted. Applying the ﬁrst-order Shanks transformation results in a new sequence eð1Þ

m which is

related to the original Tmby the following formula:

eð1Þ_m ðT Þ ¼ Tmþ1

ðTmþ2 Tmþ1ÞðTmþ1 TmÞ

ðTmþ2 Tmþ1Þ ðTmþ1 TmÞ

;

mP0: ð33Þ

The same transformation rule can be again applied on the newly formed sequence eð1Þ

m to yield a second new

sequence eð2Þ m : eð2Þ_m ðT Þ ¼ eð1Þ_mþ1 eð1Þ_mþ2 eð1Þ_mþ1 eð1Þ_mþ1 eð1Þ m eð1Þ_mþ2 eð1Þ_mþ1 eð1Þ_mþ1 eð1Þm ; mP0: ð34Þ

Expression (34) is termed the second-order Shanks transformation. The process can be continued to

Fig. 2. Sequence of instantaneous temperature proﬁles Tmat

s¼ 0:1.

Fig. 3. Sequence of instantaneous temperature proﬁles Tmat

(6)

produce still higher-order Shanks transformation se-quences. In the present case, the original temperature sequence Tm is up to m¼ 4 only; which allows for at

most the employment of the second-order Shanks transformation. Connections between the original se-quence Tmand the sequences generated by applying the

ﬁrst- and second-order Shanks transformations, eð1Þ m ðT Þ

and eð2Þ

m ðT Þ are most instructively represented by the

following triangular arrays: T0 T1 T2 eð1Þ0 ðT Þ T3 eð1Þ1 ðT Þ T4 eð1Þ2 ðT Þ e ð2Þ 0 ðT Þ:

According to Shanks, the sequence formed by elements along the diagonal line has a much better convergence property than the original one. In most cases, we may even obtain an induced convergent sequence from an originally divergent one through the Shanks transfor-mation. Notion of employing the rational fraction (33) to approximate a geometric series and the convergence property of the transformed diagonal sequence were discussed in detail by Shanks [17]. Figure 4 shows at s¼ 0:75 the sequence of temperature proﬁles constructed from the diagonal elements T0, eð1Þ0 ðT Þ, and e

ð2Þ

0 ðT Þ of the

above triangular arrays. Clearly, the convergence rate is improved in this new sequence. Figure 5 compares at s¼ 0:75 the original sequence of temperature proﬁle Tm

with that obtained from applying the second-order Shanks transformation, eð2Þ₀ ðT Þ. If we consider eð2Þ0 ðT Þ a

more accurate result, apparently, the original solution (18) ought to be expanded to a much higher order of s.

Instantaneous temperature profiles for different sur-face tension values r and Stefan numbers D are plotted in Figs. 6 and 7. All temperature profiles shown in these

plots are that from the second-order Shanks transfor-mation, eð2Þ₀ ðT Þ. In reading these plots, it is reminded that in the deﬁnition of the dimensionless temperature T (Eq. (6)), T

s Tf is a negative quantity (since it is a

freezing process). Thus the larger the value of the di-mensionless temperature T is, the lower is the physical temperature T_{. Figure 6 indicates that the physical}

temperature proﬁle T_{has a larger value everywhere in}

the solidified layer for larger surface tension value. This is not surprising on viewing the interface condition (9), which states that the action of the surface tension (under fixed Stefan number) is to increase the equilibrium temperature at the interface. The same result is found in Fig. 7 where the surface tension b is held fixed; the di-mensional temperature profile T _{has a higher value for}

smaller Stefan number D.

Fig. 4. Sequence of temperature proﬁles constructed from ap-plying the ﬁrst- and second-order Shanks transformations.

Fig. 5. Comparison of the sequence of temperature proﬁles Tm

with the more accurate result obtained from applying the sec-ond-order Shanks transformation.

Fig. 6. Instantaneous temperature proﬁles at s¼ 0:5 under diﬀerent surface tension values.

(7)

4.2. Freezing rate

Another quantity that is of interest is the freezing rate, i.e., the growth rate of the interface,

dR dt ¼ D s k0þ sk1þ s2k2þ s3k3þ s4k4þ 1 s sr ð1 sÞ2 ! : ð35Þ

As before, we may deﬁne a sequence representing dif-ferent orders of the freezing rate solutions by succes-sively adding the next higher order terms into the expression R0m¼ D s S0 m 1 s sr ð1 sÞ2 ! ; ð36Þ where S0 m¼ Pm

n¼0snkn denotes the ﬁrst m-term partial

sum of the temperature gradient kn obtained from the

series solution (18). Once again, we may apply the ﬁrst-and second-order Shanks transformations to achieve better-converged results from the above sequence R0

m. In

analogy to Eqs. (33) and (34), we have

eð1Þ m ðR0Þ ¼ R0mþ1 R0_mþ2 R0 mþ1 R0_mþ1 R0 m R0 mþ2 R0mþ1 R0 mþ1 R0m ; mP0; ð37Þ eð2Þ m ðR0Þ ¼ e ð1Þ mþ1 eð1Þ_mþ2 eð1Þ_mþ1 eð1Þ_mþ1 eð1Þ m eð1Þ_mþ2 eð1Þ_mþ1 eð1Þ_mþ1 eð1Þm ; mP0: ð38Þ

Figure 8 shows the evolution history of the freezing rate R0

m with parameters values D¼ 0:1 and r ¼ 0:5,

Ewhile plotted along is the improved result obtained from the second-order Shanks transformation eð2Þ₀ ðR0_Þ

for comparison. Convergence of the sequence R0 mis quite

good except at large s. Slow convergence of the original sequence R0

mat large time is shown in the insert of Fig. 8.

Figure 9 gives the history of the freezing rate (after Shanks transformation) under various surface tension values. As indicated by the figure, raising the surface tension at the interface expedites the freezing rate jdR=dtj. In all cases demonstrated here, the freezing rate drops sharply at the early stage of solidification, and then becomes leveling off progressively. It is found that at certain critical time the freezing rate reaches a mini-mum, after that, it rises again. Basically, the freezing rate in this problem is governed by two mechanisms. Thickening of the solidified layer reduces the tempera-ture gradient across the layer and hence reducing the growth rate of the interface; while the ever rising of the

Fig. 8. Evolution history of the freezing rate solution sequence R0

m.

Fig. 9. Evolution history of the freezing rate under various surface tension conditions.

Fig. 7. Instantaneous temperature proﬁles at s¼ 0:75 under diﬀerent Stefan number values.

(8)

equilibrium temperature at the interface as the freezing front progresses towards the center of the sphere (as given by the Gibbs–Thomson condition (9)) causes an increase in the temperature gradient at the interface and hence increasing the freezing rate. The latter effect is more pronounced near the completion of freezing when the radius of the interface is close to zero. Competition of these two opposite effects results in the observed dropping and rising of the freezing rate at the initial and final phases of solidification, respectively.

4.3. Asymptotic results of freezing rate for D 1 Figure 8 suggests that the temperature-gradient series S0

m in (36) can be truncated at OðsoÞ term, the resulting

expression for the freezing rate

dR dt R 0 0¼ D s k0 1 s " sr ð1 sÞ2 # ð39Þ

is still a good approximation, at least up to a moderate value of time s. It is recalled that the k0 in (39) is an

implicit function of D through Eq. (24). In most freezing and casting problems, the Stefan number D is usually small (D Oð102_{Þ), which enables the derivation of the}

leading-order expression of k0 for D 1. After

ex-panding the exponential and error functions at the right-hand-side of Eq. (24) for small D, we have, correct to OðDÞ, k0ﬃ ð1 þ rÞ þ D 3ð1 þ rÞ 2 : ð40Þ

On substituting (40) into Eq. (39), the analytical ap-proximation of the growth rate is obtained as

dR dt ﬃ D s D 3ð1 sÞð1 ( þ rÞ2 1 ð1 sÞð1 þ rÞ sr ð1 sÞ2 ) : ð41Þ

Depending on the relative orders of magnitude of the parameters D and s, we may further distinguish the following three cases:

(I) s D 1 (approximation 1)

In this parameter range, Eq. (41) is dominated by the leading-order term dR dt ﬃ D sð1 þ rÞ: ð42Þ (II) s D 1 (approximation 2)

In this parameter range, Eq. (41) is simpliﬁed to (correct to OðDÞ and OðsÞ)

dR dt ﬃ D s ð1 þ sÞð1 þ rÞ þ sr D 3ð1 þ rÞ 2 : ð43Þ (III) D 1 s (approximation 3)

In this parameter range, the leading-order approxi-mation of (41) is (correct to OðDÞ)

dR dt ﬃ D s 1 1 sð1 ( þ rÞ þ sr ð1 sÞ2 ) ¼ D s 1þ r s ð1 sÞ2 ( ) : ð44Þ

Eqations (42)–(44) provide in each diﬀerent range of time an asymptotic expression of the freezing rate in terms of relevant physical parameters D, r and s under the condition of small Stefan number D 1. For in-stance, to leading-order approximation, the above as-ymptotic results predict a linear dependence of the freezing rate on the surface tension parameter r for all time ranges. The three expressions are plotted in Fig. 10 for D¼ 0:01, r ¼ 0:5. In each plot, more accurate result from the second-order Shanks transformation of the original sequence R0

m is also provided for comparison.

Clearly, both results are in excellent agreement, indi-cating that Eqs. (42)–(44) are indeed good approxima-tions to the freezing rate in the respective parameters ranges. Particularly we note that in Fig. 10(c), the va-lidity of the approximation 3 under the designated pa-rameters values is shown to extend at least up to a time as large as s¼ 0:75. In fact this can be easily justiﬁed by glancing at Eq. (31), we immediately come to the con-clusion k1 OðDÞ (and so for the subsequent kn, n P 1).

Thus Eq. (44), which is derived from applying the one-term approximation to the series in (35), is theoretically the correct leading-order expansion of the freezing rate for all ﬁxed s ðs < 1Þ under the asymptotic condition D! 0.

It has been mentioned previously that the growth rate of the interface exhibits a local minimum at certain critical time sc whose value depends on the

surface tension parameter r and the Stefan number D. Again, such a relationship can be derived from the approximate form of the freezing rate, Eq. (39). The condition for a local extreme to occur in the growth rate is d ds dR dt ¼ D s2_{ð1 sÞ}3 k0ð sÞ þ 2k1 0sð1 sÞ 2rs2 ¼ 0:

The above criterion leads to the ﬁnding of the roots to the following quadratic equation:

s2 3k0 2ðk0þ rÞ

sþ k0 2ðk0þ rÞ

¼ 0:

Thus the critical time sc at which the freezing rate is

(9)

sc¼ 3k0 4ðk0þ rÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k2₀ 8rk0 =ðk0þ rÞ 2 q 4 ð45Þ

in which the solution with ‘+’ sign is to be discarded because it gives a time scgreater than 1 (i.e. greater than

the total time needed for a complete solidiﬁcation). Substituting the asymptotic expression of k0 (Eq. (40))

into Eq. (45) and expanding it for D 1, we obtain the approximate expression for the critical time sc (correct

to OðDÞ) sc¼ 1 4 3ðr h þ 1Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðr þ 1Þð9r þ 1ÞiþD 4rðr þ 1Þ ðr þ 1Þ ð9r þ 5Þ 3ð9r þ 1Þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðr þ 1Þð9r þ 1Þ p : ð46Þ Equation (46) is graphed in Fig. 11 for various sur-face tension parameter values b (b¼ rD) and Stefan number D. It is observed that critical freezing rate occurs at an earlier time for larger surface tension value and

(a) (b)

(c)

Fig. 10. Comparisons between the asymptotic results of the freezing rate and that obtained from the original series solution (with Shanks transformation) in the parameters ranges (a) s D 1, (b) s D 1, (c) D 1 s.

Fig. 11. Critical time sc for the freezing rate as a function of

(10)

smaller Stefan number – a situation corresponding to a higher equilibrium temperature (dimensional) at the solid/liquid interface (see the boundary condition (9)) and hence a higher temperature gradient at the freezing front. In all cases we have considered here, minimal freezing rate is found to develop in the ﬁrst half of the solidiﬁcation process, i.e., sc60:5. From Fig. 11, it is

also noticed that when b¼ 0; the minimal freezing rate always occurs at sc¼ 0:5, regardless the value of the

Stefan number (as long as D 1). This is already im-plied in Eq. (39) as r¼ 0 is substituted into the equa-tion; the resulting form of dR=dt has an extreme at s¼ 0:5, independent of D.

5. Concluding remarks

A small-time series expansion technique has been applied in this work to study the thermal effect of the surface tension on an inward solidification process with spherical symmetry. A nonlinear Shanks transformation has been adopted to improve the convergence rate of the series solution so as to extend its applicability to a larger time range. In summary, the present result shows that surface tension increases the equilibrium temperature at the solid/liquid interface and hence speeds up the growth rate of the freezing front. The asymptotic expression of the freezing rate has been obtained for each different time range under the condition of small Stefan number. To leading-order approximation, the results predict a linear increase of the freezing rate with the surface ten-sion parameter for all time ranges under consideration. Critical time at which the freezing rate develops a local minimum has also been determined analytically, which shows that the freezing rate reaches its minimum value more quickly for larger surface tension value at the in-terface.

The form of the Gibbs–Thomson law assumed in the present analysis is not applicable for solid/liquid inter-face with extremely small radius of curvature. Thus it is expected that the validity of the present result cannot be extended to the time when the freezing front comes very close to the center of the sphere. In practical applica-tions, the upper limit of such time depends on the numerical values of various parameters. A simple esti-mation of a typical value of this limitation time is given in the appendix. In practice, the limitation is not that stringent as it may appear. Aside from this, the present study does provide critical information regarding to the variation of the freezing rate with relevant physical parameters as the thermal effect of the surface tension at the interface is taken into consideration. It is clear from the present analysis that surface tension has nonnegli-gible influence on the freezing rate. This is particularly true at the initial stage (s 0) and close to the final phase (s 1) of solidification where the freezing rate is

shown to be augmented by a factor r=s and r=ð1 sÞ2, respectively.

Appendix A

To estimate the typical value of s beyond which the present analysis is not applicable, we start from the Gibbs–Thomson law, Eq. (1). Already, it has been stated in Section 1 that Eq. (1) holds only if e¼ 2c=LR_1.

Rewrite the condition using nondimensional parameters rand R deﬁned in the text, we obtain

e¼ rðT f T sÞ=ðRT fÞ 1: ðA:1Þ Let us take r¼ 0:1, T

f ¼ 273 K; and assume a small

degree of cooling,ðT

f TsÞ ¼ 5 K for example. Now, if

we accept that 0.05 is practically a small enough number for e, condition (A.1) then gives

s¼ 1 R 0:96;

which stands for a 96% of solidiﬁcation of the sphere. Moreover, the smaller the surface tension parameter r, the closer this value is to 1. Thus in practice, limitation on the applicability of the Gibbs–Thomson law in the form of Eq. (1) is not that restricted as it may appear. In most cases involving small surface tension, small Stefan number or small degree of cooling, we can push the limit of usefulness of the present series solution to a point not too far from a complete solidiﬁcation.

References

[1] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Clarendon Press, Oxford, 1959.

[2] J.M. Hill, One-Dimensional Stefan Problems: An Introduction, Longman Scientiﬁc and Technical, Essex, 1987.

[3] F. Kreith, F.E. Romie, A study of the thermal diﬀusion equation with boundary conditions corresponding to solidiﬁcation or melting of materials initially at the fusion temperature, Proc. Phys. Soc. 68B (1955) 277–286. [4] D. Langford, The freezing of spheres, Int. J. Heat Mass

Transfer 9 (1966) 827–836.

[5] G. Poots, On the application of integral methods to the solution of problems involving the solidiﬁcation of liquids initially at fusion temperature, Int. J. Heat Mass Transfer 5 (1962) 525–531.

[6] L.C. Tao, Generalized numerical solutions of freezing a saturated liquid in cylinders and spheres, AIChE J. 13 (1) (1967) 165–169.

[7] N.-Y. Li, Thermomechanical stresses and some asymptotic behavior in castings with spherical solidiﬁcation, J. Ther-mal Stresses 18 (1995) 165–184.

[8] J. Caldwell, C.-C. Chan, Numerical solutions of the Stefan problem by the enthalpy method and the heat balance integral method, Numer. Heat Transfer Part B 33 (1998) 99–117.

(11)

[9] C.L. Huang, Y.P. Shih, A perturbation method for spherical and cylindrical solidiﬁcation, Chem. Eng. Sci. 30 (1975) 897–906.

[10] R.I. Pedroso, G.A. Domoto, Perturbation solutions for spherical solidiﬁcation of saturated liquids, ASME J. Heat Transfer 95 (1973) 42–46.

[11] R.I. Pedroso, G.A. Domoto, Inward spherical solidiﬁca-tion–solution by the method of strained coordinates, Int. J. Heat Mass Transfer 16 (1973) 1037–1043.

[12] D.S. Riley, F.T. Smith, G. Poots, The inward solidiﬁcation of spheres and circular cylinders, Int. J. Heat Mass Transfer 17 (1974) 1507–1516.

[13] K. Stewartson, R.T. Waechter, On Stefan’s problem for spheres, Proc. R. Soc. London Ser. A 348 (1976) 415–426.

[14] G.B. Davis, J.M. Hill, A moving-boundary problem for the sphere, IMA J. Appl. Math. 29 (1982) 99–111.

[15] J.M. Hill, A. Kucera, Freezing a saturated liquid in-side a sphere, Int. J. Heat Mass Transfer 26 (1983) 1631–1637.

[16] D.T.J. Hurle (Ed.), Handbook of Crystal Growth Volume I: Fundamentals, North-Holland, Amsterdam, 1993, pp. 194–197.

[17] D. Shanks, Nonlinear transformations of divergent and slowly convergent sequences, J. Math. Phys. 34 (1955) 1–42.

[18] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Math-ematical Functions, Dover, New York, 1965.