v2_100.m

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

% v2_100.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,clg
pi=3.14159265;
plt=input('Store plots to file (Yes=1): (0)');if isempty(plt),plt=0;end;
%
%  solve 3 dof system 
%
mass=[8,0,0;0,10,0;0,0,12];
stiff=10000.*[26,-12,0;-12,22,-10;0,-10,10];    
damp=[13,-6,0;-6,11,-5;0,-5,5];    
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);
modal_mass1=psi(:,1:3).'*mass*psi(:,1:3);
modal_mass2=psi(:,4:6).'*mass*psi(:,4: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');
modal_mass1
xxx=input('Hit any key to continue');
modal_mass2
xxx=input('Hit any key to continue');
clear a,clear b,clear x,clear d;
%      Set up inversion loop
delta=1.0;
for counter=1:300;
om=(counter-1)*delta;
omega(counter)=om;
B=-mass.*om.*om + damp.*j.*om + stiff;
H=inv(B);
FRF(1,counter)=H(1,1);
FRF(2,counter)=H(1,2);
FRF(3,counter)=H(1,3);
FRF(4,counter)=H(2,2);
FRF(5,counter)=H(2,3);
FRF(6,counter)=H(3,3);
end
for counter=1:6;
clg
subplot(211),semilogy(omega,abs(FRF(counter,:)))
xlabel('Frequency (Rad/Sec)'),ylabel('Magnitude'),grid
scale=360.0/(2.0*pi);
subplot(212),plot(omega,scale.*angle(FRF(counter,:)))
xlabel('Frequency (Rad/Sec)'),ylabel('Phase'),grid
if counter==1,title('Frequency Response Function: H(1,1)'),end;
if counter==2,title('Frequency Response Function: H(1,2)'),end;
if counter==3,title('Frequency Response Function: H(1,3)'),end;
if counter==4,title('Frequency Response Function: H(2,2)'),end;
if counter==5,title('Frequency Response Function: H(2,3)'),end;
if counter==6,title('Frequency Response Function: H(3,3)'),end;
pause
end
if plt==1,print -f1 -deps v2_100a,end;