10. Underconstrained Systems

Overview

This chapter covers MechanicalSystems' SetFree and SolveFree functions, which are used to synthesize the dynamic motion of underconstrained mechanical systems. SolveFree uses an adaptive multistep integration algorithm to integrate the equations of motion built by SetFree, and returns the dynamic motion profile of the model in the form of standard Mathematica interpolating functions. SetFree and SolveFree are also used to find the equilibrium configuration of underconstrained models that are subject to conservative static loading, or the equilibrium velocity of underconstrained models subject to velocity-dependent loading.
The SetPath and SolvePath functions are used to find the motion trajectories of models with one or more nonholonomic constraints. All Mech models must first be defined as fully constrained systems with SetConstraints before they can be modified with SetFree or SetPath.

10.1 Equilibrium Configuration

This section describes the use of the SetFree and SolveFree functions to find the equilibrium configuration of an underconstrained mechanism subject to static forces. A planar model of a pendulum sliding on a massless rope is used to demonstrate the use of SetFree and SolveFree.

10.1.1 Example Mechanism

The pendulum-on-a-rope model in the following example consists of two planar bodies, a slider that moves without friction along an infinitely flexible, inextensible, massless rope, and a pendulum that hangs below the slider and is attached to the ground by a spring. Before an equilibrium analysis can be done, the model must have a full set of kinematic constraints so that it can be assembled and each body can be moved to a reasonable position to use as an initial guess for the equilibrium analysis. Some of these kinematic constraints are then dropped to make the model underconstrained.

Bodies

Three body objects are defined for the pendulum model. Note that the zeroth point, which is defined to be {0,0} by default, is used in each body.

• The ground (body 1) requires three point definitions.
P0 is the origin of the ground body and the attachment point of the left end of the rope.
P1 is the attachment point of the right end of the rope.
P2 is the attachment point of a spring that is used to pull the pendulum to one side.
• The slider (body 2) requires three point definitions.
P0 is the origin of the slider and the point where the pendulum is attached.
P1 is the point where the rope exits the left side of the slider.
P2 is the point where the rope exits the right side of the slider.
• The pendulum (body 3) requires two point definitions.
P0 is the origin of the pendulum and its center of mass.
P1 is the point where the pendulum is attached to the slider.
• Constraints

Seven constraints are used to model the pendulum-on-a-rope model.

• A RelativeDistance1 constraint is used to model the left side of the rope. Since only the total length of the rope is known, the distance from the left end of the slider to the attachment point of the rope on the ground is simply constrained to be the unknown s. The symbol s is then added to the model as a dependent variable with the Constraint function.
• Another RelativeDistance1 constraint models the right side of the rope. The length of the right side of the rope is simply the total length of the rope minus the length of the slider and the length of the left side, 5 - 1 - s.
• A Constraint function is used only to introduce the variable s, and its initial guess, into the model as a dependent variable. The Constraint has no constraint equation, so it constrains zero degrees of freedom, but it does add one degree of freedom to the model by introducing s.
• A Revolute2 constraint attaches the pendulum to the slider.
• A RelativeX1 constraint initially controls the X coordinate of the slider. This constraint is dropped for the equilibrium analysis.
• A RotationLock1 constraint initially controls the rotation of the slider. This constraint is dropped for the equilibrium analysis.
• A RotationLock1 constraint initially controls the rotation of the pendulum. This constraint is dropped for the equilibrium analysis.
• Runtime

Because there is no occurrence of time T in the model, the model can be run with the SolveMech command with no time argument. Since there are no other symbolic parameters, the model can only achieve one configuration, the one specified by the constraints.

Forces

A model must have some nonzero loads applied to it in order to have a unique equilibrium position. The pendulum-on-a-rope model has three forces applied to it: a spring force applied to the pendulum, and gravity applied to the pendulum and the slider.

10.1.2 Finding Equilibrium

To find the equilibrium configuration of a mechanism, a FreeSystem object is built that contains a new set of equations of motion with one or more constraints removed from the model.

Finding equilibrium configuration.

