Tutorial: The Cubecopter

In this page we model the movement of a simple deformable body on a reference frame moving with it and then solve the problem to find its trajectory from the point of view of an inertial frame.

The example body will be a cube attached to rotating helix, from now on comradely referred to as the cubecopter. Each particle is represented by a PointMass, and the full body is an array of point masses.

r_0 = [
      PointMass(1., [ 1.,  1.,  0.])
    , PointMass(1., [ 1., -1.,  0.])
    , PointMass(1., [-1.,  1.,  0.])
    , PointMass(1., [-1., -1.,  0.])
    , PointMass(1., [ 1.,  1., -1.])
    , PointMass(1., [ 1., -1., -1.])
    , PointMass(1., [-1.,  1., -1.])
    , PointMass(1., [-1., -1., -1.])
    , PointMass(.5, [-1.,  0.,  .5])
    , PointMass(.5, [ 1.,  0.,  .5])
    ]

For visualization purposes, the plot function from Plots.jl is overloaded to work with arrays of point masses.

using Plots: plot
plot(r_0)

It may be hard to visualize a lot of points scattered on a 3D graph. Therefore, the plot function is also overloaded to accept a parameter bodylines connecting some particles. It does not interferes on the dynamics and serves only for aiding visualization.

Below is a nasty hack for calculating the cubecopter's edges. You do not have to understand it. The important part is that the variable edges is an array of indexes representing which points should be connected.

edges = Tuple[]
for i = 1:length(r_0), j = i:length(r_0)
    if count(a -> first(a) == last(a), zip(pos(r_0[i]),pos(r_0[j]))) == 2
        push!(edges, (i,j))
    end
end

plot(r_0, bodylines=edges)

Much better, right?

Now, the array r_0 represents a stationary body. What we want is a trajectory over time, represented by an array of functions. Our chosen trajectory will be a stationary cube with an helix rotating counter-clockwise over it.

For the cube part, we represent these as constant functions.

trajectory = Function[]
for x in r_0[1:end-2]
    push!(trajectory, let x=x; t -> x;end)
end

While the helix movement is done using the function rotate over a fixed axis with an angle varying over time.

const ω = 2*π/5.0            # Angular velocity
const z_axis  = [0., 0., 1.] # Axis of rotation is orthogonal to helix

for x in r_0[end-1:end]
    push!(trajectory, let x=x; t -> rotate(x, axis=z_axis, angle=ω*t); end)
end

With this, we complete the code for the trajectory on the body frame. Let's animate it so we can see how the trajectory behaves.

using Plots: Animation, frame, gif
anime = Animation()
for t in 0.:0.1:7.
    frame(anime, plot([ x(t) for x in trajectory], bodylines=edges))
end
gif(anime)

It is now time to define the problem's model and find the trajectory on the inertial frame. First, we need to define some initial data.

tstart  = 0.                # Starting time
tend    = 5.                # Ending time
Rstart  = one(Quaternion)   # Initial rotation
L_cm    = zeros(3)          # Angular momentum as viewed from center of mass

The variables tstart and tend define respectively the starting and ending time for the dynamics. The other two variables are the initial data necessary to solve the differential equation. The term one(Quaternion) is our way to represent the identity rotation, meaning that at the starting time, the inertial frame coincides with the body frame. The variable L_cm is set to zero, meaning that the total angular momentum is zero from the point of view of the center of mass.

With this information, we are ready to define our Model.

model = Model( trajectory
             , timespan = (tstart, tend)
             , q_0  = Rstart
             , L_cm = L_cm
             )

Now that everything is setted up, finding the inertial frame trajectory is as simple as writing

solve!(model)

To visualize the final result, we can use the function plotmodel to plot and save an animation of the dynamics of both frames side by side.

plotmodel( model
         , :both
         , saveas="cubecopter.gif"
         , bodylines=edges          # Only for visualization
         , duration=5.0             # Duration in seconds
         )

In the inertial frame, the cube rotates in the opposite direction to the helix guaranteeing that the total angular momentum is conserved.