Consideremos un aro de masa m_a y radio R, con una partícula de masa m unida en un punto del mismo situados sobre un plano horizontal. En la parte izquierda de la figura, se muestra la situación inicial, en el instante t=0. En la parte derecha de la figura, se muestra la situación en el instante t: La posición angular de la partícula es θ medida desde la vertical, sus coordenadas son x e y, mientras que la posición del centro del aro es X e Y, todos dependientes del tiempo t.

Las fuerza que actúan sobre el sistema son:

• El peso (m_a+m)g=Mg que actúa en el centro de masa (c.m.) del sistema
• La reacción del plano horizontal, N
• La fuerza de rozamiento F_r. Ambas fuerzas actúan en el punto P de contacto del aro y el plano horizontal

Situamos el origen O y los ejes X e Y como se muestra en la figura

Ecuaciones del movimiento

En este apartado, deducimos las ecuaciones del movimiento de traslación del c.m. y del movimientomde rotación alrededor de un eje perpendicular al plano del aro y que pasa por el c.m.

Movimiento de traslación del centro de masas

La distancia del centro de masas (c.m.) al centro C del aro es

${\displaystyle x_{cm}=\frac{m}{m_a+m}R=\frac{m}{M}R=\gamma R}$

Como vemos en la figura, la posición del c.m. del sistema (x,y) está relacionado con la posición del centro (X,Y) del aro del siguiente modo

${\displaystyle \overrightarrow{r}\{\begin{array}{l}x=X+\gamma R\sin \theta \\ y=Y+\gamma R\cos \theta \end{array}}$

Derivando una vez y dos veces con respecto del tiempo t, obtenemos el vector velocidad y el vector aceleración, respectivamente

${\displaystyle \begin{array}{l}\overrightarrow{v}\{\begin{array}{l}\frac{dx}{dt}=\frac{dX}{dt}+\gamma R\cos \theta \frac{d\theta }{dt}\\ \frac{dy}{dt}=\frac{dY}{dt}-\gamma R\sin \theta \frac{d\theta }{dt}\end{array}\\ \overrightarrow{a}\{\begin{array}{l}\frac{d^{2}x}{d{t}^2}=\frac{d^{2}X}{d{t}^2}+\gamma R\cos \theta \frac{d^2\theta }{d{t}^2}-\gamma R\sin \theta {\left(\frac{d\theta }{dt}\right)}^2\\ \frac{d^{2}y}{d{t}^2}=\frac{d^{2}Y}{d{t}^2}-\gamma R\sin \theta \frac{d^2\theta }{d{t}^2}-\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\end{array}\end{array}}$

Cuando el aro rueda o desliza está en contacto con el plano horizontal por lo que la coordenada Y(t)=R es constante y sus derivadas respecto del tiempo t, son nulas

Aplicamos la segunda ley de Newton al c.m. del sistema en la dirección del eje X y en la dirección del eje Y son, respectivamente

${\displaystyle \begin{array}{l}\{\begin{array}{l}M\frac{d^{2}x}{d{t}^2}=F_r\\ M\frac{d^{2}y}{d{t}^2}=N-Mg\end{array}\\ \{\begin{array}{l}{F}_r=M\frac{d^{2}X}{d{t}^2}+M\gamma R\cos \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\sin \theta {\left(\frac{d\theta }{dt}\right)}^2\\ N=Mg-M\gamma R\sin \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\end{array}\end{array}}$

Movimiento de rotación alrededor de un eje que pasa por el c.m.

El momento de inercia del sistema, respecto a un eje perpendicular al plano del aro y que pasa por su centro es

${\displaystyle I_C=m{R}^2+m_a{R}^2=\left(m+m_a\right)R^2}$

El momento de inercia del sistema respecto de un eje paralelo que pasa por el c.m. es

${\displaystyle I_{cm}=I_C-\left(m+m_a\right){\left(\gamma R\right)}^2=m{R}^2+m_a{R}^2=\left(1-{\gamma }^2\right)M{R}^2}$

La ecuación del movimiento de rotación alrededor de un eje que pasa por el c.m. es

${\displaystyle \begin{array}{l}{I}_{cm}\frac{d^2\theta }{d{t}^2}=N\gamma R\sin \theta -F_r\left(R+\gamma R\cos \theta \right)\\ \left(1-{\gamma }^2\right)MR\frac{d^2\theta }{d{t}^2}=N\gamma \sin \theta -F_r\left(1+\gamma \cos \theta \right)\end{array}}$

La energía del sistema

La energía potencial del c.m. referida al plano horizontal es E_p=Mg(R+γRcosθ)

La energía cinética es la suma de la energía cinética de traslación del centro de masa y de rotación alrededor de un eje que pasa por el c.m.

${\displaystyle \begin{array}{l}{E}_k=\frac{1}{2}{I}_{cm}{\left(\frac{d\theta }{dt}\right)}^2+\frac{1}{2}M\left({\left(\frac{dx}{dt}\right)}^2+{\left(\frac{dy}{dt}\right)}^2\right)=\\ \frac{1}{2}\left(1-{\gamma }^2\right)M{R}^2{\left(\frac{d\theta }{dt}\right)}^2+\frac{1}{2}M\left({\left(\frac{dX}{dt}+\gamma R\cos \theta \frac{d\theta }{dt}\right)}^2+{\left(-\gamma R\sin \theta \frac{d\theta }{dt}\right)}^2\right)=\\ \frac{1}{2}M{R}^2{\left(\frac{d\theta }{dt}\right)}^2+\frac{1}{2}M\left({\left(\frac{dX}{dt}\right)}^2+2\gamma R\cos \theta \frac{dX}{dt}\frac{d\theta }{dt}\right)\end{array}}$

La energía total E_p+E_k es constante e igual a

${\displaystyle E=MgR\left(1+\gamma \cos \theta \right)+\frac{1}{2}M{R}^2{\left(\frac{d\theta }{dt}\right)}^2+\frac{1}{2}M\left({\left(\frac{dX}{dt}\right)}^2+2\gamma R\cos \theta \frac{dX}{dt}\frac{d\theta }{dt}\right)}$

Etapas del movimiento

Las posibles etapas del movimiento del cuerpo formado por la partícula unida al aro, son:

• Rueda sin deslizar
• Desliza, con dos posibilidades: la fuerza de rozamiento F_r>0 y F_r<0
• Salta, cuando la reacción del plano horizontal, N=0 o se detiene, cuando la velocidad es cero

Rueda sin deslizar

Cuando el cuerpo rueda sin deslizar, hay una relación entre la velocidad del centro del aro dX/dt y la velocidad angular de rotación de la partícula, dθ/dt

${\displaystyle \frac{dX}{dt}=R\frac{d\theta }{dt}}$

Cuando se cumple esta condición, la fuerza de rozamiento F_r y la reacción del plano horizontal N en el punto de contacto valen

${\displaystyle \{\begin{array}{l}{F}_r=MR\left(1+\gamma \cos \theta \right)\frac{d^2\theta }{d{t}^2}-M\gamma R\sin \theta {\left(\frac{d\theta }{dt}\right)}^2\\ N=Mg-M\gamma R\sin \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\end{array}}$

Introduciendo estas expresiones en la ecuación del movimiento de rotación, llegamos a la ecuación difrencial que nos proporciona, la posición angular θ de la partícula en función del tiempo t

${\displaystyle \frac{d^2\theta }{d{t}^2}=\frac{\gamma \sin \theta }{2R\left(1+\gamma \cos \theta \right)}\left(g+R{\left(\frac{d\theta }{dt}\right)}^2\right)}$

La expresión de la fuerza de rozamiento se simplifica

${\displaystyle F_r=M\frac{\gamma \sin \theta }{2}\left(g-R{\left(\frac{d\theta }{dt}\right)}^2\right)}$

Conservación de la energía

Cuando el aro rueda sin deslizar, la energía total es constante e igual a la energía inicial

${\displaystyle E=MR\left(1+\gamma \cos \theta \right)\left(g+R{\left(\frac{d\theta }{dt}\right)}^2\right)}$

Desliza

El aro desliza cuando la fuerza de rozamiento F_r alcanza el valor máximo (en valor absoluto), cuando se cumple que |F_r|=μ_sN, donde μ_s es el coeficiente estático. Los dos posibles casos son

• La velocidad angular de rotación es mayor que la de traslación, Rdθ/dt>dX/dt

La velocidad del punto P de contacto con el plano horizontal es v_P<0 por lo que F_r>0, tal como se muestra en la figura de la izquierda. En las ecuaciones del movimiento de traslación del c.m. ponemos F_r=μ_kN, siendo μ_k el coeficiente cinético

${\displaystyle \{\begin{array}{l}{\mu }_{k}N=M\frac{d^{2}X}{d{t}^2}+M\gamma R\cos \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\sin \theta {\left(\frac{d\theta }{dt}\right)}^2\\ N=Mg-M\gamma R\sin \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\end{array}}$

y hacemos lo mismo, en la ecuación del movimiento de rotación alrededor de un eje que pasa por el c.m.

${\displaystyle \begin{array}{l}\left(1-{\gamma }^2\right)MR\frac{d^2\theta }{d{t}^2}=N\gamma \sin \theta -{\mu }_{k}N\left(1+\gamma \cos \theta \right)\\ \left(1-{\gamma }^2\right)MR\frac{d^2\theta }{d{t}^2}=N{S}_{+}(\theta )\\ S_{+}(\theta )=\gamma \sin \theta -{\mu }_k\left(1+\gamma \cos \theta \right)\end{array}}$

La ecuación difrencial que nos proporciona, la posición angular θ de la partícula en función del tiempo t es

${\displaystyle \frac{d^2\theta }{d{t}^2}=\frac{\left(g-\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\right)S_{+}(\theta )}{R\left(1-{\gamma }^2+\gamma \sin \theta ·S_{+}(\theta )\right)}}$

La aceleración del centro del aro es

${\displaystyle \begin{array}{l}{\mu }_{k}N=M\frac{d^{2}X}{d{t}^2}+M\gamma R\cos \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\sin \theta {\left(\frac{d\theta }{dt}\right)}^2\\ \frac{d^{2}X}{d{t}^2}=\left\{\begin{array}{l}{\mu }_{k}g-\gamma R\left({\mu }_k\sin \theta +\cos \theta \right)\frac{\left(g-\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\right)S_{+}(\theta )}{R\left(1-{\gamma }^2+\gamma \sin \theta ·S_{+}(\theta )\right)}\\ +\gamma R\left(\sin \theta -{\mu }_k\cos \theta \right){\left(\frac{d\theta }{dt}\right)}^2\end{array}\right\}\end{array}}$

• La velocidad angular de rotación es menor que la de traslación, Rdθ/dt<dX/dt

La velocidad del punto P de contacto con el plano horizontal es v_P>0 por lo que F_r<0, tal como se muestra en la figura de la derecha. En las ecuaciones del movimiento de traslación del c.m. ponemos F_r=-μ_kN

${\displaystyle \{\begin{array}{l}-{\mu }_{k}N=M\frac{d^{2}X}{d{t}^2}+M\gamma R\cos \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\sin \theta {\left(\frac{d\theta }{dt}\right)}^2\\ N=Mg-M\gamma R\sin \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\end{array}}$

y hacemos lo mismo, en la ecuación del movimiento de rotación alrededor de un eje que pasa por el c.m.

${\displaystyle \begin{array}{l}\left(1-{\gamma }^2\right)MR\frac{d^2\theta }{d{t}^2}=N\gamma \sin \theta +{\mu }_{k}N\left(1+\gamma \cos \theta \right)\\ \left(1-{\gamma }^2\right)MR\frac{d^2\theta }{d{t}^2}=N{S}_{-}(\theta )\\ S_{-}(\theta )=\gamma \sin \theta +{\mu }_k\left(1+\gamma \cos \theta \right)\end{array}}$

La ecuación difrencial que nos proporciona, la posición angular θ de la partícula en función del tiempo t es

${\displaystyle \frac{d^2\theta }{d{t}^2}=\frac{\left(g-\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\right)S_{-}(\theta )}{R\left(1-{\gamma }^2+\gamma \sin \theta ·S_{-}(\theta )\right)}}$

La aceleración del centro del aro es

${\displaystyle \begin{array}{l}-{\mu }_{k}N=M\frac{d^{2}X}{d{t}^2}+M\gamma R\cos \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\sin \theta {\left(\frac{d\theta }{dt}\right)}^2\\ \frac{d^{2}X}{d{t}^2}=\left\{\begin{array}{l}-{\mu }_{k}g-\gamma R\left(-{\mu }_k\sin \theta +\cos \theta \right)\frac{\left(g-\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\right)S_{-}(\theta )}{R\left(1-{\gamma }^2+\gamma \sin \theta ·S_{-}(\theta )\right)}\\ +\gamma R\left(\sin \theta +{\mu }_k\cos \theta \right){\left(\frac{d\theta }{dt}\right)}^2\end{array}\right\}\end{array}}$

Obtenemos ecuaciones similares en ambos casos, solamente tenemos que cambiar μ_k por -μ_k

Cuando desliza la energía no se conserva ya que hay un trabajo de la fuerza de rozamiento

Salta

El aro salta, cuando la reacción N se hace cero en el instante t_j. A partir de este instante el aro deja de estar en contacto con el plano horizontal. La aceleración vertical del centro de masas del sistema y la aceleración del centro del aro son, respectivamente

${\displaystyle \begin{array}{l}M\frac{d^{2}y}{d{t}^2}=-Mg\\ \frac{d^{2}Y}{d{t}^2}=-g+\gamma R\sin \theta \frac{d^2\theta }{d{t}^2}+\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\end{array}}$

La ecuación de la dinámica de rotación alrededor del un eje que pasa por el c.m. nos dice que la aceleración angular es nula, la velocidad angular dθ/dt es constante

La aceleración del centro del aro para t>t_j es

${\displaystyle \frac{d^{2}Y}{d{t}^2}=\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2-g}$

Esta aceleración tiene que ser positiva, se tiene que cumplir en el instante t_j

${\displaystyle \gamma R\cos {\theta }_j{\left(\frac{d\theta }{dt}\right)}_j^2>g}$

Rueda, desliza, rueda, desliza y salta

Vamos a resolver un ejemplo paso a paso, para entender la complejidad de los posibles movimientos de un sistema tan simple que depende de tres parámetros (fijado el radio R del aro):

• γ, que depende del cociente entre la masa del aro y de la partícula, establece la distancia γR del centro de masa del sistema al centro del aro
• El coeficiente de rozamiento μ=μ_s=μ_k
• La velocidad angular inicial de rotación, (dθ/dt)₀

El radio del aro se ha fijado en R=1. Los valores de los parámetros son: γ=2/3, μ=0.6, el cuadrado de la velocidad angular adimensional $\frac{R}{g}{\left(\frac{d\theta }{dt}\right)}_0^2=1.1$

Situación inicial

La posición angular de la partícula en el instante t=0, es θ₀=0 y X₀=0, tal como se indica en la primera figura de esta página. Supongamos que proporcionamos al sistema una velocidad angular inicial de rotación (dθ/dt)₀. La reacción N del plano horizontal es

${\displaystyle N(0)=Mg-\gamma RM{\left(\frac{d\theta }{dt}\right)}_0^2}$

que tiene que ser positiva, lo que limita la velocidad angular inicial. Para los valores elegidos de los parámetros se cumple esta condición, g-2·1.1·g/3=0.2667·g>0

Rueda sin deslizar

Para que el aro ruede inicialmente sin deslizar, el centro del aro deberá llevar en el instante t=0, una velocidad inicial (dX/dt)₀=R(dθ/dt)₀

Para esta situación inicial θ₀=0, la fuerza de rozamiento es nula, F_r=0

El principio de conservación de la energía nos permite calcular la velocidad angular dθ/dt en función de la posición angular θ de la partícula

${\displaystyle {\left(\frac{d\theta }{dt}\right)}^2=\frac{\left(1+\gamma \cos {\theta }_0\right)}{\left(1+\gamma \cos \theta \right)}\left(\frac{g}{R}+{\left(\frac{d\theta }{dt}\right)}_0^2\right)-\frac{g}{R}}$

Conocida la expresión de la velocidad angular en función de la posición angular θ, obtenemos la aceleración angular

${\displaystyle \frac{d^2\theta }{d{t}^2}=\frac{\gamma \sin \theta }{2R\left(1+\gamma \cos \theta \right)}\left(g+R{\left(\frac{d\theta }{dt}\right)}^2\right)}$

Conocida la velocidad angular y la aceleración angular calculamos la fuerza de rozamiento F_r y la normal N

${\displaystyle \{\begin{array}{l}{F}_r=M\frac{\gamma \sin \theta }{2}\left(g-R{\left(\frac{d\theta }{dt}\right)}^2\right)\\ N=Mg-M\gamma R\sin \theta \frac{d^2\theta }{d{t}^2}-M\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\end{array}}$

La fuerza de rozamiento es positiva para θ<180° siempre que g>R(dθ/dt)² y negativa en caso contrario. Con el dato del cuadrado de la velocidad angular inicial de este ejemplo, g<R·1.1·g/R, la fuerza de rozamiento es negativa, F_r<0.

Llamaremos θ_s al ángulo para el cual la fuerza de rozamiento (en valor absoluto) |F_r| se hace igual a μN. A partir de dicha posición el aro desliza

Creamos un script de MATLAB para representar el cociente F_r/N en función de la posición angular θ y su intersección con la recta ±μ dependiendo si la fuerza de rozamiento F_r es positiva o negativa. La abscisa de dicho punto es el ángulo θ_s, que calcularemos con más precisión utilizando la función fzero de MATLAB

gamma=2/3; %distancia al c.m.
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0; %posición angular inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
disp(aDesliza*180/pi)
hold on
fplot(f,[0,2*pi])
line([0,2*pi],[-mu,-mu])
plot(aDesliza,f(aDesliza),'ro','markersize',4,'markeredgecolor',
'b','markerfacecolor','b')
hold off
set(gca,'XTick',0:pi/4:2*pi)
set(gca,'XTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
grid on
xlabel('\theta')
ylabel('f')
title('Cociente, Fr/N')

Como la partícula parte de la posicón angular θ=0, la aceleración angular es nula y la fuerza de rozamiento F_r es nula. La fuerza de rozamiento decrece, hasta que en valor absoluto se hace igual a μ_sN. A partir de la posición θ_s el aro empieza a deslizar

El aro empieza a deslizar en la posición θ_s=46.9°. Utilizamos el principio de conservación de la energía para calcular la velocidad angular final, (dθ/dt)_s

  46.9250
>> sqrt(w2(aDesliza))
ans =    3.7107

Para determinar la posición angular de la partícula θ en función del tiempo t se utiliza el procedimiento ode45 de MATLAB que resuelve la ecuación diferencial con las condiciones iniciales siguientes: en el instante t=0, θ₀=0, (dθ/dt)=(dθ/dt)₀. Como el aro rueda sin deslizar, la posición del centro del aro es X=R·θ, Y=R

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0; %posición angular inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);

x0=[0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
plot(t,x(:,1))
set(gca,'YTick',0:pi/4:2*pi)
set(gca,'YTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
grid on
xlabel('t')
ylabel('\theta')
title('Rueda sin deslizar')

El procedimiento ode45 se detiene cuando la posición angular de la partícula alcanza el ángulo θ_s=46.9°. Para ello definimos las función stop_aro_particula

function [value,isterminal,direction]=stop_aro_particula(t,x, ang)
    value=x(1)-ang; 
    isterminal=1;
    direction=0; 
end

Se alcanza la posición angular final θ_s en el instante y con velocidad

>> t(end) %instante
ans =    0.2395
>> x(end,1)*180/pi %posición
ans =   46.9250
>> x(end,2) %velocidad angular
ans =    3.7107

Confirmamos que en este instante t=0.2395 o para esta posición angular θ_s=46.9250°, se cumple que |F_r|=μN. Primero calculamos la aceleración angular

>> acel=gamma*sin(x(end,1))*(9.8+R*x(end,2)^2)
/(2*R*(1+gamma*cos(x(end,1))))
acel =    3.9433
>> Fr=R*(1+gamma*cos(x(end,1)))*acel-gamma*R*sin(x(end,1))*x(end,2)^2
Fr =   -0.9664
>> -mu*(9.8-gamma*R*sin(x(end,1))*acel-gamma*R*cos(x(end,1))*x(end,2)^2)
ans =   -0.9664

La fuerza de rozamiento F_r<0, se iguala (en valor absoluto) a μN para el ángulo θ_s=46.9°, el cuerpo formado por el aro y la partícula empieza a deslizar.

Desliza

Estamos en el segundo caso, F_r<0,

Las ecuaciones diferenciales que describen la posición angular θ de la partícula en función del tiempo t y la que describe la posición X del centro del aro en función del tiempo son

${\displaystyle \begin{array}{l}\frac{d^2\theta }{d{t}^2}=\frac{\left(g-\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\right)S_{-}(\theta )}{R\left(1-{\gamma }^2+\gamma \sin \theta ·S_{-}(\theta )\right)}\\ \frac{d^{2}X}{d{t}^2}=\left\{\begin{array}{l}-{\mu }_{k}g-\gamma R\left(-{\mu }_k\sin \theta +\cos \theta \right)\frac{\left(g-\gamma R\cos \theta {\left(\frac{d\theta }{dt}\right)}^2\right)S_{-}(\theta )}{R\left(1-{\gamma }^2+\gamma \sin \theta ·S_{-}(\theta )\right)}\\ +\gamma R\left(\sin \theta +{\mu }_k\cos \theta \right){\left(\frac{d\theta }{dt}\right)}^2\end{array}\right\}\\ S_{-}(\theta )=\gamma \sin \theta +{\mu }_k\left(1+\gamma \cos \theta \right)\end{array}}$

El movimiento de deslizar se termina, cuando se cumpla que dX/dt=Rdθ/dt, la condición de rodar sin deslizar.

Resolvemos el sistema de dos ecuaciones diferenciales, con las condiciones iniciales siguientes: en el instante t=0, θ₀=θ_s, (dθ/dt)₀=(dθ/dt)_s, X₀=R·θ_s, (dX/dt)₀=R(dθ/dt)_s

Representamos la velocidad angular (dθ/dt) de la partícula y la velocidad (dX/dt) del centro del aro en función del tiempo t.

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0; %posición angular inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);

w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
%ecuaciones para Fr<0
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))
+mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))
*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))
/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))
+gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, R, sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
hold on
plot(t,x(:,2)) 
plot(t,x(:,4)) 
hold off
%set(gca,'YTick',0:pi/4:2*pi)
%set(gca,'YTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
legend('d\theta/dt', 'dX/dt', 'location','northwest')
grid on
xlabel('t')
ylabel('d\theta/dt, dX/dt')
title('Desliza')

En el instante inicial t=0, se cumple la condición de rodar sin deslizar, ambas velocidades coinciden ya que el radio del aro R=1. La velocidad de traslación del centro del aro dX/dt se hace mayor que la velocidad angular de rotación (R·dθ/dt), la velocidad del punto de contacto P es positiva y la fuerza de rozamiento F_r es negativa

La condición de rodar sin deslizar (dX/dt)_r=R(dθ/dt)_r se restablece para un instante t_r y para una posición angular θ_r

Podría ocurrir que no se alcanzase la posición θ_r por que el aro se detuviese antes. Ahora bien, el hecho de que dθ/dt=0 no implicaría que dX/dt=0

El procedimiento ode45 de integración se detiene cuando se alcanza la condición de rodar sin deslizar, definiendo la función stop_aro_particula1

function [value,isterminal,direction]=stop_aro_particula1(t,x, R, signo)
    value=-signo;  %signo de la fuerza de rozamiento 
     if t>0.1
       value=x(4)-R*x(2); 
    end
    isterminal=1;
    direction=0; 
end

Si no ponemos la sentencia condicional, el procedimiento ode45 se detiene en el instante inicial, ya que se cumple la condición de rodar sin deslizar y por tanto, no progresaría. El valor t>0.1 es arbitrario

La posición angular final es θ_r=2.7919 (160°), comprobamos que se cumple la condición de rodar sin deslizar, tal como se muestra en la gráfica

>> t(end) %tiempo
ans =    0.3756
>> x(end,1) %posición angular
ans =    2.7919
>> x(end,3) %posición del centro del aro
ans =    2.9330
>>R*x(end,2) %velocidad angular
ans =    8.9086
>> x(end,4) %velocidad del centro del aro
ans =    8.9086

Rueda sin deslizar

El cuerpo formado por el aro y la partícula rueda sin deslizar a lo largo del plano horizontal, utilizamos el código de la primera parte de esta sección para calcular el ángulo θ_s que empieza a deslizar, partiendo de la posición angular inicial θ₀=θ_r y de la velocidad angular inicial (dθ/dt)₀=(dθ/dt)_r

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
disp(aDesliza*180/pi)

%desliza (con las ecuaciones de Fr<0)
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))
+mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))
*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, R, sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);
disp(aDesliza*180/pi)
hold on
fplot(f,[x(end,1),3*pi/2])
line([x(end,1),3*pi/2],[mu,mu])
plot(aDesliza,f(aDesliza),'ro','markersize',4,'markeredgecolor','b'
,'markerfacecolor','b')
hold off
set(gca,'XTick',0:pi/4:2*pi)
set(gca,'XTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
grid on
xlabel('\theta')
ylabel('f')
title('Cociente, Fr/N')

La fuerza de rozamiento F_r ha cambiado de negativa a positiva. Buscamos el ángulo θ_s para el cual la fuerza de rozamiento alcanza su valor máximo μ_sN

El aro empieza a deslizar en la posición θ_s=231.3°. Utilizamos el principio de conservación de la energía para calcular la velocidad angular final, (dθ/dt)_s

   231.2935
>> sqrt(w2(aDesliza))
ans =    6.8804

Para determinar la posición angular de la partícula θ en función del tiempo t, resolvemos la ecuación diferencial con las condiciones iniciales siguientes: en el instante t=0, θ₀=θ_r, (dθ/dt)₀=(dθ/dt)_r. Como el aro rueda sin deslizar, la posición del centro del aro es X=X₀+R·θ, Y=R. Donde X₀ es la posición inicial del centro del aro al terminar la etapa anterior (desliza)

Desliza

Estamos en el primer caso, F_r>0. Las ecuaciones para este caso se obtienen del segundo, simplemente, cambiando μ_k por -μ_k.

Resolvemos el sistema de dos ecuaciones diferenciales, con las condiciones iniciales siguientes: en el instante t=0, θ₀=θ_s, (dθ/dt)₀=(dθ/dt)_s, X₀=R·θ_s, (dX/dt)₀=(R·dθ/dt)_s

Representamos la velocidad angular (dθ/dt) de la partícula y la velocidad (dX/dt) del centro del aro.

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
disp(aDesliza*180/pi)

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))
*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))
/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))
+gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, R, sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0]; 
%desliza (las ecuaciones de Fr<0 con mu<0; se convierten en Fr>0 )
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))
+gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, R,sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
hold on
plot(t,R*x(:,2)) 
plot(t,R*x(:,4)) 
hold off
legend('d\theta/dt', 'dX/dt', 'location','northwest')
grid on
xlabel('t')
ylabel('d\theta/dt,dX/dt')
title('Desliza, velocidades')

En el instante inicial t=0, se cumple la condición de rodar sin deslizar, ambas velocidades coinciden. La velocidad de traslación del centro del aro dX/dt se hace menor que la velocidad angular de rotación (R·dθ/dt), la velocidad del punto de contacto P es negativa y la fuerza de rozamiento F_r es positiva

La condición de rodar sin deslizar (dX/dt)_r=R(dθ/dt)_r se restablece para un instante t_r y para una posición angular θ_r

El procedimiento ode45 de integración se detiene cuando se alcanza la condición de rodar sin deslizar, definiedo la función stop_aro_particula1

>> x(1,end)
ans =    6.8804

La posición angular final es θ_r=6.88 (394°). Sin embargo, el sistema no alcanza esta posición angular ya que la reacción N se ha hecho cero mientras deslizaba

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
disp(aDesliza*180/pi)

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))
*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, R,sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))
+gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, R,sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
N=9.8-(gamma*R*sin(x(:,1)).*(9.8-gamma*R*cos(x(:,1)).*x(:,2).^2).*
(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1)))))./(R*(1-gamma^2+gamma*sin(x(:,1)).
*(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))))-gamma*R*cos(x(:,1)).*x(:,2).^2;
plot(x(:,1),N)

grid on
set(gca,'XTick',0:pi/4:2*pi)
set(gca,'XTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
xlabel('\theta')
ylabel('N')
title('Reacción, N')

Salta

La reacción N se hace cero aproximadamente, para el ángulo θ=5.58 (319.7°). El centro de masa del aro describe una trayectoria parabólica

Para que el proceso de integración mediante ode45 se detenga cuando la reacción N se haga cero, redefinimos la función stop_aro_particula1, para incluir esta posiblidad

function [value,isterminal,direction]=stop_aro_particula1
(t,x, gamma, mu, R, signo)
    N=9.8-(gamma*R*sin(x(1))*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))))-gamma*R*cos(x(1))*x(2)^2;
    value=-signo;  
     if t>0.1
       value=[x(4)-R*x(2), N]; 
    end
    isterminal=[1, 1];
    direction=[0,-1]; 
