Sistemas resorte-masa

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

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]

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)

Reducimos el orden de las ecuaciones, haciendo el siguiente cambio de variable: (cv) Con lo que el sistema nos queda: (sistema) Aplicando los métodos, obtenemos: (métodos)

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