From 3ee30ce01951522646a4cdc1e6950ac5bccd4b05 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Fri, 7 Feb 2020 10:44:52 -0500 Subject: [PATCH] Replaced status checking for pumps & FCVs See ReleaseNotes2_3.md. --- ReleaseNotes2_3.md | 7 +++- src/hydcoeffs.c | 94 ++++++++++++++++++++-------------------------- src/hydsolver.c | 12 +++--- src/hydstatus.c | 14 +------ src/input2.c | 13 +++++-- 5 files changed, 64 insertions(+), 76 deletions(-) diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md index 1f922fc..f6b2f68 100644 --- a/ReleaseNotes2_3.md +++ b/ReleaseNotes2_3.md @@ -3,9 +3,14 @@ This document describes the changes and updates that have been made in version 2.3 of EPANET. - - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files to be opened by the toolkit. + - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files can be opened by the toolkit. - A `EN_setcurvetype` function was added to allow API clients to set a curve's type (e.g., `EN_PUMP_CURVE,` `EN_VOLUME_CURVE,` etc.). - A `EN_setvertex` function was added to allow API clients to change the coordinates of a link's vertex. - The index of a General Purpose Valve's (GPV's) head loss curve was added to the list of editable Link Properties using the symbolic constant name `EN_GPV_CURVE`. - The `EN_getlinkvalue` and `EN_setlinkvalue` functions were updated to get and set the value of `EN_GPV_CURVE`. + - For `EN_CUSTOM` type pump curves the maximum head value is now extrapolated to the y-axis intercept instead of being based on the first curve data point. Similarly, the maximum flow value is extrapolated to the x-axis intercept. + - Status checking for a pump not able to deliver enough head has been replaced by adding a penalty term to the pump's operating curve that prevents it from having negative flow (i.e., from crossing the y-axis). + - Status checking for Flow Control Valves has been eliminated by using a continuous head v. flow function. If the current flow is below the valve setting then the normal open head loss relation is used; otherwise a linear penalty function is applied to any flow in excess of the setting. Warnings are no longer issued when the valve operates fully opened at flows below the setting. + + diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index 9f4ac3d..b461946 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 10/04/2019 + Last Updated: 02/07/2020 ****************************************************************************** */ @@ -59,9 +59,9 @@ static void valvecoeff(Project *pr, int k); static void gpvcoeff(Project *pr, int k); static void pbvcoeff(Project *pr, int k); static void tcvcoeff(Project *pr, int k); +static void fcvcoeff(Project *pr, int k); static void prvcoeff(Project *pr, int k, int n1, int n2); static void psvcoeff(Project *pr, int k, int n1, int n2); -static void fcvcoeff(Project *pr, int k, int n1, int n2); void resistcoeff(Project *pr, int k) @@ -152,6 +152,8 @@ void headlosscoeffs(Project *pr) gpvcoeff(pr, k); break; case FCV: + fcvcoeff(pr, k); + break; case PRV: case PSV: if (hyd->LinkSetting[k] == MISSING) valvecoeff(pr, k); @@ -285,8 +287,8 @@ void valvecoeffs(Project *pr) ** Input: none ** Output: none ** Purpose: computes coeffs. of the linearized hydraulic eqns. -** contributed by PRVs, PSVs & FCVs whose status is -** not fixed to OPEN/CLOSED +** contributed by PRVs & PSVs whose status is not +** fixed to OPEN/CLOSED **-------------------------------------------------------------- */ { @@ -313,19 +315,8 @@ void valvecoeffs(Project *pr) n2 = link->N2; // Call valve-specific function - switch (link->Type) - { - case PRV: - prvcoeff(pr, k, n1, n2); - break; - case PSV: - psvcoeff(pr, k, n1, n2); - break; - case FCV: - fcvcoeff(pr, k, n1, n2); - break; - default: continue; - } + if (link->Type == PRV) prvcoeff(pr, k, n1, n2); + if (link->Type == PSV) psvcoeff(pr, k, n1, n2); } } @@ -701,20 +692,31 @@ void pumpcoeff(Project *pr, int k) } // Obtain reference to pump object - q = ABS(hyd->LinkFlow[k]); p = findpump(&pr->network, k); pump = &pr->network.Pump[p]; + + // Prevent negative flow + q = hyd->LinkFlow[k]; + if (q < 0.0) + { + hloss = -(SQR(setting) * pump->Hmax) + CBIG * q; + hgrad = CBIG; + hyd->P[k] = 1.0 / hgrad; + hyd->Y[k] = hloss / hgrad; + return; + } // If no pump curve treat pump as an open valve if (pump->Ptype == NOCURVE) { hyd->P[k] = 1.0 / CSMALL; - hyd->Y[k] = hyd->LinkFlow[k]; + hyd->Y[k] = q; return; } // Get pump curve coefficients for custom pump curve // (Other pump types have pre-determined coeffs.) + q = ABS(q); if (pump->Ptype == CUSTOM) { // Find intercept (h0) & slope (r) of pump curve @@ -1044,12 +1046,10 @@ void psvcoeff(Project *pr, int k, int n1, int n2) } -void fcvcoeff(Project *pr, int k, int n1, int n2) +void fcvcoeff(Project *pr, int k) /* **-------------------------------------------------------------- ** Input: k = link index -** n1 = upstream node of valve -** n2 = downstream node of valve ** Output: none ** Purpose: computes solution matrix coeffs. for flow control ** valve @@ -1059,40 +1059,28 @@ void fcvcoeff(Project *pr, int k, int n1, int n2) Hydraul *hyd = &pr->hydraul; Smatrix *sm = &hyd->smatrix; - int i, j; // Rows in solution matrix - double q; // Valve flow setting + double qset; // Valve flow setting + double flow; // Current valve flow + double hloss, hgrad; // Head loss & gradient - q = hyd->LinkSetting[k]; - i = sm->Row[n1]; - j = sm->Row[n2]; - - // If valve active, break network at valve and treat - // flow setting as external demand at upstream node - // and external supply at downstream node. - - if (hyd->LinkStatus[k] == ACTIVE) - { - hyd->Xflow[n1] -= q; - hyd->Xflow[n2] += q; - hyd->Y[k] = hyd->LinkFlow[k] - q; - sm->F[i] -= q; - sm->F[j] += q; - hyd->P[k] = 1.0 / CBIG; - sm->Aij[sm->Ndx[k]] -= hyd->P[k]; - sm->Aii[i] += hyd->P[k]; - sm->Aii[j] += hyd->P[k]; - } - - // Otherwise treat valve as an open pipe - - else + // Treat as a regular valve if status fixed or flow below setting + qset = hyd->LinkSetting[k]; + flow = hyd->LinkFlow[k]; + if (qset == MISSING || hyd->LinkStatus[k] <= CLOSED || flow < qset) { valvecoeff(pr, k); - sm->Aij[sm->Ndx[k]] -= hyd->P[k]; - sm->Aii[i] += hyd->P[k]; - sm->Aii[j] += hyd->P[k]; - sm->F[i] += (hyd->Y[k] - hyd->LinkFlow[k]); - sm->F[j] -= (hyd->Y[k] - hyd->LinkFlow[k]); + } + + // Otherwise prevent flow from exceeding the setting + else + { + hyd->LinkFlow[k] = qset; + valvecoeff(pr, k); + hloss = hyd->Y[k] / hyd->P[k] + CBIG * (flow - qset); + hgrad = CBIG; + hyd->P[k] = 1.0 / hgrad; + hyd->Y[k] = hloss / hgrad; + hyd->LinkFlow[k] = flow; } } diff --git a/src/hydsolver.c b/src/hydsolver.c index fcb8258..53057c2 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -8,7 +8,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 07/15/2019 + Last Updated: 02/07/2020 ****************************************************************************** */ @@ -111,6 +111,7 @@ int hydsolve(Project *pr, int *iter, double *relerr) maxtrials = hyd->MaxIter; if (hyd->ExtraIter > 0) maxtrials += hyd->ExtraIter; *iter = 1; + headlosscoeffs(pr); while (*iter <= maxtrials) { // Compute coefficient matrices A & F and solve A*H = F @@ -118,7 +119,6 @@ int hydsolve(Project *pr, int *iter, double *relerr) // head loss gradients, & F = flow correction terms. // Solution for H is returned in F from call to linsolve(). - headlosscoeffs(pr); matrixcoeffs(pr); errcode = linsolve(sm, net->Njuncs); @@ -138,6 +138,9 @@ int hydsolve(Project *pr, int *iter, double *relerr) } newerr = newflows(pr, &hydbal); // Update flows *relerr = newerr; + + // Compute head loss coeffs. for new flows + headlosscoeffs(pr); // Write convergence error to status report if called for if (rpt->Statflag == FULL) @@ -243,7 +246,7 @@ int badvalve(Project *pr, int n) if (n == n1 || n == n2) { t = link->Type; - if (t == PRV || t == PSV || t == FCV) + if (t == PRV || t == PSV) { if (hyd->LinkStatus[k] == ACTIVE) { @@ -253,8 +256,7 @@ int badvalve(Project *pr, int n) clocktime(rpt->Atime, time->Htime), link->ID); writeline(pr, pr->Msg); } - if (link->Type == FCV) hyd->LinkStatus[k] = XFCV; - else hyd->LinkStatus[k] = XPRESSURE; + hyd->LinkStatus[k] = XPRESSURE; return 1; } } diff --git a/src/hydstatus.c b/src/hydstatus.c index 87d45c6..76ca823 100644 --- a/src/hydstatus.c +++ b/src/hydstatus.c @@ -7,7 +7,7 @@ Description: updates hydraulic status of network elements Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 05/15/2019 +Last Updated: 02/07/2020 ****************************************************************************** */ @@ -141,18 +141,6 @@ int linkstatus(Project *pr) hyd->LinkStatus[k] = cvstatus(pr, hyd->LinkStatus[k], dh, hyd->LinkFlow[k]); } - if (link->Type == PUMP && hyd->LinkStatus[k] >= OPEN && - hyd->LinkSetting[k] > 0.0) - { - hyd->LinkStatus[k] = pumpstatus(pr, k, -dh); - } - - // Check for status changes in non-fixed FCVs - if (link->Type == FCV && hyd->LinkSetting[k] != MISSING) - { - hyd->LinkStatus[k] = fcvstatus(pr, k, status, hyd->NodeHead[n1], - hyd->NodeHead[n2]); - } // Check for flow into (out of) full (empty) tanks if (n1 > net->Njuncs || n2 > net->Njuncs) diff --git a/src/input2.c b/src/input2.c index a4d7fc8..136d297 100644 --- a/src/input2.c +++ b/src/input2.c @@ -7,7 +7,7 @@ Description: reads and interprets network data from an EPANET input file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 02/03/2020 +Last Updated: 02/07/2020 ****************************************************************************** */ @@ -429,9 +429,14 @@ int updatepumpparams(Project *pr, int pumpindex) { if (curve->Y[m] >= curve->Y[m - 1]) return 227; } - pump->Qmax = curve->X[npts - 1]; - pump->Q0 = (curve->X[0] + pump->Qmax) / 2.0; - pump->Hmax = curve->Y[0]; + pump->Q0 = (curve->X[0] + curve->X[npts-1]) / 2.0; + + // Extend curve to find Hmax (at 0 flow) and Qmax (at 0 head) + b = (curve->Y[1] - curve->Y[0]) / (curve->X[1] - curve->X[0]); + pump->Hmax = curve->Y[0] + b * curve->X[0]; + b = (curve->Y[npts-1] - curve->Y[npts-2]) / + (curve->X[npts-1] - curve->X[npts-2]); + pump->Qmax = curve->X[npts-1] - curve->Y[npts-1] / b; } // Compute shape factors & limits of power function curves