source: trunk/anuga_core/source/anuga/structures/culvert_routines.py @ 7939

Last change on this file since 7939 was 7939, checked in by steve, 12 years ago

Adding in culvert routines into structures directory

File size: 18.2 KB
Line 
1"""Collection of culvert routines for use with Culvert_flow in culvert_class
2
3This module holds various routines to determine FLOW through CULVERTS and SIMPLE BRIDGES
4
5
6
7
8Usage:
9
10   
11
12"""
13
14#NOTE:
15# Inlet control:  Delta_total_energy > inlet_specific_energy
16# Outlet control: Delta_total_energy < inlet_specific_energy
17# where total energy is (w + 0.5*v^2/g) and
18# specific energy is (h + 0.5*v^2/g)
19
20
21from math import pi, sqrt, sin, cos
22
23
24def boyd_generalised_culvert_model(inlet_depth, 
25                                   outlet_depth,
26                                   inlet_velocity,
27                                   outlet_velocity,
28                                   inlet_specific_energy, 
29                                   delta_total_energy, 
30                                   g,
31                                   culvert_length=0.0,
32                                   culvert_width=0.0,
33                                   culvert_height=0.0,
34                                   culvert_type='box',
35                                   manning=0.0,
36                                   sum_loss=0.0,
37                                   max_velocity=10.0,
38                                   log_filename=None):
39
40    """Boyd's generalisation of the US department of transportation culvert methods
41       
42        WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
43        FULL PIPE OR CRITICAL DEPTH ONLY
44        For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
45        The obtain the CORRECT velocity requires an iteration of Depth to Establish
46        the Normal Depth of flow in the pipe.
47       
48        It is proposed to provide this in a seperate routine called
49    boyd_generalised_culvert_model_complex
50       
51    The Boyd Method is based on methods described by the following:
52        1.
53        US Dept. Transportation Federal Highway Administration (1965)
54        Hydraulic Chart for Selection of Highway Culverts.
55        Hydraulic Engineering Circular No. 5 US Government Printing
56        2.
57        US Dept. Transportation Federal Highway Administration (1972)
58        Capacity charts for the Hydraulic design of highway culverts.
59        Hydraulic Engineering Circular No. 10 US Government Printing
60    These documents provide around 60 charts for various configurations of culverts and inlets.
61       
62        Note these documents have been superceded by:
63        2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
64        Which combines culvert design information previously contained in Hydraulic Engineering Circulars
65        (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
66        HEC-5 provides 20 Charts
67        HEC-10 Provides an additional 36 Charts
68        HEC-13 Discusses the Design of improved more efficient inlets
69        HDS-5 Provides 60 sets of Charts
70
71        In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
72        1987 published "Generalised Head Discharge Equations for Culverts".
73        These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
74       
75        It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
76        The additional charts cover a range of culvert shapes and inlet configurations
77       
78   
79    """
80
81    from anuga.utilities.system_tools import log_to_file
82    from anuga.config import velocity_protection
83    from anuga.utilities.numerical_tools import safe_acos as acos
84
85    local_debug ='false'
86    if inlet_depth > 0.1: #this value was 0.01:
87        if local_debug =='true':
88            log.critical('Specific E & Deltat Tot E = %s, %s'
89                         % (str(inlet_specific_energy),
90                            str(delta_total_energy)))
91            log.critical('culvert type = %s' % str(culvert_type))
92        # Water has risen above inlet
93       
94        if log_filename is not None:                           
95            s = 'Specific energy  = %f m' % inlet_specific_energy
96            log_to_file(log_filename, s)
97       
98        msg = 'Specific energy at inlet is negative'
99        assert inlet_specific_energy >= 0.0, msg
100                     
101
102        if culvert_type =='circle':
103            # Round culvert (use height as diameter)
104            diameter = culvert_height   
105
106            """
107            For a circular pipe the Boyd method reviews 3 conditions
108            1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
109            2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
110            3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
111           
112            For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
113            """
114           
115            # Calculate flows for inlet control           
116            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
117            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
118            # Note for to SUBMERGED TO OCCUR inlet_specific_energy should be > 1.2 x diameter.... Should Check !!!
119
120            if log_filename is not None:                                       
121                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
122                log_to_file(log_filename, s)
123            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
124
125            # THE LOWEST Value will Control Calcs From here
126            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
127            dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
128            dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
129            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
130            if dcrit1/diameter  > 0.85:
131                outlet_culvert_depth = dcrit2
132            else:
133                outlet_culvert_depth = dcrit1
134            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
135            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
136            if outlet_culvert_depth >= diameter:
137                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
138                flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
139                perimeter = diameter * pi
140                flow_width= diameter
141                case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
142                if local_debug == 'true':
143                    log.critical('Inlet CTRL Outlet submerged Circular '
144                                 'PIPE FULL')
145            else:
146                #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
147                alpha = acos(1-2*outlet_culvert_depth/diameter)*2
148                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
149                flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
150                flow_width= diameter*sin(alpha/2.0)
151                perimeter = alpha*diameter/2.0
152                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
153                if local_debug =='true':
154                    log.critical('INLET CTRL Culvert is open channel flow '
155                                 'we will for now assume critical depth')
156                    log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
157                                 % (str(Q), str(outlet_culvert_depth),
158                                    str(alpha)))
159            if delta_total_energy < inlet_specific_energy:  #  OUTLET CONTROL !!!!
160                # Calculate flows for outlet control
161               
162                # Determine the depth at the outlet relative to the depth of flow in the Culvert
163                if outlet_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
164                    outlet_culvert_depth=diameter
165                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
166                    perimeter = diameter * pi
167                    flow_width= diameter
168                    case = 'Outlet submerged'
169                    if local_debug =='true':
170                        log.critical('Outlet submerged')
171                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
172                    # IF  outlet_depth < diameter
173                    dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
174                    dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
175                    if dcrit1/diameter >0.85:
176                        outlet_culvert_depth= dcrit2
177                    else:
178                        outlet_culvert_depth = dcrit1
179                    if outlet_culvert_depth > diameter:
180                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
181                        flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
182                        perimeter = diameter * pi
183                        flow_width= diameter
184                        case = 'Outlet unsubmerged PIPE FULL'
185                        if local_debug =='true':
186                            log.critical('Outlet unsubmerged PIPE FULL')
187                    else:
188                        alpha = acos(1-2*outlet_culvert_depth/diameter)*2
189                        flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
190                        flow_width= diameter*sin(alpha/2.0)
191                        perimeter = alpha*diameter/2.0
192                        case = 'Outlet is open channel flow we will for now assume critical depth'
193                        if local_debug == 'true':
194                            log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
195                                         % (str(Q), str(outlet_culvert_depth),
196                                            str(alpha)))
197                            log.critical('Outlet is open channel flow we '
198                                         'will for now assume critical depth')
199            if local_debug == 'true':
200                log.critical('FLOW AREA = %s' % str(flow_area))
201                log.critical('PERIMETER = %s' % str(perimeter))
202                log.critical('Q Interim = %s' % str(Q))
203            hyd_rad = flow_area/perimeter
204
205            if log_filename is not None:
206                s = 'hydraulic radius at outlet = %f' %hyd_rad
207                log_to_file(log_filename, s)
208
209            # Outlet control velocity using tail water
210            if local_debug =='true':
211                log.critical('GOT IT ALL CALCULATING Velocity')
212                log.critical('HydRad = %s' % str(hyd_rad))
213            culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
214            Q_outlet_tailwater = flow_area * culvert_velocity
215            if local_debug =='true':
216                log.critical('VELOCITY = %s' % str(culvert_velocity))
217                log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
218            if log_filename is not None:               
219                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
220                log_to_file(log_filename, s)
221            Q = min(Q, Q_outlet_tailwater)
222            if local_debug =='true':
223                log.critical('%s,%.3f,%.3f'
224                             % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
225                log.critical('%s,%.3f,%.3f,%.3f'
226                             % ('Q and Velocity and Depth=', Q,
227                                culvert_velocity, outlet_culvert_depth))
228
229        else:
230            pass
231                #FIXME(Ole): What about inlet control?
232        # ====  END OF CODE BLOCK FOR "IF" CIRCULAR PIPE
233
234        #else....
235        if culvert_type == 'box':
236            if local_debug == 'true':
237                log.critical('BOX CULVERT')
238            # Box culvert (rectangle or square)   ========================================================================================================================
239
240            # Calculate flows for inlet control
241            height = culvert_height
242            width = culvert_width
243            flow_width=culvert_width           
244           
245            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
246            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
247
248            if log_filename is not None:                           
249                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
250                log_to_file(log_filename, s)
251
252            # FIXME(Ole): Are these functions really for inlet control?   
253            if Q_inlet_unsubmerged < Q_inlet_submerged:
254                Q = Q_inlet_unsubmerged
255                dcrit = (Q**2/g/width**2)**0.333333
256                if dcrit > height:
257                    dcrit = height
258                flow_area = width*dcrit
259                outlet_culvert_depth = dcrit
260                case = 'Inlet unsubmerged Box Acts as Weir'
261            else:   
262                Q = Q_inlet_submerged
263                flow_area = width*height
264                outlet_culvert_depth = height
265                case = 'Inlet submerged Box Acts as Orifice'                   
266
267            dcrit = (Q**2/g/width**2)**0.333333
268
269            outlet_culvert_depth = dcrit
270            if outlet_culvert_depth > height:
271                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
272                flow_area = width*height  # Cross sectional area of flow in the culvert
273                perimeter = 2*(width+height)
274                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
275            else:
276                flow_area = width * outlet_culvert_depth
277                perimeter = width+2*outlet_culvert_depth
278                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
279       
280            if delta_total_energy < inlet_specific_energy:
281                # Calculate flows for outlet control
282               
283                # Determine the depth at the outlet relative to the depth of flow in the Culvert
284                if outlet_depth > height:        # The Outlet is Submerged
285                    outlet_culvert_depth=height
286                    flow_area=width*height       # Cross sectional area of flow in the culvert
287                    perimeter=2.0*(width+height)
288                    case = 'Outlet submerged'
289                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
290                    dcrit = (Q**2/g/width**2)**0.333333
291                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
292                    if outlet_culvert_depth > height:
293                        outlet_culvert_depth=height
294                        flow_area=width*height
295                        perimeter=2.0*(width+height)
296                        case = 'Outlet is Flowing Full'
297                    else:
298                        flow_area=width*outlet_culvert_depth
299                        perimeter=(width+2.0*outlet_culvert_depth)
300                        case = 'Outlet is open channel flow'
301
302                hyd_rad = flow_area/perimeter
303
304                if log_filename is not None:                               
305                    s = 'hydraulic radius at outlet = %f' % hyd_rad
306                    log_to_file(log_filename, s)
307
308                # Outlet control velocity using tail water
309                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
310                Q_outlet_tailwater = flow_area * culvert_velocity
311
312                if log_filename is not None:                           
313                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
314                    log_to_file(log_filename, s)
315                Q = min(Q, Q_outlet_tailwater)
316            else:
317                pass
318                #FIXME(Ole): What about inlet control?
319        # ====  END OF CODE BLOCK FOR "IF" BOX
320
321               
322        # Common code for circle and box geometries ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
323        if log_filename is not None:
324            log_to_file(log_filename, 'Case: "%s"' % case)
325           
326        if log_filename is not None:                       
327            s = 'Flow Rate Control = %f' % Q
328            log_to_file(log_filename, s)
329
330        culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
331        if local_debug =='true':
332            log.critical('FLOW AREA = %s' % str(flow_area))
333            log.critical('PERIMETER = %s' % str(perimeter))
334            log.critical('Q final = %s' % str(Q))
335            log.critical('FROUDE = %s' % str(culv_froude))
336        if log_filename is not None:                           
337            s = 'Froude in Culvert = %f' % culv_froude
338            log_to_file(log_filename, s)
339
340        # Determine momentum at the outlet
341        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
342
343    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow 
344
345    else: # inlet_depth < 0.01:
346        Q = barrel_velocity = outlet_culvert_depth = 0.0
347
348    # Temporary flow limit
349    if barrel_velocity > max_velocity:
350        if log_filename is not None:                           
351            s = 'Barrel velocity was reduced from = %f m/s to %f m/s' % (barrel_velocity, max_velocity)
352            log_to_file(log_filename, s)
353       
354        barrel_velocity = max_velocity
355        Q = flow_area * barrel_velocity
356       
357       
358
359       
360    return Q, barrel_velocity, outlet_culvert_depth
361
362
Note: See TracBrowser for help on using the repository browser.