This paper establishes a multiscale design evaluation framework that integrates performance models for a thermal energy storage (TES) unit and a subsystem heat exchanger (HX). The modeling facilitates the analysis of transient input and extraction processes for the TES device which uses solid–liquid phase change to store thermal energy. We investigate sensible and latent heat transfer through the unit's matrix structure which contains phase change material (PCM) in the interstitial spacing. The heat transfer is driven by a temperature difference between fluid flow passages and the PCM matrix which experiences sensible heat transfer until it reaches the PCM fusion point; then it undergoes melting or solidification in order to receive, or reject, energy. To capture these physics, we establish a dimensionless framework to model heat transfer in the storage device much like effectiveness-number of transfer units (NTU) analysis methods for compact HX. Solution of the nondimensional governing equations is subsequently used to predict the effectiveness of the transient energy input and extraction processes. The TES is examined within the context of a larger subsystem to illustrate how a high efficiency design target can be established for specified operating conditions that correspond to a variety of applications. The general applicability of the model framework is discussed and example performance calculations are presented for the enhancement of a Rankine power plant via asynchronous cooling.

Introduction

This paper summarizes the derivation and solution of the governing equations for the energy input and extraction processes in a thermal storage device interacting with a heat exchanger (HX) that facilitates this input or extraction. Earlier analyses of phase change thermal storage performance have generally modeled specific details of heat transfer in the storage unit structure with constant boundary conditions, thereby neglecting interaction within a subsystem (e.g., see Refs. [13]). More recent explorations of phase change material (PCM) thermal storage have included the evaluation of PCM materials (e.g., see Refs. [4] and [5]) and efforts to model associate transport processes (e.g., see Refs. [68]).

Shamsundar and Srinivasan [1] look at a three-dimensional shell and tube configuration both analytically and numerically (via finite difference) in which the working fluid temperature changes axially as heat is transferred from the PCM. To generate effectiveness charts for the thermal energy storage (TES), the two authors compare the mean frozen fraction to that of the frozen fraction at the inlet for varying sizes, layouts, and Biot numbers. They indicate that the approach could be applied to estimating effectiveness for varying flow rates or inlet temperatures but that the accuracy would be sacrificed.

El-Dessouky and Al-Juwayhel [2] use a second law analysis to characterize the TES by entropy generation numbers. While the authors analyze an entire transient melt/freeze cycle, the PCM remains at the melt temperature throughout. Their case study validates the prediction that more effective TES devices should have higher Reynolds number flows, large heat transfer area, and a greater inlet temperature difference between the working fluid and PCM. Once again, these authors do not investigate beyond constant inlet temperatures.

Ismail and Goncalves [3] explore a two-dimensional model of a tube immersed in PCM. Like our analysis here, they write the energy equation in enthalpy form and couple the change of enthalpy in the PCM with that of the working fluid. Defining an appropriate control volume, the authors employ a finite difference scheme to characterize how the TES melt fraction, number of transfer units (NTU), and effectiveness are dependent on the geometry, Biot and Stefan numbers, and working fluid inlet temperature. As the papers before, these authors apply a constant inlet boundary condition without discussing the heat exchanger that might supply this.

Tay et al. [6] present a one-dimensional simplified NTU-effectiveness analysis for a shell and tube heat exchanger. The authors find that this approach is valid provided that there is a high heat transfer area. They use thermal resistances to describe the heat transfer from the fluid core through the PCM. This simplified analysis is unlike ours in many ways, including those mentioned earlier as well as in the way we determine U.

Most recently, several papers have explored natural convection in three-dimensional phase change cells. Bondareva and Sheremet [8] include both momentum and energy equations in the PCM to account for natural convection as well as conduction. The authors use a constant temperature energy input for the melting process and solve for the three-dimensional velocity and temperature fields within the cell to determine the solid–liquid interface location and, thus, the melt fraction with time. This is different from our paper, despite a similar nondimensionalization technique resulting in some overlap in dimensionless groups. The 2017 paper by Hu et al. [7] also examines natural convection, though these authors employ a three-dimensional lattice Boltzmann model as well as a particle velocity model in order to resolve the temperature field, find the phase interface, and calculate the melt fraction. These authors begin with an enthalpy-based heat equation and keep track of the velocities within the PCM to account for natural convection. Their dimensionless groups are similar to the previous paper and they perform numerical parametric studies by varying one of these.

To summarize, many papers mentioned here have applied NTU-effectiveness type modeling to TES units. These papers have tended to nondimensionalize the governing equations and then employ one of two approaches to evaluate device performance: using working fluid temperature to quantify effectiveness or comparing the local melt fraction to the inlet melt fraction of the device. In this paper, our method of quantifying effectiveness is not significantly different from prior work, but our temporally and spatially varying multiscale modeling approach, depicted in Fig. 1, is. This is novel because it facilitates exploring the relationship between two key dimensionless parameters, t* and Ntu, within the context of a subsystem that models realistic operating conditions for a TES.

While one paper in this literature review discusses the possibility of varying inlet temperature and another examines a simplified cycle, these papers do not analyze a TES device with transient boundary conditions and spatially varying initial conditions for the processes within the cycle. Furthermore, these papers consider the unit cell to the device level scale, but do not examine this within the context of a subsystem that is responsible for the delivery of heat to or from the TES. By incorporating a macroscopic lens and considering the subsystem heat exchanger, we can specify exactly how this technology could fit into and improve a power or refrigeration cycle.

We begin at the TES unit cell and model the local heat transfer between the fluid and the PCM using a mean overall heat transfer coefficient U. This is similar to the analysis of standard heat exchangers, though here, U must be averaged both spatially and temporally. Once the unit cell has been sufficiently analyzed, we move to the device level to solve the governing equations. These are cast in dimensionless form to establish a nondimensional framework for the analysis of the efficiency of the thermal energy input, energy storage, and energy retrieval processes as well as repeated cycling of the TES device. The framework is similar in some ways to the effectiveness-NTU methodology for analyzing heat exchanger performance, but is constructed to capture the inherent transience of these processes. This transience arises not only due to the time-varying nature of phase change but also due to a spatially and temporally varying working fluid temperature within the device. To capture these physics, we analyze a subsystem composed of the TES as well as a coupled heat exchanger that provides the time-varying working fluid inlet boundary condition for the extraction and charging processes. With thermodynamic effectiveness modeling in place for the entire subsystem, the TES unit can be designed to accommodate a variety of heating and cooling applications that would benefit from heat rejection load shifting.

