% Trabalho final - Vibrações Mecânicas - 2025/2
% Viga engastada-livre em translação lateral

clc
clearvars
close all

% dados dimensionais e materiais da viga
E=200e9; % Pa
ro=7850; % kg/m^3
ep=6.25e-3; % m
lg=49.8e-3;  % m
cp=580e-3; % m

% massa de motor e complementos
mm=0.245; % kg

% número de divisões
ngdl=10;

% matriz de inércia
L=cp/ngdl;
me=ro*L*lg*ep;
M=[(me/2)+mm 0 0 0 0 0 0 0 0 0; 0 me 0 0 0 0 0 0 0 0; ...
    0 0 me 0 0 0 0 0 0 0; 0 0 0 me 0 0 0 0 0 0; 0 0 0 0 me 0 0 0 0 0; ...
    0 0 0 0 0 me 0 0 0 0; 0 0 0 0 0 0 me 0 0 0; 0 0 0 0 0 0 0 me 0 0; ...
    0 0 0 0 0 0 0 0 me 0; 0 0 0 0 0 0 0 0 0 me];

% matriz de rigidez
I=lg*ep^3/12;
A=((L^3)/(6*E*I))*...
  [2000 1701 1408 1127 864 625 416 243 112 29; ...
   1701 1458 1216  980 756 550 368 216 100 26; ...
   1408 1216 1024  833 648 475 320 189  88 23; ...
   1127  980  833  686 540 400 272 162  76 20; ...
    864  756  648  540 432 325 224 135  64 17; ...
    625  550  475  400 325 250 176 108  52 14; ...
    416  368  320  272 224 176 128  81  40 11; ...
    243  216  189  162 135 108  81  54  28  8; ...
    112  100   88   76  64  52  40  28  16  5; ...
     29   26   23   20  17  14  11   8   5  2];
K=inv(A);

% matriz de amortecimento
alfa=0;
beta=0.00005;
C=alfa*M+beta*K; % amortecimento proporcional

% 1a. Parte - Cálculo de frequências naturais e modos de vibrar

% passo 1 - decomposição de Choleski
L=chol(M);
Lt=L';
Li=inv(L);
Lti=inv(Lt);

% passo 2 - matriz de rigidez normalizada pela massa
Ktil=Lti*K*Li;

% passo 3 - problema de autovalores
[V,D]=eig(Ktil);
  % frequências naturais em rad/s
omn=sqrt(diag(D));
  % ordenamento das frequências naturais
[omn,iv]=sort(omn);
  % apresentação das frequências naturais em rad/s
disp('As frequências naturais teóricas, em rad/s, são')
disp(' ')
for im=1:ngdl
  disp(['ômega_',num2str(im),' = ',num2str(omn(im)),' rad/s'])
end
disp(' ')
  % frequências naturais em Hz
fn=omn/(2*pi);

% apresentação das frequências naturais em Hz
disp('As frequências naturais teóricas, em Hz, são')
disp(' ')
for im=1:ngdl
  disp(['fn_',num2str(im),' = ',num2str(fn(im)),' Hz'])
end
disp(' ')

% passo 4 - autovetores normalizados e ordenados
P=zeros(size(V));
for il=1:length(iv)
  P(:,il)=V(:,iv(il));
end

% passo 5 - modos de vibrar em colunas
T=Li*P; 
Tt=T';
Pt=P';
Ti=Pt*L;
  % inserção de engaste
Te=[T;zeros(1,10)]; 
  
% visualização dos modos de interesse
nint=3; % número de modos de interesse
vcf=(cp:(-cp/10):0)'; % vetor de coordenadas físicas
for im=1:nint
  figure(im)
  plot(Te(:,im),vcf,'o--')
  grid on
  xlabel('amplitude modal')
  ylabel('coordenada física (m)')
  axis([-1.5 1.5 0 0.6])
  mode=['modo ',num2str(im)];
  legend(mode)
end

% % 2a. Parte - Cálculo de razões de amortecimento
%   % razões de amortecimento
% zeta=(alfa./(2*omn))+(beta*omn/2);
%   % apresentação de razões de amortecimento
% disp('As razões de amortecimento teóricas são')
% disp(' ')
% for im=1:ngdl
%   disp(['zeta_',num2str(im),' = ',num2str(zeta(im))])
% end
% disp(' ')

