¿Adónde fue la mitad de la fuerza en LBM? — Cuatro esquemas de forcing y el 1−1/(2τ)

En 1993 Shan y Chen separaron el agua del aceite sobre una malla. La receta era modesta. Todo consistía en sumar $\tau\mathbf{F}/\rho$ a la velocidad de equilibrio que entra en la colisión. En 1998 He metió la fuerza directamente en la distribución de equilibrio. En 2002 Guo señaló que las recetas anteriores dejaban un error en el tensor de esfuerzos y antepuso un $(1-\frac{1}{2\tau})$ al término de fuerza. En 2004 Kupershtokh consiguió lo mismo con la sola diferencia entre dos distribuciones de equilibrio.

Cuatro respuestas para una misma pregunta. ¿Cuál conviene usar entonces? Hoy la pregunta se le hizo directamente a un canal D2Q9. La respuesta no fue la esperada. En flujo estacionario, cualquiera de los esquemas arrojó el mismo error hasta la cuarta cifra decimal. La verdadera trampa no estaba en el momento de elegir el esquema. Este post ubica esa trampa con la derivación y mide su tamaño con números.

Se recorta el término de fuerza de la ecuación continua#

La ecuación de Boltzmann con fuerza externa gana un término.

$f$ es la función de distribución, $\boldsymbol{\xi}$ la velocidad de partícula, $\mathbf{F}$ la fuerza por unidad de volumen y $\lambda$ el tiempo de relajación.

El problema está en el tercer término. La derivada en el espacio de velocidades $\nabla_{\boldsymbol{\xi}} f$ no existe sobre la malla. Lo único disponible son nueve velocidades discretas. Por eso se reemplaza $f$ por $f^{\rm eq}$ y la derivada se trata de forma analítica.

$c_s$ es la velocidad del sonido de la malla (en D2Q9, $c_s^2=1/3$) y $\mathbf{u}$ es la velocidad del fluido.

¿Vale sustituir $f$ por $f^{\rm eq}$? En bajo Mach, sí. La parte de no equilibrio $f^{(1)}$ es de primer orden en el número de Knudsen (recorrido libre medio / longitud característica), así que al multiplicarse por la fuerza cae a un orden despreciable. A cambio, esa sustitución deja más adelante un error residual en el tensor de esfuerzos. Ahí es justo donde los esquemas se separan.

Esta expresión se proyecta sobre las velocidades discretas $\mathbf{c}_i$ y la serie de Hermite se trunca en segundo orden. De ahí sale el término de fuerza sobre la malla.

$w_i$ es el peso de la malla. La razón de truncar justo en segundo orden es clara. Los momentos que hacen falta para recuperar Navier–Stokes llegan hasta el segundo.

Conviene dejar verificados los tres momentos.

El de orden cero es la masa. La fuerza no crea masa, así que debe anularse. El de primer orden es la cantidad de movimiento; el de segundo, el tensor de esfuerzos.

La mitad que dejó la integración trapezoidal#

Ahora se integra la ecuación de Boltzmann de velocidades discretas a lo largo de las características durante $\delta t$. Integrar el término de colisión y el de fuerza con la regla del trapecio da precisión de segundo orden. A cambio, en el lado derecho se cuela $t+\delta t$.

Los superíndices $n$ y $n+1$ denotan $t$ y $t+\delta t$. Así tal cual, el esquema es implícito. Habría que resolver un sistema de ecuaciones en cada paso.

LBM borra esa implicitud con un cambio de variable.

Al reescribir en términos de $\bar f_i$ y tomar $\tau = \lambda/\delta t + 1/2$, queda la línea de siempre.

El cambio dejó dos huellas. Una es el $(1-\frac{1}{2\tau})$ que antecede a la fuerza. La otra pasa más desapercibida. Lo que se guarda en el arreglo es $\bar f_i$, no $f_i$. Por eso la cantidad de movimiento también queda corrida en esa mitad.

Las dos son una sola cosa. Salieron juntas del mismo cambio de variable. Tomar una y olvidar la otra descuadra la contabilidad. Cuánto se descuadra se mide más abajo.

