FAVAD implementation always uses US units for “m” (EN_LEAK_EXPAN)

[dapaez270]

Not sure if this is by design or if there is a missing unit conversion in leakage.c, but to replicate the FAVAD simulation results in SI units (e.g., LPS or MLD), the value for the leak expansion rate due to pressure EN_LEAK_EXPAN, known as “m”, needs to be in [10^-6 ft].

For example, consider the test case below. Only pipe C2 has an EN_LEAK_AREA = 20 and EN_LEAK_EXPAN = 0.5. The resulting pressure heads in the upstream and downstream junctions are 45.796m and 42.992m. If the units of “m” where in [10^-6 m], the leakage flow rate EN_LINK_LEAKAGE should be 4.6048Lps:

However, the current simulation gives a leakage flow rate of 3.1107Lps, which corresponds to using the previous equation but with “m” in [10^-6 ft]:

Given that most literature available on the values of “m” reports it in SI units (mainly by Kobus; e.g., Greyvenstein & Van Zyl, 2007), it seems that SI projects should use [10^-6 m] as the units for EN_LEAK_EXPAN.

`[TITLE]

[JUNCTIONS]
;;ID Elev Demand Pattern
J1 0.0000
J3 0.0000
J6 0.0000

[RESERVOIRS]
;;ID Head Pattern
R1 50.0000

[PIPES]
;;ID Node1 Node2 Length Diameter Roughness MinorLoss Status
C1 R1 J1 200.0000 100.0000 100.0000 0.0000
C2 J1 J3 200.0000 100.0000 100.0000 0.0000
C3 J3 J6 200.0000 100.0000 100.0000 0.0000

[DEMANDS]
;;Junction Demand Pattern Category
J6 5.000000

[LEAKAGE]
;;Pipe Leak Area Leak Expansion
C2 20.000000 0.500000

[ENERGY]
GLOBAL EFFIC 75.0000
DEMAND CHARGE 0.0000

[TIMES]
DURATION 0:00:00
HYDRAULIC TIMESTEP 1:00:00
QUALITY TIMESTEP 0:05:00
REPORT TIMESTEP 1:00:00
REPORT START 0:00:00
PATTERN TIMESTEP 1:00:00
PATTERN START 0:00:00
RULE TIMESTEP 0:06:00
START CLOCKTIME 0:00:00
STATISTIC NONE

[OPTIONS]
UNITS LPS
PRESSURE METERS
HEADLOSS H-W
UNBALANCED CONTINUE 10
QUALITY NONE
DEMAND MULTIPLIER 1.0000
EMITTER EXPONENT 0.5000
BACKFLOW ALLOWED YES
VISCOSITY 1.000000
DIFFUSIVITY 1.000000
SPECIFIC GRAVITY 1.000000
TRIALS 40
ACCURACY 0.00100000
TOLERANCE 0.01000000
CHECKFREQ 2
MAXCHECK 10
DAMPLIMIT 0.00000000

[REPORT]
PAGESIZE 0
STATUS NO
SUMMARY NO
ENERGY NO
MESSAGES YES
NODES NONE
LINKS NONE
Elevation NO
Demand PRECISION 2
Head PRECISION 2
Pressure PRECISION 2
Quality NO
Length NO
Diameter NO
Flow PRECISION 2
Velocity PRECISION 2
Headloss PRECISION 2
Quality NO
State NO
Setting NO
Reaction NO

[TAGS]
;;Object ID Tag

[COORDINATES]
;;ID X-Coord Y-Coord
J1 -391.000000 210.500000
J3 -240.000000 212.500000
J6 -116.000000 216.500000
R1 -466.000000 146.500000

[END]
`

[LRossman]

There are some errors in the leakage hand calculation. Instead of 32.2 for g it should be 9.81 and the leak area a should be multiplied by 1e-6 to convert it from sq mm to sq m. But regardless of that, I discovered errors in the source code in converting the leakage parameters from metric units to US units for use in the FAVAD calculations.

The units of the leakage parameters supplied to the program are fixed at sq mm for leak area and sq mm per m for leak expansion, regardless of which unit system is being used by the rest of the input data, since this is how they are typically stated in the literature and will have "reasonable" values (imagine trying to express them in units of feet). To insure that the correct conversions are made the current code needs to be modified as follows:

• In both convertunits in input1.c and EN_setlinkvalue in epanet.c, the leak area and expansion properties should be multiplied by Ucf[LENGTH] to convert them to per 100 ft of pipe when metric units are used.
• The inverse of the above factors should be used when converting the leak parameters back to their original units in EN_getlinkvalue and in the inpfile.c module.
• In leakage.c, leak area should be divided by 1.e6*SQR(MperFT) and the area expansion rate by 1e6*MperFT to convert them to feet units.

I made these code changes and ran the example through it, checking by hand that the results were correct. I did this both for the original example that used metric units (lps flow units) and for the same example using cfs flow units and therefore US units. When converting the original leak parameters to US units, they had to be multiplied by 0.3048 to make them relative to 100 ft of pipe length rather than 100 m. The results had the leakage in Pipe C2 be the same for both examples, at 1.55 lps or 0.0548 cfs.

It appears that this conversion problem is serious enough to warrant an immediate patch release (instead of a 1.55 lps leakage result the current DLL is showing a leakage of 3.11 lps for the example network).

[dapaez270]

Thanks for taking a look at this so quickly, Lew. I should have mentioned earlier that the software I used for the hand verification (SMath Studio) handles unit conversions internally when they are shown in blue (e.g., it automatically applies the 1e-6 factor when converting from mm² to m²)

I agree with the solution approach. It feels more intuitive to specify the area in mm² and the expansion rate in mm²/m, which avoids having to pass most values with an explicit 1e-6 factor. Expressing both per 100 ft or per 100 m of pipe length (depending on the project units) also makes sense (I had originally assumed they were always per 100 ft)

I’m looking forward to testing the new patch. The updated 1.55 L/s value seems reasonable, although I’d need the updated pressure head values at the junctions to verify it by hand.

[dapaez270]

I also just noticed that c_orif in leakage.c might have a small typo. The computation 0.6 * sqrt(2 * 32.2) = 4.8149766 . However, the value used in lines 158 and 293 is 4.8149866 (i.e., there seems to be a typo in the fifth decimal place).

This probably does not result in a significant error in the final computed values, though.

[LRossman]

I edited my earlier comment to reflect the actual changes I made to the leakage-related code. I also made the correction you noted above. These changes reside in a new branch named dev-leakage_units_fix. I will wait a while before submitting a PR for so that it can be tested some more. Please let me know if you see any problems with it.

[LRossman]

PR #901 fixes this issue and has been merged into dev. This issue can now be closed.