Al quitar K el agua se comprimió 1.86 veces más — un error de interfaz que refinar la malla no elimina

La ecuación de transporte de la función de color del VOF compresible arrastra un término extra: $-(\psi+K)\nabla\cdot u$. Ese $K$ se trasladó al código y se imprimió su valor. En una celda de interfaz agua–aire con $\psi=0.5$ resultó $K=-0.499947$, de modo que $\psi+K=5.3\times10^{-5}$. El coeficiente añadido para que el término fuera correcto había borrado ese término casi por completo.

Este artículo deriva $K$ desde el cierre isobárico, la misma hipótesis que emplea el método de volúmenes finitos multifásico compresible, basado en presión y totalmente acoplado, de Denner. Luego reproduce lo que ocurre en un problema de compresión 1D al anular $K$ o al invertir su signo: el volumen de líquido se reduce 1.86 y 4.23 veces más de lo que permite la física. Finalmente muestra por qué refinar la malla ocho veces elimina solo un 7% de ese error.

Por qué el VOF compresible arrastra $-(\psi+K)\nabla\cdot u$#

En el VOF incompresible la función de color simplemente viaja con el flujo. Con $\nabla\cdot u=0$ basta con advectarla. Cuando el flujo se vuelve compresible, el fluido se expande y se contrae. La función de color es una fracción volumétrica, así que debe seguir esa expansión.

El artículo lo escribe así:

Aquí $\psi$ es la fracción volumétrica de líquido, $u_i$ la velocidad y $K$ un factor de compresibilidad que depende del material. Al pasar el tercer término al otro lado aparece la forma con derivada material.

Con $K=0$ la fracción volumétrica no cambia a lo largo de una trayectoria material. Parece razonable y es incorrecto. El $\nabla\cdot u$ que ve una celda de interfaz es el valor de la mezcla, y la compresibilidad de una mezcla está dominada por la fase blanda.

Con números: a 1 atm el módulo volumétrico del agua es $\rho_l a_l^2 = 2.64\times10^9$ Pa y el del aire $\rho_g a_g^2 = 1.40\times10^5$ Pa. La razón es 18 857. La compresibilidad de la mezcla es una suma ponderada por fracción volumétrica, así que incluso con $\psi=0.5$ el término del aire se lleva prácticamente todo. La velocidad del sonido de la mezcla cae hasta 23.65 m/s: 69 veces más lenta que en el agua y 14 veces más lenta que en el aire.

En la simulación siguiente conviene mover la fracción volumétrica directamente.

Al subir el deslizador de $\psi$ apenas por encima de cero, la velocidad del sonido se desploma de 341 m/s a menos de 100 m/s. La pista inferior muestra el tiempo que tarda un pulso acústico en cruzar el mismo tramo de 2 m. Con solo un 0.1% de líquido en volumen, el carril de la mezcla ya se retrasa de forma visible respecto al del gas puro.

De dónde sale el número $K$#

$K$ sale de una hipótesis: el cierre isobárico. Las dos fases dentro de una celda de interfaz comparten una sola presión y cada una responde con su propio módulo volumétrico. Es el mismo contenido que la hipótesis local de fase única que el artículo usa en ACID.

Siguiendo una celda lagrangiana y exigiendo que se conserve la masa de la fase $k$:

$V_k$ es el volumen que ocupa esa fase y $\rho_k a_k^2$ su módulo volumétrico. Al sumar ambas fases aparece la divergencia de la mezcla.

El corchete es la compresibilidad de la mezcla, y de él sale gratis la velocidad del sonido de Wood. Ahora se deriva $\psi = V_l/V$ y se sustituyen ambos resultados:

Lo que multiplica es la diferencia de las dos compresibilidades. Si ambas fases comparten módulo volumétrico, la fracción volumétrica no se mueve en absoluto. Comparando con $\mathrm{D}\psi/\mathrm{D}t = K\nabla\cdot u$ y cancelando $\mathrm{D}p/\mathrm{D}t$ queda $K$:

El denominador contiene $1/\psi$ y $1/(1-\psi)$, que divergen en celdas puras. En el código conviene escribir la forma simplificada:

