FRF Issue of 4 DOF Frame

Question

Can AKYOL 2021-6-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/856270-frf-issue-of-4-dof-frame

评论： Mathieu NOE 2021-6-17

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 个评论
显示无隐藏无

Mathieu NOE 2021-6-16

hello

we need the input data as well (el_centro file)

Can AKYOL 2021-6-16

el_centro.xlsx

Hello, here it is sir !

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2021-6-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/856270-frf-issue-of-4-dof-frame#answer_726220

在 MATLAB Online 中打开

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 个评论
显示 3更早的评论隐藏 3更早的评论

Mathieu NOE 2021-6-17

在 MATLAB Online 中打开

Hello

I changed the excitation signal to a white noise (random) , so now we can control the sampling frequency and time length; also a white noise has constant power vs frequency, which is not the case of a "true" earthquake signal; so we apply here the same techniques as in modal analysis (excitation signals are either broadband or short pulses or chirp sine)

we can get a much prettier FRF plot now :

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));
samples = 1e4;
step = 1/2;
time = (0:samples-1)'*step;
Input = randn(samples,1);
EQ = [time Input];
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('Base Input')
xlabel('Time')
ylabel('Displacement (m)')
subplot(212),plot(time,ysol(:,1:4))
title('Floor Response')
legend('Fisrt','Second','Third','Forth')
xlabel('Time')
ylabel('Displacement (m)')
%% FRFs computations
NFFT = 2048*2;
NOVERLAP = 0.5*NFFT;
% 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('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

Can AKYOL 2021-6-17

Hello, thank you again,

Now everything is much more clear

Mathieu NOE 2021-6-17

My pleasure

read some modal analysis publications if you need to refine your skills in that matter :)

请先登录，再进行评论。

FRF Issue of 4 DOF Frame

2 个评论
显示无隐藏无

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

FRF Issue of 4 DOF Frame

2 个评论 显示 无隐藏 无

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论