An efficient iterative method for looped pipe network hydraulics

Original and improved version of the Hardy Cross iterative method with related modifications are today widely used for calculation of fluid flow through conduits in loops-like distribution networks of pipes with known node fluid consumptions. Fluid in these networks is usually natural gas for distribution in the municipalities, water in waterworks or hot water in district heating system, air in the case of ventilation systems in buildings or mines, etc. Since, the resistances in these networks depend of flow, problem is not linear like in electrical circuits, and iterative procedure must be used. In both version of the Hardy Cross method, in original and in the improved one, initial results of calculation in iteration procedure is not flow, but rather the correction of flow. Unfortunately, these corrections should be added to or subtracted from a flow calculated in previous iteration according to complicate algebraic rules. After the here presented node-loop method, final results in each of the iterations is flow directly rather than flow correction. In that way complex algebraic scheme for sign of flow correction is avoided, while the final results still remain unchanged. Numbers of required iterations for the same results are comparable with the improved Hardy Cross method.


Introduction
Since, the resistances in a network of pipes for distribution of fluids depend on flow, problem is not linear like in electric circuits and iterative procedure must be used to calculate distribution of fluid flow through pipes and distribution of pressure in the network.Usually, in a hydraulic network of pipes, consumption of fluid assigned to each node is known and stays unchanged during computation.This is also the case for the inputs in network which are also assigned to nodes and which also do not change during calculation.Further, in order to calculate flow and pressure distribution in the network of pipes, first of all, initial flow pattern through pipes in the network has to be assigned to satisfy first Kirchhoff law for each node.This means to satisfy material balance of fluid moved through network.During iterative cycles of calculation, this flow distribution will changes in order to conform second prerequisite condition govern by the second Kirchhoff law, i.e. to satisfy energy balance in each closed conduit formed by pipes in the network.In hydraulic network this energy balance is usually expressed through pressure or some of the functions in which pressure exist.While the first Kirchhoff law has to be satisfied in all iterations for each node in the network, the second Kirchhoff law has to be satisfied for each closed conduit at the end of calculation.
Usually, such as in Hardy Cross method [1] and related improved version [2], result of iterative calculation of flow distribution pattern in a hydraulic network is correction of flow [1][2][3].This correction of flow has to be added to flow calculated in the previous iteration using complex algebraic rules [3,4].This intermediate step will be eliminated, using procedure that will be shown in this paper.In that way, flow will be directly calculated in all iteration for each pipe.
All methods from this paper assume equilibrium between 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 drop where inertia force is important.Gas flow in a municipal distribution network [5], air flow in a ventilation system in buildings and mines [6], and of course water flow in waterworks [7] or district heating systems [8] and cooling systems [8] can be treated as incompressible flow since the pressure drop in these kinds of networks are minor even to compress significantly natural gas or air.The same applies to pipelines for distribution of mixed natural gas and hydrogen [9].

Loop-oriented methods; Original and improved Hardy Cross method
The Hardy Cross method [1] introduced in 1936 is the first useful procedure for the calculation of flow distribution in looped networks of pipes.Further step was made by introduction of the modification in the original Hardy Cross method in 1970 by Epp and Fowler [2].The original Hardy Cross method [1] as a sort of single adjustment method, first of all, as an intermediate step in calculation, determines correction of flow for each loop independently and then applies this corrections to compute new flow in each conduit.It is not efficient as the improved Hardy Cross method [2,3] that considers entire system simultaneously.The improved Hardy Cross method [2], still firstly as an intermediate step, determines corrections for each loop but treated all network system simultaneously, and then applies this correction to compute new flow in each conduit such as in the original version [1].It is more efficient, but also intermediate step in calculation is not eliminated.More than thirty years had to pass by before the introduction of the modification by Epp and Fowler [2] only because of matrix calculation.While use of matrix form in the original Hardy Cross method is not mandatory [1], for the improved version it is [2].In the original paper of Hardy Cross from 1936 [1], problem is not solved using any kind of matrix calculation (but also this approach can be expressed using matrix calculation with no affects on final results [7]).

