search

USEPA / Stormwater-Management-Model

Dynamic hydrology-hydraulic water quality simulation model
223 stars 173 forks source link

Incorrect Geometry for Custom Elliptical Pipes #144

Open LRossman opened 10 months ago

LRossman commented 10 months ago

As noted recently on the Open SWMM Knowledge Base, SWMM produces inconsistent results for custom sized elliptical pipes. The problem stems from it using incorrect equations for the pipe's full flow area and hydraulic radius that came from SWMM 4.4. A corrected set of equations, derived from Table 4 of the Concrete Pipe Design Manual (American Concrete Pipe Association, 2011) are:

A = 0.8117 * (Height * Width)
Rh = 0.2448 * sqrt(Height * Width)

with all values in feet.

These equations should replace those used in lines 572-573 and 600-601 in file xsect.c, i.e.,

    xsect->aFull = 0.8117 * xsect->yFull * xsect->wMax;
    xsect->rFull = 0.2448 * sqrt(xsect->yFull * xsect->wMax);
dickinsonre commented 10 months ago

The mentioned code from the 2004 version of SWMM4 (OSU version)

C======================================================================= CIM START <><><><><><><><><> CIM THE FOLLOWING HOLD TRUE FOR INPUT ONLY, ARE MODIFIED CIM TO INTERNAL NUMBERING SCHEME AT END OF INDAT1 SUBROUTINE CIM Note: In revised internal numbering scheme, cim o Conduit types 1 through 20 are reserved cim for closed conduit types that use TWNORM, ANORM, and cim HRNORM ARRAYS cim o Conduit types 21 through 50 are special types including cim rectangular pipes (22), and also includes parabolic and cim irregular sections (24) that use QCURVE arrays. cimbridge Now includes bridges type 25 cim o Conduit types 51 through 60 are reserved for equivalent cim orifice types. CIM CIM NKLASS TYPE NEWKLASS C 1 CIRCULAR PIPE NEWKLASS(1)=1 C 2 RECTANGULAR PIPE NEWKLASS(2)=21 C 3 HORSESHOE PIPE NEWKLASS(3)=2 C 4 EGGSHAPED PIPE NEWKLASS(4)=3 C 5 BASKETHANDLE PIPE NEWKLASS(5)=4 C 6 TRAPEZOIDAL CHANNEL NEWKLASS(6)=22 C 7 PARABOLIC CROSS-SECTION OR POWER C FUNCTION CROSS SECTION NEWKLASS(7)=23 C 8 IRREGULAR CROSS-SECTION NEWKLASS(8)=24 C 9 HORIZONTAL ELLIPSE C (LONG AXIS IS HORIZONTAL) NEWKLASS(9)=5 C 10 VERTICAL ELLIPSE NEWKLASS(10)=6 C 11 ARCH NEWKLASS(11)=7 C 12 BRIDGE (TO BE IMPLEMENTED) NEWKLASS(12)=25 CIM OOOOOOOOO CIM THESE ARE OLD ASSIGNMENTS AND NEW ASSIGNMENTS FOR ORIFICE CIM EQUIVALENT PIPES C 9 SIDE OUTLET CIRCULAR ORIFICE 51 C 10 BOTTOM OUTLET CIRCULAR ORIFICE 52 C 11 SIDE OUTLET RECTANGULAR ORIFICE 53 C 12 BOTTOM OUTLET RECTANGULAR ORIFICE 54 CIM OOOOOOOOOO CIM END <><><><><><><><><><><><> C======================================================================= KLASS = NKLASS(N) CIM START <><><><><><><><> CIM IF(KLASS.EQ.1.OR.KLASS.GE.9) THEN CIM 9 WAS ORIFICE, DON'T KNOW DATA YET c change this to select case SELECT CASE (KLASS) CASE (1) CIM END <><><><><><> RFULL(N) = DEEP(N)/4.0 AFULL(N) = (3.1415926/4.0)*DEEP(N)*2 WIDE(N) = DEEP(N) cim CALL SEDEPTH1 to adjust pipe if IPIPESED=1 IF (SEDEPTH(N).NE.0.0) THEN CALL SEDEPTH1(N) ELSE SEDAREA(N) = 0.0 SEDPERI(N) = 0.0 SEDRAD(N) = 0.0 ENDIF CASE (2) RFULL(N) = (WIDE(N)DEEP(N))/(2.WIDE(N)+2.0DEEP(N)) AFULL(N) = WIDE(N)DEEP(N) CASE (3) RFULL(N) = 0.25381 DEEP(N) CASE (4) RFULL(N) = 0.19311 DEEP(N) CASE (5) RFULL(N) = 0.28800DEEP(N) CASE (6) AFULL(N)=DEEP(N)(WIDE(N)+DEEP(N)/2.(STHETA(N)+SPHI(N))) IF(WIDE(N).LE.0.0) WIDE(N) = 0.01 RFULL(N)=AFULL(N)/(WIDE(N)+DEEP(N)*

