Graph-Theoretic Approach for the Dynamic Simulation of Flexible Multibody Systems

Abstract

This paper provides a general description of a variational graph-theoretic formulation for simulation of flexible multibody systems (FMSs) which includes a brief review of linear graph principles required to formulate this algorithm. The system is represented by a linear graph, in which nodes represent reference frames on flexible bodies, and edges represent components that connect these frames. The method is based on a simplistic topological approach which casts the dynamic equations of motion into a symmetrical format. To generate the equations of motion with elastic deformations, the flexible bodies are discretized using two types of finite elements. The first is a 2 node 3D beam element based on Mindlin kinematics with quadratic rotation. This element is used to discretize unidirectional bodies such as links of flexible systems. The second consists of a triangular thin shell element based on the discrete Kirchhoff criterion and can be used to discretize bidirectional bodies such as high-speed lightweight manipulators, large high precision deployable space structures, and micro/nano-electromechanical systems (MEMSs). Two flexible systems are analyzed to illustrate the performance of this new variational graph-theoretic formulation and its ability to generate directly a set of motion equations for FMS without additional user input.

1. Introduction

Stringent tolerances on mechanical components have created increasingly severe 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 predicting and improving the efficiency of an FMS design. Since the equations governing the motion of FMS are highly nonlinear and dynamically coupled, one must exploit some kind of linear graph principles [1–4] to properly define the interconnection between the bodies. By combining linear graph theory [5–8] with the principle of virtual [9–11] work and finite elements, a dynamic formulation is obtained that extends graph-theoretic (GT) modelling methods to the analysis of 3D beams and shell surfaces of FMS. The widespread interest in flexible multibody systems is evidenced by the existence of a large number of algorithms [12]. It has been shown [13–15] that for multibody systems, a graph-theoretic formulation can generate a minimal set of equations. GT-based approaches explicitly separate the linear topological equations for the entire system from the nonlinear constitutive equations for individual components, resulting in very modular and efficient algorithms.

This work differs from past publications [10, 11] since it exploits a complete new graph-theoretic algorithm that automatically generates the equations of motion for FMS, given only a topological description of the system as input. Recognizing that the choice of coordinates employed has a direct effect on the governing differential equations and therefore on the computational time to solve these nonlinear equations, this approach will automatically produce a minimal set of differential equations equal in number to the degrees of freedom. By combining linear graph theory with the principle of virtual work, a new variational Lagrangian formulation is obtained that extends graph-theoretic modelling methods to the analysis of FMS. It shows also how a Lagrangian multiplier technique can be incorporated into a flexible finite element algorithm. In order to apply formal variational graph-theoretic formulations to a multibody system, it is necessary to create a linear graph that encapsulates all the constraint reaction forces in the system. Also, past models often had some restrictions when mixing rigid and flexible bodies, especially for closed loop systems. However, since this graph-theoretic approach is based on rigid body dynamics, both rigid and flexible bodies can be present in the same model for any configuration.

To be applicable to a wide range of spatial mechanical systems containing both open and closed kinematic chains, a flexible 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. It has also been shown [16] that a proper “tree” with a set of coordinates called “branch coordinates” encompasses both sets of Cartesian and joint coordinates as special cases. Also, 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 flexible dynamic equations for any choice of spanning tree, it is necessary to introduce a dependent virtual work element.

2. System Representation by GT Modelling

Many researchers have studied the theory of graphs [17–23] and bond-graphs [24–26] 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 rigid or flexible 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 branch. 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”.

It is possible to determine the number of branches and chords in a given graph $Graph (ν, ∊)$ by applying basic theorems of graph theory [27]. Intuitively, a connected graph with ν vertices will have ( $ν - 1$ ) branches in its tree since a tree is a minimally connected graph. Consequently, by subtracting the number of branches from the total set of edges ∊, there will be ( $∊ - ν + 1$ ) chords in the cotree.

Given a graph with ν vertices and ∊ edges, the order of interconnection of a system can be summarized in a ν by ∊ incidence matrix. It is easy to construct since each edge is adjacent to exactly two vertices. The incidence matrix of $Graph (ν, ∊)$ is denoted by [κ] and is defined as follows: (1) each of the ν rows corresponds to a vertex of $Graph (ν, ∊)$ , (2) each of the ∊ columns corresponds to an edge of $Graph (ν, ∊)$ , and (3) entries $κ_{i j} = - 1$ if the ith node is the initial vertex of the ith edge, $κ_{i j} = + 1$ if the ith node is the final vertex of the ith edge, and κ_ij = 0 otherwise. It has been shown [15] that the components necessary to generate an optimum translational tree for branch coordinates should be selected in the following element order: N₁ rigid or flexible arms; N₂ position drivers (function of time); N₃ spherical, revolute, or universal joints; N₄ cylindrical or prismatic joints; N₅ rigid or deformable bodies; and N₆ force actuators.

The linear graph representation of an FMS is most easily described by means of an example. Consider the three-dimensional flexible slider-crank shown in Figure 1(a). It is obtained from a two-dimensional slider-crank (moving in the XOY plane) by rotating the joint axis at O by 45° about the global Y-axis, so that the crank OA spins to form a cone whose apex is located at O. The linear graph representation of the spatial flexible slider-crank mechanism is depicted in Figure 1(b) where the edges shown in bold comprise the spanning tree. Edge e₅₁ represents the rigid body crank and e₁₁ and e₁₂, both rigid arm elements, represent the body fixed locations of the spherical and revolute joints at O and A, respectively. Edges e₅₂ and e₁₃ model the flexible link, with flexible arm element e₁₃ locating the spherical joint e₃₁ and flexible body element e₅₂ representing the link. Edge e₅₃ represents the rigid slider body.

Figure 1:

(a) Spatial flexible slider-crank and (b) graph representation of the system.

The spherical joint at A and universal joint at B are represented by joint edges e₃₁ and e₃₂, respectively. The revolute joint O is modelled with joint edge e₃₃ and rigid arm e₁₄, which has zero length but a constant rotation matrix of the joint axis relative to the global $X Y Z$ frame. The prismatic joint at D is represented by e₄₁ for the translational motion and e₄₂ for the constant orientation of the slider relative to the ground. The input motion to the crank is represented by a motion driver e₆₁. Two edges, e₇₂ and e₇₃, corresponding to dependent virtual work elements are automatically added to the system graph to complete the virtual topological equations at joint A.

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 ${T V}$ , symbolizing quantities such as force $F$ , torque $T$ , and virtual work δW, 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 [5] that these cutset equations are obtainable by simple row operations on the incidence matrix and can be written in the following form:

