La velocidad salió un sexto con τ=0.6: omitir la media fuerza F/2 en el forcing de Guo

La misma fuerza, ¿por qué la velocidad es un sexto de la respuesta?#

Se aplicó una fuerza de volumen sinusoidal a una caja periódica. La solución analítica ocupa una línea en cualquier libro. Sin embargo, el código de lattice Boltzmann (LBM, un método que hace colisionar y desplazar funciones de distribución para resolver el flujo) devolvió una amplitud de 0.16 veces esa respuesta. Al cambiar la viscosidad, la razón también cambiaba. Con τ=1 era exactamente la mitad.

La causa fue una sola línea: faltaba sumar $\mathbf{F}/2$ en la velocidad macroscópica. Este artículo muestra por qué esa omisión produce un error de exactamente $(1-1/2\tau)$. Se mide en el flujo de Kolmogorov (flujo cortante periódico impulsado por una fuerza de volumen sinusoidal) y se cierra con una prueba de un paso que detecta el error en cualquier código LBM en unos 30 segundos.

Dentro de un paso, el momento cambia tres veces#

Un paso de LBM es colisión más transmisión. La transmisión solo traslada momento. El momento se crea únicamente en la colisión. Por eso basta con el balance de un solo nodo de fluido uniforme.

$\mathbf{m}=\sum_i f_i\mathbf{c}_i$ es el momento antes de la colisión, $\mathbf{m}'$ el momento después, $\rho\mathbf{u}$ el momento que se coloca en el equilibrio, $\tau$ el tiempo de relajación y $\mathbf{F}$ la fuerza de volumen en unidades de red.

El segundo término es la relajación. Se anula si la velocidad de equilibrio coincide con la velocidad previa a la colisión. El tercer término es el primer momento de la fuente de Guo. Lo que se ponga en $\rho\mathbf{u}$ lo decide todo.

En la primera línea, la relajación aporta el $\mathbf{F}/(2\tau)$ que falta y entra exactamente $\mathbf{F}$. En la segunda, ese complemento desaparece. Solo sobrevive el factor $(1-1/2\tau)$ de la fuente.

La simulación de abajo permite manipularlo directamente.

Al llevar τ hacia 0.51, la pendiente roja se aplana casi a cero y la parte sólida de relajación desaparece de la barra derecha.

La velocidad de medio paso: justo entre el antes y el después#

¿Por qué precisamente $\mathbf{F}/2$? Antes de la colisión el momento es $\mathbf{m}$. Tras recibir toda la fuerza es $\mathbf{m}+\mathbf{F}$. La velocidad física es el promedio en ese intervalo, es decir, el punto medio. Es la misma razón por la que funciona la regla del trapecio (integrar con el promedio de los valores extremos, con precisión de segundo orden).

Guo et al. (2002) escribieron el término fuente así:

$w_i$ son los pesos de la red, $\mathbf{c}_i$ las velocidades discretas y $c_s^2=1/3$ el cuadrado de la velocidad del sonido de la red. Cumple $\sum_i S_i\mathbf{c}_i = \mathbf{F}$ y $\sum_i S_i = 0$.

El factor $(1-1/2\tau)$ sale de la expansión de Chapman–Enskog (desarrollar la distribución en el número de Knudsen para recuperar las ecuaciones macroscópicas). La expansión de Taylor en tiempo discreto deja un residuo de la forma $-\tfrac{1}{2}S_i$. Ese residuo comparte raíz con la corrección de viscosidad de medio paso. Para cancelarlo hay que poner $\mathbf{F}/2$ en la velocidad de equilibrio y así equilibrar el primer momento. Ambos se derivan como pareja desde el principio.

Código: flujo de Kolmogorov con tres definiciones de velocidad#

Se arma una red D2Q9 como columna periódica de $1\times N$. La fuerza es $F_x = F_0\sin(ky)$ con $k=2\pi/N$. La solución analítica estacionaria es

donde $\nu$ es la viscosidad cinemática en unidades de red. No hay paredes, así que el error del rebote (bounce-back) no contamina el resultado. Solo se mide el tratamiento de la fuerza.

Se ejecutan tres casos. guo usa $\mathbf{F}/2$ en el equilibrio y en la salida. nohalf conserva la fuente de Guo pero quita $\mathbf{F}/2$ de la velocidad. rawout calcula correctamente pero reporta $\sum f_i\mathbf{c}_i/\rho$. La amplitud es el coeficiente de Fourier de $\sin(ky)$ dividido entre el valor analítico.

import math
 
# D2Q9 lattice: velocities, weights, sound speed squared
C = [(0,0),(1,0),(0,1),(-1,0),(0,-1),(1,1),(-1,1),(-1,-1),(1,-1)]
W = [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]
CS2 = 1/3
 
def f_equil(rho, ux):
    # uy = 0 everywhere: the flow is u_x(y) only
    out = []
    for (cx, cy), wi in zip(C, W):
        cu = cx*ux
        out.append(wi*rho*(1 + cu/CS2 + 0.5*(cu/CS2)**2 - 0.5*ux*ux/CS2))
    return out
 
def guo_term(ux, Fx, tau):
    out = []
    for (cx, cy), wi in zip(C, W):
        cu = cx*ux
        out.append((1 - 0.5/tau)*wi*((cx - ux)/CS2 + cu*cx/CS2**2)*Fx)
    return out
 
def kolmogorov_amp(tau, mode, N=16, F0=1e-6):
    """mode 'guo'    : u = (sum f c + F/2)/rho in equilibrium and output
       mode 'nohalf' : Guo source, but u = sum f c / rho everywhere
       mode 'rawout' : correct update, but report sum f c / rho"""
    k = 2*math.pi/N
    Fx = [F0*math.sin(k*(j + 0.5)) for j in range(N)]
    nu = (tau - 0.5)/3
    steps = int(12/(nu*k*k))                    # ~12 viscous decay times
    f = [f_equil(1.0, 0.0) for _ in range(N)]
    for _ in range(steps):
        post = []
        for j in range(N):
            fj = f[j]
            rho = sum(fj)
            mx = sum(fi*cx for fi, (cx, _) in zip(fj, C))
            ux = mx/rho if mode == 'nohalf' else (mx + 0.5*Fx[j])/rho
            feq = f_equil(rho, ux)
            src = guo_term(ux, Fx[j], tau)
            post.append([fi - (fi - fe)/tau + s for fi, fe, s in zip(fj, feq, src)])
        # streaming: only the y component moves anything (one cell in x, periodic)
        f = [[post[(j - C[i][1]) % N][i] for i in range(9)] for j in range(N)]
    amp = 0.0
    for j in range(N):
        rho = sum(f[j]); mx = sum(fi*cx for fi, (cx, _) in zip(f[j], C))
        u = mx/rho if mode in ('nohalf', 'rawout') else (mx + 0.5*Fx[j])/rho
        amp += 2*u*math.sin(k*(j + 0.5))/N
    return amp/(F0/(nu*k*k)), steps
 
print(" tau  steps   guo     nohalf  rawout   nohalf/guo  1-1/(2tau)")
for tau in [0.55, 0.6, 0.8, 1.0, 1.5]:
    r1, n = kolmogorov_amp(tau, 'guo')
    r2, _ = kolmogorov_amp(tau, 'nohalf')
    r3, _ = kolmogorov_amp(tau, 'rawout')
    print(f"{tau:4.2f} {n:6d}  {r1:.4f}  {r2:.4f}  {r3:.4f}   {r2/r1:.4f}      {1-0.5/tau:.4f}")

 tau  steps   guo     nohalf  rawout   nohalf/guo  1-1/(2tau)
0.55   4668  0.9746  0.0885  0.9734   0.0908      0.0909
0.60   2334  0.9754  0.1621  0.9729   0.1662      0.1667
0.80    778  0.9836  0.3660  0.9759   0.3721      0.3750
1.00    466  1.0001  0.4936  0.9872   0.4936      0.5000
1.50    233  1.0772  0.7010  1.0515   0.6508      0.6667

El código anterior puede ejecutarse directamente en el navegador.

Con τ=0.6, la amplitud roja de nohalf se detiene cerca de 0.16 y se pega a la línea discontinua roja $1-1/2\tau$. Al subir N a 32 o más, la verde guo se acerca a 1, mientras la razón nohalf/guo no se mueve.

Lo que dicen 0.166 y 0.973#

La columna nohalf/guo coincide con $1-1/2\tau$: 0.0908 frente a 0.0909 con τ=0.55, y 0.1662 frente a 0.1667 con τ=0.6. La conclusión del balance de momento aparece tal cual en un flujo real. En estado estacionario, la difusión viscosa equilibra el momento que entra en cada paso. Si lo que entra se escala por $(1-1/2\tau)$, la amplitud también.

Al crecer τ, la razón se aparta un poco de la teoría: 0.6508 frente a 0.6667 con τ=1.5. Con $N=16$, $k\approx0.39$, así que hay pocos nodos por longitud de onda. Ahí se cuela el error de truncamiento de la red de orden $k^2$. Por eso guo marca 1.0772 y no 1.

rawout es casi idéntico a guo. La velocidad reportada solo se retrasa en $F/(2\rho)$. El error relativo es del orden de $\nu k^2/2$, cerca de 0.1% con τ=0.55. Es un punto fácil de confundir al depurar. Omitir $\mathbf{F}/2$ en la salida produce un error pequeño. Omitirlo en el equilibrio produce uno grande.

El régimen más peligroso es el de alto número de Reynolds. Al acercar τ a 0.5 para bajar la viscosidad, la fuerza efectiva también tiende a cero. Si se aumenta la fuerza porque la velocidad no llega, un factor oculto queda en todos los resultados.

Dónde va F/2 en el pseudocódigo#

Un bucle LBM habitual sigue el orden: inicialización → cantidades macroscópicas → equilibrio → fuerza → colisión → transmisión. Ese orden es el que genera el error. Cuando se calcula $\mathbf{u}$ en la etapa macroscópica, la fuerza todavía no existe. Por eso es tentador quedarse en $\mathbf{u}=\sum f_i\mathbf{c}_i/\rho$.

La solución es reordenar un poco.

1. Calcular primero $\rho = \sum_i f_i$.
2. Calcular la fuerza $\mathbf{F}$. Si solo necesita $\rho$, como la gravedad, aquí termina.
3. Obtener la velocidad con $\rho\mathbf{u} = \sum_i f_i\mathbf{c}_i + \mathbf{F}/2$.
4. Construir $f^{eq}$ y $S_i$ con ese mismo $\mathbf{u}$. La salida también usa ese $\mathbf{u}$.

Si la fuerza depende de la velocidad (por ejemplo, un arrastre), el paso 3 se vuelve implícito. Un arrastre lineal se resuelve algebraicamente. En otros casos bastan una o dos iteraciones de punto fijo.

"Corregirlo" borrando $(1-1/2\tau)$ tampoco es la respuesta. El primer momento queda bien, pero en el segundo momento persiste un esfuerzo que depende de τ. Ese problema se trató en el artículo sobre el tercer momento del término de fuerza y en el artículo sobre el forcing de Shan–Chen por desplazamiento de velocidad.

Una prueba de un paso para encontrar F/2 en un bucle LBM propio#

Se arma una caja periódica con $\rho=1$ y $\mathbf{u}=0$. Se aplica una fuerza uniforme $\mathbf{F}$ y se ejecuta exactamente un paso. Luego se divide el cambio de $\sum f_i\mathbf{c}_i$ entre $\mathbf{F}$.

Un resultado de 1 indica que la pareja es coherente. Un resultado de $1-1/2\tau$ indica que falta $\mathbf{F}/2$ en la velocidad de equilibrio. Repetirlo con dos valores de τ da certeza. Si la razón sigue a τ, ahí está el error. Conviene revisar este único paso antes de correr cualquier caso de validación.