Numerical analysis of water melting and solidification in the interior of tubes

TECHNICAL PAPERS

Numerical analysis of water melting and solidification in the interior of tubes

S. I. S. de Souza^I; H. A. Vielmo^II

^IUniversidade Regional Integrada do Alto Uruguai e das Missões – URI; Campus de Santo Ângelo; Rua Universidade das Missões, 464; 98802-470 Santo Ângelo, RS. Brazil; [e-mail]

^IIUniversidade Federal do Rio Grande do Sul; Departamento de Engenharia Mecânica; Rua Sarmento Leite, 425; 90050-170 Porto Alegre, RS. Brazil; [e-mail]

ABSTRACT

Latent energy storage systems find applications in many engineering fields, including industrial refrigeration plants, air conditioning installations, recovery of heat in industrial processes, etc. To tackle the design of such systems, it is necessary to have correlations to account for the heat transfer during the melting and solidification of the phase change material (PCM). This work describes and analyzes the results obtained from the numerical simulation of pure water melting and solidification in the interior of tubes, which are typically present in ice banks of air conditioning systems. The shown results consider natural convection, accounting for the inversion in the water density. In the melting process, the considered initial conditions followed the classical Stefan and Neumann approach. The presented simulation results include the evolution of the phase change interface, and of the temperature, density and streamlines fields. Correlations for the Nusselt number and for the melted material volume as functions of time have been proposed.

Keywords: Phase change, melting and solidification, ice banks, finite volumes, polar geometry

Introduction

The increasing energy demand and the requirements for an environmentally sustained development have been major motivating factors for the implementation of more efficient means of energy conservation. In addition to conservation, the rational use of energy has been the focus of many studies, so that the processes are more efficient and less expensive, reducing waste, preserving the environment and the nonrenewable energy sources, and avoiding unnecessary expenses. In this sense, the employment of thermal energy storage has the objective of making the energy consume more rational. This storage is performed in three basic, distinct ways: sensible heat, latent heat and chemical. The systems that operate with sensible heat require high temperature differences with respect to the environmental, as well as a storage substance with large mass. This leads to some disadvantages, to mention the need of large spaces to accommodate the system and thick, expensive thermal insulation. The storage systems that use the latent heat present a superior amount of volumetric capacity, besides requiring smaller temperature differences with respect to the environment, which decreases the thermal insulation expenses.

According to the "ASHRAE Applications (1995)", the use of thermal storage systems may be an economically attractive approach for heating or cooling loads when they are of short time duration and occur cyclically; when there is an incompatibility with the available energy source, or its supply is limited; when the energy costs are time-dependent; and when economical incentives are provided for the use of load-shifting equipment.

In the refrigeration field, ice banks have been employed as an energy storage solution to reduce the electrical power required during the hours of peaks in the demand, shifting the consume curve to a more favorable situation. It must be noted that in many countries, the cost of industrial electricity is time dependent during the day, having higher values in the periods of high demand. It can be also mentioned the contractual demand, resulting into a fine to the user in case the energy consumption is exceeded. Aiming at economical feasibility, it has been increasingly attempted to develop systems that employ the storage of latent heat, despite its disadvantages. Among those, it can be mentioned the inherent irreversibilities due to heat transfer, the decrease in the efficiency of the thermodynamic cycle due to the decrease in the evaporation temperature in the air-conditioning installations, the increase in the physical space that is necessary to accommodate the ice bank, and uncertainties in the evaluation of the heat transfer between the PCM and the system working fluid.

Water in the solid phase has been employed as the storage substance due to its properties and characteristics: high melting latent heat (334 kJ/kg), melting temperature of 0°C at the atmospheric pressure, accessibility, low cost, no pollutant effects, and chemical stability. As the most important disadvantage characteristic, it can be pointed the supercooling phenomenon (a metastable equilibrium) during the solidification, "Chen and Lee (1998)". The ice bank consists of thermally insulated tanks that are inserted into a refrigeration loop. In some systems, the ice is formed on the outside of the tubes. In such systems, a limit must be impose to the ice thickness formed on the piping, since if it exceeds a given value there will be a decrease in the heat transfer during the discharge process, caused by the ice thickness and the flow obstruction. Examples of works that investigated the phase change occurring in the outside of the tubes are "Ho and Chen (1986)", "Zhang and Faghri (1996)", "Abugderah and Ismail (2000)", "Stampa et al. (2002)".

