Novel District Heating Systems: Methods and Simulation Results ^†

Abstract

Fifth-generation district heating and cooling (5th DHC) systems offer promising approaches to decarbonizing space heating, cooling and domestic hot water supply. By using these systems, clustered buildings combined with industrial waste heat can achieve a net-zero energy balance on a variety of time scales. Thanks to the low exergy approach, these systems are highly efficient. As part of the Smart Anergy Quarter Baden (SANBA) project, the thermal energy grid simulation tool TEGSim has been further developed and used to design an ultra-low-temperature district heating (ULTDH) network with hydraulic and thermal components fitted to the specific regional characteristics of the investigated case. Borehole thermal energy storage (BTES) used as seasonal storage ensures long-term feasibility. The annual discrepancy of input of thermal energy provided by space cooling and output of energy demanded by space heating and domestic hot water is supplied by an external low-grade industrial waste heat source. This paper presents the functionality of the simulation and shows how to interpret the findings concerning the design of all components and their interplay, energy consumption and efficiencies.

1. Introduction

In Austria in 2019, 27% (305 PJ) of overall energy consumption was used for space heating (H), cooling (C) and domestic hot water (DHW) supply, 198 PJ of that in private households alone. Approximately 46% of all people in private households live in multi-family buildings. The others live in single-family houses or semi-detached houses. Thus, about 40% (79 PJ) of the energy required for H, C and DHW is supplied by fossil fuels and about 36% (71 PJ) by ${\mathrm{CO}}_2$ neutral resources. The remaining 24% (48 PJ) are provided by electricity or district heating (DH) [1]. Depending on local conditions, the required electricity or energy used in district heating systems can be generated either in a climate-neutral manner or from fossil fuels. The Austrian government’s ambitious goal is to free H and DHW supply from coal- and oil-based fuels by 2035 and to eliminate fossil gas-based fuels by 2040 [2].

District heating systems have been used in Europe since the 14th century [3]. These systems distribute thermal energy for H and DHW between sources and sinks through liquid or steam filled pipes [4]. The required heat is often generated in a small number of combined heat and power plants (CHP). Although fossil fuel-fired power plants are highly efficient, they are clearly not a long-term option in reducing greenhouse gases [5]. In [6], Lund et al. categorize these systems as 2nd or 3rd generation DH systems, depending on whether other climate-neutral technologies such as biomass and large-scale solar power plants feed energy into the grid or not. However, due to the high temperature level of the heat supply and the associated high exergy content, such systems are still subject to high energy losses. Energy losses can be further reduced by lowering the network temperatures [6], which in turn means low-grade industrial waste heat sources and renewable energy sources can be used [7]. However, the increasing demand for cooling energy in the coming decades cannot be met with these systems [8]. New approaches include the so-called 5th Generation District Heating and Cooling (5th DHC) networks [9]. If the temperatures of the pipe network are close to the ambient temperature, 5th DHC can also be described as ULTDH networks. Despite the complexity of such smart systems, they are as yet in the early stages of development. Here, consumers and producers can extract energy for H and C from the connected 2-pipe network. From a thermodynamic perspective, cooling energy is extracted by supplying energy to the pipe network. Network participants who can both supply and extract energy are referred to as prosumers (producer + consumer). The difference in temperature between the hot and cold pipes is only about 5–10 K. The temperature level of the network is 5–30 ${}^{\circ }$C, which is very close to the ambient temperature. Highly efficient heat pumps are used to provide the required temperature level for H and DHW. The required energy is extracted from the warm pipe. In addition, energy for cooling can be extracted directly from the cold pipe and buildings can be cooled using free cooling (FC). In areas where the annual heating and cooling energy demand is approximately the same magnitude, seasonal storages can further reduce the additional energy demand from low-grade industrial waste heat.

Most sites in operation can be found in Germany and Switzerland, e.g., ETH Hönggerberg [9,10]. Currently, the largest 5th DHC network in Austria is the Viertel Zwei Plus. The dynamic interaction of all network components, i.e., heat pumps, BTES, and the ULTDH network, is crucial for the technical design [11]. In order to further encourage development, and make 5th DHC systems a competitive technology, design and simulation tools are in great demand.

1.1. State of the Art Tools

TRNSYS is a commercially available software package that has become a standard tool for simulating thermal and electrical energy systems. The software can be used for detailed thermal and electrical simulations of, for example, solar collectors, heat pumps and fuel cells. However, TRNSYS is not suitable for simulating district heating networks at the district or city level [12]. Furthermore, TRNSYS is dependent on the software libraries available for modelling the individual components.

CFD software, such as Ansys Fluent or OpenFOAM, can obtain detailed information about the pressure and temperature field in a pipe network. Although CFD tools represent the pipe network in great detail, they cannot model geothermal probes or heat pumps in large scale applications. The thermal and hydraulic behavior of buildings or geothermal probes is estimated by using boundary conditions. However, the high numerical effort often leads to unacceptably long calculation times.

Other software is mainly limited to the simulation of energy flows of geothermal probes or probe fields. The Earth Energy Designer (EED) can calculate analytical solutions for single probes or probe fields with up to 1200 probes. With the help of numerically generated g-function, different probe geometries, distances and configurations can be considered [13,14]. With FEFLOW [15], detailed numerical simulations based on finite elements of the underground can be carried out. In contrast to TEGSim, FEFLOW can simulate aquifer thermal energy storages and convective energy transports in the ground. Detailed underground simulations can also be carried out with Comsol. However, the complex interaction between the BTES and the energy distribution network cannot be taken into account.

De Carli et al. [16], for example, used TRNSYS to simulate the buildings and EED to simulate the BTES for a ULTDH network in Italy. Further comparisons between the currently available simulation tools can be found in [12,17,18].

1.2. Scope and Limitaions

