Implicit finite volume method with a posteriori limiting for transport networks

Simulating the flow of water in district heating networks requires numerical methods which are independent of the CFL condition. We develop a high order scheme for networks of advection equations allowing large time steps. With the MOOD technique, unphysical oscillations of nonsmooth solutions are avoided. In numerical tests, the applicability to real networks is shown.


Introduction
The current trend towards renewable energy sources [23] raises the importance of their efficient use. A fast and accurate simulation of their transport through the distribution network is key to the emerging optimization tasks. The most prominent examples are gas networks, the power grid and district heating networks. In district heating networks hot water is conducted from heat sources, e.g. power plants, to the consumers. As water is incompressible, the governing equation can be split into a hydraulic part, determining the flow of water, and a transport part, which models the advection of the temperature [2,17]. The hydraulic system can be solved with methods for differential algebraic systems [4,24]. In this article, we focus on constructing a suitable numerical method for the transport of the thermal energy.
For hyperbolic PDEs, a large variety of numerical methods exists [16,22]. Most of these methods are explicit, which leads to tight restrictions on the time step. In the context of district heating networks, such methods have been considered in, e.g. [9,13]. Schemes capable to handle problems with low Mach numbers [3,25] can relax the restriction due to the hydraulic wave speeds, but not those due to the velocity of the flow (≈ 1 m/s). As the lengths of the pipes in a network can vary by two orders of magnitude, short tubes can further reduce the time steps size. Local time steps on each pipe can relax this restriction [2], but bounds due to the CFL condition remain to tight for practical applications. Avoiding the time step restriction by using Lagrangian approaches, e.g. [10,18], is difficult on networks, since tracing back across the nodes is often challenging or sometimes even impossible.
Only implicit methods are completely without restrictions on the size of the time steps. Here, again central schemes can have troubles in finding suitable values at the nodes. Upwind orientated schemes further have the advantage that their systems can be solved iteratively. First order methods, like the implicit first order upwind scheme, suffer from strong numerical diffusion. A possible high order method is the implicit WENO scheme [1]. It achieves good results for moderate CFL numbers around CF L = 3 on a single segment. Applied to bigger networks, however, the arising nonlinear equations become extremely hard to solve.
In this paper, we propose a high order scheme based on a classical finite volume approach. The high order terms reduce significantly the numerical diffusion even for large time steps of CF L = 10 and above. The use of a pure upwind stencil enables a straight forward implementation of coupling conditions at the nodes and leads to an iterative solver similar to explicit schemes. It also offers the possibility to use an a posteriori limiting approach, as the values are computed in a successive order. To avoid unphysical oscillations, we use a limiting following the idea of the MOOD framework [5,6].
The paper is structured as follows: first, we introduce the mathematical model and the necessary notations. In the third chapter, schemes of 3rd and 4th order are presented and their stability is analyzed. Afterwards we define the a posteriori limiting procedure and introduce a mixed order flux by a convex combination of different orders. The numerical behavior of the scheme is analyzed with Shu's linear test and on a minimal network. Finally, we apply the scheme to a realistic district heating network and show its advantages compared to current state of the art schemes.

Model and notations
A district heating network can be described by a directed graph G = (V, E) with node set V and edge set E ⊂ V × V of ordered pairs. Following the model in [2] to each edge e, a spatially constant velocity v e is assigned. As for our numerical scheme the velocities are constant in time during a time step, we assume all velocities to be positive. This can be achieved by local transformations of edge orientations. Furthermore, we define sources S ⊂ V as a set of nodes, where Dirichlet boundary conditions can be imposed. In order to avoid some pathological configurations we further assume: The graph G is strongly connected regarding S, i.e. ∀V ∈ V : V ∈ S or ∃s ∈ S with V is reachable from s For every node V , we can identify the set of ingoing edges I V and the set of outgoing edges O V as Each edge e has a corresponding length L e and is parametrized by x ∈ [0, L e ]. As in the model in [2], on every edge e the transport of the quantity u (e.g. energy density in district heating networks) is described by an advection equation with a constant velocity v e along with suitable initial conditions. In all the nodes V ∈ V, we assume a linear mixing law of the form The first condition of (2) is the conservation of mass in the hydraulic system, while second one models a homogeneous mixing of u. At the sources Dirichlet boundary conditions at x = 0 are imposed. There are no waves entering the system at the outflow boundaries; thus, we do not prescribe any data there. In the application to district heating networks (see Section 5), the outflow nodes will be connected to consumers and their measured output acts as input data for the hydraulic solver.

