Changeset 5436


Ignore:
Timestamp:
Jun 25, 2008, 5:57:01 PM (14 years ago)
Author:
ole
Message:

Moved culvert routines into own area

Files:
1 added
1 deleted
6 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/config.py

    r5352 r5436  
    99default_smoothing_parameter = 0.001 # Default alpha for penalised
    1010                                    # least squares fitting
     11
     12velocity_protection = 1.0e-6                                     
    1113
    1214#-------------------------------------------
  • anuga_core/source/anuga/culvert_flows/culvert_routines.py

    r5433 r5436  
     1"""Collection of culvert routines for use with Culvert_flow in culvert_class
     2
     3Usage:
     4
     5   
     6
     7"""
     8
     9#   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
     10#
     11# The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
     12# Flow Is Removed at a Rate of INFLOW
     13# Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
     14#
     15# This SHould be changed to a Vertical Opening Both BOX and Circular
     16# There will be several Culvert Routines such as:
     17# CULVERT_Boyd_Channel
     18# CULVERT_Orifice_and_Weir
     19# CULVERT_Simple_FLOOR
     20# CULVERT_Simple_WALL
     21# CULVERT_Eqn_Floor
     22# CULVERT_Eqn_Wall
     23# CULVERT_Tab_Floor
     24# CULVERT_Tab_Wall
     25# BRIDGES.....
     26# NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
     27
     28#  COULD USE EPA SWMM Model !!!!
     29
     30
     31from math import pi, sqrt, sin, cos
     32
     33
     34def boyd_generalised_culvert_model(culvert, inlet, outlet, delta_Et, g):
     35
     36    """Boyd's generalisation of the US department of transportation culvert model
     37        # == The quantity of flow passing through a culvert is controlled by many factors
     38        # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
     39        # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
     40        # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled
     41    """
     42
     43    from anuga.utilities.system_tools import log_to_file
     44    from anuga.config import velocity_protection
     45    from anuga.utilities.numerical_tools import safe_acos as acos
     46
     47       
     48    Q_outlet_tailwater = 0.0
     49    inlet.rate = 0.0
     50    outlet.rate = 0.0
     51    Q_inlet_unsubmerged = 0.0
     52    Q_inlet_submerged = 0.0
     53    Q_outlet_critical_depth = 0.0
     54
     55    log_filename = culvert.log_filename
     56
     57    manning = culvert.manning
     58    sum_loss = culvert.sum_loss
     59    length = culvert.length
     60   
     61    if inlet.depth >= 0.01:
     62        # Calculate driving energy
     63        E = inlet.specific_energy
     64
     65        s = 'Driving energy  = %f m' %E
     66        log_to_file(log_filename, s)
     67       
     68        msg = 'Driving energy is negative'
     69        assert E >= 0.0, msg
     70                     
     71       
     72        # Water has risen above inlet
     73        if culvert.culvert_type == 'circle':
     74            # Round culvert
     75
     76            # Calculate flows for inlet control           
     77            diameter = culvert.diameter
     78
     79            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*E**1.63    # Inlet Ctrl Inlet Unsubmerged
     80            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**0.63      # Inlet Ctrl Inlet Submerged
     81
     82            s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
     83            log_to_file(log_filename, s)
     84
     85            case = ''
     86            if Q_inlet_unsubmerged < Q_inlet_submerged:
     87                Q = Q_inlet_unsubmerged
     88
     89                alpha = acos(1 - inlet.depth/diameter)
     90                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
     91                outlet_culvert_depth = inlet.depth
     92                width = diameter*sin(alpha)
     93                perimeter = alpha*diameter               
     94                case = 'Inlet unsubmerged'
     95            else:   
     96                Q = Q_inlet_submerged
     97                flow_area = (diameter/2)**2 * pi
     98                outlet_culvert_depth = diameter
     99                width = diameter
     100                perimeter = diameter
     101                case = 'Inlet submerged'                   
     102
     103
     104
     105            if delta_Et < E:
     106                # Calculate flows for outlet control
     107                # Determine the depth at the outlet relative to the depth of flow in the Culvert
     108
     109                if outlet.depth > diameter:        # The Outlet is Submerged
     110                    outlet_culvert_depth=diameter
     111                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
     112                    perimeter = diameter * pi
     113                    width = diameter
     114                    case = 'Outlet submerged'
     115                elif outlet.depth==0.0:
     116                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
     117                    alpha = acos(1 - inlet.depth/diameter)
     118                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
     119                    perimeter = alpha*diameter
     120                    width = diameter*sin(alpha)                   
     121
     122                    case = 'Outlet depth is zero'                       
     123                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     124                    outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
     125                    alpha = acos(1 - outlet.depth/diameter)
     126                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
     127                    perimeter = alpha*diameter
     128                    width = diameter*sin(alpha)                   
     129                    case = 'Outlet is open channel flow'
     130
     131
     132        else:
     133            # Box culvert (rectangle or square)
     134
     135            # Calculate flows for inlet control
     136            height = culvert.height
     137            width = culvert.width           
     138           
     139            Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
     140            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
     141
     142            s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
     143            log_to_file(log_filename, s)
     144
     145            case = ''
     146            if Q_inlet_unsubmerged < Q_inlet_submerged:
     147                Q = Q_inlet_unsubmerged
     148                flow_area = width*inlet.depth
     149                outlet_culvert_depth = inlet.depth
     150                perimeter=(width+2.0*inlet.depth)               
     151                case = 'Inlet unsubmerged'
     152            else:   
     153                Q = Q_inlet_submerged
     154                flow_area = width*height
     155                outlet_culvert_depth = height
     156                perimeter=2.0*(width+height)               
     157                case = 'Inlet submerged'                   
     158
     159            if delta_Et < E:
     160                # Calculate flows for outlet control
     161                # Determine the depth at the outlet relative to the depth of flow in the Culvert
     162
     163                if outlet.depth > height:        # The Outlet is Submerged
     164                    outlet_culvert_depth=height
     165                    flow_area=width*height       # Cross sectional area of flow in the culvert
     166                    perimeter=2.0*(width+height)
     167                    case = 'Outlet submerged'
     168                elif outlet.depth==0.0:
     169                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
     170                    flow_area=width*inlet.depth
     171                    perimeter=(width+2.0*inlet.depth)
     172                    case = 'Outlet depth is zero'                       
     173                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
     174                    outlet_culvert_depth=outlet.depth
     175                    flow_area=width*outlet.depth
     176                    perimeter=(width+2.0*outlet.depth)
     177                    case = 'Outlet is open channel flow'
     178                   
     179
     180
     181
     182        # Common code for rectangular and circular culvert types   
     183        hyd_rad = flow_area/perimeter
     184        s = 'hydraulic radius at outlet = %f' %hyd_rad
     185        log_to_file(log_filename, s)
     186
     187        # Outlet control velocity using tail water
     188        culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
     189        Q_outlet_tailwater = flow_area * culvert_velocity
     190
     191        s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
     192        log_to_file(log_filename, s)
     193        Q = min(Q, Q_outlet_tailwater)
     194
     195        log_to_file(log_filename, 'Case: "%s"' %case)
     196   
     197        flow_rate_control=Q
     198
     199        s = 'Flow Rate Control = %f' %flow_rate_control
     200        log_to_file(log_filename, s)
     201
     202        inlet.rate = -flow_rate_control
     203        outlet.rate = flow_rate_control               
     204           
     205        culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
     206        s = 'Froude in Culvert = %f' %culv_froude
     207        log_to_file(log_filename, s)
     208
     209        # Determine momentum at the outlet
     210        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
     211    else: #inlet.depth < 0.01:
     212        Q = barrel_velocity = outlet_culvert_depth = 0.0
     213       
     214    return Q, barrel_velocity, outlet_culvert_depth
     215
     216
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r5315 r5436  
    15271527                negative values indicate decreases.
    15281528                Rate can be None at initialisation but must be specified
    1529                 before forcting term is applied (i.e. simulation has started).
     1529                before forcing term is applied (i.e. simulation has started).
    15301530
    15311531    center [m]: Coordinates at center of flow point
     
    15651565
    15661566        # Update area if applicable
    1567         self.area = None       
     1567        self.exchange_area = None       
    15681568        if center is not None and radius is not None:
    15691569            assert len(center) == 2
     
    15711571            assert polygon is None, msg
    15721572
    1573             self.area = radius**2*pi
     1573            self.exchange_area = radius**2*pi
    15741574       
    15751575        if polygon is not None:
    1576             self.area = polygon_area(self.polygon)
     1576            self.exchange_area = polygon_area(self.polygon)
    15771577
    15781578
     
    15841584        points = domain.get_centroid_coordinates(absolute=True)
    15851585
    1586         self.indices = None
     1586        # Calculate indices in exchange area for this forcing term
     1587        self.exchange_indices = None
    15871588        if self.center is not None and self.radius is not None:
    15881589            # Inlet is circular
    15891590           
    1590             self.indices = []
     1591            self.exchange_indices = []
    15911592            for k in range(N):
    15921593                x, y = points[k,:] # Centroid
    15931594                if ((x-self.center[0])**2+(y-self.center[1])**2) < self.radius**2:
    1594                     self.indices.append(k)
     1595                    self.exchange_indices.append(k)
    15951596                   
    15961597        if self.polygon is not None:                   
    15971598            # Inlet is polygon
    1598             self.indices = inside_polygon(points, self.polygon)
     1599            self.exchange_indices = inside_polygon(points, self.polygon)
    15991600           
    16001601
     
    16211622
    16221623
    1623         if self.indices is None:
     1624        if self.exchange_indices is None:
    16241625            self.update[:] += rate
    16251626        else:
    16261627            # Brute force assignment of restricted rate
    1627             for k in self.indices:
     1628            for k in self.exchange_indices:
    16281629                self.update[k] += rate
    16291630
     
    16451646        """Return values for specified quantity restricted to opening
    16461647        """
    1647         return self.domain.quantities[self.quantity_name].get_values(indices=self.indices)
     1648        return self.domain.quantities[self.quantity_name].get_values(indices=self.exchange_indices)
    16481649   
    16491650
     
    16511652        """Set values for specified quantity restricted to opening
    16521653        """
    1653         self.domain.quantities[self.quantity_name].set_values(val, indices=self.indices)   
     1654        self.domain.quantities[self.quantity_name].set_values(val, indices=self.exchange_indices)   
    16541655
    16551656
     
    18111812       
    18121813        if callable(self.rate):
    1813             _rate = self.rate(t)/self.area
     1814            _rate = self.rate(t)/self.exchange_area
    18141815        else:
    1815             _rate = self.rate/self.area
     1816            _rate = self.rate/self.exchange_area
    18161817
    18171818        return _rate
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r5304 r5436  
    22312231        # Setup only one forcing term, constant rainfall restricted to a polygon enclosing triangle #1 (bce)
    22322232        domain.forcing_terms = []
    2233         domain.forcing_terms.append( Rainfall(domain, rate=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]]))
     2233        R = Rainfall(domain, rate=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]])
     2234
     2235        assert allclose(R.exchange_area, 1)
     2236       
     2237        domain.forcing_terms.append(R)
    22342238
    22352239        domain.compute_forcing_terms()
  • anuga_core/source/anuga/utilities/system_tools.py

    r5072 r5436  
    66import sys
    77import os
     8
     9def log_to_file(filename, s, verbose=True):
     10    """Log string to open file descriptor
     11    """
     12
     13    fid = open(filename, 'a')
     14    if verbose: print s
     15    fid.write(s + '\n')
     16    fid.close()
     17
    818
    919def get_user_name():
  • anuga_work/development/culvert_flow/culvert_boyd_channel.py

    r5429 r5436  
    5151from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
    5252from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
    53 from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
     53
     54from anuga.culvert_flows.culvert_class import Culvert_flow
     55from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
     56
     57
    5458from anuga.utilities.polygon import plot_polygons
    55 from anuga.utilities.system_tools import log_to_file
    56 
    5759
    5860
     
    6264# Setup computational domain
    6365#------------------------------------------------------------------------------
    64 
    65 
    66 # Open file for storing some specific results...
    67 fid = open('Culvert_Headwall_VarM.txt', 'w')
    6866
    6967length = 40.
     
    8684domain.set_minimum_storable_height(0.001)
    8785
    88 s='Size'+str(len(domain))
    89 log_to_file(fid, s)
    90 
    91 
    92 
    93 velocity_protection = 1.0e-5
    94 
    95 
    96 
    97 
    98 
    99 
    100 def boyd_generalised_culvert_model(inlet, outlet, delta_Et, width, height, length, sum_loss, manning, g, fid):
    101     """Boyd's generalisation of the US department of transportation culvert model
    102         # == The quantity of flow passing through a culvert is controlled by many factors
    103         # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
    104         # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
    105         # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled
    106     """
    107 
    108     from anuga.utilities.system_tools import log_to_file
    109 
    110        
    111     Q_outlet_tailwater = 0.0
    112     inlet.rate = 0.0
    113     outlet.rate = 0.0
    114     Q_inlet_unsubmerged = 0.0
    115     Q_inlet_submerged = 0.0
    116     Q_outlet_critical_depth = 0.0
    117        
    118     if inlet.depth >= 0.01:
    119         # Water has risen above inlet
    120         if width==height:  # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ???
    121             pass
    122             #Q1[t]= 0.421*g**0.5*Diam**0.87*HW**1.63    # Inlet Ctrl Inlet Unsubmerged
    123             #Q2[t]= 0.530*g**0.5*Diam**1.87*HW**0.63    # Inlet Ctrl Inlet Submerged
    124             #PipeDcrit=
    125             #Delta_E=HW-TW
    126         else:
    127             # Box culvert
    128            
    129             # Calculate flows for inlet control
    130             E = inlet.specific_energy
    131 
    132             s = 'Driving energy  = %f m' %E
    133             log_to_file(fid, s)
    134 
    135             msg = 'Driving energy is negative'
    136             assert E >= 0.0, msg
    137                      
    138             Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    139             Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    140 
    141             s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
    142             log_to_file(fid, s)
    143 
    144             case = ''
    145             if Q_inlet_unsubmerged < Q_inlet_submerged:
    146                 Q = Q_inlet_unsubmerged
    147                 flow_area = width*inlet.depth
    148                 outlet_culvert_depth = inlet.depth
    149                 case = 'Inlet unsubmerged'
    150             else:   
    151                 Q = Q_inlet_submerged
    152                 flow_area = width*height
    153                 outlet_culvert_depth = height
    154                 case = 'Inlet submerged'                   
    155 
    156             if delta_Et < E:
    157                 # Calculate flows for outlet control
    158                 # Determine the depth at the outlet relative to the depth of flow in the Culvert
    159 
    160                 if outlet.depth > height:        # The Outlet is Submerged
    161                     outlet_culvert_depth=height
    162                     flow_area=width*height       # Cross sectional area of flow in the culvert
    163                     perimeter=2.0*(width+height)
    164                     case = 'Outlet submerged'
    165                 elif outlet.depth==0.0:
    166                     outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
    167                     flow_area=width*inlet.depth
    168                     perimeter=(width+2.0*inlet.depth)
    169                     case = 'Outlet depth is zero'                       
    170                 else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    171                     outlet_culvert_depth=outlet.depth
    172                     flow_area=width*outlet.depth
    173                     perimeter=(width+2.0*outlet.depth)
    174                     case = 'Outlet is open channel flow'
    175                    
    176                 hyd_rad = flow_area/perimeter
    177                 s = 'hydraulic radius at outlet = %f' %hyd_rad
    178                 log_to_file(fid, s)
    179                
    180 
    181                 # Outlet control velocity using tail water
    182                 culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333))
    183                 Q_outlet_tailwater = flow_area * culvert_velocity
    184 
    185                 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    186                 log_to_file(fid, s)
    187 
    188                 Q = min(Q, Q_outlet_tailwater)
    189 
    190 
    191                
    192             log_to_file(fid, 'Case: "%s"' %case)
    193            
    194             flow_rate_control=Q
    195 
    196             s = 'Flow Rate Control = %f' %flow_rate_control
    197             log_to_file(fid, s)
    198 
    199             inlet.rate = -flow_rate_control
    200             outlet.rate = flow_rate_control               
    201            
    202             culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
    203             s = 'Froude in Culvert = %f' %culv_froude
    204             log_to_file(fid, s)
    205 
    206             # Determine momentum at the outlet
    207             barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
    208 
    209             return Q, barrel_velocity, outlet_culvert_depth
    21086
    21187
     
    264140
    265141
    266 
    267 #------------------------------------------------------------------------------
    268 # Setup specialised forcing terms
    269 #------------------------------------------------------------------------------
    270 
    271  #   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
    272  #
    273  # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
    274  # Flow Is Removed at a Rate of INFLOW
    275  # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
    276  #
    277  # This SHould be changed to a Vertical Opening Both BOX and Circular
    278  # There will be several Culvert Routines such as:
    279  # CULVERT_Boyd_Channel
    280  # CULVERT_Orifice_and_Weir
    281  # CULVERT_Simple_FLOOR
    282  # CULVERT_Simple_WALL
    283  # CULVERT_Eqn_Floor
    284  # CULVERT_Eqn_Wall
    285  # CULVERT_Tab_Floor
    286  # CULVERT_Tab_Wall
    287  # BRIDGES.....
    288  # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
    289  
    290  #  COULD USE EPA SWMM Model !!!!
    291  
    292  
    293  
    294 class Culvert_flow:
    295     """Culvert flow - transfer water from one hole to another
    296    
    297     Using Momentum as Calculated by Culvert Flow !!
    298     Could be Several Methods Investigated to do This !!!
    299 
    300     2008_May_08
    301     To Ole:
    302     OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on
    303     the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon
    304 
    305     The two centers are now passed on to create_polygon.
    306    
    307 
    308     Input: Two points, pipe_size (either diameter or width, height), mannings_rougness,
    309     inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
    310     top-down_blockage_factor and bottom_up_blockage_factor
    311    
    312    
    313     And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront
    314     of the culvert
    315     At the moment this script uses only Depth, later we can change it to Energy...
    316 
    317     Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity
    318     The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
    319        
    320     """
    321 
    322     def __init__(self,
    323                  domain,
    324                  label=None,
    325                  description=None,
    326                  end_point0=None,
    327                  end_point1=None,
    328                  width=None,
    329                  height=None,
    330                  manning=None,          # Mannings Roughness for Culvert
    331                  invert_level0=None,    # Invert level if not the same as the Elevation on the Domain
    332                  invert_level1=None,    # Invert level if not the same as the Elevation on the Domain
    333                  loss_exit=None,
    334                  loss_entry=None,
    335                  loss_bend=None,
    336                  loss_special=None,
    337                  blockage_topdwn=None,
    338                  blockage_bottup=None,
    339                  culvert_routine=None,
    340                  verbose=False):
    341        
    342         from Numeric import sqrt, sum
    343 
    344         # Input check
    345         if height is None: height = width
    346 
    347         # Set defaults
    348         if manning is None: manning = 0.012   # Set a Default Mannings Roughness for Pipe
    349         if loss_exit is None: loss_exit = 1.00
    350         if loss_entry is None: loss_entry = 0.50
    351         if loss_bend is None: loss_bend=0.00
    352         if loss_special is None: loss_special=0.00
    353         if blockage_topdwn is None: blockage_topdwn=0.00
    354         if blockage_bottup is None: blockage_bottup=0.00
    355         if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model       
    356        
    357 
    358         # Create the fundamental culvert polygons from POLYGON
    359         P = create_culvert_polygons(end_point0,
    360                                     end_point1,
    361                                     width=width,   
    362                                     height=height)
    363        
    364         if verbose is True:
    365             pass
    366             #plot_polygons([[end_point0, end_point1],
    367             #               P['exchange_polygon0'],
    368             #               P['exchange_polygon1'],
    369             #               P['enquiry_polygon0'],
    370             #               P['enquiry_polygon1']],
    371             #              figname='culvert_polygon_output')
    372        
    373         self.openings = []
    374         self.openings.append(Inflow(domain,
    375                                     polygon=P['exchange_polygon0']))
    376 
    377         self.openings.append(Inflow(domain,
    378                                     polygon=P['exchange_polygon1']))                                   
    379 
    380 
    381         # Assume two openings for now: Referred to as 0 and 1
    382         assert len(self.openings) == 2
    383        
    384         # Store basic geometry
    385         self.end_points = [end_point0, end_point1]
    386         self.invert_levels = [invert_level0, invert_level1]               
    387         self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
    388         self.vector = P['vector']
    389         self.distance = P['length']; assert self.distance > 0.0
    390         self.verbose = verbose
    391         self.width = width
    392         self.height = height
    393         self.last_time = 0.0       
    394         self.temp_keep_delta_t = 0.0
    395        
    396 
    397         # Store hydraulic parameters
    398         self.manning = manning
    399         self.loss_exit = loss_exit
    400         self.loss_entry = loss_entry
    401         self.loss_bend = loss_bend
    402         self.loss_special = loss_special
    403         self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
    404         self.blockage_topdwn = blockage_topdwn
    405         self.blockage_bottup = blockage_bottup
    406 
    407         # Store culvert routine
    408         self.culvert_routine = culvert_routine
    409 
    410        
    411         # Create objects to update momentum (a bit crude at this stage)
    412 
    413        
    414         xmom0 = General_forcing(domain, 'xmomentum',
    415                                 polygon=P['exchange_polygon0'])
    416 
    417         xmom1 = General_forcing(domain, 'xmomentum',
    418                                 polygon=P['exchange_polygon1'])
    419 
    420         ymom0 = General_forcing(domain, 'ymomentum',
    421                                 polygon=P['exchange_polygon0'])
    422 
    423         ymom1 = General_forcing(domain, 'ymomentum',
    424                                 polygon=P['exchange_polygon1'])
    425 
    426         self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
    427        
    428 
    429         # Print something
    430         s = 'Culvert Effective Length = %.2f m' %(self.distance)
    431         log_to_file(fid, s)
    432    
    433         s = 'Culvert Direction is %s\n' %str(self.vector)
    434         log_to_file(fid, s)       
    435        
    436     def __call__(self, domain):
    437         from anuga.utilities.numerical_tools import mean
    438         from anuga.utilities.polygon import inside_polygon
    439         from anuga.config import g, epsilon
    440         from Numeric import take, sqrt
    441 
    442          
    443         # Time stuff
    444         time = domain.get_time()
    445         delta_t = time-self.last_time
    446         s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    447         log_to_file(fid, s)
    448        
    449         msg = 'Time did not advance'
    450         if time > 0.0: assert delta_t > 0.0, msg
    451 
    452 
    453         # Get average water depths at each opening       
    454         openings = self.openings   # There are two Opening [0] and [1]
    455         for i, opening in enumerate(openings):
    456             stage = domain.quantities['stage'].get_values(location='centroids',
    457                                                           indices=opening.exchange_indices)
    458             elevation = domain.quantities['elevation'].get_values(location='centroids',
    459                                                                   indices=opening.exchange_indices)
    460 
    461             # Indices corresponding to energy enquiry field for this opening
    462             coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
    463             enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i])
    464            
    465            
    466             # Get model values for points in enquiry polygon for this opening
    467             dq = domain.quantities
    468             stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
    469             xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices)
    470             ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices)                       
    471             elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
    472             depth = stage - elevation
    473 
    474             # Compute mean values of selected quantitites in the enquiry area in front of the culvert
    475             # Epsilon handles a dry cell case
    476             ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
    477             uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
    478             v = mean(sqrt(ux**2+uy**2))      # Average velocity
    479             w = mean(stage)                  # Average stage
    480 
    481             # Store values at enquiry field
    482             opening.velocity = v
    483 
    484 
    485             # Compute mean values of selected quantitites in the exchange area in front of the culvert
    486             # Stage and velocity comes from enquiry area and elevation from exchange area
    487            
    488             # Use invert level instead of elevation if specified
    489             invert_level = self.invert_levels[i]
    490             if invert_level is not None:
    491                 z = invert_level
    492             else:
    493                 elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
    494                 z = mean(elevation)                   # Average elevation
    495                
    496             # Estimated depth above the culvert inlet
    497             d = w - z
    498 
    499             if d < 0.0:
    500                 # This is possible since w and z are taken at different locations
    501                 #msg = 'D < 0.0: %f' %d
    502                 #raise msg
    503                 d = 0.0
    504            
    505             # Ratio of depth to culvert height.
    506             # If ratio > 1 then culvert is running full
    507             ratio = d/self.height 
    508             opening.ratio = ratio
    509                
    510             # Average measures of energy in front of this opening
    511             Es = d + 0.5*v**2/g  #  Specific energy in exchange area
    512             Et = w + 0.5*v**2/g  #  Total energy in the enquiry area
    513             opening.total_energy = Et
    514             opening.specific_energy = Es           
    515            
    516             # Store current average stage and depth with each opening object
    517             opening.depth = d
    518             opening.stage = w
    519             opening.elevation = z
    520            
    521 
    522         #################  End of the FOR loop ################################################
    523 
    524            
    525         # We now need to deal with each opening individually
    526 
    527         # Determine flow direction based on total energy difference
    528         delta_Et = openings[0].total_energy - openings[1].total_energy
    529 
    530         if delta_Et > 0:
    531             #print 'Flow U/S ---> D/S'
    532             inlet=openings[0]
    533             outlet=openings[1]
    534 
    535             inlet.momentum = self.opening_momentum[0]
    536             outlet.momentum = self.opening_momentum[1]
    537 
    538         else:
    539             #print 'Flow D/S ---> U/S'
    540             inlet=openings[1]
    541             outlet=openings[0]
    542 
    543             inlet.momentum = self.opening_momentum[1]
    544             outlet.momentum = self.opening_momentum[0]
    545            
    546             delta_Et = -delta_Et
    547 
    548         msg = 'Total energy difference is negative'
    549         assert delta_Et > 0.0, msg
    550 
    551         delta_h = inlet.stage - outlet.stage
    552         delta_z = inlet.elevation - outlet.elevation
    553         culvert_slope = (delta_z/self.distance)
    554 
    555         if culvert_slope < 0.0:
    556             # Adverse gradient - flow is running uphill
    557             # Flow will be purely controlled by uphill outlet face
    558             print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
    559 
    560 
    561         s = 'Delta total energy = %.3f' %(delta_Et)
    562         log_to_file(fid, s)
    563        
    564         sum_loss = self.sum_loss
    565         manning=self.manning
    566         length = self.distance
    567         height = self.height
    568         width = self.width
    569                    
    570         Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(inlet, outlet, delta_Et, width, height, length, sum_loss, manning, g, fid)
    571 
    572         #####################################################
    573         barrel_momentum = barrel_velocity*culvert_outlet_depth
    574                
    575         s = 'Barrel velocity = %f' %barrel_velocity
    576         log_to_file(fid, s)
    577 
    578         # Compute momentum vector at outlet
    579         outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
    580            
    581         s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
    582         log_to_file(fid, s)
    583 
    584         delta_t = time - self.last_time
    585         if delta_t > 0.0:
    586             xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
    587             xmomentum_rate /= delta_t
    588                
    589             ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
    590             ymomentum_rate /= delta_t
    591                    
    592             s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
    593             log_to_file(fid, s)                   
    594         else:
    595             xmomentum_rate = ymomentum_rate = 0.0
    596 
    597 
    598         # Set momentum rates for outlet jet
    599         outlet.momentum[0].rate = xmomentum_rate
    600         outlet.momentum[1].rate = ymomentum_rate
    601 
    602         # Remember this value for next step (IMPORTANT)
    603         outlet.momentum[0].value = outlet_mom_x
    604         outlet.momentum[1].value = outlet_mom_y                   
    605 
    606         if int(domain.time*100) % 100 == 0:
    607             s = 'T=%.5f, Culvert Discharge = %.3f f'\
    608                 %(time, Q)
    609             s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
    610                  %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
    611             s += ' Momentum rate: (%.4f, %.4f)'\
    612                  %(xmomentum_rate, ymomentum_rate)                   
    613             s+='Outlet Vel= %.3f'\
    614                 %(barrel_velocity)
    615             log_to_file(fid, s)
    616            
    617                
    618 
    619 
    620        
    621         # Execute flow term for each opening
    622         # This is where Inflow objects are evaluated and update the domain
    623         for opening in self.openings:
    624             opening(domain)
    625 
    626         # Execute momentum terms
    627         # This is where Inflow objects are evaluated and update the domain
    628         outlet.momentum[0](domain)
    629         outlet.momentum[1](domain)       
    630            
    631         # Store value of time   
    632         self.last_time = time
    633 
    634 
    635 
    636 
    637142#------------------------------------------------------------------------------
    638143# Setup culvert inlets and outlets in current topography
     
    644149# Also Set the Shape & Gap Factors fo rthe Enquiry PolyGons
    645150# ALSO Allow the Invert Level to be provided by the USER
    646 culvert = Culvert_flow(domain,
     151
     152
     153culvert1 = Culvert_flow(domain,
    647154                       label='Culvert No. 1',
    648                        description=' This culvert is a test unit 1.2m Wide by 0.75m High',   
     155                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
    649156                       end_point0=[9.0, 2.5],
    650157                       end_point1=[13.0, 2.5],
     
    652159                       culvert_routine=boyd_generalised_culvert_model,       
    653160                       verbose=True)
     161
     162culvert2 = Culvert_flow(domain,
     163                       label='Culvert No. 2',
     164                       description='This culvert is a circular test with d=1.2m',   
     165                       end_point0=[9.0, 1.5],
     166                       end_point1=[30.0, 3.5],
     167                       diameter=1.20,
     168                       invert_level0=7,
     169                       culvert_routine=boyd_generalised_culvert_model,       
     170                       verbose=True)
    654171                       
    655 domain.forcing_terms.append(culvert)
     172domain.forcing_terms.append(culvert1)
     173domain.forcing_terms.append(culvert2)
    656174
    657175
     
    686204#------------------------------------------------------------------------------
    687205
    688 temp_keep_delta_t=0.0
    689206
    690207for t in domain.evolve(yieldstep = 0.1, finaltime = 20):
Note: See TracChangeset for help on using the changeset viewer.