tdofxfer.m

振动仿真 vibration simulation using matlab and ansys!

	echo off
%	tdofxfer.m	plotting frequency responses of three dof model
%	Plots tdof model poles and zeros in complex plane, user choice of
%	damping values.  Uses several different model descriptions and 
%	frequency response calculating techniques.  The model is described
%	in polynomial, transfer function and zpk forms.  Magnitude and phase
%	versus frequency are calculated using a scalar frequency "for-loop",
%	vector frequency, automatic bode plotting and bode with magnitude
%	and frequency outputs.	

	clf;

	legend off;

	subplot(1,1,1);

	clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%	"Polynomial Form, for-loop" frequency response plotting

%	assign values for masses, damping, and stiffnesses

	m1 = 1;
	m2 = 1;
	m3 = 1;
	c1 = 0;
	c2 = 0;
	k1 = 1;
	k2 = 1;

%	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,200);

%	pre-calculate the radians to degree conversion

	rad2deg = 180/pi;

%	Use a for-loop to cycle through all the frequencies, using Matlab's 
%	complex algebra capabilities to evaluate.

	for  cnt = 1:length(w)

%	define s as the imaginary operator times each frequency	

		s = j*w(cnt);

%	define the frequency responses to be evaluated

		den(cnt) = s^2*(s^4*(m1*m2*m3) + s^3*(m2*m3*c1 + m1*m3*c1 + m1*m2*c2 ...
						+ m1*m3*c2) + s^2*(m1*m3*k1 + m1*m3*k2 + m1*m2*k2 ...
						+ m2*c1*c2 + m3*c1*c2 + m1*c1*c2 + k1*m2*m3) ...
						+ s*(m3*c1*k2 + m2*c2*k1 + m1*c2*k1 + m1*c1*k2 ...
						+ m3*c2*k1 + m2*c1*k2) + (m1*k1*k2 + m2*k1*k2 + m3*k1*k2));

		z11bf(cnt) = ((m2*m3)*s^4 + (m3*c1 + m3*c2 + m2*c2)*s^3 ...
						+ (c1*c2 + m2*k2 + m3*k1 + m3*k2)*s^2 ...
						+ (c1*k2 + c2*k1)*s + (k1*k2))/den(cnt);

		z21bf(cnt) = ((m3*c1)*s^3 + (c1*c2 + m3*k1)*s^2 + (c1*k2 + c2*k1)*s ...
						+ (k1*k2))/den(cnt);

		z31bf(cnt) = ((c1*c2)*s^2 + (c1*k2 + c2*k1)*s + (k1*k2))/den(cnt);

		z22bf(cnt) = ((m1*m3)*s^4 + (m1*c2 + m3*c1)*s^3 + (m1*k2 + c1*c2 + m3*k1)*s^2 ...
						+ (c1*k2 + c2*k1)*s + (k1*k2))/den(cnt);
	
%	calculate the magnitude and phase of each frequency response

		z11bfmag(cnt) = abs(z11bf(cnt));

		z21bfmag(cnt) = abs(z21bf(cnt));

		z31bfmag(cnt) = abs(z31bf(cnt));

		z22bfmag(cnt) = abs(z22bf(cnt));

		z11bfphs(cnt) = angle(z11bf(cnt))*rad2deg;

		z21bfphs(cnt) = angle(z21bf(cnt))*rad2deg;

		z31bfphs(cnt) = angle(z31bf(cnt))*rad2deg;

		z22bfphs(cnt) = angle(z22bf(cnt))*rad2deg;

	end		%	end of for-loop

