Con tau = 0.51 la viscosidad se leyó 51 veces más alta — el medio paso que deja la LBE discreta

La viscosidad del continuo y la del lattice diferían en un factor de 51#

Mirando solo la ecuación BGK continua, el tiempo de relajación $\tau$ se liga a la viscosidad cinemática $\nu$ como $\nu = c_s^2 \tau$. Con $\tau = 0.51$, $\nu$ debería valer 0.17. Pero al meter ese mismo $\tau$ en un código D2Q9 y medir la tasa de decaimiento de un vórtice de Taylor–Green, sale 0.00333. Un factor de 51.

No es un bug ni un signo equivocado. Es el tiempo de relajación corto por exactamente $\Delta t/2$. Un camino es la expansión de Taylor a segundo orden a lo largo de una característica. El otro es la regla del punto medio sobre la integral de colisión. El deslizador de $\tau$ mueve las tres lecturas de viscosidad sobre el mismo vórtice.

La curva azul es la amplitud en el lattice de $32 \times 32$. La verde sigue $c_s^2(\tau - 1/2)$ y la naranja sigue $c_s^2 \tau$. En $\tau = 0.51$ la verde marca $0.00333$ y la naranja $0.17$.

La línea que corre sobre el lattice no es una ecuación diferencial#

Lo que realmente se ejecuta es esta única línea.

Aquí $f_i$ es la distribución en la dirección de la velocidad discreta $\mathbf{c}_i$, $\Omega_i = f_i - f_i^{\rm eq}$ es la parte de no equilibrio (cuánto se aparta la población del equilibrio) y $\Delta t$ es un paso.

El lado izquierdo no es una derivada. Es una diferencia tomada a través de un hueco finito $\Delta t$ que el código realmente salta. Convertirla en derivadas exige una expansión, y esa expansión tiene un término de segundo orden. La costumbre en CFD es quedarse con el término de primer orden y archivar el resto como "error de truncamiento". Aquí ese término de segundo orden no es un error: reescribe un coeficiente físico.

Contar órdenes muestra por qué no se puede descartar. El término de segundo orden escala como $\Delta t^2 D^2 f / 2$, mientras que el término de colisión de la derecha va dividido por $\Delta t / \tau$. En unidades de lattice, $\tau$ es un número $O(1)$ entre 0.5 y 2, y $\Delta t = 1$. Es decir, el cociente entre segundo y primer orden no se hace pequeño. Congelar en 1 lo que debía tender a cero es el punto de partida de los métodos de gas en lattice, y este término es el precio.

Qué sobrevive al expandir a segundo orden sobre una característica#

Con $D_i \equiv \partial_t + \mathbf{c}_i \cdot \nabla$, el lado izquierdo se vuelve una exponencial de $D_i$. Se conservan los términos hasta segundo orden y se divide por $\Delta t$.

El término $D_i^2$ es el que está en cuestión. Si se descarta, sale $\nu = c_s^2\tau$; si se conserva, sale otra cosa.

Apliquemos una expansión de Chapman–Enskog. Con $f_i = f_i^{(0)} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)}$, $\partial_t = \varepsilon \partial_{t_1} + \varepsilon^2 \partial_{t_2}$, $\nabla = \varepsilon \nabla_1$ y $f_i^{(0)} = f_i^{\rm eq}$, el orden $\varepsilon$ da

Reutilizar esa relación dentro del orden $\varepsilon^2$ es toda la clave. Como $D_{i1}^2 f_i^{(0)} = -\tfrac{1}{\tau} D_{i1} f_i^{(1)}$, queda

Ha aparecido un paréntesis delante de $D_{i1} f_i^{(1)}$. El esfuerzo viscoso sale del segundo momento de $f_i^{(1)}$ multiplicado por $\tau$, así que ese $\tau$ pasa a ser $\tau\left(1 - \frac{\Delta t}{2\tau}\right) = \tau - \frac{\Delta t}{2}$.

En unidades de lattice con $\Delta t = 1$ eso es $c_s^2(\tau - 1/2)$. Como $c_s^2 = 1/3$, un $\tau$ de 0.51 da $\nu = 0.00333$. Las rayas verdes de arriba son exactamente ese número.

Empujar el término sobrante dentro de la distribución#

Hay un segundo camino al mismo resultado que no necesita expansión. Se integra la BGK continua a lo largo de una característica, pero con la regla del trapecio. Llamemos $\lambda$ al tiempo de relajación.

El lado derecho contiene el nivel de tiempo futuro, así que esto es implícito. Tal como está, no corre. Entonces se cambia de variable.

Reescrito en $\bar f_i$, los términos de tiempo futuro se cancelan y desaparecen. Sustituyendo $\bar\Omega_i = \bar f_i - f_i^{\rm eq} = \left(1 + \frac{\Delta t}{2\lambda}\right)\Omega_i$ queda esto.

Completamente explícito, y carácter por carácter la línea que se usa en producción, pero solo cuando $\tau = \lambda + \Delta t/2$. O sea que el $\tau$ del código no es el tiempo de relajación físico. El físico es $\lambda = \tau - \Delta t/2$. No es un asunto de redondeo, es la factura de un cambio de variable. De dónde sale el $\Delta t/2$ pegado al tiempo de relajación está en el medio paso de la integración trapezoidal.

Releer la misma línea como una regla de punto medio#

La regla del trapecio es el caso $\theta = 1/2$ de una regla con peso $\theta$. Barrer $\theta$ muestra por qué justamente la mitad. La amplificación en un paso es $A_\theta(z) = \frac{1 - (1-\theta)z}{1 + \theta z}$ frente a un $e^{-z}$ exacto, con $z = \Delta t/\lambda$.

Al apartarse de $1/2$, el error va como $(\theta - \tfrac12)z^2$ y sube con pendiente 2. Solo en $\theta = 1/2$ baja a pendiente 3. La mitad no es una perilla: es el único valor que engrana con el $1/2$ del coeficiente de segundo orden de $e^{-z}$.

Un vórtice de Taylor–Green decaído con seis valores de tau#

Hasta aquí llega el argumento. Ahora se le pregunta al lattice directamente. Un vórtice de Taylor–Green en 2D funciona como regla graduada para $\nu$, porque su amplitud cae como $e^{-2\nu k^2 t}$.

import numpy as np
 
C = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]])
W = np.array([4/9] + [1/9]*4 + [1/36]*4)
CS2 = 1/3
 
def equilibrium_d2q9(rho, ux, uy):
    cu = C[:,0,None,None]*ux + C[:,1,None,None]*uy
    u2 = ux*ux + uy*uy
    return W[:,None,None]*rho*(1 + cu/CS2 + cu*cu/(2*CS2**2) - u2/(2*CS2))
 
def taylor_green_init(n, u0):
    k = 2*np.pi/n
    x, y = np.meshgrid(np.arange(n), np.arange(n), indexing='ij')
    ux = -u0*np.cos(k*x)*np.sin(k*y)
    uy =  u0*np.sin(k*x)*np.cos(k*y)
    return np.ones((n,n)), ux, uy, k
 
def advance_lattice(f, tau):
    rho = f.sum(axis=0)
    ux = (C[:,0,None,None]*f).sum(axis=0)/rho
    uy = (C[:,1,None,None]*f).sum(axis=0)/rho
    f = f - (f - equilibrium_d2q9(rho, ux, uy))/tau      # colision
    for i in range(9):                                    # adveccion
        f[i] = np.roll(f[i], (C[i,0], C[i,1]), axis=(0,1))
    return f, ux, uy
 
def measure_decay_rate(n, tau, u0=0.01, nstep=400, skip=100):
    rho, ux, uy, k = taylor_green_init(n, u0)
    f = equilibrium_d2q9(rho, ux, uy)
    amp = []
    for s in range(nstep):
        f, ux, uy = advance_lattice(f, tau)
        amp.append(np.sqrt((ux**2 + uy**2).mean()))
    t = np.arange(skip, nstep, dtype=float)
    a = np.log(np.array(amp[skip:]))                      # se descarta el transitorio inicial
    slope = np.polyfit(t, a, 1)[0]
    return -slope/(2*k*k), k
 
def viscosity_scan(n=32):
    print(f"{'tau':>6} {'nu_measured':>12} {'cs2*(tau-1/2)':>14} {'cs2*tau':>10} "
          f"{'err_half[%]':>11} {'err_naive[%]':>12}")
    for tau in (0.51, 0.55, 0.6, 0.8, 1.0, 1.5):
        nu, k = measure_decay_rate(n, tau)
        half = CS2*(tau - 0.5)
        naive = CS2*tau
        print(f"{tau:6.2f} {nu:12.6f} {half:14.6f} {naive:10.6f} "
              f"{100*abs(nu-half)/half:11.2f} {100*abs(nu-naive)/naive:12.2f}")
 
viscosity_scan()

   tau  nu_measured  cs2*(tau-1/2)    cs2*tau err_half[%] err_naive[%]
  0.51     0.003333       0.003333   0.170000        0.01        98.04
  0.55     0.016667       0.016667   0.183333        0.00        90.91
  0.60     0.033337       0.033333   0.200000        0.01        83.33
  0.80     0.100115       0.100000   0.266667        0.12        62.46
  1.00     0.167205       0.166667   0.333333        0.32        49.84
  1.50     0.337688       0.333333   0.500000        1.31        32.46

Los controles de $\theta$ y de $z$ recorren la misma regla que dejó la columna $\tau - 1/2$ en la tabla.

