---
title: "Rayleigh-Plateau instability"
firstLetter: "R"
publishDate: 2021-03-10
categories:
- Physics
- Fluid mechanics
- Perturbation
- Surface tension

date: 2021-03-10T09:13:22+01:00
draft: false
markup: pandoc
---

# 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 $R_0$ and surface tension $\alpha$.
In this case, the [Young-Laplace equation](/know/concept/young-laplace-law/)
states that its internal pressure
is a constant $p_i$ expressed as follows,
where $p_o$ is the exterior air pressure:

$$\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 $z$-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_\epsilon$, assumed to be small,
to the internal pressure, which we allow to vary with time and space.
We use cylindrical coordinates:

$$\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 $\vec{u} = (u_r, u_\phi, u_z)$,
which obeys the following Euler equations:

$$\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 $\vec{u}$ is so small that we can ignore
the quadratic term in the former equation, leaving:

$$\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_\epsilon$:

$$\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:

$$\begin{aligned}
    \nabla^2 p_\epsilon
    = \pdv[2]{p_\epsilon}{r} + \frac{1}{r} \pdv{p_\epsilon}{r} + \pdv[2]{p_\epsilon}{z} + \frac{1}{r^2} \pdv[2]{p_\epsilon}{\phi}
    = 0
\end{aligned}$$

Finally, we add a perturbation $R_\epsilon \ll R_0$
to the radius of the surface of the liquid column:

$$\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_\epsilon$, $p_\epsilon$ and $\vec{u}$
do not depend on $\phi$.
The Laplace equation then reduces to:

$$\begin{aligned}
    \nabla^2 p_\epsilon
    = \pdv[2]{p_\epsilon}{r} + \frac{1}{r} \pdv{p_\epsilon}{r} + \pdv[2]{p_\epsilon}{z}
    = 0
\end{aligned}$$

Before solving this, we need boundary conditions.
The radial fluid velocity $u_r$ (the $r$-component of $\vec{u}$)
at the column surface $r\!=\!R$ is the
[material derivative](/know/concept/material-derivative/) of $R_\epsilon$:

$$\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_\epsilon$
varies slowly with respect to $z$:

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

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

$$\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 $R_1$ and $R_2$ are the principal curvature radii of the column surface.
These two expressions must be equivalent,
so, by inserting the definition of $p_i = p_o + \alpha / R_0$:

$$\begin{aligned}
    p_o + \alpha \Big( \frac{1}{R_1} + \frac{1}{R_2} \Big)
    %= p_i + p_\epsilon(r\!=\!R)
    = p_o + \frac{\alpha}{R_0} + p_\epsilon(r\!=\!R)
\end{aligned}$$

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

$$\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 $R_0 + R_\epsilon$,
while the curvature along the length can be approximated
using the second $z$-derivative of $R_\epsilon$:

$$\begin{aligned}
    p_\epsilon(r\!=\!R)
    \approx \alpha \Big( \frac{1}{R_0 + R_\epsilon} - \pdv[2]{R_\epsilon}{z} \Big) - \frac{\alpha}{R_0}
\end{aligned}$$

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

$$\begin{aligned}
    p_\epsilon(r\!=\!R)
    \approx - \alpha \Big( \frac{R_\epsilon}{R_0^2 + R_\epsilon} + \pdv[2]{R_\epsilon}{z} \Big)
    \approx - \alpha \Big( \frac{R_\epsilon}{R_0^2} + \pdv[2]{R_\epsilon}{z} \Big)
\end{aligned}$$

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

$$\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_\epsilon$ yields
Bessel's modified equation of order zero:

$$\begin{aligned}
    \dv[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 $I_0$ and $K_0$.
However, because $K_0$ diverges at $r = 0$, we must set the constant $B = 0$:

$$\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_\epsilon$
gives us the following relation:

$$\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 $\vec{u}$
states that $u_r$ is given by:

$$\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 $u_r$,
we can revisit its boundary condition:

$$\begin{aligned}
    u_r(r\!=\!R)
    = - \frac{A k}{\rho \sigma} I_0'(k R)
    = \sigma R_\epsilon
\end{aligned}$$

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

$$\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 $R \approx R_0$:

$$\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 $\sigma^2 > 0$.
Since $(1 - k^2 R_0^2)$ is the only factor that can be negative,
we need $k R_0 < 1$, leading us to the **critical wavelength** $\lambda_c$:

$$\begin{aligned}
    \boxed{
        \lambda_c
        = \frac{2 \pi}{k}
        = 2 \pi R_0
    }
\end{aligned}$$

If the perturbation wavelength $\lambda$ is larger than $\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 $\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.



## References
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.