Darle su propio τ a la conducción dejó el calentamiento viscoso atrapado en Pr = 1 — el cuasi-equilibrio g* en LB de dos poblaciones

Pr cambió, ¿por qué la temperatura no?#

Considérese un canal de Couette. La pared inferior está en reposo y la superior se desliza a velocidad $U$. Ambas paredes tienen la misma temperatura. Aun así, el centro se calienta. La viscosidad convierte energía cinética en calor (calentamiento viscoso). El aumento de temperatura en el centro es exactamente proporcional al número de Prandtl (Pr, cociente entre difusión de cantidad de movimiento y difusión de calor).

Este artículo trata de cómo cambiar Pr en un modelo térmico de Boltzmann en red (LB, Lattice Boltzmann). Si la distribución se divide en dos conjuntos con tiempos de relajación propios, Pr parece fácil de ajustar. Al ejecutarlo, la conducción cambia, pero el calentamiento viscoso no se mueve de su valor de Pr = 1. Aquí se explica por qué y se presenta la solución de una línea: el cuasi-equilibrio $g^*$.

Al mover el control de Pr cambia la separación entre el reloj del esfuerzo y el de la conducción, y también a qué curva se pega el flujo de trabajo viscoso.

Con Pr = 0.71, la curva naranja (esfuerzo) y la verde discontinua (flujo de trabajo viscoso en el modelo g*) permanecen superpuestas. La roja (modelo simple de dos τ) sigue el reloj morado de la conducción y se separa del esfuerzo. Al llevar Pr a 1, las cuatro curvas se reducen a dos.

Un solo conjunto de poblaciones solo conoce Pr = 1#

El operador de colisión BGK lleva todos los momentos fuera de equilibrio hacia el equilibrio a una misma tasa $\omega$. Esfuerzo y flujo de calor se relajan con el mismo reloj. Por eso viscosidad y conductividad quedan atadas al mismo $\omega$. De ahí sale Pr = 1 en la ecuación BGK continua.

El aire tiene Pr = 0.71 y el agua cerca de 7. Un modelo fijado en Pr = 1 ni siquiera reproduce la razón de espesores de capa límite. Hay dos remedios. Uno relaja solo el momento del flujo de calor a otra tasa dentro del operador de colisión (MRT, tipo Shakhov). El otro entrega la energía a una distribución separada. Este artículo sigue el segundo camino.

Dos conjuntos: f lleva masa y cantidad de movimiento, g la energía total#

En un modelo de dos poblaciones, $f_i$ transporta densidad y cantidad de movimiento, y $g_i$ transporta la energía total $E=\rho(c_vT+u^2/2)$ ($i=0,\dots,Q-1$). Las ecuaciones cinéticas son:

$\omega_\mu$ fija la viscosidad (viscosidad cinemática $\nu$) y $\omega_\kappa$ la conducción (difusividad térmica $\alpha$). $g_i^{eq}$ es el equilibrio local y $g_i^*$ el cuasi-equilibrio.

Las dos tasas se relacionan con las difusividades así:

$c_s^2=1/3$ es el cuadrado de la velocidad del sonido de la red. Con $\tau_\mu=1/\omega_\mu$ y $\tau_\kappa=1/\omega_\kappa$, $\mathrm{Pr}=(\tau_\mu-\tfrac12)/(\tau_\kappa-\tfrac12)$.

El código fija $\tau_\mu=1$. Entonces $\tau_\mu-\tfrac12=\tfrac12$ y $\tau_\kappa=\tfrac12+(\tfrac12)/\mathrm{Pr}$. Las cuatro entradas de la tabla salen de ahí. Pr = 0.5 da $\tau_\kappa=1.5$. Pr = 0.71 suma $0.5/0.71\approx0.704$, así que $\tau_\kappa\approx1.204$. Pr = 1 da $\tau_\kappa=1$. Pr = 2 da $\tau_\kappa=0.75$. Como $\omega=1/\tau$, el calor a Pr = 2 se relaja con $\omega_\kappa\approx1.333$ y la viscosidad queda en $\omega_\mu=1$. El trabajo viscoso pertenece al reloj lento. El modelo simple de dos τ lo monta en el $\omega_\kappa$ rápido. Solo la fila Pr = 1 tiene las dos tasas iguales, y ambas columnas marcan 1.000. Si se mira solo esa fila, el error permanece callado.

El trabajo viscoso se esconde en el flujo de energía#

En este modelo de juguete, $g^{eq}$ fija solo tres momentos. Orden cero: $E$. Primer orden (flujo de energía): $(E+p)\mathbf{u}$. Segundo orden: $c_s^2E\,\mathbf{I}+E\,\mathbf{u}\mathbf{u}$. La parte fuera de equilibrio del primer momento, $\mathbf{Q}$, es el flujo de calor real.

En el primer orden de la expansión de Chapman–Enskog (que resuelve el no equilibrio para número de Knudsen pequeño), $\mathbf{Q}\approx-\tau\,\nabla\!\cdot\!\mathbf{R}^{eq}$. Aquí $\nabla(c_s^2E)$ se divide en dos partes:

El primer término es la conducción de Fourier. El segundo es la materia prima del trabajo viscoso $-\mathbf{u}\cdot\boldsymbol{\sigma}$. Tiene exactamente la forma del esfuerzo fuera de equilibrio de $f$, $\boldsymbol{\Pi}^{neq}=\mathbf{P}-\mathbf{P}^{eq}\approx-\tau_\mu\rho c_s^2(\nabla\mathbf{u}+\nabla\mathbf{u}^T)$. El problema es el $\tau$ que lo acompaña.

Si $g$ solo se arrastra hacia el equilibrio con $\omega_\kappa$, ambas partes salen con $\tau_\kappa$. El coeficiente del trabajo viscoso escala con $\alpha$ en lugar de $\nu$. La conducción es correcta, pero el calentamiento viscoso se queda en su valor de Pr = 1. Es la curva roja de la primera viz.

Cuasi-equilibrio g*: atar solo el flujo al reloj del esfuerzo#

La solución es una línea en la definición de $g^*$. Comparte todos los momentos con $g^{eq}$ salvo el flujo de energía:

$\mathbf{P}$ es el flujo de cantidad de movimiento real medido en $f$ y $\mathbf{P}^{eq}$ su valor de equilibrio.

Al tomar el primer momento de la ecuación cinética se obtiene la actualización de un paso de los dos momentos fuera de equilibrio:

$\mathbf{Q}$ se relaja con $\omega_\kappa$, así que la conducción sigue a $\alpha$. Al mismo tiempo, si $\mathbf{Q}=\mathbf{u}\cdot\boldsymbol{\Pi}^{neq}$, en el paso siguiente se mantiene $\mathbf{Q}'=\mathbf{u}\cdot\boldsymbol{\Pi}'^{neq}$. La parte de trabajo viscoso viaja con el reloj del esfuerzo $\omega_\mu$. El cálculo de Chapman–Enskog dice lo mismo: el trabajo viscoso con $\tau_\kappa$ se cancela con el término $(\omega_\kappa-\omega_\mu)$ y queda solo $\tau_\mu$.

Con $\omega_\kappa=\omega_\mu$ el término $g^*$ desaparece. El modelo se reduce a relajación simple y Pr = 1. $g^*$ es una corrección que solo actúa cuando Pr $\ne 1$. Una red térmica de una sola población dejaba Pr pegado en 1 aunque τ subiera. Es el mismo muro que las dos constantes que quedaron trabadas. Partir las poblaciones no sale de ese muro si el reloj del trabajo viscoso no se mueve con ellas. La pregunta sigue siendo qué momento se relaja a qué tasa.

Código: el mismo canal de Couette, dos tipos de g#

