Compact Model of Latent Heat Thermal Storage for Its Integration in Multi-Energy Systems

Featured Application

Dynamic simulation of Latent Heat Thermal Storage at system level; optimization of control strategies in multi-energy systems; investigation of Demand Side Management strategies.

Abstract

Nowadays, flexibility through energy storage constitutes a key feature for the optimal management of energy systems. Concerning thermal energy, Latent Heat Thermal Storage (LHTS) units are characterized by a significantly higher energy density with respect to sensible storage systems. For this reason, they represent an interesting solution where limited space is available. Nevertheless, their market development is limited by engineering issues and, most importantly, by scarce knowledge about LHTS integration in existing energy systems. This study presents a new modeling approach to quickly characterize the dynamic behavior of an LHTS unit. The thermal power released or absorbed by a LHTS module is expressed only as a function of the current and the initial state of charge. The proposed model allows simulating even partial charge and discharge processes. Results are fairly accurate when compared to a 2D finite volume model, although the computational effort is considerably lower. Summarizing, the proposed model could be used to investigate optimal LHTS control strategies at the system level. In this paper, two relevant case studies are presented: (a) the reduction of the morning thermal power peak in District Heating systems; and (b) the optimal energy supply schedule in multi-energy systems.

1. Introduction

The thermal energy storage is a key element for increasing the operational flexibility of energy systems, especially when these are supplied by renewable energy sources. When the space available for storage technologies is limited, such as in urban contexts, Latent Heat Thermal Storage (LHTS) through the melting/solidification of Phase Change Materials (PCMs) can reduce the occupied volume up to five times with respect to sensible thermal storages to reach the same energy content. Therefore, an increasing amount of attention is focused on LHTS systems in order to advance their Technology Readiness Level (TRL) and market development. Currently, LHTS technologies have scarcely been introduced to market, and their TRL is lower than seven [1], which indicates that this technology has generally been demonstrated only in operational environments.

In this framework, numerous LHTS modeling approaches have been proposed in the existing literature. LHTSs are commonly analyzed through detailed numerical models, which, in few cases, are validated by experimental results. In fact, the intrinsic physical complexity of the phase change phenomenon can be easily simulated by numerical models [2]. Most of the models are focused on LHTSs design optimization with the aim of enhancing the heat transfer between the PCM and the heat transfer fluid (HTF). As reported in [3], several enhancement techniques can be pursued. Among these, numerous works investigated the optimal sizing of HTF tube fins. For instance, Niyas et al. [4] performed a numerical study of a shell-and-tube LHTS prototype filled with PCM salts with a full 3D model. Using commercial software based on a finite element scheme, they determined the optimal number of HTF tubes and longitudinal fins per tube that minimizes the overall discharging time. They also monitored a few performance parameters of the identified configuration, such as the melting fraction, charging/discharging times, outlet HTF temperature, and energy stored. Then, results were validated with experiments [5] showing good agreement. Sciacovelli et al. [6] adopted a shape optimization strategy for Y-shaped fins with one or two bifurcations in order to enhance the performance of a shell-and-tube LHTS. A two-dimensional cross-section model, solving the underlying equations with a commercial finite volume code, was implemented. The fins optimization was performed requiring the complete discharge of the LHTS unit both in a short and long-time interval. A further application of CFD modeling to LHTS fins design optimization was presented by Pizzolato et al. [7]. They extended the investigation of the fins design through the combination of topology optimization techniques and ad hoc finite element two-dimensional and three-dimensional models. Here, only heat conduction is considered. Then, the same authors increased their model complexity by including also natural convection [8].

Esapour et al. [9] developed a detailed 3D in-house code to study the performance of an annular LHTS with multiple tubes. More specifically, they determined the effect of the arrangement and number of tubes on the liquid fraction and melting time. Then, in another work [10], the same model was also adopted to investigate the influence of variable operational parameters such as HTF inlet temperature and flow rate on the LHTS performance. Similarly, Saddegh et al. [11] compared the thermal behavior of a shell-and-tube LHTS using a 2D conduction–convection heat transfer model for two different LHTS orientations (vertical and horizontal).

