Numerical Investigation of PCM Melting in a Finned Tube Thermal Storage

Abstract

Due to their high energy storage capacity, latent heat storage units using phase change materials (PCMs) have gained considerable attention over the past three decades. The heat exchange of a PCM with the surrounding medium is managed by the thermal energy equation (solidification/melting) with different complex boundary and initial conditions. In this study, we propose to solve numerically this equation applied to a PCM by the finite difference method. To understand the storage phenomenon of solar energy in the form of latent heat in PCM, initially found under cooling at 18°C, we studied the fusion in a specific configuration corresponding to a tubular exchanger with five circular horizontal fins. In this perspective, we propose in this work a numerical investigation based on an enthalpy formulation to study the melting of a PCM in a finned heat exchanger. This numerical approach gives simultaneously the temperature distributions in the PCM storage system and temporal propagation of the melting front during the melting of the PCM when it is exposed to a hot airflow. Also, we give in this study the transient evolution of the longitudinal air temperature profiles.

1. Introduction

Every latent heat thermal energy storage (LHTES) should include three main components: an appropriate PCM in the required temperature range, a container for the stocking substance, and an appropriate carrying fluid for an effective heat transfer from the heat source to the heat storage. Moreover, the PCMs require an important heat exchange area because of its low thermal conductivity.

One method is to increase the heat transfer surface area by employing finned surfaces [1, 2, 3]. Many numerical and analytical models of PCM solidification and PCM melting in finned PCM-air exchangers have been published in order to evaluate their performance.

Rostamizadeh et al. [4] developed a numerical model of energy storage in a rectangular container of PCM based on an enthalpy formulation, and the effect of PCM thickness on temperature distribution and melting fraction was investigated. Thus, these researchers established that 5 mm is the best thickness of the PCMs. Besides, the study results show that the PCM mass and the melting time verify a linear relationship. Ismail et al. [5] analyzed a numerical and experimental study on the solidification and the melting of PCM around a vertical axially finned isothermal cylinder. The model is based upon the pure conduction mechanism of heat transfer. From the given results, it can be noticed that the fin thickness does not have any important influence on the time of solidification; indeed, the time for full solidification and the solidification rate are strongly affected by the fin length and the number of fins. The difference of temperature has a reverse impact on the solidification of PCM, and the full solidification time tends to decrease when the temperature difference increases. Seck et al. [6] studied the evolution of the melting front of a plate of paraffin (52–54) submitted to sunshine rays in order to design a less voluminous and more efficient heat storage system. The results have shown that the natural convection plays a major role in the kinetics of front propagation. Hence, to numerically reproduce the temporal evolution of the interface throughout the fusion process, the developed natural convective phenomenon should necessarily be taken into account during the melting process. El Qarnia et al. [7] have presented a numerical model of PCM melting in a rectangular cavity heated with three protruding heat sources mounted on a vertical conducting plate. They have carried out numerical investigations for studying the influence of different key parameters on the thermal performance of the PCM-based heat sink. They have developed two correlations, one for the appropriate melt fraction and the other for the maximum working time. The proposed approach can be useful for the design of PCM-based cooling systems.

The present paper presents a numerical investigation of the melting of PCM (paraffin RT27) in a finned tube thermal storage for air conditioning systems taking into account the presence of natural convection. The first part includes the presentation of a numerical model based on the continuity, momentum, and thermal energy equations, treated by the finite volume method. The main objective of this numerical approach is to investigate the temporal evolution of PCM melting and the melting front with the evolution of the liquid fraction, as well as the longitudinal profiles of the heat transfer fluid (HTF) in the duct.

2. Computational domain and mathematical formulation

In the present paper, the studied configuration is a PCM-air heat exchanger with three coaxial cylinders (Figure 1). The radius of the inner cylinder is R₀ = 8 cm, that of the middle is R₁ = 13 cm, and that of the outer is R₃ = 17 cm. The thickness of all tubes is 3 mm, and the length of the system is 1.2 m. As shown in Figure 1, the space between the inner and middle cylinders is filled with PCM (RT27). In the first-pass flow, the heat transfer fluid (HTF) passes through the central tube, and in the second pass-flow, it flows in the reversal direction in the annular space before exiting the PCM-air heat exchanger module by the same input side. The mass flow rate of air is 0.08 kg.s⁻¹. The fins are made of aluminum with a constant thickness of 3 mm. The thermophysical properties of the solid and liquid PCM are listed in Table 1.

