Short Overview of Early Developments of the Hardy Cross Type Methods for Computation of Flow Distribution in Pipe Networks

Abstract

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 that 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 modeling 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.

1. 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]. Nonlinear electrical circuits are electrical circuits containing nonlinear 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 in theory an indefinite number of flow distributions that satisfy that the conservation of mass is possible, 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, such as for gas compressible flow. However, 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 (district heating [7] and cooling networks [8]) and ventilation systems [9,10] (a related formulation is presented in Appendix).

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 $\Delta Q$ method [13]; for $\Delta p$, see [14]), which analyzes 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] that 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,18] 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.

2. Network Piping System; Flow Distribution Calculation

2.1. 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.

2.2. A Hydraulic Model

The Renouard formula; Equation (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].

$$f=\Delta {\tilde{p}}^2=p_1^2-p_2^2=4810·\frac{{\rho }_r·L·Q^{1.82}}{D^{4.82}}$$

(1)

where f is a function of pressure, ${\rho }_r$ is relative gas density (dimensionless), ${\rho }_r=0.64$, $L$ is the pipe length (m), $D$ is the pipe diameter (m), $Q$ is flow (m³/s), and $p$ is pressure (Pa).

As shown in Appendix A, other formulas are used in the case of waterworks systems [23,24] and ventilation networks [7].

Regarding the Renouard formula (Equation (1)), one has to be careful since the pressure drop function, $f$, 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 $\sqrt{\Delta {\tilde{p}}^2}=\sqrt{p_1^2-p_2^2}$ 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 $\sqrt{\Delta {\tilde{p}}^2}$ 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 $4\times {10}^5$ Pa abs, i.e., $3\times {10}^5$ 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 (Equation (2)), where the flow is treated as a variable, is used in the Hardy Cross method.

$$f^{\prime }=\frac{\partial f\left(Q\right)}{\partial Q}=1.82·4810·\frac{{\rho }_r·L·Q^{0.82}}{D^{4.82}}$$

(2)

First assumed gas flow in each pipe is listed in the third column of

Table 1. Procedure for the solution of the flow problem for the network from Figure 1 using the modified Hardy Cross method (first two iterations)—First iteration.

Iteration 1
Loop	Pipe	a$\mathit{Q}$	b$\mathit{f}={\mathit{p}}_1^2-{\mathit{p}}_2^2$	c$\left|{\mathit{f}}^{\prime }\right|$	d$\mathbf{\Delta }{\mathit{Q}}_1$	e$\mathbf{\Delta }{\mathit{Q}}_2$	f${\mathit{Q}}_1=\mathit{Q}$
I	1	−0.3342	−144518566.8	787025109.2	−0.0994		−0.4336
	7	+0.7028	+859927106.7	2226902866.0	−0.0994		+0.6034
	8	+0.3056	+306964191.0	1828124435.8	−0.0994	−0.0532 =	+0.1530
	9	+0.2778	+800657172.4	5245486154.8	−0.0994	−0.0338 =	+0.1446
	10	−0.1364	−241342976.1	3220265516.7	−0.0994	+0.0142 ‡	−0.2217
	12	−0.0167	−6238747.4	679911398.4	−0.0994	+0.0651 ‡	−0.0511
		$\Sigma$	$f_I$ = +1575448179.8	13987715480.9
II	2	−0.0026	−80628.9	56440212.4	−0.0651		−0.0677
	11	−0.1198	−14582531.0	221537615.9	−0.0651	+0.0142 ‡	−0.1707
	12	+0.0167	+6238747.4	679911398.4	−0.0651	+0.0994 $\mp$	+0.0511
		$\Sigma$	$f_{II}$ = −8424412.4	957889226.7
