**1. Introduction**

The continuous increase in the importance of the role of energy over the last few decades, as well as the rise in fuel prices and the need to limit greenhouse gas emissions, have led to a steady growth in the use of energy saving technologies and in a more effective and extensive implementation of renewable energy sources [1]. However, despite renewable energies representing one of the best alternatives to conventional sources—such as fossil or nuclear—for energy supply in most areas of the world, renewable energies are often hampered by their discontinuous nature during the day and by the actual availability of the source during the year.

As a matter of fact, energy is not always available where, when, and how it is required, and thus storage systems (both electrical and thermal) play an important role in energy managemen<sup>t</sup> to make it available accordingly. Therefore, to improve the availability of renewable energies in remote geographical areas and to overcome their intermittent nature, thermal energy storage (TES) represents a fundamental solution to increase their competitiveness. Thanks to its capability to allow a more sustainable use of available resources [2] and to decouple thermal energy generation and use, there is currently growing interest in thermal energy storage.

A number of applications for thermal energy storage can be seen in energy grids, since they allow (i) an e ffective balance in energy supply vs. demand dynamics, (ii) a decrease in heating system energy losses, reducing the number of start-up and shut-down maneuvers and the need for backup plants, (iii) an increase in deliverable capacity (heating element generation plus storage capacity), (iv) a shift in energy purchases to lower cost periods, and (v) an increase in renewable energy source exploitation.

As an example, Barbieri et al. [3] showed the storage capability to decouple electrical and thermal power production in cogeneration plants, allowing thermal energy to be stored when only the electric energy is needed and vice-versa. It is therefore apparent that thermal energy storage represents a key solution in all applications in which energy supply and demand are decoupled, giving significant advantages in terms of cost and dispatchability of the generated energy. It can also be used to support generation by conventional energy sources [4] to deal with weekly, monthly, and annual changes in energy requirements and to help peak shaving [5]. Moreover, thermal energy storage represents a fundamental technique for the optimization of overall energy conversion, transmission, and utilization processes in smart energy systems [6]. Ultimately, it has a fundamental role in energy system control strategy development [7], such as Model-in-the-Loop applications, and implementation [8,9].

Several solutions have been proposed in the literature for thermal energy storage (e.g., through sensible heat, latent heat, and thermochemical energy) [10,11] involving the most recent areas of investigation and experimentation. A detailed analysis and a performance comparison of the existing thermal energy storage systems have been performed by Sarbu and Sebarchievici [12]. Furthermore, Noro et al. [13] focused their attention on the liquid sensible and phase change material heat storage systems, pointing out their respective pros and cons. However, because of the limitations and drawbacks of latent heat and thermochemical storage solutions—e.g., large volume variations and high storage medium costs, respectively—water is still the most widely di ffused energy storage material since it is easily available, cheap, non-toxic, non-flammable, and completely harmless [14]. For these reasons, this study focuses on water sensible thermal energy storage.

A multi-node model has been developed, and the e ffects of node number and flow rate variations have been investigated alternatively. The model has been developed with a modular approach, by considering a standardized input/output causality structure for easy implementation into a library of components for the dynamic simulation of smart energy systems [6]. Unlike other simulation tools that are currently available [15,16], the present paper shows all the fundamental equations and their implementation in an explicit way allowing easy replicability and providing high-level customization features for both discretization and numerosity of inlet or outlet ports. Then the storage model behavior has been analyzed in charge and discharge conditions, taking the experimental data from the literature as a benchmark (i.e., the experimental model developed by González-Altozano et al. [17] and the experimental study carried out by Li et al. [18], respectively). Finally, the model has been applied to a real case in an operating environment.

## *Literature Review*

Accurate surveys of mathematical simulation models of TES tanks were conducted by Njoku et al. [19], and Dumont et al. [20]. It follows that these models can be divided into three main categories:


Starting from the 0D models, following the analytical approach, the storage tank can be modeled as a semi-infinite body, assuming that the inlet temperature and the mass flow rates are constant and not considering mixing or ambient heat losses. It should be noted that some attempts to relax these assumptions have been successfully performed in [21,22]. The fully-mixed hypothesis is still the simplest approach, since uniform temperature distribution is considered in the whole tank volume, assuming that all the incoming water mixes perfectly with the water already in the tank. Governing equations are based on mass and energy conservation [23].

