N-link compound pendulum simulation

/ 6 Comments

For our Mechanical Dynamics final project, a friend and I decided to generalize the equations for compound pendulums to write a simulation program that could generate animations for n-link compound pendulums with user specified lengths, masses, angular displacement and velocities. We originally wanted to write a full blown physics engine, but quickly realized that two weeks was not enough time to do something of that scale, haha.

First things first, a demonstration! The following is an simulation of a 5 link compound pendulum:

Read our final report for a more formal/comprehensive explanation of the process!
Download the code!

These were the parameters entered in the GUI.

The 5 links are identical, each with a length of 1m, with a mass of 1kg. The initial thetas are relative to the previous pendulum.

The 5 links are identical, each with a length of 1m, with a mass of 1kg. The initial thetas are relative to the previous pendulum.

Our program also automatically generates energy and position/velocity/acceleration graphs.

The energy plots. Note how the yellow line (total energy) is mostly flat. It shows that barring calculation errors, our program is conserving energy correctly.

The energy plots. Note how the potential energy (PE) added to the kinetic energy (KE) results in a mostly flat total energy (TE). It shows that barring calculation errors, our program is conserving energy correctly.

The position/velocity/acceleration graphs give an idea of how each link in the pendulum is moving. However, it's a bit messy to read...

The position/velocity/acceleration graphs give an idea of how each link in the pendulum is moving. However, it’s a bit messy to read…

Now, how did we manage to do such a cool thing?

In order to simulate a compound pendulum, we first had to fully understand the equations behind a simple pendulum. In a simple pendulum, we ignore the mass moments of inertia and only worry about gravity and tension. Take a look at my post about the equations of motion of simple pendulums for a more detailed explanation.

To move from a simple pendulum to a compound pendulum, we first have to recognize that the core difference between the two is that we must worry about the mass moments of inertia in a compound pendulum. We are no longer treating the pendulum as a mass attached to a string, but as a rectangular bar.

compound_peng

Three link compound pendulum, each red dot represents the respective center of mass for each bar in the pendulum. Each link also has its own rotation frame relative to the inertial i-j-k frame.

We have to be able to transform between each frame of reference, so by doing some trigonometry, we find that the rotation matrices are:
\begin{bmatrix} \hat{i} \\ \hat{j} \end{bmatrix} = \begin{bmatrix} sin(\theta_1) & cos(\theta_1)\\ -cos(\theta_1) & sin(\theta_1) \end{bmatrix}\begin{bmatrix} \hat{b_1} \\ \hat{b_2} \end{bmatrix}
rot1

 

 

 

and

\begin{bmatrix} \hat{c_1} \\ \hat{c_2} \end{bmatrix} = \begin{bmatrix} cos(\theta_21) & sin(\theta_21)\\ -sin(\theta_21) & cos(\theta_21) \end{bmatrix}\begin{bmatrix} \hat{b_1} \\ \hat{b_2} \end{bmatrix}
Version 2

 

 

 

 

And ecetera for any future links to this pendulum.

We know we can describe the angular momentum of each pendulum as:
\omega = \begin{bmatrix} 0 \\ 0 \\ \dot{\theta} \\ \end{bmatrix}
Where theta is the respective theta from 0 for each pendulum. There is only rotation about the third axis, the one that’s jutting out of the page. (We can verify this via the right hand rule.)

The inertial tensor for a rectangular bar it’s center of mass is:
[I]_g = \begin{bmatrix} 0 & 0 & 0\\ 0 & \dfrac{mL^2}{12} & 0\\ 0 & 0 & \dfrac{mL^2}{12} \end{bmatrix}

Where m is the mass and L is the length of the bar.

Now, using the same method as in finding equations for a simple pendulum, we do the same by finding all the accelerations and forces to solve for the \ddot{\theta} for each link. For simple pendulums, the reaction force (tension) was always aligned with along the length of the string. However, in a compound pendulum, the reaction forces aren’t always aligned so we must solve for them as well. This means, for each link the pendulum, there are three unknowns that must be solved. These are \ddot{\theta}, R_x, R_y where R_x, R_y are the reaction forces in the two orthogonal directions of the attached reference frame. The idea is illustrated below:

Free body diagram of two example links in our compound pendulum.

Free body diagram of two example links in our compound pendulum.

Since there are three unknowns per link, we’ll need three equations multiplied by however many links we have. We can get the two of the three by following the same method as with simple pendulums, by collecting the accelerations and forces and solving F=ma. The last equation we find by using the equation:

\sum M = \dfrac{d}{dt}([I]*\omega)
(i.e, the sum of the moments is equal to the derivative of the angular momentum)

A moment is equal to the force crossed by the distance over which it’s exerted:
M = \vec{r}\times\vec{F}

As a quick example, this is how that calculation would work for the last link of our pendulum where the only moment would be as a result of the reaction force connecting it to the previous pendulum (B):
\sum M_g = (-\dfrac{L_2}{2} \hat{d_1} \times - B_1 \hat{d_1}) - (\dfrac{L_2}{2} \hat{d_1} \times -B_2 \hat{d_2})
\sum M_g = \dfrac{B_2L_2}{2} \hat{d_3}

\dfrac{d}{dt}([I]_g*\omega) = \dfrac{m_2L_2^2}{12} \ddot{\theta_3} \hat{d_3}

\dfrac{m_2L_2^2}{6} \ddot{\theta_3} - B_2 = 0

Working through all this math, we observe a pattern that every pendulum on the link is basically the same and because each has a reaction force on each end, and a theta relative to the vertical 0. The pendulums are connected by the reaction forces, as one can see in the free body diagram above, where R2 is on the bottom end of pendulum 1 and the top end of pendulum 2.

We generalized this to three equations:
\dfrac{m_n{L_n}^2}{12}\ddot{\theta_n}^2= \dfrac{-R_{nx}L_n}{2}sin(\theta_n)+\dfrac{R_{ny}L_n}{2}cos(\theta_n)-\dfrac{R_{(n+1)x}L_n}{2}sin(\theta_n)+\dfrac{R_{(n+1)y}L_n}{2}cos(\theta_n)

\sum_{i=0}^{n}(-m_n L_i cos(\theta_i)\ddot{\theta_i})-\dfrac{m_n L_n}{2}cos(\theta_n)\ddot{\theta_n}-R_{nx}+R_{(n+1)x}= \sum_{i=0}^{n}-m_n L_i \dot{\theta_i}^2 sin(\theta_i)+\dfrac{m_n L_n}{2}\dot{\theta_n}^2 sin(\theta_n)

\sum_{i=0}^{n}(m_n L_i sin(\theta_i)\ddot{\theta_i})- \dfrac{m_n L_n}{2}sin(\theta_n)\ddot{\theta_n}-R_{ny}+R_{(n+1)y}= \sum_{i=0}^{n}-m_n L_i \dot{\theta_i}^2 cos(\theta_i)+\dfrac{m_n L_n}{2}\dot{\theta_n}^2 cos(\theta_n)-m_ng

Which we then used as our basis for our N-case pendulum program.

All in all, this was a super fun project. I never thought that pendulums could be so exciting, but being able to tie this all back into a coding project was pretty awesome!
-Sophie

6 Comments

  2. Brad

    Would you happen to know how to solve for the equations of motion for the same n type problem using the Lagrangian Method instead?

    Reply

    • Sophie

      Hey Brad,

      I don’t know off the top of my head what the solution will be. However, I have a couple of gut feeling directions towards a solution?

      We mainly solved the problem by solving for a triple pendulum, then realizing that all other extensions to the problem to N-case would be variations of that pendulum. Page 2-3 of our final report sort of discuss this approach when we were searching for patterns to generalize. Let me know if that helps at all. I don’t think I’ve ever learned about legrangians, so it could be a fun weekend problem for me… :9

      (Also the grammar on that report is terrible, we were pretty sleep deprived while writing it.)

      Reply

      • Brad

        Hi Sophie. Thanks for responding, I really appreciate it. Our problem is a little different than yours. In ours, we are assuming mass less rods, with a point mass on the ends of each one between the different segments of the pendulum. Would you be able to help me find a generalized form of the Kinetic Energy of the system, and to the extent, the Lagrangian of the system?

        Reply

Leave a Reply

Your email address will not be published. Required fields are marked *