Modelo de vibraciones con muelles masas (Grupo A8)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Modelo de vibraciones con muelles masas. Grupo 8-A |
| Asignatura | Ecuaciones Diferenciales |
| Curso | Curso 2015-16 |
| Autores | Daniel del Potro Gabin 1423
Humberto del Castillo 1281 Luis de la Fuente Puig 939 Ignacio del Río 1693 |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Enunciado y modelización del sistema
Muchos modelos de vibraciones se modelizan con sistemas de muelles masas. Consideramos aquí un ejemplo simple de cuatro muelles con constantes kiN/m; i= 1;2;3;4, (fuerzas restauradoras), y tres masas mi; i= 1;2;3, unidos a dos paredes que se deslizan libremente sobre un plano horizontal.
Se modeliza el sistema de ecuaciones diferenciales para el desplazamiento de 3 masas desde la posición de equilibrio. Los grados de libertad que se utilizan son x1, x2 y x3, siendo estos la distancia de cada masa a su posición de equilibrio (Se toma el desplazamiento positivo hacia la derecha).
Se desprecia la fuerza de rozamiento, por lo que se supone que deslizan sobre una superficie horizontal y lisa, y por tanto las únicas fuerzas a tener en cuenta son las ejercidas por los muelles. Según la ley de Hooke, la fuerza producida por un muelle es proporcional a su elongación.
Por tanto las fuerzas aplicadas a la primera masa serán: :
[math] F_1=-k_1x_1 [/math] :
[math] F_2=k_2(x_2-x_1) [/math]
Sobre la segunda masa: :
[math] F_3=-k_2(x_2-x_1) [/math] :
[math] F_4=-k_3x_2 [/math]
Y las fuerzas que actúan en la tercera masa: :
[math] F_5=-k_3(x_3-x_2) [/math] :
[math] F_6=-k_4x_3 [/math]
Donde k1, k2, k3 y k4 son las constantes elásticas del movimiento de cada muelle.
Aplicando la segunda Ley de Newton, [math] F=mx'' [/math], se obtiene el siguiente sistema de ecuaciones diferenciales: [math] m_1x_1''=-k_1x_1+k_2(x_2-x_1)[/math]: [math] m_2x_2''=-k_2(x_2-x_1)-k_3(x_2-x_3)[/math]: [math] m_3x_3''=-k_3(x_3-x_2)-k_4x_3[/math]
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 repentinamente,sin velocidad.Se resuelve el sistema de ecuaciones diferenciales utilizando dos métodos numéricos:
- Método de Euler modificado.
- Método de Runge-Kutta.
Primero planteamos el problema en Matlab:
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=0.8;x11=0;x22=0;x33=0;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z=[z1;z2];
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=input('Introduce instante final:');
h=input('Introduce tamaño de paso:');%Probar para h=0.1 y h=0.025
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%f=M*z son las solucines del sistema
%Selecciono la segunda mitad del vector solución, ya que contiene las
%posiciones y la parte superior las derivadas
Y lo resolvemos con ambos métodos:
%Método Euler modificado
zEM=z;
for i=1:N
K1=M*zEM(:,i);
zEM(:,i+1)=zEM(:,i)+h*(M*(zEM(:,i)+K1*h/2));
end
zEM=zEM([4 5 6],:);%Selecciono las soluciones de posición
zEM1=zEM(1,:)+2.5;%Posición masa 1
zEM2=zEM(2,:)+4;%Posición masa 2
zEM3=zEM(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(zEM1,t,zEM2,t,zEM3,t);
title('Método Euler Modificado');
hold off
%Método Runge-Kutta
zRK=z;
for j=1:N
K1=M*zRK(:,j);
K2=M*(zRK(:,j)+(1/2)*K1*h);
K3=M*(zRK(:,j)+(1/2)*K2*h);
K4=M*(zRK(:,j)+K3*h);
zRK(:,j+1)=zRK(:,j)+h/6*(K1+2*K2+2*K3+K4);
end
zRK=zRK([4 5 6],:);%Selecciono las soluciones de posición
zRK1=zRK(1,:)+2.5;%Posición masa 1
zRK2=zRK(2,:)+4;%Posición masa 2
zRK3=zRK(3,:)+8;%Posición masa 3
ylabel('Tiempo');xlabel('Posición');
hold on
subplot(1,2,2)
plot(zRK1,t,zRK2,t,zRK3,t);
title('Método Runge-Kutta orden 4');
hold off
Se obtiene la siguiente gráfica posición-tiempo para Euler-Modificado y Runge-Kutta con h=0.1:
Y lo comparamos con la gráfica obtenida con Euler-Modificado y Runge-Kutta con h=0.025:
3 EFECTO DE LA VELOCIDAD INICIAL EN LA TRAYECTORIA
Del mismo modo que en el caso anterior, las tres masas parten de una posición diferente a la del equilibrio, pero en este caso se impone una velocidad inicial. Se analiza la trayectoria de las masas en dos supuestos distintos.
3.1 Supuesto 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 izquierda, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha..
Se obtienen los desplazamientos utilizando el método del trapecio.
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=-0.5;x2=1;x3=-0.5;x11=1;x22=1;x33=1;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=20;%Tiempo final
h=0.2;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método del Trapecio
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
z(:,i+1)=inv(I-h/2*M)*(z(:,i)+h/2*(M*z(:,i)));%Despejando
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
plot(z1,t,z2,t,z3,t);
title('Método del trapecio');
hold off
Se obtiene la siguiente gráfica posición-tiempo:
3.2 Supuesto 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
%Datos iniciales
m1=2;m2=1;m3=3;%Masas
k1=4;k2=2;k3=1;k4=3;%Constante de los muelles
x1=0.5;x2=1;x3=-0.5;x11=-1;x22=0;x33=0.5;%Datos en t0
z1=[x11,x22,x33]';%Vector de las derivas en t0
z2=[x1,x2,x3]';%Vector con las posiciones iniciales
%Matrices
k=[-(k1+k2)/m1,k2/m1,0;k2/m2,-(k2+k3)/m2,k3/m2;0,k3/m3,-(k3+k4)/m3];
z0=[z1;z2];%Condiciones iniciales
%Creo una matriz 6x6 con 4 matrices 3x3
M=[zeros(size(k)),k;eye(size(k)),zeros(size(k))];
%z'=M*z es el sistema a resolver
t0=0;
tN=10;%Tiempo final
h=0.05;%Tamaño de paso
N=round((tN-t0)/h);
t=linspace(t0,tN,N+1);
%Método de Euler Modificado
z=zeros(6,N+1);
zz=z0;
z(:,1)=zz;
I=eye(size(M));
for i=1:N
K1=M*z(:,i);
z(:,i+1)=z(:,i)+h*(M*(z(:,i)+K1*h/2));
end
z=z([4 5 6],:);%Selecciono las soluciones de posición
z1=z(1,:)-2.5;%Posición masa 1
z2=z(2,:)+4;%Posición masa 2
z3=z(3,:)+8;%Posición masa 3
hold on
subplot(1,2,1);
plot(z1,t,z3,t);
title('Masas 1 y 3');
subplot(1,2,2);
plot(z1,t,z2,t,z3,t);
title('Masas 1,2 y 3');
hold off
4 SISTEMA DE UNA MASA EN MEDIO VISCOSO CON FUERZA EXTERIOR APLICADA
%Datos iniciales
m=1;%Masa
k=4;%Constante del muelle
x1=0.3;x11=0.3;%Datos en t0
x0=[x1;x11]';%Condiciones iniciales
t0=0;%Instante inicial
tN=60;%Instante final
nu=1;%Constante de amortiguamiento
h=0.2;%Tamaño de paso
N=(tN-t0)/h;
%Matriz de coeficientes
M=[-nu/m,-k/m;1,0];
%vector tiempo y matriz posicion
t=linspace(t0,tN,N+1);
x=zeros(2,N+1);
x(:,1)=x0;
%Para conocer el valor de la fuerza en cada instante t
f=[];
for i=1:N+1
f=[f;2*exp(-0.01*t(i))*sin(2*t(i))];
end
%Método del Trapecio
I=eye(size(M));
for i=1:N
x(:,i+1)=inv(I-h/2*M)*(x(:,i)+h/2*(M*x(:,i)+f(i)+f(i+1)));
end
subplot(2,1,1);
plot(t,x(2,:));%Tiempo respecto a velocidad
ylabel('Tiempo');xlabel('Velocidad');
E=k/2*(x(1,:).^2)+m/2*(x(2,:).^2);%Ecuación de la energía
subplot(2,1,2);
semilogx(t,E);
title('Cantidad de energía');
La gráfica obtenida es la siguiente:
