3.4 Solving the Rayleigh–Plesset equation

The Rayleigh–Plesset equation is a stiff nonlinear second-order ODE in the bubble radius R(t)R(t):

ρ[RR¨+32R˙2]=pB(t)p(t)2σR4μR˙R.\rho \left[R \ddot R + \frac{3}{2} \dot R^2 \right] = p_B(t) - p_\infty(t) - \frac{2 \sigma}{R} - \frac{4 \mu \dot R}{R}.

It has no general closed-form solution. Three special regimes do admit analytical treatment:

  1. Small-amplitude linear oscillation about the static equilibrium — produces a harmonic oscillator at the Minnaert frequency, with three dissipation channels (radiation, thermal, viscous).
  2. Pure inertial collapse (Rayleigh 1917) of an empty cavity in inviscid liquid — produces an analytical closed-form R(t)R(t) that diverges in wall speed as R0R \to 0.
  3. Pure inertial growth (Plesset 1949) of a bubble in a tension step that exceeds the Blake threshold — produces an asymptotically linear growth at the Rayleigh velocity 2Δp/3ρ\sqrt{2|\Delta p| / 3 \rho}.

These three regimes anchor intuition. For the full nonlinear problem under arbitrary drive p(t)p_\infty(t), we integrate numerically. This lesson develops all four cases, with the numerical solver running live for the reader to experiment with.

Small-amplitude linear oscillation: the Minnaert frequency

Linearise the Rayleigh–Plesset equation about the stable static equilibrium R=R0R = R_0. Let R(t)=R0+r(t)R(t) = R_0 + r(t) with rR0|r| \ll R_0, expand to first order in rr and r˙\dot r, and use the polytropic relation pG(R)=pG,0(R0/R)3κp_G(R) = p_{G,0}(R_0/R)^{3\kappa} to evaluate dpB/dRdp_B / dR at R0R_0:

Linearisation about static equilibrium

At equilibrium, R˙=R¨=0\dot R = \ddot R = 0 and

pB(R0)p2σ/R0=0.p_B(R_0) - p_\infty - 2\sigma/R_0 = 0.

Perturb. Inertial term to first order:

ρRR¨ρR0r¨.\rho R \ddot R \approx \rho R_0 \ddot r.

The R˙2\dot R^2 term is second-order and drops out. Right-hand side:

pB(R0+r)2σ/(R0+r)4μr˙/(R0+r)pp_B(R_0 + r) - 2\sigma/(R_0 + r) - 4\mu\dot r/(R_0 + r) - p_\infty

[pB(R0)+rdpBdRR0]2σR0(1rR0)4μr˙R0p.\approx \left[p_B(R_0) + r \frac{dp_B}{dR}\bigg|_{R_0}\right] - \frac{2\sigma}{R_0}\left(1 - \frac{r}{R_0}\right) - \frac{4\mu \dot r}{R_0} - p_\infty.

The equilibrium parts cancel by construction. With pB=pv+pG,0(R0/R)3κp_B = p_v + p_{G,0}(R_0/R)^{3\kappa}:

dpBdRR0=3κpG,0R0.\frac{dp_B}{dR}\bigg|_{R_0} = -\frac{3\kappa p_{G,0}}{R_0}.

Collecting:

ρR0r¨=3κpG,0R0r+2σR02r4μr˙R0.\rho R_0 \ddot r = -\frac{3\kappa p_{G,0}}{R_0} r + \frac{2\sigma}{R_0^2} r - \frac{4\mu \dot r}{R_0}.

Rearranging into standard damped-oscillator form r¨+2βr˙+ω02r=0\ddot r + 2\beta \dot r + \omega_0^2 r = 0:

r¨+4μρR02r˙+1ρR02(3κpG,02σR0)r=0.\ddot r + \frac{4\mu}{\rho R_0^2} \dot r + \frac{1}{\rho R_0^2}\left(3 \kappa p_{G,0} - \frac{2\sigma}{R_0}\right) r = 0.

The natural frequency is

ω02=1ρR02[3κpG,02σR0].\omega_0^2 = \frac{1}{\rho R_0^2}\left[3 \kappa p_{G,0} - \frac{2\sigma}{R_0}\right].

The natural frequency of a small-amplitude oscillation about R0R_0 is

ω0=1R03κpG,02σ/R0ρ.\boxed{\omega_0 = \frac{1}{R_0} \sqrt{\frac{3 \kappa p_{G,0} - 2\sigma/R_0}{\rho}}.}

This is the Minnaert frequency (Marcel Minnaert, 1933, who derived it for sound emission from rising bubbles in brooks). For typical air-in-water conditions with pG,01p_{G,0} \approx 1 atm, surface tension negligible compared to gas pressure, and κ=γ=1.4\kappa = \gamma = 1.4 (adiabatic, which is appropriate for the high frequencies these resonances live at):