The pendulum-on-a-rope model has three constraints that are present only to complete the kinematic model, constraints 5, 6, and 7. Each of these constraints is dropped, one at a time, and the equilibrium position is found with each new degree of freedom. The first constraint to be dropped is constraint 7, which locks the rotation of the pendulum relative to the slider. This allows the pendulum to rotate, but the slider stays locked at the center of the rope.

Next, constraint 5 is dropped, allowing the slider to move along the rope. Finally, constraint 6 is dropped, allowing the slider to rotate. The InitialGuess option is used to pass SetFree the solution to the last FreeSystem, which is named sol. This solution is a better initial guess than the defaults pulled from LastSolve[].

An option for SetFree.

10.2 Equilibrium Velocity

This section covers the use of the SetFree and SolveFree functions to find the equilibrium velocity of a mechanism that has velocity-dependent loading. A mechanism must have some loading that is a direct function of velocity for it to have a defined equilibrium velocity. Note that the centrifugal forces in a model, although they are explicitly functions of velocity, are inertial forces and do not appear in a velocity equilibrium analysis.

10.2.1 Example Mechanism

This simple planar slider-crank model is presented with little explanation because all of the techniques used in the model have been explored in previous chapters.
The one feature of interest in this model is the velocity-dependent applied loading. Two loads are applied, a driving moment on the crank and a viscous friction on the slider. Note that the driving moment applied to the crank is a function of velocity such that the moment drops to zero as velocity increases. This is done so that the velocity will not approach infinity when the piston is not moving at top-dead-center and bottom-dead-center.

10.2.2 Finding Equilibrium Velocity

To find the equilibrium velocity of a mechanism, a FreeSystem object is built that contains a new set of equations of motion with one or more constraints removed from the model.

Velocity equilibrium.

The driving constraint is dropped from the slider-crank model to give the model one degree of freedom.

To find the equilibrium velocity at a sequence of positions through a full cycle of the mechanism, the InitialCondition option is used to pass SolveFree a new set of initial conditions at each position where the equilibrium velocity is to be found. Thus, it is not necessary to create a new FreeSystem object for each solution point.

An option for SetFree and SolveFree.

10.3 Free Acceleration

10.3.1 2D Example Mechanism

A model of a two-link planar robot arm is developed to demonstrate the use of SetFree and SolveFree to find the instantaneous free acceleration of a mechanism. The robot arm has a force of unit magnitude, but varying direction, applied to the tip of the second link, while the two revolute joints in the arm are free to rotate. By finding the acceleration vector of the tip that results from a 1-unit applied force at a particular configuration of the arm, the "ellipse of mobility" can be plotted.
Of the four constraints that make up the robot arm model, two are Revolute2 constraints that model the joints between the ground and link 1, and between link 1 and link 2. The other two RelativeAngle1 constraints are present only to set up the initial position of the arm. They are dropped by SetFree to provide two unconstrained degrees of freedom. The parameter loadangle is used to change the direction of the applied load at the tip of link 2.

10.3.2 Instantaneous Free Acceleration

To find the free acceleration of an underconstrained mechanism a FreeSystem object is built that contains a new set of equations of motion with one or more constraints dropped from the model.

Instantaneous free acceleration.

The two RelativeAngle1 constraints are now dropped from the robot arm to give the model two degrees of freedom. If initial conditions are not specified, SetFree takes them from the current default initial guesses LastSolve[].

To find the free acceleration corresponding to a sequence of different values of loadangle, SetParameters could be used between each successive run of SolveFree, or the InitialCondition option can be used to pass SolveFree the new parameter values. SolveFree reads all its initial conditions from the InitialCondition option, Parameters[], and the initial conditions contained in the FreeSystem object, in that order.

10.3.3 3D Example Mechanism

