Subir la velocidad cuesta el 27% de la difusión — el flujo sobrante del modelo de convección-difusión en LBM

En lattice Boltzmann un solo tiempo de relajación fija la difusividad. $D = c_s^2(\tau - 1/2)\delta t$. Línea corta. Sin velocidad. ¿Sale el mismo $D$ con el flujo encendido? Bajo una advección uniforme $u$ cae a $D(1 - u^2/c_s^2)$. Con $u = 0.3$ se pierde un 27%. La expansión de Chapman–Enskog es donde se escapa el término, y un único término fuente es lo que lo devuelve.

Donde duele es el flujo multifásico con phase-field. Cahn–Hilliard o Allen–Cahn sobre un núcleo de lattice Boltzmann atan el espesor de interfaz a la movilidad. Una movilidad con 27% de error arrastra el espesor, y el espesor arrastra el coeficiente de tensión superficial.

La difusividad está ajustada y aun así la interfaz adelgaza#

Conviene mirarlo primero. Abajo está la ecuación de convección-difusión unidimensional

resuelta con un esquema lattice Boltzmann D1Q3. La condición inicial es una única gaussiana; la solución exacta también es una gaussiana y su anchura crece como $\sigma^2(t) = \sigma_0^2 + 2Dt$. Al empujar u de 0 a 0.30 el pico se estrecha frente a la curva punteada, y tau no mueve esa razón.

Con u = 0.30 la curva sólida (calculada) queda más aguda que la punteada (exacta). La traza de $\sigma^2$ no alcanza la pendiente punteada $2D$. Con tau en 0.6 o en 1.5 la razón de déficit se queda en la misma casilla.

D1Q3 solo tiene tres momentos que respetar#

A diferencia de un LBM que resuelve Navier–Stokes, una distribución $g_i$ construida para transporte escalar tiene menos momentos que satisfacer. El equilibrio que usó Guo en su modelo de 2009 para ecuaciones de convección-difusión no lineales es este:

donde $w_i$ son los pesos de la red, $\mathbf{c}_i$ las velocidades de la red y $c_s^2 = 1/3$ el cuadrado de la velocidad del sonido reticular. Esta distribución cumple tres condiciones de momento.

El tercero es donde esto se separa de un equilibrio de flujo: no hay término $\phi\mathbf{u}\mathbf{u}$. El término de segundo orden en la velocidad nunca se puso. Una ecuación escalar no tiene tensor de esfuerzos, así que basta la parte isótropa del segundo momento, y D2Q5 puede sustituir a D2Q9. El medio paso dentro de $\tau - 1/2$, delante de la difusividad, es el mismo resto que el desarrollo de segundo orden de la LBE discreta deja en la viscosidad.

Parece suficiente, y con $u = 0$ lo es exactamente. El problema es que el truncamiento no sale gratis.

El cuarto término que Chapman-Enskog deja atrás#

Escribamos la ecuación de lattice Boltzmann BGK con una fuente $S_i$ añadida:

Se expande $g_i = g_i^{(0)} + \varepsilon g_i^{(1)} + \varepsilon^2 g_i^{(2)}$ y se separa la derivada temporal como $\partial_t = \varepsilon\partial_{t_1} + \varepsilon^2\partial_{t_2}$. El momento de orden cero a orden $\varepsilon$ devuelve intacta la parte advectiva de la ecuación objetivo.

El primer momento de esa misma ecuación de orden $\varepsilon$ es donde se decide todo, porque fija el flujo que transporta $g^{(1)}$.

El segundo término del corchete, $c_s^2\nabla\phi$, es el flujo difusivo que se pedía. El primero viaja con él de polizón. En un esquema de flujo, el término $\rho\mathbf{u}\mathbf{u}$ del segundo momento del equilibrio lo cancela casi por completo; el modelo escalar no tiene ese término. Así que sobrevive.

Reuniendo el momento de orden cero a orden $\varepsilon^2$ se obtiene la forma final.

$D = c_s^2(\tau - 1/2)\delta t$ aparece como se esperaba. El segundo término de la derecha es un flujo que nadie encargó. Se anula en un estado estacionario genuino, con $u$ constante en el tiempo y $\phi$ ya sin cambiar. Pero mientras $\phi$ siga moviéndose —mientras el cálculo siga corriendo— $\partial_t(\phi\mathbf{u})$ no es cero.

Con $u$ uniforme ese término es difusión negativa#

Tomemos el caso más simple: $u$ constante en espacio y tiempo. Entonces $\partial_t(\phi u) = u\,\partial_t\phi$ y, a orden dominante, $\partial_t\phi \simeq -u\,\partial_x\phi$. Sustituyendo,

Un término de difusión con el signo invertido. Sumado al físico,

El déficit es un solo cuadrado. El cuadrado de la velocidad del sonido reticular es un tercio, así que la velocidad 0.3 al cuadrado da 0.09 y el triple da 0.27. Con tiempo de relajación 1, la difusividad marcada es un tercio por un medio, 0.166667. Lo que queda es 0.166667 por 0.73, es decir 0.121667. La brecha de 0.045 es la difusividad marcada por 0.27. La fila $u = 0.30$ de la salida es esa resta. A baja velocidad el cuadrado se encoge. La velocidad 0.05 al cuadrado da 0.0025, el triple da 0.0075 y el déficit se queda en 0.75%. El pico se pega al lado de la curva exacta. El ojo lo pasa por alto.

La razón no sigue al tiempo de relajación por el mismo renglón. El factor de tiempo de relajación menos un medio multiplica a la vez la difusión física y el término sobrante. Al dividir ambos, el factor se va. La razón no depende de $\tau$. Subir $\tau$ de 0.6 a 1.5 lleva $D$ de 0.033333 a 0.333333, un factor de diez. A velocidad 0.25, tres veces 0.0625 es 0.1875, y uno menos eso es 0.8125, fijo en las cuatro filas. Solo crece la brecha absoluta, de cerca de 0.006 a cerca de 0.062. No hay sitio donde enterrar el error bajo una difusividad mayor.

Cuando la velocidad alcanza la del sonido reticular, el inverso de la raíz cuadrada de un tercio, cerca de 0.577, el paréntesis llega a cero. Más allá, $D_{\text{eff}}$ es negativo. El pico se aguza solo en lugar de extenderse, y la distribución diverge.

$1 - 1/(2\tau)$ vuelve a aparecer#

Toca fijar el coeficiente de $S_i$. Para que el término espurio y el término fuente se cancelen en la ecuación anterior,

Solo hay una forma simple que cumpla esto manteniendo además $\sum_i S_i = 0$.

Ahí está otra vez $1 - 1/(2\tau)$. Este factor ocupa el mismo sitio que el límite de momento tercero que el término de fuerza tiene que respetar. En una red de tiempo discreto una fuente actúa dos veces —una a través de $g^{(1)}$ y otra a través del término de segundo orden de la expansión de Taylor— y el factor $1/2$ es el residuo de esa doble contabilidad.

¿Qué pasa si se omite el coeficiente y se pone $\sum_i \mathbf{c}_i S_i = \partial_t(\phi\mathbf{u})$? Con $\tau = 1$ se aplica exactamente el doble de corrección, y un déficit del 27% se convierte en un exceso del 27%. La magnitud del error no cambia, así que una gráfica log-log de convergencia no lo delata.

Calcular $\partial_t(\phi\mathbf{u})$ en código consiste en guardar $\phi u$ del paso anterior y hacer una diferencia hacia atrás. Un arreglo extra es todo el coste.

Midiendo $D_{\text{eff}}$ en 60 líneas de Python#