dickinsonre commented 10 months ago

Always hard to read the SWMM4 code but here ae the vertical ellipse only arfull and rfull

CIM VERTICAL ELLIPSE CASE (10) Cwch (CIM), 7/20/04. Should be .LT. for both. C IF ((DEEP(N).GT.0.0).AND.(WIDE(N).LT.0.0)) THEN IF ((DEEP(N).LT.0.0).AND.(WIDE(N).LT.0.0)) THEN DEEP(N) = - DEEP(N) WIDE(N) = - WIDE(N) AFULL(N) = 1.2692 WIDE(N) WIDE(N) RFULL(N) = 0.3061 * WIDE(N)

LRossman commented 10 months ago

I'm now having second thoughts about applying the aforementioned area and hydraulic radius formulas for custom elliptical pipes to any pair of of width (span) and height (rise) values. The formulas were derived from the following table for standard sizes of elliptical pipes (with rise and span in inches, area A in sq. ft, and hydraulic radius Rh in ft.):

Rise    Span    S/R A   Rh
14  23  1.64    1.8 0.37
19  30  1.58    3.3 0.49
22  34  1.55    4.1 0.55
24  38  1.58    5.1 0.61
27  42  1.56    6.3 0.69
29  45  1.55    7.4 0.74
32  49  1.53    8.8 0.81
34  53  1.56    10.2    0.88
38  60  1.58    12.9    0.97
43  68  1.58    16.6    1.11
48  76  1.58    20.5    1.23
53  83  1.57    24.8    1.35
58  91  1.57    29.5    1.48
63  98  1.56    34.6    1.60
68  106 1.56    40.1    1.72
72  113 1.57    46.1    1.85
77  121 1.57    52.4    1.97
82  128 1.56    59.2    2.09
87  136 1.56    66.4    2.22
92  143 1.55    74  2.34
97  151 1.56    82  2.46
106 166 1.57    99.2    2.71
116 180 1.55    118.6   2.97

Note that the ratios of span to rise are mostly the same, averaging 1.56. Thus the area and hydraulic radius formulas fitted to the product of span and rise from this table only apply to ratios of around 1.56 and not to any arbitrary pair of values. What this implies is that the dimensions that SWMM should accept for a custom elliptical pipe must comply with the 1.56 ratio. The user should only be required to supply a single value for the rise (full depth) with the span being set to 1.56 that value for horizontal elliptical pipes and 1/1.56 for vertical elliptical pipes.

cbuahin commented 9 months ago

The formula has been updated. However, longer term, we will update the GUI to either accept only rise or span and do the multiplication with 1.56 in the code.

cbuahin commented 6 months ago

The ratio is not always 1.56. I think it is best to leave it up to users to ensure they are using the correct rise and spans.

cbuahin commented 6 months ago

@LRossman I presume you reopened this because I reopened the corresponding issue on the GUI side, correct? I will post the message I posted there here as well.

I set up a quick experiment using a 100ft pipe with 10ft sub-sections and slopes of 1% and 5% to see the impact of using the constant S/R ratio of 1.56 versus the correct ratio on depth and velocities given the same flows. The results showed some differences. There were larger deviations for the smaller pipes than for the larger pipes as one would expect.

image

Also, the deviations tended to be larger when depths approached the crowns of the pipes, which is also to be expected. Finally, the vertical ellipses showed larger deviations in depth than the horizontal ellipses as changes in flow lead to relatively larger changes in depth.

image

In light of these differences, it makes sense to me to keep the status-quo. Users have the option of selecting one of the standard sizes or they could prescribe their own rise and span values which will lead to more accurate and precise calculations.

LRossman commented 6 months ago

I wanted to reopen this issue because I still believe that the equations SWMM uses for full area and hydraulic radius in custom elliptical pipes only apply to dimensions where the span (S) over the rise (R) for a horizontal ellipse is 1.56 (or more precisely 1.56363). It would be the reciprocal of that for a vertical ellipse. The original equations appearing in SWMM 4.4A were

A = 1.2692 * R * R 
Rh = 0.3061 * R

They were documented in a file named SHAPE.DOC.TXT (last updated on 7/20/2004) that was part of the 4.4A 1999 release. They would make no sense if applied to any arbitrary set of span and rise dimensions since they totally ignore the span dimension. However look what happens if we substitute R = S / (S/R) = S / 1.56363 into them:

A = 0.8117 * R * S / (S/R) = 0.8117 / 1.56363 * R * S = 0.8117 * R * S
Rh = 0.3061 * sqrt(R * S / (S/R)) = 0.3061 / sqrt(1.56363) * sqrt(R * S) = 0.2448 * sqrt(R * S)

The coefficient 0.8117 is precisely the result obtained when the area in the table of standard elliptical pipes (shown several comments above or in Table D-1 of the SWMM 5 Hydraulics Reference Manual) is fitted to the R * S entries as shown below: image A similar result holds for the 0.2448 coefficient when the tabulated hydraulic radius is fitted to the square root of R * S. Thus the original 4.4A equations are correct only for S/R = 1.56363 as are the updated ones recently added to the bug_fixes branch. Since SWMM currently accepts any pair of R and S values for a custom elliptical pipe there will be some error introduced if a different S/R ratio applies. The confusion on this issue could have been avoided when writing SWMM 5, and input for a custom shape be restricted to just its rise dimension, if the SHAPE.DOC.TXT note had explicitly stated that its formulas for custom elliptical pipes only applied to a particular S/R.

cbuahin commented 6 months ago

That makes a lot of sense to me now. So, do you have any reservations about me modifying the code (including the GUI) to enforce this ratio? On the GUI end, I could modify it to only accept the rise value and apply the ratio internally to estimate the span.

LRossman commented 6 months ago

There's the issue of backward compatibility to consider. People, especially contractors, get bent out of shape when a new release produces a different set of results than an earlier release did on their project files. So I can see several options to consider:

  1. As you proposed, use just the rise (max. depth) value for the ellipse and ignore the span if present in an input file, plus remove the corresponding Maximum Width input field in the GUI Cross-Section Editor dialog.
  2. As in option 1 but also issue a warning message when reading an input file (in both the engine and GUI) if a span value is found that doesn't satisfy the "std. ratio" when parsing the elliptical shape entry.
  3. Accept a span entry in the input file and leave the Max. Width field in the GUI, but allow the span to be 0 indicating that the "std. ratio" will be used for it. Also add a note to the Cross-Section Editor explaining this as well as the warning message mentioned in option 2 in case a non-zero span is entered.
  4. Leave everything as is but issue the warning message if span/rise doesn't match the "std. ratio".

I'm thinking that maybe option 3 might be best since it can use the correct S/R ratio but also satisfy backward compatibility with a warning that results might be inaccurate. You might want to discuss this issue on the SWMM Users Knowledge Base and ask for suggestions on how to resolve it. Maybe they will feel OK with option 1 and not care about backward compatibility.

Also note that depending on which option is selected, the Help file would need to be updated in the short run as well as the User Manual and Hydraulics Manual in the long run.

cbuahin commented 6 months ago

Number 3 makes a lot of sense to me as well. As you know, updating the user manual requires an internal review that may prolong the the next release we are planning, but I think it is well worth the effort as this bug is an important one to address in the next release. I appreciate your inputs on this.

dickinsonre commented 6 months ago

I would vote for option 4 if the user has an option to make a custom shape where the S/R can be different than 1.56. User would then use a custom shape instead of the custom option in the elliptical shape dialog.

LRossman commented 6 months ago

I think I found a better way to handle the custom elliptical cross section. That is to treat it as a true ellipse rather than use the area and hydraulic radius formulas derived from the table of standard elliptical pipes. As noted above, those formulas are only valid for S/R = 1.56 and don't represent a true ellipse (the full area formula derived from the table is A = 0.8117 * R * S while for a true ellipse it is A = 0.7854 * R * S (i.e. (PI/4)*R*S)).

If the custom elliptical shape is treated as a true ellipse then the full area (A) and hydraulic radius (Rh) formulas would be:

a = R/2
b = S/2
c = PI *(3.0*(a+b) - sqrt((3.0*a + b)*(a + 3.0*b)))  // Circumference
A = PI * a * b
Rh = A / c

After replacing the current formulas used in v5.2.4 with these I was able to eliminate the inconsistencies found by Bruno Bermedo as described in his Open SWMM Knowledge Base posting.. You can find the details in the attached zip file.

So I'm now thinking that the preferred option for handling a custom elliptical shape is:

This will not necessarily preserve backward compatibility, but since doing so would perpetuate incorrect results I think it is a reasonable course of action to take.

P.S. Since the table for standard elliptical pipes applies to shapes with S/R of 1.56 and the areas in the table don't quite conform to a true ellipse the question arises as to what shape was used to derive the geometry tables for A/Afull, Rh/Rhfull, and W/Wmax for equally spaced Y/Yfull values. By looking at the W/Wmax table (where Wmax is S) I've confirmed that it is based on a true ellipse and not whatever shape is implied by the standard pipe table. I did this by using the following formula for the chord length W of a horizontal ellipse:

W = S * sqrt(1 - (2x/R)^2)

where x is the distance between the major axis and the chord. Adapting this formula to match the contents of a W/Wmax table we get:

W/Wmax = sqrt(1 - (1 - (2*Y/Yfull)^2)

When evaluated at 25 equally spaced Y/Yfull values it produces the same numbers seen in the W_HorizEllipse array appearing in the SWMM's xsect.dat file. We can thus assume that the area and hydraulic radius tables were also derived from a true ellipse shape.

Elliptical Pipe Comparisons.zip

cbuahin commented 6 months ago

I hear you @LRossman. Definitely something to consider. However, if we are going to use actual formulas, my thinking is that this will be a good time to use formulas for all the cross-sectional shapes we have formulas for instead of the lookup tables. I can envision some performance gains from direct formulas instead of the lookups.

cbuahin commented 6 months ago

Perhaps we can do this in two steps. The first step is to just implement the warning in the next patch and implement the more comprehensive direct formula approaches in a future minor release.

LRossman commented 6 months ago

@cbuahin please note that the formulas we are discussing here only apply to determining the area and hydraulic radius of a pipe flowing full. SWMM already uses analytic formulas for these quantities. For partly full pipes SWMM uses lookup tables only for circular, elliptical and arch pipes as well as the masonry shapes. Analytical formulas don't exist for the masonry shapes. The formulas for circular pipes involve trigonometric functions which are more expensive to compute than using a lookup table. For elliptical and arch shapes there are no simple formulas for partly full area and wetted perimeter.

Also please see the updated version of my previous posting where I confirm that the geometry tables in xsect.dat for elliptical pipes are based on a true ellipse rather than whatever kind of modified elliptical shape is implied by the table in the Concrete Pipe Design Manual.

As for just issuing a warning for custom elliptical shapes that don't conform to the standard S/R ratio, it would basically have to admit that the results one obtains are incorrect -- that's not a great look. I would rather see the bug corrected.

cbuahin commented 6 months ago

@LRossman, I hear what you are saying. Let me confirm that the analytical formula does not render the A/Afull, Rh/Rhfull, and W/Wmax for equally spaced Y/Yfull values in the lookup tables incorrect on my end (not to dispute your findings but for my own sanity ☺ ) and then I will work on addressing the bug shortly.

cbuahin commented 6 months ago

For the masonry shapes, I know there are no analytical formulas but I have been experimenting with replacing the tables with with polynomial functions with some degree of success and performance improvements. Would you have any reservations about those?

LRossman commented 6 months ago

If you can show performance improvements with the polynomial functions then I say go for it.

Regarding the elliptical lookup tables, I only checked the W/Wmax table using the analytical formula for chord width for a true ellipse and got perfect agreement with it. I don't think there are similarly simple formulas for partly filled area and hydraulic radius. In particular, I think the wetted perimeter (needed for hyd. radius) requires numerical integration (I have no idea how the SWMM hyd. radius table was constructed). If you look at the area table you will see that the increments are the same both below and above the Y/Yfull = 0.5 level which leads me to believe it applies to a true ellipse.

dickinsonre commented 6 months ago

(I have no idea how the SWMM hyd. radius table was constructed) - ,it was done in the 1990's by CDM or CDM Smith. I cannot find shapes.doc but it was added on 1/96

CIM THESE NEXT ARRAYS CONTAIN DATA FOR STANDARD ELLIPTICAL AND CIM ARCH PIPE TYPES: MAJOR AXIS LENGTH, MINOR AXIS LENGTH, CIM FULL FLOW AREA, FULL FLOW HYDRAULIC RADIUS. CIM INITIALIZED IN EXTBLK.FOR CIM SEE SHAPES.DOC FOR STANDARD ELLIPTICAL AND ARCH DIMENSIONS PARAMETER(NUMELL=23) COMMON /ELLIPSE/ EMAJOR(NUMELL),EMINOR(NUMELL),EAREA(NUMELL), .ERADIUS(NUMELL) PARAMETER(NUMARCH=102) COMMON /ARCH/ AMAJOR(NUMARCH),AMINOR(NUMARCH),AAREA(NUMARCH), .ARADIUS(NUMARCH)

dickinsonre commented 6 months ago

The SWMM4 code gives me a headache. Thanks so much for bringing it to C and writing beautiful code. Lew.

dickinsonre commented 6 months ago

Mitch gave me the shapes file - maybe it will help you SHAPE.DOC.TXT