Esta versión es simplemente cero en $\psi=0$ y $\psi=1$. El signo del numerador depende por completo de a qué fase se refiere $\psi$. Copiar la fórmula publicada sin verificar esa etiqueta puede dejar el signo invertido. Más adelante se mide ese costo.

Si $\psi+K$ vale $5\times10^{-5}$, ¿se pueden eliminar ambos términos?#

Tomando el límite $\rho_l a_l^2 \gg \rho_g a_g^2$, el numerador tiende a $-\rho_l a_l^2$ y el denominador a $\rho_l a_l^2(1-\psi)$, de modo que

Físicamente es exactamente lo correcto. El agua casi no se comprime, así que su volumen solo debería advectarse. Llevar $(\psi+K)\nabla\cdot u$ a cero es precisamente ese comportamiento.

¿Por qué no omitir ambos términos entonces? Una celda de líquido puro es otro caso. Tiene $\psi=1$ y $K=0$, lo que deja $(\psi+K)\nabla\cdot u = \nabla\cdot u$. El líquido sí se comprime a su propia tasa, y ese término es el que lo entrega. Que la cancelación sea completa también depende de la razón de módulos volumétricos. Cerca de una razón de $10$ —fluidos supercríticos entre sí, o gas contra gas— con $\psi=0.5$ se obtiene $\psi+K=0.0909$. No es cero.

Cuánto cuesta realmente apretar de 1 bar a 10 bar#

Conviene medirlo en una malla euleriana. La mitad izquierda de $[0,1]$ se llena con agua y la derecha con aire. El extremo izquierdo es una pared ($u=0$); el derecho es un contorno abierto por el que entra aire puro. La presión sube linealmente de 1 bar a 10 bar, uniforme en el espacio. Se usa una ecuación de estado stiffened-gas ($\gamma_l=4.4$, $\Pi_l=6\times10^8$ Pa) y las velocidades en las caras se integran a partir de la compresibilidad de mezcla derivada arriba.

El agua nunca alcanza el contorno derecho, así que su masa se conserva. Eso hace que el volumen exacto de líquido sea simplemente $\psi_0\,\rho_l(p_0)/\rho_l(p)$: una solución de referencia gratuita.

import math
 
P0, PEND = 1.0e5, 1.0e6                    # 1 bar -> 10 bar
GAM_L, PI_L, RHO_L0 = 4.4, 6.0e8, 1000.0   # agua (stiffened gas)
GAM_G, PI_G, RHO_G0 = 1.4, 0.0,   1.2      # aire
 
def stiffened_rho(p, gam, pi, rho0):
    return rho0 * ((p + pi) / (P0 + pi)) ** (1.0 / gam)
 
def stiffened_bulk(p, gam, pi):
    return gam * (p + pi)                  # rho*a^2
 
def wood_speed(psi, p):
    rho_m = psi * stiffened_rho(p, GAM_L, PI_L, RHO_L0) \
          + (1 - psi) * stiffened_rho(p, GAM_G, PI_G, RHO_G0)
    comp  = psi / stiffened_bulk(p, GAM_L, PI_L) \
          + (1 - psi) / stiffened_bulk(p, GAM_G, PI_G)
    return math.sqrt(1.0 / (rho_m * comp)), rho_m
 
def compressibility_factor(psi, p, mode):
    bl = stiffened_bulk(p, GAM_L, PI_L)
    bg = stiffened_bulk(p, GAM_G, PI_G)
    k = (bg - bl) * psi * (1 - psi) / (bg * psi + bl * (1 - psi))
    return 0.0 if mode == 'zero' else (-k if mode == 'flip' else k)
 
