Metí la presión tal cual y una interfaz en reposo empezó a temblar — dos formas del forcing en LBM no ideal

Una gota en reposo empezó a fluir sola#

Puse la ecuación de estado de van der Waals sobre lattice Boltzmann (LBM, Lattice Boltzmann Method) y resolví un fluido bifásico. La condición inicial era una sola gota en reposo. El campo de densidad es liso y la velocidad es cero en todas partes.

Tras unos pasos, brotaron pequeñas velocidades cerca de la interfaz. Nadie la empujó y aun así el fluido fluye. A esta velocidad artificial se le llama corriente parásita (un flujo fantasma que aparece en la interfaz sin causa física).

Al acotar la causa, llegué a un punto del código. Todo se reduce a qué se mete en el término de fuerza. ¿Se mete la presión $p$ tal cual, o se mete el potencial químico $\mu$? Los libros de texto dicen que las dos son lo mismo. Este texto anota, como un libro de cuentas, hasta dónde es cierto ese "lo mismo". La respuesta es una línea: coinciden en el continuo y se separan en la red.

Al bajar temperature por debajo de 1, la isoterma se pliega. Una presión pasa a corresponder a tres densidades.

La línea horizontal es la presión de Maxwell, donde los dos lóbulos verdes tienen igual área. El hueco entre el punto azul y el punto rosa es el salto de densidad que el término de fuerza tendrá que sostener.

¿Por dónde entra la fuerza?#

LBM advecciona y colisiona funciones de distribución para recuperar las ecuaciones macroscópicas. En un gas ideal la presión que sale por sí sola es solo $c_s^2 \rho$ ($c_s$ es la velocidad del sonido en la red). La presión real de un fluido no ideal como van der Waals difiere de eso. Cubrir esa diferencia es la tarea del término de fuerza $\mathbf{F}$.

El objetivo es recuperar la siguiente ecuación de momento.

$p$ es la presión de van der Waals, $\boldsymbol{\tau}$ es el esfuerzo viscoso y el último término es el esfuerzo de Korteweg que levanta la interfaz, con $\kappa$ como su intensidad. Como el streaming aporta $c_s^2\rho$, la fuerza solo tiene que rellenar la diferencia restante.

Aquí surgen dos ramas. Una forma que usa la presión tal cual, y una forma que usa el potencial químico proveniente de la energía libre.

La primera es la forma de presión (el tema de hoy, "meter $p$ tal cual"), la segunda es la forma de energía libre. La temperatura a la que se pliega la isoterma ya se bajó en el dial de arriba. Queda anotar por qué las dos fuerzas coinciden en el continuo.

Gibbs-Duhem: la presión y el potencial químico dicen lo mismo dos veces#

Si las dos formas son iguales lo decide una sola relación. Es la relación de Gibbs–Duhem que enlaza presión y potencial químico a temperatura constante.

Sustituye esto en $\mathbf{F}_\mu$. Como $-\rho\nabla\mu = -\nabla p$, el término del potencial químico se convierte directamente en un gradiente de presión. El $c_s^2\nabla\rho$ restante engrana exactamente con el término que la forma de presión llevaba como $-\nabla(-c_s^2\rho)$. Al final $\mathbf{F}_p = \mathbf{F}_\mu$.

Es decir, el único fundamento de la afirmación "se puede meter $p$ tal cual" es Gibbs–Duhem. Van der Waals satisface esta relación de forma exacta, porque la ecuación de estado y la energía libre salen de la misma termodinámica. Es el mismo tipo de equivalencia que los esquemas de forcing de LBM confirmados por la expansión de Chapman–Enskog: cada uno con una cara distinta, y aun así todos recuperan la misma ecuación macroscópica.

Densidades de coexistencia y Gibbs-Duhem, verificadas en Python#

Las palabras solas no son fiables. Escribamos la ecuación de estado de van der Waals en unidades reducidas, hallemos las densidades de coexistencia por temperatura con el método de Newton, y luego midamos directamente el defecto de Gibbs–Duhem $\lvert \mathrm{d}p/\mathrm{d}\rho - \rho\,\mathrm{d}\mu/\mathrm{d}\rho \rvert$.

