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

Revert "Fixes zero flow issue (#164) and DW eqn. (#199)"

Browse Source 
This reverts commit 6b808c03ab.
This commit is contained in:
Lew Rossman
2018-09-11 11:04:47 -04:00
parent 6b808c03ab
commit 836d9c3668
3 changed files with 220 additions and 268 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
src/funcs.h
View File

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

455
src/hydcoeffs.c
View File

@@ -18,24 +18,21 @@ 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 = 3.14159265358979323850e+03; // 1000*PI const double A1 = 0.314159265359e04; // 1000*PI
const double A2 = 1.57079632679489661930e+03; // 500*PI const double A2 = 0.157079632679e04; // 500*PI
const double A3 = 5.02654824574366918160e+01; // 16*PI const double A3 = 0.502654824574e02; // 16*PI
const double A4 = 6.28318530717958647700e+00; // 2*PI const double A4 = 6.283185307; // 2*PI
const double A8 = 4.61841319859066668690e+00; // 5.74*(PI/4)^.9 const double A8 = 4.61841319859; // 5.74*(PI/4)^.9
const double A9 = -8.68588963806503655300e-01; // -2/ln(10) const double A9 = -8.685889638e-01; // -2/ln(10)
const double AA = -1.5634601348517065795e+00; // -2*.9*2/ln(10) const double AA = -1.5634601348; // -2*.9*2/ln(10)
const double AB = 3.28895476345399058690e-03; // 5.74/(4000^.9) const double AB = 3.28895476345e-03; // 5.74/(4000^.9)
const double AC = -5.14214965799093883760e-03; // AA*AB const double AC = -5.14214965799e-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);
//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); //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);
@@ -50,7 +47,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(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);
@@ -198,8 +195,7 @@ void linkcoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: none ** Input: none
** Output: none ** Output: none
** Purpose: computes coefficients contributed by links to the ** Purpose: computes matrix coefficients for links
** linearized system of hydraulic equations.
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
@@ -218,34 +214,32 @@ 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 nodal flow balance (X_tmp) // Update net nodal inflows (X_tmp), solution matrix (A) and RHS array (F)
// (Flow out of node is (-), flow into node is (+)) // (Use covention that 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];
// Add to off-diagonal coeff. of linear system matrix // Off-diagonal coeff.
sol->Aij[sol->Ndx[k]] -= sol->P[k]; 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) 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]);
// Update linear system coeffs. associated with end node n2 // Node n2 is junction
// ... 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]);
} }
} }
@@ -256,8 +250,8 @@ void nodecoeffs(EN_Project *pr)
**---------------------------------------------------------------- **----------------------------------------------------------------
** Input: none ** Input: none
** Output: none ** Output: none
** Purpose: completes calculation of nodal flow balance array ** Purpose: completes calculation of nodal flow imbalance (X_tmp)
** (X_tmp) & r.h.s. (F) of linearized hydraulic eqns. ** & flow correction (F) arrays
**---------------------------------------------------------------- **----------------------------------------------------------------
*/ */
{ {
@@ -267,7 +261,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 balance & add flow balance to RHS array F // flow imbalance & add imbalance 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];
@@ -281,9 +275,8 @@ void valvecoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: none ** Input: none
** Output: none ** Output: none
** Purpose: computes coeffs. of the linearized hydraulic eqns. ** Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs
** contributed by PRVs, PSVs & FCVs whose status is ** whose status is not fixed to OPEN/CLOSED
** not fixed to OPEN/CLOSED
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
@@ -332,8 +325,7 @@ void emittercoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: none ** Input: none
** Output: none ** Output: none
** Purpose: computes coeffs. of the linearized hydraulic eqns. ** Purpose: computes matrix coeffs. for emitters
** 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
@@ -342,8 +334,12 @@ void emittercoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
int i, row; int i;
double hloss, hgrad; double ke;
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;
@@ -352,57 +348,46 @@ 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.
// Find emitter head loss and gradient q = hyd->EmitterFlows[i]; // emitter flow
emitheadloss(pr, i, &hloss, &hgrad); z = ke * pow(ABS(q), hyd->Qexp); // emitter head loss
p = hyd->Qexp * z / ABS(q); // head loss gradient
// Row of solution matrix if (p < hyd->RQtol)
row = sol->Row[i]; {
p = 1.0 / hyd->RQtol;
// Addition to matrix diagonal & r.h.s }
sol->Aii[row] += 1.0 / hgrad; else
sol->F[row] += (hloss + node->El) / hgrad; {
p = 1.0 / p; // inverse head loss gradient
// Update to node flow balance }
hyd->X_tmp[i] -= hyd->EmitterFlows[i]; 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 ** Input: i = node index
** Output: hloss = head loss across node's emitter ** Output: returns change in flow at an emitter node
** hgrad = head loss gradient ** Purpose: computes flow change at an emitter node
** Purpose: computes an emitters's head loss and gradient. **--------------------------------------------------------------
**-------------------------------------------------------------
*/ */
{ {
double ke; double ke, p;
double q;
hydraulics_t *hyd = &pr->hydraulics; hydraulics_t *hyd = &pr->hydraulics;
Snode *node = &pr->network.Node[i];
// Set adjusted emitter coeff. ke = MAX(CSMALL, node->Ke);
ke = MAX(CSMALL, pr->network.Node[i].Ke); p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]), (hyd->Qexp - 1.0));
if (p < hyd->RQtol) p = 1 / hyd->RQtol;
// Use linear head loss relation for small flow else p = 1.0 / p;
q = hyd->EmitterFlows[i]; return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El));
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;
}
} }
@@ -412,8 +397,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 ** Purpose: retrieves parameters that define a pressure dependent
** dependent demand function. ** demand function.
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
@@ -442,13 +427,12 @@ void demandcoeffs(EN_Project *pr)
**-------------------------------------------------------------- **--------------------------------------------------------------
** Input: none ** Input: none
** Output: none ** Output: none
** Purpose: computes coeffs. of the linearized hydraulic eqns. ** Purpose: computes matrix coeffs. for pressure dependent demands
** contributed by pressure dependent demands.
** **
** Note: Pressure dependent demands are modelled like emitters ** 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 ** 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 ** 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/hgrad ** P = inverse head loss gradient = 1/(dh/dQ)
** Y = flow correction term = hloss / hgrad ** Y = flow correction term = h*P
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
double hloss, // Head loss double hpipe, // Normal head loss
hgrad, // Head loss gradient hml, // Minor head loss
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: hloss = CBIG*q // For closed pipe use headloss formula: h = CBIG*q
if (hyd->LinkStatus[k] <= CLOSED) if (hyd->LinkStatus[k] <= CLOSED)
{ {
sol->P[k] = 1.0 / CBIG; sol->P[k] = 1.0 / CBIG;
@@ -589,35 +575,31 @@ void pipecoeff(EN_Project *pr, int k)
return; return;
} }
q = ABS(hyd->LinkFlows[k]); // Evaluate headloss coefficients
ml = pr->network.Link[k].Km; q = ABS(hyd->LinkFlows[k]); // Absolute flow
r = pr->network.Link[k].R; ml = link->Km; // Minor loss coeff.
r = link->R; // Resistance coeff.
// Friction head loss // Use large P coefficient for small flow resistance product
if (q <= Q_CUTOFF) if ( (r+ml)*q < hyd->RQtol)
{ {
hgrad = hyd->Hexp * r * pow(Q_CUTOFF, hyd->Hexp-1.0); sol->P[k] = 1.0 / hyd->RQtol;
hloss = q * hgrad; sol->Y[k] = hyd->LinkFlows[k] / hyd->Hexp;
} return;
else
{
hgrad = hyd->Hexp * r * pow(q, hyd->Hexp - 1.0);
hloss = hgrad * q / hyd->Hexp;
} }
// 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) if (ml > 0.0)
{ {
hloss += ml * q * q; hml = ml*q*q; // Minor head loss
hgrad += 2.0 * ml * q; p += 2.0*hml; // Q*dh(Total)/dQ
} }
else hml = 0.0;
// Adjust head loss sign for flow direction p = hyd->LinkFlows[k] / p; // 1 / (dh/dQ)
hloss *= SGN(hyd->LinkFlows[k]); sol->P[k] = ABS(p);
sol->Y[k] = p*(hpipe + hml);
// P and Y coeffs.
sol->P[k] = 1.0 / hgrad;
sol->Y[k] = hloss / hgrad;
} }
@@ -636,85 +618,102 @@ 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 r = link->R; // Resistance coeff. double dfdq = 0.0;
double ml = link->Km; // Minor loss coeff. double r, r1, f, ml, p, hloss;
double e = link->Kc / link->Diam; // Relative roughness
double s = hyd->Viscos * link->Diam; // Viscosity / diameter ml = link->Km; // Minor loss coeff.
r = link->R; // Resistance coeff.
double hloss, hgrad, f, dfdq, r1; f = frictionFactor(pr,k,&dfdq); // D-W friction factor
r1 = f*r+ml;
// Compute head loss and its derivative
// ... use Hagen-Poiseuille formula for laminar flow (Re <= 2000) // Use large P coefficient for small flow resistance product
if (q <= A2 * s) if (r1*q < hyd->RQtol)
{ {
r = 16.0 * PI * s * r; sol->P[k] = 1.0/hyd->RQtol;
hloss = hyd->LinkFlows[k] * (r + ml * q); sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp;
hgrad = r + 2.0 * ml * q; return;
} }
// ... 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
sol->P[k] = 1.0 / hgrad; hloss = r1*SQR(q); // Total head loss
sol->Y[k] = hloss / hgrad; 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| ** Input: k = link index
** e = pipe roughness / diameter ** Output: returns friction factor and
** s = viscosity * pipe diameter ** replaces dfdq (derivative of f w.r.t. flow)
** 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 f; // friction factor 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;
hydraulics_t *hyd = &pr->hydraulics;
Slink *link = &pr->network.Link[k];
// For Re >= 4000 use Swamee & Jain approximation *dfdq = 0.0;
// of the Colebrook-White Formula if (link->Type > EN_PIPE)
if ( w >= A1 ) return(1.0); // Only apply to pipes
{ q = ABS(hyd->LinkFlows[k]);
y1 = A8 / pow(w, 0.9); s = hyd->Viscos * link->Diam;
y2 = e / 3.7 + y1; w = q/s; // w = Re(Pi/4)
y3 = A9 * log(y2);
f = 1.0 / (y3*y3); // For Re >= 4000 use Colebrook Formula
*dfdq = 1.8 * f * y1 * A9 / y2 / y3 / q; 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 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;
Spump *pump; double setting = hyd->LinkSetting[k];
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;
@@ -749,54 +744,35 @@ 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 (negative sign // Determine head loss coefficients.
// 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 - 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. // Adjust head loss coefficients for pump speed.
sol->P[k] = 1.0 / hgrad; h0 = SQR(setting) * pump->H0;
sol->Y[k] = hloss / hgrad; 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, // Headloss curve intercept
double h0, // Intercept of head loss curve segment q, // Abs. value of flow
r, // Slope of head loss curve segment r; // Flow resistance coeff.
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;
@@ -860,25 +835,19 @@ 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 segment of head loss curve // Otherwise utilize headloss curve
// bracketing current flow (curve index is stored // whose index is stored in K
// in valve's setting)
else else
{ {
// Index of valve's head loss curve // Find slope & intercept of headloss 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);
// Intercept and slope of curve segment containing q // Compute inverse headloss gradient (P)
curvecoeff(pr, i, q, &h0, &r); // and flow correction factor (Y).
r = MAX(r, TINY); sol->P[k] = 1.0 / MAX(r, hyd->RQtol);
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]);
} }
} }
@@ -1115,14 +1084,14 @@ void valvecoeff(EN_Project *pr, int k)
**-------------------------------------------------------------- **--------------------------------------------------------------
*/ */
{ {
double flow, q, dhdq; double p;
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];
flow = hyd->LinkFlows[k]; double 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)
@@ -1135,26 +1104,14 @@ 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)
{ {
// Adjust for very small flow p = 2.0 * link->Km * fabs(flow);
q = MAX(fabs(flow), Q_CUTOFF); if (p < hyd->RQtol) p = hyd->RQtol;
sol->P[k] = 1.0 / p;
// Gradient of minor head loss formula sol->Y[k] = flow / 2.0;
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 else
{ {
sol->P[k] = 1.0 / CSMALL; sol->P[k] = 1.0 / hyd->RQtol;
sol->Y[k] = flow; sol->Y[k] = flow;
} }
} }

22
src/hydsolver.c
View File

@@ -494,34 +494,30 @@ void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
**---------------------------------------------------------------- **----------------------------------------------------------------
*/ */
{ {
int i; double dq;
double hloss, hgrad, dh, dq; int k;
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 (i = 1; i <= net->Njuncs; i++) for (k = 1; k <= net->Njuncs; k++)
{ {
// Skip junction if it does not have an emitter // 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 // Find emitter flow change (see hydcoeffs.c)
emitheadloss(pr, i, &hloss, &hgrad); dq = emitflowchange(pr, k);
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[i]); *qsum += ABS(hyd->EmitterFlows[k]);
*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 = i; hbal->maxflownode = k;
hbal->maxflowlink = -1; hbal->maxflowlink = -1;
} }
} }