Un cilindro macizo cuya base es una elipse, rueda sin deslizar

Momento de inercia

El momento de inercia de un cilindro de base elíptica de semiejes a y b es igual al momento de inercia de una placa de forma elíptica de las mismas dimensiones respecto de un eje perpendicular a la placa y que pase por su centro.

I c = r 2 dm

Tomamos un elemento diferencial de área dx·dy, que dista r del eje de rotación, r2=x2+y2. La masa dm contenida en este elemento diferencial es

dm= m πab dx·dy

siendo πab el área de la elipse y m la masa de la placa.

La ecuación de una elipse es

x 2 a 2 + y 2 b 2 =1

Fijado x, los límites de y son

y=± b a a 2 x 2

Los límites de x son ±a. La integral doble se escribe

I c = m πab a a b a a 2 x 2 b a a 2 x 2 ( x 2 + y 2 )dx·dy

Integramos respecto de la variable y

I c = m πab a a ( x 2 y+ y 3 3 | b a a 2 x 2 b a a 2 x 2 )dx= 2m πab { b a a a x 2 a 2 x 2 dx+ 1 3 b 3 a 3 a a ( a 2 x 2 ) 3/2 dx }

Hacemos el cambio de variable, x=asint, dx=acost·dt

I c = 2m π a 2 { a 4 π/2 π/2 sin 2 t· cos 2 t·dt+ 1 3 b 2 a 2 a 4 π/2 π/2 cos 4 t·dt }= 2m π a 2 { 1 4 π/2 π/2 sin 2 ( 2t )·dt+ 1 3 b 2 a 2 ( cos 3 t·sint | π/2 π/2 +3 π/2 π/2 sin 2 t· cos 2 t·dt ) }= 2m π a 2 { 1 8 π/2 π/2 ( 1cos(4t) )·dt+ 1 3 b 2 a 2 ( 3 8 π/2 π/2 ( 1cos(4t) )·dt ) }= m 4π a 2 { ( t 1 4 sin(4t) ) | π/2 π/2 + b 2 a 2 ( t 1 4 sin(4t) ) | π/2 π/2 } I c = 1 4 m( a 2 + b 2 )

Cuando a=b, obtenemos el momento de inercia de un disco de radio a respecto de un eje perpendicular al plano del disco y que pasa por su centro

Movimiento de rodar sin deslizar sobre un plano horizontal

Sea un cilindro de base elíptica de masa m y semiejes a y b, está inicialmente en reposo sobre el plano horizontal tal como se muestra en la figura de la izquierda. La longitud del cilindro l no interviene en la descripción del movimiento

Se hace rodar el cilindro hacia la derecha, P es el punto de contacto entre la elipse y el plano horizontal. La distancia entre O y P es la longitud del arco de elipse entre O' y P. La elipse ha girado un ángulo θ, la intersección del eje mayor de la elipse y el plano horizontal forman un ángulo θ.

Las coordenadas (x, y) del centro de masa se calculan a partir del triángulo rectángulo cuya hipotenusa es r y ángulo en P, θ+φ.

x=s-r·cos(θ+φ)
y=r·sin(θ+φ)

Para describir el movimiento del cilindro de base elíptica, es necesario expresar todas las magnitudes cinemáticas en función del ángulo girado θ

Ya tenemos todos los elementos necesarios para expresar en términos del ángulo girado θ las magnitudes cinemáticas que precisamos.

Ecuación del movimiento

La energía cinética es la suma del la energía cinética correspondiente al movimiento de traslación del centro de masas y la energía cinética de rotación alrededor del eje del cilindro.

E k = 1 2 m( ( dx dt ) 2 + ( dy dt ) 2 )+ 1 2 ( 1 4 m( a 2 + b 2 ) ) ( dθ dt ) 2 E k = 1 2 m( a 4 sin 2 θ+ b 4 cos 2 θ a 2 sin 2 θ+ b 2 cos 2 θ + 1 4 ( a 2 + b 2 ) ) ( dθ dt ) 2

La energía potencial es mgy, tomando el plano horizontal como nivel cero

