Rigid Bodies#

Contents

Phantom Torque

Phantom Torque #

Here I would like to try to clarify a point mentioned in class about what Turner and Turner [Turner and Turner, 2010] call “phantom torque”. In the textbook, [Fetter and Walecka, 2003] show in (27.1) that the rate of change of angular momentum \(\vect{L}\) in an inertial frame is equal to the external torque \(\vect{\Gamma}^{(e)}\):

\[\begin{gather*} \left. \diff{\vect{L}}{t}\right|_{\text{inertial}} = \vect{\Gamma}^{(e)} \tag{27.1} \end{gather*}\]

in two cases:

About an origin fixed in some inertial frame.
About the center of mass, even if it accelerates. Note: the “inertial” in (27.1) refers to the frame not rotating about the center of mass.

The paper by Turner and Turner [Turner and Turner, 2010] points out a potential fallacy if the object is asymmetric such that .

Consider an object rotating with angle \(\theta\) and momentum of inertia \(I_p(\theta)\) about a point of contact where \(I_p'(\theta) \neq 0\). This dependence of the momentum of inertia on the coordinate is responsible for the “phantom torque”, but the actual manifestation is subtle as we now show. Note that for a rolling coin or ball which is symmetric, \(I_{p}(\theta) = I_{p}\) is constant, which is why this issue is often not discussed with many classical textbook examples.

The example raised in [Turner and Turner, 2010] is that of a half-hoop rolling without slipping along a flat surface with the motion described by the angle \(\theta\) such that \(\dot{\theta}\) is the angular velocity.

Show code cell source Hide code cell source

%%manim -v WARNING -qm RollingHoop
from manim import *
from types import SimpleNamespace


class Hoop:
    R = 4.0
    theta0 = 0.6
    r_cm = 2 / np.pi * R
    m = 1.0
    I_o = m * R ** 2
    I_cm = I_o - m * r_cm ** 2
    shift = UP

    def __init__(self, **kw):
        for k in kw:
            assert hasattr(self, k)
            setattr(self, k, kw.pop(k))
        assert not kw
        
        self.objects = self.annotations

    def get_arc(self, theta=None):
        """Return a (arc, points) for a hoop"""
        color = 'white'
        if theta is None:
            theta = self.theta0
            color = 'cyan'
            
        z_o = -self.R * theta + 0j
        z_cm = z_o - 1j * self.r_cm * np.exp(theta * 1j)
        z_p = z_o - 1j * self.R
        o, p, cm = [np.array([_z.real, _z.imag, 0]) for _z in [z_o, z_p, z_cm]]
        e1 = o + self.R * np.array([np.cos(theta), np.sin(theta), 0])

        arc = Arc(
            radius=self.R, start_angle=PI + theta, angle=PI, arc_center=o,
            color=color
        )
        arc.stroke_width = 10.0
        points = SimpleNamespace(O=o+self.shift, P=p+self.shift, C=cm+self.shift)
        arc.shift(self.shift)
        
        # These are the actual objects we display
        return (arc, points)

    @property
    def annotations(self):
        """Return a `Group` object with the hoop and all of the annotations."""
        arc, points = self.get_arc()
        
        objs = []
        P, C, O = points.P, points.C, points.O
        for _p0, _p1, _t in [
            (arc.get_end(), O, "R"),
            (C, O, r"a=\tfrac{2}{\pi}R"),
        ]:
            objs.append(BraceBetweenPoints(_p0, _p1))
            objs.append(objs[-1].get_tex(_t))
            #objs[-1].font_size = 40
        objs.extend(
            [
                LabeledDot("O", point=O),
                LabeledDot("P", point=P),
                LabeledDot("C", point=C),
            ]
        )
        
        OC = DashedLine(O, C, color='cyan')
        OP = DashedLine(O, P, color='cyan')
        theta = Angle(OP, OC, radius=0.7, dot_radius=2.0)
        theta_l = Tex(r"$\theta$")
        _m = theta.get_midpoint()
        theta_l.next_to(O + 1.5*(_m-O), 0)
        objs.extend([OC, OP, theta, theta_l])

        floor = Line((-6, -self.R, 0), (6, -self.R, 0))
        floor.shift(self.shift)
        floor.stroke_width = arc.stroke_width
        self.arc = arc
        annotations = Group(arc, floor, *objs)
        return annotations


class HoopAnimation(Animation):
    def __init__(self, hoop, **kw):
        super().__init__(hoop.get_arc(theta=hoop.theta0)[0], **kw)
        self.hoop = hoop
        
    def interpolate_mobject(self, alpha):
        # Fake dynamics
        theta = self.hoop.theta0 * np.cos(2*np.pi * alpha)
        self.mobject.points = self.hoop.get_arc(theta=theta)[0].points
        

class RollingHoop(Scene):
    def construct(self):
        hoop = Hoop()
        self.add(hoop.annotations)
        #self.play(HoopAnimation(hoop=hoop), run_time=3, rate_func=linear)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[2], line 1
----> 1 get_ipython().run_cell_magic('manim', '-v WARNING -qm RollingHoop', 'from manim import *\nfrom types import SimpleNamespace\n\n\nclass Hoop:\n    R = 4.0\n    theta0 = 0.6\n    r_cm = 2 / np.pi * R\n    m = 1.0\n    I_o = m * R ** 2\n    I_cm = I_o - m * r_cm ** 2\n    shift = UP\n\n    def __init__(self, **kw):\n        for k in kw:\n            assert hasattr(self, k)\n            setattr(self, k, kw.pop(k))\n        assert not kw\n        \n        self.objects = self.annotations\n\n    def get_arc(self, theta=None):\n        """Return a (arc, points) for a hoop"""\n        color = \'white\'\n        if theta is None:\n            theta = self.theta0\n            color = \'cyan\'\n            \n        z_o = -self.R * theta + 0j\n        z_cm = z_o - 1j * self.r_cm * np.exp(theta * 1j)\n        z_p = z_o - 1j * self.R\n        o, p, cm = [np.array([_z.real, _z.imag, 0]) for _z in [z_o, z_p, z_cm]]\n        e1 = o + self.R * np.array([np.cos(theta), np.sin(theta), 0])\n\n        arc = Arc(\n            radius=self.R, start_angle=PI + theta, angle=PI, arc_center=o,\n            color=color\n        )\n        arc.stroke_width = 10.0\n        points = SimpleNamespace(O=o+self.shift, P=p+self.shift, C=cm+self.shift)\n        arc.shift(self.shift)\n        \n        # These are the actual objects we display\n        return (arc, points)\n\n    @property\n    def annotations(self):\n        """Return a `Group` object with the hoop and all of the annotations."""\n        arc, points = self.get_arc()\n        \n        objs = []\n        P, C, O = points.P, points.C, points.O\n        for _p0, _p1, _t in [\n            (arc.get_end(), O, "R"),\n            (C, O, r"a=\\tfrac{2}{\\pi}R"),\n        ]:\n            objs.append(BraceBetweenPoints(_p0, _p1))\n            objs.append(objs[-1].get_tex(_t))\n            #objs[-1].font_size = 40\n        objs.extend(\n            [\n                LabeledDot("O", point=O),\n                LabeledDot("P", point=P),\n                LabeledDot("C", point=C),\n            ]\n        )\n        \n        OC = DashedLine(O, C, color=\'cyan\')\n        OP = DashedLine(O, P, color=\'cyan\')\n        theta = Angle(OP, OC, radius=0.7, dot_radius=2.0)\n        theta_l = Tex(r"$\\theta$")\n        _m = theta.get_midpoint()\n        theta_l.next_to(O + 1.5*(_m-O), 0)\n        objs.extend([OC, OP, theta, theta_l])\n\n        floor = Line((-6, -self.R, 0), (6, -self.R, 0))\n        floor.shift(self.shift)\n        floor.stroke_width = arc.stroke_width\n        self.arc = arc\n        annotations = Group(arc, floor, *objs)\n        return annotations\n\n\nclass HoopAnimation(Animation):\n    def __init__(self, hoop, **kw):\n        super().__init__(hoop.get_arc(theta=hoop.theta0)[0], **kw)\n        self.hoop = hoop\n        \n    def interpolate_mobject(self, alpha):\n        # Fake dynamics\n        theta = self.hoop.theta0 * np.cos(2*np.pi * alpha)\n        self.mobject.points = self.hoop.get_arc(theta=theta)[0].points\n        \n\nclass RollingHoop(Scene):\n    def construct(self):\n        hoop = Hoop()\n        self.add(hoop.annotations)\n        #self.play(HoopAnimation(hoop=hoop), run_time=3, rate_func=linear)\n')

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/IPython/core/interactiveshell.py:2493, in InteractiveShell.run_cell_magic(self, magic_name, line, cell)
   2491 with self.builtin_trap:
   2492     args = (magic_arg_s, cell)
-> 2493     result = fn(*args, **kwargs)
   2495 # The code below prevents the output from being displayed
   2496 # when using magics with decorator @output_can_be_silenced
   2497 # when the last Python token in the expression is a ';'.
   2498 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/manim/utils/ipython_magic.py:141, in ManimMagic.manim(self, line, cell, local_ns)
    139     SceneClass = local_ns[config["scene_names"][0]]
    140     scene = SceneClass(renderer=renderer)
--> 141     scene.render()
    142 finally:
    143     # Shader cache becomes invalid as the context is destroyed
    144     shader_program_cache.clear()

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/manim/scene/scene.py:223, in Scene.render(self, preview)
    221 self.setup()
    222 try:
--> 223     self.construct()
    224 except EndSceneEarlyException:
    225     pass

File <string>:97, in construct(self)

File <string>:20, in __init__(self, **kw)

File <string>:49, in annotations(self)

File <string>:35, in get_arc(self, theta)

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/manim/mobject/geometry/arc.py:314, in Arc.__init__(self, radius, start_angle, angle, num_components, arc_center, **kwargs)
    312 self.angle: float = angle
    313 self._failed_to_get_center: bool = False
--> 314 super().__init__(**kwargs)

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/manim/mobject/geometry/arc.py:101, in TipableVMobject.__init__(self, tip_length, normal_vector, tip_style, **kwargs)
     99 self.normal_vector: Vector = normal_vector
    100 self.tip_style: dict = tip_style
--> 101 super().__init__(**kwargs)

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/manim/mobject/types/vectorized_mobject.py:153, in VMobject.__init__(self, fill_color, fill_opacity, stroke_color, stroke_opacity, stroke_width, background_stroke_color, background_stroke_opacity, background_stroke_width, sheen_factor, joint_type, sheen_direction, close_new_points, pre_function_handle_to_anchor_scale_factor, make_smooth_after_applying_functions, background_image, shade_in_3d, tolerance_for_point_equality, n_points_per_cubic_curve, **kwargs)
    151 self.tolerance_for_point_equality: float = tolerance_for_point_equality
    152 self.n_points_per_cubic_curve: int = n_points_per_cubic_curve
--> 153 super().__init__(**kwargs)
    154 self.submobjects: list[VMobject]
    156 # TODO: Find where color overwrites are happening and remove the color doubling
    157 # if "color" in kwargs:
    158 #     fill_color = kwargs["color"]
    159 #     stroke_color = kwargs["color"]

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/manim/mobject/mobject.py:115, in Mobject.__init__(self, color, name, dim, target, z_index)
    113 self.updaters: list[Updater] = []
    114 self.updating_suspended = False
--> 115 self.color = ManimColor.parse(color)
    117 self.reset_points()
    118 self.generate_points()

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/manim/utils/color/core.py:653, in ManimColor.parse(cls, color, alpha)
    651 if isinstance(color, (list, tuple)):
    652     return [cls(c, alpha) for c in color]  # type: ignore
--> 653 return cls(color, alpha)

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/manim/utils/color/core.py:121, in ManimColor.__init__(self, value, alpha)
    114         self._internal_value = ManimColor._internal_from_hex_string(
    115             result.group(), alpha
    116         )
    117     else:
    118         # This is not expected to be called on module initialization time
    119         # It can be horribly slow to convert a string to a color because
    120         # it has to access the dictionary of colors and find the right color
--> 121         self._internal_value = ManimColor._internal_from_string(value)
    122 elif isinstance(value, (list, tuple, np.ndarray)):
    123     length = len(value)

File ~/checkouts/readthedocs.org/user_builds/physics-521-classical-mechanics-i/conda/latest/lib/python3.12/site-packages/manim/utils/color/core.py:357, in ManimColor._internal_from_string(name)
    355     return _all_color_dict[upper_name]._internal_value
    356 else:
--> 357     raise ValueError(f"Color {name} not found")

ValueError: Color cyan not found

Here we have half of a hoop of radius \(R\) and mass \(m\), with center of mass at \(C\) which is distance \(a=2R/\pi\) from the center of the hoop \(O\). The external torque is provided by the downward acceleration due to gravity \(g>0\). The hoop rolls without slipping along the floor with point of contact \(P\).

If you have not solved this problem before, please stop for a moment and try it with the coordinate \(\theta\) being the angle of rotation with \(\theta = 0\) the equilibrium position with \(O\), \(C\), and \(P\) vertically aligned.

Solution

For the hoop, the center of mass is located at \(a = 2R/\pi\), and the moment of inertia about the contact point \(P\) can be found by using the parallel axis theorem twice, starting from \(I_{O} = mR^2 = I_{C} + ma^2\) and \(I_{P}(\theta) = I_{C} + m(R^2 + a^2 - 2aR\cos\theta) = 2mR(R-a\cos\theta)\). This gives the following Lagrangian and solution:

\[\begin{gather*} L(\theta, \dot{\theta}) = mR(R-a\cos\theta)\dot{\theta}^2 + mga\cos\theta,\\ \diff{}{t}\pdiff{L}{\dot{\theta}} = \pdiff{L}{\theta},\\ 2mR(R-a\cos\theta)\ddot{\theta} = -mRa\sin\theta\dot{\theta}^2 - mga\sin\theta,\\ \left(1 - \frac{2}{\pi}\cos\theta\right)\ddot{\theta} = -\frac{\sin\theta}{\pi} \dot{\theta}^2 - \frac{\sin\theta}{\pi}\frac{g}{R}. \end{gather*}\]

This is (13) from [Turner and Turner, 2010].

A careful way of proceeding is to use a Lagrangian. Here we work about the instantaneously stationary point \(P\), so the kinetic energy is purely rotational:

\[\begin{split} L(\theta, \dot{\theta}) = \frac{1}{2}I_{p}(\theta)\dot{\theta}^2 - V(\theta),\\ p_\theta = L = \pdiff{L}{\dot{\theta}} = I_{p}(\theta)\dot{\theta}, \qquad \dot{L} = \pdiff{L}{\theta}. \end{split}\]

Consider carefully the terms in the last equation of motion \(\dot{L} = \partial L/\partial \theta\):

\[\begin{split} \overbrace{I_{p}'(\theta)\dot{\theta}^2}^{a} + I_{p}(\theta)\ddot{\theta} = \overbrace{\frac{1}{2}I'_{p}(\theta)\dot{\theta}^2}^{b} - V'(\theta),\\ I_{p}(\theta)\ddot{\theta} = \underbrace{- V'(\theta)}_{\Gamma^{(e)}} + \underbrace{\frac{-I'_{p}(\theta)}{2}\dot{\theta}^2}_{\text{phantom torque}}. \end{split}\]

The two terms \(a\) and \(b\) indicated give rise to the “phantom torque” \(b-a\), but notice that the point is subtle. [Turner and Turner, 2010] point out that both terms are missing if one uses the \(F = ma\) form \(\Gamma = I \dot{\omega}\), but even using (27.1) one may be misled into finding the piece \(a\) but not \(b\):

\[\begin{gather*} \left.\diff{\vect{L}}{t}\right|_{\text{inertial}} = \mat{I}\cdot \ddot{\vect{\theta}} + \underbrace{\dot{\mat{I}}\cdot \dot{\vect{\theta}}}_{a} = \vect{\Gamma}^{(e)}. \end{gather*}\]

Incomplete Discussion

This is an incomplete discussion about what goes wrong and how to fix it. One account is given in [Jensen, 2011] but, while correct, I do not find any of the four forms he discusses very intuitive.

The simple attempts fail because the reference point \(P\) may be accelerating. Instead with need the distance from \(C\) to the point \(P_0=P\) in a locally inertial frame where \(P_0\) does not depend on \(\theta\). (In this problem, \(P_0\) is fixed since the hoop does not slip against the stationary floor). In this problem, \(P_0\) lies along the floor at distance \(-R\theta_0\) from the equilibrium position where \(\theta_0 = \theta\) at the time we are discussing, but does not change. The error in blind calculation of \(\dot{I} = I_{P}'(\theta)\dot{\theta}\) is that the expression for \(I_{P}(\theta)\) does not keep this \(\theta_0\) fixed. Doing this properly, we find that the the center of mass is at

\[ x_{C} = a\sin\theta - \theta R, \qquad y_{C} = -a\cos\theta. \]

Thus, the distance to the fixed point \(P_0\) from the center of mass \(C\) is:

\[ \sqrt{\Bigl(x_{C}(\theta) + R\theta_0\Bigr)^2 + y_{C}^2(\theta)}. \]

The moment of inertia about the point \(P_0\) is thus:

\[\begin{split} I_{P_0}(\theta, \theta_0) = I_{C} + m\Bigl(x_{C}(\theta) + \theta_0 R\Bigr)^2 + my_{C}^2(\theta)\\ = I_{P}(\theta) - 2mRx_{C}(\theta)(\theta - \theta_0) - mR^2(\theta^2 - \theta_0^2). \end{split}\]

Note that \(I_{P}(\theta_0) = I_{P_0}(\theta_0)\), but the correct variation is

\[\begin{gather*} \dot{I} = \left.\diff{I_{P_0}}{t}\right|_{\text{inertial}} = \pdiff{I_{P_0}(\theta, \theta_0)}{\theta}\dot{\theta},\\ \left.\pdiff{I_{P_0}(\theta, \theta_0)}{\theta} \right|_{\theta=\theta_0}= 2m\bigl(x_{C}(\theta_0) + \theta_0R\bigr)x_{C}'(\theta_0) + 2m y_{C}(\theta_0)y_{C}'(\theta_0)\\ = I_{P}'(\theta_0) - 2mR\bigl(x_{C}(\theta_0) + \theta_0R\bigr) \end{gather*}\]

\[\begin{split} I_{P}'(\theta) = 2m\bigl(x_{C}(\theta) + \theta R\bigr)(x_{C}'(\theta) + R) + 2m y_{C}(\theta)y_{C}'(\theta)\\ \end{split}\]

Rigid Bodies

Contents

Rigid Bodies#

Phantom Torque#

Phantom Torque #