Changeset 5429


Ignore:
Timestamp:
Jun 25, 2008, 4:15:10 PM (16 years ago)
Author:
ole
Message:

factored boyd culvert routine into separate function

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/culvert_flow/culvert_boyd_channel.py

    r5428 r5429  
    5353from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
    5454from anuga.utilities.polygon import plot_polygons
     55from anuga.utilities.system_tools import log_to_file
     56
     57
    5558
    5659from math import pi,sqrt
     
    6063#------------------------------------------------------------------------------
    6164
    62 
    63 def log(fid, s):
    64     print s
    65     fid.write(s + '\n')
    66    
    6765
    6866# Open file for storing some specific results...
     
    8987
    9088s='Size'+str(len(domain))
    91 log(fid, s)
    92 
    93 velocity_protection = 1.0e-4
     89log_to_file(fid, s)
     90
     91
     92
     93velocity_protection = 1.0e-5
     94
     95
     96
     97
     98
     99
     100def 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
     210
    94211
    95212
     
    137254               z[i] += embank_hgt-embank_dnslope*(x[i] -12.0)       # Sloping D/S Face
    138255     
    139        
    140         # Constriction
    141         #if 27 < x[i] < 29 and y[i] > 3:
    142         #    z[i] += 2       
    143        
    144         # Pole
    145         #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
    146         #    z[i] += 2
    147 
    148         # HOLE For Pit at Opening[0]
    149         #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2:
    150         #  z[i]-=1
    151 
    152         # HOLE For Pit at Opening[1]
    153         #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2:
    154         #  z[i]-=1
    155 
    156256    return z
    157257
     
    237337                 blockage_topdwn=None,
    238338                 blockage_bottup=None,
     339                 culvert_routine=None,
    239340                 verbose=False):
    240341       
     
    252353        if blockage_topdwn is None: blockage_topdwn=0.00
    253354        if blockage_bottup is None: blockage_bottup=0.00
     355        if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model       
    254356       
    255357
     
    303405        self.blockage_bottup = blockage_bottup
    304406
     407        # Store culvert routine
     408        self.culvert_routine = culvert_routine
     409
    305410       
    306411        # Create objects to update momentum (a bit crude at this stage)
     
    324429        # Print something
    325430        s = 'Culvert Effective Length = %.2f m' %(self.distance)
    326         log(fid, s)
     431        log_to_file(fid, s)
    327432   
    328433        s = 'Culvert Direction is %s\n' %str(self.vector)
    329         log(fid, s)       
     434        log_to_file(fid, s)       
    330435       
    331436    def __call__(self, domain):
     
    340445        delta_t = time-self.last_time
    341446        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
    342         log(fid, s)
     447        log_to_file(fid, s)
    343448       
    344449        msg = 'Time did not advance'
     
    358463            enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i])
    359464           
    360             if self.verbose:
    361                 pass
    362                 #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\
    363                 #      %(i, len(idx), self.enquiry_polygons[i])
    364                
    365465           
    366466            # Get model values for points in enquiry polygon for this opening
     
    434534
    435535            inlet.momentum = self.opening_momentum[0]
    436             outlet.momentum = self.opening_momentum[1]           
     536            outlet.momentum = self.opening_momentum[1]
     537
    437538        else:
    438539            #print 'Flow D/S ---> U/S'
     
    459560
    460561        s = 'Delta total energy = %.3f' %(delta_Et)
    461         log(fid, s)
    462        
    463 
    464         # ====================== NOW ENTER INTO THE CULVERT EQUATIONS AS DEFINED BY BOYD GENERALISED METHOD
    465         # == The quantity of flow passing through a culvert is controlled by many factors
    466         # == 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
    467         # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
    468         # == 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
    469        
    470        
    471         Q_outlet_tailwater = 0.0
    472         inlet.rate = 0.0
    473         outlet.rate = 0.0
    474         Q_inlet_unsubmerged = 0.0
    475         Q_inlet_submerged = 0.0
    476         Q_outlet_critical_depth = 0.0
    477        
    478         if inlet.depth >= 0.01:
    479             # Water has risen above inlet
    480             if self.width==self.height:  # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ???
    481                 pass
    482                 #Q1[t]= 0.421*g**0.5*Diam**0.87*HW**1.63    # Inlet Ctrl Inlet Unsubmerged
    483                 #Q2[t]= 0.530*g**0.5*Diam**1.87*HW**0.63    # Inlet Ctrl Inlet Submerged
    484                 #PipeDcrit=
    485                 #Delta_E=HW-TW
    486             else:
    487                 # Box culvert
     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
    488574               
    489                 sum_loss=self.loss_exit+self.loss_entry+self.loss_bend+self.loss_special
    490                 manning=self.manning
    491 
    492                 # Calculate flows for inlet control
    493                 E = inlet.specific_energy
    494                 #E = min(inlet.specific_energy, delta_Et)
    495 
    496                 s = 'Driving energy  = %f m' %E
    497                 log(fid, s)
    498 
    499                 msg = 'Driving energy is negative'
    500                 assert E >= 0.0, msg
    501                          
    502                 Q_inlet_unsubmerged = 0.540*g**0.5*self.width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    503                 Q_inlet_submerged = 0.702*g**0.5*self.width*self.height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    504 
    505                 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
    506                 log(fid, s)
    507 
    508                 case = ''
    509                 if Q_inlet_unsubmerged < Q_inlet_submerged:
    510                     Q = Q_inlet_unsubmerged
    511                     flow_area = self.width*inlet.depth
    512                     outlet_culv_depth = inlet.depth
    513                     case = 'Inlet unsubmerged'
    514                 else:   
    515                     Q = Q_inlet_submerged
    516                     flow_area = self.width*self.height
    517                     outlet_culv_depth = self.height
    518                     case = 'Inlet submerged'                   
    519 
    520                 if delta_Et < Es:
    521                     # Calculate flows for outlet control
    522                     # Determine the depth at the outlet relative to the depth of flow in the Culvert
    523 
    524                     sum_loss = self.sum_loss
     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
    525588               
    526                     if outlet.depth > self.height:  # The Outlet is Submerged
    527                         outlet_culv_depth=self.height
    528                         flow_area=self.width*self.height # Cross sectional area of flow in the culvert
    529                         perimeter=2.0*(self.width+self.height)
    530                         case = 'Outlet submerged'
    531                     elif outlet.depth==0.0:
    532                         outlet_culv_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
    533                         flow_area=self.width*inlet.depth
    534                         perimeter=(self.width+2.0*inlet.depth)
    535                         case = 'Outlet depth is zero'                       
    536                     else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    537                         outlet_culv_depth=outlet.depth
    538                         flow_area=self.width*outlet.depth
    539                         perimeter=(self.width+2.0*outlet.depth)
    540                         case = 'Outlet is open channel flow'
    541                        
    542                     hyd_rad = flow_area/perimeter
    543                     s = 'hydraulic radius at outlet = %f' %hyd_rad
    544                     log(fid, s)
     589            ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
     590            ymomentum_rate /= delta_t
    545591                   
    546 
    547                     # Outlet control velocity using tail water
    548                     culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333))
    549                     Q_outlet_tailwater = flow_area * culvert_velocity
    550 
    551                     s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
    552                     log(fid, s)
    553 
    554                     Q = min(Q, Q_outlet_tailwater)
    555 
    556 
    557                    
    558                 log(fid, 'Case: "%s"' %case)
     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           
    559617               
    560                 flow_rate_control=Q
    561 
    562                 s = 'Flow Rate Control = %f' %flow_rate_control
    563                 log(fid, s)
    564 
    565                 inlet.rate = -flow_rate_control
    566                 outlet.rate = flow_rate_control               
    567                
    568                 culv_froude=sqrt(flow_rate_control**2*self.width/(g*flow_area**3))
    569                 s = 'Froude in Culvert = %f' %culv_froude
    570                 log(fid, s)
    571 
    572                
    573 
    574                 # Determine momentum at the outlet
    575                 barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
    576                 barrel_momentum = barrel_velocity*outlet_culv_depth
    577                
    578                 outlet_E_Loss= 0.5*0.5*barrel_velocity**2/g   # Ke v^2/2g
    579                 s = 'Barrel velocity = %f' %barrel_velocity
    580                 log(fid, s)
    581 
    582                 # Compute momentum vector at outlet
    583                 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
    584                
    585                 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
    586                 log(fid, s)
    587 
    588                 delta_t = time - self.last_time
    589                 if delta_t > 0.0:
    590                     xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
    591                     xmomentum_rate /= delta_t
    592                    
    593                     ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
    594                     ymomentum_rate /= delta_t
    595                    
    596                     s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
    597                     log(fid, s)                   
    598                 else:
    599                     xmomentum_rate = ymomentum_rate = 0.0
    600 
    601 
    602                 # Set momentum rates for outlet jet
    603                 outlet.momentum[0].rate = xmomentum_rate
    604                 outlet.momentum[1].rate = ymomentum_rate
    605 
    606                 # Remember this value for next step (IMPORTANT)
    607                 outlet.momentum[0].value = outlet_mom_x
    608                 outlet.momentum[1].value = outlet_mom_y                   
    609 
    610                 if int(domain.time*100) % 100 == 0:
    611                     s = 'T=%.5f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
    612                         %(time, flow_rate_control, barrel_velocity)
    613                     s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
    614                          %(outlet_culv_depth, outlet_mom_x,outlet_mom_y)
    615                     s += ' Momentum rate: (%.4f, %.4f)'\
    616                          %(xmomentum_rate, ymomentum_rate)                   
    617                     s+='Outlet Vel= %.3f'\
    618                          %(barrel_velocity)
    619                     log(fid, s)
    620          
     618
     619
    621620       
    622621        # Execute flow term for each opening
     
    650649                       end_point0=[9.0, 2.5],
    651650                       end_point1=[13.0, 2.5],
    652                        width=1.20,height=0.75,                       
     651                       width=1.20,height=0.75,
     652                       culvert_routine=boyd_generalised_culvert_model,       
    653653                       verbose=True)
    654654                       
Note: See TracChangeset for help on using the changeset viewer.