En lugar de argumentar, conviene medir. Se deja advectar una gaussiana, se ajusta por mínimos cuadrados la pendiente de crecimiento de su segundo momento $\sigma^2$, y esa pendiente es $2D_{\text{eff}}$.

import numpy as np
 
CS2 = 1.0 / 3.0
C = np.array([0, 1, -1])
W = np.array([2 / 3, 1 / 6, 1 / 6])
 
 
def d1q3_equilibrium(phi, u):
    """g_i^eq = w_i phi (1 + c_i u / cs^2) — equilibrio que solo fija tres momentos"""
    return np.stack([W[i] * phi * (1.0 + C[i] * u / CS2) for i in range(3)])
 
 
def gaussian_moments(x, phi):
    m0 = phi.sum()
    mean = (x * phi).sum() / m0
    return mean, (((x - mean) ** 2) * phi).sum() / m0
 
 
def run_cde_lbm(L, steps, tau, u, sigma0, x0, corrected):
    x = np.arange(L, dtype=float)
    phi = np.exp(-((x - x0) ** 2) / (2 * sigma0**2))
    g = d1q3_equilibrium(phi, u)
    phi_old = phi.copy()
    hist = []
    for n in range(steps + 1):
        if n % 100 == 0:
            hist.append((n, gaussian_moments(x, phi)[1]))
        src = np.zeros_like(g)
        if corrected and n > 0:
            # S_i = w_i (1 - 1/(2 tau)) c_i d_t(phi u) / cs^2,  dt = 1
            dt_phiu = (1.0 - 1.0 / (2 * tau)) * u * (phi - phi_old)
            for i in range(3):
                src[i] = W[i] * C[i] * dt_phiu / CS2
        geq = d1q3_equilibrium(phi, u)
        g = g - (g - geq) / tau + src          # colisión
        for i in range(3):
            g[i] = np.roll(g[i], C[i])         # propagación
        phi_old = phi
        phi = g.sum(axis=0)
    return np.array(hist)
 
 
def fit_diffusivity(hist):
    """lee D_eff en la pendiente de sigma^2 = sigma0^2 + 2 D_eff t"""
    return np.polyfit(hist[:, 0], hist[:, 1], 1)[0] / 2.0
 
 
L, STEPS, SIG0, X0 = 800, 1600, 10.0, 80.0
 
tau = 1.0
D = CS2 * (tau - 0.5)
print(f"tau = {tau},  D = cs^2 (tau-1/2) = {D:.6f},  cs^2 = {CS2:.6f}")
print(f"{'u':>6} {'u^2/cs^2':>9} | {'D_eff (no src)':>14} {'ratio':>7} {'1-u^2/cs^2':>11} |"
      f" {'D_eff (src)':>12} {'ratio':>7}")
for u in [0.05, 0.10, 0.20, 0.30]:
    d_raw = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, False))
    d_fix = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, True))
    print(f"{u:>6.2f} {u * u / CS2:>9.4f} | {d_raw:>14.6f} {d_raw / D:>7.4f} {1 - u * u / CS2:>11.4f} |"
          f" {d_fix:>12.6f} {d_fix / D:>7.4f}")
 
u = 0.25
print(f"\nu = {u} fixed, tau sweep   (theory: ratio = 1 - u^2/cs^2 = {1 - u * u / CS2:.4f}, tau-independent)")
print(f"{'tau':>6} {'D':>10} | {'D_eff (no src)':>14} {'ratio':>7} | {'D_eff (src)':>12} {'ratio':>7}")
for tau in [0.6, 0.8, 1.0, 1.5]:
    D = CS2 * (tau - 0.5)
    d_raw = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, False))
    d_fix = fit_diffusivity(run_cde_lbm(L, STEPS, tau, u, SIG0, X0, True))
    print(f"{tau:>6.1f} {D:>10.6f} | {d_raw:>14.6f} {d_raw / D:>7.4f} | {d_fix:>12.6f} {d_fix / D:>7.4f}")

