add an example of gyroscopic effects (#88)

* add an example of gyroscopic effects

* set explicit world.g

* explicit g again

Author: baggepinnen

docs/make.jl

@@ -35,6 +35,7 @@ makedocs(;
                  "Ropes, cables and chains" => "examples/ropes_and_cables.md",
                  "Swing" => "examples/swing.md",
                  "Bodies in space" => "examples/space.md",
+                 "Gyroscopic effects" => "examples/gyroscopic_effects.md",
              ],
              "Rotations and orientation" => "rotations.md",
              "3D rendering" => "rendering.md",

docs/src/examples/gyroscopic_effects.md

@@ -0,0 +1,58 @@
+# Gyroscopic effects
+
+![animation](gyro.gif)
+
+In this example, we demonstrate how a rotating body creates a reaction torque when its axis of rotation is changed.
+
+The system consists of a pendulum suspended in a spherical joint, a joint without any rotational constraints. The tip of the pendulum is a cylinder that is rotating around a revolute joint in its center. When the pendulum swings, the rotation axis of the rotating tip is changed, this causes the entire pendulum to rotate around the axis through the pendulum rod.
+
+
+```@example spring_mass_system
+using Multibody
+using ModelingToolkit
+using Plots
+using JuliaSimCompiler
+using OrdinaryDiffEq
+
+t = Multibody.t
+D = Differential(t)
+world = Multibody.world
+
+systems = @named begin
+    spherical = Spherical(state=true, radius=0.02, color=[1,1,0,1])
+    body1 = BodyCylinder(r = [0.25, 0, 0], diameter = 0.05, isroot=false, quat=false)
+    rot = FixedRotation(; n = [0,1,0], angle=deg2rad(45))
+    revolute = Revolute(n = [1,0,0], radius=0.06, color=[1,0,0,1])
+    trans = FixedTranslation(r = [-0.1, 0, 0])
+    body2 = BodyCylinder(r = [0.2, 0, 0], diameter = 0.1, color=[0,0,0.5,1])
+end
+
+connections = [
+    connect(world.frame_b, spherical.frame_a)
+    connect(spherical.frame_b, body1.frame_a)
+    connect(body1.frame_b, rot.frame_a)
+    connect(rot.frame_b, revolute.frame_a)
+    connect(revolute.frame_b, trans.frame_a)
+    connect(trans.frame_b, body2.frame_a)
+]
+
+@named model = ODESystem(connections, t, systems = [world; systems])
+model = complete(model)
+ssys = structural_simplify(IRSystem(model))
+
+prob = ODEProblem(ssys, [model.world.g => 9.80665, model.revolute.w => 10], (0, 5))
+
+sol = solve(prob, Rodas5P(), abstol=1e-6, reltol=1e-6)
+@assert SciMLBase.successful_retcode(sol)
+using Test # hide
+@test sol(5, idxs=collect(model.body2.r_0[1:3])) ≈ [-0.0357364, -0.188245, 0.02076935] atol=1e-3 # hide
+# plot(sol, idxs=collect(model.body2.r_0)) # hide
+
+import CairoMakie
+Multibody.render(model, sol; x=1, z=1, filename = "gyro.gif") # Use "gyro.mp4" for a video file
+nothing # hide
+```
+
+![animation](gyro.gif)
+
+Try setting `model.revolute.w => 0` and plot this variable using `plot(sol, idxs=model.revolute.w)` and you will notice that the swinging of the pendulum induces a rotation around this joint, even if it has no rotational velocity from the start.

docs/src/examples/kinematic_loops.md

@@ -61,7 +61,7 @@ connections = [
 ]
 @named fourbar = ODESystem(connections, t, systems = [world; systems])
 m = structural_simplify(IRSystem(fourbar))
-prob = ODEProblem(m, [], (0.0, 30.0))
+prob = ODEProblem(m, [world.g => 9.81], (0.0, 30.0))
 sol = solve(prob, Rodas4(autodiff=false))
 
 plot(sol, idxs = [j1.phi, j2.phi, j3.phi])
@@ -72,7 +72,6 @@ using Test
 @test SciMLBase.successful_retcode(sol)
 @test sol(sol.t[end], idxs=j3.phi) % 2pi ≈ π/2 atol=0.3 # Test the the "pendulum" is hanging almost straight down after sufficient time has passed
 @test sol(sol.t[end], idxs=j2.phi) % 2pi ≈ -π/2 atol=0.3
-
 ```
 
 

src/Multibody.jl

@@ -27,10 +27,10 @@ Create a 3D animation of a multibody system
 
 # Camera control
 The following keyword arguments are available to control the camera pose:
-_ `x = 5`
-_ `y = 0`
-_ `z = 5`
-_ `lookat = [0,0,0]`: a three-vector of coordinates indicating the point at which the camera looks.
+- `x = 2`
+- `y = 0.5`
+- `z = 2`
+- `lookat = [0,0,0]`: a three-vector of coordinates indicating the point at which the camera looks.
 - `up = [0,1,0]`: A vector indicating the direction that is up.
 """
 function render end

src/components.jl

@@ -45,7 +45,7 @@ end
     # = 12 equations 
     @named frame_b = Frame()
     @parameters n[1:3]=[0, -1, 0] [description = "gravity direction of world"]
-    @parameters g=9.81 [description = "gravitational acceleration of world"]
+    @parameters g=9.80665 [description = "gravitational acceleration of world"]
     @parameters mu=3.986004418e14 [description = "Gravity field constant [m³/s²] (default = field constant of earth)"]
     @parameters render=render
     @parameters point_gravity = point_gravity
@@ -58,7 +58,7 @@ end
 end
 
 """
-The world component is the root of all multibody models. It is a fixed frame with a parallel gravitational field and a gravity vector specified by the unit direction `world.n` (defaults to [0, -1, 0]) and magnitude `world.g` (defaults to 9.81).
+The world component is the root of all multibody models. It is a fixed frame with a parallel gravitational field and a gravity vector specified by the unit direction `world.n` (defaults to [0, -1, 0]) and magnitude `world.g` (defaults to 9.80665).
 """
 const world = World(; name = :world)
 
@@ -151,7 +151,7 @@ Fixed translation followed by a fixed rotation of `frame_b` with respect to `fra
 - `sequence`: DESCRIPTION
 - `angle`: Angle of rotation around `n`, given in radians
 """
-@component function FixedRotation(; name, r, n = [1, 0, 0], sequence = [1, 2, 3], isroot = false,
+@component function FixedRotation(; name, r=[0, 0, 0], n = [1, 0, 0], sequence = [1, 2, 3], isroot = false,
                        angle, n_x = [1, 0, 0], n_y = [0, 1, 0])
     norm(n) ≈ 1 || error("n must be a unit vector")
     @named frame_a = Frame()
@@ -581,6 +581,8 @@ end
     @structural_parameters begin
         r = [1, 0, 0]
         r_shape = [0, 0, 0]
+        isroot = false
+        quat = false
     end
 
     @parameters begin
@@ -608,6 +610,7 @@ end
         density = 7700, [
             description = "Density of cylinder (e.g., steel: 7700 .. 7900, wood : 400 .. 800)",
         ]
+        color[1:4] = purple, [description = "Color of cylinder in animations"]
     end
     begin
         radius = diameter/2
@@ -641,7 +644,7 @@ end
         frame_a = Frame()
         frame_b = Frame()
         frameTranslation = FixedTranslation(r = r)
-        body = Body(; m, r_cm, I_11 = I[1,1], I_22 = I[2,2], I_33 = I[3,3], I_21 = I[2,1], I_31 = I[3,1], I_32 = I[3,2])
+        body = Body(; m, r_cm, I_11 = I[1,1], I_22 = I[2,2], I_33 = I[3,3], I_21 = I[2,1], I_31 = I[3,1], I_32 = I[3,2], isroot, quat)
     end
 
     @equations begin

test/runtests.jl

@@ -708,7 +708,7 @@ eqs = [connect(world.frame_b, freeMotion.frame_a)
 ssys = structural_simplify(IRSystem(model))
 @test length(unknowns(ssys)) == 12
 
-prob = ODEProblem(ssys, [collect(body.w_a .=> [0, 0, 0]); collect(body.v_0 .=> [0, 0, 0]); ], (0, 10))
+prob = ODEProblem(ssys, [world.g=>9.81; collect(body.w_a .=> [0, 0, 0]); collect(body.v_0 .=> [0, 0, 0]); ], (0, 10))
 
 sol = solve(prob, Rodas4())
 doplot() && plot(sol, idxs = body.r_0[2], title = "Free falling body")
@@ -732,7 +732,7 @@ eqs = [connect(world.frame_b, freeMotion.frame_a)
 ssys = structural_simplify(IRSystem(model))
 @test length(unknowns(ssys)) == 12 
 
-prob = ODEProblem(ssys, [collect(body.w_a .=> [0, 1, 0]); collect(body.v_0 .=> [0, 0, 0]); ], (0, 10))
+prob = ODEProblem(ssys, [world.g=>9.81; collect(body.w_a .=> [0, 1, 0]); collect(body.v_0 .=> [0, 0, 0]); ], (0, 10))
 
 sol = solve(prob, Rodas4())
 doplot() && plot(sol, idxs = body.r_0[2], title = "Free falling body")
@@ -752,6 +752,7 @@ ssys = structural_simplify(IRSystem(model))
 @test length(unknowns(ssys)) == 12
 
 prob = ODEProblem(ssys, [
+    world.g=>9.81; 
     collect(body.w_a .=> [0, 1, 0]); 
     collect(body.v_0 .=> [0, 0, 0]); 
 ], (0, 10))

test/test_orientation_getters.jl

@@ -32,7 +32,7 @@ end
 model = complete(model)
 ssys = structural_simplify(IRSystem(model))
 
-prob = ODEProblem(ssys, [model.shoulder_joint.phi => 0.0, model.elbow_joint.phi => 0.1], (0, 12))
+prob = ODEProblem(ssys, [model.shoulder_joint.phi => 0.0, model.elbow_joint.phi => 0.1, model.world.g => 9.81], (0, 12))
 sol = solve(prob, Rodas4())