In an attempt to simplify the simulation of LHTS systems, Neumann et al. [12] coupled a 1D model for the heat exchanger tubes and a reduced 3D model for the PCM and fins. This model can be used to improve the geometric and operational parameters of fin-and-tubes LHTSs with a limited computational effort, although the PCM parameters should be carefully calibrated. A different simplification approach was followed by Parry et al. [13]: a 3D computational model for a shell-and-tube LHTS was calibrated with experimental results. This was successively reduced to a one-dimensional radial model by employing an effective diffusivity technique. According to the authors, this model is suited to predict performance over long-time spans. An original simple model was proposed by Tay et al. [14]. This is based on the effectiveness-number of transfer units, and it was calibrated on a tube-in-tank salt-based PCM system with radial circular fins. The correlation proposed by the authors is semi-empirical, and it can be used for sizing and optimizing LHTS systems. A further interesting simplification methodology was suggested by Johnson et al. [15]. First, they analyzed the heat transfer properties of a tube assembly cross-section through a 2D finite volume model. Then, they produced a 1D radial model in Dymola for the ensemble of PCM and fins able to yield comparable results. Finally, they coupled this latter model with a 1D axial model for the HTF again in Dymola and studied the time evolution of the thermal heat flux for different configurations.

Nevertheless, although a few studies investigated the integration of LHTS in existing heating systems, little attention has been devoted to the analysis of LHTSs operating conditions and partial charge/discharge. For instance, Colella et al. [16] analyzed the behavior of an LHTS unit in a District Heating (DH) substation through a 2D finite volume model. Xu et al. [17] studied the performance of an encapsulated LHTS unit in a residential heating system in Sweden by means of a finite element model. Johnson et al. [15] designed an LHTS for the production of high-temperature steam in a cogeneration plant using the simplified model previously mentioned. However, in order to facilitate LHTS market development, their optimal integration in energy systems is essential. Dynamic models that are fast and accurate could positively support the study of LHTSs at a system level.

In this work, a compact 0D model for a modular shell-and-tube LHTS is proposed. The model is able to predict the thermal power released or absorbed by the assembly as a function of its state of charge. Furthermore, it can represent partial charge and discharge operations. Therefore, this model can be easily deployed for simulations at a system level in order to investigate optimal control strategies or demand side management (DSM). Indeed, the proposed model is here applied to two relevant case studies regarding the integration of LHTS in DH and multi-energy systems.

2. Materials and Methods

In order to study the optimal integration of different energy technologies at the building level, a reliable representation of each energy asset is essential. In this framework, an LHTS model for the simulation of the operating conditions provides quick results while retaining an acceptable physical description. Considering the physical complexity of the phase change in PCMs and the heat transfer in shell-and-tube configurations, a two-step model simplification was adopted. First, a 2D numerical model was developed in Ansys Fluent (2020 R1). Then, a 0D mathematical description was adopted, fitting the results obtained from the former numerical model with an exponential function. The LHTS system used as a reference is a vertical shell-and-tube type, where finned tubes are crossed by water as heat transfer fluid (HTF) and its shell encloses the PCM.

2.1. D Detailed Model

Figure 1 shows the computational domain of the initial detailed model. Due to the simultaneous presence of a solid and liquid phase in the PCM, both heat conduction and natural convection are relevant heat transfer phenomena in LHTS systems. Thus, a full 3D numerical domain is generally recommended. However, the characteristics of the configuration analyzed in this study allow a substantial domain reduction. In fact, as claimed by Vogel et al. [18], compact LHTS designs hinder the development of natural convection phenomena, which are the main aspects responsible for problem asymmetries. According to [18], the Rayleigh number (Ra_W) and LHTS aspect ratio (A) can be used to determine the relevance of the wall heat flux due to natural convection in the liquid PCM. This correlation asserts that the flow is dominated by natural convection only if $\frac{R{a}_W}{A}\ge 500$. Regarding the LHTS configuration proposed in this analysis, the above-mentioned correlation yields $\frac{R{a}_W}{A}=374$. Therefore, the heat flux due to natural convection is likely to be negligible both during charging and discharging phases if compared to conduction. Hence, the computational domain is reduced to the horizontal cross-section plane. Here, the temperature gradient is expected to be larger compared to the axial direction [19], and heat conduction is dominant. The numerical domain is also further reduced, taking advantage of the symmetry created by the fins. Then, the physical problem is modeled using the finite volume method implemented in the commercial code Ansys Fluent.

