Approximate analytical and numerical solutions for a two-dimensional Stefan problem

Applied Mathematics and Computation 202 (2008) 857–869

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Approximate analytical and numerical solutions for a two-dimensional Stefan problem Faruk Yigit Department of Mechanical Engineering, King Saud University, P.O. Box 800, Riyadh 11421, Saudi Arabia

a r t i c l e

i n f o

Keywords: Stefan problem Solidification Phase-change Perturbation Heat conduction

a b s t r a c t This paper proposes an approximate analytical and a numerical solution method to a twodimensional heat conduction problem in which a liquid becomes solidified by heat transfer to a planar mold surface by using a linear perturbation method. It is assumed that the cooling rate is perturbed by a small spatially sinusoidal heat flux at the shell–mold interface. This leads to a corresponding undulation of the solidified shell thickness. Approximate analytical results are obtained for the solid/melt moving interface as a function of time and for the temperature field in the shell. The approximate analytical solution is compared with a numerical solution, and a very good agreement has been found. A limiting analytical solution in which diffusivity of the solidified shell material is assumed to be infinitely large is also obtained, and compared with the numerical predictions to establish the validity of the model and the numerical approach. It is demonstrated that solidified shell materials with higher thermal diffusivities may result in irregular growth of the shell thickness which, generally, causes cracking near the surface. Ó 2008 Elsevier Inc. All rights reserved.

1. Introduction Heat transfer problems with phase-change are very common in physics and engineering. Typical examples include the production or melting of ice, solidification of castings, and aerodynamic heating of missiles. All of these problems share the characteristic of an interface boundary which moves toward either the solid (melting) or the liquid region (solidification). In these problems, the thermal behavior is assumed to be governed by a well known partial differential equation of heat conduction. Even though this differential equation is linear, the above problems are nonlinear due to the presence of the moving interface. There is extensive research on the above and related phase-change problems. A number of studies dealing with analytical and numerical aspects of particular melting or solidification heat transfer problems have appeared in the literature [1–4]. The non-linear nature of governing equations limits the number of exact analytical solutions to such problems. Many approximate analytical techniques such as the power series expansion [3], heat balance integral [5–8], similarity [9], isotherm migration [10], and source and sink [11] methods, however, have been developed and employed to obtain approximate solutions. Another common approach is to use the enthalpy method, which introduces an enthalpy function. In this case, the heat flux condition is automatically satisfied across the phase front [12]. Crank [13], Kutluay [14] and Caldwell and Savovic [15] are applied boundary immobilization methods using the spatial coordinate transformation in order to fix the moving interface. On the other hand, various numerical methods are also developed and applied to the problem successfully [16–18]. Crank [19] presents an elaborate collection of numerical methods used for these problems. Several other approximate methods of solution have been proposed for solving the phase-change problems. For example, Lardner [20] and Prasad and Agrawal [21] obtained an approximate solution for small values of time by using Biot’s E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.03.033

858

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

Nomenclature K c k L m Q S b S s T T Tb Tm t x, y X, Y b Y erfc(x)

thermal conductivity (W/m °C) thermal capacity (m2/s2 °C) thermal diffusivity (m2 s1) latent heat of fusion (W s/kg) reciprocal wavelength (m1) mold–shell interface heat flux dimensionless solidified shell thickness (see Eq. (9)) dimensionless solidified shell thickness (see Eq. (40)) solidified shell thickness (m) temperature (°C) dimensionless temperature (see Eq. (9)) dimensionless temperature (see Eq. (41)) melting temperature (°C) time (s) Cartesian coordinates dimensionless Cartesian coordinates dimensionless Cartesian coordinate (see Eq. (40)) error function K , qk

Subscript 0, 1 subscripts used for the zeroth- and first-order processes, respectively Greek symbols q density (kg/m3) b dimensionless time (see Eq. (9)) g dimensionless time (see Eq. (40)) Q0 1/f ; modified Stefan number (see Eq. (11)) mkqL d space step size s time increment

variational principle. El-Genk and Cronenberg [22] presented an analytical solution for the freezing of a stagnant liquid in a semi-infinite slab geometry subject to a constant heat flux boundary condition, and claimed that it was the first exact solution. However, Cho and Sunderland [23] showed that it was not exact since the suggested solutions for the solidification moving front do not satisfy the heat conduction equation. They presented an approximate analytical solution for the same problem based on the idea of decoupling the boundary conditions. A common drawback of the approximate techniques followed in the above studies is that they either are limited to one-dimensional analysis or result in very complicated mathematical formulations when applied to multidimensional problems. Furthermore, most of these approximate analytical techniques have difficulty in accommodating time dependent applied surface conditions and their applicability is restricted to short times. To overcome these difficulties, numerical methods, notably finite difference [24,25] and finite element methods [26,27] have also been used in the solution of many problems. In the casting processes, solidification is affected by heat transfer from melt through the mold to the surroundings. The process therefore generally starts at the mold surface and a shell of solidified material grows from this surface towards the center of the casting. The growth of this shell is often found to be non-uniform [28], even in nominally one-dimensional plane cast surface. This non-uniform shell growth could be a response to a pre-existing non-uniformity in the heat flux drawn from the mold surface. Small spatial perturbations in heat extraction, resulting from oxide films carried on the molten material surface, isolated impurities, or mold surface topography, for example, can lead to shell growth instability during solidification of metals [28]. This is a condition where shell distortion promotes gap nucleation beneath the thinnest regions of the shell while simultaneously enhancing heat extraction beneath the thickest regions of the shell due to improved contact. If the solidification process is interrupted, periodic thickness non-uniformities at the freezing front are often observed that can have wavelengths of order of several centimeters [29]. If solidification is allowed to proceed, however, the non-uniformities tend to die out as the freezing front morphology becomes less dependent upon the mold–shell interface due to an ever-thickening shell. Such microstructures are often detrimental to subsequent forming processes and have been linked to severe ingot cracking. It should however be noted that the heat flux itself is also influenced by the contact conditions between the shell and the mold. Thus, it is possible that a small spatial perturbation in the heat flux could be unstable due to thermo-mechanical coupling. Instabilities of this kind have already been reported in many theoretical [29,30,35] and experimental studies [31–34]. In this paper, the effect of thermal diffusivity of the solidifying shell material on the evolution of the shell thickness is assessed by extending the one-dimensional solution method suggested by Cho and Sunderland [23] to a two-dimensional