E p =mgy=mg a 2 sin 2 θ+ b 2 cos 2 θ

Calculamos la lagrangiana y la ecuación del movimiento

L= E k E p = 1 2 m( a 4 sin 2 θ+ b 4 cos 2 θ a 2 sin 2 θ+ b 2 cos 2 θ + 1 4 ( a 2 + b 2 ) ) ( dθ dt ) 2 mg a 2 sin 2 θ+ b 2 cos 2 θ d dt ( L θ ˙ ) L θ =0 m( a 4 sin 2 θ+ b 4 cos 2 θ a 2 sin 2 θ+ b 2 cos 2 θ + 1 4 ( a 2 + b 2 ) ) d 2 θ d t 2 +2m sinθcosθ· a 2 b 2 ( a 2 b 2 ) ( a 2 sin 2 θ+ b 2 cos 2 θ ) 2 ( dθ dt ) 2 m sinθcosθ· a 2 b 2 ( a 2 b 2 ) ( a 2 sin 2 θ+ b 2 cos 2 θ ) 2 ( dθ dt ) 2 +mg sinθcosθ( a 2 b 2 ) a 2 sin 2 θ+ b 2 cos 2 θ =0

Simplificando, llegamos a la ecuación diferencial

( a 4 sin 2 θ+ b 4 cos 2 θ a 2 sin 2 θ+ b 2 cos 2 θ + 1 4 ( a 2 + b 2 ) ) d 2 θ d t 2 + sinθcosθ· a 2 b 2 ( a 2 b 2 ) ( a 2 sin 2 θ+ b 2 cos 2 θ ) 2 ( dθ dt ) 2 +g sinθcosθ( a 2 b 2 ) a 2 sin 2 θ+ b 2 cos 2 θ =0

Energías

La energía del cilindro que rueda sin deslizar es constante e igual a la inicial

E= E k + E p = 1 2 m( a 4 sin 2 θ+ b 4 cos 2 θ a 2 sin 2 θ+ b 2 cos 2 θ + 1 4 ( a 2 + b 2 ) ) ( dθ dt ) 2 +mg a 2 sin 2 θ+ b 2 cos 2 θ

Fuerzas

Las fuerzas que actúan sobre el cilindro son:

Aplicamos la segunda ley de Newton.

m d 2 x d t 2 = F r m d 2 y d t 2 =Nmg

Despejamos las fuerza por unidad de masa Fr/m y N/m

F r m = d dt ( a 2 sin 2 θ+ b 2 cos 2 θ dθ dt )= ( a 2 b 2 )sinθcosθ a 2 sin 2 θ+ b 2 cos 2 θ ( dθ dt ) 2 + a 2 sin 2 θ+ b 2 cos 2 θ d 2 θ d t 2 N m =g+ d 2 y d t 2 =g+ d 2 d t 2 ( a 2 sin 2 θ+ b 2 cos 2 θ )= g+ ( a 2 b 2 )sinθcosθ a 2 sin 2 θ+ b 2 cos 2 θ d 2 θ d t 2 +( a 2 b 2 ) b 2 cos 4 θ a 2 sin 4 θ ( a 2 sin 2 θ+ b 2 cos 2 θ ) 3/2 ( dθ dt ) 2

Resultados

Resolveremos la ecuación diferencial del movimiento con las siguientes condiciones iniciales: en el instante t=0, el ángulo girado es θ0 y la velocidad angular inicial ω0 =(dθ/dt)0

Representamos la energía potencial Ep/(mg) en función del ángulo θ para dos valores del semieje mayor a=1.5, 2, el semieje menor b=1

b=1;
th=110*pi/180; 
hold on
for a=[1.5,2];
    f=@(x) sqrt(a^2*sin(x).^2+b^2*cos(x).^2);
    fplot(f, [pi/2,3*pi/2])
    line([th,2*pi-th],[f(th),f(th)])
end
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'})
hold off
grid on
xlabel('\theta')
ylabel('E_p/mg')
title('Energía potencial')

