source: trunk/anuga_core/source/anuga/structures/boyd_irregular_routine.py @ 7983

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

structure classes to deal with irregular structures

File size: 9.7 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_culvert_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, structure, manning=0.0):
61
62        culvert_routine.Culvert_routine.__init__(self, culvert, manning)
63       
64        self.slice_coordinates = structure.get_slice_coordinates()
65       
66
67    def __call__(self):
68
69        self.determine_inflow()
70
71        local_debug ='false'
72       
73        if self.inflow.get_average_height() > 0.01: #this value was 0.01:
74            if local_debug =='true':
75                log.critical('Specific E & Deltat Tot E = %s, %s'
76                             % (str(self.inflow.get_average_specific_energy()),
77                                str(self.delta_total_energy)))
78                log.critical('culvert type = %s' % str(culvert_type))
79            # Water has risen above inlet
80
81            if self.log_filename is not None:
82                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
83                log_to_file(self.log_filename, s)
84
85            msg = 'Specific energy at inlet is negative'
86            assert self.inflow.get_average_specific_energy() >= 0.0, msg
87
88            for  i in len(self.slice_coordinates):
89                for j in len(self.slice_coordinates[i]):
90                   
91                    slice_invert = self.slice_coordinates[i][j][1]
92                    slice_obvert = self.slice_coordinates[i][j][2]
93                   
94                    if self.outflow.get_average_height() > slice_obvert: # The slice is Submerged
95               
96                        #if self.outflow.get_average_height() > slice_invert: # The invert is wet
97               
98                        height = self.culvert_height
99                        width = self.culvert_width
100                        flow_width = self.culvert_width
101
102                        Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
103                        Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_average_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
104
105                        # FIXME(Ole): Are these functions really for inlet control?
106                        if Q_inlet_unsubmerged < Q_inlet_submerged:
107                            Q = Q_inlet_unsubmerged
108                            dcrit = (Q**2/g/width**2)**0.333333
109                            if dcrit > height:
110                                dcrit = height
111                            flow_area = width*dcrit
112                            outlet_culvert_depth = dcrit
113                            case = 'Inlet unsubmerged Box Acts as Weir'
114                        else:
115                            Q = Q_inlet_submerged
116                            flow_area = width*height
117                            outlet_culvert_depth = height
118                            case = 'Inlet submerged Box Acts as Orifice'
119
120                        dcrit = (Q**2/g/width**2)**0.333333
121
122                        outlet_culvert_depth = dcrit
123                        if outlet_culvert_depth > height:
124                            outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
125                            flow_area = width*height  # Cross sectional area of flow in the culvert
126                            perimeter = 2*(width+height)
127                            case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
128                        else:
129                            flow_area = width * outlet_culvert_depth
130                            perimeter = width+2*outlet_culvert_depth
131                            case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
132
133                        if self.delta_total_energy < self.inflow.get_average_specific_energy():
134                            # Calculate flows for outlet control
135
136                            # Determine the depth at the outlet relative to the depth of flow in the Culvert
137                            if self.outflow.get_average_height() > height:        # The Outlet is Submerged
138                                outlet_culvert_depth=height
139                                flow_area=width*height       # Cross sectional area of flow in the culvert
140                                perimeter=2.0*(width+height)
141                                case = 'Outlet submerged'
142                            else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
143                                dcrit = (Q**2/g/width**2)**0.333333
144                                outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
145                                if outlet_culvert_depth > height:
146                                    outlet_culvert_depth=height
147                                    flow_area=width*height
148                                    perimeter=2.0*(width+height)
149                                    case = 'Outlet is Flowing Full'
150                                else:
151                                    flow_area=width*outlet_culvert_depth
152                                    perimeter=(width+2.0*outlet_culvert_depth)
153                                    case = 'Outlet is open channel flow'
154
155                            hyd_rad = flow_area/perimeter
156
157                            if self.log_filename is not None:
158                                s = 'hydraulic radius at outlet = %f' % hyd_rad
159                                log_to_file(self.log_filename, s)
160
161                            # Outlet control velocity using tail water
162                            culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
163                            Q_outlet_tailwater = flow_area * culvert_velocity
164
165                            if self.log_filename is not None:
166                                s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
167                                log_to_file(self.log_filename, s)
168                            Q = min(Q, Q_outlet_tailwater)
169                        else:
170                            pass
171                            #FIXME(Ole): What about inlet control?
172
173                        culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
174                        if local_debug =='true':
175                            log.critical('FLOW AREA = %s' % str(flow_area))
176                            log.critical('PERIMETER = %s' % str(perimeter))
177                            log.critical('Q final = %s' % str(Q))
178                            log.critical('FROUDE = %s' % str(culv_froude))
179
180                        # Determine momentum at the outlet
181                        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
182
183                    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
184
185                    else: # self.inflow.get_average_height() < 0.01:
186                        Q = barrel_velocity = outlet_culvert_depth = 0.0
187
188                    # Temporary flow limit
189                    if barrel_velocity > self.max_velocity:
190                        barrel_velocity = self.max_velocity
191                        Q = flow_area * barrel_velocity
192
193
194
195       
196
197        return Q, barrel_velocity, outlet_culvert_depth
198
199
200
Note: See TracBrowser for help on using the repository browser.