El ludión o diablillo de Descartes

El diablillo o ludión es un tubo de ensayo de longitud L, de diámetro interior d1 y de diámetro exterior d2 que se llena parcialmente de agua hasta una altura L-l0. Luego, se invierte en un recipiente grande que contiene agua, que está cerrado y conectado por una parte, a un manómetro para medir la presión del aire encerrado en su interior y por otra, a una jeringa que nos permite variar la presión P del aire del recipiente.

Estática del diablillo

En este apartado, estudiaremos el equilibrio del tubo de ensayo (diablillo) con la burbuja de aire dentro

Control de la presión

Para variar la presión del aire contenido en el recipiente se emplea una jeringa. Si el volumen total de aire del recipiente y de la jeringa es V0. Al disminuir en x el volumen de la jeringa la presión aumenta de P0 a P. Si suponemos una transformación isotérmica, tendremos que:

P(V0-x)=P0·V0

Si x<<V0 hacemos la siguiente aproximación

P= P 0 1x/ V 0 = P 0 ( 1 x V 0 ) 1 P 0 ( 1+ x V 0 )

El incremento de presión ΔP=P-P0 es proporcional a la variación del volumen x de la jeringa o bien, al desplazamiento de su émbolo.

El tubo de ensayo

Después de haber invertido el tubo de ensayo dentro del agua del recipiente, la longitud de la burbuja de aire que se ha formado en el interior del tubo invertido es l=x+z

La base del tubo se encuentra a una altura x por encima de la superficie del agua del recipiente y el agua en el tubo se encuentra a una altura z por debajo de dicha superficie tal como se indica en la primera figura.

Sobre el tubo actúan dos fuerzas: el peso y el empuje

Cuando el tubo está parcialmente sumergido x>0

Situación de equilibrio

El tubo permanece en equilibrio cuando el peso se igual al empuje o bien, cuando F=0.

Conocida la presión P determinamos x y z resolviendo un sistema de dos ecuaciones (1) y (2) con dos incógnitas.

De la segunda ecuación despejamos x

x= P 0 l 0 P+ρgz z     (3)

y la sustituimos en la primera, para despejar z, quedando una ecuación de segundo grado, az2+bz+c=0, con

a= ρ 2 g( V L +A ) b=ρP( V L +A )( ρ v ρ)Vρg c= V L ρ P 0 l 0 ( ρ v ρ)VP

Se calcula la raíz positiva de la ecuación de segundo grado

z= b+ b 2 4ac 2a

Una vez que se ha calculado z, se determina mediante la ecuación (3), la posición de equilibrio x de la parte superior del tubo de ensayo.

Presión crítica.

Si se incrementa la presión hasta un valor límite P, el tubo de ensayo se va sumergiendo en el agua hasta que la posición de la parte superior del tubo está en el origen xe=0. Si incrementamos un poco más la presión el tubo se hunde completamente.

La presión crítica P se determina poniendo x=0, en la ecuación (1) con F=0 (situación de equilibrio) y despejando P en la ecuación (2).

z= V A ( ρ v ρ 1 )P= P 0 l 0 z ρgz

El tamaño z de la burbuja es independiente de su tamaño l0 inicial, solamente depende de la geometría del tubo de ensayo, de la densidad del vidrio ρv y del agua ρ

Ejemplo

En el programa interactivo más abajo, se introduce

Se pulsa en el botón titulado Nuevo

Observaremos el tubo en la posición de equilibrio estable F=0, si es que existe. En caso contrario, un mensaje nos lo indica.

Al lado del tubo, se muestran dos vectores que representan el peso y al empuje, y en la parte superior del recipiente, se muestra el valor de la fuerza F.

Datos de los tubos utilizados, tomados del artículo citado en las referencias

Tubo Altura L (cm) Diámetro interior d1(cm) Diámetro exterior d2 (cm) Densidad ρv (g/cm3)
1 9.8 1.38 1.6 2.35
2 16.0 1.43 1.64 2.29
3 17.1 4.44 4.97 2.16
4 20.1 1.74 2.00 2.17

