Dar a cada especie su propio tau desvió la masa total un 16.8% — velocidad de corrección en LBM multicomponente

En cuanto dos especies reciben difusividades distintas, la continuidad se rompe. Una línea de velocidad de corrección convierte 0.168 en 0.0005.

El día en que la tercera especie descuadró el libro de densidades#

Cada especie que se agrega a un solver de combustión trae consigo otro juego de funciones de distribución. Cada una difunde a su propio ritmo, así que cada una recibe su propio tiempo de relajación $\tau_k$ . Hasta aquí, manual de texto. Tras unos cientos de pasos, $\sum_k Y_k$ se aleja de 1. Las fracciones másicas (la porción de la masa total que lleva cada especie) tienen que sumar 1 por definición.

Este artículo rastrea esa desviación hasta su origen con un desarrollo de Chapman–Enskog y luego la deshace con una sola línea. Una mezcla ternaria corre 2000 pasos y el déficit máximo cae de 0.168 a 0.0005. La conclusión por delante: no es un bug, es el modelo. En cuanto las especies reciben $\tau_k$ distintos, la continuidad está construida para romperse.

Tres juegos de funciones de distribución, repartidos por tarea#

Un LBM compresible multifásico multicomponente separa sus distribuciones según su función. $f_i$ lleva el campo de flujo, $g_i$ la energía, $h_{k,i}$ la especie $k$ -ésima. Las tres viven sobre las mismas velocidades de red $\mathbf{c}_i$ y resuelven la misma ecuación de Boltzmann discreta.

h_{k,i}(\mathbf{x}+\mathbf{c}_i\delta t,\;t+\delta t)-h_{k,i}(\mathbf{x},t) =-\frac{1}{\tau_k}\Big[h_{k,i}(\mathbf{x},t)-h_{k,i}^{\mathrm{eq}}(\mathbf{x},t)\Big]

$h_{k,i}$ es la distribución de la especie $k$ en la dirección $i$ , $\tau_k$ su tiempo de relajación y $\delta t$ el paso temporal.

El equilibrio solo necesita la densidad de especie $\rho_k$ y la velocidad del flujo $\mathbf{u}$ . A diferencia de $f_i$ , no lleva término de segundo orden, porque la ecuación de especies solo debe recuperar advección y difusión.

h_{k,i}^{\mathrm{eq}}=\rho_k w_i\left(1+\frac{\mathbf{c}_i\cdot\mathbf{u}}{c_s^2}\right)

$w_i$ son los pesos de red y $c_s$ la velocidad del sonido de red. La densidad de especie se lee como $\rho_k=\sum_i h_{k,i}$ y la fracción másica como $Y_k=\rho_k/\rho$ .

Al separar los deslizadores $D_A,D_B,D_C$ , las tres difusividades se abren y la suma de masa del panel inferior se aparta de 1.

D_A = 0.005 D_B = 0.100 D_C = 0.500  correction velocity

step 0 — current |sum Y_k - 1| = 0.00e+0, peak = 0.00e+0. Drag the three sliders to the same value and the lower panel flattens to round-off. Spread them apart and the error grows every step: each species relaxes with its own tau, so the three Fickian fluxes no longer cancel. Tick the checkbox to inject the correction velocity and the growth stops on the spot.

Al alinear los tres deslizadores, el panel inferior se aplana hasta el redondeo. Dejarlos separados hace que el error se acumule en cada paso. La casilla es la velocidad de corrección. Al activarla, la suma se pliega otra vez hacia 1.

Lo que Chapman–Enskog extrae de $h_k$ #

Se expande la ecuación discreta en serie de Taylor en $\delta t$ y se separa la distribución como $h_{k,i}=h_{k,i}^{\mathrm{eq}}+\epsilon h_{k,i}^{(1)}+\cdots$ . Al recolectar el momento de orden cero por órdenes aparece una ecuación de advección–difusión para una sola especie.

\partial_t \rho_k+\nabla\cdot(\rho_k\mathbf{u})=\nabla\cdot\big(D_k\nabla\rho_k\big), \qquad D_k=c_s^2\Big(\tau_k-\tfrac{1}{2}\Big)\delta t

$D_k$ es la difusividad de esa especie. El $-1/2$ proviene de la discretización, no de la física. Ese mismo $-1/2$ aparece como viscosidad en el campo de flujo y como difusividad térmica en el de energía. Al atar ambos libros contables, el número de Prandtl queda fijado a una constante — eso está en las dos constantes escondidas en el LBM térmico.

Hasta aquí todo se ve limpio. Basta girar un $\tau_k$ para que la difusividad caiga donde se quiera. Si de verdad ocurre así se mide más abajo.

