Discussion of “Jacobian Matrix for Solving Water Distribution System Equations with the Darcy-Weisbach Head-Loss Model” by Angus Simpson and Sylvan Elhay

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Discussion of “Jacobian Matrix for Solving Water Distribution System Equations with the Darcy-Weisbach Head-Loss Model” by Angus Simpson and Sylvan Elhay Dejan Brkić

The main purpose of this discussion is to recommend more appropriate formulas for the calculation of the Darcy friction factor, λ which are more suitable for computation of the derivatives in a Jacobian matrix, and also to clarify a few points in the calculation of the flowdependent Jacobian. The Darcy-Weisbach formula contains the Darcy friction factor, λ which in general is dependent on the Reynolds number and, therefore, in general on flow, i.e. λ(Re)→λ(Q). This fact can be incorporated in the Jacobian matrix [Eq. (6) in the discussed paper] as clearly shown in Simpson and Elhay (2011). On the other hand, flow, Q, in the formula for flow friction, λ, is not treated as a variable in EPANET (Computer software, Washington, DC, U.S. Environmental Protection Agency) during the computation of the Jacobian matrix (probably following the analogy with the Hazen-Williams scheme). Use of the Darcy friction factor dependent on flow λ=λ(Q) in the Jacobian matrix is more accurate but also more complicated. A certain balance between accuracy and complexity in engineering systems has to be carefully evaluated. The procedure shown in the discussed paper is more accurate mathematically and also can be more easily implemented in a computer code. This increased computing burden of the improved procedure has to be Text DB Click here to download Manuscript: Discussion DB (Simpson and Elhay).doc carefully investigated where the ratio between complexity and accuracy will decide whether this method will be used in every day engineering practice (Giustolisi et al 2011).
The pressure drop model Δp (Pa) for the Darcy-Weisbach formula Δp=8·L·ρ·Q 2 ·λ/(π 2 ·D 5 ) is used in this discussion instead of the head-loss model ΔH (m) which was used in the discussed paper where the connection between them is simple, i.e. Δp=ρ·g·ΔH. Also, the Darcy friction factor noted as λ is used as in Colebrook (1939) where the Darcy friction is noted as f in the discussed paper (one should not confuse the Darcy friction factor used in this discussion and the discussed paper with the Fanning friction factor which is numerically equal to one-fourth the Darcy friction). The Darcy friction factor λ for the turbulent zone (turbulent flow is the most common in pipe networks) can be calculated using the implicit Colebrook Eq. (16), using some of the explicit approximations to the Colebrook equation (Brkić 2011a) or using some other suitable formulas [EPANET uses the approximation by Swamee and Jain (1976) which is shown in the next formula on the far right side of the equation].
The first derivative of the Darcy-Weisbach equation (appropriate terms in the Jacobian matrix [Eq. (6) in the discussed paper]) has to be calculated using the fact that the Darcy friction factor, λ contains flow, Q, which now has to be treated as a variable, i.e., λ=λ(Q) and not as constant anymore (e.g. flow is treated as constant in Brkić (2011b) Table 2 and Eq. (2) in the original paper] is nothing more than the Darcy-Weisbach formula where the suitable Darcy friction factor λ=64/Re is included (Eq. 18). The fact that hydraulic resistance, r, from Eq. (18) is independent of flow [as shown in Eq. (2) of the discussed paper] implies that a water network in which laminar flow takes place (in all pipes) can be solved using a non-iterative procedure. On the other hand, the Darcy friction factor, λ in the laminar zone depends only on the Reynolds number, Re, which further depends on flow Q. The fact that laminar resistance, r is independent of flow, Q, is not crucial for the presented calculation [where λ=λ(Q) and r≠r(Q)]. The use of separate equations for hydraulic resistance, r or at least separate comments for the laminar and turbulent regimes in Eqs. (1-3) of the discussed paper should be avoided (better to say the parameter, λ should always be used and not r). However, in laminar, turbulent zone and also for the unstable transient zone between them (Table 2 and 4 of the discussed paper) it is certain that the Darcy friction factor, λ depends on the Reynolds number [λ=λ(Re, ε/D)] and subsequently in general on flow [where Re=Re(Q)]. This is not true for the fully turbulent regime (fully developed flow) where the Darcy friction factor, λ depends only on relative roughness ε/D and not on the Reynolds number Re and where, therefore the influence of the flow Q can be neglected).
Finally, the flow dependent Jacobian, as recommended by Simpson and Elhay (2011), and the Darcy-Weisbach scheme should be used. Furthermore, EPANET also should be revised in light of the following discussion. In the laminar zone the Darcy friction factor should be calculated using λ=64/Re. For the turbulent zone, the Swamee and Jain approximation [used in EPANET, in the discussed paper and in Spiliotis and Tsakiris (2011) Brkić (2011a).
For example, the approximation proposed by Romeo et al. (2002) with a maximal error up to 0.1345% is more accurate. Also, approximations by Buzzelli (2008) and Serghides (1984) or by Zigrang and Sylvester (1982) are very accurate with an error of no more than 0.1385%. Vatankhah and Kouchakzadeh (2008) also show very accurate approximation with a relative error up to 0.1472%. There are at least ten approximations more accurate than the Swamee-Jain equation. These very accurate approximations are sometimes complex but they can be simply included in a computer code (Matlab, Excel, etc.) without causing significant additional computational time. Also, between the laminar and turbulent zones it is more convenient to model it as a sharp jump (sudden transition to turbulence) instead of using a cubic spline (explained in the discussed paper in Table 2; also instead of a cubic spline, the simple equation λ=0.0025·Re 0.33 is more convenient). At last but not less important, Churchill (1977) also provides a single formula which covers all three cases in the discussed paper ( Table 2 in the discussed paper of Simpson and Elhay (2011), i.e. used for the laminar and turbulent zone as well as for the unstable transient zone between them as described in EPANET). The equation by Churchill (1977) perfectly fits the laminar regime and provides a reasonable transition to turbulence (cubic spline in the discussed paper). Finally it approximates the implicit Colebrook function (suitable only for the turbulent zone) at least equally as accurate as the Swamee-Jain approximation. Simpson and Elhay (2011) use Matlab and hence for the readers" convenience, Matlab code for the well known equation by Churchill (1977) is given: