Rellenar una interfaz con un fluido que no existe — Ghost Fluid Method

En la celda donde el agua se encuentra con el aire, se vierte aire falso que en realidad no está ahí. Ese fluido fantasma (un fluido ficticio extrapolado a través de la interfaz) no existe en ningún punto de la malla. Y sin embargo es justo ese fantasma el que mantiene estable el cálculo cuando una onda de choque cruza la interfaz. Lo que se copia y lo que se extrapola son trabajos distintos. El GFM original oscila donde ese reparto sale mal.

Al subir el salto de densidad la escalera se empina, y la presión sigue cruzando la línea roja discontinua. "show ghost" enciende o apaga la prolongación más allá de la interfaz.

Con las líneas discontinuas activas, cada fluido lleva su entropía al otro lado. Mover la interfaz deja la curva de presión de una sola pieza. Lo que se corta es la densidad.

Cuando dos materiales chocan en una misma celda#

Gas y líquido, o explosivo y metal: dos materiales con propiedades muy distintas se encuentran sobre la misma malla. Si se usan propiedades promediadas, la interfaz se difumina numéricamente. En una interfaz agua-aire con razón de densidad 1000:1, ese difuminado pronto se convierte en oscilaciones no físicas.

El GFM toma otro camino. No mezcla la interfaz de forma borrosa. Extiende cada fluido como un fluido falso más allá de su propia región y luego corre un solver monomaterial dos veces, como si solo existiera un material. La interfaz se rastrea como el cero de un level set (una función distancia con signo que divide el dominio) $\phi$.

Aquí $\phi$ es la función distancia con signo y su conjunto cero $\phi=0$ es la interfaz.

Qué se rompe y qué continúa al cruzar la interfaz#

Para rellenar bien el fantasma hay que saber exactamente qué es continuo al cruzar la interfaz. La presión $p$ y la velocidad normal $u_n$ son continuas ahí. En cambio, la entropía $s$ y la velocidad tangencial $u_t$ son discontinuas.

Los corchetes $[\cdot]$ denotan el salto al cruzar la interfaz. Estas cuatro líneas son todo el GFM. Las magnitudes continuas se copian; las discontinuas se extrapolan. En la figura de arriba la presión queda lisa y la densidad es el escalón, por ese reparto.

Promediar las dos fases dentro de una celda falla de otra manera. El modelo two-fluid de seis ecuaciones pasa a eigenvalores complejos cuando la razón de densidad es grande, y al refinar la malla la divergencia llega antes. El GFM no hace ese promedio. La celda contiene un solo material.

Tratar el agua como un gas — la EOS stiffened#

El agua es casi incompresible, así que la ecuación de estado de gas ideal no la puede tratar. Para encajar ambos materiales en la misma forma de Mie–Grüneisen, el GFM usa la EOS de gas stiffened.

$\rho$ es la densidad, $e$ la energía interna por unidad de masa, $\gamma$ la razón de calores específicos y $p_\infty$ una constante que codifica la rigidez del material. Un gas ideal tiene $p_\infty = 0$; el agua se ajusta a $\gamma \approx 4.4$, $p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$. La velocidad del sonido resulta:

Gracias a $p_\infty$, a igual presión la velocidad del sonido del agua sale mucho mayor que la del aire. Al unificar la EOS en una sola línea, ambos fluidos entran en el mismo solver.

Los dos toques que rellenan una celda fantasma#

Este es el procedimiento clave. Al resolver la región del fluido A, se rellenan con el fantasma de A las celdas que ocupa B. El toque se divide en dos.

1. Copiar: las magnitudes continuas —presión y velocidad normal— se toman directamente de las celdas reales de B al otro lado de la interfaz.
2. Extrapolar: las magnitudes discontinuas —entropía (y por tanto densidad) y velocidad tangencial— se extrapolan de forma isentrópica desde la celda real de A más cercana, cruzando la interfaz.

La extrapolación isentrópica conserva la entropía $s = p / \rho^\gamma$ y vuelve a despejar la densidad.

$s_A$ es la entropía del fluido A real y $p_\text{copy}$ es la presión copiada. La extrapolación suele hacerse con un método fast marching que propaga los valores rápidamente a lo largo de la normal a la interfaz. Una vez que un barrido completa el campo fantasma de A, se olvida que B existe y se resuelven las ecuaciones de Euler monomaterial de A. Se repite lo mismo para B.