In this paper, the developed simulation tool TEGSim and the simulation results for a use case are presented. The investigated use case is discussed only with respect to thermal and hydraulic aspects. An economic analysis is not presented. In contrast to commercially available CFD software, TEGSim can also simulate components such as heat pumps or BTES. Highly accurate statements can be made about the temperature field in the pipe network and the BTES with a simultaneously low numerical effort using two-and-a-half-dimensional resistance and capacity models (RCMs). Due to the reduced model complexity, the simulation of widely branched and meshed pipe networks is possible with computation times of only a few hours. The thermal calculation of the pipe network and the BTES is carried out entirely transiently without analytical g-functions. The geometry and configuration of the pipes can therefore be adapted as required. The hydraulic calculation is quasi-stationary with the help of relaxation.

In contrast to building simulation tools such as TRNSYS, no building simulations are possible. With TEGSim, it is possible to size the essential components such as heat pumps, daytime storage tanks, BTES, pipe network and circulation pumps. With regard to the complex interaction between BTES, the pipe network and prosumer, TEGSim can also size the BTES and the pipe network. This holistic approach allows the energy system to be modelled, dimensioned, simulated and analysed with the relevant components. The developed simulation tool can perform reliable simulations of any district heating and cooling networks of any size with different discretization, parameters and components.

2. Methodology

In the course of the national project SANBA (see Section 4), the DHC network simulation tool from Nagler [19] was further developed and programmed using MATLAB. The developed thermal and hydraulic simulation program TEGSim is a stand-alone tool that requires no external calculation software for the individual components apart from the programming language and the buildings’ load curves. The developed algorithm aims to make detailed predictions about the thermal and hydraulic performance of the essential elements within DHC networks over seasonal periods, in our case one year. The key elements are the pipe network, the energy transfer stations (ETS) and the BTES. The pipe network is modelled as a closed-loop, two-pipe system with the possibility of bidirectional mass and energy flow. At the ETS, prosumers can be provided with energy for H, DHW and C. Depending on the desired operation, the BTES can compensate for the annual difference between heating and cooling energy.

The necessary input to the simulation tool is the predefined topology of the pipe network and the prosumers’ heat load curves. If, as in most cases, additional energy is required to balance the annual discrepancy between heating and cooling demand, the load curve of the heat source must be given.

In the following chapters, the modelling approaches used for the main components are discussed in more detail.

2.1. Thermohydraulic Simulation

The network calculation aims to determine the thermal and hydraulic operating state at each time step. Hydraulic variables such as pressure and mass flow in the pipes, temperature and thermal losses in the network nodes are calculated. Due to the complexity of the energy distribution grid, a thermohydraulic calculation method was designed to keep the calculation time and data storage capacity as low as possible. The calculation can therefore be divided into a quasi-stationary hydraulic calculation and a transient thermal calculation.

The variables resulting from the iterative hydraulic simulation, such as pressure and mass flow, serve as boundary conditions for the following transient thermal calculation. For the selected simulation time step of one hour, the hydraulic quantities are kept constant. In the transient thermal simulation, the main task is to calculate all the temperatures in the pipe network, the BTES and the ground. All node temperatures and boundary conditions are combined in a linear system of equations and calculated simultaneously. Therefore, the BTES and pipe network temperatures of all nodes are known for every simulation time step. The simulation time step for the thermal simulation is set to 30 min per default. After the thermal calculation has been completed, the next boundary conditions are set for the hydraulic simulation. Figure 1 shows the flowchart illustrating the thermohydraulic approach for a single hydraulic time step and a user-defined number of thermal calculation steps.

The program was designed to simulate annual behavior and show the proposed concept’s long-term viability. To obtain reliable simulation results, it is necessary to consider the full range of time-dependent load profiles and boundary conditions.

2.2. Pipe Network

In this section, the applied methods for thermal as well as for hydraulic calculation are presented.

2.2.1. Hydraulic Calculation

TEGSim automatically estimates the pipe diameters in a first step according to the maximum thermal power to be transported. An adjustable maximum pressure loss per meter pipe is an efficient approach reducing the pumping work and investment costs. According to Bothe [20], the ideal pressure losses per meter pipe are in the range of 80 ^Pa/_m to 200 ^Pa/_m. Otherwise, the ideal pressure loss per meter pipe can be obtained from the manufacturer. However, since the pipes are designed for the maximum load, a manual adaptation of the pipe diameters based on technical and economic criteria such as price per diameter, installation costs and efficient flow rate is unavoidable.

The quasi-stationary hydraulic calculation was applied following existing methods and the work of Bothe [19] and Nagler [20]. The topology of the network is represented by an incidence matrix adapted from graph theory [21]. The system equations are set up according to the QH method, wherein an equation is set up for each pipe string as a function of the volume flow in the pipes Q and the pressure heads in the nodes H. However, the equations are formulated with variables of pressure and velocity. This eliminates the need for subsequent conversion of the simulation results. The Darcy–Weisbach equation is used for the friction pressure losses ${\Delta }{p}_{\mathrm{f}}$. These pipe friction losses depend on the pipe friction coefficient $\lambda$, the density of the fluid $\rho$, the pipe length L, the characteristic pipe diameter d and the flow velocity in the pipes u. Geodetic height differences ${\Delta }h$ and the resulting pressure differences $\rho g{\Delta }h$ are taken into account by extending the Darcy–Weisbach equation according to Equation (1)

$$p_1-p_2=\frac{\lambda \rho L}{2d}{u}^2\pm \rho g{\Delta }h={\Delta }{p}_{\mathrm{f}}\pm \rho g{\Delta }h.$$

(1)

The pipe string equation is computed using the linearization method according to Wood and Charles [22]. The non-linear part of Equation (1) is divided into the constant part $\frac{\lambda \rho L}{2d}\left(u^{t-1}\right)$ and the linear part $u^t$ according to Equation (2)

$${\Delta }{p}_{\mathrm{f}}=\frac{\lambda \rho L}{2d}\left(u^{t-1}\right)u^t.$$

(2)

The superscript index $t-1$ refers to the previous and the superscript index t to the current iteration step.

