v2_057.m

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

% v2_057.m   
%
% This is a script file to solve a 2 DOF system
% given the mass, damping and stiffness matrices
% in dimensionless units and plot any desired 
% frequency response function.  
%            
% Plot Formats for FRF and IRF Measurements
% Figures for UC-SDRL-CN-20-263-662, Appendix A
 
%**********************************************************************
% 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;
axis('normal');
axis([1 2 3 4]),axis;
pi=3.14159265;
alpha=5;
mass=[8,0;0,10];
stiff=[600000,-200000;-200000,200000];    
damp=alpha.*mass;
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(b,-a);
% Sort Modal Frequencies  
orig_lambda=diag(d);
[Y,I]=sort(imag(orig_lambda));
lambda=orig_lambda(I);
xx=x(:,I);
xx,lambda
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);
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);
xxx=input('Hit any key to continue');
lambda,residu
xxx=input('Hit any key to continue');
%      Calculate desired plots
t=linspace(0,1,501);
xt1=residu(1)./2*exp(lambda(1).*t);
xt2=residu(2)./2*exp(lambda(2).*t);
xt3=residu(3)./2*exp(lambda(3).*t);
xt4=residu(4)./2*exp(lambda(4).*t);
xt=xt1+xt2+xt3+xt4;
plot(t,xt)
%       Reset plot axis for same + and - limits
V=axis;
ymax=max([abs(V(3)),abs(V(4))]);V(3)=-ymax;V(4)=+ymax;axis(V);
%   
plot(t,xt)
xlabel('Time (Sec.)'),ylabel ('Amplitude'),grid
if plt==1,meta v2_057a,end;
pause
axis([1 2 3 4]),axis;
clear t; clear xt;
f=linspace(0,100,501);
om=f.*2.*pi; 
H1=residu(1)./(j.*om-lambda(1));
H2=residu(2)./(j.*om-lambda(2));
H3=residu(3)./(j.*om-lambda(3));
H4=residu(4)./(j.*om-lambda(4));
H=H1+H2+H3+H4;
plot(f,real(H))
V=axis;
ymax=max([abs(V(3)),abs(V(4))]);V(3)=-ymax;V(4)=+ymax;axis(V);
plot(f,real(H))
xlabel('Frequency (Hertz)'),ylabel('Amplitude (Real)'),grid
if plt==1,meta v2_057b,end;
pause
axis([1 2 3 4]),axis;
plot(f,imag(H))
V=axis;
ymax=max([abs(V(3)),abs(V(4))]);V(3)=-ymax;V(4)=+ymax;axis(V);
plot(f,imag(H))
xlabel('Frequency (Hertz)'),ylabel('Amplitude (Imag)'),grid
if plt==1,meta v2_057c,end;
pause
axis([1 2 3 4]),axis;
plot(f,abs(H))
xlabel('Frequency (Hertz)'),ylabel('Magnitude'),grid
if plt==1,meta v2_057d,end;
pause
axis([1 2 3 4]),axis;
plot(f,20.*log10(abs(H)))
xlabel('Frequency (Hertz)'),ylabel('Log Magnitude (dB)'),grid
if plt==1,meta v2_057e,end;
pause
plot(f,360./(2.*pi).*angle(H))
V=axis;
ymax=max([abs(V(3)),abs(V(4))]);V(3)=-ymax;V(4)=+ymax;axis(V);
plot(f,360./(2.*pi).*angle(H))
xlabel('Frequency (Hertz)'),ylabel('Phase (Degrees)'),grid
if plt==1,meta v2_057f,end;
pause
axis([1 2 3 4]),axis;
axis('square')
polar(angle(H),abs(H))
xlabel('Amplitude (Real)'),ylabel('Amplitude (Imag)'),grid
if plt==1,meta v2_057g,end;