https://www.loom.com/share/ac7b390b536b4952be0fc42714cf5934
classdef Sismo
% Se busca hacer una clase que grabe los datos de los sismos
% En la clase A es la aceleracion
% V es la velocidad
% D es el desplazamiento
% t es el tiempo
% dt es el intervalo del registro
% Unidades corresponde a las unidades del registro, se prefiere
% trabajar en [xg cm/s cm]
% Para calcular los valores y los dibujos de peak ground hay que correr
% obj.peak_ground
% Para calcular la respuesta de un oscilador hay que correr obj.SDOF
% Para calcular el espectro con la beta de Newmark hay que poner
% obj.Espectro
properties
Name
A
V
D
t
dt
Unidades
end
% methods
% function obj=Sismo()
% obj.PGA=max(abs(max(obj.A)),abs(min(obj.A)));
% end
% end
methods
% Funcion para calcular y dibujar los valores de peak ground
function [PGA,PGV,PGD]=peak_ground(obj)
[PGA_max,idx_a]=max(obj.A);
[PGA_min,idx2_a]=min(obj.A);
PGA=max(abs(PGA_max),abs(PGA_min));
[PGV_max,idx_v]=max(obj.V);
[PGV_min,idx2_v]=min(obj.V);
PGV=max(abs(PGV_max),abs(PGV_min));
[PGD_max,idx_d]=max(obj.D);
[PGD_min,idx2_d]=min(obj.D);
PGD=max(abs(PGD_max),abs(PGD_min));
subplot(3,1,1)
plot(obj.t,obj.A)
hold on
plot(obj.t(idx_a),obj.A(idx_a),'or')
text(obj.t(idx_a)+0.25*PGA,obj.A(idx_a)+0.10*PGA,['Max = ' num2str(round(obj.A(idx_a),2))])
plot(obj.t(idx2_a),obj.A(idx2_a),'or')
text(obj.t(idx2_a)+0.25*PGA,obj.A(idx2_a)-0.10*PGA,['Min = ' num2str(round(obj.A(idx2_a),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Aceleracion [' obj.Unidades ']'])
ylim([-1.30*PGA 1.30*PGA])
subplot(3,1,2)
plot(obj.t,obj.V)
hold on
plot(obj.t(idx_v),obj.V(idx_v),'or')
text(obj.t(idx_v)+0.25*PGV,obj.V(idx_v)+0.10*PGV,['Max = ' num2str(round(obj.V(idx_v),2))])
plot(obj.t(idx2_v),obj.V(idx2_v),'or')
text(obj.t(idx2_v)+0.25*PGV,obj.V(idx2_v)-0.10*PGV,['Min = ' num2str(round(obj.V(idx2_v),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Velocidad [' obj.Unidades ']'])
ylim([-1.30*PGV 1.30*PGV])
subplot(3,1,3)
plot(obj.t,obj.D)
hold on
plot(obj.t(idx_d),obj.D(idx_d),'or')
text(obj.t(idx_d)+0.25*PGD,obj.D(idx_d)+0.10*PGD,['Max = ' num2str(round(obj.D(idx_d),2))])
plot(obj.t(idx2_d),obj.D(idx2_d),'or')
text(obj.t(idx2_d)+0.25*PGD,obj.D(idx2_d)-0.10*PGD,['Min = ' num2str(round(obj.D(idx2_d),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Desplazamiento [' obj.Unidades ']'])
ylim([-1.30*PGD 1.30*PGD])
sgtitle({obj.Name},'Interpreter','none')
set(gcf,'Visible','on','position',[0 0 1920 1080])
end
end
methods
% Respuesta de un SDOF para el registro
function [xt,vt,at]=SDOF(obj)
% Funcion para calcular la respuesta en el tiempo de un SDOF
% con el registro de aceleracion
% [xt,vt,at,coef,Tabla,Check] = Newmark(m,k,xi,t,P,delta_t,x0,v0,beta,gamma)
Tn=input('Cual es el valor de Tn = ');
xi=input('Cual es el valor de xi = ');
beta=input('Cual es el valor de beta (1/4 o 1/6) = ');
gamma=input('Cual es el valor de gamma (1/2) = ');
wn=2*pi/Tn;
m=1;
k=wn^2*m;
Peff=-m*obj.A;
[xt,vt,at,~,~,~] = Newmark(m,k,xi,obj.t,Peff,obj.dt,0,0,beta,gamma);
[at_max,idx_a]=max(at);
[at_min,idx2_a]=min(at);
ao=max(abs(at_max),abs(at_min));
[vt_max,idx_v]=max(vt);
[vt_min,idx2_v]=min(vt);
vo=max(abs(vt_max),abs(vt_min));
[xt_max,idx_d]=max(xt);
[xt_min,idx2_d]=min(xt);
xo=max(abs(xt_max),abs(xt_min));
subplot(3,1,1)
plot(obj.t,xt)
hold on
plot(obj.t(idx_d),xt(idx_d),'or')
text(obj.t(idx_d)+0.25*xo,xt(idx_d)+0.10*xo,['Max = ' num2str(round(xt(idx_d),2))])
plot(obj.t(idx2_d),xt(idx2_d),'or')
text(obj.t(idx2_d)+0.25*xo,xt(idx2_d)-0.10*xo,['Min = ' num2str(round(xt(idx2_d),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Desplazamiento [' obj.Unidades ']'])
ylim([-1.30*xo 1.30*xo])
subplot(3,1,2)
plot(obj.t,vt)
hold on
plot(obj.t(idx_v),vt(idx_v),'or')
text(obj.t(idx_v)+0.25*vo,vt(idx_v)+0.10*vo,['Max = ' num2str(round(vt(idx_v),2))])
plot(obj.t(idx2_v),vt(idx2_v),'or')
text(obj.t(idx2_v)+0.25*vo,vt(idx2_v)-0.10*vo,['Min = ' num2str(round(vt(idx2_v),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Velocidad [' obj.Unidades ']'])
ylim([-1.30*vo 1.30*vo])
subplot(3,1,3)
plot(obj.t,at)
hold on
plot(obj.t(idx_a),at(idx_a),'or')
text(obj.t(idx_a)+0.25*ao,at(idx_a)+0.10*ao,['Max = ' num2str(round(at(idx_a),2))])
plot(obj.t(idx2_a),at(idx2_a),'or')
text(obj.t(idx2_a)+0.25*ao,at(idx2_a)-0.10*ao,['Min = ' num2str(round(at(idx2_a),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Aceleracion [' obj.Unidades ']'])
ylim([-1.30*ao 1.30*ao])
sgtitle({obj.Name},'Interpreter','none')
set(gcf,'Visible','on','position',[0 0 1920 1080])
end
end
methods
% Funcion para calcular y dibujar el espectro respuesta
function [Sd,Sv,Sa,Pv,Pa,Tn]=Espectro(obj)
Tn=input('Cual es el rango de Tn = ');
xi=input('Cual es el valor de xi = ');
beta=input('Cual es el valor de beta (1/4 o 1/6) = ');
gamma=input('Cual es el valor de gamma (1/2) = ');
% [Sd,Sv,Sa,Pv,Pa] = Espectro_lineal_BN(E,Tn,xi,beta,gamma)
[Sd,Sv,Sa,Pv,Pa] = Espectro_lineal_BN(obj,Tn,xi,beta,gamma);
subplot(3,1,1)
plot(Tn,Sd)
grid on
xlabel('Periodo [s]')
ylabel(['Desplazamiento [' obj.Unidades ']'])
subplot(3,1,2)
plot(Tn,Pv)
grid on
xlabel('Periodo [s]')
ylabel(['Pseudo Velocidad [' obj.Unidades ']'])
subplot(3,1,3)
plot(Tn,Pa)
grid on
xlabel('Periodo [s]')
ylabel(['Pseudo Aceleracion [' obj.Unidades ']'])
set(gcf,'Visible','on','position',[0 0 1920 1080])
end
end
end
Funcion para el calculo de la respuesta dinamica utilizando el metodo de interpolacion de la exitacion (Time Stepping Method).
[xt,vt,at,coef,Tabla] = Newmark(m,k,xi,t,P,delta_t,x0,v0,beta,gamma)
Aceleracion constante - beta=1/4 y gamma=1/2
Aceleracion lineal - beta=1/6 y gamma=1/2
function [xt,vt,at,coef,Tabla,Check] = Newmark(m,k,xi,t,P,delta_t,x0,v0,beta,gamma)
wn=(k/m)^0.5;
Tn=2*pi/wn;
% Calculos iniciales
c=xi*(2*m*wn);
P0=P(1,1);
a0=(P0-c*v0-k*x0)/(m);
a1=((1)/(beta*delta_t^2))*m+((gamma)/(beta*delta_t))*c;
a2=(1/(beta*delta_t))*m+(gamma/beta-1)*c;
a3=(1/(2*beta)-1)*m+delta_t*(gamma/(2*beta)-1)*c;
k_hat=k+a1;
coef=[a1 a2 a3 k_hat]';
% 1 2 3 4 5 6
% Tabla = [ti Pi ui vi ai pi+1]
Tabla=zeros(length(t),6);
Tabla(:,1)=t;
Tabla(:,2)=P;
Tabla(1,3)=x0;
Tabla(1,4)=v0;
Tabla(1,5)=a0;
% Tenemos que resolver la ecuacion
for i=2:length(t)
p_hat_1=Tabla(i,2)+a1*Tabla(i-1,3)+a2*Tabla(i-1,4)+a3*Tabla(i-1,5);
Tabla(i,6)=p_hat_1;
Tabla(i,3)=p_hat_1/k_hat;
Tabla(i,4)=(gamma/(beta*delta_t))*(Tabla(i,3)-Tabla(i-1,3))+(1-gamma/beta)*Tabla(i-1,4)+delta_t*(1-gamma/(2*beta))*Tabla(i-1,5);
Tabla(i,5)=(1/(beta*delta_t^2))*(Tabla(i,3)-Tabla(i-1,3))-(1/(beta*delta_t))*Tabla(i-1,4)-(1/(2*beta)-1)*Tabla(i-1,5);
end
xt=Tabla(:,3);
vt=Tabla(:,4);
at=Tabla(:,5)-P/m; % Corregimos para que nos de la aceleracion absoluta
% Verificacion de estabilidad
Lim=(1/(pi*sqrt(2)))*(1/(beta*sqrt(gamma-2*beta)));
if delta_t/Tn < Lim
Check="Ok";
else
Check="Not OK";
end
function [Sd,Sv,Sa,Pv,Pa] = Espectro_lineal_BN(E,Tn,xi,beta,gamma)
% Calculo del espectro respuesta lineal utilizando el procedimiento de la
% beta de Newmark
% [Sd,Sa,Sv,Pv,Pa] = Espectro_lineal_BN(E,Tn,xi,beta,gamma)
% Sd → Espectro relativo de desplazamientos
% Sv → Espectro relativo de velocidades
% Sa → Espectro relativo de aceleracion
% Pv → Pseudo espectro de velocidades
% Pa → Pseudo espectro de aceleraciones
% E → clase que contenga el acelerograma
% Tn → Rango de periodos evaluados
% xi → porcentaje del amortiguamiento critico
% beta y gamma → constantes del metodo
% beta = 1/4 (aceleracion constante) - 1/6 (aceleracion lineal)
% gamma = 1/2 (aceleracion constante) - 1/2 (aceleracion constante)
m=1;
wn=2*pi./Tn;
k=(wn).^2*m;
x0=0;
v0=0;
Peff=-m*E.A;
Sd=zeros(1,length(Tn));
Sv=zeros(1,length(Tn));
Sa=zeros(1,length(Tn));
for i=1:length(Tn)
[xt,vt,at] = Newmark(m,k(1,i),xi,E.t,Peff,E.dt,x0,v0,beta,gamma);
Sd(1,i)=max(xt);
Sv(1,i)=max(vt);
Sa(1,i)=max(at);
end
% Calculo de los Pseudo-Espectros
Pv=wn.*Sd;
Pa=wn.^2.*Sd;
end
function [Reg] = Leer_Registro_NGA(filename,headerlinesIn)
% Leer registros del NGA
% La funcion detecta automaticamente el dt y las unidades.
% filename = es el nombre del archivo con la extension.txt
% headerlinesIn = lee la data numerica a partir de la linea indicada + 1
Temp=importdata(filename,' ',headerlinesIn);
% Pasamos la estructura a matrices
Reg_temp=Temp.data; % Matrix con datos del acelerograma
desc=upper(Temp.textdata); % Matriz con descripcion del acelerograma, la ponemos en mayuscula para poder hacer los algoritmos de busqueda de strings
% Generamos un algoritmo para buscar el timestep del registro
pat=("0.001"|"0.005"|"0.01"|"0.02"|"0.05"|"0.0024");
idx=contains(desc,pat);
% Extraemos el valor
dt=cell2mat(extract(desc(idx),pat));
% Definimos los parametros del acelerograma dentro de la clase "Sismos"
Reg=Sismo;
Reg.Name=filename;
% Grabar el valor de dt en valor numerico
Reg.dt=str2double(dt);
% Generamos un algoritmo para leer las unidades del registro
% Generamos un algoritmo para buscar el timestep del registro
pat2=("UNITS OF G"|"UNITS: G"|"MM/SEC/SEC"|"CM/SEC/SEC");
idx2=contains(desc,pat2);
% Extraemos el valor
Reg.Unidades=cell2mat(extract(desc(idx2),pat2));
clear idx idx2 len pat pat2 dt desc
% Calculamos el tamaño de la tabla del registro [filas columnas]
len=size(Reg_temp);
% determinamos el vector de tiempo
Reg.t=(0:Reg.dt:(len(1,1)*len(1,2)-1)*Reg.dt)';
% Determinamos el tamaño del vector de aceleracion
Reg.A=zeros(len(1,1)*len(1,2),1);
% Ordenamos el vector de aceleracion
n=1;
for i=1:len(1,1)
Reg.A(n:(n-1)+len(1,2),1)=Reg_temp(i,:)';
n=n+len(1,2);
end
% Reemplazamos todos los NaN por cero
Reg.A(isnan(Reg.A))=0;
% Calculamos el vector de aceleracion y desplazamientos
Reg.V=cumtrapz(Reg.dt,Reg.A);
Reg.D=cumtrapz(Reg.dt,Reg.V);
clear filename i n len Reg_temp Temp headerlinesIn
end
function [Reg] = Leer_Registro(filename,headerlinesIn,dt,Unidades)
% Leer registros del NGA
% La funcion detecta automaticamente el dt y las unidades.
% filename = es el nombre del archivo con la extension.txt
% headerlinesIn = lee la data numerica a partir de la linea indicada + 1
% dt es el timestep del registro
% Unidades es un string con las unidades del registro
Temp=importdata(filename,' ',headerlinesIn);
% Pasamos la estructura a matrices
Reg_temp=Temp.data; % Matrix con datos del acelerograma
% Definimos los parametros del acelerograma dentro de la clase "Sismos"
Reg=Sismo;
Reg.Name=filename;
Reg.dt=dt;
Reg.Unidades=Unidades;
% Calculamos el tamaño de la tabla del registro [filas columnas]
len=size(Reg_temp);
% determinamos el vector de tiempo
Reg.t=(0:Reg.dt:(len(1,1)*len(1,2)-1)*Reg.dt)';
% Determinamos el tamaño del vector de aceleracion
Reg.A=zeros(len(1,1)*len(1,2),1);
% Ordenamos el vector de aceleracion
n=1;
for i=1:len(1,1)
Reg.A(n:(n-1)+len(1,2),1)=Reg_temp(i,:)';
n=n+len(1,2);
end
% Reemplazamos todos los NaN por cero
Reg.A(isnan(Reg.A))=0;
% Calculamos el vector de aceleracion y desplazamientos
Reg.V=cumtrapz(Reg.dt,Reg.A);
Reg.D=cumtrapz(Reg.dt,Reg.V);
end
Para llevar el registro de aceleraciones a un registro de velocidades y desplazamientos debemos realizar una integración numérica
Para organizar la información de una mejor manera se opto por realizar una clase en el Matlab.
La clase nos va a permitir grabar propiedades especificas del registro en una variable y asi mantener la información ordenada.
classdef Sismo
% Se busca hacer una clase que grabe los datos de los sismos
% En la clase A es la aceleracion
% V es la velocidad
% D es el desplazamiento
% t es el tiempo
% dt es el intervalo del registro
% Unidades corresponde a las unidades del registro, se prefiere
% trabajar en [xg cm/s cm]
% Para calcular los valores y los dibujos de peak ground hay que correr
% obj.peak_ground
% Para calcular la respuesta de un oscilador hay que correr obj.SDOF
% Para calcular el espectro con la beta de Newmark hay que poner
% obj.Espectro
properties
Name
A
V
D
t
dt
Unidades
PGA
PGV
PGD
end
methods
% Funcion para calcular y dibujar los valores de peak ground
function [PGA,PGV,PGD]=peak_ground(obj)
[PGA_max,idx_a]=max(obj.A);
[PGA_min,idx2_a]=min(obj.A);
PGA=max(abs(PGA_max),abs(PGA_min));
[PGV_max,idx_v]=max(obj.V);
[PGV_min,idx2_v]=min(obj.V);
PGV=max(abs(PGV_max),abs(PGV_min));
[PGD_max,idx_d]=max(obj.D);
[PGD_min,idx2_d]=min(obj.D);
PGD=max(abs(PGD_max),abs(PGD_min));
subplot(3,1,1)
plot(obj.t,obj.A)
hold on
plot(obj.t(idx_a),obj.A(idx_a),'or')
text(obj.t(idx_a)+0.25*PGA,obj.A(idx_a)+0.10*PGA,['Max = ' num2str(round(obj.A(idx_a),2))])
plot(obj.t(idx2_a),obj.A(idx2_a),'or')
text(obj.t(idx2_a)+0.25*PGA,obj.A(idx2_a)-0.10*PGA,['Min = ' num2str(round(obj.A(idx2_a),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Aceleracion [' obj.Unidades ']'])
ylim([-1.30*PGA 1.30*PGA])
subplot(3,1,2)
plot(obj.t,obj.V)
hold on
plot(obj.t(idx_v),obj.V(idx_v),'or')
text(obj.t(idx_v)+0.25*PGV,obj.V(idx_v)+0.10*PGV,['Max = ' num2str(round(obj.V(idx_v),2))])
plot(obj.t(idx2_v),obj.V(idx2_v),'or')
text(obj.t(idx2_v)+0.25*PGV,obj.V(idx2_v)-0.10*PGV,['Min = ' num2str(round(obj.V(idx2_v),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Velocidad [' obj.Unidades ']'])
ylim([-1.30*PGV 1.30*PGV])
subplot(3,1,3)
plot(obj.t,obj.D)
hold on
plot(obj.t(idx_d),obj.D(idx_d),'or')
text(obj.t(idx_d)+0.25*PGD,obj.D(idx_d)+0.10*PGD,['Max = ' num2str(round(obj.D(idx_d),2))])
plot(obj.t(idx2_d),obj.D(idx2_d),'or')
text(obj.t(idx2_d)+0.25*PGD,obj.D(idx2_d)-0.10*PGD,['Min = ' num2str(round(obj.D(idx2_d),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Desplazamiento [' obj.Unidades ']'])
ylim([-1.30*PGD 1.30*PGD])
sgtitle({obj.Name},'Interpreter','none')
set(gcf,'Visible','on','position',[0 0 1920 1080])
end
end
methods
function [xt,vt,at]=SDOF(obj)
% Funcion para calcular la respuesta en el tiempo de un SDOF
% con el registro de aceleracion
% [xt,vt,at,coef,Tabla,Check] = Newmark(m,k,xi,t,P,delta_t,x0,v0,beta,gamma)
Tn=input('Cual es el valor de Tn = ');
xi=input('Cual es el valor de xi = ');
beta=input('Cual es el valor de beta (1/4 o 1/6) = ');
gamma=input('Cual es el valor de gamma (1/2) = ');
wn=2*pi/Tn;
m=1;
k=wn^2*m;
Peff=-m*obj.A;
[xt,vt,at,~,~,~] = Newmark(m,k,xi,obj.t,Peff,obj.dt,0,0,beta,gamma);
[at_max,idx_a]=max(at);
[at_min,idx2_a]=min(at);
ao=max(abs(at_max),abs(at_min));
[vt_max,idx_v]=max(vt);
[vt_min,idx2_v]=min(vt);
vo=max(abs(vt_max),abs(vt_min));
[xt_max,idx_d]=max(xt);
[xt_min,idx2_d]=min(xt);
xo=max(abs(xt_max),abs(xt_min));
subplot(3,1,1)
plot(obj.t,xt)
hold on
plot(obj.t(idx_d),xt(idx_d),'or')
text(obj.t(idx_d)+0.25*xo,xt(idx_d)+0.10*xo,['Max = ' num2str(round(xt(idx_d),2))])
plot(obj.t(idx2_d),xt(idx2_d),'or')
text(obj.t(idx2_d)+0.25*xo,xt(idx2_d)-0.10*xo,['Min = ' num2str(round(xt(idx2_d),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Desplazamiento [' obj.Unidades ']'])
ylim([-1.30*xo 1.30*xo])
subplot(3,1,2)
plot(obj.t,vt)
hold on
plot(obj.t(idx_v),vt(idx_v),'or')
text(obj.t(idx_v)+0.25*vo,vt(idx_v)+0.10*vo,['Max = ' num2str(round(vt(idx_v),2))])
plot(obj.t(idx2_v),vt(idx2_v),'or')
text(obj.t(idx2_v)+0.25*vo,vt(idx2_v)-0.10*vo,['Min = ' num2str(round(vt(idx2_v),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Velocidad [' obj.Unidades ']'])
ylim([-1.30*vo 1.30*vo])
subplot(3,1,3)
plot(obj.t,at)
hold on
plot(obj.t(idx_a),at(idx_a),'or')
text(obj.t(idx_a)+0.25*ao,at(idx_a)+0.10*ao,['Max = ' num2str(round(at(idx_a),2))])
plot(obj.t(idx2_a),at(idx2_a),'or')
text(obj.t(idx2_a)+0.25*ao,at(idx2_a)-0.10*ao,['Min = ' num2str(round(at(idx2_a),2))])
grid on
xlabel('Tiempo [s]')
ylabel(['Aceleracion [' obj.Unidades ']'])
ylim([-1.30*ao 1.30*ao])
sgtitle({obj.Name},'Interpreter','none')
set(gcf,'Visible','on','position',[0 0 1920 1080])
end
end
methods
% Funcion para calcular y dibujar el espectro respuesta
function [Sd,Sv,Sa,Pv,Pa,Tn]=Espectro(obj)
Tn=input('Cual es el rango de Tn = ');
xi=input('Cual es el valor de xi = ');
beta=input('Cual es el valor de beta (1/4 o 1/6) = ');
gamma=input('Cual es el valor de gamma (1/2) = ');
% [Sd,Sv,Sa,Pv,Pa] = Espectro_lineal_BN(E,Tn,xi,beta,gamma)
[Sd,Sv,Sa,Pv,Pa] = Espectro_lineal_BN(obj,Tn,xi,beta,gamma);
subplot(3,1,1)
plot(Tn,Sd)
grid on
xlabel('Periodo [s]')
ylabel(['Desplazamiento [' obj.Unidades ']'])
subplot(3,1,2)
plot(Tn,Pv)
grid on
xlabel('Periodo [s]')
ylabel(['Pseudo Velocidad [' obj.Unidades ']'])
subplot(3,1,3)
plot(Tn,Pa)
grid on
xlabel('Periodo [s]')
ylabel(['Pseudo Aceleracion [' obj.Unidades ']'])
set(gcf,'Visible','on','position',[0 0 1920 1080])
end
end
end
A continuación se presentan los códigos utilizados para cargar los registros de aceleracion.
<aside> 💡 Este código sirve para los registros NGA
</aside>