Encapsulating the ice in tubes is one of the most employed examples of latent heat storage, consisting of a tank with a bank of horizontal tubes, across which flows a fluid that works as a thermal device that promotes the charge or discharge of the system. Among the works that deal with the process of melting and solidification in the interior of tubes, it can be pointed "Rieger et al. (1983)", "Ho and Viskanta (1984)", and "Rieger and Beer (1986)". The latter is one of the first experimental and numerical studies on phase change that employ water as the PCM. Another investigation worth mentioning is that of "Cao et al. (1991)". The encapsulations in spherical, rectangular and spiral configurations are other possibilities.

Nomenclature

A = area, m²

B = body force, m/s²

c = specific heat, J/kg

E = coefficient of the proposed function for the Kirchhoff temperature, coefficient of the Nusselt number correlation

Fo = Fourier number ,a t/r²₀

g = gravitational field, m/s²

h = enthalpy, J/kg; convection heat transfer coefficient, W/(m²K)

k = thermal conductivity, W/(mK)

L = representative grid size

= average thickness of the liquid ring, m

m = mass, kg

Nu = Nusselt number,

p = pressure, N/m²

q = heat transfer rate, W

r = radial coordinate, mesh relation, radius, m

Rac = Rayleigh number,

Ste = Stefan number,C_p(T_p - T*)/l

T = temperature, °C, K

t = time, s

u = tangential velocity, m/s

V = volume, m³

v = radial velocity, m/s

Greek Symbols

a = thermal diffusivity, m²/s

b = thermal expansion coefficient, 1/K

f = generic variable

h = 9.297173x10^-6, °C^-b

l = melting and solidification latent heat, J/kg

µ = dynamic viscosity, Ns/m²

n = kinematic viscosity, m²/s

q = angular coordinate, rad

r = density, kg/m³

Subscripts

c = relative to convection

conv relative to convective process

f = relative to melting

H = relative to hydrostatic component

ini = relative to initial

L = relative to liquid

max = relative to maximum

mist = relative to mixture

p = relative to constant pressure

r = relative to radial direction

ref = relative to reference

s = relative to solid

S = relative to solid

w = relative to wall

0 = relative to tube wall or to initial

q = relative to angular direction, angular position

Superscripts

* = melting point

b = 1.894816, constant to "Eq. 31"

Physical Model and Boundary Conditions

The focus of this work is a latent "cold storage" modulus, which consists of a bank of tubes containing pure water inside, subjected to a crossed external flow. In the present approach only one tube is considered through a symmetric, two-dimensional model, since such tubes are typically subjected to small temperature variations along its length.

The boundary conditions are now established in order to facilitate the physical model comprehension. For the thermal problem the first-kind boundary conditions are adopted, "Eq (1)", both for the melting and solidification, as seen in "Fig.1".

In the melting process T_w>T^*, and in the solidification process T_w<T^*, where T_w is the tube wall temperature and T^* is the melting temperature (in this work T^*=0 °C).

For the hydrodynamic problem, the non-slip and non-permeability boundary conditions are imposed on the surfaces of the tube and of the solid PCM,

where u and v are, respectively, the angular and radial velocities, and r₀ and r_s are the radius of the tube and of the PCM, as seen in "Fig. 1". The initial condition adopted for the thermal problem during the melting is dependent on the approach that is applied. In the so called Stefan problem,

while in the Neumann problem,

In the solidification the initial boundary condition, for the thermal problem, is:

where T(r,q) is the temperature inside the tube, T_ini is the initial temperature and the subscripts S and L stand for the solid and liquid regions, t is time, r and q are respectively the radial and angular coordinates.

For the hydrodynamic problem, the initial conditions for the melting and solidification are:

The non-slip and impermeability boundary conditions are adopted during all the situations.

Governing Equations

