.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/mechanical_engineering/dynamics/plot_01_huygens_moment_of_inertia.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via JupyterLite. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_mechanical_engineering_dynamics_plot_01_huygens_moment_of_inertia.py: Huygens's compound pendulum and the parallel-axis theorem (1673) ================================================================ Christiaan Huygens (1673) asked how long a simple pendulum must be to swing in step with a rigid body pivoted at a point. His answer, the length :math:`L = I_O/(m d)` to the *centre of oscillation*, needs the body's moment of inertia about the pivot, which Jakob Steiner later wrote in general as :math:`I_O = I_G + m d^2` (:func:`~engineeringkit.mechanical_engineering.parallel_axis`). The example checks the theorem by summing :math:`r^2\,dm` over a disk, recovers Huygens's rule that pivot and centre of oscillation can be swapped without changing the period, and finds the correction he applied to a clock's spherical bob. .. GENERATED FROM PYTHON SOURCE LINES 18-26 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from engineeringkit.integrators import integrate_ode from engineeringkit.mechanical_engineering import angular_momentum, disk_inertia, parallel_axis, solid_sphere_inertia g = 9.81 .. GENERATED FROM PYTHON SOURCE LINES 27-32 The parallel-axis theorem by direct summation --------------------------------------------- Cut a uniform disk into small squares and add up :math:`r^2\,dm` about a pivot a distance :math:`d` from the centre. .. GENERATED FROM PYTHON SOURCE LINES 32-44 .. code-block:: Python m, R = 2.0, 0.15 n = 1201 x, y = np.meshgrid(np.linspace(-R, R, n), np.linspace(-R, R, n)) inside = x**2 + y**2 <= R**2 dm = m / inside.sum() for d in (0.0, 0.05, 0.10, 0.15): summed = np.sum(((x[inside] - d) ** 2 + y[inside] ** 2) * dm) steiner = parallel_axis(disk_inertia(m, R), m, d) print(f"d = {d:.2f} m: summed I = {summed:.6f}, I_G + m d^2 = {steiner:.6f} kg m^2") assert np.isclose(summed, steiner, rtol=2e-3) .. rst-class:: sphx-glr-script-out .. code-block:: none d = 0.00 m: summed I = 0.022499, I_G + m d^2 = 0.022500 kg m^2 d = 0.05 m: summed I = 0.027499, I_G + m d^2 = 0.027500 kg m^2 d = 0.10 m: summed I = 0.042499, I_G + m d^2 = 0.042500 kg m^2 d = 0.15 m: summed I = 0.067499, I_G + m d^2 = 0.067500 kg m^2 .. GENERATED FROM PYTHON SOURCE LINES 45-52 Period of the disk pendulum and Huygens's reciprocity ----------------------------------------------------- The period :math:`2\pi\sqrt{I_O/(m g d)}` is smallest when the pivot sits at the radius of gyration :math:`k = R/\sqrt 2`. A pivot at the centre of oscillation, :math:`L - d` on the other side of the centre, gives the same period: Kater's reversible pendulum (1818) measured :math:`g` this way. .. GENERATED FROM PYTHON SOURCE LINES 52-65 .. code-block:: Python d = np.linspace(0.01, 0.3, 300) L = parallel_axis(disk_inertia(m, R), m, d) / (m * d) period = 2 * np.pi * np.sqrt(L / g) k = np.sqrt(disk_inertia(m, R) / m) d1 = 0.04 L1 = parallel_axis(disk_inertia(m, R), m, d1) / (m * d1) d2 = L1 - d1 L2 = parallel_axis(disk_inertia(m, R), m, d2) / (m * d2) print(f"pivot {d1:.3f} m: L = {L1:.4f} m; pivot {d2:.4f} m: L = {L2:.4f} m") assert np.isclose(L1, L2) assert np.isclose(d[np.argmin(period)], k, atol=d[1] - d[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none pivot 0.040 m: L = 0.3212 m; pivot 0.2812 m: L = 0.3212 m .. GENERATED FROM PYTHON SOURCE LINES 66-72 Simulated swing --------------- Integrate :math:`I_O\ddot\phi = -m g d \sin\phi` at a small amplitude and time the zero crossings. The angular momentum :math:`I_O\dot\phi` peaks as the disk passes the bottom. .. GENERATED FROM PYTHON SOURCE LINES 72-83 .. code-block:: Python I_O = parallel_axis(disk_inertia(m, R), m, d1) t = np.linspace(0, 10, 20001) swing = integrate_ode(lambda tt, s: [s[1], -m * g * d1 * np.sin(s[0]) / I_O], [0.02, 0.0], t) phi = swing.history[:, 0] crossings = t[1:][np.diff(np.sign(phi)) != 0] measured = 2 * np.mean(np.diff(crossings)) print(f"simulated period {measured:.4f} s, Huygens 2 pi sqrt(L/g) = {2 * np.pi * np.sqrt(L1 / g):.4f} s") print(f"peak angular momentum {np.max(angular_momentum(I_O, swing.history[:, 1])):.2e} kg m^2/s") assert np.isclose(measured, 2 * np.pi * np.sqrt(L1 / g), rtol=1e-3) .. rst-class:: sphx-glr-script-out .. code-block:: none simulated period 1.1371 s, Huygens 2 pi sqrt(L/g) = 1.1370 s peak angular momentum 2.84e-03 kg m^2/s .. GENERATED FROM PYTHON SOURCE LINES 84-89 The clock bob correction ------------------------ A sphere of radius :math:`a` on a light rod of length :math:`\ell` beats as a simple pendulum of length :math:`\ell + \tfrac25 a^2/\ell`. .. GENERATED FROM PYTHON SOURCE LINES 89-106 .. code-block:: Python rod, a, bob = 0.994, 0.05, 1.0 L_bob = parallel_axis(solid_sphere_inertia(bob, a), bob, rod) / (bob * rod) print(f"equivalent length {L_bob:.6f} m, correction {1e3 * (L_bob - rod):.3f} mm") assert np.isclose(L_bob - rod, 0.4 * a**2 / rod) fig, axes = plt.subplots(1, 2, figsize=(10, 3.8)) axes[0].plot(d, period) axes[0].axvline(k, color="0.5", ls="--", label="radius of gyration") axes[0].plot([d1, d2], [2 * np.pi * np.sqrt(L1 / g)] * 2, "o", label="reciprocal pivots") axes[0].set(xlabel="pivot distance from centre [m]", ylabel="period [s]", title="Disk pendulum") axes[0].legend() axes[1].plot(t, np.degrees(phi)) axes[1].set(xlabel="time [s]", ylabel="angle [deg]", title=f"Swing about the pivot at {d1} m") fig.tight_layout() plt.show() .. image-sg:: /api/gallery/mechanical_engineering/dynamics/images/sphx_glr_plot_01_huygens_moment_of_inertia_001.png :alt: Disk pendulum, Swing about the pivot at 0.04 m :srcset: /api/gallery/mechanical_engineering/dynamics/images/sphx_glr_plot_01_huygens_moment_of_inertia_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none equivalent length 0.995006 m, correction 1.006 mm .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.079 seconds) .. _sphx_glr_download_api_gallery_mechanical_engineering_dynamics_plot_01_huygens_moment_of_inertia.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../../../../lite/lab/index.html?path=api/gallery/mechanical_engineering/dynamics/plot_01_huygens_moment_of_inertia.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_01_huygens_moment_of_inertia.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_01_huygens_moment_of_inertia.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_01_huygens_moment_of_inertia.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_