Oscilador amortiguado por una fuerza proporcional al cuadrado de la velocidad

Masa unida a un muelle elástico

Las fuerzas que actúan sobre el móvil son: la fuerza elástica F=-kx, y la fuerza de rozamiento de sentido contrario a la velocidad Fr=-λv2, donde λ es una constante que depende del sistema físico particular.

Dado que v2 es positivo, λ será positiva si v>0 y negativa si v<0

La ecuación del movimiento se escribe

ma=-kx-λv2

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 +γ ( dx dt ) 2 + ω 0 2 x=0{ γ>0 dx dt >0 γ<0 dx dt <0 ω 0 2 = k m γ= λ m

Solucción numérica

Resolvemos la ecuación diferencial por procedimientos numéricos, con las siguientes condiciones iniciales: en el instante t=0, el móvil parte en reposo dx/dt=0 desde x=x0

w0=100; %frecuencia angular rad/s
g=0.7; %constante en 1/m
x0=1; %posición inicial
f=@(t,x) [x(2);-g*x(2)^2*sign(x(2))-w0*w0*x(1)]; 
[t,x]=ode45(f,[0,1],[x0,0]);
plot(t, x(:,1))
grid on
xlabel('t')
ylabel('x');
title('Oscilador amortiguado')

Solución analítica

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.

Transformamos la ecuación diferencial de segundo orden en x y t en una de primer orden en v y x

dv dt +γ v 2 + ω 0 2 x=0{ γ>0v>0 γ<0v<0 v dv dx +γ v 2 + ω 0 2 x=0 1 2 dz dx +γz= ω 0 2 xz= v 2

La solución particular es z1=Ax+B que introducimos en la ecuación diferencial para calcular los coeficientes A y B

z 1 = ω 0 2 γ ( 1 2γ x )

La solución de la homogénea es z2=C·exp(-2γx). La solución completa es z=z1+z2

z=C e 2γx ω 0 2 γ x+ ω 0 2 2 γ 2

El coeficiente C se determina a partir de las condiciones iniciales: en la posición inicial x0, el móvil parte del reposo, v=0

v 2 = ω 0 2 γ { ( x 0 1 2γ ) e 2γ(x x 0 ) ( x 1 2γ ) }

Vamos a representar la velocidad v en función de la posición x del móvil

Primera etapa

Partimos de x0>0 con velocidad inicial v=0, la partícula se mueve hacia la izquierda v<0 (γ<0) hasta que se detiene en x1. Esta posición es la raíz de la ecuación transcendente

( x 0 1 2γ ) e 2γ(x x 0 ) ( x 1 2γ )=0

Una vez calculado x1 se representa la función v(x) en el intervalo (x1, x0)

v=sgn(γ)· ω 0 1 γ { ( x 0 1 2γ ) e 2γ(x x 0 ) ( x 1 2γ ) }

donde la función sgn(γ) indica que si γ<0 entonces v<0

Segunda etapa

Ahora x1 es la posición de partida y la partícula, se mueve hacia la derecha v>0 y por tanto, γ>0. La posición x2 se calcula resolviendo la ecuación transcendente con γ cambiado de signo

Una vez calculado x2 se representa la función v(x) en el intervalo (x1, x2), donde la función sgn(γ) indica que si γ>0 entonces v>0

Se ha completado el primer periodo, se repite las dos etapas para los siguientes. El código MATLAB que calcula 10 semiperiodos es el siguiente

w0=100; %frecuencia angular rad/s
g=0.7; %constante de amortiguamiento
x0=1; %posición inicial

hold on
for i=1:10
    g=-g;
    f=@(x) (x0-1/(2*g))*exp(-2*g*(x-x0))-(x-1/(2*g));
    x1=fzero(f,-x0);
    x=linspace(x1,x0,50);
    xx=x(2:end-1);
    %los puntos de retorno x(1) y x(end) tienen velocidad cero
    y=[0,w0*sqrt(f(xx)/g),0]; 
    plot(x,sign(g)*y);
    x0=x1;   
end
hold off
grid on
xlabel('x')
ylabel('v')
title('Espacio de las fases')

Dispositivo experimental

El dispositivo experimental es un péndulo formado por una varilla de masa m y longitud L unida a un disco de masa M y radio R, perpendicular a la dirección de la velocidad, tal como se muestra en la figura.

El momento de inercia respecto de un eje perpendicular al plano de oscilación y que pasa por O es la suma de los momentos de inercia

Aplicando el teorema de Steiner

I=( 1 12 m L 2 +m ( L 2 ) 2 )+( 1 4 M R 2 +M L 2 )= 1 3 m L 2 + 1 4 M R 2 +M L 2

La posición del dentro de masas, su distancia al centro O es

h= m L 2 +ML m+M

Las fuerzas sobre el péndulo son

La variación de energía (cinética más potencial) entre los instantes t y t+dt, el péndulo se desplaza y su cambio de velocidad es , es igual al trabajo de la fuerza de rozamiento entre dichos instantes -Mr·dθ.

La relación entre estas tres magnitudes es

F r L·dθ=d E k +d E p 1 2 C D ρ f π R 2 ( ωL ) 2 L·dθ= 1 2 Id ω 2 +( m+M )ghsinθ·dθ

Dividiendo por I/2, llamando z=ω2 y definiendo los parámetros a y b del siguiente modo

b= C D ρ f π R 2 L 3 1 3 m L 2 + 1 4 M R 2 +M L 2 = C D ρ f π R 2 L 1 3 m+ 1 4 M R 2 L 2 +M a= 2( m+M )gh 1 3 m L 2 + 1 4 M R 2 +M L 2 = m+2M 1 3 m+ 1 4 M R 2 L 2 +M g L

Obtenemos la ecuación diferencial

bz·dθ=dz+asinθ·dθ dz dθ +bz=asinθ

Teniendo en cuenta que

e bθ d dθ ( z e bθ )=bz+ dz dθ

La ecuación diferencial se integra fácilmente, suponiendo que el péndulo se suelta ω=0, z=0 en el instante t=0, desde la posición θ0

e bθ d dθ ( z e bθ )=asinθ d dθ ( z e bθ )=a e bθ sinθ ω 2 e bθ =a θ 0 θ e bθ sinθ

Integrando por partes dos veces

e bθ sinθ·dθ ,{ u= e bθ ,du=b e bθ ·dθ dv=sinθ·dθ,v=cosθ = e bθ cosθ+b e bθ cosθ·dθ ,{ u= e bθ ,du=b e bθ ·dθ dv=cosθ·dθ,v=sinθ = e bθ cosθ+b{ e bθ sinθb e bθ sinθ·dθ } e bθ sinθ·dθ = e bθ ( bsinθcosθ ) 1+ b 2

El cuadrado de la velocidad angular vale

ω 2 e bθ = a e bθ ( cosθbsinθ ) 1+ b 2 | θ 0 θ ω 2 = a 1+ b 2 { cosθbsinθ e b( θ 0 θ ) ( cos θ 0 bsin θ 0 ) }

La posición de retorno ω=0, es la raíz θ1 de la ecuación trascendente

cosθbsinθ e b( θ 0 θ ) ( cos θ 0 bsin θ 0 )=0

Al cambiar el signo de la velocidad angular ω después de cada posición de retorno ω=0, debe cambiar el signo del parámetro b, con el fin de asegurar que la fuerza de rozamiento es de sentido opuesto a la velocidad

Ejemplo

Con estos datos obtenemos: a=56.9543, b=0.3190 s-2

Sea la posición inicial θ0=-π/2 (90°) en el instante t=0

  1. Se calcula la posición final de retorno θ1=1.1334 (65°), resolviendo la ecuación trascendente

  2. cosθbsinθ e b( θ 0 θ ) ( cos θ 0 bsin θ 0 )=0

  3. Se representa la función ω=ω(θ) desde θ0=-π/2 (90°) a θ1=1.1334 (65°)

  4. ω= a 1+ b 2 { cosθbsinθ e b( θ 0 θ ) ( cos θ 0 bsin θ 0 ) }

  5. Se cambia el signo de b

  6. θ0=1.1334 (65°) es la nueva posición inicial

  7. Se calcula la posición final de retorno θ1=-0.8991 (-51.5°) resolviendo la ecuación trascendente

  8. Se representa la función ω=-ω(θ) desde θ0=1.1334 (65°) a θ1=-0.8991 (-51.5°)

  9. Se cambia el signo de b

  10. θ0=-0.8991 (-51.5°) es la nueva posición inicial

Se repiten los ocho pasos...

L=38/100; %longitud varilla
m=12.8/1000; %masa varilla
R=6/100; %radio disco
M=15.2/1000; %masa disco
c_d=1.2; %coeficiente de forma
A=pi*R^2; %área del disco
rho=1.21; %densidad del aire

b=c_d*rho*A*L/(m/3+M*(R/L)^2/4+M); %parámetros
a=(m+2*M)*(9.8/L)/(m/3+M*(R/L)^2/4+M);
th_0=-pi/2; %posición de partida
hold on
for k=1:3
    f=@(x) cos(x)-b*sin(x)-(cos(th_0)-b*sin(th_0))*exp(b*(th_0-x));
    th_1=fzero(f,[0,pi/2]);
    disp(th_1*180/pi)
    g=@(x) sqrt(a*f(x)/(1+b^2));
    fplot(g,[th_0,th_1])
    th_0=th_1;
    b=-b;
    f=@(x) cos(x)-b*sin(x)-(cos(th_0)-b*sin(th_0))*exp(b*(th_0-x));
    th_1=fzero(f,[0,-pi/2]);
    disp(th_1*180/pi)
    g=@(x) -sqrt(a*f(x)/(1+b^2));
   fplot(g,[th_1,th_0])
   th_0=th_1;
   b=-b;
end
hold off
set(gca,'XTick',-pi/2:pi/6:pi/2)
set(gca,'XTickLabel',{'-\pi/2','-\pi/3','-\pi/6','0','\pi/6', '\pi/3','\pi/2'})
grid on
xlabel('\theta')
ylabel('\omega')
title('Oscilaciones amortiguadas de un péndulo')

Referencias

Carl E Mungan, Trevor C Lipscombe. Oscillations of a quadratically damped pendulum. Eur. J. Phys. 34 (2013) pp. 1243–1253