Regarding the materials, the fins and tubes are made of aluminum, while the PCM is assumed to be RT70HC manufactured by Rubitherm GmbH [20]. A summary of the PCM properties is reported in

Table 1. Phase Change Materials (PCM) physical properties.

PCM Physical Property	Value
Density (ρ) [kg/m3]	880 (constant)
Specific Heat [J/(kg*K)]	2000
Thermal Conductivity (k) [W/(m*K)]	0.2
Latent Heat [kJ/kg]	214
Solidus Temperature [°C]	69
Liquidus Temperature [°C]	71

. The enthalpy–porosity technique is used for modeling the melting and solidification process [21]. Combining this method with the hypothesis of negligible natural convection, the only governing equation to be solved is the energy equation, as expressed in Equation (1).

$$\frac{\partial }{\partial t}\left(\rho H\right)- \nabla ·\left(k\nabla T\right)=0$$

(1)

where $H$ is the sum of the sensible enthalpy and latent heat of the material.

Adiabatic boundary conditions are applied on the external face (Γ₃ in Figure 1), while symmetry is considered on Γ₂. The heat transfer with the water in the pipes is modeled with an implicit Robin boundary condition on Γ₁, as shown by Equations (2)–(4). $T_{ref}$ represents the average temperature between the water inlet ($T_{in}$) and outlet ($T_{out}$) temperatures. However, $T_{out}$ depends on the efficacy of the heat transfer between the flowing HTF and the rest of the LHTS (i.e., the ensemble of fins and PCM). Therefore, if $T_{out}$ is expressed as a function of the heat flux at the wall ($\dot{{\mathrm{q}}_{\mathrm{wall}}}$), Fluent solves Equation (2) iteratively. The heat transfer coefficient $h_{conv}$ is computed using the Dittus–Boelter correlation (5) and the Nusselt number (6). In (5), n = 0.4 for discharge or n = 0.3 for charge. This correlation was preferred over the implicit one proposed by Sieder and Tate due to its simplicity and acceptable uncertainty. Substantially, the simulated domain represents an average cross-section of the whole LHTS unit.

$$\dot{{\mathrm{q}}_{\mathrm{wall}}}=h_{conv}\left(T-T_{ref}\right)$$

(2)

$$T_{ref}=\frac{(T_{in}+T_{out})}{2}$$

(3)

$$T_{out}=T_{in}+\frac{L_{tube}\pi d_i}{G{c}_{p_w}}\dot{{\mathrm{q}}_{\mathrm{wall}}}$$

(4)

$$Nu= 0.023R{e}^{\frac{4}{5}}P{r}^{\mathrm{n}}$$

(5)

$$h_{conv}= \frac{Nu\ast k_w}{d_i}$$

(6)

The numerical model follows the same approach adopted by Sciacovelli et al. [6], who validated their 2D cross-section model for PCM solidification against experimental data. For further details, the interested reader is referred to [6]. The computational grid consists of a non-structured mesh made of 8278 cells with an average cell size of 3.5 × 10 ⁻⁴ m. The selected mesh proved to be sufficiently fine not to influence the results. A Second Order Upwind scheme is used for the spatial discretization of the energy equation together with the Least Squares Cell Based method for gradient calculation. The convergence is reached when residuals are lower than 10⁻⁸ for the energy equation. On the other hand, the transient nature of the problem is approached with a Second Order Implicit Euler method, with a time-step of 1 s. The selected value proved to be sufficiently fine not to influence the results.

2.2. D Compact Model