III	3	−0.2338	−406110098.1	3161336093.1	−0.0142		−0.2480
	4	+0.0182	+1530938.1	153093808.5	−0.0142		+0.0040
	10	+0.1364	+241342976.1	3220265516.7	−0.0142	+0.0994 $\mp$	+0.2217
	11	+0.1198	+14582531.0	221537615.9	−0.0142	+0.0651 $\mp$	+0.1707
	14	−0.0278	−21840183.8	1429824980.5	−0.0142	−0.0338$\pm$	−0.0757
		$\Sigma$	$f_{III}$ = −170493836.7	8186058014.8
IV	5	+0.0460	+7523646.2	297674697.0	+0.0338		+0.0798
	9	−0.2778	−800657172.4	5245486154.8	+0.0338	+0.0994 ‡	−0.1446
	13	+0.0278	+21840183.8	1429824980.5	+0.0338	−0.0532 =	+0.0084
	14	+0.0278	+21840183.8	1429824980.5	+0.0338	+0.0142 $\mp$	+0.0757
		$\Sigma$	$f_{IV}$ = −749453158.7	8402810812.8
V	6	+0.0182	+3479197.2	347919720.0	+0.0532		+0.0714
	8	−0.3056	−306964191.0	1828124435.8	+0.0532	+0.0994 ‡	−0.1530
	13	−0.0278	−21840183.8	1429824980.5	+0.0532	−0.0338 $\pm$	−0.0084
		$\Sigma$	$f_V$ = −325325177.5	3605869136.3
			Iteration 2
Loop	Pipe	${\mathit{Q}}_1=\mathit{Q}$	$\mathit{f}={\mathit{p}}_1^2-{\mathit{p}}_2^2$	$\left|{\mathit{f}}^{\prime }\right|$	$\mathbf{\Delta }{\mathit{Q}}_1$	$\mathbf{\Delta }{\mathit{Q}}_2$	${\mathit{Q}}_2=\mathit{Q}$
I	1	−0.4336	−232172997.6	974431560.7	−0.0058		−0.4394
	7	+0.6034	+651439280.6	1965036192.1	−0.0058		+0.5976
	8	+0.1530	+87112249.4	1036457217.8	−0.0058	−0.0178 =	+0.1294
	9	+0.1446	+243990034.4	3070921097.1	−0.0058	−0.0098 =	+0.1290
	10	−0.2217	−584137977.5	4795666298.0	−0.0058	+0.0018 ‡	−0.2257
	12	−0.0511	−47725420.6	1700518680.1	−0.0058	−2.1·10−5 $\pm$	−0.0569
		$\Sigma$	$f_I$ = +118505168.7	13543031045.9
II	2	−0.0677	−30372941.9	816962908.0	+2.1·10−5		−0.0676
	11	−0.1707	−27780459.9	296182372.8	+2.1·10−5	+0.0018 ‡	−0.1689
	12	+0.0511	+47725420.6	1700518680.1	+2.1·10−5	+0.0058 $\mp$	+0.0569
		$\Sigma$	$f_{II}$ = −10427981.2	2813663960.8
III	3	−0.2480	−451970989.4	3317464222.8	−0.0018		−0.2497
	4	+0.0040	+99061.2	44589235.4	−0.0018		+0.0023
	10	+0.2217	+584137977.5	4795666298.0	−0.0018	+0.0058 $\mp$	+0.2257
	11	+0.1707	+27780459.9	296182372.8	−0.0018	−2.1·10−5 =	+0.1689
	14	−0.0757	−135261698.0	3251481942.9	−0.0018	−0.0098 $\pm$	−0.0873
		$\Sigma$	$f_{III}$ = +24784811.3	11705384072.0
IV	5	+0.0798	+20483898.1	467437803.0	+0.0098		+0.0896
	9	−0.1446	−243990034.4	3070921097.1	+0.0098	+0.0058 ‡	−0.1290
	13	+0.0084	+2454799.0	534076127.2	+0.0098	−0.0178 =	+0.0004
	14	+0.0757	+135261698.0	3251481942.9	+0.0098	+0.0018 $\mp$	+0.0873
		$\Sigma$	$f_{IV}$ = −85789639.2	7323916970.2
