Spotlight #61 – Physics & Engineering — Archimedes, Buckling, PID Control & 3D Pendulum

I. Archimedes Principle — Buoyancy and Ship Stability

Archimedes' insight — that a submerged body is buoyed up by a force equal to the weight of fluid it displaces — is among the oldest quantitative results in physics and remains in daily use by naval architects. The upward buoyant force is:

F_b  =  ρ_fluid · g · V_submerged

For a floating body in static equilibrium, F_b exactly balances the weight W = m·g. The interesting physics begins when the vessel tilts. As a ship heels to angle θ, the underwater volume shifts to one side: the centre of buoyancy (the centroid of the displaced volume) migrates outward. The metacentre M is the point where the vertical line through the tilted centre of buoyancy intersects the vessel's centreline. If M lies above the centre of gravity G, the restoring moment rights the ship; if M lies below G, the overturning moment capsizes it.

Metacentric Height and the Restoring Moment

The vertical distance GM (metacentric height) completely characterises stability for small angles. The restoring moment is:

M_restore  =  W · GM · sin(θ)

A large positive GM (stiff ship) rights quickly but produces harsh rolling that is uncomfortable for passengers and dangerous for cargo. A small GM (tender ship) rolls sluggishly and may not recover from large waves. Passenger liners are designed to narrow tolerances; a fully-loaded container ship has a very different stability curve from the same vessel in ballast. The IMO requires ships to carry a Stability Information Booklet specifying how the metacentric height changes with loading condition — a direct application of Archimedes computed for every voyage.

The simulation computes the displaced volume numerically at each heel angle by integrating the hull geometry below the waterline. This captures the non-linear behaviour that linear metacentric theory misses: the GZ curve (righting lever vs. angle) can have a maximum, drop to zero (angle of vanishing stability), and even go negative at extreme angles where deck-edge immersion floods the hull.

II. Column Buckling — Euler’s Critical Load

A slender column under axial compression does not fail by direct crushing. Instead, at a critical load, it suddenly deflects sideways — buckles — losing its load-carrying capacity almost instantly. Leonhard Euler derived the critical load in 1744:

P_cr  =  π²EI / (KL)²

Here E is Young's modulus, I is the second moment of area of the cross-section, L is the column length, and K is the effective length factor determined by end conditions:

• K = 1.0 — both ends pinned (free to rotate, not to translate)
• K = 0.5 — both ends fixed (clamped against rotation)
• K = 0.7 — one end fixed, one pinned
• K = 2.0 — one end fixed, one free (flagpole)

The formula reveals immediately why tall, slender columns are dangerous: P_cr falls as L². Doubling the length quarters the critical load. The relevant geometric parameter is the slenderness ratio KL/r, where r = √(I/A) is the radius of gyration. For KL/r below roughly 120 (for steel), inelastic buckling governs: the material yields before the theoretical elastic buckling load is reached.

Johnson Parabola and the Transition to Yielding

The Johnson parabola approximates the transition between yielding and elastic buckling. For slenderness ratios below a critical value C_c = √(2π²E/F_y), the critical stress is:

F_cr  =  F_y · [ 1  -  (KL/r)² / (2·C_c²) ]

Above C_c, Euler's elastic formula applies. This piecewise model is the basis of steel column design in AISC specifications. The simulation colours the column by stress state — elastic (blue), inelastic (yellow), yielded (red) — so you can watch the failure mode transition as you adjust slenderness.

III. PID Controller — The Algorithm Running the World

Proportional-Integral-Derivative control is the most widely deployed algorithm in industrial automation. Estimates suggest that over 90% of control loops in process industries use PID. The control signal is:

u(t)  =  K_p · e(t)  +  K_i · ∫e(τ)dτ  +  K_d · de(t)/dt

where e(t) = r(t) - y(t) is the error between the setpoint r and measured output y. The proportional term reacts to current error. The integral term eliminates steady-state offset by accumulating past error. The derivative term anticipates future error by responding to the rate of change — it acts like a brake, damping overshoot.

Ziegler-Nichols Tuning and Anti-Windup

The classic Ziegler-Nichols method tunes PID gains experimentally: with only proportional control active, increase K_p until the system output oscillates continuously at sustained amplitude. Record this ultimate gain K_u and the ultimate period T_u. Then set:

K_p  =  0.6 · K_u
T_i  =  0.5 · T_u       (so K_i = K_p / T_i)
T_d  =  0.125 · T_u     (so K_d = K_p · T_d)

This produces a closed-loop response with roughly 25% overshoot and fast settling. For more conservative tuning, the Cohen-Coon or IMC methods are used. A practical issue always present is integral windup: when the plant saturates (a valve cannot open more than 100%), the integral term keeps accumulating error, building up a large value that then causes a prolonged overshoot when the plant comes out of saturation. Anti-windup strategies include conditional integration (stop accumulating when output is saturated) and back-calculation (subtract the saturation from the integral).

The simulation connects the PID to a second-order plant modelled as a mass-spring-damper system. You can directly observe the transition between overdamped (slow, no overshoot), underdamped (oscillating), and critically damped (fastest settling without overshoot) responses by dragging the gain sliders.

IV. 3D Pendulum — Foucault, Conical & Rosette Modes

A spherical pendulum has two degrees of freedom — polar angle θ from vertical and azimuthal angle φ around the vertical — coupled by the constraint that the bob stays on a sphere of radius L. The Lagrangian formulation, which tracks kinetic minus potential energy, yields two coupled second-order ODEs. For small amplitudes, the motion decomposes into two independent harmonic oscillators of the same frequency, meaning any initial condition produces a Lissajous-like figure in the horizontal plane: a straight line (planar pendulum), a circle (conical pendulum), or an ellipse (the general case).

Foucault Precession and Earth Rotation

When Earth's rotation is included as a Coriolis term, the apparent swing plane of a planar pendulum precesses slowly. The precession rate depends on latitude φ:

ω_F  =  Ω · sin(φ)

where Ω = 7.27 × 10⁻⁵ rad/s is Earth's rotation rate. At the North Pole (φ = 90°), the swing plane completes one full rotation per sidereal day (23h 56m). At London's latitude (51.5°N), the precession rate is Ω · sin(51.5°) ≈ 5.68 × 10⁻⁵ rad/s, giving a full rotation period of about 30.5 hours. Jean-Bernard-Léon Foucault hung a 67-metre pendulum from the dome of the Panthéon in Paris in 1851 and invited the public to watch the sand trace trace a rotating pattern on the floor — the first direct visual demonstration that Earth rotates.

The simulation time-compresses Earth's rotation by a factor of up to 1000×, making the precession visible in seconds. At mid-latitudes, the path traces the characteristic rosette pattern: successive swings are displaced by the Foucault angle, producing a star-like figure that gradually fills the plane.

V. Acoustic Lens — FDTD Wave Focusing

The Finite-Difference Time-Domain (FDTD) method solves the wave equation by discretising both space and time on a grid. For the 2D acoustic pressure field p(x,y,t) in a medium with local wave speed c(x,y):

∂²p/∂t²  =  c(x,y)² · ∇²p

The FDTD scheme uses a leapfrog discretisation: pressure and velocity are stored on staggered grids offset by half a time step. This second-order-accurate scheme is stable when the Courant-Friedrichs-Lewy (CFL) condition is satisfied: c·Δt/Δx ≤ 1/√2 in 2D. Violating CFL causes numerical instability — the simulation diverges exponentially.

Mur Absorbing Boundaries and Lens Design

Boundary conditions are critical in FDTD. Reflective (hard) walls simply enforce ∂p/∂n = 0 at the boundary. But reflections from the computational domain boundary corrupt simulations of radiation. The Mur absorbing boundary condition (1981) approximates the outgoing wave as a one-way equation and applies it at the grid edge, absorbing 95–99% of incident energy with minimal computational cost — a more modern alternative, the Perfectly Matched Layer (PML), achieves near-perfect absorption by introducing a dissipative gradient in the boundary cells.

A gradient-index acoustic lens works by varying c(x,y) across the aperture: slower wave speed at the centre, faster at the edges (or vice versa depending on lens type). Slower-propagating central rays lag behind outer rays, curving the wavefront into a convergent shape. The focal length follows the acoustic analogue of the lensmaker's equation. Fresnel zone plates achieve similar focusing by selective blocking rather than refraction: alternating annular zones of open and blocked aperture cause constructive interference at the focal point.