def squeeze_column(mode, ncell, nstep, psi0=0.5):
    """pared izquierda u=0, entra gas puro por la derecha, p(t) lineal 1->10 bar"""
    dx, dt = 1.0 / ncell, 1.0 / nstep
    psi = [1.0 if (i + 0.5) * dx < psi0 else 0.0 for i in range(ncell)]
    for n in range(nstep):
        p    = P0 + (PEND - P0) * (n * dt)
        pdot = PEND - P0
        bl   = stiffened_bulk(p, GAM_L, PI_L)
        bg   = stiffened_bulk(p, GAM_G, PI_G)
        u = [0.0] * (ncell + 1)            # compresibilidad de mezcla -> du/dx
        for i in range(ncell):
            u[i + 1] = u[i] - dx * pdot * (psi[i] / bl + (1 - psi[i]) / bg)
        nxt = psi[:]
        for i in range(ncell):
            ur, ul = u[i + 1], u[i]
            right  = psi[i + 1] if i + 1 < ncell else 0.0
            fr = ur * (psi[i] if ur > 0 else right)
            fl = (ul * (psi[i - 1] if ul > 0 else psi[i])) if i > 0 else 0.0
            src = (psi[i] + compressibility_factor(psi[i], p, mode)) * (ur - ul) / dx
            nxt[i] = min(1.0, max(0.0, psi[i] - dt * (fr - fl) / dx + dt * src))
        psi = nxt
    exact = psi0 * stiffened_rho(P0, GAM_L, PI_L, RHO_L0) \
                 / stiffened_rho(PEND, GAM_L, PI_L, RHO_L0)
    return sum(psi) * dx, exact
 
print("[1] stiffened-gas properties")
for p in (P0, PEND):
    rl = stiffened_rho(p, GAM_L, PI_L, RHO_L0)
    rg = stiffened_rho(p, GAM_G, PI_G, RHO_G0)
    print(f"  p={p/1e5:4.1f} bar  rho_l={rl:8.3f}  rho_g={rg:6.3f}"
          f"  rho_l*a_l^2={stiffened_bulk(p,GAM_L,PI_L):.3e}"
          f"  rho_g*a_g^2={stiffened_bulk(p,GAM_G,PI_G):.3e}")
 
print("\n[2] K and (psi+K) at 1 bar")
for psi in (0.10, 0.25, 0.50, 0.75, 0.90):
    k = compressibility_factor(psi, P0, 'acid')
    a, _ = wood_speed(psi, P0)
    print(f"  psi={psi:4.2f}  K={k:+9.6f}  psi+K={psi+k:+.3e}  a_wood={a:7.2f} m/s")
 
print("\n[3] liquid-volume change after 1->10 bar compression")
print("  N     exact        ACID         K=0          K flipped    | x exact")
for ncell in (50, 100, 200, 400):
    out = {}
    for mode in ('acid', 'zero', 'flip'):
        v, exact = squeeze_column(mode, ncell, 10 * ncell)
        out[mode] = v - 0.5
    d = exact - 0.5
    print(f"  {ncell:4d}  {d:.4e}  {out['acid']:.4e}  {out['zero']:.4e}  "
          f"{out['flip']:.4e} |  {out['zero']/d:5.3f}  {out['flip']/d:5.3f}")

[1] stiffened-gas properties
  p= 1.0 bar  rho_l=1000.000  rho_g= 1.200  rho_l*a_l^2=2.640e+09  rho_g*a_g^2=1.400e+05
  p=10.0 bar  rho_l=1000.341  rho_g= 6.215  rho_l*a_l^2=2.644e+09  rho_g*a_g^2=1.400e+06
 
[2] K and (psi+K) at 1 bar
  psi=0.10  K=-0.099994  psi+K=+5.891e-06  a_wood=  39.23 m/s
  psi=0.25  K=-0.249982  psi+K=+1.767e-05  a_wood=  27.28 m/s
  psi=0.50  K=-0.499947  psi+K=+5.302e-05  a_wood=  23.65 m/s
  psi=0.75  K=-0.749841  psi+K=+1.590e-04  a_wood=  27.32 m/s
  psi=0.90  K=-0.899523  psi+K=+4.770e-04  a_wood=  39.43 m/s
 
[3] liquid-volume change after 1->10 bar compression
  N     exact        ACID         K=0          K flipped    | x exact
    50  -1.7027e-04  -1.7027e-04  -3.1695e-04  -7.1967e-04 |  1.861  4.227
   100  -1.7027e-04  -1.7027e-04  -3.1559e-04  -7.0606e-04 |  1.853  4.147
   200  -1.7027e-04  -1.7027e-04  -3.1256e-04  -6.7737e-04 |  1.836  3.978
   400  -1.7027e-04  -1.7027e-04  -3.0670e-04  -6.2727e-04 |  1.801  3.684

K = 0, el cambio de signo y el deslizador de razón de módulos cambian la velocidad a la que se encoge la columna de líquido.

Con K = 0 y play, el bloque azul se encoge más rápido que la columna exacta. Al bajar la razón a $1.0$ las tres columnas se superponen, y la lectura de $\psi+K$ se aleja de $0$.

Durante un aumento de presión de diez veces, el agua pierde realmente $1.70\times10^{-4}$ de volumen, un 0.034% de su cantidad inicial. El cierre ACID reproduce ese valor con cinco decimales. Anular $K$ elimina en cambio $3.13\times10^{-4}$: 1.84 veces la compresión física.

La diferencia entre 0.034% y 0.062% suena pequeña. Pero el patrón de comparación no es el volumen total, sino la compresión física. El esquema añadió un 84% de esa cantidad de la nada. En un problema donde el cambio de volumen es la respuesta —colapso de burbuja, cavitación— esa proporción decide el resultado.

La tabla [2] es casi simétrica. En $\psi=0.10$, $K=-0.099994$, $\psi+K=5.891\times10^{-6}$ y la velocidad del sonido es $39.23$ m/s. En $\psi=0.90$, $K=-0.899523$, $\psi+K=4.770\times10^{-4}$ y la velocidad es $39.43$ m/s. Los extremos difieren $0.20$ m/s, mientras $\psi+K$ se abre unas $81$ veces. El suelo es $23.65$ m/s en $\psi=0.50$, $15.58$ m/s más lento que $39.23$, cerca de un $40\%$. La razón de módulos $2.64\times10^{9}/1.40\times10^{5}$ vale $18857$. Esa razón deja $K$ en $\psi=0.5$ junto a $-0.5$, en el $-0.499947$ impreso, así que $\psi+K=+5.302\times10^{-5}$. En la tabla [3] la razón de $K=0$ baja de $1.861$ en $N=50$ a $1.801$ en $N=400$, una caída de $0.060$. La razón con el signo invertido baja de $4.227$ a $3.684$, una caída de $0.543$. El mismo refinamiento de $8\times$, y el error de signo se reduce unas $9$ veces más que el de $K$ anulado. En $N=400$, tanto $1.801$ como $3.684$ siguen lejos de $1$. En $N=50$ la diferencia entre $3.1695\times10^{-4}$ y $1.7027\times10^{-4}$ es $1.4668\times10^{-4}$, que redondea al $1.467\times10^{-4}$ citado arriba. En $N=400$, $3.0670\times10^{-4}$ menos $1.7027\times10^{-4}$ da $1.3643\times10^{-4}$. A $10$ bar la densidad del líquido sube de $1000.000$ a $1000.341$, un aumento de $0.341$, y el gas salta de $1.200$ a $6.215$. Un $K$ ausente le cobra a la interfaz $1.86$ veces esa pequeña compresión del líquido.

¿Refinar la malla ocho veces lo hace desaparecer?#

Las dos últimas columnas de la tabla [3] responden. Pasar $N$ de 50 a 400 movió la razón de $K=0$ de 1.861 a 1.801. En valor absoluto el error pasó de $1.467\times10^{-4}$ a $1.364\times10^{-4}$: una reducción del 7%. Eso no es convergencia.

La razón es que $K=0$ no es un error de discretización. Fijar $K=0$ cambia la ecuación que se resuelve por $\mathrm{D}\psi/\mathrm{D}t=0$. El esquema converge a la solución de otra ecuación, y no hay motivo para que la diferencia se anule cuando $\Delta x\to0$.