The heat transfer from or to the encapsulated fluid originates the thermal gradients, which in turn cause the existence of gradients in the density in the interior of the fluid domain. Therefore, thermal convection may be present due to the fluid flow caused by buoyancy. For the solution of such a problem, it is necessary the knowledge of the velocity, temperature, density and pressure fields, which characterizes the coupling of the hydrodynamic and thermal problems. The equations that govern the phenomenon for the considered geometry are present below in radial coordinates. The continuity equation is given by:

where r is the density and t is de time. The momentum equations in the r and q directions are:

in which p is the pressure, and B_r and B_q are the radial and angular components of the gravitational field acceleration g. The dynamic viscosity is represented by µ. The Boussinesq approximation is adopted for the buoyancy terms in the liquid phase. The references for the temperature and density are taken from the melting point: T_ref=T^* and r_ref=r^*. From the Maxwell thermodynamic relations, one can set the state equation that relates the density to temperature and pressure for a given substance:

For small variation in the density, it is acceptable to expand the above function in the vicinity of a point of reference by means of the Taylor series to give:

which is valid for liquids. Using now the definition of the thermal expansion coefficient b, "Eq. " becomes:

inserting "Eq. " into "Eq. (8)" and "Eq. (9)" results:

where p_H is the pressure that incorporates the hydrostatic component. Since the density of the water is not a linear function of the temperature, "Fig. 2". The behavior of the volumetric expansion coefficient, with respect to the temperature, was based on the expression of "Vasseur et al. (1983)", for the water density, resulting in the following expression:

In problems of phase change, the domain can be divided into three distinct regions: solid, liquid and interface, where each region has its own characteristics. For the pure substances, the temperature is constant in the liquid-solid interface, and so it does not provide any information on the evolution of the process. To overcome this difficulty, it is considered the energy equation in the enthalpy form. Applying the model proposed by "Cao et al. (1989)", in which it is employed the concept of the Kirchhoff Temperature, described and applied in "Vielmo (1993)", results:

where h is the specific enthalpy. The remaining terms are described below:

in which c_p is the specific heat, l is the latent heat, k is the thermal conductivity.

Numerical Methodology

For the integration of the governing equations, the finite volumes method is applied, as described in "Patankar (1980)" and "Maliska (2004)".

In the equation, the diffusion coefficient depends whether the considered region is the liquid, interface or solid. The system of equations resulting from the integration of the differential equations is solved by the TDMA scheme, with a block correction. The adopted interpolation function is the Power-Law, and the coupling between the pressure and velocity fields is based on the SIMPLEC algorithm. A computational mesh formed by 40x40 volumes was employed after grid independence tests were performed. The velocities are stored in the interfaces of the control volumes. In most of the solutions, relaxation factor of 0.8 for all the variables proved to be adequate to deal with the non-linearities of the problem. Values in the order of 10^-8 were used for convergence criteria, applied to relative error of each variable between two consecutive iterations, as well as to the residual masses in each control volume. In order to obtain stable, implicit transient solutions, the time steps were varied and chosen so that the enthalpy, in the liquid domain, changed about 200 J/kg each time step (no control volumes submitted to the phase change). When there were control volumes close to the phase change region (the biggest period of time), the time steps were reduced to values that caused variation in the enthalpy in the range of approximately 10 J/kg in the liquid.

Validation of the Proposed Methodology

For the validation, the results obtained in this work were compared to the numerical and experimental results presented in "Rieger and Beer (1986)", "Fig. 3". The solid is kept fixed in the center of the tube. Figure 4 presents the solution obtained in this work for the same situation considered in the aforementioned work. As seen, there exists an excellent agreement between the two solutions.

For verification of the independence of the solutions with respect to the discretization, it was applied the Grid Convergence Index (GCI), according to "Roache (1998)", and the usual procedure for estimation and reporting of discretization error in CFD applications. The melting was simulated using the Stefan approach for the initial condition T_ini,S=0 °C, for the boundary condition T_w=10 °C and for a tube with a diameter r₀=0.032 m. The following equations were initially applied to evaluate the apparent order p of the employed method:

where L is the representative grid size, the DA_i is the area of the i^th volume, and N is the total number of volumes used. The f terms are the solutions, and the subscripts 1, 2 and 3 refer from the more refined to the less refined mesh. The three different meshes tested have 40×40, 120×120 and 360×360 control volumes. Such an approach was necessary to allow geometrical agreement between the evaluation points. The temporal agreement was also observed. The enthalpy was the variable used for the inspection of the apparent order p of the method, since this property is the key variable in the problem. The obtained results are presented in "Tab. 1".

