An entransy dissipation-based method for global optimization of district heating networks

Energy and Buildings 48 (2012) 50–60

Contents lists available at SciVerse ScienceDirect

Energy and Buildings journal homepage: www.elsevier.com/locate/enbuild

An entransy dissipation-based method for global optimization of district heating networks Yun-Chao Xu, Qun Chen ∗ Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China

a r t i c l e

i n f o

Article history: Received 30 September 2011 Received in revised form 22 December 2011 Accepted 8 January 2012 Keywords: District heating network Optimization Entransy dissipation Energy efficiency

a b s t r a c t Improving the heating performance of district heating networks (DHNs), one of the most important components in district heating systems, has significant influence on the promotion of energy efficiency. Based on the entransy theory, this paper develops an entransy dissipation-based method for global optimization of DHNs. We analyze the irreversibility of three types of heat transfer processes in DHNs: (1) the heat transfer process between the waters in primary heating networks and DHNs in substation heat exchangers; (2) the heat transfer process between supply water and indoor airs in radiators; and (3) the return water mixing process, then derive the expressions of entransy dissipation for each process, and finally deduce two different formulas of the total entransy dissipation in DHNs by the users’ demands and the structural and operating parameters of DHNs, respectively, from which we establish the mathematical relation between the users’ demands and the structural and operating parameters of DHNs. Based on this relation, two typical optimization problems of DHNs are converted into conditional extremum problems to deduce two optimization equation groups. Solving the equation groups gives the optimal structural and operating parameters for DHNs with highest energy efficiency. Finally, the method is proved valid and advantageous through the optimization for both 1-column-2-user and 2-column-4-user DHNs. © 2012 Elsevier B.V. All rights reserved.

1. Introduction District heating systems (DHSs) are widely applied in a large number of buildings to supply heat in winter. Improving the heating performance of DHSs will promote the energy efficiency and thereafter reduce the building energy consumption. Therefore, during the last several decades, several scholars developed a lot of linear program optimization methods for DHSs including the market allocation model (MARKAL) [1,2], the Brookhaven energy system optimization model (BESOM), the time-stepped energy system optimization model (TESOM) [3], and the model for optimization of dynamic energy system with time-dependent components and boundary conditions (MODEST) [4]. These models are employed in the optimization for different types of DHSs and proved valid on not only reducing system cost but also promoting energy efficiency. However, most models have a common inherent shortage that they perform well on the brief optimization of the whole systems but are not very good at the specific optimization for the component of DHSs like district heating network (DHN), one of the most important components in DHSs.

∗ Corresponding author. Tel.: +86 10 62781610; fax: +86 10 62781610. E-mail address: [email protected] (Q. Chen). 0378-7788/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.enbuild.2012.01.008

In order to improve the performance of DHNs, Bojic et al. [5,6] optimized the structural and operating parameters in a DHN based on the energy conservation equations with the goal of the least temperature loss. And Soderman and Pettersson [7] optimized a DHN with the objective of the minimum the total cost. Besides, the search method has also attracted more and more attentions. These methods aim at finding the optimal result through the comparison among all the possible combinations of unknown parameters. For instance, Baker and Sherif [8] optimized the heat transfer performance of a DHN with a conventional but low-efficiency search method. In order to promote the efficiency of optimization process, the genetic algorithm (GA) is applied in the search method for DHN optimization. Cammarata et al. [9] used the GA to design heat exchangers, water pumps and air compressors in a building to improve their efficiency. Wright et al. [10] optimized an existing DHN to minimize the running cost, and then compared the new configuration of all devices with the original one before optimization. However, they cannot globally optimize the performance of DHNs either. It is because that the global optimization has to establish the theoretical and mathematical relation of the demands of users and all the structural and operating parameters of DHN, i.e. the power consumption of every pump and the heat transfer area of each heat exchanger to minimize both the running and investment cost. Due to the lack of this relation, the search, actually trial-anderror, method has to be employed, which means to list all possible

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

Nomenclature a A cp F m g G k (kA)0

heat distribution ratio area, m2 constant pressure specific heat, J/(kg K) Lagrange function mass flux, kg/s constant of heat capacity rate, W/K heat capacity rate, W/K heat transfer coefficient, W/(m2 K) thermal conductance of the substation heat exchanger, W/K (kA)ij thermal conductance of the radiator ij, W/K q heat flux, W/m2 Q heat transfer rate, W Rh entransy dissipation-based thermal resistance, K/W T temperature, K temperature of supply water in primary heating netT works, K temperature of return water in primary heating netT works, K T temperature difference, K u constant of thermal conductance, W/K velocity vector, m/s v V volume, m3 ˚ entransy dissipation, W K  thermal conductivity, W/(m K) density, kg/m3  ˛, i , ij ,  Lagrange multipliers Subscripts b return water cold c h hot, entransy dissipation-based i inlet o outlet m after mixing maximum max min minimum s supply water water w * primary heating network’s

answers of the problem and suffer a huge amount of calculations to achieve a better, not the best solution. Actually, DHNs consist of many irreversible heat transfer processes. For heat transfer, Guo et al. [11,12] introduced a physical quantity, termed entransy, to represent the heat transfer ability of an object or a system, deduced the expression of entransy dissipation to measure the irreversibility of heat transfer, and furthermore developed the minimum entransy dissipation-based thermal resistance principle [12] to optimize all the three basic ways of heat transfer, i.e. heat conduction [12–14], heat convection [15–18], and thermal radiation [19,20]. Besides, in practical engineering applications, this principle has been also used in the optimization of heat transfer facilities [21,22], heat exchangers [23,24] and evaporative cooling systems [25,26]. The purpose of this contribution is to propose an entransy dissipation-based optimization method for global optimization of DHNs to effectively improve their performance. We divide the irreversible heat transfer processes in DHNs into three categories, derive the expression of entransy dissipation for each category, obtain the total entransy dissipation in DHNs to establish the mathematical and theoretical relation between the demands of users and

51

