.jpg.jpg)
Para resolver las ecuaciones diferenciales utilizamos uno de los metodos numéricos desarrollados para los sistemas de un grado de libertad,
function [xt,vt,at,coef,a0] = Diferencia_Central_MDOF(M,C,K,P,t,dt,x0,v0)
% Calculo de la respuesta para MDOF utilizando el procedimiento de
% diferencia central
% [xt,vt,at] = Diferencia_Central_MDOF(M,C,K,P,t,dt,u0,v0)
% M es la matriz de masas (nxn)
% C es la matriz de amortiguamientos (nxn)
% K es la matriz de rigidez (nxn)
% t es el vector de tiempo (mx1)
% dt es el intervalo de tiempo seleccionado (escalar)
% u0 son los desplazamientos iniciales de las masas (nx1)
% v0 son las velocidades iniciales de las masas (nx1)
P=[P(:,1) P]; % Replicamos la primera columna en el vector de cargas para trabajar con el instante t(-1)
a0=(P(:,2)-C*v0-K*x0)./(diag(M)); % Calculamos el vector de aceleraciones en el instante cero a artir de EQ de equilibrio
%a0=[0;0;0];
% Calculo de los coeficientes
A=M/dt^2+C/(2*dt);
A1=A^-1;
B=-K+2*M/dt^2;
C=C/(2*dt)-M/dt^2;
coef=[A B C];
% Definimos los tamanos de las matrices donde vamos a almacenar la
% informacion
[row ~]=size(M);
[len]=length(t);
xt=zeros(row,len+1);
vt=zeros(row,len+1);
at=zeros(row,len+1);
% Calculamos los valores de t(-1) y t(0)
xt(:,1)=(a0*dt^2/2)-(v0*dt)+(x0);
vt(:,1)=NaN;
at(:,1)=NaN;
xt(:,2)=x0;
vt(:,2)=v0;
at(:,2)=a0;
% Calculamos los valores de la serie paso a paso
for i=3:length(P)
% El instante n+1=i
% El instante n=i-1
% El instante n-1=i-2
xt(:,i)=(A1)*(P(:,i-1)+B*xt(:,i-1)+C*xt(:,i-2));
vt(:,i-1)=(xt(:,i)-xt(:,i-2))/(2*dt);
at(:,i-1)=(xt(:,i)-2*xt(:,i-1)+xt(:,i-2))/(dt^2);
end
xt(:,1)=[];
vt(:,1)=[];
at(:,1)=[];
MDOF Solucion de EOM acopladas
clear variables
clc
% Datos
% h = [h1 h2 ... hn]
h1=6.5; %[m] Altura de entrepiso
h2=5; %[m] Altura de entrepiso
h3=3; %[m] Altura de entrepiso
b=5; %[m] Longitud de la viga
gamma_h=2400; %[kg/m3] Peso especifico del concreto
W_lin1=156.87; %[kg/m] Peso lineal de la columna 1
W_lin2=114.03; %[kg/m] Peso lineal de la columna 2
W_lin3=77.64; %[kg/m] Peso lineal de la columna 3
Wcm=325; %[kg/m2]
e_losa=0.15; %[m]
Area=15*15; %[m2]
m_losa=e_losa*Area*gamma_h/2;
m_cm=Wcm*Area/2;
m1=2*(W_lin1*h1/2+W_lin2*h2/2);
m2=2*(W_lin2*h2/2+W_lin3*h3/2);
m3=2*(W_lin3*h3/2);
m1=m_losa+m_cm+m1;
m2=m_losa+m_cm+m2;
m3=m_losa+m_cm+m3;
E=200*10^9; %[Pa]
I1=788.710*10^6; %[mm4]
I2=457.446*10^6; %[mm4]
I3=220.325*10^6; %[mm4]
I1=I1*(1/1000)^4; %[m4]
I2=I2*(1/1000)^4; %[m4]
I3=I3*(1/1000)^4; %[m4]
k1=2*12*E*I1/h1^3;
k2=2*12*E*I2/h2^3;
k3=2*12*E*I3/h3^3;
c1=50000; %[Nm/s]
c2=50000; %[Nm/s]
c3=50000; %[Nm/s]
% Definimos la carga
F3=5000; %[N]
w=2*pi; %Frecuencia de la carga
T=2*pi/w
Calculos
S=[k1+k2 -k2 0
-k2 k2+k3 -k3
0 -k3 k3]; % Matriz de rigidez
C=[c1+c2 -c2 0
-c2 c2+c3 -c3
0 -c3 c3];
M=blkdiag(m1,m2,m3);
Procedimiento de integracion numerica
dt=0.005;
t=(0:dt:25);
P=zeros(3,length(t));
P(3,:)=F3*sin(w*t);
u0=[0 0 0]';
v0=[0,0,0]';
[xt,vt,at,~,~]=Diferencia_Central_MDOF(M,C,S,P,t,dt,u0,v0);
plot(t,xt(1,:))
hold on
plot(t,xt(2,:))
hold on
plot(t,xt(3,:))
grid on
xlabel('Tiempo [s]')
ylabel('Desplazamiento [m]')
function du2dt = odefun(t,x,params)
% Funcion para la solucion de la EQ diferencial
% 1 2 3 4 5 6 7 8 9 10 11
% params=[m1 m2 m3 k1 k2 k3 c1 c2 c3 F3 w];
c1=params(1,7);
c2=params(1,8);
c3=params(1,9);
m1=params(1,1);
m2=params(1,2);
m3=params(1,3);
k1=params(1,4);
k2=params(1,5);
k3=params(1,6);
F3=params(1,10);
w=params(1,11);
du2dt = [x(4)
x(5)
x(6)
(-(c1+c2)*x(4)+c2*x(5)-(k1+k2)*x(1)+k2*x(2))/m1
(c2*x(4)-(c2+c3)*x(5)+c3*x(6)+k2*x(1)-(k2+k3)*x(2)+k3*x(3))/m2
(F3*sin(w*t)+c3*x(5)-c3*x(6)+k3*x(2)-k3*x(3))/m3];
end
MDOF Solucion de EOM acopladas
clear variables
clc
% Datos
% h = [h1 h2 ... hn]
h1=6.5; %[m] Altura de entrepiso
h2=5; %[m] Altura de entrepiso
h3=3; %[m] Altura de entrepiso
b=5; %[m] Longitud de la viga
gamma_h=2400; %[kg/m3] Peso especifico del concreto
W_lin1=156.87; %[kg/m] Peso lineal de la columna 1
W_lin2=114.03; %[kg/m] Peso lineal de la columna 2
W_lin3=77.64; %[kg/m] Peso lineal de la columna 3
Wcm=325; %[kg/m2]
e_losa=0.15; %[m]
Area=15*15; %[m2]
m_losa=e_losa*Area*gamma_h/2;
m_cm=Wcm*Area/2;
m1=2*(W_lin1*h1/2+W_lin2*h2/2);
m2=2*(W_lin2*h2/2+W_lin3*h3/2);
m3=2*(W_lin3*h3/2);
m1=m_losa+m_cm+m1;
m2=m_losa+m_cm+m2;
m3=m_losa+m_cm+m3;
E=200*10^9; %[Pa]
I1=788.710*10^6; %[mm4]
I2=457.446*10^6; %[mm4]
I3=220.325*10^6; %[mm4]
I1=I1*(1/1000)^4; %[m4]
I2=I2*(1/1000)^4; %[m4]
I3=I3*(1/1000)^4; %[m4]
k1=2*12*E*I1/h1^3;
k2=2*12*E*I2/h2^3;
k3=2*12*E*I3/h3^3;
c1=50000; %[Nm/s]
c2=50000; %[Nm/s]
c3=50000; %[Nm/s]
% Definimos la carga
F3=5000; %[N]
w=2*pi; %Frecuencia de la carga
T=2*pi/w
Calculos
S=[k1+k2 -k2 0
-k2 k2+k3 -k3
0 -k3 k3]; % Matriz de rigidez
C=[c1+c2 -c2 0
-c2 c2+c3 -c3
0 -c3 c3];
M=blkdiag(m1,m2,m3);
Solucion con ODE45
%[t,y]=ode45(odefun,tspan,y0)
tspan=[0 25]; % Intervalo de tiempo para la solucion
y0=[0;0;0;0;0;0]; % Definimos las condiciones iniciales [u1 u2 u3 v1 v2 v3]
params=[m1 m2 m3 k1 k2 k3 c1 c2 c3 F3 w];
[t,y]=ode45(@(t,x) odefun(t,x,params),tspan,y0);
figure
plot(t,y(:,1))
hold on
plot(t,y(:,2))
hold on
plot(t,y(:,3))
xlabel('Tiempo [s]')
ylabel('Desplazamiento [m]')
set(gcf,'Visible','on','position',[0 0 1920 1080])