""" 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_class import Culvert_flow from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model from anuga.utilities.polygon import plot_polygons from math import pi,sqrt #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ 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('culv_dev_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.001) #------------------------------------------------------------------------------ # 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 floorhgt = 5 embank_hgt=10 embank_upslope=embank_hgt/5 embank_dnslope=embank_hgt/2.5 # Add bits and Pieces to topography N = len(x) for i in range(N): # Sloping Embankment Across Channel if 0.0 < x[i] < 7.51: z[i]=z[i]+5.0 if 7.5 < 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]+5.0 else: z[i] += embank_upslope*(x[i] -5.0) # Sloping Segment U/S Face if 10.0 < x[i] < 12.1: if 2.0 < y[i] < 3.0: # Cut Out Segment for Culvert (open Channel) #z[i] += z[i]+5-(x[i]-10)*2 # Sloping Channel in Embankment z[i] += embank_hgt # Flat Crest of Embankment else: z[i] += embank_hgt # 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] += embank_hgt-embank_dnslope*(x[i] -12.0) # Sloping D/S Face 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 culvert inlets and outlets in current topography #------------------------------------------------------------------------------ # Define culvert inlet and outlets # NEED TO ADD Mannings N as Fixed Value or Function # Energy Loss Coefficients as Fixed or Function # Also Set the Shape & Gap Factors fo rthe Enquiry PolyGons # ALSO Allow the Invert Level to be provided by the USER culvert1 = Culvert_flow(domain, label='Culvert No. 1', description='This culvert is a test unit 1.2m Wide by 0.75m High', end_point0=[9.0, 2.5], end_point1=[13.0, 2.5], width=1.20,height=0.75, culvert_routine=boyd_generalised_culvert_model, number_of_barrels=2, verbose=True) culvert2 = Culvert_flow(domain, label='Culvert No. 2', description='This culvert is a circular test with d=1.2m', end_point0=[9.0, 1.5], end_point1=[30.0, 3.5], diameter=1.20, invert_level0=7, culvert_routine=boyd_generalised_culvert_model, verbose=True) domain.forcing_terms.append(culvert1) domain.forcing_terms.append(culvert2) #------------------------------------------------------------------------------ # 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=20.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 = 20): pass #if int(domain.time*100) % 100 == 0: # domain.write_time() #if domain.get_time() >= 4 and tap.flow != 0.0: # print 'Turning tap off' # tap.flow = 0.0 #if domain.get_time() >= 3 and sink.flow < 0.0: # print 'Turning drain on' # sink.flow = -0.8 # 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')