PhD Dejan Brkic
email: dejanrgf@tesla.rcub.bg.ac.rs

PhD Dejan Brkić

Iterative methods for looped network pipeline calculation

Keywords: Pipeline networks, Waterworks, Natural gas distribution, Calculation methods, Flow, Hydraulic systems Pipeline networks, Water distribution system, Natural gas distribution system, Calculation methods, Flow, Hydraulic pipeline systems

Since the value of the hydraulic resistance depends on flow rate, problem of flow distribution per pipes in a gas or water distributive looped pipelines has to be solved using iterative procedure. A number of iterative methods for determining of hydraulic solution of pipeline networks, such as, Hardy Cross, Modified Hardy Cross, Node-Loop method, Modified Node method and M.M. Andrijašev method are shown in this paper. Convergence properties are compared and discussed using a simple network with three loops. In a municipal gas pipeline, natural gas can be treated as incompressible fluid. Even under this circumstance, calculation of water pipelines cannot be literary copied and applied for calculation of gas pipelines. Some differences in calculations of networks for distribution of these two fluids, i.e. water apropos natural gas are also noted.

Introduction

A number of iterative methods for determining the hydraulic solution of water and natural gas pipeline networks which take ring-like form, such as, Hardy Cross, Modified Hardy Cross (including Andrijašev method), Node-Loop method and Modified Node method are compared 2 in this paper. All presented methods assume equilibrium among pressure and friction forces in steady and incompressible flow. As a result, they cannot be successfully used in unsteady and compressible flow calculations with large pressure drops where inertia force is important.

Problem of flow in pipes and open conduits was one which had been of considerable interest to engineers for nearly 250 years. Even today, this problem is not solved definitively. The difficulty to solve the turbulent flow problems lies in the fact that the friction factor is a complex function of relative surface roughness and Reynolds number. Precisely, hydraulic resistance depends on flow rate and hence flow problem in hydraulic networks has to be solved iteratively. Similar situation is with electric current when diode is in circuit. With common resistors in electrical circuits where the electric resistances do not depend on the value of electric current in the conduit, problem is linear and no iterative procedure has to be used. So problem of flow through a single tube is already complex. Despite of it, very efficient procedures are available for solution of flow problem in a complex pipeline such as looped pipeline like waterworks or natural gas distribution network is. Most of the shown methods are based on solution of the loop equations while the node equations are used only as control of accuracy. Node-Loop method is also based on solution of the loop equations but the node equations are also used in calculation and not only for control purposes. Node method and here shown Modified Node method is based on solution of the node equations while the loop equations is used for control purposes.

In this paper, speed of convergence of presented methods are compared and discussed. This is done for one simple network with three loops, both, for water and gas distribution networks. Similar example can be done for the air ventilation systems in the buildings or mines. For ventilation problem readers can consult paper of [START_REF] Aynsley | A resistance approach to analysis of natural ventilation airflow networks[END_REF], [START_REF] Mathews | A numerical optimization procedure for complex pipe and duct network design[END_REF], [START_REF] Wang | Computer solution of three-dimensional mine ventilation networks with multiple fans and natural ventilation[END_REF], etc. While the disposal of water or gas distribution system is in single plane with certain elevation of pipes, ventilation network is almost always spatial [START_REF] Brkić | An Improvement of Hardy Cross method applied on looped spatial natural gas distribution networks[END_REF].

Hydraulics frictions and flow rates in pipes

Each pipe is connected to two nodes at its ends. In a pipe network system, pipes are the channels used to convey fluid from one location to another. Pipes are sometimes referred to as tubes or conduits, and as a part of network to as lines, edges, arcs or branches (in the graph theory, special kinds of branches are referred to as trees). The physical characteristics of a pipe include the length, inside diameter, roughness, and minor loss coefficient. When fluid is conveyed through the pipe, hydraulic energy is lost due to the friction between the moving fluid and the stationary pipe surface. This friction loss is a major energy loss in pipe flow [START_REF] Farshad | New developments in surface roughness measurements, characterization, and modeling fluid flow in pipe[END_REF][START_REF] Kumar | Parameter estimation in water distribution networks[END_REF].

Various equations were proposed to determinate the head losses due to friction, including the Darcy-Weisbach, Fanning, Chezy, Manning, Hazen-Williams and Scobey formulas. These equations relate the friction losses to physical characteristics of the pipe and various flow parameters. The Darcy-Weisbach formula for calculating of friction loss is more accurate than the Hazen-Williams. Beside this, the Hazen-Williams relation is only for water flow.

Gas flow rates and pressure drops in pipes

When a gas is forced to flow through pipes it expands to a lower pressure and changes its density. Flow-rate, i.e. pressure drop equations for condition in gas distribution networks assumes a constant density of a fluid within the pipes. This assumption applies only to incompressible, i.e. for liquids flows such as in water distribution systems for municipalities (or any other liquid, like crude oil, etc.). For the small pressure drops in typical gas distribution networks, gas density can be treated as constant, which means that gas can be treated as incompressible fluid [START_REF] Pretorius | A flow network formulation for compressible and incompressible flow[END_REF], but not as liquid flow. Assumption of gas incompressibility means that it is compressed and forced to convey through pipes, but inside the pipeline system pressure drop of already compressed gas is small and hence further changes in gas density can be neglected. Fact is that gas is actually compressed and hence that volume of gas is decreased and then such compressed volume of gas is conveying with constant density through gas distribution pipeline. So, mass of gas is constant, but volume is decreased while gas density is according to this, increased. Operate pressure for distribution gas network is 4·10 5 Pa abs i.e. 3·10 5 Pa gauge and accordingly volume of gas is decreased four times compared to volume of gas at normal conditions. Inner surface of polyethylene pipes which are almost always used in gas distribution networks are practically smooth and hence flow regime in the typical network is hydraulically "smooth" [START_REF] Sukharev | Investigation of the hydraulic resistance in polyethylene pipelines[END_REF]. For this regime is suitable Renouard"s equation adjusted for natural gas flow (1):

82 . 4 82 . 1 2 2 2 1 4810 in r n g D L Q p p C F         (1)
Regarding to Renouard"s formula has to be careful since it does not relate pressure drop but actually difference of the quadratic pressure at the input and the output of pipe. This means that C is not actually pressure drop in spite of the same unit of measurement, i.e. same unit is used for the pressure (Pa). For gas pipeline calculation is very useful fact that when C →0 this consecutive means that also C→0. Parameter C can be noted as pseudo-pressure drop.

Liquid flow in pipes

In Renoard"s equation adjusted for gas pipelines (1) friction factor is rearranged in the way to be expressed using other flow parameters and also using some thermodynamic properties of natural gas. Using formulation for Darcy friction factor in hydraulically smooth region Renouard suggests his equation for liquid flow (2):

            2 2 5 in 2 1 8 D L Q p p p F w (3)
Renouard"s equation ( 2) is based on power law. Liquid flow is better fit using some kind of logarithmic formula like Colebrook"s (4):

              in D 71 . 3 Re 51 . 2 log 2 1    (4)
Colebrook"s equation is suitable especially for flow through steel pipes [START_REF] Colebrook | Turbulent flow in pipes with particular reference to the transition region between the smooth and rough pipe laws[END_REF]. Some researchers adopt a modification of the Colebrook equation, using the 2.825 constant instead of 2.51 especially for gas flow calculation (Haaland 1983, Coelho and[START_REF] Coelho | Considerations about equations for steady state flow in natural gas pipelines[END_REF].

In the case of waterworks, pressures will also be expressed in Pa, not in meter (m) equivalents. In that way, clear comparison with gas distribution network can be done.

Pressures expressed in Pa can be very easily recalculated in heads expressed in meters (m) for water networks knowing water density. Example network is located in a single-flat area with no variation in elevation (otherwise the correction term must be used which is significant for water network).

Looped pipeline networks for distribution of fluids

A pipeline network is a collection of elements such as pipes, compressors, pumps, valves, regulators, tanks, and reservoirs interconnected in a specific way. In this article focus is on pipes. The behavior of the network is governed by the specific characteristics of the elements and how the elements are connected together. Our assumption is that pipes are connected in a smooth way, i.e. so called minor hydraulic loses are neglected. Including other elements different than pipes is a subject of sufficient diversity and complexity to merit a separate review.

The analysis of looped pipeline systems by formal algebraic procedures is very difficult if the systems are complex. Electrical models had been used in study of this problem in the time before advanced computer became available as background to support demandable numerical procedures [START_REF] Mah | Pipeline network design and synthesis[END_REF]. Here presented methods use successive corrections or better to say these methods require iterative procedure. The convergence properties of presented methods are the main subject of this article. All of such methods can be divided into two groups (1) Methods based on solution of the loop equations, and (2) Methods based on solution of the node equations.

Most of the methods used commonly in engineering practice belong to the group based on solution of the loop equations. Such method presented in this paper is Hardy Cross method [START_REF] Cross | Analysis of flow in networks of conduits or conductors[END_REF]. Method of Andrijašev is variation of Hardy Cross method (Andrijašev 1964) 1 .

1 There are some difficulties with citation of Russian methods. Both methods, Lobačev and Andrijašev, are actually from 1930"s but original papers from Soviet time are problem to be found. Many books in Russian language explain methods of Lobačev and Andrijašev but in these books reference list is not available.

Contemporary with Hardy Cross, soviet author V.G. Lobačev [START_REF] Latišenkov | Hydraulics. Gosstroizdat[END_REF] developed very similar method compared to original Hardy Cross method. Modified Hardy Cross method proposed [START_REF] Epp | Efficient code for steady flows in networks[END_REF] which considers entire system simultaneously is also sort of loop method. Node-Loop method proposed by [START_REF] Wood | Hydraulic network analysis using linear theory[END_REF] and later improved by [START_REF] Wood | Reliability of algorithms for pipe network analysis[END_REF] is combination of the loop and node oriented methods, but despite of its name is essentially belong to the group of loop methods. Node-Loop method is based on solution of the loop equations while the node equations are also involved in calculation. Only Node method proposed by [START_REF] Shamir | Water distribution systems analysis[END_REF] is real representative of node oriented methods. Node method uses idea of Hardy Cross but to solve node equations instead of loops ones. Improved Node method uses similar idea as [START_REF] Epp | Efficient code for steady flows in networks[END_REF] suggested for the improvement of original Hardy Cross method, but of course here applied to the Node method.

