El perfil de una gota de rocío

En una mañana fría cuando está rociando, observamos gotas del vapor de agua condensado en las hierbas, en las hojas de las plantas, en las telas de araña. Nos fijamos en las pequeñas gotas casi esféricas sobre la superficie de algunas hojas.

Antes de proceder al estudio de este tema complejo recordaremos un concepto geométrico

Area lateral de un cono truncado

Empezamos por el área lateral de un cono de radio r de la base y altura H

Supongamos que el cono está hecho de papel, cortamos el cono a lo largo de su generatriz y lo extendemos sobre una superficie plana, obtenemos un sector circular de radio L (generatriz) y ángulo θ=2πr/L. El área A de este sector circular es

A=πrL=πr H 2 + r 2

Sea un tronco de cono cuyas bases tiene radios b (mayor) y a (menor) y altura h

El área lateral del tronco de cono es la diferencia entre la de un cono de radio base b y altura H y la del cono de radio base a y altura H-h, tal como se muestra en la figura.

A=πb H 2 + b 2 πa ( Hh ) 2 + a 2

En los dos triángulos rectángulos, establecemos las relaciones H=bh/(b-a) y H-h=ah/(b-a). La expresión del área A se simplifica notablemente

A=π(a+b) h 2 + ( ba ) 2

Una gota de líquido sobre un plano horizontal

Cuando una gota de líquido se coloca sobre una superficie horizontal, la fuerza de la gravedad tiende a extenderla para reducir la energía potencial gravitatoria, mientras que la tensión superficial tiende a juntar el líquido para reducir la energía superficial. La forma de equilibrio de la gota será aquella que haga mínimo la suma de las energías superficial y gravitatoria. Más adelante, deduciremos las ecuaciones diferenciales que describen la forma de la gota aplicando el cálculo de variaciones.

La forma de la gota es una superficie de revolución alrededor del eje vertical Y, utilizaremos dos funciones x=x(θ) e y=y(θ) dependientes de un parámetro θ para describirla. Los datos de partida son el volumen V de la gota y la altura h de la cúspide sobre la superficie horizontal. A partir de estos datos, dibujaremos el perfil de la gota y calcularemos el ángulo de contacto θc y el radio R de la base de la gota.

Otro parámetro inicial importante es el radio de curvatura r0 de la cúspide. La línea a trazos es la circunferencia de radio r0 a la que se aproxima la cúspide, cuyo centro se encuentra en el punto de color rojo C

El volumen V de la gota y la altura h se pueden medir fácilmente, sin embargo, la curvatura r0 es muy difícil de medir, así como el radio R de la base debido a las imperfecciones inevitables de la superficie horizontal

Angulo de contacto

El ángulo de contacto, θc que forma la superficie de la gota con el plano horizontal, es independiente del volumen de la gota, como comprobaremos al final de la página.

Calculamos la presión a la altura y, debida:

p= 2γ r 0 +ρgy

Las fuerzas que mantienen el equilibrio de la superficie circular de radio x situada a la altura y son:

El equilibrio de fuerzas en la dirección vertical a la altura y, se escribe

mg+γ( 2πx )sinθ=p( π x 2 )

Si y=h, la región sombreada se extiende al volumen de la gota

ρVg+γ( 2πR )sin θ c =( 2γ r 0 +ρgh )( π R 2 )

ρ es la densidad del líquido y V el volumen de la gota

Despejamos el ángulo de contacto θc

θ c =arcsin( R r 0 +( Rh V πR ) ρg 2γ )

Equilibrio

En este apartado, vamos a calcular las ecuaciones que describen el perfil de la gota a partir de las fuerzas que la mantienen en equilibrio sobre la superficie horizontal

En la figura, mostramos una gota situada sobre una superficie plana y horizontal. Debido a la simetría cilíndrica, situamos el origen en la cúspide de la gota, con el eje Y apuntando hacia abajo y el eje X en la dirección radial

