## Abstract

This paper develops quadratic approximations to the classical models of Darcy–Weisbach and Hazen–Williams for frictional head loss in water pipes. Key elements of the technique are selecting the approximation domain and minimizing the relative error. A comparison for individual pipes over a range of diameter, roughness, and Reynolds number showed that the approximation provides accuracy consistent with the experimental error in the classical equations. Two benchmark water distribution networks are also considered. In these systems, flows and pressures computed using the approximation were consistent with simulations using the classical models. Applications of the approximation include mathematical optimization problems where polynomial expressions are desirable.

- modeling
- pipe friction
- water networks

## NOTATION

*a*- dimensional coefficient
*b*- dimensional coefficient
*f*- Darcy friction factor [–]
*g*- gravitational acceleration [m/s/s]
*h*_{f}- frictional head loss [m]
- approximation of frictional head loss [m]
*k*_{s}- roughness height [m]
*A*−*E*- constants used to compute
*a*and*b*for Hazen–Williams *C*- Hazen–Williams
*C*value *D*- pipe diameter [m]
*F*- integral of relative errors
*L*- pipe length [m]
*Q*- flow rate [m
^{3}/s] *Q*_{1}- flow rate at lower end of approximation interval [m
^{3}/s] *Q*_{2}- flow rate at upper end of approximation interval [m
^{3}/s] - Re
- Reynolds number
*V*- average fluid velocity [m/s]
*α*- constant part of Hazen–Williams formula with respect to flow
*ε*- relative error

## INTRODUCTION

The equations describing frictional head loss in water pipes are well known and widely applied. Head loss is described by the Darcy–Weisbach equation, or in certain flow regimes using empirical formulas such as those of Hazen–Williams (Liou 1998; Christensen *et al*. 2000). Today's network modeling software packages such as EPANET (Rossman 2000) implement these approaches and are widely used for planning and operation of water systems. However, applications for network models are becoming more sophisticated to include data assimilation and real-time modeling (Hutton *et al*. 2014), optimization (Broad *et al*. 2010), and inverse problems. In these applications, the properties of the classical head loss equations with respect to complexity, differentiability, and smoothness limit the techniques that can be applied.

Techniques for solving head loss equations depend on the problem studied and friction loss model. Applications in network simulation often use a Newton-type algorithm that requires a Jacobian matrix. The Jacobian can be expensive to compute particularly for the Darcy–Weisbach model (Simpson & Elhay 2011). With the Hazen–Williams equations, the Jacobian can become singular near-zero flow causing simulations to fail (Elhay & Simpson 2011). In mathematical optimization problems, conservation of energy along pipes often appears as a constraint. The challenge of using the classical equations as constraints is that the rational exponents appearing in the Hazen–Williams formula or in explicit approximations for Darcy–Weisbach may have singularities at zero flow. A quadratic approximation for head loss allows to formulate a variety of network optimization problems as bounded polynomial optimization problems, non-convex problems whose global optima can be approximated by sequences of semi-definite optimization problems, that is, a class of convex optimization problems which can be solved efficiently (Lasserre 2001).

Numerous relationships for pipe head loss have been proposed. The empirical formulas of Hazen and Williams and Manning give explicit equations, but apply over a limited range of flows. In contrast, the Darcy–Weisbach equation for head loss is accurate and generally applicable, but requires a friction factor that depends on the flow regime. For turbulent flows the factor may be computed from the implicit Colebrook–White formula, or in the limit of fully rough flow from the explicit formula of Prandtl and Karman (White 1999). Many explicit approximations of the Darcy friction factor are also available, most notably those of Swamee & Jain (1976) for popularity and Sonnad & Goudar (2006) for accuracy. A recent summary of approximations for the friction factor is given by Giustolisi *et al*. (2011).

