Coulomb’s torsion (1784)

Charles-Augustin Coulomb (1784) hung a weight on a thin wire, twisted it and timed the oscillations. The restoring torque was proportional to the angle of twist, to the fourth power of the wire’s diameter and inversely to its length, which in modern form is \(\phi = TL/(GJ)\) with \(J = \pi d^4/32\) (torsion_angle()). The same torsion balance later weighed the electric force and, in Cavendish’s hands, the Earth. The example simulates Coulomb’s pendulum and recovers the shear modulus from the period, then shows the shear stress across the wire (torsion_shear_stress()).

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

from engineeringkit.structural_engineering import torsion_angle, torsion_shear_stress

A torsion pendulum

L, d, G = 0.8, 0.5e-3, 40e9  # brass wire
J = np.pi * d**4 / 32
k = 1.0 / torsion_angle(1.0, L, G, J)  # torque per radian [N m/rad]
I_disk = 0.5 * 0.2 * 0.05**2  # 200 g disk of radius 50 mm [kg m^2]
sol = solve_ivp(lambda t, s: [s[1], -k * s[0] / I_disk], (0, 120), [0.5, 0.0], max_step=0.05, rtol=1e-9, atol=1e-12)
crossings = sol.t[1:][np.diff(np.sign(sol.y[0])) > 0]
period = float(np.mean(np.diff(crossings)))
G_measured = 4 * np.pi**2 * I_disk * L / (period**2 * J)
print(f"period {period:.2f} s, recovered G = {G_measured / 1e9:.2f} GPa")
assert np.isclose(G_measured, G, rtol=1e-3)

period 5.67 s, recovered G = 40.00 GPa

Coulomb’s fourth-power law

diameters = np.array([0.25e-3, 0.5e-3, 1e-3])
torques = 0.1 / torsion_angle(1.0, L, G, np.pi * diameters**4 / 32)  # torque for 0.1 rad
print("torque ratios for doubling the diameter:", np.round(torques[1:] / torques[:-1], 3))
assert np.allclose(torques[1:] / torques[:-1], 16)

r = np.linspace(0, d / 2, 50)
fig, axes = plt.subplots(1, 2, figsize=(10, 3.5))
axes[0].plot(sol.t, sol.y[0])
axes[0].set(xlabel="time [s]", ylabel="twist [rad]", title="Torsion pendulum")
axes[1].plot(1e3 * r, torsion_shear_stress(k * 0.5, r, J) / 1e6)
axes[1].set(xlabel="radius [mm]", ylabel="shear stress [MPa]", title="Shear across the wire at 0.5 rad")
fig.tight_layout()

plt.show()

torque ratios for doubling the diameter: [16. 16.]