Devolver la nitidez a la interfaz en cada paso — Implementando el afilado anti-difusión (IST)

Ejecuté un caso de choque–burbuja con un modelo de cinco ecuaciones. Tras 200 pasos, el borde de la burbuja de helio se había hinchado de 6 celdas a 20. Las cantidades físicas eran correctas, pero la interfaz quedó hecha una papilla. Ese es el destino del método de interfaz difusa (un enfoque que deja a propósito que la interfaz se extienda sobre un ancho finito). Nguyen et al. (2021) deshacen ese emborronamiento con una sola pasada de posprocesado después de cada paso temporal. Hoy construimos esa técnica anti-difusión desde cero y comprobamos, con nuestros propios ojos, por qué el grosor de la interfaz se bloquea en una constante.

Con numerical diffusion activo, el borde de 6 celdas se come celdas; con sharpening (IST) activo, esa extensión se detiene cerca de ε/a.

La banda amarilla es el grosor objetivo. Subir compression a la adelgaza, subir regularization ε la engrosa otra vez, y el recuento de celdas sigue ese cociente.

Información del artículo#

• Título: Numerical modeling of multiphase compressible flows with the presence of shock waves using an interface-sharpening five-equation model
• Autores: Van-Tu Nguyen, Thanh-Hoang Phan, Warn-Gyu Park (Pusan National University)
• Fuente: International Journal of Multiphase Flow, 2021, 103542
• DOI: 10.1016/j.ijmultiphaseflow.2020.103542
• Resumen en una línea: Acoplar el afilado de interfaz (IST) a un modelo compresible bifásico de cinco ecuaciones como paso de posprocesado, para mantener las interfaces a un grosor constante incluso junto a los choques.

Por qué la interfaz engorda por sí sola#

El modelo de cinco ecuaciones es un sistema casi conservativo: leyes de conservación más una función de color. La fracción de volumen $\alpha_1$ (la porción de fluido 1 dentro de una celda) obedece una ecuación de advección como esta.

Aquí $\alpha_1$ es la fracción de volumen del fluido 1, $\mathbf{u}$ es la velocidad de mezcla y el término de Kapila del lado derecho corrige la diferencia de compresibilidad entre las dos fases.

El problema empieza al resolver esto con un esquema de captura de choques. En cuanto se interpolan suavemente los valores en las caras de la celda, el escalón de $\alpha_1$ que debería seguir afilado se emborrona un poco en cada paso. Esa difusión numérica se acumula. Tras muchos pasos, la interfaz pasa de unas pocas celdas a decenas, y rasgos clave —la forma de la burbuja, las posiciones de reflexión de los choques— desaparecen por completo.

Los esquemas de alto orden como WENO reducen la difusión, pero introducen oscilaciones y son costosos en varias dimensiones. Por eso los autores giran de lado. En vez de tocar el esquema en sí, arreglan el campo de $\alpha$ por separado después de cada paso.

Anti-difusión: la ecuación de afilado#

El núcleo es la ecuación de regularización de Shukla et al. (2010). Hacemos evolucionar $\alpha$ durante unas pocas iteraciones en un pseudo-tiempo $\tau$, no en tiempo físico.

$\varepsilon$ es el coeficiente de difusión que regula el grosor, $a$ es la intensidad de compresión que aprieta la interfaz y $\hat{\mathbf{n}} = \nabla\alpha / |\nabla\alpha|$ es la dirección normal a la interfaz.

Los dos términos juegan al tira y afloja. El primero es difusión ordinaria, así que dispersa $\alpha$. El segundo, con signo opuesto, es anti-difusión: comprime la interfaz a lo largo de su normal. Gracias al factor $\alpha(1-\alpha)$, la compresión actúa solo en la interfaz ($0<\alpha<1$) y se apaga dentro de los fluidos puros ($\alpha=0$ o $1$).

Resolviendo el estado estacionario ($\partial_\tau \alpha = 0$) en una dimensión, los dos flujos se equilibran.

