Applied Mathematics and Computation 148 (2004) 321¨C329 www.elsevier.com/locate/amc

A numerical solution of the Stefan problem with a Neumann-type

boundary condition by enthalpy method

A. Esen, S. Kutluay *

Department of Mathematics, Faculty of Arts and Science, In€on€u University, 44069 Malatya, Turkey

Abstract In this paper, the enthalpy method based on suitable ?nite di?erence approximations

has been applied to the one-dimensional moving boundary problem with a Neumanntype boundary condition known as the Stefan problem. The numerical results obtained by the hopscotch technique are compared with the exact solution of the problem. It is shown that all results are found to be in very good agreement with each other, and the numerical solution displays the expected convergence to the exact one as the mesh size is re?ned. ? 2002 Elsevier Inc. All rights reserved. Keywords: Enthalpy method; Moving boundary; Stefan problem; Hopscotch technique

1. Introduction

Many important heat conduction or di?usion problems arising in science, engineering and industry involve a phase change or a moving boundary. A common example of such problems is the problem of melting of ice that was ?rst treated by Stefan [1] and after whom this kind of problems are widely referred as Stefan problems. Since the moving boundary (interface) separating liquid and solid regions of the material changes with time, its location which is not known a priori and so needs to be determined as a part of the solution.

* Corresponding author. E-mail address: skutluay@inonu.edu.tr (S. Kutluay).

0096-3003/$ - see front matter ? 2002 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(02)00846-9

322

A. Esen, S. Kutluay / Appl. Math. Comput. 148 (2004) 321¨C329

However, such problems are highly non-linear and their exact solutions seem to be inconceivable except in some very simple cases. The exact solutions of the one-dimensional Stefan problems have been surveyed by Hill [2].

In view of the di?culty in obtaining exact solution, numerical methods are far more common. As a result, in recent years, these non-linear problems have motivated considerable research by many authors. To tackle di?cult phase change problems, many numerical methods such as ?nite di?erence, ?nite element and heat balance integral methods [3¨C11] have been and are still being developed. A comprehensive review of a wide range of these methods may be found in [12].

In this paper, an enthalpy formulation based on suitable ?nite di?erence approximations has been applied to one-dimensional Stefan problem with a Neumann condition at the ?xed boundary suggested by Ho?man [13].

2. The governing equation

We consider a one-dimensional Stefan problem in a semi-in?nite half space 0 < v < 1, consisting of a material that can exist in either a liquid or solid phase. The material is initially solid at its fusion melting temperature Tf . At time s ? 0 heat energy is supplied at one end v ? 0 in the form of a constant wall boundary temperature Tw?>Tf ?. The material subsequently melt and a liquid/solid (moving) interface S?s? moves away from v ? 0 and separates the liquid region from the solid one.

In this paper, we are interested in the temperature distribution T ?v; s? in the molten region, 0 < v < S?s? and in the location of the moving interface S?s?. Sometimes we omit the arguments v and s.

The temperature T ?v; s? is governed by the heat conduction equation

qc

oT os

?

K

o2T ov2

;

0 < v < S?s?; s > 0

?1?

subject to the boundary conditions

K

oT ov

?

?h~Tw;

v ? 0; s > 0;

T ?S?s?; s? ? Tf ; v ? S?s?; s > 0

where q is the liquid density (assumed constant), c is the speci?c heat capacity, K is the liquid thermal conductivity (assumed constant), and h~ is the prescribed function of s > 0.

The location of the moving interface is governed by the heat balance

equation known as Stefan condition

A. Esen, S. Kutluay / Appl. Math. Comput. 148 (2004) 321¨C329

323

Lq dS ? ?K oT ; v ? S?s?; s > 0

?2?

ds

ov

where L is the latent heat of fusion. The term on the right hand side of Eq. (1) is called the heat ?ux which represents heat ?ow per unit area per unit time of an isothermal surface (that is, a surface of a constant temperature) in the decreasing temperature direction. Therefore, the heat ?ux is a positive quantity since the heat ?ow is prescribed in the positive v direction, 0 < v < 1.

Initially there is no liquid region which implies the condition

S?0? ? 0; s ? 0

In the above equations, the variables and heat parameters are in terms of physical units. It is convenient to non-dimensionalise the problem by introducing the change of variables

x ? v ; s ? S ; t ? js ; U ? T ? Tf ; j ? K

?3?

¡®

¡®

¡®2

Tw

qc

in which ¡® is a notional length. Under these new variables and h~ ? K exp?js=¡®2?=¡®, a dimensionless model of the problem then becomes as follows:

Heat ?ow in 0 < x < s?t?:

oU ot

?

o2U ox2

;

subject to

0 < x < s?t?; 0 < t < 1

oU ? ? exp?t?; x ? 0; t > 0 ox U ?x; t? ? 0; x ? s?t?; t > 0

Heat balance at x ? s?t?:

ds ? ? 1 oU ; dt a ox subject to

x ? s?t?; t > 0

s?0? ? 0; t ? 0

where a, the latent heat parameter known as Stefan constant, is given be a ? L=cTw and is generally taken as unity in all the numerical algorithms.

Referring to [13,14] this problem has the following exact solution for the temperature distribution U ?x; t? and the interface location s?t?, respectively,

U ?x; t? ? exp?t ? x? ? 1; 0 6 x 6 s?t?; 0 < t < 1 ?4?

s?t? ? t; t P 0

324

A. Esen, S. Kutluay / Appl. Math. Comput. 148 (2004) 321¨C329

The exact solution (4) can be used to compare the numerical solutions and can also be used to initialise the numerical scheme at some t > 0.

3. Enthalpy formulation of the problem

The enthalpy method [15] is a simple and ?exible technique for solving phase

change problems involving either melting or freezing. Instead of working en-

tirely in terms of the temperature of a material, an enthalpy function is de?ned

which represents the total heat content per unit mass of the material. The

advantage of such a reformulation is that the necessity to carefully track the

location of solid¨Cliquid interface is removed and standard numerical tech-

niques can be employed.

As the problem is stated, the non-linearity of the Stefan boundary condition

is the principal cause of di?culty when any analytic or numerical treatment of

the problem are used. To overcome this, the enthalpy formulation for the

problem is considered. The enthalpy formulation for heat di?usion moving

boundary problems is a weak formulation which is similar to that used for the

heat equation for ?nite element methods and which eliminates explicit refer-

ence to the moving boundary [16¨C18]. A function H ?T ? which is the total heat

energy per unit mass, away from the fusion temperature T ? Tf , can be de?ned as follows:

8 ><

RT

0

c dT

;

T < Tf

H

?

>: HR T1

0

6H c dT

6 H2; ? L;

T ? Tf T > Tf

?5?

where

H1

?

R Tf

0

c dT

and

H2

?

R Tf

0

c dT

? L.

Indeed, the material is solid at its melting (fusion) temperature for H ? H1

and the material is liquid at its freezing temperature for H ? H2. When

H1 < H < H2, the material is at its fusion temperature and this region is called

a mushy region in which the material is partly molten and partly solid.

Eq. (5) can be re-arranged to yield

8

< H =c;

H < cTf

T

?

:

Tf ; ?H

? L?=c;

cTf 6 H 6 cTf ? L H > cTf ? L

?6?

Substituting (6) into Eq. (1) and utilizing Eq. (2) reduce Eqs. (1) and (2) to the single equation

q

oH os

?

K

o2T ov2

;

0 < v < S?s?; s > 0

A. Esen, S. Kutluay / Appl. Math. Comput. 148 (2004) 321¨C329

325

which involves two unknown functions, enthalpy H ?T ? and temperature T ?v; s?.

Using the non-dimensional variables in Eq. (6) gives U ? ?H ? cTf ?=cTw which is the temperature distribution in the solid region. In order to distinguish the temperature in the solid region from one in the liquid region we use E in place of U . Thus, we write

E ? H ? cTf ; E < 0

?7?

cTw

Utilising Eq. (7) and using Eqs. (3) and (6) then becomes

&

U?

E ? a; 0;

E>a 06E6a

?8?

Hence the non-dimensional enthalpy formulation of Eq. (7) describing the heat ?ow becomes

oE ot

?

o2U ox2

;

0 < x < ¡®; t > 0

?9?

subject to

oU ? ? exp?t?; x ? 0; t > 0 ox

?10?

U ?¡®; t? ? 0; x ? ¡®; t > 0

where ¡® > s?t?. Thus, standard ?nite di?erence schemes may now be applied to the whole region of the problem. Since we are interested in the temperature distributions in the liquid region, the temperature U ?x; t? can easily be found by virtue of Eq. (8) after computing setting the values of E.

4. Method of solution

The explicit ?nite di?erence scheme of the problem provides a relatively straightforward expression for the determination of the unknowns temperatures at any time step, but it limits the size of the time step.

The implicit di?erence replacement is stable for all sizes of time step, thus there is no size restriction on time step. However, it requires the solution of a system of simultaneous algebraic equations.

The odd¨Ceven Hopscotch algorithm proposed by Gourlay [19] combines explicit and implicit ?nite di?erence replacements at alternate mesh points in such a way that the overall scheme does not require the solution of a system of algebraic equations, that is, it is overall explicit. The method is also unconditionally stable.

326

A. Esen, S. Kutluay / Appl. Math. Comput. 148 (2004) 321¨C329

