GCC Code Coverage Report


Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 92.3% 60 / 0 / 65
Functions: 100.0% 6 / 0 / 6
Branches: 66.7% 12 / 0 / 18

culvert.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // culvert.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 11/01/21 (Build 5.2.0)
7 // Author: L. Rossman
8 //
9 // Culvert equations for SWMM5
10 //
11 // Computes flow reduction in a culvert-type conduit due to
12 // inlet control using equations from the FHWA HEC-5 circular.
13 //
14 // Update History
15 // ==============
16 // Build 5.1.013:
17 // - C parameter corrected for Arch, Corrugated Metal, Mitered culvert.
18 //-----------------------------------------------------------------------------
19 #define _CRT_SECURE_NO_DEPRECATE
20
21 #include <math.h>
22 #include "findroot.h"
23 #include "headers.h"
24
25 //-----------------------------------------------------------------------------
26 // Constants
27 //-----------------------------------------------------------------------------
28 enum CulvertParam {FORM, K, M, C, Y};
29 static const int MAX_CULVERT_CODE = 57;
30 static const double Params[58][5] = {
31
32 // FORM K M C Y
33 //------------------------------------
34 {0.0, 0.0, 0.0, 0.0, 0.00},
35
36 //Circular concrete
37 {1.0, 0.0098, 2.00, 0.0398, 0.67}, //Square edge w/headwall
38 {1.0, 0.0018, 2.00, 0.0292, 0.74}, //Groove end w/headwall
39 {1.0, 0.0045, 2.00, 0.0317, 0.69}, //Groove end projecting
40
41 //Circular Corrugated Metal Pipe
42 {1.0, 0.0078, 2.00, 0.0379, 0.69}, //Headwall
43 {1.0, 0.0210, 1.33, 0.0463, 0.75}, //Mitered to slope
44 {1.0, 0.0340, 1.50, 0.0553, 0.54}, //Projecting
45
46 //Circular Pipe, Beveled Ring Entrance
47 {1.0, 0.0018, 2.50, 0.0300, 0.74}, //Beveled ring, 45 deg bevels
48 {1.0, 0.0018, 2.50, 0.0243, 0.83}, //Beveled ring, 33.7 deg bevels
49
50 //Rectangular Box with Flared Wingwalls
51 {1.0, 0.026, 1.0, 0.0347, 0.81}, //30-75 deg. wingwall flares
52 {1.0, 0.061, 0.75, 0.0400, 0.80}, //90 or 15 deg. wingwall flares
53 {1.0, 0.061, 0.75, 0.0423, 0.82}, //0 deg. wingwall flares (striaght sides)
54
55 //Rectanglar Box with Flared Wingwalls & Top Edge Bevel
56 {2.0, 0.510, 0.667, 0.0309, 0.80}, //45 deg. flare; 0.43D top edge bevel
57 {2.0, 0.486, 0.667, 0.0249, 0.83}, //18-33.7 deg flare; 0.083D top edge bevel
58
59 //Rectangular Box; 90-deg Headwall; Chamfered or Beveled Inlet Edges
60 {2.0, 0.515, 0.667, 0.0375, 0.79}, //chamfered 3/4-in
61 {2.0, 0.495, 0.667, 0.0314, 0.82}, //beveled 1/2-in/ft at 45 deg (1:1)
62 {2.0, 0.486, 0.667, 0.0252, 0.865}, //beveled 1-in/ft at 33.7 deg (1:1.5)
63
64 //Rectangular Box; Skewed Headwall; Chamfered or Beveled Inlet Edges
65 {2.0, 0.545, 0.667, 0.04505,0.73}, //3/4" chamfered edge, 45 deg skewed headwall
66 {2.0, 0.533, 0.667, 0.0425, 0.705}, //3/4" chamfered edge, 30 deg skewed headwall
67 {2.0, 0.522, 0.667, 0.0402, 0.68}, //3/4" chamfered edge, 15 deg skewed headwall
68 {2.0, 0.498, 0.667, 0.0327, 0.75}, //45 deg beveled edge, 10-45 deg skewed headwall
69
70 //Rectangular box, Non-offset Flared Wingwalls; 3/4" Chamfer at Top of Inlet
71 {2.0, 0.497, 0.667, 0.0339, 0.803}, //45 deg (1:1) wingwall flare
72 {2.0, 0.493, 0.667, 0.0361, 0.806}, //18.4 deg (3:1) wingwall flare
73 {2.0, 0.495, 0.667, 0.0386, 0.71}, //18.4 deg (3:1) wingwall flare, 30 deg inlet skew
74
75 //Rectangular box, Offset Flared Wingwalls, Beveled Edge at Inlet Top
76 {2.0, 0.497, 0.667, 0.0302, 0.835}, //45 deg (1:1) flare, 0.042D top edge bevel
77 {2.0, 0.495, 0.667, 0.0252, 0.881}, //33.7 deg (1.5:1) flare, 0.083D top edge bevel
78 {2.0, 0.493, 0.667, 0.0227, 0.887}, //18.4 deg (3:1) flare, 0.083D top edge bevel
79
80 // Corrugated Metal Box
81 {1.0, 0.0083, 2.00, 0.0379, 0.69}, //90 deg headwall
82 {1.0, 0.0145, 1.75, 0.0419, 0.64}, //Thick wall projecting
83 {1.0, 0.0340, 1.50, 0.0496, 0.57}, //Thin wall projecting
84
85 // Horizontal Ellipse Concrete
86 {1.0, 0.0100, 2.00, 0.0398, 0.67}, //Square edge w/headwall
87 {1.0, 0.0018, 2.50, 0.0292, 0.74}, //Grooved end w/headwall
88 {1.0, 0.0045, 2.00, 0.0317, 0.69}, //Grooved end projecting
89
90 // Vertical Ellipse Concrete
91 {1.0, 0.0100, 2.00, 0.0398, 0.67}, //Square edge w/headwall
92 {1.0, 0.0018, 2.50, 0.0292, 0.74}, //Grooved end w/headwall
93 {1.0, 0.0095, 2.00, 0.0317, 0.69}, //Grooved end projecting
94
95 // Pipe Arch, 18" Corner Radius, Corrugated Metal
96 {1.0, 0.0083, 2.00, 0.0379, 0.69}, //90 deg headwall
97 {1.0, 0.0300, 1.00, 0.0463, 0.75}, //Mitered to slope
98 {1.0, 0.0340, 1.50, 0.0496, 0.57}, //Projecting
99
100 // Pipe Arch, 18" Corner Radius, Corrugated Metal
101 {1.0, 0.0300, 1.50, 0.0496, 0.57}, //Projecting
102 {1.0, 0.0088, 2.00, 0.0368, 0.68}, //No bevels
103 {1.0, 0.0030, 2.00, 0.0269, 0.77}, //33.7 deg bevels
104
105 // Pipe Arch, 31" Corner Radius, Corrugated Metal
106 {1.0, 0.0300, 1.50, 0.0496, 0.57}, //Projecting
107 {1.0, 0.0088, 2.00, 0.0368, 0.68}, //No bevels
108 {1.0, 0.0030, 2.00, 0.0269, 0.77}, //33.7 deg. bevels
109
110 // Arch, Corrugated Metal
111 {1.0, 0.0083, 2.00, 0.0379, 0.69}, //90 deg headwall
112 {1.0, 0.0300, 1.00, 0.0473, 0.75}, //Mitered to slope
113 {1.0, 0.0340, 1.50, 0.0496, 0.57}, //Thin wall projecting
114
115 // Circular Culvert
116 {2.0, 0.534, 0.555, 0.0196, 0.90}, //Smooth tapered inlet throat
117 {2.0, 0.519, 0.640, 0.0210, 0.90}, //Rough tapered inlet throat
118
119 // Elliptical Inlet Face
120 {2.0, 0.536, 0.622, 0.0368, 0.83}, //Tapered inlet, beveled edges
121 {2.0, 0.5035,0.719, 0.0478, 0.80}, //Tapered inlet, square edges
122 {2.0, 0.547, 0.800, 0.0598, 0.75}, //Tapered inlet, thin edge projecting
123
124 // Rectangular
125 {2.0, 0.475, 0.667, 0.0179, 0.97}, //Tapered inlet throat
126
127 // Rectangular Concrete
128 {2.0, 0.560, 0.667, 0.0446, 0.85}, //Side tapered, less favorable edges
129 {2.0, 0.560, 0.667, 0.0378, 0.87}, //Side tapered, more favorable edges
130
131 // Rectangular Concrete
132 {2.0, 0.500, 0.667, 0.0446, 0.65}, //Slope tapered, less favorable edges
133 {2.0, 0.500, 0.667, 0.0378, 0.71} //Slope tapered, more favorable edges
134
135 };
136
137 //-----------------------------------------------------------------------------
138 // Culvert data structure
139 //-----------------------------------------------------------------------------
140 typedef struct
141 {
142 double yFull; // full depth of culvert (ft)
143 double scf; // slope correction factor
144 double dQdH; // Derivative of flow w.r.t. head
145 double qc; // Unsubmerged critical flow
146 double kk;
147 double mm; // Coeffs. for unsubmerged flow
148 double ad;
149 double hPlus; // Intermediate terms
150 TXsect* xsect; // Pointer to culvert cross section
151 } TCulvert;
152
153 //-----------------------------------------------------------------------------
154 // External functions (declared in funcs.h)
155 //-----------------------------------------------------------------------------
156 // double culvert_getInflow
157
158 //-----------------------------------------------------------------------------
159 // Local functions
160 //-----------------------------------------------------------------------------
161 static double getUnsubmergedFlow(int code, double h, TCulvert* culvert);
162 static double getSubmergedFlow(int code, double h, TCulvert* culvert);
163 static double getTransitionFlow(int code, double h, double h1, double h2,
164 TCulvert* culvert);
165 static double getForm1Flow(double h, TCulvert* culvert);
166 static double form1Eqn(double yc, void* p);
167 /*
168 static void report_CulvertControl(int j, double q0, double q, int condition,
169 double yRatio); //for debugging only
170 */
171
172 //=============================================================================
173
174 34053 double culvert_getInflow(int j, double q0, double h)
175 //
176 // Input: j = link index
177 // q0 = unmodified flow rate (cfs)
178 // h = upstream head (ft)
179 // Output: returns modified flow rate through culvert (cfs)
180 // Purpose: uses FHWA HEC-5 equations to find flow through inlet
181 // controlled culverts
182 //
183 {
184 int code, //culvert type code number
185 k, //conduit index
186 condition; //flow condition
187 double y, //current depth (ft)
188 y1, //unsubmerged depth limit (ft)
189 y2, //submerged depth limit (ft)
190 q; //inlet-controlled flow (cfs)
191 TCulvert culvert; //intermediate results
192
193 // --- check that we have a culvert conduit
194
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 34053 times.
34053 if ( Link[j].type != CONDUIT ) return q0;
195 34053 culvert.xsect = &Link[j].xsect;
196 34053 code = culvert.xsect->culvertCode;
197
2/4
✓ Branch 0 taken 34053 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 34053 times.
34053 if ( code <= 0 || code > MAX_CULVERT_CODE ) return q0;
198
199 // --- compute often-used variables
200 34053 k = Link[j].subIndex;
201 34053 culvert.yFull = culvert.xsect->yFull;
202 34053 culvert.ad = culvert.xsect->aFull * sqrt(culvert.yFull);
203
204 // --- slope correction factor (-7 for mitered inlets, 0.5 for others)
205
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 34053 times.
34053 switch (code)
206 {
207 case 5:
208 case 37:
209 case 46: culvert.scf = -7.0 * Conduit[k].slope; break;
210 34053 default: culvert.scf = 0.5 * Conduit[k].slope;
211 }
212
213 // --- find head relative to culvert's upstream invert
214 // (can be greater than yFull when inlet is submerged)
215 34053 y = h - (Node[Link[j].node1].invertElev + Link[j].offset1);
216
217 // --- check for submerged flow (based on FHWA criteria of Q/AD > 4)
218 34053 y2 = culvert.yFull * (16.0 * Params[code][C] + Params[code][Y] - culvert.scf);
219
2/2
✓ Branch 0 taken 6436 times.
✓ Branch 1 taken 27617 times.
34053 if ( y >= y2 )
220 {
221 6436 q = getSubmergedFlow(code, y, &culvert);
222 6436 condition = 2;
223 }
224 else
225 {
226 // --- check for unsubmerged flow (based on arbitrary limit of 0.95 full)
227 27617 y1 = 0.95 * culvert.yFull;
228
2/2
✓ Branch 0 taken 24347 times.
✓ Branch 1 taken 3270 times.
27617 if ( y <= y1 )
229 {
230 24347 q = getUnsubmergedFlow(code, y, &culvert);
231 24347 condition = 1;
232 }
233 // --- flow is in transition zone
234 else
235 {
236 3270 q = getTransitionFlow(code, y, y1, y2, &culvert);
237 3270 condition = 0;
238 }
239 }
240
241 // --- check if inlet controls and replace conduit's value of dq/dh
242
2/2
✓ Branch 0 taken 33968 times.
✓ Branch 1 taken 85 times.
34053 if ( q < q0 )
243 {
244 // --- for debugging only
245 //if ( RptFlags.controls ) report_CulvertControl(j, q0, q, condition,
246 // y / culvert.yFull);
247
248 33968 Link[j].inletControl = TRUE;
249 33968 Link[j].dqdh = culvert.dQdH;
250 33968 return q;
251 }
252 85 else return q0;
253 }
254
255 //=============================================================================
256
257 27617 double getUnsubmergedFlow(int code, double h, TCulvert* culvert)
258 //
259 // Input: code = culvert type code number
260 // h = inlet water depth above culvert invert
261 // culvert = pointer to a culvert data structure
262 // Output: returns flow rate;
263 // computes value of variable Dqdh
264 // Purpose: computes flow rate and its derivative for unsubmerged
265 // culvert inlet.
266 //
267 {
268 double arg;
269 double q;
270
271 // --- assign shared variables
272 27617 culvert->kk = Params[code][K];
273 27617 culvert->mm = Params[code][M];
274 27617 arg = h / culvert->yFull / culvert->kk;
275
276 // --- evaluate correct equation form
277
1/2
✓ Branch 0 taken 27617 times.
✗ Branch 1 not taken.
27617 if ( Params[code][FORM] == 1.0)
278 {
279 27617 q = getForm1Flow(h, culvert);
280 }
281 else q = culvert->ad * pow(arg, 1.0/culvert->mm);
282 27617 culvert->dQdH = q / h / culvert->mm;
283 27617 return q;
284 }
285
286 //=============================================================================
287
288 9706 double getSubmergedFlow(int code, double h, TCulvert* culvert)
289 //
290 // Input: code = culvert type code number
291 // h = inlet head (ft)
292 // culvert = pointer to a culvert data structure
293 // Output: returns flow rate;
294 // computes value of Dqdh
295 // Purpose: computes flow rate and its derivative for submerged
296 // culvert inlet.
297 //
298 {
299 9706 double cc = Params[code][C];
300 9706 double yy = Params[code][Y];
301 9706 double arg = (h/culvert->yFull - yy + culvert->scf) / cc ;
302 double q;
303
304
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9706 times.
9706 if ( arg <= 0.0 )
305 {
306 culvert->dQdH = 0.0;
307 return BIG;
308 }
309 9706 q = sqrt(arg) * culvert->ad;
310 9706 culvert->dQdH = 0.5 * q / arg / culvert->yFull / cc;
311 9706 return q;
312 }
313
314 //=============================================================================
315
316 3270 double getTransitionFlow(int code, double h, double h1, double h2, TCulvert* culvert)
317 //
318 // Input: code = culvert type code number
319 // h = inlet water depth above culvert invert (ft)
320 // h1 = head limit for unsubmerged condition (ft)
321 // h2 = head limit for submerged condition (ft)
322 // culvert = pointer to a culvert data structure
323 // Output: returns flow rate )cfs);
324 // computes value of Dqdh (cfs/ft)
325 // Purpose: computes flow rate and its derivative for inlet-controlled flow
326 // when inlet water depth lies in the transition range between
327 // submerged and unsubmerged conditions.
328 //
329 {
330 3270 double q1 = getUnsubmergedFlow(code, h1, culvert);
331 3270 double q2 = getSubmergedFlow(code, h2, culvert);
332 3270 double q = q1 + (q2 - q1) * (h - h1) / (h2 - h1);
333 3270 culvert->dQdH = (q2 - q1) / (h2 - h1);
334 3270 return q;
335 }
336
337 //=============================================================================
338
339 27617 double getForm1Flow(double h, TCulvert* culvert)
340 //
341 // Input: h = inlet water depth above culvert invert
342 // culvert = pointer to a culvert data structure
343 // Output: returns inlet controlled flow rate
344 // Purpose: computes inlet-controlled flow rate for unsubmerged culvert
345 // using FHWA Equation Form1.
346 //
347 // See pages 195-196 of FHWA HEC-5 (2001) for details.
348 //
349 {
350 // --- save re-used terms in culvert structure
351 27617 culvert->hPlus = h / culvert->yFull + culvert->scf;
352
353 // --- use Ridder's method to solve Equation Form 1 for critical depth
354 // between a range of 0.01h and h
355 27617 findroot_Ridder(0.01*h, h, 0.001, form1Eqn, culvert);
356
357 // --- return the flow value used in evaluating Equation Form 1
358 27617 return culvert->qc;
359 }
360
361 //=============================================================================
362
363 184453 double form1Eqn(double yc, void* p)
364 //
365 // Input: yc = critical depth
366 // p = pointer to a TCulvert object
367 // Output: returns residual error
368 // Purpose: evaluates the error in satisfying FHWA culvert Equation Form1:
369 //
370 // h/yFull + 0.5*s = yc/yFull + yh/2/yFull + K[ac/aFull*sqrt(g*yh/yFull)]^M
371 //
372 // for a given value of critical depth yc where:
373 // h = inlet depth above culvert invert
374 // s = culvert slope
375 // yFull = full depth of culvert
376 // yh = hydraulic depth at critical depth
377 // ac = flow area at critical depth
378 // g = accel. of gravity
379 // K and M = coefficients
380 //
381 {
382 double ac, wc, yh;
383 184453 TCulvert* culvert = (TCulvert *)p;
384
385 184453 ac = xsect_getAofY(culvert->xsect, yc);
386 184453 wc = xsect_getWofY(culvert->xsect, yc);
387 184453 yh = ac/wc;
388
389 184453 culvert->qc = ac * sqrt(GRAVITY * yh);
390 368906 return culvert->hPlus - yc/culvert->yFull - yh/2.0/culvert->yFull -
391 184453 culvert->kk * pow(culvert->qc/culvert->ad, culvert->mm);
392 }
393
394 //=============================================================================
395 /*
396 void report_CulvertControl(int j, double q0, double q, int condition, double yRatio)
397 //
398 // Used for debugging only
399 //
400 {
401 static char* conditionTxt[] = {"transition", "unsubmerged", "submerged"};
402 char theDate[12];
403 char theTime[9];
404 DateTime aDate = getDateTime(NewRoutingTime);
405 datetime_dateToStr(aDate, theDate);
406 datetime_timeToStr(aDate, theTime);
407 fprintf(Frpt.file,
408 "\n %11s: %8s Culvert %s flow reduced from %.3f to %.3f cfs for %s flow (%.2f).",
409 theDate, theTime, Link[j].ID, q0, q, conditionTxt[condition], yRatio);
410 }
411 */
412