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