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
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
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( 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
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.