The Spring Pendulum

The pendulum was studied scientifically by Galileo in 1602 and contemplated by ordinary folks probably since the beginning of conscious human thought.

For a simple gravity pendulum, which is a mathematically idealized pendulum, the period is (assuming angles smaller than \(\pm6^{\circ}\)):

\[T_0=2\pi\sqrt{\frac{L}{g}}\quad for\ \theta_{0}\ll0.1radians\]

Where \(L\) is the length of the pendulum, \(g\) is the acceleration due to gravity, and \(\theta_0\) is the initial angle of the mass. If you have ever taken a dynamics or a controls course you have surely derived the humble pendulum’s equations of motion. This was probably followed by the derivation of a cart with an attached pendulum. I think there is another type of pendulum that should be derived before the cart, the spring pendulum. Although it is a tad more complex, I think it is more interesting and the position vs. time plot is pretty neat.

The Spring Pendulum

The spring pendulum is a 2 degree of freedom system. This means you need two independent variables to express the location of the mass. For this demonstration I will be choosing \(L\) and \(\theta\). We also need two frames to help with describing the motion, I will use \(A\) as the stationary frame and \(B\) as the rotating frame.

Sympy Code

We will begin by importing the necessary packages. I really like the Pydy package and I urge anyone to try it out.

import sympy
from sympy import symbols, init_printing
import sympy.physics.mechanics as me
init_printing()
import matplotlib.pyplot as plt
import numpy as np
from pydy.system import System
from numpy import linspace
%matplotlib inline

We will now create the dynamic variables, the velocities, and the constant symbols. \(L_0\) is the equilibrium length of the spring.

# Create the variables
L, theta = me.dynamicsymbols('L theta')

# Create the velocities
L_dot, theta_dot = me.dynamicsymbols('L_dot theta_dot')

# Create the constants
m, k, g, t, L_0 = sympy.symbols('m k g t L_0')

Create the A frame, the B frame, and set the angular velocity of B. Because \(+Y\) is pointing downward \(+Z\) is pointing away from the reader (right-hand rule). This makes a clockwise rotation positive. In the drawing \(\theta\) is moving counterclockwise therefore, \(\dot{\theta}\) is negative.

# Create the world frame
A = me.ReferenceFrame('A')

# Create the pendulum frame
B = A.orientnew('B', 'axis', [theta, -A.z])

# Set the rotation of the pendulum frame
B.set_ang_vel(A, -theta_dot * A.z)

Now the two points are created: \(O\) the origin and \(m\) the mass point. Once the origin point has been created we can define the mass point as being a distance \(L\) from the origin. This is where frames come to the rescue. It is much simpler to define the point in its own frame and have Sympy do all the work. Although this is a simplistic example, this is a huge helper when systems grow. I will display the location of point \(m\) in the \(A\) frame to show the complexity of manually calculating the location of \(m\).

# Create the Origin
O = me.Point('O')

# Create the mass point
P = O.locatenew('P', L * B.y)

# Display the mass point location in the A frame
P.pos_from(O).express(A)

\[L \operatorname{sin}\left(\theta\right)\mathbf{\hat{a}_x} + L \operatorname{cos}\left(\theta\right)\mathbf{\hat{a}_y}\]

Now the origin’s velocity is set and the velocity of the mass point is calculated. We must use v1pt_theory because the mass point is moving away from the Origin. If this was a regular pendulum, we would use v2pt_theory. I will also display the velocity of point \(m\) in the \(A\) frame to show the work Sympy saved us.

# Set origin velocity to zero
O.set_vel(A, 0)

# Create the velocity of the mass point
P.set_vel(B, L_dot * B.y)
P.v1pt_theory(O, A, B)

\[L \dot{\theta}\mathbf{\hat{b}_x} + \dot{L}\mathbf{\hat{b}_y}\]

# Display the velocity of point m in the A frame
P.vel(A).express(A)

\[(L \dot{\theta} \operatorname{cos}\left(\theta\right) + \dot{L} \operatorname{sin}\left(\theta\right))\mathbf{\hat{a}_x} + (- L \dot{\theta} \operatorname{sin}\left(\theta\right) + \dot{L} \operatorname{cos}\left(\theta\right))\mathbf{\hat{a}_y}\]

I will be treating the mass as a particle and this is where an interesting property of the spring must discussed. The only object that has mass is point \(P\); the spring is massless. This is a little weird to think about because this is essentially an invisible force pulling on the mass. This is an example of one of the assumptions we will be making along with:

• The spring is perfect and cannot deform
• The support does not move whatsoever
• The mass is actually a point mass
• The pendulum only moves in the \(X\) and \(Y\) dimensions

pendulum = me.Particle('pend', P, m)

Now we create the forces that are acting on the particle: the spring force and gravity. We must think about how the spring force acts; when \(L\) is shorter than \(L_0\) the spring will be pushing on the mass. Since the pushing is in the positive direction of \(\mathbf{\hat{b}_y}\) the force is positive. Also just as a note: Sympy decapitalizes the frame names.

