Numerical Evaluation of the Transient Performance of Rock-Pile Seasonal Thermal Energy Storage Systems Coupled with Exhaust Heat Recovery

This study seeks to investigate the concept of using large waste rocks from mining operations as waste-heat thermal energy storage for remote arctic communities, both commercial and residential. It holds its novelty in analyzing such systems with an experimentally validated transient three-dimensional computational fluid dynamics and heat transfer model that accounts for interphase energy balance using a local thermal non-equilibrium approach. The system performance is evaluated for a wide range of distinct parameters, such as porosity between 0.2 and 0.5, fluid velocity from 0.01 to 0.07 m/s, and the aspect ratio of the bed between 1 and 1.35. It is demonstrated that the mass flow rate of the heat transfer fluid does not expressively impact the total energy storage capacity of the rock mass, but it does significantly affect the charge/discharge times. Finally, it is shown that porosity has the greatest impact on both fluid flow and heat transfer. The evaluations show that about 540 GJ can be stored on the bed with a porosity of 0.2, and about 350 GJ on the one with 0.35, while the intermediate porosity leads to a total of 450 GJ. Additionally, thermal capacity is deemed to be the most important thermophysical factor in thermal energy storage performance.


Introduction
Increasing population and economic growth in the past years has significantly amplified fossil fuel dependency and its impact on the environment [1]. This has caused industries to search for alternatives to fossil fuel and more environmentally-friendly solutions for energy generation. With the increase in usage of renewable and alternative technologies, energy storage has gained great importance [2,3]. Thermal energy storage (TES) allows for energy to be stored in the form of heat in times of low or null demand, and used later whenever demand rises [2]. Such technologies aid in the reduction of both carbon footprints and energy costs from small (residential) to big (industrial) scales, being even deemed vital in some applications [4,5]. Such technologies are mostly used in solar energy applications, mainly concentrated solar power plants (CSP) [3,[6][7][8] and solar district heating [9][10][11][12][13], but waste-heat storage systems can also take great advantage of their use [14][15][16][17]. Among the three main types of TES systems: chemical heat storage (based on thermochemical reactions); latent heat storage (which employs phase change materials); and sensible heat storage (simply based on temperature increase/decrease of the storage media) [2], the latter is currently the cheapest and easiest one to implement [5]. It has been stated that a sensible heat storage costs up to 15 CAD/kWh, while a latent heat storage needs relatively higher costs, up to 150 CAD/kWh, depending on many factors such as the size, the material, and the applied technology [18]. For that reason, it is the most promising one for immediate use in large/industrial scales [19]. That happens as it generally implements inexpensive resources as fragmented rock and air [20], and it often dispenses the need of intermediate heat exchangers [16].
It has been shown previously that rock-bed storage systems frequently are more cost-effective than similar water-based systems for heating [11]. Still, other advantages of such systems can be highlighted. They are more environmentally friendly than their counterparts, as they use naturally occurring storage media, they can usually withstand very high temperatures, and they also have a high heat transfer coefficient between phases when in use, but low conductivity of heat when air flow is not present [12]. Despite being used mostly for solar applications, due to its flexibility, packed rock-bed TES can be employed in a wide range of energy-intensive industries. A good example is the development of the application for mine sites, where there is an opportunity to use fragmented rock packed in beds as a sensible energy storage. In these storages, air is used as the heat transfer fluid (HTF) and fragmented rocks are used as the storage media, which has been shown as a favorable approach for TES systems [21][22][23][24][25].
In order to fully understand the performance of TES systems, several parameters ought to be analyzed, which are hard to attain by analytical equations. Hence, computational simulations are usually a good method to investigate and enhance such systems, by improving the many factors affecting its operation. For that reason, several works in literature have taken the task to develop such computational models. Some of these that can be highlighted focused on the heat and mass transfer phenomena occurring within packed beds used as thermal energy storage [26]. Abdel-Salam et al. established a cylindrical coordinate-based 1D numerical model, through which it was revealed that the increase of bed length improves storage capacity while not affecting the rate of charge [4]. By implementing the finite difference method in models with distinct number of sections, Arias et al. investigated the possibility of using them to represent the stratification phenomena seen in experimental tests [27]. Hänchen et al. elaborated a one-dimensional numerical model, followed by comparison validation against experimental data from small-scale testing [28]. Looking to investigate the effect of construction parameters on heat storages for industrial processes, Zanganeh et al. developed a quasi-one-dimensional transient model to study two-phase heat transfer on packed beds [16]. Barton developed a one-dimensional model for airflow within a porous rock bed medium, and used it to study the effect of different charging and discharging flow direction schedules [29]. Kuravi et al. studied a large-sized packed bed (created with bricks) circulating air as an HTF with experiments and a 1D energy model, and confirmed the viability of its fluid flow and heat transfer for solar applications [30]. Finally, Cascetta et al. worked on experimental and numerical tests, focusing on their comparison, and achieved satisfactory results on simulations that had input temperatures from experiments [31].
When it comes to design and operation of porous media storage systems, two important parameters ought to be highlighted: the rate of heat transfer between fluid and solid domains, and the power necessary to run the heat transfer fluid through the system. The former determines operating cost, while the latter is directly related to storage capacity and rate. Besides, other case-specific factors can directly impact the outcome of the system. Hence, it is crucial to fully comprehend the porous structure of the rock bed domains and how important they are for the design of heat storage systems. Some studies have focused on evaluating influential parameters. Optimizing the heat transfer rate between the fluid and solid domains has been demonstrated to be a key factor to achieve a good performance on TES systems [32]. Mertens et al. [8] and Rodat et al. [33] investigated, in detail, both numerically and experimentally, the performance of sand and rock beds, which allowed for the impact of several design and operating conditions, as temperature and HTF velocity, to be evaluated on thermal energy storage performance. Air mass flow has been shown to enhance interphase heat transfer rates, decreasing charging time [4], while some have suggested that air density would not directly affect said heat transfer [29]. It has also been stated that volumetric heat capacity would play a major role on overall system performance [28]. Regarding insulation, a thickness of about 16-20 cm was found to be enough to retain low energy losses for intermediate temperatures [16]. Studies have also shown that for short-term (diurnal) cycles, a 95% overall thermal efficiency is achievable [7]. Another study was performed to evaluate rock suitability for use on TES, regarding chemical, mechanical, and thermophysical properties, and it found that some rocks which are not thermally stable (like foliated rocks) could be unsuitable. Finally, it has been stated that there is no absolute optimum design to a rock-bed TES system, meaning every unit needs to be planned and directed towards the requirements for cost and performance of the specific application.
It is important to note that the majority of experimental data available in literature focuses on lab-scale models. Not many studies exist in literature investigating the heat transfer phenomena on large-scale rock-bed TES as the ones with potential for implementation in mining applications. Only recently has thought been put into uses of rock beds, such as cooling/heating of intake air in underground mining sites and block-caved mine ventilation systems [21,24]. It is very important to investigate pressure drop on those systems, as significant additional fan power is required to pump the air through the bed [34,35], which leads to meaningfully high costs that could affect the viability of the system.
Even though several studies have modelled heat transfer and fluid flow in a packed bed of rocks, most of them present restrictions, including a dimensional reduction to a 2D, 1D, or quasi-1D model, neglecting heat transfer in at least one of the three special dimensions. Others have employed simplifications from interphase energy balance (obtained by the local thermal non-equilibrium method (LTNE)) to a less realistic local thermal equilibrium (LTE). A few even neglect pressure differentials. This study focuses on a TES system being charged by a stream of hot air, the product of the exhaust from a local power plant. The novel approach proposed here is the use of large rocks (possibly waste rock from mining operations) as storage medium while using streams of hot gases (such as those from exhaust) as the heat source of the thermal energy storage system. The latter contains energy commonly discarded to the ambient air at establishments employing diesel-based power plants, such as the ones present on remote communities in Canada [15,17]. Not only commercial (as remote mining operations), these establishments that spread over the arctic regions in Canada commonly house first nation groups. In fact, it is estimated that about 300 remote communities of this kind exist in Canada alone [36]. It is important to note that this implementation of TES in large-scale rock beds compares to previous work done on natural heat exchange areas for mining applications [24,25], simply based on particle size of storage media and heat transfer mechanisms involved, differing completely in heat source and grade, application, charging/discharging rates, heat transfer rates, and nature of the fluid flow within the porous media. Therefore, this research holds its novelty on the investigation of the effects of several physical characteristics of a rock-bed TES charged through waste heat for space heating applications. Such a task is achieved employing a three-dimensional, transient model, based on the LTNE approach, by which heat transfer and fluid flow within the packed bed thermal storage system are studied through the use of computational fluid dynamics (CFD) models. The impact of several parameters on the performance of the large-sized, rock-bed waste-heat storage for remote communities is evaluated for several distinct scenarios, while focusing on the simultaneous impact of these with changes in system dimensions, i.e., aspect ratios (AR), which to the best knowledge of the authors, has not been performed before.

Model Description
The thermal storage system proposed here is based on a conical geometry with a truncated top. This format is selected due to the facilitation it offers for construction, since it can be easily achieved by stacking the chosen rock with trucks, complying to a determined angle of repose. Thus, this system has been referred to as a rock-pile TES [15], and is represented in Figure 1. Note that the directions of the flow and fluid inlet are inverted for charging (top to bottom) and discharging (bottom to top), to take advantage of the stratification of temperature within the HTF, as previously pointed out in literature [27].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 19 take advantage of the stratification of temperature within the HTF, as previously pointed out in literature [27].
(a) (b) As fluid flow develops within the rock bed, several heat transfer mechanisms will be present. To fully understand the operation of the system proposed, it is crucial to be attentive to the existing heat transfer mechanisms taking place. These are indicated in Figure 2. One of the main advantages of the technology presented is the fact that elevated heat transfer rates are achievable during operation due to the large total interphase heat transfer area. While in non-operating periods, with the absence of flow, the heat conduction and, consequently, dispersion and heat losses can be maintained at lower levels than other similar TES technologies. It is worth noting that, here, the solid phase is assumed to have constant thermophysical properties along the bed. In order to analyze the rock-bed TES, a volume-averaged model (VAM) was developed. Using this model, it is possible to evaluate the mass and heat transfer phenomena within the bed, while requiring lower computational power relative to other methods. This is a consequence of it using macroscale flow characteristics in the porous media to evaluate the system's outcomes. The viscous and inertial resistance coefficients needed for such a VAM model are obtained from empirical equations. Figure 3 depicts a representative elementary volume (REV), a typical unit used in the volume-averaged analysis. It is important to note that it is only suitable to use VAM for systems with As fluid flow develops within the rock bed, several heat transfer mechanisms will be present. To fully understand the operation of the system proposed, it is crucial to be attentive to the existing heat transfer mechanisms taking place. These are indicated in Figure 2. One of the main advantages of the technology presented is the fact that elevated heat transfer rates are achievable during operation due to the large total interphase heat transfer area. While in non-operating periods, with the absence of flow, the heat conduction and, consequently, dispersion and heat losses can be maintained at lower levels than other similar TES technologies. It is worth noting that, here, the solid phase is assumed to have constant thermophysical properties along the bed.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 19 take advantage of the stratification of temperature within the HTF, as previously pointed out in literature [27].
(a) (b) As fluid flow develops within the rock bed, several heat transfer mechanisms will be present. To fully understand the operation of the system proposed, it is crucial to be attentive to the existing heat transfer mechanisms taking place. These are indicated in Figure 2. One of the main advantages of the technology presented is the fact that elevated heat transfer rates are achievable during operation due to the large total interphase heat transfer area. While in non-operating periods, with the absence of flow, the heat conduction and, consequently, dispersion and heat losses can be maintained at lower levels than other similar TES technologies. It is worth noting that, here, the solid phase is assumed to have constant thermophysical properties along the bed. In order to analyze the rock-bed TES, a volume-averaged model (VAM) was developed. Using this model, it is possible to evaluate the mass and heat transfer phenomena within the bed, while requiring lower computational power relative to other methods. This is a consequence of it using macroscale flow characteristics in the porous media to evaluate the system's outcomes. The viscous and inertial resistance coefficients needed for such a VAM model are obtained from empirical equations. Figure 3 depicts a representative elementary volume (REV), a typical unit used in the volume-averaged analysis. It is important to note that it is only suitable to use VAM for systems with In order to analyze the rock-bed TES, a volume-averaged model (VAM) was developed. Using this model, it is possible to evaluate the mass and heat transfer phenomena within the bed, while requiring lower computational power relative to other methods. This is a consequence of it using macroscale flow characteristics in the porous media to evaluate the system's outcomes. The viscous and inertial resistance coefficients needed for such a VAM model are obtained from empirical equations. Figure 3 depicts a representative elementary volume (REV), a typical unit used in the volume-averaged analysis. It is important to note that it is only suitable to use VAM for systems with l ≥ 5d and l L [37,38]; in which d is the average diameter of the pore, and L and l are, respectively, the representative size of the flow domain and the REV.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 19 ≥ 5 and l ≪ L [37,38]; in which is the average diameter of the pore, and L and l are, respectively, the representative size of the flow domain and the REV. Using ANSYS Fluent 17.2, the 3D computational model was created, meshed, and finally simulated. To solve the governing equations, the semi-implicit method for pressure-linked equation (SIMPLE) was implemented. Basically, two different methods can be used to solve the equation for conservation of energy in a porous media: a local thermal equilibrium (LTE) approach, and a local thermal non-equilibrium (LTNE) approach. When the thermal properties of the solid and the fluid are not comparable and the flow conditions through the system are highly transient, the latter model should be used. In said method, since the variation of the local rate of temperature for fluid and solid are different, two separate energy equations are solved. Hence, the heat transfer phenomena for both fluid and solid phases can be examined in the proposed thermal energy storage system. Accordingly, the following mass, momentum, and energy conservation equations (respectively) are then solved by the finite-volume analysis under LTNE conditions.
where is the viscous resistance coefficient and is the inertial one, as given by Equations (5) and (6), respectively. Note that α represents the inverse of the permeability of the medium (K). Distinct porous bodies with increasing porosities ( ) from 0.2 to 0.5 and a particle size ( ) were used in this study. To find values for α and β, the commonly employed Ergun correlation [39] is used. It is important that numerical and experimental studies be performed to obtain such values. Only then should Ergun correlation be assessed for each bed, and its validity confirmed.
Thus, having in mind its TES performance, a three-dimensional VAM is suggested to evaluate the main outcomes of the system, such as pressure drop/fan power requirements, total energy Using ANSYS Fluent 17.2, the 3D computational model was created, meshed, and finally simulated. To solve the governing equations, the semi-implicit method for pressure-linked equation (SIMPLE) was implemented. Basically, two different methods can be used to solve the equation for conservation of energy in a porous media: a local thermal equilibrium (LTE) approach, and a local thermal non-equilibrium (LTNE) approach. When the thermal properties of the solid and the fluid are not comparable and the flow conditions through the system are highly transient, the latter model should be used. In said method, since the variation of the local rate of temperature for fluid and solid are different, two separate energy equations are solved. Hence, the heat transfer phenomena for both fluid and solid phases can be examined in the proposed thermal energy storage system. Accordingly, the following mass, momentum, and energy conservation equations (respectively) are then solved by the finite-volume analysis under LTNE conditions.
where α is the viscous resistance coefficient and β is the inertial one, as given by Equations (5) and (6), respectively. Note that α represents the inverse of the permeability of the medium (K). Distinct porous bodies with increasing porosities (ε) from 0.2 to 0.5 and a particle size (d) were used in this study.
To find values for α and β, the commonly employed Ergun correlation [39] is used. It is important that numerical and experimental studies be performed to obtain such values. Only then should Ergun correlation be assessed for each bed, and its validity confirmed.
Thus, having in mind its TES performance, a three-dimensional VAM is suggested to evaluate the main outcomes of the system, such as pressure drop/fan power requirements, total energy storage, charging rate, and temperature profiles of the bed, as a function of thermodynamic and physical parameters for different geometries (i.e., aspect ratios).

Heat Transfer Coefficient (HTC) in the Rock Bed
As previously mentioned, regarding the thermal energy storage systems, it is crucial to determine the convective heat transfer coefficient (HTC) between solid and fluid. Hence, it is fundamental to assume appropriate correlations to appropriately estimate the performance of such systems. Here, the correlations used are as given by Bejan [38], which are presented in Equations (7)-(9). 1

Geometry and Aspect Ratio
When designing thermal energy storage systems, a trade-off is usually encountered between heat losses and pressure losses. These factors are also usually related to the physical footprint of the system, i.e., space restrictions. However, studies show that it is possible to find an ideal aspect ratio, like the best one for a truncated cone-shaped TES system (γ = H/D upper , where H is height and D upper is upper diameter) being higher than one [26]. It is important to note that economic and technical limitations are commonly related to the mentioned challenges and directly affect the design of the system. Having all that in mind, three aspect ratios were selected for the rock bed proposed here. The ratios range from 1.0 to 1.35 and were fashioned based on a constant volume of 1830 m 3 . The truncated cone-shaped domains have three distinct heights, 10 m, 12 m, and 13.5 m, and can be seen in Figure 4a-c, respectively. All shapes have an invariable upper diameter (D lower ) of 10 m, and a lower diameter ranging from 20 to 16 m.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 19 storage, charging rate, and temperature profiles of the bed, as a function of thermodynamic and physical parameters for different geometries (i.e., aspect ratios).

Heat Transfer Coefficient (HTC) in the Rock Bed
As previously mentioned, regarding the thermal energy storage systems, it is crucial to determine the convective heat transfer coefficient (HTC) between solid and fluid. Hence, it is fundamental to assume appropriate correlations to appropriately estimate the performance of such systems. Here, the correlations used are as given by Bejan [38], which are presented in Equations (7)- (9).

Geometry and Aspect Ratio
When designing thermal energy storage systems, a trade-off is usually encountered between heat losses and pressure losses. These factors are also usually related to the physical footprint of the system, i.e., space restrictions. However, studies show that it is possible to find an ideal aspect ratio, like the best one for a truncated cone-shaped TES system ( = / , where H is height and Dupper is upper diameter) being higher than one [26]. It is important to note that economic and technical limitations are commonly related to the mentioned challenges and directly affect the design of the system. Having all that in mind, three aspect ratios were selected for the rock bed proposed here. The ratios range from 1.

Thermophysical Properties
In this study, the rock thermophysical properties were considered constant, while the air thermophysical properties were defined as functions of temperature. The main properties for the storage domain (solid material) are compiled in Table 1. As the proposed TES in this study aims to be charged with high-temperature exhaust air around 120 • C, and discharged with ambient-temperature in winter around 7 • C, it is crucial to correlate the properties of air (i.e., density, thermal conductivity, and viscosity) with temperature. Numerous equations have been introduced in literature to compute the temperature-dependent thermal properties of air. Here, the following equations (Equations (10)-(12)) were used, which are valid in a range of −73 • C-200 • C [40].
where ρ f is air density (kg/m 3 ), µ f is air viscosity (10 −6 N·s/m 2 ), k f is air thermal conductivity (W/m·K), and T is the air temperature (Kelvin). Thus, user-defined functions (UDF) written in C language were used to account for the air thermal properties as a function of temperature.

Boundary Conditions
The current study takes into account the transient fluid flow that happens within the solid domain when it is crossed by the air. In the charging phase, the bottom of the bed is defined as pressure outlet, while the top is defined as a velocity inlet. Meanwhile, during the discharge phase, the boundaries are inverted and the ambient air (i.e., cold air) is injected from the bottom of the TES system. Seeking to take advantage of the possible buoyancy effects that may happen when the cold and hot air are present inside the rock bed, the hot inlet air (at fixed 120 • C) is injected through the top of the system, whereas the cold air in the discharge stage is injected through the bottom. Besides, the lateral boundaries of the system are considered insulated (zero heat flux) and with no slip (zero velocity).

Mesh Independency Study
In order to make sure the results obtained independent of the mesh refinement for the model developed, three meshes with distinct levels of refinement were created with 8.5 × 10 4 , 1.54 × 10 5 , and 2.86 × 10 5 elements, respectively. The mesh containing about 1.54 × 10 5 elements results in a 2% deviation in terms of pressure drop and outlet temperature, compared to the mesh with around 2.86 × 10 5 elements. The deviation in results is found to be about 9% between the coarsest mesh (containing 8.5 × 10 4 elements) and the finest one (2.86 × 10 5 elements). Thus, the intermediate mesh consisting of around 1.54 × 10 5 elements was chosen for all the further numerical investigations.

Validation of the Model
An essential step of the implementation of heat transfer and fluid flow model on the prediction of real-case system performance is the validation of its results against outcomes of an experimental system. It is essential to prove that the proposed model is able to capture the real phenomena experienced by the porous media domain while used as TES. Here, the pilot-test data presented by Hänchen et al. [28] was used for experimental validation. Dimensions were rescaled, and the boundary conditions were readjusted to conform with the heat storage system proposed there. Hence, the model was run considering both the LTE and LTNE approaches. This was performed seeking the best match for the numerical model and the experimental results. The results are shown in Figure 5a for temperature distribution within the rock bed at different points in time, allowing for the validation of heat transfer, while Figure 5b depicts a comparison of pressure drop development across the bed with time, allowing validation of fluid flow. Both temperature and pressure can be seen to agree well with the experimental results. It is important to note that, although LTE and LTNE results are very similar for pressure drop, they are significantly different for temperature distribution across the bed, with LTNE being closer to the experimental ones. This is particularly true as temperatures develop with time. For that reason, the LTNE method, despite being more computationally intensive, was chosen for all the following simulations presented in this study.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 19 system. It is essential to prove that the proposed model is able to capture the real phenomena experienced by the porous media domain while used as TES. Here, the pilot-test data presented by Hänchen et al. [28] was used for experimental validation. Dimensions were rescaled, and the boundary conditions were readjusted to conform with the heat storage system proposed there. Hence, the model was run considering both the LTE and LTNE approaches. This was performed seeking the best match for the numerical model and the experimental results. The results are shown in Figure 5a for temperature distribution within the rock bed at different points in time, allowing for the validation of heat transfer, while Figure 5b

Aspect Ratios (Rock-Bed Shape)
It is expected that the dimensions of the TES system will affect the configuration of the fluid flow profiles and velocity contours within the porous media, thus impacting the ratio of heat transfer and the capacity of the system to store heat. Due to that, the distinct geometries previously presented were used throughout the study to show the impact of several parameters in the fluid flow and heat transfer of rock beds of different dimensions. Figure 6 depicts the contours for velocity inside the rock bed, in the storage period, for all the proposed aspect ratios. The selected values for inlet velocity being 0.05 m/s, the porosity 0.2, and an average particle size of 1.0 m (here considered the base case). It can be seen that having a lower aspect ratio causes higher velocity gradients to happen, leading to a lower exit velocity of air. This is reasonable, since the outlet has a greater area in that case. It should be noted that, for higher mass flows (and consequently higher Reynolds numbers), the exploitation rate, along with the energy efficiency of the system, declines [41]. The reason for this being that the residence time of the HTF inside the bed is reduced as an effect of the lower fluid velocities, which,

Aspect Ratios (Rock-Bed Shape)
It is expected that the dimensions of the TES system will affect the configuration of the fluid flow profiles and velocity contours within the porous media, thus impacting the ratio of heat transfer and the capacity of the system to store heat. Due to that, the distinct geometries previously presented were used throughout the study to show the impact of several parameters in the fluid flow and heat transfer of rock beds of different dimensions. Figure 6 depicts the contours for velocity inside the rock bed, in the storage period, for all the proposed aspect ratios. The selected values for inlet velocity being 0.05 m/s, the porosity 0.2, and an average particle size of 1.0 m (here considered the base case).
It can be seen that having a lower aspect ratio causes higher velocity gradients to happen, leading to a lower exit velocity of air. This is reasonable, since the outlet has a greater area in that case. It should be noted that, for higher mass flows (and consequently higher Reynolds numbers), the exploitation rate, along with the energy efficiency of the system, declines [41]. The reason for this being that the residence time of the HTF inside the bed is reduced as an effect of the lower fluid velocities, which, per se, decreases the opportunity for heat exchange. Therefore, it is important to note that, to improve heat exchange rate within the system, an optimum fluid velocity should be found for each application.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 19 heat exchange rate within the system, an optimum fluid velocity should be found for each application.

Effect of Porosity, Particle Size and Fluid Velocity
Simulations were performed with distinct porosities, particle diameter, and velocities for all the aspect ratios. The HTF inlet velocity was tested from 0.01 to 0.07 m/s, while the porosity was varied from 0.2 to 0.5. Results from simulations with distinct porosities and aspect ratios are shown in Figure  7a-c for the charging, and Figure 7d-f for the discharging, cycles. It can be seen that the outlet temperature for the air continuously increases during the charging period up to almost 120 °C, which makes the inlet and outlet temperature differential no more than 1 °C. It is also visible that when using greater fluid velocities (and consequently larger mass flows), temperature raise is much faster becoming practically independent from the porosity of the bed and the aspect ratio. The same can be said for the discharging cycle, in which the drop in outlet temperature is equally fast.
Having in mind the obtained outcomes, it is possible to calculate the total energy stored in the system. The evaluations show that about 540 GJ can be stored on the bed with a porosity of 0.2, and about 350 GJ on the one with 0.35, while the intermediate porosity leads to a total of 450 GJ. The idea is that this heat that was stored can be withdrawn from the TES system in a cold season, i.e., when demand is excessive, and be employed in a wide range of potential applications. Figure 8 shows that, by increasing porosity in the system, a direct decrease in total energy stored is encountered, mainly due to the reduction of mass for a domain with constant total volume. On the other hand, aspect ratio and rock particle size are not seen to have significant effects in the stored/extracted capacity, which is the reason why Figure 8 only shows results for one case of each. It is also worth noting that, even though the total energy stored in the system slightly increases for higher mass flow rates, lower velocities can provide a longer charging/discharging period for the TES system. Finally, it is possible to conclude that the total heat storing capacity of the system is virtually the same for all the distinct aspect ratios and Reynolds tested, which leads us to believe that such capacity is nearly entirely dependent on the mass of the storage media.

Effect of Porosity, Particle Size and Fluid Velocity
Simulations were performed with distinct porosities, particle diameter, and velocities for all the aspect ratios. The HTF inlet velocity was tested from 0.01 to 0.07 m/s, while the porosity was varied from 0.2 to 0.5. Results from simulations with distinct porosities and aspect ratios are shown in Figure 7a-c for the charging, and Figure 7d-f for the discharging, cycles. It can be seen that the outlet temperature for the air continuously increases during the charging period up to almost 120 • C, which makes the inlet and outlet temperature differential no more than 1 • C. It is also visible that when using greater fluid velocities (and consequently larger mass flows), temperature raise is much faster becoming practically independent from the porosity of the bed and the aspect ratio. The same can be said for the discharging cycle, in which the drop in outlet temperature is equally fast.
Having in mind the obtained outcomes, it is possible to calculate the total energy stored in the system. The evaluations show that about 540 GJ can be stored on the bed with a porosity of 0.2, and about 350 GJ on the one with 0.35, while the intermediate porosity leads to a total of 450 GJ. The idea is that this heat that was stored can be withdrawn from the TES system in a cold season, i.e., when demand is excessive, and be employed in a wide range of potential applications. Figure 8 shows that, by increasing porosity in the system, a direct decrease in total energy stored is encountered, mainly due to the reduction of mass for a domain with constant total volume. On the other hand, aspect ratio and rock particle size are not seen to have significant effects in the stored/extracted capacity, which is the reason why Figure 8 only shows results for one case of each. It is also worth noting that, even though the total energy stored in the system slightly increases for higher mass flow rates, lower velocities can provide a longer charging/discharging period for the TES system. Finally, it is possible to conclude that the total heat storing capacity of the system is virtually the same for all the distinct aspect ratios and Reynolds tested, which leads us to believe that such capacity is nearly entirely dependent on the mass of the storage media. Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 19 Charge phase Discharge phase  To evaluate the impact of these parameters on the pressure drop through the bed, the Reynolds number of the flow was varied from 100 to 750 and, through numerical calculations, the pressure drop across the system was calculated numerically. The results that can be seen in Figure 9 show the development of pressure gradient across the bed with Reynolds number. In addition, Figure 10 shows contours for maximum pressure gradient across the system against aspect ratio and inlet velocity for a constant particle size of 1.0 m. Figure 9 shows that, even though pressure drop increases with velocity for all scenarios, the growth is very subtle for higher porosities, being even an order of magnitude lower than that for a low porosity (and consequently low permeability). The same happens for the slight increase in pressure difference seen in larger aspect ratios, which only becomes significant at lower porosities. Perhaps this effect can be more clearly identified in Figure 10, as it shows the maximum pressure gradient contours growing rapidly with velocity, as opposed to γ, but being considerably low except for the lowest porosity. With that in mind, it can be understood that the permeability of the packed bed (direct consequence of the domain particle size and porosity) impacts the pressure drop in a more accentuated way than the change in shape does. It is important to mention that greater velocity of the HTF causes more significant pressure decrease across the bed; this means more power will be necessary from fans generating the flow for two aspects, to overcome the higher pressure drop and to create the larger mass flow itself.
(a) dp = 0.8 m To evaluate the impact of these parameters on the pressure drop through the bed, the Reynolds number of the flow was varied from 100 to 750 and, through numerical calculations, the pressure drop across the system was calculated numerically. The results that can be seen in Figure 9 show the development of pressure gradient across the bed with Reynolds number. In addition, Figure 10 shows contours for maximum pressure gradient across the system against aspect ratio and inlet velocity for a constant particle size of 1.0 m. Figure 9 shows that, even though pressure drop increases with velocity for all scenarios, the growth is very subtle for higher porosities, being even an order of magnitude lower than that for a low porosity (and consequently low permeability). The same happens for the slight increase in pressure difference seen in larger aspect ratios, which only becomes significant at lower porosities. Perhaps this effect can be more clearly identified in Figure 10, as it shows the maximum pressure gradient contours growing rapidly with velocity, as opposed to γ, but being considerably low except for the lowest porosity. With that in mind, it can be understood that the permeability of the packed bed (direct consequence of the domain particle size and porosity) impacts the pressure drop in a more accentuated way than the change in shape does. It is important to mention that greater velocity of the HTF causes more significant pressure decrease across the bed; this means more power will be necessary from fans generating the flow for two aspects, to overcome the higher pressure drop and to create the larger mass flow itself. To evaluate the impact of these parameters on the pressure drop through the bed, the Reynolds number of the flow was varied from 100 to 750 and, through numerical calculations, the pressure drop across the system was calculated numerically. The results that can be seen in Figure 9 show the development of pressure gradient across the bed with Reynolds number. In addition, Figure 10 shows contours for maximum pressure gradient across the system against aspect ratio and inlet velocity for a constant particle size of 1.0 m. Figure 9 shows that, even though pressure drop increases with velocity for all scenarios, the growth is very subtle for higher porosities, being even an order of magnitude lower than that for a low porosity (and consequently low permeability). The same happens for the slight increase in pressure difference seen in larger aspect ratios, which only becomes significant at lower porosities. Perhaps this effect can be more clearly identified in Figure 10, as it shows the maximum pressure gradient contours growing rapidly with velocity, as opposed to γ, but being considerably low except for the lowest porosity. With that in mind, it can be understood that the permeability of the packed bed (direct consequence of the domain particle size and porosity) impacts the pressure drop in a more accentuated way than the change in shape does. It is important to mention that greater velocity of the HTF causes more significant pressure decrease across the bed; this means more power will be necessary from fans generating the flow for two aspects, to overcome the higher pressure drop and to create the larger mass flow itself.

Effect of Rock Thermophysical Properties
Bear in mind that for a constant volume of rock (i.e., constant mass of storage material in the domain), an improvement in the natural capacity of the material to store thermal energy (i.e., ρ and cp) would lead to an enhancement of the TES system. For that reason, different rock types were investigated for the task proposed. Five types of rock were selected, and their thermophysical properties were taken from literature [42] and are given in Table 2. The results for outlet fluid temperature during the charging cycle can be seen in Figure 11. It is noticeable that the highest impacting factor is the thermal capacity of the rock (the product of density and specific heat). As shown, sandstone (lowest thermal capacity) quickly achieves its maximum temperature, while gabbro (highest thermal capacity) takes the longest time to be fully charged. This can also be noticed from the similarity between the quartzite and the granite temperature profiles, even with very distinct values for k. The thermal conductivity has a subtle effect on the charging rate, being slightly more accentuated for smaller aspect ratios. Except for that, the aspect ratio has a small impact on charging rate, mainly caused by the higher velocities of fluid within the bed (as previously mentioned in Section 3.2.1). However, its effect on the total energy storage is virtually null, as shown in Figure 12. Bear in mind that for a constant volume of rock (i.e., constant mass of storage material in the domain), an improvement in the natural capacity of the material to store thermal energy (i.e., ρ and cp) would lead to an enhancement of the TES system. For that reason, different rock types were investigated for the task proposed. Five types of rock were selected, and their thermophysical properties were taken from literature [42] and are given in Table 2. The results for outlet fluid temperature during the charging cycle can be seen in Figure 11. It is noticeable that the highest impacting factor is the thermal capacity of the rock (the product of density and specific heat). As shown, sandstone (lowest thermal capacity) quickly achieves its maximum temperature, while gabbro (highest thermal capacity) takes the longest time to be fully charged. This can also be noticed from the similarity between the quartzite and the granite temperature profiles, even with very distinct values for k. The thermal conductivity has a subtle effect on the charging rate, being slightly more accentuated for smaller aspect ratios. Except for that, the aspect ratio has a small impact on charging rate, mainly caused by the higher velocities of fluid within the bed (as previously mentioned in Section 3.2.1). However, its effect on the total energy storage is virtually null, as shown in Figure  12.  Figure 11. Outlet temperature in charging phase for distinct types of rock. Figure 11. Outlet temperature in charging phase for distinct types of rock.

Figure 12.
Total energy stored versus thermal capacity of the rocks simulated, and aspect ratio.

Carbon Emission Reduction
To further demonstrate the benefits of the system, the carbon footprint reduction introduced by its use was evaluated using emission data from literature [43,44]. Assuming the system substitutes natural gas as a heating source, which has carbon emissions of approximately 510 tons of carbon dioxide equivalent (CO2eq) per GWh, the total amount of carbon footprint reduction was estimated. It was found that for the base properties, the proposed TES (with a volume of 1830 m 3 , corresponding to 5490 tons of rock) would be able to reduce emissions in about 156 tons of CO2eq per year, for all the aspect ratios. This is equivalent to 85 tons of CO2eq/m 3 of rock per year, or 28 tons of CO2eq/ton of rock per year. It is important to note that for a system replacing coal (∼905 tons of CO2eq per GWh) or heating oil (∼730 tons of CO2eq per GWh), the decrease in carbon footprint could be even more significant, these being displayed in Table 3.

Conclusions
In this study, a three-dimensional transient model was established to study the heat transfer and fluid flow in a thermal energy storage system consisting of a truncated cone-shaped bed of packed rocks. A considerable number of simulations were performed to account for the impact of all the distinct parameters (aspect ratio, HTF velocity, bed porosity, particle size, and thermophysical properties) and their combination in performance of the system. The study shows that porous media properties of the system, such as porosity, are extremely important in the performance of the system as a TES system, regarding both charging rate, energy stored/extracted, and pressure drop. The aspect ratio however, here evaluated in a series of combinations with other parameters, was shown not to

Carbon Emission Reduction
To further demonstrate the benefits of the system, the carbon footprint reduction introduced by its use was evaluated using emission data from literature [43,44]. Assuming the system substitutes natural gas as a heating source, which has carbon emissions of approximately 510 tons of carbon dioxide equivalent (CO 2 eq) per GWh, the total amount of carbon footprint reduction was estimated. It was found that for the base properties, the proposed TES (with a volume of 1830 m 3 , corresponding to 5490 tons of rock) would be able to reduce emissions in about 156 tons of CO 2 eq per year, for all the aspect ratios. This is equivalent to 85 tons of CO 2 eq/m 3 of rock per year, or 28 tons of CO 2 eq/ton of rock per year. It is important to note that for a system replacing coal (∼905 tons of CO 2 eq per GWh) or heating oil (∼730 tons of CO 2 eq per GWh), the decrease in carbon footprint could be even more significant, these being displayed in Table 3. Table 3. Estimated carbon footprint reduction for heating using the proposed system (base properties).

Conclusions
In this study, a three-dimensional transient model was established to study the heat transfer and fluid flow in a thermal energy storage system consisting of a truncated cone-shaped bed of packed rocks. A considerable number of simulations were performed to account for the impact of all the distinct parameters (aspect ratio, HTF velocity, bed porosity, particle size, and thermophysical properties) and their combination in performance of the system. The study shows that porous media properties of the system, such as porosity, are extremely important in the performance of the system as a TES system, regarding both charging rate, energy stored/extracted, and pressure drop. The aspect ratio however, here evaluated in a series of combinations with other parameters, was shown not to have a significant effect in anything but the total pressure gradients across the system. Through the change in fluid velocity across the system, the pressure drop was also analyzed. It was shown that the pressure differential across the system quickly increases for increments in Reynolds number or a decrease in porosity. Particle size was also considered, and does not appear to significantly influence the heat transfer, having an effect mostly on permeability and, as a result, on pressure drop. A few different rock types were investigated, and it was seen that the most important thermophysical parameter for a rock bed, in terms of heat storage, is the thermal capacity. Some estimations of the environmental benefits of the system were also made, showing that it could lead to a carbon footprint reduction (from heating systems) of around 85-150 tons of CO 2 eq per cubic meter of rock per year.
Finally, such outcomes can be effectively implemented in the design and performance optimization of thermal energy storage systems based on packed rock material consisting of large particles, such as the ones found in caved zones. Using the domain with the most appropriate properties ought to bring significant savings in energy for a site implementing the technique, such as a remote (and cold) area, both commercial (such as mining operations) and non-commercial (such as first nation dwellings) in applications, such as heating for communal areas, domestic hot water provision, or even mine ventilation.