diff --git a/src/flowbalance.c b/src/flowbalance.c index 2c0efa66..aea6f199 100644 --- a/src/flowbalance.c +++ b/src/flowbalance.c @@ -73,9 +73,7 @@ void updateflowbalance(Project *pr, long hstep) flowBalance.storageDemand = 0.0; fullDemand = 0.0; - // Initialize demand deficiency & leakage loss - hyd->DeficientNodes = 0; - hyd->DemandReduction = 0.0; + // Initialize leakage loss hyd->LeakageLoss = 0.0; // Examine each junction node @@ -104,11 +102,8 @@ void updateflowbalance(Project *pr, long hstep) if (hyd->DemandModel == PDA && hyd->FullDemand[i] > 0.0) { deficit = hyd->FullDemand[i] - hyd->DemandFlow[i]; - if (deficit > TINY) - { - hyd->DeficientNodes++; + if (deficit > 0.0) flowBalance.deficitDemand += deficit; - } } } @@ -118,8 +113,8 @@ void updateflowbalance(Project *pr, long hstep) i = net->Tank[j].Node; v = hyd->NodeDemand[i]; - // For a snapshot analysis or a reservoir node - if (time->Dur == 0 || net->Tank[j].A == 0.0) + // For a reservoir node + if (net->Tank[j].A == 0.0) { if (v >= 0.0) flowBalance.totalOutflow += v; @@ -127,16 +122,16 @@ void updateflowbalance(Project *pr, long hstep) flowBalance.totalInflow += (-v); } - // For tank under extended period analysis + // For tank else flowBalance.storageDemand += v; } - // Find % demand reduction & % leakage for current period - if (fullDemand > 0.0) - hyd->DemandReduction = flowBalance.deficitDemand / fullDemand * 100.0; - if (flowBalance.totalInflow > 0.0) - hyd->LeakageLoss = flowBalance.leakageDemand / flowBalance.totalInflow * 100.0; + // Find % leakage for current period + v = flowBalance.totalInflow; + if (flowBalance.storageDemand < 0.0) v += (-flowBalance.storageDemand); + if (v > 0.0) + hyd->LeakageLoss = flowBalance.leakageDemand / v * 100.0; // Update flow balance for entire run hyd->FlowBalance.totalInflow += flowBalance.totalInflow * dt; diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index 5be8c403..670ff1ba 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -575,7 +575,7 @@ void demandcoeffs(Project *pr) for (i = 1; i <= net->Njuncs; i++) { // Skip junctions with non-positive demands - if (hyd->NodeDemand[i] <= 0.0) continue; + if (hyd->FullDemand[i] <= 0.0) continue; // Find head loss for demand outflow at node's elevation demandheadloss(pr, i, dp, n, &hloss, &hgrad); diff --git a/src/hydsolver.c b/src/hydsolver.c index d88823fb..5d76bc60 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -8,7 +8,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 06/15/2024 + Last Updated: 06/26/2024 ****************************************************************************** */ @@ -703,10 +703,15 @@ int pdaconverged(Project *pr) Hydraul *hyd = &pr->hydraul; const double QTOL = 0.0001; // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm) - int i; + 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; + // Examine each network junction for (i = 1; i <= pr->network.Njuncs; i++) { @@ -726,9 +731,20 @@ int pdaconverged(Project *pr) } // Check if demand has not converged - if (fabs(q - hyd->DemandFlow[i]) > QTOL) return 0; + if (fabs(q - hyd->DemandFlow[i]) > QTOL) + converged = 0; + + // Accumulate demand deficient node count and demand deficit + if (hyd->DemandFlow[i] + QTOL < hyd->FullDemand[i]) + { + hyd->DeficientNodes++; + totalDemand += hyd->FullDemand[i]; + totalReduction += hyd->FullDemand[i] - hyd->DemandFlow[i]; + } } - return 1; + if (totalDemand > 0.0) + hyd->DemandReduction = totalReduction / totalDemand * 100.0; + return converged; } diff --git a/src/leakage.c b/src/leakage.c index cb5cad71..a395053e 100644 --- a/src/leakage.c +++ b/src/leakage.c @@ -16,7 +16,7 @@ leaky pipes: Q = Co * L * (Ao + m * H) * sqrt(H) -where Q = leak flow rate, Co = an orifice coefficient (= 0.6*sqrt(2g)), +where Q = pipe leak flow rate, Co = an orifice coefficient (= 0.6*sqrt(2g)), L = pipe length, Ao = initial area of leak per unit of pipe length, m = change in leak area per unit of pressure head, and H = pressure head. @@ -26,7 +26,7 @@ a pipe's end node using a pair of equivalent emitters as follows: H = Cfa * Qfa^2 H = Cva * Qva^(2/3) -where Qfa = fixed area leakage rate, Qva = variable area leakage rate, +where Qfa = fixed area node leakage rate, Qva = variable area node leakage rate, Cfa = 1 / SUM(Co*(L/2)*Ao)^2, Cva = 1 / SUM(Co*(L/2)*m)^2/3, and SUM(x) is the summation of x over all pipes connected to the node. @@ -56,9 +56,9 @@ static void convert_pipe_to_node_leakage(Project *pr); static void init_node_leakage(Project *pr); static int leakage_headloss(Project* pr, int i, double *hfa, double *gfa, double *hva, double *gva); -static void eval_node_leakage(double RQtol, double q, double c, - double n, double *h, double *g); -static void add_lower_barrier(double q, double* h, double* g); +static void eval_leak_headloss(double RQtol, double q, double c, + double n, double *hloss, double *hgrad); +static void add_lower_barrier(double q, double *hloss, double *hgrad); int openleakage(Project *pr) @@ -116,7 +116,7 @@ int create_leakage_objects(Project *pr) /*------------------------------------------------------------- ** Input: none ** Output: returns an error code -** Purpose: allocates an array of Leakage objects. +** Purpose: allocates an array of node leakage objects. **------------------------------------------------------------- */ { @@ -153,9 +153,11 @@ void convert_pipe_to_node_leakage(Project *pr) Slink *link; Snode *node1; Snode *node2; + + // Orifice coeff. with conversion from sq. mm to sq. m + c_orif = 4.8149866 * 1.e-6; // Examine each link - c_orif = 4.8149866 * 1.e-6; for (i = 1; i <= net->Nlinks; i++) { // Only pipes have leakage @@ -371,8 +373,8 @@ double leakageflowchange(Project *pr, int i) Network *net = &pr->network; Hydraul *hyd = &pr->hydraul; - double hfa, gfa, hva, gva, // same as defined in leakage_solvercoeffs() - dh, dqfa, dqva; + double hfa, gfa, hva, gva; // same as defined in leakage_solvercoeffs() + double h, dqfa, dqva; // pressure head, change in leakage flows // Find the head loss and gradient of the inverted leakage // equation for both fixed and variable area leakage at the @@ -380,13 +382,13 @@ double leakageflowchange(Project *pr, int i) if (!leakage_headloss(pr, i, &hfa, &gfa, &hva, &gva)) return 0.0; // Pressure head using latest head solution - dh = hyd->NodeHead[i] - net->Node[i].El; + h = hyd->NodeHead[i] - net->Node[i].El; // GGA flow update formula for fixed area leakage dqfa = 0.0; if (gfa > 0.0) { - dqfa = (hfa - dh) / gfa * hyd->RelaxFactor; + dqfa = (hfa - h) / gfa * hyd->RelaxFactor; hyd->Leakage[i].qfa -= dqfa; } @@ -394,7 +396,7 @@ double leakageflowchange(Project *pr, int i) dqva = 0.0; if (gva > 0.0) { - dqva = (hva - dh) / gva * hyd->RelaxFactor; + dqva = (hva - h) / gva * hyd->RelaxFactor; hyd->Leakage[i].qva -= dqva; } @@ -415,10 +417,10 @@ int leakagehasconverged(Project *pr) { Network *net = &pr->network; Hydraul *hyd = &pr->hydraul; + int i; double h, qref, qtest; - const double ABSTOL = 0.0001; // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm) - const double RELTOL = 0.001; + const double QTOL = 0.0001; // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm) for (i = 1; i <= net->Njuncs; i++) { @@ -439,7 +441,7 @@ int leakagehasconverged(Project *pr) // Compare reference leakage to solution leakage qtest = hyd->Leakage[i].qfa + hyd->Leakage[i].qva; - if (fabs(qref - qtest) > ABSTOL + RELTOL * qref) return FALSE; + if (fabs(qref - qtest) > QTOL) return FALSE; } return TRUE; } @@ -460,6 +462,7 @@ int leakage_headloss(Project* pr, int i, double *hfa, double *gfa, */ { Hydraul *hyd = &pr->hydraul; + if (hyd->Leakage[i].cfa == 0.0 && hyd->Leakage[i].cva == 0.0) return FALSE; if (hyd->Leakage[i].cfa == 0.0) { @@ -467,58 +470,62 @@ int leakage_headloss(Project* pr, int i, double *hfa, double *gfa, *gfa = 0.0; } else - eval_node_leakage(hyd->RQtol, hyd->Leakage[i].qfa, hyd->Leakage[i].cfa, - 0.5, hfa, gfa); + eval_leak_headloss(hyd->RQtol, hyd->Leakage[i].qfa, hyd->Leakage[i].cfa, + 0.5, hfa, gfa); if (hyd->Leakage[i].cva == 0.0) { *hva = 0.0; *gva = 0.0; } else - eval_node_leakage(hyd->RQtol, hyd->Leakage[i].qva, hyd->Leakage[i].cva, - 1.5, hva, gva); + eval_leak_headloss(hyd->RQtol, hyd->Leakage[i].qva, hyd->Leakage[i].cva, + 1.5, hva, gva); return TRUE; } -void eval_node_leakage(double RQtol, double q, double c, double n, - double *h, double *g) +void eval_leak_headloss(double RQtol, double q, double c, double n, + double *hloss, double *hgrad) /* **-------------------------------------------------------------- ** Input: RQtol = low gradient tolerance (ft/cfs) ** q = leakage flow rate (cfs) ** c = leakage head loss coefficient ** n = leakage head loss exponent -** Output: h = leakage head loss (ft) -** g = gradient of leakage head loss (ft/cfs) +** Output: hloss = leakage head loss (ft) +** hgrad = gradient of leakage head loss (ft/cfs) ** Purpose: evaluates inverted form of leakage equation to ** compute head loss and its gradient as a function ** flow. +** +** Note: Inverted leakage equation is: +** hloss = c * q ^ (1/n) **-------------------------------------------------------------- */ { n = 1.0 / n; - *g = n * c * pow(fabs(q), n - 1.0); + *hgrad = n * c * pow(fabs(q), n - 1.0); // Use linear head loss function for small gradient - if (*g < RQtol) +/* if (*hgrad < RQtol) { - *g = RQtol / n; - *h = (*g) * q; + *hgrad = RQtol / n; + *hloss = (*hgrad) * q; } // Otherwise use normal leakage head loss function - else *h = (*g) * q / n; + else */ + *hloss = (*hgrad) * q / n; // Prevent leakage from going negative - add_lower_barrier(q, h, g); + add_lower_barrier(q, hloss, hgrad); } -void add_lower_barrier(double q, double* h, double* g) +void add_lower_barrier(double q, double* hloss, double* hgrad) /* **-------------------------------------------------------------------- ** Input: q = current flow rate -** Output: h = head loss value -** g = head loss gradient value +** Output: hloss = head loss value +** hgrad = head loss gradient value ** Purpose: adds a head loss barrier to prevent flow from falling ** below 0. **-------------------------------------------------------------------- @@ -526,6 +533,6 @@ void add_lower_barrier(double q, double* h, double* g) { double a = 1.e9 * q; double b = sqrt(a*a + 1.e-6); - *h += (a - b) / 2.; - *g += (1.e9 / 2.) * ( 1.0 - a / b); + *hloss += (a - b) / 2.; + *hgrad += (1.e9 / 2.) * ( 1.0 - a / b); }