Categories: Fluid mechanics, Perturbation, Physics, Surface tension.

Rayleigh-Plateau instability

In fluid mechanics, the Rayleigh-Plateau instability causes a column of liquid to break up due to surface tension. It is the reason why a smooth stream of water (e.g. from a tap) eventually breaks into droplets as it falls.

Consider an infinitely long cylinder of liquid with radius R0R_0 and surface tension α\alpha. In this case, the Young-Laplace equation states that its internal pressure is a constant pip_i expressed as follows, where pop_o is the exterior air pressure:

pi=po+αR0\begin{aligned} p_i = p_o + \frac{\alpha}{R_0} \end{aligned}

We assume that the liquid is at rest. Alternatively, if it is moving in the zz-direction, we can also let our coordinate system travel at the same speed. Anyway, for convenience, we neglect any motion or acceleration of the liquid column.

Next, we add a perturbation pϵp_\epsilon, assumed to be small, to the internal pressure, which we allow to vary with time and space. We use cylindrical coordinates:

p(r,ϕ,z,t)=pi+pϵ(r,ϕ,z,t)\begin{aligned} p(r, \phi, z, t) = p_i + p_\epsilon(r, \phi, z, t) \end{aligned}

This internal pressure difference will cause the liquid to start to flow. We express the flow velocity as a vector u=(ur,uϕ,uz)\vec{u} = (u_r, u_\phi, u_z), which obeys the following Euler equations:

ut+(u)u=1ρpu=0\begin{aligned} \pdv{\vec{u}}{t} + (\vec{u} \cdot \nabla) \vec{u} = - \frac{1}{\rho} \nabla p \qquad \qquad \nabla \cdot \vec{u} = 0 \end{aligned}

The latter equation states that the fluid is incompressible. We assume that u\vec{u} is so small that we can ignore the quadratic term in the former equation, leaving:

ut=1ρpϵ\begin{aligned} \pdv{\vec{u}}{t} = - \frac{1}{\rho} \nabla p_\epsilon \end{aligned}

Taking the divergence and using incompressibility yields the Laplace equation for pϵp_\epsilon:

1ρ2pϵ=t(u)=0    2pϵ=0\begin{aligned} - \frac{1}{\rho} \nabla^2 p_\epsilon = \pdv{}{t}(\nabla \cdot \vec{u}) = 0 \qquad \implies \qquad \nabla^2 p_\epsilon = 0 \end{aligned}

We write out the Laplacian in cylindrical coordinates to get the following problem:

2pϵ=2pϵr2+1rpϵr+2pϵz2+1r22pϵϕ2=0\begin{aligned} \nabla^2 p_\epsilon = \pdvn{2}{p_\epsilon}{r} + \frac{1}{r} \pdv{p_\epsilon}{r} + \pdvn{2}{p_\epsilon}{z} + \frac{1}{r^2} \pdvn{2}{p_\epsilon}{\phi} = 0 \end{aligned}

Finally, we add a perturbation RϵR0R_\epsilon \ll R_0 to the radius of the surface of the liquid column:

R(z,t)=R0+Rϵ(z,t)\begin{aligned} R(z, t) = R_0 + R_\epsilon(z, t) \end{aligned}

Note that there is no dependence on the angle ϕ\phi; the deformation is assumed to be symmetric. Imagine the cross-section of the cylinder, and convince yourself that all asymmetric deformations will be removed by surface tension, which prefers a circular shape. We thus assume that RϵR_\epsilon, pϵp_\epsilon and u\vec{u} do not depend on ϕ\phi. The Laplace equation then reduces to:

2pϵ=2pϵr2+1rpϵr+2pϵz2=0\begin{aligned} \nabla^2 p_\epsilon = \pdvn{2}{p_\epsilon}{r} + \frac{1}{r} \pdv{p_\epsilon}{r} + \pdvn{2}{p_\epsilon}{z} = 0 \end{aligned}

