[논문 리뷰] 물방울은 음파를 거의 되돌려 보낸다 — 음향 임피던스와 sharp 계면 2상 솔버

액-기 계면이 음파의 99.9%를 되돌리는 이유와 번진 계면의 가짜 음속

잠수부는 수면 위에서 부르는 소리를 듣지 못한다. 물을 통과해 들어오는 음향 에너지가 0.1% 남짓이기 때문이다. 나머지는 전부 하늘로 되돌아간다.

같은 숫자가 압축성 2상 유동 솔버 안에서도 튀어나온다. 물방울 하나를 공기 중에 띄워놓고 음파를 쏘면, 계산 격자 위에서도 파동의 99.9%는 계면에서 돌아선다. 이 극단적인 대비를 수치적으로 감당하는 것이 액체를 포함한 압축성 해석의 진짜 난관이다.

논문: A. Urbano, M. Bibal, S. Tanguy, A semi implicit compressible solver for two-phase flows of real fluids, Journal of Computational Physics (2022). ISAE-SUPAERO / IMFT. Level set + ghost fluid 기반 sharp 계면, 원시변수 반음함수 투영법, van der Waals 상태방정식.

액체가 들어오면 음향 시간스텝이 먼저 죽는다#

압축성 2상 유동을 완전 명시적으로 풀면 시간스텝이 음향 CFL에 묶인다.

\Delta t \le \frac{\Delta x}{\max(|u| + c)}

$c$ 는 음속, $\Delta x$ 는 격자 간격이다. 문제는 $c$ 가 상마다 다르다는 점이다. 공기는 347 m/s, 물은 1500 m/s다. 액체 쪽 셀 하나가 전체 시간스텝을 4배 이상 조여버린다.

그런데 정작 액체 안에서 물질이 움직이는 속도는 1 m/s 수준이다. 마하수(관성속도/음속 비)로 따지면 $10^{-3}$ 이다. 실제로 보고 싶은 현상은 계면의 변형인데, 시간스텝은 아무도 보지 않는 음파가 정한다.

기존 회피책은 액체를 비압축성으로 두는 혼합 접근이었다. 그러면 시간스텝은 풀리지만 액체 안에서 음파가 사라진다. 액적이 음파에 어떻게 반응하는지를 보려던 목적 자체가 없어진다.

이 논문은 그 대신 방정식을 대류 부분과 음향 부분으로 쪼갠다. 대류는 명시적으로, 음향은 투영법으로 음함수 처리한다. 남는 시간스텝 조건은 다음뿐이다.

\frac{1}{\Delta t} = \frac{1}{\Delta t_{\mathrm{conv}}} + \frac{1}{\Delta t_{\mu}} + \frac{1}{\Delta t_{\sigma}}

$\Delta t_{\mathrm{conv}} = \Delta x / \max\lVert \vec{u} \rVert$ 는 대류, $\Delta t_{\mu} = \rho \Delta x^{2} / 2\mu$ 는 점성, $\Delta t_{\sigma} = \tfrac{1}{2}\sqrt{\max(\rho)\Delta x^{3}/\sigma}$ 는 표면장력 조건이다. $c$ 가 어디에도 없다.

계면에서 무엇이 이어지고 무엇이 끊기는가#

sharp 계면을 쓴다는 것은 계면 $\Gamma$ 를 두께 0의 면으로 두고, 그 위에서 만족해야 할 조건을 직접 부과한다는 뜻이다. 논문이 쓰는 세 개는 이렇다.

[\vec{u}]_{\Gamma} = 0, \qquad [p]_{\Gamma} = \sigma \kappa + 2\left[\mu \frac{\partial u_n}{\partial n}\right]_{\Gamma} - \frac{2}{3}\left[\mu \nabla \cdot \vec{u}\right]_{\Gamma}, \qquad [-k \nabla T \cdot \vec{n}]_{\Gamma} = 0

