Hardy Cross method for pipe networks

Hardy Cross originally proposed a method for analysis of flow in networks of conduits or conductors in 1936. His method was the first really useful engineering method in the field of pipe network calculation. Only electrical analogs of hydraulic networks were used before the Hardy Cross method. A problem with the flow resistance versus the electrical resistance makes these electrical analog methods obsolete. The method by Hardy Cross is taught extensively at faculties and it still remains an important tool for analysis of looped pipe systems. Engineers today mostly use a modified Hardy Cross method which threats the whole looped network of pipes simultaneously (use of these methods without computers is practically impossible). A method from the Russian practice published during 1930s, which is similar to the Hardy Cross method, is described, too. Some notes from the life of Hardy Cross are also shown. Finally, an improved version of the Hardy Cross method, which significantly reduces number of iterations, is presented and discussed. Also we tested multi-point iterative methods which can be used as substitution for the Newton-Raphson approach used by Hardy Cross, but this approach didn’t reduce number of required iterations to reach the final balanced solution. Although, many new models have been developed since the time of Hardy Cross, main purpose of this paper is to illustrate the very beginning of modeling of gas and water pipe networks or ventilation systems.


Introduction
Hardy Cross solved the problem of distribution of flow in networks of pipes in his article "Analysis of flow in networks of conduits or conductors" [1] published on November 13th 1936.
Networks of pipes are nonlinear systems since the relation between flow and pressure is not linear.On the contrary, the relation between current and voltage in electrical networks with regular resistors is governed by the linear Ohm's law.Electrical circuits with diodes as well as hydraulic networks are nonlinear systems where resistance depends on current and voltage i.e. on flow and pressure, respectively [2].
Distribution of flow in a network of pipes depends on the known inputs and consumptions at all nodes, on the given geometry of pipes and topology of network.A stable state of flow in a network must satisfy Kirchhoff's laws, which represent statements of the conservation of mass and energy.Although indefinite number of flow distributions, which satisfy the conservation of mass is possible in theory, only one distribution from this set satisfies also the conservation of energy for all closed paths formed by pipes in the network.This state is unique for the given network and in and out flows [3].
Since the relation between flow and pressure is not linear, Hardy Cross used a relation between an increment of flow and an increment of pressure, which relation is linear for a given quantity of flow.If, however, the increments are fairly large, this linear relation is somewhat in error like for gas compressible flow.But if the pressure drop in pipes is minor like in municipality network for natural gas distribution, Hardy Cross method can be used without significant errors [4][5][6].It can be used also for water pipe networks and for ventilation systems [7] (related formulation is in Appendix A of this paper).
Hardy Cross method is an iterative method, i.e. the method of successive corrections [4].Lobačev and Andrijašev in 1930s, writing in Russian, offered similar methods [8,9].Probably because of the language and the political situation in Soviet Russia, Hardy Cross was not aware of Lobačev and Andrijašev contributions.
Today, engineers use the mostly improved version of Hardy Cross method (ΔQ method), which threats the whole looped network of pipes simultaneously [10].
As novel approach for the first time presented here, we tested multi-point iterative methods [11,12] which can be used as substitution for the Newton-Raphson approach used by Hardy Cross, but this approach didn't reduce number of required iterations to reach the final balanced solution.
One example of the pipe network for distribution of gas is analyzed using the original Hardy Cross method [1]; in Section 3.1, its related equivalent from Russian literature [8,9]; in Section 3.2, the improved version of the Hardy Cross method [10,13]; in Section 3.3, and finally the approach which use multi-point iterative methods instead of the commonly used Newton-Raphson method in Section 3.4.

Topology of network
The first step in solving a pipe network problem is to make a network map showing pipe diameters, lengths and connections between pipes (nodes).Sources of natural gas supply and consumption rates have to be assigned to nodes.For convenience in locating pipes, to the each pipe and the closed loop of pipes there are assigned code numbers (represented by roman numbers for loops in Figure 1).Pipes on the network periphery are common to the one loop and those in the network interior are common to two loops.Figure 1 is an example of pipe network for distribution of natural gas for consumption in households.

Figure 1. Network of pipes for natural gas distribution for domestic consumption
The next step is to write initial gas flow distribution through pipes in the network.This distribution must be done according to the first Kirchhoff's law.The choice of initial flows is not critical and the criterion should satisfy the first Kirchhoff's law for every node in the network [3].Total gas flow arriving at a node equals total gas flow that leaves the node.The same conservation law is also valid for the whole network in total (except of gas input and output nodes that cannot be changed during calculations; see consumption nodes in Figure 1).Sum of pseudo-pressure drops along any closed path must be approximately zero for the network to be in balance according to the second Kirchhoff's law.In this paper, the flow distribution, which satisfies both Kirchhoff's laws will be calculated using the Hardy Cross iterative method.