859

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

problem where a prescribed constant cooling rate is perturbed by a small spatially sinusoidal heat flux at the surface of a planar mold. We will particularly seek an approximate analytical solution for the advance of the solid/melt boundary in solidification of pure materials with finite thermal diffusivities. Furthermore, a numerical solution will be provided to check the accuracy of the proposed approximate solution. It is noteworthy that the present heat conduction model can be combined to a thermo-elastic deformation model to determine the development of the residual stresses and the contact pressure. The stability of thermo-elastic contact problem can then be examined since the perturbation method described here provides stability criterion, which models the physics of the growing solid and its deformation. Earlier models of solidification thermo-mechanics suggest that growth instability can be minimized or completely avoided by a careful choice of process parameters. The suggested model in the present work enables us to investigate how the variation of heat flux might affect the solidification rate and the thickness non-uniformity. The paper is organized as follows. In Section 2, the governing equations that need to be solved for the solidification process are summarized. Section 3 introduces dimensionless variables used in the analysis. In Section 4, the perturbation technique which is utilized to reduce the dimension of the problem is presented. Section 5 provides the methodology of the approximate solution of the problem. A numerical algorithm used in the numerical solution of the problem to check the accuracy of the approximate solution is discussed in Section 6. Section 7 demonstrates the limiting problem and its solution in order to develop and check the purely numerical solution of the general problem. Results and their implications on the solidification process are presented in Section 8. Finally, Section 9 provides a summary and conclusions of the work. 2. Formulation We consider the two-dimensional problem of a liquid initially at its melting temperature Tm, in the half space y > 0. It is assumed that the liquid is in contact with a plane mold at y = 0, on which the heat flux is a prescribed function of one space coordinate and time t. The present model is restricted to the solidification of pure materials. Therefore, solidification is assumed to occur at a distinct temperature rather than over a range of temperatures. The liquid phase being at the uniform melting temperature Tm throughout, there is no heat transfer through it; the heat released during solidification process is transferred into the mold, and raises its temperature. After a time t, the liquid solidified near the mold forms a solidified shell of thickness s(x, t). Thus, s(x, t) defines the moving interface between the solid and liquid phases as shown in Fig. 1. The prescribed heat flux is assumed to be nearly uniform in x so that the perturbation from unidirectional solidification, s = s0(t), is small. We assume that density q, thermal diffusivity k, and conductivity K, of the solid phase are constant and independent of temperature and time. The temperature of the solidified shell, T(x, y, t), must satisfy the heat conduction equation [1]: 1 oT r2 Tðx; y; tÞ ¼ ð1Þ k ot with the following initial and boundary conditions: sð0Þ ¼ 0;

ð2Þ

Tðx; s; tÞ ¼ T m ; oT os ðx; s; tÞ ¼ Lq ðx; tÞ; K oy ot oT ðx; 0; tÞ ¼ Q ðx; tÞ; K oy

ð3Þ ð4Þ ð5Þ

where y = s(x, t) defines the instantaneous position of the solid/melt boundary. Eq. (3) states that the freezing front is isothermal at the fusion temperature, while (4) defines an energy balance between the heat conducted away from the moving interface into the solidified shell and the latent heat released during solidification. Finally, Eq. (5) prescribes the heat flux drawn from the mold surface. y 2π /m Liquid

T=Tm

s1(t) Solid

s(x,t)

s0(t)

x Mold Q=Q0+Q1cos(mx) Fig. 1. Geometry of non-uniform directional solidification.

860

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

Notice that the smallness of the perturbation from unidirectional solidification implies that os/ox  1 and hence permits us to write the boundary condition (4) in terms of x and y rather than in a slightly rotated coordinate system normal to the local melt front. We consider the case in which the prescribed heat flux is of the form: Q ðx; tÞ ¼ Q 0 þ Q 1 cosðmxÞ;

Q 1  Q 0;

ð6Þ

where Q0 and Q1 are assumed to be time-independent constants. The second part of Eq. (6) defines a small sinusoidal perturbation on the constant heat flux Q0 at the shell–mold interface. Since Q1  Q0, we anticipate a corresponding small sinusoidal perturbation in the temperature field and the solidified shell thickness, which can be expressed by the equations: Tðx; y; tÞ ¼ T 0 ðy; tÞ þ T 1 ðy; tÞ cosðmxÞ;

ð7Þ

sðx; tÞ ¼ s0 ðtÞ þ s1 ðtÞ cosðmxÞ:

ð8Þ

We assume that the amplitude of this perturbation is small in comparison with its wavelength (i.e., ms1(t)  1), in which case the slope of the moving front, os/ox, is very much less than unity. 3. Dimensionless presentation Before proceeding to the solution, it is convenient to introduce the following dimensionless variables: b¼