$[\,\cdot\,]_{\Gamma}$ 는 계면을 가로지르는 점프, $\sigma$ 는 표면장력, $\kappa$ 는 곡률, $\mu$ 는 점성계수, $k$ 는 열전도도다. 속도는 이어지고, 압력은 표면장력만큼 끊기고, 열유속은 이어진다. 질량 전달은 없다고 가정하므로 온도도 연속이다( $[T]_{\Gamma}=0$ ).

여기서 눈여겨볼 것은 밀도가 목록에 없다는 점이다. 밀도는 각 상의 상태방정식이 각자 정한다. 계면 근처에서 밀도를 평균 내지 않으므로, 두 상 사이에 물리적으로 존재하지 않는 중간 상태가 생기지 않는다. 뒤에서 볼 가짜 음속 문제가 여기서 막힌다.

상태방정식으로는 3차 EoS를 쓴다.

\rho^{3} + a_{1}\rho^{2} + a_{2}\rho + a_{3} = 0

계수 $a_1, a_2, a_3$ 는 압력·온도·물성으로 정해진다. van der Waals, Peng-Robinson, Redlich-Soave-Kwong이 모두 이 꼴이다. 하나의 식으로 기체 상태와 액체 상태를 동시에 기술할 수 있다는 것이 이 선택의 이유다. 완전기체 식은 액체를 못 다루고, Tait 식은 기체를 못 다룬다.

임피던스 하나가 반사율을 전부 정한다#

계면에서 음파가 어떻게 갈라지는지는 놀랍도록 단순한 양 하나로 끝난다. 음향 임피던스다.

Z = \rho c

매질 1에서 매질 2로 파동이 들어갈 때 압력의 반사·투과 계수는 이렇게 나온다.

R = \frac{Z_2 - Z_1}{Z_2 + Z_1}, \qquad T = 1 + R = \frac{2 Z_2}{Z_2 + Z_1}

$R$ 은 반사파 압력진폭 비, $T$ 는 투과파 압력진폭 비다. 밀도만도 아니고 음속만도 아니다. 둘의 곱만 본다.

공기와 물을 넣어보자. $Z_{\text{air}} = 416$ , $Z_{\text{water}} = 1.50 \times 10^{6}$ 이다. 비는 3600이다. $R = 0.9994$ , 에너지로는 99.89%가 되돌아온다. 그리고 $T = 1.999$ — 계면에서 압력은 거의 두 배가 된다.

아래 시뮬레이션에서 직접 조작해보자.

medium 1

medium 2 →

ρ₂998 kg/m³c₂1500 m/sinterface

What to watch: set medium 1 = air, medium 2 = water. The yellow pressure pulse comes back almost intact and nearly doubles at the interface (T = 1.998), while the cyan velocity pulse flips sign and the liquid barely moves. Now press medium 2 = air — R collapses to 0 and the pulse walks straight through. Only the ratio Z₂/Z₁ matters: slide ρ₂ down and c₂ up so that the product stays put, and the picture does not change at all.

관찰 포인트는 두 가지다. medium 2를 water로 두면 노란 압력 파동은 거의 그대로 되돌아오는데 계면에서는 진폭이 두 배로 솟는다. 반면 하늘색 속도 파동은 부호가 뒤집히고 물 쪽 진폭은 $10^{-4}$ 수준으로 죽는다. 압력은 두 배, 속도는 0 — 계면이 벽처럼 굴고 있다. 그리고 $\rho_2$ 를 내리면서 $c_2$ 를 같은 비율로 올려보면 그림이 전혀 변하지 않는다. 곱만 보기 때문이다.

번진 계면은 존재하지 않는 음속을 만든다#

이제 계면을 격자 몇 칸에 걸쳐 번지게 두면 무슨 일이 벌어지는지 보자. 확산 계면(diffuse interface) 계열 기법이 필연적으로 겪는 상황이다.

한 셀 안에 기체와 액체가 섞이면 그 셀의 밀도와 압축률은 각각 이렇게 평균된다.

\rho_m = \alpha \rho_g + (1-\alpha)\rho_\ell, \qquad \frac{1}{\rho_m c_m^{2}} = \frac{\alpha}{\rho_g c_g^{2}} + \frac{1-\alpha}{\rho_\ell c_\ell^{2}}

$\alpha$ 는 기체 체적분율, 아래첨자 $g$ / $\ell$ 은 기체/액체다. Wood의 관계식이다. 이 조합이 고약한 이유는 혼합물이 기체처럼 무르면서 액체처럼 무겁기 때문이다.

$\alpha = 0.5$ 를 넣으면 $c_m = 24$ m/s가 나온다. 공기의 1/14, 물의 1/60이다. 문제에 등장하는 어떤 물질도 그 속도로 신호를 전달하지 않는다. 이 숫자는 물리가 아니라 평균의 부산물이다.

smear n8 cellsρ_liquid998 kg/m³

What to watch: at n = 0 the two lanes are the same run. Pull n up to 8 cells — a width any capturing scheme reaches within a few hundred steps — and the cyan curve digs a hole down to ~24 m/s that neither material has, so the smeared signal falls behind and arrives late. Drop ρ_liquid and the hole fills in: the dip is driven by the density ratio, not by the sound speeds. The yellow line is the bill an explicit solver pays for that hole.

n을 0에서 8칸으로 올려보자. 8칸은 어떤 포착 기법이든 수백 스텝이면 도달하는 두께다. 하늘색 곡선이 두 물질 어디에도 없는 골짜기를 파고, 아래 레이스에서 번진 신호가 그 구간에서 멈칫한다. $\rho_\ell$ 을 낮추면 골짜기가 메워진다 — 이 함정을 만드는 것은 음속 대비가 아니라 밀도 대비다.

부수 효과는 시간스텝에도 온다. 명시적 솔버라면 그 골짜기 셀에서 $\Delta t \le \Delta x / c_m$ 을 지켜야 한다. 번짐이 CFL을 조인다.

원시변수로 푼다는 선택의 대가#

논문의 정식화에서 논쟁적인 부분은 보존변수 $(\rho, \rho\vec{u}, \rho E)$ 가 아니라 원시변수 $(p, \vec{u}, T)$ 로 푼다는 결정이다. 에너지 방정식을 압력 방정식으로 바꿔 쓴다.

얻는 것은 두 가지다. 첫째, 상태방정식이 주는 압력과 반음함수 압력보정이 계산한 압력이 매 스텝 일치한다. 보존변수 정식화에서는 이 둘이 어긋나는 것이 알려진 문제다. 둘째, 열전도 항을 음함수로 넘길 수 있다. 자연대류처럼 열이 주역인 저마하수 문제에서 이건 큰 차이다.

잃는 것은 명확하다. 총에너지 보존이 나빠진다. 그래서 이 솔버는 아음속으로 적용 범위를 스스로 제한한다. 충격파를 잡을 생각이라면 다른 정식화를 써야 한다.

압력에 대한 양수성(positivity)도 보장하지 않는다. 저자들은 이걸 결함이 아니라 여지로 본다. 캐비테이션 상황의 준안정 액체는 실제로 음의 압력을 갖는다. van der Waals 같은 EoS는 음의 압력에 양의 온도를 대응시킬 수 있다.

Python으로 확인하는 반사율과 Wood 음속#

앞의 두 그림에 들어간 숫자를 직접 찍어보자.

# 공기-물 계면의 반사·투과, 그리고 번진 계면의 Wood 음속
from math import sqrt
 
AIR   = (1.2,   347.0)    # (rho, c)
WATER = (998.0, 1500.0)
 
def impedance_split(m1, m2):
    z1, z2 = m1[0] * m1[1], m2[0] * m2[1]
    R = (z2 - z1) / (z2 + z1)
    return z1, z2, R, 1.0 + R, R * R
 
def wood_speed(alpha, gas=AIR, liq=WATER):
    rho_m = alpha * gas[0] + (1.0 - alpha) * liq[0]
    inv   = alpha / (gas[0] * gas[1] ** 2) + (1.0 - alpha) / (liq[0] * liq[1] ** 2)
    return 1.0 / sqrt(rho_m * inv)
 
z1, z2, R, T, E = impedance_split(AIR, WATER)
print(f"Z_air = {z1:8.1f}   Z_water = {z2:10.1f}   Z2/Z1 = {z2/z1:6.0f}")
print(f"R = {R:.5f}   T = 1+R = {T:.5f}   reflected energy = {100*E:.3f} %")
print(f"velocity amp in water / incident = {(z1/z2)*T:.3e}")
 
print("\nalpha   c_wood [m/s]")
for a in (0.0, 0.01, 0.1, 0.5, 0.9, 0.99, 1.0):
    print(f"{a:5.2f}   {wood_speed(a):8.1f}")

Z_air =    416.4   Z_water =  1497000.0   Z2/Z1 =   3595
R = 0.99944   T = 1+R = 1.99944   reflected energy = 99.889 %
velocity amp in water / incident = 5.562e-04
 
alpha   c_wood [m/s]
 0.00     1500.0
 0.01      120.5
 0.10       40.1
 0.50       24.0
 0.90       39.9
 0.99      114.3
 1.00      347.0

$\alpha = 0.01$ 줄을 보자. 물에 기포가 1%만 섞여도 음속이 1500에서 120 m/s로 떨어진다. 실제 기포수(bubbly water)에서 관측되는 값이고, 동시에 계면을 한 칸만 번지게 해도 격자가 어떤 값을 보게 되는지를 말해준다.

정상파 속의 물방울 — 논문이 마지막에 보여준 것#

논문의 마지막 시험은 높이 $L = 8$ mm의 원통형 공동에 반지름 1 mm 물방울을 띄우고, 위쪽 경계에 압력 진동을 준다.

p(t) = p_0 + \Delta p \sin(\omega t)

$\omega = 2\pi f$ 이고 $f = c_0/(2L) \approx 21.7$ kHz는 공동의 첫 음향 모드다. $\Delta p$ 는 $p_0$ 의 2%, 즉 2026 Pa다. 중력은 0으로 둔다.

정상파가 서면 $z = L/2$ 에 압력 마디(node)가 생긴다. 압력 마디는 속도 배(antinode)다. 그 지점 속도진폭은 임피던스로 바로 추정된다.

\Delta v = \frac{\Delta p}{\rho_0 c_0} \approx \frac{2026}{416} \approx 4.9\ \text{m/s}

물방울은 이 진동 속도장 한가운데 놓인다. 압력은 액체 안에서도 전파되며, 초기 원형 액적의 내외 압력차는 Laplace 압력 $2\sigma/r = 116$ Pa로 일정하다. 이 두 숫자의 크기 차이가 결과를 지배한다. 음향 교란은 2026 Pa, 표면장력이 유지하려는 압력차는 116 Pa다. 한 주기마다 액적은 축 방향으로 늘어났다 눌렸다 하고, 진폭이 커지면 분열로 간다.

주목할 점은 이 계산이 왜 앞의 정식화를 필요로 했는가다. 액체를 비압축성으로 두면 액적 내부의 압력파가 없어진다. 계면을 번지게 두면 계면 근처 셀에서 24 m/s짜리 가짜 음속이 정상파의 마디 위치를 흔든다. sharp 계면 + 양쪽 상 모두 압축성 — 이 조합이라야 이 그림이 나온다.

기억할 점 세 줄#

계면에서 음파의 운명은 $Z = \rho c$ 하나가 정한다. 공기-물은 $Z$ 비가 3600이라 에너지의 99.9%가 반사되고, 계면 압력은 두 배로 뛰며, 액체 쪽 속도는 사실상 0이다.
계면을 격자 몇 칸에 걸쳐 번지게 두면 Wood 음속이 24 m/s까지 내려간다. 두 물질 어디에도 없는 값이며, 파동 도달 시각과 명시적 CFL을 동시에 망가뜨린다. sharp 계면은 이걸 원천 차단한다.
대류는 명시적, 음향은 투영법으로 음함수 — 이 분리 덕에 시간스텝에서 $c$ 가 빠진다. 대가는 원시변수 정식화의 총에너지 보존 저하이고, 그래서 적용 범위는 아음속이다.