Changeset 7977


Ignore:
Timestamp:
Aug 26, 2010, 5:17:24 PM (14 years ago)
Author:
habili
Message:

Refactoring of routines.

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

Legend:

Unmodified
Added
Removed

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

    r7975 r7977  
    99
    1010
    11 def boyd_box(in):
     11def boyd_box(culvert):
    1212    """Boyd's generalisation of the US department of transportation culvert methods
    1313

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

    r7976 r7977  
    22from anuga.config import g
    33import anuga.utilities.log as log
    4 import inlet
    5 import numpy as num
    6 import math
     4import box_culvert
    75
    8 class Below_interval(Exception): pass
    9 class Above_interval(Exception): pass
    10 
    11 
    12 class Generic_box_culvert:
     6class Culvert_operator:
    137    """Culvert flow - transfer water from one rectangular box to another.
    148    Sets up the geometry of problem
    … …  
    2923                 verbose=False):
    3024       
    31         # Input check
    32        
    3325        self.domain = domain
    34 
    3526        self.domain.set_fractional_step_operator(self)
    36 
    37         self.end_points = [end_point0, end_point1]
     27        end_points = [end_point0, end_point1]
    3828       
    3929        if height is None:
    … …  
    4232        self.width = width
    4333        self.height = height
    44        
    45         self.verbose=verbose
    4634       
    47         # Create the fundamental culvert polygons and create inlet objects
    48         self.create_culvert_polygons()
    49 
    50         #FIXME (SR) Put this into a foe loop to deal with more inlets
    51         self.inlets = []
    52         polygon0 = self.inlet_polygons[0]
    53         inlet0_vector = self.culvert_vector
    54         self.inlets.append(inlet.Inlet(self.domain, polygon0))
    55 
    56         polygon1 = self.inlet_polygons[1]
    57         inlet1_vector = - self.culvert_vector
    58         self.inlets.append(inlet.Inlet(self.domain, polygon1))
     35        self.culvert = box_culvert.Box_culvert(self.domain, end_points, self.width, self.height)
     36        self.inlets = self.culvert.get_inlets()
    5937   
    6038        self.print_stats()
    … …  
    6341    def __call__(self):
    6442
    65         # Time stuff
    66         time     = self.domain.get_time()
    6743        timestep = self.domain.get_timestep()
    68 
    69 
    7044
    7145        inflow  = self.inlets[0]
    … …  
    8155
    8256        delta_z = inflow.get_average_elevation() - outflow.get_average_elevation()
    83         culvert_slope = delta_z/self.culvert_length
     57        culvert_slope = delta_z/self.culvert.get_culvert_length()
    8458
    8559        # Determine controlling energy (driving head) for culvert
    … …  
    9165            driving_head = inflow.get_average_specific_energy()
    9266           
    93 
    9467        # Transfer
    95         from culvert_routines import boyd_generalised_culvert_model
     68        from culvert_routines import boyd_box, boyd_circle
    9669        Q, barrel_velocity, culvert_outlet_depth =\
    97                     boyd_generalised_culvert_model(inflow.get_average_height(),
     70                              boyd_circle(inflow.get_average_height(),
    9871                                         outflow.get_average_height(),
    9972                                         inflow.get_average_speed(),
    … …  
    10174                                         inflow.get_average_specific_energy(),
    10275                                         delta_total_energy,
    103                                          g,
    104                                          culvert_length=self.culvert_length,
     76                                         culvert_length=self.culvert.get_culvert_length(),
    10577                                         culvert_width=self.width,
    10678                                         culvert_height=self.height,
    107                                          culvert_type='box',
    10879                                         manning=0.01)
    10980
    11081        transfer_water = Q*timestep
    111 
    11282
    11383        inflow.set_heights(inflow.get_average_height() - transfer_water)
    … …  
    11585        inflow.set_ymoms(0.0)
    11686
    117 
    11887        outflow.set_heights(outflow.get_average_height() + transfer_water)
    11988        outflow.set_xmoms(0.0)
    12089        outflow.set_ymoms(0.0)
    121 
    122 
    12390
    12491    def print_stats(self):
    … …  
    141108
    142109        print '====================================='
    143 
    144 
    145 
    146 
    147 
    148     def create_culvert_polygons(self):
    149 
    150         """Create polygons at the end of a culvert inlet and outlet.
    151         At either end two polygons will be created; one for the actual flow to pass through and one a little further away
    152         for enquiring the total energy at both ends of the culvert and transferring flow.
    153         """
    154 
    155         # Calculate geometry
    156         x0, y0 = self.end_points[0]
    157         x1, y1 = self.end_points[1]
    158 
    159         dx = x1 - x0
    160         dy = y1 - y0
    161 
    162         self.culvert_vector = num.array([dx, dy])
    163         self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2))
    164         assert self.culvert_length > 0.0, 'The length of culvert is less than 0'
    165 
    166         # Unit direction vector and normal
    167         self.culvert_vector /= self.culvert_length                      # Unit vector in culvert direction
    168         self.culvert_normal = num.array([-dy, dx])/self.culvert_length  # Normal vector
    169 
    170         # Short hands
    171         w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width
    172         h = self.height*self.culvert_vector    # Vector of length=height in the
    173                              # direction of the culvert
    174 
    175         self.inlet_polygons = []
    176 
    177         # Build exchange polygon and enquiry points 0 and 1
    178         for i in [0, 1]:
    179             i0 = (2*i-1)
    180             p0 = self.end_points[i] + w
    181             p1 = self.end_points[i] - w
    182             p2 = p1 + i0*h
    183             p3 = p0 + i0*h
    184             self.inlet_polygons.append(num.array([p0, p1, p2, p3]))
    185 
    186         # Check that enquiry points are outside inlet polygons
    187         for i in [0,1]:
    188             polygon = self.inlet_polygons[i]
    189             # FIXME (SR) Probably should calculate the area of all the triangles
    190             # associated with this polygon, as there is likely to be some
    191             # inconsistency between triangles and ploygon
    192             area = polygon_area(polygon)
    193            
    194             msg = 'Polygon %s ' %(polygon)
    195             msg += ' has area = %f' % area
    196             assert area > 0.0, msg
    197    
    198110
    199111                       

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

    r7939 r7977  
    1 """Collection of culvert routines for use with Culvert_flow in culvert_class
     1"""Collection of culvert routines for use with Culvert_operator
    22
    33This module holds various routines to determine FLOW through CULVERTS and SIMPLE BRIDGES
    … …  
    1818# specific energy is (h + 0.5*v^2/g)
    1919
     20from anuga.config import velocity_protection
     21from anuga.utilities.numerical_tools import safe_acos as acos
    2022
    2123from math import pi, sqrt, sin, cos
    22 
    23 
    24 def boyd_generalised_culvert_model(inlet_depth,
    25                                    outlet_depth,
    26                                    inlet_velocity,
    27                                    outlet_velocity,
    28                                    inlet_specific_energy,
    29                                    delta_total_energy,
    30                                    g,
    31                                    culvert_length=0.0,
    32                                    culvert_width=0.0,
    33                                    culvert_height=0.0,
    34                                    culvert_type='box',
    35                                    manning=0.0,
    36                                    sum_loss=0.0,
    37                                    max_velocity=10.0,
    38                                    log_filename=None):
     24from anuga.config import g
     25
     26def 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):
    3939
    4040    """Boyd's generalisation of the US department of transportation culvert methods
    … …  
    7979    """
    8080
    81     from anuga.utilities.system_tools import log_to_file
    82     from anuga.config import velocity_protection
    83     from anuga.utilities.numerical_tools import safe_acos as acos
    84 
    8581    local_debug ='false'
    8682    if inlet_depth > 0.1: #this value was 0.01:
    … …  
    9995        assert inlet_specific_energy >= 0.0, msg
    10096                     
    101 
    102         if culvert_type =='circle':
    103             # Round culvert (use height as diameter)
    104             diameter = culvert_height   
    105 
    106             """
    107             For a circular pipe the Boyd method reviews 3 conditions
    108             1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
    109             2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
    110             3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
     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
     107            dcrit = (Q**2/g/width**2)**0.333333
     108            if dcrit > height:
     109                dcrit = height
     110            flow_area = width*dcrit
     111            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
    111134           
    112             For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
    113             """
    114            
    115             # Calculate flows for inlet control           
    116             Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
    117             Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
    118             # Note for to SUBMERGED TO OCCUR inlet_specific_energy should be > 1.2 x diameter.... Should Check !!!
    119 
    120             if log_filename is not None:                                       
    121                 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
     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:
     145                    outlet_culvert_depth=height
     146                    flow_area=width*height
     147                    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
    122158                log_to_file(log_filename, s)
    123             Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
    124 
    125             # THE LOWEST Value will Control Calcs From here
    126             # Calculate Critical Depth Based on the Adopted Flow as an Estimate
    127             dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
    128             dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
    129             # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
    130             if dcrit1/diameter  > 0.85:
    131                 outlet_culvert_depth = dcrit2
    132             else:
    133                 outlet_culvert_depth = dcrit1
    134             #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
    135             # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
    136             if outlet_culvert_depth >= diameter:
    137                 outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
    138                 flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    139                 perimeter = diameter * pi
    140                 flow_width= diameter
    141                 case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
    142                 if local_debug == 'true':
    143                     log.critical('Inlet CTRL Outlet submerged Circular '
    144                                  'PIPE FULL')
    145             else:
    146                 #alpha = acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
    147                 alpha = acos(1-2*outlet_culvert_depth/diameter)*2
    148                 #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
    149                 flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
    150                 flow_width= diameter*sin(alpha/2.0)
    151                 perimeter = alpha*diameter/2.0
    152                 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    153                 if local_debug =='true':
    154                     log.critical('INLET CTRL Culvert is open channel flow '
    155                                  'we will for now assume critical depth')
    156                     log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
    157                                  % (str(Q), str(outlet_culvert_depth),
    158                                     str(alpha)))
    159             if delta_total_energy < inlet_specific_energy:  #  OUTLET CONTROL !!!!
    160                 # Calculate flows for outlet control
    161                
    162                 # Determine the depth at the outlet relative to the depth of flow in the Culvert
    163                 if outlet_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
    164                     outlet_culvert_depth=diameter
    165                     flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    166                     perimeter = diameter * pi
    167                     flow_width= diameter
    168                     case = 'Outlet submerged'
    169                     if local_debug =='true':
    170                         log.critical('Outlet submerged')
    171                 else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    172                     # IF  outlet_depth < diameter
    173                     dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)
    174                     dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)
    175                     if dcrit1/diameter >0.85:
    176                         outlet_culvert_depth= dcrit2
    177                     else:
    178                         outlet_culvert_depth = dcrit1
    179                     if outlet_culvert_depth > diameter:
    180                         outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
    181                         flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
    182                         perimeter = diameter * pi
    183                         flow_width= diameter
    184                         case = 'Outlet unsubmerged PIPE FULL'
    185                         if local_debug =='true':
    186                             log.critical('Outlet unsubmerged PIPE FULL')
    187                     else:
    188                         alpha = acos(1-2*outlet_culvert_depth/diameter)*2
    189                         flow_area = diameter**2/8*(alpha - sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
    190                         flow_width= diameter*sin(alpha/2.0)
    191                         perimeter = alpha*diameter/2.0
    192                         case = 'Outlet is open channel flow we will for now assume critical depth'
    193                         if local_debug == 'true':
    194                             log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
    195                                          % (str(Q), str(outlet_culvert_depth),
    196                                             str(alpha)))
    197                             log.critical('Outlet is open channel flow we '
    198                                          'will for now assume critical depth')
    199             if local_debug == 'true':
    200                 log.critical('FLOW AREA = %s' % str(flow_area))
    201                 log.critical('PERIMETER = %s' % str(perimeter))
    202                 log.critical('Q Interim = %s' % str(Q))
    203             hyd_rad = flow_area/perimeter
    204 
    205             if log_filename is not None:
    206                 s = 'hydraulic radius at outlet = %f' %hyd_rad
    207                 log_to_file(log_filename, s)
    208159
    209160            # Outlet control velocity using tail water
    210             if local_debug =='true':
    211                 log.critical('GOT IT ALL CALCULATING Velocity')
    212                 log.critical('HydRad = %s' % str(hyd_rad))
    213161            culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    214162            Q_outlet_tailwater = flow_area * culvert_velocity
    215             if local_debug =='true':
    216                 log.critical('VELOCITY = %s' % str(culvert_velocity))
    217                 log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
    218             if log_filename is not None:               
    219                 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     163
     164            if log_filename is not None:                           
     165                s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
    220166                log_to_file(log_filename, s)
    221167            Q = min(Q, Q_outlet_tailwater)
    222             if local_debug =='true':
    223                 log.critical('%s,%.3f,%.3f'
    224                              % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
    225                 log.critical('%s,%.3f,%.3f,%.3f'
    226                              % ('Q and Velocity and Depth=', Q,
    227                                 culvert_velocity, outlet_culvert_depth))
    228 
    229168        else:
    230169            pass
    231                 #FIXME(Ole): What about inlet control?
    232         # ====  END OF CODE BLOCK FOR "IF" CIRCULAR PIPE
    233 
    234         #else....
    235         if culvert_type == 'box':
    236             if local_debug == 'true':
    237                 log.critical('BOX CULVERT')
    238             # Box culvert (rectangle or square)   ========================================================================================================================
    239 
    240             # Calculate flows for inlet control
    241             height = culvert_height
    242             width = culvert_width
    243             flow_width=culvert_width           
    244            
    245             Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    246             Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    247 
    248             if log_filename is not None:                           
    249                 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
    250                 log_to_file(log_filename, s)
    251 
    252             # FIXME(Ole): Are these functions really for inlet control?   
    253             if Q_inlet_unsubmerged < Q_inlet_submerged:
    254                 Q = Q_inlet_unsubmerged
    255                 dcrit = (Q**2/g/width**2)**0.333333
    256                 if dcrit > height:
    257                     dcrit = height
    258                 flow_area = width*dcrit
    259                 outlet_culvert_depth = dcrit
    260                 case = 'Inlet unsubmerged Box Acts as Weir'
    261             else:   
    262                 Q = Q_inlet_submerged
    263                 flow_area = width*height
    264                 outlet_culvert_depth = height
    265                 case = 'Inlet submerged Box Acts as Orifice'                   
    266 
    267             dcrit = (Q**2/g/width**2)**0.333333
    268 
    269             outlet_culvert_depth = dcrit
    270             if outlet_culvert_depth > height:
    271                 outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
    272                 flow_area = width*height  # Cross sectional area of flow in the culvert
    273                 perimeter = 2*(width+height)
    274                 case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
    275             else:
    276                 flow_area = width * outlet_culvert_depth
    277                 perimeter = width+2*outlet_culvert_depth
    278                 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    279        
    280             if delta_total_energy < inlet_specific_energy:
    281                 # Calculate flows for outlet control
    282                
    283                 # Determine the depth at the outlet relative to the depth of flow in the Culvert
    284                 if outlet_depth > height:        # The Outlet is Submerged
    285                     outlet_culvert_depth=height
    286                     flow_area=width*height       # Cross sectional area of flow in the culvert
    287                     perimeter=2.0*(width+height)
    288                     case = 'Outlet submerged'
    289                 else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    290                     dcrit = (Q**2/g/width**2)**0.333333
    291                     outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
    292                     if outlet_culvert_depth > height:
    293                         outlet_culvert_depth=height
    294                         flow_area=width*height
    295                         perimeter=2.0*(width+height)
    296                         case = 'Outlet is Flowing Full'
    297                     else:
    298                         flow_area=width*outlet_culvert_depth
    299                         perimeter=(width+2.0*outlet_culvert_depth)
    300                         case = 'Outlet is open channel flow'
    301 
    302                 hyd_rad = flow_area/perimeter
    303 
    304                 if log_filename is not None:                               
    305                     s = 'hydraulic radius at outlet = %f' % hyd_rad
    306                     log_to_file(log_filename, s)
    307 
    308                 # Outlet control velocity using tail water
    309                 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333))
    310                 Q_outlet_tailwater = flow_area * culvert_velocity
    311 
    312                 if log_filename is not None:                           
    313                     s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
    314                     log_to_file(log_filename, s)
    315                 Q = min(Q, Q_outlet_tailwater)
    316             else:
    317                 pass
    318                 #FIXME(Ole): What about inlet control?
    319         # ====  END OF CODE BLOCK FOR "IF" BOX
    320 
    321                
    322         # Common code for circle and box geometries ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    323         if log_filename is not None:
    324             log_to_file(log_filename, 'Case: "%s"' % case)
    325            
    326         if log_filename is not None:                       
    327             s = 'Flow Rate Control = %f' % Q
    328             log_to_file(log_filename, s)
     170            #FIXME(Ole): What about inlet control?
    329171
    330172        culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))
    … …  
    334176            log.critical('Q final = %s' % str(Q))
    335177            log.critical('FROUDE = %s' % str(culv_froude))
    336         if log_filename is not None:                           
    337             s = 'Froude in Culvert = %f' % culv_froude
    338             log_to_file(log_filename, s)
    339 
     178 
    340179        # Determine momentum at the outlet
    341180        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
    … …  
    348187    # Temporary flow limit
    349188    if barrel_velocity > max_velocity:
    350         if log_filename is not None:                           
    351             s = 'Barrel velocity was reduced from = %f m/s to %f m/s' % (barrel_velocity, max_velocity)
    352             log_to_file(log_filename, s)
    353        
    354189        barrel_velocity = max_velocity
    355190        Q = flow_area * barrel_velocity
    356191       
    357        
    358 
    359        
    360192    return Q, barrel_velocity, outlet_culvert_depth
    361193
    362194
     195def 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
     285                flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
     286                perimeter = diameter * pi
     287                flow_width= diameter
     288                case = 'Outlet submerged'
     289                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
     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 unsubmerged PIPE FULL'
     305                    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

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

    r7976 r7977  
    142142
    143143
    144 
    145144    def get_average_velocity_head(self):
    146145
    … …  
    158157
    159158
    160 
    161159    def set_heights(self,height):
    162160
    … …  
    167165
    168166        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage)
     167
    169168
    170169    def set_xmoms(self,xmom):
Note: See TracChangeset for help on using the changeset viewer.