source: anuga_work/development/culvert_flow/culvert_boyd_channel.py @ 5437

Last change on this file since 5437 was 5437, checked in by ole, 14 years ago

Finished culvert flow work with a good working solution

File size: 8.7 KB
Line 
1"""
2
3Ole Check Culvert Routine from Line 258
4
5Although it is Setup as a Culvert with the Opening presented vertically,
6for now the calculation of flow rate is assuming a horizontal hole in the
7ground (Fix this Later)
8
9MOST importantly 2 things...
101. How to use the Create Polygon Routine to enquire Depth ( or later energy)
11infront of the Culvert
12
13Done (Ole)
14
152. How to apply the Culvert velocity and thereby Momentum to the Outlet
16Ject presented at the Outlet
17
18
19
20Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction
21
22This Version CALCULATES the Culvert Velocity and Uses it to establish
23The Culvert Outlet Momentum
24
25The Aim is to define a Flow Transfer function that Simulates a Culvert
26by using the Total Available Energy to Drive the Culvert
27as Derived by determining the Difference in Total Energy
28between 2 Points, Just Up stream and Just Down Stream of the Culvert
29away from the influence of the Flow Abstraction etc..
30
31This example includes a Model Topography that shows a
32TYPICAL Headwall Configuration
33
34The aim is to change the Culvert Routine to Model more precisely the
35abstraction
36from a vertical face.
37
38The inflow must include the impact of Approach velocity.
39Similarly the Outflow has MOMENTUM Not just Up welling as in the
40Horizontal Style
41abstraction
42
43"""
44
45#------------------------------------------------------------------------------
46# Import necessary modules
47#------------------------------------------------------------------------------
48from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
49from anuga.shallow_water import Domain
50from anuga.shallow_water.shallow_water_domain import Reflective_boundary
51from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
52from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
53
54from anuga.culvert_flows.culvert_class import Culvert_flow
55from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model
56
57
58from anuga.utilities.polygon import plot_polygons
59
60
61from math import pi,sqrt
62
63#------------------------------------------------------------------------------
64# Setup computational domain
65#------------------------------------------------------------------------------
66
67length = 40.
68width = 5.
69
70#dx = dy = 1           # Resolution: Length of subdivisions on both axes
71#dx = dy = .5           # Resolution: Length of subdivisions on both axes
72dx = dy = .25           # Resolution: Length of subdivisions on both axes
73#dx = dy = .1           # Resolution: Length of subdivisions on both axes
74
75# OLE.... How do I refine the resolution... in the area where I have the Culvert Opening ???......
76#  Can I refine in a X & Y Range ???
77points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
78                                               len1=length, len2=width)
79domain = Domain(points, vertices, boundary)   
80domain.set_name('culv_dev_HW_Var_Mom')                 # Output name
81domain.set_default_order(2)
82domain.H0 = 0.01
83domain.tight_slope_limiters = True
84domain.set_minimum_storable_height(0.001)
85
86
87
88
89
90#------------------------------------------------------------------------------
91# Setup initial conditions
92#------------------------------------------------------------------------------
93
94# Define the topography (land scape)
95def topography(x, y):
96    """Set up a weir
97   
98    A culvert will connect either side of an Embankment with a Headwall type structure
99    The aim is for the Model to REALISTICALY model flow through the Culvert
100    """
101    # General Slope of Topography
102    z=-x/100
103    floorhgt = 5
104    embank_hgt=10
105    embank_upslope=embank_hgt/5
106    embank_dnslope=embank_hgt/2.5
107    # Add bits and Pieces to topography
108    N = len(x)
109    for i in range(N):
110
111        # Sloping Embankment Across Channel
112       
113        if 0.0 < x[i] < 7.51:
114            z[i]=z[i]+5.0
115        if 7.5 < x[i] < 10.1:
116            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
117               z[i]=z[i]+5.0
118            else:
119               z[i] +=  embank_upslope*(x[i] -5.0)    # Sloping Segment  U/S Face
120        if 10.0 < x[i] < 12.1:
121            if  2.0 < y[i]  < 3.0: # Cut Out Segment for Culvert (open Channel)
122                #z[i] +=  z[i]+5-(x[i]-10)*2              # Sloping Channel in Embankment
123                z[i] +=  embank_hgt                 # Flat Crest of Embankment
124            else:
125                z[i] +=  embank_hgt                 # Flat Crest of Embankment
126        if 12.0 < x[i] < 14.5:
127            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
128               z[i]=z[i]
129            else:
130               z[i] += embank_hgt-embank_dnslope*(x[i] -12.0)       # Sloping D/S Face
131     
132    return z
133
134
135domain.set_quantity('elevation', topography)  # Use function for elevation
136domain.set_quantity('friction', 0.01)         # Constant friction
137domain.set_quantity('stage',
138                    expression='elevation')   # Dry initial condition
139
140
141
142#------------------------------------------------------------------------------
143# Setup culvert inlets and outlets in current topography
144#------------------------------------------------------------------------------
145
146# Define culvert inlet and outlets
147# NEED TO ADD Mannings N as Fixed Value or Function
148# Energy Loss Coefficients as Fixed or Function
149# Also Set the Shape & Gap Factors fo rthe Enquiry PolyGons
150# ALSO Allow the Invert Level to be provided by the USER
151
152
153culvert1 = Culvert_flow(domain,
154                       label='Culvert No. 1',
155                       description='This culvert is a test unit 1.2m Wide by 0.75m High',   
156                       end_point0=[9.0, 2.5], 
157                       end_point1=[13.0, 2.5],
158                       width=1.20,height=0.75,
159                       culvert_routine=boyd_generalised_culvert_model,       
160                       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)
171                       
172domain.forcing_terms.append(culvert1)
173domain.forcing_terms.append(culvert2)
174
175
176#------------------------------------------------------------------------------
177# Setup boundary conditions
178#------------------------------------------------------------------------------
179#Bi = Dirichlet_boundary([0.5, 0.0, 0.0])          # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!!
180Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
181Br = Reflective_boundary(domain)              # Solid reflective wall
182Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
183
184domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
185
186#------------------------------------------------------------------------------
187# Setup Application of  specialised forcing terms
188#------------------------------------------------------------------------------
189
190# This is the new element implemented by Ole to allow direct input of Inflow in m^3/s
191fixed_flow = Inflow(domain,
192                    rate=20.00,
193                    center=(2.1, 2.1),
194                    radius=1.261566)   #   Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s
195
196#            flow=file_function('Q/QPMF_Rot_Sub13.tms'))        # Read Time Series in  from File
197#             flow=lambda t: min(0.01*t, 0.01942))                             # Time Varying Function   Tap turning up         
198
199domain.forcing_terms.append(fixed_flow)
200
201
202#------------------------------------------------------------------------------
203# Evolve system through time
204#------------------------------------------------------------------------------
205
206
207for t in domain.evolve(yieldstep = 0.1, finaltime = 20):
208    pass
209    #if int(domain.time*100) % 100 == 0:
210    #    domain.write_time()
211   
212    #if domain.get_time() >= 4 and tap.flow != 0.0:
213    #    print 'Turning tap off'
214    #    tap.flow = 0.0
215       
216    #if domain.get_time() >= 3 and sink.flow < 0.0:
217    #    print 'Turning drain on'
218    #    sink.flow = -0.8       
219    # Close
220
221
222
223#------------------------------------------------------------------------------
224# Query output
225#------------------------------------------------------------------------------
226
227from anuga.shallow_water.data_manager import get_flow_through_cross_section
228
229swwfilename = domain.get_name()+'.sww'  # Output name from script
230print swwfilename
231
232polyline = [[17., 0.], [17., 5.]]
233
234time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True)
235
236from pylab import ion, plot
237ion()
238plot(time, Q)
239raw_input('done')
Note: See TracBrowser for help on using the repository browser.