Movimiento sobre una cúpula en forma de catenaria

con rozamiento

En la página titulada la catenaria hemos deducido la ecuación de esta famosa curva que describe la forma que adopta un cable suspendido de los extremos de dos postes iguales. La ecuación de una catenaria en forma de cúpula con el eje Y de simetría es

y= 1 2γ ( cosh(γa)cosh(2γx) )

Representamos esta catenaria de longitud L, cuyos extremos están fijos en los puntos (-a/2,0) y (a/2,0). El parámetro γ se calcula a partir de la longitud L del cable y de la 'luz' (distancia entre extremos) a, tal como vemos en el siguiente script de MATLAB

L=1; %longitud de la catenaria
a=0.5; %'luz'
f=@(x) sinh(a*x)-L*x;
gamma=fzero(f,[0.1 100]);
f=@(x) (cosh(gamma*a)-cosh(2*gamma*x))/(2*gamma);
fplot(f,[-a/2,a/2]);
grid on
xlabel('x');
ylabel('y');
title('Catenaria')

Para este problema conviene establecer la ecuación de la trayectoria en forma paramétrica, tal como se hizo para la cúpula en forma de cicloide.

Llamamos 2γ=1/R, x=R·φ

y=R(cosh(γa)-cosh(φ))

Para φ=0, y=R por lo que cosh(γa)=2,

El valor mayor de φ se obtiene con a/2=0, , cosh(φ0)=2, φ0=1.317

>> acosh(2)
ans =    1.3170

Las ecuaciones paramétricas de la catenaria son

x=Rφ
y=R(2-cosh(φ))

Donde R y φ son dos parámetros

Cuando φ=0, x=0, y=R, es el vértice de la catenaria, cuando φ=φ0, x=R·φ0, y=0 es el extremo de la catenaria, tal como se ve en la figura

R=1;
fplot(@(x) R*x,@(x) R*(2-cosh(x)),[0,acosh(2)]);
grid on
xlabel('x');
ylabel('y');
title('Catenaria')

La pendiente de la recta tangente a la curva en la posición (x, y), tal como vemos en la figura.

tanθ= dy dx = Rsinhφ·dφ R·dφ =sinhφ

Omitimos el signo y relacionamos el parámetro φ y el ángulo (positivo) de la pendiente, sinhφ=tanθ

En el vértice φ=0, la pendiente θ=0. En el extremo, φ0=arccosh(2), tanθ0= 3 , θ0=1.0472

Calculamos la longitud del arco s, entre el vértice y el punto de coordenadas (x, y)

d s 2 =d x 2 +d y 2 =( R 2 + R 2 sinh 2 φ ) ( dφ ) 2 = R 2 cosh 2 φ· ( dφ ) 2 ds=Rcoshφ·dφ s= 0 φ Rcoshφ·dφ=Rsinhφ=Rtanθ ds= R cos 2 θ dθ

El radio de curvatura ρ de la cicloide en la posición (x, y) se calcula mediante la fórmula

ρ= ( 1+ ( dy dx ) 2 ) 3/2 d 2 y d x 2

Calculamos la derivada primera dy/dx y segunda d2y/dx2

dy dx =sinhφ d dx ( dy dx )= d dφ ( dy dx ) 1 dx dφ = coshφ R

El resultado es

ρ= ( 1+ sinh 2 φ ) 3/2 coshφ R =R cosh 2 φ=R( 1+ tan 2 θ )= R cos 2 θ

El mismo resultado se obtiene a partir de la definición de radio de curvatura, ds=ρ·dθ

ρ= ds dθ = R cos 2 θ

El radio de curvatura ρ es una cantidad positiva.

Dibujamos el arco de catenaria en el intervalo 0≤φ≤arccosh(2), calculamos la expresión del radio de curvatura, dibujamos la tangente y la normal a la catenaria en el punto correspondiente al parámetro φ=arccosh(2)/2. El extremo de la normal es el centro de curvatura (no se ve en la figura), pero puede verse desplazando la gráfica mediante la herramienta 'Pan'

