# Changeset 5429

Ignore:
Timestamp:
Jun 25, 2008, 4:15:10 PM (15 years ago)
Message:

factored boyd culvert routine into separate function

File:
1 edited

Unmodified
Removed
• ## anuga_work/development/culvert_flow/culvert_boyd_channel.py

 r5428 from anuga.culvert_flows.culvert_polygons import create_culvert_polygons from anuga.utilities.polygon import plot_polygons from anuga.utilities.system_tools import log_to_file from math import pi,sqrt #------------------------------------------------------------------------------ def log(fid, s): print s fid.write(s + '\n') # Open file for storing some specific results... s='Size'+str(len(domain)) log(fid, s) velocity_protection = 1.0e-4 log_to_file(fid, s) velocity_protection = 1.0e-5 def boyd_generalised_culvert_model(inlet, outlet, delta_Et, width, height, length, sum_loss, manning, g, fid): """Boyd's generalisation of the US department of transportation culvert model # == The quantity of flow passing through a culvert is controlled by many factors # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled """ from anuga.utilities.system_tools import log_to_file Q_outlet_tailwater = 0.0 inlet.rate = 0.0 outlet.rate = 0.0 Q_inlet_unsubmerged = 0.0 Q_inlet_submerged = 0.0 Q_outlet_critical_depth = 0.0 if inlet.depth >= 0.01: # Water has risen above inlet if width==height:  # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ??? pass #Q1[t]= 0.421*g**0.5*Diam**0.87*HW**1.63    # Inlet Ctrl Inlet Unsubmerged #Q2[t]= 0.530*g**0.5*Diam**1.87*HW**0.63    # Inlet Ctrl Inlet Submerged #PipeDcrit= #Delta_E=HW-TW else: # Box culvert # Calculate flows for inlet control E = inlet.specific_energy s = 'Driving energy  = %f m' %E log_to_file(fid, s) msg = 'Driving energy is negative' assert E >= 0.0, msg Q_inlet_unsubmerged = 0.540*g**0.5*width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged Q_inlet_submerged = 0.702*g**0.5*width*height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged) log_to_file(fid, s) case = '' if Q_inlet_unsubmerged < Q_inlet_submerged: Q = Q_inlet_unsubmerged flow_area = width*inlet.depth outlet_culvert_depth = inlet.depth case = 'Inlet unsubmerged' else: Q = Q_inlet_submerged flow_area = width*height outlet_culvert_depth = height case = 'Inlet submerged' if delta_Et < E: # Calculate flows for outlet control # Determine the depth at the outlet relative to the depth of flow in the Culvert if outlet.depth > height:        # The Outlet is Submerged outlet_culvert_depth=height flow_area=width*height       # Cross sectional area of flow in the culvert perimeter=2.0*(width+height) case = 'Outlet submerged' elif outlet.depth==0.0: outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth flow_area=width*inlet.depth perimeter=(width+2.0*inlet.depth) case = 'Outlet depth is zero' else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity outlet_culvert_depth=outlet.depth flow_area=width*outlet.depth perimeter=(width+2.0*outlet.depth) case = 'Outlet is open channel flow' hyd_rad = flow_area/perimeter s = 'hydraulic radius at outlet = %f' %hyd_rad log_to_file(fid, s) # Outlet control velocity using tail water culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*length)/hyd_rad**1.33333)) Q_outlet_tailwater = flow_area * culvert_velocity s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater log_to_file(fid, s) Q = min(Q, Q_outlet_tailwater) log_to_file(fid, 'Case: "%s"' %case) flow_rate_control=Q s = 'Flow Rate Control = %f' %flow_rate_control log_to_file(fid, s) inlet.rate = -flow_rate_control outlet.rate = flow_rate_control culv_froude=sqrt(flow_rate_control**2*width/(g*flow_area**3)) s = 'Froude in Culvert = %f' %culv_froude log_to_file(fid, s) # Determine momentum at the outlet barrel_velocity = Q/(flow_area + velocity_protection/flow_area) return Q, barrel_velocity, outlet_culvert_depth z[i] += embank_hgt-embank_dnslope*(x[i] -12.0)       # Sloping D/S Face # Constriction #if 27 < x[i] < 29 and y[i] > 3: #    z[i] += 2 # Pole #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2: #    z[i] += 2 # HOLE For Pit at Opening[0] #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2: #  z[i]-=1 # HOLE For Pit at Opening[1] #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2: #  z[i]-=1 return z blockage_topdwn=None, blockage_bottup=None, culvert_routine=None, verbose=False): if blockage_topdwn is None: blockage_topdwn=0.00 if blockage_bottup is None: blockage_bottup=0.00 if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model self.blockage_bottup = blockage_bottup # Store culvert routine self.culvert_routine = culvert_routine # Create objects to update momentum (a bit crude at this stage) # Print something s = 'Culvert Effective Length = %.2f m' %(self.distance) log(fid, s) log_to_file(fid, s) s = 'Culvert Direction is %s\n' %str(self.vector) log(fid, s) log_to_file(fid, s) def __call__(self, domain): delta_t = time-self.last_time s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) log(fid, s) log_to_file(fid, s) msg = 'Time did not advance' enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) if self.verbose: pass #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\ #      %(i, len(idx), self.enquiry_polygons[i]) # Get model values for points in enquiry polygon for this opening inlet.momentum = self.opening_momentum[0] outlet.momentum = self.opening_momentum[1] outlet.momentum = self.opening_momentum[1] else: #print 'Flow D/S ---> U/S' s = 'Delta total energy = %.3f' %(delta_Et) log(fid, s) # ====================== NOW ENTER INTO THE CULVERT EQUATIONS AS DEFINED BY BOYD GENERALISED METHOD # == The quantity of flow passing through a culvert is controlled by many factors # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled Q_outlet_tailwater = 0.0 inlet.rate = 0.0 outlet.rate = 0.0 Q_inlet_unsubmerged = 0.0 Q_inlet_submerged = 0.0 Q_outlet_critical_depth = 0.0 if inlet.depth >= 0.01: # Water has risen above inlet if self.width==self.height:  # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ??? pass #Q1[t]= 0.421*g**0.5*Diam**0.87*HW**1.63    # Inlet Ctrl Inlet Unsubmerged #Q2[t]= 0.530*g**0.5*Diam**1.87*HW**0.63    # Inlet Ctrl Inlet Submerged #PipeDcrit= #Delta_E=HW-TW else: # Box culvert log_to_file(fid, s) sum_loss = self.sum_loss manning=self.manning length = self.distance height = self.height width = self.width Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(inlet, outlet, delta_Et, width, height, length, sum_loss, manning, g, fid) ##################################################### barrel_momentum = barrel_velocity*culvert_outlet_depth sum_loss=self.loss_exit+self.loss_entry+self.loss_bend+self.loss_special manning=self.manning # Calculate flows for inlet control E = inlet.specific_energy #E = min(inlet.specific_energy, delta_Et) s = 'Driving energy  = %f m' %E log(fid, s) msg = 'Driving energy is negative' assert E >= 0.0, msg Q_inlet_unsubmerged = 0.540*g**0.5*self.width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged Q_inlet_submerged = 0.702*g**0.5*self.width*self.height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged) log(fid, s) case = '' if Q_inlet_unsubmerged < Q_inlet_submerged: Q = Q_inlet_unsubmerged flow_area = self.width*inlet.depth outlet_culv_depth = inlet.depth case = 'Inlet unsubmerged' else: Q = Q_inlet_submerged flow_area = self.width*self.height outlet_culv_depth = self.height case = 'Inlet submerged' if delta_Et < Es: # Calculate flows for outlet control # Determine the depth at the outlet relative to the depth of flow in the Culvert sum_loss = self.sum_loss s = 'Barrel velocity = %f' %barrel_velocity log_to_file(fid, s) # Compute momentum vector at outlet outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) log_to_file(fid, s) delta_t = time - self.last_time if delta_t > 0.0: xmomentum_rate = outlet_mom_x - outlet.momentum[0].value xmomentum_rate /= delta_t if outlet.depth > self.height:  # The Outlet is Submerged outlet_culv_depth=self.height flow_area=self.width*self.height # Cross sectional area of flow in the culvert perimeter=2.0*(self.width+self.height) case = 'Outlet submerged' elif outlet.depth==0.0: outlet_culv_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth flow_area=self.width*inlet.depth perimeter=(self.width+2.0*inlet.depth) case = 'Outlet depth is zero' else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity outlet_culv_depth=outlet.depth flow_area=self.width*outlet.depth perimeter=(self.width+2.0*outlet.depth) case = 'Outlet is open channel flow' hyd_rad = flow_area/perimeter s = 'hydraulic radius at outlet = %f' %hyd_rad log(fid, s) ymomentum_rate = outlet_mom_y - outlet.momentum[1].value ymomentum_rate /= delta_t # Outlet control velocity using tail water culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333)) Q_outlet_tailwater = flow_area * culvert_velocity s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater log(fid, s) Q = min(Q, Q_outlet_tailwater) log(fid, 'Case: "%s"' %case) s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate) log_to_file(fid, s) else: xmomentum_rate = ymomentum_rate = 0.0 # Set momentum rates for outlet jet outlet.momentum[0].rate = xmomentum_rate outlet.momentum[1].rate = ymomentum_rate # Remember this value for next step (IMPORTANT) outlet.momentum[0].value = outlet_mom_x outlet.momentum[1].value = outlet_mom_y if int(domain.time*100) % 100 == 0: s = 'T=%.5f, Culvert Discharge = %.3f f'\ %(time, Q) s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\ %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y) s += ' Momentum rate: (%.4f, %.4f)'\ %(xmomentum_rate, ymomentum_rate) s+='Outlet Vel= %.3f'\ %(barrel_velocity) log_to_file(fid, s) flow_rate_control=Q s = 'Flow Rate Control = %f' %flow_rate_control log(fid, s) inlet.rate = -flow_rate_control outlet.rate = flow_rate_control culv_froude=sqrt(flow_rate_control**2*self.width/(g*flow_area**3)) s = 'Froude in Culvert = %f' %culv_froude log(fid, s) # Determine momentum at the outlet barrel_velocity = Q/(flow_area + velocity_protection/flow_area) barrel_momentum = barrel_velocity*outlet_culv_depth outlet_E_Loss= 0.5*0.5*barrel_velocity**2/g   # Ke v^2/2g s = 'Barrel velocity = %f' %barrel_velocity log(fid, s) # Compute momentum vector at outlet outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) log(fid, s) delta_t = time - self.last_time if delta_t > 0.0: xmomentum_rate = outlet_mom_x - outlet.momentum[0].value xmomentum_rate /= delta_t ymomentum_rate = outlet_mom_y - outlet.momentum[1].value ymomentum_rate /= delta_t s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate) log(fid, s) else: xmomentum_rate = ymomentum_rate = 0.0 # Set momentum rates for outlet jet outlet.momentum[0].rate = xmomentum_rate outlet.momentum[1].rate = ymomentum_rate # Remember this value for next step (IMPORTANT) outlet.momentum[0].value = outlet_mom_x outlet.momentum[1].value = outlet_mom_y if int(domain.time*100) % 100 == 0: s = 'T=%.5f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\ %(time, flow_rate_control, barrel_velocity) s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\ %(outlet_culv_depth, outlet_mom_x,outlet_mom_y) s += ' Momentum rate: (%.4f, %.4f)'\ %(xmomentum_rate, ymomentum_rate) s+='Outlet Vel= %.3f'\ %(barrel_velocity) log(fid, s) # Execute flow term for each opening end_point0=[9.0, 2.5], end_point1=[13.0, 2.5], width=1.20,height=0.75, width=1.20,height=0.75, culvert_routine=boyd_generalised_culvert_model, verbose=True)
Note: See TracChangeset for help on using the changeset viewer.