Figure 1.

Studied configurations and computational domain.

		RT27
Melting interval [°C]		[27–27,5]
Latent heat [kJ.kg^−1]		179
Density [kg.m^−3]	Solid	870
	Liquid	760
Specific heat [kJ.kg^−1.°C^−1]	Solid	2.4
	Liquid	1.8
Thermal conductivity [W.m^−1.°C^−1]	Solid	0.24
	Liquid	0.15
Thermal volumetric expansion coefficient [kg.m^−3.°C^−1]		8.5 × 10^−4
Kinematic viscosity [m^2.s^−1]		1.5 × 10^−6

Table 1.

Thermophysical properties of the PCM RT27 [8].

The problem is 2D since it has a symmetry of revolution, and hence the computational domain (the hatched area in Figure 1) becomes 5 cm in width H and 1.2 m in length L.

2.1. Modeling of the heat transfer to the airflow circling in the axial duct and the annular space

The heat transfer to the HTF-governing equation is given as follows:

$${\mathrm{\rho }}_{\mathrm{air}}{\mathrm{Cp}}_{\mathrm{air}}\left(\frac{\partial {\mathrm{T}}_{\mathrm{air}}}{\partial \mathrm{t}}+w_{\mathit{air}}\frac{\partial {\mathrm{T}}_{\mathrm{air}}}{\partial \mathrm{z}}\right)={\mathrm{S}}_{\mathrm{T}}$$E1

For the first passage of air,

$${\mathrm{S}}_{\mathrm{T}}=\frac{2{\mathrm{h}}_{\mathrm{i}}\left({\mathrm{T}}_{\mathrm{PCM}}-{\mathrm{T}}_{\mathrm{air}}\right)}{{\mathrm{R}}_0}/{\mathrm{h}}_{\mathrm{i}}=166\mathrm{W}{\mathrm{m}}^{-2°}{\mathrm{C}}^{-1}$$E2

For the second passage of air,

$${\mathrm{S}}_{\mathrm{T}}=\frac{2{\mathrm{R}}_1{\mathrm{h}}_o\left({\mathrm{T}}_{\mathrm{PCM}}-{\mathrm{T}}_{\mathrm{air}}\right)}{{\mathrm{R}}_2^2-{\mathrm{R}}_1^2}/{\mathrm{h}}_{\mathrm{o}}=114\mathrm{W}{\mathrm{m}}^{-2°}{\mathrm{C}}^{-1}$$E3

2.2. Numerical modeling of heat transfer in the PCM domain

The natural convection of a transient laminar PCM flow in the PCM container has been studied numerically. Considering constant properties, except the density difference (Boussinesq approximation), and using the enthalpy method, the continuity, momentum, and thermal energy equations for the PCM can be written as follows:

2.2.1. Continuity equation

$$\frac{1}{\mathrm{r}}\frac{\partial \left(\mathrm{r}\mathrm{u}\right)}{\partial \mathrm{r}}+\frac{\partial \mathrm{w}}{\partial \mathrm{z}}=0$$E4

2.2.2. Momentum equations

The momentum equations incorporate sink terms to take into account the changing phase of the PCM. The condition that all velocities in solid regions are zero is considered using an enthalpy-porosity approach [9].

$$\frac{\mathrm{\partial U}}{\mathrm{\partial t}}=-\mathrm{div}\left(\overrightarrow{\mathrm{V}}\mathrm{U}-\mathrm{\nu }\mathrm{gra}\overrightarrow{\mathrm{d}}\mathrm{U}\right)-\frac{1}{\rho }\frac{\mathrm{\partial P}}{\mathrm{\partial r}}-\mathrm{\nu }\frac{1}{r}\left(\frac{2\partial \mathrm{V}}{\mathrm{r}\partial \mathrm{\theta }}+\frac{\mathrm{U}}{\mathrm{r}}\right)+\frac{{\mathrm{V}}^2}{\mathrm{r}}-\frac{1}{\mathrm{\rho }}\mathrm{c}\frac{{\left(1-\mathrm{f}\right)}^2}{{\mathrm{f}}^3+\mathrm{b}}\mathrm{U}$$E5