Example network with three loops is shown in the figure 1. Before calculation, maximal consumption for each node including one of more inlet nodes has to be determined. Pipe diameters and node inputs and outputs cannot be changed during the iterative procedure. Goal is to find final flow distribution for this pipeline system (simulation problem). pressure drops or pseudo-pressure drops in the case of node oriented methods. To have a better appreciation of the utility of these representations, first will be considered the laws that govern flow rates and pressure drops in a pipeline network. These are the counterparts to Kirchhoff"s laws for electrical circuits and accordingly for hydraulic networks, namely, (i) the algebraic sum of flows at each node must be zero; and (ii) the algebraic sum of pressure drops for water network (i.e. pseudo-pressure drop for gas network) around any cyclic path (loop) must be approximately zero. For a hydraulic network, fact that both laws are satisfied simultaneously means that calculation for steady state for simulation problem is finished.

Methods based on solution of the loop equations

Loops are sometimes referred also to as contours or paths. Note that contours and loops are not synonyms in the M.M. Andrijašev method.

For the loop oriented methods first Kirchhoff"s law must be satisfied for all nodes in all iterations. Second Kirchhoff"s law for each loop must be satisfied with acceptable tolerance at the end of the calculation. Or, in other words, final flows are these ones which values are not changed between two successive iterations (must be satisfied for flow in each pipe).

In a loop oriented methods, first initial flow pattern must be chosen to satisfied first Kirchhoff"s law. Endless number of flow combinations can satisfy this condition. Someone can conclude that it is crucial to start with good initial guesses. But, how does one obtain good initial guesses? And how sensitive are the methods to the initial guesses? Answer is that initial flow pattern does not have significant influence on convergence properties of observed method. These indicate that the choice of initial flows is not critical and need to be chosen only to satisfied first Kirchhoff"s law for all nodes [START_REF] Gay | The solution of pipe network problems[END_REF]. method are not actually flows Q but correction of flows ΔQ calculated for each loop. These corrections have to be added algebraically to flows from previous iteration for each pipe according to specific rules. A pipe common to two loops receives two corrections. The upper plus or minus sign (shown in tables 1-8) indicates direction of flow in that conduit in these two loops and is obtained from Q for previous iteration. The upper sign is the same as the sign in front of Q if the flow direction in each loop coincides with the assumed flow direction in the particular loop under consideration, and opposite if it does not. The lower sign is copied from the primary loop for this correction (sign from the loop where this correction is first, sign preceding the first iteration from adjacent contour for the conduit taken into consideration). The rules for sign of corrections ΔQ 2 are: (1). the algebraic operation for correction ΔQ 1 should be the opposite of its sign; i.e. add when the sign is minus. (2). the algebraic operation for corrections ΔQ 2 should be the opposite of their lower signs when their upper signs are the same as the sign in front of Q, and as indicated by their lower signs when their upper signs are opposite to the sign in front of Q. For details of sing of corrections consult paper of [START_REF] Brkić | An Improvement of Hardy Cross method applied on looped spatial natural gas distribution networks[END_REF] and article [START_REF] Corfield | Distribution design for increased demand[END_REF] from Gas Engineers Handbook.

Hardy Cross method for the pipeline calculation

Hardy Cross, American engineer, developed method in 1936, later named after him [START_REF] Cross | Analysis of flow in networks of conduits or conductors[END_REF]. Hardy Cross calculation for gas pipeline network is shown in table 1 and for water network in table 2. Only two iterations will be shown in details. Flow Q 1 calculated in first iteration become initial flow Q for second iteration. The plus or minus preceding the flow Q, indicates the direction of the conduit flow for the particular loop (Figure 2). A plus sign denotes clockwise flow in the conduit within the loop; a minus sign, counter-clockwise. A flow correction ΔQ 1 as shown in table 1 and 2 is computed for each loop. This correction must be subtracted algebraically from the assumed gas flow according to specific rules explained previously.

Correction for each loop ΔQ 1 is calculated using first derivative of pressure function C for pipes defined by Renouard"s equation for gas networks (5) and derivative of pressure drop Δp for pipes defined by Darcy-Weisbach for water networks (6) where flow Q is treated as variable:

    82 . 0 82 . 4 82 . 0 82 . 1 82 . 4 82 . 1 82 . 1 4810 82 . 1 4810 ' n in r n in r n Q R D L Q Q Q R Q D L Q Q C F                                (5)     Q R Q Q Q R Q Q Q p F                                  2 D L 16 D L 8 ' 5 2 2 5 2 2       (6)
Derivative for loops is calculated using assumed loop shown in the figure 2 with no reference to direction. For the loop I this derivative for gas network (7) and for water network (8) can be written:

          I r I I I I I D L Q D L Q D L Q Q Q C Q C Q C Q Q C F                           
 (7)           I I I I I I Q Q Q Q Q p Q p Q p Q Q p F                                  5 7 7 7 7 5 4 4 4 4 5 3 3 3 3 2 7 7 4 4 3 3 ' D L D L D L 16      (8)
In presented example loop I begins and ends in node II via pipes 3, 4, and 7.

Correction of flow (ΔQ 1 ) for each loop can be calculated (9):

' 1 F F Q    (9)
                                                                     III II I III II I III III II II I I F F F Q Q Q x Q Q Q Q Q F Q Q Q Q Q F Q Q Q Q F 8 7 6 5 5 4 2 1 7 4 3 , , , 0 0 0 , , , 0 0 0 , , (10)
Using numerical values from table 2 for the first iteration of calculation of water network matrix equation can be written (11):

                                  2070510 5040086 15928498 - 193795963 0 0 0 53910587 0 0 0 129435625 III II I Q Q Q x (11)
Then set of correction for water network calculation in the first iteration is [ΔQ I , ΔQ II , ΔQ III ] T =[-12306·10 -5 , 9349·10 -5 , 1068·10 -5 ] T which is identical as in table 2.

Modified Hardy Cross method for the pipeline calculation