the structural and operating parameters of the devices in DHNs, and finally deduce the optimization equation groups for DHNs through conditional extreme method. Moreover, the applications of this new method are demonstrated in the optimization of both a 1column-2-user and 2-column-4-user DHNs. 2. Entransy dissipation in DHNs Fig. 1 shows the sketch of a typical DHS, which consists of a primary heating network (PHN) and a district heating network (DHN), where the PHN and the DHN are connected by the substation heat exchanger. The DHN consists of a substation heat exchanger, a water pump, pipes and users’ radiators in n buildings. The building i (i = 1, 2,. . ., n) is named Bi . In each building Bi , mi stands for the amount of users in this building and all the users are numbered from 1 to mi . The user ij mentioned in the following represent the jth user in Bi . T , T , Ts and Tb stand for the temperatures of the supply and return waters in the PHN and the DHN, respectively. G* and G0 are the heat capacity rates of the water in the PHN and the DHN, respectively, where the heat capacity rate of the water in Bi is Gi  and Gi = G0 . Heat is transferred from the heat source to the building radiators through the PHN and the DHN. In the PHN, hot water flows into the substation heat exchanger to heat the water in the DHN and its temperature is decreased from T to T . Meanwhile, the water in the DHN is heated from Tb to Ts . With the driving of the pump, the heated water in the DHN flows into each building. In each building, the water flows into the radiators of each user in order and heats the indoor air of these users. After heating all the indoor airs in different buildings, the cold water mixes and returns to the substation heat exchanger. During the whole heating process, there are such three irreversible heat transfer processes in the system as (1) the heat transfer process between the waters in the substation, (2) the heat transfer process between the water and the indoor air in each radiator, and (3) the return water mixing process. Taking a 2-column-4-user DHN which contains two buildings, B1 and B2 , with two users in each for example, Fig. 2 shows the temperature variations of each working fluid versus heat flow rate q during all the heat transfer processes. Qij stands for the heat demand of user ij. The shadow trapezium ABCD represents the heat transfer between waters in the substation heat exchanger, the trapezium DILK and EFGH separately represent the heat transfer processes between the water and the indoor air in the radiators of B1 and B2 , and the triangle HIJ represents the return water mixing process between return waters from two buildings. 2.1. Entransy dissipation in the substation For a steady-state convective heat transfer process without internal heat source, the thermal energy conservation equation is: cp v · ∇ T = ∇ · (T ),

(1)

where, , cp , v,  are the density, constant pressure specific heat, fluid velocity vector and thermal conductivity, respectively. Multiplying both sides of Eq. (1) by the temperature T yields the entransy balance equation [18]:



v·∇

cp T 2 2



= ∇ · (T ∇ T ) − |∇ T |2 .

(2)

Comparing Eq. (1) with Eq. (2), there is an additional term, | ∇ T|2 , on the right hand of Eq. (2), which represents the local entransy dissipation during heat transfer. That is, the energy is conserved, whereas the entransy is not always conserved but dissipated during convective heat transfer. Integrating the local

52

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

Fig. 1. The sketch of a DHS.

entransy dissipation over the entire heat transfer domain yields the total entransy dissipation, i.e. the irreversibility measurement of heat transfer, which is expressed as:







|∇ T |2 dV =

˚=

∇ · (T ∇ T )dV −

˝

˝

v·∇



cp T 2 2

˝

 dV.

(3)

By applying Gauss’s law, Eq. (3) is simplified as [24,25]: ˚=

2 in   mcp T 2  − 2

1