Q 20 kðqLÞ2

TðY; bÞ ¼ f¼

t;

X ¼ mx;

K Tðy; tÞ; kqL

Y¼

Q¼

Q0 y; kqL

SðbÞ ¼

Q0 sðtÞ; kqL

Q ; Q0

ð9Þ ð10Þ

mkqL : Q0

ð11Þ

Note that 1/f is a modified Stefan number for the problem of prescribed interfacial heat flux [1]. In the case of prescribed constant temperature, the Stefan number is Ste ¼

cDT ; L

ð12Þ

where cDT is the sensible heat or the heat that is stored in the material due to an imposed temperature gradient. It could also be argued that modified Stefan number [35] for the present problem is Stemodified ¼

Q 0 =qmk ; L

ð13Þ

where Q0/qmk is an analogue of the sensible heat defined above. The heat conduction equation (1) and the associated initial and boundary conditions (2)–(5) then become: o2 TðY; bÞ 2

oX Sð0Þ ¼ 0;

þ

o2 TðY; bÞ f2 oY 2

¼

oTðY; bÞ f2 ob

;

ð14Þ ð15Þ

TðX; S; bÞ ¼ T m ;

ð16Þ

oT oS ðX; S; bÞ ¼ ðX; bÞ; oY ob

ð17Þ

oT ðX; 0; bÞ ¼ Q ðX; bÞ: oY

ð18Þ

4. Method of solution If the mold is at uniform temperature, there is clearly a one-dimensional solution to the problem described above, in which the solid/melt moving interface at any given time is defined by the straight line Y = S0(b) and all physical quantities are functions of Y, b, but are independent of X. If an X-dependent perturbation is to grow on this unperturbed or ‘zeroth-order’ process, its development will be influenced by the instantaneous zeroth-order quantities. The perturbation is assumed to be sufficiently small for the corresponding two-dimensional problem and will therefore be solved by a linear perturbation of the zeroth-order solution. Notice that temperatures depend on Y and b only. Suffix 0 refers to the x-independent zeroth-order process, and the problem has a simple one-dimensional solution which is called the ‘‘zeroth-order” solution. The suffix 1 refers to the amplitude of the sinusoidal perturbation or first-order process. Substituting dimensionless form of Eq. (7) into Eq. (14), and separating periodic and uniform terms in X, we obtain: o2 T 0 ðY; bÞ oY 2

¼

oT 0 ðY; bÞ ; ob

ð19Þ

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

o2 T 1 oY 2

ðY; bÞ  f2 T 1 ðY; bÞ ¼

oT 1 ðY; bÞ : ob

861

ð20Þ

Since the perturbation is small, we can expand the temperature field in the vicinity of the mean solid/melt interface position, Y = S0(b), in the form of a Taylor series, in which case boundary condition (16) can be written as " # oT 0 ðS0 ; bÞ oT 1 ðS0 ; bÞ S1 ðbÞ cosðXÞ þ    þ T 1 ðS0 ; bÞ þ S1 ðbÞ cosðXÞ þ    cosðXÞ ¼ T m : T 0 ðS0 ; bÞ þ ð21Þ oY oY Separating periodic and uniform terms and dropping second and higher order and product terms in the small quantities, T 1 and S1, we obtain the two first-order equations: T 0 ðS0 ; bÞ ¼ T m ;

ð22Þ

oT 0 ðS0 ; bÞ þ T 1 ðS0 ; bÞ ¼ 0: S1 ðbÞ oY

ð23Þ

A similar procedure applied to the boundary condition (17) also yields the two equations: dS0 ðbÞ oT 0 ¼ ðS0 ; bÞ; db oY

ð24Þ 2

dS1 ðbÞ oT 1 ðS0 ; bÞ o T 0 ðS0 ; bÞ ¼ þ S1 ðbÞ : ob oY oY 2

ð25Þ

Finally, the boundary condition (18) at the mold surface gives the two equations: oT 0 ð0; bÞ ¼ 1; oY oT 1 ð0; bÞ ¼ Q 1 : oY

ð26Þ ð27Þ

Thus, the problem is reduced to the determination of two pairs of functions T 0 ðY; bÞ, S0(b) and T 1 ðY; bÞ, S1(b) in Eqs. (19) and (20), which satisfy the boundary conditions (22), (24), (26) and (23), (25), (27), respectively. 5. Approximate solution of the problem Eqs. (20) and (21) with boundary conditions in (22)–(24) and (25)–(27) cannot be solved in closed form. Therefore, we will obtain an approximate solution for this problem as explained below. 5.1. Zeroth-order solution Derivation of the approximate solution for the zeroth-order problem outlined above is readily available in Ref. [23].1 We therefore summarize only the final results referring the reader to Ref. [23] for more details. In Ref. [23], the temperature distribution was assumed to have a similar profile to that without phase-change [1], i.e., " # rffiffiffiffiffiffiffiffiffiffiffiffiffiffi Y b þ b0 Y 2 =4ðbþb0 Þ 2 ð28Þ T 0 ðY; bÞ ¼ T m þ Y erfc pffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 þ ; e p p 2 b þ b0 where b0 = 1/p. This expression satisfies the heat conduction equation (19) and the boundary condition (26), but it cannot be made to satisfy boundary conditions (22) and (24) simultaneously, and hence does not constitute an exact solution. Imposing boundary conditions (22) alone gives: ( " #)1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi S2 ðbÞ=4ðbþb Þ dS0 ðbÞ S0 ðbÞ 0 ¼ pðb þ b0 Þ e 0 erfc pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð29Þ db 2 b þ b0 for S0(b), whereas imposing (24) alone gives: " # dS0 ðbÞ S0 ðbÞ ¼ erfc pffiffiffiffiffiffiffiffiffiffiffiffiffiffi : db 2 b þ b0