Tabla comparativa de los esquemas#

Sea $\mathbf{j}=\sum_i \bar f_i\mathbf{c}_i$ la cantidad de movimiento cruda de las distribuciones almacenadas, y $\delta\mathbf{u}=\delta t\,\mathbf{F}/\rho$ el incremento de velocidad que la fuerza aporta en un paso.

Leyendo las cinco filas en vertical aparece lo que tienen en común. Antes de que el coeficiente lo multiplique, el momento de primer orden del término de fuerza es exactamente $\mathbf{F}$ en todos. Es lo esperable, porque así se diseñaron. Lo que se separa es el momento de segundo orden, es decir, el tensor de esfuerzos.

Pero el valor objetivo que exige el cambio de variable no es $\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$ sino $(1-\frac{1}{2\tau})(\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u})$. El único que carga ese coeficiente de forma explícita es el de la última fila. Los demás arrastran una desviación de esfuerzos del tamaño de $\frac{1}{2\tau}$.

El tamaño de la desviación es el producto de fuerza por velocidad, o sea de orden $u F$. En bajo Mach, con $u$ pequeña, ese término ya es pequeño de por sí. El término $\mathbf{F}\mathbf{F}$ que dejan Shan–Chen y EDM es todavía menor, de orden $F^2$. Por eso, en problemas empujados por una gravedad débil, da casi igual cuál se use. La diferencia asoma donde la fuerza es intensa (cuando $\tau\mathbf{F}/\rho$ compite con $\mathbf{u}$) o donde la fuerza varía bruscamente en la interfaz, como en flujo multifásico.

Conviene manipular la figura de abajo.

Por más que se gire el dial de dirección de la fuerza, la barra M0 (masa) sigue pegada a cero. Una fuerza no fabrica masa, así que esa sale gratis. El objetivo de M1 (cantidad de movimiento) no es $\mathbf{F}$ sino $(1-\frac{1}{2\tau})\mathbf{F}$, porque lo que se mide es el término de fuerza ya multiplicado por el coeficiente. El $\mathbf{F}/(2\tau)$ que falta lo devuelve el equilibrio desplazado medio paso, y la comprobación es que la última fila, $\Delta j$, cae exactamente sobre $\mathbf{F}$.

M2 (esfuerzos) es la única fila donde los esquemas se separan. Al subir el deslizador de $|u|$, el M2 de plain se queda en cero y no alcanza el objetivo: lo equivocado es la forma, no el tamaño, así que ningún coeficiente lo arregla. Al empujar $\tau$ hacia 2, la desviación de Shan–Chen crece a ojos vistas. Y al apagar el interruptor del prefactor, M1, M2 y $\Delta j$ se rompen a la vez.

El Poiseuille estacionario no los distingue#

Pasemos al problema de verificación más común. Es un canal empujado por una fuerza volumétrica entre dos paredes. La solución estacionaria se conoce: una parábola.

$H$ es la altura del canal y $\nu$ la viscosidad cinemática. Con half-way bounce-back la pared queda media celda por fuera del nodo, así que $H$ coincide con el número de nodos de fluido y la coordenada del nodo $j$ es $y=j+0.5$.

Conviene manipular la simulación de abajo.

Basta con pulsar los cuatro botones de esquema. La curva turquesa no se mueve sobre la parábola punteada. Empujar $\tau$ de 0.6 a 2.0 tampoco cambia nada. Después basta con apagar uno solo de los dos interruptores de abajo. En ese instante la amplitud de la parábola cambia por completo.

Con números pasa lo mismo. Es el error relativo máximo medido tras correr 40 000 pasos con 33 nodos de fluido y $F_x = 10^{-5}$.

No hay nada que comparar entre columnas. Los valores son iguales. El error que queda no viene del esquema sino de la discretización del bounce-back.

