球形のポリトロープ圧力は、ジェットが球を割る前までの崩壊の数字です。

水中の微小気泡が星のように燃えるとき — 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 側へ押すと半径が一度大きくふくらみ、モデルの切替は崩壊後の反動の高さだけを変えます。

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

Rayleigh–Plesset 側の反動のほうが高いです。Keller–Miksis に切り替えるとその差が縮みます。コードに書かれた液体の音速は1481 m/s です。

一つの泡を支配する方程式

気泡を完全な球とみると、未知量は半径 $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$ は液体の密度、点は時間微分です。左辺は放射状に押しのけられる液体の慣性で、右辺は壁を内外へ押す圧力差です。

壁での液体圧力 $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$ は指数です。コードの指数は1.4です。半径が半分になると因子は $2^{3\times 1.4}=2^{4.2}$ 、約18倍です。この反発が崩壊を跳ね返します。

駆動は遠方圧力に乗ります。

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

コードの初期値 $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次まで入れます。

\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}

$\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)、R-ddot項を左辺に分離
    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)   # 特異な崩壊を防ぐガード
        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$ は10倍の $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_{\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}}

コードの定数 $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 のあいだへ入れると、矢印が引力から斥力へひっくり返ります。

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

半径が近いと、ここのどの周波数でも同位相です。緑は引力、赤は斥力です。約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倍です。音響補正が持っているあいだにも、この膜の圧力が陽的な時間刻みをより強く締めるかどうかが残っています。陰的表面張力が開く時間窓が、その問いの場所です。

[論文レビュー] 水中の微小気泡が星のように燃えるとき — Keller–Miksisと二次Bjerknes力