Equation (2), as well as the mass balance are combined into a linear system of equations and repeatedly solved under consideration of a relaxation step until the user-defined convergence criterion is fulfilled. The relaxation should prevent possible oscillations and ensure convergence.

2.2.2. Thermal Calculation

To model the thermal behavior of the pipe network, a 2D RCM for 1U-probes according to [23] is used, see Figure 2. The thermal conductivity and capacity of each of the different layers and materials of the probe are considered. Inside the pipes, the heat transfer coefficients are calculated according to the flow regime and the properties of the heat transfer medium, which is assumed to be incompressible. The calculation of thermal resistances within the probe is partly with empirical approaches from [23] and using the 1st order multipole method from Claesson and Javed [24]. The filling material is divided into two equal parts, each with one node in which the temperature ($T_{\mathrm{f}1}$, $T_{\mathrm{f}2}$) is calculated. The node of the borehole wall is connected to the ground using thermal resistances and capacities. The surrounding ground is modelled axially symmetrically with n concentric circular ring-shaped surfaces coupled to the borehole wall using thermal resistances and capacities. In the use case presented in Section 4, the pipes between the borehole wall have been subdivided by the outer diameter of $2.4$ m into nine circular rings of equal thickness. From a programming point of view, it makes sense to discretize each pipe with the same number of elements. The discretization was chosen in such a way that the maximum axial cell thickness does not exceed 1 m. Therefore, each pipe has been divided into 225 equidistant elements.

Based on the general transport equation, the nodes of the individual layers can be connected to a two-and-a-half-dimensional model of the probe and the surrounding ground. Within the liquid-filled pipes, heat transfer by conduction and convection is considered. In contrast, in the solid components such as the pipe, filling material and ground, heat conduction is only possible through conduction. The thermal model is rendered fully transient by using the implicit formulations of the 1st order forward differences in time and the 1st order upwind scheme in space. Axial and radial discretization was evaluated and optimized according to calculation time and accuracy. More details on the transport equations can be found in [19,23,25].

2.3. BTES

The BTES consists of several parallel or seriell connected probes. TEGSim is able to simulate CX-, 1U- and 2U-probes. In the following, only the 2U-probe will be discussed in more detail. The essential quantities of the BTES, such as the temperature field and pressure loss, are partly calculated analogously to those of the pipe network. The methods used for the hydraulic and thermal calculation are outlined more precisely below.

2.3.1. Hydraulic Calculation

The hydraulic calculation of the BTES essentially aims to determine the pressure losses. Since all probes within the BTES undergo the same mass flow, it is sufficient to calculate the pressure loss for only one probe. Analogous to the pipe network, the flow inside the BTES is also considered incompressible. The friction-induced pressure loss between the inlet and outlet of the probe is determined using Equation (3) according to Darcy–Weisbach.

$${\Delta }{p}_{\mathrm{probe}}=\frac{\rho u^2}{2}\left(\frac{\lambda L}{d}+\zeta \right).$$

(3)

The pressure loss due to directional flow change at the base of the probe with the length L is considered by using the discharge coefficient $\zeta$. The calculation of the pipe friction coefficient $\lambda$ depends on the flow regime. For laminar flows ($Re<2300$), the Hagen–Poiseuille equation is used. Iterative estimation of the pipe friction coefficient is necessary for flows in the transition region and the fully turbulent region. The Colebrook–White equation is used in the transition region and the Nikuradse equation in the fully turbulent region.

2.3.2. Thermal Calculation

The thermal calculation of the BTES is largely analogous to the pipe network. In contrast to the pipe network, however, different probes can be used. Therefore the thermal behavior of the BTES, 2D RCMs for CX-, 1U- and 2U-probes according to [23,26] were used. Furthermore, different geometry and operating modes of the BTES can be selected, e.g., 1U- or 2U-probes with diagonal flow, adjoining flow with or without Us crossing. In contrast to Nagler’s work in [19], thermophysical parameters such as density, heat capacity and thermal conductivity of the ground can be taken into account depending on depth. The model presented can thus be better adapted to local conditions. Figure 3 shows the RCM used for a 2U-probe.

The heat conduction in the ground occurs in a radial and vertical direction. The model for a 1U probe presented by Nagler in [19] takes vertical heat conduction within the transport equations into account, but the ground underneath the probe is not modelled. With decreasing probe length, the influence of the modelled ground underneath the probe increases. We therefore considered the ground below the probe in the development of TEGSim. For this purpose, the borehole was extended to the lower edge of the calculation area and connected to the surrounding ground with one node per layer. The thermophysical parameters of the extended borehole were also adapted to those of the ground. Using this procedure, the thermal behavior within the area of influence of the probes can be calculated more precisely.

The individual probes in the BTES are aligned equidistantly in x- and y-direction. The distance between the probes is typically in the range of a few meters. The neighboring probes thermally influence the probes in the centre of the BTES. By contrast, the probes at the edge are only partially influenced by the adjacent probes. A superposition of the calculation results from one probe to all other probes would violate the model assumption of identical boundary conditions. As shown in Figure 4, probes with the same number behave thermally identically. To keep the numerical effort low, it makes sense to calculate only these probes. In the case of the BTES shown with $5\times 5=25$ probes, the calculation of only 6 probes is necessary. The exact number of probes can be adjusted by removing those at the corners. For a BTES with, for example, 21 probes, the probes with the number 6 can be deleted from Figure 4.

The reciprocal influence of the quadratic probe field is taken into account with the help of Kelvin’s line source theory [27] according to Equation (4)

$$T(\tilde{r},\tilde{t})=\frac{\dot{q}}{4\pi \lambda }{\int }_{\frac{{\tilde{r}}^2}{4a\tilde{t}}}^{\infty }\frac{{\mathrm{e}}^{-u}}{u}du.$$

(4)

