Changeset 8827


Ignore:
Timestamp:
Apr 15, 2013, 8:37:33 PM (12 years ago)
Author:
steve
Message:

changes to make parallel and sequential code reuse

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

Legend:

Unmodified
Added
Removed

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

    r8825 r8827  
    11import anuga
    22import math
     3
     4
     5#=====================================================================
     6# The class
     7#=====================================================================
     8
    39
    410class Boyd_box_operator(anuga.Structure_operator):
    … …  
    8591    def discharge_routine(self):
    8692
    87         local_debug = False
     93
    8894
    8995        if self.use_velocity_head:
    … …  
    102108
    103109
     110
     111
     112        local_debug = False
     113
    104114        if self.inflow.get_enquiry_depth() > 0.01: #this value was 0.01:
     115
    105116            if local_debug:
    106117                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
    … …  
    119130                self.driving_energy = self.inflow.get_enquiry_depth()
    120131
    121             depth = self.culvert_height
    122             width = self.culvert_width
    123             flow_width = self.culvert_width
    124             # intially assume the culvert flow is controlled by the inlet
    125             # check unsubmerged and submerged condition and use Min Q
    126             # but ensure the correct flow area and wetted perimeter are used
    127             Q_inlet_unsubmerged = 0.544*anuga.g**0.5*width*self.driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    128             Q_inlet_submerged = 0.702*anuga.g**0.5*width*depth**0.89*self.driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    129 
    130             # FIXME(Ole): Are these functions really for inlet control?
    131             if Q_inlet_unsubmerged < Q_inlet_submerged:
    132                 Q = Q_inlet_unsubmerged
    133                 dcrit = (Q**2/anuga.g/width**2)**0.333333
    134                 if dcrit > depth:
    135                     dcrit = depth
    136                     flow_area = width*dcrit
    137                     perimeter= 2.0*(width+dcrit)
    138                 else: # dcrit < depth
    139                     flow_area = width*dcrit
    140                     perimeter= 2.0*dcrit+width
    141                 outlet_culvert_depth = dcrit
    142                 self.case = 'Inlet unsubmerged Box Acts as Weir'
    143             else: # Inlet Submerged but check internal culvert flow depth
    144                 Q = Q_inlet_submerged
    145                 dcrit = (Q**2/anuga.g/width**2)**0.333333
    146                 if dcrit > depth:
    147                     dcrit = depth
    148                     flow_area = width*dcrit
    149                     perimeter= 2.0*(width+dcrit)
    150                 else: # dcrit < depth
    151                     flow_area = width*dcrit
    152                     perimeter= 2.0*dcrit+width
    153                 outlet_culvert_depth = dcrit
    154                 self.case = 'Inlet submerged Box Acts as Orifice'
    155 
    156             dcrit = (Q**2/anuga.g/width**2)**0.333333
    157             # May not need this .... check if same is done above
    158             outlet_culvert_depth = dcrit
    159             if outlet_culvert_depth > depth:
    160                 outlet_culvert_depth = depth  # Once again the pipe is flowing full not partfull
    161                 flow_area = width*depth  # Cross sectional area of flow in the culvert
    162                 perimeter = 2*(width+depth)
    163                 self.case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
    164             else:
    165                 flow_area = width * outlet_culvert_depth
    166                 perimeter = width+2*outlet_culvert_depth
    167                 self.case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    168             # Initial Estimate of Flow for Outlet Control using energy slope
    169             #( may need to include Culvert Bed Slope Comparison)
    170             hyd_rad = flow_area/perimeter
    171             culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
    172             Q_outlet_tailwater = flow_area * culvert_velocity
    173            
    174            
    175             if self.delta_total_energy < self.driving_energy:
    176                 # Calculate flows for outlet control
    177 
    178                 # Determine the depth at the outlet relative to the depth of flow in the Culvert
    179                 if self.outflow.get_enquiry_depth() > depth:        # The Outlet is Submerged
    180                     outlet_culvert_depth=depth
    181                     flow_area=width*depth       # Cross sectional area of flow in the culvert
    182                     perimeter=2.0*(width+depth)
    183                     self.case = 'Outlet submerged'
    184                 else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    185                     dcrit = (Q**2/anuga.g/width**2)**0.333333
    186                     outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
    187                     if outlet_culvert_depth > depth:
    188                         outlet_culvert_depth=depth
    189                         flow_area=width*depth
    190                         perimeter=2.0*(width+depth)
    191                         self.case = 'Outlet is Flowing Full'
    192                     else:
    193                         flow_area=width*outlet_culvert_depth
    194                         perimeter=(width+2.0*outlet_culvert_depth)
    195                         self.case = 'Outlet is open channel flow'
    196 
    197                 hyd_rad = flow_area/perimeter
    198 
    199 
    200 
    201                 # Final Outlet control velocity using tail water
    202                 culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
    203                 Q_outlet_tailwater = flow_area * culvert_velocity
    204 
    205                 Q = min(Q, Q_outlet_tailwater)
    206             else:
    207                
    208                 pass
    209                 #FIXME(Ole): What about inlet control?
    210 
    211             culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
    212             if local_debug:
    213                 anuga.log.critical('FLOW AREA = %s' % str(flow_area))
    214                 anuga.log.critical('PERIMETER = %s' % str(perimeter))
    215                 anuga.log.critical('Q final = %s' % str(Q))
    216                 anuga.log.critical('FROUDE = %s' % str(culv_froude))
    217 
    218             # Determine momentum at the outlet
    219             barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
    220 
    221         # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
    222 
    223         else: # self.inflow.get_enquiry_depth() < 0.01:
    224             Q = barrel_velocity = outlet_culvert_depth = 0.0
     132
     133            Q, barrel_velocity, outlet_culvert_depth, case = \
     134                              boyd_box_function(depth               =self.culvert_height,
     135                                                width               =self.culvert_width,
     136                                                flow_width          =self.culvert_width,
     137                                                length              =self.culvert_length,
     138                                                driving_energy      =self.driving_energy,
     139                                                delta_total_energy  =self.delta_total_energy,
     140                                                outlet_enquiry_depth=self.outflow.get_enquiry_depth(),
     141                                                sum_loss            =self.sum_loss,
     142                                                manning             =self.manning)
     143        else:
     144             Q = barrel_velocity = outlet_culvert_depth = 0.0
     145             case = 'Inlet dry'
     146
     147
     148        self.case = case
     149
    225150
    226151        # Temporary flow limit
    … …  
    230155
    231156        return Q, barrel_velocity, outlet_culvert_depth
    232        
    233        
    234        
     157
     158
     159       
     160
     161#=============================================================================
     162# define separately so that can be imported in parallel code.
     163#=============================================================================
     164def boyd_box_function(depth, width, flow_width, length, driving_energy, delta_total_energy, outlet_enquiry_depth, sum_loss, manning):
     165
     166    # intially assume the culvert flow is controlled by the inlet
     167    # check unsubmerged and submerged condition and use Min Q
     168    # but ensure the correct flow area and wetted perimeter are used
     169
     170    local_debug = False
     171   
     172    Q_inlet_unsubmerged = 0.544*anuga.g**0.5*width*driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     173    Q_inlet_submerged = 0.702*anuga.g**0.5*width*depth**0.89*driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
     174
     175    # FIXME(Ole): Are these functions really for inlet control?
     176    if Q_inlet_unsubmerged < Q_inlet_submerged:
     177        Q = Q_inlet_unsubmerged
     178        dcrit = (Q**2/anuga.g/width**2)**0.333333
     179        if dcrit > depth:
     180            dcrit = depth
     181            flow_area = width*dcrit
     182            perimeter= 2.0*(width+dcrit)
     183        else: # dcrit < depth
     184            flow_area = width*dcrit
     185            perimeter= 2.0*dcrit+width
     186        outlet_culvert_depth = dcrit
     187        case = 'Inlet unsubmerged Box Acts as Weir'
     188    else: # Inlet Submerged but check internal culvert flow depth
     189        Q = Q_inlet_submerged
     190        dcrit = (Q**2/anuga.g/width**2)**0.333333
     191        if dcrit > depth:
     192            dcrit = depth
     193            flow_area = width*dcrit
     194            perimeter= 2.0*(width+dcrit)
     195        else: # dcrit < depth
     196            flow_area = width*dcrit
     197            perimeter= 2.0*dcrit+width
     198        outlet_culvert_depth = dcrit
     199        case = 'Inlet submerged Box Acts as Orifice'
     200
     201    dcrit = (Q**2/anuga.g/width**2)**0.333333
     202    # May not need this .... check if same is done above
     203    outlet_culvert_depth = dcrit
     204    if outlet_culvert_depth > depth:
     205        outlet_culvert_depth = depth  # Once again the pipe is flowing full not partfull
     206        flow_area = width*depth  # Cross sectional area of flow in the culvert
     207        perimeter = 2*(width+depth)
     208        case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
     209    else:
     210        flow_area = width * outlet_culvert_depth
     211        perimeter = width+2*outlet_culvert_depth
     212        case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     213    # Initial Estimate of Flow for Outlet Control using energy slope
     214    #( may need to include Culvert Bed Slope Comparison)
     215    hyd_rad = flow_area/perimeter
     216    culvert_velocity = math.sqrt(delta_total_energy/((sum_loss/2/anuga.g) \
     217                                                          +(manning**2*length)/hyd_rad**1.33333))
     218    Q_outlet_tailwater = flow_area * culvert_velocity
     219
     220
     221    if delta_total_energy < driving_energy:
     222        # Calculate flows for outlet control
     223
     224        # Determine the depth at the outlet relative to the depth of flow in the Culvert
     225        if outlet_enquiry_depth > depth:        # The Outlet is Submerged
     226            outlet_culvert_depth=depth
     227            flow_area=width*depth       # Cross sectional area of flow in the culvert
     228            perimeter=2.0*(width+depth)
     229            case = 'Outlet submerged'
     230        else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     231            dcrit = (Q**2/anuga.g/width**2)**0.333333
     232            outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
     233            if outlet_culvert_depth > depth:
     234                outlet_culvert_depth=depth
     235                flow_area=width*depth
     236                perimeter=2.0*(width+depth)
     237                case = 'Outlet is Flowing Full'
     238            else:
     239                flow_area=width*outlet_culvert_depth
     240                perimeter=(width+2.0*outlet_culvert_depth)
     241                case = 'Outlet is open channel flow'
     242
     243        hyd_rad = flow_area/perimeter
     244
     245        # Final Outlet control velocity using tail water
     246        culvert_velocity = math.sqrt(delta_total_energy/((sum_loss/2/anuga.g)\
     247                                                          +(manning**2*length)/hyd_rad**1.33333))
     248        Q_outlet_tailwater = flow_area * culvert_velocity
     249
     250        Q = min(Q, Q_outlet_tailwater)
     251    else:
     252
     253        pass
     254        #FIXME(Ole): What about inlet control?
     255
     256    culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
     257    if local_debug:
     258        anuga.log.critical('FLOW AREA = %s' % str(flow_area))
     259        anuga.log.critical('PERIMETER = %s' % str(perimeter))
     260        anuga.log.critical('Q final = %s' % str(Q))
     261        anuga.log.critical('FROUDE = %s' % str(culv_froude))
     262
     263    # Determine momentum at the outlet
     264    barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
     265
     266    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
     267
     268    return Q, barrel_velocity, outlet_culvert_depth, case
     269
     270
     271
     272
     273
     274
     275
     276

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

    r8825 r8827  
    22import math
    33
     4
     5#=============================================================================
     6# define separately so that can be imported in parallel code.
     7#=============================================================================
     8def discharge_routine(operator):
     9
     10    #operator.__determine_inflow_outflow()
     11
     12    local_debug = False
     13
     14    if operator.use_velocity_head:
     15        operator.delta_total_energy = operator.inlets[0].get_enquiry_total_energy() - operator.inlets[1].get_enquiry_total_energy()
     16    else:
     17        operator.delta_total_energy = operator.inlets[0].get_enquiry_stage() - operator.inlets[1].get_enquiry_stage()
     18
     19    operator.inflow  = operator.inlets[0]
     20    operator.outflow = operator.inlets[1]
     21
     22    if operator.delta_total_energy < 0:
     23        operator.inflow  = operator.inlets[1]
     24        operator.outflow = operator.inlets[0]
     25        operator.delta_total_energy = -operator.delta_total_energy
     26
     27    if operator.inflow.get_enquiry_depth() > 0.01: #this value was 0.01: Remember this needs to be compared to the Invert Lvl
     28        if local_debug:
     29            anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
     30                         % (str(operator.inflow.get_enquiry_specific_energy()),
     31                            str(operator.delta_total_energy)))
     32            anuga.log.critical('culvert type = %s' % str(culvert_type))
     33        # Water has risen above inlet
     34
     35
     36        msg = 'Specific energy at inlet is negative'
     37        assert operator.inflow.get_enquiry_specific_energy() >= 0.0, msg
     38
     39        if operator.use_velocity_head :
     40            operator.driving_energy = operator.inflow.get_enquiry_specific_energy()
     41        else:
     42            operator.driving_energy = operator.inflow.get_enquiry_depth()
     43            """
     44    For a circular pipe the Boyd method reviews 3 conditions
     45    1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
     46    2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
     47    3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
     48
     49    For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
     50    """
     51
     52    diameter = operator.culvert_diameter
     53
     54    if operator.inflow.get_average_depth() > 0.01: #this should test against invert
     55        if local_debug:
     56            anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
     57                         % (str(operator.inflow.get_average_specific_energy()),
     58                            str(operator.delta_total_energy)))
     59            anuga.log.critical('culvert type = %s' % str(culvert_type))
     60        # Water has risen above inlet
     61
     62
     63        msg = 'Specific energy at inlet is negative'
     64        assert operator.inflow.get_average_specific_energy() >= 0.0, msg
     65
     66        # Calculate flows for inlet control for circular pipe
     67        Q_inlet_unsubmerged = 0.421*anuga.g**0.5*diameter**0.87*operator.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged
     68        Q_inlet_submerged = 0.530*anuga.g**0.5*diameter**1.87*operator.inflow.get_average_specific_energy()**0.63   # Inlet Ctrl Inlet Submerged
     69        # Note for to SUBMERGED TO OCCUR operator.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
     70
     71
     72        Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
     73
     74        # THE LOWEST Value will Control Calcs From here
     75        # Calculate Critical Depth Based on the Adopted Flow as an Estimate
     76        dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
     77        dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
     78        # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
     79        if dcrit1/diameter  > 0.85:
     80            outlet_culvert_depth = dcrit2
     81        else:
     82            outlet_culvert_depth = dcrit1
     83        #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
     84        # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
     85        if outlet_culvert_depth >= diameter:
     86            outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     87            flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
     88            perimeter = diameter * math.pi
     89            flow_width= diameter
     90            operator.case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
     91            if local_debug:
     92                anuga.log.critical('Inlet CTRL Outlet submerged Circular '
     93                             'PIPE FULL')
     94        else:
     95            #alpha = anuga.acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
     96            alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
     97            #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
     98            flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     99            flow_width= diameter*math.sin(alpha/2.0)
     100            perimeter = alpha*diameter/2.0
     101            operator.case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     102            if local_debug:
     103                anuga.log.critical('INLET CTRL Culvert is open channel flow '
     104                             'we will for now assume critical depth')
     105                anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
     106                             % (str(Q), str(outlet_culvert_depth),
     107                                str(alpha)))
     108        if operator.delta_total_energy < operator.inflow.get_average_specific_energy():  #  OUTLET CONTROL !!!!
     109            # Calculate flows for outlet control
     110
     111            # Determine the depth at the outlet relative to the depth of flow in the Culvert
     112            if operator.outflow.get_average_depth() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
     113                outlet_culvert_depth=diameter
     114                flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
     115                perimeter = diameter * math.pi
     116                flow_width= diameter
     117                operator.case = 'Outlet submerged'
     118                if local_debug:
     119                    anuga.log.critical('Outlet submerged')
     120            else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     121                # IF  operator.outflow.get_average_depth() < diameter
     122                dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
     123                dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
     124                if dcrit1/diameter >0.85:
     125                    outlet_culvert_depth= dcrit2
     126                else:
     127                    outlet_culvert_depth = dcrit1
     128                if outlet_culvert_depth > diameter:
     129                    outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     130                    flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
     131                    perimeter = diameter * math.pi
     132                    flow_width= diameter
     133                    operator.case = 'Outlet unsubmerged PIPE FULL'
     134                    if local_debug:
     135                        anuga.log.critical('Outlet unsubmerged PIPE FULL')
     136                else:
     137                    alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
     138                    flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     139                    flow_width= diameter*math.sin(alpha/2.0)
     140                    perimeter = alpha*diameter/2.0
     141                    operator.case = 'Outlet is open channel flow we will for now assume critical depth'
     142                    if local_debug:
     143                        anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
     144                                     % (str(Q), str(outlet_culvert_depth),
     145                                        str(alpha)))
     146                        anuga.log.critical('Outlet is open channel flow we '
     147                                     'will for now assume critical depth')
     148        if local_debug:
     149            anuga.log.critical('FLOW AREA = %s' % str(flow_area))
     150            anuga.log.critical('PERIMETER = %s' % str(perimeter))
     151            anuga.log.critical('Q Interim = %s' % str(Q))
     152        hyd_rad = flow_area/perimeter
     153
     154
     155
     156        # Outlet control velocity using tail water
     157        if local_debug:
     158            anuga.log.critical('GOT IT ALL CALCULATING Velocity')
     159            anuga.log.critical('HydRad = %s' % str(hyd_rad))
     160        # Calculate Pipe Culvert Outlet Control Velocity.... May need initial Estimate First ??
     161
     162        culvert_velocity = math.sqrt(operator.delta_total_energy/((operator.sum_loss/2/anuga.g)+\
     163                                                                  (operator.manning**2*operator.culvert_length)/hyd_rad**1.33333))
     164        Q_outlet_tailwater = flow_area * culvert_velocity
     165
     166
     167        if local_debug:
     168            anuga.log.critical('VELOCITY = %s' % str(culvert_velocity))
     169            anuga.log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
     170
     171        Q = min(Q, Q_outlet_tailwater)
     172        if local_debug:
     173            anuga.log.critical('%s,%.3f,%.3f'
     174                         % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
     175            anuga.log.critical('%s,%.3f,%.3f,%.3f'
     176                         % ('Q and Velocity and Depth=', Q,
     177                            culvert_velocity, outlet_culvert_depth))
     178
     179        culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
     180        if local_debug:
     181            anuga.log.critical('FLOW AREA = %s' % str(flow_area))
     182            anuga.log.critical('PERIMETER = %s' % str(perimeter))
     183            anuga.log.critical('Q final = %s' % str(Q))
     184            anuga.log.critical('FROUDE = %s' % str(culv_froude))
     185
     186        # Determine momentum at the outlet
     187        barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
     188
     189    else: # operator.inflow.get_average_depth() < 0.01:
     190        Q = barrel_velocity = outlet_culvert_depth = 0.0
     191
     192    # Temporary flow limit
     193    if barrel_velocity > operator.max_velocity:
     194        barrel_velocity = operator.max_velocity
     195        Q = flow_area * barrel_velocity
     196
     197    return Q, barrel_velocity, outlet_culvert_depth
     198
     199       
     200       
     201       
     202
     203
     204
     205
     206#=====================================================================
     207# Now the class
     208#=====================================================================
    4209class Boyd_pipe_operator(anuga.Structure_operator):
    5210    """Culvert flow - transfer water from one location to another via a circular pipe culvert.
    … …  
    12217    mannings_rougness,
    13218    """
     219
     220    discharge_routine = discharge_routine
     221
    14222
    15223    def __init__(self,
    … …  
    77285        self.case = 'N/A'
    78286       
    79     def discharge_routine(self):
    80 
    81         #self.__determine_inflow_outflow()
    82 
    83         local_debug ='false'
    84        
    85         if self.use_velocity_head:
    86             self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
    87         else:
    88             self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
    89 
    90         self.inflow  = self.inlets[0]
    91         self.outflow = self.inlets[1]
    92 
    93         if self.delta_total_energy < 0:
    94             self.inflow  = self.inlets[1]
    95             self.outflow = self.inlets[0]
    96             self.delta_total_energy = -self.delta_total_energy
    97 
    98         if self.inflow.get_enquiry_depth() > 0.01: #this value was 0.01: Remember this needs to be compared to the Invert Lvl
    99             if local_debug =='true':
    100                 anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
    101                              % (str(self.inflow.get_enquiry_specific_energy()),
    102                                 str(self.delta_total_energy)))
    103                 anuga.log.critical('culvert type = %s' % str(culvert_type))
    104             # Water has risen above inlet
    105 
    106 
    107             msg = 'Specific energy at inlet is negative'
    108             assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
    109 
    110             if self.use_velocity_head :
    111                 self.driving_energy = self.inflow.get_enquiry_specific_energy()
    112             else:
    113                 self.driving_energy = self.inflow.get_enquiry_depth()
    114                 """
    115         For a circular pipe the Boyd method reviews 3 conditions
    116         1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
    117         2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
    118         3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
    119 
    120         For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
    121         """
    122 
    123         diameter = self.culvert_diameter
    124 
    125         local_debug ='false'
    126         if self.inflow.get_average_depth() > 0.01: #this should test against invert
    127             if local_debug =='true':
    128                 anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
    129                              % (str(self.inflow.get_average_specific_energy()),
    130                                 str(self.delta_total_energy)))
    131                 anuga.log.critical('culvert type = %s' % str(culvert_type))
    132             # Water has risen above inlet
    133 
    134 
    135             msg = 'Specific energy at inlet is negative'
    136             assert self.inflow.get_average_specific_energy() >= 0.0, msg
    137 
    138             # Calculate flows for inlet control for circular pipe
    139             Q_inlet_unsubmerged = 0.421*anuga.g**0.5*diameter**0.87*self.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged
    140             Q_inlet_submerged = 0.530*anuga.g**0.5*diameter**1.87*self.inflow.get_average_specific_energy()**0.63   # Inlet Ctrl Inlet Submerged
    141             # Note for to SUBMERGED TO OCCUR self.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
    142 
    143 
    144             Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
    145 
    146             # THE LOWEST Value will Control Calcs From here
    147             # Calculate Critical Depth Based on the Adopted Flow as an Estimate
    148             dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
    149             dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
    150             # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
    151             if dcrit1/diameter  > 0.85:
    152                 outlet_culvert_depth = dcrit2
    153             else:
    154                 outlet_culvert_depth = dcrit1
    155             #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
    156             # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
    157             if outlet_culvert_depth >= diameter:
    158                 outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
    159                 flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
    160                 perimeter = diameter * math.pi
    161                 flow_width= diameter
    162                 self.case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
    163                 if local_debug == 'true':
    164                     anuga.log.critical('Inlet CTRL Outlet submerged Circular '
    165                                  'PIPE FULL')
    166             else:
    167                 #alpha = anuga.acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
    168                 alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
    169                 #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
    170                 flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
    171                 flow_width= diameter*math.sin(alpha/2.0)
    172                 perimeter = alpha*diameter/2.0
    173                 self.case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
    174                 if local_debug =='true':
    175                     anuga.log.critical('INLET CTRL Culvert is open channel flow '
    176                                  'we will for now assume critical depth')
    177                     anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
    178                                  % (str(Q), str(outlet_culvert_depth),
    179                                     str(alpha)))
    180             if self.delta_total_energy < self.inflow.get_average_specific_energy():  #  OUTLET CONTROL !!!!
    181                 # Calculate flows for outlet control
    182 
    183                 # Determine the depth at the outlet relative to the depth of flow in the Culvert
    184                 if self.outflow.get_average_depth() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
    185                     outlet_culvert_depth=diameter
    186                     flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
    187                     perimeter = diameter * math.pi
    188                     flow_width= diameter
    189                     self.case = 'Outlet submerged'
    190                     if local_debug =='true':
    191                         anuga.log.critical('Outlet submerged')
    192                 else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    193                     # IF  self.outflow.get_average_depth() < diameter
    194                     dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
    195                     dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
    196                     if dcrit1/diameter >0.85:
    197                         outlet_culvert_depth= dcrit2
    198                     else:
    199                         outlet_culvert_depth = dcrit1
    200                     if outlet_culvert_depth > diameter:
    201                         outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
    202                         flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
    203                         perimeter = diameter * math.pi
    204                         flow_width= diameter
    205                         self.case = 'Outlet unsubmerged PIPE FULL'
    206                         if local_debug =='true':
    207                             anuga.log.critical('Outlet unsubmerged PIPE FULL')
    208                     else:
    209                         alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
    210                         flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
    211                         flow_width= diameter*math.sin(alpha/2.0)
    212                         perimeter = alpha*diameter/2.0
    213                         self.case = 'Outlet is open channel flow we will for now assume critical depth'
    214                         if local_debug == 'true':
    215                             anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
    216                                          % (str(Q), str(outlet_culvert_depth),
    217                                             str(alpha)))
    218                             anuga.log.critical('Outlet is open channel flow we '
    219                                          'will for now assume critical depth')
    220             if local_debug == 'true':
    221                 anuga.log.critical('FLOW AREA = %s' % str(flow_area))
    222                 anuga.log.critical('PERIMETER = %s' % str(perimeter))
    223                 anuga.log.critical('Q Interim = %s' % str(Q))
    224             hyd_rad = flow_area/perimeter
    225 
    226 
    227 
    228             # Outlet control velocity using tail water
    229             if local_debug =='true':
    230                 anuga.log.critical('GOT IT ALL CALCULATING Velocity')
    231                 anuga.log.critical('HydRad = %s' % str(hyd_rad))
    232             # Calculate Pipe Culvert Outlet Control Velocity.... May need initial Estimate First ??
    233            
    234             culvert_velocity = math.sqrt(self.delta_total_energy/((self.sum_loss/2/anuga.g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))
    235             Q_outlet_tailwater = flow_area * culvert_velocity
    236            
    237            
    238             if local_debug =='true':
    239                 anuga.log.critical('VELOCITY = %s' % str(culvert_velocity))
    240                 anuga.log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
    241 
    242             Q = min(Q, Q_outlet_tailwater)
    243             if local_debug =='true':
    244                 anuga.log.critical('%s,%.3f,%.3f'
    245                              % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
    246                 anuga.log.critical('%s,%.3f,%.3f,%.3f'
    247                              % ('Q and Velocity and Depth=', Q,
    248                                 culvert_velocity, outlet_culvert_depth))
    249 
    250             culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
    251             if local_debug =='true':
    252                 anuga.log.critical('FLOW AREA = %s' % str(flow_area))
    253                 anuga.log.critical('PERIMETER = %s' % str(perimeter))
    254                 anuga.log.critical('Q final = %s' % str(Q))
    255                 anuga.log.critical('FROUDE = %s' % str(culv_froude))
    256 
    257             # Determine momentum at the outlet
    258             barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
    259 
    260         else: # self.inflow.get_average_depth() < 0.01:
    261             Q = barrel_velocity = outlet_culvert_depth = 0.0
    262 
    263         # Temporary flow limit
    264         if barrel_velocity > self.max_velocity:
    265             barrel_velocity = self.max_velocity
    266             Q = flow_area * barrel_velocity
    267 
    268         return Q, barrel_velocity, outlet_culvert_depth
    269 
    270        
    271        
    272        
     287
Note: See TracChangeset for help on using the changeset viewer.