source: trunk/anuga_core/source/anuga/structures/boyd_box_routine.py @ 7988

Last change on this file since 7988 was 7988, checked in by habili, 14 years ago

remove debugger

File size: 8.2 KB
RevLine 
[7969]1#! /usr/bin/python
2
3# To change this template, choose Tools | Templates
4# and open the template in the editor.
5
6__author__="steve"
7__date__ ="$23/08/2010 5:18:51 PM$"
8
9
[7980]10import culvert_routine
11from anuga.config import velocity_protection
12from anuga.utilities.numerical_tools import safe_acos as acos
[7969]13
[7980]14from math import pi, sqrt, sin, cos
15from anuga.config import g
16
17
18class Boyd_box_routine(culvert_routine.Culvert_routine):
[7969]19    """Boyd's generalisation of the US department of transportation culvert methods
20
[7980]21        WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
22        FULL PIPE OR CRITICAL DEPTH ONLY
23        For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
24        The obtain the CORRECT velocity requires an iteration of Depth to Establish
25        the Normal Depth of flow in the pipe.
[7969]26
[7980]27        It is proposed to provide this in a seperate routine called
28        boyd_generalised_culvert_model_complex
[7969]29
[7980]30        The Boyd Method is based on methods described by the following:
31        1.
32        US Dept. Transportation Federal Highway Administration (1965)
33        Hydraulic Chart for Selection of Highway Culverts.
34        Hydraulic Engineering Circular No. 5 US Government Printing
35        2.
36        US Dept. Transportation Federal Highway Administration (1972)
37        Capacity charts for the Hydraulic design of highway culverts.
38        Hydraulic Engineering Circular No. 10 US Government Printing
39        These documents provide around 60 charts for various configurations of culverts and inlets.
[7969]40
[7980]41        Note these documents have been superceded by:
42        2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
43        Which combines culvert design information previously contained in Hydraulic Engineering Circulars
44        (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
45        HEC-5 provides 20 Charts
46        HEC-10 Provides an additional 36 Charts
47        HEC-13 Discusses the Design of improved more efficient inlets
48        HDS-5 Provides 60 sets of Charts
[7969]49
[7980]50        In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
51        1987 published "Generalised Head Discharge Equations for Culverts".
52        These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
[7969]53
[7980]54        It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
55        The additional charts cover a range of culvert shapes and inlet configurations
[7969]56
57
[7980]58        """
[7969]59
[7986]60    def __init__(self, culvert):
[7969]61
[7986]62        culvert_routine.Culvert_routine.__init__(self, culvert)
63       
64        self.manning = culvert.manning
[7969]65
66
67
[7980]68    def __call__(self):
[7969]69
[7980]70        self.determine_inflow()
[7969]71
[7980]72        local_debug ='false'
73       
[7984]74        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01:
[7980]75            if local_debug =='true':
76                log.critical('Specific E & Deltat Tot E = %s, %s'
[7984]77                             % (str(self.inflow.get_enquiry_specific_energy()),
[7980]78                                str(self.delta_total_energy)))
79                log.critical('culvert type = %s' % str(culvert_type))
80            # Water has risen above inlet
81
82            if self.log_filename is not None:
[7984]83                s = 'Specific energy  = %f m' % self.inflow.get_enquiry_specific_energy()
[7980]84                log_to_file(self.log_filename, s)
85
86            msg = 'Specific energy at inlet is negative'
[7984]87            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
[7980]88
89            height = self.culvert_height
90            width = self.culvert_width
91            flow_width = self.culvert_width
92
[7984]93            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_enquiry_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
94            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_enquiry_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
[7980]95
96            # FIXME(Ole): Are these functions really for inlet control?
97            if Q_inlet_unsubmerged < Q_inlet_submerged:
98                Q = Q_inlet_unsubmerged
99                dcrit = (Q**2/g/width**2)**0.333333
100                if dcrit > height:
101                    dcrit = height
102                flow_area = width*dcrit
103                outlet_culvert_depth = dcrit
104                case = 'Inlet unsubmerged Box Acts as Weir'
105            else:
106                Q = Q_inlet_submerged
107                flow_area = width*height
108                outlet_culvert_depth = height
109                case = 'Inlet submerged Box Acts as Orifice'
110
[7969]111            dcrit = (Q**2/g/width**2)**0.333333
[7980]112
113            outlet_culvert_depth = dcrit
114            if outlet_culvert_depth > height:
115                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
116                flow_area = width*height  # Cross sectional area of flow in the culvert
117                perimeter = 2*(width+height)
118                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
[7969]119            else:
[7980]120                flow_area = width * outlet_culvert_depth
121                perimeter = width+2*outlet_culvert_depth
122                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
[7969]123
[7984]124            if self.delta_total_energy < self.inflow.get_enquiry_specific_energy():
[7980]125                # Calculate flows for outlet control
[7969]126
[7980]127                # Determine the depth at the outlet relative to the depth of flow in the Culvert
[7984]128                if self.outflow.get_enquiry_height() > height:        # The Outlet is Submerged
[7980]129                    outlet_culvert_depth=height
130                    flow_area=width*height       # Cross sectional area of flow in the culvert
131                    perimeter=2.0*(width+height)
132                    case = 'Outlet submerged'
133                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
134                    dcrit = (Q**2/g/width**2)**0.333333
135                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
136                    if outlet_culvert_depth > height:
137                        outlet_culvert_depth=height
138                        flow_area=width*height
139                        perimeter=2.0*(width+height)
140                        case = 'Outlet is Flowing Full'
141                    else:
142                        flow_area=width*outlet_culvert_depth
143                        perimeter=(width+2.0*outlet_culvert_depth)
144                        case = 'Outlet is open channel flow'
[7969]145
[7980]146                hyd_rad = flow_area/perimeter
[7969]147
[7980]148                if self.log_filename is not None:
149                    s = 'hydraulic radius at outlet = %f' % hyd_rad
150                    log_to_file(self.log_filename, s)
[7969]151
[7980]152                # Outlet control velocity using tail water
153                culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
154                Q_outlet_tailwater = flow_area * culvert_velocity
[7969]155
[7980]156                if self.log_filename is not None:
157                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
158                    log_to_file(self.log_filename, s)
159                Q = min(Q, Q_outlet_tailwater)
160            else:
161                pass
162                #FIXME(Ole): What about inlet control?
[7969]163
[7980]164            culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
165            if local_debug =='true':
166                log.critical('FLOW AREA = %s' % str(flow_area))
167                log.critical('PERIMETER = %s' % str(perimeter))
168                log.critical('Q final = %s' % str(Q))
169                log.critical('FROUDE = %s' % str(culv_froude))
[7969]170
[7980]171            # Determine momentum at the outlet
172            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
[7969]173
[7980]174        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
[7969]175
[7984]176        else: # self.inflow.get_enquiry_height() < 0.01:
[7980]177            Q = barrel_velocity = outlet_culvert_depth = 0.0
[7969]178
[7980]179        # Temporary flow limit
180        if barrel_velocity > self.max_velocity:
181            barrel_velocity = self.max_velocity
182            Q = flow_area * barrel_velocity
[7969]183
184
185
[7980]186       
[7969]187
[7980]188        return Q, barrel_velocity, outlet_culvert_depth
[7969]189
[7980]190
191
Note: See TracBrowser for help on using the repository browser.