ð30Þ

Eqs. (28) and (29) define an asymptotic solution of the zeroth-order problem where the energy balance at the solidification front is relaxed, while Eqs. (28) and (30) define an approximate solution of the zeroth-order problem where temperature continuity at the solidification front is relaxed.

1 It should be noted that several other approximate analytical methods are available in the literature for the zeroth-order problem defined above (see, for example [36,37]).

862

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

It should be noted that Cho and Sunderland [23] have not been able to prove rigorously that the exact solution would lie between these two approximate solutions. However, they have stated that they made numerous calculations and comparisons that show this was to be anticipated. Furthermore, they expressed the opinion that the condition (30) would lead to a better approximation than (29). These hypotheses are supported by the numerical results presented in the next section. In particular, we found that the maximum error in S0(b) resulting from the use of Eq. (28) with condition (30) is 1%. The results of both approximations are shown in comparison with the finite difference solution of the problem in Fig. 2. 5.2. First-order solution Cho and Sunderland’s method can also be used to obtain an approximate solution to the first-order problem defined by Eq. (20) with boundary conditions (23), (25) and (27) and initial temperature zero. The first-order temperature distribution T 1 ðY; bÞ that satisfies the governing equation (20) and the boundary condition (27) can be obtained as Z b f2 ðbsÞY 2 =4ðbsÞ Q e pffiffiffiffiffiffiffiffiffiffiffi ds: ð31Þ T 1 ðY; bÞ ¼  pffiffiffi 1 pQ 0 0 bs Following Cho and Sunderland (the derivation of T 1 ðY; bÞ can be found in the Appendix), the first-order solution front S1(b) can be approximated by imposing one only of the remaining conditions (23) or (25). Imposing the condition (23) gives: ! Z b f2 ðbsÞS2 =4ðbsÞ 0 S0 Q e pffiffiffiffiffiffiffiffiffiffiffi ds ¼ 0; ð32Þ S1 erfc pffiffiffiffiffiffiffiffiffiffiffiffiffiffi  pffiffiffi 1 pQ 0 0 bs 2 b þ b0 that can be rewritten as

S1 ðbÞ ¼

pQ ffiffipQ1 0

Rb

f2 ðbsÞS2 =4ðbsÞ 0

e

0

erfc



pffiffiffiffiffiffi bs

ds

 0 pSffiffiffiffiffiffiffiffi

2

ð33Þ

;

bþb0

where we have used: ! oT 0 S0 ðS0 ; bÞ ¼ erfc pffiffiffiffiffiffiffiffiffiffiffiffiffiffi : oY 2 b þ b0

ð34Þ

Imposing the condition (25) gives: 2 dS1 ðbÞ S1 ðbÞ Q þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi eS0 =4ðbþb0 Þ ¼ pffiffiffi1 db 2 pQ 0 pðb þ b0 Þ

Z 0

b

S0

2

ðb  sÞ3=2

ef

ðbsÞS20 =4ðbsÞ

ds:

ð35Þ

Similar to the zeroth-order solution described in previous section, Eqs. (31) and (33) define an asymptotic solution of the first-order problem where the condition (25) is relaxed, while Eqs. (31) and (35) define an approximate solution of the first-order problem where the condition (23) is relaxed. The resulting approximations are compared with the numerical solution of the problem in Fig. 3 for the case of f = 0.5.

25 Soln. with Eq. (29) Numerical Soln. Soln. with Eq. (30)

20

S0(β )

15

10

5

0 0

20

40

β

60

80

100

Fig. 2. Variation of the dimensionless zeroth-order solidification front with time, b, for f = 0.5.

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

863

0.3 Soln. with Eq. (33) Numerical Soln. Soln. with Eq. (35)

S1(β)

0.2

0.1

0 0

10

20

30

40

S0(β) Fig. 3. Variation of the dimensionless first-order solidification front with the dimensionless mean shell thickness, S0(b), for f = 0.5.

6. Numerical solution To check the accuracy of the approximate solution presented in previous section Eqs. (20) and (21) with boundary conditions in (22), (24), (26) and (23), (25), (27) are solved numerically. An explicit finite difference Lagrangian scheme, previously developed in [24], is used in this work. In this case, the zerothorder solid phase 0 < Y < S0(b), is divided into a fixed number of elements N, so the space step size, dj ¼ Sj0 =N, increases with time due to the growth in S0(b). This permits the last node to be identified with the zeroth-order solid–liquid moving front at all times, but implies that the node locations move in time, necessitating the inclusion of convective terms in the updating algorithm for temperature. Thus, for example, the instantaneous zeroth-order solidified layer temperature field is represented by the temperatures at the N + 1 points Y = (i  1)d, i = 1, 2, . . ., N + 1. The increase of the shell thickness S0(b) during the next time increment s is determined from the finite difference formulation of Eq. (24). ¼ Sj0 þ Sjþ1 0

s 2dj

½3T j0Nþ1  4T j0N þ T j0N1 ;

ð36Þ

after which the solidified layer temperatures at the interior nodes i = 2, 3, . . ., N are updated using the finite difference form of heat conduction equation (19): jþ1 T jþ1 0i ¼ T 0i þ

s ðdj Þ2

½T j0iþ1  2T j0i þ T j0i1 ;

