Changeset 7980


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

Added new tests

Location:
trunk/anuga_core/source/anuga/structures
Files:
5 added
3 edited
2 moved

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
  • trunk/anuga_core/source/anuga/structures/culvert.py

    r7979 r7980  
    66import math
    77
    8 class Box_culvert:
     8class Culvert:
    99    """Culvert flow - transfer water from one rectangular box to another.
    1010    Sets up the geometry of problem
     
    3939
    4040        #FIXME (SR) Put this into a foe loop to deal with more inlets
     41
    4142        self.inlets = []
    4243        polygon0 = self.inlet_polygons[0]
    43         inlet0_vector = self.culvert_vector
    44         self.inlets.append(inlet.Inlet(self.domain, polygon0))
     44        outward_vector0 = self.culvert_vector
     45        self.inlets.append(inlet.Inlet(self.domain, polygon0, outward_vector0))
    4546
    4647        polygon1 = self.inlet_polygons[1]
    47         inlet1_vector = - self.culvert_vector
    48         self.inlets.append(inlet.Inlet(self.domain, polygon1))
    49    
     48        outward_vector1  = - self.culvert_vector
     49        self.inlets.append(inlet.Inlet(self.domain, polygon1, outward_vector1))
     50
     51
     52    def __call__(self):
     53        msg = 'Method __call__ must be overridden by Culvert subclass'
     54        raise Exception, msg
    5055
    5156    def __create_culvert_polygons(self):
     
    7378        # Short hands
    7479        w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
    75         h = self.height*self.culvert_vector    # Vector of length=height in the
     80        h = 0.5*self.height*self.culvert_vector    # Vector of length=height in the
    7681                             # direction of the culvert
    7782
  • trunk/anuga_core/source/anuga/structures/culvert_operator.py

    r7978 r7980  
    22from anuga.config import g
    33import anuga.utilities.log as log
    4 import box_culvert
    5 import culvert_routines
     4
     5from boyd_box_culvert import Boyd_box_culvert
    66
    77class Culvert_operator:
     
    3434        self.height = height
    3535       
    36         self.culvert = box_culvert.Box_culvert(self.domain, end_points, self.width, self.height)
     36        self.culvert = Boyd_box_culvert(self.domain, end_points, self.width, self.height)
     37        print self.culvert
     38        self.routine = self.culvert.routine
     39        print self.routine
    3740        self.inlets = self.culvert.get_inlets()
    3841   
     
    4346
    4447        timestep = self.domain.get_timestep()
     48       
     49        Q, barrel_speed, culvert_outlet_depth = self.routine()
    4550
    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()
     51        inflow  = self.routine.get_inflow()
     52        outflow = self.routine.get_outflow()
    5053
    51         transfer_water = Q*timestep
     54        outflow_direction = - outflow.outward_culvert_vector
    5255
    53         inflow = culvert_routine.get_inflow()
    54         outflow = culvert_routine.get_outflow()
     56        outflow_momentum_flux = barrel_speed**2*culvert_outlet_depth*outflow_direction
    5557
    56         inflow.set_heights(inflow.get_average_height() - transfer_water)
     58
     59        print Q, barrel_speed, culvert_outlet_depth, outflow_momentum_flux
     60
     61        #FIXME (SR) Check whether we need to mult/divide by inlet area
     62        inflow_transfer =  Q*timestep/inflow.get_area()
     63
     64        outflow_transfer = Q*timestep/outflow.get_area()
     65
     66
     67
     68        inflow.set_heights(inflow.get_average_height() - inflow_transfer)
     69
    5770        inflow.set_xmoms(0.0)
    5871        inflow.set_ymoms(0.0)
    5972
    60         outflow.set_heights(outflow.get_average_height() + transfer_water)
    61         outflow.set_xmoms(0.0)
    62         outflow.set_ymoms(0.0)
     73        #u = outflow.get_xvelocities()
     74        #v = outflow.get_yvelocities()
     75
     76        outflow.set_heights(outflow.get_average_height() + outflow_transfer)
     77        #outflow.set_xmoms(outflow.get_xmoms() + timestep*outflow_momentum_flux[0] )
     78        #outflow.set_ymoms(outflow.get_ymoms() + timestep*outflow_momentum_flux[1] )
    6379
    6480    def print_stats(self):
     
    7490
    7591            print 'inlet triangle indices and centres'
    76             print inlet.triangle_indices[i]
    77             print self.domain.get_centroid_coordinates()[inlet.triangle_indices[i]]
     92            print inlet.triangle_indices
     93            print self.domain.get_centroid_coordinates()[inlet.triangle_indices]
    7894       
    7995            print 'polygon'
  • trunk/anuga_core/source/anuga/structures/culvert_routine.py

    r7979 r7980  
    55from anuga.config import g
    66
    7 class Culvert_routines:
     7class Culvert_routine:
    88    """Collection of culvert routines for use with Culvert_operator
    99
     
    2323        self.inlets = culvert.inlets
    2424       
    25         self.inflow  = self.inlets[0]
    26         self.outflow = self.inlets[1]
    27        
     25       
    2826        self.culvert_length = culvert.get_culvert_length()
    2927        self.culvert_width = culvert.get_culvert_width()
     
    3432        self.log_filename = None
    3533       
    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
     34        self.determine_inflow()
    4335
    4436        #delta_z = self.self.inflow.get_average_elevation() - self.self.outflow.get_average_elevation()
     
    5345        #    driving_head = self.inflow.get_average_specific_energy()
    5446           
     47
     48
     49    def determine_inflow(self):
     50        # Determine flow direction based on total energy difference
     51
     52        self.delta_total_energy = self.inlets[0].get_average_total_energy() - self.inlets[1].get_average_total_energy()
     53
     54        self.inflow  = self.inlets[0]
     55        self.outflow = self.inlets[1]
    5556       
     57
     58        if self.delta_total_energy < 0:
     59            self.inflow  = self.inlets[1]
     60            self.outflow = self.inlets[0]
     61            self.delta_total_energy = -self.delta_total_energy
     62
     63
    5664    def get_inflow(self):
    5765       
     
    6371        return self.outflow
    6472   
    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 
    147             dcrit = (Q**2/g/width**2)**0.333333
    148 
    149             outlet_culvert_depth = dcrit
    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
    165                     outlet_culvert_depth=height
    166                     flow_area=width*height       # Cross sectional area of flow in the culvert
    167                     perimeter=2.0*(width+height)
    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
    22173
    22274
    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
    274                 flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    275                 perimeter = diameter * pi
    276                 flow_width= diameter
    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'
    289                 if local_debug =='true':
    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
    301                     flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    302                     perimeter = diameter * pi
    303                     flow_width= diameter
    304                     case = 'Outlet submerged'
    305                     if local_debug =='true':
    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
    38475
    38576
     77
     78
  • trunk/anuga_core/source/anuga/structures/inlet.py

    r7978 r7980  
    99    """
    1010
    11     def __init__(self, domain, polygon):
     11    def __init__(self, domain, polygon, outward_culvert_vector=None):
    1212
    1313        self.domain = domain
    1414        self.domain_bounding_polygon = self.domain.get_boundary_polygon()
    1515        self.polygon = polygon
     16        self.outward_culvert_vector = outward_culvert_vector
    1617
    1718        # FIXME (SR) Using get_triangle_containing_point which needs to be sped up
     
    3536                assert is_inside_polygon(point, bounding_polygon), msg
    3637
    37         self.triangle_indices = inside_polygon(domain_centroids, self.inlet_polygon)
     38
     39        self.triangle_indices = inside_polygon(domain_centroids, self.inlet_polygon, verbose=True)
    3840
    3941        if len(self.triangle_indices) == 0:
     
    4244            msg += 'specified region: %s' % region
    4345            raise Exception, msg
     46
     47
    4448
    4549
     
    6266
    6367
     68    def get_area(self):
     69
     70        return self.area
     71
     72   
    6473    def get_areas(self):
    6574       
     
    7483       
    7584    def get_average_stage(self):
    76        
    77         return num.sum(self.get_stages())/self.triangle_indices.size
    78        
     85
     86        return num.sum(self.get_stages()*self.get_areas())/self.area
    7987       
    8088    def get_elevations(self):   
     
    8492    def get_average_elevation(self):
    8593       
    86         return num.sum(self.get_elevations())/self.triangle_indices.size
     94
     95        return num.sum(self.get_elevations()*self.get_areas())/self.area
    8796   
    8897   
     
    93102       
    94103    def get_average_xmom(self):
    95        
    96         return num.sum(self.get_xmoms())/self.triangle_indices.size
     104
     105        return num.sum(self.get_xmoms()*self.get_areas())/self.area
    97106       
    98107   
     
    104113    def get_average_ymom(self):
    105114       
    106         return num.sum(self.get_ymoms())/self.triangle_indices.size
    107  
     115        return num.sum(self.get_ymoms()*self.get_areas())/self.area
    108116   
    109117    def get_heights(self):
     
    115123       
    116124       return num.sum(self.get_heights()*self.get_areas())
    117    
    118    
     125 
     126
    119127    def get_average_height(self):
    120128   
     
    124132    def get_velocities(self):
    125133       
    126             depths = self.get_stages() - self.get_elevations()
    127             u = self.get_xmoms()/(depths + velocity_protection/depths)
    128             v = self.get_ymoms()/(depths + velocity_protection/depths)
     134            heights = self.get_heights()
     135            u = self.get_xmoms()/(heights + velocity_protection/heights)
     136            v = self.get_ymoms()/(heights + velocity_protection/heights)
    129137           
    130138            return u, v
     139
     140
     141    def get_xvelocities(self):
     142
     143            heights = self.get_heights()
     144            return self.get_xmoms()/(heights + velocity_protection/heights)
     145
     146    def get_yvelocities(self):
     147
     148            heights = self.get_heights()
     149            return self.get_ymoms()/(heights + velocity_protection/heights)
    131150           
    132151           
     
    135154            u, v = self.get_velocities()
    136155           
    137             average_u = num.sum(u)/self.triangle_indices.size
    138             average_v = num.sum(v)/self.triangle_indices.size
     156            average_u = num.sum(u*self.get_areas())/self.area
     157            average_v = num.sum(v*self.get_areas())/self.area
    139158           
    140159            return math.sqrt(average_u**2 + average_v**2)
Note: See TracChangeset for help on using the changeset viewer.