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

直径5微米的空气泡遇上超声波，先胀开，下一瞬间被压到自身尺寸的十分之一。那次坍缩里，内部气体升到数千度。一杯水里的声致发光，就出在这段很短的时间上。清洗机和碎石，用的是同一次坍缩。

推着半径走的方程是 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$ 有限这件事放到一阶。

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

前面开关上看到的反弹差，就是那道声波。 $c_L=1481$ m/s 是代码里的水中声速。

直接积分坍缩

把状态降成 $[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

输出注释里的比是 $R_{\max}/R_E=3.71$ ， $R_{\min}/R_E=0.10$ 。泡先胀到平衡半径的3.7倍，再被压到十分之一。最大半径上因子反过来。 $3.71^{4.2}$ 大约是250，气体压力降到 $p_{G0}$ 的大约250分之一，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\sigma/R$ 是十倍的 $2.9\times 10^5$ Pa，放在两万个大气压旁边只是痕迹。弹回去的是气体，不是膜。

温度比用同一个指数， $(R_E/R)^{3(n_G-1)}=10^{1.2}\approx 16$ 。数千度，和室温乘上这个倍数的量级相符。压力上的“数千个大气压”和这个式子不相符。式子站在两万个大气压一侧。

积分窗口也写在调用里。 $\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 的驱动大约是这个共振的十一分之一。反相要等另一个气泡的 $f_0$ 低于60 kHz 才打开，半径得大于大约54 μm。把60 kHz 放在5 μm 和那个大泡之间，小的在共振之下，大的在共振之上。

把两个半径滑块拉开，再把驱动频率放进两者 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。

球形模型便宜，对气泡云的空隙率或声化学反应器通常够用。气泡撞上壁面或邻居、射出射流时，这个限度比那个数字更早到来。

射流一穿，两万个大气压就是球自己算出的数

$R_{\min}/R_E=0.10$ 和 $n_G=1.4$ 给出的两万个大气压，是球一直保持为球时的多方值。射流刺穿壁面，那个球已经不在了。压力留成高估，近处气泡对的坍缩该交给 DNS。

同一个5微米气泡上， $2\sigma/R_E=2.9\times 10^4$ Pa，是大气压的0.29倍。声学修正还撑得住的时候，这层膜的压力会不会把显式时间步收得更紧，仍然开放。隐式表面张力打开的时间窗就是这个问题所在的地方。

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