The proposed TES device can be used to gather energy during the day, store it, and then reject it asynchronously (at night) rather than continuously throughout the day as is standard in many applications. This is particularly useful for a number of reasons. First, asynchronous cooling removes the need for peak load production, which is typically more expensive both financially and environmentally. Second, temperature differences are greater at night when ambient air is cooler, rendering heat transfer more thermodynamically efficient (from a simple Carnot standpoint). The type of cold storage considered here can be used for asynchronous cooling in a number of applications including rural refrigeration, building air conditioning, and steam power plants.

Figure 2 depicts a thermal energy storage device paired with an external heat exchanger, one connected in open loop to a heat source and the other to a heat sink. A thermal energy storage device is connected via a closed loop that circulates through the heat exchanger. To visualize how this works, we can consider a simple cycle. Ideally, we would like to begin with a completely frozen device. To achieve this, we send a cold working fluid (m˙sink>0) from a sink (Tsink,in) through the heat exchanger to cool a counterflowing fluid entering the thermal energy storage unit. Provided that the temperature of the fluid entering the TES device is less than the melt temperature, Tm, the PCM will undergo freezing. At a later time after a quiescent storage period (m˙closed=0), a hot working fluid (m˙source>0) from a source (Tsource,in) can be chilled by sending it through the heat exchanger, delivering heat to the closed loop fluid entering the thermal storage unit. This warm working fluid will reject heat to the cold storage matrix, provided its temperature is greater than Tm, thereby melting the PCM and chilling the closed-loop working fluid. This will chill the open-loop working fluid (Tsource,out) which can be used to augment cooling in various applications. With the subsystem described earlier, the heat exchanger could be connected to a Rankine cycle or air conditioning condenser in order to precool the open loop fluid and reduce the temperature at which steam or another refrigerant is condensed. At a later time, when the temperature difference with ambient is more favorable, the energy collected from precooling could be released as heat via the same heat exchanger. In short, the idea is to use cold storage to decrease the low system temperature of the power or refrigeration cycle (by heat rejection load shifting) without burning fuel or wasting water.

Analysis Framework

Thermal Energy Storage Energy Balance.

We assume that the working fluid in the TES unit either flows through a single passage or the flow is manifolded to multiple identical passages. Turns are ignored here for the purposes of this analysis and we instead focus on a unit cell of one long passage, with the mass flow rate per passage designated as m˙closed. The unit cell, of length dz, is composed of the working fluid flow passage and the surrounding PCM section. This element includes the tube wall and fin structures that conduct heat into the PCM material.

To derive the governing equations, we start with a control volume analysis, with one control volume around a differential element of the PCM matrix and another around a differential section of the flow passage. We apply the conservation of energy to each of these control volumes, noting that the stored thermal energy must be balanced by the thermal energy transport across the control surfaces of each unit cell.

From the definition of enthalpy for a solid–liquid mixture, thermal energy in the PCM matrix can be stored in sensible and latent forms [9]. The differential change in enthalpy of the storage element is designated as dHe 
dHe=ρ¯eνdz[c¯pe(TeTref)+xehls]
(1)
where ρ¯eνdz is the mass of PCM, c¯pe(TeTref), and xehls reflects its sensible and latent thermal energy, respectively. The above expression is differentiated with respect to time, and an energy balance is written between change in stored enthalpy of the element and the convective and conductive heat transfer across control surfaces to this element [10]. We use Newton's law of cooling as the constitutive equation to describe what form the heat transfer takes across the control surface separating the PCM matrix and flow passage 
(dHe)t=Uswdz(TwTe)
(2)
where Uswdz is the conductance from the working fluid to the PCM matrix and (Tw − Te) is the driving temperature difference. Taking the time derivative of the right-hand side of Eq. (1) and dividing through by constants ρ¯eνdzc¯pe, we can group like terms to formulate a governing equation relating stored enthalpy change to conductive and convective heat transfer into or out of the unit cell 
Tet+xet(hlsc¯pe)=Uswρ¯ec¯peν(TwTe)
(3)
In Eq. (3), at a given location, either the element temperature, Te (sensible heat transfer), or the melt fraction, xe (latent heat transfer), can change with time, but not both simultaneously. Therefore, the governing differential equation (3) can be split into the following two forms for sensible and latent heat transfer: 
Tet=Uswρ¯ec¯peν(TwTe);xet=0
(4)
for Te ≠ Tm and xe = 0 or xe = 1 
xet=Uswρ¯ehlsν(TwTe);Tet=0
(5)
for Te = Tm and 0 < xe < 1.
Likewise, conservation of energy on a control volume around the working fluid (inside the flow passage) of the unit cell requires that stored energy must be balanced by advection as well as heat transfer to and from the PCM matrix 
Twt=(m˙closedρwAc)Twz+UswρwAccpw(TeTw)
(6)

These energy balances neglect conduction in the downstream direction. One can show that the ratio of streamwise conduction to transport to or from the PCM is small for TES designs of interest. In other words, the heat diffusion effect is small compared to convection and conduction normal to the flow passage walls for the configuration considered here. These coupled equations are first-order in time and space, necessitating initial conditions for temperatures and melt fraction as well as a boundary (inlet) condition for the working fluid temperature.

Thermal Energy Storage Governing Equations.

Drawing upon the previous work, we can formulate the governing energy balance equations using dimensionless numbers and use these to resolve the temperature and melt fraction fields within the TES device [11]. The transport equations are written in terms of dimensionless position, time, and temperature. Position is nondimensionalized by dividing by the length of the flow passage. Time is nondimensionalized by dividing by the residence time of the working fluid in the storage device. Temperature is nondimensionalized by dividing by the maximum temperature difference experienced in the system. This is done because the TES absorbs heat in one range and rejects it in another, dictating the minimum and maximum temperatures for the system during the cycle. With this constraint, the melt temperature, Tm, should fall somewhere between Tmin and Tmax 
θ=TeTminTmaxTmin,ϕ=TwTminTmaxTmin
(7)
 
ẑ=zL,t*=ttres,tres=ρwAcLm˙closed
(8)
Replacing dimensional parameters with their dimensionless counterparts converts the equations to 
ϕt*=ϕẑ+Ntu(θϕ)
(9)
for sensible heat transfer in the flow passage and 
θt*=NtuRwe(ϕθ);xet*=0
(10)
for θ ≠ θm and xe = 0 or xe = 1 
xet*=NtuRweStio(ϕθ);θt*=0
(11)
for θ = θm and 0 < xe < 1. For sensible and latent heat transfer in the storage element where 
Ntu=UswLm˙clcpw,Rwe=ρwcpwAcρ¯ec¯peν,Stio=c¯pe(TmaxTmin)hls

