Changeset 5434
 Timestamp:
 Jun 25, 2008, 5:01:54 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/culvert_flows/culvert_flow.py
r5433 r5434 1 2 class Culvert_flow: 3 """Culvert flow  transfer water from one hole to another 4 5 Using Momentum as Calculated by Culvert Flow !! 6 Could be Several Methods Investigated to do This !!! 7 8 2008_May_08 9 To Ole: 10 OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on 11 the two input Points (X0,Y0) and (X1,Y1)  need to be passed to create polygon 12 13 The two centers are now passed on to create_polygon. 14 15 16 Input: Two points, pipe_size (either diameter or width, height), mannings_rougness, 17 inlet/outlet energy_loss_coefficients, internal_bend_coefficent, 18 topdown_blockage_factor and bottom_up_blockage_factor 19 20 21 And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront 22 of the culvert 23 At the moment this script uses only Depth, later we can change it to Energy... 24 25 Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity 26 The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet... 27 28 """ 29 30 def __init__(self, 31 domain, 32 label=None, 33 description=None, 34 end_point0=None, 35 end_point1=None, 36 width=None, 37 height=None, 38 diameter=None, 39 manning=None, # Mannings Roughness for Culvert 40 invert_level0=None, # Invert level if not the same as the Elevation on the Domain 41 invert_level1=None, # Invert level if not the same as the Elevation on the Domain 42 loss_exit=None, 43 loss_entry=None, 44 loss_bend=None, 45 loss_special=None, 46 blockage_topdwn=None, 47 blockage_bottup=None, 48 culvert_routine=None, 49 verbose=False): 50 51 from Numeric import sqrt, sum 52 53 # Input check 54 if diameter is not None: 55 self.culvert_type = 'circle' 56 self.diameter = diameter 57 if height is not None or width is not None: 58 msg = 'Either diameter or width&height must be specified, but not both.' 59 raise Exception, msg 60 else: 61 if height is not None: 62 if width is None: 63 self.culvert_type = 'square' 64 width = height 65 else: 66 self.culvert_type = 'rectangle' 67 elif width is not None: 68 if height is None: 69 self.culvert_type = 'square' 70 height = width 71 else: 72 msg = 'Either diameter or width&height must be specified.' 73 raise Exception, msg 74 75 if height == width: 76 self.culvert_type = 'square' 77 78 self.height = height 79 self.width = width 80 81 82 assert self.culvert_type in ['circle', 'square', 'rectangle'] 83 84 # Set defaults 85 if manning is None: manning = 0.012 # Set a Default Mannings Roughness for Pipe 86 if loss_exit is None: loss_exit = 1.00 87 if loss_entry is None: loss_entry = 0.50 88 if loss_bend is None: loss_bend=0.00 89 if loss_special is None: loss_special=0.00 90 if blockage_topdwn is None: blockage_topdwn=0.00 91 if blockage_bottup is None: blockage_bottup=0.00 92 if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model 93 if label is None: label = 'culvert_flow_' + id(self) 94 95 # Open log file for storing some specific results... 96 self.log_filename = label + '.log' 97 self.label = label 98 99 100 # Create the fundamental culvert polygons from POLYGON 101 P = create_culvert_polygons(end_point0, 102 end_point1, 103 width=width, 104 height=height) 105 106 if verbose is True: 107 pass 108 #plot_polygons([[end_point0, end_point1], 109 # P['exchange_polygon0'], 110 # P['exchange_polygon1'], 111 # P['enquiry_polygon0'], 112 # P['enquiry_polygon1']], 113 # figname='culvert_polygon_output') 114 115 self.openings = [] 116 self.openings.append(Inflow(domain, 117 polygon=P['exchange_polygon0'])) 118 119 self.openings.append(Inflow(domain, 120 polygon=P['exchange_polygon1'])) 121 122 123 # Assume two openings for now: Referred to as 0 and 1 124 assert len(self.openings) == 2 125 126 # Store basic geometry 127 self.end_points = [end_point0, end_point1] 128 self.invert_levels = [invert_level0, invert_level1] 129 self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']] 130 self.vector = P['vector'] 131 self.length = P['length']; assert self.length > 0.0 132 self.verbose = verbose 133 self.last_time = 0.0 134 self.temp_keep_delta_t = 0.0 135 136 137 # Store hydraulic parameters 138 self.manning = manning 139 self.loss_exit = loss_exit 140 self.loss_entry = loss_entry 141 self.loss_bend = loss_bend 142 self.loss_special = loss_special 143 self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special 144 self.blockage_topdwn = blockage_topdwn 145 self.blockage_bottup = blockage_bottup 146 147 # Store culvert routine 148 self.culvert_routine = culvert_routine 149 150 151 # Create objects to update momentum (a bit crude at this stage) 152 153 154 xmom0 = General_forcing(domain, 'xmomentum', 155 polygon=P['exchange_polygon0']) 156 157 xmom1 = General_forcing(domain, 'xmomentum', 158 polygon=P['exchange_polygon1']) 159 160 ymom0 = General_forcing(domain, 'ymomentum', 161 polygon=P['exchange_polygon0']) 162 163 ymom1 = General_forcing(domain, 'ymomentum', 164 polygon=P['exchange_polygon1']) 165 166 self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ] 167 168 169 # Print something 170 s = 'Culvert Effective Length = %.2f m' %(self.length) 171 log_to_file(self.log_filename, s) 172 173 s = 'Culvert Direction is %s\n' %str(self.vector) 174 log_to_file(self.log_filename, s) 175 176 def __call__(self, domain): 177 from anuga.utilities.numerical_tools import mean 178 from anuga.utilities.polygon import inside_polygon 179 from anuga.config import g, epsilon 180 from Numeric import take, sqrt 181 from anuga.config import velocity_protection 182 183 184 log_filename = self.log_filename 185 186 # Time stuff 187 time = domain.get_time() 188 delta_t = timeself.last_time 189 s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) 190 log_to_file(log_filename, s) 191 192 msg = 'Time did not advance' 193 if time > 0.0: assert delta_t > 0.0, msg 194 195 196 # Get average water depths at each opening 197 openings = self.openings # There are two Opening [0] and [1] 198 for i, opening in enumerate(openings): 199 stage = domain.quantities['stage'].get_values(location='centroids', 200 indices=opening.exchange_indices) 201 elevation = domain.quantities['elevation'].get_values(location='centroids', 202 indices=opening.exchange_indices) 203 204 # Indices corresponding to energy enquiry field for this opening 205 coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y) 206 enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 207 208 209 # Get model values for points in enquiry polygon for this opening 210 dq = domain.quantities 211 stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices) 212 xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices) 213 ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices) 214 elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices) 215 depth = stage  elevation 216 217 # Compute mean values of selected quantitites in the enquiry area in front of the culvert 218 # Epsilon handles a dry cell case 219 ux = xmomentum/(depth+velocity_protection/depth) # Velocity (xdirection) 220 uy = ymomentum/(depth+velocity_protection/depth) # Velocity (ydirection) 221 v = mean(sqrt(ux**2+uy**2)) # Average velocity 222 w = mean(stage) # Average stage 223 224 # Store values at enquiry field 225 opening.velocity = v 226 227 228 # Compute mean values of selected quantitites in the exchange area in front of the culvert 229 # Stage and velocity comes from enquiry area and elevation from exchange area 230 231 # Use invert level instead of elevation if specified 232 invert_level = self.invert_levels[i] 233 if invert_level is not None: 234 z = invert_level 235 else: 236 elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices) 237 z = mean(elevation) # Average elevation 238 239 # Estimated depth above the culvert inlet 240 d = w  z 241 242 if d < 0.0: 243 # This is possible since w and z are taken at different locations 244 #msg = 'D < 0.0: %f' %d 245 #raise msg 246 d = 0.0 247 248 # Ratio of depth to culvert height. 249 # If ratio > 1 then culvert is running full 250 ratio = d/self.height 251 opening.ratio = ratio 252 253 # Average measures of energy in front of this opening 254 Es = d + 0.5*v**2/g # Specific energy in exchange area 255 Et = w + 0.5*v**2/g # Total energy in the enquiry area 256 opening.total_energy = Et 257 opening.specific_energy = Es 258 259 # Store current average stage and depth with each opening object 260 opening.depth = d 261 opening.stage = w 262 opening.elevation = z 263 264 265 ################# End of the FOR loop ################################################ 266 267 268 # We now need to deal with each opening individually 269 270 # Determine flow direction based on total energy difference 271 delta_Et = openings[0].total_energy  openings[1].total_energy 272 273 if delta_Et > 0: 274 #print 'Flow U/S > D/S' 275 inlet=openings[0] 276 outlet=openings[1] 277 278 inlet.momentum = self.opening_momentum[0] 279 outlet.momentum = self.opening_momentum[1] 280 281 else: 282 #print 'Flow D/S > U/S' 283 inlet=openings[1] 284 outlet=openings[0] 285 286 inlet.momentum = self.opening_momentum[1] 287 outlet.momentum = self.opening_momentum[0] 288 289 delta_Et = delta_Et 290 291 msg = 'Total energy difference is negative' 292 assert delta_Et > 0.0, msg 293 294 delta_h = inlet.stage  outlet.stage 295 delta_z = inlet.elevation  outlet.elevation 296 culvert_slope = (delta_z/self.length) 297 298 if culvert_slope < 0.0: 299 # Adverse gradient  flow is running uphill 300 # Flow will be purely controlled by uphill outlet face 301 print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation 302 303 304 s = 'Delta total energy = %.3f' %(delta_Et) 305 log_to_file(log_filename, s) 306 307 308 Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g) 309 ##################################################### 310 barrel_momentum = barrel_velocity*culvert_outlet_depth 311 312 s = 'Barrel velocity = %f' %barrel_velocity 313 log_to_file(log_filename, s) 314 315 # Compute momentum vector at outlet 316 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum 317 318 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) 319 log_to_file(log_filename, s) 320 321 delta_t = time  self.last_time 322 if delta_t > 0.0: 323 xmomentum_rate = outlet_mom_x  outlet.momentum[0].value 324 xmomentum_rate /= delta_t 325 326 ymomentum_rate = outlet_mom_y  outlet.momentum[1].value 327 ymomentum_rate /= delta_t 328 329 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate) 330 log_to_file(log_filename, s) 331 else: 332 xmomentum_rate = ymomentum_rate = 0.0 333 334 335 # Set momentum rates for outlet jet 336 outlet.momentum[0].rate = xmomentum_rate 337 outlet.momentum[1].rate = ymomentum_rate 338 339 # Remember this value for next step (IMPORTANT) 340 outlet.momentum[0].value = outlet_mom_x 341 outlet.momentum[1].value = outlet_mom_y 342 343 if int(domain.time*100) % 100 == 0: 344 s = 'T=%.5f, Culvert Discharge = %.3f f'\ 345 %(time, Q) 346 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\ 347 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y) 348 s += ' Momentum rate: (%.4f, %.4f)'\ 349 %(xmomentum_rate, ymomentum_rate) 350 s+='Outlet Vel= %.3f'\ 351 %(barrel_velocity) 352 log_to_file(log_filename, s) 353 354 355 356 357 358 # Execute flow term for each opening 359 # This is where Inflow objects are evaluated and update the domain 360 for opening in self.openings: 361 opening(domain) 362 363 # Execute momentum terms 364 # This is where Inflow objects are evaluated and update the domain 365 outlet.momentum[0](domain) 366 outlet.momentum[1](domain) 367 368 # Store value of time 369 self.last_time = time 370 371
Note: See TracChangeset
for help on using the changeset viewer.