Computerized generation of motion equations using variational graph-theoretic methods

Applied Mathematics and Computation 192 (2007) 135–156 www.elsevier.com/locate/amc

Computerized generation of motion equations using variational graph-theoretic methods M.J. Richard a

a,*

, J.J. McPhee b, R.J. Anderson

c

Department of Mechanical Engineering, Laval University, Quebec, Que., Canada G1K 7P4 Systems Design Engineering, University of Waterloo, Waterloo, Ont., Canada N2L 3G1 c Department of Mechanical Engineering, Queen’s University, Kingston, Canada K7L 3N6 b

Abstract Severe tolerances on mechanical components have created increasingly stringent demands on the quality of new mechanical designs. The mathematical models used to simulate the various types of mechanical systems these days need to incorporate an optimization algorithm capable of minimizing the number of matrix multiplications when deriving symbolically the equations of motion. The method is based on a simplistic topological approach which is incorporated into an efficient variational graph-theoretic process used to solve these non-linear problems. The system is represented by a linear graph, in which nodes represent reference frames on rigid bodies, and edges represent components that connect these frames. By selecting a proper spanning tree for the graph, the analyst can choose the set of coordinates appearing in the final system of equation. The procedure casts, simultaneously, the Lagrange’s equations of motion and the kinematic constraints into a symmetrical format which yields a symbolic solution. The algorithm serves as a basis for a computer program which generates the equations of motion in symbolic form, and provides the time varying response of the system. The effectiveness of this approach is demonstrated in the analysis of a spatial four-bar mechanism and an articulated semitrailer vehicle. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Simulation; Dynamics; Multibody; Computation

1. Introduction The refinement of formulation techniques for dynamic analysis of mechanisms has led to the development of general-purpose computer programs. Several have been described in the literature [1–9]. All are capable of automatically generating and numerically integrating the differential equations governing the motion of dynamic systems but they do not incorporate automatically an optimization algorithm for selecting a proper set of coordinates. The widespread interest in multibody dynamics is evidenced by the existence of a large

*

Corresponding author. E-mail address: [email protected] (M.J. Richard).

0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.02.135

136

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

number of commercial software packages. In the past 30 years, many researchers have been working on the development of algorithms for three-dimensional systems and several programs were discussed in the survey paper by Paul [1] and book by Schiehlen [10] and in recent years more researchers have been attracted to this problem and there are probably several fairly well-known methods (and computer programs based on their methods), including DRAM and ADAMS [11–18], IMP [19–22], DYMAC [1,23–26], DADS [27–30], NEWEUL [31–38], COMPAMM [39–41] and of course programs based on graph-network models [42–52]. Recently, a paper for selecting the ‘‘best’’ coordinates for a multibody system, using graph-theoretic concepts has been written [53]. This paper demonstrates that automated coordinate selection is possible and has actually been implemented in the most recent version of the commercial software package DynaFlexPro (sold by Maplesoft). This software generates symbolic equations for three-dimensional multibody (and electrical) systems. In a subsequent publication, Schmitke and McPhee [54] present a unified formulation capable of systematically generating the governing symbolic equations for multibody systems. The formulation is based on the principle of orthogonality, a powerful concept that serves as a generalization of the principle of virtual work and virtual power. In this work, the methodology is based on a graph-theoretic formulation which is different from the other methods described in the literature and presents some advantages: (i) the model uses a simple modeling framework to discretize a complex topology of interconnected rigid bodies which is a simple isomorphism of the mechanical system, (ii) the procedure may prove to be helpful due to its simplicity and fundamental application of the variational form of Lagrange’s equations of motion, (iii) in comparison with work done in the field of dynamic systems, this approach performs the minimum number of vector multiplications for the automatic generation of the equations of motion in symbolic form. To be applicable to a wide range of spatial mechanical systems containing both open and closed kinematic chains, a multibody formulation must incorporate general mathematical methods for representing both the system topology as well as the time-varying configuration. The representation of topology is naturally handled using elements of graph theory. However, selecting automatically a proper set of coordinates to represent the changing system configuration is a difficult task because it is a function of the system on hand. In a manual derivation of the equations of motion, one can select an independent set of generalized coordinates by eliminating the forces and torques corresponding to non-working kinematic constraints. In contrast, in an automated formulation, one needs to exploit a pre-defined set of coordinates that will be applicable to the configuration of all multibody systems. In past algorithms, either cartesian (absolute) or joint (relative) coordinates were used. Absolute coordinates define the position and orientation of each body in the system relative to a global reference frame, while joint coordinates define the relative motion of two adjacent bodies in a kinematic chain. The choice of coordinates is crucial for establishing a set of motion equations, and therefore, affects the computational time for solving them. Recognizing the importance of a proper set of coordinates, several researchers have investigated the use of coordinates other than absolute or joint. Fayet and Pfister [55] introduce the concept of indirect coordinates, which define the relative motion of two bodies that are not necessarily adjacent in an open-loop kinematic chain. Huston et al. [56] shows that the use of absolute rotational coordinates leads to a simplification in the equations of motion, when compared to those obtained using joint coordinates. Agrawal and Shelly [57] show the procedure for switching from one set of joint coordinates to another during a dynamic simulation. Finally, there are many formulations that use a coordinate transformation to reduce the differential algebraic equations of motion to a minimal set equal in number to the degrees of freedom. However, these formulations require the transformation equations to be either specified by the user or approximated using numerical methods. It has been shown [58] that a set of coordinates called ‘‘branch coordinates’’ encompasses both sets of cartesian and joint coordinates as special cases. Also, in a previous paper [59], the use of virtual work has been proposed and validated as a new graph-theoretic variable. In order to create a system graph that results in correct kinematic and dynamic equations for any choice of spanning tree, it is necessary to introduce a dependent virtual work element. Hence, by exploiting the principle of virtual work, it is possible to develop a variational graph-theoretic formulation in terms of branch coordinates capable of automatically generating the motion equations for rigid multibody systems. It leads to a smaller number of equations than would generally be obtained using either absolute or joint coordinates.

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

137

In the next section, our graph-theoretic method is reviewed, including a section that suggests an order for selecting trees that leads to the greatest reduction in computer multiplications for the generation of motion equations. A new constitutive or ‘‘energy terminal’’ equation for the dependent virtual work element and other components in the system will also be given. This is followed by a description of the systematic formulation procedure for three-dimensional multibody systems, as well as a symbolic computer implementation of this algorithm. Finally, two examples are presented to illustrate this new variational graph-theoretic formulation and its ability to generate directly a smaller set of motion equations without additional user input. 2. Graph-theoretic representation Many researchers have studied the theory of graphs [60–66] and bond-graphs [67–69] mainly due to the fact that among all the fields of human interest, there are few where graph theory cannot be applied to the process of analyzing or synthesizing problems. In order to extract the kinetic properties resting within lumped-parameter mechanical systems, it is convenient to discretize the system into a schematic diagram composed of nodes or vertices representing points of interconnection in the system and oriented edges identifying system elements. Combination of all the vertices delimiting the network of elements with the total set of spanning edges between appropriate nodes will result in a diagram which is a simple isomorphism of the mechanical system. One of the most appealing features of graph-theoretic methods lies in the geometric and pictorial aspect of the method. Given a spatial mechanical system, one can construct the system’s diagram by a simple mapping of the mechanism. For instance, the vertices would correspond to centres of mass for rigid bodies, points on bodies to which forces are applied or joints are connected, and a ground node that represents the origin of an inertial reference frame. Each element is represented by a line segment and each joint or connection by an appropriate point such that a user can associate the network diagram to the mechanical system in a direct fashion. The technique is very methodical and well suited for computer implementation. Much of the simplicity and efficiency of graphical methods lie in the use of a ‘‘tree’’ to assist in arranging the order of computation. By definition, a tree is a subgraph where every vertex of the graph is connected by exactly one chord. Plainly, this connotes that the subgraph is connected and contains all the nodes of a given system graph, but has no closed loops. Thus, a tree is considered a minimal connected graph in which the deletion of a single branch would separate the subgraph. The set of vectors that complement the tree are called ‘‘cotree’’. However, given a graph with m vertices and  edges, the order of interconnection of a system can be summarized in a m by  incidence matrix. It is easy to construct since each edge is adjacent to exactly two vertices. The incidence matrix of Gðm; Þ is denoted by ½j and is defined as follows: (1) each of the m rows corresponds to a vertex of G, (2) each of the  columns corresponds to an edge of G, and (3) entries jij ¼ 1 if the ith node is the initial vertex of the jth edge, jij ¼ þ1 if the ith node is the final vertex of the jth edge, and jij ¼ 0 otherwise. Note that the incidence matrix can always be compiled from inspection of the graph. It has been shown [58] that the components necessary to generate an optimum translational tree for branch coordinates should be selected in the following order: N 1: N 2: N 3: N 4: N 5: N 6:

