Changeset 7495
- Timestamp:
- Sep 7, 2009, 6:47:36 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/culvert_flows/culvert_class.py
r7323 r7495 136 136 length=None, 137 137 number_of_barrels=1, 138 number_of_smoothing_steps=2000, 138 139 trigger_depth=0.01, # Depth below which no flow happens 139 140 manning=None, # Mannings Roughness for Culvert … … 152 153 # Input check 153 154 154 if height is None: height = width 155 self.height = height 156 self.width = width 157 155 if height is None: height = width 158 156 159 157 assert number_of_barrels >= 1 160 158 assert use_velocity_head is True or use_velocity_head is False 161 159 162 msg = 'Momentum jet not yet moved to general culvert' 163 assert use_momentum_jet is False, msg 164 160 #msg = 'Momentum jet not yet moved to general culvert' 161 #assert use_momentum_jet is False, msg 162 self.use_momentum_jet = use_momentum_jet 163 165 164 self.culvert_routine = culvert_routine 166 165 self.culvert_description_filename = culvert_description_filename … … 168 167 label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename) 169 168 self.rating_curve = ensure_numeric(rating_curve) 169 170 self.height = height 171 self.width = width 172 170 173 171 174 self.domain = domain … … 292 295 self.verbose = verbose 293 296 294 297 # Circular index for flow averaging in culvert 298 self.N = N = number_of_smoothing_steps 299 self.Q_list = [0]*N 300 self.i = i 295 301 296 302 # For use with update_interval … … 298 304 self.update_interval = update_interval 299 305 306 # Create objects to update momentum (a bit crude at this stage). This is used with the momentum jet. 307 xmom0 = General_forcing(domain, 'xmomentum', 308 polygon=P['exchange_polygon0']) 309 310 xmom1 = General_forcing(domain, 'xmomentum', 311 polygon=P['exchange_polygon1']) 312 313 ymom0 = General_forcing(domain, 'ymomentum', 314 polygon=P['exchange_polygon0']) 315 316 ymom1 = General_forcing(domain, 'ymomentum', 317 polygon=P['exchange_polygon1']) 318 319 self.opening_momentum = [[xmom0, ymom0], [xmom1, ymom1]] 320 321 300 322 301 323 # Print some diagnostics to log if requested … … 402 424 def adjust_flow_for_available_water_at_inlet(self, Q, delta_t): 403 425 """Adjust Q downwards depending on available water at inlet 426 427 This is a critical step in modelling bridges and Culverts 428 the predicted flow through a structure based on an abstract 429 algorithm can at times request for water that is simply not 430 available due to any number of constrictions that limit the 431 flow approaching the structure In order to ensure that 432 there is adequate flow available certain checks are 433 required There needs to be a check using the Static Water 434 Volume sitting infront of the structure, In addition if the 435 water is moving the available water will be larger than the 436 static volume 437 438 NOTE To temporarily switch this off for Debugging purposes 439 rem out line in function def compute_rates below 404 440 """ 405 441 … … 423 459 indices=[idx]) 424 460 depth = stage-elevation 425 min_depth = min(depth.flat) 426 427 # Compute possible flow for exchange region based on 428 # triangle with smallest depth 429 max_Q = min_depth*I.exchange_area/delta_t 430 461 min_depth = min(depth.flat) # This may lead to errors if edge of area is at a higher level !!!! 462 avg_depth = mean(depth.flat) # Yes, but this one violates the conservation unit tests 463 464 465 466 # FIXME (Ole): If you want these, use log.critical() and 467 # make the statements depend on verbose 468 #print I.depth 469 #print I.velocity 470 #print self.width 471 472 # max_Q Based on Volume Calcs 473 474 475 depth_term = min_depth*I.exchange_area/delta_t 476 if min_depth < 0.2: 477 # Only add velocity term in shallow waters (< 20 cm) 478 # This is a little ad hoc, but maybe it is reasonable 479 velocity_term = self.width*min_depth*I.velocity 480 else: 481 velocity_term = 0.0 482 483 # This one takes approaching water into account 484 max_Q = max(velocity_term, depth_term) 485 486 # This one preserves Volume 487 #max_Q = depth_term 488 489 490 if self.verbose is True: 491 log.critical('Max_Q = %f' % max_Q) 492 msg = 'Width = %.2fm, Depth at inlet = %.2f m, Velocity = %.2f m/s. ' % (self.width, I.depth, I.velocity) 493 msg += 'Max Q = %.2f m^3/s' %(max_Q) 494 log.critical(msg) 495 496 if self.log_filename is not None: 497 log_to_file(self.log_filename, msg) 498 # New Procedure for assessing the flow available to the Culvert 499 # This routine uses the GET FLOW THROUGH CROSS SECTION 500 # Need to check Several Polyline however 501 # Firstly 3 sides of the exchange Poly 502 # then only the Line Directly infront of the Polygon 503 # Access polygon Points from self.inlet.polygon 504 505 # The Following computes the flow crossing over 3 sides of the exchange polygon for the structure 506 # Clearly the flow in the culvert can not be more than that flowing toward it through the exhange polygon 507 508 #q1 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][1:3]) # First Side Segment 509 #q2 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][2:]) # Second Face Segment 510 #q3 =domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'].take([3,0], axis=0)) # Third Side Segment 511 # q4 = domain.get_flow_through_cross_section([self.culvert_polygons['exchange_polygon0'][1:4]][0]) 512 #q4=max(q1,0.0)+max(q2,0.0)+max(q3,0.0) 513 # To use only the Flow crossing the 3 sides of the Exchange Polygon use the following Line Only 514 #max_Q=max(q1,q2,q3,q4) 515 # Try Simple Smoothing using Average of 2 approaches 516 #max_Q=(max(q1,q2,q3,q4)+max_Q)/2.0 431 517 # Calculate the minimum in absolute terms of 432 518 # the requsted flow and the possible flow 433 519 Q_reduced = sign(Q)*min(abs(Q), abs(max_Q)) 434 520 if self.verbose is True: 521 msg = 'Initial Q Reduced = %.2f m3/s. ' % (Q_reduced) 522 log.critical(msg) 523 524 if self.log_filename is not None: 525 log_to_file(self.log_filename, msg) 526 # Now Keep Rolling Average of Computed Discharge to Reduce / Remove Oscillations 527 # can use delta_t if we want to averageover a time frame for example 528 # N = 5.0/delta_t Will provide the average over 5 seconds 529 530 self.i=(self.i+1)%self.N 531 self.Q_list[self.i]=Q_reduced 532 Q_reduced = sum(self.Q_list)/len(self.Q_list) 533 534 if self.verbose is True: 535 msg = 'Final Q Reduced = %.2f m3/s. ' % (Q_reduced) 536 log.critical(msg) 537 538 if self.log_filename is not None: 539 log_to_file(self.log_filename, msg) 540 541 435 542 if abs(Q_reduced) < abs(Q): 436 543 msg = '%.2fs: Requested flow is ' % time 437 544 msg += 'greater than what is supported by the smallest ' 438 545 msg += 'depth at inlet exchange area:\n ' 546 msg += 'inlet exchange area: %.2f '% (I.exchange_area) 547 msg += 'velocity at inlet :%.2f '% (I.velocity) 548 msg += 'Vel* Exch Area = : %.2f '% (I.velocity*avg_depth*self.width) 439 549 msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\ 440 % ( min_depth, I.exchange_area, delta_t)550 % (avg_depth, I.exchange_area, delta_t) 441 551 msg += ' = %.2f m^3/s\n ' % Q_reduced 442 552 msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced) 553 msg += 'Note calculate max_Q from V %.2f m^3/s ' % (max_Q) 443 554 if self.verbose is True: 444 555 log.critical(msg) … … 473 584 474 585 stage = dq['stage'].get_values(location='centroids', 475 586 indices=[idx])[0] 476 587 depth = h = stage-opening.elevation 477 588 … … 482 593 ymomentum = dq['xmomentum'].get_values(location='centroids', 483 594 indices=[idx])[0] 484 595 485 596 if h > minimum_allowed_height: 486 597 u = xmomentum/(h + velocity_protection/h) … … 509 620 inlet = openings[0] 510 621 outlet = openings[1] 622 623 # FIXME: I think this whole momentum jet thing could be a bit more elegant 624 inlet.momentum = self.opening_momentum[0] 625 outlet.momentum = self.opening_momentum[1] 511 626 else: 512 627 inlet = openings[1] 513 628 outlet = openings[0] 514 629 630 inlet.momentum = self.opening_momentum[1] 631 outlet.momentum = self.opening_momentum[0] 632 515 633 delta_total_energy = -delta_total_energy 516 634 … … 604 722 # Adjust discharge for multiple barrels 605 723 Q *= self.number_of_barrels 606 607 724 725 # Adjust discharge for available water at the inlet 608 726 Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t) 609 727 … … 611 729 self.outlet.rate = Q 612 730 731 732 # Momentum jet stuff 733 if self.use_momentum_jet is True: 734 735 736 # Compute barrel momentum 737 barrel_momentum = barrel_velocity*culvert_outlet_depth 738 739 if self.log_filename is not None: 740 s = 'Barrel velocity = %f' %barrel_velocity 741 log_to_file(self.log_filename, s) 742 743 # Compute momentum vector at outlet 744 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum 745 746 if self.log_filename is not None: 747 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) 748 log_to_file(self.log_filename, s) 749 750 751 # Update momentum 752 if delta_t > 0.0: 753 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value 754 xmomentum_rate /= delta_t 755 756 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value 757 ymomentum_rate /= delta_t 758 759 if self.log_filename is not None: 760 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate) 761 log_to_file(self.log_filename, s) 762 else: 763 xmomentum_rate = ymomentum_rate = 0.0 764 765 766 # Set momentum rates for outlet jet 767 outlet.momentum[0].rate = xmomentum_rate 768 outlet.momentum[1].rate = ymomentum_rate 769 770 # Remember this value for next step (IMPORTANT) 771 outlet.momentum[0].value = outlet_mom_x 772 outlet.momentum[1].value = outlet_mom_y 773 774 if int(domain.time*100) % 100 == 0: 775 776 if self.log_filename is not None: 777 s = 'T=%.5f, Culvert Discharge = %.3f f'\ 778 %(time, Q) 779 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\ 780 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y) 781 s += ' Momentum rate: (%.4f, %.4f)'\ 782 %(xmomentum_rate, ymomentum_rate) 783 s+='Outlet Vel= %.3f'\ 784 %(barrel_velocity) 785 log_to_file(self.log_filename, s) 786 787 788 # Execute momentum terms 789 # This is where Inflow objects are evaluated and update the domain 790 self.outlet.momentum[0](domain) 791 self.outlet.momentum[1](domain) 792 793 794 613 795 # Log timeseries to file 614 796 try: … … 619 801 fid.write('%.2f, %.2f\n' %(time, Q)) 620 802 fid.close() 803 804 # Store value of time 805 self.last_time = time 806 807 808 809 810 621 811 622 812 … … 1194 1384 1195 1385 # Create objects to update momentum (a bit crude at this stage) 1196 1197 1198 1386 xmom0 = General_forcing(domain, 'xmomentum', 1199 1387 polygon=P['exchange_polygon0']) … … 1390 1578 1391 1579 # Update momentum 1392 1393 1580 if delta_t > 0.0: 1394 1581 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
Note: See TracChangeset
for help on using the changeset viewer.