V	6	+0.0714	+41857166.9	1067095933.1	+0.0178		+0.0892
	8	−0.1530	−87112249.4	1036457217.8	+0.0178	+0.0058 ‡	−0.1294
	13	−0.0084	−2454799.0	534076127.2	+0.0178	−0.0098 $\pm$	−0.0004
		$\Sigma$	$f_V$ = −47709881.5	2637629278.1

^a Pipe lengths, diameters and initial flow distribution are shown in Table 2 and Figure 1; ^b $f$ calculated using the Renouard Equation (1); ^c $f^{\prime }$ calculated using the first derivative of the Renouard, Equation (2), where flow is variable; ^d calculated using the matrix Equation (10) and entering $\Delta Q$ with the opposite sign (using the original Hardy Cross method for Iteration 1: $\Delta Q_I$ = +0.1126; $\Delta Q_{II}$ = −0.0088; $\Delta Q_{III}$ = −0.0208; $\Delta Q_{IV}$ = −0.0892; $\Delta Q_V$ = −0.0902; using the Lobačev method for Iteration 1: $\Delta Q_I$ = −0.1041; $\Delta Q_{II}$ = −0.0644; $\Delta Q_{III}$ = −0.0780; $\Delta Q_{IV}$ = +0.1069; $\Delta Q_V$ = −0.1824); ^e$\Delta Q_2$ is $\Delta Q_1$ from the adjacent loop; ^f the final calculated flow in the first iteration is used for the calculation in the second iteration, etc.; ^g if $Q$ and $Q_1$ have a different sign, this means that the flow direction is opposite to that in the previous iteration, etc. (this occurs with the flow in pipe 13 between Iteration 3 and 4).

. 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).

3. The Hardy Cross Method; Different Versions

The Hardy Cross method is presented here: 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.

3.1. 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 $\Delta Q$ has to be computed for each loop $x$ (Equation (3)).

$$\Delta Q_x={\left(\frac{\sum \pm f}{\left|f^{\prime }\right|}\right)}_x={\left(\frac{\sum \pm 4810·\frac{{\rho }_r·L·Q^{1.82}}{D^{4.82}}}{\sum \left|1.82·4810·\frac{{\rho }_r·L·Q^{0.82}}{D^{4.82}}\right|}\right)}_x$$

(3)

For the network from Figure 1, flow corrections for the first iteration in each loop can be calculated using Equation (4).

$$\begin{array}{c}\left(\left|-f_1^{\prime }\right|+\left|f_7^{\prime }\right|+\left|f_8^{\prime }\right|+\left|f_9^{\prime }\right|+\left|-f_{10}^{\prime }\right|+\left|-f_{12}^{\prime }\right|\right)·\Delta Q_I=-f_1+f_7+f_8+f_9-f_{10}-f_{12}\\ \left(\left|-f_2^{\prime }\right|+\left|-f_{11}^{\prime }\right|+\left|f_{12}^{\prime }\right|\right)·\Delta Q_{II}=-f_2-f_{11}+f_{12}\\ \left(\left|-f_3^{\prime }\right|+\left|f_4^{\prime }\right|+\left|f_{10}^{\prime }\right|+\left|f_{11}^{\prime }\right|+\left|-f_{14}^{\prime }\right|\right)·\Delta Q_{III}=-f_3+f_4+f_{10}+f_{11}-f_{14}\\ \left(\left|f_5^{\prime }\right|+\left|-f_9^{\prime }\right|+\left|f_{13}^{\prime }\right|+\left|f_{14}^{\prime }\right|\right)·\Delta Q_{IV}=f_5-f_9+f_{13}+f_{14}\\ \left(\left|f_6^{\prime }\right|+\left|-f_8^{\prime }\right|+\left|f_{13}^{\prime }\right|\right)·\Delta Q_V=f_6-f_8+f_{13}\end{array}\}$$

(4)

