Paper No. TP-012

DYNAMIC MODEL FOR SMALL DISTRICT HEATING NETWORKS IN RURAL AREAS Georg K. Bestrzynski(a), Amru Rizal Razani(b, c), Holger Janßen(b), Andrea Luke(d), Stephan Scholl(a) (a)

Institute for Chemical and Thermal Process Engineering, TU Braunschweig, 38106 Braunschweig (b) Institute for Energy and Climate Protection, Hochschule Hannover, 30459 Hannover (c) Fernwärme-Forschungsinstitut, 30966 Hannover (d) Institute of Technical Thermodynamics Kassel, Universität Kassel, 34125 Kassel Germany

ABSTRACT Wide-spread static simulation tools for urban district heating networks are not sufficient for smaller applications, such as in rural areas. Thus, designing these networks is very precarious for potential heatsuppliers. In order to match the requirements of grid operators, customers and contractors, dynamic effects are taken into account, while the system is operating. Therefore, the resulting tool has to consider storage effects of the heat-supply-network, as well as interactions between its combined parts, such as heat-sources, heat-customers, heat-storages and (again) the heat-supply-network. A concept for a dynamic thermo-hydraulic simulation-tool matching these requirements is presented. The calculation of the nodal heads is the first step towards a thorough computation of the network’s state. An extensive example for the nodal head calculation is given. Subsequently, a thermal computation is done and combined with the results from the first step. Finally, the thermo-hydraulic state of the network is calculated.

1. INTRODUCTION All around the world the consequences of the progressing climate change become more and more relevant for our daily life. Meanwhile, the economical, ecological und humanitarian extents of these climatic changes are unforeseeable. Thus, a global reduction of greenhouse gas emissions is without any alternative. Against this background, raising the climate efficiency of fossil and regenerative fuel inputs in order to meet the electrical power and heat demands of the customer is imperative. Meanwhile, the effectiveness of district heating networks is about 80 to 90%, see Rosen et al. (2005) and Shipley et al. (2008). Hence, this common and proved technology meets a key factor for a sustainable concept of energy supply. Unfortunately, heat sales in urban areas are declining despite of supreme efforts on network extensions by its operators. This is clearly the result of progressing heat insulation acts, see Wulf et al. (2012) and Yang (2013). Consequently, a completely new field of district heating network application has to be identified in order to increase the capacities of this efficient technology. Therefore, small scale CHP capacities, such as (biogas driven) combustion units or micro gas turbines, must be installed in suburban (rural) areas in order to develop this technology. This serves needs of net stability, reduces electrical grid expansions and will raise climate efficiency, see Clausen and Kahle (2012) and Yang (2013). Inherent for the attached small scale district heating networks are small cumulative heat demands and/or few attached customers, differentiating them clearly from (big) urban networks. Resulting, interactions between heat sources, heat distribution network and heat customers jeopardize an untroubled and sufficient supply. Questions concerning network structure as well as operational strategies cannot be answered yet. Thus, a flexible tool for computing and visualizing the hydraulic and thermal state of arbitrary district heating networks is compiled within the research program on “Possibilities and boundaries of district heating systems in rural areas taking into account renewable heat sources” in order to examine these topics. 4th IIR Conference on Thermophysical Properties and Transfer Processes of Refrigerants, Delft, The Netherlands, 2013 1

Paper No. TP-012

2. THERMO-HYDRAULIC MODEL OF DISTRICT HEATING NETWORKS After converting a real-life district heating network into a graph, the thermo-hydraulic computation may be realized within two consecutive steps, see Deistel (2010). The first (central) step is the hydraulic modeling of the mass-flows within the network itself, depending on the boundary conditions, such as head sources and heat demands. This step provides all currents flowing through the pipes of the network. Simultaneously, this is the foundation for computing the thermal losses within the network, see Sacoph (1992). 2.1 Hydraulic Model The computation of the hydraulic state of the district heating network has been treated within several papers and dissertations so far, for example Valdimarsson (1993), Köcher (2000), Wernstedt (2005) or Ben Hassine and Eicker (2011). Nevertheless, an extensive and clear example of a stepwise description and elaboration of the mathematical key operations is not given in literature. Unfortunately, the contents of these references need a lot of adept interpretation, thus hampering access to the topic. Therefore, an illustrative example of a hydraulic calculation will be given in the following discussion. 2.1.1 Conversion of the examined network into a graph, mathematical notations and key operations In order to map and edit any real life heating network mathematically, a conversion of this network into a graph must be done (simplified hydraulic scheme see Figure 1a). Main assumptions for this conversion are: 