La solución analítica es sencilla. La ecuación de energía estacionaria es $kT''+\mu(u')^2=0$. Con un perfil de velocidad lineal, la temperatura es una parábola:

El máximo en el centro es $\theta_{max}=\mathrm{Pr}$. El código siguiente usa D2Q9 con una sola celda periódica en $x$ y resuelve solo a lo largo de $y$. Las paredes usan la extrapolación fuera de equilibrio de Guo. Para simplificar, $f$ usa un equilibrio isotérmico y la energía se acopla en un solo sentido. En un Couette de bajo Mach ($U=0.1$) la presión es uniforme, de modo que esta simplificación no altera el resultado.

import numpy as np
 
# D2Q9, una celda periódica en x → arreglos 1D en y
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w = np.array([4/9] + [1/9]*4 + [1/36]*4)
CV, T0 = 1.0, 1.0
 
def feq_lattice(rho, ux):
    cu = cx[:, None]*ux
    return w[:, None]*rho*(1 + 3*cu + 4.5*cu**2 - 1.5*ux**2)
 
def geq_energy(E, rho, ux):
    p = rho*(E/rho - 0.5*ux**2)/CV          # p = rho*T
    qx = (E + p)*ux
    cu = cx[:, None]*ux
    return w[:, None]*(E + 3*cx[:, None]*qx + 4.5*E*cu**2 - 1.5*E*ux**2)
 
def couette_heating(Pr, quasi, N=21, U=0.1, tau_mu=1.0, steps=8000):
    om_mu = 1/tau_mu
    om_ka = 1/(0.5 + (tau_mu - 0.5)/Pr)
    rho = np.ones(N); ux = np.zeros(N)
    E = rho*CV*T0*np.ones(N)
    f = feq_lattice(rho, ux); g = geq_energy(E, rho, ux)
    uw = np.zeros(N); uw[-1] = U
    for _ in range(steps):
        rho = f.sum(0); ux = (cx[:, None]*f).sum(0)/rho; E = g.sum(0)
        fe = feq_lattice(rho, ux); ge = geq_energy(E, rho, ux)
        Pxy = ((cx*cy)[:, None]*(f - fe)).sum(0)
        f = f + om_mu*(fe - f)
        if quasi:   # g*: cuasi-equilibrio con u·(P - P_eq) sumado al flujo de energía
            gs = ge + w[:, None]*3*cy[:, None]*(ux*Pxy)
            g = g + om_mu*(ge - g) + (om_ka - om_mu)*(gs - g)
        else:       # dos τ simple: g solo se relaja al equilibrio con ω_κ
            g = g + om_ka*(ge - g)
        for k in range(9):
            f[k] = np.roll(f[k], cy[k]); g[k] = np.roll(g[k], cy[k])
        for wall, nb in ((0, 1), (N-1, N-2)):   # paredes con extrapolación fuera de equilibrio
            r = f[:, nb].sum(); un = (cx*f[:, nb]).sum()/r; En = g[:, nb].sum()
            fn = feq_lattice(np.array([r]), np.array([un]))[:, 0]
            gn = geq_energy(np.array([En]), np.array([r]), np.array([un]))[:, 0]
            Ew = r*(CV*T0 + 0.5*uw[wall]**2)
            f[:, wall] = feq_lattice(np.array([r]), np.array([uw[wall]]))[:, 0] + f[:, nb] - fn
            g[:, wall] = geq_energy(np.array([Ew]), np.array([r]), np.array([uw[wall]]))[:, 0] + g[:, nb] - gn
    rho = f.sum(0); ux = (cx[:, None]*f).sum(0)/rho
    T = (g.sum(0)/rho - 0.5*ux**2)/CV
    theta = (T - T0)/(U**2/(8*CV))
    return theta.max(), ux
 
print(" Pr   exact   2-tau(g_eq)   g* quasi-eq")
for Pr in (0.5, 0.71, 1.0, 2.0):
    a, _ = couette_heating(Pr, False)
    b, u = couette_heating(Pr, True)
    print(f"{Pr:4.2f}  {Pr:6.3f}   {a:9.3f}    {b:9.3f}")
print(f"u(mid)/U = {u[10]/0.1:.4f}")

 Pr   exact   2-tau(g_eq)   g* quasi-eq
0.50   0.500       1.000        0.498
0.71   0.710       1.000        0.710
1.00   1.000       1.000        1.000
2.00   2.000       1.000        1.997
u(mid)/U = 0.5000

Ahora el mismo código corre en el navegador. Ambos modelos avanzan lado a lado en el mismo canal.

Al mover el control de Pr, el pico verde sigue el valor analítico discontinuo Pr, mientras el rojo queda pegado en 1.000. Cambiar $U$ o $\tau_\mu$ tampoco aparta del 1 el $\theta_{max}$ rojo. Conviene probar Pr = 3, que no aparece en el bloque de salida.

Lo que dice el 1.000#

La columna del modelo simple de dos τ marca 1.000 sea cual sea Pr. No es ruido, es estructura. El coeficiente del trabajo viscoso escala con $\tau_\kappa-\tfrac12$, y la conductividad también. En su cociente, $\theta_{max}$, $\tau_\kappa$ se cancela exactamente. Para aire (Pr = 0.71) se sobrestima el aumento de temperatura entre las paredes en 41 %. Para un fluido con Pr = 7, como el agua, se subestima siete veces.

Con $U=0.1$ y $c_v=1$, el denominador es $U^2/(8c_v)=0.00125$. El $\Delta T$ real del aire es $0.71\times0.00125\approx8.9\times10^{-4}$. El modelo simple de dos τ da $1.25\times10^{-3}$. En la red eso es 1.0009 frente a 1.0013. Con dos decimales, ambos salen 1.00. A Pr = 7 el aumento real es $0.00875$ y el modelo da $0.00125$. 1.009 frente a 1.001. Solo $\theta$ separa 1 de 7.

La columna de $g^*$ marca 0.498, 0.710, 1.000 y 1.997. 0.498 queda 0.4 % por debajo de 0.500, y 1.997 queda 0.15 % por debajo de 2. 0.710 coincide con el exacto. En 21 nodos, el resto es la pared y la malla. La velocidad central $u/U=0.5000$ es idéntica porque ambos modelos comparten el mismo $f$. El error de temperatura no aparece en el campo de velocidad.

En compresible, f también recibe un cuasi-equilibrio#

El $f$ de este artículo usó un equilibrio isotérmico. En los modelos compresibles de dos poblaciones, $f$ también incorpora un cuasi-equilibrio $f^*$. Los terceros momentos de D2Q9 tienen un defecto: $c_{ix}^3=c_{ix}$. Cuando la temperatura varía o el Mach crece, ese defecto introduce un término de error en el esfuerzo. $f^*$ lleva un tensor de presión que lo corrige. $f^*$ repara el flujo de cantidad de movimiento; $g^*$ añade el trabajo viscoso al flujo de energía.

El principio es el mismo. Equilibrio y cuasi-equilibrio comparten las cantidades conservadas (densidad, cantidad de movimiento, energía total). Solo difiere el momento que necesita su propia tasa de relajación. Dónde poner el término de presión de un fluido no ideal es, al final, el sitio de ese momento.

La línea que hay que revisar cuando un código LB térmico fija τκ aparte#

Hay que abrir el paso de colisión de la distribución de energía y ver dónde multiplica $\omega_\kappa$. Si todo $g$ se arrastra hacia el equilibrio con $\omega_\kappa$, el calentamiento viscoso queda en su valor de Pr = 1. Basta una corrida de Couette para comprobarlo. Si $\theta_{max}$ no sigue a Pr, falta el término de cuasi-equilibrio.

Validar solo con $\omega_\kappa=\omega_\mu$ deja el error en silencio. Una comprobación térmica de LB necesita un caso de calentamiento viscoso con Pr $\ne 1$. La temperatura que se difunde en un fluido en reposo no pone en marcha este término. La analogía de Reynolds entre fricción y calor parte de que ese cociente queda atado a Pr. Si $\theta_{max}$ no sigue a Pr, esa premisa ya está rota en la red.