%	plot the z11, z33 frequency responses

	loglog(w,z11bfmag,'k-')
	title('Polynomial Form, for-loop z11, z33 magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude')
	grid

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

	semilogx(w,z11bfphs,'k*-')
	title('Polynomial Form, for-loop z11 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	grid

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

%	plot the z21, z12, z23, z32 frequency responses

	loglog(w,z21bfmag,'k-')
	title('Polynomial Form, for-loop z21, z12, z23, z32 magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude')
	grid

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

	semilogx(w,z21bfphs,'k*-')
	title('Polynomial Form, for-loop z21, z12, z23, z32 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	grid

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

%	plot the z13, z31 frequency responses

	loglog(w,z31bfmag,'k-')
	title('Polynomial Form, for-loop z13, z31 magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude')
	grid

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

	semilogx(w,z31bfphs,'k*-')
	title('Polynomial Form, for-loop z13, z31 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	grid

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

%	plot the z22 frequency responses

	loglog(w,z22bfmag,'k-')
	title('Polynomial Form, for-loop z22 magnitude')
	xlabel('frequency, rad/sec')
	ylabel('magnitude')
	grid

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

	semilogx(w,z22bfphs,'k*-')
	title('Polynomial Form, for-loop z22 phase')
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	grid

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

	legend off;

%	plot the four frequency responses

	loglog(w,z11bfmag,w,z21bfmag,w,z31bfmag,w,z22bfmag)
	title('Polynomial Form, for-loop z11, z21, z31 and z22 magnitude')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('magnitude')
	grid

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

	semilogx(w,z11bfphs,w,z21bfphs,w,z31bfphs,w,z22bfphs)
	title('Polynomial Form, for-loop z11, z21, z31 and z22 phase')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	grid

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

	legend off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	clear all;

	clf;

	legend off;

	subplot(1,1,1);

%	"Polynomial Form, Vector" method - using Matlab's vector capabilities instead
%	of the "for" loop.

%	assign values for masses, damping, and stiffnesses

	m1 = 1;
	m2 = 1;
	m3 = 1;
	c1 = 0;
	c2 = 0;
	k1 = 1;
	k2 = 1;

%	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,200);

%	pre-calculate the radians to degree conversion

	rad2deg = 180/pi;

%	define s as the imaginary operator times the radian frequency vector

	s = j*w;

%	define the frequency responses to be evaluated, using the "." prefix 
%	in front of each operator to indicate that each 

%	define the frequency responses to be evaluated

	den = s.^2.*(s.^4*(m1*m2*m3) + s.^3*(m2*m3*c1 + m1*m3*c1 + m1*m2*c2 ...
					+ m1*m3*c2) + s.^2*(m1*m3*k1 + m1*m3*k2 + m1*m2*k2 ...
					+ m2*c1*c2 + m3*c1*c2 + m1*c1*c2 + k1*m2*m3) ...
					+ s*(m3*c1*k2 + m2*c2*k1 + m1*c2*k1 + m1*c1*k2 ...
					+ m3*c2*k1 + m2*c1*k2) + (m1*k1*k2 + m2*k1*k2 + m3*k1*k2));

	z11bfv = ((m2*m3)*s.^4 + (m3*c1 + m3*c2 + m2*c2)*s.^3 ...
					+ (c1*c2 + m2*k2 + m3*k1 + m3*k2)*s.^2 ...
					+ (c1*k2 + c2*k1)*s + (k1*k2))./den;

	z21bfv = ((m3*c1)*s.^3 + (c1*c2 + m3*k1)*s.^2 + (c1*k2 + c2*k1)*s ...
					+ (k1*k2))./den;

	z31bfv = ((c1*c2)*s.^2 + (c1*k2 + c2*k1)*s + (k1*k2))./den;

	z22bfv = ((m1*m3)*s.^4 + (m1*c2 + m3*c1)*s.^3 + (m1*k2 + c1*c2 + m3*k1)*s.^2 ...
					+ (c1*k2 + c2*k1)*s + (k1*k2))./den;

%	calculate the magnitude and phase of each frequency response

		z11bfvmag = abs(z11bfv);

		z21bfvmag = abs(z21bfv);

		z31bfvmag = abs(z31bfv);

		z22bfvmag = abs(z22bfv);

		z11bfvphs = angle(z11bfv)*rad2deg;

		z21bfvphs = angle(z21bfv)*rad2deg;

		z31bfvphs = angle(z31bfv)*rad2deg;

		z22bfvphs = angle(z22bfv)*rad2deg;

%	plot the four frequency responses

	loglog(w,z11bfvmag,w,z21bfvmag,w,z31bfvmag,w,z22bfvmag)
	title('Polynomial Form, Vector z11, z21, z31 and z22 magnitude')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('magnitude')
	grid

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

	semilogx(w,z11bfvphs,w,z21bfvphs,w,z31bfvphs,w,z22bfvphs)
	title('Polynomial Form, Vector z11, z21, z31 and z22 phase')
	legend('z11,z33','z21,z12,z23,z32','z31,z13','z22',3)
	xlabel('frequency, rad/sec')
	ylabel('phase, deg')
	grid

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

	legend off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	clear all;

	clf;

	legend off;

	subplot(1,1,1);

%	using Matlab's automatic "bode" plotting capability, defining the transfer
%	functions in "transfer function" form by row vectors of coefficients of "s"

%	assign values for masses, damping, and stiffnesses

	m1 = 1;
	m2 = 1;
	m3 = 1;
	c1 = 0;
	c2 = 0;
	k1 = 1;
	k2 = 1;

%	define row vectors of numerator and denominator coefficients

	den = [(m1*m2*m3) (m2*m3*c1 + m1*m3*c1 + m1*m2*c2 + m1*m3*c2) ...
	       (m1*m3*k1 + m1*m3*k2 + m1*m2*k2 + m2*c1*c2 + m3*c1*c2 + ...
	        m1*c1*c2 + k1*m2*m3) ...
		   (m3*c1*k2 + m2*c2*k1 + m1*c2*k1 + m1*c1*k2 + m3*c2*k1 + m2*c1*k2) ...
		   (m1*k1*k2 + m2*k1*k2 + m3*k1*k2) 0 0];

	z11num = [(m2*m3) (m3*c1 + m3*c2 + m2*c2) (c1*c2 + m2*k2 + m3*k1 + m3*k2) ...
			  (c1*k2 + c2*k1) (k1*k2)];

	z21num = [(m3*c1) (c1*c2 + m3*k1) (c1*k2 + c2*k1) (k1*k2)];

	z31num = [(c1*c2) (c1*k2 + c2*k1) (k1*k2)];

	z22num = [(m1*m3) (m1*c2 + m3*c1) (m1*k2 + c1*c2 + m3*k1) ...
			  (c1*k2 + c2*k1) (k1*k2)];

%	the bode command with no left-hand side arguments automatically chooses
%	frequency limits and plots results

	grid on
	bode(z11num,den);

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

	bode(z21num,den);

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

	bode(z31num,den);

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

	bode(z22num,den);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	clear all;

	clf;

	legend off;

	subplot(1,1,1);

%	using Matlab's "bode" plotting capability, defining the transfer
%	functions in "transfer function" form by row vectors of coefficients of
%	"s"	and defining output vectors for magnitude and phase as well as a
%	defined range of radian frequencies

%	assign values for masses, damping, and stiffnesses

	m1 = 1;
	m2 = 1;
	m3 = 1;
	c1 = 0;
	c2 = 0;
	k1 = 1;
	k2 = 1;

%	define row vectors of numerator and denominator coefficients

	den = [(m1*m2*m3) (m2*m3*c1 + m1*m3*c1 + m1*m2*c2 + m1*m3*c2) ...
	       (m1*m3*k1 + m1*m3*k2 + m1*m2*k2 + m2*c1*c2 + m3*c1*c2 + ...
	        m1*c1*c2 + k1*m2*m3) ...
		   (m3*c1*k2 + m2*c2*k1 + m1*c2*k1 + m1*c1*k2 + m3*c2*k1 + m2*c1*k2) ...
		   (m1*k1*k2 + m2*k1*k2 + m3*k1*k2) 0 0];

	z11num = [(m2*m3) (m3*c1 + m3*c2 + m2*c2) (c1*c2 + m2*k2 + m3*k1 + m3*k2) ...
			  (c1*k2 + c2*k1) (k1*k2)];

	z21num = [(m3*c1) (c1*c2 + m3*k1) (c1*k2 + c2*k1) (k1*k2)];

	z31num = [(c1*c2) (c1*k2 + c2*k1) (k1*k2)];

	z22num = [(m1*m3) (m1*c2 + m3*c1) (m1*k2 + c1*c2 + m3*k1) ...
			  (c1*k2 + c2*k1) (k1*k2)];

%	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,200);

%	calculate the rigid-body motions for low and high frequency portions
%	of all the frequency responses, the denominator entries are vectors with
%	entries being coefficients of the "s" terms in the low and high frequency
%	asymptotes, starting with the highest power of "s" and ending with the 
%	"0"th power of "s" or the constant term

	z11num_lo = [1];

	z11den_lo = [3 0 0];			%  -1/(3*w^2)

	z11num_hi = [1];

	z11den_hi = [1 0 0];			%	-1/(w^2)

	z21num_lo = [1];

	z21den_lo = [3 0 0];			%  -1/(3*w^2)

	z21num_hi = [1];

	z21den_hi = [1 0 0 0 0];		%  -1/(3*w^4)
	  									  
	z31num_lo = [1];

	z31den_lo = [3 0 0];			%  -1/(3*w^2)

	z31num_hi = [1];

	z31den_hi = [1 0 0 0 0 0 0];	%  -1/(w^2)

	z22num_lo = [1];

	z22den_lo = [3 0 0];			%  -1/(3*w^2)

	z22num_hi = [1];

	z22den_hi = [1 0 0];			%  -1/(w^2)

%	define the "tf" models from "num, den" combinations

	z11tf = tf(z11num,den);

	z21tf = tf(z21num,den);

	z31tf = tf(z31num,den);