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

Reverted to legacy Darcy-Weisbach method

Browse Source
This commit is contained in:
Lew Rossman
2018-06-21 12:05:51 -04:00
parent 92e10a18b0
commit f5a1b9e518
2 changed files with 82 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

142
src/hydcoeffs.c
View File

@@ -41,7 +41,7 @@ static void emittercoeffs(EN_Project *pr);
static void pipecoeff(EN_Project *pr, int k); static void pipecoeff(EN_Project *pr, int k);
static void DWpipecoeff(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 pumpcoeff(EN_Project *pr, int k);
static void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r); static void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r);
@@ -192,15 +192,8 @@ void linkcoeffs(EN_Project *pr)
// Examine each link of network */ // Examine each link of network */
for (k = 1; k <= net->Nlinks; k++) for (k = 1; k <= net->Nlinks; k++)
{ {
// if (sol->P[k] == 0.0) continue; if (sol->P[k] == 0.0) continue;
link = &net->Link[k]; link = &net->Link[k];
switch (link->Type) {
case EN_PRV:
case EN_PSV:
case EN_FCV:
if (hyd->LinkSetting[k] != MISSING) continue;
}
n1 = link->N1; // Start node of link n1 = link->N1; // Start node of link
n2 = link->N2; // End node of link n2 = link->N2; // End node of link
@@ -437,88 +430,103 @@ void DWpipecoeff(EN_Project *pr, int k)
solver_t *sol = &hyd->solver; solver_t *sol = &hyd->solver;
Slink *link = &pr->network.Link[k]; Slink *link = &pr->network.Link[k];
double s = hyd->Viscos * link->Diam;
double q = ABS(hyd->LinkFlows[k]); double q = ABS(hyd->LinkFlows[k]);
double dfdq = 0.0; double dfdq = 0.0;
double e; // relative roughness height double r, r1, f, ml, p, hloss;
double r; // total resistance coeff.
double f; // friction factor
double hloss; // head loss
double hgrad; // head loss gradient
// ... use Hagen-Poiseuille formula for laminar flow (Re <= 2000) ml = link->Km; // Minor loss coeff.
if (q <= A2 * s) 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 * link->R; sol->P[k] = 1.0/hyd->RQtol;
hloss = q * (r + link->Km * q); sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp;
hgrad = r + 2.0 * link->Km * q; return;
} }
// ... use Colebrook formula for turbulent flow // Compute P and Y coefficients
else hloss = r1*SQR(q); // Total head loss
{ p = 2.0*r1*q; // |dHloss/dQ|
e = link->Kc / link->Diam; // + dfdq*r*q*q; // Ignore df/dQ term
f = frictionFactor(q, e, s, &dfdq); p = 1.0 / p;
r = f * link->R + link->Km; sol->P[k] = p;
hloss = r * q * hyd->LinkFlows[k]; sol->Y[k] = SGN(hyd->LinkFlows[k]) * hloss * p;
hgrad = (2.0 * r * q) + (dfdq * link->R * q * q);
}
// ... head loss has same sign as flow
hloss *= SGN(hyd->LinkFlows[k]);
// ... compute P and Y coeffs.
sol->P[k] = 1.0 / hgrad;
sol->Y[k] = hloss / hgrad;
} }
double frictionFactor(double q, double e, double s, double *dfdq) double frictionFactor(EN_Project *pr, int k, double *dfdq)
/* /*
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: q = flow rate (cfs) ** Input: k = link index
** e = pipe roughness / diameter ** Output: returns friction factor and
** s = viscosity * diameter ** replaces dfdq (derivative of f w.r.t. flow)
** Output: f = friction factor
** dfdq = derivative of f w.r.t. flow
** Purpose: computes Darcy-Weisbach friction factor and its ** Purpose: computes Darcy-Weisbach friction factor and its
** derivative as a function of Reynolds Number (Re). ** derivative as a function of Reynolds Number (Re).
**
** Note: Current formulas for dfdq need to be corrected
** so dfdq returned as 0.
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
double f; double q, // Abs. value of flow
double x1, x2, x3, x4, f; // Friction factor
y1, y2, y3, double x1,x2,x3,x4,
fa, fb, r; y1,y2,y3,
double w = q / s; // Re*Pi/4 fa,fb,r;
double s,w;
// ... For Re >= 4000 use Colebrook Formula 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) if (w >= A1)
{ {
y1 = A8 / pow(w, 0.9); y1 = A8/pow(w,0.9);
y2 = e / 3.7 + y1; y2 = link->Kc/(3.7*link->Diam) + y1;
y3 = A9 * log(y2); y3 = A9*log(y2);
f = 1.0 / (y3*y3); f = 1.0/SQR(y3);
*dfdq = 1.8 * f * y1 * A9 / y2 / y3 / q; /* *dfdq = (2.0+AA*y1/(y2*y3))*f; */ /* df/dq */
} }
// ... Use interpolating polynomials developed by // For Re > 2000 use Interpolation Formula
// E. Dunlop for transition flow from 2000 < Re < 4000. 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 else
{ {
y2 = e / 3.7 + AB; f = 8.0;
y3 = A9 * log(y2); *dfdq = 0.0;
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 = r * (0.032 - 3.0 * fa + 0.5 *fb);
f = x1 + r * (x2 + r * (x3 + x4));
*dfdq = (x2 + 2.0 * r * (x3 + x4)) / s / A2;
} }
return f; return(f);
} }

1
src/hydsolver.c
View File

@@ -981,6 +981,7 @@ void checkhydbalance(EN_Project *pr, Hydbalance *hbal)
hbal->maxheaderror = 0.0; hbal->maxheaderror = 0.0;
hbal->maxheadlink = 1; hbal->maxheadlink = 1;
for (k = 1; k <= net->Nlinks; k++) { for (k = 1; k <= net->Nlinks; k++) {
if (hyd->LinkStatus[k] <= CLOSED) continue;
hlosscoeff(pr, k); hlosscoeff(pr, k);
if (sol->P[k] == 0.0) continue; if (sol->P[k] == 0.0) continue;
link = &net->Link[k]; link = &net->Link[k];