source: branches/numpy/anuga/culvert_flows/culvert_routines.py @ 6533

Last change on this file since 6533 was 6533, checked in by rwilson, 15 years ago

Revert back to 6481, prior to auto-merge of trunk and numpy branch.

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