The solution domain f?x; t? : x 2 ?0; ¡®?; t 2 ?0; 1?g is discretized into cells described by the node set ?xi; ti? in which xi ? ih ?i ? 0; 1; 2; . . . ; N ?tj ? jk, h is a spatial mesh size (which does not vary with time each time step) and k is the time step. The temperature U ?xi; ti? and enthalpy E?xi; ti? are denoted by Uij and Eij, respectively.

The standard explicit odd¨Ceven hopscotch replacement for Eqs. (9) and (10)

is

Eij?1 ? Eij ? 2r?Uij?1 ? Uij ? h exp?tj??; i ? 0 Eij?1 ? Eij ? r?Uij?1 ? 2Uij ? Uij?1?; i ? 1; 2; . . . ; N ? 1

?11?

for j ? 0; 1; 2; . . . J , where r ? k=h2. Eqs. (11) are used at nodes i ? j even. At

nodes ?xi; tj? when i ? j is odd, the following standard implicit ?nite di?erence approximation of Eqs. (9) and (10) is used. A value Eeij?1 is calculated from

Ee

j?1 i

Ee

j?1 i

? ?

Eij Eij

? ?

2r?rU?Uij??ij1??111??Uhij??e11x?p;?tj?i1???;1;

i 2;

? ..

0 .;

N

?

1

?12?

for j ? 0; 1; 2; . . . ; J , where the temperatures on the right hand side of Eqs. (12)

are known from

should

Ee

j?1 i

>

a

Eqs. (11). it may be

If Eeij?1 6 deduced

a then Eij?1 6 that Eij?1 > a

a and so Uij?1 ? and therefore,

0.

However,

Eij?1 ? ?Eeij?1 ? 2ra?=?1 ? 2r?

with corresponding temperature

Uij?1 ? Eij?1 ? a:

5. Numerical result and discussion

All numerical calculations were performed in double precision arithmetic on a Pentium III processor under DOS 6.0 and using prospero FORTRAN compiler.

In order to avoid discontinuities in the initial data at t ? 0:0 when using the hopscotch technique, the numerical scheme was started at t ? 0:1 using the exact solution (4) and all results were obtained at a ?nal time tf ? 0:9 by taking a ? 1:0. The starting values for E?U ? inside the liquid region are determined from Eq. (8). It is also possible to estimate the location and speed of the moving interface by the enthalpy E?U ?.

In order to display how good the numerical solution of the problem we shall use the weighted 1-norm kek1 de?ned by

A. Esen, S. Kutluay / Appl. Math. Comput. 148 (2004) 321¨C329

327

kek1

?

1 N

X N ?1

i?1

1

?

U

Uij ?xi;

tj

?

;

e ? ?e1; e2; . . . ; eN?1?T

Tables 1 and 2 present, respectively, numerical results for the temperature distribution with the norm kek1 and the interface location for r ? 0:4. It is observed that all the results show convergence to the exact solution as the mesh size h is re?ned. The prediction of the interface location is very good. The percentage error decreasing from 0.0411% (for h ? 0:1) to 0.0001% (for h ? 0:01).

Tables 3 and 4 display, respectively, numerical results for the temperature distribution with the norm kek1 and the interface location for r ? 2:0. It may be seen from the results presented in Table 3 that the numerical scheme converges to the exact solution as the ?nite di?erence mesh h is re?ned. The estimation of the interface location is also very good whose percentage error decreases from 0.1056% (for h ? 0:1) to 0.0003% (for h ? 0:01). The purely explicit method becomes fail for r ? 2:0 since it is stable for r 6 0:5. It should also be noted that

Table 1 Values of the temperature distribution as predicted by the numerical and exact solutions at tf ? 0:9 and various mesh sizes for r ? 0:4

Explicit

h ? 0:1

x

1.50022

0.0

1.26667

0.1

1.05554

0.2

0.86450

0.3

0.69127

0.4

0.53341

0.5

0.38833

0.6

0.25320

0.7

0.12486

0.8

0.00926

0.9

Hopscotch

h ? 0:1

1.49909 1.26552 1.05438 0.86336 0.69014 0.53239 0.38746 0.25258 0.12658 0.00925

h ? 0:05

1.48121 1.24725 1.03559 0.84415 0.67101 0.51424 0.37172 0.24081 0.11822 0.00874

h ? 0:01

1.46407 1.23004 1.01828 0.82667 0.65329 0.49641 0.35448 0.22608 0.10935 0.00108

Exact solution

1.45960 1.22554 1.01375 0.82212 0.64872 0.49182 0.34986 0.22140 0.10517 0.0

kek1

0.08376

0.05139

0.01444

Table 2 Predicted location of the interface at tf ? 0:9 for r ? 0:4

Hopscotch

h

s

Error (%)