out

 ·q  )TdA, (n

(4)

 is heat flux through the where m is the mass flow rate of fluid and q boundary. In Eq. (4), the first term on the right means the net entransy flowing into the domain accompanying the fluid flow, while the second term represents the entransy transferred through the boundary of the domain due to heat conduction. Therefore, for the convective heat transfer process in the substation heat exchanger in Fig. 1, the expressions of entransy dissipation for both the hot and cold waters are respectively written as: 1 2 2 ˚h = mh cp,h (T  − T  ) − 2 ˚c =

1 mc cp,c (Tb2 − Ts2 ) + 2



→ → (− n ·− q )TdA,

(5a)

∂˝

→ → (− n ·− q )TdA,

˚a =

1 1 2 2 m c (T  − T  ) + mc cp,c (Tb2 − Ts2 ), 2 h p,h 2

(6)

where ˚a is the entransy dissipation in the substation heat exchanger. According to the energy conservation equation, Q = mh cp,h (T  − T  ) = mc cp,c (Ts − Tb ).

∂˝



Because the entransy flowing out of the hot fluid equals to that into the cold fluid at the interface between the hot and cold fluids, combining Eqs. (5a) and (5b) yields:

(5b)

∂˝

where mh and mc are the mass flow rates of the hot and cold fluids, and cp,h and cp,c are the constant pressure specific heats of hot and cold fluids, respectively.

(7)

Eq. (6) is simplified as: ˚a =

1 Q (T  − Ts + T  − Tb ), 2

(8)

where Q stands for the total heat demand of the users in the DHN. In Fig. 2, the shadow trapezium ABCD can be calculated by Eq. (8). The shadow area represents the entransy dissipation during the heat transfer in the substation heat exchanger. On the other hand, the shadow area is also calculated by the following integral equation:



˚=

Q

(Th − Tc )dq,

(9)

0

where Th and Tc stand for the temperatures of the hot and cold fluids in the heat exchanger and Q for the total heat transferred between two fluids, respectively. By analogy with electric conduction [12], the entransy dissipation-based thermal resistance (EDTR) of heat exchanger is then defined as: Rh =

˚ . Q2

(10)

Substituting the heat transfer equations for counter flow heat exchangers, Q = (kA)0

Tmax − Tmin , ln(Tmax /Tmin )

(11)

Tmax ∗ = e−(kA)0 ((1/G )−(1/G0 )) , Tmin

(12)

into Eq. (9), the entransy dissipation in the substation heat exchanger is deduced as:

Fig. 2. The T–q graph for all the heat transfer processes in the DHN.

˚a =

1 2 Q 2

1 G∗

−

1 G0

 e(kA)0 ((1/G∗ )−(1/G0 )) + 1 e(kA)0 ((1/G

∗ )−(1/G )) 0

−1

.

(13)

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

53

process. The heat transfer rate between the two return waters Qm is: Qm = G1 (Tb1 − Tb ).

(18)

Based on the energy conservation equation, Tb1 and Tb2 are calculated by: Tb1 = Ts −

Q11 + Q12 , G1

(19a)

Q21 + Q22 , G2

(19b)

and Fig. 3. Entransy dissipation of heat transfer in radiator ij.

Tb2 = Ts −

And the expression of EDTR of the substation heat exchanger is: Rh =

 1 1 G∗

2

 e(kA)0 ((1/G∗ )−(1/G0 )) + 1 1

−

G0

e(kA)0 ((1/G

∗ )−(1/G )) 0

−1

,

(14)

where (kA)0 stands for the thermal conductance of the substation heat exchanger. Tmax and Tmin stand for the maximum and minimum temperature differences between the hot and cold waters on the same side of the substation heat exchanger, respectively. 2.2. Entransy dissipation in the radiators

1 Qm (Tb1 − Tb2 ). 2

(20)

Substituting Eqs. (18), (19a), (19b) into Eq. (20) gives the entransy dissipation of the return water mixing process:

⎛

1 2

2

2

Q j=1 1j

⎜ ⎝



2 Q j=1 2j

+

G1

G2

⎞

2

−

Q2 G0

⎟ ⎠.

(21)

Similarly, the entransy dissipation during the mixing process in the DHN, ˚c , is expressed as: 1 ˚c = 2

 n Q2

Q2 − G0 Gi i

i=1

where Qi =

mi j=1



,

(22)

Qij .

2.4. The total entransy dissipation in the DHN

1 2 1 e(kA)ij /Gi + 1 Q , 2 ij Gi e(kA)ij /Gi − 1

(15)

where ˚ij represents the entransy dissipation in the radiator ij, Rij is the EDTR of the radiator, Qij stands for the heat demand of the user ij and (kA)ij stands for the thermal conductance of the radiator ij. The total entransy dissipation of Bi is the sum of the entransy dissipation in all radiators of the building. By combining Eq. (15) of all the users in Bi , the total entransy dissipation in Bi , ˚i , is calculated by: 1 1 2 e(kA)ij /Gi + 1 . Q Gi 2 ij e(kA)ij /Gi − 1 mi

˚i =

˚c =

˚c =

Fig. 3 shows the T–q graph for the heat transfer process in each radiator, where Ti stands for the indoor air temperature of users in Bi , Twi and Two stand for the temperatures of water flowing into and out of the radiators, respectively. Because the temperature of indoor air always keeps constant, the heat capacity rate of indoor air can be regarded as infinity. According to a similar analysis in Section 2.1 and based on Eq. (13), the expression of entransy dissipation in the radiator ij is derived as: ˚ij = Qij2 Rij =

respectively. Calculating the shadow triangle yields the entransy dissipation of return water mixing process:

(16)

j=1

The total entransy dissipation in the DHN is the sum of all the three parts of entransy dissipation calculated in Sections 2.1–2.3, which is expressed as: ˚ = ˚a + ˚b + ˚c =

1

G∗

1 n

+

1 2 Q 2

mi

i=1 j=1

−

1 G0

 e(kA)0 ((1/G∗ )−(1/G0 )) + 1 e(kA)0 ((1/G∗ )−(1/G0 )) − 1

 n 2 1 Q

1 e(kA)ij /Gi + 1 + Qij2 2 Gi e(kA)ij /Gi − 1 2

Q2 − Gi G0



i

i=1

= ˚(Gi , (kA)0 , (kA)ij ).

(23)

And in Fig. 2, the shadow trapezium DILK and EFGH separately represent the entransy dissipation of B1 and B2 . The total area of both trapeziums equals to the total entransy dissipation of the radiators in the DHN.

Eq. (23) indicates that the total entransy dissipation of DHN is a function of the structural and operating parameters. On the other hand, the entransy dissipation of each heat transfer process is represented by the related shadow area in Fig. 2. So the total entransy dissipation in the DHN corresponds to the total area of all the shadow parts. In Fig. 2, it is obvious that the area of the triangles DMJ and EMC is equivalent. By filling the triangle DMJ with the triangle EMC, all the shadow parts in Fig. 2 is equally converted into the shadow area in Fig. 4. Through calculating the total area of the shadow parts in Fig. 4, the total entransy dissipation in the DHN is expressed by means of the users’ demands as:

2.3. Entransy dissipation during the mixing process of return water

1 ⎝ Qij ⎠ Ti = ˚(T  , T  , Ti , Qij ). ˚ = (T  + T  )Q − 2

Similar with the analysis above, the area of shadow triangle HIJ in Fig. 2 represents the entransy dissipation of return water mixing

Eq. (24) indicates the relation between the users’ demands and the entransy dissipation in the DHN. When the users’ demands and

By adding up the entransy dissipation of all the buildings, the total entransy dissipation of the radiators in the DHN, ˚b , can be calculated as: ˚b =

n i=1

˚i =

mi n 1 1 i=1

Gi

2 j=1

Qij2

e(kA)ij /Gi + 1 e(kA)ij /Gi − 1

.

(17)

n

i=1

⎛

mi

⎞

(24)

j=1

54

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

One constraint is that the users in the same building are seriesconnected. So the relation between the temperature for the water flowing into the radiator ij and the temperature of the supply water in DHN can be expressed as: (ij)

Tw = Ts −

j−1 Qij

Gi

j=1

,

(28)

(ij)

where Tw stands for the temperature for the water flowing into the radiator ij. Similar with Eq. (8), the entransy dissipation for the radiator ij can be calculated as:

Fig. 4. Total entransy dissipation in the DHN.

˚ij = the temperatures of the supply and return waters in the primary heating network are given, the total entransy dissipation in the DHN is proved to be constant. And then combining Eqs. (23) and (24) yields the mathematical and theoretical relations between the structural and operating parameters and the users’ demands; ˚(Gi , (kA)0 , (kA)ij ) = ˚0 (T  , T  , Ti , Qij ).



(ij)

2Tw −

Qij Gi



− 2Ti

.

(29)

Combining Eqs. (15) and (31) yields: 1 1 2 1 e(kA)ij /Gi + 1 = Qij Q 2 2 ij Gi e(kA)ij /Gi − 1



(ij) 2Tw

Qij

−

Gi

 − 2Ti

.

(30)

Substituting Eq. (30) into (32) yields:

⎛

(25)

This relation mathematically connects all the unknown structural and operating parameters, i.e. the heat capacity rates of waters and the thermal conductance of heat exchangers, with the known condition, i.e. the users’ demands for all kinds of DHNs with different structures and users’ demands. Therefore, this relation makes it possible to globally optimize all the structural and operating parameters at the same time, which is obviously superior to the conventional trial-and-error method.

1 Q 2 ij

j−1 Qij

(kA)ij /Gi

1 21 e +1 1 Q = Qij ⎝2Ts − 2 2 ij Gi e(kA)ij /Gi − 1 2

j=1

Gi

⎞ −

Qij Gi

− 2Ti ⎠ .

(31)

By solving Eq. (33), new expressions of the temperature for the (ij) supply water, Ts , are obtained as: Qij 1 1 e(kA)ij /Gi + 1 Qik + + + Ti . Q 2 ij Gi e(kA)ij /Gi − 1 Gi 2Gi j−1

(ij)

Ts

=

(32)

k=1

3. Optimization models for DHNs

Then new constraints can be added into the model: (ij)

Ts = Ts , The cost of any DHN consists of two parts: one is the structural cost mainly related to the area of heat exchangers, and the other is the operating cost related to the energy consumption of pumps. Based on these two parts of DHN’s cost, the optimization objective of DHN for minimizing the DHN’s cost is mainly applied in two aspects. On one hand, as the low-energy building concept becomes more and more popular [27], the energy consumption of pumps is usually given, the minimization of the total cost of DHN can be achieved by minimizing the heat exchangers’ cost, which is corresponding to minimizing the total area of all heat exchangers. On the other hand, if the total area of all heat exchangers in DHN is given after primary design, the minimization of the total cost of DHN can be achieved by minimizing the total energy consumption of pumps. In order to solve these two problems, it is necessary to build a mathematical model to not only reflect the physical characteristic of the process, but also satisfy the constraints of all real situations. For a practical DHN, one constraint is that the supply water flowing into different buildings has the same temperature Ts which can be calculated from the expression of the entransy dissipation for the substation heat exchanger: 1 2 Q 2

1

1 − G∗ G0

 e(kA)0 ((1/G∗ )−(1/G0 )) + 1

1 = Q [(T  − Ts )+(T  − Tb )]. 2 e(kA)0 ((1/G∗ )−(1/G0 )) − 1 (26)

By solving this equation, Ts is expressed as: Ts =

1 1

1  Q − Q (T + T  ) + 2 2G0 2

G∗

−

 e(kA)0 ((1/G∗ )−(1/G0 )) + 1 1 G0

. e(kA)0 ((1/G∗ )−(1/G0 )) − 1 (27)

j= / 1.

(33)

On the other hand, the temperatures of the water flowing into different buildings are the same, which is expressed as: (11)

Ts

(k1)

= Ts

,

k = 2, 3, . . . , n.

(34)

With all the aforementioned constraints, two typical optimization problems can be converted into the conditional extremum problems as following: (1) When the total energy consumption of pumps is given, n

Gi = g,

(35)

i=1

where g is constant. The optimization objective of DHN is to minimize the total area of all heat exchangers. Because the heat transfer coefficient can be assumed constant, the optimization object is changed to minimize the total thermal conductance, which is stated as below: n i=1

Gi = g,

˚ = ˚0 ,

⎛

find min ⎝

n mi

(11)

Ts

(k1)

= Ts

,

(ij)

Ts = Ts ,

⎞ (kA)ij + (kA)0 ⎠ .

i=1 j=1

The constraints can be removed by using Lagrange multipliers method to construct a function:

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

⎛ F =⎝

mi n

⎞ (kA)ij + (kA)0 ⎠ + ˛(˚(Gi , (kA)0 , (kA)ij ) − ˚0 ) −

 n 1 i

i=1 j=1 mi n

−

 1  Q 1 − Q (T + T  ) + 2 2G0 2

ij

i=1 j=2

1

1 − G∗ G0

Q11

2G1

i=2

55

e(kA)11 /G1 + 1 e(kA)11 /G1

Q e(kA)i1 /Gi + 1 1 Q11 − i1 − Ti + T1 − Q + 2G1 2Gi i1 e(kA)i1 /Gi − 1 2Gi −1

 e(kA)0 ((1/G∗ )−(1/G0 )) + 1

1 1 e(kA)ij /Gi + 1 Qik Qij − Qij − − Ti − (kA) ((1/G )−(1/G )) (kA)ij /Gi ∗ 2 G Gi 2Gi 0 0 e −1 i e −1 j−1

 +

k=1

 n



 Gi − g

,

i=1

(36) where ˛, i , ij and  are Lagrange multipliers. The differential of F with respect to the thermal conductance of substation heat exchanger, (kA)0 , offers



⎛

n i=1

∂F = 1 − ⎝˛ + ∂(kA)0

mi

 j=2 ij

Q

⎞

 2 ⎠ Q2 1 − 1 G∗

G0

e(kA)0 ((1/G [e

∗ )−(1/G )) 0

(kA)0 ((1/G∗ )−(1/G0 ))

− 1]

= 0.

2

(37)

The differential of F with respect to the thermal conductance of each radiator, (kA)ij , yields Q2 e(kA)11 /G1 ∂F = 1 − ˛ 11 + 2 ∂(kA)11 G1 [e(kA)11 /G1 − 1]2

 n  Q11 i

e(kA)11 /G1

G12

i=2

[e(kA)11 /G1 − 1]

Q2 ∂F e(kA)11 /G1 Q e(kA)i1 /Gi = 1 − ˛ 11 − i i1 = 0, 2 2 2 ∂(kA)i1 G1 [e(kA)11 /G1 − 1] Gi [e(kA)i1 /Gi − 1]2



∂F =1− ∂(kA)ij

˛+

ij



Qij

Qij2

e(kA)ij /Gi

Gi2

[e(kA)ij /Gi − 1]

2

= 0,

2

= 0,

(38a)

i = 2, 3, ..., n,

i = 1, 2, . . . , n,

(38b)

j = 2, . . . , mi .

(38c)

The differential of F with respect to the heat capacity rate in each branch, Gi , provides



∂F =+ ∂G1

m1

˛+

⎛ +˛⎝ mi

1j



−˛⎝

⎛ ⎝

∗ )−(1/G )) 0

