v2_095.m

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

% v2_095.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.  A third DOF is
% then added to act as a spring-mass-damper to 
% reduce the magnitude of the first mode of the 
% original system.  The 3 DOF system is then
% solved and any frequency response functions
% can be plotted.
             
%**********************************************************************
% 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;
alpha=0.001
mass=[8,0;0,10];
stiff=[6000,-2000;-2000,2000];    
damp=alpha.*stiff;
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);
%     Normalize x matrix to real vectors if possible
for ii=1:4
xx(1:4,ii)=xx(1:4,ii)./xx(3,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: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);
lambda
xxx=input('Hit any key to continue');
psi
xxx=input('Hit any key to continue');
ma
xxx=input('Hit any key to continue');
%      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: (1)');if isempty(inp),inp=1;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
om=linspace(0,50,500);
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));
HH=H1+H2+H3+H4;
fig1=figure(1);
subplot(211),semilogy(om,abs(HH))
xlabel('Frequency (Rad/Sec)'),ylabel('Magnitude'),grid
title('Frequency Response Function')
subplot(212),plot(om,angle(HH))
xlabel('Frequency (Rad/Sec)'),ylabel('Phase'),grid
pause
if plt==1,print -f1 -deps v2_095a,end;
%
% now solve 3 dof system (with spring-mass-damper added)
%
pi=3.14159265;
alpha=0.001
mass=[8,0,0;0,10,0;0,0,2.50];
stiff=[6000,-2000,0;-2000,2300,-300;0,-300,300];    
damp=alpha.*stiff;
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);
x,d
xxx=input('Hit any key to continue');
%     Pick frequency roots
lambda=diag(d);
%     Normalize x matrix to real vectors if possible
for ii=1:6
x(1:6,ii)=x(1:6,ii)./x(4,ii);
end
%     Compute 'modal a' and 'modal b' matrix
ma=x.'*a*x;
mb=x.'*b*x;
%     Extract modal vectors from state-space formulation
psi(1:3,1)=x(4:6,1);
psi(1:3,2)=x(4:6,2);
psi(1:3,3)=x(4:6,3);
psi(1:3,4)=x(4:6,4);
psi(1:3,5)=x(4:6,5);
psi(1:3,6)=x(4:6,6);
%      Calculate residue matrices
A1=psi(1:3,1)*psi(1:3,1).'./ma(1,1);
A2=psi(1:3,2)*psi(1:3,2).'./ma(2,2);
A3=psi(1:3,3)*psi(1:3,3).'./ma(3,3);
A4=psi(1:3,4)*psi(1:3,4).'./ma(4,4);
A5=psi(1:3,5)*psi(1:3,5).'./ma(5,5);
A6=psi(1:3,6)*psi(1:3,6).'./ma(6,6);
residu(1) = A1(resp,inp);
residu(2) = A2(resp,inp);
residu(3) = A3(resp,inp);
residu(4) = A4(resp,inp);
residu(5) = A5(resp,inp);
residu(6) = A6(resp,inp);
A1,A2,A3,A4,A5,A6
xxx=input('Hit any key to continue');
lambda,residu
xxx=input('Hit any key to continue');
%      Calculate desired plots
om=linspace(0,50,500);
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));
H5=residu(5)./(j.*om-lambda(5));
H6=residu(6)./(j.*om-lambda(6));
H=H1+H2+H3+H4+H5+H6;
fig2=figure(2);
subplot(211),semilogy(om,abs(H),om,abs(HH))
xlabel('Frequency (Rad/Sec)'),ylabel('Magnitude'),grid
title('Frequency Response Function')
subplot(212),plot(om,angle(H),om,angle(HH))
xlabel('Frequency (Rad/Sec)'),ylabel('Phase'),grid
pause
if plt==1,print -f2 -deps v2_095b,end;
% Now compute final H22 using impedance method
M3=2.5;
C3=.3;
K3=300;
H33=-M3.*om.*om+j.*C3.*om+K3;
H33=H33./((-M3.*om.*om).*(K3+j.*om.*C3));
HHH=HH.*H33./(HH+H33);
fig3=figure(3);
subplot(211),semilogy(om,abs(H),om,abs(HHH))
xlabel('Frequency (Rad/Sec)'),ylabel('Magnitude'),grid
title('Frequency Response Function')
subplot(212),plot(om,angle(H),om,angle(HHH))
xlabel('Frequency (Rad/Sec)'),ylabel('Phase'),grid
pause
if plt==1,print -f3 -deps v2_095c,end;