An object falling slowly in a viscous fluid under the influence of gravity is subject to a drag force (Stokes drag) which varies linearly with its velocity. Its equation of motion may be written as the second order differential equation: $$ m\frac{\mathrm{d}^2z}{\mathrm{d}t^2} = -c \frac{\mathrm{d}z}{\mathrm{d}t} + mg', $$ where $z$ is the object's position as a function of time, $t$, $c$ is a drag constant which depends on the shape of the object and the fluid viscosity, and $$ g' = g\left(1-\frac{\rho_\mathrm{fluid}}{\rho_\mathrm{obj}}\right) $$ is the effective gravitational acceleration, which accounts for the buoyant force due to the fluid (density $\rho_\mathrm{fluid}$) displaced by the object (density $\rho_\mathrm{obj}$). For a small sphere of radius $r$ in a fluid of viscosity $\eta$, Stokes' law predicts $c = 6\pi\eta r$.
Consider a sphere of platinum ($\rho = 21.45\;\mathrm{g\,cm^{-3}}$) with radius 1 mm, initially at rest, falling in mercury ($\rho = 13.53\;\mathrm{g\,cm^{-3}}$, $\eta = 1.53 \times 10^{-3}\;\mathrm{Pa\,s}$). The above second-order differential equation can be solved analytically, but to integrate it numerically using odeint, it must be treated as two first-order ordinary differential equations:
\begin{align*}
\frac{\mathrm{d}z}{\mathrm{d}t} &= \dot{z}\\
\frac{\mathrm{d}^2z}{\mathrm{d}t^2} = \frac{\mathrm{d}\dot{z}}{\mathrm{d}t} &= g' - \frac{c}{m}\dot{z}
\end{align*}
In the code below, the function deriv calculates these derivatives and is passed to odeint with the intial conditions ($z=0$, $\dot{z}=0$) and a grid of time points.
# eg8-stokes-drag.py
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# Platinum sphere falling from rest in mercury.
# Acceleration due to gravity (m.s-2).
g = 9.81
# Densities (kg.m-3).
rho_Pt, rho_Hg = 21450, 13530
# Viscosity of mercury (Pa.s).
eta = 1.53e-3
# Radius and mass of the sphere.
r = 1.e-3 # radius (m)
m = 4*np.pi/3 * r**3 * rho_Pt
# Drag constant from Stokes' law.
c = 6 * np.pi * eta * r
# Effective gravitational acceleration.
gp = g * (1 - rho_Hg/rho_Pt)
def deriv(t, z):
""" Return the dz/dt and d2z/dt2. """
dz0 = z[1]
dz1 = gp - c/m * z[1]
return dz0, dz1
t0, tf = 0, 20
t_eval = np.linspace(t0, tf, 50)
# Initial conditions: z = 0, dz/dt = 0 at t = 0.
z0 = (0, 0)
# Integrate the pair of differential equations.
sol = solve_ivp(deriv, (t0, tf), z0, t_eval=t_eval)
t = sol.t
z, zdot = sol.y
plt.plot(t, zdot)
print('Estimate of terminal velocity = {:.3f} m.s-1'.format(zdot[-1]))
# Exact solution: terminal velocity vt (m.s-1) and characteristic time tau (s).
v0, vt, tau = 0, m*gp/c, m/c
print('Exact terminal velocity = {:.3f} m.s-1'.format(vt))
z = vt*t + v0*tau*(1-np.exp(-t/tau)) + vt*tau*(np.exp(-t/tau)-1)
zdot_exact = vt + (v0-vt)*np.exp(-t/tau)
plt.plot(t, zdot_exact)
plt.xlabel('$t$ /s')
plt.ylabel('$\dot{z}\;/\mathrm{m\, s^{-1}}$')
plt.show()
The plot produced by this program is shown below: the numerical and analytical results are indistinguishable at this scale but are reported to three decimal places in the output:
Estimate of terminal velocity = 11.266 m.s-1
Exact terminal velocity = 11.285 m.s-1