The need to simplify the Darcy–Weisbach formula for use in optimization problems has been recognized by a few earlier workers. Valiantzas (2008) proposes an explicit power law approximation that provides accuracy of ±4.5%. The fractional exponents in this approximation cause the same numerical difficulties as the Hazen–Williams formula. Burgschweiger *et al*. (2009) propose a globally smooth explicit approximation for pipe head loss that includes square root functions. The approximation is asymptotically consistent with classical equations, but as will be shown subsequently, may underestimate friction losses in certain situations. Gleixner *et al*. (2012) use the Prandtl–Karman relation which can also underestimate friction losses and has a singularity in the inverse derivative at zero flow.

The contributions of the present paper are quadratic approximations for friction losses in water pipes that minimize the relative error and apply to water systems parameterized by either the Darcy–Weisbach or Hazen–Williams models. The primary motivation for this approximation is to model head loss with high accuracy while providing a mathematical form, i.e., a polynomial approximation with integer powers, which is well-behaved when embedding it in constraints of optimization problems over water networks. Some examples of such problems include network design (Alperovits & Shamir 1977), valve placement (Reis *et al*. 1997), pump scheduling (Sterling & Coulbeck 1975), and valve setting (Vairavamoorthy & Lumbers 1998) among others. In particular, the work of Bragalli *et al*. (2011) on the network design problem addresses the numerical difficulties of rational exponents in the head loss function by fitting a quintic near-zero. An additional effect of the proposed approximation is a lower complexity for computing the head loss function and its derivative. The reduced complexity may also prove beneficial where hydraulic simulators make many runs guided by genetic algorithms (Savic & Walters 1997), simulated annealing (McCormick & Powell 2004), and other heuristics used for water network optimization.

## BACKGROUND

The head loss due to friction in water distribution systems is well studied and may be described by the Darcy–Weisbach equation (White 1999)
1where *f* is the Darcy friction factor; *L* is the pipe length; *D* is the diameter; *V* is the average fluid velocity; and *g* is the gravitational constant, all in consistent units.

The method for computing the Darcy friction factor varies according to whether the flow is laminar or turbulent. In the laminar regime, Poiseuille flow occurs and the friction factor is , where is the Reynolds number, and is the kinematic viscosity. Thus, under laminar conditions .

In the turbulent regime, the friction factor's behavior with respect to varies according to flow conditions as diagramed by Moody (1944). For fully rough turbulent flow, the factor is given by the classical formula of Prandtl–Karman
2where *k _{s}* is the roughness height and

*D*is the diameter. In Equation (2),

*f*is constant with respect to Re and thus . In turbulent flows that are not fully rough

*f*depends in a complex way on Re. In either case, the value of

*f*may be found graphically from the Moody diagram, or by inverting the implicit relationship of Colebrook (White 1999) 3The value of

*f*may also be estimated using one of the many explicit approximations in the literature (see Giustolisi

*et al*. (2011) for a recent summary).

In the transition from laminar to fully rough flow, the Hazen–Williams formula also applies. In SI units, friction loss based on the Hazen–Williams formula is computed from American Water Works Association (2004)
4where *h _{f}* is the head loss due to friction in meters;

*C*is the Hazen–Williams coefficient;

*D*is the inside diameter of the pipe in meters;

*L*is the length of pipe in meters; and

*Q*is the discharge in cubic meters per second. The Hazen–Williams formula applies only to water flowing in pipes 5 cm in diameter or larger and at speeds less than 3 m/s (Franzini & Finnemore 1997). A more detailed discussion of the region of application for the Hazen–Williams formula is given by Liou (1998) and discussions.

Although Equations (3) and (4) are often treated as exact, they are fitted to experimental results and therefore contain some error. Within the flow range where tests were performed, Liou (1998) reports errors from −15 to 10% in the Hazen–Williams formula depending on pipe material. On the Darcy–Weisbach side, the Colebrook equation for *f* has an average error of 11% with respect to Colebrook's data (Yoo & Singh 2005).

## QUADRATIC APPROXIMATION FOR PIPE FRICTION

