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

Last change on this file since 6126 was 6123, checked in by ole, 15 years ago

Implented general inlet/outlet control.
Works with rating curve but not Boyd.

File size: 9.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_Et, 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    log_filename = culvert.log_filename
65
66    manning = culvert.manning
67    sum_loss = culvert.sum_loss
68    length = culvert.length
69
70    if inlet.depth_trigger >= 0.01 and inlet.depth >= 0.01:
71        # Calculate driving energy
72        # FIXME(Ole): Should this be specific energy?
73        E = inlet.total_energy
74
75        s = 'Driving energy  = %f m' %E
76        log_to_file(log_filename, s)
77       
78        msg = 'Driving energy is negative'
79        assert E >= 0.0, msg
80                     
81       
82        # Water has risen above inlet
83        if culvert.culvert_type == 'circle':
84            # Round culvert
85
86            # Calculate flows for inlet control           
87            diameter = culvert.diameter
88
89            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*E**1.63    # Inlet Ctrl Inlet Unsubmerged
90            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**0.63      # Inlet Ctrl Inlet Submerged
91
92            s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
93            log_to_file(log_filename, s)
94
95            case = ''
96            if Q_inlet_unsubmerged < Q_inlet_submerged:
97                Q = Q_inlet_unsubmerged
98
99                alpha = acos(1 - inlet.depth/diameter)
100                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
101                outlet_culvert_depth = inlet.depth
102                width = diameter*sin(alpha)
103                #perimeter = alpha*diameter               
104                case = 'Inlet unsubmerged'
105            else:   
106                Q = Q_inlet_submerged
107                flow_area = (diameter/2)**2 * pi
108                outlet_culvert_depth = diameter
109                width = diameter
110                #perimeter = diameter
111                case = 'Inlet submerged'                   
112
113
114
115            if delta_Et < E:
116                # Calculate flows for outlet control
117                # Determine the depth at the outlet relative to the depth of flow in the Culvert
118
119                if outlet.depth > diameter:        # The Outlet is Submerged
120                    outlet_culvert_depth=diameter
121                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
122                    perimeter = diameter * pi
123                    width = diameter
124                    case = 'Outlet submerged'
125                elif outlet.depth==0.0: 
126                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
127                    alpha = acos(1 - inlet.depth/diameter)
128                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
129                    perimeter = alpha*diameter
130                    width = diameter*sin(alpha)                   
131
132                    case = 'Outlet depth is zero'                       
133                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
134                    outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
135                    alpha = acos(1 - outlet.depth/diameter)
136                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
137                    perimeter = alpha*diameter
138                    width = diameter*sin(alpha)                   
139                    case = 'Outlet is open channel flow'
140
141                hyd_rad = flow_area/perimeter
142                s = 'hydraulic radius at outlet = %f' %hyd_rad
143                log_to_file(log_filename, s)
144
145                # Outlet control velocity using tail water
146                culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
147                Q_outlet_tailwater = flow_area * culvert_velocity
148
149                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
150                log_to_file(log_filename, s)
151                Q = min(Q, Q_outlet_tailwater)
152                   
153
154
155        else:
156            # Box culvert (rectangle or square)
157
158            # Calculate flows for inlet control
159            height = culvert.height
160            width = culvert.width           
161           
162            Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
163            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
164
165            s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
166            log_to_file(log_filename, s)
167
168            case = ''
169            if Q_inlet_unsubmerged < Q_inlet_submerged:
170                Q = Q_inlet_unsubmerged
171                flow_area = width*inlet.depth
172                outlet_culvert_depth = inlet.depth
173                #perimeter=(width+2.0*inlet.depth)               
174                case = 'Inlet unsubmerged'
175            else:   
176                Q = Q_inlet_submerged
177                flow_area = width*height
178                outlet_culvert_depth = height
179                #perimeter=2.0*(width+height)               
180                case = 'Inlet submerged'                   
181
182            if delta_Et < E:
183                # Calculate flows for outlet control
184                # Determine the depth at the outlet relative to the depth of flow in the Culvert
185
186                if outlet.depth > height:        # The Outlet is Submerged
187                    outlet_culvert_depth=height
188                    flow_area=width*height       # Cross sectional area of flow in the culvert
189                    perimeter=2.0*(width+height)
190                    case = 'Outlet submerged'
191                elif outlet.depth==0.0: 
192                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
193                    flow_area=width*inlet.depth
194                    perimeter=(width+2.0*inlet.depth)
195                    case = 'Outlet depth is zero'                       
196                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
197                    outlet_culvert_depth=outlet.depth
198                    flow_area=width*outlet.depth
199                    perimeter=(width+2.0*outlet.depth)
200                    case = 'Outlet is open channel flow'
201
202                hyd_rad = flow_area/perimeter
203                s = 'hydraulic radius at outlet = %f' %hyd_rad
204                log_to_file(log_filename, s)
205
206                # Outlet control velocity using tail water
207                culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
208                Q_outlet_tailwater = flow_area * culvert_velocity
209
210                s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
211                log_to_file(log_filename, s)
212                Q = min(Q, Q_outlet_tailwater)
213
214
215        # Common code for circle and square geometries
216        log_to_file(log_filename, 'Case: "%s"' %case)
217        flow_rate_control=Q
218
219        s = 'Flow Rate Control = %f' %flow_rate_control
220        log_to_file(log_filename, s)
221
222        inlet.rate = -flow_rate_control
223        outlet.rate = flow_rate_control               
224           
225        culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
226        s = 'Froude in Culvert = %f' %culv_froude
227        log_to_file(log_filename, s)
228
229        # Determine momentum at the outlet
230        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
231
232
233    else: #inlet.depth < 0.01:
234        Q = barrel_velocity = outlet_culvert_depth = 0.0
235       
236    return Q, barrel_velocity, outlet_culvert_depth
237
238
Note: See TracBrowser for help on using the repository browser.