From b20859f2af0e7339c6dcc9d7e6467ceb0adfa9e2 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Tue, 25 Jun 2024 13:17:57 -0400 Subject: [PATCH] Pipe Leakage Added --- ReleaseNotes2_3.md | 19 +- include/epanet2.bas | 11 +- include/epanet2.cs | 25 ++- include/epanet2.def | 17 +- include/epanet2.h | 4 + include/epanet2.pas | 306 ++++++++++++++------------ include/epanet2.vb | 14 +- include/epanet2_2.h | 25 +++ include/epanet2_enums.h | 13 +- src/enumstxt.h | 4 +- src/epanet.c | 101 ++++++++- src/epanet2.c | 27 +++ src/flowbalance.c | 191 +++++++++++++++++ src/funcs.h | 19 ++ src/hydcoeffs.c | 5 +- src/hydraul.c | 24 ++- src/hydsolver.c | 93 +++++--- src/inpfile.c | 141 +++++++----- src/input1.c | 6 +- src/input2.c | 1 + src/input3.c | 51 ++++- src/leakage.c | 465 ++++++++++++++++++++++++++++++++++++++++ src/project.c | 11 +- src/quality.c | 1 + src/report.c | 40 +++- src/text.h | 9 +- src/types.h | 43 +++- tests/CMakeLists.txt | 1 + tests/test_leakage.cpp | 93 ++++++++ 29 files changed, 1488 insertions(+), 272 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 67680a8a..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,4 +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.h b/include/epanet2.h index 5cc6a415..18a7345f 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -214,6 +214,8 @@ extern "C" { int DLLEXPORT ENgetnodevalue(int index, int property, EN_API_FLOAT_TYPE *value); + int DLLEXPORT ENgetnodevalues(int property, EN_API_FLOAT_TYPE *value); + int DLLEXPORT ENsetnodevalue(int index, int property, EN_API_FLOAT_TYPE value); int DLLEXPORT ENsetjuncdata(int index, EN_API_FLOAT_TYPE elev, @@ -291,6 +293,8 @@ extern "C" { int DLLEXPORT ENgetlinkvalue(int index, int property, EN_API_FLOAT_TYPE *value); + int DLLEXPORT ENgetlinkvalues(int property, EN_API_FLOAT_TYPE *value); + int DLLEXPORT ENsetlinkvalue(int index, int property, EN_API_FLOAT_TYPE value); int DLLEXPORT ENsetpipedata(int index, EN_API_FLOAT_TYPE length, 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_2.h b/include/epanet2_2.h index 5a3c3af2..1094280d 100644 --- a/include/epanet2_2.h +++ b/include/epanet2_2.h @@ -881,6 +881,20 @@ typedef struct Project *EN_Project; Values are returned in units that depend on the units used for flow rate (see @ref Units). */ + + int DLLEXPORT EN_getnodevalues(EN_Project ph, int property, double *out_values); + + /** + @brief Retrieves an array of property values for all nodes. + @param ph an EPANET project handle. + @param property the property to retrieve (see @ref EN_NodeProperty). + @param[out] values an array of values for all nodes. + @return an error code. + + Values are returned in units that depend on the units used for flow rate + (see @ref Units). + */ + int DLLEXPORT EN_getnodevalue(EN_Project ph, int index, int property, double *out_value); /** @@ -1242,6 +1256,17 @@ typedef struct Project *EN_Project; */ int DLLEXPORT EN_getlinkvalue(EN_Project ph, int index, int property, double *out_value); + /** + @brief Retrieves an array of property values for all links. + @param ph an EPANET project handle. + @param property the property to retrieve (see @ref EN_LinkProperty). + @param[out] values an array of values for all links. + @return an error code. + + Values are returned in units that depend on the units used for flow rate (see @ref Units). + */ + int DLLEXPORT EN_getlinkvalues(EN_Project ph, int property, double *out_values); + /** @brief Sets a property value for a link. @param ph an EPANET project handle. 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 86b510ba..98a18df9 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -1044,6 +1044,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,7 +1867,7 @@ 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); // Actions taken when a new Junction is added @@ -2256,7 +2259,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 +2339,13 @@ 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]) / + hyd->FullDemand[index]; + if (v < TINY) v = 0.0; + v *= 100.0; break; case EN_NODE_INCONTROL: @@ -2350,6 +2355,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; @@ -2358,6 +2375,26 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val return 0; } +int DLLEXPORT EN_getnodevalues(EN_Project p, int property, double *values) +/*---------------------------------------------------------------- +** Input: property = node property code (see EN_NodeProperty) +** Output: values = array of node property values +** Returns: error code +** Purpose: retrieves an array of node property values +**---------------------------------------------------------------- +*/ +{ + int status = 0; + + for (int i = 1; i <= p->network.Nnodes; i++) + { + status = EN_getnodevalue(p, i, property, &values[i - 1]); + // if status is not 0, return the error code + if (status != 0) { return status; } + } + return 0; +} + int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double value) /*---------------------------------------------------------------- ** Input: index = node index @@ -2429,6 +2466,7 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu if (value < 0.0) return 209; if (value > 0.0) value = pow((Ucf[FLOW] / value), hyd->Qexp) / Ucf[PRESSURE]; Node[index].Ke = value; + if (hyd->EmitterFlow[index] == 0.0) hyd->EmitterFlow[index] = 1.0; break; case EN_INITQUAL: @@ -3331,6 +3369,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; @@ -3408,6 +3448,12 @@ int DLLEXPORT EN_deletelink(EN_Project p, int index, int actionCode) if (net->Valve[i].Link > index) net->Valve[i].Link -= 1; } + // Reduce the number of pipes count by one if it is a pipe. + if (linkType == PIPE) + { + net->Npipes--; + } + // Delete any pump associated with the deleted link if (linkType == PUMP) { @@ -3896,6 +3942,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 = leakage_getlinkleakage(p, index) * Ucf[FLOW]; + break; + default: return 251; } @@ -3903,6 +3961,25 @@ int DLLEXPORT EN_getlinkvalue(EN_Project p, int index, int property, double *val return 0; } +int DLLEXPORT EN_getlinkvalues(EN_Project p, int property, double *values) +/*---------------------------------------------------------------- +** Input: property = link property code (see EN_LinkProperty) +** Output: values = array of link property values +** Returns: error code +** Purpose: retrieves property values for all links +**---------------------------------------------------------------- +*/ +{ + int status = 0; + for(int i = 1; i <= p->network.Nlinks; i++) + { + status = EN_getlinkvalue(p, i, property, &values[i-1]); + // If an error occurs, return the error code + if(status != 0) { return status; } + } + return 0; +} + int DLLEXPORT EN_setlinkvalue(EN_Project p, int index, int property, double value) /*---------------------------------------------------------------- ** Input: index = link index @@ -4117,6 +4194,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/epanet2.c b/src/epanet2.c index 7d7cd50f..c3c8fa29 100644 --- a/src/epanet2.c +++ b/src/epanet2.c @@ -355,6 +355,20 @@ int DLLEXPORT ENgetnodevalue(int index, int property, EN_API_FLOAT_TYPE *value) return errcode; } +int DLLEXPORT ENgetnodevalues(int property, EN_API_FLOAT_TYPE *values) +{ + int i, errcode = 0; + EN_API_FLOAT_TYPE value; + + for (i = 1; i <= _defaultProject->network.Nnodes; i++) + { + errcode = ENgetnodevalue(i, property, &value); + values[i-1] = value; + if (errcode != 0) return errcode; + } + return 0; +} + int DLLEXPORT ENsetnodevalue(int index, int property, EN_API_FLOAT_TYPE value) { return EN_setnodevalue(_defaultProject, index, property, value); @@ -523,6 +537,19 @@ int DLLEXPORT ENgetlinkvalue(int index, int property, EN_API_FLOAT_TYPE *value) *value = (EN_API_FLOAT_TYPE)v; return errcode; } +int DLLEXPORT ENgetlinkvalues(int property, EN_API_FLOAT_TYPE *values) +{ + int i, errcode = 0; + EN_API_FLOAT_TYPE value; + + for (i = 1; i <= _defaultProject->network.Nlinks; i++) + { + errcode = ENgetlinkvalue(i, property, &value); + values[i-1] = value; + if (errcode != 0) return errcode; + } + return 0; +} int DLLEXPORT ENsetlinkvalue(int index, int property, EN_API_FLOAT_TYPE value) { diff --git a/src/flowbalance.c b/src/flowbalance.c new file mode 100644 index 00000000..a562c32d --- /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/14/2024 + ****************************************************************************** +*/ + +#include "types.h" + +// Exported functions (declared in funcs.h) +//void flowbalance_start(Project *); +//void flowbalance_update(Project *, long); +//void flowbalance_end(Project *); + +void flowbalance_start(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 flowbalance_update(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 flowbalance_end(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..119fd5e8 100755 --- a/src/funcs.h +++ b/src/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,21 @@ int savefinaloutput(Project *); int saveinpfile(Project *, const char *); +// ------- LEAKAGE.C -------------------- + +int leakage_open(Project *); +void leakage_init(Project *); +void leakage_close(Project *); + +double leakage_getlinkleakage(Project *, int); +void leakage_solvercoeffs(Project *); +double leakage_getflowchange(Project *, int); +int leakage_hasconverged(Project *); + +// ------- FLOWBALANCE.C----------------- + +void flowbalance_start(Project *); +void flowbalance_update(Project *, long); +void flowbalance_end(Project *); + #endif diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index c00d2f9f..21ac973e 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) leakage_solvercoeffs(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 7d1718b0..dd735544 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/15/2024 ****************************************************************************** */ @@ -40,6 +40,7 @@ void addenergy(Project *, long); void tanklevels(Project *, long); void resetpumpflow(Project *, int); void getallpumpsenergy(Project *); +void initleakcoeffs(Project *, int); int openhyd(Project *pr) /* @@ -63,6 +64,9 @@ int openhyd(Project *pr) // Allocate memory for hydraulic variables ERRCODE(allocmatrix(pr)); + + // Allocate memory for leakage analysis + ERRCODE(leakage_open(pr)); // Check for unconnected nodes ERRCODE(unlinked(pr)); @@ -107,13 +111,16 @@ 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; if (net->Node[i].Ke > 0.0) hyd->EmitterFlow[i] = 1.0; } + leakage_init(pr); // Initialize links for (i = 1; i <= net->Nlinks; i++) @@ -161,6 +168,9 @@ void inithyd(Project *pr, int initflag) pump->Energy.CurrentPower = 0.0; pump->Energy.CurrentEffic = 0.0; } + + // Initialize flow balance + flowbalance_start(pr); // Re-position hydraulics file if (pr->outfile.Saveflag) @@ -253,6 +263,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 + flowbalance_update(pr, hydstep); // More time remains - update current time if (time->Htime < time->Dur) @@ -267,6 +280,8 @@ int nexthyd(Project *pr, long *tstep) // No more time remains - force completion of analysis else { + flowbalance_end(pr); + if (pr->report.Statflag) writeflowbalance(pr); time->Htime++; if (pr->quality.OpenQflag) time->Qtime++; } @@ -284,8 +299,10 @@ void closehyd(Project *pr) **-------------------------------------------------------------- */ { + leakage_close(pr); freesparse(pr); freematrix(pr); + freeadjlists(&pr->network); } @@ -494,7 +511,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; @@ -1146,4 +1163,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..2886e9c3 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 = leakage_getflowchange(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 && !leakage_hasconverged(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 e23284d3..9945a2ae 100644 --- a/src/inpfile.c +++ b/src/inpfile.c @@ -1,13 +1,13 @@ /* ****************************************************************************** Project: OWA EPANET -Version: 2.2 +Version: 2.3 Module: inpfile.c Description: saves network data to an EPANET formatted text file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 04/30/2023 +Last Updated: 06/18/2024 ****************************************************************************** */ @@ -152,40 +152,46 @@ int saveinpfile(Project *pr, const char *fname) // (Leave demands for [DEMANDS] section) fprintf(f, "\n\n"); fprintf(f, s_JUNCTIONS); + fprintf(f, "\n;;%-31s\t%-12s\t%-12s\t%-31s", + "ID", "Elev", "Demand", "Pattern"); for (i = 1; i <= net->Njuncs; i++) { node = &net->Node[i]; - fprintf(f, "\n %-31s %12.4f", node->ID, node->El * pr->Ucf[ELEV]); - if (node->Comment) fprintf(f, " ;%s", node->Comment); + fprintf(f, "\n %-31s\t%-12.4f", node->ID, node->El * pr->Ucf[ELEV]); + if (node->Comment) fprintf(f, "\t;%s", node->Comment); } // Write [RESERVOIRS] section fprintf(f, "\n\n"); fprintf(f, s_RESERVOIRS); + fprintf(f, "\n;;%-31s\t%-12s\t%-31s", + "ID", "Head", "Pattern"); for (i = 1; i <= net->Ntanks; i++) { tank = &net->Tank[i]; if (tank->A == 0.0) { node = &net->Node[tank->Node]; - sprintf(s, " %-31s %12.4f", node->ID, node->El * pr->Ucf[ELEV]); - if ((j = tank->Pat) > 0) sprintf(s1, " %s", net->Pattern[j].ID); + sprintf(s, " %-31s\t%-12.4f", node->ID, node->El * pr->Ucf[ELEV]); + if ((j = tank->Pat) > 0) sprintf(s1, "%s", net->Pattern[j].ID); else strcpy(s1, " "); - fprintf(f, "\n%s %-31s", s, s1); - if (node->Comment) fprintf(f, " ;%s", node->Comment); + fprintf(f, "\n%s\t%-31s", s, s1); + if (node->Comment) fprintf(f, "\t;%s", node->Comment); } } // Write [TANKS] section fprintf(f, "\n\n"); fprintf(f, s_TANKS); + fprintf(f, "\n;;%-31s\t%-12s\t%-12s\t%-12s\t%-12s\t%-12s\t%-12s\t%-31s\t%-12s", + "ID", "Elevation", "InitLevel", "MinLevel", "MaxLevel", "Diameter", "MinVol", "VolCurve", "Overflow"); for (i = 1; i <= net->Ntanks; i++) { tank = &net->Tank[i]; if (tank->A > 0.0) { node = &net->Node[tank->Node]; - sprintf(s, " %-31s %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f", + sprintf(s, " %-31s\t%-12.4f\t%-12.4f\t%-12.4f\t%-12.4f\t%-12.4f\t%-12.4f", node->ID, node->El * pr->Ucf[ELEV], (tank->H0 - node->El) * pr->Ucf[ELEV], (tank->Hmin - node->El) * pr->Ucf[ELEV], @@ -195,15 +201,17 @@ int saveinpfile(Project *pr, const char *fname) if ((j = tank->Vcurve) > 0) sprintf(s1, "%s", net->Curve[j].ID); else if (tank->CanOverflow) strcpy(s1, "*"); else strcpy(s1, " "); - fprintf(f, "\n%s %-31s", s, s1); - if (tank->CanOverflow) fprintf(f, " YES "); - if (node->Comment) fprintf(f, " ;%s", node->Comment); + fprintf(f, "\n%s\t%-31s", s, s1); + if (tank->CanOverflow) fprintf(f, "\t%-12s", "YES"); + if (node->Comment) fprintf(f, "\t;%s", node->Comment); } } // Write [PIPES] section fprintf(f, "\n\n"); fprintf(f, s_PIPES); + fprintf(f, "\n;;%-31s\t%-31s\t%-31s\t%-12s\t%-12s\t%-12s\t%-12s\t%-6s", + "ID", "Node1", "Node2", "Length", "Diameter", "Roughness", "MinorLoss", "Status"); for (i = 1; i <= net->Nlinks; i++) { link = &net->Link[i]; @@ -214,36 +222,38 @@ int saveinpfile(Project *pr, const char *fname) if (hyd->Formflag == DW) kc = kc * pr->Ucf[ELEV] * 1000.0; km = link->Km * SQR(d) * SQR(d) / 0.02517; - sprintf(s, " %-31s %-31s %-31s %12.4f %12.4f %12.4f %12.4f", + sprintf(s, " %-31s\t%-31s\t%-31s\t%-12.4f\t%-12.4f\t%-12.4f\t%-12.4f", link->ID, net->Node[link->N1].ID, net->Node[link->N2].ID, link->Len * pr->Ucf[LENGTH], d * pr->Ucf[DIAM], kc, km); if (link->Type == CVPIPE) sprintf(s2, "CV"); else if (link->Status == CLOSED) sprintf(s2, "CLOSED"); else strcpy(s2, " "); - fprintf(f, "\n%s %-6s", s, s2); - if (link->Comment) fprintf(f, " ;%s", link->Comment); + fprintf(f, "\n%s\t%-6s", s, s2); + if (link->Comment) fprintf(f, "\t;%s", link->Comment); } } // Write [PUMPS] section fprintf(f, "\n\n"); fprintf(f, s_PUMPS); + fprintf(f, "\n;;%-31s\t%-31s\t%-31s\t%-12s", + "ID", "Node1", "Node2", "Parameters"); for (i = 1; i <= net->Npumps; i++) { n = net->Pump[i].Link; link = &net->Link[n]; pump = &net->Pump[i]; - sprintf(s, " %-31s %-31s %-31s", link->ID, net->Node[link->N1].ID, + sprintf(s, " %-31s\t%-31s\t%-31s", link->ID, net->Node[link->N1].ID, net->Node[link->N2].ID); // Pump has constant power - if (pump->Ptype == CONST_HP) sprintf(s1, " POWER %.4f", link->Km); + if (pump->Ptype == CONST_HP) sprintf(s1, "\tPOWER %.4f", link->Km); // Pump has a head curve else if ((j = pump->Hcurve) > 0) { - sprintf(s1, " HEAD %s", net->Curve[j].ID); + sprintf(s1, "\tHEAD %s", net->Curve[j].ID); } // Old format used for pump curve @@ -260,25 +270,26 @@ int saveinpfile(Project *pr, const char *fname) // Optional speed pattern if ((j = pump->Upat) > 0) { - sprintf(s1, " PATTERN %s", net->Pattern[j].ID); + sprintf(s1, "\tPATTERN %s", net->Pattern[j].ID); strcat(s, s1); } // Optional speed setting if (link->Kc != 1.0) { - sprintf(s1, " SPEED %.4f", link->Kc); + sprintf(s1, "\tSPEED %.4f", link->Kc); strcat(s, s1); } fprintf(f, "\n%s", s); - if (link->Comment) fprintf(f, " ;%s", link->Comment); - + if (link->Comment) fprintf(f, "\t;%s", link->Comment); } // Write [VALVES] section fprintf(f, "\n\n"); fprintf(f, s_VALVES); + fprintf(f, "\n;;%-31s\t%-31s\t%-31s\t%-12s\t%-6s\t%-12s\t%-12s", + "ID", "Node1", "Node2", "Diameter", "Type", "Setting", "MinorLoss"); for (i = 1; i <= net->Nvalves; i++) { n = net->Valve[i].Link; @@ -303,7 +314,7 @@ int saveinpfile(Project *pr, const char *fname) } km = link->Km * SQR(d) * SQR(d) / 0.02517; - sprintf(s, " %-31s %-31s %-31s %12.4f %5s", + sprintf(s, " %-31s\t%-31s\t%-31s\t%-12.4f\t%-6s", link->ID, net->Node[link->N1].ID, net->Node[link->N2].ID, d * pr->Ucf[DIAM], LinkTxt[link->Type]); @@ -311,22 +322,24 @@ int saveinpfile(Project *pr, const char *fname) // For GPV, setting = head curve index if (link->Type == GPV && (j = ROUND(link->Kc)) > 0) { - sprintf(s1, "%-31s %12.4f", net->Curve[j].ID, km); + sprintf(s1, "%-31s\t%-12.4f", net->Curve[j].ID, km); } // For PCV add loss curve if present else if (link->Type == PCV && (j = net->Valve[i].Curve) > 0) { - sprintf(s1, "%12.4f %12.4f %-31s", kc, km, net->Curve[j].ID); + sprintf(s1, "%-12.4f\t%-12.4f\t%-31s", kc, km, net->Curve[j].ID); } - else sprintf(s1, "%12.4f %12.4f", kc, km); - fprintf(f, "\n%s %s", s, s1); - if (link->Comment) fprintf(f, " ;%s", link->Comment); + else sprintf(s1, "%-12.4f\t%-12.4f", kc, km); + fprintf(f, "\n%s\t%s", s, s1); + if (link->Comment) fprintf(f, "\t;%s", link->Comment); } // Write [DEMANDS] section fprintf(f, "\n\n"); fprintf(f, s_DEMANDS); + fprintf(f, "\n;;%-31s\t%-14s\t%-31s\t%-31s", + "Junction", "Demand", "Pattern", "Category"); ucf = pr->Ucf[DEMAND]; for (i = 1; i <= net->Njuncs; i++) { @@ -334,11 +347,11 @@ int saveinpfile(Project *pr, const char *fname) for (demand = node->D; demand != NULL; demand = demand->next) { if (demand->Base == 0.0) continue; - sprintf(s, " %-31s %14.6f", node->ID, ucf * demand->Base); - if ((j = demand->Pat) > 0) sprintf(s1, " %-31s", net->Pattern[j].ID); + sprintf(s, " %-31s\t%-14.6f", node->ID, ucf * demand->Base); + if ((j = demand->Pat) > 0) sprintf(s1, "%-31s", net->Pattern[j].ID); else strcpy(s1, " "); - fprintf(f, "\n%s %-31s", s, s1); - if (demand->Name) fprintf(f, " ;%s", demand->Name); + fprintf(f, "\n%s\t%-31s", s, s1); + if (demand->Name) fprintf(f, "\t;%s", demand->Name); } } @@ -346,17 +359,34 @@ int saveinpfile(Project *pr, const char *fname) // Write [EMITTERS] section fprintf(f, "\n\n"); fprintf(f, s_EMITTERS); + fprintf(f, "\n;;%-31s\t%-14s", + "Junction", "Coefficient"); for (i = 1; i <= net->Njuncs; i++) { node = &net->Node[i]; if (node->Ke == 0.0) continue; ke = pr->Ucf[FLOW] / pow(pr->Ucf[PRESSURE] * node->Ke, (1.0 / hyd->Qexp)); - fprintf(f, "\n %-31s %14.6f", node->ID, ke); + 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); + fprintf(f, "\n;;%-31s\t%-12s", + "ID", "Status/Setting"); for (i = 1; i <= net->Nlinks; i++) { link = &net->Link[i]; @@ -364,7 +394,7 @@ int saveinpfile(Project *pr, const char *fname) { if (link->Status == CLOSED) { - fprintf(f, "\n %-31s %s", link->ID, StatTxt[CLOSED]); + fprintf(f, "\n %-31s\t%s", link->ID, StatTxt[CLOSED]); } // Write pump speed here for pumps with old-style pump curve input @@ -375,7 +405,7 @@ int saveinpfile(Project *pr, const char *fname) if (pump->Hcurve == 0 && pump->Ptype != CONST_HP && link->Kc != 1.0) { - fprintf(f, "\n %-31s %-.4f", link->ID, link->Kc); + fprintf(f, "\n %-31s\t%-.4f", link->ID, link->Kc); } } } @@ -385,11 +415,11 @@ int saveinpfile(Project *pr, const char *fname) { if (link->Status == OPEN) { - fprintf(f, "\n %-31s %s", link->ID, StatTxt[OPEN]); + fprintf(f, "\n %-31s\t%s", link->ID, StatTxt[OPEN]); } if (link->Status == CLOSED) { - fprintf(f, "\n%-31s %s", link->ID, StatTxt[CLOSED]); + fprintf(f, "\n%-31s\t%s", link->ID, StatTxt[CLOSED]); } } } @@ -398,26 +428,30 @@ int saveinpfile(Project *pr, const char *fname) // (Use 6 pattern factors per line) fprintf(f, "\n\n"); fprintf(f, s_PATTERNS); + fprintf(f, "\n;;%-31s\t%-12s", + "ID", "Multipliers"); for (i = 1; i <= net->Npats; i++) { if (net->Pattern[i].Comment) fprintf(f, "\n;%s", net->Pattern[i].Comment); for (j = 0; j < net->Pattern[i].Length; j++) { if (j % 6 == 0) fprintf(f, "\n %-31s", net->Pattern[i].ID); - fprintf(f, " %12.4f", net->Pattern[i].F[j]); + fprintf(f, "\t%-12.4f", net->Pattern[i].F[j]); } } // Write [CURVES] section fprintf(f, "\n\n"); fprintf(f, s_CURVES); + fprintf(f, "\n;;%-31s\t%-12s\t%-12s", + "ID", "X-Value", "Y-Value"); for (i = 1; i <= net->Ncurves; i++) { if (net->Curve[i].Comment) fprintf(f, "\n;%s", net->Curve[i].Comment); for (j = 0; j < net->Curve[i].Npts; j++) { curve = &net->Curve[i]; - fprintf(f, "\n %-31s %12.4f %12.4f", curve->ID, curve->X[j], curve->Y[j]); + fprintf(f, "\n %-31s\t%-12.4f\t%-12.4f", curve->ID, curve->X[j], curve->Y[j]); } } @@ -502,39 +536,42 @@ int saveinpfile(Project *pr, const char *fname) // (Skip nodes with default quality of 0) fprintf(f, "\n\n"); fprintf(f, s_QUALITY); + fprintf(f, "\n;;%-31s\t%-14s", "ID", "InitQual"); for (i = 1; i <= net->Nnodes; i++) { node = &net->Node[i]; if (node->C0 == 0.0) continue; - fprintf(f, "\n %-31s %14.6f", node->ID, node->C0 * pr->Ucf[QUALITY]); + fprintf(f, "\n %-31s\t%-14.6f", node->ID, node->C0 * pr->Ucf[QUALITY]); } // Write [SOURCES] section fprintf(f, "\n\n"); fprintf(f, s_SOURCES); + fprintf(f, "\n;;%-31s\t%-9s\t%-14s\t%-31s", "ID", "Type", "Quality", "Pattern"); for (i = 1; i <= net->Nnodes; i++) { node = &net->Node[i]; source = node->S; if (source == NULL) continue; - sprintf(s, " %-31s %-8s %14.6f", node->ID, SourceTxt[source->Type], + sprintf(s, " %-31s\t%-9s\t%-14.6f", node->ID, SourceTxt[source->Type], source->C0); if ((j = source->Pat) > 0) { sprintf(s1, "%s", net->Pattern[j].ID); } else strcpy(s1, ""); - fprintf(f, "\n%s %s", s, s1); + fprintf(f, "\n%s\t%s", s, s1); } // Write [MIXING] section fprintf(f, "\n\n"); fprintf(f, s_MIXING); + fprintf(f, "\n;;%-31s\t%-8s", "ID", "Model"); for (i = 1; i <= net->Ntanks; i++) { tank = &net->Tank[i]; if (tank->A == 0.0) continue; - fprintf(f, "\n %-31s %-8s %12.4f", net->Node[tank->Node].ID, + fprintf(f, "\n %-31s\t%-8s\t%12.4f", net->Node[tank->Node].ID, MixTxt[tank->MixModel], tank->V1frac); } @@ -558,6 +595,10 @@ int saveinpfile(Project *pr, const char *fname) fprintf(f, "\n ROUGHNESS CORRELATION %-.6f", qual->Rfactor); } + fprintf(f, "\n\n"); + fprintf(f, s_REACTIONS); + fprintf(f, "\n;%-9s\t%-31s\t%-12s", "Type", "Pipe/Tank", "Coefficient"); + // Pipe-specific parameters for (i = 1; i <= net->Nlinks; i++) { @@ -565,11 +606,11 @@ int saveinpfile(Project *pr, const char *fname) if (link->Type > PIPE) continue; if (link->Kb != qual->Kbulk) { - fprintf(f, "\n BULK %-31s %-.6f", link->ID, link->Kb * SECperDAY); + fprintf(f, "\n %-9s\t%-31s\t%-.6f", "BULK", link->ID, link->Kb * SECperDAY); } if (link->Kw != qual->Kwall) { - fprintf(f, "\n WALL %-31s %-.6f", link->ID, link->Kw * SECperDAY); + fprintf(f, "\n %-9s\t%-31s\t%-.6f", "WALL", link->ID, link->Kw * SECperDAY); } } @@ -580,7 +621,7 @@ int saveinpfile(Project *pr, const char *fname) if (tank->A == 0.0) continue; if (tank->Kb != qual->Kbulk) { - fprintf(f, "\n TANK %-31s %-.6f", net->Node[tank->Node].ID, + fprintf(f, "\n %-9s\t%-31s\t%-.6f", "TANK", net->Node[tank->Node].ID, tank->Kb * SECperDAY); } } @@ -682,7 +723,7 @@ int saveinpfile(Project *pr, const char *fname) fprintf(f, "\n PATTERN %s", net->Pattern[hyd->DefPat].ID); fprintf(f, "\n DEMAND MULTIPLIER %-.4f", hyd->Dmult); fprintf(f, "\n EMITTER EXPONENT %-.4f", 1.0 / hyd->Qexp); - fprintf(f, "\n EMITTER BACKFLOW %s", BackflowTxt[hyd->EmitBackFlag]); + fprintf(f, "\n BACKFLOW ALLOWED %s", BackflowTxt[hyd->EmitBackFlag]); fprintf(f, "\n VISCOSITY %-.6f", hyd->Viscos / VISCOS); fprintf(f, "\n DIFFUSIVITY %-.6f", qual->Diffus / DIFFUS); fprintf(f, "\n SPECIFIC GRAVITY %-.6f", hyd->SpGrav); @@ -793,23 +834,25 @@ int saveinpfile(Project *pr, const char *fname) // Write [COORDINATES] section fprintf(f, "\n\n"); fprintf(f, s_COORDS); + fprintf(f, "\n;;%-31s\t%-14s\t%-14s", "ID", "X-Coord", "Y-Coord"); for (i = 1; i <= net->Nnodes; i++) { node = &net->Node[i]; if (node->X == MISSING || node->Y == MISSING) continue; - fprintf(f, "\n %-31s %14.6f %14.6f", node->ID, node->X, node->Y); + fprintf(f, "\n %-31s\t%-14.6f\t%-14.6f", node->ID, node->X, node->Y); } // Write [VERTICES] section fprintf(f, "\n\n"); fprintf(f, s_VERTICES); + fprintf(f, "\n;;%-31s\t%-14s\t%-14s", "ID", "X-Coord", "Y-Coord"); for (i = 1; i <= net->Nlinks; i++) { link = &net->Link[i]; if (link->Vertices != NULL) { for (j = 0; j < link->Vertices->Npts; j++) - fprintf(f, "\n %-31s %14.6f %14.6f", + fprintf(f, "\n %-31s\t%-14.6f\t%-14.6f", link->ID, link->Vertices->X[j], link->Vertices->Y[j]); } } 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 8c10c10e..296e752e 100644 --- a/src/input3.c +++ b/src/input3.c @@ -7,7 +7,7 @@ Description: parses network data from a line of an EPANET input file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 09/28/2023 +Last Updated: 05/11/2024 ****************************************************************************** */ @@ -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) /* **-------------------------------------------------------------- @@ -1832,7 +1873,7 @@ int optionchoice(Project *pr, int n) ** UNBALANCED STOP/CONTINUE {Niter} ** PATTERN id ** DEMAND MODEL DDA/PDA -** EMITTER BACKFLOW YES/NO +** BACKFLOW ALLOWED YES/NO **-------------------------------------------------------------- */ { @@ -1956,11 +1997,11 @@ int optionchoice(Project *pr, int n) hyd->DemandModel = choice; } - // EMITTER BACKFLOW - else if (match(parser->Tok[0], w_EMITTER)) + // Emitter BACKFLOW ALLOWED + else if (match(parser->Tok[0], w_BACKFLOW)) { if (n < 2) return 0; - if (!match(parser->Tok[1], w_BACKFLOW)) return -1; + if (!match(parser->Tok[1], w_ALLOWED)) return -1; choice = findmatch(parser->Tok[2], BackflowTxt); if (choice < 0) return setError(parser, 2, 213); hyd->EmitBackFlag = choice; diff --git a/src/leakage.c b/src/leakage.c new file mode 100644 index 00000000..3752a83f --- /dev/null +++ b/src/leakage.c @@ -0,0 +1,465 @@ +/* + ****************************************************************************** + 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" + +// Exported functions (declared in funcs.h) +//int leakage_open(Project *); +//void leakage_init(Project *); +//void leakage_close(Project *); +//double leakage_getlinkleakage(Project *, int); +//void leakage_solvercoeffs(Project *); +//double leakage_getflowchange(Project *, int); +//int leakage_hasconverged(Project *); + +// Local functions +static int findnodecoeffs(Project* pr, int i, double *hfa, + double *gfa, double *hva, double *gva); +static void evalnodeleakage(double RQtol, double q, double c, + double n, double *h, double *g); +static void addlowerbarrier(double q, double* h, double* g); + +int leakage_open(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: returns an error code +** Purpose: allocates memory for nodal leakage objects +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + hyd->Leakage = (Sleakage *)calloc(net->Njuncs + 1, sizeof(Sleakage)); + if (hyd->Leakage == NULL) return 101; + else return 0; +} + +void leakage_init(Project *pr) +/*------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: initializes the coefficients used by each junction +** node's Leakage object. +**------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + + int i; + double c_area, c_expan, c_orif, len; + Slink *link; + Snode *node1; + Snode *node2; + + // Initialize contents of junction Leakage objects + 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; + } + + // Accumulate pipe leak coeffs. to their respective end nodes + hyd->HasLeakage = FALSE; + 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; + + // 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; + hyd->HasLeakage = TRUE; + + // Adjustment 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; + + // Accumulate coeffs. at 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; + } + } + + // Compute coeffs. for inverted form of the leakage formula + // at each junction node + 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 leakage_close(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 leakage_getlinkleakage(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 leakage_solvercoeffs(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes coeffs. of the linearized hydraulic eqns. +** contributed by pipe leakage. +**-------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + Hydraul *hyd = &pr->hydraul; + Smatrix *sm = &hyd->smatrix; + + int i, row; + double hfa, // fixed area head producing current fixed area leakage + gfa, // gradient of fixed area head w.r.t. leakage flow + hva, // variable area head producing current variable area leakage + gva; // gradient of variable area head w.r.t. leakage flow + + Snode* node; + + for (i = 1; i <= net->Njuncs; i++) + { + // Skip junctions that don't leak + node = &net->Node[i]; + if (!findnodecoeffs(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 leakage_getflowchange(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 (!findnodecoeffs(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 leakage_hasconverged(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns 1 if leakage calculations converged, +** 0 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 0; + } + return 1; +} + +int findnodecoeffs(Project* pr, int i, double *hfa, double *gfa, + double *hva, double *gva) +/* +**-------------------------------------------------------------- +** Input: i = node index +** Output: hfa = fixed area leak head (ft) +** gfa = gradient of fixed area head (ft/cfs) +** hva = variable area leak head (ft) +** gva = gradient of variable area head (ft/cfs) +** returns TRUE if node has leakage, FALSE otherwise +** Purpose: finds head 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 + evalnodeleakage(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 + evalnodeleakage(hyd->RQtol, hyd->Leakage[i].qva, hyd->Leakage[i].cva, + 1.5, hva, gva); + return TRUE; +} + +void evalnodeleakage(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 + addlowerbarrier(q, h, g); +} + +void addlowerbarrier(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/quality.c b/src/quality.c index 41869fe5..0739c15d 100644 --- a/src/quality.c +++ b/src/quality.c @@ -410,6 +410,7 @@ int closequal(Project *pr) FREE(qual->FlowDir); FREE(qual->SortedNodes); } + freeadjlists(&pr->network); return errcode; } diff --git a/src/report.c b/src/report.c index 31e0649e..25bb1376 100644 --- a/src/report.c +++ b/src/report.c @@ -291,7 +291,7 @@ void writesummary(Project *pr) if (qual->Qualflag == NONE || time->Dur == 0.0) sprintf(s, FMT29); else if (qual->Qualflag == CHEM) sprintf(s, FMT30, qual->ChemName); else if (qual->Qualflag == TRACE) sprintf(s, FMT31, net->Node[qual->TraceNode].ID); - else if (qual->Qualflag == AGE) printf(s, FMT32); + else if (qual->Qualflag == AGE) sprintf(s, FMT32); writeline(pr, s); if (qual->Qualflag != NONE && time->Dur > 0) { @@ -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 3f3b50ab..84b05db7 100755 --- a/src/text.h +++ b/src/text.h @@ -1,13 +1,13 @@ /* ****************************************************************************** Project: OWA EPANET - Version: 2.2 + Version: 2.3 Module: text.h Description: string constants used throughout EPANET Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 02/05/2023 + Last Updated: 06/15/2024 ****************************************************************************** */ @@ -130,6 +130,7 @@ #define w_TOLERANCE "TOLER" #define w_EMITTER "EMIT" #define w_BACKFLOW "BACK" +#define w_ALLOWED "ALLOW" #define w_PRICE "PRICE" #define w_DMNDCHARGE "DEMAN" @@ -150,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" @@ -207,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..dc2657cd --- /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/07/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.001); + +// 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()