Por qué una gota en reposo fluye sola — corrientes parásitas y tensión superficial well-balanced

La gota estaba en reposo. Las fuerzas coincidían. Al arrancar el cálculo crecieron vórtices pequeños en la superficie. Nadie empujó. El fluido se movió solo. Ese campo de velocidad falso es una corriente parásita (parasitic current). Tallois y colaboradores (2025) la apagan con una tensión superficial well-balanced.

Al bajar curvature error ε a cero, el campo de velocidad desaparece. Al subir el mismo control, los vórtices regresan sobre la interfaz.

El cuadro con ε en cero es el estado equilibrado. Al pasar a 1D · face, el mismo error pega el flujo a los ejes de la malla. multiD · nodal reparte ese mismo error en un campo mucho más débil. La energía cinética de la derecha es el número de esa diferencia.

La ley de Laplace — el salto de presión que crea la curvatura#

Una interfaz curvada, la frontera donde se encuentran dos fluidos, produce un salto de presión. Es la ley de Laplace.

$\Delta P$ es la diferencia de presión a través de la interfaz. $\sigma$ es el coeficiente de tensión superficial, fuerza por unidad de longitud. $\kappa$ es la curvatura, el inverso del radio.

En una gota 2D, $\kappa=1/R$. Un $R$ menor da un salto mayor. Un diámetro de 1 mm, es decir $R=0.5\,\mathrm{mm}$, leído como esfera, tiene $\kappa=2/R$ y $\Delta P=\sigma\cdot 2/R$. Con el $\sigma=0.072\,\mathrm{N/m}$ del agua sale $0.072\times 4000=288\,\mathrm{Pa}$, los cerca de 300 Pa de la frase habitual. El código es un cilindro 2D, así que $\kappa=1/R=2000\,\mathrm{m^{-1}}$ y el salto es $0.072/0.0005=144\,\mathrm{Pa}$. La mitad de la esfera.

La presión interior empuja hacia afuera. La tensión superficial tira hacia adentro. Cuando las dos coinciden, la gota se queda quieta. Ese estado es el equilibrio.

Por qué fluye una gota en reposo#

La computadora no clava ese equilibrio sobre las celdas.

Los métodos numéricos convierten la tensión superficial en una fuerza de volumen. Eso es CSF (Continuum Surface Force). La interfaz se difumina a lo largo de unas pocas celdas y la fuerza se reparte en esa banda. La fuerza escala con $\kappa$. La curvatura sobre una malla se tuerce con facilidad. Empujar la misma fuerza hasta el ángulo de contacto en una pared es otro cierre. Lo que se rompe aquí es anterior. La curvatura y la presión dejan de cancelarse.

Un error chico de curvatura y la fuerza de tensión superficial no coincide con el gradiente de presión. Lo que sobra empuja el fluido.

$\mathbf{r}$ es la fuerza residual. Vale cero cuando el equilibrio se sostiene. $\kappa_{\text{num}}$ es la curvatura numérica. $\nabla z$ es el gradiente de la fracción de volumen, la parte líquida de la celda.

Con $\kappa_{\text{num}}$ exacta, $\mathbf{r}=0$. Unos pocos por ciento de error y $\mathbf{r}$ ya no es cero. Ese residuo arma vórtices a lo largo de la interfaz. La ε de la figura de arriba son esos pocos por ciento.

Well-balanced — ajustar con exactitud el equilibrio discreto#

Si la discretización conserva la solución en reposo como solución, es well-balanced (preservación del equilibrio).

Tallois y colaboradores tratan la tensión superficial como un producto no conservativo (non-conservative product) y meten el término dentro del solucionador de Riemann. El gradiente de presión y la fuerza de tensión superficial, calculados con la misma regla discreta, se cancelan celda por celda. Cuando un problema que debería quedarse quieto salta en presión, el criterio de esa ruptura está en el sitio donde un término no conservativo rompe una solución en reposo. En la gota, ese término es la curvatura.

El stencil de la curvatura también parte el resultado. Un stencil 1D basado en caras solo mira las caras vecinas. El residuo se alinea con los ejes de la malla, la velocidad pega un salto en esa dirección y la interfaz se sacude. Un stencil multidimensional basado en nodos mira alrededor del nodo. El residuo se dispersa y la corriente parásita se debilita. La tensión superficial es un fenómeno multidimensional, así que el stencil también tiene que serlo. La brecha de energía cinética entre 1D · face y multiD · nodal en la figura de arriba es esa decisión.

La gota que vuelve a moverse — oscilación de Rayleigh#

Una corriente parásita es un flujo falso que hay que quitar. La tensión superficial también crea un movimiento real.

Una elipse en un gas oscila. Parte de donde la energía superficial es máxima. La tensión superficial la devuelve a un círculo. Ahí la energía cinética llega al máximo. La inercia se pasa y la elipse regresa. Sin amortiguamiento, el intercambio sigue.

El período sigue una fórmula de Rayleigh modificada.

$\omega=2\pi/T$ es la frecuencia angular. $n$ es el modo. Una elipse es $n=2$. $\rho_l,\rho_g$ son las densidades del líquido y del gas. $R$ es el radio medio en reposo.

Con $n=2$, $n^3-n=8-2=6$. Con el $\sigma=0.072$ del código y $R=5\times 10^{-4}$, $R^3=1.25\times 10^{-10}$ y $6\sigma/R^3=3.456\times 10^{9}$. $\omega^2$ es ese número dividido por $\rho_l+\rho_g$. Este texto no da las dos densidades como cifras. Un denominador mayor alarga el período. Subir $\sigma$ agranda el numerador y acorta el período.

