Changeset 5434


Ignore:
Timestamp:
Jun 25, 2008, 5:01:54 PM (14 years ago)
Author:
ole
Message:

Added culvert framework

File:
1 edited

Legend:

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

    r5433 r5434  
     1 
     2class 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    top-down_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 = time-self.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 (x-direction)
     220            uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
     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.