i ¼ 2; 3; . . . ; N;

ð37Þ

which can be corrected using convective terms through: jþ1 T jþ1 0i ¼ T 0i þ ði  1Þ

djþ1  dj dj

jþ1 ½T jþ1 0iþ1  T 0i ;

i ¼ 2; 3; . . . ; N:

ð38Þ

The solidified layer temperatures at node N + 1 remain at a constant value for all times in view of Eq. (22) and that at node 1 is updated through (26), which determines the first difference in the first element. Essentially, the same procedure is used to determine the evolution of the first-order solidified layer temperature field, using Eqs. (23), (25) and (27). The choice of an appropriate value for s is motivated by the desire for computational efficiency, while retaining acceptable numerical convergence and stability. Extensive investigations were made into the effect of both the space and time discretization to ensure that the final results are reliable. With the explicit scheme used here, the maximum time step for stability is proportional to d2 and hence the stability requirement generally places a restriction on s when good spatial accuracy is desired, necessitating very small values of d. The numerical stability of the scheme [38] is satisfied if: s d2

< 0:5;

ð39Þ

where d is the smallest dimension in the discretization. At the beginning of the process, the thickness of the solidified layer and d are very small, and therefore we need an extremely small time step that causes the numerical process to be very slow. However, as solidification progresses and the shell thickens, the time steps increase, leading to an accelerated convergence of the numerical process. In the present work, we developed the numerical algorithm that runs faster without an accompanying loss of accuracy. At the early stages of the process, the thickness of the solidified shell is very small and the temperature distribution is approximately linear. Thus, a few points are sufficient to represent the temperature field accurately. To use

864

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

fewer elements at the early stages allows us to choose larger time increments rather than using a fixed number of elements through the whole solidification process. As the solidified layer thickens, we need to increase the number of grid points for the accuracy of the process. Therefore, new grid points need to be included and the solidified shell must be rediscretized according to the new number of elements. In this case, the temperatures at the new grid points are determined by using a linear interpolation. With the algorithm described above, it is clearly not possible to start at the instant of first solidification, since at S0(b) = 0, all the nodes would coincide. Instead, we need to use an asymptotic solution of the problem at small times to provide a suitable initial condition for the numerical algorithm at finite time. Fortunately, the limiting solution given in the next section, which assumes that diffusivity of the solidified shell material is infinitely large, becomes progressively more accurate at small times, since the temperature drop across the solidified layer is small at the very beginning of the process. We can therefore start the process with a small but finite thickness, using the limiting solution given in the next section to define the initial values for the temperature fields in the solidified shell thickness. 7. The limiting solution for zero Stefan number It can be demonstrated that the solution for f ? 1 is a limiting case of the present more general theory. This simplification permits us to obtain an analytical solution in closed form. The results of this limiting case are useful in the development and checking of the numerical solution of the general problem. In physical terms, this simplifying assumption is equivalent to the statement that the solidifying shell material has zero thermal capacity. In other words, the heat diffusivity of the solidifying shell material is infinitely large. To obtain the solution in the limit as f ? 1, it is convenient to define the new dimensionless parameters2 as b ¼ my ¼ fY; Y

b SðgÞ ¼ msðtÞ ¼ fSðbÞ;

g¼

mQ 0 t ¼ fb; qL

bðY b ; gÞ ¼ mK Tðy; tÞ ¼ fTðY; bÞ: T Q0

ð40Þ ð41Þ

The governing equation for the zeroth-order temperature field given by Eq. (19) reduces to: b0 b 0 1 oT o2 T ; ¼ b2 f og oY

ð42Þ

and, in the limit f ? 1, it reduces to: b0 o2 T ¼ 0: b2 oY

ð43Þ

We can also express the boundary conditions (22), (24) and (26) in the new notation as follows: b 0 ðb S 0 ; gÞ ¼ T m ; T

b 0 ðb db S 0 ðgÞ o T S 0 ; gÞ ¼ ; b ds oY

b0 oT ð0; gÞ ¼ 1: b oY

ð44Þ

The solution of the governing equation (43) with the boundary conditions given in Eq. (44) will be obtained as b S 0 ðgÞ ¼ g; b 0 ðY b ; gÞ ¼ T bm þ Y b  g: T

ð45Þ ð46Þ

Thus, the zeroth-order temperature field is always linear in the limit f ? 1, and the moving interface advances linearly in time. However, the governing equation for the first-order temperature given by Eq. (20) reduces to: b1 b o2 T b 1 ¼ 1 oT 1 ; T b2 f og oY

ð47Þ

and, in the limit f ? 1, it reduces to: b1 o2 T b1 ¼ 0 T b2 oY

ð48Þ

with solution: b ; gÞ ¼ c1 ðgÞ sinhð Y b Þ þ c2 ðgÞ coshð Y b Þ; b 1 ðY T

ð49Þ

2 b 0 is then a function It should be noted that it would be possible to solve the entire problem in terms of these parameters. However, the temperature field T b and time g. This causes the numerical stability of the scheme is to be a function of f, which leads to of the parameter f as well as the dimensionless position Y extremely small time step size for large values of f. By contrast, T 0 depends on Y and b only [see Eqs. (19), (22), (24) and (26)], and hence the numerical stability of the scheme is independent of f. Therefore, this presentation is preferred in the main body of the work.

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

865

b , except in a limiting sense when Y b  1. To complete the solution, we express the boundary condiwhich is not linear in Y tions (23), (25) and (27) in the new notation and take the limit as f ? 1, with the result: b b 1 ðb S 0 ; gÞ; S 1 ðgÞ ¼  T