The variable $\dot{q}$ describes the radial heat flow over the borehole wall, $\lambda$ the thermal conductivity of the ground, $\tilde{r}=\frac{r}{r_b}$ the dimensionless radius with respect to the radius of the borehole wall $r_b$ and $\tilde{t}=Fo=\frac{at}{r_b^2}$ the Fourier number. In order to keep the numerical effort low despite the high quality of the results, the adjustment of the lateral temperature boundary condition of the individual probes is not carried out in every time step. Lamarche has shown in [28] that the deviation of the solution in Kelvin line source theory is small for $Fo>20$. The adjustment of the temperature boundary condition with the help of Equation (4) is therefore only carried out if the condition in Equation (5) is fulfilled

$$t_n-t_{n-1}>20\frac{r_b^{2}c\rho }{\lambda }.$$

(5)

In this case, the specific heat capacity c, the density $\rho$ and the thermal conductivity $\lambda$ refer to the surrounding ground. Through exploratory drillings at the site of the investigated use case (see Section 4), the following mean thermophysical ground parameters could be determined: $\rho =1800$ kg/m${}^3$, $\lambda =1.9$ W/m/K and $c=800$ J/kg/K. This period for this condition $t_n-t_{n-1}$ with the selected thermophysical parameters and a borehole radius of $r_b=0.075$ m is calculated as 23.68 h. In order to ensure time synchronization with the time step size of the thermal and hydraulic calculation, the temperature at the lateral boundary layer of the simulated probes is adjusted every 24 h or every 48 thermal calculation steps with a time step size of 30 min. Due to the linearity of the fundamental heat conduction equation, the temperature response of the m probes can be superpositioned according to Equation (6)

$${\Delta }T(x,y,\tilde{t})=\sum _{i=1}^m\frac{{\dot{q}}_i}{4\pi \lambda }{\int }_{\frac{{\widetilde{r_i}}^2}{4a\tilde{t}}}^{\infty }\frac{{\mathrm{e}}^{-u}}{u}du.$$

(6)

With the heat flows ${\dot{q}}_i$ of the individual probes and the thermophysical parameters of the ground, the temperature field $T(x,y,\tilde{t})$ in the entire calculation area can be determined.

2.4. Boundary Conditions

The boundary conditions of the numerical simulation can be divided into two groups: hydraulic and thermal boundary conditions.

An essential hydraulic boundary condition is the mass flow taken from or fed into the pipe network by the ETS. This mass flow is mainly determined by the load curves of the buildings and the filling levels of the thermal storages in the ETS. To fully define the mass balance of the pipe network, either the mass flow at the waste heat source or at the entry into the BTES must be specified. If several waste heat sources are connected to the pipe network, it must always be ensured that the mass flow is not defined at exactly one feed-in or extraction point. At this point, the mass balance is automatically closed based on the complex interaction of prosumer, waste heat source and the BTES.

The thermal boundary conditions of the BTES and the pipe network are formulated as Dirichlet boundary conditions. Figure 5 shows three possible options. On the one hand, a constant temperature, independent of the depth, can be imposed at the lateral edge of the calculation area (Figure 5a). On the other hand, the geothermal temperature gradient can be taken into account with and without surface weather influence (Figure 5b,c). The geothermal temperature gradient can be modeled from any number of continuous linear functions. In most cases, the geothermal temperature gradient is about 0.03 ^K/_m [29]. The near-surface weather influence and the associated hourly temperature changes extend to a depth of about 8 to 18 m, depending on the thermal conductivity of the ground [29]. The near-surface weather influence is modeled using a test reference year (TRY) dataset. The TRY dataset available to us provides hourly values of ground temperatures at depths of 0.1, 0.2 and 0.5 m. Temperatures to a depth of 0.5 m were linearly interpolated. Between 0.5 m depth and the depth of the undisturbed ground temperature, the temperatures were calculated using the one-dimensional heat conduction equation. With the geological conditions of the area we studied, the undisturbed ground temperature occurs at a depth of 15 m and is always 10 ${}^{\circ }$C throughout the year. With the help of Neumann boundary conditions at the top and Dirichlet boundary conditions at the bottom, a typical temperature curve for the topography was determined.

Due to the rotational symmetry of the probe models used, the distance between two probes in a probe field is equal to the diameter of the lateral boundary condition. Since the probe spacing is usually bound to technical or topographical conditions, the spacing of the lateral boundary condition usually cannot be freely selected. Such limitations of the model lead to inaccurate representation of the ground around the probe field. In order not to violate the Dirichlet boundary condition of the external probes, the heat flow at the edge of the computational domain must approach zero. As the ratio of volume to surface area of the geothermal energy storage increases, the influence of the violated boundary condition decreases. Therefore, increasing probe spacing together with more probes leads to smaller errors.

2.5. Energy Transfer Stations

At the ETS, prosumers can be provided with energy for H and DHW via heat pump, FC directly from the grid and active cooling (AC) via the reversible heat pump that otherwise supplies heating energy. The ETS connect the prosumers to the ULTDH network. Depending on the prosumers, their tasks are:

• Supply of hot water for H;
• Supply of hot water for DHW;
• Supply of cold water for FC or AC.

Furthermore, these components are located in the ETS: heat pumps, heat exchangers, short term water storages, circulation pumps and valves. Since not all prosumers have the same requirements in delivering energy for DHW, H and C, different types of ETS are generally necessary. However, they all consist of the same kind of components except for different numbers of heat pumps with different sizes and short term water storage sizes. Figure 6 shows the schematic hydraulic diagram of a fully equipped ETS. Depending on the requirements, various components can be left out. For example, an office or school building with minimal demand for hot water does not need a heat pump to produce DHW. The required DHW is highly efficiently provided by the use of decentralized electric hot water boilers.

The ETS modeling, especially the modeling of grid pumps, heat exchangers and short term water storages, was deliberately kept simple to reduce complexity and calculation time. A more detailed model would only affect the simulation results to a minor degree.

For the models used here, pressure losses within the ETS only occur due to pipe friction. Due to the short pipe lengths within the ETS, we made the assumption that pressure losses only occur when fluid flows through the heat exchangers. With the help of $\zeta$ values, these can be taken into account analogously to Equation (3). Therefore, the pressure loss changes the pressure boundary condition at the transfer point between the ETS and the pipe network.

