Integration of the acceleration signal to obtain cumulate displacement

14 visualizaciones (últimos 30 días)
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 comentarios
Bjorn Gustavsson
Bjorn Gustavsson el 19 de Oct. de 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)?
federico Sperati
federico Sperati el 19 de Oct. de 2021
Hi Bjorn, thanks for your answer, I need to integrate the acceleration only when it's bigger than my threshold value because only in that case my rigid block starts up.
The real problem is that i need to continue the integration until the value of the relative velocity between block and ground is greater than zero.
In this period (from the moment in which a(t)<0 till the relative velocity of the block is not null) the block is in a decelerated motion and I need to compute the linked displacement.
It's called Newmark Method

Iniciar sesión para comentar.

Respuesta aceptada

Mathieu NOE
Mathieu NOE el 25 de Oct. de 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)]);

Más respuestas (0)

Categorías

Más información sobre Simulink Functions en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by