Generalized Coordinate Modeling in the Mechanical Systems Pack

Overview

The Generalized Coordinate Package for the Mechanical Systems Pack is an extension package that allows some or all of the bodies in a Mech model to be represented by generalized coordinates. The primary utility of this technique is to increase the performance of multi-body models by reducing the total number of symbolic coordinates used.

Standard Mech models are made up of one or more independent bodies, each of which are represented by three, four, or seven Cartesian and Euler coordinates. While this representation is very general, and allows bodies to be interconnected in an unrestricted manner, it inevitably causes many more symbolic coordinates to be included in the model than are required to represent its motion.

For example, a 3D body which is attached to the ground with a revolute joint (a Revolute5 constraint in Mech3D) would normally be represented by seven coordinates. While all these coordinates are necessary to completely define the location and orientation of the body, we are really concerned with only one of them; the single coordinate representing rotation about the revolute axis. The other coordinates may be just excess baggage, conveying trivial, constant information that we do not need. The Mech Generalized Coordinate Package allows such a body to be represented by a single coordinate, supplied by the user.

Note that there are pros and cons to using generalized coordinates. Some information, such as constraint reaction forces, are lost when a constraint is eliminated from the model with generalized coordinates. This notebook contains general guidelines to help determine when to use generalized coordinates, and when to use Mech's standard representation.

Reference

Context

Mech`GenCrd2D` or Mech`GenCrd3D`

Package Version

1.2

Mathematica Version

2.2

Copyright

Copyright 1995, Dynamic Modeling

Keywords

dynamics, generalized coordinates, mechanical systems, mechanism

Requirement

The Mathematica Mechanical Systems Pack must be installed on the computer system on which the Generalized Coordinate Package is to be run.

Installation and Loading

This section covers the installation of the Mech Generalized Coordinate Package on various computer systems and loading the packages into a Mathematica session.

File Locations

The Mech Generalized Coordinate Package is contained in three standard Mathematica package files. These files must be placed in the same directory as the rest of the Mechanical Systems Pack files. While the names of the files will be truncated to eight characters on a Windows PC, all other platforms will have file names and locations similar to the following.

/usr/local/math/Applications/Mech/GenCrd2D.m
                                  GenCrd3D.m
                                  GenCrdCommon.m

A typical Unix installation.

The required addition to the $Path variable that enables Mathematica to find the Generalized Coordinate packages will already have been made by the Mechanical Systems Pack installation process. No other installation steps are required.

Loading the Packages

The Generalized Coordinate packages are loaded into Mathematica in the same manner as a standard Mathematica package; with a Get or Needs statement. The Mech packages may be loaded before loading GenCrd2D or GenCrd3D, otherwise they will be autoloaded.

Needs["Mech`GenCrd3D`"] loads all necessary components of the Mech packages that are not already loaded, and then loads the Generalized Coordinate packages.

Loading the Generalized Coordinate packages.

This loads the GenCrd2D package.

Needs["Mech`GenCrd2D`"]

Like the standard Mech packages, the Generalized Coordinate packages are not intended to be loaded twice. Thus, it is recommended that Needs be used to load the packages because Needs will not read a package if it has already been read.

The KillMech command will remove all Mech packages from a Mathematica session, including the Generalized Coordinate packages.

This removes the Mech2D or Mech3D standard and auxiliary packages.

KillMech[]

After executing KillMech nothing remains of the Mech session except whatever definitions the user had made in the Global` context, and a few undefined symbols that Mech creates in the Global` context. The Mech packages can be reloaded with Get or Needs after executing KillMech.

Generalized Coordinate Constraints

The generalized coordinate constraints are not, in fact, constraints at all. Each generalized coordinate constraint is named like a standard constraint, takes approximately the same arguments as a standard constraint, returns a SysCon constraint object that is passed to SetConstraints, and generally appears in a mechanism model in exactly the same way as the standard Mech constraints. The difference between them is best illustrated by an example.

Consider a standard Mech2D model containing a Translate2 constraint that models a translational joint between body 2 and the ground. The presence of the constraint object causes SetConstraints to do two things: three degrees of freedom (the symbols X2, Y2, and Th2) are added to the model to represent body 2, and two degrees of freedom are removed from the model by the two constraint equations in the constraint. The net effect, then, is to add one degree of freedom to the model that must be constrained by some other constraint. If the Translate2 constraint is replaced with a GenTranslate2 constraint, SetConstraints does only one thing: one degree of freedom is added to the model by adding a single symbol to the dependent variable list. The name of the symbol is specified in the GenTranslate2 constraint object.

Body 2, which was represented by three symbols, is now represented by only one. All of Mech's other constraint, load, output, and graphics functions can recognize and access body 2 just as if it were a standard body; the only difference is that the symbolic representation of all geometry on the body is a function of just one symbol, not three.

Note that generalized coordinate constraints cannot be used within the multistage constraint objects, TimeSwitch or StageSwitch. This is simply because including them in a multistage constraint would have no effect on the model; the generalized coordinate constraints contain no constraint equations, so there is nothing to switch on and off.

2D Constraints

There are only two 2D generalized coordinate constraints. The wide variety of standard 2D constraints cannot be meaningfully converted to generalized constraints; only constraints which attach exactly two bodies together such that only one degree of freedom remains between them are generally useful.

The arguments given to the generalized coordinate constraints are of the same form as the arguments given to the standard constraints. However, each of these constraints requires an additional argument to specify the name of, and an initial guess for, the new symbol that represents the added degree of freedom.

GenRevolute2[cnum, axis1, axis2, {sym, guess}] returns a constraint object that introduces a slave body into the model represented by a single generalized coordinate sym. The axis1 body is a slave to the axis2 body with one degree of rotational freedom.  The origins of the two axes are constrained to be coincident, and the counterclockwise angle from axis2 to axis1 is represented by sym.

GenTranslate2[cnum, axis1, axis2, {sym, guess}] returns a constraint object that introduces a slave body into the model represented by a single generalized coordinate sym. The axis1 body is a slave to the axis2 body with one degree of translational freedom. The two axes are constrained to be colinear, and the distance from the origin of axis2 to that of axis1, in the direction of axis2, is represented by sym.

The 2D generalized coordinate constraints.

Note that there is a clear distinction between a master and slave bodies. The slave body is the one that will be represented by the single new generalized coordinate. Thus, the slave body cannot be the ground body. The master body may be the ground body, another body represented by standard coordinates, or another body represented by generalized coordinates that was created with another generalized constraint. The master-slave relationship cannot be circular; the master-slave chain must originate from the ground, or from a standard body. Also, a slave body cannot be a slave to more than one master body.

Another restriction on the arguments to generalized coordinate constraints is that the axis objects, axis1 and axis2, must each be defined entirely on a single body. While standard constraints allow axis objects to span multiple bodies (such as a line originating on one body and ending on another) generalized coordinate constraints do not.

Example

The following example shows a standard constraint and a generalized coordinate constraint that could be used to perform exactly the same function in a model.

Here are two functionally equivalent constraints.

Needs["Mech`GenCrd2D`"]
ground = 1;
body2  = 2;
Translate2[2, Line[body2,  {0, 0}, {0, 1}],
              Line[ground, {0, 1}, {0, 2}]]
GenTranslate2[2, Line[body2,  {0, 0}, {0, 1}],
                 Line[ground, {0, 1}, {0, 2}], {newsym, 0.0}]

This removes the Mech packages.

KillMech[]

3D Constraints

There are only four 3D generalized coordinate constraints. The wide variety of standard 3D constraints cannot be meaningfully converted to generalized constraints; only constraints which attach exactly two bodies together such that only one or two degrees of freedom remain between them are generally useful.

The arguments given to the generalized coordinate constraints are of the same form as the arguments given to the standard constraints. However, each of these constraints requires additional arguments to specify the names of, and an initial guesses for, the new symbol or symbols that represent the added degrees of freedom.

GenCylindrical4[cnum, axis1, axis2, {alpha, guess}, {beta, guess}] returns a constraint object that introduces a slave body represented by two generalized coordinates alpha and beta. The axis1 body is a slave to the axis2 body with one degree of rotational freedom and one degree of translational freedom. The two axes are constrained to be colinear and the distance from the origin of axis2 to that of axis1, in the direction of axis2, is represented by alpha. The angle from the projection of the reference direction of axis2 to that of axis1, in the right-hand sense, is represented by beta.

GenOrthoRevolute4[cnum, axis1, axis2, {alpha, guess}, {beta, guess}] returns a constraint object that introduces a slave body represented by two generalized coordinates alpha and beta. The axis1 body is a slave to the axis2 body with two degrees of rotational freedom.  The origins of the two axes are constrained to be coincident, and the axes are constrained to be orthogonal.  The symbol alpha represents the angle between axis1 and the reference direction of axis2, and beta represents the angle between axis2 and the reference direction of axis1.

GenRevolute5[cnum, axis1, axis2, {sym, guess}] returns a constraint object that introduces a slave body represented by a single generalized coordinate sym. The axis1 body is a slave to the axis2 body with one degree of rotational freedom. The two axes are constrained to be parallel and their origins are constrained to be coincident. The angle from the projection of the reference direction of axis2 to that of axis1, in the right-hand sense, is represented by sym.

GenTranslate5[cnum, axis1, axis2, {sym, guess}] returns a constraint object that introduces a slave body represented by a single generalized coordinate sym.  The axis1 body is a slave to the axis2 body with one degree of translational freedom.  The two axes are constrained to be colinear, and the distance from the origin of axis2 to that of axis1, in the direction of axis2, is represented by sym.  The projections of the reference directions of the two axes are constrained to be parallel.

The 3D generalized coordinate constraints.

Note that there is a clear distinction between a master and slave bodies. The slave body is the one that will be represented by the single new generalized coordinate. Thus, the slave body cannot be the ground body. The master body may be the ground body, another body represented by standard coordinates, or another body represented by generalized coordinates that was created with another generalized constraint. The master-slave relationship cannot be circular; the master-slave chain must eventually lead back to the ground, or to a standard body.

Another restriction on the arguments to generalized coordinate constraints is that the axis objects, axis1 and axis2, must each be defined entirely on a single body. While standard constraints allow axis objects to span multiple bodies (such as a line originating on one body and ending on another) generalized coordinate constraints do not.

Example

The following example shows a standard constraint and a generalized coordinate constraint that could be used to perform exactly the same function in a model.

Here are two functionally equivalent constraints.

Needs["Mech`GenCrd3D`"]
ground = 1;
body2  = 2;
Revolute5[2, Line[body2,  {0, 0, 0}, {0, 1, 0}],
             Line[ground, {0, 1, 1}, {0, 2, 1}]]
GenRevolute5[2, Line[body2,  {0, 0, 0}, {0, 1, 0}],
                Line[ground, {0, 1, 1}, {0, 2, 1}], {newsym, 0.0}]

User-Defined Generalized Constraint

If the selection of predefined generalized coordinate constraints is not sufficient, the GenConstraint function can be used to create constraints with any number of generalized coordinates representing any geometric relationship between two bodies that can be written algebraically.

GenConstraint[cnum, point1, point2, rotation, {sym, guess}] returns a constraint object which introduces a slave body represented by generalized coordinates into the model. A point on the slave body, point1, is constrained to be coincident with a point on the master, point2. The angular orientation of the slave with respect to the master is specified by rotation. {sym, guess}, or {{sym1, guess}, ...}, are generalized coordinate symbols and their associated initial guesses which are added to the model.

The user-defined generalized coordinate constraint.

In Mech2D, rotation is a counterclockwise angle. In Mech3D, rotation is specified in the {ang, axis} format, and it can be a single rotation or a sequence of rotations of the form used by RotationLock3.

The trick to creating a valid constraint with GenConstraint is to make the pointi and rotation arguments symbolic functions of the generalized coordinate symbols symi. For example, the following constraint object would introduce a slave body that is constrained so that its origin lies on a parabola in the Y-Z plane, and its rotation is about the Z axis only. The position of the body along the parabola is represented by tau, and the rotation angle of the body about the Z axis is represented by phi.

Here is a GenConstraint example.

GenConstraint[1, Point[body2, {0, 0, 0}],
                 Point[ground, {0, tau, tau^2}],
                 {phi, {0, 0, 1}},
                 {tau, 0.0}, {phi, 0.0}]

