Empujar el marco a 0.2 redujo la viscosidad un 5.97% — el tercer momento que falta en el equilibrio de LBM

¿Deben coincidir la viscosidad del agua quieta y la del agua en movimiento?#

Supongamos que se mide el mismo fluido dos veces. Una vez dentro de una caja en reposo y otra vez mientras esa caja avanza a velocidad constante. El coeficiente de viscosidad debería salir idéntico en ambos casos. La invariancia galileana, es decir, que las leyes físicas se escriben igual en dos marcos que se mueven a velocidad relativa constante, lo exige.

Al repetir ese experimento en un solver de Boltzmann en red (LBM) con D2Q9, el segundo valor sale un 5.97% más pequeño. Refinar la malla no lo corrige. Cambiar el tiempo de relajación $\tau$ apenas mueve esa proporción. Este artículo persigue ese 5.97% hasta una sola línea de momentos de la distribución de equilibrio. La conclusión, por adelantado: no es un problema de estabilidad sino de álgebra. El conjunto de velocidades de la red no puede construir el término que falta.

Al desarrollar en Taylor la única línea que resuelve la red#

LBM avanza una sola ecuación.

Aquí $f_i$ es la población que viaja sobre la velocidad de red $\mathbf{c}_i$, $f_i^{eq}$ es el equilibrio local y $\tau$ es el tiempo de relajación.

Al desarrollar el lado izquierdo en $\Delta t$ sobrevive un término de segundo orden.

Ese término de segundo orden es el que genera el célebre $\tau - 1/2$ en el desarrollo de Chapman–Enskog, la técnica multiescala que descompone la distribución en una serie de potencias del número de Knudsen. El paso siguiente del mismo desarrollo impone una condición: para recuperar Navier–Stokes, el equilibrio debe reproducir exactamente cuatro momentos.

Los tres primeros fijan masa, cantidad de movimiento y presión. El cuarto, el momento de tercer orden, es el que reemplaza la derivada temporal del esfuerzo viscoso. Si ese falla, la continuidad queda intacta y el nivel de Euler también: solo se contamina el término viscoso.

Conviene comprobar ese tercer momento directamente en la simulación siguiente.

Al desplazar el control de velocidad, la curva ámbar (lo que exige una maxwelliana) se separa de la recta azul (lo que D2Q9 entrega de verdad). Con sweep u el barrido automático muestra que la separación sigue una potencia impar de $u$. Al activar speeds ±2 aparecen cinco barras y la brecha se cierra a cero.

$c_x^3 = c_x$, de modo que el término no puede existir#

Las componentes en $x$ de D2Q9 toman solo tres valores: $-1$, $0$ y $1$. Cada uno de ellos coincide con su propio cubo. Por lo tanto, para cualquier conjunto de poblaciones $f_i$,

se cumple de forma idéntica. No importa cómo se diseñe el equilibrio. Incluso al rellenar las nueve casillas con números aleatorios la identidad sigue en pie. El tercer momento ya está clavado al primero.

Lo que pide una maxwelliana es $\rho(u_x^3 + 3 c_s^2 u_x)$ y, como en D2Q9 $c_s^2 = 1/3$, resulta $3c_s^2 = 1$ y el objetivo se reduce a $\rho(u_x + u_x^3)$. La diferencia entre ambos es exactamente $\rho u_x^3$. No un error de ese orden de magnitud: ese término, sin resto alguno.

El rostro que el término ausente deja en la ecuación de momento#

El momento que falta cae directamente sobre el esfuerzo viscoso. Al llevar el desarrollo de Chapman–Enskog hasta el final aparece un término adicional.

Considérese ahora un flujo medio $U$ en dirección $x$ con una perturbación pequeña $u'$ encima. Al conservar solo la parte lineal en $u'$, el término queda como $-(\tau - 1/2)\Delta t\, U^2 \partial_x^2 (\rho u'_\alpha)$, que tiene exactamente la forma de un término viscoso. Es decir, el coeficiente de viscosidad en la dirección $x$ se reescribe por completo.

El error relativo es $-U^2 / c_s^2$. Nótese que $\tau$ se cancela entre numerador y denominador: por más que se ajuste la relajación, el porcentaje no cambia. Y el error crece con el cuadrado de la velocidad, que en número de Mach es simplemente $-\mathrm{Ma}^2$.

Medir el mismo vórtice desde dos marcos con Python#

Se coloca un vórtice de Taylor–Green sobre una red periódica de 64×64 y se le suma una velocidad uniforme $U_0$. Ajustando la pendiente logarítmica de la amplitud que decae se recupera el coeficiente de viscosidad. Solo cambia $U_0$; todo lo demás queda fijo.

import numpy as np
 
CX  = np.array([0, 1, 0, -1,  0, 1, -1, -1,  1], dtype=float)
CY  = np.array([0, 0, 1,  0, -1, 1,  1, -1, -1], dtype=float)
W   = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CS2 = 1.0 / 3.0
 
def f_equilibrium(rho, ux, uy):
    # equilibrio estandar, desarrollado hasta segundo orden
    cu = CX[:, None, None] * ux + CY[:, None, None] * uy
    usq = ux * ux + uy * uy
    return W[:, None, None] * rho * (1 + cu / CS2 + cu * cu / (2 * CS2**2) - usq / (2 * CS2))
 
print("[1] third moment audit:  sum f_i c_ix^3   vs   Maxwell rho*(u^3 + 3*cs2*u)")
print("   u        lattice      Maxwell        gap        gap/u^3")
for u in (0.05, 0.10, 0.20, 0.30):
    feq = f_equilibrium(np.ones((1, 1)), np.full((1, 1), u), np.zeros((1, 1)))
    lat = float((feq * (CX**3)[:, None, None]).sum())
    exact = u**3 + 3 * CS2 * u
    print(f"  {u:4.2f}  {lat:11.8f}  {exact:11.8f}  {exact-lat:11.8f}   {(exact-lat)/u**3:8.5f}")
 
rng = np.random.default_rng(7)
frand = rng.random((9, 1, 1))
d = float((frand * (CX**3)[:, None, None]).sum() - (frand * CX[:, None, None]).sum())
print(f"  any random f:  sum f c_x^3 - sum f c_x = {d:.3e}   (c_x^3 = c_x, so always 0)")
 
def taylor_green_run(U0, tau=0.8, N=64, steps=900, amp=0.04, warmup=300, sample=25):
    # el mismo vortice montado sobre una velocidad uniforme U0; la viscosidad se recupera del decaimiento
    k = 2 * np.pi / N
    x = np.arange(N)[:, None] * np.ones(N)[None, :]
    y = np.ones(N)[:, None] * np.arange(N)[None, :]
    ux = U0 - amp * np.cos(k * x) * np.sin(k * y)
    uy = amp * np.sin(k * x) * np.cos(k * y)
    rho = np.ones((N, N))
    f = f_equilibrium(rho, ux, uy)
    ts, logs = [], []
    for n in range(steps + 1):
        rho = f.sum(axis=0)
        ux = (f * CX[:, None, None]).sum(axis=0) / rho
        uy = (f * CY[:, None, None]).sum(axis=0) / rho
        if n % sample == 0:
            up, vp = ux - ux.mean(), uy - uy.mean()   # fluctuacion que queda tras quitar el flujo medio
            ts.append(n)
            logs.append(0.5 * np.log(2 * np.mean(up * up + vp * vp)))
        f += (f_equilibrium(rho, ux, uy) - f) / tau   # colision BGK
        for i in range(9):                            # propagacion
            f[i] = np.roll(np.roll(f[i], int(CX[i]), axis=0), int(CY[i]), axis=1)
    ts, logs = np.array(ts, dtype=float), np.array(logs)
    m = ts >= warmup
    slope = np.polyfit(ts[m], logs[m], 1)[0]
    return -slope / (2 * k * k)
 
print("\n[2] same vortex measured in a uniformly moving frame  (tau=0.8, nu_theory=0.100000)")
nu_th = CS2 * (0.8 - 0.5)
print("   U0      nu_eff      rel.err     rel.err/U0^2")
for U0 in (0.00, 0.05, 0.10, 0.15, 0.20):
    nu = taylor_green_run(U0)
    e = nu / nu_th - 1
    tail = f"{e/U0**2:10.4f}" if U0 > 0 else "         -"
    print(f"  {U0:4.2f}  {nu:10.7f}  {e:+9.4%}  {tail}")
 
print("\n[3] does the error depend on tau?  (U0=0.15 fixed)")
print("   tau     nu_theory    nu_eff      rel.err")
for tau in (0.6, 0.8, 1.0):
    nu = taylor_green_run(0.15, tau=tau)
    th = CS2 * (tau - 0.5)
    print(f"  {tau:4.2f}  {th:10.7f}  {nu:10.7f}  {nu/th-1:+9.4%}")

[1] third moment audit:  sum f_i c_ix^3   vs   Maxwell rho*(u^3 + 3*cs2*u)
   u        lattice      Maxwell        gap        gap/u^3
  0.05   0.05000000   0.05012500   0.00012500    1.00000
  0.10   0.10000000   0.10100000   0.00100000    1.00000
  0.20   0.20000000   0.20800000   0.00800000    1.00000
  0.30   0.30000000   0.32700000   0.02700000    1.00000
  any random f:  sum f c_x^3 - sum f c_x = 0.000e+00   (c_x^3 = c_x, so always 0)
 
[2] same vortex measured in a uniformly moving frame  (tau=0.8, nu_theory=0.100000)
   U0      nu_eff      rel.err     rel.err/U0^2
  0.00   0.1000175   +0.0175%           -
  0.05   0.0996433   -0.3567%     -1.4268
  0.10   0.0985207   -1.4793%     -1.4793
  0.15   0.0966495   -3.3505%     -1.4891
  0.20   0.0940296   -5.9704%     -1.4926
 
[3] does the error depend on tau?  (U0=0.15 fixed)
   tau     nu_theory    nu_eff      rel.err
  0.60   0.0333333   0.0321975   -3.4076%
  0.80   0.1000000   0.0966495   -3.3505%
  1.00   0.1666667   0.1611502   -3.3099%

El error viaja sobre $U^2$, no sobre $\tau$#

La última columna de [1] es 1.00000 en todas las filas. La afirmación de que el defecto es exactamente $\rho u^3$ se sostiene hasta el quinto decimal.

En [2] el marco en reposo reporta $+0.0175\%$, que marca el ruido de fondo del propio método de medición. Al subir $U_0$, el error camina hacia $-0.36\%$, $-1.48\%$, $-3.35\%$ y $-5.97\%$. La última columna divide ese error entre $U_0^2$ y, a medida que $U_0$ disminuye, se asienta en $-1.5$.

Ese $-1.5$ es la predicción de la sección anterior. El error relativo $-U^2/c_s^2$ afecta solo a la componente $x$, y el modo de Taylor–Green tiene $k_x = k_y$, de modo que su tasa de decaimiento promedia ambas direcciones. La mitad de $-U^2/c_s^2$ es $-U^2/(2c_s^2) = -1.5\,U^2$. Para $U_0 = 0.2$ la predicción da $-6.00\%$ frente a un valor medido de $-5.97\%$.

El resultado [3] pesa todavía más. Subir $\tau$ de 0.6 a 1.0 multiplica la viscosidad por cinco y el error relativo solo se desplaza de $-3.41\%$ a $-3.31\%$. Este error no se puede enterrar bajo más viscosidad.

Al subir frame velocity U, el vórtice de la derecha se desplaza por el panel mientras se desvanece más despacio que el de la izquierda. El mismo fluido y, sin embargo, el de la derecha se resiste a morir. Al mover el control de $\tau$ se observa además que el valor de error de abajo casi no se inmuta.

Tres maneras de corregirlo y lo que cuesta cada una#

Bajar la velocidad. El error va como $U^2$, así que reducir a la mitad la velocidad de red lo deja en una cuarta parte. La factura llega en número de pasos: reproducir el mismo tiempo físico exige más iteraciones. La regla empírica de LBM de "mantener el número de Mach por debajo de 0.1" suele venderse como criterio de estabilidad, pero lo que en realidad se activa primero es este término ausente.

Añadir un término de corrección. Se calcula $\partial_\gamma(\rho u_\alpha u_\beta u_\gamma)$ por diferencias finitas y se inyecta en la colisión con signo opuesto. El costo es una evaluación de gradiente más la pérdida parcial de la colisión estrictamente local que hace atractivo a LBM. Aun así, en marcos rotantes y otros problemas con flujo medio grande sigue siendo la opción más barata.

Ampliar el conjunto de velocidades. Una red multivelocidad que incluya $\pm 2$ cumple $c^3 \neq c$, con lo cual el tercer momento puede igualarse de forma exacta. El botón speeds ±2 de la primera visualización resuelve precisamente eso. Se paga con un estarcido más ancho, más memoria y un tratamiento de contorno más enredado, y además algunos pesos pueden volverse negativos, lo que vuelve a poner la estabilidad sobre la mesa.

Dónde asoma este defecto en cálculos reales#

En todo lugar donde el flujo medio sea grande. Marcos rotantes de turbomaquinaria, mallas deslizantes, turbulencia montada sobre una corriente uniforme rápida, un marco pegado a un cuerpo en movimiento. Un solver impecable en los casos de referencia estáticos que después falla el número de Reynolds efectivo justo en estos problemas está mostrando este síntoma.

En el reescalado del no equilibrio en mallas LBM refinadas la solución pasaba por recalcular $\tau$ en cada nivel, y este defecto sobrevive intacto a esa operación, porque su tamaño relativo no depende de $\tau$. Comparte además la estructura de las dos constantes escondidas en el LBM térmico: una cantidad que nadie ingresó como propiedad del material queda decidida por la estructura de la red.

El diagnóstico es barato. Basta repetir el mismo caso con un flujo medio añadido y comparar la tasa de decaimiento o el coeficiente de arrastre. Si la diferencia crece como $U^2$, el término es este. No hace falta sospechar del tratamiento de contorno, como en la prueba de paridad en la voxelización de STL.

La próxima vez que la viscosidad dependa del marco#

En Boltzmann en red la estabilidad no es lo único que limita la velocidad. El orden de momentos que el equilibrio alcanza la limita antes. D2Q9 pierde $\rho u^3$ en tercer orden y la factura llega como un error de viscosidad de $-\mathrm{Ma}^2$.

Cuando los números no cuadran, los reflejos habituales son refinar la malla y retocar $\tau$. Ninguno de los dos toca este error. Solo se reduce al bajar $u/c_s$, al devolver a mano el término ausente o al ampliar el conjunto de velocidades.