Thermal energy storage in combined cycle power plants: comparing finite volume to finite element methods

. The research in thermal energy storage (TES) systems has a long track record. However, there are several technical challenges that need to be overcome, to become omnipresent and reach their full potential. These include performance, physical size, weight and dynamic response. In many cases, it is also necessary to be able to achieve the foregoing at greater and greater scale, in terms of power and energy. One of the applications in which these challenges prevail is in the integration of a thermal energy storage with the gas turbine (GT) compressor inlet conditioning system in a combined cycle power plant. The system is intended to provide either GT cooling or heating, based on the operational strategy of the plant. As a contribution to tackle the preceding, this article describes a series of 3-dimensional (3D) numerical simulations, employing different Computational Fluid Dynamics (CFD) methods, to study the transient effects of inlet temperature and flow rate variation on the performance of an encapsulated TES with phase change materials (PCM). A sensitivity analysis is performed where the heat transfer fluid (HTF) temperature varies from -7 o C to 20 o C depending on the operating mode of the TES (charging or discharging). The flow rate ranges from 50% to 200% of the nominal inflow rate. Results show that all examined cases lead to instant thermal power above 100kWth. Moreover, increasing the flow rate leads to faster solidification and melting. The increment in each process depends on the driving temperature difference between the encapsulated PCM and the HTF inlet temperature. Lastly, the effect of the inlet temperature has a larger effect as compared to the mass flow rate on the efficiency of the heat transfer of the system.


Introduction
By the late 1950s, it was identified that to make best use of renewable energy resources with a meteorologically dependent output, a storage element to the overall system would increase the energy yield [1]. According to the public report of International Energy Agency (IEA), among the various storage systems, thermal energy storage (TES) is an enabling and crosscutting technology that can unlock potentials in various sectors [2].
For what concerns combined cycle (CC) power plants, the integration of TES acts as a counterpart to electrical storage and adds further flexibility in both thermal and electrical power generation. The capability of adding power dispatchability and overall power plant efficiency in parallel to the reduction of fossil fuels utilization can be identified as one of their most prominent advantages. Thus, they constitute a promising technology for the future energy systems.
An innovative concept where TES are integrated in real CC co-generative power plants is currently under research within the European Union's Horizon 2020 PUMP-HEAT [3]. The main objective of PUMP-HEAT is to develop and demonstrate replicable power plant layouts that are based on the integration of heat pumps (HP) and TES to un-tap the CC flexibility through low CAPEX (CAPital EXpenditure) balance of plant innovations. Additionally, the proposed innovations enable annual reduction of OPEX (OPerational EXpenditure) in the order of 3 to 5%, annually. One of the TES systems investigated in this project is applied with an integrated inlet conditioning system responsible for either lowering the inlet temperature of the GT (power augmentation) or heating it up (Minimum Environmental Load reduction), depending on the CC power plant operational strategy.
An extensive review work has been performed in [4], investigating the cooling methods used to provide GT power augmentation with particular focus on the merits of TES when combined with inlet cooling methods. Among the thermal energy storage systems, latent heat based phase change materials (PCM) show promising thermal energy storage density. Concerning the type of heat exchanger that is utilized in the selected storage system, several performance enhancement techniques have been studied in the available literature [5]. These include shell and tube configurations [6], encapsulated storages [7] and phase change slurries [8]. The materials used in these systems are organic PCMs which are carbon-based compounds, inorganic PCMs which are generally metallic and hydrated salts and eutectic PCMs of mixtures of two or more components [9], [10].
The available literature focuses on the comparison of the different mathematical models applied in packed bed storage systems. These consist of the continuous solid phase models, the Schumman's model, the single-phase models and the thermal diffusion models [11], [12]. Regarding the discretization methods, there is no other comparative study that utilizes different codes to simulate large-scale TES systems in a three-dimensional (3D) approach.
Having discussed the above, the aim of the present investigation is to numerically assess different simulation tools to study the performance of an encapsulated, packed bed TES. To achieve that end, this article describes the main differences of the selected methods (Finite Volume Method -FVM, Finite Element Method -FEM) that are utilized to solve the partial differential equations and predict the TES performance. The subsequent sections include the problem description, the methodology and lastly, the conclusions of this study. It must be noted that the scope of this work is not to evaluate which computational tool is more suitable for the research of such topics. It is rather a crosschecking effort of two independent studies performed by two institutes to evaluate the performance of the studied TES. Throughout the coordinated efforts, all of the varying parameters are compared aiming to assess the feasibility of the selected TES layout, prior to the experimental campaigns. Additionally, meshing solutions and overall approach of the problem is provided for the readers that anticipate researching on similar, computationally demanding applications. The overall goal is to serve as a contribution in the selection and design of a suitable cold storage tank that will fulfil the thermal power requirements for a validation site that is currently under construction at the University of Genoa, in Italy.