− 1 − 2(kA)0 ((1/G∗ ) − (1/G0 ))e(kA)0 ((1/G (e(kA)0 ((1/G

(e(kA)1j /G1 − 1)

G12

2(e(kA)1j /Gi − 1)

 j=2 ij



Q 2 e2(kA)0 ((1/G

2

∗ )−(1/G ) 0

2G02

Q

2Gi2

⎞ ij ⎠

j=2

2

2G1

+

(e(kA)ij /Gi − 1)

2

2

G1

i=2



2G0 2

2Gi2

+

Q1j

= 0,

2G12

∗ )−(1/G ) 0



⎠ − i

+

− 1)

(39a)

ik

Gi2

n

∗ )−(1/G ) 0

2

+

Q2



2G02

Qi1 e2(kA)i1 /Gi − 1 − 2((kA)i1 /Gi )e(kA)i1 /Gi 2Gi2

j Q k=1

2G1

2(e(kA)11 /Gi − 1)



1k G12 k=1

⎞ +

+

Q2

⎞    n Q11 Q11 e2(kA)11 /G1 − 1 − 2((kA)11 /G1 )e(kA)11 /G1 ⎠+ i + 2 2 2 2

j−1 Q

Qi2

Qij e2(kA)ij /Gi − 1 − 2((kA)ij /Gi )e(kA)ij /Gi 2Gi2

− 1)

∗ )−(1/G )) 0

