Method of lines and enthalpy method for solving moving boundary problems

Int. Comm. HeatMass Transfer, Vol.25, No. 4, pp. 531--540, 1998 Copyright © 1998 ElsevierScienceLtd Printed in the USA. All rights reserved 0735-1933/98 $19.00 + .00

Pergamon

PH S0735-1933(98)00040-2

METHOD

OF LINES AND ENTHALPY METHOD FOR SOLVING MOVING BOUNDARY PROBLEMS

S. Srinivas Shastri and R. M. Allen D e p a r t m e n t of Chemical and Process Engineering University of C a n t e r b u r y Christchurch New Zealand

(Communicated by C.L. Tien and P.F. Peterson)

ABSTRACT The enthalpy method is a simple and efficient method of solving moving boundary problems associated with melting and solidification. The numerical techniques used to solve the equations are typically explicit and implicit methods. The combination of a more effective mathematical technique, the method of lines, with the enthalpy method h~s yielded a more stable and easily proga'ammed solution which has been successfully tested against previous results. © 1998 ElsevierScienceLtd Introduction A number of methods have evolved since Stefan first formulated the melting problem (Crank [1]). The equations that describe the standard melting/solidification problem can be written ms 0Ts

dS(t) dt

~Ts - o.~ a:_~

(1)

Ot - at ,9='-

