Modelo de vibraciones con muelles y masas. Grupo 8-C
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelo de vibraciones con muelles y masas. Grupo 8-C |
| Asignatura | Ecuaciones Diferenciales |
| Curso | 2015-16 |
| Autores | Enrique Pírez, Ángel Robisco, Adrián Vera, Javier Olalde. |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas. En nuestro caso utilizamos cuatro muelles y tres masas unidas a dos paredes que se deslizan libremente sobre un plano horizontal.
Contenido
1 Modelización
El sistema consta de 3 masas m1, m2 y m3 y 4 muelles de constantes de elasticidad k1, k2, k3 y k4. Para realizar el trabajo suponemos el desplazamiento positivo a la derecha. Como el peso está equilibrado con la normal las únicas fuerzas que influyen son las recuperadoras de los muelles, puesto que no existe rozamiento ni fuerza viscosa o amortiguadora. Las fuerzas que se aplican sobre cada particula:
Sobre la primera masa: F1=-k1x1 F2=k2(x2-x1)
Sobre la segunda masa: F3=-k2(x2-x1) F4=k3(x3-x2)
Sobre la tercera masa: F5=-k3(x3-x2) F6=-k4x3
Aplicando la segunda Ley de Newton, F=ma, se obtiene el siguiente sistema de ecuaciones diferenciales:
m1a1 = -k1x1+k2(x2-x1) m2a2 = -k2(x2-x1)-k3(x2-x3) m3a3 = -k3(x3-x2)-k4x3
2 Resolución del sistema de ecuaciones
En el instante t=0, las tres masas están desplazadas 0.5,1 y 0.8 metros hacia la derecha respecto de la posición de equilibrio y se sueltan sin velocidad.
Para su resolución utilizaremos:
- Método de Runge-Kutta (4º Orden).
- Método de Euler modificado (2º Orden).
2.1 Método de Runge-Kutta
Para poder aplicar este método, se necesita reducir el orden de la ecuación, descomponiendo el sistema en seis ecuaciones de primer orden. Para ello se realizan los siguientes cambios de variable: x1’= z1 -> z1’= x1’’ x2’= z2 -> z2’= x2’’ x3’= z3 -> z3’= x3’’
Sustituyendo los datos proporcionados en el enunciado: m1=2 kg, k1=4 N/m, m2=1kg, k2=2 N/m, m3=3kg, k3=1 N/m.
Se obtiene el sistema de ecuaciones:
x1'= z1 -> z1'= x1=-3x1+x2
x2'= z2 -> z2'= x2=-3x2+2x1+x3
x3'= z3 -> z3'= x3=-(4/3)x3+(1/3)x2
Con todo ello, podemos introducir los datos en Matlab.
%Método de Runge-Kutta
%El vector generado está creado en el orden x1, x2, x3, x1', x2' y x3'
%Tamaño de paso
h=0.025;
%Condiciones iniciales
y0=[0.5 1 0.8 0 0 0;-0.5 1 -0.5 1 1 1;0.5 1 -0.5 -1 0 0.5]'; %Matriz con las condiciones iniciales
%Tiempo inicial
t0=0;
%Tiempo final
tf=20;
N=(tf-t0)/h;
t=t0:h:tf;
y=zeros(6,N+1);
A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;-3 1 0 0 0 0;2 -3 1 0 0 0;0 1/3 -4/3 0 0 0];
yy=y0(:,1);
y(:,1)=yy;
%Bucle de Runge-Kutta
for i=1:N
%Cálculo de k1=f(tn,yn)
k1=A*yy;
%Cálculo de k2=f(tn+h/2,yn+1/2*h*k1)
k2=A*(yy+1/2*h*k1);
%Cálculo de k3=f(tn+h/2,yn+1/2*h*k2)
k3=A*(yy+1/2*h*k2);
%Cálculo de k4=f(tn+h,yn+h*k3)
k4=A*(yy+h*k3);
%Vector generado
yy=yy+h/6*(k1+2*k2+2*k3+k4);
y(:,i+1)=yy;
end
%Datos de posición
y=y(1:3,:);
%Primera masa
y(1,:)=y(1,:)+ones(1,N+1)*2.5;
%Segunda masa
y(2,:)=y(2,:)+ones(1,N+1)*4;
%Tercera masa
y(3,:)=y(3,:)+ones(1,N+1)*8;
%Dibujamos el gráfico
plot(y,t,'-')
xlabel('Posicion de la masa')
ylabel('Tiempo')
title('Runge-Kutta')
legend('x1(t)','x2(t)','x3(t)');
Obtenemos la gráfica de las tres posiciones de las masas. Se observa que para tiempo 0 en el eje x se dan las posiciones iniciales (separadas 0, 0,8 y 1 metros de la posición de equilibrio).
2.2 Método de Euler modificado
Realizamos el mismo cambio de variable e introducimos los datos en Matlab.
%Método de Euler Modificado
%El vector generado está creado en el orden x1, x2, x3, x1', x2' y x3'
%Tamaño de paso
h=0.1;
%Condiciones iniciales
y0=[0.5 1 0.8 0 0 0;-0.5 1 -0.5 1 1 1;0.5 1 -0.5 -1 0 0.5]'; % Matriz con las condiciones iniciales
%Tiempo inicial
t0=0;
%Tiempo final
tf=20;
N=(tf-t0)/h;
t=t0:h:tf;
y=zeros(6,N+1);
A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;-3 1 0 0 0 0;2 -3 1 0 0 0;0 1/3 -4/3 0 0 0];
yy=y0(:,1);
y(:,1)=yy;
%Bucle de Euler modificado
for i=1:N
%Cálculo de k1=h*f(tn,yn)
k1=A*yy;
%Cálculo de k2=h*f(tn+h/2,yn+1/2*k1)
k2=A*(yy+h*k1);
%Vector generado
yy=yy+(h/2)*(k1+k2);
y(:,i+1)=yy;
end
%Datos de posición
y=y(1:3,:);
%Primera masa
y(1,:)=y(1,:)+ones(1,N+1)*2.5;
%Segunda masa
y(2,:)=y(2,:)+ones(1,N+1)*4;
%Tercera masa
y(3,:)=y(3,:)+ones(1,N+1)*8;
%Dibjamos el gráfico
plot(y,t,'-')
xlabel('Posicion de las masas')
ylabel('Tiempo')
title('Euler Modificado')
legend('x1(t)','x2(t)','x3(t)');
2.3 Comparación de ambos métodos
3 Velocidad de las masas y trayectorias
3.1 Supuesto Nº1
La primera masa se desplaza a su izquierda 0.5 metros, la segunda 1 metro a su derecha, y la tercera 0.5 metros hacia la derecha, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha.
Se obtienen los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio.
%Método del trapecio
%El vector generado está creado en el orden x1, x2, x3, x1', x2' y x3'
%Tamaño de paso
h=0.2;
%Condiciones iniciales
y0=[0.5 1 0.8 0 0 0;-0.5 1 -0.5 1 1 1;0.5 1 -0.5 -1 0 0.5]'; % Matriz con las condiciones iniciales
%Tiempo inicial
t0=0;
%Tiempo final
tf=20;
N=(tf-t0)/h;
t=t0:h:tf;
y=zeros(6,N+1);
yy=y0(:,2);
y(:,1)=yy;
A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;-3 1 0 0 0 0;2 -3 1 0 0 0;0 1/3 -4/3 0 0 0];
C=inv(eye(size(A))-h/2*A);
%Bucle del trapecio
for i=1:N
y(:,i+1)=C*(y(:,i)+h/2*(A*y(:,i)));
end
%Datos de posición
y=y(1:3,:);
%Dibujamos el gráfico
plot(y,t,'-')
xlabel('Posicion de las masas desde equilibrio')
ylabel('Tiempo')
title('Trapecial')
legend('x1(t)','x2(t)','x3(t)');
Para obtener las posiciones relativas de las masas utilizamos el método de Euler modificado.
%Método de Euler Modificado
%El vector generado está creado en el orden x1, x2, x3, x1', x2' y x3'
%Tamaño de paso
h=0.05;
%Condiciones iniciales
y0=[0.5 1 0.8 0 0 0;-0.5 1 -0.5 1 1 1;0.5 1 -0.5 -1 0 0.5]'; % Matriz con las condiciones iniciales de cada apartado
%Tiempo inicial
t0=0;
%Tiempo final
tf=10;
N=(tf-t0)/h;
t=t0:h:tf;
y=zeros(6,N+1);
A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;-3 1 0 0 0 0;2 -3 1 0 0 0;0 1/3 -4/3 0 0 0];
yy=y0(:,2);
y(:,1)=yy;
%Bucle de Euler modificado
for i=1:N
%Cálculo de k1=h*f(tn,yn)
k1=h*A*yy;
%Cálculo de k2=h*f(tn+h/2,yn+1/2*k1)
k2=h*A*(yy+1/2*k1);
%Vector generado
yy=yy+k2;
y(:,i+1)=yy;
end
%Datos de posición
y=y(1:3,:);
%Primera masa
y(1,:)=y(1,:)+ones(1,N+1)*2.5;
%Segunda masa
y(2,:)=y(2,:)+ones(1,N+1)*4;
%Tercera masa
y(3,:)=y(3,:)+ones(1,N+1)*8;
%Dibujamos los gráficos
subplot(2,1,1)
plot(y,t)
xlabel('Posicion de las masas')
ylabel('Tiempo')
title('Euler Modificado')
legend('x1(t)','x2(t)','x3(t)');
y(2,:)=[];
subplot(2,1,2)
plot(y,t)
xlabel('Posicion de la masas 1 y 3')
ylabel('Tiempo')
legend('x1(t)','x3(t)')
3.2 Supuesto Nº2
La primera masa se desplaza a su derecha 0.5 metros con una velocidad de 1 m/s en dirección hacia la izquierda, la segunda masa se suelta sin velocidad inicial con un desplazamiento de 1 metro a la derecha, y la tercera masa se desplaza a su izquierda 0.5 metros y se suelta con una velocidad de 0.5 metros hacia la derecha.
Para obtener los desplazamientos respecto a la posición de equilibrio utilizando el método del trapecio, y para hallar las posiciones relativas de las masas mediante el uso del método de Euler Modificado, el código de MATLAB para ambos casos sería el mismo que en el supuesto 1, cambiando las condiciones iniciales.
%Método del trapecio
%El vector generado está creado en el orden x1, x2, x3, x1', x2' y x3'
%Tamaño de paso
h=0.2;
%Condiciones iniciales
y0=[0.5 1 0.8 0 0 0;-0.5 1 -0.5 1 1 1;0.5 1 -0.5 -1 0 0.5]'; % Matriz con las condiciones iniciales
%Tiempo inicial
t0=0;
%Tiempo final
tf=20;
N=(tf-t0)/h;
t=t0:h:tf;
y=zeros(6,N+1);
yy=y0(:,3);
y(:,1)=yy;
A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;-3 1 0 0 0 0;2 -3 1 0 0 0;0 1/3 -4/3 0 0 0];
C=inv(eye(size(A))-h/2*A);
%Bucle del trapecio
for i=1:N
y(:,i+1)=C*(y(:,i)+h/2*(A*y(:,i)));
end
%Datos de posición
y=y(1:3,:);
%Dibujamos el gráfico
plot(y,t,'-')
xlabel('Posicion de las masas desde equilibrio')
ylabel('Tiempo')
title('Trapecial')
legend('x1(t)','x2(t)','x3(t)');
%Método de Euler Modificado
%El vector generado está creado en el orden x1, x2, x3, x1', x2' y x3'
%Tamaño de paso
h=0.05;
%Condiciones iniciales
y0=[0.5 1 0.8 0 0 0;-0.5 1 -0.5 1 1 1;0.5 1 -0.5 -1 0 0.5]'; % Matriz con las condiciones iniciales de cada apartado
%Tiempo inicial
t0=0;
%Tiempo final
tf=10;
N=(tf-t0)/h;
t=t0:h:tf;
y=zeros(6,N+1);
A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;-3 1 0 0 0 0;2 -3 1 0 0 0;0 1/3 -4/3 0 0 0];
yy=y0(:,3);
y(:,1)=yy;
%Bucle de Euler modificado
for i=1:N
%Cálculo de k1=h*f(tn,yn)
k1=h*A*yy;
%Cálculo de k2=h*f(tn+h/2,yn+1/2*k1)
k2=h*A*(yy+1/2*k1);
%Vector generado
yy=yy+k2;
y(:,i+1)=yy;
end
%Datos de posición
y=y(1:3,:);
%Primera masa
y(1,:)=y(1,:)+ones(1,N+1)*2.5;
%Segunda masa
y(2,:)=y(2,:)+ones(1,N+1)*4;
%Tercera masa
y(3,:)=y(3,:)+ones(1,N+1)*8;
%Dibujamos los gráficos
subplot(2,1,1)
plot(y,t)
xlabel('Posicion de las masas')
ylabel('Tiempo')
title('Euler Modificado')
legend('x1(t)','x2(t)','x3(t)');
y(2,:)=[];
subplot(2,1,2)
plot(y,t)
xlabel('Posicion de la masas 1 y 3')
ylabel('Tiempo')
legend('x1(t)','x3(t)')
4 Sistema de una masa en un medio viscoso aplicando una fuerza exterior
Suponemos que el sistema está constituido por un único muelle y una sola masa, sobre la cual actúa una fuerza exterior en la dirección del movimiento de la masa provocada por las vibraciones de la pared. Dicho sistema está sumergido en un medio viscoso que provoca un amortiguamiento, reduciendo la amplitud de la oscilación. Se aplica la segunda Ley de Newton al sistema y se obtiene la siguiente ecuación diferencial: mx=-kx-cx'+2exp(-0.01t)*sen(2t)
Sustituyendo en la ecuación los datos del enunciado: m=1 kg k=2 N/m μ=c=1
Se obtiene: x=-2x-x'+2exp(-0.01t)*sen(2t)
La ecuación se resuelve usando el método trapezoidal en MATLAB con h = 0.02 La fuerza exterior del enunciado como una fuerza de valor absoluto 2e^(-0.01t)·sen(2t) que actúa siempre en la dirección del movimiento, pues son vibraciones que van siempre en la dirección favorable al movimiento.
%Método trapecial
%Tiempo inicial
t0=0;
%Tiempo final
tf=60;
%Tamaño de paso
h=0.2;
t=t0:h:tf;
N=(tf-t0)/h;
y=zeros(2,N+1);
%Condiciones iniciales
y0=[0.3,0.3]';
y(:,1)=y0;
A=[0 1;-2 -1];
C=inv(eye(size(A))-h/2*A);
%Bucle método trapecial
for n=1:N
sv=y(2,n)/abs(y(2,n)); % Signo de la velocidad:con ella sabemos la direccion de la fuerza externa
y(:,n+1)=C*(y(:,n)+h/2*(A*y(:,n)+[0 sv*abs(2*exp(-0.01*t(n))*sin(2*t(n)))]'+[0 sv*abs(2*exp(-0.01*t(n+1))*sin(2*t(n+1)))]'));
end
%Dibujo de la gráfica
plot(y(1,:),t)
xlabel('Posicion de la masa')
ylabel('Tiempo')
title('Trapecial')
Con ello, procedemos al cálculo de la energía del sistema. La energía asociada a dicho sistema viene dada por la ecuación: E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2
Particularizada para los datos del enunciado: E=(1/2)x'^2+x^2
%Método del trapecio
%Tiempo inicial
t0=0;
%Tiempo final
tf=60;
%Tamaño de paso
h=0.2;
t=t0:h:tf;
N=(tf-t0)/h;
y=zeros(2,N+1);
y0=[0.3,0.3]';
y(:,1)=y0;
A=[0 1;-2 -1];
C=inv(eye(size(A))-h/2*A);
%Bucle del trapecio
for n=1:N
sv=y(2,n)/abs(y(2,n)); % Signo de la velocidad (Para saber la direccion de la fuerza externa)
y(:,n+1)=C*(y(:,n)+h/2*(A*y(:,n)+[0 sv*abs(2*exp(-0.01*t(n))*sin(2*t(n)))]'+[0 sv*abs(2*exp(-0.01*t(n+1))*sin(2*t(n+1)))]'));
end
%Cálculo de la energía
x=y(1,:);
v=y(2,:);
e=0.5*(v.*v)+(x.*x);
%Dibujamos la gráfica
plot(t,e)
xlabel('Tiempo')
ylabel('Energía')
title('Energia del sistema a lo largo del tiempo')
Este suceso se debe a la ecuación diferencial de segundo grado que se ha obtenido del sistema aplicando la segunda Ley de Newton, de la que se obtiene una solución particular y una general de la homogénea. La solución resultante es armónica, y en uno de los términos tendremos un e elevado a un término negativo, que depende del tiempo y de la constante μ, haciendo que cuanto mayor sean estos valores menor sea la solución x(t). Al ser menor la x, la energía será también cada vez menor. Desde el punto de vista de la física, es correcto, ya que cuanto mayor sea la viscosidad del medio, más se opone al movimiento y más rápido se disipará la energía.
Aquí podemos ver cómo se reduce la energía si en vez de utilizar μ=1 lo hacemos con μ=3.
