New positional control valve added

[LRossman]

This PR is in response to the recent discussion in issue #18 . In addition to code updates that implement the new PCV valve, a new unit test for a PCV was added and the input file format section of the project's documentation (input-file.dox) was updated.

[allenmlowe]

I am excited to have this new PCV valve so I tried using a simple PCV valve without a curve yet with a setting of 0.0001 to 1 to simulate a valve opening. All the results showed a headloss of 0.0 in valve 1. Only when setting is 0.0 was there a headloss. The current GUI has no way to enter a PCV curve, so a simple opening of a valve with 1 as fully open and above 0 as partially opened is the only way for now. Maybe it needs to have a default curve if there is no curve provided for it to work.

Net1_pcv.zip

[LRossman]

@allenmlowe the loss coefficient is pegged to the PCV's fully open minor loss coefficient Ko. If no valve curve is specified then a default linear curve is used which results in a partly opened loss coefficient K equal to Ko / f^2 where f is the fraction open. So if Ko is 0, as in your example, then the valve won't produce any head loss. Try using a non-zero value for Ko like 1.0. In that case, with an opening of 0.0001, the head loss at time 0 will be 15.26. As a check, change the valve type to a TCV and use a setting of 1/0.0001^2 = 100000000. The resulting head loss will be the same as found using the PCV.

[allenmlowe]

Coming from an ordinary Epanet user's point of view, a PCV will make simulating valve openings a lot simpler. However, determining the ko is not easy but the opening which can represent the new area ratio of the valve opening is simple. So a setting of 0.9 means the area is now 0.9 of the original valve area. Even without a valve curve and ko=0, because the user does not know this, ko can be assumed equal to 1000 or a default value to be determined by experts like you. I did some experiments with different ko values and setting and plotted the new flows based on different setting and k0 as show below.

I am assuming that the new reduced area will have a corresponding reduced effective diameters and correspondingly the flow reduction in the valve.

I did some modifications in pcvlosscoeff to try to reflect the changes.

double pcvlosscoeff(Project pr, int k, double s) / -------------------------------------------------------------- Input: k = link index s = valve fraction open setting Output: returns a valve loss coefficient Purpose: finds a Positional Control Valve's loss coefficient from its fraction open setting. *-------------------------------------------------------------- / { Network* net = &pr->network;

int v = findvalve(net, k);         // valve index
int c = net->Valve[v].Curve;       // Kv curve index
double d, d1;                      // valve diameters
double A1, A2;                     // Areas of full and partial opening
double kmo;                        // fully open loss coeff.
double km;                         // partly open loss coeff.
double kvr;                        // Kv / Kvo (Kvo = Kv at fully open)
double *x, *y;                     // points on kvr v. frac. open curve
int k1, k2, npts;
Scurve *curve;

// Valve has no setting so return 0
if (s == MISSING) return 0.0;

d = net->Link[k].Diam;
kmo = net->Link[k].Km;

// Valve is completely open so return its Km value
if (s >= 1.0) return kmo;

A1 = PI * d * d / 4.0; // Original Area
A2 = A1 * s; // New Area
d1 = sqrt(A2 * 4 / PI); // New diameter

// there is no minor loss coeff.
if (kmo == 0.0)
{
   const double km = 1000; // assumed km
   kmo = 0.02517 * km / (SQR(d1 * d1));
}
else
{
   km = net->Link[k].Km * SQR(d) * SQR(d) / 0.02517; // minor headloss
   kmo = 0.02517 * km / (SQR(d1 * d1));
}

// Valve has no assigned curve so assume a linear one
if (c == 0) kvr = s;

else
{
    // Valve curve data
    curve = &net->Curve[c];
    npts = curve->Npts;
    x = curve->X;            // x = frac. open
    y = curve->Y;            // y = Kv / Kvo

    // s lies below first point of curve
    if (s < x[0])
        kvr = s / x[0] * y[0];

    // s lies above last point of curve
    else if (s > x[npts-1])
    {
        k2 = npts - 1;
        kvr = (s - x[k2]) / (1. -  x[k2]) * (1. - y[k2]) + y[k2];
    }

    // Otherwise interpolate over curve segment that brackets s
    else
    {
        k2 = 0;
        while (k2 < npts && x[k2] < s) k2++;
        if (k2 == 0) k2++;
        else if (k2 == npts)  k2--;
        k1 = k2 - 1;
        kvr = (y[k2] - y[k1]) / (x[k2] - x[k1]);
        kvr = y[k1] + kvr * (s - x[k1]);
    }
}

// kvr can't be > 1 or <= 0
kvr = MIN(kvr, 1.0);
kvr = MAX(kvr, CSMALL);

// Convert from Kv ratio to minor loss coeff.
km = kmo / (kvr * kvr);
km = MIN(km, CBIG);
return km;

}

[LRossman]

@allenmlowe your loss coefficient adjustment based on area reduction is the same thing that a linear characteristic curve does. Therefore including it in the code before the curve is applied is double counting. I don't think it's a good idea to construct an effective kmowhen the user enters 0 for the fully open loss coefficient (referred to as Koin previous comments). It makes all types of partly open valves behave the same which simply isn't the case. Values for kmo for particular types of valves (gate, ball, globe, etc) are widely available in text books, web pages and from manufacturers.

[allenmlowe]

@LRossman you are right, the GUI should remind the user to enter the loss coefficient instead of assigning a default kmo when the toolkit encounters a 0 kmo. The GUI also forces that curve values should be in ascending order. However, for this PCV curve, the normal entry should be in descending order so that the Cv/Cvo can be computed right away. Should the order of the PCV curve be in ascending order?

[LRossman]

For consistency with other curve types the points should be entered in ascending order of the abscissa variable which is fraction open.

[allenmlowe]

    // s lies below first point of curve
    if (s < x[0])
    {
        kvr = s / x[0] * y[0];

Assuming data is in ascending order: if s=0.2, x[0] = 0.25, y[0] = 24 then kvr = (0..2 / 0.25) * 24 = 19.2.
But kvr is supposed to be higher than 24 since s < 0.25. This kvr value should also be converted to Cv/Cvo since Y is loss coeff or Y is already in Cv/Cvo?

Edited: This is my fault, I pre assumed that Y data was the original loss coefficient and later converted to Cv. That assumption made me decide on the descending entry and the area reduction routine. Entering the final effective (Cv/Cvo) for the Y value is so much simpler. Sorry for the confusion. Thanks for your time @LRossman .