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

Llamar a la entropía "información perdida" no explica por qué este cálculo no revienta. En el método de Lattice Boltzmann (Lattice Boltzmann Method, LBM), bajar la viscosidad empuja $\tau$ hacia $1/2$, y BGK se pasa del equilibrio. Las poblaciones caen a negativo. El LBM entrópico (Entropic LBM, ELBM) de la línea de Karlin ata ahí la H discreta a su valor previo a la colisión. El precio es una raíz no lineal por nodo.

Al mover a la derecha el deslizador de la perturbación de no equilibrio, el punto amarillo se separa de la línea punteada $\alpha=2$ y la segunda raíz de $G(\alpha)=0$ se desplaza a la izquierda.

Con una perturbación pequeña el punto queda sobre 2. El hueco entre la línea punteada y el punto amarillo es el exceso que BGK no suelta.

Cuando τ se acerca a 0.5, BGK se desmorona#

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

$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)$. $c_s^2 = 1/3$ es el cuadrado de la velocidad reticular del sonido. Un flujo de alto Reynolds, o sea de baja viscosidad, obliga a $\tau \to 1/2$, es decir $\omega \to 2$.

Cuando $\omega$ se acerca a 2, la colisión aterriza al otro lado del equilibrio. Una población que se vuelve negativa 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 las colisiones solas nunca aumentan la cantidad siguiente.

$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 baje de forma monótona es el mismo enunciado que el crecimiento monótono de la entropía. Irreversibilidad. La flecha del tiempo. La desigualdad es un certificado de estabilidad. Un esquema en el que H no puede crecer tampoco puede dejar que la amplitud se dispare, porque una amplitud sin cota tendría que arrastrar a H.

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

LBM parte la velocidad continua en nueve direcciones (D2Q9). H se discretiza en consonancia.

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

El equilibrio entrópico es la población que minimiza esta $H$ bajo los momentos conservados. La masa $\sum_i f_i = \rho$ y el momento $\sum_i \mathbf{c}_i f_i = \rho\mathbf{u}$ quedan fijos mientras H baja. Un problema de multiplicadores de Lagrange. En D2Q9 existe una forma producto cerrada.

$u_\alpha$ es una componente de velocidad y $c_{i\alpha}$ una componente de velocidad reticular (−1, 0, 1). El truncamiento de segundo orden del LBM estándar es un equilibrio polinómico cortado en el segundo orden de Hermite, y ese polinomio no es el mínimo de H. El producto de arriba sí lo es.

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

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

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

$\alpha$ no es una constante. En cada paso y en cada nodo es la raíz no trivial de

El paso avanza hacia el equilibrio solo lo justo para dejar H idéntica a su valor previo a la colisión. No queda sitio para que H crezca. $\alpha=0$ es la raíz trivial. La raíz necesaria está más allá del equilibrio. Cerca del equilibrio $\alpha \approx 2$, que es otra vez BGK. Solo una desviación grande aparta $\alpha$ de 2 y recorta la sobre-relajación.

La fórmula de la viscosidad deja de ser abstracta con números. Con $c_s^2 = 1/3$ y $\nu = c_s^2(\tau - 1/2)$, $\tau = 0.51$ deja un hueco de $0.01$, así que $\nu = 1/300 \approx 0.00333$. La tasa es $\omega = 1/0.51 \approx 1.961$. Una desviación de longitud 1 se recorre 1.961 veces y cae 0.961 más allá del equilibrio. En $\tau = 0.501$, $\omega \approx 1.996$. El hueco que queda hasta el valor espejo $\omega = 2$ es 0.004.

El $\beta = 0.95$ del listado es la misma identidad. $\tfrac{1}{2\beta}-\tfrac{1}{2} = \tfrac{1}{1.9}-\tfrac{1}{2} = \tfrac{1}{38}$, de modo que $\nu = 1/114 \approx 0.00877$. BGK con esa viscosidad tiene $\tau = 1/(2\beta) \approx 0.526$. Cerca del equilibrio, $\alpha \approx 2$ da una tasa efectiva $\alpha\beta = 1.9$ y la viscosidad coincide con ese paso BGK. Si una desviación grande baja $\alpha$ de 2, $\alpha\beta$ cae bajo 1.9 y el exceso $\alpha\beta-1$ se acorta por debajo de 0.9. La viscosidad se queda donde la puso $\beta$. Solo se recorta el exceso.

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