% % 3a. Parte - Cálculo de vibrações sob excitação harmônica
% 
% % passo 6 - matriz espectral
% Lb=Pt*Ktil*P;
% 
% % passo 7 - condições iniciais em coordenadas modais  
%   % condições iniciais em coordenadas físicas
% x0=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0];
% xp0=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0];
%   % condições iniciais em coordenadas modais
% r0=Ti*x0;
% rp0=Ti*xp0;
% 
% % passo 8 - resolução das equações pertinentes
%   % vetor de tempo
% tf=20; % tempo final (s)
% t=linspace(0,tf,2048); 
%   % forças atuantes
% F=[10; 0; 0; 0; 0; 0; 0; 0; 0; 0]; % amplitudes
% omf=omn(1); % frequência (rad/s)
% fif=0; % fase
%   % vetor de amplitudes de forças modais
% Fm=Tt*F;
%   % forças modais harmônicas
% fm=Fm*sin(omf*t+fif);
%   % visualização das forças modais
% figure(nint+1)
% subplot(2,1,1)
% plot(t,fm(1,:),'r-',t,fm(2,:),'m-',t,fm(3,:),'b-',t,fm(4,:),'c-',t,fm(5,:),'g')
% axis([0 tf -15 15])
% grid on
% xlabel('tempo (s)')
% ylabel('força modal (N)')
% title('Forças modais')
% legend('força 1','força 2','força 3','força 4','força 5')
% subplot(2,1,2)
% plot(t,fm(6,:),'r-',t,fm(7,:),'m-',t,fm(8,:),'b-',t,fm(9,:),'c-',t,fm(10,:),'g')
% axis([0 tf -15 15])
% grid on
% xlabel('tempo (s)')
% ylabel('força modal (N)')
% legend('força 6','força 7','força 8','força 9','força 10')
%   % respostas em coordenadas modais
%     % respostas permanentes
% k=omn.^2;
% e=omf./omn;
% H=1./sqrt((1-e.^2).^2+(2*zeta.*e).^2);
% R=(Fm./k).*H;
% teta=atan2((2*zeta.*e),(1-e.^2));
% tetar=teta;
% rper=zeros(length(omn),length(t));
% for ir=1:ngdl
%   rper(ir,:)=R(ir)*sin(omf*t-tetar(ir));
% end
%     % respostas transientes (amort. subcrítico, ci's nulas)
% ni=omn.*sqrt(1-zeta.^2);
% delta=zeta.*omn;
% psi=atan2((2*zeta.*sqrt(1-zeta.^2)),(2*zeta.^2-(1-e.^2)));
% rtran=zeros(length(omn),length(t));
% for ir=1:length(omn)
%   rtran(ir,:)=R(ir)*(e(ir)*exp(-delta(ir)*t)/sqrt(1-zeta(ir)^2)).*...
%               sin(ni(ir)*t+psi(ir));
% end
%     % respostas completas
% rcom=rper+rtran;
%   % visualização das respostas em coordenadas modais
% figure(nint+2)
% subplot(5,2,1)
% plot(t,rper(1,:),'r-',t,rtran(1,:),'g-',t,rcom(1,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 1 (m)');
% title('Vibrações modais sob excitação harmônica')
% legend('per.','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))]) 
% subplot(5,2,2)
% plot(t,rper(2,:),'r-',t,rtran(2,:),'g-',t,rcom(2,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 2 (m)');
% title('Vibrações modais sob excitação harmônica')
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))]) 
% subplot(5,2,3)
% plot(t,rper(3,:),'r-',t,rtran(3,:),'g-',t,rcom(3,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 3 (m)');
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))])
% subplot(5,2,4)
% plot(t,rper(4,:),'r-',t,rtran(4,:),'g-',t,rcom(4,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 4 (m)');
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))])
% subplot(5,2,5)
% plot(t,rper(5,:),'r-',t,rtran(5,:),'g-',t,rcom(5,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 5 (m)');
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))])
% subplot(5,2,6)
% plot(t,rper(6,:),'r-',t,rtran(6,:),'g-',t,rcom(6,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 6 (m)');
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))])
% subplot(5,2,7)
% plot(t,rper(7,:),'r-',t,rtran(7,:),'g-',t,rcom(7,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 7 (m)');
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))])
% subplot(5,2,8)
% plot(t,rper(8,:),'r-',t,rtran(8,:),'g-',t,rcom(8,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 8 (m)');
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))])
% subplot(5,2,9)
% plot(t,rper(9,:),'r-',t,rtran(9,:),'g-',t,rcom(9,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 9 (m)');
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))])
% subplot(5,2,10)
% plot(t,rper(10,:),'r-',t,rtran(10,:),'g-',t,rcom(10,:),'b-')
% grid on
% xlabel('tempo (s)'); ylabel('cd. mod. 10 (m)');
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))])
% 
% % passo 9 - respostas globais em coordenadas físicas
% x=T*rcom;
% %   % contribuições dos modos
% % xm1=T(:,1)*rcom(1,:); xm2=T(:,2)*rcom(2,:); xm3=T(:,3)*rcom(3,:);
% % xm4=T(:,4)*rcom(4,:); xm5=T(:,5)*rcom(5,:); xm6=T(:,6)*rcom(6,:);
% % xm7=T(:,7)*rcom(7,:); xm8=T(:,8)*rcom(8,:); xm9=T(:,9)*rcom(9,:);
% % xm10=T(:,10)*rcom(10,:);
%   % visualização das respostas em coordenadas físicas
% figure(nint+3)
% subplot(5,2,1)
% plot(t,x(1,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 1 (m)');
% title('Vibrações em coordenadas físicas')
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])
% subplot(5,2,2)
% plot(t,x(2,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 2 (m)');
% title('Vibrações em coordenadas físicas')
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])
% subplot(5,2,3)
% plot(t,x(3,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 3 (m)');
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])
% subplot(5,2,4)
% plot(t,x(4,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 4 (m)');
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])
% subplot(5,2,5)
% plot(t,x(5,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 5 (m)');
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])
% subplot(5,2,6)
% plot(t,x(6,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 6 (m)');
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])
% subplot(5,2,7)
% plot(t,x(7,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 7 (m)');
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])
% subplot(5,2,8)
% plot(t,x(8,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 8 (m)');
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])
% subplot(5,2,9)
% plot(t,x(9,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 9 (m)');
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])
% subplot(5,2,10)
% plot(t,x(10,:))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 10 (m)');
% axis([min(t) max(t) round(min(min(x'))) round(max(max(x')))])

