Lo que mata un vórtice no es la viscosidad — el 1/r de Newton y el núcleo de Lamb–Oseen

Una condición inicial habitual para verificar un solver viscoso es el vórtice libre. Se impone el campo $v_\theta = \Gamma/2\pi r$ y se avanza en el tiempo. Como hay viscosidad, el vórtice debería debilitarse poco a poco. Sin embargo, la velocidad en el radio 0.5 sigue igual hasta la cuarta cifra decimal después de cincuenta veces el tiempo de difusión inicial. El código no está mal. Este artículo explica por qué, y qué es lo que realmente mata un vórtice. La mitad de la respuesta ya estaba impresa en 1687.

A qué apuntaba el Libro II de los Principia#

Los Principia tienen tres libros. El primero trata las leyes de la fuerza y el tercero la gravitación universal. El segundo, el que casi nadie cita en física, es mecánica de fluidos. Su blanco era explícito: la cosmología de vórtices de Descartes.

Descartes explicaba las órbitas planetarias como cuerpos arrastrados por un vórtice fluido. El espacio estaría lleno de materia sutil, y esa materia formaría un enorme flujo rotatorio. Newton atacó esa imagen con argumentos de mecánica de fluidos. El movimiento de un fluido siempre lleva resistencia asociada, así que sin una fuerza externa el vórtice terminará por extinguirse.

Para sostener el argumento había que cuantificar la resistencia. Por eso en el Libro II se postula que la resistencia de un fluido es lineal en la tasa de cizalla.

Aquí $\tau$ es el esfuerzo cortante, $du/dy$ el gradiente de velocidad y la constante de proporcionalidad $\mu$ es la viscosidad. Es el primer lugar donde la viscosidad de un fluido queda definida matemáticamente. Más tarde se descubrió que son muchos más los fluidos que incumplen esa relación lineal que los que la cumplen, y de ahí que a estos últimos se les llame newtonianos.

El 1/r que sale de un solo cilindro#

La derivación de Newton sigue leyéndose como moderna. Un cilindro infinitamente largo gira a velocidad angular constante dentro de un fluido viscoso. En estado estacionario, el par transmitido hacia afuera a través de cualquier superficie cilíndrica de radio $r$ tiene que ser el mismo en todas partes. De lo contrario se acumularía momento angular en la capa intermedia.

En coordenadas cilíndricas, el esfuerzo cortante de un flujo puramente circular es $\tau_{r\theta} = \mu\, r\, \frac{d}{dr}(v_\theta/r)$. No es simplemente $\mu\, dv_\theta/dr$. La rotación rígida ($v_\theta \propto r$) no debe llevar cizalla, y esta forma cumple esa condición de manera automática. El par por unidad de longitud es el esfuerzo por el brazo $r$ por el perímetro $2\pi r$.

Al resolver esa ecuación diferencial ordinaria aparecen dos términos.

El primero es rotación rígida y el segundo es un vórtice libre. Si la frontera exterior está infinitamente lejos y allí el fluido está en reposo, entonces $A=0$. Lo que queda es $v_\theta \propto 1/r$. Ese es el resultado del Libro II, y es la misma expresión que hoy se enseña como solución estacionaria del flujo de Taylor–Couette.

El exponente que exige Kepler es 1/2#

Aquí se cierra la refutación. Si el vórtice de Descartes transporta los planetas, su perfil de velocidad tiene que reproducir la tercera ley de Kepler. El periodo orbital escala como $T \propto r^{3/2}$, de modo que la velocidad escala como $v = 2\pi r / T \propto r^{-1/2}$.

El vórtice fluido de Newton da $r^{-1}$. Kepler exige $r^{-1/2}$. Los exponentes no coinciden. Un vórtice estacionario en un fluido viscoso no puede producir órbitas planetarias.

Conviene mover el exponente a mano en el dial de abajo.

En n = −1 el radio marcado se mantiene recto: eso es rotación rígida. Al empujar hasta n = 1 el radio se enrolla con fuerza y, aun así, la barra de "viscous force" abajo a la derecha cae a verde. En el n = 0.5 de Kepler las dos barras están en rojo. Ahí está el punto de observación: sostener ese perfil exige que algo siga empujando.

Hay esfuerzo, pero no hay fuerza#

Con esto se resuelve la duda inicial. Basta sustituir $v_\theta = C\,r^{-n}$ en el término viscoso de la componente rotacional.

Ese último término $-v_\theta/r^2$ aparece por la curvatura. El laplaciano cartesiano no tiene nada equivalente. Al agrupar el paréntesis queda $n^2-1$.

Ese coeficiente es exactamente cero en $n = \pm 1$. El caso $n = -1$ es rotación rígida. El caso $n = +1$ es el vórtice libre. Es decir, sobre el vórtice libre no actúa ninguna fuerza viscosa.

No porque el esfuerzo cortante se anule. El esfuerzo vale $\tau_{r\theta} = -\mu\Gamma/\pi r^{2}$, enorme cerca del núcleo. Lo que ocurre es que el par sobre la cara interior de un elemento fluido cancela exactamente el par sobre su cara exterior. La fuerza no es el esfuerzo, es la divergencia del esfuerzo. El vórtice libre es el perfil especial cuya divergencia de esfuerzos vale cero.

La misma historia se cuenta con vorticidad. Con $\omega = \frac{1}{r}\frac{d(r v_\theta)}{dr}$ y $v_\theta = \Gamma/2\pi r$, el producto $r v_\theta$ es constante, así que $\omega = 0$ en todo punto con $r>0$. La fuerza viscosa en flujo incompresible puede escribirse como $-\nu\,\nabla\times\boldsymbol{\omega}$. Sin vorticidad no hay fuerza viscosa.

El reloj lo fija la expansión del núcleo#

Se dijo que no hay vorticidad, pero con precisión eso vale salvo en el origen. La circulación $\Gamma$ tiene que estar en alguna parte, y en el vórtice libre ideal está concentrada en una delta sobre el eje. Ahí es donde la viscosidad trabaja de verdad.

Si se concentra la circulación $\Gamma$ en el origen y se resuelve el problema de difusión, sale el vórtice de Lamb–Oseen.

La exponencial construye el núcleo. Para $r \gg \sqrt{4\nu t}$ el paréntesis vale 1 y se recupera el vórtice libre. El radio del núcleo crece como $r_c = 2.2418\sqrt{\nu t}$ y la velocidad pico cae como $t^{-1/2}$.

La cantidad clave es la circulación. Cuando $r \to \infty$ vale siempre exactamente $\Gamma$. La viscosidad reparte la vorticidad; no la elimina. Conviene alternar el interruptor de pared en el experimento de abajo.

Con la pared apagada el núcleo se hincha y el pico baja, pero la barra de circulación de abajo a la derecha no se mueve de 1.000. Al encender la pared esa misma barra empieza a caer. Lo que se manipula es $\nu$ y la pared; lo que se observa es el radio del núcleo y la diferencia entre las dos historias de circulación.

Tasas de decaimiento por radio, contadas en Python#

Toca poner números a lo dicho. Se avanza la ecuación axisimétrica de giro $\partial_t v_\theta = \nu(\partial_{rr} v_\theta + \partial_r v_\theta / r - v_\theta / r^2)$ con Euler explícito sobre 400 celdas radiales centradas.

import numpy as np
 
NU, GAMMA, R, N = 1.0e-3, 1.0, 1.0, 400
dr = R / N
r = (np.arange(N) + 0.5) * dr
 
def lamb_oseen(rr, t):
    """Velocidad tangencial del vortice de Lamb-Oseen; con t pequeno tiende al 1/r puro."""
    return GAMMA / (2 * np.pi * rr) * (1 - np.exp(-rr * rr / (4 * NU * t)))
 
def swirl_terms(v):
    """Las tres piezas del laplaciano de giro, separadas para ver la cancelacion."""
    ghost = np.concatenate(([-v[0]], v, [v[-1] * r[-1] / (r[-1] + dr)]))
    return ((ghost[2:] - 2 * ghost[1:-1] + ghost[:-2]) / dr**2,
            (ghost[2:] - ghost[:-2]) / (2 * dr) / r,
            -v / r**2)
 
