Dynamic simulation of mechanisms

Abstract

A method of dynamic simulation of mechanisms which constructs the lagrangian form of the dynamics equations as pure differential equations in terms of a minimal set of generalized coordinates. The dynamics problem is separated into two parts used in the lagrangian--a kinematic analysis and a kinetics analysis. The kinematic analysis is computed using a degrees of freedom analysis. The velocities required for the kinetic analysis are computed using screw theory. Once constructed, the lagrangian dynamics equations are solved for accelerations, applied to the mechanism, and integrated over time to simulate the dynamics of the mechanism. The kinetic and kinematic analyses are approximately linearly related to the number of components in the mechanism, allowing efficient solution of dynamics problems at interactive rates.

Description

TECHNICAL FIELD

The present invention relates to a system and method for simulating the dynamics of a mechanism, and in particular, for constructing the lagrangian form of the dynamics equation with pure differential equations permitting interactive simulation

BACKGROUND

Understanding complex mechanical devices requires the ability to simulate the behaviors of the devices, and to make reasonable generalizations about those behaviors. Mechanism simulation may be done at a purely qualitative level, a mixed numerical/qualitative level, or at a detailed numerical level followed by qualitative characterization of the results. See, Hyun-Kyung Kim. Qualitative kinematics of linkages. Technical Report UIUCDS-R-90-1603, University of Illinois, Champaign-Urbana, May 1990; L. Joskowicz and E. P. Sacks. Computational kinematics. Artificial Intelligence, 51:381-416, 1991; Andrew Gelsey. Automated Reasoning about Machines. Ph.D. thesis, Yale University, April 1990. YALEU/CSD/RR No. 785. Each of these approaches has advantages and disadvantages with regard to predictive power, generation of incomplete or impossible behaviors, and computational cost. Detailed numerical simulation provides the most accurate level of detail, but at the highest cost. Since multiple simulations are needed to extract data about trends (e.g. velocity vs. change in a parameter value), precise numerical simulation is often unattractive, particularly at the earliest stages of design.

Simple processes and mechanical systems have been described in a number of qualitative reasoning systems. See, D. S. Weld and J. de Kleer, editors. Readings in Qualitative Reasoning about Physical Systems. Morgan Kaufmann, San Mateo, California, 1990. Most systems handle only the simplest cases of nonlinearity, due to the coarse structure of the qualitative representations. More specialized qualitative reasoning, like trigonometric reasoning, can give more precise results, such as whether or not a particular link in a mechanism can rotate completely with respect to another. It cannot be used to describe specific attributes of the space curves traced by arbitrary points on the mechanism, yet this information is essential to designers.

As the description of the model is more finely discretized, more detail can be obtained from the simulation. Configuration spaces have been used to model a variety of mechanisms. The representation becomes computationally intractable for more complicated non-fixed-axis mechanisms. However, within its domain of applicability, the configuration space approach is able to deal with topology changes during the operation of the mechanism. Simple approximations to dynamic behavior, incorporating models of steady-state forces, allow a number of mechanisms to be simulated at a coarse dynamic level. See, Elisha Sacks and Leo Joskowicz. Mechanism simulation with configuration spaces and simple dynamics. Department of Computer Science Technical Report CS-TR-367-92, Princton University, Princeton, N.J., March 1992.

Detailed numerical simulation, followed by abstracting the results into qualitatively interesting regions, is the most accurate and general approach to describing dynamics in qualitative terms. However, the process can be time-consuming due to long runtimes for general-purpose simulators. Since many simulation runs may be needed to extract data about behavioral trends, this approach can be computationally infeasible.

In the matrix-based approach to dynamics, the kinematic equations and the differential equations are solved together. One approach is to use a maximally redundant set of generalized coordinates (six positional and six velocity per rigid body), and solve using sparse matrix techniques. N. Ofiandea, M. A. Chase, and D. A. Calahan. A sparsity oriented approach to the dynamic analysis and design of mechanical systems--parts I and II. Journal of Engineering for Industry, Trans. ASME Ser. B, 99:773-779, 780-784, 1977. The resulting equations are stiff, and therefore require small integration time-steps. Other approaches attempt to reduce the number of generalized coordinates in order to improve computational efficiency, and to help reduce the stiffness of the matrix. See, Edward Haug. Computer Aided Kinematics and Dynamics of Mechanical Systems, volume 1 Basic Method. Department of Mechanical Engineering, University of Iowa, Iowa City, Iowa, 1985.

General-purpose simulators usually deal with fixed topology mechanisms. A simulator with built-in capabilities for detecting collisions and other topological changes has been proposed. This simulator reformulates the equations when changes occur, and then continues dynamic simulation. See, James F. Cremer. An Architecture for General Purpose Physical System Simulation-Integrating Geometry, Dynamics, and Control. Ph.D. thesis, Cornell University, Ithaca, N.Y., April 1989. Department of Computer Science TR 89-987.

An alternative to detailed numerical solution of the mixed algebraic/differential equations is the Dyne symbolic system which generates a set of pure differential equations, which describe the system in terms of a minimal set of generalized coordinates corresponding to the system's true degrees of freedom. The Dyne system uses symbolic reasoning to derive such equations, guided by a set of algebraic transformation rules and meta-level control rules. Such systems are hard to design, maintain, and debug. However, the equations derived using Dyne are useful to designers performing sensitivity analysis at selected points in the mechanism's behavior, and for describing qualitative regions of behavior. D. R. Brown and L. J. Leifer. The role of meta-level inference in problem-solving strategy for a knowledge-based dynamics analysis aid. Journal of Mechanical Design, Trans. ASME, 113:438-445, December 1991.

Symbolic geometric solution of kinematics problems using degrees of freedom analysis has been made. See, Glenn A. Kramer. Solving Geometric Constraint Systems: A case study in kinematics. MIT Press, Cambridge, Mass., 1992. The present invention is an entirely different approach to the solution of such geometric constraint problems and avoids algebraic reformulation of the geometric entities and constraints. Instead, the present method reasons directly about the degrees of freedom of the geometric entities that are possible.

Geometric entities have degrees of freedom which allow them to vary in location or size. For example, in 3 dimensional space a general fairly rigid body has three translational and three rotational degrees of freedom. A circle with a variable radius has three translational, two rotational, and one dimensional degree of freedom. In the present application the term "configuration" of a geometric entity is considered to be translational, rotational, and dimensional degrees of freedom variables used to parameterize the geometric entity. A detailed discussion of "degrees of freedom" analysis is discussed in U.S. Pat. No. 5, 043,929 commonly assigned to applicant (incorporated by reference for all purposes).

Direct methods for dynamics can be traced to Euler, who advocated treating dynamics by partitioning the problem into two parts: kinematics and kinetics. See, R. S. Hartenberg and J. Denavit. Kinematic Synthesis of Linkages. McGraw-Hill, New York, 1964. Kinematics deals with the positions of the parts of the mechanism as constrained by geometric relationships. Relative (but not absolute) velocities and accelerations often can be calculated kinematically. Kinetics deals with how physical objects move under the effect of forces, and deals with absolute velocity, acceleration, mass, inertia.

Euler also demonstrated how an object in three-dimensional space can be moved from one position to an arbitrary second position by a combination of a single translation and a single rotation, where the rotation axis is parallel to the translation vector, resulting in a screw-like motion. At any instant in time, a body in motion may be thought of as moving about an instantaneous screw in space; this screw's position, orientation, and pitch changes over time. The theory of screws was treated in depth in a geometric manner in R. S. Ball. A Treatise on the Theory of Screws. Cambridge University Press, Cambridge, UK, 1900.

The nonpatent references cited herein are incorporated by reference for background.

It would, therefore, be a significant advance in the art if an efficient method were devised for simulating the dynamics of complex mechanisms, permitting efficient solution of dynamics problems at interactive rates.

