[논문 리뷰] 물과 공기를 반씩 섞었더니 소리가 24 m/s로 기었다 — Murrone–Guillard 5식 모델의 K∇·u 항

1930년 교과서의 한 줄짜리 식#

1930년, 영국 해군 연구소의 A. B. Wood는 음향학 교과서 한 권을 냈다. 그 책에 기포 섞인 물의 음속 식이 한 줄 들어 있다. 혼합물의 압축률은 각 상 압축률의 부피 가중 합이라는 식이다. 이 식대로라면 물과 공기를 반씩 섞은 거품의 음속은 약 24 m/s다. 물(1625 m/s)보다도, 공기(342 m/s)보다도 훨씬 느리다.

이 글은 Murrone–Guillard의 INRIA 보고서 RR-4778(2003)을 읽는다. 이 보고서는 압축성 2상 유동의 5식 축약 모델(reduced model, 7식 모델에서 이완 극한으로 줄인 모델)을 다룬다. 질문은 하나다. 체적분율 식에 붙은 $K\nabla\cdot u$ 항은 무엇을 하는가. 답을 Python 1D 계산으로 확인한다. 같은 관, 같은 펄스에 $K$를 넣고 빼서 소리의 속도를 잰다.

속도를 맞춘 뒤 압력 이완으로 식을 줄이는 길은 Saurel(2009)의 압력 이완에 있다. 오늘은 그 이완을 끝까지 보냈을 때 체적분율 식에 남는 $K$를 잰다.

α 슬라이더를 옮기면 관 A의 Wood 펄스와 관 B의 빠른 펄스가 한 화면에서 갈라진다.

α 슬라이더를 0.5에 두면 관 A의 펄스는 기어가고 관 B의 펄스는 수십 바퀴를 돈다. α를 1% 근처로 옮기면 파란 곡선이 이미 물 음속의 1/14까지 떨어진 것이 보인다.

2003년 보고서가 다시 꺼낸 질문#

Wood 식은 오래된 결과다. 문제는 수치 모델이 그 식을 스스로 재현하느냐다. 7식 Baer–Nunziato 모델은 두 상의 속도와 압력이 따로 있다. 물리적으로 정확하지만 비용이 크다. 두 상이 빠르게 같은 속도·같은 압력으로 묶인다면 변수를 줄일 수 있다.

Kapila 등(2001)은 이 극한에서 5식 모델을 얻었다. Murrone–Guillard는 같은 모델을 점근 전개(asymptotic expansion, 작은 매개변수로 해를 급수 전개하는 기법)로 다시 유도했다. 그리고 쌍곡성(hyperbolicity, 실수 파속을 갖는 성질)과 파동 구조를 정리하고, 두 가지 수치 기법을 시험했다. 식을 더 줄여 계면 평형만 남기는 선택은 4-방정식 모델이 가져간 쪽이다. 오늘은 다섯 식을 유지한 채 $K$ 한 줄만 켜고 끈다.

보고서의 그림 하나가 오늘 글의 출발점이다. 공기–물 혼합물의 음속을 공기 체적분율에 대해 그린 그림이다. 곡선은 두 개다. "5 Eq. Reduc"은 축약 모델, "5 Eq. Trans"는 체적분율을 단순 수송(transport)만 하는 모델이다. 두 곡선은 크게 벌어진다.

이완을 무한히 빠르게 하면 남는 다섯 식#

7식 모델에서 속도 이완과 압력 이완의 시간 척도를 0으로 보내면 다음이 남는다.

$\alpha_k$는 상 $k$의 체적분율, $\rho_k$는 상 밀도, $u$·$p$는 공통 속도·압력, $E$는 혼합물 총에너지다. 앞의 네 식은 보존형이다. 마지막 식만 비보존형이고, 우변의 $K\nabla\cdot u$가 오늘의 주인공이다.

K∇·u는 어디서 오는가#

압력이 평형이면 각 상은 같은 압력 변화를 겪는다. 충격파가 없는 구간에서 상 $k$는 등엔트로피로 압축된다. 그러면 $D\ln\rho_k/Dt = (\rho_k c_k^2)^{-1} Dp/Dt$다. 여기에 상 질량 보존 $D(\alpha_k\rho_k)/Dt = -\alpha_k\rho_k\nabla\cdot u$를 합친다.

두 상의 식을 더하면 $\alpha_1+\alpha_2=1$이므로 좌변이 0이다. 그러면 압력 식이 나온다.

$c_W$가 Wood 음속이다. 이를 첫 식에 다시 넣으면 $K$가 정해진다.

