DYNAMIC TEMPERATURE SIMULATION IN DISTRICT HEATING SYSTEMS IN DENMARK REGARDING PRONOUNCED TRANSIENT BEHAVIOUR

. A dynamic performance of district heating systems was analysed with an emphasis on temperature profile distortion throughout a heating system network. Therefore, a modelling approach (the so-called node method) developed at the Technical University of Denmark was applied. For comparison purposes, commercial TERMIS software was also used in this work. Typical supply conditions were investigated in the considered district heating systems in Denmark. Large and sudden temperature changes in supply temperature were depicted in the district heating system in Madumvej while the pronounced transient supply conditions of temperature wave were exhibited in the district heating system in Naestved. Time-dependent consumer data from the district heating systems were applied to compare the results obtained from modelling approaches. Based on the analysis of temperature wave propagation through the network, it was noted that temperature wave spread unevenly on different locations of the network.


Introduction
Modern district heating systems and carbon lean technologies involve the multiple usage of energy due to the possibility of redistributing heat derived from various sources such as combined heat and power plants, geothermal resources, industrial processes and other resources. Such flexibility in heat source selection facilitates the optimal usage of energy resources by varying the amount of supplied heat depending on demand, price and availability. This introduces the first source of transient regimes in district heating systems. Another reason for the occurrence of transient regimes is heat load variation in consumers due to changes in weather conditions and daily demand. As a result, operation conditions are never stable and involve pronounced thermal and hydraulic transient regimes, particularly temperature waves combined with fluctuations in temperature. These conditions were pronounced in the considered district heating systems in Denmark. Large and sudden temperature changes in supply temperature were depicted in the district heating system in Madumvej while the pronounced transient supply conditions of temperature wave were exhibited in the district heating system in Naestved.
Due to variation in the amount of supplied heat, the use of complex operational strategies becomes a necessary part for efficient system performance, which can be implemented by modelling a dynamic performance of a district heating system. The prediction of operational regimes, including thermal and hydraulic regimes in district heating systems, is required for short and long term forecasting and for performing dynamic temperature control of district heating systems.
Many papers on district heating systems and their performance have been published in recent years. Some of those have focused on district heating modelling and presented simulation results. Larsen with co-authors discussed strategies for reducing computational time and proposed new methods for reducing the size of the network (Larsen et al. , 2004. Another group of codes simulates the performance of a district heating plant in different operation modes: generating heat or power and generating heat and power simultaneously (Fonseca and Schneider 2006). The simulation results of a new alternative control approach for indirectly connected district heating substations were presented in (Gustafsson et al. 2010) where the authors performed energetic analysis and compared the results obtained using different codes. Other simulations concern the comparison between district heating and micro-combined heat and power plants (CHP) where a different degree of decentralization is taken into account (Nystedt et al. 2006). Some papers have examined particular technological aspects, e.g. the performance of geothermal district heating was illustrated in (Hepbasli and Canakci 2003) with an example of a geothermal district heating system built in Turkey. The energy and exergy losses of a district heating network were evaluated in (Westin and Lagergren 2002). Merging of energy and environmental analyses were proposed in Torchio et al. (2009). Some papers depict a steady-state model for estimating district heating hydraulics (Gabrielaitiene et al. 2002;Stevanovic et al. 2007) while other models describe simulation results concerned with the transient behaviour of DH systems. More particularly, modelling transient behaviour in several branches of a pipe network has been presented in (Gabrielaitiene et al. 2003;Pálsson et al. 1999). Modelling a district heating system has been presented in (Larsen et al. 2004) where information about time-dependent consumer behaviour is limited. In the above-mentioned work, data on timedependent heat consumption were only available for some of the consumers. For the remaining part of users, the averaged heat consumption was used. The unique aspect of this work is the application of time-dependent measured data available for all consumers in district heating systems, which enabled to represent the realistic performance of the simulation model.
Two models have been employed in this study, namely the node method developed at the Technical University of Denmark (Pálsson 2000) and commercial TERMIS software. In these models, dynamic temperature simulations are based on idealized flow conditions where changes in temperature profile are influenced by flow velocity and heat accumulation in the insulated district heating pipelines. Other methods such as the element method (Gabrielaitiene et al. 2008) and the model developed in (Stevanovic et al. 2009) have not been applied because of their similarities to the chosen methods or due to their modelling limitations.
In the present work, we focus on the analysis of system behaviour and the ways it could be represented by modelling tools for the pronounced thermal and hydraulic transient regimes, particularly temperature waves combined with temperature fluctuations. For this purpose, two district heating systems and their operational conditions were analyzed. The measured data acquired from these systems were applied for a comparison of results obtained from modelling approaches.

Case 1: Naestved District Heating System
The Naestved district heating subsystem serving l0% (i.e., 6 MW) of the total system heat loads was analysed. Time-dependent heat consumption data were measured for all consumers in the subsystem, which formed the basis of comprehensive analysis. Fig. 1 shows the distribution network of the Naestved subsystem.
The network is constructed of pre-insulated pipes with the total length of approximately 5 km. The pipe diameters range from 40 mm to 250 mm. The amount of the heat supplied to the subsystem was measured and the measurement point was treated in modelling as a heat source and therefore called so hereafter (Fig. 1). In the study, a consumer (shown as a circle in Fig. 1) was regarded as a substation, which generally served a single building. However, in some cases, the substation served a group of buildings, creating a secondary network not considered in this work.
The heat supplied to the subsystem was measured along with data on heat consumption, including power, supply and return temperature measurements on the primary side of consumers' substations. The outdoor temperature is shown in Fig. 2 along with the amount of the heat supplied to the subsystem. The Naestved district heating subsystem supplies heat to 41 consumers and includes apartment blocks and administrative and industry premises. These data were analyzed thoroughly in (Gabrielaitiene et al. 2007).
Most of the measurements were available on an hourly basis, except for five-minute readings for the six largest consumers, which represented 50% of the total heat loads. The temperature supplied to the subsystem was also measured on a five-minute basis. As can be seen from Fig. 3, the supplied temperature formed temperature wave deliberately produced to reduce flow rate during peak consumption in the morning period. The measurements of two working days in January were used.  The Madumvej district heating subsystem serves 14 consumers and provides approximately 3 MW of the total heat loads (Fig. 4). The subsystem is connected to a large district heating system through a heat exchanger, which in the model was treated as a heat source. The Madumvej distribution network is shown in Fig. 4. It is built up of pre-insulated pipes with a total length of 2.71 km. The pipe diameters range from 40 mm to 200 mm. Consumers are identified with numbers and are shown in Fig. 4. Consumers are connected to the network through heat exchangers and are equipped with a local control system causing variations in flow rate. More information about Madumvej district heating can be found in (Gabrielaitiene et al. 2006).
The data used in the presented study contained measurements at a primary side of the consumer's installations. They include the measurements of power and supply/return temperatures collected at 6-minute intervals.
The measurements of the heat and temperatures supplied to the subsystem were collected by the Supervision, Control and Data Acquisition system (SCADA). These measurements are presented in Figs 5 and 6 for the first (25-26 March) and second considered time period (29-30 March) respectively.

Modelling District Heating Systems
In this work, a node method implemented in DHSIM simulation program developed at the Technical University of Denmark and Risoe National Laboratory has been used. This program was previously used in investigations by (Pálsson et al. 1999;Bøhm et al. 2002) and is based on the so-called node model for estimating temperature dynamics in pipelines (Benonysson 1991). Commercial TERMIS software (TERMIS 2010) has also been used in the study for comparison purposes. Both models are based on the quasi-dynamic approach where temperature is estimated dynamically while flow and pressure are calculated on the basis of a static flow model. Fluid pipe flow is assumed to be one-dimensional and considered as an idealized flow. The temperature profile under such conditions is influenced by flow velocity and heat accumulation in the pipe while the impact of turbulent fluctua-tion and secondary flow appearance are neglected. Eq. 1 presents the basis for estimating temperature dynamics: where T f and T p -temperature for fluid and surroundings respectively; V -volume flow m 3 /s; A f -a cross sectional area of the pipe (m 2 ); A -an area of a pipe internal surface per unit (1m) length (m); ρ f -density (kg/m 3 ); h -the heat transfer coefficient (W/m 2 K); c p -specific heat (J/kgK).
The difference between the models lies in the estimation of heat transfer between the fluid and its surrounding. More particularly, the heat capacity of the metal pipe and heat accumulation in pipes is neglected in TERMIS software. In the node method, the overall heat transfer coefficient is based on the method of thermal resistances (Holman 2002) and the influence of two neighbouring pipes (the return and supply pipes) is accounted for. As this approach is not used in TERMIS, the heat transfer coefficient estimated by the node method was applied in TERMIS software.
To ensure similar modelling conditions, the same time step length was used in both models, which together with pipelines discretization, determined the accuracy of simulations. The discretization of pipelines was handled automatically and depended on time step length.
The dynamic behaviour of a consumer was represented by the measured time-dependent heat load and by time-dependent temperature difference between supply and return primary temperatures. The latter described the cooling ability of customer substation. The measured supply temperature at a heat source was used as model input value.

Results of Supply Temperature Simulations of Consumers for the Naestved System
In this section, the results of supply temperature simulation of consumers are described for the Naestved district heating system. Predictions about supply temperature are presented for consumers located at different distances from the heat production unit. The consumers closest to the heat source (EE041 and EE045 in Fig. 1) were located at 0.1 km and 0.2l km respectively. The following consumers were located at a greater distance from the heat source: i) consumers located at a fairly long straight pipeline; ii) consumers located at a distant pipeline with many bends and fittings. In Fig. 1, the location of two consumers in group i) are identified as EE036 and N1263 (distances 1.391 km and 1.173 km respectively) while two consumers in group ii) are labelled as EE020 and EE019 (distances 1.126 km and 1.136 km respectively). The modelling results of those consumers are shown in Figs 7-12.
The results from the node method are only shown in the figures because difference in the results obtained using two models was slight and hardly visible in the figures. To complement the figures, this obtained information was presented in the digital form in Table l. It shows the difference between the measured and predicted temperatures averaged over the period when morning temperature wave occurs (and called average error). The negative value of average error indicates that the temperature is under-predicted and vice versa. A standard deviation of this error was also shown in Table l. In Fig. 12, the results obtained using two models were presented, because they differ significantly from each other, especially within the period up to 3 hours. Figs 7-11 indicate that the trend towards predicted and measured temperature profiles is similar, except for consumers. N1263 and EE036 are located at the distant bendy pipeline (Figs 11 and 12). For those consumers the difference between prediction and actual measurements is somewhat larger in terms of time for temperature wave arrival and the time period during which temperature wave occurs at a consumer (i.e. wave spreading time).
In other words, spreading time can be defined as the time period during which temperature wave starts and ends at a given spatial location. It appears that temperature wave underwent the greatest distortion when reaching those consumers. The spreading time of temperature wave was 5 hours, which is one hour more than for consumers EE020 and EE019 (Figs 9 and 10), although they  are located at a longer distance from the heat source. Identical spreading time (4 hours) as observed for consumers 8E045 and E8041 is located closest to the heat source (Figs 7 and 8). A possible reason for large temperature wave spreading for consumers N1263 and EE036 (Figs 11 and 12) is that the wave travelled through a large number of turns, bends and other fittings in preceding pipelines. Heat transfer in such fittings is more intensive due to velocity profile distortions and secondary flow appearance. Thermal diffusion rates are then more pronounced and as a result, temperature wave spreads and propagates more quickly. These effects are not accounted for in the present models, and therefore the discrepancy between the predicted and measured temperature is greater.
Other factors that can affect temperature wave spreading such as the velocity and magnitude of its fluctuation have also been examined. It was found that the average velocity in preceding pipelines for consumers EE020, EE019 and N1263, EE036 was similar and differed on only l6% during morning wave occurrence, which is small compared to velocity differences in other parts of the network. Velocity fluctuation at consumers' substations could not affect temperature wave either, because they are not significant (less than 2 times) during the hours when temperature wave occurs for all consumers. An exception is observed for consumer EE036 where velocity increases up to 12 times during 3-5 hours (Fig. 12). It coincides with temperature wave occurrence, which is expected to increase wave spreading. Nevertheless, spreading for consumer EE036 is smaller than for consumers EE020 and EE019. Thus, spreading temperature wave is expected to be greater in pipelines with a large number of turns, bends and other fittings. The principle used for estimating minor pressure losses can also be used to account for this phenomenon. As additional components (valves, bends, tees etc.) added to the overall heat loss, they also add to the overall temperature wave spreading. An increase in temperature spreading, in principle, can be described in a dimensionless form and based on experimental data in a similar manner to express pressure drop across the component in terms of loss coefficient. This suggestion is presented in general terms and additional investigation is required for describing it fur-ther in detail. From an additional analysis of temperature wave prediction at other consumers, it was found that difference in predictions obtained by the two simulation models was insignificant. This statement corresponds to the analysed data where only relatively small and slow supply temperature increases occurred in the heat source. The supply temperature wave in the analyzed data reached a maximum after a long period (i.e., 2 hours) and the amplitude of temperature wave was fairly small and made namely 6 °C for the first and 7 °C for the second day. It is expected that in case of a large increase in supply temperature, the difference between the results obtained using the two models would be more pronounced. It is based on the fact that TERMIS neglects pipe heat capacity, thus heat accumulation in the pipe leads to the early arrival of temperature front, especially for distant pipelines. This statement is supported by a comparison of modelling results obtained with and without pipe heat capacity in Benonysson (1991). However, a more comprehensive analysis of this case was beyond the scope of this work due to a lack of experimental data.

