Add FAVA leakage model to EPANET

[LRossman]

Pipe leakage is typically modeled with an emitter at a node, using a power relationship between leakage flow rate and pressure. An issue arises in assigning meaningful values to the emitter flow coefficient and the global emitter exponent to be used. Conceptually the flow coefficient should somehow reflect the number of pipes connected to the node, the assumed number of cracks and their areas in each connecting pipe, and the pressure experienced along the pipe (as opposed to just the node). Regarding the flow exponent, empirical studies have shown that it can also vary with pressure.

The FAVAD leakage model (Cassa and van Zyl, 2014) uses the following equation to estimate the leakage rate from a single crack along a pipe:

   flow = Cd * sqrt(2g) * (Ao*h^0.5 + m*h^1.5)

where Cd is an orifice discharge coefficient (typically 0.6), Ao is the area of the crack under zero pressure conditions, h is the pressure head seen by the crack, and m is rate at which the crack area expands as pressure increases (dA/dh with experimental values between 0 and 0.001 mm depending on pipe material and crack size).

This model can easily be incorporated into EPANET. An additional set of leakage parameters would be defined for each pipe: LC1 = the area of leaks per 100 units of pipe length and LC2 = the leak expansion rate per 100 units of pipe length (i.e., LC1 = Ao * N * L / 100, LC2 = m * N * L / 100 where N is the average number of cracks per 100 length units of pipe and L is pipe length). The FAVAD equation would be evaluated using these parameters along with the average pressure head along the pipe (average of the upstream and downstream node pressures). The resulting leakage flow would then be divided evenly between the pipe's end nodes and incorporated into EPANET's hydraulic solver in similar fashion to how emitter flow is treated.

Admittedly this results in substituting the three pipe-related leakage parameters (Ao, m, and N) for the two node-related emitter parameters (flow coefficient and exponent). However it does provide a stronger physical basis for modeling pipe leaks.

[LRossman]

I've incorporated the FAVAD leakage model into a fork of the dev branch here. Since it touches a significant number of files I'd like to get some buy in on it before merging it into this repo. The updates it includes can be summarized as follows:

• A new [LEAKAGE] section has been added to the input file format where each line contains the ID name of a pipe, its leak area in sq. mm per 100 length units, and its leak expansion rate in sq. mm per unit of pressure head.
• EN_LEAK_AREA and EN_LEAK_EXPAN can be used with the functions EN_getlinkvalue and EN_setlinkvalue to retrieve and assign values for a pipe's leak area and expansion properties.
• EN_LINK_LEAKAGE can be used with EN_getlinkvalue to retrieve a pipe's leakage rate at a given point in time.
• EN_LEAKAGEFLOW can be used with EN_getnodevalue to retrieve the leakage demand generated at a node from all its connecting pipes at a given point in time.
• EN_LEAKAGELOSS can be used with EN_getstatistic to retrieve the the total leakage loss in the system at a given point in time as a percentage of total flow entering the system.
• A new Flow Balance Report has been added to end of a simulation run's Status Report that lists the various components of the system's total inflow and outflow over the simulation period. It also displays the ratio of outflow to inflow as a check on flow continuity.
• The following constants can be used with EN_getnodevalue to retrieve the components of a node's total demand at a given point in time:
  • EN_FULLDEMAND - the consumer demand requested
  • EN_DEMANDFLOW - the consumer demand delivered
  • EN_DEMANDDEFICIT - the difference between the consumer demand requested and delivered
  • EN_EMITTERFLOW - the node's emitter flow
  • EN_LEAKAGEFLOW - the node's leakage flow
  • EN_DEMAND - the sum of the node's consumer demand, emitter flow, and leakage flow

I've structured the leakage code so that it doesn't get executed if no pipes are assigned leakage parameters. Also included is a simple unit test file. As an example of the kind of results it can generate I ran the 2466 pipe EXNET-3 model under a PDA scenario where each pipe was assigned a fixed leakage area of 1 sq mm per 100 m of length and an expansion rate of 0.1 sq mm per meter of pressure head. Below is a color-coded map of the resulting leakage rate in each pipe, where the total system leakage loss was 20%. It's interesting to see how non-uniform the leakage results are and would be difficult to predict given the size of this system.

[kobusvanzyl]

Thank you Lew!

It's great seeing you incorporate the more realistic FAVAD leakage model into EPANET.

