cant_2el_guyan.m

振动仿真 vibration simulation using matlab and ansys!

	echo off
%	cant_2el_guyan.m	two element cantilever beam finite element 
%	program using guyan reduction to eliminate rotations, both Matlab 
%	and hand-calculation results are listed at each step of the program.

%	beam is 2mm wide by 20mm long by 0.2mm thick

	clear all;

%	input size of beam and material

	wbeam = 2.0;

	tbeam = 0.2;

	lbeam = 20.0;

	E = 190e6;

	density = 7.83e-6;

	num_elements = 2;

%	define length of each element, uniform lengths

	l = lbeam/num_elements;

%	define length vector for plotting, right-to-left numbering

	lvec = 0:l:lbeam;

%	define the node numbers

	n = 1:num_elements+1;

%	number the nodes for the elements

	node1 = 1:num_elements;
	
	node2 = 2:num_elements+1;

%	size the stiffness and mass matrices to have 2 times the number of nodes
%	to allow for translation and rotation dof's for each node, including built-
%	in end

	max_node1 = max(node1);

	max_node2 = max(node2);

	max_node_used = max([max_node1 max_node2]);

	mnu = max_node_used;

	k = zeros(2*max_node_used);

	m = zeros(2*max_node_used);

%	now build up the global stiffness and consistent mass matrices, element by element	

%	calculate I, area and mass per unit length of beam

	I = wbeam*tbeam^3/12;

	area = wbeam*tbeam;

	mpl = density*area;

	for  i = 1:num_elements

		dof1 = 2*node1(i)-1;
		dof2 = 2*node1(i);
		dof3 = 2*node2(i)-1;
		dof4 = 2*node2(i);

		k(dof1,dof1) = k(dof1,dof1)+(12*E*I/l^3);
		k(dof2,dof1) = k(dof2,dof1)+(6*E*I/l^2);
		k(dof3,dof1) = k(dof3,dof1)+(-12*E*I/l^3);
		k(dof4,dof1) = k(dof4,dof1)+(6*E*I/l^2);

		k(dof1,dof2) = k(dof1,dof2)+(6*E*I/l^2);
		k(dof2,dof2) = k(dof2,dof2)+(4*E*I/l);
		k(dof3,dof2) = k(dof3,dof2)+(-6*E*I/l^2);
		k(dof4,dof2) = k(dof4,dof2)+(2*E*I/l);
	
		k(dof1,dof3) = k(dof1,dof3)+(-12*E*I/l^3);
		k(dof2,dof3) = k(dof2,dof3)+(-6*E*I/l^2);
		k(dof3,dof3) = k(dof3,dof3)+(12*E*I/l^3);
		k(dof4,dof3) = k(dof4,dof3)+(-6*E*I/l^2);
	
		k(dof1,dof4) = k(dof1,dof4)+(6*E*I/l^2);
		k(dof2,dof4) = k(dof2,dof4)+(2*E*I/l);
		k(dof3,dof4) = k(dof3,dof4)+(-6*E*I/l^2);
		k(dof4,dof4) = k(dof4,dof4)+(4*E*I/l);
	
		m(dof1,dof1) = m(dof1,dof1)+(mpl/420)*(156*l);
		m(dof2,dof1) = m(dof2,dof1)+(mpl/420)*(22*l^2);
		m(dof3,dof1) = m(dof3,dof1)+(mpl/420)*(54*l);
		m(dof4,dof1) = m(dof4,dof1)+(mpl/420)*(-13*l^2);
	
		m(dof1,dof2) = m(dof1,dof2)+(mpl/420)*(22*l^2);
		m(dof2,dof2) = m(dof2,dof2)+(mpl/420)*(4*l^3);
		m(dof3,dof2) = m(dof3,dof2)+(mpl/420)*(13*l^2);
		m(dof4,dof2) = m(dof4,dof2)+(mpl/420)*(-3*l^3);
	
		m(dof1,dof3) = m(dof1,dof3)+(mpl/420)*(54*l);
		m(dof2,dof3) = m(dof2,dof3)+(mpl/420)*(13*l^2);
		m(dof3,dof3) = m(dof3,dof3)+(mpl/420)*(156*l);
		m(dof4,dof3) = m(dof4,dof3)+(mpl/420)*(-22*l^2);
	
		m(dof1,dof4) = m(dof1,dof4)+(mpl/420)*(-13*l^2);
		m(dof2,dof4) = m(dof2,dof4)+(mpl/420)*(-3*l^3);
		m(dof3,dof4) = m(dof3,dof4)+(mpl/420)*(-22*l^2);
		m(dof4,dof4) = m(dof4,dof4)+(mpl/420)*(4*l^3);
		
	end

%	now that stiffness and mass matrices are defined for all dof's, including
%	constrained dof's, need to delete rows and columns of the matrices that
%	correspond to constrained dof's, in the left-to-right case, the first two
%	rows and columns

	k(1:2,:) = [];	% translation/rotation of node 1
	k(:,1:2) = [];

	m(1:2,:) = [];
	m(:,1:2) = [];

%	compare with hand calculation

	k

	khand = E*I*[24/l^3		0		-12/l^3		6/l^2
				    0		8/l		-6/l^2		2/l
				-12/l^3    -6/l^2    12/l^3    -6/l^2
				6/l^2		2/l		-6/l^2		4/l]

	m

	mhand = (mpl/420)*[312*l		0			54*l		-13*l^2
					    0		  	8*l^3		13*l^2		-3*l^3
						54*l		13*l^2		156*l		-22*l^2
						-13*l^2		-3*l^3		-22*l^2		4*l^3]