Los datos iniciales ya fijan las densidades fantasma. A la izquierda: densidad 1, presión 1, $\gamma = 1.4$, entropía 1. Al copiar la presión 0.1, la densidad fantasma es $0.1^{1/1.4}$, cerca de 0.193, no el 0.125 real. El solver monomaterial ve una razón cercana a 5.2, no 8.

El gas derecho tiene $\gamma = 1.67$. Su entropía es $0.1 / 0.125^{1.67}$, cerca de 3.22. Copiar la presión 1 sobre esa entropía deja una densidad fantasma de unos 0.496, la mitad del estado real izquierdo. La presión cruza. La entropía se queda en el fluido dueño del fantasma.

Con $p_\infty = 0$, la velocidad del sonido es la raíz cuadrada de $\gamma p / \rho$. La izquierda da cerca de 1.183 y la derecha cerca de 1.156. El espaciado con 400 puntos es $1/399$, cerca de 0.00251. Un CFL de 0.4 da un paso cercano a 0.00085, y $0.14$ dividido por ese paso cae cerca de 165 pasos. Ese tiempo final es menos de un quinto del cruce sónico, cerca de 0.85. Una razón de densidad 8 es suave junto al 1000 del agua-aire, y 0.193 frente a 0.496 ya no coinciden con los estados reales.

Dónde oscila el GFM original#

El GFM original tiene problemas cuando un lado de la interfaz es muy rígido (stiff, como el agua, con $p_\infty$ grande). Cuando una onda de choque o detonación fuerte cruza la interfaz, el estado fantasma construido por copia y extrapolación simples discrepa de la solución de Riemann verdadera. Esa discrepancia aparece como oscilaciones de presión cerca de la interfaz, y las oscilaciones a menudo conducen a la divergencia.

En la simulación de abajo, un pulso de presión atraviesa el fluido rígido (región rosa). Cuanto más sube la stiffness ratio, mayores son las ondulaciones justo después de la interfaz.

Al llevar la stiffness ratio hasta 9, aparece con claridad una oscilación de alta frecuencia justo tras la interfaz. No desaparece al refinar la malla, porque la causa es la condición de contorno, no la malla.

El GFM modificado — fijar los valores con un problema de Riemann en la interfaz#

La solución es resolver el estado fantasma en lugar de adivinarlo. El GFM modificado (MGFM o IGFM) toma los estados reales a cada lado de la interfaz como estados izquierdo y derecho de un problema de Riemann y lo resuelve.

$W_A, W_B$ son los estados reales a cada lado, $\mathcal{R}$ es el solver de Riemann y $p^*, u^*$ son la presión y velocidad de la discontinuidad de contacto. Al asignar el mismo $p^*, u^*$ a ambos fantasmas, las condiciones de continuidad de presión y velocidad se satisfacen exactamente. Con "use modified GFM" activo, la oscilación desaparece con la misma rigidez.

Python — un tubo de choque GFM gas-gas#

Resolvamos un tubo de choque donde se encuentran dos gases ideales ($\gamma_A = 1.4$, $\gamma_B = 1.67$), con GFM de extrapolación isentrópica. Usa un flujo HLL y, en cada paso, dos barridos fantasma fusionados por el level set.

import numpy as np
 
GAMMA = (1.4, 1.67)   # razones de calores específicos de izquierda (aire) y derecha (gas monoatómico)
 
def primitive(U, g):
    rho = U[0]; u = U[1] / rho
    e = U[2] / rho - 0.5 * u * u
    p = (g - 1.0) * rho * e          # EOS stiffened, p_inf = 0 (gas ideal)
    return rho, u, p
 
def conserved(rho, u, p, g):
    E = p / (g - 1.0) + 0.5 * rho * u * u
    return np.array([rho, rho * u, E])
 
def hll(UL, UR, g):
    rL, uL, pL = primitive(UL, g); rR, uR, pR = primitive(UR, g)
    aL = np.sqrt(g * pL / rL); aR = np.sqrt(g * pR / rR)
    sL = min(uL - aL, uR - aR); sR = max(uL + aL, uR + aR)
    FL = np.array([rL * uL, rL * uL * uL + pL, uL * (UL[2] + pL)])
    FR = np.array([rR * uR, rR * uR * uR + pR, uR * (UR[2] + pR)])
    if sL >= 0: return FL
    if sR <= 0: return FR
    return (sR * FL - sL * FR + sL * sR * (UR - UL)) / (sR - sL)
 
