La entropía frena la explosión — el Lattice Boltzmann entrópico de Karlin

La entropía suele llegar etiquetada como "información perdida": la dirección en que crece el desorden, una medida de pérdida irreversible. Y sin embargo, en el método de Lattice Boltzmann (Lattice Boltzmann Method, LBM), es precisamente esa cantidad la que atrapa una simulación desbocada y la deja quieta. Este artículo trata de cómo el LBM entrópico (Entropic LBM, ELBM) de la escuela de Karlin impone el teorema H a nivel de la discretización para eliminar la divergencia en el régimen de baja viscosidad. Al terminar entenderás por qué vale la pena resolver una ecuación no lineal en cada paso por una sola línea del operador de colisión.

Cuando τ se acerca a 0.5, BGK se desmorona#

El LBM estándar usa la colisión BGK (aproximación de tiempo de relajación único). Atrae la distribución $f_i$ hacia el equilibrio $f_i^{eq}$ por un factor $\omega$ .

f_i^* = f_i + \omega\,(f_i^{eq} - f_i)

Aquí $f_i$ es la población en la dirección reticular $i$ , $f_i^{eq}$ es el equilibrio y $\omega = 1/\tau$ es la tasa de relajación.

La viscosidad la fija $\nu = c_s^2(\tau - 1/2)$ , donde $c_s^2 = 1/3$ es el cuadrado de la velocidad reticular del sonido. Para resolver flujo de alto Reynolds (baja viscosidad), se empuja $\tau \to 1/2$ , es decir $\omega \to 2$ .

El problema es que, cuando $\omega$ se acerca a 2, la colisión se pasa del equilibrio y aterriza en el lado opuesto. Las poblaciones caen a negativo, y una vez que un $f_i$ se vuelve negativo, oscila con más fuerza en el paso siguiente. La retícula explota. BGK no tiene nada que detenga esta sobre-relajación (over-relaxation).

El teorema H de Boltzmann — la flecha del tiempo#

En 1872, Boltzmann demostró que la siguiente cantidad nunca crece bajo colisiones puras.

H = \int f \ln f \, d\mathbf{c}

$f$ es la distribución de velocidades y $\mathbf{c}$ es la velocidad molecular. Se cumple $dH/dt \le 0$ , con igualdad solo en el equilibrio (la distribución de Maxwell-Boltzmann).

H es entropía negativa. Que H decrezca de forma monótona es lo mismo que decir que la entropía crece de forma monótona. Eso es la irreversibilidad: la flecha del tiempo. El punto clave es que esta desigualdad es un certificado de estabilidad. Un esquema numérico en el que H nunca crece no puede divergir, porque un crecimiento ilimitado de la amplitud exigiría que H aumentara.

La función H discreta y el equilibrio entrópico#

LBM rebana la velocidad continua en 9 direcciones (D2Q9). H se discretiza en consonancia.

H(f) = \sum_i f_i \ln\frac{f_i}{w_i}

$w_i$ son los pesos reticulares (en D2Q9: centro $4/9$ , caras $1/9$ , diagonales $1/36$ ). Es el logaritmo del cociente normalizado por el peso.

El equilibrio entrópico se define como la distribución que minimiza este $H$ bajo las restricciones de momentos conservados. Manteniendo fijas la masa $\sum_i f_i = \rho$ y el momento $\sum_i \mathbf{c}_i f_i = \rho\mathbf{u}$ , se minimiza H (un problema de multiplicadores de Lagrange). En D2Q9 existe una forma producto cerrada.

f_i^{eq} = \rho\, w_i \prod_{\alpha=x,y}\left(2 - \sqrt{1+3u_\alpha^2}\right)\left(\frac{2u_\alpha + \sqrt{1+3u_\alpha^2}}{1 - u_\alpha}\right)^{c_{i\alpha}}

$u_\alpha$ es una componente de velocidad y $c_{i\alpha}$ una componente de velocidad reticular (−1, 0, 1). A diferencia del equilibrio polinómico truncado de segundo orden del LBM estándar, este es el verdadero mínimo de H.

El estabilizador entrópico — resolver α en cada paso#

La colisión de ELBM tiene la misma forma que BGK, pero parte la tasa de relajación en dos.

f_i^* = f_i + \alpha\beta\,(f_i^{eq} - f_i)

$\beta \in (0,1)$ fija la viscosidad ( $\nu = c_s^2(\tfrac{1}{2\beta} - \tfrac{1}{2})$ , con $\beta\to1$ no viscoso) y $\alpha$ fija la estabilidad.

$\alpha$ no es una constante. En cada paso, en cada nodo de la retícula, se resuelve como la raíz no trivial de

H\big(f + \alpha(f^{eq} - f)\big) = H(f)

El significado es claro: relajar hacia el equilibrio exactamente lo que deja H idéntica a su valor previo a la colisión. Cualquier aumento de H queda estructuralmente prohibido. $\alpha=0$ es la raíz trivial; lo que buscamos es la segunda raíz más allá del equilibrio. Cerca del equilibrio $\alpha \approx 2$ , que es exactamente BGK. Solo cuando la parte de no equilibrio es grande, $\alpha$ se aparta de 2 y recorta la sobre-relajación.

Python — resolver la α entrópica por bisección#

Aquí están el equilibrio en forma producto, la función H y un paso que resuelve $\alpha$ por bisección.

import numpy as np
 
# Pesos y velocidades de la retícula D2Q9 (lattice weights & velocities)
W  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CX = np.array([0, 1, 0, -1,  0, 1, -1, -1,  1])
CY = np.array([0, 0, 1,  0, -1, 1,  1, -1, -1])
 