The relevant geometrical dimensions of the pre- and backflow are similar. Therefore the consideration of just one part of the network (pre-flow) is sufficient, see Figure 1b. Pipes are linear resistances for the mass-flow and represented by branches.

According to the resulting graph of the network, the Connectivity- (Incidence-) Matrix A, the Cutset-Matrix D and the Loop-Matrix B have to be defined/calculated. The dimension of the Connectivity-Matrix A is N x M, where N is the number of nodes (Arabic numbers: 1, 2, 3…) and M is the number of branches plus boundary elements (Roman numbers: I, II, III …). The entries anm within the Connectivity-Matrix A are: anm = 1 if the branch (pipe) m starts at node n, anm = -1 if the branch (pipe) m ends at node n or anm = 0 if both criterions above are not valid. The Connectivity-Matrix is split up into two subparts, the tree (AT) and the link part (AL) A=[AT|AL]. Examples of the division of the two subparts are given in Figure 2. No loops are found within the tree part of the matrix AT, neither over the ground (backflow), nor within the considered part of the network itself. Therefore, just two out of the three branches of the network shown in Figure 1b as well as one boundary condition form the tree part of the network.

Figure 1: a) Simplified hydraulic scheme of the district heating network b) Graph of the pre-flow network 4th IIR Conference on Thermophysical Properties and Transfer Processes of Refrigerants, Delft, The Netherlands, 2013 2

Paper No. TP-012

On the other hand, the elements of the link part AL either form a loop within the pipe structure itself (column 4 of each matrix in Figure 2) or connect the tree of the network to the ground (back-flow) once again (columns 5 & 6 of each matrix in Figure 2). A selection of graphical interpretations of Connectivity-Matrices A resulting from different interpretations of the graph itself are given in Figure 3 for the matrices Aa) & Ab). The next central step towards a hydraulic calculation is the assembling of the fundamental Cutset- Matrix D. This matrix is partitioned into two subparts D=[I|DL] again, while its dimensions are N x M. Thus, the dimensions of the link part of the Cutset-Matrix DL equals N x (M-N), where (M-N) is the number of link elements (link pipes and boundary elements) in the network. The special characteristics of this matrix are:  

Every line of the matrix D splits the whole graph into two separated parts. The separation is done by cutting just one tree element as well as several link elements of the graph (link pipes and boundary elements). Therefore, the number of fundamental cutsets equals the number of tree elements (e.g. no. of nodes).

The single entries of the Cutset-Matrix D are: dnm = 1 if the branch (pipe) m is a member of the cutset n having the same direction, dnm= -1 if the branch (pipe) m is a member of the cutset n having the opposite direction or dnm = 0 if the branch (pipe) is not a member of the cutset. It can be calculated by rearranging the Connectivty-Matrix A according to Eq. (1): (1) Occurring differences of the Cutset-Matrices Da) & Db), calculated by Eq. (1), result from different interpretations of the connectivity matrix A. Nevertheless, independent of differences in interpretations, there are general parallels. Thus, the first cutsets separate the whole network completely from the ground, whereas the other cutsets give a balance of mass-flows for each node adjacent to the link pipes (first columns of AL). The calculated Cutset-Matrices Da) & Db) are opposed to illustrations of every cutset in Figure 4. The algebraic sign of each entry becomes descriptive by considering the entries in the identity part I of the matrix first. Choosing the ground as reference system for the first fundamental cutsets, this identity entry demonstrates, that a mass-flow into the ground is assumed positive. Therefore, the entries in the fifth and sixth column are positive as well. On the other hand, by choosing the network part, all mass-flows leaving pipe structure are assumed positive, which leads to the same results concerning the algebraic sign.

Figure 2: Examples of dividing the Connectivity- (Incidence-) Matrix A into its subparts

Figure 3: Illustration of the Connectivity-Matrices Aa) & Ab) 4th IIR Conference on Thermophysical Properties and Transfer Processes of Refrigerants, Delft, The Netherlands, 2013 3

Paper No. TP-012