Hydraulic model
Renouard formula; Eq. ( 1) fits best a natural gas distribution system built with polyethylene (PVC) pipes [14,15].Computed pressure drops are always less than actual drop since the maximal consumption is occurred only during extremely severe winter days [16].
Where f is function of pressure, ρr is relative gas density (dimensionless); here ρr=0.64,L is length of pipe (m), D is diameter of pipe (m), Q is flow (m 3 /s), ΔQ is flow correction (m 3 /s), and p is pressure (Pa).
As shown in Appendix A of this paper, another formulas are used in the case of waterworks systems [17,18] or ventilation networks [7].
Regarding to the Renouard formula; Eq. ( 1) one has to be careful since the pressure drop function, f, does not relate pressure drop but actually difference of the quadratic pressure at the input and the output of pipe.This means that is not actually pressure drop in spite of the same unit of measurement, i.e. the same unit is used as for pressure (Pa).Parameter can be noted as pseudo-pressure drop.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.Operate pressure for typical distribution gas network is 4x10 5 Pa abs i.e. 3x10 5 Pa gauge and accordingly volume of gas decreases four times compared to the volume of gas at normal (or standard) conditions.Pressure in the Renouard formula is for normal (standard) conditions.
First derivative f' of the Renouard relation; Eq. ( 2) where the flow is treated as variable is used in the Hardy Cross method.

 
First assumed gas flow in each pipe is listed in the third column of Table 1.The plus or minus preceding flow indicates the direction of the flow through pipe for the particular loop [13,19].A plus sign denotes counterclockwise flow in the pipe within the loop while minus sign, clockwise.Loop direction can be chosen to be clockwise or counterclockwise (in Figure 1 all loops are counterclockwise).

Hardy Cross method; different versions
Here will be presented Hardy Cross method; original approach in Section 3.1.,version of the Hardy Cross method from Russian practice in Section 3.2., Modified Hardy Cross method in Section 3.3, and finally the method that uses multi-point iterative procedures instead the Newton-Raphson method and which can be implemented in all already mentioned methods.

Hardy Cross method; original approach
Pressure drop function for each pipe is listed in Table 1 (for initial flow pattern in the forth column).Sign in front of the pressure drop function shown in forth column is the same as for flow from the observed iteration.In the fifth column of Table 1 are listed the first derivatives of pressure drop function where flow is treated as variable.Column of the function of pressure drops is added algebraically while column of the first derivatives is added arithmetically for each loop.Flow correction ΔQ has to be computed for each loop x; Eq. (3).
For the network from Figure 1, the flow corrections for the first iteration in each loop can be calculated using Eq. ( 4).
In the second iteration, the calculated correction ΔQ has to be added algebraically to the assumed gas flow (first initial flow pattern).Further, the calculated correction ΔQ has to be subtracted algebraically from the gas flow computed in previous iteration.This means that the algebraic operation for the first correction is opposite of its sign, i.e. add when the sign is minus and opposite.A pipe common to two loops receives two corrections simultaneously.First correction is from the particular loop under consideration while the second one is from the adjacent loop which observed pipe also belong to.
The upper sign after second correction in Table 1 is plus if the flow direction in mutual pipe coincides with assumed orientation of adjacent loop, and minus if it does not (Figure 2).Lower sign is the sign in front of correction ΔQ calculated for adjacent loop (Figure 2).
Details for signs of corrections can be seen in Brkić [13] and Corfield et al. [19].
The algebraic operation for second correction should be the opposite of its lower sign when its upper sign is the same as the sign in front of flow Q, and as indicated by its lower sign when its upper sign is opposite to the sign in front of flow Q.
The calculation procedure is repeated until the net algebraic sum of pressure functions around each loop is as close to zero as the degree of precision desired demands.This also means that the calculated corrections of flow and change in calculated flow between two successive iterations is approximately zero.The pipe network is then in approximate balance and calculation after the Hardy Cross can be terminated.
In the original Hardy Cross method, the corrections for the first iteration are: ΔQI=1575448179.8

