Update hydsolver.c

This commit is contained in:
Lew Rossman
2022-02-14 11:03:22 -05:00
parent 650d8e8e53
commit 070901d6bf

View File

@@ -555,8 +555,8 @@ void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum)
dq *= hyd->RelaxFactor;
// Prevent a flow change greater than full demand
if (fabs(dq) > hyd->NodeDemand[i])
dq = 0.5 * SGN(dq) * hyd->NodeDemand[i];
if (fabs(dq) > 0.4 * hyd->NodeDemand[i])
dq = 0.4 * SGN(dq) * hyd->NodeDemand[i];
hyd->DemandFlow[i] -= dq;
// Update system flow summation
@@ -658,9 +658,11 @@ int pdaconverged(Project *pr)
{
Hydraul *hyd = &pr->hydraul;
const double TOL = 0.001;
const double QTOL = 0.0001; // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm)
int i, converged = 1;
double totalDemand = 0.0, totalReduction = 0.0;
double dp = hyd->Preq - hyd->Pmin;
double p, q, r;
hyd->DeficientNodes = 0;
hyd->DemandReduction = 0.0;
@@ -671,14 +673,24 @@ int pdaconverged(Project *pr)
// Skip nodes whose required demand is non-positive
if (hyd->NodeDemand[i] <= 0.0) continue;
// Check for negative demand flow or positive demand flow at negative pressure
if (hyd->DemandFlow[i] < -TOL) converged = 0;
if (hyd->DemandFlow[i] > TOL &&
hyd->NodeHead[i] - pr->network.Node[i].El - hyd->Pmin < -TOL)
converged = 0;
// Evaluate demand equation at current pressure solution
p = hyd->NodeHead[i] - pr->network.Node[i].El;
if (p <= hyd->Pmin)
q = 0.0;
else if (p >= hyd->Preq)
q = hyd->NodeDemand[i];
else
{
r = (p - hyd->Pmin) / dp;
q = hyd->NodeDemand[i] * pow(r, hyd->Pexp);
}
// Check if demand has not converged
if (fabs(q - hyd->DemandFlow[i]) > QTOL)
converged = 0;
// Accumulate total required demand and demand deficit
if (hyd->DemandFlow[i] + 0.0001 < hyd->NodeDemand[i])
if (hyd->DemandFlow[i] + QTOL < hyd->NodeDemand[i])
{
hyd->DeficientNodes++;
totalDemand += hyd->NodeDemand[i];