Such evaluations were carried out for three angular positions in the tube, and in each position three locations were chosen along the radius. The angles of the positions are referred to with respect to the "Fig. 1". These positions were chosen in such a way to cover the upper, intermediary and lower regions of the tube.

For the evaluation of the relative error, e, and of the GCI, "Eq's. (21) and (22)" were employed considering the apparent order p of local approximation, "Eq. (20)":

As seen in the "Tab. (1)", the maximum relative error, e, reached values lesser than 1 % in the comparison of the results obtained for the meshes 40x40 and 360x360. The CGI presents a maximum error of approximately 1.5 % for the meshes 40x40 and 120x120. One can be state, based on the results presented in "Tab. (1)", that the employed grid of 40x40 control volumes proved to be adequate for the analysis presented in this work. Also it is possible to observe the appearance of a negative value for the apparent order p, indicating an oscillatory convergence in that position (close to the center of the tube). In this position low values for the relative errors also occur. In rare situations these variations, when close to zero, are indicative that the exact solution was reached. Anyway, it does not represent unsatisfactory calculations, or parameter to invalidate the results.

In all the obtained solution, for the different meshes, it was not registered any variation in the pattern of the variables fields. The verification of the presence or not of multiple solutions was investigated through alterations in the directions of the TDMA scheme. As observed, changing the directions in all directions did not cause any modification in the solutions.

Results and Discussion

The Nusselt Number

It is first necessary to define an expression for the Nusselt number in the melting of the PCM encapsulated in the interior of the tube, and for that matter it was adopted the enclosure treatment. For the geometry considered in this work, the solution of the steady-state, one-dimensional diffusion equation, with no source term, provides the below temperature profile:

The heat transfer by diffusion, q_r , in a given position r is given by:

The convective heat transfer q_conv is written using the general formulation below:

where h_c is the heat transfer coefficient. In a situation where the PCM is partially liquid and solid, but before the convection has started, the heat rate across r₀ is therefore transferred by diffusion to the liquid up to the phase change boundary. In this situation, it can be written:

The Nusselt number is then defined for this problem as the relation q_conv/q_r0 . For the situation of pure diffusion:

The situation described by "Eq. (27)" occurs in the initial instants of the melting process. For the situation in which the natural convection takes place, all the energy that crosses the domain by diffusion is carried to the interior of the domain by convection, and the temperature profile in the liquid is dependent of the time, r and q coordinates. Employing the concept of thin wall together with "Eqs. (24) and (26)", for a two-dimensional situation, the Nusselt number can be written for a given angular position q, and time t:

Inserting the dimensionless temperature profile, T_adm, into "Eq. (28)", the local Nusselt number becomes:

The Nusselt number, thus defined, represents a relation between the existing convective heat rate and the heat rate that would exist in case the process were purely diffusive. All the previous equations provide results for a control volume in a given angular position q, and time t. In engineering, the main interest usually is the average behavior of the phenomenon. Therefore, it is interesting to determine the average behavior of the Nusselt number, which is done in the following way:

Results for the Melting Process

Next, it is presented the solutions for the wall boundary conditions of T_w=4 ºC and T_w=8 ºC, applied to a tube with radius r₀=0.032 m, and starting at the initial condition in the solid of T_ini,S=T^*=0 ºC. The thermal resistance imposed by the tube wall was neglected. Since the phenomenon is transient, solutions were collected at various instants t for the two wall boundary conditions. The streamlines in the liquid region are plotted in the same graphs of the density. This allows a better understanding of the phenomenon associated with the inversion in the density, based on the observation of the local Nusselt number, "Eq. (29)", and of the average Nusselt number, "Eq. (30)".

Firstly, "Fig. 5" shows the solutions for T_w=4 ºC. In the figure, the behavior of the local Nusselt number as a function of the angular position is presented. In the initial instants, as for t=100 s, it can be noted that the phenomenon is purely diffusive, indicating an equality between the heat rates, as described in "Eq. (27)".