Among the 1D models, the two-zone moving boundary approach is based on the assumption that the storage volume is divided into two zones with uniform temperature separated by a thermocline to prevent the hot and cold fluid from mixing and to maintain a stable temperature gradient [24]. In order to reduce the mathematical complexity of the analytical solutions, for the sake of practicability, an integral approximate approach has been successfully applied to the 1D storage tank models by Chung and Shin [25]. In the plug-flow method the tank volume is composed of a variable number of isothermal disks moving within the tank without any mixing between them [26]. Finally, the multi-node models are developed following a 1D finite-volume method that considers a uniform horizontal temperature in each layer (called a node). Therefore, a temperature gradient in the vertical direction only (i.e., stratification) and a 1D flow inside the tank [27,28] can be modeled. Using this 1D approach, González-Altozano et al. [17] have proposed a new methodology for the estimation of numerous temperature-dependent indices employed for the characterization of thermal stratification in water storage tanks (as temperature profile and thermocline thickness) during the charging phase, while Li et al. [18] have investigated the influence of several inlet structures on the stratification effectiveness and discharging performance of a tank.

Within the 2D-3D models, the zonal approach is a 3D finite-volume method with a coarse mesh based on mass and energy balances [29], instead the CFD [30,31] is a finite element method governed by mass, energy, and momentum conservation laws.

Each approach shows advantages and disadvantages depending on its specific application: the pros and cons of each model category are briefly summarized in Table 1. It is shown that the 2D-3D models allow for a more detailed simulation of physical systems compared to the 1D models, but with a higher computational burden because of the increase in dynamic states (resulting in a larger number of differential equations to be solved). Concerning the comparison between the reliability of different models, an insightful numerical study has been proposed by Cabelli [32], which states that, under particular circumstances, the temperature profile of the stratified fluid can be reasonably predicted with a simpler 1D model instead of a 2D model. Thus, in this paper specific attention has been paid to 1D models as they represent the best compromise between accuracy (which is related to the number of differential equations representing the model) and computational time (Table 1).

It is widely known that several physical phenomena occur in a water storage system that affect its performance. Among them, stratification has particularly raised the interest of researchers [27,28,33] for its strong influence on the thermal performance of plants [34,35].


**Table 1.** Pros and cons for each thermal energy storage (TES) tank modeling approach.

For instance, Van Koppen et al. [36] and Furbo and Mikkelsen [37] found that thermal stratification in solar heating systems allows a reduction in temperature at the collector inlet (which increases its efficiency) and leads to a limited operation of the auxiliary energy supply. Accordingly, stratification improves not only the water tank efficiency, but also that of the whole extended system.

Taking into account the scope of this work, the 1D approach has been chosen by focusing the attention on multi-node architecture. Even though several approaches for multi-node modeling of thermal energy storage tanks have been proposed in the literature (as reported above), in general they do not allow variations in the number and size of nodes considered in the simulation. To allow for a more flexible customization, a new multi-node model has been developed by the authors and is presented in this paper. Both the dimensions and number of nodes of the model can be set, focusing on the desired zone of the tank, by increasing the number of nodes in that specific area.

#### **2. Materials and Methods**

The model of the thermal energy storage developed by the authors falls under the multi-node 1D model category. The basic approach consists of dividing the tank into a number N of zones (called "nodes") where temperature is considered uniform (Ti). This assumption considers the vertical gradient of the temperature in the tank, while the horizontal gradient is neglected. When numbering the nodes from the top of the tank (Figure 1), the first node is the hottest and node "N" is the last and coldest.

As the developed model is physics-based, the involved equations arise from volume (1) and energy balances (2). Thus, the continuity equation and energy conservation equation are implemented for each node in the following form:

$$\frac{d\mathbf{V}}{dt} = \sum \dot{\mathbf{V}}(\mathbf{t}) = 0\tag{1}$$

$$\frac{\text{dE}}{\text{dt}} = \sum \mathbf{h}(\mathbf{t}) \cdot \dot{\mathbf{m}}(\mathbf{t}) + \sum \dot{\mathbf{Q}} \tag{2}$$

The energy Equation (2) involves enthalpy flows and heat exchanges not related to mass flows. The latter are due to heat transfer inside the tank between neighboring nodes, and to ambient losses through the tank walls because of temperature differences between the storage medium, tank walls, and the surrounding environment. The continuity and energy balance equations are implemented for each node, considering enthalpy flows, heat transfers between adjacent nodes and ambient losses. A schematic representation of the generic node is shown in Figure 2. The governing dynamic equations for the ith node are ordinary differential equations representing volume (3) and energy balances (4):

