Numerical modeling of Stefan problems

Marhl Cornput. Modelling, Vol. 14, pp. 139-144, 1990 Printed in tireat Bntain

NUMERICAL

08957177/90 $3.00 + 0.00 Pergamon Press plc

MODELING

OF STEFAN

PROBLEMS

N. S. Asaithambi

Department Mississippi

of Mathematics State

and Statistic3

University

Abstract. A numerical method for the solution of two-dimensional Stefan problems will be presented in this paper. The major features of the method are the use of boundary-fitted coordinates, and an implicit marching procedure. Illustrative examples will be given in order to demonstrate the accuracy of the method. Numerical results of the present method will be compared with those obtained previously.

Kevwords: Stefan problems; tridiagonal

moving bounda.ry;

INTRODUCTION

numerical treatments can be found in Douglas and Gallie (1955), Murray and Landis (1959), Gupta (1972, 1974), Crank and Gupta (1972), Gupta and Kumar (1980, 1981, 1984), and Asaithambi (1988). In this paper, a finitedifference method with coordinate transformation is presented. The method begins by mapping the changing physical domain to a fixed rectangular computational domain by means of a simple coordinate transformation. This approach is similar to that of Asaithambi (1987, 1988) and Gupta and Kumar (1984). The governing equation and the boundary conditions are then transformed

Certain physical problems in heat conduction can be modeled as ‘moving boundary problems’ (MBPs), often referred to as Stefan problems. Various analytical and numerical approaches for obtaining approximate solutions to Stefa.n problems have been explored. Examples of analytical methods can be found in Goodman (1964), Rasmussen (3.977), Solomon (1978), and Gupta and Kumar (1986). Approximate analytical solutions are somewhat simpler to obtain than numerical solutions. Bowever, there is no check on the accuracy of these analytical approximations. Examples of ncn 14-F

finite-differences;

systems.

139

Proc. 7rh Int. Co@ on Mathematical and Computer Modelling

140

accordingly.

An implicit

scheme

is used

to discretize the governing equations, and a marching, procedure is used to advance the solution from one time-step to the next. The accuracy of the method is illustrated by solving numerically a test problem whose exact solution is known. The method is then applied to the examples of Gupta and Kumar (198G). STATEMENT

OF

+ uYY

an initial

y = 0 at t = 0. There

interface

responding

initial

temperature

tion in the domain

on D,

(6) and D.

t > 0,

(I)

/p$ 3

0 II

II water

s”

at

2 = O,l,

(2)

u=l,

at

y=o,

(3)

u = 0,

at

y = s(z,t).

(4)

At the moving interface, an additional condition known as the Stefan condition must also hold. The Stefan condition is described by y = s(z,t).

(5)

In (2)-(5), subscripts are used to denote partial differentiation. (e.g., u1 = du/dt etc.). Usually, the melting does not start

$

_ ST=1

Figure 1 A Two-Dimensional Stefan Problem

and y = s(z, t) is the moving interface. (See Figure 1). The following boundary conditions hold:

at

0).

ice

Y = 0

St = -uy + s,u,,

by

The problem is to determine s(z,t) u(z:,y,t) for all t > 0 in the domain

0

0,

distribu-

D described

u(5, y, 0) = 1 - y/s(z,

where

U t- -

is

y = s(z, 0) and a cor-

PROBLEM

Let ((2, Y)~(z,Y) E [O, 11x io,~1) be a semi-infinite slab of ice. If the temperature at the bottom surface y = 0 is raised, the ice starts melting. An interface between water and ice forms, and starts moving up as time progresses. Letting u(z, y, t) denote the temperature at (z,y) at time t, we can describe this heat conduction process by the heat equation in its dimensionless form as follows: uI = u,,

at the surface

COORDINATE TRANSFORMATION A more general moving boundary problem is obtained by replacing the conditions (3)-(6) with at y = 0, (3a)

u = 9,

at y = 5, (4a)

u = h, s1 = -uy U=l where f(z,t),

-+ s,u,

+ f at y = s, (5a) at t = 0, (Ga)

g = g(x,t), h = h(+), f = and r = r(z,y) can be viewed as

Proc.

7th Int. ConJ on Mathematical

some forcing functions, and ~(2, t) has been abbreviated to s. Under the transformation

141

Mode&g

with At = 1/(1-l), A7 = l/(J-l), and AT chosen to satisfy appropriate accuracy criteria. The Stefan condition (12) is used to advance the moving boundary s in time, while the heat equation (8) is

E =5 17 = Y/Sk

t)

(7)

r=t the heat equation

and Computer

(1) becomes

UT = UE~ - 2Au6,

+ Bu,,,, + Cu,,,

(8)

on the domain

D’ = w?)l(l> 17)E 10,11x 10,111, where

