El caos que empieza en el fondo de una olla — Convección de Rayleigh–Bénard y las ecuaciones de Lorenz

En 1900 Henri Bénard calentó desde el fondo una capa fina de aceite y la vio partirse sola en un panal hexagonal. Nadie la removió. Dieciséis años después Lord Rayleigh escribió el umbral como un solo número. Cuarenta y siete años más tarde, Lorenz sacó el caos del mismo sistema.

El umbral es el número de Rayleigh 1708. Por debajo, conducción. Por encima, rollos. Al cortar el rollo a tres ecuaciones queda el parámetro de calentamiento ρ. Pasado 24.74, los puntos fijos ya no sujetan la trayectoria.

Al bajar el deslizador de Ra por debajo de 1708 las partículas casi se detienen. Al empujarlo por encima de 1708 se encienden rollos que giran en sentidos opuestos.

Por debajo de 1708 la temperatura solo se apila en capas horizontales. Pasado el umbral, empujar más el deslizador hace girar más rápido los mismos rollos.

Por qué un líquido en reposo empieza a moverse#

Caliente abajo, frío arriba. El fluido calentado se expande y se aligera. Lo ligero abajo, lo pesado arriba. La flotabilidad quiere voltear esa pila.

No salta de inmediato. La viscosidad frena el movimiento. La difusión térmica deja escapar el calor de una burbuja que sube y mata la flotabilidad que la elevaba. Las dos son frenos.

Si gana la flotabilidad aparece una circulación, un rollo de convección. Si ganan los frenos, el líquido se queda y el calor viaja solo por conducción. La flotabilidad también se gasta en el camino, antes de que exista un rollo. El tiro que pierde 7 Pa cuando el humo baja 1 m es ese presupuesto.

El número de Rayleigh — fuerza que voltea sobre fuerza que retiene#

Comprimida en un solo número adimensional, la pulseada es el número de Rayleigh.

$g$ es la aceleración de la gravedad, $\beta$ el coeficiente de expansión térmica (cuánto crece el volumen al subir la temperatura), $\Delta T$ la diferencia de temperatura entre arriba y abajo, $L$ el espesor de la capa, $\nu$ la viscosidad cinemática y $\alpha$ la difusividad térmica. El numerador voltea. El denominador retiene.

El término pesado es $L^3$. Al engrosar la capa, Ra crece al cubo. Una película fina de aceite rara vez convecta. El océano profundo y el manto se mueven con una diferencia de temperatura minúscula. Una olla apenas más honda acerca el umbral de golpe.

Número de Rayleigh crítico 1708 — el umbral donde se enciende la convección#

El análisis de estabilidad lineal de Rayleigh cambia el número con la frontera. Arriba y abajo libres: $\mathrm{Ra}_c = 27\pi^4/4 \approx 657.5$. Las paredes rígidas habituales en los experimentos: $\mathrm{Ra}_c \approx 1708$.

Cuando $\mathrm{Ra} < \mathrm{Ra}_c$, la viscosidad y la difusión devoran una oscilación pequeña. El fluido vuelve al reposo. En cuanto $\mathrm{Ra} > \mathrm{Ra}_c$, una oscilación de cierta longitud de onda crece sola. Esa es la semilla de un rollo.

Dividir 1708 entre 657.5 da cerca de 2.6. Mismo fluido, mismo espesor: las paredes rígidas piden una diferencia de temperatura unas 2.6 veces mayor. Pagado en espesor, la raíz cúbica es $2.6^{1/3} \approx 1.37$. Ese es el coste de cambiar una tapa de deslizamiento libre por una rígida.

La forma de un rollo y el número de Nusselt#

Justo pasado el umbral, el rollo elige una longitud de onda. El número de onda crítico es $k_c \approx 3.117/L$. Una longitud de onda es $2\pi/k$, y al dividir por $3.117$ salen unos $2.02L$. Dos rollos contiguos llenan esa longitud, así que un rollo mide cerca de $L$. "Cuadrado" es esa división. El modo que se enciende al menor Ra tiene esa longitud de onda.

Cuánto calor extra sale lo mide el número de Nusselt (transferencia total sobre conducción pura).

$q_\text{total}$ es el flujo real y $q_\text{conduction}$ el flujo sin convección. Solo conducción significa $\mathrm{Nu}=1$. Cuando la convección ya es fuerte, los experimentos quedan cerca de $\mathrm{Nu} \sim \mathrm{Ra}^{1/3}$. Leído solo ese exponente, un Ra ocho veces mayor duplica Nu. Justo en el umbral no aplica. Ahí Nu sigue cerca de 1.

La convección que Lorenz redujo a tres ecuaciones#