For t=3000 s, there are variation of the Nusselt number with the angle. It can be observed the beginning of the convective effects through a small increase on the Nusselt number for the angular positions in the upper region, 0°<q<100°. For the lower region, 100°<q<180°, there is a decrease in the value of the Nusselt number. The maximum temperature that the liquid reached in this simulation was 4 °C, while the minimum was 0 °C. The portion of liquid with the highest temperature will assume the highest density. Due to the gravitational field action, the fluid with the lowest density and temperature will be concentrated in the upper positions, and so intensifying the heat rate in that region. The fluid with the highest temperature and density will occupy the lower positions of the tube, leading to a decrease in the heat rate from the wall to the fluid. As a result, the Nusselt number presents values lower than the unity. This decrease occurs due to the assumptions taken for the derivation of "Eq. (26)", in which it was considered that the heat rate that penetrates in the boundary by diffusion is mostly transported by advection in the interior of the tube. In "Fig. 6", it can be observed the evolution of the phase change boundary as well as the behavior of the process parameters, temperature and density, which validate the aforementioned explanations. The temperature evolution demonstrates the concentration of liquid with higher temperatures in the lower regions of the domain. The aforementioned characteristics of the process intensify with the time evolution; that is, the Nusselt number increases in the upper positions and decreases in the lower positions.

Figure 5 shows the existence of two peaks in the Nusselt number for t=7000 s. For this same instant, "Fig. 6" shows that the vortex structure is complex: two clockwise (negative) vortices, while another vortex spins in the anticlockwise direction. This causes the formation of two locations with peaks in the Nusselt number. This vortex structure is also registered for t=12000 s. In the vortex boundaries, there is a larger mass displacement, which extends from the region close to the solid PCM to the tube wall, originating regions with high thermal gradients, as observed on the left side of the figures, where the temperature field is presented. In initial instants, as for t=5000 s, the first peak in the Nusselt number is due to the formation of a vortex in the upper part of the tube, and a second peak occurs in a intermediary angular position, and is formed by the PCM ring thickness variation in the liquid state. For t=18000 s, there is another vortex structure. The vortex that spins in the clockwise direction, a negative vortex, predominates, while a positive vortex, located in the lower region of the domain, does not affect significantly the heat transfer process. In this situation, there is a peak in the Nusselt number in the upper part of the tube, a region where liquid with lower temperature and density is concentrated.

Figure 7 shows the behavior of the Nusselt number along the angular position for T_w=8 °C, while "Fig. 8" presents the behaviors of the temperature and the density fields, to which streamlines are juxtaposed, for this situation. Again, the diffusive process is observed in the initial instants, t=100 s, in "Fig. 7". After some time, the convection mechanism takes place, as can be observed by the behavior of the Nusselt number. In this situation, one can observe the inversion in the behavior of the Nusselt number as compared to the previous situation. The peaks in the Nusselt number occur for angular positions q>50º. For q<50º, there is a decrease in the Nusselt number.

The difference between the two solutions can be explained by the inversion in the density for temperatures close to 4 °C. For this wall temperature, there is concentration of liquid with high temperature and lower density in the upper regions. The fluid formed from the melting of solids with lower temperatures also presents a density inferior to the maximum value, and is prone to move to the upper regions. However, as the fluid starts the motion, it mixes with fluids at higher temperature, reaching the temperature of 4 °C, a situation in which the water has the maximum density, and so it moves towards the lower positions of the tube. For the previous situation, T_w=4 °C, there was an inverse behavior, although with the same driving mechanism.

The behaviors of the density and temperature fields for T_w=8 ºC can be observed in "Fig. 8". The simulations indicate that, for situations in which T_w<4 ºC, the behaviors are similar to those found for T_w=4 ºC, while for T_w>8 ºC it is similar to T_w=8 ºC. For T_w=6 ºC, in the beginning of the process the behavior resembles the situation for T_w=4 ºC; after some time, the solution presents an inversion, with a behavior similar to that for T_w=8 ºC.

"Figs. 5, 6, 7 and 8" show the dynamic aspect of the heat transfer phenomenon, which makes it unfeasible to draw conclusions from local results to predict the heat transfer behavior. Instead, it is more adequate to consider the average behavior, as is in "Fig. 9", which shows the average Nusselt number, "Eq. (30)", for the two considered situations.

In the search for correlations that enable the evaluation of the heat transfer, and that are not prone to be affected by the "singular" behavior of the water density, it was employed an expression that has been presented in "Sasaguchi et al. (1998)". Applying it into the dimensionless momentum equations, more specifically in the buoyancy terms, leads to the following Rayleigh number:

where h = 9.297173×10^-6, is the average thickness of the liquid ring, b = 1.894816, a is the thremal diffusivity, n is the kinematic viscosity, and T_ref takes the value of T_w in the melting.

This parameter allowed the identification of a pattern behavior of the Nu_average/Ra_c ratio, as a function of the melted volume V_f to the initial volume V₀ ratio, V_f/V₀, and also as a function of Fo, as can be observed in "Fig. 10", where logarithm scales were adopted.

Such behaviors were approximated by the following relations:

where a, b, c and d are the coefficients to be determined. This allows writing an expression for the Nusselt number in the following way:

The terms E₁ and E₂ are defined as:

where the a is the thermal diffusivity.

Based on the behaviors of the presented solutions, it is possible to establish correlations for the melted volume V_f and, consequently, for the average thickness of the melted PCM:

where the terms d and c are obtained according to:

The efficacy of the proposed correlations can be evaluated in "Fig. 11", where they are compared to the results obtained with the simulation. In the values of the simulations, the percentage variations were added, allowing the visualization of the maximum existing differences between the two results. It can be observed for the Nusselt number that the maximum differences occur at the beginning of the melting process, where the diffusion is the predominant heat transfer mode. In the time periods where the convection is the dominant mode, the differences between the two solutions decrease. Indeed, in literature, convection correlations with uncertainties as high as 20% occur.

Results for the Melting Process Followed by Re-solidification

As a consequence of the variable ambient conditions, which cause variations in the thermal load, it is common the partial use of the stored energy in the ice banks along the time. Therefore, simulations considering this situation also must be carried through, allowing a better knowledge of the heat transfer mechanisms. The results for the partial melting process followed by resolidification are shown in "Figs. 12 and 13", considering that the PCM in the solid phase is initially supercooled (Neumann approach), T_ini,S=-5 ºC. For t=0 s, the tube wall was exposed to T_w=8 ºC, initiating the melting. After a period of 7760 s, the surface temperature on the external surface changes to T_w=-5.5 ºC, initiating the resolidification. The behaviors of the enthalpy and density in the domain are also presented. The dashed line represents the enthalpy in the region where the liquid reaches its maximum density, for a temperature of 4 ºC. On the left corner on the bottom, the Rayleigh number for the considered instant is presented. Together with the relative density, the streamlines are presented, enabling the visualization of the liquid behavior during the phenomenon evolution. The instants of time for which the simulations were run are presented on the bottom corner.

The behavior of the melting process is first shown for the time instant t=200 s. It can be observed that the sensible energy associated with the solid supercooling rapidly disappears, and in this instant diffusion predominates. For t=3000 s, there is no longer traces of sensible energy, and the advection mechanism is already active in the heat transfer. For t=7780 s, the boundary condition is altered so that solidification occurs. One can note a thin layer of solidified material in the external region. The liquid is then confined between the two phase change boundaries. The sensible energy associated with the liquid is released to the two boundaries. This initially reduces the advance of the solidification in the external boundary. In the inner boundary, the melting continues until the thermal equilibrium is reached. Initially, the thickness of the solidified external layer is not uniform, as for when t=7800 s. This occurs due to the fact that it is in the upper region where the liquid with the lowest density and enthalpy level concentrates. This difference is small, being the more intense the larger is the temperature T_w set to cause the melting. In the solidification, the vortices structure changes drastically, as can be observed for the instants t=7780 s, t=7800 s, and t=8000 s. As the solidification advances, the inner liquid and the solid reach thermal equilibrium, as for when t=8750 s. After this, the buoyancy effects become negligible, the velocities associated with the liquid are of such a magnitude that convection stops actuating, so that diffusion predominates again, and the liquid density becomes homogeneous. The inner solid becomes "insulated", without altering its energetic state until t=10750 s, when the two boundaries meet. From this point ahead, there is heat transfer from the inner solid by diffusion. This can be visualized for time instants later than t=12500 s.

