Short Overview of Early Developments of the Hardy Cross Type Methods for Computation of Flow Distribution in 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 flow resistance versus electrical resistance makes these electrical analog methods obsolete. The method by Hardy Cross is taught extensively at faculties, and it remains an important tool for the analysis of looped pipe systems. Engineers today mostly use a modified Hardy Cross method which considers the whole looped network of pipes simultaneously (use of these methods without computers is practically impossible). A method from a Russian practice published during the 1930s, which is similar to the Hardy Cross method, is described, too. Some notes from the work of Hardy Cross are also presented. Finally, an improved version of the Hardy Cross method, which significantly reduces the number of iterations, is presented and discussed. We also tested multi-point iterative methods, which can be used as a substitution for the Newton-Raphson approach used by Hardy Cross, but in this case this approach did not reduce the number of iterations. Although many new models have been developed since the time of Hardy Cross, the main purpose of this paper is to illustrate the very beginning of modelling of gas and water pipe networks and ventilation systems. As a novelty, a new multi-point iterative solver is introduced and compared with the standard Newton-Raphson iterative method.


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 13 November 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]. Non-linear electrical circuits are electrical circuits containing non-linear components. Nonlinear components can be resistive, capacitive, and inductive.
The 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 on network topology. A stable state of flow in a network must satisfy Kirchhoff's laws, which are statements of the conservation of mass and energy. Although there is an indefinite number of flow distributions that satisfy that the conservation of mass is possible in theory, only one distribution from this set also satisfies 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 outflows [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, as this relation is linear for the given quantity of flow. If, however, the increments are fairly large, this linear relation is somewhat in error, as for gas compressible flow. But if the pressure drop in pipes is minor, such as in a municipality network for natural gas distribution, the Hardy Cross method can be used without significant errors [4][5][6]. Moreover, the Hardy Cross method can also be used for water pipe networks (distring heating [7] and cooling networks [8]) and ventilation systems [9,10] (a related formulation is in Appendix A of this paper).
The Hardy Cross method is an iterative method, i.e. a method using successive corrections [4]. Lobačev and Andrijašev in the 1930s, writing in Russian, offered similar methods [11,12]. Probably because of the language barrier and the political situation in Soviet Russia, Hardy Cross was not aware of Lobačev and Andrijašev's contributions.
Today, engineers use the most improved version of the Hardy Cross method (the method [13]; for see [14]), which analyses the whole looped network of pipes simultaneously [15]. As a novel approach presented for the first time here, we tested multi-point iterative methods [16,17] which can be used as a substitution for the Newton-Raphson approach used by Hardy Cross. This approach, however, did not in this case reduce the 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 [11,12] in Section 3.2, the improved version of the Hardy Cross method [15,17] in Section 3.3, and finally the approach which uses multi-point iterative methods instead of the commonly used Newton-Raphson method in Section 3.4.

Topology of the 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, code numbers are assigned to each pipe and closed loop of pipes (represented by roman numbers for loops in Figure 1). Pipes on the network periphery are common to one loop and those in the network interior are common to two loops. Figure 1 is an example of a pipe network for distribution of natural gas for consumption in households. The next step is to write the initial gas flow distribution through pipes in the network. This distribution must be done according to Kirchhoff's first law. The choice of initial flows is not critical, and the criterion should satisfy Kirchhoff's first law for every node in the network [3]. The total gas flow arriving at a node equals the total gas flow that leaves that node. The same conservation law is also valid for the whole network in total (except for gas input and output nodes that cannot be changed during calculations; see consumption nodes in Figure 1). The sum of pseudo-pressure drops along any closed path must be approximately zero for the network to be in balance according to Kirchhoff's second law. In this paper, the flow distribution, which satisfies both of Kirchhoff's laws, will be calculated using the Hardy Cross iterative method.

A Hydraulic Model
The Renouard formula; Eq. (1) best fits a natural gas distribution system built with polyvinyl chloride (PVC) pipes [19,20]. Computed pressure drops are always less than the actual drop since the maximal consumption occurs only during extremely severe winter days [21,22]. (1) Where f is a function of pressure, ρr is relative gas density (dimensionless); here , is the pipe length (m), is the pipe diameter (m), is flow (m 3 /s), and is pressure (Pa). As shown in Appendix A of this paper, other formulas are used in the case of waterworks systems [23,24] and ventilation networks [7].
Regarding the Renouard formula; Eq. (1), one has to be careful since the pressure drop function, , does not relate pressure drop, but actually the difference of the quadratic pressure at the input and the output of the pipe. This means that is not actually pressure drop despite using the same unit of measurement, i.e. the same unit is used as for pressure (Pa). The parameter can be noted as a pseudo pressure drop. In fact the gas is actually compressed, and hence that volume of the gas is decreased, and then such a compressed volume of the gas is conveying with a constant density through the gas distribution pipeline. The operating pressure for a typical distribution gas network is Pa abs i.e. Pa gauge and accordingly the volume of the gas decreases four times compared to the volume of the gas in normal (or standard) conditions. Pressure in the Renouard formula is for normal (standard) conditions.
The first derivative f' of the Renouard relation; Eq. (2), where the flow is treated as a variable, is used in the Hardy Cross method.
(2) First assumed gas flow in each pipe is listed in the third column of Table 1. The plus or minus sign preceding flow indicates the direction of flow through the pipe for the particular loop [18,25]. A plus sign denotes counterclockwise flow in the pipe within the loop, while the minus sign, clockwise. The loop direction can be chosen to be clockwise or counterclockwise (in Figure 1 all loops are counterclockwise).