b 1 ðb db S 1 ðgÞ o T S 0 ; gÞ ¼ ; b dg oY

b 1 ð0; gÞ Q 1 oT ; ¼ b Q0 oY

ð50Þ

b 0 . Using Eqs. (49) and (50) and where we have used Eqs. (45) and (46) to substitute for the zeroth-order functions b S 0 and T eliminating the two arbitrary functions c1(g), c2(g), we obtain the ordinary differential equation: db S 1 ðgÞ Q1 þ tanhðgÞb S 1 ðgÞ ¼ dg Q 0 coshðgÞ

ð51Þ

for the amplitude of the perturbation in the solidification front, with solution: b S 1 ðgÞ ¼

Q 1g : Q 0 coshðgÞ

ð52Þ

The functions c1(g), c2(g) can then be found by back substitution, the resulting first-order temperature field being: b ; gÞ ¼ Q 1 fsinhð Y b Þ  ½g þ tanhð Y b Þð1  g tanhð Y b ÞÞ coshð Y b Þg: b 1 ðY T Q0

ð53Þ

8. Results and discussion Fig. 4 reports the position of the solid–melt interface at various values of the dimensionless time, b. Once S0(b) and S1(b) is found S(X, b) can be obtained using dimensionless form of Eq. (8). It indicates that shell growth will be faster at the points where the heat flux is greatest. Both approximations give the same results when b = 3 (which corresponds to S0(b) = 1.9). As the solidification progresses the shell thickness non-uniformity increases. If solidification is interrupted, periodic thickness non-uniformities at the freezing front will be observed. If solidification is allowed to proceed, however, the non-uniformities tend to die out as the freezing front morphology becomes less dependent upon the mold–shell interface due to the everthickening shell as seen in Fig. 5. It should be noted that this theoretical result is supported by many experimental observations of thermo-mechanical growth instability in casting processes (see, for example [31–34]). It is anticipated that this non-uniform shell growth could be a response to a pre-existing non-uniformity in the cooling rate. It should be noted that it is not possible to draw any definitive conclusions from the figures represented so far about the effect of thermal diffusivity on the evolution of shell thickness simply because of the definition of dimensionless variables presented in Section 3. New definitions of dimensionless variables given in Section 7, however, enable us to examine the influence of this key parameter on the process. Notice that the dimensionless parameter f stands for the thermal diffusivity of the solidified shell material since k appears only in the definition of f in the new dimensionless presentation in Section 7. In order to study the effect of thermal diffusivity on the development of shell thickness, we report Figs. 6 and 7: Fig. 6 shows the evolution of the solution of dimensionless mean shell thickness S0(g) for different values of f. The five sets of curves correspond to f = 1, 5, 10, 50, 500. There are two curves for each value of f: one curve corresponds to the numerical

3 Soln. with Eq. (33) Numerical Soln. Soln. with Eq. (35)

β =4

2

S(X,β)

β =3

β =2

1 β =1 β=0.5

0 0

2

X

4

6

Fig. 4. Growth of the solidified shell with a sinusoidal variation in the heat flux at selected dimensionless time, b, for f = 0.5.

866

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

10 β =30

S(X,β)

8

β =20 β =15

6

β =10

4 β =5

2

Soln. with Eq. (33) Numerical Soln. Soln. with Eq. (35)

0 0

2

β

4

6

Fig. 5. Growth of the solidified shell with a sinusoidal variation in the heat flux at later stages of the process for f = 0.5.

20 Limitting soln. ζ =500 ζ =50 ζ =10 ζ =5 ζ =1

16

12

S0(η)

Numerical solutions

8

4

0 0

4

8

η

12

16

20

Fig. 6. Numerical solution (3), limiting solution (3), and approximate analytical solution using Eq. (30) for the dimensionless zeroth-order solidification front, S0(g) with time, g at selected values of f.

solution (solid) and the other curve corresponds to the approximate solution obtained by using Eq. (30). The limiting solution for f ? 1 developed in Section 7 is also shown for comparison and testing the performance of the numerical algorithm for the large values of f. Increasing leads to faster growth of the mean shell thickness. For large values of f, growth of the mean shell thickness tends to be linear in time. Differences between the numerical and approximate solutions for all selected values of are negligible over the time range considered. Fig. 7 shows the numerical solution for the growth of perturbation in solidification front, S1(g), as a function of mean shell thickness S0(g) at selected values of f. The results for f > 10 exactly reduce to the limiting solution where. This limiting solution is transcribed from the analysis reported in Section 7 for comparison. As it is seen, shell material with higher diffusivity lead to higher magnitude but faster decay in perturbation in solidification. Growth of perturbation S1(g) approaches zero for all cases considered whenever mean shell thickness S0(g) becomes large enough. For all values of considered, the perturbation in shell thickness grows and reaches a maximum, and then approaches to zero whenever S0(g) = 2p (i.e., ms0(t) = 2p). A physical explanation of this behavior is provided by the consideration that when S0(g) = 2p (i.e., ms0(t) = 2p), the thickness of the mean solidified layer is comparable with the wavelength of the perturbation, in which case variations in conditions at the solidification front will have relatively little influence on those at the mold boundary. The perturbation in shell thickness is therefore frozen at zero implying that the shell growth will be planar beyond this point. As f increases, the perturbation in nominal shell thickness grows faster. This implies that the solidified shell materials with higher thermal diffusivities may result in higher probability of occurrence of irregular growth.

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

867