Observamos que la energía potencial es máxima para θ=π/2 y θ=3π/2 (posiciones de equilibrio inestables) y mínima para θ=π (posición de equilibrio estable).

Cuando la velocidad angular inicial es nula, la energía del cilindro es menor que la energía potencial máxima, el cilindro describe oscilaciones de amplitud igual a π-θ0

Cuando la excentricidad de la elipse es grande, la variación de velocidad angular es grande y el procedimiento numérico ode45 de MATLAB puede que no funcione adecuadamente. Por lo que nos restringiremos a elipses de semiejes, a≈b

Oscilaciones

Movimiento de rodar sin deslizar sobre un plano horizontal

Situamos el cilindro en la posición de equilibrio inestable θ0=π/2 y le porporcionamos una velocidad angular inicial (dθ/dt)0=0.1 rad/s, de modo que su energía inicial es superior al máximo de energía potencial. El cilindro rodará a lo largo del plano horizontal

function elipse
    b=1;
    a=1.25;
    th_0=pi/2; 
    [t,x]=ode45(@eDiferencial,[0,10],[th_0,0.1]);
    e2=1-(b/a)^2;
    xPos=a*(ellipticE(x(:,1)+pi/2,e2)-ellipticE(th_0+pi/2,e2)); 
    plot(t,xPos)
    grid on
    xlabel('t')
    ylabel('x');
    title('Abscisa x');
    
    function dr=eDiferencial(~, x)
       p1=(a^4*sin(x(1))^2+b^4*cos(x(1))^2)/(a^2*sin(x(1))^2+
b^2*cos(x(1))^2)+(a^2+b^2)/4;
       p2=a^2*b^2*(a^2-b^2)*cos(x(1))*sin(x(1))/(a^2*sin(x(1))^2+
b^2*cos(x(1))^2)^2;
       p3=9.8*(a^2-b^2)*cos(x(1))*sin(x(1))/sqrt(a^2*sin(x(1))^2+
b^2*cos(x(1))^2);
       dr=[x(2);-(p2*x(2)^2+p3)/p1];
    end
end

Actividades

Se introduce

Se observa el movimiento del cilindro de base elíptica sobre la plano horizontal. En la parte superior derecha, se proporcionan los datos de

Utilizando los botones pausa || y paso >| se puede medir el periodo de oscilación, cuando el cilindro retorna a la posición de partida su velocidad angular de rotación es cero.

El programa calcula en cada instante el tanto por ciento de error relativo en la energía o el cociente

| E E 0 E 0 |·100

donde E es la energía del sistema en cualquier instante t, y E0 es la energía inicial del sistema.

Este valor se proporciona en caracteres de color rojo en la parte superior. Su valor debe ser siempre cero, o un valor muy pequeño lo que indica que la energía del sistema permanece constante y el programa realiza los cálculos correctamente.

Movimiento de rodar sin deslizar sobre un plano inlinado

El cilindro de base elíptica rueda sin deslizar a lo largo de un plano inclinado un ángulo α. Lo único que cambia es la energía potencial

El centro de masa se encuentra por debajo del nivel cero de energía potencial, una altura, xsinα-ycosα.

Calculamos la lagrangiana y la ecuación del movimiento. Solamente tenemos que modificar el término correspondiente a la energía potencial

L= E k E p = 1 2 m( a 4 sin 2 θ+ b 4 cos 2 θ a 2 sin 2 θ+ b 2 cos 2 θ + 1 4 ( a 2 + b 2 ) ) ( dθ dt ) 2 +mg( xsinαycosα ) d dt ( L θ ˙ ) L θ =0 m( a 4 sin 2 θ+ b 4 cos 2 θ a 2 sin 2 θ+ b 2 cos 2 θ + 1 4 ( a 2 + b 2 ) ) d 2 θ d t 2 +2m sinθcosθ· a 2 b 2 ( a 2 b 2 ) ( a 2 sin 2 θ+ b 2 cos 2 θ ) 2 ( dθ dt ) 2 m sinθcosθ· a 2 b 2 ( a 2 b 2 ) ( a 2 sin 2 θ+ b 2 cos 2 θ ) 2 ( dθ dt ) 2 mg( x θ sinα y θ cosα )=0