end

Representamos la posición angular θ de la partícula y la posición X del centro del aro, en función del tiempo t

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
disp(aDesliza*180/pi)

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))
+mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))
*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R, 
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))
+mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+
cos(x(1)))*(9.8-gamma*R*cos(x(1))*x(2)^2)*
(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*
(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))))+gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R,
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
hold on
plot(t,x(:,1)) 
plot(t,x(:,3)) 
hold off
legend('\theta', 'X', 'location','northwest')
grid on
xlabel('t')
ylabel('\theta,X')
title('Desliza')

La reacción N se hace cero para la posición angular θ_j=5.6274 (322.4°). Mostramos los datos de la velocidad angular final, el desplazamiento del centro del aro en este tiempo y la velocidad final del centro del aro

>> t(end)
ans =    0.3167
>> x(end,1) %posición angular de la partícula
ans =    5.6274
>> x(end,2) %velocidad angular de la partícula
ans =    4.3067
>> x(end,3) %desplazamiento del centro del aro en este tiempo
ans =    5.3756
>> x(end,4) %velocidad del centro del aro
ans =    2.5016

Comprobamos que en el momento del salto, t_j=0.3167

${\displaystyle \gamma R\cos {\theta }_j{\left(\frac{d\theta }{dt}\right)}_j^2>g}$