In the second iteration, the calculated correction $\Delta Q$ has to be added algebraically to the assumed gas flow (the first initial flow pattern). Further, the calculated correction $\Delta Q$ 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.

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 $\Delta Q$ calculated for the adjacent loop (Figure 2).

Details of the signs of corrections were reported by 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 $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 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:

$$\Delta Q_I=\frac{1575448179.8}{13987715480.9}=+0.1126$$

$$\Delta Q_{II}=\frac{-8424412.4}{957889226.7}=-0.0088$$

$$\Delta Q_{III}=\frac{-170493836.7}{8186058014.8}=-0.0208$$

$$\Delta Q_{IV}=\frac{-749453158.7}{8402810812.8}=-0.0892$$

and

$$\Delta Q_V=\frac{-325325177.5}{3605869136.3}=-0.0902$$

3.2. 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 is 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 (Equation (5)):

$$\begin{array}{c}+\left(\left|-f_1^{\prime }\right|+\left|f_7^{\prime }\right|+\left|f_8^{\prime }\right|+\left|f_9^{\prime }\right|+\left|-f_{10}^{\prime }\right|+\left|-f_{12}^{\prime }\right|\right)·\Delta Q_I+\left|f_{12}^{\prime }\right|·\Delta Q_{II}+\left|f_{10}^{\prime }\right|·\Delta Q_{III}+\left|f_9^{\prime }\right|·\Delta Q_{IV}+\left|f_8^{\prime }\right|·\Delta Q_V=-f_1+f_7+f_8+f_9-f_{10}-f_{12}\\ \left|f_{12}^{\prime }\right|·\Delta Q_I-\left(\left|-f_2^{\prime }\right|+\left|-f_{11}^{\prime }\right|+\left|f_{12}^{\prime }\right|\right)·\Delta Q_{II}-\left|f_{11}^{\prime }\right|·\Delta Q_{III}=-f_2-f_{11}+f_{12}\\ \begin{array}{c}+\left|f_{10}^{\prime }\right|·\Delta Q_I-\left|f_{11}^{\prime }\right|·\Delta Q_{II}-\left(\left|-f_3^{\prime }\right|+\left|f_4^{\prime }\right|+\left|f_{10}^{\prime }\right|+\left|f_{11}^{\prime }\right|+\left|-f_{14}^{\prime }\right|\right)·\Delta Q_{III}-\left|f_{14}^{\prime }\right|·\Delta Q_{IV}=-f_3+f_4+f_{10}+f_{11}-f_{14}\\ +\left|f_9^{\prime }\right|·\Delta Q_I-\left|f_{14}^{\prime }\right|·\Delta Q_{III}-\left(\left|f_5^{\prime }\right|+\left|-f_9^{\prime }\right|+\left|f_{13}^{\prime }\right|+\left|f_{14}^{\prime }\right|\right)·\Delta Q_{IV}-\left|f_{13}^{\prime }\right|·\Delta Q_V=f_5-f_9+f_{13}+f_{14}\end{array}\\ +\left|f_8^{\prime }\right|·\Delta Q_I-\left|f_{13}^{\prime }\right|·\Delta Q_{IV}-\left(\left|f_6^{\prime }\right|+\left|-f_8^{\prime }\right|+\left|f_{13}^{\prime }\right|\right)·\Delta Q_V=f_6-f_8+f_{13}\end{array}\}$$

(5)

In Equation (5), signs in front of terms from the left side of the equal signs have to be determined (this is much more complex than in the Hardy Cross method). Thus, in the Lobačev method, if ${\left(\sum f\right)}_x>0$, then the sign in front of ${\left(\sum \left|f^{\prime }\right|\right)}_x$ has to be positive, and vice versa (for the first iteration, this can be seen in Table 1; $f_I=+1575448179.8>0$, $f_{II}=-8424412.4<0$, $f_{III}=-170493836.7<0$, $f_{IV}=-749453158.7<0$, $f_V=-325325177.5<0$). The sign for the other terms (these terms are sufficient in the Hardy Cross method) are determined using further rules and the scheme in Figure 3.

