v2_056.m

美国辛辛那提大学的振动教材

% v2_056.m   
%
% This is a script file to solve a 2 DOF system
% given the mass, damping and stiffness matrices
% in dimensionless units.  The output includes poles, 
% residues(modal coefficients) and time and frequency 
% domain plots of impulse and frequency response
% functions.
%
% 2 DOF Non-Proportionally Damped System
% Figures for UC-SDRL-CN-20-263-662, Chapter 5
  
%**********************************************************************
% Author: Randall J. Allemang
% Date:	18-Apr-94
% Structural Dynamics Research Lab
% University of Cincinnati
% Cincinnati, Ohio  45221-0072
% TEL:  [phone]
% FAX:  [phone]
% E-MAIL: [e-mail]
%*********************************************************************
%
clear,clg
plt=input('Store plots to file (Yes=1): (0)');if isempty(plt),plt=0;end;
pi=3.14159265;
mass=[5,0;0,10];
stiff=[4,-2;-2,6];    
damp=[6,-4;-4,5]; 
null=[0,0;0,0];
%    Form 2N x 2N state space equation.  
a=[null,mass;mass,damp];
b=[-mass,null;null,stiff];
[x,d]=eig(-inv(a)*b);
% Sort Modal Frequencies  
orig_lambda=diag(d);
[Y,I]=sort(imag(orig_lambda));
lambda=orig_lambda(I);
xx=x(:,I);
xxx=input('Hit any key to continue');
%     Compute 'modal a' and 'modal b' matrix
ma=xx.'*a*xx;
mb=xx.'*b*xx;
%     Extract modal vectors from state-space formulation
psi(1:2,1)=xx(3:4,1);
psi(1:2,2)=xx(3:4,2);
psi(1:2,3)=xx(3:4,3);
psi(1:2,4)=xx(3:4,4);
%      Calculate residue matrices
A1=psi(1:2,1)*psi(1:2,1).'./ma(1,1);
A2=psi(1:2,2)*psi(1:2,2).'./ma(2,2);
A3=psi(1:2,3)*psi(1:2,3).'./ma(3,3);
A4=psi(1:2,4)*psi(1:2,4).'./ma(4,4);
%      Formulate H(1,2) FRF as Default
resp=input('Enter response (row) DOF: (1)');if isempty(resp),resp=1;end;
inp=input('Enter input (column) DOF: (2)');if isempty(inp),inp=2;end;
residu(1) = A1(resp,inp);
residu(2) = A2(resp,inp);
residu(3) = A3(resp,inp);
residu(4) = A4(resp,inp);
A1,A2,A3,A4
xxx=input('Hit any key to continue');
lambda,residu
xxx=input('Hit any key to continue');
%      Calculate desired plots
t=linspace(0,30,500);
xt1=residu(1).*exp(lambda(1).*t) + residu(2).*exp(lambda(2).*t);
xt2=residu(3).*exp(lambda(3).*t) + residu(4).*exp(lambda(4).*t);
xt=xt1+xt2;
axis([0,30,-0.10,0.10])
plot(t,xt)
xlabel('Time (Sec.)'),ylabel('Amplitude'),grid
if plt==1,print -f1 -deps v2_056c,end;
pause
clg
plot(t,xt1)
xlabel('Time (Sec.)'),ylabel('Amplitude'),grid
if plt==1,print -f1 -deps v2_056i,end;
pause
clg
plot(t,xt2)
xlabel('Time (Sec.)'),ylabel('Amplitude'),grid
if plt==1,print -f1 -deps v2_056f,end;
pause
clg
clear t; clear xt; clear xt1; clear xt2;
f=linspace(0,2.5,500);
xf1=residu(1)./(j.*f-lambda(1)) + residu(2)./(j.*f-lambda(2));
xf2=residu(3)./(j.*f-lambda(3)) + residu(4)./(j.*f-lambda(4));
xf=xf1+xf2;
axis([0,2.5,-0.6,0.6])
plot(f,real(xf))
xlabel('Frequency (Rad/Sec)'),ylabel('Real'),grid
if plt==1,print -f1 -deps v2_056a,end;
pause
clg
plot(f,imag(xf))
xlabel('Frequency (Rad/Sec)'),ylabel('Imaginary'),grid
if plt==1,print -f1 -deps v2_056b,end;
pause
clg
plot(f,real(xf1))
xlabel('Frequency (Rad/Sec)'),ylabel('Real'),grid
if plt==1,print -f1 -deps v2_056g,end;
pause
clg
plot(f,imag(xf1))
xlabel('Frequency (Rad/Sec)'),ylabel('Imaginary'),grid
if plt==1,print -f1 -deps v2_056h,end;
pause
clg
axis([0,2.5,-0.6,0.6])
plot(f,real(xf2))
xlabel('Frequency (Rad/Sec)'),ylabel('Real'),grid
if plt==1,print -f1 -deps v2_056d,end;
pause
clg
plot(f,imag(xf2))
xlabel('Frequency (Rad/Sec)'),ylabel('Imaginary'),grid
if plt==1,print -f1 -deps v2_056e,end;
axis([1 2 3 4]);axis;