Diferencia entre revisiones de «Sistema resorte-masa»
| Línea 285: | Línea 285: | ||
hold on | hold on | ||
title('Newmark') | title('Newmark') | ||
| − | axis([0,5,0, | + | axis([0,5,0,12]) |
plot(x+2,t,'-') | plot(x+2,t,'-') | ||
plot(2,t,'-g') | plot(2,t,'-g') | ||
Revisión actual del 11:15 26 feb 2014
Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.
Contenido
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 u=-(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
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*k1u;
k2y=vv(n)+1/2*h*k1v;
k2u=-(k2+k1)/m1*(xx(n)+1/2*h*k1x)+k2/m1*(yy(n)+1/2*k1y*h);
k2v=k2/m2*(xx(n)+1/2*h*k1x)-((k2+k3)/m2)*(yy(n)+1/2*k1y*h);
% Tercera K
k3x=uu(n)+1/2*h*k2u;
k3y=vv(n)+1/2*h*k2v;
k3u=-(k2+k1)/m1*(xx(n)+1/2*h*k2x)+k2/m1*(yy(n)+1/2*k2y*h);
k3v=k2/m2*(xx(n)+1/2*h*k2x)-((k2+k3)/m2)*(yy(n)+1/2*k2y*h);
% Cuarta K
k4x=uu(n)+h*k3u;
k4y=vv(n)+h*k3v;
k4u=-(k2+k1)/m1*(xx(n)+h*k3x)+k2/m1*(yy(n)+k3y*h);
k4v=k2/m2*(xx(n)+h*k3x)-((k2+k3)/m2)*(yy(n)+k3y*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
t=t0:h:tN;
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
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
Caso 2
4 Muelle y amortiguador
En algunos casos el sistema presenta, además de la masa y el resorte, un amortiguador viscoso. El efecto de este amortiguador va reduciendo la amplitud de la oscilación hasta que la masa se estabiliza nuevamente en la posición de equilibrio.
Aplicando de nuevo la segunda ley de Newton al sistema y tomando u(t)=0 obtenemos la ecuación diferencial:
- [math]m_{1}\ddot x+c \dot x +k_{1}x=0[/math]
Tomando los parámetros anteriores::
[math]m_{1}=2~[kg][/math]: [math]k_{1}=4~[N/m][/math]: [math]\mu=c=1[/math] La ecuación se resuelve en MATLAB aplicando los metodos de Newmark y de Runge-Kutta
4.1 Método de Newmark
Al ser la ecuación de orden dos es necesario reducir el orden en uno. Aplicando el cambio de variable [math]\dot x=y[/math] resulta un sistema de ecuaciones diferenciales de orden uno:: [math]\dot x=y[/math]:
[math]m_{1}\dot y=-cy-k_{1}x[/math]
Con las condiciones iniciales [math]x_{0}=1~[m][/math] el código en MATLAB del método es:
% sistema masa-muelle con amortiguador viscoso
clear all
%tiempo
t0=0;
tN=10;
%datos
k1=4;
m1=2;
c=1;
%condiciones iniciales
x0=1;
y0=0;
%parametros newmark
beta=1/4;
gamma=1/2;
N=1000;
h=(tN-t0)/N;
%primeros valores del vector solucion
x(1)=x0;
y(1)=y0;
%vector solucion
S=[x(1) ; y(1)];
for n=1:N
A=[1+beta*(k1/m1)*h^2 , beta*(c/m1)*h^2 ; ...
h*gamma*k1/m1 , 1+h*gamma*c/m1 ];
B=[1-(1/2-beta)*k1/m1*h^2 , h-(1/2-beta)*c/m1*h^2 ; ...
-h*(1-gamma)*k1/m1 , 1-h*(1-gamma)*c/m1 ];
S=A\(B*S);
x(n+1)=S(1);
y(n+1)=S(2);
end
t=t0:h:tN;
figure(1)
hold on
title('Newmark')
axis([0,5,0,12])
plot(x+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
4.2 Método de Runge-Kutta
Para aplicar el método a la ecuación diferencial vuelve a ser necesario reducir el orden hasta uno. Aplicando el mismo cambio de variable que para Newmark el codigo de MATLAB que aproxima la solución es:
% sistema masa-muelle con amortiguador viscoso
clear all
%tiempo
t0=0;
tN=10;
%datos
k1=4;
m1=2;
c=1;
%condiciones iniciales
x0=1;
y0=0;
N=1000;
h=(tN-t0)/N;
%primeros valores del vector solucion
xx(1)=x0;
yy(1)=y0;
for n=1:N
% Primera K
k1x=yy(n);
k1y=(-c*yy(n)-k1*xx(n))/m1;
% Segunda K
k2x=yy(n)+1/2*h*k1y;
k2y=(-c*(yy(n)+1/2*k1y*h)-k1*(xx(n)+1/2*h*k1x))/m1;
%Tercera k
k3x=yy(n)+1/2*h*k2y;
k3y=(-c*(yy(n)+1/2*k2y*h)-k1*(xx(n)+1/2*h*k2x))/m1;
%Cuarta K
k4x=yy(n)+h*k3y;
k4y=(-c*(yy(n)+k3y*h)-k1*(xx(n)+h*k3x))/m1;
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);
end
% Representacion
t=t0:h:tN;
axis([0,5,0,10])
hold on
title('Runge-Kutta')
plot(xx+2,t,'-')
plot(2,t,'-g')
legend('Masa m1','Posicion de equilibrio')
hold off
4.3 Energía del sistema
En el primer caso al no existir rozamiento la energía del sistema permanece constante. Pero la presencia de un amortiguador viscoso hace que el sistema se vaya frenando y, por lo tanto, perdiendo energía. La ecuación de la energía del sistema es::
[math]E=\frac{1}{2}m_{1}\dot x ^2+ \frac{1}{2}k_{1}x^2[/math]:
4.4 Influencia del coeficiente [math]\mu[/math]
A continuación se muestran una serie de casos con valores de [math]\mu[/math] comprendidos entre 0 y 5. Las gráficas muestran las trayectorias y el decrecimiento de la energía para cada uno de los casos:
Caso 1 [math]\mu=0[/math] Caso limite sin amortiguamiento en el que la energía es constante
Caso 2 [math]\mu=0.5[/math]
Caso 3 [math]\mu=1[/math]
Caso 4 [math]\mu=2[/math]
Caso 5 [math]\mu=3[/math]
Caso 6 [math]\mu=5[/math]
La conclusión que se extrae en vista de las gráficas es que la influencia del coeficiente de amortiguamiento [math]\mu[/math] es muy clara, cuanto mayor es este coeficiente más rapidamente decrece la energía pues la masa se estabiliza antes en torno a la posición de equilibrio.
