FRF Issue of 4 DOF Frame
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Can AKYOL
el 14 de Jun. de 2021
Comentada: Mathieu NOE
el 17 de Jun. de 2021
Hi everybody, I am triying to analyze a 4 story shear building under and earthquake. I used state space formulation and obtained the displacement of each story. When I try to find FRF to observe either I can find the modal parameters from response histories or not, I can get no meaningful results. I do fft(Output)/fft(Input) but ı obtained meaningless results. Where ı do a mistake in my code ?
State space portion :
global CC EQ M
K = [10000 -5000 0 0;-5000 10000 -5000 0;0 -5000 10000 -5000;0 0 -5000 5000]; %in N/m2
M = [4000 0 0 0;0 4000 0 0;0 0 4000 0;0 0 0 4000]; % in kg
[V,WW] = eig(K,M);
W = sqrt(WW)*[1;1;1;1];
C = inv(transpose(V))*0.1*[W(1) 0 0 0;0 W(2) 0 0;0 0 W(3) 0;0 0 0 W(4)]*inv(V);
CC = [zeros(4) eye(4);-inv(M)*K -inv(M)*C];
EQ = xlsread('el_centro');
step = 0.02;
time = 0:step:max(EQ(:,1));
y0 = [0 0 0 0 0 0 0 0]; % [d1 d2 d3 d4 v1 v2 v3 v4]
[tsol,ysol] = ode23('testode_2D',time,y0);
plot(time,ysol(:,1:4))
title('El centro Floor Response')
legend('Fisrt','Second','Third','Forth')
xlabel('Time')
ylabel('Displacement (m)')
.............................................................................
function dy = testode_2D(t,y)
global CC EQ M
index = interp1(EQ(:,1),1:length(EQ(:,1)),t,'nearest');
if (index ~= 1) & (index ~= length(EQ(:,1)))
if t > EQ(index,1)
oran = (t - EQ(index,1))/0.02;
ag = (EQ(index+1,2) - EQ(index,2))*oran + EQ(index,2);
else
oran = (t - EQ(index-1,1))/0.02;
ag = (EQ(index,2) - EQ(index-1,2))*oran + EQ(index-1,2);
end
else
ag = EQ(index,2);
end
f1 = -M(1,1)*ag;
f2 = -M(2,2)*ag;
f3 = -M(3,3)*ag;
f4 = -M(4,4)*ag;
A00 = zeros(4);
FF = [A00;inv(M)]*[f1;f2;f3;f4];
dy = CC*y + FF;
FRF part :
EQ = xlsread('el_centro');
EQ(:,2) = 9.81*EQ(:,2);
EQ_floor_dat = xlsread('el_centro_dat');
fs = 50; % sampling freq
T = .02; % sampling period
N = length(EQ_floor_dat(:,1)); % number of signal data
time = EQ_floor_dat(:,1); % time vector
Output = EQ_floor_dat(:,5) ;
Input = EQ(:,2);
Y = fft(Output); % output fft
X = fft(Input); %input fft
P2 = abs(Y/N);
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(N/2))/N;
2 comentarios
Respuesta aceptada
Mathieu NOE
el 16 de Jun. de 2021
hello again
so this is the code a bit reworked and expanded
the "interesting" portion in the FRF plot is in the sub Hertz range.. so the quality of the plot is related to simulation time length. It would require to simulate on longer records to further improve the FRF plot below 1 Hz
Code :
clc
clearvars
global CC EQ M
K = [10000 -5000 0 0;-5000 10000 -5000 0;0 -5000 10000 -5000;0 0 -5000 5000]; %in N/m2
M = [4000 0 0 0;0 4000 0 0;0 0 4000 0;0 0 0 4000]; % in kg
[V,WW] = eig(K,M);
W = sqrt(WW)*[1;1;1;1];
C = inv(transpose(V))*0.1*[W(1) 0 0 0;0 W(2) 0 0;0 0 W(3) 0;0 0 0 W(4)]*inv(V);
CC = [zeros(4) eye(4);-inv(M)*K -inv(M)*C];
EQ = xlsread('el_centro.xlsx');
% step = 0.02;
% time = 0:step:max(EQ(:,1));
time = EQ(:,1);
Input = EQ(:,2);
step = mean(diff(time));
y0 = [0 0 0 0 0 0 0 0]; % [d1 d2 d3 d4 v1 v2 v3 v4]
% [tsol,ysol] = ode23('testode_2D',time,y0);
[tsol,ysol] = ode23(@testode_2D,time,y0);
%% time plot
figure(1),
subplot(211),plot(time,Input)
title('El centro Base Input')
xlabel('Time')
ylabel('Displacement (m)')
subplot(212),plot(time,ysol(:,1:4))
title('El centro Floor Response')
legend('Fisrt','Second','Third','Forth')
xlabel('Time')
ylabel('Displacement (m)')
%% FRFs computations
% NFFT = 512;
% NOVERLAP = 0.5*NFFT;
% NFFT = length(time);
NFFT = length(time); % maximal frequency resolution (but implies zero overlap)
NOVERLAP = 0;
WINDOW = hanning(NFFT);
Fs = 1/step;
for ci =1:4
[Txy(:,ci),Freq] = tfestimate(Input,ysol(:,ci),WINDOW,NOVERLAP,NFFT,Fs);
end
figure(2),
loglog(Freq,abs(Txy));
title('El centro Floor Response FRF vs. Base Input')
legend('Fisrt','Second','Third','Forth')
xlabel('Frequency (Hz)')
ylabel('Displacement ratio (m / m)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy = testode_2D(t,y)
global CC EQ M
index = interp1(EQ(:,1),1:length(EQ(:,1)),t,'nearest');
if (index ~= 1) & (index ~= length(EQ(:,1)))
if t > EQ(index,1)
oran = (t - EQ(index,1))/0.02;
ag = (EQ(index+1,2) - EQ(index,2))*oran + EQ(index,2);
else
oran = (t - EQ(index-1,1))/0.02;
ag = (EQ(index,2) - EQ(index-1,2))*oran + EQ(index-1,2);
end
else
ag = EQ(index,2);
end
f1 = -M(1,1)*ag;
f2 = -M(2,2)*ag;
f3 = -M(3,3)*ag;
f4 = -M(4,4)*ag;
A00 = zeros(4);
FF = [A00;inv(M)]*[f1;f2;f3;f4];
dy = CC*y + FF;
end
5 comentarios
Mathieu NOE
el 17 de Jun. de 2021
My pleasure
read some modal analysis publications if you need to refine your skills in that matter :)
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!