>> gamma*R*x(end,2)^2
ans =   12.3652

Descripción completa

Hemos descrito las distintas etapas del movimiento del sistema formado por el aro y la partícula en términos de la posición angular de la partícula, θ.

• Rueda sin deslizar desde 0 a 46.9°
• Desliza desde 46.9° a 160°
• Rueda sin deslizar 160° a 231.3°
• Desliza desde 231.3° a 319.7°
• Salta en la posición angular 319.7°

Reunimos las porciones de código que se han descrito en la sección anterior, para representar:

1. La posición angular θ de la partícula y la posición X del centro del aro en función del tiempo t
2. La velocidad angular R·dθ/dt de la partícula y la velocidad dX/dt del centro del aro en función la posición angular θ
3. La aceleración angular de la partícula en función de la posición angular θ
4. La fuerza de rozamiento F_r en función de la posición angular θ
5. La reacción N en función de la posición angular θ
6. La energía total del aro en función del tiempo t

Posiciones

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial

%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
x0=[0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
hold on
plot(t,x(:,1),'r')
line([t(end),t(end)],[0,x(end,1)],'lineStyle','--','color','k')
line([0,t(end)],[x(end,1),x(end,1)],'lineStyle','--','color','k')
tFin=t(end);

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))
*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))
/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))
+gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R, 
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
plot(t+tFin,x(:,1),'r')
plot(t+tFin,x(:,3),'b')
line([tFin+t(end),tFin+t(end)],[0,x(end,1)],'lineStyle','--','color','k')
line([0,tFin+t(end)],[x(end,1),x(end,1)],'lineStyle','--','color','k')
tFin=tFin+t(end);
xFin=x(end,3)-R*x(end,1);

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);
x0=[th_0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
plot(t+tFin,x(:,1),'r')
plot(t+tFin,R*x(:,1)+xFin,'b')
line([tFin+t(end),tFin+t(end)],[0,x(end,1)],'lineStyle','--','color','k')
line([0,tFin+t(end)],[x(end,1),x(end,1)],'lineStyle','--','color','k')
tFin=tFin+t(end);

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))
*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))
/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))
+gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R,
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
plot(t+tFin,x(:,1),'r')
plot(t+tFin,x(:,3)+xFin,'b')
line([tFin+t(end),tFin+t(end)],[0,x(end,1)],'lineStyle','--','color','k')
line([0,tFin+t(end)],[x(end,1),x(end,1)],'lineStyle','--','color','k')