Results of Supply Temperature Simulations of Consumers for the Madumvej system
In this section, the propagation of steep changes in supply temperature and their prediction are detailed for the Madumvej district heating system. Two parameters are introduced to describe a comparison between the predicted and measured temperature values in a systematic way. The first parameter, the average error, is estimated by finding the difference between the predicted and measured values at each time step, then averaging them over the simulation period. This parameter measured the quality of heat loss estimation. The second parameter describes a standard deviation of this error, which is a measure of how accurate time delay between a heat source and consumers is evaluated.
The above-mentioned parameters are presented in Tables 2 and 3 to predict supply temperature for all consumers. A negative value of the average error indicates that temperature is under-predicted and vice versa.
Tables 2 and 3 disclose that the average error is considerably greater for industrial consumer 114. The reason for this is low velocity regime in the consumers' intake pipe, which causes high heat losses compared to other pipes. The average velocity in this pipe was in the order of 0.02-0.04 m/s combined with velocity drop to as low as 0.001-0.01 m/s (see Figs 13 and 14) that occurred due to the small heat load for the simulated days. Flow regime in the consumers' intake pipe is very important under such extreme conditions because others, for example consumers 112 and 103 that are connected to the same branch together with consumer 114 (branch N06-N08-N09 in Fig. 4), have a smaller average error, i.e. less than 1 °C. Thus, it is more likely that low average velocities, in the order of 0.04 m/s and less, combined with velocity drop to the values very close to 0 will result in large deviations from the measurements. The relative error is greatest for consumer 114 during the period of time from 25 th to 26 th of March for both modelling approaches (see Table 2). However, for other period of time (29)(30), the difference between the approaches is not uniform and TERMIS produced somewhat larger relative errors (see Table 3).
In Figs 15-18, the propagation of steep changes in supply temperature and their prediction are detailed for the second considered time period. We consider only the cases with moderate velocity fluctuations, thus excluding the effect large fluctuations might have on temperature changes. These conditions were found in a part of the network supplying heat to four consumers, including 108, 109, 107 and 113 (Fig. 4), the velocity fluctuation of which can be observed in Figs 15-18 along with average velocity and Reynolds numbers in the consumers' intake pipe. This information is also presented for pipelines prior to the intake pipes of these consumers (Fig. 19). Figs 5-18 also indicate that the local fluctuations in velocity are moderate-maximum of 10-30% from the average value, and velocity level is fairly constant over the presented time period (from 29 th March 16:00 to 30 th March 04:00). Such fluctuations in velocity can in principal be approximated as periodical, which at relative amplitudes of 10-30%, produce a small variation in local heat flux quantities. However, this has no measurable effect on the average heat transfer (Benim et al. 2004). Consequently, it would be appropriate to conclude that fluctuations in such velocity would not affect fluid temperature significantly.
It was found that sudden and big changes in temperature at the heat source resulted in an inadequate prediction of temperature values at consumers' substations, whilst predictions about the time moments of these changes highly correlate well with the performed measurements. Two impulse-like temperature changes, considerably larger than other changes, were observed at the heat source after 17 o'clock and 22 o'clock (Fig. 16). The first impulse having a value of 16 °C between the lowest and highest impulse points introduces a deviation between the predicted and measured temperature higher than 4 °C.
The second impulse (7.2 °C) produces a deviation between the predicted and measured temperature higher than 2 °C while the typical value of maximum scattering between those values is found to be within 2 °C when changes in temperature were smaller than impulses. The prediction of the above-mentioned impulses is shown in Figs 15-19 for consumers 107, 108, 109 and 113.
Moreover, unfavourable conditions from the modelling viewpoint constitute the situations where large initial changes in temperature are present. It appears that under such conditions, the prediction of temperature values is less reliable. Fig. 18. Temperature and velocity for consumer 113 for the second considered time period Fig. 19. Measured velocity in the pipelines preceding consumers 107, 108, 109 and 113 for the second considered time period. The location of the pipelines is shown in Fig. 1

Concluding Discussion
The dynamic behaviour of two district heating systems, Naestved and Madumvej (Denmark), has been analyzed. It has been modelled by two approaches: the node method developed at the Technical University of Denmark and applying commercial TERMIS software. The measured data from the district heating systems were applied to compare the results obtained from modelling approaches.
The results received from modelling the Naestved district heating system have disclosed that difference in predictions obtained using two models was insignificant for a relatively small and slow increase in temperature observed in the measured data (changes in temperature made 6-7 °C during 2 hours). The differences between two models are pronounced when significant variations in flow rate are observed.
A comparison between measured and simulated data concludes that the prediction of temperature values deviates significantly from the measured values in the following cases: i) during relatively big and sudden changes in temperature at the heat source (e.g. an impulse of 8-16 °C produced for approximately 20 minutes) and ii) during the periods of low velocities (in the order of 0.04 m/s and less). The analysis of temperature wave propagation through the network has revealed that temperature wave spreads unevenly on different locations of the network. The achieved results indicate that discrepancies between the predicted and measured temperatures are pronounced for the consumers located at distant pipelines.