Sistema resorte-masa

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.

1 Introducción

El caso considerado está formado por tres muelles y dos masas, unidos a dos paredes separadas 9 metros, que se deslizan sin rozamiento sobre un plano horizontal.

2 Ecuaciones diferenciales

La solución del un sistema resorte-masa se obtiene de resolver las ecuaciones diferenciales obtenidas al aplicar la segunda ley de Newton a las dos masas (en adelante m1 y m2). Los grados de libertad usados son x e y, siendo estos la distancia de cada masa a su posición de equilibrio. Con este criterio las ecuaciones obtenidas son las siguientes: [math]\left\{\begin{matrix}m_{1}\ddot x=k_{2}(y-x)-k_{1}x\\m_{2}\ddot y=-k_{3}y-k_{2}(y-x)\end{matrix}\right.[/math]

3 Resolucion

Para resolver este sistema de ecuaciones diferenciales se van a realizar dos aproximaciones con métodos numéricos iterativos en MATLAB:

• Método de Newmark para sistemas de orden 2
• Método de Runge-Kutta de cuarto orden

Los valores constantes usados para la resolución son los siguientes:: [math]m_{1}=2~[kg][/math]: [math]m_{2}=1~[kg][/math]: [math]k_{1}=4~[N/m][/math]: [math]k_{2}=2~[N/m][/math]: [math]k_{3}=1~[N/m][/math]:

3.1 Método de Newmark

Para poder aplicar el metodo es necesario reducir el orden de la ecuación en uno, de manera que realizamos el cambio de variable:: [math]\dot x=u[/math]: [math]\dot y=v[/math] Una vez realizado el cambio de variable el sistema resultante es:: [math]m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y[/math]: [math]m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y[/math]: [math]\dot x=u[/math]: [math]\dot y=v[/math] Con las condiciones iniciales [math]x_{0}=1~[m][/math] e [math]y_{0}=1,5~[m][/math] el código en MATLAB del método es:

% Sistema de 2 masas y 3 muelles
clear all
 
% Intervalo de tiempo
t0=0;
tN=10;
 
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
 
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;
 
% Parametros Newmark
beta=1/4;
gamma=1/2;
 
 
N=1000; % Pasos
h=(tN-t0)/N; % Intervalo
 
% Primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
u(1)=u0;
v(1)=v0;
 
% Vector solucion S
S=[x(1);y(1);u(1);v(1)];
 
% Metodo de Newmark
 
for n=1:N
    % Matrices de coeficientes
    A=[1+beta*((k1+k2)/m1)*h^2 , -beta*(k2/m1)*h^2 , 0 , 0 ; ...
        -beta*(k2/m2)*h^2 , 1+beta*((k2+k3)/m2)*h^2 , 0 , 0 ; ...
        h*gamma*((k2+k1)/m1) , -h*gamma*(k2/m1) , 1 , 0 ; ...
        -h*gamma*(k2/m2) , h*gamma*((k2+k3)/m2) , 0 , 1];
    
    B=[1-(1/2-beta)*((k2+k1)/m1)*h^2, (1/2-beta)*(k2/m1)*h^2 , h , 0 ; ...
        (1/2-beta)*(k2/m2)*h^2 , 1-(1/2-beta)*((k2+k3)/m2)*h^2 , 0 , h ; ...
        -h*(1-gamma)*((k2+k1)/m1) , h*(1-gamma)*(k2/m1) , 1 , 0 ; ...
        h*(1-gamma)*(k2/m2) , -h*(1-gamma)*((k3+k2)/m2) , 0 , 1];

    % Sistema
    S=A\(B*S);
    x(n+1)=S(1);
    y(n+1)=S(2);
    u(n+1)=S(3);
    v(n+1)=S(4);
end


% Representacion
t=t0:h:tN;
figure(1)
axis([0,9,0,10])
hold on
title('Newmark')
plot(x+2,t,'-')
plot(y+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off

centro

3.2 Método de Runge-Kutta

Al igual que con el metodo de Newmark, para aplicar Runge-Kutta al problema es necesario reducir el orden de la ecuacion, por lo que realizamos el mismo cambio de variable y, por consiguiente, obtenemos el mismo sistema de ecuaciones:: [math]m_{1}\dot x=-(k_{1}+k_{2})x+k_{2}y[/math]: [math]m_{2}\dot v=k_{2}x-(k_{2}+k_{3})y[/math]: [math]\dot x=u[/math]: [math]\dot y=v[/math] El código en MATLAB resulta:

% Sistema de 2 masas y 3 muelles
clear all
% Intervalo de tiempo
t0=0;
tN=10;
% Datos
k1=4;
k2=2;
k3=1;
m1=2;
m2=1;
%Condiciones iniciales
x0=1;
y0=1.5;
u0=0;
v0=0;

N=1000; % Pasos
h=(tN-t0)/N; % Intervalo


% Primeros valores
xx(1)=x0;
yy(1)=y0;
uu(1)=u0;
vv(1)=v0;

% Metodo de Runge-Kutta

for n=1:N
    % Primera K
    k1x=uu(n);
    k1y=vv(n);
    k1u=-(k2+k1)/m1*xx(n)+k2/m1*yy(n);
    k1v=k2/m2*xx(n)-((k2+k3)/m2)*yy(n);
    % Segunda K
    k2x=uu(n)+1/2*h;
    k2y=vv(n)+1/2*h;
    k2u=-(k2+k1)/m1*(xx(n)+1/2*h)+k2/m1*(yy(n)+1/2*k1u*h);
    k2v=k2/m2*(xx(n)+1/2*h)-((k2+k3)/m2)*(yy(n)+1/2*k1v*h);

    % Tercera K
    k3x=uu(n)+1/2*h;
    k3y=vv(n)+1/2*h;
    k3u=-(k2+k1)/m1*(xx(n)+1/2*h)+k2/m1*(yy(n)+1/2*k2u*h);
    k3v=k2/m2*(xx(n)+1/2*h)-((k2+k3)/m2)*(yy(n)+1/2*k2v*h);

    % Cuarta K
    k4x=uu(n)+h;
    k4y=vv(n)+h;
    k4u=-(k2+k1)/m1*(xx(n)+h)+k2/m1*(yy(n)+k3u*h);
    k4v=k2/m2*(xx(n)+h)-((k2+k3)/m2)*(yy(n)+k3v*h);

    xx(n+1)=xx(n)+(h/6)*(k1x+2*k2x+2*k3x+k4x);
    yy(n+1)=yy(n)+(h/6)*(k1y+2*k2y+2*k3y+k4y);
    uu(n+1)=uu(n)+(h/6)*(k1u+2*k2u+2*k3u+k4u);
    vv(n+1)=vv(n)+(h/6)*(k1v+2*k2v+2*k3v+k4v);
end

axis([0,9,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(yy+4,t,'-r')
plot(2,t,'-g')
plot(4,t,'-g')
legend('Masa m1','Masa m2','Posicion de equilibrio')
hold off

centro

3.3 Efecto de la velocidad inicial en la trayectoria

En el caso anterior ambas masas partían de una posición distinta a la del equilibrio pero no tenían velocidad inicial impuesta. En este apartado se estudia el movimiento que describen las masas en 2 casos

• Velocidades iniciales de 1 m/s para ambas masas en el mismo sentido.
• Velocidades iniciales de 1 m/s para ambas masas en sentidos opuestos.

Estableciendo las condiciones iniciales: [math] Caso~1=\begin{cases} \dot x=1~[m/s]\\ \dot y=1~[m/s] \end{cases} [/math]:[math] Caso~2=\begin{cases} \dot x=1~[m/s] \\ \dot y=-1~[m/s] \end{cases} [/math] Las masas describen estos movimientos:
Caso 1

centro

Caso 2

centro