Problem Statement

Someone posted a question on reddit, asking for the equations of motion for a complicated mechanical system:

At the time of writing, no one actually answered the question and the comments are just variations on "it's ugly and not solvable in closed form". In practice, it's not that ugly after all, and practically no mechanical systems are solvable in closed form, but that doesn't stop us from accurately simulating the system's behavior.

So, in this document, we'll walk through the classical mechanics of how to derive and simulate this system by briefly reviewing the fundamentals of classical mechanic and working up to actual system of interest, one feature at a time.

Classical Mechanics Review

indirectly $T,$ $V,$ of a system. The advantage here is that many systems (like the one we're studying here) would be practically impossible to analyze with Newtonian mechanics. As my classical mechanics professor used to say, someone trying to describe these kinds of systems in the Newtonian way would experience "an algebraic death"!

$\mathbf q$ $\mathcal L(\mathbf q, \dot{\mathbf q}) = T - V$ $\int \mathcal L\,dt$ stationary, which gives

\frac{d}{d t} \frac{\partial L}{\partial {\dot{q}}_{i}} - \frac{\partial L}{\partial q_{i}} = 0 .

$T = \tfrac12\,\dot{\mathbf q}^\top \mathsf M(\mathbf q)\,\dot{\mathbf q}$ $V = V(\mathbf q)$ . Then

\frac{\partial L}{\partial {\dot{q}}_{i}} = \sum_{j} M_{i j} {\dot{q}}_{j}, \frac{d}{d t} \frac{\partial L}{\partial {\dot{q}}_{i}} = \sum_{j} M_{i j} {\ddot{q}}_{j} + \sum_{j, k} \frac{\partial M_{i j}}{\partial q_{k}} {\dot{q}}_{j} {\dot{q}}_{k}, \frac{\partial L}{\partial q_{i}} = \frac{1}{2} \sum_{j, k} \frac{\partial M_{j k}}{\partial q_{i}} {\dot{q}}_{j} {\dot{q}}_{k} - \frac{\partial V}{\partial q_{i}},

so the equations of motion are

M (q) \ddot{q} + h (q, \dot{q}) + \nabla V (q) = 0, h_{i} = \sum_{j, k} (\frac{\partial M_{i j}}{\partial q_{k}} - \frac{1}{2} \frac{\partial M_{j k}}{\partial q_{i}}) {\dot{q}}_{j} {\dot{q}}_{k}

$\mathsf M$ $V$ .

The rocking-arc pendulum, one body at a time

1. The rocking arc

Start simple by just analyzing how the large arc rocks back and forth.

We assume it rolls on the ground without slipping, and as it does, its potential energy is given by the height of its center of mass:

V_{1} = M g d (1 - \cos θ)

$d = \frac{4\,(R_o^3 - R_i^3)}{3\pi\,(R_o^2 - R_i^2)}$ . The kinetic energy for this system depends on the arc's own moment of inertia about its center,

I_{O} = \frac{M}{2} (R_{o}^{2} + R_{i}^{2})

$C$ $I_C = I_G + M|CG|^2$ $I_G = I_O - Md^2$ $C$ $O$ $G$ $|CO| = R_o$ $|OG| = d$ $\theta$ between them),

| C G |^{2} = R_{o}^{2} + d^{2} - 2 R_{o} d \cos θ ⟹ I_{C} = I_{O} + M R_{o}^{2} - 2 M R_{o} d \cos θ,

T_{1} = \frac{1}{2} \underset{A_{0} (θ)}{\underset{⏟}{(I_{O} + M R_{o}^{2} - 2 M R_{o} d \cos (θ))}} {\dot{θ}}^{2}

which, after differentiation gives us the equation of motion

A_{0} (θ) \ddot{θ} + M R_{o} d \sin θ {\dot{θ}}^{2} + M g d \sin θ = 0

2. Adding the rolling disk

$\mathbf r_P = \mathbf r_O + \rho\,\mathbf u(\beta)$ $\rho = R_i - r$ $\dot\varphi = (R_i\dot\theta - \rho\dot\beta)/r$ .

The potential contribution of this disk is simple, we just need to add a contribution proportional to its center of mass vertical coordinate:

V_{2} = V_{1} + m_{r} g ρ (1 - \cos β)

The kinetic contribution has terms due to its own the translational and rotational energy, as well as a term from the coupling of its position to the rocking arc.

T_{2} = T_{1} + \frac{1}{2} m_{r} | {\dot{r}}_{P} |^{2} + \frac{1}{2} κ m_{r} (R_{i} \dot{θ} - ρ \dot{β})^{2}, | {\dot{r}}_{P} |^{2} = R_{o}^{2} {\dot{θ}}^{2} - 2 R_{o} ρ \cos β \dot{θ} \dot{β} + ρ^{2} {\dot{β}}^{2}

Those two contributions update our equations of motion, giving us a new effective mass matrix

M_{2} = (\begin{matrix} A & B \\ B & D \end{matrix}), \begin{aligned} A & = A_{0} (θ) + m_{r} R_{o}^{2} + κ m_{r} R_{i}^{2} \\ B & = - m_{r} R_{o} ρ \cos β - κ m_{r} R_{i} ρ \\ D & = (1 + κ) m_{r} ρ^{2} \end{aligned}

\begin{aligned} A \ddot{θ} + B \ddot{β} + M R_{o} d \sin θ {\dot{θ}}^{2} + m_{r} R_{o} ρ \sin β {\dot{β}}^{2} + M g d \sin θ & = 0 \\ B \ddot{θ} + D \ddot{β} + m_{r} g ρ \sin β & = 0 \end{aligned}

3. Adding the hinged pendulum

$P$ $\mathbf r_\pm = \mathbf r_P + l_p\,\mathbf u(\alpha) \pm \tfrac{w_p}{2}\,\mathbf t(\alpha)$ $m = 2m_p$ $J = m\left(l_p^2 + \tfrac14 w_p^2\right)$ :

T_{3} = T_{2} + \frac{1}{2} m | {\dot{r}}_{P} |^{2} + m l_{p} \dot{α} (- R_{o} \cos α \dot{θ} + ρ \cos (β - α) \dot{β}) + \frac{1}{2} J {\dot{α}}^{2}

V_{3} = V_{2} + m g ρ (1 - \cos β) + m g l_{p} (1 - \cos α)

$\mu = m_r + m$ ,

M_{3} = (\begin{matrix} A & B & E \\ B & D & F \\ E & F & J \end{matrix}), \begin{aligned} A & = A_{0} (θ) + μ R_{o}^{2} + κ m_{r} R_{i}^{2} \\ B & = - μ R_{o} ρ \cos β - κ m_{r} R_{i} ρ \\ D & = (μ + κ m_{r}) ρ^{2} \\ E & = - m l_{p} R_{o} \cos α \\ F & = m l_{p} ρ \cos (β - α) \end{aligned}

h_{3} = (\begin{matrix} M R_{o} d \sin θ {\dot{θ}}^{2} + μ R_{o} ρ \sin β {\dot{β}}^{2} + m l_{p} R_{o} \sin α {\dot{α}}^{2} \\ m l_{p} ρ \sin (β - α) {\dot{α}}^{2} \\ - m l_{p} ρ \sin (β - α) {\dot{β}}^{2} \end{matrix}), \nabla V_{3} = (\begin{matrix} M g d \sin θ \\ μ g ρ \sin β \\ m g l_{p} \sin α \end{matrix})

M_{3} \ddot{q} + h_{3} + \nabla V_{3} = 0

4. Adding the springs

$\mathbf a_\pm = (\pm d_s, 0)$ $\mathbf B_\pm = \mathbf r_O + \mathsf Q(\theta)\,(\pm R_o\cos\beta_s,\, -R_o\sin\beta_s)$ $\mathsf Q$ $P$ $\mathbf r_Q = \mathbf r_O + R_a\mathbf u(\theta)$ . Their lengths are

L_{\pm} = | B_{\pm} - a_{\pm} |, s = \sqrt{R_{a}^{2} + ρ^{2} - 2 R_{a} ρ \cos (β - θ)},

$L_0$ $s_0$ $\mathbf q = \mathbf 0$ ),

T_{4} = T_{3}, V_{4} = V_{3} + \frac{k}{2} [(L_{+} - L_{0})^{2} + (L_{-} - L_{0})^{2} + (s - s_{0})^{2}]

\nabla V_{4} = \nabla V_{3} + (\begin{matrix} τ - f \\ f \\ 0 \end{matrix}), τ = k \sum_{\pm} (L_{\pm} - L_{0}) \frac{d L_{\pm}}{d θ}, f = k (s - s_{0}) \frac{R_{a} ρ \sin (β - θ)}{s}

M_{3} \ddot{q} + h_{3} + \nabla V_{4} = 0

Just for fun, we also plot the sum of the different kinetic and potential energy terms for the ball, arc, pendulum and springs to verify that the total energy is conserved.