- Email: [email protected]
Prediction of thermal transients in district heating systems
Energy Conversion and Management 50 (2009) 2167–2173
Contents lists available at ScienceDirect
Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman
Prediction of thermal transients in district heating systems Vladimir D. Stevanovic a,*, Branislav Zivkovic a, Sanja Prica a, Blazenka Maslovaric a, Vladan Karamarkovic b, Vojin Trkulja c a
Faculty of Mechanical Engineering, University of Belgrade, Kraljice Marije 16, 11120 Belgrade 35, Serbia Faculty of Mechanical Engineering, Dositejeva 19, 36000 Kraljevo, Serbia c JKP Beogradske Elektrane, Savski Nasip 11, 11070 Belgrade, Serbia b
a r t i c l e
i n f o
Article history: Received 3 January 2008 Received in revised form 21 September 2008 Accepted 27 April 2009 Available online 27 May 2009 Keywords: District heating system Thermal transients Simulation
a b s t r a c t A prediction of thermal transients in district heating systems is important in order to adjust in an energy efficient manner the loads of heat power plants and pump stations with dynamic consumption needs during changes of the outside air temperature, the wind intensity, the solar radiation, or system startups and shut-downs. A model and corresponding computer code are developed with the aim to simulate the thermal transients in district heating systems. They are based on the high-order accurate numerical solution of the transient energy equation, and the hydraulic prediction of pressure and fluid flow rates within the complex pipe network. Thermal transients caused by an increase and decrease of the heat power plant load are simulated for real operating conditions of the district heating system. Predicted temperature front propagations show a good agreement with data measured at three consumer substations located in different parts of the district heating network and at different distances from the heat plant. The developed computational tool provides reliable information about time periods of temperature fronts propagation and heat distribution from the heat source to consumers within the whole district heating network. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Managing a district heating system operation during temperature transients is important for the reliable supply of adequate heat amounts to consumers in an energy efficient manner. The temperature transients occur during changes of weather conditions, such as a change of temperature or wind intensity, during start-ups and shut-downs of district heating system operation, or due to the change of hot water consumption in the systems that provide sanitary water heating. The efficiency of heat supply during transients depends on the timing of heat source power according to the heat consumption, on the time periods necessary for the heat transport from the heat source plant to consumers’ substations, as well as on the heat accumulation within the distribution pipelines and heat losses to the surroundings. Due to the complexity of district heating networks, and dynamic characteristics of heat consumption and transport properties, the prediction of transient regimes is possible only with the application of numerical models and computer codes. One of the early dynamic models for the transient temperature predictions in the complex district heating systems was presented almost two decades ago [1]. The model was based on the lumped * Corresponding author. Tel.: +381 11 3370 561; fax: +381 11 3370 364. E-mail address: [email protected] (V.D. Stevanovic). 0196-8904/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.enconman.2009.04.034
parameter energy model of the district heating pipeline network. In this approach the whole pipe is described with the single control volume, the linear water temperature change along the pipe is assumed and the mean water temperature is assigned to the pipe as the simple average of the input and output water pipe temperatures. Such an approach is not able to predict the temperature wave propagation along the pipe. Hence, a large diffusion of the temperature waves propagation within the pipeline network is expected, together with a low accuracy prediction of time delays of temperature waves propagations from a heat source to consumers. Due to the complexity of pipeline networks of district heating systems, the aggregated dynamic models have been developed [2,3], which represent a number of pipeline segments with only one equivalent pipe. Aggregated models are simpler and they provide faster computer simulations, but they are not able at all to predict temperature waves propagations from heat sources to consumers within the pipeline network. A further development towards the real physical representation of temperature transients was achieved with the node method [4,5]. The principle of the node method is that it keeps trace of how long time a water mass, which in the actual time step arrives at a node, has been on its way from the last node. Based on time series for the temperature history of the different nodes the temperature of the mass at the inlet of the pipe is calculated, and on the basis of this as well as the heat loss and heat capacity of the pipe, the temperature of the mass
2168
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
Nomenclature A e g h q_ p t
cross sectional area, m2 fluid energy, J/kg gravity, m/s2 specific enthalpy, J/kg volumetric heat source, W/m3 pressure, Pa time, s, min
at the outlet of the pipe is evaluated. In this manner a new temperature for all nodes in the network can be calculated in every time step. As in the case of the lumped parameter model [1], the node method does not solve the temperature wave propagation along the pipe; hence, the temperature wave diffusion and smearing could be expected. The real physical model of the temperature waves propagation within the pipeline network of a district heating system will assume a direct numerical solving of the energy balance equation of water flow along all pipes. It is reported that such an approach is applied recently in the commercial computer code TERMIS [6,7], but no information about the numerical solving procedure is available in the open literature. The developed node method and the TERMIS code are tested for the conditions of temperature fronts propagations in parts of district heating systems [6]. At present, numerical models are applied at several district heating systems in order to optimize transient operational regimes; first results show the saving of the total fuel energy consumption of 2% at least [7,8]. Dynamics of pressure and flow changes within district heating networks are determined by the pressure waves propagation. The speed of pressure waves in the hot water of district heating systems is approximately 1000 m/s, which is 1000 times faster than the water velocity that is in the range from below 1 m/s till 2 or 3 m/s (these higher velocities are encountered in main transport district heating pipelines). The temperature variations propagate with the velocity of the water flow. Hence, the dynamics of the flow in the network are of minor importance compared with the dynamics of the temperature changes, from an operational optimization point of view. On this basis, the modeling approaches [1–6] for the prediction of temperature transients applied the steady-state model of water flow, where the pressure and flow distributions are calculated for every time instant of the numerical integration, based on the actual distribution of temperature and corresponding thermo-physical characteristics of the water in the pipeline network. This paper presents the new numerical model for the simulation and analyses of the thermal transients in the district heating system. It is based on the numerical solution of the transient energy equation with the numerical scheme of the third-order accuracy in space, and on the steady-state hydraulic solution of pressure changes and flow rates within the network. The developed method is applied for the simulation of two thermal transients in the district heating system of Zemun in Serbia, caused by an increase and decrease of the heat power plant load. The obtained numerical results are compared with data measured at three consumers substations located at different distances and in different parts of the pipeline network. The obtained good agreement shows that the developed numerical method is a reliable tool for the simulation and analyses of district heating system thermal transients. The presented modeling approach is an advantage in comparison to the lumped parameter, aggregated and node methods [1–4], since it provides the prediction of distributed temperatures along all pipes within the district heating system, and enables the accurate calculation of dynamic temperature changes without the influence of the numerical diffusion in predictions of temperature waves propagations. A com-
u x
velocity, m/s coordinate, m
Greek symbols h pipe inclination, rad q density, kg/m3 sw wall friction per channel unit volume, N/m3
parison with the recently introduced temperature dynamic model of the TERMIS code could not be done, since the TERMIS numerical model is not reported in the open literature. Also, here presented experimental data on temperature waves propagations in the complex district heating system is a contribution to the limited information about this process. 2. Model development 2.1. Transient energy balance Transient one-dimensional energy balance equation is written as
@ðqeÞ @ðqueÞ @p þ ¼ qgu sin h þ q_ @t @t @t
ð1Þ
where the total energy of fluid flow is expressed as
e¼hþ
u2 2
ð2Þ
The mass balance is introduced in the form
@ q @ðquÞ ¼0 þ @x @t
ð3Þ
and the momentum balance reads
@ðquÞ @ðqu2 Þ @p þ ¼ sw qg sin h @t @x @x
ð4Þ
Introducing Eqs. (2)–(4) into (1) it is obtained
q_ Dh 1 Dp ¼ þ usw þ Dt q Dt q
ð5Þ
where the substantial derivative is
D @ @ ¼ þu Dt @t @x
ð6Þ
Eq. (5) is applied to transient hot water flows in district heating systems, such as the system start-ups and temperature fronts propagation towards consumers, the changes of the heat source or consumption load, or the system shut-down. For such scenarios, the influence of the pressure change or the dissipation heat due to wall friction (the first and the second terms on the r.h.s. of Eq. (5)) can be neglected in comparison to the enthalpy changes due to water heating or cooling in heat power plants and consumer substations (the third term on the r.h.s. of Eq. (5)). Hence, Eq. (5) is transformed in the following form:
Dh q_ ¼ Dt q
ð7Þ
2.2. Numerical method Eq. (7) is the partial differential equation of the hyperbolic type. It is solved along the characteristic path that is determined by the
2169
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
fluid particle motion. In the space–time coordinate system the characteristic path is defined with
dt 1 ¼ dx u
ð8Þ
ui > 0 and nodes (i 1), i, (i + 1), (i + 2) in case of ui < 0. The suitability of the third degree of the LIP is demonstrated in [9]. The following calculation algorithm, suitable for the computer programming of the LIP, is applied: m X hj ; D j j¼0
where u denotes the fluid particle velocity. Material derivative D/Dt in Eq. (7) along the characteristic path defined with Eq. (8) transforms into the time derivative and the energy Eq. (7) has the following form:
hm ðxA Þ ¼ Pmþ1 ðxA Þ
dh q_ ¼ dt q
Pmþ1 ðxA Þ ¼ ðxA x0 Þ . . . ðxA xm Þ; and
ð18Þ
Dj ¼ ðxj x0 Þðxj x1 Þ . . . ðxj xj1 ÞðxA xj Þðxj xjþ1 Þ . . . ðxj xm Þ;
hi ðt þ DtÞ hA ðtÞ ¼
q_
q
xi xA ¼ uA Dt uA ui x A xi ¼ ui1 ui xi1 xi uA ui x A xi ¼ uiþ1 ui xiþ1 xi
ð10Þ
A
ð11Þ dt 1 ¼ >0 dx ui dt 1 for ¼ <0 dx ui for
ð12Þ ð13Þ
The expression for the determination of point A coordinate is derived from Eqs. (11)–(13)
u i Dt 1þa
ð19Þ
dt 1 ¼ >0 dx ui dt 1 for ¼ <0 dx ui for
2
3
4
M mþ1 jPmþ1 ðxÞj; ðm þ 1Þ! ðmþ1Þ
where M mþ1 ¼ max jh a6x6b
ðxÞj
ð20Þ ð21Þ
According to Eqs. (18), (20), and (21) the application of the LIP for the enthalpy interpolation on the grid presented in Fig. 1 gives the following truncation error (TE) for the numerical discretization of the enthalpy along the x coordinate, TE = O[(Dx)m+1]. Therefore, for the third degree LIP the TE is of the fourth order, O[(Dx)4]. Time integration of the energy Eq. (9) is performed along the characteristic path with the Euler explicit method, Eq. (10). Therefore, time integration method bears the first order TE = O(Dt). Stability analyses shows that a solution along a characteristic paths requires a time step of integration which satisfies the Courant criterion [11]
Dx ; jui j
Dt 6 min
i ¼ 2; 3; . . . ; n þ 1:
ð22Þ
ð15Þ ð16Þ
2.3. Boundary condition
For the calculation of enthalpy in node i at time level t + Dt, i.e. hi(t + Dt), from Eq. (10), the initial enthalpy value at point A, hA(t), is needed. Enthalpy hA(t) is determined with the Lagrange’s interpolation polynomial (LIP) of the third degree, where the enthalpy values at the nodes (i 2), (i 1), i and (i + 1) are used in case of
1
jhðxÞ hm ðxÞj 6
Volumetric heat flux q_ A at A is specified or, as the enthalpy (i.e. temperature) at point A is known, it may be calculated from the laws of heat transfer from the fluid, through the pipe wall and insulation, towards the surroundings. Density at point A may be determined from the equation of state in the form q = q(h, p).
ð14Þ
where
ui ui1 Dt xi xi1 uiþ1 ui Dt a¼ xiþ1 xi
j ¼ 0; 1; . . . ; m: The error of the interpolation with the LIP is given in [10]
Dt
where index i denotes a numerical node, as presented in Fig. 1, while point A is the intersection of the characteristic path defined with Eq. (8) and the x coordinate at instant t. Coordinate xA of the point A is determined using the slope of the characteristic path and the linear interpolation of velocity between nodes i and i 1 for positive flow direction u > 0, as well as i and i + 1 for negative flow direction u < 0. The following expressions are written:
a¼
where
ð9Þ
The time derivative in Eq. (9) is approximated with the finite difference and the algebraic equation is obtained
xA ¼ xi
ð17Þ
i-2
Connections of consumers with heat supply and return pipelines within a district heating system network are made through junctions of three or more pipelines, as it is illustrated in Fig. 2. The pipelines with or without pumps and heat exchangers that conduct water towards the junction are denoted with Ji
i-1
i
i+1
t+dt t t i-2
i-1 A
x Fig. 1. Numerical grid.
i
i+1
n
n+1
2170
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
I1 Ij (D’,J1)
D
(D’,Ij)
J1
Im
Ji Jn
Fig. 2. Several pipes in a junction.
(i = 1, 2, . . ., n), while the pipelines that conduct the water from the junction are denoted with Ij (j = 1, 2, . . ., m). The junction is denoted with D, while the inlets and outlets of pumps and heat exchangers are denoted with D’ and the corresponding pipe number Ji or Ij. The enthalpy change from the pump or heat exchanger inlet to outlet is denoted with DhJi or DhIj, where indices correspond to the notation of junction inflow or outflow pipe. The enthalpy change along the characteristic path that reaches point D’ at instant t + Dt in the pipeline Ji is given with
hD0 ;Ji ðt þ DtÞ ¼ hA;Ji ðtÞ þ
q_
q
Dt i ¼ 1; 2; . . . ; n:
ð23Þ
A;Ji
The enthalpy in point D is determined by the mixing of water flows in pipelines Ji (i = 1, 2, . . ., n); hence, the steady-state energy balance of water flows from points D’ in all pipelines Ji (i = 1, 2, . . ., n) to mixing point D, including enthalpy changes in heat exchangers and pumps, gives
Pn hD ¼
i¼1
qD0 ;Ji uD0 ;Ji hD0 ;Ji þ DhJi AJi Pn 0 i¼1 qD0 ;Ji uD ;Ji AJi
ð24Þ
where the denominator in Eq. (24) represents the total mass flow rate towards the junction D. The cross sectional area of pipeline Ji is denoted with AJi, while the water density at point D’ in the pipeline Ji is denoted with qD’,Ji. The energy balance for a heat exchanger or a pump in the outflowing pipeline Ij (j = 1, 2, . . ., m) gives
hD0 ;Ij ¼ hD þ DhIj ;
j ¼ 1; 2; . . . ; m
ð25Þ
2.4. Hydraulic model Thermal transients in a district heating system are slow; they last for tens of minutes or for hours. During these time periods, the water in the pipelines can be considered to be incompressible in regard to the pressure changes because the sonic velocity of pressure wave propagations is of the order of 1000 m/s (with this velocity a perturbation of pressure will be transported and attenuated within the network in much shorter time periods than the periods of the thermal transients). This assumption permits calculating the network hydraulics as a succession of steady states by utilizing the steady-state hydraulics models of network flow. The steady-state hydraulic model, which is used for the calculation of pressure and flow rates distribution within the pipeline network in this paper, is based on the loop model of the network, and it is presented in details in [12]. 3. Results and discussion The developed model is used for the computer simulations of real district heating system transients that were measured at the district heating system in Zemun, Serbia, in February 2007. The installed power of the heat plant is 65 MW. About 2/3 of the heat is delivered to consumers within the main line 3, presented in Fig. 3,
Fig. 3. The district heating system in Zemun.
which comprises more than 50 substations, mainly in resident buildings. Substations are equipped with flow control valves that adjust the hot water flow rate according to the consumers’ needs. Temperatures in the supply and return lines at the heat plant are measured, as well as in three consumer substations: SST A, SST B and SST C (denoted with 117, 108 and 206 in Fig. 3). The substations SST A and SST C are located in large residential buildings, while the SST B is located in the sport centre. Lengths of distribution pipelines from the heat plant to the substations SST A, SST B and SST C are 729 m, 1615 m, and 1718 m, respectively. Also, measured are the hot water flow rates through the substations and the total flow rate of supply water from the heat plant towards the main line. The measuring instruments of standard industrial accuracy are used (the relative error of the ultrasonic flow measurement is ±2%, while the absolute error of the temperature measurement is up to 1 °C). Fig. 4 shows the measured temperature changes in supply and return lines in case of boiler load increase. The substation SST A consumes heat only for space heating, while the SST B and SST C consume heat both for space heating and sanitary hot water preparation. The dynamic temperature changes at the substation C outlet are the result of the sanitary hot water consumption. The boiler load increases within the period from the 14th to 32nd min (Fig. 4), and consequently, temperature front propagates to the consumers with the peak temperature of 70 °C. The temperature wave arrives in the nearest substation SST A, then in SST B, and finally in the farthest SST C. The time necessary for temperature front propagation to substations is determined by the distance from the heat station and the hot water velocity. The maximum value of the temperature front that arrives in the substations is lower for approximately 2 °C. Temperature values in the return lines practically do not show the influence of the boiler load increase, which indicates the large accumulation capacity of the heated objects. The temper-
2171
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
Temperature (ºC)
a
The total flow rate from the heat plant towards the consumers presented in Fig. 3 is shown in Fig. 5. A slight volumetric flow rate increase at the beginning of the transient is caused by the supply water temperature increase, consequent water density decrease and a reduction of the total pressure drop along the supply lines. In the latter period the flow rate decreases since the control valves in the substations close due to the increased supply water temperature, and the total resistance to the water flow through the network increases. Fig. 6 shows the hot water flow rates through three substations of which measurements were performed. The
70 65 60 55 50 45 40 35 30 25 20 0
100
120
140
0.12 0.1 0.08 0.06 0.04 0.02 0 0
20
40
60 80 Time (min)
100
120
20
80 100 Time (min)
120
140
160
0.003 SST B
0.0025
20
40
60 80 Time (min)
100
120
0.002 0.0015 SST A
0.001
SST C
0.0005
140
0 0
70 65 60 55 50 45 40 35 30 25 20
20
40
60
80 100 Time (min)
120
140
160
Fig. 6. Flow rates through the substations during the boiler load increase.
70 65 o
Temperature ( C)
Temperature (ºC)
60
Fig. 5. Flow rate measured at the heat plant outlet during the boiler load increase.
70 65 60 55 50 45 40 35 30 25 20 0
d
40
140
Flow rate (m3/s)
Temperature (ºC)
60 80 Time (min)
70 65 60 55 50 45 40 35 30 25 20 0
c
40
Flow rate (m3/s)
Temperature (ºC)
b
20
0
20
40
60 80 Time (min)
100
120
140
Fig. 4. Measured supply and return temperatures during the boiler load increase (a) at the heat plant (top), (b) at the substation SST A (the second from the top), (c) at the substation SST B (the third from the top), (d) at the substation SST C (bottom).
ature increase is recorded only in SST A after 100 min, but this is caused by the total stopping of the hot water flow through the substation, since there was no further need for space heating due to the rise of the outside air temperature in the latter period of the transient.
60 55 50
Heat plant, measured SST A, simulation SST B, simulation SST C, simulation SST A, measured SST B, measured SST C, measured
45 40
0
20
40 60 Time (min)
80
100
Fig. 7. Temperature waves propagation to the consumers – comparison of the calculated and measured results for the transient of boiler load increase.
2172
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
substation SST A supplies energy only for the residential building heating, hence, after 98 min the control valve in SST A is closed and the supply water flow is stopped due to the high outside air temperature. The water flows through SST B and SST C are reduced
Temperature (ºC)
b
Temperature (ºC)
c
Temperature (ºC)
d
70 65 60 55 50 45 40 35 30 25 20 0
50
100 Time (min)
150
0
50
100 Time (min)
150
0
50
100 Time (min)
150
70 65 60 55 50 45 40 35 30 25 20
70 65 60 55 50 45 40 35 30 25 20
70 65 60 55 50 45 40 35 30 25 20
0.16 0.14 Flow rate (m3/s)
Temperature (ºC)
a
at the end of the transient due to the supply water increase. These water flow reductions are delayed in respect to the arrivals of the temperature fronts, since the action of control valves in substations is delayed according to the built in control logic. The flow rate changes in SST B and SST C are more dynamic compared to SST A since the former two substations supply heat for the preparation of the hot sanitary water besides space heating, while, as mentioned, SST A supplies energy only for space heating. The measured and calculated temperature fronts at substations SST A, SST B and SST C are compared in Fig. 7. The good agreement is obtained, despite the fact that the total supply flow rate from the heat plant is kept constant throughout the whole transient (the value of 0.09 m3/s is adopted according to Fig. 5), as well as the flow rates through the substations are kept constant (the values calibrated to the initial steady-state condition of the district heating system). The shapes and amplitudes of the calculated temperature fronts are preserved, which indicate that there is no smearing of the calculated temperature fronts due to the numerical diffusion. This feature of the numerical method is due to the applied numerical scheme of the third-order in space for the solution of energy equation. The application of the lower-order numerical scheme would lead to the front smearing, as demonstrated in [9]. The second set of the measured data shows temperature and flow rate changes at substations in case of boiler load reduction. Measurement started at 11:00 am, and the outside air temperature was +12 °C. The air temperature was increasing, so, at the end of the measurement period it was 14 °C. The measured supply and return temperatures are shown in Fig. 8. The boiler load decreases in the period from the 58th till the 100th min. The hot water temperature drop occurs in the substation SST A first, 24 min after the beginning of the boiler heat load decrease. The temperature drop in the substations SST C and SST B begins after 48 and 52 min, respectively. The supply water temperature at the heat plant drops to 42 °C after 86 min. The temperature drop front reaches the substation SST A first, and then substations SST C and SST B, after time intervals determined by the distance from the heat station and the hot water velocity. The return temperatures from the substations with the sanitary water heating drop within a few minutes after the arrival of the low temperature front in the substation supply line; while in substations where only space heating exists, this temperature drop occurs after tens of minutes (the temperature drop in the substation SST A occurs 26 min after the arrival of the cold temperature wave, the delay that is determined by the residential building transient heat behaviour). The supply water temperature increase at the heat station begins 100 min after the transient initial time (Fig. 8), while the increase of supply water temperatures at the substations occur after the same time intervals as required for the temperature drop front arrival, as described above. The total water flow rate in the supply line from the heat plant is shown in Fig. 9. The flow rate increase is observed after 106 min, and it is caused by the opening of the flow control valves in substa-
0
50
100 Time (min)
150
0.12 0.1 0.08 0.06 0.04 0.02 0 0
Fig. 8. Measured supply and return temperatures during the boiler load reduction (a) at the heat plant (top), (b) at the substation SST A (the second from the top), (c) at the substation SST B (the third from the top), (d) at the substation SST C (bottom).
50
100 Time (min)
150
Fig. 9. Flow rate measured at the heat plant outlet during the boiler load reduction.
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
0.006
Flow rate (m3/s)
0.005 0.004
SST A
0.003
SST B
0.002
SST C
0.001 0 0
50
100 Time (min)
150
Fig. 10. Flow rates through substations during the boiler load reduction.
70
2173
ing systems, as well as the results of measurements and computer simulations of two thermal transients caused by the rapid increase and decrease of the heat power plant load. The comparison of the calculated and measured temperature fronts at three consumer substations in different parts of the district heating system network and at different distances from the heat plant shows good agreement. The shapes and amplitudes of the calculated temperature fronts are preserved, which indicate that there is no smearing of the calculated fronts due to the numerical diffusion. This feature of the numerical method is due to the applied numerical scheme of the third-order in space for the solution of energy equation. The performed research is a part of the verification of the developed numerical method, which will be used for optimization of operational regimes of district heating systems, with the goal to increase the overall energy efficiency of the heat supply to consumers in transient conditions. Future research will include a modeling of the consumers’ substations dynamic behaviour and development of algorithms for the optimization of the district heating systems control in various transient conditions.
o
Temperature ( C)
65 60
Acknowledgements
55
This work was supported by the Ministry of Science of the Republic of Serbia through the National Energy Efficiency Program (Grant 242008) and the Fundamental Research Program (Grant 144022). The authors are thankful to Maja Todorovic and Radoslav Galic (University of Belgrade) and Dragan Mandic, Dragan Dragojevic and Srdjan Nikodijevic (JKP Beogradske elektrane) for their aid during measurements at the District Heating System of Zemun.
50 45 Heat plant, measured SST A, simulated SST B, simulated SST C, simulated SST A, measured SST B, measured SST C, measured
40 35 30
40
60
80
100 Time (min)
120
140
160
Fig. 11. Temperature waves propagation to the consumers – comparison of the calculated and measured results for the transient of the boiler load reduction.
tions due to the arrival of the temperature drop fronts. Hot water flow rates through the substations are shown in Fig. 10. The increase of the hot water flow in the substation SST A occurs after 92 min (10 min after the supply water temperature at the substation starts to decrease). The maximum flow rate at this substation is achieved at the 110th min. Although the supply temperature decreases till the 116th min, the further flow rate increase through SST A is not possible since the flow control valve is already fully open. Substantial changes of water flow rates in other two substations are not observed. Fig. 11 shows the measured and calculated temperature fronts at substations SST A, SST B and SST C in case of heat plant load reduction. Again, the good agreement is obtained, although the total supply flow rate from the heat plant is kept constant throughout the whole transient (the value of 0.11 m3/s is adopted according to Fig. 9), and the flow rates through the substations are kept constant (the values calibrated to the initial steady-state condition of the district heating system). Also in this case the shapes and amplitudes of the calculated temperature fronts are perfectly preserved. 4. Conclusions This paper presents the model and numerical method for the computer prediction of the thermal transients in the district heat-
References [1] Kunz VJ, Haldi PA, Sarlos G. Dynamic behaviour of district heating networks. Fernwaerme Int 1991;20:104–19. [2] Larsen HV, Palsson H, Bohm B, Ravn HF. Aggregated dynamic simulation model of district heating networks. Energy Convers Manage 2002;43:995–1019. [3] Larsen HV, Bohm B, Wigbels M. A comparison of aggregated models for simulation and operational optimization of district heating networks. Energy Convers Manage 2004;45:1119–39. [4] Benonysson A, Bøhm B, Ravn HF. Operational optimization in a district heating system. Energy Convers Manage 1995;36:297–314. [5] Park YS, Kim WT, Kim BK, Ha SK. Simple models for operational optimization, Chapter 2: state of the art. In: Bøhm B, editor. IEA district heating and cooling, annex VI, Report 2002: S1. Technical University of Denmark, Kgs. Lyngby; 2002, p. 7–10. [6] Gabrielaitiene I, Bøhm B, Sunden B. Modelling temperature dynamics of a district heating system in Naestved, Denmark – a case study. Energy Convers Manage 2007;48:78–86. [7] Iversen S, Ougaard P, Loppenthien JK. Dynamic temperature optimizationproviding instant results. EuroHeat 2006;3-II:46–9. [8] Karamarkovic, V, Solujic A, Simovic, J, Ramic B. Implementation program of energy development strategy in the field of industrial energy. In: The Society of Thermal Engineers of Serbia edition. Proceedings of the regional conference industrial energy and environmental protection in Southeast Europe, Zlatibor, Serbia, 24–28 June, 2008. [9] Stevanovic V, Jovanovic Z. A hybrid method for the numerical prediction of enthalpy transport in fluid flow. Int Commun Heat Mass Transfer 2000;27:23–34. [10] Demidovitch B, Maron I. Elements de Calcul Numerique. Moscou: Editions Mir; 1987. [11] Tannehill JC, Anderson DA, Pletcher RH. Computational fluid mechanics and heat transfer. Bristol: Taylor & Francis; 1997. [12] Stevanovic V, Prica S, Maslovaric B, Zivkovic B, Nikodijevic S. Efficient numerical method for district heating system hydraulics. Energy Convers Manage 2007;48:1536–43.
Contents lists available at ScienceDirect
Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman
Prediction of thermal transients in district heating systems Vladimir D. Stevanovic a,*, Branislav Zivkovic a, Sanja Prica a, Blazenka Maslovaric a, Vladan Karamarkovic b, Vojin Trkulja c a
Faculty of Mechanical Engineering, University of Belgrade, Kraljice Marije 16, 11120 Belgrade 35, Serbia Faculty of Mechanical Engineering, Dositejeva 19, 36000 Kraljevo, Serbia c JKP Beogradske Elektrane, Savski Nasip 11, 11070 Belgrade, Serbia b
a r t i c l e
i n f o
Article history: Received 3 January 2008 Received in revised form 21 September 2008 Accepted 27 April 2009 Available online 27 May 2009 Keywords: District heating system Thermal transients Simulation
a b s t r a c t A prediction of thermal transients in district heating systems is important in order to adjust in an energy efficient manner the loads of heat power plants and pump stations with dynamic consumption needs during changes of the outside air temperature, the wind intensity, the solar radiation, or system startups and shut-downs. A model and corresponding computer code are developed with the aim to simulate the thermal transients in district heating systems. They are based on the high-order accurate numerical solution of the transient energy equation, and the hydraulic prediction of pressure and fluid flow rates within the complex pipe network. Thermal transients caused by an increase and decrease of the heat power plant load are simulated for real operating conditions of the district heating system. Predicted temperature front propagations show a good agreement with data measured at three consumer substations located in different parts of the district heating network and at different distances from the heat plant. The developed computational tool provides reliable information about time periods of temperature fronts propagation and heat distribution from the heat source to consumers within the whole district heating network. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Managing a district heating system operation during temperature transients is important for the reliable supply of adequate heat amounts to consumers in an energy efficient manner. The temperature transients occur during changes of weather conditions, such as a change of temperature or wind intensity, during start-ups and shut-downs of district heating system operation, or due to the change of hot water consumption in the systems that provide sanitary water heating. The efficiency of heat supply during transients depends on the timing of heat source power according to the heat consumption, on the time periods necessary for the heat transport from the heat source plant to consumers’ substations, as well as on the heat accumulation within the distribution pipelines and heat losses to the surroundings. Due to the complexity of district heating networks, and dynamic characteristics of heat consumption and transport properties, the prediction of transient regimes is possible only with the application of numerical models and computer codes. One of the early dynamic models for the transient temperature predictions in the complex district heating systems was presented almost two decades ago [1]. The model was based on the lumped * Corresponding author. Tel.: +381 11 3370 561; fax: +381 11 3370 364. E-mail address: [email protected] (V.D. Stevanovic). 0196-8904/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.enconman.2009.04.034
parameter energy model of the district heating pipeline network. In this approach the whole pipe is described with the single control volume, the linear water temperature change along the pipe is assumed and the mean water temperature is assigned to the pipe as the simple average of the input and output water pipe temperatures. Such an approach is not able to predict the temperature wave propagation along the pipe. Hence, a large diffusion of the temperature waves propagation within the pipeline network is expected, together with a low accuracy prediction of time delays of temperature waves propagations from a heat source to consumers. Due to the complexity of pipeline networks of district heating systems, the aggregated dynamic models have been developed [2,3], which represent a number of pipeline segments with only one equivalent pipe. Aggregated models are simpler and they provide faster computer simulations, but they are not able at all to predict temperature waves propagations from heat sources to consumers within the pipeline network. A further development towards the real physical representation of temperature transients was achieved with the node method [4,5]. The principle of the node method is that it keeps trace of how long time a water mass, which in the actual time step arrives at a node, has been on its way from the last node. Based on time series for the temperature history of the different nodes the temperature of the mass at the inlet of the pipe is calculated, and on the basis of this as well as the heat loss and heat capacity of the pipe, the temperature of the mass
2168
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
Nomenclature A e g h q_ p t
cross sectional area, m2 fluid energy, J/kg gravity, m/s2 specific enthalpy, J/kg volumetric heat source, W/m3 pressure, Pa time, s, min
at the outlet of the pipe is evaluated. In this manner a new temperature for all nodes in the network can be calculated in every time step. As in the case of the lumped parameter model [1], the node method does not solve the temperature wave propagation along the pipe; hence, the temperature wave diffusion and smearing could be expected. The real physical model of the temperature waves propagation within the pipeline network of a district heating system will assume a direct numerical solving of the energy balance equation of water flow along all pipes. It is reported that such an approach is applied recently in the commercial computer code TERMIS [6,7], but no information about the numerical solving procedure is available in the open literature. The developed node method and the TERMIS code are tested for the conditions of temperature fronts propagations in parts of district heating systems [6]. At present, numerical models are applied at several district heating systems in order to optimize transient operational regimes; first results show the saving of the total fuel energy consumption of 2% at least [7,8]. Dynamics of pressure and flow changes within district heating networks are determined by the pressure waves propagation. The speed of pressure waves in the hot water of district heating systems is approximately 1000 m/s, which is 1000 times faster than the water velocity that is in the range from below 1 m/s till 2 or 3 m/s (these higher velocities are encountered in main transport district heating pipelines). The temperature variations propagate with the velocity of the water flow. Hence, the dynamics of the flow in the network are of minor importance compared with the dynamics of the temperature changes, from an operational optimization point of view. On this basis, the modeling approaches [1–6] for the prediction of temperature transients applied the steady-state model of water flow, where the pressure and flow distributions are calculated for every time instant of the numerical integration, based on the actual distribution of temperature and corresponding thermo-physical characteristics of the water in the pipeline network. This paper presents the new numerical model for the simulation and analyses of the thermal transients in the district heating system. It is based on the numerical solution of the transient energy equation with the numerical scheme of the third-order accuracy in space, and on the steady-state hydraulic solution of pressure changes and flow rates within the network. The developed method is applied for the simulation of two thermal transients in the district heating system of Zemun in Serbia, caused by an increase and decrease of the heat power plant load. The obtained numerical results are compared with data measured at three consumers substations located at different distances and in different parts of the pipeline network. The obtained good agreement shows that the developed numerical method is a reliable tool for the simulation and analyses of district heating system thermal transients. The presented modeling approach is an advantage in comparison to the lumped parameter, aggregated and node methods [1–4], since it provides the prediction of distributed temperatures along all pipes within the district heating system, and enables the accurate calculation of dynamic temperature changes without the influence of the numerical diffusion in predictions of temperature waves propagations. A com-
u x
velocity, m/s coordinate, m
Greek symbols h pipe inclination, rad q density, kg/m3 sw wall friction per channel unit volume, N/m3
parison with the recently introduced temperature dynamic model of the TERMIS code could not be done, since the TERMIS numerical model is not reported in the open literature. Also, here presented experimental data on temperature waves propagations in the complex district heating system is a contribution to the limited information about this process. 2. Model development 2.1. Transient energy balance Transient one-dimensional energy balance equation is written as
@ðqeÞ @ðqueÞ @p þ ¼ qgu sin h þ q_ @t @t @t
ð1Þ
where the total energy of fluid flow is expressed as
e¼hþ
u2 2
ð2Þ
The mass balance is introduced in the form
@ q @ðquÞ ¼0 þ @x @t
ð3Þ
and the momentum balance reads
@ðquÞ @ðqu2 Þ @p þ ¼ sw qg sin h @t @x @x
ð4Þ
Introducing Eqs. (2)–(4) into (1) it is obtained
q_ Dh 1 Dp ¼ þ usw þ Dt q Dt q
ð5Þ
where the substantial derivative is
D @ @ ¼ þu Dt @t @x
ð6Þ
Eq. (5) is applied to transient hot water flows in district heating systems, such as the system start-ups and temperature fronts propagation towards consumers, the changes of the heat source or consumption load, or the system shut-down. For such scenarios, the influence of the pressure change or the dissipation heat due to wall friction (the first and the second terms on the r.h.s. of Eq. (5)) can be neglected in comparison to the enthalpy changes due to water heating or cooling in heat power plants and consumer substations (the third term on the r.h.s. of Eq. (5)). Hence, Eq. (5) is transformed in the following form:
Dh q_ ¼ Dt q
ð7Þ
2.2. Numerical method Eq. (7) is the partial differential equation of the hyperbolic type. It is solved along the characteristic path that is determined by the
2169
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
fluid particle motion. In the space–time coordinate system the characteristic path is defined with
dt 1 ¼ dx u
ð8Þ
ui > 0 and nodes (i 1), i, (i + 1), (i + 2) in case of ui < 0. The suitability of the third degree of the LIP is demonstrated in [9]. The following calculation algorithm, suitable for the computer programming of the LIP, is applied: m X hj ; D j j¼0
where u denotes the fluid particle velocity. Material derivative D/Dt in Eq. (7) along the characteristic path defined with Eq. (8) transforms into the time derivative and the energy Eq. (7) has the following form:
hm ðxA Þ ¼ Pmþ1 ðxA Þ
dh q_ ¼ dt q
Pmþ1 ðxA Þ ¼ ðxA x0 Þ . . . ðxA xm Þ; and
ð18Þ
Dj ¼ ðxj x0 Þðxj x1 Þ . . . ðxj xj1 ÞðxA xj Þðxj xjþ1 Þ . . . ðxj xm Þ;
hi ðt þ DtÞ hA ðtÞ ¼
q_
q
xi xA ¼ uA Dt uA ui x A xi ¼ ui1 ui xi1 xi uA ui x A xi ¼ uiþ1 ui xiþ1 xi
ð10Þ
A
ð11Þ dt 1 ¼ >0 dx ui dt 1 for ¼ <0 dx ui for
ð12Þ ð13Þ
The expression for the determination of point A coordinate is derived from Eqs. (11)–(13)
u i Dt 1þa
ð19Þ
dt 1 ¼ >0 dx ui dt 1 for ¼ <0 dx ui for
2
3
4
M mþ1 jPmþ1 ðxÞj; ðm þ 1Þ! ðmþ1Þ
where M mþ1 ¼ max jh a6x6b
ðxÞj
ð20Þ ð21Þ
According to Eqs. (18), (20), and (21) the application of the LIP for the enthalpy interpolation on the grid presented in Fig. 1 gives the following truncation error (TE) for the numerical discretization of the enthalpy along the x coordinate, TE = O[(Dx)m+1]. Therefore, for the third degree LIP the TE is of the fourth order, O[(Dx)4]. Time integration of the energy Eq. (9) is performed along the characteristic path with the Euler explicit method, Eq. (10). Therefore, time integration method bears the first order TE = O(Dt). Stability analyses shows that a solution along a characteristic paths requires a time step of integration which satisfies the Courant criterion [11]
Dx ; jui j
Dt 6 min
i ¼ 2; 3; . . . ; n þ 1:
ð22Þ
ð15Þ ð16Þ
2.3. Boundary condition
For the calculation of enthalpy in node i at time level t + Dt, i.e. hi(t + Dt), from Eq. (10), the initial enthalpy value at point A, hA(t), is needed. Enthalpy hA(t) is determined with the Lagrange’s interpolation polynomial (LIP) of the third degree, where the enthalpy values at the nodes (i 2), (i 1), i and (i + 1) are used in case of
1
jhðxÞ hm ðxÞj 6
Volumetric heat flux q_ A at A is specified or, as the enthalpy (i.e. temperature) at point A is known, it may be calculated from the laws of heat transfer from the fluid, through the pipe wall and insulation, towards the surroundings. Density at point A may be determined from the equation of state in the form q = q(h, p).
ð14Þ
where
ui ui1 Dt xi xi1 uiþ1 ui Dt a¼ xiþ1 xi
j ¼ 0; 1; . . . ; m: The error of the interpolation with the LIP is given in [10]
Dt
where index i denotes a numerical node, as presented in Fig. 1, while point A is the intersection of the characteristic path defined with Eq. (8) and the x coordinate at instant t. Coordinate xA of the point A is determined using the slope of the characteristic path and the linear interpolation of velocity between nodes i and i 1 for positive flow direction u > 0, as well as i and i + 1 for negative flow direction u < 0. The following expressions are written:
a¼
where
ð9Þ
The time derivative in Eq. (9) is approximated with the finite difference and the algebraic equation is obtained
xA ¼ xi
ð17Þ
i-2
Connections of consumers with heat supply and return pipelines within a district heating system network are made through junctions of three or more pipelines, as it is illustrated in Fig. 2. The pipelines with or without pumps and heat exchangers that conduct water towards the junction are denoted with Ji
i-1
i
i+1
t+dt t t i-2
i-1 A
x Fig. 1. Numerical grid.
i
i+1
n
n+1
2170
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
I1 Ij (D’,J1)
D
(D’,Ij)
J1
Im
Ji Jn
Fig. 2. Several pipes in a junction.
(i = 1, 2, . . ., n), while the pipelines that conduct the water from the junction are denoted with Ij (j = 1, 2, . . ., m). The junction is denoted with D, while the inlets and outlets of pumps and heat exchangers are denoted with D’ and the corresponding pipe number Ji or Ij. The enthalpy change from the pump or heat exchanger inlet to outlet is denoted with DhJi or DhIj, where indices correspond to the notation of junction inflow or outflow pipe. The enthalpy change along the characteristic path that reaches point D’ at instant t + Dt in the pipeline Ji is given with
hD0 ;Ji ðt þ DtÞ ¼ hA;Ji ðtÞ þ
q_
q
Dt i ¼ 1; 2; . . . ; n:
ð23Þ
A;Ji
The enthalpy in point D is determined by the mixing of water flows in pipelines Ji (i = 1, 2, . . ., n); hence, the steady-state energy balance of water flows from points D’ in all pipelines Ji (i = 1, 2, . . ., n) to mixing point D, including enthalpy changes in heat exchangers and pumps, gives
Pn hD ¼
i¼1
qD0 ;Ji uD0 ;Ji hD0 ;Ji þ DhJi AJi Pn 0 i¼1 qD0 ;Ji uD ;Ji AJi
ð24Þ
where the denominator in Eq. (24) represents the total mass flow rate towards the junction D. The cross sectional area of pipeline Ji is denoted with AJi, while the water density at point D’ in the pipeline Ji is denoted with qD’,Ji. The energy balance for a heat exchanger or a pump in the outflowing pipeline Ij (j = 1, 2, . . ., m) gives
hD0 ;Ij ¼ hD þ DhIj ;
j ¼ 1; 2; . . . ; m
ð25Þ
2.4. Hydraulic model Thermal transients in a district heating system are slow; they last for tens of minutes or for hours. During these time periods, the water in the pipelines can be considered to be incompressible in regard to the pressure changes because the sonic velocity of pressure wave propagations is of the order of 1000 m/s (with this velocity a perturbation of pressure will be transported and attenuated within the network in much shorter time periods than the periods of the thermal transients). This assumption permits calculating the network hydraulics as a succession of steady states by utilizing the steady-state hydraulics models of network flow. The steady-state hydraulic model, which is used for the calculation of pressure and flow rates distribution within the pipeline network in this paper, is based on the loop model of the network, and it is presented in details in [12]. 3. Results and discussion The developed model is used for the computer simulations of real district heating system transients that were measured at the district heating system in Zemun, Serbia, in February 2007. The installed power of the heat plant is 65 MW. About 2/3 of the heat is delivered to consumers within the main line 3, presented in Fig. 3,
Fig. 3. The district heating system in Zemun.
which comprises more than 50 substations, mainly in resident buildings. Substations are equipped with flow control valves that adjust the hot water flow rate according to the consumers’ needs. Temperatures in the supply and return lines at the heat plant are measured, as well as in three consumer substations: SST A, SST B and SST C (denoted with 117, 108 and 206 in Fig. 3). The substations SST A and SST C are located in large residential buildings, while the SST B is located in the sport centre. Lengths of distribution pipelines from the heat plant to the substations SST A, SST B and SST C are 729 m, 1615 m, and 1718 m, respectively. Also, measured are the hot water flow rates through the substations and the total flow rate of supply water from the heat plant towards the main line. The measuring instruments of standard industrial accuracy are used (the relative error of the ultrasonic flow measurement is ±2%, while the absolute error of the temperature measurement is up to 1 °C). Fig. 4 shows the measured temperature changes in supply and return lines in case of boiler load increase. The substation SST A consumes heat only for space heating, while the SST B and SST C consume heat both for space heating and sanitary hot water preparation. The dynamic temperature changes at the substation C outlet are the result of the sanitary hot water consumption. The boiler load increases within the period from the 14th to 32nd min (Fig. 4), and consequently, temperature front propagates to the consumers with the peak temperature of 70 °C. The temperature wave arrives in the nearest substation SST A, then in SST B, and finally in the farthest SST C. The time necessary for temperature front propagation to substations is determined by the distance from the heat station and the hot water velocity. The maximum value of the temperature front that arrives in the substations is lower for approximately 2 °C. Temperature values in the return lines practically do not show the influence of the boiler load increase, which indicates the large accumulation capacity of the heated objects. The temper-
2171
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
Temperature (ºC)
a
The total flow rate from the heat plant towards the consumers presented in Fig. 3 is shown in Fig. 5. A slight volumetric flow rate increase at the beginning of the transient is caused by the supply water temperature increase, consequent water density decrease and a reduction of the total pressure drop along the supply lines. In the latter period the flow rate decreases since the control valves in the substations close due to the increased supply water temperature, and the total resistance to the water flow through the network increases. Fig. 6 shows the hot water flow rates through three substations of which measurements were performed. The
70 65 60 55 50 45 40 35 30 25 20 0
100
120
140
0.12 0.1 0.08 0.06 0.04 0.02 0 0
20
40
60 80 Time (min)
100
120
20
80 100 Time (min)
120
140
160
0.003 SST B
0.0025
20
40
60 80 Time (min)
100
120
0.002 0.0015 SST A
0.001
SST C
0.0005
140
0 0
70 65 60 55 50 45 40 35 30 25 20
20
40
60
80 100 Time (min)
120
140
160
Fig. 6. Flow rates through the substations during the boiler load increase.
70 65 o
Temperature ( C)
Temperature (ºC)
60
Fig. 5. Flow rate measured at the heat plant outlet during the boiler load increase.
70 65 60 55 50 45 40 35 30 25 20 0
d
40
140
Flow rate (m3/s)
Temperature (ºC)
60 80 Time (min)
70 65 60 55 50 45 40 35 30 25 20 0
c
40
Flow rate (m3/s)
Temperature (ºC)
b
20
0
20
40
60 80 Time (min)
100
120
140
Fig. 4. Measured supply and return temperatures during the boiler load increase (a) at the heat plant (top), (b) at the substation SST A (the second from the top), (c) at the substation SST B (the third from the top), (d) at the substation SST C (bottom).
ature increase is recorded only in SST A after 100 min, but this is caused by the total stopping of the hot water flow through the substation, since there was no further need for space heating due to the rise of the outside air temperature in the latter period of the transient.
60 55 50
Heat plant, measured SST A, simulation SST B, simulation SST C, simulation SST A, measured SST B, measured SST C, measured
45 40
0
20
40 60 Time (min)
80
100
Fig. 7. Temperature waves propagation to the consumers – comparison of the calculated and measured results for the transient of boiler load increase.
2172
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
substation SST A supplies energy only for the residential building heating, hence, after 98 min the control valve in SST A is closed and the supply water flow is stopped due to the high outside air temperature. The water flows through SST B and SST C are reduced
Temperature (ºC)
b
Temperature (ºC)
c
Temperature (ºC)
d
70 65 60 55 50 45 40 35 30 25 20 0
50
100 Time (min)
150
0
50
100 Time (min)
150
0
50
100 Time (min)
150
70 65 60 55 50 45 40 35 30 25 20
70 65 60 55 50 45 40 35 30 25 20
70 65 60 55 50 45 40 35 30 25 20
0.16 0.14 Flow rate (m3/s)
Temperature (ºC)
a
at the end of the transient due to the supply water increase. These water flow reductions are delayed in respect to the arrivals of the temperature fronts, since the action of control valves in substations is delayed according to the built in control logic. The flow rate changes in SST B and SST C are more dynamic compared to SST A since the former two substations supply heat for the preparation of the hot sanitary water besides space heating, while, as mentioned, SST A supplies energy only for space heating. The measured and calculated temperature fronts at substations SST A, SST B and SST C are compared in Fig. 7. The good agreement is obtained, despite the fact that the total supply flow rate from the heat plant is kept constant throughout the whole transient (the value of 0.09 m3/s is adopted according to Fig. 5), as well as the flow rates through the substations are kept constant (the values calibrated to the initial steady-state condition of the district heating system). The shapes and amplitudes of the calculated temperature fronts are preserved, which indicate that there is no smearing of the calculated temperature fronts due to the numerical diffusion. This feature of the numerical method is due to the applied numerical scheme of the third-order in space for the solution of energy equation. The application of the lower-order numerical scheme would lead to the front smearing, as demonstrated in [9]. The second set of the measured data shows temperature and flow rate changes at substations in case of boiler load reduction. Measurement started at 11:00 am, and the outside air temperature was +12 °C. The air temperature was increasing, so, at the end of the measurement period it was 14 °C. The measured supply and return temperatures are shown in Fig. 8. The boiler load decreases in the period from the 58th till the 100th min. The hot water temperature drop occurs in the substation SST A first, 24 min after the beginning of the boiler heat load decrease. The temperature drop in the substations SST C and SST B begins after 48 and 52 min, respectively. The supply water temperature at the heat plant drops to 42 °C after 86 min. The temperature drop front reaches the substation SST A first, and then substations SST C and SST B, after time intervals determined by the distance from the heat station and the hot water velocity. The return temperatures from the substations with the sanitary water heating drop within a few minutes after the arrival of the low temperature front in the substation supply line; while in substations where only space heating exists, this temperature drop occurs after tens of minutes (the temperature drop in the substation SST A occurs 26 min after the arrival of the cold temperature wave, the delay that is determined by the residential building transient heat behaviour). The supply water temperature increase at the heat station begins 100 min after the transient initial time (Fig. 8), while the increase of supply water temperatures at the substations occur after the same time intervals as required for the temperature drop front arrival, as described above. The total water flow rate in the supply line from the heat plant is shown in Fig. 9. The flow rate increase is observed after 106 min, and it is caused by the opening of the flow control valves in substa-
0
50
100 Time (min)
150
0.12 0.1 0.08 0.06 0.04 0.02 0 0
Fig. 8. Measured supply and return temperatures during the boiler load reduction (a) at the heat plant (top), (b) at the substation SST A (the second from the top), (c) at the substation SST B (the third from the top), (d) at the substation SST C (bottom).
50
100 Time (min)
150
Fig. 9. Flow rate measured at the heat plant outlet during the boiler load reduction.
V.D. Stevanovic et al. / Energy Conversion and Management 50 (2009) 2167–2173
0.006
Flow rate (m3/s)
0.005 0.004
SST A
0.003
SST B
0.002
SST C
0.001 0 0
50
100 Time (min)
150
Fig. 10. Flow rates through substations during the boiler load reduction.
70
2173
ing systems, as well as the results of measurements and computer simulations of two thermal transients caused by the rapid increase and decrease of the heat power plant load. The comparison of the calculated and measured temperature fronts at three consumer substations in different parts of the district heating system network and at different distances from the heat plant shows good agreement. The shapes and amplitudes of the calculated temperature fronts are preserved, which indicate that there is no smearing of the calculated fronts due to the numerical diffusion. This feature of the numerical method is due to the applied numerical scheme of the third-order in space for the solution of energy equation. The performed research is a part of the verification of the developed numerical method, which will be used for optimization of operational regimes of district heating systems, with the goal to increase the overall energy efficiency of the heat supply to consumers in transient conditions. Future research will include a modeling of the consumers’ substations dynamic behaviour and development of algorithms for the optimization of the district heating systems control in various transient conditions.
o
Temperature ( C)
65 60
Acknowledgements
55
This work was supported by the Ministry of Science of the Republic of Serbia through the National Energy Efficiency Program (Grant 242008) and the Fundamental Research Program (Grant 144022). The authors are thankful to Maja Todorovic and Radoslav Galic (University of Belgrade) and Dragan Mandic, Dragan Dragojevic and Srdjan Nikodijevic (JKP Beogradske elektrane) for their aid during measurements at the District Heating System of Zemun.
50 45 Heat plant, measured SST A, simulated SST B, simulated SST C, simulated SST A, measured SST B, measured SST C, measured
40 35 30
40
60
80
100 Time (min)
120
140
160
Fig. 11. Temperature waves propagation to the consumers – comparison of the calculated and measured results for the transient of the boiler load reduction.
tions due to the arrival of the temperature drop fronts. Hot water flow rates through the substations are shown in Fig. 10. The increase of the hot water flow in the substation SST A occurs after 92 min (10 min after the supply water temperature at the substation starts to decrease). The maximum flow rate at this substation is achieved at the 110th min. Although the supply temperature decreases till the 116th min, the further flow rate increase through SST A is not possible since the flow control valve is already fully open. Substantial changes of water flow rates in other two substations are not observed. Fig. 11 shows the measured and calculated temperature fronts at substations SST A, SST B and SST C in case of heat plant load reduction. Again, the good agreement is obtained, although the total supply flow rate from the heat plant is kept constant throughout the whole transient (the value of 0.11 m3/s is adopted according to Fig. 9), and the flow rates through the substations are kept constant (the values calibrated to the initial steady-state condition of the district heating system). Also in this case the shapes and amplitudes of the calculated temperature fronts are perfectly preserved. 4. Conclusions This paper presents the model and numerical method for the computer prediction of the thermal transients in the district heat-
References [1] Kunz VJ, Haldi PA, Sarlos G. Dynamic behaviour of district heating networks. Fernwaerme Int 1991;20:104–19. [2] Larsen HV, Palsson H, Bohm B, Ravn HF. Aggregated dynamic simulation model of district heating networks. Energy Convers Manage 2002;43:995–1019. [3] Larsen HV, Bohm B, Wigbels M. A comparison of aggregated models for simulation and operational optimization of district heating networks. Energy Convers Manage 2004;45:1119–39. [4] Benonysson A, Bøhm B, Ravn HF. Operational optimization in a district heating system. Energy Convers Manage 1995;36:297–314. [5] Park YS, Kim WT, Kim BK, Ha SK. Simple models for operational optimization, Chapter 2: state of the art. In: Bøhm B, editor. IEA district heating and cooling, annex VI, Report 2002: S1. Technical University of Denmark, Kgs. Lyngby; 2002, p. 7–10. [6] Gabrielaitiene I, Bøhm B, Sunden B. Modelling temperature dynamics of a district heating system in Naestved, Denmark – a case study. Energy Convers Manage 2007;48:78–86. [7] Iversen S, Ougaard P, Loppenthien JK. Dynamic temperature optimizationproviding instant results. EuroHeat 2006;3-II:46–9. [8] Karamarkovic, V, Solujic A, Simovic, J, Ramic B. Implementation program of energy development strategy in the field of industrial energy. In: The Society of Thermal Engineers of Serbia edition. Proceedings of the regional conference industrial energy and environmental protection in Southeast Europe, Zlatibor, Serbia, 24–28 June, 2008. [9] Stevanovic V, Jovanovic Z. A hybrid method for the numerical prediction of enthalpy transport in fluid flow. Int Commun Heat Mass Transfer 2000;27:23–34. [10] Demidovitch B, Maron I. Elements de Calcul Numerique. Moscou: Editions Mir; 1987. [11] Tannehill JC, Anderson DA, Pletcher RH. Computational fluid mechanics and heat transfer. Bristol: Taylor & Francis; 1997. [12] Stevanovic V, Prica S, Maslovaric B, Zivkovic B, Nikodijevic S. Efficient numerical method for district heating system hydraulics. Energy Convers Manage 2007;48:1536–43.