- Timestamp:
- Aug 30, 2010, 5:53:45 PM (14 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/structures/culvert_routine.py
r7979 r7980 5 5 from anuga.config import g 6 6 7 class Culvert_routine s:7 class Culvert_routine: 8 8 """Collection of culvert routines for use with Culvert_operator 9 9 … … 23 23 self.inlets = culvert.inlets 24 24 25 self.inflow = self.inlets[0] 26 self.outflow = self.inlets[1] 27 25 28 26 self.culvert_length = culvert.get_culvert_length() 29 27 self.culvert_width = culvert.get_culvert_width() … … 34 32 self.log_filename = None 35 33 36 # Determine flow direction based on total energy difference 37 self.delta_total_energy = self.inflow.get_average_total_energy() - self.outflow.get_average_total_energy() 38 39 if self.delta_total_energy < 0: 40 self.inflow = self.inlets[1] 41 self.outflow = self.inlets[0] 42 self.delta_total_energy = -self.delta_total_energy 34 self.determine_inflow() 43 35 44 36 #delta_z = self.self.inflow.get_average_elevation() - self.self.outflow.get_average_elevation() … … 53 45 # driving_head = self.inflow.get_average_specific_energy() 54 46 47 48 49 def determine_inflow(self): 50 # Determine flow direction based on total energy difference 51 52 self.delta_total_energy = self.inlets[0].get_average_total_energy() - self.inlets[1].get_average_total_energy() 53 54 self.inflow = self.inlets[0] 55 self.outflow = self.inlets[1] 55 56 57 58 if self.delta_total_energy < 0: 59 self.inflow = self.inlets[1] 60 self.outflow = self.inlets[0] 61 self.delta_total_energy = -self.delta_total_energy 62 63 56 64 def get_inflow(self): 57 65 … … 63 71 return self.outflow 64 72 65 66 def boyd_box(self):67 68 """Boyd's generalisation of the US department of transportation culvert methods69 70 WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER71 FULL PIPE OR CRITICAL DEPTH ONLY72 For Supercritical flow this is UNDERESTIMATING the Outlet Velocity73 The obtain the CORRECT velocity requires an iteration of Depth to Establish74 the Normal Depth of flow in the pipe.75 76 It is proposed to provide this in a seperate routine called77 boyd_generalised_culvert_model_complex78 79 The Boyd Method is based on methods described by the following:80 1.81 US Dept. Transportation Federal Highway Administration (1965)82 Hydraulic Chart for Selection of Highway Culverts.83 Hydraulic Engineering Circular No. 5 US Government Printing84 2.85 US Dept. Transportation Federal Highway Administration (1972)86 Capacity charts for the Hydraulic design of highway culverts.87 Hydraulic Engineering Circular No. 10 US Government Printing88 These documents provide around 60 charts for various configurations of culverts and inlets.89 90 Note these documents have been superceded by:91 2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),92 Which combines culvert design information previously contained in Hydraulic Engineering Circulars93 (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.94 HEC-5 provides 20 Charts95 HEC-10 Provides an additional 36 Charts96 HEC-13 Discusses the Design of improved more efficient inlets97 HDS-5 Provides 60 sets of Charts98 99 In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in100 1987 published "Generalised Head Discharge Equations for Culverts".101 These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.102 103 It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.104 The additional charts cover a range of culvert shapes and inlet configurations105 106 107 """108 109 local_debug ='false'110 if self.inflow.get_average_height() > 0.1: #this value was 0.01:111 if local_debug =='true':112 log.critical('Specific E & Deltat Tot E = %s, %s'113 % (str(self.inflow.get_average_specific_energy()),114 str(self.delta_total_energy)))115 log.critical('culvert type = %s' % str(culvert_type))116 # Water has risen above inlet117 118 if self.log_filename is not None:119 s = 'Specific energy = %f m' % self.inflow.get_average_specific_energy()120 log_to_file(self.log_filename, s)121 122 msg = 'Specific energy at inlet is negative'123 assert self.inflow.get_average_specific_energy() >= 0.0, msg124 125 height = self.culvert_height126 width = self.culvert_width127 flow_width = self.culvert_width128 129 Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged130 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 Submerged131 132 # FIXME(Ole): Are these functions really for inlet control?133 if Q_inlet_unsubmerged < Q_inlet_submerged:134 Q = Q_inlet_unsubmerged135 dcrit = (Q**2/g/width**2)**0.333333136 if dcrit > height:137 dcrit = height138 flow_area = width*dcrit139 outlet_culvert_depth = dcrit140 case = 'Inlet unsubmerged Box Acts as Weir'141 else:142 Q = Q_inlet_submerged143 flow_area = width*height144 outlet_culvert_depth = height145 case = 'Inlet submerged Box Acts as Orifice'146 147 dcrit = (Q**2/g/width**2)**0.333333148 149 outlet_culvert_depth = dcrit150 if outlet_culvert_depth > height:151 outlet_culvert_depth = height # Once again the pipe is flowing full not partfull152 flow_area = width*height # Cross sectional area of flow in the culvert153 perimeter = 2*(width+height)154 case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'155 else:156 flow_area = width * outlet_culvert_depth157 perimeter = width+2*outlet_culvert_depth158 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'159 160 if self.delta_total_energy < self.inflow.get_average_specific_energy():161 # Calculate flows for outlet control162 163 # Determine the depth at the outlet relative to the depth of flow in the Culvert164 if self.outflow.get_average_height() > height: # The Outlet is Submerged165 outlet_culvert_depth=height166 flow_area=width*height # Cross sectional area of flow in the culvert167 perimeter=2.0*(width+height)168 case = 'Outlet submerged'169 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity170 dcrit = (Q**2/g/width**2)**0.333333171 outlet_culvert_depth=dcrit # For purpose of calculation assume the outlet depth = Critical Depth172 if outlet_culvert_depth > height:173 outlet_culvert_depth=height174 flow_area=width*height175 perimeter=2.0*(width+height)176 case = 'Outlet is Flowing Full'177 else:178 flow_area=width*outlet_culvert_depth179 perimeter=(width+2.0*outlet_culvert_depth)180 case = 'Outlet is open channel flow'181 182 hyd_rad = flow_area/perimeter183 184 if self.log_filename is not None:185 s = 'hydraulic radius at outlet = %f' % hyd_rad186 log_to_file(self.log_filename, s)187 188 # Outlet control velocity using tail water189 culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))190 Q_outlet_tailwater = flow_area * culvert_velocity191 192 if self.log_filename is not None:193 s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater194 log_to_file(self.log_filename, s)195 Q = min(Q, Q_outlet_tailwater)196 else:197 pass198 #FIXME(Ole): What about inlet control?199 200 culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))201 if local_debug =='true':202 log.critical('FLOW AREA = %s' % str(flow_area))203 log.critical('PERIMETER = %s' % str(perimeter))204 log.critical('Q final = %s' % str(Q))205 log.critical('FROUDE = %s' % str(culv_froude))206 207 # Determine momentum at the outlet208 barrel_velocity = Q/(flow_area + velocity_protection/flow_area)209 210 # END CODE BLOCK for DEPTH > Required depth for CULVERT Flow211 212 else: # self.inflow.get_average_height() < 0.01:213 Q = barrel_velocity = outlet_culvert_depth = 0.0214 215 # Temporary flow limit216 if barrel_velocity > self.max_velocity:217 barrel_velocity = self.max_velocity218 Q = flow_area * barrel_velocity219 220 return Q, barrel_velocity, outlet_culvert_depth221 73 222 74 223 def boyd_circle(self):224 """225 For a circular pipe the Boyd method reviews 3 conditions226 1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)227 2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)228 3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe229 230 For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe231 """232 233 diameter = self.culvert_height234 235 local_debug ='false'236 if self.inflow.get_average_height() > 0.1: #this value was 0.01:237 if local_debug =='true':238 log.critical('Specific E & Deltat Tot E = %s, %s'239 % (str(self.inflow.get_average_specific_energy()),240 str(self.delta_total_energy)))241 log.critical('culvert type = %s' % str(culvert_type))242 # Water has risen above inlet243 244 if self.log_filename is not None:245 s = 'Specific energy = %f m' % self.inflow.get_average_specific_energy()246 log_to_file(self.log_filename, s)247 248 msg = 'Specific energy at inlet is negative'249 assert self.inflow.get_average_specific_energy() >= 0.0, msg250 251 # Calculate flows for inlet control252 Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*self.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged253 Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*self.inflow.get_average_specific_energy()**0.63 # Inlet Ctrl Inlet Submerged254 # Note for to SUBMERGED TO OCCUR self.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!255 256 if self.log_filename is not None:257 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)258 log_to_file(self.log_filename, s)259 Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)260 261 # THE LOWEST Value will Control Calcs From here262 # Calculate Critical Depth Based on the Adopted Flow as an Estimate263 dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)264 dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)265 # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as266 if dcrit1/diameter > 0.85:267 outlet_culvert_depth = dcrit2268 else:269 outlet_culvert_depth = dcrit1270 #outlet_culvert_depth = min(outlet_culvert_depth, diameter)271 # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter272 if outlet_culvert_depth >= diameter:273 outlet_culvert_depth = diameter # Once again the pipe is flowing full not partfull274 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert275 perimeter = diameter * pi276 flow_width= diameter277 case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'278 if local_debug == 'true':279 log.critical('Inlet CTRL Outlet submerged Circular '280 'PIPE FULL')281 else:282 #alpha = acos(1 - outlet_culvert_depth/diameter) # Where did this Come From ????/283 alpha = acos(1-2*outlet_culvert_depth/diameter)*2284 #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) # Pipe is Running Partly Full at the INLET WHRE did this Come From ?????285 flow_area = diameter**2/8*(alpha - sin(alpha)) # Equation from GIECK 5th Ed. Pg. B3286 flow_width= diameter*sin(alpha/2.0)287 perimeter = alpha*diameter/2.0288 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'289 if local_debug =='true':290 log.critical('INLET CTRL Culvert is open channel flow '291 'we will for now assume critical depth')292 log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'293 % (str(Q), str(outlet_culvert_depth),294 str(alpha)))295 if self.delta_total_energy < self.inflow.get_average_specific_energy(): # OUTLET CONTROL !!!!296 # Calculate flows for outlet control297 298 # Determine the depth at the outlet relative to the depth of flow in the Culvert299 if self.outflow.get_average_height() > diameter: # Outlet is submerged Assume the end of the Pipe is flowing FULL300 outlet_culvert_depth=diameter301 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert302 perimeter = diameter * pi303 flow_width= diameter304 case = 'Outlet submerged'305 if local_debug =='true':306 log.critical('Outlet submerged')307 else: # Culvert running PART FULL for PART OF ITS LENGTH Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity308 # IF self.outflow.get_average_height() < diameter309 dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)310 dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)311 if dcrit1/diameter >0.85:312 outlet_culvert_depth= dcrit2313 else:314 outlet_culvert_depth = dcrit1315 if outlet_culvert_depth > diameter:316 outlet_culvert_depth = diameter # Once again the pipe is flowing full not partfull317 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert318 perimeter = diameter * pi319 flow_width= diameter320 case = 'Outlet unsubmerged PIPE FULL'321 if local_debug =='true':322 log.critical('Outlet unsubmerged PIPE FULL')323 else:324 alpha = acos(1-2*outlet_culvert_depth/diameter)*2325 flow_area = diameter**2/8*(alpha - sin(alpha)) # Equation from GIECK 5th Ed. Pg. B3326 flow_width= diameter*sin(alpha/2.0)327 perimeter = alpha*diameter/2.0328 case = 'Outlet is open channel flow we will for now assume critical depth'329 if local_debug == 'true':330 log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'331 % (str(Q), str(outlet_culvert_depth),332 str(alpha)))333 log.critical('Outlet is open channel flow we '334 'will for now assume critical depth')335 if local_debug == 'true':336 log.critical('FLOW AREA = %s' % str(flow_area))337 log.critical('PERIMETER = %s' % str(perimeter))338 log.critical('Q Interim = %s' % str(Q))339 hyd_rad = flow_area/perimeter340 341 if self.log_filename is not None:342 s = 'hydraulic radius at outlet = %f' %hyd_rad343 log_to_file(self.log_filename, s)344 345 # Outlet control velocity using tail water346 if local_debug =='true':347 log.critical('GOT IT ALL CALCULATING Velocity')348 log.critical('HydRad = %s' % str(hyd_rad))349 culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))350 Q_outlet_tailwater = flow_area * culvert_velocity351 if local_debug =='true':352 log.critical('VELOCITY = %s' % str(culvert_velocity))353 log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))354 if self.log_filename is not None:355 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater356 log_to_file(self.log_filename, s)357 Q = min(Q, Q_outlet_tailwater)358 if local_debug =='true':359 log.critical('%s,%.3f,%.3f'360 % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))361 log.critical('%s,%.3f,%.3f,%.3f'362 % ('Q and Velocity and Depth=', Q,363 culvert_velocity, outlet_culvert_depth))364 365 culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))366 if local_debug =='true':367 log.critical('FLOW AREA = %s' % str(flow_area))368 log.critical('PERIMETER = %s' % str(perimeter))369 log.critical('Q final = %s' % str(Q))370 log.critical('FROUDE = %s' % str(culv_froude))371 372 # Determine momentum at the outlet373 barrel_velocity = Q/(flow_area + velocity_protection/flow_area)374 375 else: # self.inflow.get_average_height() < 0.01:376 Q = barrel_velocity = outlet_culvert_depth = 0.0377 378 # Temporary flow limit379 if barrel_velocity > self.max_velocity:380 barrel_velocity = self.max_velocity381 Q = flow_area * barrel_velocity382 383 return Q, barrel_velocity, outlet_culvert_depth384 75 385 76 77 78
Note: See TracChangeset
for help on using the changeset viewer.