def extrapolate_ghost(rho, u, p, mat, side):
    # rellena todo el dominio con el fluido `side`: copia p,u; extrapola densidad isentropicamente
    rg, ug, pg = rho.copy(), u.copy(), p.copy()
    g = GAMMA[side]
    idx = np.where(mat == side)[0]
    lo, hi = idx[0], idx[-1]
    ref = lo if side == 1 else hi      # celda real mas cercana a la interfaz
    s_ref = p[ref] / rho[ref] ** g     # entropia a conservar
    for i in range(len(rho)):
        if mat[i] != side:
            j = min(max(i, lo), hi)            # celda real mas cercana
            ug[i] = u[j]; pg[i] = p[j]         # copia
            rg[i] = (pg[i] / s_ref) ** (1.0 / g)   # extrapolacion isentropica
    return rg, ug, pg
 
def run_ghost_fluid_tube(N=400, t_end=0.14, cfl=0.4):
    x = np.linspace(0, 1, N); dx = x[1] - x[0]
    x0 = 0.5                                   # posicion inicial de la interfaz
    mat = np.where(x < x0, 0, 1)
    rho = np.where(x < x0, 1.0, 0.125)
    u   = np.zeros(N)
    p   = np.where(x < x0, 1.0, 0.1)
    t = 0.0
    while t < t_end:
        a = np.sqrt(np.where(mat == 0, GAMMA[0], GAMMA[1]) * p / rho)
        dt = min(cfl * dx / np.max(np.abs(u) + a), t_end - t)
        new = {}
        for side in (0, 1):                    # dos barridos fantasma
            rg, ug, pg = extrapolate_ghost(rho, u, p, mat, side)
            g = GAMMA[side]
            U = np.array([conserved(rg[i], ug[i], pg[i], g) for i in range(N)]).T
            F = np.zeros((3, N + 1))
            for i in range(1, N):
                F[:, i] = hll(U[:, i - 1], U[:, i], g)
            F[:, 0] = F[:, 1]; F[:, N] = F[:, N - 1]
            new[side] = U - dt / dx * (F[:, 1:] - F[:, :-1])
        U = np.where(mat[None, :] == 0, new[0], new[1])   # fusion por level set
        for i in range(N):
            rho[i], u[i], p[i] = primitive(U[:, i], GAMMA[mat[i]])
        x0 += dt * u[np.argmin(np.abs(x - x0))]           # movimiento de la interfaz
        mat = np.where(x < x0, 0, 1)
        t += dt
    return x, rho, u, p
 
if __name__ == "__main__":
    x, rho, u, p = run_ghost_fluid_tube()
    print(f"t=0.14:  max p = {p.max():.3f}, "
          f"densidad cerca del contacto {rho[180]:.3f} -> {rho[220]:.3f}")

La salida muestra que el salto de densidad sobrevive al cruzar la interfaz mientras la presión se mantiene suave. Dos barridos monomaterial son la columna vertebral del GFM.

Trampas que evitar al escribir el código#

Extrapolar a lo largo de los ejes de la malla mezcla la entropía en una interfaz oblicua. La dirección del barrido es solo la normal a la interfaz. Con una razón de rigidez alta, agua contra aire, la copia y la extrapolación originales casi siempre dejan ondulaciones de presión. En ese caso el fantasma se rellena desde el principio con el estado de Riemann de la interfaz. Sin reinicializar el level set, el cero puede seguir ahí mientras los valores de distancia se desordenan, y la velocidad de la interfaz se desvía.

Cuando la extrapolación que funcionó entre gases se pone sobre agua#

Dos gases ideales con $p_\infty = 0$ ($\gamma = 1.4$ y $1.67$) conservan el escalón de densidad en el contacto con una extrapolación isentrópica simple. El agua es $\gamma \approx 4.4$, $p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$. Al sumar esa constante en el numerador de la velocidad del sonido, a la misma presión, las ondulaciones del demo con stiffness ratio 9 son una versión menor del mismo desajuste. Refinar la malla las deja dentro de la condición de contorno.

Si la presión y la velocidad del fantasma siguen siendo una copia, o pasan al estado de Riemann $p^*, u^*$, se decide antes de contar el salto de densidad. Una densidad fantasma que cae por debajo de cero es un problema distinto de limitador. El limitador de flujo que preserva la positividad muestra cómo la masa puede cuadrar mientras la densidad de una celda no. Conviene seguir por ahí después de una interfaz rígida.