---
title: "Stokes' law"
firstLetter: "S"
publishDate: 2021-05-04
categories:
- Physics
- Fluid mechanics
- Fluid dynamics

date: 2021-05-04T19:28:30+02:00
draft: false
markup: pandoc
---

# Stokes' law

**Stokes' law** describes the size of the drag force $D$
at low [Reynolds number](/know/concept/reynolds-number/) $\mathrm{Re} \ll 1$
experienced by a spherical object in a steady, uniform flow at velocity $U$.


## Flow field

Imagine a sphere with radius $a$ sinking in a viscous liquid.
To model this situation, let us pretend that the sphere is fixed instead,
and the fluid comes from infinity at velocity $U$ along the $z$-axis,
flows past the sphere, and continues to infinity at the same $U$.
The Reynolds number is:

$$\begin{aligned}
    \mathrm{Re}
    = \frac{2 a U}{\nu}
\end{aligned}$$

We assume that $\mathrm{Re} \ll 1$, in which case
the incompressible [Navier-Stokes equations](/know/concept/navier-stokes-equations/)
are reduced to the **steady Stokes equations**:

$$\begin{aligned}
    \nabla p
    = \eta \nabla^2 \va{v}
    \qquad \quad
    \nabla \cdot \va{v}
    = 0
\end{aligned}$$

The goal is to solve for $p$ and $\va{v}$.
We make the following ansatz in
[spherical coordinates](/know/concept/spherical-coordinates/) $(r, \theta, \phi)$,
where $q(r)$, $f(r)$ and $g(r)$ are unknown functions:

$$\begin{gathered}
    p
    = \eta U q(r) \cos\theta
    \\
    v_r
    = U f(r) \cos\theta
    \qquad
    v_\theta
    = - U g(r) \sin\theta
    \qquad
    v_\phi
    = 0
\end{gathered}$$

The fluid hits the sphere head on,
so the solution is taken to be $\phi$-independent due to symmetry.
Note that $\theta$ is the angle to the positive $z$-axis,
which is the direction of $\va{U} = U \vu{e}_z$.
Moreover, note that $\va{U} \cdot \vu{e}_r = U \cos\theta$
and $\va{U} \cdot \vu{e}_\theta = - U \sin\theta$,
where $\vu{e}_r$ and $\vu{e}_\theta$ are basis vectors.

To begin with, we insert this ansatz into the incompressibility condition,
yielding:

$$\begin{aligned}
    0
    = \nabla \cdot \va{v}
    &= \pdv{v_r}{r} + \frac{1}{r} \pdv{v_\theta}{\theta} + \frac{2 v_r}{r} + \frac{v_\theta}{r \tan \theta}
    \\
    &= U \dv{f}{r} \cos\theta - \frac{U g}{r} \cos\theta + \frac{2 U f}{r} \cos\theta - \frac{U g}{r} \cos\theta
    \\
    &= U \cos\theta \Big( \dv{f}{r} + \frac{2}{r} f - \frac{2}{r} g \Big)
\end{aligned}$$

The parenthesized expression must be zero for all $r$,
leading us to the following relation:

$$\begin{aligned}
    g(r)
    = f + \frac{r}{2} \dv{f}{r}
\end{aligned}$$

Next, we take the divergence of the first Stokes equation,
and insert incompressibility:

$$\begin{aligned}
    \nabla^2 p
    = \eta \nabla \cdot (\nabla^2 \va{v})
    = \eta \nabla^2 (\nabla \cdot \va{v})
    = 0
\end{aligned}$$

This is simply the Laplace equation,
which is as follows for our ansatz $p(r, \theta)$:

$$\begin{aligned}
    0
    = \nabla^2 p
    &= \frac{1}{r^2} \pdv{r} \Big( r^2 \pdv{p}{r} \Big) + \frac{1}{r^2 \sin\theta} \pdv{\theta} \Big( \sin\theta \pdv{p}{\theta} \Big)
    \\
    0
    &= \frac{\eta U \cos\theta}{r^2} \dv{r} \Big( r^2 \dv{q}{r} \Big)
    - \frac{\eta U q}{r^2 \sin\theta} \pdv{\theta} \Big( \sin^2\theta \Big)
    \\
    &= \frac{\eta U \cos\theta}{r^2} \dv{r} \Big( r^2 \dv{q}{r} \Big)
    - \frac{2 \eta U q}{r^2 \sin\theta} \sin\theta \cos\theta
    \\
    &= \eta U \cos\theta \Big( \dv[2]{q}{r} + \frac{2}{r} \dv{q}{r} - \frac{2}{r^2} q \Big)