syms phi xx;
R=1;
x=R*phi;
y=R*(2-cosh(phi));
yp=diff(y,phi)/diff(x,phi);
ypp=diff(yp,phi)/diff(x,phi);
r=(1+yp^2)^(3/2)/ypp;
hold on
ezplot(x,y,[0,acosh(2)])
grid on
x1=double(subs(x,phi,acosh(2)/2));
y1=double(subs(y,phi,acosh(2)/2));
rho=double(subs(r,phi,acosh(2)/2));
m1=double(subs(yp,phi,acosh(2)/2)); %pendiente tangente
m2=-1/m1; %pendiente normal
xc=x1+rho/sqrt(m2^2+1); %centro de curvatura
yc=y1+m2*(xc-x1);
y=m2*xx+y1-m2*x1;
ezplot(y,[xc,x1]) %normal
y=m1*xx+y1-m1*x1;
ezplot(y,[x1-0.5,x1+0.5]) %tangente

axis equal
xlim([0,acosh(2)])
ylim([0,1])
hold off
xlabel('x')
ylabel('y')
title('Catenaria')

Ecuación del movimiento

Dibujamos las fuerzas que actúan sobre la partícula, cuando se encuentra en la posición (x, y).

La partícula describe un movimiento con aceleración tangencial at y aceleración normal an. Aplicando la segunda ley de Newton

Teniendo en cuenta las definiciones de aceleración tangencial at y aceleración normal an,

a t = dv dt a n = v 2 ρ

Eliminando la reacción N de la superficie en el punto de contacto, expresamos las dos ecuaciones del movimiento en una única ecuación diferencial

dv dt =gsinθμgcosθ+μ v 2 ρ

La velocidad v y la aceleración dv/dt se expresan en términos del parámetro θ, (el ángulo de la recta tangente a la trayectoria) y sus derivadas respecto del tiempo.

v= ds dt = R cos 2 θ dθ dt dv dt = 2Rsinθ cos 3 θ ( dθ dt ) 2 + R cos 2 θ d 2 θ d t 2

Se obtiene la ecuación diferencial de segundo orden

d 2 θ d t 2 +( 2tanθμ ) ( dθ dt ) 2 g R ( sinθμcosθ ) cos 2 θ=0

que se resuelve por procedimientos numéricos con las siguientes condiciones iniciales t=0, θ=0, ( dθ dt ) 0 = v 0 R . Donde v0 es la velocidad inicial. En un apartado más abajo, resolveremos utilizando MATLAB la ecuación del movimiento. En la simulación se utiliza el procedimiento de Runge-Kutta

La partícula se detiene o deja de tener contacto con la cúpula

La ecuación del movimiento se puede integrar para obtener una expresión de la velocidad v en función del parámetro θ.

dv dt =gsinθμgcosθ+μ v 2 ρ v dv ds =gsinθμgcosθ+μ v 2 cos 2 θ R 1 2 d v 2 Rdθ cos 2 θ=gsinθμgcosθ+μ v 2 cos 2 θ R

Llamando x=v2/(Rg), nos queda la ecuación diferencial

dx dθ 2μx= 2 cosθ ( tanθμ )

La solución de la ecuación diferencial se compone de dos términos: la solución de la ecuación diferencial homogénea y la particular

La solución de la ecuación diferencial homogénea es

dx dθ 2μx=0 dx x =2μ·dθ dx x = 2μ·dθ

Integrando ambos miembros obtenemos lnx=2μθ+cte, o bien,x1=D·exp(2μθ)

Buscamos una solución particular de la forma, donde u(θ) es deconocida

x 2 =D e 2μθ u(θ)

Introducimos x2 en la ecuación diferencial

2μD e 2μθ u(θ)+D e 2μθ du dθ 2μD e 2μθ u(θ)= 2 cosθ ( tanθμ ) D e 2μθ du dθ = 2 cosθ ( tanθμ ) u(θ)= 2 D 0 θ e 2μθ cosθ ( tanθμ ) dθ x 2 =2 e 2μθ 0 θ e 2μθ cosθ ( tanθμ ) dθ

La solución de la ecuación diferencial, es la suma de la solución de la homogénea y de la particular x=x1+x2

v 2 Rg =D e 2μθ +2 e 2μθ 0 θ e 2μθ cosθ ( tanθμ ) dθ

La constante D se determina a partir de las condiciones iniciales, para θ=0, v=v0

Finalmente, la ecuación que nos proporciona la velocidad v en función del parámetro θ, es

v 2 Rg = v 0 2 Rg e 2μθ +2 e 2μθ 0 θ e 2μθ cosθ ( tanθμ ) dθ

La reacción N de la cúpula se calcula conociendo v2 en la posición dada por el parámetro θ

mgcosθN=m v 2 ρ N mg =cosθ( 1( v 2 Rg )cosθ )

Si no hay rozamiento, μ=0

