Optimization of Design Variables of a Phase Change Material Storage Tank and Comparison of a 2D Implicit vs. 2D Explicit Model

: In this study, a thermal energy storage tank ﬁlled with commercial phase change material ﬂat slabs is investigated. The tank provides heat at around 15 ◦ C to the evaporator of a seasonal thermal energy storage system developed under the EU-funded project SWS-Heating. A 2D numerical model of the phase changed material storage tank based on the ﬁnite control volume approach was developed and validated with experimental data. Based on the validated model, an optimization was performed to identify the number, type and conﬁguration of slabs. The ﬁnal goal of the phase change material tank model is to be implemented into the whole generic heating system model. A trade-off between results’ accuracy and computational time of the phase change material model is needed. Therefore, a comparison between a 2D implicit and 2D explicit scheme of the model was performed. The results showed that using an explicit scheme instead of an implicit scheme with a reasonable number of nodes (15 to 25) in the heat transfer ﬂuid direction allowed a considerable decrease in the computational time (7 times for the best case) with only a slight reduction in the accuracy in terms on mean average percentage error (0.44%). Author Contributions: Conceptualization, A.C., G.Z., A.d.G.; methodology, A.C., G.Z., A.d.G., C.F.; software, A.C.; data curation, L.F.C.; validation, A.C., D.V., G.Z., A.d.G.; formal analysis, A.C., G.Z., A.d.G.; investigation, A.C., E.B., D.V.; resources, L.F.C.; writing—original


Introduction
Thermal energy storage (TES) allows the supply of thermal demand to be independent of the heat source, especially when the heat source is discontinuous, such as solar energy. Among all the TES alternatives, latent heat TES with the use of phase change materials (PCM) stands out for its high energy density and ability to provide heat at a nearly constant temperature during the phase change. PCM storage units were analyzed for very different applications, such as buildings [1,2], process heat [3], concentrated solar power (CSP) plants [4] or solar cooling [5]. The shape of the encapsulated PCM and the melting temperature of the PCM are essential variables in the design of the systems. Focusing only on HVAC applications (heating, ventilation, air conditioning), the use of flat slabs with rectangular shapes was frequently investigated in the literature, due to its simplicity of operation and its commercialization.
The use of PCM as storage material is usually complicated due to certain limitations, such as low thermal conductivity, subcooling, phase segregation, non-uniform distribution of the PCM inside the slabs or capsules, or the non-uniform composition of the PCM. Although some of these effects are very difficult to be considered in a theoretical model, a numerical simulation model can be very useful to analyze its performance before designing or experimentally testing a PCM storage unit, especially when it is coupled to another system. Several authors presented numerical models where rectangular flat slabs were of TRNSYS were compared. Crank-Nicholson was the method selected by the author because it is unconditionally stable and allowed the model to use larger time steps than other methods (Euler method). Moreover, it provided accurate results with the fewest number of iterations. Jutang et al. [16] compared four methods (explicit Euler method, implicit Euler method, implicit midpoint rule and Runge Kutta Fourth-order) with the exact solution to solve the equations of a water tank model. Focusing on the comparison between the explicit and the implicit Euler method, the implicit Euler method provided higher accuracy. Moreover, when the time step of the explicit Euler was too high, the solution became unstable. That means that the explicit method requires a small time-step [17].
In the previously reviewed literature, several authors validated numerical models of PCM storage units [2,7,13] and studied the impact of geometry (thickness of PCM, thickness of gap, etc.) and/or operational parameters (mass flow, inlet temperature, etc.) [6,11] on the results' accuracy. Furthermore, some authors compared different mathematical methods to solve the energy balance equations of TES models. However, to the best knowledge of authors, no studies have compared the uncertainties between an implicit and explicit model of a PCM storage unit and the trade-off between results accuracy and computational time. Especially when the optimization of a whole heating system or its control is performed, it is essential to identify the numerical model that minimizes the computational time while keeping good results accuracy.
In this paper, some design variables (type of slab, number of PCM slabs and their configuration) of a PCM storage unit that supplies heat to the evaporator of a seasonal TES system, based on selective water sorbent material, were optimized. For this purpose, a 2D implicit numerical model was developed and validated against experimental data of a charging process. Furthermore, to implement the PCM storage unit model into the simulation of the whole heating system, the best compromise between accuracy and computational time is needed. Therefore, an uncertainty analysis comparison between a 2D implicit and 2D explicit model was carried out.