The parameters in the nondimensional governing transport equations are dimensionless groups with physical relevance. Ntu, the number of transfer units, is used to specify dimensions and quantify the heat transfer rate associated with different designs. Rwe, the ratio of thermal capacities, can be used to select an appropriate PCM. The Stefan number, Stio, indicates whether heat transfer will be primarily sensible or latent.

The partial differential equations (PDEs) (9)(11) are solved numerically via an explicit finite difference discretization of the domain. The derivatives are replaced with forward difference algebraic expressions. For this finite difference discretization of derivatives, the initial and boundary conditions can be specified to model the cycling of the storage device. In particular, we can set initial conditions for extraction or charging to match the thermally equilibrated state after periods of storage. Furthermore, the transient inlet boundary condition can account for coupling with an external heat exchanger, and the varying temperatures and fluctuating heat transfer rates experienced there 
Att*=0:θ=θ(ẑ),xe=xe(ẑ),ϕ=ϕ(ẑ)
(12)
for 0ẑ1 
Atẑ=0:ϕ=ϕ(t*)
(13)
for t* > 0.

With any numerical method, it is imperative to ensure both stability and consistency of results. Due to the explicit nature of the forward Euler and upwind schemes, stability is only attainable below a threshold value of Δt* and Δẑ. In order to increase the time-step and grid size (thereby reducing computation time), it would be necessary to write derivatives as well as solve the governing equations differently. This can be accomplished via established methods for advective equations (e.g., Lax–Wendroff), or, the application of parabolic problem solutions to hyperbolic systems (e.g., Crank–Nicolson).

In order to quantify how well the storage device works, we define a performance metric for the extraction and charging loops. Effectiveness, εtes, is the calculated melt fraction at the last time-step to the melt fraction if the PCM is either completely melted or solidified.

Extraction effectiveness is defined as 
εtes,ext=ẑ=01xe(tend*,ẑ)ẑ=01xe,max(t*,ẑ)wherexe,max=1
(14)
where we neglect the contributions via sensible storage because, for the applications with small operating temperature differences considered here, latent storage dominates. This is consistent with the low Stefan number approximation adopted in this paper. Charging effectiveness, like extraction effectiveness, is defined as 
εtes,char=ẑ=011xe(tend*,ẑ)ẑ=011xe,min(t*,ẑ)wherexe,min=0
(15)
It follows that the effectiveness is a function of the dimensionless parameters in the model 
εtes=εtes(t*,Ntu,Rwe,Stio)
(16)

where the functional relation for εtes is embodied in the solution of the dimensionless differential equations. The overall effectiveness of the device should be taken as the minimum between extraction and charging and could be improved by adjusting working fluid mass flow rate and operation time (t*) for a given design geometry (Ntu), phase change material (Rwe), and operating conditions (Stio).

Coupled Heat Exchanger Energy Balance.

For the following analysis, we will use the well-established effectiveness-NTU method which posits that, for any given heat exchanger [10] 
(m˙cp)(ΔT)=Cmin(Thot,inTcold,in)εhx
(17)

With this equation, as well as the device modeling established earlier, we can examine the subsystem consisting of the TES as well as external heat exchanger.

Heat Exchanger During Extraction.

Conservation equation(s) 
Q˙source=(m˙cp)source(Tsource,inTsource,out)
(18)
 
Q˙source=(m˙cp)ext(Text,outText,in)
(19)
Performance equation(s) 
Q˙source=Cmin(Tsource,inText,in)εsource
(20)
Transport equation(s) 
m˙source=Q˙condcp,source(TcondTsource,out)εcond
(21)
For most relevant applications of this subsystem, the chilled open loop fluid exiting the heat exchanger would be used to condense a working fluid (e.g., steam, refrigerant) in a power or refrigeration cycle. Assuming that the goal of our subsystem is to precool a condenser, we begin by specifying the heat transfer rate in the condenser, Q˙cond, to solve for the open loop mass flow rate through the heat exchanger during extraction, as in Eq. (21). To solve this, we will have to begin by estimating Tsource,out, our target for the open loop fluid exiting the heat exchanger. If we solve Eq. (20) for Text,in, we find that 
Text,in=Tsource,inQ˙sourceCminεsource
(22)
where Cmin is either (m˙cp)source or (m˙cp)ext. We can solve Eqs. (19) and (20) for the unknown temperature, Text,out, but must specify whether Cmin is associated with the open loop or closed loop fluid 
Text,out=εsoCmin(m˙cp)extTso,in+(1εsoCmin(m˙cp)ext)Text,in
(23)

It is hard to predict ahead of time which fluid will be associated with Cmin as m˙ext is unknown. We need another functional relationship to establish a value for this flow rate. To do that, we solve the governing TES partial differential equations while iterating through m˙ext until the desired εtes is achieved. The only other unknown in the PDEs is U which is derived analytically from a Stefan type formulation [12].

Thermal Energy Storage Device Within Subsystem.

The total energy stored in the TES device can be determined in two ways. First, purely looking at the energy the PCM is capable of receiving or rejecting through latent heat transfer 
Ecap=ρpcmνLhlsεtes
(24)
The total energy stored must not exceed the energy capacity, Ecap, determined by the TES design. If the time required for the extraction or charging process is specified, knowledge of the heat transfer rate in the source heat exchanger enables us to compute the maximum amount of energy we can transfer to the TES device for a proposed operation time 
Etot=0textQ˙(t)dtQ˙sourcetext
(25)
As the goal is to eventually reject all of this as heat during asynchronous cooling, we can calculate the heat transfer rate in the sink heat exchanger as well 
Q˙sink=0tchardEdtEtottchar
(26)

This equation, in integral form, is challenging to solve for a time-varying heat transfer rate within the sink heat exchanger.

Heat Exchanger During Charging.

Conservation equation(s) 
Q˙sink=(m˙cp)sink(Tsink,outTsink,in)
(27)
 
Q˙sink=(m˙cp)char(Tchar,inTchar,out)
(28)
Performance equation(s) 
Q˙sink=Cmin(Tchar,inTsink,in)εsink
(29)
Transport equation(s) 
m˙sink=Q˙sinkcp,sink(Tsink,outTsink,in)
(30)
These relations have several unknowns, including m˙char,m˙sink,Tchar,in,Tchar,out,andTsink,out, and which fluid will have Cmin. For the charging process, we can confidently assume that Cmin is associated with the closed loop fluid (m˙cp)char such that ΔT will be greater within the TES and better drive the heat transfer necessary to refreeze the PCM. With this assumption, we can proceed by solving Eq. (29) for Tchar,in such that 
Tchar,in=Q˙sinkεsinkCmin+Tsink,in
(31)
We have justified the assumption that Cmin=(m˙cp)char, but m˙char remains unknown. Here, we can take the approach discussed previously to solve the governing TES PDEs while iterating through flow rates until εtes is achieved. In this way, we can define m˙char. In order to solve for Tchar,in, we must approximate Q˙sink. To do this, we will employ Eq. (26) and divide the total energy stored in the TES by the charging process operation time to determine an average Q˙sink. Tchar,out can be found from the conservation equation in terms of an average Q˙sink or from equating the conservation Eq. (28) and performance Eq. (29) and solving for Tchar,out 
Tch,out=εsiCmin(m˙cp)chTsi,in+(1εsiCmin(m˙cp)ch)Tch,in
(32)

