- Email: [email protected]
Solidification of a finite slab with convective cooling and shrinkage
Applied Mathematical Modelling 27 (2003) 733–762 www.elsevier.com/locate/apm
Solidification of a finite slab with convective cooling and shrinkage Zhaoji Yang 1, Mihir Sen *, Samuel Paolucci Department of Aerospace and Mechanical Engineering, University of Notre Dame, Notre Dame, IN 46556-5637, USA Received 29 April 2002; received in revised form 12 March 2003; accepted 28 March 2003
Abstract We consider the one-dimensional solidification of a finite slab of a pure substance which is initially in liquid state. Convective cooling is applied at one end of the slab while the other end is kept adiabatic. There is a precooling stage before solidification begins. A void, which has substantial overall effect on the solidification of the slab, is formed at the cooled end due to shrinkage and grows during solidification. Perturbation analysis for small Stefan numbers is used to find the interface position and temperature distributions in both solid and liquid phases. The order of magnitude of the Biot number at the cold end, Bi, determines the nature of the solution. Regular and two time-scale perturbation solutions are obtained for Bi Oð1=2 Þ and Bi Oð1Þ, respectively. The results can be extended to other orders of magnitude, as well as to slowly-varying time-dependent Bi. The results agree well with special cases in the literature, and with numerical solutions calculated using the enthalpy method. 2003 Elsevier Inc. All rights reserved.
1. Introduction Heat transfer problems involving phase change such as solidification and melting are of great importance in engineering applications [1]. The mathematical model is a system of partial differential equations of the parabolic type with moving boundaries. Difficulties are encountered in the analytical solution of these problems due to the nonlinearity of moving boundary conditions.
*
Corresponding author. Tel.: +1-574-631-5975; fax: +1-574-631-8341. E-mail address: [email protected] (M. Sen). 1 Currently at Advanced Technology Development, Ford Visteon, Electronics Technical Center, 17000 Rotunda Dr., Dearborn, MI 48121-6010. 0307-904X/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0307-904X(03)00078-7
734
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
Nomenclature Bi c h k ‘ L s t T Tf T0 Ti x Greek a b c f g h H k kv l n q r s / W
¼ hl=ks , Biot number specific heat heat transfer coefficient thermal conductivity thickness of slab latent heat of fusion solid–liquid interface position time temperature freezing temperature mold temperature initial temperature coordinate along slab symbols ¼ k=qc, thermal diffusivity ¼ 1 cð2=kv 1Þ ¼ 1 ql =qs , shrinkage factor Stefan number nondimensional coordinate for numerical method nondimensional coordinate in liquid nondimensional temperature in solid nondimensional temperature for numerical method ¼ kl =ks ¼ kv =ks ¼ al =as nondimensional coordinate in solid density ¼ s=‘, nondimensional interface position nondimensional time nondimensional temperature in liquid nondimensional enthalpy
Subscripts and superscripts l liquid s solid m matching v void outer variable precooling stage * beginning of solidification
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
735
Exact solutions are possible in similarity form only in limited cases. Many approximate and numerical methods of solution for nonsimilar cases have also been developed [2–4]. The effect of shrinkage and void formation due to solidification was studied by Wilson and Solomon [5] who obtained a similarity solution by assuming the void thermal resistance proportional to its size. Perturbation methods based on the assumption of a small Stefan number have also been used to deal with finite domain solidification problems. There are two types of problems that have been considered: one where the liquid is initially superheated above the freezing temperature and in which the temperature in both solid and liquid must be determined (two-phase), and the other where the liquid is not superheated but is always at the freezing temperature (one-phase). A regular perturbation analysis of the one-phase problem for convective and radiative cooling has been carried out by Yan and Huang [6]. The radiative boundary condition which they presented was simply linearized to convective form. They could use the interface position as an independent variable instead of time since the interface position increases monotonically with respect to time. In their conclusions, they foresaw some difficulty in applying the perturbation technique to the convective solidification of an initially superheated slab. However, their analysis did not include the effect of shrinkage. Significant progress for the two-phase problem in a finite domain was reported by Weinbaum and Jiji [7]. Cooling was by an isothermal wall below freezing temperature while the other end was either adiabatic or isothermal. Shrinkage was not included here either. The problem was solved using the method of matched asymptotics. There has been no analytical treatment, however, of the finite domain solidification problem of a two-phase region due to convective cooling. This kind of boundary condition is fairly general and includes the isothermal as a limiting case. The presence of the shrinkage void is an important factor in the solidification process. We will treat this problem using perturbation methods. The solutions thus obtained will be compared to previous results under appropriate limits and also verified by numerical computation. The geometry is taken to be one-dimensional as shown in Fig. 1. The molten material at temperature Ti is initially between x ¼ 0 and x ¼ ‘. Cooling is then applied at x ¼ 0, which could represent the surface of the mold where solidification is taking place. In the precooling stage, the heat transfer coefficient is denoted as h~, which includes effects of thermal contact between the liquid and the mold, and then convective cooling to the mold at temperature T0 . A constant value of h~ is assumed throughout the precooling stage. Eventually the temperature of the liquid at x ¼ 0 x=0
x=γs(t)
x=s(t)
x=l
Convective Cooling
Void
Liquid
Solid motion
T=T0
motion
T=Tf
x Fig. 1. Schematic of the physical problem.
Adiabatic
736
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
comes down to the freezing point and solidification begins. As the material solidifies, it shrinks and a void is formed between the solidifying material and the mold surface, as seen in experiments [8]. This is the region between x ¼ 0 and x ¼ csðtÞ, where c ¼ ð1 ql =qs Þ is the shrinkage parameter quantifying the effect of density change. It must be mentioned that in reality void formation often occurs in two different ways: through micropores within the material and what is referred to as a shrinkage gap. The proportion of one to the other depends greatly on external conditions. We assume here that it is entirely of the second kind. Note also that while we are considering a one-dimensional analysis, the result should be locally applicable to higher dimensional problems as long as the characteristic length scales in the transverse directions are much larger than c‘. We assume that the heat transfer across this void and the mold to the ambient can be represented by a heat transfer coefficient h. We will take h to be a constant [8], the same as h~. Moreover, the analysis in which it is varying slowly with time is considered in Section 6.4. The outer boundary of the liquid at x ¼ ‘ is taken to be adiabatic. In what follows, wherever numerical values are needed we will consider the solidification of pure aluminum.
2. Temperature distribution before solidification In a one-phase cooling problem, such as in [7], the material would begin to solidify as soon as cooling is applied because of the instantaneous drop of the temperature at x ¼ 0. In the present problem, however, there is a period of precooling before the solidification can actually start. The material is all in liquid form for time t < t , where t is the time necessary to bring the liquid at the boundary to the freezing temperature. Before looking at solidification, we would like to determine the temperature distribution in the liquid at t ¼ t . In the precooling stage, heat transfer in the liquid is by single-phase conduction. The governing equation and boundary conditions for the temperature distribution, Tl ðx; tÞ, are al
o2 Tl oTl ; ¼ ox2 ot
kl
oTl ð0; tÞ ¼ h~½Tl ð0; tÞ T0 ox
06x6‘
ð1Þ ð2Þ
oTl ð‘; tÞ ¼ 0 ox
ð3Þ
Tl ðx; 0Þ ¼ Ti
ð4Þ
These can be nondimensionalized to e e o/ o2 / ; ¼ os og2
06g61
e e Bi o/ e ð0; sÞ ð0; sÞ ¼ / k og
ð5Þ
ð6Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
737
e o/ ð1; sÞ ¼ 0 og
ð7Þ
e ðg; 0Þ ¼ 1 /
ð8Þ
where the nondimensional liquid temperature, distance and time are given by e ¼ Tl T0 / Ti T0 g¼
x ‘
ð9Þ ð10Þ
al t ð11Þ ‘2 respectively, and k ¼ kl =ks . The controlling parameter is the heat transfer coefficient which is nondimensionalized into a Biot number defined by s¼
~ e ¼ h‘ Bi ks
ð12Þ
It is easy to find a series solution to this conduction problem [1]; the resulting temperature distribution is e ðg; sÞ ¼ /
1 X
An expðpn2 sÞ cos½pn ð1 gÞ
ð13Þ
n¼0
where An ¼
2 sin pn # e k Bi pn 1 þ 2 e þ k2 p2 Bi n "
ð14Þ
and pn is the nth positive eigenvalue determined from the transcendental equation tan pn ¼
e Bi kpn
ð15Þ
The time s ¼ s when the cold end drops to the freezing point is determined implicitly from e ð0; s Þ ¼ Tf T0 / Ti T0 e ð0; s Þ is a parameter about initial preheating that must be specified. Note that / The temperature distribution in the liquid at the instant solidification begins is e ðg; s Þ ¼ /
1 X n¼0
An expðpn2 s Þ cos½pn ð1 gÞ
ð16Þ
ð17Þ
738
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
e it becomes flatter as Bi e decreases. Fig. 3 shows the This is shown in Fig. 2 for different Bi Bi; e We have assumed an aluminum slab with ‘ ¼ 0:577 m, k ¼ 0:500 and Bi. variation of s with Bi freezing temperature Tf ¼ 660 C. Other temperatures are Ti ¼ 710 C and T0 ¼ 640 C, thus e ð0; s Þ ¼ 0:286. Typically, a precooling heat transfer coefficient h~ ¼ 364 W/m2 K corresponds to / e ¼ 1:0. a Biot number Bi 1 Bi = 10 6 4.1 0.8
3 2
0.6 1
0.4
φ∗=0.286 0.2
0
0.2
0.4
0.6
0.8
1
e ð0; s Þ ¼ 0:286. Fig. 2. Temperature distribution in the liquid at beginning of solidification; k ¼ 0:500, /
10
10
1
0
φ∗=0.5
–1
10
φ∗=0.286
φ∗=0.1
–2
10
–3
10
0
2
4
6
8
e k ¼ 0:500. Fig. 3. Variation of preheating time s with Bi; Bi
10
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
739
3. Governing equations during solidification Once solidification begins we have to deal with both liquid and solid phases. Due to shrinkage there is a void formed which, consistent with some experiments [8], we assume to be near the cold end. The growth of the void leads to bulk motion of the solid phase, adding a convective term to the heat conduction equation in the solid. It is assumed, however, that there is no convective motion in the liquid [5]. The governing equations are thus as
o2 Ts oTs ds oTs þc ; ¼ 2 dt ox ox ot
al
o2 Tl oTl ¼ ; ox2 ot
t P t ;
t P t ; s6x6‘
cs 6 x 6 s
ð18Þ ð19Þ
for each of the two phases separated by the interface at x ¼ sðtÞ. The boundary and initial conditions are ks
oTs ðcs; tÞ ¼ h½Ts ðcs; tÞ T0 ox
ð20Þ
Ts ðs; tÞ ¼ Tf
ð21Þ
Tl ðs; tÞ ¼ Tf
ð22Þ
ks
oTs oTl ds ðs; tÞ kl ðs; tÞ ¼ ql L dt ox ox
ð23Þ
oTl ð‘; tÞ ¼ 0 ox
ð24Þ
sðt Þ ¼ 0
ð25Þ
e ðx; t Þ Tl ðx; t Þ ¼ T0 þ ðTi T0 Þ /
ð26Þ
e ðx; t Þ is calculated from Eq. (17). We introduce the liquid and solid nondimensional where / temperature variables /¼k h¼
Tl Tf Tf T0
Ts Tf Tf T0
ð27Þ ð28Þ
A Landau transformation is used to convert the moving boundary problem into one of fixed domain. The new spatial variables in the solid and liquid are n¼
x s
ð29Þ
740
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
xs ð30Þ ‘s where s is defined as before. The nondimensional position of the solid–liquid interface is r ¼ s=‘, and the thermal diffusivity ratio is l ¼ al =as . The governing equations in these variables are o2 h dr oh 2 oh ¼ l r ð31Þ ðn cÞr ; s P s ; c 6 n 6 1 2 os ds on on g¼
o2 / o/ dr o/ ¼ ð1 rÞ2 ð1 gÞð1 rÞ ; 2 og os ds og
s P s ; 0 6 g 6 1
ð32Þ
with boundary conditions hð1; sÞ ¼ 0
ð33Þ
/ð0; sÞ ¼ 0
ð34Þ
o/ ð1; sÞ ¼ 0 og oh o/ dr ð1 rÞ ð1; sÞ r ð0; sÞ ¼ rð1 rÞ lð1 cÞ on og ds
ð35Þ ð36Þ
oh ðc; sÞ ¼ rBi½hðc; sÞ þ 1 on
ð37Þ
rðs Þ ¼ 0
ð38Þ "
e ðg; s Þ / 1 /ðg; s Þ ¼ k e ð0; s Þ /
#
ð39Þ
The nondimensional parameters controlling the solidification process are ¼ cs Bi ¼
Tf T0 L
h‘ ks
ð40Þ ð41Þ
where is the Stefan number and Bi is the Biot number. The Stefan number is the ratio of sensible heat to latent heat. The Biot number represents the ratio between the conductive thermal resistance in the solid and the thermal resistance across the void to the mold. The thermal diffusion time scale in the liquid is tdiff ¼ ‘2 =al (the diffusion time in the solid is of the same order), and the void heat transfer time scale is tvoid ¼ ‘ql cl =h. They are related by ktdiff =tvoid ¼ Bi. In this paper we have chosen to nondimensionalize all times by tdiff . Note that the solidification time scale is a dependent variable given by tdiff =Bi. In addition, it is possible that there may be time scales associated with changes in e h , h and T0 . The solutions developed here will remain valid as long as they are slowly varying; more will be said about this in Section 6.4.
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
741
Bi is proportional to h, which is the heat transfer coefficient after the separation of the solidified e on the other hand is its value before solidification begins, i.e. material from the mold wall. Bi Bi, e Generally, we can before the void is formed. We will use the same constant value for Bi as for Bi Bi. e expect that BiðsÞ 6 Bi because of the increasing insulating effect of the void. The transition should e at the instant solidification starts (see the model presented in be continuous with Biðs Þ ¼ Bi Section 6.4). When Bi ! 1, the boundary condition (37) reduces to the isothermal condition analyzed by Weinbaum and Jiji [7]. Here we will investigate the problem for finite Bi. We will also assume that l, k and ð1 cÞ are all fixed by material properties and are of Oð1Þ. The Stefan number only appears in the interface condition (36). It is small in many practical applications, and we will use it as a small parameter to develop a perturbation solution. The analysis will depend on the order of magnitude of the Biot number. We will consider two different cases, Bi Oð1=2 Þ and Bi Oð1Þ, with which we can cover the entire range of Bi. The solidification problem can also be solved numerically to check the validity of the perturbation solutions. We use an enthalpy method in which the solid and liquid regions are treated as a continuum, except for an enthalpy jump at the interface [10]. To implement the numerical method, we introduce f¼
x csðtÞ ; lð1 cÞ
2
W ¼ lð1 cÞ H;
H¼
1 Ts Tf kð1 cÞ Tf T0
ð42Þ
in the solid and x f¼ ; l
W¼Hþ
lð1 cÞ ; k
H¼
Tl Tf Tf T0
ð43Þ
in the liquid. The equations are nondimensional; the convective term has also been removed. The original system of Eqs. (18) and (19) and boundary conditions (20)–(26) reduce to o2 H oW ¼ ;06f61 os of2
ð44Þ
oH ð1; sÞ ¼ 0 of
ð45Þ
oH Bi ð0; sÞ ¼ ð1 cÞBiHð0; sÞ þ of k
ð46Þ
e ðf; s Þ / 1 e ð0; s Þ /
ð47Þ
Hðf; s Þ ¼
The problem is solved using a finite difference method with second-order central difference in space and explicit Euler marching in time [9]. The numerical values for aluminum are c ¼ 0:0627, k ¼ 0:50, and l ¼ 0:582. For the temperatures assumed before, we have ¼ 0:0594.
742
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
4. Solution for Bi O (1=2 ) Let us first consider a large heat transfer coefficient corresponding to Bi Oð1=2 Þ. The diffusion time is large compared to the void heat transfer time, and the solidification process is diffusion limited. We seek perturbation solutions of the form hðn; sÞ ¼
1 X
n=2 hn ðn; sÞ
ð48Þ
n¼0
/ðg; sÞ ¼
1 X
n=2 /n ðg; sÞ
ð49Þ
n¼0
rðsÞ ¼
1 X
n=2 rn ðsÞ
ð50Þ
n¼0
We follow the standard procedure of substituting these expansions into the governing equations (31) and (32) and boundary and initial conditions (33)–(39), and balancing terms in each order. We also introduce the Oð1Þ parameter Bi0 Bi 1=2 . From the Oð1Þ terms in the interface condition (36), we have r0 ¼ 0
ð51Þ
Collecting terms of Oð1Þ in Eqs. (31), (33) and (37), we get o2 h0 ¼0 on2
ð52Þ
h0 ð1; sÞ ¼ 0
ð53Þ
oh0 ðc; sÞ ¼ r1 Bi0 ½h0 ðc; sÞ þ 1 on
ð54Þ
which gives us h0 ðn; sÞ ¼
r1 Bi0 ðn 1Þ r1 Bi0 ð1 cÞ þ 1
ð55Þ
Using this in the OðÞ balance of Eq. (36), we obtain r1 ðsÞ
dr1 1 oh0 ¼ ð1; sÞ lð1 cÞ on ds
ð56Þ
From this dr1 1 Bi0 ¼ lð1 cÞ 1 þ ð1 cÞBi0 r1 ds
ð57Þ
r1 ðs Þ ¼ 0
ð58Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
which can be solved as 1 r1 ¼ 1c
1 þ Bi0
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 s s þ 2 l Bi20
743
ð59Þ
From a balance of Oð1Þ terms in Eq. (32), and boundary conditions (34) and (35), we can write o2 /0 o/0 ¼ og2 os
ð60Þ
/0 ð0; sÞ ¼ 0
ð61Þ
o/0 ð1; sÞ ¼ 0 og "
ð62Þ
# e ðg; s Þ / 1 /0 ðg; s Þ ¼ k e ð0; s Þ /
ð63Þ
e ðg; s Þ as Oð1Þ contribution in the Here we include the initial liquid temperature distribution / perturbation expansion. By separation of variables, we obtain the solution # " 2 1 X 1 1 Cn sin n þ p2 ðs s Þ ð64Þ pg exp n þ /0 ðg; sÞ ¼ 2 2 n¼0 where Cn ¼
1 X 2k þ n þ 12 p m¼0
1 nþ p cos pm 2 " # 2 1 e ð0; s Þ n þ p2 pm2 / 2
2Am k expðpm2 s Þ
ð65Þ
Balancing the Oð1=2 Þ terms in Eqs. (31), (33) and (36), we find o2 h1 ¼0 on2
ð66Þ
h1 ð1; sÞ ¼ 0
ð67Þ
oh1 ðc; sÞ ¼ r2 Bi0 ½h0 ðc; sÞ þ 1 þ Bi0 r1 ðsÞh1 ðc; sÞ on
ð68Þ
We solve for h1 to get Bi0 r2 ½h0 ðc; sÞ þ 1 ðn 1Þ h1 ðn; sÞ ¼ 1 þ ð1 cÞBi0 r1 ðsÞ From a balance of the Oð3=2 Þ terms in the interface condition (36), we have 1 oh1 oh0 o/ dðr1 r2 Þ dr1 r21 ð1; sÞ r1 ð1; sÞ r1 0 ð0; sÞ ¼ lð1 cÞ on ds on og ds
ð69Þ
ð70Þ
744
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
Substituting previously obtained expressions for r1 ðsÞ, h0 ðn; sÞ, h1 ðn; sÞ, and /0 ðg; sÞ into this, we get dr2 ¼ AðsÞr2 þ BðsÞ ds
ð71Þ
r2 ðs Þ ¼ 0
ð72Þ
where AðsÞ ¼
Bi20 l þ 2Bi20 ðs s Þ
ð73Þ
# " 2 1 X 1 1 1 Cn n þ p2 ðs s Þ BðsÞ ¼ p exp n þ lð1 cÞ n¼0 2 2
ð74Þ
The general solution is of the form r2 ðsÞ ¼ exp
Z
s
Aðs0 Þ ds0
s
Z
s
" Bðs00 Þ exp
s
Z
s00
# Aðs0 Þ ds0 ds00
ð75Þ
s
which can be evaluated as " # Z s sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 2 2 0 1 Cn n þ p2 ðs0 s Þ ds0 1 þ Bi0 ðs s Þ exp n þ p 2 l 2 s n¼0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 ðsÞ ¼ 2 lð1 cÞ 1 þ Bi20 ðs s Þ l " # 8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 1 > > 2 > > 1 < 1 1 þ l Bi0 ðs s Þ exp n þ 2 p ðs s Þ X 1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Cn ¼ > 1 2 2 > n¼0 > nþ p lð1 cÞ 1 þ Bi0 ðs s Þ > : 2 l " 9 2 # rffiffiffi l n þ 12 p2 2 > > Bi0 exp > > rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi = 1 2Bi20 l l n þ p 1 l 2 p ffiffi ffi p ðs s Þ þ erf n þ erf þ 2 > 2 2Bi20 2Bi0 1 > > > p3=2 nþ ; 2 1 X
ð76Þ
The perturbation solution for rðsÞ has been obtained up to OðÞ. Higher order solutions are difficult to obtain in explicit analytical form because of the nonlinearity of the corresponding differential equations. Even so, we can obtain good accuracy with the first few perturbation terms as long as the parameter is small enough. That the perturbation series converges fast can be seen
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
745
4
σ1
3
1/2
σ+ε σ2 1 2
−σ 2
1
0
0
2
4
6
e ¼ Fig. 4. Comparison of first- and second-order Bi Oð1=2 Þ perturbation solutions for interface location; Bi ¼ Bi e ð0; s Þ ¼ 0:286, c ¼ 0:0627, l ¼ 0:582, ¼ 0:0594. 4:10, k ¼ 0:500, /
from comparing the terms r1 and r2 in Fig. 4. The second-order solution is only important near the beginning of solidification. One can easily verify that since r1 and r2 are Oð1Þ for Bi Oð1=2 Þ, the expansion remains uniformly valid. Fig. 5 shows temperature histories at the the two ends of the slab for Bi ¼ 1=2 . One term is used in Eq. (48) for hðc; sÞ and two in (49) for /ð1; sÞ. The perturbation solution is seen to compare well with the numerical solution. 5. Solution for Bi O(1) This corresponds to a moderate heat transfer coefficient h. Both diffusion and void time scales are of the same order. 5.1. Short-time solution Following the previous form of the perturbation expansions (48)–(50), we find that r2nþ1 ¼ 0. Actually, in both the perturbation expansion of nondimensionalized temperature and solid–liquid interface position function, the odd n terms drop out. Thus we can simply use the expansions 1 X n hn ðn; sÞ ð77Þ hðn; sÞ ¼ n¼0
/ðg; sÞ ¼
1 X n¼0
n /n ðg; sÞ
ð78Þ
746
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762 1.5
1
0.5
φ(1,τ) 0
–0.5
θ(γ,τ) –1
–1.5
0
1
2
3
4
5
6
7
8
Fig. 5. Comparison of numerical and Bi Oð1=2 Þ perturbation solutions for temperature history at boundaries; solid e ¼ 4:10, lines are from perturbation solution, and dashed lines with points correspond to numerical solution; Bi ¼ Bi e ð0; s Þ ¼ 0:286, c ¼ 0:0627, l ¼ 0:582, ¼ 0:0594. k ¼ 0:500, /
rðsÞ ¼
1 X
n rn ðsÞ
ð79Þ
n¼0
On gathering the terms of order Oð1Þ, we find that r0 ¼ 0 /0 ¼
1 X
ð80Þ Cn sin
n¼0
# " 2 1 1 nþ p2 ðs s Þ pg exp n þ 2 2
ð81Þ
and o2 h0 ¼0 on2
ð82Þ
h0 ð1; sÞ ¼ 0
ð83Þ
oh0 ðc; sÞ ¼ 0 on
ð84Þ
which means that h0 ðn; sÞ ¼ 0 The balance of OðÞ terms gives o2 h1 ¼0 on2
ð85Þ
ð86Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
747
h1 ð1; sÞ ¼ 0
ð87Þ
oh1 ðc; sÞ ¼ r1 Bi on Solving, we get
ð88Þ
h1 ðn; sÞ ¼ r1 Biðn 1Þ
ð89Þ
2
From the balance of Oð Þ terms in Eq. (36), we obtain dr1 1 oh1 o/ ¼ ð1; sÞ r1 0 ð0; sÞ r1 lð1 cÞ on ds og
ð90Þ
which gives
Thus
" )# ( 2 1 X dr1 1 1 1 Bi ¼ Cn n þ p2 ðs s Þ p exp n þ lð1 cÞ 2 2 ds n¼0
ð91Þ
r1 ðs Þ ¼ 0
ð92Þ
( " 1 X 1 C n r1 ¼ exp Biðs s Þ þ 1 lð1 cÞ n þ p n¼0 2
1 nþ 2
2
! p2 ðs s Þ
)# 1
ð93Þ
The Oð2 Þ terms in Eqs. (31), (33) and (37) give o2 h2 ¼0 on2
ð94Þ
h2 ð1; sÞ ¼ 0
ð95Þ
oh2 ðc; sÞ ¼ r1 Bih1 ðc; sÞ þ r2 Bi on which can be solved as h2 ðn; sÞ ¼ Bi2 ð1 cÞr21 þ Bir2 ðn 1Þ
ð96Þ
ð97Þ
The Oð3 Þ terms of Eq. (36) give dðr1 r2 Þ 1 oh2 oh1 o/1 o/0 2 dr1 r1 ¼ ð1; sÞ r1 ð1; sÞ r1 ð0; sÞ r2 ð0; sÞ ds lð1 cÞ on ds on og og
ð98Þ
where /1 is the first-order temperature distribution in the liquid. In principle, this can be determined from o2 /1 o/1 o/ dr1 o/0 2r1 0 ð1 gÞ ¼ 2 og os os ds og o/1 ð1; sÞ ¼ 0 og
ð99Þ ð100Þ
748
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
/1 ð0; sÞ ¼ 0
ð101Þ
/1 ðg; 0Þ ¼ 0
ð102Þ
It is difficult to solve for /1 analytically due to the series form of /0 and r1 . Notice, however, the 2 homogeneous nature of the boundary conditions and the exp½ n þ 12 p2 s rate of decay with respect to time in the /0 and dr1 =ds series in Eqs. (81) and (91). We can conclude that the time rate of change of /1 ðg; sÞ and o/1 =ogðg; sÞ will have at least an exponential decay factor for large s, although it increases from zero to a finite value near the beginning. Using Eqs. (81), (89) and (97) in Eq. (98), the equation for r2 is dr2 dr1 1 Bi 1 o/1 ¼ r1 ð0; sÞ ð103Þ þ Bi r1 ds ds 1c l lð1 cÞ og which can be integrated to give Z s Z s 1 1 Bi 1 o/1 ð0; s0 Þ ds0 r1 ðs0 Þ ds0 þ Bi r2 ðsÞ ¼ r21 2 1c l s lð1 cÞ s og
ð104Þ
As s ! 1, r1 ! Bi s where the exponential term in Eq. (93) has dropped out. As discussed at the end of the last paragraph, the last integral is also negligibly small due to its exponential decay factor. The limit of r2 with large s is r2 !
Bi3 s2 2l2 ð1 cÞ
ð105Þ
We now look at the perturbation solution we have obtained up to now. As s ! lð1 cÞ=Bi, the interface position tends to r!
Bi s Bi3 s2 2 þ Oð3 Þ lð1 cÞ 2l ð1 cÞ
ð106Þ
This is not valid near r ¼ 1 since the higher order terms become comparable to those of lower order when s lð1 cÞ=Bi. So only the early stage of solidification is described well by this expression, but not the later stage. This is thus a small-time solution. 5.2. Large-time solution To obtain the large-time solution, we will follow the steps used in Section 4 and use the expansions (48)–(50). We change the time scale to bs ¼ s. Denoting the large-time solution by b r, b h, b , the governing equations (31)–(39) become and / " # b o2 b db r ob h h 2 oh ðn cÞb r ; c6n61 ð107Þ ¼ b r obs dbs on on2 b b b b o/ o2 / dr 2 o/ ð1 gÞð1 b rÞ ; ¼ ð1 b rÞ 2 obs dbs og og with boundary conditions
06g61
ð108Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
749
b h ð1; bs Þ ¼ 0
ð109Þ
b ð0; bs Þ ¼ 0 /
ð110Þ
b o/ ð1; bs Þ ¼ 0 og " # b 1 ob h o/ db r ð1 b r Þ ð1; bs Þ b r ð0; bs Þ ¼ b r ð1 b rÞ lð1 cÞ on og dbs h i ob h ðc; bs Þ ¼ b r Bi b h ðc; bs Þ þ 1 on
ð111Þ
ð112Þ
ð113Þ
From an Oð1Þ balance in Eq. (108), we obtain b o2 / 0 ¼0 og2
ð114Þ
b ð0; sÞ ¼ 0 / 0
ð115Þ
b o/ 0 ð1; sÞ ¼ 0 og
ð116Þ
which gives b ðg; sÞ ¼ 0 / 0
ð117Þ
This solution is only valid in the later stage of solidification, not near the beginning where bs ¼ bs . Eq. (117) implies that the temperature distribution in the liquid region decays to almost the freezing temperature when the outer solution becomes valid. A balance of Oð1Þ terms in Eq. (107) is h0 o2 b ¼0 on2
ð118Þ
b h 0 ð1; sÞ ¼ 0
ð119Þ
ob h0 ðc; sÞ ¼ b r 0 Bi½b h 0 ðc; sÞ þ 1 on
ð120Þ
which can be solved to give b h 0 ðn; bs Þ ¼
Bib r0 ð1 nÞ 1 þ ð1 cÞBib r0
ð121Þ
750
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
To this order, the interfacial condition (112) is db r0 1 ob h0 b r 0 ð1 b ð1 b r0Þ ¼ ð1; bs Þ r0Þ lð1 cÞ dbs on which can be reduced to
ð122Þ
db r0 1 Bi ¼ lð1 cÞ 1 þ ð1 cÞBib r0 dbs
ð123Þ
b r 0 ðbs 0 Þ ¼ 0
ð124Þ
b 0 at Since the actual initial condition is outside the range of validity of the large-time solution, r bs ¼ bs 0 is an arbitrary constant, where bs 0 will be determined later by matching with the small-time solution. The matching also shows that the constant can be taken to be zero without loss of generality. The solution of b r 0 is s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " # 1 1 1 2ðbs bs 0 Þ b ð125Þ þ r0 ¼ þ 1c Bi Bi2 l A balance of terms of Oð1=2 Þ in Eq. (108), gives b ðg; bs Þ ¼ 0 / 1 and a balance of terms of Oð
ð126Þ 1=2
Þ in Eq. (107), gives
o2 b h1 ¼0 on2
ð127Þ
hb1 ð1; sÞ ¼ 0
ð128Þ
ob h1 ðc; sÞ ¼ b r 1 Bi½b h 0 ðc; sÞ þ 1 þ Bib r0 b h 1 ðc; sÞ on
ð129Þ
The solution is Bib r 1 ½b h 0 ðc; sÞ þ 1 b h 1 ðn; bs Þ ¼ ðn 1Þ 1 þ ð1 cÞBib r0 The interfacial condition for Oð1=2 Þ is ob h1 db r1 db r0 ð1; bs Þ ¼ lð1 cÞ b r0 b r1 on dbs dbs
ð130Þ
ð131Þ
which can be reduced to b r0
db r1 Bib h 0 ðc; sÞb r1 ¼ lð1 cÞ dbs 1 þ ð1 cÞBib r0
Using Eqs. (121), (125) and (130), we find that
ð132Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
b r1 ¼
a 1 þ ð1 cÞBir0
751
ð133Þ
where a is a constant of integration that will be determined by comparing with the corresponding small-time solution r2 after leading order matching. It can be noted that the outer solution is actually a one-phase solution with a uniform freezing temperature in the liquid. 5.3. Matching So far we have obtained small- and large-times solutions. We would now like to match these solutions in some suitable manner. To do this we propose the existence of an intermediate time, which is sm in the inner time variable and bs m in the outer time variable, at which the inner and outer solutions share common values of interface position and its first time derivative. Thus, our conditions for matching are r ðbs m Þ rðsm Þ ¼ b
ð134Þ
dr db r ðsm Þ ¼ ðbs m Þ ds dbs
ð135Þ
We emphasize that this procedure does not require matching of all derivatives as is commonly done. The leading order balance is r 0 ðbs m Þ r1 ðsm Þ ¼ b
ð136Þ
dr1 db r0 ðsm Þ ¼ ðbs m Þ ds dbs
ð137Þ
which can be written explicitly as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi# " 1 1 1 2ðbs m bs 0 Þ þ þ 2 1 c Bi Bi l " # !9 8 2 1 > > 2 > > > C exp n þ p ðs n m s Þ 1 > > > 1 < = 2 X ¼ Biðsm s Þ þ > > 1 lð1 cÞ > > n¼0 > > p nþ > > : ; 2 and
# " 2 1 X Bi 1 1 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ Bi Cn n þ p ðsm s Þ p exp n þ 2 2 n¼0 2Bi2 ðbs m bs Þ 1þ l
ð138Þ
ð139Þ
These two equations for bs m and bs 0 can be solved numerically. The functional forms of r1 ðsÞ and b r 0 to be tangent to each other at bs ¼ bs m . Fig. 6 is r 0 ðbs Þ guarantee the existence of a bs 0 for r1 and b
752
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762 10
10
10
2
τ c −τ∗
1
0
τ m−τ∗ 10
1
τ o −τ∗ 10
10
2
0
2
4
6
8
10
e ð0; s Þ ¼ 0:286, c ¼ 0:0627, l ¼ 0:582, e ¼ Bi, k ¼ 0:500, / Fig. 6. Variation of sm s , s0 s and sc s with Bi; Bi ¼ 0:0594.
calculated for the leading order terms. Only the first-order term r1 of the inner solution and the zeroth-order term b r 0 of the outer are used for the matching; sc is from Eq. (146). The figure shows sm s , s0 s and solidification completion time sc s (see Section 6.3 for definition of sc ) as functions of Bi. sm is larger than s0 , but both decrease as the heat transfer coefficient is increased. They are also much smaller than the solidification completion time. When considering the matching of the next term in the rðsÞ and b r ðbs Þ expansions, which is 1=2 order higher than the leading order terms r1 ðsÞ and b r 0 ðbs Þ respectively, we find that the constant in Eq. (133) turns out to be a ¼ 0, since the outer solution has a zero Oð3=2 Þ term. From this we know that both the next terms in rðsÞ and b r ðbs Þ expansions are OðÞ higher than their respective r 0 ðbs Þ. This means that the order of accuracy of both rðsÞ and b r ðbs Þ leading order terms r1 ðsÞ and b solutions is actually the same, being OðÞ higher than the leading term. After matching, the perturbation solution of the interface position is summarized as: 8 9 8 > > > " ( ) # > > > 2 > 1 = < X > Cn 1 > 2 > > þ Oð2 Þ; Biðs s Þ þ p ðs s Þ 1 exp n þ > < lð1 cÞ > > 1 2 > > n¼0 ; : nþ p r¼ 2 > ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s " # > > > > 1 1 1 2ðbs bs 0 Þ > > þ OðÞ; bs > bs m > : 1 c Bi þ Bi2 þ l
s 6 sm
ð140Þ
Fig. 7 shows an example of the matching of the small- and large-time solutions of the interface position. Again only the leading order solutions are shown. The forms of the solutions make it possible to always match both the position and its time derivative. Note that b r ðs0 Þ ¼ 0 corresponds to the intercept on the abscissa in this figure.
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
753
0.3
0.2
Outer solution
0.1 Matching point Inner Solution
0
τ=τm 0
το
1
2
3
4
e ð0; s Þ ¼ 0:286, c ¼ 0:0627, l ¼ 0:582, e ¼ Bi ¼ 1:0, k ¼ 0:500, / Fig. 7. Matching of inner and outer solutions; Bi ¼ 0:0594.
Since the matching principles (134) and (135) have been based only on the interface position and its velocity, we would like to check the corresponding matching of the temperature fields. r 0 ðbs m Þ in the expression (121) by r1 ðsm Þ, so that At s ¼ sm , we can replace b b h ðbs m Þ ¼
Bir1 ðsm Þ ½1 þ OðÞð1 nÞ 1 þ ð1 cÞBir1 ðsm Þ
ð141Þ
Expanding, we get b h ðbs m Þ ¼ Bir1 ðsm Þð1 nÞ þ Oð2 Þ
ð142Þ
which is identical to the inner solution hðsm Þ from Eq. (89). So the inner and outer solutions of the temperature distribution in the solid also match at bs m . Furthermore, the temperature distribution in the liquid and its time and spatial derivatives match at sm since the inner solution is exponentially small as seen from Eq. (81), and the outer solution is exactly zero as Eq. (117) shows. The derivatives oh=os, oh=on and o2 h=on2 are also matched at s ¼ sm according to (135) and the linear spatial form of (121) and (89), respectively. The largest temperature variations with time occur at the void–solid interface, and at the adiabatic liquid boundary. Fig. 8 shows the perturbation and numerical solutions of the temperature histories at these locations. /ð1; sÞ is the inner solution /0 ð1; sÞ, and hðc; sÞ is the inner h 0 ðc; sÞ when s > sm . The matching time sm solution h1 ðc; sÞ when s 6 sm and the outer solution b for the two time-scale solution is also indicated. For the history of the solid–liquid interface position, Fig. 9 shows the comparison between the numerical solutions from the enthalpy method and the perturbation solutions. The Bi ¼ 4:1 curve includes terms of the first three orders: r0 , r1 and r3 . For Bi ¼ 1 the inner solution r1 ðsÞ is used r 0 ðsÞ when s > sm . Both the regular perturbation solution of Section 4 when s 6 sm and the outer b
754
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
0.5
0.25
φ(1,τ) 0
θ(γ,τ) –0.25
τ=τ m –0.5
0
2.5
5
7.5
10
12.5
15
Fig. 8. Comparison of numerical and perturbation solutions for temperature history at boundaries: Lines are pere ð0; s Þ ¼ 0:286, e ¼ 4:10, k ¼ 0:500, / turbation solutions and points correspond to numerical solutions; Bi ¼ Bi c ¼ 0:0627, l ¼ 0:582, ¼ 0:0594.
1
0.75
Bi=4.1 Bi=1.0
0.5
0.25
0
0
2
4
6
8
10
12
14
Fig. 9. Comparison of numerical and perturbation solutions for interface positions: Lines are perturbation solutions e ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594. e c ¼ 0:0627, k ¼ 0:500, / and points correspond to numerical solutions; Bi ¼ Bi Bi,
(Bi ¼ 1=2 ¼ 4:1 is used) and the matched two time-scale solution of Section 5 (Bi ¼ 1 is used) agree well with numerical results within the expressed accuracy of the perturbation expansion.
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
755
6. Discussion 6.1. Extension of Bi range The two solutions of Sections 4 and 5 can be extended to cover the entire range of possible Bi. The Bi Oð1=2 Þ solution in Section 4 is also valid for Bi > Oð1=2 Þ. As a special case it is valid as Bi ! 1, which corresponds to isothermal cooling. In fact it can be checked that in this e ! 1, limit the solution calculated here tends to that given by Weinbaum and Jiji [7] if we take Bi Bi ! 1 and c ¼ 0. We can also consider the range of Bi between Oð1Þ and Oð1=2 Þ. The interface position computed from the Oð1=2 Þ regular perturbation or the Oð1Þ two time-scale solution are numerically similar, as shown in Fig. 10 with comparison to the numerical solution, so that either one may be used. In the figure, the regular perturbation solution includes terms of the first three orders: r0 , r1 and r2 . The matched asymptotic solution is the inner solution r1 ðsÞ when s 6 sm and the outer solution b r 0 ðsÞ when s > sm . For Bi < Oð1Þ, which corresponds to very weak cooling, the solution of Section 5 for Bi Oð1Þ remains valid. In obtaining the inner solution, the leading order balance for the boundary condition (37) could include OðÞ terms. 6.2. One-phase case Another interesting situation is when Ti ¼ Tf where the liquid is always at its freezing temperature. In this case, s ¼ Cn ¼ 0. To compare with the results of Yan and Huang [6], we take Bi Oð1Þ and c ¼ 0. The matching condition gives sm ¼ s0 ¼ 0 and the interface position is 1
0.75 Matched Asymptotic Solution
0.5
Regular Perturbation Solution
Numerical Solution
0.25
0
0
2
4
6
8
10
Fig. 10. Comparison between regular and two time-scale perturbation solutions in Oð1Þ < Bi < Oð1=2 Þ range: Lines e ð0; s Þ ¼ 0:286, e ¼ 2:0, k ¼ 0:500, / are perturbation solutions and points correspond to numerical solutions; Bi ¼ Bi c ¼ 0:0627, l ¼ 0:582, ¼ 0:0594.
756
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
1 r¼ þ Bi
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2bs þ þ OðÞ 2 Bi lð1 cÞ
ð143Þ
To the leading order this limit is identical to the solution in [6]. 6.3. Interface motion and solidification time The solutions in Sections 4 and 5 both show that the solid–liquid interface moves with an initial velocity of " # 1 X dr 1 ðs Þ ¼ Bi ð144Þ nþ pCn þ h:o:t ds lð1 cÞ 2 n¼0 where h.o.t are the higher order terms. The starting velocity is identical for all Bi, though to different orders of accuracy. From a physical point of view, this is more realistic than the infinite starting velocity predicted by the models of Wilson and Solomon [5] and Weinbaum and Jiji [7]. The subsequent interface motion after starting is much slower. Fig. 11 shows the displacement of the solidification front for different Bi. The Bi ¼ 20, 10 and 4.1 curves include terms of the first three orders: r0 , r1 and r2 . For Bi ¼ 1 and 0.5 the inner solution r1 ðsÞ is used when s 6 sm and the outer b r 0 ðsÞ when s > sm . As is intuitively obvious, the front is seen to move faster as the cooling is increased. The total time sc necessary to solidify the whole slab is of special interest in this problem. Taking r ¼ 1 in Eq. (140) for Bi Oð1Þ we find that " # 2 lð1 cÞ lð1 cÞ þ ð145Þ sc ¼ 1 þ bs 0 þ Oð1Þ Bi 2
1 Bi=20 10
4.1 1.0
0.75
0.5
0.5
0.25
0
0
5
10
15
20
25
e ð0; s Þ ¼ 0:286, c ¼ 0:0627, e k ¼ 0:500, / Fig. 11. Solid–liquid interface position for various Biot numbers, Bi ¼ Bi Bi, l ¼ 0:582, ¼ 0:0594.
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
757
where bs 0 is determined from the matching. Similarly, from the Bi Oð1=2 Þ solution, we have sc ¼
1
lð1 cÞ2 lð1 cÞ þ 1=2 D þ s þ OðÞ þ 1=2 Bi0 2
ð146Þ
where 2
sffiffiffi " " 2 # rffiffiffi# 3 n þ 12 p2 n þ 12 p l 2 6 7 erfc Bi0 exp pffiffiffi 1 6 7 Bi l 2 C 0 l X 2 6 7 Cn 61 þ D¼ ð1 cÞ 7 2 7 Bi0 n¼0 6 2 1 3=2 4 5 p nþ 2
ð147Þ
The leading orders of sc in Eqs. (145) and (146) are identical. For small , the solidification time may be explained in the following manner. We take a linear temperature profile in the solid and constant one in the liquid. The heat balance equation then gives Tf T0 ds ¼ ql L 1 ð1 cÞs dt þ h ks Integrating this equation and applying the conditions sð0Þ ¼ 0 and sðtc Þ ¼ ‘, we get ql L‘ ð1 cÞ‘L tc ¼ 1þ ðTf T0 Þh 2ks which in nondimensional terms is lð1 cÞ ð1 cÞBi 1þ sc ¼ Bi 2
ð148Þ
ð149Þ
ð150Þ
These are the first two terms of Eqs. (145) and (146). For large Bi, the limitation to interface motion is not the thermal resistance at the void–solid e ! 1, Bi ! 1, and interface but the conduction resistance within the solid. As Bi 2 sc ! lð1 cÞ =2 so that it becomes independent of Bi. The variation of sc with Bi is also graphed in Fig. 6. It may be noted that since sc is much larger than sm , the outer solution of Section 5 is the one that is applicable for most of the time. Figs. 12 and 13 show the effect of the shrinkage factor, c, on the total time taken to solidify the slab, sc . In Fig. 12 the first three orders, r0 , r1 and r2 , are used while for Fig. 13 the inner solution r 0 ðsÞ for s > sm . We see that the larger c is, the r1 ðsÞ is used for s 6 sm and the outer solution b smaller is sc . This is also indicated in the first terms of Eqs. (145) and (146). Mathematically this comes from the fact that increasing shrinkage reduces the thermal resistance of the solidified layer, so the same cooling is becoming more efficient as the solid–liquid interface advances. This incorrect prediction is due to the physically unrealistic assumption of a Biot number that is constant (as assumed by Reddy and Beckermann [8]) in spite of the growth of the void, a restriction that is lifted in the following discussion.
758
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762 1 γ=
0.2
0.1
7
0.0
62
0.0
0.75
0.5
0.25
0
0
2
4
6
8
e ¼ 4:10, Fig. 12. Shrinkage effect on solid–liquid interface position from regular perturbation solution with Bi ¼ Bi e k ¼ 0:500, / ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594. 1 γ=
0.2
0.1 7 62 0.0
0.0
0.75
0.5
0.25
0
0
4
8
12
16
e ¼ Fig. 13. Shrinkage effect on solid–liquid interface position from two time-scale perturbation solution with Bi ¼ Bi e ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594. 1:0, k ¼ 0:500, /
6.4. Time-dependent Bi So far we have assumed that Bi is a constant. However, it can be argued on physical grounds that the heat transfer coefficient must change as the void between the solid surface and the mold wall grows. The proper change will only occur if the thermal resistance of the void is not small
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
759
compared to that of the solid. The solutions derived here can also be applied to time-dependent Bi e as long as they vary slowly compared to the solidification process. and Bi Bi, As a example, let us take the heat transfer coefficient to be inversely proportional to the size of the void to compare with the results of Wilson and Solomon [5] who obtain a similarity solution for this case. Thus if we take h ¼ kv =csðtÞ, where kv is the thermal conductivity in the void, we have BiðsÞ ¼
kv crðsÞ
ð151Þ
where kv ¼ kv =ks . We note that at s ¼ 0, Bi ! 1 and thus we have s ¼ 0. Eq. (59) for large Bi reduces to vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2as t u ð152Þ s¼u t 2 1 ð1 cÞ 1 þ c kv in dimensional terms. This form is similar to that given in [5], and turns out to be the leading order of the solution obtained in [5] if the transcendental equation for the interface position given in that work is expanded in power series of 1=2 . The quantitative difference of the interface position between this solution and those in the numerical example presented there is within 2%. The model represented by Eq. (151) has the obvious drawback of predicting an infinite heat transfer coefficient at the beginning of solidification. One can improve it by adding a contact thermal resistance to remove this singularity. If heat transfer across the void is by conduction alone and natural convection and radiation are negligible, as is possible to justify in some applications of practical interest [9], the Biot number takes the form e Bi
ð153Þ e Bi 1 þ crðsÞ kv e due to small contact resistance. The above form of the time-dependent with a possibly large Bi Biot number is substituted for Bi0 into the regular perturbation solution (59) to give sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 2s þ þ ð154Þ r1 ¼ 2 2 e e b blð1 cÞ Bib Bi Bi BiðsÞ ¼
where
2 1 b¼1þc kv
ð155Þ
Fig. 14 shows the motion of the interface for different values of c calculated from r1 above. They agree reasonably well with the numerical solutions. As c increases, it takes longer for the slab to solidify and sc increases, as long as kv < 2 according to the above expression. This can be contrasted with the behavior shown in Fig. 8 which was for the constant Bi (i.e. for kv ! 1). In practice [9], kv 1 (for our case it is 2.73 · 104 ); the interface motion is very sensitive to kv , as shown in Fig. 15 where r1 is calculated from Eq. (154).
760
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762 1 0.75 0.5
γ=0.0627
0.25 0
0
2
4
1
τ−τ∗
6
8
10
6
8
10
0.75 0.5 γ=0.20
0.25 0
0
2
4
Fig. 14. Shrinkage effect on solid–liquid interface position from time varying Bi model: Lines are numerical solutions e ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594, kv ¼ 0:50. e ¼ 4:10, k ¼ 0:50, / and points correspond to analytical expressions; Bi 1 γ=0.0627 0.1 0.75
0.2
0.5
0.25
0
0
1000
2000
3000
4000
5000
6000
e ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594, kv ¼ 2:73 104 . e ¼ 4:10, k ¼ 0:50, / Fig. 15. Solid–liquid interface position, Bi
7. Conclusions The one-dimensional solidification problem in a finite domain has been solved using perturbation methods with the Stefan number as the small parameter. Initially, the liquid is at a uniform
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
761
temperature above the freezing temperature. The driving mechanism for solidification is cooling at one boundary, while the other end of the slab is adiabatic. The void that is formed due to density change during solidification is assumed to appear between the mold wall and the cold end of the solid. The two-phase problem with convective boundary condition is more complicated to handle than either the isothermal condition analyzed in [7], or the one-phase problem studied in [6]. The convective mechanism gives rise to a nondimensional parameter, the Biot number Bi, in the governing equations. The nature of the solution depends on the order of magnitude of Bi in terms of the perturbation parameter, the Stefan number . There is a precooling stage after which solidification begins. The temperature distribution in the precooling stage is obtained as a series solution. The perturbation method is used once solidification actually begins. The range of perturbation solutions of the interface position and temperature distributions in both solid and liquid that are obtained for different magnitudes of cooling are the following: 1. Oð1=2 Þ 6 Bi < 1: This corresponds to strong cooling for which a regular perturbation solution is obtained. 2. Oð1Þ < Bi < Oð1=2 Þ: This corresponds to intermediate cooling. The solid–liquid interface position computed from the solutions in cases (1) above or (3) below are not much different. Either can be used for numerical purposes. 3. 0 < Bi 6 Oð1Þ: This corresponds to weak and moderate cooling. A two time-scale solution describes the interface motion and temperature distribution. The inner solution is valid up to an intermediate time sm . From sm onwards, the outer solution takes over. The solutions approach those of Refs. [5–7] in the appropriate limits. The analysis is also validated using numerical solutions obtained by the enthalpy method. The results obtained here are valid not only for constant Biot numbers, but can also be applied to time-dependent parameters as long as they are slowly varying. This greatly extends the usefulness of the present results and enables them to be tailored to any void thermal resistance or wall temperature model.
Acknowledgement We gratefully acknowledge support of this research by A.E. Goetze of South Bend, Indiana.
References [1] [2] [3] [4] [5]
H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Oxford University Press, London, 1959. J. Crank, Free and Moving Boundary Problems, Oxford University Press, New York, 1984. L.I. Rubinstein, The Stefan problem, Trans. Math. Monographs 27 (1971) 94–181. M.A. Biot, Variational Principles in Heat Transfer, Clarendon Press, Oxford, 1970. D.G. Wilson, A.D. Solomon, A Stefan-type problem with void formation and its explicit solution, IMA J. Appl. Math. 37 (1986) 67–76.
762
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
[6] M.M. Yan, P.N.S. Huang, Perturbation solutions to phase change problem subjected to convection and radiation, ASME J. Heat Transfer 101 (1979) 96–100. [7] S. Weinbaum, L.M. Jiji, Singular perturbation theory for melting or freezing in finite domains initially not at fusion temperature, ASME J. Appl. Mech. 44 (1977) 25–30. [8] A.V. Reddy, C. Beckermann, Measurements of metal-mold interfacial coefficients during solidification of Sn and Sn–Pb alloys, Exper. Heat Transfer 6 (1993) 111–129. [9] Z. Yang, Analytical and numerical study of solidification by convective cooling with void formation, M.S. Thesis, University of Notre Dame, Notre Dame, Indiana, 1994. [10] V. Alexiades, A.D. Solomon, Mathematical Modeling of Melting and Freezing Process, Hemisphere Publ. Corp., Washington, DC, 1993.
Solidification of a finite slab with convective cooling and shrinkage Zhaoji Yang 1, Mihir Sen *, Samuel Paolucci Department of Aerospace and Mechanical Engineering, University of Notre Dame, Notre Dame, IN 46556-5637, USA Received 29 April 2002; received in revised form 12 March 2003; accepted 28 March 2003
Abstract We consider the one-dimensional solidification of a finite slab of a pure substance which is initially in liquid state. Convective cooling is applied at one end of the slab while the other end is kept adiabatic. There is a precooling stage before solidification begins. A void, which has substantial overall effect on the solidification of the slab, is formed at the cooled end due to shrinkage and grows during solidification. Perturbation analysis for small Stefan numbers is used to find the interface position and temperature distributions in both solid and liquid phases. The order of magnitude of the Biot number at the cold end, Bi, determines the nature of the solution. Regular and two time-scale perturbation solutions are obtained for Bi Oð1=2 Þ and Bi Oð1Þ, respectively. The results can be extended to other orders of magnitude, as well as to slowly-varying time-dependent Bi. The results agree well with special cases in the literature, and with numerical solutions calculated using the enthalpy method. 2003 Elsevier Inc. All rights reserved.
1. Introduction Heat transfer problems involving phase change such as solidification and melting are of great importance in engineering applications [1]. The mathematical model is a system of partial differential equations of the parabolic type with moving boundaries. Difficulties are encountered in the analytical solution of these problems due to the nonlinearity of moving boundary conditions.
*
Corresponding author. Tel.: +1-574-631-5975; fax: +1-574-631-8341. E-mail address: [email protected] (M. Sen). 1 Currently at Advanced Technology Development, Ford Visteon, Electronics Technical Center, 17000 Rotunda Dr., Dearborn, MI 48121-6010. 0307-904X/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0307-904X(03)00078-7
734
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
Nomenclature Bi c h k ‘ L s t T Tf T0 Ti x Greek a b c f g h H k kv l n q r s / W
¼ hl=ks , Biot number specific heat heat transfer coefficient thermal conductivity thickness of slab latent heat of fusion solid–liquid interface position time temperature freezing temperature mold temperature initial temperature coordinate along slab symbols ¼ k=qc, thermal diffusivity ¼ 1 cð2=kv 1Þ ¼ 1 ql =qs , shrinkage factor Stefan number nondimensional coordinate for numerical method nondimensional coordinate in liquid nondimensional temperature in solid nondimensional temperature for numerical method ¼ kl =ks ¼ kv =ks ¼ al =as nondimensional coordinate in solid density ¼ s=‘, nondimensional interface position nondimensional time nondimensional temperature in liquid nondimensional enthalpy
Subscripts and superscripts l liquid s solid m matching v void outer variable precooling stage * beginning of solidification
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
735
Exact solutions are possible in similarity form only in limited cases. Many approximate and numerical methods of solution for nonsimilar cases have also been developed [2–4]. The effect of shrinkage and void formation due to solidification was studied by Wilson and Solomon [5] who obtained a similarity solution by assuming the void thermal resistance proportional to its size. Perturbation methods based on the assumption of a small Stefan number have also been used to deal with finite domain solidification problems. There are two types of problems that have been considered: one where the liquid is initially superheated above the freezing temperature and in which the temperature in both solid and liquid must be determined (two-phase), and the other where the liquid is not superheated but is always at the freezing temperature (one-phase). A regular perturbation analysis of the one-phase problem for convective and radiative cooling has been carried out by Yan and Huang [6]. The radiative boundary condition which they presented was simply linearized to convective form. They could use the interface position as an independent variable instead of time since the interface position increases monotonically with respect to time. In their conclusions, they foresaw some difficulty in applying the perturbation technique to the convective solidification of an initially superheated slab. However, their analysis did not include the effect of shrinkage. Significant progress for the two-phase problem in a finite domain was reported by Weinbaum and Jiji [7]. Cooling was by an isothermal wall below freezing temperature while the other end was either adiabatic or isothermal. Shrinkage was not included here either. The problem was solved using the method of matched asymptotics. There has been no analytical treatment, however, of the finite domain solidification problem of a two-phase region due to convective cooling. This kind of boundary condition is fairly general and includes the isothermal as a limiting case. The presence of the shrinkage void is an important factor in the solidification process. We will treat this problem using perturbation methods. The solutions thus obtained will be compared to previous results under appropriate limits and also verified by numerical computation. The geometry is taken to be one-dimensional as shown in Fig. 1. The molten material at temperature Ti is initially between x ¼ 0 and x ¼ ‘. Cooling is then applied at x ¼ 0, which could represent the surface of the mold where solidification is taking place. In the precooling stage, the heat transfer coefficient is denoted as h~, which includes effects of thermal contact between the liquid and the mold, and then convective cooling to the mold at temperature T0 . A constant value of h~ is assumed throughout the precooling stage. Eventually the temperature of the liquid at x ¼ 0 x=0
x=γs(t)
x=s(t)
x=l
Convective Cooling
Void
Liquid
Solid motion
T=T0
motion
T=Tf
x Fig. 1. Schematic of the physical problem.
Adiabatic
736
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
comes down to the freezing point and solidification begins. As the material solidifies, it shrinks and a void is formed between the solidifying material and the mold surface, as seen in experiments [8]. This is the region between x ¼ 0 and x ¼ csðtÞ, where c ¼ ð1 ql =qs Þ is the shrinkage parameter quantifying the effect of density change. It must be mentioned that in reality void formation often occurs in two different ways: through micropores within the material and what is referred to as a shrinkage gap. The proportion of one to the other depends greatly on external conditions. We assume here that it is entirely of the second kind. Note also that while we are considering a one-dimensional analysis, the result should be locally applicable to higher dimensional problems as long as the characteristic length scales in the transverse directions are much larger than c‘. We assume that the heat transfer across this void and the mold to the ambient can be represented by a heat transfer coefficient h. We will take h to be a constant [8], the same as h~. Moreover, the analysis in which it is varying slowly with time is considered in Section 6.4. The outer boundary of the liquid at x ¼ ‘ is taken to be adiabatic. In what follows, wherever numerical values are needed we will consider the solidification of pure aluminum.
2. Temperature distribution before solidification In a one-phase cooling problem, such as in [7], the material would begin to solidify as soon as cooling is applied because of the instantaneous drop of the temperature at x ¼ 0. In the present problem, however, there is a period of precooling before the solidification can actually start. The material is all in liquid form for time t < t , where t is the time necessary to bring the liquid at the boundary to the freezing temperature. Before looking at solidification, we would like to determine the temperature distribution in the liquid at t ¼ t . In the precooling stage, heat transfer in the liquid is by single-phase conduction. The governing equation and boundary conditions for the temperature distribution, Tl ðx; tÞ, are al
o2 Tl oTl ; ¼ ox2 ot
kl
oTl ð0; tÞ ¼ h~½Tl ð0; tÞ T0 ox
06x6‘
ð1Þ ð2Þ
oTl ð‘; tÞ ¼ 0 ox
ð3Þ
Tl ðx; 0Þ ¼ Ti
ð4Þ
These can be nondimensionalized to e e o/ o2 / ; ¼ os og2
06g61
e e Bi o/ e ð0; sÞ ð0; sÞ ¼ / k og
ð5Þ
ð6Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
737
e o/ ð1; sÞ ¼ 0 og
ð7Þ
e ðg; 0Þ ¼ 1 /
ð8Þ
where the nondimensional liquid temperature, distance and time are given by e ¼ Tl T0 / Ti T0 g¼
x ‘
ð9Þ ð10Þ
al t ð11Þ ‘2 respectively, and k ¼ kl =ks . The controlling parameter is the heat transfer coefficient which is nondimensionalized into a Biot number defined by s¼
~ e ¼ h‘ Bi ks
ð12Þ
It is easy to find a series solution to this conduction problem [1]; the resulting temperature distribution is e ðg; sÞ ¼ /
1 X
An expðpn2 sÞ cos½pn ð1 gÞ
ð13Þ
n¼0
where An ¼
2 sin pn # e k Bi pn 1 þ 2 e þ k2 p2 Bi n "
ð14Þ
and pn is the nth positive eigenvalue determined from the transcendental equation tan pn ¼
e Bi kpn
ð15Þ
The time s ¼ s when the cold end drops to the freezing point is determined implicitly from e ð0; s Þ ¼ Tf T0 / Ti T0 e ð0; s Þ is a parameter about initial preheating that must be specified. Note that / The temperature distribution in the liquid at the instant solidification begins is e ðg; s Þ ¼ /
1 X n¼0
An expðpn2 s Þ cos½pn ð1 gÞ
ð16Þ
ð17Þ
738
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
e it becomes flatter as Bi e decreases. Fig. 3 shows the This is shown in Fig. 2 for different Bi Bi; e We have assumed an aluminum slab with ‘ ¼ 0:577 m, k ¼ 0:500 and Bi. variation of s with Bi freezing temperature Tf ¼ 660 C. Other temperatures are Ti ¼ 710 C and T0 ¼ 640 C, thus e ð0; s Þ ¼ 0:286. Typically, a precooling heat transfer coefficient h~ ¼ 364 W/m2 K corresponds to / e ¼ 1:0. a Biot number Bi 1 Bi = 10 6 4.1 0.8
3 2
0.6 1
0.4
φ∗=0.286 0.2
0
0.2
0.4
0.6
0.8
1
e ð0; s Þ ¼ 0:286. Fig. 2. Temperature distribution in the liquid at beginning of solidification; k ¼ 0:500, /
10
10
1
0
φ∗=0.5
–1
10
φ∗=0.286
φ∗=0.1
–2
10
–3
10
0
2
4
6
8
e k ¼ 0:500. Fig. 3. Variation of preheating time s with Bi; Bi
10
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
739
3. Governing equations during solidification Once solidification begins we have to deal with both liquid and solid phases. Due to shrinkage there is a void formed which, consistent with some experiments [8], we assume to be near the cold end. The growth of the void leads to bulk motion of the solid phase, adding a convective term to the heat conduction equation in the solid. It is assumed, however, that there is no convective motion in the liquid [5]. The governing equations are thus as
o2 Ts oTs ds oTs þc ; ¼ 2 dt ox ox ot
al
o2 Tl oTl ¼ ; ox2 ot
t P t ;
t P t ; s6x6‘
cs 6 x 6 s
ð18Þ ð19Þ
for each of the two phases separated by the interface at x ¼ sðtÞ. The boundary and initial conditions are ks
oTs ðcs; tÞ ¼ h½Ts ðcs; tÞ T0 ox
ð20Þ
Ts ðs; tÞ ¼ Tf
ð21Þ
Tl ðs; tÞ ¼ Tf
ð22Þ
ks
oTs oTl ds ðs; tÞ kl ðs; tÞ ¼ ql L dt ox ox
ð23Þ
oTl ð‘; tÞ ¼ 0 ox
ð24Þ
sðt Þ ¼ 0
ð25Þ
e ðx; t Þ Tl ðx; t Þ ¼ T0 þ ðTi T0 Þ /
ð26Þ
e ðx; t Þ is calculated from Eq. (17). We introduce the liquid and solid nondimensional where / temperature variables /¼k h¼
Tl Tf Tf T0
Ts Tf Tf T0
ð27Þ ð28Þ
A Landau transformation is used to convert the moving boundary problem into one of fixed domain. The new spatial variables in the solid and liquid are n¼
x s
ð29Þ
740
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
xs ð30Þ ‘s where s is defined as before. The nondimensional position of the solid–liquid interface is r ¼ s=‘, and the thermal diffusivity ratio is l ¼ al =as . The governing equations in these variables are o2 h dr oh 2 oh ¼ l r ð31Þ ðn cÞr ; s P s ; c 6 n 6 1 2 os ds on on g¼
o2 / o/ dr o/ ¼ ð1 rÞ2 ð1 gÞð1 rÞ ; 2 og os ds og
s P s ; 0 6 g 6 1
ð32Þ
with boundary conditions hð1; sÞ ¼ 0
ð33Þ
/ð0; sÞ ¼ 0
ð34Þ
o/ ð1; sÞ ¼ 0 og oh o/ dr ð1 rÞ ð1; sÞ r ð0; sÞ ¼ rð1 rÞ lð1 cÞ on og ds
ð35Þ ð36Þ
oh ðc; sÞ ¼ rBi½hðc; sÞ þ 1 on
ð37Þ
rðs Þ ¼ 0
ð38Þ "
e ðg; s Þ / 1 /ðg; s Þ ¼ k e ð0; s Þ /
#
ð39Þ
The nondimensional parameters controlling the solidification process are ¼ cs Bi ¼
Tf T0 L
h‘ ks
ð40Þ ð41Þ
where is the Stefan number and Bi is the Biot number. The Stefan number is the ratio of sensible heat to latent heat. The Biot number represents the ratio between the conductive thermal resistance in the solid and the thermal resistance across the void to the mold. The thermal diffusion time scale in the liquid is tdiff ¼ ‘2 =al (the diffusion time in the solid is of the same order), and the void heat transfer time scale is tvoid ¼ ‘ql cl =h. They are related by ktdiff =tvoid ¼ Bi. In this paper we have chosen to nondimensionalize all times by tdiff . Note that the solidification time scale is a dependent variable given by tdiff =Bi. In addition, it is possible that there may be time scales associated with changes in e h , h and T0 . The solutions developed here will remain valid as long as they are slowly varying; more will be said about this in Section 6.4.
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
741
Bi is proportional to h, which is the heat transfer coefficient after the separation of the solidified e on the other hand is its value before solidification begins, i.e. material from the mold wall. Bi Bi, e Generally, we can before the void is formed. We will use the same constant value for Bi as for Bi Bi. e expect that BiðsÞ 6 Bi because of the increasing insulating effect of the void. The transition should e at the instant solidification starts (see the model presented in be continuous with Biðs Þ ¼ Bi Section 6.4). When Bi ! 1, the boundary condition (37) reduces to the isothermal condition analyzed by Weinbaum and Jiji [7]. Here we will investigate the problem for finite Bi. We will also assume that l, k and ð1 cÞ are all fixed by material properties and are of Oð1Þ. The Stefan number only appears in the interface condition (36). It is small in many practical applications, and we will use it as a small parameter to develop a perturbation solution. The analysis will depend on the order of magnitude of the Biot number. We will consider two different cases, Bi Oð1=2 Þ and Bi Oð1Þ, with which we can cover the entire range of Bi. The solidification problem can also be solved numerically to check the validity of the perturbation solutions. We use an enthalpy method in which the solid and liquid regions are treated as a continuum, except for an enthalpy jump at the interface [10]. To implement the numerical method, we introduce f¼
x csðtÞ ; lð1 cÞ
2
W ¼ lð1 cÞ H;
H¼
1 Ts Tf kð1 cÞ Tf T0
ð42Þ
in the solid and x f¼ ; l
W¼Hþ
lð1 cÞ ; k
H¼
Tl Tf Tf T0
ð43Þ
in the liquid. The equations are nondimensional; the convective term has also been removed. The original system of Eqs. (18) and (19) and boundary conditions (20)–(26) reduce to o2 H oW ¼ ;06f61 os of2
ð44Þ
oH ð1; sÞ ¼ 0 of
ð45Þ
oH Bi ð0; sÞ ¼ ð1 cÞBiHð0; sÞ þ of k
ð46Þ
e ðf; s Þ / 1 e ð0; s Þ /
ð47Þ
Hðf; s Þ ¼
The problem is solved using a finite difference method with second-order central difference in space and explicit Euler marching in time [9]. The numerical values for aluminum are c ¼ 0:0627, k ¼ 0:50, and l ¼ 0:582. For the temperatures assumed before, we have ¼ 0:0594.
742
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
4. Solution for Bi O (1=2 ) Let us first consider a large heat transfer coefficient corresponding to Bi Oð1=2 Þ. The diffusion time is large compared to the void heat transfer time, and the solidification process is diffusion limited. We seek perturbation solutions of the form hðn; sÞ ¼
1 X
n=2 hn ðn; sÞ
ð48Þ
n¼0
/ðg; sÞ ¼
1 X
n=2 /n ðg; sÞ
ð49Þ
n¼0
rðsÞ ¼
1 X
n=2 rn ðsÞ
ð50Þ
n¼0
We follow the standard procedure of substituting these expansions into the governing equations (31) and (32) and boundary and initial conditions (33)–(39), and balancing terms in each order. We also introduce the Oð1Þ parameter Bi0 Bi 1=2 . From the Oð1Þ terms in the interface condition (36), we have r0 ¼ 0
ð51Þ
Collecting terms of Oð1Þ in Eqs. (31), (33) and (37), we get o2 h0 ¼0 on2
ð52Þ
h0 ð1; sÞ ¼ 0
ð53Þ
oh0 ðc; sÞ ¼ r1 Bi0 ½h0 ðc; sÞ þ 1 on
ð54Þ
which gives us h0 ðn; sÞ ¼
r1 Bi0 ðn 1Þ r1 Bi0 ð1 cÞ þ 1
ð55Þ
Using this in the OðÞ balance of Eq. (36), we obtain r1 ðsÞ
dr1 1 oh0 ¼ ð1; sÞ lð1 cÞ on ds
ð56Þ
From this dr1 1 Bi0 ¼ lð1 cÞ 1 þ ð1 cÞBi0 r1 ds
ð57Þ
r1 ðs Þ ¼ 0
ð58Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
which can be solved as 1 r1 ¼ 1c
1 þ Bi0
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 s s þ 2 l Bi20
743
ð59Þ
From a balance of Oð1Þ terms in Eq. (32), and boundary conditions (34) and (35), we can write o2 /0 o/0 ¼ og2 os
ð60Þ
/0 ð0; sÞ ¼ 0
ð61Þ
o/0 ð1; sÞ ¼ 0 og "
ð62Þ
# e ðg; s Þ / 1 /0 ðg; s Þ ¼ k e ð0; s Þ /
ð63Þ
e ðg; s Þ as Oð1Þ contribution in the Here we include the initial liquid temperature distribution / perturbation expansion. By separation of variables, we obtain the solution # " 2 1 X 1 1 Cn sin n þ p2 ðs s Þ ð64Þ pg exp n þ /0 ðg; sÞ ¼ 2 2 n¼0 where Cn ¼
1 X 2k þ n þ 12 p m¼0
1 nþ p cos pm 2 " # 2 1 e ð0; s Þ n þ p2 pm2 / 2
2Am k expðpm2 s Þ
ð65Þ
Balancing the Oð1=2 Þ terms in Eqs. (31), (33) and (36), we find o2 h1 ¼0 on2
ð66Þ
h1 ð1; sÞ ¼ 0
ð67Þ
oh1 ðc; sÞ ¼ r2 Bi0 ½h0 ðc; sÞ þ 1 þ Bi0 r1 ðsÞh1 ðc; sÞ on
ð68Þ
We solve for h1 to get Bi0 r2 ½h0 ðc; sÞ þ 1 ðn 1Þ h1 ðn; sÞ ¼ 1 þ ð1 cÞBi0 r1 ðsÞ From a balance of the Oð3=2 Þ terms in the interface condition (36), we have 1 oh1 oh0 o/ dðr1 r2 Þ dr1 r21 ð1; sÞ r1 ð1; sÞ r1 0 ð0; sÞ ¼ lð1 cÞ on ds on og ds
ð69Þ
ð70Þ
744
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
Substituting previously obtained expressions for r1 ðsÞ, h0 ðn; sÞ, h1 ðn; sÞ, and /0 ðg; sÞ into this, we get dr2 ¼ AðsÞr2 þ BðsÞ ds
ð71Þ
r2 ðs Þ ¼ 0
ð72Þ
where AðsÞ ¼
Bi20 l þ 2Bi20 ðs s Þ
ð73Þ
# " 2 1 X 1 1 1 Cn n þ p2 ðs s Þ BðsÞ ¼ p exp n þ lð1 cÞ n¼0 2 2
ð74Þ
The general solution is of the form r2 ðsÞ ¼ exp
Z
s
Aðs0 Þ ds0
s
Z
s
" Bðs00 Þ exp
s
Z
s00
# Aðs0 Þ ds0 ds00
ð75Þ
s
which can be evaluated as " # Z s sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 2 2 0 1 Cn n þ p2 ðs0 s Þ ds0 1 þ Bi0 ðs s Þ exp n þ p 2 l 2 s n¼0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r2 ðsÞ ¼ 2 lð1 cÞ 1 þ Bi20 ðs s Þ l " # 8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 1 > > 2 > > 1 < 1 1 þ l Bi0 ðs s Þ exp n þ 2 p ðs s Þ X 1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Cn ¼ > 1 2 2 > n¼0 > nþ p lð1 cÞ 1 þ Bi0 ðs s Þ > : 2 l " 9 2 # rffiffiffi l n þ 12 p2 2 > > Bi0 exp > > rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi = 1 2Bi20 l l n þ p 1 l 2 p ffiffi ffi p ðs s Þ þ erf n þ erf þ 2 > 2 2Bi20 2Bi0 1 > > > p3=2 nþ ; 2 1 X
ð76Þ
The perturbation solution for rðsÞ has been obtained up to OðÞ. Higher order solutions are difficult to obtain in explicit analytical form because of the nonlinearity of the corresponding differential equations. Even so, we can obtain good accuracy with the first few perturbation terms as long as the parameter is small enough. That the perturbation series converges fast can be seen
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
745
4
σ1
3
1/2
σ+ε σ2 1 2
−σ 2
1
0
0
2
4
6
e ¼ Fig. 4. Comparison of first- and second-order Bi Oð1=2 Þ perturbation solutions for interface location; Bi ¼ Bi e ð0; s Þ ¼ 0:286, c ¼ 0:0627, l ¼ 0:582, ¼ 0:0594. 4:10, k ¼ 0:500, /
from comparing the terms r1 and r2 in Fig. 4. The second-order solution is only important near the beginning of solidification. One can easily verify that since r1 and r2 are Oð1Þ for Bi Oð1=2 Þ, the expansion remains uniformly valid. Fig. 5 shows temperature histories at the the two ends of the slab for Bi ¼ 1=2 . One term is used in Eq. (48) for hðc; sÞ and two in (49) for /ð1; sÞ. The perturbation solution is seen to compare well with the numerical solution. 5. Solution for Bi O(1) This corresponds to a moderate heat transfer coefficient h. Both diffusion and void time scales are of the same order. 5.1. Short-time solution Following the previous form of the perturbation expansions (48)–(50), we find that r2nþ1 ¼ 0. Actually, in both the perturbation expansion of nondimensionalized temperature and solid–liquid interface position function, the odd n terms drop out. Thus we can simply use the expansions 1 X n hn ðn; sÞ ð77Þ hðn; sÞ ¼ n¼0
/ðg; sÞ ¼
1 X n¼0
n /n ðg; sÞ
ð78Þ
746
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762 1.5
1
0.5
φ(1,τ) 0
–0.5
θ(γ,τ) –1
–1.5
0
1
2
3
4
5
6
7
8
Fig. 5. Comparison of numerical and Bi Oð1=2 Þ perturbation solutions for temperature history at boundaries; solid e ¼ 4:10, lines are from perturbation solution, and dashed lines with points correspond to numerical solution; Bi ¼ Bi e ð0; s Þ ¼ 0:286, c ¼ 0:0627, l ¼ 0:582, ¼ 0:0594. k ¼ 0:500, /
rðsÞ ¼
1 X
n rn ðsÞ
ð79Þ
n¼0
On gathering the terms of order Oð1Þ, we find that r0 ¼ 0 /0 ¼
1 X
ð80Þ Cn sin
n¼0
# " 2 1 1 nþ p2 ðs s Þ pg exp n þ 2 2
ð81Þ
and o2 h0 ¼0 on2
ð82Þ
h0 ð1; sÞ ¼ 0
ð83Þ
oh0 ðc; sÞ ¼ 0 on
ð84Þ
which means that h0 ðn; sÞ ¼ 0 The balance of OðÞ terms gives o2 h1 ¼0 on2
ð85Þ
ð86Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
747
h1 ð1; sÞ ¼ 0
ð87Þ
oh1 ðc; sÞ ¼ r1 Bi on Solving, we get
ð88Þ
h1 ðn; sÞ ¼ r1 Biðn 1Þ
ð89Þ
2
From the balance of Oð Þ terms in Eq. (36), we obtain dr1 1 oh1 o/ ¼ ð1; sÞ r1 0 ð0; sÞ r1 lð1 cÞ on ds og
ð90Þ
which gives
Thus
" )# ( 2 1 X dr1 1 1 1 Bi ¼ Cn n þ p2 ðs s Þ p exp n þ lð1 cÞ 2 2 ds n¼0
ð91Þ
r1 ðs Þ ¼ 0
ð92Þ
( " 1 X 1 C n r1 ¼ exp Biðs s Þ þ 1 lð1 cÞ n þ p n¼0 2
1 nþ 2
2
! p2 ðs s Þ
)# 1
ð93Þ
The Oð2 Þ terms in Eqs. (31), (33) and (37) give o2 h2 ¼0 on2
ð94Þ
h2 ð1; sÞ ¼ 0
ð95Þ
oh2 ðc; sÞ ¼ r1 Bih1 ðc; sÞ þ r2 Bi on which can be solved as h2 ðn; sÞ ¼ Bi2 ð1 cÞr21 þ Bir2 ðn 1Þ
ð96Þ
ð97Þ
The Oð3 Þ terms of Eq. (36) give dðr1 r2 Þ 1 oh2 oh1 o/1 o/0 2 dr1 r1 ¼ ð1; sÞ r1 ð1; sÞ r1 ð0; sÞ r2 ð0; sÞ ds lð1 cÞ on ds on og og
ð98Þ
where /1 is the first-order temperature distribution in the liquid. In principle, this can be determined from o2 /1 o/1 o/ dr1 o/0 2r1 0 ð1 gÞ ¼ 2 og os os ds og o/1 ð1; sÞ ¼ 0 og
ð99Þ ð100Þ
748
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
/1 ð0; sÞ ¼ 0
ð101Þ
/1 ðg; 0Þ ¼ 0
ð102Þ
It is difficult to solve for /1 analytically due to the series form of /0 and r1 . Notice, however, the 2 homogeneous nature of the boundary conditions and the exp½ n þ 12 p2 s rate of decay with respect to time in the /0 and dr1 =ds series in Eqs. (81) and (91). We can conclude that the time rate of change of /1 ðg; sÞ and o/1 =ogðg; sÞ will have at least an exponential decay factor for large s, although it increases from zero to a finite value near the beginning. Using Eqs. (81), (89) and (97) in Eq. (98), the equation for r2 is dr2 dr1 1 Bi 1 o/1 ¼ r1 ð0; sÞ ð103Þ þ Bi r1 ds ds 1c l lð1 cÞ og which can be integrated to give Z s Z s 1 1 Bi 1 o/1 ð0; s0 Þ ds0 r1 ðs0 Þ ds0 þ Bi r2 ðsÞ ¼ r21 2 1c l s lð1 cÞ s og
ð104Þ
As s ! 1, r1 ! Bi s where the exponential term in Eq. (93) has dropped out. As discussed at the end of the last paragraph, the last integral is also negligibly small due to its exponential decay factor. The limit of r2 with large s is r2 !
Bi3 s2 2l2 ð1 cÞ
ð105Þ
We now look at the perturbation solution we have obtained up to now. As s ! lð1 cÞ=Bi, the interface position tends to r!
Bi s Bi3 s2 2 þ Oð3 Þ lð1 cÞ 2l ð1 cÞ
ð106Þ
This is not valid near r ¼ 1 since the higher order terms become comparable to those of lower order when s lð1 cÞ=Bi. So only the early stage of solidification is described well by this expression, but not the later stage. This is thus a small-time solution. 5.2. Large-time solution To obtain the large-time solution, we will follow the steps used in Section 4 and use the expansions (48)–(50). We change the time scale to bs ¼ s. Denoting the large-time solution by b r, b h, b , the governing equations (31)–(39) become and / " # b o2 b db r ob h h 2 oh ðn cÞb r ; c6n61 ð107Þ ¼ b r obs dbs on on2 b b b b o/ o2 / dr 2 o/ ð1 gÞð1 b rÞ ; ¼ ð1 b rÞ 2 obs dbs og og with boundary conditions
06g61
ð108Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
749
b h ð1; bs Þ ¼ 0
ð109Þ
b ð0; bs Þ ¼ 0 /
ð110Þ
b o/ ð1; bs Þ ¼ 0 og " # b 1 ob h o/ db r ð1 b r Þ ð1; bs Þ b r ð0; bs Þ ¼ b r ð1 b rÞ lð1 cÞ on og dbs h i ob h ðc; bs Þ ¼ b r Bi b h ðc; bs Þ þ 1 on
ð111Þ
ð112Þ
ð113Þ
From an Oð1Þ balance in Eq. (108), we obtain b o2 / 0 ¼0 og2
ð114Þ
b ð0; sÞ ¼ 0 / 0
ð115Þ
b o/ 0 ð1; sÞ ¼ 0 og
ð116Þ
which gives b ðg; sÞ ¼ 0 / 0
ð117Þ
This solution is only valid in the later stage of solidification, not near the beginning where bs ¼ bs . Eq. (117) implies that the temperature distribution in the liquid region decays to almost the freezing temperature when the outer solution becomes valid. A balance of Oð1Þ terms in Eq. (107) is h0 o2 b ¼0 on2
ð118Þ
b h 0 ð1; sÞ ¼ 0
ð119Þ
ob h0 ðc; sÞ ¼ b r 0 Bi½b h 0 ðc; sÞ þ 1 on
ð120Þ
which can be solved to give b h 0 ðn; bs Þ ¼
Bib r0 ð1 nÞ 1 þ ð1 cÞBib r0
ð121Þ
750
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
To this order, the interfacial condition (112) is db r0 1 ob h0 b r 0 ð1 b ð1 b r0Þ ¼ ð1; bs Þ r0Þ lð1 cÞ dbs on which can be reduced to
ð122Þ
db r0 1 Bi ¼ lð1 cÞ 1 þ ð1 cÞBib r0 dbs
ð123Þ
b r 0 ðbs 0 Þ ¼ 0
ð124Þ
b 0 at Since the actual initial condition is outside the range of validity of the large-time solution, r bs ¼ bs 0 is an arbitrary constant, where bs 0 will be determined later by matching with the small-time solution. The matching also shows that the constant can be taken to be zero without loss of generality. The solution of b r 0 is s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " # 1 1 1 2ðbs bs 0 Þ b ð125Þ þ r0 ¼ þ 1c Bi Bi2 l A balance of terms of Oð1=2 Þ in Eq. (108), gives b ðg; bs Þ ¼ 0 / 1 and a balance of terms of Oð
ð126Þ 1=2
Þ in Eq. (107), gives
o2 b h1 ¼0 on2
ð127Þ
hb1 ð1; sÞ ¼ 0
ð128Þ
ob h1 ðc; sÞ ¼ b r 1 Bi½b h 0 ðc; sÞ þ 1 þ Bib r0 b h 1 ðc; sÞ on
ð129Þ
The solution is Bib r 1 ½b h 0 ðc; sÞ þ 1 b h 1 ðn; bs Þ ¼ ðn 1Þ 1 þ ð1 cÞBib r0 The interfacial condition for Oð1=2 Þ is ob h1 db r1 db r0 ð1; bs Þ ¼ lð1 cÞ b r0 b r1 on dbs dbs
ð130Þ
ð131Þ
which can be reduced to b r0
db r1 Bib h 0 ðc; sÞb r1 ¼ lð1 cÞ dbs 1 þ ð1 cÞBib r0
Using Eqs. (121), (125) and (130), we find that
ð132Þ
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
b r1 ¼
a 1 þ ð1 cÞBir0
751
ð133Þ
where a is a constant of integration that will be determined by comparing with the corresponding small-time solution r2 after leading order matching. It can be noted that the outer solution is actually a one-phase solution with a uniform freezing temperature in the liquid. 5.3. Matching So far we have obtained small- and large-times solutions. We would now like to match these solutions in some suitable manner. To do this we propose the existence of an intermediate time, which is sm in the inner time variable and bs m in the outer time variable, at which the inner and outer solutions share common values of interface position and its first time derivative. Thus, our conditions for matching are r ðbs m Þ rðsm Þ ¼ b
ð134Þ
dr db r ðsm Þ ¼ ðbs m Þ ds dbs
ð135Þ
We emphasize that this procedure does not require matching of all derivatives as is commonly done. The leading order balance is r 0 ðbs m Þ r1 ðsm Þ ¼ b
ð136Þ
dr1 db r0 ðsm Þ ¼ ðbs m Þ ds dbs
ð137Þ
which can be written explicitly as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi# " 1 1 1 2ðbs m bs 0 Þ þ þ 2 1 c Bi Bi l " # !9 8 2 1 > > 2 > > > C exp n þ p ðs n m s Þ 1 > > > 1 < = 2 X ¼ Biðsm s Þ þ > > 1 lð1 cÞ > > n¼0 > > p nþ > > : ; 2 and
# " 2 1 X Bi 1 1 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ Bi Cn n þ p ðsm s Þ p exp n þ 2 2 n¼0 2Bi2 ðbs m bs Þ 1þ l
ð138Þ
ð139Þ
These two equations for bs m and bs 0 can be solved numerically. The functional forms of r1 ðsÞ and b r 0 to be tangent to each other at bs ¼ bs m . Fig. 6 is r 0 ðbs Þ guarantee the existence of a bs 0 for r1 and b
752
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762 10
10
10
2
τ c −τ∗
1
0
τ m−τ∗ 10
1
τ o −τ∗ 10
10
2
0
2
4
6
8
10
e ð0; s Þ ¼ 0:286, c ¼ 0:0627, l ¼ 0:582, e ¼ Bi, k ¼ 0:500, / Fig. 6. Variation of sm s , s0 s and sc s with Bi; Bi ¼ 0:0594.
calculated for the leading order terms. Only the first-order term r1 of the inner solution and the zeroth-order term b r 0 of the outer are used for the matching; sc is from Eq. (146). The figure shows sm s , s0 s and solidification completion time sc s (see Section 6.3 for definition of sc ) as functions of Bi. sm is larger than s0 , but both decrease as the heat transfer coefficient is increased. They are also much smaller than the solidification completion time. When considering the matching of the next term in the rðsÞ and b r ðbs Þ expansions, which is 1=2 order higher than the leading order terms r1 ðsÞ and b r 0 ðbs Þ respectively, we find that the constant in Eq. (133) turns out to be a ¼ 0, since the outer solution has a zero Oð3=2 Þ term. From this we know that both the next terms in rðsÞ and b r ðbs Þ expansions are OðÞ higher than their respective r 0 ðbs Þ. This means that the order of accuracy of both rðsÞ and b r ðbs Þ leading order terms r1 ðsÞ and b solutions is actually the same, being OðÞ higher than the leading term. After matching, the perturbation solution of the interface position is summarized as: 8 9 8 > > > " ( ) # > > > 2 > 1 = < X > Cn 1 > 2 > > þ Oð2 Þ; Biðs s Þ þ p ðs s Þ 1 exp n þ > < lð1 cÞ > > 1 2 > > n¼0 ; : nþ p r¼ 2 > ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s " # > > > > 1 1 1 2ðbs bs 0 Þ > > þ OðÞ; bs > bs m > : 1 c Bi þ Bi2 þ l
s 6 sm
ð140Þ
Fig. 7 shows an example of the matching of the small- and large-time solutions of the interface position. Again only the leading order solutions are shown. The forms of the solutions make it possible to always match both the position and its time derivative. Note that b r ðs0 Þ ¼ 0 corresponds to the intercept on the abscissa in this figure.
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
753
0.3
0.2
Outer solution
0.1 Matching point Inner Solution
0
τ=τm 0
το
1
2
3
4
e ð0; s Þ ¼ 0:286, c ¼ 0:0627, l ¼ 0:582, e ¼ Bi ¼ 1:0, k ¼ 0:500, / Fig. 7. Matching of inner and outer solutions; Bi ¼ 0:0594.
Since the matching principles (134) and (135) have been based only on the interface position and its velocity, we would like to check the corresponding matching of the temperature fields. r 0 ðbs m Þ in the expression (121) by r1 ðsm Þ, so that At s ¼ sm , we can replace b b h ðbs m Þ ¼
Bir1 ðsm Þ ½1 þ OðÞð1 nÞ 1 þ ð1 cÞBir1 ðsm Þ
ð141Þ
Expanding, we get b h ðbs m Þ ¼ Bir1 ðsm Þð1 nÞ þ Oð2 Þ
ð142Þ
which is identical to the inner solution hðsm Þ from Eq. (89). So the inner and outer solutions of the temperature distribution in the solid also match at bs m . Furthermore, the temperature distribution in the liquid and its time and spatial derivatives match at sm since the inner solution is exponentially small as seen from Eq. (81), and the outer solution is exactly zero as Eq. (117) shows. The derivatives oh=os, oh=on and o2 h=on2 are also matched at s ¼ sm according to (135) and the linear spatial form of (121) and (89), respectively. The largest temperature variations with time occur at the void–solid interface, and at the adiabatic liquid boundary. Fig. 8 shows the perturbation and numerical solutions of the temperature histories at these locations. /ð1; sÞ is the inner solution /0 ð1; sÞ, and hðc; sÞ is the inner h 0 ðc; sÞ when s > sm . The matching time sm solution h1 ðc; sÞ when s 6 sm and the outer solution b for the two time-scale solution is also indicated. For the history of the solid–liquid interface position, Fig. 9 shows the comparison between the numerical solutions from the enthalpy method and the perturbation solutions. The Bi ¼ 4:1 curve includes terms of the first three orders: r0 , r1 and r3 . For Bi ¼ 1 the inner solution r1 ðsÞ is used r 0 ðsÞ when s > sm . Both the regular perturbation solution of Section 4 when s 6 sm and the outer b
754
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
0.5
0.25
φ(1,τ) 0
θ(γ,τ) –0.25
τ=τ m –0.5
0
2.5
5
7.5
10
12.5
15
Fig. 8. Comparison of numerical and perturbation solutions for temperature history at boundaries: Lines are pere ð0; s Þ ¼ 0:286, e ¼ 4:10, k ¼ 0:500, / turbation solutions and points correspond to numerical solutions; Bi ¼ Bi c ¼ 0:0627, l ¼ 0:582, ¼ 0:0594.
1
0.75
Bi=4.1 Bi=1.0
0.5
0.25
0
0
2
4
6
8
10
12
14
Fig. 9. Comparison of numerical and perturbation solutions for interface positions: Lines are perturbation solutions e ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594. e c ¼ 0:0627, k ¼ 0:500, / and points correspond to numerical solutions; Bi ¼ Bi Bi,
(Bi ¼ 1=2 ¼ 4:1 is used) and the matched two time-scale solution of Section 5 (Bi ¼ 1 is used) agree well with numerical results within the expressed accuracy of the perturbation expansion.
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
755
6. Discussion 6.1. Extension of Bi range The two solutions of Sections 4 and 5 can be extended to cover the entire range of possible Bi. The Bi Oð1=2 Þ solution in Section 4 is also valid for Bi > Oð1=2 Þ. As a special case it is valid as Bi ! 1, which corresponds to isothermal cooling. In fact it can be checked that in this e ! 1, limit the solution calculated here tends to that given by Weinbaum and Jiji [7] if we take Bi Bi ! 1 and c ¼ 0. We can also consider the range of Bi between Oð1Þ and Oð1=2 Þ. The interface position computed from the Oð1=2 Þ regular perturbation or the Oð1Þ two time-scale solution are numerically similar, as shown in Fig. 10 with comparison to the numerical solution, so that either one may be used. In the figure, the regular perturbation solution includes terms of the first three orders: r0 , r1 and r2 . The matched asymptotic solution is the inner solution r1 ðsÞ when s 6 sm and the outer solution b r 0 ðsÞ when s > sm . For Bi < Oð1Þ, which corresponds to very weak cooling, the solution of Section 5 for Bi Oð1Þ remains valid. In obtaining the inner solution, the leading order balance for the boundary condition (37) could include OðÞ terms. 6.2. One-phase case Another interesting situation is when Ti ¼ Tf where the liquid is always at its freezing temperature. In this case, s ¼ Cn ¼ 0. To compare with the results of Yan and Huang [6], we take Bi Oð1Þ and c ¼ 0. The matching condition gives sm ¼ s0 ¼ 0 and the interface position is 1
0.75 Matched Asymptotic Solution
0.5
Regular Perturbation Solution
Numerical Solution
0.25
0
0
2
4
6
8
10
Fig. 10. Comparison between regular and two time-scale perturbation solutions in Oð1Þ < Bi < Oð1=2 Þ range: Lines e ð0; s Þ ¼ 0:286, e ¼ 2:0, k ¼ 0:500, / are perturbation solutions and points correspond to numerical solutions; Bi ¼ Bi c ¼ 0:0627, l ¼ 0:582, ¼ 0:0594.
756
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
1 r¼ þ Bi
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2bs þ þ OðÞ 2 Bi lð1 cÞ
ð143Þ
To the leading order this limit is identical to the solution in [6]. 6.3. Interface motion and solidification time The solutions in Sections 4 and 5 both show that the solid–liquid interface moves with an initial velocity of " # 1 X dr 1 ðs Þ ¼ Bi ð144Þ nþ pCn þ h:o:t ds lð1 cÞ 2 n¼0 where h.o.t are the higher order terms. The starting velocity is identical for all Bi, though to different orders of accuracy. From a physical point of view, this is more realistic than the infinite starting velocity predicted by the models of Wilson and Solomon [5] and Weinbaum and Jiji [7]. The subsequent interface motion after starting is much slower. Fig. 11 shows the displacement of the solidification front for different Bi. The Bi ¼ 20, 10 and 4.1 curves include terms of the first three orders: r0 , r1 and r2 . For Bi ¼ 1 and 0.5 the inner solution r1 ðsÞ is used when s 6 sm and the outer b r 0 ðsÞ when s > sm . As is intuitively obvious, the front is seen to move faster as the cooling is increased. The total time sc necessary to solidify the whole slab is of special interest in this problem. Taking r ¼ 1 in Eq. (140) for Bi Oð1Þ we find that " # 2 lð1 cÞ lð1 cÞ þ ð145Þ sc ¼ 1 þ bs 0 þ Oð1Þ Bi 2
1 Bi=20 10
4.1 1.0
0.75
0.5
0.5
0.25
0
0
5
10
15
20
25
e ð0; s Þ ¼ 0:286, c ¼ 0:0627, e k ¼ 0:500, / Fig. 11. Solid–liquid interface position for various Biot numbers, Bi ¼ Bi Bi, l ¼ 0:582, ¼ 0:0594.
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
757
where bs 0 is determined from the matching. Similarly, from the Bi Oð1=2 Þ solution, we have sc ¼
1
lð1 cÞ2 lð1 cÞ þ 1=2 D þ s þ OðÞ þ 1=2 Bi0 2
ð146Þ
where 2
sffiffiffi " " 2 # rffiffiffi# 3 n þ 12 p2 n þ 12 p l 2 6 7 erfc Bi0 exp pffiffiffi 1 6 7 Bi l 2 C 0 l X 2 6 7 Cn 61 þ D¼ ð1 cÞ 7 2 7 Bi0 n¼0 6 2 1 3=2 4 5 p nþ 2
ð147Þ
The leading orders of sc in Eqs. (145) and (146) are identical. For small , the solidification time may be explained in the following manner. We take a linear temperature profile in the solid and constant one in the liquid. The heat balance equation then gives Tf T0 ds ¼ ql L 1 ð1 cÞs dt þ h ks Integrating this equation and applying the conditions sð0Þ ¼ 0 and sðtc Þ ¼ ‘, we get ql L‘ ð1 cÞ‘L tc ¼ 1þ ðTf T0 Þh 2ks which in nondimensional terms is lð1 cÞ ð1 cÞBi 1þ sc ¼ Bi 2
ð148Þ
ð149Þ
ð150Þ
These are the first two terms of Eqs. (145) and (146). For large Bi, the limitation to interface motion is not the thermal resistance at the void–solid e ! 1, Bi ! 1, and interface but the conduction resistance within the solid. As Bi 2 sc ! lð1 cÞ =2 so that it becomes independent of Bi. The variation of sc with Bi is also graphed in Fig. 6. It may be noted that since sc is much larger than sm , the outer solution of Section 5 is the one that is applicable for most of the time. Figs. 12 and 13 show the effect of the shrinkage factor, c, on the total time taken to solidify the slab, sc . In Fig. 12 the first three orders, r0 , r1 and r2 , are used while for Fig. 13 the inner solution r 0 ðsÞ for s > sm . We see that the larger c is, the r1 ðsÞ is used for s 6 sm and the outer solution b smaller is sc . This is also indicated in the first terms of Eqs. (145) and (146). Mathematically this comes from the fact that increasing shrinkage reduces the thermal resistance of the solidified layer, so the same cooling is becoming more efficient as the solid–liquid interface advances. This incorrect prediction is due to the physically unrealistic assumption of a Biot number that is constant (as assumed by Reddy and Beckermann [8]) in spite of the growth of the void, a restriction that is lifted in the following discussion.
758
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762 1 γ=
0.2
0.1
7
0.0
62
0.0
0.75
0.5
0.25
0
0
2
4
6
8
e ¼ 4:10, Fig. 12. Shrinkage effect on solid–liquid interface position from regular perturbation solution with Bi ¼ Bi e k ¼ 0:500, / ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594. 1 γ=
0.2
0.1 7 62 0.0
0.0
0.75
0.5
0.25
0
0
4
8
12
16
e ¼ Fig. 13. Shrinkage effect on solid–liquid interface position from two time-scale perturbation solution with Bi ¼ Bi e ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594. 1:0, k ¼ 0:500, /
6.4. Time-dependent Bi So far we have assumed that Bi is a constant. However, it can be argued on physical grounds that the heat transfer coefficient must change as the void between the solid surface and the mold wall grows. The proper change will only occur if the thermal resistance of the void is not small
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
759
compared to that of the solid. The solutions derived here can also be applied to time-dependent Bi e as long as they vary slowly compared to the solidification process. and Bi Bi, As a example, let us take the heat transfer coefficient to be inversely proportional to the size of the void to compare with the results of Wilson and Solomon [5] who obtain a similarity solution for this case. Thus if we take h ¼ kv =csðtÞ, where kv is the thermal conductivity in the void, we have BiðsÞ ¼
kv crðsÞ
ð151Þ
where kv ¼ kv =ks . We note that at s ¼ 0, Bi ! 1 and thus we have s ¼ 0. Eq. (59) for large Bi reduces to vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2as t u ð152Þ s¼u t 2 1 ð1 cÞ 1 þ c kv in dimensional terms. This form is similar to that given in [5], and turns out to be the leading order of the solution obtained in [5] if the transcendental equation for the interface position given in that work is expanded in power series of 1=2 . The quantitative difference of the interface position between this solution and those in the numerical example presented there is within 2%. The model represented by Eq. (151) has the obvious drawback of predicting an infinite heat transfer coefficient at the beginning of solidification. One can improve it by adding a contact thermal resistance to remove this singularity. If heat transfer across the void is by conduction alone and natural convection and radiation are negligible, as is possible to justify in some applications of practical interest [9], the Biot number takes the form e Bi
ð153Þ e Bi 1 þ crðsÞ kv e due to small contact resistance. The above form of the time-dependent with a possibly large Bi Biot number is substituted for Bi0 into the regular perturbation solution (59) to give sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 2s þ þ ð154Þ r1 ¼ 2 2 e e b blð1 cÞ Bib Bi Bi BiðsÞ ¼
where
2 1 b¼1þc kv
ð155Þ
Fig. 14 shows the motion of the interface for different values of c calculated from r1 above. They agree reasonably well with the numerical solutions. As c increases, it takes longer for the slab to solidify and sc increases, as long as kv < 2 according to the above expression. This can be contrasted with the behavior shown in Fig. 8 which was for the constant Bi (i.e. for kv ! 1). In practice [9], kv 1 (for our case it is 2.73 · 104 ); the interface motion is very sensitive to kv , as shown in Fig. 15 where r1 is calculated from Eq. (154).
760
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762 1 0.75 0.5
γ=0.0627
0.25 0
0
2
4
1
τ−τ∗
6
8
10
6
8
10
0.75 0.5 γ=0.20
0.25 0
0
2
4
Fig. 14. Shrinkage effect on solid–liquid interface position from time varying Bi model: Lines are numerical solutions e ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594, kv ¼ 0:50. e ¼ 4:10, k ¼ 0:50, / and points correspond to analytical expressions; Bi 1 γ=0.0627 0.1 0.75
0.2
0.5
0.25
0
0
1000
2000
3000
4000
5000
6000
e ð0; s Þ ¼ 0:286, l ¼ 0:582, ¼ 0:0594, kv ¼ 2:73 104 . e ¼ 4:10, k ¼ 0:50, / Fig. 15. Solid–liquid interface position, Bi
7. Conclusions The one-dimensional solidification problem in a finite domain has been solved using perturbation methods with the Stefan number as the small parameter. Initially, the liquid is at a uniform
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
761
temperature above the freezing temperature. The driving mechanism for solidification is cooling at one boundary, while the other end of the slab is adiabatic. The void that is formed due to density change during solidification is assumed to appear between the mold wall and the cold end of the solid. The two-phase problem with convective boundary condition is more complicated to handle than either the isothermal condition analyzed in [7], or the one-phase problem studied in [6]. The convective mechanism gives rise to a nondimensional parameter, the Biot number Bi, in the governing equations. The nature of the solution depends on the order of magnitude of Bi in terms of the perturbation parameter, the Stefan number . There is a precooling stage after which solidification begins. The temperature distribution in the precooling stage is obtained as a series solution. The perturbation method is used once solidification actually begins. The range of perturbation solutions of the interface position and temperature distributions in both solid and liquid that are obtained for different magnitudes of cooling are the following: 1. Oð1=2 Þ 6 Bi < 1: This corresponds to strong cooling for which a regular perturbation solution is obtained. 2. Oð1Þ < Bi < Oð1=2 Þ: This corresponds to intermediate cooling. The solid–liquid interface position computed from the solutions in cases (1) above or (3) below are not much different. Either can be used for numerical purposes. 3. 0 < Bi 6 Oð1Þ: This corresponds to weak and moderate cooling. A two time-scale solution describes the interface motion and temperature distribution. The inner solution is valid up to an intermediate time sm . From sm onwards, the outer solution takes over. The solutions approach those of Refs. [5–7] in the appropriate limits. The analysis is also validated using numerical solutions obtained by the enthalpy method. The results obtained here are valid not only for constant Biot numbers, but can also be applied to time-dependent parameters as long as they are slowly varying. This greatly extends the usefulness of the present results and enables them to be tailored to any void thermal resistance or wall temperature model.
Acknowledgement We gratefully acknowledge support of this research by A.E. Goetze of South Bend, Indiana.
References [1] [2] [3] [4] [5]
H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Oxford University Press, London, 1959. J. Crank, Free and Moving Boundary Problems, Oxford University Press, New York, 1984. L.I. Rubinstein, The Stefan problem, Trans. Math. Monographs 27 (1971) 94–181. M.A. Biot, Variational Principles in Heat Transfer, Clarendon Press, Oxford, 1970. D.G. Wilson, A.D. Solomon, A Stefan-type problem with void formation and its explicit solution, IMA J. Appl. Math. 37 (1986) 67–76.
762
Z. Yang et al. / Appl. Math. Modelling 27 (2003) 733–762
[6] M.M. Yan, P.N.S. Huang, Perturbation solutions to phase change problem subjected to convection and radiation, ASME J. Heat Transfer 101 (1979) 96–100. [7] S. Weinbaum, L.M. Jiji, Singular perturbation theory for melting or freezing in finite domains initially not at fusion temperature, ASME J. Appl. Mech. 44 (1977) 25–30. [8] A.V. Reddy, C. Beckermann, Measurements of metal-mold interfacial coefficients during solidification of Sn and Sn–Pb alloys, Exper. Heat Transfer 6 (1993) 111–129. [9] Z. Yang, Analytical and numerical study of solidification by convective cooling with void formation, M.S. Thesis, University of Notre Dame, Notre Dame, Indiana, 1994. [10] V. Alexiades, A.D. Solomon, Mathematical Modeling of Melting and Freezing Process, Hemisphere Publ. Corp., Washington, DC, 1993.