La presión atmosférica es P0=101300 Pa. La densidad del mercurio (Hg) es 13.55 g/cm3

Ejemplo de presión crítica con el Tubo 2.

Calculamos la presión crítica para los estados iniciales siguientes:

Altura inicial l0 de la burbuja de aire en cm Presión crítica P (Pa) Diferencia de presión ΔP=P-P0 en mm de Hg
4 61632 -299
6 92766 -64
8 123901 170

Los datos de la última columna son los que se introducen en el control titulado Presión cuando se establece con el ratón la altura inicial l0 de la primera columna, en el programa interactivo al final de la página. La velocidad inicial deberá ser cero, en el control titulado Velocidad

L=16/100; %longitud
d_int=1.43/100; %diámetro interior
d_ext=1.64/100; %diámetro exterior
V=pi*L*(d_ext^2-d_int^2)/4;
A=pi*d_int^2/4;
rho=1000; %densidad agua
rho_g=2.29*1000; %densidad vidrio
P0=1.013e5; %presión atmosférica
for l0=(4:2:8)/100
    z=V*(rho_g/rho-1)/A;
    P=P0*l0/z-rho*9.8*z;
    disp([P,(P-P0)*1000/(13550*9.8)])
end
   1.0e+04 *    6.1632   -0.0299
   1.0e+04 *    9.2766   -0.0064
   1.0e+05 *    1.2390    0.0017

Cuando el tubo está completamente sumergido x≤0

Dinámica del diablillo

Escribimos, de nuevo, la fuerza que actúa sobre el tubo de ensayo