The Hardy Cross Method; Different Versions
Here will be presented the Hardy Cross method; the original approach in Section 3.1., a version of the Hardy Cross method from Russian practice in Section 3.2., the 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 the aforementioned methods.

The Hardy Cross Method; Original Approach
The pressure drop function for each pipe is listed in Table 1 (for initial flow pattern in the fourth column). The sign in front of the pressure drop function shown in the fourth column is the same as for flow from the observed iteration. The fifth column of Table 1 includes the first derivatives of the pressure drop function, where the flow is treated as a variable. The column of the function of pressure drops is computed algebraically, while the column of the first derivatives is estimated numerically for each loop. Flow correction has to be computed for each loop ; Eq. (3). (3) For the network from Figure 1, flow corrections for the first iteration in each loop can be calculated using Eq. (4).

(4)
In the second iteration, the calculated correction has to be added algebraically to the assumed gas flow (the first initial flow pattern). Further, the calculated correction has to be subtracted algebraically from the gas flow computed in the previous iteration. This means that the algebraic operation for the first correction is the opposite of its sign, i.e. add when the sign is minus, and vice versa. A pipe common to two loops receives two corrections simultaneously. The first correction is from the particular loop under consideration, while the second one is from the adjacent loop, which the observed pipe also belongs to.  Table 2 and Figure 1   The upper sign after the second correction in Table 1 is plus if the flow direction in the mutual pipe coincides with the assumed orientation of the adjacent loop, and minus if it does not ( Figure 2). The lower sign is the sign in front of correction calculated for the adjacent loop ( Figure 2). Details for signs of corrections can be seen in Brkić [18] and Corfield et al. [25]. The algebraic operation for the second correction should be the opposite of its lower sign when its upper sign is the same as the sign in front of flow , and as indicated by its lower sign, when its upper sign is opposite to the sign in front of flow .
The calculation procedure is repeated until the net algebraic sum of pressure functions around each loop is as close to zero as the desired degree of precision demands. This also means that the calculated corrections of flow and the change in calculated flow between two successive iterations is approximately zero. The pipe network is then in approximate balance and the calculation after the Hardy Cross can be terminated.
In the original Hardy Cross method, the corrections for the first iteration are: , , , , and .

A Version of the Hardy Cross Method from Russian Practice
As mentioned in the introduction, two Russian authors, Lobačev [11] and Andrijašev [12], proposed a similar method to Hardy Cross [1]. These two methods are also from the 1930s. 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 the answer to this question is no, for both sides. The main difference between the Hardy Cross and Andrijašev methods is that in the method of Andrijašev contours can be defined to include few loops. This strategy only complicates the situation, while the number of required iterations remains unchanged.
Further on the Andrijašev method can be seen from the example in the paper of Brkić [3].
Here the method of Lobačev will be shown in more detail.
In the Hardy Cross method, the influence of adjacent loops is neglected. The Lobačev method takes into consideration this influence; Eq. (5): (5) In the previous system of equations; Eq. (5), signs in front of terms from the left side of the equals sign have to be determined (this is much more complex than in the Hardy Cross method). So, in the Lobačev method if then the sign in front of has to be positive, and opposite (for the first iteration this can be seen in Table 1;  ,  ,  ,  , ). The sign for other terms (these terms are sufficient in the Hardy Cross method) will be determined using further rules and the scheme from Figure 3.    Figure 3 with the same coloured lines. Flow corrections ( ) shown in Figure 4 with different colours are used with the related signs in Eq. (5). They are chosen in a similar way as explained in the example from Figure 3.
So, instead of the simple equations of the original Hardy Cross method, the system of equations has to be solved in the Lobačev method; Eq. (6). (6) Underlined terms in Eq. (6) do not exist in the Hardy Cross method.
In the Lobačev method, corrections for the first iterations are , where for the first iteration is; Eq. (7).
The correction for the first loop in the first iteration is; Eq. (9).
Other corrections in the first iteration are =-0.0644, , and .
The Lobačev method is more complex compared to the original Hardy Cross method. But the number of required iterations is not reduced using the Lobačev procedure compared with the original Hardy Cross procedure.

