A Study of Thermal Energy

Y. Cao Storage Systems With Conjugate

A. Faghri
Turbulent Forced Convection
Department of Mechanical and Materials
A phase-change energy storage module 'with a turbulent transport fluid is studied.
Engineering, The forced convection due to the turbulent transport fluid is solved with the k-e
Wright State University, model and coupled with the phase-change solution in the phase-change material
Dayton, OH 45435 (PCM). The numerical method is first compared with previous investigations, then
conjugate computations for the energy storage module with different cases are
conducted. A correlation for the energy storage capacity is derived and its empirical
constants are determined by the numerical data. The correlation provides a reference
for the module design and parameter optimization.

Introduction
The study of phase-change thermal energy storage systems row and Hsu (1981) and by Shamsundar (1982) can be used,
is important due to applications in space-based power plants in which the solution for the PCM is coupled by an energy
and solar energy utilization (Viskanta, 1983). Another appli- balance over the transfer fluid by assuming a heat transfer
cation of phase-change energy storage is space-based heat-sink coefficient between the transfer fluid and the PCM. In this
systems. Since phase-change materials (PCMs) have a large technique, the heat transfer coefficient must be calculated using
latent heat, it is an efficient way to absorb the waste heat empirical correlations or other numerical solutions based on
energy in a short time when systems receive heat input and the steady-state constant heat rate or constant surface tem-
release it to space over a long period of time. perature condition. However, for transfer fluids with low
Thermal energy systems usually consist of a number of PCM Prandtl numbers, the heat transfer coefficient at the PCM-
modules with various possible configurations. Figure 1 presents transfer fluid interface changes significantly along the axial
a typical module, which consists of a hollow cylinder of phase- direction, and the boundary condition there is neither steady-
change material with a transfer fluid pumped through the state constant heat rate nor constant surface temperature. Also,
interior tube for the purpose of heat energy exchange. The axial conduction in the transfer fluid cannot be neglected for
study of these PCM modules is essential to understand the low-Pr fluids. In addition, all the results presented in the pre-
characteristics of the PCM energy storage system. The problem vious studies, except those of Shamsundar (1982) for St = 0,
by nature involves time-dependent phase-change heat transfer are only some specific cases and general trends for the module
combined with conjugate forced convection. performance, and no explicit relation for performance param-
Sparrow and Hsu (1981) numerically solved the conjugate eters related to the module design was given. This will surely
phase-change/convection problem for two-dimensional freez- limit the application of the previous studies.
ing on the outside of a coolant-carrying tube with a constant In this paper, the conservation equations for turbulent flow
heat transfer coefficient h for the coolant flow. Shamsundar were solved by applying the k-e model, which was coupled
(1982) presented an analytical closed-form solution for the with the heat conduction in the PCM. The numerical model
same problem by neglecting the heat capacities and axial con- was examined by comparing solutions for different cases with
duction effects in the PCM. The analytical results compared those of Shamsundar (1982) and Sparrow and Hsu (1981).
favorably with the numerical solution of Sparrow and Hsu Correlations for important parameters such as the energy stor-
(1981) for the case with Stefan number St = 0. In both of the age capacity were derived. The constants in the correlation
above investigations, the thickness of the PCM outside the were determined by the numerical data generated from the
circular tube was considered to be infinite, the PCM was in- conjugate numerical computation. To the authors' knowledge,
itially at its freezing temperature, and only the heat conduction this is the first effort to obtain such correlations for the PCM
in the solid region was considered. energy storage system design.
More recently, the PCM module similar to that in Fig. 1
was studied by Solomon et al. (1986), Stovall and Arimilli Mathematical Formulation
(1988), and Yimer and Adami (1989). The studies were focused
on the diffusion-controlled heat transfer in the PCM, with the Conservation Equations for the Transfer Fluid. In energy
heat transfer between the transfer fluid and the PCM calculated
using empirical correlations for the heat transfer coefficient.
Since the temperature field of the transfer fluid continues to ^ PCM ^
change as the phase-change interface progresses, the use of a
steady-state fully developed empirical heat transfer correlation !
and the assumption of a constant fluid temperature may result
in a significant error in the evaluation of the system perform- r ->
.A r
ance. To account for the temperature change in the transfer
-1
R n
X "1 W o
T-
fluid along the axial direction, the technique adopted by Spar- ln „/
^^______^

1
Contributed by the Heat Transfer Division and presented at the ASME Winter
Annual Meeting, Atlanta, Georgia, December 1-6, 1991. Manuscript received
by the Heat Transfer Division March 1991; revision received June 1992. Key-
words: Conjugate Heat Transfer, Phase-Change Phenomena, Thermal Energy
y PCM
V/////////////////////////////////////////,
- L
.J
Storage. Associate Technical Editor: J. H. Kim. Fig. 1 Schematic of PCM energy storage module

Journal of Heat Transfer NOVEMBER 1992, Vol. 114 /1019

Copyright © 1992 by ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use
 storage systems, large amounts of heat must be transferred
