Changeset 7978


Ignore:
Timestamp:
Aug 27, 2010, 2:30:30 PM (14 years ago)
Author:
habili
Message:

Refactoring of the culvert code.

Location:
trunk/anuga_core/source/anuga/structures
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/structures/culvert_operator.py

    r7977 r7978  
    33import anuga.utilities.log as log
    44import box_culvert
     5import culvert_routines
    56
    67class Culvert_operator:
    … …  
    1718    def __init__(self,
    1819                 domain,
    19                  end_point0=None,
    20                  end_point1=None,
    21                  width=None,
     20                 end_point0,
     21                 end_point1,
     22                 width,
    2223                 height=None,
    2324                 verbose=False):
    … …  
    4344        timestep = self.domain.get_timestep()
    4445
    45         inflow  = self.inlets[0]
    46         outflow = self.inlets[1]
    47 
    48         # Determine flow direction based on total energy difference
    49         delta_total_energy = inflow.get_average_total_energy() - outflow.get_average_total_energy()
    50 
    51         if delta_total_energy < 0:
    52             inflow  = self.inlets[1]
    53             outflow = self.inlets[0]
    54             delta_total_energy = -delta_total_energy
    55 
    56         delta_z = inflow.get_average_elevation() - outflow.get_average_elevation()
    57         culvert_slope = delta_z/self.culvert.get_culvert_length()
    58 
    59         # Determine controlling energy (driving head) for culvert
    60         if inflow.get_average_specific_energy() > delta_total_energy:
    61             # Outlet control
    62             driving_head = delta_total_energy
    63         else:
    64             # Inlet control
    65             driving_head = inflow.get_average_specific_energy()
    66            
    67         # Transfer
    68         from culvert_routines import boyd_box, boyd_circle
    69         Q, barrel_velocity, culvert_outlet_depth =\
    70                               boyd_circle(inflow.get_average_height(),
    71                                          outflow.get_average_height(),
    72                                          inflow.get_average_speed(),
    73                                          outflow.get_average_speed(),
    74                                          inflow.get_average_specific_energy(),
    75                                          delta_total_energy,
    76                                          culvert_length=self.culvert.get_culvert_length(),
    77                                          culvert_width=self.width,
    78                                          culvert_height=self.height,
    79                                          manning=0.01)
     46        from culvert_routines import Culvert_routines
     47        culvert_routine = culvert_routines.Culvert_routines(self.culvert)
     48       
     49        Q, barrel_velocity, culvert_outlet_depth = culvert_routine.boyd_circle()
    8050
    8151        transfer_water = Q*timestep
     52
     53        inflow = culvert_routine.get_inflow()
     54        outflow = culvert_routine.get_outflow()
    8255
    8356        inflow.set_heights(inflow.get_average_height() - transfer_water)

  • trunk/anuga_core/source/anuga/structures/culvert_routines.py

    r7977 r7978  
    1 """Collection of culvert routines for use with Culvert_operator
    2 
    3 This module holds various routines to determine FLOW through CULVERTS and SIMPLE BRIDGES
    4 
    5 
    6 
    7 
    8 Usage:
    9 
    10    
    11 
    12 """
    13 
    14 #NOTE:
    15 # Inlet control:  Delta_total_energy > inlet_specific_energy
    16 # Outlet control: Delta_total_energy < inlet_specific_energy
    17 # where total energy is (w + 0.5*v^2/g) and
    18 # specific energy is (h + 0.5*v^2/g)
    19 
    201from anuga.config import velocity_protection
    212from anuga.utilities.numerical_tools import safe_acos as acos
    … …  
    245from anuga.config import g
    256
    26 def boyd_box(inlet_depth,
    27              outlet_depth,
    28              inlet_velocity,
    29              outlet_velocity,
    30              inlet_specific_energy,
    31              delta_total_energy,
    32              culvert_length=0.0,
    33              culvert_width=0.0,
    34              culvert_height=0.0,
    35              manning=0.0,
    36              sum_loss=0.0,
    37              max_velocity=10.0,
    38              log_filename=None):
    39 
    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    
     7class Culvert_routines:
     8    """Collection of culvert routines for use with Culvert_operator
     9
     10    This module holds various routines to determine FLOW through CULVERTS and SIMPLE BRIDGES
     11
     12    Usage:
     13
     14    NOTE:
     15     Inlet control:  self.delta_total_energy > self.inflow.get_average_specific_energy()
     16     Outlet control: self.delta_total_energy < self.inflow.get_average_specific_energy()
     17     where total energy is (w + 0.5*v^2/g) and
     18     specific energy is (h + 0.5*v^2/g)
    7919    """
    8020
    81     local_debug ='false'
    82     if inlet_depth > 0.1: #this value was 0.01:
    83         if local_debug =='true':
    84             log.critical('Specific E & Deltat Tot E = %s, %s'
    85                          % (str(inlet_specific_energy),
    86                             str(delta_total_energy)))
    87             log.critical('culvert type = %s' % str(culvert_type))
    88         # Water has risen above inlet
    89        
    90         if log_filename is not None:                           
    91             s = 'Specific energy  = %f m' % inlet_specific_energy
    92             log_to_file(log_filename, s)
    93        
    94         msg = 'Specific energy at inlet is negative'
    95         assert inlet_specific_energy >= 0.0, msg
    96                      
    97         height = culvert_height
    98         width = culvert_width
    99         flow_width = culvert_width           
    100        
    101         Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    102         Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    103 
    104         # FIXME(Ole): Are these functions really for inlet control?   
    105         if Q_inlet_unsubmerged < Q_inlet_submerged:
    106             Q = Q_inlet_unsubmerged
     21    def __init__(self, culvert, manning=0.0):
     22       
     23        self.inlets = culvert.inlets
     24       
     25        self.inflow  = self.inlets[0]
     26        self.outflow = self.inlets[1]
     27       
     28        self.culvert_length = culvert.get_culvert_length()
     29        self.culvert_width = culvert.get_culvert_width()
     30        self.culvert_height = culvert.get_culvert_height()
     31        self.sum_loss = 0.0
     32        self.max_velocity = 10.0
     33        self.manning = manning
     34        self.log_filename = None
     35       
     36        # Determine flow direction based on total energy difference
     37        self.delta_total_energy = self.inflow.get_average_total_energy() - self.outflow.get_average_total_energy()
     38
     39        if self.delta_total_energy < 0:
     40            self.inflow  = self.inlets[1]
     41            self.outflow = self.inlets[0]
     42            self.delta_total_energy = -self.delta_total_energy
     43
     44        #delta_z = self.self.inflow.get_average_elevation() - self.self.outflow.get_average_elevation()
     45        #culvert_slope = delta_z/self.culvert.get_self.culvert_length()
     46
     47        # Determine controlling energy (driving head) for culvert
     48        #if self.self.inflow.get_average_specific_energy() > self.self.delta_total_energy:
     49        #    # Outlet control
     50        #    driving_head = self.delta_total_energy
     51        #else:
     52            # Inlet control
     53        #    driving_head = self.inflow.get_average_specific_energy()
     54           
     55       
     56    def get_inflow(self):
     57       
     58        return self.inflow
     59       
     60       
     61    def get_outflow(self):
     62   
     63        return self.outflow
     64   
     65       
     66    def boyd_box(self):
     67
     68        """Boyd's generalisation of the US department of transportation culvert methods
     69       
     70        WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER
     71        FULL PIPE OR CRITICAL DEPTH ONLY
     72        For Supercritical flow this is UNDERESTIMATING the Outlet Velocity
     73        The obtain the CORRECT velocity requires an iteration of Depth to Establish
     74        the Normal Depth of flow in the pipe.
     75       
     76        It is proposed to provide this in a seperate routine called
     77        boyd_generalised_culvert_model_complex
     78       
     79        The Boyd Method is based on methods described by the following:
     80        1.
     81        US Dept. Transportation Federal Highway Administration (1965)
     82        Hydraulic Chart for Selection of Highway Culverts.
     83        Hydraulic Engineering Circular No. 5 US Government Printing
     84        2.
     85        US Dept. Transportation Federal Highway Administration (1972)
     86        Capacity charts for the Hydraulic design of highway culverts.
     87        Hydraulic Engineering Circular No. 10 US Government Printing
     88        These documents provide around 60 charts for various configurations of culverts and inlets.
     89       
     90        Note these documents have been superceded by:
     91        2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),
     92        Which combines culvert design information previously contained in Hydraulic Engineering Circulars
     93        (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.
     94        HEC-5 provides 20 Charts
     95        HEC-10 Provides an additional 36 Charts
     96        HEC-13 Discusses the Design of improved more efficient inlets
     97        HDS-5 Provides 60 sets of Charts
     98
     99        In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in
     100        1987 published "Generalised Head Discharge Equations for Culverts".
     101        These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.
     102       
     103        It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.
     104        The additional charts cover a range of culvert shapes and inlet configurations
     105       
     106       
     107        """
     108
     109        local_debug ='false'
     110        if self.inflow.get_average_height() > 0.1: #this value was 0.01:
     111            if local_debug =='true':
     112                log.critical('Specific E & Deltat Tot E = %s, %s'
     113                             % (str(self.inflow.get_average_specific_energy()),
     114                                str(self.delta_total_energy)))
     115                log.critical('culvert type = %s' % str(culvert_type))
     116            # Water has risen above inlet
     117           
     118            if self.log_filename is not None:                           
     119                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
     120                log_to_file(self.log_filename, s)
     121           
     122            msg = 'Specific energy at inlet is negative'
     123            assert self.inflow.get_average_specific_energy() >= 0.0, msg
     124                         
     125            height = self.culvert_height
     126            width = self.culvert_width
     127            flow_width = self.culvert_width           
     128           
     129            Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     130            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
     131
     132            # FIXME(Ole): Are these functions really for inlet control?   
     133            if Q_inlet_unsubmerged < Q_inlet_submerged:
     134                Q = Q_inlet_unsubmerged
     135                dcrit = (Q**2/g/width**2)**0.333333
     136                if dcrit > height:
     137                    dcrit = height
     138                flow_area = width*dcrit
     139                outlet_culvert_depth = dcrit
     140                case = 'Inlet unsubmerged Box Acts as Weir'
     141            else:   
     142                Q = Q_inlet_submerged
     143                flow_area = width*height
     144                outlet_culvert_depth = height
     145                case = 'Inlet submerged Box Acts as Orifice'                   
     146
    107147            dcrit = (Q**2/g/width**2)**0.333333
    108             if dcrit > height:
    109                 dcrit = height
    110             flow_area = width*dcrit
     148
    111149            outlet_culvert_depth = dcrit
    112             case = 'Inlet unsubmerged Box Acts as Weir'
    113         else:   
    114             Q = Q_inlet_submerged
    115             flow_area = width*height
    116             outlet_culvert_depth = height
    117             case = 'Inlet submerged Box Acts as Orifice'                   
    118 
    119         dcrit = (Q**2/g/width**2)**0.333333
    120 
    121         outlet_culvert_depth = dcrit
    122         if outlet_culvert_depth > height:
    123             outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
    124             flow_area = width*height  # Cross sectional area of flow in the culvert
    125             perimeter = 2*(width+height)
    126             case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
    127         else:
    128             flow_area = width * outlet_culvert_depth
    129             perimeter = width+2*outlet_culvert_depth
    130             case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    131 
    132         if delta_total_energy < inlet_specific_energy:
    133             # Calculate flows for outlet control
    134            
    135             # Determine the depth at the outlet relative to the depth of flow in the Culvert
    136             if outlet_depth > height:        # The Outlet is Submerged
    137                 outlet_culvert_depth=height
    138                 flow_area=width*height       # Cross sectional area of flow in the culvert
    139                 perimeter=2.0*(width+height)
    140                 case = 'Outlet submerged'
    141             else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    142                 dcrit = (Q**2/g/width**2)**0.333333
    143                 outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
    144                 if outlet_culvert_depth > height:
     150            if outlet_culvert_depth > height:
     151                outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
     152                flow_area = width*height  # Cross sectional area of flow in the culvert
     153                perimeter = 2*(width+height)
     154                case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
     155            else:
     156                flow_area = width * outlet_culvert_depth
     157                perimeter = width+2*outlet_culvert_depth
     158                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     159
     160            if self.delta_total_energy < self.inflow.get_average_specific_energy():
     161                # Calculate flows for outlet control
     162               
     163                # Determine the depth at the outlet relative to the depth of flow in the Culvert
     164                if self.outflow.get_average_height() > height:        # The Outlet is Submerged
    145165                    outlet_culvert_depth=height
    146                     flow_area=width*height
     166                    flow_area=width*height       # Cross sectional area of flow in the culvert
    147167                    perimeter=2.0*(width+height)
    148                     case = 'Outlet is Flowing Full'
    149                 else:
    150                     flow_area=width*outlet_culvert_depth
    151                     perimeter=(width+2.0*outlet_culvert_depth)
    152                     case = 'Outlet is open channel flow'
    153 
    154             hyd_rad = flow_area/perimeter
    155 
    156             if log_filename is not None:                               
    157                 s = 'hydraulic radius at outlet = %f' % hyd_rad
    158                 log_to_file(log_filename, s)
    159 
    160             # Outlet control velocity using tail water
    161             culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    162             Q_outlet_tailwater = flow_area * culvert_velocity
    163 
    164             if log_filename is not None:                           
    165                 s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
    166                 log_to_file(log_filename, s)
    167             Q = min(Q, Q_outlet_tailwater)
    168         else:
    169             pass
    170             #FIXME(Ole): What about inlet control?
    171 
    172         culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
    173         if local_debug =='true':
    174             log.critical('FLOW AREA = %s' % str(flow_area))
    175             log.critical('PERIMETER = %s' % str(perimeter))
    176             log.critical('Q final = %s' % str(Q))
    177             log.critical('FROUDE = %s' % str(culv_froude))
    178  
    179         # Determine momentum at the outlet
    180         barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
    181 
    182     # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow 
    183 
    184     else: # inlet_depth < 0.01:
    185         Q = barrel_velocity = outlet_culvert_depth = 0.0
    186 
    187     # Temporary flow limit
    188     if barrel_velocity > max_velocity:
    189         barrel_velocity = max_velocity
    190         Q = flow_area * barrel_velocity
    191        
    192     return Q, barrel_velocity, outlet_culvert_depth
    193 
    194 
    195 def boyd_circle(inlet_depth,
    196                outlet_depth,
    197                inlet_velocity,
    198                outlet_velocity,
    199                inlet_specific_energy,
    200                delta_total_energy,
    201                culvert_length=0.0,
    202                culvert_width=0.0,
    203                culvert_height=0.0,
    204                manning=0.0,
    205                sum_loss=0.0,
    206                max_velocity=10.0,
    207                log_filename=None):
    208     """
    209     For a circular pipe the Boyd method reviews 3 conditions
    210     1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
    211     2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
    212     3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
    213    
    214     For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
    215     """
    216    
    217     diameter = culvert_height
    218    
    219     local_debug ='false'
    220     if inlet_depth > 0.1: #this value was 0.01:
    221         if local_debug =='true':
    222             log.critical('Specific E & Deltat Tot E = %s, %s'
    223                          % (str(inlet_specific_energy),
    224                             str(delta_total_energy)))
    225             log.critical('culvert type = %s' % str(culvert_type))
    226         # Water has risen above inlet
    227        
    228         if log_filename is not None:                           
    229             s = 'Specific energy  = %f m' % inlet_specific_energy
    230             log_to_file(log_filename, s)
    231        
    232         msg = 'Specific energy at inlet is negative'
    233         assert inlet_specific_energy >= 0.0, msg
    234                      
    235         # Calculate flows for inlet control           
    236         Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
    237         Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
    238         # Note for to SUBMERGED TO OCCUR inlet_specific_energy should be > 1.2 x diameter.... Should Check !!!
    239 
    240         if log_filename is not None:                                       
    241             s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
    242             log_to_file(log_filename, s)
    243         Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
    244 
    245         # THE LOWEST Value will Control Calcs From here
    246         # Calculate Critical Depth Based on the Adopted Flow as an Estimate
    247         dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
    248         dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
    249         # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
    250         if dcrit1/diameter  > 0.85:
    251             outlet_culvert_depth = dcrit2
    252         else:
    253             outlet_culvert_depth = dcrit1
    254         #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
    255         # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
    256         if outlet_culvert_depth >= diameter:
    257             outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
    258             flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    259             perimeter = diameter * pi
    260             flow_width= diameter
    261             case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
    262             if local_debug == 'true':
    263                 log.critical('Inlet CTRL Outlet submerged Circular '
    264                              'PIPE FULL')
    265         else:
    266             #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
    267             alpha = acos(1-2*outlet_culvert_depth/diameter)*2
    268             #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
    269             flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
    270             flow_width= diameter*sin(alpha/2.0)
    271             perimeter = alpha*diameter/2.0
    272             case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    273             if local_debug =='true':
    274                 log.critical('INLET CTRL Culvert is open channel flow '
    275                              'we will for now assume critical depth')
    276                 log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
    277                              % (str(Q), str(outlet_culvert_depth),
    278                                 str(alpha)))
    279         if delta_total_energy < inlet_specific_energy:  #  OUTLET CONTROL !!!!
    280             # Calculate flows for outlet control
    281            
    282             # Determine the depth at the outlet relative to the depth of flow in the Culvert
    283             if outlet_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
    284                 outlet_culvert_depth=diameter
     168                    case = 'Outlet submerged'
     169                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     170                    dcrit = (Q**2/g/width**2)**0.333333
     171                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
     172                    if outlet_culvert_depth > height:
     173                        outlet_culvert_depth=height
     174                        flow_area=width*height
     175                        perimeter=2.0*(width+height)
     176                        case = 'Outlet is Flowing Full'
     177                    else:
     178                        flow_area=width*outlet_culvert_depth
     179                        perimeter=(width+2.0*outlet_culvert_depth)
     180                        case = 'Outlet is open channel flow'
     181
     182                hyd_rad = flow_area/perimeter
     183
     184                if self.log_filename is not None:                               
     185                    s = 'hydraulic radius at outlet = %f' % hyd_rad
     186                    log_to_file(self.log_filename, s)
     187
     188                # Outlet control velocity using tail water
     189                culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
     190                Q_outlet_tailwater = flow_area * culvert_velocity
     191
     192                if self.log_filename is not None:                           
     193                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
     194                    log_to_file(self.log_filename, s)
     195                Q = min(Q, Q_outlet_tailwater)
     196            else:
     197                pass
     198                #FIXME(Ole): What about inlet control?
     199
     200            culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
     201            if local_debug =='true':
     202                log.critical('FLOW AREA = %s' % str(flow_area))
     203                log.critical('PERIMETER = %s' % str(perimeter))
     204                log.critical('Q final = %s' % str(Q))
     205                log.critical('FROUDE = %s' % str(culv_froude))
     206     
     207            # Determine momentum at the outlet
     208            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
     209
     210        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow 
     211
     212        else: # self.inflow.get_average_height() < 0.01:
     213            Q = barrel_velocity = outlet_culvert_depth = 0.0
     214
     215        # Temporary flow limit
     216        if barrel_velocity > self.max_velocity:
     217            barrel_velocity = self.max_velocity
     218            Q = flow_area * barrel_velocity
     219           
     220        return Q, barrel_velocity, outlet_culvert_depth
     221
     222
     223    def boyd_circle(self):
     224        """
     225        For a circular pipe the Boyd method reviews 3 conditions
     226        1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
     227        2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
     228        3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
     229       
     230        For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
     231        """
     232       
     233        diameter = self.culvert_height
     234       
     235        local_debug ='false'
     236        if self.inflow.get_average_height() > 0.1: #this value was 0.01:
     237            if local_debug =='true':
     238                log.critical('Specific E & Deltat Tot E = %s, %s'
     239                             % (str(self.inflow.get_average_specific_energy()),
     240                                str(self.delta_total_energy)))
     241                log.critical('culvert type = %s' % str(culvert_type))
     242            # Water has risen above inlet
     243           
     244            if self.log_filename is not None:                           
     245                s = 'Specific energy  = %f m' % self.inflow.get_average_specific_energy()
     246                log_to_file(self.log_filename, s)
     247           
     248            msg = 'Specific energy at inlet is negative'
     249            assert self.inflow.get_average_specific_energy() >= 0.0, msg
     250                         
     251            # Calculate flows for inlet control           
     252            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*self.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged
     253            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*self.inflow.get_average_specific_energy()**0.63   # Inlet Ctrl Inlet Submerged
     254            # Note for to SUBMERGED TO OCCUR self.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
     255
     256            if self.log_filename is not None:                                       
     257                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
     258                log_to_file(self.log_filename, s)
     259            Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
     260
     261            # THE LOWEST Value will Control Calcs From here
     262            # Calculate Critical Depth Based on the Adopted Flow as an Estimate
     263            dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
     264            dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
     265            # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
     266            if dcrit1/diameter  > 0.85:
     267                outlet_culvert_depth = dcrit2
     268            else:
     269                outlet_culvert_depth = dcrit1
     270            #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
     271            # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
     272            if outlet_culvert_depth >= diameter:
     273                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
    285274                flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    286275                perimeter = diameter * pi
    287276                flow_width= diameter
    288                 case = 'Outlet submerged'
     277                case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
     278                if local_debug == 'true':
     279                    log.critical('Inlet CTRL Outlet submerged Circular '
     280                                 'PIPE FULL')
     281            else:
     282                #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
     283                alpha = acos(1-2*outlet_culvert_depth/diameter)*2
     284                #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
     285                flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     286                flow_width= diameter*sin(alpha/2.0)
     287                perimeter = alpha*diameter/2.0
     288                case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    289289                if local_debug =='true':
    290                     log.critical('Outlet submerged')
    291             else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    292                 # IF  outlet_depth < diameter
    293                 dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
    294                 dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
    295                 if dcrit1/diameter >0.85:
    296                     outlet_culvert_depth= dcrit2
    297                 else:
    298                     outlet_culvert_depth = dcrit1
    299                 if outlet_culvert_depth > diameter:
    300                     outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     290                    log.critical('INLET CTRL Culvert is open channel flow '
     291                                 'we will for now assume critical depth')
     292                    log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
     293                                 % (str(Q), str(outlet_culvert_depth),
     294                                    str(alpha)))
     295            if self.delta_total_energy < self.inflow.get_average_specific_energy():  #  OUTLET CONTROL !!!!
     296                # Calculate flows for outlet control
     297               
     298                # Determine the depth at the outlet relative to the depth of flow in the Culvert
     299                if self.outflow.get_average_height() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
     300                    outlet_culvert_depth=diameter
    301301                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    302302                    perimeter = diameter * pi
    303303                    flow_width= diameter
    304                     case = 'Outlet unsubmerged PIPE FULL'
     304                    case = 'Outlet submerged'
    305305                    if local_debug =='true':
    306                         log.critical('Outlet unsubmerged PIPE FULL')
    307                 else:
    308                     alpha = acos(1-2*outlet_culvert_depth/diameter)*2
    309                     flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
    310                     flow_width= diameter*sin(alpha/2.0)
    311                     perimeter = alpha*diameter/2.0
    312                     case = 'Outlet is open channel flow we will for now assume critical depth'
    313                     if local_debug == 'true':
    314                         log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
    315                                      % (str(Q), str(outlet_culvert_depth),
    316                                         str(alpha)))
    317                         log.critical('Outlet is open channel flow we '
    318                                      'will for now assume critical depth')
    319         if local_debug == 'true':
    320             log.critical('FLOW AREA = %s' % str(flow_area))
    321             log.critical('PERIMETER = %s' % str(perimeter))
    322             log.critical('Q Interim = %s' % str(Q))
    323         hyd_rad = flow_area/perimeter
    324 
    325         if log_filename is not None:
    326             s = 'hydraulic radius at outlet = %f' %hyd_rad
    327             log_to_file(log_filename, s)
    328 
    329         # Outlet control velocity using tail water
    330         if local_debug =='true':
    331             log.critical('GOT IT ALL CALCULATING Velocity')
    332             log.critical('HydRad = %s' % str(hyd_rad))
    333         culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    334         Q_outlet_tailwater = flow_area * culvert_velocity
    335         if local_debug =='true':
    336             log.critical('VELOCITY = %s' % str(culvert_velocity))
    337             log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
    338         if log_filename is not None:               
    339             s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    340             log_to_file(log_filename, s)
    341         Q = min(Q, Q_outlet_tailwater)
    342         if local_debug =='true':
    343             log.critical('%s,%.3f,%.3f'
    344                          % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
    345             log.critical('%s,%.3f,%.3f,%.3f'
    346                          % ('Q and Velocity and Depth=', Q,
    347                             culvert_velocity, outlet_culvert_depth))
    348                            
    349         culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
    350         if local_debug =='true':
    351             log.critical('FLOW AREA = %s' % str(flow_area))
    352             log.critical('PERIMETER = %s' % str(perimeter))
    353             log.critical('Q final = %s' % str(Q))
    354             log.critical('FROUDE = %s' % str(culv_froude))
    355 
    356         # Determine momentum at the outlet
    357         barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
    358 
    359     else: # inlet_depth < 0.01:
    360         Q = barrel_velocity = outlet_culvert_depth = 0.0
    361 
    362     # Temporary flow limit
    363     if barrel_velocity > max_velocity:       
    364         barrel_velocity = max_velocity
    365         Q = flow_area * barrel_velocity
    366    
    367     return Q, barrel_velocity, outlet_culvert_depth
    368 
    369 
     306                        log.critical('Outlet submerged')
     307                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     308                    # IF  self.outflow.get_average_height() < diameter
     309                    dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
     310                    dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
     311                    if dcrit1/diameter >0.85:
     312                        outlet_culvert_depth= dcrit2
     313                    else:
     314                        outlet_culvert_depth = dcrit1
     315                    if outlet_culvert_depth > diameter:
     316                        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     317                        flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
     318                        perimeter = diameter * pi
     319                        flow_width= diameter
     320                        case = 'Outlet unsubmerged PIPE FULL'
     321                        if local_debug =='true':
     322                            log.critical('Outlet unsubmerged PIPE FULL')
     323                    else:
     324                        alpha = acos(1-2*outlet_culvert_depth/diameter)*2
     325                        flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     326                        flow_width= diameter*sin(alpha/2.0)
     327                        perimeter = alpha*diameter/2.0
     328                        case = 'Outlet is open channel flow we will for now assume critical depth'
     329                        if local_debug == 'true':
     330                            log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
     331                                         % (str(Q), str(outlet_culvert_depth),
     332                                            str(alpha)))
     333                            log.critical('Outlet is open channel flow we '
     334                                         'will for now assume critical depth')
     335            if local_debug == 'true':
     336                log.critical('FLOW AREA = %s' % str(flow_area))
     337                log.critical('PERIMETER = %s' % str(perimeter))
     338                log.critical('Q Interim = %s' % str(Q))
     339            hyd_rad = flow_area/perimeter
     340
     341            if self.log_filename is not None:
     342                s = 'hydraulic radius at outlet = %f' %hyd_rad
     343                log_to_file(self.log_filename, s)
     344
     345            # Outlet control velocity using tail water
     346            if local_debug =='true':
     347                log.critical('GOT IT ALL CALCULATING Velocity')
     348                log.critical('HydRad = %s' % str(hyd_rad))
     349            culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
     350            Q_outlet_tailwater = flow_area * culvert_velocity
     351            if local_debug =='true':
     352                log.critical('VELOCITY = %s' % str(culvert_velocity))
     353                log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
     354            if self.log_filename is not None:               
     355                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     356                log_to_file(self.log_filename, s)
     357            Q = min(Q, Q_outlet_tailwater)
     358            if local_debug =='true':
     359                log.critical('%s,%.3f,%.3f'
     360                             % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
     361                log.critical('%s,%.3f,%.3f,%.3f'
     362                             % ('Q and Velocity and Depth=', Q,
     363                                culvert_velocity, outlet_culvert_depth))
     364                               
     365            culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
     366            if local_debug =='true':
     367                log.critical('FLOW AREA = %s' % str(flow_area))
     368                log.critical('PERIMETER = %s' % str(perimeter))
     369                log.critical('Q final = %s' % str(Q))
     370                log.critical('FROUDE = %s' % str(culv_froude))
     371
     372            # Determine momentum at the outlet
     373            barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
     374
     375        else: # self.inflow.get_average_height() < 0.01:
     376            Q = barrel_velocity = outlet_culvert_depth = 0.0
     377
     378        # Temporary flow limit
     379        if barrel_velocity > self.max_velocity:       
     380            barrel_velocity = self.max_velocity
     381            Q = flow_area * barrel_velocity
     382       
     383        return Q, barrel_velocity, outlet_culvert_depth
     384
     385

  • trunk/anuga_core/source/anuga/structures/inlet.py

    r7977 r7978  
    1 
    21from anuga.geometry.polygon import inside_polygon, is_inside_polygon
    32from anuga.config import velocity_protection, g

  • trunk/anuga_core/source/anuga/structures/testing_culvert.py

    r7975 r7978  
    1010import anuga
    1111
    12 from anuga.structures.culvert_operator import Generic_box_culvert
     12from anuga.structures.culvert_operator import Culvert_operator
    1313                           
    1414#from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
    … …  
    8686
    8787filename=os.path.join(path, 'example_rating_curve.csv')
    88 culvert1 = Generic_box_culvert(domain,
    89                               end_point0=[9.0, 2.5],
    90                               end_point1=[13.0, 2.5],
    91                               width=1.00,
    92                               verbose=False)
    93 
     88culvert1 = Culvert_operator(domain,
     89                            end_point0=[9.0, 2.5],
     90                            end_point1=[13.0, 2.5],
     91                            width=1.00,
     92                            verbose=False)
    9493
    9594#culvert2 = Generic_box_culvert(domain,
    … …  
    132131#min_delta_w = sys.maxint
    133132#max_delta_w = -min_delta_w
    134 for t in domain.evolve(yieldstep = 1.0, finaltime = 200):
     133for t in domain.evolve(yieldstep = 1.0, finaltime = 300):
    135134    domain.write_time()
    136135
Note: See TracChangeset for help on using the changeset viewer.