hold off
set(gca,'YTick',0:pi/4:2*pi)
set(gca,'YTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
grid on
xlabel('t')
ylabel('\theta,X')
title('Aro-partícula que ruedan')

En color rojo, la posición angular de la partícula. En azul, la posición del centro del aro en función del tiempo t. Cuando el cuerpo rueda sin deslizar X=R·θ, como el radio del aro R=1, la posición angular θ de la partícula y la del centro X del aro coinciden.

Velocidades

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial

%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
x0=[0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
hold on
plot(x(:,1),x(:,2),'r')
line([x(end,1),x(end,1)],[0,x(end,2)],'lineStyle','--','color','k')
line([0,x(end,1)],[x(end,2),x(end,2)],'lineStyle','--','color','k')

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))
+gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R, 
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
plot(x(:,1),x(:,2),'r')
plot(x(:,1),x(:,4),'b')
line([x(end,1),x(end,1)],[0,x(end,2)],'lineStyle','--','color','k')
line([0,x(end,1)],[x(end,2),x(end,2)],'lineStyle','--','color','k')

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);
x0=[th_0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
plot(x(:,1),x(:,2),'r')
line([x(end,1),x(end,1)],[0,x(end,2)],'lineStyle','--','color','k')
line([0,x(end,1)],[x(end,2),x(end,2)],'lineStyle','--','color','k')

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))
+mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R,
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
plot(x(:,1),x(:,2),'r')
plot(x(:,1),x(:,4),'b')
line([x(end,1),x(end,1)],[0,x(end,2)],'lineStyle','--','color','k')
line([0,x(end,1)],[x(end,2),x(end,2)],'lineStyle','--','color','k')
hold off
grid on
set(gca,'XTick',0:pi/4:2*pi)
set(gca,'XTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
xlabel('\theta')
ylabel('d\theta/dt,dX/dt')
title('Aro-partícula que ruedan')

En color rojo, la velocidad angular de la partícula. En azul, la velocidad del centro del aro en función de la posición angular θ de la partícula. Ambas coinciden cuando el cuerpo rueda sin deslizar

Cuando rueda sin deslizar la energía E se conserva

${\displaystyle E=MR\left(1+\gamma \cos \theta \right)\left(g+R{\left(\frac{d\theta }{dt}\right)}^2\right)}$

El máximo de la velocidad angular se producirá para el ángulo θ=π (180°)

En el intervalo 46.9° - 160° el aro desliza, la velocidad del centro del aro dX/dt es mayor que la velocidad angular de rotación Rdθ/dt. La velocidad del punto P de contacto del aro con el plano horizontal es positiva, la fuerza de rozamiento F_r es negativa. En el intervalo 231.3° - 319.7° ocurre lo contrario

Aceleración

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial

%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
hold on
fplot(acel,[0,aDesliza],'r')
line([aDesliza,aDesliza],[0,N(aDesliza)],'lineStyle','--','color','k')

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R, 
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
acel=(9.8-gamma*R*cos(x(:,1)).*x(:,2).^2).*(gamma*sin(x(:,1))+
mu*(1+gamma*cos(x(:,1))))./(R*(1-gamma^2+gamma*sin(x(:,1)).*
(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))));
acel_X=-mu*9.8-gamma*R*(-mu*sin(x(:,1))+cos(x(:,1))).*
acel+gamma*R*(sin(x(:,1))+mu*cos(x(:,1))).*x(:,2).^2;
plot(x(:,1),acel,'r')
plot(x(:,1),acel_X,'b')
line([x(end,1),x(end,1)],[0,acel(end)],'lineStyle','--','color','k')

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);
fplot(acel,[x(end,1),aDesliza],'r')
line([aDesliza,aDesliza],[0,acel(aDesliza)],'lineStyle','--','color','k')

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))
+gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R,
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
acel=(9.8-gamma*R*cos(x(:,1)).*x(:,2).^2).*(gamma*sin(x(:,1))+mu*
(1+gamma*cos(x(:,1))))./(R*(1-gamma^2+gamma*sin(x(:,1)).*(gamma*sin(x(:,1))+
mu*(1+gamma*cos(x(:,1))))));
acel_X=-mu*9.8-gamma*R*(-mu*sin(x(:,1))+cos(x(:,1))).*acel+gamma*R*(sin(x(:,1))
+mu*cos(x(:,1))).*x(:,2).^2;
plot(x(:,1),acel,'r')
plot(x(:,1),acel_X,'b')
line([x(end,1),x(end,1)],[0,acel(end)],'lineStyle','--','color','k')

hold off
grid on
set(gca,'XTick',0:pi/4:2*pi)
set(gca,'XTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
xlabel('\theta')
ylabel('d^2\theta/dt^2,d^2X/dt^2 ')
title('Aro-partícula que ruedan')

En color rojo, la aceleración angular de la partícula. En azul, la aceleración del centro del aro en función de la posición angular θ de la partícula. Ambas coinciden cuendo el aro rueda sin deslizar

La primera característica que observamos es la discontinuidad de ambas funciones para θ=160°, es decir, cuando comienza a rodar después de deslizar

Consideremos en el intervalo, 46.9°-160°, el aro desliza: la velocidad del centro del aro dX/dt es al principio igual que la velocidad angular de la partícula Rdθ/dt, la primera se va haciendo más grande que la segunda, lo que implica que la aceleración del centro del aro tiene que ser mayor que la aceleración angular de la partícula. Al final del intervalo, ambas velocidades vuelven a coincidir, lo que implica que la aceleración del centro del aro tiene que ser más pequeña que la angular de la partícula

Consideremos en el intervalo, 231° - 322°, el aro desliza: la velocidad del centro del aro dX/dt es al principio igual que la velocidad angular de la partícula Rdθ/dt, luego, la velocidad del centro del aro, diminuye más rápidamente que la angular de la partícula, lo que implica que la aceleración del centro del aro tendrá que ser más negativa que la angular de la partícula, tal como se muestra en la figura

La aceleración angular presenta también un máximo acusado para θ=160°, de algo más de cuatro veces la aceleración de la gravedad, se hace cero para θ=180°

Energía

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial

%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
x0=[0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
hold on
E=R*(1+gamma*cos(x(:,1))).*(9.8+R*x(:,2).^2);
plot(t,E)
line([t(end),t(end)],[33,E(end)],'lineStyle','--','color','k')
tFin=t(end);

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))
*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R, s
ign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
E=9.8*R*(1+gamma*cos(x(:,1)))+R^2*x(:,2).^2/2+(x(:,4).^2+
2*gamma*R*(cos(x(:,1)).*x(:,4)).*x(:,2))/2;
plot(t+tFin,E)
line([tFin+t(end),tFin+t(end)],[33,E(end)],'lineStyle','--','color','k')
tFin=tFin+t(end);