Consideremos un elemento diferencial de gota comprendido entre los planos de alturas y e y+dy. El volumen de este elemento diferencial es πx2·dy y el volumen total V del líquido (supuesto incomprensible) en la gota es

V= 0 h π x 2 dy

Consideremos ahora, su superficie lateral, pintada de color amarillo, corresponde a la de un tronco de cono de bases, b=x+dx, a=x y altura h=dy. El área lateral es

dA=π( 2x+dx ) d y 2 +d x 2 =2πxds

Consideremos las fuerzas sobre dicho elemento diferencial de superficie de la gota:

En el equilibrio se ha de cumplir que

( 2γ r 0 +ρgy )2πx·dx+2πγxsinθ=2πγ(x+dx)sin(θ+dθ) ( 2γ r 0 +ρgy )2πx·dx+2πγxsinθ=2πγ(x+dx)( sinθ+cosθ·dθ ) ( 2γ r 0 +ρgy )2πx·dx=2πγxcosθ·dθ+2πγsinθ·dx

Despejando dx/dθ

dx dθ = γcosθ 2γ r 0 +ρgyγ sinθ x

Integrando esta ecuación diferencial obtenemos la función x=x(θ). La otra ecuación diferencial se obtiene teniendo en cuenta que tanθ=dy/dx

dy dθ = dy dx dx dθ = γsinθ 2γ r 0 +ρgyγ sinθ x

Integrando esta ecuación diferencial obtenemos y=y(θ). De este modo, determinaremos y representaremos el perfil de la gota

Energía mínima

Una gota colocada sobre una superficie horizontal, tiene la forma que hace mínima la suma de la energía potencial gravitatoria y la energía superficial (recuérdese que la tensión superficial, γ, es la energía por unidad de área), manteniendo constante el volumen. En este apartado, vamos a utilizar el cálculo de variaciones, para resolver una situación similar al denominado problema isoperimétrico, otra situación similar es la tratada en la página titulada La forma que adopta una cuerda bajo la acción de su peso y de la tensión superficial

La suma de estas dos energías deberá ser mínima, con la condición de que el volumen V de la gota permanezca constante

V= 0 h π x 2 dy

Formamos la función auxiliar dependiente del parámetro λ

F(y,x, x ˙ )=ρgπ x 2 y+2πγ·x 1+ ( dx dy ) 2 +λπ x 2

Dado que la función F depende de y, x y dx/dy, la ecuación de Euler-Lagrange se escribe (se ha intercambiado x por y)

F x d dy ( F x ˙ )=0

El resultado es

ρgxy+γ 1+ x ˙ 2 +λrλ{ x ˙ 2 1+ x ˙ 2 + x x ¨ ( 1+ x ˙ 2 ) 3/2 }=0 γ( 1 1+ x ˙ 2 x x ¨ ( 1+ x ˙ 2 ) 3/2 )=ρgxyλx

En vez de obtener la solución de la ecuación diferencial de segundo orden, vamos a reemplazarla por dos ecuaciones diferenciales de primer orden que nos da la solución en forma paramétrica, x=x(θ) e y=y(θ) para la forma de la gota

Expresamos la ecuación diferencial en términos del ángulo θ

x ˙ = dx dy = 1 tanθ 1+ x ˙ 2 = 1+ ( dx dy ) 2 = 1 sinθ x ¨ = d 2 x d y 2 = d dy ( 1 tanθ )= d dθ ( 1 tanθ ) dθ dy = 1 sin 2 θ dθ dy

Sustituyendo en la ecuación diferencial de segundo orden

dy dθ = γsinθ ρgyλγ sinθ x

Integrando esta ecuación diferencial obtenemos y=y(θ)

Teniendo en cuenta que

dx dθ dθ dy = 1 tanθ dx dθ = 1 tanθ dy dθ

Obtenemos la segunda ecuación diferencial

dx dθ = γcosθ ρgyλγ sinθ x

Integrando esta ecuación diferencial obtenemos x=x(θ)

Ahora bien, en la cúspide de la gota, cuando x→0, y θ→0, en el denominador tenemos una indeterminación que se resuelve aplicando la regla de L'Hopital

( dx dθ ) 0 = γcos0 λγ cos0 ( dx dθ ) 0

Despejamos la derivada dx/dθ en θ=0,

( dx dθ ) 0 = 2γ λ

Por otra parte, en la cúspide θ→0, se cumple (arco=radio·ángulo). Véase la primera figura

( dx dθ ) 0 = r 0

Donde r0 es el radio de curvatura de la cúspide. De este modo, hemos dado un significado al multiplicador, λ=-2γ/r0.

Finalmente, volvemos a obtener las dos ecuaciones diferenciales que describen el perfil de la gota situada sobre una superficie horizontal

dy dθ = γsinθ ρgy+ 2γ r 0 γ sinθ x dx dθ = γcosθ ρgy+ 2γ r 0 γ sinθ x

Perfil de la gota

Este sistema de dos ecuaciones diferenciales no tiene solución analítica, por lo que integraremos, empleando procedimientos numéricos, el sistema de dos ecuaciones diferenciales para representar el perfil de la gota, partiendo de la cúspide θ=0, x=0 e y=0. Los datos de partida son el volumen V de la gota y la altura h de su cúspide

Sistema de unidades

Conocemos la densidad ρ y la tensión superficial γ del fluido. Para el agua ρ=1000 kg/m3 y γ=72·10-3 N/m. Para que el ordenador no trabaje con números demasiado pequeños, establecemos un sistema de unidades en el que la longitud se mida en cm, la densidad ρ·1000 kg/m3 y la tensión superficial en γ·10-3 N/m

dx· 10 2 dθ = γ· 10 3 cosθ 10 3 ·9.8·ρ·y· 10 2 + 2·γ· 10 3 r 0 · 10 2 γ· 10 3 sinθ x· 10 2 dx dθ = γcosθ 980ρy+ 2·γ r 0 γ sinθ x

Una expresión similar para dy/dθ

dy dθ = γsinθ 980ρy+ 2·γ r 0 γ sinθ x

Indeterminación

La primera dificultad con la que nos tropezamos es que en la cúspide θ=0, x=0, y=0 se produce la indeterminación 0/0 en el término sinθ/x del denominador de las dos ecuaciones. Como hemos visto en el apartado anterior, para θ=0, las derivadas (dy/dθ)0=0 y (dx/dθ)0=r0

No podemos utilizar el procedimiento ode45 de MATLAB, ya que produce un error NaN en la primera iteracción. Utilizaremos el procedimiento de Runge-Kutta descrito en el apartado titulado Sistema de dos ecuaciones diferenciales de primer orden para integrar el sistema de ecuaciones. Para evitar la indeterminación, separaremos la primera iteracción del resto

%sistema de dos ecuaciones diferenciales
f=@(t,x,y) gamma*cos(t)/(2*gamma/r0+rho*980*y-gamma*sin(t)/x);
g=@(t,x,y) gamma*sin(t)/(2*gamma/r0+rho*980*y-gamma*sin(t)/x);
x(1)=0; y(1)=0; %condiciones iniciales en la cúspide

%primera iteracción
i=1;
k1=dt*r0;
l1=0;
k2=dt*f(t(i)+dt/2,x(i)+k1/2,y(i)+l1/2);
l2=dt*g(t(i)+dt/2,x(i)+k1/2,y(i)+l1/2);
k3=dt*f(t(i)+dt/2,x(i)+k2/2,y(i)+l2/2);
l3=dt*g(t(i)+dt/2,x(i)+k2/2,y(i)+l2/2);
k4=dt*f(t(i)+dt,x(i)+k3,y(i)+l3);
l4=dt*g(t(i)+dt,x(i)+k3,y(i)+l3);     
x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;
y(i+1)=y(i)+(l1+2*l2+2*l3+l4)/6;        

