Diferencia entre revisiones de «Sistemas resorte-masa»

Revisión del 02:12 5 mar 2013

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.

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[/math]: [math]\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

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