Changeset 6299


Ignore:
Timestamp:
Feb 10, 2009, 5:31:46 AM (16 years ago)
Author:
ole
Message:

Cleaned up culvert routine and made all input parameters explicit.

Location:
anuga_core/source/anuga/culvert_flows
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r6150 r6299  
    204204            log_to_file(self.log_filename, description)
    205205            log_to_file(self.log_filename, self.culvert_type)       
     206        else:
     207            self.log_filename = None
    206208
    207209
    … …  
    296298
    297299        # Print some diagnostics to log if requested
    298         if hasattr(self, 'log_filename'):
     300        if self.log_filename is not None:
    299301            s = 'Culvert Effective Length = %.2f m' %(self.length)
    300302            log_to_file(self.log_filename, s)
    … …  
    323325            if time - self.last_update > self.update_interval or time == 0.0:
    324326                update = True
    325 
    326         if hasattr(self, 'log_filename'):           
     327           
     328        if self.log_filename is not None:       
    327329            s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    328330            log_to_file(self.log_filename, s)
    … …  
    439441            if self.verbose is True:
    440442                print msg
    441             if hasattr(self, 'log_filename'):                   
     443               
     444            if self.log_filename is not None:               
    442445                log_to_file(self.log_filename, msg)
    443446       
    … …  
    524527                print '%.2fs - WARNING: Flow is running uphill.' % time
    525528           
    526         if hasattr(self, 'log_filename'):
     529        if self.log_filename is not None:
    527530            s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\
    528531                %(time, self.inlet.stage, self.outlet.stage)
    … …  
    573576                    msg += 'for culvert "%s"' % self.label
    574577                   
    575                     if hasattr(self, 'log_filename'):                   
     578                    if self.log_filename is not None:                   
    576579                        log_to_file(self.log_filename, msg)
    577580            else:
    578581                # User culvert routine
    579582                Q, barrel_velocity, culvert_outlet_depth =\
    580                     self.culvert_routine(self, delta_total_energy, g)
     583                    self.culvert_routine(inlet.depth,
     584                                         outlet.depth,
     585                                         inlet.specific_energy,
     586                                         delta_total_energy,
     587                                         g,
     588                                         culvert_length=self.length,
     589                                         culvert_width=self.width,
     590                                         culvert_height=self.height,
     591                                         culvert_type=self.culvert_type,
     592                                         manning=self.manning,
     593                                         sum_loss=self.sum_loss,
     594                                         log_filename=self.log_filename)
    581595           
    582596           
    … …  
    601615
    602616                           
    603    
     617# OBSOLETE (Except for momentum jet in Culvert_flow_energy)   
    604618class Culvert_flow_rating:
    605619    """Culvert flow - transfer water from one hole to another

  • anuga_core/source/anuga/culvert_flows/culvert_routines.py

    r6128 r6299  
    88
    99#NOTE:
    10 # Inlet control:  Delta_Et > Es at the inlet
    11 # Outlet control: Delta_Et < Es at the inlet
    12 # where Et is total energy (w + 0.5*v^2/g) and
    13 # Es is the specific energy (h + 0.5*v^2/g)
    14 
    15 
    16 
    17 #   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
    18 #
    19 # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
    20 # Flow Is Removed at a Rate of INFLOW
    21 # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
    22 #
    23 # This SHould be changed to a Vertical Opening Both BOX and Circular
    24 # There will be several Culvert Routines such as:
    25 # CULVERT_Boyd_Channel
    26 # CULVERT_Orifice_and_Weir
    27 # CULVERT_Simple_FLOOR
    28 # CULVERT_Simple_WALL
    29 # CULVERT_Eqn_Floor
    30 # CULVERT_Eqn_Wall
    31 # CULVERT_Tab_Floor
    32 # CULVERT_Tab_Wall
    33 # BRIDGES.....
    34 # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
    35 
    36 #  COULD USE EPA SWMM Model !!!!
     10# Inlet control:  Delta_total_energy > inlet_specific_energy
     11# Outlet control: Delta_total_energy < inlet_specific_energy
     12# where total energy is (w + 0.5*v^2/g) and
     13# specific energy is (h + 0.5*v^2/g)
    3714
    3815
    … …  
    4017
    4118
    42 def boyd_generalised_culvert_model(culvert, delta_total_energy, g):
    43 
    44     """Boyd's generalisation of the US department of transportation culvert model
    45         # == The quantity of flow passing through a culvert is controlled by many factors
    46         # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
    47         # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
    48         # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled
     19def boyd_generalised_culvert_model(inlet_depth,
     20                                   outlet_depth,
     21                                   inlet_specific_energy,
     22                                   delta_total_energy,
     23                                   g,
     24                                   culvert_length=0.0,
     25                                   culvert_width=0.0,
     26                                   culvert_height=0.0,
     27                                   culvert_type='box',
     28                                   manning=0.0,
     29                                   sum_loss=0.0,
     30                                   log_filename=None):
     31
     32    """Boyd's generalisation of the US department of transportation culvert
     33    model
     34     
     35    The quantity of flow passing through a culvert is controlled by many factors
     36    It could be that the culvert is controlled by the inlet only, with it
     37    being unsubmerged this is effectively equivalent to the weir Equation
     38    Else the culvert could be controlled by the inlet, with it being
     39    submerged, this is effectively the Orifice Equation
     40    Else it may be controlled by down stream conditions where depending on
     41    the down stream depth, the momentum in the culvert etc. flow is controlled
    4942    """
    5043
    … …  
    5346    from anuga.utilities.numerical_tools import safe_acos as acos
    5447
    55     inlet = culvert.inlet
    56     outlet = culvert.outlet       
    57     Q_outlet_tailwater = 0.0
    58     inlet.rate = 0.0
    59     outlet.rate = 0.0
    60     Q_inlet_unsubmerged = 0.0
    61     Q_inlet_submerged = 0.0
    62     Q_outlet_critical_depth = 0.0
    63 
    64     if hasattr(culvert, 'log_filename'):                               
    65         log_filename = culvert.log_filename
    66 
    67     manning = culvert.manning
    68     sum_loss = culvert.sum_loss
    69     length = culvert.length
    70 
    71     if inlet.depth > 0.01:
    72 
    73         E = inlet.specific_energy
    74 
    75         if hasattr(culvert, 'log_filename'):                           
    76             s = 'Specific energy  = %f m' %E
     48
     49    if inlet_depth > 0.01:
     50        # Water has risen above inlet
     51       
     52        if log_filename is not None:                           
     53            s = 'Specific energy  = %f m' % inlet_specific_energy
    7754            log_to_file(log_filename, s)
    7855       
    79         msg = 'Specific energy is negative'
    80         assert E >= 0.0, msg
     56        msg = 'Specific energy at inlet is negative'
     57        assert inlet_specific_energy >= 0.0, msg
    8158                     
    82        
    83         # Water has risen above inlet
    84         if culvert.culvert_type == 'circle':
    85             # Round culvert
     59
     60        if culvert_type == 'circle':
     61            # Round culvert (use width as diameter)
     62            diameter = culvert_width           
    8663
    8764            # Calculate flows for inlet control           
    88             diameter = culvert.diameter
    89 
    90             Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*E**1.63    # Inlet Ctrl Inlet Unsubmerged
    91             Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**0.63      # Inlet Ctrl Inlet Submerged
    92 
    93             if hasattr(culvert, 'log_filename'):
    94                 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
     65            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
     66            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
     67
     68            if log_filename is not None:                                       
     69                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
    9570                log_to_file(log_filename, s)
    9671
    97             case = ''
     72
     73            # FIXME(Ole): Are these functions really for inlet control?                   
    9874            if Q_inlet_unsubmerged < Q_inlet_submerged:
    9975                Q = Q_inlet_unsubmerged
    100 
    101                 alpha = acos(1 - inlet.depth/diameter)
     76                alpha = acos(1 - inlet_depth/diameter)
    10277                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    103                 outlet_culvert_depth = inlet.depth
    104                 width = diameter*sin(alpha)
    105                 #perimeter = alpha*diameter               
     78                outlet_culvert_depth = inlet_depth
    10679                case = 'Inlet unsubmerged'
    10780            else:   
    … …  
    10982                flow_area = (diameter/2)**2 * pi
    11083                outlet_culvert_depth = diameter
    111                 width = diameter
    112                 #perimeter = diameter
    11384                case = 'Inlet submerged'                   
    11485
    11586
    11687
    117             if delta_total_energy < E:
     88            if delta_total_energy < inlet_specific_energy:
    11889                # Calculate flows for outlet control
     90               
    11991                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    120 
    121                 if outlet.depth > diameter:        # The Outlet is Submerged
     92                if outlet_depth > diameter:        # The Outlet is Submerged
    12293                    outlet_culvert_depth=diameter
    12394                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    12495                    perimeter = diameter * pi
    125                     width = diameter
    12696                    case = 'Outlet submerged'
    127                 elif outlet.depth==0.0:
    128                     outlet_culvert_depth=inlet.depth  # For purpose of calculation assume the outlet depth = the inlet depth
    129                     alpha = acos(1 - inlet.depth/diameter)
     97                elif outlet_depth==0.0:
     98                    outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth
     99                    alpha = acos(1 - inlet_depth/diameter)
    130100                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    131101                    perimeter = alpha*diameter
    132                     width = diameter*sin(alpha)                   
    133102
    134103                    case = 'Outlet depth is zero'                       
    135104                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    136                     outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
    137                     alpha = acos(1 - outlet.depth/diameter)
     105                    outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
     106                    alpha = acos(1 - outlet_depth/diameter)
    138107                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    139108                    perimeter = alpha*diameter
    140                     width = diameter*sin(alpha)                   
    141109                    case = 'Outlet is open channel flow'
    142110
    143111                hyd_rad = flow_area/perimeter
    144112               
    145                 if hasattr(culvert, 'log_filename'):
     113                if log_filename is not None:
    146114                    s = 'hydraulic radius at outlet = %f' %hyd_rad
    147115                    log_to_file(log_filename, s)
    148116
    149117                # Outlet control velocity using tail water
    150                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
     118                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    151119                Q_outlet_tailwater = flow_area * culvert_velocity
    152120               
    153                 if hasattr(culvert, 'log_filename'):
     121                if log_filename is not None:               
    154122                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    155123                    log_to_file(log_filename, s)
    … …  
    165133
    166134            # Calculate flows for inlet control
    167             height = culvert.height
    168             width = culvert.width           
     135            height = culvert_height
     136            width = culvert_width           
    169137           
    170             Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    171             Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    172 
    173             if hasattr(culvert, 'log_filename'):
     138            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     139            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
     140
     141            if log_filename is not None:                           
    174142                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
    175143                log_to_file(log_filename, s)
    176144
    177             case = ''
     145
     146            # FIXME(Ole): Are these functions really for inlet control?   
    178147            if Q_inlet_unsubmerged < Q_inlet_submerged:
    179148                Q = Q_inlet_unsubmerged
    180                 flow_area = width*inlet.depth
    181                 outlet_culvert_depth = inlet.depth
    182                 #perimeter=(width+2.0*inlet.depth)               
     149                flow_area = width*inlet_depth
     150                outlet_culvert_depth = inlet_depth
    183151                case = 'Inlet unsubmerged'
    184152            else:   
    … …  
    186154                flow_area = width*height
    187155                outlet_culvert_depth = height
    188                 #perimeter=2.0*(width+height)               
    189156                case = 'Inlet submerged'                   
    190157
    191             if delta_total_energy < E:
     158            if delta_total_energy < inlet_specific_energy:
    192159                # Calculate flows for outlet control
     160               
    193161                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    194 
    195                 if outlet.depth > height:        # The Outlet is Submerged
     162                if outlet_depth > height:        # The Outlet is Submerged
    196163                    outlet_culvert_depth=height
    197164                    flow_area=width*height       # Cross sectional area of flow in the culvert
    198165                    perimeter=2.0*(width+height)
    199166                    case = 'Outlet submerged'
    200                 elif outlet.depth==0.0:
    201                     outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
    202                     flow_area=width*inlet.depth
    203                     perimeter=(width+2.0*inlet.depth)
     167                elif outlet_depth==0.0:
     168                    outlet_culvert_depth=inlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
     169                    flow_area=width*inlet_depth
     170                    perimeter=(width+2.0*inlet_depth)
    204171                    case = 'Outlet depth is zero'                       
    205172                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    206                     outlet_culvert_depth=outlet.depth
    207                     flow_area=width*outlet.depth
    208                     perimeter=(width+2.0*outlet.depth)
     173                    outlet_culvert_depth=outlet_depth
     174                    flow_area=width*outlet_depth
     175                    perimeter=(width+2.0*outlet_depth)
    209176                    case = 'Outlet is open channel flow'
    210177
    211178                hyd_rad = flow_area/perimeter
    212179               
    213                 if hasattr(culvert, 'log_filename'):               
    214                     s = 'hydraulic radius at outlet = %f' %hyd_rad
     180                if log_filename is not None:                               
     181                    s = 'hydraulic radius at outlet = %f' % hyd_rad
    215182                    log_to_file(log_filename, s)
    216183
    217184                # Outlet control velocity using tail water
    218                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
     185                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    219186                Q_outlet_tailwater = flow_area * culvert_velocity
    220187
    221                 if hasattr(culvert, 'log_filename'):                           
    222                     s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     188                if log_filename is not None:                           
     189                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
    223190                    log_to_file(log_filename, s)
    224191                Q = min(Q, Q_outlet_tailwater)
    … …  
    227194                #FIXME(Ole): What about inlet control?
    228195
     196               
    229197        # Common code for circle and square geometries
    230         if hasattr(culvert, 'log_filename'):                                   
    231             log_to_file(log_filename, 'Case: "%s"' %case)
     198        if log_filename is not None:
     199            log_to_file(log_filename, 'Case: "%s"' % case)
    232200           
    233         flow_rate_control=Q
    234 
    235         if hasattr(culvert, 'log_filename'):                                   
    236             s = 'Flow Rate Control = %f' %flow_rate_control
     201        if log_filename is not None:                       
     202            s = 'Flow Rate Control = %f' % Q
    237203            log_to_file(log_filename, s)
    238204
    239         inlet.rate = -flow_rate_control
    240         outlet.rate = flow_rate_control               
    241205           
    242         culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
    243        
    244         if hasattr(culvert, 'log_filename'):                           
    245             s = 'Froude in Culvert = %f' %culv_froude
     206        culv_froude=sqrt(Q**2*width/(g*flow_area**3))
     207       
     208        if log_filename is not None:                           
     209            s = 'Froude in Culvert = %f' % culv_froude
    246210            log_to_file(log_filename, s)
    247211
    … …  
    250214
    251215
    252     else: #inlet.depth < 0.01:
     216    else: # inlet_depth < 0.01:
    253217        Q = barrel_velocity = outlet_culvert_depth = 0.0
    254218       
Note: See TracChangeset for help on using the changeset viewer.