Available online at www.sciencedirect.com

AASRI Procedía

ELSEVIER

AASRI Procedía 2 (2012) 241 - 248

www.elsevier.com/locate/procedia

2012 AASRI Conference on Power and Energy Systems

Simple Mathematical Model of a Thermal Storage with PCM

Fortunato B.*9 Camporeale S.M., Torresi M.? Albano M.

Aim of this work is to characterize the thermodynamics of a thermal storage system based on the latent heat of a paraffmic Phase Change Material (PCM). The heat exchange between the heat transfer fluid and the PCM and its phase change are investigated. Under simplifying assumptions, it is shown that the governing equations are the three energy conservation equations written for the heat transfer fluid, the wall and the PCM. The PCM energy conservation equation is written in terms of enthalpy. All three parabolic differential equations are numerically solved according to a finite difference approach implemented in MATLAB language. The simulations are conducted in order to derive the temperature profiles inside the PCM and to estimate the time needed to complete the melting of the PCM. Different simulation campaigns are conducted varying the inlet temperature of the heat transfer fluid. The obtained results can be used as guidelines for more careful studies of the thermal storage system and for subsequent analysis and optimization processes.

©2012PublishedbyElsevierB.V. Selectionand/or peerreviewunder responsibilityofAmericanAppliedScience Research Institute Keywords: CFD; PCM; thermal storage system.

1. Introduction and objectives

Fossil fuel depletion and air pollution have pushed the world to face up the environmental sustainability challenge. Environmental sustainability is expressed through the use and optimization of renewable energy. In the field of energy saving, heat storage plays a very important role.

* Corresponding author. Tel.: +39-080-596-3223; fax: +39-080-596-3411. E-mail address: fortunato@poliba.it.

Politécnico di Bari, Via Re David, 200, Bari 70125, Italy

Abstract

2212-6716 © 2012 Published by Elsevier B.V. Selection and/or peer review under responsibility of American Applied Science Research Institute doi:10.1016/j.aasri.2012.09.041

In the field of power supply, energy demand and energy availability often do not coincide in time; therefore, the use of thermal storage system could be desirable. There are three main ways for storing thermal energy, i.e., in the form of: sensible, latent or thermo-chemical energy. The heat storage in the form of latent heat can be effectively achieved by means of PCMs (Phase Change Materials), i.e. materials characterized by high latent heat of fusion, which through solidification and melting are able to store and release heat.

The thermal storage, in the form of latent heat, has the important advantage of accumulating a large amount of heat in a very small volume, hence presenting a high thermal storage density. The work of Zalba et al. [1] presents an exhaustive review with regard to the classification of PCM. PCMs are generally divided into two main categories: organic (paraffinic and non-paraffinic) and inorganic (mainly hydrated salts). Organic materials have many advantages such as: congruent melting; non-corrosive behavior; high latent heat of fusion; broad range of melting temperatures (depending on the number of carbon atoms in the molecule); minimum tendency to sub-cooling; chemical stability; scarce or null toxicity. A common drawback to all organic materials is the low thermal conductivity, which determines the formation of very high thermal gradients during the transient heat transfer. The main advantage of inorganic PCMs is the high storage capacity per unit volume (about twice that of the organic compounds). Their main drawback is that they melt in an incongruent way. The thermal storage system, under investigation, is composed of a shell and tube heat exchanger. The heat exchange is performed between the Heat Transfer Fluid, which flows inside the tubes, and the PCM inside the shell.

In its generality, the problem consists of the study of conjugated heat transfer and phase change. This type of study leads to the formulation of a mathematical model that provides the definition of moving solid-liquid interfaces and therefore varying boundary conditions. The analytical solution of this kind of problems, known as the Stefan problem, is possible only for very few cases with specific geometric conditions and simple boundary conditions. Shabgard et al. [2] have formulated a mathematical model that provides the variable solid-liquid interface analysis and the natural convection in the liquid phase of PCM. Since the analytical solution of such a complex problem is very difficult, it is necessary to perform a numerical type solution to reach a sufficiently accurate solution of the heat transfer and phase change. The most commonly used methods for the analysis of the phase change is the enthalpy method ([3]-[8]) and the method of equivalent heat capacity based on the temperature. The enthalpy method consists in the formulation of the energy conservation equation in term of enthalpy and considers the solid and liquid phases as a single domain without the explicit definition of the solid-liquid interface. Trp et al. [9] have analyzed the phenomenon of heat transfer and phase change using CFD. Furthermore, Trp [10] presented a detailed mathematical and numerical analysis of the convective heat transfer between the operating fluid and the PCM. In this work the mass, momentum and energy conservation equations are discretized by means of a finite element scheme and solved simultaneously using the SIMPLER algorithm.

In this paper the unsteady heat transfer between the operating fluid and the PCM is numerically investigated using a finite difference scheme. The operating fluid is air and the PCM is organic paraffin. The convective heat transfer coefficient is computed by means of empirical relationships as functions of the Reynolds number, the Prandtl number and the geometry of the heat exchanger. Considering suitable and realistic simplifying assumptions, the governing equations of the problem are the energy conservation equations written for the HTF (Heat Transfer Fluid), the wall and the PCM. The first two equations are discretized using explicit schemes, whilst the last is discretized using a semi-implicit scheme. The phase change of the PCM is analyzed by writing the energy conservation equation in the form of enthalpy of the PCM.

The objective of this paper is to analyze the unsteady heat transfer between air and a PCM in order to track the guidelines for the definition of the thermodynamic performance and design optimization of the storage system.

2. Design of the thermal storage system

The thermal storage system under investigation consists of a "shell and tube" heat exchanger (Fig. 1-a). The tubes, which the Heat Transfer Fluid (HTF) flows inside, are arranged in a regular pattern. For simplicity, in the PCM, fictitious coaxial circular cylinders (tangent one to another) around every tube can be considered limiting the analysis to one single tube (Fig. 1-b). In the default configuration, the HTF (in this case air) enters inside each tube at Tin = 600K and pin = 1 bar with a flow velocity, v = 0.5mis. Each tube has an inner radius, Ri = 0.035m and a thickness, s = 1 mm. The radius of every singular elementary module is R0 = 0.061 m [11]. The number of elementary modules is « = 134.

Fig. 1. (a) Schematic of the thermal energy storage system and (b) elementary module [11]

The computational domain, limited to one single elementary module, can be assimilated to a dual-flow heat exchanger. The heat storage medium is an organic paraffin, whose properties are listed in Table 1.

Table 1. Thermo physical properties of the organic paraffin

Thermo physical properties Units magnitudes

