VIEWS: 1 PAGES: 18 POSTED ON: 12/17/2009 Public Domain
TR/69 DECEMBER 1976 NUMERICAL SOLUTION OF STEFAN PROBLEMS by A.B. Crowley, Department of Mathematics, Brunel University w9260458 Abstract The enthalpy method for two dimensional Stefan problems is outlined, and applied to the numerical solution of a problem involving the solidification of a square cylinder of fluid, for which experimental data have been published. The same technique is then applied to a similar problem, and the results compared with those obtained by other authors. N o me n c l a t u r e H k L s(t) T t n o n - d i me n s i o n a l e n t h a l p y ; rate of cooling at surface; non-dimensional latent heat; f r e e z i n g f r o n t i n t w o d i me n s i o n a l p r o b l e m; p o s i t i o n o f f r e e z i n g f r o n t i n o n e d i me n s i o n a l p r o b l e m; n o n - d i me n s i o n a l t e mpe r a t u r e ; n o n - d i me n s i o n a l t i me S (x, y, t) = 0 ( = x, y time x thermal diffusivity ) ; (length) 2 Cartesian coordinates; Subscripts L s liquid phase; solid phase. 1. 1. Introduction S t e f a n p r o b l e ms , t h a t i s p r o b l e ms o f h e a t c o n d u c t i o n w i t h c h a n g e o f p h a s e , have attracted much interest in the literature, and recently attention has b e e n c e n t r e d o n t h e s o l u t i o n o f mul t i - d i me n s i o n a l p r o b l e ms . A s n o e x a c t a n a l y t i c a l s o l u t i o n s o f s u c h p r o b l e ms a r e a v a i l a b l e , w e mus t i n g e n e r a l resort to numerical methods, several of which have been proposed ([1], [2], [3] ). In this situation the experimental results of Saitoh [4], for the two dimensional freezing of water, are of particular value, p r o v i d i n g a standard with which numerical results may be compared. I n t h i s p a p e r w e c o n s i d e r t w o p r o b l e ms o f t h e i n w a r d s o l i d i f i c a t i o n o f a square cylinder of liquid, initially at its freezing temperature, with different surface conditions. The first of these, in which the surface t e mpe r a t u r e i s l o w e r e d a t a c o n s t a n t r a t e , c o r r e s p o n d s t o t h e e x p e r i me n t a l configuration employed by Saitoh. The second problem, in which the s ur fa c e t e mpe r a t u r e i s l o w e r e d d i s c o n t i n u o u s l y a t t h e i n i t i a l i n s t a n t , i s o n e f o r w h i c h s e v e r a l s e t s o f n u me r i c a l r e s u l t s h a v e p r e v i o u s l y b e e n p u b l i s h e d . T h e numerical method used here is the straightforward extension to two di me ns i ons of the enthalpy method, used in one dimension by Kamenomo stskaja [5] a n d A t t h e y [ 6 ] . T h i s me t h o d i s f i r s t a p p l i e d t o t h e s o l u t i o n o f S a i t o h ' s p r o b l e m, and good agreement between the numerical results and the experimental data i s o b t a i n e d . T h e me t h o d i s t h e n a p p l i e d t o t h e s e c o n d p r o b l e m, a n d t h e r e s u l t s c o mpa r e d w i t h t h o s e p r e v i o u s l y p u b l i s h e d . T h e p h y s i c a l p r o b l e m a n d t h e e n t h a l p y r e f o r mu l a t i o n We c o n s i d e r a n i n f i n i t e l y l o n g c y l i n d e r o f s q u a r e c r o s s - s e c t i o n , f i l l e d w i t h l i q u i d a t i t s f r e e z i n g t e mpe r a t u r e . T h e l o w e r i n g o f t h e s u r f a c e t e mpe r a t u r e c a u s e s a f r o z e n l a y e r t o g r o w i n w a r d s f r o m t h e s u r f a c e . I n t h i s mo d e l t h e d e n s i t y , s p e c i f i c h e a t a n d t h e r ma l c o n d u c t i v i t y o f t h e s o l i d a r e t a k e n a s c o n s t a n t , a n d i t i s a s s u me d t h a t t h e r e a r e n o d e n s i t y c h a n g e s o n f r e e z i n g . Here there is no convection in the liquid region. I n n o n - d i me n s i o n a l v a r i a b l e s , t h e c y l i n d e r ' s c r o s s - s e c t i o n ma y b e t a k e n a s ( − 1 ≤ x ≤ 1, − 1 ≤ y ≤ 1) a n d t h e t i m e s c a l e c h o s e n i s t h a t o f c o n d u c t i o n , 2 while the freezing temperature is chosen as the origin of the n o n - d i me n s i o n a l t e mpe r a t u r e s c a l e . L e t D b e t h e d o ma i n , i n i t i a l l y zero, occupied by the frozen material, and let the freezing front b e S(x, y, t) = 0 (see fig. l). A suitable non-dimensional model is thus ∂T = ∂ 2T ∂ 2T + ∂t ∂x 2 ∂y 2 in D , (2.1) with T ≡ 0 (2.2) i n t h e l i q u i d r e g i o n , w h e r e T d e n o t e s t h e t e mpe r a t u r e a n d t t h e t i me . A t t h e f r e e z i n g front, S ( x, y, t ) = 0, T = 0, (2.3) and by conservation of energy ∂S L = − (∇T . ∇S) | s , ∂t (2,4) where L denotes the non-dimensional latent heat, and ( ∇ T. ∇ S)| s is the limit of ( ∇ T. ∇ S) as the freezing front is approached in D. In the first problem the surface temperature is lowered at a constant rate, so the boundary condition at the fixed surface is T = − kt on ( x 2 − 1) ( y 2 −1) = 0 (2.5) where k > 0. In the second problem, the temperature is lowered discontinuously. and then maintained constant, so the appropriate condition is T = − 1 on ( x 2 − 1) ( y 2 −1) = 0 , t > 0 . (2.6) 3 (2.1) - (2.5), (2.1) - (2.4) and (2.6) constitute complete descriptions of the two problems to be solved. We next outline the enthalpy reformulation of these problems, and the resulting numerical method We introduce the enthalpy, H, in order to reformulate the problem on the fixed domain (-1 ≤ x ≤ 1, -1 ≤ y ≤ 1). In non-dimensional variables the temperature is defined in terms of the enthalpy by H ≤ 0, ⎧H ⎪ T= ⎨0 <H< L, ⎪ H−L H ≥ L, ⎩ ( 2 . 1 ) , ( 2 . 2 ) t h e n b e c o me ∂ 2T ∂ 2T ∂H = + ∂t ∂x 2 ∂y 2 (2.7) w h e r e T = 0 , H = L c o r r e s p o n d s t o l i q u i d a t i t s f r e e z i n g t e mpe r a t u r e . (2.8) where H is continuous. Across the discontinuity in H, which occurs at the freezing front due to the liberation of latent heat, the conservation form of (2.8) yields ∂S ∂S = L = (∇T.∇S) L − (∇T.∇S) s ∂t ∂t which reduces immediately to (2.4) since in this case T ≡ 0 in the liquid. (H L −Hs ) Thus when the conduction equation is written in the form (2.8), with T = T(H) defined by (2.7), the conditions (2.2), (2.3) are automatically satisfied by any solution of the weak form of (2.8). Oleinik [7] showed that the problem (2.7), (2.8) with appropriate boundary and initial conditions has a unique weak solution, and it may also be shown ( [5], [6] ) that an explicit finite difference converges to this solution. Jerome [8] has shown that an implicit finite difference scheme converges to the weak solution. Thus in order to find a numerical solution, we solve the finite difference forms of (2.7), (2.8) together with the initial condition T = 0, H = L in -1 < x < 1, -1 < y <1, and the appropriate boundary condition (2.5) or (2.6). 4 By symmetry it is sufficient to work on the first quadrant of the square, ∂T ∂T with the boundary conditions = 0 on x = 0 , = 0 on y = 0 . ∂x ∂y 3. Results and discussion All the numerical results for the first problem were calculated using an explicit finite difference scheme to solve (2.7), (2.8) with the appropriate boundary conditions. Fig. 2 shows the position of the freezing front on the axis, and its x coordinate on the diagonal plotted against time. We see that the numerical results agree well with the experimental data, both on the axis and along the diagonal. During an initial time interval in which, away from the corners, the freezing front remains parallel to the fixed surface, heat conduction parallel to the boundaries is negligible. Thus during this interval the solution of the two dimensional problem near the axis is well approximated by the solution of the corresponding one dimensional problem. The exact solution of the one dimensional analogue of Saitoh's problem is not known, but in this case L is large ( ~ 160). Thus we may use a pertubation solution in inverse powers of L, which yields for the position, x = s(t), of the freezing front k t 2 ⎛ k ⎞ 3/2 (3.1) s =1− t + + .... ⎜ ⎟ L 6 ⎝L⎠ We have therefore plotted this solution for comparison with the solution of the two dimensional problem near the axis, and the agreement is excellent until t ~ 6, when the front has moved halfway to the centre. Fig. 3 shows the freezing front at various times. Figs. 4 and 5 show the position of the freezing front on the axis and its x coordinate on the diagonal as functions of time for the second problem, as calculated by the enthalpy method. Also shown are the results obtained by Allen & Severn [1], Lazaridis [2], Crank and Gupta [3]. Here the exact one dimensional solution for the position s(t) of the freezing front on the axis, s( t ) = 1 − α t , [9], where α satisfies α 2 L π αe 2 4 erf α = 1, 2 5 is plotted in Fig. 4 for comparison. The numerical solution for solidification in an infinite corner obtained by Rathjen and Jiji [10] is also shown in Fig. 5 since this is a good approximation to the solution of the present problem while the front remains parallel to the fixed surface away from the corners. The results for the first problem were calculated using an explicit finite difference scheme on a 20 × 20 spatial mesh in 185 seconds of computing time on a CDC 7600. For the second problem an explicit scheme on a 40 × 40 spatial mesh was used to obtain the results shown, in 210 seconds of computing time. However, using a 20 × 20 spatial mesh and 18 seconds computing time, very similar curves were obtained, with a difference of only 3% in the time of final solidification (t ~ 0.625) estimated from these curves. Implicit schemes were also used for the second problem, and the results were in excellent agreement with those on the same spatial mesh, but less economical on computing time. The freezing front as calculated by the enthalpy method is squarer than those from the other numerical schemes. However, as the numerical results for Saitoh's problem are a good fit to the experimental data on both the axis and the diagonal, it would appear that this shape is correct. This squarer shape of the freezing front is also in excellent agreement with the numerical solution for the solidification of an infinite corner [10]. Thus we conclude that the enthalpy reformulation provides a simple and accurate scheme for the numerical solution of multi-dimensional Stefan problems. Acknowledgements The author is most grateful to Professor Saitoh for kindly making available his experimental results. This work was supported by a grant awarded by the Science Research Council. 6 References 1. D. N. de G. Allen and R. T. Severn, The application of relaxation me t h o d s t o t h e s o l u t i o n o f n o n - e l l i p t i c p a r t i a l d i f f e r e n t i a l equations - III, Q. J. Mech. Appl. Math. 15 , 53 - 62 (1962). 2. A. Lazaridis, A numerical solution of the multi-dimensional solidification (or melting) problem. Int. J. Heat Mass Transfer 13, 1459 - 1477 (1970) . 3. J. Crank and R. S. Gupta, Isotherm migration method in two dimensions, Int. J. Heat Mass Transfer 18 , 1101 - 1107 (1975). 4. T. Saitoh, An experimental study of the cylindrical and twodimensional freezing of water with varying wall temperature, Technology Reports, Tohoku Univ. 4l , 61 - 72, (1976). 5. S. L. Kamenomostskaja, On the Stefan problem, Mat. Sb. 53 (95) , 489 - 514, (1961). 6. D. R. Atthey, A finite difference scheme for me lting problems, J. Inst. Maths. Applics. 13 , 353 - 366, (1974). 7. O. A. Oleinik, A method of solution of the general Stefan problem, Sov. Math. Dok1. 1 (2), 1350 - 1354 (1960). 8. J. W. Jerome, Nonlinear equations of evolution and a generalized St e fa n probl e m. 9. (Manuscript). H, S. Carslaw and J. C. Jaeger, Conduction of heat in solids, Chap.XI, Clarendon Press, Oxford (1959). 10. K. A. Rathjen and L. M. Jiji, Heat conduction with melting or freezing in a corner, J. Heat Transfer 93 , 101 - 109 (1971). 7 Fig. 1. Fig. 2. Sketch of the x,y plane. T he x c oordi na t e of th e freezin g fro n t o n th e ax is an d o n t h e d i a g o n a l f o r S a i t o h ' s p r o b l e m. T h e d o t t e d l i n e s h o w s the position of the front in the analogous one d i m e n s i o n a l problem (3.1). Fig. 3. Fig. 4. Position of the freezing front at various times. The x coordinate of the freezing front on the axis for the second problem. The solid line shows the results obtained from the enthalpy method, the broken line the front, s ( t ) = 1 − 1.034 problem. t , for the corresponding one dimensional Fig. 5. The x coordinate of the freezing front on the diagonal. Again the solid line shows the results obtained by the enthalpy method, while the broken line shows the solution for solidification in an infinite corner [10]. 8 9 10 11 12