Considering the LHTSs dynamics, the thermal power released or absorbed by LHTSs varies over time, as the thermal conductivity of the liquid and solid material is different. Indeed, as the solidification or liquefaction of the PCM proceeds, the power exchanged decreases, since the melting front is further from the fins. For this reason, it was assumed that the most relevant parameter influencing the thermal power input/output of an LHTS is represented by its state of charge (SOC). The SOC can be defined as the ratio between the actual energy stored inside the LHTS ($E_{stored}$) and its reference value when fully charged (7).

$$SOC= \frac{E_{stored}\left(t\right)- E_{min}}{E_{max}-E_{min}}$$

(7)

where $E_{max}$ is the energy contained in the PCM when the LHTS is fully charged. Similarly, $E_{min}$ is the energy contained in the PCM when fully discharged.

A mathematical correlation between the current state of charge of the LHTS (SOC) and its thermal power input/output was identified using the results of the 2D model as shown by Equations (8)–(10). These were obtained observing the shape of the 2D model results and searching for their best fit through the sum of exponential functions. Each expression refers to a different operational phase (i.e., charge, discharge, or idle phase).

$$\dot{Q_{dis}}= A{e}^{B\ast SO{C}_n}+C{e}^{D\ast SO{C}_n}+K\left(1-SO{C}_0\right)e^{-{\left(\frac{SO{C}_n-E}{F}\right)}^2}$$

(8)

$$\dot{Q_{chrg}}= A{e}^{B\ast SO{C}_n}+C{e}^{D\ast SO{C}_n}+K\ast SO{C}_0{e}^{-{\left(\frac{SO{C}_n-E}{F}\right)}^2}$$

(9)

$$\dot{Q_{idle}}=0$$

(10)

SOC_n is the normalized state of charge, and it is equal to $\left(\frac{SOC}{SO{C}_0}\right)$ during discharge and $\left(\frac{SOC-SO{C}_0}{1-SO{C}_0}\right)$ during charge, SOC₀ is the initial state of charge, $\dot{Q_{dis}}$ is the thermal power released by each tube during discharge, $\dot{Q_{chrg}}$ is the thermal power requested by each tube during the charging phase, and $\dot{Q_{idle}}$ is the thermal power when the LHTS is not operated. Therefore, depending on the operational phase of the LHTS, the value of the thermal power released or requested by each LHTS tube can be expressed as follows (11):

$$\dot{Q_{LHTS}}=\dot{x_1{Q}_{dis}}- \dot{x_2{Q}_{dis}}+ \dot{x_3{Q}_{idle}}$$

(11)

where $x_1,x_2,x_3 ϵ \left[0;1\right]$ and $x_1+x_2+x_3=1$.

Equations (8) and (9) were obtained fitting with exponential and Gaussian functions the LHTS thermal power and states of charge resulting from the simulations of the 2D model. The initial conditions of these simulations were, respectively, a fully charged LHTS and a fully discharged LHTS. Moreover, assuming that the thermal losses from the PCM to the HTF are negligible, the thermal power is always zero during the idle phase. Although these expressions were retrieved considering a specific LHTS design with fixed operating parameters (such as water mass flow rate and temperatures), the proposed methodology can be easily tailored to different LHTS concepts and operating conditions. Indeed, the evolution of the LHTS heat rate highlighted by several studies [12,15,22] is comparable with the one expressed by Equations (8) and (9).

3. Case Study

3.1. LHTS System Description

The LHTS system considered in this study is a vertical shell-and-tube type, with a height of 1.5 m. Each tube has an inner diameter of 19.05 mm (6/8′′) and is equipped with 16 longitudinal fins whose extension is 30 mm (thickness: 1 mm). Instead, the interaxial distance between each tube is assumed to be 91 mm. The PCM is characterized by an average phase change temperature of 70 °C, as reported in Table 1. The HTF inlet temperature is assumed to be 75 °C during the charging phase and 48 °C during the discharge. These values were selected analyzing the working temperatures of the heating system of the building described in Section 3.2. The HTF mass flow rate is 0.168 kg/s per tube. This value was obtained fixing an average HTF velocity of 0.6 m/s. Finally, the size of the LHTS system (in terms of stored energy) is not fixed a priori. As a matter of fact, the LHTS can be conceived as a modular system composed by the aggregation of several tubes ($n_{tubes}$) and their associated PCM depending on the considered application. Each unit (tube + PCM) is able to store 2637.2 kJ, considering both the sensible and latent heat content of the PCM (at the selected temperatures).