In some applications such as real-time simulation or mathematical optimization, a simpler relationship than Equation (1) with Equation (3), or an explicit approximation for *f*, that provides a sufficiently accurate representation of the underlying physics is desirable. The limiting case of fully rough flow suggests a squared function. Indeed, the head-flow curve for a pipe shows strongly quadratic behavior in the turbulent range. Thus, the following quadratic approximation for the pipe head loss curve is proposed:
5where *a* and *b* are unknown dimensional coefficients and *Q* is the volumetric flow rate. The intercept is dropped to enforce the origin as a point on the curve. The linear term is retained to provide a degree of freedom for improving the fit, and in recognition of the fact that turbulent flows outside of the fully rough range are not purely a squared function of the flow.

The derivative or reciprocal derivative of the head loss function is needed in simulation and optimization applications. Using a quadratic, the derivative expression is easily found and well behaved 6Notice that from the left and the right. The reciprocal of the derivative is non-zero provided for .

In developing an approximation function one must consider several factors:

the functional form of the approximation;

a merit function for evaluating candidates;

the relevant range of the original function;

the method for finding coefficients.

This work uses the functional form of Equation (5) whether Darcy–Weisbach or Hazen–Williams is used as an underlying model. Coefficients are selected to minimize the ‘relative error’ (MRE) over the approximation range. When the coefficients of Equation (5) are found in this way, the approximation is termed a ‘MRE quadratic’. The exact merit function, approximation range, and method of finding coefficients differ slightly between head loss models as discussed below.

### Darcy–Weisbach

The preferred method for modeling head loss on water networks is the Darcy–Weisbach law of Equation (1). If the friction factor *f* were constant the head loss curve would be purely quadratic; this is the primary motivation for developing a quadratic approximation. The fit proposed here minimizes the sum of squared relative errors over the turbulent flow range . Owing to the implicit form of Colebrook's equation for *f*, and the complexity of explicit approximations, the fit is obtained numerically by solving the least squares problem
7By using as the weight for each point, the coefficients minimize relative rather than absolute errors. A vector of flow rates having length *N* defines the approximation interval. The lower end of the approximation interval corresponds to . At the upper end of the interval, *Q _{N}* corresponds to the maximum flow rate expected in the application based on a preliminary examination of the network. A sufficiently large number of points should be used so that the values of the fitted coefficients

*a*and

*b*do not depend on the number of points.

A quadratic fit to the Darcy–Weisbach head loss curve in a sample pipe is shown in Figure 1. The sample pipe has length 1,000 m, diameter 100 mm, and roughness height 0.3 mm. Considering water with a kinematic viscosity of 1.002 × 10^{−6} m^{2}/s, a MRE quadratic approximation using 2,000 points over the range yielded coefficients and for use with flow rates in L/s. The relative error ranges from −1.9 to 8% with the largest relative errors occurring at low flow rates.

Because the quadratic fit is developed for a particular flow range, the accuracy of the approximation outside of the range is of interest. In the laminar region, head losses are linear and the quadratic approximation over the turbulent range does not provide a very good estimate (Figure 2). However, the errors noted here do not affect the flow distribution in the networks studied because head losses are low. Figure 2 also illustrates why the transition flow range is excluded from the approximation: trying to minimize relative errors through the jump in the transition range sacrifices accuracy everywhere else along the curve. On the upper end of the flow range, the shape of the relative error curve suggests good agreement with the theoretical curve beyond the range selected for approximation. In the example of Figure 1, applying the approximation at Re = 2 × 10^{5}, twice the fitted value, incurs an error of 2.2%.

Visual inspection of the Moody diagram shows that the friction factor varies over a smaller range as the relative roughness of the pipe increases. This smaller range implies that rougher pipes, such as those in aging distribution systems, are well approximated by a quadratic function. Indeed the fit improves as pipes become rougher (Figure 3(a), 3(c), and 3(e)).

Also shown in Figure 3 is the curve given by Prandtl–Karmman's formula and the approximation due to Burgschweiger *et al*. (2009). The Burgschweiger approximation is given by
8where the constants *r*, *b*, and *c* are found from pipe properties and the values of *a* and *d* are chosen to match a desired slope at . In the comparisons shown here, *a* and *d* are chosen to match the slope at which occurs under laminar flow according to the law of Hagen–Poiseuille. For pipes with roughness heights of 10 and 1 mm, both a quadratic approximation, the law of Prandtl–Karman, and Burgschweiger's approximation fall directly on the ideal curve. As a roughness height decreases, deviation of the approximations from the ideal curve becomes visible.

