Elastic Pendulum

Simulation of an elastic pendulum combining spring and pendulum dynamics.

You can also find this example as a standalone Python file in the GitHub repository.

We implement this system im pathsim through two separate ODE blocks, one for the radial and one for the angular dynamics. The connections between the blocks represent the coupling through the forces.

import numpy as np
import matplotlib.pyplot as plt

# Apply PathSim docs matplotlib style
plt.style.use('../pathsim_docs.mplstyle')

from pathsim import Simulation, Connection
from pathsim.blocks import ODE, Function, Scope
from pathsim.solvers import RKBS32

# Initial conditions
r0, vr0 = 2, 0.0
phi0, omega0 = 0.3*np.pi, 0.0

# Physical parameters
g = 9.81          # gravity [m/s^2]
l0 = 1.0          # natural spring length [m]
k = 50.0          # spring constant [N/m]
m = 1.0           # mass [kg]
c_r = 0.3         # radial damping [kg/s]
c_phi = 0.1       # angular damping [N m s]

Next, the two ODE blocks are defined through their right hand side function. For the radial component this is the following set of odes with the centrifugal force coupling and the acceleration through gravity:

\[\begin{split}\dot{r} = v_r \\\end{split}\]

\[\dot{v_r} = -\frac{k}{m} (r - l_0) - \frac{c_r}{m} v_r + \omega^2 r + g \cos(\varphi)\]

# Define the radial ODE (spring-mass-damper with coupling terms)
def rad_ode(x, u, t):
    r, vr = x
    omega, phi = u

    # radial acceleration terms
    centrifugal = r * omega**2
    spring = -(k/m) * (r - l0)
    gravity_rad = g * np.cos(phi)
    damping = -(c_r/m) * vr

    accel_r = centrifugal + spring + gravity_rad + damping

    return np.array([vr, accel_r])

rad = ODE(rad_ode, np.array([r0, vr0]))

For the angular component we have the following odes:

\[\dot{\varphi} = \omega\]

\[\dot{\omega} = - \frac{g}{r} \sin(\varphi) - \frac{2}{r} v_r \omega - \frac{c_\varphi}{r^2 m} \omega\]

# Define the angular ODE (pendulum with coupling terms)
def ang_ode(x, u, t):
    phi, omega = x
    r, vr = u

    # angular acceleration terms
    gravity_torque = -(g / r) * np.sin(phi)
    coriolis = -(2 / r) * vr * omega
    damping = -(c_phi / (m * r**2)) * omega

    accel_phi = gravity_torque + coriolis + damping

    return np.array([omega, accel_phi])

ang = ODE(ang_ode, np.array([phi0, omega0]))

For straight forward plotting and evaluation it makes sense to convert from polar to cartesian coordinates. For this, we use a Function block. It can be defined through a decorator which makes it super compact:

# Cartesian conversion
@Function
def crt(r, phi):
    return r*np.sin(phi), -r*np.cos(phi)

# Scopes for visualization
sc1 = Scope(labels=["r [m]", "vr [m/s]", "phi [rad]", "omega [rad/s]"])
sc2 = Scope(labels=["x [m]", "y [m]"], sampling_period=0.005)

Next, define the connections (Connection) between the two blocks…

blocks = [rad, ang, crt, sc1, sc2]

connections = [
    Connection(ang[1], rad[0]),           # omega -> rad input 0
    Connection(ang[0], rad[1], crt[1]),   # phi -> rad input 1
    Connection(rad[0], ang[0], crt[0]),   # r -> ang input 0
    Connection(rad[1], ang[1]),           # vr -> ang input 1
    Connection(rad[:2], sc1[:2]),         # r, vr -> scope
    Connection(ang[:2], sc1[2:4]),        # phi, omega -> scope
    Connection(crt[:2], sc2[:2])
    ]

… initialize the Simulation with the blocks and connections. The dynamics of this system can be very chaotic, so it makes sense to use an adaptive timestep ode solver such as RKBS32:

# Simulation instance
sim = Simulation(
    blocks,
    connections,
    Solver=RKBS32
    )

11:39:40 - INFO - LOGGING (log: True)
11:39:40 - INFO - BLOCKS (total: 5, dynamic: 2, static: 3, eventful: 1)
11:39:40 - INFO - GRAPH (nodes: 5, edges: 9, alg. depth: 2, loop depth: 0, runtime: 0.040ms)

# Run the simulation
sim.run(10)

11:39:40 - INFO - STARTING -> TRANSIENT (Duration: 10.00s)
11:39:40 - INFO - --------------------   1% | 0.0s<0.4s | 4560.9 it/s
11:39:40 - INFO - ####----------------  20% | 0.1s<0.3s | 4768.9 it/s
11:39:40 - INFO - ########------------  40% | 0.2s<0.3s | 4424.6 it/s
11:39:40 - INFO - ############--------  60% | 0.3s<0.2s | 4659.4 it/s
11:39:40 - INFO - ################----  80% | 0.4s<0.1s | 4060.5 it/s
11:39:40 - INFO - #################### 100% | 0.5s<--:-- | 4696.3 it/s
11:39:40 - INFO - FINISHED -> TRANSIENT (total steps: 2130, successful: 2083, runtime: 488.38 ms)

{'total_steps': 2130,
 'successful_steps': 2083,
 'runtime_ms': 488.38450699986424}

One of pathsims convenient features is that each block is physicalized. This means the Scope block records the signals internally and the data can be read and plotted through its local methods directly:

# Plot state variables
sc1.plot();

We can clearly see the mass jumping back and forth, driventhrough the contracting spring at each swing:

# Plot the trajectory

t, (x, y) = sc2.read()

fig, ax = plt.subplots()
ax.plot(0, 0, "o", c="k")
ax.plot(x, y)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_aspect(1)