LargeTESModelingToolkit: A Modelica Library for Large-scale Thermal Energy Storage Modeling and Simulation

This paper introduces the LargeTESModelingToolkit, a novel Modelica library for modeling and simulation of large-scale pit and tank thermal energy storage. This first comprehensive Modelica library in the field provides the flexibility and tools needed to develop new storage models tailored to the desired application. It also offers researchers and industrial users pre-built storage models for simulation studies to answer the relevant questions for an optimized design at storage and system level. In this paper, we present the library’s key features and structure and introduce the underlying physical and mathematical foundations and modeling approaches. More-over, we discuss the validation of the models, present the first results, and show the library’s applicability using an exemplary simulation case study.


Introduction
The integration of large-scale underground hot-water tank and pit thermal energy storage systems offers a high potential to considerably increase the share of renewable energy in future local and district energy systems.These large-scale thermal energy storage (TES) technologies can provide the flexibility needed to store volatile renewable energy sources for a few days as well as on a seasonal basis, bridging the natural gap between supply and demand (Schmidt et al. 2018).At the same time, they also offer a high economic attractiveness for storing large amounts of heat due to economies of scale and a certain flexibility in site selection due to attractive underground integration without free-standing tall structures.In contrast, the large volumes involved lead to high investment costs, which require fundamental planning at the component, storage, and system level.Experimental investigations in the design phase are limited due to the size of these storage technologies and the long time periods in question.Therefore, numerical simulations from the component to the system level are used throughout the whole design process, from the feasibility phase to the detailed design (Dahash, Ochs, Janetti, et al. 2019).
To date, for large-scale tank (TTES) and pit (PTES) thermal energy storage systems TRNSYS1 is the most widely used simulation tool for storage and system design questions as well as for scientific studies.For example, an overview of past studies and the used TRNSYS models can be found in Xiang et al. (2022).
Modelica TES models are up to now mainly focused on free-standing models (i.e., without modeling of the surrounding ground) (Leoni et al. 2020).These models are primarily used to simulate hot water storage tanks for domestic applications in small size ranges (i.e., up to a few cubic meters) and are, for instance, included in the Modelica Buildings library (Wetter et al. 2014).In recent years, dedicated TTES and PTES models have been developed incorporating modeling of the surrounding ground.Dahash, Ochs, and Tosatto (2020) demonstrated a model with cylindrical geometry and conducted a cross-comparison with generic boundary conditions between the Modelica model and a model developed in COMSOL Multiphysics 2 .Moreover, Reisenbichler et al. (2021) also developed a TTES model with cylindrical geometry in Modelica.To assess the accuracy of the developed model, a validation case study has been conducted by comparing the simulation results with real measurement data of a Danish PTES in Dronninglund, alongside a crosscomparison against other numerical models.A similar Modelica model with cylindrical geometry was developed by Fournier (2022) and also validated against the Dronninglund PTES measurement data.Recently, Kirschstein (2022) developed a Modelica PTES model with a conical geometry, which is available as part of the Modelica Solar District Heating (MoSDH) library (Formhals 2022).
Modelica's equation-based, object-oriented, and multidomain modeling approach inherently facilitates high reusability, expandability and adaptability of the produced models.Therefore, with Modelica's features and capabilities, large-scale TES modeling on the storage and system level can be taken to the next level.
However, currently available models are still scattered, limited in functionality and flexibility (e.g., in the choice of geometries), focused on specific applications, or have not yet undergone a comprehensive validation process.Consequently, our overall goal is to enable more efficient modeling and simulation of large-scale TES in multiannual dynamic system and storage simulations.Therefore, with the development of the Modelica LargeTES-ModelingToolkit (LargeTESmtk) library, we aim to provide a comprehensive toolkit for the modeling and simulation of large-scale pit and tank TES.In addition to an easy-to-use library with scientifically proven, pre-built storage models for researchers and industrial users, the library is also intended to provide the foundation and tools for the development of new storage models customized to the wanted application.
This paper focuses on the presentation of this Modelica library.We start with an overview of the main features, the implementation in Modelica, and the modeling approaches of the main models.Then we give a brief insight into the models' ongoing validation process and show the library's application in an exemplary simulation case study.

The LargeTESModelingToolkit library
This section introduces the LargeTESmtk by giving an overview of its main features and models.Afterwards, we describe the corresponding Modelica library structure.

