Numerical Simulation of Ice Melting Using the Finite Volume Method

ABSTRUCT The Aim of this paper is to investigate numerically the simulation of ice melting in one and two dimension using the cell-centered finite volume method. The mathematical model is based on the heat conduction equation associated with a fixed grid, latent heat source approach. The fully implicit time scheme is selected to represent the time discretization. The ice conductivity is chosen to be the value of the approximated conductivity at the interface between adjacent ice and water control volumes. The predicted temperature distribution, percentage melt fraction, interface location and its velocity is compared with those obtained from the exact analytical solution. A good agreement is obtained when comparing the numerical results of one dimensional temperature distribution with the analytical results


INTRODUCTION
The classical Stefan problem is a mathematical model for 'phase-change' or 'moving-boundary' engineering problems involving heat and mass transfer in materials undergoing phase change during thermal processes.Such kind of processes can be the solidification of pure metals, melting of ice, freezing of water and deep freezing of foodstuffs and so on.These materials are assumed to undergo a phase change with a continuously moving liquid-solid isothermal front 'moving boundary' that has to be tracked as part of the solution, (Carslaw and Jaeger, 1959).Owing to the nonlinear form of the thermal energy balance equation at the unknown location of this timedependent liquid-solid isothermal front 'the phase change constant temperature interface' associated with both, the absorption or release of latent heat and the abrupt discontinuity of properties.Analytical solutions are difficult to obtain except for a limited number of special cases.Therefore, approximate analytical methods are used in the analysis of phase change problems.The Picard's iterative method was used in (Witula et al., 2011).The Lie-Group shooting method was used in (Chein-Shan Liu, 2011).In addition, the heat balance integral method (HBIM) was used in (Myers, 2007).The asymptotic solution for zerophase type problem for cases having small but finite Biot numbers was used in (Naaktgeboren, 2007).On the other hand, numerical methods are used commonly.The explicit Finite Deference Method was used in (Savovic and Caldwell, 2003).The Faedo-Galerkin Finite Element Method with a Crank-Nicolson time scheme was used in (Rincon and Scardua, 2008).The numerical simulation of ice formation in one and two dimension using the cell-centered Finite Volume Method based on the latent heat source approach was investigated in (Prapainop and Maneeratana, 2004).An enthalpy transforming method with finite volume approach to investigate numerically the steady state natural convection melting process of a two-dimensional square ice cavity was adopted in (Al-Zubaidy, 2006).Several effective modeling techniques with approximate analytical and numerical methods using; Enthalpy Method, Boundary Immobilization Method (BIM), Perturbation Method, Nodal Integral Method and the (HBIM) for the solution of One-Dimensional Stefan problem have been described and compared in (Caldwell and Kwan, 2004).In addition, commercial software packages have been used to solve the 'moving boundary' problem.The melting of snow/ice and the effect of adding salt to the ice/snow has been modeled using MATLAB in (Patrick et al., 2008).The three-dimensional melting of ice around a liquid-carrying tube placed in an adiabatic rectangular cavity was investigated numerically using PHOENICS Code in (Sugawara et al., 2011).
Basically, there are two main approaches for the solution of the Stefan problem.One is the front-tracking method, where the position of the phase boundary is continuously tracked in every time step.Sadoun et al. (2012), has used the Goodman (HBIM) which explicitly tracks the motion of the isothermal liquid-solid front with an explicit finite difference method including modified boundary immobilization scheme in transformed coordinates based on front-tracking and referred to as variable space grid method (VSGM).Hence, this method does not require the discontinuity approximation for isothermal problems and it is poorly suited to multidimensional problems due to the difficulties with algorithms of implementations and large computational cost.Zhaochun et al. (2011) presented the moving interface locations and used these location coordinates as the grid points to find out the arrival time of moving interface respectively.Through this approach, the difficulty in mesh generation was avoided completely.Another approach is to use a fixed-grid method, which implicitly contains the moving interface condition within the mathematical model.This method is more flexible than the front-tracking and is suitable for multi-dimensional problems.For the second method, the latent heat is accounted for by either the temperature-based or the enthalpy-based method.
The temperature-based approach considers the temperature as the only dependent variable in the governing equation.In order to avoid the discontinuity in isothermal problems, an approximate numerical smoothing must be used and a special integration is needed to compute the latent heat.On the other hand, the enthalpy-based method is further divided into basic enthalpy, apparent heat capacity and latent heat source.In the basic enthalpy scheme, enthalpy is used as the primary variable and the temperature is calculated from a previously defined enthalpy-temperature relation.This method gives reasonably accurate results for metallic solidifying over mushy ranges, but it is complex and computationally expensive and it performs poorly for isothermal problems.In the apparent heat capacity method, the latent heat is calculated from the integration of heat capacity with respect to temperature.As the relationship between heat capacity and temperature in isothermals involves sudden changes, the zerowidth phase change interval must be approximated by a narrow range of phase change temperature.Thus, the size of time step must be small enough that this temperature range is captured in the calculation.In the latent heat source method, the latent heat is included in the source term of the governing differential equation, which is obtained from a prescribed relation between latent enthalpy and temperature.
The aim of this paper is to predict the temperature distribution, the percentage melt fraction, the liquid-solid interface location and its velocity in one and two dimension using the Finite Volume Method based on the fixed grid, latent heat source approach.The fully implicit time scheme is selected to represent the time discretization.The ice conductivity is chosen to be the value of the approximated conductivity at the liquid-solid interface between adjacent ice and water control volumes.To validate the predicted numerical results, it will be compared with those obtained from the exact analytical solution.