0.1

0.9003703

0.0411

0.05

0.9000937

0.0104

0.01

0.8999822

0.0001

Exact solution 0.9

328

A. Esen, S. Kutluay / Appl. Math. Comput. 148 (2004) 321¨C329

Table 3 Values of the temperature distribution as predicted by the numerical and exact solutions at tf ? 0:9 and various mesh sizes for r ? 2:0

Explicit

h ? 0:1

x

Method unstable 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Hopscotch

h ? 0:1

1.50235 1.26642 1.05157 0.85846 0.69746 0.55766 0.40828 0.26326 0.14908 0.04756

h ? 0:05

1.48867 1.25465 1.04185 0.84934 0.67776 0.52343 0.37981 0.24727 0.12183 0.04277

h ? 0:01

1.46450 1.23045 1.01868 0.82707 0.65369 0.49681 0.35487 0.22651 0.10969 0.00108

Exact solution

1.45960 1.22554 1.01375 0.82212 0.64872 0.49182 0.34986 0.22140 0.10517 0.0

kek1

0.11545

0.06128

0.01545

Table 4 Predicted location of the interface at tf ? 0:9 for r ? 2:0

Hopscotch

h

s

Error (%)

0.1

0.9095119

0.1056

0.05

0.9002139

0.0237

0.01

0.8999707

0.0003

Exact solution 0.9

the number of iteration corresponding to the same mesh sizes decreases when the stability parameter r increases.

In conclusion, we have demonstrated that the enthalpy approach can be used e?ectively in problems with a Neumann condition at the ?xed boundary and that the hopscotch scheme is an e?cient solver for such problems since hopscotch di?erence approximation produces the expected converging numerical solutions for r > 0.

References

[1] J. Stefan, Uber die theorie der eisbildung imbesondee uber die eisbindung im polarmeere, Ann. Phys. U. Chem. 42 (1891) 269¨C286.

[2] J. Hill, One-dimensional Stefan Problem: An Introduction, Longman Scienti?c and Technical, Essex, 1987.

[3] J. Crank, Two methods for the numerical solution of moving boundary problems in di?usion and heat ?ow, J. Mech. Appl. Math. 10 (1957) 220¨C231.

[4] J. Crank, R.S. Gupta, A moving boundary problem arising from the di?usion in absorbing tissue, J. Inst. Math. Appl. 10 (1972) 19¨C33.

A. Esen, S. Kutluay / Appl. Math. Comput. 148 (2004) 321¨C329

329

[5] S. Kutluay, A.R. Bahadir, A. O€ zdes, The numerical solution of one-phase classical Stefan problem, J. Comp. Appl. Math. 81 (1997) 35¨C44.

[6] N.S. Asaithambi, A Galerkin method for Stefan problems, Appl. Math. Comput. 52 (1992) 239¨C250.

[7] N.S. Asaithambi, A variable time step Galerkin method for a one-dimensional Stefan problem, Appl. Math. Comput. 81 (1997) 189¨C200.

[8] W.D. Finn, E. Varoglu, Finite element solution of the Stefan problem, in: J.R. Whiteman (Ed.), The Mathematics of Finite Elements and Applications, MAFELAP 1978, Academic Press, New York, 1979.

[9] T.R. Goodman, The heat balance integral and its application to problems involving a change of phase, Trans. ASME 80 (1958) 335¨C342.

[10] G.E. Bell, A re?nement of the heat balance integral method applied to a melting problem, Int. J. Heat Mass Transfer 21 (1978) 1357¨C1362.

[11] A.S. Wood, A new look at the heat balance integral method, Appl. Math. Modelling 25 (2001) 815¨C824.

[12] J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1984. [13] K.H. Ho?man, Fachreich Mathematik, Berlin Freie Universitat, 1977. [14] R.M. Furzeland, A comparative study of numerical methods for moving boundary problems,

J. Inst. Math. Appl. 26 (1980) 411¨C429. [15] M.E. Rose, A method for calculating solutions of parabolic equations with a free boundary,

Math. Comput. 14 (1960) 249¨C256. [16] D.R. Atthey, A ?nite di?erence scheme for melting problems, J. Inst. Maths. Appl. 13 (1974)

353¨C366. [17] V.R. Voller, M. Cross, Accurate solutions of moving boundary problems using the enthalpy

method, Int. J. Heat Mass Transfer 24 (1981) 545¨C556. [18] A.S. Wood, S.I.M. Ritchie, G.E. Bell, An e?cient implementation of the enthalpy method, Int.

J. Numer. Meth. Engng. 17 (1981) 301¨C305. [19] A.R. Gourlay, Hopscotch: A fast second-order partial di?erential equation solver, J. Inst.

Math. Appl. 6 (1970) 375¨C390.