물의 음속 1625 m/s가 압력 부분계에서는 18 m/s였다 — Toro–Vázquez 분할이 강성을 옮기는 자리

음속 1625 m/s짜리 물이 18 m/s로 계산됐다#

물의 음속은 대략 1625 m/s다. 그런데 이번 논문의 압력 부분계에서 1 bar 물의 파속은 18.44 m/s로 나온다. 계산 실수가 아니다. 식 하나에서 바로 나오는 값이다.

Schropff, Petitpas, Daniel의 2026년 JCP 논문이다. 제목은 "Semi-implicit method for two-phase flows in pressure disequilibrium"이다. 압력이 서로 다른 두 상을 한 속도로 푸는 IMEX(음해–양해 혼합) 기법이다. 물질을 Water나 Spinel로 바꾸면 주황 선이 거의 수직이 되고, 압력·마하 슬라이더가 $A$와 시간 간격 비를 같이 민다.

회색 점선은 원래 계의 특성선 $u\pm c$다. 주황 실선은 압력 부분계의 $\pm A$, 파란 실선은 이송 속도 $u$다. 오른쪽 막대는 저마하에서 시간 간격이 몇 배 커지는지 로그 눈금으로 보여준다. 물에서 주황 선이 수직에 가까울수록 $A/c$가 0.01 아래로 떨어진 것이다.

1993년, 압력을 대류에서 떼어낸 분할#

출발점은 Zha–Bilgen 분할(flux splitting, 플럭스를 성질별로 나눠 따로 다루는 기법)이다. 오일러 방정식의 플럭스를 두 덩어리로 나눈다. 하나는 유속으로 실려 가는 대류 부분이다. 다른 하나는 압력이 미는 부분이다.

Toro와 Vázquez는 2012년에 이 분할을 Godunov형으로 다시 썼다. 각 부분계에 리만 문제를 따로 세울 수 있게 만든 것이다. Schropff 등은 이 아이디어를 압력 비평형 2상 모델로 옮겼다. 압력 효과가 담긴 플럭스와 비보존 항을 첫 번째 부분계에 모은다. 남는 이송 항은 두 번째 부분계로 간다.

2026년, 두 상을 두 부분계로 나누다#

원시 변수는 $W=(\alpha_1,\alpha_2,\rho_1,\rho_2,p_1,p_2,u)^T$다. 압력 부분계는 다음과 같다.

$\alpha_k$는 상 $k$의 체적분율, $\rho_k$는 상 밀도, $e_k$는 비내부에너지다. $p=\sum_k\alpha_k p_k$는 혼합 압력이다.

이송 부분계에는 압력이 없다. 체적분율, 질량, 내부에너지, 운동량이 모두 $u$로 실려 간다. 그래서 고유값 7개가 전부 $u$다. 두 부분계를 더하면 원래 모델이 복원된다. 이완 항은 분석에서 뺐다.

압력 단계에서 밀도는 얼어 있다#

압력 부분계를 다시 보자. $\alpha_k$와 $\alpha_k\rho_k$가 시간에 대해 상수다. 그러니 $\rho_k$도 움직이지 않는다. 압력을 바꾸는 통로는 내부에너지 하나뿐이다. 원시 변수로 쓰면 이렇게 된다.

$\Gamma_k$는 그뤼나이젠 계수(Grüneisen parameter, 밀도를 고정하고 에너지를 넣을 때 압력이 오르는 비율)다. $e_{p,k}$는 밀도 고정 시 $\partial e_k/\partial p_k$다. 이 식과 운동량 식을 묶으면 파속 $\pm A$가 나온다.

여기서 $\rho=\sum_k\alpha_k\rho_k$는 혼합 밀도다. $A$에는 $u$가 없다. 그래서 논문은 이 파동 구조가 언제나 아음속이라고 쓴다.

강화 기체(stiffened gas) 상태방정식이면 $\Gamma_k=\gamma_k-1$이다. 그러면 논문 식 (14)가 된다.

$c_f$는 동결 음속(frozen sound speed, 상간 이완이 없을 때의 혼합 음속), $p_{\infty,k}$는 강화 기체의 강성 압력이다. 이상기체면 $A^2=c^2-p/\rho$다. 공기는 $A/c=\sqrt{(\gamma-1)/\gamma}=0.535$가 된다. 물은 사정이 다르다. $\gamma p_\infty$ 항이 $c^2$의 거의 전부를 차지하고, 그 몫이 통째로 빠진다.

네 물질의 A와 c를 한 표에 놓았다#

