Flexible Double Pendulum¶

Two very flexible beam segments connected by a hinge joint, with the root mounted to a hinge

Imports¶

In [ ]:

from jax import numpy as jnp
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

from flapjax.models.double_pendulum.double_pendulum import generate_double_pendulum
from flapjax.structure import StructureCase

from jax import numpy as jnp from matplotlib import pyplot as plt from matplotlib.animation import FuncAnimation from IPython.display import HTML from flapjax.models.double_pendulum.double_pendulum import generate_double_pendulum from flapjax.structure import StructureCase

Case parameters¶

In [ ]:

dt = 0.005  # time step length
t_end = 5.0  # total simulation time
n_tstep = int(jnp.ceil(t_end / dt)) + 1
n = 10  # number of nodes per pendulum segment

dt = 0.005 # time step length t_end = 5.0 # total simulation time n_tstep = int(jnp.ceil(t_end / dt)) + 1 n = 10 # number of nodes per pendulum segment

Create model¶

Generate the double pendulum structure. The two segments are started oriented along the x-axis, connected by a hinge that allows rotation about the y-axis.

In [ ]:

beam, root_hinge, mid_hinge = generate_double_pendulum(
    n_nodes_per_segment=n,
    length=1.0,
    hinge_axis=jnp.array([0.0, 1.0, 0.0]),
    gravity=jnp.array([0.0, 0.0, -9.81]),
    k_bending=4.0,  # very low bending stiffness
)

beam, root_hinge, mid_hinge = generate_double_pendulum( n_nodes_per_segment=n, length=1.0, hinge_axis=jnp.array([0.0, 1.0, 0.0]), gravity=jnp.array([0.0, 0.0, -9.81]), k_bending=4.0, # very low bending stiffness )

Solve dynamic system¶

Both joints are hard constraints: a GroundedHinge pins the root to its initial position, and a MultibodyHinge connects the two segments. No DOFs are prescribed — the entire structure swings freely under gravity.

In [ ]:

init = beam.reference_configuration(prescribed_dofs=()).to_dynamic()

# set the initial beam acceleration to be gravity for all dynamic nodes
# this isn't essential, but ensures that at timestep 0 the forces are balanced
g = jnp.array([0.0, 0.0, -9.81])
init.v_dot = init.v_dot.at[:, :3].set(g[None, :])

# time march solution
solution: StructureCase = beam.dynamic_solve(
    init_state=init,
    n_tstep=n_tstep,
    dt=dt,
)

init = beam.reference_configuration(prescribed_dofs=()).to_dynamic() # set the initial beam acceleration to be gravity for all dynamic nodes # this isn't essential, but ensures that at timestep 0 the forces are balanced g = jnp.array([0.0, 0.0, -9.81]) init.v_dot = init.v_dot.at[:, :3].set(g[None, :]) # time march solution solution: StructureCase = beam.dynamic_solve( init_state=init, n_tstep=n_tstep, dt=dt, )

Animate¶

Animate the pendulum motion

In [ ]:

stride = n_tstep // 200
frames = range(0, n_tstep, stride)

fig_anim, ax_anim = plt.subplots(figsize=(6, 6))
(line1,) = ax_anim.plot([], [], "-", color="tab:blue", markersize=4, lw=2)
(line2,) = ax_anim.plot([], [], "-", color="tab:orange", markersize=4, lw=2)
(hinge_dots,) = ax_anim.plot([], [], "o", color="white", markersize=7, markeredgecolor="black", markeredgewidth=1.5,
                             zorder=5)
time_text = ax_anim.text(0.02, 0.95, "", transform=ax_anim.transAxes, fontsize=10)

pad = 0.1
x_all = solution.x[:, :, 0]
z_all = solution.x[:, :, 2]
ax_anim.set_xlim(float(x_all.min()) - pad, float(x_all.max()) + pad)
ax_anim.set_ylim(float(z_all.min()) - pad, float(z_all.max()) + pad)
ax_anim.set_aspect("equal")
ax_anim.set_xlabel("x [m]")
ax_anim.set_ylabel("z [m]")
ax_anim.set_title("Double Pendulum")
ax_anim.grid(True, alpha=0.3)


def update(i_ts):
    line1.set_data(solution.x[i_ts, :n, 0], solution.x[i_ts, :n, 2])
    line2.set_data(solution.x[i_ts, n:-1, 0], solution.x[i_ts, n:-1, 2])
    hinge_dots.set_data(
        [solution.x[i_ts, 0, 0], solution.x[i_ts, n - 1, 0]],
        [solution.x[i_ts, 0, 2], solution.x[i_ts, n - 1, 2]],
    )
    time_text.set_text(f"t = {float(i_ts * dt):.3f} s")
    return line1, line2, hinge_dots, time_text


anim = FuncAnimation(fig_anim, update, frames=frames, blit=True, interval=30)
plt.close(fig_anim)
HTML(anim.to_jshtml())

stride = n_tstep // 200 frames = range(0, n_tstep, stride) fig_anim, ax_anim = plt.subplots(figsize=(6, 6)) (line1,) = ax_anim.plot([], [], "-", color="tab:blue", markersize=4, lw=2) (line2,) = ax_anim.plot([], [], "-", color="tab:orange", markersize=4, lw=2) (hinge_dots,) = ax_anim.plot([], [], "o", color="white", markersize=7, markeredgecolor="black", markeredgewidth=1.5, zorder=5) time_text = ax_anim.text(0.02, 0.95, "", transform=ax_anim.transAxes, fontsize=10) pad = 0.1 x_all = solution.x[:, :, 0] z_all = solution.x[:, :, 2] ax_anim.set_xlim(float(x_all.min()) - pad, float(x_all.max()) + pad) ax_anim.set_ylim(float(z_all.min()) - pad, float(z_all.max()) + pad) ax_anim.set_aspect("equal") ax_anim.set_xlabel("x [m]") ax_anim.set_ylabel("z [m]") ax_anim.set_title("Double Pendulum") ax_anim.grid(True, alpha=0.3) def update(i_ts): line1.set_data(solution.x[i_ts, :n, 0], solution.x[i_ts, :n, 2]) line2.set_data(solution.x[i_ts, n:-1, 0], solution.x[i_ts, n:-1, 2]) hinge_dots.set_data( [solution.x[i_ts, 0, 0], solution.x[i_ts, n - 1, 0]], [solution.x[i_ts, 0, 2], solution.x[i_ts, n - 1, 2]], ) time_text.set_text(f"t = {float(i_ts * dt):.3f} s") return line1, line2, hinge_dots, time_text anim = FuncAnimation(fig_anim, update, frames=frames, blit=True, interval=30) plt.close(fig_anim) HTML(anim.to_jshtml())