\end{aligned}$$

Again, the parenthesized expression must be zero for all $r$,
meaning it is an ODE for $q(r)$,
whose solution is straightforwardly found to be:

$$\begin{aligned}
    q(r)
    = \frac{C_3}{r^2} + C_4 r
\end{aligned}$$

Where $C_3$ and $C_4$ are linearity constants ($C_1$ and $C_2$ appear later).
The pressure is therefore:

$$\begin{aligned}
    p
    = \eta U \cos\theta \Big( \frac{C_3}{r^2} + C_4 r \Big)
\end{aligned}$$

Consequently, its gradient $\nabla p$ in spherical coordinates is as follows:

$$\begin{aligned}
    \nabla p
    = \vu{e}_r \pdv{p}{r} + \vu{e}_\theta \frac{1}{r} \pdv{p}{\theta}
    = \vu{e}_r \Big( \eta U \cos\theta \dv{q}{r} \Big) - \vu{e}_\theta \Big( \eta U \sin\theta \frac{q}{r} \Big)
\end{aligned}$$

According to the Stokes equation, this equals $\eta \nabla^2 \va{v}$.
Let us look at the $r$-component of $\nabla^2 \va{v}$:

$$\begin{aligned}
    (\nabla^2 \va{v})_r
    &= \pdv[2]{v_r}{r} + \frac{1}{r^2} \pdv[2]{v_r}{\theta} + \frac{2}{r} \pdv{v_r}{r}
    + \frac{\cot\theta}{r^2} \pdv{v_r}{\theta} - \frac{2}{r^2} \pdv{v_\theta}{\theta} - \frac{2}{r^2} v_r - \frac{2 \cot\theta}{r^2} v_\theta
    \\
    &= U \cos\theta \Big( \dv[2]{f}{r} - \frac{1}{r^2} f + \frac{2}{r} \dv{f}{r} - \frac{1}{r^2} f
    + \frac{2}{r^2} g - \frac{2}{r^2} f + \frac{2}{r^2} g \Big)
    \\
    &= U \cos\theta \Big( \dv[2]{f}{r} + \frac{2}{r} \dv{f}{r} - \frac{4}{r^2} f + \frac{4}{r^2} g \Big)
\end{aligned}$$

Substituting $g$ for the expression we found from incompressibility lets us simplify this:

$$\begin{aligned}
    \eta (\nabla^2 \va{v})_r
    &= \eta U \cos\theta \Big( \dv[2]{f}{r} + \frac{4}{r} \dv{f}{r} \Big)
\end{aligned}$$

The Stokes equation says that this must be equal to the $r$-component of $\nabla p$:

$$\begin{aligned}
    \eta U \cos\theta \Big( \dv[2]{f}{r} + \frac{4}{r} \dv{f}{r} \Big)
    = \eta U \cos\theta \Big( \!-\! \frac{2 C_3}{r^3} + C_4 \Big)
\end{aligned}$$

Where we have inserted $\dv*{q}{r}$.
Dividing out $\eta U \cos\theta$ leaves an ODE for $f(r)$,
satisfied by:

$$\begin{aligned}
    f(r)
    = C_1 + \frac{C_2}{r^3} + \frac{C_3}{r} + \frac{C_4 r^2}{10}
\end{aligned}$$

Then, thanks to our earlier relation again,
we know that $g(r)$ is as follows:

$$\begin{aligned}
    g(r)
    = C_1 - \frac{C_2}{2 r^3} + \frac{C_3}{2 r} + \frac{C_4 r^2}{5}
\end{aligned}$$

So what about $C_1$, $C_2$, $C_3$ and $C_4$?
For $r\!\to\!\infty$, we expect that $\va{v}\!\to\!\va{U}$,
meaning that $f(r)\!\to\!1$ and $g(r)\!\to\!1$.
This implies that $C_4 = 0$ and $C_1 = 1$, leaving:

$$\begin{aligned}
    f(r)
    = 1 + \frac{C_2}{r^3} + \frac{C_3}{r}
    \qquad \quad
    g(r)
    = 1 - \frac{C_2}{2 r^3} + \frac{C_3}{2 r}
\end{aligned}$$

Furthermore, the viscous *no-slip* condition demands
that $\va{v} = 0$ at the sphere's surface $r = a$, so $f(a) = g(a) = 0$ there.
Inserting $a$ into $f$ and $g$, setting them to zero,
and solving the resulting system of equations
yields $C_2 = a^3 / 2$ and $C_3 = -3 a / 2$.
Therefore the full solution is:

$$\begin{gathered}
    \boxed{
        p
        = - \frac{3 \eta U a}{2 r^2} \cos\theta
    }
    \\
    \boxed{
        v_r
        = U \cos\theta \Big( 1 + \frac{a^3}{2 r^3} - \frac{3 a}{2 r} \Big)
        \qquad
        v_\theta
        = - U \sin\theta \Big( 1 - \frac{a^3}{4 r^3} - \frac{3 a}{4 r} \Big)
    }
\end{gathered}$$


## Drag force

From the definition of [viscosity](/know/concept/viscosity/),
we know that there must be shear stresses at the sphere surface,
described by the fluid's [Cauchy stress tensor](/know/concept/cauchy-stress-tensor/) $\hat{\sigma}$.
The drag force $\va{D}$ on the surface is:

$$\begin{aligned}
    \va{D}
    = \oint \hat{\sigma} \cdot \dd{\va{S}}
    = \int_0^{2\pi} \!\!\!\! \int_0^\pi \big( \hat{\sigma} \cdot \vu{e}_r \big) \:a^2 \sin\theta \dd{\theta} \dd{\phi}
\end{aligned}$$

Where $\vu{e}_r$ is the sphere's surface normal vector.
The integrand can be expanded as follows:

$$\begin{aligned}
    \hat{\sigma} \cdot \vu{e}_r
    = \vu{e}_r \sigma_{rr} + \vu{e}_\theta \sigma_{\theta r}
\end{aligned}$$

To calculate this, we start by taking the gradient of the velocity field $\va{v}$:

$$\begin{aligned}
    \nabla\va{v}
    &= \vu{e}_r \vu{e}_r \pdv{v_r}{r} + \vu{e}_r \vu{e}_\theta \pdv{v_\theta}{r}
    + \vu{e}_\theta \vu{e}_r \Big( \frac{1}{r} \pdv{v_r}{\theta} - \frac{v_\theta}{r} \Big)
    \\
    &\qquad + \vu{e}_\theta \vu{e}_\theta \Big( \frac{1}{r} \pdv{v_\theta}{\theta} - \frac{v_r}{r} \Big)
    + \vu{e}_\phi \vu{e}_\phi \Big( \frac{v_\theta}{r \tan\theta} + \frac{v_r}{r} \Big)
\end{aligned}$$

Some of these terms are necessary to calculate the stress elements $\sigma_{rr}$ and $\sigma_{\theta r}$:

$$\begin{aligned}
    \sigma_{rr}
    &= - p + 2 \eta (\nabla\va{v})_{rr}
    = - p + 2 \eta \pdv{v_r}{r}
    \\
    &= \frac{3 \eta U a}{2 r^2} \cos\theta + 2 \eta U \cos\theta \: \Big( \!-\! \frac{3 a^3}{2 r^4} + \frac{3 a}{2 r^2} \Big)
    \\
    &= \frac{3 \eta U a}{2 r^2} \cos\theta \: \Big( 3 - 2 \frac{a^2}{r^2} \Big)
\end{aligned}$$
$$\begin{aligned}
    \sigma_{\theta r}
    &= \eta \big( (\nabla\va{v})_{\theta r} + (\nabla\va{v})_{r \theta} \big)
    = \eta \: \Big( \pdv{v_\theta}{r} + \frac{1}{r} \pdv{v_r}{\theta} - \frac{v_\theta}{r} \Big)
    \\
    &= \eta U \sin\theta \: \Big( \!-\! \frac{3 a^3}{4 r^4} - \frac{3 a}{4 r^2}
    - \frac{1}{r} - \frac{a^3}{2 r^4} + \frac{3 a}{2 r^2}
    + \frac{1}{r} - \frac{a^3}{4 r^4} - \frac{3 a}{4 r^2} \Big)
    \\
    &= - \frac{3 \eta U a^3}{2 r^4} \sin\theta
\end{aligned}$$

At the sphere's surface we set $r = a$, so these expressions reduce to the following:

$$\begin{aligned}
    \sigma_{rr}
    = \frac{3 \eta U}{2 a} \cos\theta
    \qquad \quad
    \sigma_{\theta r}
    = - \frac{3 \eta U}{2 a} \sin\theta
\end{aligned}$$

Now we can finally calculate the effective stress on the surface,
by converting the basis vectors $\vu{e}_r$ and $\vu{e}_\theta$ to Cartesian coordinates:

$$\begin{aligned}
    \hat{\sigma} \cdot \vu{e}_r
    &= \vu{e}_r \frac{3 \eta U}{2 a} \cos\theta - \vu{e}_\theta \frac{3 \eta U}{2 a} \sin\theta
    \\
    &= \Big( \vu{e}_x \sin\theta \cos\phi + \vu{e}_y \sin\theta \sin\phi + \vu{e}_z \cos\theta \Big) \frac{3 \eta U}{2 a} \cos\theta
    \\
    &\qquad - \Big( \vu{e}_x \cos\theta \cos\phi + \vu{e}_y \cos\theta \sin\phi - \vu{e}_z \sin\theta \Big) \frac{3 \eta U}{2 a} \sin\theta
    \\
    &= \Big( \vu{e}_z \cos^2\theta + \vu{e}_z \sin^2\theta \Big) \frac{3 \eta U}{2 a}
    = \vu{e}_z \frac{3 \eta U}{2 a}
\end{aligned}$$

Remarkably, the stress at every point on the sphere is purely in the $z$-direction!
This is not entirely unexpected though: symmetry cancels out all other components.

With this, we can do the integrals for $\va{D}$,
which reduce to a surface area factor $4 \pi a^2$:

$$\begin{aligned}
    \va{D}
    = \vu{e}_z \frac{3 \eta U}{2 a} \int_0^{2\pi} \!\!\!\! \int_0^\pi a^2 \sin\theta \dd{\theta} \dd{\phi}
    = \vu{e}_z \frac{3 \eta U}{2 a} 2 \pi a^2 \int_0^\pi \sin\theta \dd{\theta}
    = \vu{e}_z \: 6 \pi \eta U a
\end{aligned}$$

At last, we arrive at Stokes' law,
which simply expresses the magnitude of $\va{D}$:

$$\begin{aligned}
    \boxed{
        D
        = 6 \pi \eta U a
    }
\end{aligned}$$

To arrive at this result,
we assumed that the sphere was fixed, and the fluid was flowing past it.
We can equally well let the fluid be at rest,
with the sphere falling through it at $U$.
The force of gravity then exerts the following force $G$ on it,
subtracting [buoyancy](/know/concept/archimedes-principle/):

$$\begin{aligned}
    G
    = \frac{4 \pi a^3}{3} (\rho_s - \rho_f) g_0
\end{aligned}$$

Where $\rho_s$ and $\rho_f$ are the sphere's and fluid's densities,
and $g_0$ is the gravitational acceleration.
Since $D$ acts in the opposite sense of $G$,
after some time, they cancel out:

$$\begin{aligned}
    6 \pi \eta U a
    = \frac{4 \pi a^3}{3} (\rho_s - \rho_f) g_0
\end{aligned}$$

This is an equation for the **terminal velocity** $U_t$,
which we find to be as follows:

$$\begin{aligned}
    \boxed{
        U_t
        = \frac{2 a^2 (\rho_s - \rho_f) g_0}{9 \eta}
    }
\end{aligned}$$

The falling sphere will accelerate until $U_t$,
and then continue falling at constant speed.



## References
1.  B. Lautrup,
    *Physics of continuous matter: exotic and everyday phenomena in the macroscopic world*, 2nd edition,
    CRC Press.