This removes the Mech packages.

KillMech[]

Four-Bar Example Mechanism

To demonstrate the use of generalized coordinate constraints in a Mech2D model, a 2D model of a simple four-bar mechanism is developed. The model consists of three moving bars and one stationary bar, the ground. This same model is also used in Section 8.3 and Chapter 9 of the Mechanical Systems Pack manual, so it is presented very briefly here.

This loads the Generalized Coordinate packages and defines some useful constants.

Needs["Mech`GenCrd2D`"]
Off[General::spell1]
X = {1, 0};
Y = {0, 1};

Here is the four-bar mechanism, shown with the locus of a moving point.

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.

The ground (body 1) uses two points.
P0 is the attachment point of the drive bar.
P1 is the attachment point of the driven bar.

The drive bar (body2) uses two points.
P0 is the rotational axis of the drive bar.
P1 is the attachment point of the center bar.

The driven bar (body3) uses two points.
P0 is the rotational axis of the driven bar.
P1 is the attachment point of the center bar.

The center bar (body 4) uses three points.
P0 is the attachment point of the drive bar.
P1 is the attachment point of the driven bar.
P2 is an offset point used for tracking a point locus.

Here are the body objects for the four-bar model.

bd[1] = Body[ground = 1, PointList->{{-2, 3}}];  

bd[2] = Body[drivebar = 2, PointList->{{1, 0}}];

bd[3] = Body[drivenbar = 3, PointList->{{3, 0}},
  InitialGuess->{{-2, 3}, 0}];

bd[4] = Body[centerbar = 4, PointList->{{0, 3}, {1, 2}},
  InitialGuess->{{ 1, 0}, 0}];

The body properties are incorporated into the model.

SetBodies[bd[ground],
          bd[drivebar],
          bd[drivenbar],
          bd[centerbar]]

Constraints

Five constraints, one of which is a driving constraint, are used to model the four-bar mechanism. These constraints are simply standard Mech constraints which will be replaced with generalized coordinate constraints shortly.

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.

Here are the constraint objects for the four-bar.

cs[1] = RotationLock1[1, drivebar, 2 N[Pi] T ];
cs[2] = Revolute2[2, Point[drivebar,  0], Point[ground, 0] ];
cs[3] = Revolute2[3, Point[centerbar, 0], Point[drivebar,  1] ];
cs[4] = Revolute2[4, Point[drivenbar, 1], Point[centerbar, 1] ];
cs[5] = Revolute2[5, Point[drivenbar, 0], Point[ground, 1] ];

The constraints are incorporated into the current model.

SetConstraints[ Array[cs, 5] ]

Run-Time

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.

The model is run at T = 0.3.

SolveMech[0.3]

Here is the four-bar model at T = 0.3.

Generalized Coordinate Constraints

Now one of the standard constraints is replaced with a generalized coordinate constraint. Constraint 2, the revolute joint at the drive bar origin, is replaced first, and body 2 becomes a slave to the ground body, represented by the single coordinate alpha.

Here is the new constraint.

cs[2] = GenRevolute2[2, Axis[drivebar, 0, X],
                        Axis[ground,   0, X], {alpha, 0.0} ];
SetConstraints[ Array[cs, 5] ]

The modified model is run again at T = 0.3.

SolveMech[0.3]

Note that the symbols X2, Y2, and Th2 are no longer included in the solution rules returned by SolveMech. The location of body 2 is represented entirely by alpha, which has the same numerical value as Th2 used to have. The numerical value of alpha would be different from Th2 if the direction vectors of the two axes in the GenRevolute2 constraint were not identical (X and X).

Two more of the Revolute2 constraints can also be replaced, but not the fourth one. This is because each generalized coordinate constraint creates a slave body, and a slave body can be a slave to only one master body. Thus, there is no way to use four generalized coordinate constraints in a model with only three bodies. Any three of the four revolute constraints can be replaced, but not all of them.

Here are two more new constraints.

cs[3] = GenRevolute2[3, Axis[centerbar, 0, X], 
                        Axis[drivebar,  1, X], {beta, 0.0} ];
cs[4] = GenRevolute2[4, Axis[drivenbar, 1, X],
                        Axis[centerbar, 1, X], {gamma, 0.0} ];
SetConstraints[ Array[cs, 5] ]

The modified model is run again at T = 0.3.

SolveMech[0.3]

Now all of the Cartesian coordinates have been replaced with the generalized coordinates alpha, beta, and gamma. A quick degree of freedom count shows that this all makes sense: the three generalized coordinate constraints each add one degree of freedom to the model, while the two standard constraints Revolute2 and RotationLock1 constrain a total of three degrees of freedom.

The physical meaning of the new coordinates is quite clear: alpha is the angle between the X axis of the ground and the local x axis of the drive bar, beta is the angle between the local x axes of the drive and center bars, and gamma is the angle between the local x axes of the center and driven bars.

Output and Graphics Functions

Now that we have a solution in terms of the generalized coordinates, all of the standard Mech output and graphics functions allow it to be used in the same manner as a solution in Cartesian coordinates. For example, the location of a point is still found with the Location function and the rotation angle of a body is found with Rotation. The symbolic expressions returned by these functions reflect the new coordinate system being used, but the numerical results are equivalent.

Here are the Cartesian coordinates of point 2 on the center bar.

Location[centerbar, 2]
%/.LastSolve[]

Here is the global rotation angle of the driven bar. Note that it is numerically equal to the Th3 coordinate from the first set of solution rules.

Rotation[drivenbar]
%/.LastSolve[]

The location of point 2 on the center bar is a function of alpha and beta, but not gamma. This is because of the master-slave relationship of the generalized coordinate bodies; the center bar is a slave to the drive bar, which is a slave to the ground. Thus, the location of any point on the driven bar is a function of all three generalized coordinates, and the location of any point on the drive bar is a function of only alpha.

Mech graphics functions also operate without regard to the use of generalized coordinates. The following graphics object will generate the same image with any of the three constraint formulations used previously.

Here is the graphics object representing the four-bar model.

graph := Graphics[{Thickness[.015],PointSize[.025],
(* Ground    *) {Point[{-2, 3}], 
  Line[{{-1, -1.5}, {0, 0}, {1, -1.5}, {-2, -1.5}, {-2, 3}}]},
(* Drive Bar *) {RGBColor[1, 0, 0],
  Tangent[{0,0}, .5, Point[drivebar, 1], 0],
  Wheel[Point[drivebar, 0], 0.5, Faces->False] },
(* Driven Bar*) {RGBColor[0, 0, 1], Edge[drivenbar, {0, 1}] },
(* Center Bar*) {Vertex[centerbar, {0, 1, 2}],
                 RGBColor[0, 1, 0],Facet[centerbar, {0, 1, 2}]}},
  Frame->True,PlotRange->{{-2.4, 2.4}, {-1.6, 4.2}}, 
  AspectRatio->Automatic];

This shows the graphic at T = 0.45.

Show[{graph/.SolveMech[0.45]}];

This removes all of the Mech packages.

KillMech[]

Space-Crank Example Mechanism

To demonstrate the use of generalized coordinate constraints in a Mech3D model, a 3D model of a spatial slider-crank mechanism is developed. The model consists of three moving bodies; a slider, crank, and connecting rod. A geometrically similar model is presented in Section 3.2 of the Mechanical Systems Pack manual. The model presented here differs in that it uses a separate body to represent the connecting rod, instead of simply replacing it with a relative-distance constraint.

This loads the Generalized Coordinate packages and defines some useful constants.

Needs["Mech`GenCrd3D`"]
Off[General::spell1]
X = {1, 0, 0};
Y = {0, 1, 0};
Z = {0, 0, 1};

Here is the complete space-crank graphic.

Bodies

Four body objects are defined for the space-crank model. Note the use of the zeroth point to avoid having to define {0, 0, 0}.

The ground (body 1) has one point defined in its body object.
P1 is a point to locate the rotational axis of the crank.

The crank (body 2) has one point defined in its body object. 
P1 is the attachment point of the connecting rod. The center of the crank is at the local origin.

The connecting rod (body 3) has one point defined in its body object.
P1 is the attachment point of the connecting rod to the slider, while the connecting rod is attached to the crank at P0.

The slider (body 4) has no points in its body object, only an initial guess, because only the local origin P0 is used.

Here are the body objects for the space-crank model.

bd[1] = Body[ground = 1,
  PointList->{{ 0, 10, 12}}];

bd[2] = Body[crank = 2,
  PointList->{{0, 8, 0}},
  InitialGuess->{{0, 10, 12}, {1, 0, 0, 0}} ];

bd[3] = Body[conrod = 3,
  PointList->{{30, 0, 0}},
  InitialGuess->{{0, 18, 12}, {1, 0, 0, 0}} ];

bd[4] = Body[slider = 4,
  InitialGuess->{{20, 0, 0}, {1, 0, 0, 0}} ];

The body properties are incorporated into the model.

SetBodies[bd[ground],
          bd[crank ],
          bd[conrod],
          bd[slider]]

Constraints

Five constraints, one of which is a driving constraint, are used to model the space-crank mechanism. These constraints are simply standard Mech constraints which will be replaced with generalized coordinate constraints shortly.

A ProjectedAngle1 constraint specifies the angular relationship between a line on the ground body and a line on the crank. This driving constraint is used to rotate the crank.

A Revolute5 constraint forces the axis of the crank to be coincident with its pivot axis on the ground.

A Spherical3 constraint between the crank and the connecting rod models a ball joint.

An OrthoRevolute4 constraint between the slider and the connecting rod models a universal joint. This is used instead of a Spherical3 constraint so as to constrain the spin of the connecting rod about its own longitudinal axis.

A Translate5 constraint between the slider and the ground allows translation of the slider.

Here are the constraint objects for the space-crank.

cs[1] = ProjectedAngle1[1, Vector[crank, Y], Vector[ground, Y],
                           Vector[ground, X], 2 N[Pi] T];

cs[2] = Revolute5[2, Axis[crank, 0, X], Axis[ground, 1, X]];

cs[3] = Spherical3[3, Point[crank, 1], Point[conrod, 0]];

cs[4] = OrthoRevolute4[4, Axis[conrod, 1, Y], Axis[slider, 0, Z]];

cs[5] = Translate5[5, Axis[slider, 0, X, Y], Axis[ground, 0, X, Y]];

The constraints are incorporated into the model.

SetConstraints[ Array[cs, 5] ]

Run-Time

The model can now be run through its range of motion by varying T with the SolveMech command. The numerical value of T specifies the angle of rotation of the crank in revolutions.

The model is run at T = 0.05.

SolveMech[0.05]

Here is the slider-crank model at T = 0.05.

Generalized Coordinate Constraints

Now one of the standard constraints is replaced with a generalized coordinate constraint. Constraint 2, the revolute joint that models the rotational axis of the crank, is replaced first, and body 2 becomes a slave to the ground body, represented by the single coordinate alpha. Note that the arguments to the GenRevolute5 constraint are exactly the same as the arguments to the Revolute5 constraint, with the addition of the new symbol alpha and its initial guess.

Here is the new constraint.

cs[2] = GenRevolute5[2, Axis[crank, 0, X], Axis[ground, 1, X],
          {alpha, 0.0}];
SetConstraints[ Array[cs, 5] ]

The modified model is run again at T = 0.05.

SolveMech[0.05]

Note that all of the symbols associated with body 2 are no longer included in the solution rules returned by SolveMech. The location of body 2 is represented entirely by alpha, which is the rotation angle of the crank in full revolutions. The definition (and hence the numerical value) of alpha is dependent on the relative reference directions of the two axis objects passed to GenRevolute5. Since the reference directions were not explicitly specified in the axis objects, default reference directions were assumed by Mech (the local y axes of each body, in this case).

Only two more of the standard Mech constraints in this model can be replaced with generalized coordinate constraints. This is because each generalized coordinate constraint creates a slave body, and a slave body can be a slave to only one master body. Thus, there is no way to use more than three generalized coordinate constraints in a model with only three bodies.

Here are two more new constraints.

cs[4] = GenOrthoRevolute4[4, Axis[conrod, 1, Y],
                             Axis[slider, 0, Z],
          {beta, -1.0}, {gamma, 1.0} ];
cs[5] = GenTranslate5[5, Axis[slider, 0, X, Y],
                         Axis[ground, 0, X, Y], {delta, 20.0} ];
SetConstraints[ Array[cs, 5] ]

The modified model is run again at T = 0.05.

SolveMech[0.05]

Now the 21 Cartesian and Euler coordinates have been replaced with the four generalized coordinates alpha, beta, gamma, and delta. A quick degree of freedom count shows that this all makes sense: the three generalized coordinate constraints add a total of four degrees of freedom to the model, while the two remaining standard constraints Spherical3 and ProjectedAngle1 constrain a total of four degrees of freedom.

The physical meaning of the new coordinates is somewhat confusing: alpha is simply the angle between the global Y axis and the local y axis of the crank, and delta is the X direction distance from the global origin to the origin of the slider. beta and gamma are two angles that represent the rotation of the connecting rod with respect to the slider: beta is the altitude angle, the angle between a horizontal plane and a vector down the axis of the connecting rod, and gamma is the azimuth angle, the rotation of the connecting rod axis about the global Z axis.

Finally, to reduce the mathematical model of the space-crank to absolute minimum size, constraint 1 (the ProjectedAngle1 constraint) must be reexamined. This constraint uses a rather complicated expression relating the directions of three vectors to specify the rotation of the crank, but it is known that the rotation of the crank can be directly specified by alpha. While Mech attempts to simplify constraint expressions to some extent, the algebraic expression generated by constraint 1 is more complicated than it could be.

Here is the expression that results from constraint 1.

Constraints[1]

Constraint 1 really just sets alpha equal to 2 Pi T, which can be done more directly with the all-purpose constraint function Constraint.

Here is a simplified functional replacement for constraint 1.

Constraint[1, alpha == 2 N[Pi] T]

Output and Graphics Functions

Now that we have a solution in terms of generalized coordinates, all of the standard Mech output and graphics functions allow it to be used in the same manner as a solution in terms of Cartesian and Euler coordinates. For example, the location of a point is still found with the Location function and the spatial orientation of a body is found with Rotation. The symbolic expressions returned by these functions reflect the new coordinate system being used, but the numerical results are equivalent.

Here are the Cartesian coordinates of a point at the center of the connecting rod.

Location[conrod, {15, 0, 0}]
%/.LastSolve[]

The location of the point on the connecting rod is a function of beta, gamma, and delta, but not alpha. This is because of the master-slave relationship of the generalized coordinate bodies; the connecting rod is a slave to the slider, which is a slave to the ground. Thus, the location of any point on the connecting rod is a function of the generalized coordinates of the slider and the connecting rod, while the location of any point on the crank is a function of only alpha.

Here is the global rotation angle of the crank.

crrot = Rotation[crank]
crrot/.LastSolve[]

The symbolic expression returned by Rotation is rather complicated, considering that the correct result is known to be {alpha, {1, 0, 0}}. This is an unfortunate side effect of the general transformation that Mech uses to convert a rotation matrix to the single rotation angle and axis of rotation that it represents. The application of a couple of trig identities will reduce this expression back to the simpler form, but only in such trivial cases. In a more general case, such as the connecting rod, which has undergone two consecutive rotations relative to the ground, the symbolic expression for the angle and axis of rotation cannot be reduced significantly. 

Mech graphics functions also operate without regard to the use of generalized coordinates. The following graphics object will generate the same image with any of the three constraint formulations used in the preceding example.

Here is the graphics object representing the space-crank model.

graph := Graphics3D[ {{Thickness[.008],
  Line[{{-1,-5,0}, {-1,18, 0}, {35,18, 0}, {35,-5,0},
        {-1,-5,0}, {-1,-5,22}, {-1,18,22}, {-1,18,0}}],
  Line[{{-1, 1.8,0}, {35, 1.8,0}}],
  Line[{{-1,-1.8,0}, {35,-1.8,0}}]},
  Bar[Line[crank, 0, {-1,0,0}], 10.0, 12],
  Box[ slider, {{-4,-1.5,-2}, {4,1.5,2}} ],
  Bar[ Line[conrod, 0, 1], 1.0, 8 ]} ];

This shows the graphic at T = 0.15.

Show[{graph/.SolveMech[0.15]}];

This removes the Mech packages.

KillMech[]

Loads and Reactions

This section deals with issues related to applying loads to mechanisms modeled with generalized coordinates, and finding the reaction forces associated with those coordinates.

Applied Loads

The Generalized Coordinate Package has almost no effect on the modeling of applied Loads. The standard Mech load functions Force, Moment, Gravity, and GyroMoment are used to apply loads to a generalized coordinate slave body in exactly the same way as they are used with standard bodies. The one load function that is affected by the use of generalized coordinates is Load, which is used to apply a load directly to a specified coordinate.

Since Load references the coordinates in a model directly by their symbolic name, it can be used to apply a load directly to a generalized coordinate. Conversely, Load cannot be used to apply a load directly to a Cartesian coordinate of a slave body because a slave body doesn't have any explicit Cartesian coordinates.

For example, suppose a model has two bodies, 2 and 3, and body 3 is a slave to body 2. The two bodies are joined by a revolute joint and the symbol alpha represents the rotation of the slave, relative to the master. The Force function may be used to apply a force in the global X direction to the slave , but the Load function may not.

This loads the Generalized Coordinate Package and defines names for bodies 2 and 3.

Needs["Mech`GenCrd2D`"]
master = 2;
slave  = 3;

Here is a load object that applies a 10 unit force to the slave body.

Force[slave, Axis[slave, 0, {1, 0}], 10]

This load object does not apply the 10 unit force to the slave body, because X3 is not a coordinate in the model.

Load[X3, 10]

However, the Load function can be used to apply a moment to the slave body because the generalized coordinate alpha is an angular coordinate. A load applied directly to alpha essentially applies a counterclockwise moment to the slave body, and a clockwise moment to the master.

Here is a load object that applies a 10 unit moment to the slave, and an equal and opposite moment to the master.

Moment[slave, master, 10]

Here is a simpler load object that applies the same 10 unit couple between the slave and master.

Load[alpha, 10]

In general, the Load function should be used instead of Force or Moment whenever possible because the resulting algebraic expressions that are added to the model can be much simpler.

Constraint Reactions

While the use of generalized coordinates has little effect on the application of loads to a model, finding constraint reaction forces is another matter. When a constraint is eliminated from a model by using generalized coordinates, information regarding the reaction forces at the constraint is lost. This is one basic disadvantage of a reduced coordinate system; by algebraically eliminating coordinates from a model we lose information about those coordinates.

Thus, three restrictions are placed on models using generalized coordinates: the Reaction function cannot be used to find the reaction forces applied by a generalized coordinate constraint to any body (because it is not really a constraint), the Reaction function cannot be used to find the reaction forces applied by any constraint to a generalized coordinate slave body (because the body doesn't have Cartesian coordinates), and the Loads function cannot be used to find the loads applied to a slave body (for the same reason).

Both the Reaction and Loads functions may still be used with generalized coordinates, but they must be used to access the coordinates directly.

Reaction[cnum, sym] returns the generalized reaction applied by constraint cnum to generalized coordinate sym. The physical meaning of the returned scalar load is dependent on the nature of the generalized coordinate sym.
   
Loads[sym] returns the total external generalized load applied to generalized coordinate sym.

Loads[sym, Type->Reaction] returns the total generalized reaction load applied by all constraints to coordinate sym.

Loads[sym, Type->Dynamic] returns the total generalized inertial load applied to coordinate sym.

Alternate usages of Reaction and Loads.

The physical meaning of the quantities returned by Reaction and Loads is quite clear, in the context of the simple generalized coordinate constraints provided by this package. If the generalized coordinate sym represents an angle in radians, then the loads on the coordinate are in consistent units of moment. If the coordinate sym represents a distance, then the loads on the coordinate are in consistent units of force.

However, if abstract user-defined generalized coordinate constraints are used, the generalized loads may have little physical meaning, just like the generalized loads associated with each constraint (Lmb1, Lmb2, ...) sometimes have little physical meaning. The one meaning that always applies to a generalized load is this: a generalized load times a small variation in its associated generalized coordinate is equal to a small variation in work, in consistent energy units.

Slider-Crank Example Mechanism

To demonstrate the use of the Reaction and Loads functions with generalized coordinates a planar model of a simple slider-crank is developed. This model uses only two bodies, the slider and crank, and it is modeled with two generalized coordinates, the rotation angle of the crank and the X displacement of the slider.

This loads the Generalized Coordinate packages and defines some useful constants.

Needs["Mech`GenCrd2D`"]
Off[General::spell1]
X = {1, 0};
Y = {0, 1};

Here is the slider-crank mechanism.

Bodies

Two body objects are defined for the slider-crank model. Inertia properties are defined for each moving body so that the dynamic reaction forces can be calculated. No local point definitions are made in the body objects. All necessary point coordinates in this model are given explicitly in the constraint functions.

Here are the body objects for the slider-crank model.

ground = 1;
bd[2] = Body[crank = 2,
  Mass->200.0,
  Inertia->450.0];
bd[3] = Body[slider = 3,
  Mass->25.0,
  Inertia->0.02];

The body properties are incorporated into the model.

SetBodies[bd[crank], bd[slider]]

Constraints

Four constraints are used to model the slider-crank mechanism. Two generalized coordinate constraints are used along with two standard constraints.

A GenRevolute2 constraint defines the symbol alpha as the angular coordinate of the crank.

A RotationLock1 constraint controls the rotation of the crank.

A GenTranslate2 constraint defines the symbol beta as the displacement coordinate of the slider.

A RelativeDistance1 constraint models the connecting rod between the crank and the slider.

Here are the constraint objects for the slider-crank.

cs[1] = GenRevolute2[1, Axis[crank,  0, X],
                        Axis[ground, 0, X], {alpha, 0.0}];
cs[2] = RotationLock1[2, crank, 2 N[Pi] T];
cs[3] = GenTranslate2[3, Axis[slider, 0, X],
                         Axis[ground, 0, X], {beta, 5.0} ];
cs[4] = RelativeDistance1[4, Point[crank, {2, 0}],
                             Point[slider, 0], 3.5];

The constraints are incorporated into the current model.

SetConstraints[ Array[cs, 4] ]

Loads

A force that is a linear function of slider position is applied to the slider to simulate a simple spring. The spring has a spring constant of 10, and it is placed so that the spring force is zero when the X displacement of the slider is three units.

Here is a spring force applied to the slider.

SetLoads[
  ld[1] = Load[beta, -10 (beta - 3)] ]

The model is run at T = 0.1 with the Static option to cause the Lagrange multipliers to be calculated.

SolveMech[0.1, Solution->Static]

Reactions

Because generalized coordinates were used in this model, Mech's ability to find constraint reaction forces is limited. Normally, Reaction could be used to find the reaction force vector that the revolute joint applies to the crank, but with this model only the moment applied to the crank by constraint 2 can be found. Note that Reaction cannot directly find the moment applied to body 2, it can only find the generalized force associated with the symbol alpha. Since alpha is an angular coordinate, its associated load is a moment.

Here is the moment applied to the crank by the driving constraint.

Reaction[2, alpha]/.LastSolve[]

Of course constraint 4, the RelativeDistance1 constraint, applies an equal and opposite moment to the crank.

Here is the moment applied to the crank by the connecting rod.

Reaction[4, alpha]
%/.LastSolve[]

Constraint 4 also applies a reaction force to the slider. Again, the restraining force and torque that would be applied to the slider by a normal Translate2 constraint cannot be found. Only the generalized force associated with beta can be found, which is the X direction force on the slider.

Here is the force applied to the slider by the connecting rod.

Reaction[4, beta]/.LastSolve[]

Here are the external loads applied to the coordinates alpha and beta.

Loads[All]

The Generalized Coordinate Package automatically converts the inertia terms into the proper coordinate system so that the dynamic reaction forces may be calculated. The Solution->Dynamic option causes SolveMech to include the inertial loads and calculate the reaction forces.

Here is a plot of the torque required to drive the crank through a full cycle.

ListPlot[{alpha, Reaction[2, alpha]} /.
  SolveMech[{0,1}, 41, Solution->Dynamic],
  PlotJoined->True];

Here is the graphics function that generates the image of the slider-crank.

graph:=Graphics[{
  Thickness[.01],
  Line[{{1,-2.2},{0,0},{-1,-2.2},{5.5,-2.2},{5.5,-.25},{1.5,-.25}}],
  {Thickness[.02],RGBColor[1,0,0],
   Disk[{0,0},.5],Edge[crank, {0, {2,0}}]},
  {RGBColor[0,1,0],Thickness[.015],
   Edge[{Point[crank, {2,0}], Point[slider,0]}]},
  {RGBColor[0,0,1],Box[slider,{{-.5,-.25},{.5,.25}}]},
  PointSize[.02],Point[Location[slider,0]],Point[{0,0}],
  Thickness[.002],LocalAxes[{crank,slider}, 1]},
  AspectRatio->Automatic,PlotRange->{{-2.3,6.3},{-2.3,2.3}}]

Here is the slider-crank at T = 0.35.

Show[graph/.SolveMech[.35]];

This removes the Mech packages.

KillMech[]

Dynamic Modeling

This section covers topics related to generalized coordinate modeling of dynamic systems. Dynamic motion synthesis is often where generalized coordinates have their greatest payoff. Using a coordinate system that is closer to being a minimal and independent set can drastically reduce the size of the mass matrix and centrifugal force vector, resulting in reduced run-times for the numerical integration routine.

The reduced mass matrix, in terms of generalized coordinates, is often useful in its own right for certain types of analysis. Examples of various mass matrix manipulations are presented in this section, including the conversion of an open-loop mass matrix to end-point coordinates, and the reduction of a closed loop mass matrix to minimal coordinates.

Two-Link Manipulator Example

It is often useful to find the mass matrix of a model with respect to an alternative coordinate system. The Generalized Coordinate Package automatically generates the model's mass matrix in terms of the specified coordinates, and converting the mass matrix to another consistent coordinate system is a straight-forward process.

An open-loop two-link manipulator model is developed to demonstrate the conversion of a mass matrix from one coordinate system to another. The model has two moving bodies, link 1 and link 2, which are represented by two generalized coordinates, alpha and beta. The alpha coordinate represents the rotation angle of link 1 relative to the ground and beta represents the rotation of link2 relative to link 1.

This loads the Generalized Coordinate packages and defines some useful constants.

Needs["Mech`GenCrd2D`"]
Off[General::spell1]
X = {1, 0};
Y = {0, 1};

Here is a graphic of the two-link manipulator model.

Bodies

Two body objects are used in the manipulator model. Each body object is used only to define inertia properties for the body, no local point definitions are made. Thus the ground (body 1) needs no body object because its inertia is immaterial. Note that the lengths of the links are left as parameters L1 and L2, causing the resulting mass matrix to be completely general.

Here are the body objects for the manipulator model.

ground = 1;
link1  = 2;
link2  = 3;
bd[2] = Body[link1, Mass->1, Inertia->3, Centroid->{L1/2, 0}];
bd[3] = Body[link2, Mass->1, Inertia->3, Centroid->{L2/2, 0}];

The inertia properties are incorporated into the model.

SetBodies[ bd[link1], bd[link2] ]

Constraints

Four constraint objects are used to model the manipulator. Two constraints are generalized coordinate constraints which introduce the two degrees of freedom, alpha and beta, and the other two constraints simply constrain these two degrees of freedom.

A GenRevolute2 constraint sets up the relationship between link1 and the ground.

A GenRevolute2 constraint sets up the relationship between link2 and link1.

Two Constraint objects constrain alpha and beta to be equal to the parameters q1 and q2.

Note that L1 appears explicitly in constraint 3 to locate the tip of link1, where it is attached to link2.

Here are the constraint objects for the manipulator.

cs[1]=GenRevolute2[1, Axis[link1,  0, X],
                      Axis[ground, 0, X], {alpha, 0}];
cs[2]=Constraint[2, alpha == q1];
cs[3]=GenRevolute2[3, Axis[link2, 0      , X],
                      Axis[link1, {L1, 0}, X], {beta, 0}];
cs[4]=Constraint[4, beta  == q2];

The constraints are incorporated into the current model.

SetConstraints[ Array[cs, 4] ]

Mass Matrix

The kinematics of this model are truly trivial. The entire constraint set is nothing more than a pair of expressions setting alpha and beta equal to q1 and q2. There is really no point in "running" the model at all, in the normal sense in which other Mech models are "run" because the configuration of the model can be specified directly by alpha and beta. The usefulness of the model lies in the other definitions that have been made by Mech for inertia properties and coordinate transformations.

Here is the constraint vector.

Constraints[All]

Here is the global location of the tip of link 2.

link2tip = Location[link2, {L2, 0}]

Here is the 2x2 system mass matrix.

m = Simplify[ MassMatrix[All] ]

End Point Mobility

Now we wish to find the ellipse of mobility of the end point of link 2. The ellipse of mobility of an equivalent two-link manipulator model was found in Section 10.3 of the Mechanical Systems Pack manual. In that example, a one unit force was applied to the end point of link 2 while varying the direction of the force through 360 degrees, and the resulting free acceleration vector of the endpoint was found with the SetFree and SolveFree functions.

This section presents a more elegant method of finding the ellipse of mobility directly from the mass matrix. The method presented in Section 10.3 is much more general than this one in that this method (as it is presented here) is restricted to finding the mobility of the end point at zero velocity, and is restricted to models that have the same total number of coordinates as the number of physical degrees of freedom of the end point.

The basic trick here is to convert the mass matrix from the generalized (joint angle) coordinates to endpoint coordinates. The mass matrix m, in the context of f=m.a, transforms accelerations of the generalized coordinates a, which are the relative angular accelerations of each of the joints, into forces on the generalized coordinates f, which are moments applied at the joints. We need to have a matrix P in P.F=A, which transforms a force vector F applied to the endpoint into the linear acceleration vector A of the end point. First, we need a transformation matrix J from generalized coordinates to endpoint coordinates such that A = J.a. Then P=J.Inverse[m].Transpose[J].

The matrix J is the Jacobian of the global coordinates of the endpoint.

Here is the transformation matrix J.

J = Outer[D, link2tip, {alpha, beta}]

Here is the inverse end point mass matrix P, a rather large symbolic matrix.

P = J.Inverse[m].Transpose[J];

We can now find the acceleration of the endpoint that results from any force vector at any particular configuration simply by multiplying P times the force vector.

Here is the acceleration vector resulting from an applied force {0, 2}.

config = {alpha->N[Pi]/6, beta->2 N[Pi]/3, L1->6, L2->6};
P . {0, 2} /. config

Now finding the ellipse of mobility is simply a matter of applying a unit force vector that revolves through a full circle and plotting the components of the resulting acceleration.

Here is the ellipse of mobility.

ParametricPlot[Evaluate[P . PolarToXY[{1, ang}] /. config],
  {ang, 0, 2 Pi},
  Epilog->{Line[{{0, 0}, P . X /. config}]} ];

Four-Bar Example

Another useful manipulation of the mass matrix is to find a minimal mass matrix for a model with a non-minimal coordinate system. While this can theoretically be done with or without generalized coordinates, it often becomes intractable when modeling in Cartesian coordinates because of the need to invert a large symbolic matrix.

The four-bar model that was developed previously in this notebook is used again here to demonstrate the formulation of a minimal mass matrix. In this case, the minimal mass matrix is a 1x1 matrix because the four-bar model has only one degree of freedom. Since the four-bar is modeled with three generalized coordinates alpha, beta, and gamma, the default mass matrix generated by Mech is a 3x3 matrix; the task at hand is to reduce it to 1x1.

Since the Generalized Coordinate Package was loaded into Mathematica in the last section, it does not need to be loaded again here. Note, however, that inertia properties are included in the body objects of this version of the four-bar.

This clears all data in the current Mech model.

ClearMech[]

Bodies and Constraints

The entire four-bar model is redefined without explanation here, because its kinematics are completely documented in a previous section.

Here are the body objects for the four-bar model.

SetBodies[
  bd[1] = Body[ground = 1,
    PointList -> {{-2, 3}}],

  bd[2] = Body[drivebar = 2,
    PointList -> {{1, 0}},
    Mass      -> 5,
    Inertia   -> 1],

  bd[3] = Body[drivenbar = 3,
    PointList    -> {{3, 0}},
    Mass         -> 2,
    Inertia      -> 4,
    InitialGuess -> {{-2, 3}, 0}],

  bd[4] = Body[centerbar = 4,
    PointList->{{0, 3}, {1, 2}},
    Mass         -> 3,
    Inertia      -> 12,
    InitialGuess->{{ 1, 0}, 0}] ]

Here are the generalized coordinate constraint objects.

SetConstraints[
  cs[1] = RotationLock1[1, drivebar, 2 N[Pi] T ],
  cs[2] = GenRevolute2[2, Axis[drivebar, 0, X],
                          Axis[ground,   0, X], {alpha, 0.0} ],
  cs[3] = GenRevolute2[3, Axis[centerbar, 0, X], 
                          Axis[drivebar,  1, X], {beta, 0.0} ],
  cs[4] = GenRevolute2[4, Axis[drivenbar, 1, X],
                          Axis[centerbar, 1, X], {gamma, 0.0} ],
  cs[5] = Revolute2[5, Point[drivenbar, 0], Point[ground, 1] ] ]

Mass Matrix

We intend to transform the 3x3 system mass matrix M into a 1x1 mass matrix m where f = m.alphadd and f is the generalized force in the direction of alpha, which is simply a counterclockwise moment applied to the drive bar. This transformation is of the form m = Transpose[J] . M . J.

The necessary transformation matrix J converts motions in all three generalized coordinates into motions of the single coordinate alpha. This matrix J is created by inverting the part of the model's Jacobian associated with beta and gamma and multiplying it by the part associated with alpha. The resulting 2x1 matrix converts alpha only into beta and gamma, so it must be prepended with a 1 to convert alpha into all three coordinates.

Only the part of the Jacobian associated with constraint 5 is included in the calculation because the driving constraint (constraint 1) must be dropped for the model to have more than zero degrees of freedom, and the other 3 generalized coordinate constraints contribute no constraint equations.

Here is the 3x1 transformation matrix J.

J = Prepend[
  -Inverse[Jacobian[5, {beta, gamma}]] . Jacobian[5, alpha], {1}];

Here is the 3x3 system mass matrix.

M = Simplify[ MassMatrix[All] ]

Here is the 1x1 mass matrix in terms of the coordinate alpha.

m = Transpose[J] . M . J;

Since m is a 1x1 matrix it may as well be a scalar, and the acceleration that results from applying a moment to the drive bar is simply alphadd = f/m. Remember that m is still a function of alpha, beta, and gamma, and not just any set of alpha, beta and gamma will do (unlike the two-link manipulator in which any values of the coordinates could be specified). SolveMech must be used to find a set of coordinates that represents a valid configuration of the four-bar model.

The 1x1 matrix m is converted to a scalar.

m = m[[1,1]];

Here is the configuration of the four-bar at 1/8 turn of the drive bar.

config = SolveMech[1/8]

Here is the acceleration that results from a 10 unit moment applied to the drive bar.

a = 10/m /. config

Dynamic Motion

The constrained dynamic motion profile of the four-bar can be found with the SetFree and SolveFree functions, as described in Chapter 10 of the Mechanical Systems Pack manual. In this example, the drive bar is given an initial angular velocity of 2 pi radians/second, a condition that is inherent in the definition of constraint 1. Constraint 1 is dropped from the model, and the free motion of the underconstrained system is found.

The model is run at T = 0 to make initial conditions available to SetFree.

SolveMech[0, Solution->Velocity];

Here is the FreeSystem object that contains the dynamical equations of motion, with constraint 1 dropped.

fsys = SetFree[1, Solution->Dynamic]

The FreeSystem object fsys is integrated in the time domain from T = 0.0 to 2.0 with SolveFree. All initial conditions for the integration are taken from the most recent solution returned by SolveMech.

Here is the time domain solution.

sol = SolveFree[fsys, 1.0]

Here is a plot of the angular velocity of the drive bar through approximately one full revolution of the drive bar.

Plot[Evaluate[D[alpha/.sol, T]], {T, 0., 1.}];

This removes all of the Mech packages.

KillMech[]