$\rho_k c_k^2$는 상 $k$의 체적탄성률(bulk modulus, 압축에 버티는 정도)이다. 물은 공기보다 약 1.9만 배 단단하다. 그래서 공기를 상 1로 두면 $K\approx\alpha_1$이다. 혼합물이 압축되면 그 부피 변화를 거의 전부 공기가 떠안는다는 뜻이다.

Wood 음속: 압축률은 부피로, 관성은 밀도로#

Wood 식이 느린 이유는 두 성질이 서로 다른 상에서 오기 때문이다. 압축률은 부드러운 공기가 정한다. 관성(밀도)은 무거운 물이 정한다. 공기 50%면 압축률은 순수 공기의 절반쯤이다. 밀도는 순수 공기의 약 400배다. 음속은 $\sqrt{1/(\text{밀도}\times\text{압축률})}$이므로 공기 음속의 $1/\sqrt{200}$ 근처로 떨어진다.

$K=0$으로 두면 이야기가 달라진다. 체적분율은 유체를 따라 흘러가기만 한다. 이때 혼합물은 하나의 stiffened gas(물처럼 단단한 액체를 근사하는 상태방정식)처럼 행동한다. 음속은 $\rho\xi c^2 = \sum_k \alpha_k\xi_k\rho_k c_k^2$, $\xi_k = 1/(\gamma_k-1)$로 주어진다. 물의 큰 체적탄성률이 에너지 가중으로 섞여 들어온다. 그래서 공기가 반이어도 소리가 빠르다.

코드: 같은 관에 두 모델을 걸었다#

상태방정식은 공기 $\gamma=1.4,\ \pi=0$, 물 $\gamma=4.4,\ \pi=6\times10^8$ Pa다. 1 m 관에 공기 50% 혼합물을 채운다. 중앙에 진폭 0.1%의 가우스 압력 펄스를 놓는다. 오른쪽 펄스가 0.3 m를 갈 시간 동안 계산하고 마루 위치로 속도를 잰다. 보존 변수에는 Rusanov 플럭스, 체적분율에는 중앙차분 + 같은 크기의 인공 확산을 쓴다.

import math
 
