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

Las corrientes parásitas que nacen del error de curvatura y la técnica well-balanced que las elimina

Una gota en reposo no debería moverse. Las fuerzas están perfectamente equilibradas. Sin embargo, al encender la simulación, pequeños remolinos brotan lentamente sobre la superficie de la gota. Nadie la empuja y aun así el fluido fluye. A este flujo fantasma se le llama corriente parásita (parasitic current: campo de velocidad falso que surge del error numérico). Tallois y colaboradores (2025) atacan este problema con una técnica de tensión superficial well-balanced. Hoy seguimos el rastro de por qué una gota fluye sola y de cómo se la detiene.

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

Una interfaz curvada (frontera donde se encuentran dos fluidos distintos) genera un salto de presión. Esto es la ley de Laplace.

\Delta P = \sigma \kappa

$\Delta P$ es la diferencia de presión entre el interior y el exterior de la interfaz. $\sigma$ es el coeficiente de tensión superficial (fuerza por unidad de longitud). $\kappa$ es la curvatura de la interfaz (el inverso del radio).

Para una gota en dos dimensiones, $\kappa = 1/R$ . Cuanto menor es el radio $R$ , mayor es el salto de presión. En una gota de 1 mm de diámetro, el interior está unos 300 Pa más alto. Esa diferencia de presión es la fuerza que compacta la gota en forma esférica.

Lo clave aquí es que se trata de un equilibrio. La alta presión interior empuja hacia afuera y la tensión superficial tira hacia adentro. Ambas fuerzas son exactamente iguales. Por eso la gota permanece quieta.

Por qué fluye una gota en reposo#

El problema está en que la computadora no logra ajustar ese equilibrio con exactitud.

Las técnicas numéricas convierten la tensión superficial en una fuerza volumétrica. A esto se le llama CSF (Continuum Surface Force: método que difunde la interfaz a lo largo de unas pocas celdas y reparte la fuerza sobre esa banda). La magnitud de la fuerza es proporcional a la curvatura $\kappa$ . Sin embargo, calcular la curvatura con exactitud sobre una malla es difícil.

Basta con calcular la curvatura con un pequeño error para que la fuerza de tensión superficial deje de coincidir con el gradiente de presión. La fuerza restante empuja al fluido.

\mathbf{r} = \nabla P - \sigma \kappa_{\text{num}} \nabla z

$\mathbf{r}$ es la fuerza residual (cero si hay equilibrio). $\kappa_{\text{num}}$ es la curvatura calculada numéricamente. $\nabla z$ es el gradiente de la fracción de volumen (proporción de líquido dentro de la celda).

Si la curvatura es exacta, $\mathbf{r} = 0$ y la gota está en silencio. Pero basta un error de unos pocos por ciento en $\kappa_{\text{num}}$ para que $\mathbf{r} \neq 0$ . Ese residuo crea remolinos a lo largo de la interfaz. Eso es la corriente parásita.

Well-balanced — ajustar con exactitud el equilibrio discreto#

La clave de la solución es una propiedad llamada well-balanced (preservación del equilibrio). Consiste en diseñar las ecuaciones discretizadas de modo que preserven con exactitud la solución en reposo.

Tallois y colaboradores tratan la tensión superficial como un producto no conservativo (non-conservative product). Y colocan este término directamente dentro del solucionador de Riemann (herramienta numérica que resuelve las ondas en la interfaz). Al calcular el gradiente de presión y la fuerza de tensión superficial con la misma regla discreta, ambos se cancelan con exactitud celda por celda.

La forma de calcular la curvatura también importa. El artículo calcula la curvatura con un stencil basado en nodos (node) y no en caras (face). ¿Por qué?

1D (basado en caras): solo mira las caras vecinas. El residuo se alinea con los ejes de la malla. El campo de velocidad rebota en la dirección de la malla y desestabiliza la interfaz.
Multidimensional (basado en nodos): mira todo el entorno del nodo. El residuo se dispersa de forma suave. La corriente parásita es mucho más débil.

La tensión superficial es, en esencia, un fenómeno multidimensional. Por eso conviene un stencil verdaderamente multidimensional.

Probémoslo directamente en la simulación de abajo. Al bajar curvature error ε a cero, el campo de velocidad desaparece. Ese es el estado well-balanced.

curvature error ε: 0.35 surface tension σ: 0.072 N/mdroplet radius R: 0.50 mm

Laplace jump ΔP = σκ = σ/R = 144.0 Pa

kinetic energy Σ½|u|² = 0.00e+0

Al subir ε, crecen remolinos en la interfaz. Al elegir 1D · face, el flujo se alinea con los ejes de la malla y se amplifica con fuerza. Al cambiar a multiD · nodal, con el mismo error el flujo es mucho más débil y suave. Observemos cuánto difiere el valor de energía cinética a la derecha entre ambos casos.

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

La corriente parásita es un flujo falso que hay que eliminar. Pero también existe un movimiento real que crea la tensión superficial.

Si colocamos una gota elíptica dentro de un gas, oscila. Parte de una elipse con energía superficial máxima. La tensión superficial la devuelve a un círculo. En ese momento la energía cinética es máxima. La inercia se pasa de largo y vuelve a formar una elipse. Este vaivén, en condiciones ideales, se prolonga para siempre.

El período de oscilación sigue la fórmula de Rayleigh modificada.

\omega^2 = (n^3 - n)\, \frac{\sigma}{(\rho_l + \rho_g) R^3}

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

Aquí well-balanced vuelve a ser importante. Idealmente, la gota debería oscilar sin amortiguamiento. Pero una técnica con mucha difusión recorta rápidamente la amplitud. La gota colapsa a un círculo en apenas unos pocos períodos. El artículo suprime este amortiguamiento numérico con una corrección low-Mach y mantiene la oscilación durante varios períodos.

Probémoslo directamente abajo.

surface tension σ: 340 N/mliquid density ρ_l: 100 kg/m³numerical damping: 0.40

Rayleigh period T = 2π/ω = 81.2 ms

ω = √[(n³−n)σ / ((ρ_l+ρ_g)R³)] , n = 2

Al dejar numerical damping en cero, la amplitud se mantiene. Esta es la técnica ideal sin difusión. Al subir el valor, la gota se asienta rápidamente en un círculo. Al aumentar σ o reducir ρ_l, el período $T$ se acorta. Tal como dice la fórmula.

El error de equilibrio visto con Python#

Comprobemos directamente cómo el error de curvatura da origen a la corriente parásita. Tomamos una sección radial de la gota. Difuminamos la fracción de volumen de forma suave con tanh. Y calculamos 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)))

El resultado de la ejecución es este.

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

Con la curvatura exacta, el residuo está al nivel del error de máquina. La energía cinética es prácticamente cero. Con un error de apenas 3 % en la curvatura, el residuo explota. La energía cinética salta 24 órdenes de magnitud. Esa es la verdadera identidad de la corriente parásita que vimos en pantalla.

Para recordar#

La ley de Laplace $\Delta P = \sigma \kappa$ es el equilibrio de una gota en reposo. La presión interior y la tensión superficial se cancelan con exactitud.
La corriente parásita es el campo de velocidad falso que surge cuando el error de curvatura rompe esa cancelación. Basta un error del 3 % en la curvatura para que la energía cinética explote.
La técnica well-balanced calcula la presión y la tensión superficial con la misma regla discreta y preserva el equilibrio con exactitud. El stencil multidimensional basado en nodos evita la inestabilidad de alineación con la malla.