Approximate analytical methods for multi-dimensional stefan problems

COMPUTER METHODS NORTH-HOLLAND

IN APPLIED

MECHANICS

AND ENGINEERING

56 (1986) 127-138

APPROXIMATE ANALYTICAL METHODS FOR MULTI-DIMENSIONAL STEFAN PROBLEMS Radhey S. GUPTA

and Ambreesh

Department of Mathematics, University of Roorkee,

KUMAR Roorkee 247672, India

Received 23 July 1984

Two approximate analytical methods are presented for solving two-dimensional Stefan problems. One of them selects one-dimensional quadratic temperature profiles in one variable at the preselected values of the other variable. The second one is based on piecewise linear profiles. The first method is extended for solving a three-dimensional problem also. The numerical results for some selected problems are obtained by using the present methods as well as some methods due to other authors. The merits of the proposed methods are discussed.

1. Introduction There exist a great many transient problems in heat conduction and diffusion whose domain of reference changes its shape and size with time. Such problems are most commonly known as ‘moving-boundary problems’ (MBPs). Since contours of the moving boundary have to be determined along with the solution, these problems are not amenable to analytical methods. The most common approach for solving MBPs, though not trivial, is through finite difference methods. Examples of some of these methods may be found in [l-5]. The other line of thought is to look for some kinds of approximate analytical solutions which are much easier to manipulate. The chief exponent of one of these methods is Goodman [6] who introduced his ‘integral method’ for solving heat conduction (diffusion) problems including the moving boundaries. Basically, the method assumes a global temperature profile, comprising unknown (time-dependent) parameters which are determined by satisfying the prescribed boundary conditions. One parameter, defined as ‘penetration distance’, is determined by equating the total amount of heat at any time t under the assumed profile to that of the true solution. In spite of the handicap that there is no check on the accuracy of the solution and, also, that there is no foolproof device to improve the accuracy, the approximate analytical methods are still much sought for, because of their simple forms. Hpwever, in order to achieve a better accuracy in one-dimensional MBPs, Bell [7,8] uses piecewise linear functions instead of a global temperature profile. This approach, which scans a smaller subdomain, is likely to give better results. In the higher dimensions, Poots [9] considers freezing of a square prism and assumes a two-dimensional surface profile. Riley and Duck [lo], while considering solidification of a cuboid, use a three-dimensional surface profile. Although important information may be extracted from these solutions, specially for small times t, the simplicity of the ‘integral method’ is unfortunately lost. The expressions for the two- or three-dimensional profiles as obtained in [9, lo] are very cumbersome to work with. To make matters simpler, Rasmussen [ll] solves a 00457825/86/$3.50

@ 1986, Elsevier Science Publishers

B.V (North-Holland)

128

R.S. Gupta, A. Kumar, Approximate analysis for Stefan problems

two-dimensional problem choosing a set of one-dimensional profiles instead of one single two-dimensional global profile. For example, referring to the problem cited in [ll], he fits straight lines in the y-direction at preselected points on the x-axis (details are given in Section 3). But the linear profiles, as expected, can hardly suit to any practical problem. In the present paper we suggest a method for selecting quadratic profiles which can cover a much wider class of problems. An aIternative procedure is also developed for using linear profiles in a piecewise manner. It is shown for already tested problems that the methods suggested in the present paper work very well when the other methods cannot be employed or give erratic results.

2. The ~wo-d~e~io~a~

problem

If u(x, y, t) denotes the temperature at a point (x, y) in a two-dimensional space domain D at any time t, the heat conduction equation in its nondimensional farm may be written as uI=u,,+uYY

inL),

t>O,

(2-I)

where ut, u,,, and uYy are the partial derivatives of u. Let us assume that the domain D is bounded by the fixed boundaries x = 0, x = 1, and y = 0 and the moving boundary y = s(x, t), Further, let D={x,

y]O=%Gl,

O~y~s(x,t)}

be occupied by water with the following boundary

conditions:

Ux=O,

x=0

U==l,

y=O,