Element Element Element Element Element Element

representing representing representing representing representing representing

rigid arms; position drivers; spherical and revolute joints; cylindrical and prismatic joints; rigid bodies; force actuators.

Let us examine the spatial four-bar mechanism represented in Fig. 1. The system consists of only two rigid bodies, a crank (body 1) of mass m1 and a coupler link (body 2) of mass m2 connected with a revolute joint located at the center of mass of body 1 at one unit from the reference frame. The length of body 2 is l2 with its centroid located at the center. A link with unit length and negligible mass (no inertia) connects, with spherical joints, the end of body 2 and point (0, 2, 0) on the y-axis as shown in Fig. 1. An x01 y 01 z01 frame is fixed in the first link at centroid m1, with its x01 axis along the first bar and its z01 axis parallel to the global Z-axis. Gravity is

138

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156 Z

(0, 2, 0) Y

T 1

1

z'2

z'1

y'1 2

α1

X

α2

m1

y'2

m2 x'2

x'1

Fig. 1. Spatial four-bar mechanism.

acting in the negative-Z direction. The variables a1 and a2 represent angular displacements of the two revolute joints. Finally, a driving torque T causes the crank to oscillate. The linear graph representation of the spatial four-bar mechanism is depicted in Fig. 2. Edges e11, e12 and e13 correspond to rigid-arm elements which represent the position and orientation of a reference frame at the distal end of the link. Edges e31, e32 and e33, e34 represent revolute joints and spherical joints, respectively. Edges e51 and e52 represent the two rigid bodies relative to the inertial frame (datum node O) and edge T is the driving torque. Edge e35 represents a constrained rigid link and edge e61 corresponds to the gravitational force. A set of edges, e71–e77 not shown, corresponding to dependent virtual work elements [70] are automatically added to the system graph to generate the correct virtual topological equations at every vertex. It is well known that the general motion of a rigid body can be decomposed into independent translational and rotational motions. The linear graph in Fig. 2 represents strictly the translational motion of elements in the system. A separate graph, not shown, is therefore needed to represent the rotational motion of the four-bar mechanism components. It is obvious that the rotational graph is simpler than the translational one since the angular velocity is a characteristic of the whole body and not an individual point. Thus, all points on a rigid

T

0 e21 e31 e34 e11

e51

e52

e61 e35

e12

e13 e33

e32 Fig. 2. Graph of the spatial four-bar mechanism.

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

139

body are assembled to a single node, thereby eliminating all body-fixed vectors. The only elements appearing in the rotational graph are the two bodies and the joints. Unlike the traditional approach to dynamics in which Newton’s second law is the basic postulate, in the graphical method there are two equally important fundamental postulates: the vertex and circuit postulates. The vertex postulate states that the algebraic sum of through-variables fTVg, symbolizing quantities such as force F, torque T and virtual work dW , corresponding to all the vectors incident with any vertex of the graph is, identically, zero. Essentially, this is recognizable as the dynamic force-balance law or d’Alembert’s principle which requires that force summation upon each body, including d’Alembert’s inertial variable, must be equal to zero. However, in this form it is more general, since it applies to any physical system; for example, in electrical systems, it is equivalent to Kirchhoff’s current law. A cutset isolates a part of the system and is equivalent to the construction of a free-body-diagram. From graph theory it can be shown [43] that these cutset equations are obtainable by simple row operations on the incidence matrix, and can be written in the form: ½U t

DfTVg ¼ 0;

ð1Þ

where ½Dij  is a sub-matrix containing +1, 1 and 0 depending whether element Nj is incident to Ni and submatrix Ut is a unit matrix associated tree elements. The column matrix fTVg represents the through-variables (Fi, Ti or dW i ) applied by element Ni. Consider the four-bar mechanism, the vertex equation for the node representing the end point of rigid arm e12 is TV32 þ TV73  TV12 ¼ 0

ð2Þ

and if we use virtual work as through variables, dW 32 þ dW 73  dW 12 ¼ 0;

ð3Þ

where dW 32 ¼0, since a frictionless joint cannot add or remove energy from the system. Another governing postulate is the circuit postulate which states that the algebraic sum of across-variables fAVg, representing those quantities such as displacements r, velocities v, accelerations v_ , virtual angular dis_ corresponding to all the vectors included in placement dh, angular velocities x and angular accelerations x, any circuit is, identically, zero. Basically, this postulate represents the geometrical relations guiding the motion of mechanical systems. To be precise, each circuit equation alludes to a closed vector polygon respecting the geometric fit or compatibility law of rigid body dynamics. From graph theory, it can be shown that these kinematic constraint equations can be obtained from the cutset equations and are usually written in the form: ½E

U c fAVg ¼ 0;

ð4Þ

where sub-matrix ½Eji  specifies which element Ni is included in the closed loop Nj and sub-matrix Uc is a unit matrix associated with cotree elements. The column matrix fAVg represents the across-variables (ri, vi, v_ i , dhi , xi and x_ i ) for element Ni. Note that the submatrices [D] and [E] are related by the principle of orthogonality [43,54] (not well-known in dynamics): ½E ¼ ½D

T

T

or ½D ¼ ½E :

ð5Þ

Since there is one independent circuit equation for each chord in the graph, consider the circuit equation for edge e52, AV52  AV12  AV32  AV11  AV31 ¼ 0:

ð6Þ

Naturally, each across-variable may be replaced with the virtual displacement dr or even the virtual rotation dh since these vectors obey the rules for vector quantities and may be so treated. However, since the acrossvariables cannot be substituted for finite rotations because they do not obey the commutative law of vector addition, one must define an equivalent circuit equation in terms of rotation transformation matrices ½Rð33Þ associated with each element. For example, the corresponding equation for edge e52 can be expressed as T

T

T

T

½R52 ½R12  ½R32  ½R11  ½R31  ¼ ½1;

ð7Þ

140

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

where [1] is a unit matrix. Note that the additive nature of the vector circuit equation is replaced by a multiplicative expression in its rotational counterpart. Making use of the orthogonal nature of rotation transformation matrices, ½R52  can then be isolated as a branch transformation ½R52  ¼ ½R31 ½R11 ½R32 ½R12 :

ð8Þ

Each of the rotation matrices will be functions of the rotation variables for the corresponding element. Having chosen a complete set of branch coordinates for a system with j degrees of freedom (DOF), this procedure may specify j branch coordinates as primary coordinates. The remaining position variables are called secondary coordinates and may be computed from the branch coordinates. In this formulation, the primary variables are the tree across and cotree through variables and the secondary variables are the tree through and cotree across variables. Note that the secondary variables can be expressed as functions of the primary variables using the chord and branch transformation equations (1) and (4). For a closed-loop system, like the four-bar mechanism, one may select a subset of the angles to serve as branch coordinates. Note that the associated constraint equations frequently cannot be solved explicitly for the remaining secondary coordinates. Even when they can be found, the expressions for the secondary variables in terms of primary coordinates will be very cumbersome for most closed-loop systems, and it will be difficult to write the transformation equations. Therefore, it seldom pays, for computational purposes, to work exclusively with generalized coordinates when dealing with closed-loop systems and hence, the retention of secondary variables can be very advantageous. 3. Virtual work terminal equations The basic premise of the graphical approach is that each system element is modelled separately by defining its characteristics, independent of the other system elements, in the form of a ‘‘virtual terminal’’ (or constitutive) equation. Associated with every edge is one or more terminal equations that define the generalized force Q of an element corresponding to its branch coordinate q. These terminal equations are written in terms of the system through and across variables. In this formulation, a translational and rotational network is required in order to ensure consistency in the initial theory. Hence, these equations are functions of virtual displacements and virtual work. In this formulation a rigid-arm element is used to represent the relative position and orientation of two reference frames on the same body. This second reference frame may be needed to locate a point of application for some other element. The terminal equations, in terms of virtual displacements, for a rigid-arm element are dh1 ¼ 0; dr1 ¼ dh5  r1 ;

ð9Þ ð10Þ

where r1 is the rigid-arm position vector which is function of rotational body-fixed frame h5 since the direction of r1 varies as the rigid body rotates. If the virtual work is the work done by specified forces on virtual displacements which are consistent with the constraints, with all other coordinates being kept constant, then the virtual work of force F is dW ¼ drT F:

ð11Þ

Now, let us consider a system whose kinematics are parameterized by a branch coordinate vector q where r ¼ rðqÞ, then a virtual displacement dr is related to a branch coordinate variation dq by dr ¼ rq dq;

ð12Þ

where the subscript q indicates a partial derivative with respect to vector q. Then the virtual work of a force may be obtained by substituting Eq. (12) in Eq. (11), dW ¼ dqT rTq F:

ð13Þ

Eq. (13) may be interpreted as the scalar product of a branch variation dq and its generalized force Q, dW ¼ dqT Q;

ð14Þ

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

141

where the generalized force is defined as Q  rTq F. In this formulation, one may separate forces into applied forces FA and constraint forces FC such that, dW ¼ dqT ðQA þ QC Þ ¼ dW A þ dW C :

ð15Þ

If branch coordinates are consistent with constraints, then q is a vector with independent coordinates. Note that if constraints that act on the system are workless and if the branch coordinates q must satisfy constraint equations of the form Uðq; tÞ ¼ 0;

ð16Þ

they are dependent. Thus branch coordinate variations must satisfy Uq dq ¼ 0;

ð17Þ

where Uq is the Jacobian matrix of the constraint equation (16). If constraint forces are workless, it can be proven that dW C ¼ dqT QC ¼ dqT UTq k;

ð18Þ

where the Lagrange multiplier k represents vector-generalized constraint forces. Up to now, we have assumed that the reaction forces due to constraints neither produce nor consume work. This restriction has prohibited us from analyzing systems with non-rigid links (like elastic bodies which as already been treated [70]) or systems with energy-dissipating devices such as viscous or Coulomb dampers. The effect of such constraints is to produce certain pairs of internal forces which may be treated exactly as we would any other applied forces. Hence, the new form of terminal equations for applied forces has been written in the first part of Eq. (15). For elements containing physical constraints the terminal equation (18) is introduced. Care must be taken when applying the constraint generalized force because the physical constraints do no work and its Lagrange multiplier form must be introduced in the initial cutset equations at the beginning of the substitution procedure. Spring and gravitational forces belong to the wider class of so-called conservative forces for which the generalized force can be efficiently calculated from the potential energy. Thus, if forces that act on a system can be expressed as the gradient of a scalar function of generalized coordinates, the virtual work may be written as the negative of a variation in potential energy V and the new terminal equation for conservative forces becomes dW cons ¼ dqT Qcons ¼ dqT V Tq :

ð19Þ

Now, let us consider a moving system defined by a set of branch coordinates. If between these branch coordinates and time t there exists some relation of the form of Eq. (16), it is said that the system is moving under constraint. This means that the functions of constraints are geometric or kinematic conditions which restrain the possibilities of motion of the system. For a spatial mechanical system there is a limited number of types of functions of constraint, represented by the joints between the bodies. Fig. 3 presents the four most common types of ideal joints which can be found in multibody systems. A spherical joint shown in Fig. 3a is defined by the condition that the center of the ball coincides with the center of the socket. The terminal equations are, then, three scalar constraint equations that restricts the relative positions between the bodies, ½ ex

ey

ez T dr3 ¼ ½ 0

0

0 T ;

ð20Þ

where ex ; ey ; ez are unit vectors and dr3 represents the relative virtual displacement between the bodies. Essentially, the three constraining equation (20) impose that there is no relative displacements at the joint in the x  y  z directions. A revolute joint, shown in Fig. 3b, is constructed with bearings that allow relative rotation about a common axis in a pair of bodies, but precludes relative translation along this axis. The five terminal constrained equations are ½ ex

ey

ez T dr3 ¼ ½ 0

0

0 T ;

½ ex

ey T dh3 ¼ ½ 0

0 T ;

where dh3 represents the relative virtual angular rotation between the bodies.

ð21Þ

142

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

Fig. 3. Ideal kinematic joints (a) spherical, (b) revolute, (c) cylindrical and (d) prismatic.

A cylindrical joint between a pair of bodies is shown in Fig. 3c. It permits relative translation and relative rotation between bodies about a common axis. The four terminal constrained equations are ½ ex

T

ey  dr4 ¼ ½ 0

T

0 ;

½ ex

T

ey  dh4 ¼ ½ 0

T

ð22Þ

0 :

A prismatic joint, shown in Fig. 3d, allows relative translation along a common axis between a pair of bodies, but precludes relative rotation about this axis. The five terminal constrained equations are ½ ex

T

ey  dr4 ¼ ½ 0

T

0 ;

½ ex

ey

T

ez  dh4 ¼ ½ 0 0

T

0 :

ð23Þ

Let us now consider rigid body elements N5. The graph for this type of element consists of an edge e5 from the datum node to a local reference frame on the body. The use of the variational form of Lagrange’s equation relies upon a correct formulation of the kinetic energy of the multibody system in terms of branch coordinates. If the system is composed of nb rigid bodies, the kinetic energy of the ith body is defined as 1 1 T i ¼ mi vTi vi þ x0T ½I 0 x0 ; ð24Þ 2 2 i i i where mi is the mass of body i, vi is the velocity vector of the centre of mass of the body in inertial space, ½I 0i  is the inertia tensor of the body expressed in a body-fixed 0 coordinate system and x0i is the angular velocity vector of the body in inertial space and expressed in the inertia 0 tensor coordinate system. Note that absolute or inertial-space velocities must be used although they may be expressed in any convenient (inertial or non-inertial) coordinate system. Then, the virtual work terminal equation for the ith body becomes (   T ) d oT oT i i dW 5i ¼ dqT : ð25Þ  dt oq_ oq It is assumed at this point that the velocities {vi} and {xi} have been found in symbolic form for each body component in terms of branch coordinates and their derivatives. In order to apply the variational form of

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156 i Lagrange’s equation, the terms dtd ðoT Þ and oq_ entiation of Eq. (24),  oT i  T 0 0 ¼ mi vi viq þ x0T i ½I i xiq oq

oT i oq

143

will be required for each branch coordinates. By partial differ-

and similarly,  oT i  T 0 0 ¼ mi vi viq_ þ x0T i ½I i xiq_ : oq_ The time derivative of Eq. (27) is     d oT i 0 0 0T 0 _ 0iq_ : ¼ mi v_ Ti viq_ þ mi vTi v_ iq_ þ x_ 0T i ½I i xiq_ þ xi ½I i x dt oq_ Substituting Eqs. (26) and (28) into the variational form of Lagrange’s equation (25) gives h iT x0 0 0 T v 0T 0 dW 5i ¼ dqT mi v_ Ti viq_ þ x_ 0T ½I x þ m v P þ x ½I P ; i i i=q i i iq_ i i i=q

ð26Þ

ð27Þ

ð28Þ

ð29Þ

where Pvi=q ¼ v_ iq_  viq ;

ð30Þ

0 Pxi=q

ð31Þ

¼

x_ 0iq_



x0iq :

The quantity Pvi=q can be shown to be equal to zero. Since ri ¼ ri ðq; tÞ, we have vi ¼ r_ i ¼ riq q_ þ rit :

ð32Þ

From Eq. (32), form the partial derivative of vi with respect to q_ viq_ ¼ riq :

ð33Þ

Take the time derivatives of both sides of Eq. (33) to get v_ iq_ ¼ r_ iq ¼ viq :

ð34Þ

Substitute Eq. (34) into Eq. (30) to show that Pvi=q ¼ 0. When the angular velocity vector is the exact derivative 0 of another vector function of branch coordinates and time, the preceding argument will show that Pxi=q ¼ 0. This is always the case for problems formulated in one or two-dimensional space where angular velocities are simply the time derivatives of angular coordinates. In problems which must be formulated in three dimen0 sions, finite rotations cannot be represented as vectors and Pxi=q is generally non-zero. The virtual equation (29) is now simplified to h iT x0 0 0 0T 0 dW 5i ¼ dqT mi v_ Ti viq_ þ x_ 0T ½I x þ x ½I P : ð35Þ i i iq_ i i i=q The right hand side of Eq. (35) can be separated into terms containing branch accelerations €q and terms con_ This is necessary in order to place the coefficients of the terms €q into taining products of branch velocities q_ q. an augmented mass matrix c M i and the product terms into a generalized force QK i . Write the time-derivatives of the velocity vectors as q þ v_ pi ; v_ i ¼ viq_ € q þ x_ 0p x_ i ¼ x0iq_ € i :

ð36Þ ð37Þ

The terms superscripted p contain all products of generalized velocities. Then, exploit equations (36) and (37) to rewrite Eq. (35) into the final terminal equation for rigid bodies as follows h iT dW 5i ¼ dqT c q þ QK Mi€ ; ð38Þ i

144

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

where   0 0 c ½I x M i ¼ mi vTiq_ viq_ þ x0T _ _ iq i iq

ð39Þ

  x0 0 0 0T 0 _ pT _ 0pT QK i viq_ þ x i ½I i xiq_ þ xi ½I i Pi=q : i ¼ mi v

ð40Þ

and

The elements of the coefficient matrix c M i are the coefficient of €q appearing in the differential equation corresponding to the branch coordinate q. The forcing vector QK i shown in Eq. (40) is made up of all the terms in the equations of motion which do not contain second derivatives. This vector contains part of the contribution of the kinetic energy to the equations of motion. 4. Matrix formulation of constrained motion equations Together, the topological and terminal equations form a necessary and sufficient set of motion equations for multibody systems. However, an efficient approach to this problem consists in reducing the number of equations that need to be solved simultaneously by using branch transformation equations for a tree selection. The first step consists in defining the problem by creating a proper spanning tree and generating the branch transformation equations for all the elements in the cotree. These transformations will be used to replace the across-variables for these cotree elements as a function of branch coordinates q. Depending on the topology of the mechanical system and the specified tree, the branch coordinates may not be independent quantities. If the number of coordinates n is greater than the degrees of freedom f, then (c ¼ n  f ) constraint equations are required to express the dependency between coordinates. The constraint equations are obtained directly from the joints and motion drivers in the cotree, by projecting their circuit equations onto the joint reaction space [59]. Upon substitution of the branch transformation equations into these circuit equations, the constraint equations are obtained in the form of Eq. (16) which constitutes a set of c non-linear algebraic equations. Differentiating Eq. (16), using the chain rule of differentiation, yields the velocity constraint equations, Uq q_ ¼ Ut : Differentiating again this result with respect to time yields the c constraint acceleration equations h i _ q q_ þ 2Uqt q_ þ Utt  K: Uq € q ¼  ðUq qÞ

ð41Þ

ð42Þ

Alternatively, the velocity and acceleration equations (41) and (42) can easily be obtained by symbolic differentiation of the position-level constraint equation (16). Since the total kinetic energy of the system of rigid bodies is the sum of the individual kinetic energies, the algorithm must incorporate the sum of all bodies in the formulation using a global cutset equation. Then, the variational equations of motion for multibody systems may be obtained from the tree joint element which connects the whole system of bodies to the ground body through a path consisting entirely of branches. Substituting the general form of virtual work terminal equations for each element in the cutset equation for this tree body constraint or joint, one gets n o dqT c M€ q þ QK þ QC þ Qcons  QA ¼ 0: ð43Þ This variational equation of motion holds for arbitrary virtual displacement dq, so it is equivalent to the Lagrange multiplier form of constrained equations of motion. In general, the n branch coordinates are not independent, but are related by the c kinematic constraint equation (16). Thus, these constraint acceleration equation (42) must be appended to the set of dynamic equation (43), giving (n þ c) equations to solve for branch coordinates q and the c reaction loads k. Substituting Eq. (18) into Eq. (43) and combining with Eq. (42) finally yields the classical system differential–algebraic equations of motion in matrix form

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

"

c M Uq

) # ( € q Qtotal ; ¼ k K 0

UTq

145

ð44Þ

where c M is a symmetric augmented (n  n) mass matrix, Uq is a (c  n) constraint Jacobian matrix and the total generalized force Qtotal ¼ QA  QK  Qcons . Using Eq. (39) it is easy to form symbolically the coefficient matrix c M, element by element, from partial derivatives of the velocity vectors. Non-linearities have been preserved throughout the formulation as each element may contain branch coordinates and quantities which vary with time. Note that the global matrix coefficient is symmetric. This requires only that the inertia tensor be symmetric, which it is. The off-diagonal mass sub-matrices represent coupling effects between adjacent bodies. The fundamental form of these equations and physical properties of kinetic energy guarantee that a unique solution of the constrained equations of motion exists. 5. Automatic generation of motion equations Once the topology and parameters for a mechanical system has been defined, the translational and rotational graphs, with proper trees, can be generated automatically. At this point, the equations governing the motion of the mechanical system can be systematically formulated by following the five steps of the procedure outlined below. Let us consider the example of the spatial four-bar mechanism described earlier. Step 1: Across-variables function of branch coordinates In the spatial four-bar mechanism, a proper tree has been highlighted in bold in Fig. 2. From the linear graph, the following branch transformation equations are automatically generated, for cotree edge T, dhT ¼ dh31  a1

ð45Þ

for body edge e51 dh51 ¼ dh31 þ dh11  a1  q1

ð46Þ

and Eq. (6), for body edge e52 dh52 ¼ dh31 þ dh11 þ dh32 þ dh12  a1 þ a2  q1 þ q2 ;

ð47Þ

where the terminal equations of rigid-arm elements are equal to zero. Hence, the branch coordinates consists of a1 and a2, the coordinates associated with tree joint elements e31 and e32, respectively. Note that kinematic terminal equations, such as r31 ¼ 0. Using the branch transformation equations, all across-variables can be written as explicit functions of q ¼ ½a1 ; a2 T . Once the branch coordinates have been selected, the computerized algorithm registers the initial parameters, masses m1, m2 and inertia tensors ½I 01  ¼ diagð0; m31 ; m31 Þ, ½I 02  ¼ diagðm32 ; 0; m32 Þ. Now, from the circuit equations associated with edges e51 and e52, r51 ¼ r31 þ r11 ¼ R1 r011

ð48Þ

r52 ¼ r31 þ r11 þ r32 þ r12 ¼ R1 r011 þ R12 r012 ;

ð49Þ

and

where r01i is a (3  1) column matrix containing the constant component of r1i for (i ¼ 1; 2) in a local frame attached to the crank and R1 and R12 are (3  3) transformation matrices from the x01 y 01 z01 frame to the XYZ frame and from x02 y 02 z02 frame to the x01 y 01 z01 frame 2 3 2 3 cos a1  sin a1 0 1 0 0 6 7 6 7 ð50Þ cos a1 0 5; R12 ¼ 4 0 sin a2 cos a2 5 : R1 ¼ 4 sin a1 0 0 1 0  cos a2 sin a2 Hence, the positions of the centroid of the two bodies can be automatically generated

146

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

9 8 > = < cos a1 > r51  r1 ¼ sin a1 ; > > ; : 0

9 8 > = < cos a1  ðl2 =2Þ sin a1 sin a2 > r52  r2 ¼ sin a1 þ ðl2 =2Þ cos a1 sin a2 : > > ; : ðl2 =2Þ cos a2

ð51Þ

In a similar fashion the computerized version also needs to generate the velocities and angular velocities 9 9 8 8 > > = = < ½ sin a1  ðl2 =2Þ cos a1 sin a2 a_ 1  ½ðl2 =2Þ sin a1 cos a2 a_ 2 > <  sin a1 a_ 1 > ð52Þ fv1 g ¼ cos a1 a_ 1 ; fv2 g ¼ ½cos a1  ðl2 =2Þ sin a1 sin a2 a_ 1 þ ½ðl2 =2Þ cos a1 cos a2 a_ 2 ; > > > > ; ; : : 0 ðl2 =2Þ sin a2 a_ 2 9 8 9 8 a_ 2 > > > = = <0> < x01 ¼ ð53Þ 0 ; x02 ¼ RT12 ðx01 þ x012 Þ ¼  cos a2 a_ 1 : > > > ; ; : > : a_ 1 sin a2 a_ 1

Step 2: Velocity derivatives Once the translational and rotational velocities have been derived in terms of branch coordinates, then the algorithm calculates the non-zero partial derivatives, which are 9 8 > = <  sin a1  ðl2 =2Þ cos a1 sin a2 > v2q_ 1 ¼ ; ð54Þ cos a1  ðl2 =2Þ sin a1 sin a2 > > ; : 0 9 9 8 8 > > = = < ðl2 =2Þ sin a1 cos a2 > <  sin a1 > ; ð55Þ v1q_ 1 ¼ cos a1 ; v2q_ 2 ¼ ðl2 =2Þ cos a1 cos a2 > > > > ; ; : : 0 ðl2 =2Þ sin a2 9 8 8 9 8 9 0 > > > > = = = < <0> <1> ð56Þ x01q_ 1 ¼ 0 ; x02q_ 1 ¼  cos a2 ; x02q_ 2 ¼ 0 : > > > > ; ; ; : : > : > 1 0 sin a2 Also, the algorithm requires the acceleration terms for products in generalized velocities for each body in the system 9 9 8 8 2 0 > > > = = < <  cos a1 a_ 1 > p 0p 2 v_ 1 ¼  sin a1 a_ 1 ; x_ 2 ¼ cos a2 a_ 1 a_ 2 ; ð57Þ > > > > ; ; : : sin a2 a_ 1 a_ 2 0 9 8 2 2 > = < ½ cos a1 þ ðl2 =2Þ sin a1 sin a2 a_ 1  l2 cos a1 cos a2 a_ 1 a_ 2 þ ðl2 =2Þ sin a1 sin a2 a_ 2 > ð58Þ v_ p2 ¼ ½ sin a1  ðl2 =2Þ cos a1 sin a2 a_ 21  l2 sin a1 cos a2 a_ 1 a_ 2  ðl2 =2Þ cos a1 sin a2 a_ 22 : > > ; : 2 ðl2 =2Þ cos a2 a_ 2 To complete the virtual work terminal equations for the nb bodies in terms of branch coordinates, the algorithm needs to substitute these derivatives in Eqs. (39) and (40). Step 3: Generation of virtual cutset equations and mass matrix The only tree joint element that connects the system of bodies back to the ground, through a path consisting entirely of branches, is the revolute joint e31. Writing the cutset equation for this tree body element of virtual work dW T þ dW 61 þ dW 71 þ dW 31 þ dW 51 þ dW 72 þ dW 73 þ dW 52 þ dW 74 þ dW 75 þ dW 35 ¼ 0:

ð59Þ

Note that dW 31 and dW 34 are both zero by virtue of their terminal equation and that dW 72 and dW 73 cancel out, as well as dW 74 and dW 75 . Since the ideal joint e31 can do no work when connected to the ground body, the virtual work of the corresponding element dW 71 must also equal zero. Thus, the system virtual work, with their corresponding terminal equations, reduces to

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

dW T þ dW 61 þ dW 51 þ dW 52 þ dW 35 ¼ 0; n o c q þ QK þ QC þ Qcons þ dW T ¼ 0; q þ QK M 51 € dqT c 51 þ M 52 € 52 35 61

147

ð60Þ ð61Þ

which includes the contribution of all physical elements capable of doing work. Note that elements e35 and e61 represent the constraining link and the gravitational force, respectively. Once the branch derivatives are substituted into Eq. (61), one obtains an expression for the system virtual work entirely in terms of branch coordinates and Lagrange multipliers (da1 ; da2 ; k). Naturally, the complete augmented mass matrix and the kinetic forcing vector is obtained by summing Eqs. (39) and (40), respectively, over both bodies c M¼c M 51 þ c M 52 and K K K Q ¼ Q51 þ Q52 . At this point, the algorithm can generate the mass matrix, element by element, using Eq. (39), and summing for both bodies " #

2  3l2 þ 4 sin2 a2 4m1 b M 11 ¼ þ m2 1 þ ; ð62Þ 3 12 b 12 ¼ m2 ðl2 =2Þ cos a2 ¼ M 21 ; M   2 b 22 ¼ m2 1 þ l2 cos2 a2 þ sin2 a2 : M 3 4

ð63Þ ð64Þ

Step 4: Computation of the forcing vectors Now that the inertial part of the problem is complete and, after eliminating all the secondary variables, the T automated algorithm may generate symbolically the two generalized kinetic forces QK ¼ ½QK 1 ; QK 2  ,   0pT 0 x0 0 0T 0 _ QK 1 ¼  m2 v_ pT v þ x ½I x þ x ½I P ð65Þ 2 2q_ 1 2 2=1 ; 2 2q_ 1 2 2 where 0

Px2=1 ¼ x_ 02q_ 1  x02q1 ¼

8 > < > :

0

9 > =

sin a2 a_ 1 a_ 2 ; > ; cos a2 a_ 1 a_ 2

ð66Þ

hence QK 1 ¼ m2 sin a2

 2    l2 3l2 þ 4 cos a2 a_ 1 a_ 2 a_ 2  6 2

ð67Þ

and   x0 0T 0 v þ x ½I P QK 2 ¼  m2 v_ pT 2 q _ 2 2 2=2 ; 2 2 where 0

ð68Þ

9 0 > = ¼ sin a2 a_ 1 ; > > ; : cos a2 a_ 1

ð69Þ

 3l22 þ 4 sin a2 cos a2 a_ 21 : 12

ð70Þ

Px2=2 ¼ x_ 02q_ 2  x02q2

8 > <

hence QK 2 ¼ m2



The potential energy and generalized applied force T are necessary to complete the formulation. Gravitational forces in this example are conservative and can be generated from the potential energy V, V ¼ m2 g

l2 ð1  cos a2 Þ: 2

ð71Þ

148

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

To satisfy Eq. (19), one also needs the partial derivative of V with respect to branch coordinate a2. l2 V q2 ¼ m2 g sin a2 : ð72Þ 2 Step 5: Constraint circuit equations If a multibody system has no virtual constraining elements dW C in the dynamical equations of motion generated in the preceding step, then this step 5 is not required to solve for the dynamic response of the system. However, for the spatial four-bar mechanism, there exists a constraining element e35 which is a rigid-couplerlink. For a mechanical system with closed kinematic chains, one can generally define branch coordinates that are fewer in number than a set of joint coordinates for the same system. Nevertheless, the n branch coordinates will still be dependent variables, related by c non-linear algebraic equations of constraint. These are obtained by first generating the circuit equation (4) for active and passive kinematic constraints. By projecting these circuit equations onto the corresponding reaction space [59], the across-variables for cotree constraint elements are automatically eliminated. The resulting constraint equations are in terms of the branch coordinates. Note that one algebraic equation is obtained for each reaction force or torque appearing in the dynamic equations. Since the spatial four-bar mechanism has one degree of freedom and two branch coordinates, one constraint equation is required. This equation is obtained from the circuit cotree link e35 equation, r35  r33 þ r13  r12  r32  r11  r31 þ r21  r34 ¼ 0:

ð73Þ 35

After eliminating elements equal to zero and projecting the resulting equation onto the unit vector e , ðr35 Þ  e35 þ ðr13  r12  r11 þ r21 Þ  e35 ¼ 0:

ð74Þ

Since the magnitude of the constraining link jr35 j ¼ 1, with

½r11 þ ðr12  r13 Þ  r21   ½r11 þ ðr12  r13 Þ  r21  ¼ 1

ð75Þ

9 8 9 9 8 8 0 > > > = = > <0> = < < cos a1 >  2 r11 þ ðr12  r13 Þ  r21 ¼ sin a1 þ R1 l2 sin a2 > > > > ; ; > : > ; : : 0 0 l2 cos a2

ð76Þ

and the constraint equation is U  l22  4l2 cos a1 sin a2  4 sin a1 þ 4 ¼ 0: Thus dq must satisfy Uq dq ¼ ½4l2 sin a1 sin a2  4 cos a1 ; 4l2 cos a1 cos a2 dq:

ð77Þ ð78Þ

By the Lagrange multiplier theorem, there exists a Lagrange multiplier k such that augmented variational equation of motion holds for all dq. The acceleration constraint equation that completes the system equations of motion is Uq € q ¼ ð4l2 cos a1 sin a2 þ 4 sin a1 Þa_ 21  8l2 sin a1 cos a2 a_ 1 a_ 2  4l2 cos a1 sin a2 a_ 22 ¼ K: Thus the differential–algebraic equations of motion for the spatial four-bar mechanism are 2 3 h i 8 9 ð3l22 þ4Þ sin2 a2 4m1 þ m 1 þ ðl =2Þ cos a 4l sin a sin a  4 cos a m 2 2 2 2 2 1 2 1 >€ 12 = < a1 > 6 3 7 h i 6 7 2 l € 2 a 6 7 1 2 2 2 m2 ðl2 =2Þ cos a2 m2 3 þ 4 cos a2 þ sin a2 4l2 cos a1 cos a2 4 5> ; : > k 4l2 sin a1 sin a2  4 cos a1 4l2 cos a1 cos a2 1 9 8 K1 > = < T þQ > ¼ QK 2  V q2 ; > > ; : K which can be solved by numerical integration to obtain the time varying response of the system.

ð79Þ

ð80Þ

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

149

In this second example, the equations of motion will be derived for a planar articulated semi-trailer system which consists of two rigid bodies. As shown in Fig. 4, the car is a conventional vehicle of mass mc, with three degrees of freedom and supports the front end of a trailer of mass mt on a platform or ‘‘fifth-wheel’’ which has freedom to rotate in yaw and the rear end of the trailer is mounted on a single axle. As two units are involved, four degrees of freedom are required to completely describe the system. The forces and moments acting on the system are shown in Fig. 5. Since this general formulation exploits kinetic energy, one may choose to represent the inertial-space branch coordinates with respect to a moving reference frame attached to the car. Even though, absolute across-variables must be used, velocities may be expressed in a convenient non-inertial coordinate system. A right hand set of moving unit vectors ex ; ey ; ez are defined so that; ex is aligned with the forward direction of the car, ey is perpendicular to ex and directed to the right and ez is directed vertically downward. All vectors and tensors defined in this problem are expressed with respect to these unit vectors. The tree for the linear graph representation of the semi-trailer is shown in bold in Fig. 4, where the rigidarm elements to the axles and all the driving forces, traced on Fig. 4, have not been drawn on the graph for clarity. The two rigid bodies are represented by edges e51 and e52 while rigid-arm elements are represented by edges e11 and e12. Note that the non-inertial reference frame is moving and rotating with the vehicle. For this open-loop system, both rigid-arm elements can be selected into the spanning tree. From this choice of tree T selection, one obtains the set of branch coordinates q ¼ ½x; y; w; wt  . Since these coordinates are independent, the dynamic equations of motion will consists of four ordinary differential equations and there are no constraint equations or corresponding Lagrange multiplier, hence, step 5 is no longer needed for this formulation. Once the branch coordinates have been selected, the algorithm begins with the identification of system parameters for the mathematical model, masses m1 ¼ mc , m2 ¼ mt and inertia tensors ½I 01  ¼ diagð0; 0; I c Þ, ½I 02  ¼ diagð0; 0; I t Þ. From the circuit equation associated with edge e52, r52 ¼ r51 þ r11 þ r31  r12 ¼ r51 þ r11  R1 r012

ð81Þ

δf a

b

d1

Y

X

e11

ψ e51

e31 δr ψt

e52

e12

Fig. 4. Graph of the schematic top view of a car-trailer model.

d2

c

150

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

PACx

δf Xf Yf Tf XAC TAC

YAC

δr

Xr

Yr

Tr

PATx

XAT

YAT Yt

TAT Tt

Xt

Fig. 5. Forces and torques acting on a car-trailer model (X = longitudinal force, Y = lateral force, T = torque,  ) aerodynamic reference point) (subscripts; f = front tire, r = rear tire, t = trailer tire, A = aerodynamic, C = car, T = trailer).

and in terms of branch coordinates, 9 8 9 8 9 8 9 8 > = > = > = > = < d 1 > < d 2 cos wt > < x  d 1  d 2 cos wt >  d 2 sin wt ¼ : r2 ¼ y þ 0 y  d 2 sin wt > > > > ; > ; > ; > ; : > : : : 0 0 0 0

ð82Þ

Care must be taken when differentiating the body positions since the reference frame is non-inertial. When the axes of reference are translating and rotating with the mechanical system, a moving frame operator for all vectors must be applied such as   d½  d½  þ x05  ½ : ð83Þ ¼ dt dt x0 y 0 z0 Hence, the velocity vector of the trailer, with respect to moving axis, becomes 9 9 8 8 _ _ > = = > < wðy  d 2 sin wt Þ > < x_ þ d 2 wt sin wt > _  d 1  d 2 cos w Þ ; v2 ¼ y_  d 2 w_ t cos wt þ wðx t > > > ; ; > : : 0 0

ð84Þ

where naturally, the second term of Eq. (84) takes into account the rotational effect of the reference frame. Note that x ¼ y ¼ 0 since the moving frame is superimposed over the car’s reference frame. A consequence of this is that the properties of the across-variables, measured across an element, from the position to the velocity level, no longer holds for problems with rotating reference frames. Nonetheless, the linear graph generated by the algorithm is still valid. Once the linear graph has been constructed, relatively simple expressions are generated, in Appendix, by the symbolic program. The equations of motion for the car-trailer system can finally be automatically generated

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

151

by systematically applying the six steps formulation provided in Appendix, leading to the four differential equations of motion 2 38 9 mc þ mt 0 mt d 2 sin wt mt d 2 sin wt €x > > > > = < €y > 6 7> 0 m þ m m ðd þ d cos w Þ m d cos w c t t 1 2 t 2 6 7 t t 6 7 €> 4 mt d 2 sin wt mt d 1  mt d 2 cos wt I c þ I t þ mt d 21 þ mt ðd 22 þ 2d 1 d 2 cos wt Þ I t þ mt d 22 þ mt d 1 d 2 cos wt 5> w > > > ; :€ > 2 2 mt d 2 sin wt mt d 2 cos wt I t þ mt ðd 2 þ d 1 d 2 cos wt Þ I t þ mt d 2 wt 9 8 K1 Q þ QA1 > > > > > = < QK 2 þ QA2 > : ¼ > > QK 3 þ QA3 > > > > ; : K4 Q þ QA4 ð85Þ Numerical techniques can then be easily exploited to investigate the behaviour of the solution of these nonlinear differential equations. The two examples presented in this section illustrate a generally applicable method for deriving differential equations of motion for multibody systems, providing that independent branch coordinates may be found. In the semi-trailer example, it was possible to define independent branch coordinates that lead to non-linear ordinary differential equations of motions. In the spatial four-bar mechanism, a constraint equation function of branch coordinates was retained and incorporated in the differential-algebraic equations of motion, using the Lagrange multiplier theorem. Even in these relatively elementary examples, extremely non-linear differential equations are obtained symbolically, precluding practical generation of analytical solutions. Numerical integration techniques are therefore, generally, needed to solve these equations of motion. By exploiting the symbolic C++ language, steps 2–5 of this formulation have been implemented into a computer program called VARNET (for VARiational NETwork). Given only a spanning tree with the velocity expressions of bodies in the system, this program automatically generates the symbolic equations of motion. Since the selection of a proper tree only requires an elementary knowledge of graph theory, the objective consists in choosing an optimal joint tree to keep the number of branch coordinates and constraint equations to a minimum. Unfortunately, this symbolic software is not fully automated, since it is used mainly for academic purposes, because the branch coordinates are picked and optimized by the user for each particular multibody system. 6. Conclusion By combining the mechanical system topology with the Lagrange variational constitutive equations, a new systematic formulation procedure has been introduced and used to describe the time-varying configuration of spatial multibody systems. This method assembles automatically the governing equations of motion in a symmetrical format where the structure and organization of the mass matrix parallels that of structural finite element mass and stiffness matrices, which are also derived using variational methods. The systematic nature of the variational approach, for both dynamics and structural mechanics, follows from the additive property of energies in the variational formulation. The classical approach to the analysis of spatial dynamic systems consists in deriving the equations of motion in all system generalized coordinates. As seen in the past by other methods, these equations are so complicated that a symbolic representation is always difficult and, hence, solutions must be constructed numerically. Given this undesirable situation, it is important to consider methods for the formulation of motion equations that have attractive computational properties such as the variational graph-theoretic method presented in this work. For open-loop systems, a joint tree results in independent branch coordinates and a set of reduced ordinary differential equations can be generated. However, if the graph of a multibody system has closed loops, a tree structure is formed by mathematically cutting the constraining elements or joints yielding the constraint circuit equations. The kinematics formulation of step 5 may then be developed for the resulting spanning tree, so it

152

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

neglects momentarily the effect of kinematic constraints between other bodies. While the variational equation of motion is still valid, it holds only for branch coordinate variations dq that are consistent with the constraint. It is, therefore, necessary to introduce Lagrange multipliers associated with the constraint equations. In this new formulation, graph-theoretic methods are used to obtain the velocities of the bodies in the multibody system. Then, the algorithm generates automatically the governing equations of motion by following steps 2–5. The main advantage of this method resides in the fact that vector multiplications, particularly for spatial systems, have been reduced considerably since the variational Lagrange formulation begins its substitution process at the velocity level. A consequence of this approach implies that it can also analyze systems with holonomic or non-holonomic constraints. Furthermore, the method can choose branch coordinates that can be defined with respect to an inertial or non-inertial reference frame, thereby reducing the computational time for obtaining the motion equations expressed in a moving reference frame. Finally, another advantage of this symbolic approach consists in eliminating mechanical singularities that often arise with a purely numerical approach. The systematic nature of this method will be extended, in the future, to incorporate the control analysis of automated guided vehicle systems and mechatronic systems. Prior to the advent of symbolic computation methods, evaluation of these types of systems was oppressive. This symbolic computation program, however, carries out the algebraic manipulations required to form executable computer code to analyse these systems and makes evaluation of such systems practical and computationally attractive. Acknowledgement Financial support of this research by the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. Appendix In this appendix, the analytical expressions for quantities necessary for the variational Lagrangian derivation of the equations of motion are derived. This section contains all the vectors and coefficients generated by the symbolic algorithm in four sequential steps to obtain the differential equation (85) of motion for the articulated semi-trailer system. Step 1: Across-variables function of branch coordinates 8 9 8 9
ð86Þ

The velocities of the car and trailer masses can be expressed as a vector in the x  y  z moving coordinate system 8 9 8 9 x_ þ d 2 ðw_ þ w_ t Þ sin wt < x_ = < = v1 ¼ y_ ; v2 ¼ y_  d 1 w_  d 2 ðw_ þ w_ t Þ cos wt ; : ; : ; 0 8 9 8 9 0 ð87Þ <0= < 0 = x01 ¼ 0 ; x02 ¼ : 0 :_; :_ ; w w þ w_ t Step 2: Velocity derivatives Only non-zero derivatives are shown in this section 8 9 8 9 > > = = <1> <0> v1q_ 1 ¼ v2q_ 1 ¼ 0 ; v1q_ 2 ¼ v2q_ 2 ¼ 1 ; > > ; ; : > : > 0 0

ð88Þ

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

8 <

9 8 9 d 2 sin wt = < d 2 sin wt = v2q_ 3 ¼ d 1  d 2 cos wt ; v2q_ 4 ¼ d 2 cos wt ; : ; : ; 0 0 8 9 <0= 0 0 0 x1q_ 3 ¼ x2q_ 3 ¼ x2q_ 4 ¼ 0 ; : ; 1 9 9 8 8 _2 _ _ _ 2 _ > > > = < y_ w = < y_ w þ d 1 w þ d 2 ðw þ wt Þ cos wt > ; v_ p2 ¼ : v_ p1 ¼ x_ w_ _ w_ þ d 2 ðw_ þ w_ t Þ2 sin wt x > > > > ; ; : : 0 0

153

ð89Þ

ð90Þ

ð91Þ

Step 3: Generation of virtual cutset equations and mass matrix The elements of the mass matrix, defined by summing equation (39) for both bodies, can now be found to be M 11 M 12 M 13 M 23 M 24 M 33 M 34 M 44

¼ M 22 ¼ mc þ mt ; ¼ M 21 ¼ 0; ¼ M 31 ¼ M 14 ¼ M 41 ¼ mt d 2 sin wt ; ¼ M 32 ¼ mt ðd 1 þ d 2 cos wt Þ; ¼ M 42 ¼ mt d 2 cos wt ; ¼ I c þ I t þ mt ðd 21 þ d 22 þ 2d 1 d 2 cos wt Þ; ¼ M 43 ¼ I t þ mt ðd 22 þ d 1 d 2 cos wt Þ; ¼ I t þ mt d 22 :

ð92Þ ð93Þ ð94Þ ð95Þ ð96Þ ð97Þ ð98Þ ð99Þ

Step 4: Computation of the forcing vectors The kinetic generalized forces can now be found for each generalized coordinates by summing Eq. (40) for both bodies QK 1 ¼ ðmc þ mt Þ_y w_  mt d 1 w_ 2  mt d 2 ðw_ þ w_ t Þ2 cos wt ; ð100Þ 2 K2 _ _ _ Q ¼ ðmc þ mt Þ_xw  mt d 2 ðw þ wt Þ sin w ; ð101Þ t

K3

Q

K4

Q

_ 2 sin w þ mt x_ wðd _ 1 þ d 2 cos w Þ þ mt d 1 d 2 sin w ðw_ 2 þ 2w_ w_ t Þ; ¼ mt y_ wd t t t t 2 _ _ _ ¼ mt y_ wd 2 sin w þ mt x_ wd 2 cos w  mt d 1 d 2 w sin w : t

t

t

ð102Þ ð103Þ

Aj

Finally, the generalized applied forces Q corresponding to the virtual displacement dqj , with all other coordinates being kept constant can be calculated, QA1 ¼ X f cos df þ X r cos dr þ X t cos wt  Y f sin df  Y r sin dr  Y t sin wt þ X Ac þ X At cos wt  Y At sin wt ; A2 Q ¼ X f sin df þ X r sin dr þ X t sin wt þ Y f cos df þ Y r cos dr þ Y t cos wt þ Y Ac þ X At sin wt  Y At cos wt ; A3 Q ¼ aðX f sin df þ Y f cos df Þ  bðX r sin dr þ Y r cos dr Þ  Y Ac P Acx þ ðd 2 þ P Atx Þ sin wt ðX At cos wt  Y At sin wt Þ  fd 1 þ ðd 2 þ P Atx Þ cos wt gðX At sin wt þ Y At cos wt Þ þ ðd 2 þ cÞ sin wt ðX t cos wt  Y t sin wt Þ  fd 1 þ ðd 2 þ cÞ cos wt gðX t sin wt þ Y t cos wt Þ þ T f þ T Ac þ T r þ T At þ T t

ð104Þ ð105Þ

ð106Þ

by simplifying, QA3 ¼ aX f sin df  bX r sin dr  d 1 X t sin wt þ aY f cos df  bY r cos dr  ðd 2 þ c þ d 1 cos wt ÞY t  Y Ac P Acx  d 1 X At sin wt  fðd 2 þ P Atx Þ þ d 1 cos wt gY At þ T f þ T Ac þ T r þ T At þ T t

ð107Þ

154

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

and QA4 ¼ ðd 2 þ P Atx Þ sin wt ðX At cos wt  Y At sin wt Þ  ðd 2 þ P Atx Þ cos wt ðX At sin wt  Y At cos wt Þ þ ðd 2 þ cÞ sin wt ðX t cos wt  Y t sin wt Þ  ðd 2 þ cÞ cos wt ðX t sin wt þ Y t cos wt Þ þ T At þ T t  T YB ð108Þ by simplifying, QA4 ¼ ðd 2 þ cÞY t  ðd 2 þ P Atx ÞY At þ T At þ T t  T YB :

ð109Þ

References [1] B. Paul, Analytical dynamics of mechanisms – a computer oriented overview, Mechanism and Machine Theory 10 (1975) 481. [2] J. Unda, J. Garcia de Jalon, F. Losantos, R. Enparantza, A comparative study on some different formulations of the dynamic equations of constrained mechanical systems, Journal of Mechanisms, Transmissions and Automation in Design 109 (1987) 466. [3] R.E. Roberson, Computer-oriented dynamic modeling of spacecraft: historical evolution of Eulerian multibody formalisms since 1750, in: 28th International Astronautical Congress, Prag, IAF Paper 77-A11, 1977. [4] R. Schwertassek, R.E. Roberson, A perspective on computer oriented multibody dynamic formalisms and their implementations, in: Bianchi, Schiehlen (Eds.), Proceedings IUTAM/IFToMM Symposium on Dynamics of Multibody Systems, Springer-Verlag, Udine, Italy, 1986, p. 261. [5] E.J. Haug, Elements and methods of computational dynamics, in: E.J. Haug (Ed.), Computer Aided Analysis and Optimization of Mechanical System Dynamics, Springer-Verlag, 1984, p. 3. [6] B. Paul, Computer oriented analytical dynamics of machinery, in: E.J. Haug (Ed.), Computer Aided Analysis and Optimization of Mechanical System Dynamics, Springer-Verlag, 1984, p. 41. [7] J. Wittenburg, Dynamics of multibody systems – a brief review, in: 39th International Astronautical Federation Conference, Bangalore, India, 1988, p. 89. [8] G.C. Andrews, A brief survey of self-formulating simulation programs for multi-body dynamic systems, in: IASTED International Symposium on Simulation and Modelling, Orlando, FL, 1980. [9] D. Benson, Performance analysis of four computational dynamics algorithms, in: Proceedings of the 2nd International Computer Engineering Conference, San Diego, CA, 1982. [10] W. Schiehlen, Multibody Systems Handbook, Springer-Verlag, 1990. [11] D.A. Smith, Reaction forces and impact in generalized two-dimensional mechanical dynamic systems, Ph.D. thesis, University of Michigan, Ann Arbor, MI, 1971. [12] N. Orlandea, Node-analogous, sparsity-oriented methods for simulation of mechanical dynamic systems, Ph.D. thesis, University of Michigan, Ann Arbor, MI, 1973. [13] N. Orlandea, M.A. Chace, D.A. Calahan, A sparsity-oriented approach to the dynamic analysis and design of mechanical systems – Part 1, Journal of Engineering for Industry 99 (1977) 773. [14] N. Orlandea, D.A. Calahan, M.A. Chace, A sparsity-oriented approach to the dynamic analysis and design of mechanical systems – Part 2, Journal of Engineering for Industry 99 (1977) 780. [15] M.A. Chace, Y.O. Bayazitoglu, Development and application of a generalized D’Alembert force for multi-freedom mechanical systems, Journal of Engineering for Industry 93 (1971) 317. [16] M.A. Chace, Methods and experience in computer aided design of large displacement mechanical systems, in: E.J. Haug (Ed.), Computer Aided Analysis and Optimization of Mechanical System Dynamics, Springer-Verlag, Berlin, 1984, p. 233. [17] M.A. Chace, Using DRAM and ADAMS programs to simulate machinery, Vehicles, Agricultural Engineering, November, p. 18 and December, p. 16, 1978. [18] ADAMS Users Manual, Mechanical Dynamics Incorporated, 2301 Commonwealth Blvd., Ann Arbor, MI 48105, November 1989. [19] P.N. Sheth, A digital computer based simulation procedure for multiple degree of freedom mechanical systems with geometric constraints, Ph.D. thesis, University of Wisconsin, Madison, WI, 1972. [20] P.N. Sheth, J.J. Uicker, IMP (integrated mechanisms program), a computer-aided design analysis system for mechanisms and linkage, Journal of Engineering for Industry 94 (1972) 454. [21] P.N. Sheth, J.J. Uicker, A generalized symbolic notation for mechanisms, Journal of Engineering for Industry 93 (1971) 102. [22] P.N. Sheth, T.M. Hodges, J.J. Uicker, Matrix analysis method for direct and multiple contact multibody systems, Journal of Mechanisms, Transmissions and Automation in Design 112 (2) (1990) 145. [23] B. Paul, Kinematics and Dynamics of Planar Machinery, Prentice-Hall, 1979. [24] B. Paul, D. Krajcinovic, Computer analysis of machines with planar motion – 1. Kinematics, Journal of Applied Mechanics 37 (1970) 697. [25] B. Paul, D. Krajcinovic, Computer analysis of machines with planar motion – 2. Dynamics, Journal of Applied Mechanics 37 (1970) 703. [26] B. Paul, Dynamic analysis of machinery via program DYMAC, SAE Paper 770049, Society of Automotive for Engineers, Warrendale, PA, 1977.

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

155

[27] R.A. Wehage, E.J. Haug, Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic systems, Journal of Mechanical Design 104 (1982) 247. [28] M.W. Dubetz, J.G. Kuhl, E.J. Haug, Network implementation of real-time dynamic simulation with interactive animated graphics, in: 1988 ASME Design Automation Conference, Kissimmee, FL, 1988, p. 519. [29] DADS Users Manual, Computer Aided Design Software Incorporated, P.O. Box 203, Oakdale, Iowa 52319, 1990. [30] T.W. Park, E.J. Haug, Hybrid numerical integration method for machine dynamic simulation, Journal of Mechanisms, Transmissions and Automation in Design 108 (2) (1986) 211. [31] W.O. Schiehlen, E.J. Kreuzer, Symbolic computerized derivation of equations of motion, in: K. Magnus (Ed.), Dynamics of Multibody Systems, Springer-Verlag, 1978. [32] W.O. Schiehlen, Computer generation of equations of motion, in: E.J. Haug (Ed.), Computer Aided Analysis and Optimization of Mechanical System Dynamics, Springer-Verlag, Berlin, 1984, p. 83. [33] W.O. Schiehlen, Dynamics of complex multibody systems, Solid Mechanics Archives 9 (2) (1984) 159. [34] W.O. Schiehlen, Computational dynamics of multibody systems, in: Proceedings of the 1990 ASME International Computers in Engineering Conference and Exposition, Boston, MA, 1990, p. 519. [35] E.J. Kreuzer, W.O. Schiehlen, Equations of motion and equations of stress for robots and manipulators; an application of the NEWEUL formalism, in: 5th CISM/IFToMM Symposium, Udine, Italy, 1985, p. 79. [36] E.J. Kreuzer, W.O. Schiehlen, Generation of symbolic equations of motion for complex spacecraft using formalism NEWEUL, in: Proceedings of the AAS/AIAA Astrodynamics Conference, Lake Placid, NY, 1984, p. 21. [37] E.J. Kreuzer, W.O. Schiehlen, Computerized generation of symbolic equations of motion for spacecraft, Journal of Guidance, Control and Dynamics 8 (2) (1985) 284. [38] W.O. Schiehlen, Modeling of complex vehicle systems, in: 8th IAVSD-IUTAM Symposium on Dynamics of Vehicles on Roads and Tracks, vol. 12, 1983, p. 12. [39] J. Unda, A. Avello, J.M. Jimenez, J. Garcia de Jalon, COMPAMM – a program for the dynamic analysis of multi-rigid-body systems, in: Proceedings of the 4th International Conference on Engineering Software IV, London, England, 1985, p. 16-3. [40] J. Garcia de Jalon, J. Unda, A. Avello, J.M. Jimenez, Dynamic analysis of three-dimensional mechanisms in ‘‘natural’’ coordinates, Journal of Mechanisms, Transmissions and Automation in Design 109 (1987) 460. [41] J. Garcia de Jalon, J. Unda, A. Avello, Natural coordinates for the computer analysis of multibody systems, Journal of Computer Methods in Applied Mechanics and Engineering 56 (3) (1986) 309. [42] G.C. Andrews, H.K. Kesavan, The vector-network model: a new approach to vector dynamics, Mechanism and Machine Theory 10 (1975) 57. [43] G.C. Andrews, The vector-network model: a topological approach to mechanics, Ph.D. thesis, University of Waterloo, Waterloo, Ont., 1971. [44] G.C. Andrews, A general re-statement of the laws of dynamics based on graph theory, in: Problem Analysis in Science and Engineering, Academic Press, 1977. [45] G.C. Andrews, H.K. Kesavan, Simulation of multibody systems using the vector-network model, in: Dynamics of Multibody Systems, IUTAM Symposium, Munich, Germany, 1977. [46] G.C. Andrews, Dynamics using vector-network techniques, Course Lecture-Notes, Course No. ME-730, Department of Mechanical Engineering, University of Waterloo, Waterloo, Ont., 1975. [47] M.J. Richard, R.J. Anderson, Dynamic simulation of three-dimensional rigid body systems using vector-network techniques, International Journal in Mathematics and Computer Simulation 26 (1984) 289. [48] M.J. Richard, Dynamic simulation of constrained three-dimensional multi-body systems using vector-network techniques, Ph.D. thesis, Queen’s University, Kingston, Ont., 1985. [49] G.C. Andrews, M.J. Richard, R.J. Anderson, A general vector-network formulation for dynamic systems with kinematic constraints, Journal of Mechanism and Machine Theory 23 (3) (1988) 243. [50] M.J. Richard, R.J. Anderson, G.C. Andrews, Generalized vector-network formulation for the dynamic simulation of multibody systems, Journal of Dynamic Systems, Measurement and Control 108 (1986) 322. [51] M.J. Richard, R.J. Anderson, G.C. Andrews, The vector-network method for the modelling of mechanical systems, International Journal of Mathematics and Computers in Simulation 31 (1990) 565. [52] M.J. Richard, Dynamic simulation of multi-body mechanical systems using the vector-network model, Transaction of the Canadian Society of Mechanical Engineering 12 (1) (1988) 21. [53] M. Leger, J. McPhee, Selection of modelling coordinates for forward dynamic multibody simulations, Multibody System Dynamics, submitted for publication. [54] C. Schmitke, J. McPhee, Using linear graph theory and the principle of orthogonality to model multibody, multi-domain systems, Advanced Engineering Informatics, submitted for publication. [55] M. Fayet, F. Pfister, Analysis of multibody systems with indirect coordinates and global inertia tensors, European Journal of Mechanics 13 (3) (1994) 431. [56] R.L. Huston, Y.S. Liu, C. Liu, Use of absolute coordinates in computational multibody dynamics, Computers and Structures 52 (1) (1994) 17. [57] S.K. Agrawal, M.P. Shelly, Dynamic simulation and choice of generalized coordinates for a free-floating closed-chain planar manipulator, Mechanism and Machine Theory 28 (1993) 615. [58] J.J. McPhee, Automatic generation of motion equations for planar mechanical systems using the new set of branch coordinates, Mechanism and Machine Theory 33 (1998) 805.

156

M.J. Richard et al. / Applied Mathematics and Computation 192 (2007) 135–156

[59] P. Shi, J.J. McPhee, On the use of virtual work in a graph-theoretic formulation for multibody dynamics, in: ASME Design Engineering Technical Conference, Sacramento, CA DETC97/VIB-4199, 1997. [60] M. Behzad, G. Chartrand, Introduction to the Theory of Graphs, Allyn and Bacon, 1971. [61] N. Christofides, Graph Theory, An Algorithm Approach, Academic Press, New York, 1975. [62] R.E. Roberson, The path matrix of a graph, its construction and its use in evaluating certain products, Journal of Computer Methods in Applied Mechanics and Engineering 42 (1984) 47. [63] K. Arczewski, Application of graph theory to the mathematical modelling of a class of rigid body systems, Journal of the Franklin Institute 327 (2) (1990) 209. [64] G. Baciu, J.C.K. Chou, H.K. Kesavan, Constrained multibody systems: graph-theoretic Newton–Euler formulation, IEEE Transactions on Systems, Man and Cybernetics 20 (5) (1990) 1025. [65] K. Arczewski, Application of graph theory to the determination of kinetic energy of rigid body systems, Journal of the Franklin Institute 324 (3) (1987) 351. [66] J.C. Chou, H.K. Kesavan, K. Singhal, Dynamics of 3-D isolated rigid-body systems: graph-theoretic models, Journal of Mechanism and Machine Theory 21 (3) (1986) 261. [67] A.M. Bos, Modelling multibody systems in terms of multibond graphs with application to a motorcycle, Dissertation Twente University, Enschede, The Netherlands, 1986. [68] Y. Hu, Applications of bond graphs and vector bond graphs to rigid body dynamics, Journal of China Textile University, English Edition 5 (4) (1988) 67. [69] D.L. Margolis, Bond graphs for automated simulation and control of nonlinear vehicle systems, in: Future Transportation Technology Conference and Exposition, Seattle WA, SAE paper no. 871558, 1987, 8 p. [70] P. Shi, J.J. McPhee, Dynamics of flexible multibody systems using virtual work and linear graph theory, Multibody System Dynamics 4 (2000) 355.