0.2 Limiting Soln. ζ =10 ζ =0.5 ζ =0.3 ζ =0.1

0.16

S1(η)

0.12

0.08

0.04

0

0

2

4

S0(η)

6

8

10

Fig. 7. Perturbation in solidification front, S1(g), as a function of nominal shell thickness, S0(g), for selected values of f.

9. Summary and conclusions An approximate analytical and numerical method for solving a two-dimensional Stefan problem on a planar mold have been developed using a linear perturbation method. A one-dimensional solidification process is perturbed by the superposition of a small spatial sinusoidal variation in cooling rate on the semi-infinite mold surface. Approximate analytical solutions have been obtained for the advance of the solid/melt boundary and the temperature distribution in the shell. The approximate analytical solution for the advance of the solid/melt boundary is compared with a numerical solution, and a very good agreement has been obtained. The effect of thermal diffusivity of the solidified shell material on the development of the shell thickness was investigated. It was shown that solidified shell materials with higher thermal diffusivities may result in irregular growth which, generally, causes cracking near the surface. We have also discussed the limiting solution for zero Stefan number which was useful in the development and checking of the numerical solution of the general problem. The present model developed in this paper can be extended to incorporate the thermal stress field which is coupled with the temperature field along the mold–shell interface through a pressure-dependent thermal contact resistance [29,30]. In this case, changes in the local mold–shell interface pressure will be reflected in the local casting thickness as non-uniformities (i.e., cellular undulations), which is not captured in the physics of the present model. These macroscopic undulations represent growth instability where certain regions of the solidification front grow faster than others. The result is a significant thickness variation in the casting which can have a detrimental effect on the quality of the final cast product and the casting process itself. Non-uniform growth can be accompanied by air gap nucleation which leads to a loss of heat transfer and separation of metallurgical constituents, promotes surface cracking during subsequent process-intensive metalworking operations. The present model does not account for these important phenomena and hence predictions made about the irregular growth of the shell here are clearly restrictive. It is anticipated that the fully coupled version of the present model may provide conclusive demonstration of the growth instability mechanism observed in many casting experiments (see, for example [31–34]).

Appendix. Derivation of Eq. (31) We suppose the governing equation (20) with the mold boundary condition (27) has a solution that is piecewise smooth and absolutely integrable on (0, 1) including its first and second-order derivatives. We consider the Fourier cosine transform of T 1 ðY; bÞ,  Z 1 1 T 1c ðz; bÞ ¼ pffiffiffi T 1 ðY; bÞ cosðzYÞdY: ðA1Þ p 0 Differentiating Eq. (A1) with respect to b under the integral sign and using the differential equation (20), we will have: !1   Z 1 oT 1c 1 oT 1 1 ðz; bÞ ¼ pffiffiffi cosðzYÞ þ zT 1 sinðzYÞ  f2 z2 pffiffiffi T 1 ðY; bÞ cosðzYÞdY : ðA2Þ ob oY p p 0 0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ¼T 1c ðz;bÞ

868

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

ðY;bÞ Under the assumption that T 1 ðY; bÞ and oT 1oY vanish as Y ? 1, and taking note of the boundary conditions, this leads to the non-homogeneous differential equation:

oT 1c 1 Q ðz; bÞ þ f2 z2 T 1c ðz; bÞ ¼  pffiffiffi 1 : ob p Q0

ðA3Þ

The solution of Eq. (A3) that, by initial condition T 1 ðY; 0Þ ¼ 0, satisfies the condition T 1c ðz; 0Þ ¼ 0 is given by  Z b Q 2 2 T 1c ðz; bÞ ¼  pffiffiffi 1 ef z ðbsÞ ds: pQ 0 0 Hence, by taking the inverse transform, we obtain: Z b 1 T 1 ðY; bÞ ¼ pffiffiffi T 1c ðz; bÞ cosðzYÞdz: p 0 Substitution of Eq. (A4) into Eq. (A5) gives: Z 1 Z b Q 2 2 T 1 ðY; bÞ ¼  1 ef z ðbsÞ cosðzYÞds dz: pQ 0 0 0

ðA4Þ

ðA5Þ

ðA6Þ

If we formally interchange the order of integration and use the table of Fourier cosine transform Eq. (A6) simplifies to the form: Q T 1 ðY; bÞ ¼  pffiffiffi 1 pQ 0

Z

b 0

2

ef

ðbsÞY 2 =4ðbsÞ

pffiffiffiffiffiffiffiffiffiffiffi bs

ds:

ðA7Þ

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]

[33]