Version of Hardy Cross method from Russian practice
As mentioned in introduction, two Russian authors, Lobačev [8] and Andrijašev [9], proposed similar method as Hardy Cross [1].These two methods are also from the 1930's.It is not clear if Hardy Cross had been aware of the contribution of these two authors from Soviet Russia and vice versa, but most probably answer to this question is no, for both sides.The main difference between Hardy Cross and Andrijašev method is that in method of Andrijašev contours can be defined to include few loops.This strategy only complicated situation while the number of required iteration remains unchanged.
Further on Andrijašev method can be seen from the example in paper of Brkić [3].
Here it will be shown the method of Lobačev in more details.
In Hardy Cross method, influence of adjacent loops is neglected.Lobačev method takes into consideration this influence; Eq. ( 5): In previous system of equations; Eq. ( 5), sings in front of terms from the left side of equal sign have to be determined (this is much more complex than in the Hardy Cross method).So, in Lobačev method if (Σf)x>0 then sign in front of (Σ|f'|)x has to be positive, and opposite (for the first iteration this can be seen in table 1; fI=+1575448179.8>0,fII=-8424412.4<0,fIII=-170493836.7<0,fIV=-749453158.7<0,fV=-325325177.5<0).Sign for other terms (these terms are sufficient in the Hardy Cross method) will be determined using further rules and scheme from Figure 3. From Figure 3 it can be seen that if (Σf)x>0 and if the assumed flow coincides with the loop direction, then the sing of flow in adjacent pipe is negative and if the flow does not coincide with the loop direction then the sing of flow in the adjacent pipe is positive.And opposite, if (Σf)x<0 and if the assumed flow coincides with the loop direction, then the sing of flow in the adjacent pipe is positive and if the flow does not coincide with the loop direction then the sing of flow in the adjacent pipe is negative.This procedure determines the sign in the front of flow correction (ΔQ) which are shown in Figure 3 with black letters (also in Figure 4 for our example network of pipes).If (Σf)x from adjacent loop is positive while loop direction and assumed flow do not coincide, flow correction from adjacent loop changes its sign and opposite if (Σf)x from adjacent loop is positive while loop direction and assumed flow coincide, flow correction from adjacent loop does not change its sign.If (Σf)x from adjacent loop is negative while loop direction and assumed flow do not coincide, flow correction from adjacent loop does not change its sign and opposite if (Σf)x from adjacent loop is negative while loop direction and assumed flow do not coincide, flow correction from adjacent loop changes its sign.These four parameters are connected in Figure 3 with the same colored lines.Flow corrections (ΔQ) shown in Figure 4 with different colors are used with the related signs in Eq. ( 5).They are chosen in a similar way as explained in example from Figure 3.
So, instead simple equations as in the original Hardy Cross method, in Lobačev method, the system of Eqs.(6) has to be solved.
Underlined terms in Eqs.(6) do not exist in the Hardy Cross method.
In the Lobačev method, corrections for the first iterations are ΔQx=Δ(ΔQx)/Δ, where Δ for the first iteration is; Eq. (7) Correction for the first loop in the first iteration is; Eq. ( 9).  -0.1041 10 3.97 10 4.14 - Other corrections in the first iteration are ΔQII=-0.0644,ΔQIII=-0.0780,ΔQIV=+0.1069 and ΔQV=-0.1824.The Lobačev method is more complex compared to the original Hardy Cross method.Number of required iterations is not reduced using Lobačev procedure compared with the original Hardy Cross procedure.

Modified Hardy Cross method
The Hardy Cross method can be noted in matrix form.Gas distribution network from Figure 1 has five independent loops; Eq. (10).
Same results for the corrections are by using Eq. ( 4) for each particular loop in the network and Eq. ( 10) using matrix calculation.Epp and Fowler [10] improved original Hardy Cross method [1] by replacing some of the zeroes in non-diagonal terms of Eq. (10).For example, if pipe 8 is mutual for loop I and V, first derivative of pressure drop function for the observed pipe where flow is treated as variable will be put with the negative sign in the first column and the fifth row and also in the fifth column and the first row; Eq. (11).
In the modified Hardy Cross method, corrections for the first iterations are (12); where solutions are listed in Table 1.This procedure reduces number of iterations required for solution of the problem significantly (Figure 5).  1. Pipe diameters and lengths, as well as first assumed and final, calculated flow distribution for the network in balance are shown in Table 2.
Gas velocity in network is small (can be up to 10-15 m/s).Network can be subject of diameter optimization (as in [4]) which can be done also by using Hardy Cross method (diameter correction ΔD should be calculated for known and locked up flow where first derivative of Renouard function have to be calculated for diameter as variable).Network should stay as is if gas consumption in distance future which cannot be estimated now, are planned to be detached on nodes 5, 6, 8 and 10 (then pipes 4 and 13 will be useful in future with increased flow of gas).
Some similar examples but in case of water flow can be seen in [20].Optimization of pipe diameters in water distributive pipe network using this method can be seen in [6].1), c calculated to satisfy first Kirchhoff's law for all nodes and second Kirchhoff's law for all closed path formed by pipes (full errors in figure 1), d sign minus means that direction of flow is opposite then in initial pattern for assumed flows