f0=ω02π12πR03γpatmρw3.26 Hz⋅mR0.f_0 = \frac{\omega_0}{2\pi} \approx \frac{1}{2\pi R_0}\sqrt{\frac{3 \gamma p_{atm}}{\rho_w}} \approx \frac{3.26\ \text{Hz·m}}{R_0}.

A 1 mm bubble resonates at 3.26 kHz. A 10 μm bubble (the size of medical ultrasound contrast agents) resonates at 326 kHz. A 100 nm bubble at 32.6 MHz. The Minnaert relation is the central pitch-vs-size relation in bubble acoustics and underlies essentially every applied-bubble technology.

The viscous damping in the linearised equation is small for water bubbles larger than a few μm. Two other damping channels — radiation of acoustic energy to infinity, and thermal damping inside the gas — typically dominate. Their full treatment is reserved for a later chapter on driven oscillating bubbles, not yet drafted.

Rayleigh inertial collapse (1917)

Take an empty cavity (pB=0p_B = 0) in an inviscid surface-tension-free liquid at uniform ambient pressure p>0p_\infty > 0. The Rayleigh–Plesset equation reduces to

ρ[RR¨+32R˙2]=p.\rho \left[R \ddot R + \frac{3}{2} \dot R^2\right] = -p_\infty.

Multiply by R˙\dot R and recognise both sides as time derivatives:

ρddt[12R3R˙2]=pR2R˙3/3=pddt[R3]/3.\rho \frac{d}{dt}\left[\frac{1}{2} R^3 \dot R^2\right] = -p_\infty \cdot R^2 \dot R \cdot 3 / 3 = -p_\infty \frac{d}{dt}\left[R^3\right] / 3.

Wait — let me be more careful. The left-hand side of the R–P equation, multiplied by R˙\dot R:

ρR˙[RR¨+32R˙2]=ρddt[12R2R˙2R]=ρddt[12RR2R˙2].\rho \dot R \left[R \ddot R + \frac{3}{2} \dot R^2 \right] = \rho \frac{d}{dt}\left[\frac{1}{2} R^2 \dot R^2 \cdot R\right] = \rho \frac{d}{dt}\left[\frac{1}{2} R \cdot R^2 \dot R^2\right].

The cleanest form is to recognise that ddt(R3R˙2)=3R2R˙R˙2+2R3R˙R¨=R2R˙(3R˙2+2RR¨)\frac{d}{dt}(R^3 \dot R^2) = 3 R^2 \dot R \cdot \dot R^2 + 2 R^3 \dot R \ddot R = R^2 \dot R (3 \dot R^2 + 2 R \ddot R), so

ρR˙[RR¨+(3/2)R˙2]=ρ2R2R˙2R21R2R˙ddt(R3R˙2R1/2)\rho \dot R [R \ddot R + (3/2)\dot R^2] = \frac{\rho}{2} R^2 \dot R \cdot \frac{2}{R^2} \cdot \frac{1}{R^2 \dot R} \cdot \frac{d}{dt}(R^3 \dot R^2 \cdot R^{-1}/2)

This is getting awkward. Let me proceed differently.

Rayleigh's first integral

Multiply the R–P equation ρ[RR¨+(3/2)R˙2]=p\rho[R\ddot R + (3/2)\dot R^2] = -p_\infty through by R2R˙R^2 \dot R:

ρR3R˙R¨+32ρR2R˙3=pR2R˙.\rho R^3 \dot R \ddot R + \frac{3}{2} \rho R^2 \dot R^3 = -p_\infty R^2 \dot R.

The left side is 12ρddt(R3R˙2)\frac{1}{2} \rho \frac{d}{dt}(R^3 \dot R^2):

ddt(R3R˙2)=3R2R˙R˙2+2R3R˙R¨=R2R˙(3R˙2+2RR¨),\frac{d}{dt}(R^3 \dot R^2) = 3 R^2 \dot R \dot R^2 + 2 R^3 \dot R \ddot R = R^2 \dot R(3\dot R^2 + 2 R \ddot R),

so

12ρddt(R3R˙2)=12ρR2R˙(3R˙2+2RR¨)=ρ(R3R˙R¨+32R2R˙3),\frac{1}{2} \rho \frac{d}{dt}(R^3 \dot R^2) = \frac{1}{2}\rho R^2 \dot R(3 \dot R^2 + 2 R \ddot R) = \rho(R^3 \dot R \ddot R + \tfrac{3}{2} R^2 \dot R^3),

which is exactly the left-hand side. The right-hand side is

pR2R˙=p3ddt(R3).-p_\infty R^2 \dot R = -\frac{p_\infty}{3}\frac{d}{dt}(R^3).

Integrating in time, with initial condition R(0)=RmaxR(0) = R_{\max}, R˙(0)=0\dot R(0) = 0:

12ρR3R˙2=p3(Rmax3R3),\frac{1}{2}\rho R^3 \dot R^2 = \frac{p_\infty}{3}\left(R_{\max}^3 - R^3\right),

so