Overview
Figure 1 provides an overview of the library.Displayed in the center is the basic structure of each storage model, consisting of the two main models for the fluid and ground domain, surrounded by an extract of available model configuration options.Along with other features, the library should provide a wide range of model configuration options in terms of geometry (e.g., cylindrical, conical, or hybrid geometries), heat transfer mechanisms (e.g., pure convection or combined convection and radiation) or ground properties (e.g., uniform ground or specific ground layers) that can be tailored specifically to the wanted application and level of detail.
Underground TTES are typically surrounded by vertical retaining walls to withstand the horizontal loads (e.g., lateral earth pressure of the enclosing ground or hydrostatic water pressure) and have cylindrical or cuboidal geometries.PTES, on the other hand, typically have conical or pyramidal geometries with slopped walls and hence do not require retaining walls to accommodate the horizontal loads (Pauschinger et al. 2020).Both PTES and TTES can be either fully buried (below original ground level) or partly buried with the excavated soil as embankments.In addition, the integration of large-scale TTES in district heating grids as free-standing steel tanks for mainly shortterm heat storage is widely used (Dahash, Ochs, Janetti, et al. 2019).As shown in Figure 1, the library provides all the necessary building blocks to model the mentioned storage types and geometry options.
Eventually, the LargeTESmtk library should offer the following key features and benefits: • High adaptability, extensibility, and reusability of the models and sub-models • Large portfolio of configuration options for initial model generation and later customizations (e.g., concerning geometry, ground properties, heat transfer mechanisms) • Adaptable level of detail, enabling an applicationoriented adjustment between accuracy and calculation performance • Broad range of application of the models from preliminary design to detailed system and storage design studies • Simple integration and coupling with other relevant system components (e.g., solar thermal systems) for holistic investigations at system level The developed models with the LargeTESmtk are to be applied in parameter studies, sensitivity, and technoeconomic analyses for optimized design on storage and system level.This may include addressing important storage design questions regarding the volume, geometry, insulation quality of the cover, side walls, and bottom, or the number and position of inlets and outlets (i.e., diffusers).In addition, for instance, the investigation of long-term effects (e.g., the development of storage performance in the first years of operation during the heat-up of the surrounding ground) of different system integration concepts (e.g., post-heating concepts via large heat pumps) or storage operation strategies is possible.Finally, the developed models can be used in conjunction with appropriate case studies and methods to obtain relevant techno-economic key performance indicators for decision-making and project planning, such as storage efficiency, thermal losses, stratification quality, investment costs, levelized cost of storage and heat, primary energy consumption, CO 2 emissions or savings.
It is also necessary to point out current limitations of the library models.Mainly to reduce the numerical effort, the ground domain model is restricted to axisymmetric geometries.Thus, non-axisymmetric strorage geometries (e.g., pyramids or cuboids) cannot be modeled directly, but are implemented by corresponding parametrizations and geometry transformations.Furthermore, this approach leads to limitations in the modeling of occurring groundwater flows.

Modelica library structure
Figure 2a shows the top-level library structure of the LargeTESmtk.The Modelica library reflects the models and features described in the overview in subsection 2.1.The structure is oriented around the Modelica Standard Library (MSL) (Modelica Association 2020) and therefore contains common packages such as Utilities, Types, Icons or BaseClasses.The Examples package is intended for models demonstrating the application of the library.The TankTES and PitTES packages contain basic prebuilt storage models for different fluid domain geometries (e.g., cylindrical or conical geometries) and construction types (free-standing, partly buried, fully buried) as shown in Figure 2b.Also included are Validation sub-packages providing validation examples for the individual storage models.The HybridTES package is intended to include storage models for hybrid geometries (e.g., a combination of a cylinder at the top and a truncated cone at the bottom).The StorageComponents package contains the models from which the pre-built models are assembled and contains various fluid and ground domain models for different geometries and component models for covers or side walls.As such, this package also contains the core components for building new storage models.

