- Email: [email protected]
Function method for dynamic temperature simulation of district heating network
Applied Thermal Engineering 123 (2017) 682–688
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Function method for dynamic temperature simulation of district heating network Jinfu Zheng, Zhigang Zhou, Jianing Zhao ⇑, Jinda Wang School of Municipal and Environmental Engineering, Harbin Institute of Technology, Harbin 150090, China
h i g h l i g h t s A new physical method for dynamic temperature simulation of DHN was proposed. Dynamic thermal conditions can be solved analytically and efficiently. Method was validated at large and quick temperature changes and different flows. Simulations of new method is agreed better with measurements than node method. Compared with node method, the calculation time of new method is reduced by 37%.
a r t i c l e
i n f o
Article history: Received 8 February 2017 Revised 26 April 2017 Accepted 15 May 2017 Available online 24 May 2017 Keywords: District heating network Dynamic temperature simulation Function method Delay time Relative attenuation degree Calculation time
a b s t r a c t A thermal dynamic characteristic of district heating network was analyzed with an emphasis on simulating dynamic temperature distribution. Therefore, a new physical method (namely function method) and corresponding computer codes for dynamic temperature simulation were proposed. The function method takes into account various factors, such as flow time, heat capacity of the pipe and heat loss, in one single step by using Fourier series expansion as the mathematical basis to obtain the analytical solution of transient energy equation. Further, two key parameters, delay time and relative attenuation degree, were described to represent the time delay and heat loss respectively and were used to obtain the temperature distribution of the network. Moreover, the function method is validated by applying it to a real district heating network at different temperature change stages. For comparison purpose, the node method was also used, and the simulated supply temperatures at substations were compared with measurements. Comparisons of both methods were also performed considering different fluid velocities at heat source using numerical results. The results showed that the function method could be recommended to simulate the dynamic temperature of district heating network accurately and quickly for efficient system performance. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction District heating (DH) is represented by heat networks used to deliver heat from heat source (HS) to consumers to satisfy space heating [1]. Nowadays, heat produced by various sources in a DH system is widely used for energy efficiency and the heat production proportion of these sources vary over time depending on demand, price and availability, which results in heat redistribution. Additionally, renewable energy sources, such as wind, solar and geothermal energy, are gradually used to supply heat for con-
⇑ Corresponding author at: Box 2644, 73 Huanghe Road, Nangang District, Harbin Institute of Technology, Harbin 150090, China. E-mail address: [email protected] (J. Zhao). http://dx.doi.org/10.1016/j.applthermaleng.2017.05.083 1359-4311/Ó 2017 Elsevier Ltd. All rights reserved.
sumers of DH system, which changes irregularly due to weather conditions [2–4]. Furthermore, heat load variation for consumers due to changes in weather condition and daily demand also intensifies the process of heat redistribution [5,6]. As a result of abovementioned factors, operation conditions are never stable and static temperature simulation for a network [7] cannot satisfy dynamic operation conditions. Therefore, modeling a dynamic characteristic of a DH system becomes necessary for efficient system performance, such as determining optimal start-up time of heat plants or pumps to reduce pumping cost [8,9] and determining optimal supply temperature at HS to reduce heat loss cost [10,11]. In the published papers, there are mainly two models, statistical model and physical model, to simulate the dynamic temperature distribution of DH network. The statistical models, based on
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
683
Nomenclature A c cp C ext C int F Fp L M N n n0 tðs; 0Þ tðs; xÞ tðs; xÞ te t p ðs; xÞ
amplitude (m) specific heat of water (J/(kg °C)) specific heat of pipe wall (J/(kg °C)) exterior perimeter of pipe wall (m) internal perimeter of pipe wall (m) cross-sectional area of flow channel (m2) cross-sectional area of pipe wall (m2) pipe length (m) the maximum value of the second derivative of the function the number of Fourier series expansion number of collected data n0 ¼ 0:5n pipe inlet temperature (°C) temperature distribution of water (°C) average temperature of water (°C) environment temperature(°C) temperature distribution of pipe wall(°C)
standard transfer function models or neural networks, are generally simpler to update and easier to operate for state estimation of DH networks when compared to physical models but require availability of measurements for model estimation and validation for different operation strategies [12,13]. In addition, the estimation and tracking of time delays in the networks played an important role in statistical models and were reported to be the main source of problems, especially if temperatures were changed abruptly [14]. As for the physical models, a full physical modeling of the network is involved without introducing artificial changes to the structure of the network. The physical models mainly contain the element method firstly presented in [15], node method firstly presented in [15] and characteristic method [16]. Compared with node method, both the element method and characteristic method divide the pipe into many discrete ‘‘elements” or ‘‘nodes” and in each calculation step, every ‘‘element” or ‘‘node” is required to be calculated, which tend to be computationally intensive. Additionally, according to published literatures [14,15], the element method was found to be inferior to node method, both with respect to accuracy and stability. Node method is most widely used for modeling thermal dynamic characteristics of DH networks and the accuracy of node method has been validated in several real DH networks by comparing simulations and measurements of supply temperature at substations [5,6,17]. In node method, the temperature of outlet node needs three steps for establishing: firstly, estimating preliminarily the temperature of inlet node by taking flow time into account; secondly, taking into account the heat capacity of the pipe; finally, taking heat loss into account. It is not simple enough to obtain the outlet temperature based on the inlet temperature and will result in more time consumption for simulation. However, there have been few published literatures to optimize the node method or introduce a new method to reduce the time consumption for simulation with significant accuracy in recent years, which play an important role in online operating system and efficient system performance. Additionally, the accuracy validation of the dynamic methods mentioned above was not under condition of large and quick temperature changes occurring in the HS, which usually results in extreme deviations between measured and predicted temperatures. Furthermore, the validation methods were carried out at a fixed velocity only, not related to the validations at different fluid
u x
a b
g h
r r2E u x q qp s
water velocity (m/s) distance from the pipe inlet (m) convective heat transfer coefficient between flow and pipe inner wall (W/(m2 °C)) convective heat transfer coefficient between environment and pipe outer wall (W/(m °C)) delay time (h) initial phase (rad) error caused by substitution for inlet temperature mean squared error of temperature test data relative attenuation degree angular frequency (rad/h) density of water (kg/m3) density of pipe wall (kg/m3) time (h)
Subscripts i item of Fourier series expansion
velocities. Therefore, the aim of this paper was to solve those problems and make sufficient contributions to the knowledge. 2. Dynamic model of heating network In the process of heat transmission and distribution from HS to consumers, there are time delay and heat loss, which have a great influence on temperature prediction [18,19]. In addition, heating network consists of numerous pipes in which time delay and heat losses occur. Consequently, research on temperature prediction is mainly focused on time delay and heat loss of pipe [19]. 2.1. Modeling assumptions Generally, in order to investigate the main thermal dynamic characteristic of DH network, the following assumptions are used in physical models [20]: 1. The steady-state hydraulic condition is considered. Even though the volume flows vary over the simulation time, the pseudodynamic method is used, which means that the simulation time will be divided into many time steps and a new flow situation is determined by steady-state flow calculation at every time step. 2. The heat storage of insulation, mantle and ground outside the pipe is ignored. As for the insulation layer, the thermal conductivity should not be greater than 0.12 W/(m K), and the density should not be greater than 200 kg/m3 [21]. As a result, the heat storage of insulation is small and the effect of fluctuation of ground temperature on water temperature is small, so they can be ignored. 3. The fluid is treated as an ideal one. 4. The impact of hydraulic dispersion, thermal diffusion and axial heat transmission are neglected. These assumptions were used in most literatures about dynamic temperature simulation for their rationality [5,6,16,17]. 2.2. Modeling of district heating network The heat is transferred from the fluid to the pipe wall. Simultaneously, the heat is transferred from pipe wall to environment
684
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
after thermal conduction of insulation layer and soil, as well as thermal convection between ground and environment (see Fig. 1). Both the heat transfers result in temperature variation of the fluid and pipe wall. So the energy balance equations of the pipe wall and fluid are written as Eqs. (1) and (2).
@t p ðs; xÞ ¼ Ha ½tðs; xÞ tp ðs; xÞ þ Hb ½te tp ðs; xÞ @s
ð1Þ
@tðs; xÞ @tðs; xÞ þu ¼ H½t p ðs; xÞ tðs; xÞ @s @x
ð2Þ
where H ¼ aC int =ðFcqÞ, Ha ¼ aC int =ðF p cp qp Þ. In order to solve Eqs. (1) and (2) analytically, we need to know the form of analytical solution. After the heat is transmitted along a pipe, there is rightward shift (caused by time delay) and downward shift (caused by heat loss) from inlet temperature curve to outlet temperature curve (see Fig. 4), which can be illustrated by trigonometric functions: phase variation represents the rightward shift and amplitude variation represents the downward shift. Additionally, the pipe inlet temperature tðs; 0Þ satisfy the Dirichlet conditions [22] in general, which means that tðs; 0Þ can be expanded into Fourier series, as described in Eq. (3) [23,24].
tðs; 0Þ ¼
N X Ai sinðxi s þ hi Þ þ tðs; 0Þ
ð3Þ
i¼1
Based on analysis of the previous paragraph, the outlet temperature tðs; xÞ also varies in the form of sum of simple sine waves. The difference between the sine waves is that the phase of tðs; xÞ decreases and amplitude of tðs; xÞ decreases in the form of exponent relatively, and both of the decreases should be the function of length from entrance [24,25]. So tðs; xÞ can be expressed as follows:
tðs; xÞ ¼
N X Ai expðmi xÞ sinðxi s þ ni x þ hi Þ þ tðs; xÞ
ð4Þ
i¼1
Note that Eq. (1) contains environment temperature t e ; but it is not considered in the solution form described in Eq. (4). Mathematically, when Eq. (4) is substituted into Eq. (1) and Eq. (2), the uncertain parameters m, n cannot be solved. Physically, environment temperature has an effect on the heat loss to some extent [26]. In order to make the model solvable with accuracy, Eq. (4) is equivalently transformed into Eq. (5).
tðs; xÞ ¼
" # N X Ai expðmi xÞ sinðxi s þ ni x þ hi Þ þ te i¼1
þ ½tðs; xÞ t e
ð5Þ PN
Eq. (5) contains unsteady-state term i¼1 Ai expðmi xÞ sinðxi s þ ni x þ hi Þ þ t e (denoted by ~tðs; xÞ in the next contents) and steady-state term tðs; xÞ t e : In the two terms, only ~tðs; xÞ has an impact on the dynamic parameters of DH network, so
~tðs; xÞ is substituted into Eqs. (1) and (2). The temperature distribution of ~tðs; xÞ which is x meters from the inlet of pipe is obtained, as given by Eq. (6):
~tðs; xÞ ¼
" !# N X x Ha ðHa þ Hb Þ 1 Ai exp H u x2i þ ðHa þ Hb Þ2 i¼1 " !# x HHa þ te sin xi s þ hi xi 1þ u x2i þ ðHa þ Hb Þ2
Based on Eqs. (3) and (6), two parameters, delay time and relative attenuation degree, were defined to represent the thermal dynamic characteristics of the DH systems. The delay time (g) refers to the time it takes for the water temperature to change from the inlet temperature to the outlet temperature. The relative attenuation degree (u) refers to the radio of outlet temperature to the inlet temperature. In Eqs. (3) and (6), there are summation symbols and one-toone correspondence for each item between Eqs. (3) and (6), which means that Ai expðmi xÞ sinðxi s þ ni x þ hi Þ is only obtained by Ai sinðxi s þ hi Þ. As a result, each item has one delay time (gi ) and one relative attenuation degree (ui ), so we average them to obtain practical delay time (g) and relative attenuation degree (u) representing the thermal dynamic characteristics of the DH systems, as described in Eqs. (7) and (8). N 1X L Ha H b 1þ g¼ 2 N i¼1 u x þ Ha þ Hb 2
! ð7Þ
i
" !# N 1X L Ha ðHa þ Hb Þ 1 u¼ exp H N i¼1 u x2i þ ðHa þ Hb Þ2
ð8Þ
In Eqs. (7) and (8), the number of Fourier series expansion N is a major parameter, so the minimum error-analysis method is used to derive the parameter N. The inlet temperature is decomposed into the sum of a set of sine functions (Fourier series) based on measured temperature. This substitution will result in errors in the following aspects: (1) the continuous temperature field is replaced with discrete field; (2) the Fourier series is truncated into finite terms; (3) the measured data contains system errors which cannot be eliminated artificially. From aspects (1) and (2) it can be found that more is the number of Fourier series expansions considered, more accurate the substitution for inlet temperature is, but the aspect (3) causes it to be inversed. Therefore, it does not mean that more N is better and there exists an optimal number N to minimize the error [27,28]. The objective function used to minimize the error is given in Eq. (9) and derivation of the objective function can be obtained from literatures [27,29].
" min rðNÞ ¼ n 4M
2
1 1 þ n0 N
2
# þ 2Nr
2 E
ð9Þ
Through the studies mentioned in previous contents, the delay time (g) and relative attenuation degree (u) were obtained, taking into consideration various factors: water temperature, fluid velocity, pipe diameter, overall heat transfer coefficient, pipe length, outer temperature, number of Fourier series expansion and so on. Therefore, the pipe outlet temperature can be obtained in one single step:
t½ðs þ gÞ; L ¼ utðs; 0Þ þ t e ðsÞð1 uÞ Fig. 1. Schematic diagram of directly buried pipe.
ð6Þ
ð10Þ
Therefore, based on the Eqs. (1)–(10), the procedure of the calculation can be clarified as follows:
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
685
Fig. 4. The measured temperature at HS and substations.
Fig. 2. The district heating system in Changchun.
1. Firstly, the measured supply temperature at HS is used to obtain the Fourier series of inlet temperature of pipes and the optimal number (N) by Eq.(9). 2. Secondly, N and related constitutive parameters are substituted into Eqs. (7) and (8) to obtain parameter g and u for every pipe. 3. Thirdly, substituting g, u and known inlet temperature (usually, the known inlet temperature is the supply temperature at the HS) into Eq. (10) to obtain the pipe outlet temperature, which is also the inlet temperature of connection pipes [6]. Therefore, the temperature distribution of the entire network can be obtained.
Fig. 5. The measured velocity and environment temperature at HS.
of accuracy and calculation time for simulation. The procedure of node method can be easily found in the published papers [14,17]. 3.1. DH system
3. Model validations The function method was applied to simulate dynamic supply temperature in substations at three different temperature change stages in a real district heating network in Changchun, a northeastern city in China. Simultaneously, the node method was also chosen for providing a comparison to the function method in terms
In the DH system, the installed power of heat plant is about 406 MW and serves more than 2.9 million square meters of building area. Fig. 2 shows the distribution network of the DH system. The network is constructed of pre-insulated pipes with total length of approximately 6.6 km. The pipe diameters range from 80 mm to 1000 mm.
Fig. 3. Schematic view of measurement site.
686
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
Fig. 6. The measured and simulated supply temperature at SU1.
Fig. 8. The measured and simulated supply temperature at SU3.
Fig. 7. The measured and simulated supply temperature at SU2.
Fig. 9. The measured and simulated supply temperature at SU4.
3.2. Operational strategy In entire heating season, the DH system is operated by a simple strategy, ‘‘change the flow mass at different stages” [30], which means that the supply temperature is regulated by environment temperature, besides, the circulation flow remains constant over a period of time, which could be one or several months. Additionally, the primary side at substations is not equipped with automatic control device, and do not consume heat for sanitary hot water consumption, which means that when the circulation flow at HS remains constant, the steady-state hydraulic condition is obtained. Therefore, the operational strategy will result in a substantial temperature variation with steady-state hydraulic condition, which is advantageous in studying the thermal dynamic characteristics of DH network (see Figs. 4 and 5). Fig. 10. The measured and simulated supply temperature at SU5.
3.3. Measurements Measurements and simulations of supply temperature are presented for substations located at different distances from HS (see Fig. 2): SU1–SU5 (distances 1121 m, 1881 m, 3499 m, 4566 m and 6581 m from HS, respectively). SU5 is the most remote substation in the DH system with most branches, turns and junctions. Temperature and flow measurements were performed at supply pipe for HS and substations (see Fig. 3). A thermoelectric
thermometer (with accuracy of 0.2 °C) and an ultrasonic flow meter (with accuracy of 0.1 t/h) were used to measure temperature and flow and therefore obtain them precisely. The measurements were carried out every 10 min on 15th January 2015, which is during the coldest period in winter 2015. The measured supply temperature at HS and substations are shown in Fig. 4. The measured velocity and environment temperature at HS are shown in Fig. 5.
687
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688 Table 1 Average error and standard deviation for predicted supply temperature at substations for time period of 6:00–21:00. Methods
SU1 Aver. error (Std. dev.) [°C]
SU2 Aver. error (Std. dev.) [°C]
SU3 Aver. error (Std. dev.) [°C]
SU4 Aver. error (Std. dev.) [°C]
SU5 Aver. error (Std. dev.) [°C]
Average value Aver. error (Std. dev.) [°C]
Function method Node method
0.14 (0.46) 0.23 (0.48)
0.14 (0.96) 0.31 (0.99)
0.23 (1.43) 0.41 (1.48)
0.16 (1.66) 0.49 (1.71)
0.41 (1.74) 0.67 (1.79)
0.22 (1.33) 0.42 (1.38)
Table 2 Average error and standard deviation for predicted supply temperature at substations for different temperature change stages. Methods
Stages
SU1 Aver. error (Std. dev.) [°C]
SU2 Aver. error (Std. dev.) [°C]
SU3 Aver. error (Std. dev.) [°C]
SU4 Aver. error (Std. dev.) [°C]
SU5 Aver. error (Std. dev.) [°C]
Average value Aver. error (Std. dev.) [°C]
Function method
Stage 1 Stage 2 Stage 3
0.13 (0.49) 0.01 (0.69) 0.25 (0.29)
0.73 (1.23) 0.75 (1.52) 0.17 (0.28)
1.26 (1.59) 0.88 (1.56) 0.31 (0.37)
1.45 (1.73) 1.68 (2.63) 0.36 (0.37)
2.09 (2.18) 1.91 (2.97) 0.55 (0.58)
1.13 (1.55) 1.05 (2.05) 0.33 (0.39)
Node method
Stage 1 Stage 2 Stage 3
0.18 (0.52) 0.11 (0.66) 0.36 (0.36)
0.88 (1.25) 0.58 (1.43) 0.35 (0.34)
1.32 (2.00) 0.60 (1.43) 0.52 (0.57)
1.64 (2.25) 1.25 (2.31) 0.73 (0.71)
2.46 (3.10) 1.48 (2.69) 1.02 (0.98)
1.30 (2.03) 0.76 (1.85) 0.60 (0.64)
3.4. Results and discussion In order to highlight the difference at stages of large and quick temperature changes, the results of the supply temperature at substations obtained by both methods are presented for a time period of 6:00–21:00 together with the measurements, as shown in Figs. 6–10. As seen in Figs. 6–10, the results obtained by both models agrees well with the measured data and the difference between both models was small. Therefore, average error (aver. error) and standard deviation (std. dev) were introduced to compare the accuracy of the models for the time period of 6:00–21:00, as presented in Table 1. Table 1 discloses that average error and standard deviation of function method were less than that of node method. Therefore, we can conclude that function method is more accurate than node method for simulating the supply temperature at substations in this case. In addition, the delay times obtained by both models are almost same, which indicates that substituting flow time of water (used in node method) for temperature wave spread time (used in function method) will not result in any significant deviation of delay time. As also seen in Figs. 6–10, the magnitude of deviation between the measured and predicted temperatures is not uniform during the simulation period. Therefore, comparisons were also carried out at different temperature change stages, namely rapid dropping stage (stage 1), rapid rising stage (stage 2) and relatively stable stage (stage 3), as presented in Table 2 for a period of 3 h, which is about the time taken for temperature to reach minimum at rapid-dropping stage, the time taken for temperature to reach maximum at rapid-rising stage or the time for relatively stable stage after rapid-rising stage. From Table 2, it can be seen that average error and standard deviation of both methods at the three stages tend to be larger with the increase in distance from HS. The same phenomenon can be seen in Figs. 6–10 and Table 1. A possible reason is that the effects of branches, turns, junctions and other fittings are not taken into account in the present models [5]. Simultaneously, the average error and standard deviation of both methods at rapid dropping stage and rapid rising stage are greater than that at relatively stable stage. This is mainly affected by the delay time: a little bit of delay time deviation has a great
Fig. 11. Simulated supply temperature of SU4 at different fluid velocities.
effect on stage of large and quick temperature changes. The accuracy at relatively stable stage is mainly affected by the relative attenuation degree and it is not as intense as the effect of delay time on stage of large temperature changes. As also seen in Table 2, the function method is more accurate than the node method at the rapid dropping stage and relatively stable stage, while it is opposite at rapid rising stage. This phenomenon is caused by several factors, such as delay time, relative attenuation degree and the trend of temperature changes. Taking Fig. 10 as an example, on the basis of overestimation of delay time and relative attenuation degree, at rapid dropping stage, average error and standard deviation increase with the increase of relative attenuation degree; on the contrary, average error and standard deviation decrease with the increase of relative attenuation degree at rapid rising stage. Other situations require a specific analysis. Additionally, the calculation time of function method for simulation was also compared with node method. The result indicates that both methods can run on a computer with normal configuration, but the function method is 37% faster than the node method on the same computer configuration. The optimization of calculation time is mainly due to reduction of the calculation steps (‘‘three steps” used in node method) by taking account of various factors in one single step (Eq. (10)).
688
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
References
Fig. 12. Simulated supply temperature of SU5 at different fluid velocities.
Comparisons were also performed considering different fluid velocities at HS using numerical results. The simulated supply temperature at SU4–SU5 are shown in Figs. 11 and 12, respectively. As seen in Fig. 12, higher the fluid velocity is, faster the temperature wave propagates and less the heat loss is. This is because the increased heat loss caused by strengthened thermal convection between water and pipe is much less than the reduced heat loss due to reduction of the propagation time. It can also be seen that irrespective of higher or lower velocity, the simulation results of both methods are similar. Therefore, the analysis at fluid velocity of 1.17 m/s is also applicable for different fluid velocities and we can conclude that the proposed method can also perform well when volume flows in the network are varying by using pseudodynamic method. 4. Conclusions The main conclusions of this study were summarized as follows: 1. For function method, the average error and standard deviation over the simulation time are 0.22 °C and 0.42 °C respectively, both of which are smaller than that of node method. 2. The function method is more accurate than the node method at the rapid dropping stage and relatively stable stage, while it is opposite at rapid rising stage. 3. Compared with node method, the calculation time of function method is reduced by approximately 37%. 4. The proposed method can also perform well when volume flows in the network are varying. Therefore, the function method could be recommended to simulate the dynamic temperature of district heating network accurately and efficiently. Acknowledgements This work was supported by the National Key Technologies P&D Program of China (No. 2015BAA01B01) and the Key Projects of Natural Science Foundation of Heilongjiang Province of China (No. ZD2016010).
[1] G. Comodi, M. Lorenzetti, D. Salvi, A. Arteconi, Criticalities of district heating in Southern Europe: lesson learned from a CHP-DH in Central Italy, Appl. Therm. Eng. 112 (2016) 649–659. [2] Z. Pan, Q. Guo, H. Sun, Feasible region method based integrated heat and electricity dispatch considering building thermal inertia, Appl. Energy (2016). [3] I.B. Hassine, U. Eicker, Impact of load structure variation and solar thermal energy integration on an existing district heating network, Appl. Therm. Eng. 50 (2013) 1437–1446. [4] S.A. Kyriakis, P.L. Younger, Towards the increased utilization of geothermal energy in a district heating network through the use of a heat storage, Appl. Therm. Eng. 94 (2015) 99–110. [5] I. Gabrielaitiene, B. Bøhm, B. Sunden, Modelling temperature dynamics of a district heating system in Naestved, Denmark—a case study, Energy Convers. Manage. 48 (2007) 78–86. [6] I. Gabrielaitiene, B. Bøhm, B. Sunden, Dynamic temperature simulation in district heating systems in Denmark regarding pronounced transient behaviour, J. Civil Eng. Manage. 17 (2011) 79–87. [7] M. Kuosa, K. Kontu, T. Mäkilä, M. Lampinen, R. Lahdelma, Static study of traditional and ring networks and the use of mass flow control in district heating applications, Appl. Therm. Eng. 54 (2013) 450–459. [8] X. Sheng, D. Lin, Energy saving analyses on the reconstruction project in district heating system with distributed variable speed pumps, Appl. Therm. Eng. 101 (2016) 432–445. [9] P. Jie, N. Zhu, D. Li, Operation optimization of existing district heating systems, Appl. Therm. Eng. 78 (2015) 278–288. [10] A. Colmenar-Santos, E. Rosales-Asensio, D. Borge-Diez, E. Collado-Fernández, Evaluation of the cost of using power plant reject heat in low-temperature district heating and cooling networks, Appl. Energy 162 (2016) 892–907. [11] S. Lim, S. Park, H. Chung, M. Kim, Y.J. Baik, S. Shin, Dynamic modeling of building heat network system using Simulink, Appl. Therm. Eng. 84 (2015) 375–389. [12] O.P. Palsson, Stochastic modeling, control and optimization of district heating systems, 1993. [13] H. Zhao, Analysis, Modelling and Operational Optimization of District Heating Systems, Technical University of Denmark, 1995. [14] H. Palsson, Analysis of numerical methods for simulating temperature dynamics in district heating pipes, in: International Symposium on District Heating and Cooling Simulation, 2010, pp. 219–224. [15] A. Benonysson, Dynamic modelling and operational optimization of district heating systems, Laboratory of Heating and Air Conditioning, Technical University of Denmark, 1991. [16] V.D. Stevanovic, B. Zivkovic, S. Prica, B. Maslovaric, V. Karamarkovic, V. Trkulja, Prediction of thermal transients in district heating systems, Energy Convers. Manage.ent 50 (2009) 2167–2173. [17] I. Gabrielaitiene, B. Bøhm, B. Sunden, Evaluation of approaches for modeling temperature wave propagation in district heating pipelines, Heat Transfer Eng. 29 (2008) 45–56. [18] H. Zhao, Optimum operation of a CHP-based district heating system by mathematical modelling, Fuel Energy Abstr. (1996) 296. [19] L. Dobos, J. Jäschke, J. Abonyi, S. Skogestad, Dynamic model and control of heat exchanger networks for district heating, Hung. J. Ind. Chem. 37 (2009) 37–49. [20] H.V. Larsen, M. Lucht, M. Wigbels, Simple models for operational optimization, in: IEA District heating and cooling, Technical University of Denmark, 2002. [21] MOHURD, Design code for city heating network, in: Vol. CJJ34 - 2010, China Architecture & Building Press, Beijing, 2010. [22] Z. Zhang, Engineering Mathematics Analysis, China Higher Education Press, 2008. [23] Q. Bing, Studies on the dynamic thermal characteristics of district heating systems, Tsinghua University, 2004, pp. 162. [24] P. Jie, Z. Tian, S. Yuan, N. Zhu, Modeling the dynamic characteristics of a district heating network, Energy 39 (2012) 126–134. [25] J. Duquette, A. Rowe, P. Wild, Thermal performance of a steady state physical pipe model for simulating district heating grids with variable flow, Appl. Energy 178 (2016) 383–393. [26] K. Çomaklı, B. Yüksel, Ö. Çomaklı, Evaluation of energy and exergy losses in district heating network, Appl. Therm. Eng. 24 (2004) 1009–1017. [27] L. Hua, Y. Wang, Numerical Integration and Application, Science Press, Beijing, 1963. [28] J.A. Simmons, Dynamical Prediction: Some Results from Operational Forecasting and Research Experiments, World Meteorological Organization, 1983. [29] A. Xiu, D. Liao, Study on the determination of optimal Fourier series, Acta Meteorol. Sin. 49 (1991) 308–320. [30] P. He, G. Sun, F. Wang, Heating Engineering, fourth ed., China Architecture & Building Press, Beijing, 2009.
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Research Paper
Function method for dynamic temperature simulation of district heating network Jinfu Zheng, Zhigang Zhou, Jianing Zhao ⇑, Jinda Wang School of Municipal and Environmental Engineering, Harbin Institute of Technology, Harbin 150090, China
h i g h l i g h t s A new physical method for dynamic temperature simulation of DHN was proposed. Dynamic thermal conditions can be solved analytically and efficiently. Method was validated at large and quick temperature changes and different flows. Simulations of new method is agreed better with measurements than node method. Compared with node method, the calculation time of new method is reduced by 37%.
a r t i c l e
i n f o
Article history: Received 8 February 2017 Revised 26 April 2017 Accepted 15 May 2017 Available online 24 May 2017 Keywords: District heating network Dynamic temperature simulation Function method Delay time Relative attenuation degree Calculation time
a b s t r a c t A thermal dynamic characteristic of district heating network was analyzed with an emphasis on simulating dynamic temperature distribution. Therefore, a new physical method (namely function method) and corresponding computer codes for dynamic temperature simulation were proposed. The function method takes into account various factors, such as flow time, heat capacity of the pipe and heat loss, in one single step by using Fourier series expansion as the mathematical basis to obtain the analytical solution of transient energy equation. Further, two key parameters, delay time and relative attenuation degree, were described to represent the time delay and heat loss respectively and were used to obtain the temperature distribution of the network. Moreover, the function method is validated by applying it to a real district heating network at different temperature change stages. For comparison purpose, the node method was also used, and the simulated supply temperatures at substations were compared with measurements. Comparisons of both methods were also performed considering different fluid velocities at heat source using numerical results. The results showed that the function method could be recommended to simulate the dynamic temperature of district heating network accurately and quickly for efficient system performance. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction District heating (DH) is represented by heat networks used to deliver heat from heat source (HS) to consumers to satisfy space heating [1]. Nowadays, heat produced by various sources in a DH system is widely used for energy efficiency and the heat production proportion of these sources vary over time depending on demand, price and availability, which results in heat redistribution. Additionally, renewable energy sources, such as wind, solar and geothermal energy, are gradually used to supply heat for con-
⇑ Corresponding author at: Box 2644, 73 Huanghe Road, Nangang District, Harbin Institute of Technology, Harbin 150090, China. E-mail address: [email protected] (J. Zhao). http://dx.doi.org/10.1016/j.applthermaleng.2017.05.083 1359-4311/Ó 2017 Elsevier Ltd. All rights reserved.
sumers of DH system, which changes irregularly due to weather conditions [2–4]. Furthermore, heat load variation for consumers due to changes in weather condition and daily demand also intensifies the process of heat redistribution [5,6]. As a result of abovementioned factors, operation conditions are never stable and static temperature simulation for a network [7] cannot satisfy dynamic operation conditions. Therefore, modeling a dynamic characteristic of a DH system becomes necessary for efficient system performance, such as determining optimal start-up time of heat plants or pumps to reduce pumping cost [8,9] and determining optimal supply temperature at HS to reduce heat loss cost [10,11]. In the published papers, there are mainly two models, statistical model and physical model, to simulate the dynamic temperature distribution of DH network. The statistical models, based on
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
683
Nomenclature A c cp C ext C int F Fp L M N n n0 tðs; 0Þ tðs; xÞ tðs; xÞ te t p ðs; xÞ
amplitude (m) specific heat of water (J/(kg °C)) specific heat of pipe wall (J/(kg °C)) exterior perimeter of pipe wall (m) internal perimeter of pipe wall (m) cross-sectional area of flow channel (m2) cross-sectional area of pipe wall (m2) pipe length (m) the maximum value of the second derivative of the function the number of Fourier series expansion number of collected data n0 ¼ 0:5n pipe inlet temperature (°C) temperature distribution of water (°C) average temperature of water (°C) environment temperature(°C) temperature distribution of pipe wall(°C)
standard transfer function models or neural networks, are generally simpler to update and easier to operate for state estimation of DH networks when compared to physical models but require availability of measurements for model estimation and validation for different operation strategies [12,13]. In addition, the estimation and tracking of time delays in the networks played an important role in statistical models and were reported to be the main source of problems, especially if temperatures were changed abruptly [14]. As for the physical models, a full physical modeling of the network is involved without introducing artificial changes to the structure of the network. The physical models mainly contain the element method firstly presented in [15], node method firstly presented in [15] and characteristic method [16]. Compared with node method, both the element method and characteristic method divide the pipe into many discrete ‘‘elements” or ‘‘nodes” and in each calculation step, every ‘‘element” or ‘‘node” is required to be calculated, which tend to be computationally intensive. Additionally, according to published literatures [14,15], the element method was found to be inferior to node method, both with respect to accuracy and stability. Node method is most widely used for modeling thermal dynamic characteristics of DH networks and the accuracy of node method has been validated in several real DH networks by comparing simulations and measurements of supply temperature at substations [5,6,17]. In node method, the temperature of outlet node needs three steps for establishing: firstly, estimating preliminarily the temperature of inlet node by taking flow time into account; secondly, taking into account the heat capacity of the pipe; finally, taking heat loss into account. It is not simple enough to obtain the outlet temperature based on the inlet temperature and will result in more time consumption for simulation. However, there have been few published literatures to optimize the node method or introduce a new method to reduce the time consumption for simulation with significant accuracy in recent years, which play an important role in online operating system and efficient system performance. Additionally, the accuracy validation of the dynamic methods mentioned above was not under condition of large and quick temperature changes occurring in the HS, which usually results in extreme deviations between measured and predicted temperatures. Furthermore, the validation methods were carried out at a fixed velocity only, not related to the validations at different fluid
u x
a b
g h
r r2E u x q qp s
water velocity (m/s) distance from the pipe inlet (m) convective heat transfer coefficient between flow and pipe inner wall (W/(m2 °C)) convective heat transfer coefficient between environment and pipe outer wall (W/(m °C)) delay time (h) initial phase (rad) error caused by substitution for inlet temperature mean squared error of temperature test data relative attenuation degree angular frequency (rad/h) density of water (kg/m3) density of pipe wall (kg/m3) time (h)
Subscripts i item of Fourier series expansion
velocities. Therefore, the aim of this paper was to solve those problems and make sufficient contributions to the knowledge. 2. Dynamic model of heating network In the process of heat transmission and distribution from HS to consumers, there are time delay and heat loss, which have a great influence on temperature prediction [18,19]. In addition, heating network consists of numerous pipes in which time delay and heat losses occur. Consequently, research on temperature prediction is mainly focused on time delay and heat loss of pipe [19]. 2.1. Modeling assumptions Generally, in order to investigate the main thermal dynamic characteristic of DH network, the following assumptions are used in physical models [20]: 1. The steady-state hydraulic condition is considered. Even though the volume flows vary over the simulation time, the pseudodynamic method is used, which means that the simulation time will be divided into many time steps and a new flow situation is determined by steady-state flow calculation at every time step. 2. The heat storage of insulation, mantle and ground outside the pipe is ignored. As for the insulation layer, the thermal conductivity should not be greater than 0.12 W/(m K), and the density should not be greater than 200 kg/m3 [21]. As a result, the heat storage of insulation is small and the effect of fluctuation of ground temperature on water temperature is small, so they can be ignored. 3. The fluid is treated as an ideal one. 4. The impact of hydraulic dispersion, thermal diffusion and axial heat transmission are neglected. These assumptions were used in most literatures about dynamic temperature simulation for their rationality [5,6,16,17]. 2.2. Modeling of district heating network The heat is transferred from the fluid to the pipe wall. Simultaneously, the heat is transferred from pipe wall to environment
684
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
after thermal conduction of insulation layer and soil, as well as thermal convection between ground and environment (see Fig. 1). Both the heat transfers result in temperature variation of the fluid and pipe wall. So the energy balance equations of the pipe wall and fluid are written as Eqs. (1) and (2).
@t p ðs; xÞ ¼ Ha ½tðs; xÞ tp ðs; xÞ þ Hb ½te tp ðs; xÞ @s
ð1Þ
@tðs; xÞ @tðs; xÞ þu ¼ H½t p ðs; xÞ tðs; xÞ @s @x
ð2Þ
where H ¼ aC int =ðFcqÞ, Ha ¼ aC int =ðF p cp qp Þ. In order to solve Eqs. (1) and (2) analytically, we need to know the form of analytical solution. After the heat is transmitted along a pipe, there is rightward shift (caused by time delay) and downward shift (caused by heat loss) from inlet temperature curve to outlet temperature curve (see Fig. 4), which can be illustrated by trigonometric functions: phase variation represents the rightward shift and amplitude variation represents the downward shift. Additionally, the pipe inlet temperature tðs; 0Þ satisfy the Dirichlet conditions [22] in general, which means that tðs; 0Þ can be expanded into Fourier series, as described in Eq. (3) [23,24].
tðs; 0Þ ¼
N X Ai sinðxi s þ hi Þ þ tðs; 0Þ
ð3Þ
i¼1
Based on analysis of the previous paragraph, the outlet temperature tðs; xÞ also varies in the form of sum of simple sine waves. The difference between the sine waves is that the phase of tðs; xÞ decreases and amplitude of tðs; xÞ decreases in the form of exponent relatively, and both of the decreases should be the function of length from entrance [24,25]. So tðs; xÞ can be expressed as follows:
tðs; xÞ ¼
N X Ai expðmi xÞ sinðxi s þ ni x þ hi Þ þ tðs; xÞ
ð4Þ
i¼1
Note that Eq. (1) contains environment temperature t e ; but it is not considered in the solution form described in Eq. (4). Mathematically, when Eq. (4) is substituted into Eq. (1) and Eq. (2), the uncertain parameters m, n cannot be solved. Physically, environment temperature has an effect on the heat loss to some extent [26]. In order to make the model solvable with accuracy, Eq. (4) is equivalently transformed into Eq. (5).
tðs; xÞ ¼
" # N X Ai expðmi xÞ sinðxi s þ ni x þ hi Þ þ te i¼1
þ ½tðs; xÞ t e
ð5Þ PN
Eq. (5) contains unsteady-state term i¼1 Ai expðmi xÞ sinðxi s þ ni x þ hi Þ þ t e (denoted by ~tðs; xÞ in the next contents) and steady-state term tðs; xÞ t e : In the two terms, only ~tðs; xÞ has an impact on the dynamic parameters of DH network, so
~tðs; xÞ is substituted into Eqs. (1) and (2). The temperature distribution of ~tðs; xÞ which is x meters from the inlet of pipe is obtained, as given by Eq. (6):
~tðs; xÞ ¼
" !# N X x Ha ðHa þ Hb Þ 1 Ai exp H u x2i þ ðHa þ Hb Þ2 i¼1 " !# x HHa þ te sin xi s þ hi xi 1þ u x2i þ ðHa þ Hb Þ2
Based on Eqs. (3) and (6), two parameters, delay time and relative attenuation degree, were defined to represent the thermal dynamic characteristics of the DH systems. The delay time (g) refers to the time it takes for the water temperature to change from the inlet temperature to the outlet temperature. The relative attenuation degree (u) refers to the radio of outlet temperature to the inlet temperature. In Eqs. (3) and (6), there are summation symbols and one-toone correspondence for each item between Eqs. (3) and (6), which means that Ai expðmi xÞ sinðxi s þ ni x þ hi Þ is only obtained by Ai sinðxi s þ hi Þ. As a result, each item has one delay time (gi ) and one relative attenuation degree (ui ), so we average them to obtain practical delay time (g) and relative attenuation degree (u) representing the thermal dynamic characteristics of the DH systems, as described in Eqs. (7) and (8). N 1X L Ha H b 1þ g¼ 2 N i¼1 u x þ Ha þ Hb 2
! ð7Þ
i
" !# N 1X L Ha ðHa þ Hb Þ 1 u¼ exp H N i¼1 u x2i þ ðHa þ Hb Þ2
ð8Þ
In Eqs. (7) and (8), the number of Fourier series expansion N is a major parameter, so the minimum error-analysis method is used to derive the parameter N. The inlet temperature is decomposed into the sum of a set of sine functions (Fourier series) based on measured temperature. This substitution will result in errors in the following aspects: (1) the continuous temperature field is replaced with discrete field; (2) the Fourier series is truncated into finite terms; (3) the measured data contains system errors which cannot be eliminated artificially. From aspects (1) and (2) it can be found that more is the number of Fourier series expansions considered, more accurate the substitution for inlet temperature is, but the aspect (3) causes it to be inversed. Therefore, it does not mean that more N is better and there exists an optimal number N to minimize the error [27,28]. The objective function used to minimize the error is given in Eq. (9) and derivation of the objective function can be obtained from literatures [27,29].
" min rðNÞ ¼ n 4M
2
1 1 þ n0 N
2
# þ 2Nr
2 E
ð9Þ
Through the studies mentioned in previous contents, the delay time (g) and relative attenuation degree (u) were obtained, taking into consideration various factors: water temperature, fluid velocity, pipe diameter, overall heat transfer coefficient, pipe length, outer temperature, number of Fourier series expansion and so on. Therefore, the pipe outlet temperature can be obtained in one single step:
t½ðs þ gÞ; L ¼ utðs; 0Þ þ t e ðsÞð1 uÞ Fig. 1. Schematic diagram of directly buried pipe.
ð6Þ
ð10Þ
Therefore, based on the Eqs. (1)–(10), the procedure of the calculation can be clarified as follows:
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
685
Fig. 4. The measured temperature at HS and substations.
Fig. 2. The district heating system in Changchun.
1. Firstly, the measured supply temperature at HS is used to obtain the Fourier series of inlet temperature of pipes and the optimal number (N) by Eq.(9). 2. Secondly, N and related constitutive parameters are substituted into Eqs. (7) and (8) to obtain parameter g and u for every pipe. 3. Thirdly, substituting g, u and known inlet temperature (usually, the known inlet temperature is the supply temperature at the HS) into Eq. (10) to obtain the pipe outlet temperature, which is also the inlet temperature of connection pipes [6]. Therefore, the temperature distribution of the entire network can be obtained.
Fig. 5. The measured velocity and environment temperature at HS.
of accuracy and calculation time for simulation. The procedure of node method can be easily found in the published papers [14,17]. 3.1. DH system
3. Model validations The function method was applied to simulate dynamic supply temperature in substations at three different temperature change stages in a real district heating network in Changchun, a northeastern city in China. Simultaneously, the node method was also chosen for providing a comparison to the function method in terms
In the DH system, the installed power of heat plant is about 406 MW and serves more than 2.9 million square meters of building area. Fig. 2 shows the distribution network of the DH system. The network is constructed of pre-insulated pipes with total length of approximately 6.6 km. The pipe diameters range from 80 mm to 1000 mm.
Fig. 3. Schematic view of measurement site.
686
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
Fig. 6. The measured and simulated supply temperature at SU1.
Fig. 8. The measured and simulated supply temperature at SU3.
Fig. 7. The measured and simulated supply temperature at SU2.
Fig. 9. The measured and simulated supply temperature at SU4.
3.2. Operational strategy In entire heating season, the DH system is operated by a simple strategy, ‘‘change the flow mass at different stages” [30], which means that the supply temperature is regulated by environment temperature, besides, the circulation flow remains constant over a period of time, which could be one or several months. Additionally, the primary side at substations is not equipped with automatic control device, and do not consume heat for sanitary hot water consumption, which means that when the circulation flow at HS remains constant, the steady-state hydraulic condition is obtained. Therefore, the operational strategy will result in a substantial temperature variation with steady-state hydraulic condition, which is advantageous in studying the thermal dynamic characteristics of DH network (see Figs. 4 and 5). Fig. 10. The measured and simulated supply temperature at SU5.
3.3. Measurements Measurements and simulations of supply temperature are presented for substations located at different distances from HS (see Fig. 2): SU1–SU5 (distances 1121 m, 1881 m, 3499 m, 4566 m and 6581 m from HS, respectively). SU5 is the most remote substation in the DH system with most branches, turns and junctions. Temperature and flow measurements were performed at supply pipe for HS and substations (see Fig. 3). A thermoelectric
thermometer (with accuracy of 0.2 °C) and an ultrasonic flow meter (with accuracy of 0.1 t/h) were used to measure temperature and flow and therefore obtain them precisely. The measurements were carried out every 10 min on 15th January 2015, which is during the coldest period in winter 2015. The measured supply temperature at HS and substations are shown in Fig. 4. The measured velocity and environment temperature at HS are shown in Fig. 5.
687
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688 Table 1 Average error and standard deviation for predicted supply temperature at substations for time period of 6:00–21:00. Methods
SU1 Aver. error (Std. dev.) [°C]
SU2 Aver. error (Std. dev.) [°C]
SU3 Aver. error (Std. dev.) [°C]
SU4 Aver. error (Std. dev.) [°C]
SU5 Aver. error (Std. dev.) [°C]
Average value Aver. error (Std. dev.) [°C]
Function method Node method
0.14 (0.46) 0.23 (0.48)
0.14 (0.96) 0.31 (0.99)
0.23 (1.43) 0.41 (1.48)
0.16 (1.66) 0.49 (1.71)
0.41 (1.74) 0.67 (1.79)
0.22 (1.33) 0.42 (1.38)
Table 2 Average error and standard deviation for predicted supply temperature at substations for different temperature change stages. Methods
Stages
SU1 Aver. error (Std. dev.) [°C]
SU2 Aver. error (Std. dev.) [°C]
SU3 Aver. error (Std. dev.) [°C]
SU4 Aver. error (Std. dev.) [°C]
SU5 Aver. error (Std. dev.) [°C]
Average value Aver. error (Std. dev.) [°C]
Function method
Stage 1 Stage 2 Stage 3
0.13 (0.49) 0.01 (0.69) 0.25 (0.29)
0.73 (1.23) 0.75 (1.52) 0.17 (0.28)
1.26 (1.59) 0.88 (1.56) 0.31 (0.37)
1.45 (1.73) 1.68 (2.63) 0.36 (0.37)
2.09 (2.18) 1.91 (2.97) 0.55 (0.58)
1.13 (1.55) 1.05 (2.05) 0.33 (0.39)
Node method
Stage 1 Stage 2 Stage 3
0.18 (0.52) 0.11 (0.66) 0.36 (0.36)
0.88 (1.25) 0.58 (1.43) 0.35 (0.34)
1.32 (2.00) 0.60 (1.43) 0.52 (0.57)
1.64 (2.25) 1.25 (2.31) 0.73 (0.71)
2.46 (3.10) 1.48 (2.69) 1.02 (0.98)
1.30 (2.03) 0.76 (1.85) 0.60 (0.64)
3.4. Results and discussion In order to highlight the difference at stages of large and quick temperature changes, the results of the supply temperature at substations obtained by both methods are presented for a time period of 6:00–21:00 together with the measurements, as shown in Figs. 6–10. As seen in Figs. 6–10, the results obtained by both models agrees well with the measured data and the difference between both models was small. Therefore, average error (aver. error) and standard deviation (std. dev) were introduced to compare the accuracy of the models for the time period of 6:00–21:00, as presented in Table 1. Table 1 discloses that average error and standard deviation of function method were less than that of node method. Therefore, we can conclude that function method is more accurate than node method for simulating the supply temperature at substations in this case. In addition, the delay times obtained by both models are almost same, which indicates that substituting flow time of water (used in node method) for temperature wave spread time (used in function method) will not result in any significant deviation of delay time. As also seen in Figs. 6–10, the magnitude of deviation between the measured and predicted temperatures is not uniform during the simulation period. Therefore, comparisons were also carried out at different temperature change stages, namely rapid dropping stage (stage 1), rapid rising stage (stage 2) and relatively stable stage (stage 3), as presented in Table 2 for a period of 3 h, which is about the time taken for temperature to reach minimum at rapid-dropping stage, the time taken for temperature to reach maximum at rapid-rising stage or the time for relatively stable stage after rapid-rising stage. From Table 2, it can be seen that average error and standard deviation of both methods at the three stages tend to be larger with the increase in distance from HS. The same phenomenon can be seen in Figs. 6–10 and Table 1. A possible reason is that the effects of branches, turns, junctions and other fittings are not taken into account in the present models [5]. Simultaneously, the average error and standard deviation of both methods at rapid dropping stage and rapid rising stage are greater than that at relatively stable stage. This is mainly affected by the delay time: a little bit of delay time deviation has a great
Fig. 11. Simulated supply temperature of SU4 at different fluid velocities.
effect on stage of large and quick temperature changes. The accuracy at relatively stable stage is mainly affected by the relative attenuation degree and it is not as intense as the effect of delay time on stage of large temperature changes. As also seen in Table 2, the function method is more accurate than the node method at the rapid dropping stage and relatively stable stage, while it is opposite at rapid rising stage. This phenomenon is caused by several factors, such as delay time, relative attenuation degree and the trend of temperature changes. Taking Fig. 10 as an example, on the basis of overestimation of delay time and relative attenuation degree, at rapid dropping stage, average error and standard deviation increase with the increase of relative attenuation degree; on the contrary, average error and standard deviation decrease with the increase of relative attenuation degree at rapid rising stage. Other situations require a specific analysis. Additionally, the calculation time of function method for simulation was also compared with node method. The result indicates that both methods can run on a computer with normal configuration, but the function method is 37% faster than the node method on the same computer configuration. The optimization of calculation time is mainly due to reduction of the calculation steps (‘‘three steps” used in node method) by taking account of various factors in one single step (Eq. (10)).
688
J. Zheng et al. / Applied Thermal Engineering 123 (2017) 682–688
References
Fig. 12. Simulated supply temperature of SU5 at different fluid velocities.
Comparisons were also performed considering different fluid velocities at HS using numerical results. The simulated supply temperature at SU4–SU5 are shown in Figs. 11 and 12, respectively. As seen in Fig. 12, higher the fluid velocity is, faster the temperature wave propagates and less the heat loss is. This is because the increased heat loss caused by strengthened thermal convection between water and pipe is much less than the reduced heat loss due to reduction of the propagation time. It can also be seen that irrespective of higher or lower velocity, the simulation results of both methods are similar. Therefore, the analysis at fluid velocity of 1.17 m/s is also applicable for different fluid velocities and we can conclude that the proposed method can also perform well when volume flows in the network are varying by using pseudodynamic method. 4. Conclusions The main conclusions of this study were summarized as follows: 1. For function method, the average error and standard deviation over the simulation time are 0.22 °C and 0.42 °C respectively, both of which are smaller than that of node method. 2. The function method is more accurate than the node method at the rapid dropping stage and relatively stable stage, while it is opposite at rapid rising stage. 3. Compared with node method, the calculation time of function method is reduced by approximately 37%. 4. The proposed method can also perform well when volume flows in the network are varying. Therefore, the function method could be recommended to simulate the dynamic temperature of district heating network accurately and efficiently. Acknowledgements This work was supported by the National Key Technologies P&D Program of China (No. 2015BAA01B01) and the Key Projects of Natural Science Foundation of Heilongjiang Province of China (No. ZD2016010).
[1] G. Comodi, M. Lorenzetti, D. Salvi, A. Arteconi, Criticalities of district heating in Southern Europe: lesson learned from a CHP-DH in Central Italy, Appl. Therm. Eng. 112 (2016) 649–659. [2] Z. Pan, Q. Guo, H. Sun, Feasible region method based integrated heat and electricity dispatch considering building thermal inertia, Appl. Energy (2016). [3] I.B. Hassine, U. Eicker, Impact of load structure variation and solar thermal energy integration on an existing district heating network, Appl. Therm. Eng. 50 (2013) 1437–1446. [4] S.A. Kyriakis, P.L. Younger, Towards the increased utilization of geothermal energy in a district heating network through the use of a heat storage, Appl. Therm. Eng. 94 (2015) 99–110. [5] I. Gabrielaitiene, B. Bøhm, B. Sunden, Modelling temperature dynamics of a district heating system in Naestved, Denmark—a case study, Energy Convers. Manage. 48 (2007) 78–86. [6] I. Gabrielaitiene, B. Bøhm, B. Sunden, Dynamic temperature simulation in district heating systems in Denmark regarding pronounced transient behaviour, J. Civil Eng. Manage. 17 (2011) 79–87. [7] M. Kuosa, K. Kontu, T. Mäkilä, M. Lampinen, R. Lahdelma, Static study of traditional and ring networks and the use of mass flow control in district heating applications, Appl. Therm. Eng. 54 (2013) 450–459. [8] X. Sheng, D. Lin, Energy saving analyses on the reconstruction project in district heating system with distributed variable speed pumps, Appl. Therm. Eng. 101 (2016) 432–445. [9] P. Jie, N. Zhu, D. Li, Operation optimization of existing district heating systems, Appl. Therm. Eng. 78 (2015) 278–288. [10] A. Colmenar-Santos, E. Rosales-Asensio, D. Borge-Diez, E. Collado-Fernández, Evaluation of the cost of using power plant reject heat in low-temperature district heating and cooling networks, Appl. Energy 162 (2016) 892–907. [11] S. Lim, S. Park, H. Chung, M. Kim, Y.J. Baik, S. Shin, Dynamic modeling of building heat network system using Simulink, Appl. Therm. Eng. 84 (2015) 375–389. [12] O.P. Palsson, Stochastic modeling, control and optimization of district heating systems, 1993. [13] H. Zhao, Analysis, Modelling and Operational Optimization of District Heating Systems, Technical University of Denmark, 1995. [14] H. Palsson, Analysis of numerical methods for simulating temperature dynamics in district heating pipes, in: International Symposium on District Heating and Cooling Simulation, 2010, pp. 219–224. [15] A. Benonysson, Dynamic modelling and operational optimization of district heating systems, Laboratory of Heating and Air Conditioning, Technical University of Denmark, 1991. [16] V.D. Stevanovic, B. Zivkovic, S. Prica, B. Maslovaric, V. Karamarkovic, V. Trkulja, Prediction of thermal transients in district heating systems, Energy Convers. Manage.ent 50 (2009) 2167–2173. [17] I. Gabrielaitiene, B. Bøhm, B. Sunden, Evaluation of approaches for modeling temperature wave propagation in district heating pipelines, Heat Transfer Eng. 29 (2008) 45–56. [18] H. Zhao, Optimum operation of a CHP-based district heating system by mathematical modelling, Fuel Energy Abstr. (1996) 296. [19] L. Dobos, J. Jäschke, J. Abonyi, S. Skogestad, Dynamic model and control of heat exchanger networks for district heating, Hung. J. Ind. Chem. 37 (2009) 37–49. [20] H.V. Larsen, M. Lucht, M. Wigbels, Simple models for operational optimization, in: IEA District heating and cooling, Technical University of Denmark, 2002. [21] MOHURD, Design code for city heating network, in: Vol. CJJ34 - 2010, China Architecture & Building Press, Beijing, 2010. [22] Z. Zhang, Engineering Mathematics Analysis, China Higher Education Press, 2008. [23] Q. Bing, Studies on the dynamic thermal characteristics of district heating systems, Tsinghua University, 2004, pp. 162. [24] P. Jie, Z. Tian, S. Yuan, N. Zhu, Modeling the dynamic characteristics of a district heating network, Energy 39 (2012) 126–134. [25] J. Duquette, A. Rowe, P. Wild, Thermal performance of a steady state physical pipe model for simulating district heating grids with variable flow, Appl. Energy 178 (2016) 383–393. [26] K. Çomaklı, B. Yüksel, Ö. Çomaklı, Evaluation of energy and exergy losses in district heating network, Appl. Therm. Eng. 24 (2004) 1009–1017. [27] L. Hua, Y. Wang, Numerical Integration and Application, Science Press, Beijing, 1963. [28] J.A. Simmons, Dynamical Prediction: Some Results from Operational Forecasting and Research Experiments, World Meteorological Organization, 1983. [29] A. Xiu, D. Liao, Study on the determination of optimal Fourier series, Acta Meteorol. Sin. 49 (1991) 308–320. [30] P. He, G. Sun, F. Wang, Heating Engineering, fourth ed., China Architecture & Building Press, Beijing, 2009.