# stiffened gas: 상 1 = 공기, 상 2 = 물
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):
    """alpha를 수송만 할 때(K = 0): 혼합물은 하나의 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                        # Rusanov 속도 = 그 모델 자신의 음속
    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)]   # 기울기 0 유령 셀
    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

위 코드를 브라우저에서 직접 돌려보자. 같은 Rusanov 계산이 두 관에서 한 스텝씩 진행된다.

기본값(N = 400, α = 0.5)에서 Play를 누르면 표의 measured가 23.55와 742.19로 끝난다. N을 100으로 줄이면 두 비율이 함께 내려가지만 둘의 31배 차이는 그대로다. α를 옮기면 출력 표에 없는 혼합비의 펄스 속도를 직접 잴 수 있다.

23.55와 742.19가 말하는 것#

두 계산은 같은 격자, 같은 플럭스, 같은 초기 펄스를 쓴다. 다른 것은 $K$ 한 줄뿐이다. 그런데 펄스 속도는 23.55 m/s와 742.19 m/s로 31.5배 갈린다. 각 값은 해당 모델의 이론 음속과 0.4% 안에서 맞는다. 남은 0.4%는 1차 Rusanov 확산이 마루를 뒤로 끄는 몫이다.

즉 수치 기법이 틀린 것이 아니다. 두 모델이 서로 다른 물리를 풀고 있다. $K=0$ 모델은 이 펄스를 정확하게, 그러나 틀린 속도로 전한다. 격자를 아무리 조여도 742가 23으로 내려오지 않는다. $K$를 지우면 계면 근처 물이 과하게 압축되는 일이 같은 항의 다른 얼굴이다. 오늘 숫자는 균질한 혼합 영역에서 그 항이 남기는 흔적이다.

표의 두 열을 $\alpha = 0.5$에서 나누면 $745.3/23.6 \approx 31.6$이다. 펄스 측정은 $742.19/23.55 \approx 31.5$로 거의 같다. 각 측정은 이론의 0.996배다. 23.65에서 23.55를 빼면 0.10 m/s, 745.30에서 742.19를 빼면 3.11 m/s다. 절대 오차는 빠른 쪽이 크지만 비율은 같다. 1차 Rusanov가 마루를 같은 비율로 뒤로 끈 것이다. 스텝 수는 316과 315다. 각 실행이 자기 음속으로 CFL을 잡아서, 느린 소리도 빠른 소리도 같은 걸음 수를 쓴다. 차이는 걸린 시간이고, 그 비는 음속 비와 같다.

$K$도 출력에 있다. $\alpha = 0.5$에서 $K = 0.499947$이다. 두 체적탄성률이 무한히 벌어지면 이 식의 $K$는 0.5가 된다. 인쇄된 값은 0.5보다 $5.3\times10^{-5}$ 작다. 물이 공기보다 약 1.9만 배 단단하다는 말과 맞다. $\alpha = 0.01$에서 Wood 음속은 118.6 m/s, 순수 물 1624.9 m/s의 약 1/14다. $\alpha = 0.0001$이면 956.6 m/s로, 1624.9에서 668 m/s가 이미 빠진다. 같은 자리의 $K=0$ 열은 1624.3 m/s다. 기포가 음속을 깎는 일은 $K$가 있을 때만 일어난다.

최소 23.65 m/s는 $\alpha = 0.501$이다. 0.5의 23.6 m/s와 사실상 같다. 체적탄성률 비가 크면 최소가 0.5에 붙는다. $\alpha = 0.99$에서 $K=0$은 540.9 m/s로, 순수 공기 341.6 m/s보다 199 m/s 빠르다. 물 1%의 $\pi$가 공기 쪽으로 샌 결과다. 같은 $\alpha$의 Wood 값은 112.4 m/s라, 순수 공기보다 오히려 느리다.

공기 1%에서 99%까지, 두 곡선이 갈라지는 자리#

왼쪽 열은 공기 한 줌에 민감하다. 0.01%만 있어도 1625 m/s가 957 m/s가 되고, 1%면 119 m/s다. 최소가 $\alpha\approx0.5$에 붙는 것도 체적탄성률 차이가 크기 때문이다.

오른쪽 열은 공기 99%에서 순수 공기보다 빠르다. 물 1%의 $\pi$가 공기 쪽으로 샌 결과이고, 보고서의 "5 Eq. Trans"가 축약 모델에서 벗어나는 자리도 여기다. 물방울이 음파를 거의 되돌리는 이야기는 음향 임피던스와 sharp 계면에서 따로 봤다. 오늘의 두 곡선은 계면 반사 이전에, 균질 혼합물이 어떤 음속을 갖는지를 묻는다.

보고서가 실험과 맞춰 본 것#

보고서는 음속 곡선에서 멈추지 않는다. 공기–물 충격관에서 두 모델의 상 밀도와 내부에너지를 나란히 그린다. 그리고 에폭시–스피넬 분말 혼합물의 충격 Hugoniot를 그린다. 이는 충격파 속도 $U_s$와 입자 속도 $U_p$의 관계다. 그림에는 5식 축약 모델, 7식 모델, 실험점이 함께 실려 있다.

이 비교가 중요한 이유가 있다. $K\nabla\cdot u$는 비보존 항이라 충격파 점프 조건이 모델만으로 유일하게 정해지지 않는다. 수치 점성의 형태가 충격파 뒤 상태를 바꿀 수 있다. 그래서 축약 모델이 7식 모델과 실험에 얼마나 가까운지를 그림으로 확인하는 것이다. 약한 파동에서는 오늘 계산처럼 차이가 사라진다. 강한 충격에서는 이 검증이 별도로 필요하다.

23 m/s가 아니면 그 5식 코드는 K를 버린 것이다#

균질 혼합물 $\alpha=0.5$에 작은 압력 펄스를 놓고 속도를 잰다. 23 m/s 근처면 $K$가 살아 있다. 745 m/s 근처, 곧 표의 $K=0$ 이론값이면 체적분율은 수송만 되고 있다. CFL에 넣는 음속도 같이 본다. Wood 음속 코드에 동결 음속을 넣으면 시간 간격이 수십 배 보수적이 되고, $K=0$ 모델에 23.6 m/s를 넣으면 CFL을 넘긴다. 그 강성을 음향과 수송 중 어디로 보낼지는 Peluchon의 음향-수송 IMEX가 고른 자리와 겹친다.

공기 0.1%만 돼도 Wood 음속은 365 m/s로, 물 1625 m/s의 1/4 아래다. 캐비테이션 초반이나 기포가 들어오는 펌프에서 이 숫자가 압력파 도착 시간을 정한다. 그 칸이 1600 m/s로 찍혀 있으면 모델이 $K$ 없는 쪽으로 가 있다. 0.01%에서도 이미 957 m/s까지 떨어지는데, 순수 물 음속을 그대로 쓰면 도착 시간이 틀린다.