tdof_non_prop_damped.m

vibration simulation using ansys and matlab 一书中的程序

	echo off
%	tdof_non_prop_damped.m	non-proportionally damped tdof model
%	This code is used to develop an understanding of the results 
%	of Matlab's eigenvalue analysis and complex modes.

	clf;

	legend off;

	subplot(1,1,1);

	clear all;

%	define the values of masses, springs, dampers

	m1 = 1;
	m2 = 1;
	m3 = 1;

	k1 = 1;
	k2 = 1;

%	define arbitrary damping values

   	c1 = input('input value for c1, default 0.1, ...  ');
	
	if  (isempty(c1))
		c1 = 0.1;
	else
	end

	c2 = input('input value for c1, default 0.2, ...  ');
	
	if  (isempty(c2))
		c2 = 0.2;
	else
	end

%	define the system matrix, aphys, in physical coordinates

	aphys = [  0      1         0			 0		   0	  0
     		-k1/m1 -c1/m1     k1/m1        c1/m1       0	  0	
	     	   0      0         0            1         0      0
	     	 k1/m2  c1/m2  -(k1+k2)/m2  -(c1+c2)/m2  k2/m2  c2/m2
		       0      0         0            0         0      1
		       0      0       k2/m3        c2/m3    -k2/m3 -c2/m3];

%	solve for the eigenvalues of the system matrix

	[xm,lambda] = eig(aphys)

%	take the diagonal elements of the generalized eigenvalue matrix lambda

	lambdad = diag(lambda)

%	now reorder the eigenvalues and eigenvectors from low to high frequency, 
%	keeping track of how the eigenvalues are ordered in reorder the 
%	eigenvectors to match, using indexhz

	[lambdaorder,indexhz] = sort(abs(imag(lambdad)));

	for  cnt = 1:length(lambdad)

		lambdao(cnt,1) = lambdad(indexhz(cnt));	%  reorder eigenvalues

		xmo(:,cnt) = xm(:,indexhz(cnt));		%  reorder eigenvector columns

	end

%	now normalize the eigenvectors with respect to the position of mass 1, which
%	will be set to 1.0

	for  cnt = 1:length(lambdad)

		xmon1(:,cnt) = xmo(:,cnt)/xmo(1,cnt);

	end

%	now calculate the magnitude and phase angle of each of the eigenvector
%	entries

	for  row = 1:length(lambdad)

		for  col = 1:length(lambdad)

			xmon1mag(row,col) = abs(xmon1(row,col));

			xmon1ang(row,col) = (180/pi)*angle(xmon1(row,col));

		end

	end

	lambdao

	xmo

	xmon1

	xmon1mag

	xmon1ang

%	calculate the percentage of critical damping for each mode

	zeta1 = 0
	
	theta2 = atan(real(lambdao(3))/imag(lambdao(3)));
	zeta2 = abs(sin(theta2))

	theta3 = atan(real(lambdao(5))/imag(lambdao(5)));
	zeta3 = abs(sin(theta3))

	plot(lambda,'k*')
	grid on
	axis([-3 1 -2 2])
	axis('square')
	title('Damped Eigenvalues')
	xlabel('real')
	ylabel('imaginary')
	text(real(lambdao(3))-1.2,imag(lambdao(3))+0.1,['zeta = ',num2str(zeta2)])	
	text(real(lambdao(5))-1.2,imag(lambdao(5))+0.1,['zeta = ',num2str(zeta3)])	
	disp('execution paused to display figure, "enter" to continue'); pause