논문 표 1의 강화 기체 계수(공기, 물, 에폭시, 스피넬)로 $A$와 $c$를 계산했다. 같은 코드에서 공기의 압력 부분계만 떼어내 펄스 실험도 돌렸다. 양해 Lax–Friedrichs와 음해 후진 오일러를 같은 Courant 수로 비교한다. 마지막은 저마하 시간 간격 비다.

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

Courant 슬라이더를 0.5에서 20으로 밀고 양해와 음해를 바꾸면, 표의 exp_peak와 imp_peak가 다시 계산된다.

0.9까지는 둘 다 파속이 180 m/s 근처에 있다. 1.5를 넘기면 양해 봉우리가 240 Pa로 튀고, 음해는 20에서도 9.82 Pa까지 납작해질 뿐 숫자가 터지지 않는다.

1 bar 물에서 $A/c$는 0.01135다. 에폭시는 0.00333, 스피넬은 0.00052까지 내려간다. 압력을 100 bar로 올리면 물의 $A$는 정확히 10배가 된다. $A^2$가 $p$에 비례하기 때문이다. 반면 $c$는 1624.9에서 1638.3으로 0.8%만 움직인다. 공기는 압력과 무관하게 0.53452로 고정이다.

물–공기 혼합물도 볼 만하다. 물 체적분율 0.5에서 $A$는 19.48 m/s다. 0.99일 때의 18.45 m/s와 거의 같다. 혼합 밀도와 분자 $\sum\alpha_k p_k\Gamma_k$를 모두 물이 정하기 때문이다. 공기가 99%가 되어야 62.00 m/s로 뛴다.

1 bar에서 100 bar로 올리면 물의 $A$는 18.44에서 184.39가 된다. 정확히 10배다. $A^2$가 $p$에 비례해서 압력이 100배가 되면 파속은 10배가 된다. 같은 구간에서 $c$는 1624.9에서 1638.3으로 0.8%만 움직인다. 에폭시의 $A$도 10.99에서 109.85로 10배다. $c$는 3296.7에서 3299.8에 머문다. 스피넬은 $A$가 4.14에서 41.37, $c$는 7941.3과 7941.6이다. $A/c$는 0.00052에서 0.00521로 커져도 0.01을 넘지 않는다.

공기는 1 bar와 100 bar가 같은 줄이다. $c=341.6$, $A=182.57$, $A/c=0.53452$. 이상 기체라 밀도를 압력에 비례시켰고, 온도가 고정되면 $c$도 고정된다. $\sqrt{(\gamma-1)/\gamma}=\sqrt{0.4/1.4}$가 0.535로 적힌 이유와 같다.

혼합물은 물이 조금만 있어도 $A$가 물 값으로 붙는다. $\alpha_w=0.01$에서 $\rho=11.2$, $A=62.00$이다. $\alpha_w=0.50$이면 $\rho=500.6$, $A=19.48$. $\alpha_w=0.99$이면 $\rho=990.0$, $A=18.45$다. 순수 물 18.44와 0.01 m/s 차이다. 분모의 혼합 밀도와 분자의 $\sum \alpha_k p_k \Gamma_k$를 물이 같이 가져간다. 펄스 표의 양해 봉우리는 Courant 0.9에서 42.37 Pa, 1.5에서 240.14 Pa, 5.0에서 102438.55 Pa이다. 음해는 같은 세 칸에서 30.52, 25.76, 16.02 Pa로 오히려 줄어든다. Courant 20의 음해 봉우리 9.82 Pa은 오른쪽 파의 기준 50 Pa의 1/5이고, 측정 속도는 129.3 m/s까지 내려간다.

강성은 사라지지 않고 단계 사이로 옮겨간다#

그렇다면 물의 1625 m/s는 어디로 갔나. 두 부분계 어디에도 $u\pm c$는 없다. 압력 부분계는 $\pm A$, 이송 부분계는 $u$뿐이다.

원래 음속은 두 단계를 번갈아 돌릴 때 돌아온다. 이송 단계에서 $\rho_k$가 바뀐다. 상태방정식이 그 변화를 $\rho c^2$ 비율로 압력 변화로 바꾼다. 그 압력이 다음 압력 단계의 입력이 된다. 물의 강성은 부분계 고유값에서 빠져 단계 사이의 상태방정식 재계산으로 옮겨간 것이다. 음향을 한 단계, 수송을 다음 단계로 떼어 강성을 옮기는 풀이는 음향과 수송을 단계로 가른 IMEX에도 있다.