2.5.1. Heat Pumps

The temperature difference between the inlet and outlet temperature of the heat pumps was set to a user-defined value of 4 K. Due to the non-constant inlet temperature to the heat pumps, the temperature dependence of the coefficient of performance (COP) must be taken into account. On the one hand, temperature-dependent COP values can be accounted for according to Equation (7)

$$CO{P}_{\mathrm{real}}={\eta }_{\mathrm{HP}}\frac{T_{\mathrm{h}}}{T_{\mathrm{h}}-T_{\mathrm{c}}}.$$

(7)

The COP values are calculated depending on the heat pump efficiency ${\eta }_{\mathrm{HP}}$, the temperature of the heat source $T_{\mathrm{c}}$ and the temperature of the heat sink $T_{\mathrm{h}}$. The problem with this procedure is that the performance of genuine heat pumps cannot be taken into account. Therefore, on the other hand, manufacturer-specific models with given COP values and performance classes can be implemented. Thus, the temperature dependence of the COP values of the heat pumps was taken into account using COP values specified by the manufacturer for two temperature levels (e.g., 35 ${}^{\circ }$C for H and 63 ${}^{\circ }$C for DHW). TEGSim can generically determine the number of heat pumps required by specifying specific maximum and minimum loads out of a set of predefined heat pumps. However, it has been shown that the assumed source temperature for the prosumers (e.g., 15 ${}^{\circ }$C) varies only slightly. A more detailed calculation of the COP values requires more computing time, while the influence on the resulting electrical energy demand is negligible.

2.5.2. Heat Exchangers

The heat exchangers in the ETS have been modeled in such a way that a user-defined temperature difference of, for example, 5 K can always be achieved. This temperature difference can always be achieved by adjusting the flow rate at the heat exchangers depending on the load. Detailed modeling of the heat exchangers did not make sense due to long calculation times and negligible improvement in accuracy.

2.5.3. Short Term Water Storages

In a first approach, TEGSim estimates the capacity of the storage tanks in such a way that the heat pumps only have to be designed for 80% of the peak load. This guarantees an economically rational dimensioning of manufacturer-specific devices. Furthermore, TEGSim selects the storage sizes within predefined storage volumes. This makes it possible to implement off-the-shelf storages without iteratively adapting the storage size to manufacturer-specific sizes. The water storages are modeled with two temperature levels, a one-dimensional thermocline and a charging level. Evaluating the results, it was found that the water storages are frequently charged and discharged and idle charging cycle times are in the range of a few hours, which justifies the simple model.

2.5.4. Circulation Pumps

Each of the ETS has several circulation pumps depending on its configuration. The pumps are controlled so that a user-defined temperature difference of, e.g., 4 K at the heat exchangers can always be ensured. The pressure losses of the components within the ETS are taken into account using manufacturer-specific zeta values. The dimensioning of the circulation pumps is conducted in post-processing based on the pressure differences and the required mass flows.

3. Validation

The validation of the presented simulation tool is not possible with standard software due to the highly individualized ETS. Due to missing measured values, we were not able to validate the overall model. Since the ETS only define the temperature difference between feed and return pipe and the mass flow, it is sufficient to validate only the pipe network and BTES. The methodology of the hydraulic pipe network calculation was validated by Bothe in [20] with the commercial software PSS^®SINCAL. The methodology for the thermal pipe network calculation is analogous to that for 1U probes. Bauer validated the simulation results in [23] for one probe using measured data from a thermal response test. The validation of the BTES and the resulting reciprocal influence of the probes was carried out with the commercial software FEFLOW. Thus, the finite element program can simulate heat transport processes in the underground and model geothermal heating systems [15].

Numerical simulations were performed with TEGSim and FEFLOW using identical and simplified initial and boundary conditions. A homogenous subsurface model with a depth-averaged effective thermal conductivity of 1.75 W/(mK) and a constant temperature at the lateral boundary of 15 ${}^{\circ }$C were chosen. Within the BTES there are 349 2U-probes with a distance of 4 m in x- and y-direction, see Section 4. The water-carrying pipes are connected in parallel and have a dimension of 40 × 3.7 mm corresponding to current standards. The borehole diameter was chosen as 150 mm. Furthermore, the same simplified input data, mass flow of the fluid and the feed temperature into the BTES were used. The test function of the mass flow rate is shown in Figure 7. The mass flow was chosen so that the flow within the pipes is always fully turbulent. The corresponding Re numbers are $Re=1.03\times {10}^4$ and $Re=1.72\times {10}^4$.

Figure 8 shows the inlet temperature into the BTES and the resulting mean outlet temperatures. The inlet temperature and the mass flow rate were defined in such a way that each of the three different inlet temperatures correlates with two different mass flow rates. The profile of the inlet temperature was defined according to the actually occurring conditions. Thus, six different operating conditions can be simulated. Throughout the validation process, the resulting mean outlet temperatures of the BTES were compared for each time step. The values of the thermal resistances of the FEFLOW model were adjusted to keep the mean annual deviation at less than 0.65%.

The relative deviation for the FEFLOW simulation is shown in Figure 9. More significant deviations can only be observed when input values change abruptly. The highest resulting local deviations are between 10 and 60%, but they decrease again to less than 1% after only a few minutes. A possible reason is the different time discretization for solving the systems of both models. TEGSim uses constant time intervals of one hour, whereas FEFLOW uses an automatic time-step control scheme, a predictor–corrector scheme.The time-step control scheme ensures that after long periods with constant entry conditions, the abrupt changes of these conditions are detected with a slight delay. A few seconds after sudden changes, the time step width decreases strongly until it becomes longer again when the entry conditions remain the same. Nevertheless, the minor overall deviations show a good correspondence of the simulation results and the validation’s success.

4. Use Case