Numerical Model (2D)
A numerical model was developed to simulate the charge and discharge process of the PCM storage unit. The modeled PCM storage unit (see Figure 1) consists of a rectangular tank filled with water and horizontal PCM slabs stacked one above another with a gap between them to allow the HTF flow through the tank. At both the inlet and outlet of the tank, there is a column of water, in which a diffuser is arranged to distribute the flow of water through the entire tank. The numerical model allowed to simulate the heat transfer in the PCM slabs, in the HTF flowing inside the channels between slabs, and in the two columns of water. The numerical model of the PCM tank was developed in Python [18]. The following assumptions were made: • The density and thermal conductivity of the PCM were considered constant in each of the phases and also during the phase change.

•
The phase change was taken into account through an equivalent specific heat capac- The following assumptions were made: • The density and thermal conductivity of the PCM were considered constant in each of the phases and also during the phase change.

•
The phase change was taken into account through an equivalent specific heat capacity [19]. During the phase change, the temperature dependence on the PCM specific heat capacity had a trapezium shape centred on the melting temperature of the PCM with an upper base of the trapezium of 1 • C. It was an analogous methodology to the one presented by Farid [20], who used a triangle shape distribution, which was found to be successful for modelling heat transfer in PCM.

•
The PCM was considered homogeneous and isotropic. • All the PCM was initially at the same temperature. The initial temperature of the HTF inside the tank and of the plastic container has been also set to the same temperature. • Thermal losses to the environment were not considered. • All the slabs had symmetrical thermal behaviour with respect to the horizontal plane that crosses their centre.

