4. Velocity and Acceleration

Overview

This chapter covers MechanicalSystems' velocity and acceleration solution blocks, and issues related to the representation of angular orientation and velocity in spatial systems. Mech automatically differentiates the kinematic equations with respect to time and can solve for the velocity and acceleration of each body in a model without any further input from the user. Output functions are provided for utilizing the higher-order solution data, and definitions for the derivatives of Mech dependent variables are automatically made in the Global` context so that expressions containing Mech dependent variables can be properly differentiated.

4.1 Velocity and Acceleration Solutions

MechanicalSystems calculates all velocities and accelerations with respect to the time variable (the symbol T by default). As was shown in Chapter 2, the use of the time variable in a model is optional, if the model is to be used to generate a location solution only. However, a model must have some explicit functional dependence on the time variable in order to have a nontrivial solution to the higher-order (velocity and acceleration) equations.

4.1.1 2D Example Mechanism

To demonstrate the use of the Mech velocity and acceleration solutions, a 2D model is developed that can be measured. A 3D model is treated in the same manner as 2D with regard to higher-order solutions. The model is be a pull-rod and rocker rear suspension linkage for a car, consisting of three moving bodies: the chassis, the carrier, and the rocker. A real suspension linkage would have three more moving bodies: the upper and lower A-arms, and the pull-rod between the wheel carrier and the rocker.
The input to the model is the chassis moving vertically relative to the ground. This model could have been reduced in size by making the ground body the chassis, and simply moving the wheel carrier relative to the chassis, but the chosen method allows more flexibility and makes for better graphics.
The wheel carrier is attached to the chassis by two A-arms. These A-arms are not modeled as independent bodies because their function can be represented by a pair of relative distance constraints between the chassis and the wheel carrier.
The rocker pivots about an axis on the chassis and is attached to the wheel carrier by the pull-rod, a connecting link in tension. This pull-rod is also modeled with a relative distance constraint, not a separate body.
The shock absorber is connected between the other end of the rocker and the chassis. The damper is not actually part of the kinematic model at all, but the length of the damper as a function of wheel motion is a primary point of interest for any suspension linkage. Thus, the shock absorber exists only as a pair of mounting points on the chassis and the rocker, the separation of which can be measured as the car is bounced.

Here is a graphic of the 2D pull-rod suspension model.

Bodies

Four body objects are defined for the pull-rod suspension model.

Here are each of the bodies in the pull-rod model.

This loads the Modeler2D package and defines names for each of the bodies in the model.

In[1]:=

Here are the four body objects for the pull-rod model.

In[6]:=

In[7]:=

In[8]:=

In[9]:=

The body objects must be incorporated into the current model with SetBodies.

In[10]:=

Constraints

Seven constraints, one of which is a driving constraint with explicit dependence on time, are used to model the pull-rod suspension.

Here are the seven constraint objects for the pull-rod model.

In[11]:=

The constraint objects must be incorporated into the current model with SetConstraints.

In[19]:=