Discretisation Methods
Finite volume and finite element methods are numerical methods based on spatial and temporal discretization of mathematical model equations. The time is usually discretized with some time-stepping schemes for the ordinary differential equations. The FVM solver used in this study is the ANSYS Fluent 18.2 and the FEM solver is the COMSOL Multiphysics 5.4.
The FVM is a special finite difference formulation and its characteristic that distinguishes it from the FEM technique is the integration of the governing equations of fluid flow over the total number of the finite control volumes of the solution domain. The resulting statements express the exact conservation of relevant properties (mass, momentum, energy) for each cell. ANSYS Fluent typically uses the cell-centered finite volume approach which solves conservation equations at the cell centers and uses interpolation to express the variable values φ at the element centroid in terms of the nodal values at the control volume surface φf [13]. A beneficial feature of the method is that its base operation leads to local conservation of the net flux for each cell and is also providing a discretization stabilization for flows dominated by convection.
The FEM divides the computational domain into a finite number of subdomains, the socalled elements. The variables φ within each element are interpolated using local polynomials in terms of the values φj at a set of nodal points j (the corners of each element) in a way that continuity of the solution is guaranteed across the element sides [14]. A summary of this method can be described by three steps: Discretisation, interpolation and assembling of the equations. One distinct advantage of this method is that it allows to easily handle complex arbitrary geometries as it can be easily applied using irregular grids of various shapes. This also implies that it is a method providing the highest possible accuracy on coarse grids.

Phase Change Simulation
The FVM solver is based on the enthalpy-porosity technique used to model the solidification and melting processes. In this technique, the melt interface is not tracked explicitly. Instead, liquid fraction indicating the fraction of the cell volume that is in liquid form, is associated with each cell in the domain. The liquid fraction is computed at each iteration, based on an enthalpy balance. The mushy zone is a region in which the liquid fraction lies between 0 and 1. The mushy zone is modelled as a "pseudo" porous medium in which the porosity decreases from 1 to 0 as the material solidifies. When the material has fully solidified in a cell, the porosity becomes zero and hence the velocities also drop to zero [15].
On the other hand, the FEM solver is based on the enthalpy method, to track the moving boundary location and account for the phase change process. In this method, enthalpy is tracked as a function of temperature of the discretized elements by its ability to account for the latent heat lost/absorbed during the solidification/melting process.

Geometry
The cylindrical, vertically oriented TES tank is a packed bed configuration with encapsulated PCMs (Figure 1). The apparent advantages of the encapsulation are the increased heat transfer they provide and the impeding of mixing with the heat transfer fluid (HTF). The selected encapsulation morphology is an ellipsoid. This shape is characterized by a slightly increased heat transfer rate compared to the respective spherical capsule, due to its increased surface to volume ratio. Since the total volume of the tank is 5m 3 with a diameter of 1.9m., it is reasonable to down-scale the geometry to accelerate the 3D transient phase change simulation. Therefore, the model is reduced to a single module problem ( Figure 2) including 40 capsules in a staggered layout. The selection of this layout (instead of an aligned layout [16]) is deemed more realistic and is closed to a freely packed experimental setup where capsules are stacked inside the tank. At this point, it should be noted that the results of the downscaled single module are upscaled to the real geometry, by taking into consideration the contribution of the total module number which is derived from the available TES space.