•
The mass flow was divided evenly between each HTF channel. • Natural convection was neglected inside the liquid phase of the PCM due to the thin, flat nature of the slabs and their horizontal layout [21]. • Constant C p of the HTF was considered since according to Zsembinszki et al. [6] the density and specific heat of the HTF could be considered constant to reduce the computational time as they have less importance in the mean absolute percentage error (MAPE).
The convective heat transfer coefficient h HTF between the HTF and the PCM slab was calculated using the correlation for the average Nusselt number presented by Kakac et al. [22] for laminar flow in a rectangular duct, given the Reynolds number and the Prandtl number at the HTF bulk temperature, the ratio of the tube length to hydraulic diameter and the ratio of the thickness of the duct to its width. The correlation can be solved for a constant temperature wall or for a constant heat flux. For this study, the correlation for a constant temperature wall was considered, since as Zsembinszki et al. [6] pointed out, that condition is more suitable when a phase change occurs. At every time-step, an average convective heat transfer coefficient h HTF based on the assumed bulk HTF properties was used. Avoiding the calculation of a different h HTF for each node results in a reduction of the computational time with very little effect on the results [6].
Considering these assumptions, a fully implicit finite control volume method in the two-dimensional Cartesian coordinates system [23] was used for solving the energy balance equations. The iterative method Gauss-Seidel was used to solve the set of equations. As it can be seen in Figure 2, a two-dimensional heat transfer in the x and y direction (Q x and Q y ) was considered inside of each PCM node (T PCM (x,y)) and each container wall node (T w (x,y)), and one-dimensional heat transfer was considered for the HTF (T HTF (x)). Assuming symmetry, half of the water channel and half of the PCM slab were simulated. Heat transfer in the PCM and container wall is governed by the heat conduction equation: Heat transfer during the circulation of HTF through the channels is governed by the following equation: The set of governing Equations (1) and (2) is completed by the following initial and boundary conditions ( Figure 2): The set of governing Equations (1) and (2) is completed by the following initial and boundary conditions ( = 0 = 0 due to symmetry (8) = + + = 0 due to the symmetry (10) Finally, to calculate the instantaneous heat transfer rate (HTR) during the charging process of the whole PCM tank, the power released by the HTF was used (Equation (11)), understood as a temperature gradient between inlet and outlet temperature of the HTF. It is important to highlight that this HTR was not just the latent energy stored by the PCM slabs, but also the sensible energy stored by the PCM and the HTF:

Experimental Set-up
The experimental set-up used to validate the numerical model was located in the laboratory of GREiA research group at the University of Lleida (Spain). The main component of the set-up consisted of a rectangular tank (Figure 3a), as explained in the previous section, filled with water and 10 horizontal PCM slabs stacked one above another with a Finally, to calculate the instantaneous heat transfer rate (HTR) during the charging process of the whole PCM tank, the power released by the HTF was used (Equation (11)), understood as a temperature gradient between inlet and outlet temperature of the HTF. It is important to highlight that this HTR was not just the latent energy stored by the PCM slabs, but also the sensible energy stored by the PCM and the HTF: m tank C p,HTF (T in,HTF − T out,HTF )

Experimental Set-Up
The experimental set-up used to validate the numerical model was located in the laboratory of GREiA research group at the University of Lleida (Spain). The main component of the set-up consisted of a rectangular tank (Figure 3a), as explained in the previous section, filled with water and 10 horizontal PCM slabs stacked one above another with a gap of 10 mm between them to allow the HTF flow through the tank. The slabs used for the experiments were commercial Flat-ICE slabs (Figure 3b) provided by PCM Products Ltd. containing PlusICE 15 [24] as PCM, which is a hydrated salt. Each Flat-ICE slab had dimensions of 0.5 m × 0.25 m × 0.035 m. The container of the PCM slabs was rigid and made of HDPE 2.5 mm thick. The dimensions L 1 , L 2 and L 3 (see Figure 1) of the tank are respectively 0.06 m, 0.56 m and 0.062 m. Pure water was used as HTF. The desired inlet temperature of the HTF was obtained by a compression cooling unit (Zanotti model GCU2030ED01B) of 5 kW of maximum cooling power, and a thermal bath with two immersed thermostats (OVAN TH100E and JP SELECTA) of 3 kW of total thermal power. Nine calibrated Pt-100 class B temperatures sensors with accuracy ±0.3 • C were attached to the outer part of three PCM slabs at different levels: 3 at the top, 3 at the middle and 3 at  Figure 4) to analyze the PCM phase change process and estimate when the phase change is occurring. Moreover, to measure the power delivered by the HTF, a flow meter and two temperature sensors were placed at the inlet and outlet of the HTF circuit. A data acquisition system was used to collect all monitored variables at 10 s time interval.
made of HDPE 2.5 mm thick. The dimensions L1, L2 and L3 (see Figure 1) of the tank are respectively 0.06 m, 0.56 m and 0.062 m. Pure water was used as HTF. The desired inlet temperature of the HTF was obtained by a compression cooling unit (Zanotti model GCU2030ED01B) of 5 kW of maximum cooling power, and a thermal bath with two immersed thermostats (OVAN TH100E and JP SELECTA) of 3 kW of total thermal power. Nine calibrated Pt-100 class B temperatures sensors with accuracy ±0.3 °C were attached to the outer part of three PCM slabs at different levels: 3 at the top, 3 at the middle and 3 at the bottom (see Figure 4) to analyze the PCM phase change process and estimate when the phase change is occurring. Moreover, to measure the power delivered by the HTF, a flow meter and two temperature sensors were placed at the inlet and outlet of the HTF circuit. A data acquisition system was used to collect all monitored variables at 10 s time interval.
(a) (b)  The physical properties of the PCM and plastic container are shown in Tables 1 and  2, respectively. For this case study, a temperature range of four degrees Celsius (lower base of trapezium) during the phase change was considered. The value of latent heat presented in Table 1 was obtained as the average value of the latent temperature range (±2 °C from the melting temperature) given by the manufacturer. The same methodology was used for the calculation of the solid and liquid specific heat capacity average values in their corresponding temperature range (solid and liquid sensible range). The experiment used to validate the model consisted of a charging process at a total mass flow rate of 240 kg/h at an HTF inlet temperature around 25 °C. To perform the charging process, HTF made of HDPE 2.5 mm thick. The dimensions L1, L2 and L3 (see Figure 1) of the tank are respectively 0.06 m, 0.56 m and 0.062 m. Pure water was used as HTF. The desired inlet temperature of the HTF was obtained by a compression cooling unit (Zanotti model GCU2030ED01B) of 5 kW of maximum cooling power, and a thermal bath with two immersed thermostats (OVAN TH100E and JP SELECTA) of 3 kW of total thermal power. Nine calibrated Pt-100 class B temperatures sensors with accuracy ±0.3 °C were attached to the outer part of three PCM slabs at different levels: 3 at the top, 3 at the middle and 3 at the bottom (see Figure 4) to analyze the PCM phase change process and estimate when the phase change is occurring. Moreover, to measure the power delivered by the HTF, a flow meter and two temperature sensors were placed at the inlet and outlet of the HTF circuit. A data acquisition system was used to collect all monitored variables at 10 s time interval.
(a) (b)  The physical properties of the PCM and plastic container are shown in Tables 1 and  2, respectively. For this case study, a temperature range of four degrees Celsius (lower base of trapezium) during the phase change was considered. The value of latent heat presented in Table 1 was obtained as the average value of the latent temperature range (±2 °C from the melting temperature) given by the manufacturer. The same methodology was used for the calculation of the solid and liquid specific heat capacity average values in their corresponding temperature range (solid and liquid sensible range). The experiment used to validate the model consisted of a charging process at a total mass flow rate of 240 kg/h at an HTF inlet temperature around 25 °C. To perform the charging process, HTF The physical properties of the PCM and plastic container are shown in Tables 1 and 2, respectively. For this case study, a temperature range of four degrees Celsius (lower base of trapezium) during the phase change was considered. The value of latent heat presented in Table 1 was obtained as the average value of the latent temperature range (±2 • C from the melting temperature) given by the manufacturer. The same methodology was used for the calculation of the solid and liquid specific heat capacity average values in their corresponding temperature range (solid and liquid sensible range). The experiment used to validate the model consisted of a charging process at a total mass flow rate of 240 kg/h at an HTF inlet temperature around 25 • C. To perform the charging process, HTF was circulated through the PCM tank until the average temperature of all the sensors inside the tank reached a temperature of 24 • C. The temperature of the entire PCM tank at the beginning of the charging experiment was 5 • C. Table 1. Thermophysical properties of PlusICE S15 for the melting process [24].

