Oscilaciones del líquido contenido en dos vasos comunicantes

Consideremos dos recipientes de sección S1 y S2, comunicados por su parte inferior a través de un tubo de longitud d y pequeña sección S. Se cierra la llave de paso del tubo horizontal y se vierte un líquido en los recipientes hasta una altura h01 y h02 tal como se muestra en la figura

Cuando se abre la llave de paso y se comunican los dos recipientes, el líquido fluye a través del tubo desde un recipiente al otro, hasta que al cabo de un tiempo se alcanza una altura de equilibrio h

Como el volumen total de líquido no ha cambiado, desde la situación inicial a la de equilibrio, tendemos que

S1h01+S2h02=S1h+S2h

En esta página, se considera que el líquido que se vierte en los recipientes es ideal, el estado de equilibrio no se alcanza. El líquido pasa de un recipiente hacia el otro y vicersa, oscilando entre alturas constantes, h02 y 2h-h02 en el recipiente de la derecha.

En un instante dado t, el nivel del líquido por debajo de la altura de equilibrio en el primer vaso es x1 y en el segundo, el nivel está por encima x2. Como el volumen de líquido no ha cambiado

S1x1=S2x2

Si v1 es la velocidad del fluido en el primer recipiente, v2 en el segundo y u en el tubo que comunica ambos recipientes. Se cumplirá por la ecuación de continuidad que

S1v1=S2v2=Su

Energías

La energía cinética total es

E k = 1 2 ρ S 2 { ( 1+ S 2 S 1 )h+( 1 S 2 2 S 1 2 ) x 2 + S 2 S d } v 2 2

La posición del centro de masas del recipiente izquierdo y derecho, se señalan por un punto de color rojo en la figura. La posición del centro de masas del líquido contenido en los dos recipientes es

x c = S 1 (h x 1 )( h x 1 2 )+ S 2 (h+ x 2 )( h+ x 2 2 ) S 1 (h x 1 )+ S 2 (h+ x 2 ) = 1 2 ( h+ S 2 S 1 x 2 2 h )

La energía potencial del centro de masas es

E p =ρ( S 1 + S 2 )h·g· x c = 1 2 ρ( S 1 + S 2 )g( h 2 + S 2 S 1 x 2 2 )

Ecuación del movimiento

La lagrangiana L=Ek-Ep es

L= 1 2 ρ S 2 { ( 1+ S 2 S 1 )h+( 1 S 2 2 S 1 2 ) x 2 + S 2 S d } ( d x 2 dt ) 2 1 2 ρ S 2 ( 1+ S 1 S 2 )g( h 2 + S 2 S 1 x 2 2 )

Para simplificar la notación, denominaremos y=x2. La ecuación del movimiento es

d dt ( L y ˙ ) L y =0 d dt ( 2{ ( 1+ S 2 S 1 )h+( 1 S 2 2 S 1 2 )y+ S 2 S d } dy dt )( 1 S 2 2 S 1 2 ) ( dy dt ) 2 +2( 1+ S 1 S 2 )g S 2 S 1 y=0 2{ ( 1+ S 2 S 1 )h+( 1 S 2 2 S 1 2 )y+ S 2 S d } d 2 y d t 2 +( 1 S 2 2 S 1 2 ) ( dy dt ) 2 +2( 1+ S 2 S 1 )gy=0 { h+ S 1 S 2 S( S 1 + S 2 ) d+( 1 S 2 S 1 )y } d 2 y d t 2 + 1 2 ( 1 S 2 S 1 ) ( dy dt ) 2 +gy=0 (A+By) d 2 y d t 2 + 1 2 B ( dy dt ) 2 +gy=0

Resolveremos la ecuación del movimiento por procedimientos numéricos con la siguientes condiciones iniciales: en el instante t=0, la altura del fluido en el recipiente de la derecha por encima de la de equilibrio es y0=x02=h02-h (véase la primera figura). La velocidad inicial del líquido en los recipientes es cero, dy/dt=0

Oscilaciones armónicas

Cuando las secciones de los recipientes son iguales, S1=S2. La ecuación del movimiento es

{ h+ S 2 2S d } d 2 y d t 2 +gy=0

Se trata de la ecuación diferencial de las oscilaciones libres de frecuencia angular ω0 y periodo P=2π/ω0

ω 0 2 = g h+ S 2 2S d

Energía

Hemos calculado, la energía cinética y potencial del sistema formado por los dos recipientes y el tubo horizontal que los comunica

La energía inicial es solamente potencial

E 0 = 1 2 ρ S 2 ( 1+ S 1 S 2 )g( h 2 + S 2 S 1 x 02 2 )

La energía en el instante t es

E= 1 2 ρ S 2 { ( 1+ S 2 S 1 )h+( 1 S 2 2 S 1 2 ) x 2 + S 2 S d } v 2 2 + 1 2 ρ S 2 ( 1+ S 1 S 2 )g( h 2 + S 2 S 1 x 2 2 )

Si el líquido es ideal, la energía se mantiene constante

{ h+ S 1 S 2 S( S 1 + S 2 ) d+( 1 S 2 S 1 ) x 2 } v 2 2 =g( x 02 2 x 2 2 ) ( A+By ) ( dy dt ) 2 =g( y 0 2 y 2 )

Periodo

Despejamos dt e integramos

dt= A+By g( y 0 2 y 2 ) dy t= B g y 0 y C+x ( y 0 x)( y 0 +x) dx

Donde C=A/B. Buscamos la solución a la integral

b u xc ( ax )(xb) dx=2 ac ·E(θ,k)2 (au)(ub) uc θ=arcsin (ac)(ub) (ab)(uc) ,k= ab ac

E(θ,k), se denomina integral elíptica incompleta de segunda especie. Identificamos: u=y, c=-C=-A/B, a=y0, b=-y0

t=2 B g { y 0 +C ·E(θ,k) ( y 0 y)(y+ y 0 ) y+C } θ=arcsin ( y 0 +C)(y+ y 0 ) 2 y 0 (y+C) ,k= 2 y 0 y 0 +C

El periodo es el doble del tiempo que tarda el nivel del líquido en el recipiente de la derecha, en subir desde h-y0 a h+y0. Cuando el límite superior de la integral y=y0

θ=arcsin(1)= π 2 ,k= 2 y 0 y 0 +C P=2 B g y 0 y 0 C+x ( y 0 x)( y 0 +x) dx=4 B g y 0 +C ·E(k)

E(k), es una integral elíptica completa

Solución numérica

Consideremos el sistema formado por dos vasos comunicantes

En el instante inicial t=0, los niveles del líquido en cada uno de los recipientes es

Calculamos la altura h de equilibrio

h= r 1 2 h 1 + r 2 2 h 2 r 1 2 + r 2 2

Se deberá cumplir que 2h-h2≥0. en caso contrario se modifica h2

Calculamos el periodo P y resolvemos la ecuación diferencial por el procedimiento numérico ode45 de MATLAB. Las condiciones iniciales son: t=0, y0=x02=h2-h, en reposo, dy/dt=0

El nivel del líquido en el recipiente de la derecha en el instante t es h2=h+x2=h+y

r2=5; % fijo
r1=10; %mayor que r2
h2=20; %mayor que h1
h1=10; %fijo
d=10; %tubo que comunica los recipientes
rd=0.2; 

%altura de equilibrio
h=(r1^2*h1+r2^2*h2)/(r1^2+r2^2);
if 2*h-h2<0
    disp('Diminuir la altura inicial')
    return;
end
A=h+r1^2*r2^2*d/(rd^2*(r1^2+r2^2));
B=1-r2^2/r1^2;
%periodo
if r1==r2
    P=2*pi*sqrt((h+r2^2*d/(2*rd^2))/9.8);
else
    k2=2*(h2-h)/(h2-h+A/B);
    [K,E]=ellipke(k2);
   P=4*sqrt(B/9.8)*sqrt(h2-h+A/B)*E;
end

f=@(t,x) [x(2);(-9.8*x(1)-B*x(2)^2/2)/(A+B*x(1))]; 
[t,x]=ode45(f,[0,P],[h2-h,0]);
plot(t,x(:,1)+h)
line([0,P],[h,h],'lineStyle','--','color','k')
grid on
xlabel('t (s)')
ylabel('h_2 (cm)');
title('Vasos comunicantes ')

La línea a trazos señala la altura h de equilibrio. El periodo P en s, vale

P =  142.0929

Comprobamos que la energía se mantiene constante.

( A+By ) ( dy dt ) 2 =g( y 0 2 y 2 )

En el lado izquierdo, tenemos el cambio de energía cinética, en el derecho, el cambio de energía potencial, la suma debe anularse, ΔEk+ΔEp=0

>> (A+B*x(:,1)).*x(:,2).^2-9.8*((h2-h)^2-x(:,1).^2)
ans =
         0
    0.0000
   -0.0000
   ...
   -0.1115
   -0.1116

Solución analítica

La ecuación

t=2 B g { y 0 +C ·E(θ,k) ( y 0 y)(y+ y 0 ) y+C } θ=arcsin ( y 0 +C)(y+ y 0 ) 2 y 0 (y+C) ,k= 2 y 0 y 0 +C

nos proporciona la ecuación implícita t=f(y). Damos valores a y en el intervalo -y0yy0, obteniendo el nivel de líquido en el recipiente derecho h+y en función del tiempo t durante medio periodo de la oscilación

Comprobamos que la solución analítica coincide con la solución utilizando procedimientos numéricos

r2=5; %cm fijo
r1=10; %mayor que r2
h2=20; %mayor que h1
h1=10; %fijo
d=10; %tubo que comunica los recipientes
rd=0.2; 

%altura de equilibrio
h=(r1^2*h1+r2^2*h2)/(r1^2+r2^2);
if 2*h-h2<0
    disp('Diminuir la altura inicial')
    return;
end
A=h+r1^2*r2^2*d/(rd^2*(r1^2+r2^2));
B=1-r2^2/r1^2;

%periodo
if r1==r2 %armónico
    P=2*pi*sqrt((h+r2^2*d/(2*rd^2))/9.8);
else
    k2=h2/(h2+A/B);
    [K,E]=ellipke(k2);
    P=4*sqrt(B/9.8)*sqrt(h2+A/B)*E;
end

%solución analítica
if r1==r2 %armónico
    w0=sqrt(9.8/(h+r2^2*d/(2*rd^2)));
    t=linspace(0,pi/w0,100);
    y=-(h2-h)*cos(w0*t);
else
    y=linspace(-(h2-h),(h2-h),100);
    k2=2*(h2-h)/(h2-h+A/B);
    th=asin(sqrt((h2-h+A/B)*(y+h2-h)./(2*(h2-h)*(y+A/B))));
    E=ellipticE(th,k2);
    t=2*sqrt(B/9.8)*(sqrt((h2-h)+A/B)*E-sqrt((h2-h-y).*(y+h2-h)./(y+A/B)));
end

hold on
plot(t,y+h)
plot(t+P/2,-y+h)
line([0,P],[h,h],'lineStyle','--','color','k')

%solución numérica
f=@(t,x) [x(2);(-9.8*x(1)-B*x(2)^2/2)/(A+B*x(1))]; 
[t,x]=ode45(f,[0,P],[-h2+h,0]);
plot(t,x(:,1)+h)

hold off
grid on
xlabel('t (s)')
ylabel('h_2 (cm)');
title('Vasos comunicantes ')

Actividades

Se introduce

En el instante inicial t=0, los niveles del líquido en cada uno de los recipientes son

Calculamos la altura h de equilibrio. Si 2h-h2<0, un mensaje nos invita a modificar la altura h2 o el radio r1


Referencias

I.S. Gradshteyn, I.M. Ryzhik. Table of Integrals, Series, and Products. Seventh Edition. Elsevier (2007). 3.141, n° 16, pág. 267,

A.R. Kazachkov, V.A. Lykah, K.A. Minakova, E.S. Syrkin, O.Y. Tkachenko. Liquid nonlinear oscillatios in the U-tube system. Proceedings of the 5th International Conference on Nonlinear Dynamics, ND-KhPI2016. September 27-30, 2016, Kharkov, Ukraine