Application of Third-Order Schemes to Improve the Convergence of the Hardy Cross Method in Pipe Network Analysis

In this study, twenty-two new mathematical schemes with third-order of convergence are gathered from the literature and applied to pipe network analysis. The presented methods were classified into one-step, two-step, and three-step schemes based on the number of hypothetical discharges utilized in solving pipe networks. The performances of these new methods and Hardy Cross method were compared by solving a sample pipe network considering four different scenarios (92 cases). The results show that the one-step methods improve the rate of convergence of the Hardy Cross method in 10 out of 24 cases (41%), while this improvement was found to be 39 out of 56 cases (69.64%) and 5 out of 8 cases (62.5%) for the two-step and three-step methods, respectively. This obviously indicates that the modified schemes, particularly the three-step methods, improve the performance of the original loop corrector method by taking lower number of iterations with the compensation of relatively more computational efforts.


Introduction
Water distribution networks (WDNs) are a critical infrastructures that affect the daily life of each and everyone who consumes purified water. In essence, management and design of WDNs require the knowledge of flow and pressure fields within WDNs. Such knowledge is exclusively provided through pipe network analysis. This analysis is mainly aimed at solving a set of first-order nonlinear differential algebraic equations with iterative methods. Hardy Cross method, which exploits the Newton-Raphson method, is one of the methods available for this purpose. However, it has several disadvantages to be employed as an adequate hydraulic solvers: (1) It requires an initial guess that conserves continuity equations in the very first iteration in advance of analyzing WDNs [1]. (2) It has a slow rate of convergence, and (3) coding of this method may be harder than other available methods [2]. These shortcomings confine the application of Hardy Cross method to solving WDNs [3].
Despite the challenges mentioned, several studies were conducted in literature about this method in a bit to improve its performance. Lopes [1] implemented the Hardy Cross method for solving pipe networks. Demir et al. [4] proposed a spreadsheet to apply a modified version of the Hardy Cross method for time-dependent simulation of WDNs. Moosavian and Jaefarzadeh [5] suggested a modified Hardy Cross method for analyzing WDNs. Baugh and Liu [6] investigated a general characterization of the Hardy Cross method. Niazkar and Afzali [7] linked MATLAB with MS Excel to solve WDNs under steady-state conditions using three dischargebased methods including the Hardy Cross method. Türkkan et al. [8] studied the implementation of the Hardy Cross method by the programming language C# for teaching calculations of WDNs. Gokyay [9] presented an MS Excel Visual Basic for Applications program for design of closed-loop pipe networks using the Hardy Cross method.
Although the Hardy Cross method is considered one of the oldest approaches for the analysis of pipe networks, it is not being obsoleted because of its advantages [6]. Not only this method is adaptable with handy calculations, but also it can be easily implemented in spreadsheets [4]. Moreover, the Hardy Cross method solves water networks without forming matrix equations [10]. Nonetheless, the convergence of this method is low in comparison to some of matrix-based approaches [11,12]. Therefore, one of the main challenges of this method is its relatively low rate of convergence.
In this study, new modified variants of the loop corrector method with third-order rate of convergence were applied for solving WDNs. For this purpose, twenty-two schemes were gathered from the literature to improve the rate of convergence of the loop corrector method. The implementations of these new modified versions are elaborated in details for a simple water network using four scenarios. The performance of the original Hardy Cross method was compared with the proposed modified versions. The results demonstrate that some of the recommended modifications improve the rate of convergence of the Hardy Cross method significantly.

Method and Materials
2.1. Analysis of Water Distribution Network. Analysis of WDN may be conducted under three different conditions [13]: (1) steady-state condition, (2) extended-period simulation, and (3) transient condition. The rate of change of the state variables of water networks with time is the main difference among these conditions. In the first condition, the temporal variation of state variables is not taken into account whereas abrupt changes of state variables are simulated in the third condition. Unlike these two extreme conditions, the time of network analysis is divided into several intervals in which the steady-state condition governs in each interval in the second condition. Therefore, any modification to the analyzing WDN under the steady-state condition can be simply extended to the extended-period simulation. In this research, the water network analysis is conducted under the steady-state condition.
The governing equations in the steady-state condition comprise (1) water continuity equation and (2) energy equation. The unknown state variables of a pipe network, i.e., pipe flow rate (Q) and nodal hydraulic head (h), are determined by solving these two equations [7]. The iteration-based methods of water distribution networks can be classified as two main approaches [10]: (1) matrix-based and (2) nonmatrix-based methods. The former ones form a matrix equation comprising the governing equations and attempt to solve these equations simultaneously whereas the latter deals with solving equations separately in two steps: First, linear continuity equations are satisfied, and second, nonlinear energy equations are solved iteratively. The former method includes the linear theory, Newton-Raphson, finite element method, and gradient global algorithm while the loop corrector method is the latter approach [14].
The analysis of water network is nothing but solving two governing equations. The first equation balances water for each node (Equation (1)) while the second equation satisfies the energy in the network for each loop (Equation (2)).
where n 1 is the number of pipes attached to a typical node whose water demand is Q d , n 2 is the number of pipes in a typical loop, K is the pipe coefficient, and n is the power. The two last parameters depend on the resistance equation adopted.
To be more specific, K and n are equal to 8f L/gπ 2 D 5 and 2, respectively, if the Darcy-Weisback (D-W) equation is utilized as the resistance equation. In this relation, L is pipe length, f is the D-W friction factor, g is the gravity acceleration, π is the pi number, and D denotes the pipe diameter. In calculations conducted in this study, the D-W equation was used, and f was calculated by the explicit Swamee-Jain's relation, i.e., f = 1:325/½ln ððε/3:71DÞ + ð5:74/Re 0:9 ÞÞ 2 , where ε/ D is relative roughness and Re is the Reynolds number [15].

