OpenWaterAnalytics
Calculate the flow from/to a surface model
Some recent discussions on the OpenSWMM mailing list shows that there is some interest in linking SWMM to a 2D surface domain. This as been done by some commercial packages, and by researchers as well [1, 2]. I achieved this in Itzï by using Cython, but it could be much cleaner (and possibly faster) if this computation is done directly in the SWMM model.
I believe this would be a nice addition to SWMM. I have an adea on how it could be implemented. I'll detail it is subsequent posts.
See related issues: #110
[1] Seyoum, Solomon Dagnachew, Zoran Vojinovic, Roland K Price, and Sutat Weesakul. 2012. “Coupled 1D and Noninertia 2D Flood Inundation Model for Simulation of Urban Flooding.” Journal of Hydraulic Engineering 138 (1). American Society of Civil Engineers:23–34. [2] Leandro, J., and R. Martins. 2016. “A Methodology for Linking 2D Overland Flow Models with the Sewer Network Model SWMM 5.1 Based on Dynamic Link Libraries.” Water Science and Technology 73 (12):3017–26. https://doi.org/10.2166/wst.2016.171.
I think that it could be implemented in the following way, by using orifice and weir equations to calculate the flow:
A function that loop through the nodes, calling another function for each node:
- if the node is not linked, continue
- determine the type of equation to use, according to the relative heads in the node and the domain
- calculate and apply the lateral flow
To calculate the flow, the main function needs to be aware of two values from the surface model:
- The surface area on which it applies (i.e, raster cell, computing element, etc.)
- The water depth on the surface (considering that the ground surface elevation is equal to the crest elevation)
This could be done in two way:
- Either the loop function takes as an argument a AoS of the same size as the SWMM Node array, or
- the Node struct contain those values, which are set with another function
I think that the first option my be the best, as the AoS could contain the calculated flow, as a feedback to the surface model.
The TNode struct should have new values added:
- The linkage type (int), from an Enum that could be:
- Not linked
- Linked, no flow
- Orifice
- Free weir
- Submerged weir
- Orifice coefficient
- Free weir coefficient
- Submerged weir coefficient
A more advanced way could be to link to another "manhole cover" object, that would determine those coefficients and the manhole shape. Those values influence the actual flow
My knowledge of the SWMM code base is not very deep, so I'd be glad to have the opinion of other developers on that.
It might be adequate to use ExtInflow for the calculated flow, as it could transfer pollutants. However, could ExtInflow be negative in case of outflow?
Some theory on the use of weir and orifice for the flow calculation are found here: http://documents.irevues.inist.fr/bitstream/handle/2042/25250/0465_239chen.pdf
I would propose to simplify the solution described in this paper by using only the orifice equation when the node is overflowing.
This solution as been proven with comparison to a physical model:
Rubinato, Matteo, Ricardo Martins, Georges Kesserwani, Jorge Leandro, Slobodan Djordjević, and James Shucksmith. 2017. “Experimental Calibration and Validation of Sewer/surface Flow Exchange Equations in Steady and Unsteady Flow Conditions.” Journal of Hydrology 552 (September):421–32. https://doi.org/10.1016/j.jhydrol.2017.06.024.
I would expect that there would be a benefit to having BOTH the Weir and orifice equation applying, as described in many texts such as the 1956 John Hopkins University "The Design of Storm Water Inlets", whereby the opening operates as a weir, until it is drowned and then acts as an orifice. For this to occur an effective weir length L is required and an effective clear opening area A (accounting for bars etc.)
Noting that at times there may be two sets of these required for say a Grated Opening plus a Kerb Inlet.
The Alternative approach is a RATING Curve.
Note also that application of a Blockage factor or functional would also be very useful.
I would expect that there would be a benefit to having BOTH the Weir and orifice equation applying, as described in many texts such as the 1956 John Hopkins University "The Design of Storm Water Inlets", whereby the opening operates as a weir, until it is drowned and then acts as an orifice. For this to occur an effective weir length L is required and an effective clear opening area A (accounting for bars etc.)
What I propose:
- When the surface drains into the sewer, use free weir, submerged weir or orifice, depending on the relative water heads and crest elevation (see the PDF I shared)
- When the drainage is overflowing, orifice.
However, different authors use different combinations, and it does not seems that there is a definite answer to that.
For the weir length and opening area, I was thinking that we could start with a simplified model that:
- consider the opening area equal to the manhole surface area.
- deduce the weir width from the manhole area, assuming it is circular.
However, in reality the lid is usually smaller than the actual manhole, and a more complex system might be useful. It is what I was referring to when proposing the addition of a "manhole cover" object. This object could contain the following information:
- cover surface area
- circumference of the opening
- orifice and weir coefficients
- ...
Noting that at times there may be two sets of these required for say a Grated Opening plus a Kerb Inlet.
Of course, the reality is more complex, and there are way to represent it more finely. For example, see this article about drainage and overflow from various inlets, and taking into account the lifting force of the manhole cover:
Chen, A. S., Leandro, J., & Djordjević, S. (2015). Modelling sewer discharge via displacement of manhole covers during flood events using 1D/2D SIPSON/P-DWave dual drainage simulations. Urban Water Journal, 9006(August), 1–11. https://doi.org/10.1080/1573062X.2015.1041991
However, implementing a solution like this might be much more complex, and we have to evaluate the its actual benefit. For a software like SWMM that is usually used on a large scale (city-wide), I'm not sure if it is worth the pain. In the real world, having complete information of the main drainage network is already quite rare. I don't know if many utilities have informations about the inlets locations and numbers, yet alone their types. There might be some use cases where those modelling details are useful, but for now it looks like a first world problem to me.
Note also that application of a Blockage factor or functional would also be very useful.
That would be indeed interesting. Do you have any literature recommendation on that matter?
In Australia there are a few papers on Culvert Blockage (more so than pits) however many local authorities require a blockage factor to be applied when undertaking analysis. I will try and track down resources. For example Wollongong City Council specifies 50% and 80% blockage in:
http://www.wollongong.nsw.gov.au/council/governance/Policies/Wollongong%20DCP%202009%20Chapter%20E14%20-%20Stormwater%20Management.pdf
Similarly:
http://www.nillumbik.vic.gov.au/files/assets/public/council/council-publications/drainage_design_guidelines_january_2013.pdf
However the notion of using a Factor on the calculated flow (Q) for the pit inlet capacity is useful, for both analysing the impact of blockage, but also testing the system with say a Factor of 2 to double the capacity (to easily assess impact of doubling the pit capacity).
@RudyFrom3 It seems to me that those documents are more like general design guidelines that call for a "safety factor", more than the actual occurring clogging, and even less how the blockage build-up in time during an event. We could potentially add a inlet efficiency coefficient, but it might be redundant with the weir and orifice coefficients. Furthermore, clogging is not the only source of inlet inefficiency [1]. Accurately estimating the actual flow passing through an inlet is a complex subject and it does not seems settled yet among the scientific and engineering community
[1] Inefficiency of storm water inlets as a source of urban floods J. Despotovic, J. Plavsic, N. Stefanovic, D. Pavlovic Water Science and Technology Jan 2005, 51 (2) 139-145;
As for representing multiple inlets, creating multiple smaller junctions (see #112) might be the easiest way to do it, instead of setting complex multiple-linking elements.
I would have thought that the NODE simply be given some surface OPENING characteristics that define how FLOW from the surface enters the NODE. The NODE is simply a small storage that is part of the the internal transfer between links, and allows for additional surface flow to enter through the OPENING that is on the surface. Typically the opening would act as a weir initially and then at some higher depth act as an orifice. (Or it could be described by a rating curve Depth Versus Flow)
@RudyFrom3 I agree with you. That is why I was proposing a a "manhole cover" object. See my post above.
However, in reality the lid is usually smaller than the actual manhole, and a more complex system might be useful. It is what I was referring to when proposing the addition of a "manhole cover" object. This object could contain the following information:
cover surface area circumference of the opening orifice and weir coefficients ...
But this does not negate the need of a way to set the manhole (not opening) area individually.
The only reason I think it would be useful for the Manhole Lid Object to have Multiple Child Objects that represent openings is as in the following example:
Where the Grate has a Weir Length initially (Pink Lines) then an Effective Area (The openings in the grate) and the Kerb Inlet has initially a weir length (green line) followed by a vertical orifice opening being the yellow outline....
So the calculation for discharge would involve:
Qgrate(low) = CLD^1.5
Qkerb(low) =CLD^1.5
Qgrate(high) = Q = Cd A\sqrt{2gD}
Qkerb(high) = Q = Cd A\sqrt{2gD} (Or better for large lateral apertures?)
Etc....
@RudyFrom3 Maybe one Node object could have multiple nodeOpenning objects to take into account those cases?
But then again, this type of approach implies substantial changes in the swmm code.
As a first step, we could consider that the nodes are circular and have one opening of the size of the node surfaceArea. It is far from perfect and it is a strong simplification, but the implementation would be much easier. I think I could implement this rather quickly, as I already did something very similar in another software. This would allow us to test it and work on the interface with our respective 2d models.
Then, later, we could work on adding the nodeOpenning objects.
But, the nodeOpenning object struct would have the advantage of structuring better the data. Without it, we will have orifice/weir coefficients inside the Node struct, which is not ideal.
Another thing to think about is what happen to the overflow and ponding calculations that are already present in swmm.
I'm working on an implementation. It might take some time before it is functional.
You should also worry about what happens when a manhole surcharges as it shifts to a new node solution when the depth in the node is above the ycrown of the highest connecting link to the node. It might be cleaner to just have storage nodes connect to a 2D mesh.
@dickinsonre Thank you for your comment.
I think you refer to what is done in setNodeDepth()?
Why is a surcharged node treated differently? If it is for stability, what would make a storage more stable?
A quick answer and a better answer is to took at the EPA Hydraulics manual for SWMM5, Lew Rossman does a great job of presenting HOW the Surcharge solution is different than the normal node continuity equation. Here is a link as well as on OPENSWMM.ORG https://www.openswmm.org/Topic/4341/surcharged-conduits-and-phantom-storage talks about problems with tunnels and the surcharge algorithm. If i may make a suggestion, before you make changes to the SWMM5 code you should try to understand how it works fully, but to answer your last question. A storage nodes has a solution that is very different in many aspects than a node that is surcharged and a node that is not surcharged. If you would look at the SWMM5 dynamic wave code it will be more clear as to HOW it is different.
Thinking about this more there might be a simple flag that can be added to make the solution for storage or surcharged node transparent to the user. A key engine parameter is called yCrown (the crown of the highest connected pipe to a node), The program shifts to the surcharge algorithm at 0.96 * yCrown. If you change the code so that a 2D node has a yCrown above the rim elevation of the node then the surcharge algothim will not be used for node but the normal node continuity equation with the default surface area.
Regarding the weir vs orifice, i would suggest a weir only for the following reason Note: Weir and Orifice Flow Equations for a Weir in SWMM 5 If you use a weir in SWMM 5 then two flow equations are used
-
The weir uses the weir flow equation when the head at the weir is between the invert elevation of the weir and the crown of the weir and -
An orifice equation when the head is above the weir crown or the weir is submerged
Image is shown below https://swmm5.org/2013/08/04/weir-and-orifice-flow-equations-for-a-weir-in-swmm-5/
@dickinsonre Thank you ery much for your comments.
A storage nodes has a solution that is very different in many aspects than a node that is surcharged and a node that is not surcharged. If you would look at the SWMM5 dynamic wave code it will be more clear as to HOW it is different.
I read the relevant part in the swmm hydraulic manual. I think I understand HOW the codes that update H for surcharged node and storage nodes are different. What I still fail to understand is WHY it is different. In dynamic wave routing, a node will always have a minSurfaceArea and a volume. So it is a kind of simplified storage. AFAIU, a storage node can surcharge (i.e all connected links are full), but will always use the under-relaxation scheme for updating its water depth, never the surcharge scheme. Why all types nodes are not doing the same?
Thinking about this more there might be a simple flag that can be added to make the solution for storage or surcharged node transparent to the user. A key engine parameter is called yCrown (the crown of the highest connected pipe to a node), The program shifts to the surcharge algorithm at 0.96 * yCrown. If you change the code so that a 2D node has a yCrown above the rim elevation of the node then the surcharge algothim will not be used for node but the normal node continuity equation with the default surface area.
That is indeed a neat idea. The drawback is that those nodes will never be reported as surcharged, I guess.
The other option would be to create another class of node, but it will imply large changes in the engine.
Regarding the weir vs orifice, i would suggest a weir only for the following reason Note: Weir and Orifice Flow Equations for a Weir in SWMM 5 If you use a weir in SWMM 5 then two flow equations are used
The weir uses the weir flow equation when the head at the weir is between the invert elevation of the weir and the crown of the weir and An orifice equation when the head is above the weir crown or the weir is submergedImage is shown below https://swmm5.org/2013/08/04/weir-and-orifice-flow-equations-for-a-weir-in-swmm-5/
I was thinking of simply hard-coding the weir and orifice equations in the coupling function. But it could make sense to re-use the swmm weir/orifice equations.
It is always better to use the existing code. You can automatically make a new orifice and/or weir in the engine to communicate flows. Key code from SWMM5 for the handling of the calculation of depth in junctions, storage nodes etc
void setNodeDepth(int i, double dt)
//
// Input: i = node index
// dt = time step (sec)
// Output: none
// Purpose: sets depth at non-outfall node after current time step.
//
{
int canPond; // TRUE if node can pond overflows
int isPonded; // TRUE if node is currently ponded
double dQ; // inflow minus outflow at node (cfs)
double dV; // change in node volume (ft3)
double dy; // change in node depth (ft)
double yMax; // max. depth at node (ft)
double yOld; // node depth at previous time step (ft)
double yLast; // previous node depth (ft)
double yNew; // new node depth (ft)
double yCrown; // depth to node crown (ft)
double surfArea; // node surface area (ft2)
double denom; // denominator term
double corr; // correction factor
double f=0.0; // relative surcharge depth
// --- see if node can pond water above it
canPond = (AllowPonding && Node[i].pondedArea > 0.0);
isPonded = (canPond && Node[i].newDepth > Node[i].fullDepth);
// --- initialize values
yCrown = Node[i].crownElev - Node[i].invertElev;
yOld = Node[i].oldDepth;
yLast = Node[i].newDepth;
Node[i].overflow = 0.0;
surfArea = Xnode[i].newSurfArea;
// --- determine average net flow volume into node over the time step
dQ = Node[i].inflow - Node[i].outflow;
dV = 0.5 * (Node[i].oldNetInflow + dQ) * dt;
**# // --- if node not surcharged or the node is a storage node , base depth change on surface area**
if ( yLast <= yCrown || Node[i].type == STORAGE || isPonded )
{
dy = dV / surfArea;
yNew = yOld + dy;
// --- save non-ponded surface area for use in surcharge algorithm //(5.1.002)
if ( !isPonded ) Xnode[i].oldSurfArea = surfArea; //(5.1.002)
// --- apply under-relaxation to new depth estimate
if ( Steps > 0 )
{
yNew = (1.0 - Omega) * yLast + Omega * yNew;
}
// --- don't allow a ponded node to drop much below full depth
if ( isPonded && yNew < Node[i].fullDepth )
yNew = Node[i].fullDepth - FUDGE;
}
// --- if node surcharged, base depth change on dqdh
// NOTE: depth change is w.r.t depth from previous
// iteration; also, do not apply under-relaxation.
else
{
// --- apply correction factor for upstream terminal nodes
corr = 1.0;
if ( Node[i].degree < 0 ) corr = 0.6;
// --- allow surface area from last non-surcharged condition
// to influence dqdh if depth close to crown depth
denom = Xnode[i].sumdqdh;
if ( yLast < 1.25 * yCrown )
{
f = (yLast - yCrown) / yCrown;
denom += (Xnode[i].oldSurfArea/dt -
Xnode[i].sumdqdh) * exp(-15.0 * f);
}
// --- compute new estimate of node depth
if ( denom == 0.0 ) dy = 0.0;
else dy = corr * dQ / denom;
yNew = yLast + dy;
if ( yNew < yCrown ) yNew = yCrown - FUDGE;
// --- don't allow a newly ponded node to rise much above full depth
if ( canPond && yNew > Node[i].fullDepth )
yNew = Node[i].fullDepth + FUDGE;
}
// --- depth cannot be negative
if ( yNew < 0 ) yNew = 0.0;
// --- determine max. non-flooded depth
yMax = Node[i].fullDepth;
if ( canPond == FALSE ) yMax += Node[i].surDepth;
// --- find flooded depth & volume
if ( yNew > yMax )
{
yNew = getFloodedDepth(i, canPond, dV, yNew, yMax, dt);
}
else Node[i].newVolume = node_getVolume(i, yNew);
// --- compute change in depth w.r.t. time
Xnode[i].dYdT = fabs(yNew - yOld) / dt;
based on #124 being closed off, I am assuming that the weir/orifice connection to SWMM is now operational and that the stumbling block for 2D models now is only on how to link the exchange of information to and from SWMM? So to drive the Weir/Orifice a value of depth is required that is representative of the depth over a pit Top. For Pits that Surcharge, the Flow from SWMM needs to be passed back to the 2D model. There is the third required component.... Once the depth is passed, the Weir/Orifice calculation identifies a Flow, which may or may not be excepted by the connecting pipe link. So a secondary calculation is required in SWMM to identify that the Pipe has the capacity to except the total flow (that coming from pipes upstream if any, and the Pit Top Flow). This calculation may inform that the PitTop flow is in fact limited to something less than that calculated by the PitTop. It is important to do this check as this is the flow that needs to be removed from the 2D model, and also to ensure Volumetric accounting is accurate. That is, no water is lost, because of differences in Pipe & Pit capacities! So to Summarise there are three pieces of information required and a further check:
- Water Depth in the 2D domain
- Calculated PitTop FLOW
- Confirmed Pit Top Flow (Adjusted based on pipe Capacity) (Flow taken from the 2D domain)
- CHECK if Flow Surcharged in the case that the upstream pipe has greater capacity than a downstream pipe and the PitTop inflow is negative... (Surcharging)
@RudyFrom3
Once the depth is passed, the Weir/Orifice calculation identifies a Flow, which may or may not be excepted by the connecting pipe link. So a secondary calculation is required in SWMM to identify that the Pipe has the capacity to except the total flow (that coming from pipes upstream if any, and the Pit Top Flow)
When using the dynamic solver, a SWMM node is a storage. All the flows that enter a node(from pipes, user inflow, catchment, etc.) are added to/subtracted from the amount of water stored in the node. This in turn changes the water surface elevation in the node. The latter is used to calculate the flow in the pipes connected to the node at the next time step. The mass balance is maintained at all time.
Now, the added interchange code can add or remove quite a bit of mass from one time-step to another, which in turn could induce oscillations. Some techniques are implemented in the code to limit the potential instabilities. See especially the coupling_findNodeInflow function. First the exchange flow cannot change sign from one time-step to another. Second, there is a flow limiter that prevents negative depth in the surface model.
In my PhD thesis there are more details on how this was implemented in Itzï, and the influence of those techniques on the coupled model stability. The sections of interest are 2.2 and 3.2.
Hi all, just a desperate Plea for a Status Update ?? What needs to be completed for there to be a Generic available method to link SWMM to 2D models ? I am at a point where I am willing to pay some one if there needs to be coding done ? Where to from here ?
@RudyFrom3,
the feature "as-is" is already sitting in the feature-2dflood feature branch. The biggest thing that needs to happen at this point are unit/reg tests and documentation. I'm not sure of the status of it - but if you are to begin using it, you should log bugs and other issues that you have with it. When a contributor brings this 2D link implementation to the finish line, the we can merge it all the way up to the develop branch. But nothing should be stopping anyone from using it "as-is" today.
Does this make sense?
I'm particularly interesting in @dickinsonre's thoughts on this as well
Ok... is there any where a small example of how it is to operate ? - How to pass Depth from the 2D model to SWMM - How to check the SWMM state (Pipes already full ? or able to accept Water) - Report Flow accepted by SWMM back to 2D model - Pass Flow from 'End of Pipe' to the 2D model
Once I can determine how this is done I'd be very happy to set up various tests....
I will comment on this, Rudy, later in the week. It might take me all werk @RudyFrom3 due to just getting back to my normal work schedule after a long month of traveling
@bemcdonnell I guess the best documentation would be a code example. Maybe I should make a pyswmm PR for the code I have? I really want to implement that new functionality in Itzï, but I did not find the time to do it yet. When I eventually do so, I can write a more generic documentation at the same time.
@RudyFrom3 I try to do a synthesis:
First, in the nodes you want to couple you have to create openings with swmm_setNodeOpening(), and set the couplingArea with swmm_setNodeParam().
As well, you have to make sure that you allow ponding in nodes, otherwise the node will never overflow (that is still something that could be improved, but @dickinsonre had a good idea how to do it.).
Then, at every timestep:
swmm_setNodeParam() to set the new overlandDepth (and couplingArea if changed).
swmm_getNodeResult() to get couplingInflow to remove from the surface model.
That is pretty much it. At every SWMM timestep, SWMM will calculate couplingInflow using overlandDepth and the openings parameters. couplingInflow is counted by SWMM as part of lateral inflows.