여기서부터는 논문에 쓰인 문장이 아니라 내 해석이다. 압력 단계 하나만 떼면 $\Delta t\le\Delta x/A$로 양해 처리가 가능해 보인다. 하지만 합성된 스킴의 안정성은 그 부등식이 보장하지 않는다. 논문이 $A$가 작은 물에서도 압력 단계를 무조건 음해로 푸는 것은 이 맥락에서 자연스럽다.

압력 단계를 음해로 풀면 펄스는 버티고 납작해진다#

표의 두 번째 블록이 펄스 실험이다. 양해 Lax–Friedrichs는 Courant 0.9까지 버티다가 1.5에서 무너진다. 1.5에서 최대 압력이 240 Pa로 초기값을 넘었다. 5.0에서는 10만 Pa을 넘었다. 음해는 20에서도 발산하지 않는다. 대신 최댓값이 9.82 Pa로 정답 50 Pa의 1/5로 납작해졌다.

이 감쇠가 IMEX의 계약 조건이다. 저마하 문제에서 음향파의 진폭은 관심 대상이 아니다. 압력장이 속도장에 맞춰 즉시 따라오기만 하면 된다. 음향파를 정확히 추적해야 하는 문제라면 이 거래는 손해다.

논문은 음해 압력 단계에서 상별 압력 방정식을 세운다. $\Gamma_k$가 압력의 함수라 식이 비선형이다. 이를 Picard 반복(일부 계수를 직전 반복값으로 고정해 선형화하는 반복법)으로 푼다. 공간 이산화는 조밀 스텐실이라 행렬이 좁은 띠 행렬이 된다. 위 코드의 thomas_solve가 1D에서 그 띠 행렬 하나를 푸는 셈이다.

저마하에서 Δt가 1001배가 되는 산수#

양해 스킴의 시간 간격은 $\Delta x/(|u|+c)$에 묶인다. IMEX는 이송 단계만 양해라 $\Delta x/|u|$로 풀린다. 논문은 $\Delta t_{IMEX}=0.5\,\Delta x/|\sigma|$를 쓴다. 같은 Courant 수를 가정하면 비율은 다음과 같다.

$M$은 마하수다. $M=0.3$이면 4.3배, $M=0.001$이면 1001배다. 다만 이것은 단계당 비용을 같다고 본 상한이다. 음해 단계는 Picard 반복마다 선형계를 푼다. 논문 결론도 같은 말을 한다. 저마하 문제에서는 이득이 크고, 충격파관 문제에서는 이득이 작거나 없을 수 있다.

비보존 항이 남긴 에너지를 되돌리는 마지막 단계#

분할에는 뒷정리가 하나 따른다. 상별 내부에너지 식에는 $\alpha_k p_k\nabla\cdot u$ 같은 비보존 항이 있다. 이 항은 이산화하면 총에너지를 정확히 보존하지 못한다.

그래서 논문은 보정 연산자를 세 번째 단계로 둔다. 질량, 운동량, 혼합 총에너지, 체적분율은 예측값을 그대로 쓴다. 상별 내부에너지에만 $\alpha_k^n\,\delta\varepsilon$을 더한다. 혼합 총에너지와 맞도록 차이를 나눠 주는 것이다. 압력 이완은 이 IMEX 단계가 끝난 뒤 별도로 수행한다.

물의 $A$만 보고 압력 단계를 양해로 두면 합성 스킴이 안 버틴다#

스피넬의 $A/c$는 0.00052, 에폭시는 0.00333, 1 bar 물은 0.01135다. 셋 다 0.01 근처이거나 그 아래다. 압력 단계만 보면 $\Delta t\le\Delta x/A$로 양해가 되어 보인다. 그 부등식은 부분계 이야기다. 합성된 스킴의 강성은 $A$가 아니라 $c$에 있다. 1 bar 물에서 18.44와 1624.9는 다른 숫자다. $M=0.3$이면 시간 간격 비는 4.3배, $M=0.001$이면 1001배다. 1D에서 $1+1/M$이 Picard 반복 비용을 넘을 때만 이 비가 이득으로 남는다.

Courant 20의 후방 오일러는 이 출력에서 펄스 봉우리를 50 Pa의 1/5인 9.82 Pa까지 깎는다. 음향 진폭이 결과에 들어가면 그 숫자를 적어야 한다. 봉우리가 이렇게 눌리는 현상은 음해 압력 단계에서 $L_\infty$가 납작해지는 경우에서 같은 숫자로 다시 만난다. 진폭을 해석에 쓰면 시간 간격을 음향에 맞추거나, 후방 오일러보다 높은 차수의 음해 적분으로 간다.