How to use 'trapz' command in a foor loop

5 visualizaciones (últimos 30 días)
Sumit Saha
Sumit Saha el 9 de Mayo de 2021
Comentada: Star Strider el 11 de Mayo de 2021
main file:
% Analysis of response spectra for given ground motion acceleration data
clear all
close all
clc
% Supershear & Subshear Earthquake Data
load Histories/IND.X0kY-1k.CXY.sema;% Ground-motion acceleration data
load Histories/IND.X0kY-2k.CXY.sema;% Ground-motion acceleration data
ground_acc(:,:,1) = IND_X0kY_1k_CXY;
ground_acc(:,:,2) = IND_X0kY_2k_CXY;
T_D = 8 ; % duration of earthquake in secs
dt = 0.002; % dt = sampling interval in seconds
xi = 0.02; % xi = ratio of critical damping
T_series_EQ = [0.002:dt:T_D]'; % Time column for EQ data
g =9.81 ; % accleration due to gravity in m/sec^2
% time period in vector form for response spectrum plot
Time_Period_Spectrum = [0.01:0.01:10];
for j=1:2
[PSA(:,:,j), PSV(:,:,j), SD(:,:,j), SA(:,:,j), SV(:,:,j), OUT(:,:,j),vel(:,:,j)] = responseSpectra(xi, Time_Period_Spectrum, ground_acc(:,:,j), dt);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Different Intensity Measure Calculation for EQ DATA %%%%%%%%%%%%%%%%
%%% Reference : Ground motion intensity measures for seismic probabilistic risk analysis
%% Peak Based Intensity Measure Calculation
for p=1:2
Peak_Ground_Velocity(p) = max(abs(vel(:,:,p)));
Peak_Ground_Accleration(p) = max(abs(ground_acc(:,:,p)));
%%% Pseudo-spectral accleration at fundamental period
T_f_structure = 1.42; % Fundamental Period of structure in sec
PSA_TF(p) = interp1 (Time_Period_Spectrum,PSA(:,:,p),T_f_structure);
%%% Spectral Intensity (related to kinetic energy stored in the structure During Earthquake)
S_v(:,p) = PSV(:,:,p)';
T = Time_Period_Spectrum';
Trange_SI = (T >= 0.1) & (T <= 2.5); %-period interval of Structures
I_H(p) = trapz(T(Trange_SI), S_v(:,p,(Trange_SI))); %Spectral_Intensity
end
Function file: given below
How to use Trapz command in case of loop for multiple eq data?
  13 comentarios
Star Strider
Star Strider el 10 de Mayo de 2021
In the context of the code, what are ‘1.txt’ and ‘2.txt’?
They are vectors, however everything else seem to be matrices.
Sumit Saha
Sumit Saha el 10 de Mayo de 2021
@Star delete these two lines...write load 1.txt && load 2.txt
% Supershear & Subshear Earthquake Data
load Histories/IND.X0kY-1k.CXY.sema;% Ground-motion acceleration data
load Histories/IND.X0kY-2k.CXY.sema;% Ground-motion acceleration data

Iniciar sesión para comentar.

Respuesta aceptada

