Ignore:
Timestamp:
Aug 30, 2010, 5:53:45 PM (12 years ago)
Author:
steve
Message:

Added new tests

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/structures/boyd_box_routine.py

    r7977 r7980  
    88
    99
     10import culvert_routine
     11from anuga.config import velocity_protection
     12from anuga.utilities.numerical_tools import safe_acos as acos
    1013
    11 def boyd_box(culvert):
     14from math import pi, sqrt, sin, cos
     15from anuga.config import g
     16
     17
     18class Boyd_box_routine(culvert_routine.Culvert_routine):
    1219    """Boyd's generalisation of the US department of transportation culvert methods
    1320
    14         WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
    15         FULL PIPE OR CRITICAL DEPTH ONLY
    16         For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
    17         The obtain the CORRECT velocity requires an iteration of Depth to Establish
    18         the Normal Depth of flow in the pipe.
     21        WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
     22        FULL PIPE OR CRITICAL DEPTH ONLY
     23        For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
     24        The obtain the CORRECT velocity requires an iteration of Depth to Establish
     25        the Normal Depth of flow in the pipe.
    1926
    20         It is proposed to provide this in a seperate routine called
    21     boyd_generalised_culvert_model_complex
     27        It is proposed to provide this in a seperate routine called
     28        boyd_generalised_culvert_model_complex
    2229
    23     The Boyd Method is based on methods described by the following:
    24         1.
    25         US Dept. Transportation Federal Highway Administration (1965)
    26         Hydraulic Chart for Selection of Highway Culverts.
    27         Hydraulic Engineering Circular No. 5 US Government Printing
    28         2.
    29         US Dept. Transportation Federal Highway Administration (1972)
    30         Capacity charts for the Hydraulic design of highway culverts.
    31         Hydraulic Engineering Circular No. 10 US Government Printing
    32     These documents provide around 60 charts for various configurations of culverts and inlets.
     30        The Boyd Method is based on methods described by the following:
     31        1.
     32        US Dept. Transportation Federal Highway Administration (1965)
     33        Hydraulic Chart for Selection of Highway Culverts.
     34        Hydraulic Engineering Circular No. 5 US Government Printing
     35        2.
     36        US Dept. Transportation Federal Highway Administration (1972)
     37        Capacity charts for the Hydraulic design of highway culverts.
     38        Hydraulic Engineering Circular No. 10 US Government Printing
     39        These documents provide around 60 charts for various configurations of culverts and inlets.
    3340
    34         Note these documents have been superceded by:
    35         2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
    36         Which combines culvert design information previously contained in Hydraulic Engineering Circulars
    37         (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
    38         HEC-5 provides 20 Charts
    39         HEC-10 Provides an additional 36 Charts
    40         HEC-13 Discusses the Design of improved more efficient inlets
    41         HDS-5 Provides 60 sets of Charts
     41        Note these documents have been superceded by:
     42        2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
     43        Which combines culvert design information previously contained in Hydraulic Engineering Circulars
     44        (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
     45        HEC-5 provides 20 Charts
     46        HEC-10 Provides an additional 36 Charts
     47        HEC-13 Discusses the Design of improved more efficient inlets
     48        HDS-5 Provides 60 sets of Charts
    4249
    43         In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
    44         1987 published "Generalised Head Discharge Equations for Culverts".
    45         These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
     50        In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
     51        1987 published "Generalised Head Discharge Equations for Culverts".
     52        These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
    4653
    47         It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
    48         The additional charts cover a range of culvert shapes and inlet configurations
    49 
    50     """
    51 
    52     # Calculate flows for inflow control
    53 
    54     Q_inflow_unsubmerged = 0.540*g**0.5*width*inflow_specific_energy**1.50 # Flow based on inflow Ctrl inflow Unsubmerged
    55     Q_inflow_submerged = 0.702*g**0.5*width*height**0.89*inflow_specific_energy**0.61  # Flow based on inflow Ctrl inflow Submerged
    56 
    57     if log_filename is not None:
    58         s = 'Q_inflow_unsubmerged = %.6f, Q_inflow_submerged = %.6f' %(Q_inflow_unsubmerged, Q_inflow_submerged)
    59         log_to_file(log_filename, s)
    60 
    61     # FIXME(Ole): Are these functions really for inflow control?
    62     if Q_inflow_unsubmerged < Q_inflow_submerged:
    63         Q = Q_inflow_unsubmerged
    64         dcrit = (Q**2/g/width**2)**0.333333
    65         if dcrit > height:
    66             dcrit = height
    67         flow_area = width*dcrit
    68         outflow_culvert_depth = dcrit
    69         case = 'inflow unsubmerged Box Acts as Weir'
    70     else:
    71         Q = Q_inflow_submerged
    72         flow_area = width*height
    73         outflow_culvert_depth = height
    74         case = 'inflow submerged Box Acts as Orifice'
    75 
    76     dcrit = (Q**2/g/width**2)**0.333333
    77 
    78     outflow_culvert_depth = dcrit
    79     if outflow_culvert_depth > height:
    80         outflow_culvert_depth = height  # Once again the pipe is flowing full not partfull
    81         flow_area = width*height  # Cross sectional area of flow in the culvert
    82         perimeter = 2*(width+height)
    83         case = 'inflow CTRL outflow unsubmerged PIPE PART FULL'
    84     else:
    85         flow_area = width * outflow_culvert_depth
    86         perimeter = width+2*outflow_culvert_depth
    87         case = 'inflow CTRL Culvert is open channel flow we will for now assume critical depth'
    88 
    89     if delta_total_energy < inflow_specific_energy:
    90         # Calculate flows for outflow control
    91 
    92         # Determine the depth at the outflow relative to the depth of flow in the Culvert
    93         if outflow_depth > height:        # The outflow is Submerged
    94             outflow_culvert_depth=height
    95             flow_area=width*height       # Cross sectional area of flow in the culvert
    96             perimeter=2.0*(width+height)
    97             case = 'outflow submerged'
    98         else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    99             dcrit = (Q**2/g/width**2)**0.333333
    100             outflow_culvert_depth=dcrit   # For purpose of calculation assume the outflow depth = Critical Depth
    101             if outflow_culvert_depth > height:
    102                 outflow_culvert_depth=height
    103                 flow_area=width*height
    104                 perimeter=2.0*(width+height)
    105                 case = 'outflow is Flowing Full'
    106             else:
    107                 flow_area=width*outflow_culvert_depth
    108                 perimeter=(width+2.0*outflow_culvert_depth)
    109                 case = 'outflow is open channel flow'
    110 
    111         hyd_rad = flow_area/perimeter
    112 
    113         if log_filename is not None:
    114             s = 'hydraulic radius at outflow = %f' % hyd_rad
    115             log_to_file(log_filename, s)
    116 
    117         # outflow control velocity using tail water
    118         culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    119         Q_outflow_tailwater = flow_area * culvert_velocity
    120 
    121         if log_filename is not None:
    122             s = 'Q_outflow_tailwater = %.6f' % Q_outflow_tailwater
    123             log_to_file(log_filename, s)
    124         Q = min(Q, Q_outflow_tailwater)
    125 
    126     return Q
     54        It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
     55        The additional charts cover a range of culvert shapes and inlet configurations
    12756
    12857
    129 if __name__ == "__main__":
     58        """
     59
     60    def __init__(self, culvert, manning=0.0):
     61
     62        culvert_routine.Culvert_routine.__init__(self, culvert, manning)
    13063
    13164
    132     g=9.81
    133     culvert_slope=0.1  # Downward
    13465
    135     inlet_depth=2.0
    136     outlet_depth=0.0
     66    def __call__(self):
    13767
    138     inlet_velocity=0.0,
    139     outlet_velocity=0.0,
     68        self.determine_inflow()
    14069
    141     culvert_length=4.0
    142     culvert_width=1.2
    143     culvert_height=0.75
     70        local_debug ='false'
     71       
     72        if self.inflow.get_average_height() > 0.01: #this value was 0.01:
     73            if local_debug =='true':
     74                log.critical('Specific E & Deltat Tot E = %s, %s'
     75                             % (str(self.inflow.get_average_specific_energy()),
     76                                str(self.delta_total_energy)))
     77                log.critical('culvert type = %s' % str(culvert_type))
     78            # Water has risen above inlet
    14479
    145     culvert_type='box'
    146     manning=0.013
    147     sum_loss=0.0
     80            if self.log_filename is not None:
     81                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
     82                log_to_file(self.log_filename, s)
    14883
    149     inlet_specific_energy=inlet_depth #+0.5*v**2/g
    150     z_in = 0.0
    151     z_out = -culvert_length*culvert_slope/100
    152     E_in = z_in+inlet_depth # +
    153     E_out = z_out+outlet_depth # +
    154     delta_total_energy = E_in-E_out
     84            msg = 'Specific energy at inlet is negative'
     85            assert self.inflow.get_average_specific_energy() >= 0.0, msg
    15586
    156     Q = boyd_box(culvert_height, culvert_width, culvert_width, inlet_specific_energy)
     87            height = self.culvert_height
     88            width = self.culvert_width
     89            flow_width = self.culvert_width
    15790
    158     print 'Q ',Q
     91            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     92            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_average_specific_energy()**0.61  # Flow based on Inlet Ctrl Inlet Submerged
     93
     94            # FIXME(Ole): Are these functions really for inlet control?
     95            if Q_inlet_unsubmerged < Q_inlet_submerged:
     96                Q = Q_inlet_unsubmerged
     97                dcrit = (Q**2/g/width**2)**0.333333
     98                if dcrit > height:
     99                    dcrit = height
     100                flow_area = width*dcrit
     101                outlet_culvert_depth = dcrit
     102                case = 'Inlet unsubmerged Box Acts as Weir'
     103            else:
     104                Q = Q_inlet_submerged
     105                flow_area = width*height
     106                outlet_culvert_depth = height
     107                case = 'Inlet submerged Box Acts as Orifice'
     108
     109            dcrit = (Q**2/g/width**2)**0.333333
     110
     111            outlet_culvert_depth = dcrit
     112            if outlet_culvert_depth > height:
     113                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
     114                flow_area = width*height  # Cross sectional area of flow in the culvert
     115                perimeter = 2*(width+height)
     116                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
     117            else:
     118                flow_area = width * outlet_culvert_depth
     119                perimeter = width+2*outlet_culvert_depth
     120                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     121
     122            if self.delta_total_energy < self.inflow.get_average_specific_energy():
     123                # Calculate flows for outlet control
     124
     125                # Determine the depth at the outlet relative to the depth of flow in the Culvert
     126                if self.outflow.get_average_height() > height:        # The Outlet is Submerged
     127                    outlet_culvert_depth=height
     128                    flow_area=width*height       # Cross sectional area of flow in the culvert
     129                    perimeter=2.0*(width+height)
     130                    case = 'Outlet submerged'
     131                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     132                    dcrit = (Q**2/g/width**2)**0.333333
     133                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
     134                    if outlet_culvert_depth > height:
     135                        outlet_culvert_depth=height
     136                        flow_area=width*height
     137                        perimeter=2.0*(width+height)
     138                        case = 'Outlet is Flowing Full'
     139                    else:
     140                        flow_area=width*outlet_culvert_depth
     141                        perimeter=(width+2.0*outlet_culvert_depth)
     142                        case = 'Outlet is open channel flow'
     143
     144                hyd_rad = flow_area/perimeter
     145
     146                if self.log_filename is not None:
     147                    s = 'hydraulic radius at outlet = %f' % hyd_rad
     148                    log_to_file(self.log_filename, s)
     149
     150                # Outlet control velocity using tail water
     151                culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
     152                Q_outlet_tailwater = flow_area * culvert_velocity
     153
     154                if self.log_filename is not None:
     155                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
     156                    log_to_file(self.log_filename, s)
     157                Q = min(Q, Q_outlet_tailwater)
     158            else:
     159                pass
     160                #FIXME(Ole): What about inlet control?
     161
     162            culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
     163            if local_debug =='true':
     164                log.critical('FLOW AREA = %s' % str(flow_area))
     165                log.critical('PERIMETER = %s' % str(perimeter))
     166                log.critical('Q final = %s' % str(Q))
     167                log.critical('FROUDE = %s' % str(culv_froude))
     168
     169            # Determine momentum at the outlet
     170            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
     171
     172        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
     173
     174        else: # self.inflow.get_average_height() < 0.01:
     175            Q = barrel_velocity = outlet_culvert_depth = 0.0
     176
     177        # Temporary flow limit
     178        if barrel_velocity > self.max_velocity:
     179            barrel_velocity = self.max_velocity
     180            Q = flow_area * barrel_velocity
     181
     182
     183
     184       
     185
     186        return Q, barrel_velocity, outlet_culvert_depth
     187
     188
     189
Note: See TracChangeset for help on using the changeset viewer.