La salida:

tau = 1.0,  D = cs^2 (tau-1/2) = 0.166667,  cs^2 = 0.333333
     u  u^2/cs^2 | D_eff (no src)   ratio  1-u^2/cs^2 |  D_eff (src)   ratio
  0.05    0.0075 |       0.165417  0.9925      0.9925 |     0.166666  1.0000
  0.10    0.0300 |       0.161667  0.9700      0.9700 |     0.166666  1.0000
  0.20    0.1200 |       0.146667  0.8800      0.8800 |     0.166663  1.0000
  0.30    0.2700 |       0.121667  0.7300      0.7300 |     0.166658  0.9999
 
u = 0.25 fixed, tau sweep   (theory: ratio = 1 - u^2/cs^2 = 0.8125, tau-independent)
   tau          D | D_eff (no src)   ratio |  D_eff (src)   ratio
   0.6   0.033333 |       0.027096  0.8129 |     0.033345  1.0004
   0.8   0.100000 |       0.081258  0.8126 |     0.100006  1.0001
   1.0   0.166667 |       0.135417  0.8125 |     0.166661  1.0000
   1.5   0.333333 |       0.270794  0.8124 |     0.333275  0.9998

Al subir u solo crece el flujo fantasma rosa, como $u^2$. Al llevar alpha de 0 a 1 el término fuente borra ese ancho.

Con alpha = 0 el lóbulo rosa pone el flujo físico azul cabeza abajo. Con alpha = 1 el verde cae sobre el rosa y la suma ámbar vuelve a la curva azul. Es la columna con la fuente activa, donde ratio se pega a 1.0000.

En la primera tabla la columna ratio coincide con 1-u^2/cs^2 hasta la cuarta cifra decimal. La predicción no es una estimación: es el orden dominante exacto. Al activar el término fuente, las cuatro velocidades vuelven a 1.0000.

La segunda tabla muerde más fuerte. Barrer $\tau$ de 0.6 a 1.5 cambia $D$ en un factor de diez, y la razón de déficit apenas se mueve de 0.8129 a 0.8124. Intentar enterrar el error bajo una difusividad mayor no funciona: si $D$ se multiplica por diez, la cantidad que desaparece también.

El presupuesto de velocidad reticular ya estaba gastado#

$u^2/c_s^2$ es el cuadrado del número de Mach reticular. Con $u = 0.1$, $0.01 \times 3 = 0.03$: se va un 3% de la difusión. Con $u = 0.2$, $0.04 \times 3 = 0.12$, un 12%. En 0.05 o por debajo la pérdida es 0.75% y se esconde en el error de malla. La costumbre de atar $u$ cerca de 0.1 por la compresibilidad todavía no ha visto la factura del escalar. Si un espesor de interfaz o un número de Sherwood sale como cifra, el término fuente se enciende en cuanto $u$ pasa de 0.1. Por debajo, ahorrar un arreglo deja en pie el 0.9925 de la tabla.

En phase-field no alcanza con fiarse de $u^2/c_s^2$. El término sobrante es el $\partial_t(\phi\mathbf{u})$ completo, y el cuadrado vale solo cuando $u$ es uniforme y estacionaria. Cerca de la interfaz $u$ salta y $\partial_t u$ no es cero, así que también sobrevive la pieza $\phi\,\partial_t u$. $\partial_t(\phi\mathbf{u})$ entra como diferencia hacia atrás de $\phi u$ del paso anterior. El $\tau$ dentro de $1 - 1/(2\tau)$ es el tiempo de relajación del primer momento. Cómo cae la viscosidad al desplazar el marco de velocidad es la misma clase de factura que el momento tercero que le falta al equilibrio. Si la aritmética de la movilidad cuadra y la interfaz igual adelgaza, $\tau$ no era el fallo. El flujo se llevó primero el presupuesto de difusión.