%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);
x0=[th_0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
E=R*(1+gamma*cos(x(:,1))).*(9.8+R*x(:,2).^2);
plot(t+tFin,E)
line([tFin+t(end),tFin+t(end)],[33,E(end)],'lineStyle','--','color','k')
tFin=tFin+t(end);

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R,
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
E=9.8*R*(1+gamma*cos(x(:,1)))+R^2*x(:,2).^2/2+(x(:,4).^2+
2*gamma*R*(cos(x(:,1)).*x(:,4)).*x(:,2))/2;
plot(t+tFin,E)
line([tFin+t(end),tFin+t(end)],[33,E(end)],'lineStyle','--','color','k')

hold off
grid on
xlabel('t')
ylabel('E')
title('Aro-partícula que ruedan')

La energía se conserva cuando el cuerpo rueda sin deslizar y disminuye cuando desliza

La reacción del plano

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial

%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
hold on
fplot(N,[0,aDesliza])
line([aDesliza,aDesliza],[0,N(aDesliza)],'lineStyle','--','color','k')

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R, 
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
N=9.8-gamma*R*sin(x(:,1)).*(9.8-gamma*R*cos(x(:,1)).*x(:,2).^2).*
(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))./(R*(1-gamma^2+gamma*sin(x(:,1)).
*(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))))-gamma*R*cos(x(:,1)).*x(:,2).^2;
plot(x(:,1),N)
line([x(end,1),x(end,1)],[0,N(end)],'lineStyle','--','color','k')

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);
fplot(N,[x(end,1),aDesliza])
line([aDesliza,aDesliza],[0,N(aDesliza)],'lineStyle','--','color','k')

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R,
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
N=9.8-gamma*R*sin(x(:,1)).*(9.8-gamma*R*cos(x(:,1)).*x(:,2).^2).*
(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))./(R*(1-gamma^2+
gamma*sin(x(:,1)).*(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))))-
gamma*R*cos(x(:,1)).*x(:,2).^2;
plot(x(:,1),N)
line([x(end,1),x(end,1)],[0,N(end)],'lineStyle','--','color','k')

hold off
grid on
set(gca,'XTick',0:pi/4:2*pi)
set(gca,'XTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
xlabel('\theta')
ylabel('N')
title('Aro-partícula que ruedan')

La reacción del plano horizontal presenta un máximo acusado para θ=180°

Fuerza de rozamiento

gamma=2/3;
R=1;  %radio 
mu=0.6; %coeficiente estático
w0=sqrt(1.1*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial

%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
hold on
fplot(Fr,[0,aDesliza])
line([aDesliza,aDesliza],[0,Fr(aDesliza)],'lineStyle','--','color','k')

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R, 
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
N=9.8-gamma*R*sin(x(:,1)).*(9.8-gamma*R*cos(x(:,1)).*x(:,2).^2).*
(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))./(R*(1-gamma^2+
gamma*sin(x(:,1)).*(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))))-
gamma*R*cos(x(:,1)).*x(:,2).^2;
plot(x(:,1),-mu*N)
line([x(end,1),x(end,1)],[0,-mu*N(end)],'lineStyle','--','color','k')

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) R*(1+gamma*cos(x)).*acel(x)-gamma*R*sin(x).*w2(x);
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);
fplot(Fr,[x(end,1),aDesliza])
line([aDesliza,aDesliza],[0,Fr(aDesliza)],'lineStyle','--','color','k')

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R,
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
N=9.8-gamma*R*sin(x(:,1)).*(9.8-gamma*R*cos(x(:,1)).*x(:,2).^2).*
(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))./(R*(1-gamma^2+
gamma*sin(x(:,1)).*(gamma*sin(x(:,1))+mu*(1+gamma*cos(x(:,1))))))-
gamma*R*cos(x(:,1)).*x(:,2).^2;
plot(x(:,1),-mu*N)
line([x(end,1),x(end,1)],[0,-mu*N(end)],'lineStyle','--','color','k')

