From b0796f346a4a6de9718a6ed4097dc8699c2fdc81 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Wed, 26 Jun 2024 11:15:01 -0400 Subject: [PATCH 1/5] Adding leakage model --- ReleaseNotes2_3.md | 17 +- include/epanet2.bas | 11 +- include/epanet2.cs | 25 +- include/epanet2.def | 17 +- include/epanet2.pas | 306 ++++++++++++----------- include/epanet2.vb | 14 +- include/epanet2_enums.h | 13 +- src/enumstxt.h | 4 +- src/epanet.c | 66 ++++- src/flowbalance.c | 191 +++++++++++++++ src/funcs.h | 19 +- src/hydcoeffs.c | 5 +- src/hydraul.c | 19 +- src/hydsolver.c | 93 ++++--- src/inpfile.c | 17 +- src/input1.c | 6 +- src/input2.c | 1 + src/input3.c | 41 ++++ src/leakage.c | 531 ++++++++++++++++++++++++++++++++++++++++ src/project.c | 11 +- src/report.c | 38 +++ src/text.h | 6 +- src/types.h | 43 +++- tests/CMakeLists.txt | 1 + tests/test_leakage.cpp | 93 +++++++ 25 files changed, 1366 insertions(+), 222 deletions(-) create mode 100644 src/flowbalance.c create mode 100644 src/leakage.c create mode 100644 tests/test_leakage.cpp diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md index ea5e96de..e561ee81 100644 --- a/ReleaseNotes2_3.md +++ b/ReleaseNotes2_3.md @@ -4,6 +4,7 @@ This document describes the changes and updates that have been made in version 2.3 of EPANET. - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files could be opened by the toolkit. + - A `EN_setcurvetype` function was added to allow API clients to set a curve's type (e.g., `EN_PUMP_CURVE,` `EN_VOLUME_CURVE,` etc.). - A `EN_setvertex` function was added to allow API clients to change the coordinates of a single link vertex. - The indices of a General Purpose Valve (GPV) and a Positional Control Valve (PCV) were added to the list of editable Link Properties using the symbolic constant names `EN_GPV_CURVE` and `EN_PCV_CURVE`, respectively. @@ -52,6 +53,20 @@ This document describes the changes and updates that have been made in version 2 - Setting the demand multiplier within the `[DEMANDS]` section of INP has been depreciated, please use `DEMAND MULTIPLIER` inside `[OPTIONS]` instead. - `EN_PRESS_UNITS` can now be used with `EN_getoption` and `EN_setoption` to get or set the pressure unit used in EPANET. - Continuous barrier functions were added to constrain emitter flows to allowable values. -- The `EN_openx` function has been added to enable the opening of input files with formatting errors through the API. This allows users to continue using toolkit functions even when such errors are present. +- The `EN_openX` function has been added to enable the opening of input files with formatting errors through the API. This allows users to continue using toolkit functions even when such errors are present. - The `EN_getnodesvalues` and `EN_getlinksvalues` were added to retrieve a property value for all nodes or links in the network. - Fixed a bug in EN_setnodevalue with EN_EMITTER option that could cause NaN results. +- Support has been added for FAVAD (Fixed And Variable Area Discharge) modeling of pipe leaks: + - 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 \ No newline at end of file diff --git a/include/epanet2.bas b/include/epanet2.bas index 6d707591..6ddb0ada 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -5,7 +5,7 @@ Attribute VB_Name = "Module1" 'Declarations of functions in the EPANET PROGRAMMERs TOOLKIT '(EPANET2.DLL) -'Last updated on 09/28/2023 +'Last updated on 06/23/2024 ' These are codes used by the DLL functions Public Const EN_ELEVATION = 0 ' Node parameters @@ -38,6 +38,9 @@ Public Const EN_CANOVERFLOW = 26 Public Const EN_DEMANDDEFICIT = 27 Public Const EN_NODE_INCONTROL = 28 Public Const EN_EMITTERFLOW = 29 +Public Const EN_LEAKAGEFLOW = 30 +Public Const EN_DEMANDFLOW = 31 +Public Const EN_FULLDEMAND = 32 Public Const EN_DIAMETER = 0 ' Link parameters Public Const EN_LENGTH = 1 @@ -65,6 +68,9 @@ Public Const EN_PUMP_EPAT = 22 Public Const EN_LINK_INCONTROL = 23 Public Const EN_GPV_CURVE = 24 Public Const EN_PCV_CURVE = 25 +Public Const EN_LEAK_AREA = 26 +Public Const EN_LEAK_EXPAN = 27 +Public Const EN_LINK_LEAKAGE = 28 Public Const EN_DURATION = 0 ' Time parameters Public Const EN_HYDSTEP = 1 @@ -90,6 +96,7 @@ Public Const EN_MAXFLOWCHANGE = 3 Public Const EN_MASSBALANCE = 4 Public Const EN_DEFICIENTNODES = 5 Public Const EN_DEMANDREDUCTION = 6 +Public Const EN_LEAKAGELOSS = 7 Public Const EN_NODE = 0 ' Component types Public Const EN_LINK = 1 @@ -350,6 +357,7 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENsettankdata Lib "epanet2.dll" (ByVal index As Long, ByVal elev As Single, ByVal initlvl As Single, ByVal minlvl As Single, ByVal maxlvl As Single, ByVal diam As Single, ByVal minvol As Single, ByVal volcurve As String) As Long Declare Function ENgetcoord Lib "epanet2.dll" (ByVal index As Long, x As Double, y As Double) As Long Declare Function ENsetcoord Lib "epanet2.dll" (ByVal index As Long, ByVal x As Double, ByVal y As Double) As Long + Declare Function ENgetnodevalues Lib "epanet2.dll" (ByVal property as Long, values as Any) As Long 'Nodal Demand Functions Declare Function ENgetdemandmodel Lib "epanet2.dll" (type_ As Long, pmin As Single, preq As Single, pexp As Single) As Long @@ -382,6 +390,7 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENgetvertex Lib "epanet2.dll" (ByVal index As Long, ByVal vertex As Long, x As Double, y As Double) As Long Declare Function ENsetvertex Lib "epanet2.dll" (ByVal index As Long, ByVal vertex As Long, ByVal x As Double, ByVal y As Double) As Long Declare Function ENsetvertices Lib "epanet2.dll" (ByVal index As Long, xCoords As Any, yCoords As Any, ByVal count As Long) As Long + Declare Function ENgetlinkvalues Lib "epanet2.dll" (ByVal property as Long, values as Any) As Long 'Pump Functions Declare Function ENgetheadcurveindex Lib "epanet2.dll" (ByVal linkIndex As Long, curveIndex As Long) As Long diff --git a/include/epanet2.cs b/include/epanet2.cs index 57668bf7..0327c4d8 100644 --- a/include/epanet2.cs +++ b/include/epanet2.cs @@ -3,7 +3,7 @@ using System.Runtime.InteropServices; //epanet2.cs[By Oscar Vegas] -//Last updated on 09/28/2023 +//Last updated on 06/23/2024 //Declarations of functions in the EPANET PROGRAMMERs TOOLKIT //(EPANET2.DLL) for use with C# @@ -50,6 +50,9 @@ public static class Epanet public const int EN_DEMANDDEFICIT = 27; public const int EN_NODE_INCONTROL = 28; public const int EN_EMITTERFLOW = 29; + public const int EN_LEAKAGEFLOW = 30; + public const int EN_DEMANDFLOW = 31; + public const int EN_FULLDEMAND = 32; public const int EN_DIAMETER = 0; //Link parameters public const int EN_LENGTH = 1; @@ -78,6 +81,9 @@ public static class Epanet public const int EN_LINK_INCONTROL = 23; public const int EN_GPV_CURVE = 24; public const int EN_PCV_CURVE = 25; + public const int EN_LEAK_AREA = 26; + public const int EN_LEAK_EXPAN = 27; + public const int EN_LINK_LEAKAGE = 28; public const int EN_DURATION = 0; //Time parameters public const int EN_HYDSTEP = 1; @@ -102,6 +108,7 @@ public static class Epanet public const int EN_MASSBALANCE = 4; public const int EN_DEFICIENTNODES = 5; public const int EN_DEMANDREDUCTION = 6; + public const int EN_LEAKAGELOSS = 7; public const int EN_NODE = 0; //Component types public const int EN_LINK = 1; @@ -390,6 +397,9 @@ public static class Epanet [DllImport(EPANETDLL, EntryPoint = "ENgetresultindex")] public static extern int ENgetresultindex(int type, int index, ref int value); + [DllImport(EPANETDLL, EntryPoint = "ENtimetonextevent")] + public static extern int ENtimetonextevent(ref int eventType, ref long duration, ref int elementIndex); + //Analysis Options Functions [DllImport(EPANETDLL, EntryPoint = "ENgetoption")] @@ -440,10 +450,10 @@ public static class Epanet public static extern int ENgetnodetype(int index, ref int nodeType); [DllImport(EPANETDLL, EntryPoint = "ENgetnodevalue")] - public static extern int ENgetnodevalue(int index, int paramcode, ref float value); + public static extern int ENgetnodevalue(int index, int param, ref float value); [DllImport(EPANETDLL, EntryPoint = "ENsetnodevalue")] - public static extern int ENsetnodevalue(int index, int code, float value); + public static extern int ENsetnodevalue(int index, int param, float value); [DllImport(EPANETDLL, EntryPoint = "ENsetjuncdata")] public static extern int ENsetjuncdata(int index, float elev, float dmnd, string dmndpat); @@ -457,6 +467,8 @@ public static class Epanet [DllImport(EPANETDLL, EntryPoint = "ENsetcoord")] public static extern int ENsetcoord(int index, double x, double y); + [DllImport(EPANETDLL, EntryPoint = "ENgetnodevalues")] + public static extern int ENgetnodevalues(int param, ref float values); //Nodal Demand Functions [DllImport(EPANETDLL, EntryPoint = "ENgetdemandmodel")] @@ -525,10 +537,10 @@ public static class Epanet public static extern int ENsetlinknodes(int index, int node1, int node2); [DllImport(EPANETDLL, EntryPoint = "ENgetlinkvalue")] - public static extern int ENgetlinkvalue(int index, int code, ref float value); + public static extern int ENgetlinkvalue(int index, int param, ref float value); [DllImport(EPANETDLL, EntryPoint = "ENsetlinkvalue")] - public static extern int ENsetlinkvalue(int index, int code, float value); + public static extern int ENsetlinkvalue(int index, int param, float value); [DllImport(EPANETDLL, EntryPoint = "ENsetpipedata")] public static extern int ENsetpipedata(int index, float length, float diam, float rough, float mloss); @@ -542,6 +554,9 @@ public static class Epanet [DllImport(EPANETDLL, EntryPoint = "ENsetvertices")] public static extern int ENsetvertices(int index, ref double[] x, ref double[] y, int count); + [DllImport(EPANETDLL, EntryPoint = "ENgetlinkvalues")] + public static extern int ENgetlinkvalues(int param, ref float values); + //Pump Functions [DllImport(EPANETDLL, EntryPoint = "ENgetheadcurveindex")] diff --git a/include/epanet2.def b/include/epanet2.def index f5372af4..d3472c10 100644 --- a/include/epanet2.def +++ b/include/epanet2.def @@ -25,6 +25,7 @@ EXPORTS ENgetbasedemand = _ENgetbasedemand@12 ENgetcomment = _ENgetcomment@12 ENgetcontrol = _ENgetcontrol@24 + ENgetcontrolenabled = _ENgetcontrolenabled@8 ENgetcoord = _ENgetcoord@12 ENgetcount = _ENgetcount@8 ENgetcurve = _ENgetcurve@20 @@ -44,13 +45,14 @@ EXPORTS ENgetlinkid = _ENgetlinkid@8 ENgetlinkindex = _ENgetlinkindex@8 ENgetlinknodes = _ENgetlinknodes@12 - ENsetlinknodes = _ENsetlinknodes@12 ENgetlinktype = _ENgetlinktype@8 ENgetlinkvalue = _ENgetlinkvalue@12 + ENgetlinkvalues = _ENgetlinkvalues@8 ENgetnodeid = _ENgetnodeid@8 ENgetnodeindex = _ENgetnodeindex@8 ENgetnodetype = _ENgetnodetype@8 - ENgetnodevalue = _ENgetnodevalue@12 + ENgetnodevalue = _ENgetnodevalue@12 + ENgetnodevalues = _ENgetnodevalues@8 ENgetnumdemands = _ENgetnumdemands@8 ENgetoption = _ENgetoption@8 ENgetpatternid = _ENgetpatternid@8 @@ -63,6 +65,7 @@ EXPORTS ENgetqualtype = _ENgetqualtype@8 ENgetresultindex = _ENgetresultindex@12 ENgetrule = _ENgetrule@20 + ENgetruleenabled = _ENgetruleenabled@8 ENgetruleID = _ENgetruleID@8 ENgetstatistic = _ENgetstatistic@8 ENgetthenaction = _ENgetthenaction@20 @@ -79,6 +82,7 @@ EXPORTS ENopen = _ENopen@12 ENopenH = _ENopenH@0 ENopenQ = _ENopenQ@0 + ENopenX = _ENopenX@12 ENreport = _ENreport@0 ENresetreport = _ENresetreport@0 ENrunH = _ENrunH@4 @@ -89,6 +93,7 @@ EXPORTS ENsetbasedemand = _ENsetbasedemand@12 ENsetcomment = _ENsetcomment@12 ENsetcontrol = _ENsetcontrol@24 + ENsetcontrolenabled = _ENsetcontrolenabled@8 ENsetcoord = _ENsetcoord@20 ENsetcurve = _ENsetcurve@16 ENsetcurveid = _ENsetcurveid@8 @@ -118,6 +123,7 @@ EXPORTS ENsetpremisevalue = _ENsetpremisevalue@12 ENsetqualtype = _ENsetqualtype@16 ENsetreport = _ENsetreport@4 + ENsetruleenabled = _ENsetruleenabled@8 ENsetrulepriority = _ENsetrulepriority@8 ENsetstatusreport = _ENsetstatusreport@4 ENsettankdata = _ENsettankdata@32 @@ -129,11 +135,6 @@ EXPORTS ENsolveH = _ENsolveH@0 ENsolveQ = _ENsolveQ@0 ENstepQ = _ENstepQ@4 + ENtimetonextevent = _ENtimetonextevent@12 ENusehydfile = _ENusehydfile@4 ENwriteline = _ENwriteline@4 - ENtimetonextevent = _ENtimetonextevent@12 - ENopenX = _ENopenX@12 - ENgetcontrolenabled = _ENgetcontrolenabled@8 - ENsetcontrolenabled = _ENsetcontrolenabled@8 - ENgetruleenabled = _ENgetruleenabled@8 - ENsetruleenabled = _ENsetruleenabled@8 diff --git a/include/epanet2.pas b/include/epanet2.pas index 4ce4a03c..6f2e8dae 100644 --- a/include/epanet2.pas +++ b/include/epanet2.pas @@ -3,7 +3,7 @@ { Declarations of imported procedures from the EPANET PROGRAMMERs TOOLKIT } { (EPANET2.DLL) } -{Last updated on 09/28/2023} +{Last updated on 06/06/2024} interface @@ -46,6 +46,9 @@ interface EN_DEMANDDEFICIT = 27; EN_NODE_INCONTROL = 28; EN_EMITTERFLOW = 29; + EN_LEAKAGEFLOW = 30; + EN_DEMANDFLOW = 31; + EN_FULLDEMAND = 32; EN_DIAMETER = 0; { Link parameters } EN_LENGTH = 1; @@ -71,8 +74,11 @@ interface EN_PUMP_ECOST = 21; EN_PUMP_EPAT = 22; EN_LINK_INCONTROL = 23; - EN_GPV_CURVE = 24; - EN_PCV_CURVE = 25; + EN_GPV_CURVE = 24; + EN_PCV_CURVE = 25; + EN_LEAK_AREA = 26; + EN_LEAK_EXPAN = 27; + EN_LINK_LEAKAGE = 28; EN_DURATION = 0; { Time parameters } EN_HYDSTEP = 1; @@ -98,6 +104,7 @@ interface EN_MASSBALANCE = 4; EN_DEFICIENTNODES = 5; EN_DEMANDREDUCTION = 6; + EN_LEAKAGELOSS = 7; EN_NODE = 0; { Component Types } EN_LINK = 1; @@ -275,180 +282,199 @@ interface EN_TRUE = 1; { boolean true } +{$MACRO ON} + {$ifdef MSWINDOWS} - EpanetLib = 'epanet2.dll'; + EpanetLib = 'epanet2.dll'; + {$DEFINE cdecl := stdcall} +{$endif} + +{$ifdef UNIX} + {$ifdef DARWIN} + EpanetLib = 'libepanet2.dylib'; + {$linklib libepanet2} + {$else} + EpanetLib = 'libepanet2.so'; + {$endif} +{$endif} + +{$ifdef UNIX} + {$DEFINE TimeType:=Int64} {$else} - EpanetLib = 'libepanet2.so'; + {$DEFINE TimeType:=Integer} {$endif} {Project Functions} - function ENepanet(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar; F4: Pointer): Integer; stdcall; external EpanetLib; - function ENinit(F2: PAnsiChar; F3: PAnsiChar; UnitsType: Integer; HeadlossType: Integer): Integer; stdcall; external EpanetLib; - function ENopen(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENopenX(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetcount(Code: Integer; var Count: Integer): Integer; stdcall; external EpanetLib; - function ENgettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsaveinpfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENclose: Integer; stdcall; external EpanetLib; + function ENepanet(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar; F4: Pointer): Integer; cdecl; external EpanetLib; + function ENinit(F2: PAnsiChar; F3: PAnsiChar; UnitsType: Integer; HeadlossType: Integer): Integer; cdecl; external EpanetLib; + function ENopen(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENopenX(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetcount(Code: Integer; var Count: Integer): Integer; cdecl; external EpanetLib; + function ENgettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsaveinpfile(F: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENclose: Integer; cdecl; external EpanetLib; {Hydraulic Functions} - function ENsolveH: Integer; stdcall; external EpanetLib; - function ENsaveH: Integer; stdcall; external EpanetLib; - function ENopenH: Integer; stdcall; external EpanetLib; - function ENinitH(SaveFlag: Integer): Integer; stdcall; external EpanetLib; - function ENrunH(var T: Integer): Integer; stdcall; external EpanetLib; - function ENnextH(var Tstep: Integer): Integer; stdcall; external EpanetLib; - function ENcloseH: Integer; stdcall; external EpanetLib; - function ENsavehydfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENusehydfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsolveH: Integer; cdecl; external EpanetLib; + function ENsaveH: Integer; cdecl; external EpanetLib; + function ENopenH: Integer; cdecl; external EpanetLib; + function ENinitH(SaveFlag: Integer): Integer; cdecl; external EpanetLib; + function ENrunH(var T: TimeType): Integer; cdecl; external EpanetLib; + function ENnextH(var Tstep: TimeType): Integer; cdecl; external EpanetLib; + function ENcloseH: Integer; cdecl; external EpanetLib; + function ENsavehydfile(F: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENusehydfile(F: PAnsiChar): Integer; cdecl; external EpanetLib; {Quality Functions} - function ENsolveQ: Integer; stdcall; external EpanetLib; - function ENopenQ: Integer; stdcall; external EpanetLib; - function ENinitQ(SaveFlag: Integer): Integer; stdcall; external EpanetLib; - function ENrunQ(var T: Integer): Integer; stdcall; external EpanetLib; - function ENnextQ(var Tstep: Integer): Integer; stdcall; external EpanetLib; - function ENstepQ(var Tleft: Integer): Integer; stdcall; external EpanetLib; - function ENcloseQ: Integer; stdcall; external EpanetLib; + function ENsolveQ: Integer; cdecl; external EpanetLib; + function ENopenQ: Integer; cdecl; external EpanetLib; + function ENinitQ(SaveFlag: Integer): Integer; cdecl; external EpanetLib; + function ENrunQ(var T: TimeType): Integer; cdecl; external EpanetLib; + function ENnextQ(var Tstep: TimeType): Integer; cdecl; external EpanetLib; + function ENstepQ(var Tleft: TimeType): Integer; cdecl; external EpanetLib; + function ENcloseQ: Integer; cdecl; external EpanetLib; {Reporting Functions} - function ENwriteline(S: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENreport: Integer; stdcall; external EpanetLib; - function ENcopyreport(F: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENclearreport: Integer; stdcall; external EpanetLib; - function ENresetreport: Integer; stdcall; external EpanetLib; - function ENsetreport(S: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetstatusreport(Code: Integer): Integer; stdcall; external EpanetLib; - function ENgetversion(var Version: Integer): Integer; stdcall; external EpanetLib; - function ENgeterror(Errcode: Integer; Errmsg: PAnsiChar; MaxLen: Integer): Integer; stdcall; external EpanetLib; - function ENgetstatistic(StatType: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENgetresultindex(Code: Integer; Index: Integer; var Value: Integer): Integer; stdcall; external EpanetLib; + function ENwriteline(S: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENreport: Integer; cdecl; external EpanetLib; + function ENcopyreport(F: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENclearreport: Integer; cdecl; external EpanetLib; + function ENresetreport: Integer; cdecl; external EpanetLib; + function ENsetreport(S: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetstatusreport(Code: Integer): Integer; cdecl; external EpanetLib; + function ENgetversion(var Version: Integer): Integer; cdecl; external EpanetLib; + function ENgeterror(Errcode: Integer; Errmsg: PAnsiChar; MaxLen: Integer): Integer; cdecl; external EpanetLib; + function ENgetstatistic(StatType: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENgetresultindex(Code: Integer; Index: Integer; var Value: Integer): Integer; cdecl; external EpanetLib; + function ENtimetonextevent(var EventType: Integer; var Duration: TimeType; var ElementIndex: Integer): Integer; cdecl; external EpanetLib; {Analysis Options Functions} - function ENgetoption(Code: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetoption(Code: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENgetflowunits(var Code: Integer): Integer; stdcall; external EpanetLib; - function ENsetflowunits(Code: Integer): Integer; stdcall; external EpanetLib; - function ENgettimeparam(Code: Integer; var Value: Integer): Integer; stdcall; external EpanetLib; - function ENsettimeparam(Code: Integer; Value: Integer): Integer; stdcall; external EpanetLib; - function ENgetqualinfo(var QualType: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; var TraceNode: Integer): Integer; stdcall; external EpanetLib; - function ENgetqualtype(var QualCode: Integer; var TraceNode: Integer): Integer; stdcall; external EpanetLib; - function ENsetqualtype(QualCode: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; TraceNodeID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetoption(Code: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetoption(Code: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENgetflowunits(var Code: Integer): Integer; cdecl; external EpanetLib; + function ENsetflowunits(Code: Integer): Integer; cdecl; external EpanetLib; + function ENgettimeparam(Code: Integer; var Value: TimeType): Integer; cdecl; external EpanetLib; + function ENsettimeparam(Code: Integer; Value: TimeType): Integer; cdecl; external EpanetLib; + function ENgetqualinfo(var QualType: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; var TraceNode: Integer): Integer; cdecl; external EpanetLib; + function ENgetqualtype(var QualCode: Integer; var TraceNode: Integer): Integer; cdecl; external EpanetLib; + function ENsetqualtype(QualCode: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; TraceNodeID: PAnsiChar): Integer; cdecl; external EpanetLib; {Node Functions} - function ENaddnode(ID: PAnsiChar; NodeType: Integer; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENdeletenode(Index: Integer; ActionCode: Integer): Integer; stdcall; external EpanetLib; - function ENgetnodeindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetnodeid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetnodeid(Index: Integer; NewID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetnodetype(Index: Integer; var Code: Integer): Integer; stdcall; external EpanetLib; - function ENgetnodevalue(Index: Integer; Code: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetnodevalue(Index: Integer; Code: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENsetjuncdata(Index: Integer; Elev: Single; Dmnd: Single; DmndPat: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsettankdata(Index: Integer; Elev, InitLvl, MinLvl, MaxLvl, Diam, MinVol: Single; VolCurve: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetcoord(Index: Integer; var X: Double; var Y: Double): Integer; stdcall; external EpanetLib; - function ENsetcoord(Index: Integer; X: Double; Y: Double): Integer; stdcall; external EpanetLib; + function ENaddnode(ID: PAnsiChar; NodeType: Integer; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENdeletenode(Index: Integer; ActionCode: Integer): Integer; cdecl; external EpanetLib; + function ENgetnodeindex(ID: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetnodeid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetnodeid(Index: Integer; NewID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetnodetype(Index: Integer; var Code: Integer): Integer; cdecl; external EpanetLib; + function ENgetnodevalue(Index: Integer; Code: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetnodevalue(Index: Integer; Code: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENsetjuncdata(Index: Integer; Elev: Single; Dmnd: Single; DmndPat: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsettankdata(Index: Integer; Elev, InitLvl, MinLvl, MaxLvl, Diam, MinVol: Single; VolCurve: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetcoord(Index: Integer; var X: Double; var Y: Double): Integer; cdecl; external EpanetLib; + function ENsetcoord(Index: Integer; X: Double; Y: Double): Integer; cdecl; external EpanetLib; + function ENgetnodevalues(Code: Integer; var X: array of Single): Integer; cdecl; external EpanetLib; {Demand Functions} - function ENgetdemandmodel(var Model: Integer; var Pmin: Single; var Preq: Single; var Pexp: Single): Integer; stdcall; external EpanetLib; - function ENsetdemandmodel(Model: Integer; Pmin: Single; Preq: Single; Pexp: Single): Integer; stdcall; external EpanetLib; - function ENgetnumdemands(NodeIndex: Integer; var NumDemands: Integer): Integer; stdcall; external EpanetLib; - function ENadddemand(NodeIndex: Integer; BaseDemand: Single; PatName: PAnsiChar; DemandName: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENdeletedemand(NodeIndex: Integer; DemandIndex: Integer): Integer; stdcall; external EpanetLib; - function ENgetdemandindex(NodeIndex: Integer; DemandName: PAnsiChar; var DemandIndex: Integer): Integer; stdcall; external EpanetLib; - function ENgetbasedemand(NodeIndex: Integer; DemandIndex: Integer; var BaseDemand: Single): Integer; stdcall; external EpanetLib; - function ENsetbasedemand(NodeIndex: Integer; DemandIndex: Integer; BaseDemand: Single): Integer; stdcall; external EpanetLib; - function ENgetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; var PatIndex: Integer): Integer; stdcall; external EpanetLib; - function ENsetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; PatIndex: Integer): Integer; stdcall; external EpanetLib; - function ENgetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetdemandmodel(var Model: Integer; var Pmin: Single; var Preq: Single; var Pexp: Single): Integer; cdecl; external EpanetLib; + function ENsetdemandmodel(Model: Integer; Pmin: Single; Preq: Single; Pexp: Single): Integer; cdecl; external EpanetLib; + function ENgetnumdemands(NodeIndex: Integer; var NumDemands: Integer): Integer; cdecl; external EpanetLib; + function ENadddemand(NodeIndex: Integer; BaseDemand: Single; PatName: PAnsiChar; DemandName: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENdeletedemand(NodeIndex: Integer; DemandIndex: Integer): Integer; cdecl; external EpanetLib; + function ENgetdemandindex(NodeIndex: Integer; DemandName: PAnsiChar; var DemandIndex: Integer): Integer; cdecl; external EpanetLib; + function ENgetbasedemand(NodeIndex: Integer; DemandIndex: Integer; var BaseDemand: Single): Integer; cdecl; external EpanetLib; + function ENsetbasedemand(NodeIndex: Integer; DemandIndex: Integer; BaseDemand: Single): Integer; cdecl; external EpanetLib; + function ENgetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; var PatIndex: Integer): Integer; cdecl; external EpanetLib; + function ENsetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; PatIndex: Integer): Integer; cdecl; external EpanetLib; + function ENgetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; cdecl; external EpanetLib; {Link Functions} - function ENaddlink(ID: PAnsiChar; LinkType: Integer; FromNode: PAnsiChar; ToNode: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENdeletelink(Index: Integer; ActionCode: Integer): Integer; stdcall; external EpanetLib; - function ENgetlinkindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetlinkid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetlinkid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetlinktype(Index: Integer; var Code: Integer): Integer; stdcall; external EpanetLib; - function ENsetlinktype(var Index: Integer; LinkType: Integer; ActionCode: Integer): Integer; stdcall; external EpanetLib; - function ENgetlinknodes(Index: Integer; var Node1: Integer; var Node2: Integer): Integer; stdcall; external EpanetLib; - function ENsetlinknodes(Index: Integer; Node1: Integer; Node2: Integer): Integer; stdcall; external EpanetLib; - function ENgetlinkvalue(Index: Integer; Code: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetlinkvalue(Index: Integer; Code: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENsetpipedata(Index: Integer; Length: Single; Diam: Single; Rough: Single; Mloss:Single): Integer; stdcall; external EpanetLib; - - function ENgetvertexcount(Index: Integer; var Count: Integer): Integer; stdcall; external EpanetLib; - function ENgetvertex(Index: Integer; Vertex: Integer; var X: Double; var Y: Double): Integer; stdcall; external EpanetLib; - function ENsetvertex(Index: Integer; Vertex: Integer; X: Double; Y: Double): Integer; stdcall; external EpanetLib; - function ENsetvertices(Index: Integer; var X: Double; var Y: Double; Count: Integer): Integer; stdcall; external EpanetLib; + function ENaddlink(ID: PAnsiChar; LinkType: Integer; FromNode: PAnsiChar; ToNode: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENdeletelink(Index: Integer; ActionCode: Integer): Integer; cdecl; external EpanetLib; + function ENgetlinkindex(ID: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetlinkid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetlinkid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetlinktype(Index: Integer; var Code: Integer): Integer; cdecl; external EpanetLib; + function ENsetlinktype(var Index: Integer; LinkType: Integer; ActionCode: Integer): Integer; cdecl; external EpanetLib; + function ENgetlinknodes(Index: Integer; var Node1: Integer; var Node2: Integer): Integer; cdecl; external EpanetLib; + function ENsetlinknodes(Index: Integer; Node1: Integer; Node2: Integer): Integer; cdecl; external EpanetLib; + function ENgetlinkvalue(Index: Integer; Code: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetlinkvalue(Index: Integer; Code: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENsetpipedata(Index: Integer; Length: Single; Diam: Single; Rough: Single; Mloss:Single): Integer; cdecl; external EpanetLib; + function ENgetlinkvalues(Code: Integer; var X: array of Single): Integer; cdecl; external EpanetLib; + + function ENgetvertexcount(Index: Integer; var Count: Integer): Integer; cdecl; external EpanetLib; + function ENgetvertex(Index: Integer; Vertex: Integer; var X: Double; var Y: Double): Integer; cdecl; external EpanetLib; + function ENsetvertex(Index: Integer; Vertex: Integer; X: Double; Y: Double): Integer; cdecl; external EpanetLib; + function ENsetvertices(Index: Integer; var X: Double; var Y: Double; Count: Integer): Integer; cdecl; external EpanetLib; {Pump Functions} - function ENgetpumptype(LinkIndex: Integer; var PumpType: Integer): Integer; stdcall; external EpanetLib; - function ENgetheadcurveindex(LinkIndex: Integer; var CurveIndex: Integer): Integer; stdcall; external EpanetLib; - function ENsetheadcurveindex(LinkIndex: Integer; CurveIndex: Integer): Integer; stdcall; external EpanetLib; + function ENgetpumptype(LinkIndex: Integer; var PumpType: Integer): Integer; cdecl; external EpanetLib; + function ENgetheadcurveindex(LinkIndex: Integer; var CurveIndex: Integer): Integer; cdecl; external EpanetLib; + function ENsetheadcurveindex(LinkIndex: Integer; CurveIndex: Integer): Integer; cdecl; external EpanetLib; {Pattern Functions} - function ENaddpattern(ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENdeletepattern(Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetpatternindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetpatternid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetpatternid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetpatternlen(Index: Integer; var Len: Integer): Integer; stdcall; external EpanetLib; - function ENgetpatternvalue(Index: Integer; Period: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetpatternvalue(Index: Integer; Period: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENgetaveragepatternvalue(Index: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetpattern(Index: Integer; var F: Single; N: Integer): Integer; stdcall; external EpanetLib; + function ENaddpattern(ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENdeletepattern(Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetpatternindex(ID: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetpatternid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetpatternid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetpatternlen(Index: Integer; var Len: Integer): Integer; cdecl; external EpanetLib; + function ENgetpatternvalue(Index: Integer; Period: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetpatternvalue(Index: Integer; Period: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENgetaveragepatternvalue(Index: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetpattern(Index: Integer; var F: Single; N: Integer): Integer; cdecl; external EpanetLib; {Curve Functions} - function ENaddcurve(ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENdeletecurve(Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetcurveindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetcurveid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetcurveid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetcurvelen(Index: Integer; var Len: Integer): Integer; stdcall; external EpanetLib; - function ENgetcurvetype(Index: Integer; var CurveType: Integer): Integer; stdcall; external EpanetLib; - function ENsetcurvetype(Index: Integer; CurveType: Integer): Integer; stdcall; external EpanetLib; - function ENgetcurvevalue(CurveIndex: Integer; PointIndex: Integer; var X: Single; var Y: Single): Integer; stdcall; external EpanetLib; - function ENsetcurvevalue(CurveIndex: Integer; PointIndex: Integer; X: Single; Y: Single): Integer; stdcall; external EpanetLib; - function ENgetcurve(Index: Integer; ID: PAnsiChar; var N: Integer; var X: Single; var Y: Single): Integer; stdcall; external EpanetLib; - function ENsetcurve(Index: Integer; var X: Single; var Y: Single; N: Integer): Integer; stdcall; external EpanetLib; + function ENaddcurve(ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENdeletecurve(Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetcurveindex(ID: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetcurveid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetcurveid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetcurvelen(Index: Integer; var Len: Integer): Integer; cdecl; external EpanetLib; + function ENgetcurvetype(Index: Integer; var CurveType: Integer): Integer; cdecl; external EpanetLib; + function ENsetcurvetype(Index: Integer; CurveType: Integer): Integer; cdecl; external EpanetLib; + function ENgetcurvevalue(CurveIndex: Integer; PointIndex: Integer; var X: Single; var Y: Single): Integer; cdecl; external EpanetLib; + function ENsetcurvevalue(CurveIndex: Integer; PointIndex: Integer; X: Single; Y: Single): Integer; cdecl; external EpanetLib; + function ENgetcurve(Index: Integer; ID: PAnsiChar; var N: Integer; var X: Single; var Y: Single): Integer; cdecl; external EpanetLib; + function ENsetcurve(Index: Integer; var X: Single; var Y: Single; N: Integer): Integer; cdecl; external EpanetLib; {Control Functions} - function ENaddcontrol(Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENdeletecontrol(Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetcontrol(Index: Integer; var Ctype: Integer; var Link: Integer; var Setting: Single; var Node: Integer; var Level: Single): Integer; stdcall; external EpanetLib; - function ENsetcontrol(Index: Integer; Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single): Integer; stdcall; external EpanetLib; - function ENgetcontrolenabled(Index: Integer; out_enabled: Integer): Integer; stdcall; external EpanetLib; - function ENsetcontrolenabled(Index: Integer; var enabled: Integer): Integer; stdcall; external EpanetLib; + function ENaddcontrol(Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENdeletecontrol(Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetcontrol(Index: Integer; var Ctype: Integer; var Link: Integer; var Setting: Single; var Node: Integer; var Level: Single): Integer; cdecl; external EpanetLib; + function ENsetcontrol(Index: Integer; Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single): Integer; cdecl; external EpanetLib; + function ENgetcontrolenabled(Index: Integer; out_enabled: Integer): Integer; cdecl; external EpanetLib; + function ENsetcontrolenabled(Index: Integer; var enabled: Integer): Integer; cdecl; external EpanetLib; {Rule-Based Control Functions} - function ENaddrule(Rule: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENdeleterule(Index: Integer): Integer; stdcall; external EpanetLib; + function ENaddrule(Rule: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENdeleterule(Index: Integer): Integer; cdecl; external EpanetLib; function ENgetrule(Index: Integer; var Npremises: Integer; var NthenActions: Integer; - var NelseActions: Integer; var Priority: Single): Integer; stdcall; external EpanetLib; - function ENgetruleID(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetrulepriority(Index: Integer; Priority: Single): Integer; stdcall; external EpanetLib; + var NelseActions: Integer; var Priority: Single): Integer; cdecl; external EpanetLib; + function ENgetruleID(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetrulepriority(Index: Integer; Priority: Single): Integer; cdecl; external EpanetLib; function ENgetpremise(RuleIndex: Integer; PremiseIndex: Integer; var LogOp: Integer; var ObjType: Integer; var ObjIndex: Integer; var Param: Integer; var RelOp: Integer; - var Status: Integer; var Value: Single): Integer; stdcall; external EpanetLib; + var Status: Integer; var Value: Single): Integer; cdecl; external EpanetLib; function ENsetpremise(RuleIndex: Integer; PremiseIndex: Integer; LogOp: Integer; ObjType: Integer; - ObjIndex: Integer; Param: Integer; RelOp: Integer; Status: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENsetpremiseindex(RuleIndex: Integer; PremiseIndex: Integer; ObjIndex: Integer): Integer; stdcall; external EpanetLib; - function ENsetpremisestatus(RuleIndex: Integer; PremiseIndex: Integer; Status: Integer): Integer; stdcall; external EpanetLib; - function ENsetpremisevalue(RuleIndex: Integer; PremiseIndex: Integer; Value: Single): Integer; stdcall; external EpanetLib; + ObjIndex: Integer; Param: Integer; RelOp: Integer; Status: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENsetpremiseindex(RuleIndex: Integer; PremiseIndex: Integer; ObjIndex: Integer): Integer; cdecl; external EpanetLib; + function ENsetpremisestatus(RuleIndex: Integer; PremiseIndex: Integer; Status: Integer): Integer; cdecl; external EpanetLib; + function ENsetpremisevalue(RuleIndex: Integer; PremiseIndex: Integer; Value: Single): Integer; cdecl; external EpanetLib; function ENgetthenaction(RuleIndex: Integer; ActionIndex: Integer; var LinkIndex: Integer; - var Status: Integer; var Setting: Single): Integer; stdcall; external EpanetLib; + var Status: Integer; var Setting: Single): Integer; cdecl; external EpanetLib; function ENsetthenaction(RuleIndex: Integer; ActionIndex: Integer; LinkIndex: Integer; - Status: Integer; Setting: Single): Integer; stdcall; external EpanetLib; + Status: Integer; Setting: Single): Integer; cdecl; external EpanetLib; function ENgetelseaction(RuleIndex: Integer; ActionIndex: Integer; var LinkIndex: Integer; - var Status: Integer; var Setting: Single): Integer; stdcall; external EpanetLib; + var Status: Integer; var Setting: Single): Integer; cdecl; external EpanetLib; function ENsetelseaction(RuleIndex: Integer; ActionIndex: Integer; LinkIndex: Integer; - Status: Integer; Setting: Single): Integer; stdcall; external EpanetLib; - function ENgetruleenabled(Index: Integer; out_enabled: Integer): Integer; stdcall; external EpanetLib; - function ENsetruleenabled(Index: Integer; var enabled: Integer): Integer; stdcall; external EpanetLib; + Status: Integer; Setting: Single): Integer; cdecl; external EpanetLib; + function ENgetruleenabled(Index: Integer; out_enabled: Integer): Integer; cdecl; external EpanetLib; + function ENsetruleenabled(Index: Integer; var enabled: Integer): Integer; cdecl; external EpanetLib; implementation diff --git a/include/epanet2.vb b/include/epanet2.vb index 73e46511..a7b31799 100644 --- a/include/epanet2.vb +++ b/include/epanet2.vb @@ -4,7 +4,7 @@ 'Declarations of functions in the EPANET PROGRAMMERs TOOLKIT '(EPANET2.DLL) for use with VB.Net. -'Last updated on 09/28/2023 +'Last updated on 06/23/2024 Imports System.Runtime.InteropServices Imports System.Text @@ -42,6 +42,9 @@ Public Const EN_CANOVERFLOW = 26 Public Const EN_DEMANDDEFICIT = 27 Public Const EN_NODE_INCONTROL = 28 Public Const EN_EMITTERFLOW = 29 +Public Const EN_LEAKAGEFLOW = 30 +Public Const EN_DEMANDFLOW = 31 +Public Const EN_FULLDEMAND = 32 Public Const EN_DIAMETER = 0 ' Link parameters Public Const EN_LENGTH = 1 @@ -69,6 +72,9 @@ Public Const EN_PUMP_EPAT = 22 Public Const EN_LINK_INCONTROL = 23 Public Const EN_GPV_CURVE = 24 Public Const EN_PCV_CURVE = 25 +Public Const EN_LEAK_AREA = 26 +Public Const EN_LEAK_EXPAN = 27 +Public Const EN_LINK_LEAKAGE = 28 Public Const EN_DURATION = 0 ' Time parameters Public Const EN_HYDSTEP = 1 @@ -93,6 +99,7 @@ Public Const EN_MAXFLOWCHANGE = 3 Public Const EN_MASSBALANCE = 4 Public Const EN_DEFICIENTNODES = 5 Public Const EN_DEMANDREDUCTION = 6 +Public Const EN_LEAKAGELOSS = 7 Public Const EN_NODE = 0 ' Component types Public Const EN_LINK = 1 @@ -309,6 +316,7 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENgeterror Lib "epanet2.dll" (ByVal errcode As Int32, ByVal errmsg As String, ByVal maxLen As Int32) As Int32 Declare Function ENgetstatistic Lib "epanet2.dll" (ByVal type_ As Int32, ByRef value As Single) As Int32 Declare Function ENgetresultindex Lib "epanet2.dll" (ByVal type_ As Int32, ByVal index As Int32, ByRef value As Int32) As Int32 + Declare Function ENtimetonextevent Lib "epanet2.dll" (eventType As Int32, duration As Int32, elementIndex As Int32) As Int32 'Analysis Options Functions Declare Function ENgetoption Lib "epanet2.dll" (ByVal option As Int32, value As Single) As Int32 @@ -334,7 +342,8 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENsettankdata Lib "epanet2.dll" (ByVal index As Int32, ByVal elev As Single, ByVal initlvl As Single, ByVal minlvl As Single, ByVal maxlvl As Single, ByVal diam As Single, ByVal minvol As Single, ByVal volcurve As String) As Int32 Declare Function ENgetcoord Lib "epanet2.dll" (ByVal index As Int32, x As Double, y As Double) As Int32 Declare Function ENsetcoord Lib "epanet2.dll" (ByVal index As Int32, ByVal x As Double, ByVal y As Double) As Int32 - + Declare Function ENgetnodevalues Lib "epanet2.dll" (ByVal property as Int32, values as Any) As Int32 + 'Nodal Demand Functions Declare Function ENgetdemandmodel Lib "epanet2.dll" (type_ As Int32, pmin As Single, preq As Single, pexp As Single) As Int32 Declare Function ENsetdemandmodel Lib "epanet2.dll" (ByVal type_ As Int32, ByVal pmin As Single, ByVal preq As Single, ByVal pexp As Single) As Int32 @@ -366,6 +375,7 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENgetvertex Lib "epanet2.dll" (ByVal index As Int32, ByVal vertex As Int32, x As Double, y As Double) As Int32 Declare Function ENsetvertex Lib "epanet2.dll" (ByVal index As Int32, ByVal vertex As Int32, ByVal x As Double, ByVal y As Double) As Int32 Declare Function ENsetvertices Lib "epanet2.dll" (ByVal index As Int32, xCoords As Any, yCoords As Any, ByVal count As Int32) As Int32 + Declare Function ENgetlinkvalues Lib "epanet2.dll" (ByVal property as Int32, values as Any) As Int32 'Pump Functions Declare Function ENgetheadcurveindex Lib "epanet2.dll" (ByVal linkIndex As Int32, curveIndex As Int32) As Int32 diff --git a/include/epanet2_enums.h b/include/epanet2_enums.h index 557f48fe..4fc7755e 100644 --- a/include/epanet2_enums.h +++ b/include/epanet2_enums.h @@ -65,7 +65,10 @@ typedef enum { EN_CANOVERFLOW = 26, //!< Tank can overflow (= 1) or not (= 0) EN_DEMANDDEFICIT = 27,//!< Amount that full demand is reduced under PDA (read only) EN_NODE_INCONTROL = 28, //!< Is present in any simple or rule-based control (= 1) or not (= 0) - EN_EMITTERFLOW = 29 //!< Current emitter flow (read only) + EN_EMITTERFLOW = 29, //!< Current emitter flow (read only) + EN_LEAKAGEFLOW = 30, //!< Current leakage flow (read only) + EN_DEMANDFLOW = 31, //!< Current consumer demand delivered (read only) + EN_FULLDEMAND = 32 //!< Current consumer demand requested (read only) } EN_NodeProperty; /// Link properties @@ -99,7 +102,10 @@ typedef enum { EN_PUMP_EPAT = 22, //!< Pump energy price time pattern index EN_LINK_INCONTROL = 23, //!< Is present in any simple or rule-based control (= 1) or not (= 0) EN_GPV_CURVE = 24, //!< GPV head loss v. flow curve index - EN_PCV_CURVE = 25 //!< PCV loss coeff. curve index + EN_PCV_CURVE = 25, //!< PCV loss coeff. curve index + EN_LEAK_AREA = 26, //!< Pipe leak area (sq mm per 100 length units) + EN_LEAK_EXPAN = 27, //!< Leak expansion rate (sq mm per unit of pressure head) + EN_LINK_LEAKAGE = 28 //!< Current leakage rate (read only) } EN_LinkProperty; /// Time parameters @@ -152,7 +158,8 @@ typedef enum { EN_MAXFLOWCHANGE = 3, //!< Largest flow change in links EN_MASSBALANCE = 4, //!< Cumulative water quality mass balance ratio EN_DEFICIENTNODES = 5, //!< Number of pressure deficient nodes - EN_DEMANDREDUCTION = 6 //!< % demand reduction at pressure deficient nodes + EN_DEMANDREDUCTION = 6, //!< % demand reduction at pressure deficient nodes + EN_LEAKAGELOSS = 7 //!< % flow lost to system leakage } EN_AnalysisStatistic; /// Types of network objects diff --git a/src/enumstxt.h b/src/enumstxt.h index dbb4c74c..65962a6d 100755 --- a/src/enumstxt.h +++ b/src/enumstxt.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 02/05/2023 + Last Updated: 06/24/2024 ****************************************************************************** */ @@ -127,7 +127,7 @@ char *SectTxt[] = {s_TITLE, s_JUNCTIONS, s_RESERVOIRS, s_REACTIONS, s_MIXING, s_REPORT, s_TIMES, s_OPTIONS, s_COORDS, s_VERTICES, s_LABELS, s_BACKDROP, - s_TAGS, s_END, + s_TAGS, s_LEAKAGE, s_END, NULL}; char *Fldname[] = {t_ELEV, t_DEMAND, t_HEAD, diff --git a/src/epanet.c b/src/epanet.c index 991024cd..7e68c9d3 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 09/28/2023 + Last Updated: 06/26/2024 ****************************************************************************** */ @@ -510,6 +510,10 @@ int DLLEXPORT EN_initH(EN_Project p, int initFlag) return errcode; } } + + // Open pipe leakage modeling system + errcode = openleakage(p); + if (errcode) return errcode; // Initialize hydraulics solver inithyd(p, fflag); @@ -564,7 +568,11 @@ int DLLEXPORT EN_closeH(EN_Project p) */ { if (!p->Openflag) return 102; - if (p->hydraul.OpenHflag) closehyd(p); + if (p->hydraul.OpenHflag) + { + closeleakage(p); + closehyd(p); + } p->hydraul.OpenHflag = FALSE; return 0; } @@ -1044,6 +1052,9 @@ int DLLEXPORT EN_getstatistic(EN_Project p, int type, double *value) case EN_DEMANDREDUCTION: *value = p->hydraul.DemandReduction; break; + case EN_LEAKAGELOSS: + *value = p->hydraul.LeakageLoss; + break; case EN_MASSBALANCE: *value = p->quality.MassBalance.ratio; break; @@ -1864,8 +1875,10 @@ int DLLEXPORT EN_addnode(EN_Project p, const char *id, int nodeType, int *index) hyd->NodeDemand = (double *)realloc(hyd->NodeDemand, size); qual->NodeQual = (double *)realloc(qual->NodeQual, size); hyd->NodeHead = (double *)realloc(hyd->NodeHead, size); - hyd->DemandFlow = (double *)realloc(hyd->DemandFlow, size); + hyd->FullDemand = (double *)realloc(hyd->FullDemand, size); hyd->EmitterFlow = (double *)realloc(hyd->EmitterFlow, size); + hyd->LeakageFlow = (double *)realloc(hyd->LeakageFlow, size); + hyd->DemandFlow = (double *)realloc(hyd->DemandFlow, size); // Actions taken when a new Junction is added if (nodeType == EN_JUNCTION) @@ -2256,7 +2269,7 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val Ucf[VOLUME]; break; - case EN_DEMAND: + case EN_DEMAND: // Consumer Demand + Emitter Flow + Leakage Flow v = hyd->NodeDemand[index] * Ucf[FLOW]; break; @@ -2336,11 +2349,10 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val case EN_DEMANDDEFICIT: if (index > nJuncs) return 0; - // After an analysis, DemandFlow contains node's required demand - // while NodeDemand contains delivered demand + emitter flow - if (hyd->DemandFlow[index] < 0.0) return 0; - v = (hyd->DemandFlow[index] - - (hyd->NodeDemand[index] - hyd->EmitterFlow[index])) * Ucf[FLOW]; + // FullDemand contains node's required consumer demand + // while DemandFlow contains delivered consumer demand + if (hyd->FullDemand[index] <= 0.0) return 0; + v = (hyd->FullDemand[index] - hyd->DemandFlow[index]) * Ucf[FLOW]; break; case EN_NODE_INCONTROL: @@ -2350,6 +2362,18 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val case EN_EMITTERFLOW: v = hyd->EmitterFlow[index] * Ucf[FLOW]; break; + + case EN_LEAKAGEFLOW: + v = hyd->LeakageFlow[index] * Ucf[FLOW]; + break; + + case EN_DEMANDFLOW: // Consumer demand delivered + v = hyd->DemandFlow[index] * Ucf[FLOW]; + break; + + case EN_FULLDEMAND: // Consumer demand requested + v = hyd->FullDemand[index] * Ucf[FLOW]; + break; default: return 251; @@ -3352,6 +3376,8 @@ int DLLEXPORT EN_addlink(EN_Project p, const char *id, int linkType, } link->Kb = 0; link->Kw = 0; + link->LeakArea = 0; + link->LeakExpan = 0; link->R = 0; link->Rc = 0; link->Rpt = 0; @@ -3923,6 +3949,18 @@ int DLLEXPORT EN_getlinkvalue(EN_Project p, int index, int property, double *val v = (double)incontrols(p, LINK, index); break; + case EN_LEAK_AREA: + v = Link[index].LeakArea * Ucf[LENGTH]; + break; + + case EN_LEAK_EXPAN: + v = Link[index].LeakExpan * Ucf[LENGTH]; + break; + + case EN_LINK_LEAKAGE: + v = findlinkleakage(p, index) * Ucf[FLOW]; + break; + default: return 251; } @@ -4163,6 +4201,16 @@ int DLLEXPORT EN_setlinkvalue(EN_Project p, int index, int property, double valu } break; + case EN_LEAK_AREA: // leak area in sq mm per 100 pipe length units + if (value < 0.0) return 211; + Link[index].LeakArea = value / Ucf[LENGTH]; + break; + + case EN_LEAK_EXPAN: // leak area expansion slope (sq mm per unit of head) + if (value < 0.0) return 211; + Link[index].LeakExpan = value / Ucf[LENGTH]; + break; + default: return 251; } diff --git a/src/flowbalance.c b/src/flowbalance.c new file mode 100644 index 00000000..2c0efa66 --- /dev/null +++ b/src/flowbalance.c @@ -0,0 +1,191 @@ +/* + ****************************************************************************** + Project: OWA EPANET + Version: 2.3 + Module: flowbalance.c + Description: computes components of network's flow balance + Authors: see AUTHORS + Copyright: see AUTHORS + License: see LICENSE + Last Updated: 06/26/2024 + ****************************************************************************** +*/ + +#include "types.h" + +// Exported functions (declared in funcs.h) +//void startflowbalance(Project *); +//void updateflowbalance(Project *, long); +//void endflowbalance(Project *); + +void startflowbalance(Project *pr) +/* +**------------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: initializes components of the network's flow balance. +**------------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + hyd->FlowBalance.totalInflow = 0.0; + hyd->FlowBalance.totalOutflow = 0.0; + hyd->FlowBalance.consumerDemand = 0.0; + hyd->FlowBalance.emitterDemand = 0.0; + hyd->FlowBalance.leakageDemand = 0.0; + hyd->FlowBalance.deficitDemand = 0.0; + hyd->FlowBalance.storageDemand = 0.0; + hyd->FlowBalance.ratio = 0.0; +} + +void updateflowbalance(Project *pr, long hstep) +/* +**------------------------------------------------------------------- +** Input: hstep = time step (sec) +** Output: none +** Purpose: updates components of the system flow balance. +**------------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + Times *time = &pr->times; + + int i, j; + double v, dt, deficit, fullDemand; + SflowBalance flowBalance; + + // Determine current time interval in seconds + if (time->Dur == 0) dt = 1.0; + else if (time->Htime < time->Dur) + { + dt = (double) hstep; + } + else return; + + // Initialize local flow balance + flowBalance.totalInflow = 0.0; + flowBalance.totalOutflow = 0.0; + flowBalance.consumerDemand = 0.0; + flowBalance.emitterDemand = 0.0; + flowBalance.leakageDemand = 0.0; + flowBalance.deficitDemand = 0.0; + flowBalance.storageDemand = 0.0; + fullDemand = 0.0; + + // Initialize demand deficiency & leakage loss + hyd->DeficientNodes = 0; + hyd->DemandReduction = 0.0; + hyd->LeakageLoss = 0.0; + + // Examine each junction node + for (i = 1; i <= net->Njuncs; i++) + { + // Accumulate consumer demand flow + v = hyd->DemandFlow[i]; + if (v < 0.0) + flowBalance.totalInflow += (-v); + else + { + fullDemand += hyd->FullDemand[i]; + flowBalance.consumerDemand += v; + flowBalance.totalOutflow += v; + } + + // Accumulate emitter and leakage flow + v = hyd->EmitterFlow[i]; + flowBalance.emitterDemand += v; + flowBalance.totalOutflow += v; + v = hyd->LeakageFlow[i]; + flowBalance.leakageDemand += v; + flowBalance.totalOutflow += v; + + // Accumulate demand deficit flow + if (hyd->DemandModel == PDA && hyd->FullDemand[i] > 0.0) + { + deficit = hyd->FullDemand[i] - hyd->DemandFlow[i]; + if (deficit > TINY) + { + hyd->DeficientNodes++; + flowBalance.deficitDemand += deficit; + } + } + } + + // Examine each tank/reservoir node + for (j = 1; j <= net->Ntanks; j++) + { + 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) + { + if (v >= 0.0) + flowBalance.totalOutflow += v; + else + flowBalance.totalInflow += (-v); + } + + // For tank under extended period analysis + 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; + + // Update flow balance for entire run + hyd->FlowBalance.totalInflow += flowBalance.totalInflow * dt; + hyd->FlowBalance.totalOutflow += flowBalance.totalOutflow * dt; + hyd->FlowBalance.consumerDemand += flowBalance.consumerDemand * dt; + hyd->FlowBalance.emitterDemand += flowBalance.emitterDemand * dt; + hyd->FlowBalance.leakageDemand += flowBalance.leakageDemand * dt; + hyd->FlowBalance.deficitDemand += flowBalance.deficitDemand * dt; + hyd->FlowBalance.storageDemand += flowBalance.storageDemand * dt; +} + +void endflowbalance(Project *pr) +/* +**------------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: finalizes components of the system flow balance. +**------------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + Times *time = &pr->times; + + double seconds, qin, qout, qstor, r; + + if (time->Htime > 0) + seconds = time->Htime; + else + seconds = 1.0; + hyd->FlowBalance.totalInflow /= seconds; + hyd->FlowBalance.totalOutflow /= seconds; + hyd->FlowBalance.consumerDemand /= seconds; + hyd->FlowBalance.emitterDemand /= seconds; + hyd->FlowBalance.leakageDemand /= seconds; + hyd->FlowBalance.deficitDemand /= seconds; + hyd->FlowBalance.storageDemand /= seconds; + + qin = hyd->FlowBalance.totalInflow; + qout = hyd->FlowBalance.totalOutflow; + qstor = hyd->FlowBalance.storageDemand; + if (qstor > 0.0) + qout += qstor; + else + qin -= qstor; + if (qin == qout) + r = 1.0; + else if (qin > 0.0) + r = qout / qin; + else + r = 0.0; + hyd->FlowBalance.ratio = r; +} diff --git a/src/funcs.h b/src/funcs.h index 8910c361..87a6be24 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 09/28/2023 + Last Updated: 06/26/2024 ****************************************************************************** */ #ifndef FUNCS_H @@ -100,6 +100,7 @@ int controldata(Project *); int energydata(Project *); int sourcedata(Project *); int emitterdata(Project *); +int leakagedata(Project *); int qualdata(Project *); int reactdata(Project *); int mixingdata(Project *); @@ -142,6 +143,7 @@ void writecontrolaction(Project *, int, int); void writeruleaction(Project *, int, char *); int writehydwarn(Project *, int,double); void writehyderr(Project *, int); +void writeflowbalance(Project *); void writemassbalance(Project *); void writetime(Project *, char *); char *clocktime(char *, long); @@ -195,4 +197,19 @@ int savefinaloutput(Project *); int saveinpfile(Project *, const char *); +// ------- LEAKAGE.C -------------------- + +int openleakage(Project *); +void closeleakage(Project *); +double findlinkleakage(Project *, int); +void leakagecoeffs(Project *); +double leakageflowchange(Project *, int); +int leakagehasconverged(Project *); + +// ------- FLOWBALANCE.C----------------- + +void startflowbalance(Project *); +void updateflowbalance(Project *, long); +void endflowbalance(Project *); + #endif diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index c00d2f9f..5be8c403 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 03/29/2023 + Last Updated: 06/15/2024 ****************************************************************************** */ @@ -310,6 +310,7 @@ void matrixcoeffs(Project *pr) linkcoeffs(pr); emittercoeffs(pr); demandcoeffs(pr); + if (hyd->HasLeakage) leakagecoeffs(pr); // Update nodal flow balances with demands and add onto r.h.s. coeffs. nodecoeffs(pr); @@ -606,7 +607,7 @@ void demandheadloss(Project *pr, int i, double dp, double n, Hydraul *hyd = &pr->hydraul; double d = hyd->DemandFlow[i]; - double dfull = hyd->NodeDemand[i]; + double dfull = hyd->FullDemand[i]; double r = d / dfull; // Evaluate inverted demand function diff --git a/src/hydraul.c b/src/hydraul.c index 39769212..ae40ad63 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 09/28/2023 + Last Updated: 06/26/2024 ****************************************************************************** */ @@ -63,7 +63,7 @@ int openhyd(Project *pr) // Allocate memory for hydraulic variables ERRCODE(allocmatrix(pr)); - + // Check for unconnected nodes ERRCODE(unlinked(pr)); @@ -107,8 +107,10 @@ void inithyd(Project *pr, int initflag) hyd->OldStatus[net->Nlinks+i] = TEMPCLOSED; } - // Initialize emitter flows + // Initialize node outflows + memset(hyd->DemandFlow,0,(net->Nnodes+1)*sizeof(double)); memset(hyd->EmitterFlow,0,(net->Nnodes+1)*sizeof(double)); + memset(hyd->LeakageFlow,0,(net->Nnodes+1)*sizeof(double)); for (i = 1; i <= net->Nnodes; i++) { net->Node[i].ResultIndex = i; @@ -161,6 +163,9 @@ void inithyd(Project *pr, int initflag) pump->Energy.CurrentPower = 0.0; pump->Energy.CurrentEffic = 0.0; } + + // Initialize flow balance + startflowbalance(pr); // Re-position hydraulics file if (pr->outfile.Saveflag) @@ -253,6 +258,9 @@ int nexthyd(Project *pr, long *tstep) // Accumulate pumping energy if (time->Dur == 0) addenergy(pr,0); else if (time->Htime < time->Dur) addenergy(pr,hydstep); + + // Update flow balance + updateflowbalance(pr, hydstep); // More time remains - update current time if (time->Htime < time->Dur) @@ -267,6 +275,8 @@ int nexthyd(Project *pr, long *tstep) // No more time remains - force completion of analysis else { + endflowbalance(pr); + if (pr->report.Statflag) writeflowbalance(pr); time->Htime++; if (pr->quality.OpenQflag) time->Qtime++; } @@ -495,7 +505,7 @@ void demands(Project *pr) if (djunc > 0.0) hyd->Dsystem += djunc; sum += djunc; } - hyd->NodeDemand[i] = sum; + hyd->FullDemand[i] = sum; // Initialize pressure dependent demand hyd->DemandFlow[i] = sum; @@ -1147,4 +1157,3 @@ void resetpumpflow(Project *pr, int i) if (pump->Ptype == CONST_HP) pr->hydraul.LinkFlow[i] = pump->Q0; } - diff --git a/src/hydsolver.c b/src/hydsolver.c index 25c6752c..d88823fb 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -8,7 +8,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 02/14/2022 + Last Updated: 06/15/2024 ****************************************************************************** */ @@ -47,6 +47,7 @@ static double newflows(Project *, Hydbalance *); static void newlinkflows(Project *, Hydbalance *, double *, double *); static void newemitterflows(Project *, Hydbalance *, double *, double *); static void newdemandflows(Project *, Hydbalance *, double *, double *); +static void newleakageflows(Project *, Hydbalance *, double *, double *); static void checkhydbalance(Project *, Hydbalance *); static int hasconverged(Project *, double *, Hydbalance *); @@ -93,7 +94,6 @@ int hydsolve(Project *pr, int *iter, double *relerr) int valveChange; // Valve status change flag int statChange; // Non-valve status change flag Hydbalance hydbal; // Hydraulic balance errors - double fullDemand; // Full demand for a node (cfs) // Initialize status checking & relaxation factor nextcheck = hyd->CheckFreq; @@ -195,12 +195,12 @@ int hydsolve(Project *pr, int *iter, double *relerr) errcode = 110; } - // Store actual junction outflow in NodeDemand & full demand in DemandFlow + // Save total outflow (NodeDemand) at each junction for (i = 1; i <= net->Njuncs; i++) { - fullDemand = hyd->NodeDemand[i]; - hyd->NodeDemand[i] = hyd->DemandFlow[i] + hyd->EmitterFlow[i]; - hyd->DemandFlow[i] = fullDemand; + hyd->NodeDemand[i] = hyd->DemandFlow[i] + + hyd->EmitterFlow[i] + + hyd->LeakageFlow[i]; } // Save convergence info @@ -381,6 +381,7 @@ double newflows(Project *pr, Hydbalance *hbal) newlinkflows(pr, hbal, &qsum, &dqsum); newemitterflows(pr, hbal, &qsum, &dqsum); newdemandflows(pr, hbal, &qsum, &dqsum); + if (hyd->HasLeakage) newleakageflows(pr, hbal, &qsum, &dqsum); // Return ratio of total flow corrections to total flow if (qsum > hyd->Hacc) return (dqsum / qsum); @@ -514,6 +515,45 @@ void newemitterflows(Project *pr, Hydbalance *hbal, double *qsum, } +void newleakageflows(Project *pr, Hydbalance *hbal, double *qsum, + double *dqsum) +/* +**---------------------------------------------------------------- +** Input: hbal = ptr. to hydraulic balance information +** qsum = sum of current system flows +** dqsum = sum of system flow changes +** Output: updates hbal, qsum and dqsum +** Purpose: updates nodal leakage flows after new nodal heads computed +**---------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + int i; + double dq; + + for (i = 1; i <= net->Njuncs; i++) + { + // Update leakage flow at node i + dq = leakageflowchange(pr, i); + if (dq == 0.0) continue; + + // Update system flow summation + *qsum += ABS(hyd->LeakageFlow[i]); + *dqsum += ABS(dq); + + // Update identity of element with max. flow change + if (ABS(dq) > hbal->maxflowchange) + { + hbal->maxflowchange = ABS(dq); + hbal->maxflownode = i; + hbal->maxflowlink = -1; + } + } +} + + void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum) /* **---------------------------------------------------------------- @@ -546,7 +586,7 @@ void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum) for (i = 1; i <= net->Njuncs; i++) { // Skip junctions with no positive demand - if (hyd->NodeDemand[i] <= 0.0) continue; + if (hyd->FullDemand[i] <= 0.0) continue; // Find change in demand flow (see hydcoeffs.c) demandheadloss(pr, i, dp, n, &hloss, &hgrad); @@ -555,8 +595,8 @@ void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum) dq *= hyd->RelaxFactor; // Prevent a flow change greater than full demand - if (fabs(dq) > 0.4 * hyd->NodeDemand[i]) - dq = 0.4 * SGN(dq) * hyd->NodeDemand[i]; + if (fabs(dq) > 0.4 * hyd->FullDemand[i]) + dq = 0.4 * SGN(dq) * hyd->FullDemand[i]; hyd->DemandFlow[i] -= dq; // Update system flow summation @@ -641,11 +681,15 @@ int hasconverged(Project *pr, double *relerr, Hydbalance *hbal) if (hyd->FlowChangeLimit > 0.0 && hbal->maxflowchange > hyd->FlowChangeLimit) return 0; + // Check for node leakage convergence + if (hyd->HasLeakage && !leakagehasconverged(pr)) return 0; + // Check for pressure driven analysis convergence if (hyd->DemandModel == PDA) return pdaconverged(pr); return 1; } + int pdaconverged(Project *pr) /* **-------------------------------------------------------------- @@ -659,47 +703,32 @@ int pdaconverged(Project *pr) Hydraul *hyd = &pr->hydraul; const double QTOL = 0.0001; // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm) - int i, converged = 1; - double totalDemand = 0.0, totalReduction = 0.0; + int i; double dp = hyd->Preq - hyd->Pmin; double p, q, r; - hyd->DeficientNodes = 0; - hyd->DemandReduction = 0.0; - - // Add up number of junctions with demand deficits + // Examine each network junction for (i = 1; i <= pr->network.Njuncs; i++) { // Skip nodes whose required demand is non-positive - if (hyd->NodeDemand[i] <= 0.0) continue; + if (hyd->FullDemand[i] <= 0.0) continue; // Evaluate demand equation at current pressure solution p = hyd->NodeHead[i] - pr->network.Node[i].El; if (p <= hyd->Pmin) q = 0.0; else if (p >= hyd->Preq) - q = hyd->NodeDemand[i]; + q = hyd->FullDemand[i]; else { r = (p - hyd->Pmin) / dp; - q = hyd->NodeDemand[i] * pow(r, hyd->Pexp); + q = hyd->FullDemand[i] * pow(r, hyd->Pexp); } // Check if demand has not converged - if (fabs(q - hyd->DemandFlow[i]) > QTOL) - converged = 0; - - // Accumulate total required demand and demand deficit - if (hyd->DemandFlow[i] + QTOL < hyd->NodeDemand[i]) - { - hyd->DeficientNodes++; - totalDemand += hyd->NodeDemand[i]; - totalReduction += hyd->NodeDemand[i] - hyd->DemandFlow[i]; - } - } - if (totalDemand > 0.0) - hyd->DemandReduction = totalReduction / totalDemand * 100.0; - return converged; + if (fabs(q - hyd->DemandFlow[i]) > QTOL) return 0; + } + return 1; } diff --git a/src/inpfile.c b/src/inpfile.c index ec786b26..9945a2ae 100644 --- a/src/inpfile.c +++ b/src/inpfile.c @@ -7,7 +7,7 @@ Description: saves network data to an EPANET formatted text file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 05/11/2024 +Last Updated: 06/18/2024 ****************************************************************************** */ @@ -369,6 +369,19 @@ int saveinpfile(Project *pr, const char *fname) fprintf(f, "\n %-31s\t%-14.6f", node->ID, ke); } + // Write [LEAKAGE] section + fprintf(f, "\n\n"); + fprintf(f, s_LEAKAGE); + fprintf(f, "\n;;%-31s\t%-14s\t%-14s", + "Pipe", "Leak Area", "Leak Expansion"); + for (i = 1; i <= net->Nlinks; i++) + { + link = &net->Link[i]; + if (link->LeakArea == 0.0 && link->LeakExpan == 0.0) continue; + fprintf(f, "\n %-31s %14.6f %14.6f", link->ID, + link->LeakArea * pr->Ucf[LENGTH], link->LeakExpan * pr->Ucf[LENGTH]); + } + // Write [STATUS] section fprintf(f, "\n\n"); fprintf(f, s_STATUS); @@ -584,7 +597,7 @@ int saveinpfile(Project *pr, const char *fname) fprintf(f, "\n\n"); fprintf(f, s_REACTIONS); - fprintf(f, "\n;;%-9s\t%-31s\t%-12s", "Type", "Pipe/Tank", "Coefficient"); + fprintf(f, "\n;%-9s\t%-31s\t%-12s", "Type", "Pipe/Tank", "Coefficient"); // Pipe-specific parameters for (i = 1; i <= net->Nlinks; i++) diff --git a/src/input1.c b/src/input1.c index ee40d1ec..1074decd 100644 --- a/src/input1.c +++ b/src/input1.c @@ -7,7 +7,7 @@ Description: retrieves network data from an EPANET input file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 09/28/2023 +Last Updated: 06/15/2024 ****************************************************************************** */ @@ -593,6 +593,10 @@ void convertunits(Project *pr) // Convert units on reaction coeffs. link->Kb /= SECperDAY; link->Kw /= SECperDAY; + + // Convert leakage parameters + link->LeakArea /= pr->Ucf[LENGTH]; + link->LeakExpan /= pr->Ucf[LENGTH]; } else if (link->Type == PUMP) diff --git a/src/input2.c b/src/input2.c index c5af2bae..f1ce8bfe 100644 --- a/src/input2.c +++ b/src/input2.c @@ -310,6 +310,7 @@ int newline(Project *pr, int sect, char *line) else return 0; case _SOURCES: return (sourcedata(pr)); case _EMITTERS: return (emitterdata(pr)); + case _LEAKAGE: return (leakagedata(pr)); case _QUALITY: return (qualdata(pr)); case _STATUS: return (statusdata(pr)); case _ROUGHNESS: return (0); diff --git a/src/input3.c b/src/input3.c index 1d4f03a7..296e752e 100644 --- a/src/input3.c +++ b/src/input3.c @@ -387,6 +387,8 @@ int pipedata(Project *pr) link->Km = 0.0; link->Kb = MISSING; link->Kw = MISSING; + link->LeakArea = 0.0; + link->LeakExpan = 0.0; link->Type = PIPE; link->Status = OPEN; link->Rpt = 0; @@ -494,6 +496,8 @@ int pumpdata(Project *pr) link->Km = 0.0; link->Kb = 0.0; link->Kw = 0.0; + link->LeakArea = 0.0; + link->LeakExpan = 0.0; link->Type = PUMP; link->Status = OPEN; link->Rpt = 0; @@ -613,6 +617,8 @@ int valvedata(Project *pr) link->Km = 0.0; link->Kb = 0.0; link->Kw = 0.0; + link->LeakArea = 0.0; + link->LeakExpan = 0.0; link->Type = type; link->Status = ACTIVE; link->Rpt = 0; @@ -1118,6 +1124,41 @@ int emitterdata(Project *pr) return 0; } +int leakagedata(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns error code +** Purpose: processes link leakage data +** Format: +** [LEAKAGE] +** link C1 C2 +**-------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Parser *parser = &pr->parser; + + int j, // Link index + n; // # data items + double c1, c2; // Flow coeff. + + // Check that link exists & is a pipe + n = parser->Ntokens; + if (n < 3) return 201; + if ((j = findlink(net, parser->Tok[0])) == 0) return setError(parser, 0, 203); + if (net->Link[j].Type > PIPE) return 0; + + // Parse leakage coeffs. + if (!getfloat(parser->Tok[1], &c1)) return setError(parser, 1, 202); + if (c1 < 0.0) return setError(parser, 1, 209); + if (!getfloat(parser->Tok[2], &c2)) return setError(parser, 2, 202); + if (c2 < 0.0) return setError(parser, 1, 209); + net->Link[j].LeakArea = c1; + net->Link[j].LeakExpan = c2; + return 0; +} + int qualdata(Project *pr) /* **-------------------------------------------------------------- diff --git a/src/leakage.c b/src/leakage.c new file mode 100644 index 00000000..cb5cad71 --- /dev/null +++ b/src/leakage.c @@ -0,0 +1,531 @@ +/* + ****************************************************************************** + Project: OWA EPANET + Version: 2.3 + Module: leakage.c + Description: models additional nodal demands due to pipe leaks. + Authors: see AUTHORS + Copyright: see AUTHORS + License: see LICENSE + Last Updated: 06/14/2024 + ****************************************************************************** +*/ +/* +This module uses the FAVAD (Fixed and Variable Discharge) equation to model +leaky pipes: + + Q = Co * L * (Ao + m * H) * sqrt(H) + +where Q = 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. + +The inverted form of this equation is used to model the leakage demand from +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, +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. + +In implementing this model, the pipe property "LeakArea" represents Ao in +sq. mm per 100 units of pipe length and "LeakExpan" represents m in sq. mm +per unit of pressure head. + +*/ +#include +#include + +#include "types.h" +#include "funcs.h" + +// Exported functions (declared in funcs.h) +//int openleakage(Project *); +//void closeleakage(Project *); +//double findlinkleakage(Project *, int); +//void leakagecoeffs(Project *); +//double leakageflowchange(Project *, int); +//int leakagehasconverged(Project *); + +// Local functions +static int check_for_leakage(Project *pr); +static int create_leakage_objects(Project *pr); +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); + + +int openleakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: returns an error code +** Purpose: opens the pipe leakage modeling system +**------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + + int err; + + // Check if project includes leakage + closeleakage(pr); + hyd->HasLeakage = check_for_leakage(pr); + if (!hyd->HasLeakage) return 0; + + // Allocate memory for leakage data objects + err = create_leakage_objects(pr); + if (err > 0) return err; + + // Convert pipe leakage coeffs. to node coeffs. + convert_pipe_to_node_leakage(pr); + init_node_leakage(pr); + return 0; +} + + +int check_for_leakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: returns TRUE if any pipes can leak, FALSE otherwise +** Purpose: checks if any pipes can leak. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + int i; + Slink *link; + + for (i = 1; i <= net->Nlinks; i++) + { + // Only pipes have leakage + link = &net->Link[i]; + if (link->Type > PIPE) continue; + if (link->LeakArea > 0.0 || link->LeakExpan > 0.0) return TRUE; + } + return FALSE; +} + + +int create_leakage_objects(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: returns an error code +** Purpose: allocates an array of Leakage objects. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + int i; + + hyd->Leakage = (Sleakage *)calloc(net->Njuncs + 1, sizeof(Sleakage)); + if (hyd->Leakage == NULL) return 101; + for (i = 1; i <= net->Njuncs; i++) + { + hyd->Leakage[i].cfa = 0.0; + hyd->Leakage[i].cva = 0.0; + hyd->Leakage[i].qfa = 0.0; + hyd->Leakage[i].qva = 0.0; + } + return 0; +} + +void convert_pipe_to_node_leakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: converts pipe leakage parameters into node leakage +** coefficents. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + int i; + double c_area, c_expan, c_orif, len; + Slink *link; + Snode *node1; + Snode *node2; + + // Examine each link + c_orif = 4.8149866 * 1.e-6; + for (i = 1; i <= net->Nlinks; i++) + { + // Only pipes have leakage + link = &net->Link[i]; + if (link->Type > PIPE) continue; + + // Ignore leakage in a pipe connecting two tanks or + // reservoirs (since those nodes don't have demands) + node1 = &net->Node[link->N1]; + node2 = &net->Node[link->N2]; + if (node1->Type != JUNCTION && node2->Type != JUNCTION) continue; + + // Get pipe's fixed and variable area leak coeffs. + if (link->LeakArea == 0.0 && link->LeakExpan == 0.0) continue; + c_area = c_orif * link->LeakArea / SQR(MperFT); + c_expan = c_orif * link->LeakExpan; + + // Adjust for number of 100-ft pipe sections + len = link->Len * pr->Ucf[LENGTH] / 100.; + if (node1->Type == JUNCTION && node2->Type == JUNCTION) + { + len *= 0.5; + } + c_area *= len; + c_expan *= len; + + // Add these coeffs. to pipe's end nodes + if (node1->Type == JUNCTION) + { + hyd->Leakage[link->N1].cfa += c_area; + hyd->Leakage[link->N1].cva += c_expan; + } + if (node2->Type == JUNCTION) + { + hyd->Leakage[link->N2].cfa += c_area; + hyd->Leakage[link->N2].cva += c_expan; + } + } +} + +void init_node_leakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: initializes node leakage coeffs. and flows. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + int i; + double c_area, c_expan; + + for (i = 1; i <= net->Njuncs; i++) + { + // Coeff. for fixed area leakage + c_area = hyd->Leakage[i].cfa; + if (c_area > 0.0) + hyd->Leakage[i].cfa = 1.0 / (c_area * c_area); + else + hyd->Leakage[i].cfa = 0.0; + + // Coeff. for variable area leakage + c_expan = hyd->Leakage[i].cva; + if (c_expan > 0.0) + hyd->Leakage[i].cva = 1.0 / pow(c_expan, 2./3.); + else + hyd->Leakage[i].cva = 0.0; + + // Initialize leakage flow to a non-zero value (as required by + // the hydraulic solver) + if (hyd->Leakage[i].cfa > 0.0) + hyd->Leakage[i].qfa = 0.001; + if (hyd->Leakage[i].cva > 0.0) + hyd->Leakage[i].qva = 0.001; + } +} + +void closeleakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: frees memory for nodal leakage objects. +**------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + if (hyd->Leakage) free(hyd->Leakage); + hyd->Leakage = NULL; + hyd->HasLeakage = FALSE; +} + +double findlinkleakage(Project *pr, int i) +/*------------------------------------------------------------- +** Input: i = link index +** Output: returns link leakage flow (cfs) +** Purpose: computes leakage flow from link i at current +** hydraulic solution. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + Smatrix *sm = &hyd->smatrix; + Slink *link = &net->Link[i]; + + int n1, n2; + double h1, h2, hsqrt, a, m, c, len, q1, q2; + + // Only pipes can leak + link = &net->Link[i]; + if (link->Type > PIPE) return 0.0; + + // No leakage if area & expansion are 0 + if (link->LeakArea == 0.0 && link->LeakExpan == 0.0) return 0.0; + + // No leakage if link's end nodes are both fixed grade + n1 = link->N1; + n2 = link->N2; + if (n1 > net->Njuncs && n2 > net->Njuncs) return 0.0; + + // Pressure head of end nodes + h1 = hyd->NodeHead[n1] - net->Node[n1].El; + h1 = MAX(h1, 0.0); + h2 = hyd->NodeHead[n2] - net->Node[n2].El; + h2 = MAX(h2, 0.0); + + // Pipe leak parameters converted to feet + a = link->LeakArea / SQR(MperFT); + m = link->LeakExpan; + len = link->Len * pr->Ucf[LENGTH] / 100.; // # 100 ft pipe lengths + c = 4.8149866 * len / 2.0 * 1.e-6; + + // Leakage from 1st half of pipe connected to node n1 + q1 = 0.0; + if (n1 <= net->Njuncs) + { + hsqrt = sqrt(h1); + q1 = c * (a + m * h1) * hsqrt; + } + + // Leakage from 2nd half of pipe connected to node n2 + q2 = 0.0; + if (n2 <= net->Njuncs) + { + hsqrt = sqrt(h2); + q2 = c * (a + m * h2) * hsqrt; + } + + // Adjust leakage flows to account for one node being fixed grade + if (q2 == 0.0) q1 *= 2.0; + if (q1 == 0.0) q2 *= 2.0; + return q1 + q2; +} + +void leakagecoeffs(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes coeffs. of the linearized hydraulic eqns. +** contributed by node leakages. +**-------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + Smatrix *sm = &hyd->smatrix; + + int i, row; + double hfa, // head loss producing current fixed area leakage + gfa, // gradient of fixed area head loss + hva, // head loss producing current variable area leakage + gva; // gradient of variable area head loss + + Snode* node; + + for (i = 1; i <= net->Njuncs; i++) + { + // Skip junctions that don't leak + node = &net->Node[i]; + if (!leakage_headloss(pr, i, &hfa, &gfa, &hva, &gva)) continue; + + // Addition to matrix diagonal & r.h.s + row = sm->Row[i]; + if (gfa > 0.0) + { + sm->Aii[row] += 1.0 / gfa; + sm->F[row] += (hfa + node->El) / gfa; + } + if (gva > 0.0) + { + sm->Aii[row] += 1.0 / gva; + sm->F[row] += (hva + node->El) / gva; + } + + // Update node's flow excess (inflow - outflow) + hyd->Xflow[i] -= (hyd->Leakage[i].qfa + hyd->Leakage[i].qva); + } +} + +double leakageflowchange(Project *pr, int i) +/* +**-------------------------------------------------------------- +** Input: i = node index +** Output: returns change in leakage flow rate +** Purpose: finds new leakage flow rate at a node after new +** heads are computed by the hydraulic solver. +**-------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + double hfa, gfa, hva, gva, // same as defined in leakage_solvercoeffs() + dh, dqfa, dqva; + + // 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; + + // GGA flow update formula for fixed area leakage + dqfa = 0.0; + if (gfa > 0.0) + { + dqfa = (hfa - dh) / 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; + hyd->Leakage[i].qva -= dqva; + } + + // New leakage flow at the node + hyd->LeakageFlow[i] = hyd->Leakage[i].qfa + hyd->Leakage[i].qva; + return dqfa + dqva; +} + +int leakagehasconverged(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns TRUE if leakage calculations converged, +** FALSE if not +** Purpose: checks if leakage calculations have converged. +**-------------------------------------------------------------- +*/ +{ + 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; + + for (i = 1; i <= net->Njuncs; i++) + { + // Skip junctions that don't leak + if (hyd->Leakage[i].cfa == 0 && hyd->Leakage[i].cva == 0) continue; + + // Evaluate node's pressure head + h = hyd->NodeHead[i] - net->Node[i].El; + + // Directly compute a reference leakage at this pressure head + qref = 0.0; + // Contribution from pipes with fixed area leaks + if (hyd->Leakage[i].cfa > 0.0) + qref = sqrt(h / hyd->Leakage[i].cfa); + // Contribution from pipes with variable area leaks + if (hyd->Leakage[i].cva > 0.0) + qref += pow((h / hyd->Leakage[i].cva), 1.5); + + // Compare reference leakage to solution leakage + qtest = hyd->Leakage[i].qfa + hyd->Leakage[i].qva; + if (fabs(qref - qtest) > ABSTOL + RELTOL * qref) return FALSE; + } + return TRUE; +} + +int leakage_headloss(Project* pr, int i, double *hfa, double *gfa, + double *hva, double *gva) +/* +**-------------------------------------------------------------- +** Input: i = node index +** Output: hfa = fixed area leak head loss (ft) +** gfa = gradient of fixed area head loss (ft/cfs) +** hva = variable area leak head loss (ft) +** gva = gradient of variable area head loss (ft/cfs) +** returns TRUE if node has leakage, FALSE otherwise +** Purpose: finds head loss and its gradient for a node's +** leakage as a function of leakage flow. +**-------------------------------------------------------------- +*/ +{ + 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); + 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); + return TRUE; +} + +void eval_node_leakage(double RQtol, double q, double c, double n, + double *h, double *g) +/* +**-------------------------------------------------------------- +** 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) +** Purpose: evaluates inverted form of leakage equation to +** compute head loss and its gradient as a function +** flow. +**-------------------------------------------------------------- +*/ +{ + n = 1.0 / n; + *g = n * c * pow(fabs(q), n - 1.0); + + // Use linear head loss function for small gradient + if (*g < RQtol) + { + *g = RQtol / n; + *h = (*g) * q; + } + + // Otherwise use normal leakage head loss function + else *h = (*g) * q / n; + + // Prevent leakage from going negative + add_lower_barrier(q, h, g); +} + +void add_lower_barrier(double q, double* h, double* g) +/* +**-------------------------------------------------------------------- +** Input: q = current flow rate +** Output: h = head loss value +** g = 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); +} diff --git a/src/project.c b/src/project.c index 2104d060..80c76979 100644 --- a/src/project.c +++ b/src/project.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 09/28/2023 + Last Updated: 06/24/2024 ****************************************************************************** */ @@ -328,8 +328,11 @@ void initpointers(Project *pr) pr->hydraul.P = NULL; pr->hydraul.Y = NULL; pr->hydraul.Xflow = NULL; + pr->hydraul.FullDemand = NULL; pr->hydraul.DemandFlow = NULL; pr->hydraul.EmitterFlow = NULL; + pr->hydraul.LeakageFlow = NULL; + pr->hydraul.Leakage = NULL; pr->quality.NodeQual = NULL; pr->quality.PipeRateCoeff = NULL; @@ -392,14 +395,18 @@ int allocdata(Project *pr) pr->hydraul.NodeDemand = (double *)calloc(n, sizeof(double)); pr->hydraul.NodeHead = (double *)calloc(n, sizeof(double)); pr->quality.NodeQual = (double *)calloc(n, sizeof(double)); + pr->hydraul.FullDemand = (double *)calloc(n, sizeof(double)); pr->hydraul.DemandFlow = (double *)calloc(n, sizeof(double)); pr->hydraul.EmitterFlow = (double *)calloc(n, sizeof(double)); + pr->hydraul.LeakageFlow = (double *)calloc(n, sizeof(double)); ERRCODE(MEMCHECK(pr->network.Node)); ERRCODE(MEMCHECK(pr->hydraul.NodeDemand)); ERRCODE(MEMCHECK(pr->hydraul.NodeHead)); ERRCODE(MEMCHECK(pr->quality.NodeQual)); + ERRCODE(MEMCHECK(pr->hydraul.FullDemand)); ERRCODE(MEMCHECK(pr->hydraul.DemandFlow)); ERRCODE(MEMCHECK(pr->hydraul.EmitterFlow)); + ERRCODE(MEMCHECK(pr->hydraul.LeakageFlow)); } // Allocate memory for network links @@ -471,8 +478,10 @@ void freedata(Project *pr) free(pr->hydraul.LinkFlow); free(pr->hydraul.LinkSetting); free(pr->hydraul.LinkStatus); + free(pr->hydraul.FullDemand); free(pr->hydraul.DemandFlow); free(pr->hydraul.EmitterFlow); + free(pr->hydraul.LeakageFlow); free(pr->quality.NodeQual); // Free memory used for nodal adjacency lists diff --git a/src/report.c b/src/report.c index ffdc2d22..25bb1376 100644 --- a/src/report.c +++ b/src/report.c @@ -422,6 +422,44 @@ void writehydstat(Project *pr, int iter, double relerr) writeline(pr, " "); } +void writeflowbalance(Project *pr) +/* +**------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: writes hydraulic flow balance ratio to report file. +**------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + Report *rpt = &pr->report; + char s1[MAXMSG+1]; + double ucf = pr->Ucf[FLOW]; + + snprintf(s1, MAXMSG, "Hydraulic Flow Balance (%s)", rpt->Field[DEMAND].Units); + writeline(pr, s1); + snprintf(s1, MAXMSG, "================================"); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Total Inflow: %12.3f", hyd->FlowBalance.totalInflow*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Consumer Demand: %12.3f", hyd->FlowBalance.consumerDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Demand Deficit: %12.3f", hyd->FlowBalance.deficitDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Emitter Flow: %12.3f", hyd->FlowBalance.emitterDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Leakage Flow: %12.3f", hyd->FlowBalance.leakageDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Total Outflow: %12.3f", hyd->FlowBalance.totalOutflow*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Storage Flow: %12.3f", hyd->FlowBalance.storageDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Flow Ratio: %12.3f", hyd->FlowBalance.ratio); + writeline(pr, s1); + snprintf(s1, MAXMSG, "================================\n"); + writeline(pr, s1); +} + void writemassbalance(Project *pr) /* **------------------------------------------------------------- diff --git a/src/text.h b/src/text.h index 1a78573b..84b05db7 100755 --- a/src/text.h +++ b/src/text.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 05/11/2024 + Last Updated: 06/15/2024 ****************************************************************************** */ @@ -151,6 +151,9 @@ #define w_REQUIRED "REQ" #define w_EXPONENT "EXP" +#define w_AREA "AREA" +#define w_EXPAN "EXPAN" + #define w_SECONDS "SEC" #define w_MINUTES "MIN" #define w_HOURS "HOU" @@ -208,6 +211,7 @@ #define s_DEMANDS "[DEMANDS]" #define s_SOURCES "[SOURCES]" #define s_EMITTERS "[EMITTERS]" +#define s_LEAKAGE "[LEAKAGE]" #define s_PATTERNS "[PATTERNS]" #define s_CURVES "[CURVES]" #define s_QUALITY "[QUALITY]" diff --git a/src/types.h b/src/types.h index ff7c3599..30fd64e8 100755 --- a/src/types.h +++ b/src/types.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 07/17/2023 + Last Updated: 06/15/2024 ****************************************************************************** */ @@ -294,7 +294,7 @@ typedef enum { _VALVES, _CONTROLS, _RULES, _DEMANDS, _SOURCES, _EMITTERS, _PATTERNS, _CURVES, _QUALITY, _STATUS, _ROUGHNESS, _ENERGY, _REACTIONS, _MIXING, _REPORT, _TIMES, _OPTIONS, - _COORDS, _VERTICES, _LABELS, _BACKDROP, _TAGS, _END + _COORDS, _VERTICES, _LABELS, _BACKDROP, _TAGS, _LEAKAGE, _END } SectionType; typedef enum { @@ -413,6 +413,8 @@ typedef struct // Link Object double Kw; // wall react. coef. double R; // flow resistance double Rc; // reaction coeff. + double LeakArea; // leak area (sq mm per 100 pipe length units + double LeakExpan; // leak expansion (sq mm per unit of head) LinkType Type; // link type StatusType Status; // initial status Pvertices Vertices; // internal vertex coordinates @@ -549,6 +551,26 @@ typedef struct // Mass Balance Components double ratio; // ratio of mass added to mass lost } SmassBalance; +typedef struct +{ + double totalInflow; + double totalOutflow; + double consumerDemand; + double emitterDemand; + double leakageDemand; + double deficitDemand; + double storageDemand; + double ratio; +} SflowBalance; + +typedef struct // Node Leakage Object +{ + double qfa; // fixed area leakage flow + double qva; // variable area leakage flow + double cfa; // fixed area leakage coeff. + double cva; // variable area leakage coeff. +} Sleakage; + /* ------------------------------------------------------ Wrapper Data Structures @@ -712,9 +734,11 @@ typedef struct { double *NodeHead, // Node hydraulic heads - *NodeDemand, // Node demand + emitter flows - *DemandFlow, // Work array of demand flows - *EmitterFlow, // Emitter outflows + *NodeDemand, // Node total demand (consumer + emitter + leakage) + *FullDemand, // Required consumer demand + *DemandFlow, // Demand flow from nodes + *EmitterFlow, // Emitter flow from nodes + *LeakageFlow, // Leakage flow from nodes *LinkFlow, // Link flows *LinkSetting, // Link settings Htol, // Hydraulic head tolerance @@ -741,6 +765,7 @@ typedef struct { MaxHeadError, // Max. error for link head loss MaxFlowChange, // Max. change in link flow DemandReduction, // % demand reduction at pressure deficient nodes + LeakageLoss, // % system leakage loss RelaxFactor, // Relaxation factor for flow updating *P, // Inverse of head loss derivatives *Y, // Flow correction factors @@ -759,12 +784,18 @@ typedef struct { MaxCheck, // Hydraulic trials limit on status checks OpenHflag, // Hydraulic system opened flag Haltflag, // Flag to halt simulation - DeficientNodes; // Number of pressure deficient nodes + DeficientNodes, // Number of pressure deficient nodes + HasLeakage; // TRUE if project has non-zero leakage parameters + + Sleakage *Leakage; // Array of node leakage parameters StatusType *LinkStatus, // Link status *OldStatus; // Previous link/tank status + SflowBalance + FlowBalance; // Flow balance components + Smatrix smatrix; // Sparse matrix storage } Hydraul; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d90167f7..744a58aa 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -41,6 +41,7 @@ set(toolkit_test_srcs test_pda.cpp test_valve.cpp test_units.cpp + test_leakage.cpp ) add_executable(test_toolkit ${toolkit_test_srcs}) diff --git a/tests/test_leakage.cpp b/tests/test_leakage.cpp new file mode 100644 index 00000000..4a1c7161 --- /dev/null +++ b/tests/test_leakage.cpp @@ -0,0 +1,93 @@ +/* + ****************************************************************************** + Project: OWA EPANET + Version: 2.3 + Module: test_leakage.cpp + Description: Tests EPANET toolkit api functions + Authors: see AUTHORS + Copyright: see AUTHORS + License: see LICENSE + Last Updated: 06/26/2024 + ****************************************************************************** +*/ + +/* + Tests Pipe Leakage Feature +*/ + +#include + +#include "test_toolkit.hpp" + +BOOST_AUTO_TEST_SUITE (test_leakage) + +BOOST_AUTO_TEST_CASE(test_leakage_model) + +//#include +//#include +//#include "epanet2_2.h" + +//int main() +{ + int error = 0; + int Pipe21, Junc21, Junc22; + double pipe21Leak, junc21Leak, junc22Leak; + EN_Project ph = NULL; + + error = EN_createproject(&ph); + BOOST_REQUIRE(error == 0); + error = EN_open(ph, DATA_PATH_NET1, DATA_PATH_RPT, ""); +// error = EN_open(ph, "Net1.inp", "Net1.rpt", ""); + BOOST_REQUIRE(error == 0); + + // single period analysis + error = EN_settimeparam(ph, EN_DURATION, 0); + BOOST_REQUIRE(error == 0); + + // Get index of Pipe 21 + error = EN_getlinkindex(ph, "21", &Pipe21); + BOOST_REQUIRE(error == 0); + + // Set Pipe21 leak area to 1.0 sq mm per 100 ft of pipe + // and its expansion rate to 0.1 sq mm per ft of head + error = EN_setlinkvalue(ph, Pipe21, EN_LEAK_AREA, 1.0); + BOOST_REQUIRE(error == 0); + error = EN_setlinkvalue(ph, Pipe21, EN_LEAK_EXPAN, 0.1); + BOOST_REQUIRE(error == 0); + + // Solve for hydraulics + error = EN_solveH(ph); + BOOST_REQUIRE(error == 0); + + // Compute Pipe 21 leakage flow using the FAVAD formula + // Note: we can't just sum the leak rates at both end nodes + // together since in general the nodes can have leakage + // contributed by other connecting pipes. + error = EN_getlinkvalue(ph, Pipe21, EN_LINK_LEAKAGE, &pipe21Leak); + BOOST_REQUIRE(error == 0); +// printf("\n Pipe leakage flow: %.4f", pipe21Leak); + + // Retrieve leakage flow at end nodes + // Note: In this case all of the leakage at these nodes is from Pipe 21. + error = EN_getnodeindex(ph, "21", &Junc21); + BOOST_REQUIRE(error == 0); + error = EN_getnodeindex(ph, "22", &Junc22); + BOOST_REQUIRE(error == 0); + error = EN_getnodevalue(ph, Junc21, EN_LEAKAGEFLOW, &junc21Leak); + BOOST_REQUIRE(error == 0); + error = EN_getnodevalue(ph, Junc22, EN_LEAKAGEFLOW, &junc22Leak); + BOOST_REQUIRE(error == 0); + + // Check that the sum of the node leakages equals the pipe leakage + //printf("\n Node leakage flow: %.4f\n", junc21Leak + junc22Leak); + BOOST_REQUIRE(abs(pipe21Leak - (junc21Leak+junc22Leak)) < 0.01); + +// Clean up + error = EN_close(ph); + BOOST_REQUIRE(error == 0); + error = EN_deleteproject(ph); + BOOST_REQUIRE(error == 0); +// return 0; +} + +BOOST_AUTO_TEST_SUITE_END() From 037ca41af6f2c92d702fecf4a5fedcbc43b3354f Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Wed, 26 Jun 2024 11:34:19 -0400 Subject: [PATCH 2/5] Adding Pipe Leakage Modeling --- ReleaseNotes2_3.md | 17 +- include/epanet2.bas | 11 +- include/epanet2.cs | 25 +- include/epanet2.def | 17 +- include/epanet2.pas | 306 ++++++++++++----------- include/epanet2.vb | 14 +- include/epanet2_enums.h | 13 +- src/enumstxt.h | 4 +- src/epanet.c | 66 ++++- src/flowbalance.c | 191 +++++++++++++++ src/funcs.h | 19 +- src/hydcoeffs.c | 5 +- src/hydraul.c | 19 +- src/hydsolver.c | 93 ++++--- src/inpfile.c | 17 +- src/input1.c | 6 +- src/input2.c | 1 + src/input3.c | 41 ++++ src/leakage.c | 531 ++++++++++++++++++++++++++++++++++++++++ src/project.c | 11 +- src/report.c | 38 +++ src/text.h | 6 +- src/types.h | 43 +++- tests/CMakeLists.txt | 1 + tests/test_leakage.cpp | 93 +++++++ 25 files changed, 1366 insertions(+), 222 deletions(-) create mode 100644 src/flowbalance.c create mode 100644 src/leakage.c create mode 100644 tests/test_leakage.cpp diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md index ea5e96de..e561ee81 100644 --- a/ReleaseNotes2_3.md +++ b/ReleaseNotes2_3.md @@ -4,6 +4,7 @@ This document describes the changes and updates that have been made in version 2.3 of EPANET. - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files could be opened by the toolkit. + - A `EN_setcurvetype` function was added to allow API clients to set a curve's type (e.g., `EN_PUMP_CURVE,` `EN_VOLUME_CURVE,` etc.). - A `EN_setvertex` function was added to allow API clients to change the coordinates of a single link vertex. - The indices of a General Purpose Valve (GPV) and a Positional Control Valve (PCV) were added to the list of editable Link Properties using the symbolic constant names `EN_GPV_CURVE` and `EN_PCV_CURVE`, respectively. @@ -52,6 +53,20 @@ This document describes the changes and updates that have been made in version 2 - Setting the demand multiplier within the `[DEMANDS]` section of INP has been depreciated, please use `DEMAND MULTIPLIER` inside `[OPTIONS]` instead. - `EN_PRESS_UNITS` can now be used with `EN_getoption` and `EN_setoption` to get or set the pressure unit used in EPANET. - Continuous barrier functions were added to constrain emitter flows to allowable values. -- The `EN_openx` function has been added to enable the opening of input files with formatting errors through the API. This allows users to continue using toolkit functions even when such errors are present. +- The `EN_openX` function has been added to enable the opening of input files with formatting errors through the API. This allows users to continue using toolkit functions even when such errors are present. - The `EN_getnodesvalues` and `EN_getlinksvalues` were added to retrieve a property value for all nodes or links in the network. - Fixed a bug in EN_setnodevalue with EN_EMITTER option that could cause NaN results. +- Support has been added for FAVAD (Fixed And Variable Area Discharge) modeling of pipe leaks: + - 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 \ No newline at end of file diff --git a/include/epanet2.bas b/include/epanet2.bas index 6d707591..6ddb0ada 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -5,7 +5,7 @@ Attribute VB_Name = "Module1" 'Declarations of functions in the EPANET PROGRAMMERs TOOLKIT '(EPANET2.DLL) -'Last updated on 09/28/2023 +'Last updated on 06/23/2024 ' These are codes used by the DLL functions Public Const EN_ELEVATION = 0 ' Node parameters @@ -38,6 +38,9 @@ Public Const EN_CANOVERFLOW = 26 Public Const EN_DEMANDDEFICIT = 27 Public Const EN_NODE_INCONTROL = 28 Public Const EN_EMITTERFLOW = 29 +Public Const EN_LEAKAGEFLOW = 30 +Public Const EN_DEMANDFLOW = 31 +Public Const EN_FULLDEMAND = 32 Public Const EN_DIAMETER = 0 ' Link parameters Public Const EN_LENGTH = 1 @@ -65,6 +68,9 @@ Public Const EN_PUMP_EPAT = 22 Public Const EN_LINK_INCONTROL = 23 Public Const EN_GPV_CURVE = 24 Public Const EN_PCV_CURVE = 25 +Public Const EN_LEAK_AREA = 26 +Public Const EN_LEAK_EXPAN = 27 +Public Const EN_LINK_LEAKAGE = 28 Public Const EN_DURATION = 0 ' Time parameters Public Const EN_HYDSTEP = 1 @@ -90,6 +96,7 @@ Public Const EN_MAXFLOWCHANGE = 3 Public Const EN_MASSBALANCE = 4 Public Const EN_DEFICIENTNODES = 5 Public Const EN_DEMANDREDUCTION = 6 +Public Const EN_LEAKAGELOSS = 7 Public Const EN_NODE = 0 ' Component types Public Const EN_LINK = 1 @@ -350,6 +357,7 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENsettankdata Lib "epanet2.dll" (ByVal index As Long, ByVal elev As Single, ByVal initlvl As Single, ByVal minlvl As Single, ByVal maxlvl As Single, ByVal diam As Single, ByVal minvol As Single, ByVal volcurve As String) As Long Declare Function ENgetcoord Lib "epanet2.dll" (ByVal index As Long, x As Double, y As Double) As Long Declare Function ENsetcoord Lib "epanet2.dll" (ByVal index As Long, ByVal x As Double, ByVal y As Double) As Long + Declare Function ENgetnodevalues Lib "epanet2.dll" (ByVal property as Long, values as Any) As Long 'Nodal Demand Functions Declare Function ENgetdemandmodel Lib "epanet2.dll" (type_ As Long, pmin As Single, preq As Single, pexp As Single) As Long @@ -382,6 +390,7 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENgetvertex Lib "epanet2.dll" (ByVal index As Long, ByVal vertex As Long, x As Double, y As Double) As Long Declare Function ENsetvertex Lib "epanet2.dll" (ByVal index As Long, ByVal vertex As Long, ByVal x As Double, ByVal y As Double) As Long Declare Function ENsetvertices Lib "epanet2.dll" (ByVal index As Long, xCoords As Any, yCoords As Any, ByVal count As Long) As Long + Declare Function ENgetlinkvalues Lib "epanet2.dll" (ByVal property as Long, values as Any) As Long 'Pump Functions Declare Function ENgetheadcurveindex Lib "epanet2.dll" (ByVal linkIndex As Long, curveIndex As Long) As Long diff --git a/include/epanet2.cs b/include/epanet2.cs index 57668bf7..0327c4d8 100644 --- a/include/epanet2.cs +++ b/include/epanet2.cs @@ -3,7 +3,7 @@ using System.Runtime.InteropServices; //epanet2.cs[By Oscar Vegas] -//Last updated on 09/28/2023 +//Last updated on 06/23/2024 //Declarations of functions in the EPANET PROGRAMMERs TOOLKIT //(EPANET2.DLL) for use with C# @@ -50,6 +50,9 @@ public static class Epanet public const int EN_DEMANDDEFICIT = 27; public const int EN_NODE_INCONTROL = 28; public const int EN_EMITTERFLOW = 29; + public const int EN_LEAKAGEFLOW = 30; + public const int EN_DEMANDFLOW = 31; + public const int EN_FULLDEMAND = 32; public const int EN_DIAMETER = 0; //Link parameters public const int EN_LENGTH = 1; @@ -78,6 +81,9 @@ public static class Epanet public const int EN_LINK_INCONTROL = 23; public const int EN_GPV_CURVE = 24; public const int EN_PCV_CURVE = 25; + public const int EN_LEAK_AREA = 26; + public const int EN_LEAK_EXPAN = 27; + public const int EN_LINK_LEAKAGE = 28; public const int EN_DURATION = 0; //Time parameters public const int EN_HYDSTEP = 1; @@ -102,6 +108,7 @@ public static class Epanet public const int EN_MASSBALANCE = 4; public const int EN_DEFICIENTNODES = 5; public const int EN_DEMANDREDUCTION = 6; + public const int EN_LEAKAGELOSS = 7; public const int EN_NODE = 0; //Component types public const int EN_LINK = 1; @@ -390,6 +397,9 @@ public static class Epanet [DllImport(EPANETDLL, EntryPoint = "ENgetresultindex")] public static extern int ENgetresultindex(int type, int index, ref int value); + [DllImport(EPANETDLL, EntryPoint = "ENtimetonextevent")] + public static extern int ENtimetonextevent(ref int eventType, ref long duration, ref int elementIndex); + //Analysis Options Functions [DllImport(EPANETDLL, EntryPoint = "ENgetoption")] @@ -440,10 +450,10 @@ public static class Epanet public static extern int ENgetnodetype(int index, ref int nodeType); [DllImport(EPANETDLL, EntryPoint = "ENgetnodevalue")] - public static extern int ENgetnodevalue(int index, int paramcode, ref float value); + public static extern int ENgetnodevalue(int index, int param, ref float value); [DllImport(EPANETDLL, EntryPoint = "ENsetnodevalue")] - public static extern int ENsetnodevalue(int index, int code, float value); + public static extern int ENsetnodevalue(int index, int param, float value); [DllImport(EPANETDLL, EntryPoint = "ENsetjuncdata")] public static extern int ENsetjuncdata(int index, float elev, float dmnd, string dmndpat); @@ -457,6 +467,8 @@ public static class Epanet [DllImport(EPANETDLL, EntryPoint = "ENsetcoord")] public static extern int ENsetcoord(int index, double x, double y); + [DllImport(EPANETDLL, EntryPoint = "ENgetnodevalues")] + public static extern int ENgetnodevalues(int param, ref float values); //Nodal Demand Functions [DllImport(EPANETDLL, EntryPoint = "ENgetdemandmodel")] @@ -525,10 +537,10 @@ public static class Epanet public static extern int ENsetlinknodes(int index, int node1, int node2); [DllImport(EPANETDLL, EntryPoint = "ENgetlinkvalue")] - public static extern int ENgetlinkvalue(int index, int code, ref float value); + public static extern int ENgetlinkvalue(int index, int param, ref float value); [DllImport(EPANETDLL, EntryPoint = "ENsetlinkvalue")] - public static extern int ENsetlinkvalue(int index, int code, float value); + public static extern int ENsetlinkvalue(int index, int param, float value); [DllImport(EPANETDLL, EntryPoint = "ENsetpipedata")] public static extern int ENsetpipedata(int index, float length, float diam, float rough, float mloss); @@ -542,6 +554,9 @@ public static class Epanet [DllImport(EPANETDLL, EntryPoint = "ENsetvertices")] public static extern int ENsetvertices(int index, ref double[] x, ref double[] y, int count); + [DllImport(EPANETDLL, EntryPoint = "ENgetlinkvalues")] + public static extern int ENgetlinkvalues(int param, ref float values); + //Pump Functions [DllImport(EPANETDLL, EntryPoint = "ENgetheadcurveindex")] diff --git a/include/epanet2.def b/include/epanet2.def index f5372af4..d3472c10 100644 --- a/include/epanet2.def +++ b/include/epanet2.def @@ -25,6 +25,7 @@ EXPORTS ENgetbasedemand = _ENgetbasedemand@12 ENgetcomment = _ENgetcomment@12 ENgetcontrol = _ENgetcontrol@24 + ENgetcontrolenabled = _ENgetcontrolenabled@8 ENgetcoord = _ENgetcoord@12 ENgetcount = _ENgetcount@8 ENgetcurve = _ENgetcurve@20 @@ -44,13 +45,14 @@ EXPORTS ENgetlinkid = _ENgetlinkid@8 ENgetlinkindex = _ENgetlinkindex@8 ENgetlinknodes = _ENgetlinknodes@12 - ENsetlinknodes = _ENsetlinknodes@12 ENgetlinktype = _ENgetlinktype@8 ENgetlinkvalue = _ENgetlinkvalue@12 + ENgetlinkvalues = _ENgetlinkvalues@8 ENgetnodeid = _ENgetnodeid@8 ENgetnodeindex = _ENgetnodeindex@8 ENgetnodetype = _ENgetnodetype@8 - ENgetnodevalue = _ENgetnodevalue@12 + ENgetnodevalue = _ENgetnodevalue@12 + ENgetnodevalues = _ENgetnodevalues@8 ENgetnumdemands = _ENgetnumdemands@8 ENgetoption = _ENgetoption@8 ENgetpatternid = _ENgetpatternid@8 @@ -63,6 +65,7 @@ EXPORTS ENgetqualtype = _ENgetqualtype@8 ENgetresultindex = _ENgetresultindex@12 ENgetrule = _ENgetrule@20 + ENgetruleenabled = _ENgetruleenabled@8 ENgetruleID = _ENgetruleID@8 ENgetstatistic = _ENgetstatistic@8 ENgetthenaction = _ENgetthenaction@20 @@ -79,6 +82,7 @@ EXPORTS ENopen = _ENopen@12 ENopenH = _ENopenH@0 ENopenQ = _ENopenQ@0 + ENopenX = _ENopenX@12 ENreport = _ENreport@0 ENresetreport = _ENresetreport@0 ENrunH = _ENrunH@4 @@ -89,6 +93,7 @@ EXPORTS ENsetbasedemand = _ENsetbasedemand@12 ENsetcomment = _ENsetcomment@12 ENsetcontrol = _ENsetcontrol@24 + ENsetcontrolenabled = _ENsetcontrolenabled@8 ENsetcoord = _ENsetcoord@20 ENsetcurve = _ENsetcurve@16 ENsetcurveid = _ENsetcurveid@8 @@ -118,6 +123,7 @@ EXPORTS ENsetpremisevalue = _ENsetpremisevalue@12 ENsetqualtype = _ENsetqualtype@16 ENsetreport = _ENsetreport@4 + ENsetruleenabled = _ENsetruleenabled@8 ENsetrulepriority = _ENsetrulepriority@8 ENsetstatusreport = _ENsetstatusreport@4 ENsettankdata = _ENsettankdata@32 @@ -129,11 +135,6 @@ EXPORTS ENsolveH = _ENsolveH@0 ENsolveQ = _ENsolveQ@0 ENstepQ = _ENstepQ@4 + ENtimetonextevent = _ENtimetonextevent@12 ENusehydfile = _ENusehydfile@4 ENwriteline = _ENwriteline@4 - ENtimetonextevent = _ENtimetonextevent@12 - ENopenX = _ENopenX@12 - ENgetcontrolenabled = _ENgetcontrolenabled@8 - ENsetcontrolenabled = _ENsetcontrolenabled@8 - ENgetruleenabled = _ENgetruleenabled@8 - ENsetruleenabled = _ENsetruleenabled@8 diff --git a/include/epanet2.pas b/include/epanet2.pas index 4ce4a03c..6f2e8dae 100644 --- a/include/epanet2.pas +++ b/include/epanet2.pas @@ -3,7 +3,7 @@ { Declarations of imported procedures from the EPANET PROGRAMMERs TOOLKIT } { (EPANET2.DLL) } -{Last updated on 09/28/2023} +{Last updated on 06/06/2024} interface @@ -46,6 +46,9 @@ interface EN_DEMANDDEFICIT = 27; EN_NODE_INCONTROL = 28; EN_EMITTERFLOW = 29; + EN_LEAKAGEFLOW = 30; + EN_DEMANDFLOW = 31; + EN_FULLDEMAND = 32; EN_DIAMETER = 0; { Link parameters } EN_LENGTH = 1; @@ -71,8 +74,11 @@ interface EN_PUMP_ECOST = 21; EN_PUMP_EPAT = 22; EN_LINK_INCONTROL = 23; - EN_GPV_CURVE = 24; - EN_PCV_CURVE = 25; + EN_GPV_CURVE = 24; + EN_PCV_CURVE = 25; + EN_LEAK_AREA = 26; + EN_LEAK_EXPAN = 27; + EN_LINK_LEAKAGE = 28; EN_DURATION = 0; { Time parameters } EN_HYDSTEP = 1; @@ -98,6 +104,7 @@ interface EN_MASSBALANCE = 4; EN_DEFICIENTNODES = 5; EN_DEMANDREDUCTION = 6; + EN_LEAKAGELOSS = 7; EN_NODE = 0; { Component Types } EN_LINK = 1; @@ -275,180 +282,199 @@ interface EN_TRUE = 1; { boolean true } +{$MACRO ON} + {$ifdef MSWINDOWS} - EpanetLib = 'epanet2.dll'; + EpanetLib = 'epanet2.dll'; + {$DEFINE cdecl := stdcall} +{$endif} + +{$ifdef UNIX} + {$ifdef DARWIN} + EpanetLib = 'libepanet2.dylib'; + {$linklib libepanet2} + {$else} + EpanetLib = 'libepanet2.so'; + {$endif} +{$endif} + +{$ifdef UNIX} + {$DEFINE TimeType:=Int64} {$else} - EpanetLib = 'libepanet2.so'; + {$DEFINE TimeType:=Integer} {$endif} {Project Functions} - function ENepanet(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar; F4: Pointer): Integer; stdcall; external EpanetLib; - function ENinit(F2: PAnsiChar; F3: PAnsiChar; UnitsType: Integer; HeadlossType: Integer): Integer; stdcall; external EpanetLib; - function ENopen(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENopenX(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetcount(Code: Integer; var Count: Integer): Integer; stdcall; external EpanetLib; - function ENgettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsaveinpfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENclose: Integer; stdcall; external EpanetLib; + function ENepanet(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar; F4: Pointer): Integer; cdecl; external EpanetLib; + function ENinit(F2: PAnsiChar; F3: PAnsiChar; UnitsType: Integer; HeadlossType: Integer): Integer; cdecl; external EpanetLib; + function ENopen(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENopenX(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetcount(Code: Integer; var Count: Integer): Integer; cdecl; external EpanetLib; + function ENgettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsaveinpfile(F: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENclose: Integer; cdecl; external EpanetLib; {Hydraulic Functions} - function ENsolveH: Integer; stdcall; external EpanetLib; - function ENsaveH: Integer; stdcall; external EpanetLib; - function ENopenH: Integer; stdcall; external EpanetLib; - function ENinitH(SaveFlag: Integer): Integer; stdcall; external EpanetLib; - function ENrunH(var T: Integer): Integer; stdcall; external EpanetLib; - function ENnextH(var Tstep: Integer): Integer; stdcall; external EpanetLib; - function ENcloseH: Integer; stdcall; external EpanetLib; - function ENsavehydfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENusehydfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsolveH: Integer; cdecl; external EpanetLib; + function ENsaveH: Integer; cdecl; external EpanetLib; + function ENopenH: Integer; cdecl; external EpanetLib; + function ENinitH(SaveFlag: Integer): Integer; cdecl; external EpanetLib; + function ENrunH(var T: TimeType): Integer; cdecl; external EpanetLib; + function ENnextH(var Tstep: TimeType): Integer; cdecl; external EpanetLib; + function ENcloseH: Integer; cdecl; external EpanetLib; + function ENsavehydfile(F: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENusehydfile(F: PAnsiChar): Integer; cdecl; external EpanetLib; {Quality Functions} - function ENsolveQ: Integer; stdcall; external EpanetLib; - function ENopenQ: Integer; stdcall; external EpanetLib; - function ENinitQ(SaveFlag: Integer): Integer; stdcall; external EpanetLib; - function ENrunQ(var T: Integer): Integer; stdcall; external EpanetLib; - function ENnextQ(var Tstep: Integer): Integer; stdcall; external EpanetLib; - function ENstepQ(var Tleft: Integer): Integer; stdcall; external EpanetLib; - function ENcloseQ: Integer; stdcall; external EpanetLib; + function ENsolveQ: Integer; cdecl; external EpanetLib; + function ENopenQ: Integer; cdecl; external EpanetLib; + function ENinitQ(SaveFlag: Integer): Integer; cdecl; external EpanetLib; + function ENrunQ(var T: TimeType): Integer; cdecl; external EpanetLib; + function ENnextQ(var Tstep: TimeType): Integer; cdecl; external EpanetLib; + function ENstepQ(var Tleft: TimeType): Integer; cdecl; external EpanetLib; + function ENcloseQ: Integer; cdecl; external EpanetLib; {Reporting Functions} - function ENwriteline(S: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENreport: Integer; stdcall; external EpanetLib; - function ENcopyreport(F: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENclearreport: Integer; stdcall; external EpanetLib; - function ENresetreport: Integer; stdcall; external EpanetLib; - function ENsetreport(S: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetstatusreport(Code: Integer): Integer; stdcall; external EpanetLib; - function ENgetversion(var Version: Integer): Integer; stdcall; external EpanetLib; - function ENgeterror(Errcode: Integer; Errmsg: PAnsiChar; MaxLen: Integer): Integer; stdcall; external EpanetLib; - function ENgetstatistic(StatType: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENgetresultindex(Code: Integer; Index: Integer; var Value: Integer): Integer; stdcall; external EpanetLib; + function ENwriteline(S: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENreport: Integer; cdecl; external EpanetLib; + function ENcopyreport(F: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENclearreport: Integer; cdecl; external EpanetLib; + function ENresetreport: Integer; cdecl; external EpanetLib; + function ENsetreport(S: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetstatusreport(Code: Integer): Integer; cdecl; external EpanetLib; + function ENgetversion(var Version: Integer): Integer; cdecl; external EpanetLib; + function ENgeterror(Errcode: Integer; Errmsg: PAnsiChar; MaxLen: Integer): Integer; cdecl; external EpanetLib; + function ENgetstatistic(StatType: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENgetresultindex(Code: Integer; Index: Integer; var Value: Integer): Integer; cdecl; external EpanetLib; + function ENtimetonextevent(var EventType: Integer; var Duration: TimeType; var ElementIndex: Integer): Integer; cdecl; external EpanetLib; {Analysis Options Functions} - function ENgetoption(Code: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetoption(Code: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENgetflowunits(var Code: Integer): Integer; stdcall; external EpanetLib; - function ENsetflowunits(Code: Integer): Integer; stdcall; external EpanetLib; - function ENgettimeparam(Code: Integer; var Value: Integer): Integer; stdcall; external EpanetLib; - function ENsettimeparam(Code: Integer; Value: Integer): Integer; stdcall; external EpanetLib; - function ENgetqualinfo(var QualType: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; var TraceNode: Integer): Integer; stdcall; external EpanetLib; - function ENgetqualtype(var QualCode: Integer; var TraceNode: Integer): Integer; stdcall; external EpanetLib; - function ENsetqualtype(QualCode: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; TraceNodeID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetoption(Code: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetoption(Code: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENgetflowunits(var Code: Integer): Integer; cdecl; external EpanetLib; + function ENsetflowunits(Code: Integer): Integer; cdecl; external EpanetLib; + function ENgettimeparam(Code: Integer; var Value: TimeType): Integer; cdecl; external EpanetLib; + function ENsettimeparam(Code: Integer; Value: TimeType): Integer; cdecl; external EpanetLib; + function ENgetqualinfo(var QualType: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; var TraceNode: Integer): Integer; cdecl; external EpanetLib; + function ENgetqualtype(var QualCode: Integer; var TraceNode: Integer): Integer; cdecl; external EpanetLib; + function ENsetqualtype(QualCode: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; TraceNodeID: PAnsiChar): Integer; cdecl; external EpanetLib; {Node Functions} - function ENaddnode(ID: PAnsiChar; NodeType: Integer; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENdeletenode(Index: Integer; ActionCode: Integer): Integer; stdcall; external EpanetLib; - function ENgetnodeindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetnodeid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetnodeid(Index: Integer; NewID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetnodetype(Index: Integer; var Code: Integer): Integer; stdcall; external EpanetLib; - function ENgetnodevalue(Index: Integer; Code: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetnodevalue(Index: Integer; Code: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENsetjuncdata(Index: Integer; Elev: Single; Dmnd: Single; DmndPat: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsettankdata(Index: Integer; Elev, InitLvl, MinLvl, MaxLvl, Diam, MinVol: Single; VolCurve: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetcoord(Index: Integer; var X: Double; var Y: Double): Integer; stdcall; external EpanetLib; - function ENsetcoord(Index: Integer; X: Double; Y: Double): Integer; stdcall; external EpanetLib; + function ENaddnode(ID: PAnsiChar; NodeType: Integer; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENdeletenode(Index: Integer; ActionCode: Integer): Integer; cdecl; external EpanetLib; + function ENgetnodeindex(ID: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetnodeid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetnodeid(Index: Integer; NewID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetnodetype(Index: Integer; var Code: Integer): Integer; cdecl; external EpanetLib; + function ENgetnodevalue(Index: Integer; Code: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetnodevalue(Index: Integer; Code: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENsetjuncdata(Index: Integer; Elev: Single; Dmnd: Single; DmndPat: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsettankdata(Index: Integer; Elev, InitLvl, MinLvl, MaxLvl, Diam, MinVol: Single; VolCurve: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetcoord(Index: Integer; var X: Double; var Y: Double): Integer; cdecl; external EpanetLib; + function ENsetcoord(Index: Integer; X: Double; Y: Double): Integer; cdecl; external EpanetLib; + function ENgetnodevalues(Code: Integer; var X: array of Single): Integer; cdecl; external EpanetLib; {Demand Functions} - function ENgetdemandmodel(var Model: Integer; var Pmin: Single; var Preq: Single; var Pexp: Single): Integer; stdcall; external EpanetLib; - function ENsetdemandmodel(Model: Integer; Pmin: Single; Preq: Single; Pexp: Single): Integer; stdcall; external EpanetLib; - function ENgetnumdemands(NodeIndex: Integer; var NumDemands: Integer): Integer; stdcall; external EpanetLib; - function ENadddemand(NodeIndex: Integer; BaseDemand: Single; PatName: PAnsiChar; DemandName: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENdeletedemand(NodeIndex: Integer; DemandIndex: Integer): Integer; stdcall; external EpanetLib; - function ENgetdemandindex(NodeIndex: Integer; DemandName: PAnsiChar; var DemandIndex: Integer): Integer; stdcall; external EpanetLib; - function ENgetbasedemand(NodeIndex: Integer; DemandIndex: Integer; var BaseDemand: Single): Integer; stdcall; external EpanetLib; - function ENsetbasedemand(NodeIndex: Integer; DemandIndex: Integer; BaseDemand: Single): Integer; stdcall; external EpanetLib; - function ENgetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; var PatIndex: Integer): Integer; stdcall; external EpanetLib; - function ENsetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; PatIndex: Integer): Integer; stdcall; external EpanetLib; - function ENgetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetdemandmodel(var Model: Integer; var Pmin: Single; var Preq: Single; var Pexp: Single): Integer; cdecl; external EpanetLib; + function ENsetdemandmodel(Model: Integer; Pmin: Single; Preq: Single; Pexp: Single): Integer; cdecl; external EpanetLib; + function ENgetnumdemands(NodeIndex: Integer; var NumDemands: Integer): Integer; cdecl; external EpanetLib; + function ENadddemand(NodeIndex: Integer; BaseDemand: Single; PatName: PAnsiChar; DemandName: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENdeletedemand(NodeIndex: Integer; DemandIndex: Integer): Integer; cdecl; external EpanetLib; + function ENgetdemandindex(NodeIndex: Integer; DemandName: PAnsiChar; var DemandIndex: Integer): Integer; cdecl; external EpanetLib; + function ENgetbasedemand(NodeIndex: Integer; DemandIndex: Integer; var BaseDemand: Single): Integer; cdecl; external EpanetLib; + function ENsetbasedemand(NodeIndex: Integer; DemandIndex: Integer; BaseDemand: Single): Integer; cdecl; external EpanetLib; + function ENgetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; var PatIndex: Integer): Integer; cdecl; external EpanetLib; + function ENsetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; PatIndex: Integer): Integer; cdecl; external EpanetLib; + function ENgetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; cdecl; external EpanetLib; {Link Functions} - function ENaddlink(ID: PAnsiChar; LinkType: Integer; FromNode: PAnsiChar; ToNode: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENdeletelink(Index: Integer; ActionCode: Integer): Integer; stdcall; external EpanetLib; - function ENgetlinkindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetlinkid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetlinkid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetlinktype(Index: Integer; var Code: Integer): Integer; stdcall; external EpanetLib; - function ENsetlinktype(var Index: Integer; LinkType: Integer; ActionCode: Integer): Integer; stdcall; external EpanetLib; - function ENgetlinknodes(Index: Integer; var Node1: Integer; var Node2: Integer): Integer; stdcall; external EpanetLib; - function ENsetlinknodes(Index: Integer; Node1: Integer; Node2: Integer): Integer; stdcall; external EpanetLib; - function ENgetlinkvalue(Index: Integer; Code: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetlinkvalue(Index: Integer; Code: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENsetpipedata(Index: Integer; Length: Single; Diam: Single; Rough: Single; Mloss:Single): Integer; stdcall; external EpanetLib; - - function ENgetvertexcount(Index: Integer; var Count: Integer): Integer; stdcall; external EpanetLib; - function ENgetvertex(Index: Integer; Vertex: Integer; var X: Double; var Y: Double): Integer; stdcall; external EpanetLib; - function ENsetvertex(Index: Integer; Vertex: Integer; X: Double; Y: Double): Integer; stdcall; external EpanetLib; - function ENsetvertices(Index: Integer; var X: Double; var Y: Double; Count: Integer): Integer; stdcall; external EpanetLib; + function ENaddlink(ID: PAnsiChar; LinkType: Integer; FromNode: PAnsiChar; ToNode: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENdeletelink(Index: Integer; ActionCode: Integer): Integer; cdecl; external EpanetLib; + function ENgetlinkindex(ID: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetlinkid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetlinkid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetlinktype(Index: Integer; var Code: Integer): Integer; cdecl; external EpanetLib; + function ENsetlinktype(var Index: Integer; LinkType: Integer; ActionCode: Integer): Integer; cdecl; external EpanetLib; + function ENgetlinknodes(Index: Integer; var Node1: Integer; var Node2: Integer): Integer; cdecl; external EpanetLib; + function ENsetlinknodes(Index: Integer; Node1: Integer; Node2: Integer): Integer; cdecl; external EpanetLib; + function ENgetlinkvalue(Index: Integer; Code: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetlinkvalue(Index: Integer; Code: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENsetpipedata(Index: Integer; Length: Single; Diam: Single; Rough: Single; Mloss:Single): Integer; cdecl; external EpanetLib; + function ENgetlinkvalues(Code: Integer; var X: array of Single): Integer; cdecl; external EpanetLib; + + function ENgetvertexcount(Index: Integer; var Count: Integer): Integer; cdecl; external EpanetLib; + function ENgetvertex(Index: Integer; Vertex: Integer; var X: Double; var Y: Double): Integer; cdecl; external EpanetLib; + function ENsetvertex(Index: Integer; Vertex: Integer; X: Double; Y: Double): Integer; cdecl; external EpanetLib; + function ENsetvertices(Index: Integer; var X: Double; var Y: Double; Count: Integer): Integer; cdecl; external EpanetLib; {Pump Functions} - function ENgetpumptype(LinkIndex: Integer; var PumpType: Integer): Integer; stdcall; external EpanetLib; - function ENgetheadcurveindex(LinkIndex: Integer; var CurveIndex: Integer): Integer; stdcall; external EpanetLib; - function ENsetheadcurveindex(LinkIndex: Integer; CurveIndex: Integer): Integer; stdcall; external EpanetLib; + function ENgetpumptype(LinkIndex: Integer; var PumpType: Integer): Integer; cdecl; external EpanetLib; + function ENgetheadcurveindex(LinkIndex: Integer; var CurveIndex: Integer): Integer; cdecl; external EpanetLib; + function ENsetheadcurveindex(LinkIndex: Integer; CurveIndex: Integer): Integer; cdecl; external EpanetLib; {Pattern Functions} - function ENaddpattern(ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENdeletepattern(Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetpatternindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetpatternid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetpatternid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetpatternlen(Index: Integer; var Len: Integer): Integer; stdcall; external EpanetLib; - function ENgetpatternvalue(Index: Integer; Period: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetpatternvalue(Index: Integer; Period: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENgetaveragepatternvalue(Index: Integer; var Value: Single): Integer; stdcall; external EpanetLib; - function ENsetpattern(Index: Integer; var F: Single; N: Integer): Integer; stdcall; external EpanetLib; + function ENaddpattern(ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENdeletepattern(Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetpatternindex(ID: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetpatternid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetpatternid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetpatternlen(Index: Integer; var Len: Integer): Integer; cdecl; external EpanetLib; + function ENgetpatternvalue(Index: Integer; Period: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetpatternvalue(Index: Integer; Period: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENgetaveragepatternvalue(Index: Integer; var Value: Single): Integer; cdecl; external EpanetLib; + function ENsetpattern(Index: Integer; var F: Single; N: Integer): Integer; cdecl; external EpanetLib; {Curve Functions} - function ENaddcurve(ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENdeletecurve(Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetcurveindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetcurveid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetcurveid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENgetcurvelen(Index: Integer; var Len: Integer): Integer; stdcall; external EpanetLib; - function ENgetcurvetype(Index: Integer; var CurveType: Integer): Integer; stdcall; external EpanetLib; - function ENsetcurvetype(Index: Integer; CurveType: Integer): Integer; stdcall; external EpanetLib; - function ENgetcurvevalue(CurveIndex: Integer; PointIndex: Integer; var X: Single; var Y: Single): Integer; stdcall; external EpanetLib; - function ENsetcurvevalue(CurveIndex: Integer; PointIndex: Integer; X: Single; Y: Single): Integer; stdcall; external EpanetLib; - function ENgetcurve(Index: Integer; ID: PAnsiChar; var N: Integer; var X: Single; var Y: Single): Integer; stdcall; external EpanetLib; - function ENsetcurve(Index: Integer; var X: Single; var Y: Single; N: Integer): Integer; stdcall; external EpanetLib; + function ENaddcurve(ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENdeletecurve(Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetcurveindex(ID: PAnsiChar; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetcurveid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetcurveid(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENgetcurvelen(Index: Integer; var Len: Integer): Integer; cdecl; external EpanetLib; + function ENgetcurvetype(Index: Integer; var CurveType: Integer): Integer; cdecl; external EpanetLib; + function ENsetcurvetype(Index: Integer; CurveType: Integer): Integer; cdecl; external EpanetLib; + function ENgetcurvevalue(CurveIndex: Integer; PointIndex: Integer; var X: Single; var Y: Single): Integer; cdecl; external EpanetLib; + function ENsetcurvevalue(CurveIndex: Integer; PointIndex: Integer; X: Single; Y: Single): Integer; cdecl; external EpanetLib; + function ENgetcurve(Index: Integer; ID: PAnsiChar; var N: Integer; var X: Single; var Y: Single): Integer; cdecl; external EpanetLib; + function ENsetcurve(Index: Integer; var X: Single; var Y: Single; N: Integer): Integer; cdecl; external EpanetLib; {Control Functions} - function ENaddcontrol(Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single; var Index: Integer): Integer; stdcall; external EpanetLib; - function ENdeletecontrol(Index: Integer): Integer; stdcall; external EpanetLib; - function ENgetcontrol(Index: Integer; var Ctype: Integer; var Link: Integer; var Setting: Single; var Node: Integer; var Level: Single): Integer; stdcall; external EpanetLib; - function ENsetcontrol(Index: Integer; Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single): Integer; stdcall; external EpanetLib; - function ENgetcontrolenabled(Index: Integer; out_enabled: Integer): Integer; stdcall; external EpanetLib; - function ENsetcontrolenabled(Index: Integer; var enabled: Integer): Integer; stdcall; external EpanetLib; + function ENaddcontrol(Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single; var Index: Integer): Integer; cdecl; external EpanetLib; + function ENdeletecontrol(Index: Integer): Integer; cdecl; external EpanetLib; + function ENgetcontrol(Index: Integer; var Ctype: Integer; var Link: Integer; var Setting: Single; var Node: Integer; var Level: Single): Integer; cdecl; external EpanetLib; + function ENsetcontrol(Index: Integer; Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single): Integer; cdecl; external EpanetLib; + function ENgetcontrolenabled(Index: Integer; out_enabled: Integer): Integer; cdecl; external EpanetLib; + function ENsetcontrolenabled(Index: Integer; var enabled: Integer): Integer; cdecl; external EpanetLib; {Rule-Based Control Functions} - function ENaddrule(Rule: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENdeleterule(Index: Integer): Integer; stdcall; external EpanetLib; + function ENaddrule(Rule: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENdeleterule(Index: Integer): Integer; cdecl; external EpanetLib; function ENgetrule(Index: Integer; var Npremises: Integer; var NthenActions: Integer; - var NelseActions: Integer; var Priority: Single): Integer; stdcall; external EpanetLib; - function ENgetruleID(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; - function ENsetrulepriority(Index: Integer; Priority: Single): Integer; stdcall; external EpanetLib; + var NelseActions: Integer; var Priority: Single): Integer; cdecl; external EpanetLib; + function ENgetruleID(Index: Integer; ID: PAnsiChar): Integer; cdecl; external EpanetLib; + function ENsetrulepriority(Index: Integer; Priority: Single): Integer; cdecl; external EpanetLib; function ENgetpremise(RuleIndex: Integer; PremiseIndex: Integer; var LogOp: Integer; var ObjType: Integer; var ObjIndex: Integer; var Param: Integer; var RelOp: Integer; - var Status: Integer; var Value: Single): Integer; stdcall; external EpanetLib; + var Status: Integer; var Value: Single): Integer; cdecl; external EpanetLib; function ENsetpremise(RuleIndex: Integer; PremiseIndex: Integer; LogOp: Integer; ObjType: Integer; - ObjIndex: Integer; Param: Integer; RelOp: Integer; Status: Integer; Value: Single): Integer; stdcall; external EpanetLib; - function ENsetpremiseindex(RuleIndex: Integer; PremiseIndex: Integer; ObjIndex: Integer): Integer; stdcall; external EpanetLib; - function ENsetpremisestatus(RuleIndex: Integer; PremiseIndex: Integer; Status: Integer): Integer; stdcall; external EpanetLib; - function ENsetpremisevalue(RuleIndex: Integer; PremiseIndex: Integer; Value: Single): Integer; stdcall; external EpanetLib; + ObjIndex: Integer; Param: Integer; RelOp: Integer; Status: Integer; Value: Single): Integer; cdecl; external EpanetLib; + function ENsetpremiseindex(RuleIndex: Integer; PremiseIndex: Integer; ObjIndex: Integer): Integer; cdecl; external EpanetLib; + function ENsetpremisestatus(RuleIndex: Integer; PremiseIndex: Integer; Status: Integer): Integer; cdecl; external EpanetLib; + function ENsetpremisevalue(RuleIndex: Integer; PremiseIndex: Integer; Value: Single): Integer; cdecl; external EpanetLib; function ENgetthenaction(RuleIndex: Integer; ActionIndex: Integer; var LinkIndex: Integer; - var Status: Integer; var Setting: Single): Integer; stdcall; external EpanetLib; + var Status: Integer; var Setting: Single): Integer; cdecl; external EpanetLib; function ENsetthenaction(RuleIndex: Integer; ActionIndex: Integer; LinkIndex: Integer; - Status: Integer; Setting: Single): Integer; stdcall; external EpanetLib; + Status: Integer; Setting: Single): Integer; cdecl; external EpanetLib; function ENgetelseaction(RuleIndex: Integer; ActionIndex: Integer; var LinkIndex: Integer; - var Status: Integer; var Setting: Single): Integer; stdcall; external EpanetLib; + var Status: Integer; var Setting: Single): Integer; cdecl; external EpanetLib; function ENsetelseaction(RuleIndex: Integer; ActionIndex: Integer; LinkIndex: Integer; - Status: Integer; Setting: Single): Integer; stdcall; external EpanetLib; - function ENgetruleenabled(Index: Integer; out_enabled: Integer): Integer; stdcall; external EpanetLib; - function ENsetruleenabled(Index: Integer; var enabled: Integer): Integer; stdcall; external EpanetLib; + Status: Integer; Setting: Single): Integer; cdecl; external EpanetLib; + function ENgetruleenabled(Index: Integer; out_enabled: Integer): Integer; cdecl; external EpanetLib; + function ENsetruleenabled(Index: Integer; var enabled: Integer): Integer; cdecl; external EpanetLib; implementation diff --git a/include/epanet2.vb b/include/epanet2.vb index 73e46511..a7b31799 100644 --- a/include/epanet2.vb +++ b/include/epanet2.vb @@ -4,7 +4,7 @@ 'Declarations of functions in the EPANET PROGRAMMERs TOOLKIT '(EPANET2.DLL) for use with VB.Net. -'Last updated on 09/28/2023 +'Last updated on 06/23/2024 Imports System.Runtime.InteropServices Imports System.Text @@ -42,6 +42,9 @@ Public Const EN_CANOVERFLOW = 26 Public Const EN_DEMANDDEFICIT = 27 Public Const EN_NODE_INCONTROL = 28 Public Const EN_EMITTERFLOW = 29 +Public Const EN_LEAKAGEFLOW = 30 +Public Const EN_DEMANDFLOW = 31 +Public Const EN_FULLDEMAND = 32 Public Const EN_DIAMETER = 0 ' Link parameters Public Const EN_LENGTH = 1 @@ -69,6 +72,9 @@ Public Const EN_PUMP_EPAT = 22 Public Const EN_LINK_INCONTROL = 23 Public Const EN_GPV_CURVE = 24 Public Const EN_PCV_CURVE = 25 +Public Const EN_LEAK_AREA = 26 +Public Const EN_LEAK_EXPAN = 27 +Public Const EN_LINK_LEAKAGE = 28 Public Const EN_DURATION = 0 ' Time parameters Public Const EN_HYDSTEP = 1 @@ -93,6 +99,7 @@ Public Const EN_MAXFLOWCHANGE = 3 Public Const EN_MASSBALANCE = 4 Public Const EN_DEFICIENTNODES = 5 Public Const EN_DEMANDREDUCTION = 6 +Public Const EN_LEAKAGELOSS = 7 Public Const EN_NODE = 0 ' Component types Public Const EN_LINK = 1 @@ -309,6 +316,7 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENgeterror Lib "epanet2.dll" (ByVal errcode As Int32, ByVal errmsg As String, ByVal maxLen As Int32) As Int32 Declare Function ENgetstatistic Lib "epanet2.dll" (ByVal type_ As Int32, ByRef value As Single) As Int32 Declare Function ENgetresultindex Lib "epanet2.dll" (ByVal type_ As Int32, ByVal index As Int32, ByRef value As Int32) As Int32 + Declare Function ENtimetonextevent Lib "epanet2.dll" (eventType As Int32, duration As Int32, elementIndex As Int32) As Int32 'Analysis Options Functions Declare Function ENgetoption Lib "epanet2.dll" (ByVal option As Int32, value As Single) As Int32 @@ -334,7 +342,8 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENsettankdata Lib "epanet2.dll" (ByVal index As Int32, ByVal elev As Single, ByVal initlvl As Single, ByVal minlvl As Single, ByVal maxlvl As Single, ByVal diam As Single, ByVal minvol As Single, ByVal volcurve As String) As Int32 Declare Function ENgetcoord Lib "epanet2.dll" (ByVal index As Int32, x As Double, y As Double) As Int32 Declare Function ENsetcoord Lib "epanet2.dll" (ByVal index As Int32, ByVal x As Double, ByVal y As Double) As Int32 - + Declare Function ENgetnodevalues Lib "epanet2.dll" (ByVal property as Int32, values as Any) As Int32 + 'Nodal Demand Functions Declare Function ENgetdemandmodel Lib "epanet2.dll" (type_ As Int32, pmin As Single, preq As Single, pexp As Single) As Int32 Declare Function ENsetdemandmodel Lib "epanet2.dll" (ByVal type_ As Int32, ByVal pmin As Single, ByVal preq As Single, ByVal pexp As Single) As Int32 @@ -366,6 +375,7 @@ Public Const EN_TRUE = 1 ' boolean true Declare Function ENgetvertex Lib "epanet2.dll" (ByVal index As Int32, ByVal vertex As Int32, x As Double, y As Double) As Int32 Declare Function ENsetvertex Lib "epanet2.dll" (ByVal index As Int32, ByVal vertex As Int32, ByVal x As Double, ByVal y As Double) As Int32 Declare Function ENsetvertices Lib "epanet2.dll" (ByVal index As Int32, xCoords As Any, yCoords As Any, ByVal count As Int32) As Int32 + Declare Function ENgetlinkvalues Lib "epanet2.dll" (ByVal property as Int32, values as Any) As Int32 'Pump Functions Declare Function ENgetheadcurveindex Lib "epanet2.dll" (ByVal linkIndex As Int32, curveIndex As Int32) As Int32 diff --git a/include/epanet2_enums.h b/include/epanet2_enums.h index 557f48fe..4fc7755e 100644 --- a/include/epanet2_enums.h +++ b/include/epanet2_enums.h @@ -65,7 +65,10 @@ typedef enum { EN_CANOVERFLOW = 26, //!< Tank can overflow (= 1) or not (= 0) EN_DEMANDDEFICIT = 27,//!< Amount that full demand is reduced under PDA (read only) EN_NODE_INCONTROL = 28, //!< Is present in any simple or rule-based control (= 1) or not (= 0) - EN_EMITTERFLOW = 29 //!< Current emitter flow (read only) + EN_EMITTERFLOW = 29, //!< Current emitter flow (read only) + EN_LEAKAGEFLOW = 30, //!< Current leakage flow (read only) + EN_DEMANDFLOW = 31, //!< Current consumer demand delivered (read only) + EN_FULLDEMAND = 32 //!< Current consumer demand requested (read only) } EN_NodeProperty; /// Link properties @@ -99,7 +102,10 @@ typedef enum { EN_PUMP_EPAT = 22, //!< Pump energy price time pattern index EN_LINK_INCONTROL = 23, //!< Is present in any simple or rule-based control (= 1) or not (= 0) EN_GPV_CURVE = 24, //!< GPV head loss v. flow curve index - EN_PCV_CURVE = 25 //!< PCV loss coeff. curve index + EN_PCV_CURVE = 25, //!< PCV loss coeff. curve index + EN_LEAK_AREA = 26, //!< Pipe leak area (sq mm per 100 length units) + EN_LEAK_EXPAN = 27, //!< Leak expansion rate (sq mm per unit of pressure head) + EN_LINK_LEAKAGE = 28 //!< Current leakage rate (read only) } EN_LinkProperty; /// Time parameters @@ -152,7 +158,8 @@ typedef enum { EN_MAXFLOWCHANGE = 3, //!< Largest flow change in links EN_MASSBALANCE = 4, //!< Cumulative water quality mass balance ratio EN_DEFICIENTNODES = 5, //!< Number of pressure deficient nodes - EN_DEMANDREDUCTION = 6 //!< % demand reduction at pressure deficient nodes + EN_DEMANDREDUCTION = 6, //!< % demand reduction at pressure deficient nodes + EN_LEAKAGELOSS = 7 //!< % flow lost to system leakage } EN_AnalysisStatistic; /// Types of network objects diff --git a/src/enumstxt.h b/src/enumstxt.h index dbb4c74c..65962a6d 100755 --- a/src/enumstxt.h +++ b/src/enumstxt.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 02/05/2023 + Last Updated: 06/24/2024 ****************************************************************************** */ @@ -127,7 +127,7 @@ char *SectTxt[] = {s_TITLE, s_JUNCTIONS, s_RESERVOIRS, s_REACTIONS, s_MIXING, s_REPORT, s_TIMES, s_OPTIONS, s_COORDS, s_VERTICES, s_LABELS, s_BACKDROP, - s_TAGS, s_END, + s_TAGS, s_LEAKAGE, s_END, NULL}; char *Fldname[] = {t_ELEV, t_DEMAND, t_HEAD, diff --git a/src/epanet.c b/src/epanet.c index 991024cd..7e68c9d3 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 09/28/2023 + Last Updated: 06/26/2024 ****************************************************************************** */ @@ -510,6 +510,10 @@ int DLLEXPORT EN_initH(EN_Project p, int initFlag) return errcode; } } + + // Open pipe leakage modeling system + errcode = openleakage(p); + if (errcode) return errcode; // Initialize hydraulics solver inithyd(p, fflag); @@ -564,7 +568,11 @@ int DLLEXPORT EN_closeH(EN_Project p) */ { if (!p->Openflag) return 102; - if (p->hydraul.OpenHflag) closehyd(p); + if (p->hydraul.OpenHflag) + { + closeleakage(p); + closehyd(p); + } p->hydraul.OpenHflag = FALSE; return 0; } @@ -1044,6 +1052,9 @@ int DLLEXPORT EN_getstatistic(EN_Project p, int type, double *value) case EN_DEMANDREDUCTION: *value = p->hydraul.DemandReduction; break; + case EN_LEAKAGELOSS: + *value = p->hydraul.LeakageLoss; + break; case EN_MASSBALANCE: *value = p->quality.MassBalance.ratio; break; @@ -1864,8 +1875,10 @@ int DLLEXPORT EN_addnode(EN_Project p, const char *id, int nodeType, int *index) hyd->NodeDemand = (double *)realloc(hyd->NodeDemand, size); qual->NodeQual = (double *)realloc(qual->NodeQual, size); hyd->NodeHead = (double *)realloc(hyd->NodeHead, size); - hyd->DemandFlow = (double *)realloc(hyd->DemandFlow, size); + hyd->FullDemand = (double *)realloc(hyd->FullDemand, size); hyd->EmitterFlow = (double *)realloc(hyd->EmitterFlow, size); + hyd->LeakageFlow = (double *)realloc(hyd->LeakageFlow, size); + hyd->DemandFlow = (double *)realloc(hyd->DemandFlow, size); // Actions taken when a new Junction is added if (nodeType == EN_JUNCTION) @@ -2256,7 +2269,7 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val Ucf[VOLUME]; break; - case EN_DEMAND: + case EN_DEMAND: // Consumer Demand + Emitter Flow + Leakage Flow v = hyd->NodeDemand[index] * Ucf[FLOW]; break; @@ -2336,11 +2349,10 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val case EN_DEMANDDEFICIT: if (index > nJuncs) return 0; - // After an analysis, DemandFlow contains node's required demand - // while NodeDemand contains delivered demand + emitter flow - if (hyd->DemandFlow[index] < 0.0) return 0; - v = (hyd->DemandFlow[index] - - (hyd->NodeDemand[index] - hyd->EmitterFlow[index])) * Ucf[FLOW]; + // FullDemand contains node's required consumer demand + // while DemandFlow contains delivered consumer demand + if (hyd->FullDemand[index] <= 0.0) return 0; + v = (hyd->FullDemand[index] - hyd->DemandFlow[index]) * Ucf[FLOW]; break; case EN_NODE_INCONTROL: @@ -2350,6 +2362,18 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val case EN_EMITTERFLOW: v = hyd->EmitterFlow[index] * Ucf[FLOW]; break; + + case EN_LEAKAGEFLOW: + v = hyd->LeakageFlow[index] * Ucf[FLOW]; + break; + + case EN_DEMANDFLOW: // Consumer demand delivered + v = hyd->DemandFlow[index] * Ucf[FLOW]; + break; + + case EN_FULLDEMAND: // Consumer demand requested + v = hyd->FullDemand[index] * Ucf[FLOW]; + break; default: return 251; @@ -3352,6 +3376,8 @@ int DLLEXPORT EN_addlink(EN_Project p, const char *id, int linkType, } link->Kb = 0; link->Kw = 0; + link->LeakArea = 0; + link->LeakExpan = 0; link->R = 0; link->Rc = 0; link->Rpt = 0; @@ -3923,6 +3949,18 @@ int DLLEXPORT EN_getlinkvalue(EN_Project p, int index, int property, double *val v = (double)incontrols(p, LINK, index); break; + case EN_LEAK_AREA: + v = Link[index].LeakArea * Ucf[LENGTH]; + break; + + case EN_LEAK_EXPAN: + v = Link[index].LeakExpan * Ucf[LENGTH]; + break; + + case EN_LINK_LEAKAGE: + v = findlinkleakage(p, index) * Ucf[FLOW]; + break; + default: return 251; } @@ -4163,6 +4201,16 @@ int DLLEXPORT EN_setlinkvalue(EN_Project p, int index, int property, double valu } break; + case EN_LEAK_AREA: // leak area in sq mm per 100 pipe length units + if (value < 0.0) return 211; + Link[index].LeakArea = value / Ucf[LENGTH]; + break; + + case EN_LEAK_EXPAN: // leak area expansion slope (sq mm per unit of head) + if (value < 0.0) return 211; + Link[index].LeakExpan = value / Ucf[LENGTH]; + break; + default: return 251; } diff --git a/src/flowbalance.c b/src/flowbalance.c new file mode 100644 index 00000000..2c0efa66 --- /dev/null +++ b/src/flowbalance.c @@ -0,0 +1,191 @@ +/* + ****************************************************************************** + Project: OWA EPANET + Version: 2.3 + Module: flowbalance.c + Description: computes components of network's flow balance + Authors: see AUTHORS + Copyright: see AUTHORS + License: see LICENSE + Last Updated: 06/26/2024 + ****************************************************************************** +*/ + +#include "types.h" + +// Exported functions (declared in funcs.h) +//void startflowbalance(Project *); +//void updateflowbalance(Project *, long); +//void endflowbalance(Project *); + +void startflowbalance(Project *pr) +/* +**------------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: initializes components of the network's flow balance. +**------------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + hyd->FlowBalance.totalInflow = 0.0; + hyd->FlowBalance.totalOutflow = 0.0; + hyd->FlowBalance.consumerDemand = 0.0; + hyd->FlowBalance.emitterDemand = 0.0; + hyd->FlowBalance.leakageDemand = 0.0; + hyd->FlowBalance.deficitDemand = 0.0; + hyd->FlowBalance.storageDemand = 0.0; + hyd->FlowBalance.ratio = 0.0; +} + +void updateflowbalance(Project *pr, long hstep) +/* +**------------------------------------------------------------------- +** Input: hstep = time step (sec) +** Output: none +** Purpose: updates components of the system flow balance. +**------------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + Times *time = &pr->times; + + int i, j; + double v, dt, deficit, fullDemand; + SflowBalance flowBalance; + + // Determine current time interval in seconds + if (time->Dur == 0) dt = 1.0; + else if (time->Htime < time->Dur) + { + dt = (double) hstep; + } + else return; + + // Initialize local flow balance + flowBalance.totalInflow = 0.0; + flowBalance.totalOutflow = 0.0; + flowBalance.consumerDemand = 0.0; + flowBalance.emitterDemand = 0.0; + flowBalance.leakageDemand = 0.0; + flowBalance.deficitDemand = 0.0; + flowBalance.storageDemand = 0.0; + fullDemand = 0.0; + + // Initialize demand deficiency & leakage loss + hyd->DeficientNodes = 0; + hyd->DemandReduction = 0.0; + hyd->LeakageLoss = 0.0; + + // Examine each junction node + for (i = 1; i <= net->Njuncs; i++) + { + // Accumulate consumer demand flow + v = hyd->DemandFlow[i]; + if (v < 0.0) + flowBalance.totalInflow += (-v); + else + { + fullDemand += hyd->FullDemand[i]; + flowBalance.consumerDemand += v; + flowBalance.totalOutflow += v; + } + + // Accumulate emitter and leakage flow + v = hyd->EmitterFlow[i]; + flowBalance.emitterDemand += v; + flowBalance.totalOutflow += v; + v = hyd->LeakageFlow[i]; + flowBalance.leakageDemand += v; + flowBalance.totalOutflow += v; + + // Accumulate demand deficit flow + if (hyd->DemandModel == PDA && hyd->FullDemand[i] > 0.0) + { + deficit = hyd->FullDemand[i] - hyd->DemandFlow[i]; + if (deficit > TINY) + { + hyd->DeficientNodes++; + flowBalance.deficitDemand += deficit; + } + } + } + + // Examine each tank/reservoir node + for (j = 1; j <= net->Ntanks; j++) + { + 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) + { + if (v >= 0.0) + flowBalance.totalOutflow += v; + else + flowBalance.totalInflow += (-v); + } + + // For tank under extended period analysis + 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; + + // Update flow balance for entire run + hyd->FlowBalance.totalInflow += flowBalance.totalInflow * dt; + hyd->FlowBalance.totalOutflow += flowBalance.totalOutflow * dt; + hyd->FlowBalance.consumerDemand += flowBalance.consumerDemand * dt; + hyd->FlowBalance.emitterDemand += flowBalance.emitterDemand * dt; + hyd->FlowBalance.leakageDemand += flowBalance.leakageDemand * dt; + hyd->FlowBalance.deficitDemand += flowBalance.deficitDemand * dt; + hyd->FlowBalance.storageDemand += flowBalance.storageDemand * dt; +} + +void endflowbalance(Project *pr) +/* +**------------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: finalizes components of the system flow balance. +**------------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + Times *time = &pr->times; + + double seconds, qin, qout, qstor, r; + + if (time->Htime > 0) + seconds = time->Htime; + else + seconds = 1.0; + hyd->FlowBalance.totalInflow /= seconds; + hyd->FlowBalance.totalOutflow /= seconds; + hyd->FlowBalance.consumerDemand /= seconds; + hyd->FlowBalance.emitterDemand /= seconds; + hyd->FlowBalance.leakageDemand /= seconds; + hyd->FlowBalance.deficitDemand /= seconds; + hyd->FlowBalance.storageDemand /= seconds; + + qin = hyd->FlowBalance.totalInflow; + qout = hyd->FlowBalance.totalOutflow; + qstor = hyd->FlowBalance.storageDemand; + if (qstor > 0.0) + qout += qstor; + else + qin -= qstor; + if (qin == qout) + r = 1.0; + else if (qin > 0.0) + r = qout / qin; + else + r = 0.0; + hyd->FlowBalance.ratio = r; +} diff --git a/src/funcs.h b/src/funcs.h index 8910c361..87a6be24 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 09/28/2023 + Last Updated: 06/26/2024 ****************************************************************************** */ #ifndef FUNCS_H @@ -100,6 +100,7 @@ int controldata(Project *); int energydata(Project *); int sourcedata(Project *); int emitterdata(Project *); +int leakagedata(Project *); int qualdata(Project *); int reactdata(Project *); int mixingdata(Project *); @@ -142,6 +143,7 @@ void writecontrolaction(Project *, int, int); void writeruleaction(Project *, int, char *); int writehydwarn(Project *, int,double); void writehyderr(Project *, int); +void writeflowbalance(Project *); void writemassbalance(Project *); void writetime(Project *, char *); char *clocktime(char *, long); @@ -195,4 +197,19 @@ int savefinaloutput(Project *); int saveinpfile(Project *, const char *); +// ------- LEAKAGE.C -------------------- + +int openleakage(Project *); +void closeleakage(Project *); +double findlinkleakage(Project *, int); +void leakagecoeffs(Project *); +double leakageflowchange(Project *, int); +int leakagehasconverged(Project *); + +// ------- FLOWBALANCE.C----------------- + +void startflowbalance(Project *); +void updateflowbalance(Project *, long); +void endflowbalance(Project *); + #endif diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index c00d2f9f..5be8c403 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 03/29/2023 + Last Updated: 06/15/2024 ****************************************************************************** */ @@ -310,6 +310,7 @@ void matrixcoeffs(Project *pr) linkcoeffs(pr); emittercoeffs(pr); demandcoeffs(pr); + if (hyd->HasLeakage) leakagecoeffs(pr); // Update nodal flow balances with demands and add onto r.h.s. coeffs. nodecoeffs(pr); @@ -606,7 +607,7 @@ void demandheadloss(Project *pr, int i, double dp, double n, Hydraul *hyd = &pr->hydraul; double d = hyd->DemandFlow[i]; - double dfull = hyd->NodeDemand[i]; + double dfull = hyd->FullDemand[i]; double r = d / dfull; // Evaluate inverted demand function diff --git a/src/hydraul.c b/src/hydraul.c index 39769212..ae40ad63 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 09/28/2023 + Last Updated: 06/26/2024 ****************************************************************************** */ @@ -63,7 +63,7 @@ int openhyd(Project *pr) // Allocate memory for hydraulic variables ERRCODE(allocmatrix(pr)); - + // Check for unconnected nodes ERRCODE(unlinked(pr)); @@ -107,8 +107,10 @@ void inithyd(Project *pr, int initflag) hyd->OldStatus[net->Nlinks+i] = TEMPCLOSED; } - // Initialize emitter flows + // Initialize node outflows + memset(hyd->DemandFlow,0,(net->Nnodes+1)*sizeof(double)); memset(hyd->EmitterFlow,0,(net->Nnodes+1)*sizeof(double)); + memset(hyd->LeakageFlow,0,(net->Nnodes+1)*sizeof(double)); for (i = 1; i <= net->Nnodes; i++) { net->Node[i].ResultIndex = i; @@ -161,6 +163,9 @@ void inithyd(Project *pr, int initflag) pump->Energy.CurrentPower = 0.0; pump->Energy.CurrentEffic = 0.0; } + + // Initialize flow balance + startflowbalance(pr); // Re-position hydraulics file if (pr->outfile.Saveflag) @@ -253,6 +258,9 @@ int nexthyd(Project *pr, long *tstep) // Accumulate pumping energy if (time->Dur == 0) addenergy(pr,0); else if (time->Htime < time->Dur) addenergy(pr,hydstep); + + // Update flow balance + updateflowbalance(pr, hydstep); // More time remains - update current time if (time->Htime < time->Dur) @@ -267,6 +275,8 @@ int nexthyd(Project *pr, long *tstep) // No more time remains - force completion of analysis else { + endflowbalance(pr); + if (pr->report.Statflag) writeflowbalance(pr); time->Htime++; if (pr->quality.OpenQflag) time->Qtime++; } @@ -495,7 +505,7 @@ void demands(Project *pr) if (djunc > 0.0) hyd->Dsystem += djunc; sum += djunc; } - hyd->NodeDemand[i] = sum; + hyd->FullDemand[i] = sum; // Initialize pressure dependent demand hyd->DemandFlow[i] = sum; @@ -1147,4 +1157,3 @@ void resetpumpflow(Project *pr, int i) if (pump->Ptype == CONST_HP) pr->hydraul.LinkFlow[i] = pump->Q0; } - diff --git a/src/hydsolver.c b/src/hydsolver.c index 25c6752c..d88823fb 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -8,7 +8,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 02/14/2022 + Last Updated: 06/15/2024 ****************************************************************************** */ @@ -47,6 +47,7 @@ static double newflows(Project *, Hydbalance *); static void newlinkflows(Project *, Hydbalance *, double *, double *); static void newemitterflows(Project *, Hydbalance *, double *, double *); static void newdemandflows(Project *, Hydbalance *, double *, double *); +static void newleakageflows(Project *, Hydbalance *, double *, double *); static void checkhydbalance(Project *, Hydbalance *); static int hasconverged(Project *, double *, Hydbalance *); @@ -93,7 +94,6 @@ int hydsolve(Project *pr, int *iter, double *relerr) int valveChange; // Valve status change flag int statChange; // Non-valve status change flag Hydbalance hydbal; // Hydraulic balance errors - double fullDemand; // Full demand for a node (cfs) // Initialize status checking & relaxation factor nextcheck = hyd->CheckFreq; @@ -195,12 +195,12 @@ int hydsolve(Project *pr, int *iter, double *relerr) errcode = 110; } - // Store actual junction outflow in NodeDemand & full demand in DemandFlow + // Save total outflow (NodeDemand) at each junction for (i = 1; i <= net->Njuncs; i++) { - fullDemand = hyd->NodeDemand[i]; - hyd->NodeDemand[i] = hyd->DemandFlow[i] + hyd->EmitterFlow[i]; - hyd->DemandFlow[i] = fullDemand; + hyd->NodeDemand[i] = hyd->DemandFlow[i] + + hyd->EmitterFlow[i] + + hyd->LeakageFlow[i]; } // Save convergence info @@ -381,6 +381,7 @@ double newflows(Project *pr, Hydbalance *hbal) newlinkflows(pr, hbal, &qsum, &dqsum); newemitterflows(pr, hbal, &qsum, &dqsum); newdemandflows(pr, hbal, &qsum, &dqsum); + if (hyd->HasLeakage) newleakageflows(pr, hbal, &qsum, &dqsum); // Return ratio of total flow corrections to total flow if (qsum > hyd->Hacc) return (dqsum / qsum); @@ -514,6 +515,45 @@ void newemitterflows(Project *pr, Hydbalance *hbal, double *qsum, } +void newleakageflows(Project *pr, Hydbalance *hbal, double *qsum, + double *dqsum) +/* +**---------------------------------------------------------------- +** Input: hbal = ptr. to hydraulic balance information +** qsum = sum of current system flows +** dqsum = sum of system flow changes +** Output: updates hbal, qsum and dqsum +** Purpose: updates nodal leakage flows after new nodal heads computed +**---------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + int i; + double dq; + + for (i = 1; i <= net->Njuncs; i++) + { + // Update leakage flow at node i + dq = leakageflowchange(pr, i); + if (dq == 0.0) continue; + + // Update system flow summation + *qsum += ABS(hyd->LeakageFlow[i]); + *dqsum += ABS(dq); + + // Update identity of element with max. flow change + if (ABS(dq) > hbal->maxflowchange) + { + hbal->maxflowchange = ABS(dq); + hbal->maxflownode = i; + hbal->maxflowlink = -1; + } + } +} + + void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum) /* **---------------------------------------------------------------- @@ -546,7 +586,7 @@ void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum) for (i = 1; i <= net->Njuncs; i++) { // Skip junctions with no positive demand - if (hyd->NodeDemand[i] <= 0.0) continue; + if (hyd->FullDemand[i] <= 0.0) continue; // Find change in demand flow (see hydcoeffs.c) demandheadloss(pr, i, dp, n, &hloss, &hgrad); @@ -555,8 +595,8 @@ void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum) dq *= hyd->RelaxFactor; // Prevent a flow change greater than full demand - if (fabs(dq) > 0.4 * hyd->NodeDemand[i]) - dq = 0.4 * SGN(dq) * hyd->NodeDemand[i]; + if (fabs(dq) > 0.4 * hyd->FullDemand[i]) + dq = 0.4 * SGN(dq) * hyd->FullDemand[i]; hyd->DemandFlow[i] -= dq; // Update system flow summation @@ -641,11 +681,15 @@ int hasconverged(Project *pr, double *relerr, Hydbalance *hbal) if (hyd->FlowChangeLimit > 0.0 && hbal->maxflowchange > hyd->FlowChangeLimit) return 0; + // Check for node leakage convergence + if (hyd->HasLeakage && !leakagehasconverged(pr)) return 0; + // Check for pressure driven analysis convergence if (hyd->DemandModel == PDA) return pdaconverged(pr); return 1; } + int pdaconverged(Project *pr) /* **-------------------------------------------------------------- @@ -659,47 +703,32 @@ int pdaconverged(Project *pr) Hydraul *hyd = &pr->hydraul; const double QTOL = 0.0001; // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm) - int i, converged = 1; - double totalDemand = 0.0, totalReduction = 0.0; + int i; double dp = hyd->Preq - hyd->Pmin; double p, q, r; - hyd->DeficientNodes = 0; - hyd->DemandReduction = 0.0; - - // Add up number of junctions with demand deficits + // Examine each network junction for (i = 1; i <= pr->network.Njuncs; i++) { // Skip nodes whose required demand is non-positive - if (hyd->NodeDemand[i] <= 0.0) continue; + if (hyd->FullDemand[i] <= 0.0) continue; // Evaluate demand equation at current pressure solution p = hyd->NodeHead[i] - pr->network.Node[i].El; if (p <= hyd->Pmin) q = 0.0; else if (p >= hyd->Preq) - q = hyd->NodeDemand[i]; + q = hyd->FullDemand[i]; else { r = (p - hyd->Pmin) / dp; - q = hyd->NodeDemand[i] * pow(r, hyd->Pexp); + q = hyd->FullDemand[i] * pow(r, hyd->Pexp); } // Check if demand has not converged - if (fabs(q - hyd->DemandFlow[i]) > QTOL) - converged = 0; - - // Accumulate total required demand and demand deficit - if (hyd->DemandFlow[i] + QTOL < hyd->NodeDemand[i]) - { - hyd->DeficientNodes++; - totalDemand += hyd->NodeDemand[i]; - totalReduction += hyd->NodeDemand[i] - hyd->DemandFlow[i]; - } - } - if (totalDemand > 0.0) - hyd->DemandReduction = totalReduction / totalDemand * 100.0; - return converged; + if (fabs(q - hyd->DemandFlow[i]) > QTOL) return 0; + } + return 1; } diff --git a/src/inpfile.c b/src/inpfile.c index ec786b26..9945a2ae 100644 --- a/src/inpfile.c +++ b/src/inpfile.c @@ -7,7 +7,7 @@ Description: saves network data to an EPANET formatted text file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 05/11/2024 +Last Updated: 06/18/2024 ****************************************************************************** */ @@ -369,6 +369,19 @@ int saveinpfile(Project *pr, const char *fname) fprintf(f, "\n %-31s\t%-14.6f", node->ID, ke); } + // Write [LEAKAGE] section + fprintf(f, "\n\n"); + fprintf(f, s_LEAKAGE); + fprintf(f, "\n;;%-31s\t%-14s\t%-14s", + "Pipe", "Leak Area", "Leak Expansion"); + for (i = 1; i <= net->Nlinks; i++) + { + link = &net->Link[i]; + if (link->LeakArea == 0.0 && link->LeakExpan == 0.0) continue; + fprintf(f, "\n %-31s %14.6f %14.6f", link->ID, + link->LeakArea * pr->Ucf[LENGTH], link->LeakExpan * pr->Ucf[LENGTH]); + } + // Write [STATUS] section fprintf(f, "\n\n"); fprintf(f, s_STATUS); @@ -584,7 +597,7 @@ int saveinpfile(Project *pr, const char *fname) fprintf(f, "\n\n"); fprintf(f, s_REACTIONS); - fprintf(f, "\n;;%-9s\t%-31s\t%-12s", "Type", "Pipe/Tank", "Coefficient"); + fprintf(f, "\n;%-9s\t%-31s\t%-12s", "Type", "Pipe/Tank", "Coefficient"); // Pipe-specific parameters for (i = 1; i <= net->Nlinks; i++) diff --git a/src/input1.c b/src/input1.c index ee40d1ec..1074decd 100644 --- a/src/input1.c +++ b/src/input1.c @@ -7,7 +7,7 @@ Description: retrieves network data from an EPANET input file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 09/28/2023 +Last Updated: 06/15/2024 ****************************************************************************** */ @@ -593,6 +593,10 @@ void convertunits(Project *pr) // Convert units on reaction coeffs. link->Kb /= SECperDAY; link->Kw /= SECperDAY; + + // Convert leakage parameters + link->LeakArea /= pr->Ucf[LENGTH]; + link->LeakExpan /= pr->Ucf[LENGTH]; } else if (link->Type == PUMP) diff --git a/src/input2.c b/src/input2.c index c5af2bae..f1ce8bfe 100644 --- a/src/input2.c +++ b/src/input2.c @@ -310,6 +310,7 @@ int newline(Project *pr, int sect, char *line) else return 0; case _SOURCES: return (sourcedata(pr)); case _EMITTERS: return (emitterdata(pr)); + case _LEAKAGE: return (leakagedata(pr)); case _QUALITY: return (qualdata(pr)); case _STATUS: return (statusdata(pr)); case _ROUGHNESS: return (0); diff --git a/src/input3.c b/src/input3.c index 1d4f03a7..296e752e 100644 --- a/src/input3.c +++ b/src/input3.c @@ -387,6 +387,8 @@ int pipedata(Project *pr) link->Km = 0.0; link->Kb = MISSING; link->Kw = MISSING; + link->LeakArea = 0.0; + link->LeakExpan = 0.0; link->Type = PIPE; link->Status = OPEN; link->Rpt = 0; @@ -494,6 +496,8 @@ int pumpdata(Project *pr) link->Km = 0.0; link->Kb = 0.0; link->Kw = 0.0; + link->LeakArea = 0.0; + link->LeakExpan = 0.0; link->Type = PUMP; link->Status = OPEN; link->Rpt = 0; @@ -613,6 +617,8 @@ int valvedata(Project *pr) link->Km = 0.0; link->Kb = 0.0; link->Kw = 0.0; + link->LeakArea = 0.0; + link->LeakExpan = 0.0; link->Type = type; link->Status = ACTIVE; link->Rpt = 0; @@ -1118,6 +1124,41 @@ int emitterdata(Project *pr) return 0; } +int leakagedata(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns error code +** Purpose: processes link leakage data +** Format: +** [LEAKAGE] +** link C1 C2 +**-------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Parser *parser = &pr->parser; + + int j, // Link index + n; // # data items + double c1, c2; // Flow coeff. + + // Check that link exists & is a pipe + n = parser->Ntokens; + if (n < 3) return 201; + if ((j = findlink(net, parser->Tok[0])) == 0) return setError(parser, 0, 203); + if (net->Link[j].Type > PIPE) return 0; + + // Parse leakage coeffs. + if (!getfloat(parser->Tok[1], &c1)) return setError(parser, 1, 202); + if (c1 < 0.0) return setError(parser, 1, 209); + if (!getfloat(parser->Tok[2], &c2)) return setError(parser, 2, 202); + if (c2 < 0.0) return setError(parser, 1, 209); + net->Link[j].LeakArea = c1; + net->Link[j].LeakExpan = c2; + return 0; +} + int qualdata(Project *pr) /* **-------------------------------------------------------------- diff --git a/src/leakage.c b/src/leakage.c new file mode 100644 index 00000000..cb5cad71 --- /dev/null +++ b/src/leakage.c @@ -0,0 +1,531 @@ +/* + ****************************************************************************** + Project: OWA EPANET + Version: 2.3 + Module: leakage.c + Description: models additional nodal demands due to pipe leaks. + Authors: see AUTHORS + Copyright: see AUTHORS + License: see LICENSE + Last Updated: 06/14/2024 + ****************************************************************************** +*/ +/* +This module uses the FAVAD (Fixed and Variable Discharge) equation to model +leaky pipes: + + Q = Co * L * (Ao + m * H) * sqrt(H) + +where Q = 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. + +The inverted form of this equation is used to model the leakage demand from +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, +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. + +In implementing this model, the pipe property "LeakArea" represents Ao in +sq. mm per 100 units of pipe length and "LeakExpan" represents m in sq. mm +per unit of pressure head. + +*/ +#include +#include + +#include "types.h" +#include "funcs.h" + +// Exported functions (declared in funcs.h) +//int openleakage(Project *); +//void closeleakage(Project *); +//double findlinkleakage(Project *, int); +//void leakagecoeffs(Project *); +//double leakageflowchange(Project *, int); +//int leakagehasconverged(Project *); + +// Local functions +static int check_for_leakage(Project *pr); +static int create_leakage_objects(Project *pr); +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); + + +int openleakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: returns an error code +** Purpose: opens the pipe leakage modeling system +**------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + + int err; + + // Check if project includes leakage + closeleakage(pr); + hyd->HasLeakage = check_for_leakage(pr); + if (!hyd->HasLeakage) return 0; + + // Allocate memory for leakage data objects + err = create_leakage_objects(pr); + if (err > 0) return err; + + // Convert pipe leakage coeffs. to node coeffs. + convert_pipe_to_node_leakage(pr); + init_node_leakage(pr); + return 0; +} + + +int check_for_leakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: returns TRUE if any pipes can leak, FALSE otherwise +** Purpose: checks if any pipes can leak. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + int i; + Slink *link; + + for (i = 1; i <= net->Nlinks; i++) + { + // Only pipes have leakage + link = &net->Link[i]; + if (link->Type > PIPE) continue; + if (link->LeakArea > 0.0 || link->LeakExpan > 0.0) return TRUE; + } + return FALSE; +} + + +int create_leakage_objects(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: returns an error code +** Purpose: allocates an array of Leakage objects. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + int i; + + hyd->Leakage = (Sleakage *)calloc(net->Njuncs + 1, sizeof(Sleakage)); + if (hyd->Leakage == NULL) return 101; + for (i = 1; i <= net->Njuncs; i++) + { + hyd->Leakage[i].cfa = 0.0; + hyd->Leakage[i].cva = 0.0; + hyd->Leakage[i].qfa = 0.0; + hyd->Leakage[i].qva = 0.0; + } + return 0; +} + +void convert_pipe_to_node_leakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: converts pipe leakage parameters into node leakage +** coefficents. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + int i; + double c_area, c_expan, c_orif, len; + Slink *link; + Snode *node1; + Snode *node2; + + // Examine each link + c_orif = 4.8149866 * 1.e-6; + for (i = 1; i <= net->Nlinks; i++) + { + // Only pipes have leakage + link = &net->Link[i]; + if (link->Type > PIPE) continue; + + // Ignore leakage in a pipe connecting two tanks or + // reservoirs (since those nodes don't have demands) + node1 = &net->Node[link->N1]; + node2 = &net->Node[link->N2]; + if (node1->Type != JUNCTION && node2->Type != JUNCTION) continue; + + // Get pipe's fixed and variable area leak coeffs. + if (link->LeakArea == 0.0 && link->LeakExpan == 0.0) continue; + c_area = c_orif * link->LeakArea / SQR(MperFT); + c_expan = c_orif * link->LeakExpan; + + // Adjust for number of 100-ft pipe sections + len = link->Len * pr->Ucf[LENGTH] / 100.; + if (node1->Type == JUNCTION && node2->Type == JUNCTION) + { + len *= 0.5; + } + c_area *= len; + c_expan *= len; + + // Add these coeffs. to pipe's end nodes + if (node1->Type == JUNCTION) + { + hyd->Leakage[link->N1].cfa += c_area; + hyd->Leakage[link->N1].cva += c_expan; + } + if (node2->Type == JUNCTION) + { + hyd->Leakage[link->N2].cfa += c_area; + hyd->Leakage[link->N2].cva += c_expan; + } + } +} + +void init_node_leakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: initializes node leakage coeffs. and flows. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + int i; + double c_area, c_expan; + + for (i = 1; i <= net->Njuncs; i++) + { + // Coeff. for fixed area leakage + c_area = hyd->Leakage[i].cfa; + if (c_area > 0.0) + hyd->Leakage[i].cfa = 1.0 / (c_area * c_area); + else + hyd->Leakage[i].cfa = 0.0; + + // Coeff. for variable area leakage + c_expan = hyd->Leakage[i].cva; + if (c_expan > 0.0) + hyd->Leakage[i].cva = 1.0 / pow(c_expan, 2./3.); + else + hyd->Leakage[i].cva = 0.0; + + // Initialize leakage flow to a non-zero value (as required by + // the hydraulic solver) + if (hyd->Leakage[i].cfa > 0.0) + hyd->Leakage[i].qfa = 0.001; + if (hyd->Leakage[i].cva > 0.0) + hyd->Leakage[i].qva = 0.001; + } +} + +void closeleakage(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: frees memory for nodal leakage objects. +**------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + if (hyd->Leakage) free(hyd->Leakage); + hyd->Leakage = NULL; + hyd->HasLeakage = FALSE; +} + +double findlinkleakage(Project *pr, int i) +/*------------------------------------------------------------- +** Input: i = link index +** Output: returns link leakage flow (cfs) +** Purpose: computes leakage flow from link i at current +** hydraulic solution. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + Smatrix *sm = &hyd->smatrix; + Slink *link = &net->Link[i]; + + int n1, n2; + double h1, h2, hsqrt, a, m, c, len, q1, q2; + + // Only pipes can leak + link = &net->Link[i]; + if (link->Type > PIPE) return 0.0; + + // No leakage if area & expansion are 0 + if (link->LeakArea == 0.0 && link->LeakExpan == 0.0) return 0.0; + + // No leakage if link's end nodes are both fixed grade + n1 = link->N1; + n2 = link->N2; + if (n1 > net->Njuncs && n2 > net->Njuncs) return 0.0; + + // Pressure head of end nodes + h1 = hyd->NodeHead[n1] - net->Node[n1].El; + h1 = MAX(h1, 0.0); + h2 = hyd->NodeHead[n2] - net->Node[n2].El; + h2 = MAX(h2, 0.0); + + // Pipe leak parameters converted to feet + a = link->LeakArea / SQR(MperFT); + m = link->LeakExpan; + len = link->Len * pr->Ucf[LENGTH] / 100.; // # 100 ft pipe lengths + c = 4.8149866 * len / 2.0 * 1.e-6; + + // Leakage from 1st half of pipe connected to node n1 + q1 = 0.0; + if (n1 <= net->Njuncs) + { + hsqrt = sqrt(h1); + q1 = c * (a + m * h1) * hsqrt; + } + + // Leakage from 2nd half of pipe connected to node n2 + q2 = 0.0; + if (n2 <= net->Njuncs) + { + hsqrt = sqrt(h2); + q2 = c * (a + m * h2) * hsqrt; + } + + // Adjust leakage flows to account for one node being fixed grade + if (q2 == 0.0) q1 *= 2.0; + if (q1 == 0.0) q2 *= 2.0; + return q1 + q2; +} + +void leakagecoeffs(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes coeffs. of the linearized hydraulic eqns. +** contributed by node leakages. +**-------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + Smatrix *sm = &hyd->smatrix; + + int i, row; + double hfa, // head loss producing current fixed area leakage + gfa, // gradient of fixed area head loss + hva, // head loss producing current variable area leakage + gva; // gradient of variable area head loss + + Snode* node; + + for (i = 1; i <= net->Njuncs; i++) + { + // Skip junctions that don't leak + node = &net->Node[i]; + if (!leakage_headloss(pr, i, &hfa, &gfa, &hva, &gva)) continue; + + // Addition to matrix diagonal & r.h.s + row = sm->Row[i]; + if (gfa > 0.0) + { + sm->Aii[row] += 1.0 / gfa; + sm->F[row] += (hfa + node->El) / gfa; + } + if (gva > 0.0) + { + sm->Aii[row] += 1.0 / gva; + sm->F[row] += (hva + node->El) / gva; + } + + // Update node's flow excess (inflow - outflow) + hyd->Xflow[i] -= (hyd->Leakage[i].qfa + hyd->Leakage[i].qva); + } +} + +double leakageflowchange(Project *pr, int i) +/* +**-------------------------------------------------------------- +** Input: i = node index +** Output: returns change in leakage flow rate +** Purpose: finds new leakage flow rate at a node after new +** heads are computed by the hydraulic solver. +**-------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + double hfa, gfa, hva, gva, // same as defined in leakage_solvercoeffs() + dh, dqfa, dqva; + + // 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; + + // GGA flow update formula for fixed area leakage + dqfa = 0.0; + if (gfa > 0.0) + { + dqfa = (hfa - dh) / 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; + hyd->Leakage[i].qva -= dqva; + } + + // New leakage flow at the node + hyd->LeakageFlow[i] = hyd->Leakage[i].qfa + hyd->Leakage[i].qva; + return dqfa + dqva; +} + +int leakagehasconverged(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns TRUE if leakage calculations converged, +** FALSE if not +** Purpose: checks if leakage calculations have converged. +**-------------------------------------------------------------- +*/ +{ + 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; + + for (i = 1; i <= net->Njuncs; i++) + { + // Skip junctions that don't leak + if (hyd->Leakage[i].cfa == 0 && hyd->Leakage[i].cva == 0) continue; + + // Evaluate node's pressure head + h = hyd->NodeHead[i] - net->Node[i].El; + + // Directly compute a reference leakage at this pressure head + qref = 0.0; + // Contribution from pipes with fixed area leaks + if (hyd->Leakage[i].cfa > 0.0) + qref = sqrt(h / hyd->Leakage[i].cfa); + // Contribution from pipes with variable area leaks + if (hyd->Leakage[i].cva > 0.0) + qref += pow((h / hyd->Leakage[i].cva), 1.5); + + // Compare reference leakage to solution leakage + qtest = hyd->Leakage[i].qfa + hyd->Leakage[i].qva; + if (fabs(qref - qtest) > ABSTOL + RELTOL * qref) return FALSE; + } + return TRUE; +} + +int leakage_headloss(Project* pr, int i, double *hfa, double *gfa, + double *hva, double *gva) +/* +**-------------------------------------------------------------- +** Input: i = node index +** Output: hfa = fixed area leak head loss (ft) +** gfa = gradient of fixed area head loss (ft/cfs) +** hva = variable area leak head loss (ft) +** gva = gradient of variable area head loss (ft/cfs) +** returns TRUE if node has leakage, FALSE otherwise +** Purpose: finds head loss and its gradient for a node's +** leakage as a function of leakage flow. +**-------------------------------------------------------------- +*/ +{ + 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); + 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); + return TRUE; +} + +void eval_node_leakage(double RQtol, double q, double c, double n, + double *h, double *g) +/* +**-------------------------------------------------------------- +** 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) +** Purpose: evaluates inverted form of leakage equation to +** compute head loss and its gradient as a function +** flow. +**-------------------------------------------------------------- +*/ +{ + n = 1.0 / n; + *g = n * c * pow(fabs(q), n - 1.0); + + // Use linear head loss function for small gradient + if (*g < RQtol) + { + *g = RQtol / n; + *h = (*g) * q; + } + + // Otherwise use normal leakage head loss function + else *h = (*g) * q / n; + + // Prevent leakage from going negative + add_lower_barrier(q, h, g); +} + +void add_lower_barrier(double q, double* h, double* g) +/* +**-------------------------------------------------------------------- +** Input: q = current flow rate +** Output: h = head loss value +** g = 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); +} diff --git a/src/project.c b/src/project.c index 2104d060..80c76979 100644 --- a/src/project.c +++ b/src/project.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 09/28/2023 + Last Updated: 06/24/2024 ****************************************************************************** */ @@ -328,8 +328,11 @@ void initpointers(Project *pr) pr->hydraul.P = NULL; pr->hydraul.Y = NULL; pr->hydraul.Xflow = NULL; + pr->hydraul.FullDemand = NULL; pr->hydraul.DemandFlow = NULL; pr->hydraul.EmitterFlow = NULL; + pr->hydraul.LeakageFlow = NULL; + pr->hydraul.Leakage = NULL; pr->quality.NodeQual = NULL; pr->quality.PipeRateCoeff = NULL; @@ -392,14 +395,18 @@ int allocdata(Project *pr) pr->hydraul.NodeDemand = (double *)calloc(n, sizeof(double)); pr->hydraul.NodeHead = (double *)calloc(n, sizeof(double)); pr->quality.NodeQual = (double *)calloc(n, sizeof(double)); + pr->hydraul.FullDemand = (double *)calloc(n, sizeof(double)); pr->hydraul.DemandFlow = (double *)calloc(n, sizeof(double)); pr->hydraul.EmitterFlow = (double *)calloc(n, sizeof(double)); + pr->hydraul.LeakageFlow = (double *)calloc(n, sizeof(double)); ERRCODE(MEMCHECK(pr->network.Node)); ERRCODE(MEMCHECK(pr->hydraul.NodeDemand)); ERRCODE(MEMCHECK(pr->hydraul.NodeHead)); ERRCODE(MEMCHECK(pr->quality.NodeQual)); + ERRCODE(MEMCHECK(pr->hydraul.FullDemand)); ERRCODE(MEMCHECK(pr->hydraul.DemandFlow)); ERRCODE(MEMCHECK(pr->hydraul.EmitterFlow)); + ERRCODE(MEMCHECK(pr->hydraul.LeakageFlow)); } // Allocate memory for network links @@ -471,8 +478,10 @@ void freedata(Project *pr) free(pr->hydraul.LinkFlow); free(pr->hydraul.LinkSetting); free(pr->hydraul.LinkStatus); + free(pr->hydraul.FullDemand); free(pr->hydraul.DemandFlow); free(pr->hydraul.EmitterFlow); + free(pr->hydraul.LeakageFlow); free(pr->quality.NodeQual); // Free memory used for nodal adjacency lists diff --git a/src/report.c b/src/report.c index ffdc2d22..25bb1376 100644 --- a/src/report.c +++ b/src/report.c @@ -422,6 +422,44 @@ void writehydstat(Project *pr, int iter, double relerr) writeline(pr, " "); } +void writeflowbalance(Project *pr) +/* +**------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: writes hydraulic flow balance ratio to report file. +**------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + Report *rpt = &pr->report; + char s1[MAXMSG+1]; + double ucf = pr->Ucf[FLOW]; + + snprintf(s1, MAXMSG, "Hydraulic Flow Balance (%s)", rpt->Field[DEMAND].Units); + writeline(pr, s1); + snprintf(s1, MAXMSG, "================================"); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Total Inflow: %12.3f", hyd->FlowBalance.totalInflow*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Consumer Demand: %12.3f", hyd->FlowBalance.consumerDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Demand Deficit: %12.3f", hyd->FlowBalance.deficitDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Emitter Flow: %12.3f", hyd->FlowBalance.emitterDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Leakage Flow: %12.3f", hyd->FlowBalance.leakageDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Total Outflow: %12.3f", hyd->FlowBalance.totalOutflow*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Storage Flow: %12.3f", hyd->FlowBalance.storageDemand*ucf); + writeline(pr, s1); + snprintf(s1, MAXMSG, "Flow Ratio: %12.3f", hyd->FlowBalance.ratio); + writeline(pr, s1); + snprintf(s1, MAXMSG, "================================\n"); + writeline(pr, s1); +} + void writemassbalance(Project *pr) /* **------------------------------------------------------------- diff --git a/src/text.h b/src/text.h index 1a78573b..84b05db7 100755 --- a/src/text.h +++ b/src/text.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 05/11/2024 + Last Updated: 06/15/2024 ****************************************************************************** */ @@ -151,6 +151,9 @@ #define w_REQUIRED "REQ" #define w_EXPONENT "EXP" +#define w_AREA "AREA" +#define w_EXPAN "EXPAN" + #define w_SECONDS "SEC" #define w_MINUTES "MIN" #define w_HOURS "HOU" @@ -208,6 +211,7 @@ #define s_DEMANDS "[DEMANDS]" #define s_SOURCES "[SOURCES]" #define s_EMITTERS "[EMITTERS]" +#define s_LEAKAGE "[LEAKAGE]" #define s_PATTERNS "[PATTERNS]" #define s_CURVES "[CURVES]" #define s_QUALITY "[QUALITY]" diff --git a/src/types.h b/src/types.h index ff7c3599..30fd64e8 100755 --- a/src/types.h +++ b/src/types.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 07/17/2023 + Last Updated: 06/15/2024 ****************************************************************************** */ @@ -294,7 +294,7 @@ typedef enum { _VALVES, _CONTROLS, _RULES, _DEMANDS, _SOURCES, _EMITTERS, _PATTERNS, _CURVES, _QUALITY, _STATUS, _ROUGHNESS, _ENERGY, _REACTIONS, _MIXING, _REPORT, _TIMES, _OPTIONS, - _COORDS, _VERTICES, _LABELS, _BACKDROP, _TAGS, _END + _COORDS, _VERTICES, _LABELS, _BACKDROP, _TAGS, _LEAKAGE, _END } SectionType; typedef enum { @@ -413,6 +413,8 @@ typedef struct // Link Object double Kw; // wall react. coef. double R; // flow resistance double Rc; // reaction coeff. + double LeakArea; // leak area (sq mm per 100 pipe length units + double LeakExpan; // leak expansion (sq mm per unit of head) LinkType Type; // link type StatusType Status; // initial status Pvertices Vertices; // internal vertex coordinates @@ -549,6 +551,26 @@ typedef struct // Mass Balance Components double ratio; // ratio of mass added to mass lost } SmassBalance; +typedef struct +{ + double totalInflow; + double totalOutflow; + double consumerDemand; + double emitterDemand; + double leakageDemand; + double deficitDemand; + double storageDemand; + double ratio; +} SflowBalance; + +typedef struct // Node Leakage Object +{ + double qfa; // fixed area leakage flow + double qva; // variable area leakage flow + double cfa; // fixed area leakage coeff. + double cva; // variable area leakage coeff. +} Sleakage; + /* ------------------------------------------------------ Wrapper Data Structures @@ -712,9 +734,11 @@ typedef struct { double *NodeHead, // Node hydraulic heads - *NodeDemand, // Node demand + emitter flows - *DemandFlow, // Work array of demand flows - *EmitterFlow, // Emitter outflows + *NodeDemand, // Node total demand (consumer + emitter + leakage) + *FullDemand, // Required consumer demand + *DemandFlow, // Demand flow from nodes + *EmitterFlow, // Emitter flow from nodes + *LeakageFlow, // Leakage flow from nodes *LinkFlow, // Link flows *LinkSetting, // Link settings Htol, // Hydraulic head tolerance @@ -741,6 +765,7 @@ typedef struct { MaxHeadError, // Max. error for link head loss MaxFlowChange, // Max. change in link flow DemandReduction, // % demand reduction at pressure deficient nodes + LeakageLoss, // % system leakage loss RelaxFactor, // Relaxation factor for flow updating *P, // Inverse of head loss derivatives *Y, // Flow correction factors @@ -759,12 +784,18 @@ typedef struct { MaxCheck, // Hydraulic trials limit on status checks OpenHflag, // Hydraulic system opened flag Haltflag, // Flag to halt simulation - DeficientNodes; // Number of pressure deficient nodes + DeficientNodes, // Number of pressure deficient nodes + HasLeakage; // TRUE if project has non-zero leakage parameters + + Sleakage *Leakage; // Array of node leakage parameters StatusType *LinkStatus, // Link status *OldStatus; // Previous link/tank status + SflowBalance + FlowBalance; // Flow balance components + Smatrix smatrix; // Sparse matrix storage } Hydraul; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d90167f7..744a58aa 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -41,6 +41,7 @@ set(toolkit_test_srcs test_pda.cpp test_valve.cpp test_units.cpp + test_leakage.cpp ) add_executable(test_toolkit ${toolkit_test_srcs}) diff --git a/tests/test_leakage.cpp b/tests/test_leakage.cpp new file mode 100644 index 00000000..4a1c7161 --- /dev/null +++ b/tests/test_leakage.cpp @@ -0,0 +1,93 @@ +/* + ****************************************************************************** + Project: OWA EPANET + Version: 2.3 + Module: test_leakage.cpp + Description: Tests EPANET toolkit api functions + Authors: see AUTHORS + Copyright: see AUTHORS + License: see LICENSE + Last Updated: 06/26/2024 + ****************************************************************************** +*/ + +/* + Tests Pipe Leakage Feature +*/ + +#include + +#include "test_toolkit.hpp" + +BOOST_AUTO_TEST_SUITE (test_leakage) + +BOOST_AUTO_TEST_CASE(test_leakage_model) + +//#include +//#include +//#include "epanet2_2.h" + +//int main() +{ + int error = 0; + int Pipe21, Junc21, Junc22; + double pipe21Leak, junc21Leak, junc22Leak; + EN_Project ph = NULL; + + error = EN_createproject(&ph); + BOOST_REQUIRE(error == 0); + error = EN_open(ph, DATA_PATH_NET1, DATA_PATH_RPT, ""); +// error = EN_open(ph, "Net1.inp", "Net1.rpt", ""); + BOOST_REQUIRE(error == 0); + + // single period analysis + error = EN_settimeparam(ph, EN_DURATION, 0); + BOOST_REQUIRE(error == 0); + + // Get index of Pipe 21 + error = EN_getlinkindex(ph, "21", &Pipe21); + BOOST_REQUIRE(error == 0); + + // Set Pipe21 leak area to 1.0 sq mm per 100 ft of pipe + // and its expansion rate to 0.1 sq mm per ft of head + error = EN_setlinkvalue(ph, Pipe21, EN_LEAK_AREA, 1.0); + BOOST_REQUIRE(error == 0); + error = EN_setlinkvalue(ph, Pipe21, EN_LEAK_EXPAN, 0.1); + BOOST_REQUIRE(error == 0); + + // Solve for hydraulics + error = EN_solveH(ph); + BOOST_REQUIRE(error == 0); + + // Compute Pipe 21 leakage flow using the FAVAD formula + // Note: we can't just sum the leak rates at both end nodes + // together since in general the nodes can have leakage + // contributed by other connecting pipes. + error = EN_getlinkvalue(ph, Pipe21, EN_LINK_LEAKAGE, &pipe21Leak); + BOOST_REQUIRE(error == 0); +// printf("\n Pipe leakage flow: %.4f", pipe21Leak); + + // Retrieve leakage flow at end nodes + // Note: In this case all of the leakage at these nodes is from Pipe 21. + error = EN_getnodeindex(ph, "21", &Junc21); + BOOST_REQUIRE(error == 0); + error = EN_getnodeindex(ph, "22", &Junc22); + BOOST_REQUIRE(error == 0); + error = EN_getnodevalue(ph, Junc21, EN_LEAKAGEFLOW, &junc21Leak); + BOOST_REQUIRE(error == 0); + error = EN_getnodevalue(ph, Junc22, EN_LEAKAGEFLOW, &junc22Leak); + BOOST_REQUIRE(error == 0); + + // Check that the sum of the node leakages equals the pipe leakage + //printf("\n Node leakage flow: %.4f\n", junc21Leak + junc22Leak); + BOOST_REQUIRE(abs(pipe21Leak - (junc21Leak+junc22Leak)) < 0.01); + +// Clean up + error = EN_close(ph); + BOOST_REQUIRE(error == 0); + error = EN_deleteproject(ph); + BOOST_REQUIRE(error == 0); +// return 0; +} + +BOOST_AUTO_TEST_SUITE_END() From 30a3fcb4b981d0a776123da028fe1f314bbd321f Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Wed, 26 Jun 2024 11:47:35 -0400 Subject: [PATCH 3/5] debug test_pda.cpp --- tests/test_pda.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test_pda.cpp b/tests/test_pda.cpp index 6091def7..7e951c0d 100644 --- a/tests/test_pda.cpp +++ b/tests/test_pda.cpp @@ -79,6 +79,9 @@ BOOST_AUTO_TEST_CASE(test_pda_model) error = EN_getnodeindex(ph, (char *)"21", &index); BOOST_REQUIRE(error == 0); error = EN_getnodevalue(ph, index, EN_DEMANDDEFICIT, &reduction); + + printf("\nreduction = %f", reduction); + BOOST_REQUIRE(error == 0); BOOST_REQUIRE(abs(reduction - 413.67) < 0.01); From 68b73a14f1a20daa292da8d40554ef2c40b74baf Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Wed, 26 Jun 2024 17:50:24 -0400 Subject: [PATCH 4/5] Fixes for PDA --- src/flowbalance.c | 25 ++++++++++--------------- src/hydsolver.c | 26 +++++++++++++++++++++----- tests/test_pda.cpp | 3 --- 3 files changed, 31 insertions(+), 23 deletions(-) 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/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/tests/test_pda.cpp b/tests/test_pda.cpp index 7e951c0d..6091def7 100644 --- a/tests/test_pda.cpp +++ b/tests/test_pda.cpp @@ -79,9 +79,6 @@ BOOST_AUTO_TEST_CASE(test_pda_model) error = EN_getnodeindex(ph, (char *)"21", &index); BOOST_REQUIRE(error == 0); error = EN_getnodevalue(ph, index, EN_DEMANDDEFICIT, &reduction); - - printf("\nreduction = %f", reduction); - BOOST_REQUIRE(error == 0); BOOST_REQUIRE(abs(reduction - 413.67) < 0.01); From 6089b93a51223dc8458e68355b5ba97913595b14 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Tue, 2 Jul 2024 20:19:45 -0400 Subject: [PATCH 5/5] Fix refactoring error in hydcoeffs.c --- src/hydcoeffs.c | 2 +- src/leakage.c | 75 +++++++++++++++++++++++++++---------------------- 2 files changed, 42 insertions(+), 35 deletions(-) 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/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); }