El agua con sonido a 1625 m/s apareció a 18 m/s en el subsistema de presión — adónde lleva la rigidez el splitting de Toro–Vázquez

El agua a 1625 m/s, calculada a 18 m/s#

El sonido viaja en el agua a unos 1625 m/s. Sin embargo, el subsistema de presión de este artículo le asigna al agua a 1 bar una velocidad de onda de 18.44 m/s. No es un error de cálculo. Sale directamente de una sola fórmula.

El artículo es el de Schropff, Petitpas y Daniel en JCP 2026, "Semi-implicit method for two-phase flows in pressure disequilibrium". Es un esquema IMEX (mezcla implícita–explícita) para dos fases con presiones distintas y una sola velocidad. Al pasar el material a Water o Spinel las líneas naranjas quedan casi verticales; los controles de presión y de Mach mueven $A$ y el cociente de pasos de tiempo a la vez.

Las líneas grises discontinuas son las características del sistema completo $u\pm c$. Las líneas naranjas son las $\pm A$ del subsistema de presión, y la azul es la velocidad de advección $u$. Las barras de la derecha muestran, en escala logarítmica, cuánto crece el paso de tiempo a bajo Mach. Cuanto más vertical queda la línea naranja en el agua, más ha bajado $A/c$ de 0.01.

1993: separar la presión de la convección#

El punto de partida es el splitting de Zha–Bilgen (flux splitting: dividir el flujo según su papel físico y tratar cada parte por separado). Corta el flujo de Euler en dos piezas. Una es la parte convectiva, transportada por la velocidad. La otra es la parte que empuja la presión.

En 2012, Toro y Vázquez reescribieron este splitting en forma de Godunov. Así se puede plantear un problema de Riemann para cada subsistema por separado. Schropff y colaboradores llevan la idea a un modelo bifásico con desequilibrio de presiones. Los flujos y términos no conservativos ligados a la presión van al primer subsistema. Los términos de transporte restantes van al segundo.

2026: dos fases, dos subsistemas#

Las variables primitivas son $W=(\alpha_1,\alpha_2,\rho_1,\rho_2,p_1,p_2,u)^T$. El subsistema de presión es:

$\alpha_k$ es la fracción volumétrica de la fase $k$, $\rho_k$ su densidad y $e_k$ su energía interna específica. $p=\sum_k\alpha_k p_k$ es la presión de mezcla.

El subsistema de advección no contiene presión. Fracciones volumétricas, masa, energía interna y cantidad de movimiento se transportan con $u$. Por eso los siete autovalores valen $u$. Sumar ambos subsistemas recupera el modelo original. Los términos de relajación se omiten en el análisis.

Durante el paso de presión la densidad queda congelada#

Volviendo al subsistema de presión, $\alpha_k$ y $\alpha_k\rho_k$ son constantes en el tiempo. Entonces $\rho_k$ tampoco cambia. La única vía para cambiar la presión es la energía interna. En forma primitiva:

$\Gamma_k$ es el parámetro de Grüneisen (cuánto sube la presión al añadir energía a densidad fija). $e_{p,k}$ es $\partial e_k/\partial p_k$ a densidad fija. Combinando esto con la ecuación de cantidad de movimiento aparecen ondas a $\pm A$:

Aquí $\rho=\sum_k\alpha_k\rho_k$ es la densidad de mezcla. $A$ no contiene $u$. Por eso el artículo afirma que esta estructura de ondas es siempre subsónica.

Con la ecuación de estado de gas rigidizado (stiffened gas), $\Gamma_k=\gamma_k-1$. Eso da la ecuación (14) del artículo:

$c_f$ es la velocidad del sonido congelada (frozen sound speed, la de la mezcla sin relajación entre fases) y $p_{\infty,k}$ la presión de rigidez del gas rigidizado. Para un gas ideal, $A^2=c^2-p/\rho$, y el aire da $A/c=\sqrt{(\gamma-1)/\gamma}=0.535$. El agua es otra historia. El término $\gamma p_\infty$ ocupa casi todo $c^2$, y esa parte desaparece por completo.

A y c de cuatro materiales en una sola tabla#

Con los parámetros de gas rigidizado de la tabla 1 del artículo (aire, agua, epoxi, espinela) se calculan $A$ y $c$. El mismo código aísla el subsistema de presión del aire y ejecuta una prueba de pulso. Compara Lax–Friedrichs explícito con Euler implícito hacia atrás al mismo número de Courant. El último bloque es la razón de pasos de tiempo a bajo Mach.

import math
 
# Stiffened-gas parameters (gamma, p_inf [Pa]) from the paper's Table 1
MATS = {"Air": (1.4, 0.0), "Water": (4.4, 6e8), "Epoxy": (2.43, 5.3e9), "Spinel": (1.62, 141e9)}
RHO = {"Air": 1.2, "Water": 1000.0, "Epoxy": 1185.0, "Spinel": 3622.0}
 
