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

잠수부는 수면 위에서 부르는 소리를 거의 듣지 못한다. 물로 들어가는 음향 에너지는 0.1% 남짓이고, 나머지는 하늘로 돌아간다. 공기 속 물방울에 음파를 쏘면 격자 위에서도 99.9%가 계면에서 돌아선다. 액체를 넣은 압축성 계산의 난관이 이 대비다.

medium 2를 water로 두면 노란 압력파는 거의 되돌아오고, 계면에서 진폭이 두 배가 된다. $\rho_2$를 내리는 만큼 $c_2$를 올리면 그림은 그대로다.

하늘색 속도파는 부호가 뒤집히고, 물 쪽 진폭은 $10^{-4}$ 근처로 죽는다. 압력은 두 배, 속도는 거의 0이다. 계면이 벽처럼 구는 장면이다.

“

논문: 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에 묶인다.

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

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

예전 우회는 액체를 비압축으로 두는 혼합이었다. 시간스텝은 풀리지만 액체 안의 음파가 사라진다. 액적이 소리에 어떻게 반응하는지를 보려던 목적이 같이 사라진다.

이 논문은 방정식을 대류와 음향으로 가른다. 대류는 명시적, 음향은 투영법으로 음함수다. 남는 조건은 다음뿐이다.

$\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인 면으로 두고, 그 위에서 지킬 조건을 직접 건다. 논문의 세 조건은 이렇다.

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

목록에 밀도가 없다. 밀도는 각 상의 상태방정식이 따로 정한다. 계면 근처에서 밀도를 평균 내지 않으니, 두 상 사이에 없는 중간 상태가 생기지 않는다. 뒤의 가짜 음속이 여기서 막힌다.

상태방정식은 3차다.

계수 $a_1, a_2, a_3$는 압력·온도·물성으로 정해진다. van der Waals, Peng-Robinson, Redlich-Soave-Kwong이 모두 이 꼴이다. 기체와 액체를 한 식이 같이 담는다. 완전기체는 액체를 못하고, Tait 식은 기체를 못한다.

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

계면에서 파동이 갈라지는 방식은 음향 임피던스 하나로 끝난다.

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

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

공기와 물을 넣으면 $Z_{\text{air}} = 416$, $Z_{\text{water}} = 1.50 \times 10^{6}$, 비는 3600이다. $R = 0.9994$, 에너지로는 99.89%가 돌아온다. $T = 1.999$라 계면 압력은 거의 두 배다. 앞의 그림에서 $\rho_2$와 $c_2$를 반대로 밀어도 안 변한 이유가 이 곱이다.

계면이 돌려보낸 파는 계산 영역 출구에서 한 번 더 만난다. 출구가 그것을 삼키는지 다시 튕기는지는 NSCBC의 $\sigma$가 정한다.

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

계면을 격자 몇 칸에 걸쳐 번지게 두면, 확산 계면 기법이 늘 만나는 상황이 된다.

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

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

$\alpha = 0.5$이면 $c_m = 24$ m/s다. 공기의 1/14, 물의 1/60이다. 문제 안의 어떤 물질도 그 속도로 신호를 보내지 않는다. 평균이 만든 숫자다. 번진 폭을 손으로 늘리는 경주는 출력 표 바로 뒤에 둔다.

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

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

얻는 것은 압력의 일치와 열전도의 음함수 처리다. 상태방정식이 돌려준 압력과 반음함수 보정이 계산한 압력이 매 스텝 같다. 보존변수로 풀면 이 둘이 어긋난다고 알려져 있다. 열전도 항도 음함수로 넘길 수 있어서, 자연대류처럼 열이 주역인 저마하수에서 차이가 크다.

잃는 것은 총에너지 보존이다. 그래서 솔버는 아음속으로 범위를 스스로 자른다. 충격파를 잡을 생각이면 다른 정식화가 필요하다.

압력의 양수성도 보장하지 않는다. 저자들은 이것을 결함보다 여지로 본다. 캐비테이션의 준안정 액체는 실제로 음의 압력을 갖고, van der Waals 같은 상태방정식은 음의 압력에 양의 온도를 대응시킬 수 있다.

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

n 슬라이더를 0에서 8칸으로 올리면 하늘색 곡선이 두 물질 어디에도 없는 골짜기를 판다. $\rho_\ell$을 내리면 골짜기가 메워진다. 출력의 $\alpha = 0.01$ 행, 120.5 m/s가 그 골짜기의 입구다. 밀도 대비를 줄이면 음속 대비를 그대로 둬도 골이 사라진다.

반사 에너지는 $R^{2}$이다. 출력의 $R = 0.99944$를 제곱하면 99.889%로, 적힌 퍼센트와 같다. 빠져나간 쪽은 0.111%라 도입의 0.1%와 같은 자리다. 음속 비 1500/347은 4.32인데 임피던스 비는 3595다. 밀도 비 998/1.2에 음속 비를 곱하면 같은 3595가 나온다. 그 벽은 음속이 아니라 밀도가 세운다.

$\alpha = 0.01$에서 이미 1500이 120.5로 떨어진다. 12.4배다. $\alpha = 0.5$의 24 m/s까지는 거기서 다시 약 5배다. 피해의 대부분은 첫 1%에 있다. 명시적 시간스텝은 그 셀의 음속에 묶이므로, 한 칸만 번져도 열두 배 조여진다.

정상파는 2026 Pa 대 Laplace 116 Pa, 비는 17.5다. 속도 진폭 4.9 m/s는 공기 음속 347 m/s의 1.4%라 마하 0.014다. 앞에서 말한 액체 대류 1 m/s보다 다섯 배 가깝다. $f = 347/(2 \times 0.008) = 21688$ Hz로, 적힌 21.7 kHz와 맞다. 액체를 비압축으로 두면 이 17.5배 교란이 액적 안으로 들어가지 못한다. 계면을 한 칸만 섞으면 120.5 m/s 셀이 21.7 kHz 마디의 자리를 민다.

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

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

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

$\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)다. 그 속도 진폭은 임피던스로 추정된다.

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

이 그림이 앞의 정식화를 필요로 한 이유다. 액체를 비압축으로 두면 액적 내부의 압력파가 없어진다. 계면을 번지게 두면 근처 셀의 24 m/s 가짜 음속이 마디 위치를 흔든다. sharp 계면에 양쪽 상 모두 압축성이어야 이 그림이 나온다.

한 칸의 번짐이 마디를 밀어 버리는 선택#

$\alpha = 0.5$의 24 m/s는 계면을 번진 채로 두면 솔버가 밟는 숫자다. 그 음속을 식 안에 남기는 5방정식 쪽은 Murrone–Guillard의 Wood 음속이 같은 24 m/s를 일부러 끌고 간다.

물방울 정상파에서 그 24 m/s를 남기면 마디가 버티는가, 아니면 한 칸의 120.5 m/s를 지워야 21.7 kHz가 제자리에 서는가?