import numpy as np
 
A, B, R = 9.0 / 8.0, 1.0 / 3.0, 1.0   # van der Waals, unidades reducidas (rho_c=1, T_c=1)
 
def p_eos(rho, T):                     # presión de van der Waals
    return rho * R * T / (1.0 - B * rho) - A * rho * rho
 
def mu_eos(rho, T):                    # potencial químico mu = df/drho
    return R * T * (np.log(rho / (1.0 - B * rho)) + B * rho / (1.0 - B * rho)) - 2.0 * A * rho
 
def maxwell(T):                        # regla de igual área: densidades donde p y mu coinciden en las dos fases
    x = np.array([0.30, 1.90]); h = 1e-8
    for _ in range(80):
        f = np.array([p_eos(x[0], T) - p_eos(x[1], T), mu_eos(x[0], T) - mu_eos(x[1], T)])
        J = np.empty((2, 2))
        for k in range(2):
            y = x.copy(); y[k] += h
            g = np.array([p_eos(y[0], T) - p_eos(y[1], T), mu_eos(y[0], T) - mu_eos(y[1], T)])
            J[:, k] = (g - f) / h
        x -= np.linalg.solve(J, f)
    return x[0], x[1]
 
# Gibbs-Duhem: dp = rho d(mu). El único fundamento para meter la presión tal cual.
print("T/Tc  rho_vap  rho_liq |  max| dp/drho - rho*dmu/drho |")
for T in (0.95, 0.90, 0.85):
    rv, rl = maxwell(T)
    r = np.linspace(rv, rl, 400)
    dp = np.gradient(p_eos(r, T), r)
    dmu = np.gradient(mu_eos(r, T), r)
    err = np.abs(dp - r * dmu).max()
    print(f"{T:4.2f}  {rv:7.4f}  {rl:7.4f} |  {err:.2e}")

T/Tc  rho_vap  rho_liq |  max| dp/drho - rho*dmu/drho |
0.95   0.5790   1.4617 |  2.94e-04
0.90   0.4257   1.6573 |  9.44e-04
0.85   0.3197   1.8071 |  1.98e-03

El defecto está en el nivel de $10^{-3}$, y hasta eso se debe a que las derivadas se midieron con diferencias finitas. Cuanto más baja la temperatura, más se ensancha el hueco entre las densidades de coexistencia. Gibbs–Duhem se cumple de forma exacta en el continuo. Hasta aquí la forma de presión y la forma de energía libre son completamente idénticas.

Conviene restar las columnas de coexistencia. En $T/T_c=0.95$ el vapor es $0.5790$ y el líquido $1.4617$, un hueco de $0.8827$. En $0.90$ el hueco es $1.6573-0.4257=1.2316$. En $0.85$ es $1.8071-0.3197=1.4874$. Bajar la temperatura $0.10$ estira el salto de densidad de $0.8827$ a $1.4874$, cerca de $1.7$ veces. En el mismo tramo el defecto crece por $2.94\times10^{-4}$, $9.44\times10^{-4}$ y $1.98\times10^{-3}$. La identidad $\mathrm{d}p=\rho\,\mathrm{d}\mu$ no se rompió. Un np.gradient de 400 puntos solo gasta más error en una isoterma más empinada. En $T=0.90$ la razón líquido/vapor es $1.6573/0.4257\approx 3.89$. Esa razón es el salto de densidad que el forcing tiene que sostener.

Al superponer las dos líneas de fuerza, solo un lugar discrepa. $-\nabla(p-c_s^2\rho)$ en $\mathbf{F}_p$ y $-\rho\nabla\mu+c_s^2\nabla\rho$ en $\mathbf{F}_\mu$ coinciden cuando vale $\mathrm{d}p=\rho\,\mathrm{d}\mu$. El término de Korteweg $\kappa\rho\nabla(\nabla^2\rho)$ está escrito igual en ambos lados, así que no abre hueco. $c_s^2=1/3$ es la corrección de la velocidad del sonido en la red, ya anotada, no un empuje físico sobre la interfaz. En el libro continuo las dos formas coinciden, y esa frase termina aquí. En la red la frase se abre por el $1.60\times10^{-2}$ que imprime la tabla posterior.

