直接积分超声场中气泡的半径振荡与坍缩，以及气泡对的相互作用

当水中微气泡像恒星般燃烧 — Keller–Miksis 方程与二次 Bjerknes 力#

一个直径 5 微米的气泡遇到超声波会发生什么。气泡先膨胀，下一瞬间被压缩到自身尺寸的十分之一。在这次短暂的坍缩中，内部气体被加热到数千度。一杯水里迸出如星光般的闪光，声致发光（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$ 有限这一事实纳入到一阶。

\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 是零阶，Keller–Miksis 是一阶，Gilmore 模型是二阶。

在下面的模拟中亲自操作一下吧。交替开启两种模型，施加相同的驱动。

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

把驱动振幅提高到 1.2 个大气压，气泡会大幅膨胀后再尖锐地坍缩。Rayleigh–Plesset 把反弹画得过大。切换到 Keller–Miksis，坍缩后的反弹会明显减弱。以声波形式逃逸的能量正是这么多。

直接积分坍缩

把方程降为状态 $[R,\dot{R}]$ 的一阶系统，用四阶 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 的二阶导项孤立到左边
    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 倍，再被压缩到十分之一。在最小半径处，气体压力可达数千个大气压，温度达数千度。这正是产生声致发光的条件。

气泡之间相互推拉

论文真正的主题不是单个气泡，而是气泡对。一个气泡振荡时，会向自身周围的液体辐射压力波。在不可压缩假设下，第 $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 把液体声速纳入一阶，诚实地压低了坍缩反弹。
气泡对通过辐射压力波耦合，相位差决定二次 Bjerknes 力的符号。驱动频率在两个共振之间就排斥，否则就吸引。
球形模型对缓慢坍缩很精确，但对产生射流的剧烈近距离坍缩，DNS 才是答案。

[论文评述] 当水中微气泡像恒星般燃烧 — Keller–Miksis 方程与二次 Bjerknes 力