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

Update hydcoeffs.c

Browse Source
This commit is contained in:
Lew Rossman
2021-04-20 09:30:11 -04:00
parent 9224d67a00
commit 13416942f8
1 changed files with 46 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

53
src/hydcoeffs.c
View File

@@ -7,7 +7,7 @@
Authors: see AUTHORS Authors: see AUTHORS
Copyright: see AUTHORS Copyright: see AUTHORS
License: see LICENSE License: see LICENSE
Last Updated: 10/04/2019 Last Updated: 4/20/2021
****************************************************************************** ******************************************************************************
*/ */
@@ -271,11 +271,31 @@ void nodecoeffs(Project *pr)
// For junction nodes, subtract demand flow from net // For junction nodes, subtract demand flow from net
// flow excess & add flow excess to RHS array F // flow excess & add flow excess to RHS array F
for (i = 1; i <= net->Njuncs; i++)
if (hyd->DemandModel == DDA) //ET-PDA
{ {
hyd->Xflow[i] -= hyd->DemandFlow[i]; for (i = 1; i <= net->Njuncs; i++)
sm->F[sm->Row[i]] += hyd->Xflow[i]; {
hyd->Xflow[i] -= hyd->DemandFlow[i];
sm->F[sm->Row[i]] += hyd->Xflow[i];
}
} }
else
{
for (i = 1; i <= net->Njuncs; i++)
{
hyd->DemandFlow[i] = MAX(hyd->DemandFlow[i], 0.0);
hyd->DemandFlow[i] = MIN(hyd->DemandFlow[i], hyd->NodeDemand[i]);
switch (net->Node[i].demandState)
{
case NEGATIVE_DEMAND:
case FULL_DEMAND:
hyd->DemandFlow[i] = hyd->NodeDemand[i];
hyd->Xflow[i] -= hyd->NodeDemand[i];
}
sm->F[sm->Row[i]] += hyd->Xflow[i];
}
}
} }
@@ -444,6 +464,20 @@ void demandcoeffs(Project *pr)
{ {
// Skip junctions with non-positive demands // Skip junctions with non-positive demands
if (hyd->NodeDemand[i] <= 0.0) continue; if (hyd->NodeDemand[i] <= 0.0) continue;
net->Node[i].demandState = PARTIAL_DEMAND; //ET-PDA
if (hyd->DemandFlow[i] <= 0 &&
hyd->NodeHead[i] <= (net->Node[i].El + hyd->Pmin))
{
net->Node[i].demandState = ZERO_DEMAND;
continue;
}
else if (MAX(0,hyd->DemandFlow[i]) * hyd->NodeHead[i] >=
hyd->NodeDemand[i] * (net->Node[i].El + hyd->Preq))
{
net->Node[i].demandState = FULL_DEMAND;
continue;
}
// Find head loss for demand outflow at node's elevation // Find head loss for demand outflow at node's elevation
demandheadloss(pr, i, dp, n, &hloss, &hgrad); demandheadloss(pr, i, dp, n, &hloss, &hgrad);
@@ -474,10 +508,14 @@ void demandheadloss(Project *pr, int i, double dp, double n,
{ {
Hydraul *hyd = &pr->hydraul; Hydraul *hyd = &pr->hydraul;
double d = hyd->DemandFlow[i]; double d = MAX(hyd->DemandFlow[i], 0.0001); //ET-PDA
double dfull = hyd->NodeDemand[i]; double dfull = hyd->NodeDemand[i];
double r = d / dfull; double r = d / dfull;
*hgrad = n * dp * pow(r, n - 1.0) / dfull;
*hloss = (*hgrad) * d / n;
/* //ET-PDA
// Use lower barrier function for negative demand // Use lower barrier function for negative demand
if (r <= 0) if (r <= 0)
{ {
@@ -486,7 +524,7 @@ void demandheadloss(Project *pr, int i, double dp, double n,
} }
// Use power head loss function for demand less than full // Use power head loss function for demand less than full
else if (r <= 1.0) else if (r <= 1.0) //(2.2.1)
{ {
*hgrad = n * dp * pow(r, n - 1.0) / dfull; *hgrad = n * dp * pow(r, n - 1.0) / dfull;
// ... use linear function for very small gradient // ... use linear function for very small gradient
@@ -504,6 +542,7 @@ void demandheadloss(Project *pr, int i, double dp, double n,
*hgrad = CBIG; *hgrad = CBIG;
*hloss = dp + CBIG * (d - dfull); *hloss = dp + CBIG * (d - dfull);
} }
*/
} }