$$\frac{\mathrm{\partial W}}{\mathrm{\partial t}}=-\mathrm{div}\left(\overrightarrow{\mathrm{V}}\mathrm{W}-\mathrm{\nu }\mathrm{gra}\overrightarrow{\mathrm{d}}\mathrm{W}\right)-\frac{1}{\rho }\frac{\mathrm{\partial P}}{\mathrm{\partial z}}+\mathrm{g}\mathrm{\beta }\left(\mathrm{T}-{\mathrm{T}}_0\right)-\frac{1}{\rho }\mathrm{c}\frac{{\left(1-\mathrm{f}\right)}^2}{{\mathrm{f}}^3+\mathrm{b}}\mathrm{W}$$E6

where S_u is constituted by two terms: the first of which corresponds to the melting sink term, and the second is the thermo-convective generation term. Furthermore, S_w is the extinction of “W” component due to the melting of PCM.

2.2.3. Energy equation

$$\rho \left(\mathrm{T}\right){\mathrm{C}}_{\mathrm{p}}\left(\mathrm{T}\right)\frac{\mathrm{\partial T}}{\mathrm{\partial t}}=-\mathrm{div}\left(\rho \left(\mathrm{T}\right){\mathrm{C}}_p\left(\mathrm{T}\right)\overrightarrow{\mathrm{V}}\mathrm{T}-\mathrm{\lambda }\left(\mathrm{T}\right)\mathrm{g}\overrightarrow{\mathrm{r}}\mathrm{ad}\mathrm{T}\right)$$E7

where [10]

$$\rho \left(\mathrm{T}\right)\mathrm{Cp}\left(\mathrm{T}\right)=\left\{\begin{array}{l}{\mathrm{\rho }}_{\mathrm{s}}{\mathrm{Cp}}_{\mathrm{s}}\mathrm{T}\le {\mathrm{T}}_{\mathrm{sd}}\\ {\rho }_{\mathrm{l}}{\mathrm{Cp}}_{\mathrm{l}}\mathrm{T}\ge {\mathrm{T}}_{\mathrm{liq}}\\ \frac{{\rho }_{\mathrm{mix}}\mathrm{L}}{{\mathrm{T}}_{\mathrm{liq}}-{\mathrm{T}}_{\mathrm{sd}}}{\mathrm{T}}_{\mathrm{sd}}\le \mathrm{T}\le {\mathrm{T}}_{\mathrm{liq}}/{\rho }_{\mathrm{mix}}=\mathrm{f}{\rho }_{\mathrm{l}}+\left(1-\mathrm{f}\right){\mathrm{\rho }}_{\mathrm{s}}\end{array}\right.$$E8

$$\mathrm{\lambda }\left(\mathrm{T}\right)=\left\{\begin{array}{l}{\mathrm{\lambda }}_{\mathrm{s}}\mathrm{T}\le {\mathrm{T}}_{\mathrm{sd}}\\ {\mathrm{\lambda }}_l\mathrm{T}\ge {\mathrm{T}}_{\mathrm{liq}}\\ \frac{{\mathrm{\lambda }}_{\mathrm{l}}-{\mathrm{\lambda }}_{\mathrm{s}}}{{\mathrm{T}}_{\mathrm{liq}}-{\mathrm{T}}_{\mathrm{sd}}}\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{liq}}\right)+{\mathrm{\lambda }}_{\mathrm{l}}{\mathrm{T}}_{\mathrm{sd}}\le \mathrm{T}\le {\mathrm{T}}_{\mathrm{liq}}\end{array}\right.$$E9

The liquid fraction “f” can be analyzed by the following equation:

$$\mathrm{f}=\frac{\mathrm{T}-{\mathrm{T}}_{\mathrm{sd}}}{{\mathrm{T}}_{\mathrm{liq}}-{\mathrm{T}}_{\mathrm{sd}}}$$E10

2.3. Initial and boundary conditions