%segunda, tercera, n iteracciones
for i=2:n
        k1=dt*f(t(i),x(i),y(i));
        l1=dt*g(t(i),x(i),y(i));
        k2=dt*f(t(i)+dt/2,x(i)+k1/2,y(i)+l1/2);
        l2=dt*g(t(i)+dt/2,x(i)+k1/2,y(i)+l1/2);
        k3=dt*f(t(i)+dt/2,x(i)+k2/2,y(i)+l2/2);
        l3=dt*g(t(i)+dt/2,x(i)+k2/2,y(i)+l2/2);
        k4=dt*f(t(i)+dt,x(i)+k3,y(i)+l3);
        l4=dt*g(t(i)+dt,x(i)+k3,y(i)+l3);
        x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;
        y(i+1)=y(i)+(l1+2*l2+2*l3+l4)/6;
        
        if y(i+1)>=h %completa la integración
            break;
        end
end

El procedimiento numérico termina en la iteracción i+1 cuando la ordenada y alcanza la posición h del plano horizontal, o altura de la gota. El número n de iteracciones deberá ser lo suficientemente grande para asegurarnos de que alcanzamos esta posición

Radio de curvatura de la cúspide

El sistema de dos ecuaciones diferenciales nos pide el radio de curvatura r0 de la cúspide, pero este dato es muy difícil de obtener observando la gota sobre la superficie horizontal. Hacemos que r0=h y calculamos el volumen Vc de la gota empleando el procedimiento numérico denominado fórmula del trapecio

V c = 0 h π x 2 dy

...
vol=0;
for j=1:i
     vol=vol+(y(j+1)-y(j))*(pi*x(j)^2+pi*x(j+1)^2)/2;
end
disp(['volumen: ', num2str(vol)]);

El proceso se repite hasta que se alcance la precisión deseada.

Tenemos una ecuación transcedente V(r0, h)-V=0, cuya solución buscamos aplicando el procedimiento de Newton-Raphson. Dado que no conocemos la función f(x)=V(x)-V=0 de forma explícita ni su derivada, el procedimiento a aplicar es

x n+1 = x n δ V( x n )V V( x n +δ)V( x n )

donde δ es una cantidad pequeña

El proceso de convergencia se termina en un número N de pasos o cuando se alcanza la precisión deseada

| x n+1 x n x n+1 |<ε

...
r0=h; %valor inicial; %cm
%rutina de convergencia hacia el valor r0 que hace que el volumen de la
%gota se aproxime a V
delta=0.0001;
N=20; %número de pasos
for i=1:N
    [R, th_c, V0, x, y,k]=forma_gota_ode(rho, gamma, h, r0);
    [R, th_c, Vc, x, y,k]=forma_gota_ode(rho, gamma, h, r0+delta);
    r1=r0-delta*(V0-V)/(Vc-V0);
    if abs((r1-r0)/r1)<0.001
        break;
    end
    r0=r1;
end
if(i==N)
    disp('No se alcanza la precisión deseada')
end

Radio base y ángulo de contacto

Se termina el proceso de integración de las dos ecuaciones diferenciales, cuando la ordenada y alcanza la posición h del plano horizontal, el procedimiento nos proporciona el radio R de la base y el ángulo de contacto

...
disp([y(i+1),x(i+1),t(i+1)*180/pi]) %altura, radio base, ángulo de contacto

Obtenemos, de forma alternativa, el ángulo de contacto mediante su expresión en el nuevo sistema de unidades establecido

θ c =arcsin( R r 0 +( Rh V πR ) ρg 2γ ) θ c =arcsin( R· 10 2 r 0 · 10 2 +( Rh· 10 4 V· 10 6 πR· 10 2 ) 10 3 ρ9.8 2γ· 10 3 ) θ c =arcsin( R r 0 +( Rh V πR ) 980ρ 2γ )

