cantbeam_guyan.m

MATLAB PROGRAMS FOR "VIBRATION SIMULATION USING MATLAB AND ANSYS" All the M-files which are liste

	echo off
%	cantbeam_guyan.m	cantilever beam finite element program, 
%	selectable number of elements.  Solves for eigenvalues and 
%	eigenvectors of a cantilever with user-defined dimensions, 
%	material properties, number of elements and number of mode shapes 
%	to plot.  Guyan reduction is an option.  A 10 element beam is used 
%	as an example.  Default beam is 2mm wide by 20mm long by 0.2mm thick.

	clf;

	clear all;

	inp = input('Input "1" to enter beam dimensions, "Enter" to use default ...  ');

	if (isempty(inp))
   		inp = 0;
	else
	end

	if  inp == 0

   		wbeam = 2.0
   		tbeam = 0.2
   		lbeam = 20.0
   		E = 190e6
   		density = 7.83e-6

	else

%	input size of beam and material

	wbeam = input('Input width of beam, default 2mm, ...  ');

	if (isempty(wbeam))
   		wbeam = 2.0;
	else
	end

	tbeam = input('Input thickness of beam, default 0.2mm, ...  ');

	if (isempty(tbeam))
   		tbeam = 0.2;
	else
	end

	lbeam = input('Input length of beam, default 20mm, ...  ');

	if (isempty(lbeam))
   		lbeam = 20.0;
	else
	end

	E = input('Input modulus of material, mN/mm^2, default stainless steel 190e6 ...  ');
	
	if (isempty(E))
   		E = 190e6;
	else
	end

	density = input('Input density of material, Kg/mm^3, default stainless steel 7.83e-6 ...  ');

	if (isempty(density))
   		density = 7.83e-6;
	else
	end

	end

%	input number of elements
	
	num_elements = input('Input number of elements for beam, minimum 2, default 10 ...  ');

	if (isempty(num_elements))
   		num_elements = 10;
	else
	end

%	define whether or not to do Guyan Reduction

	guyan = input('enter "1" to do Guyan elimination of rotations, "enter" to not do Guyan ...  ');

	if (isempty(guyan))
   		guyan = 0;
	else
	end

	if  guyan == 0

		num_plot_max = 2*num_elements;

		num_plot_default = num_elements;

	else

		num_plot_max = num_elements;

		num_plot_default = num_elements;

	end

	num_plot = input(['enter the number of modes to plot, max ',num2str(num_plot_max),', default ',num2str(num_plot_default),' ...  ']);

	if (isempty(num_plot))
   		num_plot = 9;
	else
	end

%	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*mnu);

	m = zeros(2*mnu);

%	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) = [];

	if  guyan == 1

%	Guyan Reduction - reduce out the rotation dof's, leaving displacement dof's
%	re-order the matrices

%	re-order the columns of k

	kr = zeros(2*(mnu-1));

	krr = zeros(2*(mnu-1));

%   rearrange columns, rotation and then displacement dof's

	mkrcolcnt = 0;

	for  mkcolcnt = 2:2:2*(mnu-1)

		mkrcolcnt = mkrcolcnt + 1;

		kr(:,mkrcolcnt) = k(:,mkcolcnt);

		mr(:,mkrcolcnt) = m(:,mkcolcnt);

	end

	mkrcolcnt = num_elements;

	for  mkcolcnt = 1:2:2*(mnu-1)

		mkrcolcnt = mkrcolcnt + 1;

		kr(:,mkrcolcnt) = k(:,mkcolcnt);

		mr(:,mkrcolcnt) = m(:,mkcolcnt);

	end

%   rearrange rows, rotation and then displacement dof's

	mkrrowcnt = 0;

	for  mkrowcnt = 2:2:2*(mnu-1)

		mkrrowcnt = mkrrowcnt + 1;

		krr(mkrrowcnt,:) = kr(mkrowcnt,:);

		mrr(mkrrowcnt,:) = mr(mkrowcnt,:);

	end

	mkrrowcnt = num_elements;

	for  mkrowcnt = 1:2:2*(mnu-1)

		mkrrowcnt = mkrrowcnt + 1;

		krr(mkrrowcnt,:) = kr(mkrowcnt,:);

		mrr(mkrrowcnt,:) = mr(mkrowcnt,:);

	end

%	define sub-matrices and transformation matrix T

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

	kab = krr(1:num_elements,num_elements+1:2*num_elements);

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

%	calculate reduced mass and stiffness matrices

	kbb = T'*krr*T
	
	mbb = T'*mrr*T

	else

	kbb = k;

	mbb = m;

	end

%	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;

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

%	divide the negative of the stiffness matrix by the mass matrix

	ksm = inv(mbb)*(-kbb);

%	now expand to state space size
%	fill out unit values in mass and stiffness matrices

	for  row = 1:2:ssdof

		aud(row,row+1) = 1;

	end