Meshing
For the generation of the computational grid, two different mesh strategies are employed for each individual study. A coarse mesh is used with the FEM solver while a fine mesh is used with the FVM solver ( Figure 3). The coarse mesh amounts to approximately 260.000 cells with average mesh quality of 0.63 and cell size of 4.5 to 18 millimeters. As already stated earlier in this section, the FEM formulation provides the best accuracy for coarse grids. On the other hand, the finer mesh consists of an order of magnitude higher amount of cells (~2.5 million elements) with an average skewness of 0.35 (values close to 1 indicate poor quality).
The coarse mesh allows for faster solutions which can be evaluated and allow for model refinement. The fine mesh allows to capture the flow field in a more accurate way. Both methods attempt to predict the TES performance in an independent manner. This is achieved by examining the effect of the mass flow rate over the outlet temperature of the system through both methods. The goal is to evaluate the physical phenomena that are captured by each method and assess which model over or underpredicts the actual behavior of the TES.  additional grids were also generated. The first one with 5M elements and the second one with 1.25M elements. The selected parameter for the grid independence study was the outlet temperature which is one of the actual performance parameters that are evaluated within this study. The initial grid size was deemed appropriate (2.5M) and all the results presented in the following sections correspond to the 2.5M case. Moreover, the time-step size was defined as a function of the liquid fraction of the PCM. For the purposes of this study, three time step sizes were selected : 0.5, 1 and 2 seconds. Following the results shown in Figure 4, the time step size of 0.5 seconds was selected. The difference of the solution between the selected time steps is lower than 2%. In order to mesh the single-module in a staggered configuration with the capsules aligned in the way described above, some challenges occured in both grids. These are associated to highly skewed elements that are present at the junction of capsule shells (Figure 5a). To avoid these elements that introduce considerable numerical errors, a specific geometrical treatment is applied (Figure 5b) that only introduces a minor systematic error, thus resulting in a discretized domain consisting of acceptable elements in terms of orthogonality and cell skewness. Moreover, this approach leads to acceptable solution convergence meeting the chosen criteria of continuity <10 -4 , momentum <10 -6 , energy <10 -6 .

Initial conditions and specific heat capacity (cp) treatment
The performance for both charging and discharging operations of the TES are evaluated herein, with mass flow rates of 50% nominal, nominal and 200% nominal conditions. The different values of mass flow examined in this sensitivity analysis are selected for the expected operating conditions. The PCM thermophysical properties were obtained from T-History methods of a proposed PCM in the PUMP-HEAT project [3] with peak phase change temperatures at around 5 o C.
This part describes how the PCM specific heat capacity (cp) is pre-treated in the modelling part. In order to integrate the cp values into the commercial software, particular data treatment is required. A normal (Gaussian) distribution is applied to the specific heat to smoothen out negative values due to supercooling and reduce superseding error from the experimental measurements. The resulting cp distributions for the two processes (melting and freezing) are shown in Figure 6.

Figure 6. cp graphs for cooling (charging) and heating (discharging)
In general, most materials do not melt at a constant temperature, but over a temperature range due to impurities and because of blending composition, such as paraffins and salt hydrates [17]. Though the range is small, it may significantly influence the performance of the system. In this investigation, the melting of the PCM selected is taken as Tsol=3.3 o C and Tliq=6.3 o C. The selection of the aforementioned temperatures is based both on the experimental acquisition of the PCM properties which leads to the derived cp graphs and also, on the observed hysteresis which amounts to approximately 1 o C. This discrepancy is expected to lead to a difference on the charge/discharge capacity if the operating temperature is not wide enough to cover the full freezing/melting range. As a result, it is expected that this would have an effect on the round-trip efficiency (ratio of energy stored to energy retrieved). Other assumptions are that the PCM is homogeneous and isotropic, freezing/melting occur congruently, buoyancy force is negligible, flow is laminar, the thermal conductivity and density of PCM are constant in solid and liquid phases and capsule thermal resistance is negligible.
The following cases are tested, utilizing both FEM and FVM numerical models and employing the initial conditions shown in Table 1. Two different discharging modes and one charging mode are examined based on the different operations of the power plant and hence, the different requirements from the TES system. The two discharging modes correspond to the minimum and maximum power required by the tank. The HTF selected is an aqueous solution of ethylene glycol (35% wt.).

