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

Last change on this file since 7977 was 7977, checked in by habili, 12 years ago

Refactoring of routines.

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