[\begin{bmatrix} U_{t} & D \end{bmatrix}] {T V} = 0,

(1)

where $[D_{i j}]$ is a submatrix containing +1, −1, and 0 depending on whether element N_j is incident to N_i and submatrix U_t is a unit matrix associated with tree elements. The column matrix ${T V}$ represents the through-variables ( $F_{i}$ , $T_{i}$ , or δW_i) applied by element N_i. Considering the 3D slider-crank mechanism, the vertex equation for the node representing the end point of rigid arm e₃₃ is

{T V}_{33} + {T V}_{61} + {T V}_{51} + {T V}_{72} - {T V}_{31} = 0,

(2)

and if we use virtual work as through variables,

δ W_{33} + δ W_{61} + δ W_{51} + δ W_{72} - δ W_{31} = 0,

(3)

where $δ W_{33} = δ W_{42} = 0$ , since a frictionless joint cannot add or remove energy from the system. Furthermore, the sum of the virtual works done by the spherical joint on the crank and connecting rod must also be zero, $δ W_{72} + δ W_{73} = 0$ .

Another governing postulate is the circuit postulate which states that the algebraic sum of across-variables ${A V}$ , representing those quantities such as displacements $r$ , velocities $v$ , accelerations $\dot{v}$ , virtual angular displacement $δ θ$ , angular velocities ω, and angular accelerations $\dot{ω}$ , corresponding to all the vectors included in 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 following form:

[\begin{bmatrix} E & U_{c} \end{bmatrix}] {A V} = 0,

(4)

where submatrix $[E_{j i}]$ specifies which element N_i is included in the closed loop N_j and sub-matrix U_c is a unit matrix associated with cotree elements. The column matrix ${A V}$ represents the across-variables ( $r_{i}$ , $v_{i}$ , ${\dot{v}}_{i}$ , $δ θ_{i}$ , $ω_{i}$ , and ${\dot{ω}}_{i}$ ) for element N_i. Note that the submatrices [D] and [E] are related by the principle of orthogonality [5, 28] (not well known in dynamics):

[E] = - {[D]}^{T} or [D] = - {[E]}^{T} .

(5)

Since there is one independent circuit equation for each chord in the graph, consider the circuit equation for edge e₃₁:

\begin{matrix} {A V}_{31} - {A V}_{13} - {A V}_{32} - {A V}_{42} - {A V}_{41} \\ + {A V}_{14} + {A V}_{33} - {A V}_{12} + {A V}_{11} = 0 . \end{matrix}

(6)

Naturally, in this circuit equation, $r_{42} = r_{14} = r_{33} = r_{31} = r_{32} = 0$ . Each across-variable in (6) may be replaced with the virtual displacement $δ r$ or even the virtual rotation $δ θ$ since these vectors obey the rules for vector quantities and may be so treated. However, since the across-variables 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_{(3 ×3)} associated with each element. For example, the corresponding circuit equation for edge e₃₁ can be expressed as

[R_{31}] {[R_{13}]}^{T} {[R_{32}]}^{T} {[R_{42}]}^{T} {[R_{41}]}^{T} [R_{14}] [R_{33}] {[R_{12}]}^{T} [R_{11}] = [1],

(7)

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, $[R_{31}]$ can then be isolated as a branch transformation:

[R_{31}] = {[R_{11}]}^{T} [R_{12}] {[R_{33}]}^{T} {[R_{14}]}^{T} [R_{41}] [R_{42}] [R_{32}] [R_{13}] .

(8)

Each of the rotation matrices will be functions of the rotation variables for the corresponding element.

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

\begin{gathered} δ θ_{1} = 0, \\ δ r_{1} = δ θ_{5} \times r_{1}, \end{gathered}

(9)

where $r_{1}$ is the rigid-arm position vector which is a function of rotational body-fixed frame $θ_{5}$ since the direction of $r_{1}$ 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

δ W = δ r^{T} F .

(10)

Now, let us consider a system whose kinematics are parameterized by a branch coordinate vector $q$ where $r = r (q)$ ; then a virtual displacement $δ r$ is related to a branch coordinate variation $δ q$ by

δ r = r_{q} δ q,

(11)

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 (11) in (10):

δ W = δ q^{T} r_{q}^{T} F .

(12)

Equation (12) may be interpreted as the scalar product of a branch variation $δ q$ and its generalized force $Q$ :

δ W = δ q^{T} Q,

(13)

where the generalized force is defined as $Q \equiv r_{q}^{T} F$ . In this formulation, one may separate forces into applied forces $F^{A}$ and constraint forces $F^{C}$ such that

δ W = δ q^{T} (Q^{A} + Q^{C}) = δ W^{A} + δ W^{C} .

(14)

If branch coordinates are consistent with constraints, then $q$ is a vector with independent coordinates.

Up to now, we developed the equations of motion of a rigid component. However, rarely do we simulate the motion of a single body. This restriction has prohibited us from analyzing multibody systems which is defined as being a group of interconnected rigid and deformable components, each of which may undergo large translational and rotational motions.

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 restrictive relation, 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 multibody systems, the connectivity between different bodies is often assured by kinematic joints (such as spherical, revolute, cylindrical, and/or prismatic joints) [8]. These joints introduce geometrical constraints which impose specific relative motions between bodies. For a spatial mechanical system there is a limited number of types of functions of constraint, represented by the joints between the bodies.

Note that if constraints that act on the system are workless and if the branch coordinates $q$ must satisfy constraint equations of the form

Φ (q, t) = 0,

(15)

they are dependent. Thus branch coordinate variations must satisfy

Φ_{q} δ q = 0,

(16)

where $Φ_{q}$ is the Jacobian matrix of the constraint equation (15). If constraint forces are workless, it can be proven that

δ W^{C} = δ q^{T} Q^{C} = δ q^{T} Φ_{q}^{T} λ,