In the solidification, the time and intensity of the natural convection actuation are larger for situation where the initial sensible energy, or equivalently the initial liquid temperature (T_ini,L), is large and when the T_w temperature is low. To evaluate the influence of the natural convection in the solidification process, it was plotted the behavior of the T_mist, defined in "Eq. 42" with m_L being the liquid mass.

In both situations, it was considered the same physical model, and the same boundary and initial conditions, the difference being solely on the solution method. The first case was dealt considering the solidification process to be purely diffusive; in the other case, the convection mechanism was taken into account.

The observations show a faster decrease on T_mist in the simulation that took the natural convection into account, and this trend remained the same for the other simulations having different boundary and initial conditions. The difference between the solutions narrowed with the decrease of the liquid sensible energy in the beginning of the process. The time needed for the temperature homogenization in the interior of the liquid domain also decreased for the same reason; the decrease on T_ini,L leads to the same effects.

The influences of the natural convection, and of the sensible energy existing in the beginning of the process on the PCM solidification time, can be seen in "Fig. 14" (right side). The results for the two different solidification situations are exposed. For the first and second cases, it was set T_ini,L=8 °C and 16 °C. For both cases, T_w=-5.5 °C and the solutions were obtained for processes with pure diffusion and with natural convection.

Observing the left size of "Fig. 14", one realizes the proximity between the purely diffusive and natural convection solutions, since the steeper decrease on T_mist when convection is taken into account is not sufficient to cause noticeable differences between them (solidification time). This behavior can be justified by the sensible energy being much smaller than the latent energy (that is, low Stefan numbers). The differences in the values of T_mist /T_ini,L between the two observed cases in the left side are consequences of the fluid movement in the convective process, causing a faster homogenization of the temperature in the liquid. The amount of sensible energy in the liquid has small influence on the solidification speed, for temperatures around T_ini,L=8 °C, so that there is a coincidence between the solutions with natural convection and pure diffusion.

For higher temperatures, as for T_ini,L=16 °C, the solidification speed is smaller for the initial instants in the convective process, V_f/V₀>0.6. After the temperature homogenization is reached, for V_f/V₀>0.4, there is an inversion in the solidification speed, and the convective process becomes faster. In the end of the process, the solidified volume is identical for the two cases, and the time needed for this to occur, t_max, was similar for all the simulations.

Conclusions

Melting initiates by purely diffusive heat transfer. In the evolution of the process, the space occupied by the liquid phase increases. This causes the increase in the Rayleigh number, intensifying the natural convection, which in many cases can become the dominant heat transfer mechanism. The solidification process initiates in the presence of natural convection, since the larger portion of the domain is in the liquid phase. As time progresses, the thermal gradients and the spaces occupied by the liquid diminish, decreasing the Rayleigh number and causing, from a certain point in time, the diffusion to become the dominant heat transfer mechanism, as seen in "Fig. 13" for Rac=63. As expected, the dominant mechanism of heat transfer in the interior of the duct changes with time and the boundary conditions.

The vortices structure, which varies with time and the boundary conditions, drastically interferes in the change of the heat transfer mechanism. The degree of complexity of the heat transfer phenomenon, in the considered problems, is associated with the "singular" behavior of the water density. This behavior is responsible for the variation in the vortices structure and in the isothermals in the interior of the tube. In the melting process, this structure "absorbs" the energy on the wall and then "releases" it to the solid, intensifying the exchange mechanism in specific positions, where the fluid with low enthalpy concentrates. The vortices structures have influence on the phase change boundary. As long as the diffusive process is dominant, the solid maintains a cylindrical configuration; however, as the vortices structure intensifies, the geometry of the phase change boundary "deforms". In the initial instants of the resolidification, the liquid sensible energy is released to the two phase change boundaries. The external phase change boundary maintains a homogeneous geometry, for the liquid sensible energy is rapidly consumed, establishing a homogeneous density and minimizing the effect of natural convection.

In the solidification, the assumption of pure diffusion leads to results that are close to the solutions that consider natural convection, for the engineering oriented cases that were considered in this work. In the melting, a larger discrepancy between the two solutions was observed.