The last remaining unknowns are Tsink,out and m˙sink, but it quickly becomes evident that we have no other relation to use in order to define values for either of these variables. In order to proceed, we will do as we did before and specify Tsink,out. This is more arbitrary than specifying Tsource,out as a desired temperature for the condenser inlet.

In specifying Tsink,out, we have several necessary conditions we have to meet. First, we must ensure that whatever choice of Tsink,out, (m˙cp)sink should remain greater than (m˙cp)char. We can think about our temperature constraints from a conservative, simplified perspective 
ΔTsink,min=Tmmax(Tsink,in)
(33)
For the heat exchanger to function as desired throughout the charging process, the open loop fluid temperature difference, (Tsink,out − Tsink,in), must be less than ΔT at all times. This is guaranteed to be the case if (Tsink,out − Tsink,in) is less than ΔTsink,min. Furthermore, we can assert that the following inequality should hold: 
Tsink,in<Tchar,out<Tsink,out<Tchar,in
(34)

We can use any guess for Tsink,out to set m˙sink so long as these requirements are met.

With relations in place for heat transfer rates, temperatures, and mass flow rates in the coupled heat exchanger and TES, we can proceed with a nondimensional formulation for the subsystem.

Subsystem Governing Equations.

Conservation equation(s) 
ϕclosed,outϕclosed,in=Q˙hx(m˙cp)closed(TmaxTmin)
(35)
Performance equation(s) 
εhx(ϕopen,inϕclosed,in)=Q˙hxCmin(TmaxTmin)
(36)
The dimensionless formulation established for the thermal energy storage device uses high and low system temperatures to nondimensionalize both the working fluid and element temperatures. We can also use this for the external heat exchanger during transient operation. In a realistic device, external temperatures would fluctuate throughout the day. These are input into the subsystem model as Topen,in in order to determine temperatures in the external heat exchanger, including the inlet boundary conditions to the TES device, Tclosed,out. This is no different than the analysis done for the source and sink heat exchanger above except that it is generalized to be applicable to either. Tclosed,out can be subsequently nondimensionalized to serve as ϕbc = ϕcl,out 
ϕcl,out=Cmin(m˙cp)clεhxϕop,in+(1Cmin(m˙cp)clεhx)ϕcl,in
(37)
Equation (37) is applicable to either heat exchanger and serves as a transient boundary condition for the inlet to the TES if Topen,in is a function of time. Alternatively, we could start with a nondimensional form and combine Eqs. (35) and (36) and solve for εhx 
εhx=(m˙cp)closedCmin(ϕclosed,outϕclosed,inϕopen,inϕclosed,in)
(38)
By similar heat exchanger analysis, the effectiveness of the TES can be written as 
εtes=ϕclosed,outϕclosed,inϕclosed,outθm
(39)
Here, the problem of specifying Cmin is averted by only having one fluid flowing through the TES. Equation (39) is applicable if the device is primarily undergoing latent heat transfer, as is the case for low Stefan number designs of interest here. These are of interest because the subsystem is intended for operation in a typical daytime ambient temperature range, which might encompass a maximum temperature difference of 20 °C, containing the fusion temperature somewhere between. As is shown in the case study example, this relatively small temperature difference renders the Stefan number much less than one which is appropriate for such TES subsystems. Equation (39) can be solved for ϕclosed,in 
ϕclosed,in=(1εtes)ϕclosed,out+εtesθm
(40)
This result can be substituted into Eq. (38) and ϕclosed,out can be found 
ϕcl,out=εhxϕop,inCmin+((m˙cp)clεhxCmin)εtesθmεhxCminεhxεtesCmin+εtes(m˙cp)cl
(41)

As above in Eq. (37), ϕbc = ϕclosed,out. Here, in Eq. (41), the boundary condition takes a very different form and incorporates information about the TES. Knowing which form to use comes down to the numerical scheme employed as well as an understanding of the assumptions in place. Equation (37) requires knowledge of both ϕopen,in and ϕclosed,in, the fluid temperature exiting the TES device. This equation makes sense to use with finite difference schemes in which the temperatures are defined throughout the device at all time steps. Conversely, Eq. (41) requires knowledge of ϕopen,in and θm, both of which are inputs to the programming. While this is simpler to use numerically, it should be employed with caution; this equation is only applicable during latent heat transfer and cannot be used if the PCM near the inlet is not at its melt temperature.

Performance at Short Times and Validation Via Closed Form Solution

We consider solution of the TES performance equations for a TES unit starting uniformly at temperature T0 = Te = Tw = Tm and initial melt fraction xe = xe,0. The equations will be nondimensionalized with Tmin = min(T0, Twi) and Tmax = max(T0, Twi) substituted into the previous definition for θ, ϕ, and Stio. Because of the initial conditions, θ = θm = 0 everywhere, and at all times of interest here, modifying the differential equation for ϕ to 
dϕdẑ=Ntuϕ
(42)
Note that before the transient starts, ϕ = θm = 0 everywhere. During the first residence time (ttres or t* ≤ 1), a wave-like step change in ϕ propagates downstream. We can integrate the above differential equation from the inlet to the local z location of interest noting that the definition of ϕ dictates that ϕ = 1 at ẑ=0. Integration of the equation leads to 
ϕ=eNtuẑ
(43)

This is the solution for ϕ for all times beyond the step-change front passage (as long as θ = θm = 0 everywhere). During the transient, each segment will only begin transferring heat and changing the melt fraction after the front passes when t*ẑ.

For these conditions, the energy balance equation for the PCM matrix reduces to 
xet*=±NtuRweStioϕ
(44)

where (+) applies to extraction and (−) applies to charging. This sign convention is used to properly account for the increasing or decreasing PCM melt fraction over time.

