GitHub

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

prescribed motion

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])
Example block output

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]")
Example block output