Melting/solidification Temperature [K\ 300.7 Latent heat of fusion [kJ/kg] 206 Thermal conductivity solid/liquid [W/mK\ 0.18/0.19 Specific heatsolid/liquid [kJ/kgK] 1.8/2.4

Density_[kg/m3]_789/750

3. Mathematical and numerical model

The problem of phase change of the paraffin conjugated to conduction in the tube wall and to the forced

convection of the HTF is unsteady. The simplifying assumptions, which the mathematical model relies on, for the description of the physics of the problem, are as follows: the PCM is homogeneous and isotropic; the thermo physical properties of the PCM are constant with temperature; natural convection in the liquid phase of the PCM is neglected; azimuthal temperature are everywhere negligible (2D solution).

The effect of natural convection in PCM liquid phase is widely neglected in modeling PCM thermal system (as done by Trp [11]). The main reason is that, considering only the conduction mode of heat transfer in the PCM, the analysis of thermal systems is highly simplified. Furthermore, with Rayleigh numbers lower than 1' 106, this hypothesis leads to an acceptable margin of error [12].

3.1. Equations

Considering these assumptions, the governing equations can be written in a cylindrical reference system. They are: the conservation of mass and energy and the momentum balance equations for the HTF; the energy conservation equations for the walls and the PCM. Regarding the PCM, temperature is a function of the enthalpy (Fig. 2).

Fig. 2. Isothermal phase transition

In order to determine the PCM state (solid, liquid or undergoing phase change), the liquid fraction, x = milmtot, is introduced, where is the total mass of the PCM, whereas ms and mt are the masses of the PCM in solid and liquid states, respectively. The liquid fraction can be also written as / = (H - cs Tm)/q, hence:

' ^0 H < csTm

0 <X<1 csTm < H <csTm + q. (1)

. X>1 H >c,Tm + q

The system of equations, described above, requires the solution of the flow field, the heat transfer and the phase change simultaneously, using finite element technique [13]. Additionally, the HTF flow is assumed inviscid, incompressible and ID [14]. These three hypotheses implies that it is possible to decouple the equations of motion from the equations of energy and to ignore the continuity and momentum equations of the HTF. The main heat exchange mechanism within the HTF is convection. Hence it is possible to model the heat transfer within the HTF through the definition of a convective heat transfer coefficient h. Known the Reynolds number, Re = 1150, the Prandtl number, Pr = 0.727, and the characteristic size, Du of the heat exchanger, it is possible to calculate the Nusselt number, Nu = 1.86 (Re Pr DJL)1/3 [15], [16]. Knowing the Nusselt number, the coefficient of convective heat exchange can be determined as h = Nu A/Dj.

Considering the tube cross section, Ahtf, the energy conservation equation for the HTF becomes:

n ^ dlk. --Ch'ftilh'f ar*'/ + hnD• IT -T Ph'f Mf at - Ahlf dx + Ahlf T -

htf = -htf'"-h,f + nji T _ T ) ^2)

The thermal resistance of the wall is assumed to be negligible, so that assumption arises that the equation of energy conservation of the wall is:

P. - nT _t; ) + T )

where Aw is the cross section area of the wall, and and kpcm are the thermal conductivity of the wall and the PCM, respectively, and Awi and Avcml the side surfaces of the wall and the PCM, respectively.

The conduction heat transfer within the PCM in the axial direction is neglected. Based on this assumption, the energy conservation equation in the first cell of PCM adjacent to the wall becomes:

C u A1 k A' ^

W W + Pcm Pcm

\ pcm w

On the other hand the energy conservation equation for the generic cell of the PCM becomes:

dHpcm _ kpcm i d ( dHp

r dr^ dr

3.2. Initial and boundary conditions

The physics of the problem is governed by the differential equations 2-5. They are parabolic differential equations that need, in order to be solved, initial (Table 2) and boundary (Table 3) conditions.

Table 2. Initial conditions, t = 0s

Radial coordinate Axial coordinate Initial condition

0<R<R¡ Om<x<lm Tuf=Tln Ri<R<Re Om<x<lm Tw = (Thtf+ Tpcm)/2 Re<R<R0_Om<x<lm_Tpcm = 293K

Table 3. Boundary conditions, t > Os

Radial coordinate Axial coordinate Boundary condition

0<R<R¡ x = Om Thtf=

R¡<R<Re x = Om dTJdx = 0

Re<R<R0 x = Om dTpcJdx = 0

R¡<R<Re x = lm dTJdx = 0

Re<R<R0 x = lm dTpcJdx = 0

R=R¡ Om<x< lm AhtjdThtJdr = A„ dTJdr

R=Re Om<x< lm dTJdr = Xpcm dTpcJdr

R=Ro Om<x< lm dTpcJdr = 0

3.3. Numerical Solution

The unsteady 2D axisymmetric problem described by the differential equations 2-5 is numerically solved using a finite difference scheme. Being i the generic axial node, the spatial derivative is discretized by means of a first order forward finite difference in the first point and a first order backward finite difference for all the other points. Eq.2 and Eq.3 are discretized in time using the forward Euler first-order accurate scheme. Being m the time step, Eq.6 for the HTF and Eq.7 for the wall are obtained:

Thtf (ra +1, i)-Thtf {mj) 1

Phtf chtf

d Thtf (m,i) hnD¡

(t (ra, i)-Thtf (m, i $

T-(+^M - -P- ('■))+i ¥

A' Aw Pw C„, A,,

kpcmA'pcm 1 ((m,i)—Tw (m,i) )

sw pw cw

In the energy conservation equation for the PCM, the second order spatial derivative, indicating the thermal diffusion in the radial direction (/' index), is discretized according to a centered finite difference scheme with second order accuracy. In the time discretization, the semi-implicit, first order accurate, Crank Nicholson scheme is adopted:

Hpcm (m +1)-Hpcm (m, i) l

1 a ( a'H,

dr 1 dr

pcm pcm

i a ( ah„,

dr 1 dr

The computational domain is discretized using 21 nodes in the axial direction and only 5 nodes in the radial direction. This is a rather coarse grid but the main idea is to develop a fast code capable to rapidly provide an estimation of the time scale of the melting process. The grid can be easily refined increasing, however the computational costs. Once Ax and Ar are been defined, the time step At is computed in such a way to respect the stability of the numerical calculation. In particular At = 0.1 minfAx, Ar)/(2vin).

The numerical solution is based on the computation of the values of the fluid and wall temperatures for every z'-th node at the new time as a function of the values measured at the previous time according to the explicit Euler scheme. These temperature values are set as boundary conditions for the solution of the PCM energy conservation equation. Observing the spatial discretization of the PCM energy conservation equation, it is observed that the value of the enthalpy in they'-th node depends on the enthalpy value assumed in nodesy-1 andy+1. In this way, for every z'-th axial node there is a tridiagonal block system.

The numerical solution of this system is based on the Thomas algorithm, which consists in transforming the coefficient matrix into a tridiagonal upper triangular matrix and applying a process of back substitution [17]. The simulations are carried out until the difference between the temperature of the HTF, the wall and the PCM is at least equal to 10"5.

4. Results

The solution of the above equations allows us to determine the profiles of temperature for each node and at each time in the PCM.

The PCM absorbs the heat transferred from the fluid through the wall. From the radial temperature distribution it is possible to observe how, inside the PCM, a thermal gradient is generated, whose slope decreases as time increases. The high thermal gradients depend on the low thermal conductivity of the PCM and its large thermal inertia (high specific heat).

Furthermore, the increase in temperature over time becomes slower with increasing radial abscissa (Fig. 3-a). In fact, when the heat is transferred from the HTF to the first cell of the PCM, a portion of it is accumulated whereas the remaining portion is transmitted to the adjacent cells far away, hence receiving only a small amount of the heat flux. Moreover, as the heat is transferred to the PCM, the enthalpy of the cell at immediate contact with the wall reaches its incipient melting condition and, until the entire latent heat of fusion is not transferred, the temperature remains constant. As a result, the thermal gradients between the various cells are lacking and there is a further slowing of the increase in temperature.

Fig. 3-b shows that the time needed to complete the melting process increases moving radially outward. Furthermore the melting rate slowdowns also in the axial direction.

The same considerations can be performed for different values of the HTF inlet temperature. The temperature distributions in the axial position with i = 11, Fig. 3-c and Fig. 3-d, show how lowering the inlet temperature from 600K to 475K (Fig. 3-c and Fig. 3-d) and 350K respectively (Fig. 3-e and Fig. 3-f), the trend in the radial temperature distribution is preserved but the time needed to reach a regime configuration increases.

(a) T ~ 600K v - ft. 5 m/s

Ub) T-600K v - 0. .1 in s I

\(c) T = 475K

■■ ft 5 in s

\(d) T = 475K v = 0.5 m/s

-.lOifiiii —: — i li

-it-1+1-2 h —o—: it

'•it

(e) T ~ 350K v - ft. 5 m/s

60« = ?«

f2S FOII 4"5 Jrll

425 400

375 .<?» 325 300 275

15 30 45 (¡0

<10 lOi'^liO

\(f) T = 3>0K v - 0.5

Fig. 3 Temperature distributions for different radial positions (left) and time histories of the PCM temperature (right) at section / = 11

5. Conclusions

In this work, the study of heat exchange between a Heat Transfer Fluid (HTF) and a Phase Change Material (PCM) has been analyzed. The solution of the temperature field inside the PCM has been obtained by

numerically solving the three differential equations governing the physical problem. The first two differential equations related to the energy conservation of the Heat Transfer Fluid and the wall have been discretized in time with an explicit scheme, while the PCM energy conservation equation has been discretized with a semi-implicit Crank-Nicholson scheme, in order to have bigger time steps, due to the better numerical stability behavior. The numerical model has proved to be stable and able to give the problem solution very fast (about lA hour of simulation with a standard PC for 2 physical hours of the transient problem).

The numerical simulations have demonstrated that the diffusion of heat in the PCM is very slow: the temperature gradients within the PCM are very low and its melting occurs very slowly along the radial direction.

The main advantage of the simple theoretical model proposed in this work is the great computational simplicity with consequent high calculation speed, obtaining results that accurately show the unsteady behavior of the accumulation system here investigated.

References

[1] B. Zalba, J.M. Marin, L.F. Cabeza, H. Mehling, "Review on thermal energy storage with phase change materials, heat transfer analysis and applications", Applied Thermal Engineering, 2003

[2] H. Shabgard, T.L. Bergman, N. Sharifi, A. Faghiri, "High Temperature latent heat thermal energy storage using heat pipes", International Journal of Heat and Mass Transfer, 2010

[3] D.J. Morrison, and S.I. Abdel Khalik, "Effects of Phase-change Energy storage on the performance of air -Based and liquid-Based Solar Heating Systems", Solar Energy, 1978

[4] C. Bellecci, M. Conti, "Transient behaviour analysis of a latent heat thermal storage module" Int. J. Heat Mass Trans. 36, 3851-3857, 1993

[5] M. Lacroix, "Numerical simulation of a shell-and-tube latent heat thermal energy storage unit", Sol. Energy 50, 357-367, 1993

[6] C. Bellecci, M. Conti, "Phase change thermal storage: transient behaviour analysis of a solar receiver/storage module using the enthalpy method", Int. J. Heat Mass Transfer 36, 2157-2163, 1993

[7] C. Bellecci, M. Conti, "Latent heat thermal storage for solar dynamic power generation", Sol. Energy 51, 169-173, 1993

[8] M. Lacroix, "Study of the heat transfer behaviour of a latent heat thermal energy storage unit with a finned tube", Int. J. Heat Mass Trans. 36, 2083-2092, 1993.

[9] A. Trp, B. Frankovic, K. Lenic, "An analysis of phase change heat transfer in a solar thermal energy store", in: Proceedings of ISES 1999 Solar World Congress, vol. Ill, Jerusalem, Israel, 1999

[10] A. Trp, "A thermodynamic analysis of heat storage in a latent heat storage unit", PhD thesis, Faculty of Engineering University of Rijeka, Croatia (in Croatian), 2002

[11] A. Trp, "An experimental and numerical investigation of heat transfer during technical grade paraffin melting and solidification in a shell-and-tube latent thermal energy storage unit", Solar Energy 79 (2005) 648-660.

[12] E.M. Alawadhi, "Phase change process with free convection in a circular enclosure: numerical simulations", Computers & Fluids, 33 (2004) 1335-1348

[13] S.V. Patankar, "Numerical Heat Transfer and Fluid Flow" Hemisphere Publishing Corporation, New York, 1980

[14] J. Banaszek, R. Domanaski, M. Rebow, F. El-Sagier, "Numerical analysis of the paraffin wax-air spiral thermal energy storage unit", Applied Thermal Engineering 20 (2000) 323-354

[15] Y Cengel, "Termodinámica e trasmissione del calore"

[16] H.Schlichting, "Boundary- layer Theory"

[17] C.A.J. Fletcher, "Computational Techniques for Fluid Dynamics" Vol 1 e Vol 2.