Node-oriented methods
Two years before modification of the original Hardy Cross method, Shamir and Howard in 1968 [10] reformulated original method to solve node equations and not any more loop equations like in the original Hardy Cross method [1].The node equations expressed in the node method in terms of unknown pressure in nodes [11].Methods based on node equations are less reliable which means that the single adjustment methods based on idea from the original Hardy Cross method (but here adjust for nodes) must be employed with caution.Idea for these node-oriented methods is simple knowing principle of loop-oriented method developed by Hardy Cross [1].In a loop-oriented method, energy distribution for all closed paths in a network governs by the second Kirchhoff law will be always satisfied, while material balance for all nodes in a network governs by the first Kirchhoff law will be balanced in an iterative procedure.Similar principle applies as in the original Hardy Cross method, but only with opposite approach.Still, as an intermediate step, correction of pressure has to be calculated [12][13][14] (in the original method by Hardy Cross this is correction of flow [15][16][17]), and then after that, pressure as a final result of iteration has to be calculated using complex algebraic rules.Pressure can be expressed in different quantities, such lengths of water elevation or similar.

Node-loop oriented method
After the development of the loop-oriented and node-oriented methods, and after introduction of matrix calculus, all necessary tools are available, i.e. matrix form of loop method and matrix form of node method to unite both, the loop and the node equations in matrix form which has a result completely new and innovated method [18,19].This transformation makes possible direct calculation of final flow in each of the iterations, and not the correction of flow like in methods mentioned before (Figure 1).Unfortunately, as already explained these corrections of flow calculated after previous methods should be added to or subtracted from flow (or pressure in the node method) calculated in previous iterations according to complicated algebraic rules [3].Brkić (2009) [3], Gay & Middleton (1971) [38] See Brkić (2009) [3], Epp & Fowler (1970) [2] See Brkić (2009)  So, the main strength of the node-loop method introduced in 1972 by Wood and Charles [18] for waterworks calculation does not reflect in noticeably reduced number of iteration compared to the modified Hardy Cross method.Main advantage of this method is in the capability to solve directly the pipe flow rate rather than flow correction.The method uses a linear head loss term which allows a network of n pipes to be described by a set of n linear equations which can be solved simultaneously for the flow distribution.Wood and Rayes in 1981 introduced improvement in the node-loop method [19].Here will be shown improved version of this method rearranged for gas flow and for water flow in terms of pressure distribution rather than head distribution (which quantity is express in length; such as elevation of water).

Some literary overview of the existed methods for calculation of flow distribution in a looped network of pipes
Excellent example of calculation of looped natural gas distribution network after original Hardy Cross method can be found in Gas Engineers Handbook from 1974 [4].Already mentioned algebraic rules for correction of flow calculated as an intermediate step in iterative procedure that can be used for both versions of Hardy Cross method can be found in this reference book [4] (and also for the node-oriented method but where correction of pressure is calculated as an intermediate step rather than correction of flow).This algebraic rules were further additionally and developed in Brkić [3].Same spatial gas network as shown in Brkić [3] will be also used here for calculation after the node-loop method.Same topology of the network with same diameter will be used here for calculation of water flow as comparisons of the results obtained for liquid flow.
Excellent book in this issue, but only for waterworks calculation by Boulos et al. [20] can be recommended for further reading.In this book, unfortunately, the Hazen-Williams equation, an obsolete relation is used to correlate only water flow, pressure drops in pipes and hydraulics frictions.
Further, for details on natural ventilation airflow networks one can consult paper of Aynsley [6].There is no space here to calculate separately air ventilation network, but readers interested in this matter can make this in a very effective way, according to natural gas and water flow calculation shown in this paper.Specific details on airflow resistances are also given in Aynsley [6].
Also, Todini and Pilati [21] for water networks and Hamam and Brameller [22] for gas networks wrote conservation of energy for each pipe and as result beside of flow correction in each pipe, pressure drop also can be simultaneously calculated.This method is also known as hybrid or gradient approach.Some comparisons of available methods for pipeline network calculations can be found in Mah [23], Mah and Shacham [24], Mah and Lin [25], etc.To compare calculation of water networks using the Hazen-Williams equation and approach with pseudo-loops consult book of Boulous et al [20].Lopes [26] also deals with the program for the Hardy Cross solution of the piping networks.Shown kind of problems today can be solved very easily using MS Excel [27,28].
The first computer solutions of network problems were done on analog computers where electrical elements are used to simulate pipe networks [29].Today, this approach is obsolete.Also, today natural gas is mostly distributed in cities, but earlier it was gas derived from coal [30].