La razón está en el lugar por donde entra el momento de segundo orden. Este flujo es estacionario y unidireccional. La velocidad se reduce a $u_x(y)$ y la fuerza solo tiene componente $x$. El término que separa a los esquemas es $\mathbf{u}\mathbf{F}$, es decir, la componente $xx$. Pero lo que gobierna de verdad el balance de cantidad de movimiento de este flujo es el esfuerzo cortante $xy$. La desviación $xx$ que genera la fuerza no tiene por dónde entrar en el balance.

También pesa que $\partial_t\mathbf{u}=0$. No hay margen para que el error del término de fuerza se acumule en el tiempo. Para ver la diferencia entre esquemas hay que ir a flujo no estacionario o a una fuerza que varíe en el espacio. Con esta prueba, al menos, no se pueden separar los cinco.

No es una mala noticia. Es información que conviene tener a la hora de elegir el problema de verificación. Significa que aprobar Poiseuille no permite afirmar que la implementación del forcing sea correcta.

Contar la mitad dos veces agranda la fuerza#

¿Qué es lo que sí se rompe entonces? Las dos cosas que antes se dijeron inseparables: la corrección de media fuerza en la velocidad y el coeficiente que antecede al término de fuerza.

Contemos la cantidad de movimiento que realmente se inyecta en un paso. Al sumar $\delta t\mathbf{F}/2\rho$ a la velocidad de equilibrio, la colisión transporta cada paso una fracción $1/\tau$ de eso. A ello se añade la parte que aporta directamente el término de fuerza. Llamando $s\in\{0,1\}$ al estado del interruptor de la corrección de media fuerza y $g$ al coeficiente que antecede a la fuerza,

Esto debe igualar a $\delta t F$. La condición es una sola: $g = 1 - s/(2\tau)$. Los dos interruptores no son independientes.

Al desajustarlos, la amplitud se desvía en una proporción predecible. Estos son los cocientes de amplitud medidos dejando fijo el término de fuerza de Guo y cambiando solo los dos interruptores.

Medición y predicción coinciden hasta la tercera cifra decimal. Con $\tau=0.6$, activar solo la corrección de media fuerza agranda la fuerza un 83%. Activar solo el coeficiente la reduce a una sexta parte.

No es un problema sutil de precisión. Es exactamente lo que pasa al copiar tal cual el $F_i$ de Guo de un artículo y dejar como velocidad el $\mathbf{j}/\rho$ del código viejo. Si además, al ver una viscosidad rara, se empieza a tocar $\tau$, el pozo se hace más hondo. Cada vez que cambia $\tau$, la razón del error se mueve con él.

Hay una forma de reconocerlo por los síntomas. Si al duplicar la malla la razón del error no baja, no se trata de error de discretización. Si la razón converge hacia 1.5 conforme $\tau$ se acerca a 1 y crece como si divergiera conforme $\tau$ baja hacia 0.5, esa es la huella de $1+\frac{1}{2\tau}$. Y si al aumentar la fuerza la razón se mantiene igual, el problema es de contabilidad y no de no linealidad. Cuando los tres síntomas coinciden, en vez de sospechar del esquema conviene abrir primero la definición de la velocidad.

Python — se desajustan los dos interruptores#

Este es un solver de canal D2Q9 con los dos interruptores expuestos como argumentos. Como el flujo es homogéneo en $x$, se conserva una sola columna. Los cinco esquemas de la tabla anterior se obtienen cambiando únicamente source_terms y la velocidad que entra en el equilibrio, así que aquí queda fijo el término de fuerza de Guo.

import numpy as np
 
EX = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
EY = 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])
OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])
CS2 = 1.0 / 3.0
 
 
def lattice_equilibrium(rho, ux, uy):
    feq = np.empty((9, rho.size))
    usq = ux**2 + uy**2
    for i in range(9):
        eu = EX[i] * ux + EY[i] * uy
        feq[i] = W[i] * rho * (1 + eu/CS2 + eu**2/(2*CS2**2) - usq/(2*CS2))
    return feq
 
 
