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

Updates to fix zero flow (#164) and D-W eqn. (#199) issues

Browse Source
This commit is contained in:
Lew Rossman
2018-09-11 13:15:15 -04:00
parent 836d9c3668
commit dab7be8446
3 changed files with 274 additions and 221 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
src/funcs.h
View File

@@ -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 resistcoeff(EN_Project *pr, int k); /* Finds pipe flow resistance */
void headlosscoeffs(EN_Project *pr); // Finds link head loss coeffs. void headlosscoeffs(EN_Project *pr); // Finds link head loss coeffs.
void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. */ void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. */
double emitflowchange(EN_Project *pr, int i); /* Change in emitter outflow */ void emitheadloss(EN_Project *pr, int, // Finds emitter head loss
double demandflowchange(EN_Project *pr, int i, // Change in demand outflow double *, double *);
double dp, double n); double demandflowchange(EN_Project *pr, int, // Change in demand outflow
void demandparams(EN_Project *pr, double *dp, // PDA function parameters double, double);
double *n); void demandparams(EN_Project *pr, double *, // PDA function parameters
double *);
/* ----------- SMATRIX.C ---------------*/ /* ----------- SMATRIX.C ---------------*/
int createsparse(EN_Project *pr); /* Creates sparse matrix */ int createsparse(EN_Project *pr); /* Creates sparse matrix */

462
src/hydcoeffs.c
View File

@@ -18,21 +18,24 @@ HYDCOEFFS.C -- hydraulic coefficients for the EPANET Program
#include "funcs.h" #include "funcs.h"
// Constants used for computing Darcy-Weisbach friction factor // Constants used for computing Darcy-Weisbach friction factor
const double A1 = 0.314159265359e04; // 1000*PI const double A1 = 3.14159265358979323850e+03; // 1000*PI
const double A2 = 0.157079632679e04; // 500*PI const double A2 = 1.57079632679489661930e+03; // 500*PI
const double A3 = 0.502654824574e02; // 16*PI const double A3 = 5.02654824574366918160e+01; // 16*PI
const double A4 = 6.283185307; // 2*PI const double A4 = 6.28318530717958647700e+00; // 2*PI
const double A8 = 4.61841319859; // 5.74*(PI/4)^.9 const double A8 = 4.61841319859066668690e+00; // 5.74*(PI/4)^.9
const double A9 = -8.685889638e-01; // -2/ln(10) const double A9 = -8.68588963806503655300e-01; // -2/ln(10)
const double AA = -1.5634601348; // -2*.9*2/ln(10) const double AA = -1.5634601348517065795e+00; // -2*.9*2/ln(10)
const double AB = 3.28895476345e-03; // 5.74/(4000^.9) const double AB = 3.28895476345399058690e-03; // 5.74/(4000^.9)
const double AC = -5.14214965799e-03; // AA*AB 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 // External functions
//void resistcoeff(EN_Project *pr, int k); //void resistcoeff(EN_Project *pr, int k);
//void headlosscoeffs(EN_Project *pr); //void headlosscoeffs(EN_Project *pr);
//void matrixcoeffs(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); //double demandflowchange(EN_Project *pr, int i, double dp, double n);
//void demandparams(EN_Project *pr, 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 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(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 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);
@@ -195,7 +198,8 @@ void linkcoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: none ** Input: none
** Output: 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 n1 = link->N1; // Start node of link
n2 = link->N2; // End node of link n2 = link->N2; // End node of link
// Update net nodal inflows (X_tmp), solution matrix (A) and RHS array (F) // Update nodal flow balance (X_tmp)
// (Use covention that flow out of node is (-), flow into node is (+)) // (Flow out of node is (-), flow into node is (+))
hyd->X_tmp[n1] -= hyd->LinkFlows[k]; hyd->X_tmp[n1] -= hyd->LinkFlows[k];
hyd->X_tmp[n2] += 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]; 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) if (n1 <= net->Njuncs)
{ {
sol->Aii[sol->Row[n1]] += sol->P[k]; // Diagonal coeff. sol->Aii[sol->Row[n1]] += sol->P[k]; // Diagonal coeff.
sol->F[sol->Row[n1]] += sol->Y[k]; // RHS 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]); 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) if (n2 <= net->Njuncs)
{ {
sol->Aii[sol->Row[n2]] += sol->P[k]; // Diagonal coeff. sol->Aii[sol->Row[n2]] += sol->P[k]; // Diagonal coeff.
sol->F[sol->Row[n2]] -= sol->Y[k]; // RHS 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]); else sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]);
} }
} }
@@ -250,8 +256,8 @@ void nodecoeffs(EN_Project *pr)
**---------------------------------------------------------------- **----------------------------------------------------------------
** Input: none ** Input: none
** Output: none ** Output: none
** Purpose: completes calculation of nodal flow imbalance (X_tmp) ** Purpose: completes calculation of nodal flow balance array
** & flow correction (F) arrays ** (X_tmp) & r.h.s. (F) of linearized hydraulic eqns.
**---------------------------------------------------------------- **----------------------------------------------------------------
*/ */
{ {
@@ -261,7 +267,7 @@ void nodecoeffs(EN_Project *pr)
EN_Network *net = &pr->network; EN_Network *net = &pr->network;
// For junction nodes, subtract demand flow from net // 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++) for (i = 1; i <= net->Njuncs; i++)
{ {
hyd->X_tmp[i] -= hyd->DemandFlows[i]; hyd->X_tmp[i] -= hyd->DemandFlows[i];
@@ -275,8 +281,9 @@ void valvecoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: none ** Input: none
** Output: none ** Output: none
** Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs ** Purpose: computes coeffs. of the linearized hydraulic eqns.
** whose status is not fixed to OPEN/CLOSED ** contributed by PRVs, PSVs & FCVs whose status is
** not fixed to OPEN/CLOSED
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
@@ -325,7 +332,8 @@ void emittercoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: none ** Input: none
** Output: 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 ** Note: Emitters consist of a fictitious pipe connected to
** a fictitious reservoir whose elevation equals that ** a fictitious reservoir whose elevation equals that
@@ -334,12 +342,8 @@ void emittercoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
int i; int i, row;
double ke; double hloss, hgrad;
double p;
double q;
double y;
double z;
hydraulics_t *hyd = &pr->hydraulics; hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver; solver_t *sol = &hyd->solver;
@@ -348,46 +352,57 @@ void emittercoeffs(EN_Project *pr)
for (i = 1; i <= net->Njuncs; i++) for (i = 1; i <= net->Njuncs; i++)
{ {
// Skip junctions without emitters
node = &net->Node[i]; node = &net->Node[i];
if (node->Ke == 0.0) continue; if (node->Ke == 0.0) continue;
ke = MAX(CSMALL, node->Ke); // emitter coeff.
q = hyd->EmitterFlows[i]; // emitter flow // Find emitter head loss and gradient
z = ke * pow(ABS(q), hyd->Qexp); // emitter head loss emitheadloss(pr, i, &hloss, &hgrad);
p = hyd->Qexp * z / ABS(q); // head loss gradient
if (p < hyd->RQtol) // Row of solution matrix
{ row = sol->Row[i];
p = 1.0 / hyd->RQtol;
} // Addition to matrix diagonal & r.h.s
else sol->Aii[row] += 1.0 / hgrad;
{ sol->F[row] += (hloss + node->El) / hgrad;
p = 1.0 / p; // inverse head loss gradient
} // Update to node flow balance
y = SGN(q)*z*p; // head loss / gradient hyd->X_tmp[i] -= hyd->EmitterFlows[i];
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
} }
} }
double emitflowchange(EN_Project *pr, int i) void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad)
/* /*
**-------------------------------------------------------------- **-------------------------------------------------------------
** Input: i = node index ** Input: i = node index
** Output: returns change in flow at an emitter node ** Output: hloss = head loss across node's emitter
** Purpose: computes flow change at an emitter node ** 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; hydraulics_t *hyd = &pr->hydraulics;
Snode *node = &pr->network.Node[i];
ke = MAX(CSMALL, node->Ke); // Set adjusted emitter coeff.
p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]), (hyd->Qexp - 1.0)); ke = MAX(CSMALL, pr->network.Node[i].Ke);
if (p < hyd->RQtol) p = 1 / hyd->RQtol;
else p = 1.0 / p; // Use linear head loss relation for small flow
return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El)); 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 ** Input: none
** Output: dp = pressure range over which demands can vary ** Output: dp = pressure range over which demands can vary
** n = exponent in head loss v. demand function ** n = exponent in head loss v. demand function
** Purpose: retrieves parameters that define a pressure dependent ** Purpose: retrieves parameters that define a pressure
** demand function. ** dependent demand function.
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
@@ -427,12 +442,13 @@ void demandcoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: none ** Input: none
** Output: 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 ** 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 ** 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 // Use linear head loss function for near zero demand
else if (r < EPS) else if (r < EPS)
{ {
*hgrad = dp * pow(EPS, n - 1.0) / dfull; *hgrad = dp * pow(EPS, n) / dfull / EPS;
*hloss = (*hgrad) * d; *hloss = (*hgrad) * d;
} }
@@ -542,25 +558,23 @@ void pipecoeff(EN_Project *pr, int k)
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: k = link index ** Input: k = link index
** Output: none ** 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) ** P = inverse head loss gradient = 1/hgrad
** Y = flow correction term = h*P ** Y = flow correction term = hloss / hgrad
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
double hpipe, // Normal head loss double hloss, // Head loss
hml, // Minor head loss hgrad, // Head loss gradient
ml, // Minor loss coeff. ml, // Minor loss coeff.
p, // q*(dh/dq)
q, // Abs. value of flow q, // Abs. value of flow
r; // Resistance coeff. r; // Resistance coeff.
hydraulics_t *hyd = &pr->hydraulics; hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver; 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) if (hyd->LinkStatus[k] <= CLOSED)
{ {
sol->P[k] = 1.0 / CBIG; sol->P[k] = 1.0 / CBIG;
@@ -575,31 +589,37 @@ void pipecoeff(EN_Project *pr, int k)
return; return;
} }
// Evaluate headloss coefficients q = ABS(hyd->LinkFlows[k]);
q = ABS(hyd->LinkFlows[k]); // Absolute flow ml = pr->network.Link[k].Km;
ml = link->Km; // Minor loss coeff. r = pr->network.Link[k].R;
r = link->R; // Resistance coeff.
// Use large P coefficient for small flow resistance product // Friction head loss
if ( (r+ml)*q < hyd->RQtol) // ... use linear relation for small flows
if (q <= Q_CUTOFF)
{ {
sol->P[k] = 1.0 / hyd->RQtol; hgrad = r * pow(Q_CUTOFF, hyd->Hexp) / Q_CUTOFF;
sol->Y[k] = hyd->LinkFlows[k] / hyd->Hexp; hloss = hgrad * q;
return; }
// ... 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 // Contribution of minor head loss
hpipe = r*pow(q, hyd->Hexp); // Friction head loss
p = hyd->Hexp*hpipe; // Q*dh(friction)/dQ
if (ml > 0.0) if (ml > 0.0)
{ {
hml = ml*q*q; // Minor head loss hloss += ml * q * q;
p += 2.0*hml; // Q*dh(Total)/dQ hgrad += 2.0 * ml * q;
} }
else hml = 0.0;
p = hyd->LinkFlows[k] / p; // 1 / (dh/dQ) // Adjust head loss sign for flow direction
sol->P[k] = ABS(p); hloss *= SGN(hyd->LinkFlows[k]);
sol->Y[k] = p*(hpipe + hml);
// 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]; Slink *link = &pr->network.Link[k];
double q = ABS(hyd->LinkFlows[k]); double q = ABS(hyd->LinkFlows[k]);
double dfdq = 0.0; double r = link->R; // Resistance coeff.
double r, r1, f, ml, p, hloss; double ml = link->Km; // Minor loss coeff.
double e = link->Kc / link->Diam; // Relative roughness
ml = link->Km; // Minor loss coeff. double s = hyd->Viscos * link->Diam; // Viscosity / diameter
r = link->R; // Resistance coeff.
f = frictionFactor(pr,k,&dfdq); // D-W friction factor double hloss, hgrad, f, dfdq, r1;
r1 = f*r+ml;
// Compute head loss and its derivative
// Use large P coefficient for small flow resistance product // ... use Hagen-Poiseuille formula for laminar flow (Re <= 2000)
if (r1*q < hyd->RQtol) if (q <= A2 * s)
{ {
sol->P[k] = 1.0/hyd->RQtol; r = 16.0 * PI * s * r;
sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp; hloss = hyd->LinkFlows[k] * (r + ml * q);
return; 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 // Compute P and Y coefficients
hloss = r1*SQR(q); // Total head loss sol->P[k] = 1.0 / hgrad;
p = 2.0*r1*q; // |dHloss/dQ| sol->Y[k] = hloss / hgrad;
// + 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(EN_Project *pr, int k, double *dfdq) double frictionFactor(double q, double e, double s, double *dfdq)
/* /*
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: k = link index ** Input: q = |pipe flow|
** Output: returns friction factor and ** e = pipe roughness / diameter
** replaces dfdq (derivative of f w.r.t. flow) ** 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 ** 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 q, // Abs. value of flow double f; // friction factor
f; // Friction factor double x1, x2, x3, x4,
double x1,x2,x3,x4, y1, y2, y3,
y1,y2,y3, fa, fb, r;
fa,fb,r; double w = q / s; // Re*Pi/4
double s,w;
hydraulics_t *hyd = &pr->hydraulics;
Slink *link = &pr->network.Link[k];
*dfdq = 0.0; // For Re >= 4000 use Swamee & Jain approximation
if (link->Type > EN_PIPE) // of the Colebrook-White Formula
return(1.0); // Only apply to pipes if ( w >= A1 )
q = ABS(hyd->LinkFlows[k]); {
s = hyd->Viscos * link->Diam; y1 = A8 / pow(w, 0.9);
w = q/s; // w = Re(Pi/4) y2 = e / 3.7 + y1;
y3 = A9 * log(y2);
// For Re >= 4000 use Colebrook Formula f = 1.0 / (y3*y3);
if (w >= A1) *dfdq = 1.8 * f * y1 * A9 / y2 / y3 / q;
{ }
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;
} }
@@ -730,13 +733,17 @@ void pumpcoeff(EN_Project *pr, int k)
double h0, // Shutoff head double h0, // Shutoff head
q, // Abs. value of flow q, // Abs. value of flow
r, // Flow resistance coeff. 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; hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver; 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 // Use high resistance pipe if pump closed or cannot deliver head
setting = hyd->LinkSetting[k];
if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0)
{ {
sol->P[k] = 1.0 / CBIG; sol->P[k] = 1.0 / CBIG;
@@ -744,35 +751,54 @@ void pumpcoeff(EN_Project *pr, int k)
return; return;
} }
// Obtain reference to pump object & its speed setting
q = ABS(hyd->LinkFlows[k]); q = ABS(hyd->LinkFlows[k]);
q = MAX(q, TINY);
// Obtain reference to pump object
p = findpump(&pr->network, k); p = findpump(&pr->network, k);
pump = &pr->network.Pump[p]; 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) if (pump->Ptype == CUSTOM)
{ {
// Find intercept (h0) & slope (r) of pump curve // Find intercept (h0) & slope (r) of pump curve
// line segment which contains speed-adjusted flow. // line segment which contains speed-adjusted flow.
curvecoeff(pr, pump->Hcurve, q / setting, &h0, &r); 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->H0 = -h0;
pump->R = -r; pump->R = -r;
pump->N = 1.0; 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) / 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;
}
} }
// Adjust head loss coefficients for pump speed. // P and Y coeffs.
h0 = SQR(setting) * pump->H0; sol->P[k] = 1.0 / hgrad;
n = pump->N; sol->Y[k] = hloss / hgrad;
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;
} }
@@ -825,9 +851,10 @@ void gpvcoeff(EN_Project *pr, int k)
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
double h0, // Headloss curve intercept int i;
q, // Abs. value of flow double h0, // Intercept of head loss curve segment
r; // Flow resistance coeff. r, // Slope of head loss curve segment
q; // Abs. value of flow
hydraulics_t *hyd = &pr->hydraulics; hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver; solver_t *sol = &hyd->solver;
@@ -835,19 +862,25 @@ void gpvcoeff(EN_Project *pr, int k)
// Treat as a pipe if valve closed // Treat as a pipe if valve closed
if (hyd->LinkStatus[k] == CLOSED) valvecoeff(pr, k); if (hyd->LinkStatus[k] == CLOSED) valvecoeff(pr, k);
// Otherwise utilize headloss curve // Otherwise utilize segment of head loss curve
// whose index is stored in K // bracketing current flow (curve index is stored
// in valve's setting)
else 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 = ABS(hyd->LinkFlows[k]);
q = MAX(q, TINY); q = MAX(q, TINY);
curvecoeff(pr, (int)ROUND(hyd->LinkSetting[k]), q, &h0, &r);
// Compute inverse headloss gradient (P) // Intercept and slope of curve segment containing q
// and flow correction factor (Y). curvecoeff(pr, i, q, &h0, &r);
sol->P[k] = 1.0 / MAX(r, hyd->RQtol); r = MAX(r, TINY);
sol->Y[k] = sol->P[k] * (h0 + r*q) * SGN(hyd->LinkFlows[k]);
// 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; EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics; hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver; solver_t *sol = &hyd->solver;
Slink *link = &net->Link[k]; Slink *link = &net->Link[k];
double flow = hyd->LinkFlows[k]; flow = hyd->LinkFlows[k];
// Valve is closed. Use a very small matrix coeff. // Valve is closed. Use a very small matrix coeff.
if (hyd->LinkStatus[k] <= CLOSED) if (hyd->LinkStatus[k] <= CLOSED)
@@ -1104,14 +1137,29 @@ void valvecoeff(EN_Project *pr, int k)
// Account for any minor headloss through the valve // Account for any minor headloss through the valve
if (link->Km > 0.0) if (link->Km > 0.0)
{ {
p = 2.0 * link->Km * fabs(flow); // Adjust for very small flow
if (p < hyd->RQtol) p = hyd->RQtol; q = fabs(flow);
sol->P[k] = 1.0 / p; if (q <= Q_CUTOFF)
sol->Y[k] = flow / 2.0; {
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 else
{ {
sol->P[k] = 1.0 / hyd->RQtol; sol->P[k] = 1.0 / CSMALL;
sol->Y[k] = flow; sol->Y[k] = flow;
} }
} }

22
src/hydsolver.c
View File

@@ -494,30 +494,34 @@ void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
**---------------------------------------------------------------- **----------------------------------------------------------------
*/ */
{ {
double dq; int i;
int k; double hloss, hgrad, dh, dq;
EN_Network *net = &pr->network; EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics; hydraulics_t *hyd = &pr->hydraulics;
// Examine each network junction // 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 // 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) // Find emitter head loss and gradient
dq = emitflowchange(pr, k); emitheadloss(pr, i, &hloss, &hgrad);
hyd->EmitterFlows[k] -= dq;
// Find emitter flow change
dh = hyd->NodeHead[i] - net->Node[i].El;
dq = (hloss - dh) / hgrad;
hyd->EmitterFlows[i] -= dq;
// Update system flow summation // Update system flow summation
*qsum += ABS(hyd->EmitterFlows[k]); *qsum += ABS(hyd->EmitterFlows[i]);
*dqsum += ABS(dq); *dqsum += ABS(dq);
// Update identity of element with max. flow change // Update identity of element with max. flow change
if (ABS(dq) > hbal->maxflowchange) if (ABS(dq) > hbal->maxflowchange)
{ {
hbal->maxflowchange = ABS(dq); hbal->maxflowchange = ABS(dq);
hbal->maxflownode = k; hbal->maxflownode = i;
hbal->maxflowlink = -1; hbal->maxflowlink = -1;
} }
} }