Completar el tercer momento con honestidad multiplicó el error por 2.8 — lo que un término de fuerza realmente debe

El tercer momento de un término de fuerza llega a la ecuación de momento solo en orden Burnett, y una corrección tensorial que la red concede a medias resulta peor que ninguna.

¿Hasta qué momento debe ajustarse un término de fuerza?#

Al meter gravedad o tensión superficial en un código LBM se ajustan tres líneas. Momento cero nulo, primer momento $\mathbf{F}$ , segundo momento $\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$ . ¿Y el tercero? El término de fuerza de la ecuación de Boltzmann continua también tiene un tercer momento bien definido. Ajustar ese valor sobre la red parece una ganancia de precisión gratuita.

Los números salieron al revés. Al completar el tercer momento con lo que exige la teoría cinética, el error pasó de 0.053% a 0.150%. Los deslizadores $u_x$ y $u_y$ son la velocidad. Al subirlos, las filas de momento que coincidían se abren. La casilla cierra solo las filas amarillas.

ux = 0.22 uy = 0.18 F angle = 0 deg  fill the third moment

Drag ux, uy toward zero and every row lines up. Push them out and the two yellow rows open a gap that the checkbox can close, while the red row cannot be touched at all: on D2Q9 cx^3 is the same vector as cx. Largest remaining gap, relative to the biggest moment on screen: 0.00%.

Al devolver la velocidad hacia cero, las filas vuelven a coincidir. La fila roja sigue abierta aunque los deslizadores recorran todo el rango.

El tercer momento que exige la ecuación de Boltzmann#

Con una fuerza volumétrica $\mathbf{F}$ , el término de fuerza de la ecuación de Boltzmann es una derivada en el espacio de velocidades.

\Omega^F = -\frac{F_\lambda}{\rho}\,\frac{\partial f^{eq}}{\partial \xi_\lambda}

Aquí $\xi$ es la velocidad molecular y $\rho$ la densidad. Una integración por partes cierra el tercer momento.

\int \xi_\alpha \xi_\beta \xi_\gamma\, \Omega^F \, d\boldsymbol{\xi} = \frac{F_\alpha}{\rho} M_{\beta\gamma} + \frac{F_\beta}{\rho} M_{\alpha\gamma} + \frac{F_\gamma}{\rho} M_{\alpha\beta}

con $M_{\alpha\beta}=\int f^{eq}\xi_\alpha\xi_\beta\,d\boldsymbol{\xi} = \rho(u_\alpha u_\beta + c_s^2\delta_{\alpha\beta})$ . Así, el tercer momento que se pide se escribe de esta forma.

\Phi_{\alpha\beta\gamma} = F_\alpha\!\left(u_\beta u_\gamma + c_s^2\delta_{\beta\gamma}\right) + F_\beta\!\left(u_\alpha u_\gamma + c_s^2\delta_{\alpha\gamma}\right) + F_\gamma\!\left(u_\alpha u_\beta + c_s^2\delta_{\alpha\beta}\right)

La pieza $c_s^2$ y la pieza $uuF$ conviven ahí. La primera no depende de la velocidad; la segunda es cuadrática en ella.

Lo que Guo forcing suministra de verdad#

El término de fuerza sobre la red suele escribirse en la forma de Guo. Quitando el factor de medio paso, al corchete lo llamamos $\hat{S}_i$ .

\hat{S}_i = w_i\left[\frac{\mathbf{c}_i\cdot\mathbf{F}}{c_s^2} + \frac{\left(u_\alpha F_\beta + u_\beta F_\alpha\right)\left(c_{i\alpha}c_{i\beta} - c_s^2\delta_{\alpha\beta}\right)}{2\,c_s^4}\right]

$w_i$ son los pesos de la red y $\mathbf{c}_i$ las velocidades discretas. El $S_i$ real lleva además un factor $(1-1/2\tau)$ . Por qué hace falta ese factor se trató aparte en el medio paso que deja la LBE discreta.

Al calcular el tercer momento de este $\hat{S}_i$ con la isotropía de cuarto orden de D2Q9, solo queda lo siguiente.

\sum_i \hat{S}_i c_{i\alpha}c_{i\beta}c_{i\gamma} = c_s^2\left(F_\alpha\delta_{\beta\gamma} + F_\beta\delta_{\alpha\gamma} + F_\gamma\delta_{\alpha\beta}\right)

La pieza $c_s^2$ es exacta. De los términos $uuF$ no hay ninguno. El segundo corchete arma un quinto momento, que en D2Q9 se anula.

Eso es lo que abrieron los deslizadores de velocidad: dos filas amarillas se mueven y la fila roja sigue bloqueada.

$c_x^3 = c_x$ : los dos componentes que la red bloqueó#

La fila roja es la red misma. Las componentes de velocidad de D2Q9 solo valen $-1, 0, 1$ . Entonces $c_{ix}^3 = c_{ix}$ en toda dirección. Los dos vectores son el mismo vector.

\sum_i \hat{S}_i c_{ix}^3 = \sum_i \hat{S}_i c_{ix} = F_x

$\Phi_{xxx}$ no es un valor elegible. Queda fijado en el mismo instante en que el primer momento se fija a $F_x$ . La exigencia es $3F_x(u_x^2+c_s^2) = F_x + 3F_x u_x^2$ , de modo que $3F_x u_x^2$ es un déficit permanente.

En cambio $\Phi_{xxy}$ y $\Phi_{xyy}$ son independientes. De nueve grados de libertad, los momentos de orden 0 a 2 toman seis y quedan tres. Así que de los cuatro componentes, dos son libres y dos están bloqueados. Solo se puede arreglar la mitad de un tensor.

Un término que no transporta masa ni momento acaba moviendo momento#

Se puede construir un vector de corrección $v_i$ que cambie solo $\Phi_{xyy}$ , en una cantidad $d$ . Sus momentos de orden 0, 1 y 2 se anulan, así que en el sitio no toca densidad, ni momento, ni esfuerzo. En D2Q9 ese vector tiene esta forma.

v_i = \left(0,\,-\tfrac12,\,0,\,\tfrac12,\,0,\,\tfrac14,\,-\tfrac14,\,-\tfrac14,\,\tfrac14\right)

Un paso de propagación cambia la historia. Cada componente vuela hacia un vecino distinto. Si se recoge solo el momento en $x$ , queda lo siguiente.

\Delta(\rho u_x)(y) = \tfrac{1}{2}\left[\,d(y-1) - 2\,d(y) + d(y+1)\,\right]

Una diferencia central de segundo orden. Es exactamente nula cuando $d$ es uniforme y crece con el cuadrado del número de onda cuando $d$ oscila. Con una forma uniforme el cambio de celda a celda desaparece. Al acortar una senoide, ese cambio crece con el cuadrado del número de onda. La suma de las celdas sigue en cero.

Un canal con el patrón de fuerza fijo y un flujo transversal#

La verificación necesita un problema con solución exacta. Se fija una fuerza volumétrica senoidal en $x$ dentro de una caja periódica y se sopla un flujo transversal uniforme $V$ en $y$ . Paredes porosas que inyectan por un lado y aspiran por el otro.

\nu \frac{d^2 u_x}{dy^2} - V\frac{d u_x}{dy} + \frac{a_0}{\rho}\sin(ky) = 0

$\nu$ es la viscosidad cinemática, $k$ el número de onda y $a_0$ la amplitud de la fuerza. Un solo modo da forma cerrada.

u_x(y) = C_1\sin(ky) - \frac{V C_1}{\nu k}\cos(ky), \qquad C_1 = \frac{a_0/\rho}{\nu k^2 + V^2/\nu}

El flujo transversal hace dos cosas. Recorta la amplitud y retrasa la fase en $\arctan\!\left(V/\nu k\right)$ . A la vez $u_y = V$ sigue vivo, de modo que $d = F_x u_y^2$ oscila en el espacio. Ese es el interruptor de la fuga anterior.

Cuatro combinaciones sobre la misma red#

Red periódica de 16 nodos, $\tau = 1$ , $a_0 = 4\times10^{-3}$ , cuatro corridas. Guo estándar, tercer momento de la fuerza completado, tercer momento del equilibrio completado, y ambos.

import numpy as np
 
CX  = np.array([0, 1, 0,-1, 0, 1,-1,-1, 1], float)
CY  = np.array([0, 0, 1, 0,-1, 1, 1,-1,-1], float)
WT  = 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 hermite3_vectors():
    """Momentos 0 a 2 nulos, tercer momento unitario solo en xxy / xyy."""
    rows = np.vstack([np.ones(9), CX, CY, CX*CX, CX*CY, CY*CY, CX*CX*CY, CX*CY*CY])
    pinv = np.linalg.pinv(rows)
    return pinv @ np.eye(8)[6], pinv @ np.eye(8)[7]
 
VXXY, VXYY = hermite3_vectors()
 
def feq_second_order(rho, ux, uy):
    cu = CX[:, None, None]*ux + CY[:, None, None]*uy
    return WT[:, None, None]*rho*(1 + cu/CS2 + cu*cu/(2*CS2**2)
                                  - (ux*ux + uy*uy)/(2*CS2))
 
def force_populations(ux, uy, fx, fy, tau):
    """Guo forcing: primer momento F, segundo momento uF+Fu."""
    cF  = CX[:, None, None]*fx + CY[:, None, None]*fy
    gxx, gyy, gxy = ux*fx, uy*fy, 0.5*(ux*fy + uy*fx)
    quad = (gxx*(CX[:, None, None]**2 - CS2) + gyy*(CY[:, None, None]**2 - CS2)
            + 2*gxy*(CX*CY)[:, None, None])
    return WT[:, None, None]*(1 - 1/(2*tau))*(cF/CS2 + quad/CS2**2)
 
def analytic_wave(ny, nu, V, a0):
    """Solucion estacionaria periodica de nu u'' - V u' + a0 sin(ky) = 0."""
    k = 2*np.pi/ny
    y = np.arange(ny, dtype=float)
    c1 = a0/(nu*k*k + V*V/nu)
    return c1*np.sin(k*y) - (V*c1/(nu*k))*np.cos(k*y)
 
def run_cross_flow(ny, tau, V, a0, cubic_eq, cubic_force, nstep=6000, tol=1e-12):
    nx = 3
    k  = 2*np.pi/ny
    fx = np.tile(a0*np.sin(k*np.arange(ny)), (nx, 1))
    fy = np.zeros((nx, ny))
    rho = np.ones((nx, ny)); ux = np.zeros((nx, ny)); uy = np.full((nx, ny), V)
    f = feq_second_order(rho, ux, uy)
    half, prev = 1 - 1/(2*tau), ux.copy()
    for s in range(nstep):
        rho = f.sum(0)
        ux  = ((CX[:, None, None]*f).sum(0) + 0.5*fx)/rho
        uy  = ((CY[:, None, None]*f).sum(0) + 0.5*fy)/rho
        feq = feq_second_order(rho, ux, uy)
        if cubic_eq:                      # completa el defecto cubico rho*u*u*u del equilibrio
            feq = feq + VXXY[:, None, None]*(rho*ux*ux*uy) \
                      + VXYY[:, None, None]*(rho*ux*uy*uy)
        src = force_populations(ux, uy, fx, fy, tau)
        if cubic_force:                   # completa el defecto cubico F*u*u del termino de fuerza
            src = src + half*(VXXY[:, None, None]*(2*fx*ux*uy + fy*ux*ux)
                              + VXYY[:, None, None]*(fx*uy*uy + 2*fy*ux*uy))
        f = f - (f - feq)/tau + src
        for i in range(9):
            f[i] = np.roll(f[i], (int(CX[i]), int(CY[i])), axis=(0, 1))
        if s % 1000 == 999:
            if np.max(np.abs(ux - prev)) < tol:
                break
            prev = ux.copy()
    return ux[0]
 
print("D2Q9 lock check:  max|cx^3 - cx| =", float(np.abs(CX**3 - CX).max()))
for tag, mom in (("xxy", CX*CX*CY), ("xyy", CX*CY*CY), ("xxx", CX**3), ("xy", CX*CY)):
    print(f"  correction vector v_xxy :  sum v_i c_{tag:3s} = {float(mom @ VXXY): .6f}")
 
NY, TAU = 16, 1.0
NU, A0  = CS2*(TAU - 0.5), 0.004
print(f"\nny={NY}  tau={TAU}  nu={NU:.5f}  k={2*np.pi/NY:.4f}  a0={A0}")
print("\n   V    |u|max   Guo       +force    +feq3     +both")
for V in (0.05, 0.10, 0.20):
    ref = analytic_wave(NY, NU, V, A0)
    errs = []
    for ce, cf in ((0, 0), (0, 1), (1, 0), (1, 1)):
        u = run_cross_flow(NY, TAU, V, A0, bool(ce), bool(cf))
        errs.append(100*np.linalg.norm(u - ref)/np.linalg.norm(ref))
    print(f"{V:5.2f}  {np.abs(ref).max():.5f}  " + "  ".join(f"{e:7.4f}%" for e in errs))

D2Q9 lock check:  max|cx^3 - cx| = 0.0
  correction vector v_xxy :  sum v_i c_xxy =  1.000000
  correction vector v_xxy :  sum v_i c_xyy = -0.000000
  correction vector v_xxy :  sum v_i c_xxx =  0.000000
  correction vector v_xxy :  sum v_i c_xy  = -0.000000
 
ny=16  tau=1.0  nu=0.16667  k=0.3927  a0=0.004
 
   V    |u|max   Guo       +force    +feq3     +both
 0.05  0.12258   0.5993%   0.5915%   0.0181%   0.0143%
 0.10  0.08374   1.6422%   1.6207%   0.0295%   0.0375%
 0.20  0.04826   3.7307%   3.6840%   0.0531%   0.1503%

La separación entre +feq3 y +both en esa tabla crece al acortar la longitud de onda. El interruptor de forma de onda y el deslizador de longitud son las barras de esa fuga.

d(y) = nodes per wave = 8 amplitude = 1.00  stream

Switch d(y) to uniform: the dots still fly, the green bars all collapse to zero. Switch back to sine and shorten the wavelength: the bars grow like the second difference. Peak delivered momentum 0.00e+0, total over the column 0 (exactly), predicted peak 0.5*|d''| = 0.00e+0.

En uniform las barras verdes quedan en cero. Al acortar la longitud en sine, las barras crecen con el cuadrado del número de onda y la suma de las celdas sigue en cero.

La corrección que convirtió 0.053% en 0.150%#

Leyendo la tabla fila por fila aparece el orden. Los números son para $V=0.20$ .

Combinación	Error $L_2$	frente al estándar
Guo estándar	3.7307%	—
Tercer momento de la fuerza	3.6840%	1.3% mejor
Tercer momento del equilibrio	0.0531%	70 veces mejor
Ambos	0.1503%	2.8 veces peor que solo el equilibrio

Completar el tercer momento de la fuerza rinde 1.3% como máximo. Completar el del equilibrio recorta el error setenta veces. Y al poner los dos vuelve a empeorar. Los tres números cuentan la misma historia.

La razón está en el orden del desarrollo de Chapman-Enskog. El segundo momento de la fuerza entra directo en la ecuación de esfuerzo de $O(\varepsilon)$ . El tercer momento de la fuerza entra en la ecuación de tercer momento de $O(\varepsilon)$ , baja de ahí al esfuerzo y del esfuerzo al momento. Llega a la ecuación de momento en $O(\varepsilon^3)$ , es decir, en orden Burnett. Para una LBM isoterma que apunta a Navier-Stokes, ese término no es algo que haya que ajustar.

El tercer momento del equilibrio, en cambio, entra directo en el esfuerzo de $O(\varepsilon)$ . Por eso domina el defecto cúbico del equilibrio.

El camino por el que la corrección fabricó error#

¿Por qué entonces completar el tercer momento de la fuerza resulta dañino? La fórmula de fuga anterior lo responde. Como $d = F_x u_y^2 = a_0 V^2 \sin(ky)$ , la propagación fabrica una fuerza volumétrica espuria de este tamaño.

\tfrac{1}{2}\left(1-\frac{1}{2\tau}\right) V^2 k^2 a_0

Con $V=0.2$ , $k=0.3927$ , $\tau=1$ resulta $0.154\%$ de $a_0$ . La diferencia máxima medida en el código entre las dos soluciones fue $0.156\%$ de $|u|_{\max}$ . Los dos números coinciden. La red entregó con toda fidelidad un término que la ecuación de momento nunca pidió.

La dependencia $V^2 k^2$ también se ve tal cual. Subiendo $V$ por $0.05 \to 0.10 \to 0.20$ , la diferencia entre ambas soluciones fue $0.010\% \to 0.039\% \to 0.156\%$ . Cuatro veces cada vez.

Dónde el tercer momento empieza a importar de verdad#

Nada de esto vuelve inútil la derivación de $\Phi_{\alpha\beta\gamma}$ . Hay lugares que la necesitan.

En LBM térmica o compresible que resuelve una ecuación de energía, el tercer momento es el flujo de calor. Allí el tercer momento de la fuerza entra directo en la ecuación de energía de $O(\varepsilon^2)$ . No hace falta llegar a Burnett. La LBM multifásica que construye interfases con una ecuación de estado no ideal es parecida. La fuerza es grande y espacialmente abrupta. Cuando $k$ se acerca a uno, $V^2k^2$ ya no es pequeño.

Con $V=0.20$ , Guo estándar da 3.7307%. Completar solo el tercer momento de la fuerza da 3.6840%. La diferencia es 0.0467 puntos porcentuales, cerca del 1.3% de 3.7307. Solo el equilibrio da 0.0531%. Al dividir 3.7307 entre 0.0531 el factor ronda 70. Ambos juntos dan 0.1503%. Al dividir 0.1503 entre 0.0531 sale cerca de 2.83, el factor 2.8 de la tabla. Con $V=0.05$ , ambos juntos es 0.0143%, menor que el equilibrio solo en 0.0181%. Con $V=0.10$ el orden se invierte, 0.0375% contra 0.0295%. Con $V=0.20$ la separación es 0.0972 puntos porcentuales.

Al meter $V=0.2$ , $k=0.3927$ , $\tau=1$ en la fuga, $1-1/(2\tau)$ vale $1/2$ . Por el $1/2$ de delante queda $1/4$ . $V^2=0.04$ y $k^2$ es cerca de $0.1542$ . Un cuarto por 0.04 por 0.1542 es 0.001542, el 0.154% de $a_0$ . El código midió 0.156% de $|u|_{\max}$ . Cabe en 0.002 puntos porcentuales. Doblar $V$ de 0.05 a 0.10 y a 0.20 multiplica $V^2$ por cuatro, y las separaciones 0.010%, 0.039%, 0.156% hacen lo mismo. $|u|_{\max}$ baja por 0.12258, 0.08374, 0.04826. El denominador de $C_1$ es $\nu k^2+V^2/\nu$ con $\nu=0.16667$ , así que un $V$ mayor recorta la amplitud.

Dónde va la mitad que falta de la fuerza es lo que separa la comparación de Guo, Shan–Chen y EDM sobre el mismo factor $(1-1/2\tau)$ . Antes de sumar la línea, toca contar la ecuación objetivo y el orden en que llega el término.

Por qué no se puede completar a medias un término que la red no lleva#

La lección real de este experimento está antes del conteo de órdenes. $\Phi_{\alpha\beta\gamma}$ es un tensor. D2Q9 concede $xxy$ y $xyy$ , y bloquea $xxx$ y $yyy$ . Completar la mitad de un tensor equivale a insertar una corrección con la isotropía rota. Errores que se estaban cancelando pueden asomar, y aquí asomaron.

Al encontrar un término que la red no puede transportar, hay dos opciones. Ampliar el conjunto de velocidades hasta poseer ese momento (D2Q17, D2Q37), o no tocarlo en absoluto. No hay intermedio. Antes de sumar una línea a un término de fuerza, conviene contar si el resto de esa línea existe en la red.