# source:anuga_core/source/anuga/culvert_flows/culvert_routines.py@6299

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

Cleaned up culvert routine and made all input parameters explicit.

File size: 9.8 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_total_energy > inlet_specific_energy
11# Outlet control: Delta_total_energy < inlet_specific_energy
12# where total energy is (w + 0.5*v^2/g) and
13# specific energy is (h + 0.5*v^2/g)
14
15
16from math import pi, sqrt, sin, cos
17
18
19def boyd_generalised_culvert_model(inlet_depth,
20                                   outlet_depth,
21                                   inlet_specific_energy,
22                                   delta_total_energy,
23                                   g,
24                                   culvert_length=0.0,
25                                   culvert_width=0.0,
26                                   culvert_height=0.0,
27                                   culvert_type='box',
28                                   manning=0.0,
29                                   sum_loss=0.0,
30                                   log_filename=None):
31
32    """Boyd's generalisation of the US department of transportation culvert
33    model
34
35    The quantity of flow passing through a culvert is controlled by many factors
36    It could be that the culvert is controlled by the inlet only, with it
37    being unsubmerged this is effectively equivalent to the weir Equation
38    Else the culvert could be controlled by the inlet, with it being
39    submerged, this is effectively the Orifice Equation
40    Else it may be controlled by down stream conditions where depending on
41    the down stream depth, the momentum in the culvert etc. flow is controlled
42    """
43
44    from anuga.utilities.system_tools import log_to_file
45    from anuga.config import velocity_protection
46    from anuga.utilities.numerical_tools import safe_acos as acos
47
48
49    if inlet_depth > 0.01:
50        # Water has risen above inlet
51
52        if log_filename is not None:
53            s = 'Specific energy  = %f m' % inlet_specific_energy
54            log_to_file(log_filename, s)
55
56        msg = 'Specific energy at inlet is negative'
57        assert inlet_specific_energy >= 0.0, msg
58
59
60        if culvert_type == 'circle':
61            # Round culvert (use width as diameter)
62            diameter = culvert_width
63
64            # Calculate flows for inlet control
65            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
66            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
67
68            if log_filename is not None:
69                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
70                log_to_file(log_filename, s)
71
72
73            # FIXME(Ole): Are these functions really for inlet control?
74            if Q_inlet_unsubmerged < Q_inlet_submerged:
75                Q = Q_inlet_unsubmerged
76                alpha = acos(1 - inlet_depth/diameter)
77                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
78                outlet_culvert_depth = inlet_depth
79                case = 'Inlet unsubmerged'
80            else:
81                Q = Q_inlet_submerged
82                flow_area = (diameter/2)**2 * pi
83                outlet_culvert_depth = diameter
84                case = 'Inlet submerged'
85
86
87
88            if delta_total_energy < inlet_specific_energy:
89                # Calculate flows for outlet control
90
91                # Determine the depth at the outlet relative to the depth of flow in the Culvert
92                if outlet_depth > diameter:        # The Outlet is Submerged
93                    outlet_culvert_depth=diameter
94                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
95                    perimeter = diameter * pi
96                    case = 'Outlet submerged'
97                elif outlet_depth==0.0:
98                    outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth
99                    alpha = acos(1 - inlet_depth/diameter)
100                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
101                    perimeter = alpha*diameter
102
103                    case = 'Outlet depth is zero'
104                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
105                    outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
106                    alpha = acos(1 - outlet_depth/diameter)
107                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
108                    perimeter = alpha*diameter
109                    case = 'Outlet is open channel flow'
110
112
113                if log_filename is not None:
115                    log_to_file(log_filename, s)
116
117                # Outlet control velocity using tail water
119                Q_outlet_tailwater = flow_area * culvert_velocity
120
121                if log_filename is not None:
122                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
123                    log_to_file(log_filename, s)
124
125                Q = min(Q, Q_outlet_tailwater)
126            else:
127                pass
128                #FIXME(Ole): What about inlet control?
129
130
131        else:
132            # Box culvert (rectangle or square)
133
134            # Calculate flows for inlet control
135            height = culvert_height
136            width = culvert_width
137
138            Q_inlet_unsubmerged = 0.540*g**0.5*width*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
139            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
140
141            if log_filename is not None:
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
146            # FIXME(Ole): Are these functions really for inlet control?
147            if Q_inlet_unsubmerged < Q_inlet_submerged:
148                Q = Q_inlet_unsubmerged
149                flow_area = width*inlet_depth
150                outlet_culvert_depth = inlet_depth
151                case = 'Inlet unsubmerged'
152            else:
153                Q = Q_inlet_submerged
154                flow_area = width*height
155                outlet_culvert_depth = height
156                case = 'Inlet submerged'
157
158            if delta_total_energy < inlet_specific_energy:
159                # Calculate flows for outlet control
160
161                # Determine the depth at the outlet relative to the depth of flow in the Culvert
162                if outlet_depth > height:        # The Outlet is Submerged
163                    outlet_culvert_depth=height
164                    flow_area=width*height       # Cross sectional area of flow in the culvert
165                    perimeter=2.0*(width+height)
166                    case = 'Outlet submerged'
167                elif outlet_depth==0.0:
168                    outlet_culvert_depth=inlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
169                    flow_area=width*inlet_depth
170                    perimeter=(width+2.0*inlet_depth)
171                    case = 'Outlet depth is zero'
172                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
173                    outlet_culvert_depth=outlet_depth
174                    flow_area=width*outlet_depth
175                    perimeter=(width+2.0*outlet_depth)
176                    case = 'Outlet is open channel flow'
177
179
180                if log_filename is not None:
182                    log_to_file(log_filename, s)
183
184                # Outlet control velocity using tail water
186                Q_outlet_tailwater = flow_area * culvert_velocity
187
188                if log_filename is not None:
189                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
190                    log_to_file(log_filename, s)
191                Q = min(Q, Q_outlet_tailwater)
192            else:
193                pass
194                #FIXME(Ole): What about inlet control?
195
196
197        # Common code for circle and square geometries
198        if log_filename is not None:
199            log_to_file(log_filename, 'Case: "%s"' % case)
200
201        if log_filename is not None:
202            s = 'Flow Rate Control = %f' % Q
203            log_to_file(log_filename, s)
204
205
206        culv_froude=sqrt(Q**2*width/(g*flow_area**3))
207
208        if log_filename is not None:
209            s = 'Froude in Culvert = %f' % culv_froude
210            log_to_file(log_filename, s)
211
212        # Determine momentum at the outlet
213        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
214
215
216    else: # inlet_depth < 0.01:
217        Q = barrel_velocity = outlet_culvert_depth = 0.0
218
219    return Q, barrel_velocity, outlet_culvert_depth
220
221
Note: See TracBrowser for help on using the repository browser.