Solo cambió τ y cambió la densidad de un gas en reposo — el esfuerzo (τ−½)² del forcing de Shan–Chen por desplazamiento de velocidad

Cambió la viscosidad y cambió un fluido que no se mueve#

Un gas en reposo no sabe cuál es su viscosidad. Con velocidad nula, el esfuerzo viscoso es nulo. Por eso un estado de reposo donde fuerza y presión se equilibran no debería depender de la viscosidad. En el método de Boltzmann en red (LB, Lattice Boltzmann), la viscosidad la fija el tiempo de relajación $\tau$. Cambiar $\tau$ no debería alterar, entonces, el perfil de densidad en reposo.

Hay una familia de esquemas de forcing que rompe esto. Es el forcing por desplazamiento de velocidad (velocity shift) que usan los modelos multifásicos de Shan–Chen. Misma fuerza, misma red, solo $\tau$ de 0.8 a 3, y el perfil de densidad se mueve. Este artículo compara tres esquemas de forcing lado a lado y localiza el origen de la diferencia. El resultado cabe en una línea por esquema: el coeficiente del esfuerzo que cada uno añade en silencio al flujo de momento. Guo da 0, Shan–Chen da $(\tau-\tfrac12)^2$ y EDM da $\tfrac14$.

Criterio 1: ¿hacia dónde se desplaza la velocidad de equilibrio?#

Los tres esquemas añaden un momento $F\,\delta t$ por paso. Difieren en qué distribución lo transporta. Abajo, $u'=\sum_i f_i c_i/\rho$ es la velocidad leída directamente de las poblaciones. En los tres, la velocidad física se define como $u = u' + F/(2\rho)$ (corrección de medio paso).

En reposo $u=0$, así que $u'=-F/(2\rho)$. La velocidad que Shan–Chen introduce en su equilibrio queda

$u_{eq}$ es la velocidad con la que se construye el equilibrio, $\tau$ el tiempo de relajación adimensional y $F$ la fuerza por unidad de volumen. El gas está quieto, pero su equilibrio se mueve con una velocidad proporcional a $\tau$.

El segundo momento del equilibrio (el flujo de momento) es una parábola en esa velocidad:

$c_s^2=1/3$ es el cuadrado de la velocidad del sonido de la red. Cualquier $u_{eq}$ distinto de cero eleva el flujo en $\rho u_{eq}^2$. La simulación de abajo permite manipularlo directamente. Un cursor recorre la columna y la fuerza local $F(y)$ coloca el punto de equilibrio de cada esquema sobre la parábola.

Al llevar $\tau$ desde cerca de 0.5 hasta 3, el punto naranja (Shan–Chen) sube por la parábola mientras el verde (Guo) se queda en el fondo. Los puntos morados (EDM) se ubican en $\pm F/(2\rho)$ y mantienen la misma altura $F^2/(4\rho)$ sin importar $\tau$.

Criterio 2: ¿sobrevive esa elevación hasta el estado de reposo?#

Un equilibrio elevado no implica de inmediato un esfuerzo macroscópico elevado. Lo que transporta momento es el segundo momento de las poblaciones post-colisión, $\Pi^*$. Igualando los segundos momentos antes y después de la colisión en estado estacionario se obtiene

donde $S_\Pi$ es el segundo momento de la fuente adicional. La fuente de Guo tiene segundo momento nulo en $u=0$. La de EDM da $\rho(u'+F/\rho)^2-\rho u'^2 = 0$. Así, los tres esquemas difieren solo por la velocidad colocada en el equilibrio. El balance de fuerzas en la columna en reposo pasa a ser

$b$ es el coeficiente del esfuerzo espurio. La hidrostática del continuo corresponde a $b=0$.

¿Por qué este término no aparece en las deducciones de libro? La expansión estándar de Chapman–Enskog (una expansión multiescala que ordena las poblaciones por potencias de un parámetro pequeño $\epsilon$) trata $F$ como de primer orden. Entonces $F^2/\rho$ es un término $\epsilon^2$ y se descarta al nivel de Navier–Stokes. Cuando su coeficiente crece como $\tau^2$, ya nada garantiza que sea pequeño. El artículo citado en las notas fuente, Revised Chapman–Enskog analysis for a class of forcing schemes in the lattice Boltzmann method, revisa justamente esa asignación de orden del término de fuerza. Si la fuerza debe entrar como gradiente de presión o como fuerza de cuerpo se trató en las dos formas del forcing en LBM no ideal. Aquí la pregunta es cómo se carga la misma fuerza sobre las poblaciones.

Código: una columna periódica, tres forcings, cinco valores de τ#

Una columna periódica D1Q3 (velocidades 0, ±1) de 32 celdas recibe $F = G\sin(ky)$. Parte con densidad uniforme. Cuando la viscosidad amortigua las ondas sonoras, queda un estado de reposo. Sobre la densidad convergida, un ajuste por mínimos cuadrados de $c_s^2\rho = \lambda P(y) - bF^2/\rho + c$ mide $b$. $P(y)=-\tfrac{G}{k}\cos(ky)$ es el potencial hidrostático con $P'=F$. $\lambda$ absorbe el factor de amplitud que introduce la discretización de la red.

import math
 
# Red D1Q3: velocidades 0, +1, -1 ; cs^2 = 1/3
C = (0.0, 1.0, -1.0)
W = (2/3, 1/6, 1/6)
CS2 = 1/3
 
def d1q3_eq(rho, u):
    return [W[i]*rho*(1 + C[i]*u/CS2 + 0.5*(C[i]*u)**2/CS2**2 - 0.5*u*u/CS2)
            for i in range(3)]
 
def collide_node(fn, F, tau, scheme):
    rho = sum(fn)
    mom = fn[1] - fn[2]
    if scheme == "guo":                 # Guo: u = (sum f c + F/2)/rho, fuente por (1 - 1/2tau)
        u = (mom + 0.5*F)/rho
        feq = d1q3_eq(rho, u)
        out = []
        for i in range(3):
            S = W[i]*((C[i] - u)/CS2 + C[i]*u*C[i]/CS2**2)*F
            out.append(fn[i] - (fn[i] - feq[i])/tau + (1 - 0.5/tau)*S)
        return out
    if scheme == "sc":                  # Shan-Chen: desplaza la velocidad de equilibrio en tau F / rho
        feq = d1q3_eq(rho, mom/rho + tau*F/rho)
        return [fn[i] - (fn[i] - feq[i])/tau for i in range(3)]
    up = mom/rho                        # EDM: suma feq(u'+F/rho) - feq(u')
    feq = d1q3_eq(rho, up)
    fsh = d1q3_eq(rho, up + F/rho)
    return [fn[i] - (fn[i] - feq[i])/tau + fsh[i] - feq[i] for i in range(3)]
 
def run_column(scheme, tau, G, N=32, tol=1e-12):
    k = 2*math.pi/N
    F = [G*math.sin(k*(j + 0.5)) for j in range(N)]
    f = [d1q3_eq(1.0, 0.0) for _ in range(N)]
    for it in range(1, 100001):
        post = [collide_node(f[j], F[j], tau, scheme) for j in range(N)]
        new = [[post[j][0], post[j-1][1], post[(j+1) % N][2]] for j in range(N)]
        if it % 100 == 0 and max(abs(sum(new[j]) - sum(f[j])) for j in range(N)) < tol:
            return F, [sum(fn) for fn in new], it
        f = new
    return F, [sum(fn) for fn in f], it
 
def stress_fit(F, rho, G):
    # ajuste  cs^2 rho = lam*P(y) - b*F^2/rho + c,  P' = F  (hidrostatico: b = 0)
    N = len(F); k = 2*math.pi/N
    X = [[F[j]**2/rho[j], -G/k*math.cos(k*(j + 0.5)), 1.0] for j in range(N)]
    r = [CS2*x for x in rho]
    A = [[sum(X[n][a]*X[n][b] for n in range(N)) for b in range(3)] for a in range(3)]
    v = [sum(X[n][a]*r[n] for n in range(N)) for a in range(3)]
    for i in range(3):
        for j in range(i + 1, 3):
            m = A[j][i]/A[i][i]
            A[j] = [A[j][c] - m*A[i][c] for c in range(3)]
            v[j] -= m*v[i]
    x = [0.0]*3
    for i in (2, 1, 0):
        x[i] = (v[i] - sum(A[i][c]*x[c] for c in range(i + 1, 3)))/A[i][i]
    return -x[0], x[1]
 
G = 0.02
print(f"periodic column N=32, F = {G} sin(ky), rho0 = 1")
print(" tau | b_Guo  b_SC   (tau-1/2)^2  b_EDM | lam   | max|rho_SC-rho_Guo|  iters")
for tau in [0.8, 1.0, 1.5, 2.0, 3.0]:
    out = {s: run_column(s, tau, G) for s in ("guo", "sc", "edm")}
    b = {s: stress_fit(out[s][0], out[s][1], G) for s in out}
    drho = max(abs(p - q) for p, q in zip(out["sc"][1], out["guo"][1]))
    bg = round(b["guo"][0], 3) + 0.0
    print(f"{tau:4.1f} | {bg:5.3f} {b['sc'][0]:6.3f}   {(tau-0.5)**2:6.3f}    {b['edm'][0]:5.3f} |"
          f" {b['guo'][1]:.4f} |  {drho:.2e}            {out['sc'][2]}")

periodic column N=32, F = 0.02 sin(ky), rho0 = 1
 tau | b_Guo  b_SC   (tau-1/2)^2  b_EDM | lam   | max|rho_SC-rho_Guo|  iters
 0.8 | 0.000  0.090    0.090    0.250 | 0.9968 |  5.49e-05            6100
 1.0 | 0.000  0.250    0.250    0.250 | 0.9968 |  1.53e-04            3600
 1.5 | 0.000  1.000    1.000    0.250 | 0.9968 |  6.11e-04            1900
 2.0 | 0.000  2.250    2.250    0.250 | 0.9968 |  1.38e-03            1300
 3.0 | 0.000  6.250    6.250    0.250 | 0.9968 |  3.84e-03            800

Lo que dicen 0, (τ−½)² y 1/4#

El código de arriba puede ejecutarse directamente en el navegador. Las tres columnas avanzan a la vez bajo la misma fuerza y $b$ se reajusta en cada cuadro.

Con los valores por defecto $\tau=1.5$, $G=0.02$, la curva naranja se fija en 1.000, la morada en 0.250 y la verde en 0.000: la tercera fila del bloque de salida. Al subir $\tau$ a 3 y pulsar Reset, solo la naranja sube, hasta 6.25. Cambiar $G$ reescala los puntos como $G^2$, pero $b$ no se mueve.

Lo primero que hay que revisar en la salida es la columna $b_{SC}$ frente a la columna $(\tau-\tfrac12)^2$. Las cinco filas coinciden hasta el tercer decimal. No se perdió nada en la deducción. El 0.000 de Guo también importa: el equilibrio hidrostático se cumple hasta el nivel $F^2$.

En $\tau=1$, Shan–Chen y EDM marcan ambos 0.250. No es casualidad. Con $\tau=1$, las poblaciones post-colisión son el propio equilibrio. Shan–Chen se reduce a $f^{eq}(u'+F/\rho)$ y EDM a $f^{eq}(u') + f^{eq}(u'+F/\rho) - f^{eq}(u')$, que es lo mismo. Los dos esquemas se encuentran en un solo punto.

$\lambda=0.9968$ no tiene que ver con el forcing. Es el factor de discretización que una red de 32 celdas aplica a la amplitud del seno, y coincide con $1-k^2/12$. Tiende a 1 al refinar la malla. El ajuste está planteado para que $b$ quede separado de ese factor.

La diferencia de densidad parece pequeña: $3.84\times10^{-3}$ en $\tau=3$, el 1.3% de la amplitud de densidad 0.30. Pero esa proporción crece en proporción a $G$. El efecto principal escala con $F$ y el esfuerzo espurio con $F^2$.

Los tres esquemas en una tabla#

EDM no depende de $\tau$, pero se desvía en una constante $\tfrac14$, lo que facilita corregirlo. Shan–Chen solo es pequeño cerca de $\tau=0.5$, justo donde LB es menos estable.

En una interfaz multifásica, F es grande#

La fuerza de este experimento era suave y débil. En los modelos multifásicos de pseudopotencial es distinta. La fuerza de interacción $F=-\mathcal{G}\psi\nabla\psi$ ($\mathcal{G}$ es la intensidad de interacción y $\psi$ la densidad efectiva) es grande y se concentra en pocas celdas alrededor de la interfaz. Allí, $bF^2/\rho$ puede alcanzar el mismo orden que la presión de la ecuación de estado.

Por eso, una interfaz plana construida con forcing de Shan–Chen tiene densidades de coexistencia (la densidad de cada fase en equilibrio) que se desplazan con $\tau$. Se toca $\tau$ para cambiar la viscosidad y cambia la densidad del líquido. El experimento de la columna aísla ese efecto sin interfaz alguna. El esfuerzo espurio entra tal cual en la condición de equilibrio mecánico de la interfaz.

Algunos métodos usan este término a propósito. Si $b$ se trata como parámetro ajustable, se convierte en una perilla para acercar las densidades de coexistencia a la consistencia termodinámica: se suma lo necesario sobre el $\tfrac14$ de EDM o el 0 de Guo. En cualquier caso, primero hay que saber qué $b$ tiene el propio código. Otro término oculto hallado con la misma maquinaria de Chapman–Enskog aparece en el cuasi-equilibrio g* del LB de dos poblaciones.

La prueba de 32 celdas antes de confiar en un forcing nuevo#

Al incorporar un esquema de forcing nuevo, conviene no empezar por un caso multifásico. En una interfaz, las causas se mezclan. Primero se corre la columna periódica de arriba. Se mide $b$ con dos valores de $\tau$. Si $b$ cambia con $\tau$, hay un desplazamiento de velocidad escondido en el equilibrio. Si $b$ es una constante distinta de cero, esa constante es lo que se sumará a la presión de la interfaz. Bastan 32 celdas y unos miles de pasos.