From 6b808c03ab23f042777352de03eface41cda9cae Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Tue, 11 Sep 2018 09:48:01 -0400 Subject: [PATCH 1/6] Fixes zero flow issue (#164) and DW eqn. (#199) --- src/funcs.h | 11 +- src/hydcoeffs.c | 457 ++++++++++++++++++++++++++---------------------- src/hydsolver.c | 22 ++- 3 files changed, 269 insertions(+), 221 deletions(-) diff --git a/src/funcs.h b/src/funcs.h index b961a82b..52b2d70c 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -178,11 +178,12 @@ int hydsolve(EN_Project *pr, int *,double *); /* Solves network equations void resistcoeff(EN_Project *pr, int k); /* Finds pipe flow resistance */ void headlosscoeffs(EN_Project *pr); // Finds link head loss coeffs. void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. */ -double emitflowchange(EN_Project *pr, int i); /* Change in emitter outflow */ -double demandflowchange(EN_Project *pr, int i, // Change in demand outflow - double dp, double n); -void demandparams(EN_Project *pr, double *dp, // PDA function parameters - double *n); +void emitheadloss(EN_Project *pr, int, // Finds emitter head loss + double *, double *); +double demandflowchange(EN_Project *pr, int, // Change in demand outflow + double, double); +void demandparams(EN_Project *pr, double *, // PDA function parameters + double *); /* ----------- SMATRIX.C ---------------*/ int createsparse(EN_Project *pr); /* Creates sparse matrix */ diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index 00b0b820..84b67c87 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -18,21 +18,24 @@ HYDCOEFFS.C -- hydraulic coefficients for the EPANET Program #include "funcs.h" // Constants used for computing Darcy-Weisbach friction factor -const double A1 = 0.314159265359e04; // 1000*PI -const double A2 = 0.157079632679e04; // 500*PI -const double A3 = 0.502654824574e02; // 16*PI -const double A4 = 6.283185307; // 2*PI -const double A8 = 4.61841319859; // 5.74*(PI/4)^.9 -const double A9 = -8.685889638e-01; // -2/ln(10) -const double AA = -1.5634601348; // -2*.9*2/ln(10) -const double AB = 3.28895476345e-03; // 5.74/(4000^.9) -const double AC = -5.14214965799e-03; // AA*AB +const double A1 = 3.14159265358979323850e+03; // 1000*PI +const double A2 = 1.57079632679489661930e+03; // 500*PI +const double A3 = 5.02654824574366918160e+01; // 16*PI +const double A4 = 6.28318530717958647700e+00; // 2*PI +const double A8 = 4.61841319859066668690e+00; // 5.74*(PI/4)^.9 +const double A9 = -8.68588963806503655300e-01; // -2/ln(10) +const double AA = -1.5634601348517065795e+00; // -2*.9*2/ln(10) +const double AB = 3.28895476345399058690e-03; // 5.74/(4000^.9) +const double AC = -5.14214965799093883760e-03; // AA*AB + +// Cutoff flow for using linear head loss relation +const double Q_CUTOFF = 1.0e-5; // External functions //void resistcoeff(EN_Project *pr, int k); //void headlosscoeffs(EN_Project *pr); //void matrixcoeffs(EN_Project *pr); -//double emitflowchange(EN_Project *pr, int i); +//void emitheadloss(EN_Project *pr, int i, double *hloss, double *dhdq); //double demandflowchange(EN_Project *pr, int i, double dp, double n); //void demandparams(EN_Project *pr, double *dp, double *n); @@ -47,7 +50,7 @@ static void demandheadloss(double d, double dfull, double dp, static void pipecoeff(EN_Project *pr, int k); static void DWpipecoeff(EN_Project *pr, int k); -static double frictionFactor(EN_Project *pr, int k, double *dfdq); +static double frictionFactor(double q, double e, double s, double *dfdq); static void pumpcoeff(EN_Project *pr, int k); static void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r); @@ -195,7 +198,8 @@ void linkcoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes matrix coefficients for links +** Purpose: computes coefficients contributed by links to the +** linearized system of hydraulic equations. **-------------------------------------------------------------- */ { @@ -214,32 +218,34 @@ void linkcoeffs(EN_Project *pr) n1 = link->N1; // Start node of link n2 = link->N2; // End node of link - // Update net nodal inflows (X_tmp), solution matrix (A) and RHS array (F) - // (Use covention that flow out of node is (-), flow into node is (+)) + // Update nodal flow balance (X_tmp) + // (Flow out of node is (-), flow into node is (+)) hyd->X_tmp[n1] -= hyd->LinkFlows[k]; hyd->X_tmp[n2] += hyd->LinkFlows[k]; - // Off-diagonal coeff. + // Add to off-diagonal coeff. of linear system matrix sol->Aij[sol->Ndx[k]] -= sol->P[k]; - // Node n1 is junction + // Update linear system coeffs. associated with start node n1 + // ... node n1 is junction if (n1 <= net->Njuncs) { sol->Aii[sol->Row[n1]] += sol->P[k]; // Diagonal coeff. sol->F[sol->Row[n1]] += sol->Y[k]; // RHS coeff. } - // Node n1 is a tank/reservoir + // ... node n1 is a tank/reservoir else sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]); - // Node n2 is junction + // Update linear system coeffs. associated with end node n2 + // ... node n2 is junction if (n2 <= net->Njuncs) { sol->Aii[sol->Row[n2]] += sol->P[k]; // Diagonal coeff. sol->F[sol->Row[n2]] -= sol->Y[k]; // RHS coeff. } - // Node n2 is a tank/reservoir + // ... node n2 is a tank/reservoir else sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]); } } @@ -250,8 +256,8 @@ void nodecoeffs(EN_Project *pr) **---------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: completes calculation of nodal flow imbalance (X_tmp) -** & flow correction (F) arrays +** Purpose: completes calculation of nodal flow balance array +** (X_tmp) & r.h.s. (F) of linearized hydraulic eqns. **---------------------------------------------------------------- */ { @@ -261,7 +267,7 @@ void nodecoeffs(EN_Project *pr) EN_Network *net = &pr->network; // For junction nodes, subtract demand flow from net - // flow imbalance & add imbalance to RHS array F. + // flow balance & add flow balance to RHS array F for (i = 1; i <= net->Njuncs; i++) { hyd->X_tmp[i] -= hyd->DemandFlows[i]; @@ -275,8 +281,9 @@ void valvecoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs -** whose status is not fixed to OPEN/CLOSED +** Purpose: computes coeffs. of the linearized hydraulic eqns. +** contributed by PRVs, PSVs & FCVs whose status is +** not fixed to OPEN/CLOSED **-------------------------------------------------------------- */ { @@ -325,7 +332,8 @@ void emittercoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes matrix coeffs. for emitters +** Purpose: computes coeffs. of the linearized hydraulic eqns. +** contributed by emitters. ** ** Note: Emitters consist of a fictitious pipe connected to ** a fictitious reservoir whose elevation equals that @@ -334,12 +342,8 @@ void emittercoeffs(EN_Project *pr) **-------------------------------------------------------------- */ { - int i; - double ke; - double p; - double q; - double y; - double z; + int i, row; + double hloss, hgrad; hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; @@ -348,46 +352,57 @@ void emittercoeffs(EN_Project *pr) for (i = 1; i <= net->Njuncs; i++) { + // Skip junctions without emitters node = &net->Node[i]; if (node->Ke == 0.0) continue; - ke = MAX(CSMALL, node->Ke); // emitter coeff. - q = hyd->EmitterFlows[i]; // emitter flow - z = ke * pow(ABS(q), hyd->Qexp); // emitter head loss - p = hyd->Qexp * z / ABS(q); // head loss gradient - if (p < hyd->RQtol) - { - p = 1.0 / hyd->RQtol; - } - else - { - p = 1.0 / p; // inverse head loss gradient - } - y = SGN(q)*z*p; // head loss / gradient - sol->Aii[sol->Row[i]] += p; // addition to main diagonal - sol->F[sol->Row[i]] += y + p * node->El; // addition to r.h.s. - hyd->X_tmp[i] -= q; // addition to net node inflow + + // Find emitter head loss and gradient + emitheadloss(pr, i, &hloss, &hgrad); + + // Row of solution matrix + row = sol->Row[i]; + + // Addition to matrix diagonal & r.h.s + sol->Aii[row] += 1.0 / hgrad; + sol->F[row] += (hloss + node->El) / hgrad; + + // Update to node flow balance + hyd->X_tmp[i] -= hyd->EmitterFlows[i]; } } -double emitflowchange(EN_Project *pr, int i) +void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad) /* -**-------------------------------------------------------------- +**------------------------------------------------------------- ** Input: i = node index -** Output: returns change in flow at an emitter node -** Purpose: computes flow change at an emitter node -**-------------------------------------------------------------- +** Output: hloss = head loss across node's emitter +** hgrad = head loss gradient +** Purpose: computes an emitters's head loss and gradient. +**------------------------------------------------------------- */ { - double ke, p; + double ke; + double q; hydraulics_t *hyd = &pr->hydraulics; - Snode *node = &pr->network.Node[i]; - ke = MAX(CSMALL, node->Ke); - p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]), (hyd->Qexp - 1.0)); - if (p < hyd->RQtol) p = 1 / hyd->RQtol; - else p = 1.0 / p; - return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El)); + // Set adjusted emitter coeff. + ke = MAX(CSMALL, pr->network.Node[i].Ke); + + // Use linear head loss relation for small flow + q = hyd->EmitterFlows[i]; + if (fabs(q) <= Q_CUTOFF) + { + *hgrad = hyd->Qexp * ke * pow(Q_CUTOFF, hyd->Qexp - 1.0); + *hloss = (*hgrad) * q; + } + + // Otherwise use normal emitter function + else + { + *hgrad = hyd->Qexp * ke * pow(fabs(q), hyd->Qexp - 1.0); + *hloss = (*hgrad) * q / hyd->Qexp; + } } @@ -397,8 +412,8 @@ void demandparams(EN_Project *pr, double *dp, double *n) ** Input: none ** Output: dp = pressure range over which demands can vary ** n = exponent in head loss v. demand function -** Purpose: retrieves parameters that define a pressure dependent -** demand function. +** Purpose: retrieves parameters that define a pressure +** dependent demand function. **-------------------------------------------------------------- */ { @@ -427,12 +442,13 @@ void demandcoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes matrix coeffs. for pressure dependent demands +** Purpose: computes coeffs. of the linearized hydraulic eqns. +** contributed by pressure dependent demands. ** ** Note: Pressure dependent demands are modelled like emitters -** with Hloss = Pserv * (D / Dfull)^(1/Pexp) +** with Hloss = Preq * (D / Dfull)^(1/Pexp) ** where D (actual demand) is zero for negative pressure -** and is Dfull above pressure Pserv. +** and is Dfull above pressure Preq. **-------------------------------------------------------------- */ { @@ -542,25 +558,23 @@ void pipecoeff(EN_Project *pr, int k) **-------------------------------------------------------------- ** Input: k = link index ** Output: none -** Purpose: computes P & Y coefficients for pipe k +** Purpose: computes P & Y coefficients for pipe k. ** -** P = inverse head loss gradient = 1/(dh/dQ) -** Y = flow correction term = h*P +** P = inverse head loss gradient = 1/hgrad +** Y = flow correction term = hloss / hgrad **-------------------------------------------------------------- */ { - double hpipe, // Normal head loss - hml, // Minor head loss + double hloss, // Head loss + hgrad, // Head loss gradient ml, // Minor loss coeff. - p, // q*(dh/dq) q, // Abs. value of flow r; // Resistance coeff. hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; - Slink *link = &pr->network.Link[k]; - // For closed pipe use headloss formula: h = CBIG*q + // For closed pipe use headloss formula: hloss = CBIG*q if (hyd->LinkStatus[k] <= CLOSED) { sol->P[k] = 1.0 / CBIG; @@ -575,31 +589,35 @@ void pipecoeff(EN_Project *pr, int k) return; } - // Evaluate headloss coefficients - q = ABS(hyd->LinkFlows[k]); // Absolute flow - ml = link->Km; // Minor loss coeff. - r = link->R; // Resistance coeff. + q = ABS(hyd->LinkFlows[k]); + ml = pr->network.Link[k].Km; + r = pr->network.Link[k].R; - // Use large P coefficient for small flow resistance product - if ( (r+ml)*q < hyd->RQtol) + // Friction head loss + if (q <= Q_CUTOFF) { - sol->P[k] = 1.0 / hyd->RQtol; - sol->Y[k] = hyd->LinkFlows[k] / hyd->Hexp; - return; + hgrad = hyd->Hexp * r * pow(Q_CUTOFF, hyd->Hexp-1.0); + hloss = q * hgrad; + } + else + { + hgrad = hyd->Hexp * r * pow(q, hyd->Hexp - 1.0); + hloss = hgrad * q / hyd->Hexp; } - // Compute P and Y coefficients - hpipe = r*pow(q, hyd->Hexp); // Friction head loss - p = hyd->Hexp*hpipe; // Q*dh(friction)/dQ + // Contribution of minor head loss if (ml > 0.0) { - hml = ml*q*q; // Minor head loss - p += 2.0*hml; // Q*dh(Total)/dQ + hloss += ml * q * q; + hgrad += 2.0 * ml * q; } - else hml = 0.0; - p = hyd->LinkFlows[k] / p; // 1 / (dh/dQ) - sol->P[k] = ABS(p); - sol->Y[k] = p*(hpipe + hml); + + // Adjust head loss sign for flow direction + hloss *= SGN(hyd->LinkFlows[k]); + + // P and Y coeffs. + sol->P[k] = 1.0 / hgrad; + sol->Y[k] = hloss / hgrad; } @@ -618,102 +636,85 @@ void DWpipecoeff(EN_Project *pr, int k) Slink *link = &pr->network.Link[k]; double q = ABS(hyd->LinkFlows[k]); - double dfdq = 0.0; - double r, r1, f, ml, p, hloss; - - ml = link->Km; // Minor loss coeff. - r = link->R; // Resistance coeff. - f = frictionFactor(pr,k,&dfdq); // D-W friction factor - r1 = f*r+ml; - - // Use large P coefficient for small flow resistance product - if (r1*q < hyd->RQtol) + double r = link->R; // Resistance coeff. + double ml = link->Km; // Minor loss coeff. + double e = link->Kc / link->Diam; // Relative roughness + double s = hyd->Viscos * link->Diam; // Viscosity / diameter + + double hloss, hgrad, f, dfdq, r1; + + // Compute head loss and its derivative + // ... use Hagen-Poiseuille formula for laminar flow (Re <= 2000) + if (q <= A2 * s) { - sol->P[k] = 1.0/hyd->RQtol; - sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp; - return; + r = 16.0 * PI * s * r; + hloss = hyd->LinkFlows[k] * (r + ml * q); + hgrad = r + 2.0 * ml * q; } - + + // ... otherwise use Darcy-Weisbach formula with friction factor + else + { + dfdq = 0.0; + f = frictionFactor(q, e, s, &dfdq); + r1 = f * r + ml; + hloss = r1 * q * hyd->LinkFlows[k]; + hgrad = (2.0 * r1 * q) + (dfdq * r * q * q); + } + // Compute P and Y coefficients - hloss = r1*SQR(q); // Total head loss - p = 2.0*r1*q; // |dHloss/dQ| - // + dfdq*r*q*q; // Ignore df/dQ term - p = 1.0 / p; - sol->P[k] = p; - sol->Y[k] = SGN(hyd->LinkFlows[k]) * hloss * p; + sol->P[k] = 1.0 / hgrad; + sol->Y[k] = hloss / hgrad; } -double frictionFactor(EN_Project *pr, int k, double *dfdq) +double frictionFactor(double q, double e, double s, double *dfdq) /* **-------------------------------------------------------------- -** Input: k = link index -** Output: returns friction factor and -** replaces dfdq (derivative of f w.r.t. flow) +** Input: q = |pipe flow| +** e = pipe roughness / diameter +** s = viscosity * pipe diameter +** Output: dfdq = derivative of friction factor w.r.t. flow +** Returns: pipe's friction factor ** Purpose: computes Darcy-Weisbach friction factor and its ** derivative as a function of Reynolds Number (Re). -** -** Note: Current formulas for dfdq need to be corrected -** so dfdq returned as 0. **-------------------------------------------------------------- */ { - double q, // Abs. value of flow - f; // Friction factor - double x1,x2,x3,x4, - y1,y2,y3, - fa,fb,r; - double s,w; - - hydraulics_t *hyd = &pr->hydraulics; - Slink *link = &pr->network.Link[k]; - - *dfdq = 0.0; - if (link->Type > EN_PIPE) - return(1.0); // Only apply to pipes - q = ABS(hyd->LinkFlows[k]); - s = hyd->Viscos * link->Diam; - w = q/s; // w = Re(Pi/4) - - // For Re >= 4000 use Colebrook Formula - if (w >= A1) - { - y1 = A8/pow(w,0.9); - y2 = link->Kc/(3.7*link->Diam) + y1; - y3 = A9*log(y2); - f = 1.0/SQR(y3); - /* *dfdq = (2.0+AA*y1/(y2*y3))*f; */ /* df/dq */ - } - - // For Re > 2000 use Interpolation Formula - else if (w > A2) - { - y2 = link->Kc/(3.7*link->Diam) + AB; - y3 = A9*log(y2); - fa = 1.0/SQR(y3); - fb = (2.0+AC/(y2*y3))*fa; - r = w/A2; - x1 = 7.0*fa - fb; - x2 = 0.128 - 17.0*fa + 2.5*fb; - x3 = -0.128 + 13.0*fa - (fb+fb); - x4 = r*(0.032 - 3.0*fa + 0.5*fb); - f = x1 + r*(x2 + r*(x3+x4)); - /* *dfdq = (x1 + x1 + r*(3.0*x2 + r*(4.0*x3 + 5.0*x4))); */ - } - - // For Re > 8 (Laminar Flow) use Hagen-Poiseuille Formula - else if (w > A4) - { - f = A3*s/q; // 16 * PI * Viscos * Diam / Flow - /* *dfdq = A3*s; */ - } - else - { - f = 8.0; - *dfdq = 0.0; - } - return(f); + double f; // friction factor + double x1, x2, x3, x4, + y1, y2, y3, + fa, fb, r; + double w = q / s; // Re*Pi/4 + + // For Re >= 4000 use Swamee & Jain approximation + // of the Colebrook-White Formula + if ( w >= A1 ) + { + y1 = A8 / pow(w, 0.9); + y2 = e / 3.7 + y1; + y3 = A9 * log(y2); + f = 1.0 / (y3*y3); + *dfdq = 1.8 * f * y1 * A9 / y2 / y3 / q; + } + // Use interpolating polynomials developed by + // E. Dunlop for transition flow from 2000 < Re < 4000. + else + { + y2 = e / 3.7 + AB; + y3 = A9 * log(y2); + fa = 1.0 / (y3*y3); + fb = (2.0 + AC / (y2*y3)) * fa; + r = w / A2; + x1 = 7.0 * fa - fb; + x2 = 0.128 - 17.0 * fa + 2.5 * fb; + x3 = -0.128 + 13.0 * fa - (fb + fb); + x4 = 0.032 - 3.0 * fa + 0.5 *fb; + f = x1 + r * (x2 + r * (x3 + r * x4)); + *dfdq = (x2 + r * (2.0 * x3 + r * 3.0 * x4)) / s / A2; + } + return f; } @@ -730,13 +731,17 @@ void pumpcoeff(EN_Project *pr, int k) double h0, // Shutoff head q, // Abs. value of flow r, // Flow resistance coeff. - n; // Flow exponent coeff. + n, // Flow exponent coeff. + setting, // Pump speed setting + hloss, // Head loss across pump + hgrad; // Head loss gradient + hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; - double setting = hyd->LinkSetting[k]; - Spump *pump; + Spump *pump; // Use high resistance pipe if pump closed or cannot deliver head + setting = hyd->LinkSetting[k]; if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) { sol->P[k] = 1.0 / CBIG; @@ -744,35 +749,54 @@ void pumpcoeff(EN_Project *pr, int k) return; } + // Obtain reference to pump object & its speed setting q = ABS(hyd->LinkFlows[k]); - q = MAX(q, TINY); - - // Obtain reference to pump object p = findpump(&pr->network, k); pump = &pr->network.Pump[p]; - // Get pump curve coefficients for custom pump curve. + // Get pump curve coefficients for custom pump curve + // (Other pump types have pre-determined coeffs.) if (pump->Ptype == CUSTOM) { // Find intercept (h0) & slope (r) of pump curve // line segment which contains speed-adjusted flow. curvecoeff(pr, pump->Hcurve, q / setting, &h0, &r); - // Determine head loss coefficients. + // Determine head loss coefficients (negative sign + // converts from pump curve's head gain to head loss) pump->H0 = -h0; pump->R = -r; pump->N = 1.0; - } - // Adjust head loss coefficients for pump speed. - h0 = SQR(setting) * pump->H0; - n = pump->N; - r = pump->R * pow(setting, 2.0 - n); - if (n != 1.0) r = n * r * pow(q, n - 1.0); + // Compute head loss and its gradient + hgrad = pump->R * setting ; + hloss = pump->H0 * SQR(setting) + hgrad * hyd->LinkFlows[k]; + } + else + { + // Adjust head loss coefficients for pump speed + h0 = SQR(setting) * pump->H0; + n = pump->N; + r = pump->R * pow(setting, 2.0 - n); + + // Compute head loss and its gradient + // ... use linear approx. to pump curve for small flows + if (q <= Q_CUTOFF) + { + hgrad = r * pow(Q_CUTOFF, n - 1.0); + hloss = h0 + hgrad * hyd->LinkFlows[k]; + } + // ... use original pump curve for normal flows + else + { + hgrad = n * r * pow(q, n - 1.0); + hloss = h0 + hgrad * hyd->LinkFlows[k] / n; + } + } - // Compute inverse headloss gradient (P) and flow correction factor (Y) - sol->P[k] = 1.0 / MAX(r, hyd->RQtol); - sol->Y[k] = hyd->LinkFlows[k] / n + sol->P[k] * h0; + // P and Y coeffs. + sol->P[k] = 1.0 / hgrad; + sol->Y[k] = hloss / hgrad; } @@ -825,9 +849,10 @@ void gpvcoeff(EN_Project *pr, int k) **-------------------------------------------------------------- */ { - double h0, // Headloss curve intercept - q, // Abs. value of flow - r; // Flow resistance coeff. + int i; + double h0, // Intercept of head loss curve segment + r, // Slope of head loss curve segment + q; // Abs. value of flow hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; @@ -835,19 +860,25 @@ void gpvcoeff(EN_Project *pr, int k) // Treat as a pipe if valve closed if (hyd->LinkStatus[k] == CLOSED) valvecoeff(pr, k); - // Otherwise utilize headloss curve - // whose index is stored in K + // Otherwise utilize segment of head loss curve + // bracketing current flow (curve index is stored + // in valve's setting) else { - // Find slope & intercept of headloss curve. + // Index of valve's head loss curve + i = (int)ROUND(hyd->LinkSetting[k]); + + // Adjusted flow rate q = ABS(hyd->LinkFlows[k]); q = MAX(q, TINY); - curvecoeff(pr, (int)ROUND(hyd->LinkSetting[k]), q, &h0, &r); - // Compute inverse headloss gradient (P) - // and flow correction factor (Y). - sol->P[k] = 1.0 / MAX(r, hyd->RQtol); - sol->Y[k] = sol->P[k] * (h0 + r*q) * SGN(hyd->LinkFlows[k]); + // Intercept and slope of curve segment containing q + curvecoeff(pr, i, q, &h0, &r); + r = MAX(r, TINY); + + // Resulting P and Y coeffs. + sol->P[k] = 1.0 / r; + sol->Y[k] = (h0 / r + q) * SGN(hyd->LinkFlows[k]); } } @@ -1084,14 +1115,14 @@ void valvecoeff(EN_Project *pr, int k) **-------------------------------------------------------------- */ { - double p; + double flow, q, dhdq; EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; Slink *link = &net->Link[k]; - double flow = hyd->LinkFlows[k]; + flow = hyd->LinkFlows[k]; // Valve is closed. Use a very small matrix coeff. if (hyd->LinkStatus[k] <= CLOSED) @@ -1104,14 +1135,26 @@ void valvecoeff(EN_Project *pr, int k) // Account for any minor headloss through the valve if (link->Km > 0.0) { - p = 2.0 * link->Km * fabs(flow); - if (p < hyd->RQtol) p = hyd->RQtol; - sol->P[k] = 1.0 / p; - sol->Y[k] = flow / 2.0; + // Adjust for very small flow + q = MAX(fabs(flow), Q_CUTOFF); + + // Gradient of minor head loss formula + dhdq = 2.0 * link->Km * q; + sol->P[k] = 1.0 / dhdq; + + // Y = head loss / gradient + // Use linear head loss relation for low flow + if (q == Q_CUTOFF) sol->Y[k] = flow; + + // Otherwise use normal minor loss formula + else sol->Y[k] = flow / 2.0; } + + // If no minor loss coeff. specified use a + // low resistance linear head loss relation else { - sol->P[k] = 1.0 / hyd->RQtol; + sol->P[k] = 1.0 / CSMALL; sol->Y[k] = flow; } } diff --git a/src/hydsolver.c b/src/hydsolver.c index 7e6f0d65..63ed951c 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -494,30 +494,34 @@ void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum, **---------------------------------------------------------------- */ { - double dq; - int k; + int i; + double hloss, hgrad, dh, dq; EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; // Examine each network junction - for (k = 1; k <= net->Njuncs; k++) + for (i = 1; i <= net->Njuncs; i++) { // Skip junction if it does not have an emitter - if (net->Node[k].Ke == 0.0) continue; + if (net->Node[i].Ke == 0.0) continue; - // Find emitter flow change (see hydcoeffs.c) - dq = emitflowchange(pr, k); - hyd->EmitterFlows[k] -= dq; + // Find emitter head loss and gradient + emitheadloss(pr, i, &hloss, &hgrad); + + // Find emitter flow change + dh = hyd->NodeHead[i] - net->Node[i].El; + dq = (hloss - dh) / hgrad; + hyd->EmitterFlows[i] -= dq; // Update system flow summation - *qsum += ABS(hyd->EmitterFlows[k]); + *qsum += ABS(hyd->EmitterFlows[i]); *dqsum += ABS(dq); // Update identity of element with max. flow change if (ABS(dq) > hbal->maxflowchange) { hbal->maxflowchange = ABS(dq); - hbal->maxflownode = k; + hbal->maxflownode = i; hbal->maxflowlink = -1; } } From 836d9c36682b7345d95e110ddfb0770a36a43c11 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Tue, 11 Sep 2018 11:04:47 -0400 Subject: [PATCH 2/6] Revert "Fixes zero flow issue (#164) and DW eqn. (#199)" This reverts commit 6b808c03ab23f042777352de03eface41cda9cae. --- src/funcs.h | 11 +- src/hydcoeffs.c | 457 ++++++++++++++++++++++-------------------------- src/hydsolver.c | 22 +-- 3 files changed, 221 insertions(+), 269 deletions(-) diff --git a/src/funcs.h b/src/funcs.h index 52b2d70c..b961a82b 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -178,12 +178,11 @@ int hydsolve(EN_Project *pr, int *,double *); /* Solves network equations void resistcoeff(EN_Project *pr, int k); /* Finds pipe flow resistance */ void headlosscoeffs(EN_Project *pr); // Finds link head loss coeffs. void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. */ -void emitheadloss(EN_Project *pr, int, // Finds emitter head loss - double *, double *); -double demandflowchange(EN_Project *pr, int, // Change in demand outflow - double, double); -void demandparams(EN_Project *pr, double *, // PDA function parameters - double *); +double emitflowchange(EN_Project *pr, int i); /* Change in emitter outflow */ +double demandflowchange(EN_Project *pr, int i, // Change in demand outflow + double dp, double n); +void demandparams(EN_Project *pr, double *dp, // PDA function parameters + double *n); /* ----------- SMATRIX.C ---------------*/ int createsparse(EN_Project *pr); /* Creates sparse matrix */ diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index 84b67c87..00b0b820 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -18,24 +18,21 @@ HYDCOEFFS.C -- hydraulic coefficients for the EPANET Program #include "funcs.h" // Constants used for computing Darcy-Weisbach friction factor -const double A1 = 3.14159265358979323850e+03; // 1000*PI -const double A2 = 1.57079632679489661930e+03; // 500*PI -const double A3 = 5.02654824574366918160e+01; // 16*PI -const double A4 = 6.28318530717958647700e+00; // 2*PI -const double A8 = 4.61841319859066668690e+00; // 5.74*(PI/4)^.9 -const double A9 = -8.68588963806503655300e-01; // -2/ln(10) -const double AA = -1.5634601348517065795e+00; // -2*.9*2/ln(10) -const double AB = 3.28895476345399058690e-03; // 5.74/(4000^.9) -const double AC = -5.14214965799093883760e-03; // AA*AB - -// Cutoff flow for using linear head loss relation -const double Q_CUTOFF = 1.0e-5; +const double A1 = 0.314159265359e04; // 1000*PI +const double A2 = 0.157079632679e04; // 500*PI +const double A3 = 0.502654824574e02; // 16*PI +const double A4 = 6.283185307; // 2*PI +const double A8 = 4.61841319859; // 5.74*(PI/4)^.9 +const double A9 = -8.685889638e-01; // -2/ln(10) +const double AA = -1.5634601348; // -2*.9*2/ln(10) +const double AB = 3.28895476345e-03; // 5.74/(4000^.9) +const double AC = -5.14214965799e-03; // AA*AB // External functions //void resistcoeff(EN_Project *pr, int k); //void headlosscoeffs(EN_Project *pr); //void matrixcoeffs(EN_Project *pr); -//void emitheadloss(EN_Project *pr, int i, double *hloss, double *dhdq); +//double emitflowchange(EN_Project *pr, int i); //double demandflowchange(EN_Project *pr, int i, double dp, double n); //void demandparams(EN_Project *pr, double *dp, double *n); @@ -50,7 +47,7 @@ static void demandheadloss(double d, double dfull, double dp, static void pipecoeff(EN_Project *pr, int k); static void DWpipecoeff(EN_Project *pr, int k); -static double frictionFactor(double q, double e, double s, double *dfdq); +static double frictionFactor(EN_Project *pr, int k, double *dfdq); static void pumpcoeff(EN_Project *pr, int k); static void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r); @@ -198,8 +195,7 @@ void linkcoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes coefficients contributed by links to the -** linearized system of hydraulic equations. +** Purpose: computes matrix coefficients for links **-------------------------------------------------------------- */ { @@ -218,34 +214,32 @@ void linkcoeffs(EN_Project *pr) n1 = link->N1; // Start node of link n2 = link->N2; // End node of link - // Update nodal flow balance (X_tmp) - // (Flow out of node is (-), flow into node is (+)) + // Update net nodal inflows (X_tmp), solution matrix (A) and RHS array (F) + // (Use covention that flow out of node is (-), flow into node is (+)) hyd->X_tmp[n1] -= hyd->LinkFlows[k]; hyd->X_tmp[n2] += hyd->LinkFlows[k]; - // Add to off-diagonal coeff. of linear system matrix + // Off-diagonal coeff. sol->Aij[sol->Ndx[k]] -= sol->P[k]; - // Update linear system coeffs. associated with start node n1 - // ... node n1 is junction + // Node n1 is junction if (n1 <= net->Njuncs) { sol->Aii[sol->Row[n1]] += sol->P[k]; // Diagonal coeff. sol->F[sol->Row[n1]] += sol->Y[k]; // RHS coeff. } - // ... node n1 is a tank/reservoir + // Node n1 is a tank/reservoir else sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]); - // Update linear system coeffs. associated with end node n2 - // ... node n2 is junction + // Node n2 is junction if (n2 <= net->Njuncs) { sol->Aii[sol->Row[n2]] += sol->P[k]; // Diagonal coeff. sol->F[sol->Row[n2]] -= sol->Y[k]; // RHS coeff. } - // ... node n2 is a tank/reservoir + // Node n2 is a tank/reservoir else sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]); } } @@ -256,8 +250,8 @@ void nodecoeffs(EN_Project *pr) **---------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: completes calculation of nodal flow balance array -** (X_tmp) & r.h.s. (F) of linearized hydraulic eqns. +** Purpose: completes calculation of nodal flow imbalance (X_tmp) +** & flow correction (F) arrays **---------------------------------------------------------------- */ { @@ -267,7 +261,7 @@ void nodecoeffs(EN_Project *pr) EN_Network *net = &pr->network; // For junction nodes, subtract demand flow from net - // flow balance & add flow balance to RHS array F + // flow imbalance & add imbalance to RHS array F. for (i = 1; i <= net->Njuncs; i++) { hyd->X_tmp[i] -= hyd->DemandFlows[i]; @@ -281,9 +275,8 @@ void valvecoeffs(EN_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 +** Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs +** whose status is not fixed to OPEN/CLOSED **-------------------------------------------------------------- */ { @@ -332,8 +325,7 @@ void emittercoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes coeffs. of the linearized hydraulic eqns. -** contributed by emitters. +** Purpose: computes matrix coeffs. for emitters ** ** Note: Emitters consist of a fictitious pipe connected to ** a fictitious reservoir whose elevation equals that @@ -342,8 +334,12 @@ void emittercoeffs(EN_Project *pr) **-------------------------------------------------------------- */ { - int i, row; - double hloss, hgrad; + int i; + double ke; + double p; + double q; + double y; + double z; hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; @@ -352,57 +348,46 @@ void emittercoeffs(EN_Project *pr) for (i = 1; i <= net->Njuncs; i++) { - // Skip junctions without emitters node = &net->Node[i]; if (node->Ke == 0.0) continue; - - // Find emitter head loss and gradient - emitheadloss(pr, i, &hloss, &hgrad); - - // Row of solution matrix - row = sol->Row[i]; - - // Addition to matrix diagonal & r.h.s - sol->Aii[row] += 1.0 / hgrad; - sol->F[row] += (hloss + node->El) / hgrad; - - // Update to node flow balance - hyd->X_tmp[i] -= hyd->EmitterFlows[i]; + ke = MAX(CSMALL, node->Ke); // emitter coeff. + q = hyd->EmitterFlows[i]; // emitter flow + z = ke * pow(ABS(q), hyd->Qexp); // emitter head loss + p = hyd->Qexp * z / ABS(q); // head loss gradient + if (p < hyd->RQtol) + { + p = 1.0 / hyd->RQtol; + } + else + { + p = 1.0 / p; // inverse head loss gradient + } + y = SGN(q)*z*p; // head loss / gradient + sol->Aii[sol->Row[i]] += p; // addition to main diagonal + sol->F[sol->Row[i]] += y + p * node->El; // addition to r.h.s. + hyd->X_tmp[i] -= q; // addition to net node inflow } } -void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad) +double emitflowchange(EN_Project *pr, int i) /* -**------------------------------------------------------------- +**-------------------------------------------------------------- ** Input: i = node index -** Output: hloss = head loss across node's emitter -** hgrad = head loss gradient -** Purpose: computes an emitters's head loss and gradient. -**------------------------------------------------------------- +** Output: returns change in flow at an emitter node +** Purpose: computes flow change at an emitter node +**-------------------------------------------------------------- */ { - double ke; - double q; + double ke, p; hydraulics_t *hyd = &pr->hydraulics; + Snode *node = &pr->network.Node[i]; - // Set adjusted emitter coeff. - ke = MAX(CSMALL, pr->network.Node[i].Ke); - - // Use linear head loss relation for small flow - q = hyd->EmitterFlows[i]; - if (fabs(q) <= Q_CUTOFF) - { - *hgrad = hyd->Qexp * ke * pow(Q_CUTOFF, hyd->Qexp - 1.0); - *hloss = (*hgrad) * q; - } - - // Otherwise use normal emitter function - else - { - *hgrad = hyd->Qexp * ke * pow(fabs(q), hyd->Qexp - 1.0); - *hloss = (*hgrad) * q / hyd->Qexp; - } + ke = MAX(CSMALL, node->Ke); + p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]), (hyd->Qexp - 1.0)); + if (p < hyd->RQtol) p = 1 / hyd->RQtol; + else p = 1.0 / p; + return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El)); } @@ -412,8 +397,8 @@ void demandparams(EN_Project *pr, double *dp, double *n) ** Input: none ** Output: dp = pressure range over which demands can vary ** n = exponent in head loss v. demand function -** Purpose: retrieves parameters that define a pressure -** dependent demand function. +** Purpose: retrieves parameters that define a pressure dependent +** demand function. **-------------------------------------------------------------- */ { @@ -442,13 +427,12 @@ void demandcoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes coeffs. of the linearized hydraulic eqns. -** contributed by pressure dependent demands. +** Purpose: computes matrix coeffs. for pressure dependent demands ** ** Note: Pressure dependent demands are modelled like emitters -** with Hloss = Preq * (D / Dfull)^(1/Pexp) +** with Hloss = Pserv * (D / Dfull)^(1/Pexp) ** where D (actual demand) is zero for negative pressure -** and is Dfull above pressure Preq. +** and is Dfull above pressure Pserv. **-------------------------------------------------------------- */ { @@ -558,23 +542,25 @@ void pipecoeff(EN_Project *pr, int k) **-------------------------------------------------------------- ** Input: k = link index ** Output: none -** Purpose: computes P & Y coefficients for pipe k. +** Purpose: computes P & Y coefficients for pipe k ** -** P = inverse head loss gradient = 1/hgrad -** Y = flow correction term = hloss / hgrad +** P = inverse head loss gradient = 1/(dh/dQ) +** Y = flow correction term = h*P **-------------------------------------------------------------- */ { - double hloss, // Head loss - hgrad, // Head loss gradient + double hpipe, // Normal head loss + hml, // Minor head loss ml, // Minor loss coeff. + p, // q*(dh/dq) q, // Abs. value of flow r; // Resistance coeff. hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; + Slink *link = &pr->network.Link[k]; - // For closed pipe use headloss formula: hloss = CBIG*q + // For closed pipe use headloss formula: h = CBIG*q if (hyd->LinkStatus[k] <= CLOSED) { sol->P[k] = 1.0 / CBIG; @@ -589,35 +575,31 @@ void pipecoeff(EN_Project *pr, int k) return; } - q = ABS(hyd->LinkFlows[k]); - ml = pr->network.Link[k].Km; - r = pr->network.Link[k].R; + // Evaluate headloss coefficients + q = ABS(hyd->LinkFlows[k]); // Absolute flow + ml = link->Km; // Minor loss coeff. + r = link->R; // Resistance coeff. - // Friction head loss - if (q <= Q_CUTOFF) - { - hgrad = hyd->Hexp * r * pow(Q_CUTOFF, hyd->Hexp-1.0); - hloss = q * hgrad; - } - else + // Use large P coefficient for small flow resistance product + if ( (r+ml)*q < hyd->RQtol) { - hgrad = hyd->Hexp * r * pow(q, hyd->Hexp - 1.0); - hloss = hgrad * q / hyd->Hexp; + sol->P[k] = 1.0 / hyd->RQtol; + sol->Y[k] = hyd->LinkFlows[k] / hyd->Hexp; + return; } - // Contribution of minor head loss + // Compute P and Y coefficients + hpipe = r*pow(q, hyd->Hexp); // Friction head loss + p = hyd->Hexp*hpipe; // Q*dh(friction)/dQ if (ml > 0.0) { - hloss += ml * q * q; - hgrad += 2.0 * ml * q; + hml = ml*q*q; // Minor head loss + p += 2.0*hml; // Q*dh(Total)/dQ } - - // Adjust head loss sign for flow direction - hloss *= SGN(hyd->LinkFlows[k]); - - // P and Y coeffs. - sol->P[k] = 1.0 / hgrad; - sol->Y[k] = hloss / hgrad; + else hml = 0.0; + p = hyd->LinkFlows[k] / p; // 1 / (dh/dQ) + sol->P[k] = ABS(p); + sol->Y[k] = p*(hpipe + hml); } @@ -636,85 +618,102 @@ void DWpipecoeff(EN_Project *pr, int k) Slink *link = &pr->network.Link[k]; double q = ABS(hyd->LinkFlows[k]); - double r = link->R; // Resistance coeff. - double ml = link->Km; // Minor loss coeff. - double e = link->Kc / link->Diam; // Relative roughness - double s = hyd->Viscos * link->Diam; // Viscosity / diameter - - double hloss, hgrad, f, dfdq, r1; - - // Compute head loss and its derivative - // ... use Hagen-Poiseuille formula for laminar flow (Re <= 2000) - if (q <= A2 * s) + double dfdq = 0.0; + double r, r1, f, ml, p, hloss; + + ml = link->Km; // Minor loss coeff. + r = link->R; // Resistance coeff. + f = frictionFactor(pr,k,&dfdq); // D-W friction factor + r1 = f*r+ml; + + // Use large P coefficient for small flow resistance product + if (r1*q < hyd->RQtol) { - r = 16.0 * PI * s * r; - hloss = hyd->LinkFlows[k] * (r + ml * q); - hgrad = r + 2.0 * ml * q; - } - - // ... otherwise use Darcy-Weisbach formula with friction factor - else - { - dfdq = 0.0; - f = frictionFactor(q, e, s, &dfdq); - r1 = f * r + ml; - hloss = r1 * q * hyd->LinkFlows[k]; - hgrad = (2.0 * r1 * q) + (dfdq * r * q * q); + sol->P[k] = 1.0/hyd->RQtol; + sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp; + return; } - + // Compute P and Y coefficients - sol->P[k] = 1.0 / hgrad; - sol->Y[k] = hloss / hgrad; + hloss = r1*SQR(q); // Total head loss + p = 2.0*r1*q; // |dHloss/dQ| + // + dfdq*r*q*q; // Ignore df/dQ term + p = 1.0 / p; + sol->P[k] = p; + sol->Y[k] = SGN(hyd->LinkFlows[k]) * hloss * p; } -double frictionFactor(double q, double e, double s, double *dfdq) +double frictionFactor(EN_Project *pr, int k, double *dfdq) /* **-------------------------------------------------------------- -** Input: q = |pipe flow| -** e = pipe roughness / diameter -** s = viscosity * pipe diameter -** Output: dfdq = derivative of friction factor w.r.t. flow -** Returns: pipe's friction factor +** Input: k = link index +** Output: returns friction factor and +** replaces dfdq (derivative of f w.r.t. flow) ** Purpose: computes Darcy-Weisbach friction factor and its ** derivative as a function of Reynolds Number (Re). +** +** Note: Current formulas for dfdq need to be corrected +** so dfdq returned as 0. **-------------------------------------------------------------- */ { - double f; // friction factor - double x1, x2, x3, x4, - y1, y2, y3, - fa, fb, r; - double w = q / s; // Re*Pi/4 - - // For Re >= 4000 use Swamee & Jain approximation - // of the Colebrook-White Formula - if ( w >= A1 ) - { - y1 = A8 / pow(w, 0.9); - y2 = e / 3.7 + y1; - y3 = A9 * log(y2); - f = 1.0 / (y3*y3); - *dfdq = 1.8 * f * y1 * A9 / y2 / y3 / q; - } + double q, // Abs. value of flow + f; // Friction factor + double x1,x2,x3,x4, + y1,y2,y3, + fa,fb,r; + double s,w; + + hydraulics_t *hyd = &pr->hydraulics; + Slink *link = &pr->network.Link[k]; + + *dfdq = 0.0; + if (link->Type > EN_PIPE) + return(1.0); // Only apply to pipes + q = ABS(hyd->LinkFlows[k]); + s = hyd->Viscos * link->Diam; + w = q/s; // w = Re(Pi/4) + + // For Re >= 4000 use Colebrook Formula + if (w >= A1) + { + y1 = A8/pow(w,0.9); + y2 = link->Kc/(3.7*link->Diam) + y1; + y3 = A9*log(y2); + f = 1.0/SQR(y3); + /* *dfdq = (2.0+AA*y1/(y2*y3))*f; */ /* df/dq */ + } + + // For Re > 2000 use Interpolation Formula + else if (w > A2) + { + y2 = link->Kc/(3.7*link->Diam) + AB; + y3 = A9*log(y2); + fa = 1.0/SQR(y3); + fb = (2.0+AC/(y2*y3))*fa; + r = w/A2; + x1 = 7.0*fa - fb; + x2 = 0.128 - 17.0*fa + 2.5*fb; + x3 = -0.128 + 13.0*fa - (fb+fb); + x4 = r*(0.032 - 3.0*fa + 0.5*fb); + f = x1 + r*(x2 + r*(x3+x4)); + /* *dfdq = (x1 + x1 + r*(3.0*x2 + r*(4.0*x3 + 5.0*x4))); */ + } + + // For Re > 8 (Laminar Flow) use Hagen-Poiseuille Formula + else if (w > A4) + { + f = A3*s/q; // 16 * PI * Viscos * Diam / Flow + /* *dfdq = A3*s; */ + } + else + { + f = 8.0; + *dfdq = 0.0; + } + return(f); - // Use interpolating polynomials developed by - // E. Dunlop for transition flow from 2000 < Re < 4000. - else - { - y2 = e / 3.7 + AB; - y3 = A9 * log(y2); - fa = 1.0 / (y3*y3); - fb = (2.0 + AC / (y2*y3)) * fa; - r = w / A2; - x1 = 7.0 * fa - fb; - x2 = 0.128 - 17.0 * fa + 2.5 * fb; - x3 = -0.128 + 13.0 * fa - (fb + fb); - x4 = 0.032 - 3.0 * fa + 0.5 *fb; - f = x1 + r * (x2 + r * (x3 + r * x4)); - *dfdq = (x2 + r * (2.0 * x3 + r * 3.0 * x4)) / s / A2; - } - return f; } @@ -731,17 +730,13 @@ void pumpcoeff(EN_Project *pr, int k) double h0, // Shutoff head q, // Abs. value of flow r, // Flow resistance coeff. - n, // Flow exponent coeff. - setting, // Pump speed setting - hloss, // Head loss across pump - hgrad; // Head loss gradient - + n; // Flow exponent coeff. hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; - Spump *pump; + double setting = hyd->LinkSetting[k]; + Spump *pump; // Use high resistance pipe if pump closed or cannot deliver head - setting = hyd->LinkSetting[k]; if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) { sol->P[k] = 1.0 / CBIG; @@ -749,54 +744,35 @@ void pumpcoeff(EN_Project *pr, int k) return; } - // Obtain reference to pump object & its speed setting q = ABS(hyd->LinkFlows[k]); + q = MAX(q, TINY); + + // Obtain reference to pump object p = findpump(&pr->network, k); pump = &pr->network.Pump[p]; - // Get pump curve coefficients for custom pump curve - // (Other pump types have pre-determined coeffs.) + // Get pump curve coefficients for custom pump curve. if (pump->Ptype == CUSTOM) { // Find intercept (h0) & slope (r) of pump curve // line segment which contains speed-adjusted flow. curvecoeff(pr, pump->Hcurve, q / setting, &h0, &r); - // Determine head loss coefficients (negative sign - // converts from pump curve's head gain to head loss) + // Determine head loss coefficients. pump->H0 = -h0; pump->R = -r; pump->N = 1.0; - - // Compute head loss and its gradient - hgrad = pump->R * setting ; - hloss = pump->H0 * SQR(setting) + hgrad * hyd->LinkFlows[k]; - } - else - { - // Adjust head loss coefficients for pump speed - h0 = SQR(setting) * pump->H0; - n = pump->N; - r = pump->R * pow(setting, 2.0 - n); - - // Compute head loss and its gradient - // ... use linear approx. to pump curve for small flows - if (q <= Q_CUTOFF) - { - hgrad = r * pow(Q_CUTOFF, n - 1.0); - hloss = h0 + hgrad * hyd->LinkFlows[k]; - } - // ... use original pump curve for normal flows - else - { - hgrad = n * r * pow(q, n - 1.0); - hloss = h0 + hgrad * hyd->LinkFlows[k] / n; - } } - // P and Y coeffs. - sol->P[k] = 1.0 / hgrad; - sol->Y[k] = hloss / hgrad; + // Adjust head loss coefficients for pump speed. + h0 = SQR(setting) * pump->H0; + n = pump->N; + r = pump->R * pow(setting, 2.0 - n); + if (n != 1.0) r = n * r * pow(q, n - 1.0); + + // Compute inverse headloss gradient (P) and flow correction factor (Y) + sol->P[k] = 1.0 / MAX(r, hyd->RQtol); + sol->Y[k] = hyd->LinkFlows[k] / n + sol->P[k] * h0; } @@ -849,10 +825,9 @@ void gpvcoeff(EN_Project *pr, int k) **-------------------------------------------------------------- */ { - int i; - double h0, // Intercept of head loss curve segment - r, // Slope of head loss curve segment - q; // Abs. value of flow + double h0, // Headloss curve intercept + q, // Abs. value of flow + r; // Flow resistance coeff. hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; @@ -860,25 +835,19 @@ void gpvcoeff(EN_Project *pr, int k) // Treat as a pipe if valve closed if (hyd->LinkStatus[k] == CLOSED) valvecoeff(pr, k); - // Otherwise utilize segment of head loss curve - // bracketing current flow (curve index is stored - // in valve's setting) + // Otherwise utilize headloss curve + // whose index is stored in K else { - // Index of valve's head loss curve - i = (int)ROUND(hyd->LinkSetting[k]); - - // Adjusted flow rate + // Find slope & intercept of headloss curve. q = ABS(hyd->LinkFlows[k]); q = MAX(q, TINY); + curvecoeff(pr, (int)ROUND(hyd->LinkSetting[k]), q, &h0, &r); - // Intercept and slope of curve segment containing q - curvecoeff(pr, i, q, &h0, &r); - r = MAX(r, TINY); - - // Resulting P and Y coeffs. - sol->P[k] = 1.0 / r; - sol->Y[k] = (h0 / r + q) * SGN(hyd->LinkFlows[k]); + // Compute inverse headloss gradient (P) + // and flow correction factor (Y). + sol->P[k] = 1.0 / MAX(r, hyd->RQtol); + sol->Y[k] = sol->P[k] * (h0 + r*q) * SGN(hyd->LinkFlows[k]); } } @@ -1115,14 +1084,14 @@ void valvecoeff(EN_Project *pr, int k) **-------------------------------------------------------------- */ { - double flow, q, dhdq; + double p; EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; Slink *link = &net->Link[k]; - flow = hyd->LinkFlows[k]; + double flow = hyd->LinkFlows[k]; // Valve is closed. Use a very small matrix coeff. if (hyd->LinkStatus[k] <= CLOSED) @@ -1135,26 +1104,14 @@ void valvecoeff(EN_Project *pr, int k) // Account for any minor headloss through the valve if (link->Km > 0.0) { - // Adjust for very small flow - q = MAX(fabs(flow), Q_CUTOFF); - - // Gradient of minor head loss formula - dhdq = 2.0 * link->Km * q; - sol->P[k] = 1.0 / dhdq; - - // Y = head loss / gradient - // Use linear head loss relation for low flow - if (q == Q_CUTOFF) sol->Y[k] = flow; - - // Otherwise use normal minor loss formula - else sol->Y[k] = flow / 2.0; + p = 2.0 * link->Km * fabs(flow); + if (p < hyd->RQtol) p = hyd->RQtol; + sol->P[k] = 1.0 / p; + sol->Y[k] = flow / 2.0; } - - // If no minor loss coeff. specified use a - // low resistance linear head loss relation else { - sol->P[k] = 1.0 / CSMALL; + sol->P[k] = 1.0 / hyd->RQtol; sol->Y[k] = flow; } } diff --git a/src/hydsolver.c b/src/hydsolver.c index 63ed951c..7e6f0d65 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -494,34 +494,30 @@ void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum, **---------------------------------------------------------------- */ { - int i; - double hloss, hgrad, dh, dq; + double dq; + int k; EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; // Examine each network junction - for (i = 1; i <= net->Njuncs; i++) + for (k = 1; k <= net->Njuncs; k++) { // Skip junction if it does not have an emitter - if (net->Node[i].Ke == 0.0) continue; + if (net->Node[k].Ke == 0.0) continue; - // Find emitter head loss and gradient - emitheadloss(pr, i, &hloss, &hgrad); - - // Find emitter flow change - dh = hyd->NodeHead[i] - net->Node[i].El; - dq = (hloss - dh) / hgrad; - hyd->EmitterFlows[i] -= dq; + // Find emitter flow change (see hydcoeffs.c) + dq = emitflowchange(pr, k); + hyd->EmitterFlows[k] -= dq; // Update system flow summation - *qsum += ABS(hyd->EmitterFlows[i]); + *qsum += ABS(hyd->EmitterFlows[k]); *dqsum += ABS(dq); // Update identity of element with max. flow change if (ABS(dq) > hbal->maxflowchange) { hbal->maxflowchange = ABS(dq); - hbal->maxflownode = i; + hbal->maxflownode = k; hbal->maxflowlink = -1; } } From dab7be844690da62a9e84f73b0eb59b2e0d9f415 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Tue, 11 Sep 2018 13:15:15 -0400 Subject: [PATCH 3/6] Updates to fix zero flow (#164) and D-W eqn. (#199) issues --- src/funcs.h | 11 +- src/hydcoeffs.c | 464 ++++++++++++++++++++++++++---------------------- src/hydsolver.c | 22 ++- 3 files changed, 275 insertions(+), 222 deletions(-) diff --git a/src/funcs.h b/src/funcs.h index b961a82b..52b2d70c 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -178,11 +178,12 @@ int hydsolve(EN_Project *pr, int *,double *); /* Solves network equations void resistcoeff(EN_Project *pr, int k); /* Finds pipe flow resistance */ void headlosscoeffs(EN_Project *pr); // Finds link head loss coeffs. void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. */ -double emitflowchange(EN_Project *pr, int i); /* Change in emitter outflow */ -double demandflowchange(EN_Project *pr, int i, // Change in demand outflow - double dp, double n); -void demandparams(EN_Project *pr, double *dp, // PDA function parameters - double *n); +void emitheadloss(EN_Project *pr, int, // Finds emitter head loss + double *, double *); +double demandflowchange(EN_Project *pr, int, // Change in demand outflow + double, double); +void demandparams(EN_Project *pr, double *, // PDA function parameters + double *); /* ----------- SMATRIX.C ---------------*/ int createsparse(EN_Project *pr); /* Creates sparse matrix */ diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index 00b0b820..22a8b2dc 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -18,21 +18,24 @@ HYDCOEFFS.C -- hydraulic coefficients for the EPANET Program #include "funcs.h" // Constants used for computing Darcy-Weisbach friction factor -const double A1 = 0.314159265359e04; // 1000*PI -const double A2 = 0.157079632679e04; // 500*PI -const double A3 = 0.502654824574e02; // 16*PI -const double A4 = 6.283185307; // 2*PI -const double A8 = 4.61841319859; // 5.74*(PI/4)^.9 -const double A9 = -8.685889638e-01; // -2/ln(10) -const double AA = -1.5634601348; // -2*.9*2/ln(10) -const double AB = 3.28895476345e-03; // 5.74/(4000^.9) -const double AC = -5.14214965799e-03; // AA*AB +const double A1 = 3.14159265358979323850e+03; // 1000*PI +const double A2 = 1.57079632679489661930e+03; // 500*PI +const double A3 = 5.02654824574366918160e+01; // 16*PI +const double A4 = 6.28318530717958647700e+00; // 2*PI +const double A8 = 4.61841319859066668690e+00; // 5.74*(PI/4)^.9 +const double A9 = -8.68588963806503655300e-01; // -2/ln(10) +const double AA = -1.5634601348517065795e+00; // -2*.9*2/ln(10) +const double AB = 3.28895476345399058690e-03; // 5.74/(4000^.9) +const double AC = -5.14214965799093883760e-03; // AA*AB + +// Cutoff flow for using linear head loss relation +const double Q_CUTOFF = 1.0e-5; // External functions //void resistcoeff(EN_Project *pr, int k); //void headlosscoeffs(EN_Project *pr); //void matrixcoeffs(EN_Project *pr); -//double emitflowchange(EN_Project *pr, int i); +//void emitheadloss(EN_Project *pr, int i, double *hloss, double *dhdq); //double demandflowchange(EN_Project *pr, int i, double dp, double n); //void demandparams(EN_Project *pr, double *dp, double *n); @@ -47,7 +50,7 @@ static void demandheadloss(double d, double dfull, double dp, static void pipecoeff(EN_Project *pr, int k); static void DWpipecoeff(EN_Project *pr, int k); -static double frictionFactor(EN_Project *pr, int k, double *dfdq); +static double frictionFactor(double q, double e, double s, double *dfdq); static void pumpcoeff(EN_Project *pr, int k); static void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r); @@ -195,7 +198,8 @@ void linkcoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes matrix coefficients for links +** Purpose: computes coefficients contributed by links to the +** linearized system of hydraulic equations. **-------------------------------------------------------------- */ { @@ -214,32 +218,34 @@ void linkcoeffs(EN_Project *pr) n1 = link->N1; // Start node of link n2 = link->N2; // End node of link - // Update net nodal inflows (X_tmp), solution matrix (A) and RHS array (F) - // (Use covention that flow out of node is (-), flow into node is (+)) + // Update nodal flow balance (X_tmp) + // (Flow out of node is (-), flow into node is (+)) hyd->X_tmp[n1] -= hyd->LinkFlows[k]; hyd->X_tmp[n2] += hyd->LinkFlows[k]; - // Off-diagonal coeff. + // Add to off-diagonal coeff. of linear system matrix sol->Aij[sol->Ndx[k]] -= sol->P[k]; - // Node n1 is junction + // Update linear system coeffs. associated with start node n1 + // ... node n1 is junction if (n1 <= net->Njuncs) { sol->Aii[sol->Row[n1]] += sol->P[k]; // Diagonal coeff. sol->F[sol->Row[n1]] += sol->Y[k]; // RHS coeff. } - // Node n1 is a tank/reservoir + // ... node n1 is a tank/reservoir else sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]); - // Node n2 is junction + // Update linear system coeffs. associated with end node n2 + // ... node n2 is junction if (n2 <= net->Njuncs) { sol->Aii[sol->Row[n2]] += sol->P[k]; // Diagonal coeff. sol->F[sol->Row[n2]] -= sol->Y[k]; // RHS coeff. } - // Node n2 is a tank/reservoir + // ... node n2 is a tank/reservoir else sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]); } } @@ -250,8 +256,8 @@ void nodecoeffs(EN_Project *pr) **---------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: completes calculation of nodal flow imbalance (X_tmp) -** & flow correction (F) arrays +** Purpose: completes calculation of nodal flow balance array +** (X_tmp) & r.h.s. (F) of linearized hydraulic eqns. **---------------------------------------------------------------- */ { @@ -261,7 +267,7 @@ void nodecoeffs(EN_Project *pr) EN_Network *net = &pr->network; // For junction nodes, subtract demand flow from net - // flow imbalance & add imbalance to RHS array F. + // flow balance & add flow balance to RHS array F for (i = 1; i <= net->Njuncs; i++) { hyd->X_tmp[i] -= hyd->DemandFlows[i]; @@ -275,8 +281,9 @@ void valvecoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs -** whose status is not fixed to OPEN/CLOSED +** Purpose: computes coeffs. of the linearized hydraulic eqns. +** contributed by PRVs, PSVs & FCVs whose status is +** not fixed to OPEN/CLOSED **-------------------------------------------------------------- */ { @@ -325,7 +332,8 @@ void emittercoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes matrix coeffs. for emitters +** Purpose: computes coeffs. of the linearized hydraulic eqns. +** contributed by emitters. ** ** Note: Emitters consist of a fictitious pipe connected to ** a fictitious reservoir whose elevation equals that @@ -334,12 +342,8 @@ void emittercoeffs(EN_Project *pr) **-------------------------------------------------------------- */ { - int i; - double ke; - double p; - double q; - double y; - double z; + int i, row; + double hloss, hgrad; hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; @@ -348,46 +352,57 @@ void emittercoeffs(EN_Project *pr) for (i = 1; i <= net->Njuncs; i++) { + // Skip junctions without emitters node = &net->Node[i]; if (node->Ke == 0.0) continue; - ke = MAX(CSMALL, node->Ke); // emitter coeff. - q = hyd->EmitterFlows[i]; // emitter flow - z = ke * pow(ABS(q), hyd->Qexp); // emitter head loss - p = hyd->Qexp * z / ABS(q); // head loss gradient - if (p < hyd->RQtol) - { - p = 1.0 / hyd->RQtol; - } - else - { - p = 1.0 / p; // inverse head loss gradient - } - y = SGN(q)*z*p; // head loss / gradient - sol->Aii[sol->Row[i]] += p; // addition to main diagonal - sol->F[sol->Row[i]] += y + p * node->El; // addition to r.h.s. - hyd->X_tmp[i] -= q; // addition to net node inflow + + // Find emitter head loss and gradient + emitheadloss(pr, i, &hloss, &hgrad); + + // Row of solution matrix + row = sol->Row[i]; + + // Addition to matrix diagonal & r.h.s + sol->Aii[row] += 1.0 / hgrad; + sol->F[row] += (hloss + node->El) / hgrad; + + // Update to node flow balance + hyd->X_tmp[i] -= hyd->EmitterFlows[i]; } } -double emitflowchange(EN_Project *pr, int i) +void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad) /* -**-------------------------------------------------------------- +**------------------------------------------------------------- ** Input: i = node index -** Output: returns change in flow at an emitter node -** Purpose: computes flow change at an emitter node -**-------------------------------------------------------------- +** Output: hloss = head loss across node's emitter +** hgrad = head loss gradient +** Purpose: computes an emitters's head loss and gradient. +**------------------------------------------------------------- */ { - double ke, p; + double ke; + double q; hydraulics_t *hyd = &pr->hydraulics; - Snode *node = &pr->network.Node[i]; - ke = MAX(CSMALL, node->Ke); - p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]), (hyd->Qexp - 1.0)); - if (p < hyd->RQtol) p = 1 / hyd->RQtol; - else p = 1.0 / p; - return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El)); + // Set adjusted emitter coeff. + ke = MAX(CSMALL, pr->network.Node[i].Ke); + + // Use linear head loss relation for small flow + q = hyd->EmitterFlows[i]; + if (fabs(q) <= Q_CUTOFF) + { + *hgrad = ke * pow(Q_CUTOFF, hyd->Qexp) / Q_CUTOFF; + *hloss = (*hgrad) * q; + } + + // Otherwise use normal emitter function + else + { + *hgrad = hyd->Qexp * ke * pow(fabs(q), hyd->Qexp - 1.0); + *hloss = (*hgrad) * q / hyd->Qexp; + } } @@ -397,8 +412,8 @@ void demandparams(EN_Project *pr, double *dp, double *n) ** Input: none ** Output: dp = pressure range over which demands can vary ** n = exponent in head loss v. demand function -** Purpose: retrieves parameters that define a pressure dependent -** demand function. +** Purpose: retrieves parameters that define a pressure +** dependent demand function. **-------------------------------------------------------------- */ { @@ -427,12 +442,13 @@ void demandcoeffs(EN_Project *pr) **-------------------------------------------------------------- ** Input: none ** Output: none -** Purpose: computes matrix coeffs. for pressure dependent demands +** Purpose: computes coeffs. of the linearized hydraulic eqns. +** contributed by pressure dependent demands. ** ** Note: Pressure dependent demands are modelled like emitters -** with Hloss = Pserv * (D / Dfull)^(1/Pexp) +** with Hloss = Preq * (D / Dfull)^(1/Pexp) ** where D (actual demand) is zero for negative pressure -** and is Dfull above pressure Pserv. +** and is Dfull above pressure Preq. **-------------------------------------------------------------- */ { @@ -524,7 +540,7 @@ void demandheadloss(double d, double dfull, double dp, double n, // Use linear head loss function for near zero demand else if (r < EPS) { - *hgrad = dp * pow(EPS, n - 1.0) / dfull; + *hgrad = dp * pow(EPS, n) / dfull / EPS; *hloss = (*hgrad) * d; } @@ -542,25 +558,23 @@ void pipecoeff(EN_Project *pr, int k) **-------------------------------------------------------------- ** Input: k = link index ** Output: none -** Purpose: computes P & Y coefficients for pipe k +** Purpose: computes P & Y coefficients for pipe k. ** -** P = inverse head loss gradient = 1/(dh/dQ) -** Y = flow correction term = h*P +** P = inverse head loss gradient = 1/hgrad +** Y = flow correction term = hloss / hgrad **-------------------------------------------------------------- */ { - double hpipe, // Normal head loss - hml, // Minor head loss + double hloss, // Head loss + hgrad, // Head loss gradient ml, // Minor loss coeff. - p, // q*(dh/dq) q, // Abs. value of flow r; // Resistance coeff. hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; - Slink *link = &pr->network.Link[k]; - // For closed pipe use headloss formula: h = CBIG*q + // For closed pipe use headloss formula: hloss = CBIG*q if (hyd->LinkStatus[k] <= CLOSED) { sol->P[k] = 1.0 / CBIG; @@ -575,31 +589,37 @@ void pipecoeff(EN_Project *pr, int k) return; } - // Evaluate headloss coefficients - q = ABS(hyd->LinkFlows[k]); // Absolute flow - ml = link->Km; // Minor loss coeff. - r = link->R; // Resistance coeff. + q = ABS(hyd->LinkFlows[k]); + ml = pr->network.Link[k].Km; + r = pr->network.Link[k].R; - // Use large P coefficient for small flow resistance product - if ( (r+ml)*q < hyd->RQtol) + // Friction head loss + // ... use linear relation for small flows + if (q <= Q_CUTOFF) { - sol->P[k] = 1.0 / hyd->RQtol; - sol->Y[k] = hyd->LinkFlows[k] / hyd->Hexp; - return; + hgrad = r * pow(Q_CUTOFF, hyd->Hexp) / Q_CUTOFF; + hloss = hgrad * q; + } + // ... use original formula for other flows + else + { + hgrad = hyd->Hexp * r * pow(q, hyd->Hexp - 1.0); + hloss = hgrad * q / hyd->Hexp; } - // Compute P and Y coefficients - hpipe = r*pow(q, hyd->Hexp); // Friction head loss - p = hyd->Hexp*hpipe; // Q*dh(friction)/dQ + // Contribution of minor head loss if (ml > 0.0) { - hml = ml*q*q; // Minor head loss - p += 2.0*hml; // Q*dh(Total)/dQ + hloss += ml * q * q; + hgrad += 2.0 * ml * q; } - else hml = 0.0; - p = hyd->LinkFlows[k] / p; // 1 / (dh/dQ) - sol->P[k] = ABS(p); - sol->Y[k] = p*(hpipe + hml); + + // Adjust head loss sign for flow direction + hloss *= SGN(hyd->LinkFlows[k]); + + // P and Y coeffs. + sol->P[k] = 1.0 / hgrad; + sol->Y[k] = hloss / hgrad; } @@ -618,102 +638,85 @@ void DWpipecoeff(EN_Project *pr, int k) Slink *link = &pr->network.Link[k]; double q = ABS(hyd->LinkFlows[k]); - double dfdq = 0.0; - double r, r1, f, ml, p, hloss; - - ml = link->Km; // Minor loss coeff. - r = link->R; // Resistance coeff. - f = frictionFactor(pr,k,&dfdq); // D-W friction factor - r1 = f*r+ml; - - // Use large P coefficient for small flow resistance product - if (r1*q < hyd->RQtol) + double r = link->R; // Resistance coeff. + double ml = link->Km; // Minor loss coeff. + double e = link->Kc / link->Diam; // Relative roughness + double s = hyd->Viscos * link->Diam; // Viscosity / diameter + + double hloss, hgrad, f, dfdq, r1; + + // Compute head loss and its derivative + // ... use Hagen-Poiseuille formula for laminar flow (Re <= 2000) + if (q <= A2 * s) { - sol->P[k] = 1.0/hyd->RQtol; - sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp; - return; + r = 16.0 * PI * s * r; + hloss = hyd->LinkFlows[k] * (r + ml * q); + hgrad = r + 2.0 * ml * q; } - + + // ... otherwise use Darcy-Weisbach formula with friction factor + else + { + dfdq = 0.0; + f = frictionFactor(q, e, s, &dfdq); + r1 = f * r + ml; + hloss = r1 * q * hyd->LinkFlows[k]; + hgrad = (2.0 * r1 * q) + (dfdq * r * q * q); + } + // Compute P and Y coefficients - hloss = r1*SQR(q); // Total head loss - p = 2.0*r1*q; // |dHloss/dQ| - // + dfdq*r*q*q; // Ignore df/dQ term - p = 1.0 / p; - sol->P[k] = p; - sol->Y[k] = SGN(hyd->LinkFlows[k]) * hloss * p; + sol->P[k] = 1.0 / hgrad; + sol->Y[k] = hloss / hgrad; } -double frictionFactor(EN_Project *pr, int k, double *dfdq) +double frictionFactor(double q, double e, double s, double *dfdq) /* **-------------------------------------------------------------- -** Input: k = link index -** Output: returns friction factor and -** replaces dfdq (derivative of f w.r.t. flow) +** Input: q = |pipe flow| +** e = pipe roughness / diameter +** s = viscosity * pipe diameter +** Output: dfdq = derivative of friction factor w.r.t. flow +** Returns: pipe's friction factor ** Purpose: computes Darcy-Weisbach friction factor and its ** derivative as a function of Reynolds Number (Re). -** -** Note: Current formulas for dfdq need to be corrected -** so dfdq returned as 0. **-------------------------------------------------------------- */ { - double q, // Abs. value of flow - f; // Friction factor - double x1,x2,x3,x4, - y1,y2,y3, - fa,fb,r; - double s,w; - - hydraulics_t *hyd = &pr->hydraulics; - Slink *link = &pr->network.Link[k]; - - *dfdq = 0.0; - if (link->Type > EN_PIPE) - return(1.0); // Only apply to pipes - q = ABS(hyd->LinkFlows[k]); - s = hyd->Viscos * link->Diam; - w = q/s; // w = Re(Pi/4) - - // For Re >= 4000 use Colebrook Formula - if (w >= A1) - { - y1 = A8/pow(w,0.9); - y2 = link->Kc/(3.7*link->Diam) + y1; - y3 = A9*log(y2); - f = 1.0/SQR(y3); - /* *dfdq = (2.0+AA*y1/(y2*y3))*f; */ /* df/dq */ - } - - // For Re > 2000 use Interpolation Formula - else if (w > A2) - { - y2 = link->Kc/(3.7*link->Diam) + AB; - y3 = A9*log(y2); - fa = 1.0/SQR(y3); - fb = (2.0+AC/(y2*y3))*fa; - r = w/A2; - x1 = 7.0*fa - fb; - x2 = 0.128 - 17.0*fa + 2.5*fb; - x3 = -0.128 + 13.0*fa - (fb+fb); - x4 = r*(0.032 - 3.0*fa + 0.5*fb); - f = x1 + r*(x2 + r*(x3+x4)); - /* *dfdq = (x1 + x1 + r*(3.0*x2 + r*(4.0*x3 + 5.0*x4))); */ - } - - // For Re > 8 (Laminar Flow) use Hagen-Poiseuille Formula - else if (w > A4) - { - f = A3*s/q; // 16 * PI * Viscos * Diam / Flow - /* *dfdq = A3*s; */ - } - else - { - f = 8.0; - *dfdq = 0.0; - } - return(f); + double f; // friction factor + double x1, x2, x3, x4, + y1, y2, y3, + fa, fb, r; + double w = q / s; // Re*Pi/4 + + // For Re >= 4000 use Swamee & Jain approximation + // of the Colebrook-White Formula + if ( w >= A1 ) + { + y1 = A8 / pow(w, 0.9); + y2 = e / 3.7 + y1; + y3 = A9 * log(y2); + f = 1.0 / (y3*y3); + *dfdq = 1.8 * f * y1 * A9 / y2 / y3 / q; + } + // Use interpolating polynomials developed by + // E. Dunlop for transition flow from 2000 < Re < 4000. + else + { + y2 = e / 3.7 + AB; + y3 = A9 * log(y2); + fa = 1.0 / (y3*y3); + fb = (2.0 + AC / (y2*y3)) * fa; + r = w / A2; + x1 = 7.0 * fa - fb; + x2 = 0.128 - 17.0 * fa + 2.5 * fb; + x3 = -0.128 + 13.0 * fa - (fb + fb); + x4 = 0.032 - 3.0 * fa + 0.5 *fb; + f = x1 + r * (x2 + r * (x3 + r * x4)); + *dfdq = (x2 + r * (2.0 * x3 + r * 3.0 * x4)) / s / A2; + } + return f; } @@ -730,13 +733,17 @@ void pumpcoeff(EN_Project *pr, int k) double h0, // Shutoff head q, // Abs. value of flow r, // Flow resistance coeff. - n; // Flow exponent coeff. + n, // Flow exponent coeff. + setting, // Pump speed setting + hloss, // Head loss across pump + hgrad; // Head loss gradient + hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; - double setting = hyd->LinkSetting[k]; - Spump *pump; + Spump *pump; // Use high resistance pipe if pump closed or cannot deliver head + setting = hyd->LinkSetting[k]; if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) { sol->P[k] = 1.0 / CBIG; @@ -744,35 +751,54 @@ void pumpcoeff(EN_Project *pr, int k) return; } + // Obtain reference to pump object & its speed setting q = ABS(hyd->LinkFlows[k]); - q = MAX(q, TINY); - - // Obtain reference to pump object p = findpump(&pr->network, k); pump = &pr->network.Pump[p]; - // Get pump curve coefficients for custom pump curve. + // Get pump curve coefficients for custom pump curve + // (Other pump types have pre-determined coeffs.) if (pump->Ptype == CUSTOM) { // Find intercept (h0) & slope (r) of pump curve // line segment which contains speed-adjusted flow. curvecoeff(pr, pump->Hcurve, q / setting, &h0, &r); - // Determine head loss coefficients. + // Determine head loss coefficients (negative sign + // converts from pump curve's head gain to head loss) pump->H0 = -h0; pump->R = -r; pump->N = 1.0; - } - // Adjust head loss coefficients for pump speed. - h0 = SQR(setting) * pump->H0; - n = pump->N; - r = pump->R * pow(setting, 2.0 - n); - if (n != 1.0) r = n * r * pow(q, n - 1.0); + // Compute head loss and its gradient + hgrad = pump->R * setting ; + hloss = pump->H0 * SQR(setting) + hgrad * hyd->LinkFlows[k]; + } + else + { + // Adjust head loss coefficients for pump speed + h0 = SQR(setting) * pump->H0; + n = pump->N; + r = pump->R * pow(setting, 2.0 - n); + + // Compute head loss and its gradient + // ... use linear approx. to pump curve for small flows + if (q <= Q_CUTOFF) + { + hgrad = r * pow(Q_CUTOFF, n) / Q_CUTOFF; + hloss = h0 + hgrad * hyd->LinkFlows[k]; + } + // ... use original pump curve for normal flows + else + { + hgrad = n * r * pow(q, n - 1.0); + hloss = h0 + hgrad * hyd->LinkFlows[k] / n; + } + } - // Compute inverse headloss gradient (P) and flow correction factor (Y) - sol->P[k] = 1.0 / MAX(r, hyd->RQtol); - sol->Y[k] = hyd->LinkFlows[k] / n + sol->P[k] * h0; + // P and Y coeffs. + sol->P[k] = 1.0 / hgrad; + sol->Y[k] = hloss / hgrad; } @@ -825,9 +851,10 @@ void gpvcoeff(EN_Project *pr, int k) **-------------------------------------------------------------- */ { - double h0, // Headloss curve intercept - q, // Abs. value of flow - r; // Flow resistance coeff. + int i; + double h0, // Intercept of head loss curve segment + r, // Slope of head loss curve segment + q; // Abs. value of flow hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; @@ -835,19 +862,25 @@ void gpvcoeff(EN_Project *pr, int k) // Treat as a pipe if valve closed if (hyd->LinkStatus[k] == CLOSED) valvecoeff(pr, k); - // Otherwise utilize headloss curve - // whose index is stored in K + // Otherwise utilize segment of head loss curve + // bracketing current flow (curve index is stored + // in valve's setting) else { - // Find slope & intercept of headloss curve. + // Index of valve's head loss curve + i = (int)ROUND(hyd->LinkSetting[k]); + + // Adjusted flow rate q = ABS(hyd->LinkFlows[k]); q = MAX(q, TINY); - curvecoeff(pr, (int)ROUND(hyd->LinkSetting[k]), q, &h0, &r); - // Compute inverse headloss gradient (P) - // and flow correction factor (Y). - sol->P[k] = 1.0 / MAX(r, hyd->RQtol); - sol->Y[k] = sol->P[k] * (h0 + r*q) * SGN(hyd->LinkFlows[k]); + // Intercept and slope of curve segment containing q + curvecoeff(pr, i, q, &h0, &r); + r = MAX(r, TINY); + + // Resulting P and Y coeffs. + sol->P[k] = 1.0 / r; + sol->Y[k] = (h0 / r + q) * SGN(hyd->LinkFlows[k]); } } @@ -1084,14 +1117,14 @@ void valvecoeff(EN_Project *pr, int k) **-------------------------------------------------------------- */ { - double p; + double flow, q, y, hgrad; EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; Slink *link = &net->Link[k]; - double flow = hyd->LinkFlows[k]; + flow = hyd->LinkFlows[k]; // Valve is closed. Use a very small matrix coeff. if (hyd->LinkStatus[k] <= CLOSED) @@ -1104,14 +1137,29 @@ void valvecoeff(EN_Project *pr, int k) // Account for any minor headloss through the valve if (link->Km > 0.0) { - p = 2.0 * link->Km * fabs(flow); - if (p < hyd->RQtol) p = hyd->RQtol; - sol->P[k] = 1.0 / p; - sol->Y[k] = flow / 2.0; + // Adjust for very small flow + q = fabs(flow); + if (q <= Q_CUTOFF) + { + hgrad = link->Km * Q_CUTOFF; + y = flow; + } + else + { + hgrad = 2.0 * link->Km * q; + y = flow / 2.0; + } + + // P and Y coeffs. + sol->P[k] = 1.0 / hgrad; + sol->Y[k] = y; } + + // If no minor loss coeff. specified use a + // low resistance linear head loss relation else { - sol->P[k] = 1.0 / hyd->RQtol; + sol->P[k] = 1.0 / CSMALL; sol->Y[k] = flow; } } diff --git a/src/hydsolver.c b/src/hydsolver.c index 7e6f0d65..63ed951c 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -494,30 +494,34 @@ void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum, **---------------------------------------------------------------- */ { - double dq; - int k; + int i; + double hloss, hgrad, dh, dq; EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; // Examine each network junction - for (k = 1; k <= net->Njuncs; k++) + for (i = 1; i <= net->Njuncs; i++) { // Skip junction if it does not have an emitter - if (net->Node[k].Ke == 0.0) continue; + if (net->Node[i].Ke == 0.0) continue; - // Find emitter flow change (see hydcoeffs.c) - dq = emitflowchange(pr, k); - hyd->EmitterFlows[k] -= dq; + // Find emitter head loss and gradient + emitheadloss(pr, i, &hloss, &hgrad); + + // Find emitter flow change + dh = hyd->NodeHead[i] - net->Node[i].El; + dq = (hloss - dh) / hgrad; + hyd->EmitterFlows[i] -= dq; // Update system flow summation - *qsum += ABS(hyd->EmitterFlows[k]); + *qsum += ABS(hyd->EmitterFlows[i]); *dqsum += ABS(dq); // Update identity of element with max. flow change if (ABS(dq) > hbal->maxflowchange) { hbal->maxflowchange = ABS(dq); - hbal->maxflownode = k; + hbal->maxflownode = i; hbal->maxflowlink = -1; } } From e45a23c4efeb666a203428a1453efe7f30444550 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Wed, 12 Sep 2018 10:38:01 -0400 Subject: [PATCH 4/6] More robust method for handling zero/low flows (#164) --- src/hydcoeffs.c | 40 ++++++++++++++++++++++++++-------------- src/types.h | 3 ++- 2 files changed, 28 insertions(+), 15 deletions(-) diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index 22a8b2dc..ed0c258b 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -28,9 +28,6 @@ const double AA = -1.5634601348517065795e+00; // -2*.9*2/ln(10) const double AB = 3.28895476345399058690e-03; // 5.74/(4000^.9) const double AC = -5.14214965799093883760e-03; // AA*AB -// Cutoff flow for using linear head loss relation -const double Q_CUTOFF = 1.0e-5; - // External functions //void resistcoeff(EN_Project *pr, int k); //void headlosscoeffs(EN_Project *pr); @@ -74,10 +71,12 @@ void resistcoeff(EN_Project *pr, int k) */ { double e, d, L; - EN_Network *net = &pr->network; + + EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; Slink *link = &net->Link[k]; + link->Qa = 0.0; switch (link->Type) { // ... Link is a pipe. Compute resistance based on headloss formula. @@ -91,14 +90,20 @@ void resistcoeff(EN_Project *pr, int k) switch (hyd->Formflag) { case HW: + // ... resistance coeff. 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; case DW: link->R = L / 2.0 / 32.2 / d / SQR(PI * SQR(d) / 4.0); break; case CM: + // ... resistance coeff. link->R = SQR(4.0 * e / (1.49 * PI * SQR(d))) * pow((d / 4.0), -1.333) * L; + // ... flow below which linear head loss applies + link->Qa = hyd->RQtol / 2.0 / link->R; } break; @@ -384,16 +389,20 @@ void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad) { double ke; double q; + double qa; hydraulics_t *hyd = &pr->hydraulics; // Set adjusted emitter coeff. ke = MAX(CSMALL, pr->network.Node[i].Ke); + // Find flow below which head loss is linear + qa = pow(hyd->RQtol / ke / hyd->Qexp, 1.0 / (hyd->Qexp - 1.0)); + // Use linear head loss relation for small flow q = hyd->EmitterFlows[i]; - if (fabs(q) <= Q_CUTOFF) + if (fabs(q) <= qa) { - *hgrad = ke * pow(Q_CUTOFF, hyd->Qexp) / Q_CUTOFF; + *hgrad = hyd->RQtol; *hloss = (*hgrad) * q; } @@ -595,9 +604,9 @@ void pipecoeff(EN_Project *pr, int k) // Friction head loss // ... use linear relation for small flows - if (q <= Q_CUTOFF) + if (q <= pr->network.Link[k].Qa) { - hgrad = r * pow(Q_CUTOFF, hyd->Hexp) / Q_CUTOFF; + hgrad = hyd->RQtol; hloss = hgrad * q; } // ... use original formula for other flows @@ -735,6 +744,7 @@ void pumpcoeff(EN_Project *pr, int k) r, // Flow resistance coeff. n, // Flow exponent coeff. setting, // Pump speed setting + qa, // Flow limit for linear head loss hloss, // Head loss across pump hgrad; // Head loss gradient @@ -783,9 +793,10 @@ void pumpcoeff(EN_Project *pr, int k) // Compute head loss and its gradient // ... use linear approx. to pump curve for small flows - if (q <= Q_CUTOFF) + qa = pow(hyd->RQtol / n / r, 1.0 / (n - 1.0)); + if (q <= qa) { - hgrad = r * pow(Q_CUTOFF, n) / Q_CUTOFF; + hgrad = hyd->RQtol; hloss = h0 + hgrad * hyd->LinkFlows[k]; } // ... use original pump curve for normal flows @@ -1117,8 +1128,8 @@ void valvecoeff(EN_Project *pr, int k) **-------------------------------------------------------------- */ { - double flow, q, y, hgrad; - + double flow, q, y, qa, hgrad; + EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; @@ -1139,9 +1150,10 @@ void valvecoeff(EN_Project *pr, int k) { // Adjust for very small flow q = fabs(flow); - if (q <= Q_CUTOFF) + qa = hyd->RQtol / 2.0 / link->Km; + if (q <= qa) { - hgrad = link->Km * Q_CUTOFF; + hgrad = hyd->RQtol; y = flow; } else diff --git a/src/types.h b/src/types.h index 0a2d26ca..61c2cdb2 100755 --- a/src/types.h +++ b/src/types.h @@ -410,7 +410,8 @@ typedef struct /* LINK OBJECT */ double Kb; /* Bulk react. coeff */ double Kw; /* Wall react. coeff */ double R; /* Flow resistance */ - double Rc; /* Reaction cal */ + double Rc; /* Reaction coeff. */ + double Qa; // Low flow limit EN_LinkType Type; /* Link type */ StatType Stat; /* Initial status */ char Rpt; /* Reporting flag */ From 2a86252a1042a34b0f809cc25423cf3ca357b28b Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Thu, 13 Sep 2018 16:20:52 -0400 Subject: [PATCH 5/6] Updating benchmark --- tools/before-test.cmd | 4 ++-- tools/nrtest-epanet/report_diff.py | 14 ++++++++++---- tools/run-nrtest.cmd | 2 +- 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/tools/before-test.cmd b/tools/before-test.cmd index 0f64cd96..3e79017b 100644 --- a/tools/before-test.cmd +++ b/tools/before-test.cmd @@ -25,8 +25,8 @@ set SCRIPT_HOME=%~dp0 set TEST_HOME=%~1 -set EXAMPLES_VER=1.0.2-dev.2 -set BENCHMARK_VER=220dev2 +set EXAMPLES_VER=1.0.2-dev.3 +set BENCHMARK_VER=220dev3 set TESTFILES_URL=https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v%EXAMPLES_VER%.zip diff --git a/tools/nrtest-epanet/report_diff.py b/tools/nrtest-epanet/report_diff.py index b3f49e35..b82d746c 100644 --- a/tools/nrtest-epanet/report_diff.py +++ b/tools/nrtest-epanet/report_diff.py @@ -51,16 +51,22 @@ def _log_relative_error(q, c): ''' diff = np.subtract(q, c) tmp_c = np.copy(c) + # If ref value is small compute absolute error - tmp_c[np.fabs(tmp_c) < 1.0e-6] = 1.0 - + tmp_c[np.fabs(tmp_c) < 1.0e-6] = 1.0 re = np.fabs(diff)/np.fabs(tmp_c) + # If re is tiny set lre to number of digits re[re < 1.0e-7] = 1.0e-7 # If re is very large set lre to zero re[re > 2.0] = 1.0 - return np.negative(np.log10(re)) + lre = np.negative(np.log10(re)) + + # If lre is negative set to zero + lre[lre < 1.0] = 0.0 + + return lre def _print_diff(idx, lre, test, ref): @@ -71,7 +77,7 @@ def _print_diff(idx, lre, test, ref): diff_val = (test_val - ref_val) lre_val = (lre[idx[0]]) - print("Idx: %s\nSut: %f Ref: %f Diff: %f LRE: %f\n" + print("Idx: %s\nSut: %e Ref: %e Diff: %e LRE: %.2f\n" % (idx_val, test_val, ref_val, diff_val, lre_val)) diff --git a/tools/run-nrtest.cmd b/tools/run-nrtest.cmd index 8efb44bd..2fc2ef25 100644 --- a/tools/run-nrtest.cmd +++ b/tools/run-nrtest.cmd @@ -18,7 +18,7 @@ setlocal set NRTEST_SCRIPT_PATH=%~1 set TEST_SUITE_PATH=%~2 -set BENCHMARK_VER=220dev2 +set BENCHMARK_VER=220dev3 set NRTEST_EXECUTE_CMD=python %NRTEST_SCRIPT_PATH%\nrtest execute From d164293d8d8fc3e572c871136a3087692c777d9a Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Thu, 13 Sep 2018 16:48:45 -0400 Subject: [PATCH 6/6] Triggering appveyor build --- README.md | 1 - 1 file changed, 1 deletion(-) diff --git a/README.md b/README.md index 2c845383..85bf93e3 100755 --- a/README.md +++ b/README.md @@ -18,4 +18,3 @@ A step-by-step tutorial on how to contribute to EPANET using GitHub is also [ava __Note:__ This repository is not affiliated with, or endorsed by, the USEPA. For the last "official" release of EPANET (2.00.12 UI and Toolkit) please go to the [EPA's GitHub repo](https://github.com/USEPA/Water-Distribution-Network-Model) or [the USEPA website](http://www2.epa.gov/water-research/epanet). It is also not the graphical user interface version. This is the hydraulic and water quality solver engine. However, if you are interested in extending EPANET for academic, personal, or commercial use, then you've come to the right place. For community discussion, FAQ, and roadmapping of the project, go to the [Community Forum](http://community.wateranalytics.org/category/epanet). -