""" Ole Check Culvert Routine from Line 258 Although it is Setup as a Culvert with the Opening presented vertically, for now the calculation of flow rate is assuming a horizontal hole in the ground (Fix this Later) MOST importantly 2 things... 1. How to use the Create Polygon Routine to enquire Depth ( or later energy) infront of the Culvert Done (Ole) 2. How to apply the Culvert velocity and thereby Momentum to the Outlet Ject presented at the Outlet Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction This Version CALCULATES the Culvert Velocity and Uses it to establish The Culvert Outlet Momentum The Aim is to define a Flow Transfer function that Simulates a Culvert by using the Total Available Energy to Drive the Culvert as Derived by determining the Difference in Total Energy between 2 Points, Just Up stream and Just Down Stream of the Culvert away from the influence of the Flow Abstraction etc.. This example includes a Model Topography that shows a TYPICAL Headwall Configuration The aim is to change the Culvert Routine to Model more precisely the abstraction from a vertical face. The inflow must include the impact of Approach velocity. Similarly the Outflow has MOMENTUM Not just Up welling as in the Horizontal Style abstraction """ #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross from anuga.shallow_water import Domain from anuga.shallow_water.shallow_water_domain import Reflective_boundary from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing from anuga.culvert_flows.culvert_polygons import create_culvert_polygons from anuga.utilities.polygon import plot_polygons from math import pi,sqrt #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ # Open file for storing some specific results... fid = open('Culvert_Headwall_VarM.txt', 'w') length = 40. width = 5. #dx = dy = 1 # Resolution: Length of subdivisions on both axes #dx = dy = .5 # Resolution: Length of subdivisions on both axes dx = dy = .25 # Resolution: Length of subdivisions on both axes #dx = dy = .1 # Resolution: Length of subdivisions on both axes # OLE.... How do I refine the resolution... in the area where I have the Culvert Opening ???...... # Can I refine in a X & Y Range ??? points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), len1=length, len2=width) domain = Domain(points, vertices, boundary) domain.set_name('culvert_HW_Var_Mom') # Output name domain.set_default_order(2) domain.H0 = 0.01 domain.tight_slope_limiters = True domain.set_minimum_storable_height(0.02) print 'Size', len(domain) #------------------------------------------------------------------------------ # Setup initial conditions #------------------------------------------------------------------------------ # Define the topography (land scape) def topography(x, y): """Set up a weir A culvert will connect either side of an Embankment with a Headwall type structure The aim is for the Model to REALISTICALY model flow through the Culvert """ # General Slope of Topography z=-x/100 # Add bits and Pieces to topography N = len(x) for i in range(N): # Sloping Embankment Across Channel if 5.0 < x[i] < 10.1: if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: # Cut Out Segment for Culvert FACE z[i]=z[i] else: z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face if 10.0 < x[i] < 12.1: z[i] += 2.5 # Flat Crest of Embankment if 12.0 < x[i] < 14.5: if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: # Cut Out Segment for Culvert FACE z[i]=z[i] else: z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face # Constriction #if 27 < x[i] < 29 and y[i] > 3: # z[i] += 2 # Pole #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2: # z[i] += 2 # HOLE For Pit at Opening[0] #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2: # z[i]-=1 # HOLE For Pit at Opening[1] #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2: # z[i]-=1 return z domain.set_quantity('elevation', topography) # Use function for elevation domain.set_quantity('friction', 0.01) # Constant friction domain.set_quantity('stage', expression='elevation') # Dry initial condition #------------------------------------------------------------------------------ # Setup specialised forcing terms #------------------------------------------------------------------------------ # NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet) # # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed # Flow Is Removed at a Rate of INFLOW # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges # # This SHould be changed to a Verical Opening Both BOX and Circular # There will be several Culvert Routines such as: # CULVERT_Simple_FLOOR # CULVERT_Simple_WALL # CULVERT_Eqn_Floor # CULVERT_Eqn_Wall # CULVERT_Tab_Floor # CULVERT_Tab_Wall # BRIDGES..... # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!! # COULD USE EPA SWMM Model !!!! class Culvert_flow: """Culvert flow - transfer water from one hole to another Using Momentum as Calculated by Culvert Flow !! Could be Several Methods Investigated to do This !!! 2008_May_08 To Ole: OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon The two centers are now passed on to create_polygon. Input: Two points, pipe_size (either diameter or width, height), mannings_rougness, inlet/outlet energy_loss_coefficients, internal_bend_coefficent, top-down_blockage_factor and bottom_up_blockage_factor And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront of the culvert At the moment this script uses only Depth, later we can change it to Energy... Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet... """ def __init__(self, domain, end_point0=None, end_point1=None, width=None, height=None, verbose=False): from Numeric import sqrt, sum if height is None: height = width # Create the fundamental culvert polygons P = create_culvert_polygons(end_point0, end_point1, width=width, height=height) if verbose is True: pass #plot_polygons([[end_point0, end_point1], # P['exchange_polygon0'], # P['exchange_polygon1'], # P['enquiry_polygon0'], # P['enquiry_polygon1']], # figname='culvert_polygon_output') self.openings = [] self.openings.append(Inflow(domain, polygon=P['exchange_polygon0'])) self.openings.append(Inflow(domain, polygon=P['exchange_polygon1'])) # Assume two openings for now: Referred to as 0 and 1 assert len(self.openings) == 2 # Store basic geometry self.end_points = [end_point0, end_point1] self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']] self.vector = P['vector'] self.distance = P['length'] self.verbose = verbose self.width = width self.height = height self.last_time = 0.0 # Create objects to update momentum (a bit crude at this stage) self.xmom_forcing0 = General_forcing(domain, 'xmomentum', polygon=P['exchange_polygon0']) self.xmom_forcing1 = General_forcing(domain, 'xmomentum', polygon=P['exchange_polygon1']) self.ymom_forcing0 = General_forcing(domain, 'ymomentum', polygon=P['exchange_polygon0']) self.ymom_forcing1 = General_forcing(domain, 'ymomentum', polygon=P['exchange_polygon1']) # Print something print 'Culvert Effective Length =', self.distance print 'Culvert Slope is Delta Z / Dist' print 'Culvert Direction is ', self.vector print 'Point 1m Up Stream is X,Y =' print 'Point 1m Down Stream is X,Y =' def __call__(self, domain): from anuga.utilities.numerical_tools import mean from anuga.utilities.polygon import inside_polygon from anuga.config import g, epsilon from Numeric import take, sqrt # Get average water depths at each opening time = domain.get_time() print 'Time = %f, delta_t = %f' %(time, time-self.last_time) openings = self.openings for i, opening in enumerate(openings): # Quantities corresponding to fluid exchange field for this opening stage_o = domain.quantities['stage'].get_values(location='centroids', indices=opening.indices) elevation_o = domain.quantities['elevation'].get_values(location='centroids', indices=opening.indices) xmomentum_o = domain.quantities['xmomentum'].get_values(location='centroids', indices=opening.indices) ymomentum_o = domain.quantities['xmomentum'].get_values(location='centroids', indices=opening.indices) # Compute mean velocity in the exchange area in front of the culvert (taking zero depths into account) depth_o = stage_o - elevation_o ux_o = xmomentum_o/(depth_o+epsilon) uy_o = ymomentum_o/(depth_o+epsilon) v_o = mean(sqrt(ux_o**2+uy_o**2)) d_o = mean(depth_o) w_o = mean(stage_o) z_o = mean(elevation_o) self.mean_xmomentum_o = mean(xmomentum_o) self.mean_ymomentum_o = mean(ymomentum_o) # Indices corresponding to energy enquiry field for this opening coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y) idx = inside_polygon(coordinates, self.enquiry_polygons[i]) if self.verbose: pass #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\ # %(i, len(idx), self.enquiry_polygons[i]) # Get average model values for points in # enquiry polygon for this opening dq = domain.quantities stage = dq['stage'].get_values(location='centroids', indices=idx) xmomentum = dq['xmomentum'].get_values(location='centroids', indices=idx) ymomentum = dq['ymomentum'].get_values(location='centroids', indices=idx) elevation = dq['elevation'].get_values(location='centroids', indices=idx) depth = stage - elevation # Compute mean velocity in the area in front of the culvert (taking zero depths into account) ux = xmomentum/(depth+epsilon) uy = ymomentum/(depth+epsilon) v = mean(sqrt(ux**2+uy**2)) d = mean(depth) w = mean(stage) z = mean(elevation) self.mean_xmomentum = mean(xmomentum) self.mean_ymomentum = mean(ymomentum) # Ratio of depth to culvert height ratio = d/(2*self.height) if ratio > 1.0: # Assume culvert is running full & ratio = 1.0 # under pressure. Note this is usually ~ 1.35 opening.ratio = ratio # Average measure of total energy (D + V^2/2g) in enquiry field in front of this opening E = d + 0.5*v**2/g # RUDY - please check this opening.energy = E # Store current average stage and depth at enquiry field with each opening object opening.depth = d opening.stage = w opening.elevation = z # Now we are done calculating energy etc for each opening. # At this point we can work on the transfer functions etc. # Handy values (all calculated at enquiry polygon - # if you need them at exchange polygons we can easily do that.) delta_h = openings[1].stage - openings[0].stage #delta_h = openings[1].depth - openings[0].depth delta_z = openings[1].elevation - openings[0].elevation # Ideas..... if openings[0].depth > 0 and openings[1].depth > 0: flow_rate = delta_h else: flow_rate = 0.0 if openings[0].depth > self.height: # This could be usefull mayby #print 'Inlet has been overflowed' pass if openings[1].depth > self.height: # This could be usefull mayby #print 'Outlet has been overflowed' pass if delta_h > 0: # Water Level U/S is higher than DS flow_rate_orifice = 0.6*openings[0].area*(2*g*delta_h)**0.5 # Orifice Eqn Q= cA(2gh)^0.5 flow_rate_weir = 1.69*(openings[0].area)*openings[0].depth**1.5 # WEIR Eqn Q= CLH^1.5 # Does Weir or Orifice Control Flow Rate ? if flow_rate_weir > flow_rate_orifice: flow_rate_control = flow_rate_orifice else: flow_rate_control = flow_rate_weir elif delta_h < 0: # Water Level D/S is higher than US # That is reverse flow in culvert flow_rate_orifice = 0.6*openings[0].area*(2*g*-delta_h)**0.5 # Orifice Eqn Q= cA(2gh)^0.5 flow_rate_weir = 1.69*openings[0].area*openings[0].depth**1.5 # WEIR Eqn Q= CLH^1.5 # Does Weir or Orifice Control Flow Rate ? if flow_rate_weir > flow_rate_orifice: flow_rate_control = flow_rate_orifice else: flow_rate_control = flow_rate_weir if openings[0].depth < epsilon: # Do nothing because No water over Opening That is Set Flow to Zero!! self.openings[0].rate = 0 self.openings[1].rate = 0 elif openings[0].depth > 0: # Only calculate fllow if there is some depth over the inlet #if delta_h > 0: if 1: # FIXME (OLE): We never get in here. Why? #print 'We got water at inlet and dh > 0' # Flow will go from opening 0 to opening 1 #- that is abstract from [0] and add to [1] self.openings[0].rate = -flow_rate_control self.openings[1].rate = flow_rate_control # Add jet in the form of absolute momentum to opening 1 #speed = 20.0 # m/s # This Should be Based on the VELOCITY in the Culvert # Based on the Flow Depth in the Culvert for part full # or Flow Jet of the Pressurised Culvert if FUll # If Part Full Flow Calculate Part Full Depth # Based on Depth Calculate Area... the Vel = flow_rate_control / Area # Need to Break Velocity into X & Y Components #self.openings[1].set_quantity_values(delta_h*speed, 'xmomentum') # Previous Calculated Depth/ Culvert Heigh Ratio Use it to Determine Velocity !!!! FIX LAter #if (ratio*self.area)== 0: # Don't Really need this already established water depth here # outlet_vel=0.0 #else: # FIXME (Ole): Shouldn't this be openings[1]?? outlet_vel =(flow_rate_control/(openings[0].ratio*openings[0].area)) outlet_dep = 2.0*openings[0].depth*openings[0].ratio #????? outlet_mom = outlet_vel*outlet_dep # Eventually will Need Momentum in X & Y Components based on the orientation of # the culvert from X0,Y0, X1, Y1 from Create Polygon Routine # YES - use self.vector which is a unit vector in the direction of the culvert. outlet_mom_x, outlet_mom_y = self.vector * outlet_mom #print 'Directional momentum', outlet_mom_x, outlet_mom_y # Update the momentum forcing terms with directional momentum at the outlet delta_t = time - self.last_time if delta_t > 0.0: xmomentum_rate = outlet_mom_x - self.mean_xmomentum if xmomentum_rate > 0: xmomentum_rate /= delta_t else: xmomentum_rate = 0.0 ymomentum_rate = outlet_mom_y - self.mean_ymomentum if ymomentum_rate > 0: ymomentum_rate /= delta_t else: ymomentum_rate = 0.0 else: xmomentum_rate = ymomentum_rate = 0.0 # Set momentum rates for outlet jet self.xmom_forcing1.rate = xmomentum_rate self.ymom_forcing1.rate = ymomentum_rate # Remember this value for next step (IMPORTANT) self.xmom_forcing1.value = outlet_mom_x self.ymom_forcing1.value = outlet_mom_y if int(time*100) % 100 == 0: s = 'T=%.3f, Culvert Discharge = %.3f Culv. Vel. %0.3f'\ %(time, flow_rate_control, outlet_vel) s += ' Depth= %0.3f\n Outlet Momentum = (%0.3f, %0.3f)\n'\ %(outlet_dep, outlet_mom_x, outlet_mom_y) s += ' Avg Momentum at opening = (%0.3f, %0.3f)\n'\ %(self.mean_xmomentum_o, self.mean_ymomentum_o) s += ' Avg Momentum in enquiry = (%0.3f, %0.3f)\n'\ %(self.mean_xmomentum, self.mean_ymomentum) s += ' Momentum rate: (%.4f, %.4f)'\ %(xmomentum_rate, ymomentum_rate) fid.write(s) print s else: # Opening 1 has the greatest depth. Therefore Reverse the Flow !!! # Flow will go from opening 1 to opening 0, That is Abstract from [1] and add to [0] self.openings[0].rate = flow_rate_control # Else it will be Orifice Flow (Going US) self.openings[1].rate = -flow_rate_control # Second Else.... if water at outlet before at inlet !!! self.openings[1].rate = 10 # Execute flow term for each opening # This is where Inflow objects are evaluated and update the domain for opening in self.openings: opening(domain) # Execute momentum terms # This is where Inflow objects are evaluated and update the domain self.xmom_forcing0(domain) self.ymom_forcing0(domain) self.xmom_forcing1(domain) self.ymom_forcing1(domain) # Print out flows every 1 seconds if int(time*100) % 100 == 0: s = 'Time %.3f\n' %time s += ' Opening[0]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\ %(openings[0].depth, openings[0].area, openings[0].energy, openings[0].ratio) s += ' Opening[1]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\ %(openings[1].depth, openings[1].area, openings[1].energy, openings[1].ratio) s += ' Distance=%.2f, W=%.3f, O=%.3f, C=%.3f\n'\ %(self.distance, flow_rate_weir, flow_rate_orifice, flow_rate_control) print s fid.write(s) # Store value of time self.last_time = time #------------------------------------------------------------------------------ # Setup culvert inlets and outlets in current topography #------------------------------------------------------------------------------ # Define culvert inlet and outlets culvert = Culvert_flow(domain, end_point0=[9.0, 2.5], end_point1=[13.0, 2.5], width=1.00, verbose=True) domain.forcing_terms.append(culvert) #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ #Bi = Dirichlet_boundary([0.5, 0.0, 0.0]) # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!! Bi = Dirichlet_boundary([0.0, 0.0, 0.0]) # Inflow based on Flow Depth and Approaching Momentum !!! Br = Reflective_boundary(domain) # Solid reflective wall Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br}) #------------------------------------------------------------------------------ # Setup Application of specialised forcing terms #------------------------------------------------------------------------------ # This is the new element implemented by Ole to allow direct input of Inflow in m^3/s fixed_flow = Inflow(domain, rate=6.00, center=(2.1, 2.1), radius=1.261566) # Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s # flow=file_function('Q/QPMF_Rot_Sub13.tms')) # Read Time Series in from File # flow=lambda t: min(0.01*t, 0.01942)) # Time Varying Function Tap turning up domain.forcing_terms.append(fixed_flow) #------------------------------------------------------------------------------ # Evolve system through time #------------------------------------------------------------------------------ for t in domain.evolve(yieldstep = 0.1, finaltime = 10): pass #if int(domain.time*100) % 100 == 0: # domain.write_time() #if domain.get_time() >= 4 and tap.rate != 0.0: # print 'Turning tap off' # tap.rate = 0.0 #if domain.get_time() >= 3 and sink.rate < 0.0: # print 'Turning drain on' # sink.rate = -0.8 # Close fid.close() #------------------------------------------------------------------------------ # Query output #------------------------------------------------------------------------------ from anuga.shallow_water.data_manager import get_flow_through_cross_section swwfilename = domain.get_name()+'.sww' # Output name from script print swwfilename polyline = [[17., 0.], [17., 5.]] time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True) from pylab import ion, plot ion() plot(time, Q) raw_input('done')