Context Navigation

source: anuga_core/source/anuga/culvert_flows/culvert_class.py @ 5437

Last change on this file since 5437 was 5437, checked in by ole, 15 years ago
Finished culvert flow work with a good working solution
File size: 15.8 KB

Rev	Line
[5435]	1	from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
	2	from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
	3	from anuga.utilities.system_tools import log_to_file
	4
[5434]	5	class Culvert_flow:
	6	"""Culvert flow - transfer water from one hole to another
	7
	8	Using Momentum as Calculated by Culvert Flow !!
	9	Could be Several Methods Investigated to do This !!!
	10
	11	2008_May_08
	12	To Ole:
	13	OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on
	14	the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon
	15
	16	The two centers are now passed on to create_polygon.
	17
	18
	19	Input: Two points, pipe_size (either diameter or width, height), mannings_rougness,
	20	inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
	21	top-down_blockage_factor and bottom_up_blockage_factor
	22
	23
	24	And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront
	25	of the culvert
	26	At the moment this script uses only Depth, later we can change it to Energy...
	27
	28	Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity
	29	The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
	30
	31	"""
	32
	33	def __init__(self,
	34	domain,
	35	label=None,
	36	description=None,
	37	end_point0=None,
	38	end_point1=None,
	39	width=None,
	40	height=None,
	41	diameter=None,
	42	manning=None, # Mannings Roughness for Culvert
	43	invert_level0=None, # Invert level if not the same as the Elevation on the Domain
	44	invert_level1=None, # Invert level if not the same as the Elevation on the Domain
	45	loss_exit=None,
	46	loss_entry=None,
	47	loss_bend=None,
	48	loss_special=None,
	49	blockage_topdwn=None,
	50	blockage_bottup=None,
	51	culvert_routine=None,
	52	verbose=False):
	53
	54	from Numeric import sqrt, sum
	55
	56	# Input check
	57	if diameter is not None:
	58	self.culvert_type = 'circle'
	59	self.diameter = diameter
	60	if height is not None or width is not None:
	61	msg = 'Either diameter or width&height must be specified, but not both.'
	62	raise Exception, msg
	63	else:
	64	if height is not None:
	65	if width is None:
	66	self.culvert_type = 'square'
	67	width = height
	68	else:
	69	self.culvert_type = 'rectangle'
	70	elif width is not None:
	71	if height is None:
	72	self.culvert_type = 'square'
	73	height = width
	74	else:
	75	msg = 'Either diameter or width&height must be specified.'
	76	raise Exception, msg
	77
	78	if height == width:
	79	self.culvert_type = 'square'
	80
	81	self.height = height
	82	self.width = width
	83
	84
	85	assert self.culvert_type in ['circle', 'square', 'rectangle']
	86
	87	# Set defaults
	88	if manning is None: manning = 0.012 # Set a Default Mannings Roughness for Pipe
	89	if loss_exit is None: loss_exit = 1.00
	90	if loss_entry is None: loss_entry = 0.50
	91	if loss_bend is None: loss_bend=0.00
	92	if loss_special is None: loss_special=0.00
	93	if blockage_topdwn is None: blockage_topdwn=0.00
	94	if blockage_bottup is None: blockage_bottup=0.00
	95	if culvert_routine is None: culvert_routine=boyd_generalised_culvert_model
[5437]	96	if label is None: label = 'culvert_flow'
	97	label += '_' + str(id(self))
[5434]	98
	99	# Open log file for storing some specific results...
	100	self.log_filename = label + '.log'
	101	self.label = label
	102
[5437]	103	# Print something
	104	log_to_file(self.log_filename, self.label)
	105	log_to_file(self.log_filename, description)
	106	log_to_file(self.log_filename, self.culvert_type)
[5434]	107
[5437]	108
	109	# Create the fundamental culvert polygons from POLYGON
	110	if self.culvert_type == 'circle':
	111	# Redefine width and height for use with create_culvert_polygons
	112	width = height = diameter
	113