used to advance the termperature. We will use the subscripts i, j and the superscript 72 to denote the dependence of a function on E, r), and r respectively. For instance, Ulj M U([i,qj,Tn), and S? M s([i, 7,). The computational problem is then to determine .sr and u~j for all n 2 1 given sf and up,j. A marching procedure can be used to advance the solution from sn and u!‘. to srfl respectiveli The dkretization and (8) leads to systems

A = vcls, B = (1 + v2s;)/s2, c = rl (+/s

+ (s&)2

The condition UC =

(2) gets transformed

96

sr

=

to (9)

(3a) and (4a) become

4,

u = h(t,r), The Stefan forms to

*

at t = 0,l

Au,,,

and the conditions

u=

- (SE/S)0

condition

-pu, + spt

at 7j = 0,

(IO)

at 77= 1.

(11)

at 77 =

1 trans-

+ f(b),

(12)

ci = (i - l)A[, =

1 - (j - l)Aq,

rn = nor,

al-

Note that the equation (12) is in the form s, = F(s, SE, 21,us, uV) for some F. Discretization along 7 using a backward difference yields s n+l

=

Sn

+

&J’“+’

(13)

METHOD

A rectangular grid is imposed on the computational domain so that a grid point is described by (E;, qj, 7,) where

qj

of nonlinear

gebraic equations in the unknowns ST+’ and u:;‘. Hence, the solution method turns out to be iterative. Starting with an initial (‘guess” for the position sl+r of the moving boundary for each [’ at r,+i , equation (8) could be solved to obtain u in the interior. With the approximation to u in the interior, the Stefan condition is used to obtain the next iterate for .sn+l L +

where F”+l denotes evaluating F at rn+r. Spatial discretization (along 0 results in a system of the form

where P = (1 + $)/s.

NUMERICAL

and usj” of (i2)

i = 1,2,-e.

,I,

j = 1,2,.-a

,J,

n = O,l.,...

,

for i = 2,3,... , I - 1, when a centered difference is used for the spatial derivatives. Conditions (9) and (12) are combined to yield

Proc. 7th Int. Cony. on Mathematical and Computer Modelling

142

at [ = 0,l

so that

their

discretization

at the left and right ends also results in equations similar to (14). Let ~y+l,(~) denote the k ‘I’ iterate for s.r-e’. Newton’s method can be used to obtain ‘i

n+r,(k+I)

=

S~+LG9

+

Asi

= -‘Hi,

proceeds

ai

(17) PI

1, where

by

U~+l,(k+l) = u[J”~(~) + Aui,j (21) I,i in which Aui,j is obtained by solving PiAui-r,j

+ oiAui,j

+ riAui+r,j

=

=

(l +

AT/(U)“,

z/J1 + ~2

=

= 2;..

= Ki,j (22) ,J-1,

2Bi,jP2), A~/(ATJ)~,

Ki,j

=

and w is the SOR acceleration parameter. The boundary conditions at t = 0,l and n = 0,l will supply the necessary equations corresponding to i, j = -wSi,j,

ai = 1 - (~tl)l”P,“‘A~/sl”) bi = jjX(U[)lsl

scheme

for i = 2,. . . ,I-1 andj where pi = yi = -/or,

biAsi-l+aiAsi+ciAsi+* ,I--

the line-SOR

setting

(16)

by solving for Asi in

for i = 2,3,...

Then

- x(s~)~+‘(u~)~+l/s~+‘,

1, i = I, and j = J. Once again, the linear systems to be solved for each j are

ci = -bi, and X = AT/A<. Corresponding 1 and i = I, we have from (Is), al = 1 - Ar(~,,)~+‘/(s~+~)~, aI = 1 - Ar(~,#‘/(s;+‘)~.

to i =

(18)

With cr = by = 0, the 1inea.r system (17) can be extended to i = 1 and i = I as well. Note that the resulting system of linear equations is a kidiagonal system, and can be solved for Asi very efficiently. The discretization of (8) yields a linear system of equations for the ui,j when the coefficients A, B, and C are evaluated using the most recent iterates available for the si. An iterative technique such as the line successive overrela.xation, (or simply, line-S OR) can be used to solve this linear system. Suppose (8) is expressed as Cu = 0. Let Cau

= 0

(19)

be the corresponding discretization using difference formulas that are second order accurate in space. Unless the u has converged, Cau will not be zero. So, let (C*U);,j

=

Si,j.

(20)

tridiagonal, and hence very efficiently.

can be handled

Now, our numerical method can be described by the following algorithm: Algorithm. [To solve a Stefan in two dimensions]

problem

(1) [initialize] (2)

(3)

72 t 0. r +-- 0. initialize s and u for T,,. [start loop] k t-- 0. “guess” u and s for 7,$1. [boundary] Solve (17) to obtain sn+l,(k+l)

(4) [interior] (5) it;

Perform one sweep of lineSOR on (19). [exit?] Exit if happy; else k t k + 1; go to step (3). [advance] n c n + 1; T c 7,. [stop?] Stop if r, > rmax; else go to step (2). I

Note that the boundary and the interior iterations have been coupled in steps (3) and (4) so that a complete solution of u in the interior is not necessary to obtain the next iterate for s. The stopping criterion used in place of the “if happy” condition in step (5) above was that the relative change in the iterates be less than a prescribed tolerance E.

143

Proc. 7th Inr. Co@ on Mathematical and Computer Modelring

NUMERICAL A test problem

RESULTS

was generated

with exact

solution s(Q)

= (1+

cosz)cost,

U(2, y, t) = &--n2)t--y which satisfies

(1) and The functions g, h, f, propriately chosen so (3a)-(6a) as well. The

AC! 0.1 0.05 0.025

were superior

mussen

(1977).

in s

error( 10M4) 4.11 1.03 0.52

to those

II. Interface

of Ras-

at z = 0

(23)

(2) automatically. and r can be apthat (23) satisfies present numerical

I. Errors

results

TABLE

cos 7r2,

method was tested on this problem for various grid spacings. It was observed that the method is consistent and second order accurate in space. The error in s(z, t) at a fixed time t was measured using the /a-norm. The following errors were recorded corresponding to t = 0.05, with AT = 0.001 and I = J = 11,21,41. TABLE

Gupta and Kumar (1986) have not provided any evidence showing why their

_

Suppose the errors in s are 0 ((A[)p). Let ei denote the error in s corresponding to Ati. Let A&, A&, and A& be such that A&/A& = A[s/A& = V. Then the order of accuracy p can be estimated using

From the results of Table I which correspond to v = 2, p M 2.5 indicating second order a.ccuracy. Next, the numerical method was tested on the physical problem described by (l)-(6), starting with initial profiles provided by Gupta and Kumar (1986). It is clear from Table II that the present method reproduces their results thus validating their method. However,

t

present

0.05 0.10 0.15 0.20 0.25

0.4147 0.5020 0.5761 0.6417 0.7011

0.30 0.35 0.40

0.7557 0.8065 0.8542

previous 0.4144 0.5016 0.5756 0.6410 0.7002 0.7547 0.8055 0.8532

Table II shows the position of the interface at 2 = 0 at various times as computed by the present method and the previous (Method I of Gupta and Kumar (1986)) method corresponding to the profile s(2,O)

= 0.5 - 0.2 cos +71x.

The results of the present method agree with those produced by the Method I of Gupta and Kumar (1986) for all the other initial profiles assumed as well. CONCLUSIONS A finite-difference marching method for the solution of two-dimensional Stefan problems has been presented. It employs boundary-fitted coordinates and a backward difference scheme for the numerical solution. The consistency and the order of accuracy of the method were demonstrated by testing its performance on a test problem with a known exact solution. The method reproduced the results on a set of problems which were reported earlier in the literature by Gupta and Kumar (1986).

Proc. 7th Int. Conf on Mathematical and Computer Modelling

Gupta, R. S., and D. Kumar (1984). Variable time-step method with co-

REFERENCES Asaithambi, N. S. (1987). Computation of free-Surface flows, J. Comput. Phys. 73, 380-394.

ordinate Meths. 91-103.

transformation, Appl.

Comput.

mech.

Engrg.

44,

Asaithambi, N. S. (1988). On A variable time-step method for the one dimenCompu t. sional Stefan problem, Meths. Appl. Mech. Engrg. 7l, 1-13.

1986). ApGupta, R. S., and A. Kumar proximate analytical so‘I utions for multi-dimensional Stefan problems, Comput. Meths. Appl. Mech. Engrg. 56, 127-138.

Crank, J., and R. S. Gupta (1972). A method for solving moving bound-

Murray, W. D., and F. Landis (1959). Numerical and machine solutions of

ary problems in heat flow using cubic splines or polynomials, J. Inst.

involving

Math.

Heat Transfer 8l, 106-112.

Appl.

a,

296-304.

Dou las, J., and T. M. Gallie (1955). On tB e numerical integration of a parabolic differential equation subject to a moving boundary condition, Duke Math. J. 22, 557-570. Goodman, T. R. (1964). Application of integral methods to transient nonlinear heat transfer, in: T. F. Irvine, eds., AdJr., and J. P. Hartnett, vances in Heat Transfer, Academic Press 1, 51-122. Gupta, R. S. (1972). Some numerical and integral methods for solving moving boundary problems in heat flow and diffusion, Ph. D. Thesis, Brunei University, Uxbridge, U. K. Gupta, R. S. (1974). Moving grid method without interpolations, Comput.

Meths. 152.

Appl.

Mech.

Engrg.

4, 143-

Gupta, R. S., and D. Kumar (1980). A modified variable time-step method for the one-dimensional Stefan problem, Comput. Meths. Appf. Mech. Engrg. 23, 101-109. Gupta, R. S., and D. Kumsr (1981). Variable step methods for one dimensional Stefan problems with mixed boundary conditions, Internat. J. Hea.t Transfer 24, 251-259.

transient

heat

conduction

melting

or

problems

freezing,

J.

Rasmussen H. (1977). An a proximate method for solving two- 1 imensional Stefan problems, Lett. Heat Mass Tra.nsfer 4, 273-277. al, eds. Solomon, A. D. et. Moving Boundary Problems, mic Press, New York.

(1978). A cade-