Document Actions
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: 513-556-2725
% FAX: 513-556-3390
% E-MAIL: randy.allemang@uc.edu
%*********************************************************************
%
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;