Cuando una microburbuja arde como una estrella bajo el agua — Keller–Miksis y la fuerza secundaria de Bjerknes#

Una burbuja de aire de 5 μm encuentra una onda de ultrasonido, se hincha y en el instante siguiente queda aplastada a un décimo de su tamaño. En ese colapso el gas interior llega a miles de grados. La sonoluminiscencia en un vaso de agua sale de ese intervalo corto. Un baño de limpieza y un litotriptor usan el mismo colapso.

Keller–Miksis es la ecuación que empuja el radio. Siguiendo a Nagy y Hegedűs (2025), se integra una burbuja y después la fuerza secundaria de Bjerknes entre dos.

“

Artículo: 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

Al llevar el deslizador de amplitud hacia 1,2 bar el radio se hincha una vez, con fuerza; el interruptor de modelo solo cambia la altura del rebote tras el colapso.

El rebote de Rayleigh–Plesset es el más alto. En Keller–Miksis la diferencia se reduce. La velocidad del sonido del líquido escrita en el código es 1481 m/s.

La ecuación que gobierna una sola burbuja#

Si la burbuja es una esfera, la única incógnita es el radio $R(t)$. Al integrar la cantidad de movimiento en un líquido incompresible aparece Rayleigh–Plesset.

$\rho_L$ es la densidad del líquido. Los puntos son derivadas temporales. El lado izquierdo es la inercia del líquido empujado en dirección radial. El derecho es la diferencia de presión sobre la pared.

La presión líquida en la pared, $p_L$, tiene tres términos.

$p_G$ es la presión del gas interior, $\sigma$ la tensión superficial, $\mu_L$ la viscosidad del líquido. El segundo término cierra la burbuja. El tercero se opone al movimiento de la pared.

Un cierre politrópico cierra el gas.

$R_E$ es el radio de equilibrio, $p_0$ la presión ambiente, $n_G$ el exponente. El código usa 1,4. Al partir el radio a la mitad el factor es $2^{3\times 1{,}4}=2^{4{,}2}$, cerca de 18. Ese rebote devuelve el colapso hacia afuera.

El forzamiento viaja en la presión lejana.

Con los valores por defecto del código, $p_A=1{,}2\times 10^5$ Pa y $p_0=1{,}0\times 10^5$ Pa, el suelo de $p_\infty$ es $p_0-p_A=-2\times 10^4$ Pa. En el semiciclo en que la presión absoluta entra en negativo, una ley de aire politrópico no sustituye a una cavidad de vapor. El número de cavitación definido desde la presión de vapor es la escala de ese intervalo. La burbuja se hincha con la presión baja y queda aplastada con la presión alta.

El líquido también se comprime — la corrección de Keller–Miksis#

El punto débil de Rayleigh–Plesset es el líquido incompresible. Cuando la velocidad de la pared en el colapso se acerca a la del sonido, esa hipótesis se rompe. Keller–Miksis conserva una velocidad del sonido finita en la pared, $c_L$, hasta primer orden.

Cuando $\dot{R}/c_L$ es pequeño los paréntesis se recogen en 1 y vuelve Rayleigh–Plesset. Al crecer la velocidad de la pared, esos términos irradian la energía del colapso como sonido y contienen la amplitud. Nagy y Hegedűs ordenan los modelos por el número de Mach de la pared. Rayleigh–Plesset es de orden cero, Keller–Miksis de orden uno, Gilmore de orden dos.

La diferencia de rebote en el interruptor de arriba es ese sonido. $c_L=1481$ m/s es la velocidad del sonido en el agua dentro del código.

Integración directa del colapso#

El estado baja a $[R,\dot{R}]$ y se integra con Runge–Kutta de cuarto orden. El colapso es tan agudo que el paso es de un nanosegundo.

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), con el término R-ddot aislado del lado izquierdo
    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)   # protege contra el colapso singular
        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

Los comentarios dejan $R_{\max}/R_E=3{,}71$ y $R_{\min}/R_E=0{,}10$. La burbuja llega a 3,7 radios de equilibrio y después a un décimo. En el radio máximo el factor corre al revés. $3{,}71^{4{,}2}$ vale cerca de 250, así que la presión del gas baja a cerca de 1/250 de $p_{G0}$, unos 500 Pa. La burbuja está casi vacía, y el colapso es la presión exterior cayendo en ese volumen vacío.