3.2. Model Application to Thermal Energy Networks

Among the strategies to improve the operational management of DH networks, the adoption of distributed thermal storage systems at the building level is gaining interest as a measure to shave peak demands [23]. Therefore, the first case study considered in this article is the reduction of the peak thermal request from a substation of the DH network of Torino (Italy), which is characterized by a climate zone E. This substation belongs to a building known as the Energy Center, which is a research center owned by Politecnico di Torino. Thanks to its monitoring system, the building energy data are available with a frequency of 15 min. For this case study, a representative day of January 2020 is taken as a reference. As can be seen in Figure 2, the building heat demand has a common shape characterized by a morning peak (from 6:30 to 9 am) and a rather steady thermal request during the remainder of the day. During the peak hours, the amount of heat requested to the DH network is 3100 MJ, which represents the 27.1% of the overall heat demand within the same day.

Overall, an LHTS system composed of 300 units (791 MJ) is needed in order to achieve an energy demand reduction of 25.5% during the peak. However, due to the high power ceased or absorbed by each LHTS unit at the beginning of its operation, the overall storage system is subdivided into six bundles made of 50 units (tube + PCM) as depicted in Figure 3. These bundles are sequentially activated with a delay of few minutes both during charge (in the night) and discharge (in the morning peak). Moreover, each bundle follows the subsequent duty cycle: it is charged up to SOC = 0.97 and discharged down to SOC = 0.02. Further partial charge and discharge cycles are not considered for this case study. This is done since the building heat demand is rather constant for the rest of the day.

3.3. Model Application to Multi-Energy Systems

This section shows the application of the 0D model for the achievement of the best operation in a multi-energy system context, as schematized in Figure 4. The multi-energy system is conceived to supply heat, electricity, and cold to a block of flats with 15 dwellings. The production system is characterized by a group of several technologies. In this context, the adoption of a thermal storage provides significant benefits to improve the system flexibility. In particular, an LHTS can be extremely advantageous in buildings, since it requires small installation volumes. In order to easily model the time-dependent behavior of the LHTS, the 0D model is adopted.

The most convenient technologies to be operated are selected by a Mixed Integer Non-Linear Programming optimization tool that includes the 0D LHTS model. The optimization algorithm allows achieving the best operation for the multi-energy system. Therefore, its output is the optimal time evolution of (a) the power from the production/conversion technologies; (b) the power purchase from the grid; and (c) the released/absorbed power in the storages.

The optimization model includes constraints such as (a) the overall produced and purchased energy vectors must be equal to the loads and (b) the energy production of each technology must not exceed its maximum capacity (the same for the conversion technologies and the storages). Moreover, additional constraints are set to model the LHTS. The maximum thermal power that can be absorbed/released by the LHTS is imposed at each time step in relation to the outcome of the LHTS 0D model.

Concerning the energy production/conversion, both traditional and innovative technologies are considered to be installed. Among the traditional technologies that can be adopted are a heat-only boiler (HOB) for heat production, an electric chiller (EHP) for cold production, and the grid connection for the electricity/heat supply. Among the innovative technologies are a micro-cogeneration system (for combined heat and power production) fed by natural gas, an absorption chiller, fuel cells, photovoltaics (PV), thermal solar, and wind turbine. Furthermore, a series of storages (hot, cold, and electric) provide a significant level of flexibility for the multi-energy system.

4. Results and Discussion

In order to assess the validity of the 0D model, this section compares the results with the ones obtained from the 2D model. Afterwards, the impact of an LHTS system on the operational strategies of two relevant case studies is quantified through the application of the compact 0D model. In the former case, it shows how the LHTS allows achieving a significant reduction of the thermal power peak in a District Heating (DH) substation. In the latter case, it is used to include the latent thermal storage in an optimizer that finds the optimal performance of a multi-energy system.

