Diferencia entre revisiones de «Sistema de muelles y masas (Grupo 3B)»

Revisión del 12:25 4 mar 2014

Muchos modelos de vibraciones se modelizan con sistemas de muelles y masas.

                                      [[Archivo: ]]

1 MODELIZACIÓN DEL SISTEMA DE ECUACIONES

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 k_1, k_2, k_3 y k_4 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 Runge-Kutta de cuarto orden.
  
• Método de Euler modificiado de segundo orden.
  

2.1 Método 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: [math] x_1’= z_1[/math]: [math] z_1’= x_1’’[/math]: [math] x_2’= z_2[/math]: [math] z_2’= x_2’’[/math]: [math] x_3’= z_3[/math]: [math] z_3’= x_3’’[/math]:

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: [math] x_1’= z_1[/math]: [math] z_1’= x_1’’=-3x_1+x_2[/math]: [math] x_2’= z_2[/math]: [math] z_2’= x_2’’=-3x_2+2x_1+x_3[/math]: [math] x_3’= z_3[/math]: [math] z_3’= x_3’’=-(4/3)x_3+(1/3)x_2[/math]:

El código en MATLAB del método sería:

% Las 6 filas que uso corresponden a x_1, x_2, x_3, x_1', x_2' y x_3' en ese orden
clear all
h=0.025;
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
t0=0;
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); % Cambiar el 1 por un 2 o un 3 si se quiere ver el resultado de los otros 2 apartados,
%             aunque el primero es el unico que se hace con RK4
 
y(:,1)=yy;
 
for i=1:N
% Calculamos k1=f(tn,yn)
k1=A*yy;
% Calculamos k2=f(tn+h/2,yn+1/2*h*k1)
k2=A*(yy+1/2*h*k1);
% Calculamos k3=f(tn+h/2,yn+1/2*h*k2)
k3=A*(yy+1/2*h*k2);
% Calculamos k4=f(tn+h,yn+h*k3)
k4=A*(yy+h*k3);
% Meto los datos en el vector
yy=yy+h/6*(k1+2*k2+2*k3+k4);
y(:,i+1)=yy;
end
 
y=y(1:3,:); % Para coger solo los datos de la posicion
y(1,:)=y(1,:)+ones(1,N+1)*2.5;
y(2,:)=y(2,:)+ones(1,N+1)*4;
y(3,:)=y(3,:)+ones(1,N+1)*8;
 
 
plot(y,t,':')
xlabel('Posicion de la masa')
ylabel('Tiempo')
legend('x1(t)','x2(t)','x3(t)')

Se obtiene la siguiente gráfica posición-tiempo:

Aproximación numérica del sistema de ecuaciones por Runge-Kutta

2.2 Método Euler Modificado

Al igual que con el método de Runge-Kutta, es necesario reducir el orden de la ecuación por lo que realizamos el mismo cambio de variable y por tanto se obtiene el mismo sistema de ecuaciones.

El código en MATLAB sería:

% Resolver Sistema de ecuaciones con Euler Modificado

% Las 6 filas que usamos corresponden a x1, x2, x3, x1', x2' y x3' en ese orden

clear all
h=0.1;
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
t0=0;
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); % Cambiar el 1 por un 2 o un 3 si se quiere ver el resultado de los otros 2 apartados,
%             aunque el primero es el unico que se hace con RK4

y(:,1)=yy;


for i=1:N
% Calculamos k1=h*f(tn,yn)
k1=A*yy;
% Calculamos k2=h*f(tn+h/2,yn+1/2*k1)
k2=A*(yy+h*k1);
% Meto los datos en el vector
yy=yy+(h/2)*(k1+k2);
y(:,i+1)=yy;
end

y=y(1:3,:); % Para coger solo los datos de la posicion

y(1,:)=y(1,:)+ones(1,N+1)*2.5;
y(2,:)=y(2,:)+ones(1,N+1)*4;
y(3,:)=y(3,:)+ones(1,N+1)*8;

plot(y,t)
xlabel('Posicion de las masas')
ylabel('Tiempo')
legend('x1(t)','x2(t)','x3(t)')

Se obtiene la siguiente gráfica posición-tiempo:

Aproximación numérica del sistema de ecuaciones por Euler Modificado

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 derecha, imponiendo a las tres masas una velocidad de 1 m/s hacia la derecha..

Se obtienen los desplazamientos utilizando el método del trapecio.

% Resolver Sistema de ecuaciones con el metodo trapezoidal
 
% Las 6 filas que usamos corresponden a x1, x2, x3, x1', x2' y x3' en ese orden
 
clear all
h=0.2;
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 los apartados 2, 3a y 3b
t0=0;
tf=20;
N=(tf-t0)/h;
t=t0:h:tf;
y=zeros(6,N+1);
 
yy=y0(:,2); % Cambiar el 1 por un 2 o un 3 si se quiere ver el resultado de los otros 2 apartados
 
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);
for i=1:N
y(:,i+1)=C*(y(:,i)+h/2*(A*y(:,i)));
end
 
y=y(1:3,:); % Para coger solo los datos de la posicion
 
plot(y,t,':')
xlabel('Posicion de las masas respecto al equilibrio')
ylabel('Tiempo')
legend('x1(t)','x2(t)','x3(t)')

Se obtiene la siguiente gráfica posición-tiempo:

Aproximación numérica del sistema de ecuaciones por el método del trapecio

Las posiciones relativas de las masas se obtienen utilizando el método de Euler modificado.

% Resolver Sistema de ecuaciones con Euler Modificado
 
% Las 6 filas que usamos corresponden a x1, x2, x3, x1', x2' y x3' en ese orden
 
clear all
h=0.05;
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
t0=0;
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); % Cambiar el 1 por un 2 o un 3 si se quiere ver el resultado de los otros 2 apartados,
%             aunque el primero es el unico que se hace con RK4
 
y(:,1)=yy;
 
 
for i=1:N
% Calculamos k1=h*f(tn,yn)
k1=h*A*yy;
% Calculamos k2=h*f(tn+h/2,yn+1/2*k1)
k2=h*A*(yy+1/2*k1);
% Meto los datos en el vector
yy=yy+k2;
y(:,i+1)=yy;
end
 
y=y(1:3,:); % Para coger solo los datos de la posicion
 
y(1,:)=y(1,:)+ones(1,N+1)*2.5;
y(2,:)=y(2,:)+ones(1,N+1)*4;
y(3,:)=y(3,:)+ones(1,N+1)*8;
subplot(2,1,1)
plot(y,t)
xlabel('Posicion de las masas')
ylabel('Tiempo')
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)')

Se obtiene la siguiente gráfica posición-tiempo:

Aproximación numérica del sistema de ecuaciones por Euler Modificado

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

% Las condiciones iniciales se introducen en el programa con la siguiente matriz:
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]';

% En el supuesto 1 las condiciones iniciales se corresponden con la segunda columna de la matriz y0
yy=y0(:,2);

% Para el supuesto número 2 utilizamos la tercera columna de la matriz
yy=y0(:,3);

Aproximación numérica del sistema de ecuaciones por el método del trapecio

Aproximación numérica del sistema de ecuaciones por Euler Modificado

4 SISTEMA DE UNA MASA EN MEDIO VISCOSO CON FUERZA EXTERIOR APLICADA

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: :

[math]mx''=-kx-cx'+2exp(-0.01t)sen(2t)[/math]

Sustituyendo en la ecuación los datos del enunciado: : [math]m=1 kg[/math] : [math]k=2 N/m[/math] : [math]μ=c=1[/math]

Se obtiene: : [math]x''=-2x-x'+2exp(-0.01t)sen(2t)[/math]

La ecuación se resuelve usando el método trapezoidal en MATLAB con h = 0.02

clear all
clf
t0=0;
tf=60;
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);

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

plot(y(1,:),t)
xlabel('Posicion de la masa')
ylabel('Tiempo')

Aproximación numérica del sistema de ecuaciones por el método del trapecio

La energía asociada al sistema viene dada por la ecuación: : [math]E=Ec+Ep=(1/2)mx'^2+(1/2)kx^2[/math] Particularizando para los datos del enunciado: : [math]E=(1/2)x'^2+x^2[/math]

clear all
clf
t0=0;
tf=60;
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);

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

x=y(1,:);
v=y(2,:);
e=0.5*(v.*v)+(x.*x);
plot(t,e)
xlabel('Tiempo')
ylabel('Energía')

Aproximación numérica del sistema de ecuaciones por el método del trapecio