[Revisión de artículo] Mitad agua, mitad aire y el sonido se arrastra a 24 m/s — El término K∇·u del modelo de cinco ecuaciones de Murrone–Guillard

Una fórmula de una línea en un libro de 1930#

En 1930, A. B. Wood, del laboratorio de investigación del Almirantazgo británico, publicó un libro de acústica. Dentro aparece una fórmula de una sola línea para la velocidad del sonido en agua con burbujas. Dice que la compresibilidad de una mezcla es la suma de las compresibilidades de las fases, ponderada por volumen. Con una espuma mitad aire y mitad agua, el resultado es de unos 24 m/s. Es mucho más lento que el agua (1625 m/s) y mucho más lento que el aire (342 m/s).

Esta entrada revisa el informe INRIA RR-4778 (2003) de Murrone y Guillard. El informe estudia un modelo reducido de cinco ecuaciones (obtenido del modelo de siete ecuaciones en el límite de relajación instantánea) para flujo bifásico compresible. Una sola pregunta guía el texto: ¿qué hace realmente el término $K\nabla\cdot u$ de la ecuación de fracción volumétrica? La respuesta sale de un cálculo 1D en Python. Mismo tubo, mismo pulso, con y sin $K$, y se mide la velocidad del sonido.

Reducir las ecuaciones igualando la velocidad y luego relajando la presión está en la revisión de la relajación de presión de Saurel (2009). Aquí se mide la $K$ que queda en la ecuación de fracción volumétrica cuando esa relajación se lleva hasta el final.

El deslizador α separa en la misma pantalla el pulso de Wood del tubo A y el pulso rápido del tubo B.

Con el deslizador α en 0,5, el pulso del tubo A avanza a paso lento mientras el del tubo B da decenas de vueltas. Al mover α cerca del 1%, la curva azul ya cayó a 1/14 de la velocidad del sonido en el agua.

La pregunta que reabrió el informe de 2003#

La fórmula de Wood es un resultado antiguo. La cuestión es si un modelo numérico la reproduce por sí solo. El modelo de siete ecuaciones de Baer–Nunziato mantiene velocidades y presiones separadas para ambas fases. Es físicamente rico, pero costoso. Si las fases se igualan rápido en velocidad y presión, se pueden eliminar variables.

Kapila et al. (2001) obtuvieron un modelo de cinco ecuaciones en ese límite. Murrone y Guillard lo volvieron a derivar con un desarrollo asintótico (expandir la solución en potencias de un parámetro pequeño). Después analizaron su hiperbolicidad (velocidades de onda reales y finitas) y su estructura de ondas, y probaron dos esquemas numéricos. Seguir recortando hasta dejar solo el equilibrio de interfaz es la vía del modelo de cuatro ecuaciones. Aquí se mantienen las cinco ecuaciones y solo se enciende o se apaga la línea de $K$.

Una figura del informe es el punto de partida. Muestra la velocidad del sonido de la mezcla aire–agua frente a la fracción volumétrica de aire. Hay dos curvas. "5 Eq. Reduc" es el modelo reducido. "5 Eq. Trans" es un modelo que solo transporta la fracción volumétrica. Las dos curvas quedan muy separadas.

Cinco ecuaciones tras una relajación infinitamente rápida#

Al llevar a cero los tiempos de relajación de velocidad y presión del modelo de siete ecuaciones, queda esto:

$\alpha_k$ es la fracción volumétrica de la fase $k$, $\rho_k$ su densidad, $u$ y $p$ la velocidad y presión comunes, $E$ la energía total de la mezcla. Las cuatro primeras ecuaciones son conservativas. Solo la última no lo es, y su lado derecho $K\nabla\cdot u$ es el protagonista de esta entrada.

De dónde sale K∇·u#

En equilibrio de presión, ambas fases sufren el mismo cambio de presión. Lejos de los choques, la fase $k$ se comprime de forma isentrópica, así que $D\ln\rho_k/Dt = (\rho_k c_k^2)^{-1} Dp/Dt$. Esto se combina con la conservación de masa de cada fase, $D(\alpha_k\rho_k)/Dt = -\alpha_k\rho_k\nabla\cdot u$:

Al sumar las ecuaciones de las dos fases, como $\alpha_1+\alpha_2=1$, el lado izquierdo se anula. Aparece entonces una ecuación para la presión:

$c_W$ es la velocidad del sonido de Wood. Al sustituirla de vuelta, $K$ queda determinado:

$\rho_k c_k^2$ es el módulo volumétrico (bulk modulus, la resistencia a la compresión) de la fase $k$. El agua es unas 19 000 veces más rígida que el aire. Con el aire como fase 1, $K\approx\alpha_1$. Cuando la mezcla se comprime, el aire absorbe casi todo el cambio de volumen.

Velocidad de Wood: compresibilidad por volumen, inercia por densidad#

La velocidad de Wood es baja porque dos propiedades vienen de fases distintas. El aire, blando, fija la compresibilidad. El agua, pesada, fija la inercia. Con 50% de aire, la compresibilidad es cerca de la mitad de la del aire puro. La densidad es unas 400 veces la del aire. Como la velocidad del sonido escala con $\sqrt{1/(\text{densidad}\times\text{compresibilidad})}$, cae a cerca de $1/\sqrt{200}$ de la del aire.

Con $K=0$ la historia cambia. La fracción volumétrica solo viaja con el flujo. La mezcla se comporta entonces como un único stiffened gas (una ecuación de estado que aproxima líquidos rígidos como el agua). Su velocidad del sonido cumple $\rho\xi c^2 = \sum_k \alpha_k\xi_k\rho_k c_k^2$ con $\xi_k = 1/(\gamma_k-1)$. El enorme módulo volumétrico del agua se cuela a través de una ponderación energética. Por eso el sonido sigue siendo rápido aun con 50% de aire.

Código: dos modelos en el mismo tubo#

Ecuaciones de estado: aire $\gamma=1.4,\ \pi=0$; agua $\gamma=4.4,\ \pi=6\times10^8$ Pa. Se llena un tubo de 1 m con una mezcla de 50% de aire. En el centro se coloca un pulso gaussiano de presión con amplitud de 0,1%. Se calcula hasta que el pulso que viaja a la derecha debería haber recorrido 0,3 m, y la velocidad se mide con la posición del pico. Las variables conservativas usan flujo de Rusanov. La fracción volumétrica usa diferencias centradas más una difusión artificial del mismo tamaño.

import math
 
# stiffened gas: fase 1 = aire, fase 2 = agua
G = (1.4, 4.4)
PI = (0.0, 6.0e8)
P0 = 1.0e5
RHO0 = (1.2, 1000.0)
 
def wood_speed(a1):
    rho = a1 * RHO0[0] + (1 - a1) * RHO0[1]
    comp = a1 / (G[0] * (P0 + PI[0])) + (1 - a1) / (G[1] * (P0 + PI[1]))
    return 1.0 / math.sqrt(rho * comp)
 
def isobaric_speed(a1):
    """solo se transporta alpha (K = 0): la mezcla es un único stiffened gas"""
    rho = a1 * RHO0[0] + (1 - a1) * RHO0[1]
    xi = a1 / (G[0] - 1) + (1 - a1) / (G[1] - 1)
    gpi = a1 * G[0] * PI[0] / (G[0] - 1) + (1 - a1) * G[1] * PI[1] / (G[1] - 1)
    gam = 1 + 1 / xi
    return math.sqrt(gam * (P0 + gpi / xi / gam) / rho)
 
def mix_pressure(a1, rhoe):
    a2 = 1 - a1
    num = rhoe - a1 * G[0] * PI[0] / (G[0] - 1) - a2 * G[1] * PI[1] / (G[1] - 1)
    return num / (a1 / (G[0] - 1) + a2 / (G[1] - 1))
 
def kapila_K(a1, p):
    k1, k2 = G[0] * (p + PI[0]), G[1] * (p + PI[1])
    a2 = 1 - a1
    return a1 * a2 * (k2 - k1) / (a1 * k2 + a2 * k1)
 
def pulse_speed(a_bg, use_K, n=400, travel=0.3, amp=1e-3):
    dx = 1.0 / n
    x = [(i + 0.5) * dx for i in range(n)]
    a1 = [a_bg] * n
    p = [P0 * (1 + amp * math.exp(-((xi - 0.5) / 0.03) ** 2)) for xi in x]
    m1 = [a_bg * RHO0[0]] * n
    m2 = [(1 - a_bg) * RHO0[1]] * n
    mom = [0.0] * n
    E = [a_bg * (pi + G[0] * PI[0]) / (G[0] - 1) + (1 - a_bg) * (pi + G[1] * PI[1]) / (G[1] - 1) for pi in p]
    c = wood_speed(a_bg) if use_K else isobaric_speed(a_bg)
    s = 1.05 * c                        # velocidad de Rusanov = la velocidad del sonido del propio modelo
    t_end = travel / c
    steps = math.ceil(t_end / (0.4 * dx / s))
    dt = t_end / steps
    g = lambda q, i: q[min(max(i, 0), n - 1)]   # celdas fantasma de gradiente nulo
    for _ in range(steps):
        u = [mom[i] / (m1[i] + m2[i]) for i in range(n)]
        pr = [mix_pressure(a1[i], E[i] - 0.5 * mom[i] * u[i]) for i in range(n)]
        U = (m1, m2, mom, E)
        F = ([m1[i] * u[i] for i in range(n)], [m2[i] * u[i] for i in range(n)],
             [mom[i] * u[i] + pr[i] for i in range(n)], [(E[i] + pr[i]) * u[i] for i in range(n)])
        new = []
        for q, f in zip(U, F):
            fl = [0.5 * (g(f, i - 1) + g(f, i)) - 0.5 * s * (g(q, i) - g(q, i - 1)) for i in range(n + 1)]
            new.append([q[i] - dt / dx * (fl[i + 1] - fl[i]) for i in range(n)])
        na = []
        for i in range(n):
            K = kapila_K(a1[i], pr[i]) if use_K else 0.0
            dudx = (g(u, i + 1) - g(u, i - 1)) / (2 * dx)
            dadx = (g(a1, i + 1) - g(a1, i - 1)) / (2 * dx)
            lap = g(a1, i + 1) - 2 * a1[i] + g(a1, i - 1)
            na.append(a1[i] + dt * (-u[i] * dadx + K * dudx + s * lap / (2 * dx)))
        m1, m2, mom, E = new
        a1 = na
    pr = [mix_pressure(a1[i], E[i] - 0.5 * mom[i] ** 2 / (m1[i] + m2[i])) for i in range(n)]
    ipk = max((i for i in range(n) if x[i] > 0.6), key=lambda i: pr[i])
    return c, (x[ipk] - 0.5) / t_end, steps
 
print("alpha_air  wood[m/s]  K=0[m/s]")
for a in [0.0, 1e-4, 1e-3, 0.01, 0.1, 0.5, 0.9, 0.99, 1.0]:
    print(f"{a:9.4f}  {wood_speed(a):9.1f}  {isobaric_speed(a):8.1f}")
amin = min((k / 1000 for k in range(1, 1000)), key=wood_speed)
print(f"Wood minimum: {wood_speed(amin):.2f} m/s at alpha_air = {amin:.3f}")
print(f"K at alpha_air=0.5: {kapila_K(0.5, P0):.6f}")
print("\n1D pulse, alpha_air = 0.5, N = 400")
for useK in (True, False):
    c, meas, st = pulse_speed(0.5, useK)
    tag = "with K" if useK else "K = 0 "
    print(f"{tag}  theory {c:7.2f}  measured {meas:7.2f}  ratio {meas / c:.3f}  steps {st}")

alpha_air  wood[m/s]  K=0[m/s]
   0.0000     1624.9    1624.9
   0.0001      956.6    1624.3
   0.0010      364.8    1618.9
   0.0100      118.6    1567.2
   0.1000       39.4    1228.3
   0.5000       23.6     745.3
   0.9000       39.2     581.7
   0.9900      112.4     540.9
   1.0000      341.6     341.6
Wood minimum: 23.65 m/s at alpha_air = 0.501
K at alpha_air=0.5: 0.499947
 
1D pulse, alpha_air = 0.5, N = 400
with K  theory   23.65  measured   23.55  ratio 0.996  steps 316
K = 0   theory  745.30  measured  742.19  ratio 0.996  steps 315

El código de arriba se puede correr directamente en el navegador. El mismo cálculo de Rusanov avanza paso a paso en ambos tubos.

Con los valores por defecto (N = 400, α = 0,5), al pulsar Play la columna measured termina en 23,55 y 742,19. Al bajar N a 100 ambos cocientes caen juntos, pero la brecha de 31× entre modelos se mantiene. Moviendo α se miden pulsos en proporciones de mezcla que la tabla de salida no incluye.

Lo que dicen 23,55 y 742,19#

Ambas corridas comparten malla, flujo y pulso inicial. Solo cambia la línea de $K$. Aun así, las velocidades del pulso se separan por un factor de 31,5: 23,55 m/s frente a 742,19 m/s. Cada una coincide con la velocidad teórica de su modelo dentro de 0,4%. Ese 0,4% es la difusión de Rusanov de primer orden que arrastra el pico hacia atrás.

Así que el esquema no está mal. Los dos modelos resuelven físicas distintas. El modelo con $K=0$ transporta el pulso con precisión, pero a la velocidad equivocada. Ningún refinamiento de malla baja 742 a 23. Borrar $K$ y sobrecomprimir el agua cerca de una interfaz es la otra cara del mismo término. Los números de hoy son la huella de ese término dentro de una mezcla uniforme.

