Changeset 7073
 Timestamp:
 May 25, 2009, 2:58:05 PM (11 years ago)
 Location:
 anuga_core/source/anuga/culvert_flows
 Files:

 4 added
 1 deleted
 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/culvert_flows/culvert_routines.py
r6948 r7073 1 1 """Collection of culvert routines for use with Culvert_flow in culvert_class 2 3 This module holds various riutines to determine FLOW through CULVERTS and SIMPLE BRIDGES 4 5 6 2 7 3 8 Usage: … … 33 38 log_filename=None): 34 39 35 """Boyd's generalisation of the US department of transportation culvert 36 model 37 38 The quantity of flow passing through a culvert is controlled by many factors 39 It could be that the culvert is controlled by the inlet only, with it 40 being unsubmerged this is effectively equivalent to the weir Equation 41 Else the culvert could be controlled by the inlet, with it being 42 submerged, this is effectively the Orifice Equation 43 Else it may be controlled by down stream conditions where depending on 44 the down stream depth, the momentum in the culvert etc. flow is controlled 40 """Boyd's generalisation of the US department of transportation culvert methods 41 42 WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER 43 FULL PIPE OR CRITICAL DEPTH ONLY 44 For Supercritical flow this is UNDERESTIMATING the Outlet Velocity 45 The obtain the CORRECT velocity requires an iteration of Depth to Establish 46 the Normal Depth of flow in the pipe. 47 48 It is proposed to provide this in a seperate routine called 49 boyd_generalised_culvert_model_complex 50 51 The Boyd Method is based on methods described by the following: 52 1. 53 US Dept. Transportation Federal Highway Administration (1965) 54 Hydraulic Chart for Selection of Highway Culverts. 55 Hydraulic Engineering Circular No. 5 US Government Printing 56 2. 57 US Dept. Transportation Federal Highway Administration (1972) 58 Capacity charts for the Hydraulic design of highway culverts. 59 Hydraulic Engineering Circular No. 10 US Government Printing 60 These documents provide around 60 charts for various configurations of culverts and inlets. 61 62 Note these documents have been superceded by: 63 2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS5), 64 Which combines culvert design information previously contained in Hydraulic Engineering Circulars 65 (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information. 66 HEC5 provides 20 Charts 67 HEC10 Provides an additional 36 Charts 68 HEC13 Discusses the Design of improved more efficient inlets 69 HDS5 Provides 60 sets of Charts 70 71 In 1985 Professor Michael Boyd Published "HeadDischarge Relations for Culverts", and in 72 1987 published "Generalised Head Discharge Equations for Culverts". 73 These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations. 74 75 It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done. 76 The additional charts cover a range of culvert shapes and inlet configurations 77 78 45 79 """ 46 80 … … 48 82 from anuga.config import velocity_protection 49 83 from anuga.utilities.numerical_tools import safe_acos as acos 50 51 84 print "STARTING..BOYD................." 85 local_debug ='false' 52 86 if inlet_depth > 0.1: #this value was 0.01: 87 if local_debug =='true': 88 print 'Specific E & Deltat Tot E = ',inlet_specific_energy,delta_total_energy 89 print 'culvert typ = ',culvert_type 53 90 # Water has risen above inlet 54 91 … … 61 98 62 99 63 if culvert_type == 'circle': 64 # Round culvert (use width as diameter) 65 diameter = culvert_width 66 100 if culvert_type =='circle': 101 # Round culvert (use height as diameter) 102 diameter = culvert_height 103 104 """ 105 For a circular pipe the Boyd method reviews 3 conditions 106 1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet) 107 2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice) 108 3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe 109 110 For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe 111 """ 112 67 113 # Calculate flows for inlet control 68 114 Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged 69 115 Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63 # Inlet Ctrl Inlet Submerged 116 # Note for to SUBMERGED TO OCCUR inlet_specific_energy should be > 1.2 x diameter.... Should Check !!! 70 117 71 118 if log_filename is not None: 72 119 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged) 73 120 log_to_file(log_filename, s) 74 75 76 # FIXME(Ole): Are these functions really for inlet control? 77 if Q_inlet_unsubmerged < Q_inlet_submerged: 78 Q = Q_inlet_unsubmerged 79 alpha = acos(1  inlet_depth/diameter) 80 flow_area = diameter**2 * (alpha  sin(alpha)*cos(alpha)) 81 outlet_culvert_depth = inlet_depth 82 case = 'Inlet unsubmerged' 83 else: 84 Q = Q_inlet_submerged 85 flow_area = (diameter/2)**2 * pi 86 outlet_culvert_depth = diameter 87 case = 'Inlet submerged' 88 89 90 91 if delta_total_energy < inlet_specific_energy: 121 Q = min(Q_inlet_unsubmerged, Q_inlet_submerged) 122 123 # THE LOWEST Value will Control Calcs From here 124 # Calculate Critical Depth Based on the Adopted Flow as an Estimate 125 dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75) 126 dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95) 127 # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as 128 if dcrit1/diameter > 0.85: 129 outlet_culvert_depth = dcrit2 130 else: 131 outlet_culvert_depth = dcrit1 132 #outlet_culvert_depth = min(outlet_culvert_depth, diameter) 133 # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter 134 if outlet_culvert_depth >= diameter: 135 outlet_culvert_depth = diameter # Once again the pipe is flowing full not partfull 136 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert 137 perimeter = diameter * pi 138 flow_width= diameter 139 case = 'Inlet CTRL Outlet submerged Circular PIPE FULL' 140 if local_debug =='true': 141 print 'Inlet CTRL Outlet submerged Circular PIPE FULL' 142 else: 143 #alpha = acos(1  outlet_culvert_depth/diameter) # Where did this Come From ????/ 144 alpha = acos(12*outlet_culvert_depth/diameter)*2 145 #flow_area = diameter**2 * (alpha  sin(alpha)*cos(alpha)) # Pipe is Running Partly Full at the INLET WHRE did this Come From ????? 146 flow_area = diameter**2/8*(alpha  sin(alpha)) # Equation from GIECK 5th Ed. Pg. B3 147 flow_width= diameter*sin(alpha/2.0) 148 perimeter = alpha*diameter/2.0 149 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth' 150 if local_debug =='true': 151 print 'INLET CTRL Culvert is open channel flow we will for now assume critical depth' 152 print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,' ',alpha 153 if delta_total_energy < inlet_specific_energy: # OUTLET CONTROL !!!! 92 154 # Calculate flows for outlet control 93 155 94 156 # Determine the depth at the outlet relative to the depth of flow in the Culvert 95 if outlet_depth > diameter: # The Outlet is Submerged157 if outlet_depth > diameter: # Outlet is submerged Assume the end of the Pipe is flowing FULL 96 158 outlet_culvert_depth=diameter 97 159 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert 98 160 perimeter = diameter * pi 161 flow_width= diameter 99 162 case = 'Outlet submerged' 100 elif outlet_depth==0.0: 101 outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth 102 alpha = acos(1  inlet_depth/diameter) 103 flow_area = diameter**2 * (alpha  sin(alpha)*cos(alpha)) 104 perimeter = alpha*diameter 105 106 case = 'Outlet depth is zero' 107 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 108 outlet_culvert_depth=outlet_depth # For purpose of calculation assume the outlet depth = the inlet depth 109 alpha = acos(1  outlet_depth/diameter) 110 flow_area = diameter**2 * (alpha  sin(alpha)*cos(alpha)) 111 perimeter = alpha*diameter 112 case = 'Outlet is open channel flow' 113 114 hyd_rad = flow_area/perimeter 115 116 if log_filename is not None: 117 s = 'hydraulic radius at outlet = %f' %hyd_rad 118 log_to_file(log_filename, s) 119 120 # Outlet control velocity using tail water 121 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 122 Q_outlet_tailwater = flow_area * culvert_velocity 123 124 if log_filename is not None: 125 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater 126 log_to_file(log_filename, s) 127 128 Q = min(Q, Q_outlet_tailwater) 129 else: 130 pass 163 if local_debug =='true': 164 print 'Outlet submerged' 165 else: # Culvert running PART FULL for PART OF ITS LENGTH Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 166 # IF outlet_depth < diameter 167 dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75) 168 dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95) 169 if dcrit1/diameter >0.85: 170 outlet_culvert_depth= dcrit2 171 else: 172 outlet_culvert_depth = dcrit1 173 if outlet_culvert_depth > diameter: 174 outlet_culvert_depth = diameter # Once again the pipe is flowing full not partfull 175 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert 176 perimeter = diameter * pi 177 flow_width= diameter 178 case = 'Outlet unsubmerged PIPE FULL' 179 if local_debug =='true': 180 print 'Outlet unsubmerged PIPE FULL' 181 else: 182 alpha = acos(12*outlet_culvert_depth/diameter)*2 183 flow_area = diameter**2/8*(alpha  sin(alpha)) # Equation from GIECK 5th Ed. Pg. B3 184 flow_width= diameter*sin(alpha/2.0) 185 perimeter = alpha*diameter/2.0 186 case = 'Outlet is open channel flow we will for now assume critical depth' 187 if local_debug =='true': 188 print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,' ',alpha 189 print 'Outlet is open channel flow we will for now assume critical depth' 190 if local_debug =='true': 191 print 'FLOW AREA = ',flow_area 192 print 'PERIMETER = ',perimeter 193 print 'Q Interim = ',Q 194 hyd_rad = flow_area/perimeter 195 196 if log_filename is not None: 197 s = 'hydraulic radius at outlet = %f' %hyd_rad 198 log_to_file(log_filename, s) 199 200 # Outlet control velocity using tail water 201 if local_debug =='true': 202 print 'GOT IT ALL CALCULATING Velocity' 203 print 'HydRad = ',hyd_rad 204 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 205 Q_outlet_tailwater = flow_area * culvert_velocity 206 if local_debug =='true': 207 print 'VELOCITY = ',culvert_velocity 208 print 'Outlet Ctrl Q = ',Q_outlet_tailwater 209 if log_filename is not None: 210 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater 211 log_to_file(log_filename, s) 212 Q = min(Q, Q_outlet_tailwater) 213 if local_debug =='true': 214 print ('%s,%.3f,%.3f' %('dcrit 1 , dcit2 =',dcrit1,dcrit2)) 215 print ('%s,%.3f,%.3f,%.3f' %('Q and Velocity and Depth=',Q,culvert_velocity,outlet_culvert_depth)) 216 217 else: 218 pass 131 219 #FIXME(Ole): What about inlet control? 132 133 134 else: 135 # Box culvert (rectangle or square) 220 # ==== END OF CODE BLOCK FOR "IF" CIRCULAR PIPE 221 222 #else.... 223 if culvert_type =='box': 224 if local_debug =='true': 225 print 'BOX CULVERT' 226 # Box culvert (rectangle or square) ======================================================================================================================== 136 227 137 228 # Calculate flows for inlet control 138 229 height = culvert_height 139 width = culvert_width 230 width = culvert_width 231 flow_width=culvert_width 140 232 141 233 Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged … … 146 238 log_to_file(log_filename, s) 147 239 148 149 240 # FIXME(Ole): Are these functions really for inlet control? 150 241 if Q_inlet_unsubmerged < Q_inlet_submerged: 151 242 Q = Q_inlet_unsubmerged 152 flow_area = width*inlet_depth 153 outlet_culvert_depth = inlet_depth 154 case = 'Inlet unsubmerged' 243 dcrit = (Q**2/g/width**2)**0.333333 244 if dcrit > height: 245 dcrit = height 246 flow_area = width*dcrit 247 outlet_culvert_depth = dcrit 248 case = 'Inlet unsubmerged Box Acts as Weir' 155 249 else: 156 250 Q = Q_inlet_submerged 157 251 flow_area = width*height 158 252 outlet_culvert_depth = height 159 case = 'Inlet submerged' 160 253 case = 'Inlet submerged Box Acts as Orifice' 254 255 dcrit = (Q**2/g/width**2)**0.333333 256 257 outlet_culvert_depth = dcrit 258 if outlet_culvert_depth > height: 259 outlet_culvert_depth = height # Once again the pipe is flowing full not partfull 260 flow_area = width*height # Cross sectional area of flow in the culvert 261 perimeter = 2*(width+height) 262 case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL' 263 else: 264 flow_area = width * outlet_culvert_depth 265 perimeter = width+2*outlet_culvert_depth 266 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth' 267 161 268 if delta_total_energy < inlet_specific_energy: 162 269 # Calculate flows for outlet control … … 168 275 perimeter=2.0*(width+height) 169 276 case = 'Outlet submerged' 170 elif outlet_depth==0.0:171 outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth172 flow_area=width*inlet_depth173 perimeter=(width+2.0*inlet_depth)174 case = 'Outlet depth is zero'175 277 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 176 outlet_culvert_depth=outlet_depth 177 flow_area=width*outlet_depth 178 perimeter=(width+2.0*outlet_depth) 179 case = 'Outlet is open channel flow' 278 dcrit = (Q**2/g/width**2)**0.333333 279 outlet_culvert_depth=dcrit # For purpose of calculation assume the outlet depth = Critical Depth 280 if outlet_culvert_depth > height: 281 outlet_culvert_depth=height 282 flow_area=width*height 283 perimeter=2.0*(width+height) 284 case = 'Outlet is Flowing Full' 285 else: 286 flow_area=width*outlet_culvert_depth 287 perimeter=(width+2.0*outlet_culvert_depth) 288 case = 'Outlet is open channel flow' 180 289 181 290 hyd_rad = flow_area/perimeter 182 291 183 292 if log_filename is not None: 184 293 s = 'hydraulic radius at outlet = %f' % hyd_rad … … 196 305 pass 197 306 #FIXME(Ole): What about inlet control? 307 # ==== END OF CODE BLOCK FOR "IF" BOX 198 308 199 309 200 # Common code for circle and square geometries310 # Common code for circle and box geometries ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 201 311 if log_filename is not None: 202 312 log_to_file(log_filename, 'Case: "%s"' % case) … … 205 315 s = 'Flow Rate Control = %f' % Q 206 316 log_to_file(log_filename, s) 207 208 209 culv_froude=sqrt(Q**2*width/(g*flow_area**3)) 210 317 culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3)) 318 if local_debug =='true': 319 print 'FLOW AREA = ',flow_area 320 print 'PERIMETER = ',perimeter 321 print 'Q final = ',Q 322 print 'FROUDE = ',culv_froude 211 323 if log_filename is not None: 212 324 s = 'Froude in Culvert = %f' % culv_froude … … 215 327 # Determine momentum at the outlet 216 328 barrel_velocity = Q/(flow_area + velocity_protection/flow_area) 217 329 # END CODE BLOCK for DEPTH > Required depth for CULVERT Flow 218 330 219 331 else: # inlet_depth < 0.01:
Note: See TracChangeset
for help on using the changeset viewer.