Ignore:
Timestamp:
Mar 17, 2009, 4:02:54 PM (15 years ago)
Author:
rwilson
Message:

Revert back to 6481, prior to auto-merge of trunk and numpy branch.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/numpy/anuga/culvert_flows/culvert_routines.py

    r6517 r6533  
    88
    99#NOTE:
    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)
     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 !!!!
    1437
    1538
    … …  
    1740
    1841
    19 def 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
     42def 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
    4249    """
    4350
    … …  
    4653    from anuga.utilities.numerical_tools import safe_acos as acos
    4754
    48     if inlet_depth > 0.01:
     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
     77            log_to_file(log_filename, s)
     78       
     79        msg = 'Specific energy is negative'
     80        assert E >= 0.0, msg
     81                     
     82       
    4983        # Water has risen above inlet
    50        
    51         if log_filename is not None:                           
    52             s = 'Specific energy  = %f m' % inlet_specific_energy
    53             log_to_file(log_filename, s)
    54        
    55         msg = 'Specific energy at inlet is negative'
    56         assert inlet_specific_energy >= 0.0, msg
    57                      
    58         if culvert_type == 'circle':
    59             # Round culvert (use width as diameter)
    60             diameter = culvert_width           
     84        if culvert.culvert_type == 'circle':
     85            # Round culvert
    6186
    6287            # Calculate flows for inlet control           
    63             Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
    64             Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
    65 
    66             if log_filename is not None:                                       
    67                 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
     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)
    6895                log_to_file(log_filename, s)
    6996
    70             # FIXME(Ole): Are these functions really for inlet control?                   
     97            case = ''
    7198            if Q_inlet_unsubmerged < Q_inlet_submerged:
    7299                Q = Q_inlet_unsubmerged
    73100
    74                 alpha = acos(1 - inlet_depth/diameter)
     101                alpha = acos(1 - inlet.depth/diameter)
    75102                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    76                 outlet_culvert_depth = inlet_depth
     103                outlet_culvert_depth = inlet.depth
     104                width = diameter*sin(alpha)
     105                #perimeter = alpha*diameter               
    77106                case = 'Inlet unsubmerged'
    78107            else:   
    … …  
    80109                flow_area = (diameter/2)**2 * pi
    81110                outlet_culvert_depth = diameter
     111                width = diameter
     112                #perimeter = diameter
    82113                case = 'Inlet submerged'                   
    83114
    84115
    85116
    86             if delta_total_energy < inlet_specific_energy:
     117            if delta_total_energy < E:
    87118                # Calculate flows for outlet control
    88119                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    89120
    90                 if outlet_depth > diameter:        # The Outlet is Submerged
     121                if outlet.depth > diameter:        # The Outlet is Submerged
    91122                    outlet_culvert_depth=diameter
    92123                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    93124                    perimeter = diameter * pi
     125                    width = diameter
    94126                    case = 'Outlet submerged'
    95                 elif outlet_depth==0.0:
    96                     outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth
    97                     alpha = acos(1 - inlet_depth/diameter)
     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)
    98130                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    99131                    perimeter = alpha*diameter
     132                    width = diameter*sin(alpha)                   
    100133
    101134                    case = 'Outlet depth is zero'                       
    102135                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    103                     outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
    104                     alpha = acos(1 - outlet_depth/diameter)
     136                    outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
     137                    alpha = acos(1 - outlet.depth/diameter)
    105138                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    106139                    perimeter = alpha*diameter
     140                    width = diameter*sin(alpha)                   
    107141                    case = 'Outlet is open channel flow'
    108142
    109143                hyd_rad = flow_area/perimeter
    110144               
    111                 if log_filename is not None:
     145                if hasattr(culvert, 'log_filename'):
    112146                    s = 'hydraulic radius at outlet = %f' %hyd_rad
    113147                    log_to_file(log_filename, s)
    114148
    115149                # Outlet control velocity using tail water
    116                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
     150                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
    117151                Q_outlet_tailwater = flow_area * culvert_velocity
    118152               
    119                 if log_filename is not None:               
     153                if hasattr(culvert, 'log_filename'):
    120154                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    121155                    log_to_file(log_filename, s)
    … …  
    131165
    132166            # Calculate flows for inlet control
    133             height = culvert_height
    134             width = culvert_width           
     167            height = culvert.height
     168            width = culvert.width           
    135169           
    136             Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    137             Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    138 
    139             if log_filename is not None:                           
     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'):
    140174                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
    141175                log_to_file(log_filename, s)
    142176
    143             # FIXME(Ole): Are these functions really for inlet control?   
     177            case = ''
    144178            if Q_inlet_unsubmerged < Q_inlet_submerged:
    145179                Q = Q_inlet_unsubmerged
    146                 flow_area = width*inlet_depth
    147                 outlet_culvert_depth = inlet_depth
     180                flow_area = width*inlet.depth
     181                outlet_culvert_depth = inlet.depth
     182                #perimeter=(width+2.0*inlet.depth)               
    148183                case = 'Inlet unsubmerged'
    149184            else:   
    … …  
    154189                case = 'Inlet submerged'                   
    155190
    156             if delta_total_energy < inlet_specific_energy:
     191            if delta_total_energy < E:
    157192                # Calculate flows for outlet control
    158193                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    159194
    160                 if outlet_depth > height:        # The Outlet is Submerged
     195                if outlet.depth > height:        # The Outlet is Submerged
    161196                    outlet_culvert_depth=height
    162197                    flow_area=width*height       # Cross sectional area of flow in the culvert
    163198                    perimeter=2.0*(width+height)
    164199                    case = 'Outlet submerged'
    165                 elif outlet_depth==0.0:
    166                     outlet_culvert_depth=inlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
    167                     flow_area=width*inlet_depth
    168                     perimeter=(width+2.0*inlet_depth)
     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)
    169204                    case = 'Outlet depth is zero'                       
    170205                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    171                     outlet_culvert_depth=outlet_depth
    172                     flow_area=width*outlet_depth
    173                     perimeter=(width+2.0*outlet_depth)
     206                    outlet_culvert_depth=outlet.depth
     207                    flow_area=width*outlet.depth
     208                    perimeter=(width+2.0*outlet.depth)
    174209                    case = 'Outlet is open channel flow'
    175210
    176211                hyd_rad = flow_area/perimeter
    177212               
    178                 if log_filename is not None:                               
    179                     s = 'hydraulic radius at outlet = %f' % hyd_rad
     213                if hasattr(culvert, 'log_filename'):               
     214                    s = 'hydraulic radius at outlet = %f' %hyd_rad
    180215                    log_to_file(log_filename, s)
    181216
    182217                # Outlet control velocity using tail water
    183                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*culvert_length)/hyd_rad**1.33333))
     218                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
    184219                Q_outlet_tailwater = flow_area * culvert_velocity
    185220
    186                 if log_filename is not None:                           
    187                     s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
     221                if hasattr(culvert, 'log_filename'):                           
     222                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    188223                    log_to_file(log_filename, s)
    189224                Q = min(Q, Q_outlet_tailwater)
    … …  
    193228
    194229        # Common code for circle and square geometries
    195         if log_filename is not None:
    196             log_to_file(log_filename, 'Case: "%s"' % case)
     230        if hasattr(culvert, 'log_filename'):                                   
     231            log_to_file(log_filename, 'Case: "%s"' %case)
    197232           
    198         if log_filename is not None:                       
    199             s = 'Flow Rate Control = %f' % Q
     233        flow_rate_control=Q
     234
     235        if hasattr(culvert, 'log_filename'):                                   
     236            s = 'Flow Rate Control = %f' %flow_rate_control
    200237            log_to_file(log_filename, s)
    201238
    202         culv_froude=sqrt(Q**2*width/(g*flow_area**3))
    203        
    204         if log_filename is not None:                           
    205             s = 'Froude in Culvert = %f' % culv_froude
     239        inlet.rate = -flow_rate_control
     240        outlet.rate = flow_rate_control               
     241           
     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
    206246            log_to_file(log_filename, s)
    207247
    … …  
    210250
    211251
    212     else: # inlet_depth < 0.01:
     252    else: #inlet.depth < 0.01:
    213253        Q = barrel_velocity = outlet_culvert_depth = 0.0
    214254       
Note: See TracChangeset for help on using the changeset viewer.