As shown in Figure 3 if ${\left(\sum f\right)}_x>0$ and if the assumed flow coincides with the loop direction, then the sign of flow in the adjacent pipe is negative and, if the flow does not coincide with the loop direction, then the sign of flow in the adjacent pipe is positive. Conversely, if ${\left(\sum f\right)}_x<0$ and if the assumed flow coincides with the loop direction, then the sign of flow in the adjacent pipe is positive and, if the flow does not coincide with the loop direction, then the sign of flow in the adjacent pipe is negative. This procedure determines the signs in the front of the flow corrections ($\Delta Q$), which are shown in Figure 3 with black letters (and also in Figure 4 for our pipe network example).

If ${\left(\sum f\right)}_x$ from the adjacent loop is positive, while the loop direction and assumed flow do not coincide, the flow correction from the adjacent loop changes its sign, and, conversely, if ${\left(\sum f\right)}_x$ from the adjacent loop is positive, while the loop direction and assumed flow coincide, the flow correction from the adjacent loop does not change its sign. If ${\left(\sum f\right)}_x$ from the adjacent loop is negative, while the loop direction and assumed flow do not coincide, the flow correction from the adjacent loop does not change its sign, and, conversely, if ${\left(\sum f\right)}_x$ from the adjacent loop is negative, while the loop direction and assumed flow do not coincide, the flow correction from the adjacent loop changes its sign. These four parameters are connected in Figure 3 with the same colored lines. Flow corrections ($\Delta Q$) shown in Figure 4 with different colors are used with the related signs in Equation (5). They are chosen in a similar way as explained in the example in Figure 3.

Thus, instead of the simple equations of the original Hardy Cross method, the system of equations has to be solved in the Lobačev method (Equation (6)).

$$\begin{array}{c}+13987715480.9·\Delta Q_I\underset{\_}{+679911398.4·\Delta Q_{II}+3220265516.7·\Delta Q_{III}+5245486154.8·\Delta Q_{IV}+1828124435.8·\Delta Q_V}=+1575448179.8\\ \underset{\_}{+679911398.4·\Delta Q_I}-957889226.7·\Delta Q_{II}\underset{\_}{-221537615.9·\Delta Q_{III}}=-8424412.4\\ \underset{\_}{+3220265516.7·\Delta Q_I}\underset{\_}{-221537615.9·\Delta Q_{II}}-8186058014.8·\Delta Q_{III}\underset{\_}{-1429824980.5·\Delta Q_{IV}}=-170493836.7\\ \underset{\_}{+5245486154.8·\Delta Q_I-1429824980.5·\Delta Q_{III}}-8402810812.8·\Delta Q_{IV}\underset{\_}{-1429824980.5·\Delta Q_V}=-749453158.7\\ \underset{\_}{+1828124435.8·\Delta Q_I-1429824980.5·\Delta Q_{IV}}-3605869136.3·\Delta Q_V=-325325177.5\end{array}\}$$

(6)

Underlined terms in Equation (6) do not exist in the Hardy Cross method.

In the Lobačev method, corrections for the first iterations are $\Delta Q_x=\frac{\Delta \left(\Delta Q_x\right)}{\Delta }$, where $\Delta$ for the first iteration is (Equation (7)):

$$\Delta =\left|\begin{array}{ccc}\begin{array}{c}13987715480.9\\ 679911398.4\\ \begin{array}{c}3220265516.7\\ 5245486154.8\\ 1828124435.8\end{array}\end{array}& \begin{array}{c}679911398.4\\ -957889226.7\\ \begin{array}{c}-221537615.9\\ 0\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}3220265516.7\\ -221537615.9\\ \begin{array}{c}-8186058014.8\\ -1429824980.5\\ 0\end{array}\end{array}& \begin{array}{c}5245486154.8\\ 0\\ \begin{array}{c}-1429824980.5\\ -8402810812.8\\ -1429824980.5\end{array}\end{array}& \begin{array}{c}1828124435.8\\ 0\\ \begin{array}{c}0\\ -1429824980.5\\ -3605869136.3\end{array}\end{array}\end{array}\end{array}\right|=+3.97·{10}^{+48}$$

