Changeset 5777
- Timestamp:
- Sep 20, 2008, 12:18:52 AM (16 years ago)
- Location:
- anuga_core/source/anuga/culvert_flows
- Files:
-
- 2 edited
-
Test_Culvert_Flat_Water_Lev.py (modified) (1 diff)
-
culvert_class.py (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py
r5774 r5777 37 37 width = 5. 38 38 39 #dx = dy = 1 # Resolution: Length of subdivisions on both axes39 dx = dy = 1 # Resolution: Length of subdivisions on both axes 40 40 #dx = dy = .5 # Resolution: Length of subdivisions on both axes 41 dx = dy = .5 # Resolution: Length of subdivisions on both axes41 #dx = dy = .5 # Resolution: Length of subdivisions on both axes 42 42 #dx = dy = .1 # Resolution: Length of subdivisions on both axes 43 43 -
anuga_core/source/anuga/culvert_flows/culvert_class.py
r5586 r5777 16 16 2008_May_08 17 17 To Ole: 18 OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on 19 the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon 18 OK so here we need to get the Polygon Creating code to create 19 polygons for the culvert Based on 20 the two input Points (X0,Y0) and (X1,Y1) - need to be passed 21 to create polygon 20 22 21 23 The two centers are now passed on to create_polygon. 22 24 23 25 24 Input: Two points, pipe_size (either diameter or width, height), mannings_rougness, 26 Input: Two points, pipe_size (either diameter or width, height), 27 mannings_rougness, 25 28 inlet/outlet energy_loss_coefficients, internal_bend_coefficent, 26 29 top-down_blockage_factor and bottom_up_blockage_factor 27 30 28 31 29 And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront 30 of the culvert 31 At the moment this script uses only Depth, later we can change it to Energy... 32 33 Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity 32 And the Delta H enquiry should be change from Openings in line 412 33 to the enquiry Polygons infront of the culvert 34 At the moment this script uses only Depth, later we can change it to 35 Energy... 36 37 Once we have Delta H can calculate a Flow Rate and from Flow Rate 38 an Outlet Velocity 34 39 The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet... 35 40 41 Invert levels are optional. If left out they will default to the 42 elevation at the opening. 43 36 44 """ 37 45 … … 46 54 diameter=None, 47 55 manning=None, # Mannings Roughness for Culvert 48 invert_level0=None, # Invert level if not the same as the Elevation on the Domain49 invert_level1=None, # Invert level if not the same as the Elevation on the Domain56 invert_level0=None, # Invert level at opening 0 57 invert_level1=None, # Invert level at opening 1 50 58 loss_exit=None, 51 59 loss_entry=None, … … 65 73 self.diameter = diameter 66 74 if height is not None or width is not None: 67 msg = 'Either diameter or width&height must be specified, but not both.' 75 msg = 'Either diameter or width&height must be specified, ' 76 msg += 'but not both.' 68 77 raise Exception, msg 69 78 else: … … 96 105 97 106 # Set defaults 98 if manning is None: manning = 0.012 # Set a Default Mannings Roughness for Pipe107 if manning is None: manning = 0.012 # Default roughness for pipe 99 108 if loss_exit is None: loss_exit = 1.00 100 109 if loss_entry is None: loss_entry = 0.50 … … 103 112 if blockage_topdwn is None: blockage_topdwn=0.00 104 113 if blockage_bottup is None: blockage_bottup=0.00 105 if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model 114 if culvert_routine is None: 115 culvert_routine=boyd_generalised_culvert_model 116 106 117 if label is None: label = 'culvert_flow' 107 118 label += '_' + str(id(self)) 108 109 119 self.label = label 110 120 … … 174 184 self.end_points = [end_point0, end_point1] 175 185 self.invert_levels = [invert_level0, invert_level1] 176 self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']] 186 #self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']] 187 self.enquiry_polylines = [P['enquiry_polygon0'][:2], 188 P['enquiry_polygon1'][:2]] 177 189 self.vector = P['vector'] 178 190 self.length = P['length']; assert self.length > 0.0 … … 245 257 dq = domain.quantities 246 258 247 stage = dq['stage'].get_values(location='centroids',248 indices=opening.exchange_indices)249 elevation = dq['elevation'].get_values(location='centroids',250 indices=opening.exchange_indices)251 252 # Indices corresponding to energy enquiry field for this opening253 coordinates = domain.get_centroid_coordinates(absolute=True) # Get all centroid points (x,y)254 enquiry_indices = inside_polygon(coordinates,255 self.enquiry_polygons[i])256 257 if len(enquiry_indices) == 0:258 msg = 'No triangles have been identified in specified region: %s' %str(self.enquiry_polygons[i])259 raise Exception, msg260 261 # Get model values for points in enquiry polygon for this opening262 stage = dq['stage'].get_values(location='centroids',263 indices=enquiry_indices)264 xmomentum = dq['xmomentum'].get_values(location='centroids',265 indices=enquiry_indices)266 ymomentum = dq['ymomentum'].get_values(location='centroids',267 indices=enquiry_indices)268 elevation = dq['elevation'].get_values(location='centroids',269 indices=enquiry_indices)270 depth = stage - elevation271 272 # Compute mean values of selected quantitites in the enquiry area in front of the culvert273 # Epsilon handles a dry cell case274 ux = xmomentum/(depth+velocity_protection/depth) # Velocity (x-direction)275 uy = ymomentum/(depth+velocity_protection/depth) # Velocity (y-direction)276 #print 'Velocity in culvert:', ux, uy, depth, xmomentum, ymomentum277 v = mean(sqrt(ux**2+uy**2)) # Average velocity278 w = mean(stage) # Average stage279 280 # Store values at enquiry field281 opening.velocity = v282 283 284 259 # Compute mean values of selected quantitites in the exchange area in front of the culvert 285 260 # Stage and velocity comes from enquiry area and elevation from exchange area 286 261 262 stage = dq['stage'].get_values(location='centroids', 263 indices=opening.exchange_indices) 264 w = mean(stage) # Average stage 265 287 266 # Use invert level instead of elevation if specified 288 267 invert_level = self.invert_levels[i] … … 291 270 else: 292 271 elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices) 293 z = mean(elevation) 272 z = mean(elevation) # Average elevation 294 273 295 274 # Estimated depth above the culvert inlet … … 302 281 303 282 304 305 # Depth at exchange area used to trigger calculations306 stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)307 elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)308 depth = stage - elevation309 d_trigger = mean(depth)310 311 312 313 283 # Ratio of depth to culvert height. 314 284 # If ratio > 1 then culvert is running full … … 319 289 opening.ratio = ratio 320 290 291 321 292 # Average measures of energy in front of this opening 322 Es = d + 0.5*v**2/g # Specific energy in exchange area 323 Et = w + 0.5*v**2/g # Total energy in the enquiry area 324 opening.total_energy = Et 325 opening.specific_energy = Es 326 293 polyline = self.enquiry_polylines[i] 294 #print 't = %.4f, opening=%d,' %(domain.time, i), 295 opening.total_energy = domain.get_energy_through_cross_section(polyline, 296 kind='total') 297 #print 'Et = %.3f m' %opening.total_energy 298 327 299 # Store current average stage and depth with each opening object 328 300 opening.depth = d 329 opening.depth_trigger = d _trigger301 opening.depth_trigger = d # Use this for now 330 302 opening.stage = w 331 303 opening.elevation = z … … 342 314 if delta_Et > 0: 343 315 #print 'Flow U/S ---> D/S' 344 inlet =openings[0]345 outlet =openings[1]316 inlet = openings[0] 317 outlet = openings[1] 346 318 347 319 inlet.momentum = self.opening_momentum[0] … … 350 322 else: 351 323 #print 'Flow D/S ---> U/S' 352 inlet =openings[1]353 outlet =openings[0]324 inlet = openings[1] 325 outlet = openings[0] 354 326 355 327 inlet.momentum = self.opening_momentum[1]
Note: See TracChangeset
for help on using the changeset viewer.
