Analyze System Dynamics in MATLAB

This example shows how to analyze the dynamics of a four-bar linkage system using the MATLAB classes of Simscape Multibody. The example system is a crank-rocker mechanism that has four links measuring 5 cm, 8 cm, 10 cm, and 12 cm in length.

The image shows the initial configuration of the system, where the angle between the crank link and the base link is 60 degrees. The crank link drives the system with an angular velocity of 10 degrees per second. Also, an angular acceleration of 5 deg/s² propels the crank link, while a load of 50 N acts on the distal end of the rocker link. This force aligns along the x-axis of the world frame.

Create Four-Bar System

To construct the four-bar multibody system, open the example and use the custom function fourbar with the specified link lengths. For detailed instructions on building a four-bar system in MATLAB, see Model a Multibody System in MATLAB.

import simscape.multibody.* simscape.Value simscape.op.*
[fb,op] = fourBar(simscape.Value(12,"cm"), simscape.Value(10,"cm"), simscape.Value(5,"cm"), simscape.Value(8,"cm"));

The four-bar system operates within the X-Y plane, with gravity acting in the -y direction.

Before performing analyses on a multibody system, compile the associated multibody object.

cfb = compile(fb);

To set the initial joint primitive targets, create an operating point object, jointOP.

jointOP = simscape.op.OperatingPoint;
jointOP("Bottom_Left_Joint/Rz/q") = simscape.op.Target(60,"deg","High");
jointOP("Bottom_Right_Joint/Rz/q") = simscape.op.Target(90,"deg","Low");
jointOP("Bottom_Left_Joint/Rz/w") = simscape.op.Target(10,"deg/s","High");

To verify that the four-bar system meets all the specified targets, compute the state of the system by using the computeState method.

state = computeState(cfb,jointOP)

state = 
  State:

  Status: Valid

  Assembly diagnostics:
  x
    Bottom_Left_Joint
        'Bottom_Left_Joint' successfully assembled
        Rz
            High priority position target +60 (deg) achieved
            High priority velocity target +10 (deg/s) achieved
    Top_Left_Joint
        'Top_Left_Joint' successfully assembled
        Rz
            Free position value: +38.485 (deg)
            Free velocity value: +12.7749 (deg/s)
    Bottom_Right_Joint
        'Bottom_Right_Joint' successfully assembled
        Rz
            Low priority position target +90 (deg) not achieved; actual value: +91.4095 (deg)
            Free velocity value: +4.14183 (deg/s)
    Top_Right_Joint
        'Top_Right_Joint' successfully assembled
        Rz
            Free position value: +69.8945 (deg)
            Free velocity value: +6.91669 (deg/s)

The results indicate that the system meets all the specified targets.

To visualize the system, use the visualize method:

visualize(cfb,computeState(cfb,jointOP),"vizFourBar")

Specify Joint Internal Mechanics

In real-life applications, the joints in a mechanical system are not ideal and often experience energy losses. To account for these phenomena in the system, specify the spring and damping properties by using the subclasses of the simscape.multibody.JointForceLaw class.

To represent the torsional spring-damper law, create a simscape.multibody.TorsionalSpringDamper object.

tsd = simscape.multibody.TorsionalSpringDamper;

To set the equilibrium position and stiffness of the spring, use a simscape.Value object.

tsd.EquilibriumPosition = simscape.Value(90,"deg");
tsd.SpringStiffness = simscape.Value(0.5,"N*cm/deg")

tsd = 
  TorsionalSpringDamper with properties:

    EquilibriumPosition: 90 (deg)
        SpringStiffness: 0.5000 (N*cm/deg)
     DampingCoefficient: 0 (N*m*s/deg)

To specify the damping coefficient of the damper, use a simscape.Value object.

tsd.DampingCoefficient = simscape.Value(0.1,"N*cm/(rev/s)")

tsd = 
  TorsionalSpringDamper with properties:

    EquilibriumPosition: 90 (deg)
        SpringStiffness: 0.5000 (N*cm/deg)
     DampingCoefficient: 0.1000 (N*cm*s/rev)

For simplicity, apply the same spring-damper law to all joint primitives.

Bottom_Left_Joint.Rz.ForceLaws = tsd;
Top_Left_Joint.Rz.ForceLaws = tsd;
Top_Right_Joint.Rz.ForceLaws = tsd;
Bottom_Right_Joint.Rz.ForceLaws = tsd;

Model Joint Actuation and External Loads

As shown in the system diagram, the crank link drives the system at an angular acceleration of 5 deg/s^2, while a load of 50 N acts on the distal end of the rocker link. To construct the acceleration and load for the system, use the simscape.multibody.JointAccelerationDictionary and simscape.multibody.ExternalForceTorqueDictionary classes.

To pair and store the joint acceleration with the associated joint primitive, use the dict_accel dictionary.

dict_accel = simscape.multibody.JointAccelerationDictionary;
accel_revolute = simscape.multibody.RevolutePrimitiveAcceleration(simscape.Value(5,"deg/s^2"));
dict_accel("Bottom_Left_Joint/Rz") = accel_revolute;

By default, all joint primitives in the system have zero actuation torque. However, the bottom left joint undergoes an angular acceleration. To maintain balance during dynamics computation, you must specify one joint primitive to automatically compute torque. For example, you can specify the joint primitive between the connector link and the rocker link as automatically computed by using the simscape.multibody.JointActuationDictionary class.

dict_t = simscape.multibody.JointActuationDictionary;
torque = simscape.multibody.RevolutePrimitiveActuationTorque("Computed");
dict_t("Top_Right_Joint/Rz") = torque;

To pair and store the load with the associated frame connector, use the dict_f dictionary. The load aligns along the x-axis of the world frame.

dict_f = simscape.multibody.ExternalForceTorqueDictionary;
force = simscape.multibody.ExternalForce(simscape.Value([50 0 0],"N"),"World");
dict_f("Right_Link/pos_end/f") = force;

Compute Dynamics

To compute the dynamics of the four-bar system using the given state and actuation inputs, use the computeDynamics method.

dynamicsResult = computeDynamics(cfb,state,dict_accel,dict_t,dict_f);

You can now query the positions, velocities, and accelerations of every joint primitive in the system. For example, to obtain the angular velocity and acceleration of the joint primitive between the base link and the rocker link, use the jointPrimitiveVelocity and jointPrimitiveAcceleration methods.

vel_bottom_right = jointPrimitiveVelocity(cfb,"Bottom_Right_Joint/Rz",state)

vel_bottom_right = 
    4.1418 (deg/s)

accel_bottom_right = jointPrimitiveAcceleration(cfb,"Bottom_Right_Joint/Rz",dynamicsResult)

accel_bottom_right = 
    3.0495 (deg/s^2)