Diferencia entre revisiones de «Sistemas resorte-masa»

Revisión actual del 00:30 3 jun 2013

Es el trabajo propuesto en la asignatura de Ecuaciones Diferenciales y realizado por el grupo 18, que plantea un sistema de dos masas y tres muelles.

1 Sistemas resorte-masa

Consideremos un sistema formado por dos masas ancladas a la pared por los muelles de constantes k1 y k3, y unidos entre ellos por otro de constante k2, tal y como se indica en el dibujo.

Sistema formado por dos masas y tres muelles, dos de ellos anclados a las paredes laterales

El objetivo será deducir las ecuaciones del movimiento de cada masa. Para ello, establecemos dos parámetros “x” e “y” que nos indican los desplazamientos de las masas respecto de su posición de equilibrio respectivamente. (NOTA: Tomaremos el desplazamiento positivo hacia la derecha). Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que:

[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]

2 Aproximación de la posición mediante métodos numéricos

Para aplicar tanto Runge-Kutta como Newmark necesitamos un sistema de ecuaciones donde sólo haya derivadas primeras, para ello hacemos los siguientes cambios de variable:

[math]\dot x=u\\\dot y=v[/math]

Y así obtenemos el sistema: [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)\\\dot x=u\\\dot y=v\end{matrix}\right.[/math]

En estas ecuaciones sustituimos los datos que plantean en el problema:

[math]m_1= 2kg\\m_2= 1kg\\k_1= 4N/m\\k_2= 2N/m\\k_3= 1N/m[/math]

Y para poder resolver el problema de valor inicial, se toman los valores también dados en el enunciado que se muestran a continuación: [math]x(0)= 1\\y(0)= 1,5\\u(0)= 0\\v(0)= 0[/math]

2.1 Método Runge-Kutta

Para poder aplicar el método necesitamos pasar el sistema de ecuaciones a uno matricial, obteniendo así M.

Sabiendo que el método de Runge-Kutta consiste en: [math]y_0\\y_{n+1} =y_n+h/6(k_1+2k_2+2k_3+k_4)\\k_1=f(t_n,y_n)\\k_2=f(t_n+1/2h,y_n+1/2k_1h)\\k_3=f(t_n+1/2h,y_n+1/2k_2h)\\k_4=f(t_n+h,y_n+k_3h)[/math] En nuestro, el programa queda de la siguiente manera:

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;
 
t0=0; tN=10; N=400; 
h=(tN-t0)/N;
t=t0:h:tN;
y(:,1)=[1;0;1.5;0]; %Condiciones iniciales x0=1; u0=0; y0=1.5; w0=0
 
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; %Matriz del sistema
 
%Aplicación del método Runge-Kutta
for n=1:N
    k1=M*y(:,n);
    k2=M*(y(:,n)+1/2*k1*h);
    k3=M*(y(:,n)+1/2*k2*h);
    k4=M*(y(:,n)+k3*h);
    y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
 
figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off
 
%Para visualizar mejor el comportamiento de ambas masas, planteamos la
%siguiente gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black');
plot(y(2,:),t,'green');
plot(y(3,:),t,'black');
plot(y(4,:),t,'green');
hold off

Así obtendremos la siguiente gráfica de la posición y la velocidad de ambas masas a lo largo del tiempo:

Las posiciones se muestran en negro y las velocidades en verde

2.2 Método de Newmark

Una vez introducidas las dos variables “u” y “v” que serán igual a las velocidades (derivada primera de x e y). A continuación aplicamos las fórmulas del método de Newmark dadas para cada variable, obteniendo 4 ecuaciones en las que solo aparece la primera derivada.:

[math]x_{n+1} = x_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\y_{n+1} = y_n + hz_n + h^2(βf_n+1 + (1/2 - β)f_n)\\u_{n+1} = u_n + h(γf_(n+1) + (1 - γ)f_n)\\v_{n+1} = v_n + h(γf_(n+1) + (1 - γ)f_n)[/math]

Ahora ordenamos matricialmente los términos obteniendo las matrices A y B correspondientes a los términos n+1 y n respectivamente. Estas matrices se ven reflejadas en el programa de Matlab que se muestra a continuación.

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
% Newmark parameters
beta=1/4;
gamma=1/2;
N=400; % number of steps
h=(tN-t0)/N; % mesh size
Y=[1;1.5;0;0]; %Condiciones iniciales x0=1; y0=1.5; u0=0; w0=0;
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]

B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
for n=1:N
Y(:,n+1)=inv(A)*B*Y(:,n);
end

% Plot the numerical approximation
figure(1) %desplazamientos
t=t0:h:tN;
hold on
plot(t,Y(1,:)); %x
plot(t,Y(2,:),'black'); %y
hold off

figure(2) %velocidades
hold on
plot(t,Y(3,:),'green'); %u=x'
plot(t,Y(4,:),'red');   %w=y'
hold off

Así pues obtendremos las gráficas de las posiciones y velocidades de las partículas

Masa 1 en color azul, masa 2 en negro

La velocidad de la masa 1 se muestra en rojo y la de la masa 2 en verde

