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