The Modified Hardy Cross Method
The Hardy Cross method can be noted in matrix form. The gas distribution network from Figure 1 has five independent loops; Eq. (10). (10) Eq. (4) provides for each particular loop in the network the same corrections as Eq. (10) using matrix calculation. Epp and Fowler [15] improved the original Hardy Cross method [1] by replacing some of the zeroes in the non-diagonal terms of Eq. (10). For example, if pipe 8 is mutual for loop I and V, the first derivative of the pressure drop function for the observed pipe, where flow treated as a variable, will be put with a negative sign in the first column and the fifth row, and also in the fifth column and the first row; Eq. (11). (11) In the modified Hardy Cross method, corrections for the first iterations are (12); where solutions are listed in Table 1. (12) This procedure significantly reduces the number of iterations required for the solution of the problem ( Figure 5).  Table 1. Pipe diameters and lengths, as well as the first, assumed, and the final calculated flow distributions for the network in balance are shown in Table 2.
The gas velocity in the network is small (can be up to 10-15 m/s). The network can be the subject of diameter optimization (as in [4]), which can also be done by using the Hardy Cross method (diameter correction should be calculated for known and locked flow, where the first derivative of the Renouard function has to be calculated for diameter as a variable). The network should stay unchanged, even if planned gas consumption on nodes 5, 6, 8 and 10 increases, as pipes 4 and 13 will be useful thanks to increased gas flow.
Similar examples, but for water flow, can be seen in [26]. Optimization of pipe diameters in a water distributive pipe network using the same approach can be seen in [6].  Figure 1 (flows are for normal pressure conditions; real pressure in the network is Pa abs i.e. Pa) b chosen to satisfy Kirchhoff's first law for all nodes (dash arrows in Figure 1) c calculated to satisfy Kirchhoff's first law for all nodes and Kirchhoff's second law for all closed path formed by pipes (full errors in Figure 1) d the minus sign means that the direction of flow is opposite to the initial pattern for assumed flows

The Multi-Point Iterative Hardy Cross method
The here described multipoint method can substitute the Newton-Raphson iterative procedure used in all the above described methods. Recently, we successfully used the here presented multipoint method for acceleration of the iterative solution of the Colebrook equation for flow friction modelling [16,17]. On the contrary, for the gas network example of Figure 1, the multipoint method requires the same number of iterations as the original Newton-Raphson procedure. For the test, we used the three-point method from Džunić et al. [27]. Flow corrections from Eqs. This iterative process can continue, as the formula from the third iteration is used also for iterations i=4, 5, 6, 7 and so forth. This procedure should be done for all loops in the network separately (in our case for I, II, III, IV and V). However, in order to simplify calculations, derivative-free methods can be used [28,29].

Conclusions
Hardy Cross simplified mathematical modelling of complex problems in structural and hydraulic engineering long before the computer age. Moment distributions in indeterminate concrete structures described with differential equations were too complex for the time before computers. These finding from structural analysis, Hardy Cross later applied to balancing of flow in pipe networks. He revolutionized how the profession addressed complicated problems. Today, in engineering practice, the modified Hardy Cross method proposed by Epp and Fowler [15] is used rather than the original version of the Hardy Cross method [1]. Methods proposed by Hamam and Brameller [30], and those by Wood and Charles [31], and Wood and Rayes [32] are used in common practice [33], too. Moreover, the node oriented method proposed by Shamir and Howard [34] is also based on the Hardy Cross method.
Professional engineers use a different kind of looped pipeline in professional software [35], but even today, engineers invoke the name of Hardy Cross with awe. When petroleum and natural gas or civil engineers have to figure out what is happening in looped piping systems [36], they inevitably turn to what is generally known as the Hardy Cross method. The original Hardy Cross method is still extensively used for teaching and learning purpose [6]. 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 in our case.
The view of Hardy Cross was that engineers lived in the real world with real problems and that it was their job to come up with answers to questions in design tasks, even if initial approximations were involved. After Hardy Cross, the essential idea which he wished to present involves no mathematical relations except the simplest arithmetic.
For example, ruptures of pipes with leakage can be detected using the Hardy Cross method because every single-point disturbance affects the general distribution of flow and pressure [37,38]. This paper has the purpose of illustrating the very beginning of modelling of gas or water pipe networks. As noted by Todini and Rossman [39], many new models have been developed since the time of Hardy Cross.
Some details about the life and work of Hardy Cross are given in Appendix B.

Appendix A: Hydraulic models for water pipe networks and for ventilation systems
To relate pressure [40] [1]. His name is also famous in the field of structural engineering [53][54][55]. He developed the moment distribution method for statically indeterminate structures in 1932 [56]. This method has been superseded by more powerful procedures, but still, the moment distribution method made possible the efficient and safe design of many reinforced concrete buildings for the duration of an entire generation. Furthermore, the solution of the here discussed pipe network problems was a by-product of his explorations in structural analysis. Later, Hardy Cross was Chair of the Department of Civil Engineering at Yale, from 1937 until the early 1950s.
In 1922, Konstantin A. Čališev, emigrant from Soviet Russia, writing in Serbo-Croatian, offered a method of solving the slope deflection equations by successive approximations [57][58][59][60]. Hardy Cross improved Kališev's method as noted in [49]: "It was Hardy Cross's genius that he recognized he could bypass adjusting rotations to get to the moment balance at each and every node."

Nomenclature
The following symbols are used in this paper: relative gas density (-); here