### Hazen–Williams

The Hazen–Williams formula for pipe head loss applies over a limited range; outside this range it does not reflect the underlying physics and has considerable errors compared to experimental results (Liou 1998; Franzini & Finnemore 1997). These deficiencies notwithstanding, the equation remains in wide use because the explicit formulation is easy to compute. Since existing models of real networks often use Hazen–Williams, a quadratic approximation is provided here.

A quadratic approximation to the Hazen–Williams formula can be found numerically as described above for Darcy–Weisbach. However, with its explicit functional form a quadratic approximation to the Hazen–Williams formula can also be obtained analytically by minimizing the integral of relative errors
9where and specify the range of interest and *α* is the part of Equation (4) that does not vary with flow rate . The integral of Equation (9) is found and the resulting form is differentiated by *a* and *b*. These differentials are set equal to zero, giving a 2 × 2 system of linear equations. Solving the system yields expressions for *a* and *b* as functions of the approximation interval
10where
11A quadratic approximation to the Hazen–Williams curve for a sample pipe is shown in Figure 4 along with the relative error curve. Several points of interest are noted:

the quadratic fit provides errors within 10% over the fitted range of flows;

since the relative error is minimized, there is a singularity in the relative error at zero flow;

the shape of the relative error curve suggests a method for choosing the approximation interval.

Choosing the approximation interval is a key step because the values of *a* and *b* depend on the interval. The value of should be chosen to reflect the maximum flow likely to be encountered in the application. A value can be chosen through investigation of the network or selection of a maximum velocity or power dissipation in each pipe. Once a value for is established, may be selected so that the minimum of the relative error curve sits at an error tolerance. The location of the minimum relative error is found from calculus
12With Equation (12) and an error tolerance, ε_{tol}, is found at the root of
13This approach to setting the approximation interval is conservative in that the worst underestimate of head losses is specified. Head losses which are overestimated occur only at the extreme ends of the interval and relative errors are largest for small flows. A quadratic approximation to the Hazen–Williams head loss curve is shown for a range of diameters, flows, and roughnesses in Figure 3(b), 3(d), and 3(f).

## COMPARISON ON PIPE NETWORKS

A numerical model was created to study the accuracy of the MRE quadratic approximation for networks of pipes. The model solves the set of nonlinear equations for steady-state network hydraulics using the global gradient algorithm of Todini & Pilati (1988). For the MRE quadratic, Jacobian terms are given by Equation (6). For the Darcy–Weisbach model, the Jacobian matrix given by Simpson & Elhay (2011) is used. Note that the coefficients of the MRE quadratic approximation are pre-computed, and so are treated as constants by the simulator. Because the goal here was a comparison of accuracies, the implementation was not optimized for speed. In the following subsections, the simulation model is applied on two networks from the literature.

### Twenty-five node network

The first network considered here is the benchmark network of Sterling & Bargiela (1984) that has also been studied by Jowitt & Xu (1990), Reis *et al*. (1997), Vairavamoorthy & Lumbers (1998), and Nicolini & Zovatto (2009) among others. The network layout is shown in Figure 5. With three sources, 22 junctions and 37 links, the network is highly connected. Head loss parameterization is given using Hazen–Williams coefficients.

Quadratic approximations for head loss in each pipe were developed using the methods given above with parameter corresponding to a velocity of 2 m/s and . For convenience of other workers, network information including coefficients of a quadratic approximation is provided in Tables 1 and 2.

The accuracy of the quadratic approximation was assessed by comparing flows and pressures from a simulation using the approximation to the solution obtained from Hazen–Williams coefficients (Figure 6). Differences between pressures obtained by the models were small, with the largest relative error of 0.62% or 0.24 m at node 14. In the flow solution relative errors had a median of 0% and relative errors in the middle 95% of the data ranged from −12 to 20%. In the remaining 5% of pipes absolute errors ranged from −1.9 to 0.48 L/s. As expected, high relative errors occurred at low flows and had a negligible impact on the overall mass balance.

