Numerical Model for the Analysis of Thermal Transients in District Heating Networks

. As District Heating (DH) networks are experiencing an evolution towards the so-called 4th generation, there is a need to update the currently used models to take into account the ever-increasing complexity of this technology. Indeed, to further improve the reduction in energy consumption and carbon-dioxide emissions, a wide range of technologies and management strategies are being introduced within district heating, such as a large exploitation of Renewable Energy Sources (RES). As a consequence, thermal transients assume a major importance, posing the need to redefine the relevant physical parameters and to develop a model which accurately describes their behaviour. In this framework, this paper proposes a quantitative analysis of the influence of the pipe heat-capacity on the model. Moreover, an equivalent-model, which is able to take into account the two heat capacities of steel and water in just one equation, is proposed and compared with two commonly used approaches. One of the features of the proposed model is the suitability for application to large networks. To prove its capabilities, an application to the Turin district heating network, which is among the largest systems in Europe, is proposed. Results show significant improvements in terms of accuracy over computational time ratio.


Introduction
Nowadays, District Heating (DH) technology is regarded as a core element in the costeffective decarbonization of the European energy system [1].Although DH infrastructures are well-established in Europe since the 1970s [2], in recent years they are facing multiple challenges in order to reach a future renewable non-fossil energy system.Hence, the new developments have led to the conceptualization of 4th generation district heating (4GDH) [3,4].
In order to reach the 4GDH, numerical models are needed to simulate and optimize the configuration of existing and planned district heating networks.In the literature, two main approaches have been adopted for district heating modelling: black-box models and physical models [5].The former are statistical models based on standard transfer function models or neural networks.These methods suffer from low accuracy in the time-delay estimation, especially when temperatures are changed abruptly [5].By contrast, the latter address the physical description of all the relevant components of a network and are preferred in case of large and quick temperature variations or when the estimation of physical parameters within the network is particularly relevant [6].
Among the physical models, many different approaches can be found.The element method and the node method were firstly introduced by Benonysson [7] and later used in many other contributions [5,8].Other popular methods are the characteristic method [9,10], the plug-flow model [11] and the finite volume method [6,10], which is also used in this contribution.
In recent years, together with the evolution towards the fourth-generation systems, revisiting the models has become a necessity.Indeed, these systems will have new features, including lower operating temperatures and lower grid losses, the possibility to recycle heat from low-temperature sources and integrate renewable heat sources and the ability to be part of an integrated smart energy system [3].These new aspects will deeply influence the thermal behavior of the network, leading to the need of a major focus on the prediction of thermal transients.In this framework, the heat capacity of the steel pipe is among the parameters which affect the thermal response of the network.This parameter is explicitly taken into account by various models in literature [8,[12][13][14], while other authors do not consider its influence [10][11]15].Therefore, there is a need to determine its physical relevance.
Hence, the purpose of this study is to provide a quantitative analysis of the influence of the pipe heat-capacity on the thermal response of the network.This goal is reached by a comparison of the two most used approaches, which consist in a) exclusively considering the water heat capacity, using a 1-equation model for the resolution of the energy equation; b) taking into account both water and steel heat capacities, by means of a 2-equation model.Also, a third model is used to take on board the steel heat capacity in a 1-equation model by means of an approximation.
This paper is organized as follows: after this brief introduction, the methodology is presented; then, the results of the different approaches for the pure advection problem in a pipe are discussed and an application to the Turin district heating network, which is among the largest systems in Europe, is presented; finally, the last section is devoted to draw conclusions.

