%Programa análise Modal
%Dados e entrada:
clc
clear all
% Unidades em SI
%% Calculo do Keq
Db = 4.5e-3; %Diametro bambu
Ib = (pi*Db^4)/64 %Momento de inercia do bambu
l = 0.25; %comprimento da vareta de bambu (entre as placas)
Eb = 9.02E9; %Modulo de elasticidade do bambu
n = 8; %Numero de varetas dO bambu
Keq = n*((12*Eb*Ib)/(l^3)); %Rigidez equivalente (viga engastada)
M1 = 0.335; %Massa do primeiro andar
M2 = 0.330; %Massa do segundo andar
M3 = 0.334; %Massa do terceiro andar
M = [M1 0 0; 0 M2 0; 0 0 M3]; %Matriz de massa
K = [2*Keq -Keq 0; -Keq 2*Keq -Keq; 0 -Keq Keq;]; %matriz de rigidez
X0 = [0; 0; 0]; %Desl. Inicial
V0 = [0; 0; 0]; %Veloc. Inicial
%% 1 - M121inv
M12 = chol(M)';
M12inv = M12^(-1);
%% 2 - Matril Ktil
Ktil = M12inv*K*M12inv;
%% 3,4 - Problemas de autovalores
[Vt Vl] = eig(Ktil);
omega = sqrt(diag(Vl));
%% 5 - Matriz S e Sinv
S = M12inv*Vt;
Sinv = Vt'*M12;
%% 6 - Condições inicais nas coordenas nodais
r0 = Sinv*X0;
rv0 = Sinv*V0;
%% 7 - Solução do problema nas coordenas nodais
t0 = 0; %tempo incial
tf = 10; %tempo final
dt = tf/100; %Delta tempo
t = t0:dt:tf;
% Amplitude e angulo phi:
for i=1:length(X0)
A(i) = sqrt(r0(i)^2 + (rv0(i)/omega(i))^2);
phi(i) = atan2 (r0(i)*omega(i), rv0(i));
end