SUMMARY

The present invention provides an intelligent computational methodology for performing accurate dynamic analysis of rigid-body mechanisms in an efficient, interactive manner. This methodology allows fast evaluation of design alternatives, and provides information for sensitivity analysis and force analysis. The approach of the present invention separates the problem into kinematics and kinetics, as first proposed by Euler. The kinematic analysis is efficiently computed using: degrees of freedom analysis. Velocity analysis, required for kinetics, is performed using screw theory, formulated as another geometric constraint problem. This analysis allows easy derivation of the Lagrangian form of the dynamics equations in terms of a minimal set of generalized coordinates, thereby resulting in efficient and stable computation methods. In addition, by adjusting the time-step, coarse grained dynamic analysis is possible in an even more efficient manner.

Broadly speaking, the method of dynamic simulation of the present invention automatically formulates a representation (preferably lagrangian) of the dynamics of the mechanism. Kinematic information of the mechanism is determined and used as a potential energy term in the lagrangian representation. The kinetic coenergy term of the lagrangian representation is determined by finding the instantaneous screw axes of at least some of the mechanism components and calculating the angular and linear velocities of at least some of the mechanism components. The angular and linear velocities are ratios of distances of components center of mass to respective screw axes. The dynamic behavior of the mechanism is simulated over time, given starting values of mechanism position and velocities, by solving the lagrangian representation for accelerations and integrating over time.

In a preferred form, the kinematic information is determined as a function of the generalized position of the mechanism components using a degrees of freedom analysis. Preferably, the kinematic information is determined by calculating for one component the change in position of the mass center as a function of the change of the generalized position coordinates for the mechanism. With a known change of position of a component, the change in position of the other components' center of mass can be determined as a function of the ratio of the linear velocities of the components' center of masses.

Preferably the lagrangian representation is a differential equation of the form ##EQU1## where .XI. is the sum of any nonconservative forces; in the case of a conservative force, such as gravity, .XI. is equal to 0. In this lagrangian representation, the first term is the kinetic coenergy term and the second term is the potential (i.e. kinematic) energy term.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow chart describing the method of the present invention;

FIG. 2 is a schematic depiction of a four-bar linkage mechanism;

FIG. 3 is the four-bar linkage mechanism of FIG. 2, illustrating a solution technique for deriving the instantaneous center of rotation for link L.sub.2 ;

FIG. 4 is a schematic describing the method for calculating the derivative of an instantaneous center; and

FIG. 5 geometrically depicts an alternative form for calculating the derivative of an instantaneous center.

DESCRIPTION OF A PREFERRED EMBODIMENT

1 Overview

The present invention formulates dynamic simulation into kinematic (or geometric) and kinetic (or force related) constraints, which allows efficient solution of dynamics problems at interactive rates. The method uses a geometric constraint engine, based on degrees of freedom analysis, for the kinematic solution with extensions to incorporate reasoning about dynamics. Additionally, the method calculates velocity ratios, formulated as a problem in geometry analogizing body displacements to screw axes. These calculations, in conjunction with the kinematic solution, allow determination of all kinetic and potential energy terms required for dynamic simulation. The constructed system of dynamics equations are pure differential equations in terms of a minimal set of generalized coordinates. The formulation results in small, dense matrices, rather than large, sparse ones. In principle, large time-steps may be used for coarse dynamic behavior, with smaller time-steps yielding better approximations to the true behavior, at the cost of more computation time.

The rigid bodies (or components) of the mechanism are thought of as moving about an instantaneous screws in space. With the instantaneous screws known for each body in a rigid-body mechanism, all relative velocities in the system are related by ratios of distances from the appropriate screws. Specifying one absolute velocity then allows finding all absolute velocities. The velocities provide the kinetic coenergy terms in the lagrangian formulation of dynamics, and the kinematic information provides the potential energy terms.

All terms required in the lagrangian are attainable from simple geometric constructions. This leads to the following method for constructing the lagrangian and using it in dynamic analysis as illustrated in FIG. 1.

1. Calculate the kinematic information in terms of the generalized position coordinates. 2. Find the instantaneous screw axes by geometric construction. 3. Calculate angular and linear velocities as ratios of distances between points on the mechanism and the screw axes. 4. Calculate the time derivatives of the screw axes. 5. Use the position and screw information to construct the lagrangian directly. 6. Solve the lagrangian for accelerations. 7. Integrate over time.

Steps 1-6 generate a plan which is iteratively reused to produce new values for accelerations which are integrated over time (step 7), resulting in dynamic simulation.

2 Example

FIG. 2 describes a four-bar planar linkage as an example of a simple rigid-body mechanism. One of the links is Founded, or fixed to the global reference frame, and is not shown. The moving links (i.e., components or bodies) have lengths l.sub.1, l.sub.2, and l.sub.3, and are modeled with masses m.sub.1, m.sub.2, and m.sub.3 and rotational inertias I.sub.1, I.sub.2, and I.sub.3. The revolute joints j.sub.i and j.sub.4 connect the cranks (links l.sub.2, l.sub.3) to the ground; the remaining two joints j.sub.2, j.sub.3 connect the cranks (links l.sub.1, l.sub.3) to the coupler (link l.sub.2). The linkage has one degree of freedom; the position of all links are determined fully by the crank angle .THETA., and the velocities by the first derivative .THETA. of with respect to time .THETA..

The example problem is an initial values problem. Given initial values for .THETA. and .THETA., find the behavior of the linkage over time. In this problem, gravity (a conservative force) exerts a downward force, as illustrated, and the joints are assumed to be frictionless. Other acceleration forces may, of course, be applied to the mechanism.

The first portions of the method require: 1 ) calculating the: kinematic information as a function of the generalized position coordinates; 2) finding the instantaneous screw axes by geometric construction; and 3) calculating angular and linear velocities as ratios of distances of mass centers from the screw axes. The first item is covered by kinematic analysis U.S. Pat. No. 5,043,929 (incorporated by reference).

In planar mechanisms, the axes of all instantaneous screws are normal to the plane, and the pitch is always zero. In this specialized case, the instantaneous screw is known as the instantaneous center of motion. FIG. 3 illustrates how the instantaneous centers and velocities are calculated. The instantaneous centers of the cranks (links l.sub.1, l.sub.3) are the fixed joints (j.sub.i for link l.sub.4 and j.sub.4 for link l.sub.3); these centers do not change over time.

The instantaneous center for link l.sub.2 is computed geometrically. Since joint j.sub.2 is on the link l.sub.2, it must rotate about the instantaneous center for link l.sub.2. However, since j.sub.2 is also on link l.sub.1, it must rotate about j.sub.1, and hence, its linear velocity must be perpendicular to the line describing link l.sub.1. Therefore, the instantaneous center must lie on a line through j.sub.1 and j.sub.2. Using a similar argument for the motion of joint j3, the instantaneous center for link l.sub.2 is found by intersecting the two lines describing l.sub.1, l.sub.3. Link l.sub.2 rotates about the intersection point with angular velocity .omega..sub.2, which is yet to be determined. Similar constructions exist for general screws. See, R. S. Ball. A Treatise on the Theory of Screws. Cambridge University Press, Cambridge, UK, 1900.

The angular velocity of link l.sub.1, .omega..sub.1, is specified to be .THETA., the first derivative of the crank angle .THETA. with respect to time. The angular velocity .omega..sub.2 of link l.sub.2 is found by equating the linear velocities at joint j.sub.2 : .omega..sub.1 l.sub.1 =.omega..sub.2 r.sub.1. In similar fashion, the remaining angular velocities are calculated:

.omega..sub.1 =.THETA. .omega..sub.2 =-.omega..sub.1 (l.sub.1 /r.sub.1)=-.THETA..sub.1 (l.sub.1 /r.sub.1) .omega..sub.3 =-.omega..sub.2 (r.sub.3 /l.sub.3)=.THETA..sub.1 (l.sub.1 r.sub.3 /l.sub.3 r.sub.1)

The directions of the linear velocities are derived from the angular velocities .omega..sub.1, .omega..sub.2, .omega..sub.3. The magnitudes of the linear velocities (assuming each center of mass is at the center of the link) follow directly from the angular velocities:

.vertline.v.sub.1 .vertline.=.omega..sub.1 (l.sub.1 /2) =.THETA..sub.1 (l.sub.1 /2) .vertline.v.sub.2 .vertline.=.omega..sub.2 r.sub.2 =.THETA..sub.1 (l.sub.1 r.sub.2 /r.sub.1) .vertline.v.sub.3 .vertline.=.omega..sub.3 (l.sub.3 /2)=.THETA..sub.1 (l.sub.1 r.sub.3 /2r.sub.1)

Thus, all angular and linear velocities are calculated as simple ratios of distances between instantaneous centers and points on the mechanism. For a more compact notation, the ratios are written as .omega..sub.i =.THETA.k.sub.i, and .vertline.v.sub.i .vertline.=.THETA.j.sub.i.

The lagrangian of the system is generated as follows (assuming a positive y axis moving upward in FIG. 3, and with v.sub.i =.vertline.v.sub.i .vertline.. ##EQU2##

Here, I.sub.i is the moment of inertia for mass link i, and g is the acceleration of gravity. An alternative representation is the lagrangian differential equation describing the motion of the linkage as: ##EQU3##

The second term in this lagrangian differential equation is the potential energy term, and is calculated as follows: ##EQU4##

The change in height (y) for each point mass is: ##EQU5##

where .delta.p.sub.1 /.delta..THETA. is the change in position of point mass m.sub.i with respect to a .THETA.. This vector is always in the direction of v.sub.i. The value of .delta.p.sub.1 /.delta..THETA. is calculated from .THETA. as follows: ##EQU6##

The remaining position change magnitudes are related to each other directly as the magnitude of the linear velocities, which have already been computed geometrically: ##EQU7##

The first term in the lagrangian differential equation is the kinetic coenergy term, and is calculated as follows. The derivative with respect to .THETA. is: ##EQU8##

The time derivative of this quantity is: ##EQU9##

The first term of this time derivative involves the lengths of geometric entities already constructed. Evaluating the second term involves understanding how the instantaneous centers move over time. In this example, the only instantaneous center that moves over time is the center for link l.sub.2.

Using the chain rule, finding the time derivatives of the j.sub.i 's (as well as the k.sub.i 's) is reduced to finding the derivatives of these quantities with respect to .THETA.: ##EQU10##

In the general case, the differential dj.sub.i is the sum of the partials. Since this is a one variable problem, (i.e. .THETA.) the total and partial differentials are ,equivalent.

The partial derivative .delta.j.sub.i /.delta..THETA. may be found by differential analysis of the construction lines that are intersected to form the instant center of rotation. FIG. 4 illustrates the process. The figure shows two lines, labeled 1 and 3, which intersect in point p.sub.2 ; they are rotated to new positions 1' and 3' respectively, where they intersect in a new point, distance .delta.p.sub.2 from p.sub.2. The lines are rotated by amounts .delta..THETA..sub.1 and .delta..THETA..sub.3, respectively.

In the case where p.sub.1 and p.sub.3 are the grounded joints of the four-bar linkage, and point p.sub.2 is the instantaneous center of link 2, the length of line segment p.sub.1 p.sub.2 is l.sub.1 +r.sub.1, and the length of line segment p.sub.3 p.sub.2 is l.sub.3 +r.sub.3. Since l.sub.1 and l.sub.3 are parameters of the linkage, and constant, the value of .delta.p.sub.2 yields r.sub.1 and r.sub.3 trivially. The only remaining detail is to determine .delta..THETA..sub.3 in terms of .delta..THETA..sub.1, so only one independent quantity is used when taking the derivative. The ratio of the angular velocities provides this information: ##EQU11##

In FIG. 4, the lines have been positioned to simplify the equations. Since line 1 passes through the origin, its equation is y=xtan.THETA..sub.1. Line 3 passes the x axis at x=d, so the equation is y=(x-d) tan.THETA..sub.3. The intersection point is then: ##EQU12## The derivatives with respect to .THETA..sub.1 are calculated using the known relationship between .delta..THETA..sub.1 and .delta..THETA..sub.3 : ##EQU13## Then, ##EQU14## is found as the vector sum of the two above quantifies.

The quantity .delta.j.sub.i /.delta..THETA. may alternatively be found through geometric sensitivity analysis, involving simple geometric constructions. FIG. 5 illustrates the technique.

Consider three points p.sub.1, p.sub.2, and p.sub.3 where points p.sub.1 and p.sub.3 are fixed in space. If the line through p.sub.1 p.sub.2 is rotated through a small displacement .delta..THETA. and p.sub.2 is constrained on line p.sub.3 p.sub.2, it will move in direction .delta.p.sub.2,1. If, on the other hand, line p.sub.3 p.sub.2 is rotated by a small displacement .delta..THETA..sub.3, p.sub.2 will move in direction .delta.p.sub.2,3. For a composite change of both .delta..THETA., and .delta..THETA..sub.3 the movement of point p.sub.2 will be the, vector sum of the two independent movements (for infinitesimal displacements, the quadrangle on which the three vectors lie becomes a parallelogram).

In the case where p.sub.1 and p.sub.3 are the grounded joints of the four bar linkage, and point p.sub.2 is the instantaneous center of link l.sub.1, the length p.sub.1 p.sub.2 is l.sub.1 +r.sub.1 and is l.sub.1 +r.sub.3. Since l.sub.1 and l.sub.3 are parameters of the linkage, and constant, the .delta.p.sub.2 value of yields .delta.r.sub.1 and .delta.r.sub.3 trivially. The only remaining detail is to determine .delta..THETA..sub.3 in terms of .delta..THETA..sub.1 so only one independent quantity is used when taking the derivative. The ratio of the angular velocities provides this information: ##EQU15##

Many types of integration schemes may be used, but since the system of equations is purely differential, with no algebraic equations, stiffness should not be an overriding concern. Traditional integration schemes, such as adaptive step size Runge-Kutta and predictor-corrector methods, are believed to be preferred.

If integration uses large time-steps, the effect should be to have "approximate" dynamics. There is no danger of the mechanism "flying apart," since the kinematic constraints are not considered in the time integration. There is a trade-off between computation time, accuracy, and interactiveness.

3 Discussion

Generating the plan of geometric constructions to find velocities need only be done once for a given mechanism. After that, the plan may be reused during each step of the dynamic simulation. This approach of plan generation and reuse is fully explained in U.S. Pat. No. 5,043,929.

Since there are as many instantaneous screws as there are bodies in the mechanism, evaluation of the kinetic coenergy terms of the lagrangian takes time linearly proportional to the size of the mechanism. Kinematic analysis, necessary for the potential energy terms, is O(n log n), but typically linear in n, where n is the number of bodies in the mechanism.

Solving the lagrangian for the accelerations requires inverting a matrix of size d, the number of true degrees of freedom in the mechanism. This contrasts with the standard matrix-based approaches, where the matrix to be inverted is of size proportional to the number of bodies in the mechanism. In the worst case of a mechanism comprised exclusively of open chains, the number of true degrees of freedom will be proportional to the number of bodies; however, the absolute number of generalized coordinates being considered will still be less using the geometric algorithm, since the kinematic constraints are already eliminated.

Beneficially, the geometric constructions required for the velocity analysis can all be performed in a highly stylized fashion. If the gains in computational efficiency are comparable with the gains in kinematic simulation, the speedup in dynamic simulation could be substantial, affording interactive simulation speeds for many complex mechanisms.

While only mechanisms with fixed topology are currently simulated, known techniques should permit simulation of mechanisms with changeable topology.

Besides their use in formulating the dynamics equations, the geometric sensitivities can also be used for force analysis at selected points in the mechanism's trajectory, and for kinematic analysis of velocity ratios. These are quantities that could be optimized in a design at interactive rates.

If the dynamic behaviors can be simulated efficiently enough, it may be possible to make multi-dimensional "maps" of simulated behavior as a function of design parameter values. Such behavioral maps could be of substantial benefit in the design and debugging of complex mechanical devices.

Claims

1. In a digital computer system a method of dynamic simulation of a mechanism having a plurality of components, comprising the steps of:

a. determining kinematic information of the mechanism;

b. finding the instantaneous screw axes of at least some of the mechanism components;

c. calculating the angular and linear velocities of at least some of the mechanism components as ratios of distances of component center of mass to respective screw axes;

d. determining kinetic coenergy terms of the mechanism using the angular and linear velocities;

e. formulating a representation of the dynamics of the mechanism using the kinetic coenergy terms, and the kinematic information as the potential energy term;

f. by integrating the accelerations over time; and

g. simulating the dynamic behavior of the mechanism over time, given the starting values of position and velocities.

2. The method of claim 1, the representation of the dynamics being a lagrangian representation.

3. The method of claim 1, the representation being a lagrangian representation of the form: ##EQU16## where the first term is the kinetic coenergy term, the second term is the potential energy term, and .XI. is the sum of any nonconservative forces.

4. The method of claim 1, the kinematic information being determined as a function of generalized position coordinates describing the mechanism components.

5. The method of claim 4, including determining the potential energy term using the kinematic information including calculating for one component the change in position of the mass center as a function of the change of the generalized position coordinates for the mechanism.

6. The method of claim 4, including determining the potential energy term using the kinematic information where the change in position for a component is: ##EQU17## .DELTA.P.sub.1 =Change in magnitude of position of component 1's center of mass relative to a change in some generalized position coordinate.

v.sub.1 =linear velocity of component 1's center of mass

v.sub.2 =linear velocity of component 2's center of mass

.DELTA.P.sub.2 =known change in magnitude of position of component 2.

7. The method of claim 3, where the kinetic coenergy term is determined as follows, the derivative with respect to .THETA. is: ##EQU18## and the time derivative of this quantity is: ##EQU19## where the j's and k's are ratios of distances to instantaneous screw axes.

8. The method of claim 7, where the time derivatives of the j's and k's, ##EQU20## are reduced to the derivatives of the j's and k's with respect to position .THETA.: ##EQU21##

9. The method of claim 1, including compiling steps (a-g) into a plan, said plan being iteratively reused during simulation.

10. In a digital computer, a system for simulating the dynamics of a rigid-body mechanism, comprising:

a kinematic simulator for determining a potenial energy term for the mechanism using a degrees of freedom analysis in terms of the generalized position coordinates of the mechanism;

means for determining a kinetic coenergy term using information on the instantaneous screw axes for at least some of the bodies of the mechanism;

a lagrangian constructor for formulating a lagrangian representation of the dynamics of the mechanism using the potential energy term and kinetic coenergy term;

a plan generator which solves the lagrangian representation for accelerations to build a plan which simulates the dynamics when integrated over time; and

a video display for showing said simulation of the dynamics of the mechanism over time.

11. The system according to claim 10, including an integrator for iteratively reusing the plan to produce new accelerations which are integrated over time, simulating the mechanism dynamics over time.

12. The system of claim 10, the kinetic coenergy term determining means operable for calculating angular and linear velocities as ratios of distances between points on the mechanism and the screw axes.

13. The system of claim 12, the kinetic coenergy term determining means operable for calculating the time derivatives of the screw axes.