Hydraulics resistance of a single pipe
Source-issue that cause problem with the calculation of hydraulic networks is non-constant value of hydraulic resistance when fluid convey through pipe.On the other hand, electrical resistance of a wire or a resistor has a constant value which has a consequence, non-iterative calculation of electrical circuits.To establish relation between flow rate of natural gas through a single pipe and related pressure drop, the Renouard equation for gas flow will be used and in that case (1) [25].Using that approach, resistance will not be calculated at all since Renouard's equation relates pressure and flow rate using other properties, parameters and quantities to connect these two variables.On the other hand, for the calculation of hydraulic resistance in a single pipe, well known Colebrook equation will be used [26] (which is also iterative and which caused also some problems [33][34][35]) where pressure drop is calculated using Darcy-Weisbach equation.Finally, for calculation of air-flow through ventilation system, one can consult Aynsley [6], as already mentioned before.
The Hazen-Williams equation, which is used in here recommended book of Boulos et al. [20], is useless for calculation of gas flow.Introduced in the early 1900s, the Hazen-Williams equation determines pipe friction head loss for water, requiring a single roughness coefficient (roughness is also very important parameter also in Darcy-Weisbach scheme for calculation [36]).Unfortunately even for water it may produce errors as high as ±40% when applied outside a limited and somewhat controversial range of the Reynolds numbers, pipe diameters and coefficients.Not only inaccurate the Hazen-Williams equation is conceptually incorrect [37].
In this paper the focus is on pipes, and other parts of systems are not examined.Furthermore, in a water or gas distribution system, the pipe friction head losses usually predominate and other minor losses can be ordinarily neglected without serious errors [38][39][40][41].

Topology of looped pipe system
First of all, maximal consumption per each node including one or more inlet nodes has to be determined (red in Figure 2).These parameters are looked up during the calculation.Further, initial guess of flow per conduits has to be assigned to satisfy first Kirchhoff's law and in that way chosen values are to be used for first iteration [3].Final flows do not depend on first assumed flows per pipes (countless initial flow pattern can satisfy first Kirchhoff's law and all of them equally can be used with the same final results [3,38]).After the iteration procedure is completed, and if the value of gas or water flow velocity for all conduits are bellow standard values, calculated flows become flow distribution per pipes for maximal possible consumptions per nodes.Further, pressure per all nodes (can be heads in case of water) can be calculated.Whole network can be supplied by gas or water 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 i.e. water consumption of households.Disposal of households is along the network's conduits, and only their consumption is to be assigned to nodes.Main purpose of the method is to calculate flow pattern per pipes and pressure pattern per nodes for the maximal load of the network 1 .First assumed flow pattern has to be chosen to satisfy first Kirchhoff's law (continuity of flow) which means that algebraic sum of flows per each node must be zero exactly.On the other hand, second Kirchhoff's law (continuity of potential), which means that algebraic sum of pressure drops per each contour must be approximately zero at the end of iterative procedure.Procedure can be interrupted when algebraic sum per all nodes become approximately zero, or when flows per pipes are not changed in calculation after two successive iterations.
One spatial fluid distribution network of pipelines will be examined as example (Figure 2).Polyethylene pipes (PVC) are used in the example shown in this paper.

