Changeset 7073


Ignore:
Timestamp:
May 25, 2009, 2:58:05 PM (11 years ago)
Author:
ole
Message:

Updated Boyd culvert flow routine with new tests from Rudy van Drie.

Location:
anuga_core/source/anuga/culvert_flows
Files:
4 added
1 deleted
1 edited

Legend:

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

    r6948 r7073  
    11"""Collection of culvert routines for use with Culvert_flow in culvert_class
     2
     3This module holds various riutines to determine FLOW through CULVERTS and SIMPLE BRIDGES
     4
     5
     6
    27
    38Usage:
     
    3338                                   log_filename=None):
    3439
    35     """Boyd's generalisation of the US department of transportation culvert
    36     model
    37      
    38     The quantity of flow passing through a culvert is controlled by many factors
    39     It could be that the culvert is controlled by the inlet only, with it
    40     being unsubmerged this is effectively equivalent to the weir Equation
    41     Else the culvert could be controlled by the inlet, with it being
    42     submerged, this is effectively the Orifice Equation
    43     Else it may be controlled by down stream conditions where depending on
    44     the down stream depth, the momentum in the culvert etc. flow is controlled
     40    """Boyd's generalisation of the US department of transportation culvert methods
     41       
     42        WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
     43        FULL PIPE OR CRITICAL DEPTH ONLY
     44        For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
     45        The obtain the CORRECT velocity requires an iteration of Depth to Establish
     46        the Normal Depth of flow in the pipe.
     47       
     48        It is proposed to provide this in a seperate routine called
     49    boyd_generalised_culvert_model_complex
     50       
     51    The Boyd Method is based on methods described by the following:
     52        1.
     53        US Dept. Transportation Federal Highway Administration (1965)
     54        Hydraulic Chart for Selection of Highway Culverts.
     55        Hydraulic Engineering Circular No. 5 US Government Printing
     56        2.
     57        US Dept. Transportation Federal Highway Administration (1972)
     58        Capacity charts for the Hydraulic design of highway culverts.
     59        Hydraulic Engineering Circular No. 10 US Government Printing
     60    These documents provide around 60 charts for various configurations of culverts and inlets.
     61       
     62        Note these documents have been superceded by:
     63        2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
     64        Which combines culvert design information previously contained in Hydraulic Engineering Circulars
     65        (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
     66        HEC-5 provides 20 Charts
     67        HEC-10 Provides an additional 36 Charts
     68        HEC-13 Discusses the Design of improved more efficient inlets
     69        HDS-5 Provides 60 sets of Charts
     70
     71        In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
     72        1987 published "Generalised Head Discharge Equations for Culverts".
     73        These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
     74       
     75        It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
     76        The additional charts cover a range of culvert shapes and inlet configurations
     77       
     78   
    4579    """
    4680
     
    4882    from anuga.config import velocity_protection
    4983    from anuga.utilities.numerical_tools import safe_acos as acos
    50 
    51 
     84    print "STARTING..BOYD................."
     85    local_debug ='false'
    5286    if inlet_depth > 0.1: #this value was 0.01:
     87        if local_debug =='true':
     88            print 'Specific E & Deltat Tot E = ',inlet_specific_energy,delta_total_energy
     89            print 'culvert typ = ',culvert_type
    5390        # Water has risen above inlet
    5491       
     
    6198                     
    6299
    63         if culvert_type == 'circle':
    64             # Round culvert (use width as diameter)
    65             diameter = culvert_width           
    66 
     100        if culvert_type =='circle':
     101            # Round culvert (use height as diameter)
     102            diameter = culvert_height   
     103
     104            """
     105            For a circular pipe the Boyd method reviews 3 conditions
     106            1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
     107            2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
     108            3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
     109           
     110            For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
     111            """
     112           
    67113            # Calculate flows for inlet control           
    68114            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
    69115            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
     116            # Note for to SUBMERGED TO OCCUR inlet_specific_energy should be > 1.2 x diameter.... Should Check !!!
    70117
    71118            if log_filename is not None:                                       
    72119                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
    73120                log_to_file(log_filename, s)
    74 
    75 
    76             # FIXME(Ole): Are these functions really for inlet control?                   
    77             if Q_inlet_unsubmerged < Q_inlet_submerged:
    78                 Q = Q_inlet_unsubmerged
    79                 alpha = acos(1 - inlet_depth/diameter)
    80                 flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    81                 outlet_culvert_depth = inlet_depth
    82                 case = 'Inlet unsubmerged'
    83             else:   
    84                 Q = Q_inlet_submerged
    85                 flow_area = (diameter/2)**2 * pi
    86                 outlet_culvert_depth = diameter
    87                 case = 'Inlet submerged'                   
    88 
    89 
    90 
    91             if delta_total_energy < inlet_specific_energy:
     121            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
     122
     123            # THE LOWEST Value will Control Calcs From here
     124            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
     125            dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
     126            dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
     127            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
     128            if dcrit1/diameter  > 0.85:
     129                outlet_culvert_depth = dcrit2
     130            else:
     131                outlet_culvert_depth = dcrit1
     132            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
     133            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
     134            if outlet_culvert_depth >= diameter:
     135                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     136                flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
     137                perimeter = diameter * pi
     138                flow_width= diameter
     139                case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
     140                if local_debug =='true':
     141                    print 'Inlet CTRL Outlet submerged Circular PIPE FULL'
     142            else:
     143                #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
     144                alpha = acos(1-2*outlet_culvert_depth/diameter)*2
     145                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
     146                flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     147                flow_width= diameter*sin(alpha/2.0)
     148                perimeter = alpha*diameter/2.0
     149                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     150                if local_debug =='true':
     151                    print 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     152                    print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,'  ',alpha
     153            if delta_total_energy < inlet_specific_energy:  #  OUTLET CONTROL !!!!
    92154                # Calculate flows for outlet control
    93155               
    94156                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    95                 if outlet_depth > diameter:        # The Outlet is Submerged
     157                if outlet_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
    96158                    outlet_culvert_depth=diameter
    97159                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    98160                    perimeter = diameter * pi
     161                    flow_width= diameter
    99162                    case = 'Outlet submerged'
    100                 elif outlet_depth==0.0:
    101                     outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth
    102                     alpha = acos(1 - inlet_depth/diameter)
    103                     flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    104                     perimeter = alpha*diameter
    105 
    106                     case = 'Outlet depth is zero'                       
    107                 else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    108                     outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
    109                     alpha = acos(1 - outlet_depth/diameter)
    110                     flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
    111                     perimeter = alpha*diameter
    112                     case = 'Outlet is open channel flow'
    113 
    114                 hyd_rad = flow_area/perimeter
    115                
    116                 if log_filename is not None:
    117                     s = 'hydraulic radius at outlet = %f' %hyd_rad
    118                     log_to_file(log_filename, s)
    119 
    120                 # Outlet control velocity using tail water
    121                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    122                 Q_outlet_tailwater = flow_area * culvert_velocity
    123                
    124                 if log_filename is not None:               
    125                     s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    126                     log_to_file(log_filename, s)
    127                    
    128                 Q = min(Q, Q_outlet_tailwater)
    129             else:
    130                 pass
     163                    if local_debug =='true':
     164                        print 'Outlet submerged'
     165                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     166                    # IF  outlet_depth < diameter
     167                    dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
     168                    dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
     169                    if dcrit1/diameter >0.85:
     170                        outlet_culvert_depth= dcrit2
     171                    else:
     172                        outlet_culvert_depth = dcrit1
     173                    if outlet_culvert_depth > diameter:
     174                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     175                        flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
     176                        perimeter = diameter * pi
     177                        flow_width= diameter
     178                        case = 'Outlet unsubmerged PIPE FULL'
     179                        if local_debug =='true':
     180                            print  'Outlet unsubmerged PIPE FULL'
     181                    else:
     182                        alpha = acos(1-2*outlet_culvert_depth/diameter)*2
     183                        flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     184                        flow_width= diameter*sin(alpha/2.0)
     185                        perimeter = alpha*diameter/2.0
     186                        case = 'Outlet is open channel flow we will for now assume critical depth'
     187                        if local_debug =='true':
     188                            print 'Q Outlet Depth and ALPHA = ',Q,' ',outlet_culvert_depth,'  ',alpha
     189                            print 'Outlet is open channel flow we will for now assume critical depth'
     190            if local_debug =='true':
     191                print 'FLOW AREA = ',flow_area
     192                print 'PERIMETER = ',perimeter
     193                print 'Q Interim = ',Q
     194            hyd_rad = flow_area/perimeter
     195
     196            if log_filename is not None:
     197                s = 'hydraulic radius at outlet = %f' %hyd_rad
     198                log_to_file(log_filename, s)
     199
     200            # Outlet control velocity using tail water
     201            if local_debug =='true':
     202                print 'GOT IT ALL CALCULATING Velocity'
     203                print 'HydRad = ',hyd_rad     
     204            culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
     205            Q_outlet_tailwater = flow_area * culvert_velocity
     206            if local_debug =='true':
     207                print 'VELOCITY = ',culvert_velocity
     208                print 'Outlet Ctrl Q = ',Q_outlet_tailwater
     209            if log_filename is not None:               
     210                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     211                log_to_file(log_filename, s)
     212            Q = min(Q, Q_outlet_tailwater)
     213            if local_debug =='true':
     214                print ('%s,%.3f,%.3f' %('dcrit 1 , dcit2 =',dcrit1,dcrit2))
     215                print ('%s,%.3f,%.3f,%.3f' %('Q and Velocity and Depth=',Q,culvert_velocity,outlet_culvert_depth))
     216
     217        else:
     218            pass
    131219                #FIXME(Ole): What about inlet control?
    132 
    133 
    134         else:
    135             # Box culvert (rectangle or square)
     220        # ====  END OF CODE BLOCK FOR "IF" CIRCULAR PIPE
     221
     222        #else....
     223        if culvert_type =='box':
     224            if local_debug =='true':
     225                print 'BOX CULVERT'
     226            # Box culvert (rectangle or square)   ========================================================================================================================
    136227
    137228            # Calculate flows for inlet control
    138229            height = culvert_height
    139             width = culvert_width           
     230            width = culvert_width
     231            flow_width=culvert_width           
    140232           
    141233            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     
    146238                log_to_file(log_filename, s)
    147239
    148 
    149240            # FIXME(Ole): Are these functions really for inlet control?   
    150241            if Q_inlet_unsubmerged < Q_inlet_submerged:
    151242                Q = Q_inlet_unsubmerged
    152                 flow_area = width*inlet_depth
    153                 outlet_culvert_depth = inlet_depth
    154                 case = 'Inlet unsubmerged'
     243                dcrit = (Q**2/g/width**2)**0.333333
     244                if dcrit > height:
     245                    dcrit = height
     246                flow_area = width*dcrit
     247                outlet_culvert_depth = dcrit
     248                case = 'Inlet unsubmerged Box Acts as Weir'
    155249            else:   
    156250                Q = Q_inlet_submerged
    157251                flow_area = width*height
    158252                outlet_culvert_depth = height
    159                 case = 'Inlet submerged'                   
    160 
     253                case = 'Inlet submerged Box Acts as Orifice'                   
     254
     255            dcrit = (Q**2/g/width**2)**0.333333
     256
     257            outlet_culvert_depth = dcrit
     258            if outlet_culvert_depth > height:
     259                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
     260                flow_area = width*height  # Cross sectional area of flow in the culvert
     261                perimeter = 2*(width+height)
     262                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
     263            else:
     264                flow_area = width * outlet_culvert_depth
     265                perimeter = width+2*outlet_culvert_depth
     266                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     267       
    161268            if delta_total_energy < inlet_specific_energy:
    162269                # Calculate flows for outlet control
     
    168275                    perimeter=2.0*(width+height)
    169276                    case = 'Outlet submerged'
    170                 elif outlet_depth==0.0:
    171                     outlet_culvert_depth=inlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
    172                     flow_area=width*inlet_depth
    173                     perimeter=(width+2.0*inlet_depth)
    174                     case = 'Outlet depth is zero'                       
    175277                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    176                     outlet_culvert_depth=outlet_depth
    177                     flow_area=width*outlet_depth
    178                     perimeter=(width+2.0*outlet_depth)
    179                     case = 'Outlet is open channel flow'
     278                    dcrit = (Q**2/g/width**2)**0.333333
     279                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
     280                    if outlet_culvert_depth > height:
     281                        outlet_culvert_depth=height
     282                        flow_area=width*height
     283                        perimeter=2.0*(width+height)
     284                        case = 'Outlet is Flowing Full'
     285                    else:
     286                        flow_area=width*outlet_culvert_depth
     287                        perimeter=(width+2.0*outlet_culvert_depth)
     288                        case = 'Outlet is open channel flow'
    180289
    181290                hyd_rad = flow_area/perimeter
    182                 
     291           
    183292                if log_filename is not None:                               
    184293                    s = 'hydraulic radius at outlet = %f' % hyd_rad
     
    196305                pass
    197306                #FIXME(Ole): What about inlet control?
     307        # ====  END OF CODE BLOCK FOR "IF" BOX
    198308
    199309               
    200         # Common code for circle and square geometries
     310        # Common code for circle and box geometries ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    201311        if log_filename is not None:
    202312            log_to_file(log_filename, 'Case: "%s"' % case)
     
    205315            s = 'Flow Rate Control = %f' % Q
    206316            log_to_file(log_filename, s)
    207 
    208            
    209         culv_froude=sqrt(Q**2*width/(g*flow_area**3))
    210        
     317        culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
     318        if local_debug =='true':
     319            print 'FLOW AREA = ',flow_area
     320            print 'PERIMETER = ',perimeter
     321            print 'Q final = ',Q
     322            print 'FROUDE = ',culv_froude
    211323        if log_filename is not None:                           
    212324            s = 'Froude in Culvert = %f' % culv_froude
     
    215327        # Determine momentum at the outlet
    216328        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
    217 
     329    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow 
    218330
    219331    else: # inlet_depth < 0.01:
Note: See TracChangeset for help on using the changeset viewer.