Una gota ideal oscila sin amortiguamiento. Un método muy difusivo recorta la amplitud hasta un círculo en unos pocos períodos. El artículo contiene ese amortiguamiento numérico con una corrección de bajo Mach.

El error de equilibrio visto con Python#

En un corte radial, se difumina la fracción de volumen con una tanh y se calcula la fuerza residual.

import numpy as np
 
sigma = 0.072            # N/m, agua-aire
R = 0.5e-3               # radio de la gota (m)
kappa_exact = 1.0 / R    # curvatura de cilindro 2D
 
def volume_fraction(r, radius, width):
    # fraccion de liquido con la interfaz difuminada (0=gas, 1=liquido)
    return 0.5 * (1.0 - np.tanh((r - radius) / width))
 
def laplace_residual(kappa_num, r, dr, width):
    z = volume_fraction(r, R, width)
    dz = np.gradient(z, dr)                         # d z / d r
    f_st = sigma * kappa_num * dz                   # fuerza de tension superficial CSF
    p = np.cumsum(sigma * kappa_exact * dz) * dr    # presion de equilibrio de Laplace
    dp = np.gradient(p, dr)
    return dp - f_st                                # fuerza residual (0 si es well-balanced)
 
def parasitic_energy(residual, dt=1e-6, steps=200):
    # la fuerza restante acelera el fluido: u <- u + dt * residual
    u = np.zeros_like(residual)
    for _ in range(steps):
        u += dt * residual
    return 0.5 * np.sum(u * u)
 
r = np.linspace(0.2e-3, 0.8e-3, 400)
dr = r[1] - r[0]
width = 3 * dr
 
res_ok  = laplace_residual(kappa_exact,        r, dr, width)   # curvatura exacta
res_bad = laplace_residual(kappa_exact * 1.03, r, dr, width)   # error del 3%
 
print("curvatura exacta : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_ok).max(),  parasitic_energy(res_ok)))
print("error del 3%%  : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_bad).max(), parasitic_energy(res_bad)))

La corrida imprime esto.

curvatura exacta : max|r| = 4.1e-10,  KE = 2.3e-19
error del 3% : max|r| = 4.8e+05,  KE = 1.7e+05

La curvatura exacta tiene un residuo máximo de $4.1\times 10^{-10}$. La energía cinética $2.3\times 10^{-19}$ es prácticamente cero. El caso del 3%, $\kappa$ multiplicada por 1.03, tiene un residuo máximo de $4.8\times 10^{5}$ y una energía cinética de $1.7\times 10^{5}$. La razón es $1.7\times 10^{5}/2.3\times 10^{-19}\approx 7.4\times 10^{23}$, cerca de 23.9 órdenes. Esa es la línea de las 24 órdenes.

Ese $4.8\times 10^{5}$ es el 3% del salto de presión, convertido en gradiente. $\sigma\cdot 0.03\cdot\kappa=0.072\times 60=4.32\,\mathrm{Pa}$ es el error de presión que queda en la interfaz. Sumar $dt=10^{-6}$ doscientas veces da $2\times 10^{-4}$. Por el residuo máximo, la escala de velocidad es $2\times 10^{-4}\times 4.8\times 10^{5}=96$. La fuerza que deja un 3% apila esa escala de velocidad sobre la interfaz.

Con numerical damping en cero, la elipse conserva la amplitud. Al subir el valor, el colapso a un círculo se acelera. Al aumentar σ o reducir ρ_l, el período $T$ se acorta. El numerador de arriba es $6\sigma/R^3$, así que duplicar $\sigma$ duplica ese numerador y multiplica $\omega$ por $\sqrt{2}$.

El radio de 0.5 mm cae en el medio de las muestras. El código coloca 400 puntos de 0.2 mm a 0.8 mm. El tramo mide 0.6 mm en 399 huecos, así que un hueco mide cerca de 1.50 μm. La interfaz ocupa tres huecos, cerca de 4.51 μm. El radio es unas 111 veces ese ancho, y la banda difusa queda delgada junto a la gota.

La curvatura de esfera 2/R, con la misma tensión superficial, da un salto de 288 Pa, y el 3% de ese salto son 8.64 Pa. El cilindro del código tiene la mitad de curvatura, así que la presión que sobra también es la mitad, 4.32 Pa. Meter la fórmula de la esfera en este residuo bidimensional lee la fuerza del error al doble.

Duplicar el radio multiplica el cubo por 8, así que la frecuencia angular baja por la raíz de 8, cerca de 2.83, y el período crece en ese factor. Duplicar la tensión superficial solo multiplica la frecuencia por la raíz de 2, cerca de 1.41. El radio empuja el período con más fuerza. La suma de densidades sigue sin ser un número. El período en segundos no se cierra aquí.

Dejar el 3% y solo agrandar el paso de tiempo sigue siendo un fallo#

Si la gota en reposo no está quieta, la curvatura se ajusta antes de medir un período. Mientras el residuo $4.8\times 10^{5}$ sigue vivo, un $\Delta t$ mayor solo apila más rápido la escala de velocidad 96. Well-balanced es lo que baja ese residuo a $4.1\times 10^{-10}$ y da derecho a hablar del paso de tiempo. Dejar el stencil en las caras devuelve el flujo axial y retira ese derecho.

El límite capilar del paso de tiempo sigue ahí después de que muere el residuo. Hasta dónde se puede empujar $\Delta t$ no lo responde este residuo. Qué ancho abre de verdad el paso de tiempo cuando la tensión superficial se trata de forma implícita carga el fallo siguiente.