- Email: [email protected]
A numerical analysis of freezing and melting with convection
COMPUTER METHODS IN APPLIED MECHANICS NORTH-HOLLAND PUBLISHING COMPANY
A NUMERICAL
AND ENGINEERING
ANALYSIS OF FREEZING CONVECTION
28 (1981) 275-284
AND MELTING
WITH
K. MORGAN Department
of Civil Engineering,
University of Wales, Swansea
SA2 8PP, United Kingdom
Received 17 November 1980 Revised manuscript received 26 March 1981
An explicit finite element algorithm is presented for the solution of the basic equations describing combined conductive and convective transfer of heat in materials which may undergo a liquid/solid change of phase. The method is used to analyse the influence of natural convection on the process of freezing and melting in a cylindrical thermal cavity. The results produced in the case of freezing are shown to agree qualitatively with experimental observations.
1. Introduction The numerical modelling of heat transfer problems involving a solid/liquid change of phase has been an area which has attracted much attention [l]. Problems of this type are currently being studied to investigate the possibilities of thermal energy storage via the phase change process and there are also many practical phase change problems which are of importance in the more traditional steel and nuclear industries. In the great majority of cases, the numerical models which have been developed are based on heat conduction theory and the effects of convection are neglected, whereas it is known that convection plays a dominant role in many situations of practical interest. This has been confirmed by a recent series of experiments [2,3] which shows the importance of natural convection in the processes of melting and freezing. Studies which have attempted to take account of convection have generally either included simplifying assumptions about the geometry of the region under consideration [4,5] or used empirical heat transfer coefficients [6]. A significant advance was provided by the work of Gartling [7] who removed the need for restrictions of this type. By using the finite element method he produced a technique which could be applied to two-dimensional regions of arbitrary shape, and the full Boussinesq equations for coupled laminar flow with heat transfer were solved in the fluid. A quadratic finite element representation for temperature was found to be necessary for accurate modelling of the phase change process by the enthalpy method [8] with the implicit method of time-stepping which was adopted. The method presented here retains the full generality of the work of Gartling [7] by solving the same equations via an explicit finite element technique. Accurate modelling of phase change problems can then be obtained by using linear elements to represent the temperature field [9]. An important additional feature of the method is that it could be used, in an extended form, to analyse similar phase change problems in three-dimensional regions of arbitrary shape where the use of implicit methods to solve the fully coupled system of equations could prove to be impractical [lo]. 00457825/81/0000-0000/$02.75
@ 1981 North-Holland
K. Morgan, Numerical analysis of freezing and melting with convection
276
2. The governing equations An algorithm is to be developed which is capable of analysing problems in which there exists initially a liquid (solid) region which may change phase to solid (liquid), as time proceeds, due to heat transfer processes. Heat transfer in a fluid is governed by the basic equations expressing conservation of mass, momentum and energy. If the Boussinesq approximation [ll] is made and the summation convention adopted, then in two dimensions the mass conservation equation can be written in the form au,/ax, = 0
(i = 1,2)
(1)
where u1 and u2 are the velocity components in the direction of the Cartesian coordinate axes x1 and x2 respectively and where the x2 axis is assumed to be vertically upwards. If p is the fluid density, g is the acceleration due to gravity and p is the coefficient of conservation equations may be expressed volume expansion of the fluid, then the momentum as
= *-axj
pgP(T - T,,,)2
(2) I
where CTij
=
-p&j
+
/L
(Z+$!)
(3)
and p=p’+pgxz
(4)
and Tret is a reference Here p’ is the pressure, p is the viscosity, T is the temperature temperature for which buoyancy forces are zero. Conservation of energy in the fluid then requires satisfaction of the equation
(5) where c is the specific heat and k is the thermal conductivity. When a phase change from liquid to solid occurs, any accompanying change of density will be neglected and the heat transfer in the solid is then governed by (5) with the velocity components Uj set to zero.
3. Modelling the phase change process In constructing a solution to a problem involving a liquid/solid change of phase, a possible approach would be to accurately track the position of the phase boundary [l, 12,131 and then to solve equations (l-5) in the fluid region and the appropriate form of (5) in the solid region.
K. Morgan, Numerical analysis of freezing and melting with convection
277
These solutions would have to satisfy the boundary conditions on temperature, velocity and heat removal at the phase change surface [14]. However, the complexity of front tracking methods in two dimensions led to the adoption in the present work of the alternative approach in which the phase change process is modelled by a variant of the enthalpy method [S]. In this method the phase change is assumed to occur over a temperature range and the associated latent heat effect is handled by increasing suitably the specific heat in this range. Thus if the phase change is assumed to occur over the temperature interval [T,, T,], where T, is the solidus temperature and T, the liquidus temperature, then the specific heat cA used in the calculation is defined by CA = c + [H(T
- Ts)-
H(T
- T,)]L/AT
(6)
where H denotes the Heaviside function, L is the latent heat and AT = T, - T, is the phase change interval. When this method is applied to the analysis of pure materials, in which the phase change occurs at a specified temperature (i.e. T, = T,), a phase change interval AT (# 0) must be assumed. Then it has been demonstrated [7,15,16] that good results can be obtained for problems involving conduction provided that the right choice is made for the value of AT. However, for materials in which the phase change does indeed occur over a reasonable temperature range this problem does not arise and the actual physical values of T, and T, can then be used successfully [17]. The enthalpy method is used here in conjunction with a fixed finite element mesh. Gartling [7] noted that such an approach was well suited to the computation of the fluid motion except for the problem of the application of the no-slip bounda~ condition at the phase boundary. Exact application of this boundary condition of zero-fluid velocity at the phase change surface will not be possible as the fluid/solid interface will, in general, not coincide with the nodes of the finite element mesh. Gartling [7] overcame this problem by using a smearing approach in which the viscosity was greatly increased in the phase change region, but this method was found to be unsuccessful within the context of the explicit solution algorithm to be adopted here. Instead, the problem was overcome by adopting the simpler, but perhaps less accurate, approach of reducing the nodal velocity to zero whenever the nodal temperature lies below the liquidus. 4. Finite element solution procedure The approach outlined in the previous section for dealing with a liquid/solid change of phase was implemented in the computer code CONDIF, developed at J.R.C. Ispra by J. Donea et al. [18] for the analysis of transient convective-conductive heat transfer problems. Space discretisation is performed via the finite element method while time-stepping is achieved by the use of the fractional step method 1181. This method was originally used by Chorin [19] in the finite difference solution of incompressible flow problems. The solution procedure adopted in the CONDIF code has been widely reported in the literature [18,20,21], and a brief description of the salient points will therefore suffice here. The spatial discretisation is accomplished via the use of 4-noded isoparametric e’tements
278
K. Morgan, Numerical analysis of freezing and melting with convection
with associated bilinear velocity and temperature fields and elementwise constant pressure. Application of the standard method of weighted residuals [22] to equations (l)-(5) then yields a matrix equation system of the form
cu=o,
(7)
du
Mz=fa+fd+fb+
cp,
dT Madt=ga+g,+g,.
(9)
Here C is the pressure matrix, M is the mass matrix, MA the capacitance subscripts a, d and b denote the vectors arising from the advection terms, the and the boundary conditions and body force terms respectively. With the solution known everywhere at time tn, the fractional step method determine the solution of (7) and (8) at the next time level t,+* (= t, + At). In intermediate velocity field u* is calculated so as to satisfy only a discretised without the pressure term, i.e., U
*= u”+AtM-‘(f:+f:+f:).
matrix and the diffusion terms is then used to this method an version of (8)
(10)
The required velocity and pressure fields, untl and p”+l, are then obtained by adding to u* the dynamical effect of pressure, determined so as to ensure that the incompressibility condition remains satisfied, i.e., AtC’M-‘Cp”+’ U “+l =
= - C’U* ,
u* + AtM-‘Cp””
.
(11)
(12)
Problems associated with the application of velocity boundary conditions ensure that the matrix C’M-‘C in (11) splits into two independent submatrices, each of which is singular. It has been argued [23] that this may be due to the fact that incorrect pressure equations are produced in the boundary elements if the conditions on prescribed tangential velocity are strictly applied. As a consequence, the boundary conditions on the tangential component of velocity are satisfied by u”+’ in a weak sense only, thus ensuring that the pressure equations have the correct form in the boundary elements and enabling (11) to be solved, upon specification of a single reference value for pressure. With the velocity and pressure at time so determined, the solution is completed by evaluating the temperature at this time from g)‘in the form T “+I= T” + AtMi’(g:
+ g: + g:) .
(13)
In explicit schemes of this type, stability considerations restrict the size of the time step At which may be used, and it is therefore computationally advantageous to replace the consistent
K. Morgan, Numerical analysis of freezing and melting with convection
279
mass and capacitance matrices in equations (lOb(13) by the lumped ~diagona1) matrices obtained by adding all terms of each row of the consistent matrices and placing the result on the diagonal. A detailed analysis of the effects of this approximation and the determination of the resulting stability requirements has been made by Donea et al. [18]. Their results are used here and then, at any stage in the process. The time step At must satisfy the relation
PcIfr12 PIhI IhI ’ 4~
4k
’ uihi
(14)
where L = (hl, hz) is the element side vector. If a temperature dependent specific heat c is used in (6), the lumped capacitance matrix has to be recalculated at each time step. This can be avoided if the specific heat in the solid and liquid regions can be assumed separately constant (though not necessarily equal). Then recomputation is required only for the elements currently undergoing a change of phase, and the computational efficiency of the scheme is improved. The efficiency of the scheme is not affected by the inclusion of a temperature dependent thermal conductivity, k, in (5) since the diffusion matrix gd in equation (9) is already recalculated at each time step.
5. Freezing in a thermal cavity
This finite element scheme has been shown [9] to be capable of providing accurate solutions to phase change problems involving only conduction, but it is difficult to assess the accuracy of its predictions for problems involving convection due to the scarcity of suitable experimental data. However, in a recent paper Sparrow et al. [3] describe results obtained in a series of experiments designed to investigate the effect of natural convection on the process of freezing. In these experiments a liquid is placed at temperature To in a cylindrical container whose ends are thermally insulated and whose outer surface is held fixed at the temperature To. Freezing commences when a cooled inner cylinder at temperature Ti is inserted into the liquid. A detailed analysis of these experiments was not possible as the author did not have access to the full material property data. It was therefore decided to perform calculations in a similar geometry and subject to similar boundary conditions with a view to comparing the general experimental and numerical conclusions. The dimensions of the region analysed and the type of boundary conditions applied are shown diagrammatically in Fig. 1. The region was covered by a mesh of 100 square elements and boundary conditions of zero velocity were applied on the walls of the region at all times. Although somewhat unrealistic physically, the material properties were taken to be constant and equal in both phases and the values adopted for the problem parameters were defined by p=
1.0,
k =
c=
1.0,
p = 0.003 )
Tre3,= cl ,
0.001 , (13
L = 5.0 .
The phase change in the material was assumed to occur over a temperature
range defined by
280
K. Morgan,
Numerical
analysis of freezing
\’
01
Fig. 1. Region
analysed
T, = -0.6,
and temperature
and melting with convection
’
‘\ar
’
’
’
azzo
boundary
T, = -0.65
’
’
’
’
11
conditions applied. Initially T = TO everywhere.
(15b)
and, in the problems to be considered here, such a range should certainly ensure that the enthalpy method takes correct account of the latent heat effect [17]. In the initial calculations, the value of g was adjusted so as to give a Rayleigh number of 3000 when To- Ti = 1, as previous experience with the model [21] had shown that the coarse mesh adopted would then give a good representation of the flow. The computer runs performed are detailed in Table 1. In Run A the liquid is initially at a temperature only just above the liquidus value, and the effects of convection in this case are negligible. The frozen region is then in the form of a perfect cylinder at each instant and this is illustrated in Fig. 2 which shows the progress of the liquidus to time 3000 when the computation ceased. In the remaining Runs B-D the liquid is superheated and the table shows how natural convection slows and ultimately stops the freezing process. Plots of the isotherms and streamlines at steady state conditions are shown in Fig. 3. Close examination of the numerical results produces the same conclusions as those reached by analysing experimental observations [3] viz. increasing To at a fixed value of Ti:,,or decreasing -Ti at a fixed value of To, decreases the
Table 1 Details of freezing calculations Run
TO
Ti
Time for steady state
A B
-0.59
c
0 2
-1 -1 -2 -2
3000’ 648 1578 299
D
0
‘Denotes end of computation conditions not attained.
time; steady state
281
K. Morgan, Numerical analysis of freezing and melting with convection
Run B
Run C
(al Fig. 2. Position of the liquidus at various times in Run
Run II
lb)
Fig. 3. Temperature (a) and stream function,(b) at steady state conditions for Runs B, C and D. The solidified region is shaded.
A.
time at which steady state is attained and decreases the volume of liquid which is eventually frozen. The effect of the magnitude of the latent heat was investigated by repeating Run B firstly with L = 0 and then with L = 10. It was found in all three cases that the distribution of temperature and stream function at steady state was virtually identical but that the time at
RunE
(a)
Fig. 4. Temperature shaded.
Run F
(a) and stream function (b) at steady state conditions
for Runs E and F. The solidified region is
K. Morgan,
282
Numerical
analysis
of freezing
and melting with convection
which steady state conditions were reached increased with latent heat from 414 with L = 0 to 918 with L = 10. In Runs E and F, Run B was repeated with the value of g increased by a factor of ten (Rayleigh number = 30 000) and by a factor of fifty (Rayleigh number = 150 000) respectively. The effect of convection is then significantly increased and a finer mesh was necessary in Run F at the highest Rayleigh number. Again the effect of the increased convection is to reduce the final amount of frozen material and also to reduce the time at which steady state conditions result. The steady state was attained at time 249 in Run E and at time 184 in Run F. The distribution of temperature and stream function at steady state for these two cases is shown in Fig. 4.
6. Melting in a thermal cavity To further demonstrate the capabilities of the model, the inverse problem of melting in a cylindrical thermal cavity was also considered. The geometry of the previous section is used but now the initial temperature of the region is taken to be below the solidus, and melting is initiated by increasing the temperature of the inner wall. The values given in (15) were again Table 2 Details of melting calculations Run
TO
Ti
Time for steady state
G H
-0.66 -1
2 2
1900 2609
RunG
Run H
Run
(al
I
RunJ
Fig. 5. Temperature (a) and stream function (b) at steady state conditions for Runs G, H, I and J. The solidified region is shaded.
K. Morgan, Numerical analysis of freezing and melting with convection
283
used for the thermal properties but T,, was put equal to the liquidus T,. Two basic Runs (G and H) were made as detailed in Table 2 and the effect of increasing the amount of convection was demonstrated by repeating Run H with the value of g increased by a factor of ten (Run I) and by a factor of fifty (Run J). As the effect of convection is increased in this manner, the steady state is reached in a shorter time (in Run I at time 945 and at time 436 in Run J) and the amount of material remaining frozen increases. For this complete set of calculations, the temperature and stream function distributions at steady-state are shown in Fig. 5.
7. Conchdons An explicit finite element method for analysing problems involving melting and freezing in the presence of convection has been described. The method has been shown to be capable of producing results which agree qualitatively with experimental observations, but a detailed analysis of the accuracy of the method is still not possible because of the scarcity of suitable experimental results. The explicit nature of the method means that it should be possible in the future to extend the method to the analysis of three-dimensional phase change problems involving convection.
The work reported in this paper was performed while the author was a Visiting Scientist in the Applied Mechanics Division of the Joint Research Centre of the European Communities, Ispra, Italy. The author wishes to thank the director of the J.R.C. for providing the financial support which made the visit possible and the staff of the Applied Mechanics Division, particularly Dr. J. Donea, for their help and encouragement.
References [l] J. Crank, How to deal with moving boundaries in thermal problems, in: R.W. Lewis, K. Morgan and O.C. Zienkiewicz, eds., Numerical Methods in Heat Transfer (Wiley, Chichester, 1981). [2] E.M. Sparrow, R.R. Schmidt and J.W. Ramsey, Experiments on the role of natural convection in the melting of solids, J. Heat Transfer 100 (1978) 11-16. [3] E.M. Sparrow, J.W. Ramsey and R.G. Kemink, Freezing controlled by natural convection, J. Heat Transfer 101 (1979) 578-584. [4] C.A. Lapuduta and W.K. Mueller, The effect of buoyancy on the formation of a solid deposit freezing on a vertical surface, Internat. J. Heat Mass Transfer 13 (1970) 13-25. [5] E.M. Sparrow, S.V. Patankar and S. Ramadhyani, Analysis of melting in the presence of natural convection in the melt region, J. Heat Transfer 99 (1977) 520-526. [6] J. Szekely and P.S. Chhabra, The effect of natural convection on the shape and movement of the melt-solid interface in the controlled solidi~cation of lead, Metallurg. Trans. (1970) 1195-1203. [7] D.K. Gartling, Finite element analysis of convective heat transfer problems with change of phase, in: K. Morgan, C. Taylor and C.A. Brebbia, eds., Computer Methods in Fluids (Pentech, London, 1980). [8] D.R. Atthey, A finite difference scheme for melting problems, J. Inst. Math. Appl. 13 (1974) 353-366.
284
K. Morgan, Numerical
analysis of freezing and melting with convection
[9] K. Morgan, An explicit method for convection problems involving a change of phase, SMiRT-6, Paris, 1981, Paper B278. [lo] P.M. Gresho, R.L. Lee and R.L. Sani, On the time-dependent solution of the incompressible Navier-Stokes equations in two and three dimensions, in: C. Taylor and K. Morgan, eds., Recent Advances in Numerical Methods in Fluids (Pineridge, Swansea, 1980). [ll] J.C. Heinrich, P.S. Marshall and O.C. Zienkiewicz, Penalty function solution of coupled convective and conductive heat transfer, in: C. Taylor, K. Morgan and CA. Brebbia, eds., Proc. Conf. Numer. Meths. in Laminar and Turbulent Flow (Pentech, London, 1978). [12] G.H. Meyer, The numerical solution of multidimensional Stefan problems-a survey, in: D.G. Wilson, A.D. Solomon and P.T. Boggs, eds., Moving Boundary Problems (Academic Press, New York, 1978). [13] K. O’Neill and D.R. Lynch, A finite element solution of freezing problems using a continuously deforming coordinate system, in: R.W. Lewis, K. Morgan and O.C. Zienkiewicz, eds., Numerical Methods in Heat Transfer (Wiley, Chichester, 1981). [14] H.S. Carslaw and J.C. Jaeger, Conduction of heat in solids (Oxford University Press, London, 2nd ed., 1959). [15] G. Comini, S. Del Guidice, R.W. Lewis and O.C. Zienkiewicz, Finite element solution of non-linear heat conduction equations with special reference to phase change, Internat. J. Numer. Meths. Engrg. 8 (1974) 613-624. [16] K. Morgan, R.W. Lewis and O.C. Zienkiewicz, An improved algorithm for heat conduction problems with phase change, Internat. J. Numer. Meths. Engrg. 13 (1978) 1191-1195. [17] K. Morgan, R.W. Lewis and H.R. Thomas, Finite element modelling of drying stresses in timber and cooling stresses in cast metal, in: M.J. O’Carroll et al., eds., Modelling and Simulation in Practice, Vol. 2 (Emjoc, Northallerton, 1980). [18] J. Donea, S. Giuliani, H. Lava1 and L. Quartapelle, Solution of the unsteady Navier-Stokes equations by a finite element projection method, in: C. Taylor and K. Morgan, eds., Recent Advances in Numerical Methods in Fluids, Vol. 2: Computational Techniques in Transient and Turbulent Flow (Pineridge, Swansea, 1981). [19] A.J. Chorin, Numerical solution of the Navier-Stokes equations, Math. Comput. 22 (1968) 745-762. [20] J. Donea and S. Giuliani, The computer code CONDIF for transient convective-conductive heat transfer, Ispra Report EUR 6822/I EN, 1980. [21] J. Donea, S. Giuliani and H. Lava], Accurate explicit finite element schemes for convective
A NUMERICAL
AND ENGINEERING
ANALYSIS OF FREEZING CONVECTION
28 (1981) 275-284
AND MELTING
WITH
K. MORGAN Department
of Civil Engineering,
University of Wales, Swansea
SA2 8PP, United Kingdom
Received 17 November 1980 Revised manuscript received 26 March 1981
An explicit finite element algorithm is presented for the solution of the basic equations describing combined conductive and convective transfer of heat in materials which may undergo a liquid/solid change of phase. The method is used to analyse the influence of natural convection on the process of freezing and melting in a cylindrical thermal cavity. The results produced in the case of freezing are shown to agree qualitatively with experimental observations.
1. Introduction The numerical modelling of heat transfer problems involving a solid/liquid change of phase has been an area which has attracted much attention [l]. Problems of this type are currently being studied to investigate the possibilities of thermal energy storage via the phase change process and there are also many practical phase change problems which are of importance in the more traditional steel and nuclear industries. In the great majority of cases, the numerical models which have been developed are based on heat conduction theory and the effects of convection are neglected, whereas it is known that convection plays a dominant role in many situations of practical interest. This has been confirmed by a recent series of experiments [2,3] which shows the importance of natural convection in the processes of melting and freezing. Studies which have attempted to take account of convection have generally either included simplifying assumptions about the geometry of the region under consideration [4,5] or used empirical heat transfer coefficients [6]. A significant advance was provided by the work of Gartling [7] who removed the need for restrictions of this type. By using the finite element method he produced a technique which could be applied to two-dimensional regions of arbitrary shape, and the full Boussinesq equations for coupled laminar flow with heat transfer were solved in the fluid. A quadratic finite element representation for temperature was found to be necessary for accurate modelling of the phase change process by the enthalpy method [8] with the implicit method of time-stepping which was adopted. The method presented here retains the full generality of the work of Gartling [7] by solving the same equations via an explicit finite element technique. Accurate modelling of phase change problems can then be obtained by using linear elements to represent the temperature field [9]. An important additional feature of the method is that it could be used, in an extended form, to analyse similar phase change problems in three-dimensional regions of arbitrary shape where the use of implicit methods to solve the fully coupled system of equations could prove to be impractical [lo]. 00457825/81/0000-0000/$02.75
@ 1981 North-Holland
K. Morgan, Numerical analysis of freezing and melting with convection
276
2. The governing equations An algorithm is to be developed which is capable of analysing problems in which there exists initially a liquid (solid) region which may change phase to solid (liquid), as time proceeds, due to heat transfer processes. Heat transfer in a fluid is governed by the basic equations expressing conservation of mass, momentum and energy. If the Boussinesq approximation [ll] is made and the summation convention adopted, then in two dimensions the mass conservation equation can be written in the form au,/ax, = 0
(i = 1,2)
(1)
where u1 and u2 are the velocity components in the direction of the Cartesian coordinate axes x1 and x2 respectively and where the x2 axis is assumed to be vertically upwards. If p is the fluid density, g is the acceleration due to gravity and p is the coefficient of conservation equations may be expressed volume expansion of the fluid, then the momentum as
= *-axj
pgP(T - T,,,)2
(2) I
where CTij
=
-p&j
+
/L
(Z+$!)
(3)
and p=p’+pgxz
(4)
and Tret is a reference Here p’ is the pressure, p is the viscosity, T is the temperature temperature for which buoyancy forces are zero. Conservation of energy in the fluid then requires satisfaction of the equation
(5) where c is the specific heat and k is the thermal conductivity. When a phase change from liquid to solid occurs, any accompanying change of density will be neglected and the heat transfer in the solid is then governed by (5) with the velocity components Uj set to zero.
3. Modelling the phase change process In constructing a solution to a problem involving a liquid/solid change of phase, a possible approach would be to accurately track the position of the phase boundary [l, 12,131 and then to solve equations (l-5) in the fluid region and the appropriate form of (5) in the solid region.
K. Morgan, Numerical analysis of freezing and melting with convection
277
These solutions would have to satisfy the boundary conditions on temperature, velocity and heat removal at the phase change surface [14]. However, the complexity of front tracking methods in two dimensions led to the adoption in the present work of the alternative approach in which the phase change process is modelled by a variant of the enthalpy method [S]. In this method the phase change is assumed to occur over a temperature range and the associated latent heat effect is handled by increasing suitably the specific heat in this range. Thus if the phase change is assumed to occur over the temperature interval [T,, T,], where T, is the solidus temperature and T, the liquidus temperature, then the specific heat cA used in the calculation is defined by CA = c + [H(T
- Ts)-
H(T
- T,)]L/AT
(6)
where H denotes the Heaviside function, L is the latent heat and AT = T, - T, is the phase change interval. When this method is applied to the analysis of pure materials, in which the phase change occurs at a specified temperature (i.e. T, = T,), a phase change interval AT (# 0) must be assumed. Then it has been demonstrated [7,15,16] that good results can be obtained for problems involving conduction provided that the right choice is made for the value of AT. However, for materials in which the phase change does indeed occur over a reasonable temperature range this problem does not arise and the actual physical values of T, and T, can then be used successfully [17]. The enthalpy method is used here in conjunction with a fixed finite element mesh. Gartling [7] noted that such an approach was well suited to the computation of the fluid motion except for the problem of the application of the no-slip bounda~ condition at the phase boundary. Exact application of this boundary condition of zero-fluid velocity at the phase change surface will not be possible as the fluid/solid interface will, in general, not coincide with the nodes of the finite element mesh. Gartling [7] overcame this problem by using a smearing approach in which the viscosity was greatly increased in the phase change region, but this method was found to be unsuccessful within the context of the explicit solution algorithm to be adopted here. Instead, the problem was overcome by adopting the simpler, but perhaps less accurate, approach of reducing the nodal velocity to zero whenever the nodal temperature lies below the liquidus. 4. Finite element solution procedure The approach outlined in the previous section for dealing with a liquid/solid change of phase was implemented in the computer code CONDIF, developed at J.R.C. Ispra by J. Donea et al. [18] for the analysis of transient convective-conductive heat transfer problems. Space discretisation is performed via the finite element method while time-stepping is achieved by the use of the fractional step method 1181. This method was originally used by Chorin [19] in the finite difference solution of incompressible flow problems. The solution procedure adopted in the CONDIF code has been widely reported in the literature [18,20,21], and a brief description of the salient points will therefore suffice here. The spatial discretisation is accomplished via the use of 4-noded isoparametric e’tements
278
K. Morgan, Numerical analysis of freezing and melting with convection
with associated bilinear velocity and temperature fields and elementwise constant pressure. Application of the standard method of weighted residuals [22] to equations (l)-(5) then yields a matrix equation system of the form
cu=o,
(7)
du
Mz=fa+fd+fb+
cp,
dT Madt=ga+g,+g,.
(9)
Here C is the pressure matrix, M is the mass matrix, MA the capacitance subscripts a, d and b denote the vectors arising from the advection terms, the and the boundary conditions and body force terms respectively. With the solution known everywhere at time tn, the fractional step method determine the solution of (7) and (8) at the next time level t,+* (= t, + At). In intermediate velocity field u* is calculated so as to satisfy only a discretised without the pressure term, i.e., U
*= u”+AtM-‘(f:+f:+f:).
matrix and the diffusion terms is then used to this method an version of (8)
(10)
The required velocity and pressure fields, untl and p”+l, are then obtained by adding to u* the dynamical effect of pressure, determined so as to ensure that the incompressibility condition remains satisfied, i.e., AtC’M-‘Cp”+’ U “+l =
= - C’U* ,
u* + AtM-‘Cp””
.
(11)
(12)
Problems associated with the application of velocity boundary conditions ensure that the matrix C’M-‘C in (11) splits into two independent submatrices, each of which is singular. It has been argued [23] that this may be due to the fact that incorrect pressure equations are produced in the boundary elements if the conditions on prescribed tangential velocity are strictly applied. As a consequence, the boundary conditions on the tangential component of velocity are satisfied by u”+’ in a weak sense only, thus ensuring that the pressure equations have the correct form in the boundary elements and enabling (11) to be solved, upon specification of a single reference value for pressure. With the velocity and pressure at time so determined, the solution is completed by evaluating the temperature at this time from g)‘in the form T “+I= T” + AtMi’(g:
+ g: + g:) .
(13)
In explicit schemes of this type, stability considerations restrict the size of the time step At which may be used, and it is therefore computationally advantageous to replace the consistent
K. Morgan, Numerical analysis of freezing and melting with convection
279
mass and capacitance matrices in equations (lOb(13) by the lumped ~diagona1) matrices obtained by adding all terms of each row of the consistent matrices and placing the result on the diagonal. A detailed analysis of the effects of this approximation and the determination of the resulting stability requirements has been made by Donea et al. [18]. Their results are used here and then, at any stage in the process. The time step At must satisfy the relation
PcIfr12 PIhI IhI ’ 4~
4k
’ uihi
(14)
where L = (hl, hz) is the element side vector. If a temperature dependent specific heat c is used in (6), the lumped capacitance matrix has to be recalculated at each time step. This can be avoided if the specific heat in the solid and liquid regions can be assumed separately constant (though not necessarily equal). Then recomputation is required only for the elements currently undergoing a change of phase, and the computational efficiency of the scheme is improved. The efficiency of the scheme is not affected by the inclusion of a temperature dependent thermal conductivity, k, in (5) since the diffusion matrix gd in equation (9) is already recalculated at each time step.
5. Freezing in a thermal cavity
This finite element scheme has been shown [9] to be capable of providing accurate solutions to phase change problems involving only conduction, but it is difficult to assess the accuracy of its predictions for problems involving convection due to the scarcity of suitable experimental data. However, in a recent paper Sparrow et al. [3] describe results obtained in a series of experiments designed to investigate the effect of natural convection on the process of freezing. In these experiments a liquid is placed at temperature To in a cylindrical container whose ends are thermally insulated and whose outer surface is held fixed at the temperature To. Freezing commences when a cooled inner cylinder at temperature Ti is inserted into the liquid. A detailed analysis of these experiments was not possible as the author did not have access to the full material property data. It was therefore decided to perform calculations in a similar geometry and subject to similar boundary conditions with a view to comparing the general experimental and numerical conclusions. The dimensions of the region analysed and the type of boundary conditions applied are shown diagrammatically in Fig. 1. The region was covered by a mesh of 100 square elements and boundary conditions of zero velocity were applied on the walls of the region at all times. Although somewhat unrealistic physically, the material properties were taken to be constant and equal in both phases and the values adopted for the problem parameters were defined by p=
1.0,
k =
c=
1.0,
p = 0.003 )
Tre3,= cl ,
0.001 , (13
L = 5.0 .
The phase change in the material was assumed to occur over a temperature
range defined by
280
K. Morgan,
Numerical
analysis of freezing
\’
01
Fig. 1. Region
analysed
T, = -0.6,
and temperature
and melting with convection
’
‘\ar
’
’
’
azzo
boundary
T, = -0.65
’
’
’
’
11
conditions applied. Initially T = TO everywhere.
(15b)
and, in the problems to be considered here, such a range should certainly ensure that the enthalpy method takes correct account of the latent heat effect [17]. In the initial calculations, the value of g was adjusted so as to give a Rayleigh number of 3000 when To- Ti = 1, as previous experience with the model [21] had shown that the coarse mesh adopted would then give a good representation of the flow. The computer runs performed are detailed in Table 1. In Run A the liquid is initially at a temperature only just above the liquidus value, and the effects of convection in this case are negligible. The frozen region is then in the form of a perfect cylinder at each instant and this is illustrated in Fig. 2 which shows the progress of the liquidus to time 3000 when the computation ceased. In the remaining Runs B-D the liquid is superheated and the table shows how natural convection slows and ultimately stops the freezing process. Plots of the isotherms and streamlines at steady state conditions are shown in Fig. 3. Close examination of the numerical results produces the same conclusions as those reached by analysing experimental observations [3] viz. increasing To at a fixed value of Ti:,,or decreasing -Ti at a fixed value of To, decreases the
Table 1 Details of freezing calculations Run
TO
Ti
Time for steady state
A B
-0.59
c
0 2
-1 -1 -2 -2
3000’ 648 1578 299
D
0
‘Denotes end of computation conditions not attained.
time; steady state
281
K. Morgan, Numerical analysis of freezing and melting with convection
Run B
Run C
(al Fig. 2. Position of the liquidus at various times in Run
Run II
lb)
Fig. 3. Temperature (a) and stream function,(b) at steady state conditions for Runs B, C and D. The solidified region is shaded.
A.
time at which steady state is attained and decreases the volume of liquid which is eventually frozen. The effect of the magnitude of the latent heat was investigated by repeating Run B firstly with L = 0 and then with L = 10. It was found in all three cases that the distribution of temperature and stream function at steady state was virtually identical but that the time at
RunE
(a)
Fig. 4. Temperature shaded.
Run F
(a) and stream function (b) at steady state conditions
for Runs E and F. The solidified region is
K. Morgan,
282
Numerical
analysis
of freezing
and melting with convection
which steady state conditions were reached increased with latent heat from 414 with L = 0 to 918 with L = 10. In Runs E and F, Run B was repeated with the value of g increased by a factor of ten (Rayleigh number = 30 000) and by a factor of fifty (Rayleigh number = 150 000) respectively. The effect of convection is then significantly increased and a finer mesh was necessary in Run F at the highest Rayleigh number. Again the effect of the increased convection is to reduce the final amount of frozen material and also to reduce the time at which steady state conditions result. The steady state was attained at time 249 in Run E and at time 184 in Run F. The distribution of temperature and stream function at steady state for these two cases is shown in Fig. 4.
6. Melting in a thermal cavity To further demonstrate the capabilities of the model, the inverse problem of melting in a cylindrical thermal cavity was also considered. The geometry of the previous section is used but now the initial temperature of the region is taken to be below the solidus, and melting is initiated by increasing the temperature of the inner wall. The values given in (15) were again Table 2 Details of melting calculations Run
TO
Ti
Time for steady state
G H
-0.66 -1
2 2
1900 2609
RunG
Run H
Run
(al
I
RunJ
Fig. 5. Temperature (a) and stream function (b) at steady state conditions for Runs G, H, I and J. The solidified region is shaded.
K. Morgan, Numerical analysis of freezing and melting with convection
283
used for the thermal properties but T,, was put equal to the liquidus T,. Two basic Runs (G and H) were made as detailed in Table 2 and the effect of increasing the amount of convection was demonstrated by repeating Run H with the value of g increased by a factor of ten (Run I) and by a factor of fifty (Run J). As the effect of convection is increased in this manner, the steady state is reached in a shorter time (in Run I at time 945 and at time 436 in Run J) and the amount of material remaining frozen increases. For this complete set of calculations, the temperature and stream function distributions at steady-state are shown in Fig. 5.
7. Conchdons An explicit finite element method for analysing problems involving melting and freezing in the presence of convection has been described. The method has been shown to be capable of producing results which agree qualitatively with experimental observations, but a detailed analysis of the accuracy of the method is still not possible because of the scarcity of suitable experimental results. The explicit nature of the method means that it should be possible in the future to extend the method to the analysis of three-dimensional phase change problems involving convection.
The work reported in this paper was performed while the author was a Visiting Scientist in the Applied Mechanics Division of the Joint Research Centre of the European Communities, Ispra, Italy. The author wishes to thank the director of the J.R.C. for providing the financial support which made the visit possible and the staff of the Applied Mechanics Division, particularly Dr. J. Donea, for their help and encouragement.
References [l] J. Crank, How to deal with moving boundaries in thermal problems, in: R.W. Lewis, K. Morgan and O.C. Zienkiewicz, eds., Numerical Methods in Heat Transfer (Wiley, Chichester, 1981). [2] E.M. Sparrow, R.R. Schmidt and J.W. Ramsey, Experiments on the role of natural convection in the melting of solids, J. Heat Transfer 100 (1978) 11-16. [3] E.M. Sparrow, J.W. Ramsey and R.G. Kemink, Freezing controlled by natural convection, J. Heat Transfer 101 (1979) 578-584. [4] C.A. Lapuduta and W.K. Mueller, The effect of buoyancy on the formation of a solid deposit freezing on a vertical surface, Internat. J. Heat Mass Transfer 13 (1970) 13-25. [5] E.M. Sparrow, S.V. Patankar and S. Ramadhyani, Analysis of melting in the presence of natural convection in the melt region, J. Heat Transfer 99 (1977) 520-526. [6] J. Szekely and P.S. Chhabra, The effect of natural convection on the shape and movement of the melt-solid interface in the controlled solidi~cation of lead, Metallurg. Trans. (1970) 1195-1203. [7] D.K. Gartling, Finite element analysis of convective heat transfer problems with change of phase, in: K. Morgan, C. Taylor and C.A. Brebbia, eds., Computer Methods in Fluids (Pentech, London, 1980). [8] D.R. Atthey, A finite difference scheme for melting problems, J. Inst. Math. Appl. 13 (1974) 353-366.
284
K. Morgan, Numerical
analysis of freezing and melting with convection
[9] K. Morgan, An explicit method for convection problems involving a change of phase, SMiRT-6, Paris, 1981, Paper B278. [lo] P.M. Gresho, R.L. Lee and R.L. Sani, On the time-dependent solution of the incompressible Navier-Stokes equations in two and three dimensions, in: C. Taylor and K. Morgan, eds., Recent Advances in Numerical Methods in Fluids (Pineridge, Swansea, 1980). [ll] J.C. Heinrich, P.S. Marshall and O.C. Zienkiewicz, Penalty function solution of coupled convective and conductive heat transfer, in: C. Taylor, K. Morgan and CA. Brebbia, eds., Proc. Conf. Numer. Meths. in Laminar and Turbulent Flow (Pentech, London, 1978). [12] G.H. Meyer, The numerical solution of multidimensional Stefan problems-a survey, in: D.G. Wilson, A.D. Solomon and P.T. Boggs, eds., Moving Boundary Problems (Academic Press, New York, 1978). [13] K. O’Neill and D.R. Lynch, A finite element solution of freezing problems using a continuously deforming coordinate system, in: R.W. Lewis, K. Morgan and O.C. Zienkiewicz, eds., Numerical Methods in Heat Transfer (Wiley, Chichester, 1981). [14] H.S. Carslaw and J.C. Jaeger, Conduction of heat in solids (Oxford University Press, London, 2nd ed., 1959). [15] G. Comini, S. Del Guidice, R.W. Lewis and O.C. Zienkiewicz, Finite element solution of non-linear heat conduction equations with special reference to phase change, Internat. J. Numer. Meths. Engrg. 8 (1974) 613-624. [16] K. Morgan, R.W. Lewis and O.C. Zienkiewicz, An improved algorithm for heat conduction problems with phase change, Internat. J. Numer. Meths. Engrg. 13 (1978) 1191-1195. [17] K. Morgan, R.W. Lewis and H.R. Thomas, Finite element modelling of drying stresses in timber and cooling stresses in cast metal, in: M.J. O’Carroll et al., eds., Modelling and Simulation in Practice, Vol. 2 (Emjoc, Northallerton, 1980). [18] J. Donea, S. Giuliani, H. Lava1 and L. Quartapelle, Solution of the unsteady Navier-Stokes equations by a finite element projection method, in: C. Taylor and K. Morgan, eds., Recent Advances in Numerical Methods in Fluids, Vol. 2: Computational Techniques in Transient and Turbulent Flow (Pineridge, Swansea, 1981). [19] A.J. Chorin, Numerical solution of the Navier-Stokes equations, Math. Comput. 22 (1968) 745-762. [20] J. Donea and S. Giuliani, The computer code CONDIF for transient convective-conductive heat transfer, Ispra Report EUR 6822/I EN, 1980. [21] J. Donea, S. Giuliani and H. Lava], Accurate explicit finite element schemes for convective