Oscilaciones forzadas en un sistema formado por partículas unidas por muelles (II)

Fuerza periódica

La solución de estas ecuaciones diferenciales es la suma de la solución de la ecuación homogénea que depende de las condiciones iniciales y la solución particular. El primer término describe el estado transitorio que desaparece al cabo de cierto tiempo (teóricamente infinito) si hay rozamiento y el segundo, el estado estacionario. En esta página, supondremos que las condiciones iniciales son: en el instante t=0, x0i=0 (las partículas parten del origen) con velocidad inicial nula v0i=0.

Primero, consideraremos que no hay rozamiento que es el caso más sencillo y posteriormente, le añadiremos el rozamiento.

Sin rozamiento

En este apartado, estudiaremos el comprortamiento de un sistema bajo la acción de fuerzas, sin tener en cuenta el rozamiento.

En la figura, representamos un sistema de dos partículas de masas m1 y m2 unidas por tres muelles de constantes k1, k2 y k3. Sobre la primera partícula actúa una fuerza dependiente del tiempo F1(t) y sobre la segunda la fuerza F2(t). Las ecuaciones del movimiento son

m 1 d 2 x 1 d t 2 = k 1 x 1 + k 2 ( x 2 x 1 )+ F 1 (t) m 2 d 2 x 2 d t 2 = k 3 x 2 k 2 ( x 2 x 1 )+ F 2 (t)

En forma matricial escribimos

( m 1 0 0 m 2 )( d 2 x 1 d t 2 d 2 x 2 d t 2 )+( k 1 + k 2 k 2 k 2 k 2 + k 3 )( x 1 x 2 )=( F 1 (t) F 2 (t) ) M x ¨ +Kx=F

En la página anterior, calculamos los valores propios y los vectores propios de la matriz M-1·K, donde M es la matriz diagonal de las masas y K de las constantes elásticas de los muelles.

>> [V,D]=eig(inv(M)*K)

La matriz diagonal D contiene en su diagonal principal los cuadrados de las frecuencias de los modos normales de vibración, ω1 y ω2.

Los vectores columna de la matriz V son los vectores propios correspondientes a cada uno de los valores propios. Modificamos la matriz V multiplicando cada vector propio por un factor de escala, de modo que se cumpla

M g = ( V ) T ·M·V=( 1 0 0 1 ) K g = ( V ) T ·K·V=( ω 1 2 0 0 ω 2 2 )

Definimos un nuevo vector u(t) de modo que

x(t)=V·u(t) x ¨ (t)=V· u ¨ (t)

La ecuación diferencial del movimiento en forma matricial se transforma en

M·V· u ¨ (t)+KV·u(t)=F ( V ) T ·M·V· u ¨ (t)+ ( V ) T ·KV·u(t)= ( V ) T ·F M g · u ¨ (t)+ K g ·u(t)=f(t) M g =( 1 0 0 1 ) K g =( ω 1 2 0 0 ω 2 2 )

El comportamiento del sistema se describe mediante un sistema de dos ecuaciones diferenciales desacopladas

