Next Article in Journal
Evaluation of Hybrid Constructed Wetland Performance and Reuse of Treated Wastewater in Agricultural Irrigation
Previous Article in Journal
Hydrochemical Zoning and Chemical Evolution of the Deep Upper Jurassic Thermal Groundwater Reservoir Using Water Chemical and Environmental Isotope Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Head Formulation for the Steady-State Analysis of Water Distribution Systems Using an Explicit and Exact Expression of the Colebrook–White Equation

Faculty of Civil and Environmental Engineering, Technion-Israel Institute of Technology, Haifa 32000, Israel
*
Author to whom correspondence should be addressed.
Passed away.
Water 2021, 13(9), 1163; https://doi.org/10.3390/w13091163
Submission received: 13 March 2021 / Revised: 16 April 2021 / Accepted: 19 April 2021 / Published: 22 April 2021
(This article belongs to the Section Hydraulics and Hydrodynamics)

Abstract

:
Steady-state demand-driven water distribution system (WDS) solution is the bedrock for much research conducted in the field related to WDSs. WDSs are modeled using the Darcy–Weisbach equation with the Swamee–Jain equation. However, the Swamee–Jain equation approximates the Colebrook–White equation, errors of which are within 1% for ϵ / D [ 10 6 , 10 2 ] and R e [ 5000 , 10 8 ] . A formulation is presented for the solution of WDSs using the Colebrook–White equation. The correctness and efficacy of the head formulation have been demonstrated by applying it to six WDSs with the number of pipes ranges from 454 to 157,044 and the number of nodes ranges from 443 to 150,630. The addition of a physically and fundamentally more accurate WDS solution method can improve the quality of the results achieved in both academic research and industrial application, such as contamination source identification, water hammer analysis, WDS network calibration, sensor placement, and least-cost design and operation of WDSs.

1. Introduction

Water distribution systems (WDSs) are essential infrastructures of every city and town, the purpose of which is to satisfy the water requirements of the population, of agriculture and for industry with the required quality and quantity. The hydraulic steady-state solutions of WDSs (the flows and heads) is the foundation of many, if not all, WDS academic research and industry application. Therefore, the speed and accuracy of the hydraulic simulation model that is used to find the steady-state of WDSs underpin the quality of the research outputs and industry applications. The quest for a solution method for finding the steady-state solution of a looped WDS can date back to 1936. Since then, the research community diverged into three main branches: (1) loop-based methods, (2) null-space method, and (3) range space methods.
The loop-based methods use the loop energy equations and continuity equations to model the demand-driven steady-state of WDSs. The first loop-based methods and the first WDS solution method is the Hardy–Cross method [1]. The Hardy–Cross method is an iterative manual method that uses successive approximation to solve the above nonlinear system of equations in which a set of initial flows that satisfies the mass conservation equations is successively corrected until the stopping test has been met. The loop identification and the requirement that the initial guess of flows must satisfy continuity are the factors that affect the performance of the Hardy–Cross method. The performance of the loop-based methods is dependent on the sum of the length of all identified loop. The Shortest Cycle Basis [2] is the best set of loops that minimizes the time that is required to execute any loop-based method. This is because the use of the shortest cycle basis can achieve the minimum number of non-zeros in the key matrix. Ref. [3] explored the use of minimum loops in the Newton–Raphson loop flows method and showed the time used to identify the shortest cycle basis is 100 to 10,000 times the time used to solve the network. Ref. [4] proposed two algorithms to select a set of network loops to achieve a highly sparse matrix. Although a smaller number of non-zeros in the Schur complement was reported in Creaco and Franchini [3], the substantial improvement in terms of the efficiency reported by Alvarruiz et al. [4] suggests the latter algorithm is the better practical choice. The overhead to identify the shortest cycle basis is clearly the bottleneck in the loop-based method even though the loop identification algorithm is only required to execute once for a given network topology.
The null-space methods, methods that operate in the subspace defined by a null space that is orthogonal to the column of the unknown-head node-arc incidences matrix, partitions the network into a spanning tree acyclic graph and the complementary chord tree edges. The addition of any chord tree pipe to the spanning tree graph creates a loop. The null-space method, in the context of WDS solution methods, is a special case of the loop-based methods. This is because the cycle basis created by any spanning tree graph and chord edges is a subset of the set of all cycle bases of the loop-based method, whereas the cycle basis of the loop-based method (particularly the shortest cycle basis) cannot be expressed as a combination of a spanning tree graph and the complementary chord tree edges. By sacrificing the generality, the null-space method requires fewer computation resources to find a combination of spanning tree graph and chord tree edges. Rahal [5] proposed a co-tree flows formulation in which it is necessary to (1) identify the associated circulating graph; (2) determine the demands that are to be carried by the spanning tree branches; (3) find the associated chain of branches closing a circuit for each co-tree chord; and (4) compute pseudo link head losses. Later, Elhay et al. [6] exploits the relationship between the co-tree flows and spanning tree flow by applying the Schilders’ factorization [7] to permute the A1 matrix into a lower triangular square block at the top, representing a spanning tree, and a rectangular block below, representing the corresponding co-tree.
The range space methods include a collection of methods that operate in the subspace defined by the rows of the unknown-head node-arc incidences matrix. The global gradient algorithm [8] is the most widely used WDS solution method, which solves for the flows and the heads in WDSs simultaneously by exploiting the block structure of the Jacobian matrix to reduce the size of the key matrix in the linearization of the Newton method. The graph matrix partitioning algorithm (GMPA) [9] exploits the linear relationship between the flows in the internal trees and the flows in the superlinks to speed up the solution process of the GGA.
In addition to the different solution methods developed for the simulation of the WDS steady-state, network partitioning using graph theory has been another active research avenue. Network partitioning identifies sections of a WDS network that are hydraulically independent. The first graph structure that is exploited is the forest component of a WDS. A tree in a graph is any two vertices are connected by exactly one edge. Most WDSs have trees, the collections of which are called forests. Simpson et al. [10] proposed a forest-core partitioning algorithm to partition a WDS graph into its linear forest component and nonlinear core component. The flows in the forest pipes can be computed a priori and the heads in the forest nodes can be computed a posteriori by a linear process. Qiu et al. [11] proposed a bridge-block partitioning algorithm to partition a WDS graph into several linear independent blocks and bridges. The steady-state solution of one block of a WDS can be computed independently from other WDS blocks.
EPANET2 [12] is one of the most widely used WDS simulation packages, in which the GGA is used to provide steady-state demand-driven solution of WDSs. The code for EPANET 2 is in the public domain, allowing many studies to be conducted. These include the least-cost design and operation of WDSs, sensor placement, chlorine decay models calibration, contamination event detection, network vulnerability analysis, cyber-attack detection, network decontamination, and many others. Any inaccuracies that existed in the WDS solution method used will be inherited and sometimes exacerbated. One of the accuracy-related problems is the use of the Swamee–Jain equation [13], an approximation of the implicit Colebrook–White equation [14], a general equation that can be used when the Reynolds number is greater than or equal to 4000, to calculate the friction factor. It is reported in Swamee and Jain [13] that the errors involved in friction factor are within 1% for ϵ / D [ 10 6 , 10 2 ] and when R e [ 5000 , 10 8 ] . All WDS solution methods that are currently available (1) ignored the above 1% error and (2) extended the applicability of the Swamee–Jain equation to calculate the friction factor when R e 4000 and all values of ϵ / D . This is mainly due to the solutions of the inexplicit Colebrook–White equation are required many times for multiple iterations, which is time-consuming.
In this paper, an iterative solution method that finds the steady-state demand-driven WDS solution using the Colebrook–White equation for the turbulent flow regime is proposed. This is achieved by (1) expressing the Colebrook–White equation, an implicit function, as an explicit, exact, and differentiable function to describe the friction factor in the turbulent flow regime, (2) using the Hagen-Poiseuille equation to describe the friction factor in the laminar flow regime, and (3) using cubic interpolation to fit a curve between the turbulent and the laminar flow regimes to describe the transitional flow regime. It is important to note that unlike the Swamee–Jain equation, the explicit expression of the Colebrook–White equation agree at all points, the use of which eliminates the inaccuracy that is associated with the Swamee–Jain equation. It is shown in this paper that the different steady-state solutions of WDSs are observed which can be a critical problem for some research areas (such as the least-cost design and operation, contamination event detection, network vulnerability analysis, sensor placement, district meter area, etc.).
This paper is organized as follows. A description of the general WDS demand-driven steady-state problem is given in the next section. This is followed by Section 3 the derivation of the proposed head formulation. Section 4 gives an algorithmic description of the proposed head formulation, followed by the validation of the proposed friction factor equation against the Colebrook–White equation. Six case study networks are then described in Section 6, the results of which are discussed in the next section. The last section offers some conclusions.