hold off
grid on
set(gca,'XTick',0:pi/4:2*pi)
set(gca,'XTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
xlabel('\theta')
ylabel('F_r')
title('Aro-partícula que ruedan')

Como en el gráfico de la aceleración, la primera característica que observamos es la discontinuidad de ambas funciones para θ=160°, es decir, cuando comienza a rodar después de deslizar

La fuerza de rozamiento primero es negativa y después positiva. Cuando el cuerpo desliza la fuerza de rozamiento es igual a μN. Cuando el cuerpo rueda sin deslizar, la fuerza de rozamiento se obtiene de las ecuaciones del movimiento, se hace cero para θ=180°, cambiando de signo, de negativa a positiva

Otros ejemplos

Da una vuelta completa

Los valores de los parámetros son: γ=3/4, μ=0.4, el cuadrado de la velocidad angular adimensional $\frac{R}{g}{\left(\frac{d\theta }{dt}\right)}_0^2=0.5$

En primer lugar, representamos el cociente F_r/N, para que dado el coeficiente estático, μ, conozcamos el intervalo en el que se encuentra el ángulo θ_s raíz de la ecuación trascendente |F_r/N|=μ y se lo suministremos a la función fzero de MATLAB

gamma=3/4;
R=1;  %radio 
w0=sqrt(0.5*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial

%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) gamma*sin(x).*(9.8-R*w2(x))/2;
f=@(x) Fr(x)./N(x);

mu=0.4; %coeficiente estático
g=@(x) f(x)+mu;
hold on
fplot(f,[0,2*pi])
line([0,2*pi],[-mu,-mu],'color','b')
aDesliza=fzero(g,[pi/3,pi/2]);
plot(aDesliza,f(aDesliza),'bo','markersize',4,'markeredgecolor',
'b','markerfacecolor','b')

mu=0.1; %coeficiente estático
g=@(x) f(x)-mu;
line([0,2*pi],[mu,mu],'color','r')
aDesliza=fzero(g,[0.1,pi/4]);
plot(aDesliza,f(aDesliza),'ro','markersize',4,'markeredgecolor',
'r','markerfacecolor','r')
hold off
set(gca,'XTick',0:pi/4:2*pi)
set(gca,'XTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
grid on
xlabel('\theta')
ylabel('f')
title('Cociente, Fr/N')

Por ejemplo, para μ=0.4, la raíz se encuentra en la intersección de la función F_r/N y la recta -μ (en color azul). Para μ=0.1, la raíz se encuentra en la intersección entre F_r/N y la recta μ (en color rojo). En el primer caso, la fuerza de rozamiento es negativa y en el segundo, positiva. Como vemos en la gráfica tenemos que elegir cuidadosamente el intervalo en el que se encuentra la raíz para que fzero pueda calcularla

Creamos un script para representar la velocidad angular de la partícula dθ/dt y la velocidad del centro del aro dX/dt en función de la posición angular θ de la partícula en las distintas etapas del movimiento del cuerpo. El lector puede ampliar los scripts de la sección anterior para representar otras magnitudes en función de la posición θ o del tiempo t

gamma=3/4;
R=1;  %radio 
mu=0.4; %coeficiente estático
w0=sqrt(0.5*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial

%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) gamma*sin(x).*(9.8-R*w2(x))/2;
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
x0=[0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
hold on
plot(x(:,1),x(:,2),'r')
line([x(end,1),x(end,1)],[0,x(end,2)],'lineStyle','--','color','k')
line([0,x(end,1)],[x(end,2),x(end,2)],'lineStyle','--','color','k')
disp(aDesliza*180/pi)

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R, 
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
plot(x(:,1),x(:,2),'r')
plot(x(:,1),x(:,4),'b')
line([x(end,1),x(end,1)],[0,x(end,2)],'lineStyle','--','color','k')
line([0,x(end,1)],[x(end,2),x(end,2)],'lineStyle','--','color','k')
disp(x(end,1)*180/pi)

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) gamma*sin(x).*(9.8-R*w2(x))/2;
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);
x0=[th_0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
plot(x(:,1),x(:,2),'r')
line([x(end,1),x(end,1)],[0,x(end,2)],'lineStyle','--','color','k')
line([0,x(end,1)],[x(end,2),x(end,2)],'lineStyle','--','color','k')
disp(aDesliza*180/pi)

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, 
R,sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
plot(x(:,1),x(:,2),'r')
plot(x(:,1),x(:,4),'b')
line([x(end,1),x(end,1)],[0,x(end,2)],'lineStyle','--','color','k')
line([0,x(end,1)],[x(end,2),x(end,2)],'lineStyle','--','color','k')
disp(x(end,1)*180/pi)

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) gamma*sin(x).*(9.8-R*w2(x))/2;
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=2*pi;
%aDesliza=fzero(g,[x(end,1),2*pi]);
x0=[th_0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
plot(x(:,1),x(:,2),'r')
line([x(end,1),x(end,1)],[0,x(end,2)],'lineStyle','--','color','k')
line([0,x(end,1)],[x(end,2),x(end,2)],'lineStyle','--','color','k')
disp(aDesliza*180/pi)


hold off
grid on
set(gca,'XTick',0:pi/4:2*pi)
set(gca,'XTickLabel',{'0','\pi/4','\pi/2','3\pi/4','\pi','5\pi/4',
'3\pi/2','7\pi/4','2*\pi'})
xlabel('\theta')
ylabel('d\theta/dt,dX/dt')
title('Aro-partícula que ruedan')

  78.9025
  168.8163
  218.1270
  324.5845
  360

Observamos las distintas etapas del movimiento del cuerpo en la representación gráfica de las velocidades en función de la posición angular θ de la partícula.

• Rueda sin deslizar, desde 0 hasta 78.9°
• Desliza, desde 78.9° hasta 168.8°
• Rueda sin deslizar, desde 168.8° hasta 218.1°
• Desliza, desde 218.1° hasta 324.6°
• Rueda sin deslizar, desde 324.6° hasta 360°

Da una vuelta completa, empezando en la posición θ₀=0, con velocidad angular inicial (dθ/dt)₀ igual a

>> x(end,2)
ans =    2.1024

El movimiento completo ha comenzado con una velocidad angular inicial de 2.21 y ha terminado con una velocidad angular de 2.10

Comprobamos que la energía se mantiene constante cuando rueda sin deslizar y disminuye cuando desliza

gamma=3/4;
R=1;  %radio 
mu=0.4; %coeficiente estático
w0=sqrt(0.5*9.8/R); %velocidad angular inicial
th_0=0;  %posición inicial

%rueda sin deslizar
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x) gamma*sin(x).*(9.8-R*w2(x))/2;
f=@(x) Fr(x)./N(x);
g=@(x) f(x)+mu;
aDesliza=fzero(g,pi/4);
x0=[0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
hold on
E=R*(1+gamma*cos(x(:,1))).*(9.8+R*x(:,2).^2);
plot(t,E)
line([t(end),t(end)],[24,E(end)],'lineStyle','--','color','k')
tFin=t(end);

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))*
(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R, 
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
E=9.8*R*(1+gamma*cos(x(:,1)))+R^2*x(:,2).^2/2+(x(:,4).^2+
2*gamma*R*(cos(x(:,1)).*x(:,4)).*x(:,2))/2;
plot(t+tFin,E)
line([tFin+t(end),tFin+t(end)],[24,E(end)],'lineStyle','--','color','k')
tFin=tFin+t(end);


%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x)  gamma*sin(x).*(9.8-R*w2(x))/2;
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=fzero(g,[x(end,1),3*pi/2]);
x0=[th_0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
E=R*(1+gamma*cos(x(:,1))).*(9.8+R*x(:,2).^2);
plot(t+tFin,E)
line([tFin+t(end),tFin+t(end)],[24,E(end)],'lineStyle','--','color','k')
tFin=tFin+t(end);

%desliza
w0=sqrt(w2(aDesliza));
x0=[aDesliza,w0,R*aDesliza,R*w0];
mu=-mu;
f=@(t,x) [x(2); (9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1))))/(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+
mu*(1+gamma*cos(x(1)))))); x(4);-mu*9.8-gamma*R*(-mu*sin(x(1))+cos(x(1)))
*(9.8-gamma*R*cos(x(1))*x(2)^2)*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))/
(R*(1-gamma^2+gamma*sin(x(1))*(gamma*sin(x(1))+mu*(1+gamma*cos(x(1))))))+
gamma*R*(sin(x(1))+mu*cos(x(1)))*x(2)^2 ]; 
opts=odeset('events',@(t,x) stop_aro_particula1(t,x, gamma, mu, R,
sign(Fr(aDesliza))));
[t,x]=ode45(f,[0,50],x0,opts);
E=9.8*R*(1+gamma*cos(x(:,1)))+R^2*x(:,2).^2/2+(x(:,4).^2+
2*gamma*R*(cos(x(:,1)).*x(:,4)).*x(:,2))/2;
plot(t+tFin,E)
line([tFin+t(end),tFin+t(end)],[24,E(end)],'lineStyle','--','color','k')
tFin=tFin+t(end);

%rueda sin deslizar
w0=x(end,2); %velocidad angular inicial
th_0=x(end,1);  %posición inicial
w2=@(x) (1+gamma*cos(th_0))*(9.8/R+w0^2)./(1+gamma*cos(x))-9.8/R;
acel=@(x) gamma*sin(x).*(9.8+R*w2(x))./(2*R*(1+gamma*cos(x)));
N=@(x) 9.8-gamma*R*sin(x).*acel(x)-gamma*R*cos(x).*w2(x);
Fr=@(x)  gamma*sin(x).*(9.8-R*w2(x))/2;
f=@(x) Fr(x)./N(x);
g=@(x) f(x)-mu;
aDesliza=2*pi;
x0=[th_0,w0];
f=@(t,x) [x(2); gamma*sin(x(1))*(9.8+R*x(2)^2)/(2*R*(1+gamma*cos(x(1))))]; 
opts=odeset('events',@(t,x) stop_aro_particula(t,x,aDesliza));
[t,x]=ode45(f,[0,50],x0,opts);
E=R*(1+gamma*cos(x(:,1))).*(9.8+R*x(:,2).^2);
plot(t+tFin,E)
line([tFin+t(end),tFin+t(end)],[24,E(end)],'lineStyle','--','color','k')

hold off
grid on
xlabel('t')
ylabel('E')
title('Aro-partícula que ruedan')

El procedimiento ode45 no calcula correctamente el tercer intervalo, cuya energía debería ser constante

Referencias

W.F.D. Theron, M. F. Maritz. The amazing variety of motions of a loaded hoop. Mathematical and Computer Modelling, 47 (2008) 1077-1088