Changeset 5437


Ignore:
Timestamp:
Jun 25, 2008, 7:10:21 PM (14 years ago)
Author:
ole
Message:

Finished culvert flow work with a good working solution

Files:
3 edited

Legend:

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

    r5435 r5437  
    9494        if blockage_bottup is None: blockage_bottup=0.00
    9595        if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model
    96         if label is None: label = 'culvert_flow_' + id(self)
     96        if label is None: label = 'culvert_flow'
     97        label += '_' + str(id(self))
    9798       
    9899        # Open log file for storing some specific results...
     
    100101        self.label = label
    101102
    102 
    103         # Create the fundamental culvert polygons from POLYGON
     103        # Print something
     104        log_to_file(self.log_filename, self.label)       
     105        log_to_file(self.log_filename, description)
     106        log_to_file(self.log_filename, self.culvert_type)       
     107
     108
     109        # Create the fundamental culvert polygons from POLYGON
     110        if self.culvert_type == 'circle':
     111            # Redefine width and height for use with create_culvert_polygons
     112            width = height = diameter
     113       
    104114        P = create_culvert_polygons(end_point0,
    105115                                    end_point1,
     
    135145        self.verbose = verbose
    136146        self.last_time = 0.0       
    137         self.temp_keep_delta_t = 0.0
    138147       
    139148
     
    239248                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
    240249                z = mean(elevation)                   # Average elevation
    241                
     250
    242251            # Estimated depth above the culvert inlet
    243             d = w - z
    244 
     252            d = w - z  # Used for calculations involving depth
    245253            if d < 0.0:
    246254                # This is possible since w and z are taken at different locations
     
    249257                d = 0.0
    250258           
     259
     260
     261            # Depth at exchange area used to trigger calculations
     262            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
     263            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
     264            depth = stage - elevation
     265            d_trigger = mean(depth)
     266
     267
     268
    251269            # Ratio of depth to culvert height.
    252270            # If ratio > 1 then culvert is running full
    253             ratio = d/self.height 
     271            if self.culvert_type == 'circle':
     272                ratio = d/self.diameter
     273            else:   
     274                ratio = d/self.height 
    254275            opening.ratio = ratio
    255276               
     
    260281            opening.specific_energy = Es           
    261282           
    262             # Store current average stage and depth with each opening object 
     283            # Store current average stage and depth with each opening object
    263284            opening.depth = d
     285            opening.depth_trigger = d_trigger           
    264286            opening.stage = w
    265287            opening.elevation = z
  • anuga_core/source/anuga/culvert_flows/culvert_routines.py

    r5436 r5437  
    5858    sum_loss = culvert.sum_loss
    5959    length = culvert.length
    60    
    61     if inlet.depth >= 0.01:
     60
     61    if inlet.depth_trigger >= 0.01 and inlet.depth >= 0.01:
    6262        # Calculate driving energy
    6363        E = inlet.specific_energy
     
    9191                outlet_culvert_depth = inlet.depth
    9292                width = diameter*sin(alpha)
    93                 perimeter = alpha*diameter               
     93                #perimeter = alpha*diameter               
    9494                case = 'Inlet unsubmerged'
    9595            else:   
     
    9898                outlet_culvert_depth = diameter
    9999                width = diameter
    100                 perimeter = diameter
     100                #perimeter = diameter
    101101                case = 'Inlet submerged'                   
    102102
     
    129129                    case = 'Outlet is open channel flow'
    130130
     131                hyd_rad = flow_area/perimeter
     132                s = 'hydraulic radius at outlet = %f' %hyd_rad
     133                log_to_file(log_filename, s)
     134
     135                # Outlet control velocity using tail water
     136                culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
     137                Q_outlet_tailwater = flow_area * culvert_velocity
     138
     139                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     140                log_to_file(log_filename, s)
     141                Q = min(Q, Q_outlet_tailwater)
     142                   
     143
    131144
    132145        else:
     
    148161                flow_area = width*inlet.depth
    149162                outlet_culvert_depth = inlet.depth
    150                 perimeter=(width+2.0*inlet.depth)               
     163                #perimeter=(width+2.0*inlet.depth)               
    151164                case = 'Inlet unsubmerged'
    152165            else:   
     
    154167                flow_area = width*height
    155168                outlet_culvert_depth = height
    156                 perimeter=2.0*(width+height)               
     169                #perimeter=2.0*(width+height)               
    157170                case = 'Inlet submerged'                   
    158171
     
    176189                    perimeter=(width+2.0*outlet.depth)
    177190                    case = 'Outlet is open channel flow'
    178                    
    179 
    180 
    181 
    182         # Common code for rectangular and circular culvert types   
    183         hyd_rad = flow_area/perimeter
    184         s = 'hydraulic radius at outlet = %f' %hyd_rad
    185         log_to_file(log_filename, s)
    186 
    187         # Outlet control velocity using tail water
    188         culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
    189         Q_outlet_tailwater = flow_area * culvert_velocity
    190 
    191         s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    192         log_to_file(log_filename, s)
    193         Q = min(Q, Q_outlet_tailwater)
    194 
     191
     192                hyd_rad = flow_area/perimeter
     193                s = 'hydraulic radius at outlet = %f' %hyd_rad
     194                log_to_file(log_filename, s)
     195
     196                # Outlet control velocity using tail water
     197                culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
     198                Q_outlet_tailwater = flow_area * culvert_velocity
     199
     200                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     201                log_to_file(log_filename, s)
     202                Q = min(Q, Q_outlet_tailwater)
     203
     204
     205        # Common code for circle and square geometries
    195206        log_to_file(log_filename, 'Case: "%s"' %case)
    196    
    197207        flow_rate_control=Q
    198208
     
    209219        # Determine momentum at the outlet
    210220        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
     221
     222
    211223    else: #inlet.depth < 0.01:
    212224        Q = barrel_velocity = outlet_culvert_depth = 0.0
  • anuga_work/development/culvert_flow/culvert_boyd_channel.py

    r5436 r5437  
    219219    # Close
    220220
    221 fid.close()
    222221
    223222
Note: See TracChangeset for help on using the changeset viewer.