Methodology
The current investigation involved the simulation of the hydraulic and thermal behavior of a district heating network.The simulation is intended to accurately reproduce the evolution of pressures, mass-flow rates and temperatures within the whole system.
To this aim, a pseudo-dynamic thermo-fluid dynamic model was used.This means that the hydraulic problem, expressed by the conservation equation of mass and momentum, is simulated with steady-state conditions.In this case, the negligibility of the unsteady term is allowed by the rapidity of the fluid-dynamic perturbations, which are transferred to the whole network in a period of time of few seconds, much smaller than the time step adopted for calculation.In contrast, temperature perturbations travel at the fluid velocity and could take a long time to be propagated in the whole network.For this reason, the energy conservation equation needs to be solved dynamically.
The problem was treated as one-dimensional and the complex structure of the network was described by means of the graph theory [16].Hence, each pipe was treated as a branch which connects two nodes, corresponding to the inlet and the outlet sections.The network topology was described by means of the incidence matrix A, which has as many rows as the number of nodes () and as many columns as the number of branches ().The general element   is equal to 1 if the i-th node is the inlet node of the j-th branch, -1 if the i-th node is the outlet node of the j-th branch and 0 if the i-th node and the j-th branch are unrelated.
The finite volume method [17] was used to solve the problem.In particular, the continuity and energy equations were integrated over control volumes including each junction node and half of the branches entering or exiting that node.However, in the case of the energy equation, more control volumes were used to reduce the size of the computational mesh and consequently the effect of numerical diffusivity [17].On the other hand, the integration of momentum equation was performed over control volumes including a branch and the two delimiting nodes.
The hydraulic problem can be expressed as follows: Where Eq. ( 1) and Eq. ( 2) are respectively the conservation equations of mass and momentum.While   takes into account the viscous forces,  1 represents the source term and also accounts for the effect of local fluid dynamic resistance due to valves or junctions and the effects of pressure rise due to pumps.The integration of continuity and momentum equation brought to the following formulation of the hydraulic problem: where the unknown terms are the arrays  (length: NB) and  (length: NN), which respectively represent the mass flow rate in each pipe/branch and the pressure in each node.The known terms are: the incidence matrix  (size: NN × NB); the array   (length: NN), which contains the mass flow rates injected in or extracted from the system; the fluid dynamic conductance matrix  (size: NB × NB), accounting for the pressure losses; the vector  (length: NB) that includes the pressure rise due to pumps.Details about the strategy adopted for the solution of the hydraulic problem can be found in Sciacovelli et al. [18].On the other hand, the energy equation was studied in three different forms.All of them assume incompressible fluid and constant specific heat; moreover, the axial thermal conduction and the heat capacity of insulation and ground are neglected in all the cases.The three formulations can be summarized as follows: a) In the first formulation, the heat capacity of the steel pipe was neglected, as done by [10][11]15].This means that just the fluid heat capacity was considered.Hence, the thermal problem was expressed by means of one PDE, reported in Eq. ( 5).
a) With the second formulation, the heat capacity of the pipe was included in the analysis.To do that, the solution of two coupled PDEs -reported in Eq. ( 6) -was needed to solve the thermal problem.This approach to the thermal problem is the standard way the pipe thermal capacity is taken account by various models in literature [8,[12][13][14].Obviously, it also represents the most accurate formulation among the ones proposed in this paper.b) A third formulation was used to include the heat capacity of the pipe in an equivalent one-equation model.This was made possible by assuming that the pipe immediately reaches the thermal equilibrium with the water.Hence, the temperature of the pipe was supposed equal to the temperature of the fluid (the thermal resistance between the two bodies was supposed equal to zero).This can be considered as a reasonable hypothesis since the heat transfer between the water and the pipe is very high.The resulting PDE, reported in Eq. ( 7), differs from Eq. ( 5) in that it considers a greater thermal inertia for the fluid, taking into account both the water and the fluid heat capacities. ( To solve the thermal problem, the various equations were integrated over their corresponding control volumes.The Upwind Differencing Scheme [17] was used to relate boundary and nodal values of temperature.Note that, as previously said, for the solution of the energy equation the number of control volumes was increased to reduce the effect of numerical diffusivity.Therefore,  and  are greater than the ones used for the hydraulic problem.At this stage, the matrix form of the problem -given in Eq. ( 8) -was obtained.
The size and the meaning of the terms in Eq. ( 8) are different according to the model formulation adopted.In the cases (a) and (c) the unknown vector  accounts for the nodal values of water temperature and has a length corresponding to .The mass matrix  is  ×  and considers in the first case just the water thermal inertia, while in the second case it involves both water and steel thermal inertia.The stiffness matrix  and the vector  accounting for the losses are the same in both the cases, and their size is respectively  ×  and  × 1.The size of the problem is doubled when formulation (b) is adopted.Under these circumstances, the unknown vector  includes both water and steel temperatures and has a length equal to 2 .The thermal inertia of water and steel are considered in different rows of the mass matrix  that in this case is 2  × 2 .Also  and  are respectively 2  × 2  and 2  × 1.
In order to test and compare the different formulations, some simulations concerning an adiabatic pure advection problem were performed over a pipe.The pipe under analysis was 1 km long and had internal diameter, external diameter and insulation thickness respectively equal to 50 mm, 125 mm and 29 mm.The initial temperature of the water in the pipe was supposed to be 20 °C.Then, a mass flow rate at 120 °C was injected at the inlet section.Its velocity was imposed equal to 0.3 m/s.A transient of 80 minutes was simulated by means of the three thermal models described above.In this case, the simulations were conducted with very fine spatial and temporal discretization (Δ = 0.01  and Δ = 0.1 ), in order to minimize the numerical error.
Afterwards, the two fastest approaches, i.e. models (a) and (c), were compared in the case of a real distribution network, belonging to the Turin district heating system, which is among the largest in Europe.The selected distribution network is represented in Figure 1.It connects 11 buildings to the main pipeline.The thermo-fluid dynamic simulation was performed both on the supply and return lines.This allowed evaluating the heat flux request of the distribution network at the barycenter (BCT -i.e. the connection point between the distribution network and the transport network), expressed by Eq. ( 9): Φ , () =  , ()  (  −  , ()) where  , represents the total mass flow rate that flows in the network,   is the supply temperature -equal to 120 °C in the proposed application -and  , is the temperature of the water exiting the distribution network, evaluated by means of the model.The simulation was performed for a typical winter day with the data made available by IREN.

Results and discussion
The aim of this paper was to provide a quantitative analysis of the influence of the steel pipe heat capacity on the thermal transients of district heating networks.Indeed, although this parameter is explicitly considered by various models in literature [8,[12][13][14], other authors do not include it in their models [10][11]15].Therefore, there is a need to identify its physical relevance, in line with the growing interest towards an accurate prediction of thermal transients, which will become crucial for 4 th generation district heating.
To this purpose, some numerical tests were conducted using different approaches.In a first case, the steel heat capacity was ignored.Secondly, it was taken into account by introducing a second equation in the model.Finally, an equivalent model which considers both the heat capacities of water and steel in a single equation, by means of an approximation regarding the thermal resistance between the two, was proposed.The three approaches, which are fully detailed in the methodology section, have been used to carry out a numerical simulation of a pure advection transient problem over a single pipe (length = 1 km, internal diameter = 50 mm, external diameter = 125 mm, insulation thickness = 29 mm, fluid velocity = 0.3 m/s, initial temperature = 20 °C, inlet temperature = 120 °C).The results of this analysis are shown in Figure 2, Figure 3 and Table 1.Explanation of these results is given in the following.
Figure 2 shows the temperature distribution over the pipe after 40 minutes.By comparing the results obtained using methods (a) and (b), it turns out that the consideration of the heat capacity noticeably changes the solution of the problem, since it is responsible for the cooldown of a significant portion of the mass-flow rate, whose hydraulic front is located at x = 720 m.Therefore, the graph seems to suggest that is worth to not neglect the pipe heat capacity contribution.On the other hand, the solution obtained with approach (c) differs from solution (b), which is the most accurate one, just by the fact that it is not able to reproduce the smoothing effect.However, generally speaking, it could be observed that this error could counteract the truncation error introduced by artificial diffusivity when the discretization is not fine enough.
In Figure 3, the temperature evolution of the control volume located at x = 0.5 km is illustrated.It is possible to notice that the difference between (a) and (b) consists in a slowdown of the thermal perturbation.While in the first case the temperature of the analyzed control volume rises to 120 °C after 27.8 minutes, in the second case the temperature gradually increases from 15 °C to 120 °C in the time range between 35 and 60 minutes.Approach (c) is able to take into account the time delay due to the greater thermal inertia and shifts the temperature change at t = 47.2 min (with a 19.4 minutes delay with respect to the model without pipe thermal inertia).Again, it can be observed that the graduality given by approach (b) is lost, since the temperature change is sharp.However, as previously mentioned, this error can be considered as acceptable.
Obviously, the magnitude of the delay of the temperature perturbation depends on the volume of the steel of the pipe considered, and, more in detail, on its ratio with respect to the water volume.This delay (for a given control volume CVP located at xP) can be derived as: (10) where    ,() and    ,() are the time it takes for the control volume CVP to reach the target temperature respectively with approach (a) and approach (c),   is the hydraulic velocity of water and  ℎ, is an equivalent thermal velocity, which can be defined as follows: ℎ, =           +        (11) As a consequence, the time delay linearly depends on the volumes ratio, and quadratically depends on the diameters ratio, as shown in Figure 4.In the specific case, the equivalent thermal velocity assumes the value of about 0.18 m/s, with a hydraulic velocity of 0.3 m/s.The computational time of the three approaches is compared in Table 1.Inspection of the table clearly indicates that the computational cost of method (b) is drastically higher than the other two approaches.For this reason, despite it is the most accurate one, it is inapplicable for large and real networks.On the other hand, the times required by approaches (a) and (c) are of the same order of magnitude (which is more than 97% less than the one required for approach (b)).
Overall, this analysis suggests that the preferred model would be method (c) since it adequately approximates the effect of the thermal capacity of the pipe without overloading the simulation.To sum up, it can be considered as the best trade-off between performance and accuracy.Fig. 4. Time delay of the temperature perturbation at x = 0.5 km produced by approach (c) with respect to the solution given by approach (a), as a function of the ratio between the external and the internal diameter of the pipe.The circle in red indicates the case of the pipe studied in this paper.Table 1.Computational time of the three methods for the solution of a 80-minutes-lasting transient pure advection problem in a pipe 1 km long, with dx = 0.01 m and dt = 0.1 s.All simulations were run on a PC laptop with a total of 16 GB of memory and an Intel i7-8565U CPU @ 1.80 GHz.The software MATLAB® was used and the function adopted for the solution of the linear systems was mldivide with its default options.

Pipe heat capacity
Equivalent approach

Computational time [s]
(a) 1 Finally, a winter application to a distribution network belonging to the Turin district heating network is analyzed.Results are displayed in Figure 5.In detail, Figure 5(a) represents the temperature of the mass-flow rate coming back to the barycenter, while Figure 5(b) illustrates the thermal load of the whole distribution network.Both the figures report the daily evolution obtained using approaches (a) and (c).It could be inferred that the consideration of the thermal capacity of the pipe modifies both the temporal profiles.In particular, the thermal delay due to the increased heat capacity is particularly evident in the first hours of the day, during the heating phase of the system.Despite the two curves do not present dramatic differences (the relative error is in all the cases less than 5%), it is worth to observe that the modification of the thermal response of the network may cause important implications in some optimization applications.An example could be the determination of the optimal set of anticipation to apply in a demand-side management strategy for district heating systems.In this kind of application, the results strongly depend on the temperature dynamics within the network.Therefore, differences in the temperature profile like the ones shown in Figure 5(a) could significantly influence the determination of the optimal solution, suggesting a worthy subject for future works.

Conclusions
In this paper, a numerical model for the analysis of thermal transients in district heating network is proposed.To reach this goal, a quantitative analysis of the influence of the heat capacity of the pipe on the thermal behavior of the network is carried out.The scientific interest towards this theme arises from the growing importance that the thermal transients are assuming in the simulation of district heating networks, related to their evolution towards the 4 th generation.Moreover, a literature review shows that while some models do take into account this parameter, other authors neglect its influence.Hence, it is crucial to determine whether is fundamental to consider its influence in a proper thermo-fluid dynamic model.
Three different approaches are proposed and compared in this work.The first one neglects the heat capacity of the pipe, as done by [10][11]15].In the second, both the heat capacities of water and pipe are considered.This is obtained by adding a further equation in the model.Therefore, in this case the thermal problem is made up by two coupled partial differential equations.This second model is the most accurate one and represents the standard way the pipe heat capacity is taken into account by many models in literature [8,[12][13][14].Finally, a third model takes into account water and steel heat capacities in just one equation, by approximating the water and steel temperatures as equal.
From a straightforward analysis of a pure advection problem in a pipe, it appears evident that the heat capacity of the pipe has a significant influence on the thermal response of the system, in that it is responsible for a time delay in temperature propagation, which is particularly evident in the start-up heating phase considered.Hence, it is worth to include this parameter in the simulation models.On the other hand, it is also shown that the approximation obtained by including the steel heat capacity in a one-equation model is affordable, since the temperature delay is reproduced in a sufficiently accurate way.Moreover, the use of a one-equation model noticeably reduced the computational effort required for the simulation.The computational time is 97.5 % less than the one needed for the two-equations model, which would be unusable for real and large networks.
Finally, an application of the equivalent one-equation model to the Turin district heating network illustrates the changes in temperature dynamics which can be caught thanks to the consideration of the steel pipe thermal capacity in the model.
Overall, the results obtained from the analysis bring to the conclusion that the equivalent one-equation model is the best trade-off between accuracy and computational cost.Hence, it is suggested as reference for the development of future numerical models.

Fig. 1 .
Fig. 1.Schematic representation of the topology of the distribution network under analysis.

Fig. 2 .
Fig. 2. Temperature distributions of water -obtained with model (a), (b) and (c) -and steel -obtained with model (b) -for a pure advection transient problem in a pipe (length = 1 km, internal diameter = 50 mm, external diameter = 125 mm, insulation thickness = 29 mm, fluid velocity = 0.3 m/s, initial temperature = 20 °C, inlet temperature = 120 °C) after 40 minutes.Note that the hydraulic front is located at x = 720 m.

Fig. 5 .
Fig. 5. Temperature of the mass-flow rate coming back to the barycenter (a) and thermal load of the distribution network (b) depicted in Figure 1, obtained with the simulation of a 24-hours thermal transient for a typical winter day.The black line represents the solution obtained without considering the steel pipe heat capacity; the red line represents the solution obtained with the equivalent model which considers both the water and the steel pipe thermal capacities. +