Changeset 5428


Ignore:
Timestamp:
Jun 25, 2008, 3:32:05 PM (15 years ago)
Author:
ole
Message:

First workable version of culvert flows based on Boyd's generalisation of the US department of transportation's culvert model. A momentum 'jet' has been added.
The ANUGA version has the potential improvement of working with the total energy difference and the specific energy at the inlet rather than just the head water.

File:
1 edited

Legend:

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

    r5427 r5428  
    6060#------------------------------------------------------------------------------
    6161
     62
     63def log(fid, s):
     64    print s
     65    fid.write(s + '\n')
     66   
    6267
    6368# Open file for storing some specific results...
     
    8186domain.H0 = 0.01
    8287domain.tight_slope_limiters = True
    83 domain.set_minimum_storable_height(0.02)
     88domain.set_minimum_storable_height(0.001)
    8489
    8590s='Size'+str(len(domain))
    86 fid.write(s)
    87 print s
     91log(fid, s)
     92
     93velocity_protection = 1.0e-4
    8894
    8995
     
    171177 # This SHould be changed to a Vertical Opening Both BOX and Circular
    172178 # There will be several Culvert Routines such as:
     179 # CULVERT_Boyd_Channel
     180 # CULVERT_Orifice_and_Weir
    173181 # CULVERT_Simple_FLOOR
    174182 # CULVERT_Simple_WALL
     
    221229                 height=None,
    222230                 manning=None,          # Mannings Roughness for Culvert
    223                  inv_point0=None,       # Invert level if not the same as the Elevation on the Domain
    224                  inv_point1=None,       # Invert level if not the same as the Elevation on the Domain
     231                 invert_level0=None,    # Invert level if not the same as the Elevation on the Domain
     232                 invert_level1=None,    # Invert level if not the same as the Elevation on the Domain
    225233                 loss_exit=None,
    226234                 loss_entry=None,
     
    233241        from Numeric import sqrt, sum
    234242
     243        # Input check
    235244        if height is None: height = width
     245
     246        # Set defaults
    236247        if manning is None: manning = 0.012   # Set a Default Mannings Roughness for Pipe
    237248        if loss_exit is None: loss_exit = 1.00
     
    241252        if blockage_topdwn is None: blockage_topdwn=0.00
    242253        if blockage_bottup is None: blockage_bottup=0.00
    243         #sum_loss=loss_exit+loss_entry+loss_bend+loss_special    # Note if losses are a function can't do this !!!!
     254       
    244255
    245256        # Create the fundamental culvert polygons from POLYGON
     
    270281       
    271282        # Store basic geometry
    272         self.end_points = [end_point0, end_point1]       
     283        self.end_points = [end_point0, end_point1]
     284        self.invert_levels = [invert_level0, invert_level1]               
    273285        self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
    274286        self.vector = P['vector']
    275         self.distance = P['length']
     287        self.distance = P['length']; assert self.distance > 0.0
    276288        self.verbose = verbose
    277289        self.width = width
    278290        self.height = height
    279291        self.last_time = 0.0       
    280         self.temp_keep_delta_t=0.0
     292        self.temp_keep_delta_t = 0.0
     293       
     294
     295        # Store hydraulic parameters
     296        self.manning = manning
     297        self.loss_exit = loss_exit
     298        self.loss_entry = loss_entry
     299        self.loss_bend = loss_bend
     300        self.loss_special = loss_special
     301        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
     302        self.blockage_topdwn = blockage_topdwn
     303        self.blockage_bottup = blockage_bottup
     304
     305       
    281306        # Create objects to update momentum (a bit crude at this stage)
    282         self.xmom_forcing0 = General_forcing(domain,
    283                                              'xmomentum',
    284                                              polygon=P['exchange_polygon0'])
    285 
    286         self.xmom_forcing1 = General_forcing(domain,
    287                                              'xmomentum',
    288                                              polygon=P['exchange_polygon1'])
    289 
    290         self.ymom_forcing0 = General_forcing(domain,
    291                                              'ymomentum',
    292                                              polygon=P['exchange_polygon0'])
    293 
    294         self.ymom_forcing1 = General_forcing(domain,
    295                                              'ymomentum',
    296                                              polygon=P['exchange_polygon1'])               
     307
     308       
     309        xmom0 = General_forcing(domain, 'xmomentum',
     310                                polygon=P['exchange_polygon0'])
     311
     312        xmom1 = General_forcing(domain, 'xmomentum',
     313                                polygon=P['exchange_polygon1'])
     314
     315        ymom0 = General_forcing(domain, 'ymomentum',
     316                                polygon=P['exchange_polygon0'])
     317
     318        ymom1 = General_forcing(domain, 'ymomentum',
     319                                polygon=P['exchange_polygon1'])
     320
     321        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
     322       
    297323
    298324        # Print something
    299        
    300         s= 'Culvert Effective Length ='+str( self.distance)
    301         print s
    302         fid.write(s)
     325        s = 'Culvert Effective Length = %.2f m' %(self.distance)
     326        log(fid, s)
    303327   
    304         s= 'Culvert Direction is '+str( self.vector)
    305         print s
    306         fid.write(s)
    307         s= 'Point 1m Up Stream is X,Y ='
    308         print s
    309         fid.write(s)
    310         s= 'Point 1m Down Stream is X,Y ='
    311         print s
    312         fid.write(s)
    313        
    314 
     328        s = 'Culvert Direction is %s\n' %str(self.vector)
     329        log(fid, s)       
    315330       
    316331    def __call__(self, domain):
     
    319334        from anuga.config import g, epsilon
    320335        from Numeric import take, sqrt
    321        
    322         #  OLE    I think this routine gets Called twice every time Step, so that each time the second time around Delta_t =0 and therefore Momentum is Set to ZERO
    323                 # Why does this NEED to be calle d twice ???
    324 
    325 
    326         # Get average water depths at each opening
    327 
     336
     337         
     338        # Time stuff
    328339        time = domain.get_time()
    329         openings = self.openings   # There are two Opening  [0] and [1]   so this happens TWICE Right ???????????
     340        delta_t = time-self.last_time
     341        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
     342        log(fid, s)
     343       
     344        msg = 'Time did not advance'
     345        if time > 0.0: assert delta_t > 0.0, msg
     346
     347
     348        # Get average water depths at each opening       
     349        openings = self.openings   # There are two Opening [0] and [1]
    330350        for i, opening in enumerate(openings):
    331351            stage = domain.quantities['stage'].get_values(location='centroids',
    332                                                           indices=opening.indices)
     352                                                          indices=opening.exchange_indices)
    333353            elevation = domain.quantities['elevation'].get_values(location='centroids',
    334                                                                   indices=opening.indices)
     354                                                                  indices=opening.exchange_indices)
    335355
    336356            # Indices corresponding to energy enquiry field for this opening
    337357            coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
    338             idx = inside_polygon(coordinates, self.enquiry_polygons[i])
     358            enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i])
    339359           
    340360            if self.verbose:
     
    344364               
    345365           
    346             # Get average model values for points in ?????? WHAT IF THE DOMAIN IS DRY ???????
    347             # enquiry polygon for this opening
     366            # Get model values for points in enquiry polygon for this opening
    348367            dq = domain.quantities
    349             stage = dq['stage'].get_values(location='centroids', indices=idx)
    350             xmomentum = dq['xmomentum'].get_values(location='centroids', indices=idx)
    351             ymomentum = dq['ymomentum'].get_values(location='centroids', indices=idx)                       
    352             elevation = dq['elevation'].get_values(location='centroids', indices=idx)
    353             depth = stage - elevation           
    354 
    355             # Compute mean velocity in the area in front of the culvert (taking zero depths into account)
    356             ux = xmomentum/(depth+epsilon)   # Does epsilon.... handle a Dry Cell Case ????????
    357             uy = ymomentum/(depth+epsilon)
    358             v = mean(sqrt(ux**2+uy**2))
    359             d = mean(depth)
    360             w = mean(stage)
    361             z = mean(elevation)
    362            
    363             # Ratio of depth to culvert height
    364             #ratio = d/(2*self.height)  # Is Height Now REAL Height or Radius   if Radius then Need 2 else NOT !!
    365             ratio = d/self.height  # Use this if Height is NOT Radius!
    366             if ratio > 1.0:    # Assume culvert is running full &
    367                 ratio = 1.0    # under pressure. Note this is usually ~ 1.35
     368            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
     369            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices)
     370            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices)                       
     371            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
     372            depth = stage - elevation
     373
     374            # Compute mean values of selected quantitites in the enquiry area in front of the culvert
     375            # Epsilon handles a dry cell case
     376            ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
     377            uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
     378            v = mean(sqrt(ux**2+uy**2))      # Average velocity
     379            w = mean(stage)                  # Average stage
     380
     381            # Store values at enquiry field
     382            opening.velocity = v
     383
     384
     385            # Compute mean values of selected quantitites in the exchange area in front of the culvert
     386            # Stage and velocity comes from enquiry area and elevation from exchange area
     387           
     388            # Use invert level instead of elevation if specified
     389            invert_level = self.invert_levels[i]
     390            if invert_level is not None:
     391                z = invert_level
     392            else:
     393                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
     394                z = mean(elevation)                   # Average elevation
     395               
     396            # Estimated depth above the culvert inlet
     397            d = w - z
     398
     399            if d < 0.0:
     400                # This is possible since w and z are taken at different locations
     401                #msg = 'D < 0.0: %f' %d
     402                #raise msg
     403                d = 0.0
     404           
     405            # Ratio of depth to culvert height.
     406            # If ratio > 1 then culvert is running full
     407            ratio = d/self.height 
    368408            opening.ratio = ratio
    369409               
    370 
    371             # Average measure of total energy (D + V^2/2g) in enquiry field in front of this opening
    372             Es = d + 0.5*v**2/g  #  SPECIFIC ENERGY
    373             Et = w + 0.5*v**2/g  #  TOTAL ENERGY
    374             opening.energy = Et
    375             opening.vel=v
    376            
    377 
    378             # Store current average stage and depth at enquiry field with each opening object
     410            # Average measures of energy in front of this opening
     411            Es = d + 0.5*v**2/g  #  Specific energy in exchange area
     412            Et = w + 0.5*v**2/g  #  Total energy in the enquiry area
     413            opening.total_energy = Et
     414            opening.specific_energy = Es           
     415           
     416            # Store current average stage and depth with each opening object
    379417            opening.depth = d
    380418            opening.stage = w
    381419            opening.elevation = z
    382 
    383         #  End of the FOR  LOOP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    384         # we now need to deal with each opening
    385 
    386         # Now we are done calculating energy etc for each opening.
    387         # At this point we can work on the transfer functions etc.
    388        
    389         # Handy values (all calculated at enquiry polygon -
    390         # if you need them at exchange polygons we can easily do that.)
    391 
    392 
    393         #Then we need delta_h or better delta_E to be used to drive water through culvert      IS THIS FROM EXCHANGE or ENQUIRY POLYGON ?????? Need to Know !!!!
    394         #delta_h = openings[0].energy - openings[1].energy
    395         delta_h = openings[0].stage - openings[1].stage
    396         delta_Et = openings[0].energy - openings[1].energy
    397         delta_z = openings[1].elevation - openings[0].elevation
    398         culvert_slope=(delta_z/self.distance)      # COULD USE CULVERT SLOPE To Determine Internal Flow Velocity That is Not Inlet Or Outlet Velocity
    399         s= 'Delta_Energy = %.3f'%(delta_Et)
    400         print s
    401         fid.write(s)
    402        
    403         # Ideas.....  GET RID OF THIS  !!!!!!!!!!!!!!!!
    404         if openings[0].depth > 0 and openings[1].depth > 0:
    405             flow_rate = delta_h
    406         else:
    407             flow_rate = 0.0
    408 
    409         if openings[0].depth > self.height:
    410             # This could be usefull mayby
    411             #print 'Inlet has been overflowed'
    412             pass
    413 
    414         if openings[1].depth > self.height:
    415             # This could be usefull mayby
    416             #print 'Outlet has been overflowed'           
    417             pass
    418         #    ====================== NOW ENTER INTO THE CULVERT EQUATIONS AS DEFINED BY BOYD GENERALISED METHOD
    419         #  ==    The quantity of flow passing through a culvert is controlled by many factors
    420         # ==     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
    421         # ==    Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
    422         # ==    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
    423        
     420           
     421
     422        #################  End of the FOR loop ################################################
     423
     424           
     425        # We now need to deal with each opening individually
     426
     427        # Determine flow direction based on total energy difference
     428        delta_Et = openings[0].total_energy - openings[1].total_energy
     429
    424430        if delta_Et > 0:
    425431            #print 'Flow U/S ---> D/S'
    426432            inlet=openings[0]
    427433            outlet=openings[1]
     434
     435            inlet.momentum = self.opening_momentum[0]
     436            outlet.momentum = self.opening_momentum[1]           
    428437        else:
    429438            #print 'Flow D/S ---> U/S'
     
    431440            outlet=openings[0]
    432441
    433         if inlet.depth < 0.01:   #epsilon:      10mm Minium should be OK....for NOW
    434             # Do nothing because No water over Opening  That is Set Flow to Zero!!
    435             #print 'NO water over Inlet yet!'
    436             self.openings[0].rate = 0
    437             self.openings[1].rate = 0
    438             Qb1=0.0
    439             Qb2=0.0
    440             Qb3tw=0.0
    441             Qb3dc=0.0
    442 
    443 
    444         else: #For PIPES
    445             print 'Water over Inlet NOW!'
     442            inlet.momentum = self.opening_momentum[1]
     443            outlet.momentum = self.opening_momentum[0]
     444           
     445            delta_Et = -delta_Et
     446
     447        msg = 'Total energy difference is negative'
     448        assert delta_Et > 0.0, msg
     449
     450        delta_h = inlet.stage - outlet.stage
     451        delta_z = inlet.elevation - outlet.elevation
     452        culvert_slope = (delta_z/self.distance)
     453
     454        if culvert_slope < 0.0:
     455            # Adverse gradient - flow is running uphill
     456            # Flow will be purely controlled by uphill outlet face
     457            print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
     458
     459
     460        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
    446480            if self.width==self.height:  # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ???
    447481                pass
     
    450484                #PipeDcrit=
    451485                #Delta_E=HW-TW
    452             else:     #FOR BOX CULVERT
    453                         # Why does it not KNOW these Losses from __INIT__  ????
    454                 #sum_loss=loss_exit+loss_entry+loss_bend+loss_special    # Note if losses are a function can't do this !!!!
    455                 sum_loss=1.5
    456                 manning=0.012     # THis should come from Culvert data !!!
    457                 print 'Depth over inlet =',inlet.depth
    458                 Qb1=0.540*g**0.5*self.width*inlet.depth**1.50                    # Flow based on Inlet Ctrl Inlet Unsubmerged
    459                 Qb2=0.702*g**0.5*self.width*self.height**0.89*inlet.depth**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    460                 print'Qb1= %.6f, Qb2 = %.6f'%(Qb1,Qb2)
    461                 Qb_Inlet=min(Qb1,Qb2)
    462                 # NOW FOR OUTLET CONTROL
    463                 # Need to determine the depth at the outlet relative to the depth of flow in the Culvert
    464 #   THIS IS WHERE THE TROUBLE STARTS
    465                 # USUALLY YOU NEED TO ITERATE THECALCULATION OF Q as itis Reliant on Knowing Dcrit which is reliant on Q
    466                 # Instead we will calculate Q assuming TW and then Q again using Dcrit assuming the Qtw
    467                 #If TW >H then Assume Dcrit =H else Assume Dcrit= TW
    468                 if outlet.depth > self.height:  # The Outlet is Submerged
    469                     outlet_culv_depth=self.height
    470                     flow_area=self.width*self.height
    471                     perimeter=2.0*(self.width+self.height)
    472                     hyd_rad=flow_area/perimeter
    473                 elif outlet.depth==0.0:
    474                     outlet_culv_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
    475                     flow_area=self.width*inlet.depth
    476                     perimeter=(self.width+2.0*inlet.depth)
    477                     hyd_rad=flow_area/perimeter
    478                 else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    479                     outlet_culv_depth=outlet.depth
    480                     flow_area=self.width*outlet.depth
    481                     perimeter=(self.width+2.0*outlet.depth)
    482                     hyd_rad=flow_area/perimeter
    483                 print 'hydrad at outlet= ',hyd_rad
    484                 Vel3=sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333))  # Outlet Ctrl usingTW as Control Depth Either TW or H
    485                 Qb3tw=flow_area*Vel3
    486                 Qb3dc=Qhold=Qb_Inlet
    487                 Qcheck=10.0
    488                 while Qcheck > 0.001:
    489                     Dcrit = 0.4672*(Qb3dc/self.width)**0.66666   # This assumes that Flow is at Critical Depth  Not Alsways True...
    490                                         # Better to Use Mannings Equation to Estimate Q and iterate Depth in the Culvert to Match the Energy Loss
    491                     print 'Dcrit =',Dcrit
    492                     if Dcrit > self.height:  #
    493                         flow_area=self.width*self.height
     486            else:
     487                # Box culvert
     488               
     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
     525               
     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
    494529                        perimeter=2.0*(self.width+self.height)
    495                         hyd_rad=flow_area/perimeter
    496                     else:
    497                         flow_area=self.width*Dcrit
    498                         perimeter=(self.width+2.0*Dcrit)
    499                         hyd_rad=flow_area/perimeter
    500                     print 'Adjusted hydrad at outlet= ',hyd_rad
    501                
    502                     Qb3dc=flow_area*sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333))  # Outlet Ctrl using Dcrit as Control Depth Either Dcrit or H
    503                     Qcheck =abs(Qhold-Qb3dc)
    504                     Qhold=Qb3dc
    505                     print 'Qcheck Qb3dc  ',Qcheck,Qb3dc
    506                 flow_rate_control=min(Qb1,Qb2) #
    507                 print 'OutletCtrl Qbtw and Qbdc= ',Qb3tw,Qb3dc
    508                 print 'Flow Rate Control = ',flow_rate_control
    509                                
    510                                 # NOW UPDATE Flows on the DOMAIN once again have to check direction of FLOW
    511                 if delta_Et > 0:
    512                     print 'Flow U/S ---> D/S'
    513                     self.openings[0].rate=-flow_rate_control
    514                     self.openings[1].rate=flow_rate_control
    515                 else:
    516                     print 'Flow D/S ---> U/S'
    517                     self.openings[1].rate=-flow_rate_control
    518                     self.openings[0].rate=flow_rate_control     # Flow Rates have been set  NOW  Add jet in the form of absolute momentum to outlet
    519                                        
    520                                        
     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)
     545                   
     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)
     559               
     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               
    521568                culv_froude=sqrt(flow_rate_control**2*self.width/(g*flow_area**3))
    522                 print 'Froude in Culvert = ',culv_froude
    523 
    524                 #+++++++++++ NOW DETERMINE MOMENTUM AT THE OUTLET !!!! +++++++++++++++
    525 
    526 
    527                 # This Should be Based on the VELOCITY in the Culvert
    528                 # Based on the Flow Depth in the Culvert for part full  So need to determine actual flow depth in the culvert based on slope !!!
    529                 # or Flow Jet of the Pressurised Culvert if FUll
    530                 # If Part Full Flow Calculate Part Full Depth
    531                 # Based on Depth Calculate Area... the Vel = flow_rate_control / Area
    532                 # Need to Break Velocity into X & Y Components
    533                 #self.openings[1].set_quantity_values(delta_h*speed, 'xmomentum')
    534                 # Previous Calculated Depth/ Culvert Heigh Ratio Use it to Determine Velocity !!!! FIX LAter
    535                
    536                 #if (ratio*self.area)== 0: # Don't Really need this already established water depth here
    537                 #   outlet_vel=0.0
    538                 #else:
    539 
    540 
    541                 # FIXME (Ole): Shouldn't this be openings[1]??
    542                 # This attempts to resolve the Velocity in the Culvert at the Entrance and Outlet.... Again if flow reversal Occurs  then need to swap END !!!
    543                 #inlet_vel =(flow_rate_control/(openings[0].ratio*openings[0].area))
    544                 #inlet_dep = openings[0].depth*openings[0].ratio #?????
    545                 #inlet_mom = outlet_vel*outlet_dep
    546                 #inlet_E_Loss= 0.5*0.5*inlet_vel**2/g
    547                 #  NEED to Calculate Critical depth inthe Culvert for Box & Pipe Sections.. If Greater than Depth on the Down Stream Flood Plain use Dcrit
    548                 outlet_vel =(flow_rate_control/flow_area)
    549                 outlet_dep = flow_area/self.width #?????
    550                 outlet_mom = outlet_vel*outlet_dep
    551                 outlet_E_Loss= 0.5*0.5*outlet_vel**2/g   # Ke v^2/2g
    552                 print 'Outlet Velocity =',outlet_vel
    553 
    554                 #barrel_vel=(inlet_vel+outlet_vel)/2
    555                 # use this to determine Barrel Loss
    556                
    557                 # Eventually will Need Momentum in X & Y Components based on the orientation of
    558                 # the culvert from X0,Y0, X1, Y1 from Create Polygon Routine
    559                 # YES - use self.vector which is a unit vector in the direction of the culvert.
    560                                 # Will also need to Check Normal Depth For a Steep Culvert using the Mannings Equation to determine the Flow Depth in the Culvert this requires Iteration and a Simpsons Rule Approach to Refine Depth
    561                                 #    So that Qcalculated = Qestimated....
    562                    
    563                 outlet_mom_x, outlet_mom_y = self.vector * outlet_mom
    564                 #print 'Directional momentum', outlet_mom_x, outlet_mom_y
    565                 xmomentum_rate = outlet_mom_x
    566                 ymomentum_rate = outlet_mom_y
    567                 print 'X Y MOM ',xmomentum_rate,ymomentum_rate
    568                 # Update the momentum forcing terms with directional momentum at the outlet
    569                                 # >>>>>>>>>>>>>>OLE Passing through this the 2nd time ALWAYS sets delta_t to 0 so that Momentum is also set to 0.0    This is NOT RIGHT !!!!  THIS SEGMENT OF CODE NEEDS ONLY TO BE RUn FOR THE OUTLET !!!
    570                                 # ONLY NEEDS TO BE CALLED ONCE for the CULVERT not TWICE for Each Opening as a check of each opening is DONE HERe any WAY !!!
    571                                
    572                                 # I have tried to trick it by introducing this term temp_keep_Delta_t...... but this is really bad !!!      How do we avoid this Situation ????
     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
    573588                delta_t = time - self.last_time
    574                 if delta_t==0:
    575                     delta_t=self.temp_keep_delta_t
    576                     xmomentum_rate=self.xmom_forcing1.rate
    577                     ymomentum_rate=self.ymom_forcing1.rate
    578                 self.temp_keep_delta_t=delta_t
    579                 print 'time and delta t = ',time,delta_t
    580589                if delta_t > 0.0:
    581                     print 'delta_t Calc MOM'
    582                     xmomentum_rate = outlet_mom_x - self.xmom_forcing1.value
     590                    xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
    583591                    xmomentum_rate /= delta_t
    584592                   
    585                     ymomentum_rate = outlet_mom_y - self.ymom_forcing1.value
     593                    ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
    586594                    ymomentum_rate /= delta_t
    587                     print 'X Y MOM_RATE ',xmomentum_rate,ymomentum_rate
    588                     #ymomentum_rate = 3
    589                 #else:
    590                 #    xmomentum_rate = ymomentum_rate = 0.0
     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
    591601
    592602                # Set momentum rates for outlet jet
    593                 self.xmom_forcing1.rate = xmomentum_rate
    594                 self.ymom_forcing1.rate = ymomentum_rate
    595 
    596                 # Remember this value for next step (IMPORTANT)   
    597                 self.xmom_forcing1.value = outlet_mom_x
    598                 self.ymom_forcing1.value = outlet_mom_y                   
    599 
     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                   
    600609
    601610                if int(domain.time*100) % 100 == 0:
    602611                    s = 'T=%.5f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
    603                         %(time, flow_rate_control, outlet_vel)
     612                        %(time, flow_rate_control, barrel_velocity)
    604613                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
    605                          %(outlet_dep, outlet_mom_x,outlet_mom_y)
     614                         %(outlet_culv_depth, outlet_mom_x,outlet_mom_y)
    606615                    s += ' Momentum rate: (%.4f, %.4f)'\
    607616                         %(xmomentum_rate, ymomentum_rate)                   
    608                     s+='Outlet Vel= %.3f Outlet Dep = %.3f'\
    609                          %(outlet_vel,outlet_dep)
    610                     fid.write(s)
    611                     print s
     617                    s+='Outlet Vel= %.3f'\
     618                         %(barrel_velocity)
     619                    log(fid, s)
    612620         
    613621       
     
    619627        # Execute momentum terms
    620628        # This is where Inflow objects are evaluated and update the domain
    621         self.xmom_forcing0(domain)
    622         self.ymom_forcing0(domain)               
    623         self.xmom_forcing1(domain)
    624         self.ymom_forcing1(domain)       
    625 
    626            
    627 
    628         # Print out flows every 1 seconds
    629         if int(time*100) % 100 == 0:
    630             s = 'Time %.6f\n' %time
    631             s += '    Opening[0]: d=%.3f, v=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
    632                  %(openings[0].depth,
    633                    openings[0].vel,
    634                    openings[0].area,
    635                    openings[0].energy,
    636                    openings[0].ratio)                   
    637             s += '    Opening[1]: d=%.3f, v=%.3f,  A=%.3f, E=%.3f, r=%.3f\n'\
    638                  %(openings[1].depth,
    639                    openings[1].vel,
    640                    openings[1].area,
    641                    openings[1].energy,
    642                    openings[1].ratio)                 
    643             s += '    Distance=%.2f, Qb1=%.3f, Qb2=%.3f, Qb3tw=%.3f,Qb3dc=%.3f\n'\
    644                  %(self.distance,
    645                    Qb1,
    646                    Qb2,
    647                    Qb3tw,Qb3dc)
    648            
    649             print s
    650             fid.write(s)
    651 
     629        outlet.momentum[0](domain)
     630        outlet.momentum[1](domain)       
     631           
    652632        # Store value of time   
    653633        self.last_time = time
    654         #self.temp_keep_delta_t=delta_t
     634
     635
     636
     637
    655638#------------------------------------------------------------------------------
    656639# Setup culvert inlets and outlets in current topography
Note: See TracChangeset for help on using the changeset viewer.