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

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

Back-merge from Numeric trunk to numpy branch.

File size: 10.3 KB
RevLine 
[5436]1"""Collection of culvert routines for use with Culvert_flow in culvert_class
2
3Usage:
4
5   
6
7"""
8
[6121]9#NOTE:
[6553]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)
[6121]14
15
[5436]16from math import pi, sqrt, sin, cos
17
18
[6553]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,
[6689]30                                   max_velocity=10.0,
[6553]31                                   log_filename=None):
[5436]32
[6553]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
[5436]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
[6533]49
[6689]50    if inlet_depth > 0.1: #this value was 0.01:
[6553]51        # Water has risen above inlet
52       
53        if log_filename is not None:                           
54            s = 'Specific energy  = %f m' % inlet_specific_energy
[6128]55            log_to_file(log_filename, s)
[5436]56       
[6553]57        msg = 'Specific energy at inlet is negative'
58        assert inlet_specific_energy >= 0.0, msg
[5436]59                     
60
[6553]61        if culvert_type == 'circle':
62            # Round culvert (use width as diameter)
63            diameter = culvert_width           
64
[5436]65            # Calculate flows for inlet control           
[6553]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
[5436]68
[6553]69            if log_filename is not None:                                       
70                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)
[6128]71                log_to_file(log_filename, s)
[5436]72
[6553]73
74            # FIXME(Ole): Are these functions really for inlet control?                   
[5436]75            if Q_inlet_unsubmerged < Q_inlet_submerged:
76                Q = Q_inlet_unsubmerged
[6553]77                alpha = acos(1 - inlet_depth/diameter)
[5436]78                flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
[6553]79                outlet_culvert_depth = inlet_depth
[5436]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
[6553]89            if delta_total_energy < inlet_specific_energy:
[5436]90                # Calculate flows for outlet control
[6553]91               
[5436]92                # Determine the depth at the outlet relative to the depth of flow in the Culvert
[6553]93                if outlet_depth > diameter:        # The Outlet is Submerged
[5436]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'
[6553]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)
[5436]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
[6553]106                    outlet_culvert_depth=outlet_depth   # For purpose of calculation assume the outlet depth = the inlet depth
107                    alpha = acos(1 - outlet_depth/diameter)
[5436]108                    flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha))
109                    perimeter = alpha*diameter
110                    case = 'Outlet is open channel flow'
111
[5437]112                hyd_rad = flow_area/perimeter
[6128]113               
[6553]114                if log_filename is not None:
[6128]115                    s = 'hydraulic radius at outlet = %f' %hyd_rad
116                    log_to_file(log_filename, s)
[5436]117
[5437]118                # Outlet control velocity using tail water
[6689]119                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
[5437]120                Q_outlet_tailwater = flow_area * culvert_velocity
[6128]121               
[6553]122                if log_filename is not None:               
[6128]123                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
124                    log_to_file(log_filename, s)
125                   
[5437]126                Q = min(Q, Q_outlet_tailwater)
[6128]127            else:
128                pass
129                #FIXME(Ole): What about inlet control?
[5437]130
131
[5436]132        else:
133            # Box culvert (rectangle or square)
134
135            # Calculate flows for inlet control
[6553]136            height = culvert_height
137            width = culvert_width           
[5436]138           
[6553]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
[5436]141
[6553]142            if log_filename is not None:                           
[6128]143                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
144                log_to_file(log_filename, s)
[5436]145
[6553]146
147            # FIXME(Ole): Are these functions really for inlet control?   
[5436]148            if Q_inlet_unsubmerged < Q_inlet_submerged:
149                Q = Q_inlet_unsubmerged
[6553]150                flow_area = width*inlet_depth
151                outlet_culvert_depth = inlet_depth
[5436]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
[6553]159            if delta_total_energy < inlet_specific_energy:
[5436]160                # Calculate flows for outlet control
[6553]161               
[5436]162                # Determine the depth at the outlet relative to the depth of flow in the Culvert
[6553]163                if outlet_depth > height:        # The Outlet is Submerged
[5436]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'
[6553]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)
[5436]172                    case = 'Outlet depth is zero'                       
173                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
[6553]174                    outlet_culvert_depth=outlet_depth
175                    flow_area=width*outlet_depth
176                    perimeter=(width+2.0*outlet_depth)
[5436]177                    case = 'Outlet is open channel flow'
178
[5437]179                hyd_rad = flow_area/perimeter
[6128]180               
[6553]181                if log_filename is not None:                               
182                    s = 'hydraulic radius at outlet = %f' % hyd_rad
[6128]183                    log_to_file(log_filename, s)
[5436]184
[5437]185                # Outlet control velocity using tail water
[6689]186                culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 
[5437]187                Q_outlet_tailwater = flow_area * culvert_velocity
[5436]188
[6553]189                if log_filename is not None:                           
190                    s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater
[6128]191                    log_to_file(log_filename, s)
[5437]192                Q = min(Q, Q_outlet_tailwater)
[6128]193            else:
194                pass
195                #FIXME(Ole): What about inlet control?
[5436]196
[6553]197               
[5437]198        # Common code for circle and square geometries
[6553]199        if log_filename is not None:
200            log_to_file(log_filename, 'Case: "%s"' % case)
[6128]201           
[6553]202        if log_filename is not None:                       
203            s = 'Flow Rate Control = %f' % Q
[6128]204            log_to_file(log_filename, s)
[5436]205
[6533]206           
[6553]207        culv_froude=sqrt(Q**2*width/(g*flow_area**3))
[6128]208       
[6553]209        if log_filename is not None:                           
210            s = 'Froude in Culvert = %f' % culv_froude
[6128]211            log_to_file(log_filename, s)
[5436]212
213        # Determine momentum at the outlet
214        barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
[5437]215
216
[6553]217    else: # inlet_depth < 0.01:
[5436]218        Q = barrel_velocity = outlet_culvert_depth = 0.0
[6689]219    # Temporary flow limit
220    if barrel_velocity > max_velocity:
221        if log_filename is not None:                           
222            s = 'Barrel velocity was reduced from = %f m/s to %f m/s' % (barrel_velocity, max_velocity)
223            log_to_file(log_filename, s)
[5436]224       
[6689]225        barrel_velocity = max_velocity
226        Q = flow_area * barrel_velocity
227       
228       
229
230       
[5436]231    return Q, barrel_velocity, outlet_culvert_depth
232
233
Note: See TracBrowser for help on using the repository browser.