Quadratic approximations for pipe friction

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³/s]

Q₁

flow rate at lower end of approximation interval [m³/s]

Q₂

flow rate at upper end of approximation interval [m³/s]

Re

Reynolds number

V

average fluid velocity [m/s]

α

constant part of Hazen–Williams formula with respect to flow

ε

relative error

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.

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

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

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.

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⁻⁶ m²/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⁵, 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)).

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

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⁻⁸ in seven (quadratic) or eight (Hazen–Williams) iterations.

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⁻⁸. 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.

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.