For this small network, the MRE quadratic approximation showed similar computational performance to the conventional Hazen–Williams model. On a Windows laptop with 2.2 GHz processor and 8 Gb RAM, simulations using the quadratic converged in 5 ms (average of 100 runs) and simulations using the Hazen–Williams coefficients converged in 6 ms. Both simulations converged to a tolerance of 10^{−8} in seven (quadratic) or eight (Hazen–Williams) iterations.

### EXNET

To illustrate the applicability of the MRE quadratic approximation on another network, the EXNET water system first mentioned by Farmani *et al*. (2005) is also considered (Figure 7). This system is based on a real network and is suggested as a benchmark by University of Exeter Centre for Water Systems (2012). Head losses are modeled by the Darcy–Weisbach method. For the purposes of the comparison here, the valves in the network file available online were converted to pipes of length 10 m. A quadratic approximation for the head loss in each pipe was created by choosing at and based on 2 m/s.

Flows and pressures calculated using the quadratic approximation compared favorably to the EPANET solution (Figure 8). In the pressure solution, errors ranged from −4.4 to 0.63% and averaged −0.46%. In the flow solution, relative errors had a median of zero and relative errors in the middle 95% of the data ranged from −11 to 3.6%. In the remaining 5% of pipes, flow rates ranged from −15 to 71 L/s and absolute errors ranged from −9.7 to 1.3 L/s. As seen in the figure, these errors did not affect the overall mass balance.

From a computational perspective, the quadratic approximation showed similar performance to the Darcy–Weisbach model. One hundred simulations converged in an average 15 s for the quadratic model and 17 s for the Darcy–Weisbach model. The quadratic model required 10 iterations and the Darcy–Weisbach model required 11 iterations to converge to a tolerance of 10^{−8}. Profiling of the code indicated that the difference in time was attributable both to the computation of the head loss terms and to the additional iteration for the Darcy–Weisbach model. These results suggest benefits of the quadratic approximation reside not in the computational speed but in the functional form.

## CONCLUSIONS

This paper has proposed a quadratic approximation for pipe head loss curves which minimizes the relative approximation error. The approximation applies to pipes and networks parameterized by either the Darcy–Weisbach or Hazen–Williams models. By focusing on the relative rather than absolute error, the approximation stays consistent over the range of possible flows. The accuracy of the approximation depends on the range of pipe flow rates and the pipe roughness. Because rougher pipes develop fully turbulent flow at lower Reynolds numbers, the friction factor varies over a smaller range and a quadratic approximation is more accurate. This result means that the approximation is especially well suited to older water systems where pipes tend to be rougher.

Evaluations of the MRE quadratic approximation proposed here included comparison with other approximations and use in a simulation model. For the Darcy–Weisbach parameterization, MRE quadratic approximations were compared to other essentially quadratic forms. For rough pipes, all approximations show good agreement with the theoretical curve. As pipes become smoother, the MRE quadratic shows closer agreement with the ideal curve. The MRE quadratic approximation was also implemented in a network simulation model and tested on two networks from the literature. Results showed that the MRE quadratic yielded a flow and pressure distribution consistent with results obtained by simulating with the original equations.

In the networks simulated here, the MRE quadratic provided a slight but insignificant improvement in simulation speed compared to the classical equations. Speedup due to the MRE quadratic may become more pronounced for larger networks, especially where thousands of model runs are used to explore the search space of a combinatorial problem such as network expansion or pump scheduling.

Pipe friction models based on quadratic forms open up several new avenues for the theory and practice of water network analysis. Some operational and planning problems can now be formulated as polynomial optimization problems. When additionally exploiting the sparse network structure of an urban water network, this may allow advantage to be taken of sparse semi-definite optimization techniques to approximate the globally optimal solutions of these operational and planning problems.

- First received 30 June 2014.
- Accepted in revised form 14 November 2014.

- © IWA Publishing 2015

Sign-up for alerts