Numerical Simulation Of Ice Melting Using Engineering Essay

Published: Last Edited:

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

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 food-stuffs 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 Ref. [1]. Owing to the nonlinear form of the thermal energy balance equation at the unknown location of this time-dependent 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 Ref. [2]. The Lie-group shooting method was used in Ref. [3]. In addition, the heat balance integral method (HBIM) was used in Ref. [4]. The asymptotic solution for zero-phase type problem for cases having small but finite Biot numbers was used in Ref. [5]. On the other hand, numerical methods are used commonly. The explicit Finite Deference Method was used in Ref. [6]. The Faedo-Galerkin Finite Element Method with a Crank-Nicolson time scheme was used in Ref. [7]. Reference [8] investigated the simulation of ice formation in one and two dimension using the cell-centered Finite Volume Method based on the latent heat source approach. Reference [9] proposed 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. Reference [10] have described and compared 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. In addition, commercial software packages have been used to solve the 'moving boundary' problem. Reference [11] has used MATLAB to model both; the melting of snow/ice and the effect of adding salt to the ice/snow. Reference [12] investigated numerically the three-dimensional melting of ice around a liquid-carrying tube placed in an adiabatic rectangular cavity using PHOENICS Code.

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. Reference [13] 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 multi-dimensional problems due to the difficulties with algorithms of implementations and large computational cost. Reference [14] 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 zero-width 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 temporal 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 is initially at a uniform temperature, Fig. (1-a). This body is suddenly subjected with a constant temperature and , Fig. (1-b). It is required to find: The temperature distribution of the solid phase (Ice), and liquid phase (Water), , the Melt fraction, , the Liquid-solid Interface location, , and its velocity, . To solve this classical two-phase Stefan problem, the following assumptions are considered Ref. [15, 16]:

Fig. 1 Mathematical Model

Thermal conductivity of Solid Thermal conductivity of Liquid.

Thermal Diffusivity of Solid Thermal Diffusivity of Liquid.

Specific heat of Solid Specific heat of Liquid.

Neglect the volume variation for both phases by assuming a constant and equal density for both phases.

The thermal properties of a single phase are independent of the temperature.

The initial temperature is constant.

A planer and sharp surface (interface) is separating the phases at a constant phase change temperature .

The Latent heat of fusion is constant.

Neglect the surface tension and curvature effects at the interface.

No convection or radiation heat transfer at the boundaries, only pure and isotropic conduction is considered, and gravity effect is neglected.

No internal heat generation.

The governing equation with its Initial and Boundary Conditions are:

(3)

(4)

(5)

(6)

The energy balance at the interface yields;

The general solution of the temperature distribution for both the liquid and solid phases is as follows:-

Where is the root of the following transcendental equation:-

And the interface velocity is as follows;

And the melt fraction is as follows;

Where L is a finite length of the domain needed to validate the numerical results of the One-Dimensional model.

FINITE VOLUME 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 for Two-Dimensional domain. A general nodal point is identified by P and its neighbors. In one-Dimensional domain, the nodes to the west and east of P are identified by W and E respectively. The west side boundary of the control volume is referred to by 'w' and the east side of the control volume is referred by 'e', Fig. 2. In Two-Dimensional domain, the nodes to the west, east, south and north of P are identified by W, E, S and N respectively. The side boundaries of the control volume is referred to by; 'w', 'e', 's' and 'n' for the west, east, south, and north sides respectively, Fig. 3. The time domain is divided into a number of time steps of size . Variables at the previous time level are indicated by the superscript (o, Old). In contrast, the variables at a new time level are not superscripted Ref. [17].

Fig. 2 A typical One-Dimensional control volume (shaded area)

Fig. 3 A typical Two-Dimensional control volume (shaded area)

Discretization of One-Dimensional model

A 0.1-m slab has initial temperature . The boundary condition at one end of the slab is subjected to a constant temperature while the other end is kept at the value of the initial temperature as illustrated in Fig. 4. Uniform control volumes with size are used.

The law of energy conservation is combined with the Fourier's law to give a local form as:

where are the enthalpy, temperature and thermal conductivity respectively.

Single phase Formulation:

Substitute (14) into (13), yields;

Multiply (15) by, , then integrate over the control volume faces, yields;

Using fully implicit scheme;

Two-phase Formulation:

Substitute (18) into (13), yields;

Multiply (19) by, , then integrate over the control volume faces, yields;

Using fully implicit scheme;

Each of the boundary conditions is substituted; into (17) for discretizing the Single Phase One-Dimensional model and into (21) for discretizing the Two-Phase One-Dimensional model respectively. Then divide each of the resulting equation by and rearrange yields:

Where are listed in Table I.

Fig. 4 One-Dimensional Model

TABLE I GLOBAL DISCRETIZATION OF ONE-DIMENSIONAL MODEL.

0

0

0

0

0

0

East Boundary

0

Discretization of Two-Dimensional model

Consider a region with area is initially at a uniform temperature. 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 and are used. The governing Two-Dimensional energy equation is:

Single phase Formulation:

Substitute (14) into (23), yields;

Multiply (24) by, , then integrate over the control volume faces, yields;

Using fully implicit scheme;

Two-phase Formulation:

Substitute (18) into (23), yields;

Multiply (27) by, , then integrate over the control volume faces, yields;

Using fully implicit scheme;

Each of the boundary conditions is substituted; into (26) for discretizing the Single Phase Two-Dimensional model and into (29) for discretizing the Two-Phase Two-Dimensional model respectively. Then divide each of the resulting equation by and rearrange yields:

Where are listed in Table II.

Solver

Equation (22) and (30) are solved using the Tri-Diagonal Matrix Algorithm TDMA, Ref. [17].

TABLE II GLOBAL DISCRETIZATION OF TWO-DIMENSIONAL MODEL.

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

Corner D

0

0

0

Fig. 5 Two-Dimensional Model

RESULTS AND DISCUSSION

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 shown in TABLE III. A constant phase change (Melting) temperature . The Latent heat for fusion is constant kJ/kg. 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 equivalent to with a time step of was obtained with a maximum error of melt fraction of less than 3% as shown in FIG. 6.

Fig. 7 shows a maximum percentage error less than 0.78% when the numerical results of temperature distribution is compared with the analytical results at time intervals of 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.

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 shown in TABLE III. A constant phase change (Melting) temperature. The Latent heat for fusion is constant kJ/kg. 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 are used to model the domain with a time step of

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 One-dimensional 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 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.

TABLE III ONE DIMENSION MODEL MATERIAL PROPERTIES.

Liquid

k (W/m.K)

2.22

0.556

C (kJ/kg.K)

1.762

4.22

1000

1000

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 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.

Fig. 6 One-Dimensional GIT Test.

Fig. 7 One dimensional temperature profiles at 1, 5, 20 and 40 min.

Fig. 8 One Dimensional interface location as function of time.

Fig. 9 One Dimension melt fraction as function of time.

Fig. 10 One Dimension interface velocity as function of time.

Fig. 11 Predicted temperature profiles at the diagonal line BD compared with the One-dimension analytical solution at 1, 5, 20 and 40 min.

Fig. 12 Predicted velocity for a particle at the interface moving along the diagonal line BD compared with the One-dimension analytical solution at 1, 5, 20 and 40 min.

1 min.

Fig. 14 Two dimension temperature profile.

5 min.

Fig. 13 Two dimension melt fraction as function of time.

20 min.

40 min.

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.

ACKNOWLEDGMENT

The author would like to thank his wife Dr. Maysoon Basheer Abid (Lecturer at The University of Baghdad, College of Engineering, Water Resources Engineering) for her infinite support and valuable advice.