Changeset 5429
- Timestamp:
- Jun 25, 2008, 4:15:10 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/culvert_flow/culvert_boyd_channel.py
r5428 r5429 53 53 from anuga.culvert_flows.culvert_polygons import create_culvert_polygons 54 54 from anuga.utilities.polygon import plot_polygons 55 from anuga.utilities.system_tools import log_to_file 56 57 55 58 56 59 from math import pi,sqrt … … 60 63 #------------------------------------------------------------------------------ 61 64 62 63 def log(fid, s):64 print s65 fid.write(s + '\n')66 67 65 68 66 # Open file for storing some specific results... … … 89 87 90 88 s='Size'+str(len(domain)) 91 log(fid, s) 92 93 velocity_protection = 1.0e-4 89 log_to_file(fid, s) 90 91 92 93 velocity_protection = 1.0e-5 94 95 96 97 98 99 100 def boyd_generalised_culvert_model(inlet, outlet, delta_Et, width, height, length, sum_loss, manning, g, fid): 101 """Boyd's generalisation of the US department of transportation culvert model 102 # == The quantity of flow passing through a culvert is controlled by many factors 103 # == 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 104 # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation 105 # == 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 106 """ 107 108 from anuga.utilities.system_tools import log_to_file 109 110 111 Q_outlet_tailwater = 0.0 112 inlet.rate = 0.0 113 outlet.rate = 0.0 114 Q_inlet_unsubmerged = 0.0 115 Q_inlet_submerged = 0.0 116 Q_outlet_critical_depth = 0.0 117 118 if inlet.depth >= 0.01: 119 # Water has risen above inlet 120 if width==height: # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ??? 121 pass 122 #Q1[t]= 0.421*g**0.5*Diam**0.87*HW**1.63 # Inlet Ctrl Inlet Unsubmerged 123 #Q2[t]= 0.530*g**0.5*Diam**1.87*HW**0.63 # Inlet Ctrl Inlet Submerged 124 #PipeDcrit= 125 #Delta_E=HW-TW 126 else: 127 # Box culvert 128 129 # Calculate flows for inlet control 130 E = inlet.specific_energy 131 132 s = 'Driving energy = %f m' %E 133 log_to_file(fid, s) 134 135 msg = 'Driving energy is negative' 136 assert E >= 0.0, msg 137 138 Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged 139 Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61 # Flow based on Inlet Ctrl Inlet Submerged 140 141 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged) 142 log_to_file(fid, s) 143 144 case = '' 145 if Q_inlet_unsubmerged < Q_inlet_submerged: 146 Q = Q_inlet_unsubmerged 147 flow_area = width*inlet.depth 148 outlet_culvert_depth = inlet.depth 149 case = 'Inlet unsubmerged' 150 else: 151 Q = Q_inlet_submerged 152 flow_area = width*height 153 outlet_culvert_depth = height 154 case = 'Inlet submerged' 155 156 if delta_Et < E: 157 # Calculate flows for outlet control 158 # Determine the depth at the outlet relative to the depth of flow in the Culvert 159 160 if outlet.depth > height: # The Outlet is Submerged 161 outlet_culvert_depth=height 162 flow_area=width*height # Cross sectional area of flow in the culvert 163 perimeter=2.0*(width+height) 164 case = 'Outlet submerged' 165 elif outlet.depth==0.0: 166 outlet_culvert_depth=inlet.depth # For purpose of calculation assume the outlet depth = the inlet depth 167 flow_area=width*inlet.depth 168 perimeter=(width+2.0*inlet.depth) 169 case = 'Outlet depth is zero' 170 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 171 outlet_culvert_depth=outlet.depth 172 flow_area=width*outlet.depth 173 perimeter=(width+2.0*outlet.depth) 174 case = 'Outlet is open channel flow' 175 176 hyd_rad = flow_area/perimeter 177 s = 'hydraulic radius at outlet = %f' %hyd_rad 178 log_to_file(fid, s) 179 180 181 # Outlet control velocity using tail water 182 culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 183 Q_outlet_tailwater = flow_area * culvert_velocity 184 185 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater 186 log_to_file(fid, s) 187 188 Q = min(Q, Q_outlet_tailwater) 189 190 191 192 log_to_file(fid, 'Case: "%s"' %case) 193 194 flow_rate_control=Q 195 196 s = 'Flow Rate Control = %f' %flow_rate_control 197 log_to_file(fid, s) 198 199 inlet.rate = -flow_rate_control 200 outlet.rate = flow_rate_control 201 202 culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3)) 203 s = 'Froude in Culvert = %f' %culv_froude 204 log_to_file(fid, s) 205 206 # Determine momentum at the outlet 207 barrel_velocity = Q/(flow_area + velocity_protection/flow_area) 208 209 return Q, barrel_velocity, outlet_culvert_depth 210 94 211 95 212 … … 137 254 z[i] += embank_hgt-embank_dnslope*(x[i] -12.0) # Sloping D/S Face 138 255 139 140 # Constriction141 #if 27 < x[i] < 29 and y[i] > 3:142 # z[i] += 2143 144 # Pole145 #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:146 # z[i] += 2147 148 # HOLE For Pit at Opening[0]149 #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2:150 # z[i]-=1151 152 # HOLE For Pit at Opening[1]153 #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2:154 # z[i]-=1155 156 256 return z 157 257 … … 237 337 blockage_topdwn=None, 238 338 blockage_bottup=None, 339 culvert_routine=None, 239 340 verbose=False): 240 341 … … 252 353 if blockage_topdwn is None: blockage_topdwn=0.00 253 354 if blockage_bottup is None: blockage_bottup=0.00 355 if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model 254 356 255 357 … … 303 405 self.blockage_bottup = blockage_bottup 304 406 407 # Store culvert routine 408 self.culvert_routine = culvert_routine 409 305 410 306 411 # Create objects to update momentum (a bit crude at this stage) … … 324 429 # Print something 325 430 s = 'Culvert Effective Length = %.2f m' %(self.distance) 326 log (fid, s)431 log_to_file(fid, s) 327 432 328 433 s = 'Culvert Direction is %s\n' %str(self.vector) 329 log (fid, s)434 log_to_file(fid, s) 330 435 331 436 def __call__(self, domain): … … 340 445 delta_t = time-self.last_time 341 446 s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) 342 log (fid, s)447 log_to_file(fid, s) 343 448 344 449 msg = 'Time did not advance' … … 358 463 enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 359 464 360 if self.verbose:361 pass362 #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\363 # %(i, len(idx), self.enquiry_polygons[i])364 365 465 366 466 # Get model values for points in enquiry polygon for this opening … … 434 534 435 535 inlet.momentum = self.opening_momentum[0] 436 outlet.momentum = self.opening_momentum[1] 536 outlet.momentum = self.opening_momentum[1] 537 437 538 else: 438 539 #print 'Flow D/S ---> U/S' … … 459 560 460 561 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 480 if self.width==self.height: # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ??? 481 pass 482 #Q1[t]= 0.421*g**0.5*Diam**0.87*HW**1.63 # Inlet Ctrl Inlet Unsubmerged 483 #Q2[t]= 0.530*g**0.5*Diam**1.87*HW**0.63 # Inlet Ctrl Inlet Submerged 484 #PipeDcrit= 485 #Delta_E=HW-TW 486 else: 487 # Box culvert 562 log_to_file(fid, s) 563 564 sum_loss = self.sum_loss 565 manning=self.manning 566 length = self.distance 567 height = self.height 568 width = self.width 569 570 Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(inlet, outlet, delta_Et, width, height, length, sum_loss, manning, g, fid) 571 572 ##################################################### 573 barrel_momentum = barrel_velocity*culvert_outlet_depth 488 574 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 575 s = 'Barrel velocity = %f' %barrel_velocity 576 log_to_file(fid, s) 577 578 # Compute momentum vector at outlet 579 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum 580 581 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) 582 log_to_file(fid, s) 583 584 delta_t = time - self.last_time 585 if delta_t > 0.0: 586 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value 587 xmomentum_rate /= delta_t 525 588 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 529 perimeter=2.0*(self.width+self.height) 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) 589 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value 590 ymomentum_rate /= delta_t 545 591 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) 592 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate) 593 log_to_file(fid, s) 594 else: 595 xmomentum_rate = ymomentum_rate = 0.0 596 597 598 # Set momentum rates for outlet jet 599 outlet.momentum[0].rate = xmomentum_rate 600 outlet.momentum[1].rate = ymomentum_rate 601 602 # Remember this value for next step (IMPORTANT) 603 outlet.momentum[0].value = outlet_mom_x 604 outlet.momentum[1].value = outlet_mom_y 605 606 if int(domain.time*100) % 100 == 0: 607 s = 'T=%.5f, Culvert Discharge = %.3f f'\ 608 %(time, Q) 609 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\ 610 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y) 611 s += ' Momentum rate: (%.4f, %.4f)'\ 612 %(xmomentum_rate, ymomentum_rate) 613 s+='Outlet Vel= %.3f'\ 614 %(barrel_velocity) 615 log_to_file(fid, s) 616 559 617 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 568 culv_froude=sqrt(flow_rate_control**2*self.width/(g*flow_area**3)) 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 588 delta_t = time - self.last_time 589 if delta_t > 0.0: 590 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value 591 xmomentum_rate /= delta_t 592 593 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value 594 ymomentum_rate /= delta_t 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 601 602 # Set momentum rates for outlet jet 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 609 610 if int(domain.time*100) % 100 == 0: 611 s = 'T=%.5f, Culvert Discharge = %.3f Culv. Vel. %0.3f'\ 612 %(time, flow_rate_control, barrel_velocity) 613 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\ 614 %(outlet_culv_depth, outlet_mom_x,outlet_mom_y) 615 s += ' Momentum rate: (%.4f, %.4f)'\ 616 %(xmomentum_rate, ymomentum_rate) 617 s+='Outlet Vel= %.3f'\ 618 %(barrel_velocity) 619 log(fid, s) 620 618 619 621 620 622 621 # Execute flow term for each opening … … 650 649 end_point0=[9.0, 2.5], 651 650 end_point1=[13.0, 2.5], 652 width=1.20,height=0.75, 651 width=1.20,height=0.75, 652 culvert_routine=boyd_generalised_culvert_model, 653 653 verbose=True) 654 654
Note: See TracChangeset
for help on using the changeset viewer.