def split_speeds(gam, pinf, p, rho):
    c = math.sqrt(gam * (p + pinf) / rho)     # full sound speed
    A = math.sqrt((gam - 1.0) * p / rho)       # pressure-subsystem speed, Eq. (14)
    return c, A
 
def thomas_solve(rhs, beta):
    # (I - beta*L) x = rhs, compact Laplacian, zero-gradient ends
    n = len(rhs)
    cp, dp = [0.0] * n, [0.0] * n
    b0 = 1 + beta
    cp[0], dp[0] = -beta / b0, rhs[0] / b0
    for i in range(1, n):
        b = 1 + beta if i == n - 1 else 1 + 2 * beta
        m = b + beta * cp[i - 1]
        cp[i] = -beta / m
        dp[i] = (rhs[i] + beta * dp[i - 1]) / m
    x = [0.0] * n
    x[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x
 
def ddx(f, dx):
    n = len(f)
    return [(f[min(i + 1, n - 1)] - f[max(i - 1, 0)]) / (2 * dx) for i in range(n)]
 
def pulse_run(courant, implicit, n=200, travel=0.3):
    # linearised pressure subsystem (air, 1 bar): p_t + K u_x = 0, u_t + p_x / rho = 0
    rho, K = 1.2, 0.4 * 1e5
    A = math.sqrt(K / rho)
    dx = 1.0 / n
    x = [(i + 0.5) * dx for i in range(n)]
    p = [100.0 * math.exp(-((xi - 0.5) / 0.04) ** 2) for xi in x]
    u = [0.0] * n
    dt = courant * dx / A
    steps = math.ceil(travel / A / dt)
    for _ in range(steps):
        if implicit:
            ux = ddx(u, dx)
            p = thomas_solve([p[i] - dt * K * ux[i] for i in range(n)], dt * dt * K / rho / dx ** 2)
            px = ddx(p, dx)
            u = [u[i] - dt / rho * px[i] for i in range(n)]
        else:  # Lax-Friedrichs
            P = [p[0]] + p + [p[-1]]; U = [u[0]] + u + [u[-1]]
            p = [0.5 * (P[i] + P[i + 2]) - dt * K * (U[i + 2] - U[i]) / (2 * dx) for i in range(n)]
            u = [0.5 * (U[i] + U[i + 2]) - dt / rho * (P[i + 2] - P[i]) / (2 * dx) for i in range(n)]
        if max(abs(v) for v in p) > 1e9:
            return steps, float("inf"), float("nan")
    right = max(range(n // 2, n), key=lambda i: p[i])
    speed = (x[right] - 0.5) / (steps * dt)
    return steps, max(abs(v) for v in p[n // 2:]), speed
 
print("material  p[bar]     c[m/s]     A[m/s]      A/c")
for name, (g, pinf) in MATS.items():
    for pbar in (1.0, 100.0):
        rho = RHO[name] * (pbar if name == "Air" else 1.0)  # isothermal air
        c, A = split_speeds(g, pinf, pbar * 1e5, rho)
        print(f"{name:7s} {pbar:7.0f} {c:10.1f} {A:10.2f} {A / c:9.5f}")
 
print("\nwater/air mixture, 1 bar, Eq. (13)")
for aw in (0.01, 0.5, 0.99):
    rho = aw * 1000.0 + (1 - aw) * 1.2
    A = math.sqrt((aw * 1e5 * 3.4 + (1 - aw) * 1e5 * 0.4) / rho)
    print(f"alpha_w={aw:4.2f}  rho={rho:7.1f}  A={A:7.2f} m/s")
 
A_air = math.sqrt(0.4e5 / 1.2)
print(f"\npressure pulse, air 1 bar, A = {A_air:.1f} m/s, c = {math.sqrt(1.4e5 / 1.2):.1f} m/s")
print("Courant_A steps  exp_peak  exp_speed  imp_peak  imp_speed")
for cr in (0.5, 0.9, 1.5, 5.0, 20.0):
    s, pe, ve = pulse_run(cr, False)
    _, pi, vi = pulse_run(cr, True)
    print(f"{cr:7.1f} {s:6d} {pe:9.2f} {ve:9.1f} {pi:9.2f} {vi:9.1f}")
 
print("\nlow-Mach air: dt ratio at equal Courant number = (|u|+c)/|u|")
for M in (0.3, 0.1, 0.01, 0.001):
    print(f"M={M:6.3f}  ratio={(1 + M) / M:9.1f}")

material  p[bar]     c[m/s]     A[m/s]      A/c
Air           1      341.6     182.57   0.53452
Air         100      341.6     182.57   0.53452
Water         1     1624.9      18.44   0.01135
Water       100     1638.3     184.39   0.11255
Epoxy         1     3296.7      10.99   0.00333
Epoxy       100     3299.8     109.85   0.03329
Spinel        1     7941.3       4.14   0.00052
Spinel      100     7941.6      41.37   0.00521
 
water/air mixture, 1 bar, Eq. (13)
alpha_w=0.01  rho=   11.2  A=  62.00 m/s
alpha_w=0.50  rho=  500.6  A=  19.48 m/s
alpha_w=0.99  rho=  990.0  A=  18.45 m/s
 
pressure pulse, air 1 bar, A = 182.6 m/s, c = 341.6 m/s
Courant_A steps  exp_peak  exp_speed  imp_peak  imp_speed
    0.5    120     25.59     184.1     35.90     181.1
    0.9     67     42.37     183.2     30.52     180.2
    1.5     40    240.14     129.3     25.76     178.0
    5.0     12 102438.55      13.7     16.02     168.9
   20.0      3   3336.76      32.0      9.82     129.3
 
low-Mach air: dt ratio at equal Courant number = (|u|+c)/|u|
M= 0.300  ratio=      4.3
M= 0.100  ratio=     11.0
M= 0.010  ratio=    101.0
M= 0.001  ratio=   1001.0

Al mover el control de Courant de 0.5 a 20 y cambiar entre explícito e implícito, exp_peak e imp_peak de la tabla se vuelven a calcular en pantalla.

Hasta 0.9 las dos velocidades quedan cerca de 180 m/s. Pasado 1.5 el pico explícito salta a 240 Pa. En 20 el implícito solo se aplana, hasta 9.82 Pa, y el número sigue finito.

Para agua a 1 bar, $A/c$ vale 0.01135. El epoxi baja a 0.00333 y la espinela a 0.00052. Si la presión sube a 100 bar, la $A$ del agua crece exactamente 10 veces, porque $A^2$ es proporcional a $p$. En cambio $c$ solo pasa de 1624.9 a 1638.3, un 0.8%. El aire se queda en 0.53452 sin importar la presión.

La mezcla agua–aire también merece atención. Con fracción volumétrica de agua 0.5, $A$ es 19.48 m/s. Casi igual a los 18.45 m/s con 0.99. El agua fija tanto la densidad de mezcla como el numerador $\sum\alpha_k p_k\Gamma_k$. Solo con 99% de aire salta $A$ a 62.00 m/s.

Al subir la presión de 1 bar a 100 bar, la $A$ del agua pasa de 18.44 a 184.39. Exactamente 10 veces. Como $A^2$ escala con $p$, una presión 100 veces mayor es una velocidad de onda 10 veces mayor. En ese mismo tramo $c$ solo va de 1624.9 a 1638.3, el 0.8% que ya está en la página. La $A$ del epoxi también se multiplica por 10, de 10.99 a 109.85, y $c$ queda en 3296.7 y luego en 3299.8. En la espinela $A$ pasa de 4.14 a 41.37 y $c$ de 7941.3 a 7941.6. $A/c$ crece de 0.00052 a 0.00521 y sigue por debajo de 0.01.

El aire es la misma línea a 1 bar y a 100 bar. $c=341.6$, $A=182.57$, $A/c=0.53452$. La densidad del gas ideal se escaló con la presión, así que la temperatura, y por tanto $c$, queda fija. Por eso $\sqrt{(\gamma-1)/\gamma}=\sqrt{0.4/1.4}$ aparece escrito como 0.535.

Basta un poco de agua para clavar $A$ en el valor del agua. Con $\alpha_w=0.01$, $\rho=11.2$ y $A=62.00$. Con $\alpha_w=0.50$, $\rho=500.6$ y $A=19.48$. Con $\alpha_w=0.99$, $\rho=990.0$ y $A=18.45$, a 0.01 m/s del agua pura en 18.44. El agua se lleva a la vez la densidad de mezcla del denominador y $\sum \alpha_k p_k \Gamma_k$ del numerador. En la tabla del pulso los picos explícitos son 42.37 Pa en Courant 0.9, 240.14 Pa en 1.5 y 102438.55 Pa en 5.0. Los picos implícitos en esas mismas tres celdas bajan: 30.52, 25.76, 16.02 Pa. En Courant 20 el pico implícito es 9.82 Pa, un quinto de los 50 Pa de referencia de la onda a la derecha, y la velocidad medida baja a 129.3 m/s.

La rigidez no desaparece: se traslada entre los pasos#

Entonces, ¿adónde fueron los 1625 m/s del agua? Ningún subsistema contiene $u\pm c$. El de presión tiene $\pm A$; el de advección, solo $u$.

La velocidad del sonido real reaparece cuando los dos pasos se alternan. En el paso de advección cambia $\rho_k$. La ecuación de estado convierte ese cambio en un cambio de presión a razón $\rho c^2$. Esa presión es la entrada del siguiente paso de presión. La rigidez del agua sale de los autovalores de los subsistemas y pasa al recálculo de la ecuación de estado entre pasos. Separar la acústica en un paso y el transporte en el siguiente es también el camino de un IMEX que parte acústica y transporte.

Lo que sigue es una interpretación propia, no una afirmación del artículo. Visto aislado, el paso de presión parece admitir un tratamiento explícito con $\Delta t\le\Delta x/A$. Pero esa desigualdad no garantiza la estabilidad del esquema compuesto. En ese contexto tiene sentido que el artículo resuelva siempre el paso de presión de forma implícita, incluso para el agua, donde $A$ es diminuta.

Un paso de presión implícito mantiene vivo el pulso, pero lo aplana#

El segundo bloque de la salida es la prueba de pulso. Lax–Friedrichs explícito aguanta hasta Courant 0.9 y se rompe en 1.5. En 1.5 la presión máxima llegó a 240 Pa, por encima del valor inicial. En 5.0 superó los 100 000 Pa. El esquema implícito no diverge ni siquiera en 20. A cambio, su pico se aplana a 9.82 Pa, una quinta parte de los 50 Pa exactos.

Ese amortiguamiento es el trato que propone IMEX. En problemas de bajo Mach, la amplitud acústica no es lo que interesa. Basta con que el campo de presión siga al de velocidad de inmediato. Si el problema exige acústica precisa, el trato sale caro.

En el paso implícito de presión, el artículo plantea una ecuación de presión por fase. Como $\Gamma_k$ depende de la presión, la ecuación es no lineal. Se resuelve con iteración de Picard (una linealización que congela algunos coeficientes en el valor de la iteración anterior). La discretización espacial usa un esténcil compacto, así que la matriz es de banda estrecha. La función thomas_solve del código de arriba resuelve justamente una de esas matrices de banda en 1D.

La aritmética de un paso de tiempo 1001 veces mayor a bajo Mach#

El paso de tiempo de un esquema explícito está atado a $\Delta x/(|u|+c)$. En IMEX solo el paso de advección es explícito, así que se relaja a $\Delta x/|u|$. El artículo usa $\Delta t_{IMEX}=0.5\,\Delta x/|\sigma|$. Suponiendo el mismo número de Courant, la razón es:

$M$ es el número de Mach. Con $M=0.3$ son 4.3 veces; con $M=0.001$, 1001 veces. Pero es una cota superior que supone igual costo por paso. El paso implícito resuelve un sistema lineal en cada iteración de Picard. La conclusión del artículo dice lo mismo. La ganancia es grande en problemas de bajo Mach y pequeña o nula en tubos de choque.

Un último paso para devolver la energía que pierden los términos no conservativos#

El splitting deja una tarea de limpieza. Las ecuaciones de energía interna por fase llevan términos no conservativos como $\alpha_k p_k\nabla\cdot u$. Al discretizarlos, la energía total no se conserva exactamente.

Por eso el artículo añade un operador de corrección como tercer paso. Masa, cantidad de movimiento, energía total de mezcla y fracción volumétrica conservan sus valores predichos. Solo a las energías internas por fase se les suma $\alpha_k^n\,\delta\varepsilon$. Así se reparte la diferencia para que cuadre con la energía total de mezcla. La relajación de presiones se ejecuta aparte, al terminar el paso IMEX.

Elegir el paso de presión explícito solo porque la $A$ del agua es pequeña no sostiene el esquema compuesto#

La $A/c$ de la espinela es 0.00052, la del epoxi 0.00333 y la del agua a 1 bar 0.01135. Las tres quedan en 0.01 o por debajo. El paso de presión, aislado, parece explícito bajo $\Delta t\le\Delta x/A$. Esa desigualdad habla solo del subsistema. La rigidez del esquema compuesto está en $c$, no en $A$. A 1 bar, 18.44 y 1624.9 son números distintos. Con $M=0.3$ la razón de pasos de tiempo es 4.3 veces; con $M=0.001$, 1001 veces. En 1D esa razón queda como ganancia solo cuando $1+1/M$ supera el costo de las iteraciones de Picard.

Euler hacia atrás con Courant 20 recorta el pico del pulso, en esta salida, a 9.82 Pa, un quinto de 50 Pa. Si la amplitud acústica entra en el resultado, ese número tiene que quedar escrito. El mismo aplanamiento de un pico $L_\infty$ bajo un paso de presión implícito aparece en un paso de presión implícito que aplana el pico $L_\infty$. Cuando la amplitud forma parte de la respuesta, el paso de tiempo se ajusta a la acústica, o se deja Euler hacia atrás por un integrador implícito de mayor orden.