Dynamic Free-Flying Structure¶

A very flexible beam subject to a time-dependent tip load at its base, causing large elastic deformation and rigid body motion. This script calculates the time history of the structural response.

Imports¶

In [ ]:

from jax import numpy as jnp
from matplotlib import pyplot as plt

from flapjax.models.flying_spaghetti.flying_spaghetti import generate_flying_spaghetti
from flapjax.structure import StructureCase

from jax import numpy as jnp from matplotlib import pyplot as plt from flapjax.models.flying_spaghetti.flying_spaghetti import generate_flying_spaghetti from flapjax.structure import StructureCase

Case parametrisation¶

In [ ]:

n_nodes = 21  # number of beam nodes
dt = 0.01  # time step length
t_end = 10.0  # time to run simulation for

n_nodes = 21 # number of beam nodes dt = 0.01 # time step length t_end = 10.0 # time to run simulation for

Create case¶

Create the data structure for the flying spaghetti model. This also generates an array for the external forces which is used in the canonical flying spaghetti case in "On the Dynamics in Space of Rods Undergoing Large Motions - A Geometrically-exact Approach", by Simo and Vu-Quoc, 1988.

For dynamic problems, we prescribe a spectral radius in range [0, 1] to control the damping of high-frequency content in the structural response in the time integrator. A value of 1 introduces no damping, while a value of 0 introduces considerable damping. Larger values are more likely to suffer from convergence issues, while smaller values may damp out the response too much. A value of 0.7 was found to be suitable for this case.

In [ ]:

# deduce number of timesteps from simulation time and time step
n_tstep = int(jnp.ceil(t_end / dt)) + 1

# time array for the simulation
t = jnp.arange(n_tstep) * dt

# create case and get external forces
# time-dependent forces f(t) have shape [n_tstep, n_nodes, 6]
struct, f_dead_2d, _ = generate_flying_spaghetti(
    n_nodes=n_nodes, t=t, spectral_radius=0.7
)

# deduce number of timesteps from simulation time and time step n_tstep = int(jnp.ceil(t_end / dt)) + 1 # time array for the simulation t = jnp.arange(n_tstep) * dt # create case and get external forces # time-dependent forces f(t) have shape [n_tstep, n_nodes, 6] struct, f_dead_2d, _ = generate_flying_spaghetti( n_nodes=n_nodes, t=t, spectral_radius=0.7 )

Solve dynamic system¶

Solve the structural dynamics problem using the dynamic_solve method. This takes the following inputs:

• init_state: The initial state of the system, for instance if initialising a dynamic system at a deformed state. If set to None, the initial state is set to the undeformed state.
• n_tstep: The number of time steps to solve for.
• dt: The time step length.
• f_ext_follower: The time-dependent nodal follower forces, [n_tstep, n_nodes, 6]. If set to None, no follower forces are applied.
• f_ext_dead: The time-dependent nodal dead forces, [n_tstep, n_nodes, 6]. If set to None, no dead forces are applied.
• f_ext_aero: The time-dependent nodal aerodynamic forces, [n_tstep, n_nodes, 6]. If set to None, no aerodynamic forces are applied. This is designed for use in coupled simulations with an aerodynamic solver.
• prescribed_dofs: An array of prescribed degrees of freedom index. If set to None, no degrees of freedom are prescribed, which is done in this free-flying case. This routine returns a Structure object which contains the time history of the structural response, which includes properties such as the nodal coordinates, orientations, velocities, and accelerations at each time step.

In [ ]:

solution: StructureCase = struct.dynamic_solve(
    init_state=None,
    n_tstep=n_tstep,
    dt=dt,
    f_ext_follower=None,
    f_ext_dead=f_dead_2d,
    f_ext_aero=None,
    prescribed_dofs=(),
)

solution: StructureCase = struct.dynamic_solve( init_state=None, n_tstep=n_tstep, dt=dt, f_ext_follower=None, f_ext_dead=f_dead_2d, f_ext_aero=None, prescribed_dofs=(), )

Plot results¶

In [ ]:

i_ts_plot = [0, 200, 400, 600, 800, 1000]  # extract the coordinates at these timestep index

# nodal x and z coordinates
x_out = solution.x[i_ts_plot, :, 0]  # [6, n_nodes]
z_out = solution.x[i_ts_plot, :, 2]  # [6, n_nodes]

# create plot over time
fig, ax = plt.subplots(figsize=(10, 4))
for i_snapshot in range(len(i_ts_plot)):
    ax.plot(x_out[i_snapshot, ...], z_out[i_snapshot, ...])  # plot beam
ax.set_xlabel("x")
ax.set_ylabel("z")
ax.set_title("Flying Spaghetti Trajectory")
ax.axis("equal")
plt.show()

i_ts_plot = [0, 200, 400, 600, 800, 1000] # extract the coordinates at these timestep index # nodal x and z coordinates x_out = solution.x[i_ts_plot, :, 0] # [6, n_nodes] z_out = solution.x[i_ts_plot, :, 2] # [6, n_nodes] # create plot over time fig, ax = plt.subplots(figsize=(10, 4)) for i_snapshot in range(len(i_ts_plot)): ax.plot(x_out[i_snapshot, ...], z_out[i_snapshot, ...]) # plot beam ax.set_xlabel("x") ax.set_ylabel("z") ax.set_title("Flying Spaghetti Trajectory") ax.axis("equal") plt.show()

Plot to Paraview¶

We can plot the 3D structure to Paraview using the plot method of the Structure object. This has the n_interp argument for interpolating to create mid-element nodes, as well as the index argument for selecting which time steps to plot. In this case, we plot every 10th timestep by computing the required timestep index. Two types of output file are written:

• .vtu files for each time step, which can be opened in Paraview to visualise the 3D structure.
• A .pvd file which contains the time stamps of the simulation, which can be opened in Paraview to animate the structure over time.

In [ ]:

pvd_path = solution.plot(directory="./flying_spaghetti_output", n_interp=3, index=jnp.arange(0, n_tstep, 10))

pvd_path = solution.plot(directory="./flying_spaghetti_output", n_interp=3, index=jnp.arange(0, n_tstep, 10))