Numerical method
For problems of type (1), many numerical methods are available [16,22]. For the application of district heating networks, fast schemes with time steps corresponding to CFL numbers about 20 are required, as large time horizons on networks with largely varying pipe length are considered. Furthermore, we seek for high order methods to obtain accurate results even with these large time steps. In the following, we focus on finite volume methods of order three and four and incorporate an a posteriori limiter in order to eliminate the unphysical oscillations.
When constructing schemes for conservation laws on networks, central stencils pose some difficulties. In general, the coupling conditions are not injective and therefore data in the downwind direction can not be reconstructed at the nodes. Thus, we focus on purely upwind based schemes. Since the reconstruction across the nodes is difficult even in the upwind direction, we aim for stencils as small as possible.
On every edge of length L, we assume a uniform spacial discretization (x i ) i=1,...,N with grid size x. The global time step is denoted by t and the local CFL number is c = v t x . For the sake of simplicity, we will drop the edge index in the formulation of the scheme and lower indices will refer to the spacial discretization. Thus, the classical update formula of finite volume schemes for the cell averages u n+1 i at time t n+1 reads where F i+ 1 2 are the numerical fluxes at the interfaces.

Scheme of order 4
For a scheme of order P , we need to process for the numerical flux at least d F = P values of u and d = P +2 values for the final update. Thus, for a scheme of order 4, the most compact upwind directed stencil for the flux is given by (u n i , u n+1 i , u n i−1 , u n+1 i−1 ). As indicated in Fig. 1, when solving (1), these values are transported to the vertical time axis at x i+ 1 2 along the characteristics. At the interface, we perform a conservative interpolation to obtain a polynomial p of order 3 with the constraints The corresponding numerical flux is computed by F i+ 1 2 = t n+1 t n vp(τ )dτ and has the form Thus, the final scheme reads This method can be obtained via various interpretations, e.g. a slightly different formulation of the same scheme can be found in [21], named CTCS. It is written in finite difference form with a central stencil, but with the shift (u n i−2 , u n i−1 , u n i ) → (u n i−1 , u n i , u n i+1 ) the scheme (3) can be transformed to the CTCS.

Scheme of order 3
Similarly, we construct a scheme of order 3. As stencil, we just use (u n i , u n+1 i , u n+1 i−1 ). The numerical flux resulting from the second order polynomial interpolation is and the resulting scheme of order 3 Here, we emphasize that the choice of the stencil is not arbitrary. In the following sections, we will discuss some aspects of the stability of the above schemes and possible alternative choices, e.g. Table 1.

Stability
A central basic property of a numerical method is its stability. In the following, we discuss the linear stability of the above schemes. They both can be written in the general form Au n+1 = Bu n .
First, we investigate, as commonly done, the case of periodic boundary conditions. The stability conditions for this case are easily verified, as the resulting matrices are normal. For the case with Dirichlet data, those conditions are no longer sufficient but only necessary for stability (see, e.g. [14,19]).
x 1 x 1 x 1 Unstable combinations are not listed