Model Validation
The accuracy of the numerical model and its validation was evaluated with experimental data of a charging process of the PCM tank, presented in Section 2.2. Since the initial temperature in the tank was 5 • C and the experiment was stopped when the average temperature of all PCM slabs temperature sensors reached 24 • C, a temperature gradient of 20 • C between initial and final PCM state was analyzed.
Since the experimental data were collected every 10 s, this time step was also considered for the computational simulation. Moreover, for the validation of the numerical model, the following assumptions were considered: 25 nodes in the HTF direction, 4 nodes in the plastic wall, and 9 nodes in the PCM slab. Both the number of nodes and the time-step showed grid independence.
To assess the uncertainty of the measurements used to validate the model, a sensitivity analysis of the inlet and outlet temperature and heat transfer rate was performed. The uncertainty of the measured temperatures was calculated using the Equation (12). The uncertainty of the heat transfer rate was calculated applying the Equation (13) [25] to the density and specific heat. The uncertainty of the HTF thermophysical properties and heat transfer rate was estimated at each registered time step, and then the mean value was computed. The uncertainty values obtained were: ±0.201 • C for the inlet temperature, ±0.199 • C for the outlet temperature and ±0.054 kW for the heat transfer rate. Figure 5 presents the outlet temperature and HTR of the HTF (Equation (11)), respectively, for both experimental and numerical results. As shown in the figures, there is a good match between experimental and simulated profiles for both the HTF outlet temperature and HTR. In the early period of melting, the model predicted slightly higher values of HTF than those obtained from the experiment. The temperature of PCM slabs and HTF were initially at around 5 • C. However, the vertical column of HTF at the outlet of the tank could still have some nodes at a slightly higher temperature (just 3 sensors were placed in the vertical column of water). This may be the reason for that disagreement. On the other hand, at the end of the experiment, the model predicted lower values of heat transfer rate. This is probably because at that stage, the PCM temperatures were above ambient temperature, and the PCM tank was absorbing a small amount of energy from the ambient. Neglecting the early period of the charging process, when the error was considerably higher, the error between the numerical model results and the experimental values was ±1%. pose.

Design Optimization of the System for Real-Scale Application
A 2D implicit numerical model of a PCM tank was developed and validated to be used to design the optimal configuration of the PCM tank for the particular application of the generic heating system of the EU-funded project SWS-Heating (GA 764025). During its operation, the PCM tank is first charged by a solar thermal field to store heat at lowgrade temperature (around 15 °C). Afterwards, the PCM tank provides the required heat to the evaporator of the seasonal sorption tank during its discharge ( Figure 6). Therefore, the PCM tank must be able to provide the heat power required by the evaporator. To obtain the optimal configuration of the theoretical PCM tank, the number of slabs per column and number of columns in series (meaning that the total length of the TES would be the slab length multiplied by the number of columns in series) were optimized based on the discharging process of the PCM tank. Moreover, two different types of PCM slabs containing the same PCM but with two different thickness were analyzed in the optimization: the Flat-ICE slab [24], which was already used in the experimental study, and the Thin-ICE slab [24]. The physical properties of the PCM slabs and the parameters used in the optimization are presented in Table 3.
The validity of the optimal configuration of the PCM tank was analyzed in terms of the mean absolute percentage error (MAPE), which provides information about the relative error. In our case, the interest lies in the evaluation of the MAPE based on the delivered heat power of the tank ( ) with respect to the desired delivered heat power to the evaporator (see Equation (14)). The application of this PCM tank is providing heat to the evaporator of a seasonal sorption TES, which is part of a novel complex generic heating system that provides space heating and domestic hot water for residential applications. The mathematical model will be used to analyze the performance of the whole system and optimize its control. Hence, the numerical model presented in this study is considered enough accurate for this purpose.

