8. Inverse Dynamics

## Overview

This chapter covers MechanicalSystems basic dynamic modeling capabilities for fully constrained models. If mass properties are defined for the bodies in a model, Mech can calculate the inertial loads on the model and the resulting reaction forces. Because the inertial loads are a function of the velocity and acceleration of the bodies, a model must have some dependence on time in order to have nonzero inertial loads. The forces that result from an inverse dynamic analysis are represented in the same manner as those of a static analysis, in the form of Lagrange multipliers (generalized forces). The free response dynamics of underconstrained systems is covered in Chapter 10.

# 8.1 Inertia Properties

In order for any inertial loads to exist in a mechanism model, the bodies of the mechanism must have nonzero inertia properties, mass and inertia. The inertia properties of the bodies in a model are defined in Mech body objects, and then added to the current model with the SetBodies command.

## 8.1.1 Body Options

Mech body objects are used to define specific properties of a body in a mechanism model. Properties are defined with options to the Body function.

The body object builder.

Two of the options for Body that are used for basic kinematic models are PointList and InitialGuess, which are discussed in Section 2.1. The remaining three options for Body are used to define the inertia properties of a body.

### 2D

2D inertia property options for Body.

### 3D

3D inertia property options for Body.

Executing the Body function does not, in itself, have any affect on the current Mech model. The returned SysBody data objects must be added to the model with the SetBodies command.

The body object processing function.

## 8.1.2 Examples

To use the 2D example mechanism of Section 1.2 in a dynamic analysis, inertia properties would be added with body objects and the SetBodies function. The crankshaft-piston mechanism used three Modeler2D bodies: the ground, crankshaft, and piston.

For this example, it is assumed that the constraints of the crankshaft-piston model have already been defined as they were in Section 1.2.
There is no use in adding inertia properties to the ground body, as it cannot move, and therefore can sustain no inertial loads. The inertia properties of the crankshaft and piston are defined and added to the model.

Because the default setting for all options for Body is Automatic, which leaves the current values of properties unchanged, a single inertia property can be changed without affecting the rest of the model. An example is changing the location of the centroid of the piston without changing any other body properties.

Each of the inertia property options for Body can be interrogated to find the current settings of a body property.

Inertia property functions.

## 8.1.3 Compound Inertia Properties

The centroidal inertia properties of a 3D body in a mechanism model must be determined before the dynamic reaction forces in a model can be found. For complex 3D bodies, the values of the moments and products of inertia can be quite difficult to determine manually. This problem is alleviated by the availability of 3D solid modeling systems that can generally calculate the inertia properties of any solid body that they are capable of modeling.
If the inertia properties of a body must be found manually, the body can be subdivided into several elementary geometrical components such as spheres, cylinders, or plates. Modeler3D provides several functions to aid in manipulating the inertia properties of a body or multiple bodies. These functions can make the task of determining the inertia properties of a slightly complex body more manageable.

Utilities for manipulating 3D inertia properties.

### Inertia Transformations

The following example takes the inertia properties of a triangular plate that are calculated in a noncentroidal reference frame and convert them to a centroidal frame. The plate is an isosceles right triangle with two one-unit sides.

The inertia properties of the plate, in the shown reference frame, can be integrated directly or found in an engineering or mathematics handbook.

The location of the centroid of the plate can also be found by direct integration, knowing that the mass of the plate is 1/2 because the density of the plate is one unit mass per unit area.

To use the ParallelAxis function to transform noncentroidal inertia properties into a centroidal reference frame, instead of vice versa, the vector from the centroid to the origin of the noncentroidal frame is specified in the centroidal reference frame coordinates, and a negative mass is specified. Note that the sign of the centroid vector is irrelevant, but it must be given in the correct reference frame if there is any relative rotation between the frames.
ParallelAxis returns the centroidal inertias {Ixx, Iyy, ...} that are reduced in magnitude relative to the noncentroidal properties.

The principal axes of inertia of the plate can be found with the PrincipalAxes function. Note that the list of the three principal axes of inertia forms the rotation matrix from the local reference frame of the inertia properties to the principal reference frame.

### Composite Bodies

The CompositeInertia function is used to add the inertia properties of several bodies together and return the centroidal inertia properties of their superset. In the following example, the inertia properties of the triangular plate are added to the inertia properties of a cylinder that is attached to the plate and rotated through a 30-degree angle about the local x axis, as shown below. The radius of the cylinder is 0.25 and the length is 1 unit. The bottom surface of the cylinder is coplanar with the local x axis of the plate.

The centroidal inertia properties of a cylinder in an axial reference frame can be found in an engineering manual. These properties have been determined and are given in the following example.

CompositeInertia can now be used to combine the inertia properties of the plate and the cylinder into one set of centroidal inertia properties. Note that the optional rotation argument is used for the cylinder since the reference frame in which its inertia properties were calculated was not aligned with the local reference frame.

### 2D Inertia Functions

Modeler2D also provides two inertia transformation functions that parallel those in Modeler3D, although the calculations they perform are trivial by comparison.

Utilities for manipulating 2D inertia properties.

# 8.2 Kinematic and Dynamic Solutions

## 8.2.1 Dynamics Options

Once the kinematic constraints and the inertia properties of a body have been defined, the dynamic forces in a mechanism can be determined. The Solution option to SolveMech is used to tell Mech to calculate the generalized reaction forces in the model and to include the forces due to inertia.

An option for SolveMech.

The purpose of the Kinematic setting is to allow for velocity-dependent loading, such as damping, while still not including the dynamic loading. Note that using the Static setting does not necessarily produce an error in a model that has velocity-dependent loading, because Static sets the velocities to zero. This condition may result in an incidental error, however, if the velocity-dependent terms happen to blow up at zero velocity.
The Dynamic setting calculates the mechanism accelerations and then adds the inertial loads due to the acceleration of each body to the model's applied load vector before calculating the generalized reaction forces. Note that the Dynamic setting must be used to include the centrifugal forces, even though they depend on velocity and not acceleration.

## 8.2.2 3D Example Mechanism

A Modeler3D example model is developed to demonstrate a simple gyroscopic inertial loading condition. The model is a flywheel with an angular velocity vector with a changing direction.
The model consists of only one moving body, the flywheel. The flywheel spins about an axis that is attached to the ground at one point. The axis of rotation of the flywheel is then rotated about another nonparallel axis affixed to the ground. Thus, the direction of the axis of rotation of the flywheel changes with time.

The Method -> Angular option is given to SetSymbols before building any other part of the flywheel model. This is done so that the angular velocity and angular acceleration of the flywheel is calculated directly, instead of being calculated in terms of derivatives of Euler parameters. Note that this is not a necessary step to obtain the dynamic forces.

### Inertia

The flywheel is modeled by four one-unit point masses, symmetrically placed about the axis of rotation. The Modeler3D CompositeInertia function can be used, trivially, to find the net inertia properties of the flywheel.

### Bodies

One body object is defined for the flywheel model. The ground body does not have a body object because it has no inertia or required point definitions. The body object for the flywheel contains the inertia properties and the locations of each of the four point masses.

### Constraints

Three constraints, two of which are functionally dependent on time, are required to model the flywheel mechanism.
A Spherical3 constraint places the center of the flywheel at the global origin.
A Parallel2 constraint forces the axis of rotation of the flywheel to be parallel with a specified vector on the ground body. The direction of the axis vector on the ground body is explicitly a function of time T. This constraint effectively causes the axis of rotation of the flywheel to change direction with time.
A ProjectedAngle1 constraint rotates the flywheel about its axis. Note the choice of the vector that was used on the ground body. The vector in the Y direction on the ground body is always normal to the axis of rotation of the flywheel. If a vector in the X direction was chosen instead, the constraint would become singular at T = 0.25. Note also that the flywheel is spinning 10 times more rapidly about its symmetry axis than its axis is rotating with respect to the ground.

### Runtime

Because of the presence of T in the constraints, the model can be run through its intended range of motion by varying T directly with the SolveMech command.

No forces are applied to the model with SetLoads because only the inertial loading of the flywheel is being considered. All of the inertia properties have already been defined, so all that must be done to find the dynamic solution is to request it of SolveMech with the solution option.

The angular velocity and acceleration of the flywheel can be read directly from the solution rules returned by SolveMech, or Modeler3D output functions can be used. Note that the angular acceleration is orthogonal to the angular velocity.

The moment that is applied to the flywheel to enforce the specified rotation is the reaction force to constraint 2, the constraint that controls the direction of the axis of the flywheel.
Note that the reaction to constraint 3 is zero. No moment is required to spin the flywheel about its primary axis, only to change the direction of the axis.

As the position of the flywheel is advanced in time, the components of the angular velocity change, and the components of the angular acceleration and the applied moment follow.

# 8.3 Friction and Damping

This section demonstrates various techniques of applying friction and damping forces to a mechanism model with MechanicalSystems.

## 8.3.1 Types of Friction

Loads that can be applied to a Mech model fall into four basic categories with considerable overlap between them.
An applied load is a load in which the magnitude and/or direction depends on the configuration of the model, on time, or is constant. Applied loads may or may not be conservative. Loads of this type are discussed throughout Chapter
7.
A damping load is a nonconservative load in which the magnitude and/or direction depends on the velocity of the model. These loads are modeled by giving an expression for the load magnitude or direction that is functionally dependent on velocity terms, and then using the Kinematic or Dynamic option settings for SolveMech.
A frictional load is a nonconservative load in which the magnitude and/or direction depends on the reaction forces applied by the constraints. These loads are modeled by giving an expression for the load magnitude or direction that is functionally dependent on the Lagrange multipliers, usually by using the Reaction function to find the reaction force at a constraint.
An inertial load is a load in which the magnitude and/or direction depends on the linear and angular acceleration and angular velocity of bodies in the model. Inertial loads are applied to the model automatically by specifying inertia properties for each body, and then using the Dynamic option setting for SolveMech, as discussed in Section 8.2.

### Damping

Damping loads are usually applied to a model in a relatively simple manner. To apply linear damping to a point on a moving body, a force is applied to the body that is inversely proportional to the velocity of the point. The magnitude of the force is proportional to the magnitude of the velocity of the point (relative to whatever the point is in contact with) and the direction vector of the force is opposite the velocity vector of the point.

The velocity of a point.

The following example shows how damping with a damping constant of 4.0 would be applied to a point at local coordinates {1, 0, 0} on body 2 in a 3D model.

The line of action of the force passes through Point[2, {1,0,0}], which is not dependent on velocity, and points in the direction of the velocity of the point. The negation of the velocity vector causes the applied force to oppose the motion of the point, and the Magnitude -> Relative setting causes the magnitude of the force to be 4.0 times the magnitude of the velocity vector, instead of a constant magnitude of 4.0.
To model point damping that occurs at an interface between two moving bodies, the velocity of the point to which damping is applied must be known relative to the body that the point is in contact with. The Mech RelativeVelocity function returns this vector quantity.

The velocity of a point relative to another body.

The following example shows how damping that exists between body 2 and body 3, with a damping constant of 4.0, would be applied to body 2 at local coordinates {0, 1, 0} in a 3D model.

Note that the preceding example only applies the frictional force to body 2. An equal and opposite force should also be applied to body 3.

### Friction

A Coulombic friction load is functionally dependent on reaction forces in the model. Such a load is usually applied by making the magnitude of the load a function of reaction forces by using the Mech Reaction function.

The reaction load from a constraint.

Coulombic friction forces may also be dependent on velocity terms to establish the direction of the frictional force.
The following example shows how the friction of a point sliding on a plane might be modeled. The model uses a PointOnPlane1 constraint to place a point on body 2 in contact with a plane on body 3, and a coefficient of friction between the point and the plane is 0.3. The frictional force would be modeled in the general case by using the velocity vector of the point, relative to the plane, to define the direction of the force, and the reaction to the constraint times the coefficient of friction to define the magnitude.
The second constraint in the following SetConstraints call is just a dummy constraint to fill out the needed 12 degrees of freedom.

The optional second argument of Force is used here to apply a force to body 2 and an equal and opposite force to body 3.

## 8.3.2 Damping Example

To demonstrate the application of damping to a Mech model, a 2D model of a simple four-bar linkage is developed. The model consists of three moving bodies: the drive bar, the driven bar, and the center bar. The kinematic input to the model is the rotation of the drive bar.
Damping is first applied to the model at a point on the center bar that traces out an eccentric path in the 2D plane. Coulombic friction is then applied to the rotational axis of the drive bar.

### Bodies

Four body objects are required for the four-bar model. Each bar is attached to another bar at its local origin and at a point at the other end of the bar on the local X or Y axis.

### Constraints

Five constraints, one of which is a driving constraint, are required to model the four-bar mechanism. A RotationLock1 constraint is used as a driving constraint to rotate the drive bar. Four Revolute2 constraints model the four pivot points in the four-bar mechanism.

A damping force is applied to point 2 on the center bar relative to the ground. The direction of the frictional force opposes the global velocity of point 2, and the coefficient of damping is 5.0.

### Runtime

Because of the presence of T in the driving constraint cs[1], the model can be run through its intended range of motion by varying T directly with the SolveMech command. The numerical value of T specifies the angle of rotation of the drive bar in revolutions.

Because the loading in this model is dependent on velocity terms, the Solution -> Kinematic option must be used with SolveMech to solve for the reaction forces. If Solution -> Static were used, SolveMech would assume that all velocities were zero, thus, the applied damping would have zero magnitude and zero reactions would result.

### Plots

To generate a plot of the torque required to turn the drive bar, the four-bar model is run from T = 0 to T = 1 in 21 steps, a full turn of the drive bar, while requesting that SolveMech interpolate the results.

A parametric plot of the x-y reaction at the axis of the drive bar can be easily created. The plot ranges only from T = 0 to T = 0.95 so the start and end points are apparent.

A very simple Coulombic friction load is applied to the rotational axis of the drive bar, with respect to the ground. The load is a pure moment resisting the direction of motion, the magnitude of which is proportional to the magnitude of the reaction vector at the drive bar axis. A coefficient of friction of 0.5 is assumed and is applied to a journal bearing with radius 0.4.

### More Plots

The drive torque function that was developed earlier can be plotted again with the new solution rules. Note that the required drive torque is slightly higher.

## 8.3.3 Friction Example

To demonstrate the application of Coulombic friction to a Mech mechanism model, a 2D model of a slider on an inclined ramp is developed. The model consists of one moving body, the slider.
As the slider is pushed along its track on the ground body, it travels along the initial horizontal part of the track until the tip of the slider encounters the upward turn of the track. The tip then progresses up the inclined ramp while the heel of the slider stays on the horizontal section. Coulombic friction is modeled at the two sliding interfaces at the tip and the heel of the slider. The Modeler2D TimeSwitch function is used to switch between the different friction forces as the tip of the slider passes from one section of the ramp to the next.

### Bodies

Two body objects are defined for the slider-ramp model.

• The ramp (body 1) requires five point definitions.
P1 is the origin of the ramp.
P2 is a point directly to the right of the origin to define the lower flat of the ramp.
P3 is the center point of the curved section of the ramp.
P4 and P5 are two points that define the upper flat of the ramp.
• The slider (body 2) requires two local point definitions.
P1 is the center of the heel of the slider, which is the local origin.
P2 is the center of the tip of the slider.
• Note that points 3, 4, and 5 on the ramp are defined as functions of two user-defined variables radius and incline, so that the radius at the curve of the ramp and the angle of the incline can be changed easily.
Because a gravitational force will be applied to the slider, the mass and centroid location are defined in the slider's body object.

### Constraints

Three constraints, one of which is a driving constraint, are required to model the slider-ramp mechanism.

• A RelativeX1 constraint controls the horizontal position of the heel of the slider. This constraint is a driving constraint because it is functionally dependent on the time variable T.
• A RelativeY1 constraint keeps the heel of the slider in contact with the ramp.
• A TimeSwitch constraint keeps the tip of the slider in contact with the ramp. The TimeSwitch constraint has three stages:
1. A PointOnLine1 constraint keeps the tip of the slider in contact with the horizontal section of the ramp.
2. A RelativeDistance1 constraint keeps the tip of the slider in contact with the curved section of the ramp.
3. A PointOnLine1 constraint keeps the tip of the slider in contact with the inclined section of the ramp.
• The value of the initial switch time determined by inspection is 0.5. However, the value of the second switch time is more difficult to determine analytically, so it is left as the symbol swtime and the FindSwitchTime function is used to find a numeric value.

### Runtime

Before the model can be run the undefined symbols, incline, radius, and swtime, must be defined. \$MechVariables returns the names of any symbols that have been recognized as being part of the current model.

### Switch Time

Now the exact value of swtime can be found with FindSwitchTime. Because the slider passes tangentially from the curved section of the ramp onto the inclined section (which implies that there is no discontinuity in the velocity of the slider at the transition), the TangentConstraint option must be used with FindSwitchTime. This option tells FindSwitchTime to look for the point of velocity equality instead of just location equality.

The goal of this example is to determine the force required to push the slider around the curved section of the ramp and up the incline, with gravity pulling down on the slider, and friction between the slider and the ramp at both ends.
First, gravity is applied without any friction forces. A plot of the force required to push the slider is made from a list of solutions returned by SolveMech.

Now a friction force is added to the left end of the slider. This is fairly straightforward, as the line of action of the force is constant throughout the motion of the slider. The force is applied to the slider and the line of action is along the surface of the ramp. The magnitude of the force is equal to the coefficient of friction times the normal force at the left end of the slider, which is the Y reaction to constraint 2.

Now a friction force is added to the right end of the slider. This is rather complicated, as the line of action of the friction force changes as the slider passes from one section of the ramp to another.
The friction force is applied to the slider, but the line of action is along the surface of the ramp. On the two flat sections of the ramp, the friction force is applied in the same manner as the force on the left end of the slider. However, the line of action of the force around the curve on the ramp has to be developed directly with vector algebra.
Although the line of action changes as the slider moves through the different stages of the ramp, the magnitude of the normal force is always equal to the coefficient of friction times the magnitude of the reaction force of constraint 3, even though constraint 3 is changing from one constraint type to another.
TimeSwitch is used to switch between the forces in exactly the same manner as it is used to switch constraints, without the constraint number argument. Since the forcing functions switch at the same times as the constraints, the same switch times can be used.

SolveMech is called with intermediate times in between the beginning and end times that are explicitly included in the points in time that are solved for. This ensures that a data point is included exactly at the peak of the force plot.