The final step of preparing the network in order to start a hydraulic calculation is the definition of the LoopMatrix B=[BT|I], which again can be obtained from the Cutset-Matrix D and therefore from the Connectivity-Matrix A following Eq. (2): (2) The dimensions of the tree part of the Loop-Matrix BT equals (M-N) x N and its entries are: bnm = 1 if the branch (pipe) m is a member of the loop n having the same direction, bnm = -1 if the branch (pipe) m is a member of the loop n having the opposite direction or bnm = 0 if the branch (pipe) is not a member of the cutset. The calculated Loop-Matrices Ba) & Bb) and illustrations of every fundamental loop represented by each line are given in Figure 5. The single loops are formed by one link element, considered within the identity part I of the Loop Matrix B, together with several tree elements, considered within BT. Furthermore, the first lines of the Loop-Matrices represent “internal” loops of the network without applying any boundary element (the ground), whereas the subsequent lines incorporate the ground. Thus, for our examples, the first lines give “internal” loops, while the other fundamental loops apply one link boundary element each. Meanwhile, the interpretation of the direction/algebraic sign is clearly connected to the originally defined direction of the branches (pipes). For each branch being a member of the loop, the entries are positive as long as the direction of the pipe is the same as the direction of the considered loop defined in the identity part I (e.g. link element) of the matrix. Otherwise the entry is negative, cf. Figure 5.

Figure 4: Resulting Cutset Matrices Da) & Db) opposed to their graphical interpretations

Figure 5: Loop Matrices Ba) & Bb) opposed to their graphical interpretations

2.1.2 The Laws of Kirchhoff Linear systems of equations are the result of combining Kirchhoff’s Current and Voltage Laws with the defined Connectivity-Matrix A, the calculated Cutset-Matrix D and the Loop-Matrix B, see Weißgerber (2009). Intuitively the correctness of the Current Law for the whole network is clear, as the sum of entering and leaving mass-flows for the whole network must be equal to zero (Eq. 3a & 3b). Furthermore, by adding another virtual node for the ground, which is not shown in any Figures above, the Current Law tells us that the sum of mass-flows entering and leaving one single node must be zero. ⃗⃗

⃗⃗

(3a; 3b)

Finally, the Voltage Law is applied for the calculated loop matrix B. This matrix represents all fundamental loops of the graph taking the ground into account as well: ⃗

(4)

The dimensions of the vectors ( ⃗⃗ & ⃗ ) are equal to the number of branches within the network. The given linear systems however are coupled by the vectorised Darcy-Weisbach equation, see VDI-Gesellschaft (2010) for example: 4th IIR Conference on Thermophysical Properties and Transfer Processes of Refrigerants, Delft, The Netherlands, 2013 4

Paper No. TP-012

̇

̇

(5)

For each pipe i, its length Li, diameter di and mass-flow ̇ has to be defined. The modulus of the mass-flow considers that there might be negative pressure drops along a pipe, as we arbitrarily choose a flow direction within each pipe in the beginning, in order to convert the network into a graph. This is essential for the calculation of networks with a “non-ground” link (branch IV for the Connectivity-Matrix Aa)). Finally, the friction factor f has to be calculated as well, e.g. by the Darcy-Weißbach or Swamee-Jain correlation, see VDI-Gesellschaft (2010). As this essential parameter for the pressure drop cannot be correlated to any other common parameter characterizing a surface (e.g. surface roughness), a fixed value of the friction factor f is chosen for our given example. 2.1.3 Calculation of the hydraulic state An exemplary calculation of the hydraulic state of the network is done for the Connectivity-Matrix Aa). A system of equations results from Eq. (3b), by converting the heat demands of the customers into a mass-flow ̇ ̇ ). The conversion of the heat demand into a mass-flow boundary condition is done by applying ( ̇ the first law of thermodynamics, see again VDI-Gesellschaft (2010). ⃗⃗

̇

]