Qué queda cuando los $\tau_k$ difieren#

El problema no aparece en una ecuación de especie sino en su suma. Se suma la ecuación sobre $k$ y se usa $\sum_k\rho_k=\rho$ .

\partial_t \rho+\nabla\cdot(\rho\mathbf{u})=\nabla\cdot\left(\sum_k D_k\nabla\rho_k\right)

El lado izquierdo es la continuidad. El derecho debe anularse. En general no lo hace.

\sum_k D_k\nabla\rho_k=\rho\sum_k D_k\nabla Y_k \neq 0

Solo cuando todos los $D_k$ son iguales se cumple $\sum_k D_k\nabla Y_k=D\nabla\sum_k Y_k=D\nabla(1)=0$ . Dar a cada especie su propio $\tau_k$ es exactamente el acto de romper esa cancelación. Escrito con flujos queda aún más claro: los flujos de Fick $\mathbf{J}_k=-\rho D_k\nabla Y_k$ no suman cero.

Físicamente es un hecho conocido. La ley de Fick es una aproximación de mezcla binaria, y con tres o más especies el cierre correcto son las relaciones de Stefan–Maxwell. El LBM hereda el déficit completo al transcribir esa aproximación como un $\tau_k$ independiente por especie.

Dos mil pasos de una mezcla ternaria, y el cuadre de cuentas#

Coloqué tres juegos de poblaciones de especie sobre una red D1Q3. La mezcla está en reposo: $\rho=1$ , $\mathbf{u}=0$ . No ocurre nada salvo difusión, de modo que $\sum_k\rho_k$ debe quedarse en 1 para siempre. Primero un modo seno amortiguado verifica si $D_k$ realmente vale $c_s^2(\tau_k-1/2)$ , y después se mide el déficit con varias combinaciones de $D_k$ .

import numpy as np
 
N, CS2 = 256, 1.0/3.0
W = np.array([2.0/3.0, 1.0/6.0, 1.0/6.0])   # pesos D1Q3: c = 0, +1, -1
C = np.array([0, 1, -1])
 
def feq_species(rho_k, u):
    return rho_k[None, :] * W[:, None] * (1.0 + C[:, None] * u[None, :] / CS2)
 
def correction_velocity(rho_k, D):
    dY = np.array([(np.roll(r, -1) - np.roll(r, 1)) * 0.5 for r in rho_k])
    return np.sum(D[:, None] * dY, axis=0)           # V_c = sum_k D_k grad Y_k
 
def lbm_species_step(h, D, use_correction):
    rho_k = h.sum(axis=1)                            # (K, N)
    u = correction_velocity(rho_k, D) if use_correction else np.zeros(N)
    tau = D / CS2 + 0.5
    for k in range(len(D)):
        h[k] += (feq_species(rho_k[k], u) - h[k]) / tau[k]
        for i in range(3):
            h[k, i] = np.roll(h[k, i], C[i])
    return h
 
def diffusion_leak_run(D, steps, use_correction):
    x = np.arange(N)
    r1 = 0.30 + 0.20 * np.sin(2 * np.pi * x / N)
    r2 = 0.45 + 0.15 * np.cos(4 * np.pi * x / N)
    r3 = 1.0 - r1 - r2                               # al inicio la suma es exactamente 1
    h = np.array([feq_species(r, np.zeros(N)) for r in [r1, r2, r3]])
    worst = 0.0
    for _ in range(steps):
        h = lbm_species_step(h, D, use_correction)
        worst = max(worst, np.abs(h.sum(axis=(0, 1)) - 1.0).max())
    return worst
 
def measure_D_eff(tau, steps=400):
    x = np.arange(N)
    amp0 = 0.01
    rho = 0.5 + amp0 * np.sin(2 * np.pi * x / N)     # inyectar un modo seno y medir su decaimiento
    h = feq_species(rho, np.zeros(N))
    for _ in range(steps):
        r = h.sum(axis=0)
        h += (feq_species(r, np.zeros(N)) - h) / tau
        for i in range(3):
            h[i] = np.roll(h[i], C[i])
    amp = np.abs(np.fft.rfft(h.sum(axis=0))[1]) * 2.0 / N
    k2 = (2 * np.pi / N) ** 2
    return -np.log(amp / amp0) / (k2 * steps)
 
print("tau    D_theory=cs^2(tau-1/2)   D_measured   error[%]")
for tau in (0.55, 0.80, 1.20, 2.00):
    dth, dm = CS2 * (tau - 0.5), measure_D_eff(tau)
    print(f"{tau:.2f}   {dth:.6f}                {dm:.6f}     {100*(dm-dth)/dth:+.3f}")
 
