source: trunk/anuga_core/source/anuga/structures/boyd_circle_routine.py @ 7980

Last change on this file since 7980 was 7980, checked in by steve, 12 years ago

Added new tests

File size: 11.9 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__ ="$30/08/2010 10:15:08 AM$"
8
9import culvert_routine
10from anuga.config import velocity_protection
11from anuga.utilities.numerical_tools import safe_acos as acos
12
13from math import pi, sqrt, sin, cos
14from anuga.config import g
15
16
17class Boyd_box_routine(culvert_routine.Culvert_routine):
18    """Boyd's generalisation of the US department of transportation culvert methods
19
20        WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
21        FULL PIPE OR CRITICAL DEPTH ONLY
22        For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
23        The obtain the CORRECT velocity requires an iteration of Depth to Establish
24        the Normal Depth of flow in the pipe.
25
26        It is proposed to provide this in a seperate routine called
27        boyd_generalised_culvert_model_complex
28
29        The Boyd Method is based on methods described by the following:
30        1.
31        US Dept. Transportation Federal Highway Administration (1965)
32        Hydraulic Chart for Selection of Highway Culverts.
33        Hydraulic Engineering Circular No. 5 US Government Printing
34        2.
35        US Dept. Transportation Federal Highway Administration (1972)
36        Capacity charts for the Hydraulic design of highway culverts.
37        Hydraulic Engineering Circular No. 10 US Government Printing
38        These documents provide around 60 charts for various configurations of culverts and inlets.
39
40        Note these documents have been superceded by:
41        2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
42        Which combines culvert design information previously contained in Hydraulic Engineering Circulars
43        (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
44        HEC-5 provides 20 Charts
45        HEC-10 Provides an additional 36 Charts
46        HEC-13 Discusses the Design of improved more efficient inlets
47        HDS-5 Provides 60 sets of Charts
48
49        In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
50        1987 published "Generalised Head Discharge Equations for Culverts".
51        These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
52
53        It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
54        The additional charts cover a range of culvert shapes and inlet configurations
55
56
57        """
58
59    def __init__(self):
60
61        Culvert_routine.__init__(self)
62
63
64
65    def __call__(self):
66
67        """
68        For a circular pipe the Boyd method reviews 3 conditions
69        1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
70        2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
71        3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
72
73        For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
74        """
75
76        diameter = self.culvert_height
77
78        local_debug ='false'
79        if self.inflow.get_average_height() > 0.1: #this value was 0.01:
80            if local_debug =='true':
81                log.critical('Specific E & Deltat Tot E = %s, %s'
82                             % (str(self.inflow.get_average_specific_energy()),
83                                str(self.delta_total_energy)))
84                log.critical('culvert type = %s' % str(culvert_type))
85            # Water has risen above inlet
86
87            if self.log_filename is not None:
88                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
89                log_to_file(self.log_filename, s)
90
91            msg = 'Specific energy at inlet is negative'
92            assert self.inflow.get_average_specific_energy() >= 0.0, msg
93
94            # Calculate flows for inlet control
95            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*self.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged
96            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*self.inflow.get_average_specific_energy()**0.63   # Inlet Ctrl Inlet Submerged
97            # Note for to SUBMERGED TO OCCUR self.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
98
99            if self.log_filename is not None:
100                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
101                log_to_file(self.log_filename, s)
102            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
103
104            # THE LOWEST Value will Control Calcs From here
105            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
106            dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
107            dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
108            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
109            if dcrit1/diameter  > 0.85:
110                outlet_culvert_depth = dcrit2
111            else:
112                outlet_culvert_depth = dcrit1
113            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
114            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
115            if outlet_culvert_depth >= diameter:
116                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
117                flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
118                perimeter = diameter * pi
119                flow_width= diameter
120                case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
121                if local_debug == 'true':
122                    log.critical('Inlet CTRL Outlet submerged Circular '
123                                 'PIPE FULL')
124            else:
125                #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
126                alpha = acos(1-2*outlet_culvert_depth/diameter)*2
127                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
128                flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
129                flow_width= diameter*sin(alpha/2.0)
130                perimeter = alpha*diameter/2.0
131                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
132                if local_debug =='true':
133                    log.critical('INLET CTRL Culvert is open channel flow '
134                                 'we will for now assume critical depth')
135                    log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
136                                 % (str(Q), str(outlet_culvert_depth),
137                                    str(alpha)))
138            if self.delta_total_energy < self.inflow.get_average_specific_energy():  #  OUTLET CONTROL !!!!
139                # Calculate flows for outlet control
140
141                # Determine the depth at the outlet relative to the depth of flow in the Culvert
142                if self.outflow.get_average_height() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
143                    outlet_culvert_depth=diameter
144                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
145                    perimeter = diameter * pi
146                    flow_width= diameter
147                    case = 'Outlet submerged'
148                    if local_debug =='true':
149                        log.critical('Outlet submerged')
150                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
151                    # IF  self.outflow.get_average_height() < diameter
152                    dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
153                    dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
154                    if dcrit1/diameter >0.85:
155                        outlet_culvert_depth= dcrit2
156                    else:
157                        outlet_culvert_depth = dcrit1
158                    if outlet_culvert_depth > diameter:
159                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
160                        flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
161                        perimeter = diameter * pi
162                        flow_width= diameter
163                        case = 'Outlet unsubmerged PIPE FULL'
164                        if local_debug =='true':
165                            log.critical('Outlet unsubmerged PIPE FULL')
166                    else:
167                        alpha = acos(1-2*outlet_culvert_depth/diameter)*2
168                        flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
169                        flow_width= diameter*sin(alpha/2.0)
170                        perimeter = alpha*diameter/2.0
171                        case = 'Outlet is open channel flow we will for now assume critical depth'
172                        if local_debug == 'true':
173                            log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
174                                         % (str(Q), str(outlet_culvert_depth),
175                                            str(alpha)))
176                            log.critical('Outlet is open channel flow we '
177                                         'will for now assume critical depth')
178            if local_debug == 'true':
179                log.critical('FLOW AREA = %s' % str(flow_area))
180                log.critical('PERIMETER = %s' % str(perimeter))
181                log.critical('Q Interim = %s' % str(Q))
182            hyd_rad = flow_area/perimeter
183
184            if self.log_filename is not None:
185                s = 'hydraulic radius at outlet = %f' %hyd_rad
186                log_to_file(self.log_filename, s)
187
188            # Outlet control velocity using tail water
189            if local_debug =='true':
190                log.critical('GOT IT ALL CALCULATING Velocity')
191                log.critical('HydRad = %s' % str(hyd_rad))
192            culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
193            Q_outlet_tailwater = flow_area * culvert_velocity
194            if local_debug =='true':
195                log.critical('VELOCITY = %s' % str(culvert_velocity))
196                log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
197            if self.log_filename is not None:
198                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
199                log_to_file(self.log_filename, s)
200            Q = min(Q, Q_outlet_tailwater)
201            if local_debug =='true':
202                log.critical('%s,%.3f,%.3f'
203                             % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
204                log.critical('%s,%.3f,%.3f,%.3f'
205                             % ('Q and Velocity and Depth=', Q,
206                                culvert_velocity, outlet_culvert_depth))
207
208            culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
209            if local_debug =='true':
210                log.critical('FLOW AREA = %s' % str(flow_area))
211                log.critical('PERIMETER = %s' % str(perimeter))
212                log.critical('Q final = %s' % str(Q))
213                log.critical('FROUDE = %s' % str(culv_froude))
214
215            # Determine momentum at the outlet
216            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
217
218        else: # self.inflow.get_average_height() < 0.01:
219            Q = barrel_velocity = outlet_culvert_depth = 0.0
220
221        # Temporary flow limit
222        if barrel_velocity > self.max_velocity:
223            barrel_velocity = self.max_velocity
224            Q = flow_area * barrel_velocity
225
226        return Q, barrel_velocity, outlet_culvert_depth
227
228
229
Note: See TracBrowser for help on using the repository browser.