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

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

Multibarrel culvert functionality.

File size: 8.8 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                       number_of_barrels=2,
161                       verbose=True)
162
163culvert2 = Culvert_flow(domain,
164                       label='Culvert No. 2',
165                       description='This culvert is a circular test with d=1.2m',   
166                       end_point0=[9.0, 1.5], 
167                       end_point1=[30.0, 3.5],
168                       diameter=1.20,
169                       invert_level0=7,
170                       culvert_routine=boyd_generalised_culvert_model,       
171                       verbose=True)
172                       
173domain.forcing_terms.append(culvert1)
174domain.forcing_terms.append(culvert2)
175
176
177#------------------------------------------------------------------------------
178# Setup boundary conditions
179#------------------------------------------------------------------------------
180#Bi = Dirichlet_boundary([0.5, 0.0, 0.0])          # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!!
181Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
182Br = Reflective_boundary(domain)              # Solid reflective wall
183Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
184
185domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
186
187#------------------------------------------------------------------------------
188# Setup Application of  specialised forcing terms
189#------------------------------------------------------------------------------
190
191# This is the new element implemented by Ole to allow direct input of Inflow in m^3/s
192fixed_flow = Inflow(domain,
193                    rate=20.00,
194                    center=(2.1, 2.1),
195                    radius=1.261566)   #   Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s
196
197#            flow=file_function('Q/QPMF_Rot_Sub13.tms'))        # Read Time Series in  from File
198#             flow=lambda t: min(0.01*t, 0.01942))                             # Time Varying Function   Tap turning up         
199
200domain.forcing_terms.append(fixed_flow)
201
202
203#------------------------------------------------------------------------------
204# Evolve system through time
205#------------------------------------------------------------------------------
206
207
208for t in domain.evolve(yieldstep = 0.1, finaltime = 20):
209    pass
210    #if int(domain.time*100) % 100 == 0:
211    #    domain.write_time()
212   
213    #if domain.get_time() >= 4 and tap.flow != 0.0:
214    #    print 'Turning tap off'
215    #    tap.flow = 0.0
216       
217    #if domain.get_time() >= 3 and sink.flow < 0.0:
218    #    print 'Turning drain on'
219    #    sink.flow = -0.8       
220    # Close
221
222
223
224#------------------------------------------------------------------------------
225# Query output
226#------------------------------------------------------------------------------
227
228from anuga.shallow_water.data_manager import get_flow_through_cross_section
229
230swwfilename = domain.get_name()+'.sww'  # Output name from script
231print swwfilename
232
233polyline = [[17., 0.], [17., 5.]]
234
235time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True)
236
237from pylab import ion, plot
238ion()
239plot(time, Q)
240raw_input('done')
Note: See TracBrowser for help on using the repository browser.