cleaned up code, added chamfer

flexure/flexure_stiffness.py

@@ -5,6 +5,8 @@ import argparse
 from pyframe3dd.frame3dd import write_frame3dd_file, read_lowest_mode, read_frame3dd_displacements, compute_mass
 from pyframe3dd.util import magnitudes, close
 import subprocess
+import matplotlib.pyplot as plt
+plt.style.use('bmh')
 
 def plot_connections(nodes,beamsets):
     #for debug only, this is slow!
@@ -45,49 +47,61 @@ def clean_up_frame3dd(filename):
 def build(args):
 	#return nodes,rods as numpy arrays
 	dxy = args.attach_radius/sqrt(2)
-
-	if args.flexure_type == 'mirrored':
-		nodes = args.l*array([[0,0,0],[0,1,0],[-1,1,0],[1,0,0],[1,-1,0]]) #one side of flexure plate
-		beams = array([[0,1],[1,2],[0,3],[3,4]])
-		nodes += array([dxy,dxy,0]) #offset by attachment radius
-
-	elif args.flexure_type == 'cyclic':
-		nodes = args.l*array([[0,0,0],[0,1,0],[-1,1,0]]) + array([dxy,dxy,0])
-		beams = array([[0,1],[1,2],[3,4]])
-		nodes = vstack((nodes, args.l*array([[1,1,0],[1,0,0]]) + array([dxy,-dxy,0])))
-
-	nodes = vstack((nodes,array([[-n[0],-n[1],0] for n in nodes]))) #append reflection
-	beams = vstack((beams, beams + 5))
+	nodes = array([[dxy,dxy,0],[-dxy,dxy,0],[-dxy,-dxy,0],[dxy,-dxy,0]])
+	solid_beams = array([ [0,1],[1,2],[2,3],[3,0],[0,2],[1,3] ])
 	z_os = array([0,0,.5*args.sep])
 	nodes = vstack(( nodes+z_os, nodes-z_os ))
-	beams = vstack((beams, beams + 10))
-
-	solid_nodes = array([ [dxy,-dxy,.5*args.sep],[-dxy,dxy,.5*args.sep],[dxy,-dxy,-.5*args.sep],[-dxy,dxy,-.5*args.sep] ])
-	nodes = vstack((nodes, solid_nodes))
-	solid_beams = array([
-		[0,5],[0,10],[5,15],[10,15],#[0,15],[5,10],
-		[0,20],[0,21],[5,20],[5,21],[20,21],
-		[10,22],[10,23],[15,22],[15,23],[22,23],
-		[20,22],[21,23],#[20,23],[21,22],
-		[0,22],[0,23],[5,22],[5,23],
-		[10,20],[10,21],[15,20],[15,21]
-		])
+	solid_beams = vstack((solid_beams, solid_beams + 4))
+	solid_beams = vstack((solid_beams, array([ 
+		[0,4],[1,5],[2,6],[3,7],
+		[0,5],[1,6],[2,7],[3,4],
+		[0,7],[1,4],[2,5],[3,6] ])))
+	#sensor nodes
+	nodes = vstack((nodes, array([ [args.sensor_radius,0,0],[0,args.sensor_radius,0],[-args.sensor_radius,0,0],[0,-args.sensor_radius,0] ])))
+	solid_beams = vstack((solid_beams, array([ [0,8],[3,8],[0,9],[1,9],[1,10],[2,10],[2,11],[3,11] ]) ))
+	solid_beams = vstack((solid_beams, array([ [4,8],[7,8],[4,9],[5,9],[5,10],[6,10],[6,11],[7,11] ]) ))
+
 	if args.flexure_type == 'cyclic':
-		beams = vstack((beams, array([[4,20],[9,21],[14,22],[19,23]])))
+		if args.chamfer > 0:
+			l = args.l; ch = args.chamfer
+			#chamfer the flexure	
+			flexure_nodes = array([[0,l-ch,0],[-ch,l,0],[-l,l,0]]) + array([dxy,dxy,0])
+			flexure_nodes = vstack(( flexure_nodes, array([[l-ch,0,0],[l,ch,0],[l,l,0]]) + array([dxy,-dxy,0]) ))
+			#append reflection
+			flexure_nodes = vstack((flexure_nodes,array([[-n[0],-n[1],0] for n in flexure_nodes])))
+			flexure_beams = array([[0,12],[12,13],[13,14],[3,15],[15,16],[16,17],[2,18],[18,19],[19,20],[1,21],[21,22],[22,23]])
+			#append both plates
+			flexure_nodes = vstack(( flexure_nodes + z_os, flexure_nodes - z_os))
+			flexure_beams = vstack(( flexure_beams, array([[4,24],[24,25],[25,26],[7,27],[27,28],[28,29],[6,30],[30,31],[31,32],[5,33],[33,34],[34,35]]) ))
+			fixed_nodes = [14,17,20,23,26,29,32,35]
+		else:
+			flexure_nodes = args.l*array([[0,1,0],[-1,1,0]]) + array([dxy,dxy,0])
+			flexure_nodes = vstack(( flexure_nodes, args.l*array([[1,0,0],[1,1,0]]) + array([dxy,-dxy,0]) ))
+			#append reflection
+			flexure_nodes = vstack((flexure_nodes,array([[-n[0],-n[1],0] for n in flexure_nodes])))
+			flexure_beams = array([[0,12],[12,13],[3,14],[14,15],[2,16],[16,17],[1,18],[18,19]])
+			#append both plates
+			flexure_nodes = vstack(( flexure_nodes + z_os, flexure_nodes - z_os))
+			flexure_beams = vstack(( flexure_beams, array([[4,20],[20,21],[7,22],[22,23],[6,24],[24,25],[5,26],[26,27]]) ))
+			fixed_nodes = [13,15,17,19,21,23,25,27]
+	elif args.flexure_type == 'mirrored':
+		flexure_nodes = args.l*array([[0,1,0],[-1,1,0]]) + array([dxy,dxy,0])
+		flexure_nodes = vstack(( flexure_nodes, args.l*array([[1,0,0],[1,-1,0]]) + array([dxy,dxy,0]) ))
+		#append reflection
+		flexure_nodes = vstack((flexure_nodes,array([[-n[0],-n[1],0] for n in flexure_nodes])))
+		flexure_beams = array([[0,12],[12,13],[0,14],[14,15],[2,16],[16,17],[2,18],[18,19]])
+		#append both plates
+		flexure_nodes = vstack(( flexure_nodes + z_os, flexure_nodes - z_os))
+		flexure_beams = vstack(( flexure_beams, array([[4,20],[20,21],[4,22],[22,23],[6,24],[24,25],[6,26],[26,27]]) ))
+		fixed_nodes = [13,15,17,19,21,23,25,27]
 
-	nodes = vstack((nodes, array([ [args.sensor_radius,0,0],[0,args.sensor_radius,0],[-args.sensor_radius,0,0],[0,-args.sensor_radius,0] ])))
-	solid_beams = vstack(( solid_beams, array([ 
-		[24,0],[24,10],[24,20],[24,22],
-		[25,0],[25,10],[25,21],[25,23],
-		[26,5],[26,15],[26,21],[26,23],
-		[27,5],[27,15],[27,20],[27,22]
-		])))
 
-	return nodes, beams, solid_beams
+	nodes = vstack((nodes, flexure_nodes))
+	return nodes, flexure_beams, solid_beams, fixed_nodes
 
 def run_simulation(args):
 	#set up simulation
-	nodes,beams,solid_beams = build(args)
+	nodes,beams,solid_beams,fixed_nodes = build(args)
 	global_args = {
 		'n_modes':args.n_modes,'length_scaling':args.length_scaling,'exagerration':10,
 		'zoom_scale':2.,'node_radius':zeros(shape(nodes)[0]),
@@ -98,15 +112,14 @@ def run_simulation(args):
 	(beams,{'E':args.E,'nu':args.nu,'rho':args.rho,'cross_section':'rectangular','d2':args.w,'d1':args.t,'roll':0.,'loads':[],'beam_divisions':args.bd,'prestresses':[]}),
 	(solid_beams,{'E':10*args.E,'nu':args.nu,'rho':args.rho,'cross_section':'rectangular','d1':.003,'d2':.003,'roll':0.,'loads':[],'beam_divisions':1,'prestresses':[]})
 	]
-	if args.flexure_type == 'mirrored':
-		fixed_nodes = [2,4,7,9,12,14,17,19]
-	elif args.flexure_type == 'cyclic':
-		fixed_nodes = [2,3,7,8,12,13,17,18]
 	constraints = [{'node':node,'DOF':dof,'value':0} for dof in [0,1,2,3,4,5] for node in fixed_nodes]
 
 	
-	loaded_nodes = [0,5,10,15,20,21,22,23]
-	sensor_nodes = [24,25,26,27]
+	#loaded_nodes = [0,5,10,15,20,21,22,23]
+	#sensor_nodes = [24,25,26,27]
+	loaded_nodes = range(8)
+	sensor_nodes = [8,9,10,11]
+
 	results = []
 
 	for force_dof in [0,1,2]:
@@ -125,7 +138,7 @@ def run_simulation(args):
 		loads = [{'node':n,'DOF':torque_dof,'value':torque_force if nodes[n][2]>0 else -torque_force} for n in loaded_nodes]
 		run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads)
 		disps = read_frame3dd_displacements(global_args['frame3dd_filename'])
-		moving_sensor_nodes = [24,26] if torque_dof==0 else [25,27]
+		moving_sensor_nodes = [8,10] if torque_dof==0 else [9,11]
 
 		#print disps[sensor_nodes]
 		#axis = (array([0,0,0]), array([0,-1,0]) if torque_dof==0 else array([1,0,0]) )
@@ -159,45 +172,9 @@ def run_simulation(args):
 	#todo: plot displacements vs. design parameters
 	return results
 
-def find_stability_threshold(args):
-	#out loop of simulations to determine the buckling load
-	lower = 0 #lower bound
-	upper = 10*args.force_res #initial upper bound before bracketing
-	bracketed=False
-	#actually not necessary, but fun to have the unloaded frequency
-	args.force = lower
-	res = run_simulation(args)	
-	freqs = [res['fundamental_frequency']]
-	forces = [args.force]
-
-	i = 0
-	while not bracketed:
-		print lower,upper,bracketed,res['fundamental_frequency']
-		args.force = upper
-		res = run_simulation(args); i += 1
-		if res['fundamental_frequency']<0:
-			bracketed=True
-		else:
-			freqs.append(res['fundamental_frequency'])
-			forces.append(args.force)
-			lower = upper
-			upper = 2*upper
-	while (upper-lower > args.force_res):
-		print lower,upper,bracketed
-		args.force = .5*(upper+lower)
-		res = run_simulation(args); i += 1
-		if res['fundamental_frequency']>0:
-			freqs.append(res['fundamental_frequency'])
-			forces.append(args.force)
-			lower = .5*(upper+lower)
-		else:
-			upper = .5*(upper+lower)
-	return forces,freqs,res
-
-
 if __name__ == '__main__':
 	parser = argparse.ArgumentParser()
-	parser.add_argument('-M','--mode',choices=('simulate','search', 'visualize'), required=True)
+	parser.add_argument('-M','--mode',choices=('simulate','graph', 'visualize'), required=True)
 	parser.add_argument('-flexure_type','--flexure_type',choices=('cyclic','mirrored'), required=True)
 	parser.add_argument('-Q','--quiet',action='store_true',help='Whether to suppress frame3dd output')
 	parser.add_argument("-f","--force",  type=double, default=.1, help="force to apply (N)")
@@ -208,8 +185,9 @@ if __name__ == '__main__':
 	parser.add_argument("-t","--t", type=double, default=.0023, help="thickness of flexure material (m)")
 	parser.add_argument("-l","--l", type=double, default=.0068, help="length of flexure segment (m)")
 	parser.add_argument("-attach_radius","--attach_radius", type=double, default=.0043, help="distance from z axis to flexure attachment (m)")
-	parser.add_argument("-sep","--sep", type=double, default=.025, help="flexure plate z separation (m)")
+	parser.add_argument("-sep","--sep", type=double, default=.0185, help="flexure plate z separation (m)")
 	parser.add_argument("-sensor_radius","--sensor_radius", type=double, default=.012, help="distance from rotation axis to sensor (m)")
+	parser.add_argument("-chamfer","--chamfer", type=double, default=.001, help="chamfer length for flexure (m), zero for no chamfer")
 
 	parser.add_argument("-bd","--bd", type=int, default=1, help='how many divisions for each rod, useful in buckling analysis')
 	parser.add_argument("-E","--E", type=double, default=70e9, help="Young's Modulus of laminate")
@@ -220,21 +198,83 @@ if __name__ == '__main__':
 	parser.add_argument("-ls","--length_scaling", type=double, default=1., help="Scale factor to keep numbers commesurate")
 	args = parser.parse_args()
 
-	if args.mode=='search':
-		forces,freqs,last_res = find_stability_threshold(args)
-		print "Fundamental frequency: %.3f Hz"%(freqs[-1])
-		print "Critical force: %.3f N"%(forces[-1])
-		print "Critical stress: %.3f MPa"%(last_res['stress']/1e6)
+	#if args.mode=='search':
+	#	forces,freqs,last_res = find_stability_threshold(args)
+	#	print "Fundamental frequency: %.3f Hz"%(freqs[-1])
+	#	print "Critical force: %.3f N"%(forces[-1])
+	#	print "Critical stress: %.3f MPa"%(last_res['stress']/1e6)
+	if args.mode == 'graph':
+		ws = linspace(.000, .9*args.l, 10)
+		res = {}
+		for wi in ws:
+			#args.w = wi
+			args.chamfer = wi
+			res[wi] = run_simulation(args)
+
+		X = [1e6*res[wi][0]['displacement'] for wi in ws]
+		Y = [1e6*res[wi][1]['displacement'] for wi in ws]
+		Z = [1e6*res[wi][2]['displacement'] for wi in ws]
+		rX = [1e6*res[wi][3]['displacement'] for wi in ws]
+		rY = [1e6*res[wi][4]['displacement'] for wi in ws]
+		rZ = [1e6*res[wi][5]['displacement'] for wi in ws]
+		print ws,X,Y,Z
+		plt.plot( 1e3*ws, X, label='X' )
+		plt.plot( 1e3*ws, Y, label='Y' )
+		plt.plot( 1e3*ws, Z, label='Z' )
+		plt.plot( 1e3*ws, rX, label='rX' )
+		plt.plot( 1e3*ws, rY, label='rY' )
+		plt.plot( 1e3*ws, rZ, label='rZ' )
+		plt.ylabel('displacement at sensor (microns)')
+		plt.xlabel('chamfer width (mm)')
+		plt.xlim([1e3*ws[0],1e3*ws[-1]])
+		plt.legend(loc='upper right')
+		plt.show()
 	elif args.mode=='simulate':
 		res = run_simulation(args)
 		print res
 		#print "Fundamental frequency: %.3f Hz"%res['fundamental_frequency']
 		#print "Stress: %.3f MPa"%(res['stress']/1e6)
 	elif args.mode=='visualize':
-		nodes,rods,solid_beams = build(args)
+		nodes,rods,solid_beams,fixed_nodes = build(args)
 		plot_connections(nodes,[rods,solid_beams])
 	else:
 		assert(0) #should not be here
 
 
 
+'''
+def find_stability_threshold(args):
+	#out loop of simulations to determine the buckling load
+	lower = 0 #lower bound
+	upper = 10*args.force_res #initial upper bound before bracketing
+	bracketed=False
+	#actually not necessary, but fun to have the unloaded frequency
+	args.force = lower
+	res = run_simulation(args)	
+	freqs = [res['fundamental_frequency']]
+	forces = [args.force]
+
+	i = 0
+	while not bracketed:
+		print lower,upper,bracketed,res['fundamental_frequency']
+		args.force = upper
+		res = run_simulation(args); i += 1
+		if res['fundamental_frequency']<0:
+			bracketed=True
+		else:
+			freqs.append(res['fundamental_frequency'])
+			forces.append(args.force)
+			lower = upper
+			upper = 2*upper
+	while (upper-lower > args.force_res):
+		print lower,upper,bracketed
+		args.force = .5*(upper+lower)
+		res = run_simulation(args); i += 1
+		if res['fundamental_frequency']>0:
+			freqs.append(res['fundamental_frequency'])
+			forces.append(args.force)
+			lower = .5*(upper+lower)
+		else:
+			upper = .5*(upper+lower)
+	return forces,freqs,res
+'''