La región azul es el área bajo $e^{-s/\lambda}$. Las rayas rosas son la altura única con la que la regla $\theta$ la sustituye. Al salir de $\theta = 1/2$ la pendiente del error pasa de 3 a 2, y el $1.31\%$ de $\tau = 1.5$ es el resto que esa pendiente ya no cancela.

La columna de $\tau - 1/2$ se desvía entre 0.01 y 1.31 por ciento, el nivel de error de alto orden que deja la discretización del lattice. La columna de $\tau$ se desvía entre 32 y 98 por ciento. Esa magnitud ya no se puede llamar error.

La amplitud inicial se fija en 0.01 para borrar la compresibilidad. Con un número de Mach de orden 0.02, el término de tercer orden truncado de $f^{\rm eq}$ no se filtra en la tasa de decaimiento. Descartar los primeros 100 pasos cumple el mismo propósito. Al inicializar con la distribución de equilibrio, $f^{(1)}$ arranca en cero, de modo que la parte de no equilibrio necesita unas decenas de pasos para asentarse. Si se ajusta la pendiente incluyendo esa ventana, la viscosidad sale más pequeña de lo real.

También vale leer cómo el error del lado $\tau - 1/2$ crece despacio con $\tau$. En $\tau = 1.5$ llega a 1.31 por ciento, porque los términos $O(\Delta t^2)$ descartados en la expansión crecen con $\tau$. La corrección de medio paso no es la última palabra: es la historia hasta segundo orden.

Los cocientes $\tau/(\tau - 0.5)$ de la tabla no son un adorno. En $\tau = 0.51$ sale $0.51/0.01 = 51$, y $0.170000/0.003333$ repite ese 51. En $\tau = 0.55$ es $11$, en $0.60$ es $6$, en $1.00$ es $2$ y en $1.50$ es $1.5$. El error ingenuo de $\tau = 0.55$, $90.91\%$, es $1 - 0.016667/0.183333 = 1 - 1/11$. En $\tau = 0.80$ la medida $0.100115$ se aparta de $0.100000$ en $0.000115/0.100000 = 0.115\%$, que la tabla redondea a $0.12\%$. En $\tau = 1.50$, $(0.337688 - 0.333333)/0.333333$ da el $1.31\%$. La amplitud $0.01$, los $100$ pasos descartados y la malla $32$ son los mismos en toda la columna. Los errores ingenuos que quedan en la última columna son $98.04$, $90.91$, $83.33$, $62.46$, $49.84$ y $32.46$.

La frontera y el término de fuerza vuelven a facturar el mismo medio paso#

Este medio paso no aparece una sola vez.

Es la razón por la que bounce-back coloca la pared a mitad de camino entre dos nodos y no sobre un nodo: la reflexión ocurre a mitad de paso. Cuando $\tau$ se aleja de 0.5, esa posición efectiva de pared empieza a derivar con $\tau$, y eso es la conocida pérdida de precisión de segundo orden.

Es también la razón por la que el término de fuerza tipo Guo lleva un factor $\left(1 - \frac{1}{2\tau}\right)$. La fuerza se integra sobre la misma característica con la misma regla del trapecio, y como la transformación a $\bar f$ actúa solo sobre $\Omega_i$, queda un coeficiente suelto. Si falta ese paréntesis, un fluido que debería estar quieto se desplaza lentamente, o un perfil hidrostático bajo gravedad se inclina un poco.

Es la razón por la que reducir a la mitad el espaciado del lattice no significa duplicar $\tau$. Mantener $\nu$ fijo exige resolver $\tau_{\rm fine} - 1/2 = 2(\tau_{\rm coarse} - 1/2)$.

Cuando $\tau$ varía nodo a nodo, este medio paso es el piso. Como $\nu \to 0$ cuando $\tau \to 0.5$, cualquier problema donde la viscosidad caiga localmente hacia cero — fluidos no newtonianos, por ejemplo — choca con un muro construido en el propio lattice. Delante de ese muro el nodo pide un $\tau$ que la colisión ya no puede pagar.

0.5 no es una cota inferior, es el origen de coordenadas#

Memorizar $\tau > 0.5$ como condición de estabilidad es quedarse con la consecuencia. La causa es que el tiempo de relajación físico vale $\tau - 1/2$, y un tiempo de relajación negativo no es física. En $\tau = 0.51$ ese hueco es $0.01$, y $0.51/0.01 = 51$.

Lo que entra al código es $\nu$, no $\tau$. Se fija $\tau = 1/2 + \nu/c_s^2$ en un solo punto y después no se toca $\tau$ a mano. El factor de Guo $\left(1 - \frac{1}{2\tau}\right)$ obedece la misma cuenta, y el detalle de esa línea está en la línea que debe respetar un término de fuerza. Si el decaimiento se frena por un factor de 51, el código está usando $\tau$ donde corresponde $\nu$.