v2_099.m

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

% v2_099.m   
%
% This is a script file to solve a 3 DOF system
% given the mass, damping and stiffness matrices
% in dimensionless units and plot any desired 
% frequency response function.  
%            
% Solution for FRF by Inverse of [B] 
             
%**********************************************************************
% 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,clf,close all;
pi=3.14159265;
plt=input('Store plots to file (Yes=1): (0)');if isempty(plt),plt=0;end;
%
%  solve 3 dof system 
%
mass=[10,0,0;0,14,0;0,0,12];
stiff=[5000,-3000,0;-3000,5500,-2500;0,-2500,2500];    
damp=[50,-30,0;-30,55,-25;0,-25,25];    
null=[0,0,0;0,0,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);
%     Normalize x matrix to real vectors if possible
for ii=1:6
xx(1:6,ii)=xx(1:6,ii)./xx(4,ii);
end
%     Compute 'modal a' and 'modal b' matrix
ma=xx.'*a*xx;
mb=xx.'*b*xx;
%     Extract modal vectors from state-space formulation
psi(1:3,1)=xx(4:6,1);
psi(1:3,2)=xx(4:6,2);
psi(1:3,3)=xx(4:6,3);
psi(1:3,4)=xx(4:6,4);
psi(1:3,5)=xx(4:6,5);
psi(1:3,6)=xx(4:6,6);
lambda
xxx=input('Hit any key to continue');
psi
xxx=input('Hit any key to continue');
ma
xxx=input('Hit any key to continue');
clear a,clear b,clear x,clear d;
 
omega=linspace(0,100,501);
for ii=1:501;
FRF(:,:,ii)=inv(-mass.*omega(ii).*omega(ii) + damp.*j.*omega(ii) + stiff);
end;

for Ni=1:3;
for No=1:3;
H=squeeze(FRF(No,Ni,:));
figure;
subplot(211),semilogy(omega,abs(H))
xlabel('Frequency (Rad/Sec)'),ylabel('Magnitude'),grid
scale=360.0/(2.0*pi);
subplot(212),plot(omega,scale.*angle(H))
xlabel('Frequency (Rad/Sec)'),ylabel('Phase'),grid
title('Frequency Response Function: H(x,x)');
pause
end
end
if plt==1,print -f1 -deps v2_099a,end;
pause