H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Oxford University Press, London, England, 1956. M.N. Ozisik, Heat Conduction, Wiley, New York, 1980. J. Hill, One-Dimensional Stefan Problem: An Introduction, Longman Scientific and Technical, Essex, 1987. E. Javierre, C. Vuik, F.J. Vermolen, S. van der Zwaag, A comparison of numerical models for one-dimensional Stefan problems, J. Comput. Appl. Math. 192 (2006) 445–459. N. Sadoun, E.K. Si-Ahmed, A new analytical expression of the freezing constant in the Stefan problem with initial superheat, Numer. Meth. Therm. Probl. 9 (1995) 843–854. J. Mennig, M.N. Ozisik, Coupled integral equation approach for solving melting or solidification, Int. J. Heat Mass Transfer 28 (1985) 1481–1485. F. Mosally, A.S. Wood, A. Al-Fhaid, An exponential heat balance integral method, Appl. Math. Comput. 130 (2002) 87–100. N. Sadoun, E.K. Si-Ahmed, P. Colinet, On the integral method for the one-phase Stefan problem with time-dependent boundary conditions, Appl. Math. Modell. 30 (2006) 531–544. S.R. Coriell, G.B. McFadden, R.F. Sekerka, W.J. Boettinger, Multiple similarity for solidification and melting, J. Cryst. Growth 191 (1998) 573–585. J. Crank, R. Gupta, Isotherm migration method in two dimensions, Int. J. Heat Mass Transfer 18 (1975) 1101–1117. C.S. Keung, The use of sources and sinks in solving two-dimensional heat conduction problems with change of phase in arbitrary domains, Ph.D. Dissertation, Columbia University, New York, 1980. J. Caldwell, C-C. Chan, Spherical solidification by the enthalpy method and the heat balance integral method, Appl. Math. Modell. 24 (2000) 45–53. J. Crank, Two methods for the numerical solution of moving boundary problems in diffusion and heat flow, Quart. J. Appl. Mech. X (2) (1957) 220–231. S. Kutluay, Numerical schemes for one-dimensional Stefan-like problems with a forcing term, Appl. Math. Comput. 168 (2005) 1159–1168. J. Caldwell, S. Savovic, Numerical solution of Stefan problem by variable space grid method and boundary immobilisation method, J. Math. Sci. 13 (2002) 67–79. N. Zabaras, S. Mukherjee, An analysis of solidification problems by the boundary element method, Int. J. Numer. Meth. Eng. 24 (1987) 1879–1900. J.S. Hsiao, An efficient algorithm for finite-difference analyses of heat transfer with melting and solidification, Numer. Heat Transfer 8 (1985) 653–666. S. Jana, S. Ray, F. Durst, A numerical method to compute solidification and melting processes, Appl. Math. Modell. 31 (2007) 93–119. J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1984. T.J. Lardner, Biot’s variational principle in heat conduction, AIAA J. 1 (1963) 196–206. A. Prasad, H.C. Agrawal, Biot’s variational principle for a Stefan problem, AIAA J. 10 (1971) 325–327. M. El-Genk, A.W. Cronenberg, Solidification in a semi-infinite region with boundary conditions of the second kind: an exact solution, Lett. Heat Mass Transfer 6 (1979) 321–327. S.H. Cho, J.E. Sunderland, Approximate temperature distribution for phase change of a semi-infinite body, ASME J. Heat Transfer 103 (1981) 401–403. F. Yigit, Perturbation solution for solidification of pure metals on a sinusoidal mold surface, Int. J. Heat Mass Transfer 50 (2007) 2624–2633. S. Savovic, J. Caldwell, Finite difference solution of one-dimensional stefan problem with periodic boundary conditions, Int. J. Heat Mass Transfer 46 (2003) 2911–2916. R.T. Tenchev, J.A. Mackenzie, T.J. Scanlon, M.T. Stickland, Finite element moving mesh analysis of phase change problems with natural convection, Int. J. Heat Fluid Flow 26 (2005) 597–612. W. Rolph III, K.-J. Bathe, An efficient algorithm for analysis of nonlinear heat transfer with phase changes, Int. J. Numer. Meth. Eng. 18 (1982) 119–134. O. Richmond, N.C. Huang, Interface instability during unidirectional solidification of a pure metal, in: Proceedings of the Sixth Canadian Congress of Applied Mechanics, Vancouver, 1977, pp. 453–454. F. Yigit, J.R. Barber, Effect of Stefan number on thermoelastic instabilities in unidirectional solidification, Int. J. Mech. Sci. 36 (8) (1994) 707–723. F. Yigit, Existence of critical wavelength for gap nucleation in solidification on a rigid mold, ASME J. Appl. Mech. 71 (2004) 96–108. H. Murakami, M. Suzuki, T. Kitagawa, S. Miyahara, Control of uneven solidified shell formation of hypo-peritectic carbon steels in continuous casting mold, J. Iron Steel Inst. Jpn. 78 (1992) 105–112. D.A. Weirauch Jr., A. Giron, The early stages of aluminium solidification in the presence of a moving meniscus, in: Proceedings on the Integration of Material, Process and Product Design – A Conference Dedicated to the 70th Birthday of Owen Richmond, A.A. Balkema Publishers, Rotterdam, Netherlands, 1998, pp. 183–191. P.N. Anyalebechi, Undulatory solid shell growth of aluminum alloy 3003 as a function of the wavelength of a grooved mold surface topography, in: P.N. Anyalebechi (Ed.), Materials Processing Fundamentals, TMS (The Minerals, Metals and Materials Society), 2007, pp. 31–47.

F. Yigit / Applied Mathematics and Computation 202 (2008) 857–869

869

[34] P.N. Anyalebechi, Ungrooved mold surface topography effects on cast subsurface microstructure, in: P.N. Anyalebechi (Ed.), Materials Processing Fundamentals, TMS (The Minerals, Metals and Materials Society), 2007, pp. 49–62. [35] F. Yigit, L.G. Hector Jr., Solidification of pure metal with finite thermal capacitance on a sinusoidal mold surface, J. Therm. Stress. 25 (2002) 663–690. [36] T.R. Goodman, The heat-balance integral and its application to problems involving a change of phase, ASME J. Heat Transfer 80 (1958) 335–341. [37] G.W. Evans II, E. Isaacson, J.K.L. MacDonald, Stefan like problems, Quart. Appl. Math. 8 (1950) 312–319. [38] D.A. Anderson, J.C. Tannehill, R.H. Pletcher, Computational Fluid Mechanics and Heat Transfer, Hemisphere Publishing Corporation, New York, 1984.