Equilibrio en forma producto, la función H y una bisección para $\alpha$.

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 $\alpha$ impresa queda cerca de 2 y no es exactamente 2. Ese hueco es la corrección que mantiene H fija.

Las dos líneas de la perturbación son un cambio de momentos. Sumar $0.03$ a $f_1$ y restar $0.03$ de $f_3$ deja intacta la suma de masa. Las velocidades reticulares de esos enlaces son $+1$ y $-1$, así que el momento en $x$ sube $0.03-(-0.03)=0.06$. Desde $\rho=1$, $u_x=0.05$, el estado justo tras el golpe tiene $u_x=0.11$. El factor producto es $(2u+\sqrt{1+3u^2})/(1-u)$: cerca de 1.162 en $u=0.05$, cerca de 1.391 en $u=0.11$, y cerca de 2.467 a velocidad adimensional 0.3, donde $1-u$ vale solo 0.7. La banda $|u| \lesssim 0.1$ es donde ese factor todavía es suave.

El tope de la bisección es una restricción de signo. Una componente con $d_i$ negativo llega a cero en $\alpha=-f_i/d_i$, y el menor de esos valores es $a_{\max}$. La búsqueda se corta en $\min(0.999\,a_{\max},\,4)$. Una raíz más allá de 4, o una que cambia el signo de una población, no entra como candidata.

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

La curva de arriba es esta $G$. El punto amarillo es la raíz no trivial $\alpha^*$. Una perturbación mayor la despega de 2. BGK clava ese sitio en 2. Si la diferencia se acumula paso a paso, la retícula diverge.

Interactivo — bloquear el instante en que H crece#

Al pasar el mismo estado de no equilibrio por las dos reglas de colisión, solo las poblaciones que se volvieron negativas pasan a barras rojas.

El lado entrópico se queda en $H_{ent} = H_{pre}$. Así está construido. En el lado BGK, en cuanto crece la desviación, $H_{BGK}$ cruza la H previa a la colisión. Ahí se enciende la alerta roja. Esa subida de H es donde una retícula de baja viscosidad empieza a romperse.

Tres cosas con las que choqué al implementar#

Una resolución no lineal en cada nodo pesa antes que el término de viscosidad. Cerca del equilibrio, cuando $\|f - f^{eq}\|$ es pequeño, se fija $\alpha \approx 2$ y la bisección se enciende solo donde la desviación cruza un umbral. El costo baja a cerca de una décima parte.

El equilibrio producto solo está definido para $|u| < 1$. El denominador lleva $1 - u$. Pasada una velocidad adimensional de 0.3, el tensor de presión empieza a torcerse según la dirección. Incluso ELBM queda más seguro dentro de $|u| \lesssim 0.1$.

Si el equilibrio sigue siendo el polinomio truncado de segundo orden, desaparece la garantía de mínimo de H. El estabilizador rinde del todo solo cuando el propio equilibrio es el producto entrópico. Atornillar $\alpha$ sobre un polinomio deja la mitad.

Subir solo β con α clavada en 2 suelta la cota de H#

Saltar a $\alpha \approx 2$ cerca del equilibrio es un recorte de costo para $\|f - f^{eq}\|$ pequeño. Empujar $\beta$ hacia 1 y congelar ese atajo en todos los nodos deja la tasa efectiva $\alpha\beta$ otra vez sobre 2. Vuelve el exceso hasta la imagen especular, y la H discreta puede superar su valor previo a la colisión. Bajó la viscosidad. La raíz no se resolvió.

Una colisión puede mantener H plana y llegar tarde igual, si el streaming importa una población que ya es negativa en un nodo vecino. Cuando el signo se rompe en la advección, queda abierto si aún hace falta un limitador de flujo que preserva la positividad o si la raíz entrópica alcanza por sí sola.