Substituting for ϕ using Eq. (43) and integrating this equation with respect to time using the initial condition of xe = xe,0 at a given ẑ location for time t*=ẑ yields 
xexe,0=±NtuRweStio(t*ẑ)eNtuẑ
(45)
Thus, for this either process, the solutions are 
ϕ=0;xe=xe,0
(46)
before the step change front passes (t*<ẑ) and 
ϕ=eNtuẑ;xe=xe,0±NtuRweStio(t*ẑ)eNtuẑ
(47)
after the step change front passes (t*ẑ). For the melting transient, these solutions are valid until xe = 1 at the inlet location, ẑ=0. Setting xe = 1 and ẑ=0 and solving for t* yields the t* value tend* at which the validity of the model ends for extraction 
tend*=1xe,0NtuRweStio
(48)
For the freezing transient, the heat flow is in the opposite direction (from the PCM to the fluid). During charging, these solutions are valid until xe = 0 at the inlet location, ẑ=0. Setting xe = 0 and ẑ=0 and solving for t* yields the t* value tend* at which the validity of the model ends 
tend*=xe,0NtuRweStio
(49)

For these short times of interest, we can compare the analytical and computational solutions in order to validate the numerical methods employed in this paper.

The working fluid temperature, ϕ, is fairly consistent between the numerical (approximated) and the analytical (exact) solution as seen in Fig. 3. The two differ in treatment of the step change front. The exact solution exhibits a vertical drop in temperature, as the information from upstream has not yet reached the location downstream. The numerical scheme cannot handle such discontinuities and smears the solution around the shock. This is purely an artifact of the numerics and can be improved by using a finer grid, or employing a different numerical scheme while solving the differential equations.

Despite the error introduced in ϕ while attempting to capture the physics of the step change front, the solution of xe is compelling. As evidenced in Fig. 4, the numerical and analytical solutions are almost indistinguishable from each other. This comparison with the closed form analytical solution provides necessary validation of our numerical method as well as understanding of the physics governing the TES device performance.

Parametric Studies

As mentioned earlier in the derivation of the governing equations, the values of dimensionless groups t*, Ntu, Rwe, and Stio, are embodied in the effectiveness prediction for a potential TES unit. To better understand this functional relationship, we can examine the parameters individually and establish desired ranges of these for high effectiveness designs.

Variation of Ntu.

In examining these dimensionless groups, we will start with Ntu which incorporates the overall heat transfer coefficient, U, in its definition. The analytical framework is set up for a constant Ntu based on a mean U¯ value, which must be averaged both spatially and temporally as described in a previous paper [13]. With an average U¯ determined, we can look at how its corresponding dimensionless counterpart, Ntu, impacts effectiveness.

For specified values of the storage design parameters Rwe, Stio, and θm, computations can be done for different combinations of Ntu and dimensionless termination time, tend*, to determine the resulting effectiveness of the extraction and charging processes. Figure 5 illustrates the results of multiple solutions to define the dependence of the effectiveness on tend* and Ntu.

From this figure, we can draw conclusions on what a reasonable value for Ntu should be in order to achieve a target effectiveness (for a specified Rwe and Stio). For high efficiency TES devices, the recommended range of Ntu is between 1.0 and 12.0 for a reasonable operation time (tend*>15) and conditions (Rwe = 1, Stio = 0.1) as device operation below this threshold cannot achieve desired effectiveness targets. As Ntu=UswL/m˙wcpw, the recommendation that Ntu should be one or greater is a reflection that the heat transfer to/from the PCM matrix must equal or outweigh the working fluid's capacity to advect the energy along the flow passage.

In most design scenarios, there is a fixed time window for operation and an impetus to make the TES as small as possible while still accomplishing heat transfer during that time. This inevitably leads to a design trade-off between tend* and Ntu. If the process can take longer, the task can be done with a smaller heat exchanger, or conversely, if the TES device is larger, the process can take less time. In order to go about examining the design space for an effective TES device, there are a couple of insights that can be gained from Fig. 5. There are several regimes which are of great interest. At low Ntu (Ntu ≤ 1), the only way to achieve high performance is to greatly increase tend*. To do this would entail increasing the total operation time or decreasing the residence time in the flow passage. Lowering the residence time could be done by selecting a lower density fluid, decreasing the channel cross-sectional area, decreasing the length of the flow passage, or increasing the mass flow rate. At midrange Ntu (1.0 < Ntu < 12.0), there is some room to decrease tend* while increasing Ntu and vice versa. Ntu can be increased by enhancing the overall heat transfer, increasing the heat transfer area, and decreasing the mass flow rate or specific heat of the working fluid. At high Ntu (Ntu ≥ 12), we notice diminishing returns. Any further increase in Ntu will not make much difference, and the only relevant parameter to adjust is tend*. In this range of Ntu, once tend* is increased to the point of achieving a desired effectiveness value, any further increase in tend* is superfluous.

Variation of Rwe.

For specified values of the storage design parameters Ntu, Stio, and θm, computations can be done for different combinations of Rwe and dimensionless termination time, tend*, to determine the resulting effectiveness of the extraction and charging processes. Figure 6 illustrates the results of multiple solutions to define the dependence of the effectiveness on tend* and Rwe.

From this figure, we can draw conclusions on what a reasonable value for Rwe should be in order to achieve a target effectiveness (for a specified Ntu and Stio). For high efficiency TES devices, the recommended range of Rwe is between 0.4 and 1.8 for a reasonable operation time (tend*>15) and conditions (Ntu = 10, Stio = 0.1) as device operation below this threshold cannot achieve desired effectiveness targets. As defined, Rwe=ρwcpwAc/ρ¯ec¯peν, so the recommendation that Rwe be somewhere near unity is a reflection that the working fluid and PCM should have a similar energy capacity.

Variation of Stio.

For specified values of the storage design parameters Ntu, Rwe, and θm, computations can be done for different combinations of Stio and dimensionless termination time, tend*, to determine the resulting effectiveness of the extraction and charging processes. Figure 7 illustrates the results of multiple solutions to define the dependence of the effectiveness on tend* and Stio.

From this figure, we can draw conclusions on what a reasonable value for Stio should be in order to achieve a target effectiveness (for a specified Ntu and Rwe). For high efficiency TES devices, the recommended range of Stio is between 0.06 and 0.20 for a reasonable operation time (tend*>15) and conditions (Ntu = 10, Rwe = 1) as device operation below this threshold cannot achieve desired effectiveness targets. The Stefan number is written as Stio=c¯pe(TmaxTmin)/hls. This dimensionless group reflects the ratio of sensible to latent heat transfer. As such, it makes sense that for the thermal energy storage we aim to accomplish, employing phase change in a constrained temperature range, the latent heat transfer and resulting effectiveness should be high, rendering Stio less than one, consistent with the low Stefan number approximation adopted here.

Key Takeaways From Parametric Studies.