Modified Hardy Cross method (somewhere called improved) is also known as the simultaneous loop adjustment method. As seen in figure 1, several loops have common pipes, so corrections to these loops will cause energy losses around more than one loop. In figure 1, pipe 4 belongs to two loops (loop I and II), pipe 7 to loop I and III, and finally pipe 5 to II and III. Modified Hardy Cross method is a sort of Newton-Raphson method used to solve unknown flow corrections taking into consideration whole system simultaneously. [START_REF] Epp | Efficient code for steady flows in networks[END_REF] gave idea for this approach. To increase efficiency of Hardy Cross method zeros from non-diagonal term in matrix equation (10) will be replaced to include influence of pipes mutual with adjacent loop (12).

                                                                                                                          
F F F Q Q Q x Q Q Q Q Q F Q Q F Q Q F Q Q F Q Q Q Q Q F Q Q F Q Q F Q Q F Q Q Q Q F 8 7 6 5 5 7 5 5 4 2 1 4 7 4 7 4 3 , , , , , , , , (12)
Presented matrix is symmetrical; for example (13):

    I III III I Q Q F Q Q F        ) ( ) ( 7 7 (13)
13 This is because pipe 7 is mutual for two adjacent loops (loop I and loop III). Non-diagonal terms have always opposite sign than diagonal. Spatial networks common for ventilation systems in buildings or mines are exceptions [START_REF] Brkić | An Improvement of Hardy Cross method applied on looped spatial natural gas distribution networks[END_REF].

Using numerical values from table 2 for the first iteration of calculation of water network matrix equation can be written ( 14):

                                  20705102 50400861 159284985 - 193795963 429677 - 70157780 - 429677 - 53910587 4978573 - 70157780 - 4978573 - 5 1294356252 III II I Q Q Q x (14)
First two iterations from our example using the modified Hardy Cross method is shown in table 3 for gas network and in table 4 for water network. Table 3. Calculation after the modified Hardy Cross method for example gas network Table 4. Calculation after the modified Hardy Cross method for example water network

Modified Andrijašev method for the pipeline calculation

Andrijašev method can be used in the formulation as in the original Hardy Cross method and as in Modified Hardy Cross method. Here it will be given in notation as improved method because this approach shows better convergence performance (example for two iteration of gas network calculation is shown in table 5 and for water network in table 6). It can be noted that some pipes in table 1-4 received only one correction per iteration (for example pipe 3 in contour I). This means that pipe 3 belongs only to one loop. In method of M.M. Andrijašev contours can be defined to include few loops. Thus, contours can be defined in other way and then each pipe in the network belongs to two contours (Figure 3). Loop is not synonym with contour in M.M. Andrijašev method as in Hardy Cross approach. Andrijašev"s contour will be marked with special sign ( O ). 3) starting and ending in node I via pipes 4, 5. 6, 8, 3, contour II O (green circuit in figure 3) starting and ending in referent node via pipes 1, 6, 8, 7, 4, 2, and finally contour III O (blue circuit in figure 3) starting and ending in referent node via pipes 1, 5, 7, 3, 2.

Matrix formulation of this method for example gas network can be written as (15):

                                                                                                                         o III o II o I o III o II o I o III o III o II o III o I o III o III o II o II o II o I o II o III o I o II o I o I o I F F F Q Q Q x Q Q Q Q Q Q F Q Q Q Q F Q Q Q F Q Q Q Q F Q Q Q Q Q Q Q F Q Q Q Q F Q Q Q F Q Q Q Q F Q Q Q Q Q Q F 7 5 3 2 1 7 2 1 5 3 7 2 1 8 7 6 4 2 1 8 6 4 5 3 8 6 4 8 6 5 4 3 , , , , , , , , , , , , , , , , , , , , , , , (15)
Here has to be very careful because non-diagonal terms are not always negative as in modified Hardy Cross method (16): which is used for calculation. Furthermore, calculation of single-flat pipeline networks using M.M. Andrijašev method is equally complex as calculation of spatial, multidimensional network, using Modified Hardy Cross method [START_REF] Brkić | An Improvement of Hardy Cross method applied on looped spatial natural gas distribution networks[END_REF].

                                 

Methods based on solution of the node equations

Nodes are sometimes also referred to as junctions, points or vertices. [START_REF] Shamir | Water distribution systems analysis[END_REF] introduced first node oriented method. For the node oriented methods second Kirchhoff"s law must be satisfied for each loop in all iterations. First Kirchhoff"s law for each node must be satisfied with acceptable tolerance at the end of the calculation. Here will be presented only one method. Improvement of presented method is done according to the same idea as used to improve original Hardy Cross method but here applied to the matrix solution for the node equations.

Node method

Pipe equations in previous text were expressed as Δp=F(Q) for waterworks, or C= 2 2 2 1 p p  =F(Q) for gas networks. These relations can be rewritten in form as Q=f(Δp) for waterworks or Q=f(C) for gas networks. After that, Renouard"s equation (1) can be rearranged (17):

  82 . 1 1 82 . 4 82 . 1 1 82 . 4 2 2 2 1 4810 4810 ) (                           r in r in g L D C L D p p C Q f   (17)
Similar, can be done for Darcy-Weisbach equation ( 18):

                        L 8 D L 8 D ) ( 2 5 in 2 5 in 2 1 p p p p Q f w (18)
In Node method, (C) for each pipe in a gas network and (Δp) for each pipe in a water network has to be assumed, not flows. These assumed functions of pressure must be chosen to satisfy second Kirchhoff"s law (Figure 4). First Kirchhoff"s law will be fulfilled with demanded tolerance at the end of calculation. For gas network, correction of C, noted as ΔC must be calculated, and for water network, correction of Δp, here noted as Δ Δp . Same algebraic rules as for loop oriented methods are valid. Correction ΔC or Δ Δp will be calculated using first derivative of Q=f(Δp) for waterworks and Q=f(C) for gas networks, where Δp, i.e. C is treated as variable (19):

' f f x    (19)
In previous equation x is Δp or C. Example network adjusted for node oriented calculation is shown in figure 4 (red letters for gas network and blue for water network).

Figure 4. Example of pipeline gas network with three loops adjusted for node oriented methods Calculation of the network using non-improved Node method is not shown here.

Modified Node method

Modified Node method is also referred to as the simultaneous node adjustment method. For the improvement of Node method, same idea as for the improvement of the Hardy Cross method is used. One node must be chosen to be omitted from the calculation shown in table 7 and 8, because linear independency among node equations is preserved in that way [START_REF] Mathews | A numerical optimization procedure for complex pipe and duct network design[END_REF]Köhler 1995, Mah and[START_REF] Mah | Pipeline network design and synthesis[END_REF]. More details about graph theory of networks will be explained in further text (section 3.4). Gas network calculation after Modified Node method is shown in table 7:

Table 7. Calculation after Modified Node method for example gas network with three loops Corrections in the table 7 and 8 are calculated using matrix equation ( 20):

                                                                                                                                        V IV III II I Q Q Q Q Q x x x x x x x x f x f x f x f x x x f x f x f x x x f x f x f x f x f x x x f x f x f x f x x x f
' 2 3 ' 2 4 ' 1 3 ' 1 4 3 2 ' 1 Q Q Q Q Q Q Q Q Q Q Q Q Q Q , 0 0 , , 0 0 0 , , 0 , , 0 0 , , (20)
When non-diagonal terms in the first matrix (21) are equalized with zero, results are equal as for Node method in basic form (19). For the first iteration, numerical values for the matrix equation are extracted from table 6 (21): For the node oriented calculation applied for water networks, pressure function (C) must be replaced with pressure drop (Δp). Further analogy is clear (Table 8).

                                                                          0.
Table 8. Calculation after Modified Node method for example water network with three loops

Node-Loop method

Conditionally, Node-Loop method can be sorted in the group of methods based on solution of the loop equations according to the previous discussion, but better solution is to be treated as combination of the loop and node oriented methods as its name unambiguously suggest.

Node-Loop method is also known as the flow adjustment method. [START_REF] Wood | Hydraulic network analysis using linear theory[END_REF] developed the flow adjustment method by coupling the loop equations with the node equations. [START_REF] Wood | Reliability of algorithms for pipe network analysis[END_REF] improved this method. Rather than solve for loop corrections, in this method, conservation of energy around a loop is written directly in the terms of the pipe flow rates. Final result after this method is not flow correction, but even better flow itself.

Pipeline network with three loops and six nodes is used as example in this article. Graph has X branches (pipes) and Y nodes where in our example, X=8 and Y=6. Graph with Y nodes (in our case 6) has Y-1 independent nodes (in our case 5) and X-Y+1 independent loops (in our case 3). Tree is a set of connected branches chosen to connect all nodes, but not to make any closed path (not to form a loop). Branches, which do not belong to a tree, are links (number of links are X-Y+1). Loops in the network are formed using pipes from tree and one more chosen among the link pipes. Number of the loops is determined by number of links. Network from example has six nodes and five independent nodes. One node can be omitted from calculation while no information on the topology in that way will be lost. Rows in the node matrix with all node included are not linearly independent. To obtain linear independence any row of the node matrix has to be omitted. For example, pipe 6 is between node IV and V, and reasonable assumption is that if node IV is output node for flow through pipe 6, then node V must be input node for flow through this pipe. In our example, node VI will be omitted. First Kirchhoff's law for the initial flow pattern shown in figure 2 can be written using set of equations ( 22

                          Q Q Q node Q Q Q node Q Q Q Q node Q Q Q Q node Q Q Q Q node Q Q Q Q node VI VI V V IV IV III III II II I I (22)
Node VI will be omitted from the node matrix to assure linear independency of the rows as shown in figure 4 (23):

                                 1 0 1 0 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 1 1 1 0 N (23)
First Kirchhoff's law must be fulfilled in all iterations for all nodes. Second Kirchhoff's law for the initial water flow pattern shown in the figure 1 can be written using set of equations ( 24):

III III II II I I F F F F F loop F F F F F loop F F F F loop              8 7 6 5 5 4 2 1 7 4 3 ~ (24)
Second Kirchhoff's law for the initial flow pattern shown in the figure 1 also can be noted in matrix form using loop matrix (25):

                        1 1 1 1 0 0 0 0 0 0 0 1 1 0 1 1 0 1 0 0 1 1 0 0 L (25)
Second Kirchhoff's law must be fulfilled for all loops at the end of calculation with demanded accuracy (i.e. F I →0, F II →0 and F III →0).

In Node-Loop method these two matrices become one with some modifications. The nodes and the loops equations shown in previous text here will be united in one coherent system by coupling these two set of equations. To introduce matrix calculation, the node-loop matrix

[NL], matrix of calculated flow in observed iteration [Q], and [V] matrix will be defined (26):

                                                               ' 8 8 ' 7 7 ' 6 6 ' 5 5 ' 5 5 ' 4 4 ' 2 2 ' 1 1 ' 7 7 ' 4 4 ' 3 3 F F F F F F F F F F F F F F F F F F F F F F F F F Q Q Q Q Q V III II I V IV III II I (26)
Sign minus preceding some of the flows Q in matrix [V] means that this particular Q is not consumption (sing minus represent inlet of fluid). Node-Loop matrix [NL] can be defined using node matrix, loop matrix and first derivative of Renouard"s function for gas pipes or of Darcy-Weisbach for water networks, as follows ( 27):

                                                                 ' 8 ' 7 ' 6 ' 5 ' 5 ' 4 ' 2 ' 1 ' 7 ' 4 ' 3 1 1 1 1 0 0 0 0 0 0 0 1 1 0 1 1 0 1 0 0 1 1 0 0 1 0 1 0 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 1 1 1 0 F F F F F F F F F F F NL (27)
Further, vector [Q] of the unknown flows can be calculated in the first iteration (28). Here will be used values from table 1 for gas network (29) and from table 2 for water network

      V x NL inv Q Q Q Q Q Q Q Q Q                             8
0 1 0 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 1 1 1 0 8 7 6 5 4 3 2 1 Q Q Q Q Q Q Q Q x (29)                                                                                     
1 0 1 0 0 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 1 1 1 0 8 7 6 5 4 3 2 1 Q Q Q Q Q Q Q Q x (30)
For gas network last three rows are calculated as follows using values from table 1, Excellent book for waterworks calculation by [START_REF] Boulos | Comprehensive water distribution systems analysis handbook for engineers and planners[END_REF] can be recommended for further reading. In this book, authors instead to omit one node in the node matrix to preserve linear independency of rows in this matrix introduce one pseudo-loop in loop matrix. This procedure is not practical because at least two nodes with equal pressure must be found in the network. This is not always possible. Further in that way the node-loop matrix has two additional rows which could be avoided. [START_REF] Mathews | A numerical optimization procedure for complex pipe and duct network design[END_REF] in his discussion use simplest way, i.e. they omit one row.

Comparison of solution techniques for looped piping networks

Final flows are unique after all presented methods, and will be listed in table 9, both for water and for gas network. Table 9. Final flows for network presented in this paper Each method has advantages and shortcomings. Convergence performances will be compared for all presented methods in figure 5. Note that Modified Node method cannot be compared literary because initial values cannot be equalized. In all other methods initial patterns are given in the form of flows and equalized, while in Node method initial pattern is in the form of pressures (better to say function of pressure). Best way to compare water and gas distribution network is to compare velocity of gas and water through pipes. For water, this can be done using (31):

    2 4 in w D Q v (31)
Velocity of gaseous fluids depends on the pressure in pipe since they are compressible (32):

      2 4 in n n g D p Q p v (32)
Velocities for water and gas for calculated flows through the pipes in our example is listed in table 10:

Table 10: Velocities for water and gas for calculated flows from example network

Conclusion

Comparison between analyzed methods was carried out, taking as a criterion the number of iteration to achieve final results. Modified Hardy Cross method, Modified M.M. Andrijahshev method and Node-Loop method have equal performances according to above adopted criterion.

For more complex networks, using Node-Loop method, number of required iterations is smaller even compared with Modified Hardy Cross method. Among these three methods, Node-Loop method is superior because it does not require complex numerical scheme for algebraic addition of corrections in each of iterations. In Node-Loop method, final results of each of the iterations are flows directly and not correction of flows. Modified Andrijašev method is more complex compared with Modified Hardy Cross method but without any improvement in the properties of the convergence. Node method has the worst performance of convergence, but this method is different in its approach compared with all other methods shown in this paper. Node method cannot be rejected easily based only on calculation shown in this paper. Hardy Cross method in original form has historical value and should be replaced with Modified Hardy Cross method, or even better with Node-Loop method.

Node-Loop method, presented among others in the text, is powerful numerical procedure for calculation of flows in looped fluid distribution networks. Main advantage is that flow in each pipe can be calculated directly, which is not possible after other available methods. In other methods, results of calculation are flow corrections which have to be added to flows calculated in previous iteration using complex algebraic rules. Node-Loop method is recommended to be used.

In the real network, consumers are not concentrated in the nodes. This can cause two-way flow in some pipes [START_REF] Brkić | An Improvement of Hardy Cross method applied on looped spatial natural gas distribution networks[END_REF]. This can cause disturbance in convergence properties of certain method. In such case, method should be changed. Some details on convergence properties can be found in the paper of [START_REF] Mah | Pipeline network calculations using sparse computation techniques[END_REF], [START_REF] Mah | Comparison of modified Newton"s methods[END_REF] and [START_REF] Altman | Convergence of Newton method nonlinear network analysis[END_REF].

Simulation problem today can be solved using different software [START_REF] Huddleston | Water distribution network analysis using Excel[END_REF][START_REF] Lopes | Implementation of the Hardy-Cross method for the solution of piping networks[END_REF]) in which can be implemented shown methods. Method for solution of pipe equations proposed by [START_REF] Hamam | Hybrid method for the solution of piping networks[END_REF] for gas networks and [START_REF] Todini | A gradient method for the analysis of pipe networks[END_REF] for water networks is available, but not shown here.

Acknowledgements