In the national project SANBA, a ULTDH grid was investigated to be applied at the former military base Martinek–Kaserne using the TEGSim program. This area of 40 hectares in the city of Baden south of Vienna was abandoned in 2014. The consortium developed three different urban mix-use scenarios: MINI, MIDI and MAXI [30]. In scenario MINI, only existing buildings with a total gross floor area (GFA) of 65,591 m${}^2$ are renovated. In scenario MIDI, the existing buildings will be renovated. Besides this, apartments, a school and a supermarket with a total GFA of 94,000 m${}^2$ are considered. In Scenrio MAXI, in addition to the renovated existing buildings, further new buildings with a total GFA of 142,000 m${}^2$ are considered. Both the existing buildings and the new buildings were divided into different usage categories. Depending on the category, corresponding requirements for the provision of energy for H, C and DHW result.

Table 1. Classification of usage and whether energy is needed for H, C and DHW. Categories marked with * are not part of the presented MIDI Scenario.

Usage	DHW	H	C
office		x	x
education		x	x
archive *		x
residential	x	x	x
gastro, event, shop		x	x
supermarket *		x	x
hotel *	x	x	x

summarizes the categories used for the three scenarios. For example, the usage archive, which occurs only in scenario MINI, requires energy only for H to keep buildings non-freezing in winter. In contrast, the usage residential and hotel require energy for H, C, and DHW. All usage categories without DHW supply are equipped with decentralized electric water boilers.

In this paper, the use of TEGSim is based on one of the investigated scenarios–the so-called MIDI scenario, which was estimated as very likely be realized, are discussed. The findings from this scenario show the methods and results obtained by TEGSim. Figure 10 shows the area of interest of scenario MIDI from the project SANBA. The black lines represent the annular thermal energy grid itself. The black circles each represent a prosumer and their ETS. The low-grade industrial waste heat is provided by a neighboring dairy plant.

4.1. System Design and Heat Input

The uninsulated pipes in the pipe network are made of polyethylene and installed in the ground at a depth of 1.4 m. Due to the low fluid temperatures close to the ambient temperature in the pipes, insulation of the pipes is not necessary from a technical point of view and not desireable from an economic perspective. In total, 4458 m of supply and return pipes (3031 m ring line, 1427 m branch line) are considered for the ULTDH network. The dimension of the ring line is 280 × 25.4 mm. The dimensions of the branch lines depend on the load and range between 160 × 14.6 mm and 280 × 25.4 mm. Due to the high pipe volume of 157 m${}^3$ and moderate minimum temperatures inside the network way above 0 ${}^{\circ }$C, water is used as the heat transfer medium.

Table 2. Categorisation of prosumers and equipment of ETS: Gross floor area, capacity and number of heat pumps for H/AC and DHW, volume of short term water storages.

Prosumer	Usage	GFA	${\mathit{P}}^{\mathbf{H}/\mathbf{AC}}$	${\mathit{n}}^{\mathbf{H}/\mathbf{AC}}$	${\mathit{V}}^{\mathbf{H}/\mathbf{AC}}$	${\mathit{P}}^{\mathbf{DHW}}$	${\mathit{n}}^{\mathbf{DHW}}$	${\mathit{V}}^{\mathbf{DHW}}$
-	-	m${}^2$	kW	-	m${}^3$	kW	-	m${}^3$
P1	education	3383	122.5	1	5	76.4	1	1
P2	residential	30,000	81.5	6	15	198.2	1	1
P3	education	13,716	152.8	2	9	120.1	1	1
P4	residential	22,000	152.8	3	14	120.1	1	1
P5	residential	15,000	152.8	2	9	120.1	1	1
P6	education	27,000	152.8	7	26	-	-	-
P7	education	3730	122.5	2	6	-	-	-
P8	gastro, event, shop	4445	152.8	2	7	-	-	-
P9	gastro, event, shop	4445	122.5	2	6	-	-	-
P10	office	12,251	81.5	4	10	-	-	-
P11	office	12,264	122.5	3	11	-	-	-
P12	office	11,357	152.8	2	9	-	-	-

lists the different prosumers, the intended usage of the associated buildings and the corresponding GFA. Furthermore, the performance of the individual heat pumps for H/AC ($P^{\mathrm{H}/\mathrm{AC}}$) and DHW ($P^{\mathrm{DHW}}$) as well as their quantity ($n^{\mathrm{H}/\mathrm{AC},\mathrm{DHW}}$) is shown. The number of heat pumps varies depending on the load profiles. One heat pump is sufficient to meet the DHW demand in each case. Buildings that are used as schools, offices, or gastro, event, shop do not have heat pumps to provide DHW. The demand for DHW in these buildings is covered exclusively by decentralized electric hot water boilers. Overall, there are twelve ETS and more than 30 heat pumps at a cumulative peak load of more than 4.4 MW_th.

For a seasonally balanced system, a couple of parameters play a crucial role. The most critical parameters are the size of the BTES, which means the number of probes, and the low-grade waste heat supplied from the dairy plant. Optimizing these two values for a balanced BTES is always possible. Finally, it is crucial for energy efficiency to keep temperatures low in summer to ensure FC instead of AC via a heat pump. FC is only possible up to 18 ${}^{\circ }$C supply temperature of the heat exchangers. The highest return temperature for FC is set to 23 ${}^{\circ }$C, considering the plate heat exchangers’ pinch point of 5 K between the cooling surfaces in the buildings and the pipe network.

After determining the minimal probe number, a parameter study was performed to find the adjustable input parameters’ optimal values. In this scenario, 349 probes with a depth of 180 m are required in the BTES. Due to spatial conditions, a probe spacing of 4 meters in x- and y-direction was chosen. The probes are designed as parallel flow 2U probes with a dimension of 40 × 3.7 mm. The diagonal distance between the two U-tubes is 98.995 mm and the borehole diameter is 150 mm.

