JIANG/EPANET
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #537 from OpenWaterAnalytics/lrossman-hydcoeffs_refactor

Browse Source 
Refactor of hydcoeffs.c
This commit is contained in:
Lew Rossman
2019-10-08 10:17:52 -04:00
committed by GitHub
parent ad139a4f08 6e80fd8d99
commit 5601973d67
3 changed files with 66 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
ReleaseNotes2_2.md
View File

@@ -94,7 +94,7 @@ EPANET's original node re-ordering scheme has been replaced by the more efficien
EPANET's hydraulic solver can generate an ill-conditioned solution matrix when pipe flows approach zero unless some adjustment is made (i.e., as a pipe's flow approaches 0 its head loss gradient also approaches 0 causing its reciprocal, which is used to form the solution matrix's coefficients, to approach infinity). EPANET 2.0 used an arbitrary cutoff on head loss gradient to prevent it from becoming 0. This approach doesn't allow a pipe to follow any head loss v. flow relation in the region below the cutoff and can produce incorrect solutions for some networks (see [Estrada et al., 2009](https://ascelibrary.org/doi/full/10.1061/%28ASCE%29IR.1943-4774.0000100)). EPANET's hydraulic solver can generate an ill-conditioned solution matrix when pipe flows approach zero unless some adjustment is made (i.e., as a pipe's flow approaches 0 its head loss gradient also approaches 0 causing its reciprocal, which is used to form the solution matrix's coefficients, to approach infinity). EPANET 2.0 used an arbitrary cutoff on head loss gradient to prevent it from becoming 0. This approach doesn't allow a pipe to follow any head loss v. flow relation in the region below the cutoff and can produce incorrect solutions for some networks (see [Estrada et al., 2009](https://ascelibrary.org/doi/full/10.1061/%28ASCE%29IR.1943-4774.0000100)).
The hydraulic solver has been modified to use a linear head loss v. flow relation for flows approaching zero. For the Darcy-Weisbach equation, the linear Hagen-Poiseuille formula is used for laminar flow where the Reynolds Number is <= 2000. For the Hazen-Williams and Chezy-Manning equations, a flow limit is established for each pipe, equal to the flow that produces the EPANET 2 gradient cutoff. For flows below this a linear head loss relation is used whose gradient always equals the cutoff. EPANET 2.2 is now able to correctly solve the examples presented in Estrada et al. (2009) as well as those in [Gorev et al., (2013)](https://ascelibrary.org/doi/10.1061/%28ASCE%29HY.1943-7900.0000694) and [Elhay and Simpson (2011)](https://ascelibrary.org/doi/10.1061/%28ASCE%29HY.1943-7900.0000411). The hydraulic solver has been modified to use a linear head loss v. flow relation for flows approaching zero. For the Darcy-Weisbach equation, the linear Hagen-Poiseuille formula is used for laminar flow where the Reynolds Number is <= 2000. For the Hazen-Williams and Chezy-Manning equations, if the head loss gradient at a given flow is below the EPANET 2.0 gradient cutoff then a linear head loss relation is used whose slope equals the cutoff. EPANET 2.2 is now able to correctly solve the examples presented in Estrada et al. (2009) as well as those in [Gorev et al., (2013)](https://ascelibrary.org/doi/10.1061/%28ASCE%29HY.1943-7900.0000694) and [Elhay and Simpson (2011)](https://ascelibrary.org/doi/10.1061/%28ASCE%29HY.1943-7900.0000411).
## Pressure Dependent Demands ## Pressure Dependent Demands

124
src/hydcoeffs.c
View File

@@ -7,7 +7,7 @@
Authors: see AUTHORS Authors: see AUTHORS
Copyright: see AUTHORS Copyright: see AUTHORS
License: see LICENSE License: see LICENSE
Last Updated: 07/08/2019 Last Updated: 10/04/2019
****************************************************************************** ******************************************************************************
*/ */
@@ -79,7 +79,6 @@ void resistcoeff(Project *pr, int k)
double e, d, L; double e, d, L;
Slink *link = &net->Link[k]; Slink *link = &net->Link[k];
link->Qa = 0.0;
switch (link->Type) { switch (link->Type) {
// ... Link is a pipe. Compute resistance based on headloss formula. // ... Link is a pipe. Compute resistance based on headloss formula.
@@ -93,20 +92,14 @@ void resistcoeff(Project *pr, int k)
switch (hyd->Formflag) switch (hyd->Formflag)
{ {
case HW: case HW:
// ... resistance coeff.
link->R = 4.727 * L / pow(e, hyd->Hexp) / pow(d, 4.871); link->R = 4.727 * L / pow(e, hyd->Hexp) / pow(d, 4.871);
// ... flow below which linear head loss applies
link->Qa = pow(hyd->RQtol / hyd->Hexp / link->R, 1.17371);
break; break;
case DW: case DW:
link->R = L / 2.0 / 32.2 / d / SQR(PI * SQR(d) / 4.0); link->R = L / 2.0 / 32.2 / d / SQR(PI * SQR(d) / 4.0);
break; break;
case CM: case CM:
// ... resistance coeff.
link->R = SQR(4.0 * e / (1.49 * PI * SQR(d))) * link->R = SQR(4.0 * e / (1.49 * PI * SQR(d))) *
pow((d / 4.0), -1.333) * L; pow((d / 4.0), -1.333) * L;
// ... flow below which linear head loss applies
link->Qa = hyd->RQtol / 2.0 / link->R;
} }
break; break;
@@ -396,28 +389,23 @@ void emitterheadloss(Project *pr, int i, double *hloss, double *hgrad)
double ke; double ke;
double q; double q;
double qa;
// Set adjusted emitter coeff. // Set adjusted emitter coeff.
ke = MAX(CSMALL, pr->network.Node[i].Ke); ke = MAX(CSMALL, pr->network.Node[i].Ke);
// Find flow below which head loss is linear // Compute gradient of head loss through emitter
qa = pow(hyd->RQtol / ke / hyd->Qexp, 1.0 / (hyd->Qexp - 1.0));
// Use linear head loss relation for small flow
q = hyd->EmitterFlow[i]; q = hyd->EmitterFlow[i];
if (fabs(q) <= qa) *hgrad = hyd->Qexp * ke * pow(fabs(q), hyd->Qexp - 1.0);
// Use linear head loss function for small gradient
if (*hgrad < hyd->RQtol)
{ {
*hgrad = hyd->RQtol; *hgrad = hyd->RQtol;
*hloss = (*hgrad) * q; *hloss = (*hgrad) * q;
} }
// Otherwise use normal emitter function // Otherwise use normal emitter head loss function
else else *hloss = (*hgrad) * q / hyd->Qexp;
{
*hgrad = hyd->Qexp * ke * pow(fabs(q), hyd->Qexp - 1.0);
*hloss = (*hgrad) * q / hyd->Qexp;
}
} }
@@ -486,7 +474,6 @@ void demandheadloss(Project *pr, int i, double dp, double n,
{ {
Hydraul *hyd = &pr->hydraul; Hydraul *hyd = &pr->hydraul;
const double EPS = 0.01;
double d = hyd->DemandFlow[i]; double d = hyd->DemandFlow[i];
double dfull = hyd->NodeDemand[i]; double dfull = hyd->NodeDemand[i];
double r = d / dfull; double r = d / dfull;
@@ -498,18 +485,17 @@ void demandheadloss(Project *pr, int i, double dp, double n,
*hloss = CBIG * d; *hloss = CBIG * d;
} }
// Use linear head loss function for near zero demand
else if (r < EPS)
{
*hgrad = dp * pow(EPS, n) / dfull / EPS;
*hloss = (*hgrad) * d;
}
// Use power head loss function for demand less than full // Use power head loss function for demand less than full
else if (r < 1.0) else if (r < 1.0)
{ {
*hgrad = n * dp * pow(r, n - 1.0) / dfull; *hgrad = n * dp * pow(r, n - 1.0) / dfull;
*hloss = (*hgrad) * d / n; // ... use linear function for very small gradient
if (*hgrad < hyd->RQtol)
{
*hgrad = hyd->RQtol;
*hloss = (*hgrad) * d;
}
else *hloss = (*hgrad) * d / n;
} }
// Use upper barrier function for demand above full value // Use upper barrier function for demand above full value
@@ -560,19 +546,18 @@ void pipecoeff(Project *pr, int k)
ml = pr->network.Link[k].Km; ml = pr->network.Link[k].Km;
r = pr->network.Link[k].R; r = pr->network.Link[k].R;
// Friction head loss // Friction head loss gradient
// ... use linear relation for small flows hgrad = hyd->Hexp * r * pow(q, hyd->Hexp - 1.0);
if (q <= pr->network.Link[k].Qa)
// Friction head loss:
// ... use linear function for very small gradient
if (hgrad < hyd->RQtol)
{ {
hgrad = hyd->RQtol; hgrad = hyd->RQtol;
hloss = hgrad * q; hloss = hgrad * q;
} }
// ... use original formula for other flows // ... otherwise use original formula
else else hloss = hgrad * q / hyd->Hexp;
{
hgrad = hyd->Hexp * r * pow(q, hyd->Hexp - 1.0);
hloss = hgrad * q / hyd->Hexp;
}
// Contribution of minor head loss // Contribution of minor head loss
if (ml > 0.0) if (ml > 0.0)
@@ -702,7 +687,6 @@ void pumpcoeff(Project *pr, int k)
r, // Flow resistance coeff. r, // Flow resistance coeff.
n, // Flow exponent coeff. n, // Flow exponent coeff.
setting, // Pump speed setting setting, // Pump speed setting
qa, // Flow limit for linear head loss
hloss, // Head loss across pump hloss, // Head loss across pump
hgrad; // Head loss gradient hgrad; // Head loss gradient
Spump *pump; Spump *pump;
@@ -743,7 +727,7 @@ void pumpcoeff(Project *pr, int k)
pump->R = -r; pump->R = -r;
pump->N = 1.0; pump->N = 1.0;
// Compute head loss and its gradient // Compute head loss and its gradient (with speed adjustment)
hgrad = pump->R * setting ; hgrad = pump->R * setting ;
hloss = pump->H0 * SQR(setting) + hgrad * hyd->LinkFlow[k]; hloss = pump->H0 * SQR(setting) + hgrad * hyd->LinkFlow[k];
} }
@@ -755,24 +739,43 @@ void pumpcoeff(Project *pr, int k)
if (ABS(n - 1.0) < TINY) n = 1.0; if (ABS(n - 1.0) < TINY) n = 1.0;
r = pump->R * pow(setting, 2.0 - n); r = pump->R * pow(setting, 2.0 - n);
// Constant HP pump
if (pump->Ptype == CONST_HP)
{
// ... compute pump curve's gradient
hgrad = -r / q / q;
// ... use linear curve if gradient too large or too small
if (hgrad > CBIG)
{
hgrad = CBIG;
hloss = -hgrad * hyd->LinkFlow[k];
}
else if (hgrad < hyd->RQtol)
{
hgrad = hyd->RQtol;
hloss = -hgrad * hyd->LinkFlow[k];
}
// ... otherwise compute head loss from pump curve
else
{
hloss = r / hyd->LinkFlow[k];
}
}
// Compute head loss and its gradient // Compute head loss and its gradient
// ... pump curve is nonlinear // ... pump curve is nonlinear
if (n != 1.0) else if (n != 1.0)
{ {
// ... use linear approx. to pump curve for small flows // ... compute pump curve's gradient
if (pump->Ptype == CONST_HP) qa = -CBIG; hgrad = n * r * pow(q, n - 1.0);
else qa = pow(hyd->RQtol / n / r, 1.0 / (n - 1.0)); // ... use linear pump curve if gradient too small
if (q <= qa) if (hgrad < hyd->RQtol)
{ {
hgrad = hyd->RQtol; hgrad = hyd->RQtol;
hloss = h0 + hgrad * hyd->LinkFlow[k]; hloss = h0 + hgrad * hyd->LinkFlow[k];
} }
// ... use original pump curve for normal flows // ... otherwise compute head loss from pump curve
else else hloss = h0 + hgrad * hyd->LinkFlow[k] / n;
{
hgrad = n * r * pow(q, n - 1.0);
hloss = h0 + hgrad * hyd->LinkFlow[k] / n;
}
} }
// ... pump curve is linear // ... pump curve is linear
else else
@@ -1107,7 +1110,7 @@ void valvecoeff(Project *pr, int k)
Hydraul *hyd = &pr->hydraul; Hydraul *hyd = &pr->hydraul;
Slink *link = &pr->network.Link[k]; Slink *link = &pr->network.Link[k];
double flow, q, y, qa, hgrad; double flow, q, hloss, hgrad;
flow = hyd->LinkFlow[k]; flow = hyd->LinkFlow[k];
@@ -1122,23 +1125,20 @@ void valvecoeff(Project *pr, int k)
// Account for any minor headloss through the valve // Account for any minor headloss through the valve
if (link->Km > 0.0) if (link->Km > 0.0)
{ {
// Adjust for very small flow
q = fabs(flow); q = fabs(flow);
qa = hyd->RQtol / 2.0 / link->Km; hgrad = 2.0 * link->Km * q;
if (q <= qa)
// Guard against too small a head loss gradient
if (hgrad < hyd->RQtol)
{ {
hgrad = hyd->RQtol; hgrad = hyd->RQtol;
y = flow; hloss = flow * hgrad;
}
else
{
hgrad = 2.0 * link->Km * q;
y = flow / 2.0;
} }
else hloss = flow * hgrad / 2.0;
// P and Y coeffs. // P and Y coeffs.
hyd->P[k] = 1.0 / hgrad; hyd->P[k] = 1.0 / hgrad;
hyd->Y[k] = y; hyd->Y[k] = hloss / hgrad;
} }
// If no minor loss coeff. specified use a // If no minor loss coeff. specified use a

1
src/types.h
View File

@@ -394,7 +394,6 @@ typedef struct // Link Object
double Kw; // wall react. coef. double Kw; // wall react. coef.
double R; // flow resistance double R; // flow resistance
double Rc; // reaction coeff. double Rc; // reaction coeff.
double Qa; // low flow limit
LinkType Type; // link type LinkType Type; // link type
StatusType Status; // initial status StatusType Status; // initial status
int Rpt; // reporting flag int Rpt; // reporting flag