between the P C M and the transfer fluid. Under this operating 0»O = V-(pVe + qe) + Se (6)
dt
condition, the flow rate of the transfer fluid must be high,
and the flow is turbulent in most cases. In this study, the where
incompressible Reynolds equations and the standard k-e model Sk=G- pe
were applied to the transport fluid. The turbulent time-aver-
aged continuity, momentum, and energy equations in a two-
dimensional cylindrical coordinate system can be written as: Se^-(ClG-C2pe)
k
13 „ dU
(1)
r dr dx qk=-TkVk, Tk = ^
Ok
d(pV) d 1d , dp
ot dx r dr dr
1 9 dV\ d
q e =-r £ ve, r£ = ^
+ r or
T r + /*eff ., I +S„ (2)
^Tr dx

dt dx r dr dx
2
Mi^Mi
The standard turbulence model constants are (Launder and
+
(dU dV"
y dr dx/

Spalding, 1974):
d_U d_
/>eff / W dx + SU (3)
r dr dr dx C„ = 0.09, C, = 1.44, C2=1.92, o>=1.0, CT£=1.3 (7)

dT d 1d Id dt Boundary Conditions for the Transfer Fluid:

rT„
dt• +dx
—.(UT) +r- dr
— (rVT) =r- dr dr
Inlet plane: x = 0

dx
r e f f f) (4) 0<r<Rr. U=UiH, T=Tin (8)
Outlet plane: x = L
where
r, dU 9V dT
dU\ 13 dV 11 0<r<Rf. — = 0, — = n0, — = 0 (9)
1 + ft rn, dr rkp dx dx dx
*"" r dx dr r dr r dr
For the turbulent flow and heat transfer near the fluid-wall
d ( art 1d dV a
k interface, the wall functions are applied (Launder and Spald-
+ Tdx \-, P
r dr
^•Yx ing, 1974; Djilali et al., 1989).
!*<;{[= P +1*1
Wall shear stress:
1 eff = « + T T -
^ = - In (y+PE) (10)
Pr,
Hi = pCJ^/e is the turbulent eddy viscosity, k is the kinetic where
energy of turbulence, and e is the dissipation rate of turbulence.
The differential equations conserving the transport of k and
UT = ^Jp~, y+p = ^ r , £ = 9 . 0 .
e are

(pk) = V-(pVk + qk) + Sk (5) This expression is used to evaluate the near-wall diffusion terms
dt in the momentum equations.

h = heat transfer coeffi-

c = specific heat, J/(kg-K) cient, W/(m 2 -K)
Cm = specific heat of mushy k = kinetic energy of turbu- Qm = energy storage den-
phase = (cs + ci)/2, 3/ lence, m 2 / s 2 sity = Q,/M, J/kg
(kg-K) K = thermal conductivity, qjw = heat flux at the wall-
cP = specific heat of P C M , W/(m-K) PCM interface, W/m 2
J/(kg-K) L = length of P C M module, R = phase-change interface
C\, C2, C^ = turbulence model con- m radius, m
stants M = total mass of the P C M Ri = inside radius of the
c = coefficient
J/(kg-K)
in E q . (14), module, kg pipe, m
Nu/ = Nusselt n u m b e r R0 = outside radius of the
C, C" = correlation constants = 2R,h/Kf PCM module, m
D = diameter of the circular p = pressure, N/m 2 Rw = radius of the wall-PCM
cylinder = 2 Rh m Pr/ = fluid Prandtl number interface, m
E = constant in near-wall = v/a Re/ = fluid Reynolds number
description of velocity Pr, = turbulent fluid Prandtl = UJD/v
profile number S = term in Eq. (14), J/kg
8i, gs = integrals in Eq. (19) Q = heat energy, J S2 = characteristic area of
G = source term in Eqs. (5) Q, = total energy stored or the annular PCM
and (6) energy storage capacity, = TT(R20-RJ)L/(R0
H = latent heat, J / k g J -Rd, m2

1020 / Vol. 114, NOVEMBER 1992 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 Turbulent Kinetic Energy: The turbulent kinetic energy k
is neither convected nor diffused into the control volume across P = 9.0
Pr,
the wall. A zero normal velocity assures the former, while • » &

setting the boundary condition dk/dy I „ = 0 guarantees the lat- and qfw should be coupled with the PCM solution.
ter. Also,
Conservation Equations for the PCM. The phase-change
In (£>;) model applied is the efficient temperature-based fixed-grid for-
(11)
«-yP mulation that was developed by Cao and Faghri (1990). This
is used to evaluate the dissipation term - p e in Eq. (5). model has the advantage of eliminating the time step and grid
size limitations that are normally encountered in other fixed-
grid methods.
Dissipation Rate: At a node p adjacent to the wall in the
d ( p C 7 ) _ i a (KldT) L(K^L d(pS)
transport fluid, e is evaluated under the local equilibrium con- (14)
dition as dt r dr \ dr 1 dx \ dx dt
where

(T<Tm-dT)

Cm +
H
cm-- 2df ^T>»-ST^T^Tm + 5T)

c, (T>Tm + 8T)
c&T- Tm) (T<T,„-bT)
H H
5(7) = c„, + ^.J T„, + cm5T+ - {Tm-8T<T<Tn + 57)

-CiTm + c£T+H (T>Tm + oT)

Ks (T<Tm-8T

K(T)-- Ks+(K,-Ks)(T-Tm + bT)/2oT {Tm-bT<T<Tm + bT)

K, (T>T„, + 8T)

Equation (14) is also applicable to heat transfer in the con-

(12) tainer wall with C = cm K = Kw, and S = 0. In addition,
i<yP
the possibility of natural convection occurring in the PCM was
Temperature: Given a heat flux qf„, the interface temper- neglected due to the primary space application concerned. The
ature can be obtained from initial condition for the system and boundary conditions for
the PCM and pipe wall are as follows:
(Tfv-Tp)pcjCxt% 1
= Pr, In (£>;) + /> (13)
Qfw Initial Conditions: t = 0
where Entire domain: 0<x<L, 0<r<Ro

Ok, 0€ == turbulence model con-

St = Stefan number yP = distance from the wall stant
= cp{Tin~Tm)/H for T == Fourier number = apt/
to the adjacent node p,
melting or cp(Tm m S2, or shear stress,
- Tin)/H for freezing a = molecular thermal dif- N/m 2
T = temperature, K fusivity, m 2 /s Subscripts
Tfm = fluid mixed-mean tem- r = diffusion coefficient eff == effective
perature, K 8W = thickness of container / == transfer fluid
T,„ = melting or freezing tem- wall, m fw --= fluid-wall interface
perature, K 25T° = phase-change tempera- i == initial condition
t = time, s ture range or mushy in == inlet
5 == solid PCM, or sensible
uin == inlet velocity, m/s
mean velocity in the
phase range, K
e = dissipation rate of tur- heat energy
um pipe, m/s bulence, m 2 /s 3 / == liquid PCM or latent
Up = velocity at node p, m/s K = von Karman constant heat
U, V = time-averaged turbulent /* = molecular viscosity, m -= mushy phase or mean
velocities, m/s kg/(m-s) value
VP = total PCM volume fit = turbulent viscosity, P --= PCM or grid point near
of the module kg/(m-s) fluid-wall interface in
= ir(R20-RhL,m3 V = kinematic viscosity, the transfer fluid
V = velocity vector m 2 /s t --= turbulent
3 w -= container wall
x, r = coordinate directions P = density, kg/m

Journal of Heat Transfer NOVEMBER 1992, Vol. 114 /1021

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 T=T, (15) - C lfw

Boundary Conditions: t > 0

-+f ^~
Outer Wall: 0<x<L, r = R0 Wal
dT
(16) Fluid-Wall
3r=° T,„ - @ 0 #—' —•©• Interface
At both ends: x = 0 and L, Rj-<r<R0
& © m ® ®- •o—L

^=0 (17)
dx
PCM-wall interface: 0<x<L, R = RW
Transport Fluid
dT dT
K„ =KW Fig. 2 Schematic of grid arrangement around the fluid-wall interface
dr dr

Fluid-wall interface: 0<x<L, r = Rj For the present problem, however, the transport fluid is
T=Tfw (18) turbulent so this approach is not applicable. The reason is that
the temperature gradient in the transfer fluid is very steep at
The transient conjugate analysis of Cao and Faghri (1991)
the fluid-wall interface. In this numerical model, the transfer
for laminar flow showed that the velocity reached the steady
fluid and PCM-wall regions were solved separately and cou-
state very quickly from the initial zero velocity, while the tem-
pled via the wall function. Figure 2 represents the grid ar-
perature counterpart continued to change as the phase-change
rangement around the fluid-wall interface. The two-domain
front moved in the PCM. Therefore, the velocity transient
problem was solved sequentially and iteratively. During the
period can be ignored without losing generality in the numerical
iterations at each time step, the PCM-wall region was first
modeling. Also, for most engineering problems the inlet ve-
solved with the T/w from the previous iteration as a boundary
locity can be assumed to be in the fully developed turbulent
condition. The Qf„ was then estimated by the relation:
state because of the connecting pipeline between the PCM
system and the pump. These assumptions simplified the model dT
and reduced computation time significantly. (21)
The key module parameters used to evaluate the perform-
ance of a thermal energy storage system are the energy storage using the temperatures at the grids in the pipe wall close to
capacity or total energy stored/extracted, Q, (J), and the energy the interface. This q/„ was then applied in the turbulent wall
storage density Q,„ = Q,/M (J/kg), where M i s the total mass function, Eq. (13), to solve for the temperature in the transfer
of the PCM module. fluid. The interface temperature TfW was recalculated from the
transfer fluid solution and used as a boundary condition for
The total energy stored Q, consists of two parts, i.e., the
the PCM-wall region in the next iteration. The whole iteration
latent heat energy and the sensible heat energy:
procedure was repeated until a converged solution was ob-
tained for each time step. The converged results were assumed
Q< = Qi + Qs = PHT \ (R\x, t) - RJ)dx to be reached when the maximum relative change of temper-
ature between the consecutive iterations was less than 5 x
10" 5 . Different grid sizes for the same problem were tested
+ 2irp cJT- Ti)rdrdx = pHirg, + 2irpgs (19) and it was demonstrated that the numerical solutions were
essentially independent of the grid size used. A change from
25 x 29 to 50 x 56, for example, resulted in changes in the
where g, = \ {R.\x, t)-RJ)dx,gs = \ ° [ cp(T-T,) rdrdx, energy storage capacity and the outlet transport fluid temper-
ature of about 2 and 1 percent, respectively. The grid size and
R(x, t) is the radius of the melting interface at time t, and 7} time step are given in the next section for the case presented.
is the initial temperature of the PCM system.
Results and Discussion
Numerical Procedure and Fluid-Wall Interface Treat- The computer codes for both the turbulent model of the
ment transport fluid and the phase-change model of the PCM in a
The problem specified in the last section was solved by ap- cylindrical coordinate system have been checked against two
plying the control-volume finite-difference approach described correlations and a numerical solution in the literature, re-
by Patankar (1980, 1988). The discretization equations were spectively, to validate the numerical accuracy. Figure 3 shows
obtained by applying the conservation laws over a finite size the comparison of the present numerical solution and the ex-
control volume surrounding the grid node and integrating over perimental correlations for the fully developed turbulent forced
the control volume. The SIMPLE algorithm was used to solve convection Nusselt number with Pr = 0.706. Both correlations
the velocity and pressure for the transport fluid. Since the fluid shown in the figure were considered to have a high degree of
is incompressible, they are decoupled from the temperature accuracy in representing a large variety of turbulent experi-
solution. For the temperature solution, Cao and Faghri (1991) mental data (Kays and Crawford, 1980; Bhatti and Shah, 1987).
applied the harmonic mean of the thermal conductivity at the As can be seen, the agreement between the numerical solution
fluid-PCM interface for the laminar flow such as (for a uni- and the correlations is good. Figure 4 compares the present
form grid size) numerical phase-change radius with the one-dimensional nu-
merical result of Lunardini (1981). The geometry of the phase-
dT dT 2KpKf dT change system studied was a circular cylinder surrounded by
-Kf (20)
dr r=R-~Kf+Kpdr an infinite, homogeneous medium. The agreement between the
two numerical solutions is very good. More information con-
and solved the transport fluid and PCM as a single-domain cerning the PCM model verification can be found in the paper
problem. by Cao and Faghri (1990).

1022/Vol. 114, NOVEMBER 1992 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 1400 0.8
1000
- U6.I5 (Shomiiundor.1982)
0.7 = MIU S ^ « _ e n i solution) .,
^_730.6 s (Sr.omstmdor.1982) ,
= 7~30.6 s (Present solution)
0.6 =.1461.1 5 (Shamsundor.1982) ,

= 2922.3 s (ShomsiJndor.1982^
=_2<m.~Sj ( P r e s e t Solution) _
= 4_38314 s_(Sh_oms_Ufid_ar,J982i
=_4_383.4 5 (Prssenl solution) '
<f 0.4- = 5844.6 s (5homstindar,1982)
= 5844.6 s (Prosunl saluiion)

0.3
Re, = 4.77*10'
O Present Numerical solution 0 . 2 -I Pr, - 0.706
A Correlation by Sieicher and Rouse (1975) St = 0.0
< T h -T,)/(T - T m ) = 1.0
Correlation by Gnielinski (1976)
1 K,7K = 0.
_ _ L/0 = 50
0.0
0
x/R,

w Fig. 6 Comparisons of phase-change fronts at different times for the

10 100 1000 2000
conjugate one-region phase-change problem (St = 0.0)
Re*10"3
Fig. 3 Comparison of numerical solution and the experimental corre-
lations for fully developed turbulent Nusselt number with Pr = 0.706 0.7 -
q r l/R.'=0.0?5(Sporro w and Hsu)
o r t/R.'=0.025'{Presenl solution)
nJ/R i '=b^05(Spofr O w and Hsu)
0.6 a r l ^ , ' = 0.0"5(Pr.sinl solullon)
1
• . W - - * ' ! * ™ •>"" ""•>
o 0 l/R i 2 =OJ(P^senl_sol i ^lioni
0.5 O p t / R ^ q ^ S p o r r o w and Hsu)
nJ/R^q^Pr^.nljolulid^L .
O p t / j ^ O ^ f S p o r r o w ^ a n d Hsu) _
O Present Numerical Solution 0.4 ^ l / » i ' = " 0 , j ( P r . « n l salulian)
4- QrlA?=0.4fSporr0wandHsu)
A Numerical result by Lunordini (198f n l/ft ( ' = 0.4(Presenl solulipn)
0.3

Re, = 4.77>104
0.2 Pr. = 0.706
3 - SI = 1.

0. 1 ' " = 0.899

L/Di ' 50

0.0
150
x/U

Fig. 7 Comparisons of phase-change fronts at different times for the

conjugate one-region problem (St = 1.0)

0.02 100
(at/R,2)*St
Fig. 4 Comparision of phase-change radius of a cylinder surrounded
of the container wall on the heat transfer in the PCM is very
by an infinite PCM medium small. The same conclusion was also drawn by Cao and Faghri
(1991) with laminar transfer fluids. For this reason, the thick-
ness of the container wall is considered to be negligibly small
ouu - in the following analysis for simplicity.
Before proceeding to the analysis of the energy storage mod-
ule in Fig. 1, the numerical method is examined by simulating
600 - ° <5.../R: = 0-0 the one-region conjugate phase-change problem presented by
Sparrow and Hsu (1981). Figure 6 compares the present so-
lutions for the thickness of the frozen layer with the closed-
400 - ^ ^ " form solution by Shamsundar (1982) using the Newton-Raph-
/ ^ son/Secant method. The agreement between the two solutions
Re, = 4.77*t04
/</ Pr, = 0.706 is very good except at the very early operation time in which
200 -I = 1.0 the transient behavior of the transfer fluid is very important,
L/D = 50
which was not taken into account in the closed-form solution
of Shamsundar (1982). In the analysis of Shamsundar (1982),
0 - 1 t the sensible heat and axial conduction in the PCM were ne-
0.0 0.1 0.2 0.3 0.4 0.: glected, and the initial system temperature was assumed to be
VA2 at the freezing temperature. The transfer fluid temperature
Fig. 5 The effect of the pipe wall on the energy storage capacity was obtained by making an energy balance over the PCM and
the transfer fluid with an assumed heat transfer coefficient.
The solution is a function of Biot number, hR/K, and Stanton
The effect of the pipe wall on the performance of the energy number, h/pUmCj. The case in which the comparison was made
storage module was examined by running cases with a wall of •corresponds to hRj/K = 5 . 0 and h/pUmCf = 0.003. In the
zero thickness, bw/R; = 0.0, and with a fairly thick wall, b„/ present computation, the operational parameters such as Re/
R, = 0.08. The total energies extracted from the PCM for the and K/Kp were determined based on these Biot and Stanton
two cases are plotted in Fig. 5. In the early period, the total numbers, in which the heat transfer coefficient h was obtained
energy extracted from the PCM with a thick wall was slightly via the correlation by Gnielinski (1976) for a given Prandtl
less than that with a wall of zero thickness. Over time, however, number. It can be seen from the comparison that little error
the total energies extracted for the two cases are almost the was introduced by using a constant h for the gaseous transfer
same. Since the pipe wall has a much higher thermal conduc- fluid with relatively high Prandtl number except in the region
tivity than that of the PCM, which provides a significant heat very close to the entrance. Figure 7 compares the present so-
path between the transfer fluid and the PCM, the influence lution and the numerical solution by Sparrow and Hsu (1981)

Journal of Heat Transfer NOVEMBER 1992, Vol. 114 /1023

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 0 . 8 "I ;
i r l/R.'=O.Q25{Spdrrow and Hsu),

mwzzm

x/Ri

Fig. 8 Comparisons of phase-change fronts at different times for the Fig. 10 Mixed-mean fluid temperature along axial direction at different
conjugate one-region problem with lower Prandtl number time periods for the melting problem

1 -r* - ^ ,
problem in which both the solid and liquid regions need to be
•O O T = 1.4233-10-'
considered.
The numerical computation for the thermal energy storage
module shown in Fig. 1 was first made with a melting problem.
The phase-change material studied is LiH, which is considered
to be one of the most promising candidates for the space-based
energy storage system. The system was initially at a temperature
Tj lower than the melting temperature Tm = 962 K of the
PCM. The grid size used in the calculation was 51 (axial) x
[20 (radial fluid) + 50 (radial PCM)] and a time step At = 8
min. The phase-change temperature range was taken to be 87"
= 0.5 to simulate phase change at a single temperature. It has
0 70 140 210 280
been shown by Cao and Faghri (1990) that a moderately small
x/D 57" is enough to simulate the phase-change problem occurring
Fig. 9 Melting front locations at different time periods
at a single temperature.
Figure 9 shows melting front locations along the axial di-
rection for different time periods. The melting interface moves
for St = 1, with the other parameters being the same. The in the radial direction with time. However, the interface moves
comparison was also made with the gaseous transfer fluid and faster upstream due to the temperature drop of the transport
the same trend was observed. Figure 8 compares the present fluid along the axial direction. Although the melting interface
solution and the numerical results by Sparrow and Hsu (1981) curves are sloped along the axial direction, the curves are still
for the same case but with a relatively low Prandtl number, relatively flat owing to the high flow rate and the strong tur-
Pr = 0.05. For the present numerical result, the operational bulent heat transfer between the transport fluid and the PCM.
parameters Re/ and K/Kp were calculated based on Biot num- Figure 10 shows the mixed-mean transport fluid temperature
ber hRj/K = 5, Stanton number h/pUmcf = 0.003, and Pr = along the axial direction at different times, which is defined
0.05 using the correlation by Notter and Sleicher (Bhatti and as
Shah, 1987) for the constant heat rate. The difference between
T
the two solutions is significant. For transfer fluids with low M = -^7T \ ' UTrdr. (22)
Prandtl numbers, the heat transfer coefficient changes signif- KjUm J0
icantly in the axial direction, and the axial conduction and The fluid temperature is initially at the starting system tem-
developing flow effects are dominant along the pipe length. perature T„ which is lower than Tm. As operation began, the
The present numerical results show that the fully developed outlet temperature rose very quickly. After the initial sharp
flow was not reached even at the outlet of the pipe flow, x/ increase, the outlet temperature maintained a continuous but
Rj = 100. As a result, a significant error may be introduced very slow increasing trend for the rest of the operational period.
by employing a steady-state fully developed heat transfer coef- Figure 11 shows the radial temperature distribution at an axial
ficient and neglecting the axial conduction in the transfer fluid. location of x/D = 100 for different time periods. Two different
In this case, the more comprehensive method adopted by the regions (the PCM and the transport fluid) are also indicated
present study is needed to obtain accurate solutions. Another in the figure, noting that the pipe wall between the transport
aspect that needs to be mentioned is that the technique adopted fluid and the PCM has been assumed to be negligibly small.
by Sparrow and Hsu (1981) and by Shamsundar (1982) deals The phase-change front at different times is the intersection
with the phase change on the outside of a coolant-carrying of (T- Tm)/(Tir, - Tm) = 0 and the corresponding temperature
tube surrounded by an infinite PCM with the initial temper- curves. The largest radial temperature drop in the transport
ature being the freezing temperature. In this case, only the fluid occurs at the fluid-PCM interface due to the turbulent
heat conduction in the solid region needs to be considered and nature of the transport fluid.
the sensible heat may only play a minor role when the thickness The numerical computation for the energy storage module
of the frozen layer is relatively small ((i? — i?,)/i?;<0.7). For was then made with a freezing problem. The phase-change
the energy storage module presented in Fig. 1, however, the material considered is octadecane with an initial temperature
initial temperature may be much lower or higher than the T\ higher than its freezing temperature Tm = 301.2 K. Figures
phase-change temperature, and in many cases, the sensible heat 12-14 show the freezing fronts along the axial direction, mixed-
can constitute more than 50 percent of the total energy storage mean transport fluid temperature along the axial direction,
capacity. Obviously, the problem by nature is a two-region and the radial temperature distribution at x/D = 25 for dif-

1024/Vol. 114, NOVEMBER 1992 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 ?tr=~::"^^rr:rrzr*
V T.5.0.255,J0;'
A r- 1.423*10-*
A T = 4.0B*1Q-4
+ T = 4.033'1Q-J
X T = 10.2*10-' -t- T - 13.*1Q-*
X T = 2~3".7*1Q-<
0.5 O T - 36.2*10-4

Re, = 1.1*105 Re, = S.'IO4

Pr, = 0.66 Trv=TJ706~
St = 0.616 Si = 0.425
(T i „-T i )/(T i n -T m )
= 2.93 = 1.3
c,/c, = 0.736
\ \ c /c. = 0.864
-0.5 • K./K„ = 0.106
K ( A p = 0.114 K V = 0.4
K./K;= 0-75 (R -R,)/L = 0.014
(R -R,)/L = 0.005 L/D = 50
L/D = 200 -PCM-
i Fluid — Fluid -
-PCM_-
— I — 0.4 O.i 1.2 1.6 2.4
0.5 1.5 2.5 r/R,
r/R, Fig. 14 Radial temperature distribution at x/D = 25 and different time
Fig. 11 Radial temperature distribution at x/D = 1 0 0 for different time periods for the freezing problem
periods for the melting problem

O T J 1.02*10-
A T = 4.34*10-'
+ T = 13.*1Q-*
X T = 2~3~7*1Q-'
O T = 36.2*10"'

I 0.6
Re. = 5.«104
Pr, = 0.706
SI = 0.425
0.4- (T,„-T,)/(T,n-Tm)
= 1.3
c,/c. = 0.864
K,/KD = 0.106
0.2 K,/K'= 0.4 O Q,. («n-R,)A =
(R„-R,)A = 0-014
L/D = 50. : 0.736 + 9u ' r.RjZkf.
= 0.114 X 0^
K/K = 0.75 O 0,. (>VR,)/L =
10 20 30 40 50 60 70 „ , - T , ) / ( T , „ - T j L/D i 200
0 "2.93 '
x/D
400 800
Fig. 12 Freezing fronts along axial direction at different time periods
t (min.)
Fig. 15 Energy storage capacity and density with different geometry
parameters for the melting problem

Re, = 5.*104
Pr, = 0.706 the storage capacity and the density are time dependent. At
St = 0.425
Or,„-T,)/(r,„-Tj very early times, the energy storage capacities for the three
= 1.3 PCM modules with different (R0-Rj)/L are almost the same.
c /c. = 0.864
K * = 0.106
K X = 0.4
For longer times, the storage capacity of the module with
(R „-R:)/L = 0.014 smaller PCM thickness R0 - R, is lower than those with larger
o T = 0.255*10 L/D = 50
T = 2.04*10-"
thicknesses. At about 400 min, the PCM encapsulated in the
T == 6.S9*10"
4 module with (R0-R^/L - 3 x 10~ 3 was almost completely
T -= 28.3*10-*
4
melted, while the heat storage capacities of the other two mod-
T -= 42.8*10"
0.6 ules continue to increase since the phase change still takes place
in the PCM. For the energy storage density, the module with
small R0 - R, has a higher value in most cases due to the smaller
0.5
20 30 40 70 total PCM in the module. For the module design, optimization
x/D is required, and sometimes a tradeoff is needed among different
Fig. 13 Mixed-mean fluid temperature along axial direction at different
module parameters. It should be pointed out that the results
time periods for the freezing problem shown in Fig. 15 only involve variation of a single parameter
(R0~Ri)/L. In the real design process, many parameters need
to be considered simultaneously. It is important to obtain some
ferent time periods. The general trends of the three graphs are explicit expressions of Q, and Qm of the purpose of design and
similar to those of the melting problem presented earlier. It optimization. Since any purely analytical expressions are not
should be pointed out that in Figs. 13 and 14, the temperature available, an attempt to obtain some semi-empirical relations
distributions are presented in dimensionless form. While the was made in this paper.
trends of the dimensionless temperature distributions are sim- In order to evaluate the integrals g/ and gs in Eq. (19) and
ilar for freezing and melting, the trends are opposite for the to obtain a correlation frame for Q,, consider a simplified
dimensional temperature profiles. This is due to the fact that quasi-stationary model as described below. A hollow PCM
the phase-change temperature is higher than the inlet temper- cylinder of inner radius Rt is initially solid at its melting tem-
ature for freezing and lower for melting. perature. At the surface r = Rh heat is input via convective
Figure 15 shows the energy storage capacity and density at heat transfer with a film coefficient h from a transfer fluid at
different times for the melting problem described earlier with ambient temperature 7}„. In this case, a uniform phase-change
different values of the geometry parameter (R0-Ri)/L. Both front R(t) exists, separating liquid PCM (r<R(f)) from solid

Journal of Heat Transfer NOVEMBER 1992, Vol. 114 /1025

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 lO^i Stefan number is St = cp (Tin - T„,)/H, and the Fourier number
for the P C M is r = tap/S2.
With the expression for R2 - R2 above and a constant specific
heat, we can evaluate the integrals in Q,:

g,= \ (R2-R2)dx=(R2-R2)L

I5»
if
- / *
® Present numerical data
Present correlation
&= f ° t
jR. J0
cp(T~ Tm)rdrdx=Cp ^

In the above integration, we have assumed that the liquid

I
z
d L ^-Rj)/2

0 . 1 -./ • •
P C M has an average temperature of {Tin+ Tm)/2. Therefore,
Q, = pHirgi + 2Trpgs

0.02
0.4 100
= pirL(R2 -R2) H+-*(Tin-T„)
^''"'ft-^'
Fig. 16 Comparison of numerical data and the proposed correlation
= prLHS2 -£ Nu/rSt ( 1 + f St
(r>R(t)). In the liquid region there is a radially symmetric,
time-dependent temperature function T(r, t), satisfying the

* § Nu^stU
L(R20-R2)T:
heat equation = pHL st
Rn — Ri
dT_ad_ / dT
(23)
dt~ r dr \ dr
with boundary conditions
T{R{t),t)=Tm (24) pHL(R20-R2)ir

dT R0-R>\-l K>
MdR Nu/rSt 1+-St (29)
-Kp — (R{t),t) (25) K,

At the inner radius, Rh The dimensionaless energy storage capacity Q,/

pHL{R20 - R2) -K = Q,/MHis the ratio of the total energy stored
d_
Kp-JT^m^hiT^D-TJ (26) at time t to the total latent energy avialable for the P C M
dr module.
Similarly, the dimensionless energy storage density can be
An analytical approximation to T(r, t) and R(t) can be o b -
written as
tained based o n the quasi-stationary approximation a p p r o a c h
(Solomon et al., 1986; L o n d o n and Seban, 1943), by replacing
the heat equation, Eq. (23), with the steady-state equation
Qui.
H
Rn-R
f) N ^ S t ( 1 + , S t (30)

— I r —— 1 = 0 , while retaining the remaining conditions (24)- It should be pointed out that the above expressions have
been obtained based on a very simplified model. The purpose
(26). The solution for the above problem is of the derivation is to obtain an expression frame for the
correlation. The above parameters such as R(f) and h could
Kpt(Tin-Tm) R2-R2 R2 In (R/R) only be considered as average or effective parameters. The
pH 4 2 whole expression must be modified and the final version be
determined by the conjugate numerical data.
= ||0K2-tf/) (27) In the simplified model, the initial temperature 7} of the
P C M was assumed to be the melting temperature Tm, which
Since the magnitude of R: is usually on the order of 10 m m , is not the case in general for engineering applications. There-
the term fore, a correction factor
R2-R}_R2ITL(R/R,) T - T
C 1
m 1 >
4 2
is much smaller t h a n the first term, was included in the expression. With additional exponent m o d -
Kpt(Tj„-Tm) ifications for the rest of the dimensionless n u m b e r s , the cor-
relation can be written as
pH
T — JT K, Ro-Ri
and can usually be neglected. Therefore, C" 1
in / c rl
MH~ K, f f
T -T
R2_RiJ.Riht(Tin-Tm)
T^SKI + O.SSOF (31)
pH
where the constants C", a, b, c, d, e, f, and g are to be
or
determined by the numerical solutions. In the above expres-
R2-R} = S2§LNufStT (28) sion, we have expressed the effective Nusselt number Nuy as
N u / ~ R e / P r } , which is generally acceptable for engineering
where the characteristic area of the annular PCM is S = Vp/ applications.
(R0-Ri), Vp is t h e total P C M volume of t h e module, t h e Numerical computations for the P C M energy storage m o d -
Nusselt number for the transfer fluid is N u / = IRjh/Kf, the ule with different parameters were m a d e extensively, and about

1026/Vol. 114, NOVEMBER 1992 Transactions of the ASME

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 150 numerical data points, most for high-temperature energy (1982) and the numerical solution of Sparrow and Hsu (1981)
storage systems, were taken from the numerical solutions. These provided an efficient way to predict system performance. For
numerical data were used in conjunction with a regression two-region conjugate phase-change problems or PCM modules
analysis to determine the constants in Eq. (31). The expression with lower Prandtl number transfer fluids, the present method
thus obtained is was needed to obtain accurate solutions. The performance of
R0 - Rj the energy storage module in Fig. 1 was simulated systemat-
_2L. = 0.203 Tm-T, K,
ically for different parameters, and correlations for energy
HM~ •* in ' •* n. KP/ storage capacity and energy storage density were obtained using
the conjugate numerical data. These correlations provide guid-
riy.0.089-r"-
0.54r
4
[St(l+0.5St)] u ' U J (32)
C/
ance for module design and optimization.
Figure 16 shows the comparison of the numerical data and
the correlations, Eq. (32). For most of the data points, the
differences between the numerical data and the correlation are References
within 25 percent. Bhatti, M. S., and Shah, R. K., 1987, "Turbulent and Transition Flow Con-
vective Heat Transfer in Ducts," in: Handbook of Single-Phase Convective Heal
The correlation is applicable within Transfer, S. Kakac et al., eds.
Cao, Y., and Faghri, A., 1990, " A Temperature Transforming Model With
Qt
0-20 < T 7 T < 3.0 a Fixed Grid Numerical Methodology for Phase-Change Problems Including
MH Natural Convection," ASME JOURNAL OF HEAT TRANSFER, Vol. 112, pp. 812—
816.
which covers the most useful operating ranges. Cao, Y., and Faghri, A., 1991, "Performance Characteristics of a Thermal
Other parameters covered by the numerical data are Energy Storage Module: A Transient PCM/Forced Convection Conjugate Anal-
ysis," Int. J. Heat Mass Transfer, Vol. 34, No. 1, pp. 93-101.
K, T, Djilali, N., Gartshore, I., and Salcudean, M., 1989, "Calculation of Con-
0.05: <0.72; 1.05 < <3.78;
~K, T —T vective Heat Transfer in Recirculating Turbulent Flow Using Various Near-Wall
1
in x
m Turbulence Models," Numerical Heat Transfer, Part A, Vol. 16, pp. 189-212.
Ro-R Ganic, E. N., Hartnett, J. P., and Rohsenow, W. M., 1985, "Basic Concepts
0.0014< <0.02 0.1<St<1.34 of Heat Transfer," Handbook of Heat Transfer Fundamentals, W. M. Rho-
L senow et al., eds., pp. 15-25.
0.1<Pry<1.0 5 x l 0 4 < R e y < 1.1x10" Gnielinski, V., 1976, "New Equations for Heat and Mass Transfer in Tur-
bulent Pipe and Channel Flow," Int. Chem. Eng., Vol. 16, pp. 359-368.
The expression for energy storage density Q„„ following the Kays, W. M., and Crawford, M. E., 1980, Convective Heat and Mass Transfer,
result for Qt, can be written as 2nd ed., McGraw-Hill, New York.
Launder, B. E., and Spalding, D. B., 1974, "The Numerical Computation
/ s 0.907 / \ 0.029 / /?\-°-327 of Turbulent Flow," Comp. Meth. Appl. Mech. Eng., Vol. 3, pp. 263-289.
1 1 J K K
y»' '- =r,0.203
™-> I in- i\ I ±L\ l o~ i\ Re0-308 London, A. L., and Seban, R. A., 1943, Transactions of the ASME, Vol.
H K, 65, p. 771.
Lunardini, V. J., 1981, "Phase Change Around a Circular Cylinder," ASME
JOURNAL OF HEAT TRANSFER, Vol. 103, pp. 589-600.
TJM [St(l+0.5St)f' U J (33)
•if
Prf Patankar, S. V., 1980, Numerical Heat Transfer and Fluid Flow, McGraw-
which is applicable within the same parameter ranges covered Hill, New York.
Patankar, S. V., 1988, "Elliptic Systems: Finite-Difference Method I , " in:
by the numerical data. For phase-change materials with dif- Handbook of Numerical Heat Transfer, W. J. Minkowycz et al., eds., pp. 215-
ferent values of properties between the liquid and the solid 290.
phases, such as thermal conductivity and specific heat, the Shamsundar, N., 1982, "Formulae for Freezing Outside a Circular Tube With
parameters in Eqs. (32) and (33) can be evaluated based on Axial Variation of Coolant Temperature," Int. J. Heat Mass Transfer, Vol. 25,
No. 10, pp. 1614-1616.
the averaged values over the two phases. Although most of Sleicher, C. A., and Rouse, M. W., 1975, " A Convenient Correlation for
the numerical data used to obtain Eqs. (32) and (33) were taken Heat Transfer to Constant and Variable Property Fluids in Turbulent Pipe
from the melting solutions of high-temperature storage mod- Flow," Int. J. Heat Mass Transfer, Vol. 18, pp. 667-683.
ules, Eqs. (32) and (33) are equally applicable to freezing prob- Solomon, A. D., Morris, M. D., Martin, J., and Olszewski, M., 1986, "The
Development of a Simulation Code for a Latent Heat Thermal Energy Storage
lems of high-temperature storage modules with St = System in a Space Station," Technical Report ORNL-6213.
cp(Tm-Tin)/H. Sparrow, E. M., and Hsu, C. F., 1981, "Analysis of Two-Dimensional Freez-
ing on the Outside of a Coolant-Carrying Tube," Int. J. Heat Mass Transfer,
Vol. 24, Vol. 8, pp. 1345-1357.
Stovall, T. K., and Arimilli, R. V., 1988, "Transient Thermal Analysis of
Three Fast-Changing Latent Heat Storage Configurations for a Space-Based
Conclusions Power System," Proc. 23rd Intersociety Energy Conversion Engineering Con-
A numerical procedure solving the transport fluid and the ference, pp. 171-177.
Viskanta, R., 1983, "Phase-Change Heat Transfer," in: Solar Heat Storage:
PCM as a conjugate problem was presented. It has been shown Latent Heat Material, G. A. Lane, ed., Vol. 1, pp. 153-222.
that for one-region conjugate phase-change problems with the Yimer, B., and Adami, M., 1989, "Parametric Study and Optimization of
initial temperature close to the freezing temperature and gas- Phase Change Thermal Storage Systems," Proc. 1989 ASME National Heat
eous transfer fluids, the closed-form solution of Shamsundar Transfer Conf., HTD-Vol. 109, pp. 1-8.

Journal of Heat Transfer NOVEMBER 1992, Vol. 114 /1027

Downloaded From: http://heattransfer.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use