Theorem 1 (Stability of fourth order scheme) The fourth order method defined in (3) is unconditionally stable.
Proof For the 4th order scheme the matrix A ∈ R N×N is given by The stability of systems with circulant matrices is expressed by the condition where λ A,k and λ B,k are the eigenvalues of A and B respectively. For circulant matrices, they have the form with the frequencies ω = e 2πi N , see [7]. For the six point stencil, we obtain for A Thus, all eigenvalues of the solution operator A −1 B have modulus one and therefore the 4th order scheme is unconditionally stable.
For the scheme of order 3, evaluation of the stability condition is more involved, as the symmetry in A and B is missing. Unfortunately, we are not able to give an analytical proof; therefore, we evaluate (5) numerically to indicate the stability region of the scheme. In Fig. 2, the modulus of the eigenvalues is shown for different CFL numbers and N = 50.
We observe, that for c < 1 the modulus is larger than 1 and thus the scheme is instable. For c > 1 all eigenvalues are smaller or equal than 1, which indicates stability.
For the following, we emphasize that the stability of a scheme can depend on the boundary conditions imposed. Therefore, we now discuss the stability of the two schemes in case of Dirichlet boundary conditions. The layout of the matrices slightly changes, i.e. compared to (4) the upper right block is eliminated and the matrixÃ has the formÃ The new form ofB is analogue. These matrices are no longer normal, which implies that the error amplification factor is now measured by the largest singular value of A −1B , which is not related to its eigenvalues. As we did not find any analytical expressions, we rely on a numerical evaluation of the stability regions. In Fig. 3, the largest singular value is shown for different matrix sizes N and different CFL numbers.
We observe that both schemes are unstable for c < 1, but stable for c > 1. Note that for the scheme of order 4 this differs from the case of periodic boundary conditions. Since for c < 1 the size of the largest singular value increases with N, we do expect this instability to hold for any choice of discretization. Thus, both schemes are unstable for c < 1. But this is a mild drawback, as in this case a classical explicit method can be used.
Furthermore, we applied this analysis on all possible upwind biased combinations of stencils up to order 4. The coefficients of all these methods can be constructed as outlined in Section 3.1. We omit these details for those methods, if not needed for the final hybrid scheme. The results of the individual stability analysis are summarized in Table 1.
Up to order 2 it is possible to find schemes which are globally stable. For order 3 there is only one stable choice available. In most cases the choice of boundary conditions is not affecting the stability, but for the scheme of order 4 and for a, not very intuitive, first order scheme. The unconditionally stable first order scheme is the classical implicit upwind method, that is used as the fallback option of the limiting procedure (see Section 4). Among the two second order schemes is one obtained by the trapezoidal rule. Both are legitimate methods, but we focus on even higher orders to keep the signals on long traveling distances.

Coupling
In order to apply the above methods in a network, suitable numerical coupling procedures are needed at the nodes. Following the derivation using polynomials in time at the interfaces, we can directly extend the schemes using the coupling conditions (2). At a node V , we denote by p e the temporal polynomials of all ingoing edges e ∈ I V . Thus, the temporal evolution of the output is given by The flux over the first cell of the outgoing edges is then With this flux the first cell can be updated. For the construction of the flux at the next interface, we furthermore need to fill the required ghost cells at the left boundary. Since we already have a polynomial representation of the mass in time for the leftmost cell interface, we can reconstruct the cell values as they were derived in Section 3.1. We set With these ghost cells, all cells along the edge can be updated.

A posteriori limiting
High order schemes tend to show spurious oscillations in the presence of discontinuities or sharp gradients. The oscillations can be suppressed by applying an appropriate limiter. There are different possible approaches for the limiting of implicit schemes.
In [1] implicit WENO schemes are proposed, where a spacial WENO reconstruction is combined with an implicit Runge-Kutta scheme. They achieve acceptable results for moderate CFL numbers around C = 3, but the scheme involves large nonlinear systems that become harder to solve, the larger C gets. Another possibility is the use of flux corrected transport techniques [15,20]. They include the transport with a monotone scheme and a anti-diffusive correction by some high order term. Especially for large time steps, this ansatz suffers from the large diffusion of the first order solution and also involves solving large nonlinear systems.
We therefore want to make use of the recent a posteriori limiting technique, benefiting from the special structure of our problem. Originally designed for multidimensional finite element schemes, the MOOD paradigm [5,6] can be adapted to our requirements while eliminating most disadvantages of the options above.
The key idea of a posteriori limiting is to begin with the computation of a high order solution. Only if that solution appears to be unsuitable, the limiting is activated and a lower order solution is computed. Whether a solution is suitable can be checked with some (possibly problem specific) criteria.
The basic procedure is the following: Physical Admissibility Detection (PAD): Most importantly, the solution has to fulfill some basic admissibility criteria in order to ensure that the numerical code does not crash and the variables stay in a meaningful regime.
A classical physical admissibility condition is positivity of the solution. For our application on district heating networks, we incorporate checks for lower and upper bounds of the energy density. The energy must not fall below the return temperature in order to have a well defined consumer model for example. Possible additional conditions are the check for NaN or INF, indicating that something with the solution went wrong.
While the PAD can detect oscillations overshooting the total range, it does not find the ones near steep gradients with moderate values. Those have to be found by numerically analyzing the computed solution. Using a neighborhood of the current cell and evaluating curvature data, we can distinguish local oscillations from the smooth solution as follows.
Extrema Detector (ED): If does not hold, cell u i is a local extremum. This can be either a true extremum of the solution, a local oscillation or some artifact due to round off errors. If we do not find an extremum, the computed value is considered valid. In presence of an extremum in order to identify the right case, we need to look at the local curvatures next to the cell.
We define a discrete approximation to the local curvature (e.g. by central difference). The local curvature indicators are defined as Depending on these indicators, we can distinguish the three cases: with SD = 1 2 . The check can be schematically drawn as in Fig. 4. Whenever the check finds 'invalid' cells, the order of the corresponding flux is reduced. Then, the check is performed again with the new cell data. We use a cascade of 3 schemes with the orders 4 → 3 → 1, where the high order schemes are as described above and the first order scheme is the implicit upwind scheme, which always produces valid results.