def swirl_operator(v, outer):
    """nu * (v_rr + v_r/r - v/r^2): todo lo que queda de la viscosidad en giro puro."""
    g_out = -v[-1] if outer == 'wall' else v[-1] * r[-1] / (r[-1] + dr)
    w = np.concatenate(([-v[0]], v, [g_out]))
    v_rr = (w[2:] - 2 * w[1:-1] + w[:-2]) / dr**2
    v_r = (w[2:] - w[:-2]) / (2 * dr)
    return NU * (v_rr + v_r / r - v / r**2)
 
def march_swirl(v, t_end, outer):
    dt = 0.2 * dr * dr / NU
    for _ in range(int(round(t_end / dt))):
        v = v + dt * swirl_operator(v, outer)
    return v
 
def circulation_at(v, radius):
    return 2 * np.pi * radius * np.interp(radius, r, v)
 
# 1. el vortice 1/r puro: tres terminos grandes que se cancelan
a, b, c = (np.interp(0.10, r, x) for x in swirl_terms(GAMMA / (2 * np.pi * r)))
print(f"free vortex at r=0.10:  v_rr={a:+8.2f}  v_r/r={b:+8.2f}  -v/r^2={c:+8.2f}  sum={a+b+c:+.2e}")
print(f"                        shear stress tau_rtheta = {-NU * GAMMA / (np.pi * 0.10**2):+.4f}")
 
# 2. avanzar un vortice real y contrastarlo con la solucion analitica
t0, t1 = 0.02, 1.0
v0 = lamb_oseen(r, t0)
v1 = march_swirl(v0, t1 - t0, 'free')
print(f"\nmarched {t0} -> {t1} s   L-inf vs Lamb-Oseen = {np.max(np.abs(v1 - lamb_oseen(r, t1))):.1e}")
print("    r    v(0.02)   v(1.00)    change")
for x in (0.02, 0.05, 0.15, 0.50, 0.95):
    p, q = np.interp(x, r, v0), np.interp(x, r, v1)
    print(f"{x:5.2f} {p:9.4f} {q:9.4f} {100 * (q - p) / p:+8.1f}%")
 
# 3. el nucleo crece como sqrt(nu t), el pico cae como 1/sqrt(t), la circulacion aguanta
print("\n    t   r_core   v_peak  v_peak*sqrt(t)  Gamma(0.9)")
v, tc = v0.copy(), t0
for t in (0.05, 0.20, 0.50, 1.00):
    v = march_swirl(v, t - tc, 'free'); tc = t
    i = int(np.argmax(v))
    print(f"{t:5.2f} {r[i]:8.4f} {v[i]:8.4f} {v[i] * np.sqrt(t):13.4f} {circulation_at(v, 0.9):11.4f}")
 
# 4. con una pared en r = R el mismo vortice muere
print("\n    t   Gamma(0.9)  unbounded    walled")
vw, tc = v0.copy(), t0
for t in (1.0, 20.0, 70.0, 200.0):
    vw = march_swirl(vw, t - tc, 'wall'); tc = t
    print(f"{t:8.1f}     {circulation_at(lamb_oseen(r, t), 0.9):9.4f} {circulation_at(vw, 0.9):9.4f}")
print(f"\nslowest walled mode:  R^2/(nu*j11^2) = {R**2 / (NU * 3.8317**2):.1f} s")

free vortex at r=0.10:  v_rr= +318.81  v_r/r= -159.40  -v/r^2= -159.30  sum=+9.98e-02
                        shear stress tau_rtheta = -0.0318
 
