source: branches/numpy/anuga/culvert_flows/culvert_routines.py @ 7035

Last change on this file since 7035 was 7035, checked in by rwilson, 15 years ago

Back-merge from Numeric trunk.

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