%	Guyan Reduction - reduce out the rotation dof's, leaving displacement dof's
%	re-order the matrices, rotations in row/column 2 and 4 become 1 and 2 and
%   translations in row/column 1 and 3 become 3 and 4

%	re-order the k columns

	kr = [k(:,2) k(:,4) k(:,1) k(:,3)];

%	re-order the rows

	krr = [kr(2,:)
		   kr(4,:)
		   kr(1,:)
		   kr(3,:)]

	krrhand = E*I*[8/l		2/l		0		-6/l^2
				   2/l		4/l		6/l^2	-6/l^2
				   0		6/l^2	24/l^3	-12/l^3
				  -6/l^2   -6/l^2  -12/l^3   12/l^3]		

	krr - krrhand

%	re-order the m columns

	mr = [m(:,2) m(:,4) m(:,1) m(:,3)];

%	re-order the rows

	mrr = [mr(2,:)
		   mr(4,:)
		   mr(1,:)
		   mr(3,:)]

	mrrhand = (mpl/420)*[8*l^3		-3*l^3			0		13*l^2
						-3*l^3		4*l^3		-13*l^2	   -22*l^2
						0			-13*l^2		  312*l		54*l
						13*l^2		-22*l^2		  54*l		156*l]
				
%	define sub-matrices and transformation matrix T

	kaa = krr(1:2,1:2);

	kab = krr(1:2,3:4);

	T = [-inv(kaa)*kab
			eye(2,2)]

	Thand = [6/(14*l)	6/(14*l)
		   -24/(14*l)  18/(14*l)
		        1		   0
				0		   1]
				
%	calculate reduced mass and stiffness matrices

	kbb = T'*krr*T
	
	kbbhand = (E*I/(14*l^3))*[192   -60
							  -60    24]

	mbb = T'*mrr*T

	mbbhand = mpl*l*[1528/1715  	241/1372
					  241/1372		471/1715]

%	now convert to state-space form for Matlab calculations

%	define the number of dof for state-space version, 2 times dof left after
%	removing constrained dof's

	[dof,dof] = size(kbb);

%	define the sizes of mass and stiffness matrices for state-space

	ssdof = 2*dof;

	kss = zeros(ssdof);		% creates a ssdof x ssdof null matrix

	mss = zeros(ssdof);

%	fill out unit values in mass and stiffness matrices

	for  row = 1:2:ssdof

		kss(row,row+1) = -1;

		mss(row,row) = 1;

	end

%	fill out mass and stiffness terms from m and k

	for  row = 2:2:ssdof

		for  col = 2:2:ssdof

			kss(row,col-1) = kbb(row/2,col/2);

			mss(row,col) = mbb(row/2,col/2);

		end

	end

	kss

	ksshand1= (E*I/(14*l^3))*[0		0		0		0
							 192    0      -60      0
							  0		0		0		0
							 -60    0       24      0];

	ksshand = ksshand1 + [0 -1 0 0
						  0 0 0 0
						  0 0 0 -1
						  0 0 0 0]

	mss

	msshand1 = mpl*l*[0		0		0		0
			          0 1528/1715	0    241/1372
				      0		0		0		0
				      0  241/1372   0    471/1715];

	msshand = msshand1 + [1 0 0 0
						  0 0 0 0
						  0 0 1 0
						  0 0 0 0]

%	take the inverse of mss

	mssinv = inv(mss)

	mssinvhand = [1 				0					0					0
				  0     263760/(205367*mpl*l)			0     -168700/(205367*mpl*l)
				  0					0					1					0
				  0    -168700/(205367*mpl*l)			0      855680/(205367*mpl*l)]

%	premultipy -kss by the inverse of mss to get the system matrix

	a = mssinv*(-kss)

	ahand = [		0			     1				0					0
	  -4340280*E*I/(205367*mpl*l^4)	 0  1419600*E*I/(205367*mpl*l^4)	0
					0				 0				0			   		1
	   5980800*E*I/(205367*mpl*l^4)  0 -2189800*E*I/(205367*mpl*l^4)  	0]
	   
%	calculate the eigenvalues/eigenvectors of the undamped matrix for plotting
%	and for calculating the damping matrix c

	[evec1,evalu] = eig(a);

	evalud = diag(evalu);

	evaludhz = evalud/(2*pi);

%	now reorder the eigenvalues and eigenvectors from low to high freq	

	[evalorder,indexhz] = sort(abs(evalud));

	for  cnt = 1:length(evalud)

		eval(cnt,1) = evalud(indexhz(cnt));

		evalhzr(cnt,1) = round(evaludhz(indexhz(cnt)));

		evec(:,cnt) = evec1(:,indexhz(cnt));

	end

%	hand calculation of resonant frequencies

	f1 = (1/(2*pi))*(2/205367)*(1/(mpl*l^2))*sqrt(43127070)*sqrt(I*mpl*E*(3887-20*sqrt(34178)));

	f2 = (1/(2*pi))*(2/205367)*(1/(mpl*l^2))*sqrt(43127070)*sqrt(I*mpl*E*(3887+20*sqrt(34178)));

	eig_matlab = sort(evaludhz(1:2:4))
	
	eighand = [f1; f2]