Design Optimization of the System for Real-Scale Application
A 2D implicit numerical model of a PCM tank was developed and validated to be used to design the optimal configuration of the PCM tank for the particular application of Energies 2021, 14, 2605 9 of 15 the generic heating system of the EU-funded project SWS-Heating (GA 764025). During its operation, the PCM tank is first charged by a solar thermal field to store heat at low-grade temperature (around 15 • C). Afterwards, the PCM tank provides the required heat to the evaporator of the seasonal sorption tank during its discharge ( Figure 6). Therefore, the PCM tank must be able to provide the heat power required by the evaporator. which to search must be defined. The objective function consisted of minimizing the mean absolute percentage error (MAPE) based on HTR. Due to the high price of PCM and the requirements of the system, a maximum PCM mass between 180 and 220 kg was considered. To study all the possible combinations in the optimization space that fits the mass restriction, the range values for the number of slabs per column and number of columns in series were selected (shown in Table 3). The combinations that did not fulfill the mass restriction were discarded by the optimization algorithm.
The evaporator of the sorption tank allows inlet HTF coming from the PCM tank at a temperature range between 5 and 15 °C, at a mass flow rate of 600 kg/h. For the optimization of the PCM tank design, the HTF inlet temperature to the PCM tank (equal to the HTF outlet temperature from the evaporator) was set at 9 °C. The nominal heat power required by the evaporator ( ) used in the calculation of the objective function (Equation (14)) was 2.1 kW. However, this heat power demand can vary between 1 and 3 kW during the operation of the seasonal sorption tank, because its value depends on the inlet HTF coming from the PCM tank to the evaporators and on the HTF temperature coming from the water tank of the system ( Figure 6). Hence, to avoid optimization results where MAPE was minimized, but the operation power was out of the desired range from 1 to 3 kW, an additional constrain in the optimization was considered: at least half of the time of the discharging process, the HTR must be between ±1 kW with respect to the average desired HTR. Furthermore, to speed up the optimization time, a constant Nusselt number to calculate the convective heat transfer coefficient was used. The constant Nusselt value (7.54) provided by Çengel [29] for a fully developed laminar flow in tubes with a rectangular section with an infinite value of the ratio width/thickness (valid for both Flat-ICE and Thin-ICE) was used. This constant value presented a slight difference (error of 2.5% with respect to the Flat-ICE) compared to the correlation presented by Kakac et al. [22] calculated at every time-step.

