¿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. La derivación ubica la trampa. Los números miden el tamaño.

El dial de momentos cambia el esquema, $\tau$, $|u|$ y el prefactor, y las barras M0, M1 y M2 se mueven con esos controles.

Por más que gire el dial, M0 se queda en cero. Al subir $|u|$, solo el M2 de plain pierde el objetivo. Al apagar el prefactor, M1 y $\Delta j$ se rompen en el mismo paso.

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. El medio paso que la regla del trapecio deja en la viscosidad sale de esa misma raíz. 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.

El objetivo de M1 no es $\mathbf{F}$ sino $(1-\frac{1}{2\tau})\mathbf{F}$, porque la barra mide el término ya multiplicado por el coeficiente. El $\mathbf{F}/(2\tau)$ que falta vuelve con el equilibrio desplazado medio paso, y $\Delta j$ cae sobre $\mathbf{F}$. Solo M2 separa los esquemas. Plain tiene la forma equivocada, así que ningún coeficiente lo arregla. Al empujar $\tau$ hacia 2, la desviación de Shan–Chen crece.

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$.

Cambiar de esquema deja quieta la amplitud de la parábola. Error relativo máximo tras 40 000 pasos, 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.

Con $\tau=0.6$, $2\tau=1.2$ y $1/(2\tau)=5/6$, o sea 0.8333. Ambos interruptores activos significan $s=1$ y $g=1-1/(2\tau)=1/6$. Al dividir el momento inyectado en un paso por $\delta t F$ queda $(1/\tau)\cdot(s/2)+g$. Sustituyendo, $(1/0.6)\cdot 0.5 + 1/6 = 5/6 + 1/6 = 1$. Coincide con el objetivo. La columna both, 0.9991, es ese 1 menos solo lo que se lleva la malla.

Con solo la corrección, $g$ queda fijo en 1. El mismo lugar pasa a $5/6+1=1.8333$. La tabla da 1.8316. La diferencia 0.0017 es cerca del 0.09% de la predicción. Con solo el coeficiente, $s=0$ y $g=1/6$, así que la inyección es 0.1667. El valor medido 0.1666 se encuentra en la cuarta cifra.

Con $\tau=1$, $1/(2\tau)=1/2$. Solo la corrección predice 1.5 y solo el coeficiente 0.5. La tabla dice 1.5002 y 0.5005. Ambos activos dan 1.0003, un 0.03% por encima de 1. Con $\tau=1.8$, $1/(2\tau)=1/3.6=0.2778$. Las predicciones 1.2778 y 0.7222 quedan junto a 1.2867 y 0.7280. $(1.2867-1.2778)/1.2778$ es cerca de 0.007, un 0.7%. El error relativo del esquema a ese $\tau$ es 0.7358%. Mismo tamaño.

La tabla de error relativo se lee igual. 0.0875% con $\tau=0.6$, 0.0306% con $\tau=1.0$ y 0.7358% con $\tau=1.8$. Con $\nu=c_s^2(\tau-1/2)$ y $c_s^2=1/3$, las viscosidades son $1/30$, $1/6$ y $13/30$. Los cocientes de $(\tau-1/2)$ van de 0.1 a 0.5 a 1.3, o sea 1 a 5 a 13. El pico es $u(H/2)=F_x H^2/(8\nu)$, de modo que con $F_x=10^{-5}$ y $H=33$ la altura va al revés de $\nu$. El pico con $\tau=1.8$ es $1/13$ del pico con $\tau=0.6$. El porcentaje crece ahí. Las cinco casillas siguen compartiendo un solo número.

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)

En el laboratorio del canal, alinear los botones de esquema, $\tau$, el medio paso y el prefactor con las columnas both, shift_only y gain_only de la salida de arriba.

Los cuatro esquemas dejan la curva turquesa sobre la parábola punteada. Al apagar un solo interruptor, con $\tau=0.6$ la amplitud salta hacia 1.83 o hacia 0.17.

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.

Con solo el medio paso, la fuerza queda en 1.83#

Con $\tau=0.6$, shift_only marca 1.8316. La predicción $1+\frac{1}{2\tau}=1.8333$ coincide hasta la tercera cifra. Copiar el término de Guo y dejar la velocidad en $\mathbf{j}/\rho$ cae justo en esa casilla: corrección activa, prefactor apagado. Duplicar la malla no baja el cociente. Llevar $\tau$ a 1 deja el cociente cerca de 1.5. Bajar $\tau$ hacia 0.5 lo agranda. Es la huella de $1+\frac{1}{2\tau}$.

El momento de primer orden de los cinco esquemas es $\mathbf{F}$. La separación está en el segundo, y el objetivo lleva $(1-\frac{1}{2\tau})$. Un canal estacionario no tiene por dónde meter ese error $xx$ en el balance de cortante, así que las cinco columnas comparten 0.0875%, 0.0306% y 0.7358%. M2 importa cuando el flujo deja de ser estacionario o la fuerza salta en la interfaz. Llenar también el momento de tercer orden es otra cuenta, en el límite de tercer momento del término de fuerza.