En la red discreta las dos se separan#

El problema es la red. No hay garantía de que el $\mathrm{d}p = \rho\,\mathrm{d}\mu$ del continuo se cumpla bajo la derivación discreta. Un $\nabla p$ con diferencia centrada y $\rho\,\nabla\mu$ se desvían entre sí donde la densidad se dobla de forma brusca, como en una interfaz.

Resolví una interfaz plana en reposo con la condición de Euler–Lagrange, y luego calculé las dos formas de fuerza directamente sobre ella.

import numpy as np
 
A, B, R = 9.0 / 8.0, 1.0 / 3.0, 1.0
CS2, KAPPA, T = 1.0 / 3.0, 0.02, 0.90
 
def p_eos(rho):  return rho * R * T / (1.0 - B * rho) - A * rho * rho
def mu_eos(rho): return R * T * (np.log(rho / (1.0 - B * rho)) + B * rho / (1.0 - B * rho)) - 2.0 * A * rho
def dmu(rho):    return R * T * (1.0 / (rho * (1.0 - B * rho)) + B / (1.0 - B * rho) ** 2) - 2.0 * A
def diff1(a, dx): return (np.roll(a, -1) - np.roll(a, 1)) / (2.0 * dx)
def lap(a, dx):   return (np.roll(a, -1) - 2.0 * a + np.roll(a, 1)) / dx ** 2
 
def maxwell():
    x = np.array([0.30, 1.90]); h = 1e-8
    for _ in range(80):
        f = np.array([p_eos(x[0]) - p_eos(x[1]), mu_eos(x[0]) - mu_eos(x[1])])
        J = np.empty((2, 2))
        for k in range(2):
            y = x.copy(); y[k] += h
            g = np.array([p_eos(y[0]) - p_eos(y[1]), mu_eos(y[0]) - mu_eos(y[1])])
            J[:, k] = (g - f) / h
        x -= np.linalg.solve(J, f)
    return x[0], x[1]
 
rv, rl = maxwell()
mu_co = mu_eos(np.array([rv]))[0]
 
# (1) Resuelve la interfaz plana en reposo y compara las dos formas de forcing
NX = 240
xs = np.arange(NX)
rho = 0.5 * (rl + rv) + 0.5 * (rl - rv) * (np.tanh((xs - NX / 4) / 6.0) - np.tanh((xs - 3 * NX / 4) / 6.0) - 1.0)
for _ in range(6000):                                    # relaja el residuo de Euler-Lagrange
    rho -= 0.15 * (mu_eos(rho) - KAPPA * lap(rho, 1.0) - mu_co) / dmu(rho)
 
Gp  = -diff1(p_eos(rho) - CS2 * rho, 1.0) + KAPPA * rho * diff1(lap(rho, 1.0), 1.0) - diff1(CS2 * rho, 1.0)
Gmu = -rho * diff1(mu_eos(rho) - KAPPA * lap(rho, 1.0), 1.0) + CS2 * diff1(rho, 1.0) - diff1(CS2 * rho, 1.0)
gd  = np.abs(diff1(p_eos(rho), 1.0) - rho * diff1(mu_eos(rho), 1.0)).max()
 
print(f"coexistence         rho_vap = {rv:.4f}   rho_liq = {rl:.4f}")
print(f"free-energy form    max|G_mu|        = {np.abs(Gmu).max():.2e}   (well-balanced)")
print(f"pressure form       max|G_p|         = {np.abs(Gp).max():.2e}   (spurious force)")
print(f"gap between forms   max|G_p - G_mu|  = {np.abs(Gp - Gmu).max():.2e}")
print(f"discrete Gibbs-Duhem defect          = {gd:.2e}   <- the gap, exactly")
 