[5434]	114	P = create_culvert_polygons(end_point0,
	115	end_point1,
	116	width=width,
	117	height=height)
	118
	119	if verbose is True:
	120	pass
	121	#plot_polygons([[end_point0, end_point1],
	122	# P['exchange_polygon0'],
	123	# P['exchange_polygon1'],
	124	# P['enquiry_polygon0'],
	125	# P['enquiry_polygon1']],
	126	# figname='culvert_polygon_output')
	127
	128	self.openings = []
	129	self.openings.append(Inflow(domain,
	130	polygon=P['exchange_polygon0']))
	131
	132	self.openings.append(Inflow(domain,
	133	polygon=P['exchange_polygon1']))
	134
	135
	136	# Assume two openings for now: Referred to as 0 and 1
	137	assert len(self.openings) == 2
	138
	139	# Store basic geometry
	140	self.end_points = [end_point0, end_point1]
	141	self.invert_levels = [invert_level0, invert_level1]
	142	self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
	143	self.vector = P['vector']
	144	self.length = P['length']; assert self.length > 0.0
	145	self.verbose = verbose
	146	self.last_time = 0.0
	147
	148
	149	# Store hydraulic parameters
	150	self.manning = manning
	151	self.loss_exit = loss_exit
	152	self.loss_entry = loss_entry
	153	self.loss_bend = loss_bend
	154	self.loss_special = loss_special
	155	self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
	156	self.blockage_topdwn = blockage_topdwn
	157	self.blockage_bottup = blockage_bottup
	158
	159	# Store culvert routine
	160	self.culvert_routine = culvert_routine
	161
	162
	163	# Create objects to update momentum (a bit crude at this stage)
	164
	165
	166	xmom0 = General_forcing(domain, 'xmomentum',
	167	polygon=P['exchange_polygon0'])
	168
	169	xmom1 = General_forcing(domain, 'xmomentum',
	170	polygon=P['exchange_polygon1'])
	171
	172	ymom0 = General_forcing(domain, 'ymomentum',
	173	polygon=P['exchange_polygon0'])
	174
	175	ymom1 = General_forcing(domain, 'ymomentum',
	176	polygon=P['exchange_polygon1'])
	177
	178	self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
	179
	180
	181	# Print something
	182	s = 'Culvert Effective Length = %.2f m' %(self.length)
	183	log_to_file(self.log_filename, s)
	184
	185	s = 'Culvert Direction is %s\n' %str(self.vector)
	186	log_to_file(self.log_filename, s)
	187
	188	def __call__(self, domain):
	189	from anuga.utilities.numerical_tools import mean
	190	from anuga.utilities.polygon import inside_polygon
	191	from anuga.config import g, epsilon
	192	from Numeric import take, sqrt
	193	from anuga.config import velocity_protection
	194
	195
	196	log_filename = self.log_filename
	197
	198	# Time stuff
	199	time = domain.get_time()
	200	delta_t = time-self.last_time
	201	s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
	202	log_to_file(log_filename, s)
	203
	204	msg = 'Time did not advance'
	205	if time > 0.0: assert delta_t > 0.0, msg
	206
	207
	208	# Get average water depths at each opening
	209	openings = self.openings # There are two Opening [0] and [1]
	210	for i, opening in enumerate(openings):
	211	stage = domain.quantities['stage'].get_values(location='centroids',
	212	indices=opening.exchange_indices)
	213	elevation = domain.quantities['elevation'].get_values(location='centroids',
	214	indices=opening.exchange_indices)
	215
	216	# Indices corresponding to energy enquiry field for this opening
	217	coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
	218	enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i])
	219
	220
	221	# Get model values for points in enquiry polygon for this opening
	222	dq = domain.quantities
	223	stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
	224	xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices)
	225	ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices)
	226	elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
	227	depth = stage - elevation
	228
	229	# Compute mean values of selected quantitites in the enquiry area in front of the culvert
	230	# Epsilon handles a dry cell case
	231	ux = xmomentum/(depth+velocity_protection/depth) # Velocity (x-direction)
	232	uy = ymomentum/(depth+velocity_protection/depth) # Velocity (y-direction)
	233	v = mean(sqrt(ux2+uy2)) # Average velocity
	234	w = mean(stage) # Average stage
	235
	236	# Store values at enquiry field
	237	opening.velocity = v
	238
	239
	240	# Compute mean values of selected quantitites in the exchange area in front of the culvert
	241	# Stage and velocity comes from enquiry area and elevation from exchange area
	242
	243	# Use invert level instead of elevation if specified
	244	invert_level = self.invert_levels[i]
	245	if invert_level is not None:
	246	z = invert_level
	247	else:
	248	elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
	249	z = mean(elevation) # Average elevation
[5437]	250
[5434]	251	# Estimated depth above the culvert inlet
[5437]	252	d = w - z # Used for calculations involving depth
[5434]	253	if d < 0.0:
	254	# This is possible since w and z are taken at different locations
	255	#msg = 'D < 0.0: %f' %d
	256	#raise msg
	257	d = 0.0
	258