[

̇

̇ ̇

̇

̇ ̇

̇

̇

(6)

), the first & ̇ By inserting the boundary flow conditions at the customer’s site ( ̇ line of Eq. (6) is solved, giving the entering mass-flow into the system. The negative algebraic sign is obviously correct considering the original and arbitrarily chosen direction of ̇ in Figure 1b. Furthermore, ̇ ) within the network. just two more equations are given for three unknown mass-flows ( ̇ , ̇

) has to be done for each link of the network, giving a first mass-flow Thus, an initial guess ( ̇ distribution for the rest of the network. This also gives a pressure drop ΔhLinkDarcy related to the initially guessed mass-flow according to Eq. (5). Besides, as all mass-flows within the rest of the network are known now as well, the pressure drops within the remaining pipes of the network can be calculated with the same equation. In combination with a nodal head at the head/heat source (100m of hydrostatic pressure), this gives a pressure distribution for all nodes of the network. This gives a pressure drop along the link ΔhLinkNodal as well. Unfortunately, this pressure drop along the link will (most likely) not be the same as the previously calculated pressure drop ΔhLinkDarcy. In other words, the initial guess/previous value for the mass-flow is wrong and must be adjusted. Finally, this adjustment can be done by linearizing the mass-flow around the actual state of flow, see Eq. (7). ̇

̇ ̇

̇

̇

(7)

This linearization gives hints concerning quality and quantity (direction and extent) of adjustment. The single ̇ ̇ are calculated by using Eq. (5) again. The combination of the ̇ ̇ & pressures pressure-drop-differences ΔhLinkNodal - ΔhLinkNodal and the linearization, gives a fixed value for the mass-flow within the link for the next iteration according to Eq. (8), as long as the residual drops below a defined value: ̇

̇

(8)

The result of the several consecutive steps of iteration concluding the hydraulic calculations, is given in, Table 1. In order to simplify the model and to make it more illustrative and intuitively comprehensible, all pipes have the same geometrical dimensions (diameter di=53mm and length Li=1000m). Furthermore, the friction factor f, whose calculation itself is rather sophisticated, is a fixed and representative value of f=0.004. Thus, after four iterations, the correction of the mass-flow through the link pipe becomes negligible. 4th IIR Conference on Thermophysical Properties and Transfer Processes of Refrigerants, Delft, The Netherlands, 2013 5

Paper No. TP-012

̇

Table 1: Stepwise iteration of the hydraulic calculation for Iterations

2

3

Correction [kg/s] ̇

5

1

100

ΔhLinkNodal32

12,95 m

III

3

2

79.76

ΔhLinkDarcy32

0. 81 m

IV

1 (guess)

3

92.71

Δhlin

14.57 m/(kg/s)

II

4.17

1

100

ΔhLinkNodal32

2.15 m

III

3.83

2

85.95

ΔhLinkDarcy32

2.72 m

IV

1.833

3

88.10

Δhlin

15.92 m/(kg/s)

II

4.20

1

100

ΔhLinkNodal32

2.61 m

III

3.79

2

85.71

ΔhLinkDarcy32

2.6172 m

IV

1.798

3

88.32

Δhlin

15.86 m/(kg/s)

II 1

̇

̇

0.833

1.833

̇

̇

-0.035

̇

1.798 ̇

-6.36E-5

1.798

2.2 Coupling of the thermal model Based on the results of the hydraulic calculation, a simple thermal model is applied in order to illustrate the thermal state of the network. For this purpose, the model presented by Jarfelt and Persson (2005) seems to be most promising, see the visualization in Figure 6. Furthermore, this model seems to give realistic values for heat losses of singular pipes, comparing these to data given by manufacturers of twin pipes. Due to the brevity of this paper, just the fundamental assumptions for the calculations according to this model are given, without presenting the underlying system of equations at all.     

The generally relevant geometry of a twin pipe system, cf. Figure 6, is considered within the model. The outdoor temperature influences the heat loss of the system, which is not deep under the ground. The internal heat transfer from the fluid to the insulation is ideal, e.g. independent of the mass-flow. The thermal resistance of the pipe material (steel) is neglected. The backflow influences the pre-flow by absorbing thermal energy. This virtual heat transfer represents the asymmetric temperature distribution around the twin pipe system.

Especially the last assumption seems to be very well customized concerning the aim of the calculations reducing the computational efforts to a minimum. The results of the thermal calculation as well as its boundary conditions are given in Table 2. The temperature calculated for the 2nd node is a mixture temperature of the mass- flows coming from node one and node three. It is calculated by assuming a constant heat capacity cp of the medium (water). The temperature drop inside each pipe is calculated by an exponential equation, see Glück (1985). However, discharging the calculated heat loss at the end of each pipe at once, leads to tolerable deviations (