Sobrevivir al régimen de baja viscosidad donde BGK explota — Colisión MRT en Lattice Boltzmann

Dos modelos de colisión relajan hacia el mismo equilibrio en el mismo retículo, y aun así uno gira con suavidad mientras el otro se rompe en un patrón de tablero de ajedrez. La diferencia es una sola elección. BGK relaja las nueve funciones de distribución con un único tiempo de relajación, mientras que MRT (Multiple Relaxation Time, tiempo de relajación múltiple) las transforma en momentos y relaja cada uno a su propio ritmo. Este artículo rastrea por qué esa elección decide la estabilidad, usando el método de lattice Boltzmann D2Q9. Movemos las distribuciones a momentos con una matriz de transformación $M$, ajustamos por separado los grados de libertad que nada tienen que ver con la viscosidad, y observamos con un vórtice de Taylor–Green cómo MRT sostiene el régimen de baja viscosidad donde BGK explota.

La debilidad de relajarlo todo en un solo compás#

El método de lattice Boltzmann (LBM) repite dos pasos: las funciones de distribución $f_i$ se desplazan por el retículo (streaming) y, en cada punto, relajan hacia el equilibrio (collision). La colisión BGK más común se escribe como

donde $f_i$ es la distribución a lo largo de la velocidad $\mathbf{c}_i$, $\tau$ es el tiempo de relajación y $f_i^{\text{eq}}$ es el equilibrio local. La viscosidad queda fijada solo por $\tau$.

El problema empieza aquí. Para bajar la viscosidad (subir el número de Reynolds) hay que acercar $\tau$ a $1/2$. Pero cuando $\tau \to 1/2$, los modos de alto orden que nada tienen que ver con la viscosidad apenas relajan y recorren el retículo. Cuando esos modos crecen cerca del límite de resolución de la malla, surge una oscilación de tablero y la simulación diverge. En BGK, la perilla de la viscosidad y la de la estabilidad están atadas al mismo mango.

Mirar la colisión en el espacio de momentos#

La idea central es simple. En lugar de relajar las nueve funciones de distribución directamente, primero se las transforma en momentos con significado físico y luego se relaja. Los momentos son combinaciones lineales de las distribuciones.

Aquí $\mathbf{f} = (f_0,\dots,f_8)^\top$ es el vector de distribuciones, $M$ es la matriz de transformación $9\times 9$ y $\mathbf{m}$ es el vector de momentos. La colisión relaja hacia los momentos de equilibrio $\mathbf{m}^{\text{eq}}$ en el espacio de momentos y luego regresa al espacio de velocidades.

$S = \mathrm{diag}(s_0,\dots,s_8)$ es la matriz diagonal con las tasas de relajación por momento. Si todos los $s_i$ se fijan al mismo valor $1/\tau$, esta expresión vuelve exactamente a BGK. MRT no hace más que resolver esas entradas diagonales momento por momento.

Los nueve momentos de D2Q9 y la matriz de transformación M#

En el retículo D2Q9 (dos dimensiones, nueve velocidades) los momentos son estos nueve: densidad $\rho$, energía $e$, energía al cuadrado $\varepsilon$, momento $j_x, j_y$, flujo de energía $q_x, q_y$ y los tipo tensión $p_{xx}, p_{xy}$.

Entre ellos, $\rho, j_x, j_y$ son cantidades conservadas que la colisión no altera, así que sus tasas de relajación se fijan en $s=0$. Los otros seis son las perillas ajustables. La matriz de transformación estándar de Lallemand–Luo se construye con filas ortogonales, por lo que la inversa es simplemente la transpuesta dividida por las normas de fila.

Los momentos de equilibrio cierran solo con densidad y momento.

Aquí $j_x = \rho u_x$, $j_y = \rho u_y$ son las componentes del momento. El equilibrio de los momentos de tensión $p_{xx}, p_{xy}$ se acopla a los gradientes de velocidad, así que sus tasas de relajación fijan la viscosidad.

Una tasa distinta por momento — grados de libertad ajenos a la viscosidad#

Quien fija la viscosidad de corte es la tasa de relajación de los momentos de tensión, $s_\nu$.

La clave es que las demás tasas $s_e$ (volumétrica), $s_\varepsilon$ y $s_q$ (flujo de energía) no afectan en nada a la viscosidad. Son parámetros libres. Aun empujando $s_\nu$ cerca de $2$ para baja viscosidad, los modos de alto orden pueden mantenerse aparte en la banda estable de $1.5\sim1.8$. Las dos perillas que BGK mantenía atadas se separan.

Ajusta las tasas de relajación directamente en las barras de abajo.

Sube $s_\nu$ hasta $1.99$ y la viscosidad $\nu$ cae casi a $0$, pero aún así puedes mantener $s_e, s_q$ en la banda estable. Observa también que las barras de los momentos conservados ($\rho, j_x, j_y$) quedan ancladas en cero.

Colisión en el espacio de momentos, streaming en el de velocidades#

El flujo del algoritmo es el siguiente.

1. Calcular las magnitudes macroscópicas $\rho, u_x, u_y$ a partir de $f_i$
2. Transformar a momentos con $\mathbf{m} = M\mathbf{f}$
3. Calcular los momentos de equilibrio $\mathbf{m}^{\text{eq}}$ (a partir de densidad y momento)
4. Relajar por momento: $m_a \leftarrow m_a - s_a(m_a - m_a^{\text{eq}})$
5. Restaurar al espacio de velocidades con $\mathbf{f}^{*} = M^{-1}\mathbf{m}$
6. Streaming: mover $f_i^{*}$ al vecino a lo largo de $\mathbf{c}_i$
7. Aplicar condiciones de contorno y volver al paso 1

Solo la colisión ocurre en el espacio de momentos; el streaming transcurre en el espacio de velocidades como siempre. El costo adicional es apenas dos productos de matrices $9\times 9$ por celda.

Python — Medir el decaimiento viscoso con un vórtice de Taylor–Green#

Para la verificación usamos un vórtice de Taylor–Green. Este flujo tiene una solución exacta de decaimiento viscoso, así que se puede contrastar la implementación con números. La energía cinética debe decaer como $\exp(-4\nu k^2 t)$, con $k$ el número de onda del vórtice.

import numpy as np
 
# Retículo D2Q9
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
 
# Matriz de transformacion Lallemand-Luo (filas: rho,e,eps,jx,qx,jy,qy,pxx,pxy)
M = np.array([
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [-4,-1,-1,-1,-1, 2, 2, 2, 2],
    [ 4,-2,-2,-2,-2, 1, 1, 1, 1],
    [ 0, 1, 0,-1, 0, 1,-1,-1, 1],
    [ 0,-2, 0, 2, 0, 1,-1,-1, 1],
    [ 0, 0, 1, 0,-1, 1, 1,-1,-1],
    [ 0, 0,-2, 0, 2, 1, 1,-1,-1],
    [ 0, 1,-1, 1,-1, 0, 0, 0, 0],
    [ 0, 0, 0, 0, 0, 1,-1, 1,-1],
], dtype=float)
Minv = np.linalg.inv(M)
 
def d2q9_equilibrium(rho, ux, uy):
    eq = np.zeros((9,) + rho.shape)
    for i in range(9):
        cu = cx[i] * ux + cy[i] * uy
        eq[i] = w[i] * rho * (1 + 3*cu + 4.5*cu**2 - 1.5*(ux**2 + uy**2))
    return eq
 
def relaxation_rates(nu):
    s_nu = 1.0 / (3*nu + 0.5)          # la tasa que fija la viscosidad
    # conservados (rho,jx,jy)=0, modos de alto orden fijos en la banda estable
    return np.array([0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu])
 
def taylor_green_init(N, U0):
    x = np.arange(N)
    X, Y = np.meshgrid(x, x, indexing='ij')
    k = 2*np.pi / N
    ux = -U0 * np.cos(k*X) * np.sin(k*Y)
    uy =  U0 * np.sin(k*X) * np.cos(k*Y)
    rho = 1 - 0.75 * U0**2 * (np.cos(2*k*X) + np.cos(2*k*Y))
    return d2q9_equilibrium(rho, ux, uy), k
 
def macros(f):
    rho = f.sum(axis=0)
    ux = (cx[:, None, None] * f).sum(axis=0) / rho
    uy = (cy[:, None, None] * f).sum(axis=0) / rho
    return rho, ux, uy
 
def mrt_collide(f, s):
    rho, ux, uy = macros(f)
    jx, jy = rho*ux, rho*uy
    m = np.tensordot(M, f, axes=([1], [0]))        # al espacio de momentos
    meq = np.zeros_like(m)
    meq[0] = rho
    meq[1] = -2*rho + 3*(jx**2 + jy**2)
    meq[2] = rho - 3*(jx**2 + jy**2)
    meq[3], meq[4] = jx, -jx
    meq[5], meq[6] = jy, -jy
    meq[7] = jx**2 - jy**2
    meq[8] = jx*jy
    m -= s[:, None, None] * (m - meq)               # relajacion por momento
    return np.tensordot(Minv, m, axes=([1], [0]))    # de vuelta al espacio de velocidades
 
def stream(f):
    for i in range(9):
        f[i] = np.roll(f[i], (cx[i], cy[i]), axis=(0, 1))
    return f
 
def run_mrt_lbm(N=64, nu=0.004, U0=0.04, steps=2000):
    f, k = taylor_green_init(N, U0)
    s = relaxation_rates(nu)
    e0 = None
    for t in range(steps):
        f = mrt_collide(f, s)
        f = stream(f)
        if t % 400 == 0:
            _, ux, uy = macros(f)
            E = np.mean(ux**2 + uy**2)
            e0 = E if e0 is None else e0
            print(f"t={t:5d}  KE/KE0={E/e0:.4f}  analytic={np.exp(-4*nu*k**2*t):.4f}")
 
run_mrt_lbm()

La salida muestra que el decaimiento numérico coincide con la curva analítica $\exp(-4\nu k^2 t)$ hasta dos decimales. Al correr BGK en el mismo retículo y la misma $\nu$, diverge primero a una viscosidad más baja.

BGK y MRT lado a lado al mismo Re#

En la simulación de abajo, baja el deslizador de viscosidad y cambia el modelo de colisión con el botón BGK/MRT.

Baja la viscosidad hasta cerca de $0.002$ y déjalo en BGK: un patrón rugoso invade el retículo y aparece "diverged". Cambia a MRT con la misma viscosidad y haz reset, y el vórtice decae con suavidad. Aunque ambos modelos apuntan al mismo número de Reynolds, sus límites de estabilidad se distinguen de un vistazo.

La próxima vez que toques lattice Boltzmann#

MRT no es magia. Solo separa físicamente la perilla que fija la viscosidad ($s_\nu$) de las perillas que fijan la estabilidad ($s_e, s_q$). Esa separación amplía el rango de operación del lattice Boltzmann a baja viscosidad y alto Reynolds.

Quedan tres cosas. Primera: BGK explota a baja viscosidad porque los modos de alto orden ajenos a la viscosidad se frenan junto con ella. Segunda: MRT mueve solo la colisión al espacio de momentos para ajustar esos modos aparte, a un costo de dos productos de matrices por celda. Tercera: siempre puedes verificar la implementación contra el decaimiento $\exp(-4\nu k^2 t)$ de un vórtice de Taylor–Green.