pp  F'=   CQ Q   c ΔQ 1 = ' F F d ΔQ 2 e Q 1 =Q b C= 22 12 pp  F'=   CQ Q   c ΔQ 1 =
pp  F'=   CQ Q   c ΔQ 1 d ΔQ 2 e Q 1 =Q b C= 22 12 pp  F'=   CQ Q   c ΔQ 1 d ΔQ 2 Q 2 =Q I -19444·10 -5 -1264933339 11839776055 +15114·10 -5 - -4331·10 -5 -82226369 3455539181 +2078·10 -5 - -2253·10 -5 +2778·10 -5 +20357137 1333799622 +15114·10 -5 +7937·10 -5  +25828·10 -5 +1178151884 8301891430 +2078·10 -5 -5572·10 -5 = +22263·10 -5 -30556·10 -5 -2399620963 14293015047 +15114·10 -5 -4288·10 -5  -19730·10 -5 -1082435403 9985057887 +2078·10 -5 -4154·10 -5 ‡ -17620·10 -5 Σ F I =-3644197166 27466590725 F I =+13490112 21742488499 II +27778·10 -5 +1344982710 8812326713 -7937·10 -5 - +19841·10 -5 +729037716 6687432899 +5643·10 -5 - +25484·10 -5 -27778·10 -5 -200615476 1314432601 -7937·10 -5 - -35715·10 -5 -316967670 1615251716 +5643·10 -5 - -30072·10 -5 -2778·10 -5 -20357137 1333799622 -7937·10 -5 -15114·10 -5  -25828·10 -5
-1178151884 8301891430

+5643·10 -5 +1128·10 -5  -22263·10 -5 +2778·10 -5 +1828426 119798452 -7937·10 -5 -4288·10 -5 = -9447·10 -5f -16966070 326858346.3 +5643·10 -5 -4154·10 -5 ‡ -3772·10 -5 Σ F II =+1125838522 11580357388 F II =-783047908 16931434392 III -2778·10 -5 -1828425 119798452 +4288·10 -5 +7937·10 -5 ‡ +9447·10 -5f +16966070 326858346.3 -32·10 -5 -5572·10 -5 = +3772·10 -5 +2778·10 -5 +65604940 4298435730 +4288·10 -5 - +7066·10 -5 +358812846 9242405949 -32·10 -5 - +7034·10 -5 +30556·10 -5
+2399620963 14293015047 +4288·10 -5 -15114·10 -5 = +19730·10 -5 +1082435403 9985057887 -32·10 -5 +1128·10 -5 = +17620·10 -5 -16667·10 -5 -2096864105 22897756035 +4288·10 -5 --12379·10 -5 -1220331619 17942054130 -32·10 -5 --12411·10 -5 Σ F III =+366533372 41609005264 F III =+237882701 37496376312 Pipe lengths and diameters are shown in figure 1 and initial flow pattern a in figure 2; 3 iterations are enough to achieve results shown in table 9 (this is for ΣF→0) b this is F calculated using (1), c using ( 12) where (Δp) is replaced with (C), d ΔQ 2 is ΔQ 1 from adjacent loop, e final calculated flow in the first iteration is used for the calculation in the second iteration, f opposite flow direction than in previous iteration

Q b Δp=p 1 -p 2 F'=   pQ Q   c ΔQ 1 d ΔQ 2 e Q 1 =Q b Δp=p 1 -p 2 F'=   pQ Q   c ΔQ 1 d ΔQ 2 Q 2 =Q I 3 -19444·10 -5 -5279095 54299272 +14195·10 -5 - -5250·10 -5 -416024 +2505·10 -5 - -2744·10 -5 4 +2778·10 -5 +69146 4978573 +14195·10 -5 +8006·10 -5  +24978·10 -5 +4800902 +2505·10 -5 -5312·10 -5 = +22172·10 -5 7 -30556·10 -5 -10718549 70157780 +14195·10 -5 -4053·10 -5  -20413·10 -5 -4840373 +2505·10 -5 -107·10 -5  -18015·10 -5 Σ F I =-15928498 129435625 F I =-455494 II 1 +27778·10 -5 +5919850 42622924 -8006·10 -5 - +19772·10 -5 +3030676 +5312·10 -5 - +25084·10 -5 2 -27778·10 -5 -816585 5879413 -8006·10 -5 - -35784·10 -5 -1342790 +5312·10 -5 - -30472·10 -5 4 -2778·10 -5 -69146 4978573 -8006·10 -5 -14195·10 -5  -24978·10 -5 -4800902 +5312·10 -5 -2505·10 -5  -22172·10 -5 5 +2778·10 -5 +5967 429677 -8006·10 -5 -4053·10 -5 = -9281·10 -5 -57430 +5312·10 -5 -107·10 -5  -4076·10 -5 Σ F II =+5040086 53910587 F II =-3170446 III 5 -2778·10 -5 -5967 429677 +4053·10 -5 +8006·10 -5 ‡ +9281·10 -5 +57430 +107·10 -5 -5312·10 -5 = +4076·10 -5 6 +2778·10 -5 +232186 16717415 +4053·10 -5 - +6830·10 -5 +1312311 +107·10 -5 - +6937·10 -5 7 +30556·10 -5 +10718549 70157780 +4053·10 -5 -14195·10 -5 = +20413·10 -5 +4840373 +107·10 -5 -2505·10 -5 = +18015·10 -5 8 -16667·10 -5 -8874257 106491089 +4053·10 -5 - -12614·10 -5 -5136681 +107·10 -5 - -12507·10 -5 Σ F III =+2070510
193795963 F III =+1073434 Pipe lengths and diameters are shown in figure 1 and initial flow pattern a in figure 2; 3 iterations are enough to achieve results shown in table 9 (this is for ΣF→0) b this is F calculated using (3), c using (12), d ΔQ 2 is ΔQ 1 from adjacent loop, e final calculated flow in the first iteration is used for the calculation in the second iteration, f opposite flow direction than in previous iteration 17), b F'=|∂Q(C)/∂C|, c ΔC 1 after eq. ( 20), d ΔC 2 is ΔC 1 from adjacent node, e final calculated pressure function in the first iteration is used for the calculation in the second iteration, f opposite flow direction than in previous iteration 2409·10 -10 -10320·10 -5 654·10 -10 Pipe lengths and diameters are shown in figure 1; see figure 4 for initial pattern (blue letters); 9 iterations are enough to achieve results shown in table 9 (this is for Σf→0) a using (18); to calculate friction factor λ, Reynolds number has to be calculated and for this velocity have to be chosen (this velocity does not have effect on final results, here is chosen extremely large velocity 100 m/s), b F'=|∂Q(Δp)/∂Δp|, c Δ Δp1 after eq. ( 20), d Δ Δp2 is Δ Δp1 from adjacent node, e final calculated pressure function in the first iteration is used for the calculation in the second iteration, f opposite flow direction than in previous iteration

pp  F'=   CQ Q   1 co Q  2 do Q  e Q 1 =Q b C= 22 12 pp  F'=   CQ Q   1 co Q  2 do

  [START_REF] Ekinci | An optimization strategy for water distribution networks[END_REF]) can be found very easily using Darcy-Weisbach equation (3):

Figure 1 .

 1 Figure 1. Example of pipeline network

  figure 2.

Figure 2 .

 2 Figure 2. Example network initial parameters prepared for loop oriented calculation

Figure 3 .

 3 Figure 3. Contours for method of M.M. Andrijašev calculation

  in a front of flow Q in the matrix [Q] means that calculated flow direction is opposite compared to shown one in the previous iteration (or in the figure 1 in our case for the first calculated values of flows compared with initials flow pattern). If all values of pressure drops sums calculated after (24) are not approximate zero with reasonable accuracy, calculation has to be repeated using values calculated in previous iteration. At the end of calculation calculated set of flows [Q] stays unchanged.

( 30 )

 30 as example. These values are valid for the first iteration. First five rows (first matrix) are from node equation, and next three is from loop equation but multiplied with first derivate marked in tables as F":

  027·119798452+0.027·4298435730+0.305·14293015047+(-0.166·22897756035)) After first iteration for gas network vector of flows is [0.198409265, 0.357146291, 0.043307855, 0.25828288, -0.094469817, 0.07065686, 0.197298048, 0.123787585] T . Minus in front of flow in pipe 5 means: change assumed flow direction from previous iteration. After first iteration for water network vector of flows is [0.197719798, 0.357835758, 0.052496097, 0.249784106, -0.092806697, 0.068304272, 0.204133702, 0.126140172] T . Flows are expressed in m 3 /s.

Figure 5 .

 5 Figure 5. Comparisons of convergence for presented methods (gas network)

Figure 2 .

 2 Figure 1. Example of pipeline network

Table 1

 1

. Hardy Cross calculation for example gas network Table

2

. Hardy Cross calculation for example water network In the original Hardy Cross method, each loop correction is determined independently of other loops. Original Hardy Cross method is a sort of Newton-Raphson method but used to solve each single loop equation solely, one by one. Hence, Hardy Cross method is also known as the single contour adjustment method. In matrix form, original Hardy Cross method for the example network for water distribution from figure 1 can be noted as (10):

Table 5 .

 5 Calculation after Modified M.M. Andrijašev method for example gas network

	4 4048956589	3 2586239214	3 1171997760		Q	o	-	3277663792
								I
	3 2586239214	8 5294976574	-	4166255733	x	Q	o	1492371894	(16)
								II
	3 1171997760	-	4166255733	8 3637934886		Q	o	-	2518358643
							III
	For example term in first row and second column from (16) is
	25862392143=4298435730+22897756035-1333799622 according to table 5. Same value has
	term in second row, first column, etc. Presented matrix is symmetrical. Similar can be done
	for water network (table 6).

Table 6 .

 6 Calculation after Modified M.M. Andrijašev method for example water network M.M. Andrijašev method does not improve convergence properties compared with Modified Hardy Cross method. It only complicates calculation of non-diagonal terms in the matrix

Table 1 .

 1 The Hardy Cross calculation for example gas network

	Iteration 1

Table 3 .

 3 Calculation after the modified Hardy Cross method for example gas network

	Iteration 1

table 3 DB Click here to download table: Table 3 DB.docTable 4 .

 34 Calculation after the modified Hardy Cross method for example water network

	Iteration 1

table 4 DB Click here to download table: Table 4 DB.docTable 5 .

 45 Calculation after Modified M.M. Andrijašev method for example gas network

	Iteration 1

Table 7 .

 7 Calculation after Modified Node method for example gas network with three loops Pipe lengths and diameters are shown in figure1; see figure4for initial pattern (red letters); 26 iterations are enough to achieve results shown in table 9 (this is for Σf→0)

	Iteration 1

