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.

Al subir el control de velocidad desde cero, la curva ámbar de la maxwelliana se separa de la recta azul de la red.

sweep u muestra que la separación sigue una potencia impar de $u$. Con 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$.

Conviene dejar la aritmética escrita. Con $c_s^2 = 1/3$ y $U = 0.2$, $U^2 = 0.04$ y $U^2/c_s^2 = 0.12$, así que la viscosidad en $x$ sola cae un 12%. Taylor–Green tiene $k_x = k_y$, y la mitad de ese error es 6%. La fila $U_0 = 0.20$ mide $\nu = 0.0940296$ frente a $0.100000$, o sea $-5.9704\%$. Entre 6.00 y 5.97 hay 0.03 puntos porcentuales.

Mover $\tau$ no cambia el sitio. En $\tau = 0.6$ la teoría da $c_s^2 \times 0.1 = 0.0333333$. En $\tau = 1.0$ da $c_s^2 \times 0.5 = 0.1666667$, cinco veces más. El error relativo en $U_0 = 0.15$ solo pasa de $-3.4076\%$ a $-3.3099\%$. Al meter $U^2 = 0.0225$ en $-U^2/(2c_s^2)$ sale $-0.03375$, es decir $-3.375\%$. Las tres filas se pegan a ese número. El álgebra que produce $\tau - 1/2$ está en la expansión de segundo orden de la LBE discreta. Aquí el factor se cancela arriba y abajo.

La columna de auditoría ya guarda el defecto. En $u = 0.20$ la suma de la red es $0.20000000$, el lado maxwelliano es $0.20800000$ y la brecha es $0.00800000$. Como $0.2^3 = 0.008$, el cociente vale 1.00000. En $u = 0.05$ la brecha es $0.00012500$ y el cociente sigue en 1. En $u = 0.30$ la brecha $0.02700000$ iguala a $0.3^3$. Poblaciones al azar dejan una diferencia de $0.000\mathrm{e}{+}00$. Bajar el marco y no poder guardar $\rho u^3$ son el mismo término.

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%

Al subir frame velocity U, el vórtice de la derecha se desplaza y su decaimiento se vuelve más lento.

El mismo fluido se desvanece más despacio a la derecha. Al mover el control de $\tau$, el error de abajo casi no cambia. Ahí están el $-3.41\%$ y el $-3.31\%$ de la tabla.

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.

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.

Recalcular $\tau$ en cada nivel deja este defecto intacto, porque el error relativo no depende de $\tau$. Una cantidad que nadie ingresó como propiedad queda fijada 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. El tratamiento de contorno no es el primer sospechoso.

Bajar $u/c_s$ o ampliar las velocidades hasta $\pm 2$#

D2Q9 pierde $\rho u^3$ en tercer orden. En $U_0 = 0.2$ la factura es $-5.97\%$. Refinar la malla, o subir $\tau$ de 0.6 a 1.0, apenas mueve esa proporción.

Hay tres cortes que sí la reducen. Al bajar $u/c_s$ el error cae con $U^2$. También se puede devolver $\partial_\gamma(\rho u_\alpha u_\beta u_\gamma)$ a la colisión. O ampliar el conjunto de velocidades. Por qué tres velocidades quedan trabadas en $c^3 = c$ sigue en las raíces de Hermite que eligen las velocidades de la red. Sin cambiar esa elección, un marco en 0.2 vuelve a recortar la viscosidad un 5.97%.