For the velocity components, the nonslip boundary conditions on all solid walls are imposed.

The thermal limiting condition imposed on the interfaces between the PCM envelope and the air evolving in the tubular exchanger results from the application of the continuity rule. Thus, a conducto-convective mixed condition is applied on both sides of the PCM container. During the first passage of air, the heat transfer between the PCM and the HTF through the wall of the PCM container is written as follows:

$${\mathrm{h}}_{\mathrm{i}}\left({\mathrm{T}}_{\mathrm{PCM}}-{\mathrm{T}}_{\mathrm{air}}\left(\mathrm{z}\right)\right)=\lambda {\left.\frac{\mathrm{\partial T}}{\mathrm{\partial r}}\right|}_{\mathrm{r}={\mathrm{R}}_0}$$E11

During the second air passage,

$${\mathrm{h}}_{\mathrm{o}}\left({\mathrm{T}}_{\mathrm{PCM}}-{\mathrm{T}}_{\mathrm{air}}\left(\mathrm{z}\right)\right)=-\lambda {\left.\frac{\mathrm{\partial T}}{\mathrm{\partial r}}\right|}_{\mathrm{r}={\mathrm{R}}_1}$$E12

As initial conditions, it is assumed in the charging mode (melting) that the storage unit and the air in the duct are initially at a uniform temperature of 45°C and the temperature of the inlet air is 18°C.

2.4. Resolution method

In the present study, a finite volume method for the numerical solution of 2D unsteady natural convection flow and energy Eqs. (5)–(7) is used [1]. Thus, the domain including the PCM is divided into a convenient number of control volumes (NR = 40; NZ = 300).

The staggered grid method is employed to accurately control mass conservation and heat transfer on all control volumes, and an exponentially fitted spatial discretization scheme is used. Furthermore, the equation is solved by the alternating direction implicit (ADI) technique, in which the Crank-Nicolson scheme is used for the time discretization, and the SIMPLER algorithm for the treatment of the pressure-velocity coupling.

Numerical results are obtained by developing a specific code using the FORTRAN language.

3. Validation of the numerical approach

The performance of the proposed method is verified with the experimental data performed by Longeon et al. [11]. This was done by simulating, under the same operating conditions, the melting process of the paraffin RT35 used as PCM, with the thermophysical properties as listed in Table 2.

Melting point [°C]		35
Latent heat [kJ.kg^−1]		157
Density [kg.m^−3]	Solid	880
	Liquid	760
Specific heat [kJ.kg^−1.K^−1]	Solid	1.8
	Liquid	2.4
Thermal conductivity [W.m^−1.K^−1]		0.2
Kinematic viscosity [m^2.s^−1]		3.3 × 10^–6

Table 2.

Thermophysical properties of RT35 [11].

The PCM storage system is composed of two concentric cylinders (Figure 2): the first one with an inner diameter of 44 mm is made of Plexiglas, and the other with an inner diameter of 15 mm is made of stainless steel. The length of the whole system is 400 mm. The HTF flows into the inner tube, and the PCM fills up the annular space. As it is shown in Figure 2, the thermocouples are distributed on section B.

Figure 2.

Scheme of the experimental loop and thermocouple positions in the test section [11].

Figure 3.

Comparison of the experimental results and our numerical solution.

It is found that the temperature evolution does not depend on angular position, indicating that the heat exchange can be modeled with a 2D approach. In Figure 3, which gives the temperature evolution with time at two radial positions (a and b), our numerical results are compared with those determined experimentally by Longeon et al. [11]. A good agreement between numerical and experimental data can be noted. In fact, in both positions, the same shape tendency with a rapid decrease of the temperature at the beginning followed by a slow rating decrease characterizing the phase change process can be noted. The average deviation between our numerical simulation and the experimental data of Longeon et al. [11] is approximately 5%. This deviation is very acceptable for this type of study [12].

4. Results

4.1. Temporal evolution of temperature and velocity fields

Figures 4 and 5 present the distinctive temporal evolution of the distribution of the temperature and current lines. The analysis of the thermal fields has shown that the temperature is higher in the inner part near the central tube and the fins. The peripheral part is also warm, but to a lesser degree, because this surface is exposed to less warm air, since it has lost some of its sensible heat during the first passage. In the early stages of the fusion, very little free space exists for the flow. Hence, the homologous mappings given the current lines and the isotherms have revealed that the liquefaction front is parallel to the tubes and the fins, which confirms the conduction predominance.

