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

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

Made get_intersecting_segments cacheable and used that
with get_energy_through_cross_section and similar functions.

This helped the culvert routine along heaps, but I still believe
the interpolation could be made much more efficient.
File size: 9.4 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_trigger >= 0.01 and inlet.depth >= 0.01:
62        # Calculate driving energy
63        E = inlet.total_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                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/2*g)+(manning**2*length)/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.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
153            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**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/2*g)+(manning**2*length)/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**2*width/(g*flow_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.