def h_function(f):
    # función H discreta: H = Σ f_i ln(f_i / w_i)  (versión reticular de la H de Boltzmann)
    return np.sum(f * np.log(f / W))
 
def entropic_equilibrium(rho, ux, uy):
    # equilibrio entrópico exacto en forma producto de Ansumali–Karlin
    def factor(u):
        s = np.sqrt(1 + 3 * u * u)
        return (2 - s), (2 * u + s) / (1 - u)
    bx, rx = factor(ux)
    by, ry = factor(uy)
    return rho * W * bx * by * (rx ** CX) * (ry ** CY)
 
def entropic_alpha(f, feq):
    # raíz no trivial α* de H(f + α(feq − f)) = H(f), por bisección
    d = feq - f
    H0 = h_function(f)
    G = lambda a: h_function(f + a * d) - H0
    amax = np.min(np.where(d < 0, -f / d, np.inf))   # cota de positividad
    hi = min(0.999 * amax, 4.0) if np.isfinite(amax) else 4.0
    if G(hi) <= 0:
        return 2.0                                   # valor seguro cerca del equilibrio
    lo = 1.0
    while G(lo) > 0 and lo > 1e-3:
        lo *= 0.5
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        if G(mid) > 0: hi = mid
        else:          lo = mid
    return 0.5 * (lo + hi)
 
# --- un nodo, un paso ---
feq = entropic_equilibrium(1.0, 0.05, 0.0)
f = feq.copy()
f[1] += 0.03; f[3] -= 0.03            # no equilibrio traído por el streaming (non-equilibrium)
rho = f.sum(); ux = (CX * f).sum() / rho; uy = (CY * f).sum() / rho
feq = entropic_equilibrium(rho, ux, uy)
 
alpha = entropic_alpha(f, feq)
beta  = 0.95                          # ν = c_s²(1/(2β) − 1/2), β→1 es baja viscosidad
f_post = f + alpha * beta * (feq - f)
 
print(f"α entrópica  = {alpha:.4f}   (valor fijo de BGK 2.0)")
print(f"H antes      = {h_function(f):.8f}")
print(f"H después    = {h_function(f_post):.8f}")
print(f"pob. mínima  = {f_post.min():.3e}  (> 0)")

La salida muestra una $\alpha$ cercana a 2 pero no exactamente 2. Ese pequeño hueco es la corrección que conserva H.

Interactivo — la segunda raíz de G(α)#

Manipula la simulación de abajo. Aumenta la perturbación de no equilibrio y observa cómo se curva la gráfica.

Non-equilibrium perturbation: 0.060

The yellow dot is the non-trivial root α* where the discrete H-function returns to its pre-collision value. Near equilibrium α* ≈ 2, which is exactly standard BGK. Push the perturbation up and α* drifts away from 2 — that gap is the entropic correction BGK ignores.

El punto amarillo es la $\alpha^*$ que buscamos. Para perturbaciones pequeñas se pega justo sobre la línea punteada ( $\alpha=2$ ). Empuja el deslizador a la derecha y $\alpha^*$ se despega de 2. BGK ignora este desplazamiento y siempre usa 2. Cuando el hueco se acumula, la retícula diverge.

Interactivo — bloquear el instante en que H crece#

El mismo estado de no equilibrio, procesado por dos reglas de colisión en paralelo. Las barras rojas son poblaciones que se volvieron negativas.

Non-equilibrium shear: 0.070

Both panels show post-collision populations. The entropic collision lands exactly on the equal-entropy point, so H_ent = H_pre. Fixed-α BGK misses it; once the non-equilibrium part grows, BGK over-shoots and H actually increases — a violated discrete H-theorem, which is the seed of low-viscosity blow-up.

El lado entrópico siempre cumple $H_{ent} = H_{pre}$ : así está construido. En el lado BGK, en cuanto la parte de no equilibrio crece, llega un instante en que $H_{BGK}$ supera al valor previo a la colisión. Ahí se enciende la alerta roja. Ese aumento de H es justo la semilla que revienta la retícula a baja viscosidad.

Tres cosas con las que choqué al implementar#

Primera, una resolución no lineal en cada nodo es cara. Por eso, en la práctica se fija $\alpha \approx 2$ cerca del equilibrio (cuando $\|f - f^{eq}\|$ es pequeño) y se activa la bisección solo en los nodos donde la parte de no equilibrio cruza un umbral. El costo cae a una décima parte aproximadamente.

Segunda, el equilibrio en forma producto solo está definido para $|u| < 1$ : hay un $1 - u$ en el denominador. Pasada una velocidad adimensional de 0.3, el tensor de presión empieza a torcerse de forma anisótropa, así que incluso ELBM es más seguro en el régimen $|u| \lesssim 0.1$ .

Tercera, si conservas el equilibrio polinómico truncado de segundo orden, la garantía de minimización de H se rompe. Para obtener el beneficio completo del estabilizador entrópico, el propio equilibrio debe ser la forma producto entrópica. Atornillar $\alpha$ sobre un equilibrio polinómico te da solo la mitad del efecto.

Resumen en tres líneas#

BGK explota a baja viscosidad porque la sobre-relajación aumenta la H discreta. ELBM lo evita resolviendo, en cada paso, la $\alpha$ que hace que H tras la colisión sea igual a H antes de ella.
$\alpha$ es la raíz no trivial de $H(f + \alpha(f^{eq}-f)) = H(f)$ , y cerca del equilibrio $\alpha \approx 2$ recupera BGK automáticamente.
La viscosidad la lleva $\beta$ y la estabilidad la lleva $\alpha$ . El costo es una resolución no lineal por paso, pero el atajo cerca del equilibrio lo mantiene asequible en la práctica.