print()
print("D1,D2,D3            correction   max|sum(rho_k)-1|")
for D in (np.array([0.05, 0.05, 0.05]),
          np.array([0.05, 0.10, 0.02]),
          np.array([0.01, 0.10, 0.30]),
          np.array([0.005, 0.10, 0.50])):
    for corr in (False, True):
        print(f"{D[0]:.3f},{D[1]:.3f},{D[2]:.3f}     {str(corr):5s}        "
              f"{diffusion_leak_run(D, 2000, corr):.3e}")

tau    D_theory=cs^2(tau-1/2)   D_measured   error[%]
0.55   0.016667                0.016874     +1.243
0.80   0.100000                0.100137     +0.137
1.20   0.233333                0.233122     -0.091
2.00   0.500000                0.498135     -0.373
 
D1,D2,D3            correction   max|sum(rho_k)-1|
0.050,0.050,0.050     False        1.865e-13
0.050,0.050,0.050     True         2.091e-13
0.050,0.100,0.020     False        5.050e-02
0.050,0.100,0.020     True         5.201e-05
0.010,0.100,0.300     False        1.155e-01
0.010,0.100,0.300     True         1.354e-04
0.005,0.100,0.500     False        1.683e-01
0.005,0.100,0.500     True         5.215e-04

Al separar un deslizador $D_k$ , cambia la longitud de la flecha roja resultante.

D_A = 0.005 D_B = 0.100 D_C = 0.500 dY_A/dx = 0.40 dY_B/dx = -0.25  correction velocity

sum_k J_k = 0.00e+0, V_c = 0.0000. With three different D_k the red resultant almost never lands on zero, whatever the gradients — that residue is exactly the mass the solver invents each step. Switch the correction on and the green segments rebalance the three arrows so the resultant collapses to zero for every slider setting.

Con cualquier combinación de $D_k$ y gradientes, la flecha resultante roja casi nunca cae en cero. Al marcar la casilla aparecen los segmentos verdes y la resultante se pliega hasta cero exacto.

Conviene leer primero la tabla de arriba. $D_k=c_s^2(\tau_k-1/2)$ se cumple al nivel del 0.1%. El 1.24% en $\tau=0.55$ es el error de red que crece a medida que $\tau$ se acerca a $1/2$ . La relación en sí es confiable.

La tabla de abajo es el punto. Con las tres difusividades en 0.05, el déficit vale $1.9\times10^{-13}$ , redondeo de doble precisión. Separarlas por un factor 2.5 ya da $5\times10^{-2}$ . Separarlas por 100 lleva a 0.168. Con $\rho=1$ , eso es el 16.8% de la masa total.

Al meter $c_s^2=1/3$ en $D_k=c_s^2(\tau_k-1/2)$ sale la columna teórica de arriba. Con $\tau=0.80$ , $(0.30)/3=0.100000$ . La medida es 0.100137, error $+0.137\%$ . Con $\tau=2.00$ , $(1.50)/3=0.500000$ , medida 0.498135, error $-0.373\%$ . El signo cambia. Con $\tau=0.55$ , $(0.05)/3=0.016667$ frente a una medida de 0.016874, $+1.243\%$ . Con $\tau=1.20$ la teoría es 0.233333, la medida 0.233122, $-0.091\%$ . Cuanto más lejos de $1/2$ , más chico el error porcentual.

Los cocientes del déficit también salen de la tabla inferior. Para $D=(0.050,\,0.100,\,0.020)$ , $5.050\times10^{-2}$ sobre $5.201\times10^{-5}$ da 971. Para $D=(0.005,\,0.100,\,0.500)$ , $1.683\times10^{-1}$ sobre $5.215\times10^{-4}$ da 323. Un factor 100 entre difusividades compra una corrección menor que un factor 2.5. Con los tres valores en 0.05 el déficit es $1.865\times10^{-13}$ , y con la corrección activa queda en $2.091\times10^{-13}$ . La cancelación ya vale, así que la velocidad de corrección no tiene trabajo. En 2000 pasos con $\rho=1$ , 0.168 es el 16.8% de la masa, y $5.215\times10^{-4}$ es el residuo de ese hueco.

La línea que convirtió 0.168 en 0.0005#

El arreglo se conoce desde hace mucho. Se suma a todas las especies la misma velocidad de corrección $\mathbf{V}_c$ , forzando que los flujos sumen cero.

\mathbf{J}_k=-\rho D_k\nabla Y_k+\rho Y_k\mathbf{V}_c, \qquad \mathbf{V}_c=\sum_k D_k\nabla Y_k