[5437]	259
	260
	261	# Depth at exchange area used to trigger calculations
	262	stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
	263	elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
	264	depth = stage - elevation
	265	d_trigger = mean(depth)
	266
	267
	268
[5434]	269	# Ratio of depth to culvert height.
	270	# If ratio > 1 then culvert is running full
[5437]	271	if self.culvert_type == 'circle':
	272	ratio = d/self.diameter
	273	else:
	274	ratio = d/self.height
[5434]	275	opening.ratio = ratio
	276
	277	# Average measures of energy in front of this opening
	278	Es = d + 0.5v*2/g # Specific energy in exchange area
	279	Et = w + 0.5v*2/g # Total energy in the enquiry area
	280	opening.total_energy = Et
	281	opening.specific_energy = Es
	282
[5437]	283	# Store current average stage and depth with each opening object
[5434]	284	opening.depth = d
[5437]	285	opening.depth_trigger = d_trigger
[5434]	286	opening.stage = w
	287	opening.elevation = z
	288
	289
	290	################# End of the FOR loop ################################################
	291
	292
	293	# We now need to deal with each opening individually
	294
	295	# Determine flow direction based on total energy difference
	296	delta_Et = openings[0].total_energy - openings[1].total_energy
	297
	298	if delta_Et > 0:
	299	#print 'Flow U/S ---> D/S'
	300	inlet=openings[0]
	301	outlet=openings[1]
	302
	303	inlet.momentum = self.opening_momentum[0]
	304	outlet.momentum = self.opening_momentum[1]
	305
	306	else:
	307	#print 'Flow D/S ---> U/S'
	308	inlet=openings[1]
	309	outlet=openings[0]
	310
	311	inlet.momentum = self.opening_momentum[1]
	312	outlet.momentum = self.opening_momentum[0]
	313
	314	delta_Et = -delta_Et
	315
	316	msg = 'Total energy difference is negative'
	317	assert delta_Et > 0.0, msg
	318
	319	delta_h = inlet.stage - outlet.stage
	320	delta_z = inlet.elevation - outlet.elevation
	321	culvert_slope = (delta_z/self.length)
	322
	323	if culvert_slope < 0.0:
	324	# Adverse gradient - flow is running uphill
	325	# Flow will be purely controlled by uphill outlet face
	326	print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
	327
	328
	329	s = 'Delta total energy = %.3f' %(delta_Et)
	330	log_to_file(log_filename, s)
	331
	332
	333	Q, barrel_velocity, culvert_outlet_depth = self.culvert_routine(self, inlet, outlet, delta_Et, g)
	334	#####################################################
	335	barrel_momentum = barrel_velocity*culvert_outlet_depth
	336
	337	s = 'Barrel velocity = %f' %barrel_velocity
	338	log_to_file(log_filename, s)
	339
	340	# Compute momentum vector at outlet
	341	outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
	342
	343	s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
	344	log_to_file(log_filename, s)
	345
	346	delta_t = time - self.last_time
	347	if delta_t > 0.0:
	348	xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
	349	xmomentum_rate /= delta_t
	350
	351	ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
	352	ymomentum_rate /= delta_t
	353
	354	s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
	355	log_to_file(log_filename, s)
	356	else:
	357	xmomentum_rate = ymomentum_rate = 0.0
	358
	359
	360	# Set momentum rates for outlet jet
	361	outlet.momentum[0].rate = xmomentum_rate
	362	outlet.momentum[1].rate = ymomentum_rate
	363
	364	# Remember this value for next step (IMPORTANT)
	365	outlet.momentum[0].value = outlet_mom_x
	366	outlet.momentum[1].value = outlet_mom_y
	367
	368	if int(domain.time*100) % 100 == 0:
	369	s = 'T=%.5f, Culvert Discharge = %.3f f'\
	370	%(time, Q)
	371	s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\
	372	%(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)
	373	s += ' Momentum rate: (%.4f, %.4f)'\
	374	%(xmomentum_rate, ymomentum_rate)
	375	s+='Outlet Vel= %.3f'\
	376	%(barrel_velocity)
	377	log_to_file(log_filename, s)
	378
	379
	380
	381
	382
	383	# Execute flow term for each opening
	384	# This is where Inflow objects are evaluated and update the domain
	385	for opening in self.openings:
	386	opening(domain)
	387
	388	# Execute momentum terms
	389	# This is where Inflow objects are evaluated and update the domain
	390	outlet.momentum[0](domain)
	391	outlet.momentum[1](domain)
	392
	393	# Store value of time
	394	self.last_time = time
	395
	396

Note: See TracBrowser for help on using the repository browser.

Download in other formats: