물속 미세기포가 별처럼 타오를 때 — Keller–Miksis와 이차 Bjerknes 힘#

지름 5 μm 공기 방울이 초음파를 만나면, 부풀었다가 다음 순간에 크기의 10분의 1로 짓눌린다. 그 붕괴에서 내부 기체가 수천 도까지 달아오른다. 물컵 안의 소노루미네선스(기포 발광)가 그 짧은 구간에서 나온다. 세척기와 결석 파쇄도 같은 붕괴다.

반경을 미는 식이 Keller–Miksis다. Nagy와 Hegedűs(2025)를 따라 방울 하나를 적분하고, 이어서 두 방울의 이차 Bjerknes 힘까지 간다.

“

논문: D. Nagy, F. Hegedűs, "Assessing the accuracy of the coupled-spherical-bubble approach for bubble pairs in an acoustic field", Ultrasonics Sonochemistry 123 (2025) 107651. DOI: 10.1016/j.ultsonch.2025.107651

진폭 슬라이더를 1.2 bar 쪽으로 밀면 반경이 한 번 크게 부풀고, 모델 토글은 붕괴 뒤 반동의 높이만 바꾼다.

Rayleigh–Plesset 쪽 반동이 더 높다. Keller–Miksis로 넘기면 그 차이가 줄어드는데, 코드에 적힌 액체 음속은 1481 m/s다.

방울 하나를 지배하는 방정식#

기포를 완전한 구로 보면 미지수는 반경 $R(t)$ 하나다. 주변 액체를 비압축으로 두고 운동량을 적분하면 Rayleigh–Plesset이 나온다.

$\rho_L$ 은 액체 밀도, 점 표기는 시간 미분. 좌변은 방사상으로 밀려나는 액체의 관성이고, 우변은 벽을 안팎으로 미는 압력 차다.

벽의 액체 압력 $p_L$ 은 세 항이다.

$p_G$ 는 내부 기체 압력, $\sigma$ 는 표면장력, $\mu_L$ 은 액체 점성. 표면장력 항은 방울을 오므리고, 점성 항은 벽 운동을 거스른다.

기체는 폴리트로픽으로 닫는다.

$R_E$ 는 평형 반경, $p_0$ 는 주위 압력, $n_G$ 는 지수. 코드의 지수는 1.4다. 반경이 절반이 되면 인자는 $2^{3\times 1.4}=2^{4.2}$, 약 18배다. 이 반발이 붕괴를 튕겨 낸다.

구동은 원거리 압력에 실린다.

코드 기본값 $p_A=1.2\times 10^5$ Pa, $p_0=1.0\times 10^5$ Pa 이면 $p_\infty$ 의 바닥은 $p_0-p_A=-2\times 10^4$ Pa 다. 절대압이 음수로 들어가는 반주기에는, 폴리트로픽 공기 법칙이 수증기 공동을 대신하지 못한다. 증기압으로 정의하는 캐비테이션 수가 그 구간의 기준이다. 압력이 낮은 반주기에 방울이 부풀고, 높은 반주기에 짓눌린다.

액체도 압축된다 — Keller–Miksis 보정#

Rayleigh–Plesset의 약한 자리는 액체가 비압축이라는 가정이다. 붕괴 순간 벽 속도가 음속에 가까워지면 가정이 깨진다. Keller–Miksis는 벽의 음속 $c_L$ 이 유한하다는 사실을 1차까지 넣는다.

$\dot{R}/c_L$ 이 작으면 괄호가 1로 모여 Rayleigh–Plesset로 돌아간다. 벽 속도가 커지면 이 항들이 붕괴 에너지를 음파로 내보내 진폭을 누른다. Nagy와 Hegedűs는 벽 마하수로 차수를 가른다. Rayleigh–Plesset은 0차, Keller–Miksis는 1차, Gilmore는 2차다.

앞 토글에서 본 반동 차이가 그 음파다. $c_L=1481$ m/s 가 코드의 물 음속이다.

붕괴를 직접 적분한다#

상태를 $[R,\dot{R}]$ 로 낮추고 4차 Runge–Kutta로 적분한다. 붕괴가 뾰족해 시간 스텝은 나노초다.

import numpy as np
 
P0, RHO, SIGMA, MU, C, KAPPA = 1.0e5, 998.0, 0.0725, 1.0e-3, 1481.0, 1.4
 
def bubble_rhs(y, t, RE, pA, f):
    R, V = y
    w = 2 * np.pi * f
    pG0 = P0 + 2 * SIGMA / RE
    pG = pG0 * (RE / R) ** (3 * KAPPA)
    pL = pG - 2 * SIGMA / R - 4 * MU * V / R
    pInf = P0 - pA * np.sin(w * t)
    # d/dt(pL - pInf), with the R-ddot term isolated to the left side
    dP = (-3 * KAPPA * pG * V / R + 2 * SIGMA * V / R**2
          + 4 * MU * V**2 / R**2 + pA * w * np.cos(w * t))
    num = (-1.5 * (1 - V / (3 * C)) * V**2
           + (1 + V / C) * (pL - pInf) / RHO
           + R * dP / (RHO * C))
    den = (1 - V / C) * R + 4 * MU / (RHO * C)
    return np.array([V, num / den])
 
def integrate_km(RE=5e-6, pA=1.2e5, f=60e3, dt=5e-10, steps=200000):
    y = np.array([RE, 0.0])
    t, Rlog = 0.0, []
    for _ in range(steps):
        k1 = bubble_rhs(y, t, RE, pA, f)
        k2 = bubble_rhs(y + 0.5 * dt * k1, t + 0.5 * dt, RE, pA, f)
        k3 = bubble_rhs(y + 0.5 * dt * k2, t + 0.5 * dt, RE, pA, f)
        k4 = bubble_rhs(y + dt * k3, t + dt, RE, pA, f)
        y = y + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        y[0] = max(y[0], 0.1 * RE)   # guard against singular collapse
        t += dt
        Rlog.append(y[0])
    return np.array(Rlog)
 
R = integrate_km()
print(f"R_max / R_E = {R.max() / 5e-6:.2f}")
print(f"R_min / R_E = {R.min() / 5e-6:.2f}")
# R_max / R_E = 3.71
# R_min / R_E = 0.10

출력 주석의 비는 $R_{\max}/R_E=3.71$, $R_{\min}/R_E=0.10$ 이다. 평형 반경의 3.7배까지 부풀었다가 10분의 1로 짓눌린다. 최대 반경에서는 인자가 거꾸로다. $3.71^{4.2}$ 는 약 250이라, 기체 압력은 $p_{G0}$ 의 약 250분의 1, 500 Pa 안팎까지 내려간다. 방울은 거의 비어 있고, 붕괴는 바깥 압력이 그 빈 공간으로 떨어지는 일이다.

그 0.10에 코드 상수를 넣으면 압력이 닫힌다. $\sigma=0.0725$, $R_E=5\times 10^{-6}$ m 이면 $2\sigma/R_E=2.9\times 10^4$ Pa 다. $p_{G0}=p_0+2\sigma/R_E=1.29\times 10^5$ Pa, 대기압의 약 1.29배다. 최소 반경에서 부피비는 $(R_E/R)^{3 n_G}=10^{4.2}\approx 1.6\times 10^4$ 다. 기체 압력은 $1.29\times 10^5\times 1.6\times 10^4\approx 2\times 10^9$ Pa, 약 2만 기압이다. 같은 자리의 표면장력 $2\sigma/R$ 은 열 배인 $2.9\times 10^5$ Pa 이라, 2만 기압 옆에서는 흔적이다. 튕겨 내는 힘은 막이 아니라 기체다.

온도 비는 같은 지수로 $(R_E/R)^{3(n_G-1)}=10^{1.2}\approx 16$ 이다. 수천 도라는 말은 상온에 이 배율을 얹은 자릿수와 맞다. 압력의 "수천 기압"은 이 식과 안 맞는다. 식은 2만 기압 쪽이다.