% % 4a. Parte - Cálculo de receptância de interesse
% 
%   % vetor de frequências
% nf=2048; % número de frequências
% fx=linspace(0,400,nf); % Hz
% omx=2*pi*fx; % rad/s
%   % elemento da matriz de receptâncias
%     % seleção do elemento de interesse
% es=1; % vibração
% er=9; % excitação
%     % montagem da função receptância selecionada
% Hsrx=zeros(ngdl,nf);
% Hsr=zeros(1,nf);
% for ix=1:nf
%   for im=1:ngdl
%     Hsrx(im,ix)=(T(es,im)*T(er,im))/...
%                 ((omn(im)^2-omx(ix)^2)+1i*(2*zeta(im)*omn(im)*omx(ix)));
%   end
%   Hsr(ix)=sum(Hsrx(:,ix));
% end
%     % visualização da função receptância selecionada
% ifrf=[num2str(es),num2str(er)];
% mfrf=['Esp. Amplitude de H_{',ifrf,'} (dB re 1m/N)'];
% afrf=['Esp. Fase de H_{',ifrf,'} (rad)'];
% tfrf=['Receptância H_{',ifrf,'}'];
% figure(nint+4)
% subplot(2,1,1)
% plot(fx,20*log10(abs(Hsr)))
% % axis([...])
% grid on
% xlabel('frequência (Hz)')
% ylabel(mfrf)
% title(tfrf)
% % hold on
% % plot(fx,20*log10(abs(Hsrx(1,:))),fx,20*log10(abs(Hsrx(2,:))),...
% %      fx,20*log10(abs(Hsrx(3,:))),fx,20*log10(abs(Hsrx(3,:))),...
% %      fx,20*log10(abs(Hsrx(5,:))),fx,20*log10(abs(Hsrx(6,:))),...
% %      fx,20*log10(abs(Hsrx(7,:))),fx,20*log10(abs(Hsrx(8,:))),...
% %      fx,20*log10(abs(Hsrx(9,:))),fx,20*log10(abs(Hsrx(10,:))))
% % legend('total','modo 1','modo 2','modo 3','modo 4','modo 5',...
% %                'modo 6','modo 7','modo 8','modo 9','modo 10')
% % hold off
% subplot(2,1,2)
% plot(fx,angle(Hsr))
% % axis([...])
% grid on
% xlabel('frequência (Hz)')
% ylabel(afrf)
% % hold on
% % plot(fx,angle(Hsrx(1,:)),fx,angle(Hsrx(2,:)),fx,angle(Hsrx(3,:)),...
% %      fx,angle(Hsrx(4,:)),fx,angle(Hsrx(5,:)),fx,angle(Hsrx(6,:)),...
% %      fx,angle(Hsrx(7,:)),fx,angle(Hsrx(8,:)),fx,angle(Hsrx(9,:)),...
% %      fx,angle(Hsrx(10,:)))
% % legend('total','modo 1','modo 2','modo 3','modo 4','modo 5',...
% %                'modo 6','modo 7','modo 8','modo 9','modo 10')
% % hold off