%	fill out mass and stiffness terms from m and k

	for  row = 2:2:ssdof

		for  col = 2:2:ssdof

			aud(row,col-1) = ksm(row/2,col/2);

		end

	end

%	calculate the eigenvalues/eigenvectors of the undamped matrix for plotting
%	and for calculating the damping matrix c

	[evec1,evalu] = eig(aud);

	evalud = diag(evalu);

	evaludhz = evalud/(2*pi);

	num_modes = length(evalud)/2;

%	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

%	now check for any imaginary eigenvectors and convert to real

	for  cnt = 1:length(evalud)

		if  (imag(evec(1,cnt)) & imag(evec(3,cnt)) & imag(evec(5,cnt))) ~= 0

			evec(:,cnt) = imag(evec(:,cnt));

		else

		end

	end

	if  guyan == 0

%	now separate the displacement and rotations in the eigenvectors
%	for plotting mode shapes

	evec_disp = zeros(ceil(dof/2),ssdof);
	
	rownew = 0;

	for  row = 1:4:ssdof

		rownew = rownew+1;

		evec_disp(rownew,:) = evec(row,:);
	
	end

	evec_rotation = zeros(ceil(dof/2),ssdof);

	rownew = 0;
	
	for  row = 3:4:ssdof
		
		rownew = rownew+1;

		evec_rotation(rownew,:) = evec(row,:);

	end

	else

	evec_disp = zeros(ceil(dof/4),ssdof);
	
	rownew = 0;

	for  row = 1:2:ssdof

		rownew = rownew+1;

		evec_disp(rownew,:) = evec(row,:);
	
	end

	end

%	normalize the displacement eigenvectors wrt one for plotting

	for  col = 1:ssdof

		evec_disp(:,col) = evec_disp(:,col)/max(abs(real(evec_disp(:,col))));

		if  evec_disp(floor(dof/2),col) >= 0

			evec_disp(:,col) = -evec_disp(:,col);

		else
		end

	end

%	list eigenvalues, hz

	format long	e

	evaludhz_list = sort(evaludhz(1:2:2*num_modes))

	format short

%	list displacement (not velocity) eigenvectors

	evec_disp(:,1:2:2*num_plot)

	if  guyan == 0

%	plot mode shapes

  	for  mode_cnt = 1:num_plot

  		evec_cnt = 2*mode_cnt -1;

  		plot(lvec,[0; evec_disp(:,evec_cnt)],'ko-')
  		title(['Cantilever Beam, Mode ', ...
  			num2str(mode_cnt),': ',num2str(abs(evalhzr(evec_cnt))),' hz']);
  		xlabel('Distance From Built-In End')
  		ylabel('Normalized Y-Displacement')
  		axis([0 lbeam -1.5 1.5])
  		grid on
	
	disp('execution paused to display figure, "enter" to continue'); pause

  	end

	else

%	plot mode shapes, Guyan Reduced

  	for  mode_cnt = 1:num_plot

  		evec_cnt = 2*mode_cnt -1;

  		plot(lvec,[0; evec_disp(:,evec_cnt)],'ko-')
  		title(['Cantilever Beam, Mode ', ...
  			num2str(mode_cnt),': ',num2str(abs(evalhzr(evec_cnt))),' hz']);
  		xlabel('Distance From Built-In End')
  		ylabel('Normalized Y-Displacement')
  		axis([0 lbeam -1.5 1.5])
  		grid on
	
	disp('execution paused to display figure, "enter" to continue'); pause

  	end

	end

%	normalization with respect to mass on a filled (not diagonal) mass matrix

%	calculate the displacement (displacement and rotation) eigenvectors
%	to be used for the modal model eigenvectors

	xm = zeros(dof);
	
	col = 0;

	for  mode = 1:2:ssdof

		col = col + 1;

		row = 0;

		for	 ndof = 1:2:ssdof

			row = row + 1;

			xm(row,col) = evec(ndof,mode);

		end

	end

%	normalize with respect to mass

	for  mode = 1:dof

		xn(:,mode) = xm(:,mode)/sqrt(xm(:,mode)'*mbb*xm(:,mode));

	end

%	calculate the normalized mass and stiffness matrices for checking

	mm = xn'*mbb*xn;

	km = xn'*kbb*xn;

%	check that the sqrt of diagonal elements of km are eigenvalues

	p = (diag(km)).^0.5;

	row = 0;

	for cnt = 1:2:ssdof

		row = row + 1;

		evalrad(row) = abs((eval(cnt)));

	end

	[p evalrad']/(2*pi)

	evalhz = evalrad/(2*pi);

	semilogy(evalhz)
	title('Resonant Frequencies, Hz')
	xlabel('Mode Number')
	ylabel('Frequency, hz')
	grid
	disp('execution paused to display figure, "enter" to continue'); pause