Example:Timoshenko Locking

From Nemesis
Jump to: navigation, search
u=[]
r=[2,3,4,5,10,20,50,100,500,1000]
for beam in ['EULER','TIM2F','TIM2R','TIM3F','TIM3R']:
	t=[]
	for ratio in r:
		print 'Solving for type %s and r=%4i.'%(beam,ratio)
		domain.dim(2)
		material.uniElastic(1,1000,0.)
		section.rect(1,1.0,20./ratio)

		# nodes
		node.add(1,0.,0.)
		for i in range(20): 
			node.add(2*i+2,i+0.5, 0.)
			node.add(2*i+3,i+1.0, 0.)
		
		# beams
		for i in range(20):
			if beam=='EULER':
				element.beam2e(i+1,2*i+1,2*i+3,1,1)
			elif beam=='TIM2F':
				element.beam2t(i+1,2*i+1,2*i+3,1,1,2)
			elif beam=='TIM2R':
				element.beam2t(i+1,2*i+1,2*i+3,1,1,1)
			elif beam=='TIM3F':
				element.beam3t(i+1,2*i+1,2*i+3,2*i+2,1,1,3)
			elif beam=='TIM3R':
				element.beam3t(i+1,2*i+1,2*i+3,2*i+2,1,1,3)
		
		node.fix(1,1)
		node.fix(1,2)
		node.fix(1,6)

		lc.define(1)
		load.node(21,2,1.0)

		analysis.static()
		analysis.run(1,1)

		t.append(node.data(41)['disp'][1])
		domain.clear()
	u.append(t)

# Results
for i in range(1,5):
	for j in range(len(r)):
		u[i][j]=u[i][j]/u[0][j]
for j in range(len(r)):
	u[0][j]=1.0
print 
print "====================================================="
print "|Ratio    |Tip Displacement (normalized)            |"
print "-----------------------------------------------------"
print "|         |   EULER |TIM2F | TIM2R | TIM3F |  TIM3R |"
print "====================================================="
for i in range(len(r)):
	print '|%9.3f|'%(r[i]),
	for j in range(5):
		print ' %6.4f'%(u[j][i]),
	print '|'
print "====================================================="


import pylab
from numpy import array
x=array(r)
y1=array(u[0])
y2=array(u[1])
y3=array(u[2])
y4=array(u[3])
y5=array(u[4])
pylab.semilogx(x,y1,'k',x,y2,'r:o',x,y3,'r-',x,y4,'g:o',x,y5,'g-')
pylab.legend(('Euler-Bernulli', 'Timoshenko 2 nodes, 2 Gauss-points', 'Timoshenko 2 nodes, 1 Gauss-point','Timoshenko 3 nodes, 3 Gauss-points', 'Timoshenko 3 nodes, 2 Gauss-points'),
       'upper right', shadow=True)
pylab.xlabel('log (L/h)')
pylab.ylabel('w/wEB')
pylab.grid(True)
pylab.show()
Personal tools