marched 0.02 -> 1.0 s   L-inf vs Lamb-Oseen = 7.7e-04
    r    v(0.02)   v(1.00)    change
 0.02    7.9233    0.7562    -90.5%
 0.05    3.1851    1.4782    -53.6%
 0.15    1.0611    1.0573     -0.4%
 0.50    0.3183    0.3183     +0.0%
 0.95    0.1675    0.1675     +0.0%
 
    t   r_core   v_peak  v_peak*sqrt(t)  Gamma(0.9)
 0.05   0.0163   7.1699        1.6032      1.0000
 0.20   0.0312   3.5880        1.6046      1.0000
 0.50   0.0513   2.2697        1.6049      1.0000
 1.00   0.0713   1.6057        1.6057      1.0000
 
    t   Gamma(0.9)  unbounded    walled
     1.0        1.0000    0.9773
    20.0        1.0000    0.4182
    70.0        0.9446    0.2186
   200.0        0.6367    0.0342
 
slowest walled mode:  R^2/(nu*j11^2) = 68.1 s

El primer bloque muestra la cancelación. Tres términos del orden de 300 suman 0.1. En términos relativos eso es $3\times10^{-4}$, y corresponde al error de discretización de la diferencia de segundo orden. Analíticamente la suma es exactamente cero. El esfuerzo cortante en el mismo punto no lo es.

El segundo bloque responde la duda inicial. En $r = 0.5$ y $r = 0.95$ el cambio es del 0.0%, tras cincuenta veces el tiempo inicial. Solo $r = 0.02$ perdió un 90%. La viscosidad trabajó en el núcleo y en ningún otro sitio.

Las dos últimas columnas del tercer bloque confirman el escalado. v_peak*sqrt(t) se mantiene entre 1.603 y 1.606, un margen del 0.2%. El radio del núcleo coincide: $2.2418\sqrt{\nu t} = 0.0709$ frente a 0.0713 medido. Durante todo ese rato la circulación marca 1.0000.

Con una pared alrededor, ese mismo vórtice muere#

El cuarto bloque es la conclusión. En el caso no acotado, la circulación en $r=0.9$ todavía marca 0.6367 en $t = 200$. Eso no es decaimiento: es el núcleo que ya se hinchó hasta ese radio. Medida en un radio mayor, sigue siendo 1.

El caso con pared marca 0.0342. Se fue más del 95%. Lo que hace la pared es sacar momento angular del sistema. El tiempo de decaimiento lo fija el modo más lento. El modo que cumple $v_\theta = 0$ en la pared y regularidad en el eje es $J_1(j_{1,1} r/R)$, con constante de tiempo $R^2/(\nu\, j_{1,1}^2) = 68.1$ segundos. La corrida marca 0.2186 en $t = 70$, lo cual encaja bien.

Así que la frase de Newton admite una enmienda. Lo que mata un vórtice no es la viscosidad por sí sola. Un vórtice muere cuando hay viscosidad y una frontera por donde pueda escapar el momento angular. En un dominio no acotado la viscosidad solo reparte vorticidad. Para refutar a Descartes eso bastaba igual: con universo finito o infinito, un planeta montado en un vórtice no reproduce el exponente de Kepler.

Antes de escribir el próximo caso de prueba con vórtices#

De aquí salen tres consecuencias prácticas.

No conviene usar un vórtice libre para verificar precisión. Fuera del núcleo el término viscoso es idénticamente cero, de modo que cualquier error allí es error del esquema convectivo y no de la discretización viscosa. Para verificar el término viscoso hay que resolver en malla un núcleo de Lamb–Oseen y medir el crecimiento de $r_c(t)$.

Sí conviene usarlo para medir amortiguamiento numérico. Al revés, esa misma propiedad es un diagnóstico excelente. Se siembra un vórtice libre y toda caída de velocidad fuera del núcleo es disipación numérica del código, no física. Los esquemas de la familia upwind quedan al descubierto de inmediato.

El tamaño del dominio cambia la respuesta. En problemas donde el vórtice debe sobrevivir mucho tiempo — vórtices de punta de ala, estelas de rotor — una frontera exterior demasiado cercana se comporta como pared. Antes de fijar el dominio conviene comparar $R^2/(\nu\, j_{1,1}^2)$ con el tiempo físico de interés. Si la corrida usa modelo de turbulencia, en el lugar de $\nu$ va $\nu_t$ y esa constante de tiempo se acorta en órdenes de magnitud.