Movimiento de un péndulo bajo la acción de una fuerza de rozamiento proporcional al cuadrado de la velocidad

Estudiamos ahora, el péndulo simple cuyo comportamiento difiere del oscilador consistente en una masa unida a un muelle elástico

Cuando un péndulo simple oscila en el aire, es una suposición razonable considerar que la fuerza de rozamiento que actúa sobre la partícula sea proporcional al cuadrado de la velocidad.

El péndulo describe una trayectoria circular, un arco de una circunferencia de radio l (la longitud del péndulo).

Las fuerzas que actúan sobre la partícula de masa m son:

Descomponemos el peso en la acción simultánea de dos componentes, mg·sinθ en la dirección tangencial y mg·cosθ en la dirección radial.

La aceleración tangencial de la partícula es at=dv/dt.

La relación entre la aceleración tangencial at y la aceleración angular α es at=α·l=l·d2θ/dt2. La ecuación del movimiento se escribe en forma de ecuación diferencial

d 2 θ d t 2 + g l sinθ+γl ( dθ dt ) 2 =0{ γ>0 dθ dt >0 γ<0 dθ dt <0

Solució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 dθ/dt=0 desde θ=θ0

w0=100; %frecuencia angular rad/s
k=0.7; % gamma·l, constante de amortiguamiento
x0=pi/3; %ángulo inicial

f=@(t,x) [x(2);-k*x(2)^2*sign(x(2))-w0^2*sin(x(1))]; 
[t,x]=ode45(f,[0,1],[x0,0]);
plot(t, x(:,1))
grid on
xlabel('t')
ylabel('\theta');
title('Oscilador amortiguado')

Un resultado similar al del oscilador lineal

La mayor riqueza de comportamientos se produce para la condición inicial siguiente: en la posición de partida θ0, la velocidad de la partícula es v0, que trataremos en el siguiente sección

Solución analítica

Transformamos la ecuación diferencial de segundo orden en θ y t en una de primer orden en v=l·dθ/dt y θ donde v es la velocidad de la partícula

1 l dv dt + g l sinθ+γ v 2 l =0 dv dt = dv dθ dθ dt = v l dv dθ = 1 2l d v 2 dθ d v 2 dt +2γl v 2 =2glsinθ

Llamando z=v2, la solución particular es z1=Asinθ+Bcosθ que introducimos en la ecuación diferencial para calcular los coeficientes A y B

z 1 = 2gl 1+4 γ 2 l 2 ( cosθ2γlsinθ )

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

v 2 =Cexp( 2γlθ )+ 2gl 1+4 γ 2 l 2 ( cosθ2γlsinθ )

El coeficiente C se determina a partir de las condiciones iniciales: en la posición inicial θ0, la partícula lleva una velocidad v0.

v 0 2 =Cexp( 2γl θ 0 )+ 2gl 1+4 γ 2 l 2 ( cos θ 0 2γlsin θ 0 )

La velocidad v de la partícula en la posición angular θ es

v 2 = v 0 2 exp( 2γl(θ θ 0 ) )+ 2gl 1+4 γ 2 l 2 { cosθ2γlsinθ( cos θ 0 2γlsin θ 0 )exp( 2γl(θ θ 0 ) ) }

La tensión del hilo es

T mg =2 v 0 2 2gl exp( 2γl(θ θ 0 ) )+ 2 1+4 γ 2 l 2 { cosθ2γlsinθ( cos θ 0 2γlsin θ 0 )exp( 2γl(θ θ 0 ) ) }+cosθ

Si no hay rozamiento γ=0,

v 2 = v 0 2 +2gl( cosθcos θ 0 )

Ecuación que se obtiene aplicando el principio de conservación de la energía

El estudio de los distintos casos es similar al de la página titulada Rozamiento en el bucle

Posición angular final, θ<π/2

Cuando la posición angular final θ<π/2, es el caso más simple de analizar ya que la tensión del hilo T es siempre positiva

Calculamos la velocidad inicial v0,π/2 que deberíamos impartir a la partícula en la posición inicial θ0 para que llegue justamente a la posición θ=π/2, (con velocidad nula)

0= v 0,π/2 2 exp( 2γl( π 2 θ 0 ) )+ 2gl 1+4 γ 2 l 2 { 2γl( cos θ 0 2γlsin θ 0 )exp( 2γl( π 2 θ 0 ) ) } v 0,π/2 2 = 2gl 1+4 γ 2 l 2 { 2γlexp( γl(π2 θ 0 ) )+cos θ 0 2γlsin θ 0 }

Una expresión más simple se obtiene cuando la posición inicial θ0=0, el péndulo parte del origen

v 0,π/2 2 2gl = 1 1+4 γ 2 l 2 { 2γlexp( γlπ )+1 }

Si la velocidad inicial v0 es menor, entonces la posición final del péndulo θ1<π/2, que vamos a calcular resolviendo la ecuación transcendente por procedimiento fzero de MATLAB

0= v 0 2 exp( 2γl(θ θ 0 ) )+ 2gl 1+4 γ 2 l 2 { cosθ2γlsinθ( cos θ 0 2γlsin θ 0 )exp( 2γl(θ θ 0 ) ) } v 0 2 2gl ( 1+4 γ 2 l 2 )+( cosθ2γlsinθ )exp( 2γl(θ θ 0 ) )( cos θ 0 2γlsin θ 0 )=0

Supongamos que el péndulo parte del origen θ0=0, con velocidad v0 menor que la máxima v0, π/2. Sea k=γl=0.25, la velocidad máxima v0,π/2=1.2951 en unidades 2gl . Si el péndulo parte del origen con velocidad v0=0.25, llega a la posición final θ1=70.38°

k=0.25; %rozamiento
v_max=sqrt((1+2*k*exp(k*pi))/(1+4*k^2)); %velocidad inicial para pi/2
disp(v_max)
v0=1; %velocidad inicial
f=@(x) v0^2*(1+4*k^2)+(cos(x)-2*k*sin(x))*exp(2*k*x)-1;
th_1=fzero(f, [0, pi/2]); %posición final
disp(th_1*180/pi)
    1.2951
   70.3839

Representamos la posición final θ1 en función de la velocidad inicial v 0 2gl , para varios valores de k=γl. La posición inicial θ0=0 (el péndulo parte del origen)

hold on
for k=[0,0.05,0.1,0.25] %rozamiento
    v_max=sqrt((1+2*k*exp(k*pi))/(1+4*k^2)); %velocidad inicial para pi/2
    vv=linspace(0,v_max,100); 
    th=zeros(1,length(vv));
    th(end)=pi/2;
    i=1;
    for v0=vv(1:end-1)
        f=@(x) v0^2*(1+4*k^2)+(cos(x)-2*k*sin(x))*exp(2*k*x)-1;
        th(i)=fzero(f, [0, pi/2]);
        i=i+1;
    end
    plot(vv,th, 'displayName',num2str(k)),
end
hold off
set(gca,'YTick',0:pi/12:pi/2)
set(gca,'YTickLabel',{'0','\pi/12','\pi/6','\pi/4','\pi/3', '5\pi/12','\pi/2'})
xlabel('v_0')
ylabel('\theta')
grid on
legend('-DynamicLegend','location','best')
title('Posición angular final')

Posición angular final, π/2<θ

Para que la partícula pase por la posición más alta de su trayectoria θ=π, debe llevar una velocidad mínima. Aplicando la ecuación de la dinámica del movimiento circular

m v 2 l =mg+T

La velocidad mínima, se obtiene cuando la tensión del hilo T es cero

v= gl

La velocidad inicial v0,π que deberá llevar la partícula para que llegue a la posición final θ=π con esta velocidad mínima v2=gl es

gl= v 0,π 2 exp( 2γl(π θ 0 ) )+ 2gl 1+4 γ 2 l 2 { 1( cos θ 0 2γlsin θ 0 )exp( 2γl(π θ 0 ) ) } v 0,π 2 = 2gl 1+4 γ 2 l 2 { cos θ 0 2γlsin θ 0 +exp( 2γl(π θ 0 ) ) }+glexp( 2γl(π θ 0 ) )

Si la velocidad inicial de la partícula es inferior a v0,π pero superior a v0,π/2, la posición angular final θ1>π/2 pero θ1<π. La posición final se obtiene cuando la tensión T del hilo se hace cero o bien, cuando la velocidad de la partícula vale (véase la figura)

m v 2 l =mgcos θ 1 v 2 =glcos θ 1

A partir de esta posición, la partícula describe una trayectoria parabólica

La posición final θ1 se calcula resolviendo la ecuación transcendente por el procedimiento fzero de MATLAB

v 0 2 gl exp( 2γl(θ θ 0 ) )+ 2 1+4 γ 2 l 2 { cosθ2γlsinθ( cos θ 0 2γlsin θ 0 )exp( 2γl(θ θ 0 ) ) }+cosθ=0 v 0 2 2gl = 1 1+4 γ 2 l 2 { cos θ 0 2γlsin θ 0 ( cosθ2γlsinθ )exp( 2γl(θ θ 0 ) ) } 1 2 cosθexp( 2γl(θ θ 0 ) )

Supongamos que el péndulo parte del origen θ0=0, con velocidad v0 mayor que v0,π/2 y menor que v0,π. Sea k=γl=0.25, la velocidad mínima v0,π/2=1.2951 y la máxima v0,π=2.6559 en unidades 2gl . Si el péndulo parte del origen con velocidad v0=1.8, llega a la posición final θ1=114.31°

k=0.25; %rozamiento
v_min=sqrt((1+2*k*exp(k*pi))/(1+4*k^2)); %velocidad inicial para pi/2
v_max=sqrt((1+exp(2*k*pi))/(1+4*k^2)+exp(2*k*pi)/2); %velocidad inicial para pi
disp([v_min,v_max])
v0=1.8;
f=@(x) (1-(cos(x)-2*k*sin(x))*exp(2*k*x))/(1+4*k^2)-exp(2*k*x)*cos(x)/2-v0^2;
th_1=fzero(f, [pi/2, pi]);
disp(th_1*180/pi) 
    1.2951    2.6559
  114.3149

Representamos la posición final θ1 en función de la velocidad inicial v 0 2gl , para varios valores de k=γl. La posición inicial θ0=0 (el péndulo parte del origen)

hold on
for k=[0,0.05,0.1,0.25]
    v_min=sqrt((1+2*k*exp(k*pi))/(1+4*k^2)); %velocidad inicial para  pi/2
    v_max=sqrt((1+exp(2*k*pi))/(1+4*k^2)+exp(2*k*pi)/2); %v. inicial para  pi
    vv=linspace(v_min,v_max,100); 
    th=zeros(1,length(vv));
    th(1)=pi/2;
    th(end)=pi;
    i=2;
    for v0=vv(2:end-1)
        f=@(x) (1-(cos(x)-2*k*sin(x))*exp(2*k*x))/(1+4*k^2)-exp(2*k*x)*cos(x)/2-
v0^2;
        th(i)=fzero(f, [pi/2, pi]);
        i=i+1;
    end
    plot(vv,th, 'displayName',num2str(k)),
end
hold off
set(gca,'YTick',pi/2:pi/12:pi)
set(gca,'YTickLabel',{'\pi/2','7\pi/12','2\pi/3','3\pi/4','5\pi/6', '11\pi/12',
'\pi'})
xlabel('v_0')
ylabel('\theta')
grid on
legend('-DynamicLegend','location','best')
title('Posición angular final')

Unimos las dos porciones de código, para representar la posición final θ1 en función de la velocidad inicial para varios valores del parámetro k: 0 (negro), 0.05 (azul), 0.1 (rojo), 0.25 (verde)

hold on
colores=['k','b','r','g'];
indice=1;
for k=[0,0.05,0.1,0.25]
    v_max=sqrt((1+2*k*exp(k*pi))/(1+4*k^2)); %velocidad inicial para pi/2
    vv=linspace(0,v_max,100); 
    th=zeros(1,length(vv));
    th(end)=pi/2;
    i=1;
    for v0=vv(1:end-1)
        f=@(x) v0^2*(1+4*k^2)+(cos(x)-2*k*sin(x))*exp(2*k*x)-1;
        th(i)=fzero(f, [0, pi/2]);
        i=i+1;
    end
    plot(vv,th, colores(indice))

    v_min=sqrt((1+2*k*exp(k*pi))/(1+4*k^2)); %velocidad inicial para  pi/2
    v_max=sqrt((1+exp(2*k*pi))/(1+4*k^2)+exp(2*k*pi)/2); %v. inicial para  pi
    vv=linspace(v_min,v_max,100); 
    th=zeros(1,length(vv));
    th(1)=pi/2;
    th(end)=pi;
    i=2;
    for v0=vv(2:end-1)
        f=@(x) (1-(cos(x)-2*k*sin(x))*exp(2*k*x))/(1+4*k^2)-exp(2*k*x)*cos(x)/2-
v0^2;
        th(i)=fzero(f, [pi/2, pi]);
        i=i+1;
    end
    plot(vv,th, colores(indice))
    indice=indice+1;
end
hold off
set(gca,'YTick',0:pi/6:pi)
set(gca,'YTickLabel',{'0','\pi/6','\pi/3','\pi/2','2\pi/3','5\pi/6','\pi'})
xlabel('v_0')
ylabel('\theta')
grid on
title('Posición angular final')

Posición angular final, π<θ<2π

Si la velocidad inicial v0 es mayor que v0,π, la tensión del hilo T en la posición más elevada θ=π es positiva y el movimiento del péndulo se extiende más allá, hasta que la tensión T se anule en el intervalo π<θ<3π/2.

Por ejemplo, para k=0.5, v0,π=4.8622. Representamos la tensión T/(mg) del hilo en función del ángulo θ en el intervalo (π, 3π/2) para una velocidad inicial un poco mayor, v0=5.07

k=0.5; %rozamiento
v_max=sqrt((1+exp(2*k*pi))/(1+4*k^2)+exp(2*k*pi)/2); %velocidad inicial para pi
disp(v_max)
v0=5.07; %velocidad inicial
f=@(x) 2*v0^2*exp(-2*k*x)+2*(cos(x)-2*k*sin(x)-exp(-2*k*x))/(1+4*k^2)+cos(x);
hold on
fplot(f, [pi/2,3*pi/2])
line([2*pi/3,3*pi/2],[0,0],'color','k')
plot(198*pi/180,0,'o','markersize',3,'markeredgecolor','r','markerfacecolor','r')
line([198*pi/180,198*pi/180],[-1,0],'color','k','lineStyle','--')
hold off
set(gca,'XTick',pi/2:pi/6:3*pi/2)
set(gca,'XTickLabel',{'\pi/2','2\pi/3','5\pi/6','\pi','7\pi/6', '4\pi/3','3\pi/2'})
xlabel('\theta')
ylabel('T/mg')
grid on
ylim([-1,1])
title('Tensión del hilo')

4.8622

Observamos que la tensión del hilo es positiva, toma un valor próximo a cero para un ángulo cercano a θc=198°. Para este ángulo crítico, la tensión T y su derivada dT/dθ (pendiente de la recta tangente) se anulan simultáneamente, tal como apreciamos en la figura. El péndulo continuará su movimiento hasta alcanzar el punto más bajo de su trayectoria θ=2π

Para hacer los cálculos más sencillos, supongamos que partimos del origen θ0=0. La tensión del hilo es

T mg =2 v 0 2 2gl exp( 2γlθ )+ 2 1+4 γ 2 l 2 { cosθ2γlsinθexp( 2γlθ ) }+cosθ

Su derivada es

1 mg ( dT dθ )=2 v 0 2 2gl ( 2γl )exp( 2γlθ )+ 2 1+4 γ 2 l 2 { sinθ2γlcosθ+2γlexp( 2γlθ ) }sinθ

Resolvemos el sistema de dos ecuaciones con dos incógnitas, el ángulo crítico θc y la velocidad inicial crítica v0,2π

{ 2 v 0 2 2gl exp( 2γlθ )+ 2 1+4 γ 2 l 2 { cosθ2γlsinθexp( 2γlθ ) }+cosθ=0 2 v 0 2 2gl ( 2γl )exp( 2γlθ )+ 2 1+4 γ 2 l 2 { sinθ2γlcosθ+2γlexp( 2γlθ ) }sinθ=0

Eliminamos v0 y llegamos a la ecuación

tan θ c = 2 3 k

que es un ángulo del tercer cuadrante. Una vez calculado este ángulo θc despejamos la velocidad inicial crítica, v0,2π, tomando la primera del sistema de dos ecuaciones

v 0,2π 2 2gl = 1 1+4 γ 2 l 2 { 1( cos θ c 2γlsin θ c )exp( 2γl θ c ) } 1 2 cos θ c exp( 2γl θ c )

>> k=0.5;
>> th_c=pi+atan(2*k/3) %ángulo del tercer cuadrante
th_c =    3.4633

>> v_c=sqrt(((2*k*sin(th_c)-cos(th_c))*exp(2*k*th_c)+1)/(1+4*k^2)-
exp(2*k*th_c)*cos(th_c)/2)
v_c =    5.0732 %velocidad inicial crítica

>> f(3.4633) %tensión del hilo
ans =  -2.5339e-05

Posición angular final, θ≥2π

Cuando la velocidad inicial en el origen θ0=0, es v0v0,2π el péndulo completa una vuelta θ=2π y regresa a la posición de partida. La nueva velocidad inicial de la partícula es

v 2 = v 0 2 exp( 2γl·2π )+ 2gl 1+4 γ 2 l 2 { 1exp( 2γl·2π ) }

Con esta velocidad inicial, hay que volver a analizar el movimiento del péndulo para el valor dado de k=γ·l

Oscilaciones amortiguadas

El caso más importante se produce cuando la posición final θ1<π/2. El primero tratado en esta página.

  1. El péndulo se mueve en el sentido contrario a las agujas del reloj, v>0.

  2. Como ya se ha analizado, la ecuación diferencial que describe el movimiento desde θ0 con velocidad v0 hasta θ1 en reposo es

    d v 2 dt +2γl v 2 =2glsinθ

    La solución, con la condición inicial: en la posición θ0 la velocidad inicial es v0

    v 2 = v 0 2 exp( 2γl(θ θ 0 ) )+ 2gl 1+4 γ 2 l 2 { cosθ2γlsinθ( cos θ 0 2γlsin θ 0 )exp( 2γl(θ θ 0 ) ) }

    El péndulo se detiene para el ángulo θ1, solución de la ecuación transcendente

    v 0 2 2gl + 1 1+4 γ 2 l 2 { ( cosθ2γlsinθ )exp( 2γl(θ θ 0 ) )( cos θ 0 2γlsin θ 0 ) }=0

  3. El péndulo se mueve en el sentido de las agujas del reloj, v<0.

  4. Cuando el péndulo inicial el camino de vuelta, la velocidad cambia de signo v<0 y la fuerza de rozamiento también, la ecuación del movimiento es

    d v 2 dt 2γl v 2 =2glsinθ

    Llamando z=v2, la solución particular es z1=Asinθ+Bcosθ que introducimos en la ecuación diferencial para calcular los coeficientes A y B

    z 1 = 2gl 1+4 γ 2 l 2 ( cosθ+2γlsinθ )

    La solución de la homogénea es z2=C·exp(2γlθ). La solución completa es v2=z1+z2

    v 2 =Cexp( 2γlθ )+ 2gl 1+4 γ 2 l 2 ( cosθ+2γlsinθ )

    Es la misma solución deducida anteriormente, sustituyendo γ por -γ

    El coeficiente C se determina a partir de las condiciones iniciales: en la posición inicial θ1, el móvil parte del reposo.

    C= 2gl 1+4 γ 2 l 2 ( cos θ 1 +2γlsin θ 1 )exp( 2γl θ 1 ) v 2 = 2gl 1+4 γ 2 l 2 { cosθ+2γlsinθ( cos θ 1 +2γlsin θ 1 )exp( 2γl(θ θ 1 ) ) }

    El péndulo se detiene para el ángulo θ2, solución de la ecuación transcendente

    cosθ+2γlsinθ( cos θ 1 +2γlsin θ 1 )exp( 2γl(θ θ 1 ) )=0

  5. El péndulo se mueve en el sentido contrario a las agujas del reloj, v>0.

  6. El péndulo parte en reposo desde la posición θ2, su velocidad es v>0,

    v 2 = 2gl 1+4 γ 2 l 2 { cosθ2γlsinθ( cos θ 2 2γlsin θ 2 )exp( 2γl(θ θ 2 ) ) }

    El péndulo se detiene para el ángulo θ3, solución de la ecuación transcendente

    cosθ2γlsinθ( cos θ 2 2γlsin θ 2 )exp( 2γl(θ θ 2 ) )=0

  7. El péndulo se mueve en el sentido de las agujas del reloj, v<0.

  8. y así, sucesivamente

Representamos el movimiento del péndulo en el espacio de las fases para k=0.5. El péndulo parte del origen θ0=0, con velocidad inicial v0<v0,π/2=1.2951

k=0.25; %rozamiento
v_max=sqrt((1+2*k*exp(k*pi))/(1+4*k^2)); %velocidad inicial para pi/2
disp(v_max)
v0=1.2; %velocidad inicial
f=@(x)  v0^2*(1+4*k^2)+(cos(x)-2*k*sin(x))*exp(2*k*x)-1;
th_1=fzero(f,[0,pi/2]); %posición angular final
th=linspace(0,th_1,100);
xx=th(2:end-1);
v=@(x) sqrt(v0^2*exp(-2*k*x)+(cos(x)-2*k*sin(x)-exp(-2*k*x))/(1+4*k^2));
y=[v0,v(xx),0]; 
hold on
plot(th,y);
for i=1:10
    f=@(x) cos(x)+2*k*sin(x)-(cos(th_1)+2*k*sin(th_1))*exp(2*k*(x-th_1));
    th_2=fzero(f,[0,-th_1]);
    th=linspace(th_2,th_1,100);
    xx=th(2:end-1);
    %los puntos de retorno x(1) y x(end) tienen velocidad cero
    v=@(x) sqrt((cos(x)+2*k*sin(x)-(cos(th_1)+2*k*sin(th_1))*
exp(2*k*(x-th_1)))/(1+4*k^2));
    y=[0,v(xx),0]; 
    plot(th,-sign(k)*y);
    th_1=th_2;
    k=-k;
end
hold off
grid on
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'})
xlabel('\theta')
ylabel('v')
title('Espacio de las fases')

1.2951

Actividades

Se introduce

Se pulsa el botón titulado Nuevo y se observa el movimiento del péngulo

En la parte izquierda, se proporcionanlos datos del tiempo t, la posición angular θ, en grados, la velocidad angular ω=dθ/dt, y la tensión del hilo T

Cuando la velocidad v es cero, el cuerpo se detiene. Cuando la tensión del hilo T es nula, el cuerpo describe una trayectoria parabólica.

En la parte derecha, se representa el balance energético: la energía cinética (en color rojo), la energía potencial (en color azul) y el trabajo de la fuerza de rozamiento (en color negro)


Referencias

Marko V Lubarda, Vlado A Lubarda. An analysis of pendulum motion in the presence of quadratic and linear drag. Eur. J. Phys. 42 (2021) 055014