OpenWaterAnalytics / EPANET

The Water Distribution System Hydraulic and Water Quality Analysis Toolkit
MIT License
279 stars 204 forks source link

Incorrect coefficient used for linear head loss approximation at low flow #662

Closed LRossman closed 2 years ago

LRossman commented 2 years ago

To prevent its Newton-Raphson Jacobian matrix from becoming ill-conditioned when solving network hydraulics, EPANET replaces the nonlinear Hazen-Williams (and Chezy-Manning) head loss equation with a linear approximation whenever pipe flow is below the flow q* that produces a head loss gradient of RQtol(equal to 1e-7). From 0 to q* the head loss should equal (RQtol / Hexp) * q, where Hexpis the exponent used in the H-W (or C-M) equation. However the 1/Hexp factor is missing in the current code (see line 556 in hydcoeffs.c). This same bug occurs for emitter head loss (missing a 1/Qexp factor) and for pressure dependent demand head loss (missing a Pexp factor).

The derivation of the (RQtol / Hexp) coefficient for low pipe flow is as follows. The H-W head loss equation can be expressed as:

    hloss = R * q^Hexp

where hlossis head loss, Ris a resistance factor, q is flow rate and Hexpis 1.852. Its derivative with respect to flow, hgrad, is:

    hgrad = Hexp * R * q^(Hexp - 1) = Hexp * hloss / q

At the flow q* where hgradequals RQtolwe have:

    hloss* / q* = RQtol / Hexp

Using this point to define a linear head loss relation from 0 to q* yields RQtol / Hexp as the coefficient (or slope) for this relation.