Prescribed movements
The movement of a frame can be prescribed using any of the components
The motion of joints Revolute
and Prismatic
can also be prescribed by using axisflange = true
and attaching a ModelingToolkitStandardLibrary.Rotational.Position
or ModelingToolkitStandardLibrary.TranslationalModelica.Position
to the axis flange of the joint.
Prescribed motions can be useful to, e.g.,
- Use the 3D rendering capabilities to visualize the movement of a mechanism moving in a prescribed way without simulating the system.
- Simplify simulations by prescribing the movement of a part of the mechanism.
In the example below, we prescribe the motion of a suspension system using Pose
. We load the Rotations
package in order to get access to the constructor RotXYZ
which allows us to specify a rotation matrix using Euler angles. The component prescribing the motion is
body_upright = Pose(; r = [0, 0.17 + 0.1sin(t), 0], R = RotXYZ(0, 0.5sin(t), 0))
This makes the body_upright
frame move in the vertical direction and rotate around the y-axis. Notice how the reference translation r
and the reference orientation R
are symbolic expressions that depend on time t
. To make a frame follow recorded data, you may use DataInterpolations.jl to create interpolation objects that when accessed with t
return the desired position and orientation.
using Multibody
using Multibody.Rotations # To specify orientations using Euler angles
using ModelingToolkit
using Plots
using OrdinaryDiffEq
using LinearAlgebra
using JuliaSimCompiler
using Test
t = Multibody.t
D = Differential(t)
n = [1, 0, 0]
AB = 146.5 / 1000
BC = 233.84 / 1000
CD = 228.60 / 1000
DA = 221.43 / 1000
BP = 129.03 / 1000
DE = 310.31 / 1000
t5 = 19.84 |> deg2rad
@mtkmodel QuarterCarSuspension begin
@structural_parameters begin
spring = true
(jc = [0.5, 0.5, 0.5, 0.7])#, [description = "Joint color"]
mirror = false
end
@parameters begin
cs = 4000, [description = "Damping constant [Ns/m]"]
ks = 44000, [description = "Spring constant [N/m]"]
rod_radius = 0.02
jr = 0.03, [description = "Radius of revolute joint"]
end
begin
dir = mirror ? -1 : 1
rRod1_ia = AB*normalize([0, -0.1, 0.2dir])
rRod2_ib = BC*normalize([0, 0.2, 0dir])
end
@components begin
r123 = JointRRR(n_a = n*dir, n_b = n*dir, rRod1_ia, rRod2_ib, rod_radius=0.018, rod_color=jc)
r2 = Revolute(; n=n*dir, radius=jr, color=jc)
b1 = FixedTranslation(radius = rod_radius, r = CD*normalize([0, -0.1, 0.3dir])) # CD
chassis = FixedTranslation(r = DA*normalize([0, 0.2, 0.2*sin(t5)*dir]), render=false)
chassis_frame = Frame()
if spring
springdamper = SpringDamperParallel(c = ks, d = cs, s_unstretched = 1.3*BC, radius=rod_radius, num_windings=10)
end
if spring
spring_mount_F = FixedTranslation(r = 0.7*CD*normalize([0, -0.1, 0.3dir]), render=false)
end
if spring
spring_mount_E = FixedTranslation(r = 1.3DA*normalize([0, 0.2, 0.2*sin(t5)*dir]), render=true)
end
end
begin
A = chassis.frame_b
D = chassis.frame_a
end
@equations begin
# Main loop
connect(A, r123.frame_a)
connect(r123.frame_b, b1.frame_b)
connect(b1.frame_a, r2.frame_b)
connect(r2.frame_a, D)
# Spring damper
if spring
connect(springdamper.frame_b, spring_mount_E.frame_b)
connect(b1.frame_a, spring_mount_F.frame_a)
connect(D, spring_mount_E.frame_a)
connect(springdamper.frame_a, spring_mount_F.frame_b)
end
connect(chassis_frame, chassis.frame_a)
end
end
@mtkmodel ExcitedWheelAssembly begin
@parameters begin
rod_radius = 0.02
end
@components begin
chassis_frame = Frame()
suspension = QuarterCarSuspension(; rod_radius)
wheel = SlippingWheel(
radius = 0.2,
m = 15,
I_axis = 0.06,
I_long = 0.12,
x0 = 0.0,
z0 = 0.0,
der_angles = [0, 0, 0],
iscut = true,
)
end
@equations begin
connect(wheel.frame_a, suspension.r123.frame_ib)
connect(chassis_frame, suspension.chassis_frame)
end
end
@mtkmodel SuspensionWithExcitationAndMass begin
@parameters begin
ms = 1500/4, [description = "Mass of the car [kg]"]
rod_radius = 0.02
end
@components begin
world = World()
mass = Body(m=ms, r_cm = 0.5DA*normalize([0, 0.2, 0.2*sin(t5)]))
excited_suspension = ExcitedWheelAssembly(; rod_radius)
prescribed_motion = Pose(; r = [0, 0.1 + 0.1sin(t), 0], R = RotXYZ(0, 0.5sin(t), 0))
end
@equations begin
connect(prescribed_motion.frame_b, excited_suspension.chassis_frame, mass.frame_a)
end
end
@named model = SuspensionWithExcitationAndMass()
model = complete(model)
ssys = structural_simplify(multibody(model))
display([unknowns(ssys) diag(ssys.mass_matrix)])
defs = [
model.excited_suspension.suspension.ks => 30*44000
model.excited_suspension.suspension.cs => 30*4000
model.excited_suspension.suspension.r2.phi => -0.6031*(1)
]
prob = ODEProblem(ssys, defs, (0, 2π))
sol = solve(prob, Rodas5P(autodiff=false), initializealg = BrownFullBasicInit())
@test SciMLBase.successful_retcode(sol)
Multibody.render(model, sol, show_axis=false, x=-0.8, y=0.7, z=0.1, lookat=[0,0.1,0.0], filename="prescribed_motion.gif") # Video
┌ Warning: Rosenbrock methods on equations without differential states do not bound the error on interpolations.
└ @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/55UVY/src/solve.jl:108
Even though we formulated an ODEProblem
and called solve
, we do not actually perform any simulation here! If we look at the mass matrix of the system
ssys.mass_matrix
6×6 LinearAlgebra.Diagonal{Int64, Vector{Int64}}:
0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 0
we see that it is all zeros. This means that there are no differential equations at all in the system, only algebraic equations. The solver will thus only solve for the algebraic variables using a nonlinear root finder. In general, prescribing the value of some state variables removes the need for the solver to solve for them, which can be useful for simplifying simulations. Using the "simulation" above, we can use the solution object to, e.g., find the compression of the spring and the forces acting on the ground over time.
plot(sol, idxs=[model.excited_suspension.suspension.springdamper.s, -model.excited_suspension.suspension.springdamper.f, model.excited_suspension.wheel.wheeljoint.f_n], labels=["Spring length [m]" "Spring force [N] " "Normal force [N]"], layout=(2,1), sp=[1 2 2])
Here, we see that the total spring force and the normal force acting on the ground are not equal, this is due to the spring not applying force only in the vertical direction. We can also compute the slip velocity, the velocity with which the contact between the wheel and the ground is sliding along the ground due to the prescribed motion.
wj = model.excited_suspension.wheel.wheeljoint
plot(sol, idxs=[wj.v_slip, wj.v_slip_long, wj.v_slip_lat], labels=["Slip velocity magnitude" "Longitudinal slip velocity" "Lateral slip velocity"], ylabel="Velocity [m/s]")