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.

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\u02d9sink>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 $(m\u02d9closed=0)$, a hot working fluid $(m\u02d9source>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*, thereby melting the PCM and chilling the closed-loop working fluid. This will chill the open-loop working fluid (

_{m}*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.

## 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\u02d9closed$. 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.

*dH*

_{e}*x*

_{e}h_{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

*Us*is the conductance from the working fluid to the PCM matrix and (

_{w}dz*T*−

_{w}*T*) is the driving temperature difference. Taking the time derivative of the right-hand side of Eq. (1) and dividing through by constants $\rho \xafe\nu \u2032dzc\xafpe$, 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

_{e}*T*(sensible heat transfer), or the melt fraction,

_{e}*x*(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:

_{e}*T*≠

_{e}*T*and

_{m}*x*= 0 or

_{e}*x*= 1

_{e}*T*=

_{e}*T*and 0 <

_{m}*x*< 1.

_{e}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.

*T*, should fall somewhere between

_{m}*T*

_{min}and

*T*

_{max}

*θ*≠

*θ*and

_{m}*x*= 0 or

_{e}*x*= 1

_{e}*θ*=

*θ*and 0 <

_{m}*x*< 1. For sensible and latent heat transfer in the storage element where

_{e}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.

*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 $\Delta z\u0302$. 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.

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.

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.

*T*

_{source,out}, our target for the open loop fluid exiting the heat exchanger. If we solve Eq. (20) for

*T*

_{ext,in}, we find that

*C*

_{min}is either $(m\u02d9cp)source$ or $(m\u02d9cp)ext$. We can solve Eqs. (19) and (20) for the unknown temperature,

*T*

_{ext,out}, but must specify whether

*C*

_{min}is associated with the open loop or closed loop fluid

It is hard to predict ahead of time which fluid will be associated with *C*_{min} as $m\u02d9ext$ 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\u02d9ext$ 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.

*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

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.

*C*

_{min}. For the charging process, we can confidently assume that

*C*

_{min}is associated with the closed loop fluid $(m\u02d9cp)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

*T*

_{char,in}such that

*ε*

_{tes}is achieved. In this way, we can define $m\u02d9char$. In order to solve for

*T*

_{char,in}, we must approximate $Q\u02d9sink$. 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\u02d9sink$.

*T*

_{char,out}can be found from the conservation equation in terms of an average $Q\u02d9sink$ or from equating the conservation Eq. (28) and performance Eq. (29) and solving for

*T*

_{char,out}

The last remaining unknowns are *T*_{sink,out} and $m\u02d9sink$, 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.

*T*

_{sink,out}, we have several necessary conditions we have to meet. First, we must ensure that whatever choice of

*T*

_{sink,out}, $(m\u02d9cp)sink$ should remain greater than $(m\u02d9cp)char$. We can think about our temperature constraints from a conservative, simplified perspective

*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:

We can use any guess for *T*_{sink,out} to set $m\u02d9sink$ 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.

*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}

*T*

_{open,in}is a function of time. Alternatively, we could start with a nondimensional form and combine Eqs. (35) and (36) and solve for

*ε*

_{hx}

*C*

_{min}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,out}can be found

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

*T*

_{0}=

*T*=

_{e}*T*=

_{w}*T*and initial melt fraction

_{m}*x*=

_{e}*x*

_{e}_{,0}. The equations will be nondimensionalized with

*T*

_{min}= min(

*T*

_{0},

*T*

_{wi}) and

*T*

_{max}= max(

*T*

_{0},

*T*

_{wi}) substituted into the previous definition for

*θ*,

*ϕ*, and St

_{io}. Because of the initial conditions,

*θ*=

*θ*= 0 everywhere, and at all times of interest here, modifying the differential equation for

_{m}*ϕ*to

*ϕ*=

*θ*= 0 everywhere. During the first residence time (

_{m}*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 $z\u0302=0$. Integration of the equation leads to

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*\u2265z\u0302$.

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.

*ϕ*using Eq. (43) and integrating this equation with respect to time using the initial condition of

*x*=

_{e}*x*

_{e}_{,0}at a given $z\u0302$ location for time $t*=z\u0302$ yields

*x*= 1 at the inlet location, $z\u0302=0$. Setting

_{e}*x*= 1 and $z\u0302=0$ and solving for

_{e}*t** yields the

*t** value $tend*$ at which the validity of the model ends for extraction

*x*= 0 at the inlet location, $z\u0302=0$. Setting

_{e}*x*= 0 and $z\u0302=0$ and solving for

_{e}*t** yields the

*t** value $tend*$ at which the validity of the model ends

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

## 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 $U\xaf$ value, which must be averaged both spatially and temporally as described in a previous paper [13]. With an average $U\xaf$ 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, $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

*N*

_{tu}.

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 $(tend*>15)$ and conditions (*R*_{we} = 1, St_{io} = 0.1) as device operation below this threshold cannot achieve desired effectiveness targets. As $Ntu=UswL/m\u02d9wcpw$, 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 $tend*$ 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 $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 *N*_{tu} (1.0 < *N*_{tu} < 12.0), there is some room to decrease $tend*$ 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 $tend*$. In this range of *N*_{tu}, once $tend*$ is increased to the point of achieving a desired effectiveness value, any further increase in $tend*$ 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, $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

*R*

_{we}.

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 ($tend*>15$) and conditions (*N*_{tu} = 10, St_{io} = 0.1) as device operation below this threshold cannot achieve desired effectiveness targets. As defined, $Rwe=\rho wcpwAc/\rho \xafec\xafpe\nu \u2032$, 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, $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 St

_{io}.

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 $(tend*>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 $Stio=c\xafpe(Tmax\u2212Tmin)/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 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 $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, *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^{−5} m^{2}/s) and the working fluid (∼10^{−7} m^{2}/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*^{2}/*α* (∼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.

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, $m\u02d9ext$ | 2990 | kg/s |

Charging mass flow rate, $m\u02d9char$ | 1975 | kg/s |

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, $m\u02d9ext$ | 2990 | kg/s |

Charging mass flow rate, $m\u02d9char$ | 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.

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 $tend*$ | 58.99 |

Dimensionless charging time $tchar*$ | 34.66 |

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 $tend*$ | 58.99 |

Dimensionless charging time $tchar*$ | 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 $\Delta z\u0302$ 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 $\Delta z\u0302$ 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 *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.

*T*

_{open,out}

where *T*_{open,in} is known, *T*_{closed,in} is calculated from the TES PDEs, $m\u02d9open$ is determined from an average $Q\u02d9hx$, and $m\u02d9closed$ is found via iteration for a prescribed effectiveness. If *C*_{min} is $(m\u02d9cp)closed$, then the computation for $m\u02d9open$ requires knowledge of *T*_{open,out} as in Eq. (21). To calculate an appropriate constant flowrate, $m\u02d9open$, 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*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.

_{m}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 $m\u02d9closed$ 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 *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 $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

*A*=_{c}cross-sectional area of flow passage

*c*_{pw}=working fluid specific heat

*C*_{min}=minimum heat capacity rate in external HX

- $c\xafpe$ =
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

- $k\xafe$ =
effective thermal conductivity of differential element

*L*=length of flow passage

- $m\u02d9closed$ =
working fluid closed loop mass flow rate through TES; may be written as $m\u02d9ext$ or $m\u02d9char$

- $m\u02d9open$ =
working fluid open loop mass flow rate through external HX; may be written as $m\u02d9source$ or $m\u02d9sink$

- $Q\u02d9$ =
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*_{0}=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

- $\nu \u2032$ =
PCM matrix volume per unit flow length

*ρ*_{pcm}=density of PCM

*ρ*=_{w}working fluid density

- $\rho \xafe$ =
effective density of TES differential element

*ϕ*=dimensionless working fluid temperature