R˙2=2p3ρ[(RmaxR)31].\dot R^2 = \frac{2 p_\infty}{3 \rho}\left[\left(\frac{R_{\max}}{R}\right)^3 - 1\right].

The collapse velocity R˙|\dot R| diverges as R0R \to 0: R˙R3/2\dot R \sim -R^{-3/2}. The collapse time can be obtained by integrating dR/R˙=dtdR / \dot R = -dt from R=RmaxR = R_{\max} to R=0R = 0:

τRayleigh=Rmax3ρ2p0.91470.915Rmaxρp.\tau_\text{Rayleigh} = R_{\max} \sqrt{\frac{3 \rho}{2 p_\infty}} \cdot 0.9147 \approx 0.915\, R_{\max} \sqrt{\frac{\rho}{p_\infty}}.

For a 1 mm bubble collapsing at p=1p_\infty = 1 atm in water, this is about 90 μs. For a 10 μm bubble at the same pressure, about 0.9 μs. The collapse is fast: characteristic acoustic timescales are not long compared to it, so the assumption of incompressible liquid breaks down. Real collapses radiate acoustic energy and emit shock waves at the moment of minimum radius — a regime that belongs to a later chapter on bubble collapse, not yet drafted.

Inertial growth — Plesset’s tension step

Now consider the opposite problem: a small bubble in a liquid suddenly subjected to a tension p<0p_\infty < 0 that exceeds the Blake threshold. At the instant the tension is applied, the bubble’s gas pressure plus vapour pressure exceeds the ambient pressure plus surface-tension term, and the bubble begins to grow.

In the limit where the gas-pressure term becomes negligible (the bubble has grown large enough that K/R30K/R^3 \to 0) and surface tension is also negligible, the R–P equation reduces to

ρ[RR¨+(3/2)R˙2]=p=p.\rho[R\ddot R + (3/2)\dot R^2] = -p_\infty = |p_\infty|.

The same first-integral analysis as above gives

R˙2=2p3ρ[1(R0R)3]2p3ρasRR0.\dot R^2 = \frac{2 |p_\infty|}{3 \rho}\left[1 - \left(\frac{R_0}{R}\right)^3\right] \to \frac{2 |p_\infty|}{3 \rho} \quad \text{as} \quad R \gg R_0.

The bubble grows at the asymptotic Rayleigh growth velocity

R˙=2p3ρ.\dot R_\infty = \sqrt{\frac{2 |p_\infty|}{3 \rho}}.

For p=1p_\infty = -1 atm in water, R˙8.2\dot R_\infty \approx 8.2 m/s — a 100 μm bubble grown in ~12 μs. Inertial cavitation growth is fast — the bubble explodes outward — but slower than inertial collapse, where wall speeds reach hundreds of m/s before the cavity vanishes.

Numerical integration of the full equation

For arbitrary drive p(t)p_\infty(t) and finite values of pBp_B, σ\sigma, μ\mu, no closed form exists. Standard practice is to recast the second-order ODE as two coupled first-order equations,

R˙=U,U˙=1R[pBp2σ/R4μU/Rρ32U2],\dot R = U, \qquad \dot U = \frac{1}{R}\left[\frac{p_B - p_\infty - 2\sigma/R - 4\mu U/R}{\rho} - \frac{3}{2} U^2\right],

and integrate with an adaptive Runge–Kutta scheme (RK4 with step control, RK45, or implicit methods for stiff regions during collapse).

R₀ = 10 μm0.085.2170.4255.6340.9bubble radius (μm)10p_∞ (atm)0.010.020.030.040.0time (μs)peak radius309.9 μm(31.0 × R₀)min radius10.00 μm(1.00 × R₀)gas pressureat R₀: 1.12 atmterms✓ surface tension✓ viscous damping
drive:

The Rayleigh-Plesset equation ρ[RR̈ + (3/2)Ṙ²] = (p_g + p_v) − p_∞(t) − 2σ/R − 4μṘ/R is a stiff nonlinear ODE in the bubble radius R(t). At equilibrium the gas pressure balances atmospheric plus surface tension; perturb the ambient pressure and the bubble responds. Step drives below −1 atm produce explosive growth (the bubble "cavitates"); sinusoidal drives below the Blake threshold produce stable oscillations and above it produce transient inertial collapse. Lithotripter pulses (positive spike + negative tail) drive a brief overpressure followed by intense tension, producing a characteristic growth-and-collapse signature that emits a shock wave at minimum radius. Toggle surface tension and viscosity to see their relative roles: surface tension dominates the equilibrium for small R; viscous damping is negligible for water bubbles above a few μm.

The solver above integrates the full Rayleigh–Plesset equation under three drive options:

What’s left to develop

The Rayleigh–Plesset equation, in its full nonlinear form, is the workhorse of the rest of the field. The natural next directions are:

In each of these, the Rayleigh–Plesset equation appears as the same compact second-order ODE we have just derived and solved here, but interpreted under different boundary conditions and approximations. The chapters that develop them are still to be written.