(2.3)

u=o,

y = s(x, t) .

(2.42

and x=1,

(2.2)

The portion D’ = {x, y IOs x c 1, y 3 s(x, t)) is assumed to be occupied by ice at 0 “C. It may be noted that the size of D increases with time while L)’ dimi~shes (see Fig. 1). an extra condition emerges at the ice/water (in general In addition to (2.2)-(2.4), solid/liquid) interface due to latent heat and given by s, = -uy + u,s, ,

y =

s(x, 0 ,

(2.5)

where s, and s, are partial derivatives of s. If the melting were to start from the face y = 0 at zero time, the problem would have been straightforward. Instead, there is a condition imposed that the initial temperature distribution is given by

u(x, Y> = 1 -

Yk(X)

7

where g(x) is an arbitrary function such that s(x, 0) = g(x).

P-6)

R.S. Gupta, A. Kumar, Approximate analysis for Stefan problems

129

ICE

WATER

Fig. 1. Moving interface at any time t.

We wish to find the solution of (2.1) in the time-varying domain D subject to conditions (2.2)-(2.6) which obviously calls for the determination of the position of the interface y = s(x, t) in order to define D.

3. The integral method with one-Dimensions profiles According to the conventional ‘integral method’ of Goodman [6] one shall have to choose a two-dimensional temperature profile U(X, y, t) which satisfies the prescribed boundary conditions. However, the selection of such a profile requires some intuition and the resulting form, as stated in Section 1, may be very complex to deal with. Instead of choosing a ho-dimensional global profile, Rasmussen [ll] proposes a finite number of one-dimensional linear profiles, one for each prescribed value of x. Let us illustrate this method in brief for the present problem. The linear profile satisfying conditions (2.3) and (2.4) may be written as 24 =

1

- y/s(x, t) ,

for a particular x .

(3.1)

130

R.S. Gupta, A. Kumar, Approximate analysis for Stefan problems

In order to determine s, (2.1) is integrated with respect to y from 0 to s (in accordance with the ‘integral method’, see [ll]). After using (2.6) and approximating u by (3.1), the following equation is obtained: s, =

4 - s,, + 3 * (1 /s)

)

(3.2)

from which the position of the interface can be determined step-by-step manner.

at any time t for a given x in a

4. The present methods As the choice of a linear profile is a very crude approximation for the temperature distribution in the y-direction at all times, the method suggested in the preceding section may give rise to very crude results. We suggest here two methods which are likely to give better results-one assumes quadratic profiles and the other piecewise linear functions. They are discussed in detail in Sections 4.1 and 4.2. 4.1. Method I: Quadratic profile Let us assume that the temperature

profile for constant x, satisfying (2.3) and (2.4), is given

by

u=(l-f)+d(1-;),

(4.1)

where (Yis an unknown parameter and s, as before, is the position of the solid/liquid interface (for a given x). For determining cr we derive an extra condition by differentiating U(X, y, t) = 0 on the interface y = s(x, t). For x = constant, it gives UYS,+ u, = 0 )

(4.2)

which, after using the basic equation

(2.1) and the condition

uY(u,s, - uY) + u,, + uYY= 0 Further,

