Categories: Fluid dynamics, Fluid mechanics, Physics.

Stokes’ law

Stokes’ law describes the size of the drag force DD at low Reynolds number Re1\mathrm{Re} \ll 1 experienced by a spherical object in a steady, uniform flow at velocity UU.

Flow field

Imagine a sphere with radius aa 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 UU along the zz-axis, flows past the sphere, and continues to infinity at the same UU. The Reynolds number is:

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

We assume that Re1\mathrm{Re} \ll 1, in which case the incompressible Navier-Stokes equations are reduced to the steady Stokes equations:

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

The goal is to solve for pp and v\va{v}. We make the following ansatz in spherical coordinates (r,θ,ϕ)(r, \theta, \phi), where q(r)q(r), f(r)f(r) and g(r)g(r) are unknown functions:

p=ηUq(r)cosθvr=Uf(r)cosθvθ=Ug(r)sinθvϕ=0\begin{gathered} p = \eta U q(r) \cos\theta \\ v_r = U f(r) \cos\theta \qquad \quad v_\theta = - U g(r) \sin\theta \qquad \quad 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 zz-axis, which is the direction of U=Ue^z\va{U} = U \vu{e}_z. Moreover, note that Ue^r=Ucosθ\va{U} \cdot \vu{e}_r = U \cos\theta and Ue^θ=Usinθ\va{U} \cdot \vu{e}_\theta = - U \sin\theta, where e^r\vu{e}_r and e^θ\vu{e}_\theta are basis vectors.

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

0=v=vrr+1rvθθ+2vrr+vθrtanθ=UdfdrcosθUgrcosθ+2UfrcosθUgrcosθ=U(dfdr+2rf2rg)cosθ\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 \Big( \dv{f}{r} + \frac{2}{r} f - \frac{2}{r} g \Big) \cos\theta \end{aligned}

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

g(r)=f+r2dfdr\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:

2p=η(2v)=η2(v)=0\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,θ)p(r, \theta):

0=2p=1r2r(r2pr)+1r2sinθθ(pθsinθ)0=ηUcosθr2ddr(r2dqdr)ηUqr2sinθθ(sin2θ)=ηUr2ddr(r2dqdr)cosθ2ηUqr2sinθsinθcosθ=ηU(d2qdr2+2rdqdr2r2q)cosθ\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( \pdv{p}{\theta} \sin\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}{r^2} \dv{}{r}\Big( r^2 \dv{q}{r} \Big) \cos\theta - \frac{2 \eta U q}{r^2 \sin\theta} \sin\theta \cos\theta \\ &= \eta U \Big( \dvn{2}{q}{r} + \frac{2}{r} \dv{q}{r} - \frac{2}{r^2} q \Big) \cos\theta \end{aligned}

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

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

Where C3C_3 and C4C_4 are linearity constants (C1C_1 and C2C_2 appear later). The pressure is therefore:

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

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

p=e^rpr+e^θ1rpθ=e^r(ηUdqdrcosθ)e^θ(ηUqrsinθ)\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 \dv{q}{r} \cos\theta \Big) - \vu{e}_\theta \Big( \eta U \frac{q}{r} \sin\theta \Big) \end{aligned}

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