2.3 Comparación

A pesar de la diferente disposición de las gráficas en el método de Newmark y en el de Runge-Kutta, es apreciable con la propia vista que el resultado obtenido por ambos métodos es muy parecido, aun así hemos planteado un código de Matlab de tal forma que se superpongan ambas gráficas. Al realizar esto se puede apreciar que hasta que no ampliamos la imagen hasta obtener escalas centesimales no se parecia la variación, por lo que ambos métodos dan resultados muy similares. Además, con el fin de asegurarnos de realizar una comparación precisa, se ha planteado en el código un cociente entre las posiciones obtenidas según ambos métodos, obteniéndose como resultado siempre la unidad.

%Resolucion del sistema x'=u; u=-3x+y; y'=w; w=2x+3y;
 
t0=0; tN=10; N=400; 
h=(tN-t0)/N;
t=t0:h:tN;
 
%RUNGE-KUTTA
y=[1;0;1.5;0];
 
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0]; 
 
for n=1:N
    k1=M*y(:,n);
    k2=M*(y(:,n)+1/2*k1*h);
    k3=M*(y(:,n)+1/2*k2*h);
    k4=M*(y(:,n)+k3*h);
    y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
 
%NEWMARK
beta=1/4;
gamma=1/2;
S=[1;1.5;0;0]; 
 
A=[1+3*(h^2)/4,-(h^2)/4,0,0;-(h^2)/2,1+3*(h^2)/4,0,0;3*h/2,-h/2,1,0;-h,3*h/2,0,1]
B=[1-3*(h^2)/4,(h^2)/4,h,0;(h^2)/2,1-3*(h^2)/4,0,h;(-3)*h/2,h/2,1,0;h,-3*h/2,0,1]
 
for n=1:N
    S(:,n+1)= A\(B*S(:,n));
end
 
%COMPARACION DE DESPLAZAMIENTOS
figure(1)
hold on
y(3,:)=y(3,:)+4
S(2,:)=S(2,:)+4
plot(y(1,:),t,'black');
plot(y(3,:),t,'black');
plot(S(1,:),t,'red');
plot(S(2,:),t,'red');
hold off
 
figure(2)
z=y(1,:)/S(1,:);
plot(t,z);
%como podemos ver el la gráfica que nos divide las x de rk con las de
%newmark, el cociente es 1 o extremadamente próximo a 1, lo que quiere
%decir que la diferencia entre ambos métodos en prácticamente nula.
 
%NOTA:se podría hacer lo mismo comparando las velocidades y saldría
%lo mismo.

Al activar el código aquí se planteado se obtienen las siguientes gráficas:

Posiciones de ambas masas obtenidas por los dos métodos superpuestas

Al hacer zoom en la gráfica se acaba apreciando una separación entre ellas

Cociente entre las posiciones de la masa 1 según ambos métodos

3 Posiciones de las masas en función del tiempo para los siguientes casos

3.1 Partiendo del equilibrio con Vi=1m/s en el mismo sentido para ambas masas

Lo único que se alterará en el problema de valor inicial serán los valores iniciales dados para u y v, derivadas primeras de x e y.:

[math]u(0)= 1m/s\\v(0)= 1m/s[/math]

El código de Matlab aplicado es el mismo que en el apartado dos (ambos son válidos), con la unica diferencia de la variación de los valores iniciales, por lo que no se mostrará de nuevo. A continuación se mostrará el utilizado siguiendo el método Runge-Kutta:

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;1;1.5;1]; %Condiciones iniciales x0=1; u0=1; y0=1.5; w0=1;
M=[0,1,0,0;-3,0,1,0; 0,0,0,1; 2,0,-3,0];
for n=1:N
    k1=M*y(:,n);
    k2=M*(y(:,n)+1/2*k1*h);
    k3=M*(y(:,n)+1/2*k2*h);
    k4=M*(y(:,n)+k3*h);
    y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end
figure(1)
hold on
plot(t,y(1,:)); %x en azul
plot(t,y(2,:),'green'); %u en verde
plot(t,y(3,:),'black'); %y en negro
plot(t,y(4,:),'red'); %w en rojo
hold off
%Gráfica alternativa.
figure(2)
hold on
y(3,:)=y(3,:)+4
y(4,:)=y(4,:)+4
plot(y(1,:),t,'black'); %posición de m1
plot(y(2,:),t,'green'); %velocidad de m1
plot(y(3,:),t,'black'); %posición de m2
plot(y(4,:),t,'green'); %velocidad de m2
hold off

Las gráfica que define el movimiento en este caso es la siguiente:

Las posiciones se muestran en negro y las velocidades en verde

3.2 Partiendo del equilibrio con Vi=1m/s en sentido opuesto para ambas masas

Este segundo caso es muy similar al anterior, con la diferencia de que la velocidad de la segunda masa es de sentido contrario, así pues los valores iniciales de las velocidades serán las siguientes:

[math]u(0)= 1m/s\\v(0)= -1m/s[/math]

La gráfica que define la posición y el movimiento de ambas partículas en este caso es la que se muestra a continuación:

Las posiciones se muestran en negro y las velocidades en verde

4 Sistema de una masa en medio viscoso

Suponemos ahora que el sistema tiene sólo una masa. Es decir, desenganchamos el muelle de la segunda masa. Suponemos también que el sistema está sumergido en un medio viscoso que provoca un amortiguamiento en el comportamiento del muelle proporcional a la velocidad de la masa, con coeficiente μ= 1. En este caso m=2 k=4

Sistema con amortiguamiento

Aplicando las ecuaciones de la mecánica clásica (Newton- Euler), tenemos que: [math]-C\dot x-Kx=m\ddot x\\\ddot x=-2x-(C/2)\dot x[/math]

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable: [math]\dot x=u[/math] Con lo que el sistema nos queda: [math]\dot x=u\\\dot u=-2x-(C/2)u[/math]

A continuación, aplicando los métodos de Runge-Kutta y Newmark, obtenemos:

4.1 Método Runge-Kutta

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
M=[0,1;-1,-1/2];
for n=1:N
    k1=M*y(:,n);
    k2=M*(y(:,n)+1/2*k1*h);
    k3=M*(y(:,n)+1/2*k2*h);
    k4=M*(y(:,n)+k3*h);
    y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end

figure(4)
hold on
plot(t,y(1,:),'red'); %posición de la partícula
plot(t,y(2,:),'green'); %velocidad de la partícula
hold off
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v)); %energía mecánica de la partícula
figure(5)
plot(t,E)

4.2 Método de Newmark

% Solve second order ODE with Newmark method
clear all
t0=0;
tN=10;
y0=1;
z0=0;

% Newmark parameters
beta=1/4;
gamma=1/2;
N=100; % number of steps
h=(tN-t0)/N; % mesh size
y(1)=y0;
z(1)=z0;
Y=[y(1);z(1)];
for n=1:N
    A=[1+(h^2)/4, (h^2)/8; h/2, 1+h/4];
    B=[1-(h^2)/4, h-(h^2)/8; ...
        -h/2, 1-h/4];
    Y=A\(B*Y);
    y(n+1)=Y(1);
    z(n+1)=Y(2);
end

% Plot the numerical approximation
figure(1)
t=t0:h:tN;
E=1/2*4*(y.*y)+1/2*2*(z.*z);
hold on
plot(t,y,'green');
plot(t,z)
hold off
figure(2)
plot(t,E)

A partir de ambos métodos se obtienen semejantes gráficas de la velocidad y posición que es la siguiente:

Se muestra la posición de la partícula en rojo y la velocidad en verde a lo largo del tiempo

4.3 Energía

Se muestra una gráfica de la energía de la partícula a lo largo del tiempo en escala logarítmica.

Energía mecánica de la partícula a lo largo del tiempo

4.4 Influencia de μ

En este último apartado del trabajo comprobaremos el efecto que tiene la constante del amortiguamiento sobre la posción, la velocidad y la energía. Para ello creamos un nuevo programa en el que creamos un bucle de forma que cambie la matriz en función de la cte. de amortguamiento y un nuevo bucle dentro de este para aproximar la ecuación diferencial por Runge-Kutta. Dentro del primer bucle y fuera y del segundo ploteamos las gráficas de la posición, la velocidad y la energía, obteniendo las tres para cada valor que da el bucle. Finalmente, comparamos las gráficas llegando a la conclusión de que al aumentar el amortiguamiento se reduce las oscilaciones, la masa frena con mayor rapidez y la energía se disipa antes al disminuir más la velocidad en el mismo tiempo obteniendo una pérdida energía cinética, puesto que es proporcional a la velocidad.

t0=0; tN=10; N=400;
h=(tN-t0)/N;
t=t0:h:tN;
y=[1;0];
%lo que vamos a hacer es crear un bucle que me va a coger 10 valores de mu,
%desde 1 hasta 5, aumentando de 0.5 en 0.5
g=5;
for i=1:0.5:g
M=[0,1;-2,-i/2];
for n=1:N
    k1=M*y(:,n);
    k2=M*(y(:,n)+1/2*k1*h);
    k3=M*(y(:,n)+1/2*k2*h);
    k4=M*(y(:,n)+k3*h);
   y(:,n+1)=y(:,n)+h/6*(k1+2*k2+2*k3+k4);
end 
figure(1)
hold on
subplot(g*2,3,1+(((i*2)-1)*3)) %el desplazamiento lo colocamos en la primera columna
plot(t,y(1,:),'red');
subplot(g*2,3,2+(((i*2)-1)*3))%la velocidad en la segunda columna
plot(t,y(2,:),'green');
x=y(1,:);
v=y(2,:);
E=(1/2*4*(x.*x))+(1/2*2*(v.*v));
subplot(g*2,3,3+(((i*2)-1)*3))%la energia en la tercera
plot(t,E)
end

Se muestran a continuación los cambios que produce la μ en la posición, la velocidad y la energía:

Desplazamiento, velocidad y energía aumentando μ de 1 a 5 en intervalos de 0,5