초음파 속 기포의 반경 진동과 붕괴, 기포쌍 상호작용을 직접 적분한다

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

지름 5마이크로미터짜리 공기 방울 하나가 초음파를 만나면 어떤 일이 벌어질까. 방울은 부풀었다가, 다음 순간 자기 크기의 10분의 1로 짓눌린다. 그 짧은 붕괴에서 내부 기체는 수천 도까지 달아오른다. 물 한 컵 안에서 별빛 같은 섬광이 튀는 소노루미네선스(sonoluminescence·기포 발광)가 여기서 나온다. 세척기가 안경을 닦고, 초음파가 결석을 부수는 것도 같은 붕괴가 한 일이다.

이 폭력적인 반경 진동을 지배하는 방정식이 Keller–Miksis 방정식이다. 이 글은 Nagy와 Hegedűs(2025)의 논문을 따라, 기포 하나의 반경 동역학을 세우고 Python으로 직접 적분한다. 그다음 기포 두 개가 서로 밀고 당기는 이차 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

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

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

\rho_L\left(R\ddot{R} + \frac{3}{2}\dot{R}^2\right) = p_L(R,t) - p_\infty(t)

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

벽에서의 액체 압력 $p_L$ 은 세 항으로 갈라진다.

p_L(R,t) = p_G(t) - \frac{2\sigma}{R} - \frac{4\mu_L\dot{R}}{R}

$p_G$ 는 기포 안 기체 압력, $\sigma$ 는 표면장력, $\mu_L$ 은 액체 점성이다. 둘째 항은 방울을 오므리는 표면장력, 셋째 항은 벽 운동을 거스르는 점성 저항이다.

기체는 단열에 가깝게 압축된다고 보면, 폴리트로픽 관계로 닫힌다.

p_G(t) = \left(p_0 + \frac{2\sigma}{R_E}\right)\left(\frac{R_E}{R}\right)^{3 n_G}

$R_E$ 는 평형 반경, $p_0$ 는 주위 압력, $n_G$ 는 기체의 폴리트로픽 지수다. 반경이 절반이 되면 기체 압력은 $2^{3n_G}\approx 11$ 배로 치솟는다. 이 급격한 반발이 붕괴를 튕겨내는 힘이다.

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

p_\infty(t) = p_0 - p_A\sin(2\pi f t)

$p_A$ 는 초음파 진폭, $f$ 는 주파수다. 압력이 낮아지는 반주기에 방울이 부풀고, 높아지는 반주기에 짓눌린다.

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

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

\left(1-\frac{\dot{R}}{c_L}\right)R\ddot{R} + \frac{3}{2}\left(1-\frac{\dot{R}}{3c_L}\right)\dot{R}^2 = \left(1+\frac{\dot{R}}{c_L}+\frac{R}{c_L}\frac{\mathrm{d}}{\mathrm{d}t}\right)\frac{p_L-p_\infty}{\rho_L}

$c_L$ 은 액체 음속이다. $\dot{R}/c_L$ 이 작으면 괄호가 모두 1로 수렴해 Rayleigh–Plesset로 되돌아간다. 벽 속도가 커질수록 이 항들이 붕괴 에너지를 음파로 내보내며 진폭을 눌러 준다. Nagy와 Hegedűs는 정확도를 벽 마하수 기준으로 등급화한다. Rayleigh–Plesset은 0차, Keller–Miksis는 1차, Gilmore 모델은 2차다.

아래 시뮬레이션에서 직접 조작해보자. 두 모델을 번갈아 켜며 같은 구동을 준다.

Driving amplitude p_A: 1.20 atmFrequency f: 60 kHzEquilibrium radius R_E: 5 um

구동 진폭을 1.2기압까지 올리면 방울이 크게 부풀었다가 날카롭게 붕괴한다. Rayleigh–Plesset은 반동을 과하게 그린다. Keller–Miksis로 바꾸면 붕괴 후 반동이 눈에 띄게 잦아든다. 음파로 빠져나간 에너지가 그만큼이다.

붕괴를 직접 적분한다#

방정식을 상태 $[R,\dot{R}]$ 의 1계 시스템으로 낮추고 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

평형 반경의 3.7배까지 부풀었다가 10분의 1로 짓눌린다. 최소 반경에서 기체 압력은 수천 기압, 온도는 수천 도에 이른다. 소노루미네선스가 나오는 조건이다.