Discrete Maximum Principle (DMP):
An important criterion for most numerical schemes is the satisfaction of a discrete maximum principle. This is a very powerful tool for the detection of spurious oscillations, as the high order overshoots would violate the DMP and can therefore be eliminated. For the application of implicit schemes on transport equations, however, we do not have a discrete maximum principle. The long travel distance of cell values makes it impossible to state bounds for the solution. This is why we have to drop the DMP check from our procedure.

Extrapolation and computation of downwind data for MOOD check
The input data required by the MOOD check consists of at least the cell values (u n+1 i−1 , u n+1 i , u n+1 i+1 ) in order to perform the ED step. That means additionally to the computed candidate solution u n+1 i we have to provide the next cell value as well. While inside the domain, we have access to u n i+1 and the computation of the new cell value at i + 1 is possible. Note here, that this new data has to be computed with the same order as the value to be checked. While at the boundary, i.e. at the last cell of any edge in the network, however, we do not have the data required to compute another cell value. In this case, we have to extrapolate u n+1 i+1 only using values inside the domain. Note that also here we require an accuracy of the same order than used before. The numerical flux functions we use are for the 4th order and for the 3rd order.

Convex combination of high order fluxes
The classical MOOD strategy determines the order of the spacial reconstruction polynomial for a given cell. Therefore there is a hard switch of the corresponding fluxes, when the order is decreased. Since we are using the strategy directly on the fluxes, a continuous transition between the different orders is possible. Whenever the order for a given cell has to be reduced from order 4 to order 3 and the 3rd order solution is suitable, we might as well use a flux of the form The parameter α can then be chosen in a way, such that α = max{α|u n+1 i satisfies the MOOD check}.
Unfortunately, this optimization problem is discontinuous, as the MOOD check will only return TRUE or FALSE. The parameter α here is determined by a bisection method. The incorporation of sufficient criteria for the check result into the optimization might accelerate the process. Altogether, we have constructed a hybrid scheme taking advantage of the monotonicity of the first order upwind as well as the good approximation quality of the high order ones.

Application to district heating networks
We designed this scheme for the application of district heating networks. A district heating network can be modeled by a graph describing the connection between the different components. The edges include the pipes, consumers and the source. In the nodes, the coupling conditions between the different components are formulated. In [2] and [17] the modeling of such networks is described in more detail. Here, we want to restrict ourselves to the solution of the energy advection problem. The system consists of a network of advection equations ∂ t e + v∂ x e = 0, as introduced in (1) and (2). The transported quantity is the energy density e and we will ignore source terms modeling thermal losses for now. As described in [2], the velocity for the system is solved separately in every time step by an equation of the form g(v, e, Q) = 0, where Q contains the demand of the consumers of the network. The energies at the outflow nodes of the network are measured (see Section 2) and the new velocities are set to fulfill the demand of all consumers given the measured output. This kind of system can efficiently be solved, see [11,13]. To conclude, a splitting technique is used, where in every time step, we get an actualized velocity depending on the old energy densities and current consumption. This results in a first order approximation of the full system. As the largest numerical errors are due to diffusive effects of the energy transport, the overall error is still much lower than a classical first order scheme for the full system.

Numerical examples
In this section, we analyze the constructed schemes in four different test cases. First, a classical test case for schemes for hyperbolic conservation laws, Shu's linear test [12], is performed and the convergence is investigated numerically. The third test case is a small example of a network, where during the simulation velocities change their sign. An analytic solution of this case is given in [17]. The last test is a part of an existing network geometry with realistic parameters and boundary conditions. We can show that in all cases the new scheme is superior to a classical first order scheme.