La interfaz converge a una forma de tangente hiperbólica, y su grosor escala como $\varepsilon/a$. En otras palabras, el grosor queda fijado de antemano. Esta es la propiedad de "grosor constante" que los autores subrayan. Cuando la interfaz mantiene siempre el mismo número de celdas, los choques y las discontinuidades de contacto pueden capturarse de forma estable.

En 1D: el tira y afloja entre emborronado y afilado#

El mismo tira y afloja, escrito, es el script de abajo. Imita el emborronado que deja el transporte y luego corre la ecuación de afilado en pseudo-tiempo.

import numpy as np
 
N, dx = 201, 1.0 / 200
x = np.linspace(0, 1, N)
 
def band_width(a, lo=0.05, hi=0.95):
    # número de celdas contadas como interfaz (medida del grosor)
    return int(np.sum((a > lo) & (a < hi)))
 
def smear_once(a, D=0.16):
    # el esquema de transporte engordando la interfaz en cada paso
    lap = np.zeros_like(a)
    lap[1:-1] = a[2:] - 2 * a[1:-1] + a[:-2]
    return a + D * lap
 
def sharpen_sweep(a, comp_a, eps, iters=200):
    # regularización anti-difusión en forma de flujo (Shukla 2010, ecs. (38)-(39))
    dtau = 0.9 * min(dx * dx / (2 * eps), dx / comp_a)
    for _ in range(iters):
        af = 0.5 * (a[:-1] + a[1:])          # alpha centrada en la cara
        dA = a[1:] - a[:-1]
        s = np.sign(dA)
        J = comp_a * af * (1 - af) * s - eps * dA / dx   # flujo en la cara
        a[1:-1] = np.clip(a[1:-1] - dtau / dx * (J[1:] - J[:-1]), 0, 1)
    return a
 
# empezar desde un escalón afilado -> dejar que se emborrone 20 veces
a = 0.5 * (1 + np.tanh((x - 0.5) / 0.012))
for _ in range(20):
    a = smear_once(a)
print("banda tras emborronar:", band_width(a), "celdas")     # -> más ancha
 
a = sharpen_sweep(a, comp_a=1.0, eps=0.004)
print("banda tras afilar:", band_width(a), "celdas")         # -> se encoge al objetivo

Tras el emborronado la banda queda ancha y, tras el afilado, cae a las pocas celdas que fija $\varepsilon/a$. El esquema no se toca. sharpen_sweep solo recibe el arreglo. No sabe qué solver de Riemann lo llenó.

Los valores por defecto son comp_a=1.0 y eps=0.004, así que la escala de grosor $\varepsilon/a$ vale 0.004. El factor de pendiente del perfil estacionario es $a/(2\varepsilon)=125$. La tanh con la que arranca el código tiene denominador 0.012, de modo que el objetivo es más empinado que el escalón inicial. $dx=1/200=0.005$. El tope del pseudo-tiempo es el menor entre $dx^2/(2\varepsilon)=0.000025/0.008=0.003125$ y $dx/a=0.005$. Gana el límite difusivo. Multiplicado por 0.9, $d\tau$ vale 0.0028125. Por encima de ese tope la anti-difusión cambia de signo y la interfaz oscila.

La pieza de compresión del flujo en la cara escala con $\alpha(1-\alpha)$. Ese producto llega a 0.25 en $\alpha=0.5$ y es cero en un fluido puro, así que allí la compresión se apaga. El contador de banda solo suma celdas entre 0.05 y 0.95. Doblar solo $\varepsilon$ dobla el grosor objetivo; doblar solo $a$ lo parte por dos. Los 200 pasos del arranque son tiempo físico. Pasar de 6 celdas a 20 es un factor de unos 3.3. D=0.16 dentro de smear_once e iters=200 dentro del afilado son capas distintas. Las 1 a 3 pasadas por paso del artículo suponen que el paso físico apenas difunde. El término de Kapila $\alpha_1 K\nabla\cdot\mathbf{u}$ empuja la fracción de volumen cuando un choque dispara la divergencia. Esa fuente no es difusión numérica, y el IST no la borra. A diferencia de la relajación de presión separada en el procedimiento de Saurel (2009), esta pasada solo devuelve el número de celdas.

