Consideremos un puente líquido entre dos placas planas paralelas separadas una distancia h. Las placas están hechas del mismo material para que el ángulo de contacto θ_c sea el mismo y la distancia del eje Z al punto de contacto x_c sea también el mismo

En la figura, se representa el perfil z=z(x) que al girar alrededor del eje Z genera la supericie de revolución que describe el puente líquido

Por simetría, solamente tenemos que describir la función z=z(x) entre x_m y x_c

Calculamos el área de la superficie de revolución

${\displaystyle dA=2\pi x·ds=2\pi x\sqrt{1+{\left(\frac{dz}{dx}\right)}^2}dx}$

La energía superficial (multiplicamos por dos, debido a la simetría) es

${\displaystyle E_s=2\gamma {\displaystyle \underset{x_m}{\overset{x_c}{\int }}2\pi x·\sqrt{1+{\left(\frac{dz}{dx}\right)}^2}dx}=4\gamma \pi {\displaystyle \underset{x_m}{\overset{x_c}{\int }}x·\sqrt{1+{\left(\frac{dz}{dx}\right)}^2}dx}}$

Volumen de líquido

${\displaystyle dV=\pi x^{2}dz=\pi x^2\frac{dz}{dx}dx}$

Por simetría, el volumen entre x_m y x_c se multiplica por dos

${\displaystyle V=2\pi {\displaystyle \underset{x_m}{\overset{x_c}{\int }}{x}^2\frac{dz}{dx}dx}}$

Tenemos que hacer mínima la energía superficial E_s manteniendo constante el volumen V

Vamos a resolver este problema de extremo condicionado. Creamos la función auxiliar F, dependiente de un parámetro λ

${\displaystyle F\left(x,z,\dot{z}\right)=4\gamma \pi x\sqrt{1+{\dot{z}}^2}+2\pi \lambda x^2\dot{z},\dot{z}=\frac{dz}{dx}}$

Aplicamos la ecuación Euler-Lagrange

${\displaystyle \frac{\partial F}{\partial z}-\frac{d}{dx}\left(\frac{\partial F}{\partial \dot{z}}\right)=0}$

La función F no depende de z

${\displaystyle \begin{array}{l}-\frac{d}{dx}\left(\frac{4\pi \gamma x\dot{z}}{\sqrt{1+{\dot{z}}^2}}+2\pi \lambda x^2\right)=0\\ \frac{Rx\frac{dz}{dx}}{\sqrt{1+{\left(\frac{dz}{dx}\right)}^2}}+x^2=C,R=\frac{2\gamma }{\lambda }\end{array}}$

donde C es una constante.

${\displaystyle \begin{array}{l}{\frac{dz}{dx}|}_{x=x_m}=\infty ,R{x}_m+x_m^2=C\\ {\frac{dz}{dx}|}_{x=x_c}=\tan {\theta }_c,R{x}_c\sin {\theta }_c+x_c^2=C\end{array}}$

Eliminamos C de este sistema de dos ecuaciones y despejamos R

${\displaystyle R=\frac{x_c^2-x_m^2}{x_m-x_c\sin {\theta }_c}}$

Como vemos en la figura x_c>x_m. El parámetro R puede ser positivo, negativo o infinito.

Elevamos al cuadrado y despejamos dz/dx

${\displaystyle \begin{array}{l}\frac{{\left(Rx\frac{dz}{dx}\right)}^2}{1+{\left(\frac{dz}{dx}\right)}^2}={\left(C-x^2\right)}^2\\ {\left(\frac{dz}{dx}\right)}^2=\frac{{\left(C-x^2\right)}^2}{R^2{x}^2-{\left(C-x^2\right)}^2}\\ \frac{dz}{dx}=\pm \frac{C-x^2}{\sqrt{R^2{x}^2-{\left(C-x^2\right)}^2}}\\ \frac{dz}{dx}=\pm \frac{R{x}_m+x_m^2-x^2}{\sqrt{R^2{x}^2-{\left(R{x}_m+x_m^2-x^2\right)}^2}}\end{array}}$

Como vemos en la figura, en el intervalo x_m≤x≤x_c, la pendiente dz/dx es positiva.

${\displaystyle {\frac{dz}{dx}|}_{x=x_c}=\pm \frac{R{x}_m+x_m^2-x_c^2}{\sqrt{R^2{x}_c^2-{\left(R{x}_m+x_m^2-x_c^2\right)}^2}}=\pm \frac{R{x}_c\sin {\theta }_c}{\sqrt{R^2{x}_c^2-{\left(R{x}_m+x_m^2-x_c^2\right)}^2}}>0}$

• Si R>0, entonces se toma el signo +
• Si R<0, entonces se toma el signo -

En esta página, estudiamos el primer caso

${\displaystyle \frac{dz}{dx}=\frac{R{x}_m+x_m^2-x^2}{\sqrt{R^2{x}^2-{\left(R{x}_m+x_m^2-x_c^2\right)}^2}}}$

Se define el parámetro k

${\displaystyle \begin{array}{l}k=\frac{R+x_m}{x_m}>1\\ \frac{dz}{dx}=\frac{\left(k-1\right)x_m^2+x_m^2-x^2}{\sqrt{{\left(k-1\right)}^2{x}_m^2{x}^2-{\left(\left(k-1\right)x_m^2+x_m^2-x^2\right)}^2}}=\frac{k{x}_m^2-x^2}{\sqrt{{\left(k-1\right)}^2{x}_m^2{x}^2-{\left(k{x}_m^2-x^2\right)}^2}}\\ \frac{dz}{dx}=\frac{k{x}_m^2-x^2}{\sqrt{k^2{x}_m^2{x}^2+x_m^2{x}^2-k^2{x}_m^4-x^4}}\end{array}}$

Hacemos el cambio de variable

${\displaystyle \begin{array}{l}x=k{x}_{m}u,dx=k{x}_m·du\\ \frac{dz}{k{x}_m·du}=\frac{k{x}_m^2-k^2{x}_m^2{u}^2}{\sqrt{k^4{x}_m^4{u}^2+x_m^4{k}^2{u}^2-k^2{x}_m^4-k^4{x}_m^4{u}^4}}\\ \frac{dz}{du}=\frac{k{x}_m\left(1-k{u}^2\right)}{\sqrt{\left(k^2{u}^2-1\right)\left(1-u^2\right)}}=\frac{x_m\left(1-k{u}^2\right)}{\sqrt{\left(u^2-\frac{1}{k^2}\right)\left(1-u^2\right)}}\end{array}}$

Solución numérica

El perfil del puente líquido viene dado por la función z=z(x)

${\displaystyle z=x_m{\displaystyle \underset{1/k}{\overset{\frac{x}{k{x}_m}}{\int }}\frac{1-k{u}^2}{\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}}du},k=\frac{{\left(\frac{x_c}{x_m}\right)}^2-1}{1-\left(\frac{x_c}{x_m}\right)\sin {\theta }_c}+1>1}$

Para que k>1, el cociente (x_c/x_m)<1/sinθ_c

Para x=x_c, z=h, la mitad de la separación entre las placas

${\displaystyle h=x_m{\displaystyle \underset{1/k}{\overset{\frac{x_c}{k{x}_m}}{\int }}\frac{1-k{u}^2}{\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}}du}}$

Volumen de líquido

${\displaystyle V=2\pi {\displaystyle \underset{x_m}{\overset{x_c}{\int }}{x}^2\frac{k{x}_m^2-x^2}{\sqrt{k^2{x}_m^2{x}^2+x_m^2{x}^2-k^2{x}_m^4-x^4}}dx}}$

Hacemos el cambio de variable

${\displaystyle \begin{array}{l}x=k{x}_{m}u,dx=k{x}_{m}du\\ V=2\pi {\displaystyle \underset{\frac{1}{k}}{\overset{\frac{x_c}{k{x}_m}}{\int }}{k}^2{x}_m^2{u}^2\frac{k{x}_m^2-k^2{x}_m^2{u}^2}{\sqrt{k^4{x}_m^4{u}^2+x_m^4{k}^2{u}^2-k^2{x}_m^4-k^4{x}_m^4{u}^4}}k{x}_{m}du}\\ V=2\pi k^3{x}_m^3{\displaystyle \underset{\frac{1}{k}}{\overset{\frac{x_c}{k{x}_m}}{\int }}\frac{u^2\left(1-k{u}^2\right)}{\sqrt{k^2{u}^2+u^2-1-k^2{u}^4}}du}\\ V=2\pi k^2{x}_m^3{\displaystyle \underset{\frac{1}{k}}{\overset{\frac{x_c}{k{x}_m}}{\int }}\frac{u^2\left(1-k{u}^2\right)}{\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}}du}\end{array}}$

El cociente h³/V es una función de ξ=x_c/x_m

${\displaystyle \begin{array}{l}r=\frac{h}{V^{\frac{1}{3}}}\\ h^3\left(\xi \right)-r^{3}V\left(\xi \right)=\mathrm{0,}\xi =\frac{x_c}{x_m}>1\end{array}}$

Conocido el ángulo de contacto θ_c. Fijamos el valor $r=\frac{h}{V^{1/3}}$, determinamos ξ=x_c/x_m resolviendo la ecuación transcendente

Fijado el valor del volumen V=1, determinamos x_m

Resultados

• Angulo de contacto, θ_c=π/6 (30°)
• El parámtro, r=0.5

th=pi/6; %ángulo de contacto
kk=@(x) 1+(x^2-1)/(1-x*sin(th)); %x=xc/xm
n=1;
rel=0.5; %relación h^3/V
xx=linspace(1.1,1/sin(th)-0.01,50);
h=zeros(1, length(xx));
V=zeros(1, length(xx));
for x=xx
    k=kk(x);
    f=@(u) (1-k*u.^2)./sqrt((1-u.^2).*(u.^2-1/k^2));
    h(n)=integral(f,1/k,x/k);
    g=@(u) (1-k*u.^2).*u.^2./sqrt((1-u.^2).*(u.^2-1/k^2));
    V(n)=2*pi*(k^2)*integral(g,1/k,x/k);
    n=n+1;
end
yy=rel^3*V-h.^3;
plot(xx,yy) 
grid on
xlabel('\xi')
ylabel('f(\xi)')
title('Raíz')

La función transcendente tiene una raíz proxima a ξ=1.82. Creamos una función raiz que nos aproxime a la raíz mediante una interpolación lineal

La ecuación de la recta que pasa por dos puntos es

${\displaystyle y=\frac{y_2-y_1}{x_2-x_1}x+\frac{y_1{x}_2-y_2{x}_1}{x_2-x_1}}$

La abcisa x que hace y=0 es

${\displaystyle x=\frac{y_2{x}_1-y_1{x}_2}{y_2-y_1}}$

   function r = raiz(xx, yy)
         indice=find(yy(1:end-1).*yy(2:end)<0);
        x1=xx(indice); y1=yy(indice);
        x2=xx(indice+1); y2=yy(indice+1);
        r=(x1*y2-x2*y1)/(y2-y1); %interpolación lineal
     end

La primera tarea es buscar el intervalo en el que la función cambia de signo, después calcular la raíz mediante interpolación lineal

Cuando esperamos más de una raíz, el código de la función raiz es el siguiente

 function r = raiz(xx, yy)
        indices=find(yy(1:end-1).*yy(2:end)<0);
        r=zeros(1,length(indices));
        for j=1:length(indices)
            x1=xx(indices(j)); y1=yy(indices(j));
            x2=xx(indices(j)+1); y2=yy(indices(j)+1);
            r(j)=(x1*y2-x2*y1)/(y2-y1); %interpolación lineal
        end
    end

Una vez que tenemos la raíz ξ=x_c/x_m, calculamos x_m sabiendo que el volumen V=1, después x_c=ξ·x_m y representamos el perfil z=z(x). El código completo es

function pelicula_2
    th=pi/6; %ángulo de contacto
    kk=@(x) 1+(x^2-1)/(1-x*sin(th)); %x=xc/xm
    n=1;
    rel=0.5;
    xx=linspace(1.1,1/sin(th)-0.01,50);
    h=zeros(1, length(xx));
    V=zeros(1, length(xx));
    for x=xx
        k=kk(x);
        f=@(u) (1-k*u.^2)./sqrt((1-u.^2).*(u.^2-1/k^2));
        h(n)=integral(f,1/k,x/k);
        g=@(u) (1-k*u.^2).*u.^2./sqrt((1-u.^2).*(u.^2-1/k^2));
        V(n)=2*pi*(k^2)*integral(g,1/k,x/k);
        n=n+1;
    end
    yy=rel^3*V-h.^3;
%     plot(xx,yy) %para probar si hay raíces
    rr=raiz(xx,yy);
    xr=rr(1); %primera raíz
    k=kk(xr);
    g=@(u) (1-k*u.^2).*u.^2./sqrt((1-u.^2).*(u.^2-1/k^2));
    Vr=2*pi*(k^2)*integral(g,1/k,xr/k);
    xm=1/Vr^(1/3);
    xc=xm*xr;
    disp([xm,xc])
    f=@(u) (1-k*u.^2)./sqrt((1-u.^2).*(u.^2-1/k^2));
    xx=linspace(xm+0.01,xc-0.01,50);
    z=zeros(1, length(xx));
    n=1;
    for x=xx
        z(n)=xm*integral(f,1/k,x/(xm*k));
        n=n+1;
    end
     hold on
     h1=z(end);
    fill([xm,xx,xc,-xc,-fliplr(xx),-xm, -xx,-xc,xc,fliplr(xx)],
[0,z,h1,h1,fliplr(z),0,-z,-h1,-h1,-flip(z)],'c')
    plot([xm,xx],[0,z],'b')
    plot([xm,xx],[0,-z],'b')
    plot([-xm,-xx],[0,z],'b')
    plot([-xm,-xx],[0,-z],'b')
    line([-xc,xc],[z(end),z(end)],'lineWidth',1.5,'color','k')
    line([-xc,xc],[-z(end),-z(end)],'lineWidth',1.5,'color','k')
    hold off
    axis equal
    grid on
    xlabel('x')
    ylabel('z')
    title('Perfil')

    function r = raiz(xx, yy)
        indices=find(yy(1:end-1).*yy(2:end)<0);
        r=zeros(1,length(indices));
        for j=1:length(indices)
            x1=xx(indices(j)); y1=yy(indices(j));
            x2=xx(indices(j)+1); y2=yy(indices(j)+1);
            r(j)=(x1*y2-x2*y1)/(y2-y1); %interpolación lineal
        end
    end
end

Los valores x_m=0.4421 y x_c=0.8037

    0.4421    0.8037

Girando alrededor del eje de simetría Z, obtenemos una superficie de revolución, que es el perfil del puente líquido entre dos planos paralelos en asusencia de gravedad

Dos configuraciones

Para ángulos de contacto más pequeños, por ejemplo θ_c=8°, encontramos dos raíces

• Para la primera raíz ξ=2.9121, obtenemos el siguiente perfil



Los valores x_m=0.3615 y x_c=1.0526

    0.3615    1.0526

• Para la segunda raíz ξ=5.9568, obtenemos el siguiente perfil



Los valores x_m=0.2102 y x_c=1.2523

    0.2102    1.2523

Relación volumen/separación

Podemos cambiar el factor $r=\frac{h}{V^{1/3}}$, por ejemplo r=0.2 manteniendo el ángulo de contacto θ_c=π/6

Los valores x_m=0.8573 y x_c=0.9764

0.8573    0.9764

Cuando r>0.53 o r<0.17 aproximadamente, no se encuentran raíces en el intervalo especificado

Solución analítica

Obtenemos la función z=z(x), integrando desde x_m a x o bien, desde 1/k a u=x/(kx_m)

${\displaystyle \begin{array}{l}z=x_m{\displaystyle \underset{1/k}{\overset{u}{\int }}\frac{1-k{u}^2}{\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}}du}=\\ x_m\left\{{\displaystyle \underset{1/k}{\overset{u}{\int }}\frac{du}{\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}}-k{\displaystyle \underset{1/k}{\overset{u}{\int }}\frac{u^2}{\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}}du}}\right\}\end{array}}$

Buscamos el resultado en una tabla de integrales

${\displaystyle \begin{array}{l}{\displaystyle \underset{b}{\overset{u}{\int }}\frac{dx}{\sqrt{\left(a^2-x^2\right)\left(x^2-b^2\right)}}}=\frac{1}{a}F\left(\kappa ,q\right),a\ge u>b>0\\ {\displaystyle \underset{b}{\overset{u}{\int }}\frac{x^{2}dx}{\sqrt{\left(a^2-x^2\right)\left(x^2-b^2\right)}}=}aE\left(\kappa ,q\right)-\frac{1}{u}\sqrt{\left(a^2-u^2\right)\left(u^2-b^2\right)},a\ge u>b>0\\ q=\frac{\sqrt{a^2-b^2}}{a},\kappa =\arcsin \frac{a}{u}\sqrt{\frac{u^2-b^2}{a^2-b^2}}\end{array}}$

El resultado es

${\displaystyle \begin{array}{l}\{\begin{array}{l}a=1\\ b=\frac{1}{k}<1\end{array}\\ z=x_m\left\{F\left(\kappa ,q\right)-k\left(E\left(\kappa ,q\right)-\frac{1}{u}\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}\right)\right\}\\ k=\frac{{\left(\frac{x_c}{x_m}\right)}^2-1}{1-\left(\frac{x_c}{x_m}\right)\sin {\theta }_c}+\mathrm{1,}q=\frac{1}{k}\sqrt{k^2-1},\\ \kappa =\arcsin \frac{1}{u}\sqrt{\frac{k^2{u}^2-1}{k^2-1}}=\arcsin \frac{1}{\frac{x}{k{x}_m}}\sqrt{\frac{k^2{\left(\frac{x}{k{x}_m}\right)}^2-1}{k^2-1}}=\arcsin k\sqrt{\frac{1-{\left(\frac{x_m}{x}\right)}^2}{k^2-1}}\\ z=x_m\left\{F\left(\kappa ,q\right)-k\left(E\left(\kappa ,q\right)-\sqrt{\left(1-\frac{1}{k^2}{\left(\frac{x}{x_m}\right)}^2\right)\left(1-{\left(\frac{x_m}{x}\right)}^2\right)}\right)\right\}\end{array}}$

Comprobamos que z(x_m)=0. Véase Integrales elípticas

La separación entre las placas es 2h. Tenemos que h=z(x_c)

${\displaystyle \begin{array}{l}h=x_m\left\{F\left(\kappa ,q\right)-k\left(E\left(\kappa ,q\right)-\sqrt{\left(1-\frac{1}{k^2}{\left(\frac{x_c}{x_m}\right)}^2\right)\left(1-{\left(\frac{x_m}{x_c}\right)}^2\right)}\right)\right\}\\ q=\frac{1}{k}\sqrt{k^2-1},\kappa =\arcsin k\sqrt{\frac{1-{\left(\frac{x_m}{x_c}\right)}^2}{k^2-1}}\end{array}}$

Volumen de líquido

En el apartado anterior, llegamos a la expresión

${\displaystyle \begin{array}{l}V=2\pi k^2{x}_m^3{\displaystyle \underset{\frac{1}{k}}{\overset{\frac{x_c}{k{x}_m}}{\int }}\frac{u^2\left(1-k{u}^2\right)}{\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}}du}=\\ 2\pi k^2{x}_m^3\left\{{\displaystyle \underset{\frac{1}{k}}{\overset{\frac{x_c}{k{x}_m}}{\int }}\frac{u^2}{\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}}du}-k{\displaystyle \underset{\frac{1}{k}}{\overset{\frac{x_c}{k{x}_m}}{\int }}\frac{u^4}{\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}}du}\right\}\end{array}}$

Buscamos el resultado en una tabla de integrales

${\displaystyle \begin{array}{l}{\displaystyle \underset{b}{\overset{u}{\int }}\frac{x^{4}dx}{\sqrt{\left(1-x^2\right)\left(x^2-b^2\right)}}=\frac{a}{3}}\left\{2\left(a^2+b^2\right)E\left(\kappa ,q\right)-b^{2}F\left(\kappa ,q\right)\right\}-\frac{u^2+2a^2+2b^2}{3u}\sqrt{\left(a^2-u^2\right)\left(u^2-b^2\right)},a\ge u>b>0\\ q=\frac{\sqrt{a^2-b^2}}{a},\kappa =\arcsin \frac{a}{u}\sqrt{\frac{u^2-b^2}{a^2-b^2}}\end{array}}$

La integral vale

${\displaystyle \begin{array}{l}\{\begin{array}{l}a=1\\ b=\frac{1}{k}<1\end{array}\\ \frac{1}{3}\left\{2\left(1+\frac{1}{k^2}\right)E\left(\kappa ,q\right)-\frac{1}{k^2}F\left(\kappa ,q\right)\right\}-\frac{u^2+2+2\frac{1}{k^2}}{3u}\sqrt{\left(1-u^2\right)\left(u^2-\frac{1}{k^2}\right)}\\ \frac{1}{3}\left\{2\left(1+\frac{1}{k^2}\right)E\left(\kappa ,q\right)-\frac{1}{k^2}F\left(\kappa ,q\right)\right\}-\frac{{\left(\frac{x_c}{k{x}_m}\right)}^2+2+2\frac{1}{k^2}}{3\frac{x_c}{k{x}_m}}\sqrt{\left(1-{\left(\frac{x_c}{k{x}_m}\right)}^2\right)\left({\left(\frac{x_c}{k{x}_m}\right)}^2-\frac{1}{k^2}\right)}\\ \frac{1}{3}\left\{2\left(1+\frac{1}{k^2}\right)E\left(\kappa ,q\right)-\frac{1}{k^2}F\left(\kappa ,q\right)\right\}-\frac{{\left(\frac{x_c}{x_m}\right)}^2+2k^2+2}{3k^2}\sqrt{\left(1-\frac{1}{k^2}{\left(\frac{x_c}{x_m}\right)}^2\right)\left(1-{\left(\frac{x_m}{x_v}\right)}^2\right)}\end{array}}$

El volumen V es

${\displaystyle \begin{array}{l}V=2\pi k^2{x}_m^3\left\{\begin{array}{l}E\left(\kappa ,q\right)-\sqrt{\left(1-\frac{1}{k^2}{\left(\frac{x_c}{x_m}\right)}^2\right)\left(1-{\left(\frac{x_m}{x_c}\right)}^2\right)}\\ -\frac{1}{3}\left\{2k\left(1+\frac{1}{k^2}\right)E\left(\kappa ,q\right)-\frac{1}{k}F\left(\kappa ,q\right)\right\}+\frac{{\left(\frac{x_c}{x_m}\right)}^2+2k^2+2}{3k}\sqrt{\left(1-\frac{1}{k^2}{\left(\frac{x_c}{x_m}\right)}^2\right)\left(1-{\left(\frac{x_m}{x_c}\right)}^2\right)}\end{array}\right\}\\ k=\frac{{\left(\frac{x_c}{x_m}\right)}^2-1}{1-\left(\frac{x_c}{x_m}\right)\sin {\theta }_c}+\mathrm{1,}q=\frac{1}{k}\sqrt{k^2-1},\kappa =\arcsin k\sqrt{\frac{1-{\left(\frac{x_m}{x_c}\right)}^2}{k^2-1}}\end{array}}$

Conocido el ángulo de contacto θ_c. Fijamos el valor $r=\frac{h}{V^{1/3}}$, determinamos ξ=x_c/x_m resolviendo la ecuación transcendente

Fijado el valor del volumen V=1, determinamos x_m

Resultados

• Angulo de contacto, θ_c=π/6 (30°)
• El parámtro, r=0.5

th=pi/6; %ángulo de contacto
k=@(x) 1+(x.^2-1)./(1-x*sin(th));
q2=@(x) (k(x).^2-1)./k(x).^2;
K=@(x) asin(k(x).*sqrt((1-1./x.^2)./(k(x).^2-1)));
rel=0.5;
h=@(x) ellipticF(K(x),q2(x))-k(x).*(ellipticE(K(x),q2(x))-
sqrt((1-(x.^2)./k(x).^2).*(1-1./x.^2)));
V=@(x) (2*pi*k(x).^2).*(ellipticE(K(x),q2(x))-sqrt((1-(x.^2).
/k(x).^2).*(1-1./x.^2))-(2*k(x).*(1+1./k(x).^2).*ellipticE(K(x),q2(x))
-ellipticF(K(x),q2(x))./k(x))/3+((x.^2+2*k(x).^2+2).*sqrt((1-(x.^2)./k(x).^2)
.*(1-1./x.^2)))./(3*k(x)));
f=@(x) rel^3*V(x)-h(x).^3; 
fplot(f,[1.05,1/sin(th)])
xr=fzero(f, [1.05,1/sin(th)-0.01]); %cociente xc/xm
grid on
xlabel('\xi')
ylabel('f(\xi)')
title('Raíz')

La raíz buscada es ξ=1.8182, coincide aproximadamente, con la obtenida mediante procedimientos numéricos

>> xr
xr =    1.8182

Representamos la función z=z(x) entre x_m y x_c

th=pi/6; %ángulo de contacto
k=@(x) 1+(x.^2-1)./(1-x*sin(th));
q2=@(x) (k(x).^2-1)./k(x).^2;
K=@(x) asin(k(x).*sqrt((1-1./x.^2)./(k(x).^2-1)));
rel=0.5;
h=@(x) ellipticF(K(x),q2(x))-k(x).*(ellipticE(K(x),q2(x))-
sqrt((1-(x.^2)./k(x).^2).*(1-1./x.^2)));
V=@(x) (2*pi*k(x).^2).*(ellipticE(K(x),q2(x))-sqrt((1-(x.^2)./k(x).^2).
*(1-1./x.^2))-(2*k(x).*(1+1./k(x).^2).*ellipticE(K(x),q2(x))-
ellipticF(K(x),q2(x))./k(x))/3+((x.^2+2*k(x).^2+2).*sqrt((1-(x.^2)./k(x).^2).
*(1-1./x.^2)))./(3*k(x)));
f=@(x) rel^3*V(x)-h(x).^3; 
% fplot(f,[1.05,1/sin(th)])
xr=fzero(f, [1.05,1/sin(th)-0.01]); %cociente xc/xm
xm=1/V(xr)^(1/3);
xc=xm*xr;
kr=k(xr);
qr2=(kr^2-1)/kr^2;
K=@(x) asin(kr*sqrt((1-(xm./x).^2)/(kr^2-1)));
z=@(x) xm*(ellipticF(K(x),qr2)-kr.*(ellipticE(K(x),qr2)-sqrt((1-(x/xm).^2/kr^2).
*(1-(xm./x).^2))));

hold on
fplot(z, [xm,xc])
hold off
axis equal
grid on
xlabel('x')
ylabel('z')
title('Perfil')

Los valores x_m=0.4421 y x_c=0.8037

>> xm,xc
xm =    0.4421
xc =    0.8037

th=pi/6; %ángulo de contacto
k=@(x) 1+(x.^2-1)./(1-x*sin(th));
q2=@(x) (k(x).^2-1)./k(x).^2;
K=@(x) asin(k(x).*sqrt((1-1./x.^2)./(k(x).^2-1)));
rel=0.5;
h=@(x) ellipticF(K(x),q2(x))-k(x).*(ellipticE(K(x),q2(x))-sqrt((1-(x.^2).
/k(x).^2).*(1-1./x.^2)));
V=@(x) (2*pi*k(x).^2).*(ellipticE(K(x),q2(x))-sqrt((1-(x.^2)./k(x).^2).*
(1-1./x.^2))-(2*k(x).*(1+1./k(x).^2).*ellipticE(K(x),q2(x))-ellipticF(K(x),q2(x)).
/k(x))/3+((x.^2+2*k(x).^2+2).*sqrt((1-(x.^2)./k(x).^2).*(1-1./x.^2)))./(3*k(x)));
f=@(x) rel^3*V(x)-h(x).^3; 
% fplot(f,[1.05,1/sin(th)])
xr=fzero(f, [1.05,1/sin(th)-0.01]); %cociente xc/xm
xm=1/V(xr)^(1/3);
xc=xm*xr;
kr=k(xr);
qr2=(kr^2-1)/kr^2;
K=@(x) asin(kr*sqrt((1-(xm./x).^2)/(kr^2-1)));
z=@(x) xm*(ellipticF(K(x),qr2)-kr.*(ellipticE(K(x),qr2)-sqrt((1-(x/xm).^2/kr^2).*
(1-(xm./x).^2))));

 hold on
 h1=z(xc);
 fp=fplot(z, [xm,xc],'b');
  fill([xm,fp.XData,xc,-xc,-fliplr(fp.XData),-xm, -fp.XData,-xc,xc,
fliplr(fp.XData)],
[0,fp.YData,h1,h1,fliplr(fp.YData),0,-fp.YData,-h1,-h1,-flip(fp.YData)],'c')
 plot(fp.XData,fp.YData,'b')
 plot(fp.XData,-fp.YData,'b')
 plot(-fp.XData,fp.YData,'b')
  plot(-fp.XData,-fp.YData,'b')
hold off
axis equal
grid on
xlabel('x')
ylabel('z')
title('Perfil')

Girando alrededor del eje de simetría Z, obtenemos una superficie de revolución, que es el perfil del puente líquido entre dos planos paralelos en asusencia de gravedad

Representación en tres dimensiones del puente líquido

th=pi/6; %ángulo de contacto
k=@(x) 1+(x.^2-1)./(1-x*sin(th));
q2=@(x) (k(x).^2-1)./k(x).^2;
K=@(x) asin(k(x).*sqrt((1-1./x.^2)./(k(x).^2-1)));
rel=0.5;
h=@(x) ellipticF(K(x),q2(x))-k(x).*(ellipticE(K(x),q2(x))-
sqrt((1-(x.^2)./k(x).^2).*(1-1./x.^2)));
V=@(x) (2*pi*k(x).^2).*(ellipticE(K(x),q2(x))-sqrt((1-(x.^2)./k(x).^2).
*(1-1./x.^2))-(2*k(x).*(1+1./k(x).^2).*ellipticE(K(x),q2(x))-
ellipticF(K(x),q2(x))./k(x))/3+((x.^2+2*k(x).^2+2).*sqrt((1-(x.^2)./k(x).^2).
*(1-1./x.^2)))./(3*k(x)));
f=@(x) rel^3*V(x)-h(x).^3; 
% fplot(f,[1.05,1/sin(th)])
xr=fzero(f, [1.05,1/sin(th)-0.01]); %cociente xc/xm
xm=1/V(xr)^(1/3);
xc=xm*xr;
kr=k(xr);
qr2=(kr^2-1)/kr^2;
K=@(x) asin(kr*sqrt((1-(xm./x).^2)/(kr^2-1)));
z=@(x) xm*(ellipticF(K(x),qr2)-kr.*(ellipticE(K(x),qr2)-
sqrt((1-(x/xm).^2/kr^2).*(1-(xm./x).^2))));

r=linspace(xm,xc,50);
phi=linspace(0,2*pi,30);
[R,Phi]=meshgrid(r,phi);
X=R.*cos(Phi);
Y=R.*sin(Phi);
Z=z(R);
hold on
surfl(X,Y,Z)
surfl(X,Y,-Z)
shading interp
colormap(gray);
hold off
grid on
xlabel('x')
ylabel('y')
zlabel('z')
title('Perfil')
view(105,22)

Dos configuraciones

Para ángulos de contacto más pequeños, por ejemplo θ_c=8°, encontramos dos raíces

La primera raíz está en el intervalo [1,4] y la segunda en el intervalo [5,7]. Calculamos el valor de las raíces mediante el procedimiento fzero de MATLAB en los intervalos apropiados

• Para la primera raíz ξ= 2.9123, obtenemos el siguiente perfil



Los valores x_m=0.3614 y x_c=1.0526

>> xm,xc
xm =    0.3614
xc =    1.0526

• Para la segunda raíz ξ= 5.9578, obtenemos el siguiente perfil



Los valores x_m=0.2102 y x_c=1.2524

>> xm,xc
xm =    0.2102
xc =    1.2524

Para decidir cuál es la configuración real, calculamos la energía superficial

${\displaystyle \begin{array}{l}{E}_s=4\gamma \pi {\displaystyle \underset{x_m}{\overset{x_c}{\int }}x·\sqrt{1+{\left(\frac{dz}{dx}\right)}^2}dx}\\ \frac{dz}{dx}=\frac{dz}{k{x}_m·du}=\frac{1}{k}\frac{1-k{u}^2}{\sqrt{\left(u^2-\frac{1}{k^2}\right)\left(1-u^2\right)}}=\frac{1-k{u}^2}{\sqrt{\left(k^2{u}^2-1\right)\left(1-u^2\right)}}\\ E_s=4\gamma \pi {\displaystyle \underset{\frac{1}{k}}{\overset{\frac{x_c}{k{x}_m}}{\int }}k{x}_{m}u·\sqrt{1+{\left(\frac{1-k{u}^2}{\sqrt{\left(k^2{u}^2-1\right)\left(1-u^2\right)}}\right)}^2}k{x}_{m}du}=\\ 4\gamma \pi k^2\left(k-1\right)x_m^2{\displaystyle \underset{\frac{1}{k}}{\overset{\frac{x_c}{k{x}_m}}{\int }}\frac{u^2}{\sqrt{\left(k^2{u}^2-1\right)\left(1-u^2\right)}}du}=4\gamma \pi k{x}_m^2\left(k-1\right){\displaystyle \underset{\frac{1}{k}}{\overset{\frac{x_c}{k{x}_m}}{\int }}\frac{u^2}{\sqrt{\left(u^2-\frac{1}{k^2}\right)\left(1-u^2\right)}}du}\\ E_s=4\gamma \pi k\left(k-1\right)x_m^2\left\{E\left(\kappa ,q\right)-\sqrt{\left(1-\frac{1}{k^2}{\left(\frac{x_c}{x_m}\right)}^2\right)\left(1-{\left(\frac{x_m}{x_c}\right)}^2\right)}\right\}\\ k=\frac{{\left(\frac{x_c}{x_m}\right)}^2-1}{1-\left(\frac{x_c}{x_m}\right)\sin {\theta }_c}+\mathrm{1,}q=\frac{1}{k}\sqrt{k^2-1},\kappa =\arcsin k\sqrt{\frac{1-{\left(\frac{x_m}{x_c}\right)}^2}{k^2-1}}\end{array}}$

th=8*pi/180; %ángulo de contacto
k=@(x) 1+(x.^2-1)./(1-x*sin(th));
q2=@(x) (k(x).^2-1)./k(x).^2;
K=@(x) asin(k(x).*sqrt((1-1./x.^2)./(k(x).^2-1)));
rel=0.5;
V=@(x) (2*pi*k(x).^2).*(ellipticE(K(x),q2(x))-sqrt((1-(x.^2)./k(x).^2)
.*(1-1./x.^2))-(2*k(x).*(1+1./k(x).^2).*ellipticE(K(x),q2(x))-
ellipticF(K(x),q2(x))./k(x))/3+((x.^2+2*k(x).^2+2).*sqrt((1-(x.^2)./
k(x).^2).*(1-1./x.^2)))./(3*k(x)));
for xr=[2.9123,5.9578]
    xm=1/V(xr)^(1/3);
    xc=xm*xr;
    kr=k(xr);
    K=@(x) asin(kr*sqrt((1-(xm/x)^2)/(kr^2-1)));
    qr2=(kr^2-1)/kr^2;
    Es=kr*(kr-1)*xm^2*(ellipticE(K(xc),qr2)-sqrt((1-(xc/xm)^2/kr^2)*
(1-(xm/xc)^2)));
    disp([xr,Es])
end

    2.9123    0.5936
    5.9578    0.8238

La energía superficial de la segunda configuración es la mayor. La primera configuración ξ=2.9123, tiene la energía superficial mínima

Relación volumen/separación

Podemos cambiar el factor $r=\frac{h}{V^{1/3}}$, por ejemplo r=0.2 manteniendo el ángulo de contacto θ_c=π/6

Los valores x_m=0.8564 y x_c=0.9758

>> xm,xc
xm =    0.8564
xc =    0.9758

Los valores de x_m y x_c coinciden prácticamente con los calculados mediante procedimientos numéricos

Referencias

I. S. Gradshteyn, I. M. Ryzhik. Table of Integrals, Series, and Products.. Eighth Edition. Elsevier. 2015

• pp. 282, 3.152, n° 9
• pp. 283, 3.153, n° 7
• pp. 284, 3.154, n° 7

YU Feng-jun, JU Lin, ZHANG Xi-wei, TANG Zhen-jie, TIAN Jun-long. Theoretical study on space liquid bridge experiment — another teaching case caused by space class. College Physics. 2023, 42 (05):