Hay dos posibles ángulos suplementarios, θc y π-θc

R=x(i+1);  %radio de la base
th_c=asin(R/r0+(R*h-V/(pi*R))*980*rho/(2*gamma));
disp((pi-th_c)*180/pi) %ángulo suplementario

El código

En primer lugar, proporcionamos los datos de:

gamma=72; %tensión superficial
rho=1; %g/cm^3
%unidades para que x e y se midan en cm
V=0.04; %cm^3
h=0.3; %cm
r0=h; %cm

Creamos una función forma_gota_ode, a la que pasamos la densidad ρ, la la tensión superficial γ, la altura de la cúspide sobre el plano horizontal h, y el radio de curvatura de la cúspide r0.

La función devuelve:

Como el perfil de la gota es simétrico respecto del eje vertical Y, se calcula solamente la parte derecha

function [R, th_c, Vc, x, y, k] =forma_gota_ode(rho, gamma, h, r0)
    %sistema de dos ecuaciones diferenciales
    f=@(t,x,y) gamma*cos(t)/(2*gamma/r0+rho*980*y-gamma*sin(t)/x);
    g=@(t,x,y) gamma*sin(t)/(2*gamma/r0+rho*980*y-gamma*sin(t)/x);
    n=31415; %número máximo de iteracciones
    dt=pi/n; %paso
    t=0:dt:pi;
%reserva memoria para n+1 elementos de los vectores x e y
    x=zeros(n+1,1); 
    y=zeros(n+1,1);
    x(1)=0; y(1)=0; %condiciones iniciales
    
    %primera iteracción
    i=1;
    k1=dt*r0;
    l1=0;
    k2=dt*f(t(i)+dt/2,x(i)+k1/2,y(i)+l1/2);
    l2=dt*g(t(i)+dt/2,x(i)+k1/2,y(i)+l1/2);
    k3=dt*f(t(i)+dt/2,x(i)+k2/2,y(i)+l2/2);
    l3=dt*g(t(i)+dt/2,x(i)+k2/2,y(i)+l2/2);
    k4=dt*f(t(i)+dt,x(i)+k3,y(i)+l3);
    l4=dt*g(t(i)+dt,x(i)+k3,y(i)+l3);     
    x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;
    y(i+1)=y(i)+(l1+2*l2+2*l3+l4)/6;  
    %resto de las iteracciones
    for i=2:n
        k1=dt*f(t(i),x(i),y(i));
        l1=dt*g(t(i),x(i),y(i));
        k2=dt*f(t(i)+dt/2,x(i)+k1/2,y(i)+l1/2);
        l2=dt*g(t(i)+dt/2,x(i)+k1/2,y(i)+l1/2);
        k3=dt*f(t(i)+dt/2,x(i)+k2/2,y(i)+l2/2);
        l3=dt*g(t(i)+dt/2,x(i)+k2/2,y(i)+l2/2);
        k4=dt*f(t(i)+dt,x(i)+k3,y(i)+l3);
        l4=dt*g(t(i)+dt,x(i)+k3,y(i)+l3);
        x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;
        y(i+1)=y(i)+(l1+2*l2+2*l3+l4)/6;
        %final del proceso de integración numérica
        if y(i+1)>=h 
            break;
        end
    end
    k=i+1; %iteracción final
    R=x(k); %radio de la base de la gota
    th_c=t(k); %ángulo de contacto
  
  %volumen calculado de la gota 
    Vc=0;
    for j=1:k-1
        Vc=Vc+(y(j+1)-y(j))*(pi*x(j)^2+pi*x(j+1)^2)/2;
    end
 end

Comparamos el volumen calculado Vc con el volumen dado V y modificamos el valor de r0, repetimos el proceso hasta alcanzar la precisión deseada. Una vez alcanzada, representmos el perfil de la gota y proporcionamos los parametros más importantes: ángulo de contacto θc, radio de la base R, radio de curvatura de la cúspide r0