Con ese $\mathbf{V}_c$ queda $\sum_k\mathbf{J}_k=-\rho\sum_k D_k\nabla Y_k+\rho\mathbf{V}_c\sum_k Y_k=0$ . Usa directamente la restricción $\sum_k Y_k=1$ .

En LBM esto nunca se convierte en un término fuente. Basta cambiar la velocidad dentro del equilibrio de $\mathbf{u}$ a $\mathbf{u}+\mathbf{V}_c$ : la línea correction_velocity y la rama use_correction de arriba son todo el cambio. El flujo de más que aparece al subir la velocidad vive en el término sobrante de un LBM de convección–difusión.

La tabla dice lo mismo. Un déficit de $5.05\times10^{-2}$ pasa a $5.20\times10^{-5}$ , y $1.68\times10^{-1}$ pasa a $5.22\times10^{-4}$ : factores de 971 y 323.

Por qué la corrección aún no llega a cero#

Hay un detalle que salta a la vista. Con difusividades iguales el déficit es $1.9\times10^{-13}$ , pero con la corrección activa se estanca cerca de $5\times10^{-4}$ . Eso no es precisión de máquina.

La razón es que hay dos operadores discretos en juego. correction_velocity construye $\nabla Y_k$ con una diferencia central. El flujo difusivo que el LBM produce realmente es el primer momento de la parte de no equilibrio $h_{k,i}^{(1)}$ , es decir, otro gradiente discreto definido por el streaming. Ambos discrepan en $O(\delta x^2)$ , y esa discrepancia sobrevive en vez de cancelarse.

Para eliminarla conviene tomar la velocidad de corrección del mismo momento de no equilibrio, algo como $\mathbf{V}_c=\frac{1}{\rho}\sum_k\big(1-\frac{1}{2\tau_k}\big)\sum_i\mathbf{c}_i h_{k,i}^{(1)}$ . La diferencia central también va un paso atrasada, lo que deja un residuo por el mismo motivo. La idea general — que una restricción discreta transcrita tal cual del continuo no se satisface de forma exacta — es el mismo tipo de hueco que el defecto del momento cúbico del equilibrio LBM.

En la práctica, eso sí, $5\times10^{-4}$ y $1.7\times10^{-1}$ son cosas muy distintas. El primero no se acumula y baja al refinar la malla. El segundo crece con el tiempo.

Dónde asoma este término en un código de combustión#

Los números de Lewis en combustión metano–aire se dispersan entre 0.3 y 1.8. H y H $_2$ son extremadamente rápidos; los hidrocarburos grandes, lentos. En difusividad eso es un rango de más de veinte. La tercera fila de la tabla cae justo ahí.

El déficit no se queda en un error de densidad. Las tasas de reacción dependen de $Y_k$ de forma exponencial. Si la suma pasa de 1, la tasa se sobreestima, ese calor vuelve a empujar la densidad y la corrida diverge. Muchos códigos renormalizan dividiendo $Y_k$ por la suma en cada paso, pero eso no elimina el déficit: lo reparte de manera arbitraria entre especies. Qué especie paga la cuenta depende del índice del nodo.

El remedio más común en producción es designar la especie más pesada como especie de relleno (bath species) y calcular $Y_{N_2}=1-\sum_{k\neq N_2}Y_k$ de forma directa. La suma queda exactamente en 1, pero todo el déficit se apila sobre el N $_2$ . La velocidad de corrección reparte esa carga en proporción a $Y_k$ .

Volcar el déficit sobre N $_2$ cuadra el libro en la especie equivocada#

Al dar a cada especie su propio $\tau_k$ se abren tres caminos. Dejar el déficit donde está, dividir por la suma y cargarlo sobre una especie, o repartirlo con una velocidad de corrección en proporción a $Y_k$ . Dejarlo no es un camino. Cargarlo sobre una especie solo se sostiene cuando esa especie de relleno pesa mucho, combustión diluida con $Y_{N_2}\approx 0.7$ . El valor por defecto es la velocidad de corrección. El costo es un juego de gradientes por paso y una velocidad cambiada en el equilibrio.

El rango de $\tau$ se cierra al mismo tiempo. La tabla muestra un error de difusividad de $+1.243\%$ en $\tau=0.55$ y de $-0.373\%$ en $\tau=2.00$ . El signo cambia. Ese piso es el mismo sitio que el medio paso que lee mal la viscosidad por un factor 51 en $\tau=0.51$ . A lo largo del eje de especies el piso se tiende una vez más. Cuando un registro imprime $\sum_k Y_k=1.0003$ , conviene alinear primero los $D_k$ . Si difieren, la causa ya está en esa línea.