v 2 Rg = v 0 2 Rg +2 0 θ sinθ cos 2 θ dθ v 2 Rg = v 0 2 Rg +2( 1 cosθ 1 )

La misma expresión se obtiene aplicando el principio de conservación de la energía.

1 2 m v 0 2 +mgR= 1 2 m v 2 +mgy v 2 = v 0 2 +2g( RR(2coshφ) ) v 2 Rg = v 0 2 Rg +2( 1+ sinh 2 φ 1 ) v 2 Rg = v 0 2 Rg +2( 1+ tan 2 θ 1 ) v 2 Rg = v 0 2 Rg +2( 1 cosθ 1 )

La partícula, deja de tener contacto con la cúpula para el valor del parámetro θs tal que N=0.

mgcosθN=m v 2 ρ mgcosθ=m v 2 cos 2 θ R 1 cos θ s =2 v 0 2 Rg

Si la velocidad inicial v0=0, θs=π/3 (60°). En el caso de una cúpula semiesférica θs=arccos(2/3) (48°), en el caso de una cúpula en forma de cicloide, θs=π/4 (45°). La partícula permanece en el vértice θ=0, que es una posición de equilibrio inestable, hasta que una pequeña perturbación la lleva a deslizar a lo largo de la cúpula, tal como se ha descrito en la página Movimiento sobre una cúpula semiesférica, sin rozamiento

Una vez alcanzada la posición θs, la partícula describe un movimiento parabólico. Dicha posición (x,y) se calcula mediante la relación entre la pendiente θ y el parámetro φ y ecuaciones paramétricas de la catenaria

sinhφ=tanθ
x=Rφ
y=R(2-cosh(φ))=R(2-1/cosθ)

Representamos la reacción N/mg en función del parámetro θ hasta el ángulo crítico θs, posición en la que se anula N=0

R=1;
v0=2; %velocidad inicial
v2=@(x) v0^2+2*R*9.8*(1./cos(x)-1); %cuadrado velocidad
N=@(x) cos(x).*(1-v2(x).*cos(x)/(9.8*R)); %normal

th_c=acos(1/(2-v0^2/(R*9.8))); %ángulo crítico N=0
fplot(N,[0,th_c])
grid on
xlabel('\theta')
ylabel('N/mg')
title('Normal')

Si hay rozamiento μ≠0, pueden ocurrir dos casos:

Sea una cúpula en forma de catenaria de R=1, la velocidad inicial en el vértice, θ=0, es v0=2 m/s, el coeficiente de rozamiento es μ=0.2. Representamos v2/(Rg) y N en función del parámetro θ

R=1;
mu=0.2; %coeficiente de rozamiento
v0=2; %velocidad inicial
f=@(x) (exp(-2*mu*x).*(tan(x)-mu))./cos(x);
v2=@(x) exp(2*mu*x).*(v0^2/(R*9.8)+2*integral(f,0,x)); %cuadrado velocidad
N=@(x) cos(x).*(1-v2(x).*cos(x)); %normal

subplot(2,1,1)
fplot(v2,[0,1])
grid on
ylim([0,1.5])
xlabel('\theta')
ylabel('v^2/Rg')
title('Velocidad')

subplot(2,1,2)
fplot(N,[0,1])
grid on
ylim([-0.5,1])
xlabel('\theta')
ylabel('N/mg')
title('Normal')

La velocidad disminuye y luego aumenta. La reacción N se hace cero para un valor crítico del parámetro θc, que calculamos mediante la función fzero de MATLAB. Recuérdese que la valor mayor de θ es tanθ0= 3 , θ0=1.0472

>> fzero(N,[0,atan(sqrt(3))])
ans =    0.9758

Cambiamos el coeficiente de rozamiento a μ=0.4 y la velocidad inicial v0=1. Representamos v2/(Rg) en función del parámetro θ

R=1;
mu=0.4; %coeficiente de rozamiento
v0=1; %velocidad inicial
f=@(x) (exp(-2*mu*x).*(tan(x)-mu))./cos(x);
v2=@(x) exp(2*mu*x).*(v0^2/(R*9.8)+2*integral(f,0,x)); %cuadrado velocidad

fplot(v2,[0,1.0472])
ylim([-0.5,1])
grid on
xlabel('\theta')
ylabel('v^2/Rg')
title('Velocidad')

La partícula se detiene v=0, en la posición dada por el parámetro θc que calculamos utilizando la función fzero de MATLAB.

>> fzero(v2,[0,0.2])
ans =    0.1725