기포는 서로 밀고 당긴다#

논문의 진짜 주제는 기포 하나가 아니라 기포쌍이다. 한 방울이 진동하면 자기 주변 액체로 압력파를 방사한다. 비압축 가정에서 $j$ 번 방울이 거리 $r$ 에 만드는 압력은 부피 가속도에 비례한다.

p_{\mathrm{ac},j}(r,t) = \frac{\rho_L}{r}\left(2\dot{R}_j^2 R_j + R_j^2\ddot{R}_j\right)

이 압력이 옆 방울의 $p_\infty$ 에 더해져 두 반경 방정식을 엮는다. 결합의 결과가 이차 Bjerknes 힘(진동하는 두 기포 사이의 시간 평균 힘)이다. 방향은 두 방울의 위상차가 정한다.

F_B \;\propto\; -\frac{\langle \dot{V}_1\,\dot{V}_2\rangle}{4\pi D^2}

$V_i$ 는 각 기포 부피, $D$ 는 중심 간 거리다. 두 방울이 같은 위상으로 뛰면 서로 끌어당기고, 반대 위상이면 밀어낸다. 위상은 구동 주파수와 각 방울의 공명 주파수 관계로 갈린다. 여기서 공명은 Minnaert 진동수다.

f_0 = \frac{1}{2\pi R_E}\sqrt{\frac{3 n_G p_0}{\rho_L}}

작은 방울일수록 $f_0$ 가 높다. 구동 주파수가 두 방울의 공명 사이에 끼면, 한 방울은 공명 아래(동위상)로, 다른 방울은 공명 위(반위상)로 응답해 둘이 반대 위상이 된다. 그때 두 방울은 서로 밀어낸다.

아래에서 두 기포의 반경과 구동 주파수를 바꿔 보자.

Bubble 1 radius: 40 umBubble 2 radius: 70 umDriving frequency f: 50 kHzSeparation D: 6 x (R1+R2)

두 반경이 비슷하면 어느 주파수에서든 같은 위상으로 응답해 초록 화살표(끌림)가 뜬다. 한 방울을 크게, 다른 방울을 작게 벌린 뒤 구동 주파수를 둘의 공명 사이로 밀어 넣으면 위상차가 벌어지며 빨간 화살표(밀침)로 뒤집힌다.

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

Keller–Miksis도 Gilmore도 대전제는 "구형"이다. Nagy와 Hegedűs는 이 가정이 언제 깨지는지를 ALPACA 다상 유동 솔버의 직접수치모사(DNS)와 견줘 짚는다. 세 가지 결론이 나온다.

첫째, 고립된 방울이 완만히 붕괴할 때는 구형 모델이 놀랍도록 정확하다. 완벽한 구가 아니어도 붕괴 압력을 잘 맞힌다. 둘째, 벽 마하수가 1에 가까워지는 격렬한 붕괴에서는 Gilmore가 Keller–Miksis보다 DNS에 가깝다. 액체 압축성을 상태방정식으로 더 정직하게 다루기 때문이다. 셋째, 기포쌍이 가까이서 강하게 붕괴하면 제트(jet·한쪽으로 뚫고 들어가는 액체 분사)가 생긴다. 구형이 깨지고, 구형 모델은 내부 압력을 과대평가한다. 이때는 DNS가 필요하다.

즉 구형 모델은 값싸고 대체로 쓸 만한 근사다. 기포 구름의 공극률이나 소노반응기 성능을 예측하기엔 충분하다. 다만 방울이 벽이나 이웃과 부딪혀 제트를 낼 때는 한계를 분명히 안다.

기억할 점#

Rayleigh–Plesset은 구형 기포의 기본 모형이고, Keller–Miksis는 액체 음속을 1차로 넣어 붕괴 반동을 정직하게 눌러 준다.
기포쌍은 방사 압력파로 엮이고, 위상차가 이차 Bjerknes 힘의 부호를 정한다. 구동이 두 공명 사이면 밀치고, 아니면 끌린다.
구형 모델은 완만한 붕괴엔 정확하지만, 제트가 생기는 격렬한 근접 붕괴엔 DNS가 답이다.

[논문 리뷰] 물속 미세기포가 별처럼 타오를 때 — Keller–Miksis와 이차 Bjerknes 힘