Comparison between 2D Implicit and 2D Explicit Model
The computational time of the simulation of the discharge process of the PCM tank using the 2D implicit model with a 10 s time-step lasts more than 300 s due to the high number of nodes and iterative nature of the implicit approach.
In future studies, the implementation of smart control policies to optimize the control of the SWS-Heating system will be analyzed. If smart control algorithms are based on reinforcement learning, a learning process of the controller is required. The learning process is an iterative process whose duration depends on several factors, such as the com- To obtain the optimal configuration of the theoretical PCM tank, the number of slabs per column and number of columns in series (meaning that the total length of the TES would be the slab length multiplied by the number of columns in series) were optimized based on the discharging process of the PCM tank. Moreover, two different types of PCM slabs containing the same PCM but with two different thickness were analyzed in the optimization: the Flat-ICE slab [24], which was already used in the experimental study, and the Thin-ICE slab [24]. The physical properties of the PCM slabs and the parameters used in the optimization are presented in Table 3. The validity of the optimal configuration of the PCM tank was analyzed in terms of the mean absolute percentage error (MAPE), which provides information about the relative error. In our case, the interest lies in the evaluation of the MAPE based on the delivered heat power of the tank (Q sim ) with respect to the desired delivered heat power to the evaporator Q re f (see Equation (14)).
The optimization of the parameters was carried out using the library Hyperopt [26], which uses different stochastic search algorithms to find the best scalar value into a space. In this study, the Tree-structured Parzen Estimators (TPE) algorithm [27,28] was chosen. For the hyperparametric optimization, a maximum of 180 iterations were set, which means that on each iteration a new combination of the number of slabs and number of columns in series was studied. Within the TPE algorithm, a sample of 1000 candidates was set, and the first 20 iterations were considered randomly. Gamma value of 20% was used, which means that 20% of the best observations are used to estimate the next set of parameters. To perform the optimization, the objective function to minimize and the space over which to search must be defined. The objective function consisted of minimizing the mean absolute percentage error (MAPE) based on HTR.
Due to the high price of PCM and the requirements of the system, a maximum PCM mass between 180 and 220 kg was considered. To study all the possible combinations in the optimization space that fits the mass restriction, the range values for the number of slabs per column and number of columns in series were selected (shown in Table 3). The combinations that did not fulfill the mass restriction were discarded by the optimization algorithm.
The evaporator of the sorption tank allows inlet HTF coming from the PCM tank at a temperature range between 5 and 15 • C, at a mass flow rate of 600 kg/h. For the optimization of the PCM tank design, the HTF inlet temperature to the PCM tank (equal to the HTF outlet temperature from the evaporator) was set at 9 • C. The nominal heat power required by the evaporator (Q t re f ) used in the calculation of the objective function (Equation (14)) was 2.1 kW. However, this heat power demand can vary between 1 and 3 kW during the operation of the seasonal sorption tank, because its value depends on the inlet HTF coming from the PCM tank to the evaporators and on the HTF temperature coming from the water tank of the system ( Figure 6). Hence, to avoid optimization results where MAPE was minimized, but the operation power was out of the desired range from 1 to 3 kW, an additional constrain in the optimization was considered: at least half of the time of the discharging process, the HTR must be between ±1 kW with respect to the average desired HTR.
Furthermore, to speed up the optimization time, a constant Nusselt number to calculate the convective heat transfer coefficient was used. The constant Nusselt value (7.54) provided by Çengel [29] for a fully developed laminar flow in tubes with a rectangular section with an infinite value of the ratio width/thickness (valid for both Flat-ICE and Thin-ICE) was used. This constant value presented a slight difference (error of 2.5% with respect to the Flat-ICE) compared to the correlation presented by Kakac et al. [22] calculated at every time-step.

Comparison between 2D Implicit and 2D Explicit Model
The computational time of the simulation of the discharge process of the PCM tank using the 2D implicit model with a 10 s time-step lasts more than 300 s due to the high number of nodes and iterative nature of the implicit approach.
In future studies, the implementation of smart control policies to optimize the control of the SWS-Heating system will be analyzed. If smart control algorithms are based on reinforcement learning, a learning process of the controller is required. The learning process is an iterative process whose duration depends on several factors, such as the computational time of the system simulation (in this case the whole heating system), the control time slot (i.e., 10 or 15 min) or the duration of every test (i.e., 12 days or 1 whole year). During the learning process, thousands of simulations are frequently required to identify the optimal control. To reduce the computational time during the control optimization of an energy system, different solutions were presented in the literature. For instance, Romani et al. [30] used an explicit model, and de Gracia et al. [31] increased the control time slot from 10 to 15 min to reduce the computational effort.
The simulation time reduction of every element of the whole system is essential to minimize the learning process of its smart controller. For this reason, the computational effort of the whole simulation, and therefore of the PCM tank, in each time-step, must be reduced to the maximum extent, while keeping a considerable accuracy in the model results. In this section, with the goal to identify the trade-off between reliable model and low computational time, a 2D explicit model was developed and compared against the 2D implicit model validated in Section 2. nodes in the x-axis (N X ), wall nodes and PCM nodes were also analyzed. The 2D explicit model was based on the same equations presented in Section 2.3, but solved in an explicit scheme. For this reason, a low simulation time-step must be used to allow the convergence of the model. For the convergence of the 2D explicit model, a maximum time-step of 1 s was considered.
The simulations were done using a computer with an 8 GB RAM and an AMD FX(tm)-6100 Six-Core Processor 3.30 GHz.

Optimization of Design Variables of the PCM Tank
During the first period of the discharge process, due to the thermal inertia of the HTF column at the entrance and exit of the tank, there is a high temperature gradient between the inlet and the outlet temperature of the HTF. Hence, the heat transfer rate during the first minutes of the discharge is very high compared to the rest of the process. For this reason, the first 15 min of the discharge process were neglected for the calculation of the objective function, MAPE.
Two optimizations processes were launched, one for the Flat-ICE and one for the Thin-ICE. The minimization of the objective function was obtained for a tank composed of 33 slabs of Flat-ICE (MAPE = 31.26% with respect to 2.1 kW). The optimization for a PCM storage tank composed of Thin-ICE slabs could not find any value, which minimized the MAPE and fulfilled the two proposed constraints (range of heat power and mass limitation). For comparison with the optimized Flat-ICE configuration, a second optimization for the Thin-ICE slabs was launched. This time, the heat power constraint was neglected to let the optimization provide an optimum result. The optimal configuration (MAPE = 39%) was found for two columns of PCM slabs in series with 30 Thin-ICE slabs inside each tank. The optimized configurations for both types of slabs are shown in Figure 7.  Figure 8 presents the results of the optimal PCM tank configuration, which co sponds to the Flat-ICE and its comparison with the Thin-ICE. Since Thin-ICE has a thin thickness and less mass per slabs compared with the Flat-ICE, which means that for same tank mass, there is more heat transfer area between the PCM and the HTF, the T ICE is able to release the heat at higher power with the same amount of PCM. Never less, for the studied application, where heat power is required between 1 and 3 kW, Flat-ICE fulfilled better the requirements. The energy stored in the Flat-ICE slabs was leased at a lower heat power for a longer time compared to the Thin-ICE.   Figure 8 presents the results of the optimal PCM tank configuration, which corresponds to the Flat-ICE and its comparison with the Thin-ICE. Since Thin-ICE has a thinner thickness and less mass per slabs compared with the Flat-ICE, which means that for the same tank mass, there is more heat transfer area between the PCM and the HTF, the Thin-ICE is able to release the heat at higher power with the same amount of PCM. Nevertheless, for the studied application, where heat power is required between 1 and 3 kW, the Flat-ICE fulfilled better the requirements. The energy stored in the Flat-ICE slabs was released at a lower heat power for a longer time compared to the Thin-ICE. Figure 8 presents the results of the optimal PCM tank configuration, which corresponds to the Flat-ICE and its comparison with the Thin-ICE. Since Thin-ICE has a thinner thickness and less mass per slabs compared with the Flat-ICE, which means that for the same tank mass, there is more heat transfer area between the PCM and the HTF, the Thin-ICE is able to release the heat at higher power with the same amount of PCM. Nevertheless, for the studied application, where heat power is required between 1 and 3 kW, the Flat-ICE fulfilled better the requirements. The energy stored in the Flat-ICE slabs was released at a lower heat power for a longer time compared to the Thin-ICE.

Comparison of Implicit and Explicit 2D Schemes during Discharge
To compare the 2D implicit and 2D explicit schemes in the developed numerical model, the discharge process of the optimized PCM storage unit composed of 33 Flat-ICE slabs was used. The implicit PCM model with 25 nodes in the x-axis (Nx), four nodes of wall (N_wall) and nine nodes of PCM (N_PCM) with a time-step of 10 s was considered as reference (the same parameters as the validated model). To demonstrate the time-step grid independence of the model, the results for the same number of nodes with a timestep of 1 s were calculated, obtaining a MAPE of 0.3% with respect to the reference. The results of both implicit and explicit schemes can be seen in Table 4.