%	calculate the motions of the three masses for all three modes - damped case

	t = 0:.12:15;

	sigma11 = real(lambdao(1));		%  sigma for first eigenvalue for mode 1
	omega11 = imag(lambdao(1));		%  omega for first eigenvalue for mode 1

	sigma12 = real(lambdao(2));		%  sigma for second eigenvalue for mode 1
	omega12 = imag(lambdao(2));		%  omega for second eigenvalue for mode 1

	sigma21 = real(lambdao(3));		%  sigma for first eigenvalue for mode 2
	omega21 = imag(lambdao(3));		%  omega for first eigenvalue for mode 2

	sigma22 = real(lambdao(4));		%  sigma for second eigenvalue for mode 2
	omega22 = imag(lambdao(4));		%  omega for second eigenvalue for mode 2

	sigma31 = real(lambdao(5));		%  sigma for first eigenvalue for mode 3
	omega31 = imag(lambdao(5));		%  omega for first eigenvalue for mode 3

	sigma32 = real(lambdao(6));		%  sigma for second eigenvalue for mode 3
	omega32 = imag(lambdao(6));		%  omega for second eigenvalue for mode 3

%	motion of three masses for mode 1

	z111 = exp(sigma11*t).*(exp(i*omega11*t)*xmon1(1,1));		%  mass 1
	z112 = exp(sigma12*t).*(exp(i*omega12*t)*xmon1(1,2));		%  mass 1

	z121 = exp(sigma11*t).*(exp(i*omega11*t)*xmon1(3,1));		%  mass 2
	z122 = exp(sigma12*t).*(exp(i*omega12*t)*xmon1(3,2));		%  mass 2

	z131 = exp(sigma11*t).*(exp(i*omega11*t)*xmon1(5,1));		%  mass 3
	z132 = exp(sigma12*t).*(exp(i*omega12*t)*xmon1(5,2));		%  mass 3

%	motion of three masses for mode 2

	z211 = exp(sigma21*t).*(exp(i*omega21*t)*xmon1(1,3));		%  mass 1
	z212 = exp(sigma22*t).*(exp(i*omega22*t)*xmon1(1,4));		%  mass 1

	z221 = exp(sigma21*t).*(exp(i*omega21*t)*xmon1(3,3));		%  mass 2
	z222 = exp(sigma22*t).*(exp(i*omega22*t)*xmon1(3,4));		%  mass 2

	z231 = exp(sigma21*t).*(exp(i*omega21*t)*xmon1(5,3));		%  mass 3
	z232 = exp(sigma22*t).*(exp(i*omega22*t)*xmon1(5,4));		%  mass 3

%	motion of three masses for mode 3

	z311 = exp(sigma31*t).*(exp(i*omega31*t)*xmon1(1,5));		%  mass 1
	z312 = exp(sigma32*t).*(exp(i*omega32*t)*xmon1(1,6));		%  mass 1

	z321 = exp(sigma31*t).*(exp(i*omega31*t)*xmon1(3,5));		%  mass 2
	z322 = exp(sigma32*t).*(exp(i*omega32*t)*xmon1(3,6));		%  mass 2

	z331 = exp(sigma31*t).*(exp(i*omega31*t)*xmon1(5,5));		%  mass 3
	z332 = exp(sigma32*t).*(exp(i*omega32*t)*xmon1(5,6));		%  mass 3

%	plot real and imaginary motions of each mass for the two complex conjugate
%   eigenvectors of mode 2

	plot(t,real(z211),'k-',t,real(z212),'k+-',t,imag(z211),'k.-',t,imag(z212),'ko-')
	title('non-prop damped real and imag for z1, mode 2')
	legend('real','real','imag','imag')
	xlabel('time, sec')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z221),'k-',t,real(z222),'k+-',t,imag(z221),'k.-',t,imag(z222),'ko-')
	title('non-prop damped real and imag for z2 mode 2')
	legend('real','real','imag','imag')
	xlabel('time, sec')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z231),'k-',t,real(z232),'k+-',t,imag(z231),'k.-',t,imag(z232),'ko-')
	title('non-prop damped real and imag for z3 mode 2')
	legend('real','real','imag','imag')
	xlabel('time, sec')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z211+z212),'k-',t,real(z221+z222),'k+-',t,real(z231+z232),'k.-')
	title('non-prop damped, z1, z2, z3 mode 2')
	legend('mass 1','mass 2','mass 3')
	xlabel('time, sec')
	axis([0 max(t) -2 2])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

