Reynolds and Sommerfeld: hydrodynamic lubrication (1886-1904)

Beauchamp Tower (1883) drilled an oil hole in the cap of a test bearing and found the oil pushing the plug out: the film carried the load on a pressure far above atmospheric. Osborne Reynolds (1886) explained it. A journal running off-centre drags oil into a converging wedge, and in a thin film the pressure obeys

\[\frac{d}{d\theta}\left(h^3\frac{dp}{d\theta}\right) = 6\mu\omega r^2\frac{dh}{d\theta}, \qquad h = c(1 + \varepsilon\cos\theta).\]

Arnold Sommerfeld (1904) integrated it in closed form for a long bearing. The load, viscosity, speed and geometry then enter only through the Sommerfeld number \(S = (r/c)^2\mu N/P\) (sommerfeld_number()), which fixes the eccentricity (long_bearing_eccentricity()). The example integrates Reynolds’s equation numerically and recovers Sommerfeld’s load and friction formulas.

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import cumulative_trapezoid, trapezoid

from engineeringkit.mechanical_engineering import long_bearing_eccentricity, petroff_friction_torque, sommerfeld_number

Integrate Reynolds’s equation round the journal

Integrating once gives \(dp/d\theta = 6\mu\omega r^2 (h - h^*)/h^3\). The constant \(h^*\) makes the pressure return to its starting value after one turn. The film force and the shear on the journal follow by integrating round the circumference, per metre of bearing length.

viscosity, omega, r, L, c = 0.04, 2 * np.pi * 40, 0.025, 1.0, 30e-6
theta = np.linspace(0, 2 * np.pi, 20001)


def film(eps):
    h = c * (1 + eps * np.cos(theta))
    h_star = trapezoid(h**-2, theta) / trapezoid(h**-3, theta)
    dp = 6 * viscosity * omega * r**2 * (h - h_star) / h**3
    p = cumulative_trapezoid(dp, theta, initial=0)
    p -= p.mean()
    load = np.hypot(trapezoid(p * np.cos(theta), theta), trapezoid(p * np.sin(theta), theta)) * r * L
    torque = trapezoid(viscosity * omega * r / h + h / (2 * r) * dp, theta) * r**2 * L
    return p, load, torque


for eps in (0.2, 0.5, 0.8):
    p, load, torque = film(eps)
    W_sommerfeld = 12 * np.pi * viscosity * omega * r**3 * L * eps / (c**2 * (2 + eps**2) * np.sqrt(1 - eps**2))
    T_sommerfeld = 4 * np.pi * viscosity * omega * r**3 * L * (1 + 2 * eps**2) / (c * (2 + eps**2) * np.sqrt(1 - eps**2))
    S = sommerfeld_number(viscosity, omega, load, r, L, c)
    print(f"eps = {eps}: load {load / 1e3:7.1f} kN per metre (Sommerfeld {W_sommerfeld / 1e3:7.1f})")
    print(f"          S = {S:.4f} -> eps = {long_bearing_eccentricity(S):.4f}")
    assert np.isclose(load, W_sommerfeld, rtol=1e-6) and np.isclose(torque, T_sommerfeld, rtol=1e-6)
    assert np.isclose(long_bearing_eccentricity(S), eps)

eps = 0.2: load   658.4 kN per metre (Sommerfeld   658.4)
          S = 0.0844 -> eps = 0.2000
eps = 0.5: load  1688.4 kN per metre (Sommerfeld  1688.4)
          S = 0.0329 -> eps = 0.5000
eps = 0.8: load  3323.1 kN per metre (Sommerfeld  3323.1)
          S = 0.0167 -> eps = 0.8000

Lightly loaded bearings tend to Petroff’s

As the Sommerfeld number grows the journal moves to the centre and the friction torque falls to Petroff’s value.

S = np.geomspace(0.01, 10, 60)
eps = long_bearing_eccentricity(S)
torque_ratio = (1 + 2 * eps**2) * 2 / ((2 + eps**2) * np.sqrt(1 - eps**2))
_, _, torque = film(0.01)
print(f"eps = 0.01: torque / Petroff = {torque / petroff_friction_torque(viscosity, r, L, c, omega):.6f}")
assert np.isclose(torque, petroff_friction_torque(viscosity, r, L, c, omega), rtol=1e-3)

fig, axes = plt.subplots(1, 2, figsize=(10, 3.8))
for e in (0.2, 0.5, 0.8):
    axes[0].plot(np.degrees(theta), film(e)[0] / 1e6, label=f"eps = {e}")
axes[0].set(xlabel="angle from widest gap [deg]", ylabel="film pressure [MPa]", title="Full Sommerfeld pressure")
axes[0].legend()
axes[1].semilogx(S, eps, label="eccentricity ratio")
axes[1].semilogx(S, torque_ratio, label="torque / Petroff torque")
axes[1].set(xlabel="Sommerfeld number S", title="Long journal bearing")
axes[1].legend()
fig.tight_layout()

plt.show()

eps = 0.01: torque / Petroff = 1.000200