2. General WDS Demand-Driven Steady-State Problem

2.1. Definitions and Notation

Consider a water distribution system that contains n p pipes, n j junctions, and n r fixed-head nodes.
The i -th node of the network has two properties: its nodal demand d i and its elevation head z i . Let h = h 1 , h 2 , , h n j T denote the vector of unknown heads, d = d 1 , d 2 , , d n j T denote the vector of nodal demands, e l = e l 1 , e l 2 , e l n r T denote the vector of fixed-head elevations.
The p-th pipe of the network can be characterized by its diameter D p , length L p , flows q p , and Hazen–William coefficient C p for Hazen–William head loss model or roughness height ϵ p for Darcy–Weisbach head loss model. Let q = q 1 , q 2 , , q n p T denote the vector of unknown flows, C = C 1 , C 2 , , C n p T denote the vector of Hazen–William coefficients, ϵ = ϵ 1 , ϵ 2 , , ϵ n p T denote the vector of roughness heights, L = L 1 , L 2 , , L n p T denote the vector of pipe lengths, D = D 1 , D 2 , , D n p T denote the vector of pipe diameters,
The matrix A 1 is the full-rank, unknown-head, node-arc incidence matrix. The matrix A 2 is the fixed-head node-arc incidence matrix. The head loss exponent n is assumed to be dependent only on the head loss model: n = 2 for the Darcy–Weisbach head loss model and n = 1.852 for the Hazen-Williams head loss model. The head loss within the pipe p is modeled by h f p = r p q p | q p | n 1 . Denote by G q R n p × n p , a diagonal square matrix with elements G p p = r p | q p | n 1 for j = 1 , 2 , . . . . n p . Denote by F q R n p × n p , a diagonal square matrix where the p-th element on its diagonal F p p = d d q p G p p q p .

2.2. System of Equations

The steady-state flows and heads in a WDS system are modeled by the demand-driven model (DDM) continuity Equation (1) and the energy conservation Equation (2):
A 1 T q d = O
h f A 1 h A 2 e l = O ,
which can be expressed as
G ( q ) A 1 A 1 T O q h A 2 e l d = 0 ,
where its Jacobian matrix used in the solution process is
J = F ( q ) A 1 A 2 T O
and it is sometimes referred to as a nonlinear saddle point problem [15].
This nonlinear system is often solved by the Newton method, in which q ( m + 1 ) and h ( m + 1 ) are repeatedly computed from q ( m ) and h ( m ) by
F ( m ) ( q ( m ) ) A 1 A 1 T O q ( m + 1 ) q ( m ) h ( m + 1 ) h ( m ) = G ( m ) q ( m ) A 1 h ( m ) A 2 e l A 1 T q ( m ) d ,
until the relative differences | | q ( m + 1 ) q ( m ) | | | | q ( m + 1 ) | | and | | h ( m + 1 ) h ( m ) | | | | h ( m + 1 ) | | are sufficiently small.

3. Derivation of the Head Formulation

Consider the continuity equations in Equation (1) and a vector of unknown heads, h, for the n j nodes in a network.
Assume the flows in pipe p can be expressed as a function of the head loss in pipe p:
q p = f ( [ h f ] p ) .
Let a * , j R n p × 1 denote the j-th column of the A 1 matrix, the continuity equation can be rewritten as:
f c = f c 1 f c 2 f c j f c n j = a * , 1 T q + d 1 a * , 2 T q + d 2 a * , j T q + d j a * , n j T q + d n j = a * , 1 T f ( h f ) + d 1 a * , 2 T f ( h f ) + d 2 a * , j T f ( h f ) + d j a * , n j T f ( h f ) + d n j
is the mass balance equation expressed for every node in a WDS.
The partial derivative of f c ( q 1 , q 2 , , q n p ) with respect to h f can be expressed as:
J = f c h f = f c 1 h 1 f c 1 h 2 f c 1 h n j f c 2 h 1 f c 2 h 2 f c 2 h n j f c n j h 1 f c n j h 2 f c n j h n j
in which the partial derivative of f c j ( q 1 , q 2 , , q n p ) with respect to h f p can be expressed as:
f c j h j = i = 1 n p f c j q i q i h j .
As a result, Equation (7) can be written as:
J = f c 1 q 1 f c 1 q 2 f c 1 q n p f c 2 q 1 f c 2 q 2 f c 2 q n p f c n j q 1 f c n j q 2 f c n j q n p q 1 h 1 q 1 h 2 q 1 h n j q 2 h 1 q 2 h 2 q 2 h n j q n p h 1 q n p h 2 q n p h n j + d 1 h 1 0 0 0 d 2 h 2 0 0 0 d n j h n ) j
where the first matrix of the matrix multiplication above can be written as A 1 T q / q , which is A 1 T , and the partial derivative of q p ( h f p ) (Equation (6)) with respect to h j can be expressed as:
q p h j = q p h f p h f p h j .
Substitute Equation (10) into Equation (9), we get:
J = A 1 T q 1 h f 1 h f 1 h 1 q 1 h f 1 h f 1 h 2 q 1 h f 1 h f 1 h n j q 2 h f 2 h f 2 h 1 q 2 h f 2 h f 2 h 2 q 2 h f 2 h f 2 h n j q n p h f n p h f n p h 1 q n p h f n p h f n p h 2 q n p h f n p h f n p h n j + d h = A 1 T q 1 h f 1 0 0 0 q 2 h f 2 0 0 0 q n p h f n p h f 1 h 1 h f 1 h 2 h f 1 h n j h f 2 h 1 h f 2 h 2 h f 2 h n j h f n p h 1 h f n p h 2 h f n p h n j + d h
where the third matrix in the matrix multiplication above can be expressed as A 1 h + A 2 e l / h , which is A 1 . Denote by W R n p × n p , a diagonal square matrix where the p-th element on its diagonal [ W ] p p = f ( h f p ) . The Jacobian matrix in Equation (11) can be expressed as:
J = A 1 T W A 1 + d h
Omitting term d h for the demand-driven analysis, this new WDS solution method can be solved by the Newton method in which d h ( m + 1 ) are repeatedly computed from h ( m ) using:
J h ( m + 1 ) h ( m ) = A 1 T q d
Moving h ( m ) to the right-hand-side of Equation (13), we get:
J h ( m + 1 ) = A 1 T q ( m ) d + A 1 T W h f ( m ) A 2 e l
Equation (6) is now derived for Hazen-Williams and Darcy–Weisbach head loss models in the next subsections.

3.1. Hazen-Williams Equation

The Hazen-Williams Equation head loss equation is
h f p = 10.67 L p C p 1.852 D p 4.8704 q p | q p | 0.852
and can be rewritten as:
q p = f H W ( h f [ ) = C p D p 2.63 3.6 L p 0.54 h f p | h f p | 0.46
and the partial derivative of q p ( h f p ) (Equation (16)) with respect to h f p can be expressed as:
f H W ( h f p ) = 0.54 C p D p 2.63 3.6 L p 0.54 | h f p | 0.46

3.2. Darcy–Weisbach Equation

The Darcy–Weisbach head loss equation is:
h f p = 8 f p L p π 2 g D p 5 | q p | q p
Equation (18) can be rearranged into
f p = π 2 g D p 5 8 L p | q p | q p h f p
The friction factor for turbulent flow regime ( R e 4000 ) can be calculated by using the Colebrook–White equation:
1 f p = 2 log ϵ p 3.7 D p + 2.51 R e p f p
in which the Reynolds number is
R e p = 4 | q p | π D p ν .
Substitute Equations (19) and (21) into the RHS of Equation (20), we get:
f p = 0.25 log 2 ϵ p 3.7 D p + 2.51 2 ν L p 0.5 2 | h f p | 0.5 g 0.5 D p 1.5
Denote by S p = h f p / L p , the hydraulic gradient in pipe p. The friction factor in pipe p can be expressed as:
f p = 0.25 log 2 ϵ 3.7 D p + 2.51 2 ν 2 g 0.5 D p 1.5 | S p | 0.5
The friction factor for laminar flow regime ( R e 2000 ) can be calculated by using the Hagen-Poiseuille equation:
| q p | = 16 π D p ν f p
Substitute Equation (24) into Equation (18), we get:
f p = 2048 υ 2 g D p 3 | S p | 1
The friction factor for transitional flow regime ( 2000 R e 4000 ) can be calculated by using a cubic interpolation between the laminar flow regime and the turbulent flow regime using:
f = ( c 0 + S ( c 1 + S ( c 2 + S c 3 ) ) )
in which S = ( | S p | S p 2000 ) / ( S p 4000 S p 2000 ) , c 0 = 0.032 , c 1 = 0.032 , c 2 = 0.032 + 3 f p 4000 f p 4000 , c 3 = 0.032 2 f p 4000 + f p 4000 , f p 4000 is the friction factor for pipe p when the Reynolds number is 4000, and f p 4000 is the derivative of friction factor when the Reynolds number is 4000.
Let a p = 2.51 2 υ 2 g 0.5 D p 1.5 , we can express:
f p = 2048 υ 2 g D p 3 | S p | 1 | S p | S p 2000 ( c 0 + S ( c 1 + S ( c 2 + S c 3 ) ) ) S p 2000 | S p | S p 4000 0.25 ln 2 ( 10 ) ln 2 ϵ 3.7 D + a p | S p | 0.5 | S p | S p 4000
Let b p = π g 0.5 D p 2.5 2 2 , the flows in pipe p as a function of the head loss in pipe p (Equation (6)) for the Darcy–Weisbach head loss model can be expressed as
q p ( S p ) = b p S p | S p | 0.5 f p 0.5
where its derivative can be expressed as:
q ( S p ) = f D W ( S p ) = 0.5 b p ( | S p | f p ) 0.5 1 S p f p d f p d S p
in which
d f p d S p = 2048 υ 2 g D P 3 ( S P | S P | ) 1 | S P | S p 2000 c 1 + 2 S c 2 + 3 S 2 c 3 ( S 4000 S 2000 ) S p | S p | 1 S p 2000 | S p | S p 4000 a p ln 2 ( 10 ) 4 ln 3 ( ϵ 3.7 D + a p | S p | 0.5 ) ( ϵ 3.7 D | S p | 0.5 + a p ) S p | S p | S p 4000 ,
and d S p d h f p = L 1 . It is important to note that special care must be taken to calculate the q ( h f ) value for laminar flows as the pipe head loss will be cancelled when calculating the pipe flow:
q ( S p ) = b p 2048 υ 2 g D p 3 0.5 S p
as a result, the derivative of the pipe flow with respect to the pipe head loss in the laminar flow regime can be expressed:
q ( S P ) = b p 2048 υ 2 g D p 3 0.5
which is a constant value. Therefore, a pipe with zero head loss will not cause any numerical problems as a result of the zero head loss.
Finally, S p 2000 and S p 4000 need to be identified for pipe p for all pipes. S p 2000 can be easily identified as f 2000 = 64 / 2000 = 0.032 for any pipes and q p 2000 = 500 π D p ν . Therefore,
S p 2000 = 64000 ν 2 g D p 3 .
The determination of S p 4000 involves the solution of f p 4000 using the implicit Colebrook–White equation for pipe p when R e = 4000 and substitute the value of f p 4000 in Equation (2). This is only time when the direct solution of the implicit Colebrook–White Equation is required.

4. Head Formulation Algorithm

The steps of the proposed head formulation are described in Algorithm 1. The proposed head formulation, a single-phase formulation, within the iterative phase (between lines 4 and 18 in Algorithm 1 is similar in terms of the computational intensive when compared to the global gradient algorithm, a two-phase formulation.
Meanwhile, the overhead of the proposed head formulation, particularly the computation of the S p 4000 that requires the solution of the Colebrook–White equation for the Darcy–Weisbach head loss model, can increase the computation burden of the algorithm. However, as the manufacturing limitation, there is only a limit number of pipes that is available for the construction of WDSs. Therefore, the number of combinations of the pipe roughness heights and pipe diameters is limited, so does the number of distinct values of S p 4000 . This is also true for the least-cost design problem of a WDS. S p 4000 value is only required to be computed and stored for each of the commercially available pipes a priori. This stored value can then be retrieved during the optimization phase.
Moreover, Equation (16) is just an inverse function of Equation (15). As a result, the d q d h f required in the proposed head formulation is the inverse of the d h f d q required in the GGA. Therefore, the Schur Complement in the GGA, A 1 T F 1 A 1 , is the same as the Jacobian matrix shown in Equation (12) as W = F 1 . Thus, if the proposed head formulation and the GGA starts with the same initial starting points when applied to a WDS with Hazen-Williams head loss model, they will also have the same iterative solutions and the same final solution.
Algorithm 1: Head formulation
Water 13 01163 i001

5. Validation of the Proposed Friction Factor Equation

The proposed friction factor equation is validated in this section. The proposed friction factor equation and the Swamee–Jain friction factor equation will be used to compute friction factors for pipes with the value of ϵ = (1 mm, 0.1 mm, 0.01 mm, 0.001 mm), D = 100 mm and Reynolds number ranging from 4000 to 10 8 , and the computed friction factor value in turbulent flow regime will be plotted against the value of the friction factor computed by using the Colebrook–White equation. As can be seen from Figure 1, the differences are observed between the friction factors computed by using the Swamee–Jain friction factor equation and that computed by using the Colebrook–White equation as shown in Swamee and Jain [13], whereas the friction factors computed by using the proposed friction factor equation are the same as that computed by using the Colebrook–White equation.

6. Case Studies

The proposed head formulation was implemented in WDSLib [16] and has been applied to six WDSs. The networks used here were assigned network identifiers from N 1 to N 6 (see details in Table 1):
  • Network identifier N 1 is assigned to Balerma Network, first introduced by Reca and Martínez [17], is comprised of 454 pipes, 443 junctions, and four reservoirs.
  • Network identifier N 2 is assigned to Richmond Network, first introduced by Van Zyl et al. [18], is comprised of 834 pipes, 848 junctions, and eight reservoirs.
  • Network identifier N 3 is assigned to exnet Network, first introduced by Farmani et al. [19], is comprised of 2465 pipes, 1890 junctions, and three reservoirs. This is important to note that valves are replaced by pipes.
  • Network identifier N 4 is assigned to the large Network used in Sitzenfrei et al. [20]. This network has 4021 pipes, 3557 junctions and one reservoir.
  • Network identifier N 5 is assigned to Network 2 of the Battle of Network Sensors competition [21]. Network N 5 is comprised of 14,830 pipes, 12,523 junctions, and seven reservoirs. This is important to note that valves and pumps are replaced by pipes and demand patterns are removed.
  • Network identifier N 6 is assigned to virtRome Network, first introduced by [20], is comprised of 157,044 pipes, 150,630 junctions, and four reservoirs.
All case studies were performed on an Intel Core i7-7700 running at 3.6 GHz with four cores in C++ under IEEE-standard double-precision floating arithmetic [22] with machine epsilon ϵ m a c h = 2.204 × 10 16 .
The infinity norm of the relative head differences h m + 1 h m h m + 1 will be used as convergence test. Without a WDS solution method that can be used as a benchmark, the friction factor residual in Equation (34), the continuity equation residual in Equation (35), and the energy equation residual in Equation (36) were used to validate the correctness of the proposed head formulation.
σ f ( m ) = 1 f ( m ) + 2 log ϵ 3.7 D + 2.51 R e ( m ) f ( m )
σ c ( m ) = A 1 T q ( m ) + d ( m )
σ E ( m ) = h f ( m ) + A 1 h ( m ) + A 2 e l ( m )

7. Results and Discussion

Figure 2 shows the convergence of the proposed iterative head formulation. Networks N 1 and N 2 are both relatively small network, both of which have a small number of loops (11 loops for network N 1 and 14 loops for network N 2 ). The stopping test and the continuity residual have been met for both smaller size networks. However, convergence problem has been observed for the medium size networks (networks N 3 and N 4 ) and large size networks (networks N 5 and N 6 ). It is worth noting that the relatively head differences use in the stopping test has a better convergence property that the continuity residual. This is particularly pronounced in network N 2 as shown in Figure 2b. The stopping test has been met at seven iteration, whereas the continuity residual is 10 3 . This significant difference is also observed in Figure 2c–f for networks N 3 N 6 .
Using junction 776 and reservoir E in network N 3 as an example, the head iterates of junction 776 is observed to be oscillating around the elevation head of reservoir E as shown in Figure 3a. The relative head difference of the nodal head at junction 775 is 3.94 × 10 10 at 6th iteration. However, this small perturbation in the nodal head caused the flow direction reversed. This flow direction reversal happened three more times while the value of the nodal heads converges to the true value. This is because the pipe (Pipe 1865) that connects junction 776 and reservoir E is 1 m in length and 0.999 m in diameter with 0.15 mm roughness height, which means the resistance factor of pipe 1865 is a very small value. This can also be seen from Figure 3a as the head loss between reservoir E and junction 776 is 10 8 .
Once a damping factor of 0.67 has been applied, the head convergence of this particular node is more well-behaved as a faster convergence is achieved and no head oscillation has been observed as can be seen from Figure 3b. Please note that in Figure 3b and in Figure 4 the damped Newton’s method was used. Within this scheme, the derivative within the Newton method is multiplied by a damping factor between zero and one for accelerating convergence. On the one hand, the application of the damping factor caused a slower head convergence of the node that is well-behaved before its application. This can be seen from Figure 4a,b as a slower rate of convergence is observed when compared to that from Figure 2a,b. On the other hand, the application of the damping factor guaranteed convergence for networks N 3 N 6 as shown in Figure 4c–f when compared to before the application of the damping factor as shown in Figure 2c–f.
The convergence properties, including the relative head differences, continuity residuals, energy residuals, and Colebrook–White equation residuals, at the final iteration of the proposed head formulation after applying the damping factor is shown in Table 2.
Table 3 presents the detailed timing results of pre-processing, iterative phase and post-processing operations of the GGA and proposed head formulation applied to the six case study networks. The total wall-clock time required to apply the proposed head formulation is higher than that required by the GGA for all six case study networks. This is because the proposed head formulation requires more time to perform the iterative phase while the GGA and the proposed head formulation require a similar time to execute the pre-processing and post-processing operations. In addition, a similar amount of per-iteration runtime is required by the GGA and the proposed head formulation, which means the longer wall-clock time required by the proposed head formulation is the result of higher number of iterations required due the damping factor applied.
Table 4 shows the detailed statistics of the absolute head differences between the GGA and the proposed head formulation applied to each of the six case study networks. This head differences in the results of optimization problems can cause constraints violation of the optimal solution identified. Take the Balerma network N 1 , one of the well explored network, as an example. This problem is first introduced by [17]. Two least-cost designs presented by [23,24], both of which using the GGA formulation to find the steady-state solution in which the Swamee–Jain equation is used to model the pipe friction factors, have found to have nodal pressure violations when using the proposed head formulation as shown in Table 5, in which the Colebrook–White equation is used to model the pipe friction factors. The above two least-cost designs of network N 3 are both infeasible as the Swamee–Jain equation approximates the Colebrook–White equation. This is because (1) the errors involved in Swamee–Jain friction factor are within 1% for ϵ / D [ 10 6 , 10 2 ] and R e [ 5000 , 10 8 ] as reported in [13], whereas GGA has extended the applicability of the Swamee–Jain equation to all values of ϵ / D and R e 4000 and (2) the errors involved in Swamee–Jain friction factor have been ignored in all optimization problems.
In addition to the differences in the head solutions, the use of friction factor produced by the Colebrook–White equation in the proposed head formulation also produces different flow results when compared to that produced by the GGA as shown in Table 6 and the spatial distribution of the pipes with different flows is shown in Figure 5. The flows of the pipes in the looped component found by using the proposed head formulation for each of the six case study networks are different from that found by using the GGA, whereas the flows of the pipes in the forest component found by using the head formulation for each of the six case study networks are the same as that found by using the GGA. In addition, the number of pipes with flow direction reversal is ranging from six pipes in network N 2 to 730 pipes in network N 6 .
In addition to the single-period steady-state WDS simulation where the boundary conditions (pumps and tanks) are fixed, errors caused by the difference between the Colebrook–White equation and the Swamee–Jain equation can accumulate as the boundary condition in the extended-period simulation. The error accumulation is manifested as the different tank level at each time step, the time when an operation activated by the trigger level, and the pump operating points. Net3 in EPANET is used as an example to demonstrate the error accumulation described above. It is important to note that the Net3, which is a simple WDS with 92 nodes, two reservoirs, three tanks, two pumps, and 117 pipes, in EPANET has been converted from the Hazen-William head loss model to the Darcy–Weisbach head loss model. The model difference between the GGA and the head formulation is relatively small between 00:00 to 05:00, which is also observed in network N 1 . As can be seen from Figure 6, however, the tank levels start to diverge as the tank level rises. This is because the insignificant model errors for small network can accumulate over a period of time in the form of boundary conditions. It is also worth noting that when the tank operation mode changes, from filling to discharging or from discharging to filling, the tank level differences between the two model start to narrowing, because the cancelling of errors, but never reaching zero. The maximum tank level difference is 0.074 m for tank 1, 0.12 m for tank 2 and 0.37 m for tank 3. Due to the error build-up in the boundary conditions, the differences between GGA and the head formulation flow results are relatively small with no flow direction reversal before the mid-day as can be seen from Figure 7a. After 12:00, however, significant differences between GGA and the head formulation flow results are observed with the reversal of flow direction as can be seen from Figure 7c–e. Finally, flow result differences start to decrease as the tank level differences narrow. The GGA and the head formulation start with the same tank levels and end up with different tank levels over the 24 h simulation period. This boundary condition differences will keep accumulating over a longer simulation period.

8. Possible Applications

The main improvement of using the proposed head formulation is the use of the Colebrook–White equation addresses a fundamental error that is associated with the use of Swamee–Jain equation, an approximation of the Colebrook–White equation, in all existing WDS solution methods. It is shown in the last section that both flow and head results of the head formulation are different from that of the global gradient algorithm. The differences in flow and head results can significantly affect several research areas.

8.1. Water Quality Simulation

Water quality is one of the most important water distribution system research subject. Water quality is simulated using the chemical transportation equation (advection) and the chemical decay equation, both of which are functions of pipe flows and time for both the contaminant and the disinfectant. As a result, the differences in flow results between the proposed head formulation and the GGA, particularly the flow reversal, can significantly affect the evaluation of the chemical transportation equation. This can be a critical problem for the identification of the contamination sources, network decontamination, disinfection by-product modeling, water quality sensor placement, and others.
In addition, the water quality modeling is temporal in nature, therefore, the flow results of the extended-period simulation are normally used in the chemical transportation equation. As can be seen from Figure 7 that the differences in flow results have no correlation and the pipes with flow reversal are not always the same. As a result, there is another level of error accumulation in water quality simulation on top of the extended-period simulation as opposed to the water quality simulation using the steady-state simulation where only the initial conditions are different. These two levels of error accumulation are compounded so that the errors in the water quality modeling can be amplified.

8.2. Water Hammer Analysis

Water hammer analysis is another research area that is affected by the flows, heads, and head losses in pipes and nodes of water distribution systems. The governing equations for the water hammer analysis are the unsteady momentum equation and the unsteady continuity equation. Both unsteady equations are functions of heads, flows, and the direction of flows. The reversal of the flow direction in a pipe will significantly alter the result of water hammer analysis in a looped water distribution system.

8.3. WDS Network Calibration

The proposed head formulation can also improve the quality of the WDS calibration models. Water distribution system calibration is a process of comparing model results with field data and making the appropriate adjustments so that both results agree, and it is usually applied to the estimation of pipe roughness values and nodal demands [25]. All existing calibration studies compared the ’measurement data’, which is generated using EPANET instead of using the real measurement data, against the existing model results, both of which uses the Swamee–Jain equation to model the turbulent flow regime. As a result, the errors existed in the Swamee–Jain approximation is inherited by the calibrated water distribution system. A better network calibration model can be produced by using the proposed head formulation as the differences between the Colebrook–White equation and the experimental data are smaller than that between the Swamee–Jain equation and the experimental data.

8.4. WDS Optimization and Operation

The small head differences can invalidate all existing optimal solutions of water distribution system modeled by Darcy–Weisbach head loss equation. Examples of the pressure violations in the optimal solutions found in [23,24] for the Balerma network using the global gradient algorithm using Swamee–Jain equation are shown in Table 5, both of which become infeasible if the proposed Colebrook–White equation is used.

9. Conclusions

This paper presents an efficient iterative head formulation for the steady-state demand-driven solution of water distribution systems for both Darcy–Weisbach and Hazen-Williams head loss model. When the Hazen-Williams head loss model is used, the proposed solution method produces the same final and iterative flow and head solutions if the same initial guess is used as the global gradient algorithm. When the Darcy–Weisbach head loss model is used, an exact and explicit expression of the inexplicit Colebrook–White friction factor equation is proposed in this study. A cubic interpolation between this explicit expression of the inexplicit Colebrook-White equation for the turbulent flow regime and the Hagen-Poiseuille equation for the laminar flow regime is generated to describe the friction factor in the transitional flow regime.
The main features of the proposed head formulation of the steady-state demand-driven WDS simulation include:
(1)
friction factor for the turbulent flow regime can be calculated using an explicit and exact expression of the Colebrook–White equation without the need for an iterative method;
(2)
the use of the proposed head formulation can significantly improve the accuracy of the steady-state demand-driven solution of WDSs when compared to the GGA. This is because the proposed head formulation uses an explicit and exact expression of the Colebrook–White equation to calculate friction factor for the turbulent flow regime as opposed to the Swamee–Jain equation, an approximation of the Colebrook–White equation, used in previous WDS solution methods.
The efficacy of the proposed head formulation has been demonstrated by applying it to six case study networks, the results of which have been validated using the continuity residuals, energy residuals, and Colebrook–White residuals. It should also be observed that the proposed method could be selected for analytical (e.g., Di Nucci [26]) and for numerical applications (e.g., Pasculli [27]), and in particular for improving the Wall Function.Differences between the proposed head formulation and the GGA have been observed in both the flows and heads. On the one hand, the flow differences, particularly the flow direction reversal, can be a critical problem for some research areas such as water quality simulation, contamination source identification, water hammer analysis, WDS network calibration, sensor placement, and network clustering. On the other hand, the heads differences can cause pressure violations in the WDS least-cost design when the GGA is used to perform the steady-state demand-driven analysis. It is important to note that damping factor has been applied to the Newton method to achieve a more stable convergence. The choice of initial guess, another important factor in the Newton method convergence, and inclusion of valves, pumps, and other network elements will be interesting avenues for future research.

Author Contributions

Conceptualization, M.Q.; Data curation, M.Q.; Formal analysis, M.Q.; Funding acquisition, A.O.; Methodology, M.Q.; Project administration, A.O.; Supervision, A.O.; Validation, M.Q.; Visualization, M.Q.; Writing—original draft, M.Q.; Writing—review and editing, M.Q. and A.O. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the ISRAEL SCIENCE FOUNDATION (grant No. 555/18).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

This research was supported by the ISRAEL SCIENCE FOUNDATION (grant No. 555/18).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Cross, H. Analysis of Flow in Networks of Conduits or Conductors; Engineering Experiment Station Bulletin, University of Illinois: Champaign, IL, USA, 1936. [Google Scholar]
  2. Horton, J.D. A polynomial-time algorithm to find the shortest cycle basis of a graph. SIAM J. Comput. 1987, 16, 358–366. [Google Scholar] [CrossRef]
  3. Creaco, E.; Franchini, M. Comparison of Newton-Raphson global and loop algorithms for water distribution network resolution. J. Hydraul. Eng. 2013, 140, 313–321. [Google Scholar] [CrossRef]
  4. Alvarruiz, F.; Martínez-Alzamora, F.; Vidal, A. Improving the efficiency of the loop method for the simulation of water distribution systems. J. Water Resour. Plan. Manag. 2015, 141, 04015019. [Google Scholar] [CrossRef] [Green Version]
  5. Rahal, H. A co-tree flows formulation for steady state in water distribution networks. Adv. Eng. Softw. 1995, 22, 169–178. [Google Scholar] [CrossRef]
  6. Elhay, S.; Simpson, A.R.; Deuerlein, J.; Alexander, B.; Schilders, W.H. Reformulated co-tree flows method competitive with the global gradient algorithm for solving water distribution system equations. J. Water Resour. Plan. Manag. 2014, 140, 04014040. [Google Scholar] [CrossRef] [Green Version]
  7. Schilders, W.H. Solution of indefinite linear systems using an LQ decomposition for the linear constraints. Linear Algebra Its Appl. 2009, 431, 381–395. [Google Scholar] [CrossRef] [Green Version]
  8. Todini, E.; Pilati, S. A gradient algorithm for the analysis of pipe networks. In Computer Applications in Water Supply: Volume 1—Systems Analysis and Simulation; Research Studies Press Ltd.: Somerset, UK, 1988; pp. 1–20. [Google Scholar]
  9. Elhay, S.; Deuerlein, J.; Piller, O.; Simpson, A.R. Graph Partitioning in the Analysis of Pressure Dependent Water Distribution Systems. J. Water Resour. Plan. Manag. 2018, 144, 04018011. [Google Scholar] [CrossRef] [Green Version]
  10. Simpson, A.R.; Elhay, S.; Alexander, B. Forest-core partitioning algorithm for speeding up analysis of water distribution systems. J. Water Resour. Plan. Manag. 2012, 140, 435–443. [Google Scholar] [CrossRef] [Green Version]
  11. Qiu, M.; Simpson, A.; Elhay, S.; Alexander, B. Bridge-Block Partitioning Algorithm for Speeding Up Analysis of Water Distribution Systems. J. Water Resour. Plan. Manag. 2019, 145, 04019036. [Google Scholar] [CrossRef]
  12. Rossman, L.A. EPANET 2 Users Manual, US Environmental Protection Agency; Water Supply and Water Resources Division, National Risk Management Research Laboratory: Cincinnati, OH, USA, 2000; Volume 45268. [Google Scholar]
  13. Swamee, P.K.; Jain, A.K. Explicit equations for pipe-flow problems. J. Hydraul. Div. 1976, 102, 657–664. [Google Scholar] [CrossRef]
  14. Colebrook, C.F.; Blench, T.; Chatley, H.; Essex, E.; Finniecome, J.; Lacey, G.; Williamson, J.; Macdonald, G. Correspondence. turbulent flow in pipes, with particular reference to the transition region between the smooth and rough pipe laws. (includes plates). J. Inst. Civ. Eng. 1939, 12, 393–422. [Google Scholar] [CrossRef]
  15. Benzi, M.; Golub, G.; Liesen, J. Numerical solution of saddle point problems. Acta. Num. 2005, 14, 1–137. [Google Scholar] [CrossRef] [Green Version]
  16. Qiu, M.; Alexander, B.; Simpson, A.R.; Elhay, S. A software tool for assessing the performance of and implementing water distribution system solution methods. Environ. Model. Softw. 2019, 112, 52–69. [Google Scholar] [CrossRef]
  17. Reca, J.; Martínez, J. Genetic algorithms for the design of looped irrigation water distribution networks. Water Resour. Res. 2006, 42. [Google Scholar] [CrossRef]
  18. Van Zyl, J.E.; Savic, D.A.; Walters, G.A. Operational optimization of water distribution systems using a hybrid genetic algorithm. J. Water Resour. Plan. Manag. 2004, 130, 160–170. [Google Scholar] [CrossRef]
  19. Farmani, R.; Savic, D.A.; Walters, G.A. Exnet benchmark problem for multi-objective optimization of large water systems. Modelling and Control for Participatory Planning and Managing Water Systems. In Proceedings of the IFAC Workshop, Venice, Italy, 29 September–1 October 2004. [Google Scholar]
  20. Sitzenfrei, R.; Wang, Q.; Kapelan, Z.; SaviÄ, D. Using Complex Network Analysis for Optimization of Water Distribution Networks. Water Resour. Res. 2020, 56, e2020WR027929. [Google Scholar] [CrossRef]
  21. Ostfeld, A.; Uber, J.G.; Salomons, E.; Berry, J.W.; Hart, W.E.; Phillips, C.A.; Watson, J.P.; Dorini, G.; Jonkergouw, P.; Kapelan, Z.; et al. The battle of the water sensor networks (BWSN): A design challenge for engineers and algorithms. J. Water Resour. Plan. Manag. 2008, 134, 556–568. [Google Scholar] [CrossRef] [Green Version]
  22. IEEE. IEEE Standard for Floating-Point Arithmetic. IEEE Std 754–2019 (Revision of IEEE 754-2008). 2019, pp. 1–84. Available online: https://betalearn.circuitverse.org/docs/binary-algebra/ieee-std-754.html (accessed on 22 April 2021).
  23. Tolson, B.A.; Asadzadeh, M.; Maier, H.R.; Zecchin, A. Hybrid discrete dynamically dimensioned search (HD-DDS) algorithm for water distribution system design optimization. Water Resour. Res. 2009, 45. [Google Scholar] [CrossRef]
  24. Barlow, E.; Tanyimboh, T.T. Multiobjective memetic algorithm applied to the optimisation of water distribution systems. Water Resour. Manag. 2014, 28, 2229–2242. [Google Scholar] [CrossRef] [Green Version]
  25. Ostfeld, A.; Salomons, E.; Ormsbee, L.; Uber, J.G.; Bros, C.M.; Kalungi, P.; Burd, R.; Zazula-Coetzee, B.; Belrain, T.; Kang, D.; et al. Battle of the water calibration networks. J. Water Resour. Plan. Manag. 2012, 138, 523–532. [Google Scholar] [CrossRef] [Green Version]
  26. Di Nucci, C.; Petrilli, M.; Spena, A.R. Unsteady friction and viscoelasticity in pipe fluid transients. J. Hydraul. Res. 2011, 49. [Google Scholar] [CrossRef]
  27. Pasculli, A. FD-FEM 2D Modelling of a local water flow. Some numerical results. Alp. Mediterr. Quat. 2008, 21, 215–228. [Google Scholar]
Figure 1. The friction factors calculated for the turbulent flow region using the proposed equation in Equation (27) and Swamee–Jain equation for: (a) ϵ / D = 10 2 , (b) ϵ / D = 10 3 , (c) ϵ / D = 10 4 , and (d) ϵ / D = 10 5 .
Figure 1. The friction factors calculated for the turbulent flow region using the proposed equation in Equation (27) and Swamee–Jain equation for: (a) ϵ / D = 10 2 , (b) ϵ / D = 10 3 , (c) ϵ / D = 10 4 , and (d) ϵ / D = 10 5 .
Water 13 01163 g001
Figure 2. The infinity norm of the relative head difference h m + 1 h m h m + 1 and infinity norm of the continuity residuals after applied the proposed head formulation to each of the six case study networks N 1 N 6 using the Newton method as shown in Figure 2, at sub-figures (ef).
Figure 2. The infinity norm of the relative head difference h m + 1 h m h m + 1 and infinity norm of the continuity residuals after applied the proposed head formulation to each of the six case study networks N 1 N 6 using the Newton method as shown in Figure 2, at sub-figures (ef).
Water 13 01163 g002
Figure 3. The convergence property of the proposed head formulation applied to the Balerma Network: (a) the head iterates at Junction 776 in Balerma Network and source elevation at Reservoir E using the Newton method; (b) the head iterates at Junction 776 in Balerma Network and source elevation at Reservoir E using the damped Newton method.
Figure 3. The convergence property of the proposed head formulation applied to the Balerma Network: (a) the head iterates at Junction 776 in Balerma Network and source elevation at Reservoir E using the Newton method; (b) the head iterates at Junction 776 in Balerma Network and source elevation at Reservoir E using the damped Newton method.
Water 13 01163 g003
Figure 4. The infinity norm of the relative head difference h m + 1 h m h m + 1 and infinity norm of the continuity residuals after applied the proposed head formulation to each of the six case study networks N 1 N 6 using the damped Newton method as shown in Figure 4, at sub-figures (ef).
Figure 4. The infinity norm of the relative head difference h m + 1 h m h m + 1 and infinity norm of the continuity residuals after applied the proposed head formulation to each of the six case study networks N 1 N 6 using the damped Newton method as shown in Figure 4, at sub-figures (ef).
Water 13 01163 g004
Figure 5. The spatial distribution of the different levels of the relative differences between flow results of the GGA and the proposed head formulation applied to (a) network N 2 , (b) network N 3 , (c) network N 4 , and (d) network N 5 .
Figure 5. The spatial distribution of the different levels of the relative differences between flow results of the GGA and the proposed head formulation applied to (a) network N 2 , (b) network N 3 , (c) network N 4 , and (d) network N 5 .
Water 13 01163 g005
Figure 6. The differences between the tank level found by applying the head formulation and the GGA to the network Net3 in EPANET for the network three tanks, as shown in sub-figures (ac).
Figure 6. The differences between the tank level found by applying the head formulation and the GGA to the network Net3 in EPANET for the network three tanks, as shown in sub-figures (ac).
Water 13 01163 g006
Figure 7. The relative differences between the flows in each pipe of the EPANET example 3 found by the GGA and the proposed head formulation over 24 h extended-period simulation: (a) temporal variation over the 24 h extended-period simulation, (b) spatial distribution of different levels of relative difference of pipe flows at t = 09:00, (c) spatial distribution of different levels of relative difference of pipe flows at t = 13:00, (d) spatial distribution of different levels of relative difference of pipe flows at t = 17:00, and , (e) spatial distribution of different level of relative difference of pipe flows at t = 18:00.
Figure 7. The relative differences between the flows in each pipe of the EPANET example 3 found by the GGA and the proposed head formulation over 24 h extended-period simulation: (a) temporal variation over the 24 h extended-period simulation, (b) spatial distribution of different levels of relative difference of pipe flows at t = 09:00, (c) spatial distribution of different levels of relative difference of pipe flows at t = 13:00, (d) spatial distribution of different levels of relative difference of pipe flows at t = 17:00, and , (e) spatial distribution of different level of relative difference of pipe flows at t = 18:00.
Water 13 01163 g007
Table 1. Benchmark networks summary.
Table 1. Benchmark networks summary.
NetworkNo. of PipesNo. of NodesNo. of SourcesNo. of Forest Pipes
N 1 4544434288
N 2 9348488361
N 3 246518903429
N 4 4021355711566
N 5 14,83012,52372932
N 6 157,044150,630445,736
Table 2. Convergence property of the proposed head formulation applied to the six case study networks.
Table 2. Convergence property of the proposed head formulation applied to the six case study networks.
NetworkNum. of Iter. h m + 1 h m + 1 h m + 1 σ C σ E σ f
N 1 16 7 × 10 7 1 × 10 7 1 × 10 15 1 × 10 15
N 2 11 7 × 10 7 9 × 10 6 1 × 10 15 2 × 10 15
N 3 18 7 × 10 7 2 × 10 6 1 × 10 15 1 × 10 15
N 4 14 9 × 10 7 5 × 10 6 1 × 10 15 1 × 10 15
N 5 18 9 × 10 7 1 × 10 6 1 × 10 15 1 × 10 15
N 6 15 9 × 10 7 2 × 10 6 1 × 10 15 3 × 10 15
Table 3. The number of iterations and the wall-clock time (second) required to perform the pre-processing, iterative phase, and the post-processing operations for the GGA and the proposed head formulation applied to each of the six case study networks (the number in the bracket indicates the per-iteration run time required to execute the iterative phase of the Newton method).
Table 3. The number of iterations and the wall-clock time (second) required to perform the pre-processing, iterative phase, and the post-processing operations for the GGA and the proposed head formulation applied to each of the six case study networks (the number in the bracket indicates the per-iteration run time required to execute the iterative phase of the Newton method).
NetworkMethodsNum. of Iter.Preproc.Iterative PhasePostProc.Total
N 1 GGA60.0030.003 (0.0005)0.0010.007
Head160.0040.009 (0.0005)0.0010.014
N 2 GGA60.0050.006 (0.001)0.0010.012
Head110.0070.013 (0.001)0.0010.021
N 3 GGA90.0110.025 (0.003)0.0010.037
Head180.0140.053 (0.003)0.0010.068
N 4 GGA60.0320.027 (0.005)0.0010.060
Head140.0240.070 (0.005)0.0010.095
N 5 GGA80.0620.143 (0.018)0.0020.207
Head180.0890.367 (0.02)0.0010.457
N 6 GGA81.931.78 (0.22)0.023.73
Head152.143.61 (0.26)0.0065.75
Table 4. The detailed statistics of the absolute head differences (m) between the GGA and the proposed head formulation, | h G G A h h e a d | , applied to six case study networks.
Table 4. The detailed statistics of the absolute head differences (m) between the GGA and the proposed head formulation, | h G G A h h e a d | , applied to six case study networks.
NetworkMinMeanMedianMaxstd.dev
N 1 1.00e-48.98e-27.89e-23.07e-15.99e-2
N 2 0.003.12e-32.70e-39.90e-32.14e-3
N 3 5.00e-42.04e-12.12e-13.09e-15.27e-2
N 4 1.30e-21.20e-11.19e-12.37e-15.13e-2
N 5 4.65e-25.88e-25.88e-28.15e-23.72e-3
N 6 3.00e-32.31e-21.61e-29.22e-21.62e-2
Table 5. The nodal pressure violation in some of the solutions of least-cost design of network N 1 using the proposed head formulation.
Table 5. The nodal pressure violation in some of the solutions of least-cost design of network N 1 using the proposed head formulation.
(a) 1,940,923 design found in Tolson et al. [23]
JunctionElevationTotal HeadPressure Deficit
Junc 323.9043.6219.722.76e-1
Junc 5921.0040.9219.928.18e-2
Junc 15156.5076.3919.891.14e-1
Junc 23387.17107.0819.919.08e-2
Junc 27073.7093.5919.891.07e-1
Junc 28175.0094.9219.928.46e-2
Junc 33275.7095.6019.909.69e-2
Junc 35980.70100.6219.928.25e-2
Junc 36380.50100.4119.919.23e-2
Junc 39456.4076.3019.901.01e-1
Junc 39780.50100.4119.918.62e-2
Junc 39880.50100.4019.909.84e-2
Junc 40181.70101.7019.999.00e-4
(b) 1,920,656 design found in Barlow and Tanyimboh [24]
JunctionElevationTotal HeadPressure HeadPressure Deficit
Junc 323.9043.7219.821.81e-1
Junc 13560.0079.9519.955.12e-2
Junc 15045.0064.8219.821.76e-1
Junc 15255.0074.9919.992.38e-2
Junc 20195.00114.9519.955.03e-2
Junc 23387.17107.1119.945.58e-2
Junc 28175.0094.9619.964.09e-2
Junc 35980.70100.5919.891.07e-1
Junc 36380.50100.4119.918.83e-2
Junc 37469.5089.4119.918.75e-2
Junc 39780.50100.4119.919.10e-2
Junc 39880.50100.3919.891.12e-1
Junc 40181.79101.6919.895.70e-3
Junc 179,00160.0079.9619.964.00e-2
Table 6. The number of pipes in different bins of the relative difference between the Flow results of the GGA and the proposed head formulation applied to six case study networks.
Table 6. The number of pipes in different bins of the relative difference between the Flow results of the GGA and the proposed head formulation applied to six case study networks.
N 1 N 2 N 3 N 4 N 5 N 6
Same302 (66.52%)371 (39.72%)438 (17.77%)1567 (38.97%)2951 (19.90%)45,788 (29.16%)
0–1%146 (32.16%)390 (41.76%)858 (34.81%)2312 (57.50%)8113 (54.71%)93,859 (59.77%)
1–10%(1.32%)138 (14.78%)671 (27.22%)106 (2.64%)1926 (12.99%)11,210 (7.14%)
10–20%0 (0.00%)15 (1.61%)193 (7.83%)4 (0.10%)295 (1.99%)1684 (1.07%)
20–30%0 (0.00%)2 (0.21%)85 (3.45%)8 (0.20%)179 (1.21%)940 (0.60%)
30–40%0 (0.00%)2 (0.21%)43 (1.74%)1 (0.02%)125 (0.84%)648 (0.41%)
40–50%0 (0.00%)0 (0%)26 (1.05%)0 (0.00%)117 (0.79%)453 (0.29%)
>50%0 (0.00%)16 (1.71%)151 (6.13%)23 (0.57%)1124 (7.58%)2462 (1.57%)
Flow reversal0 (0.00%)7 (0.75%)45 (1.83%)12 (0.30%)448 (3.02%)730 (0.46%)
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Qiu, M.; Ostfeld, A. A Head Formulation for the Steady-State Analysis of Water Distribution Systems Using an Explicit and Exact Expression of the Colebrook–White Equation. Water 2021, 13, 1163. https://doi.org/10.3390/w13091163

AMA Style

Qiu M, Ostfeld A. A Head Formulation for the Steady-State Analysis of Water Distribution Systems Using an Explicit and Exact Expression of the Colebrook–White Equation. Water. 2021; 13(9):1163. https://doi.org/10.3390/w13091163

Chicago/Turabian Style

Qiu, Mengning, and Avi Ostfeld. 2021. "A Head Formulation for the Steady-State Analysis of Water Distribution Systems Using an Explicit and Exact Expression of the Colebrook–White Equation" Water 13, no. 9: 1163. https://doi.org/10.3390/w13091163

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop