Changeset 7495


Ignore:
Timestamp:
Sep 7, 2009, 6:47:36 PM (11 years ago)
Author:
ole
Message:

Work on culvert flow adjustment and momentum jets

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r7323 r7495  
    136136                 length=None,
    137137                 number_of_barrels=1,
     138                 number_of_smoothing_steps=2000,
    138139                 trigger_depth=0.01, # Depth below which no flow happens
    139140                 manning=None,          # Mannings Roughness for Culvert
     
    152153        # Input check
    153154       
    154         if height is None: height = width       
    155         self.height = height
    156         self.width = width
    157 
     155        if height is None: height = width
    158156       
    159157        assert number_of_barrels >= 1
    160158        assert use_velocity_head is True or use_velocity_head is False
    161159       
    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 
    165164        self.culvert_routine = culvert_routine       
    166165        self.culvert_description_filename = culvert_description_filename
     
    168167            label, type, width, height, length, number_of_barrels, description, rating_curve = read_culvert_description(culvert_description_filename)
    169168            self.rating_curve = ensure_numeric(rating_curve)           
     169
     170        self.height = height
     171        self.width = width
     172
    170173       
    171174        self.domain = domain
     
    292295        self.verbose = verbose
    293296
    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
    295301       
    296302        # For use with update_interval                       
     
    298304        self.update_interval = update_interval
    299305       
     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
    300322
    301323        # Print some diagnostics to log if requested
     
    402424    def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
    403425        """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
    404440        """
    405441   
     
    423459                                               indices=[idx])       
    424460        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
    431517        # Calculate the minimum in absolute terms of
    432518        # the requsted flow and the possible flow
    433519        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
    435542        if abs(Q_reduced) < abs(Q):
    436543            msg = '%.2fs: Requested flow is ' % time
    437544            msg += 'greater than what is supported by the smallest '
    438545            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)
    439549            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)
    441551            msg += ' = %.2f m^3/s\n        ' % Q_reduced
    442552            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)
    443554            if self.verbose is True:
    444555                log.critical(msg)
     
    473584           
    474585            stage = dq['stage'].get_values(location='centroids',
    475                                                indices=[idx])[0]
     586                                           indices=[idx])[0]
    476587            depth = h = stage-opening.elevation
    477588                                                           
     
    482593            ymomentum = dq['xmomentum'].get_values(location='centroids',
    483594                                                   indices=[idx])[0]
    484            
     595
    485596            if h > minimum_allowed_height:
    486597                u = xmomentum/(h + velocity_protection/h)
     
    509620            inlet = openings[0]
    510621            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]
    511626        else:
    512627            inlet = openings[1]
    513628            outlet = openings[0]
    514629           
     630            inlet.momentum = self.opening_momentum[1]
     631            outlet.momentum = self.opening_momentum[0]
     632
    515633            delta_total_energy = -delta_total_energy
    516634
     
    604722        # Adjust discharge for multiple barrels
    605723        Q *= self.number_of_barrels
    606        
    607 
     724
     725        # Adjust discharge for available water at the inlet
    608726        Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
    609727       
     
    611729        self.outlet.rate = Q
    612730
     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           
    613795        # Log timeseries to file
    614796        try:
     
    619801            fid.write('%.2f, %.2f\n' %(time, Q))
    620802            fid.close()
     803
     804        # Store value of time
     805        self.last_time = time
     806
     807
     808           
     809
     810
    621811
    622812                           
     
    11941384       
    11951385        # Create objects to update momentum (a bit crude at this stage)
    1196 
    1197        
    11981386        xmom0 = General_forcing(domain, 'xmomentum',
    11991387                                polygon=P['exchange_polygon0'])
     
    13901578
    13911579            # Update momentum       
    1392 
    13931580            if delta_t > 0.0:
    13941581                xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
Note: See TracChangeset for help on using the changeset viewer.