MATHEMATICAL MODEL
Consider a semi-infinite ice body with This body is suddenly subjected ( ) with a constant temperature o T and ( ) .To solve this classical two-phase Stefan problem, the following assumptions are considered, Ozisik (1993) and Alexiades and Solomon (1993): 1.The thermal properties of solid are not equal to that of liquid.2. Neglect the volume variation for both phases by assuming a constant and equal density ( ) ρ for both phases.3. The thermal properties of a single phase are independent of the temperature.
4. The initial temperature is constant ( ) 6.The Latent heat of fusion ( ) SL h is constant.7. Neglect the surface tension and curvature effects at the interface.
8. No convection or radiation heat transfer at the boundaries, only pure and isotropic conduction is considered, and gravity effect is neglected.9.No internal heat generation.

ONE DIMENSIONAL CLOSED FORM SOLUTION
The governing equation with its Initial and Boundary Conditions are: ) The energy balance at the interface yields; The general solution of the temperature distribution for both the liquid and solid phases is given as follows, Carslaw and Jaeger (1959): Where ( ) δ is the root of the following transcendental equation:- ( ) And the melt fraction ( ) r M is as follows; ( ) Where L is a finite length of the domain needed to validate the numerical results of the One-Dimensional model.

Grid Generation
In the Finite Volume Method, the first step is to divide the domain into a number of discrete control volumes; ( ) for One-Dimensional domain and ( )

Two Phase Formulation ( )
m T T ≥

:
The enthalpy method-latent heat source approach is used to predict the solution of the Two Phase problem, as follows: Substitute Eq. ( 18) into Eq.( 13), yields; Multiply Eq. ( 19) by ( ) dt dV . ,( ) and then integrate over the control volume faces, yields;  ( ) Each of the boundary conditions is substituted; into Eq.( 17) for discretizing the Single Phase One-Dimensional model and into Eq.( 21) for discretizing the Two-Phase One-Dimensional model respectively.Then divide each of the resulting equation by t ∆ and rearrange yields: Where ( ) are listed in Table 1.

Discretization of Two Dimensional Model (Case Study 1)
Consider a . This region is suddenly subjected at the boundaries with a constant temperature . Due to summitry, only one-fourth of the total area (shaded area) is modeled as shown in Fig. 5 with no heat flux across the surfaces AD and CD.Uniform control volumes with size x ∆ and y ∆ are used.The governing Two-Dimensional energy equation is:

:
Substitute Eq. ( 14) into Eq.( 23), yields; Multiply Eq. ( 24) by ( ) The enthalpy method-latent heat source approach is used to predict the solution of the Two Phase problem, as follows: Substitute Eq. ( 18) into Eq.( 23), yields; Multiply Eq. ( 27) by ( ) dt dV . ,where ( ) , and then integrate over the control volume faces, yields; Each of the boundary conditions is substituted; into Eq.( 26) for discretizing the Single Phase Two-Dimensional model and into Eq.( 29) for discretizing the Two-Phase Two-Dimensional model respectively.Then divide each of the resulting equation by t ∆ and rearrange yields: Where ( ) are listed in Table 2.

One Dimensional Model
The One-Dimensional model with the chosen values of initial and boundary conditions is shown in Fig. 4. Selected material properties for solid and liquid phases are given in Table 3 .The value of Stefan Number equal (0.25).As the convection of the liquid across cell faces and the induced stress due to the expansion of the ice are not included in the present work, the density of the ice is approximated to be equal to that of water to ensure the conservation of mass for each control volume.
The first step to validate the numerical solution is to choose the grid and time sizes that are adequate for obtaining the minimum error of results.This task is accomplished through a Gird Independency Test GIT.A grid size of with a time step of .sec 1 = ∆t was obtained with a maximum error of melt fraction of less than 3% as shown in Fig. 6.Fig. 7 shows a good agreement when the numerical results of temperature distribution is compared with the analytical results at time intervals of .min 40 20 , 5 , 1 and t = using a fully implicit time scheme.The ice conductivity is chosen to be the value of the approximated conductivity at the interface between adjacent ice and water control volumes.The temperature distributions show that ice cells heat up faster due to a high diffusivity and a high temperature gradient of almost a linear shape.Once a control volume is melted, its temperature increase will be slow.This clearly illustrates that the rate of heat transfer is predominantly controlled by the position of the melting front as illustrated in Fig. 8 which indicates the percentage of melt fraction as shown in Fig. 9.As ice is a good insulator, the melting front advances at a decelerating rate as shown in Fig. 10.

Case Study 1:Two Dimensional Model
The Two-Dimensional model with the chosen values of initial and boundary conditions is shown in Fig. 5. Selected material properties for solid and liquid phases are given in Table 3 To validate the numerical solution of the Two-Dimensional model, the length of the diagonal line BD is assumed to be equal to that of the Onedimensional domain (L=0.1 m).Fig. 11 shows the comparisons between the predicted temperature distribution along the diagonal line BD and the One-Dimensional analytical solution with a fully implicit time scheme and the ice conductivity is chosen to be the value of the approximated conductivity at the interface between adjacent ice and water control volumes at time values of .min 40 20 , 5 , 1 and t = It is clear that the predicted temperature distribution along the diagonal line BD exceeds the One-Dimensional analytical solution due to the area of the exposed corner of the Two-Dimensional domain is approximately twice of that for the One-Dimensional domain.This clearly illustrates that the rate of heat transfer is predominantly controlled by the position of the melting front with a decelerating rate at the beginning of melting process then ending with an accelerated rate due to a decrease in the exposed area of the solid phase of the Two-Dimensional domain, as shown in Fig. 12.The percentage of melt fraction is shown in Fig. 13.Fig. 14 shows the Two-dimensional temperature distribution of simulation with a fully implicit time scheme and the ice conductivity is chosen to be the value of the approximated conductivity at the interface between adjacent ice and water control volumes at time values of .min 40 20 , 5 , 1 and t = The temperature distribution clearly show heat transfer through both exposed edges, causing a parabolic temperature profile instead of a linear one as in the One-Dimensional temperature distribution.

Case Study 2:Two Dimensional Model
A Two-Dimensional model of an industrial ice block is simulated with a developed version of the computer program using a body fitted coordinate system.It is subjected to a convection heat transfer boundary condition.Due to summitry, only onefourth of the total area (shaded area) is modeled as shown in Fig. 15.The value of the free convection heat transfer coefficient between ambient air and outer surface of the ice block is ( ) , and the value of the forced convection heat transfer coefficient between the inner surface of the ice block and the flowing fluid is ( ) Holman (1992).The .The Latent heat for fusion is constant kg kJ h SL / 338 = .As the induced stress due to the expansion of the ice is not included in the present work, the density of the ice is approximated to be equal to that of water to ensure the conservation of mass for each control volume.A grid size of ( 20X40) is used to model the domain with a time step of .sec 1 = ∆t as shown in Fig. 16.
The percentage of melt fraction is shown in Fig. 17.
Fig. 18 shows the Two-dimensional temperature distribution of simulation with a fully implicit time scheme and the ice conductivity is chosen to be the value of the approximated conductivity at the interface between adjacent ice and water control volumes at time values of .sec 12 10 , The temperature distribution clearly show that ice melting of the Two-Dimensional model (case study 2) is faster than the model of (case study 1).

CONCLUSIONS
The modeling of ice melting in one and two dimension using the cell-centered finite volume method of a fully implicit time scheme associated with a fixed grid, latent heat source approach is successfully performed.As the melting front should be of a one control volume thickness, this dominating restriction controls both the time interval and the grid sizes.
The temperature distributions show that ice cells heat up faster with a temperature gradient of almost a linear shape.Once an ice cell is melted, its temperature increase will be slow.This clearly illustrates that the rate of heat transfer is predominantly controlled by the position of the melting front which advances at a decelerating rate at the beginning of melting process then ending with an accelerated rate due to a decrease in the exposed area of the solid phase of the Two-Dimensional domain.
. (1).It is required to find: The temperature distribution of the solid phase (Ice), Fig. 4. Uniform control volumes with size x ∆ are used.The governing One-Dimensional energy equation is:

.
The value of Stefan Number equal (0.25).As the convection of the liquid across cell faces and the induced stress due to the expansion of the ice are not included in the present work, the density of the ice is approximated to be equal to that of water to ensure the conservation of mass for each control volume.A grid size of ( to model the domain with a time step of .sec 1 = ∆t

Fig. 12
Fig. 12 Predicted velocity (Case Study 1) Fig. 13 Case Study 1: Two dimensional melt for a particle at the interface moving along fraction as function of time.the diagonal line BD compared with the one dimensional analytical solution as function of time.

Table 1 Global discretization of one dimensional model.
ρInternal nodesw