The evaluation of the waste heat potential showed that four main sources of waste heat are available: hot waste water, compressed air generation, refrigeration and steam generation. Since the available waste heat exceeds the required energy by far, waste heat is used exclusively from refrigeration. The measured supply temperature at the chillers connected to the refrigeration production varies according to the season between 32 and 36 ${}^{\circ }$C. The measured temperature difference between supply and return is about 4–5 K. By using ULTDH networks, even this temperature level, which is very close to the ambient temperature, can be used. The average waste heat capacity of the chillers comes to 1692 kW. If all the available low-grade waste heat is fed in, the BTES and pipe network temperature will increase to such an extent that FC is no longer possible. Therefore, only a fraction of the low-grade waste heat potential from the chillers is required. However, by feeding the waste heat into the ULTDH network, the dairy plant cannot eliminate the chillers. On the one hand, because not all the available energy can be used, and on the other hand, because the continuous production of dairy products requires high reliability of the cooling systems.

Figure 11 shows the actual power taken from the dairy plant. During the cooling season, approximately from hour 3400 to 6200, the available power taken is considerably lower. As a result, the temperature in the BTES and the pipe network can be reduced to such an extent that mainly FC can be used.

4.2. Energy Demand

A distinction is made between electrical and thermal energy demand. The thermal energy demand is subdivided according to DHW, H and C. These fluctuating energy demands serve as input for TEGSim. They can be created with any building simulation program. In this case, the software EDSL Tas was used. The load profiles of the buildings were calculated with a time step width of one hour. Figure 12 shows the annual energy demand for C, H and DHW. The thermal energy demand for H is $6223.5$ MWh_th/a, for DHW it is $1842.6$ MWh_th/a and for C it is $-2814$ MWh_th/a. In contrast to the demand for H and C, the demand for DHW is almost constant throughout the year. In addition, the sorted annual load curves (ALC) for H and C are plotted in Figure 12 (dotted lines). The absolute value of the slope of the ALC for the cooling demand decreases faster at high values than the slope for the heating demand. The reason for this is that in contrast to energy for H, energy for C must be provided much faster to achieve the necessary living comfort and prevent overheating on hot summer days. Therefore, the maximum cooling power (2.89 MW_th) is approximately the same as the maximum power needed for H (2.97 MW_th). Contrary to the other thermal energy demands, energy for C has to be supplied by charging the seasonal storage.

The electricity demand is subdivided into two parts. On the one hand, the electricity demand for users and residents of the buildings was calculated using standard load profiles according to [31] at 9.44 GWh_el/a. This electricity demand does not include the electricity required to operate the ULTHC network and to power the heat pumps. This energy demand is only intended to provide a rough estimation and will not be discussed any further. On the other hand, the electricity demand for the operation of the ULTDH network (circulation pumps and heat pumps) results from the simulation of the holistic system. The electricity demand of the circulation pumps is linearly interpolated based on five operating points of the pump maps. The database is made up of discrete operating states of different pumps from a selected manufacturer. Based on the maximum pump head and the mass flow, the simulation determines the electrical power demand for each calculated time step.

5. Results and Discussion

5.1. Dynamic Behavior BTES

Since the simulation calculates temperatures, pressures and velocities at every node, much information can be found in the trends apart from highest and lowest values. For example, Figure 13 shows the feed temperature into the BTES in red and the mean outlet temperature in black for one year. From hour 1 to 1400, it can be seen that the outlet temperature is higher than the feed temperature, which means discharging during winter. The opposite is true in summer when the BTES is charged. It is easy to see that the mean temperature rises in summer and reaches its highest value. In winter, the BTES cools down considerably. Since the cooling load dominates in summer and the excess heat from the buildings is stored in the BTES, the low-grade industrial waste heat can therefore be reduced to a minimum. The horizontal line in Figure 13 at a temperature of 23 ${}^{\circ }$C represents the boundary between FC and AC. If the return temperature to the pipe network is less than 23 ${}^{\circ }$C, the required cooling demand can be covered highly efficiently by FC. At higher return temperatures, the heat can still be extracted from the rooms to be cooled, but the upper temperature limit of 23 ${}^{\circ }$C ensures the necessary living comfort.

The BTES is charged a few times during the heating period at hour 750, 1400, 7200 and 7800, as shown in a positive sign of the mass flow in Figure 14. The changing sign of the mass flow is due to the complex interplay of low heat demands and high amounts of energy fed into the system by the supplier. The transitional seasons of spring and autumn are marked by frequent changes in charging/discharging and varying feed temperatures. The highest mass flows are needed in summer when cooling demand and grid temperatures are high.

The aggregated energy of the BTES shows a smooth trend near the ideal sinus curve, see Figure 15. A negative slope means that the BTES is discharged and a positive slope means that it is charged. The BTES is therefore charged during the cooling season between hours 1800 and 6800 and discharged during the rest of the time in the heating season. The BTES can only be used as seasonal thermal energy storage over a long period if the same amount of energy is discharged during the heating season as is charged during the cooling seasonm which is well met in the use case.

5.2. Share of AC and FC

A maximum permissible indoor temperature of 26 ${}^{\circ }$C was assumed when simulating the cooling demand of the buildings. In order to dampen any load peaks, the maximum permissible indoor temperature may be exceeded during 8 h per year.

Table 3. Prosumer cooling demand, duration of cooling period and whether it is provided by AC or FC.

Prosumer	${\mathit{E}}^{\mathit{C}}$	${\mathit{e}}^{\mathit{C}}$	${\mathit{p}}^{\mathit{C}}$	${\mathit{t}}^{\mathit{FC}}$	${\mathit{t}}^{\mathit{AC}}$	${\mathit{t}}^{\mathit{C}}$
-	MWhth/a	kWhth/a/m${}^2$	%	h	h	h
P1	21.33	6.30	0.8	1401	1147	2548
P2	466.01	15.53	16.6	1047	1955	3002
P3	71.44	5.21	2.5	606	1877	2483
P4	581.63	26.44	20.7	2365	1543	3908
P5	219.18	14.61	7.8	841	2146	2987
P6	734.90	27.22	26.1	1747	241	1988
P7	72.75	19.50	2.6	1196	195	1391
P8	177.27	39.88	6.3	2346	163	2509
P9	155.48	34.98	5.5	2343	164	2507
P10	78.74	6.43	2.8	1204	301	1505
P11	177.27	14.45	6.3	1192	312	1504
P12	58.01	5.11	2.1	1238	258	1496