# Create damping forces for spring and rotation
spring = k * (L_0-L) * B.y
gravity = m * g * A.y
forces = spring + gravity
forces

\[k \left(L_{0} - L\right)\mathbf{\hat{b}_y} + g m\mathbf{\hat{a}_y}\] Now we create the KanesMethod object by suppling the inertial frame, the coordinates, the generalized velocities, and the kinematic differential equations.

kane = me.KanesMethod(A,
                      q_ind=[L, theta],
                      u_ind=[L_dot, theta_dot],
                      kd_eqs=[L_dot - L.diff(t),
                      theta_dot - theta.diff(t)])

Now we must create \(Fr\) and \(Fr^{*}\)

fr, frstar = kane.kanes_equations([(P, forces)], [pendulum])

The hard work is done and the equations of motion have been created so lets see what they look like.

M = kane.mass_matrix_full
f = kane.forcing_full
M.inv() * f

\[\left[\begin{matrix}\operatorname{\dot{L}}{\left (t \right )}\\dot{\theta}{\left (t \right )}\\frac{1}{m} \left(g m \cos{\left (\theta{\left (t \right )} \right )} + k \left(L_{0} - L{\left (t \right )}\right) + m L{\left (t \right )} \dot{\theta}^{2}{\left (t \right )}\right)\\frac{1}{m L^{2}{\left (t \right )}} \left(- g m L{\left (t \right )} \sin{\left (\theta{\left (t \right )} \right )} - 2 m L{\left (t \right )} \operatorname{\dot{L}}{\left (t \right )} \dot{\theta}{\left (t \right )}\right)\end{matrix}\right]\]

I will now be using PyDy to solve the ODE and plotting the data. You must first define the pydy.system with your KanesMethod object.

sys = System(kane)

Now define your constants and initial conditions.

sys.constants = {m:1.0,
                 g:9.81,
                 k:100.0,
                 L_0:1.0}
sys.initial_conditions = {L:1.0,
                      theta:np.deg2rad(55)}

Now create the time array, which is the length of time you want the simulation to run for. Also, though it’s not really necessary here, I create the ode function with cython. This makes the integration much faster for larger systems. If you do not have the cython package just omit that line of code and it will be solved using lamdify, the default solver. If you want more information on the backends check here.

runtime = 6
# For 30fps I want 30 datapoints per second so runtime is multiplied by 30
sys.times = linspace(0.0, runtime, runtime*30)
sys.generate_ode_function(generator='cython')
resp = sys.integrate()

sim_time = sys.times
fig = plt.figure(0, figsize=(12,5))
fig.add_subplot(121)
plt.plot(sim_time, resp[:,0], label='Unshaped')
plt.title(r'L', fontsize=18)
plt.xlabel('Time', fontsize=18)
plt.ylabel('Meters', fontsize=18)
fig.add_subplot(122)
plt.plot(sim_time, np.rad2deg(resp[:,1]), label='Unshaped')
plt.title(r'$\theta$', fontsize=18)
plt.xlabel('Time', fontsize=18)
plt.ylabel('Degrees', fontsize=18)
plt.tight_layout()

fig = plt.figure(0, figsize=(12,7))
plt.plot(resp[:,0]*np.sin(resp[:,1]),
         resp[:,0]*np.cos(resp[:,1]), label='Unshaped')
plt.gca().invert_yaxis()
plt.axes().set_aspect('equal', 'datalim')
plt.xlabel('Horizontal Motion')
plt.ylabel('Vertical Motion')

Going Further

There is a slight problem with this model, specifically due to the springs. When modeling with cables, wire, or rope there is somewhat of a breakdown in the spring assumption. Cables cannot push which leads to an interesting problem when constructing the equations of motion. I do not have the knowledge to solve such a problem by hand but SymPy can. By setting the forces as functions we can use a boolean operator to stop the spring from pushing. This will make the model behave more like a piece of string. The two functions will be (I also want damping to not push):

# L is the equilibrium length of the spring
def K(L1):
    k = 100
    L = 1
    return k * (L1 >= L)

def C(L1):
    c = 0.5
    L = 1
    return c * (L1 >= L)

Reading the function: if the supplied L1 is larger or equal to L then the function will return k or c. If L1 was smaller the function will return 0. Next we must alter the forces applied to include the functions:

# Create damping forces for spring and rotation
damping = -L_dot * C(L) * B.y
spring = K(L) * (L_0-L) * B.y
gravity = m * g * A.y
forces = damping + spring + gravity
forces

\[(100 \left(L_{0} - L\right) \left(L \geq 1\right) - 0.5 L_{dot} \left(L \geq 1\right))\mathbf{\hat{b}_y} + g m\mathbf{\hat{a}_y}\]

Finally in the sys.constants we no longer need to set k and c because they are defined in the functions.

The plots are slightly different as would be expected since there is not much spring compression and no rotational damping. Now the responses will be drastically different if the spring is first compressed by setting L = 0.5.

The no push model falls straight down much like what would be expected in real life. The push model behaves like a compressed spring and vibrates wildly as it swings.

This is a small sample of what Sympy is capable of.