Free motions
This example demonstrates how a free-floating Body
can be simulated. The body is attached to the world through a FreeMotion
joint, i.e., a joint that imposes no constraints. The joint is required to add the appropriate relative states between the world and the body. We choose enforceStates = true
and isroot = true
in the FreeMotion
constructor.
using Multibody
using ModelingToolkit
using Plots
using SymbolicIR
using OrdinaryDiffEq
t = Multibody.t
D = Differential(t)
world = Multibody.world
@named freeMotion = FreeMotion(enforceStates = true, isroot = true)
@named body = Body(m = 1)
eqs = [connect(world.frame_b, freeMotion.frame_a)
connect(freeMotion.frame_b, body.frame_a)]
@named model = ODESystem(eqs, t,
systems = [world;
freeMotion;
body])
ssys = structural_simplify(IRSystem(model))
prob = ODEProblem(ssys, [
D.(freeMotion.r_rel_a) .=> randn();
D.(D.(freeMotion.r_rel_a)) .=> randn();
D.(freeMotion.phi) .=> randn();
D.(D.(freeMotion.phi)) .=> randn();
D.(body.w_a) .=> randn();
], (0, 10))
sol = solve(prob, Rodas4())
plot(sol, idxs = body.r_0[2], title="Free falling body")
# Plot analytical solution
tvec = 0:0.1:sol.t[end]
plot!(tvec, -9.81/2 .* tvec .^ 2, lab="Analytical solution")
The figure indicates that the body is falling freely, experiencing a constant acceleration of -9.81 m/s² in the $y$ direction, corresponding to the gravity parameters of the world
:
Model world with 12 (18) equations
States (9):
(frame_b₊r_0(t))[1] [defaults to 0.0]: Position vector directed from the origin of the world frame to the connector frame origin, resolved in world frame
(frame_b₊r_0(t))[2] [defaults to 0.0]: Position vector directed from the origin of the world frame to the connector frame origin, resolved in world frame
(frame_b₊r_0(t))[3] [defaults to 0.0]: Position vector directed from the origin of the world frame to the connector frame origin, resolved in world frame
(frame_b₊f(t))[1] [defaults to 0.0]: Cut force resolved in connector frame
(frame_b₊f(t))[2] [defaults to 0.0]: Cut force resolved in connector frame
(frame_b₊f(t))[3] [defaults to 0.0]: Cut force resolved in connector frame
(frame_b₊tau(t))[1] [defaults to 0.0]: Cut torque resolved in connector frame
(frame_b₊tau(t))[2] [defaults to 0.0]: Cut torque resolved in connector frame
(frame_b₊tau(t))[3] [defaults to 0.0]: Cut torque resolved in connector frame
Parameters (4):
n[1] [defaults to 0]: gravity direction of world
n[2] [defaults to -1]: gravity direction of world
n[3] [defaults to 0]: gravity direction of world
g [defaults to 9.81]: gravitational acceleration of world