Results and discussion
The results presented below describe the modeling process followed that predicts an estimated behavior of the storage tank, prior to the final design and validation of the calculations with experimental testing. The results of the two models provide an approximation on the outlet temperature of the TES tank as a function of the HTF mass flow, through an ordinary least squares (OLS) estimation method. Throughout this approach, a simple regression line was fitted by the minimization of the squared sum of the calculated residuals. The results that are presented below will be readjusted based on experimental results. The case depicted in Figure 7 represents the TES outlet temperature for the minimum power discharging case. This leads to thermal power estimation of the examined systems, as a function of time. The final results that emerge out of this analysis are presented in the following figures. The first case study (Figure 7) considers the temperature profile of the ethylene-glycol HTF for the discharging case at the time when ambient conditions impose the minimum temperature difference scenario (minimum thermal power). The temperature profile of the model shows a phase change temperature plateau in all the three cases. At 50% of the nominal flow rate (mdot/2) the velocities of the fluid inside the tank are low. It is observed that the HTF outlet temperature is nearly constant (within 1 o C variation) for more than thirty minutes of operation which is a desirable aspect of the application as this favors the stabilization of the GT inlet temperature. The second studied case considers the maximum power. The higher PCM-HTF temperature difference results in higher rates of heat transfer between the PCM and HTF. Once again, the lowest mass flow rate (mdot/2) results in a longer phase change in terms of temporal distribution and at a more constant temperature range ( Figure 8).
Similar results are obtained in the charging case (Figure 9) when the TES acts Figure 9. Numerically predicted outlet temperature distribution during charging for the examined mass flow rates as the heat source and the PCM solidifies. The range of the latent heat transfer region in this case is smaller as compared to the two discharging cases. The difference is associated with the specific heat capacity (cp) integration in the models.
Moreover, Figure 10 shows the liquid fraction of the packed bed as a function of time during charging and discharging, for different mass flow rates. As expected, the time for complete melting and freezing decreases with increasing mass flow rates due to increased velocities and hence the increased heat transfer rate from the capsules. It is also observed that the complete solidification time is slightly longer than the melting time. This is because of the overall heat transfer coefficient which is lower during solidification than the melting process and due to subdued natural convection in freezing. Additionally, the full phase change is realized within two hours with all three HTF flow rates.

Figure 10. Temporal variation of liquid fraction of the packed bed during discharge (max_power) and charge processes
The inlet temperature has a larger effect as compared to the mass flow rate on the heat transfer rate of the system. This can be observed in Figure 11 showing the liquid fraction of the packed bed system during the discharging process of the minimum power scenario. The complete melting requires one more hour of discharging as compared to the maximum power case.
The thermal power during all processes are shown in Figure 12. The upper two diagrams correspond to discharging and the lower diagram to the charging case respectively. The charged/extracted thermal power rate is calculated with the change in temperature of the HTF from inlet to outlet, with the following equation: The accumulated charged and discharged energy capacity is calculated by the following equation within the tn limit of integration: The results show that the thermal power that can be attained through the studied configurations is more than 100kWth. The peak values are as expected at the beginning of Figure 11. Temporal variation of liquid fraction of the packed bed during minimum power discharge process each process when the highest temperature difference exists between the inlet and the storage which is however reduced progressively with time. These values are then flattened out near the end when reaching full charge/discharge capacity. As far as the maximum power scenario is considered, it is observed that more than 75kWth can be achieved even for the lowest mass flow rates.  Figure 13 demonstrates the amount of energy capacity discharged and stored during the two operations respectively. After two hours of operation when the PCM has completely melted or solidified, the storage system cannot continue its main operation as its energy reserve is depleted. It can also be observed that the energy is not linearly proportional to the mass flow rate.
An additional comment that can be derived based on this analysis at this point considers the round trip efficiency of the TES system. Based on the discussion earlier in this works considering the PCM properties and the observed hysteresis, a difference on the charge/discharge capacity is expected if the operating temperature is not wide enough to cover the full freezing/melting range. The comparison of the charging and discharging energy capacities cannot be performed in an absolute manner due to the different operating conditions studied for each mode. For instance, the difference between HTF inlet temperature and PCM phase change temperature is 3 o C at minimum power discharging case while the same absolute difference is 12.27 o C at the single charging case studied. However, it is estimated that the round trip efficiency in these systems will be higher than 80%.

Conclusions
This work investigates the transient effects of inlet temperature and heat transfer fluid flow rate variation, on the performance of an encapsulated thermal energy storage employing phase change materials. Different computational softwares have been utilized, each based on independent simulation approaches.
Both finite volume and finite element methods have given comparable results to the performed studies, validating thus the obtained conclusions which are outlined as follows: • The lowest mass flow rates result in a nearly constant outlet temperature for 15 to 30 minutes of operation. • All examined cases lead to instant thermal power above 100kWth and an average thermal power of more than 45kW. • The charging and discharging time have been shown to reach to approximately 2 hours of operation. Increasing the flow rate leads to faster solidification and melting. The increment in each process depends on the driving temperature difference between the encapsulated PCM and the HTF inlet temperature • The inlet temperature has a larger effect as compared to the mass flow rate on the heat transfer performance of the system