Figure 4.

Temporal evolution of the temperature field and with corresponding current lines.

Figure 5.

Temporal evolution of temperature fields in heat storage mode (five circular fins).

As a result of the free convection movements which have become increasingly strong in molten paraffin, the volume is mainly extended over the entire height of the warmer central tube. While mounting along the air-heated tube, the molten paraffin progressively accumulates the heat of the central tube and comes against the fins before diverting toward the interior of the domain. Then, it descends along the solid-liquid interface by gradually releasing the accumulated heat, till reaching the bottom of the liquefied PCM cavity. The same phenomenon is observed on the temperature profiles near the outer tube. The observation of the temperature evolution in the figure reveals that the conduction phenomenon predominates on the opposite side of the tubes.

For all heat storage times, a comparison of the temperature levels in the zones located between two consecutive internal fins shows higher temperature levels in the upper zones.

Figure 6 shows the longitudinal profiles of the air evolving in the exchanger. In storage mode, the air entering the exchanger at 45°C is cooled in the first passage and then in the second passage by providing its PCM-sensible heat, initially at 18°C. This gradually leads to:

• preheating the PCM to 18°C at its melting point (27°C);

• liquefying the PCM in the temperature range 27–27.5°C; and

• overheating the liquefied PCM.

Figure 6.

Temporal evolution of the longitudinal air profile in the duct.

These exchange natures are presented in Figure 7, in which the temperature evolution values of the air leaving the exchanger as a function of time are plotted. Indeed, in the time interval between 5 and 40 min, an almost constant output temperature can be observed, indicating a heat supply that is equally constant during the phase of state change.

Figure 7.

Temporal evolution of the air temperature at the outlet air duct.

4.2. Temporal evolution of the melting front

Figure 8 shows the temporal evolution of the melting front during the energy storage. After a very short time (about 1 min), we have already observed that the fusion initiates at the wall of the central tube, heated to a temperature higher than that of the material liquefaction by the HTF.

Figure 8.

Temporal evolution of the melting front of the storage mode.

During the PCM melting, three successive transfer regimes can be noted in the transmission mode of the heat: in the first moments, the zone occupied by the liquid phase has a very small width, so there is practically no convection heat transfer. Consequently, the liquid-solid interface moves parallel to the tubes and the fins under the unique action of conduction. At this stage, the morphology of the melting/fusion front is presented as a flat surface. Then, when the width of the liquid phase gives birth to convection, the regime is called mixed and corresponds to a combination of natural convection and conduction. This results in a higher melting rate in the solid part; the shape of the interface is thus modified. Finally, when the liquid phase has become wide enough to allow a full development of the convective movements, the convective regime is reached, in which the main part of the heat transfer inside the liquid phase is due to the natural convection. The solid-liquid interface then is slightly deformed, and the natural convection will be the main motor for the enhancement of the melting front inside the PCM enclosure.

These transfer regimes are clearly seen in Figure 9, in which the evolution of the liquid fraction is shown, indicating the melted volume according to the total volume, as a function of time. Initially, there was a low growth rate corresponding to the conductive regime. Between roughly 20 min and 30 min, the growth rate becomes more important, underlining the contribution of thermo-convective transfer. Starting from 35 min, the liquefaction rate becomes more important, highlighting the dominance of the thermo-convective effects. The total liquefaction ends after 45 min.

Figure 9.

Transient evolution of liquid fraction.

5. Conclusion

In the present work, we have studied the melting of the PCM (Paraffin RT27) in a tubular PCM-air heat exchanger. We have concluded that in the storage mode, the complete PCM melting time was found to be substantially shorter (about one-fourth of the time) compared to that recorded for solidification during the destocking [10]. Indeed, in contrast to what happens during solidification involving an increase in the thermal resistance to the wall due to the formation of a solid PCM layer, the initiation of liquefaction on the side of the hot walls promotes the formation of thermo-convective loops which are intensified during the fusion.