Note that the explicit symbol that is used to represent time in Modeler2D (T by default) may be changed. The SymbolBasis option for SetSymbols is a nested list of strings that determines the basis for all of the Global` symbols created by Modeler2D. Redefining the first element of SymbolBasis changes the symbol that is recognized as the time variable.

SymbolBasis is the basis for all Modeler2D global symbols.

In[20]:=

Out[20]=

Runtime

Before the model can be run, the symbol rodlength, which determines the length of the pull-rod, must be defined. Because of the presence of T in the driving constraint (constraint 1), the model can be run through its intended range of motion by varying T directly with the first argument to the SolveMech command. The way that this model is defined, the numerical value of T specifies the distance between the chassis and the ground.

The symbol rodlength must be defined.

In[21]:=

Now the model is run at T = 2.5.

In[22]:=

Out[22]=

Here is the pull-rod suspension model at T = 2.5.

4.1.2 The Solution Option

When a Mech model is constructed with a functional dependence on the time variable T its time-domain motion is fully defined. Thus, nothing further needs to be done to generate velocity and acceleration solutions other than to ask for them. This is done with the Solution option for SolveMech.

An option for SolveMech.

The Modeler2D pull-rod suspension model defined in Section 4.1.1 is used to demonstrate the velocity and acceleration solutions. The Solution -> Velocity option causes SolveMech to calculate the location and velocity of each body in the model.

This finds the location and velocity of each body at T = 3.0.

In[23]:=

Out[23]=

Mech creates a new set of symbols to represent the velocities and accelerations by appending a lowercase d or dd to each of the symbols that were used to represent body locations. Thus, while {X2, Y2} represents the location of the origin of body 2, {X2d, Y2d} represents the velocity of the origin of body 2.
Note that in this case, {X2d, Y2d} is equal to {0, 1}. This is expected since the driving constraint forces the height of the chassis to be T; therefore, its vertical velocity is equal to one.
The Solution -> Acceleration option causes SolveMech to calculate the location, velocity, and acceleration of each body in the model.

This finds the location, velocity, and acceleration of each body at T = 3.0.

In[24]:=

Out[24]=

The location, velocity, and acceleration of the rocker (body 4) reflect the fact that its origin is stationary on the chassis, but its rotation is controlled in a nonlinear manner by the suspension linkage.

Here is the rocker location, velocity, and acceleration.

In[25]:=

Out[25]//MatrixForm=

When SetConstraints is called to build the mathematical model, all of the expressions required for the velocity and acceleration solutions are not immediately built. Some parts of the mathematical model are built only when they are needed for the first time by SolveMech. Thus, the first time that SolveMech is called requesting a velocity or acceleration solution, it may take much longer to return a result than in any subsequent calls.
SetConstraints accepts the BuildMech option to force the modeling equations to be built to any level immediately, instead of waiting until they are needed.

An option for SetConstraints.

BuildMech can be used to build a model's velocity and acceleration equations before all of the symbolic parameters in the model have been defined, so that the parameters will remain variable in the model. If undefined parameters are present in a model, SolveMech cannot be called to seek a solution, but BuildMech can still build the required equations.

Undocumented Graphics Generation

In[29]:=

In[30]:=

In[31]:=

In[32]:=

In[33]:=

In[34]:=

In[35]:=

In[36]:=

In[39]:=

In[40]:=

In[43]:=

In[44]:=

In[45]:=

In[46]:=

In[47]:=

In[48]:=

4.1.3 3D Example Mechanism

A Modeler3D example model is developed to demonstrate the different angular velocity and acceleration formulations used by Modeler3D. This model is of a MacPherson strut automotive front suspension.
The model consists of only two moving bodies: the chassis and the wheel carrier. The wheel carrier is attached to the chassis by a lower A-arm, a MacPherson strut, and a steering tie rod that controls the steer angle of the wheel carrier. The A-arm, strut, and tie rod are not modeled as independent bodies in the model. They are each modeled as constraints that entirely represent the function they perform in the suspension. The input to the model is the chassis being moved vertically relative to the ground.
A MacPherson strut is attached to the chassis of a car by a ball joint, but it is attached to the wheel carrier in a cantilever fashion, not by another ball joint such as would be used in a more conventional double A-arm suspension. Thus, the MacPherson strut not only provides vertical support for the wheel carrier, it also provides the axis about which the wheel carrier pivots for steering. Since the strut cannot pivot with respect to the wheel carrier, it is modeled simply by forcing a ray on the wheel carrier to pass through a point on the chassis (the upper mounting point of the strut).

This loads the Modeler3D package.

In[26]:=

Here is a 3D MacPherson strut front suspension model graphic.

Bodies

Three body objects are defined for the MacPherson strut front suspension model.

Names for each of the body numbers in the model are defined.

In[27]:=

Here are the three body objects used in the MacPherson strut model.

In[30]:=

In[31]:=

In[32]:=

The body objects are incorporated into the current model.

In[33]:=

Constraints

Eight constraints, including a driving constraint that is explicitly dependent on the time variable, are required to model the MacPherson strut front suspension.

Here are the eight constraint objects for the MacPherson strut model.

In[34]:=

The constraints are incorporated into the current model.

In[42]:=

Note that the explicit symbol used to represent time in Modeler3D (T by default) may be changed. The SymbolBasis option for SetSymbols is a nested list of strings that determines the basis for all of the symbols created by Modeler3D in the Global` context. Changing the first element of SymbolBasis changes the symbol that is recognized as the time variable.

SymbolBasis specifies the basis for Modeler3D automatic symbols.

In[43]:=

Out[43]=

Runtime

Because of the presence of T in the driving constraint (constraint 3) the model can be run through its intended range of motion by varying T directly with the first argument to the SolveMech command. The way that this model is defined, the numerical value of T specifies the distance between the chassis and the ground.

Now the model is run at T = 2.5.

In[44]:=

Out[44]=

Here is the MacPherson strut front suspension image at T = 2.5.

4.1.4 3D Euler and Angular Formulations

Velocity and acceleration solutions are calculated in a 3D model in the same manner as in a 2D model. However, there are three options in a 3D model as to how the higher-order data is represented.
Since the angular orientation of Modeler3D bodies is specified by four Euler generalized parameters, {Eon, Ein, Ejn, Ekn}, one way of representing the angular velocity and acceleration of a body in 3-space is with the first and second derivatives of the Euler parameters, {Eond, ...} and {Eondd, ...}. Another more familiar way of representing angular velocity and acceleration is with the angular velocity and acceleration vectors, {omegaX, omegaY, omegaZ} and {alphaX, alphaY, alphaZ}, in either global coordinates or local coordinates attached to each individual body. Modeler3D can formulate the higher-order solutions in any of these representations, each of which requires that different sets of equations be built by the SetConstraints function. Which formulation is used is specified with the Method and Coordinates options for SetSymbols.

A system set-up function.

The Method option for SetSymbols may also be used with Modeler2D to specify angular coordinates or a degenerate version of Euler coordinates, although the Euler coordinates are seldom used in 2D models. Specifying Method -> Euler in a Modeler2D model causes the angular orientation of each body to be represented by a pair of coordinates {Ejn, Ekn} = {cos(n), sin(n)} instead of a single coordinate, n. Velocity and acceleration are also represented by the coordinate pairs {Ejnd, Eknd} and {Ejndd, Ekndd}. The inclusion of this feature in Modeler2D is largely academic; specifying Method -> Euler causes four equations to be generated for each body, instead of only three, and usually slows down the solution phase. There are two advantages afforded by 2D degenerate Euler coordinates: no trigonometric functions are generated in the constraint expressions so all of the constraints are purely algebraic, and the second derivative of a 2D rotation matrix is very simple which leads to faster solutions in some dynamical systems.

2D/3D

An option for SetSymbols.

3D

An option for SetSymbols.

The MacPherson strut front model defined in Section 4.1.3 is used to demonstrate the 3D velocity and acceleration solutions. The Solution -> Velocity option causes SolveMech to calculate the location and velocity of each body in the model.

This finds the location and velocity of each body at T = 2.5.

In[45]:=

Out[45]=

Note that the angular velocity of each body is returned in the form of derivatives of Euler parameters. This is because the constraint equations were built with all default options for SetSymbols. If the model is rebuilt with the Method -> Angular option for SetSymbols, the first derivatives of the Euler parameters will be replaced by angular velocities.

This rebuilds the constraints with angular velocity formulation.

In[46]:=

Now find the velocity solution at T = 2.0.

In[48]:=

Out[48]=

The angular velocity of the chassis in global coordinates is represented by the symbols {x2, y2, z2}.
The Solution -> Acceleration option causes SolveMech to calculate the location, velocity, and acceleration of each body in the model. Note that it takes more time to execute the following command the first time than it will on all subsequent executions. This is because the acceleration equations must be built the first time through.

