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

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

added self.inflow to self.get_enquiry_height()
File size: 8.3 KB
Line 
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
10import culvert_routine
11from anuga.config import velocity_protection
12from anuga.utilities.numerical_tools import safe_acos as acos
13
14from math import pi, sqrt, sin, cos
15from anuga.config import g
16
17
18class Boyd_box_routine(culvert_routine.Culvert_routine):
19    """Boyd's generalisation of the US department of transportation culvert methods
20
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.
26
27        It is proposed to provide this in a seperate routine called
28        boyd_generalised_culvert_model_complex
29
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.
40
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
49
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.
53
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
56
57
58        """
59
60    def __init__(self, culvert):
61
62        culvert_routine.Culvert_routine.__init__(self, culvert)
63       
64        self.manning = culvert.manning
65
66
67
68    def __call__(self):
69
70        self.determine_inflow()
71
72        local_debug ='false'
73       
74        if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01:
75            if local_debug =='true':
76                log.critical('Specific E & Deltat Tot E = %s, %s'
77                             % (str(self.inflow.get_enquiry_specific_energy()),
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:
83                s = 'Specific energy  = %f m' % self.inflow.get_enquiry_specific_energy()
84                log_to_file(self.log_filename, s)
85
86            msg = 'Specific energy at inlet is negative'
87            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
88
89            if self.use_velocity_head :
90                driving_energy = self.inflow.get_enquiry_specific_energy()
91            else:
92                driving_energy = self.inflow.get_enquiry_height)_
93
94            height = self.culvert_height
95            width = self.culvert_width
96            flow_width = self.culvert_width
97
98            Q_inlet_unsubmerged = 0.540*g**0.5*width*driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
99            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
100
101            # FIXME(Ole): Are these functions really for inlet control?
102            if Q_inlet_unsubmerged < Q_inlet_submerged:
103                Q = Q_inlet_unsubmerged
104                dcrit = (Q**2/g/width**2)**0.333333
105                if dcrit > height:
106                    dcrit = height
107                flow_area = width*dcrit
108                outlet_culvert_depth = dcrit
109                case = 'Inlet unsubmerged Box Acts as Weir'
110            else:
111                Q = Q_inlet_submerged
112                flow_area = width*height
113                outlet_culvert_depth = height
114                case = 'Inlet submerged Box Acts as Orifice'
115
116            dcrit = (Q**2/g/width**2)**0.333333
117
118            outlet_culvert_depth = dcrit
119            if outlet_culvert_depth > height:
120                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
121                flow_area = width*height  # Cross sectional area of flow in the culvert
122                perimeter = 2*(width+height)
123                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
124            else:
125                flow_area = width * outlet_culvert_depth
126                perimeter = width+2*outlet_culvert_depth
127                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
128
129            if self.delta_total_energy < driving_energy:
130                # Calculate flows for outlet control
131
132                # Determine the depth at the outlet relative to the depth of flow in the Culvert
133                if self.outflow.get_enquiry_height() > height:        # The Outlet is Submerged
134                    outlet_culvert_depth=height
135                    flow_area=width*height       # Cross sectional area of flow in the culvert
136                    perimeter=2.0*(width+height)
137                    case = 'Outlet submerged'
138                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
139                    dcrit = (Q**2/g/width**2)**0.333333
140                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
141                    if outlet_culvert_depth > height:
142                        outlet_culvert_depth=height
143                        flow_area=width*height
144                        perimeter=2.0*(width+height)
145                        case = 'Outlet is Flowing Full'
146                    else:
147                        flow_area=width*outlet_culvert_depth
148                        perimeter=(width+2.0*outlet_culvert_depth)
149                        case = 'Outlet is open channel flow'
150
151                hyd_rad = flow_area/perimeter
152
153                if self.log_filename is not None:
154                    s = 'hydraulic radius at outlet = %f' % hyd_rad
155                    log_to_file(self.log_filename, s)
156
157                # Outlet control velocity using tail water
158                culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
159                Q_outlet_tailwater = flow_area * culvert_velocity
160
161                if self.log_filename is not None:
162                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
163                    log_to_file(self.log_filename, s)
164                Q = min(Q, Q_outlet_tailwater)
165            else:
166                pass
167                #FIXME(Ole): What about inlet control?
168
169            culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
170            if local_debug =='true':
171                log.critical('FLOW AREA = %s' % str(flow_area))
172                log.critical('PERIMETER = %s' % str(perimeter))
173                log.critical('Q final = %s' % str(Q))
174                log.critical('FROUDE = %s' % str(culv_froude))
175
176            # Determine momentum at the outlet
177            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
178
179        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
180
181        else: # self.inflow.get_enquiry_height() < 0.01:
182            Q = barrel_velocity = outlet_culvert_depth = 0.0
183
184        # Temporary flow limit
185        if barrel_velocity > self.max_velocity:
186            barrel_velocity = self.max_velocity
187            Q = flow_area * barrel_velocity
188
189
190
191       
192
193        return Q, barrel_velocity, outlet_culvert_depth
194
195
196
Note: See TracBrowser for help on using the repository browser.