(17)

where the Lagrange multiplier $λ$ represents vector-generalized constraint forces.

Let us now consider rigid body elements N₅. The graph for this type of element consists of an edge e₅ 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 n_b rigid bodies, the kinetic energy of the ith body is defined as

T_{i} = \frac{1}{2} m_{i} v_{i}^{T} v_{i} + \frac{1}{2} ω_{i}^{' T} [I_{i}^{'}] ω_{i}^{'},

(18)

where m_i is the mass of body i, $v_{i}$ is the velocity vector of the centre of mass of the body in inertial space, $[I_{i}^{'}]$ is the inertia tensor of the body expressed in a body-fixed $^{'}$ coordinate system, and $ω_{i}^{'}$ is the angular velocity vector of the body in inertial space and expressed in the inertia $^{'}$ tensor coordinate system. Note that absolute or inertial-space velocities must be used although they may be expressed in any convenient (inertial or noninertial) coordinate system. Then, the virtual work terminal equation for the ith body becomes

δ W_{5 i} = δ q^{T} {{[\frac{d}{d t} (\frac{\partial T_{i}}{\partial \dot{q}}) - \frac{\partial T_{i}}{\partial q}]}^{T}} .

(19)

It is assumed at this point that the velocities ${v_{i}}$ and ${ω_{i}}$ 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 Lagrange's equation, the terms $(d / d t) (\partial T_{i} / \partial \dot{q})$ and $\partial T_{i} / \partial q$ will be required for each branch coordinate. By partial differentiation of (18),

\frac{\partial T_{i}}{\partial q} = (m_{i} v_{i}^{T} v_{i q} + ω_{i}^{'}^{T} [I_{i}^{'}] ω_{i q}^{'}),

(20)

and similarly,

\frac{\partial T_{i}}{\partial \dot{q}} = (m_{i} v_{i}^{T} v_{i \dot{q}} + ω_{i}^{'}^{T} [I_{i}^{'}] ω_{i \dot{q}}^{'}) .

(21)

The time derivative of (21) is

\begin{matrix} \frac{d}{d t} (\frac{\partial T_{i}}{\partial \dot{q}}) \\ = (m_{i} {\dot{v}}_{i}^{T} v_{i \dot{q}} + m_{i} v_{i}^{T} {\dot{v}}_{i \dot{q}} + {\dot{ω}}_{i}^{'}^{T} [I_{i}^{'}] ω_{i \dot{q}}^{'} + ω_{i}^{'}^{T} [I_{i}^{'}] {\dot{ω}}_{i \dot{q}}^{'}) . \end{matrix}

(22)

Substituting (20) and (22) into the variational form of Lagrange's equation (19) gives

\begin{matrix} δ W_{5 i} \\ = δ q^{T} {{[m_{i} {\dot{v}}_{i}^{T} v_{i \dot{q}} + {\dot{ω}}_{i}^{'}^{T} [I_{i}^{'}] ω_{i \dot{q}}^{'} + m_{i} v_{i}^{T} P_{i / q}^{v} + ω_{i}^{'}^{T} [I_{i}^{'}] P_{i / q}^{ω^{'}}]}^{T}}, \end{matrix}

(23)

where

\begin{gathered} P_{i / q}^{v} = {\dot{v}}_{i \dot{q}} - v_{i q}, \\ P_{i / q}^{ω^{'}} = {\dot{ω}}_{i \dot{q}}^{'} - ω_{i q}^{'} . \end{gathered}

(24)

The quantity $P_{i / q}^{v}$ can be shown to be equal to zero. Since $r_{i} = r_{i} (q, t)$ , we have $v_{i} = {\dot{r}}_{i} = r_{i q} \dot{q} + r_{i t}$ . Then, the partial derivative of $v_{i}$ with respect to $\dot{q}$ is $v_{i \dot{q}} = r_{i q}$ . Take the time derivatives of both sides of the last equation to get ${\dot{v}}_{i \dot{q}} = {\dot{r}}_{i q} = v_{i q}$ , then $P_{i / q}^{v} = 0$ . When the angular velocity vector is the exact derivative of another vector function of branch coordinates and time, the preceding argument will show that $P_{i / 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 dimensions, finite rotations cannot be represented as vectors and $P_{i / q}^{ω^{'}}$ is generally nonzero.

The virtual equation (23) is now simplified to

δ W_{5 i} = δ q^{T} {{[m_{i} {\dot{v}}_{i}^{T} v_{i \dot{q}} + {\dot{ω}}_{i}^{'}^{T} [I_{i}^{'}] ω_{i \dot{q}}^{'} + ω_{i}^{'}^{T} [I_{i}^{'}] P_{i / q}^{ω^{'}}]}^{T}} .

(25)

The right-hand side of (25) can be separated into terms containing branch accelerations $\ddot{q}$ and terms containing products of branch velocities $\dot{q} \dot{q}$ . This is necessary in order to place the coefficients of the terms $\ddot{q}$ into an augmented mass matrix ${\hat{M}}_{i}$ and the product terms into a generalized force $Q_{i}^{K}$ . Write the time-derivatives of the velocity vectors as

\begin{gathered} {\dot{v}}_{i} = v_{i \dot{q}} \ddot{q} + {\dot{v}}_{i}^{p}, \\ {\dot{ω}}_{i} = ω_{i \dot{q}}^{'} \ddot{q} + {\dot{ω}}_{i}^{'}^{p} . \end{gathered}

(26)

The terms superscripted $^{p}$ contain all products of generalized velocities. Then, exploit (26) to rewrite (25) into the final terminal equation for rigid bodies as follows:

δ W_{5 i} = δ q^{T} {{[{\hat{M}}_{i} \ddot{q} + Q_{i}^{K}]}^{T}},

(27)

where

\begin{gathered} {\hat{M}}_{i} = (m_{i} v_{i \dot{q}}^{T} v_{i \dot{q}} + ω_{i}^{'}_{\dot{q}}^{T} [I_{i}^{'}] ω_{i \dot{q}}^{'}), \end{gathered}

(28)

\begin{gathered} \begin{matrix} Q_{i}^{K} = (m_{i} {\dot{v}}_{i}^{p T} v_{i \dot{q}} + {\dot{ω}}_{i}^{'}^{p T} [I_{i}^{'}] ω_{i \dot{q}}^{'} + ω_{i}^{'}^{T} [I_{i}^{'}] P_{i / q}^{ω^{'}}) \end{matrix} . \end{gathered}

(29)

The elements of the coefficient matrix ${\hat{M}}_{i}$ are the coefficient of $\ddot{q}$ appearing in the differential equation corresponding to the branch coordinate $q$ . The forcing vector $Q_{i}^{K}$ shown in (29) 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.

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

δ q^{T} {\hat{M} \ddot{q} + Q^{K} + Q^{C} + Q^{cons} - Q^{A}} = 0 .

(30)

This variational equation of motion holds for arbitrary virtual displacement $δ q$ , so it is equivalent to the Lagrange multiplier form of constrained equations of motion.

Together, the cutset, circuit, and terminal equations form a necessary and sufficient set of motion equations for determining the time response of 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 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 [16]. Upon substitution of the branch transformation equations into these circuit equations, the constraint equations are obtained in the form of (15) which constitutes a set of c nonlinear algebraic equations. Differentiating twice (15), using the chain rule of differentiation, yields the c constraint acceleration equations:

Φ_{q} \ddot{q} = - [{(Φ_{q} \dot{q})}_{q} \dot{q} + 2 Φ_{q t} \dot{q} + Φ_{t t}] \equiv Λ .

(31)

In general, the n branch coordinates are not independent but are related by the c kinematic constraint equations (15). Thus, these constraint acceleration equations must be appended to the set of dynamic equations (30), giving ( $n + c$ ) equations to solve for branch coordinates $q$ and the c reaction loads $λ$ . Substituting (17) into (30) and combining with (31) finally yields the classical system differential-algebraic equations of motion for rigid multibody systems in matrix form,

[\begin{bmatrix} \hat{M} & Φ_{q}^{T} \\ Φ_{q} & 0 \end{bmatrix}] {\begin{Bmatrix} \ddot{q} \\ λ \end{Bmatrix}} = {\begin{Bmatrix} Q^{total} \\ Λ \end{Bmatrix}},

(32)

where $\hat{M}$ is a symmetric augmented ( $n \times n$ ) mass matrix, $Φ_{q}$ is a ( $c \times n$ ) constraint Jacobian matrix, and the total generalized force is $Q^{total} = Q^{A} - Q^{K} - Q^{cons}$ . Using (28) it is easy to form symbolically the coefficient matrix $\hat{M}$ , element by element, from partial derivatives of the velocity vectors. Nonlinearities 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 submatrices 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.

4. Flexible Multibody Systems (FMSs)

This variational graph-theoretical approach is based on the flexible models developed by Shabana [29] and Tennich [30]. In this formulation a flexible arm element, as shown in Figure 2, 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. A flexible arm element is defined as being made up of a continuum of particles that can move relative to one another. Since actual bodies are never perfectly rigid, small deformation effects are often added and can have great influence on the motion of mechanisms that are made up of multiple bodies. A flexible arm element can be located by defining the position vector $R$ of the origin of its body-fixed frame $ℜ^{'} (x^{'}, y^{'}, z^{'})$ , relative to the global reference frame $ℜ (X Y Z)$ . The global location of an arbitrary point p on a flexible body can be described as

r = R + Π s^{'} = R + Π (s_{o}^{'} + s_{f}^{'}),

(33)

where $Π$ represents an Euler transformation matrix [31] from the structure $ℜ^{'} (x^{'}, y^{'}, z^{'})$ coordinate system to the global reference frame $ℜ (X Y Z)$ . The vector $s_{o}^{'}$ is the initial undeformed position of point p with respect to the body-fixed reference frame, $s_{f}^{'}$ represents the flexible displacement of point p, and $s^{'} (= s_{o}^{'} + s_{f}^{'})$ represents a flexible arm element N₁.

Figure 2:

Flexible arm element.

Since the body-fixed frame $ℜ^{'} (x^{'}, y^{'}, z^{'})$ translates and rotates relative to the global frame, the vector $R$ and Euler transformation matrix $Π$ are functions of time. Both sides of (33) may be differentiated twice with respect to time to obtain the acceleration equation:

\ddot{r} = \ddot{R} + Π ({\ddot{s}}^{'} + \tilde{\dot{ω}} s^{'} + \tilde{ω} \tilde{ω} s^{'} + 2 \tilde{ω} {\dot{s}}^{'}),

(34)

where $\dot{ω}$ and ${\ddot{s}}^{'}$ are, respectively, the angular acceleration of the body and the second derivative of the elastic displacement of point p and $\tilde{ω}$ represents a 3×3 skew-symmetric matrix which performs a cross-product multiplication and represents the angular velocity vector of the body-fixed frame ℜ′.

To describe the deformation of a flexible body, one can discretize the structure into finite elements. If a reference frame $ℜ (e_{1}, e_{2}, e_{3})$ is attached to the jth element, the displacement of point p on a flexible body can be given by

s^{'} = N^{j} (s_{o n}^{'} + s_{f n}^{'}) = N^{j} s_{n}^{'} _{}^{j},

(35)

where $s_{o n}^{'}$ and $s_{f n}^{'}$ are the nondeformed and the nodal viscoelastic displacement of point p and $s_{n}^{'}^{j}$ is the total displacement nodal vector of the jth element. Note that $N^{j}$ is a finite element spatial interpolation matrix from the jth element to the body frame ℜ′ function of a transformation matrix from the frame $ℜ (e_{1}, e_{2}, e_{3})$ to the body reference frame $ℜ (\overset{̅}{x}, \overset{̅}{y}, \overset{̅}{z})$ .

Figures 3(a) and 3(b) represent the 3D beam and the triangular shell elements, respectively. For the 3D beam element, the two spatial interpolation function is given

N_{1}^{j e} = \frac{1}{2} (1 - ξ), N_{2}^{j e} = \frac{1}{2} (1 + ξ),

(36)

where ξ is a variable measured along each beam element, with $ξ = - 1$ at node 1 and ξ = 1 at node 2. For the triangular shell element, the three spatial interpolation functions are given by

N_{1}^{j e} = 1 - ξ - η; N_{2}^{j e} = ξ; N_{3}^{j e} = η,

(37)

where ξ and η are variables measured along the triangular element.

Figure 3:

(a) Linear finite element for a flexible beam and (b) triangular finite element for a flexible shell.

Vectors $e_{1}^{j}$ , $e_{2}^{j}$ , and $e_{3}^{j}$ , in Figure 3, form a reference surface base for the jth element which is independent of the nodal numbering. Hence, the transfer matrix between the frame $ℜ (e_{1}, e_{2}, e_{3})$ and the body reference frame $ℜ (\overset{̅}{x}, \overset{̅}{y}, \overset{̅}{z})$ can be written as Q^j = (e₁,e₂,e₃)^j. Then, the diagonal transfer matrix $T^{j}$ between the element local frame ℜQ^j = (e₁,e₂,e₃)^j and the deformable body frame $ℜ (\overset{̅}{x}, \overset{̅}{y}, \overset{̅}{z})$ can be written with six rotation tensors $Q^{j}$ (of dimension 3 × 3).

Finally, the acceleration expression of point p on the flexible body can then be rewritten from the discretized form of $s^{'}$ , from (34):

\ddot{r} = \ddot{R} - Π {\tilde{s}}^{'} Λ \ddot{θ} + Π N^{j} {\ddot{s}}_{n}^{'} ^{j} + Π Ω

(38)

with

Ω = \tilde{ω} \tilde{ω} N^{j} s_{n}^{'} ^{j} + 2 \tilde{ω} N^{j} {\dot{s}}_{n}^{'} ^{j},

(39)

where Λ represents an Euler transformation matrix for the angular velocity of the body with respect to the body reference frame ℜ′. The generalized coordinate vector $q$ for a flexible body can then be assembled from the position vector $R$ , orientation $θ$ of the origin of the body reference frame ℜ′, and the nodal coordinate vector $s_{n}^{'}$ with q = {R, θ s′_n}^T.

Consider a volume element $d v = d x d y d z$ near point p on a deformable body. The constitutive virtual work equation of this volume element dv can be separated in three different virtual work components [9]:

δ W_{5} = δ W_{5}^{inertial} - δ W_{5}^{forces} + δ W_{5}^{internal} = 0

(40)

with

\begin{gathered} δ W_{5}^{inertial} = \int_{V} δ r^{T} ρ \ddot{r} d v, \end{gathered}

(41)

\begin{matrix} δ W_{5}^{forces} = \int_{V} δ r^{T} f_{v} d v, \end{matrix}

(42)

\begin{gathered} δ W_{5}^{internal} = \int_{V} δ ϵ^{T} σ d v, \end{gathered}

(43)

where $δ ∊$ is the column matrix of varied strain components in the local frame; ( $ρ \ddot{r}$ ) is the inertial force per unit volume; $σ$ contains the corresponding internal stress components; $f_{v}$ is the body volume forces including gravity.

The virtual work of all inertial forces, for a deformable body, can be written, by using (41) for $\ddot{r}$ , under the following general form:

δ W_{5}^{inertial} = δ q^{T} M \ddot{q} - δ q^{T} Q^{q},

(44)

where $M$ and $Q^{q}$ represent the global mass matrix and the quadratic velocity vector including Coriolis term, respectively. The global mass matrix of the jth finite element of the flexible body can be assembled from the elementary mass matrices:

\begin{matrix} M^{j} & = \int_{V^{j}} ρ^{j} [\begin{bmatrix} I & - Π {\tilde{s}}^{'} Λ & Π N^{j} \\ Λ^{T} {\tilde{s}}^{'}^{T} {\tilde{s}}^{'} Λ & - Λ^{T} {\tilde{s}}^{'}^{T} N^{j} \\ Sym . & {N^{j}}^{T} N^{j} \end{bmatrix}] d v \\ = {[\begin{bmatrix} M_{R R} & M_{R θ} & M_{R f} \\ M_{θ θ} & M_{θ f} \\ Sym . & M_{f f} \end{bmatrix}]}^{j}, \end{matrix}

(45)

where $ρ^{j}$ and $V^{j}$ are, respectively, the mass density and the volume of the jth element. The elements of the mass matrix $M$ are the coefficient of $\ddot{q}$ appearing in the differential equation corresponding to the branch coordinate $q$ . Using (41), it is easy to assemble the coefficient matrix $M$ , element by element, from the transformation matrices between reference frames. The formulation of the time-invariant matrices $M_{R R}$ and $M_{f f}$ is straight-forward. The matrix $M_{R R}$ is a diagonal matrix whose elements are equal to the mass of the element. The matrix $M_{f f}$ is the conventional mass matrix that arises in any finite element analysis. Nonlinearities have been preserved throughout the formulation as each element may contain branch coordinates and quantities which vary with time. The off-diagonal mass sub-matrices represent coupling effects between translational, rotational, and flexible elements. The matrices $M_{R f}$ and $M_{θ f}$ represent the inertia coupling between gross body motion and small body deformation. These matrices, in addition to the inertia tensor $M_{θ θ}$ , are implicitly time dependent since they are functions of the body generalized coordinates.

In a similar fashion, the global quadratic velocity vectors can be assembled for the jth finite element of the flexible body:

{\begin{Bmatrix} Q_{R}^{q} \\ Q_{θ}^{q} \\ Q_{f}^{q} \end{Bmatrix}}^{j} = - \int_{V^{j}} ρ^{j} {\begin{Bmatrix} Π Ω \\ Λ^{T} {\tilde{s}}^{'}^{T} Ω \\ {N^{j}}^{T} Ω \end{Bmatrix}} d v .

(46)

The inertial forcing vector $Q^{q}$ is made up of all the terms in the equations of motion which do not contain second derivatives.

The virtual work of body forces $f_{v}$ for the jth finite element can also be written at once:

\begin{matrix} {(δ W_{6}^{forces})}^{j} & = \int_{V^{j}} δ r^{T} f_{v}^{j} d v \\ = [δ R^{T} F_{R}^{j} δ θ^{T} F_{θ}^{j} δ {s_{n}^{'} ^{j}}^{T} F_{f}^{j}] . \end{matrix}

(47)

Finally, the virtual work for internal constraints of the jth finite element was defined earlier:

{(δ W_{5}^{internal})}^{j} = \int_{V^{j}} δ {ϵ^{j}}^{T} σ^{j} d v,

(48)

where ${∊^{j}}^{T}$ and $σ^{j}$ represent, respectively, Cauchy's deformation and strain vectors. For a viscous elastic material governed by the Kelvin-Voigt model [32, 33], the behaviour law between stress and strain is established as

σ^{j} = H^{j} ϵ^{j} + G^{j} {\dot{ϵ}}^{j} = H^{j} Ξ^{j} s_{n}^{'} ^{j} + G^{j} Ξ^{j} {\dot{s}}_{n}^{'} ^{j},

(49)

where $H^{j}$ and $G^{j}$ are, respectively, the elastic and viscous tensors for this behaviour law and are function of Young's modulus and Poisson's coefficient. The matrix $Ξ^{j}$ represents spatial interpolation deformation functions. Hence, the internal virtual work for the jth finite element can be written under the following form:

{(δ W_{5}^{internal})}^{j} = δ {s_{n}^{'} ^{j}}^{T} K_{f f}^{j} s_{n}^{'} ^{j} + δ {s_{n}^{'} ^{j}}^{T} C_{f f}^{j} {\dot{s}}_{n}^{'} ^{j}

(50)

with

\begin{gathered} K_{f f}^{j} = \int_{V^{j}} {[Ξ^{T} H Ξ]}^{j} d v, \end{gathered}

(51)

\begin{gathered} C_{f f}^{j} = \int_{V^{j}} {[Ξ^{T} G Ξ]}^{j} d v, \end{gathered}

(52)

Equations (51) and (52) are, respectively, the stiffness matrix and the viscous damping matrix of the jth element. With the body kinetic and strain energy in hand, the virtual work can be exploited to generate the system equations of motion of a single flexible body. Hence, the general form of virtual work terminal equation for each body $δ W_{5 i}$ , required in the cutset equation, can be written as

δ q^{T} {M \ddot{q} + C \dot{q} + K q - F - Q} = 0 .

(53)

This variational equation holds for arbitrary virtual displacement $δ q$ , so the terms in parentheses are the well-known equations of motion in standard form. Hence, for each flexible body in the system, the translational, rotational, and viscous-elastic equations of motion can now be assembled into a partitioned matrix formulation:

\begin{matrix} [\begin{bmatrix} M_{R R} & M_{R θ} & M_{R f} \\ M_{θ θ} & M_{θ f} \\ Sym . & M_{f f} \end{bmatrix}] {\begin{Bmatrix} \ddot{R} \\ \ddot{θ} \\ {\ddot{s}}_{n}^{'} \end{Bmatrix}} + [\begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & C_{f f} \end{bmatrix}] \\ \times {\begin{Bmatrix} \dot{R} \\ \dot{θ} \\ s_{n}^{'} \end{Bmatrix}} + [\begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & K_{f f} \end{bmatrix}] {\begin{Bmatrix} R \\ θ \\ s_{n}^{'} \end{Bmatrix}} = {\begin{Bmatrix} {F Q}_{R}^{q} \\ {F Q}_{θ}^{q} \\ {F Q}_{f}^{q} \end{Bmatrix}}, \end{matrix}

(54)

where $F Q_{i}^{q} = F_{i} + Q_{i}^{q}$ for ( $i = R, θ, f$ ).

Since the total virtual work of the system of flexible bodies is the sum of the individual virtual work, the algorithm must incorporate the sum of all bodies in the formulation using a global cutset equation. Then, the variational equations of motion for flexible 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. To get the contribution of all physical components to the system, the terminal virtual work equations are written for all bodies in the tree. By summing the cutset equation, similar to (3), for all tree flexible bodies, the contribution of all physical components to the system virtual work equations is captured. The dynamic form of these expressions and fundamental properties of virtual work guarantee that a unique solution of the flexible set of equations of motion exists. For multibody systems, we must finally add, to (54), the kinematic joints which are represented by a (6×6) Boolean matrix defining the translational and rotational connexions between the bodies [30].

By exploiting GT methods and virtual work principles, this formulation has been implemented into a computer program called FLEXNET (for FLEXible NETwork). Given only a spanning tree with the terminal expressions for deformable bodies in the system, this program automatically generates the 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. Since the equations of motion for deformable bodies are function of the finite element discretization, all the constant tensors that are function of finite elements have been identified and generated by a finite element preprocessor. Hence, the preprocessor generates look-up tables that can be exploited when needed during the dynamic simulation of FMS. To enforce the constraints at the position and velocity levels, an energy algorithm proposed by Bauchau et al. [34] has been implemented.

A similar symbolic software that fully exploits this graph-theoretic approach in multidomain modelling is MapleSim, commercialized by Maplesoft. This GT formulation has been extended in MapleSim5 [35] to the analysis of mechatronic and other multidisciplinary systems such as mechanical, electrical, thermal, signal/control, and hydraulic systems which can all be naturally combined in a model diagram similar the one presented in Section 2.

5. Examples

In this section, the dynamic analyses of two different FMS are presented. Once the topology and parameters for an FMS have been defined, the translational and rotational graphs, with proper trees, can be generated automatically. Due to the systematic nature of GT methods, the previous formulation was encoded with relative ease into a general computer program. Exploiting conventional GT methods, the cutset and circuit equations are automatically generated from the given topology and for the trees selected by the user. The terminal equations developed in the preceding section are contained in a library of modelling components that can be easily updated to include new components. From the projected cutset and circuit equations, the set of motion equations governing the dynamic response of a given FMS is automatically assembled that provides insight into its structural motion.

5.1. Spatial Flexible Slider-Crank Mechanism

Let us consider the example of the spatial flexible slider-crank mechanism described earlier in Section 2. For this FMS, shown in Figure 1(a), a proper tree has been highlighted in bold in Figure 1(b). This example has been previously analyzed by Shi and McPhee with MapleSim [9] and Jonker [36]. This three-dimensional mechanism has a spherical joint at A and a universal joint at B. Crank OA and flexible link AB are initially aligned with the global X-axis. The rigid crank has a length of 0.15 m and is rotated at a constant rate ω = 150 rad/s. The coupler link is flexible with a length of 0.3 m, a circular cross-sectional area of $2.826 \times 1 0^{- 5}$ m², an area moment of inertia of $6.361 \times 1 0^{- 11}$ m⁴, a Young's modulus of $2.0 \times 1 0^{11}$ N/m², and a mass density of $7.87 \times 1 0^{3}$ kg/m³. The rigid slider has a mass of $3.34 \times 1 0^{- 2}$ kg, which is half that of the flexible link. The moment of inertia of the crank OA and the link AB are $I_{OA}^{'} = diag (0.0001,0.001,0.001)$ kg · m² and $I_{AB}^{'} = diag (0.0001,0.002,0.002)$ kg · m².

The coupler flexible link is assumed rectilinear at the initial state. The flexible link is discretized by twenty (NEL = 20) linear 3D beam elements. The deflection $μ = s / L$ at the midpoint on the link is measured perpendicular to the center of the link. The two first vibration modes are considered in this simulation. Shown in Figures 4(a) and 4(b) are the numerical results for the relative deflection μ of the mid-point on the link, projected onto the y′ and z′ axes of the local reference frame ℜ′, and plotted against the crank angle θ₃₃. The results are in good agreement with those obtained by Jonker and the software MapleSim5.

Figure 4:

Deflection along axis (a) $Y_{2}^{'}$ and (b) $Z_{2}^{'}$ .

5.2. Deployment of Two Flexible Panels

As a final example, consider the unfolding of the spatial structure, drawn in Figure 5(a), composed of two flexible panels (1) and (2) attached on a rigid base. The linear graph of this flexible structure is shown in Figure 5(b). Each panel of dimension $2 m \times 0.003 m \times 2 m$ is made of steel with Young's modulus of $2.1 \times 1 0^{11} N / m^{2}$ , Poisson's coefficient of 0.3, and a mass density of 7800 kg/m³. The edges comprising the tree are traced in bold. The deformable plates are represented by edges e₅₁ and e₅₂. The revolute joint e₃₁ connects panel (1) to the ground and revolute joint e₃₂ connects panel (2) to panel (1). Each panel is discretized by 32 triangular elements as shown on Figure 6. The flexible panel elements are represented by edges e₁₁ and e₁₂. By imposing a symmetrical meshing throughout the panels, this avoids the generation of torsion deformations outside of their respective plane. To assure good alignment of the nodes on which are fixed the revolute joints, rigid beam elements are conveniently pasted along the edges joining these nodes.

Figure 5:

(a) Deployment of the panels and (b) graph representation of the system.

Figure 6:

Mesh topology of the panels.

The first four frequency modes of vibration of each panel have been calculated as 0.60, 1.47, 3.71, and 5.53 Hz. The complete deployment of the flexible panels has been achieved in T = 14 seconds where angles ψ₁ and ψ₂ go from 0° to 90° and from 180° to 0°, respectively. The imposed drivers at the articulations of the structure have zero velocity and zero acceleration at the beginning and end of the simulation. The angular drivers e₆₁ and e₆₂ are given by

\begin{gathered} ψ_{1} (t) = \frac{π}{2 T} [t - \frac{T}{2 π} \sin (\frac{2 π t}{T})], \\ ψ_{2} (t) = π - \frac{π}{T} [t - \frac{T}{2 π} \sin (\frac{2 π t}{T})], \end{gathered}

(55)

where $ψ_{1} (t) = π / 2$ and $ψ_{2} (t) = 0$ when $t \geq T$ .

Simulations for the flexible panels were performed during t = 20 seconds. Results are plotted from the time of release of the panels (at t = 0 s) through complete deployment (at t = 14 s) and rebound effects of the panels (until t = 20 s). Since MapleSim5 is restricted to the simulation of flexible links only, Figure 7 provides a comparison with the software ADAMS/Flex [37] for the deflection of the centers of the two panels considering the first four vibration modes.

Figure 7:

Transversal deflection of panels 1 and 2 centers.

6. Conclusion

By combining the mechanical system topology with the variational virtual work constitutive equations, a new systematic graph-theoretic formulation has been introduced and used to describe the time-varying configuration of spatial FMS. 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 GT self-formulating algorithm presented in this work illustrates a general applicable method for deriving the differential equations of motion for FMS. In the slider-crank mechanism, a constraint equation function of branch coordinates was retained and incorporated in the differential equations of motion. In the two-panel example, it was possible to define independent branch coordinates that lead to nonlinear differential equations of motions. Even in these relatively elementary examples, extremely nonlinear differential equations are obtained numerically, precluding practical generation of analytical solutions. Numerical integration techniques are therefore, generally, needed to solve these equations of motion.

For open-loop systems, like the deployment of panels, 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, like in the case of the slider-crank mechanism, a tree structure is formed by mathematically cutting the constraining elements or joints yielding the constraint circuit equations. A kinematic formulation may then be developed for the resulting spanning tree, so it 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 $δ q$ that are consistent with the constraint. It is, therefore, necessary to introduce the equation associated with the physical constraint.

In order to derive the equations of multibody systems with structural flexibility, the deformable bodies were discretized with two simple finite elements (3D beam and triangular shell) with $C^{°}$ continuity. It has been shown, in this work, that accurate results can be obtained by using very simple finite elements of continuity $C^{°}$ . However, the use of finite elements naturally generates a large set of elasticity equations. To limit the number of coordinates, the equations of elastic deformations were projected onto the modal base of the deformable bodies to reduce the dimension of the differential algebraic equations.

Footnotes

Acknowledgment

Financial support of this paper by the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged.

References

Tutte

W. T.

, Graph Theory, Cambridge University Press, New York, NY, USA, 2001.

Gross

J. L.

and Yellen

, Graph Theory and Its Applications, Chapman & Hall/CRC, Boca Raton, Fla, USA, 2nd edition, 2005.

Even

, Graph Algorithms, Computer Science Press, 1979.

Koenig

H. E.

and Blackwell

W. S.

, “Linear graph theory—a fundamental engineering discipline,” IRE Transaction on Education, vol. 3, pp. 42–62, 1960.

Andrews

G. C.

, The vector-network model: a topological approach to mechanics, Ph.D. thesis, University of Waterloo, Waterloo, Ontario, Canada, 1971.

Koenig

H. E.

Tokad

, and Kesavan

H. K.

, Analysis of Discrete Physical Systems, McGraw-Hill, 1967.

Andrews

G. C.

and Kesavan

H. K.

, “The vector-network model: a new approach to vector dynamics,” Mechanism and Machine Theory, vol. 10, no. 1, pp. 57–75, 1975.

Richard

M. J.

, Dynamic simulation of constrained three dimensional multibody systems using vector network techniques, Ph.D. thesis, Queen's University, Kingston, Ontario, Canada, 1985.

Shi

and McPhee

, “Dynamics of flexible multibody systems using virtual work and linear graph theory,” Multibody System Dynamics, vol. 4, no. 4, pp. 355–381, 2000.

10.

Shi

McPhee

, and Heppler

G. R.

, “A deformation field for Euler-Bernoulli beams with applications to flexible multibody dynamics,” Multibody System Dynamics, vol. 5, no. 1, pp. 79–104, 2001.

11.

Richard

M. J.

Bouazara

, and Therien

J. N.

, “Analysis of multibody systems with flexible plates using variational graph-theoretic methods,” Multibody System Dynamics, vol. 25, no. 1, pp. 43–63, 2011.

12.

Wasfy

T. M.

and Noor

A. K.

, “Computational strategies for flexible multibody systems,” Applied Mechanics Reviews, vol. 56, no. 6, pp. 553–613, 2003.

13.

McPhee

J. J.

and Redmond

S. M.

, “Modelling multibody systems with indirect coordinates,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 50–51, pp. 6942–6957, 2006.

14.

Richard

M. J.

Huang

M. Z.

, and Bouazara

, “Computer aided analysis and optimal design of mechanical systems using vector-network techniques,” Applied Mathematics and Computation, vol. 157, no. 1, pp. 175–200, 2004.

15.

Shi

and McPhee

, “On the use of virtual work in a graph-theoretic formulation for multibody dynamics,” in Proceedings of the ASME Design Engineering Technical Conference, Sacramento, Caif, USA, 1997, DETC97/VIB-4199.

16.

McPhee

J. J.

, “Automatic generation of motion equations for planar mechanical systems using the new set of branch coordinates,” Mechanism and Machine Theory, vol. 33, no. 6, pp. 805–823, 1998.

17.

Behzad

and Chartrand

, Introduction to the Theory of Graphs, Allyn and Bacon, 1971.

18.

Christofides

, Graph Theory, An Algorithmic Approach, Academic Press, New York, NY, USA, 1975.

19.

Roberson

R. E.

, “The path matrix of a graph, its construction and its use in evaluating certain products,” Computer Methods in Applied Mechanics and Engineering, vol. 42, no. 1, pp. 47–56, 1984.

20.

Arczewski

, “Application of graph theory to the mathematical modelling of a class of rigid body systems,” Journal of the Franklin Institute, vol. 327, no. 2, pp. 209–223, 1990.

21.

Baciu

Chou

J. C. K.

, and Kesavan

H. K.

, “Constrained multibody systems: graph-theoretic Newton-Euler formulation,” IEEE Transactions on Systems, Man and Cybernetics, vol. 20, no. 5, pp. 1025–1048, 1990.

22.

Arczewski

, “Application of graph theory to the determination of kinetic energy of rigid body systems,” Journal of the Franklin Institute, vol. 324, no. 3, pp. 351–367, 1987.

23.

Chou

J. C. K.

Kesavan

H. K.

, and Singhal

, “Dynamics of 3-D isolated rigid-body systems: graph-theoretic models,” Mechanism and Machine Theory, vol. 21, no. 3, pp. 261–272, 1986.

24.

Bos

A. M.

, Modelling Multibody Systems in Terms of Multibond Graphs with Application to a Motorcycle, Dissertation Twente University, Enschede, The Netherlands, 1986.

25.

, “Applications of bond graphs and vector bond graphs to rigid body dynamics,” Journal of China Textile University, vol. 5, no. 4, p. 67, 1988.

26.

Margolis

D. L.

, “Bond graphs for automated simulation and control of nonlinear vehicle systems,” in Proceedings of the Future Transportation Technology Conference and Exposition, p. 8, Seattle, Wash, USA, 1987, SAE paper no. 871558.

27.

Andrews

G. C.

, A General Re-statement of the Laws of Dynamics Based on Graph Theory, Problem Analysis in Science and Engineering, Academic Press, 1977.

28.

Schmitke

and McPhee

, “Using linear graph theory and the principle of orthogonality to model multibody, multi-domain systems,” Advanced Engineering Informatics, vol. 22, no. 2, pp. 147–160, 2008.

29.

Shabana

A. A.

, “Transient analysis of flexible multi-body systems. Part I: dynamics of flexible bodies,” Computer Methods in Applied Mechanics and Engineering, vol. 54, no. 1, pp. 75–91, 1986.

30.

Tennich

, Dynamique de systèmes multi-corps exibles, une approche générale, Ph.D. dissertation, Laval University, Québec, Canda, 1994.

31.

Wittenburg

, Dynamics of Systems of Rigid Bodies, Teubner, Stuttgart, Germany, 1977.

32.

Flugge

, Viscoelasticity, Blaisdell, New York, NY, USA, 1967.

33.

Christensen

R. M.

, Theory of Viscoelasticity: An Introduction, Academic Press, New York, NY, USA, 1975.

34.

Bauchau

O. A.

Damilano

, and Theron

N. J.

, “Numerical integration of non-linear elastic multi-body systems,” International Journal for Numerical Methods in Engineering, vol. 38, no. 16, pp. 2727–2751, 1995.

35.

MapleSim5, High-Performance Multi-Domain Modeling and Simulation, 2011, http://www.maplesoft.com/products/maplesim/index.aspx.

36.

Jonker

, “A finite element dynamic analysis of spatial mechanisms with flexible links,” Computer Methods in Applied Mechanics and Engineering, vol. 76, no. 1, pp. 17–40, 1989.

37.

ADAMS, MSC/Software Simulating Reality, Delivering Certainty, 2011, http://www.mscsoftware.com/Products/Modeling-Solutions/Default.aspx.