Before solving this, we need boundary conditions. The radial fluid velocity uru_r (the rr-component of u\vec{u}) at the column surface r ⁣= ⁣Rr\!=\!R is the material derivative of RϵR_\epsilon:

ur(r ⁣= ⁣R)=DRϵDt=Rϵt+uz(r ⁣= ⁣R)Rϵz\begin{aligned} u_r(r\!=\!R) = \frac{\mathrm{D} R_\epsilon}{\mathrm{D} t} = \pdv{R_\epsilon}{t} + u_z(r\!=\!R) \pdv{R_\epsilon}{z} \end{aligned}

We linearize this by assuming that the deformation RϵR_\epsilon varies slowly with respect to zz:

ur(r ⁣= ⁣R)Rϵt\begin{aligned} u_r(r\!=\!R) \approx \pdv{R_\epsilon}{t} \end{aligned}

Meanwhile, we can write the boundary condition of the pressure pp in two ways, respectively from the Young-Laplace equation and the definition of the perturbation pϵp_\epsilon:

p(r ⁣= ⁣R)=po+α(1R1+1R2)p(r ⁣= ⁣R)=pi+pϵ(r ⁣= ⁣R)\begin{aligned} p(r\!=\!R) = p_o + \alpha \Big( \frac{1}{R_1} + \frac{1}{R_2} \Big) \qquad \quad p(r\!=\!R) = p_i + p_\epsilon(r\!=\!R) \end{aligned}

Where R1R_1 and R2R_2 are the principal curvature radii of the column surface. These two expressions must be equivalent, so, by inserting the definition of pi=po+α/R0p_i = p_o + \alpha / R_0:

po+α(1R1+1R2)=po+αR0+pϵ(r ⁣= ⁣R)\begin{aligned} p_o + \alpha \Big( \frac{1}{R_1} + \frac{1}{R_2} \Big) = p_o + \frac{\alpha}{R_0} + p_\epsilon(r\!=\!R) \end{aligned}

Isolating this equation for pϵp_\epsilon yields the desired boundary condition:

pϵ(r ⁣= ⁣R)=α(1R1+1R2)αR0\begin{aligned} p_\epsilon(r\!=\!R) = \alpha \Big( \frac{1}{R_1} + \frac{1}{R_2} \Big) - \frac{\alpha}{R_0} \end{aligned}

The principal radius around the circumference is R0+RϵR_0 + R_\epsilon, while the curvature along the length can be approximated using the second zz-derivative of RϵR_\epsilon:

pϵ(r ⁣= ⁣R)α(1R0+Rϵ2Rϵz2)αR0\begin{aligned} p_\epsilon(r\!=\!R) \approx \alpha \Big( \frac{1}{R_0 + R_\epsilon} - \pdvn{2}{R_\epsilon}{z} \Big) - \frac{\alpha}{R_0} \end{aligned}

This can be simplified a bit by using the assumption that RϵR_\epsilon is small:

pϵ(r ⁣= ⁣R)α(RϵR02+Rϵ+2Rϵz2)α(RϵR02+2Rϵz2)\begin{aligned} p_\epsilon(r\!=\!R) \approx - \alpha \Big( \frac{R_\epsilon}{R_0^2 + R_\epsilon} + \pdvn{2}{R_\epsilon}{z} \Big) \approx - \alpha \Big( \frac{R_\epsilon}{R_0^2} + \pdvn{2}{R_\epsilon}{z} \Big) \end{aligned}

At last, we have all the necessary boundary condition. We now make the following ansatz, where kk is the wavenumber and σ\sigma describes exponential growth or decay:

u(r,z,t)=u(r)exp(σt)cos(kz)pϵ(r,z,t)=pϵ(r)exp(σt)cos(kz)Rϵ(z,t)=Rϵexp(σt)cos(kz)\begin{aligned} \vec{u}(r, z, t) &= \vec{u}(r) \exp(\sigma t) \cos(k z) \\ p_\epsilon(r, z, t) &= p_\epsilon(r) \exp(\sigma t) \cos(k z) \\ R_\epsilon(z, t) &= R_\epsilon \exp(\sigma t) \cos(k z) \end{aligned}