4.1. Comparison between 0D and 2D Models

The coefficients resulting from the fitting procedure described in Section 2 are reported in

Table 2. Fitting coefficients in Equations (8) and (9).

Coefficient	Discharge Phase	Charge Phase
A	0.1752 [kW]	3.353 [kW]
B	3.112 [-]	−45.93 [-]
C	−0.2078 [kW]	1.337 [kW]
D	−0.9345 [-]	−3.606 [-]
K	1.758 [kW]	0.3296 [kW]
E	0.5518 [-]	0.5908 [-]
F	0.3442 [-]	0.4197 [-]

. Equations (8) and (9) approximate very well the 2D power discharge and charge curves when the initial conditions are, respectively, SOC₀ = 1 and SOC₀ = 0 (Figure 5 and Figure 6). The former fitting curve is characterized by a correlation coefficient R² = 0.995, while the latter has R² = 0.992. Overall, the minimum value is R² = 0.939, while the maximum standard deviation is 0.138 kW.

The correspondence between the LHTS thermal power and its state of charge is not univocal. As a matter of fact, the LHTS thermal power depends both on the current state of charge (SOC) and on the initial state of charge (SOC₀). However, SOC₀ affects the results only in case a new charging or discharging phase starts. For instance, if there are two consecutive charging or discharging phases (only interrupted by an idle period), the value of SOC₀ should not be updated at the beginning of the second charge/discharge simulation. Instead, the previous value must be retained. This issue is due to the fact that when a partial charge is performed, only the PCM close to the fins becomes liquid, while the PCM far from the fins remains solid. If the subsequent operational phase is a discharge, a peak of thermal power occurs at the beginning of this latter process because the first layer of PCM encountered in the heat propagation is liquid (Figure 5).

A similar consideration is valid when a charge is performed, starting from a partial discharge of the system (Figure 6). The 0D model for the LHTS discharge thermal power (8) is slightly less accurate when SOC₀ is much smaller than 1 (i.e., when a discharge is preceded by a very short charging phase). However, it is rather unlikely to discharge a thermal storage starting from such a low content of energy. Similarly, the 0D model for the LHTS charge thermal power (9) is less accurate when the SOC₀ is much larger than 0 (i.e., when a charge is preceded by a very short discharge phase). As a matter of fact, the 2D curves in Figure 5 are both shifted and contracted with respect to one other, while the 2D curves in Figure 6 are mainly shifted. This dissimilarity might be due to the fact that the HTF inlet temperature during the charging phase (75 °C) is much closer to the PCM average phase change temperature (70 °C) compared to the HTF inlet temperature during discharge (48 °C). In fact, as indicated by [24], the HTF inlet temperature affects the LHTS thermal power.

4.2. Distributed LHTS in DH Networks

The 0D LHTS model is here applied to the reduction of the thermal power peak that is observed in DH networks at the beginning of their daily operation in some areas of the world (such as the Mediterranean area). In order to achieve this goal, the LHTS system follows the subsequent duty cycle, as detailed in Figure 7:

• The charging phase of each LHTS unit lasts 3 h; thus, the whole LHTS bundle is charged between 1.20 and 5.40 am (due to the imposed delay between the activation of each unit);
• The discharging phase of each LHTS unit lasts 1.5 h; thus, the whole LHTS bundle is discharged between 6.30 am and 9.20 am (due to the imposed delay between the activation of each unit).

Although all the energy absorbed by the LHTS during the charging phase is released in the discharging phase, the duration of these two processes is remarkably different. Since the HTF mass flow rate is constant, the thermal power is affected only by the temperature difference between the incoming HTF (48 °C or 75 °C) and the PCM average phase change temperature (70 °C) in both the operational phases. Consequently, the discharging phase is faster because the temperature difference is four times higher compared to the charging phase. Moreover, the spikes registered in the LHTS thermal power does not constitute a desirable feature from the point of view of the building heating system management. However, this highly unsteady behavior could be easily changed by a finer regulation of the HTF mass flow rate [16,17] (that is out of the goals of the present analysis). Furthermore, as detailed in Section 2.2, a perfect LHTS insulation is assumed. Consequently, heat losses are negligible when the storage is fully charged.