Star Strider
Star Strider el 10 de Mayo de 2021
The problem is that ‘S_v’ has only one non-singleton dimension (i.e. it is a vector).
That is throwing the error.
Changing the assignment to —
I_H(p) = trapz(T(Trange_SI), S_v((Trange_SI))); %Spectral_Intensity
no longer throws the error, however the values fo ‘I_H’ aare both identical.
I have no idea what you are doing, so I will let you sort that out.
I ran this in the Online Run Feature—
ground_acc(:,:,1) = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/613135/1.txt');
ground_acc(:,:,2) = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/613140/2.txt');
T_D = 8 ; % duration of earthquake in secs
dt = 0.002; % dt = sampling interval in seconds
xi = 0.02; % xi = ratio of critical damping
T_series_EQ = [0.002:dt:T_D]'; % Time column for EQ data
g =9.81 ; % accleration due to gravity in m/sec^2
% time period in vector form for response spectrum plot
Time_Period_Spectrum = [0.01:0.01:10];
for j=1:2
[PSA(:,:,j), PSV(:,:,j), SD(:,:,j), SA(:,:,j), SV(:,:,j), OUT(:,:,j),vel(:,:,j)] = responseSpectra(xi, Time_Period_Spectrum, ground_acc(:,:,j), dt);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Different Intensity Measure Calculation for EQ DATA %%%%%%%%%%%%%%%%
%%% Reference : Ground motion intensity measures for seismic probabilistic risk analysis
%% Peak Based Intensity Measure Calculation
for p=1:2
Peak_Ground_Velocity(p) = max(abs(vel(:,:,p)));
Peak_Ground_Accleration(p) = max(abs(ground_acc(:,:,p)));
%%% Pseudo-spectral accleration at fundamental period
T_f_structure = 1.42; % Fundamental Period of structure in sec
PSA_TF(p) = interp1 (Time_Period_Spectrum,PSA(:,:,p),T_f_structure);
%%% Spectral Intensity (related to kinetic energy stored in the structure During Earthquake)
S_v(:,p) = PSV(:,:,p)';
T = Time_Period_Spectrum';
Trange_SI = (T >= 0.1) & (T <= 2.5); %-period interval of Structures
% Qp = p
% QTrange_SI = Trange_SI
% QsizeT = size(T)
% QsizeS_v = size(S_v)
% QT = T(Trange_SI)
% % QS_v = S_v(:,p,(Trange_SI))
% QS_v = S_v((Trange_SI))
% % I_H(p) = trapz(T(Trange_SI), S_v(:,p,(Trange_SI))); %Spectral_Intensity
I_H(p) = trapz(T(Trange_SI), S_v((Trange_SI))); %Spectral_Intensity
QI_H = sprintf('I_H(%d) = %23.15E',p,I_H(p))
end
QI_H = 'I_H(1) = 4.013137327351848E-01'
QI_H = 'I_H(2) = 4.013137327351848E-01'
function [PSA, PSV, SD, SA, SV, OUT,vel] = responseSpectra(xi, Time_Period_Spectrum, ground_acc, dt)
vel = cumtrapz(ground_acc)*dt;
disp = cumtrapz(vel)*dt;
% Spectral solution
for i = 1:length(Time_Period_Spectrum)
omega_w = 2*pi/Time_Period_Spectrum(i);
C = 2*xi*omega_w;
K = omega_w^2;
y(:,1) = [0;0];
A = [0 1; -K -C];
Ae = expm(A*dt); AeB = A\(Ae-eye(2))*[0;1];
for k = 2:numel(ground_acc)
y(:,k) = Ae*y(:,k-1) + AeB*ground_acc(k);
end
displ = (y(1,:))'; % Relative displacement vector (cm)
veloc = (y(2,:))'; % Relative velocity (cm/s)
foverm = omega_w^2*displ; % Lateral resisting force over mass (cm/s^2)
absacc = -2*xi*omega_w*veloc-foverm; % Absolute acceleration from equilibrium (cm/s^2)
% Extract peak values
displ_max(i) = max(abs(displ)); % Spectral relative displacement (cm)
veloc_max(i) = max(abs(veloc)); % Spectral relative velocity (cm/s)
absacc_max(i) = max(abs(absacc)); % Spectral absolute acceleration (cm/s^2)
foverm_max(i) = max(abs(foverm)); % Spectral value of lateral resisting force over mass (cm/s^2)
pseudo_acc_max(i) = displ_max(i)*omega_w^2; % Pseudo spectral acceleration (cm/s)
pseudo_veloc_max(i) = displ_max(i)*omega_w; % Pseudo spectral velocity (cm/s)
PSA(i) = pseudo_acc_max(i); % PSA (cm/s^2)
SA(i) = absacc_max(i); % SA (cm/s^2)
PSV(i) = pseudo_veloc_max(i); % PSV (cm/s)
SV(i) = veloc_max(i); % SV (cm/s)
SD(i) = displ_max(i); % SD (cm)
% Time series of acceleration, velocity and displacement response of SDF system
OUT.acc(:,i) = absacc;
OUT.vel(:,i) = veloc;
OUT.disp(:,i) = displ;
end
end
I kept the debugging steps in the code, although commented-out. Delete them if you no longer need them, or un-comment or change them if there are still problems and you want to see the interim results.
  6 comentarios
Sumit Saha
Sumit Saha el 10 de Mayo de 2021
@Star Strider thanks a lot.
Star Strider
Star Strider el 11 de Mayo de 2021
As always, my pleasure!

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Earthquake Engineering en Help Center y File Exchange.

Productos


Versión

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by