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

Last change on this file since 5436 was 5436, checked in by ole, 16 years ago

Moved culvert routines into own area

File size: 8.8 KB
Line 
1"""Collection of culvert routines for use with Culvert_flow in culvert_class
2
3Usage:
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
31from math import pi, sqrt, sin, cos
32
33
34def 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 >= 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.421*g**0.5*diameter**0.87*E**1.63    # Inlet Ctrl Inlet Unsubmerged
80            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**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
132        else:
133            # Box culvert (rectangle or square)
134
135            # Calculate flows for inlet control
136            height = culvert.height
137            width = culvert.width           
138           
139            Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
140            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
141
142            s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
143            log_to_file(log_filename, s)
144
145            case = ''
146            if Q_inlet_unsubmerged < Q_inlet_submerged:
147                Q = Q_inlet_unsubmerged
148                flow_area = width*inlet.depth
149                outlet_culvert_depth = inlet.depth
150                perimeter=(width+2.0*inlet.depth)               
151                case = 'Inlet unsubmerged'
152            else:   
153                Q = Q_inlet_submerged
154                flow_area = width*height
155                outlet_culvert_depth = height
156                perimeter=2.0*(width+height)               
157                case = 'Inlet submerged'                   
158
159            if delta_Et < E:
160                # Calculate flows for outlet control
161                # Determine the depth at the outlet relative to the depth of flow in the Culvert
162
163                if outlet.depth > height:        # The Outlet is Submerged
164                    outlet_culvert_depth=height
165                    flow_area=width*height       # Cross sectional area of flow in the culvert
166                    perimeter=2.0*(width+height)
167                    case = 'Outlet submerged'
168                elif outlet.depth==0.0: 
169                    outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
170                    flow_area=width*inlet.depth
171                    perimeter=(width+2.0*inlet.depth)
172                    case = 'Outlet depth is zero'                       
173                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
174                    outlet_culvert_depth=outlet.depth
175                    flow_area=width*outlet.depth
176                    perimeter=(width+2.0*outlet.depth)
177                    case = 'Outlet is open channel flow'
178                   
179
180
181
182        # Common code for rectangular and circular culvert types   
183        hyd_rad = flow_area/perimeter
184        s = 'hydraulic radius at outlet = %f' %hyd_rad
185        log_to_file(log_filename, s)
186
187        # Outlet control velocity using tail water
188        culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) 
189        Q_outlet_tailwater = flow_area * culvert_velocity
190
191        s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
192        log_to_file(log_filename, s)
193        Q = min(Q, Q_outlet_tailwater)
194
195        log_to_file(log_filename, 'Case: "%s"' %case)
196   
197        flow_rate_control=Q
198
199        s = 'Flow Rate Control = %f' %flow_rate_control
200        log_to_file(log_filename, s)
201
202        inlet.rate = -flow_rate_control
203        outlet.rate = flow_rate_control               
204           
205        culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3))
206        s = 'Froude in Culvert = %f' %culv_froude
207        log_to_file(log_filename, s)
208
209        # Determine momentum at the outlet
210        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
211    else: #inlet.depth < 0.01:
212        Q = barrel_velocity = outlet_culvert_depth = 0.0
213       
214    return Q, barrel_velocity, outlet_culvert_depth
215
216
Note: See TracBrowser for help on using the repository browser.