These parametric studies are particularly valuable in defining an appropriate design space for the TES. From a practical perspective, it is important to parse through the suggested ranges and make design decisions. The ratio of thermal capacities, Rwe, is set early on by selecting materials for the working fluid and PCM and is constrained by the properties of both of these. The Stefan number, Stio, is defined by the operating conditions of the system and is constrained by the ambient temperatures surrounding the subsystem. The least constrained terms in designing a TES device are Ntu and tend*. Optimizing effectiveness is a trade-off between these two parameters.

Case Study Example

To provide concreteness, the effectiveness modeling of the thermal energy storage will be used to predict device performance. The case study considers a TES unit with a nominal capacity of 1.5 TJ, which would be about the right size to provide asynchronous cooling for a 100 MW power plant. In order to delve into a detailed analysis of the TES unit, operating conditions for the device must first be specified. We prescribe a target effectiveness for the TES as well as desired times for operation. We assume that the extraction process, removing heat from the working fluid and melting the PCM, would occur during a 10 h period. For our preliminary study, we consider that the extraction process begins with the TES completely frozen. The initial temperature of the PCM during extraction, T0,ext, is assumed to be at the melting temperature.

We assume that the charging process of refreezing the PCM would occur during a 9 h period following a 3 h period of storage. The initial temperature of the PCM during charging, T0,char, is close to the final temperature distribution after extraction Tf,ext(z), but slightly shifted due to thermal equilibration during storage. This distribution exists because the working fluid and the element temperature vary along the flow passage as heat is transferred between the working fluid and the PCM. Only where latent heat transfer is taking place can we expect the element temperature to be Tm.

However, during the extraction process, the PCM near the inlet undergoes and completes phase change more quickly than the PCM further downstream. After complete melting, this PCM near the inlet still accepts heat in the form of sensible heat transfer raising its temperature.

When the extraction process stops, energy is stored in the device until a more favorable temperature difference with ambient can be achieved. During storage, the same set of governing equations is solved, with the mass flow rate set to zero, eliminating the first-order advection term from the working fluid equation. The temperature and melt fraction distributions change slightly due to the driving normal temperature gradient between the PCM matrix and the working fluid. As thermal diffusivities of both the PCM (∼10−5 m2/s) and the working fluid (∼10−7 m2/s) are quite small, and the storage time is relatively brief, axial conduction along the passage is neglected as in the governing partial differential equations. If axial conduction were included, the thermal diffusivities of the working fluid and the PCM matrix would serve as coefficients to their respective second-order temperature diffusion terms. Considering the larger of these, it seems that the relevant time scale for axial conduction in the PCM matrix is something like L2/α (∼16,500 s). This is roughly double the storage time considered here (∼7200–10,800 s). If the materials in the TES unit have higher thermal diffusivities, this modeling framework could be modified to include the effects of axial conduction. Retaining the original set of governing equations for storage, T0,c will be a function of the axial coordinate, z, as reflected in Table 1. The working fluid inlet temperatures for extraction and charging must also be defined. These temperatures are the boundary conditions for the TES device. These selected baseline conditions are enumerated in Table 1.

For this case study, the thermal storage device is assumed to be a plate-fin heat exchanger with alternating layers of phase change material and rectangular working fluid flow passages on the liquid side. For the case study considered here, we carefully selected our working fluid and phase change material. A 50/50 mixture of ethylene glycol and water is ideal in this system because of its applicability in a wide range of operating temperatures. This working fluid has also been combined with anticorrosion agents for many years to reduce fouling in heat exchangers. Lithium nitrate trihydrate (LNT) is an ideal candidate for phase change material due to its high energy density. This PCM also has desirable thermal conductivity and cost. Furthermore, it is chemically nonreactive and stable. The melting temperature of the PCM was taken to be 30 °C, which is possible with a specific variation of LNT [14]. Heat transfer into the PCM can be enhanced using metal structures; effective thermal properties of the matrix containing LNT and aluminum fins are used in the model.

Dimensional parameters for the case study are used to calculate the dimensionless variables in Table 2. Some variables, including Ntu, differ for the extraction and charging processes because the mass flow rate and the overall heat transfer coefficient differ. Others such as Rwe and Stio remain consistent. If the dimensionless variables are changed, the simulation will naturally lead to different results.

Extraction/Melting.

During cold extraction, heat is delivered to initially frozen PCM in the storage device. The operation time and the mass flow rate through the channel are chosen to give rise to dimensionless parameters that make it possible to achieve a prescribed TES effectiveness. The dimensionless parameter values for extraction from Table 2 were input into the three governing partial differential equations (9)(11) to determine melt fraction and temperature throughout the element material.

Computed solutions of the governing equations were generated using an explicit finite difference scheme with forward time and upwind spatial finite difference representations for derivatives in the equations [15]. Reported performance calculation results are for a time-step Δt* of 0.0025 and a spatial mesh Δẑ of 0.005. Reducing these by a factor of 5 produced less than 1% change in the computed effectiveness, indicating that the solutions are not sensitive to the choices of Δẑ and Δt* at this threshold or that the solution is convergent [16].

For the dimensionless parameter values listed in Table 2, the spatial variations of dimensionless working fluid and storage element temperatures, and melt fraction are shown for 3 moments in time in Fig. 8(a). As the transient proceeds, more of the storage raises in temperature toward the inlet working fluid temperature, and an increasing fraction of the PCM is melted. In the first snapshot, the dimensionless working fluid temperature enters the TES at a hot working fluid inlet temperature. As it travels through the device, it rejects heat through the channel walls to the PCM and exits at a lower temperature. The dimensionless element temperature begins to increase due to the hot working fluid. The melt fraction has begun its ascent from a solid toward a liquid state. At a later time, the hot working fluid has heated the element temperature up as well as almost melted the PCM to its dimensionless liquid state of 1.

Charging/Freezing.

During the cold charging operation, the cooling working fluid loop is activated to refreeze the PCM. The set of partial differential equations (9)(11) is also solved for charging the TES device. The dimensionless parameter values calculated from the baseline case for charging are listed in Table 2.

Figure 8(b) presents three snapshots in time of the cold storage charging process. In the first snapshot, the dimensionless working fluid temperature enters the TES at a cold working fluid inlet temperature. As it travels through the device, it accepts heat through the channel walls from the PCM and exits at a higher temperature. The dimensionless element temperature begins to decrease due to the cool working fluid. The melt fraction has begun its descent from a saturated liquid toward a solid state. At a later time, the cool working fluid has almost cooled the element temperature down as well as frozen the PCM toward its dimensionless solid state of 0.

Subsystem Full Cycle.

Moving forward, we are able to fully take into account the cycling of the device. We can visualize the paired extraction and charging processes by looking at the entire cycle as shown in Fig. 9.