As shown in Figure 8, the time delay between the activation of each LHTS unit strongly affects the maximum power requested by the building to the DH network, whose reference value is 385 kW for this case study. Considering the LHTS system configuration analyzed, the peak thermal power of the DH supply is minimized when the activation of each LHTS unit is delayed by 16 min with respect to its preceding unit during the discharging phase. On the contrary, the minimization of the peak DH supply during the charge of the whole LHTS system would yield the obvious solution of charging each LHTS unit separately. However, this situation would not be compatible with the available time in the night. For this reason, an activation delay of 16 min was retained also for the charging phase.

Finally, Figure 9 shows the DH supply curve when the LHTS system presented in this study is deployed. For consistency with the data of the building energy demand, the DH supply curve describes the average thermal power over 15-min intervals. The morning peak reduction is evident from the figure. The thermal energy supplied by the DH between 6.30 and 9 am is reduced by 25.5%, while the maximum value of the DH thermal power is decreased by 18%. This result is achieved by simply shifting the supplied thermal energy in the night hours. This shifting could also be extensively advantageous for DH users in case of different tariffs during the daytime and nighttime. As an example, in Turin, the heat during the night is sold at a price that is up to 50% lower than the daily price [25]. Furthermore, the overall volume occupied by the LHTS system proposed in this study amounts to 3.4 m³, while a sensible water storage tank would require 7 m³ for the same energy content (assuming also the same temperature difference 75–48 °C). Therefore, this feature constitutes a great advantage when the thermal energy storage is located in the small technical rooms available in the buildings.

4.3. LHTS in Multi-Energy Systems

This section reports the results achieved in the analysis of the Multi-Energy System. The technologies activated for supplying the electricity load are the micro-cogeneration (CHP) system and the photovoltaics (PV). The extra demand is supplied purchasing the electricity from the grid; the adoption of the batteries makes the production more flexible. The electricity consumption is due to both the building load and the electricity absorbed by the electric heat pump and the electric storage. Concerning the cooling load, this is supplied through an electric heat pump (EHP). As far as heating is concerned, Figure 10 includes the power evolution of the technologies selected by the optimizer. On the left of the figure, the production evolutions are reported; on the right, consumptions are reported. The production consists in the sum of the power due to the production/conversion technologies operated, the energy purchased from the grid, and the storage discharging phase. The consumption consists in the sum of the building load, the energy sold to the grid, the storage charging phase, and the consumption of the other technologies such as the heat consumed by the absorption heat pump. The sum of the production and consumption must be equal. As a result, the heat load is mainly provided by the micro-cogeneration system, the heat-only boiler, and the thermal solar. The LHTS allows storing energy in the valley to supply it during the peak loads.

5. Conclusions

In this study, a simple 0D model for a modular shell-and-tube LHTS is presented. The model was obtained considering the main quantities affecting the charging/discharging evolution and through a fitting of the LHTS thermal power achieved with a more detailed CFD model. The compact model is obtained parameterizing the 2D model results as a function of the LHTS initial states of charge in order to account for partial charge and discharge. Results show that the 0D model is extremely accurate when compared with the 2D model, as the minimum value of the correlation coefficient (R²) is 0.939, and the maximum standard deviation is 0.138 kW.

Overall, thanks to its simplicity and low computational cost, the 0D model can be easily used for simulations at the system level. It allows investigating the performance of numerous control strategies in energy systems equipped with an LHTS, even in case the LHTS is only partially charged or discharged. In this work, two specific relevant applications of the compact model are considered. The former regards the thermal power peak reduction in DH substations, while the latter concerns the optimal operational strategy of a set of production/conversion/storage technologies (included LHTS) in a multi-energy system. However, this study describes the functioning of the LHTS at the system level only on the basis of numerical results. Therefore, future work should test the LHTS also in a real case application to better quantify the model uncertainties.