def source_terms(rho, ux, uy, fx):
    """Término de fuerza de Guo, antes de la ganancia (1 - 1/2tau)."""
    src = np.empty((9, rho.size))
    for i in range(9):
        eu = EX[i]*ux + EY[i]*uy
        src[i] = W[i] * ((EX[i] - ux)/CS2 + eu*EX[i]/CS2**2) * fx
    return src
 
 
def run_forced_channel(tau, half_shift, prefactor, ny=33, fx=1.0e-5, steps=40000):
    rho = np.ones(ny)
    f = lattice_equilibrium(rho, np.zeros(ny), np.zeros(ny))
    gain = (1.0 - 1.0/(2*tau)) if prefactor else 1.0
    for _ in range(steps):
        rho = f.sum(axis=0)
        jx = (f*EX[:, None]).sum(axis=0)
        jy = (f*EY[:, None]).sum(axis=0)
        ux = (jx + (0.5*fx if half_shift else 0.0)) / rho   # interruptor 1
        uy = jy / rho
        fpost = f - (f - lattice_equilibrium(rho, ux, uy))/tau \
                + gain*source_terms(rho, ux, uy, fx)        # interruptor 2
        for i in range(9):                                   # propagación en y
            f[i] = fpost[i] if EY[i] == 0 else np.roll(fpost[i], EY[i])
        for i in range(9):                                   # half-way bounce-back
            if EY[i] > 0:
                f[i, 0] = fpost[OPP[i], 0]
            elif EY[i] < 0:
                f[i, -1] = fpost[OPP[i], -1]
    jx = (f*EX[:, None]).sum(axis=0)
    return (jx + 0.5*fx) / f.sum(axis=0)      # velocidad física, siempre desplazada
 
 
NY, FX = 33, 1.0e-5
y = np.arange(NY) + 0.5
for tau in (0.6, 1.0, 1.8):
    exact = FX / (2*CS2*(tau - 0.5)) * y * (NY - y)
    ok = run_forced_channel(tau, True, True).max() / exact.max()
    m1 = run_forced_channel(tau, True, False).max() / exact.max()
    m2 = run_forced_channel(tau, False, True).max() / exact.max()
    print(f"tau={tau:4.2f}  both={ok:.4f}  shift_only={m1:.4f} (pred {1+1/(2*tau):.4f})"
          f"  gain_only={m2:.4f} (pred {1-1/(2*tau):.4f})")

La salida es la siguiente.

tau=0.60  both=0.9991  shift_only=1.8316 (pred 1.8333)  gain_only=0.1666 (pred 0.1667)
tau=1.00  both=1.0003  shift_only=1.5002 (pred 1.5000)  gain_only=0.5005 (pred 0.5000)
tau=1.80  both=1.0074  shift_only=1.2867 (pred 1.2778)  gain_only=0.7280 (pred 0.7222)

Con $\tau=1.8$ la desviación respecto de la predicción es del 0.7%, y en las mismas condiciones el error de discretización propio del esquema es del 0.74%. Son del mismo tamaño. Es decir, la desviación viene de la malla y no de la contabilidad.

Lo que conviene recordar#

• El $(1-\frac{1}{2\tau})$ delante del término de fuerza y el $+\delta t\mathbf{F}/2\rho$ de la velocidad son un par que salió junto del cambio de variable de la integración trapezoidal. Usar solo uno de los dos reescala la fuerza en un factor $1\pm\frac{1}{2\tau}$. Con $\tau=0.6$ eso es un 83% de exceso.
• El Poiseuille estacionario y unidireccional no distingue los esquemas de forcing. Los cinco dan el mismo error hasta la cuarta cifra. La comparación entre esquemas hay que hacerla en flujo no estacionario o con una fuerza no homogénea.
• En cualquiera de los esquemas, el momento de primer orden del propio término de fuerza está ajustado a $\mathbf{F}$. La diferencia está en el momento de segundo orden, es decir, en el tensor de esfuerzos, y el valor objetivo lleva dentro el $(1-\frac{1}{2\tau})$. Al revisar un esquema nuevo, empezar por M2 es lo más rápido.