provides an overview of prosumer cooling requirements. Compared to the existing buildings with no more than two floors, the newly constructed buildings (P2, P4, P5 and P6) are in some cases up to five floors high. The large GFA are reflected in the high cooling demand. The proposed newly constructed buildings in Scnario MIDI alone account for $71.2$% of the total cooling demand. The specific cooling demand depends, amongst other factors, on the geographical location, the building structure, the ratio of window area to façade area and the use of the buildings. In the category gastro, event, shop, the specific cooling demand is on average $37.43$ kWh/a/m${}^2$ the highest. Specific cooling demand is usually higher in new buildings than in existing buildings. The reason for this is the higher ratio of window to façade area in new buildings from scenario MIDI, 10% compared to 30%.

Whether a prosumer is cooled with AC or with FC depends exclusively on the supply temperature in the pipe network. As shown in Figure 16, the prosumers P1 to P5 are supplied with FC up to a maximum of 60%. In contrast, prosumers P6 to P12 are primarily cooled by FC to more than 80%. The reason for this is, on the one hand, the generally longer cooling period of prosumers P1 to P5, see Table 3. Towards the end of the cooling season, the BTES is already too warm to be able to cool buildings with FC. On the other hand, in the area where prosumers with very high cooling demands (P2, P4, P5) are located, the water in the pipe system is heated rapidly. Therefore, prosumers with small cooling demands (P1, P3) in the vicinity of prosumers with high cooling demands can no longer be sufficiently cooled with FC, see Figure 10. From this, it can be concluded that prosumers with small cooling loads next to prosumers with high cooling loads can be cooled less effectively with FC. This often results in financial disadvantages for prosumers with small cooling loads due to higher electricity costs for the reversible operation of heat pumps. In contrast, however, the position of the prosumer cannot be changed in most cases. Furthermore, no business models for the operation of ULTDH networks have been developed in the context of this work.

5.3. Energy Flows

Together with other parameters, such as configuration and number of ETS, network topology, size of the BTES and waste heat costs, these results allow an economic evaluation of the investigated energy supply system. Figure 17 shows the electrical and thermal energy flows of the investigated use case in daily resolution. Negative bars indicate the demand flows of the overall system and include the energy required for H, DHW and BTES charging. Positive bars indicate the supply flows of energy by discharging the BTES, feeding the low-grade waste heat, C and the electrical energy for the heat pumps and circulation pumps. The magnitudes of the supply and demand flows are identical except for the thermal losses. It is easy to see that the energy required during the heating season is mainly covered by discharging the BTES and by the additional energy fed in by the supplier. During the cooling period, the required cooling energy is almost exclusively provided by charging the BTES. The energy supplied by the dairy plant is reduced to a minimum during this period. Furthermore, it is easy to see that the electrical energy required increases on days with high cooling loads by switching on the reversible heat pumps.

Figure 18 shows the annual demand and supply flows to operate the case-specific ULTDH network. The left side shows the supply flows and how the required energy for H, DHW, BTES charge and the thermal losses is covered. The right side shows the demand flows and how the supplied energy is used. The sum of the energy flows of the left and the right side are identical and amount to 11,154.2 MWh/a each. The energy balance of the whole system is thus fulfilled in any case. In contrast to Figure 17, the electrical energy is divided between the heat pumps for H, DHW, C and the demand of the circulation pumps. Around $2795.8$ MWh_th/a is supplied to the grid by the waste heat producer. The thermal energy required for H and DHW ($8066.1$ MWh_th/a) is more than twice the amount of the required energy for C (2814 MWh_th/a). The electricity required to operate the ULTDH comes down to $2507.9$ MWh_el/a representing 22.5% of total demand/supply flows.

Since the BTES is part of the ground, no specific storage capacity can be assigned. Similar to the geothermal storage in private households with only one or two probes, it is theoretically possible to continuously discharge the storage. However, the temperature of the ground and therefore the withdrawal capacity will decrease over the years until saturation occurs. Lower ground temperatures and the resulting lower supply temperatures of the connected heat pumps decrease the COP value, see Equation (7). To keep the system efficient and to ensure the BTES long term feasibility, almost the same amount of energy must be charged and discharged periodically. As a result, the temperature of the ground remains at a reasonable level in the long term. About $3030.6$ MWh_th/a, which is a share of 27% of the total demand/supply flows, is charged and discharged within the BTES in a yearly cycle. The BTES is not entirely balanced to zero. It is discharged more than it is charged by 6 MWh_th/a.

Furthermore, it is possible to discuss the calculated heat losses and reduce them by varying pipe parameters, such as the material, diameter, or insulation. Due to the short residual times of the water in the pipe network, the thermal losses are shallow and about $51.5$ MWh_th/a, representing a share of 0.5% of total demand/supply flows. The main reason for the exceptionally low thermal losses are fluid temperatures in the pipe network that are close to the ambient temperature.

6. Conclusions

The MATLAB simulation tool TEGSim is able to give detailed information about specific use cases of ULTDH networks. In contrast to other available software, all crucial thermal and hydraulic components of the network such as BTES, pipe network and ETS are calculated simultaneously. A further adavantage is the partially holistic approach to component design. The most important models of the simulation, that of the BTES and the pipe network, were validated using FEFLOW and PSS^®SINCAL, respectively.

An additional feature is the partly holistic approach to components’ design. The results of TEGSim, exemplified in the use case of the national funded project SANBA, are the design of the grid and its components, the necessary auxiliary heat input for a balanced system, and the dynamic behavior and energy flows for the simulated time period. One significant advantage of ULTDH networks is evident in the results: the use of free cooling lowers the consumption of electrical energy significantly compared to conventional cooling technologies.

These results show that TEGSim can be used as a very effective tool for decarbonizing space heating, cooling and domestic hot water supply.