Methods
This section explains the main modeling approaches of the fluid and ground domain models applied in the LargeTESmtk.For both models, the finite difference method (FDM) is used to solve the governing basic partial differential equations (PDE).The so-called "methodof-lines" is applied to replace the spatial derivative terms in the PDEs with algebraic approximations (finite differences), leading to a system of ordinary differential equations (ODE) and differential algebraic equations (DAE), including only time-dependent functions, which can be solved with the common ODE-solvers in Modelica-based simulation environments (Fritzson 2015).Furthermore, instead of deriving the corresponding finite difference formulations directly from the basic PDEs, we will use the energy balance approach here since it is considered more intuitive (Çengel and Ghajar 2015).The energy balance approach is based on subdividing the respective calculation domain into a sufficient number of volume elements (control volumes) and subsequently forming energy and mass balances for each element.Accordingly, in the middle of each volume element are the nodal points (nodes) at which the temperatures are to be determined.

Fluid domain model
Applying the energy balance method, the fluid region is subdivided along its axial direction into equidistant volume elements with a uniform temperature per element.Then, the energy balances are established and the resulting ODEs for each element are solved.This approach is widely used in dynamic system simulation tools and in the literature often referred to as 1D multi-node model or approach (Untrau et al. 2023).Since we assume that the storage fluid is incompressible and the storage is always fully filled, the formation of the mass balance can be omitted (Powell and Edgar 2013).Figure 3 shows a schematic representation of the fluid domain modeling approach with the respective geometrical parameters and energy flows.The energy balance for each volume element [n] can be expressed as the sum of all incoming and outgoing enthalpy flows Ḣ and heat flow rates Q equal to the change in internal energy U or, in case of constant thermophysical properties of the storage fluid (water), the temperature change T of the element with time t: where ρ and c p are the density and the specific heat capacity of the storage fluid.
The following equations describe the calculation of the necessary geometrical parameters for a conical fluid domain geometry: (5) The top and bottom surface areas of the volume element The internal enthalpy flows Ḣ[n] and Ḣ[n−1] result from the interaction between the adjacent upper and lower volume elements and depend on whether a downward flow (typically during charging) or an upward flow (typically during discharging) in the storage occurs.Additional enthalpy flows Ḣin [n] or Ḣout[n] may occur if the respective fluid element is also used for the external volume flows for charging and discharging the storage: The sum of all heat flow rates of each element results from the heat conduction between the adjacent elements, the heat losses to the surroundings, and a buoyancyinduced heat flow rate that may arise: The prevailing heat flow rates due to heat conduction between the adjacent elements Q[n] and Q[n−1] are calculated with the thermal conductivity of the storage fluid k, the heat transfer area A, the distance between the fluid nodes ∆z [n] and the corresponding temperature difference ∆T .The lateral heat losses Qs[n] are derived with the overall heat transfer coefficient U s[n] , composed of the inner convective heat transfer coefficient and the thermal resistance of the wall, the heat transfer area A s[n] and the temperature difference between the fluid element T [n] and the surrounding ground T s,g [n] : Similar to Qs[n] , additional heat losses to the top Qt and the bottom Qb occur for the first and the last fluid element.
A buoyancy model is applied to account for buoyancyinduced natural convection in the storage when temperature inversion occurs (i.e., a higher fluid layer has a lower temperature than the layer below).Instead of the buoyancy-induced volume flow into the adjacent fluid layer above that occurs in reality, this volume flow is emulated by adding a corresponding heat flow rate Qbuo[n] to the fluid element and can, for instance, be expressed as (Wetter et al. 2014): where k buo is a proportionality constant since a heat flow rate is used instead of a volume flow, and τ is a time constant that determines how fast the temperature compensation between the fluid layers occurs.The Fluid.Storage.Stratified model of the Modelica IBPSA library (IBPSA 2018) served as the basis for the Modelica implementation of the fluid domain model in the LargeTESmtk.However, as described above, multiple extensions and adjustments to the model (e.g., extension to conical geometry, coupling with ground model) were made.

Ground domain model
The basic mathematical description and the governing equations for the ground domain model follow the partial differential equations of two-dimensional transient heat conduction in cylindrical coordinates with constant thermophysical properties.Again, we use the energy balance approach to derive the corresponding ordinary differential equations for each element, which are then solved.Figure 4 shows a schematic representation of the ground domain modeling approach with the respective geometrical parameters and energy flows.The energy balance for each volume element [m, n] can be expressed as the heat flow rates Q into the element from the top, bottom, left and right surface equal to the change in internal energy U or, in case of constant thermophysical properties, the temperature change T of the element with time t: where ρ and c p are the density and the specific heat capacity of the ground.The volume of the individual ground elements results from: with the respective geometrical parameters illustrated in Figure 4.The heat flow rates follow the basic equation Q = G • ∆T , where G is the thermal conductance of the ground, the product of the thermal conductivity k and the corresponding geometrical relationship, and ∆T the temperature difference between the adjacent volume element and