...
%gráficos
hold on
fill([-fliplr(x(1:k));0],[-fliplr(y(1:k));-h],'c')
fill([0;x(1:k)],[-h;-y(1:k)],'c')
plot([-fliplr(x(1:k)), x(1:k)],[-fliplr(y(1:k)),-y(1:k)],'b','lineWidth',1.5)
text(0,-h/2,num2str(r0)) %radio curvatura
text(R/2,-h+0.02,num2str(R)) %radio base
line([-R-0.05, R+0.05],[-h,-h],'lineWidth',1.5,'color','k')
line([R,R+0.05*cos(pi-th_c)], [-h, -h+0.05*sin(pi-th_c)],'color','r')
text(R,-h-0.01,num2str(th_c*180/pi)) %ángulo de contacto
hold off
axis equal
xlabel('x')
ylabel('y')
title('Forma de la gota de agua')

Cambiando el volumen V y la altura h, probamos que el ángulo de contacto θc es independiente del volumen

Volumen medido, VAltura medida, h
0.040.300
0.030.280
0.020.250
0.010.205

Con el mismo el volumen V=0.04, cambiamos la altura h=0.195 cm de la cúspide de la gota. El resultado es

gamma=72; %tensión superficial
rho=1; %g/cm^3
%unidades para que x e y se midan en cm
V=0.04; %cm^3
h=0.195; %cm
r0=h; %valor inicial; %cm
....

Cambiando el volumen V y la altura h, probamos que el ángulo de contacto θc es independiente del volumen

Volumen medido, VAltura medida, h
0.040.195
0.030.180
0.020.160
0.010.130

El programa completo consta de la función forma_gota_ode (en color rosa, más arriba) y del script

gamma=72; %tensión superficial
rho=1; %g/cm^3
%unidades para que x e y se midan en cm
V=0.04; %cm^3
h=0.3; %cm
r0=h; %valor inicial; %cm

%rutina de convergencia hacia el valor r0 que hace que el volumen de la
%gota se aproxime a V
delta=0.0001;
N=20; %número de pasos
for i=1:N
    [R, th_c, V0, x, y,k]=forma_gota_ode(rho, gamma, h, r0);
    [R, th_c, Vc, x, y,k]=forma_gota_ode(rho, gamma, h, r0+delta);
    r1=r0-delta*(V0-V)/(Vc-V0);
    if abs((r1-r0)/r1)<0.001
        break;
    end
    r0=r1;
end
if(i==N)
    disp('No se alcanza la precisión deseada')
end
   
%gráficos
hold on
fill([-fliplr(x(1:k));0],[-fliplr(y(1:k));-h],'c')
fill([0;x(1:k)],[-h;-y(1:k)],'c')
plot([-fliplr(x(1:k)), x(1:k)],[-fliplr(y(1:k)),-y(1:k)],'b','lineWidth',1.5)
text(0,-h/2,num2str(r0))
text(R/2,-h+0.02,num2str(R))
line([-R-0.05, R+0.05],[-h,-h],'lineWidth',1.5,'color','k')
line([R,R+0.05*cos(pi-th_c)], [-h, -h+0.05*sin(pi-th_c)],'color','r')
text(R,-h-0.01,num2str(th_c*180/pi))
hold off
axis equal
xlabel('x')
ylabel('y')
title('Forma de la gota de agua')

Referencias

F. Behroozi, P. S. Behroozi. Reliable determination of contact angle from the height and volume of sessile drops. Am. J. Phys. 87 (1) January 2019, pp. 28-32

Feredoon Behroozi, Hilliard K. Macomber, Jack A. Dostal, Cyrus H. Behroozi. The profile of a dew drop. Am. J. Phys. 64 (9) September 1996, pp. 1120-1125