(7)

while ΔQ_x for the first iteration is (Equation (8)).

$$\Delta \left(\Delta Q_I\right)=\left|\begin{array}{ccc}\begin{array}{c}1575448179.8\\ -8424412.4\\ \begin{array}{c}-170493836.7\\ -749453158.7\\ -325325177.5\end{array}\end{array}& \begin{array}{c}679911398.4\\ -957889226.7\\ \begin{array}{c}-221537615.9\\ 0\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}3220265516.7\\ -221537615.9\\ \begin{array}{c}-8186058014.8\\ -1429824980.5\\ 0\end{array}\end{array}& \begin{array}{c}5245486154.8\\ 0\\ \begin{array}{c}-1429824980.5\\ -8402810812.8\\ -1429824980.5\end{array}\end{array}& \begin{array}{c}1828124435.8\\ 0\\ \begin{array}{c}0\\ -1429824980.5\\ -3605869136.3\end{array}\end{array}\end{array}\end{array}\right|=-4.14·{10}^{+47}$$

(8)

The correction for the first loop in the first iteration is (Equation (9)).

$$\Delta Q_I=\frac{\Delta \left(\Delta Q_I\right)}{\Delta }=\frac{-4.14·{10}^{+47}}{+3.97·{10}^{+48}}=-0.1041$$

(9)

Other corrections in the first iteration are $\Delta Q_{II}$ = −0.0644, $\Delta Q_{III}=-0.0780$, $\Delta Q_{IV}=+0.1069$ and $\Delta Q_V=-0.1824$.

The Lobačev method is more complex compared to the original Hardy Cross method. However, the number of required iterations is not reduced using the Lobačev procedure compared with the original Hardy Cross procedure.

3.3. The Modified Hardy Cross Method

The Hardy Cross method can be noted in matrix form. The gas distribution network in Figure 1 has five independent loops (Equation (10)).

$$\left[\begin{array}{ccc}\begin{array}{c}\sum \left|f_I^{\prime }\right|\\ 0\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ \sum \left|f_{II}^{\prime }\right|\\ \begin{array}{c}0\\ 0\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}0\\ 0\\ \begin{array}{c}\sum \left|f_{III}^{\prime }\right|\\ 0\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ \sum \left|f_{IV}^{\prime }\right|\\ 0\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}0\\ 0\\ \sum \left|f_V^{\prime }\right|\end{array}\end{array}\end{array}\end{array}\right]x\left[\begin{array}{c}\Delta Q_I\\ \Delta Q_{II}\\ \begin{array}{c}\Delta Q_{III}\\ \Delta Q_{IV}\\ \Delta Q_V\end{array}\end{array}\right]=\left[\begin{array}{c}\sum \pm f_I\\ \sum \pm f_{II}\\ \begin{array}{c}\sum \pm f_{III}\\ \sum \pm f_{IV}\\ \sum \pm f_V\end{array}\end{array}\right]$$

(10)