a using (

table 7 DB Click here to download table: Table 7 DB.docTable 8 .

 78 Calculation after Modified Node method for example water network with three loops

table 8 DB Click here to download table: Table 8 DB.docTable 9 .

 89 Final flows for network presented in this paper .73 94.86 802.87 -146.23 a 248.50 643.36 451.50 Gas 913.72 1086.28 82.01 804.27 -137.86 a 251.58 633.60 448.42 a sing minus means flow direction opposite than first assumed in figure 2

	Flows (m 3 /h)
	Pipe	1	2	3	4	5	6	7	8
	Water	902.27 1097

table 9 DB Click here to download table: Table 9 DB.doc

 9

Table 10 :

 10 Velocities for water and gas for calculated flows from example network

	Velocity (m/s)
	Pipe	1	2	3	4	5	6	7	8
	Water	21.0	15.4	2.2	18.6	-1.7 a	9.4	14.9	13.2
	Gas	5.3	3.8	0.5	4.7	-0.4 a	2.4	3.7	3.3
	a sing minus means flow direction opposite than first assumed in figure 2

table 10 DB Click here to download table: Table 10 DB.doc

-10 -11 -11 -11 -11 -10 -11

List of tables:

Table 1. Hardy Cross calculation for example gas network Table 2. Hardy Cross calculation for example water network Table 3. Calculation after the modified Hardy Cross method for example gas network Table 4. Calculation after the modified Hardy Cross method for example water network Table 5. Calculation after Modified M.M. Andrijašev method for example gas network Table 6. Calculation after Modified M.M. Andrijašev method for example water network Table 7. Calculation after Modified Node method for example gas network with three loops Table 8. Calculation after Modified Node method for example water network with three loops

-1264933339 11839776055 +13268·10 -5 --6177·10 -5 -156904917 -1128·10 -5 --7305·10 -5 4 +2778·10 -5 +20357137 1333799622 +13268·10 -5 +9722·10 -5  +25767·10 -5 +1173109335 -1128·10 -5 -5572·10 -5 = +19067·10 -5

+2399620963 14293015047 -881·10 -5 -13268·10 -5 = +16407·10 -5 +773797561 +4154·10 -5 +1128·10 -5  +21689·10 -5 8 -16667·10 -5 -2096864105 22897756035 -881·10 -5 --17548·10 -5 -2302927589 +4154·10 -5 --13394·10 -5 Σ F III =+366533372 41609005264 F III =-1488793175 Pipe lengths and diameters are shown in figure 1 and initial flow pattern a in figure 2; 6 iterations are enough to achieve results shown in table 9 (this is for ΣF→0) b this is F calculated using (1), c also using (9), d ΔQ 2 is ΔQ 1 from adjacent loop, e final calculated flow in the first iteration is used for the calculation in the second iteration, f opposite flow direction than in previous iteration

193795963 F III =-6463424 164880214 Pipe lengths and diameters are shown in figure 1 and initial flow pattern a in figure 2; 6 iterations are enough to achieve results shown in table 9 (this is for ΣF→0) b this is F calculated using (3), c also using ( 9) or (10), d ΔQ 2 is ΔQ 1 from adjacent loop, e final calculated flow in the first iteration is used for the calculation in the second iteration, f opposite flow direction than in previous iteration

-316967670 +3876·10 -5 +1767·10 -5 = -30072·10 -5 Σ F°=-2518358644 36379348868 F°=-769557797 Pipe lengths and diameters are shown in figure 1, initial flow pattern in figure 3 and numerical values for initial flows a in figure 2; 3 iterations are enough to achieve results shown in table 9 (this is for ΣF→0) b this is F calculated using (1), c using (15),

adjacent contour, e final calculated flow in the first iteration is used for the calculation in the second iteration, f opposite flow direction than in previous iteration -1342790 +3855·10 -5 +1457·10 -5 = -30472·10 -5 Σ F°=-10888413 173389067 F°=-3625941 Pipe lengths and diameters are shown in figure 1, initial flow pattern in figure 3 and numerical values for initial flows a in figure 2; 3 iterations are enough to achieve results shown in table 9 (this is for ΣF→0) b this is F calculated using (3), c using (15),


Dejan Mr

M.Sc Brkić

Dejan Brkić

An Improvement of Hardy Cross Method Applied on Looped Spatial Natural Gas Distribution Networks

Keywords: Hardy Cross method, Natural gas, Flow, Pipeline network Hardy Cross method, Natural gas, Flow, Pipeline network

An improvement of Hardy Cross method applied on looped spatial natural gas distribution networks

Nomenclature

Q a -gas flow reduced to absolute pressure (m 3 /s), Q-gas flow reduced to standard value of pressure and temperature (m 3 /s) v-gas flow velocity (m/s), A c -surface value of cross section (m 2 ), -mark complementary value of Bull's variable (if Bull's variable u has value 1, then its complementary value ū is 0, and opposite)

Text
C-auxiliary Bull's variable (1 or 0)

Introduction

For single source, branching networks, only the reach equations need to be solved (e.g. using

Hazen-Williams for water or Renouard for natural gas). In this case, optimization by linear programming could yield the best alternative network, both hydraulically and economically [START_REF] Mcclure | Linear-programming offers way to optimize pipeline analysis[END_REF].

Critical path analysis for optimizing branched pipe networks and implementation of the finite element method can be done for branched networks. For looped networks, however, techniques that are more powerful are required. Algorithms have been developed to ensure convergence of the iterative procedures.

Today, Hardy Cross method [START_REF] Cross | Analysis of flow in networks of conduits or conductors[END_REF] is very often used for optimization of gas distribution networks with loops of conduits. Hardy Cross 1 developed a numerical method for calculating flow and pressure distribution in a looped network. This method also has been widely used in modeling of waterworks with ring-like structures of conduits in municipalities. Hardy Cross method assumes an equilibrium between pressure and friction forces in steady and incompressible flow. As a result, it cannot be successfully used in unsteady and compressible flow calculations with large pressure drop where inertia force is important. Original and improved Hardy Cross method, are methods of successive iterative corrections, but for the first one, corrections are calculated for each contour separately, while for the second one, corrections for all contours in each iteration is calculated simultaneously using the Newton-Raphson numerical procedure. Hardy Cross method is one approach to solve the loop equations. Original method first determines corrections for each loop independently and applies this correction to compute new flow in each conduit. It is not efficient compared to improved Hardy Cross method that considers entire system simultaneously. Simultaneous method is used by Epp and Fowler [START_REF] Epp | Efficient code for steady flows in networks[END_REF], but only for looped waterworks systems. 1 Hardy Cross (Nansemond County, Virginia 1885 -Virginia Beach, Virginia 1959), American engineer, assistant professor of civil engineering at Brown University, professor of structural engineering at the University of Illinois, Chair of the Department of Civil Engineering at Yale. He had also developed a method for analysing indeterminate structures that minimized the inconveniences and risks involved in the use and development of reinforced concrete.

Some methods developed by Russian authors are similar with original Hardy Cross method.

Contemporary with Hardy Cross, soviet author V.G. Lobachev [START_REF] Latysenkov | Hydraulics[END_REF] was being developed very similar method compared to original Hardy Cross method. Andrijashev method [START_REF] Andriyashev | Hydraulics calculation of water distribution networks[END_REF] was very often being used in Russia during the soviet era. According to this method, contour and loop are not synonyms (contours for calculations has to be chosen to include few loops and only by exception one).

Using the loop equations to represent conservation of energy, Wood and Charles [START_REF] Wood | Hydraulic network analysis using linear theory[END_REF] developed a linear theory (flow adjustment) method by coupling the loop equations with node equations.

Convergence characteristics of linear theory are later improved by Wood and Rayes [START_REF] Wood | Reliability of algorithms for pipe network analysis[END_REF].

