Ignore:
Timestamp:
Apr 16, 2013, 8:04:09 PM (12 years ago)
Author:
steve
Message:

Fixed some problems with dem2pts which turned up in merewether case study

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

Legend:

Unmodified
Added
Removed

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

    r8828 r8832  
    128128
    129129
    130             Q, barrel_velocity, outlet_culvert_depth, case = \
    131                               boyd_box_function(depth               =self.culvert_height,
    132                                                 width               =self.culvert_width,
     130            Q, barrel_velocity, outlet_culvert_depth, flow_area, case = \
     131                              boyd_box_function(width               =self.culvert_width,
     132                                                depth               =self.culvert_height,
    133133                                                flow_width          =self.culvert_width,
    134134                                                length              =self.culvert_length,
    … …  
    159159# define separately so that can be imported in parallel code.
    160160#=============================================================================
    161 def boyd_box_function(depth, width, flow_width, length, driving_energy, delta_total_energy, outlet_enquiry_depth, sum_loss, manning):
     161def boyd_box_function(width, depth, flow_width, length, driving_energy, delta_total_energy, outlet_enquiry_depth, sum_loss, manning):
    162162
    163163    # intially assume the culvert flow is controlled by the inlet
    … …  
    263263    # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
    264264
    265     return Q, barrel_velocity, outlet_culvert_depth, case
    266 
    267 
    268 
    269 
    270 
    271 
    272 
    273 
     265    return Q, barrel_velocity, outlet_culvert_depth, flow_area, case
     266
     267
     268
     269
     270
     271
     272
     273

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

    r8827 r8832  
    33
    44
    5 #=============================================================================
    6 # define separately so that can be imported in parallel code.
    7 #=============================================================================
    8 def 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        
    2025
    2036
    … …  
    21720    mannings_rougness,
    21821    """
    219 
    220     discharge_routine = discharge_routine
    22122
    22223
    … …  
    28687       
    28788
     89
     90    def discharge_routine(self):
     91
     92        local_debug = False
     93
     94        if self.use_velocity_head:
     95            self.delta_total_energy = \
     96                 self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
     97        else:
     98            self.delta_total_energy = \
     99                 self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
     100
     101        self.inflow  = self.inlets[0]
     102        self.outflow = self.inlets[1]
     103
     104        if self.delta_total_energy < 0:
     105            self.inflow  = self.inlets[1]
     106            self.outflow = self.inlets[0]
     107            self.delta_total_energy = -self.delta_total_energy
     108
     109
     110        if self.inflow.get_enquiry_depth() > 0.01: #this value was 0.01:
     111
     112            if local_debug:
     113                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
     114                             % (str(self.inflow.get_enquiry_specific_energy()),
     115                                str(self.delta_total_energy)))
     116                anuga.log.critical('culvert type = %s' % str(culvert_type))
     117            # Water has risen above inlet
     118
     119
     120            msg = 'Specific energy at inlet is negative'
     121            assert self.inflow.get_enquiry_specific_energy() >= 0.0, msg
     122
     123            if self.use_velocity_head :
     124                self.driving_energy = self.inflow.get_enquiry_specific_energy()
     125            else:
     126                self.driving_energy = self.inflow.get_enquiry_depth()
     127
     128
     129
     130            Q, barrel_velocity, outlet_culvert_depth, flow_area, case = \
     131                              boyd_pipe_function(depth               =self.inflow.get_enquiry_depth(),
     132                                                 diameter            =self.culvert_diameter,
     133                                                 length              =self.culvert_length,
     134                                                 driving_energy      =self.driving_energy,
     135                                                 delta_total_energy  =self.delta_total_energy,
     136                                                 outlet_enquiry_depth=self.outflow.get_enquiry_depth(),
     137                                                 sum_loss            =self.sum_loss,
     138                                                 manning             =self.manning)
     139        else:
     140             Q = barrel_velocity = outlet_culvert_depth = 0.0
     141             case = 'Inlet dry'
     142
     143
     144        self.case = case
     145
     146
     147        # Temporary flow limit
     148        if barrel_velocity > self.max_velocity:
     149            barrel_velocity = self.max_velocity
     150            Q = flow_area * barrel_velocity
     151
     152        return Q, barrel_velocity, outlet_culvert_depth
     153
     154#=============================================================================
     155# define separately so that can be imported in parallel code.
     156#=============================================================================
     157def boyd_pipe_function(depth, diameter, length, driving_energy, delta_total_energy, outlet_enquiry_depth, sum_loss, manning):
     158
     159
     160    local_debug = False
     161
     162
     163    """
     164    For a circular pipe the Boyd method reviews 3 conditions
     165    1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)
     166    2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)
     167    3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe
     168
     169    For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe
     170    """
     171
     172    # Calculate flows for inlet control for circular pipe
     173    Q_inlet_unsubmerged = 0.421*anuga.g**0.5*diameter**0.87*driving_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
     174    Q_inlet_submerged = 0.530*anuga.g**0.5*diameter**1.87*driving_energy**0.63   # Inlet Ctrl Inlet Submerged
     175    # Note for to SUBMERGED TO OCCUR operator.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!
     176
     177
     178    Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)
     179
     180    # THE LOWEST Value will Control Calcs From here
     181    # Calculate Critical Depth Based on the Adopted Flow as an Estimate
     182    dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
     183    dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
     184    # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as
     185    if dcrit1/diameter  > 0.85:
     186        outlet_culvert_depth = dcrit2
     187    else:
     188        outlet_culvert_depth = dcrit1
     189    #outlet_culvert_depth = min(outlet_culvert_depth, diameter)
     190    # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter
     191    if outlet_culvert_depth >= diameter:
     192        outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     193        flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
     194        perimeter = diameter * math.pi
     195        flow_width= diameter
     196        case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'
     197        if local_debug:
     198            anuga.log.critical('Inlet CTRL Outlet submerged Circular '
     199                         'PIPE FULL')
     200    else:
     201        #alpha = anuga.acos(1 - outlet_culvert_depth/diameter)    # Where did this Come From ????/
     202        alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
     203        #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))        # Pipe is Running Partly Full at the INLET   WHRE did this Come From ?????
     204        flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     205        flow_width= diameter*math.sin(alpha/2.0)
     206        perimeter = alpha*diameter/2.0
     207        case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'
     208        if local_debug:
     209            anuga.log.critical('INLET CTRL Culvert is open channel flow '
     210                         'we will for now assume critical depth')
     211            anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
     212                         % (str(Q), str(outlet_culvert_depth),
     213                            str(alpha)))
     214    if delta_total_energy < driving_energy:  #  OUTLET CONTROL !!!!
     215        # Calculate flows for outlet control
     216
     217        # Determine the depth at the outlet relative to the depth of flow in the Culvert
     218        if outlet_enquiry_depth > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
     219            outlet_culvert_depth=diameter
     220            flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
     221            perimeter = diameter * math.pi
     222            flow_width= diameter
     223            case = 'Outlet submerged'
     224            if local_debug:
     225                anuga.log.critical('Outlet submerged')
     226        else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     227            # IF  operator.outflow.get_average_depth() < diameter
     228            dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
     229            dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
     230            if dcrit1/diameter >0.85:
     231                outlet_culvert_depth= dcrit2
     232            else:
     233                outlet_culvert_depth = dcrit1
     234            if outlet_culvert_depth > diameter:
     235                outlet_culvert_depth = diameter  # Once again the pipe is flowing full not partfull
     236                flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
     237                perimeter = diameter * math.pi
     238                flow_width= diameter
     239                case = 'Outlet unsubmerged PIPE FULL'
     240                if local_debug:
     241                    anuga.log.critical('Outlet unsubmerged PIPE FULL')
     242            else:
     243                alpha = anuga.acos(1-2*outlet_culvert_depth/diameter)*2
     244                flow_area = diameter**2/8*(alpha - math.sin(alpha))   # Equation from  GIECK 5th Ed. Pg. B3
     245                flow_width= diameter*math.sin(alpha/2.0)
     246                perimeter = alpha*diameter/2.0
     247                case = 'Outlet is open channel flow we will for now assume critical depth'
     248                if local_debug:
     249                    anuga.log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'
     250                                 % (str(Q), str(outlet_culvert_depth),
     251                                    str(alpha)))
     252                    anuga.log.critical('Outlet is open channel flow we '
     253                                 'will for now assume critical depth')
     254    if local_debug:
     255        anuga.log.critical('FLOW AREA = %s' % str(flow_area))
     256        anuga.log.critical('PERIMETER = %s' % str(perimeter))
     257        anuga.log.critical('Q Interim = %s' % str(Q))
     258    hyd_rad = flow_area/perimeter
     259
     260
     261
     262    # Outlet control velocity using tail water
     263    if local_debug:
     264        anuga.log.critical('GOT IT ALL CALCULATING Velocity')
     265        anuga.log.critical('HydRad = %s' % str(hyd_rad))
     266    # Calculate Pipe Culvert Outlet Control Velocity.... May need initial Estimate First ??
     267
     268    culvert_velocity = math.sqrt(delta_total_energy/((sum_loss/2/anuga.g)+\
     269                                                              (manning**2*length)/hyd_rad**1.33333))
     270    Q_outlet_tailwater = flow_area * culvert_velocity
     271
     272
     273    if local_debug:
     274        anuga.log.critical('VELOCITY = %s' % str(culvert_velocity))
     275        anuga.log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))
     276
     277    Q = min(Q, Q_outlet_tailwater)
     278    if local_debug:
     279        anuga.log.critical('%s,%.3f,%.3f'
     280                     % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))
     281        anuga.log.critical('%s,%.3f,%.3f,%.3f'
     282                     % ('Q and Velocity and Depth=', Q,
     283                        culvert_velocity, outlet_culvert_depth))
     284
     285    culv_froude=math.sqrt(Q**2*flow_width/(anuga.g*flow_area**3))
     286    if local_debug:
     287        anuga.log.critical('FLOW AREA = %s' % str(flow_area))
     288        anuga.log.critical('PERIMETER = %s' % str(perimeter))
     289        anuga.log.critical('Q final = %s' % str(Q))
     290        anuga.log.critical('FROUDE = %s' % str(culv_froude))
     291
     292    # Determine momentum at the outlet
     293    barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
     294
     295
     296    return Q, barrel_velocity, outlet_culvert_depth, flow_area, case
     297
     298
     299
     300
Note: See TracChangeset for help on using the changeset viewer.