% % 5a. Parte - Apresentação de inertância e receptância experimentais
% 
%   % inertância experimental
%     % nomeação de arquivo de dados
% arqdf='grupo1_fr.xlsx'; 
%     % leitura de arquivo de dados
% dadosf=xlsread(arqdf,'C1:ST10');
%     % atribuição de dados
% fe=dadosf(1,:); % frequência (Hz)
% % cor=dadosf(2,:);  % coerência 
% Ie_1r=dadosf(3,:);  % inertância 1, parte real
% Ie_1i=dadosf(4,:);  % inertância 1, parte imaginária
% % Ie_1a=dadosf(5,:);  % inertância 1, amplitude
% % Ie_1f=dadosf(6,:);  % inertância 1, fase
% Ie_2r=dadosf(7,:);  % inertância 2, parte real
% Ie_2i=dadosf(8,:);  % inertância 2, parte imaginária
% Ie_2a=dadosf(9,:);  % inertância 2, amplitude
% Ie_2f=dadosf(10,:); % inertância 2, fase
%     % montagem de inertância experimental
% Iexp=Ie_2r+i*Ie_2i;
%   % receptância experimental
% fe2=fe; % frequência em receptância (Hz)
% fe2(1)=1e-6; % substituição de frequência inicial (Hz)
% Hexp=Iexp./(-(2*pi*fe2).^2);
%   % visualização de FRFs experimentais
% figure(nint+5)
% subplot(2,1,1)
% mafrf=['Espectro de Amplitude de I_{',ifrf,'} (dB re 1m/Ns^{2})'];
% tafrf=['Inertância Experimental I_{',ifrf,'}'];
% plot(fe,20*log10(abs(Iexp)),'r-')
% axis([5 40 -40 20])
% grid on
% xlabel('frequência (Hz)')
% ylabel(mafrf)
% title(tafrf)
% subplot(2,1,2)
% mafrf=['Espectro de Amplitude de H_{',ifrf,'} (dB re 1m/Ns^{2})'];
% tafrf=['Receptância Experimental H_{',ifrf,'}'];
% plot(fe2,20*log10(abs(Hexp)),'b-')
% axis([5 40 -120 -50])
% grid on
% xlabel('frequência (Hz)')
% ylabel(mafrf)
% title(tafrf)

% % 6a. Parte - Comparação entre receptâncias teórica e experimental
% 
%   % visualização de receptâncias
% mafrf=['Espectro de Amplitude de H_{',ifrf,'} (dB re 1m/Ns^{2})'];
% tafrf=['Receptância H_{',ifrf,'}'];
% figure(nint+6)
% plot(fx,20*log10(abs(Hsr)),'c-')
% hold on
% plot(fe2,20*log10(abs(Hexp)),'b-')
% axis([5 40 -120 -50])
% grid on
% xlabel('frequência (Hz)')
% ylabel(mafrf)
% title(tafrf)
% legend('teórica','experimental')
% hold off

% % 7a. Parte - Cálculo de vibração via Runge-Kutta
% % (ainda em desenvolvimento!)
% 
% % definição de tempo e condições iniciais
%   % limites de tempo
% ti=0; % instante inicial
% tf=20; % instante final
%   % vetor de tempo
% clear t; % uso de vetor de tempo distinto
% dt=tf/1024; % passo de tempo
% t=ti:dt:tf; % novo vetor de tempo
%   % condições iniciais
% ci=zeros(2*ngdl,1);
%   
% % resolução de equação matricial via Runge-Kutta
% [t,xn]=ode45(@(t,u)veh10(t,u,M,C,K,ngdl),t,ci);
% 
% % visualização de solução
% figure(nint+7)
% subplot(5,2,1)
% plot(t,xn(:,1))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 1 (m)');
% title('Vibrações em coordenadas físicas')
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])
% subplot(5,2,2)
% plot(t,xn(:,2))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 2 (m)');
% title('Vibrações em coordenadas físicas')
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])
% subplot(5,2,3)
% plot(t,xn(:,3))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 3 (m)');
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])
% subplot(5,2,4)
% plot(t,xn(:,4))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 4 (m)');
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])
% subplot(5,2,5)
% plot(t,xn(:,5))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 5 (m)');
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])
% subplot(5,2,6)
% plot(t,xn(:,6))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 6 (m)');
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])
% subplot(5,2,7)
% plot(t,xn(:,7))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 7 (m)');
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])
% subplot(5,2,8)
% plot(t,xn(:,8))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 8 (m)');
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])
% subplot(5,2,9)
% plot(t,xn(:,9))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 9 (m)');
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])
% subplot(5,2,10)
% plot(t,xn(:,10))
% grid on
% xlabel('tempo (s)'); ylabel('vib. pos. 10 (m)');
% axis([min(t) max(t) round(min(min(xn(:,1:10)))) round(max(max(xn(:,1:10))))])