source: branches/numpy/anuga/culvert_flows/culvert_routines.py @ 7184

Last change on this file since 7184 was 7176, checked in by rwilson, 16 years ago

Back-merge from Numeric to numpy.
File size: 17.5 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            print 'Specific E & Deltat Tot E = ',inlet_specific_energy,delta_total_energy
89            print 'culvert typ = ',culvert_type
90        # Water has risen above inlet
91       
92        if log_filename is not None:                           
93            s = 'Specific energy  = %f m' % inlet_specific_energy
94            log_to_file(log_filename, s)
95       
96        msg = 'Specific energy at inlet is negative'
97        assert inlet_specific_energy >= 0.0, msg
98                     
99
100        if culvert_type =='circle':
101            # Round culvert (use height as diameter)
102            diameter = culvert_height   
103
104            """
105            For a circular pipe the Boyd method reviews 3 conditions
106            1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
107            2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
108            3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
109           
110            For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
111            """
112           
113            # Calculate flows for inlet control           
114            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
115            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
116            # Note for to SUBMERGED TO OCCUR inlet_specific_energy should be > 1.2 x diameter.... Should Check !!!
117
118            if log_filename is not None:                                       
119                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
120                log_to_file(log_filename, s)
121            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
122
123            # THE LOWEST Value will Control Calcs From here
124            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
125            dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
126            dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
127            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
128            if dcrit1/diameter  > 0.85:
129                outlet_culvert_depth = dcrit2
130            else:
131                outlet_culvert_depth = dcrit1
132            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
133            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
134            if outlet_culvert_depth >= diameter:
135                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
136                flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
137                perimeter = diameter * pi
138                flow_width= diameter
139                case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
140                if local_debug =='true':
141                    print 'Inlet CTRL Outlet submerged Circular PIPE FULL'
142            else:
143                #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
144                alpha = acos(1-2*outlet_culvert_depth/diameter)*2
145                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
146                flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
147                flow_width= diameter*sin(alpha/2.0)
148                perimeter = alpha*diameter/2.0
149                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
150                if local_debug =='true':
151                    print 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
152                    print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,'  ',alpha
153            if delta_total_energy < inlet_specific_energy:  #  OUTLET CONTROL !!!!
154                # Calculate flows for outlet control
155               
156                # Determine the depth at the outlet relative to the depth of flow in the Culvert
157                if outlet_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
158                    outlet_culvert_depth=diameter
159                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
160                    perimeter = diameter * pi
161                    flow_width= diameter
162                    case = 'Outlet submerged'
163                    if local_debug =='true':
164                        print 'Outlet submerged'
165                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
166                    # IF  outlet_depth < diameter
167                    dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
168                    dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
169                    if dcrit1/diameter >0.85:
170                        outlet_culvert_depth= dcrit2
171                    else:
172                        outlet_culvert_depth = dcrit1
173                    if outlet_culvert_depth > diameter:
174                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
175                        flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
176                        perimeter = diameter * pi
177                        flow_width= diameter
178                        case = 'Outlet unsubmerged PIPE FULL'
179                        if local_debug =='true':
180                            print  'Outlet unsubmerged PIPE FULL'
181                    else:
182                        alpha = acos(1-2*outlet_culvert_depth/diameter)*2
183                        flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
184                        flow_width= diameter*sin(alpha/2.0)
185                        perimeter = alpha*diameter/2.0
186                        case = 'Outlet is open channel flow we will for now assume critical depth'
187                        if local_debug =='true':
188                            print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,'  ',alpha
189                            print 'Outlet is open channel flow we will for now assume critical depth'
190            if local_debug =='true':
191                print 'FLOW AREA = ',flow_area
192                print 'PERIMETER = ',perimeter
193                print 'Q Interim = ',Q
194            hyd_rad = flow_area/perimeter
195
196            if log_filename is not None:
197                s = 'hydraulic radius at outlet = %f' %hyd_rad
198                log_to_file(log_filename, s)
199
200            # Outlet control velocity using tail water
201            if local_debug =='true':
202                print 'GOT IT ALL CALCULATING Velocity'
203                print 'HydRad = ',hyd_rad     
204            culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
205            Q_outlet_tailwater = flow_area * culvert_velocity
206            if local_debug =='true':
207                print 'VELOCITY = ',culvert_velocity
208                print 'Outlet Ctrl Q = ',Q_outlet_tailwater
209            if log_filename is not None:               
210                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
211                log_to_file(log_filename, s)
212            Q = min(Q, Q_outlet_tailwater)
213            if local_debug =='true':
214                print ('%s,%.3f,%.3f' %('dcrit 1 , dcit2 =',dcrit1,dcrit2))
215                print ('%s,%.3f,%.3f,%.3f' %('Q and Velocity and Depth=',Q,culvert_velocity,outlet_culvert_depth))
216
217        else:
218            pass
219                #FIXME(Ole): What about inlet control?
220        # ====  END OF CODE BLOCK FOR "IF" CIRCULAR PIPE
221
222        #else....
223        if culvert_type =='box':
224            if local_debug =='true':
225                print 'BOX CULVERT'
226            # Box culvert (rectangle or square)   ========================================================================================================================
227
228            # Calculate flows for inlet control
229            height = culvert_height
230            width = culvert_width
231            flow_width=culvert_width           
232           
233            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
234            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
235
236            if log_filename is not None:                           
237                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
238                log_to_file(log_filename, s)
239
240            # FIXME(Ole): Are these functions really for inlet control?   
241            if Q_inlet_unsubmerged < Q_inlet_submerged:
242                Q = Q_inlet_unsubmerged
243                dcrit = (Q**2/g/width**2)**0.333333
244                if dcrit > height:
245                    dcrit = height
246                flow_area = width*dcrit
247                outlet_culvert_depth = dcrit
248                case = 'Inlet unsubmerged Box Acts as Weir'
249            else:   
250                Q = Q_inlet_submerged
251                flow_area = width*height
252                outlet_culvert_depth = height
253                case = 'Inlet submerged Box Acts as Orifice'                   
254
255            dcrit = (Q**2/g/width**2)**0.333333
256
257            outlet_culvert_depth = dcrit
258            if outlet_culvert_depth > height:
259                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
260                flow_area = width*height  # Cross sectional area of flow in the culvert
261                perimeter = 2*(width+height)
262                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
263            else:
264                flow_area = width * outlet_culvert_depth
265                perimeter = width+2*outlet_culvert_depth
266                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
267       
268            if delta_total_energy < inlet_specific_energy:
269                # Calculate flows for outlet control
270               
271                # Determine the depth at the outlet relative to the depth of flow in the Culvert
272                if outlet_depth > height:        # The Outlet is Submerged
273                    outlet_culvert_depth=height
274                    flow_area=width*height       # Cross sectional area of flow in the culvert
275                    perimeter=2.0*(width+height)
276                    case = 'Outlet submerged'
277                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
278                    dcrit = (Q**2/g/width**2)**0.333333
279                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
280                    if outlet_culvert_depth > height:
281                        outlet_culvert_depth=height
282                        flow_area=width*height
283                        perimeter=2.0*(width+height)
284                        case = 'Outlet is Flowing Full'
285                    else:
286                        flow_area=width*outlet_culvert_depth
287                        perimeter=(width+2.0*outlet_culvert_depth)
288                        case = 'Outlet is open channel flow'
289
290                hyd_rad = flow_area/perimeter
291
292                if log_filename is not None:                               
293                    s = 'hydraulic radius at outlet = %f' % hyd_rad
294                    log_to_file(log_filename, s)
295
296                # Outlet control velocity using tail water
297                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
298                Q_outlet_tailwater = flow_area * culvert_velocity
299
300                if log_filename is not None:                           
301                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
302                    log_to_file(log_filename, s)
303                Q = min(Q, Q_outlet_tailwater)
304            else:
305                pass
306                #FIXME(Ole): What about inlet control?
307        # ====  END OF CODE BLOCK FOR "IF" BOX
308
309               
310        # Common code for circle and box geometries ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
311        if log_filename is not None:
312            log_to_file(log_filename, 'Case: "%s"' % case)
313           
314        if log_filename is not None:                       
315            s = 'Flow Rate Control = %f' % Q
316            log_to_file(log_filename, s)
317
318        culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
319        if local_debug =='true':
320            print 'FLOW AREA = ',flow_area
321            print 'PERIMETER = ',perimeter
322            print 'Q final = ',Q
323            print 'FROUDE = ',culv_froude
324        if log_filename is not None:                           
325            s = 'Froude in Culvert = %f' % culv_froude
326            log_to_file(log_filename, s)
327
328        # Determine momentum at the outlet
329        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
330
331    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow 
332
333    else: # inlet_depth < 0.01:
334        Q = barrel_velocity = outlet_culvert_depth = 0.0
335
336    # Temporary flow limit
337    if barrel_velocity > max_velocity:
338        if log_filename is not None:                           
339            s = 'Barrel velocity was reduced from = %f m/s to %f m/s' % (barrel_velocity, max_velocity)
340            log_to_file(log_filename, s)
341       
342        barrel_velocity = max_velocity
343        Q = flow_area * barrel_velocity
344       
345       
346
347       
348    return Q, barrel_velocity, outlet_culvert_depth
349
350
Note: See TracBrowser for help on using the repository browser.