tdof_modal_xfer.m

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

	echo off
%	tdof_modal_xfer.m	plotting modal transfer functions of three dof model
%	Calculates and plots the four distinct frequency responses and
%	the individual modal contributions.	

	clf;

	legend off;

	subplot(1,1,1);

	clear all;

%	Define a vector of frequencies to use, radians/sec.  The logspace command uses
%	the log10 value as limits, i.e. -1 is 10^-1 = 0.1 rad/sec, and 1 is
%	10^1 = 10 rad/sec.  The 200 defines 200 frequency points.

	w = logspace(-1,1,150);

%	calculate the rigid-body motions for low and high frequency portions
%	of all the transfer functions

%	z11, output 1 due to force 1 transfer functions	

	z111num = 1/3;

	z111den = [1 0 0];

	z112num = 1/2;

	z112den = [1 0 1];

	z113num = 1/6;

	z113den = [1 0 3];

	[z111mag,z111phs] = bode(z111num,z111den,w);	

	[z112mag,z112phs] = bode(z112num,z112den,w);	

	[z113mag,z113phs] = bode(z113num,z113den,w);	

	if  abs(z111phs(1)) >= 10
	
		z111text = '(-)';
		
	else
	
		z111text = '(+)';
		
	end				 

	if  abs(z112phs(1)) >= 10
	
		z112text = '(-)';
		
	else
	
		z112text = '(+)';
		
	end				 

	if  abs(z113phs(1)) >= 10
	
		z113text = '(-)';
		
	else
	
		z113text = '(+)';
		
	end				 

	z111magdb = 20*log10(z111mag);

	z112magdb = 20*log10(z112mag);

	z113magdb = 20*log10(z113mag);

%	calculate the complete transfer function

	z11 = ((1/3)./((j*w).^2) + ((1/2)./((j*w).^2 + 1)) + ((1/6)./((j*w).^2 + 3)));

	z11mag = abs(z11);

	z11magdb = 20*log10(z11mag);

	z11phs = 180*angle(z11)/pi ;

%	truncate peaks for microsoft word plotting of expanded linear scale

	z11plotmag = z11mag;

	z111plotmag = z111mag;

	z112plotmag = z112mag;

	z113plotmag = z113mag;


	for  cnt = 1:length(z11mag)

		if  z11plotmag(cnt) >= 3.0
		
			z11plotmag(cnt) = 3.0;
			
		end
		
		if  z111plotmag(cnt) >= 3.0
		
			z111plotmag(cnt) = 3.0;
			
		end
		
		if  z112plotmag(cnt) >= 3.0
		
			z112plotmag(cnt) = 3.0;
			
		end
		
		if  z113plotmag(cnt) >= 3.0
		
			z113plotmag(cnt) = 3.0;
			
		end
		
	end								

%	plot the three modal contribution transfer functions and the total using
%   log magnitude versus frequency

	loglog(w,z111mag,'k+-',w,z112mag,'kx-',w,z113mag,'k.-',w,z11mag,'k-')
	title('Transfer Functions - z111, z112, z113 and z11 magnitude')
	legend('z111-1st Mode','z112-2nd Mode','z113-3rd Mode','z11-Total')
	text(.11,1.2*z111mag(1),z111text)
	text(.11,1.2*z112mag(1),z112text)
	text(.11,1.2*z113mag(1),z113text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	grid

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

%	plot the four transfer functions using db

	semilogx(w,z111magdb,'k+-',w,z112magdb,'kx-',w,z113magdb,'k.-',w,z11magdb,'k-')
	title('Transfer Function - z111, z112, z113 and z11 Magnitude')
	legend('z111-1st Mode','z112-2nd Mode','z113-3rd Mode','z11-Total')
	text(.11,2+z111magdb(1),z111text)
	text(.11,2+z112magdb(1),z112text)
	text(.11,2+z113magdb(1),z113text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude, db')
	grid

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

%	plot the four transfer functions using a linear magnitude scale so that
%	the amplitudes can be added directly

	semilogx(w,z111mag,'k+-',w,z112mag,'kx-',w,z113mag,'k.-',w,z11mag,'k-')
	title('Transfer Function - z111, z112, z113 and z11 Linear Magnitude')
	legend('z111-1st Mode','z112-2nd Mode','z113-3rd Mode','z11-Total')
	text(.11,1.0*z111mag(1),z111text)
	text(.11,1.1*z112mag(1),z112text)
	text(.11,1.1*z113mag(1),z113text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	grid

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

	semilogx(w,z111plotmag,'k+-',w,z112plotmag,'kx-',w,z113plotmag,'k.-',w,z11plotmag,'k-')
	title('Transfer Function - z111, z112, z113 and z11 Linear Magnitude')
	legend('z111-1st Mode','z112-2nd Mode','z113-3rd Mode','z11-Total')
	text(.11,1.0*z111mag(1),z111text)
	text(.11,1.1*z112mag(1),z112text)
	text(.11,1.1*z113mag(1),z113text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	axis([.1 10 0 3]);
	grid

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

%	plot phase

	semilogx(w,z111phs,'k+-',w,z112phs,'kx-',w,z113phs,'k.-',w,z11phs,'k-')
	title('Transfer Function - z111, z112, z113 and z11 Phase')
	legend('z111-1st Mode','z112-2nd Mode','z113-3rd Mode','z11-Total')
	xlabel('Frequency, rad/sec')
	ylabel('Phase, Deg')
	grid

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

	legend off;

%	z21, output 2 due to force 1 transfer functions	

	z211num = 1/3;

	z211den = [1 0 0];

%	z212 contribution = 0

	z213num = -1/3;

	z213den = [1 0 3];

	[z211mag,z211phs] = bode(z211num,z211den,w);	

	[z213mag,z213phs] = bode(z213num,z213den,w);	

	if  abs(z211phs(1)) >= 10
	
		z211text = '(-)';
		
	else
	
		z211text = '(+)';
		
	end				 

	if  abs(z213phs(1)) >= 10
	
		z213text = '(-)';
		
	else
	
		z213text = '(+)';
		
	end				 

	z211magdb = 20*log10(z211mag);

	z213magdb = 20*log10(z213mag);

%	calculate the complete transfer function

	z21 = ((1/3)./((j*w).^2) - ((1/3)./((j*w).^2 + 3)));

	z21mag = abs(z21);

	z21magdb = 20*log10(z21mag);

	z21phs = z211phs + z213phs;

%	plot the three modal contribution transfer functions and the total using
%   log magnitude versus frequency

	loglog(w,z211mag,'k+-',w,z213mag,'kx-',w,z21mag,'k-')
	title('Transfer Functions - z211, z213 and z21 magnitude')
	legend('z211-1st Mode','z213-3rd Mode','z21-Total')
	text(.11,1.2*z211mag(1),z211text)
	text(.11,1.2*z213mag(1),z213text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	grid

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

%	plot the four transfer functions using db

	semilogx(w,z211magdb,'k+-',w,z213magdb,'kx-',w,z21magdb,'k-')
	title('Transfer Function - z211, z213 and z21 Magnitude')
	legend('z211-1st Mode','z213-3rd Mode','z21-Total')
	text(.11,2+z211magdb(1),z211text)
	text(.11,2+z213magdb(1),z213text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude, db')
	grid

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

%	plot the four transfer functions using a linear magnitude scale so that
%	the amplitudes can be added directly

	semilogx(w,z211mag,'k+-',w,z213mag,'kx-',w,z21mag,'k-')
	title('Transfer Function - z211, z213 and z21 Linear Magnitude')
	legend('z211-1st Mode','z213-3rd Mode','z21-Total')
	text(.11,1.0*z211mag(1),z211text)
	text(.11,1.1*z213mag(1),z213text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	grid

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

	semilogx(w,z211mag,'k+-',w,z213mag,'kx-',w,z21mag,'k-')
	title('Transfer Function - z211, z213 and z21 Linear Magnitude')
	legend('z211-1st Mode','z213-3rd Mode','z21-Total')
	text(.11,1.0*z211mag(1),z211text)
	text(.11,1.1*z213mag(1),z213text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	axis([.1 10 0 3]);
	grid

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

%	plot phase

	semilogx(w,z211phs,'k+-',w,z213phs,'kx-',w,z21phs,'k-')
	title('Transfer Function - z211, z213 and z21 Phase')
	legend('z211-1st Mode','z213-3rd Mode','z21-Total')
	xlabel('Frequency, rad/sec')
	ylabel('Phase, Deg')
	grid

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

	legend off;

%	z31, output 3 due to force 1 transfer functions	

	z311num = 1/3;

	z311den = [1 0 0];

	z312num = -1/2;

	z312den = [1 0 1];

	z313num = 1/6;

	z313den = [1 0 3];

	[z311mag,z311phs] = bode(z311num,z311den,w);	

	[z312mag,z312phs] = bode(z312num,z312den,w);	

	[z313mag,z313phs] = bode(z313num,z313den,w);	

	if  abs(z311phs(1)) >= 10
	
		z311text = '(-)';
		
	else
	
		z311text = '(+)';
		
	end				 

	if  abs(z312phs(1)) >= 10
	
		z312text = '(-)';
		
	else
	
		z312text = '(+)';
		
	end				 

	if  abs(z313phs(1)) >= 10
	
		z313text = '(-)';
		
	else
	
		z313text = '(+)';
		
	end				 

	z311magdb = 20*log10(z311mag);

	z312magdb = 20*log10(z312mag);

	z313magdb = 20*log10(z313mag);

%	calculate the complete transfer function

	z31 = ((1/3)./((j*w).^2) - ((1/2)./((j*w).^2 + 1)) + ((1/6)./((j*w).^2 + 3)));

	z31mag = abs(z31);

	z31magdb = 20*log10(z31mag);

	z31phs = z311phs + z312phs + z313phs;

%	plot the three modal contribution transfer functions and the total using
%   log magnitude versus frequency

	loglog(w,z311mag,'k+-',w,z312mag,'kx-',w,z313mag,'k.-',w,z31mag,'k-')
	title('Transfer Functions - z311, z312, z313 and z31 magnitude')
	legend('z311-1st Mode','z312-2nd Mode','z313-3rd Mode','z31-Total')
	text(.11,1.2*z311mag(1),z311text)
	text(.11,1.2*z312mag(1),z312text)
	text(.11,1.2*z313mag(1),z313text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	grid

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

%	plot the four transfer functions using db

	semilogx(w,z311magdb,'k+-',w,z312magdb,'kx-',w,z313magdb,'k.-',w,z31magdb,'k-')
	title('Transfer Function - z311, z312, z313 and z31 Magnitude')
	legend('z311-1st Mode','z312-2nd Mode','z313-3rd Mode','z31-Total')
	text(.11,2+z311magdb(1),z311text)
	text(.11,2+z312magdb(1),z312text)
	text(.11,2+z313magdb(1),z313text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude, db')
	grid

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

%	plot the four transfer functions using a linear magnitude scale so that
%	the amplitudes can be added directly

	semilogx(w,z311mag,'k+-',w,z312mag,'kx-',w,z313mag,'k.-',w,z31mag,'k-')
	title('Transfer Function - z311, z312, z313 and z31 Linear Magnitude')
	legend('z311-1st Mode','z312-2nd Mode','z313-3rd Mode','z31-Total')
	text(.11,1.0*z311mag(1),z311text)
	text(.11,1.1*z312mag(1),z312text)
	text(.11,1.1*z313mag(1),z313text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	grid

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

	semilogx(w,z311mag,'k+-',w,z312mag,'kx-',w,z313mag,'k.-',w,z31mag,'k-')
	title('Transfer Function - z311, z312, z313 and z31 Linear Magnitude')
	legend('z311-1st Mode','z312-2nd Mode','z313-3rd Mode','z31-Total')
	text(.11,1.0*z311mag(1),z311text)
	text(.11,1.1*z312mag(1),z312text)
	text(.11,1.1*z313mag(1),z313text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	axis([.1 10 0 3]);
	grid

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

%	plot phase

	semilogx(w,z311phs,'k+-',w,z312phs,'kx-',w,z313phs,'k.-',w,z31phs,'k-')
	title('Transfer Function - z311, z312, z313 and z31 Phase')
	legend('z311-1st Mode','z312-2nd Mode','z313-3rd Mode','z31-Total')
	xlabel('Frequency, rad/sec')
	ylabel('Phase, Deg')
	grid

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

	legend off;

%	z22, output 2 due to force 2 transfer functions	

	z221num = 1/3;

	z221den = [1 0 0];

%	z222 contribution = 0

	z223num = 2/3;

	z223den = [1 0 3];

	[z221mag,z221phs] = bode(z221num,z221den,w);	

	[z223mag,z223phs] = bode(z223num,z223den,w);	

	if  abs(z221phs(1)) >= 10
	
		z221text = '(-)';
		
	else
	
		z221text = '(+)';
		
	end				 

	if  abs(z223phs(1)) >= 10
	
		z223text = '(-)';
		
	else
	
		z223text = '(+)';
		
	end				 

	z221magdb = 20*log10(z221mag);

	z223magdb = 20*log10(z223mag);

%	calculate the complete transfer function

	z22 = ((1/3)./((j*w).^2) + ((2/3)./((j*w).^2 + 3)));

	z22mag = abs(z22);

	z22magdb = 20*log10(z22mag);

	z22phs = z221phs + z223phs;

%	plot the three modal contribution transfer functions and the total using
%   log magnitude versus frequency

	loglog(w,z221mag,'kx-',w,z223mag,'k.-',w,z22mag,'k-')
	title('Transfer Functions - z221, z223 and z22 magnitude')
	legend('z221-1st Mode','z223-3rd Mode','z22-Total')
	text(.11,1.2*z221mag(1),z221text)
	text(.11,1.2*z223mag(1),z223text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	grid

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

%	plot the four transfer functions using db

	semilogx(w,z221magdb,'kx-',w,z223magdb,'k.-',w,z22magdb,'k-')
	title('Transfer Function - z221, z223 and z22 Magnitude')
	legend('z221-1st Mode','z223-3rd Mode','z22-Total')
	text(.11,2+z221magdb(1),z221text)
	text(.11,2+z223magdb(1),z223text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude, db')
	grid

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

%	plot the four transfer functions using a linear magnitude scale so that
%	the amplitudes can be added directly

	semilogx(w,z221mag,'kx-',w,z223mag,'k.-',w,z22mag,'k-')
	title('Transfer Function - z221, z223 and z22 Linear Magnitude')
	legend('z221-1st Mode','z223-3rd Mode','z22-Total')
	text(.11,1.0*z221mag(1),z221text)
	text(.11,1.1*z223mag(1),z223text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	grid

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

	semilogx(w,z221mag,'kx-',w,z223mag,'k.-',w,z22mag,'k-')
	title('Transfer Function - z221, z223 and z22 Linear Magnitude')
	legend('z221-1st Mode','z223-3rd Mode','z22-Total')
	text(.11,1.0*z221mag(1),z221text)
	text(.11,1.1*z223mag(1),z223text)
	xlabel('Frequency, rad/sec')
	ylabel('Magnitude')
	axis([.1 10 0 3]);
	grid

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

%	plot phase

	semilogx(w,z221phs,'kx-',w,z223phs,'k.-',w,z22phs,'k-')
	title('Transfer Function - z221, z223 and z22 Phase')
	legend('z221-1st Mode','z223-3rd Mode','z22-Total')
	xlabel('Frequency, rad/sec')
	ylabel('Phase, Deg')
	grid

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

	legend off;