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

Mismo retículo, mismo equilibrio. Uno gira. El otro se rompe en un tablero. BGK empuja las nueve distribuciones con un único tiempo de relajación. MRT (Multiple Relaxation Time, tiempo de relajación múltiple) las pasa a momentos y empuja cada uno a su propio ritmo.

Al llevar la barra $s_\nu$ hasta 1.99, el número de viscosidad cae hacia 0. Las barras $s_e$ y $s_q$ se quedan donde estaban.

Las barras conservadas $\rho$, $j_x$, $j_y$ siguen en cero. BGK ata todas las barras a un solo número. Ese nudo es el sitio del tablero a baja viscosidad.

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#

No se relajan las nueve distribuciones tal como están. Primero se transforman en momentos con significado físico y después 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$.

Las otras tasas $s_e$ (volumétrica), $s_\varepsilon$ y $s_q$ (flujo de energía) no entran en la viscosidad. Son libres. Aunque $s_\nu$ se pegue a $2$, los modos de alto orden pueden quedarse en la banda estable $1.5\sim1.8$. Lo que BGK ataba a una sola perilla se separa aquí.

La viscosidad por defecto del listado es $\nu=0.004$. Entonces $s_\nu = 1/(3\nu + 1/2)$, de modo que $3\times 0.004 + 0.5 = 0.512$ y $s_\nu = 1/0.512 = 1.953125$. Al llevar la barra a 1.99, $\nu = (1/3)(1/1.99 - 1/2)$. Restar 0.5 de $1/1.99 \approx 0.5025$ y dividir por 3 deja $\nu \approx 8.4\times 10^{-4}$. Cerca de cero, no cero. A esa viscosidad BGK también empuja $s_e$, $s_\varepsilon$ y $s_q$ hasta 1.99. El MRT del listado fija esos tres sitios en $1.64$, $1.54$ y $1.7$. El vector que corre es [0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu], así que los únicos ceros son los tres sitios conservados.

Los momentos de equilibrio se cierran solo con densidad y momento porque la distribución ya se cortó a orden bajo. Ese corte está en Maxwell–Boltzmann dentro de nueve velocidades.

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()

El valor analítico que el print pone al lado de la energía cinética es $\exp(-4\nu k^2 t)$. Los valores por defecto son $N=64$, $\nu=0.004$, $k=2\pi/N$, así que $4\nu k^2 = 0.016\,(\pi/32)^2 \approx 1.54\times 10^{-4}$. En $t=400$ el exponente vale cerca de $-0.0617$ y el factor cerca de 0.940. En $t=1600$ queda cerca de $-0.247$ y el factor cerca de 0.781. Con range(2000) el último $t$ impreso es 1600. Coincidir a dos decimales es comparar ese factor con la energía cinética numérica. BGK en el mismo retículo y la misma $\nu$ diverge antes, a una viscosidad más baja.

BGK y MRT lado a lado al mismo Re#

El control de viscosidad frena el decaimiento. El botón BGK/MRT decide si el retículo aguanta con esa misma $\nu$.

Al bajar la viscosidad hasta cerca de $0.002$ y dejar BGK, un patrón rugoso invade el retículo y aparece "diverged". Al pasar a MRT con la misma viscosidad y reiniciar, el vórtice decae con suavidad. Con $\nu=0.002$, $s_\nu=1/(3\times 0.002+0.5)=1/0.506\approx 1.976$. BGK empuja también los modos altos a ese número. Los $1.64$, $1.54$ y $1.7$ del listado se quedan al lado.

Dejar los modos fantasma cuando s_ν se acerca a 2#

MRT solo separa la perilla de la viscosidad $s_\nu$ de las perillas de estabilidad $s_e$ y $s_q$. Copiar $s_\nu=1.953$ de $\nu=0.004$ a todos los modos es BGK. Dejar $1.64$, $1.54$ y $1.7$ mientras solo sube $s_\nu$ es MRT. El costo extra son dos productos $9\times 9$ por celda.

Sin esa separación, un Taylor–Green de baja viscosidad muere en tablero. Otro camino para frenar la misma escapada es un teorema H dentro de la colisión. La comprobación del decaimiento sigue siendo la línea $\exp(-4\nu k^2 t)$.