Now find the acceleration of the bodies at T = 2.0.

In[49]:=

Out[49]=

Since angular acceleration is the first derivative of angular velocity, the angular acceleration vector of the carrier, is represented by the derivatives of the angular velocity {x3d, y3d, z3d}.
The model is now rebuilt with the Coordinates -> Local option so that angular velocity is represented in local, instead of global, coordinates.

This rebuilds the constraints with a local angular velocity formulation.

In[50]:=

Now find the velocity solution at T = 2.0.

In[52]:=

Out[52]=

Note that the local angular velocity of body 3 is only slightly different from the global angular velocity. This is because the angular orientation of the wheel carrier is very close to being aligned with the global coordinate system.

Undocumented Graphics Generation

In[81]:=

In[82]:=

In[83]:=

In[84]:=

In[85]:=

4.2 Velocity and Acceleration Output

Now that the velocity and acceleration of each mechanism body has been calculated, useful output may be generated by developing expressions that are functions of Mech velocity and acceleration variables, and subjecting those expressions to the standard solution rules produced by SolveMech.

4.2.1 3D Example Mechanism

The 3D MacPherson strut suspension model is used to demonstrate various methods of producing velocity and acceleration output with Mech. The model is redefined as follows in abbreviated form.

The complete MacPherson strut suspension model follows.

In[1]:=

In[3]:=

In[4]:=

Here the model is run at T = 2.5 to prove that it still works.

In[5]:=

Out[5]=

4.2.2 Output Functions

While Mech provides a wide array of output functions that return various quantities related to body location or orientation, only a small subset of those output functions are available in versions that are related to body velocities or accelerations. The most basic of those are the two higher-order analogies of the Location function. For each of the following functions, the point argument can be a Mech point object or just the arguments to a Mech point object.

Velocity and acceleration output functions.

The velocity of the steering tie rod attachment point on the wheel carrier is very near zero in the global reference frame because the chassis is allowed to move up and down while the bottom of the tire (the carrier body) remains in contact with the ground. Thus the motion of the carrier with respect to the ground is just a small tipping about the bottom of the tire.

This calculates the velocity of a point.

In[6]:=

Out[6]=

A related velocity function.

The velocity of the steering tie rod attachment point with respect to the chassis has a much larger vertical component, while the longitudinal and transverse components remain the same.

The velocity of a point, relative to another body, is shown here.

In[7]:=

Out[7]=

Two functions provide the higher-order analogies to the Rotation function.

2D

2D angular velocity and acceleration functions.

3D

3D angular velocity and acceleration functions.

Note that all Mech output functions formulate expressions in terms of derivatives of Euler parameters, or global angular velocities and accelerations, or local angular velocities and accelerations. Which form is returned depends on the settings of the Method and Coordinates options for SetSymbols that were in effect the last time that the model was built.
For example, the current model was built with the default Method -> Euler setting in SetSymbols; therefore, all output functions will return expressions in terms of derivatives of Euler parameters that can be replaced by a solution rule.

Here is the angular velocity of the wheel carrier.

In[8]:=

Out[8]=

Out[9]=

The Method -> Angular option is now used with SetSymbols. This changes the formulation of the Modeler3D output functions, but the numbers remain the same.

The model is rebuilt with the Angular option setting.

In[10]:=

The angular velocity of the wheel carrier is the same.

In[12]:=

Out[12]=

Out[13]=

Again, the internal representation of the angular velocity can be changed to local coordinates, but the numerical results from the output functions remain the same.

The model is rebuilt again with local angular velocity representation.

In[14]:=

The angular velocity of the wheel carrier is still the same.

In[16]:=

Out[16]=

Out[17]=

The only other output functions that are defined by Mech are the higher-order analogies to the Distance function. These two functions are defined largely for academic reasons, because they can be reproduced simply by differentiating the Distance function, as is shown in Section 4.2.3.