Comparison of Implicit and Explicit 2D Schemes during Discharge
To compare the 2D implicit and 2D explicit schemes in the developed numerical model, the discharge process of the optimized PCM storage unit composed of 33 Flat-ICE slabs was used. The implicit PCM model with 25 nodes in the x-axis (N x ), four nodes of wall (N_wall) and nine nodes of PCM (N_PCM) with a time-step of 10 s was considered as reference (the same parameters as the validated model). To demonstrate the time-step grid independence of the model, the results for the same number of nodes with a time-step of 1 s were calculated, obtaining a MAPE of 0.3% with respect to the reference. The results of both implicit and explicit schemes can be seen in Table 4. With respect to the implicit scheme, the iterative process between PCM nodes, especially during the phase change, requires high computational time. Therefore, if the number of N x nodes was kept constant at 25, but the number of nodes in the y-axis (N_PCM and N_wall) was slightly reduced from nine to seven and from four to three respectively, the computational time is reduced to almost half while maintaining a high accuracy (MAPE = 0.08%). In the same way, reducing from 25 to 15 (40%) the number of nodes in the x-axis (N x ) allowed decreasing the computational time to less than half with a MAPE of around 1.5%. If the number of nodes in both x-axis and y-axis was reduced (N X = 15, N_wall = 3, N_PCM = 7), the accuracy of the results stayed below 1.4% with much lower computational times (less than 4 times the reference time). Nevertheless, if we further reduced the N X nodes to 5, the results present low accuracy with respect to the reference (MAPE around 9.5%).
The advantage of the explicit scheme was the considerable reduction in computation times due to the lack of iterations in its scheme of resolution. For example, the explicit model with 5 N X nodes and a time-step of 1 s followed the same pattern with respect to MAPE than its implicit counterpart, but with much less computational time compared with the implicit model (6.9 s vs. 137.7 s). In the same way, the parameter combination of N X = 25, N_wall = 4, N_PCM = 9 was run with the explicit model with 1 s time-step reduced with respect to the reference (10 s time-step) from 427.7 to 43.7 s, with a slight difference of results accuracy (MAPE = 0.38%). If additionally, the number of nodes in the computationally demanding y-axis is reduced (N_wall = 3, N_PCM = 7), the computational time decreased even more (33.6 s) with still a high results accuracy (MAPE = 0.44%).
Since the goal of the analysis was to identify a model which provides the best trade-off between low computational time and results accuracy, the value (TIME × MAPE)/100% was analyzed. Based on this parameter, the best trade-off was obtained with the explicit scheme with the combination N X = 25, N_wall = 3, N_PCM = 7, for which (TIME × MAPE)/ 100% = 0.15. The ideal power range required by the evaporator of the sorption TES is between 1 and 3 kW. The reference model stayed 56.8% of the discharge time in this range, while the model with the combination N X = 25, N_wall = 3, N_PCM = 7 stays 56.7% of the time. These results supported the decision to use the explicit scheme. If a faster model is required, the explicit model with the parameter combination N = 15, N_wall = 3, N_PCM = 7 can be used, which decreases the computational time around 30% with reasonable accuracy of the model (MAPE of 1%).

Conclusions
In this paper, a two-dimensional numerical model of a storage unit using encapsulated flat PCM slabs was developed and validated against experimental data. The finite control volume approach was used to solve the energy balance equations to calculate the HTR from the PCM to the HTF and the HTF outlet temperature during the charging and discharging processes. The presented model can be used to analyze the thermal performance of PCM storage units with flat slabs, since it showed good accuracy (error of ±1% in the range of interest).
The PCM tank presented in this study was intended to provide heat to the evaporator of a seasonal sorption TES system based on selective water-sorbent materials. To obtain optimal design variables of the PCM tank that provided the desired heat powers required by the evaporator, an optimization based on the validated model of the number of slabs and their configuration into the PCM tank was performed. Two different slabs (Flat-ICE and Thin-ICE) with two different thickness (3.5 mm and 1.7 mm, respectively) were analysed. The objective function of the optimization consisted in minimizing the MAPE of the HTR during the discharging process of the tank. The optimal configuration was obtained for 33 slabs of the Flat-ICE. Due to the greater thickness of the Flat-ICE compared to the Thin-ICE, it has less heat transfer surface per kg of mass, and therefore the desired thermal power between 1 and 3 kW can be provided for a longer period.
The PCM storage unit presented in this study is part of a complex heating system. The performance analysis of the whole system and the optimization of its control will be performed through computational simulations in a future study. Computationally heavy system simulation models can cause too long computational times during the system control optimization. Hence, in order to obtain a trade-off between computational times and results accuracy of the PCM model, a comparison between the discharging process using a 2D implicit and a 2D explicit scheme for a different number of nodes was carried out. The results showed that it was not worth using an implicit scheme with a high number of nodes (N x = 25, N_wall = 4 and N_PCM = 9), since very accurate outcomes could be obtained (MAPE = 0.38%) with the same number of nodes solved with an explicit scheme, reducing the computational time from 428 s to 44 s even with a smaller time-step of 1 s. Furthermore, to study the compromise between computational time and results accuracy, the value (TIME × MAPE)/100% was analyzed. The minimum value among the cases studied was also obtained for the explicit model with the nodes combination of N x = 25, N_wall = 3 and N_PCM = 7 and a time-step of 1 s. This means that reducing the number of nodes in the y-axis to a minimum of N_wall = 3 and N_PCM = 7, while keeping a high number of nodes in the x-axis, was beneficial in terms of computational effort (33 s) and results accuracy (MAPE of 0.44%).
To sum up, the integration of the PCM tank model into the whole heating system using an explicit scheme will reduce the computational time (from 237 to 34 sec for the best case) to allow faster control optimization while keeping good results accuracy (MAPE of 0.44%).

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
The authors would like to thank the Catalan Government for the quality accreditation given to their research group GREiA (2017 SGR 1537). GREiA is a certified agent TECNIO in the category of technology developers from the Government of Catalonia.

Conflicts of Interest:
The authors declare no conflict of interest.