La velocidad inicial en el vértice, θ=0, es v0=2 m/s, el coeficiente de rozamiento es μ=0.5. Representamos v2/(Rg) y N en función del parámetro θ

R=1;
mu=0.5; %coeficiente de rozamiento
v0=2; %velocidad inicial
f=@(x) (exp(-2*mu*x).*(tan(x)-mu))./cos(x);
v2=@(x) exp(2*mu*x).*(v0^2/(R*9.8)+2*integral(f,0,x)); %cuadrado velocidad
N=@(x) cos(x).*(1-v2(x).*cos(x)); %normal

subplot(2,1,1)
fplot(v2,[0,atan(sqrt(3))])
grid on
ylim([0,1.5])
xlabel('\theta')
ylabel('v^2/Rg')
title('Velocidad')

subplot(2,1,2)
fplot(N,[0,atan(sqrt(3))])
grid on
ylim([0,1.5])
xlabel('\theta')
ylabel('N/mg')
title('Normal')

Ni la velocidad v ni la normal N se hacen cero para 0≤θθ0

Velocidad inicial crítica y ángulo crítico

Fijamos el coeficiente de rozamiento, por ejemplo, μ=0.4, y vamos cambiando la velocidad inicial v0 con el que lanzamos la partícula en la posición θ=0.

mu=0.4; %coeficiente de rozamiento
R=1;
hold on
for v0=[1.1,1.3,1.5] %velocidades iniciales
    f=@(x) (exp(-2*mu*x).*(tan(x)-mu))./cos(x);
    v2=@(x) exp(2*mu*x).*(v0^2/(R*9.8)+2*integral(f,0,x)); %cuadrado velocidad
    fplot(v2,[0,1.0472],'displayName',num2str(v0))
end
hold off
legend('-DynamicLegend','location','northeast')
line([atan(mu),atan(mu)],[-0.05,0.2],'linestyle','--');
ylim([-0.05,0.6])
grid on
xlabel('\theta')
ylabel('v^2/Rg')
title('Velocidades')

En la figura, se representa el cuadrado de la velocidad de la partícula v2/(Rg) en función del parámetro θ, para varias velocidades iniciales v0=1.1, 1.3, 1.5 m/s. Se ha fijado el parámetro R=1.0 m

Existe una velocidad inicial crítica v0,c para el cual la velocidad v de la partícula presenta un mínimo v=0 en el ángulo θc denominado ángulo crítico.

v 2 ( θ c )=0 ( d v 2 dθ ) θ c =0

Tenemos un sistema de dos ecuaciones, haciendo operaciones y simplificando

