FRF Issue of 4 DOF Frame
1 view (last 30 days)
Show older comments
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 Comments
Accepted Answer
Mathieu NOE
on 16 Jun 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 Comments
Mathieu NOE
on 17 Jun 2021
My pleasure
read some modal analysis publications if you need to refine your skills in that matter :)
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!