超音波中の気泡の半径振動と崩壊、気泡対の相互作用を直接積分します

水中の微小気泡が星のように燃えるとき — 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)、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

平衡半径の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力