Integration of the acceleration signal to obtain cumulate displacement
13 views (last 30 days)
Show older comments
federico Sperati
on 18 Oct 2021
Commented: Mathieu NOE
on 5 Nov 2021
Hi everyone,
I have to integrate twice in time an accelerogram of an earthquake to obtain the diagram of the cumulative displacements of a rigid block. In particular my signal has to be analyzed only from a threshold value of the acceleration (a_y) which sets the block in motion. Also, the displacement only needs to be evaluated as long as the relative velocity between the ground and the rigid body is not zero. I have seen that the cumtrapz command can be used, however, I cannot solve my problem. Can anyone help me?
That's what i shoud re-create:
That's my actual code... i think that there are some problem with drifting factors.
clc; clear all; close all;
data=load("accelerogramma.txt");%carico i dati
time=data(:,1);%tempo in secondi
T=length(time);
acceleration=data(:,2);%accelerazione m//s2
figure
plot(time,acceleration);xlim([0 max(time)]);%grafico input
ay=0.09;
ay_plot=time*0+ay;
hold on
plot(time,ay_plot);xlim([0 max(time)]);
plot(time,zeros(1,T));xlim([0 max(time)]);
dt=abs((time(1,1)-time(2,1)));
velocity=cumtrapz(time,acceleration);
displacement=cumtrapz(time,velocity);
figure
subplot(3,1,1)
plot(time,acceleration,time,zeros(1,T));xlim([0 max(time)]);
subplot(3,1,2)
plot(time,velocity,time,zeros(1,T));xlim([0 max(time)]);
subplot(3,1,3)
plot(time,displacement,time,zeros(1,T));xlim([0 max(time)]);
acc_rel=acceleration-ay;
vel_rel=cumtrapz(time,acc_rel);
for i=1:T
if vel_rel(i)<0
vel_rel(i)=0;
else
end
end
figure
plot(time,vel_rel)
hold on
plot(time,acc_rel)
disp_rel=cumtrapz(time,vel_rel);
figure
plot(time,disp_rel);
2 Comments
Bjorn Gustavsson
on 19 Oct 2021
Are you sure you should only take accelerations larger than some absolute threshold into account? To me (not earthquake-experienced) it seems possible that you should start integrating from when the acceleration exceeds the threshold until the earthquake "quiets down" and things settle again (once the block starts sliding around it will do so until Earth stops moving about and friction stops it)?
Accepted Answer
Mathieu NOE
on 25 Oct 2021
hello Federico
check my code version :
clc; clear all; close all;
data=load("accelerogramma.txt");%carico i dati
time=data(:,1);%tempo in secondi
samples=length(time);
acceleration=data(:,2);%accelerazione m//s2
figure
plot(time,acceleration);xlim([0 max(time)]);%grafico input
ay=0.09;
ay_plot=time*0+ay;
hold on
plot(time,ay_plot);xlim([0 max(time)]);
plot(time,zeros(1,samples));xlim([0 max(time)]);
dt=mean(diff(time));
ind = find(acceleration>ay);
ind = ind(1);
acc_r = zeros(1,samples)- ay; % relative acceleration
velocity = zeros(1,samples);
a = 1;
for ci = ind:samples
acc_r(ci) = acceleration(ci) - ay;% relative acceleration
velocity(ci) = velocity(ci-1) + a*dt/2*(acc_r(ci) + acc_r(ci-1)); % iterative trapz integration
if velocity(ci)<0 % a = 0 desactivate the integration process
a = 0;
velocity(ci) = 0;
end
if acc_r(ci)>0 % a = 1 reactivate the integration process
a = 1;
end
end
displacement=cumtrapz(time,velocity);
figure
subplot(3,1,1)
plot(time,acceleration,time,zeros(1,samples));xlim([0 max(time)]);
subplot(3,1,2)
plot(time,velocity,time,zeros(1,samples));xlim([0 max(time)]);
subplot(3,1,3)
plot(time,displacement,time,zeros(1,samples));xlim([0 max(time)]);
More Answers (0)
See Also
Categories
Find more on Simulink Functions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!