Figure 2. Spatial gas/water distribution network with loops -example
The first step in solving a problem is to make a network map showing pipe sizes and lengths, connections between pipes (nodes), and sources of supply.For convenience in locating pipes, assign each loop and each pipe a code number.Some of the pipes are mutual to one loop and some to two 1 Problem can be treated as inverse, i.e. flow per pipes assigned in the first iteration is not only initial pattern see (17).This flow pattern is not variable in further calculation.Instead of flows per pipe which are now constants, pipe diameters become variables, and according to this approach, optimized pipes' diameters in the network are the final result of calculation (see section 8 of this paper).or even three contours (i.e.pipe 12 belongs to the loops II, IV, V).Special cases may occur in which two pipes cross each other but are not connected (like pipes 6 and 15), resulting in certain pipes being common to three or more loops.The distribution network then becomes three-dimensional (rare for gas with exception of maybe some chemical engineering facilities, water networks or district heating system, and usually for airflow networks).For example, loop V consists of conduits 15, 9, 10, via 11, and 12. Gas/water flow into the network from a source on the left side is 7000 m 3 /h, and points of delivery are at junctions of pipes (nodes), with the red arrows pointing to volumes delivered (node consumption).Summation of these deliveries equals 7000 m 3 /h.Assumed gas flows and their directions are indicated by black arrows near the pipes (Figure 2).

Topology equations for the observed looped network of pipes
After the network map with its pipes and loop numbers and delivery and supply data has been prepared, mathematical description of the network can be done.To introduce matrix form in calculation, it is necessary to represent distribution network from Figure 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 where in Figure 2, X = 15 and Y = 11).Graph with n nodes (in our case 11) has Y-1 independent nodes (in our case 10) and X-Y+1 independent loops (in our case 5).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.In graph, one node is referent 2 (in Figure 2 referent node is XI) and all others are so called dependent nodes.

Loop equations
The Renouard equation (1) will be used for calculation of pressure drop in pipes in the case of natural gas distribution [31]. (1) Regarding to Renouard formula (1) one 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 conduit.This means that is not actually pressure drop in spite of the same unit of measurement, i.e. same unit is used as for pressure (Pa).Parameter rather can be noted as pseudo-pressure drop.Fact that when this consecutive means that also is very useful for calculation of gas pipeline with loops.So, notation for pseudo-pressure drop is ambiguous [3] (only or with appropriate index should be used instead of ).
First derivative of previous relation where the flow is treated as variable is (2): (2) 2 In other approach, with no referent node, one pseudo-loop must be introduced [39].This is very complicated and should be avoided.
The Colebrook-White equation (3) will be used for calculation the Darcy's friction factor in the case of water distribution [32].The Colebrook-White equation is implicit in friction factor, and here it is solved using MS Excel.
(3) Friction factor calculated after Colebrook's relation will be incorporated into the Darcy-Weisbach relation to calculate pressure drop in water network (4).
(4) Similar as for the gas lines, first derivate of previous relation where the flow is treated as variable is (5): Then, according to previous, for the gas network from figure 2, set of loop equation can be written as ( 6): (6) Previous relations can be noted in matrix form as (7): , Or for waterworks or district heating systems from figure 2 can be noted as (8): i.e. in matrix form for water distribution (9)., In the left matrix of the relations ( 7) and ( 9) rows represent loops and columns represent pipes.These relations are matrix reformulation of the second Kirchhoff's law.The sign for the term relates if the assumed flow is clockwise (1) or counter-clockwise (-1) relative to the loop.

