Oscilaciones amortiguadas

Fuerza de rozamiento proporcional a la velocidad

La experiencia nos muestra que la amplitud de un cuerpo vibrante tal como un resorte o un péndulo, decrece gradualmente hasta que se detiene.

Para explicar el amortiguamiento, supondremos que además de la fuerza elástica F=-kx, actúa otra fuerza opuesta a la velocidad Fr=-λv, donde λ es una constante que depende del sistema físico particular. Todo cuerpo que se mueve en el seno de un fluido viscoso en régimen laminar experimenta una fuerza de rozamiento proporcional a la velocidad y de sentido contrario a ésta.

La ecuación del movimiento se escribe

ma=-kx-λv

Expresamos la ecuación del movimiento en forma de ecuación diferencial, teniendo en cuenta que la aceleración es la derivada segunda de la posición x, y la velocidad es la derivada primera de x.

d 2 x d t 2 +2γ dx dt + ω 0 2 x=0 ω 0 2 = k m 2γ= λ m

Las raíces de la ecuación característica son:

s 2 +2γs+ ω 0 2 =0{ s 1 =γ+ γ 2 ω 0 2 s 1 =γ γ 2 ω 0 2

Las soluciones de la ecuación diferencial son

Oscilaciones amortiguadas

La posición x y velocidad dx/dt en función del tiempo t son

x=exp( γt ){ Asin(ωt)+Bcos(ωt) } dx dt =γexp( γt ){ Asin(ωt)+Bcos(ωt) }+ωexp( γt ){ Acos(ωt)Bsin(ωt) }

Los coeficientes A y B se determinan a partir de las condiciones iniciales: posición inicial x0 y velocidad inicial v0 en el instante t=0

t=0{ x 0 =B v 0 =γB+ωA x=exp(γt)( v 0 +γ x 0 ω sin(ωt)+ x 0 cos(ωt) )

Tomamos γ=7 y ω0=100. Las condiciones iniciales son: x0=5, v0=0

x0=5; %condiciones iniciales
v0=0;
g=7; %constante de amortiguamiento
w0=100; %frecuencia natural
w=sqrt(w0^2-g^2);
x=@(t) (x0*cos(w*t)+(v0+g*x0)*sin(w*t)/w).*exp(-g*t);
fplot(x,[0,0.7])
title('Oscilación amortiguada')
xlabel('t')
ylabel('x');
grid on

Otra forma alternativa de expresarlo es

x= A 0 exp(γt)sin( ωt+φ ) A 0 = A 2 + B 2 = v 0 2 +2 v 0 γ x 0 + x 0 2 ω 0 2 ω 0 2 γ 2 tanφ= B A = x 0 ω 0 2 γ 2 v 0 +γ x 0

Sea una oscilación amortiguada de frecuencia angular propia ω0=100 rad/s y cuya constante de amortiguamiento γ=7.0 s-1. Sabiendo que la partícula parte de la posición x0=5 con velocidad inicial nula, v0=0, escribimos la ecuación de la oscilación amortiguada.

La frecuencia angular de la oscilación amortiguada ω es

ω = 100 2 7 2 = 99.75 rad/s

A0=5.01
tanφ=14.25

La ecuación de la oscilación amortiguada es

x=5.01·exp(-7t)·sin(99.75t+1.5)

Representamos la posición x y velocidad v en función del tiempo t, señalando los puntos donde la velocidad es nula, o el desplazamiento x es máximo o mínimo

x=Aexp( γt )sin( ωt+φ ) v= dx dt =Aexp( γt ){ γsin( ωt+φ )+ωcos( ωt+φ ) } v=0, t n = arctan( ω/γ )+nπφ ω

w0=100; %frecuencia propia
g=7; %constante amortiguamiento
x0=5; %condiciones iniciales
v0=0; 
A0=sqrt((v0^2+2*v0*g*x0+x0^2*w0^2)/(w0^2-g^2)); %amplitud
phi=atan(x0*sqrt(w0^2-g^2)/(v0+g*x0)); %fase inicial
w=sqrt(w0^2-g^2);

subplot(2,1,1)
x=@(t) A0*exp(-g*t).*sin(w*t+phi);
hold on
fplot(x,[0,0.5])
for n=0:15
    t0=(atan(w/g)+n*pi-phi)/w;
    plot(t0,x(t0),'o','markersize',3,'markeredgecolor','r','markerfacecolor','r')
end
hold off
grid on
xlabel('t')
ylabel('x')
title('Posición')

subplot(2,1,2)
v=@(t) A0*exp(-g*t).*(-g*sin(w*t+phi)+w*cos(w*t+phi));
hold on
fplot(v,[0,0.5])
for n=0:15
    t0=(atan(w/g)+n*pi-phi)/w;
    plot(t0,0,'o','markersize',3,'markeredgecolor','r','markerfacecolor','r')
end
hold off
grid on
xlabel('t')
ylabel('v')
title('Velocidad')

Determinamos los puntos de intersección de la amplitud A0exp(-γt) con la posición x

sin(ωt+φ)=±1, t n ' = π/2 +nπφ ω

Representamos la posición x en función del tiempo t, la amplitud A0exp(-γt), los instantes t'n en los que ambas curvas son tangentes (puntos de color azul) y los instantes tn en los que el desplazamiento x es máximo o mínimo (puntos de color rojo). Para apreciar la diferencia entre ambos, se ha establecido el valor de la constante de amortiguamiento en γ=30.

w0=100; %frecuencia propia
g=30; %constante amortiguamiento
x0=5; %condiciones iniciales
v0=0; 
A0=sqrt((v0^2+2*v0*g*x0+x0^2*w0^2)/(w0^2-g^2));
phi=atan(x0*sqrt(w0^2-g^2)/(v0+g*x0));
w=sqrt(w0^2-g^2);
x=@(t) A0*exp(-g*t).*sin(w*t+phi);
hold on
fplot(x,[0,0.2])
fplot(@(t) A0*exp(-g*t),[0,0.2])
fplot(@(t) -A0*exp(-g*t),[0,0.2])
for n=0:4
    t0=(atan(w/g)+n*pi-phi)/w;
    plot(t0,x(t0),'o','markersize',3,'markeredgecolor','r','markerfacecolor','r')
    t1=(pi/2+n*pi-phi)/w;
    plot(t1,x(t1),'o','markersize',3,'markeredgecolor','b','markerfacecolor','b')
end

hold off
grid on
xlabel('t')
ylabel('x')
title('oscilaciones amortiguadas')

La amplitud A0exp(-γt) no pasa por los máximos, pero muy cerca de los mismos para la mayor parte de las situaciones en las que el rozamiento es pequeño, γ<<ω0. La diferencia t'n-tn es muy pequeña

t n ' t n = π 2 arctan( ω γ ) ω

La energía del oscilador amortiguado

La energía de la partícula que describe una oscilación amortiguada es la suma de la energía cinética de la partícula y de la energía potencial del muelle elástico deformado.

E= 1 2 m v 2 + 1 2 k x 2 = 1 2 m v 2 + 1 2 m ω 0 2 x 2

Introducimos las expresiones de la posición x y de la velocidad v de la partícula en función del tiempo t.

E= 1 2 m ω 0 2 A 2 e 2γt 1 2 mγω A 2 e 2γt sin( 2(ωt+φ) )

Si la constante de amortiguamiento γ es pequeña, como hemos visto en el ejemplo del apartado anterior ω0≈ω

E= 1 2 m ω 0 2 A 2 e 2γt ( 1 γ ω 0 sin( 2( ω 0 t+φ) ) )

La energía decrece exponencialmente con el tiempo, pero con una pequeña ondulación debida al segundo término entre paréntesis.

Añadimos al script anterior las siguientes líneas para representar la energía del oscilador amortiguado en función del tiempo. La energía del oscilador decrece rápidamente con el tiempo.

%energía
v=diff(x,t);
e=0.5*v^2+0.5*100^2*x^2;  %la masa m es un factor de escala m=1
ee=subs(e,{g w0 x0 v0},{7 100 5 0})
figure
ezplot(ee,[0 0.7])
grid on
title('Energía')

La energía perdida hasta el instante t a causa de la fuerza de rozamiento se calcula mediante la integral

0 t (λv )v·dt=2mγ 0 t v 2 dt

que será igual a la diferencia entre la energía del oscilador en el instante t y la energía inicial del oscilador en el instante t=0

Comprobamos con MATLAB

clear
syms t g w0 x0 v0;
w=sqrt(w0^2-g^2);
x=exp(-g*t)*((v0+g*x0)*sin(w*t)/w+x0*cos(w*t)); %posición
v=diff(x,t);  %velocidad
Ei=v0^2/2+w0^2*x0^2/2; %energía inicial del oscilador
Ef=v^2/2+w0^2*x^2/2; %energía final 
DE=Ef-Ei; %diferencia de energías
 %características del oscilador, masa m=1
E1=subs(DE,{g w0 x0 v0},{7 100 5 0});
Ep=-int(2*g*v^2,t,0,t); %trabajo de la fuerza de rozamiento
E2=subs(Ep,{g w0 x0 v0},{7 100 5 0});
subs(E1,t,0.7) %valor para el instante t=0.7 s
subs(E2,t,0.7)

Corremos el script en la ventana de comandos

ans = -1.2499e+005
ans = -1.2499e+005

Valores medios

La posición y velocidad de un oscilador libre, una partícula de masa m unida a un muelle de constante elástica k son

x=Asin( ω 0 t+φ ) ω 0 2 = k m v= dx dt =A ω 0 cos( ω 0 t+φ )

Las energías cinética y potencial son, respectivamente

E p = 1 2 k x 2 = 1 2 k A 2 sin 2 ( ω 0 t+φ ) E k = 1 2 m v 2 = 1 2 m A 2 ω 0 2 cos 2 ( ω 0 t+φ )

La suma de las energías cinética y potencial es constante e independiente del tiempo t

E= 1 2 k x 2 + 1 2 m v 2 = 1 2 m ω 0 2 A 2 sin 2 ( ω 0 t+φ )+ 1 2 m ω 0 2 A 2 cos 2 ( ω 0 t+φ )= 1 2 m ω 0 2 A 2

Se denomina valor medio de una función periódica de periodo P a la integral

f(t) = 1 P 0 P f(t)·dt

fplot(@(t) [sin(t), sin(t).^2],[0,2*pi])
grid on
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'})
legend('sin(t)', 'sin(t)^2', 'location','southwest')
xlabel('t')
ylabel('x')
title('Funciones periódicas')

La función sin2(ω0t+φ) es una función periódica de periodo P=π/ω0

sin 2 ( ω 0 t+φ ) = ω 0 π 0 π/ ω 0 sin 2 ( ω 0 t+φ )·dt= ω 0 2π 0 π/ ω 0 { 1cos( 2 ω 0 t+2φ ) }dt= ω 0 2π { t 1 2 ω 0 sin( 2 ω 0 t+2φ ) } 0 π/ ω 0 = ω 0 2π { π ω 0 1 2 ω 0 sin( 2π+2φ )+ 1 2 ω 0 sin( 2φ ) }= 1 2

Del mismo modo

cos 2 ( ω 0 t+φ ) = ω 0 π 0 π/ ω 0 cos 2 ( ω 0 t+φ )·dt= ω 0 2π 0 π/ ω 0 { 1+cos( 2 ω 0 t+2φ ) }dt= ω 0 2π { t+ 1 2 ω 0 sin( 2 ω 0 t+2φ ) } 0 π/ ω 0 = ω 0 2π { π ω 0 + 1 2 ω 0 sin( 2π+2φ ) 1 2 ω 0 sin( 2φ ) }= 1 2

El valor medio de la energía cinética y potencial son, respectivamente

1 2 m v 2 = E 2 1 2 k x 2 = E 2

La ecuación diferencial del oscilador amortiguado es

m d 2 x d t 2 +λ dx dt +kx=0

Multiplicando por v=dx/dt

mv dv dt +kx dx dt =λ v 2 d dt ( 1 2 m v 2 + 1 2 k x 2 )=λ v 2 dE dt =λ v 2

Tomando valores medios

d E dt = 2λ m 1 2 m v 2 d E dt = 2λ m 1 2 E d E E = λ m dt

Integrando entre 0 y t

E = E 0 exp( 2γt )

El valor medio de la energía en cada periodo, decrece exponencialmente con el tiempo.

Oscilaciones críticas (γ=ω0)

La posición x y velocidad dx/dt en función del tiempo t son

x={ A+Bt }exp( γt ) dx dt =Bexp( γt )γ( A+Bt )exp( γt )

Los coeficientes A y B se determinan a partir de las condiciones iniciales: posición inicial x0 y velocidad inicial v0 en el instante t=0

t=0{ x 0 =A v 0 =BγA

La posición x en función del tiempo t es

x={ x 0 +( v 0 +γ x 0 )t }exp( γt )

Tomamos g=w0=100. Las condiciones iniciales son: x0=5, v0=0

x0=5; %condiciones iniciales
v0=0;
g=100; %constante de amortiguamiento
x=@(t) (x0+(v0+g*x0)*t).*exp(-g*t);
fplot(x,[0,0.02*pi])
title('Oscilación crítica')
xlabel('t')
ylabel('x');
grid on

Oscilaciones sobreamortiguadas (γ>ω0)

La posición x y velocidad dx/dt en función del tiempo t son

x=exp( γt ){ Asinh(ωt)+Bcosh(ωt) } dx dt =γexp( γt ){ Asinh(ωt)+Bcosh(ωt) }+ωexp( γt ){ Acosh(ωt)+Bsinh(ωt) }

Los coeficientes A y B se determinan a partir de las condiciones iniciales: posición inicial x0 y velocidad inicial v0 en el instante t=0

t=0{ x 0 =B v 0 =γB+ωA

La posición x en función del tiempo t es

x=exp( γt ){ x 0 cosh(ωt)+( v 0 +γ x 0 ω )sinh(ωt) }

Tomamos γ=120 y ω0=100. Las condiciones iniciales son: x0=5, v0=0

x0=5; %condiciones iniciales
v0=0;
g=120; %constante de amortiguamiento
w0=100; %frecuencia natural
w=sqrt(g^2-w0^2);
x=@(t) (x0*cosh(w*t)+(v0+g*x0)*sinh(w*t)/w).*exp(-g*t);
fplot(x,[0,0.02*pi])
title('Oscilación sobreamortiguada')
xlabel('t')
ylabel('x');
grid on

Solución numérica

Expresamos la ecuación diferencial de segundo orden que describe las oscilaciones amortiguadas en forma de dos ecuaciones diferenciales de primer orden, para resolverlas utilizando la función ode45 de MATLAB

d 2 x d t 2 +2γ dx dt + ω 0 2 x=0 { dx dt =v dv dt =2γv ω 0 2 x

w0=100; %frecuencia angular propia
g=7; %rozamiento, gamma,
x0=[5,0];%posición inicial, x0, velocidad inicial, v0
tf=[0,0.3*pi]; %intervalo de tiempo

f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)]; 
tspan=[0 tf];
[t,x]=ode45(f,tf,x0);
%plot(t,x(:,1)) %tiempo, posición
plot(x(:,1),x(:,2)) %posición, velocidad (espacio de las fases)
xlabel('x')
ylabel('v');
title('oscilador amortiguado')
grid on

Con la sentencia anulada plot(t,x(:,1)), obtenemos una gráfica similar a la primera figura de esta página

El espacio de las fases nos muestra otra perspectiva del comportamiento del oscilador. Se representa el momento lineal (o la velocidad) v en el eje vertical, y la posición del móvil x en el eje horizontal. La trayectoria es una espiral que converge hacia el origen.

Actividades

Se introduce

Se pulsa el botón titulado Nuevo.

Probar con los siguientes valores de la constante de amortiguamiento γ : 5 (amortiguadas), 100 (críticas), 110 (sobreamortiguadas).


Medida de la constante de amortiguamiento

El decrecimiento logarítmico δ nos permite calcular fácilmente la constante de amortiguamiento γ, cuando su valor es pequeño comparado con la frecuencia natural del oscilador ω0.

A 1 A 2 = exp(γt) exp( γ(t+P) ) =exp(γP) δ=ln A 1 A 2 =γP= 2πγ ω 0 2 γ 2 =2π ( γ/ ω 0 ) 1 ( γ/ ω 0 ) 2 γ ω 0 δ2π γ ω 0

En una experiencia hemos medido los máximos desplazamientos

δ=ln(A1/A2)=2πγ/ω0, como ω0=100 rad/s, entonces γ=7.02.

En una experiencia de laboratorio con el péndulo de Pohl se han tomado los siguientes datos de la amplitud de una oscilación amortiguada

t (s) A
0 20
0.905 15.8
1.810 13.3
2.715 10.5
3.620 8.7
4.525 6.8
5.430 5.7
6.335 4.4
7.240 3.4
8.145 2.8

Si la amplitud disminuye exponencialmente con el tiempo, A=A0exp(-γt), o bien, lnA=lnA0-γt que es la ecuación de una línea recta

Representamos ln(A) en función del tiempo t y realizamos un ajuste lineal de los datos. La pendiente de la recta será la constante γ de amortiguamiento.

t=0.905*(0:9); %tiempo
A=[20,15.8,13.3,10.5,8.7,6.8,5.7,4.4,3.4,2.8]; %amplitud
plot(t,log(A),'ro','markersize',4,'markerfacecolor','r')
grid on
xlabel('t')
ylabel('log(A')
title('Oscilaciones amortiguadas')

En la ventana gráfica de MATLAB elegimos Tools/Basic fitting, activamos la casilla linear y obtenemos para la pendiente de la recta que mejor ajusta a los datos experimentales p1=-0.24148. Por tanto, la constante γ=0.24 s-1