Two more higher-order output functions.

There is a slight difference between using the velocity output functions Velocity, DDistanceDT, and so on, and simply differentiating their first-order counterparts Location and Distance. The first- and second-order output functions use the Mathematica D operator to differentiate the local coordinates of each point, thereby assuming that the local coordinates are constant unless they are explicitly functions of time T.

Here are the length, and the first and second derivatives of the length of the strut.

In[18]:=

Out[18]=

4.2.3 Differentiating Expressions

Although Mech provides only a small number of predefined output functions that relate to velocity and acceleration, any expression that is written in terms of Mech dependent variables can be differentiated with Dt because definitions have been automatically made for Dt[X, T], where X is any Mech dependent variable.

Now the model is rebuilt with the Method -> Euler setting.

In[19]:=

Here is the rate of change of the camber angle.

In[21]:=

Out[21]=

Out[22]=

This generates a numerical result at T = 2.5.

In[23]:=

Out[23]=

Note that Dt automatically replaced derivatives of dependent variables with the appropriate symbols. This form of differentiation is valid even when crossing from Euler to angular coordinates.
Consider the following example. With the Method -> Euler setting for SetSymbols, differentiating the seven Modeler3D location coordinates associated with body 3 with respect to T has the apparent effect of appending a d to the name of each variable. These new symbols are explicitly included in the solution rules returned by SolveMech.
However, with the Method -> Angular option set, differentiating the seven Modeler3D symbols with respect to T produces different symbolic results. This shows the functional relationship between the Euler parameters and the angular velocity. After replacement with a SolveMech solution rule, the results are numerically identical.

This shows the derivatives of dependent variables.

In[24]:=

Out[24]=

Out[25]=

Out[26]=

Now the model is rebuilt with the Method -> Angular option.

In[27]:=

The derivatives of dependent variables have changed form, but their values remain the same.

In[29]:=

Out[29]=

Out[30]=

Out[31]=

The same holds true for the local coordinate formulation.

In[32]:=

Out[34]=

Out[35]=

Out[36]=

Note that the only difference between the derivatives of the Euler parameters in terms of global or local coordinates is that a few of the signs are swapped.

4.2.4 Euler's Angles

As shown in Sections 4.2.2 and 4.2.3, Euler parameters and their derivatives are directly related to the angular velocity vector of a spatial body, and a direct conversion between angular velocity and the derivatives of Euler parameters is possible. However, there is no vector representation of spatial angular orientation so other representations must be used.
One method of representing the angular orientation of a spatial body is to define a sequence of three finite rotations about three fixed axes that take a body from one orientation to another. These three rotations are called Euler's angles. Angular orientations represented by Euler's angles are readily converted to and from Euler parameters with a combination of Modeler3D output functions. Mech does not provide a single function to make this conversion, because there are many possible combinations of successive rotations that can be used to represent the angular orientation of a body.
The following example shows the general procedure required to convert an angular orientation from a sequence of finite rotations to a set of Euler parameters. In this particular case, the sequence of rotations are the most commonly used Euler's angles.

Here are the three rotation angles, and vector constants for each axis.

In[37]:=

A single rotation matrix that combines the three successive rotations is created by multiplying the rotation matrices associated with each.

In[43]:=

Out[43]//MatrixForm=

The EulerParameters function is used to convert the rotation matrix to a set of Euler parameters.

In[44]:=

Out[44]=

The conversion of Euler parameters to a sequence of finite rotations is not such a straightforward process. The general procedure is to measure the angle between an appropriate pair of axes located in the global and local reference frames. The following example shows how to convert the Euler parameters eulpars back to three Euler's angles.

Since Mech output functions that return expressions in terms of the symbolic Euler parameters are used, a list of replacement rules is needed.

In[45]:=

Out[45]=

Here is the second Euler's angle theta. It is the included angle between the two Z axes in the global and local reference frames.

In[46]:=

Out[46]=