The following mechanism is a simple gyroscope that, dynamically, is not a very simple mechanism. The gyroscope consists of three moving bodies: the base, the gimbal, and the rotor. The base is constrained to rotate about a vertical axis affixed to the ground, while the gimbal rotates about a horizontal axis affixed to the base, and the rotor is constrained to rotate about a fixed axis on the gimbal that is orthogonal to the gimbal's horizontal rotational axis on the base.
The model is defined so that the origins of each body are always coincident, and so that each body is initially aligned with the global coordinate system, therefore no initial guesses are required because the default initial guesses are exact. The two Orthogonal1 constraints 2 and 4 are present only so that the model can be assembled in a fixed state. They will be dropped to find the free response. A single parameter freq is included in the model to specify the angular speed of the rotor, relative to the gimbal. The gyroscope is modeled in local, angular coordinates so that the angular velocities of interest appear directly in the solution rules.
All of the bodies are massless except the rotor. This makes the results much easier to interpret. The moments of inertia of the rotor are representative of a thin disk.
Only one load is applied to the model, a moment applied to the base in the Z direction. Initially this moment has no effect because the base is prevented from rotating by constraint 2. However, when the constraints that prevent the base and gimbal from rotating are dropped, the orthogonal reaction to the moment applied to the base can be observed at the rotor.

10.3.4 Instantaneous Free Acceleration

To find the free acceleration of the gyroscope, one or more constraints must be dropped. The constraint that prevents rotation of the base, constraint 2, is dropped first. The rotor cannot precess in this state because the gimbal is still locked to the base, but the reaction torques applied to the base that accelerate the rotor can be measured.
No initial conditions are explicitly given. Initial guesses are taken from the default initial guesses that were set by the last run of SolveMech.

To make this model a little more interesting, the base is given a nonzero angular velocity with the InitialCondition option for SolveFree. This induces gyroscopic torques in the base constraint. To be consistent with the constraints the Z components of the initial angular velocity of each body must be identical.

Now the FreeSystem is rebuilt with both the base and the gimbal unlocked. The same initial velocity for the base, gimbal, and rotor is used so that the gyroscope precesses.

10.3.5 Gyroscopic Moments

Gyroscopic forces and moments modeled as per the example in Section 10.3.4 work well for single point analyses. However, if the time-domain free response is desired, such a modeling technique has a serious computational problem. While the rotor is spinning rapidly, the Euler parameters that represent the angular orientation of the rotor are changing sinusoidally and very rapidly. To numerically integrate such a system, the integrator must faithfully track very many periods of the sinusoidal Euler parameters to cover only a few periods of the motion of the more interesting parts of the gyroscope. If the period of the rotor's local rotation is much, much shorter than the period of interest, the integration becomes computationally infeasible.
The way around this problem is to introduce a special forcing function that mimics the gyroscopic torque of a spinning axially symmetric body and apply it to the rotor. This allows the rotor to be rotationally static.

A function to simulate gyroscopic moments.

GyroMoment has two important caveats attached to its use: the gyroscopic motion that it models is that of a body rotating about a principal axis of inertia and the two principle moments not associated with the axis of rotation must be equal. For most practical purposes, this means that the body is axially symmetric, and it is rotating about its axis of symmetry.
To demonstrate the use of GyroMoment, the gyroscope model that was presented in Section 10.3.3 is modified to remove the large angular velocity of the rotor and replace it with a GyroMoment load.

The constraint reaction forces that result from using this method are identical to the reaction forces that result from the normal procedure, without the use of GyroMoment. The true values of alpha and omega for the simulated rotor are reproduced by the following calculation.

The equations of motion that are produced by this method are essentially equivalent to the equations of motion without the spinning rotor simulation, with one important difference, these equations of motion can be easily integrated to obtain the time-domain motion of the entire gyroscope assembly without an overwhelming computational burden.

10.4 Dynamic Motion Synthesis

10.4.1 Numerical Integration

Numerical integration with SolveFree.

There are several issues that are specifically associated with the numerical solution of the constrained equations of motion generated by MechanicalSystems. The equations of motion are a nonminimal differential-algebraic system of equations, that may or may not be linear in the highest-order terms.
One important issue is constraint violation. The equations of motion do not include the actual kinematic constraints, they only include the second time derivatives of the constraint equations. Because of this, small errors in the numerical integration of the equations of motion can accumulate into large violations of the kinematic constraints and their first derivatives. To deal with this problem, the equations can be algebraically reduced to a minimal set, but this is not generally possible with the type of constraints allowed by Mech. Alternatively, the numerical values of the dependent variables can be corrected as the integration progresses so as to continuously satisfy the constraints. This method is implemented in SolveFree.
Another issue is the transformation of angular velocity and acceleration in 3D models. If a model is built with angular coordinates (the Method -> Angular option in SetSymbols) the first and second derivatives of angular orientation are represented by the angular velocity, omega, and the angular acceleration, alpha. Unfortunately, omega is not the first derivative of the Euler parameters used to represent angular orientation. Thus, a coordinate transformation must occur at each time step to allow the Euler parameters to be integrated. SolveFree handles this transformation automatically, so models built in any coordinate system can be integrated.
The equations of motion are not necessarily linear in their highest-order terms, depending on the form of the applied loads. SolveFree automatically switches to an iterative solution method for nonlinear systems that allows them to be solved, albeit slowly.

Initialization Options

The only option that is concerned with the initialization of the integrator is the InitialCondition option. All other information needed to integrate the equations of motion is contained in the FreeSystem object. This option is most useful for repeatedly integrating the same FreeSystem with different sets of initial conditions without having to rerun SetFree.

An initialization option for SolveFree.

The list of rules required by InitialCondition is of the form returned by SolveMech. Using SolveMech to generate the initial conditions (which is essentially what is done by the default setting) guarantees that the initial conditions constitute a valid solution to the constraint equations. If initial conditions that do not satisfy the constraints are passed to SolveFree, the integration will proceed without difficulty, producing results that satisfy the differential equations of motion, but do not satisfy the kinematic constraints.
Note that a sufficient set of initial conditions for integrating a free acceleration system is returned by SolveMech[time, Solution -> Velocity], while a sufficient set of initial conditions for integrating an equilibrium velocity system is returned by SolveMech[time].

Integration Options

The following options are used to control the actual integration algorithm used by SolveFree. The default settings for these options are intended to produce reasonable results in the shortest possible time with a well-behaved system. Some tweaking may be required to get the best result from models with large accelerations or force discontinuities.

Options for controlling the SolveFree integrator.

More options for controlling the integrator.

The error estimation used by SolveFree makes an assumption that the approximate magnitude of each of the dependent variables is near to 1. While this assumption is certainly valid for Euler parameters that are between -1 and 1, and the 2D angular coordinate that has units of radians, it is up to the user to choose a system of units such that the linear dimensions of the model are not extremely large or small. While the error estimation routine does not look at the Lagrange multipliers (1, 2, ...) it is best to keep the units of force and mass such that the Lagrange multipliers are not extremely large or small, because other parts of the numerical solution method begin to produce significant errors when dealing with a wide range of numerical magnitudes.
When step size switching is enabled, SolveFree basically limits the maximum error in any dependent variable at a single time step to be less than the MaxError setting. If this is exceeded at a time step, the integrator backs up, halves the step size, and proceeds forward again. This error limit seems rather loose. If the same variable exhibited the maximum signed error at each time step, the cumulative error in that variable after 1000 time steps would be 0.1! In reality, the large errors occur at discrete points in the motion of a model, so the cumulative error is normally smaller.
The MaxError option can also be used to explicitly specify the error magnitude limits at which the time step size is halved, doubled, or when a corrective versus predictive step is taken. MathError -> {a, b, c} specifies that the time step size will be doubled if the sum of the maximum errors in the last four steps drops below c, it will be halved if the error in the current step goes above a, and an Adams-Moulton corrective step will be taken if the error is between a and b. Thus, MaxError -> {0.0001, 0.00001, 0.000001} is equivalent to the default setting of MaxError -> 0.0001.

Iterative Solution Options

The following options are used only to control the iterative solution block used for nonlinear models. These options have no effect on linear models.