substitution (r2(1+sZ)+

at the interface (2.9,

on y = s(x, t) .

becomes (4.3)

of u from (4.1) into (4.3) results in a quadratic equation (Y(4--s’s,,

+ 6s;) + (1 - s - s,, + 3s;) = 0 ,

(4.4)

the solution of which provides the value of (Y. It may be noted that cx is a function of t, which should satisfy (Y(O)= 0 in conformity with (2.6). In order to establish a relationship between s and t, instead of using the heat balance integral, we match the condition at the interface which, after employing the assumed profile, gives

131

R.S. Gupta, A. Kumar, Approximate analysis for Stefan problems

SZ)/S .

s, = (1+ a)(1 +

(4.9

Equation (4.5) is solved numerically interface s at different times.

in a step-by-step

manner

to give the position

of the

4.2. Method II: Piecewise linear profile We subdivide the temperature temperatures as ui=l-j/N,

range [0, l] into N equal intervals,

and denote

j=O(l)N,

(4.6)

where u,, = 1 and uN = 0 are the extreme temperatures. be subdivided into M equal intervals, such that

Further, let the space interval 0 s x s 1

i = O(l)M .

xi = ilM ,

(4.7)

The position of each isotherm uj at a given xi is denoted by yi, j* By integrating (2.1) over each subregion [y,, j, yi, j+l]’ j=O(l)N YI,,+I u, dy =

I Yt. ,

various

Yr,,+1 u,, dy + [u,]~:~;+’,

I Yi, ,

i = O(l)M

- 1, we get

(4.8)

.

This can also be written as YL,+I

I

d

at

u

Yr. ,

=-

dy - Uj+l_Vi,j+l + Uj)ii, j

d2 Yr, dX2 Yh

I

+

(U,)i,

,+l

dy - (U,)i, j+IY:l, j+l + (U,)i, jY:, j -,Uj+lY’i, j+l + ‘jy:I j

u

I

j+l

-

(‘Y)i,

(4.9)

j ’

where j, is the partial derivative with respect to t and y’ and y” with respect to x. Approximating u by a linear profile

u=uj+

‘i+l Yi,j+l

uj

_y,

z,,

(4.10)

(Y-Yi,j)

in each region [yi, j’ yi, j+l]’ i = O(l)M, and substituting partial differential equations: b(U, - u,)Jji, 1 = 4t”fJ - ul)Y: lC”j

- uj+*)(Yi,

it in (4.9), gives the following set of

1 - C”_r>i, IYl, 1 + C”y>i, 1 - (‘y>i,O

j+l + 3i, j) = lC”j

- ‘j+l)(Y:(,

+ C”x)i, jY;

j +

j+l + Yl j) - (‘x)i,

(‘y>i,

j+l

-

(u,)~,

(4.11)

7

j,

j+lYz’, j+l

j-= l(l)N

- 2, (4.12)

132

R.S. Gupta, A. Kumar, Approximate analysis for Stefan problems

(1 + 4%1)&v + ~~N-lL+,

= au,-,(Yy,,

+ yl:N-1)+ (“x>i,N-lYi’, N-1

-

C”y>i,

N-1

for’ = ‘C1jM *

3

(4.13)

It may be noted that conditions (2.3)-(2.5) h ave been used in deriving (4.11)-(4.13). The above set of partial differential equations can be solved numerically provided the x-derivative of U, i.e. u,, is known. In the present case, this is done by the linear interpolation as follows. Let us denote by ui, i the temperature at x = xi_ 1 corresponding to the y-value at the point (xi, uj) and by u,: j at x = xi+, for the same value of y. Then we can write

(‘x)i,

(4.14)

j =

where u_ j and u+ j are computed -

‘j+l

u,;j=uj+

Yi+l, j+l

from the following interpolation

formulae:

uj

- Yi+l.

j

(Yi, j-Yi+l,j)

and Uj Ui, j =

uj-l

‘j-1

+ Yi-1,

j-

Yi-1,

j-l

(Yi, j-Yi-1.

j-1).

It may be noted that the y-derivative of u, i.e. u,, is approximated of gradients of two adjacent profiles, i.e.,

by taking the average

(4.15)

(4.16) (4.17) The solution of the above set of PDEs provides the movement of various isotherms with time including the interface ( yi. ,,, in the present problem).

5. Numerical

results and discussion

In order to make a comparative study, we take the same initial position of the interface as chosen by Rasmussen [ll]. He considers three cases, viz. : (i) (ii)

g(x) = 1 - 0.2e-25”2 , g(x) = 0.5 - 0.2 cos $nx ,

133

R.S. Gupta, A. Kumar, Approximate analysis for Stefan problems

(iii)

g(x)=O.5-0.2cos7Fx,

OSxGl.

In both methods we have employed forward difference in the time direction and central difference in the space direction. The time and space steps are chosen to be At = 0.001 and Ax = 0.1, respectively, in both methods. In Method II, the temperature range is subdivided into 5 equal intervals, which implies that Au = 0.2. The results obtained from Methods I and II are given in Table 1. As can be seen they are quite close to each other. The corresponding values from Rasmussen [ll] are also given. Although we are unable to prove that our results are better than those of Rasmussen, we feel that they must be so, because there is a good agreement between the results from Method I and Method II (see also Fig. 2). Let us consider yet another initial interface profile given by g(x) = 2 + cos ITX)

0c x d

1)

which proves that the method of Rasmussen gives inaccurate results while Method I works very Table 1 Positions of the interface along x = 0 an9 x = 1 at various times. (A) S(X, 0) = 0.5 - 0.2 cos fnx; 0.5 - 0.2 cos 7rx; (C) S(X, 0) = 1 - 0.2e-25” x=1

x=0 t

(B) s(x, 0) =

Rasmussen [ 1l]

Method I

Method II

Rasmussen [ 1l]

Method I

Method II

(A) 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40

0.4022 0.4832 0.5523 0.6135 0.6689 0.7197 0.7668 0.8110

0.4144 0.5016 0.5756 0.6410 0.7002 0.7547 0.8055 0.8532

0.4084 0.4938 0.5664 0.6308 0.6892 0.7429 0.7928 0.8396

0.5252 0.5721 0.6196 0.6657 0.7099 0.7523 0.7929 0.8320

0.5533 0.6085 0.6623 0.7135 0.7621 0.8083 0.8524 0.8945

0.5460 0.5990 0.6509 0.7003 0.7472 0.7918 0.8344 0.8752

(B) 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40

0.4197 0.5111 0.5873 0.6534 0.7122 0.7654 0.8142 0.8594

0.4251 0.5193 0.5983 0.6674 0.7296 0.7864 0.8390 0.8881

0.4244 0.5218 0.6034 0.6741 0.7369 0.7938 0.8459 0.8942

0.7179 0.7409 0.7674 0.7961 0.8263 0.8573 0.8885 0.9198

0.7459 0.7861 0.8251 0.8634 0.9009 0.9375 0.9732 0.0081

0.7371 0.7745 0.8116 0.8481 0.8838 0.9188 0.9530 0.9865

(C) 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40

0.9133 0.9742 1.0211 1.0616 1.0984 1.1327 1.1651 1.1960

0.8913 0.9593 1.0147 1.00629 1.1064 1.1466 1.1843 1.2199

0.9324 0.9940 1.0421 1.0840 1.1223 1.1582 1.1923 1.2249

1.0328 1.0639 1.0931 1.1205 1.1470 1.1729 1.1984 1.2236

1.0421 1.0801 1.1158 1.1497 1.1824 1.2140 1.2446 1.2744

1.0360 1.0707 1.1040 1.1362 1.1672 1.1972 1.2263 1.2547

R.S. Gupta, A. Kumar, Approximate analysis for S&fan problems

134

I

METHOD

0

-II

a.2 -

o-1 -

01

0

,

I

0.05

0.1

I

0.15

1

0.2

I

0.25

k

0.3

I

0.35

1

t---c Fig. 2. Interface position versus time graph g(x) = 0.5 - 0.2 cos TX.

well. For comparison we have solved the same problem by the IMM of Crank and Gupta 1121, for which the results are given in Table 2. It may be observed that the results obtained by Method I agree with those of IMM, while the results due to Rasmussen are far out (see Fig. 3). It may be stressed that Rasmussen assumes that the temperature distribution at zero time as well as at subsequent times remains linear, while in the present method we can start from a quadratic temperature profiie in y. 6. The three-dimensional problem In this section we extend Method problems. The problem, communicated

I, presented earlier, for solving three-dimensional to us by Meyer [13], is as follows:

R.S. Gupta, A. Kumar, Approximate analysis for Stefan problems

135

Table 2 Positions of the interface along x = 0 and x = 1 at various times with initial interface profile g(x) = 2 + cos rrx ,X=1

x=0 t

[ 111

Rasmussen

Method I

IMM [12]

[ 1l]

Rasmussen

Method I

IMM [12]

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5

3.0000 2.3461 2.3478 2.4605 2.5881 2.7130 2.8330 2.9483 3.0593 3.1664 3.2699 3.3704

3.0000 3.0140 3.0426 3.1078 3.2059 3.3113 3.4178 3.5224 3.6246 3.7241 3.8211 3.9157

3.0000 3.0770 3.1518 3.2363 3.3281 3.4245 3.5233 3.6231 3.7228 3.8218 3.9197 -

1.0000 1.9864 2.2821 2.4487 2.5859 2.7126 2.8329 2.9483 3.0593 3.1664 3.2699 3.3704

1.0000 2.5071 2.8598 3.0458 3.1840 3.3038 3.4152 3.5215 3.6243 3.7240 3.8210 3.9157

1.0000 1.8298 2.3275 2.6488 2.8909 3.0895 3.2609 3.4138 3.5534 3.6828 3.8044 3.9089

Total time

8.981

5.956

5.416

8.981

5.956

5.600

ALbNG

I

0.5

x : I

I

I.0

1

1.5

-

PRESENT

-----_

I MM

-.-.-

RASMUSSEN

I

4

2.0

2.5

METHOD-I [12] [1 I]

1

3 .o

I

3.5

L.0

t-

Fig. 3. Interface

position versus time graph g(x) = 2 + cos ITX.

I

L.5

5.0

136

R. S. Gupta , A. Kumar, Approximate analysis for Stefan problems

Consider a box B defined by B={x,y,zIO
yG1,

(6.1)

o~zzo.4},

filled with ice at its melting temperature u = 0. All the four vertical sides and the top of the box are insulated. The bottom face z = 0 is subjected to a Dirichlet condition. Due to this, the ice starts melting and the ice/water interface starts moving. Mathematically, we have to solve the following heat conduction equation, expressed in nondimensional form, in the water region 0, Uf = %* +

uyy+ uzz in

L? ,

(6.2)

where 0=(x,

y,z~o~x,y~l,o~z~s(x,

y,t)} *

z = s(x, y, t) denotes the position of the interface at the specified values of x and y at any time t. The corresponding boundary conditions are

u,=o,

x=0

and x=1,

r>O,

(6.3)

z&,=0,

y=O

and y=l,

t>O,

(6.4)

u,=O,

z=o.4,

t>O,

(6.5)

u=U,(X,y,t)={18x(1-x)y(l-y)}25t+0.5f,

u=O

atz=s(x,y,t),

The latent heat condition

z=o,

t>O,

(6.6)

r>O.

(6.7)

at the interface may be written as

St=-uu,+u,s,+us

Y Y'

2 =

(x, y, t) , t > 0.

(6.8)

A quadratic temperature profile in z satisfying the boundary conditions (6.6) and (6.7) for a given x and y may be written as (6.9) As before, cy is determined u,(-u,

Substituting

+ s,u,

+

by deriving the following extra condition

syuy) + u,, + uyy +

u,, = 0

7

at the interface:

z = 4x, Yy 1)*

(6.10)

from (6.9) into (6.10), we obtain the following quadratic equation for ty:

$(I

+ s; + s;> + a[2 + 2U,(l + 3; + s:> - @,, + &) + 4X + sZ)l

+ ui( 1 + s; -I-s;> - u,(s * s,, + s * syy - 2s: - 2s;) = 0 9 the solution of which provides cy.

(6.11)

R.S. Gupta, A. Kumar, Approximate analysis for Stefan problems

137

Table 3 Positions of the interface along the diagonal x = y at various times: upper entry corresponds to Method I and the lower to Meyer [13] x=y

t

0.0

0.1

0.2

0.3

0.4

0.5

0.05

0.0365 -

0.0365 -

0.0365 -

0.0367 -

0.0883 -

0.1663 -

0.10

0.0711 0.06

0.0711 0.06

0.0711 0.06

0.0724 0.0620

0.1760 0.1417

0.2730 0.1867

0.15

0.1057

0.1057 -

0.1057 -

0.1107 -

0.2708 -

0.3691 -

0.20

0.1401 0.1211

0.1401 0.1211

0.1404 0.1214

0.1539 0.1386

0.3661 0.2374

0.4000 0.2758

0.25

0.1742 -

0.1742 -

0.1756 -

0.2017 -

0.4000 -

0.4000

0.30

0.2081 0.1669

0.2083 0.1670

0.2120 0.1684

0.2477 0.2002

0.4000 0.2909

0.4000 0.3241

0.35

0.2419 -

0.2425 -

0.2492 -

0.2887 -

0.4000 -

0.4000 -

0.40

0.2758 0.2028

0.2770 0.2030

0.2863 0.2059

0.3249 0.2461

0.4000 0.3285

0.4000 0.3583

0.45

0.3099

0.3118 -

0.3228 -

0.3573 -

0.4000 -

0.4000 -

0.50

0.3445 -

0.3469 -

0.3582 -

0.3867 -

0.4000 -

0.4000 -

By matching the condition at the interface using (6.9) and (6.11), we arrive at an implicit relationship in s and t which is given by s, =

(u. + a)(1 + S; + S$S

(6.12)

.

The above equation, which is similar to (4.5), time in a step-by-step manner. The numerical results are obtained by choosing The relevant positions of the interface at various Table 3. The corresponding figures from Meyer comparison. It is observed that the deviation proceeds.

provides the position of the interface at any the grid sizes Ax = Ay = times along the diagonal [13], wherever available, between the two results

0.1 and At = 0.001. x = y are shown in are also given for increases as time

138

R.S. Gupta, A. Kumar, Approximate analysis for Stefan problems

References [l] H.G. Landau, Heat conduction in melting solids, Quart. Appl. Math. 8 (1950) 81-94. [2] J. Crank, Two methods for the numerical solution of moving boundary problems in diffusion and heat flow, J. Mech. Appl. Math. 10 (1957) 220-231. [3] W.D. Murray and F. Landis, Numerical and machine solutions of transient heat conduction problems involving melting or freezing, J. Heat Transfer 81 (1959) 106-112. [4] J. Crank and R.S. Gupta, A moving boundary problem arising from the diffusion in absorbing tissue, J. Inst. Math. Appl. 10 (1972) 19-33. [S] R.S. Gupta and D. Kumar, A modified variable time step method for the one-dimensional Stefan problem, Comput. Meths. Appl. Mech. Engrg. 23 (1980) 101-109. [6] T.R. Goodman, Application of integral methods to transient nonlinear heat transfer, in: T.F. Irvine, Jr. and J.P. Hartnett, eds., Advances in Heat Transfer, Vol. 1 (Academic Press, New York, 1964) 51-122. [7] G.E. Bell, A refinement of the heat balance integral method applied to a melting problem, Intemat. J. Heat Mass Transfer 21 (1978) 1357-1362. [8] G.E. Bell, Solidification of a liquid about a cylindrical pipe, Internat. J. Heat Mass Transfer 22 (1979) 1681-1686. [9] G. Poots, An approximate treatment of a heat conduction problem involving a two-dimensional solidification front, Internat. J. Heat Mass Transfer 5 (1962) 339-348. [lo] D.S. Riley and P.W. Duck, Application of the heat-balance integral method to the freezing of a cuboid, Internat. J. Heat Mass Transfer 20 (1977) 294-296. [ll] H. Rasmussen, An approximate method for solving two-dimensional Stefan problems, Lett. Heat Mass Transfer 4 (1977) 273-277. [12] J. Crank and R.S. Gupta, Isotherm migration method in two dimensions, Internat. J. Heat Mass Transfer 18 (1975) 1101-1107. [13] G.H. Meyer, School of Mathematics, Georgia Institute of Technology, Atlanta, GA, personal communication.