('2)

h'~ ffl~ Ps L /9-=

531

kt OTj Pt L a z

(3)

532

S. Srinivas Shastri and R.M. Allen

Vol. 25, No. 4

where T i s the temperature, k thermal conductivity, (~ thermal diffusivity, L latent heat of fusion, S position of the interface, p density, t time and z spatial coordinate. Subscripts s and I refer to solid and liquid phases respectively. Equation 3 is the result of the boundary condition at the interface ~'~ = T~ = 7;.,~,

.t

: =

s(t)

(4)

The boundary conditions are T(t, o) -

aT(t, 0) o~-

+ ~t~(t)

0T(t, l) T(t, 0 =

O---2--- + 72(t)

(5)

0
0 < t <

(6)

oo

and the initial condition

~(0,z) = T0(z)

(7)

0<:
where "y is some function of time and l the length under consideration.

The Problem

The solution of the problem is straightforward except for the determination of the interface. Once the position of the boundary is fixed then the solution of the differential equations is relatively simple.

Melting or solidification may occur over a t.emperature range. Problems of this type are more realistic.

The Methods

Most numerical techniques employed to solve the differential equations associated with moving boundary problems involve approximations by finite differences (Crank [1]).

M e t h o d o f I n v a r i a n t Imbedding This method was developed by Meyer [2]. The one phase Stefan problem is approximated by a sequence of free boundary problems for ordinary differential equations. These equations axe solved by converting them to initial value problel~s by applying the Riccati transformation.

The one phase Stefan problem is

02T

LT =- ~

OT

+ a(t, z) ~ - b(t, z ) T - c(t, z)

~T

= f(t, -)

(s)

Vol. 25, No. 4

METHODS FOR MOVING BOUNDARY PROBLEMS

533

with boundary condition a.s in Equation 5 and initial condition T(0, z ) = T 0 ( z )

0
S(0) > 0

(9)

Neglecting the mathematical manipulations involved, the above equations can be converted to a series of ordinary differential equations, and an algebraic equation. The algebraic equation is solved iteratively to determine the position of the interface.

Variable T i m e Step G u p t a and Kumar [3] have reviewed some variable time step methods, and have developed the modified variable time step method. The modification is in the accuracy check which is made at the interface rather than at the fixed boundaries.

In variable time step methods, the melting (solidifying) interface is allowed to move in steps of only one grid point, and the time is determined iteratively. This method presupposes that the interface is at the fusion temperature. While this is correct for pure substances, it is not the case of mixtures or alloys. Thus, in spite of the ease of implementation, the applicability of this method is limited.

Enthalpy Method The enthalpy method also uses a finite difference approximation to solve the partial differential equation. If enthalpy becomes the dependent variable, then Equations 1 to 4 can be replaced by

01-1 Ot

k 02T p c9:~ + p

(10)

where H is the enthalpy and ¢ is the body heating term. This equation is valid over the entire solution domain, in both the solid and liquid phases a.s well ~s at the interface. The above equation can be solved by explicit means, or any other suitable method. The position of the interface is estimated by linear interpolation between two times by interpreting the enthalpy (Voller and Cross [4J).

Voller and Cross [4] discuss an implicit method where the interface is allowed to move only one grid point at a time, and the time required for this is determined iteratively. In this implicit method, the temperature T is replaced by a Heaviside function of the form

T = F(H)

(1l)

F ( H ) = H/C H < L

(12)

where

534

S. Srinivas Shastri and R.M. Allen

Vol. 25, No. 4

1:(11) = Tm~u O <_ H < L F( H) -

(H-L) C

H>

L

(13) (14)

and C is the specific heat. The implicit equation itself would have the form H j+l = ; 2 : 2 F ( H J + I )

(15)

F i ( H ) = [F(Hi-1) - 2F(H~) + F(Hi+I)]

(16)

where

The non linearity makes this equation difficult to solve.

The above set of equations can be rewritten in standard non linear parabolic form

OH

ka2V(H) = p

Oz ~

e~ + P

(,r)

and the method of lines technique applied.

M e t h o d o f Lines a n d E n t h a l p y M e t h o d The combination of the numerical techniques for solving parabolic partial differential equations by the method of lines and the enthalpy method h~s not been previously exploited. The method of lines technique converts the partial differential equation into a set of simultaneous ordinary differential equations (Hyman [5]). These equations can be quite stiff, possibly one of the reasons in the method of lines not being a popular technique (Polls and Goodson [61).

The mathematical solver CVODE (Hindmarsh and Cohen [7]) has been used to implement the solution. CVODE is a solver for stiff and non stiff initial value problems for systems of ordinary differentia.1 equations. CVODE integrates the ordinary differential equation over a.n interval in time. Depending on user option, CVODE (in N O R M A L mode) integrates from its current internal time value to a point, at or beyond t o u t , the desired time, and interpolates to return the value of the dependent variable at tout. For stiff problems, C V O D E uses the backward differentiation formula: and the Newlon iteration method to solve the linear system of equations.

A welding problem (Voller and Cross [4]) has been solved by the above method, and the results are compared with the original work. This example illustrates the application of the enthalpy method with the method of lines to a typical situation.

Vol. 25, No. 4

METHODS FOR MOVING BOUNDARY PROBLEMS

535

o.~G Yk\\ 0.54

""

"

. ~

o..s I I

\

,. so,idos ~-,~".,

\':L

o..21-co, o..o|r;

10

~ p

12

2 CenTre

!

14

~_____. ~" ~ 16

18

,

20

PosiTion o'f'-solidus,cenl"re,liquidus boundaries,

~I

22

I() 4 M

FIG. 1 Movement of Boundary - Voller and Cross [4]

The Welding Problem

The spot welding of two large steel plates is considered. The weld length is 2.2ram and the initial condition is T(z, 0) = 0 and boundary conditions are 0T --=0

Oz OT

z=l

o_=TM : = o

The body heating term is ¢ ----[197"/48 + 50} × 10s T < = 800°C ¢ = [ ( T - 800)/16 + 1100/3] x 10s T > 800°C

TABLE 1 M e t a l P r o p e r t i e s &: I n t e g r a t i o n S t e p s i z e s specific heat thermal conductivity density latent heat melting point melt range, 2~ dz dt

6203/kg°C 24W/m°C 7800kg/m 3 271 x 10:~J/kg 1530°C 20°C 10-4m [41 2.2 x 10-4s [4]

536

S, Srinivas Shastri and R.M. Allen

Vol. 25, No. 4

0.70 "[he Weldlng Problem

1: Problem 4.4

0.60 ? E

0.50 solidus ~ 0,4-0

I I I I l l l l

Illl

fill

Ilql

e Ill]

II1111111

IIII

II

5.00E-004 1,00E-003 1.50E-003 2.00E-003 2.50E-005 position of the interface (m) FIG. 2 Movemeat of Boundary - thLs work

0.70 \

\\ ~

"[he WeldingProblem 1981:Problem4.4

Q.60 6q

0.50 E

\ 0.4-0

~

"----___~nt~e ldus

0 30 5.0E-004 1.0E-003 1.5E-003 2.0E-003 2.5E--003 position of the interface (m) FIC. 3 Movement of Boundary - time step is lO-3s

Vol. 25, No. 4

METHODS FOR MOVING BOUNDARY PROBLEMS

2400

/ 2 . 1 0 -z /,I.5.z6

(b)

z

2000

t800

' 10-3

S -

71./

• O. 5 10-3

1400

E ~. io00 ,2oo ~_ 800

600 400 200 0.20

0.40 Time,

0.60 $

J

0.80

FIG. 4 T e m p e r a t u r e History - Voller and Cross (1981)

5000 Temperature Profiles 2.0 * 1 0 4 ~

~

2500

or02000

10-3

21500 2 ©

O -~ m

E© 1ooo

500

0 0.0

I

0,20

I

0.40 t i m e (s)

I

0.60

FIG. 5 T e m p e r a t u r e History - this work

0.80

537

538

S. Srinivas Shastri and R.M. Allen

Vol. 25, No. 4

Results and Discussion

T h e invariant imbedding and variable time step methods have been briefly described to highlight the limitations of other approaches compared to the enthalpy method.

These approaches are

restricted to those substances t h a t have a fixed fusion temperature. However, Meyer [8] has applied the m e t h o d of invariant imbedding to a solidificatioa problem.

This m e t h o d is m a t h e m a t i c a l l y

complicated as well as numerically intense.

T h e enthalpy m e t h o d has done away with the need to determine the exact location of the interface for the solution to evolve. Instead the interface position is determined from the enthalpy by linear interpolation. It is this feature t h a t is of importance to engineering problems. T h e usefulness of the enthalpy m e t h o d is the ease with which it deals with 'mushy problems'. In the ca.se of this class of problems, Equations 12, 13 and 14 become continuous, with the latent heat becoming a function of t e m p e r a t u r e over the fusion range. T h e determination of the solidus and liquidus positions is also by interpolation.

T h e results from our solution compare favourably with those of Voller and Cross.

Figure 2

shows the results of this work using the same time and spatial steps as Voller and Cross (1981). T h e irregularities in the liquidus curve are due to the stringent tolerance used (10 .6 J / k g ) in the present work. These disappear when the time step is made smaller (10-3s) from 2 x 10 2s (Figure 3). This suggests t h a t the equations are stiff. Furthermore the overall accuracy is improved by choosing snch a stringent tolerance criterion and a very small time step; as the time step tends to 0, the actual differential equation is approached.

T h e t e m p e r a t u r e histories also agree (Figures 4 and 5). It is observed t h a t the plateau in Figure 5 extends over a longer time interval t h a n in Figure 4 because our simulation predicts a thicker mushy zone, and this acts aus a larger heat reservoir. The response to heat i n p u t is therefore more sluggish.

Conclusions

T h e m e t h o d of lines in combination with the enthalpy m e t h o d is a powerful tool to solve moving b o u n d a r y problems t h a t occur in solidifcation and melting. T h e availability of stiff solvers like C V O D E has made it e~.siel to tackle difficult problems with fewer c o m p u t a t i o n a l complications. T h e ability to use very small time steps coupled wit h stringent aecm'acy checks make this method

Vol. 25, No. 4

METHODS FOR MOVING BOUNDARY PROBLEMS

539

favourable. Tile nature of the partial differential equation presents the interesting possibility of applying optimal control techniques to such problems. Acknowledgements

S. S. Sha.stri wishes to acknowledge the NZODA Postgraduate Scholarship and the New Zealand Aluminium Smelters Postgraduate Scholarship in Process Control for this work. Nomenclature C H k kl ks L l T

specific heat enthalpy thermal conductivity liquid thermal conductivity solid thermal conductivity latent heat length temperature liquid temperature Tmett fusion temperature Ts solid temperature t time z spatial coordinate S(t) position of interface al liquid thermal diffusivity as solid thermal diffusivity half melting range ¢ power dermity p density of region Pl liquid deI~sity Ps solid density References

1. Crank. J., Free and Moving Boundary Problems, Clarendon Press, (1984) 2. Meyer, G.H. Numer. Math, l_fi, 248-267 (1970) 3. Gupta, R.S. and Kumar, D.. Int. J. Heat Mass Transfer, 24, 251-259 (1981) 4. Volter, V. and Cross, M., Int. J. Heat Mass Transfer, 24, 545-556 (1981)

540

S. Srinivas Shastri and R.M. Allen

Vol. 25, No. 4

5. llyman, J. M., Method of lines solution of partial differential equations, National Technical Information Service, Springfield, Va., (1977) 6. Dolis, M. P. and Coodson, R. E., Proc. IEEE, 6._3_4(I), 45-61 (1976) 7. Cohen, S. D. and Ilindmarsh, A. C; CVODE user gaJide, (1994) 8. Meyer, G.H., Int. J. Heat Ma.ss Transfer, 2__4,778-781 (1981)

Received August 8, 1997