F={ ρgAz+ρgV(1x/L) ρ v gVx>0 ρgAz+ρgV ρ v gVx0

Relacionamos z y x mediante la ley de Boyle

P 0 l 0 ={ (P+ρgz)(z+x)x>0 (P+ρg(zx))zx0

Tenemos que despejar z de las ecuaciones de segundo grado

ρg z 2 +(P+ρgx)z+Px P 0 l 0 =0x>0 ρg z 2 +(Pρgx)z P 0 l 0 =0x0 z= b+ b 2 4ac 2a a=ρg   b={ P+ρgx Pρgx c={ Px P 0 l 0 x>0 P 0 l 0 x0

La fuerza F es ahora, una función de la posición x de la base del tubo de ensayo invertido

F=ρgV( 1 x L ) ρ v gV+ A 2 ( P+ρgx+ ( ρgxP ) 2 +4ρg P 0 l 0 )

La fuerza F depende de la posición x, se trata de una fuerza conservativa cuya energía potencial se calcula del siguiente modo

x e x F dx= E p ( x e ) E p (x)

Donde xe es la posición de equilibrio al cual le vamos a asignar una energía potencial cero. Establecemos por tanto, en esta posición el nivel cero de energía potencial Ep(xe)=0.

Para x>0 tenemos que hacer la integral

E 1 (x)= F(x)dx+ C 1 E 1 (x)= ρg 2 ( A 2 + V L ) x 2 +( AP 2 +V( ρ v ρ) )x A 2 f(x)+ C 1

donde C1 es una constante de integración que se determina a partir de la condición de que E1(xe)=0,  f(x) es una función que se define más abajo.

La fuerza es

F=( ρ v ρ)gV+ A 2 ( Pρgx+ ( ρgxP ) 2 +4ρg P 0 l 0 )

Para calcular la energía potencial, tenemos que hacer la integral

E 2 (x)= F(x)dx+ C 2 E 2 (x)= Aρg 4 x 2 +( AP 2 +V( ρ v ρ) )x A 2 f(x)+ C 2

La fuerza F es una función discontinua en x=0, pero la energía potencial Ep(x) es una función continua de modo que, en x=0 se tiene que cumplir que E1(0)=E2(0), esta condición determina el valor de la constante C2 de integración.

La función f(x) es una integral que no es inmediata

f(x)= ( ρgxP ) 2 +4ρg P 0 l 0 ·dx

Esta integral se puede escribir de la forma

4 P 0 l 0 t 2 +1 ·dtt=( ρgxP 4ρg P 0 l 0 )

Se descompone en la suma de dos integrales

t 2 +1 dt= t 2 +1 t 2 +1 dt = t·t t 2 +1 dt+ dt t 2 +1 = t t 2 +1 +ln( t+ t 2 +1 ) 2

El primer término se integra por partes, y el segundo es una integral inmediata. El resultado final es

f(x)=( x 2 P 2ρg ) ( ρgxP ) 2 +4ρg P 0 l 0 + 2 P 0 l 0 ln( ρgxP+ ( ρgxP ) 2 +4ρg P 0 l 0 )

Una de las posibles representaciones gráficas de Ep(x) se muestra en la figura.

Si en la posición de equilibrio, xe le damos al tubo una velocidad inicial v0, la energía total del móvil es la energía cinética ya que la energía potencial es nula.

E= 1 2 m v 0 2 = 1 2 ( ρ v V) v 0 2

La energía E total es constante en todos los puntos de la trayectoria. Si E es menor que el máximo local de la energía potencial, el tubo que sale de la posición de equilibrio xe oscila entre dos posiciones, x1 y x2 determinadas por las raíces de la ecuación trascendente Ep(x)-E=0, abscisas de los puntos de intersección de la curva de la energía potencial y la recta horizontal Ep(x)=E, tal como se indica en al figura.

Cuando la función energía potencial Ep(x) no tiene mínimo local o posición de equilibrio estable, el tubo sale de la posición x=0 con velocidad v0 y se hunde en el recipiente hasta que llega al fondo.

Se simula el movimiento del diablillo de Descartes resolviendo numéricamente, la ecuación diferencial del movimiento, con las condiciones iniciales especificadas.

Cálculo con MATLAB

Elegimos el tubo número 2. La representación gráfica de la energía potencial Ep(x) para varios valores de la presión P es

L=16/100; %longitud
l0=7/100; %longitud vacía
d_int=1.43/100; %diámetro interior
d_ext=1.64/100; %diámetro exterior
V=pi*L*(d_ext^2-d_int^2)/4;
A=pi*d_int^2/4;
rho=1000; %densidad agua
rho_g=2.29*1000; %densidad vidrio
P0=1.013e5; %presión atmosférica
hold on
for P=(1:0.02:1.08)*P0
    a=rho^2*9.8*(V/L+A);
    b=rho*P*(V/L+A)-(rho_g-rho)*V*rho*9.8;
    c=-V*rho*P0*l0/L-(rho_g-rho)*V*P;
    z=(-b+sqrt(b^2-4*a*c))/(2*a);
    xe=P0*l0/(P+z*rho*9.8)-z; %equilibrio
    f=@(x) -(P/(2*rho*9.8)-x/2).*sqrt((P-rho*9.8*x).^2+4*rho*9.8*P0*l0)
+2*P0*l0*log(sqrt((P-rho*9.8*x).^2+4*rho*9.8*P0*l0)-(P-rho*9.8*x));
    U1=@(x) rho*9.8*(A/2+V/L)*x.^2/2+(A*P/2+V*(rho_g-rho)*9.8)*x-A*f(x)/2;
    U2=@(x) -A*rho*9.8*x.^2/4+(A*P/2+V*(rho_g-rho)*9.8)*x-A*f(x)/2;
    C1=-U1(xe);
    C2=-U2(xe);
    U=@(x) (U1(x)+C1).*(x>0)+(U2(x)+C2).*(x<=0);
    fplot(U,[-1.2,0.1], 'displayName',num2str(P/P0))
end
hold off
grid on
xlabel('x')
ylim([-0.002,0.004])
legend('-DynamicLegend','location','best')
ylabel('E_p(x)')
title('Energía potencial')

Movimiento del diablillo de energía E entre dos posiciones para una presión P=P0

L=16/100; %longitud
l0=7/100; %longitud vacía
d_int=1.43/100; %diámetro interior
d_ext=1.64/100; %diámetro exterior
V=pi*L*(d_ext^2-d_int^2)/4;
A=pi*d_int^2/4;
rho=1000; %densidad agua
rho_g=2.29*1000; %densidad vidrio
P0=1.013e5; %presión atmosférica
P=P0; %presión
a=rho^2*9.8*(V/L+A);
b=rho*P*(V/L+A)-(rho_g-rho)*V*rho*9.8;
c=-V*rho*P0*l0/L-(rho_g-rho)*V*P;
z=(-b+sqrt(b^2-4*a*c))/(2*a);
xe=P0*l0/(P+z*rho*9.8)-z; %equilibrio
f=@(x) -(P/(2*rho*9.8)-x/2).*sqrt((P-rho*9.8*x).^2+4*rho*9.8*P0*l0)+
2*P0*l0*log(sqrt((P-rho*9.8*x).^2+4*rho*9.8*P0*l0)-(P-rho*9.8*x));
U1=@(x) rho*9.8*(A/2+V/L)*x.^2/2+(A*P/2+V*(rho_g-rho)*9.8)*x-A*f(x)/2;
U2=@(x) -A*rho*9.8*x.^2/4+(A*P/2+V*(rho_g-rho)*9.8)*x-A*f(x)/2;
C1=-U1(xe);
C2=-U2(xe);
U=@(x) (U1(x)+C1).*(x>0)+(U2(x)+C2).*(x<=0);
fplot(U,[-1.2,0.1])

E=0.002; %energía
ff=@(x) U(x)-E;
line([-1.2,0.1],[E,E],'color','k')
x1=fzero(ff,[xe,0.1]); %posiciones de retorno
x2=fzero(ff,[-0.7,xe]);
line([x1,x1],[0,U(x1)],'color','k','lineStyle','--')
line([x2,x2],[0,U(x2)],'color','k','lineStyle','--')
line([x1,x2],[0,0],'color','r')
ylim([0,0.003])
grid on
xlabel('x')
ylabel('E_p(x)')
title('Energía potencial')

El diablillo de energía E=0.002 se puede mover entre las dos posiciones x1 y x2 en los extremos del segmento de color rojo. Estas dos posiciones son las raíces de la ecuación transcendente Ep(x)-E=0 y se obtienen mediante la función fzero de MATLAB

El diablillo de la situación anterior, de masa m, parte de la posición de equilibrio xe con velocidad v 0 = 2E/m . La energía potencial en la posición de equilibrio Ep(xe)=0 es nula

function ludion_1
    L=16/100; %longitud
    l0=7/100; %longitud vacía
    d_int=1.43/100; %diámetro interior
    d_ext=1.64/100; %diámetro exterior
    V=pi*L*(d_ext^2-d_int^2)/4;
    A=pi*d_int^2/4;
    rho=1000; %densidad agua
    rho_g=2.29*1000; %densidad vidrio
    P0=1.013e5; %presión atmosférica
    P=P0;
    m=18.6/1000; %masa
 
    a=rho^2*9.8*(V/L+A);
    b=rho*P*(V/L+A)-(rho_g-rho)*V*rho*9.8;
    c=-V*rho*P0*l0/L-(rho_g-rho)*V*P;
    z=(-b+sqrt(b^2-4*a*c))/(2*a);
    xe=P0*l0/(P+z*rho*9.8)-z; %equilibrio
    E=0.002; %energía
    v=sqrt(2*E/m); %velocidad
    [t,x]=ode45(@ludion_ode45,[0,10],[xe,-v]);
    plot(t,x(:,1))
    grid on
    xlabel('t')
    ylabel('x')
    title('Movimiento')
   
    function dr=ludion_ode45(~, x)
        a=rho*9.8;
        b=P+abs(x(1))*rho*9.8;
        c=P*x(1)-P0*l0;
        if x(1)<=0
            c=-P0*l0;
        end
        xi=(-b+sqrt(b^2-4*a*c))/(2*a);
        temp=A*xi*rho*9.8+V*(1-x(1)/L)*rho*9.8-V*rho_g*9.8;%-0.01*x(2);
        if x(1)<=0
            temp=A*xi*rho*9.8+V*(rho-rho_g)*9.8;%-0.01*x(2);
        end
        dr=[x(2); temp/m];
    end
end

Observamos que el procedimiento ode45 no resuelve bien la ecuación diferencial, debido a los cambios bruscos de velocidad cerca del origen x=0. Sin embargo, el procedimiento de Runge-Kutta empleado en la simulación, más abajo, resuelve bastante bien este sistema conservativo

Añadimos un rozamiento proporcional a la velocidad, λv, debido a las bajas velocidades con las que se mueve el diablillo en el seno del fluido. En la simulación, más abajo, supondremos que no hay rozamiento

function ludion_2
    L=16/100; %longitud
    l0=7/100; %longitud vacía
    d_int=1.43/100; %diámetro interior
    d_ext=1.64/100; %diámetro exterior
    V=pi*L*(d_ext^2-d_int^2)/4;
    A=pi*d_int^2/4;
    rho=1000; %densidad agua
    rho_g=2.29*1000; %densidad vidrio
    P0=1.013e5; %presión atmosférica
    P=P0;
    m=18.6/1000; %masa
 
    a=rho^2*9.8*(V/L+A);
    b=rho*P*(V/L+A)-(rho_g-rho)*V*rho*9.8;
    c=-V*rho*P0*l0/L-(rho_g-rho)*V*P;
    z=(-b+sqrt(b^2-4*a*c))/(2*a);
    xe=P0*l0/(P+z*rho*9.8)-z; %equilibrio
    E=0.002;
    v=sqrt(2*E/m); %velocidad
    [t,x]=ode45(@ludion_ode45,[0,10],[xe,-v]);
    plot(t,x(:,1))

    grid on
    xlabel('t')
    ylabel('x')
    title('Movimiento')
   
    function dr=ludion_ode45(~, x)
        a=rho*9.8;
        b=P+abs(x(1))*rho*9.8;
        c=P*x(1)-P0*l0;
        if x(1)<=0
            c=-P0*l0;
        end
        xi=(-b+sqrt(b^2-4*a*c))/(2*a);
        temp=A*xi*rho*9.8+V*(1-x(1)/L)*rho*9.8-V*rho_g*9.8-0.01*x(2);
        if x(1)<=0
            temp=A*xi*rho*9.8+V*(rho-rho_g)*9.8-0.01*x(2);
        end
        dr=[x(2); temp/m];
    end
end

Actividades

Se introduce

Se pulsa en el botón titulado Nuevo

Se pulsa en el botón titulado

Observamos el movimiento del tubo, oscilatorio, si la energía total E del tubo es menor que la del máximo local de la energía potencial. Se representa a la izquierda, la curva Ep(x), la energía total E, y los puntos x1 y x2 de retorno.

En el caso de que E sea mayor que el máximo local, el tubo se hunde en el recipiente y llega hasta el fondo. Lo mismo ocurre si la curva de energía potencial no tiene mínimo ni máximo local, es decir, la presión del aire contenido en el recipiente es superior a la crítica.

Si no hay posición de equilibrio estable, modificamos la presión o bien, modificamos el tamaño inicial l0 de la burbuja de aire actuando con el puntero del ratón en el círculo de color rojo.

Para comenzar una nueva “experiencia” se pulsa en el botón titulado Nuevo. Se utiliza los botones pausa || y paso a paso >| para medir la amplitud de la oscilación y el periodo



Referencias

Güémez J, Fiolhais C, Fiolhais M. The Cartesian diver and the fold catastrophe. Am. J. Phys. 70 (7) July 2002, pp. 710-714.