(2v)r=2vrr2+1r22vrθ2+2rvrr+cotθr2vrθ2r2vθθ2r2vr2cotθr2vθ=U(d2fdr21r2f+2rdfdr1r2f+2r2g2r2f+2r2g)cosθ=U(d2fdr2+2rdfdr4r2f+4r2g)cosθ\begin{aligned} (\nabla^2 \va{v})_r &= \pdvn{2}{v_r}{r} + \frac{1}{r^2} \pdvn{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 \Big( \dvn{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) \cos\theta \\ &= U \Big( \dvn{2}{f}{r} + \frac{2}{r} \dv{f}{r} - \frac{4}{r^2} f + \frac{4}{r^2} g \Big) \cos\theta \end{aligned}

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

η(2v)r=ηU(d2fdr2+4rdfdr)cosθ\begin{aligned} \eta (\nabla^2 \va{v})_r &= \eta U \Big( \dvn{2}{f}{r} + \frac{4}{r} \dv{f}{r} \Big) \cos\theta \end{aligned}

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

ηU(d2fdr2+4rdfdr)cosθ=ηU( ⁣ ⁣2C3r3+C4)cosθ\begin{aligned} \eta U \Big( \dvn{2}{f}{r} + \frac{4}{r} \dv{f}{r} \Big) \cos\theta = \eta U \Big( \!-\! \frac{2 C_3}{r^3} + C_4 \Big) \cos\theta \end{aligned}

Where we have inserted dq/dr\idv{q}{r}. Dividing out ηUcosθ\eta U \cos\theta leaves an ODE for f(r)f(r), satisfied by:

f(r)=C1+C2r3+C3r+C4r210\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)g(r) is as follows:

g(r)=C1C22r3+C32r+C4r25\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 C1C_1, C2C_2, C3C_3 and C4C_4? For r ⁣ ⁣r\!\to\!\infty, we expect that v ⁣ ⁣U\va{v}\!\to\!\va{U}, meaning that f(r) ⁣ ⁣1f(r)\!\to\!1 and g(r) ⁣ ⁣1g(r)\!\to\!1. This implies that C4=0C_4 = 0 and C1=1C_1 = 1, leaving:

f(r)=1+C2r3+C3rg(r)=1C22r3+C32r\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 v=0\va{v} = 0 at the sphere’s surface r=ar = a, so f(a)=g(a)=0f(a) = g(a) = 0 there. Inserting aa into ff and gg, setting them to zero, and solving the resulting system of equations yields C2=a3/2C_2 = a^3 / 2 and C3=3a/2C_3 = -3 a / 2. Therefore the full solution is:

p=3ηUa2r2cosθvr=U(1+a32r33a2r)cosθvθ=U(1a34r33a4r)sinθ\begin{gathered} \boxed{ p = - \frac{3 \eta U a}{2 r^2} \cos\theta } \\ \boxed{ v_r = U \Big( 1 + \frac{a^3}{2 r^3} - \frac{3 a}{2 r} \Big) \cos\theta \qquad v_\theta = - U \Big( 1 - \frac{a^3}{4 r^3} - \frac{3 a}{4 r} \Big) \sin\theta } \end{gathered}

Drag force

From the definition of viscosity, we know that there must be shear stresses at the sphere surface, described by the fluid’s Cauchy stress tensor σ^\hat{\sigma}. The drag force D\va{D} on the surface is:

D=σ^dS=02π ⁣ ⁣ ⁣ ⁣0π(σ^e^r)a2sinθdθdϕ\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 e^r\vu{e}_r is the sphere’s surface normal vector. The integrand can be expanded as follows:

σ^e^r=e^rσrr+e^θσθr\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 v\va{v}:

v=e^re^rvrr+e^re^θvθr+e^θe^r(1rvrθvθr)+e^θe^θ(1rvθθvrr)+e^ϕe^ϕ(vθrtanθ+vrr)\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 σrr\sigma_{rr} and σθr\sigma_{\theta r}:

σrr=p+2η(v)rr=p+2ηvrr=3ηUa2r2cosθ+2ηUcosθ( ⁣ ⁣3a32r4+3a2r2)=3ηUa2r2cosθ(32a2r2)\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} σθr=η((v)θr+(v)rθ)=η(vθr+1rvrθvθr)=ηUsinθ( ⁣ ⁣3a34r43a4r21ra32r4+3a2r2+1ra34r43a4r2)=3ηUa32r4sinθ\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=ar = a, so these expressions reduce to the following:

σrr=3ηU2acosθσθr=3ηU2asinθ\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 e^r\vu{e}_r and e^θ\vu{e}_\theta to Cartesian coordinates:

σ^e^r=e^r3ηU2acosθe^θ3ηU2asinθ=(e^xsinθcosϕ+e^ysinθsinϕ+e^zcosθ)3ηU2acosθ(e^xcosθcosϕ+e^ycosθsinϕe^zsinθ)3ηU2asinθ=(e^zcos2θ+e^zsin2θ)3ηU2a=e^z3ηU2a\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 zz-direction! This is not entirely unexpected though: symmetry cancels out all other components.

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

D=e^z3ηU2a02π ⁣ ⁣ ⁣ ⁣0πa2sinθdθdϕ=e^z3ηU2a2πa20πsinθdθ=e^z6πηUa\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 D\va{D}:

D=6πηUa\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 UU. The force of gravity then exerts the following force GG on it, subtracting buoyancy:

G=4πa33(ρsρf)g0\begin{aligned} G = \frac{4 \pi a^3}{3} (\rho_s - \rho_f) g_0 \end{aligned}

Where ρs\rho_s and ρf\rho_f are the sphere’s and fluid’s densities, and g0g_0 is the gravitational acceleration. Since DD acts in the opposite sense of GG, after some time, they cancel out:

6πηUa=4πa33(ρsρf)g0\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 UtU_t, which we find to be as follows:

Ut=2a2(ρsρf)g09η\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 UtU_t, and then continue falling at constant speed.


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