{ d 2 u 1 d t 2 + ω 1 2 u 1 = f 1 (t) d 2 u 2 d t 2 + ω 2 2 u 2 = f 2 (t)

La solución de cada una de estas ecuaciones diferenciales es la suma de su homogénea que depende de las condiciones iniciales y de la solución particular que depende de la expresión de cada una de las fuerzas f1(t) y f2(t).

Una vez que obtenemos las soluciones en términos de las coordendas u1(t) y u2(t) se vuelven a transformar en coordenadas físicas x1(t) y x2(t), mediante x=V·u.

Ejemplo

Supongamos un sistema de dos partículas de masas m1=2 y m2=1, unidas por muelles de constantes k1=6, k2=3 y k3=0. Sobre la segunda partícula actúa una fuerza F2(t)=2·cos(ωf t), sobre la primera no actúa ninguna fuerza F1(t)=0.

Las ecuaciones del movimiento en forma matricial se escriben

( 2 0 0 1 )( d 2 x 1 d t 2 d 2 x 2 d t 2 )+( 9 3 3 3 )( x 1 x 2 )=( 0 2cos( ω f t) ) M=( 2 0 0 1 )K=( 9 3 3 3 )

Calculamos los valores propios (cuadrados de las frecuencias de los modos normales de vibración) y los vectores propios correspondientes a cada valor propio de la matriz M-1·K

>> M=sym('[2,0;0,1]');
>> K=sym('[9,-3;-3,3]');
>> [V,D]=eig(inv(M)*K)
V = 
[ 1/2, -1]
[   1,  1]
D =
[ 3/2, 0]
[   0, 6]

Los cuadrados de las frecuencias de los modos normales de vibración son los elementos de la diagonal de la matriz D.

ω 1 2 = 3 2 ω 2 2 =6

Creamos una nueva matriz V, multiplicando los vectores propios por un factor de escala de modo que

M g = ( V ) T ·M·V=( 1 0 0 1 ) K g = ( V ) T ·K·V=( ω 1 2 0 0 ω 2 2 )

>> X1=V(:,1);
>> X2=V(:,2);
>> r=X1'*M*X1
r =3/2
>> X1=X1/sqrt(r);
>> r=X2'*M*X2
r =3
>> X2=X2/sqrt(r);
>> V=[X1,X2]
V =
[ 6^(1/2)/6, -3^(1/2)/3]
[ 6^(1/2)/3,  3^(1/2)/3]
>> Mg=V'*M*V
Mg =
[ 1, 0]
[ 0, 1]
>> Kg=V'*K*V
Kg =
[ 3/2, 0]
[   0, 6] 

La nueva matriz V es

V=( 6 6 3 3 6 3 3 3 )

En el sistema de coordenadas u, que hemos definido mediante la transformación x=V·u, las fuerzas f1(t) y f2(t) se expresan

>> syms t wf;
>> F=[0;2*cos(wf*t)];
>> f=V'*F
f =
 (2*2^(1/2)*3^(1/2)*cos(t*wf))/3
         (2*3^(1/2)*cos(t*wf))/3

f(t)= ( V ) T ·F=( 6 3 F 2 (t) 3 3 F 2 (t) )=( 2 6 3 cos( ω f t) 2 3 3 cos( ω f t) )

Las ecuaciones diferenciales del movimiento se desacoplan y su solución particular se puede calcular fácilmente

M g u ¨ + K g u=f(t) M g =( 1 0 0 1 ) K g =( ω 1 2 0 0 ω 2 2 )=( 3 2 0 0 6 ) { d 2 u 1 d t 2 + ω 1 2 u 1 = 2 6 3 cos( ω f t) d 2 u 2 d t 2 + ω 2 2 u 2 = 2 3 3 cos( ω f t)

La solución particular es u=Acos(ωf·t)+Bsin(ωf·t), pero B=0, como podemos comprobar en oscilaciones forzadas

u 1 (t)= 2 6 3 1 ω 1 2 ω f 2 cos( ω f t) u 2 (t)= 2 3 3 1 ω 2 2 ω f 2 cos( ω f t)

Calculamos las posiciones x1(t) y x2(t) de cada una de las partículas mediante la transformación x=V·u

( x 1 (t) x 2 (t) )=( 6 6 3 3 6 3 3 3 )( u 1 (t) u 2 (t) ) x 1 (t)= 2 3 1 ω 1 2 ω f 2 cos( ω f t) 2 3 1 ω 2 2 ω f 2 cos( ω f t) x 2 (t)= 4 3 1 ω 1 2 ω f 2 cos( ω f t)+ 2 3 1 ω 2 2 ω f 2 cos( ω f t) ω 1 = 3 2 ω 2 = 6

w1=sym('sqrt(3/2)');
w2=sym(sqrt(6)');
wf=sym('1');  %frecuencia de la fuerza oscilante
syms t;

x1=(2/3)*cos(wf*t)/(w1^2-wf^2)-(2/3)*cos(wf*t)/(w2^2-wf^2);
x2=(4/3)*cos(wf*t)/(w1^2-wf^2)+(2/3)*cos(wf*t)/(w2^2-wf^2);

hold on
h=ezplot(x1,[0,10]);
set(h,'color','b');
h=ezplot(x2,[0,10]);
set(h,'color','r');
hold off
title('Dos osciladores acoplados, forzados')
xlabel('t')
ylabel('x1, x2')
grid on

Más abajo representaremos x1(t) y x2(t) para varias frecuencias ωf.

Amplitud de las oscilaciones en función de la frecuencia ωf

Las dos partículas describen Movimientos Armónicos Simples x=Acos(ωf·t) cuyas amplitudes son

A 1 = 2 3 1 ω 1 2 ω f 2 2 3 1 ω 2 2 ω f 2 = 6 ( 32 ω f 2 )( 6 ω f 2 ) A 2 = 4 3 1 ω 1 2 ω f 2 + 2 3 1 ω 2 2 ω f 2 = 2( 92 ω f 2 ) ( 32 ω f 2 )( 6 ω f 2 )

Representamos las amplitudes A1y A2 en función de la frecuencia de la fuerza oscilante ωf.

wf=linspace(0,4,400);
A1=6./((3-2*wf.^2).*(6-wf.^2));
A2=2*(9-2*wf.^2)./((3-2*wf.^2).*(6-wf.^2));

A1=A1.*(abs(A1)<10); %limita el máximo y mínimo de A1
A2=A2.*(abs(A2)<10);
hold on
plot(wf,A1,'b')
plot(wf,A2,'r')
grid on
title('Amplitudes de los osciladores')
ylabel('A1,A2')
xlabel('\omega_f')
hold off

Supongamos un sistema bajo la acción de fuerzas del tipo Fi(t)=F0i·cos(ωf·t).

En el espacio u, en el estado estacionario, la soluciones particulares de las ecuaciones diferenciales desacopladas tienen la forma

ui=Auicos(ωf t), i=1,2

Obtendremos los valores de Au haciendo que cumpla cada una de las ecuaciones diferenciales

Mediante la transformación x=V·u, obtenemos las posiciones x1(t) y x2(t) de cada una de las dos partículas en función del tiempo

xi=Axicos(ωf t), i=1,2

Elaboraremos un script más general, para representar la respuesta en amplitud en el estado estacionario de las oscilaciones de las partículas de un sistema bajo la acción de fuerzas del tipo Fi(t)=F0i·cos(ωf t).

syms wf;
M=sym('[2,0;0,1]');
K=sym('[9,-3;-3,3]');

[V,D]=eig(inv(M)*K); %valores y vectores propios
w=diag(sqrt(D)) %vector de frecuencias propias
n=length(w);
for i=1:n
    r=V(:,i)'*M*V(:,i);
    V(:,i)=V(:,i)/sqrt(r);
end
%vector fuerza
F0=sym('[0;2]');
f0=V'*F0;

%Amplitudes
for i=1:n
    Au(i)=f0(i)/(w(i)^2-wf^2);
end
Ax=V*Au';

%representación gráfica
color=['b','r','g'];
hold on
for i=1:n
    h=ezplot(Ax(i),[0,4]);
    set(h,'color',color(i))
end
title('Respuesta en amplitud')
ylabel('A1,A2')
xlabel('\omega_f')
grid on
hold off

En esta gráfica vemos que las amplitudes A1(color azul) y A2 (color rojo) se hacen muy grandes para las frecuencias ω1=1.2 y ω2=2.4. La amplitud A2 se anula para la frecuencia ωf=2.1. Hay una frecuencia situada entre 1.5 y 2 para la cual A1=A2 . Igualamos ambas amplitudes (basta igualar los numeradores) y obtenemos la frecuencia buscada ω f = 3

Posición de las partículas en función del tiempo

Elaboramos un script que nos calcule las posiciones x1(t) y x2(t) de cada una de las partículas para distintos valores de la frecuencia de la fuerza oscilante ωf.

syms t wf;
M=sym('[2,0;0,1]');
K=sym('[9,-3;-3,3]');

[V,D]=eig(inv(M)*K)
w=diag(sqrt(D)); %vector de frecuencias propias
%calcula una nueva matriz V
for i=1:length(w)
    r=V(:,i)'*M*V(:,i);
    V(:,i)=V(:,i)/sqrt(r);
end
%vector fuerza
F=sym('[0;2*cos(wf*t)]');
f=V'*F;

%solución particular de las ecuaciones del movimiento desacopladas
for i=1:length(w)
    u(i)=f(i)/(w(i)^2-wf^2);
end
x=V*u';
x=subs(x,wf,sym('1'));  %cambiar la frecuencia de la fuerza oscilante

%representación gráfica
color=['b','r','g'];
hold on
for i=1:length(w)
    h=ezplot(x(i),[0,12]);
    set(h,'color',color(i))
end
title('Dos osciladores acoplados, forzados')
ylabel('x1,x2')
xlabel('t')
grid on
hold off

Respuesta a una fuerza constante:

....
    x=subs(x,wf,sym('0'));
....

Como apreciamos en las ecuaciones y en la gráfica x1(t)=1/3 y x2(t)=1

Respuesta a una fuerza, con ωf=1, ωf1 y ωf<ω2.

....
    x=subs(x,wf,sym('1'));
....

Para ωf=1, las amplitudes de estos dos MAS son, respectivamente: A1=6/5 y A2=14/5

Respuesta a una fuerza, con ω f = 3 2 , ω1<ωf<ω2.

F 1 (t)=0 F 2 (t)=2cos( 3 2 t )

....
    x=subs(x,wf,sym('3/sqrt(2)'));
....

Como vemos la segunda partícula no se mueve x2(t). A esta situación se denomina antiresonancia. Calculamos esta frecuencia haciendo que la amplitud A2=0

A 2 = 2( 92 ω f 2 ) ( 32 ω f 2 )( 6 ω f 2 ) A 2 =0 ω f = 3 2

Las amplitudes de estos dos MAS son, respectivamente: A1=2/3 y A2=0

Respuesta a una fuerza, con ωf=3, ωf1 y ωf>ω2.

....
    x=subs(x,wf,sym('3'));
....

Para ωf=3, las amplitudes de estos dos MAS son, respectivamente: A1=2/15 y A2=2/5. Las dos partículas oscilan en oposición de fase.