Q1 2

(e(kA)0 ((1/G

(e(kA)ij /Gi − 1)

∗ )−(1/G )) 0

− 1 − 2(kA)0 ((1/G∗ ) − (1/G0 )e(kA)0 ((1/G

mi Qij2 e2(kA)ij /Gi − 1 − 2((kA)ij /Gi )e(kA)ij /Gi

mi

+

2

Q1j e2(kA)1j /G1 − 1 − 2((kA)1j /G1 )e(kA)1j /G1

mi

j=1

−

2G12

˛+

⎛

Q 2 e2(kA)0 ((1/G 2G02

Q



j=2

∂F =+ ∂Gi



mi 2 Q1j e2(kA)1j /G1 − 1 − 2((kA)1j /G1 )e(kA)1j /G1 j=1

−

 j=2 1j

(e(kA)i1 /Gi − 1)

2

+

Qi1



2Gi2

 +

Qij 2Gi2

= 0,

i= / 1,

(39b)

n

Eqs. (25), (33)–(35) and (37)–(39) have 2 i=1 mi + n + 2 unknown parameters as well as 2 i=1 mi + n + 2 equations. Therefore, solving the equations will obtain the values of all the unknown parameters especially the optimal arrangements of the thermal conductance and heat capacity rates to minimize the total thermal conductance of all the heat exchangers. (2) When the total thermal conductance of all heat exchangers is given, mi n

(kA)ij + (kA)0 = u

i=1 j=1

where u is constant.

(40)

56

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

The optimization objective of DHN is to minimize the heat capacity rates of waters which can be stated below: mi n

(kA)ij + (kA)0 = u,

Ts11

˚ = ˚0 ,

j1 Ts ,

=

Ts =

(ij) Ts ,

find min

 n  Gi

i=1 j=1

.

i=1

The constraints can also be removed by using Lagrange multipliers method to construct a function: F =

 n  Gi

+ ˛(˚(Gi , (kA)0 , (kA)ij ) − ˚0 ) −

i=1

−

mi n

Q 1 1  − Q (T + T  ) + 2 2G0 2

ij

⎛

+ ⎝

n mi

i

1

1 − G∗ G0

e(kA)11 /G1 + 1

Q11

2G1

i=2



i=1 j=2

 n 1

e(kA)11 /G1 − 1

+

Q 1 e(kA)i1 /Gi + 1 Q11 − i1 − Ti +T − Q 2G1 2Gi i1 e(kA)i1 /Gi − 1 2Gi

 e(kA)0 ((1/G∗ )−(1/G0 )) + 1

Qij 1 1 e(kA)ij /Gi + 1 Qik − Qij − − − Ti 2 Gi e(kA)ij /Gi − 1 Gi 2Gi e(kA)0 ((1/G∗ )−(1/G0 )) − 1 j−1

 

k=1

⎞ (kA)ij + (kA)0 − u⎠

(41)

i=1 j=1

where ˛, i , ij and  are Lagrange multipliers. The differential of F with respect to the thermal conductance of substation heat exchanger, (kA)0 , offers



⎛

n i=1

∂F =  − ⎝˛ + ∂(kA)0

mi

 j=2 ij

Q

⎞

 2 ⎠ Q2 1 − 1 G∗

G0

e(kA)0 ((1/G [e

∗ )−(1/G )) 0

(kA)0 ((1/G∗ )−(1/G0 ))

− 1]

=0

2

(42)

The differential of F with respect to the thermal conductance of each radiator, (kA)ij , yields Q2 e(kA)11 /G1 ∂F =  − ˛ 11 + 2 ∂(kA)11 G1 [e(kA)11 /G1 − 1]2

 n  Q11 i

e(kA)11 /G1

G12

i=2

[e(kA)11 /G1 − 1]

Q2 ∂F e(kA)11 /G1 Q e(kA)i1 /Gi =  − ˛ 11 − i i1 = 0, 2 2 2 ∂(kA)i1 G1 [e(kA)11 /G1 − 1] Gi [e(kA)i1 /Gi − 1]2



∂F =− ∂(kA)ij

˛+

ij



Qij2

e(kA)ij /Gi

Gi2 [e(kA)ij /Gi − 1]2

Qij

= 0,

2

=0

(43a)

i = 2, 3, ..., n

i = 1, 2, . . . , n,

(43b)

j = 2, . . . , mi .

(43c)

The differential of F with respect to the heat capacity rate in each branch, Gi , provides



∂F = 1+ ∂G1

m1

˛+

⎛ +˛⎝ mi

1j



−˛⎝

−

⎝

∗ )−(1/G )) 0

− 1 − 22(kA)0 ((1/G∗ ) − (1/G0 ))e(kA)0 ((1/G (e(kA)0 ((1/G

(e(kA)1j /G1 − 1)

mi j=2

+

2

Q1j e2(kA)1j /G1 − 1 − 2((kA)1j /G1 )e(kA)1j /G1 G12

2(e(kA)1j /G1 − 1)

 j=2 ij



Q 2 e2(kA)0 ((1/G

2

∗ )−(1/G )) 0

2G02

Q

2Gi2

⎞ ij ⎠

2

(e(kA)ij /Gi − 1)

2

2

2G1

i=2

j−1 Q

Q1j

1k G12 k=1

Q2



2G02

G1

Qi2 2Gi2

⎠ − i

+

j Q

+



ik Gi2 k=1

n

2G1

2(e(kA)11 /G1 − 1)

 = 0,

2G12

∗ )−(1/G )) 0

⎞ +

+

⎞    n Q11 e2(kA)11 /G1 − 1 − 2((kA)11 /G1 )e(kA)11 /G1 Q11 ⎠+  + i 2 2 2 2

+

Qij e2(kA)ij /Gi − 1 − 2((kA)ij /Gi )e(kA)ij /Gi 2Gi2

− 1)

∗ )−(1/G )) 0

Q12

(e(kA)0 ((1/G

(e(kA)ij /Gi − 1)

∗ )−(1/G )) 0

(44a)

− 1 − 2(kA)0 ((1/G∗ ) − (1/G0 ))e(kA)0 ((1/G

mi Qij2 e2(kA)ij /Gi − 1 − 2((kA)ij /Gi )e(kA)ij /Gi j=1

⎛

2G12

mi

˛+

⎛

Q 2 e2(kA)0 ((1/G 2G02

Q



j=2

∂F = 1+ ∂Gi



mi 2 Q1j e2(kA)1j /G1 − 1 − 2((kA)1j /G1 )e(kA)1j /G1 j=1

−

 j=2 1j

− 1)

∗ )−(1/G )) 0

2

+

Q2



2G02

Qi1 e2(kA)i1 /Gi − 1 − 2((kA)i1 /Gi )e(kA)i1 /Gi 2Gi2

(e(kA)i1 /Gi − 1)

2

+

Qi1



2Gi2

 +

Qij 2Gi2

= 0,

i= / 1

(44b)

n

Similarly, Eqs. (25), (33), (34), (40) and (42)–(44) also have 2 i=1 mi + n + 2 unknown parameters as well as 2 i=1 mi + n + 2 equations. Therefore, through the solution of optimization equations, we can get the optimal configuration of all the unknown parameters especially the optimal arrangements of thermal conductance and heat capacity rates to minimize the total heat capacity rate of waters.

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

57

Table 1 Optimized results with entransy dissipation optimization method. Parameters (W/K) Optimization results

(kA)0 619.1

(kA)1 315.5

(kA)2 408.0

In summary, the above optimization equations indicate it is convenient and advantageous to optimize the heating performance of DHNs theoretically with the entransy dissipation-based method. Furthermore, different from conventional methods, this method does not need to give all the possible combinations of unknown parameters which lead to huge calculation amount but low accuracy. 4. Applications of the entransy dissipation-based method 4.1. Optimization of a 1-building-2-user DHN For the 1-column-2-user DHN which contain one building with two users in Fig. 1, the heat capacity rate of the supply water in the DHN, G0 , keeps 1400 W/K. The temperatures of the supply and return waters, T and T , in the primary heating network are assumed 393 K and 333 K [28,29], respectively, and the temperatures of the indoor airs for both users are 292 K. Q1 and Q2 represent the heat demands of users 1 and 2, respectively. Q1 is supposed to be a constant, 10 kW. When user 2 has the same heat demand Q2 as user 1, Q1 , Table 1 lists the optimized results, where (kA)1 and (kA)2 represent the thermal conductance of radiators for users 1 and 2, respectively. As listed in Table 1, (kA)2 is larger than (kA)1 , which means that the downstream user should have a larger radiator than that of the upstream user when the heat demands are the same. When the heat demand of user 2, Q2 , varies from 5 kW to 20 kW, the optimized results are showed in Fig. 5. Fig. 5(a) shows the optimized total thermal conductance versus the variation of Q2 . (kA)ED represents the optimized total thermal conductance of heat exchangers with different heat demand Q2 . As Q2 increases, the total thermal conductance increases as well, which means that the heat exchangers with larger area can deliver more heat. Fig. 5(b) shows the optimized temperature variation of the supply water versus Q2 . When Q2 increases, the optimized temperature of the supply water also increases. As the temperature of supply water increases, the temperature difference between supply water and indoor air increases which means more heat can be delivered to users. If we use traditional method to design the DHN with the temperature of the supply water given as 313 K, the optimization results are shown in Fig. 5(a), where (kA)tr represents the designed thermal conductance by traditional method. From Fig. 5(a), comparing the results of two methods, it is obvious that the optimized thermal conductance of heat exchangers in DHN obtained by entransy dissipation-based method obtains is less than that by the traditional method. 4.2. Optimization of a 2-column-4-user DHN The temperatures of supply and return waters T and T in the primary heating network in Fig. 1 are 393 K and 333 K, respectively. The users 11 and 12 in B1 have the same indoor air temperature and heat demand, 293 K and 15 kW, while the user 21 and 22 in B2 have the same indoor air temperature, 288 K, but different heat demands, 18 kW and 12 kW, respectively. In this case, two typical optimization problems come out. First, when the total thermal conductance of all heat exchangers keeps constant, 3000 W/K, Table 2 lists the optimized thermal conductance and heat capacity rate for water of

Fig. 5. Optimization results with the variation of Q2 . (a) Total thermal conductance variation versus Q2 . (b) The temperature of supply water variation versus Q2 .

each heat exchanger. Second, when the total heat capacity rate is a constant, 2000 W/K, the optimized results are listed in Table 3. Through the entransy dissipation-based optimization method, we get the optimal configuration of the thermal conductance of heat exchangers, i.e. (kA)0 , (kA)11 , (kA)12 , (kA)21 and (kA)22 , and the heat capacity rates of the waters flowing into each building, G1 and G2 . Among all the design parameters, there are two independent parameters during design processes. If G1 is exactly fixed with the optimized value in Table 3, varying (kA)11 yields the variation of the total thermal conductance of heat exchangers in the DHN. The total thermal conductance variation versus (kA)11 is showed and represented by G1,ED in Fig. 6, where it is quite obvious that the total thermal conductance of heat exchangers reaches its minimum when (kA)11 equals to the optimized value in Table 3, (kA)11 = 389.7 W/K. On the other hand, if (kA)0 is exactly fixed with the optimized value in Table 3, the total thermal conductance variation versus G1 is showed in Fig. 7, where the total thermal conductance of heat exchangers reaches its minimum when G1 equals to the optimized value in Table 3, G1 = 1051.5 W/K. In actual design processes, it is almost impossible to fix either the thermal conductance or the heat capacity rates exactly with the optimal values. If G1 is respectively fixed with the values deviating ±10% from the optimized one listed in Table 3, Fig. 6 gives the total thermal conductance of heat exchangers variation versus (kA)11 , where the minimums of the total thermal conductance in both situations are apparently larger than the value obtained by

58

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

Table 2 The optimized result with total thermal conductance of all heat exchangers of 3000 W/K. Parameters (W/K) Optimized results

(kA)0 1697.0

(kA)11 236.0

(kA)12 285.1

(kA)21 319.3

(kA)22 462.6

G1 1374.8

G2 491.5

(kA)12 625.9

(kA)21 440.7

(kA)22 475.5

G1 1051.5

G2 948.5

Table 3 The optimized result with the total heat capacity rate of 2000 W/K. Parameters (W/K) Optimized results

(kA)0 1828.3

(kA)11 389.7

Fig. 6. The total thermal conductance of heat exchangers variation versus (kA)11 with fixed G1 .

entransy dissipation-based method. Similarly, the other two situations where (kA)0 is fixed with values deviating ±5% from the optimized value in Table 3 are calculated, which are showed in Fig. 7. In Fig. 7, the minimums of the total thermal conductance in both situations are apparently larger than the value obtained by entransy dissipation-based method. These calculations and comparisons prove the superior of new optimization method than other traditional methods. If the total heat demand in B1 is equivalent to that in B2 , and both users in B1 have the same heat demand as 10 kW, the DHN is optimized when the ratio of the heat demands a = Q21 /(Q21 + Q22 ) varies from 0.1 to 0.9. The total heat capacity rate of the water in the DHN keeps the same, 2000 W/K. Fig. 8 shows the optimization results with the variation of a. In this case, the total heat demand

Fig. 7. The total thermal conductance of heat exchangers variation versus G1 with fixed (kA)0 .

of the users in each building is given, so it is generally believed that the variations of heat demands among users in one building can only cause the change of the area of the radiators in this building but have no influence on the radiators in the other buildings. However, Fig. 8(a) shows the variations of the thermal conductance of radiators, (kA)1 and (kA)2 , in B1 and B2 , respectively, versus the ratio of the heat demands a. As the ratio of the heat demands, a, increasing, the total thermal conductance of the radiators in B2 decreases, while that in B1 increases. The variation of the thermal conductance for the radiators in B1 and B2 proves that the thermal conductance of radiators in the DHN has to be adjusted when the heat distribution in one building changes. On the other hand, it is always believed that the total thermal conductance of all heat exchangers will keep constant as long as the total heat demand of all users in the DHN is given. However, Fig. 8(b) shows that the total thermal conductance of all heat exchangers in the DHN should change slightly when the heat ratio of heat demands in B2 varies, which proves that the variation of the heat demands among users of one building in the DHN leads to the change of the total thermal conductance for all the heat exchangers in the DHN. Fig. 9 shows the ratio between the optimized heat capacity rates of the waters in B1 and B2 versus the ratio of heat demands, a. When

Fig. 8. Optimization results with the variation of a. (a) Total thermal conductance of the radiators in B1 and B2 variation versus a. (b) Total thermal conductance in DHN variation versus a.

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

59

Fig. 11. The thermal conductance variation versus Q2 /Q1 . Fig. 9. The variation of the ratio between G2 and G1 versus a.

the ratio of the heat demands, a, increases, the ratio between the optimized heat capacity rates of the waters in B1 and B2 decreases, which means the heat demand in one building does have influence on the distribution for the heat capacity rates of the supply waters among all the buildings in the DHN. For the same DHN, if the heat demands of all users are the same as 15 kW and the temperature T2 increases from 288 K to 300 K while T1 keeps constant, 293 K, the optimization results are showed in Fig. 10. Fig. 10(a) shows the variation of the thermal conductance of the substation heat exchanger and the radiators in B1 and B2 , where (kA)0 , (kA)1 and (kA)2 represent the optimized thermal conductance of the substation heat exchanger and the radiators in

B1 and B2 , respectively. As T2 increases, both (kA)0 and (kA)2 get larger whereas (kA)1 decreases. It indicates that the increase of the indoor air temperature in a building will lead to the need for larger radiators in this building and a larger substation heat exchanger in the DHN to deliver the same heat. Fig. 10(b) shows the variation of the ratio between the heat capacity rates of water in B1 and B2 versus T2 . As T2 increases, the ratio of heat capacity rates of the waters in B1 and B2 , G2 /G1 , increases, which means the distribution of heat capacity rates of the supply waters in the DHN can also be affected by the heat demand of any user in the DHN. If Q11 = Q12 , Q21 = Q22 and Q1 keeps 10 kW, Fig. 11 shows the variation of (kA)0 , (kA)1 and (kA)2 due to the increase of Q2 /Q1 . As Q2 /Q1 increases, both (kA)0 and (kA)2 get larger while (kA)1 decreases. It indicates that in this situation, the increase of heat demand for the users in one building causes the increase of the thermal conductance of the radiators in this building and the substation heat exchanger in the DHN. The decrease of (kA)1 indicates that although heat demand of a building keep constant, the optimized thermal conductance of the radiators in this building should vary as well. 5. Conclusion

Fig. 10. Optimization results with the variation of T2 . (a) The thermal conductance variation versus T2 . (b) The ratio G2 /G1 variation versus T2 .

A theoretically global optimization method for DHNs is developed in this paper based on the physical quantities, i.e. entransy, entransy dissipation and entransy dissipation-based thermal resistance. Three typical irreversible heat transfer processes in DHNs are analyzed and measured by the entransy dissipation and then the total entransy dissipation of DHNs is calculated by two different ways. One is about structural and operating parameters in DHNs and the other is about the users’ demands in DHNs. Therefore, by the different expressions of entransy dissipation, the mathematical relation between all the structural and operating parameters and users’ demands in DHNs is established theoretically. Compared with traditional methods, the entransy dissipationbased optimization method can theoretically achieve the global optimization of DHN by directly establishing the relation between all the unknown parameters and users’ demands. Based on the relation, two optimization problems, i.e. minimizing total thermal conductance of heat exchangers with total energy consumption of pumps given, minimizing total energy consumption of pumps with total thermal conductance of heat exchangers given, can be mathematically converted into conditional extremum problems which offer optimization equations. By solving the equations, the optimized configuration of all the structural and operating parameters in DHN can be optimized essentially. This method is proved valid and convenient through the optimization of two simple DHNs. The optimization results indicate

60

Y.-C. Xu, Q. Chen / Energy and Buildings 48 (2012) 50–60

that the entransy dissipation-based optimization method can lead to optimization results better than the traditional methods. Most importantly, it shows that the variation of the heat distribution or the temperature of indoor air in any building in DHN can change the optimized thermal conductance of heat exchangers and the water distribution among all the buildings in the whole DHN, which denies the generally believed opinions. Acknowledgment The present work is supported by the National Natural Science Foundation of China (Grant No. 51006060). References [1] L.G. Fishbone, H. Abilock, MARKAL, a linear-programming model for energysystems analysis – technical description of the BNL version, International Journal of Energy Research 5 (4) (1981) 353–375. [2] M. Kleemann, D. Wilde, Intertemporal capacity expansion models, Energy 15 (7–8) (1990) 549–560. [3] A.S. Kydes, Flow models, Energy 15 (7–8) (1990) 561–571. [4] D. Henning, MODEST – an energy-system optimisation model applicable to local utilities and countries, Energy 22 (12) (1997) 1135–1150. [5] M. Bojic, N. Trifunovic, Linear programming optimization of heat distribution in a district-heating system by valve adjustments and substation retrofit, Building and Environment 35 (2) (2000) 151–159. [6] M. Bojic, N. Trifunovic, S.I. Gustafsson, Mixed 0-1 sequential linear programming optimization of heat distribution in a district-heating system, Energy and Buildings 32 (3) (2000) 309–317. [7] J. Soderman, F. Pettersson, Structural and operational optimisation of distributed energy systems, Applied Thermal Engineering 26 (13) (2006) 1400–1408. [8] D.K. Baker, S.A. Sherif, Heat transfer optimization of a district heating system using search methods, International Journal of Energy Research 21 (3) (1997) 233–252. [9] G. Cammarata, A. Fichera, L. Marletta, Using genetic algorithms and the exergonomic approach to optimize district heating networks, Journal of Energy Resources Technology – Transactions of the ASME 120 (3) (1998) 241–246. [10] J.A. Wright, H.A. Loosemore, R. Farmani, Optimization of building thermal design and control by multi-criterion genetic algorithm, Energy and Buildings 34 (9) (2002) 959–972. [11] Z.Y. Guo, X.G. Cheng, Z.Z. Xia, Least dissipation principle of heat transport potential capacity and its application in heat conduction optimization, Chinese Science Bulletin 48 (4) (2003) 406–410.

[12] Z.Y. Guo, H.Y. Zhu, X.G. Liang, Entransy – a physical quantity describing heat transfer ability, International Journal of Heat and Mass Transfer 50 (13–14) (2007) 2545–2556. [13] Q. Chen, Z.Y. Guo, Entransy theory and its application to heat transfer analyses in porous media, International Journal of Nonlinear Sciences and Numerical Simulation 11 (1) (2010) 11–22. [14] Q. Chen, H.Y. Zhu, N. Pan, Z.Y. Guo, An alternative criterion in heat transfer optimization, Proceedings of the Royal Society A: Mathematical Physical and Engineering Sciences 467 (2128) (2011) 1012–1028. [15] J.A. Meng, X.G. Liang, Z.X. Li, Field synergy optimization and enhanced heat transfer by multi-longitudinal vortexes flow in tube, International Journal of Heat and Mass Transfer 48 (16) (2005) 3331–3337. [16] Q. Chen, J.X. Ren, J.A. Meng, Field synergy equation for turbulent heat transfer and its application, International Journal of Heat and Mass Transfer 50 (25–26) (2007) 5334–5339. [17] Q. Chen, J.X. Ren, Generalized thermal resistance for convective heat transfer and its relation to entransy dissipation, Chinese Science Bulletin 53 (23) (2008) 3753–3761. [18] Q. Chen, M.R. Wang, N. Pan, Z.Y. Guo, Optimization principles for convective heat transfer, Energy 34 (9) (2009) 1199–1206. [19] J. Wu, X.G. Liang, Application of entransy dissipation extremum principle in radiative heat transfer optimization, Science in China Series E: Technological Sciences 51 (8) (2008) 1306–1314. [20] X.T. Cheng, X.G. Liang, Entransy flux of thermal radiation and its application to enclosures with opaque surfaces, International Journal of Heat and Mass Transfer 54 (1–3) (2011) 269–278. [21] J.A. Meng, X.G. Liang, Z.J. Chen, Z.X. Li, Experimental study on convective heat transfer in alternating elliptical axis tubes, Experimental Thermal and Fluid Science 29 (4) (2005) 457–465. [22] X.W. Li, J.A. Meng, Z.Y. Guo, Turbulent flow and heat transfer in discrete double inclined ribs tube, International Journal of Heat and Mass Transfer 52 (3–4) (2009) 962–970. [23] L. Chen, Q. Chen, Z. Li, Z.Y. Guo, Optimization for a heat exchanger couple based on the minimum thermal resistance principle, International Journal of Heat and Mass Transfer 52 (21–22) (2009) 4778–4784. [24] X.B. Liu, J.A. Meng, Z.Y. Guo, Entropy generation extremum and entransy dissipation extremum for heat exchanger optimization, Chinese Science Bulletin 54 (6) (2009) 943–947. [25] Q. Chen, K.D. Yang, M.R. Wang, N. Pan, Z.Y. Guo, A new approach to analysis and optimization of evaporative cooling system I: theory, Energy 35 (6) (2010) 2448–2454. [26] Q. Chen, N. Pan, Z.Y. Guo, A new approach to analysis and optimization of evaporative cooling system II: applications, Energy 36 (5) (2011) 2890–2898. [27] E. Abel, Low-energy buildings, Energy and Buildings 21 (3) (1994) 169–174. [28] D. Kut, R.M.E. Dianmant, District Heating and Cooling for Energy Conservation, The Architectural Press, London, 1981. [29] I.E. Agency, District Heating and Combined Heat and Power Systems, OECD Publications and Information Center, Washington, DC, 1983.