This is justified by the fact that we can Fourier-expand any perturbation; this ansatz is simply the dominant term of the resulting series.

Inserting this into the Laplace equation for pϵp_\epsilon yields Bessel’s modified equation of order zero:

d2pϵdr2+1rdpϵdrk2pϵ=0\begin{aligned} \dvn{2}{p_\epsilon}{r} + \frac{1}{r} \dv{p_\epsilon}{r} - k^2 p_\epsilon = 0 \end{aligned}

This has well-known solutions: the modified Bessel functions I0I_0 and K0K_0. However, because K0K_0 diverges at r=0r = 0, we must set the constant B=0B = 0:

pϵ(r)=AI0(kr)+BK0(kr)=AI0(kr)\begin{aligned} p_\epsilon(r) = A I_0(kr) + B K_0(kr) = A I_0(kr) \end{aligned}

Inserting the ansatz into the boundary condition for pϵp_\epsilon gives us the following relation:

pϵ(r ⁣= ⁣R)=αRϵ(1R02+k2)=AI0(kR)\begin{aligned} p_\epsilon(r\!=\!R) = - \alpha R_\epsilon \Big( \frac{1}{R_0^2} + k^2 \Big) = A I_0(k R) \end{aligned}

Meanwhile, the linearized Euler equation governing u\vec{u} states that uru_r is given by:

σur=1ρdpϵdr=AkρI0(kr)\begin{aligned} \sigma u_r = - \frac{1}{\rho} \dv{p_\epsilon}{r} = - \frac{A k}{\rho} I_0'(kr) \end{aligned}

Now that we have an expression for uru_r, we can revisit its boundary condition:

ur(r ⁣= ⁣R)=AkρσI0(kR)=σRϵ\begin{aligned} u_r(r\!=\!R) = - \frac{A k}{\rho \sigma} I_0'(k R) = \sigma R_\epsilon \end{aligned}

Isolating this for RϵR_\epsilon and inserting it into the boundary condition for pϵp_\epsilon yields:

pϵ(r ⁣= ⁣R)=AI0(kR)=α(1R02+k2)(Akρσ2I0(kR))\begin{aligned} p_\epsilon(r\!=\!R) = A I_0(kR) = \alpha \Big( \frac{1}{R_0^2} + k^2 \Big) \Big( \frac{A k}{\rho \sigma^2} I_0'(k R) \Big) \end{aligned}

Isolating this for the exponential growth/decay parameter σ\sigma gives us the desired result, where we have also used the fact that RR0R \approx R_0:

σ2=αkρR02(1k2R02)I0(kR0)I0(kR0)\begin{aligned} \sigma^2 = \frac{\alpha k}{\rho R_0^2} (1 - k^2 R_0^2) \frac{I_0'(kR_0)}{I_0(kR_0)} \end{aligned}

To get exponential growth (i.e. instability), we need σ2>0\sigma^2 > 0. Since (1k2R02)(1 - k^2 R_0^2) is the only factor that can be negative, we need kR0<1k R_0 < 1, leading us to the critical wavelength λc\lambda_c:

λc=2πk=2πR0\begin{aligned} \boxed{ \lambda_c = \frac{2 \pi}{k} = 2 \pi R_0 } \end{aligned}

If the perturbation wavelength λ\lambda is larger than λc\lambda_c, surface tension creates a higher pressure in the narrower sections compared to the wider ones, thereby pumping the liquid into the bulges, further increasing their size until they become droplets.

Else, if λ<λc\lambda < \lambda_c, the tighter curvatures dominate the action of surface tension, which will then try to smoothen the surface by shrinking the bulges and widening the constrictions. In other words, the liquid column is stable in this case.


  1. B. Lautrup, Physics of continuous matter: exotic and everyday phenomena in the macroscopic world, 2nd edition, CRC Press.
  2. T. Bohr, A. Andersen, The Rayleigh-Plateau instability of a liquid column, 2020, unpublished.