# (2) Refina dx sobre la misma interfaz física y el defecto converge a 0 rápidamente
print("\ncells/interface |  Gibbs-Duhem defect  order")
prev = None
for n in (10, 20, 40, 80):
    L = 40.0; N = int(L * n / 10)
    z = np.linspace(-L / 2, L / 2, N, endpoint=False); dx = z[1] - z[0]
    r = 0.5 * (rl + rv) - 0.5 * (rl - rv) * np.tanh(z / (0.1 * n))
    d = np.abs(diff1(p_eos(r), dx) - r * diff1(mu_eos(r), dx))[N // 4:3 * N // 4].max()
    order = "" if prev is None else f"{np.log(prev / d) / np.log(2.0):5.2f}"
    print(f"{n:9d}       |  {d:.3e}          {order}")
    prev = d

coexistence         rho_vap = 0.4257   rho_liq = 1.6573
free-energy form    max|G_mu|        = 2.02e-16   (well-balanced)
pressure form       max|G_p|         = 1.60e-02   (spurious force)
gap between forms   max|G_p - G_mu|  = 1.60e-02
discrete Gibbs-Duhem defect          = 1.60e-02   <- the gap, exactly
 
cells/interface |  Gibbs-Duhem defect  order
       10       |  1.140e-02          
       20       |  1.428e-03           3.00
       40       |  5.115e-05           4.80
       80       |  1.611e-06           4.99

Al subir resolution, la interfaz abarca más celdas y la altura del pico rojo cae con ellas.

El pico rojo es la fuerza extra que la forma de presión apila sobre la de energía libre, $G_p - G_\mu$. La línea verde discontinua es cero. Con más celdas, el pico se sienta. La fuerza fantasma es un resto de esa diferencia, no de la física.

La forma de energía libre queda en $2.02\times10^{-16}$ sobre la interfaz. La de presión deja $1.60\times10^{-2}$. La diferencia coincide con el defecto discreto de Gibbs–Duhem hasta el decimal impreso. $1.60\times10^{-2}/2.02\times10^{-16}\approx 7.9\times10^{13}$. El redondeo de máquina y la fuerza fantasma comparten una línea de la tabla.

La tabla de malla es el mismo libro. Diez celdas por interfaz dan un defecto de $1.140\times10^{-2}$; veinte dan $1.428\times10^{-3}$, orden observado $3.00$. Cuarenta dan $5.115\times10^{-5}$ (orden $4.80$) y ochenta $1.611\times10^{-6}$ (orden $4.99$). Pasar de 10 a 80 celdas, un factor ocho, recorta el defecto en $1.140\times10^{-2}/1.611\times10^{-6}\approx 7.1\times10^{3}$. Lo que la diferencia centrada deja escapar de $\mathrm{d}p=\rho\,\mathrm{d}\mu$ en la pendiente de la interfaz es la brecha entre las dos fuerzas. Esa fuerza empuja la interfaz en reposo.

Dejar la forma de presión en una interfaz delgada pasa la factura#

La forma de energía libre ($-\rho\nabla\mu$) queda balanceada en la interfaz por definición. Incluso en una red gruesa la fuerza fantasma se pega a $2.02\times10^{-16}$. Si la presión se queda, $\nabla p$ no se diferencia a gusto: tiene que ser la misma diferencia que $\rho\,\nabla\mu$. Solo así vive Gibbs–Duhem discreto y baja el $1.60\times10^{-2}$ de la tabla.

Con la interfaz resuelta en 4 o 5 celdas o más, el defecto de la forma de presión cae hasta orden $4.99$. En la práctica la interfaz anda por 3 celdas. Ahí la fuerza fantasma sigue el cuadrado del salto de densidad. Un término de fuerza que cumple la identidad continua y deja vacíos los momentos discretos conserva el mismo residuo; ese libro está en el límite del tercer momento del término de fuerza. La comodidad de "usar $p$ tal cual" se cobra en espesor de interfaz.