En 1963 Edward Lorenz dejó solo los modos más gruesos de temperatura y flujo. Un truncamiento de Galerkin. Las ecuaciones en derivadas parciales caen a tres ordinarias.

$x$ es la intensidad de circulación del rollo, $y$ el contraste de temperatura entre columnas que suben y bajan, $z$ la distorsión del perfil vertical de temperatura. $\sigma$ es el número de Prandtl (difusividad de momento sobre difusividad térmica), $\beta$ una constante geométrica y $\rho$ es $\mathrm{Ra}/\mathrm{Ra}_c$.

Subir $\rho$ es calentar la olla con más fuerza. Para $\rho<1$ toda trayectoria va al origen. Conducción. Para $1<\rho<24.74$ se asienta en un punto fijo fuera del origen. Un rollo estacionario. Pasado $\rho>24.74$ no se queda en ningún punto fijo y deambula entre dos alas. Caos.

Un estado estacionario es donde las tres derivadas valen cero. Si $\sigma \neq 0$, entonces $y=x$. Si $x \neq 0$, entonces $z=\rho-1$. La ecuación que queda es $x^2 = \beta(\rho-1)$. Con el $\beta = 8/3$ del código y $\rho = 15$, $z=14$ y $x^2 = 112/3 \approx 37.3$. La fila impresa del medio, steady roll, significa que las últimas 4000 muestras de $x$ abarcan menos de 0.1. Con $\rho=28$, $z=27$ y $x^2=72$, así que las dos alas están ahí. El chaos impreso significa que esa cola no pasó la puerta de 0.1.

Python — integrando las tres caras de la convección#

Runge–Kutta de cuarto orden. $\sigma=10$, $\beta=8/3$, paso 0.01, 20000 pasos. La cola son las últimas 4000 muestras.

import numpy as np
 
SIGMA, BETA = 10.0, 8.0 / 3.0
 
def lorenz_deriv(state, rho):
    x, y, z = state
    return np.array([
        SIGMA * (y - x),        # cambio en la intensidad de circulacion del rollo
        x * (rho - z) - y,      # cambio en el contraste ascenso/descenso
        x * y - BETA * z,       # cambio en la distorsion vertical de temperatura
    ])
 
def integrate_rk4(rho, dt=0.01, n=20000, s0=(0.6, 0.4, 0.2)):
    s = np.array(s0, dtype=float)
    traj = np.empty((n, 3))
    for i in range(n):
        k1 = lorenz_deriv(s, rho)
        k2 = lorenz_deriv(s + 0.5 * dt * k1, rho)
        k3 = lorenz_deriv(s + 0.5 * dt * k2, rho)
        k4 = lorenz_deriv(s + dt * k3, rho)
        s = s + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
        traj[i] = s
    return traj
 
def convection_regime(rho):
    tail = integrate_rk4(rho)[-4000:]          # cola tras descartar el transitorio
    spread = tail[:, 0].max() - tail[:, 0].min()
    if abs(tail[-1, 0]) < 1e-3 and rho < 1.0:
        return "conduction"                     # converge al origen
    if spread < 0.1:
        return "steady roll"                    # converge a un punto fijo
    return "chaos"                              # oscila entre dos alas
 
for rho in (0.8, 15.0, 28.0):
    print(f"rho={rho:5.1f} -> {convection_regime(rho)}")

Vuelven tres filas.

rho=  0.8 -> conduction
rho= 15.0 -> steady roll
rho= 28.0 -> chaos

Con el deslizador de ρ por debajo de 24.74 la trayectoria cae en un punto fijo amarillo. Al cruzarlo se abren las alas de mariposa. Reset desplaza un poco el inicio. Dos caminos que al principio se superponen se separan enseguida.

Dar por terminado el trabajo al cruzar 1708 es un fallo#

Doblar el espesor multiplica Ra por ocho, por $L^3$. Una capa que acababa de tocar 1708 queda entonces en $\rho = \mathrm{Ra}/\mathrm{Ra}_c = 8$. La ventana de rollo estacionario de Lorenz va de 1 a 24.74. Ocho está dentro. El caos quiere ρ en 24.74, lo que a la misma diferencia de temperatura exige un espesor de $24.74^{1/3} \approx 2.9$ veces. Llamar caos al rollo recién pasado el umbral se equivoca por ese trecho.

Ese trecho se construye subiendo $\Delta T$ con $\beta$ fijo. Cuando la diferencia de temperatura pasa de 30 K, ¿falla antes esa constante? Donde Boussinesq se rompe a partir de esa temperatura muestra qué número, 1708 o 24.74, pierde antes el sentido.