Correlations for the Nusselt number, and for the melted volume as functions of time were also obtained for practical use in engineering.

Paper accepted May, 2005. Technical Editor: Clóvis R. Maliska.

• Abugderah, M. M., and Ismail, K. A. R., 2000, "Low Temperature Applications of a Phase Change Thermal Storage System Performance", VIII ENCIT/MERCOFRIO 2000, Porto Alegre, RS, Brazil.
• ASHRAE Handbook, 1995, "HVAC Applications", Atlanta, GA, USA.
• Cao, Y., Faghri, A. and Chang, W. S., 1989, "A Numerical Analysis of Stefan Problems for Generalized Multi-Dimensional Phase-Change Structures Using the Enthalpy Transforming Model, Int. J. Heat Mass Transfer", vol. 32, n. 5, pp. 1289-1298.
• Cao, Y., Faghri, A. and Juhasz, A., 1991, "A PCM/Forced Convection Conjugate Transient Analysis of Energy Storage Systems With Annular and Countercurrent Flows", Journal of Heat Transfer, vol. 113, pp. 37-42.
• Chen, Sih-Li and Lee Tzong-Shing, 1998, "A Study of Supercooling Phenomenon an Freezing Probability of Water Inside Horizontal Cylinders", Int. J. of Heat and Mass Transfer, vol. 41, n° 4-5, pp. 769783.
• Ho, C. J. and Chen, S., 1986, "Numerical Simulation of Melting of Ice Around a Horizontal Cylinder", Int. J. of Heat and Mass Transfer, vol. 29, n° 9, pp. 1359-1369.
• Ho, C. J. and Viskanta, R., 1984, "Heat Transfer During Inward Melting in a Horizontal Tube", Int. J. Heat Mass Transfer, vol. 27, n. 5, pp. 705-716.
• Maliska, C. R., 2004, "Transferência de Calor e Mecânica dos Fluidos Computacional", 2a Ed, Livros Técnicos e Científicos Editora SA, Rio de Janeiro, RJ, Brasil.
• Patankar, S. V., 1980, "Numerical Heat Transfer and Fluid Flow", McGraw-Hill, New York, USA.
• Rieger, H. and Beer, H., 1986, "The Melting Process of Ice Inside a Horizontal Cylinder: Effects of Density Anomaly", Transactions of the ASME, vol. 108, pp. 166-173.
• Rieger, H., Projahn, U., Bareiss, M. and Beer, H., 1983, "Heat Transfer During Melting Inside a Horizontal Tube", Transactions of the ASME, vol. 105, pp. 226-234.
• Roache, P. J., 1998, "Verification and Validation in Computational Science and Engineering", Hermosa Publishers, New Mexico, USA.
• Sasaguchi, K., Kuwabara, K., Kusano, K. and Kitagawa, H., 1998, "Transient Cooling of Water Around a Cylinder in a Rectangular Cavity A Numerical Analysis of the Effect of the Position of the Cylinder", Int. Journal of Heat Transfer, vol. 41, pp. 3149-3156.
• Stampa, C. S., Nieckele, A. O. and Braga, S. L., 2002, "Um Estudo Numérico do Crescimento da Camada de Gelo do Lado Externo de um Tubo Vertical", IX Congresso Brasileiro de Engenharia e Ciências Térmicas, Paper CIT02-0595.
• Vasseur, P., Robillard L. and Shekar, B. C., 1983, "Natural Convection Heat Transfer of Water Within a Horizontal Cylindrical Annulus With Density Inversion Effects", Journal of Heat Transfer, vol. 105, pp. 117-123.
• Vielmo, H. A., 1993, "Simulação Numérica da Transferência de Calor e Massa na Solidificação de Ligas Binárias", Doctorate Dissertation, Federal University of Santa Catarina, Florianópolis, SC., Brasil, 193 p.
• Zhang, Y. and Faghri, A., 1996, "Semi-Analytical Solution of Thermal Energy Storage System with Conjugate Laminar Forced Convection", Int. J. of Heat and Mass Transfer, vol. 39, n° 4, pp. 717724.

Publication Dates

• Publication in this collection
  31 Aug 2005
• Date of issue
  June 2005