Modified linear theory solves directly for the conduits flow rates rather than the loop equations approach of Hardy Cross method. Shamir and Howard [START_REF] Shamir | Water distribution systems analysis[END_REF] solved node equations instead loop equations using the Newton-Raphson method. After the nodal heads are computed, they computed the conduits flow rates. Previous methods solve for the conduits flows or nodal heads separately then use conservation of energy to determine the other set of unknowns. Haman and Brameller [START_REF] Hamam | Hybrid Method for the Solution of Piping Networks[END_REF], and, Todini and Pilati [START_REF] Todini | A gradient method for the analysis of pipe networks[END_REF] devised a method to solve for flows and heads simultaneously. Here, each conduit equation is written to include, both, the conduit flows and nodal heads. In addition, although the number of equations is larger than the other methods, the algorithm does not require defining loops. Further procedure is developed by Patankar [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. He developed a finite volume procedure to solve for Navier-Stokes equations in a structured coordinate system. Since the publication of the original paper in 1972, there have been several

developments reported to improve the numerical performance of the original algorithm. Datta

and Majumdar [START_REF] Datta | Flow distribution in parallel and reverse flow manifolds[END_REF] used this solution algorithm to develop a calculation procedure for manifold flow systems.

Examples of calculation of looped natural gas distribution network after original Hardy Cross method can be found in handbooks [START_REF] Corfield | Distribution Design for Increased Demand, Chapter 9[END_REF] and similar calculation of real gas network is shown in paper of Manojlović et al. [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF]. But, deeper improvement of Hardy Cross method is only shown in case of looped waterworks systems [START_REF] Boulos | Comprehensive water distribution systems analysis handbook for engineers and planners[END_REF]. Analysis of looped gas distribution networks is shown in paper of Osiadacz [START_REF] Osiadacz | Simulation and analysis of gas networks[END_REF] while in many papers are shown methods for calculation of waterworks [START_REF] Arsene | Modelling and simulation of water systems based on loop equations[END_REF][START_REF] Collins | Solving the pipe network analysis problem using optimization techniques[END_REF][START_REF] Johnson | Finite-element method for water distribution networks[END_REF][START_REF] Gay | The solution of pipe network problems[END_REF][START_REF] Chiplunkar | Analysis of looped water distribution networks[END_REF][START_REF] Gupta | Optimization of water distribution system[END_REF][START_REF] Kessler | Analysis of the linear programming gradient method for optimal design of water supply networks[END_REF][START_REF] Varma | Optimal design of water distribution systems using an NLP method[END_REF][START_REF] Chenoweth | Pipe network analysis[END_REF][START_REF] Eiger | Optimal design of water distribution networks[END_REF][START_REF] Basha | Analysis of water distribution systems using a perturbation method[END_REF][START_REF] Todini | Looped water distribution networks design using a resilience index based heuristic approach[END_REF][START_REF] Ahuja | Some recent advances in network flows[END_REF]. Very interesting application of two ant colony optimisation algorithms to water distribution system optimisation is shown in the paper of Zecchin et al. [START_REF] Zecchin | Application of two ant colony optimisation algorithms to water distribution system optimisation[END_REF].

Concept of distribution network

Hardy Cross method is powerful toll for calculation of looped gas distribution network in settlements without limitation factors, such as: number of conduits per contours, number of loops, number of nodes or number of input nodes. In Hardy Cross calculation, previously, has to be determinated maximal consumption per each node (Q output ), and one or more inlet nodes (Table 1). These parameters are looked up. Now, initial guess of flow per conduits has to be assigned (Table 2), and in that way chosen values are to be used for first iteration. After the iteration procedure is completed, and if the value of gas flow velocity for all conduits are bellow standard values, calculated flows 2 become flow distribution per conduits for maximal possible consumptions per nodes. Further, pressure per all nodes can be calculated. Whole network can be supplied by gas from one or more points (nodes). Distribution network must be design for largest consumption assigned to nodes of networks chosen to satisfy larges possible gas consumption of households in condition of very severe winters (Table 1) 3 . Disposal of households is along the network's conduits, and only their consumption is to be assigned to nodes for the purpose of calculation according to Hardy Cross. Task of Hardy Cross method is to calculate gas flow distribution per conduits for fixed looked up maximal natural gas consumption per nodes and inputs in network. If the flow speed in some conduits are above standard values after the calculation is finished, these conduits must be chosen with larger diameters, and whole iteration procedure have to be repeated 4 . Hardy Cross method is only suitable for calculation of looped networks, not for branch-like networks. Before start of iterative procedure, consumption of gas per node, input nodes, their spatial disposal and length of conduits must be chosen and looked up. Then, first guess of flow for all conduits must be chosen according to first Kirchhoff's law for each of the node nodes. Next, diameters of all conduits have to be chosen (after eq. ( 1); see I in Table 2, or not; see II in Table 2).

π v Q 4 D 4 π D v A v Q a in 2 in c a          (1)
In previous equation [START_REF] Mcclure | Linear-programming offers way to optimize pipeline analysis[END_REF], flow is reduced to absolute value of pressure, and if flow is given for standard value is to be reduced using p a •Q a =p n •Q n . In distribution network pressure value is usually p a =4•105 Pa abs., or 3•10 5 Pa gauge pressure. First larger standard diameter than calculated after eq. ( 1) has to be chosen (du~Du) from the tables of standard polyethylene pipes.

Each equation for determination of pressure drop in distribution gas network [START_REF] Coelho | Considerations about equations for steady state flow in natural gas pipelines[END_REF] (e.g. in the technical literature the Renouard equation for conditions of pressure values in gas distribution networks comes written in an explicit form in terms of the pressure difference) (2) 5 :

82 . 4 in 82 . 1 r 2 1 2 2 D Q L 4810 p p       (2)
can be written in form Δp 2 =R•Q n , where is: n=1.82, and hydraulics resistance R is to be written

as R=4810•ρ r •L•D in -4.82 .
Basics assumption is to be satisfied for Hardy Cross calculation:

-Algebraic sum of flows per each node must be zero exactly (first Kirchhoff's law-continuity of flow), -Algebraic sum of pressure drops per each contour must be approximately zero at the end of iterative procedure (second Kirchhoff's law -continuity of potential).

It is usually desired to determine the total loss of pressure or voltage between inlet and outlet. If a single conductor connected these two points, the loss of head for given flow could be computed directly from the relation between flow and head loss. In a network, however, this loss depends on the distribution of the flow in the system. If such distribution is known, the drop of potential in each conductor can be determined directly, and the total drop found as the sum of the drops along any path connecting inlet and outlet, the total drop being of course the same whatever path is chosen. If, however, the relation is not linear as for gas flow (2), serious difficulties arise in solving the equations. This problem must be solved during iterative procedure. 5 whereas in the technical literature it appears as: Δp 2 =4088•ρ r 0.82

•L•Q 1.82 •D in -4.82
In general, systems for distributing natural gas in cities may, for purposes of analysis, be considered as in a single plane. In other cases, as, for example, the distribution may take place in several planes, with interconnection between the planar systems of distribution (Fig. 1). This type of problem presents no especially new features except that distribution must be made in circuits closed by the risers as well as in the circuits which lie in a plane. It will be noted that in such problems any conduit may lie in only one circuit (an outside conduit in a floor) or in two circuits, three circuits, or even in four circuits (two floor circuits and two riser circuits). The total change in flow in the conduit is the sum of the changes in all the circuits of which it is a member.

All networks with three or more dimensions can be noted as spatial. Spatial network is each network with at least one conduit mutual for three or more contours. In fig. 1,networks A) can be reduced to single plane problem, but case C) is three-dimensional problem as in our example in fig 2 6 .

Mathematical description of network

The first step in solving a problem is to make a network map showing conduit sizes and lengths, connections between conduits (nodes), and sources of supply. For convenience in locating conduits, assign each contour and each main a code number. Conduits on the network periphery are common to one contour and those in the network interior are common to two contours.

Special cases may occur in which two conduits cross each other but are not connected, resulting in certain conduits being common to three or more loops. The distribution network then becomes three-dimensional rather than two-dimensional. for example, by starting at the right end of conduit 15 via conduit 8 and returning to the left end of conduit 15 via conduit 2 or 3.

Gas distribution system, composed of fifteen conduits (Fig. 2), has been analyzed by the Hardy Cross method (original and improved version) to determine the individual conduit flow rates and pressure drops. Gas flow into the network from a source on the left side is 7000 m 3 /h (Table 1.).

Points of delivery are at junctions of conduits, with the arrows pointing to volumes delivered (summation of these deliveries equals 7000 m 3 /h). Assumed gas flow and its direction, also indicated by an arrow (near pipe in fig. 2).

After the network map with its conduit and contour numbers and delivery and supply data has been prepared, the next step is to assume a flow pattern in the network (Table 2). This may be done by starting at sources with volumes of gas delivered into the system, and distributing these volumes through the conduits until they have been allocated to the various delivery points. The flows thus assumed are entered next to their respective conduits, with arrows to indicate direction. The total gas flow arriving at a junction must equal the total gas flow leaving it (first Kirchhoff's law). The assumed flow pattern will approximate the correct flow pattern if consideration is given to the relative flow capacity of various network conduits.

To introduce matrix form in calculation, it is necessary to represent distribution network (Fig. 2) as a graph according to Euler's theorem from mineralogy (number of polyhedral angles and edges of minerals). Graph has X branches and Y nodes (in Fig. 2

                                                                                                                                       output XI output X output IX output VIII output VII output VI output V output IV output III output II Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q 15 14
0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 (3)
In (3); 1 means that particular conduit has input in observed node, -1 if not and 0 if this particular conduit do not belong to observed node. If network has two or more points of supply, that can be included in eqs. (3) as Q i but with negative sign. Equation ( 3) can be written as (4) in develop form:

nodes Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q output XI
                                                     (4)
For independent loops can be written set of X-Y+1 independent equations of energy continuity for network shown in fig. 2. can be written as (5):

0 * 1 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 2 15 2 3 2 2 2 1                                          p p p p  (5)
In ( 5); 1 means that flow direction in conduit coincides with assumed direction of observed contour, -1 if not, and 0 if conduit do not belong to observed contour. Equation (3) can be written as (4) in develop form:

contours F F F F F 0 D Q L D Q L D Q L D Q L D Q L 4810 p p p p p 0 D Q L D Q L D Q L D Q L 4810 p p p p 0 D Q L D Q L D Q L D Q L D Q L 4810 p p p p p 0 D Q L D Q L D Q L D Q L D Q L 4810 p p p p p 0 D Q L D Q L D Q L D Q L
                                                                                                                                                                            (6)
Equations [START_REF] Epp | Efficient code for steady flows in networks[END_REF][START_REF] Latysenkov | Hydraulics[END_REF][START_REF] Andriyashev | Hydraulics calculation of water distribution networks[END_REF][START_REF] Wood | Hydraulic network analysis using linear theory[END_REF] are for initial guess of gas flow rate I from Table 2 andFig

General solution after Hardy Cross method

Two types of methods based on Hardy Cross's idea are shown for the solution of loop equations:

-Hardy Cross method; successive substitution method (single loop adjustment method) -Modified Hardy Cross method; simultaneous loop solution method (Newton-Raphson method)

Hardy Cross's successive substitution method (single loop adjusment method)

Equations of flow (for conduits i which belong to contour j) for all contours (each particular contour is marked by j) can be written as follows [START_REF] Wood | Reliability of algorithms for pipe network analysis[END_REF]:

      j n i i j Q R Q F    ( 7
)
The basic idea of the Hardy-Cross method is that conservation of mass at each node can be established initially. This means we must first assume an initial guess of flows in every pipe element before starting the pressure drop calculation. For any pipe in which Q 0 is assumed to be the initial flow rate, eq. ( 7) can be estimated using a Taylor series expansion as follows (8):

                j j i j m i m i m j i m j i Q Q F Q Q Q F Q F            ! 1 ) 1 ( ) ( 1 (8)
Because of second Kirchof's law in late iterations; F(Q) (m) →0, and

  j m i m i j j Q Q Q 1               0 1 1           ) ( m Q j j i j m Q Q F Q F ( 9
)
For contour (loop) j, where conduits i belong to loop j (10):

                  1 1 1 1 1                   m j n i i m j n i i Q j j i m j i j Q R n Q R Q Q F Q F m (10)
For Renouard equation ( 2): n=1,82. Equation [START_REF] Todini | A gradient method for the analysis of pipe networks[END_REF] in matrix form can be writen as [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]:

                                                                                                                              1 m V 1 m IV 1 m III 1 m II 1 m I m V m IV m III m II m I V 1 - m IV 1 - m III 1 - m II 1 - m I 1 - m F F F F F ΔQ F 0 0 0 0 0 ΔQ F 0 0 0 0 0 ΔQ F 0 0 0 0 0 ΔQ F 0 0 0 0 0 ΔQ F V IV III II I (11)
Detail example for calculation according Hardy Cross's successive substitution method (single loop adjusment method) for network shown in fig. 2 are given in table 4.

Modified Hardy Cross method (Newton-Raphson simulaneous method)

The Newton-Raphson method is a numerical method that can solve a set of equations [START_REF] Wood | Reliability of algorithms for pipe network analysis[END_REF] simultaneously. Convergence of Newton method nonlinear network analysis is studied by Altman and Boulos [START_REF] Altman | Convergence of Newton method nonlinear network analysis[END_REF]. The method is particularly convenient for solving differentiable equations when the value of the desired unknown parameters is known approximately. Using the Taylor's series expansion, a first-order approximation can be written as [START_REF] Datta | Flow distribution in parallel and reverse flow manifolds[END_REF]:

                                                                                                                                                                                                                                                    1 m V 1 m IV 1 m III 1 m II 1 m I m V m IV m III m II m I V 1 - m IV 1 - m III 1 - m II 1 - m I 1 - m V 1 - m IV 1 - m III 1 - m II 1 - m I 1 - m V 1 - m IV 1 - m III 1 - m II 1 - m I 1 - m V 1 - m IV 1 - m III 1 - m II 1 - m I 1 - m V 1 - m IV 1 - m III 1 - m II 1 - m I 1 - m F F F F F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F ΔQ F V V V V V IV IV IV IV IV III III III III III II II II II II I I I I I (12)
First matrix in [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] and ( 12) is a Jacobian matrix of first derivates of loop equations evaluated at Q (m-1) of all the functions for (12) and for only observed loop in [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. Jacobian matrix is square (number of contours) and symmetric. The rows in Jacobian matrix are corresponding to the loop equations and the columns correlated the loop corrections. The diagonal terms of Jacobian matrix are the sum of first derivates of the conduit equation in particular contour. The difference between successive substitution method and simultaneous solution method is that some of nondiagonal terms are non-zero in simultaneous solution method.

For the network shown in Fig. 2 Jacobian matrix is for the first iteration (13):

                                                                                                                                                      V
m ΔQ Q , Q , Q , Q , Q F ΔQ Q - F ΔQ Q Q F ΔQ Q Q F 0 ΔQ Q - F ΔQ Q , Q , Q , Q F 0 ΔQ Q - F ΔQ Q - F ΔQ Q Q F 0 ΔQ Q , Q , Q , Q , Q F ΔQ Q - F 0 ΔQ Q Q F ΔQ Q - F ΔQ Q - F ΔQ Q , Q , Q , Q , Q F ΔQ Q - F 0 ΔQ Q - F 0 ΔQ Q - F ΔQ Q , Q , Q , Q F J (13)
Numbers in [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF] are from first iteration shown in Table 4 (but note that Table 4 is for single loop adjustment method and only input parameters for first iteration have equal numerical values): [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF] Note that in eq. ( 14); 559997875 = 233875825 + 326122050, and 179192287 = 14916225 + 164276062 (see Table 4).

                      3045914618

Rules for determination of algebraic signs preceding corection of flow

These rules are to be applied for both version of method (shown in chapter 4).

Contours and conduits numbers are listed in the first and second columns of the table 4, respectively. Diameters and lengths of conduits are listed in the third and forth column of table 4.

The assumed gas flow in each conduit for iteration 1 (shown in Table 2) is listed in the fifth column in table 4. The plus or minus preceding the flow, Q, indicates the direction of the conduit flow for the particular contour. A plus sign denotes clocwisese flow in the conduit within the contour; a minus sign, counterclockwise.

The first computation-resistance to change in gas flow in each conduit, R•Q 0,82 -is listed in table 4 (R is according to (2). The coefficient n, which in this case equals 1.82 -from Renouard equation ( 2), has become n-1=0,82 (first derivate). The pressure drop in each main, R•Q 1.82 , is listed, and carries the same sign as the gas flow. Column R•Q 0,82 is added arithmetically for each contour. Column R•Q 1.82 is added algebraically for each contour. A flow correction, Δ, is computed for each contour - [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] for original method and (12) for improved. This correction must be subtracted algebraically from the assumed gas flow. A conduit common to two loops receives two corrections, and a main common to three or more contours receives three or more corrections.

Correction Δ 1 is from the particular contour under consideration. Corrections Δ 2 and Δ 3 are from the second and third contours to which a conduit belongs. The upper plus or minus sign shown indicates direction of flow in that conduit in these two contours and is obtained from Q for previous iteration. The upper sign is the same as the sign in front of Q if the flow direction in each contour coincides with the assumed flow direction in the particular contour under consideration, and opposite if it does not.

The flow, Q, and corrections are totaled across to obtain the Q listed under next iteration according to the following rules:

1. The algebraic operation for correction 1 should be the opposite of its sign; i.e., add when the sign is minus.

2. The algebraic operation for corrections 2 and 3 should be the opposite of their lower signs when their upper signs are the same as the sign in front of Q, and as indicated by their lower signs when their upper signs are opposite to the sign in front of Q.

Rules for determination of preceding sign for flow corrections is shown in Table 3. Rules shown in table 3 can be presented by following logical equation; Bull's logic [START_REF] Boulos | Comprehensive water distribution systems analysis handbook for engineers and planners[END_REF] and ( 16):

D G K D G K D G K D G K A             (15)
D C D C A     (16)
Where in ( 16) is auxiliary relation for [START_REF] Arsene | Modelling and simulation of water systems based on loop equations[END_REF]:

G K G K C     (17)
Some possible logical circuits according previous equations ( 15) and ( 16) are shown in Fig. 3.

Computation according to these rules is by an algebraic subtraction of the flow correction terms.

The iteration procedure is repeated until the net pressure drop around each loop, is as close to zero as the degree of precision desired demands. The network is then in approximate balance.

Pressure drops in conduits along a path from the node of lowest pressure to the supply source are summed to obtain the total pressure drop in the network. When the gas flow in a conduit is in the same direction as the path taken between two nodes, there is a pressure loss. When the flow is opposite to the direction of the path, there is a pressure gain. Since the network is in approximate balance, the total pressure drop should be computed along several paths and averaged to obtain a better value.

Results and identification of possible problems

Detail calculation (first two iterations -1 st guess, network in Fig. 2) after original Hardy Cross method (shown in subchapter 4.1) is given in Table 4. Corrections are calculated using eq. ( 15).

Calculation after improved Hardy Cross method is not given here in such detailed table, but calculation is done [START_REF] Datta | Flow distribution in parallel and reverse flow manifolds[END_REF] and compared graphically in figs. 4-6 with (11). Set of corrections calculated using eq. ( 12) for first iteration (1 st guess, network in Fig. 2) after modified Cross method (shown in subchapter 4.2) are Δ 1 =-44/325, Δ 2 =-7/93, Δ 3 =-48/577, Δ 4 =135/862, and Δ 5 =103/538. Modified Hardy Cross method has better convergence performance for all contours (approximately 3-5 increased speed of convergence in our case) for both guess -1 st and 2 nd .

Same conclusions can be done for convergence of flow corrections (shown in fig. 5 for contour I and in fig. 6 for contour IV).

Example of symmetric network is good example to solve some misunderstandings (initial and final flow pattern are shown in fig. 7, node consumption in Table 1; 3 rd guess, initial flow in Table 2, and final results in Table 5). This network is selected very carefully because under above-mentioned circumstances, in conduit 6 two-way flow must be expected (Fig. 8) [START_REF] Kaluđerčić | Problem of two-way supplied pipes in a gas network[END_REF]. But in Hardy Cross calculation that kind of flow in conduit is forbidden and cannot be calculated.

Anyway, after the calculation of this network is finished, value of flow in conduit 6 for 3 rd guess is 0 m 3 /s. But real consumers are located between node XI and node VII. That implies that in real network, in conduit 6 some value of flow must be expected (two-way supplied pipes) or some household will be left gasless. In some rare cases, convergence of Hardy Cross method can be spoiled [START_REF] Altman | Convergence of Newton method nonlinear network analysis[END_REF]. Recommendations for these cases is to changed method for calculation, i.e. original Hardy Cross instead of improved Hardy Cross, or opposite, or to applied some of other available methods. Even, if the node consumption is satisfied, in some regimes of exploitation of network, some households can be felt lack of gas. These households are place in the middle of conduits, between the nodes. In two-way supplied conduits this case can be occurred (Fig. 9). Also, in our case of symmetric network, convergence is stabile as for 1 st and 2 nd guess. Under some special cases, modified Hardy Cross method has not always better convergences characteristics in comparisons to original Hardy Cross method.

Conclusions

Approach of Hardy Cross was extremely practical. His view was that engineers lived in a real world with real problems and that it was their job to come up with answers to questions in design even if approximations were involved. Hardy Cross method procedure can give good results when designing a looped gas-pipeline network of composite structure. According to the price and velocity limits, the optimal design can be predicted. This paper addresses to the problem of construction of networks for distribution of natural gas in the cities and with subject to all the practical requirements for the engineers charged with design and/or analysis of such system. This paper is especially addressed to those engineers willing to understand and interpret the results of calculation properly and to make good engineering decision based on this subject. While the spatial natural gas networks with loops are maybe purely hypothetic, this kind of networks, but with some adjustment can find application in calculation of mines ventilation systems. a constant in calculation (network -Fig. 2) b constant in calculation (network -Fig. 7) c Flow into the network is sum of consumptions per nodes a must be chosen to satisfy rule after first Kirchhoff's law for all nodes (exactly) b initial guess means: distribution of flows for first iteration c chosen after eq. ( 1) -(1 st guess, network Fig. 2) d chosen to satisfy only mandatory first Kirchhoff's law for all nodes, but not chosen after eq. ( 1) -(2 nd guess, network Fig. 2) e flow pattern chosen for symmetric node consumption (see Table 1) -(3 rd guess, network Fig. 7) 1.82 and R•Q con 0.82 are not shown explicitly as in iteration 1 b Δ 1 must be added with opposite preceding sign, Δ 2 and Δ 3 must be added with adopted preceding sign (according to rule shown in table 3) c first assumed (initial) guess (flows per all conduits must be chosen to satisfy first Kirchhoff's law for all nodes -sum of pressure drops *** per all loops according to second Kirchhoff's law must be (approximate) zero when the network is in balance (here after 15 iterations) -see a calculation is over when second Kirchhoff's law is approximately satisfied (after n iterations) for all contours (loops) -algebraic sum of pressure drop per conduit is approximately equal zero. Note that final flow is not depend on first assumed gas flow and chosen type of methods, b must be under standard values [START_REF] Manojlović | Optimized design of a gas-distribution pipeline network[END_REF] -The velocity limits are 6 m/s for the pipes of small diameter (up to 0,09 m) and 12 m/s for the pipes of large diameter (up to 0,225 m) (if not must be changed diameter of conduit and must be repeated whole calculation) c symmetric network, Fig. 7 d velocity limit is exceeded in conduit 3 (must be increased diameter of conduit and must be repeated whole calculation)

a c Q (1) R•Q 1,82 (2) R•Q ,0,82 Flow correction Δ b Q Flow correction Δ b Q 1 2 3 1 2

DB ver 3 D

 3 in -inner diameter of conduit (m), π-Ludolph's number (=3.14159), p 2 -pressure at conduit exit (Pa), p 1 -pressure at conduit entrance (Pa), ρ r -relative gas density (-), L-conduit length (m), R-hydraulic resistance, flow resistance (m/m n ) (in electric networks equivalent Ω) n-flow equation exponent (n=1,82 in Renouard equation), m-mark observed iteration, X-number of conduits (pipes, main, manifolds) Y-number of nodes i, j, k-counter Δ-flow correction (m 3 /s) A, K, G, D-Bull's variables (1 or 0)

Figure 2 7

 2 Figure 27 is an example of a three-dimensional network because conduit 15 is not connected to

  2.

Figure 1 .Figure 2 .Figure 3 .Figure 4 . 2 )Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .

 1234256789 Figure 1. Transformations of pseudo-spatial network A) and spatial network C)

Figure 1 DB ver 3

 13 Figure 1 DB ver 3 Click here to download high resolution image

  Fig 4 d change of flow direction (opposite than in previous iteration -opposite upper sing in Δ 2 and Δ 3 )

  : X=15, Y=11). Graph with n nodes has Y-1 independent nodes and X-Y+1 independent loops. Tree is a set of connected branches chosen to connect all nodes, but not to closed any closed path (not to form loop) -in fig.2 e.g. conduits 13, 11, 10, 9, 15, 3, 4, 1, 5, 7 or other combination. Branches, which do not

belong to a tree, are links (number of links are X-Y+1). Number of independent loops in network are formed using tree conduits and one of the links conduit). So, number of loops are determined by number of links. In graph, one node is referent (in Fig.

2

referent node is I) and all others are so called dependent nodes. In example from fig.

2

for 1 st guess referent node is I (3). So first row in matrix is for node 2, second for node 3, etc., and last row is for node Y-1 (i.e. XI).

Table captions Table 1

 captions1 Constant node outflow

Table 2

 2 First assumed flowsTable3Combinations for choose of final algebraic sign of second rang correction or higher Table4Calculation of spatial natural gas distributive network of conduits with loops after original Hardy Cross method (example from fig.2)-1 st initial guess Table5Final (calculated) flows in the network for both methods

Table 1

 1

	Constant node outflows
	Consumption a Asymmetric 1 st and 2 nd guess Node m 3 /h m 3 /s	b Symmetric -3 rd guess Node m 3 /h m 3 /s
	I	60	1/60	I'	60	1/60
	II	2100	7/12	II'	2300	23/36
	III	170	17/360	III'	185	37/720
	IV	90	1/40	IV'	90	1/40
	V	200	1/18	V'	185	37/720
	VI	2500	25/36	VI'	2300	23/36
	VII	300	1/12	VII'	290	29/360
	VIII	170	17/360	VIII'	225	1/16
	IX	850	17/72	IX'	850	17/72
	X	280	7/90	X'	225	1/16
	XI Node m 3 /h 280 c Flow into the network 7/90 XI' m 3 /s Node	290 m 3 /h	29/360 m 3 /s
	I	-7000 -1 17/18 I'	-7000 -1 17/18

Table 1 DB ver 3

 1

Table 2

 2

	First assumed flows
	a First assumed flows per conduit -m 3 /s b initial guess of gas flow rate
	conduit	c I	d II	conduit	e III
	1	1/18	1/36	1'	7/36
	2	5/72	1/12	2'	1/9
	3	17/30	47/180	3'	47/180
	4	23/36	5/9	4'	5/6
	5	7/90	1/120	5'	61/240
	6	1/72	31/180	6'	13/144
	7	1/120	7/72	7'	19/72
	8	7/180	1/20	8'	29/144
	9	41/360	7/60	9'	1/48
	10	13/360	7/36	10'	1/12
	11	1/18	4/9	11'	61/240
	12	1/12	1/12	12'	1/9
	13	1/36	5/12	13'	7/36
	14	13/18	1 1/9	14'	5/6
	15	7/18	5/72	15'	1/72

Table 2 DB ver 3

 2

Table 4

 4 Calculation of spatial natural gas distributive network of conduits with loops after original Hardy Cross method (example from fig.2)-1 st initial guess

				Iteration 1
	*	**	D in	L

  ** 29403 *-mark Loop i.e. Contour, **-mark Conduit i.e. Pipe (Main), ***-pressure drop per Loop i.e. Contour (Pa) a calculation of R•Q con

																				3
	I	1 2	0.4064 0.3048	100 100	+1/18 -5/72	114959 -690438	2069265 9942302	+35/123 +35/123		… -6/619±			… …	+251/738 +91/443 d	-15/104 ... -15/104 -4/73=	... ...	+19/97 +5/782
		3	0.1524	100	-17/30	-889949040 1570498307 +35/123		+53/219 ‡	…	-31/773	-15/104 -89/743±	...	-59/194
		4	0.3048	100	+23/36	39193885	61346951	+35/123		…					…	+193/209	-15/104 ...	…	+60/77
					∑ -851330634 1643856824 *** 29178		I		82 , 1 	1643856824 851330634 			123 35	*** 8805		I		295343529 77528086 82 , 1 		104 15
	II	5 6	0.1524 0.3048	100 200	+7/90 -1/72	23969880 -73795	308184165 5313266	+6/619 +6/619		… +1/560 ‡		… …	+37/423 -1/415	+4/73 +4/73	... -27/539±	... ...	+34/239 +1/437 d
		11 12 2	0.1524 0.1524 0.3048	100 100 100	-1/18 -1/12 +5/72	-12993101 -27176838 690438	233875825 326122050 9942302	+6/619 +6/619 +6/619		-43/303  +53/219 ‡ -35/123=	… -114/731  …	-89/441 +2/161 d -91/443 d	+4/73 +4/73 +4/73	-13/392  ... -89/743= -13/392 ‡ +15/104 ‡ ...	-20/111 -3/35 d -5/782
					∑ -15583417 *** 3948	883437609		II			82 , 1 	883437609 15583417 			619 6	*** 10506		II		82 , 1 	1106807745 110372041 			73 4
	III	7 8 9 10 6	0.1524 0.1524 0.3048 0.1524 0.3048	100 100 100 100 200	+1/120 -7/180 +41/360 1698792 411338 -6788773 +13/360 5932191 +1/72 73795	49360570 174568437 14916225 164276062 5313266	-1/560 -1/560 -1/560 -1/560 -1/560		… … -114/731 ‡ -114/731 ‡ -6/619=	... ... ... ... ...	+4/611 -29/713 -31/707 d -9/74 d +1/415	+27/539 ... +27/539 ... +27/539 -13/392  ... ... ... +27/539 -13/392  ... +27/539 -4/73= ...	+29/512 +5/531 d -7/260 -78/745 -1/437 d
					∑ 1327344 *** 1152	408434560		III			408434560 1327344 82 , 1 			560 1	*** 7841		III		82 , 1 	674361136 61479037 			539 27
	IV	3 12 13	0.1524 0.1524 0.1524	100 100 100	+17/30 +1/12 -1/36	889949040 27176838 -3679919	1570498307 -53/219 326122050 -53/219 132477076 -53/219		-35/123= -6/619= …	... +114/731  ...	+31/773 -2/161 d -242/897	+89/743 +15/104  ... +89/743 -4/73± +13/392 ‡ +89/743 ... ...	+59/194 +3/35 d -3/20
		14	0.4064	100	-13/18	-12243919	16953118	-53/219		…					...	-647/671	+89/743 ...	...	-38/45
					∑ 901202040 *** 30020	2046050552		IV		2046050552 901202040 82 , 1 		219 53	*** 15651		IV		82 , 1 	1123546748 244942051 			743 89
	V	15 9 10	0.1524 0.3048 0.1524	200 100 100	+7/18 +41/360 1698792 897059511 +13/360 5932191	2306724456 -114/731 … 14916225 -114/731 -1/560 ‡ 164276062 -114/731 -1/560 ‡			... ... ...	+157/674 -31/707 d -9/74 d	-13/392 ... -13/392 +27/539= ... ... -13/392 +27/539= ...	+177/886 -7/260 -78/745
		11 12	0.1524 0.1524	100 100	-1/18 -1/12	-12993101 -27176838	233875825 326122050	-114/731 +6/619= -114/731 +6/619=		... +53/219 ‡	-89/441 +2/161 d	-13/392 +4/73 = -13/392 +4/73 ±	... -89/743=	-20/111 -3/35 d
					∑ 864520555	3045914618		V			304591461 864520555 82 , 1 		731 114	*** 12786		V		2708713524 163492506 82 , 1 		392 13

*

Table 4 DB ver 3

 4

Table 5

 5 Final (calculated) flows in the network for both methods 1 st and 2 nd guess 3 rd guess d

	conduit D in (m)	a Final flow m 3 /s	m 3 /h	b Velocity check m/s	c conduit D in (m) m 3 /s a Final flow	m 3 /h	b Velocity check m/s
	1	0.4064	334/979 1228.19	2.63	1'	0.1524 85/421	726.84	11.07
	2 3	0.3048 0.1524	13/129 82/539	362.80 547.68	1.38 8.34	2' 3'	0.1524 1/29 0.1524 81/329	124.14 886.32	1.89 13.50 d
	4	0.3048	551/596 3328.19 12.67	4'	0.3048 639/760	3026.84 11.52
	5	0.1524	130/673 695.39	10.59	5'	0.1524 180/973	665.98	10.14
	6	0.3048	1/71	50.73	0.19	6'	0.3048 0	0	0
	7	0.1524	9/94	344.66	5.25	7'	0.1524 40/383	375.98	5.73
	8	0.1524	18/371	174.66	2.66	8'	0.3048 32/763	150.98	0.57
	9	0.3048	22/687	115.28	0.44	9'	0.3048 32/763	150.98	0.57
	10	0.1524	28/255	395.28	6.02	10' 0.1524 40/383	375.98	5.73
	11	0.1524	106/611 624.55	9.51	11' 0.1524 180/973	665.98	10.14
	12	0.1524	17/235	260.43	3.97	12' 0.1524 1/29	124.14	1.89
	13	0.1524	76/485	564.13	8.59	13' 0.1524 85/421	726.84	11.07
	14	0.4064	223/262 3064.13 6.56	14' 0.3048 639/760	3026.84 11.52
	15	0.1524	7/45	560.05	8.53	15' 0.1524 58/381	548.03	8.35

Table 5 DB ver 3

 5

for example networks shown in Fig.2

see conduit 3' in Table

Table 3

Combinations for choose of final algebraic sign of second rang correction or higher

indicates the direction of the conduit flow for particular loop (contour) -a plus sign denotes clockwise flow in the conduit within the loop; a minus sign, counterclockwise b second rang correction is from the second contour to which a conduit belongs, third rang correction is from the third contour to which a conduit belongs, etc. c the upper plus or minus sign shown indicates direction of flown in that conduit in these two contours and is obtained from Q for previous iteration. The upper sign is the same as the sign in front Q if the flow direction in the particular contour under consideration, and opposite if it does not. d the lower sign is copied from the primary contour for this correction (sign from the contour where this correction is first, sign preceding the first iteration from adjacent contour for the conduit taken into consideration) e logical zero is equivalent to sign (-) minus f logical one is equivalent to sign (+) plus