A Physics-based Scaling of the Charging Rate in Latent Heat Thermal Energy Storage Devices
Kedar Prashant Shete, S. M. de Bruyn Kops, Dragoljub, Kosanovic
AA Physics-based Scaling of the Charging Rate in Latent HeatThermal Energy Storage Devices
Kedar Prashant Shete a, ∗ , S. M. de Bruyn Kops , Dragoljub (Beka) Kosanovic a Mechanical and Industrial Engineering, 219 Engineering Laboratory, University of Massachusetts, 160Governors Drive, Amherst, MA 01003-2210
Abstract
Thermal energy storage (TES) is increasingly recognized as an essential component of ef-ficient Combined Heat and Power (CHP), Concentrated Solar Power (CSP), Heating Ven-tilation and Air Conditioning (HVAC), and refrigeration as it reduces peak demand whilehelping to manage intermittent availability of energy (e.g., from solar or wind). Latent HeatThermal Energy Storage (LHTES) is a viable option because of its high energy storage den-sity. Parametric analysis of LHTES heat exchangers have been focused on obtaining datawith laminar flow in the phase changing fluid and then fitting a functional form, such asa power law or polynomial, to those data. While this approach can produce an accuratecorrelation applicable within the range of data used for its creation, it does not reveal detailsabout the underlying physics. In this paper we present a parametric framework to analyzeLHTES devices by identifying all relevant fluid parameters and corresponding dimensionlessnumbers. We present 64 simulations of an LHTES device using the finite volume methodat four values of the Grashof, Prandtl and Reynolds numbers in the phase change material(PCM) and heat transfer fluid (HTF). We observe that with sufficient energy available in theHTF, the effects of the HTF Reynolds number and Prandtl number on the heat transfer rateare negligible. Under these conditions, we propose a time scale for the variation of energystored (or melt fraction) of the LHTES device based on the Fourier number(
F o ), Grashofnumber( Gr p ) and Prandtl number( P r p ) and observe a Gr p and P r (1 / p dependency. We alsoidentify two distinct regions in the variation of the melt fraction with time, namely, thelinear and the asymptotic region. The linear region is characterized by constant and highheat transfer rates, making it the relevant region for operating an energy storage device. Wealso predict the critical value of the melt fraction at the transition between the two regions.From these analyses, we draw some conclusions regarding the design procedure for LHTESdevices. ∗ Corresponding author
Email addresses: [email protected], [email protected] (Kedar Prashant Shete ), [email protected] (S. M. de Bruyn Kops), [email protected] (Dragoljub (Beka) Kosanovic)
Preprint submitted to Applied Energy September 23, 2020 a r X i v : . [ phy s i c s . f l u - dyn ] S e p . Introduction Thermal energy storage (TES) is increasingly recognized as an essential component ofefficient combined heat and power (CHP), concentrated solar power (CSP), heating ven-tilation and air conditioning (HVAC), and refrigeration as it reduces peak demand whilehelping to manage intermittent availability of energy (e.g., from solar or wind). As discussedin more detail below, it has the potential to reduce energy consumption and reduce pollutiongeneration by making existing technologies more efficient and by enabling the integration ofrenewable energy sources with minimum energy curtailment.Given the thermo-physical properties of a heat storage material, it is straightforward tocompute the amount of that material required to store a given amount of heat. The chal-lenge is in designing a physical device that enables sufficiently high heat transfer rate for apractical system. If, for example, a TES is to be coupled with a CHP plant, the TES mustbe able to store and release heat at the time scale of the transients in the CHP system.Designing a TES system to meet this requirement is difficult because a very large numberof parameters affect the heat transfer rate including the properties of the working fluids, thefluid dynamical regimes of those fluids when the system is operating, the geometry of theheat exchanger and storage device, and the operating conditions for the entire system. Inthis paper, we present an approach for dealing with this complexity that consists of system-atically defining the relevant dimensionless parameters and then writing the relationshipsbetween these parameters based on physical understanding derived from theory and fromthe literature about TES systems. Of course this approach is not unique to this paper, butwe apply it the specific case of latent heat thermal energy storage (LHTES) to demonstratehow the approach can introduce physical understanding into relationships between param-eters that have typically been studied empirically and, thereby, simplify the overall designprocess.In the remainder of this section we review the motivation for studying TES and, inparticular, applications that significantly benefit from LHTES. We then review some of thefundamental studies in LHTES that provide the physical understanding necessary for ourapproach. In § § §
4. We demonstrate our methodology in § § § Applications that are being improved significantly with thermal energy storage includeconcentrating solar power (CSP) plants; Denholm et al. [8] report round-trip efficiencies closeto 100% when energy from CSP’s is stored as thermal energy rather than electrical energy.They also report that “cold storage” enables extremely high efficiency of cooling systems byshifting demand to off-peak hours. Nithyanandam and Pitchumani [21], based on their studyof charging and discharging cycles of a LHTES heat exchanger, emphasize the importance of2HTES for the effective functioning of CSP. Performance of cogeneration power plants alsoimproves when they are combined with thermal energy storage [12, 19]. Venkitaraj et al. [25]investigate experimentally the use of nano-particle enhanced LHTES for waste heat recoveryfrom IC engines and observe up to 18% increase in the energy savings.In addition to improvements in energy efficiency, thermal energy storage can reduce emis-sion of pollutants. For example, Li et al. [17] calculate the effect of a LHTES system usedto recover waste heat from a heavy duty diesel engine and conclude a potential 40% im-provement in engine warm up time during which the engine produces suboptimal emissions.Arbabzadeh et al. [2] report the huge potential impact of energy storage on decarbonizationof electricity production by allowing electricity usage for heating and cooling to be synchro-nized with when renewable energy is available. Specifically, they conclude that, for the stateof California, thermal energy storage can result in an 18% reduction in carbon dioxide emis-sion and a 21% reduction in renewable energy curtailment, that is, the reduction of outputof a renewable resource below what it could have otherwise produced.
Thermal energy storage can be classified into three major types: sensible heat storage,latent heat storage and thermochemical energy storage. For the applications discussed in thepreceding paragraphs, a desired characteristics of TES include: • High volumetric energystorage density • Heat recovery at constant temperature • Low cost • Fast heat transferrate. LHTES has inherent advantages over other TES systems with respect to high storagedensity and heat transfer at constant temperature. High storage density, in turn, tends tolead to lower cost. Thus, LHTES would appear to be a very attractive option for improvingthe energy efficiency and reducing emmisions of a variety of types of power plants andengines. Indeed, Mongibello et al. [20] study two different types of thermal energy storagefor residential micro-CHP systems and conclude that LHTES is preferred over sensible energystorage (such as hot water) in terms of cost and size. They also conclude that further analysisshould be made, including of the long-term performance and degradation of these systemsover time, in order to assess the convenience of using them for thermal energy storage. Joharet al. [15] implement a LHTES system within a micro-CHP plant and shows LHTES canbe a viable option. They note, though, that improved design procedures and performancemodeling of phase change heat exchangers are essential.The last characteristic in the list above, fast heat transfer rate, is the motivation forthe research reported in this paper. Heat transfer rate is determined primarily by the fluiddynamics and geometry of the heat exchanger rather than specifically by the storage mech-anism, with turbulent flow over large surface areas leading to high heat transfer rates. Asreviewed in § Given the latent heat of fusion of a phase changing material (PCM), it is relativelysimple to calculate the amount of PCM that a LHTES system needs in order to store aspecific amount of energy. The challenge is in designing a system with the required heat ransfer rate , which as evident from previous studies depends on a number of geometrical,material and operating parameters. Given the complexity of the problem, it is common forindividual research studies to focus on a subset of the parameters affecting heat transfer rate.An important first step is to begin with simplified governing equations for heat transfer, forexample, neglecting convective heat transfer in the PCM [3, 6, 7, 24, 27]. Natural convection,however, is a key component of accurately modeling energy storage rates [3]. Heat exchange geometry is a crucial factor affecting the heat transfer rate of LHTES.Geometry parameters that have been studied include the inner and outer diameters in anannular geometry with PCM in the annulus and HTF in the inner pipe [6, 16], HTF pipe wallthickness [7] and diameter of the HTF pipe [27]. Adding fins in the PCM has been shown toimprove charging rates, stored energy and melting front depth [5, 16, 27]. Bhagat et al. [5]conduct an optimization study of fin height, fin thickness and number of fins using ANSYSFluent and laboratory scale heat exchanger data and conclude that for a given percentageof fin material/metal inside the heat exchanger, a higher number of thinner fins lead tobetter heat transfer. The overall configuration of the LHTES is also an important factor,and various configurations including a single HTF pipe inside an annular PCM container,multiple HTF tubes inside a PCM pipe, PCM modules floating inside an HTF container anddirect contact between HTF and PCM have been studied [11]. The orientation of the devicealso affects its performance and has been studied by Kalapala and Devanuri [16].
The thermophysical properties of the HTF and PCM such as thermal conductivity andspecific heat capacity are also important parameters affecting the performance of LHTES[6, 7, 10, 11, 27]. Gasia et al. [11] conclude that an increase in specific heat capacity ofHTF of 4.9 times and in thermal conductivity of HTF of 3 times improves the chargingtimes by 44 %. Farid et al. [10] note the importance of material properties by observingthat materials such as paraffins have moderate energy storage density and low cost, butalso have low thermal conductivity, which affects their utility as energy storage materials.Hydrated salts, on the other hand have larger thermal conductivity and large energy storagecapacity, but their use is affected by other material properties like supercooling and phasesegregation. They conclude that the melting point is the most important characteristic inselecting a phase change material and point out the importance of creating materials thathave an adjustable melting point.
System operating parameters such HTF mass flow rate and temperature have a dominanteffect on the LHTES performance because it is the HTF that determines the maximum rateat which energy can be exchanged with the PCM. The effect of HTF mass flow rate andHTF temperature has been studied by a number of authors [3, 16, 24]. For example, a studyconducted in terms of dimensionless parameters is that of Teamah et al. [24], which was aparametric numerical finite difference analysis of total heat transfer gain in an cylindrical4ank with encapsulated PCM. The parameters they studied are the HTF Reynolds number inthe range 20 < Re < . < Ste < .
4, 0 . < ( ρC p ) ∗ < . ρC p ) ∗ is the ratio of effective thermal capacity (Density × Specific Heat Capacity) ofthe PCM to HTF, 0 . < θ m < . θ m is the ratio of the difference between the PCMmelting and HTF inlet temperature to the difference between the HTF inlet temperature andthe starting temperature of the system, and the Fourier Number F o which is non-dimensionalcharging time. They obtained a dependency of
F oRe . θ m and Ste ( ρC p ) ∗ for the total energygain and concluded that the dependency of Re . originates from the turbulent convectioncoefficient correlation used within their finite difference calculation. Understanding the effectof individual parameters on the performance and quantifying their importance relative toother parameters will greatly support the design process for LHTES devices [10]. It isadvantageous to have dimensionless results instead of purely experimental data pertainingto just one device [3]. From the foregoing review, it is apparent that the foundation has been laid for under-standing the individual factors affecting heat transfer rate in LHTES systems. Less progresshas been made on combining these individual factors to form a complete set of relevant di-mensionless parameters suitable for robust modeling and design guidance for creating LHTESsystems having sufficiently high heat transfer rates for commercial applications. We beginour study in § § §
4, wediscuss the equations used and validation of our simulations. In § § § § §
6, we examine the underlying convection physics and proposescaling laws for heat transfer rate as a function of the Reynolds number in HTF, Grashofnumber in PCM and Prandtl numbers in both the HTF and PCM. Some conclusions aboutthe scaling obtained and the reason causing these regions are presented in §
2. Parametrizing the problem
A variety of configurations exist for LHTES systems, but they have certain commonelements. Typical LHTES devices consist of a heat exchanger with a heat transfer fluid(HTF), such as oil, pumped across one side of a solid interface and a PCM driven by naturalconvection on the other side. Starting from the solid state in the PCM, introduction of heatto the system via the HTF melts some of the PCM and buoyancy begins to drive flow. Threefactors quantify the practical performace of an energy storage/LHTES device: the chargingrate, the discharging rate and the storage capacity. In a LHTES device, the storage capacity5 omenclature
Symbol Description Subscript Description η Melt fraction f Heat transfer fluid (HTF) ρ Density p Phase change material (PCM) µ Viscosity t HTF tube β Volumetric expansioncoefficient c PCM container T Temperature i Inner u Velocity o Outer M Mass in Inlet α Thermal diffusivity mean
Mean ν Kinematic viscosity f r
HTF Reynolds number q Heat transfer rate into controlvolume f p
HTF Prandtl number h Mean heat transfer coefficient pg PCM Grashof number D Diameter pp PCM Prandtl number A mush Mushy zone constant τ Dimensionless time λ Liquid fraction - - τ Generic dimensionless time - -is very simple to predict because it is directly proportional to the mass of the of PCM inthe system. The charging and discharging rates are more difficult to predict because, asreviewed in the previous section, they depend on the geometry of the heat transfer surface,the thermophysical propoerties of the fluids, and the operating conditions of the entiresystem. Here we consider only the charging rate because, while the discharging rate may bedifferent, the approach to parameterizing the modeling both rates is the same.A common approach to modeling the charging rate is to fit an assumed function toexperimental or numerical data. For example, Rathod and Jyotirmay [22] use polynomialregression to describe the melting time as a function of the Reynolds number in the HTF,the Stefan number of the PCM and the ratio of initial temperature of the PCM and inlettemperature of the HTF. Diarce et al. [9] assume a product of power-law relationship to fit theFourier number as a function of the Biot number, the Stefan number and two dimensionlesstemperature constants. This approach can produce effective correlations over the range ofdata used to produce them but offer limited physical insight to enable predicting the heattransfer rate outside the range that was measured.
In commercial applications, the geometry of the heat exchanger is such that turbulentflow of the PCM can be expected unless the melted fraction is extremely small. Turbulentflow studies in LHTES systems are limited by the practical size of laboratory experiments6nd current limitations in computing capability. Therefore, we procede using a physics-based approach to hypothesize the correct functional forms for the relationships betweendimensionless flow parameters. This approach begins with identifying the dimensional sys-tem parameters expected to be important for the performance of LHTES systems. Theseare tabulated in Table 1.We note that the number of parameters affecting LHTES system performance is extremelylarge. For example, the macroscopic geometry of the device can be quite complicated and,e.g., microsopic geometry of the heat transfer surfaces is a topic unto itself. Here we haveassumed an annular geometry of smooth-walled tubes with the HTF in the inner tube andPCM in the annulus.
Given the very large number of parameters in Table 1 and are narrowing of the focusof this paper to a simple geometry, we procede using only the 16 parameters in the tablethat are indicated by boxes. To further simplify the problem, we define the mean surfacetemperatures as the average temperature over that surface. For example, T mci is the averagetemperature over the inner surface of the inner boundary of the PCM container (diameter D ci ). Next we apply the Buckingham Pi theorem to determine the minimum number ofdimensionless groups given the dimensional parameters in Table 1 and the assumption thatmass, length, time, and temperature are independent dimensions. This leads us to expect12 dimensionless parameters. Given the expected number of groups and the well-establisheddefinitions of many of them, we arrive at the dimensionless groups in Table 3. In the foregoing analysis we have sought the minimum number of dimensionless groupswhile recognizing that there are multiple ways to define these groups. In particular, it isuseful to consider that the rate of change of the melt fraction is related by conservation ofenergy to the heat transfer rate with the assumption of isothermal heat transfer. The HTFtransports energy into the system, which is then transferred to other system components.Let the total heat transfer rate to the system be denoted by q f , as given in Table 2, whichis equal to the heat transfer out of the control volume HTF. This heat is then distributedbetween the PCM( q p ), HTF tube( q t ) and PCM container( q c ). There will be a transient asthe temperatures of the PCM and heat transfer surfaces adjust to the melting point of thePCM. Once this transient is finished, q t and q c are expected to be small compared to q p dueto the high volume and the high heat capacity of the PCM. q p can be further split into twocomponents, the sensible heating rate q ps which causes temperature rise in the PCM andthe latent heating rate q pl which causes melting of the PCM. q ps is typically much smallerthan q pl due to reasons similar as above; the latent heat capacity of the PCM L is a coupleof orders of magnitude higher than the sensible heat capacity Cp p . Assuming that the PCMcontainer is well insulated, the heat transfer rate balance can be written as q f = q p + q t + q c = q pl + q ps + q t + q c (1)7roperties HTF PCM HTF Tube PCM ContainerDensity ρ f ρ p ρ t ρ c Specific Heat Capacity Cp f Cp p Cp t Cp c Viscosity µ f µ p - -Thermal Conductivity k f k p k t k c Volumetric ExpansionCoefficient β f β p β t β c Inlet Temperature T in - - -Initial Temperature T i Freezing Temperature - T solidus - -Melting Temperature - T liquidus - -Latent Heat Capacity - L - -Time t Time for Solidificationwith Under-cooling - ∆ t s - -Average InletVelocity/Average Velocity u f - - -Length - - l t l c Initial Mass M f M p M t M c Diameters - - D t D ci , D co Container to FluidInterface Area A f A p - -Derived parametersMean Surface Temperature( (cid:82) A T dA/A ) - - - T mci Mean Heat TransferCoefficient ( h ) h f h p - -Mean Melting Temperature(( T solidus + T liquidus ) /
2) - T mean - - Table 1: Dimensional parameters that affect heat transfer rate of LHTES devices. Boxed parameters havebeen used for the Buckingham Pi analysis. k/ρC p ) α f α p - -Kinematic Viscosity ( µ/ρ ) ν f ν p - -Mean volume temperature( (cid:82) V T dV /V ) T f T p - -Heat transfer rate outof/into control volume q f q p q t q c Table 2: Other derived parameters
Out of these, q pl is of particular interest as it represents high quality energy available at afixed temperature. The integral of q pl from the onset of melting to the current time is relatedto the melted fraction of PCM ( η ) by (2), where η is defined as the mass of melted PCM tothe total mass of PCM. q pl = M p L dηdt (2) q pl or η can be written as a function of all the parameters in Table 1. Depending on howmany of those we vary for our simulations, we get a corresponding number of dimensionlessnumbers. This is discussed further in §
3. Numerical simulations approach
To investigate the relationships between the dimensionless groups described in Table 3,we seek benchmark simulations that are free from models. In practice, some modeling isinherent in simulations, starting with the continuum approximation, which omits moleculareffects inherent in the phase change process. Our approach is to limit the modeling in thesimulations to: 1. The HTF is incompressible and Newtonian. 2. The initial temperatureof the entire unit is uniform and the PCM is in the solid phase 3. The thermophysicalproperties of the liquid HTF, the PCM and the container are constant except for the densityof the PCM. 4. The density changes in the PCM and their scaling height are small so thatthe non-hydrostatic Boussinesq approximation is applicable. 5. The kinetic and thermalenergies of the PCM are decoupled. 6. The equations of motion for the liquid and solidphases of the PCM are coupled using the approach of Voller and Prakash [26]. To make thesimulations more tractable, only laminar flow of the HTF and PCM are considered so thatthe axisymmetric equations of motion are applicable.The PCM flow is assumed to satisfy the non-hydrostatic Boussinesq assumptions forconservation of mass and momentum which can be written in cylindrical coordinates as1 r ∂ ( ru r ) ∂r + ∂ ( u z ) ∂z = 0 (3a)9 a b l e : N o n - D i m e n s i o n a l g r o up s a ff ec t i n g h e a tt r a n s f e rr a t e s i n L H T E Sd e v i ce s . T h e i nd e p e nd e n t g r o up s c o rr e s p o nd i n g t o p a r a m e t e r s i n h a v e b ee nb o x e d . O t h e r d e p e nd e n t g r o up s li k e t h e R a y l e i g hnu m b e r h a v e b ee n m e n t i o n e ddu e t o t h e i r i m p o r t a n ce i n li t e r a t u r e . N u m b e r s H T F P C M H T F T ub e P C M C o n t a i n e r R e y n o l d s nu m b e r ( R e ) ρ f D t u f / µ f --- F o u r i e r nu m b e r ( F o ) - α p t / ( D c o − D c i ) -- P r a nd t l nu m b e r ( P r ) C p f µ f / k f C p p µ p / k p -- P ´ ec l e t nu m b e r ( P e ) D t u f / α f --- G r a s h o f nu m b e r ( G r ) - g β p ( D c o − D c i ) ( T m c i − T m e a n ) / ν p -- R a y l e i g hnu m b e r ( R a ) - g β p ( D c o − D c i ) ( T m c i − T m e a n ) / ν p α p -- A s p ec tr a t i o ( A R ) -- l t / D t l c / ( D c o − D c i ) S t e f a nnu m b e r ( S t e ) - C p p ( T m c i − T m e a n ) / L -- B i o t nu m b e r( B i ) -- h f ( D c i − D t ) / k t - N u ss e l t nu m b e r ( N u ) h f D t / k f h p D c i / k p -- M e l t i n g t o H e a t T r a n s f e r T i m e s c a l e - ∆ t s / t (cid:117) -- (cid:18) ∂u r ∂t + u r ∂u r ∂r + u z ∂u r ∂z (cid:19) = − ∂p ∗ ∂r + 2 µ ∂ u r ∂r + µ ∂∂z (cid:18) ∂u r ∂z + ∂u z ∂r (cid:19) + 2 µr (cid:18) ∂u r ∂r − u r r (cid:19) + S r (3b) ρ (cid:18) ∂u z ∂t + u r ∂u z ∂r + u z ∂u z ∂z (cid:19) = − ∂p ∗ ∂z + 2 µ ∂ u z ∂z + µ ∂∂r (cid:18) ∂u z ∂r + ∂u r ∂z (cid:19) + µr (cid:18) ∂u r ∂z + ∂u z ∂r (cid:19) + S z + S b (3c)Here, the force terms S r , S z are the momentum sinks used by the melting/solidificationmodel of Voller and Prakash [26] and are given as S r = A mush (1 − λ ) ( λ + (cid:15) ) u r , S z = A mush (1 − λ ) ( λ + (cid:15) ) u z (4)where A mush is the mushy zone constant of Voller and Prakash [26]. In their model, theliquid fraction λ is calculated as λ = T liquidus − TT liquidus − T solidus . (5)The source term S b is the buoyancy force given by ρ β ( T − T ) g , where T and ρ are thereference temperature and reference density used for the Boussinesq approximation and β isthe coefficient of thermal expansion. Within the Boussinesq approximation, viscous heatingof the fluid is taken to be negligible and so the thermal and mechanical energy equationsdecouple. The mechanical energy equation can be derived by taking the dot product ofvelocity and momentum. The thermal energy equation can be written in terms of enthalpyor temperature. ANSYS Fluent, which is the code used for simulations, uses the enthalpyform of the equation, given as ρ (cid:18) ∂h∂t + u r ∂h∂r + u z ∂h∂z (cid:19) = ρ ˙ q g + 1 r ∂∂r (cid:18) kr ∂T∂r (cid:19) + ∂∂z (cid:18) k ∂T∂z (cid:19) + S e (6)where S e = ∂ρ ∆ H∂t + ∇ · ( ρ(cid:126)u ∆ H ) . (7)Here ∆ H = Lλ (8)is the latent heat enthalpy change for a material volume of PCM.The work associated with the momentum terms S r + S z is S (cid:48) e = (cid:126)u · ( S r ˆ r + S z ˆ z ) (9)Due to the small velocities, S (cid:48) e (cid:28) S e and is neglected.11 . Simulations The simulation geometry is chosen for validation against the laboratory results of Longeonet al. [18]. The physical configuration is shown in Figure 1a. Due to the fact that thecylinder is oriented vertically and that the flow regime is laminar, it can be assumed thatthe flow is axially symmetric and the equations of motion in § (a) Experimental setup by Longeon et al.[18] (b)
Our computational domain and grid.Figure 1: Simulation geometry based on the laboratory experiments of Longeon et al. [18]. Panel (a) showstheir experimental setup, reprinted with their permission. the computational geometry are in Table 4a. The structured grid used for the simulations isshown as Figure 1b. The properties of the PCM, given in Table 4c, are matched to those inLongeon et al. [18], with the exception of the sensible specific heat capacity and the density,which are different between the solid to liquid in the experiments but in the simulations areset to average values shown in Table 4b. The properties of stainless steel in the simulationsare density 8030 kg/m , specific heat capacity 502 . J/kgK and thermal conductivity of16 . W/mK . The simulations are conducted using the finite volume code ANSYS Fluent. The sim-ulation parameters are in Tables 4b and 4d. In the phase-change model (4), the constant12 a) Geometry of simulation domain forvalidation
Parameter Value UnitHTF tubeOuter Radius 10 mmInner Radius 7.5 mmLength 400 mmPCM containerInner Radius 22 mmLength 400 mm (b) Properties of PCM used for simulation
Property Value Unit ρ kg/m L kJ/kgC p kJ/kg.Kµ kg/m.sβ /KT solidus ◦ CT liquidus ◦ Ck W/m.K . (c) Properties of PCM RT35 Rubitherm asreported by [18] Property Value Unit ρ kg/m L kJ/kgC p kJ/kg.Kµ kg/m.sβ /KT m ◦ Ck W/m.K (d) Properties of HTF Water
Property Value Unit ρ kg/m C p kJ/kg.Kµ kg/m.sk W/m.K
Table 4: Simulation parameters used in Longeon et al. [18] A mush defining the mushy zone is taken to be 100,000 and it is observed that the solution isnot strongly dependent on this value.The HTF inlet boundary condition is defined to be a uniform velocity of 0 . m/s andwith static temperature 53 ◦ C. The HTF outlet boundary condition is constant gauge pressureof 0
P a . The internal walls of the tube and PCM container are conjugate heat transferinternal boundaries with no slip. Heat transfer between the outer walls of PCM containerand the room is ignored due to the low temperature differences between the heat transfermedium and ambient conditions. Thus, the outer walls are adiabatic with no slip. Thevertical axis of the heat transfer fluid tube is defined to be a symmetry boundary conditionso that the simulations are axisymmetric.The mass and momentum equations are solved using the pressure based solver with theSIMPLE algorithm used for the pressure velocity coupling. Pressure is discretized using thePRESTO scheme [23]. The momentum and energy equations are discretized using secondorder upwind schemes. The evolution in time is first order implicit, as it is sufficient for mostproblems [1]. The solution is initialized with zero velocity in all directions and an ambienttemperature of 23 ◦ C. The highest velocity in the domain is expected in the HTF and is twicethe mean velocity, for a fully developed flow, which is 0 . m/s .13 Time Step (s)Grid Size (nodes, cells) 0.1 0.013446, 3170 2% 1%4696, 4356 1% 1%27996, 26733 1% 0
Table 5: Grid sensitivity for melt fraction η , difference relative to nodes= 27996 and ∆ t = 0 . Sensitivity of the solutions to grid resolution and time step size are examined by varyingthe time step size by two orders of magnitude and the number of finite volumes in the gridby a factor of approximately eight. The important variable for energy storage is the meltedfraction of the PCM, given as η . We perform simulations with grid sizes 3446, 4696 and27996 and time-step sizes of 0 . .
01 seconds. Note that our smallest grid size andlargest time-step are the same order of magnitude as the grid size of 9000 and time-step of0 . t ) sizes. The results show that with increasing spatialand temporal resolution, the curves approach the results for the finest grid of 27996 and thefinest time-step size of 0 .
01 seconds, and the difference between the intermediate resolutionof 4696 and 0 .
01 seconds and the finest resolution is negligible. Table 5 shows the maximumof the absolute error(as percentage) in η in reference to the finest resolution case, and we seethat the error reduces by less than 1% beyond the intermediate resolution of 4696 and 0 . t = 0 .
01 seconds is sufficientto obtain grid insensitive results.Figure 3 shows the comparison of temperature at a specific location D obtained fromsimulations and measured experimentally by Longeon et al. [18]. We see good agreementbetween the shapes of the experimental and numerical data curves.
Given the results of the validation experiments, we conclude that the simulation techniqueis adequate. After validating our simulation procedure, we proceeded to our parametricstudy. We vary the parameters u f , T in , k p and k f for the geometry and setup discussed inLongeon et al. [18], which we use for validating our simulations described in § u f is flowvelocity for HTF and is the easiest to change through the use of a pump. T in is the inlettemperature of the HTF, which depends on the system from which we are extracting energy. k p and k f are dependent on material properties and additive enhancements and we havemoderate control over them. Table 6 shows the endpoint values of parameters changed andtable 7 shows the corresponding nominal dimensionless numbers.For example, for the four parameter case described here, q pl = q pl ( u f , k f , k p , T in ) (10)14 igure 2: Grid Sensitivity, Average Melt Fraction in PCMFigure 3: Comparison of measured temperature at Point D[18] with simulation for validating simulationprocedure f k f k p T in Re f P r f P r p Gr p Table 6: Physical parameters corresponding to values of dimensionless numbers in table 7. There are atotal of 64 simulations created by varying each number in table 7 independently. For reasons of space, onlythe endpoint cases (numbering 16) have been shown here. Re f P r f P r p Gr p Re f P r f P r p Gr p Re f P r f P r p Gr p P r p Gr p P r p Gr p Table 7: List of dimensionless numbers in parameter space and their values under study. Since theparameter space is four-dimensional, the total number of simulations are 64. q pl . Itsdimensionless equivalent is η , the melt fraction, which we shall use henceforth for presentingresults.
5. The structure of the melt fraction curve η ( t ) In this section, we discuss some observations about the structure of η prior to looking atthe scaled results in section 6. We observed two distinct regions for η , specifically a linear andan asymptotic region, and name the melt fraction at the transition between these regions asthe transition melt fraction denoted by η critical . The time at which η critical occurs is denotedby τ critical ; see section 6 for further discussion about obtaining a dimensionless time τ from t .We elaborate in § § Based on Newton’s law of cooling, the heat transfer rate q f is set by the heat transfercoefficient h f and the temperature difference T in − T mci . Between the onset of melting andthe transition melt fraction η critical , the temperature T mci is roughly constant due providedthat transport of heat by the HTF does not limit the heat transfer into the PCM. If thevariation is not significant compared to the total temperature difference T in − T mci , T mci andthe difference can be considered to be constant. Indeed, for constant wall temperature withinternal laminar flow, the non-dimensional heat transfer coefficient, which is the Nusseltnumber, is constant with value equal to 3.66 [13, eq. 8.55]. Thus, q f is expected to be aconstant in this temperature region.In the PCM, q p , which is comparable to q f , is proportional to the temperature difference T mci − T mean , which is also constant. After an initial transient, q pl is the major componentof q p . The melt fraction η , which is proportional to integral of q pl as shown in (2), isexpected be linear with time. At η critical , the quantity of solid PCM gets small such thatthe characteristic temperature difference in the PCM is T mci − T p , where T p is the meantemperature of the PCM and is approximately equal to the far field temperature. T p isrising inverse-exponentially, which results in the asymptotic behavior of η , as explained inthe following section. Since the purpose of this section is to analze behavior rather than predicting data fromfirst principles, we shall use simplified notation to obtain uncluttered equations. Termsexpected to be constant have been grouped into numbered constants for brevity. Let themass of solid PCM at the time η reaches η critical be m pcm , and let its surface area be a pcm . Letthe heat transfer coefficient on the solid liquid interface be h pcm . Figure 4 shows a cartoonrepresentation of the variables of interest.After η reaches a critical fraction η critical , the mass of solid PCM is small and thecharacteristic temperature difference is closer to T mci − T p rather than T mci − T mean . Theconfiguraton is shown in cartoon form in Figure 4 along with the notation used in thefollowing discussion In this regime, there is limited contact area between the liquid and the17 igure 4: Cartoon figure showing the variables for the asymptotic model. The interface between m pcm andthe liquid is an arbitrarily drawn curve. solid and so most of the heat transferred from the HTF to the PCM raises the temperatureof the PCM. With this approximation, q p = M p Cp p dT p dt = h p A p ( T mci − T p ) ⇒ d ( T mci − T p ) dt = − h p A p M p Cp p ( T mci − T p ) ⇒ d ( T mci − T p ) T mci − T p = − c (cid:48) dt ⇒ ln ( T mci − T p ) = − c (cid:48) t + c ⇒ T p = T mci − c e ( − c (cid:48) t ) (11)In short, the PCM acts as a lumped capacitance Incropera et al. [13, eq. 5.8a] becausethe mass of the solid PCM is insufficient to affect T p . Since the purpose of the foregoing isto arrive at the expected functional form rather than a numerically exact model, we havecombined terms that are approximately constant into the coefficients the coefficients c , c , c (cid:48) . The prime notation denotes constants that carry forward into the final expression givenin (14).Assuming that the remaining solid PCM is at melting temperature and there is no sig-nificant sensible heating of the residual solid, the heat transfer to the solid PCM is q pl = L dm pcm dt = h pcm A pcm ( T p − T mean ) ⇒ d (1 − η ) dt = h pcm A pcm LM p ( T p − T mean ) . (12)The area A pcm depends on the mass of solid PCM and can be calculated if the the shapeof m pcm and its density is known. Since the PCM is close to the melting temperature, the18ensity can be considered to be a constant. At constant density, if m pcm is a sphere, A pcm isproportional to m / pcm . In general, A pcm is proportional to m γpcm where γ is some real numberless than 1, expected to be constant if the melting front geometry and density do not changein the duration of the exponential melting region. Substituting this and (11) into (12) yields d (1 − η ) dt = c h pcm m ppcm LM p ( T p − T mean ) ⇒ d (1 − η ) dt = c h pcm (1 − η ) γ LM − γp ( T p − T mean ) ⇒ d (1 − η )(1 − η ) γ = c h pcm LM − γp ( T p − T mean ) dt = c (cid:16) ( T mci − T mean ) − c e ( − c (cid:48) t ) (cid:17) dt = c ( T mci − T mean ) dt − c c e ( − c (cid:48) t ) dt ⇒ (1 − η ) γ +1 γ + 1 = c ( T mci − T mean ) t + c e ( − c (cid:48) t ) + c ⇒ − η = (cid:16) c (cid:48) + c (cid:48) ( T mci − T mean ) t + c (cid:48) e ( − c (cid:48) t ) (cid:17) γ +1 (13)From this we conclude that the function form for η ( t ) is η = 1 − (cid:16) c (cid:48) + c (cid:48) ( T mci − T mean ) t + c (cid:48) e ( − c (cid:48) t ) (cid:17) γ +1 . (14)Again, we are interested in the form of the equation and the numbered coefficients collectconstant terms that would make the form more difficult to read if included in full.Figures 5 and 6 show the melt fraction and temperature contours when melting hasreached η critical , for cases with different Gr p and P r p values. The cases have vastly differentoperating parameters, but we can see that there are similarities in the melt fraction profiles,for example, the shape of the remaining PCM, which has been identified in § η critical might beuniversal for a given device, at least in the range of dimensionless numbers studied.From the arguments in § § η to vary linearly in time when η issmall and to vary asymptotically with time when η is large with the transition between thetwo regimes being η critical . This conclusion is based entirely on physical reasoning. In thefollowing section, we apply this physical reasoning to the simulation data to understand thetime evolution of η and, in particular, how this time scales with the dimensional quantitiesin Table 1, and to determine the empirical value of η critical , which is of practical importance.
6. Application of methodology to understand flow physics
In the previous sections, we develop our approach for defining the dimensionless groupsimportant for describing a simple LHTES systems such that they can be related to the19 igure 5: Plot of liquid fraction at time when η critical is reached for two Grashof numbers at Re f P r f P r p
2. The corresponding temperature plot is shown to the right. Notice the similarity between theshape of the solid PCM for the two different cases. The figures are best viewed in conjunction with thegeometry and grid shown in figure 1b.Figure 6: Plot of liquid fraction at time when η critical is reached for two Grashof numbers at Re f P r f P r p
4. The corresponding temperature plot is shown to the right. Notice the similarity between theshape of the solid PCM for the two different cases and the cases from figure 5 η as a function oftime in the multi-dimensional space defined Gr p , P r p , Re f and P r f . The melt fraction η hasbeen defined previously as a function of dimensional time. For consistency, we denote theequivalent of η that accepts a generic dimensionless time τ as an argument, by η τ . The goalis to find η τ = η τ ( τ ) with the τ defined in terms of a physically relevant time scale. If we aresuccessful in this then the curves for η τ ( τ ) corresponding to different values of the parameterbeing varied will collapse to a single curve. In this section, we attempt to define a suitable τ based on the flow physics and observations from literature. As described in 2, we use theBuckingham-Pi theorem to organize our approach. A general dimensionless timescale canbe defined as follows - τ = τ ( F o p , Gr p , P r p , Re f , P r f ) (15)where F o p is the PCM Fourier number defined in table 3.The simulation data base consists of 64 cases with parameters tabulated in Table 7. Thesimulations span a four-dimensional parameter space defined by Re f , P r f , P r p and Gr p witha high and low value for the first two and 4 values each for the rest. As mentioned in § Re f and P r f , as long as there is sufficient heat flowinginto the HTF domain. To confirm this, we perform additional simulations with two moreintermediate values of Re f = 796 and Re f = 1443 at P r f = 42, Gr p = 8 . × and P r p = 57, the results of which are shown in figure 12. Similarly, we conduct simulations oftwo intermediate values of P r f = 8 and P r f = 6 for the case with Re f = 2090, P r p = 57and Gr p = 8 . × , results of which are shown in figure 13. Given the expected anddemonstrated weak effect of Re f and P r f , we fix their values and apply the methodology from §
5, beginning with a cut through the parameter space along the plane defined by Re f = 2090and P r f = 5, that is, a particular set of HTF parameters. Based on the reasoning in § §
1, we expect time to scalewith τ | Re f ,P r f = constant = P r / p Gr p F o p (16)for fixed Re f and P r f indicated by the subscripts Re f , P r f = constant . For brevity, we alsointroduce a shorthand notation for a one dimensional slice through τ , where all parametersexcept one are kept constant. For example, if all parameters except the PCM Grashofnumber Gr p were held constant, τ would be given as τ | Re f ,P r f ,P r p = constant ( Gr p ) = τ pg (17)where the subscript p denotes PCM and the additional subscript g denotes that the Grashofnumber is the variable in question. To verify the hypothesis of (16), the data are plottedwith this and other time scalings in Figure 7. In figure 8, it is observed that Gr p F o p collapsesthe data to a single curve provided that P r p is constant but from figure 9 it is apparent thatthe collapse also occurs for multiple values of P r p . Similarly, figure 10 shows that P r / p F o p collapses the data to a single curve provided that Gr p is constant. These two relationships21re combined in figure 7 and time is scaled according to (16) to collapse to almost a singlecurve all 16 cases having the same values for Re f and P r f .A question that cannot be answered with the approach in § η critical .This value is needed in order to inform whether the linear or asymptotic scaling of η withtime is appropriate. From observing figure 7, we estimate η critical ≈ .
9. The existence ofthe linear and inverse-exponential regions is further supported by figure 14, which clearlyshows the collapsed curves diverging as straight lines on the onset of the inverse-exponentialregion, as expected on a log-linear axes. Additionally, figures 13, 8 and 9 show that η τ transitions to non-linear around η critical . We have noted in several places in this paperthat the simulations are limited to laminar flow in the PCM whereas practical systems mayemploy turbulent flow. We know the reason, though, why the value of η critical will dependon the hydrodynamic regime of the PCM. In the remaining parts of this section, we discussthe scaling obtained and its practical implications. Gr p and P r p : Discussion and implications In this section, we discuss the reasons for the scaling obtained in figure 7 and the im-plications of this scaling for the design and operation of LHTES devices. Figures 8 and 9shows plots for various cases from table 7 and we can see that the different cases collapseto one curve when the time is scaled with the PCM Grashof number, which agrees withour predictions in 5.1. This is consistent with what Jany and Bejan [14] observed by scalinganalysis for mixed conduction-convection flow regimes in an enclosure with laminar flow. Onreaching η critical , the curve changes shape from linear to asymptotic, as predicted in 5.2. Inorder to further confirm our hypothesis from 5.1, we plot the measured Grashof number infigure 11, which indicates that the PCM container wall temperature T mci is indeed constantfor the cases under consideration. Figure 10 further confirms P r / p scaling obtained in figure7. This scaling corresponds to the P r / scaling observed in laminar forced convection withuniform heat flux. As explained in section 5.1, both the temperature difference T mci − T mean and the heat transfer rate q pl are constant for majority of the time, as demonstrated by thelinearity of the melt fraction curve. Due to the Boussinesq approximation, the simulatedflow conserves volume, and a downward movement of volume must be matched by an up-ward movement. Thus, even though we cannot explain the scaling entirely, we note thatthe conditions in the PCM match those given in the analysis of Bejan [4, Eq 2.121], whichpredicts a P r / dependency.This implies that the temperature difference T in − T mean is the most important parameterfor obtaining fast heat transfer, and should be maximized. The HTF inlet temperature isconstrained by the application being studied. Thus, the temperature difference may bemaximized by picking a PCM with a lower mean melting temperature. However, increasingthis difference corresponds to a loss in quality of heat stored. The temperature T mean is alsoexpected to be an important parameter for heat transfer during discharging of the device,as it shall affect the discharging heat transfer rate. Thus, it is desirable to find an optimizedvalue of T mean that maximizes the quality stored energy and the charging and dischargingrates. 22 igure 7: The melt fraction vs dimensionless time on the parameter plane with Re f = 2090 and P r f = 5.All curves collapse, indicating that the timescale defined in (16) is appropriate if Re f is high enough and P r f is low enough, or if u f and k f are both high enough to ensure sufficient heat flow. The markers areplotted to distinguish the curves. igure 8: The scaling of melt fraction with theGrashof number on the parameter plane with Re f = 2090, P r f = 5, P r p = 14 Figure 9: Analogous to figure 8, except P r p = 6.Figure 10: Melt Fraction η as a function of P r p atconstant Gr p = 199248, P r f = 5 and Re f = 2090.The figure shows the perfect scaling with P r p / Figure 11: Measured Gr p as a function ofnon-dimensional time τ , at constant Re f = 2090, P r f = 5 and P r p = 6. The plateau shows that theassumptions from § Re f atconstant Gr p = 83020, P r f = 42 and P r p = 57. Theplots show that the heat transfer rate does notimprove much past Re = 800. Figure 13: Melt Fraction as a function of four P r f atconstant Gr p = 83020, P r p = 57 and Re f = 2090 igure 14: Dimensionless temperature plotted as 1 − θ at constant Re f = 2090, P r f = 5 and P r p = 6. Thefigure shows the existence of both linear and inverse-exponential regions as described in § § .2. Effects of Re f and P r f : Discussion and implications Figure 12 shows η τ for different values of Re f . Increasing Re f reduces the boundarylayer thickness in the HTF which increases q f . We conclude that beyond Re f = 800, there isno significant improvement in q f . Since the range of Re f presented here are in the laminarregion, another possibility is that we would get further enhancement in q f with turbulent flowin the HTF pipe, which reduces the boundary layer thickness further. However, increasingthe HTF thermal conductivity k f thus reducing the HTF Prandtl number P r f also reducesthe boundary layer thickness. Figure 10 shows that q f does not change significantly byincreasing the thermal conductivity. This indicates that for the range of Gr p presented here,at Re f = 800, there is sufficient flow of energy into the HTF domain. Thus, it is desirable tomaximize the energy available in the HTF region by selecting a HTF with sufficiently highconductivity and pumping it with sufficient velocity to remove the weak dependency of theheat transfer rate on Re f and P r f . As mentioned in §
2, the fluid velocity is a parameterthat can be easily controlled.
7. Conclusions
The parametric performance modeling of LHTES devices is essential for their effectiveuse. In this paper, we present a framework to analyse LHTES devices and apply it to a typicalshell and tube heat exchanger geometry. Out of all the parameters listed in table 1, we pickthe fundamental operating parameters u f , T in and two fundamental material parameters k p , k f and study their effect on the heat transfer rate q pl and the melt fraction η by conducting amatrix of 64 simulations. We observe that the melt fraction scales with the PCM Grashofnumber Gr p and the PCM Prandtl number P r (1 / p provided that there is sufficient energyprovided by the HTF. No significant scaling is observed for the HTF Reynolds number Re f and HTF Prandtl number P r f in the range studied and we conclude that these parametersdo not matter provided that the heat transfer rate is not limited by the HTF.The form of η versus time as the PCM melts has a linear region and a nonlinear regionwith the separation between them defined by a critical melt fraction η critical ≈ η critical for cases with vastly different parameters areobserved to be similar in shape, which suggests a universality for the critical melt fraction η critical . Based on this, we make the following conclusions about the design process forLHTES devices that shall enable the maximization of heat transfer rates.1. The HTF velocity and thermal conductivity have weak effects on the heat transfer rate,even at moderate values, provided that the HTF does not limit the overall availabilityof energy. As noted in §
6, the velocity and the choice of HTF fluid are somewhateasier to customize than the PCM parameters, and the HTF velocity is limited only byconsiderations of optimizing pumping power and reducing pipe wear. This is termed asthe ‘sufficient’ condition, and is indicated by the HTF tube walls approaching constant26emperature. The values Re f = 800 and P r f = 8 are found to be sufficient for thegeometry studied here.2. The variation of melt fraction η (which is a measure of stored latent energy) withtime consists of linear and asymptotic regions. The linear region is characterized by aconstant and higher heat transfer rate, which makes it the relevant region for operatingthe heat exchanger as an energy storage device. The critical melt fraction η critical denotes the transition between these regions, and the device should only be operatedupto that value. For the current geometry, the value of η critical is 0.9.3. Given ‘sufficient’ conditions in the HTF, the energy stored is given by the correlation aP r / p Gr p F o p where a is a constant governed by the particular choice of geometry.This correlation, applicable in the linear region, can also stated as a dimensionlesstimescale given in (16).4. The effect of the PCM Grashof number Gr p is much stronger than the PCM Prandtlnumber P r p . In terms of selecting PCM materials and operating parameters, thisindicates that varying the melting point and/or HTF inlet temperature has a strongereffect on the heat transfer rate than enhancing the thermal conductivity of the PCM.However, if the charging and discharging HTF temperatures are fixed (this is expected,since they are governed by the application), increasing the charging Gr p reduces thedischarging Gr p . The melting point of the PCM should be optimized in order tosatisfy both charging and discharging conditions. Hence, finding materials for whichthe melting point can be varied, with means such as additives or chemical composition,is indicated to be an important area for further research. This conclusion agrees withthat of the the review presented by Farid et al. [10]5. The HTF Prandtl number is a parameter that can be used to eliminate limiting HTFeffects on heat transfer. Thus, enhancement of HTF conductivity through additives isindicated as a future research subject. Acknowledgments.
This work was sponsored by the US Department of Energy grants DE-EE00007708 and DE-EE00008277. 27 eferences [1] ANSYS Fluent. Fluent User Guide. In
ANSYS FLUENT , volume 123, pages 407–408.2011.[2] Maryam Arbabzadeh, Ramteen Sioshansi, Jeremiah X Johnson, and Gregory A Ke-oleian. The role of energy storage in deep decarbonization of electricity production inCalifornia.
Nature Communications , (2019):1–35, 2019. ISSN 2041-1723.[3] Bechiri, Mohammed and Mansouri, Kacem. Analytical solution of heat transfer in ashell-and-tube latent thermal energy storage system.
Renewable Energy , 74:825–838,2015. ISSN 09601481.[4] Adrian Bejan.
Convection Heat Transfer, Fourth Edition . 03 2013.[5] Kunal Bhagat, Mohit Prabhakar, and Sandip K. Saha. Estimation of thermal perfor-mance and design optimization of finned multitube latent heat thermal energy storage.
Journal of Energy Storage , 2018. ISSN 2352152X.[6] Y Cao and A Faghri. Performance characteristics of a thermal energy storage module :a transient PCM/forced convection conjugate analysis.
Int. J. Heat and Mass Transfer ,pages 93–101, 1991.[7] Y Cao and A Faghri. A Study of Thermal Energy Storage Systems With ConjugateTurbulent Forced Convection.
Journal of Heat Transfer , pages 1019–1027, 1992.[8] P. Denholm, E. Ela, B. Kirby, and M. Milligan. Role of Energy Storage with Renew-able Electricity Generation. Technical Report January, National Renewable EnergyLaboratory, 2010.[9] G. Diarce, Campos-Celador, J. M. Sala, and A. Garc´ıa-Romero. A novel correlationfor the direct determination of the discharging time of plate-based latent heat ther-mal energy storage systems.
Applied Thermal Engineering , 129:521–534, 2018. ISSN13594311.[10] Mohammed M. Farid, Amar M. Khudhair, Siddique Ali K. Razack, and Said Al-Hallaj.A review on phase change energy storage: Materials and applications, 2004. ISSN01968904.[11] Jaume Gasia, Jan Diriken, Malcolm Bourke, Johan Van Bael, and Luisa F. Cabeza.Comparative study of the thermal performance of four different shell-and-tube heatexchangers used as latent heat thermal energy storage systems.
Renewable Energy ,2017. ISSN 18790682.[12] Kang Hu, Lei Chen, Qun Chen, Xiao-Hai Wang, Jun Qi, Fei Xu, and Yong Min. Phase-change heat storage installation in combined heat and power plants for integration ofrenewable energy sources into power system.
Energy , 124:640–651, 2017. ISSN 03605442.2813] Frank P Incropera, David P DeWitt, Adrienne S Lavine, and Theodore L Bergman.
Fundamentals of heat and mass transfer . John Wiley & Sons, 2011.[14] Peter Jany and Adrian Bejan. Scaling theory of melting with natural convection inan enclosure.
International Journal of Heat and Mass Transfer , 31(6):1221–1235, 1988.ISSN 00179310.[15] Dheeraj Kishor Johar, Dilip Sharma, Shyam Lal Soni, Pradeep K. Gupta, and RahulGoyal. Experimental investigation and exergy analysis on thermal storage integratedmicro-cogeneration system.
Energy Conversion and Management , 131:127–134, 2017.ISSN 01968904.[16] Lokesh Kalapala and Jaya Krishna Devanuri. Influence of operational and design pa-rameters on the performance of a pcm based heat exchanger for thermal energy storagea review.
Journal of Energy Storage , 20:497 – 519, 2018. ISSN 2352-152X.[17] Jun Li, Chin Hong Tam, and Guohong Tian. Investigation of an HD Engine ThermalStorage System.
Energy Procedia , 105(0):4110–4115, 2017. ISSN 18766102.[18] Martin Longeon, Ad`ele Soupart, Jean Fran¸cois Fourmigu´e, Arnaud Bruch, and PhilippeMarty. Experimental and numerical study of annular PCM storage in the presence ofnatural convection.
Applied Energy , 2013. ISSN 03062619.[19] Benjamin McDaniel and Dragoljub Kosanovic. Modeling of combined heat and powerplant performance with seasonal thermal energy storage.
Journal of Energy Storage , 7:13–23, 2016.[20] L. Mongibello, M. Capezzuto, and G. Graditi. Technical and cost analyses of twodifferent heat storage systems for residential micro-CHP plants.
Applied Thermal En-gineering , 71(2):636–642, 2014. ISSN 13594311.[21] K. Nithyanandam and R. Pitchumani. Computational Modeling of Dynamic Responseof a Latent Thermal Energy Storage System With Embedded Heat Pipes.
Journal ofSolar Energy Engineering , 2013. ISSN 0199-6231.[22] K. Rathod and Banerjee Jyotirmay. Development of Correlation for Melting Time ofPhase Change Material in Latent Heat Storage Unit. In
Energy Procedia , pages 2125–2130, 2015. ISBN 1876-6102.[23] H. Shmueli, G. Ziskind, and R. Letan. Melting in a vertical cylindrical tube: Numericalinvestigation and comparison with experiments.
International Journal of Heat and MassTransfer , 2010. ISSN 00179310.[24] Hebat-Allah M. Teamah, Marilyn F. Lightstone, and James S. Cotton. NumericalInvestigation and Nondimensional Analysis of the Dynamic Performance of a ThermalEnergy Storage System Containing Phase Change Materials and Liquid Water.
Journalof Solar Energy Engineering , 2016. ISSN 0199-6231.2925] K. P. Venkitaraj, S. Suresh, and Arjun Venugopal. Experimental study on the ther-mal performance of nano enhanced pentaerythritol in IC engine exhaust heat recoveryapplication.
Applied Thermal Engineering , 137(October 2017):461–474, 2018. ISSN13594311.[26] V. R. Voller and C. Prakash. A fixed grid numerical modelling methodology forconvection-diffusion mushy region phase-change problems.
International Journal of Heatand Mass Transfer , 1987. ISSN 00179310.[27] B. Yimer and M. Adami. Parametric study of phase change thermal energy storagesystems for space application.