Multiscale Transient Modeling of Latent Energy Storage for Asynchronous Cooling

Helmns, Andrea; Carey, Van P.

doi:10.1115/1.4039460

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. [1–3]). 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. [6–8]).

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 N_tu, within the context of a subsystem that models realistic operating conditions for a TES.

Fig. 1

View large Download slide

Schematic depicting multiscale nature of problem ranging from the subsystem including the external heat exchanger to the TES device to the unit cell differential element (containing fins, PCM, and flow passage)

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 $({\dot{m}}_{sink} > 0)$ from a sink (T_sink,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, T_m, the PCM will undergo freezing. At a later time after a quiescent storage period $({\dot{m}}_{closed} = 0)$ ⁠, a hot working fluid $({\dot{m}}_{source} > 0)$ from a source (T_source,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 T_m, thereby melting the PCM and chilling the closed-loop working fluid. This will chill the open-loop working fluid (T_source,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.

Fig. 2

View large Download slide

Thermal energy storage coupled with an external heat exchanger for cold extraction (left) and cold charging (right)

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 ${\dot{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 dH_e

d H_{e} = {\bar{ρ}}_{e} ν' d z [{\bar{c}}_{pe} (T_{e} - T_{ref}) + x_{e} h_{ls}]

(1)

where

{\bar{ρ}}_{e} ν' d z

is the mass of PCM,

{\bar{c}}_{pe} (T_{e} - T_{ref})

⁠, and x_eh_ls 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

\frac{\partial (d H_{e})}{\partial t} = U s_{w} d z (T_{w} - T_{e})

(2)

where Us_wdz is the conductance from the working fluid to the PCM matrix and (T_w − T_e) is the driving temperature difference. Taking the time derivative of the right-hand side of Eq. and dividing through by constants

{\bar{ρ}}_{e} ν' d z {\bar{c}}_{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

\frac{\partial T_{e}}{\partial t} + \frac{\partial x_{e}}{\partial t} (\frac{h_{ls}}{{\bar{c}}_{pe}}) = \frac{U s_{w}}{{\bar{ρ}}_{e} {\bar{c}}_{pe} ν'} (T_{w} - T_{e})

(3)

In Eq. , at a given location, either the element temperature, T_e (sensible heat transfer), or the melt fraction, x_e (latent heat transfer), can change with time, but not both simultaneously. Therefore, the governing differential equation can be split into the following two forms for sensible and latent heat transfer:

\frac{\partial T_{e}}{\partial t} = \frac{U s_{w}}{{\bar{ρ}}_{e} {\bar{c}}_{pe} ν'} (T_{w} - T_{e}); \frac{\partial x_{e}}{\partial t} = 0

(4)

for T_e ≠ T_m and x_e = 0 or x_e = 1

\frac{\partial x_{e}}{\partial t} = \frac{U s_{w}}{{\bar{ρ}}_{e} h_{ls} ν'} (T_{w} - T_{e}); \frac{\partial T_{e}}{\partial t} = 0

(5)

for T_e = T_m and 0 < x_e < 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

\frac{\partial T_{w}}{\partial t} = - (\frac{{\dot{m}}_{closed}}{ρ_{w} A_{c}}) \frac{\partial T_{w}}{\partial z} + \frac{U s_{w}}{ρ_{w} A_{c} c_{pw}} (T_{e} - T_{w})

(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, T_m, should fall somewhere between T_min and T_max

θ = \frac{T_{e} - T_{min}}{T_{max} - T_{min}}, ϕ = \frac{T_{w} - T_{min}}{T_{max} - T_{min}}

(7)

\hat{z} = \frac{z}{L}, t^{*} = \frac{t}{t_{res}}, t_{res} = \frac{ρ_{w} A_{c} L}{{\dot{m}}_{closed}}

(8)

Replacing dimensional parameters with their dimensionless counterparts converts the equations to

\frac{\partial ϕ}{\partial t^{*}} = - \frac{\partial ϕ}{\partial \hat{z}} + N_{tu} (θ - ϕ)

(9)

for sensible heat transfer in the flow passage and

\frac{\partial θ}{\partial t^{*}} = N_{tu} R_{we} (ϕ - θ); \frac{\partial x_{e}}{\partial t^{*}} = 0

(10)

for θ ≠ θ_m and x_e = 0 or x_e = 1

\frac{\partial x_{e}}{\partial t^{*}} = N_{tu} R_{we} {St}_{io} (ϕ - θ); \frac{\partial θ}{\partial t^{*}} = 0

(11)

for θ = θ_m and 0 < x_e < 1. For sensible and latent heat transfer in the storage element where

N_{tu} = \frac{U s_{w} L}{{\dot{m}}_{cl} c_{pw}}, R_{we} = \frac{ρ_{w} c_{pw} A_{c}}{{\bar{ρ}}_{e} {\bar{c}}_{pe} ν'}, {St}_{io} = \frac{{\bar{c}}_{pe} (T_{max} - T_{min})}{h_{ls}}

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

The partial differential equations (PDEs) – 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

At t^{*} = 0 : θ = θ (\hat{z}), x_{e} = x_{e} (\hat{z}), ϕ = ϕ (\hat{z})

(12)

for

0 \leq \hat{z} \leq 1

At \hat{z} = 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 $Δ \hat{z}$ ⁠. 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} = \frac{\sum_{\hat{z} = 0}^{1} x_{e} (t_{end}^{*}, \hat{z})}{\sum_{\hat{z} = 0}^{1} x_{e, max} (t^{*}, \hat{z})} where x_{e, 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} = \frac{\sum_{\hat{z} = 0}^{1} 1 - x_{e} (t_{end}^{*}, \hat{z})}{\sum_{\hat{z} = 0}^{1} 1 - x_{e, min} (t^{*}, \hat{z})} where x_{e, min} = 0

(15)

It follows that the effectiveness is a function of the dimensionless parameters in the model

ε_{tes} = ε_{tes} (t^{*}, N_{tu}, R_{we}, {St}_{io})

(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 (N_tu), phase change material (R_we), and operating conditions (St_io).

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]

(\dot{m} c_{p}) (Δ T) = C_{min} (T_{hot, in} - T_{cold, 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)

{\dot{Q}}_{source} = {(\dot{m} c_{p})}_{source} (T_{source, in} - T_{source, out})

(18)

{\dot{Q}}_{source} = {(\dot{m} c_{p})}_{ext} (T_{ext, out} - T_{ext, in})

(19)

Performance equation(s)

{\dot{Q}}_{source} = C_{min} (T_{source, in} - T_{ext, in}) ε_{source}

(20)

Transport equation(s)

{\dot{m}}_{source} = \frac{{\dot{Q}}_{cond}}{c_{p, source} (T_{cond} - T_{source, 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,

{\dot{Q}}_{cond}

⁠, to solve for the open loop mass flow rate through the heat exchanger during extraction, as in Eq. . To solve this, we will have to begin by estimating T_source,out, our target for the open loop fluid exiting the heat exchanger. If we solve Eq. for T_ext,in, we find that

T_{ext, in} = T_{source, in} - \frac{{\dot{Q}}_{source}}{C_{min} ε_{source}}

(22)

where C_min is either

{(\dot{m} c_{p})}_{source}

or

{(\dot{m} c_{p})}_{ext}

⁠. We can solve Eqs. and for the unknown temperature, T_ext,out, but must specify whether C_min is associated with the open loop or closed loop fluid

T_{ext, out} = ε_{so} \frac{C_{min}}{{(\dot{m} c_{p})}_{ext}} T_{so, in} + (1 - ε_{so} \frac{C_{min}}{{(\dot{m} c_{p})}_{ext}}) T_{ext, in}

(23)

It is hard to predict ahead of time which fluid will be associated with C_min as ${\dot{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 ${\dot{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

E_{cap} = ρ_{pcm} ν' L h_{ls} ε_{tes}

(24)

The total energy stored must not exceed the energy capacity, E_cap, 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

E_{tot} = \int_{0}^{t_{ext}} \dot{Q} (t) d t \approx {\dot{Q}}_{source} t_{ext}

(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

{\dot{Q}}_{sink} = \int_{0}^{t_{char}} \frac{d E}{d t} \approx \frac{E_{tot}}{t_{char}}

(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)

{\dot{Q}}_{sink} = {(\dot{m} c_{p})}_{sink} (T_{sink, out} - T_{sink, in})

(27)

{\dot{Q}}_{sink} = {(\dot{m} c_{p})}_{char} (T_{char, in} - T_{char, out})

(28)

Performance equation(s)

{\dot{Q}}_{sink} = C_{min} (T_{char, in} - T_{sink, in}) ε_{sink}

(29)

Transport equation(s)

{\dot{m}}_{sink} = \frac{{\dot{Q}}_{sink}}{c_{p, sink} (T_{sink, out} - T_{sink, in})}

(30)

These relations have several unknowns, including

{\dot{m}}_{char}, {\dot{m}}_{sink}, T_{char, in}, T_{char, out}, and T_{sink, out}

⁠, and which fluid will have C_min. For the charging process, we can confidently assume that C_min is associated with the closed loop fluid

{(\dot{m} c_{p})}_{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. for T_char,in such that

T_{char, in} = \frac{{\dot{Q}}_{sink}}{ε_{sink} C_{min}} + T_{sink, in}

(31)

We have justified the assumption that

C_{min} = {(\dot{m} c_{p})}_{char}

⁠, but

{\dot{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

{\dot{m}}_{char}

⁠. In order to solve for T_char,in, we must approximate

{\dot{Q}}_{sink}

⁠. To do this, we will employ Eq. and divide the total energy stored in the TES by the charging process operation time to determine an average

{\dot{Q}}_{sink}

⁠. T_char,out can be found from the conservation equation in terms of an average

{\dot{Q}}_{sink}

or from equating the conservation Eq. and performance Eq. and solving for T_char,out

T_{ch, out} = ε_{si} \frac{C_{min}}{{(\dot{m} c_{p})}_{ch}} T_{si, in} + (1 - ε_{si} \frac{C_{min}}{{(\dot{m} c_{p})}_{ch}}) T_{ch, in}

(32)

The last remaining unknowns are T_sink,out and ${\dot{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 T_sink,out. This is more arbitrary than specifying T_source,out as a desired temperature for the condenser inlet.

In specifying T_sink,out, we have several necessary conditions we have to meet. First, we must ensure that whatever choice of T_sink,out,

{(\dot{m} c_{p})}_{sink}

should remain greater than

{(\dot{m} c_{p})}_{char}

⁠. We can think about our temperature constraints from a conservative, simplified perspective

Δ T_{sink, min} = T_{m} - \max (T_{sink, in})

(33)

For the heat exchanger to function as desired throughout the charging process, the open loop fluid temperature difference, (T_sink,out − T_sink,in), must be less than ΔT at all times. This is guaranteed to be the case if (T_sink,out − T_sink,in) is less than ΔT_sink,min. Furthermore, we can assert that the following inequality should hold:

T_{sink, in} < T_{char, out} < T_{sink, out} < T_{char, in}

(34)

We can use any guess for T_sink,out to set ${\dot{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} = \frac{{\dot{Q}}_{hx}}{{(\dot{m} c_{p})}_{closed} (T_{max} - T_{min})}

(35)

Performance equation(s)

ε_{hx} (ϕ_{open, in} - ϕ_{closed, in}) = \frac{{\dot{Q}}_{hx}}{C_{min} (T_{max} - T_{min})}

(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 T_open,in in order to determine temperatures in the external heat exchanger, including the inlet boundary conditions to the TES device, T_closed,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. T_closed,out can be subsequently nondimensionalized to serve as ϕ_bc = ϕ_cl,out

ϕ_{cl, out} = \frac{C_{min}}{{(\dot{m} c_{p})}_{cl}} ε_{hx} ϕ_{op, in} + (1 - \frac{C_{min}}{{(\dot{m} c_{p})}_{cl}} ε_{hx}) ϕ_{cl, in}

(37)

Equation is applicable to either heat exchanger and serves as a transient boundary condition for the inlet to the TES if T_open,in is a function of time. Alternatively, we could start with a nondimensional form and combine Eqs. and and solve for ε_hx

ε_{hx} = \frac{{(\dot{m} c_{p})}_{closed}}{C_{min}} (\frac{ϕ_{closed, out} - ϕ_{closed, in}}{ϕ_{open, in} - ϕ_{closed, in}})

(38)

By similar heat exchanger analysis, the effectiveness of the TES can be written as

ε_{tes} = \frac{ϕ_{closed, out} - ϕ_{closed, in}}{ϕ_{closed, out} - θ_{m}}

(39)

Here, the problem of specifying C_min is averted by only having one fluid flowing through the TES. Equation 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 can be solved for ϕ_closed,in

ϕ_{closed, in} = (1 - ε_{tes}) ϕ_{closed, out} + ε_{tes} θ_{m}

(40)

This result can be substituted into Eq. and ϕ_closed,out can be found

ϕ_{cl, out} = \frac{ε_{hx} ϕ_{op, in} C_{min} + ({(\dot{m} c_{p})}_{cl} - ε_{hx} C_{min}) ε_{tes} θ_{m}}{ε_{hx} C_{min} - ε_{hx} ε_{tes} C_{min} + ε_{tes} {(\dot{m} c_{p})}_{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 T₀ = T_e = T_w = T_m and initial melt fraction x_e = x_e_,0. The equations will be nondimensionalized with T_min = min(T₀, T_wi) and T_max = max(T₀, T_wi) substituted into the previous definition for θ, ϕ, and St_io. Because of the initial conditions, θ = θ_m = 0 everywhere, and at all times of interest here, modifying the differential equation for ϕ to

\frac{d ϕ}{d \hat{z}} = - N_{tu} ϕ

(42)

Note that before the transient starts, ϕ = θ_m = 0 everywhere. During the first residence time (t ≤ t_res 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

\hat{z} = 0

⁠. Integration of the equation leads to

ϕ = e^{- N_{tu} \hat{z}}

(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^{*} \geq \hat{z}$ ⁠.

For these conditions, the energy balance equation for the PCM matrix reduces to

\frac{\partial x_{e}}{\partial t^{*}} = \pm N_{tu} R_{we} {St}_{io} ϕ

(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. and integrating this equation with respect to time using the initial condition of x_e = x_e_,0 at a given

\hat{z}

location for time

t^{*} = \hat{z}

yields

x_{e} - x_{e, 0} = \pm N_{tu} R_{we} {St}_{io} (t^{*} - \hat{z}) e^{- N_{tu} \hat{z}}

(45)

Thus, for this either process, the solutions are

ϕ = 0; x_{e} = x_{e, 0}

(46)

before the step change front passes

(t^{*} < \hat{z})

and

ϕ = e^{- N_{tu} \hat{z}}; x_{e} = x_{e, 0} \pm N_{tu} R_{we} {St}_{io} (t^{*} - \hat{z}) e^{- N_{tu} \hat{z}}

(47)

after the step change front passes

(t^{*} \geq \hat{z})

⁠. For the melting transient, these solutions are valid until x_e = 1 at the inlet location,

\hat{z} = 0

⁠. Setting x_e = 1 and

\hat{z} = 0

and solving for t* yields the t* value

t_{end}^{*}

at which the validity of the model ends for extraction

t_{end}^{*} = \frac{1 - x_{e, 0}}{N_{tu} R_{we} {St}_{io}}

(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 x_e = 0 at the inlet location,

\hat{z} = 0

⁠. Setting x_e = 0 and

\hat{z} = 0

and solving for t* yields the t* value

t_{end}^{*}

at which the validity of the model ends

t_{end}^{*} = \frac{x_{e, 0}}{N_{tu} R_{we} {St}_{io}}

(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.

Fig. 3

View large Download slide

Comparison of analytical and numerical ϕ during extraction (N_tu = 52.88, R_we = 0.47, and St_io = 0.05)

Despite the error introduced in ϕ while attempting to capture the physics of the step change front, the solution of x_e 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.

Fig. 4

View large Download slide

Comparison of analytical and numerical x_e during extraction (N_tu = 52.88, R_we = 0.47, and St_io = 0.05)

Parametric Studies

As mentioned earlier in the derivation of the governing equations, the values of dimensionless groups t*, N_tu, R_we, and St_io, 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 N_tu.

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

For specified values of the storage design parameters R_we, St_io, and θ_m, computations can be done for different combinations of N_tu and dimensionless termination time, $t_{end}^{*}$ ⁠, 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 $t_{end}^{*}$ and N_tu.

Fig. 5

View large Download slide

Predicted effectiveness variations with $t_{end}^{*}$ and N_tu for R_we = 1, St_io = 0.1, and θ_m = 0.5, and extraction starting from a completely frozen state or charging starting from a completely melted state with constant thermal properties

From this figure, we can draw conclusions on what a reasonable value for N_tu should be in order to achieve a target effectiveness (for a specified R_we and St_io). For high efficiency TES devices, the recommended range of N_tu is between 1.0 and 12.0 for a reasonable operation time $(t_{end}^{*} > 15)$ and conditions (R_we = 1, St_io = 0.1) as device operation below this threshold cannot achieve desired effectiveness targets. As $N_{tu} = U s_{w} L / {\dot{m}}_{w} c_{pw}$ ⁠, the recommendation that N_tu 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 $t_{end}^{*}$ and N_tu. 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 N_tu (N_tu ≤ 1), the only way to achieve high performance is to greatly increase $t_{end}^{*}$ ⁠. 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 N_tu (1.0 < N_tu < 12.0), there is some room to decrease $t_{end}^{*}$ while increasing N_tu and vice versa. N_tu 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 N_tu (N_tu ≥ 12), we notice diminishing returns. Any further increase in N_tu will not make much difference, and the only relevant parameter to adjust is $t_{end}^{*}$ ⁠. In this range of N_tu, once $t_{end}^{*}$ is increased to the point of achieving a desired effectiveness value, any further increase in $t_{end}^{*}$ is superfluous.

Variation of R_we.

For specified values of the storage design parameters N_tu, St_io, and θ_m, computations can be done for different combinations of R_we and dimensionless termination time, $t_{end}^{*}$ ⁠, 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 $t_{end}^{*}$ and R_we.

Fig. 6

View large Download slide

Predicted effectiveness variations with $t_{end}^{*}$ and R_we for N_tu = 10, St_io = 0.1, and θ_m = 0.5, and extraction starting from a completely frozen state or charging starting from a completely melted state with constant thermal properties

From this figure, we can draw conclusions on what a reasonable value for R_we should be in order to achieve a target effectiveness (for a specified N_tu and St_io). For high efficiency TES devices, the recommended range of R_we is between 0.4 and 1.8 for a reasonable operation time (⁠ $t_{end}^{*} > 15$ ⁠) and conditions (N_tu = 10, St_io = 0.1) as device operation below this threshold cannot achieve desired effectiveness targets. As defined, $R_{we} = ρ_{w} c_{pw} A_{c} / {\bar{ρ}}_{e} {\bar{c}}_{pe} ν'$ ⁠, so the recommendation that R_we be somewhere near unity is a reflection that the working fluid and PCM should have a similar energy capacity.

Variation of St_io.

For specified values of the storage design parameters N_tu, R_we, and θ_m, computations can be done for different combinations of St_io and dimensionless termination time, $t_{end}^{*}$ ⁠, 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 $t_{end}^{*}$ and St_io.

Fig. 7

View large Download slide

Predicted effectiveness variations with $t_{end}^{*}$ and St_io for N_tu = 10, R_we = 1, and θ_m = 0.5, and extraction starting from a completely frozen state or charging starting from a completely melted state with constant thermal properties

From this figure, we can draw conclusions on what a reasonable value for St_io should be in order to achieve a target effectiveness (for a specified N_tu and R_we). For high efficiency TES devices, the recommended range of St_io is between 0.06 and 0.20 for a reasonable operation time $(t_{end}^{*} > 15)$ and conditions (N_tu = 10, R_we = 1) as device operation below this threshold cannot achieve desired effectiveness targets. The Stefan number is written as ${St}_{io} = {\bar{c}}_{pe} (T_{max} - T_{min}) / h_{ls}$ ⁠. 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 St_io 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, R_we, 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, St_io, 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 N_tu and $t_{end}^{*}$ ⁠. 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, T_0,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, T_0,char, is close to the final temperature distribution after extraction T_f_,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 T_m.

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⁻⁵ m²/s) and the working fluid (∼10⁻⁷ m²/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 L²/α (∼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, T_0,_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.

Table 1

Thermal energy storage case study operating conditions

Minimum storage effectiveness, ε_tes

0.95

TES nominal capacity, E_tot

1.5

TJ

Extraction time, t_ext

36,000

s

Charging time, t_char

32,400

s

Initial extraction temperature, T_0,ext

30

°C

Initial charging temperature, T_0,char

f(z)

°C

Hot working fluid inlet temp., T_wi,ext

f(t)

°C

Cold working fluid inlet temp., T_wi,char

f(t)

°C

Extraction mass flow rate,

{\dot{m}}_{ext}

2990

kg/s

Charging mass flow rate,

{\dot{m}}_{char}

1975

kg/s

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 N_tu, differ for the extraction and charging processes because the mass flow rate and the overall heat transfer coefficient differ. Others such as R_we and St_io remain consistent. If the dimensionless variables are changed, the simulation will naturally lead to different results.

Table 2

Thermal energy storage case study dimensionless variables

Dimensionless N_tu,ext

36.84

Dimensionless N_tu,char

37.81

Dimensionless R_we

0.47

Dimensionless St_io

0.13

Dimensionless extraction time

t_{end}^{*}

58.99

Dimensionless charging time

t_{char}^{*}

34.66

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 $Δ \hat{z}$ 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 $Δ \hat{z}$ 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.

Fig. 8

$Dimensionless fluid and element temperature (ϕ, θ) and melt fraction (xe) profiles through space and time with case study parameters from Table 2. (a) Extraction process and (b) charging process.$

View large Download slide

Dimensionless fluid and element temperature (ϕ, θ) and melt fraction (x_e) profiles through space and time with case study parameters from Table 2. (a) Extraction process and (b) charging process.

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.

Fig. 9

View large Download slide

Case study cycle consisting of precooling, storage, night cooling, and storage once more

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 T_m, 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.

Fig. 10

View large Download slide

24 h cycling of a 1.5 TJ TES device and subsystem with case study parameters from Tables 1 and 2

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 T_open,out

T_{op, out} = T_{op, in} (1 - \frac{C_{min}}{{(\dot{m} c_{p})}_{op}} ε_{hx}) + T_{cl, in} \frac{C_{min}}{{(\dot{m} c_{p})}_{op}} ε_{hx}

(50)

where T_open,in is known, T_closed,in is calculated from the TES PDEs, ${\dot{m}}_{open}$ is determined from an average ${\dot{Q}}_{hx}$ ⁠, and ${\dot{m}}_{closed}$ is found via iteration for a prescribed effectiveness. If C_min is ${(\dot{m} c_{p})}_{closed}$ ⁠, then the computation for ${\dot{m}}_{open}$ requires knowledge of T_open,out as in Eq. (21). To calculate an appropriate constant flowrate, ${\dot{m}}_{open}$ ⁠, we initially approximate T_open,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 T_m during extraction and lower than T_m 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. R_we 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. St_io 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 St_io between ∼0.06 and 0.20. N_tu corresponds most directly to heat exchanger effectiveness as it encapsulates transport parameters including ${\dot{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 $t_{end}^{*}$ and N_tu. For different N_tu, R_we, and St_io, the time required to complete the prescribed amount of latent heat transfer varies. Once the required values of N_tu, R_we, St_io, and $t_{end}^{*}$ 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

A_c =: cross-sectional area of flow passage
c_pw =: working fluid specific heat
C_min =: minimum heat capacity rate in external HX
${\bar{c}}_{pe}$ =: effective specific heat of differential element
E_cap =: total energy capacity of the PCM in the TES device
E_tot =: total energy stored in the TES device
h =: convective heat transfer coefficient
h_chan =: height of working fluid flow passage
h_ls =: latent heat of fusion of PCM
h_pcm =: height of PCM matrix section
H_e =: enthalpy of differential element
k_w =: working fluid thermal conductivity
${\bar{k}}_{e}$ =: effective thermal conductivity of differential element
L =: length of flow passage
${\dot{m}}_{closed}$ =: working fluid closed loop mass flow rate through TES; may be written as ${\dot{m}}_{ext}$ or ${\dot{m}}_{char}$
${\dot{m}}_{open}$ =: working fluid open loop mass flow rate through external HX; may be written as ${\dot{m}}_{source}$ or ${\dot{m}}_{sink}$
$\dot{Q}$ =: heat transfer rate in the external heat exchanger
s_w =: wetted perimeter of flow passage
t_char =: time elapsed to end of charging process
t_ext =: time elapsed to end of extraction process
T_closed =: TES closed loop working fluid temperature; may be written as T_ext or T_char
T_e =: TES device matrix element temperature
T_m =: melting temperature of PCM
T_max =: maximum temperature in the system
T_min =: minimum temperature in the system
T_open =: external HX open loop working fluid temperature; may be written as T_source or T_sink
T_w =: TES device working fluid temperature
T₀ =: initial storage temperature for extraction or charging
U =: overall heat transfer coefficient between bulk working fluid and PCM matrix element
w_chan =: width of flow passage
x_e =: 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
${\bar{ρ}}_{e}$ =: effective density of TES differential element
ϕ =: dimensionless working fluid temperature

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

3.

Ismail

,

K.

, and

Goncalves

,

M.

,

1999

, “

Thermal Performance of a Pcm Storage Unit

,”

Energy Convers. Manage.

,

40

(

2

), pp.

115

–

138

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

5.

Pielichowska

,

K.

, and

Pielichowski

,

K.

,

2014

, “

Phase Change Materials for Thermal Energy Storage

,”

Prog. Mater. Sci.

,

65

, pp.

67

–

123

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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

.

Google Scholar

Crossref

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.

Multiscale Transient Modeling of Latent Energy Storage for Asynchronous Cooling

Introduction

Analysis Framework

Thermal Energy Storage Energy Balance.

Thermal Energy Storage Governing Equations.

Coupled Heat Exchanger Energy Balance.

Heat Exchanger During Extraction.

Thermal Energy Storage Device Within Subsystem.

Heat Exchanger During Charging.

Subsystem Governing Equations.

Performance at Short Times and Validation Via Closed Form Solution

Parametric Studies

Variation of N_tu.

Variation of R_we.

Variation of St_io.

Key Takeaways From Parametric Studies.

Case Study Example

Extraction/Melting.

Charging/Freezing.

Subsystem Full Cycle.

Conclusions

Acknowledgment

Funding Data

Nomenclature

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Multiscale Transient Modeling of Latent Energy Storage for Asynchronous Cooling

Introduction

Analysis Framework

Thermal Energy Storage Energy Balance.

Thermal Energy Storage Governing Equations.

Coupled Heat Exchanger Energy Balance.

Heat Exchanger During Extraction.

Thermal Energy Storage Device Within Subsystem.

Heat Exchanger During Charging.

Subsystem Governing Equations.

Performance at Short Times and Validation Via Closed Form Solution

Parametric Studies

Variation of Ntu.

Variation of Rwe.

Variation of Stio.

Key Takeaways From Parametric Studies.

Case Study Example

Extraction/Melting.

Charging/Freezing.

Subsystem Full Cycle.

Conclusions

Acknowledgment

Funding Data

Nomenclature

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only

Variation of N_tu.

Variation of R_we.

Variation of St_io.