Merge pull request #5 from LRossman/leakage

Adding Pipe Leakage Modeling

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,25 +113,25 @@ 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;
             else
                 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;

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);

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;
 }
 
 

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,30 +373,30 @@ 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
     // current leakage flow rates
     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;
     }
 
     // GGA flow update formula for variable area leakage
     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,72 +462,77 @@ 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)
     {
         *hfa = 0.0;
         *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.
 **--------------------------------------------------------------------
 */
 {
     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);
 }