Multi-point iterative Hardy Cross method
Here described multipoint method can substitute the Newton-Raphson iterative procedure used in all above described methods.It is developed used experience for solving the Colebrook equation for flow friction [11,12].Unfortunately used in Hardy Cross method, the multi-point methods do not show high performances (the solution is reached after the same number of iterations as using the Newton-Raphson procedure).
For test we used the three-point method Džunić et al. [21].Flow corrections ΔQI from Eqs. ( 10) and (11) from the first loop I should be calculated using three-point procedure; Eq. ( 13): This means; in the first iteration i, calculation of ΔQI need to follow the procedure as already explained in standard version of the method, while for the second i+1 and the third iteration i+2, the procedure is given with Eq. ( 13).For the forth iteration, the counter i restarts; i=1.This should be done for all loops in the network separately (in our case for I, II, III, IV and V).

Conclusions
Essentially what Hardy Cross did was to simplify the monumental mathematical task of calculating innumerable equations to solve complex problems in the fields of structural and hydraulic engineering long before the computer age.He revolutionized how the profession addressed complicated problems.Today, in engineering practice, the modified Hardy Cross method proposed by Epp and Fowler [10] is used rather than the original version of Hardy Cross method [1].Also, methods proposed by Hamam and Brameller [22], and those by Wood, and Charles [23], and Wood and Rayes [24] are in common practice.Also, node oriented method proposed by Shamir and Howard [25] is also a sort of Hardy Cross method.
Professional engineers use different kind of looped pipeline professional software [26], but even today, engineers invoke name of Hardy Cross with awe.When petroleum and natural gas or civil engineers have to figure out what was happening in looped piping systems [27], they inevitably turned to what is generally known as the Hardy Cross method.Original Hardy Cross is still extensively used for teaching and learning purpose [6].This method is even today constantly being improved.Here we introduced in the Hardy Cross method the multi-point iterative approach instead of the Newton-Raphson iterative approach, but it does not affect the number of required iterations to reach the final solution.
View of Hardy Cross 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.After Hardy Cross, essential idea which he wished to present involves no mathematical relations except the simplest arithmetic.
Ruptures of pipes with leakage can be detected using the Hardy Cross method because every single-point disturbances affects the general distribution of flow and pressure [28,29].
This paper has a purpose to illustrate the very beginning of modeling of gas or water pipe networks.As noted by Todini and Rossman [30], many new models have been developed since the time of Hardy Cross.

Figure 2 .
Figure 2. Rules for upper and lower sign (correction from adjacent loop; second correction)

Figure 3 .
Figure 3. Rules for terms from Lobačev equations which do not exist in Hardy Cross method

Figure 4 .
Figure 4. Rules for terms from Lobačev equations which do not exist in the Hardy Cross method applied for the network from Figure 1

Figure 5 .
Figure 5. Number of required iteration for solution using original vs. improved Hardy Cross method ) flow friction factor (dimensionless) Re Reynolds number (dimensionless) relative roughness of inner pipe surface (dimensionless) flow discharge coefficient (dimensionless) A area of ventilation opening (m 2 ) π Ludolph number; π≈3.1415

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 28 January 2019Table 1 .
Procedure for solution of flow problem for network from figure 1 using modified Hardy Cross method (first two iterations) -First iteration pipe lengths, diameters and initial flow distribution are shown in Table2and Figure1, b f calculated using a f final calculated flow in the first iteration is used for the calculation in the second iteration, etc., g if Q and Q1 have different sign, this means that flow direction is opposite than in previous iteration, etc (this will be with flow in pipe 13 between iteration 3 and 4).Preprints (www.

Table 2 .
Pipe diameters and lengths, flows and velocities of gas within pipes a network from figure1(flows are for normal pressure conditions; real pressure in network is 4x105Pa abs i.e. 3x10 5 Pa), b chosen to satisfy first Kirchhoff's law for all nodes (dash arrows in figure