# Changeset 5436

Ignore:
Timestamp:
Jun 25, 2008, 5:57:01 PM (15 years ago)
Message:

Moved culvert routines into own area

Files:
1 deleted
6 edited

Unmodified
Removed
• ## anuga_core/source/anuga/config.py

 r5352 default_smoothing_parameter = 0.001 # Default alpha for penalised # least squares fitting velocity_protection = 1.0e-6 #-------------------------------------------
• ## anuga_core/source/anuga/culvert_flows/culvert_routines.py

 r5433 """Collection of culvert routines for use with Culvert_flow in culvert_class Usage: """ #   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet) # # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed # Flow Is Removed at a Rate of INFLOW # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges # # This SHould be changed to a Vertical Opening Both BOX and Circular # There will be several Culvert Routines such as: # CULVERT_Boyd_Channel # CULVERT_Orifice_and_Weir # CULVERT_Simple_FLOOR # CULVERT_Simple_WALL # CULVERT_Eqn_Floor # CULVERT_Eqn_Wall # CULVERT_Tab_Floor # CULVERT_Tab_Wall # BRIDGES..... # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!! #  COULD USE EPA SWMM Model !!!! from math import pi, sqrt, sin, cos def boyd_generalised_culvert_model(culvert, inlet, outlet, delta_Et, g): """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 from anuga.config import velocity_protection from anuga.utilities.numerical_tools import safe_acos as acos 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 log_filename = culvert.log_filename manning = culvert.manning sum_loss = culvert.sum_loss length = culvert.length if inlet.depth >= 0.01: # Calculate driving energy E = inlet.specific_energy s = 'Driving energy  = %f m' %E log_to_file(log_filename, s) msg = 'Driving energy is negative' assert E >= 0.0, msg # Water has risen above inlet if culvert.culvert_type == 'circle': # Round culvert # Calculate flows for inlet control diameter = culvert.diameter Q_inlet_unsubmerged = 0.421*g**0.5*diameter**0.87*E**1.63    # Inlet Ctrl Inlet Unsubmerged Q_inlet_submerged = 0.530*g**0.5*diameter**1.87*E**0.63      # Inlet Ctrl Inlet Submerged s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged) log_to_file(log_filename, s) case = '' if Q_inlet_unsubmerged < Q_inlet_submerged: Q = Q_inlet_unsubmerged alpha = acos(1 - inlet.depth/diameter) flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) outlet_culvert_depth = inlet.depth width = diameter*sin(alpha) perimeter = alpha*diameter case = 'Inlet unsubmerged' else: Q = Q_inlet_submerged flow_area = (diameter/2)**2 * pi outlet_culvert_depth = diameter width = diameter perimeter = diameter 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 > diameter:        # The Outlet is Submerged outlet_culvert_depth=diameter flow_area = (diameter/2)**2 * pi  # Cross sectional area of flow in the culvert perimeter = diameter * pi width = diameter case = 'Outlet submerged' elif outlet.depth==0.0: outlet_culvert_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth alpha = acos(1 - inlet.depth/diameter) flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) perimeter = alpha*diameter width = diameter*sin(alpha) case = 'Outlet depth is zero' else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity outlet_culvert_depth=outlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth alpha = acos(1 - outlet.depth/diameter) flow_area = diameter**2 * (alpha - sin(alpha)*cos(alpha)) perimeter = alpha*diameter width = diameter*sin(alpha) case = 'Outlet is open channel flow' else: # Box culvert (rectangle or square) # Calculate flows for inlet control height = culvert.height width = culvert.width 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(log_filename, s) case = '' if Q_inlet_unsubmerged < Q_inlet_submerged: Q = Q_inlet_unsubmerged flow_area = width*inlet.depth outlet_culvert_depth = inlet.depth perimeter=(width+2.0*inlet.depth) case = 'Inlet unsubmerged' else: Q = Q_inlet_submerged flow_area = width*height outlet_culvert_depth = height perimeter=2.0*(width+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' # Common code for rectangular and circular culvert types hyd_rad = flow_area/perimeter s = 'hydraulic radius at outlet = %f' %hyd_rad log_to_file(log_filename, 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(log_filename, s) Q = min(Q, Q_outlet_tailwater) log_to_file(log_filename, 'Case: "%s"' %case) flow_rate_control=Q s = 'Flow Rate Control = %f' %flow_rate_control log_to_file(log_filename, 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(log_filename, s) # Determine momentum at the outlet barrel_velocity = Q/(flow_area + velocity_protection/flow_area) else: #inlet.depth < 0.01: Q = barrel_velocity = outlet_culvert_depth = 0.0 return Q, barrel_velocity, outlet_culvert_depth
• ## anuga_core/source/anuga/shallow_water/shallow_water_domain.py

 r5315 negative values indicate decreases. Rate can be None at initialisation but must be specified before forcting term is applied (i.e. simulation has started). before forcing term is applied (i.e. simulation has started). center [m]: Coordinates at center of flow point # Update area if applicable self.area = None self.exchange_area = None if center is not None and radius is not None: assert len(center) == 2 assert polygon is None, msg self.area = radius**2*pi self.exchange_area = radius**2*pi if polygon is not None: self.area = polygon_area(self.polygon) self.exchange_area = polygon_area(self.polygon) points = domain.get_centroid_coordinates(absolute=True) self.indices = None # Calculate indices in exchange area for this forcing term self.exchange_indices = None if self.center is not None and self.radius is not None: # Inlet is circular self.indices = [] self.exchange_indices = [] for k in range(N): x, y = points[k,:] # Centroid if ((x-self.center)**2+(y-self.center)**2) < self.radius**2: self.indices.append(k) self.exchange_indices.append(k) if self.polygon is not None: # Inlet is polygon self.indices = inside_polygon(points, self.polygon) self.exchange_indices = inside_polygon(points, self.polygon) if self.indices is None: if self.exchange_indices is None: self.update[:] += rate else: # Brute force assignment of restricted rate for k in self.indices: for k in self.exchange_indices: self.update[k] += rate """Return values for specified quantity restricted to opening """ return self.domain.quantities[self.quantity_name].get_values(indices=self.indices) return self.domain.quantities[self.quantity_name].get_values(indices=self.exchange_indices) """Set values for specified quantity restricted to opening """ self.domain.quantities[self.quantity_name].set_values(val, indices=self.indices) self.domain.quantities[self.quantity_name].set_values(val, indices=self.exchange_indices) if callable(self.rate): _rate = self.rate(t)/self.area _rate = self.rate(t)/self.exchange_area else: _rate = self.rate/self.area _rate = self.rate/self.exchange_area return _rate
• ## anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

 r5304 # Setup only one forcing term, constant rainfall restricted to a polygon enclosing triangle #1 (bce) domain.forcing_terms = [] domain.forcing_terms.append( Rainfall(domain, rate=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]])) R = Rainfall(domain, rate=2.0, polygon = [[1,1], [2,1], [2,2], [1,2]]) assert allclose(R.exchange_area, 1) domain.forcing_terms.append(R) domain.compute_forcing_terms()
• ## anuga_core/source/anuga/utilities/system_tools.py

 r5072 import sys import os def log_to_file(filename, s, verbose=True): """Log string to open file descriptor """ fid = open(filename, 'a') if verbose: print s fid.write(s + '\n') fid.close() def get_user_name():