Figure 9, depicting the 24 h cycle schematic, has a lot of information to digest. First, we can look at the top left box: extraction (1). Here, hot ambient air passes through a heat exchanger referred to as the precooler which cools air by transferring heat to a closed loop working fluid. This closed loop working fluid flows through the TES, transfers heat to the PCM matrix, and re-emerges to continue cooling the flowing air. This cooled dry air is subsequently used to condense steam in the condensate manifold. The precooling/extraction process continues for the desired duration specified by plant operators, transferring heat to melt the amount of PCM prescribed by the TES effectiveness. After this period of extraction, we move into the next process indicated by the top right box: storage (2). During this time, no fluid flows through the TES device. Instead, the device thermally equilibrates such that the fluid and the adjacent PCM eventually achieve the same temperature. If this temperature equilibration causes the PCM to reach its melting point, it will undergo further melting as dictated by the driving temperature difference. After this period of storage and equilibration, the PCM in the TES device needs to be frozen again for future use in precooling. To accomplish this, we move into the charging process, depicted in the bottom right box (3). At this time, the sink heat exchanger is employed in order to transfer heat from the PCM matrix within the thermal storage to the atmosphere. A cool open loop fluid passes through the sink heat exchanger, removing heat from the closed loop fluid. This chilled fluid enters the TES device, and, because it flows through at a lower temperature than Tm, enables the PCM to begin freezing. This closed loop fluid exits the TES warmer than when it enters, but is subsequently rechilled by the open loop fluid in the night cooler. This process continues for a specified amount of time, which should be enough to return the TES to its desired starting point, or melt fraction distribution, for later precooling to take place. Before that happens, the device undergoes another period of storage, shown in the bottom left (4). Once again, the device thermally equilibrates and is ready to recommence the cycle.

Further cycling of the device should lead to a steady-state operation with a more symmetric distribution of average melt fraction through the paired extraction and charging processes. The final charging condition will serve as the initial extraction condition after storage and vice versa. We specify operating conditions so that the device remains in this window, melting and freezing the fraction of PCM according to its effectiveness target.

For the applications of interest discussed in the introduction, we are interested in device modeling that accounts for the time-varying temperatures and boundary conditions encountered in real systems. In order to understand what happens during this asynchronous cooling cycle, we can look at Fig. 10. This figure contains the relevant results obtained by solving the TES device and subsystem governing equations.

The open loop working fluid (air) inlet temperatures were approximated as parabolic in nature, reflecting the variation we might see during a typical day. The open loop air exit temperatures were not specified. Instead, we solve the heat exchanger performance equation for Topen,out 
Top,out=Top,in(1Cmin(m˙cp)opεhx)+Tcl,inCmin(m˙cp)opεhx
(50)

where Topen,in is known, Tclosed,in is calculated from the TES PDEs, m˙open is determined from an average Q˙hx, and m˙closed is found via iteration for a prescribed effectiveness. If Cmin is (m˙cp)closed, then the computation for m˙open requires knowledge of Topen,out as in Eq. (21). To calculate an appropriate constant flowrate, m˙open, we initially approximate Topen,out and finally solve for it here.

We find that the air temperatures depicted in Fig. 10 satisfy the requirement that they are higher than Tm during extraction and lower than Tm during charging. This ensures that the heat transfer occurs in the appropriate direction for both processes—melting the PCM during the day and freezing the PCM overnight. This is reflected by the spatially averaged melt fraction throughout the entire 24 h cycle.

Most importantly, this figure highlights the usefulness of thermal energy storage for these applications. Load shifting, storing energy during the day and subsequently rejecting it at night, takes advantage of favorable temperature differences that are a function of a time-varying ambient temperature. Furthermore, the lower air inlet temperature to the steam condenser improves the Rankine cycle efficiency, which is highly beneficial both financially and environmentally.

Conclusions

The dimensionless effectiveness-NTU analysis framework developed here for the thermal storage and coupled heat exchanger can be a useful tool for design optimization of the TES unit and the associated subsystems for asynchronous cooling and other thermal storage applications involving transient storage and retrieval processes. This formulation also defines the key dimensionless parameters that dictate performance and facilitates the analysis to define optimal ranges of these parameters.

After deriving the relevant conservation and performance equations, we examined various parameters of a thermal energy storage device to determine an effective and efficient design space. We used these parametric studies to define the required combination of dimensionless parameters for a high performing device.

We indicate which dimensionless parameters affect performance as well as the relationship between them. Rwe relates the thermal capacities of the working fluid to typical phase change materials. This parameter can be optimized between ∼0.4 and 1.8 by selecting a working fluid with ideal thermal properties as well as a low cost phase change material with high energy density. Stio relates sensible to latent heat transfer in the phase change material. This can be varied by adjusting the initial condition of the TES unit to incorporate greater subcooling or superheating. This study, instead, focused primarily on demonstrating the usefulness of this device by defining a performance metric, εtes, that quantified how much latent heat transfer via phase change occurred rendering Stio between ∼0.06 and 0.20. Ntu corresponds most directly to heat exchanger effectiveness as it encapsulates transport parameters including m˙closed and U; it is optimal between ∼1.0 and 12.0. It is apparent that there is a threshold dimensionless time, t*, in order to reach prescribed effectiveness targets. Optimizing design of a TES device and subsystem often requires a trade-off between tend* and Ntu. For different Ntu, Rwe, and Stio, the time required to complete the prescribed amount of latent heat transfer varies. Once the required values of Ntu, Rwe, Stio, and tend* are defined, detailed design of the thermal storage unit can be focused on the goal of establishing a flow passage and PCM packaging design that achieves the required values of these parameters. This framework thus can be used to optimize the thermal energy storage device for repeated cycling in load shifting applications.

In order to accomplish heat rejection load shifting, we examined the TES device in the context of a subsystem, paired to a heat exchanger that is used to input and later reject energy from the TES. In the case study examined in this paper, we demonstrated that the model has the capability to predict performance for time-varying working fluid inlet temperature to the heat exchanger and thus the TES. Unlike earlier models, the time-varying inlet temperature case is significantly more complicated because it introduces time-varying boundary conditions for the differential equations in the model. Our numerical approach enables us to handle this challenge well. Furthermore, we can confidently apply this numerical method having validated it via comparison to an exact analytical solution. The computational framework, described in-depth in this paper, can serve as a valuable tool to model and optimize time-varying thermal energy storage subsystems for a multitude of power and refrigeration applications.

Acknowledgment

Support for this research was provided by the ARID program of ARPA-E. We would also like to thank the reviewers for their helpful suggestions for improving this paper.

Funding Data

  • Advanced Research Projects Agency (DE-AR0000577).

Nomenclature

     
  • Ac =

    cross-sectional area of flow passage

  •  
  • cpw =

    working fluid specific heat

  •  
  • Cmin =

    minimum heat capacity rate in external HX

  •  
  • c¯pe =

    effective specific heat of differential element

  •  
  • Ecap =

    total energy capacity of the PCM in the TES device

  •  
  • Etot =

    total energy stored in the TES device

  •  
  • h =

    convective heat transfer coefficient

  •  
  • hchan =

    height of working fluid flow passage

  •  
  • hls =

    latent heat of fusion of PCM

  •  
  • hpcm =

    height of PCM matrix section

  •  
  • He =

    enthalpy of differential element

  •  
  • kw =

    working fluid thermal conductivity

  •  
  • k¯e =

    effective thermal conductivity of differential element

  •  
  • L =

    length of flow passage

  •  
  • m˙closed =

    working fluid closed loop mass flow rate through TES; may be written as m˙ext or m˙char

  •  
  • m˙open =

    working fluid open loop mass flow rate through external HX; may be written as m˙source or m˙sink

  •  
  • Q˙ =

    heat transfer rate in the external heat exchanger

  •  
  • sw =

    wetted perimeter of flow passage

  •  
  • tchar =

    time elapsed to end of charging process

  •  
  • text =

    time elapsed to end of extraction process

  •  
  • Tclosed =

    TES closed loop working fluid temperature; may be written as Text or Tchar

  •  
  • Te =

    TES device matrix element temperature

  •  
  • Tm =

    melting temperature of PCM

  •  
  • Tmax =

    maximum temperature in the system

  •  
  • Tmin =

    minimum temperature in the system

  •  
  • Topen =

    external HX open loop working fluid temperature; may be written as Tsource or Tsink

  •  
  • Tw =

    TES device working fluid temperature

  •  
  • T0 =

    initial storage temperature for extraction or charging

  •  
  • U =

    overall heat transfer coefficient between bulk working fluid and PCM matrix element

  •  
  • wchan =

    width of flow passage

  •  
  • xe =

    melt fraction of PCM in matrix element

  •  
  • εcond =

    effectiveness of precooled condenser

  •  
  • εhx =

    effectiveness of external heat exchanger; may be written as εsource or εsink

  •  
  • εtes =

    effectiveness of TES

  •  
  • θ =

    dimensionless matrix element temperature

  •  
  • μw =

    working fluid viscosity

  •  
  • ν =

    PCM matrix volume per unit flow length

  •  
  • ρpcm =

    density of PCM

  •  
  • ρw =

    working fluid density

  •  
  • ρ¯e =

    effective density of TES differential element

  •  
  • ϕ =

    dimensionless working fluid temperature

References

References
1.
Shamsundar
,
N.
, and
Srinivasan
,
R.
,
1980
, “
Effectiveness-NTU Charts for Heat Recovery From Latent Heat Storage Units
,”
ASME J. Sol. Energy Eng.
,
102
(
4
), pp.
263
271
.
2.
El-Dessouky
,
H.
, and
Al-Juwayhel
,
F.
,
1997
, “
Effectiveness of a Thermal Energy Storage System Using Phase-Change Materials
,”
Energy Convers. Manage.
,
38
(
6
), pp.
601
617
.
3.
Ismail
,
K.
, and
Goncalves
,
M.
,
1999
, “
Thermal Performance of a Pcm Storage Unit
,”
Energy Convers. Manage.
,
40
(
2
), pp.
115
138
.
4.
Sharma
,
A.
,
Tyagi
,
V.
,
Chen
,
C.
, and
Buddhi
,
D.
,
2009
, “
Review on Thermal Energy Storage With Phase Change Materials and Applications
,”
Renewable Sustainable Energy Rev.
,
13
(
2
), pp.
318
345
.
5.
Pielichowska
,
K.
, and
Pielichowski
,
K.
,
2014
, “
Phase Change Materials for Thermal Energy Storage
,”
Prog. Mater. Sci.
,
65
, pp.
67
123
.
6.
Tay
,
N.
,
Belusko
,
M.
, and
Bruno
,
F.
,
2012
, “
An Effectiveness-NTU Technique for Characterising Tube-in-Tank Phase Change Thermal Energy Storage Systems
,”
Appl. Energy
,
91
(
1
), pp.
309
319
.
7.
Hu
,
Y.
,
Li
,
D.
,
Shu
,
S.
, and
Niu
,
X.
,
2017
, “
Lattice Boltzmann Simulation for Three-Dimensional Natural Convection With Solid-Liquid Phase Change
,”
Int. J. Heat Mass Transfer
,
113
, pp.
1168
1178
.
8.
Bondareva
,
N. S.
, and
Sheremet
,
M. A.
,
2017
, “
3D Natural Convection Melting in a Cubical Cavity With a Heat Source
,”
Int. J. Therm. Sci.
,
115
, pp.
43
53
.
9.
Yunus
,
A. C.
, and
Michael
,
A. B.
,
2006
,
Thermodynamics: An Engineering Approach
,
McGraw-Hill
,
New York
.
10.
Bergman
,
T. L.
,
Lavine
,
A. S.
,
Incropera
,
F. P.
, and
DeWitt
,
D. P.
,
2011
,
Fundamentals of Heat and Mass Transfer
,
Wiley
, Hoboken, NJ.
11.
Helmns
,
A.
, and
Carey
,
V. P.
,
2016
, “
Modeling of Heat Transfer and Energy Efficiency Performance of Transient Cold Storage in Phase Change Thermal Storage Components
,”
ASME
Paper No. HT2016-7237.
12.
Lunardini
,
V. J.
,
1981
,
Heat Transfer in Cold Climates
,
Van Nostrand Reinhold Company
, New York.
13.
Helmns
,
A.
, and
Carey
,
V. P.
,
2017
, “
Modeling of Intramatrix Heat Transfer in Thermal Energy Storage for Asynchronous Cooling
,”
ASME
Paper No. HT2017-4870.
14.
Shamberger
,
P. J.
, and
Reid
,
T.
,
2012
, “
Thermophysical Properties of Lithium Nitrate Trihydrate From (253 to 353) k
,”
J. Chem. Eng. Data
,
57
(
5
), pp.
1404
1411
.
15.
Jaluria
,
Y.
,
2012
,
Computer Methods for Engineering With MATLAB Applications
,
Taylor & Francis, Inc
, Boca Raton, FL.
16.
LeVeque
,
R. J.
,
2007
,
Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems
,
SIAM
, Philadelphia, PA.