Equation (4) provides for each particular loop in the network the same corrections as Equation (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 Equation (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 (Equation (11)).

$$\left[\begin{array}{ccccc}\sum \left|f_I^{\prime }\right|& -f_{12}^{\prime }& -f_{10}^{\prime }& -f_9^{\prime }& -f_8^{\prime }\\ -f_{12}^{\prime }& \sum \left|f_{II}^{\prime }\right|& -f_{11}^{\prime }& 0& 0\\ -f_{10}^{\prime }& -f_{11}^{\prime }& \sum \left|f_{III}^{\prime }\right|& -f_{14}^{\prime }& 0\\ -f_9^{\prime }& 0& -f_{14}^{\prime }& \sum \left|f_{IV}^{\prime }\right|& -f_{13}^{\prime }\\ -f_8^{\prime }& 0& 0& -f_{13}^{\prime }& \sum \left|f_V^{\prime }\right|\end{array}\right]x\left[\begin{array}{c}\Delta Q_I\\ \Delta Q_{II}\\ \begin{array}{c}\Delta Q_{III}\\ \Delta Q_{IV}\\ \Delta Q_V\end{array}\end{array}\right]=\left[\begin{array}{c}\sum \pm f_I\\ \sum \pm f_{II}\\ \begin{array}{c}\sum \pm f_{III}\\ \sum \pm f_{IV}\\ \sum \pm f_V\end{array}\end{array}\right]$$

(11)

In the modified Hardy Cross method, corrections for the first iterations are shown in Equation (12), and the solutions are listed in Table 1.

$$\left[\begin{array}{ccc}\begin{array}{c}+13987715480.9\\ -679911398.4\\ \begin{array}{c}-3220265516.7\\ -5245486154.8\\ -1828124435.8\end{array}\end{array}& \begin{array}{c}-679911398.4\\ +957889226.7\\ \begin{array}{c}-221537615.9\\ 0\\ 0\end{array}\end{array}& \begin{array}{ccc}\begin{array}{c}-3220265516.7\\ -221537615.9\\ \begin{array}{c}+8186058014.8\\ -1429824980.5\\ 0\end{array}\end{array}& \begin{array}{c}-5245486154.8\\ 0\\ \begin{array}{c}-1429824980.5\\ +8402810812.8\\ -1429824980.5\end{array}\end{array}& \begin{array}{c}-1828124435.8\\ 0\\ \begin{array}{c}0\\ -1429824980.5\\ +3605869136.3\end{array}\end{array}\end{array}\end{array}\right]x\left[\begin{array}{c}\Delta Q_I\\ \Delta Q_{II}\\ \begin{array}{c}\Delta Q_{III}\\ \Delta Q_{IV}\\ \Delta Q_V\end{array}\end{array}\right]=\left[\begin{array}{c}+1575448179.8\\ -8424412.4\\ \begin{array}{c}-170493836.7\\ -749453158.7\\ -325325177.5\end{array}\end{array}\right]$$

(12)

This procedure significantly reduces the number of iterations required for the solution of the problem (Figure 5).

The first two iterations for the network in Figure 1 are shown in 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. Pipe diameters and lengths, flows, and velocities of gas within pipes.

a Pipe Number	Diameter (m)	Length (m)	b Assumed Flows (m3/h)	c Calculated Flows (m3/h)	Gas Velocity (m/s)
1	0.305	1127.8	1203.2	1583.6	1.5
2	0.203	609.6	9.2	245.2	0.5
3	0.203	853.4	841.6	899.7	1.9
4	0.203	335.3	65.6	7.5	0.01
5	0.203	304.8	165.6	320.2	0.7
6	0.203	762.0	65.6	322.7	0.7
7	0.203	243.8	2530.0	2149.6	4.6
8	0.203	396.2	1100.0	462.4	1.0
9	0.152	304.8	1000.0	465.0	1.8
10	0.152	335.3	491.2	813.5	3.1
11	0.254	304.8	431.2	609.1	0.8
12	0.152	396.2	60.0	204.8	0.8
13	0.152	548.6	100.0	d−2.6	−0.009
14	0.152	548.6	100.0	312.7	1.2

^a Network from Figure 1 (flows are for normal pressure conditions; real pressure in the network is $4\times {10}^5$ Pa abs, i.e., $3\times {10}^5$ 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 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 $\Delta D$ 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].

3.4. 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 this multipoint method for acceleration of the iterative solution of the Colebrook equation for flow friction modeling [16,17]. On the contrary, for the gas network example in 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 $\Delta Q_I$ from Equations (10) and (11) from the first loop I should be calculated using the three-point procedure (Equation (13)):

$$\{\begin{array}{l}\begin{array}{l}\Delta Q_I=-\frac{{\left(\sum \pm f_I\right)}_i}{{\left(\sum \left|f_I^{\prime }\right|\right)}_i}, where i=1;\hfill \end{array}\\ \Delta Q_I=-\frac{{\left(\sum \pm f_I\right)}_i}{{\left(\sum \pm f_I\right)}_i-2·{\left(\sum \pm f_I\right)}_{i+1}}·\frac{{\left(\sum \pm f_I\right)}_{i+1}}{{\left(\sum \left|f_I^{\prime }\right|\right)}_i}, where i=2;\\ \Delta Q_I=-\frac{{\left(\sum \pm f_x\right)}_{i+2}}{{\left(\sum \left|f_I^{\prime }\right|\right)}_i·\left[1-2·\frac{{\left(\sum \pm f_I\right)}_{i+1}}{{\left(\sum \pm f_I\right)}_i}-{\left(\frac{{\left(\sum \pm f_I\right)}_{i+1}}{{\left(\sum \pm f_I\right)}_i}\right)}^2\right]·\left[1-\frac{{\left(\sum \pm f_I\right)}_{i+2}}{{\left(\sum \pm f_I\right)}_{i+1}}\right]·\left[1-2·\frac{{\left(\sum \pm f_I\right)}_{i+2}}{{\left(\sum \pm f_I\right)}_i}\right]}, where i>2;\end{array}$$

(13)

Formulas of flow corrections $\Delta Q_I$ depend on the counter i. The algorithm starts from i = 1, in which the multipoint method is the same as the original Newton–Raphson procedure (Equation (13a)):

$$\Delta Q_I=-\frac{{\left(\sum \pm f_I\right)}_i}{{\left(\sum \left|f_I^{\prime }\right|\right)}_i}$$

(13a)

In the second iteration, i = 2, flow corrections $\Delta Q_I$ have a little bit more complicated form (Equation (13b)):

$$\Delta Q_I=-\frac{{\left(\sum \pm f_I\right)}_i}{{\left(\sum \pm f_I\right)}_i-2·{\left(\sum \pm f_I\right)}_{i+1}}·\frac{{\left(\sum \pm f_I\right)}_{i+1}}{{\left(\sum \left|f_I^{\prime }\right|\right)}_i}$$

(13b)

The symbol ${\left(\sum \pm f_I\right)}_i$ represents stored values from the first iteration, whereas ${\left(\sum \pm f_I\right)}_{i+1}$ represents values from the second iteration.

For the third iteration, i = 3, flow corrections $\Delta Q_I$ have the most complicated form (Equation (13c)):

$$\Delta Q_I=-\frac{{\left(\sum \pm f_x\right)}_{i+2}}{{\left(\sum \left|f_I^{\prime }\right|\right)}_i·\left[1-2·\frac{{\left(\sum \pm f_I\right)}_{i+1}}{{\left(\sum \pm f_I\right)}_i}-{\left(\frac{{\left(\sum \pm f_I\right)}_{i+1}}{{\left(\sum \pm f_I\right)}_i}\right)}^2\right]·\left[1-\frac{{\left(\sum \pm f_I\right)}_{i+2}}{{\left(\sum \pm f_I\right)}_{i+1}}\right]·\left[1-2·\frac{{\left(\sum \pm f_I\right)}_{i+2}}{{\left(\sum \pm f_I\right)}_i}\right]}$$

(13c)

This iterative process can continue, as the formula from the third iteration is also used for iterations i = 4, 5, 6, 7, etc. This procedure should be done for all loops in the network separately (in our case for I, II, III, IV and V). However, to simplify calculations, derivative-free methods can be used [28,29].

4. Conclusions

Hardy Cross simplified mathematical modeling 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. Hardy Cross later applied these finding from structural analysis 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], as well as 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 into 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 modeling 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.