Categories: Fluid dynamics, Fluid mechanics, Physics.

Stokes’ law

Stokes’ law describes the size of the drag force DD at low Reynolds number Re≪1\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 Re≪1\mathrm{Re} \ll 1, in which case the incompressible Navier-Stokes equations are reduced to the steady Stokes equations:

∇p=η∇2v⃗∇⋅v⃗=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 U⃗⋅e^r=Ucos⁡θ\va{U} \cdot \vu{e}_r = U \cos\theta and U⃗⋅e^θ=−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⃗=∂vr∂r+1r∂vθ∂θ+2vrr+vθrtan⁡θ=Udfdrcos⁡θ−Ugrcos⁡θ+2Ufrcos⁡θ−Ugrcos⁡θ=U(dfdr+2rf−2rg)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=1r2∂∂r(r2∂p∂r)+1r2sin⁡θ∂∂θ(∂p∂θsin⁡θ)0=ηUcos⁡θr2ddr(r2dqdr)−ηUqr2sin⁡θ∂∂θ(sin⁡2θ)=ηUr2ddr(r2dqdr)cos⁡θ−2ηUqr2sin⁡θsin⁡θcos⁡θ=ηU(d2qdr2+2rdqdr−2r2q)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^r∂p∂r+e^θ1r∂p∂θ=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=∂2vr∂r2+1r2∂2vr∂θ2+2r∂vr∂r+cot⁡θr2∂vr∂θ−2r2∂vθ∂θ−2r2vr−2cot⁡θr2vθ=U(d2fdr2−1r2f+2rdfdr−1r2f+2r2g−2r2f+2r2g)cos⁡θ=U(d2fdr2+2rdfdr−4r2f+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)=C1−C22r3+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)=1−C22r3+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+a32r3−3a2r)cos⁡θvθ=−U(1−a34r3−3a4r)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^r∂vr∂r+e^re^θ∂vθ∂r+e^θe^r(1r∂vr∂θ−vθr)+e^θe^θ(1r∂vθ∂θ−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η∂vr∂r=3ηUa2r2cos⁡θ+2ηUcos⁡θ ( ⁣− ⁣3a32r4+3a2r2)=3ηUa2r2cos⁡θ (3−2a2r2)\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+1r∂vr∂θ−vθr)=ηUsin⁡θ ( ⁣− ⁣3a34r4−3a4r2−1r−a32r4+3a2r2+1r−a34r4−3a4r2)=−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^zcos⁡2θ+e^zsin⁡2θ)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ηU2a∫02π ⁣ ⁣ ⁣ ⁣∫0πa2sin⁡θdθdϕ=e^z3ηU2a2πa2∫0πsin⁡θdθ=e^z 6πη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.

References

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