%	plot subplots for notes

	subplot(2,2,1)
	plot(t,real(z211),'k-',t,real(z212),'k+',t,imag(z211),'k.-',t,imag(z212),'ko-')
	title('non-prop damped real and imag for z1, mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -1 1])
	grid on

	subplot(2,2,2)
	plot(t,real(z221),'k-',t,real(z222),'k+',t,imag(z221),'k.-',t,imag(z222),'ko-')
	title('non-prop damped real and imag for z2 mode 2')
	legend('real','real','imag','imag')
	axis([0 max(t) -1 1])
	grid on

	subplot(2,2,3)
	plot(t,real(z231),'k-',t,real(z232),'k+',t,imag(z231),'k.-',t,imag(z232),'ko-')
	title('non-prop damped real and imag for z3 mode 2')
	legend('real','real','imag','imag')
	xlabel('time, sec')
	axis([0 max(t) -1 1])
	grid on

	subplot(2,2,4)
	plot(t,real(z211+z212),'k-',t,real(z221+z222),'k+-',t,real(z231+z232),'k.-')
	title('non-prop damped, z1, z2, z3 mode 2')
	legend('mass 1','mass 2','mass 3')
	grid on
	xlabel('time, sec')
	axis([0 max(t) -2 2])

	disp('execution paused to display figure, "enter" to continue'); pause

	subplot(1,1,1)

%	plot mode 3	

	plot(t,real(z311),'k-',t,real(z312),'k+-',t,imag(z311),'k.-',t,imag(z312),'ko-')
	title('non-prop damped real and imag for z1, mode 3')
	legend('real','real','imag','imag')
	xlabel('time, sec')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z321),'k-',t,real(z322),'k+-',t,imag(z321),'k.-',t,imag(z322),'ko-')
	title('non-prop damped real and imag for z2 mode 3')
	legend('real','real','imag','imag')
	xlabel('time, sec')
	axis([0 max(t) -2 2])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z331),'k-',t,real(z332),'k+-',t,imag(z331),'k.-',t,imag(z332),'ko-')
	title('non-prop damped real and imag for z3 mode 3')
	legend('real','real','imag','imag')
	xlabel('time, sec')
	axis([0 max(t) -1 1])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

	plot(t,real(z311+z312),'k-',t,real(z321+z322),'k+-',t,real(z331+z332),'k.-')
	title('non-prop damped, z1, z2, z3 mode 3')
	legend('mass 1','mass 2','mass 3')
	xlabel('time, sec')
	axis([0 max(t) -4 4])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause

%	plot subplots for notes

	subplot(2,2,1)
	plot(t,real(z311),'k-',t,real(z312),'k+-',t,imag(z311),'k.-',t,imag(z312),'ko-')
	title('non-prop damped real and imag for z1, mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -1 1])
	grid on

	subplot(2,2,2)
	plot(t,real(z321),'k-',t,real(z322),'k+-',t,imag(z321),'k.-',t,imag(z322),'ko-')
	title('non-prop damped real and imag for z2 mode 3')
	legend('real','real','imag','imag')
	axis([0 max(t) -2 2])
	grid on

	subplot(2,2,3)
	plot(t,real(z331),'k-',t,real(z332),'k+-',t,imag(z331),'k.-',t,imag(z332),'ko-')
	title('non-prop damped real and imag for z3 mode 3')
	legend('real','real','imag','imag')
	xlabel('time, sec')
	axis([0 max(t) -1 1])
	grid on

	subplot(2,2,4)
	plot(t,real(z311+z312),'k-',t,real(z321+z322),'k+-',t,real(z331+z332),'k.-')
	title('non-prop damped, z1, z2, z3 mode 3')
	legend('mass 1','mass 2','mass 3')
	xlabel('time, sec')
	axis([0 max(t) -4 4])
	grid on

	disp('execution paused to display figure, "enter" to continue'); pause