Iterative Methods for Looped Network Pipeline Calculation

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.

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.

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 (Farshad et al. 2001, Kumar 2010.
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 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 4 networks, gas density can be treated as constant, which means that gas can be treated as incompressible fluid (Pretorius et al. 2008), 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" (Sukharev et al 2005) 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. 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   5 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):

Liquid flow in pipes
Then pressure drop (Ekinci and Konak 2009) can be found very easily using Darcy-Weisbach equation (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): Colebrook"s equation is suitable especially for flow through steel pipes (Colebrook 1939).
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 andPinho 2007).
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.

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 (Mah and Shacham 1978). 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 (Cross 1936). 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.  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65   7 Contemporary with Hardy Cross, soviet author V.G. Lobačev (Latišenkov and Lobačev 1956) developed very similar method compared to original Hardy Cross method. Modified Hardy Cross method proposed Epp and Fowler (1970) which considers entire system simultaneously is also sort of loop method. Node-Loop method proposed by Wood and Charles (1972) and later improved by Wood and Rayes (1981) 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 Shamir and Howard (1968) 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 Epp and Fowler (1970) 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 8 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 (Gay and Middleton 1971). It would be equally satisfactory to generate the initial distribution within the computer program. Initial flow pattern for our example network is shown in the figure 2. Loops are also defined in the figure 2.   Corfield et al (1974) from Gas Engineers Handbook.

Hardy Cross method for the pipeline calculation
Hardy Cross, American engineer, developed method in 1936, later named after him (Cross 1936 10 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) 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:  In presented example loop I begins and ends in node II via pipes 3, 4, and 7.

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. Epp and Fowler (1970) 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). III   II   I   III   II   I   III   III   II   III   I   III   III   II   II   II   I   II   III   I   II   I   I   I Presented matrix is symmetrical; for example (13): )   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 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 (Brkić 2009).
Using numerical values from table 2 for the first iteration of calculation of water network matrix equation can be written (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 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 14 contour in M.M. Andrijašev method as in Hardy Cross approach. Andrijašev"s contour will be marked with special sign ( O ). Matrix formulation of this method for example gas network can be written as (15): Here has to be very careful because non-diagonal terms are not always negative as in modified Hardy Cross method (16):

Methods based on solution of the node equations
Nodes are sometimes also referred to as junctions, points or vertices.
Shamir and Howard (1968) 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) Similar, can be done for Darcy-Weisbach equation (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): 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). Calculation of the network using non-improved Node method is not shown here.

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.
Wood and Charles (1972)  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   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 Node VI will be omitted from the node matrix to assure linear independency of the rows as 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): 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) Boulos et al. (2006) can be recommended for further reading. In this book, authors instead to omit one node in the node matrix to preserve 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 23 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. Mathews and Köhler (1995) 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).

Figure 5. Comparisons of convergence for presented methods (gas network)
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): Velocity of gaseous fluids depends on the pressure in pipe since they are compressible (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 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   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 25 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 (Brkić 2009). 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 Mah (1974), Mah and Lin (1980) and Altman and Boulos (1995).