Shu's linear test
The scheme is tested on a 1D domain, with transport velocity v = 1 and a demanding initial condition [12]. Due to the iterative structure of the scheme, we cannot use it with periodic boundary conditions (see Section 3.3). Therefore, we extend the spacial domain to from [0, 1] to [0, 2] and initialize the new values with 0. A 0-Dirichlet boundary condition is prescribed at x = 0.
In Fig. 5, the solution is shown for the time t = 1 with hybrid scheme (red) and the unlimited first, third and fourth order schemes (dashed) with n = 200 and a CFL number of c = 2.5. Due to the prescribed boundary condition, all relevant information has traveled to the interval [1,2], which is displayed in the figures.
Compared to the exact solution (black), we can see that the first order solution looses any contour. The unlimited high order schemes, while being more accurate, show oscillations near the steep gradients of the exact solution. These oscillations can be removed by the proposed a posteriori limiting strategy of the hybrid scheme.
In Fig. 6, we show the behavior of our scheme for different CFL numbers c. In all cases, we can observe the cancellation of overshoots in the whole domain. As the 4th order scheme is exact for c = 2, we can expect accurate results for c = 2.5. Even there, we notice that the upwind scheme is smearing out almost all the information from the initial condition. As expected, the numerical diffusion rises with larger time steps, but even for c = 10, we can distinguish the four peaks and the resulting heights produced by the hybrid scheme are a good compromise between fourth and third order approximations.
In Fig. 7, the activation of the limiter is displayed for the test with n = 200 and CF L = 5. The data appears to be coarse, as the corresponding time step is t = 0.025. The colors are according to the local orders o. They are defined as where h and l are the local high and low orders and α is the convex combination parameter from Section 4.2. This leads to a continuous range of o from 1 to 4. We can see that the limiter mostly is active near the sharp gradients of the solution, moving with the features through the domain. In total, only about 26% of the cells have a local order less than four. This means that in most of the cells, only the high order scheme is computed and no comparison to different methods is needed. This rate can be even better, but at the right of the last feature, we observe a  larger domain of reduced orders. Here, the exact solution is constantly zero, but the numerical solution is building a plateau of marginally negative values. This can be reduced by omitting the strict positivity bound, as discussed in Section 6.3.

Convergence analysis
In another test case, we provide smooth initial conditions and check the theoretical order of convergence of the scheme. As for smooth data and small enough discretization the MOOD step should not trigger a recalculation, the scheme fully uses the high order fluxes and therefore achieves 4th order accuracy. We use an initial condition of the form and transport it with v = 1 and CF L = 5 until t end = 0.5. The convergence rate of 3 different variations of our proposed scheme is shown in Fig. 8. With 'hybrid3' we denote a variation, where we only use the 3rd and 1st order fluxes as high and low order respectively, without application of the convex combination. Furthermore, with 'hybrid4' we indicate the use of the full cascade 4 → 3 → 1 but also without convex combination. Finally, 'hybrid' shows the results for the full scheme described above. For a course grid, a mild limiting is still active. As the grid is refined, we can see that for smooth data all the schemes perform at their theoretical order of convergence.

Limitations of the scheme
For explicit methods, the MOOD approach guarantees certain properties by design. However, with the implicit method, there are cases, where the solution with the limiter applied still produces overshoots. This is due to the fact that the flux orders are chosen sequentially throughout the domain. It is possible, that a choice at some point in the domain may locally fulfill the smoothness test but produces problems a few cells further downstream. The limited window of available data for the check makes it impossible to predict these cases. A possible way out would be an iteration process with backtracking elements to eliminate these cases. In practical applications, however, this might imply a huge growth in computation times. For our application on district heating networks, it is not essential to avoid all small overshoots inside the range of admissible temperature values. When the total temperature however gets below a minimal value, the system is not solvable anymore and in these cases, such a mechanism has to be implemented. In realistic settings however this does not happen as near such points, the resulting velocities and pressures would tent to infinity. In real applications, there always is a margin between true solution and the set of solutions that are not admissible.