Validation
We are going to present in this paper only a part and excerpts of the validation and cross-comparison studies that have already been carried out.For details, please refer to the respective publications mentioned below.Furthermore, not all of the studies conducted have been published yet, but we will give a brief insight into the ongoing work.
Mainly at the beginning of the development of the models and the library itself, certain sub-models were compared with analytical solutions (e.g., steady-state onedimensional heat conduction) and with available similar models of other Modelica libraries.These examples are not presented here, but are included in the respective Validation sub-packages in the Modelica library.
To validate and assess the accuracy of the developed TTES model with a cylindrical fluid geometry, a validation case study was conducted comparing the simulation results with real measurement data of the Danish PTES in Dronninglund (Reisenbichler et al. 2021).Figure 5 shows an excerpt of this study with the comparison between the simulated and measured storage temperatures.In summary, the validation case study revealed that the storage temperatures as well as the charged and discharged energies (with deviations in the range of 1%) could be accurately represented compared to the measurement data.Somewhat higher deviations from the measured data (in the range of 10%) were only seen in the simulated total thermal losses, particularly in the side and bottom thermal losses.Presumably, this is due to the differences between the cylindrical model geometry and the actual pyramidal geometry of the real storage, despite adjusting certain model parameters to the geometry of the real storage (e.g., assuming constant thermal conductance values of the top, side and bottom surfaces between both geometries and corresponding adjustment of the overall heat transfer coefficients).This deviation in thermal losses must be considered in detailed design studies, in particular when the model cannot accurately represent the actual storage geometry.
Furthermore, the cylindrical TTES model of the LargeTESmtk was part of a cross-comparison study of various large-scale TES models from different simulation tools (COMSOL Multiphysics, TRNSYS, Modelica/Dy-5.Comparison between simulated and measured storage temperatures of the PTES in Dronninglund for the year 2015 in daily resolution; a) at the three inlet and outlet diffuser heights; and b) at the top, between top and middle, and middle and bottom diffuser (Reisenbichler et al. 2021).In the corresponding colored boxes, the coefficient of determination values (R 2 ) for the entire year are shown.mola3 and MATLAB/Simulink4 ) (Ochs et al. 2022).In this study, all models were examined in a scenario with generic boundary conditions, in which the system was neglected and emulated by a simplified charging and discharging profile, and with different insulation levels (e.g., insulated and non-insulated).The results were compared in terms of storage and ground temperatures, charged and discharged energy, and thermal losses.Overall, a good agreement between the LargeTESmtk TTES model and the other models with respect to the compared temperatures as well as energy values could be demonstrated.
The same study was extended to a cross-comparison with PTES models.
In addition, a similar crosscomparison study is conducted in the course of the IEA ES TCP Task 39 "Large Thermal Energy Storages for District Heating"5 .In both studies, the LargeTESmtk PTES model with conical geometry is included and the results show good agreement with the other involved models.However, the results of both studies have not yet been published.
Further validation studies of the models from the LargeTESmtk (including further comparisons with measurement data from real TES systems) are in progress.A main application area of the storage models is the integration in (multi-annual) system simulations.Thereby, the interaction of the storage together with other system components (e.g., large heat pumps, district heating systems) is evaluated on system level.For example, the storage models of the LargeTESmtk have already been used in case studies of large-scale TES with volumes up to millions of cubic meters integrated into DH systems for storing heat from geothermal and solar thermal plants (O'Donovan 2020).
However, in this paper, we demonstrate the application of the library models through an exemplary simulation case study on storage level.The study compares the performance of tank and pit TES with cylindrical and conical fluid geometry under different operation modes ranging from short-term to seasonal storage operation.For this purpose, we will use the prebuilt storage models TTESCylinderFullyBuried and PTESConeFullyBuried of the LargeTESmtk.

Case study description
As the focus is on the storage level, we neglect other system components such as heat supply units and heat consumers (e.g., the district heating system).Instead of the system components, simplified generic charging and discharging profiles are applied for the different operation modes.Yet, the operation modes are based on real storage application scenarios.The seasonal operation is based on the application of the storage in a solar district heating (SDH) system and the short-term operation on the application for the optimization of a combined heat and power (CHP) plant (Pauschinger et al. 2020).
The number of nominal storage cycles per year resulting from the different operation modes is specified in Table 1.Each full storage cycle (t cycle ) consists of a charging phase t ch , a storage phase (TES in the charged state) t store , a discharging phase t dis and an idle phase (TES in the discharged state) t idle .During the charging phase, the inlet volume flow rate Vch,in and temperature T ch,in at the top diffuser and, during the discharge phase, Vdis,in and T dis,in at the bottom diffuser (implemented as step functions) are applied.T ch,in and T dis,in are assumed to be 95 °C and 55°C, across all operation modes.
Table 2 shows the applied model parameters in terms of storage dimensions, fluid, ground, and cover properties, and the applied heat transfer coefficients (HTC).The thermophysical properties of the storage fluid (water) are assumed to be constant at a temperature of 75 °C based on the prevailing storage temperatures (Kretzschmar and Wagner 2019).The ground properties were derived from general values for unconsolidated rocks of gravel, sand, and clay/silt, in both dry and water-saturated conditions (Verein Deutscher Ingenieure 2000).Both investigated storage types are assumed to be insulated only at the top (with a slight extension of the insulation layer beyond the storage edges), while the side walls and bottom remain uninsulated.The applied properties of the storage cover are based on the currently deployed cover constructions (after their revision) of the Danish PTES in Marstal and Dronninglund (Bobach 2020).Only convective heat transfer is considered for the cover and ground surface with the ambient air.The corresponding ambient temperatures of a typical meteorological year (TMY) for the time period between 2015 and 2020 for a generic Central European location (in this case, Vienna) were obtained from the PVGIS online tool (Huld, Müller, and Gambardella 2012).Radiative heat transfer mechanisms are not considered.
To ensure a quasi-stationary storage operation and to neglect the effect of the heat-up phase, the simulation time is five years and the results are only evaluated for the last simulation year.Accordingly, the extent of the ground domain is chosen with a distance of 50 m in axial and radial direction from the fluid domain.Thus, the ground domain is sufficiently large so that the outer nodes remain unaffected by the fluid domain and adiabatic boundary conditions can be applied for the lateral and bottom ground domain boundaries.Moreover, we shifted the simulation start time to the beginning of May to start directly with a charging phase for the seasonal operation.Consequently, the ambient temperature profile is also considered with a corresponding time shift.

Results and discussion
Various definitions of energy and exergy efficiencies (e.g., with or without consideration of the difference between internal energy content or specifically for the application on seasonal storage) are used and discussed in the literature (Dahash, Ochs, Janetti, et al. 2019;Sifnaios et al. 2022).However, mainly to give a first indication, this study uses the simple definition of the annual storage efficiency η as the annual discharged energy Q dis divided by the annual charged energy Q ch : As identical charging and discharging profiles are continuously repeated for each nominal storage cycle and the simulation time is five years, the influence of the change in internal energy content is very small.Furthermore, this definition of storage efficiency allows for a comparison across all operational modes and can be interpreted simply as the ratio of energy recovered to energy stored, regardless of the observation period between the longest and shortest nominal storage cycle duration.
Figure 6 shows the resulting storage efficiencies depending on the operation mode (i.e., number of nominal storage cycles) ranging from short-term to seasonal storage operation.The expected trend that a higher number of storage cycles leads to higher storage efficiency is clearly evident for both storage types.Thus, the short-term operation with a number of 120 nominal storage cycles shows the highest efficiencies with values above 99%.With a lower number of storage cycles, the storage utilization decreases and the actual storage phases become longer, resulting in the lowest storage efficiency for the seasonal operation, whereby the storage efficiency only falls below 90% starting from a number of four storage cycles.The lowest efficiency is about 58% for the PTES case, which is in the range of the measured storage efficiencies for the PTES in Marstal (Sifnaios et al. 2022).
As expected, the TTES generally performs better than the PTES.One main reason for that is the lower surfaceto-volume ratio of the cylindrical geometry compared to the conical geometry.At high storage cycle numbers, there is only a slight difference between the storage types because of the short storage phases and the high storage utilization.However, for seasonal storage operations, the difference between TTES and PTES is considerable at about 17%.
As the computational performance is an important factor for the applicability of the models, especially for parameter studies with numerous simulation runs, a first in-LargeTESModelingToolkit: A Modelica Library for Large-scale Thermal Energy Storage Modeling and Simulation sight of the needed calculation time will be given here.Since the calculation time depends on many factors (e.g., model discretization, number of other system components, solver settings), detailed analyses are the subject of further studies.The simulations required approx.30 minutes for each single TTES case and 45 minutes for each PTES case6 .Considering the rather long simulation time of five years per case, these calculation times are considered to be within an acceptable range for more extensive parameter studies.
It is important to mention that the main purpose of this rather small simulation case study was to show the application of the storage models of the library.Storage parameters such as volume, height-to-diameter ratio, ground properties or insulation level, which can have a high impact on the storage efficiency, were not evaluated.However, this simulation case study can serve as a basis for these broader parameter studies.Furthermore, with the provided models and flexibility of the LargeTESmtk, the study may also be extended to other storage geometries or combined with additional system components for evaluations on system level.

Conclusion and outlook
This study introduced the LargeTESmtk, a Modelica library for the modeling and simulation of large-scale pit and tank TES at storage and system level.This comprehensive Modelica library provides pre-built storage models for multi-annual dynamic system simulations for researchers and industrial users.In addition, the provided sub-models and the wide range of model configuration options (e.g., in terms of geometry, heat transfer mechanisms, or ground properties) allow high flexibility in the modeling process and facilitate the development of new models specifically tailored to the desired application.
The library models are to be applied in simulation studies to address relevant storage design questions (e.g., regarding the proper storage geometry or insulation quality) or to investigate different system concepts to achieve an optimized design at storage and system level.These simulation studies also form the basis for techno-economic evaluations to obtain the relevant key performance indicators, such as storage efficiency and levelized cost of heat, for decision-making and project planning.
The accuracy of selected library models has been shown with excerpts from completed and ongoing validation case studies, including measurement data of real TES, and cross-comparison studies with other storage models.Further studies in this regard are underway.
To demonstrate the application of the LargeTESmtk, an exemplary simulation case study was conducted comparing the performance of pit and tank TES under different operation modes ranging from short-term to seasonal stor-age operation.As expected, the simulation case study revealed that the storage efficiency drops from above 99% for the short-term operation to around 58% for seasonal operation for the PTES case and that TTES generally performs better than the PTES.With the simulation study, the good applicability of the models with reasonable calculation times suitable for large parameter studies was shown.Moreover, the LargeTESmtk provides the flexibility and models to extend the study to other storage geometries or studies on system level.
The Modelica library is currently still undergoing a continuous development process.As shown in this paper, many features and models are already available, but some of the presented models and features (e.g., hybrid geometries) still need to be added and subjected to validation studies.Furthermore, it is intended to make the library available to everyone soon.So far, Dymola has been this library's main development and simulation environment.However, it is planned to test and ensure compatibility of the library with other simulation environments such as OpenModelica7 .
A t[n] and A b[n] are obtained by inserting the respective coordinates z t[n] and z b[n] in Equation 4. A similar approach can be used for other fluid domain geometries.For instance, simply choosing an angle of α = 90 • will result in a cylindrical fluid domain.The sum of all enthalpy flows for each volume element is the result of the induced volume flows during charging and discharging of the storage, where each enthalpy flow follows the basic equation Ḣ = V • ρ • c p • T , with the corresponding volume flow V and temperature T :

Cover/ground surface Heat transfer: Cover, walls, bottom Diffusers Ground properties Flexible ground layers Uniform ground Pure thermal resistance Thermal resistance & capacitance Convection and radiation Pure convection Variable number and position Cylinder/Cuboid Cone/Pyramid Hybrid geometries Free-standing Partly buried Fully buried
Figure 1.Overview and features of the LargeTESmtk.

Table 1 .
Investigated operation modes of the simulation case study.

Table 2 .
Figure 6.Annual storage efficiencies ranging from short-term to seasonal storage operation.Applied model parameters of the simulation case study.