• ## anuga_work/development/culvert_flow/culvert_boyd_channel.py

 r5429 from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing from anuga.culvert_flows.culvert_polygons import create_culvert_polygons from anuga.culvert_flows.culvert_class import Culvert_flow from anuga.culvert_flows.culvert_routines import boyd_generalised_culvert_model from anuga.utilities.polygon import plot_polygons from anuga.utilities.system_tools import log_to_file # Setup computational domain #------------------------------------------------------------------------------ # Open file for storing some specific results... fid = open('Culvert_Headwall_VarM.txt', 'w') length = 40. domain.set_minimum_storable_height(0.001) s='Size'+str(len(domain)) 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 #------------------------------------------------------------------------------ # Setup specialised forcing terms #------------------------------------------------------------------------------ #   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet) # # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed # Flow Is Removed at a Rate of INFLOW # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges # # This SHould be changed to a Vertical Opening Both BOX and Circular # There will be several Culvert Routines such as: # CULVERT_Boyd_Channel # CULVERT_Orifice_and_Weir # CULVERT_Simple_FLOOR # CULVERT_Simple_WALL # CULVERT_Eqn_Floor # CULVERT_Eqn_Wall # CULVERT_Tab_Floor # CULVERT_Tab_Wall # BRIDGES..... # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!! #  COULD USE EPA SWMM Model !!!! class Culvert_flow: """Culvert flow - transfer water from one hole to another Using Momentum as Calculated by Culvert Flow !! Could be Several Methods Investigated to do This !!! 2008_May_08 To Ole: OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon The two centers are now passed on to create_polygon. Input: Two points, pipe_size (either diameter or width, height), mannings_rougness, inlet/outlet energy_loss_coefficients, internal_bend_coefficent, top-down_blockage_factor and bottom_up_blockage_factor And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront of the culvert At the moment this script uses only Depth, later we can change it to Energy... Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet... """ def __init__(self, domain, label=None, description=None, end_point0=None, end_point1=None, width=None, height=None, manning=None,          # Mannings Roughness for Culvert invert_level0=None,    # Invert level if not the same as the Elevation on the Domain invert_level1=None,    # Invert level if not the same as the Elevation on the Domain loss_exit=None, loss_entry=None, loss_bend=None, loss_special=None, blockage_topdwn=None, blockage_bottup=None, culvert_routine=None, verbose=False): from Numeric import sqrt, sum # Input check if height is None: height = width # Set defaults if manning is None: manning = 0.012   # Set a Default Mannings Roughness for Pipe if loss_exit is None: loss_exit = 1.00 if loss_entry is None: loss_entry = 0.50 if loss_bend is None: loss_bend=0.00 if loss_special is None: loss_special=0.00 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 # Create the fundamental culvert polygons from POLYGON P = create_culvert_polygons(end_point0, end_point1, width=width, height=height) if verbose is True: pass #plot_polygons([[end_point0, end_point1], #               P['exchange_polygon0'], #               P['exchange_polygon1'], #               P['enquiry_polygon0'], #               P['enquiry_polygon1']], #              figname='culvert_polygon_output') self.openings = [] self.openings.append(Inflow(domain, polygon=P['exchange_polygon0'])) self.openings.append(Inflow(domain, polygon=P['exchange_polygon1'])) # Assume two openings for now: Referred to as 0 and 1 assert len(self.openings) == 2 # Store basic geometry self.end_points = [end_point0, end_point1] self.invert_levels = [invert_level0, invert_level1] self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']] self.vector = P['vector'] self.distance = P['length']; assert self.distance > 0.0 self.verbose = verbose self.width = width self.height = height self.last_time = 0.0 self.temp_keep_delta_t = 0.0 # Store hydraulic parameters self.manning = manning self.loss_exit = loss_exit self.loss_entry = loss_entry self.loss_bend = loss_bend self.loss_special = loss_special self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special self.blockage_topdwn = blockage_topdwn self.blockage_bottup = blockage_bottup # Store culvert routine self.culvert_routine = culvert_routine # Create objects to update momentum (a bit crude at this stage) xmom0 = General_forcing(domain, 'xmomentum', polygon=P['exchange_polygon0']) xmom1 = General_forcing(domain, 'xmomentum', polygon=P['exchange_polygon1']) ymom0 = General_forcing(domain, 'ymomentum', polygon=P['exchange_polygon0']) ymom1 = General_forcing(domain, 'ymomentum', polygon=P['exchange_polygon1']) self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ] # Print something s = 'Culvert Effective Length = %.2f m' %(self.distance) log_to_file(fid, s) s = 'Culvert Direction is %s\n' %str(self.vector) log_to_file(fid, s) def __call__(self, domain): from anuga.utilities.numerical_tools import mean from anuga.utilities.polygon import inside_polygon from anuga.config import g, epsilon from Numeric import take, sqrt # Time stuff time = domain.get_time() delta_t = time-self.last_time s = '\nTime = %.2f, delta_t = %f' %(time, delta_t) log_to_file(fid, s) msg = 'Time did not advance' if time > 0.0: assert delta_t > 0.0, msg # Get average water depths at each opening openings = self.openings   # There are two Opening  and  for i, opening in enumerate(openings): stage = domain.quantities['stage'].get_values(location='centroids', indices=opening.exchange_indices) elevation = domain.quantities['elevation'].get_values(location='centroids', indices=opening.exchange_indices) # Indices corresponding to energy enquiry field for this opening coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y) enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) # Get model values for points in enquiry polygon for this opening dq = domain.quantities stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices) xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices) ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices) elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices) depth = stage - elevation # Compute mean values of selected quantitites in the enquiry area in front of the culvert # Epsilon handles a dry cell case ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction) uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction) v = mean(sqrt(ux**2+uy**2))      # Average velocity w = mean(stage)                  # Average stage # Store values at enquiry field opening.velocity = v # Compute mean values of selected quantitites in the exchange area in front of the culvert # Stage and velocity comes from enquiry area and elevation from exchange area # Use invert level instead of elevation if specified invert_level = self.invert_levels[i] if invert_level is not None: z = invert_level else: elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices) z = mean(elevation)                   # Average elevation # Estimated depth above the culvert inlet d = w - z if d < 0.0: # This is possible since w and z are taken at different locations #msg = 'D < 0.0: %f' %d #raise msg d = 0.0 # Ratio of depth to culvert height. # If ratio > 1 then culvert is running full ratio = d/self.height opening.ratio = ratio # Average measures of energy in front of this opening Es = d + 0.5*v**2/g  #  Specific energy in exchange area Et = w + 0.5*v**2/g  #  Total energy in the enquiry area opening.total_energy = Et opening.specific_energy = Es # Store current average stage and depth with each opening object opening.depth = d opening.stage = w opening.elevation = z #################  End of the FOR loop ################################################ # We now need to deal with each opening individually # Determine flow direction based on total energy difference delta_Et = openings.total_energy - openings.total_energy if delta_Et > 0: #print 'Flow U/S ---> D/S' inlet=openings outlet=openings inlet.momentum = self.opening_momentum outlet.momentum = self.opening_momentum else: #print 'Flow D/S ---> U/S' inlet=openings outlet=openings inlet.momentum = self.opening_momentum outlet.momentum = self.opening_momentum delta_Et = -delta_Et msg = 'Total energy difference is negative' assert delta_Et > 0.0, msg delta_h = inlet.stage - outlet.stage delta_z = inlet.elevation - outlet.elevation culvert_slope = (delta_z/self.distance) if culvert_slope < 0.0: # Adverse gradient - flow is running uphill # Flow will be purely controlled by uphill outlet face print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation s = 'Delta total energy = %.3f' %(delta_Et) 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 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.value xmomentum_rate /= delta_t ymomentum_rate = outlet_mom_y - outlet.momentum.value ymomentum_rate /= delta_t 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.rate = xmomentum_rate outlet.momentum.rate = ymomentum_rate # Remember this value for next step (IMPORTANT) outlet.momentum.value = outlet_mom_x outlet.momentum.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) # Execute flow term for each opening # This is where Inflow objects are evaluated and update the domain for opening in self.openings: opening(domain) # Execute momentum terms # This is where Inflow objects are evaluated and update the domain outlet.momentum(domain) outlet.momentum(domain) # Store value of time self.last_time = time #------------------------------------------------------------------------------ # Setup culvert inlets and outlets in current topography # Also Set the Shape & Gap Factors fo rthe Enquiry PolyGons # ALSO Allow the Invert Level to be provided by the USER culvert = Culvert_flow(domain, culvert1 = Culvert_flow(domain, label='Culvert No. 1', description=' This culvert is a test unit 1.2m Wide by 0.75m High', description='This culvert is a test unit 1.2m Wide by 0.75m High', end_point0=[9.0, 2.5], end_point1=[13.0, 2.5], culvert_routine=boyd_generalised_culvert_model, verbose=True) culvert2 = Culvert_flow(domain, label='Culvert No. 2', description='This culvert is a circular test with d=1.2m', end_point0=[9.0, 1.5], end_point1=[30.0, 3.5], diameter=1.20, invert_level0=7, culvert_routine=boyd_generalised_culvert_model, verbose=True) domain.forcing_terms.append(culvert) domain.forcing_terms.append(culvert1) domain.forcing_terms.append(culvert2) #------------------------------------------------------------------------------ temp_keep_delta_t=0.0 for t in domain.evolve(yieldstep = 0.1, finaltime = 20):
Note: See TracChangeset for help on using the changeset viewer.