Al dividir las dos columnas en $\alpha = 0.5$ sale $745.3/23.6 \approx 31.6$. La medición del pulso es $742.19/23.55 \approx 31.5$, casi igual. Cada medición es 0.996 de la teoría. Restar deja 0.10 m/s ($23.65-23.55$) y 3.11 m/s ($745.30-742.19$). El error absoluto es mayor del lado rápido. El cociente no. Rusanov de primer orden arrastra la cresta hacia atrás en la misma fracción. Los pasos son 316 y 315. Cada corrida fija su CFL con su propia velocidad del sonido, así que la onda lenta y la rápida usan el mismo número de pasos. Lo que cambia es el tiempo, en la misma razón que las velocidades del sonido.

$K$ también está en la salida. En $\alpha = 0.5$, $K = 0.499947$. Si los dos módulos volumétricos se separan sin límite, la $K$ de esta fórmula llega a 0.5. El valor impreso queda $5.3\times10^{-5}$ por debajo de 0.5, en línea con que el agua es unas 19 000 veces más rígida que el aire. En $\alpha = 0.01$ la velocidad de Wood es 118.6 m/s, cerca de 1/14 del agua pura a 1624.9 m/s. En $\alpha = 0.0001$ ya es 956.6 m/s, 668 m/s por debajo de 1624.9. La columna $K=0$ en ese mismo sitio marca 1624.3 m/s. Las burbujas recortan la velocidad del sonido solo cuando $K$ está presente.

Del 1% al 99% de aire: dónde se separan las curvas#

La columna izquierda es sensible a una pizca de aire. Con 0,01% la velocidad baja de 1625 m/s a 957 m/s, y con 1% queda en 119 m/s. El mínimo se pega cerca de $\alpha\approx0.5$ porque el contraste de módulos volumétricos es grande.

La columna derecha, con 99% de aire, es más rápida que el aire puro. El 1% de agua filtró su $\pi$ hacia el gas, y ahí también la curva "5 Eq. Trans" del informe se aparta del modelo reducido. Una gota que devuelve casi toda la onda es otro relato, en la entrada sobre impedancia acústica e interfaz sharp. Las dos curvas de hoy preguntan, antes de cualquier reflexión en una interfaz, qué velocidad del sonido tiene de verdad una mezcla uniforme.

Lo que el informe contrasta con experimentos#

El informe no se detiene en las curvas de velocidad del sonido. En un tubo de choque aire–agua grafica lado a lado las densidades de fase y energías internas de ambos modelos. Luego grafica la Hugoniot de choque de una mezcla de polvo epoxi–espinela, la relación entre velocidad del choque $U_s$ y velocidad de partícula $U_p$. La figura superpone el modelo reducido de cinco ecuaciones, el de siete ecuaciones y puntos experimentales.

Esa comparación importa por una razón. $K\nabla\cdot u$ es no conservativo, así que el modelo por sí solo no fija de forma única las condiciones de salto del choque. La forma de la viscosidad numérica puede desplazar el estado detrás del choque. De ahí la necesidad de comprobar, figura por figura, cuánto se acerca el modelo reducido al de siete ecuaciones y al experimento. Para ondas débiles, como en el cálculo de hoy, la diferencia desaparece. Para choques fuertes hace falta una validación propia.

Si el pulso no ronda los 23 m/s, el código dejó K fuera#

Se coloca un pulso pequeño de presión en una mezcla uniforme con $\alpha=0.5$ y se cronometra. Cerca de 23 m/s, $K$ sigue activo. Cerca de 745 m/s, el valor teórico de la columna $K=0$, la fracción volumétrica solo viaja arrastrada. La velocidad del sonido que entra en el CFL va en la misma prueba. Una velocidad congelada en un código que debería correr a la de Wood vuelve el paso decenas de veces demasiado conservador, y meter 23,6 m/s en un modelo con $K=0$ rompe el CFL. A dónde va esa rigidez, acústica o transporte, coincide con el corte del IMEX acústico-transporte de Peluchon.

Con solo 0,1% de aire, la velocidad de Wood ya es 365 m/s, por debajo de un cuarto de los 1625 m/s del agua. En cavitación incipiente, o en una bomba que ingiere burbujas, ese número fija la llegada de la onda de presión. Una celda que todavía imprime 1600 m/s se fue al lado sin $K$. Incluso con 0,01% la velocidad ya bajó a 957 m/s, así que conservar el valor del agua pura desfasa la llegada.