Ignore:
Timestamp:
Aug 25, 2010, 4:35:40 PM (14 years ago)
Author:
steve
Message:

Got a working culvert!

Location:
trunk/anuga_core/source/anuga/structures
Files:
3 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/structures/boyd_box_routine.py

    r7969 r7975  
    99
    1010
    11 def boyd_box(height, width, flow_width, inflow_specific_energy):
     11def boyd_box(in):
    1212    """Boyd's generalisation of the US department of transportation culvert methods
    1313

  • trunk/anuga_core/source/anuga/structures/culvert_operator.py

    r7974 r7975  
    1313import anuga.utilities.log as log
    1414
    15 import Inlet
     15import inlet
    1616
    1717import numpy as num
    … …  
    6868        enquiry_pt0 = self.enquiry_points[0]
    6969        inlet0_vector = self.culvert_vector
    70         self.inlets.append(Inlet.Inlet(self.domain, polygon0, enquiry_pt0, inlet0_vector))
     70        self.inlets.append(inlet.Inlet(self.domain, polygon0, enquiry_pt0, inlet0_vector))
    7171
    7272        polygon1 = self.inlet_polygons[1]
    7373        enquiry_pt1 = self.enquiry_points[1]
    7474        inlet1_vector = - self.culvert_vector
    75         self.inlets.append(Inlet.Inlet(self.domain, polygon1, enquiry_pt1, inlet1_vector))
     75        self.inlets.append(inlet.Inlet(self.domain, polygon1, enquiry_pt1, inlet1_vector))
    7676 
    77         self.areas = self.domain.areas
    78         self.stages = self.domain.quantities['stage'].centroid_values
    79         self.elevations = self.domain.quantities['elevation'].centroid_values
    80         self.xmoms = self.domain.quantities['xmomentum'].centroid_values
    81         self.ymoms = self.domain.quantities['ymomentum'].centroid_values
     77
    8278   
    8379        self.print_stats()
    … …  
    9086        timestep = self.domain.get_timestep()
    9187
    92         inlet0 = self.inlets[0]
    93         inlet1 = self.inlets[1]
    94 
    95         # Aliases to cell indices
    96         inlet0_indices = inlet0.triangle_indices
    97         inlet1_indices = inlet1.triangle_indices
    98 
    99         # Inlet0 averages
    100         inlet0_heights = inlet0.get_heights()
    101         inlet0_areas = inlet0.get_areas()
    102         inlet0_water_volume = num.sum(inlet0_heights*inlet0_areas)
    103         average_inlet0_height = inlet0_water_volume/inlet0.area
    104 
    105         # Inlet1 averages
    106         inlet1_heights = inlet1.get_heights()
    107         inlet1_areas = inlet1.get_areas()
    108         inlet1_water_volume = num.sum(inlet1_heights*inlet1_areas)
    109         average_inlet1_height = inlet1_water_volume/inlet1.area
     88
     89
     90        inflow  = self.inlets[0]
     91        outflow = self.inlets[1]
     92
     93        # Determine flow direction based on total energy difference
     94        delta_total_energy = inflow.get_average_total_energy() - outflow.get_average_total_energy()
     95
     96        if delta_total_energy < 0:
     97            inflow  = self.inlets[1]
     98            outflow = self.inlets[0]
     99            delta_total_energy = -delta_total_energy
     100
     101        delta_z = inflow.get_average_elevation() - outflow.get_average_elevation()
     102        culvert_slope = delta_z/self.culvert_length
     103
     104        # Determine controlling energy (driving head) for culvert
     105        if inflow.get_average_specific_energy() > delta_total_energy:
     106            # Outlet control
     107            driving_head = delta_total_energy
     108        else:
     109            # Inlet control
     110            driving_head = inflow.get_average_specific_energy()
     111           
    110112
    111113        # Transfer
    112         transfer_water = timestep*inlet0_water_volume
    113        
    114         self.stages.put(inlet0_indices, inlet0.get_elevations() + average_inlet0_height - transfer_water)
    115         self.xmoms.put(inlet0_indices, 0.0)
    116         self.ymoms.put(inlet0_indices, 0.0)
    117 
    118         self.stages.put(inlet1_indices, inlet1.get_elevations() + average_inlet1_height + transfer_water)
    119         self.xmoms.put(inlet1_indices, 0.0)
    120         self.ymoms.put(inlet1_indices, 0.0)
     114        from culvert_routines import boyd_generalised_culvert_model
     115        Q, barrel_velocity, culvert_outlet_depth =\
     116                    boyd_generalised_culvert_model(inflow.get_average_height(),
     117                                         outflow.get_average_height(),
     118                                         inflow.get_average_speed(),
     119                                         outflow.get_average_speed(),
     120                                         inflow.get_average_specific_energy(),
     121                                         delta_total_energy,
     122                                         g,
     123                                         culvert_length=self.culvert_length,
     124                                         culvert_width=self.width,
     125                                         culvert_height=self.height,
     126                                         culvert_type='box',
     127                                         manning=0.01)
     128
     129        transfer_water = Q*timestep
     130
     131
     132        inflow.set_heights(inflow.get_average_height() - transfer_water)
     133        inflow.set_xmoms(0.0)
     134        inflow.set_ymoms(0.0)
     135
     136
     137        outflow.set_heights(outflow.get_average_height() + transfer_water)
     138        outflow.set_xmoms(0.0)
     139        outflow.set_ymoms(0.0)
     140
    121141
    122142
    … …  
    147167
    148168
    149     def set_store_hydrograph_discharge(self, filename=None):
    150 
    151         if filename is None:
    152             self.filename = 'culvert_discharge_hydrograph'
    153         else:
    154             self.filename = filename
    155 
    156         self.discharge_hydrograph = True
    157        
    158         self.timeseries_filename = self.filename + '_timeseries.csv'
    159         fid = open(self.timeseries_filename, 'w')
    160         fid.write('time, discharge\n')
    161         fid.close()
     169
    162170
    163171
    … …  
    223231
    224232   
    225     def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
    226         """Adjust Q downwards depending on available water at inlet
    227 
    228            This is a critical step in modelling bridges and Culverts
    229            the predicted flow through a structure based on an abstract
    230            algorithm can at times request for water that is simply not
    231            available due to any number of constrictions that limit the
    232            flow approaching the structure In order to ensure that
    233            there is adequate flow available certain checks are
    234            required There needs to be a check using the Static Water
    235            Volume sitting infront of the structure, In addition if the
    236            water is moving the available water will be larger than the
    237            static volume
    238            
    239            NOTE To temporarily switch this off for Debugging purposes
    240            rem out line in function def compute_rates below
    241         """
    242    
    243         if delta_t < epsilon:
    244             # No need to adjust if time step is very small or zero
    245             # In this case the possible flow will be very large
    246             # anyway.
    247             return Q
    248        
    249         # Short hands
    250         domain = self.domain       
    251         dq = domain.quantities               
    252         time = domain.get_time()
    253         I = self.inlet
    254         idx = I.exchange_indices   
    255 
    256         # Find triangle with the smallest depth
    257         stage = dq['stage'].get_values(location='centroids',
    258                                                indices=[idx])
    259         elevation = dq['elevation'].get_values(location='centroids',
    260                                                indices=[idx])       
    261         depth = stage-elevation
    262         min_depth = min(depth.flat)  # This may lead to errors if edge of area is at a higher level !!!!
    263         avg_depth = mean(depth.flat) # Yes, but this one violates the conservation unit tests
    264 
    265 
    266 
    267         # FIXME (Ole): If you want these, use log.critical() and
    268         # make the statements depend on verbose
    269         #print I.depth
    270         #print I.velocity
    271         #print self.width
    272 
    273         # max_Q Based on Volume Calcs
    274 
    275 
    276         depth_term = min_depth*I.exchange_area/delta_t
    277         if min_depth < 0.2:
    278             # Only add velocity term in shallow waters (< 20 cm)
    279             # This is a little ad hoc, but maybe it is reasonable
    280             velocity_term = self.width*min_depth*I.velocity
    281         else:
    282             velocity_term = 0.0
    283 
    284         # This one takes approaching water into account   
    285         max_Q = max(velocity_term, depth_term)
    286 
    287         # This one preserves Volume
    288         #max_Q = depth_term
    289 
    290 
    291         if self.verbose is True:
    292             log.critical('Max_Q = %f' % max_Q)           
    293             msg = 'Width = %.2fm, Depth at inlet = %.2f m, Velocity = %.2f m/s.      ' % (self.width, I.depth, I.velocity)
    294             msg += 'Max Q = %.2f m^3/s' %(max_Q)
    295             log.critical(msg)
    296 
    297         if self.log_filename is not None:               
    298             log_to_file(self.log_filename, msg)
    299         # New Procedure for assessing the flow available to the Culvert
    300         # This routine uses the GET FLOW THROUGH CROSS SECTION
    301         #   Need to check Several Polyline however
    302         # Firstly 3 sides of the exchange Poly
    303         # then only the Line Directly infront of the Polygon
    304         # Access polygon Points from   self.inlet.polygon
    305      
    306         #  The Following computes the flow crossing over 3 sides of the exchange polygon for the structure
    307         # Clearly the flow in the culvert can not be more than that flowing toward it through the exhange polygon
    308        
    309         #q1 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][1:3]) # First Side Segment
    310         #q2 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][2:])   # Second Face Segment
    311         #q3 =domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'].take([3,0], axis=0)) # Third Side Segment
    312         # q4 = domain.get_flow_through_cross_section([self.culvert_polygons['exchange_polygon0'][1:4]][0])
    313         #q4=max(q1,0.0)+max(q2,0.0)+max(q3,0.0)
    314         # To use only the Flow crossing the 3 sides of the Exchange Polygon use the following Line Only
    315         #max_Q=max(q1,q2,q3,q4)
    316         # Try Simple Smoothing using Average of 2 approaches
    317         #max_Q=(max(q1,q2,q3,q4)+max_Q)/2.0
    318         # Calculate the minimum in absolute terms of
    319         # the requsted flow and the possible flow
    320         Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
    321         if self.verbose is True:
    322             msg = 'Initial Q Reduced = %.2f m3/s.      ' % (Q_reduced)
    323             log.critical(msg)
    324 
    325         if self.log_filename is not None:               
    326             log_to_file(self.log_filename, msg)
    327         # Now Keep Rolling Average of Computed Discharge to Reduce / Remove Oscillations
    328         #  can use delta_t if we want to averageover a time frame for example
    329         # N = 5.0/delta_t  Will provide the average over 5 seconds
    330 
    331         self.i=(self.i+1)%self.N
    332         self.Q_list[self.i]=Q_reduced
    333         Q_reduced = sum(self.Q_list)/len(self.Q_list)
    334 
    335         if self.verbose is True:
    336             msg = 'Final Q Reduced = %.2f m3/s.      ' % (Q_reduced)
    337             log.critical(msg)
    338 
    339         if self.log_filename is not None:               
    340             log_to_file(self.log_filename, msg)
    341 
    342 
    343         if abs(Q_reduced) < abs(Q):
    344             msg = '%.2fs: Requested flow is ' % time
    345             msg += 'greater than what is supported by the smallest '
    346             msg += 'depth at inlet exchange area:\n        '
    347             msg += 'inlet exchange area: %.2f '% (I.exchange_area)
    348             msg += 'velocity at inlet :%.2f '% (I.velocity)
    349             msg += 'Vel* Exch Area = : %.2f '% (I.velocity*avg_depth*self.width)
    350             msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\
    351                 % (avg_depth, I.exchange_area, delta_t)
    352             msg += ' = %.2f m^3/s\n        ' % Q_reduced
    353             msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
    354             msg += 'Note calculate max_Q from V %.2f m^3/s ' % (max_Q)
    355             if self.verbose is True:
    356                 log.critical(msg)
    357                
    358             if self.log_filename is not None:               
    359                 log_to_file(self.log_filename, msg)
    360        
    361         return Q_reduced   
     233
    362234                       
    363            
    364     def compute_rates(self, delta_t):
    365         """Compute new rates for inlet and outlet
    366         """
    367         # Time stuff
    368         time = self.domain.get_time()
    369         self.last_update = time
    370 
    371         if hasattr(self, 'log_filename'):
    372             log_filename = self.log_filename
    373            
    374         # Compute stage, energy and velocity at the
    375         # enquiry points at each end of the culvert
    376         total_energies = []
    377         specific_energies = []
    378         velocities = []
    379        
    380         for i, inlet in enumerate(self.inlets):
    381            
    382             stage = inlet.get_average_stage()
    383             depth = stage - inlet.get_average_elevation()
    384                                                                        
    385             if depth > minimum_allowed_height:
    386                 u, v = inlet.get_average_velocities
    387             else:
    388                 u = 0.0
    389                 v = 0.0
    390                
    391             uv2 =  u**2 + v**2   
    392                
    393             if self.use_velocity_head is True:
    394                 velocity_head = 0.5*uv2/g   
    395             else:
    396                 velocity_head = 0.0
    397            
    398             inlet.set_total_energy(velocity_head + stage)
    399             inlet.set_specific_energy(velocity_head + depth)
    400 
    401         # We now need to deal with each opening individually
    402         # Determine flow direction based on total energy difference
    403         delta_total_energy = inlets[0].get_total_energy() - inlets[1].get_total_energy()
    404        
    405         if delta_total_energy >= 0:
    406             inlet = inlets[0]
    407             outlet = inlets[1]
    408         else:
    409             inlet = inlets[1]
    410             outlet = inlets[2]
    411            
    412             #msg = 'Total energy difference is negative'
    413             #assert delta_total_energy >= 0.0, msg
    414 
    415         # Recompute slope and issue warning if flow is uphill
    416         # These values do not enter the computation
    417         delta_z = inlet.get_average_elevation() - outlet.get_average_elevation()
    418         culvert_slope = (delta_z/self.culvert_length)
    419         if culvert_slope < 0.0:
    420             # Adverse gradient - flow is running uphill
    421             # Flow will be purely controlled by uphill outlet face
    422             if self.verbose is True:
    423                 log.critical('%.2fs - WARNING: Flow is running uphill.' % time)
    424            
    425         if self.log_filename is not None:
    426             s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
    427                 %(time, self.inlet.stage, self.outlet.stage)
    428             log_to_file(self.log_filename, s)
    429             s = 'Delta total energy = %.3f' %(delta_total_energy)
    430             log_to_file(log_filename, s)
    431            
    432         # Determine controlling energy (driving head) for culvert
    433         if inlet.get_specific_energy() > delta_total_energy:
    434             # Outlet control
    435             driving_head = delta_total_energy
    436         else:
    437             # Inlet control
    438             driving_head = inlet.get_specific_energy()
    439            
    440         inlet_depth = inlet.get_average_stage() - inlet.get_average_elevation()
    441         outlet_depth = outlet.get_average_stage() - outlet.get_average_elevation()
    442        
    443         if inlet_depth <= self.trigger_depth:
    444             Q = 0.0
    445         else:
    446             # Calculate discharge for one barrel and
    447             # set inlet.rate and outlet.rate
    448            
    449             if self.culvert_description_filename is not None:
    450                 try:
    451                     Q = interpolate_linearly(driving_head,
    452                                              self.rating_curve[:,0],
    453                                              self.rating_curve[:,1])
    454                 except Below_interval, e:
    455                     Q = self.rating_curve[0,1]             
    456                     msg = '%.2fs: ' % time
    457                     msg += 'Delta head smaller than rating curve minimum: '
    458                     msg += str(e)
    459                     msg += '\n        '
    460                     msg += 'I will use minimum discharge %.2f m^3/s ' % Q
    461                     msg += 'for culvert "%s"' % self.label
    462                    
    463                     if hasattr(self, 'log_filename'):                   
    464                         log_to_file(self.log_filename, msg)
    465                 except Above_interval, e:
    466                     Q = self.rating_curve[-1,1]         
    467                     msg = '%.2fs: ' % time                 
    468                     msg += 'Delta head greater than rating curve maximum: '
    469                     msg += str(e)
    470                     msg += '\n        '
    471                     msg += 'I will use maximum discharge %.2f m^3/s ' % Q
    472                     msg += 'for culvert "%s"' % self.label
    473                    
    474                     if self.log_filename is not None:                   
    475                         log_to_file(self.log_filename, msg)
    476             else:
    477                 # User culvert routine
    478                 Q, barrel_velocity, culvert_outlet_depth =\
    479                     self.culvert_routine(inlet_depth,
    480                                          outlet_depth,
    481                                          inlet.get_average_velocity(),
    482                                          outlet.get_average_velocity(),
    483                                          inlet.get_specific_energy(),
    484                                          delta_total_energy,
    485                                          g,
    486                                          culvert_length=self.length,
    487                                          culvert_width=self.width,
    488                                          culvert_height=self.height,
    489                                          culvert_type=self.culvert_type,
    490                                          manning=self.manning,
    491                                          sum_loss=self.sum_loss,
    492                                          log_filename=self.log_filename)
    493            
    494         # Adjust discharge for multiple barrels
    495         Q *= self.number_of_barrels
    496 
    497         # Adjust discharge for available water at the inlet
    498         Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
    499        
    500         self.inlet.rate = -Q
    501         self.outlet.rate = Q
    502 
    503         # Momentum jet stuff
    504         if self.use_momentum_jet is True:
    505 
    506             # Compute barrel momentum
    507             barrel_momentum = barrel_velocity*culvert_outlet_depth
    508 
    509             if self.log_filename is not None:                                   
    510                 s = 'Barrel velocity = %f' %barrel_velocity
    511                 log_to_file(self.log_filename, s)
    512 
    513             # Compute momentum vector at outlet
    514             outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
    515                
    516             if self.log_filename is not None:               
    517                 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
    518                 log_to_file(self.log_filename, s)
    519 
    520             # Update momentum       
    521             if delta_t > 0.0:
    522                 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
    523                 xmomentum_rate /= delta_t
    524                    
    525                 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
    526                 ymomentum_rate /= delta_t
    527                        
    528                 if self.log_filename is not None:               
    529                     s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
    530                     log_to_file(self.log_filename, s)                   
    531             else:
    532                 xmomentum_rate = ymomentum_rate = 0.0
    533 
    534 
    535             # Set momentum rates for outlet jet
    536             outlet.momentum[0].rate = xmomentum_rate
    537             outlet.momentum[1].rate = ymomentum_rate
    538 
    539             # Remember this value for next step (IMPORTANT)
    540             outlet.momentum[0].value = outlet_mom_x
    541             outlet.momentum[1].value = outlet_mom_y                   
    542 
    543             if int(domain.time*100) % 100 == 0:
    544 
    545                 if self.log_filename is not None:               
    546                     s = 'T=%.5f, Culvert Discharge = %.3f f'\
    547                         %(time, Q)
    548                     s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
    549                         %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
    550                     s += ' Momentum rate: (%.4f, %.4f)'\
    551                         %(xmomentum_rate, ymomentum_rate)                   
    552                     s+='Outlet Vel= %.3f'\
    553                         %(barrel_velocity)
    554                     log_to_file(self.log_filename, s)
    555 
    556 
    557             # Execute momentum terms
    558             # This is where Inflow objects are evaluated and update the domain
    559                 self.outlet.momentum[0](domain)
    560                 self.outlet.momentum[1](domain)       
    561            
    562 
    563            
    564         # Log timeseries to file
    565         try:
    566             fid = open(self.timeseries_filename, 'a')       
    567         except:
    568             pass
    569         else:   
    570             fid.write('%.2f, %.2f\n' %(time, Q))
    571             fid.close()
    572 
    573         # Store value of time
    574         self.last_time = time
    575 
    576 
    577235# FIXME(Ole): Write in C and reuse this function by similar code
    578236# in interpolate.py

  • trunk/anuga_core/source/anuga/structures/inlet.py

    r7974 r7975  
     1
    12from anuga.geometry.polygon import inside_polygon, is_inside_polygon
    2 from anuga.config import velocity_protection 
     3from anuga.config import velocity_protection, g
    34import math
    45
    … …  
    117118    def get_total_water_volume(self):
    118119       
    119        return num.sum(self.get_heights*self.get_areas())
     120       return num.sum(self.get_heights()*self.get_areas())
    120121   
    121122   
    … …  
    134135           
    135136           
    136     def get_average_velocity(self):
     137    def get_average_speed(self):
    137138 
    138139            u, v = self.get_velocities()
    … …  
    141142            average_v = num.sum(v)/self.triangle_indices.size
    142143           
    143             return math.sqrt(average_u**2 + average_v**2) 
     144            return math.sqrt(average_u**2 + average_v**2)
    144145
    145146
    146     def set_total_energy(self, total_energy):
     147
     148    def get_average_velocity_head(self):
     149
     150        return 0.5*self.get_average_speed()**2/g
     151
     152
     153    def get_average_total_energy(self):
    147154       
    148         self.total_energy = total_energy
     155        return self.get_average_velocity_head() + self.get_average_stage()
    149156       
     157   
     158    def get_average_specific_energy(self):
    150159       
    151     def get_total_energy(self):
    152        
    153         return self.total_energy
    154        
    155        
    156     def set_specific_energy(self, specific_energy):
    157        
    158         self.specific_energy = specific_energy
     160        return self.get_average_velocity_head() + self.get_average_height()
    159161
    160    
    161     def get_specific_energy(self):
    162        
    163         return self.specific_energy
    164        
    165        
     162
     163
     164    def set_heights(self,height):
     165
     166        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + height)
     167
     168
     169    def set_stages(self,stage):
     170
     171        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage)
     172
     173    def set_xmoms(self,xmom):
     174
     175        self.xmoms=self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom)
     176
     177
     178    def set_ymoms(self,ymom):
     179
     180        self.xmoms=self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom)
     181
     182
     183    def set_elevations(self,elevation):
     184
     185        self.domain.quantities['elevation'].centroid_values.put(self.triangle_indices, elevation)

  • trunk/anuga_core/source/anuga/structures/testing_culvert.py

    r7968 r7975  
    4242domain.set_name('Test_culvert')                 # Output name
    4343domain.set_default_order(2)
     44#domain.set_beta(1.5)
    4445
    4546
    … …  
    110111
    111112## Inflow based on Flow Depth and Approaching Momentum
    112 Bi = anuga.Dirichlet_boundary([1.0, 0.0, 0.0])
     113Bi = anuga.Dirichlet_boundary([2.0, 0.0, 0.0])
    113114Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
    114115#Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
    … …  
    131132#min_delta_w = sys.maxint
    132133#max_delta_w = -min_delta_w
    133 for t in domain.evolve(yieldstep = 1.0, finaltime = 100):
     134for t in domain.evolve(yieldstep = 1.0, finaltime = 200):
    134135    domain.write_time()
     136
     137    if domain.get_time() > 150.5 and domain.get_time() < 151.5 :
     138        Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
     139        domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br})
     140
    135141    #delta_w = culvert.inlet.stage - culvert.outlet.stage
    136142   
Note: See TracChangeset for help on using the changeset viewer.