El 7% restante proviene de que la banda de interfaz difundida cambia ligeramente de forma con la malla. El error equivale a la integral de $\psi\,\nabla\cdot u$ a través de esa banda, y la difusión numérica del esquema upwind de primer orden ensancha la banda en número de celdas al refinar. El espesor físico de la banda se mantiene casi fijo, así que la integral también. Cuánto puede avanzar esa banda en un paso lo fija aparte el techo de CFL de la advección de interfaz VOF.

Esa propiedad es una mala noticia y una buena. Un estudio de convergencia de malla no detecta este fallo. A la inversa, si dos mallas dejan la razón plana, esa misma planitud indica que el cierre está mal.

Qué aspecto tiene un signo invertido#

La tercera columna es la corrida con el signo invertido. Con $K\to-K$ el factor pasa a $\psi+K \to 2\psi$ y la compresión salta a entre 3.68 y 4.23 veces el valor físico. También se invierte la dirección: la fracción de líquido baja durante la compresión, así que en pantalla una burbuja parece crecer.

El problema es que esto también produce una imagen verosímil. La interfaz sigue suave, la fracción volumétrica se mantiene dentro de $[0,1]$, nada diverge. Verificar los límites monofásicos tampoco lo detecta, porque $K$ vale cero en $\psi\to0$ y $\psi\to1$ con ambos signos.

Así que al transcribir la fórmula publicada conviene fijar en un comentario a qué fase corresponde $\psi$ y luego bloquear el signo con la prueba de compresión 1D anterior. Con 50 celdas corre en un segundo.

Dónde reaparece la misma hipótesis, ahora en la entalpía#

El ACID del artículo (Acoustically-Conservative Interface Discretisation) aplica este cierre isobárico a la densidad y a la entalpía por separado. La densidad es la mitad fácil. Interpolar linealmente las densidades parciales da

y eso es matemáticamente idéntico al cierre isobárico. La conservación de masa y de momento se siguen de inmediato.

La entalpía es distinta. $h$ es una variable primaria que se resuelve dentro del sistema lineal, así que su interpolación no se puede sustituir como la de la densidad. En su lugar el método calcula una entalpía de cara objetivo que lleva incorporada la hipótesis ACID.

$U$ y $D$ son las celdas aguas arriba y aguas abajo, $L_f$ es el peso de interpolación y $\xi_f$ un limitador. La diferencia respecto al valor actual se suma después al lado derecho de la ecuación de energía como corrección diferida.

$K$ y $\delta h_f$ son dos caras de una misma hipótesis. En la ecuación de fracción volumétrica cierra algebraicamente y colapsa en un único coeficiente; en la entalpía, variable primaria, no cierra y queda diferida en un término fuente explícito. La mitad diferida tiene que converger dentro de la iteración, así que se refleja en el número de iteraciones no lineales. Cerrar en una sola presión compartida es el final de un camino que resuelve dos presiones y luego las relaja a una. La diferencia es que este cierre pliega esa relajación en un solo $K$.

Lo que permitió a Denner soltar el solver de Riemann#

Un método basado en densidad mantiene nítida la discontinuidad de contacto resolviendo un problema de Riemann en la interfaz. El método basado en presión totalmente acoplado saca $\rho$ de las variables de solución y amarra $u$, $p$ y $h$ en un único sistema lineal. El desacoplamiento presión–velocidad propio de una malla colocada se contiene con MWI. Con razones de densidad extremas, la densidad de cara se pondera con una media armónica. Otra vía para conservar el contacto es una reconstrucción que elige la amplitud menor, MUSCL-THINC-BVD.

Soltar el solver de Riemann también suelta la maquinaria que se encargaba de la física de la interfaz. Ese hueco hay que rellenarlo a mano, con hipótesis de cierre. $K$ y $\delta h_f$ son el relleno.

Avanzar implícitamente por encima de una velocidad del sonido de mezcla de 23.65 m/s sin difundir la discontinuidad de contacto exige que el cierre sea exacto. Cuando la acústica pasa a implícita, la factura del CFL se muda al transporte. Esa separación es la de un desacople IMEX entre acústica y transporte. El número 1.86 es una medición de cuánto trabajo está haciendo ese cierre.