Changeset 5428
- Timestamp:
- Jun 25, 2008, 3:32:05 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/culvert_flow/culvert_boyd_channel.py
r5427 r5428 60 60 #------------------------------------------------------------------------------ 61 61 62 63 def log(fid, s): 64 print s 65 fid.write(s + '\n') 66 62 67 63 68 # Open file for storing some specific results... … … 81 86 domain.H0 = 0.01 82 87 domain.tight_slope_limiters = True 83 domain.set_minimum_storable_height(0.0 2)88 domain.set_minimum_storable_height(0.001) 84 89 85 90 s='Size'+str(len(domain)) 86 fid.write(s) 87 print s 91 log(fid, s) 92 93 velocity_protection = 1.0e-4 88 94 89 95 … … 171 177 # This SHould be changed to a Vertical Opening Both BOX and Circular 172 178 # There will be several Culvert Routines such as: 179 # CULVERT_Boyd_Channel 180 # CULVERT_Orifice_and_Weir 173 181 # CULVERT_Simple_FLOOR 174 182 # CULVERT_Simple_WALL … … 221 229 height=None, 222 230 manning=None, # Mannings Roughness for Culvert 223 inv _point0=None,# Invert level if not the same as the Elevation on the Domain224 inv _point1=None,# Invert level if not the same as the Elevation on the Domain231 invert_level0=None, # Invert level if not the same as the Elevation on the Domain 232 invert_level1=None, # Invert level if not the same as the Elevation on the Domain 225 233 loss_exit=None, 226 234 loss_entry=None, … … 233 241 from Numeric import sqrt, sum 234 242 243 # Input check 235 244 if height is None: height = width 245 246 # Set defaults 236 247 if manning is None: manning = 0.012 # Set a Default Mannings Roughness for Pipe 237 248 if loss_exit is None: loss_exit = 1.00 … … 241 252 if blockage_topdwn is None: blockage_topdwn=0.00 242 253 if blockage_bottup is None: blockage_bottup=0.00 243 #sum_loss=loss_exit+loss_entry+loss_bend+loss_special # Note if losses are a function can't do this !!!!254 244 255 245 256 # Create the fundamental culvert polygons from POLYGON … … 270 281 271 282 # Store basic geometry 272 self.end_points = [end_point0, end_point1] 283 self.end_points = [end_point0, end_point1] 284 self.invert_levels = [invert_level0, invert_level1] 273 285 self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']] 274 286 self.vector = P['vector'] 275 self.distance = P['length'] 287 self.distance = P['length']; assert self.distance > 0.0 276 288 self.verbose = verbose 277 289 self.width = width 278 290 self.height = height 279 291 self.last_time = 0.0 280 self.temp_keep_delta_t=0.0 292 self.temp_keep_delta_t = 0.0 293 294 295 # Store hydraulic parameters 296 self.manning = manning 297 self.loss_exit = loss_exit 298 self.loss_entry = loss_entry 299 self.loss_bend = loss_bend 300 self.loss_special = loss_special 301 self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special 302 self.blockage_topdwn = blockage_topdwn 303 self.blockage_bottup = blockage_bottup 304 305 281 306 # Create objects to update momentum (a bit crude at this stage) 282 self.xmom_forcing0 = General_forcing(domain, 283 'xmomentum', 284 polygon=P['exchange_polygon0']) 285 286 self.xmom_forcing1 = General_forcing(domain, 287 'xmomentum', 288 polygon=P['exchange_polygon1']) 289 290 self.ymom_forcing0 = General_forcing(domain, 291 'ymomentum', 292 polygon=P['exchange_polygon0']) 293 294 self.ymom_forcing1 = General_forcing(domain, 295 'ymomentum', 296 polygon=P['exchange_polygon1']) 307 308 309 xmom0 = General_forcing(domain, 'xmomentum', 310 polygon=P['exchange_polygon0']) 311 312 xmom1 = General_forcing(domain, 'xmomentum', 313 polygon=P['exchange_polygon1']) 314 315 ymom0 = General_forcing(domain, 'ymomentum', 316 polygon=P['exchange_polygon0']) 317 318 ymom1 = General_forcing(domain, 'ymomentum', 319 polygon=P['exchange_polygon1']) 320 321 self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ] 322 297 323 298 324 # Print something 299 300 s= 'Culvert Effective Length ='+str( self.distance) 301 print s 302 fid.write(s) 325 s = 'Culvert Effective Length = %.2f m' %(self.distance) 326 log(fid, s) 303 327 304 s= 'Culvert Direction is '+str( self.vector) 305 print s 306 fid.write(s) 307 s= 'Point 1m Up Stream is X,Y =' 308 print s 309 fid.write(s) 310 s= 'Point 1m Down Stream is X,Y =' 311 print s 312 fid.write(s) 313 314 328 s = 'Culvert Direction is %s\n' %str(self.vector) 329 log(fid, s) 315 330 316 331 def __call__(self, domain): … … 319 334 from anuga.config import g, epsilon 320 335 from Numeric import take, sqrt 321 322 # OLE I think this routine gets Called twice every time Step, so that each time the second time around Delta_t =0 and therefore Momentum is Set to ZERO 323 # Why does this NEED to be calle d twice ??? 324 325 326 # Get average water depths at each opening 327 336 337 338 # Time stuff 328 339 time = domain.get_time() 329 openings = self.openings # There are two Opening [0] and [1] so this happens TWICE Right ??????????? 340 delta_t = time-self.last_time 341 s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) 342 log(fid, s) 343 344 msg = 'Time did not advance' 345 if time > 0.0: assert delta_t > 0.0, msg 346 347 348 # Get average water depths at each opening 349 openings = self.openings # There are two Opening [0] and [1] 330 350 for i, opening in enumerate(openings): 331 351 stage = domain.quantities['stage'].get_values(location='centroids', 332 indices=opening. indices)352 indices=opening.exchange_indices) 333 353 elevation = domain.quantities['elevation'].get_values(location='centroids', 334 indices=opening. indices)354 indices=opening.exchange_indices) 335 355 336 356 # Indices corresponding to energy enquiry field for this opening 337 357 coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y) 338 idx= inside_polygon(coordinates, self.enquiry_polygons[i])358 enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 339 359 340 360 if self.verbose: … … 344 364 345 365 346 # Get average model values for points in ?????? WHAT IF THE DOMAIN IS DRY ??????? 347 # enquiry polygon for this opening 366 # Get model values for points in enquiry polygon for this opening 348 367 dq = domain.quantities 349 stage = dq['stage'].get_values(location='centroids', indices=idx) 350 xmomentum = dq['xmomentum'].get_values(location='centroids', indices=idx) 351 ymomentum = dq['ymomentum'].get_values(location='centroids', indices=idx) 352 elevation = dq['elevation'].get_values(location='centroids', indices=idx) 353 depth = stage - elevation 354 355 # Compute mean velocity in the area in front of the culvert (taking zero depths into account) 356 ux = xmomentum/(depth+epsilon) # Does epsilon.... handle a Dry Cell Case ???????? 357 uy = ymomentum/(depth+epsilon) 358 v = mean(sqrt(ux**2+uy**2)) 359 d = mean(depth) 360 w = mean(stage) 361 z = mean(elevation) 362 363 # Ratio of depth to culvert height 364 #ratio = d/(2*self.height) # Is Height Now REAL Height or Radius if Radius then Need 2 else NOT !! 365 ratio = d/self.height # Use this if Height is NOT Radius! 366 if ratio > 1.0: # Assume culvert is running full & 367 ratio = 1.0 # under pressure. Note this is usually ~ 1.35 368 stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices) 369 xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices) 370 ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices) 371 elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices) 372 depth = stage - elevation 373 374 # Compute mean values of selected quantitites in the enquiry area in front of the culvert 375 # Epsilon handles a dry cell case 376 ux = xmomentum/(depth+velocity_protection/depth) # Velocity (x-direction) 377 uy = ymomentum/(depth+velocity_protection/depth) # Velocity (y-direction) 378 v = mean(sqrt(ux**2+uy**2)) # Average velocity 379 w = mean(stage) # Average stage 380 381 # Store values at enquiry field 382 opening.velocity = v 383 384 385 # Compute mean values of selected quantitites in the exchange area in front of the culvert 386 # Stage and velocity comes from enquiry area and elevation from exchange area 387 388 # Use invert level instead of elevation if specified 389 invert_level = self.invert_levels[i] 390 if invert_level is not None: 391 z = invert_level 392 else: 393 elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices) 394 z = mean(elevation) # Average elevation 395 396 # Estimated depth above the culvert inlet 397 d = w - z 398 399 if d < 0.0: 400 # This is possible since w and z are taken at different locations 401 #msg = 'D < 0.0: %f' %d 402 #raise msg 403 d = 0.0 404 405 # Ratio of depth to culvert height. 406 # If ratio > 1 then culvert is running full 407 ratio = d/self.height 368 408 opening.ratio = ratio 369 409 370 371 # Average measure of total energy (D + V^2/2g) in enquiry field in front of this opening 372 Es = d + 0.5*v**2/g # SPECIFIC ENERGY 373 Et = w + 0.5*v**2/g # TOTAL ENERGY 374 opening.energy = Et 375 opening.vel=v 376 377 378 # Store current average stage and depth at enquiry field with each opening object 410 # Average measures of energy in front of this opening 411 Es = d + 0.5*v**2/g # Specific energy in exchange area 412 Et = w + 0.5*v**2/g # Total energy in the enquiry area 413 opening.total_energy = Et 414 opening.specific_energy = Es 415 416 # Store current average stage and depth with each opening object 379 417 opening.depth = d 380 418 opening.stage = w 381 419 opening.elevation = z 382 383 # End of the FOR LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 384 # we now need to deal with each opening 385 386 # Now we are done calculating energy etc for each opening. 387 # At this point we can work on the transfer functions etc. 388 389 # Handy values (all calculated at enquiry polygon - 390 # if you need them at exchange polygons we can easily do that.) 391 392 393 #Then we need delta_h or better delta_E to be used to drive water through culvert IS THIS FROM EXCHANGE or ENQUIRY POLYGON ?????? Need to Know !!!! 394 #delta_h = openings[0].energy - openings[1].energy 395 delta_h = openings[0].stage - openings[1].stage 396 delta_Et = openings[0].energy - openings[1].energy 397 delta_z = openings[1].elevation - openings[0].elevation 398 culvert_slope=(delta_z/self.distance) # COULD USE CULVERT SLOPE To Determine Internal Flow Velocity That is Not Inlet Or Outlet Velocity 399 s= 'Delta_Energy = %.3f'%(delta_Et) 400 print s 401 fid.write(s) 402 403 # Ideas..... GET RID OF THIS !!!!!!!!!!!!!!!! 404 if openings[0].depth > 0 and openings[1].depth > 0: 405 flow_rate = delta_h 406 else: 407 flow_rate = 0.0 408 409 if openings[0].depth > self.height: 410 # This could be usefull mayby 411 #print 'Inlet has been overflowed' 412 pass 413 414 if openings[1].depth > self.height: 415 # This could be usefull mayby 416 #print 'Outlet has been overflowed' 417 pass 418 # ====================== NOW ENTER INTO THE CULVERT EQUATIONS AS DEFINED BY BOYD GENERALISED METHOD 419 # == The quantity of flow passing through a culvert is controlled by many factors 420 # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation 421 # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation 422 # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled 423 420 421 422 ################# End of the FOR loop ################################################ 423 424 425 # We now need to deal with each opening individually 426 427 # Determine flow direction based on total energy difference 428 delta_Et = openings[0].total_energy - openings[1].total_energy 429 424 430 if delta_Et > 0: 425 431 #print 'Flow U/S ---> D/S' 426 432 inlet=openings[0] 427 433 outlet=openings[1] 434 435 inlet.momentum = self.opening_momentum[0] 436 outlet.momentum = self.opening_momentum[1] 428 437 else: 429 438 #print 'Flow D/S ---> U/S' … … 431 440 outlet=openings[0] 432 441 433 if inlet.depth < 0.01: #epsilon: 10mm Minium should be OK....for NOW 434 # Do nothing because No water over Opening That is Set Flow to Zero!! 435 #print 'NO water over Inlet yet!' 436 self.openings[0].rate = 0 437 self.openings[1].rate = 0 438 Qb1=0.0 439 Qb2=0.0 440 Qb3tw=0.0 441 Qb3dc=0.0 442 443 444 else: #For PIPES 445 print 'Water over Inlet NOW!' 442 inlet.momentum = self.opening_momentum[1] 443 outlet.momentum = self.opening_momentum[0] 444 445 delta_Et = -delta_Et 446 447 msg = 'Total energy difference is negative' 448 assert delta_Et > 0.0, msg 449 450 delta_h = inlet.stage - outlet.stage 451 delta_z = inlet.elevation - outlet.elevation 452 culvert_slope = (delta_z/self.distance) 453 454 if culvert_slope < 0.0: 455 # Adverse gradient - flow is running uphill 456 # Flow will be purely controlled by uphill outlet face 457 print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation 458 459 460 s = 'Delta total energy = %.3f' %(delta_Et) 461 log(fid, s) 462 463 464 # ====================== NOW ENTER INTO THE CULVERT EQUATIONS AS DEFINED BY BOYD GENERALISED METHOD 465 # == The quantity of flow passing through a culvert is controlled by many factors 466 # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation 467 # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation 468 # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled 469 470 471 Q_outlet_tailwater = 0.0 472 inlet.rate = 0.0 473 outlet.rate = 0.0 474 Q_inlet_unsubmerged = 0.0 475 Q_inlet_submerged = 0.0 476 Q_outlet_critical_depth = 0.0 477 478 if inlet.depth >= 0.01: 479 # Water has risen above inlet 446 480 if self.width==self.height: # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ??? 447 481 pass … … 450 484 #PipeDcrit= 451 485 #Delta_E=HW-TW 452 else: #FOR BOX CULVERT 453 # Why does it not KNOW these Losses from __INIT__ ???? 454 #sum_loss=loss_exit+loss_entry+loss_bend+loss_special # Note if losses are a function can't do this !!!! 455 sum_loss=1.5 456 manning=0.012 # THis should come from Culvert data !!! 457 print 'Depth over inlet =',inlet.depth 458 Qb1=0.540*g**0.5*self.width*inlet.depth**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged 459 Qb2=0.702*g**0.5*self.width*self.height**0.89*inlet.depth**0.61 # Flow based on Inlet Ctrl Inlet Submerged 460 print'Qb1= %.6f, Qb2 = %.6f'%(Qb1,Qb2) 461 Qb_Inlet=min(Qb1,Qb2) 462 # NOW FOR OUTLET CONTROL 463 # Need to determine the depth at the outlet relative to the depth of flow in the Culvert 464 # THIS IS WHERE THE TROUBLE STARTS 465 # USUALLY YOU NEED TO ITERATE THECALCULATION OF Q as itis Reliant on Knowing Dcrit which is reliant on Q 466 # Instead we will calculate Q assuming TW and then Q again using Dcrit assuming the Qtw 467 #If TW >H then Assume Dcrit =H else Assume Dcrit= TW 468 if outlet.depth > self.height: # The Outlet is Submerged 469 outlet_culv_depth=self.height 470 flow_area=self.width*self.height 471 perimeter=2.0*(self.width+self.height) 472 hyd_rad=flow_area/perimeter 473 elif outlet.depth==0.0: 474 outlet_culv_depth=inlet.depth # For purpose of calculation assume the outlet depth = the inlet depth 475 flow_area=self.width*inlet.depth 476 perimeter=(self.width+2.0*inlet.depth) 477 hyd_rad=flow_area/perimeter 478 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 479 outlet_culv_depth=outlet.depth 480 flow_area=self.width*outlet.depth 481 perimeter=(self.width+2.0*outlet.depth) 482 hyd_rad=flow_area/perimeter 483 print 'hydrad at outlet= ',hyd_rad 484 Vel3=sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333)) # Outlet Ctrl usingTW as Control Depth Either TW or H 485 Qb3tw=flow_area*Vel3 486 Qb3dc=Qhold=Qb_Inlet 487 Qcheck=10.0 488 while Qcheck > 0.001: 489 Dcrit = 0.4672*(Qb3dc/self.width)**0.66666 # This assumes that Flow is at Critical Depth Not Alsways True... 490 # Better to Use Mannings Equation to Estimate Q and iterate Depth in the Culvert to Match the Energy Loss 491 print 'Dcrit =',Dcrit 492 if Dcrit > self.height: # 493 flow_area=self.width*self.height 486 else: 487 # Box culvert 488 489 sum_loss=self.loss_exit+self.loss_entry+self.loss_bend+self.loss_special 490 manning=self.manning 491 492 # Calculate flows for inlet control 493 E = inlet.specific_energy 494 #E = min(inlet.specific_energy, delta_Et) 495 496 s = 'Driving energy = %f m' %E 497 log(fid, s) 498 499 msg = 'Driving energy is negative' 500 assert E >= 0.0, msg 501 502 Q_inlet_unsubmerged = 0.540*g**0.5*self.width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged 503 Q_inlet_submerged = 0.702*g**0.5*self.width*self.height**0.89*E**0.61 # Flow based on Inlet Ctrl Inlet Submerged 504 505 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged) 506 log(fid, s) 507 508 case = '' 509 if Q_inlet_unsubmerged < Q_inlet_submerged: 510 Q = Q_inlet_unsubmerged 511 flow_area = self.width*inlet.depth 512 outlet_culv_depth = inlet.depth 513 case = 'Inlet unsubmerged' 514 else: 515 Q = Q_inlet_submerged 516 flow_area = self.width*self.height 517 outlet_culv_depth = self.height 518 case = 'Inlet submerged' 519 520 if delta_Et < Es: 521 # Calculate flows for outlet control 522 # Determine the depth at the outlet relative to the depth of flow in the Culvert 523 524 sum_loss = self.sum_loss 525 526 if outlet.depth > self.height: # The Outlet is Submerged 527 outlet_culv_depth=self.height 528 flow_area=self.width*self.height # Cross sectional area of flow in the culvert 494 529 perimeter=2.0*(self.width+self.height) 495 hyd_rad=flow_area/perimeter 496 else: 497 flow_area=self.width*Dcrit 498 perimeter=(self.width+2.0*Dcrit) 499 hyd_rad=flow_area/perimeter 500 print 'Adjusted hydrad at outlet= ',hyd_rad 501 502 Qb3dc=flow_area*sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333)) # Outlet Ctrl using Dcrit as Control Depth Either Dcrit or H 503 Qcheck =abs(Qhold-Qb3dc) 504 Qhold=Qb3dc 505 print 'Qcheck Qb3dc ',Qcheck,Qb3dc 506 flow_rate_control=min(Qb1,Qb2) # 507 print 'OutletCtrl Qbtw and Qbdc= ',Qb3tw,Qb3dc 508 print 'Flow Rate Control = ',flow_rate_control 509 510 # NOW UPDATE Flows on the DOMAIN once again have to check direction of FLOW 511 if delta_Et > 0: 512 print 'Flow U/S ---> D/S' 513 self.openings[0].rate=-flow_rate_control 514 self.openings[1].rate=flow_rate_control 515 else: 516 print 'Flow D/S ---> U/S' 517 self.openings[1].rate=-flow_rate_control 518 self.openings[0].rate=flow_rate_control # Flow Rates have been set NOW Add jet in the form of absolute momentum to outlet 519 520 530 case = 'Outlet submerged' 531 elif outlet.depth==0.0: 532 outlet_culv_depth=inlet.depth # For purpose of calculation assume the outlet depth = the inlet depth 533 flow_area=self.width*inlet.depth 534 perimeter=(self.width+2.0*inlet.depth) 535 case = 'Outlet depth is zero' 536 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 537 outlet_culv_depth=outlet.depth 538 flow_area=self.width*outlet.depth 539 perimeter=(self.width+2.0*outlet.depth) 540 case = 'Outlet is open channel flow' 541 542 hyd_rad = flow_area/perimeter 543 s = 'hydraulic radius at outlet = %f' %hyd_rad 544 log(fid, s) 545 546 547 # Outlet control velocity using tail water 548 culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333)) 549 Q_outlet_tailwater = flow_area * culvert_velocity 550 551 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater 552 log(fid, s) 553 554 Q = min(Q, Q_outlet_tailwater) 555 556 557 558 log(fid, 'Case: "%s"' %case) 559 560 flow_rate_control=Q 561 562 s = 'Flow Rate Control = %f' %flow_rate_control 563 log(fid, s) 564 565 inlet.rate = -flow_rate_control 566 outlet.rate = flow_rate_control 567 521 568 culv_froude=sqrt(flow_rate_control**2*self.width/(g*flow_area**3)) 522 print 'Froude in Culvert = ',culv_froude 523 524 #+++++++++++ NOW DETERMINE MOMENTUM AT THE OUTLET !!!! +++++++++++++++ 525 526 527 # This Should be Based on the VELOCITY in the Culvert 528 # Based on the Flow Depth in the Culvert for part full So need to determine actual flow depth in the culvert based on slope !!! 529 # or Flow Jet of the Pressurised Culvert if FUll 530 # If Part Full Flow Calculate Part Full Depth 531 # Based on Depth Calculate Area... the Vel = flow_rate_control / Area 532 # Need to Break Velocity into X & Y Components 533 #self.openings[1].set_quantity_values(delta_h*speed, 'xmomentum') 534 # Previous Calculated Depth/ Culvert Heigh Ratio Use it to Determine Velocity !!!! FIX LAter 535 536 #if (ratio*self.area)== 0: # Don't Really need this already established water depth here 537 # outlet_vel=0.0 538 #else: 539 540 541 # FIXME (Ole): Shouldn't this be openings[1]?? 542 # This attempts to resolve the Velocity in the Culvert at the Entrance and Outlet.... Again if flow reversal Occurs then need to swap END !!! 543 #inlet_vel =(flow_rate_control/(openings[0].ratio*openings[0].area)) 544 #inlet_dep = openings[0].depth*openings[0].ratio #????? 545 #inlet_mom = outlet_vel*outlet_dep 546 #inlet_E_Loss= 0.5*0.5*inlet_vel**2/g 547 # NEED to Calculate Critical depth inthe Culvert for Box & Pipe Sections.. If Greater than Depth on the Down Stream Flood Plain use Dcrit 548 outlet_vel =(flow_rate_control/flow_area) 549 outlet_dep = flow_area/self.width #????? 550 outlet_mom = outlet_vel*outlet_dep 551 outlet_E_Loss= 0.5*0.5*outlet_vel**2/g # Ke v^2/2g 552 print 'Outlet Velocity =',outlet_vel 553 554 #barrel_vel=(inlet_vel+outlet_vel)/2 555 # use this to determine Barrel Loss 556 557 # Eventually will Need Momentum in X & Y Components based on the orientation of 558 # the culvert from X0,Y0, X1, Y1 from Create Polygon Routine 559 # YES - use self.vector which is a unit vector in the direction of the culvert. 560 # Will also need to Check Normal Depth For a Steep Culvert using the Mannings Equation to determine the Flow Depth in the Culvert this requires Iteration and a Simpsons Rule Approach to Refine Depth 561 # So that Qcalculated = Qestimated.... 562 563 outlet_mom_x, outlet_mom_y = self.vector * outlet_mom 564 #print 'Directional momentum', outlet_mom_x, outlet_mom_y 565 xmomentum_rate = outlet_mom_x 566 ymomentum_rate = outlet_mom_y 567 print 'X Y MOM ',xmomentum_rate,ymomentum_rate 568 # Update the momentum forcing terms with directional momentum at the outlet 569 # >>>>>>>>>>>>>>OLE Passing through this the 2nd time ALWAYS sets delta_t to 0 so that Momentum is also set to 0.0 This is NOT RIGHT !!!! THIS SEGMENT OF CODE NEEDS ONLY TO BE RUn FOR THE OUTLET !!! 570 # ONLY NEEDS TO BE CALLED ONCE for the CULVERT not TWICE for Each Opening as a check of each opening is DONE HERe any WAY !!! 571 572 # I have tried to trick it by introducing this term temp_keep_Delta_t...... but this is really bad !!! How do we avoid this Situation ???? 569 s = 'Froude in Culvert = %f' %culv_froude 570 log(fid, s) 571 572 573 574 # Determine momentum at the outlet 575 barrel_velocity = Q/(flow_area + velocity_protection/flow_area) 576 barrel_momentum = barrel_velocity*outlet_culv_depth 577 578 outlet_E_Loss= 0.5*0.5*barrel_velocity**2/g # Ke v^2/2g 579 s = 'Barrel velocity = %f' %barrel_velocity 580 log(fid, s) 581 582 # Compute momentum vector at outlet 583 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum 584 585 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) 586 log(fid, s) 587 573 588 delta_t = time - self.last_time 574 if delta_t==0:575 delta_t=self.temp_keep_delta_t576 xmomentum_rate=self.xmom_forcing1.rate577 ymomentum_rate=self.ymom_forcing1.rate578 self.temp_keep_delta_t=delta_t579 print 'time and delta t = ',time,delta_t580 589 if delta_t > 0.0: 581 print 'delta_t Calc MOM' 582 xmomentum_rate = outlet_mom_x - self.xmom_forcing1.value 590 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value 583 591 xmomentum_rate /= delta_t 584 592 585 ymomentum_rate = outlet_mom_y - self.ymom_forcing1.value593 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value 586 594 ymomentum_rate /= delta_t 587 print 'X Y MOM_RATE ',xmomentum_rate,ymomentum_rate 588 #ymomentum_rate = 3 589 #else: 590 # xmomentum_rate = ymomentum_rate = 0.0 595 596 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate) 597 log(fid, s) 598 else: 599 xmomentum_rate = ymomentum_rate = 0.0 600 591 601 592 602 # Set momentum rates for outlet jet 593 self.xmom_forcing1.rate = xmomentum_rate 594 self.ymom_forcing1.rate = ymomentum_rate 595 596 # Remember this value for next step (IMPORTANT) 597 self.xmom_forcing1.value = outlet_mom_x 598 self.ymom_forcing1.value = outlet_mom_y 599 603 outlet.momentum[0].rate = xmomentum_rate 604 outlet.momentum[1].rate = ymomentum_rate 605 606 # Remember this value for next step (IMPORTANT) 607 outlet.momentum[0].value = outlet_mom_x 608 outlet.momentum[1].value = outlet_mom_y 600 609 601 610 if int(domain.time*100) % 100 == 0: 602 611 s = 'T=%.5f, Culvert Discharge = %.3f Culv. Vel. %0.3f'\ 603 %(time, flow_rate_control, outlet_vel)612 %(time, flow_rate_control, barrel_velocity) 604 613 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\ 605 %(outlet_ dep, outlet_mom_x,outlet_mom_y)614 %(outlet_culv_depth, outlet_mom_x,outlet_mom_y) 606 615 s += ' Momentum rate: (%.4f, %.4f)'\ 607 616 %(xmomentum_rate, ymomentum_rate) 608 s+='Outlet Vel= %.3f Outlet Dep = %.3f'\ 609 %(outlet_vel,outlet_dep) 610 fid.write(s) 611 print s 617 s+='Outlet Vel= %.3f'\ 618 %(barrel_velocity) 619 log(fid, s) 612 620 613 621 … … 619 627 # Execute momentum terms 620 628 # This is where Inflow objects are evaluated and update the domain 621 self.xmom_forcing0(domain) 622 self.ymom_forcing0(domain) 623 self.xmom_forcing1(domain) 624 self.ymom_forcing1(domain) 625 626 627 628 # Print out flows every 1 seconds 629 if int(time*100) % 100 == 0: 630 s = 'Time %.6f\n' %time 631 s += ' Opening[0]: d=%.3f, v=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\ 632 %(openings[0].depth, 633 openings[0].vel, 634 openings[0].area, 635 openings[0].energy, 636 openings[0].ratio) 637 s += ' Opening[1]: d=%.3f, v=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\ 638 %(openings[1].depth, 639 openings[1].vel, 640 openings[1].area, 641 openings[1].energy, 642 openings[1].ratio) 643 s += ' Distance=%.2f, Qb1=%.3f, Qb2=%.3f, Qb3tw=%.3f,Qb3dc=%.3f\n'\ 644 %(self.distance, 645 Qb1, 646 Qb2, 647 Qb3tw,Qb3dc) 648 649 print s 650 fid.write(s) 651 629 outlet.momentum[0](domain) 630 outlet.momentum[1](domain) 631 652 632 # Store value of time 653 633 self.last_time = time 654 #self.temp_keep_delta_t=delta_t 634 635 636 637 655 638 #------------------------------------------------------------------------------ 656 639 # Setup culvert inlets and outlets in current topography
Note: See TracChangeset
for help on using the changeset viewer.
