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

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

Added culvert test info from Petar. One tests doesn't pass, but values are in the ballpark.

Also found error in calculations. Changed /2*g to /2/g


File size: 10.3 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                                   max_velocity=10.0,
31                                   log_filename=None):
32
33    """Boyd's generalisation of the US department of transportation culvert
34    model
35     
36    The quantity of flow passing through a culvert is controlled by many factors
37    It could be that the culvert is controlled by the inlet only, with it
38    being unsubmerged this is effectively equivalent to the weir Equation
39    Else the culvert could be controlled by the inlet, with it being
40    submerged, this is effectively the Orifice Equation
41    Else it may be controlled by down stream conditions where depending on
42    the down stream depth, the momentum in the culvert etc. flow is controlled
43    """
44
45    from anuga.utilities.system_tools import log_to_file
46    from anuga.config import velocity_protection
47    from anuga.utilities.numerical_tools import safe_acos as acos
48
49
50    if inlet_depth > 0.1: #this value was 0.01:
51        # Water has risen above inlet
52       
53        if log_filename is not None:                           
54            s = 'Specific energy  = %f m' % inlet_specific_energy
55            log_to_file(log_filename, s)
56       
57        msg = 'Specific energy at inlet is negative'
58        assert inlet_specific_energy >= 0.0, msg
59                     
60
61        if culvert_type == 'circle':
62            # Round culvert (use width as diameter)
63            diameter = culvert_width           
64
65            # Calculate flows for inlet control           
66            Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*inlet_specific_energy**1.63 # Inlet Ctrl Inlet Unsubmerged
67            Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*inlet_specific_energy**0.63   # Inlet Ctrl Inlet Submerged
68
69            if log_filename is not None:                                       
70                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
71                log_to_file(log_filename, s)
72
73
74            # FIXME(Ole): Are these functions really for inlet control?                   
75            if Q_inlet_unsubmerged < Q_inlet_submerged:
76                Q = Q_inlet_unsubmerged
77                alpha = acos(1 - inlet_depth/diameter)
78                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
79                outlet_culvert_depth = inlet_depth
80                case = 'Inlet unsubmerged'
81            else:   
82                Q = Q_inlet_submerged
83                flow_area = (diameter/2)**2 * pi
84                outlet_culvert_depth = diameter
85                case = 'Inlet submerged'                   
86
87
88
89            if delta_total_energy < inlet_specific_energy:
90                # Calculate flows for outlet control
91               
92                # Determine the depth at the outlet relative to the depth of flow in the Culvert
93                if outlet_depth > diameter:        # The Outlet is Submerged
94                    outlet_culvert_depth=diameter
95                    flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert
96                    perimeter = diameter * pi
97                    case = 'Outlet submerged'
98                elif outlet_depth==0.0: 
99                    outlet_culvert_depth=inlet_depth # For purpose of calculation assume the outlet depth = the inlet depth
100                    alpha = acos(1 - inlet_depth/diameter)
101                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
102                    perimeter = alpha*diameter
103
104                    case = 'Outlet depth is zero'                       
105                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
106                    outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
107                    alpha = acos(1 - outlet_depth/diameter)
108                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
109                    perimeter = alpha*diameter
110                    case = 'Outlet is open channel flow'
111
112                hyd_rad = flow_area/perimeter
113               
114                if log_filename is not None:
115                    s = 'hydraulic radius at outlet = %f' %hyd_rad
116                    log_to_file(log_filename, s)
117
118                # Outlet control velocity using tail water
119                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
120                Q_outlet_tailwater = flow_area * culvert_velocity
121               
122                if log_filename is not None:               
123                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
124                    log_to_file(log_filename, s)
125                   
126                Q = min(Q, Q_outlet_tailwater)
127            else:
128                pass
129                #FIXME(Ole): What about inlet control?
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*inlet_specific_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
140            Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*inlet_specific_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
141
142            if log_filename is not None:                           
143                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
144                log_to_file(log_filename, s)
145
146
147            # FIXME(Ole): Are these functions really for inlet control?   
148            if Q_inlet_unsubmerged < Q_inlet_submerged:
149                Q = Q_inlet_unsubmerged
150                flow_area = width*inlet_depth
151                outlet_culvert_depth = inlet_depth
152                case = 'Inlet unsubmerged'
153            else:   
154                Q = Q_inlet_submerged
155                flow_area = width*height
156                outlet_culvert_depth = height
157                case = 'Inlet submerged'                   
158
159            if delta_total_energy < inlet_specific_energy:
160                # Calculate flows for outlet control
161               
162                # Determine the depth at the outlet relative to the depth of flow in the Culvert
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                hyd_rad = flow_area/perimeter
180               
181                if log_filename is not None:                               
182                    s = 'hydraulic radius at outlet = %f' % hyd_rad
183                    log_to_file(log_filename, s)
184
185                # Outlet control velocity using tail water
186                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
187                Q_outlet_tailwater = flow_area * culvert_velocity
188
189                if log_filename is not None:                           
190                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
191                    log_to_file(log_filename, s)
192                Q = min(Q, Q_outlet_tailwater)
193            else:
194                pass
195                #FIXME(Ole): What about inlet control?
196
197               
198        # Common code for circle and square geometries
199        if log_filename is not None:
200            log_to_file(log_filename, 'Case: "%s"' % case)
201           
202        if log_filename is not None:                       
203            s = 'Flow Rate Control = %f' % Q
204            log_to_file(log_filename, s)
205
206           
207        culv_froude=sqrt(Q**2*width/(g*flow_area**3))
208       
209        if log_filename is not None:                           
210            s = 'Froude in Culvert = %f' % culv_froude
211            log_to_file(log_filename, s)
212
213        # Determine momentum at the outlet
214        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
215
216
217    else: # inlet_depth < 0.01:
218        Q = barrel_velocity = outlet_culvert_depth = 0.0
219
220    # Temporary flow limit
221    if barrel_velocity > max_velocity:
222        if log_filename is not None:                           
223            s = 'Barrel velocity was reduced from = %f m/s to %f m/s' % (barrel_velocity, max_velocity)
224            log_to_file(log_filename, s)
225       
226        barrel_velocity = max_velocity
227        Q = flow_area * barrel_velocity
228       
229       
230
231       
232    return Q, barrel_velocity, outlet_culvert_depth
233
234
Note: See TracBrowser for help on using the repository browser.