Changeset 7980
- Timestamp:
- Aug 30, 2010, 5:53:45 PM (14 years ago)
- Location:
- trunk/anuga_core/source/anuga/structures
- Files:
-
- 5 added
- 3 edited
- 2 moved
-
boyd_box_culvert.py (added)
-
boyd_box_routine.py (modified) (1 diff)
-
boyd_circle_routine.py (added)
-
culvert.py (moved) (moved from trunk/anuga_core/source/anuga/structures/box_culvert.py) (3 diffs)
-
culvert_operator.py (modified) (4 diffs)
-
culvert_routine.py (moved) (moved from trunk/anuga_core/source/anuga/structures/culvert_routines.py) (5 diffs)
-
inlet.py (modified) (11 diffs)
-
testing_open_slot_wide_bridge.py (added)
-
testing_wide_bridge.py (added)
-
testing_wide_bridge_old.py (added)
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/structures/boyd_box_routine.py
r7977 r7980 8 8 9 9 10 import culvert_routine 11 from anuga.config import velocity_protection 12 from anuga.utilities.numerical_tools import safe_acos as acos 10 13 11 def boyd_box(culvert): 14 from math import pi, sqrt, sin, cos 15 from anuga.config import g 16 17 18 class Boyd_box_routine(culvert_routine.Culvert_routine): 12 19 """Boyd's generalisation of the US department of transportation culvert methods 13 20 14 15 16 17 18 21 WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER 22 FULL PIPE OR CRITICAL DEPTH ONLY 23 For Supercritical flow this is UNDERESTIMATING the Outlet Velocity 24 The obtain the CORRECT velocity requires an iteration of Depth to Establish 25 the Normal Depth of flow in the pipe. 19 26 20 21 boyd_generalised_culvert_model_complex27 It is proposed to provide this in a seperate routine called 28 boyd_generalised_culvert_model_complex 22 29 23 The Boyd Method is based on methods described by the following:24 25 26 27 28 29 30 31 32 These documents provide around 60 charts for various configurations of culverts and inlets.30 The Boyd Method is based on methods described by the following: 31 1. 32 US Dept. Transportation Federal Highway Administration (1965) 33 Hydraulic Chart for Selection of Highway Culverts. 34 Hydraulic Engineering Circular No. 5 US Government Printing 35 2. 36 US Dept. Transportation Federal Highway Administration (1972) 37 Capacity charts for the Hydraulic design of highway culverts. 38 Hydraulic Engineering Circular No. 10 US Government Printing 39 These documents provide around 60 charts for various configurations of culverts and inlets. 33 40 34 35 36 37 38 39 40 41 41 Note these documents have been superceded by: 42 2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5), 43 Which combines culvert design information previously contained in Hydraulic Engineering Circulars 44 (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information. 45 HEC-5 provides 20 Charts 46 HEC-10 Provides an additional 36 Charts 47 HEC-13 Discusses the Design of improved more efficient inlets 48 HDS-5 Provides 60 sets of Charts 42 49 43 44 45 50 In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in 51 1987 published "Generalised Head Discharge Equations for Culverts". 52 These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations. 46 53 47 It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done. 48 The additional charts cover a range of culvert shapes and inlet configurations 49 50 """ 51 52 # Calculate flows for inflow control 53 54 Q_inflow_unsubmerged = 0.540*g**0.5*width*inflow_specific_energy**1.50 # Flow based on inflow Ctrl inflow Unsubmerged 55 Q_inflow_submerged = 0.702*g**0.5*width*height**0.89*inflow_specific_energy**0.61 # Flow based on inflow Ctrl inflow Submerged 56 57 if log_filename is not None: 58 s = 'Q_inflow_unsubmerged = %.6f, Q_inflow_submerged = %.6f' %(Q_inflow_unsubmerged, Q_inflow_submerged) 59 log_to_file(log_filename, s) 60 61 # FIXME(Ole): Are these functions really for inflow control? 62 if Q_inflow_unsubmerged < Q_inflow_submerged: 63 Q = Q_inflow_unsubmerged 64 dcrit = (Q**2/g/width**2)**0.333333 65 if dcrit > height: 66 dcrit = height 67 flow_area = width*dcrit 68 outflow_culvert_depth = dcrit 69 case = 'inflow unsubmerged Box Acts as Weir' 70 else: 71 Q = Q_inflow_submerged 72 flow_area = width*height 73 outflow_culvert_depth = height 74 case = 'inflow submerged Box Acts as Orifice' 75 76 dcrit = (Q**2/g/width**2)**0.333333 77 78 outflow_culvert_depth = dcrit 79 if outflow_culvert_depth > height: 80 outflow_culvert_depth = height # Once again the pipe is flowing full not partfull 81 flow_area = width*height # Cross sectional area of flow in the culvert 82 perimeter = 2*(width+height) 83 case = 'inflow CTRL outflow unsubmerged PIPE PART FULL' 84 else: 85 flow_area = width * outflow_culvert_depth 86 perimeter = width+2*outflow_culvert_depth 87 case = 'inflow CTRL Culvert is open channel flow we will for now assume critical depth' 88 89 if delta_total_energy < inflow_specific_energy: 90 # Calculate flows for outflow control 91 92 # Determine the depth at the outflow relative to the depth of flow in the Culvert 93 if outflow_depth > height: # The outflow is Submerged 94 outflow_culvert_depth=height 95 flow_area=width*height # Cross sectional area of flow in the culvert 96 perimeter=2.0*(width+height) 97 case = 'outflow submerged' 98 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 99 dcrit = (Q**2/g/width**2)**0.333333 100 outflow_culvert_depth=dcrit # For purpose of calculation assume the outflow depth = Critical Depth 101 if outflow_culvert_depth > height: 102 outflow_culvert_depth=height 103 flow_area=width*height 104 perimeter=2.0*(width+height) 105 case = 'outflow is Flowing Full' 106 else: 107 flow_area=width*outflow_culvert_depth 108 perimeter=(width+2.0*outflow_culvert_depth) 109 case = 'outflow is open channel flow' 110 111 hyd_rad = flow_area/perimeter 112 113 if log_filename is not None: 114 s = 'hydraulic radius at outflow = %f' % hyd_rad 115 log_to_file(log_filename, s) 116 117 # outflow control velocity using tail water 118 culvert_velocity = sqrt(delta_total_energy/((sum_loss/2/g)+(manning**2*culvert_length)/hyd_rad**1.33333)) 119 Q_outflow_tailwater = flow_area * culvert_velocity 120 121 if log_filename is not None: 122 s = 'Q_outflow_tailwater = %.6f' % Q_outflow_tailwater 123 log_to_file(log_filename, s) 124 Q = min(Q, Q_outflow_tailwater) 125 126 return Q 54 It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done. 55 The additional charts cover a range of culvert shapes and inlet configurations 127 56 128 57 129 if __name__ == "__main__": 58 """ 59 60 def __init__(self, culvert, manning=0.0): 61 62 culvert_routine.Culvert_routine.__init__(self, culvert, manning) 130 63 131 64 132 g=9.81133 culvert_slope=0.1 # Downward134 65 135 inlet_depth=2.0 136 outlet_depth=0.0 66 def __call__(self): 137 67 138 inlet_velocity=0.0, 139 outlet_velocity=0.0, 68 self.determine_inflow() 140 69 141 culvert_length=4.0 142 culvert_width=1.2 143 culvert_height=0.75 70 local_debug ='false' 71 72 if self.inflow.get_average_height() > 0.01: #this value was 0.01: 73 if local_debug =='true': 74 log.critical('Specific E & Deltat Tot E = %s, %s' 75 % (str(self.inflow.get_average_specific_energy()), 76 str(self.delta_total_energy))) 77 log.critical('culvert type = %s' % str(culvert_type)) 78 # Water has risen above inlet 144 79 145 culvert_type='box'146 manning=0.013147 sum_loss=0.080 if self.log_filename is not None: 81 s = 'Specific energy = %f m' % self.inflow.get_average_specific_energy() 82 log_to_file(self.log_filename, s) 148 83 149 inlet_specific_energy=inlet_depth #+0.5*v**2/g 150 z_in = 0.0 151 z_out = -culvert_length*culvert_slope/100 152 E_in = z_in+inlet_depth # + 153 E_out = z_out+outlet_depth # + 154 delta_total_energy = E_in-E_out 84 msg = 'Specific energy at inlet is negative' 85 assert self.inflow.get_average_specific_energy() >= 0.0, msg 155 86 156 Q = boyd_box(culvert_height, culvert_width, culvert_width, inlet_specific_energy) 87 height = self.culvert_height 88 width = self.culvert_width 89 flow_width = self.culvert_width 157 90 158 print 'Q ',Q 91 Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged 92 Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_average_specific_energy()**0.61 # Flow based on Inlet Ctrl Inlet Submerged 93 94 # FIXME(Ole): Are these functions really for inlet control? 95 if Q_inlet_unsubmerged < Q_inlet_submerged: 96 Q = Q_inlet_unsubmerged 97 dcrit = (Q**2/g/width**2)**0.333333 98 if dcrit > height: 99 dcrit = height 100 flow_area = width*dcrit 101 outlet_culvert_depth = dcrit 102 case = 'Inlet unsubmerged Box Acts as Weir' 103 else: 104 Q = Q_inlet_submerged 105 flow_area = width*height 106 outlet_culvert_depth = height 107 case = 'Inlet submerged Box Acts as Orifice' 108 109 dcrit = (Q**2/g/width**2)**0.333333 110 111 outlet_culvert_depth = dcrit 112 if outlet_culvert_depth > height: 113 outlet_culvert_depth = height # Once again the pipe is flowing full not partfull 114 flow_area = width*height # Cross sectional area of flow in the culvert 115 perimeter = 2*(width+height) 116 case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL' 117 else: 118 flow_area = width * outlet_culvert_depth 119 perimeter = width+2*outlet_culvert_depth 120 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth' 121 122 if self.delta_total_energy < self.inflow.get_average_specific_energy(): 123 # Calculate flows for outlet control 124 125 # Determine the depth at the outlet relative to the depth of flow in the Culvert 126 if self.outflow.get_average_height() > height: # The Outlet is Submerged 127 outlet_culvert_depth=height 128 flow_area=width*height # Cross sectional area of flow in the culvert 129 perimeter=2.0*(width+height) 130 case = 'Outlet submerged' 131 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity 132 dcrit = (Q**2/g/width**2)**0.333333 133 outlet_culvert_depth=dcrit # For purpose of calculation assume the outlet depth = Critical Depth 134 if outlet_culvert_depth > height: 135 outlet_culvert_depth=height 136 flow_area=width*height 137 perimeter=2.0*(width+height) 138 case = 'Outlet is Flowing Full' 139 else: 140 flow_area=width*outlet_culvert_depth 141 perimeter=(width+2.0*outlet_culvert_depth) 142 case = 'Outlet is open channel flow' 143 144 hyd_rad = flow_area/perimeter 145 146 if self.log_filename is not None: 147 s = 'hydraulic radius at outlet = %f' % hyd_rad 148 log_to_file(self.log_filename, s) 149 150 # Outlet control velocity using tail water 151 culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333)) 152 Q_outlet_tailwater = flow_area * culvert_velocity 153 154 if self.log_filename is not None: 155 s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater 156 log_to_file(self.log_filename, s) 157 Q = min(Q, Q_outlet_tailwater) 158 else: 159 pass 160 #FIXME(Ole): What about inlet control? 161 162 culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3)) 163 if local_debug =='true': 164 log.critical('FLOW AREA = %s' % str(flow_area)) 165 log.critical('PERIMETER = %s' % str(perimeter)) 166 log.critical('Q final = %s' % str(Q)) 167 log.critical('FROUDE = %s' % str(culv_froude)) 168 169 # Determine momentum at the outlet 170 barrel_velocity = Q/(flow_area + velocity_protection/flow_area) 171 172 # END CODE BLOCK for DEPTH > Required depth for CULVERT Flow 173 174 else: # self.inflow.get_average_height() < 0.01: 175 Q = barrel_velocity = outlet_culvert_depth = 0.0 176 177 # Temporary flow limit 178 if barrel_velocity > self.max_velocity: 179 barrel_velocity = self.max_velocity 180 Q = flow_area * barrel_velocity 181 182 183 184 185 186 return Q, barrel_velocity, outlet_culvert_depth 187 188 189 -
trunk/anuga_core/source/anuga/structures/culvert.py
r7979 r7980 6 6 import math 7 7 8 class Box_culvert:8 class Culvert: 9 9 """Culvert flow - transfer water from one rectangular box to another. 10 10 Sets up the geometry of problem … … 39 39 40 40 #FIXME (SR) Put this into a foe loop to deal with more inlets 41 41 42 self.inlets = [] 42 43 polygon0 = self.inlet_polygons[0] 43 inlet0_vector= self.culvert_vector44 self.inlets.append(inlet.Inlet(self.domain, polygon0 ))44 outward_vector0 = self.culvert_vector 45 self.inlets.append(inlet.Inlet(self.domain, polygon0, outward_vector0)) 45 46 46 47 polygon1 = self.inlet_polygons[1] 47 inlet1_vector = - self.culvert_vector 48 self.inlets.append(inlet.Inlet(self.domain, polygon1)) 49 48 outward_vector1 = - self.culvert_vector 49 self.inlets.append(inlet.Inlet(self.domain, polygon1, outward_vector1)) 50 51 52 def __call__(self): 53 msg = 'Method __call__ must be overridden by Culvert subclass' 54 raise Exception, msg 50 55 51 56 def __create_culvert_polygons(self): … … 73 78 # Short hands 74 79 w = 0.5*self.width*self.culvert_normal # Perpendicular vector of 1/2 width 75 h = self.height*self.culvert_vector # Vector of length=height in the80 h = 0.5*self.height*self.culvert_vector # Vector of length=height in the 76 81 # direction of the culvert 77 82 -
trunk/anuga_core/source/anuga/structures/culvert_operator.py
r7978 r7980 2 2 from anuga.config import g 3 3 import anuga.utilities.log as log 4 import box_culvert 5 import culvert_routines 4 5 from boyd_box_culvert import Boyd_box_culvert 6 6 7 7 class Culvert_operator: … … 34 34 self.height = height 35 35 36 self.culvert = box_culvert.Box_culvert(self.domain, end_points, self.width, self.height) 36 self.culvert = Boyd_box_culvert(self.domain, end_points, self.width, self.height) 37 print self.culvert 38 self.routine = self.culvert.routine 39 print self.routine 37 40 self.inlets = self.culvert.get_inlets() 38 41 … … 43 46 44 47 timestep = self.domain.get_timestep() 48 49 Q, barrel_speed, culvert_outlet_depth = self.routine() 45 50 46 from culvert_routines import Culvert_routines 47 culvert_routine = culvert_routines.Culvert_routines(self.culvert) 48 49 Q, barrel_velocity, culvert_outlet_depth = culvert_routine.boyd_circle() 51 inflow = self.routine.get_inflow() 52 outflow = self.routine.get_outflow() 50 53 51 transfer_water = Q*timestep54 outflow_direction = - outflow.outward_culvert_vector 52 55 53 inflow = culvert_routine.get_inflow() 54 outflow = culvert_routine.get_outflow() 56 outflow_momentum_flux = barrel_speed**2*culvert_outlet_depth*outflow_direction 55 57 56 inflow.set_heights(inflow.get_average_height() - transfer_water) 58 59 print Q, barrel_speed, culvert_outlet_depth, outflow_momentum_flux 60 61 #FIXME (SR) Check whether we need to mult/divide by inlet area 62 inflow_transfer = Q*timestep/inflow.get_area() 63 64 outflow_transfer = Q*timestep/outflow.get_area() 65 66 67 68 inflow.set_heights(inflow.get_average_height() - inflow_transfer) 69 57 70 inflow.set_xmoms(0.0) 58 71 inflow.set_ymoms(0.0) 59 72 60 outflow.set_heights(outflow.get_average_height() + transfer_water) 61 outflow.set_xmoms(0.0) 62 outflow.set_ymoms(0.0) 73 #u = outflow.get_xvelocities() 74 #v = outflow.get_yvelocities() 75 76 outflow.set_heights(outflow.get_average_height() + outflow_transfer) 77 #outflow.set_xmoms(outflow.get_xmoms() + timestep*outflow_momentum_flux[0] ) 78 #outflow.set_ymoms(outflow.get_ymoms() + timestep*outflow_momentum_flux[1] ) 63 79 64 80 def print_stats(self): … … 74 90 75 91 print 'inlet triangle indices and centres' 76 print inlet.triangle_indices [i]77 print self.domain.get_centroid_coordinates()[inlet.triangle_indices [i]]92 print inlet.triangle_indices 93 print self.domain.get_centroid_coordinates()[inlet.triangle_indices] 78 94 79 95 print 'polygon' -
trunk/anuga_core/source/anuga/structures/culvert_routine.py
r7979 r7980 5 5 from anuga.config import g 6 6 7 class Culvert_routine s:7 class Culvert_routine: 8 8 """Collection of culvert routines for use with Culvert_operator 9 9 … … 23 23 self.inlets = culvert.inlets 24 24 25 self.inflow = self.inlets[0] 26 self.outflow = self.inlets[1] 27 25 28 26 self.culvert_length = culvert.get_culvert_length() 29 27 self.culvert_width = culvert.get_culvert_width() … … 34 32 self.log_filename = None 35 33 36 # Determine flow direction based on total energy difference 37 self.delta_total_energy = self.inflow.get_average_total_energy() - self.outflow.get_average_total_energy() 38 39 if self.delta_total_energy < 0: 40 self.inflow = self.inlets[1] 41 self.outflow = self.inlets[0] 42 self.delta_total_energy = -self.delta_total_energy 34 self.determine_inflow() 43 35 44 36 #delta_z = self.self.inflow.get_average_elevation() - self.self.outflow.get_average_elevation() … … 53 45 # driving_head = self.inflow.get_average_specific_energy() 54 46 47 48 49 def determine_inflow(self): 50 # Determine flow direction based on total energy difference 51 52 self.delta_total_energy = self.inlets[0].get_average_total_energy() - self.inlets[1].get_average_total_energy() 53 54 self.inflow = self.inlets[0] 55 self.outflow = self.inlets[1] 55 56 57 58 if self.delta_total_energy < 0: 59 self.inflow = self.inlets[1] 60 self.outflow = self.inlets[0] 61 self.delta_total_energy = -self.delta_total_energy 62 63 56 64 def get_inflow(self): 57 65 … … 63 71 return self.outflow 64 72 65 66 def boyd_box(self):67 68 """Boyd's generalisation of the US department of transportation culvert methods69 70 WARNING THIS IS A SIMPLISTIC APPROACH and OUTLET VELOCITIES ARE LIMITED TO EITHER71 FULL PIPE OR CRITICAL DEPTH ONLY72 For Supercritical flow this is UNDERESTIMATING the Outlet Velocity73 The obtain the CORRECT velocity requires an iteration of Depth to Establish74 the Normal Depth of flow in the pipe.75 76 It is proposed to provide this in a seperate routine called77 boyd_generalised_culvert_model_complex78 79 The Boyd Method is based on methods described by the following:80 1.81 US Dept. Transportation Federal Highway Administration (1965)82 Hydraulic Chart for Selection of Highway Culverts.83 Hydraulic Engineering Circular No. 5 US Government Printing84 2.85 US Dept. Transportation Federal Highway Administration (1972)86 Capacity charts for the Hydraulic design of highway culverts.87 Hydraulic Engineering Circular No. 10 US Government Printing88 These documents provide around 60 charts for various configurations of culverts and inlets.89 90 Note these documents have been superceded by:91 2005 Hydraulic Design of Highway Culverts, Hydraulic Design Series No. 5 (HDS-5),92 Which combines culvert design information previously contained in Hydraulic Engineering Circulars93 (HEC) No. 5, No. 10, and No. 13 with hydrologic, storage routing, and special culvert design information.94 HEC-5 provides 20 Charts95 HEC-10 Provides an additional 36 Charts96 HEC-13 Discusses the Design of improved more efficient inlets97 HDS-5 Provides 60 sets of Charts98 99 In 1985 Professor Michael Boyd Published "Head-Discharge Relations for Culverts", and in100 1987 published "Generalised Head Discharge Equations for Culverts".101 These papers reviewed the previous work by the US DOT and provided a simplistic approach for 3 configurations.102 103 It may be possible to extend the same approach for additional charts in the original work, but to date this has not been done.104 The additional charts cover a range of culvert shapes and inlet configurations105 106 107 """108 109 local_debug ='false'110 if self.inflow.get_average_height() > 0.1: #this value was 0.01:111 if local_debug =='true':112 log.critical('Specific E & Deltat Tot E = %s, %s'113 % (str(self.inflow.get_average_specific_energy()),114 str(self.delta_total_energy)))115 log.critical('culvert type = %s' % str(culvert_type))116 # Water has risen above inlet117 118 if self.log_filename is not None:119 s = 'Specific energy = %f m' % self.inflow.get_average_specific_energy()120 log_to_file(self.log_filename, s)121 122 msg = 'Specific energy at inlet is negative'123 assert self.inflow.get_average_specific_energy() >= 0.0, msg124 125 height = self.culvert_height126 width = self.culvert_width127 flow_width = self.culvert_width128 129 Q_inlet_unsubmerged = 0.540*g**0.5*width*self.inflow.get_average_specific_energy()**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged130 Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*self.inflow.get_average_specific_energy()**0.61 # Flow based on Inlet Ctrl Inlet Submerged131 132 # FIXME(Ole): Are these functions really for inlet control?133 if Q_inlet_unsubmerged < Q_inlet_submerged:134 Q = Q_inlet_unsubmerged135 dcrit = (Q**2/g/width**2)**0.333333136 if dcrit > height:137 dcrit = height138 flow_area = width*dcrit139 outlet_culvert_depth = dcrit140 case = 'Inlet unsubmerged Box Acts as Weir'141 else:142 Q = Q_inlet_submerged143 flow_area = width*height144 outlet_culvert_depth = height145 case = 'Inlet submerged Box Acts as Orifice'146 147 dcrit = (Q**2/g/width**2)**0.333333148 149 outlet_culvert_depth = dcrit150 if outlet_culvert_depth > height:151 outlet_culvert_depth = height # Once again the pipe is flowing full not partfull152 flow_area = width*height # Cross sectional area of flow in the culvert153 perimeter = 2*(width+height)154 case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'155 else:156 flow_area = width * outlet_culvert_depth157 perimeter = width+2*outlet_culvert_depth158 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'159 160 if self.delta_total_energy < self.inflow.get_average_specific_energy():161 # Calculate flows for outlet control162 163 # Determine the depth at the outlet relative to the depth of flow in the Culvert164 if self.outflow.get_average_height() > height: # The Outlet is Submerged165 outlet_culvert_depth=height166 flow_area=width*height # Cross sectional area of flow in the culvert167 perimeter=2.0*(width+height)168 case = 'Outlet submerged'169 else: # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity170 dcrit = (Q**2/g/width**2)**0.333333171 outlet_culvert_depth=dcrit # For purpose of calculation assume the outlet depth = Critical Depth172 if outlet_culvert_depth > height:173 outlet_culvert_depth=height174 flow_area=width*height175 perimeter=2.0*(width+height)176 case = 'Outlet is Flowing Full'177 else:178 flow_area=width*outlet_culvert_depth179 perimeter=(width+2.0*outlet_culvert_depth)180 case = 'Outlet is open channel flow'181 182 hyd_rad = flow_area/perimeter183 184 if self.log_filename is not None:185 s = 'hydraulic radius at outlet = %f' % hyd_rad186 log_to_file(self.log_filename, s)187 188 # Outlet control velocity using tail water189 culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))190 Q_outlet_tailwater = flow_area * culvert_velocity191 192 if self.log_filename is not None:193 s = 'Q_outlet_tailwater = %.6f' % Q_outlet_tailwater194 log_to_file(self.log_filename, s)195 Q = min(Q, Q_outlet_tailwater)196 else:197 pass198 #FIXME(Ole): What about inlet control?199 200 culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))201 if local_debug =='true':202 log.critical('FLOW AREA = %s' % str(flow_area))203 log.critical('PERIMETER = %s' % str(perimeter))204 log.critical('Q final = %s' % str(Q))205 log.critical('FROUDE = %s' % str(culv_froude))206 207 # Determine momentum at the outlet208 barrel_velocity = Q/(flow_area + velocity_protection/flow_area)209 210 # END CODE BLOCK for DEPTH > Required depth for CULVERT Flow211 212 else: # self.inflow.get_average_height() < 0.01:213 Q = barrel_velocity = outlet_culvert_depth = 0.0214 215 # Temporary flow limit216 if barrel_velocity > self.max_velocity:217 barrel_velocity = self.max_velocity218 Q = flow_area * barrel_velocity219 220 return Q, barrel_velocity, outlet_culvert_depth221 73 222 74 223 def boyd_circle(self):224 """225 For a circular pipe the Boyd method reviews 3 conditions226 1. Whether the Pipe Inlet is Unsubmerged (acting as weir flow into the inlet)227 2. Whether the Pipe Inlet is Fully Submerged (acting as an Orifice)228 3. Whether the energy loss in the pipe results in the Pipe being controlled by Channel Flow of the Pipe229 230 For these conditions we also would like to assess the pipe flow characteristics as it leaves the pipe231 """232 233 diameter = self.culvert_height234 235 local_debug ='false'236 if self.inflow.get_average_height() > 0.1: #this value was 0.01:237 if local_debug =='true':238 log.critical('Specific E & Deltat Tot E = %s, %s'239 % (str(self.inflow.get_average_specific_energy()),240 str(self.delta_total_energy)))241 log.critical('culvert type = %s' % str(culvert_type))242 # Water has risen above inlet243 244 if self.log_filename is not None:245 s = 'Specific energy = %f m' % self.inflow.get_average_specific_energy()246 log_to_file(self.log_filename, s)247 248 msg = 'Specific energy at inlet is negative'249 assert self.inflow.get_average_specific_energy() >= 0.0, msg250 251 # Calculate flows for inlet control252 Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*self.inflow.get_average_specific_energy()**1.63 # Inlet Ctrl Inlet Unsubmerged253 Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*self.inflow.get_average_specific_energy()**0.63 # Inlet Ctrl Inlet Submerged254 # Note for to SUBMERGED TO OCCUR self.inflow.get_average_specific_energy() should be > 1.2 x diameter.... Should Check !!!255 256 if self.log_filename is not None:257 s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' % (Q_inlet_unsubmerged, Q_inlet_submerged)258 log_to_file(self.log_filename, s)259 Q = min(Q_inlet_unsubmerged, Q_inlet_submerged)260 261 # THE LOWEST Value will Control Calcs From here262 # Calculate Critical Depth Based on the Adopted Flow as an Estimate263 dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)264 dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)265 # From Boyd Paper ESTIMATE of Dcrit has 2 criteria as266 if dcrit1/diameter > 0.85:267 outlet_culvert_depth = dcrit2268 else:269 outlet_culvert_depth = dcrit1270 #outlet_culvert_depth = min(outlet_culvert_depth, diameter)271 # Now determine Hydraulic Radius Parameters Area & Wetted Perimeter272 if outlet_culvert_depth >= diameter:273 outlet_culvert_depth = diameter # Once again the pipe is flowing full not partfull274 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert275 perimeter = diameter * pi276 flow_width= diameter277 case = 'Inlet CTRL Outlet submerged Circular PIPE FULL'278 if local_debug == 'true':279 log.critical('Inlet CTRL Outlet submerged Circular '280 'PIPE FULL')281 else:282 #alpha = acos(1 - outlet_culvert_depth/diameter) # Where did this Come From ????/283 alpha = acos(1-2*outlet_culvert_depth/diameter)*2284 #flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) # Pipe is Running Partly Full at the INLET WHRE did this Come From ?????285 flow_area = diameter**2/8*(alpha - sin(alpha)) # Equation from GIECK 5th Ed. Pg. B3286 flow_width= diameter*sin(alpha/2.0)287 perimeter = alpha*diameter/2.0288 case = 'INLET CTRL Culvert is open channel flow we will for now assume critical depth'289 if local_debug =='true':290 log.critical('INLET CTRL Culvert is open channel flow '291 'we will for now assume critical depth')292 log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'293 % (str(Q), str(outlet_culvert_depth),294 str(alpha)))295 if self.delta_total_energy < self.inflow.get_average_specific_energy(): # OUTLET CONTROL !!!!296 # Calculate flows for outlet control297 298 # Determine the depth at the outlet relative to the depth of flow in the Culvert299 if self.outflow.get_average_height() > diameter: # Outlet is submerged Assume the end of the Pipe is flowing FULL300 outlet_culvert_depth=diameter301 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert302 perimeter = diameter * pi303 flow_width= diameter304 case = 'Outlet submerged'305 if local_debug =='true':306 log.critical('Outlet submerged')307 else: # Culvert running PART FULL for PART OF ITS LENGTH Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity308 # IF self.outflow.get_average_height() < diameter309 dcrit1 = diameter/1.26*(Q/g**0.5*diameter**2.5)**(1/3.75)310 dcrit2 = diameter/0.95*(Q/g**0.5*diameter**2.5)**(1/1.95)311 if dcrit1/diameter >0.85:312 outlet_culvert_depth= dcrit2313 else:314 outlet_culvert_depth = dcrit1315 if outlet_culvert_depth > diameter:316 outlet_culvert_depth = diameter # Once again the pipe is flowing full not partfull317 flow_area = (diameter/2)**2 * pi # Cross sectional area of flow in the culvert318 perimeter = diameter * pi319 flow_width= diameter320 case = 'Outlet unsubmerged PIPE FULL'321 if local_debug =='true':322 log.critical('Outlet unsubmerged PIPE FULL')323 else:324 alpha = acos(1-2*outlet_culvert_depth/diameter)*2325 flow_area = diameter**2/8*(alpha - sin(alpha)) # Equation from GIECK 5th Ed. Pg. B3326 flow_width= diameter*sin(alpha/2.0)327 perimeter = alpha*diameter/2.0328 case = 'Outlet is open channel flow we will for now assume critical depth'329 if local_debug == 'true':330 log.critical('Q Outlet Depth and ALPHA = %s, %s, %s'331 % (str(Q), str(outlet_culvert_depth),332 str(alpha)))333 log.critical('Outlet is open channel flow we '334 'will for now assume critical depth')335 if local_debug == 'true':336 log.critical('FLOW AREA = %s' % str(flow_area))337 log.critical('PERIMETER = %s' % str(perimeter))338 log.critical('Q Interim = %s' % str(Q))339 hyd_rad = flow_area/perimeter340 341 if self.log_filename is not None:342 s = 'hydraulic radius at outlet = %f' %hyd_rad343 log_to_file(self.log_filename, s)344 345 # Outlet control velocity using tail water346 if local_debug =='true':347 log.critical('GOT IT ALL CALCULATING Velocity')348 log.critical('HydRad = %s' % str(hyd_rad))349 culvert_velocity = sqrt(self.delta_total_energy/((self.sum_loss/2/g)+(self.manning**2*self.culvert_length)/hyd_rad**1.33333))350 Q_outlet_tailwater = flow_area * culvert_velocity351 if local_debug =='true':352 log.critical('VELOCITY = %s' % str(culvert_velocity))353 log.critical('Outlet Ctrl Q = %s' % str(Q_outlet_tailwater))354 if self.log_filename is not None:355 s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater356 log_to_file(self.log_filename, s)357 Q = min(Q, Q_outlet_tailwater)358 if local_debug =='true':359 log.critical('%s,%.3f,%.3f'360 % ('dcrit 1 , dcit2 =',dcrit1,dcrit2))361 log.critical('%s,%.3f,%.3f,%.3f'362 % ('Q and Velocity and Depth=', Q,363 culvert_velocity, outlet_culvert_depth))364 365 culv_froude=sqrt(Q**2*flow_width/(g*flow_area**3))366 if local_debug =='true':367 log.critical('FLOW AREA = %s' % str(flow_area))368 log.critical('PERIMETER = %s' % str(perimeter))369 log.critical('Q final = %s' % str(Q))370 log.critical('FROUDE = %s' % str(culv_froude))371 372 # Determine momentum at the outlet373 barrel_velocity = Q/(flow_area + velocity_protection/flow_area)374 375 else: # self.inflow.get_average_height() < 0.01:376 Q = barrel_velocity = outlet_culvert_depth = 0.0377 378 # Temporary flow limit379 if barrel_velocity > self.max_velocity:380 barrel_velocity = self.max_velocity381 Q = flow_area * barrel_velocity382 383 return Q, barrel_velocity, outlet_culvert_depth384 75 385 76 77 78 -
trunk/anuga_core/source/anuga/structures/inlet.py
r7978 r7980 9 9 """ 10 10 11 def __init__(self, domain, polygon ):11 def __init__(self, domain, polygon, outward_culvert_vector=None): 12 12 13 13 self.domain = domain 14 14 self.domain_bounding_polygon = self.domain.get_boundary_polygon() 15 15 self.polygon = polygon 16 self.outward_culvert_vector = outward_culvert_vector 16 17 17 18 # FIXME (SR) Using get_triangle_containing_point which needs to be sped up … … 35 36 assert is_inside_polygon(point, bounding_polygon), msg 36 37 37 self.triangle_indices = inside_polygon(domain_centroids, self.inlet_polygon) 38 39 self.triangle_indices = inside_polygon(domain_centroids, self.inlet_polygon, verbose=True) 38 40 39 41 if len(self.triangle_indices) == 0: … … 42 44 msg += 'specified region: %s' % region 43 45 raise Exception, msg 46 47 44 48 45 49 … … 62 66 63 67 68 def get_area(self): 69 70 return self.area 71 72 64 73 def get_areas(self): 65 74 … … 74 83 75 84 def get_average_stage(self): 76 77 return num.sum(self.get_stages())/self.triangle_indices.size 78 85 86 return num.sum(self.get_stages()*self.get_areas())/self.area 79 87 80 88 def get_elevations(self): … … 84 92 def get_average_elevation(self): 85 93 86 return num.sum(self.get_elevations())/self.triangle_indices.size 94 95 return num.sum(self.get_elevations()*self.get_areas())/self.area 87 96 88 97 … … 93 102 94 103 def get_average_xmom(self): 95 96 return num.sum(self.get_xmoms() )/self.triangle_indices.size104 105 return num.sum(self.get_xmoms()*self.get_areas())/self.area 97 106 98 107 … … 104 113 def get_average_ymom(self): 105 114 106 return num.sum(self.get_ymoms())/self.triangle_indices.size 107 115 return num.sum(self.get_ymoms()*self.get_areas())/self.area 108 116 109 117 def get_heights(self): … … 115 123 116 124 return num.sum(self.get_heights()*self.get_areas()) 117 118 125 126 119 127 def get_average_height(self): 120 128 … … 124 132 def get_velocities(self): 125 133 126 depths = self.get_stages() - self.get_elevations()127 u = self.get_xmoms()/( depths + velocity_protection/depths)128 v = self.get_ymoms()/( depths + velocity_protection/depths)134 heights = self.get_heights() 135 u = self.get_xmoms()/(heights + velocity_protection/heights) 136 v = self.get_ymoms()/(heights + velocity_protection/heights) 129 137 130 138 return u, v 139 140 141 def get_xvelocities(self): 142 143 heights = self.get_heights() 144 return self.get_xmoms()/(heights + velocity_protection/heights) 145 146 def get_yvelocities(self): 147 148 heights = self.get_heights() 149 return self.get_ymoms()/(heights + velocity_protection/heights) 131 150 132 151 … … 135 154 u, v = self.get_velocities() 136 155 137 average_u = num.sum(u )/self.triangle_indices.size138 average_v = num.sum(v )/self.triangle_indices.size156 average_u = num.sum(u*self.get_areas())/self.area 157 average_v = num.sum(v*self.get_areas())/self.area 139 158 140 159 return math.sqrt(average_u**2 + average_v**2)
Note: See TracChangeset
for help on using the changeset viewer.