Your implementation allows the user to specify leak properties per 100-length units, which will be useful when distributing leakage proportional to pipe length (possibly considering material, diameter, etc.). However, this is not a realistic representation of leakage. We recently analysed leak dimensions and found that leak lengths/diameters typically adhere to log-normal distributions. Thus, most leaks will be small, with a few large leaks contributing most of the leakage flow.

We may not know the distribution of leaks in a given system, but we can use stochastic analysis to investigate a range of scenarios based on known leak property distributions. Many researchers investigate leak detection by allocating large leaks to particular nodes, with or without 'background' leaks distributed through the system.

I think there are some strong arguments to associate leaks with nodes rather than pipes in the EPANET model:

1. Leaks are demands just like emitters and fixed demands (which are allocated to nodes).
2. It will give the user direct control over how leaks are allocated throughout the system, thus simplifying the process of allocating single large leaks or a given distribution of leak sizes.
3. Uniform leakage can be achieved with relative ease by estimating the leakage properties of each pipe and then splitting these between the start and end nodes.

[kobusvanzyl]

We are starting to get a pretty good understanding of the head-area slope m of individual leaks. Here are some suggestions:

Smaller diameter steel and iron
m is very small (typically between -0.02 mm2/m and +0.02 mm2/m) and can be assumed zero (i.e. orifice flow).

Round holes in any material
m is very small (typically between -0.02 mm2/m and +0.02 mm2/m) and can be assumed zero (i.e. orifice flow).

Circumferential cracks
m range is wider han for round holes, but is still pretty small (typically between -0.5 mm2/m and +0.5 mm2/m).
m is often negative (crack area reduces with increasing pressure).
It is reasonable to assume m = 0.

Longitudinal cracks
m is positive and can be very large. The largest m we've measured is 170 mm2/m for a 200 mm longitudinal slit in a 50 mm diameter class PN6 PVC pipe.
Li et al. (2022) recently demonstrated that the physics-based Tada-Paris equation provides good estimates for m.

Reference
Yuanzhe Li , Jinliang Gao , Chang Shen , Yu Guan and Wei Wang (2022), Water Research, 221, 118692.

[LRossman]

Thanks Kobus for the feedback. I just wanted to clarify that the new leakage option converts the FAVAD parameters supplied to pipes into an equivalent set of coefficients for a pair of virtual emitters (one with a 0.5 exponent and the other with a 1.5 exponent) placed at their end nodes. So "under the hood" leakage is still being computed using a nodal emitter approach. Please also note that one doesn't have to assume that all pipes have the same leak areas and expansion rates. Different values can be assigned on a pipe-by-pipe basis. I chose to normalize the FAVAD parameters to 100 units of pipe length to make their values easier to work with (e.g., a typical length of a delivered pipe section is 20 to 40 feet, so 100 feet of pipe would have 4 joints which could be possible sources of leaks in addition to any cracks).

It's interesting to compare results from the pipe-based leakage approach to that of the node-based approach where, for lack of any better information on the location and size of leak areas, all pipes (or nodes) are assigned the same leakage parameters. The figures below show computed leakage rates at the nodes of the same EXNET-3 model used earlier with a pressure dependent demand model between 20 to 40 m. The top figure is for the FAVAD model with each pipe assigned a fixed leakage area of 1 sq mm per 100 m of length and an expansion rate of 0.1 sq mm per meter of pressure head. This results in a system leakage loss of 15%. The bottom figure is for the nodal emitter approach where each node is assigned an emitter coefficient of 0.05 which also produces a 15% system loss. The nodal leakage approach is not identifying the leakage hot spots the way that the pipe leakage method does (because it doesn't take into account the number of pipes and their lengths that connect to a node).

If as you say that leakage is usually dominated by a few large cracks then these can still be accounted for in EPANET's FAVAD model by simply assigning those particular pipes an appropriate set of leak area coefficients.

[kobusvanzyl]

Thanks, Lew. Using two virtual emitters, one with a 0.5 exponent and the other with a 1.5 exponent, is exactly right. The 0.5 exponent term describes the (orifice) leakage flow through the leak opening as if fixed. The 1.5 exponent term calculates the leakage flow through the expanded part of the leak opening.

I see your point that assigning the leak properties to pipes makes better sense from a user perspective and allows for a more realistic distribution across the network. Makes complete sense.