Simplificando, llegamos a la ecuación diferencial

( a 4 sin 2 θ+ b 4 cos 2 θ a 2 sin 2 θ+ b 2 cos 2 θ + 1 4 ( a 2 + b 2 ) ) d 2 θ d t 2 + sinθcosθ· a 2 b 2 ( a 2 b 2 ) ( a 2 sin 2 θ+ b 2 cos 2 θ ) 2 ( dθ dt ) 2 g( sinα a 2 sin 2 θ+ b 2 cos 2 θ sinθcosθ( a 2 b 2 ) a 2 sin 2 θ+ b 2 cos 2 θ cosα )=0

Cuando la pendiente α=0, obtenemos la ecuación del movimiento del cilindro sobre el plano horizontal

Resolvemos la ecuación diferencial por el procedimiento numérico ode45 con las siguientes condiciones inciales: en el instante t=0, la posición incial es θ0=π/2, (el eje mayor perpendicular al plano inclinado). Se suelta el cilindro en esta posición y baja rodando sin deslizar a lo largo del plano inclinado

Utilizaremos un cilindro cuya base es una elipse de semiejes a=2.5 y b=2 cm que rueda sobre un plano inclinado α=3.5°. Representamos la velocidad angular dθ/dt de rotación alrededor del eje del cilindro en función del tiempo

function elipse_1
    a=0.025;
    b=0.02;
    aRampa=3.5*pi/180;
    th_0=pi/2;   
    [t,x]=ode45(@eDiferencial,[0,1.5],[th_0,0]);
    yPos=sqrt(a^2*sin(x(:,1)).^2+b^2*cos(x(:,1)).^2);
    e2=1-(b/a)^2;
    xPos=a*(ellipticE(x(:,1)+pi/2,e2)-ellipticE(th_0+pi/2,e2)); 
    plot(t,x(:,2))
    grid on
    xlabel('t')
    ylabel('x');
    title('velocidad angular');
    m1=(a^4*sin(x(:,1)).^2+b^4*cos(x(:,1)).^2)./(a^2*sin(x(:,1)).^2+
b^2*cos(x(:,1)).^2)+(a^2+b^2)/4;
    Energia=m1.*x(:,2).^2/2-9.8*(xPos*sin(aRampa)-yPos*cos(aRampa));
    disp(Energia)
     
    function dr=eDiferencial(~, x)
       p1=(a^4*sin(x(1))^2+b^4*cos(x(1))^2)/(a^2*sin(x(1))^2+
b^2*cos(x(1))^2)+(a^2+b^2)/4;
       p2=a^2*b^2*(a^2-b^2)*cos(x(1))*sin(x(1))/(a^2*sin(x(1))^2+
b^2*cos(x(1))^2)^2;
       p3=9.8*(a^2-b^2)*cos(x(1))*sin(x(1))*cos(aRampa)/sqrt(a^2*sin(x(1))^2+
b^2*cos(x(1))^2)-9.8*sin(aRampa)*sqrt(a^2*sin(x(1))^2+b^2*cos(x(1))^2);
       dr=[x(2);-(p2*x(2)^2+p3)/p1];
    end
end

Comprobamos la conservación de la energía

    0.2445
    0.2445
....
    0.2472
    0.2490

Representamos la posición x del eje del cilindro en función del tiempo

Con una inclinación mayor, el procedimiento ode45 de MATLAB ya no responde adecuadamente, sin embargo, el procedimiento de Runge-Kutta utilizado en la simulación más abajo, sigue funcionando correctamente.

Actividades

Se introduce

Observamos el movimiento del cilindro rodando a lo largo del plano inclinado, en la parte superior derecha se proporcionan los datos de


Referencias

Johan Lindén, Kjell-Mikael Källman, Markus Lindberg. The rolling elliptical cylinder. Am. J. Phys. 89 (4), April 2021, pp. 358-364