Triangle network
The first network is forming a simple triangle. It contains three pipes and two consumers. For this network, there is an analytical solution available, which produces a flow reversal and therefore a discontinuity in energy while having continuous boundary-and initial conditions. The details about the solution can be found in [17].
The geometry and initial flow directions are shown in Fig. 9. The source is located at the left and the two consumers are connected to the two nodes at the right side respectively. A time dependent consumer behavior, where the situation continuously changes from the upper consumer demanding more energy than the lower one to a situation vice versa causes a change in flow direction on the vertical pipe to the right. The initial condition is outlined by the brightness on the edges, brighter meaning higher temperature. The results of the temperature evolution of the vertical edge at the lower consumption node are shown in Fig. 10. The general behavior of the exact solution can be explained as follows: starting from the initial condition, the rising infeed temperature enters the pipe. When the change of flow direction is reached, the temperature decreases again, as the values further down that pipe return to the node they have passed before. At the same time, at the upper node of that edge, new, hot water enters the pipe, inducing a discontinuity traveling downwards. At some point this discontinuity reaches the lower end of that edge. Thereafter, we have a linear behavior again.
In Fig. 10, we can see that the hybrid scheme is able to resolve the jump much sharper than the implicit upwind scheme. A significant increase in accuracy of the jump height and the left side value are apparent. For this example, the exact velocity description in time is used for all the schemes. When the advection solver is furthermore coupled to the flow solving algorithm, a precise estimate for the temperatures at the consumers is even more important as it will influence the flow velocities for the next time steps.

Street network
The street network is a part of a real district heating network. It contains 162 pipes and 32 consumers. The total pipe length is 1672 m with the largest pipe measuring 80 m and the shortest having 0.06 m. Even though not being obvious from the topology plot Fig. 11, it does contain one cycle near the label C 2 . In the figure, we can see a temperature front traveling through the network. For the pipes, red color means hotter water, blue is cooler water. The dots represent the consumers with the color coding the mass flow from high (red) to low (blue). The solution shown here was solved by an implicit upwind scheme, to demonstrate the smearing of the sharp front inside the network. The initial condition in the network was T 0 = 80 • C and the boundary condition at the source is set to We want to focus especially on the two consumers marked by red circles, C 1 and C 2 . In Fig. 12, we see the temperature signal in time at the rightmost node (C 1 ) of the network. The resolution for all schemes is x = 0.5m and t = 5s. Since we do not have an analytic description of the solution, the reference is computed by the upwind scheme with 50 times finer resolution ( x = 0.01m and t = 0.1s). We can see that the plain fourth order scheme does not produce stable results due to the extremely high oscillations in front of the shock. The upwind scheme heavily smears out the discontinuity. The hybrid scheme however accurately resolves position and height of the traveling front.  The advantage of the hybrid scheme becomes even more evident when looking at the solution at the consumer furthest away from the source (C 2 ). The more accurate we can compute the temperatures at the consumers, the better the resulting velocities for the next time steps. Large errors in the temperature may therefore lead to completely different flow fields and consequently different looking solutions. This is even more important when there are circles in the network, that can possibly change their flow direction. Such a case can be seen in Fig. 13. Several flow changes and different mixing ratios lead to a rather complicated temperature signal entering the consumer C 2 . We can see that compared to the reference solution, the upwind scheme loses most of the dynamics. The solution of the plain fourth order scheme is omitted, since it too oscillatory to give any meaningful contribution. The hybrid scheme however accurately resolves all the kinks and waves we observe in the reference solution. For a precise simulation of such networks, it seems evident that high order approaches like this are much more powerful than a straight forward first order approach.

Conclusion
We constructed a high order implicit scheme for district heating networks. The used stencil is upwind orientated, as at the coupling points the data cannot be reconstructed properly, due to the non-injective mixing at the nodes. Furthermore, the upwind orientation allows that the solution can be computed successively without solving large systems of equations. The oscillatory behavior of the high order is massively reduced by an a posteriori limiting approach. Even though overshoots cannot always be avoided by this ansatz, the results show a big improvement of the achieved solutions. Compared to a first order approach the constructed scheme achieves a higher precision on coarser grids and larger time steps.
Funding Open Access funding enabled and organized by Projekt DEAL. The authors acknowledge the financial support by the Federal Ministry of Education and Research of Germany in the framework of the project EiFer: Energieeffizienz durch intelligente Fernwärmenetze (grant number 05M18AMB).

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.