Out[47]=

Here is the first Euler's angle phi. It is the angle between the local Z axis and the global -Y axis, as viewed looking down onto the global X-Y plane.

In[48]:=

Out[48]=

Out[49]=

Here is the last Euler's angle psi. It is the angle between the local Y axis and the global Z axis, as viewed looking down onto the local X-Y plane.

In[50]:=

Out[50]=

Out[51]=

Any angular orientation can be represented by a single finite rotation about a fixed axis. This single rotation angle and axis are directly related to the Euler parameters so a direct conversion can be made with the Rotation function.

Here is the single rotation about a fixed axis that is equivalent to the three successive rotations phi, theta, and psi.

In[52]:=

Out[52]=

Some readers may be more familiar with the use of quaternions to represent the orientations of spatial bodies. Quaternions are four-dimensional complex quantities, with associated rules for multiplication, that allow them to be used much like 3D rotation matrices. The numerical values of the four elements of the quaternion qw + qx i + qy j + qz k are equivalent to the numerical values of the four Euler parameters {eo, ei, ej, ek} representing the same orientation. Quaternions and Euler parameters are notationally and conceptually quite different, but they are computationally equivalent.

4.2.5 2D Example Mechanism

The 2D pull-rod suspension model that was presented in Section 4.1 is used to demonstrate various methods of producing velocity and acceleration output with Mech. The model is redefined as follows in abbreviated form.

Here is the complete pull-rod suspension model.

In[53]:=

In[54]:=

In[55]:=

Here the model is run at T = 2.5 to prove that it still works.

In[57]:=

Out[58]=

4.2.6 Interpolation

MechanicalSystems provides a function that enhances the functionality of the Mathematica Interpolation function when used with a Mech model. This function is described in Section 3.2, but this subsection covers its capability of automatically utilizing higher-order data to obtain higher accuracy the interpolated results.

Interpolation utilities.

To demonstrate the use of TimeInterpolate, we must have an expression to interpolate and a list of rules to interpolate it over.

shocklength is the length of the damper in the pull-rod suspension model.

In[59]:=

Here is a list of three location solutions.

In[60]:=

The following plot shows a course representation of the shock length versus time. Only three data points were calculated by the previous call to SolveMech, T = 3, 1.5, and 0, so the ListPlot graphic consists of only two straight line segments.

Here is a plot of the shock length versus time.

In[61]:=

Out[61]=

If we use TimeInterpolate to interpolate among the three points, the curve looks much smoother, but it is still a very poor representation, because there is no way for the InterpolatingFunction to know that the curvature increases as T approaches zero.

Here we interpolate the shock length function with a default InterpolationOrder that is too high.

In[62]:=

Out[62]=

Here is a plot of the zeroth-order representation of shocklength.

In[63]:=

Out[63]=

A new table of three solution rules is calculated, but this time velocity and acceleration data is included. To avoid recalculating all of the location solution points, SolveMech is passed the list of location solutions that have already been generated.

This generates a new list of solution rules using the existing rules as initial guesses.

In[64]:=

The higher-order slope and curvature information contained in acctab is automatically included in the new InterpolatingFunction object returned by TimeInterpolate.

Again, we interpolate the shocklength function.

In[65]:=

Here is a plot of the second-order representation of shocklength.

In[66]:=

Out[66]=

Here are the three plots shown together.

In[67]:=

Out[67]=

Interpolating Solutions

Another method of applying interpolation is to directly interpolate the rules returned by SolveMech. SolveMech accepts the symbol Interpolation as an option name to tell SolveMech to interpolate the rules it returns. Note that SolveMech does not interpolate variables that are constant.

An option for SolveMech.

Because the Solution -> Acceleration option is given in the following example, the InterpolatingFunction objects returned contain slope and curvature information from the velocity and acceleration solutions.

Here is a set of interpolated solution rules.

In[68]:=

Out[69]=

Here is another plot of shocklength.

In[70]:=

Out[70]=