$$
\dot{\mathbf{V}}\_{\text{vert}\_{\text{out}}} = \frac{\dot{\mathbf{m}}\_{\text{in}}}{\rho(\mathbf{T}\_{\text{in}})} + \dot{\mathbf{V}}\_{\text{vert}\_{\text{in}}} - \frac{\dot{\mathbf{m}}\_{\text{out}}}{\rho(\mathbf{T}\_{\text{i}})} \tag{3}
$$

$$\begin{array}{l} \mathbf{M} \cdot \mathbf{c} \cdot \frac{\mathbf{d} \mathbf{T}\_{i}}{\mathbf{d} \mathbf{t}} = \sum \dot{\mathbf{m}}\_{\text{in}} \cdot \mathbf{c} \cdot \mathbf{T}\_{\text{in}} - \sum \dot{\mathbf{m}}\_{\text{out}} \cdot \mathbf{c} \cdot \mathbf{T}\_{\text{i}} \\\ -\dot{\mathbf{V}}\_{\text{vert}\_{\text{out}}} \cdot \mathbf{p} (\mathbf{T}\_{\text{vert}\_{\text{out}}}) \cdot \mathbf{c} \cdot \mathbf{T}\_{\text{vert}\_{\text{out}}} + \dot{\mathbf{V}}\_{\text{vert}\_{\text{in}}} \cdot \mathbf{p} (\mathbf{T}\_{\text{vert}\_{\text{in}}}) \cdot \mathbf{c} \cdot \mathbf{T}\_{\text{vert}\_{\text{in}}} \\\ +\mathbf{k} \cdot (\mathbf{T}\_{\text{i}-1} - \mathbf{T}\_{\text{i}}) - \mathbf{k} \cdot (\mathbf{T}\_{\text{i}} - \mathbf{T}\_{\text{i}+1}) - \mathbf{U} \cdot \mathbf{A} \cdot (\mathbf{T}\_{\text{i}} - \mathbf{T}\_{\text{amb}}) \end{array} \tag{4}$$

The water density varies with the temperature, and it is thus assumed constant in each node. It must also be specified that . Vvertin and . Vvertout are positive when water is flowing downwards. Accordingly, the sign of Tvertin and Tvertout is defined by the following relationships:

$$\mathbf{T\_{vert\_{out}}} = \begin{cases} \begin{array}{ll} \mathbf{T\_{i}} & \text{if } \dot{\mathbf{V}}\_{\text{vert\_{out}}} > 0 \\ \mathbf{T\_{i+1}} & \text{if } \dot{\mathbf{V}}\_{\text{vert\_{out}}} < 0 \end{array} \end{cases} \tag{5}$$

$$\mathbf{T\_{vert\_{in}}} = \begin{cases} \mathbf{T\_{i-1}} & \text{if } \dot{\mathbf{V}}\_{\text{vert\_{in}}} > 0 \\ \mathbf{T\_{i}} & \text{if } \dot{\mathbf{V}}\_{\text{vert\_{in}}} < 0 \end{cases} \tag{6}$$

The energy conservation equation was applied to each node in order to predict its time-temperature history taking into account the thermal losses both to the surroundings and to adjacent nodes because of the convection and conduction. The latter is represented in Equation (4) by the terms k·(Ti−<sup>1</sup> − Ti) and k·(Ti − Ti+<sup>1</sup>), which can be appointed as "pseudo-conduction terms" since they are used to represent the thermal exchange that would take place between bordering nodes by convection.

**Figure 1.** Liquid storage tank schematic representation for the multi-node modeling approach.

**Figure 2.** Generic node schematic representation.

This is a simplified way to consider these contributions (as suggested in [17,28]), but it has been noted by a preliminary analysis, which is not reported herein for the sake of brevity, that their influence on the global energy balance could be neglected. For this reason, this is an approach which is commonly used when dealing with multi-node models [15,16]. The terms . Vvertout ·<sup>ρ</sup>(Tvertout)·c·Tvertout and . Vvertin ·<sup>ρ</sup>(Tvertin )·c·Tvertin represent the energy associated with the incoming and outgoing flows among neighboring nodes and can be considered as "transport terms."
