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.

Al empujar el dial del exponente hasta n = 1, la barra de "viscous force" abajo a la derecha cae a verde. Ese es el lado del vórtice libre.

En el n = 0.5 de Kepler las dos barras se quedan en rojo. Sostener ese perfil exige un empuje continuo.

La solución del cilindro de 1687 ya dice la mitad. El reloj lo llevan el núcleo y la pared.

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 término de delante es rotación rígida y el de detrás 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. Un vórtice de Rankine que junta rotación forzada y vórtice libre en un mismo círculo está armado con esos dos términos.

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. En el dial del principio, n = −1 deja el radio recto, y n = 1 lo enrolla mientras solo la barra de fuerza viscosa pasa a verde.

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.

En el exponente de Kepler $n=1/2$, $n^2-1=1/4-1=-3/4$. No es cero. El coeficiente desaparece solo en $n=1$. Los tres términos que salen más adelante en r=0.10, 318.81, −159.40 y −159.30, suman apenas $+9.98\times 10^{-2}$ por ese cero. Frente a 318.81 eso es cerca de $3\times 10^{-4}$, el resto de una diferencia de segundo orden. El esfuerzo cortante ahí, −0.0318, no es cero. El esfuerzo está. Solo falta su divergencia.

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 referencia es la circulación. Cuando $r \to \infty$ vale siempre exactamente $\Gamma$. La viscosidad reparte la vorticidad; no la elimina. Sin pared, esa $\Gamma$ se queda dentro del dominio.

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

Al cambiar ν y el interruptor de pared, el radio del núcleo se mueve junto con la barra de circulación de abajo a la derecha.

Con la pared apagada el pico baja y la barra de circulación sigue pegada a 1.000. Con la pared encendida esa barra cae, y la circulación con pared 0.0342 en t=200 es el final.

El bloque inicial 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 bloque de tiempo 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 columnas del pico 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.

Llevar t de 0.02 a 1.00 es un factor de 50. El 0.3183 en r=0.50 coincide con $\Gamma/(2\pi r)=1/(2\pi\cdot 0.5)$, y el cambio es +0.0%. Solo se mueve r=0.02, un −90.5%. El radio del núcleo crece de 0.0163 en t=0.05 a 0.0713 en t=1.00. Con $\nu=10^{-3}$ y t=1, $2.2418\sqrt{\nu t}$ da 0.0709. La diferencia es 0.0004, menor que el espaciamiento dr=1/400=0.0025. v_peak*sqrt(t) queda entre 1.6032 y 1.6057, dentro del 0.2%.

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

El bloque de la pared 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. Queda solo el 3.4% de la circulación inicial de 1. Junto al caso no acotado de 0.6367 en el mismo instante, la pared se llevó mucho más. 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. El $j_{1,1}=3.8317$ del código, con $\nu=10^{-3}$ y $R=1$, da $3.8317^2\approx 14.68$, y $1/(10^{-3}\times 14.68)\approx 68.1$. 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.

Una caída de velocidad fuera del núcleo no es viscosidad#

Usar el vórtice libre fuera del núcleo para verificar precisión no examina la discretización viscosa. El término viscoso ahí es analíticamente cero. El error de esa zona es del esquema convectivo. Para ver el término viscoso se cubre un núcleo de Lamb–Oseen y se mide el crecimiento de $r_c(t)=2.2418\sqrt{\nu t}$. Con $\nu=10^{-3}$ y t=1 eso vale 0.0709, que coincide con el radio de núcleo impreso, 0.0713.

Al revés, una caída de velocidad fuera del núcleo es disipación numérica, no la física. Los esquemas upwind quedan al descubierto ahí. En un cálculo donde el vórtice tiene que sobrevivir, una frontera exterior demasiado cercana actúa como pared. Antes de fijar el dominio se compara $R^2/(\nu\, j_{1,1}^2)=68.1$ segundos con el tiempo físico. Con un modelo de turbulencia, $\nu_t$ ocupa el lugar de $\nu$ y esa constante de tiempo se acorta mucho. Es el mismo juicio que el número de Reynolds, que intercambia los papeles de la viscosidad y la inercia. Con viscosidad sola y sin salida para el momento angular, el vórtice no muere.