Al meter las constantes del código en ese 0,10 la presión se cierra. Con $\sigma=0{,}0725$ y $R_E=5\times 10^{-6}$ m, $2\sigma/R_E=2{,}9\times 10^4$ Pa. Entonces $p_{G0}=p_0+2\sigma/R_E=1{,}29\times 10^5$ Pa, cerca de 1,29 atmósferas. En el radio mínimo la razón de volumen es $(R_E/R)^{3 n_G}=10^{4{,}2}\approx 1{,}6\times 10^4$. La presión del gas es $1{,}29\times 10^5\times 1{,}6\times 10^4\approx 2\times 10^9$ Pa, unas 20 000 atm. La tensión superficial en ese mismo radio, $2\sigma/R$, es diez veces la de equilibrio, $2{,}9\times 10^5$ Pa, y sigue siendo un rastro junto a 20 000 atm. El rebote es el gas, no la película.

La razón de temperatura con el mismo exponente es $(R_E/R)^{3(n_G-1)}=10^{1{,}2}\approx 16$. «Miles de grados» coincide con un arranque a temperatura ambiente por ese factor. «Unos pocos miles de atmósferas» no coincide con la fórmula de presión. La fórmula queda cerca de 20 000 atm.

La ventana de integración está escrita en la llamada. $\Delta t=5\times 10^{-10}$ s y 200 000 pasos cubren $10^{-4}$ s. A 60 kHz el periodo es $1/60000$ s, así que la ventana son exactamente seis ciclos.

Las burbujas se empujan y se atraen entre sí#

El artículo trata de un par. Una burbuja que vibra radia una onda de presión. Bajo una hipótesis incompresible, la presión de la burbuja $j$ a una distancia $r$ escala con la aceleración del volumen.

Esa presión se suma al $p_\infty$ de la vecina y acopla las dos ecuaciones del radio. El promedio temporal es la fuerza secundaria de Bjerknes.

$V_i$ es el volumen, $D$ la distancia entre centros. En fase, se atraen. En oposición, se repelen. La fase sigue a la frecuencia de forzamiento frente a la frecuencia de Minnaert de cada burbuja.

Con los valores del código $n_G=1{,}4$, $p_0=10^5$ Pa y $\rho_L=998$, $\sqrt{3 n_G p_0/\rho_L}\approx 20{,}5$. En $R_E=5$ μm, $f_0\approx 6{,}53\times 10^5$ Hz, cerca de 650 kHz. Un forzamiento de 60 kHz es cerca de un onceavo de esa resonancia. La fase opuesta se abre cuando la $f_0$ de la otra burbuja queda por debajo de 60 kHz, lo que pide un radio por encima de unos 54 μm. Poner 60 kHz entre una burbuja de 5 μm y esa más grande deja a la pequeña bajo la resonancia y a la grande sobre ella.

Al separar los dos deslizadores de radio y colocar el forzamiento entre sus frecuencias de Minnaert, la flecha pasa de atracción a repulsión.

Radios parecidos siguen en fase en cualquier frecuencia de aquí. Verde es atracción, rojo es repulsión. Solo una compañera por encima de unos 54 μm deja que 60 kHz caiga entre las dos resonancias.

Dónde falla el modelo esférico#

Keller–Miksis y Gilmore suponen una esfera. Nagy y Hegedűs enfrentan esa hipótesis a la simulación numérica directa con ALPACA.

En una burbuja aislada con un colapso suave el modelo esférico queda sorprendentemente cerca. La presión del colapso sigue el caso incluso cuando la burbuja no es una esfera perfecta.

Cuando el número de Mach de la pared se acerca a 1, Gilmore se queda más cerca del DNS que Keller–Miksis. Empuja la compresibilidad del líquido a través de una ecuación de estado.

Cuando un par colapsa con fuerza a corta distancia aparece un chorro, líquido que atraviesa un lado. La esfera se rompe y un modelo esférico exagera la presión interior. Ese caso es de DNS.

El modelo esférico es barato, y suele bastar para la fracción de vacío de una nube o para un sonorreactor. Cuando una burbuja choca con una pared o con una vecina y forma un chorro, ese límite llega antes que el número.

Un chorro, y las 20 000 atm eran el número de la esfera#

Las 20 000 atm de $R_{\min}/R_E=0{,}10$ y $n_G=1{,}4$ son el valor politrópico mientras la burbuja sigue siendo una esfera. Cuando un chorro atraviesa la pared, esa esfera ya no está. La presión queda como una sobreestimación, y el colapso de un par cercano es un problema de DNS.

En la misma burbuja de 5 μm, $2\sigma/R_E=2{,}9\times 10^4$ Pa, 0,29 de una atmósfera. Incluso mientras la corrección acústica se sostiene, sigue abierto si esta presión de la película aprieta un paso explícito más que Keller–Miksis. La ventana de tiempo que abre la tensión superficial implícita es el sitio de esa pregunta.