Qué significa mantener el grosor constante#

Adelgazar la interfaz no es lo mismo que fijar su grosor a una constante. Aquí es exactamente donde los autores trazan una línea frente a enfoques anteriores (Tiwari et al. 2013). Añadir una función de afilado al término fuente reduce el error de difusión, pero el grosor vaga en el espacio y el tiempo. Fijar $\varepsilon/a$ como paso de posprocesado, en cambio, produce el mismo grosor en todas partes.

En flujo compresible hay una capa más. Si arreglas $\alpha$, la densidad y la energía de mezcla deben cambiar con ella para que la termodinámica sea consistente. Por eso los autores no corrigen solo $\alpha$; redistribuyen todo el vector de variables conservativas ($\alpha_1\rho_1$, $\alpha_2\rho_2$, cantidad de movimiento, energía) según las reglas de mezcla. Gracias a esta "regularización consistente con la mezcla", el equilibrio de velocidad, presión y temperatura en la interfaz no se rompe. La prueba de advección de interfaz pura en 1D del artículo muestra justo esto: cuando el agua y el aire se mueven a la misma velocidad, activar el afilado mantiene intactas la presión y la temperatura.

En 2D: mantener redonda la interfaz de la burbuja#

Funcionó en 1D, así que sigue una burbuja en 2D. Las figuras de choque–burbuja del artículo sostienen un círculo de pocas celdas de grosor. Con numerical diffusion activo el borde se vuelve borroso; con sharpening activo y compression a más alto, el borde vuelve a un anillo nítido.

El anillo turquesa se aprieta y el centro de la burbuja se queda. Ni crece ni encoge. Lo que baja es solo el número de celdas del anillo. La compresión sigue $\hat{\mathbf{n}}$, así que la curvatura no se arruga.

Con qué me topé al reproducirlo#

En el fluido puro $|\nabla\alpha|\to 0$ y la normal divide por cero. El artículo lo deja pasar porque $a\,\alpha(1-\alpha)$ es cero allí, pero sin un número pequeño en el denominador el código escribe NaN. Un $d\tau$ grande hace oscilar la interfaz. Hace falta un tope tipo CFL. Las 1 a 3 pasadas por paso del artículo suponen que el paso físico apenas difunde. Una malla gruesa o un choque fuerte puede necesitar más.

Una cosa más. La elegancia del "posprocesado" tiene un precio. El afilado cambia $\alpha$ independientemente de las ecuaciones de gobierno, así que en ese instante la advección no es una solución de la EDP original. A cambio de grosor constante, la conservación local cerca de la interfaz puede desviarse un poco. Los autores lo minimizan con la redistribución consistente con la mezcla, pero no es exactamente cero. El interfoam de OpenFOAM usa un término de compresión de la misma familia ($\nabla\cdot(\alpha(1-\alpha)\mathbf{u}_r)$), lo que sirve de vara de medir útil al portar esta técnica a un solver compresible.

Sin bloquear ε/a el número de celdas vuelve a flotar#

Un afilado en el término fuente reduce el error de difusión y aun así deja que el grosor vague en el espacio y el tiempo. Fijar $\varepsilon/a$ después del paso es lo que mantiene el mismo número de celdas junto al choque. Corregir solo $\alpha$, sin repartir de nuevo las variables conservativas por las reglas de mezcla, rompe primero el equilibrio de presión y temperatura en esa celda.

Si $d\tau$ supera $dx^2/(2\varepsilon)$, la anti-difusión hace oscilar la interfaz en vez de apretarla. El precio de un grosor constante es que, en ese instante, la advección no es solución de la EDP original. La deriva local de conservación es del mismo tipo que los términos no conservativos que atrapa el criterio de Abgrall.