Context Navigation

source: anuga_core/source/anuga/culvert_flows/culvert_routines.py @ 5453

Last change on this file since 5453 was 5437, checked in by ole, 16 years ago
Finished culvert flow work with a good working solution
File size: 9.5 KB

Line
1	"""Collection of culvert routines for use with Culvert_flow in culvert_class
2
3	Usage:
4
5
6
7	"""
8
9	# NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
10	#
11	# The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
12	# Flow Is Removed at a Rate of INFLOW
13	# Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
14	#
15	# This SHould be changed to a Vertical Opening Both BOX and Circular
16	# There will be several Culvert Routines such as:
17	# CULVERT_Boyd_Channel
18	# CULVERT_Orifice_and_Weir
19	# CULVERT_Simple_FLOOR
20	# CULVERT_Simple_WALL
21	# CULVERT_Eqn_Floor
22	# CULVERT_Eqn_Wall
23	# CULVERT_Tab_Floor
24	# CULVERT_Tab_Wall
25	# BRIDGES.....
26	# NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
27
28	# COULD USE EPA SWMM Model !!!!
29
30
31	from math import pi, sqrt, sin, cos
32
33
34	def boyd_generalised_culvert_model(culvert, inlet, outlet, delta_Et, g):
35
36	"""Boyd's generalisation of the US department of transportation culvert model
37	# == The quantity of flow passing through a culvert is controlled by many factors
38	# == 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
39	# == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
40	# == 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
41	"""
42
43	from anuga.utilities.system_tools import log_to_file
44	from anuga.config import velocity_protection
45	from anuga.utilities.numerical_tools import safe_acos as acos
46
47
48	Q_outlet_tailwater = 0.0
49	inlet.rate = 0.0
50	outlet.rate = 0.0
51	Q_inlet_unsubmerged = 0.0
52	Q_inlet_submerged = 0.0
53	Q_outlet_critical_depth = 0.0
54
55	log_filename = culvert.log_filename
56
57	manning = culvert.manning
58	sum_loss = culvert.sum_loss
59	length = culvert.length
60
61	if inlet.depth_trigger >= 0.01 and inlet.depth >= 0.01:
62	# Calculate driving energy
63	E = inlet.specific_energy
64
65	s = 'Driving energy = %f m' %E
66	log_to_file(log_filename, s)
67
68	msg = 'Driving energy is negative'
69	assert E >= 0.0, msg
70
71
72	# Water has risen above inlet
73	if culvert.culvert_type == 'circle':
74	# Round culvert
75
76	# Calculate flows for inlet control
77	diameter = culvert.diameter
78
79	Q_inlet_unsubmerged = 0.421g0.5diameter*0.87E**1.63 # Inlet Ctrl Inlet Unsubmerged
80	Q_inlet_submerged = 0.530g0.5diameter*1.87E**0.63 # Inlet Ctrl Inlet Submerged
81
82	s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
83	log_to_file(log_filename, s)
84
85	case = ''
86	if Q_inlet_unsubmerged < Q_inlet_submerged:
87	Q = Q_inlet_unsubmerged
88
89	alpha = acos(1 - inlet.depth/diameter)
90	flow_area = diameter*2 (alpha - sin(alpha)*cos(alpha))
91	outlet_culvert_depth = inlet.depth
92	width = diameter*sin(alpha)
93	#perimeter = alpha*diameter
94	case = 'Inlet unsubmerged'
95	else:
96	Q = Q_inlet_submerged
97	flow_area = (diameter/2)*2 pi
98	outlet_culvert_depth = diameter
99	width = diameter
100	#perimeter = diameter
101	case = 'Inlet submerged'
102
103
104
105	if delta_Et < E:
106	# Calculate flows for outlet control
107	# Determine the depth at the outlet relative to the depth of flow in the Culvert
108
109	if outlet.depth > diameter: # The Outlet is Submerged
110	outlet_culvert_depth=diameter
111	flow_area = (diameter/2)*2 pi # Cross sectional area of flow in the culvert
112	perimeter = diameter * pi
113	width = diameter
114	case = 'Outlet submerged'
115	elif outlet.depth==0.0:
116	outlet_culvert_depth=inlet.depth # For purpose of calculation assume the outlet depth = the inlet depth
117	alpha = acos(1 - inlet.depth/diameter)
118	flow_area = diameter*2 (alpha - sin(alpha)*cos(alpha))
119	perimeter = alpha*diameter
120	width = diameter*sin(alpha)
121
122	case = 'Outlet depth is zero'
123	else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
124	outlet_culvert_depth=outlet.depth # For purpose of calculation assume the outlet depth = the inlet depth
125	alpha = acos(1 - outlet.depth/diameter)
126	flow_area = diameter*2 (alpha - sin(alpha)*cos(alpha))
127	perimeter = alpha*diameter
128	width = diameter*sin(alpha)
129	case = 'Outlet is open channel flow'
130
131	hyd_rad = flow_area/perimeter
132	s = 'hydraulic radius at outlet = %f' %hyd_rad
133	log_to_file(log_filename, s)
134
135	# Outlet control velocity using tail water
136	culvert_velocity = sqrt(delta_Et/((sum_loss/2g)+(manning2length)/hyd_rad**1.33333))
137	Q_outlet_tailwater = flow_area * culvert_velocity
138
139	s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
140	log_to_file(log_filename, s)
141	Q = min(Q, Q_outlet_tailwater)
142
143
144
145	else:
146	# Box culvert (rectangle or square)
147
148	# Calculate flows for inlet control
149	height = culvert.height
150	width = culvert.width
151
152	Q_inlet_unsubmerged = 0.540g0.5widthE*1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
153	Q_inlet_submerged = 0.702g0.5widthheight0.89E**0.61 # Flow based on Inlet Ctrl Inlet Submerged
154
155	s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
156	log_to_file(log_filename, s)
157
158	case = ''
159	if Q_inlet_unsubmerged < Q_inlet_submerged:
160	Q = Q_inlet_unsubmerged
161	flow_area = width*inlet.depth
162	outlet_culvert_depth = inlet.depth
163	#perimeter=(width+2.0*inlet.depth)
164	case = 'Inlet unsubmerged'
165	else:
166	Q = Q_inlet_submerged
167	flow_area = width*height
168	outlet_culvert_depth = height
169	#perimeter=2.0*(width+height)
170	case = 'Inlet submerged'
171
172	if delta_Et < E:
173	# Calculate flows for outlet control
174	# Determine the depth at the outlet relative to the depth of flow in the Culvert
175
176	if outlet.depth > height: # The Outlet is Submerged
177	outlet_culvert_depth=height
178	flow_area=width*height # Cross sectional area of flow in the culvert
179	perimeter=2.0*(width+height)
180	case = 'Outlet submerged'
181	elif outlet.depth==0.0:
182	outlet_culvert_depth=inlet.depth # For purpose of calculation assume the outlet depth = the inlet depth
183	flow_area=width*inlet.depth
184	perimeter=(width+2.0*inlet.depth)
185	case = 'Outlet depth is zero'
186	else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
187	outlet_culvert_depth=outlet.depth
188	flow_area=width*outlet.depth
189	perimeter=(width+2.0*outlet.depth)
190	case = 'Outlet is open channel flow'
191
192	hyd_rad = flow_area/perimeter
193	s = 'hydraulic radius at outlet = %f' %hyd_rad
194	log_to_file(log_filename, s)
195
196	# Outlet control velocity using tail water
197	culvert_velocity = sqrt(delta_Et/((sum_loss/2g)+(manning2length)/hyd_rad**1.33333))
198	Q_outlet_tailwater = flow_area * culvert_velocity
199
200	s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
201	log_to_file(log_filename, s)
202	Q = min(Q, Q_outlet_tailwater)
203
204
205	# Common code for circle and square geometries
206	log_to_file(log_filename, 'Case: "%s"' %case)
207	flow_rate_control=Q
208
209	s = 'Flow Rate Control = %f' %flow_rate_control
210	log_to_file(log_filename, s)
211
212	inlet.rate = -flow_rate_control
213	outlet.rate = flow_rate_control
214
215	culv_froude=sqrt(flow_rate_control*2width/(gflow_area*3))
216	s = 'Froude in Culvert = %f' %culv_froude
217	log_to_file(log_filename, s)
218
219	# Determine momentum at the outlet
220	barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
221
222
223	else: #inlet.depth < 0.01:
224	Q = barrel_velocity = outlet_culvert_depth = 0.0
225
226	return Q, barrel_velocity, outlet_culvert_depth
227
228

Note: See TracBrowser for help on using the repository browser.

Download in other formats: