Newer
Older
#!/usr/bin/env python
from __future__ import division
from numpy import *
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!
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure()
cmap = plt.cm.get_cmap('Dark2', len(beamsets)+1)
ax = fig.gca(projection='3d')
for i,beamset in enumerate(beamsets):
for seg in beamset:
ax.plot(nodes[seg,0], nodes[seg,1], nodes[seg,2], c=cmap(i))
plt.show()
def get_rotation_from_nodes(nodes,axis,disps,targets):
base = asarray(axis[0]);
a = axis[1]-base;
a = a/magnitudes(a) #set up axis coordinates
x = nodes[targets]-base
d = x - dot(x, a)[...,None]*a
y = [dot(a,b) for a,b in zip(disps[targets,:3], cross(d,a))]
return arctan2(y,magnitudes(d)),magnitudes(d)
def run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads):
write_frame3dd_file(nodes,global_args,beam_sets,constraints,loads)
cmd = ["frame3dd", "-i",global_args['frame3dd_filename']+'.csv']
if args.quiet: cmd.append("-q")
print ' '.join(cmd)
subprocess.call(cmd)
def clean_up_frame3dd(filename):
#Delete files generated by frame3dd
files = [filename+end for end in ["_out.csv",".csv.out",".csv.plt",".csv.if01",".csv"]]
subprocess.call(["rm"]+files)
def build(args):
#return nodes,rods as numpy arrays
dxy = args.attach_radius/sqrt(2)
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] ])
nodes = vstack(( nodes+z_os, nodes-z_os ))
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] ]) ))
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
if args.flexure_type == 'cyclic':
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, flexure_nodes))
return nodes, flexure_beams, solid_beams, fixed_nodes
def run_simulation(args):
#set up simulation
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]),
'frame3dd_filename':args.base_filename+"_frame3dd"
}
clean_up_frame3dd(global_args['frame3dd_filename'])
beam_sets = [
(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':[]})
]
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 = range(8)
sensor_nodes = [8,9,10,11]
for force_dof in [0,1,2]:
loads = [{'node':n,'DOF':force_dof,'value':args.force/len(loaded_nodes)} for n in loaded_nodes]
run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads)
disps = read_frame3dd_displacements(global_args['frame3dd_filename'])
force_disp = average(disps[loaded_nodes,force_dof])
#print "Degree of freedom: %d"%force_dof
#print "Force applied: %.1f N"%(args.force)
#print "Displacement at sensor: %.1f microns"%(force_disp*1e6)
results.append( {'dof':force_dof, 'force/torque':args.force, 'displacement':force_disp} )
#sys.exit(0)
#Tx, Ty
torque_force = args.torque/(.5*args.sep)/len(loaded_nodes) # F = torque / dist
for torque_dof in [0,1]:
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 = [8,10] if torque_dof==0 else [9,11]
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#print disps[sensor_nodes]
#axis = (array([0,0,0]), array([0,-1,0]) if torque_dof==0 else array([1,0,0]) )
#rots,ds = get_rotation_from_nodes(nodes,axis,disps/global_args['length_scaling'],loaded_nodes)
#print average(rots)
#print "Degree of freedom: %d"%torque_dof
#print "Torque applied: %.1f Nm"%(args.torque)
#print "Torque force: %.2f N"%torque_force
#print "Displacement at sensor: %.1f microns"%(sin(average(rots))*args.sensor_radius*1e6)
#print "Displacement at sensor: %.1f microns"%( average(abs(disps[moving_sensor_nodes,2]))*1e6 )
results.append( {'dof':torque_dof+3, 'force/torque':args.torque, 'displacement':average(abs(disps[moving_sensor_nodes,2]))} )
#Tz
torque_force = args.torque/args.attach_radius/len(loaded_nodes) # F = torque / dist
def node_to_force(n):
return array([-nodes[n,1],nodes[n,0],0]) / sqrt(nodes[n,0]**2 + nodes[n,1]**2)
loads = [{'node':n,'DOF':i,'value':torque_force*node_to_force(n)[i]} for n in loaded_nodes for i in [0,1]]
run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads)
disps = read_frame3dd_displacements(global_args['frame3dd_filename'])
#axis = (array([0,0,0]), array([0,0,1]))
#rots,ds = get_rotation_from_nodes(nodes,axis,disps/global_args['length_scaling'],loaded_nodes)
#print average(rots)
#print "Degree of freedom: %d"%2
#print "Torque applied: %.1f Nm"%(args.torque)
#print "Torque force: %.2f N"%torque_force
#print "Displacement at sensor: %.1f microns"%(sin(average(rots))*args.sensor_radius*1e6)
#print "Displacement at sensor: %.1f microns"%( average(magnitudes(disps[sensor_nodes,:3]))*1e6 )
results.append( {'dof':5, 'force/torque':args.torque, 'displacement':average(magnitudes(disps[sensor_nodes,:3]))} )
return results
if __name__ == '__main__':
parser = argparse.ArgumentParser()
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)")
parser.add_argument("-torque","--torque", type=double, default=1., help="torque to apply (Nm)")
#parser.add_argument("-fr","--force_res", type=double, default=.01, help="Final resolution of force for search mode")
parser.add_argument("-w","--w", type=double, default=.0005, help="width of flexure (m)")
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=.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")
parser.add_argument("-nu","--nu", type=double, default=.33, help="Poisson Ratio")
parser.add_argument("-base_filename","--base_filename", default='buckle', help="Base filename for segments and frame3dd")
parser.add_argument("-rho","--rho",type=double,default=2700.,help='density of beam material, kg/m^3')
parser.add_argument("-n_modes","--n_modes",type=int,default=0,help='number of dynamic modes to compute')
parser.add_argument("-ls","--length_scaling", type=double, default=1., help="Scale factor to keep numbers commesurate")
args = parser.parse_args()
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
#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 "Fundamental frequency: %.3f Hz"%res['fundamental_frequency']
#print "Stress: %.3f MPa"%(res['stress']/1e6)
elif args.mode=='visualize':
nodes,rods,solid_beams,fixed_nodes = build(args)
plot_connections(nodes,[rods,solid_beams])
else:
assert(0) #should not be here
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
'''
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
'''