Node equations
For all nodes in the network from figure 2, relations after the first Kirchhoff's law can be noted as (10):   (10) Or in matrix form as (11) where in the first matrix rows represents nodes excluding referent node 3 .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.No information on the topology in that way will be lost [22].
First row corresponds to the first node, etc.Last row is for node 10 from figure 2, since the node 11 is chosen to be referent node and therefore must be omitted from the matrix.For example, node 1 has connection with other nodes via pipes 3, 4 and 14, and for first assumed flow pattern, all flows are from node 1 via connected pipes to other nodes.Therefore, terms 3, 4, and 14 in first row are -1.Other pipes are not connected with node 1, and therefore all other terms in the first row of node matrix are 0.
Note that there is no difference in cases of water apropos gas calculation when the node equations are observed.

Network calculation according to the node-loop method
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.This method will be examined in details for the network shown in Figure 2.This network will be treated as natural gas network in the sections 7.1 and as water network in 7.2.This approach also gives good insight into the differences which can be occurred in the cases of distribution of liquids apropos gaseous fluids.

The node-loop calculation of gas networks
First iteration for the gas calculation for the network from Figure 2 is shown in Table 1.If sign of calculated flow is negative, this means that flow direction from previous iteration must be changed, otherwise, sing must be remained unchanged.In Table 1, loop and the pipes numbers are listed in the first and the second column, respectively.Pipe length expressed in meters is listed in the third column, and assumed gas flow in each pipe expressed in m 3 /s is listed in the fourth column.The 1 or -1 in fifth column indicates sing preceding flow in the fourth column.The plus or minus preceding the flow, Q, indicates the direction of the pipe flow for the particular loop.A plus sign denotes clockwise flow in the pipe within the loop, a minus sign counterclockwise.All these assumption will not be changed also in the case of waterworks or district heating system calculation.13) and ( 14) To introduce matrix calculation, the node-loop matrix [NL], matrix of calculated flow in observed iteration [Q], and [V] matrix in the right side of (12) will be defined.
First ten rows in the NL (13) matrix are from node matrix (11), and next five rows are loop matrix (7 and 9).These five rows from the loop matrix are multiplied by first derivate of pressure drop function (2) from Table 1 for gas 4 (column F').
First ten rows in matrix [V] are node consumption 5 , and the rest five terms are from Table 1 (14).
Sign minus in front of some term means that sing preceding this term from the previous iteration must be changed.
Five iterations are enough for the calculation of gas network from Figure 2. Calculated flows for these first five iterations will be listed in Table 2. a First assumed flows per pipes chosen after the first Kirchhoff's law (black letters in figure 2) b Values in iterations 5 are equal as in iteration 4, stopping criterion is fulfilled c Gas velocity (10-15 m/s recommended); υ=(4·p·Q)/(δ 2 ·π); where p= pn/pa=0.25 and pa is absolute pressure of gas in the pipeline, here pa=400 000 Pa, and pn=normal pressure ~100 000 Pa, p=400 000 Pa/100 000 Pa=1/4 Flow direction is changed in pipe 2, 6, 8, 9 and 10 (opposite than assumed in first assumed flows).Note that the velocities in the last column of Table 2 are listed.Gas pressure in the network is circa 4x10 5 Pa abs.Flow velocity per pipes is not balanced, somewhere is too small somewhere is too high.Whole problem can be treaded now as inverse by fixing the flows per pipes and optimized pipe diameters as noted in section 4.This can be done using here presented node-loop method, Hardy Cross or similar available methods (note that different pressure in the gas apropos water network causes different values of speed of gas compared to speed of water; last column in Table 2 and 4, respectively).

The node-loop calculation of waterworks or district heating systems
Similar as for gas networks, network from Figure 2 will be used for water distribution calculation (Table 3).The calculated flows listed in Table 4 are slightly different than for the gas flow calculation.

A note on optimization problem
Renouard formula (1) for condition in gas distribution networks assumes a constant density of a fluid within the conduits.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.Assumption of gas incompressibility means that it is compressed and forced to convey through conduits, but inside the pipeline system pressure drop of already compressed gas is minor 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 typical distribution gas network is 4x10 5 Pa abs i.e. 3x10 5 Pa gauge and accordingly volume of gas is decreased four times compared to volume of gas at normal (standard) conditions.But operate pressure for gas distribution network can be lower (this case is valid for network in paper of Brkić [3]).This is not typical for natural gas distributive networks.This was common practice in obsolete systems for distribution of city gas derived from coal [30].So, flow in Renourad formula (1) adjusted for natural gas is usually expressed in normal (standard) conditions.Consequence is that if flows in previous paper of Brkić [3] are expressed in their real (compressed) values and if these real values are numerically equalized with values expressed for normal (standard) conditions, this means that operate pressure in gas network is normal (standard).Otherwise, velocities in previous paper of Brkić [3] have to be corrected.Velocities in previous paper of Brkić [3] are calculated to be comparable with the procedure shown in Manojlović et al [40] where calculation of gas distribution network in Serbian town Kragujevac is discussed.In Manojlović et al [40], flows are expressed in their real values and not for normal or standard conditions of pressure as common practice is or this network is calculated to work with lower pressure typical for gasses derived from coal.Second assumption can be rejected as less possible because in the part of Serbia south of rivers Sava and Danube where Kragujevac is situated, such gas was never used and especially not in 1990's.Some comments about that issue was also shown in Brkić [5].So, to avoid any further ambiguity, conclusion is that all flows previous paper of Brkić [3] are expressed in their real (compressed) values while operate pressure at the inputs of shown networks is normal (standard).
If these values of flows are noted for normal (standard) conditions of pressure as common practice is (Table 2), while operate pressure is 4x10 5 Pa abs i.e. 3x10 5 Pa gauge, velocities of gas are different than those in previous paper of Brkić [3] while flows remain unchanged.
Velocities in Table 2 are calculated using ( 16): Now, for such values of flows, diameters of conduits are too large and in such case Hardy Cross method [1] as well as improved Hardy Cross method [2,3] can be used for optimization of diameters of conduits shown in Figure 2. In a problem of optimization of pipe diameters, in Renouard formula (1), flow is not any more treated as variable (17) while correction Δ is now correction of diameters.(17) Ambiguity related to pressure conditions in a gas distributive network can cause very different and large consequences in an interpretation of calculated results.
Similar analogy regarding to water networks is clear (18): Diameters of conduits in presented gas pipeline should be optimized while diameters in water network are in an accepted tolerance.

Conclusions
Here presented the node-loop method is powerful numerical procedure for calculation of flows or diameters as inverse problems in looped fluid distribution networks.Main advantages is that flow in each pipe can be calculated directly, which is not possible after Hardy Cross and improved Hardy Cross methods (Figure 3).Similar numbers of iterations are necessary to achieve demanded accuracy in calculation as in the modified Hardy Cross method6 (Figure 4).The hydraulic computations involved in designing water or gas distribution systems can be only approximated, as it is impossible to consider all the factors affecting loss of head in a complicated network of pipes.
The here presented methods can be easily readapted for detection of a position of leakage in a pipe network [42,43].

Figure 3 .
Figure 3. Main conceptual difference between the Hardy Cross method (original and improved) and the node-loop method

Figure 4 .
Figure 4. Comparison of the convergence performances for the Hardy Cross methods, original and improved, and the node-loop method

Table 1 .
Node-loop analysis for the gas network from Figure1.
b letters used in (

Table 2 .
First five iteration for gas network from Figure1 -example

Table 3 .
(5)e-loop analysis for the water network from Figure2.from Figure2but expressed in m 3 /s, b letters used in (13) and (14), c Reynolds number; dynamic water viscosity 0.00089 Pa•s, d Relative roughness; absolute roughness ε=0.00002 m for PVC pipes, e Friction factor (3) calculated using MS Excel, f Pressure drop in pipe (4), g see(5) a

Table 4 .
First seven iteration for water network from Figure2-example, Flow in m 3 /h a First assumed flows per pipes chosen after the first Kirchhoff's law (black letters in figure2), b Values in iterations 7 are equal as in iteration 6, stopping criterion is fulfilled