Options for the Newton-Rhapson solution block in SolveFree.

Free acceleration models are always linear in the acceleration variables (X2dd, Y3dd, ...) but they may be nonlinear in the Lagrange multipliers (1, 2, ... ) if there are loads applied that are algebraically dependent on reaction forces, such as friction. Equilibrium velocity models can be nonlinear in the Lagrange multipliers, and they can also be nonlinear in velocity variables (X2d, Y3d, ...) if the applied loads are nonlinear functions of velocity. The second element in a FreeSystem object tells whether the model is linear or not. FreeSystem[...][[2]] returns True for nonlinear models.

Output Control Options

The following options are used to determine what data is returned by SolveFree, and what format is used to returned it.

Output options for SolveFree.

The StepSize data returned by the Diagnostics setting can be useful to help determine what should be chosen for a constant step size, if one is desired. The MaxError data included by Diagnostics is the measure of error used by the adaptive time step adjustment algorithm to decide when to change the step size.

10.4.2 Planar Example Mechanism

The following example uses the pendulum-on-a-rope model that was developed in Section 10.1. The pendulum is released from a constrained zero velocity position, that is not in static equilibrium, and the free response is found.
The spring force that was used to pull the pendulum off to one side is removed from the model, but gravity is retained. Also, nonzero moments of inertia are defined for each body, since the moments of inertia are relevant to the dynamics of the system, while they were not relevant to the equilibrium analysis.

The pendulum-on-a-rope model uses seven constraint objects. The constraints that fix the initial position of the model have been altered so that the pendulum may be released from an offset position.

10.4.3 Free Response

To find the free acceleration of the pendulum, a FreeSystem object is built that contains a new set of equations of motion with constraints 5, 6, and 7 dropped from the model. Note that no initial conditions are yet available for the velocity of the bodies in the model. SetFree sets all of the initial velocities to zero if none are specified with the InitialCondition option.

Since the pendulum model is a conservative system, the sum of the kinetic and potential energy in the system should remain constant. This fact can be used to measure the accuracy of the numerical integrator.

A special output function.

The kinetic energy of each body is given by the BodyEnergy function, and the potential energy is simply the gravitational potential energy mass G altitude. The following plot shows that the numerical error in the integration has allowed the total energy of the model to increase by about 0.04% over 4.5 seconds.

Tightening the error budget by a factor of 100 by specifying the MaxError -> 0.000001 option for SolveFree will reduce the total energy rise to only 0.0007%, over the same time period, with an associated doubling of run time. Using the ConstraintCorrection -> True option, with the default error budget, will reduce the total energy rise to 0.002% with a slightly smaller increase in run time. In this case, either constraint correction or reducing step size (which is the result of decreasing the error budget) gives similar returns on computer time investment.
Changing the FitDegree setting can also affect the integration error. Setting FitDegree to Quadratic for this example will increase the net energy rise to 0.08%, while setting FitDegree to Quartic will reverse the sign of the error, changing the net energy rise to - 0.04%. Specifying a higher FitDegree setting often produces unexpected results because the automatic time step size adjustment works to counteract the higher degree setting. The integrator measures a smaller error in each step so it chooses a larger step size, resulting in the same net error.

10.4.4 Gyroscopic Example

The following example integrates the motion of the gyroscope model that was presented in Section 10.3. The model uses the GyroMoment function, which was also presented in Section 10.3, to reduce the numerical burden of integrating the motion of the rotor. The angular local coordinate system is used, which requires SolveFree to use a coordinate transformation during the integration.

10.4.5 Free Response

To find the unconstrained motion of the gyroscope, a FreeSystem object is built that contains a new set of equations of motion with constraints 2, 4, and 6 dropped from the model. Constraint 6, the constraint that locks the rotation of the rotor to that of the gimbal, is dropped in this simulation because the gross motion of the rotor results in some relative angular acceleration between these two bodies. If constraint 6 was not dropped, the model would represent an enforced constant angular velocity between the rotor and the gimbal, as if the rotor had a constant velocity motor attaching it to the gimbal.
The constant moment that was applied to the base is retained, so the model is expected to gain energy with the passage of time.

Note that the Rotation function was not used in the preceding example to resolve the rotation angle of the base from the Euler parameters, but ProjectedAngle was used instead. This was done because the small numerical errors in the integration caused one of the Euler parameters to drift slightly above 1.0. The Euler parameters are not meaningful for absolute values greater than 1, so the Rotation function would return an error. The ProjectedAngle function uses vector algebra, instead of a mathematical identity, and it does not care if the values of the Euler parameters stray above 1.
Since the gyroscope model is a conservative system (if the applied moment is recognized as a conservative force) the sum of the kinetic and potential energy in the model should be constant. The potential energy stored by the constant applied moment is the work done by the gyroscope on the applied moment, which is the negative of the magnitude of the applied moment times the is rotation of the base. The following plot shows that the numerical error in the integration has allowed the total energy of the model to decrease by about 2% of the maximum kinetic energy over the integration period.

Tightening the error budget by a factor of 10 with MaxError -> 0.00001 will reduce the total energy rise to only 0.12% over the same period, with an associated increase in run time. Using the ConstraintCorrection -> True option, with the default error budget, will have virtually no effect on the error in the total energy. This is because the error stems from the deviation in the Euler parameters that are very close to 1, an error on which constraint correction has little effect. Using the FitDegree -> Quartic option, with the default error budget, will reduce the total energy rise to 0.71%, with only a small increase in run time.

10.4.6 Reducing Gyroscopic Models

There is one more technique that can be applied to 3D gyroscopic models to reduce their size and thereby reduce their run time. This technique can be very useful to models with more than one rotor.

A function for modifying an inertia tensor.

The restrictions that apply to the use of GyroFilter are identical to the restrictions that apply to the use of GyroMoment. The body represented by the inertia properties must be axially symmetric. The purpose of GyroFilter is to make it possible to add a zero-torque spinning body to a model without actually adding a body to the model at all. To demonstrate, the rotor body is deleted from the previously defined gyroscope model, and replaced with a special load and inertia properties on the gimbal. This time, the model is built in Euler coordinates, just to show that the results are the same.

The model is now functionally equivalent to the previously defined model. The model can be integrated with SolveFree, after dropping constraints 2 and 4. Constraint 6 was previously dropped and is no longer part of the model.

10.5 Nonholonomic Systems

Nonholonomic systems, or systems that have a kinematic constraint that is an explicit function of velocity, are dealt with using the SetPath and SolvePath functions. SetPath is used to take a normally constrained kinematic model and replace one of its constraints with a constraint that is a function of mechanism velocities. SolvePath then integrates the resulting first-order system to obtain a time history of the motion of the mechanism.

10.5.1 Example Mechanism

The following planar example is probably one of the more common nonholonomic systems, a caster, such as a wheel on a shopping cart. The caster is modeled in two dimensions by viewing the caster from above, like looking down through the shopping cart. The steady state motion of a caster is for the wheel to trail directly behind the pivot point of the caster. What is modeled here is the motion of the caster when the direction of travel of the pivot is reversed. The caster swings out to one side of the pivot and then settles in behind the pivot again.
The wheel of the caster is modeled as a single nonholonomic constraint that constrains the direction of travel of the wheel to be parallel to the direction that the wheel is pointed. In vector algebra terms the velocity vector of the wheel is perpendicular to a vector through the wheel's axis of rotation. Normal Mech constraint functions can be used to model first-order relationships such as this by making at least one of the geometry specifications in the constraint a function of velocity.
Before the velocity constraint can be put in place by SetPath the model must be built in a normal kinematically constrained way. This requires a dummy constraint that will be replaced by SetPath later. The dummy constraint is a RotationLock1 constraint to lock the rotation of the caster relative to the ground. The other constraint in the model, constraint 1, simply constrains the pivot point of the caster to move steadily along the global X axis.

10.5.2 Nonholonomic Constraints

The nonholonomic system solving functions.