- Timestamp:
- Aug 30, 2010, 5:53:45 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/structures/boyd_box_routine.py
r7977 r7980 8 8 9 9 10 import culvert_routine 11 from anuga.config import velocity_protection 12 from anuga.utilities.numerical_tools import safe_acos as acos 10 13 11 def boyd_box(culvert): 14 from math import pi, sqrt, sin, cos 15 from anuga.config import g 16 17 18 class Boyd_box_routine(culvert_routine.Culvert_routine): 12 19 """Boyd's generalisation of the US department of transportation culvert methods 13 20 14 15 16 17 18 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. 19 26 20 21 boyd_generalised_culvert_model_complex27 It is proposed to provide this in a seperate routine called 28 boyd_generalised_culvert_model_complex 22 29 23 The Boyd Method is based on methods described by the following:24 25 26 27 28 29 30 31 32 These documents provide around 60 charts for various configurations of culverts and inlets.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. 33 40 34 35 36 37 38 39 40 41 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 42 49 43 44 45 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. 46 53 47 It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done. 48 The additional charts cover a range of culvert shapes and inlet configurations 49 50 """ 51 52 # Calculate flows for inflow control 53 54 Q_inflow_unsubmerged = 0.540*g**0.5*width*inflow_specific_energy**1.50 # Flow based on inflow Ctrl inflow Unsubmerged 55 Q_inflow_submerged = 0.702*g**0.5*width*height**0.89*inflow_specific_energy**0.61 # Flow based on inflow Ctrl inflow Submerged 56 57 if log_filename is not None: 58 s = 'Q_inflow_unsubmerged = %.6f, Q_inflow_submerged = %.6f' %(Q_inflow_unsubmerged, Q_inflow_submerged) 59 log_to_file(log_filename, s) 60 61 # FIXME(Ole): Are these functions really for inflow control? 62 if Q_inflow_unsubmerged < Q_inflow_submerged: 63 Q = Q_inflow_unsubmerged 64 dcrit = (Q**2/g/width**2)**0.333333 65 if dcrit > height: 66 dcrit = height 67 flow_area = width*dcrit 68 outflow_culvert_depth = dcrit 69 case = 'inflow unsubmerged Box Acts as Weir' 70 else: 71 Q = Q_inflow_submerged 72 flow_area = width*height 73 outflow_culvert_depth = height 74 case = 'inflow submerged Box Acts as Orifice' 75 76 dcrit = (Q**2/g/width**2)**0.333333 77 78 outflow_culvert_depth = dcrit 79 if outflow_culvert_depth > height: 80 outflow_culvert_depth = height # Once again the pipe is flowing full not partfull 81 flow_area = width*height # Cross sectional area of flow in the culvert 82 perimeter = 2*(width+height) 83 case = 'inflow CTRL outflow unsubmerged PIPE PART FULL' 84 else: 85 flow_area = width * outflow_culvert_depth 86 perimeter = width+2*outflow_culvert_depth 87 case = 'inflow CTRL Culvert is open channel flow we will for now assume critical depth' 88 89 if delta_total_energy < inflow_specific_energy: 90 # Calculate flows for outflow control 91 92 # Determine the depth at the outflow relative to the depth of flow in the Culvert 93 if outflow_depth > height: # The outflow is Submerged 94 outflow_culvert_depth=height 95 flow_area=width*height # Cross sectional area of flow in the culvert 96 perimeter=2.0*(width+height) 97 case = 'outflow submerged' 98 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 99 dcrit = (Q**2/g/width**2)**0.333333 100 outflow_culvert_depth=dcrit # For purpose of calculation assume the outflow depth = Critical Depth 101 if outflow_culvert_depth > height: 102 outflow_culvert_depth=height 103 flow_area=width*height 104 perimeter=2.0*(width+height) 105 case = 'outflow is Flowing Full' 106 else: 107 flow_area=width*outflow_culvert_depth 108 perimeter=(width+2.0*outflow_culvert_depth) 109 case = 'outflow is open channel flow' 110 111 hyd_rad = flow_area/perimeter 112 113 if log_filename is not None: 114 s = 'hydraulic radius at outflow = %f' % hyd_rad 115 log_to_file(log_filename, s) 116 117 # outflow control velocity using tail water 118 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 119 Q_outflow_tailwater = flow_area * culvert_velocity 120 121 if log_filename is not None: 122 s = 'Q_outflow_tailwater = %.6f' % Q_outflow_tailwater 123 log_to_file(log_filename, s) 124 Q = min(Q, Q_outflow_tailwater) 125 126 return Q 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 127 56 128 57 129 if __name__ == "__main__": 58 """ 59 60 def __init__(self, culvert, manning=0.0): 61 62 culvert_routine.Culvert_routine.__init__(self, culvert, manning) 130 63 131 64 132 g=9.81133 culvert_slope=0.1 # Downward134 65 135 inlet_depth=2.0 136 outlet_depth=0.0 66 def __call__(self): 137 67 138 inlet_velocity=0.0, 139 outlet_velocity=0.0, 68 self.determine_inflow() 140 69 141 culvert_length=4.0 142 culvert_width=1.2 143 culvert_height=0.75 70 local_debug ='false' 71 72 if self.inflow.get_average_height() > 0.01: #this value was 0.01: 73 if local_debug =='true': 74 log.critical('Specific E & Deltat Tot E = %s, %s' 75 % (str(self.inflow.get_average_specific_energy()), 76 str(self.delta_total_energy))) 77 log.critical('culvert type = %s' % str(culvert_type)) 78 # Water has risen above inlet 144 79 145 culvert_type='box'146 manning=0.013147 sum_loss=0.080 if self.log_filename is not None: 81 s = 'Specific energy = %f m' % self.inflow.get_average_specific_energy() 82 log_to_file(self.log_filename, s) 148 83 149 inlet_specific_energy=inlet_depth #+0.5*v**2/g 150 z_in = 0.0 151 z_out = -culvert_length*culvert_slope/100 152 E_in = z_in+inlet_depth # + 153 E_out = z_out+outlet_depth # + 154 delta_total_energy = E_in-E_out 84 msg = 'Specific energy at inlet is negative' 85 assert self.inflow.get_average_specific_energy() >= 0.0, msg 155 86 156 Q = boyd_box(culvert_height, culvert_width, culvert_width, inlet_specific_energy) 87 height = self.culvert_height 88 width = self.culvert_width 89 flow_width = self.culvert_width 157 90 158 print 'Q ',Q 91 Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged 92 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 93 94 # FIXME(Ole): Are these functions really for inlet control? 95 if Q_inlet_unsubmerged < Q_inlet_submerged: 96 Q = Q_inlet_unsubmerged 97 dcrit = (Q**2/g/width**2)**0.333333 98 if dcrit > height: 99 dcrit = height 100 flow_area = width*dcrit 101 outlet_culvert_depth = dcrit 102 case = 'Inlet unsubmerged Box Acts as Weir' 103 else: 104 Q = Q_inlet_submerged 105 flow_area = width*height 106 outlet_culvert_depth = height 107 case = 'Inlet submerged Box Acts as Orifice' 108 109 dcrit = (Q**2/g/width**2)**0.333333 110 111 outlet_culvert_depth = dcrit 112 if outlet_culvert_depth > height: 113 outlet_culvert_depth = height # Once again the pipe is flowing full not partfull 114 flow_area = width*height # Cross sectional area of flow in the culvert 115 perimeter = 2*(width+height) 116 case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL' 117 else: 118 flow_area = width * outlet_culvert_depth 119 perimeter = width+2*outlet_culvert_depth 120 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth' 121 122 if self.delta_total_energy < self.inflow.get_average_specific_energy(): 123 # Calculate flows for outlet control 124 125 # Determine the depth at the outlet relative to the depth of flow in the Culvert 126 if self.outflow.get_average_height() > height: # The Outlet is Submerged 127 outlet_culvert_depth=height 128 flow_area=width*height # Cross sectional area of flow in the culvert 129 perimeter=2.0*(width+height) 130 case = 'Outlet submerged' 131 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 132 dcrit = (Q**2/g/width**2)**0.333333 133 outlet_culvert_depth=dcrit # For purpose of calculation assume the outlet depth = Critical Depth 134 if outlet_culvert_depth > height: 135 outlet_culvert_depth=height 136 flow_area=width*height 137 perimeter=2.0*(width+height) 138 case = 'Outlet is Flowing Full' 139 else: 140 flow_area=width*outlet_culvert_depth 141 perimeter=(width+2.0*outlet_culvert_depth) 142 case = 'Outlet is open channel flow' 143 144 hyd_rad = flow_area/perimeter 145 146 if self.log_filename is not None: 147 s = 'hydraulic radius at outlet = %f' % hyd_rad 148 log_to_file(self.log_filename, s) 149 150 # Outlet control velocity using tail water 151 culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333)) 152 Q_outlet_tailwater = flow_area * culvert_velocity 153 154 if self.log_filename is not None: 155 s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater 156 log_to_file(self.log_filename, s) 157 Q = min(Q, Q_outlet_tailwater) 158 else: 159 pass 160 #FIXME(Ole): What about inlet control? 161 162 culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3)) 163 if local_debug =='true': 164 log.critical('FLOW AREA = %s' % str(flow_area)) 165 log.critical('PERIMETER = %s' % str(perimeter)) 166 log.critical('Q final = %s' % str(Q)) 167 log.critical('FROUDE = %s' % str(culv_froude)) 168 169 # Determine momentum at the outlet 170 barrel_velocity = Q/(flow_area + velocity_protection/flow_area) 171 172 # END CODE BLOCK for DEPTH > Required depth for CULVERT Flow 173 174 else: # self.inflow.get_average_height() < 0.01: 175 Q = barrel_velocity = outlet_culvert_depth = 0.0 176 177 # Temporary flow limit 178 if barrel_velocity > self.max_velocity: 179 barrel_velocity = self.max_velocity 180 Q = flow_area * barrel_velocity 181 182 183 184 185 186 return Q, barrel_velocity, outlet_culvert_depth 187 188 189
Note: See TracChangeset
for help on using the changeset viewer.