적분 창도 코드에 적혀 있다. $\Delta t=5\times 10^{-10}$ s, 20만 스텝이면 총 시간은 $10^{-4}$ s 다. 구동 60 kHz의 주기는 $1/60000$ s 라, 창은 정확히 6주기다.

기포는 서로 밀고 당긴다#

논문의 대상은 기포쌍이다. 한 방울이 진동하면 압력파를 낸다. 비압축에서 $j$ 번 방울이 거리 $r$ 에 만드는 압력은 부피 가속도에 비례한다.

이 압력이 옆 방울의 $p_\infty$ 에 더해져 두 반경 식을 엮는다. 시간 평균이 이차 Bjerknes 힘이다.

$V_i$ 는 부피, $D$ 는 중심 거리. 같은 위상이면 끌고, 반대 위상이면 민다. 위상은 구동 주파수와 각 방울의 Minnaert 진동수가 가른다.

코드 상수 $n_G=1.4$, $p_0=10^5$ Pa, $\rho_L=998$ 이면 $\sqrt{3 n_G p_0/\rho_L}\approx 20.5$ 다. $R_E=5$ μm 에서 $f_0\approx 6.53\times 10^5$ Hz, 약 650 kHz 다. 구동 60 kHz는 그 공명의 약 11분의 1이다. 반대 위상은 다른 방울의 $f_0$ 가 60 kHz 아래일 때 열린다. 그 반경은 약 54 μm보다 커야 한다. 5 μm와 그 큰 방울 사이에 60 kHz를 두면, 작은 쪽은 공명 아래, 큰 쪽은 공명 위다.

두 반경 슬라이더를 벌리고 구동 주파수를 둘의 Minnaert 사이로 넣으면, 화살표가 끌림에서 밀침으로 뒤집힌다.

반경이 비슷하면 어느 주파수에서든 같은 위상이다. 초록은 끌림, 빨강은 밀침. 54 μm를 넘는 짝에서만 60 kHz가 그 사이를 가른다.

구형 모델은 어디서 무너지는가#

Keller–Miksis도 Gilmore도 대전제는 구형이다. Nagy와 Hegedűs는 ALPACA 직접수치모사와 견줘 가정이 깨지는 자리를 짚는다.

고립된 방울이 완만히 붕괴할 때 구형 모델은 놀랍도록 맞다. 완전한 구가 아니어도 붕괴 압력이 따라온다.

벽 마하수가 1에 가까운 격렬한 붕괴에서는 Gilmore가 Keller–Miksis보다 DNS에 가깝다. 액체 압축성을 상태방정식으로 더 밀기 때문이다.

기포쌍이 가까이서 세게 붕괴하면 제트(한쪽으로 뚫고 들어가는 액체)가 생긴다. 구가 깨지고, 구형 모델은 내부 압력을 과대평가한다. 그때는 DNS다.

구형 모델은 싸고, 기포 구름의 공극률이나 소노반응기에는 대개 충분하다. 방울이 벽이나 이웃과 부딪혀 제트를 낼 때는 그 한계가 숫자보다 먼저다.

제트가 뚫리면 2만 기압은 구가 만든 숫자다#

$R_{\min}/R_E=0.10$ 과 $n_G=1.4$ 가 주는 2만 기압은, 구가 끝까지 구로 남았을 때의 폴리트로픽 값이다. 제트가 벽을 뚫으면 그 구는 이미 없다. 압력은 과대평가로 남고, 가까운 쌍의 붕괴는 DNS로 넘기는 쪽이 맞다.

같은 5 μm에서 $2\sigma/R_E=2.9\times 10^4$ Pa 은 대기압의 0.29배다. 음향 보정이 버티는 동안에도 이 막의 압력이 명시 시간스텝을 더 세게 조이는지가 남아 있다. 암시적 표면장력이 여는 시간 창이 그 질문의 자리다.