{ v 0 2 Rg +2 0 θ c e 2μθ cosθ ( tanθμ ) dθ=0 2μ v 0 2 Rg e 2μ θ c +4μ e 2μ θ c 0 θ c e 2μθ cosθ ( tanθμ ) dθ+2 tan θ c μ cos θ c =0

llegamos a la siguiente relación

tanθc=μ

Calculamos la velocidad inicial crítica v0,c en θ=0 para la cual la partícula se detiene en θc.

v 0c 2 Rg =2 0 arctanμ e 2μθ cosθ ( tanθμ ) dθ

Representamos el cuadrado de la velocidad inicial crítica, en función del coeficiente de rozamiento μ

R=1;
roza=linspace(0,4,100);
Vc=zeros(0,length(roza));
i=1;
for mu=roza
    f=@(x) (exp(-2*mu*x).*(tan(x)-mu))./cos(x);  
    Vc(i)=-2*integral(f,0,atan(mu)); %cuadrado velocidad
    i=i+1;
end
plot(roza,Vc)
grid on
xlabel('\mu')
ylabel('v^2/Rg')
title('Velocidad crítica')

La curva divide el plano (μ, v0) en dos regiones. Fijado el valor de μ

Balance energético

La energía de la partícula en la posición inicial es

Cuando la partícula se encuentra en la posición definida por el parámetro es θ.

La fuerza de rozamiento Fr=μN, tiene la misma dirección (tangencial) que el desplazamiento, el arco ds, pero de sentido contrario. El trabajo realizado por la fuerza de rozamiento es.

W nc = 0 θ μN·ds=μm 0 θ ( gcosθ v 2 ρ )ds= μm 0 θ ( gcosθ v 2 cos 2 θ R ) R cos 2 θ ·dθ=μmgR 0 θ ( 1 cosθ v 2 Rg )dθ

El trabajo realizado por la fuerza de rozamiento Wnc<0 es igual a la diferencia entre la energía final menos la energía inicial.

W nc =( 1 2 m v 2 +mgR( 1 cosθ 1 ) )( 1 2 m v 0 2 +mgR )

Resulta

v 2 Rg v 0 2 Rg =4 2 cosθ 2μ 0 θ ( 1 cosθ v 2 Rg )dθ x x 0 =4 2 cosθ 2μ 0 θ ( 1 cosθ x )dθ

Derivando con respecto a θ

dx dθ = 2sinθ cos 2 θ 2μ( 1 cosθ x ) dx dθ 2μx= 2 cosθ ( tanθμ )

La misma ecuación diferencial que obtuvimos por dinámica

Solución numérica de la ecuación del movimiento

Primero, calculamos el ángulo límite para el cual la partícula se detiene v=0, o bien, la reacción de la cúpula N=0. Después, resolvemos la ecuación diferencial del movimiento que se ha deducido al principio de esta página.

mu=0.4; %coeficiente de rozamiento
v0=2; %velocidad inicial para theta=0,
R=1;

f=@(x) (exp(-2*mu*x).*(tan(x)-mu))./cos(x);
Vc=sqrt(-2*R*9.8*integral(f,0,atan(mu))); % velocidad crítica
v2=@(x) exp(2*mu*x).*(v0^2/(R*9.8)+2*integral(f,0,x)); %cuadrado velocidad
N=@(x) cos(x).*(1-v2(x).*cos(x)); %normal/mg

if v0<Vc
    aLimite=fzero(v2,atan(mu));
    opts=odeset('events',@(t,x) stop_cupula_1(t,x));
    disp('velocidad nula')
else
    th_c=acos(1/(2-v0^2/(R*9.8)));
    aLimite=fzero(N,th_c);
    opts=odeset('events',@(t,x) stop_cupula(t,x,aLimite));
    disp('normal nula')
end

x0=[0, v0/R];
f=@(t,x) [x(2);-(2*tan(x(1))-mu)*x(2)^2+9.8*(sin(x(1))-
mu*cos(x(1)))*cos(x(1))^2/R];
tspan=[0,20]; 
[t,x]=ode45(f,tspan,x0,opts);

plot(t,x(:,1))
grid on
xlabel('t')
ylabel('\theta')
title('Catenaria con rozamiento')

Definimos dos funciones para que el proceso de integración se detenga cuando N sea cero o cuando la derivada dθ/dt sea cero.

function [value,isterminal,direction]=stop_cupula(~,x,xFin)
    value=xFin-x(1); 
    isterminal=1;
    direction=0; 
end
function [value,isterminal,direction]=stop_cupula_1(~,x)
    value=x(2); 
    isterminal=1;
    direction=0; 
end

Cuando μ=0.2 y v0=2. La partícula deja de tener contacto con la cúpula, N=0 cuando θ=0.9758, tal como vemos en la gráfica (ordenda Y).

normal nula
>> aLimite
aLimite =    0.9758

Cambiamos la velocidad inicial v0=1.1, la partícula se detiene dθ/dt=0, para θ=0.2366.

velocidad nula
>> aLimite
aLimite =    0.2366

Actividades

Se introduce

Se pulsa el botón titulado Nuevo

Observamos el movimiento de la partícula deslizando sobre la cúpula. Sobre la partícula se dibujan las fuerzas: peso mg, reacción N, y fuerza de rozamiento Fr.

En la parte superior derecha, se proporcionan los datos de

Cuando la reacción N se hace cero, se muestran los datos de la a velocidad v en la posición θc en la que la partícula deja de tener contacto con la cúpula.

El círculo situado en la parte superior izquierda representa la energía total de la partícula, la porción de color rojo representa la energía cinética, y la porción azul, la energía potencial. Observamos que la energía potencial se va transformando en energía cinética, pero la suma de los valores de ambas clases de energía no se mantiene constante a lo largo de la trayectoria de la partícula si hay rozamiento. El trabajo de la fuerza de rozamiento viene indicado por la porción negra del círculo de mayor radio.

Mediante una línea roja a trazos se señala el valor límite del parámetro θ, para el cual la reacción N=0 o la partícula se detiene, su velocidad v=0


Referencias

Felipe González-Cataldo, Gonzalo Gutiérrez. Julio M. Yáñez. Sliding down an arbitrary curve in the presence of friction. Am. J. Phys. 85 (2) February 2017, 108-114