Original Loop Corrector
Method. The Hardy Cross method, as a nonmatrix-based approach, initiates the analysis procedure with a set of pipe flow rates with presumed directions which satisfies the continuity equations. In this regard, the energy balance in each loop is exclusively considered to find a loop correction (ΔQ) to refine the initial guess (Equation (3)).
In Equation (3), f ðQ n Þ is the energy equation, which may give either a negative, positive, or zero value; while f ′ ðQ n Þ is the first derivative of the energy equation with positive values.
The energy function and its first derivative are evaluated using initial guesses for pipe flow rates in the first iteration. In the next iteration(s), these functions are evaluated using the corrected values of the pipe flow rates which were achieved right at the end of previous iteration. Equation (3) is applied for each loop to find the loop correction in each loop. It should be noted that the energy function, i.e., f ðQ n Þ, is a summation of energy flowing through all pipes in a typical loop.
As the loop correction is determined for each loop, the discharge values of the pipes in that loop are refined by algebraically adding this loop corrector to each pipe flow rate. If the water in a typical pipe is assumed to flow in the direction of the loop correction, the loop correction is added; while in case the pipe flow rate and loop corrector have opposite directions, the loop corrector values are subtracted from the pipe flow rate. Moreover, if a pipe is a part of two adjacent loops, the loop corrector of both loops should be utilized to refine the corresponding flow rate. This refinement procedure of discharge values in pipes will be proceeded until the loop corrector values become negligible which practically means two things: (1) The pipe flow rates will no longer refine if the procedure goes on, and (2) the numerator of Equation (3), i.e., the energy function, approaches enough close to zero which means that the energy equation is acceptably satisfied.

Advances in Mathematical Physics
Since Equation (3) is the truncated Taylor series, the Newton-Raphson, the rate of convergence of the original Hardy Cross method is second-order. This rate of convergence is relatively low in comparison to other methods and can be considered a major shortcoming of this method. In order to improve this approach, modified versions of Hardy Cross methods were proposed which has third-order of convergence in this research.

New Variants of Loop Corrector
Methods with Third-Order Convergence. In light of obtaining more precise results in relatively fewer numbers of iterations, the loop corrector method is improved by utilizing higher order of convergence schemes. In this regard, twenty-two third-order convergence Newton-Raphson schemes are selected from the literature. In order to enhance the rate of convergence of the loop corrector method, these schemes are applied as the substitution of Equation (3). In this study, these modified algorithms are classified into three groups: (1) one-step scheme, (2) two-step schemes, and (3) three-step schemes. The first group only uses one discharge, while the second and third groups utilize one and two additional discharges, so-called the hypothetical flow rates, in each loop calculation. The relations of the hypothetical flow rates used in the twenty-two algorithms are presented in the following equations: Q r = 0:5Q m + 0:5Q n , ð8Þ where Q m , Q m1 , Q m2 , Q m3 , Q r , and Q r1 are hypothetical pipe flow rates exclusively computed to calculate the loop corrector. For better clarification, the step-by-step process of solving WDNs using one of the modified schemes (Amat et al.'s algorithm) is depicted in Figure 1. As shown, this two-step scheme uses one hypothetical discharge (Q i m ), and consequently, flow velocity (V i m ), the Reynolds number (Re i m ), D-W friction factor (f i m ), and pipe coefficient (K i m ) are calculated for Q i m . In other words, these four variables are required to be computed for each hypothetical discharge. This requirement inevitably adds more computational efforts to the analysis of WDNs in comparison with the original Hardy Cross method when formulas of the modified schemes contain one or two hypothetical discharges, while this may be compensated with more accuracy achieved by the thirdorder convergence algorithms.
The modified versions of the loop corrector method are presented in the following: (1) Schemes with one discharge per pipe for each loop calculation (one-step modified schemes): (i) Chebyshev's algorithm [16]: This one-step algorithm uses the evaluations of one energy function, i.e., f ðQ n Þ, using the pipe flow rates in the nth iteration (Q n ) and one first and one second derivatives of energy equation, i.e., f ' ðQ n Þ, f " ðQ n Þ, respectively. This algorithm is shown below: (ii) Halley's algorithm [17]: This algorithm uses the evaluations of one energy function and one first and one second derivatives of energy equation: (iii) Abbasbandy's algorithm [18]: Abbasbandy's algorithm computes one energy function and one first and one second derivatives of energy equation: (iv) Chun's first algorithm [19]: This one-step algorithm computes the loop corrector using one function and one first-derivative energy equation evaluation: (v) Chun and Kim's first algorithm [20]: This algorithm comprises one function, one first and one second-derivative evaluations (vi) Chun and Kim's third algorithm [20]: Chun and Kim's third algorithm has one function, one 3 Advances in Mathematical Physics first-derivative, and one second-derivative function evaluations (2) Schemes with two-discharge per pipe for each loop calculation (two-step modified schemes): (i) Weerakoon and Fernando's algorithm [21]: This two-step algorithm utilizes the evaluations of one energy function and two first derivatives of energy equation: (ii) Frontini and Sormani's algorithm [22]: This algorithm evaluates one energy function and two first derivatives of energy equation in two steps: (iii) Amat et al.'s algorithm [23]: This two-step algorithm calculates the loop corrector by evaluating two energy functions and one first derivative of energy equation: (iv) Özban's algorithm [24]: This algorithm comprises the evaluation of one energy function and two first-derivatives of energy equation in two steps: (v) Kou  Select a set of initial guesses for discharges that satisfy continuity equations at all nodes & j = 1 P = Number of pipes in the n th loop (vi) Chun's second algorithm [19]: Chun's second algorithm utilizes two functions and two firstderivative energy equation evaluations in two steps: (vii) Chun's third algorithm [19]: This algorithm calculates loop corrector by evaluating two energy function and one first-derivative energy equation in tow steps: (x) Darvishi and Barati's algorithm [27]: This twostep algorithm has two function evaluations and one first-derivative evaluation: (xi) Zhou's algorithm [28]: This α-power mean Newton's algorithm computes the loop corrector in two steps while it uses a free parameter (α ) as the exponent of the first derivative of energy equations. In this study, the magnitude of this free parameter is assumed to be -1: where sign is the sign function (xiii) Chun and Kim's second algorithm [20]: This algorithm calculates loop corrector in tow steps as below: (xiv) Chun and Kim's fourth algorithm [20]: The calculation of loop corrector using this algorithm requires the evaluations of one function and two first derivatives: (3) Schemes with three-discharge per pipe for each loop calculation (three-step modified schemes): (i) Jisheng et al.'s third algorithm [26]: This algorithm computes the loop corrector using one function evaluation and two first-derivative energy equation evaluations: (ii) Zavalani's algorithm [30]: Zavalani's algorithm is a three-step algorithm presented in the following: 2.4. Comparison of the Modified Algorithms with Third-Order of Convergence. The twenty-two schemes presented in the previous section are compared in Table 1 in respect to the number of functions and derivatives required to be evaluated. According to this table, most of these methods are two-step algorithms, which use two discharges for each pipe in each loop calculation. Obviously, when an additional discharge for each pipe is required in a modified scheme, it demands to compute V, Re, f , and K for the new discharge values, as shown in Figure 1. This certainly requires much more calculations, which bring about a trade-off between accuracy and computational efforts.

Results and Discussion
A sample pipe network was solved using twenty-three methods, which include twenty-two modified algorithms and the original Hardy Cross method, to compare    Advances in Mathematical Physics performances of different schemes. This network, which is adopted from the literature [2,14], consists of two loops, seven pipes, and six nodes. Figure 2 depicts the layout, pipe characteristics, water demands at nodal points, and the reservoir hydraulic head connected to this pipe network.
Each and every algorithm was coded in MATLAB and MS Excel to solve the sample pipe networks for four different scenarios (92 cases overall), while applications of these programs were previously recommended for implementing numerical modeling [31][32][33] and in particular for pipe network analysis [2,7,9]. The detail of the four different scenarios is presented in the following: The detail results achieved for solving the sample pipe network using 23 methods for four scenarios are presented in the following: 3.1. Results of the One-Step Modified Schemes. The performance of the six on-step methods in solving WDNs is compared with that of the original Hardy Cross method in Figure 3. As shown, these methods solved the sample pipe network in three iterations for the first and second scenarios, whereas the third and fourth scenarios demand four, five, or six iterations of running these methods. This implies that the initial guess has an impact on the number of iterations required for solving a typical WDN using loop corrector   Advances in Mathematical Physics methods. To be more specific, the closeness of initial guess and the final solution would reduce the number of iterations in the process of analyzing WDNs. Furthermore, the performances of all schemes seem to be similar in the first two scenarios. However, based on results of the third and fourth scenarios shown in Figure 3, Chun's first algorithm was performed as the original Hardy Cross method, while the rest of one-step schemes improve the performance of the latter. Finally, Figure 3 shows that the one-step methods reach the final solutions using a lower number of iterations in 10 out of 24 cases (41%). Figure 4 compares the computation time for solving the sample pipe network using six one-step loop corrector methods and the original Hardy Cross method. As shown, the lowest computation time in each scenario belongs to the Hardy Cross method, whereas the highest computation time was achieved by Abbasbandy's algorithm for four scenarios. By considering both Figures 3 and 4, it is observed that Chebyshev's algorithm and Halley's algorithm performed better than other one-step methods based on the iteration number and computation time. Although Chun's first algorithm gained relatively low computation time in Figure 4, it utilized the same iteration numbers as the Hardy Cross method (6 and 5 iterations) for solving the third and fourth scenarios in Figure 3, respectively.

Results of the Two-
Step Modified Schemes. The numbers of iterations required by the two-step schemes and the original Hardy Cross method to solve the sample WDN are depicted in Figure 4 for the four scenarios. As shown, a few two-step algorithms improved solving the sample network by satisfying the stopping criterion by lower number of iterations. To be more specific, eight algorithms took 2 iterations for the first scenario, whereas the original Hardy Cross method and other two-step loop corrector methods required three iterations.    obtained by the Hardy Cross method from 5 to 4. According to Figure 5, Chun and Kim's second and fourth algorithms yielded to the lowest number of iterations by taking into account the four scenarios. In summary, Figure 5 demonstrates that the two-step methods improved the rate of convergence of the Hardy Cross method in 39 out of 56 cases (69.64%), while Chun's second algorithm performed worse than the latter in the second scenario.
The computation time required for solving the sample pipe network using the two-step algorithms and the Hardy Cross method is compared in Figure 6. As shown, the lowest computation time was obtained by the Hardy Cross method, whereas Chun's second algorithm (for the first three scenarios) and Ham    Additionally, these three methods used the same iteration numbers in the second scenario. Also, Figure 7 obviously indicates that both three-step methods improve the convergence rate of the original Hardy Cross method in the third and fourth scenarios. To be more precise, they solved the third scenario in four iterations, whereas the original Hardy Cross method took six iterations. Moreover, Zavalani's algorithm and Jisheng et al.'s third algorithm analyzed the fourth scenario in four iterations, while five iterations were required for the original Hardy Cross method to reach the final solutions of the fourth scenario. Based on Figure 7, the threestep methods took a lower number of iterations in 5 out of 8 cases (62.5%). Therefore, many cases considered demonstrate that the proposed loop corrector methods improve the rate of convergence of the Hardy Cross method. Figure 8 shows the computation time of the Hardy Cross method and the three-step schemes for solving the sample WDN. According to Figure 8

Conclusions
In this paper, twenty-two new loop corrector methods with third-order of convergence are proposed. Despite the original loop corrector method, i.e., Hardy Cross method, these new methods theoretically have one higher order of convergence. A sample water network was analyzed using four scenarios (92 cases overall) for comparing the performance of these new versions of loop corrector methods with the original Hardy Cross method. The results indicate that the closeness of initial guesses to the final solutions was found to be an important factor in number of iterations required to solve a typical WDS. However, considering different scenarios reveals that one-step, two-step, and three-step schemes improve the rate of convergence of Hardy Cross method by 41%, 69.64%, and 62.5%, respectively. Additionally, one of the two-step methods, Chun's third algorithm, was found to solve the sample network in more number of iterations that the Hardy Cross method for one out of four scenarios. Moreover, based on comparing the iteration number and computation time of the four scenarios of the sample pipe network, Chebyshev's algorithm and Halley's algorithm performed better than other one-step methods, while Chun and Kim's fourth algorithm and Zavalani's algorithm outperformed other two-step and three-step methods, respectively. Finally, the improvement obtained by applying the modified schemes demonstrates that the rate of convergence of the loop corrector method can be considerably increased by adopting these schemes.

Data Availability
The pipe network data used in this study are reported in the manuscript.