Simulation of Flexible Multibody Systems Using Linear Graph Theory

This chapter provides a general description of a variational graph-theoretic formulation for simulation of flexible multibody systems (FMS) 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. 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 3-D 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, deployable space structures and micro-nano electro-mechanical systems (MEMS).


Introduction
This chapter provides a general description of a variational graph-theoretic formulation for simulation of flexible multibody systems (FMS) 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.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 3-D 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, deployable space structures and micro-nano electro-mechanical systems (MEMS).
Realistic dynamic simulation of industrial mechanisms that requires tracking accuracy at high operational speed is becoming increasingly important for engineers.Hence, to accurately describe such motions, the effects of flexibility and damping must be included in the dynamic model.Since the equations governing the motion of FMS are highly non-linear and dynamically coupled, one must exploit some kind of linear graph principles (Andrews, 1971;Behzad & Chartrand, 1971;Christofides, 1975;Even, 1979;Koenig & Blackwell, 1960) to properly define the interconnection between the bodies.By combining linear graph theory (Andrews & Kesavan, 1975;Koenig et al., 1967;Richard, 1985) with the principle of virtual work (Richard et al., 2011;Shi & McPhee, 2000;Shi et al., 2001) and finite elements, a dynamic formulation is obtained that extends graph-theoretic (GT) modelling methods to the analysis of 3-D beams and shell surfaces of FMS.The widespread interest in flexible multibody systems (FMS) is evidenced by the existence of a large number of algorithms (Wasfy & Noor, 2003).It has been shown (McPhee & Redmond, 2006;Richard et al., 2004;Shi & McPhee, 1997) 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 non-linear constitutive equations for individual components, resulting in very modular and efficient algorithms.
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 (McPhee, 1998) 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.Hence, by combining linear graph theory and 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 FMS.

System representation by linear graph
Many researchers have studied the theory of graphs (Arczewski, 1990;Baciu et al., 1990;Chou et al., 1986;Roberson, 1984) and bond-graphs (Bos, 1986;Hu, 1988) 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 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 chord.This connotes that the subgraph is connected and contains all the nodes of a given system graph, but has no closed loops.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 has been shown (McPhee, 1998) that the components necessary to generate an optimum tree for branch coordinates should be selected in the following order: N 1 rigid or flexible arms; N 2 position drivers (function of time); N 3 spherical or revolute joints; N 4 cylindrical or prismatic joints; N 5 rigid or deformable bodies, N 6 force actuators and N 7 virtual work elements.
The linear graph representation of a mechanical system is most easily described by means of an example.Consider the simple planar dynamic system, shown in figure 1(a), which consists of a body with center of gravity located at C.G. acted upon by two springs and a dashpot.Assuming negligible weight, the vector-network diagram, including the body traced in dashed lines to make the network easier to identify, is depicted in figure 1(b) and portrays a linear graph with six vertices and eight edges.
The edge e 1 is the inertial displacement vector representing the center of gravity, e 2 and e 3 define rigid or flexible arm elements, e 4 and e 5 specify displacement drivers while e 6 and e 7 Fig. 1. a) Simple planar mechanical system and b) graph representation of the system represent the springs and e 8 the dashpot.All vectors such as e 1 , e 4 and e 5 , the properties of which are established from an inertial frame, must emanate from the single ground node (datum) since that node is the only absolute fixed reference in the diagram, as compared to e 2 and e 3 which are defined relative to the body.Visibly, edges e 6 , e 7 and e 8 span their respective two points of interconnection in the system.
It is possible to determine the number of branches and chords in a given graph Graph(ν, ǫ) by applying basic theorems of graph theory (Andrews, 1977).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 κ ij = −1 if the i th node is the initial vertex of the j th edge, κ ij =+ 1i ft h ei th node is the final vertex of the j th edge, and κ ij = 0 otherwise.All columns contain exactly two non-zero entries and (ν − 2) zero entries.The incidence matrix can always be compiled from inspection of the graph.As an illustration, the graph of figure 1(b) has the following incidence matrix: From the incidence matrix, one can generate a cutset matrix.A cutset is a disconnecting set of edges such that after removal of that set of edges, the graph is divided into two or more components.A fundamental cutset (f-cutset) is considered a minimal cutset.A f-cutset reduces a connected graph into exactly two components.The selection of a tree branch will always identify one f-cutset associated to the branch and if a single cotree chord of the f-cutset is re-introduced into the graph, it unites the two residual subgraphs into a connected graph.As opposed to trees which represent a minimal set of edges which connect all the vertices of Graph(ν, ǫ), a f-cutset is a minimal set of edges which disconnect some vertices from others.
Since a f-cutset exists for each tree branch, there will be (ν − 1) f-cutsets in a graph with ν vertices.Therefore, all f-cutsets can be assembled in a (ν − 1) by ǫ cutset matrix [D] identifying all cotree terminal graphs acting through each vertex.For example, the cutset matrix of the diagram sketched in figure 1(b) can be obtained by simple row operations (Andrews, 1971) performed on the (ν − 1) rows of the incidence matrix.In this case, rows (2) and (3) are added to row (1) leading to the five branches cutset matrix, tree cotree where [U t ]i sa( ν − 1) by (ν − 1) unit sub-matrix associated with tree branches and [D]i s a( ν − 1) by (ǫ − ν + 1) sub-matrix associated with cotree chords.Since a f-cutset isolates a part of the system, its application to dynamic systems will become obvious when solving the relationship among forces which requires the construction of a "free-body diagram".
From the cutset matrix, one can generate a circuit matrix.A circuit is a connected subgraph in which exactly two edges are incident with each vertex.This concept will be molded into graph-networks through the use of a fundamental-circuit (f-circuit) which is defined as a subgraph which contains a single cotree chord, forming a closed chain, with tree branches.Following this definition, each chord of the cotree will form a f-circuit, thereby producing an independent set of closed chains since at least one edge will not be found in any other circuit.Earlier, the exact number of cotree chords in a Graph(ν, ǫ) was found to be (ǫ − ν + 1).

352
New Frontiers in Graph Theory www.intechopen.com Combining this result with the fact that for each chord there is a corresponding f-circuit, the total number of independent circuits in a Graph(ν, ǫ) is (ǫ − ν + 1).The total set of f-circuits can be generated in a (ǫ − ν + 1) by ǫ circuit matrix [E] specifying each inertial terminal graph compatible with each circuit.The circuit matrix can be obtained by applying another basic theorem (Andrews, 1971) of graph theory which proves that the cutset and circuit matrices are orthogonal.This orthogonal relationship states that the scalar product of the cutset matrix and circuit matrix must vanish.Hence, in matrix notation, the submatrices [D] and [E] are related by the negative transpose principle of orthogonality (Andrews & Kesavan, 1975;Richard, 1985) (not well-known in dynamics): This liaison regulates the entire structure of this GT formulation.This principle can be demonstrated by using the sample graph of figure 1(b).Relation (3) can now be exploited to automatically transform the cutset matrix into the three chords circuit matrix, tree cotree where [U c ]i sa( ǫ − ν + 1) by (ǫ − ν + 1) unit sub-matrix associated with cotree chords and [E]isa(ǫ − ν + 1) by (ν − 1) sub-matrix associated with tree branches.In mechanical systems, the limitations to the freedom of movement of a body are specified by certain compatibility criteria which can be extricated from the circuit matrix.This matrix will prove to be a necessary mathematical tool in the resolution of dynamic systems since it provides some internal information relevant to the geometry of contact between terminal graphs.
At this point, one must introduce a physical meaning to these matrices.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 {TV}, 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.It has been shown earlier that these cutset equations are obtainable by simple row operations on the incidence matrix, and can be written in the form: where [D ij ] is a sub-matrix containing +1, -1 and 0 depending whether element N j is incident to N i and sub-matrix U t is a unit matrix associated tree elements.The column matrix {TV} represents the through-variables (F i , T i or δW i ) applied by element N i .For the example depicted in figure 1(a), by combining equations ( 2) and ( 5), the vertex equation for node 1 representing the center of mass of the body is In this case, conventional through variables, such as force vectors, are introduced in equation ( 6) (where F 1 = −m r represents the inertial force).
The other governing postulate is the circuit postulate which states that the algebraic sum of across-variables {AV }, representing those quantities such as displacements r, velocities v, accelerations v, virtual angular displacement δθ, angular velocities ω and angular accelerations ω, 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 form: where sub-matrix [E ji ] 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 {AV } represents the across-variables (r i , v i , vi , δθ i , ω i and ωi ) for element N i .Since there is one independent circuit equation for each chord in the graph, the circuit equation for edge e 8 is where r represents the translational displacement of the corresponding element.
Consider the planar four-bar mechanism shown in figure 2(a).The system consists of three moving bodies (plus one fixed link) and four revolute joints.The link OA that is connected to the power source is called the input link (or crank).A driving torque causes the crank to rotate with angular velocity ω.The output link connects the moving pivot B to the ground pivot C. The coupler link connects the two moving pivots, A and B. The linear graph representation of the planar flexible four-bar linkage is depicted in figure 2(b) where the edges shown in bold comprise the spanning tree.
In this GT model, nodes (or vertices) are used to represent reference frames in the system while edges (lines) represent physical elements that connects these frames.Edges e 51 , e 52 and e 53 represent the bodies in the system and e 11 , e 12 , both rigid arm elements, representing the body-fixed locations of the revolute joints at O and A, respectively.Edges e 13 and e 14 model the flexible links.The four revolute joints are modelled with joint edges e 3i (i = 1, ..., 4).Edge e 21 represents a position driver and e 61 is the driving torque.Finally, seven edges e 7i (i = 1, ..., 7) corresponding to dependent virtual work elements are automatically added to the system graph to complete the virtual topological equations.
Since we can also use virtual work δW as a through variable, the cutset equation for the tree joint element e 34 , for example, is where each term corresponds directly to a physical component.work elements should automatically be included in the diagram model.Note that δW 33 =0, since a frictionless joint cannot add or remove energy from a system.Furthermore, the sum of the virtual works done by joint B must also be zero, hence, δW 74 + δW 75 = 0.This fact can be exploited in this GT formulation to eliminate all joint reactions from the constitutive equations.Later, it will be shown that joint reactions can be retained in the GT algorithm by means of Lagrange multipliers.
The second set of topological equations, the circuit equations, can be generated automatically from the cutset equations.As an example, the circuit equation for cotree edge e 33 is r 33 − r 14 − r 34 − r 21 + r 31 − r 11 + r 12 + r 32 − r 13 = 0 (10) Clearly, this linear equation represents the vector loop closure condition for the circuit containing e 33 where r 31 = r 32 = r 33 = r 34 = 0. Thus, the circuit equation ( 7) represent a set of linearly independent equations, one for each chord in the cotree.

Virtual work terminal equations for rigid bodies
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 where r 1 is the rigid-arm position vector which is 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 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 (14 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 equation ( 14) in equation ( 13), Equation ( 13) may be interpreted as the scalar product of a branch variation δq and its generalized force where the generalized force is defined as Q ≡ r T q F. In this formulation, one may separate forces into applied forces F A and constraint forces F C such that, 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 they are dependent.Thus branch coordinate variations must satisfy where Φ q is the Jacobian matrix of the constraint equation ( 18).If constraint forces are workless, it can be proven that where the Lagrange multiplier λ 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 will be treated in the next section) 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 equation ( 17).For elements containing physical constraints the terminal equation ( 20) 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 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 equation ( 18), 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.Figure 3 presents the four most common types of ideal joints which can be found in multibody systems.
A spherical joint shown in figure 3(a) is defined by the condition that the center of the ball coincides with the center of the socket.The terminal equations are, then, 3 scalar constraint equations that restricts the relative positions between the bodies, [ e x e y e z ] T δr 3 =[000 ] T where e x , e y , e z are unit vectors and δr 3 represents the relative virtual displacement between the bodies.Essentially, the three constraining equations ( 22) impose that there is no relative displacements at the joint in the x − y − z directions.
A revolute joint, shown in figure 3(b), 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 At this point, one can represent these kinematic joints in a (6 × 6) Boolean matrix, where E k t and E k r are two diagonal (3 × 3) tensors which define the translational and rotational connexions between bodies.
Let us now consider rigid body elements N 5 .The graph for this type of element consists of an edge e 5 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 i th body is defined as 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 non-inertial) coordinate system.Then, the virtual work terminal equation for the i th body becomes 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 dt ∂T i ∂ q and ∂T i ∂q will be required for each branch coordinates.By partial differentiation of equation ( 27), and similarly, The time derivative of equation ( 30) is d dt Substituting equations ( 29) and (31) into the variational form of Lagrange's equation ( 28) gives where The quantity P v i/q can be shown to be equal to zero.Since r i = r i (q, t), we have From equation ( 35), form the partial derivative of v i with respect to q v i q = r iq (36) Take the time derivatives of both sides of equation ( 36) to get Substitute equation (37) into equation ( 33) to show that P v i/q = 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 non-zero.
The virtual equation ( 32) is now simplified to The right hand side of equation ( 38) can be separated into terms containing branch accelerations q and terms containing products of branch velocities qq.This is necessary in order to place the coefficients of the terms q into an augmented mass matrix Mi and the product terms into a generalized force Q K i .Write the time-derivatives of the velocity vectors as vi = v i q q + vp i (39 The terms superscripted p contain all products of generalized velocities.Then, exploit equations ( 39) and ( 40) to rewrite equation ( 38) into the final terminal equation for rigid bodies as follows where and The elements of the coefficient matrix Mi are the coefficient of q appearing in the differential equation corresponding to the branch coordinate q.The forcing vector Q K i shown in equation ( 43) 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.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 This variational equation of motion holds for arbitrary virtual displacement δq,s oi ti s 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 function of branch coordinates.
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 (McPhee, 1998).Upon substitution of the branch transformation equations these circuit equations, the constraint equations are obtained in the form of equation ( 18) which constitutes a set of c nonlinear algebraic equations.Differentiating twice equation ( 18), using the chain rule of differentiation, yields the c constraint acceleration equations Φ q q = − (Φ q q) q q + 2Φ qt q + Φ tt ≡ Λ (45) In general, the n branch coordinates are not independent, but are related by the c kinematic constraint equations ( 18).Thus, these constraint acceleration equations must be appended to the set of dynamic equations (44), giving (n + c) equations to solve for branch coordinates q and the c reaction loads λ.Substituting equation ( 18) into equation ( 44) and combining with equation ( 45) finally yields the classical system differential-algebraic equations of motion for rigid multibody systems in matrix form, where M is a symmetric augmented (n × n) mass matrix, Φ q isa(c × n) constraint Jacobian matrix and the total generalized force 42) it is easy to form symbolically the coefficient matrix 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 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.

Flexible multibody systems (FMS)
This variational graph-theoretical approach is based on the flexible models developped by (Shabana, 1986) and (Tennich, 1994).In this formulation a flexible arm element, as shown in figure 4, 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 The global location of an arbitrary point P on a flexible body can be described as where Π represents an Euler parameter transformation matrix (Wittenburg, 1977) from the structure ℜ ′ (X ′ , Y ′ , Z ′ ) coordinate system to the global reference frame ℜ(XYZ).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 eq.( 47) may be differentiated twice with respect to time to obtain the acceleration equation, where ω and s′ are, respectively, angular acceleration of the body and the second derivative of the elastic displacement of point p and ω 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 j th element, the displacement of point p on a flexible body can be given by with where s ′ on and s ′ fn are the non-deformed and the nodal visco-elastic displacement of point p and s ′j n is the total displacement nodal vector of the j th element.Note that N j is a finite element spatial interpolation function from the j th element to the body frame ′ where Q j represents a transformation matrix from the frame ℜ(e 1 , e 2 , e 3 ) to the body reference frame ℜ( x, ȳ, z) and T j is a transfert matrix assembled with the Q j matrices.
Figures 5(a) and 5(b) represent the 3-D beam and the triangular shell elements, respectively.For the 3-D beam element, the two spatial interpolation function is given where ξ is an adimensional variable mesured 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 where ξ and η are adimensional variables mesured along the triangular element.5, form a reference surface base for the j th element which is independent of the nodal numbering.Hence, the transfert matrix between the frame ℜ(e 1 , e 2 , e 3 ) to the body reference frame ℜ( x, ȳ, z) can be written as Q j =(e 1 , e 2 , e 3 ) j .Then, the 363 Simulation of Flexible Multibody Systems Using Linear Graph Theory www.intechopen.comdiagonal transfert matrix T j between the element local frame ℜ(e 1 , e 2 , e 3 ) j and the deformable body frame ℜ( x, ȳ, 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 eq.( 48), where ω = Λ θ with Λ representing 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 dv = dx dy dz 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 (Shi et al., 2001;Tennich, 1994), where δǫ is the column matrix of varied strain components in the local frame; (ρ r) is the inertial force per unit volume; σ contains the corresponding internal stress components and f v is the body volume forces including gravity.
For a deformable body, the virtual work of all inertial forces can be written from equation ( 54) under the following general form, 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 j th finite element of the flexible body can be assembled from the elementary mass matrices, where ρ j and V j are, respectively, the mass density and the volume of the j th element.The elements of the mass matrix M are the coefficient of q appearing in the differential equation corresponding to the branch coordinate q.Using equation ( 53), 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 RR and M ff is straight-forward.The matrix M RR is a diagonal matrix whose elements are equal to the mass of the element.The matrix M ff is the conventional mass matrix that arises in any finite element analysis.Non-linearities 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 Rf 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 j th finite element of the flexible body, 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 j th finite element can also be written at once, Finally, the virtual work for internal constraints of the j th finite element was defined earlier, where ǫ j T and σ j represent, respectively, Cauchy's deformation and strain vectors.For a viscous elastic material governed by the Kelvin-Voigt model (Christensen, 1975;Flugge, 1967), the behaviour law between stress and strain is established as, 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 j th finite element can be written under the following form, Equations ( 64) and ( 65) are, respectively, the stiffness matrix and the viscous damping matrix of the j th 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 5i , required in the cutset equation, can be written 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, where 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 eq.( 5), for all tree flexible bodies, the contribution of all physical components to the system virtual work equations are 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.
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 is greater than the degrees of freedom, then 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, similar to eq.( 7), onto the joint reaction space.Then, if kinematic joint k establishes a translational connexion between points a and b located on bodies i and j, the following vectorial constraint expression must be respected, where r i and r j are the position vectors of the fixation points of the kinematic joint k on bodies i and j while r ′k k represents the length of the joint.In a similar fashion, the rotational connexion 366 New Frontiers in Graph Theory www.intechopen.comat kinematic joint k between points a and b can be written, where ω i , ω j and ω k k are angular velocities of bodies i and j and kinematic joint k, expressed in the joint reference frame.By introducing the tensor E k , obtained in equation ( 26), to eliminate all superflous constraints and substituting the relations for the derivatives of the transformation matrices, the translational and rotational geometrical constraint equations ( 68) and ( 69) can be derived twice with respect to time to obtain, with Ω = 2Π ω ṡ′ n + Π ω ωs ′ .Together, the topological and terminal equations form a necessary and sufficient set of motion equations for FMS.The variational equations of motion for FMS 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 a complete set of equations of motion for the FMS.In general, the branch coordinates are not independent, but are related by the kinematic constraint equations.Thus, these constraint acceleration equations (70) must be appended to the set of dynamic equations ( 67), giving all the equations to solve for branch coordinates and finally yields the classical system differential-algebraic equations of motion for FMS.

Examples
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 a good 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., 1995) has also been implemented.
A similar symbolic software that fully exploits this graph-theoretic approach in multi-domain modelling is MapleSim, commercialized by Maplesoft.This GT formulation has been extended in the second version of the software (MapleSimII, 2011) 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. However, this software is limited to the simulation of rigid bodies and flexible beams only.

Flexible four-bar mechanism
Let us consider the example of the planar flexible four-bar mechanism described earlier in section 2. For this FMS, shown in figure 2(a), a proper tree has been highlighted in bold in figure 2(b).This example has been previously analyzed by (Khulief, 1992) and was also analyzed with the software MapleSimII.The rigid crank OA has a length of 0.3 m and is driven at a constant angular speed of ω = 210 rad/s.The flexible coupler link AB and the flexible follower BC are discretized by eight linear 3-D beam elements with bending moments of inertia equal to 3.112 × 10 −8 m 4 and have cross-sectional areas of 1.767 × 10 −4 m 2 .Flexible links AB and BC are made of steel with modulus of elasticity of 2.0 × 10 11 N/m 2 with mass densities of 7800 kg/m 3 .The rigid crank is assumed horizontal at the initial state.The at the mid-point of the links are mesured perpendicular to the initial links.The first two vibration modes are considered in this simulation.Shown in the figures 6(a) and 6(b) are the numerical results for the relative deflection of the mid-point of the links plotted against the crank angle θ 31 .
Once the topology and parameters for the FMS has been defined, the translational and rotational graphs 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.The terminal equations developped in the preceeding 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.The results are in good agreement with those obtained by (Khulief, 1992) and the software (MapleSimII, 2011).

Deployment of two flexible panels
As a final example, consider the unfolding of the spatial structure, drawn in figure 7(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 7(b).Each panel of dimension 2m × 0.003m × 2m is made of steel with a Young's modulus of 2.1 × 10 11 N/m 2 , Poisson's coefficient of 0.3 and a mass density of 7800 kg/m 3 .edges comprising the tree are traced in bold.The deformable plates are represented by e 51 and e 52 .The revolute joint e 31 connects panel (1) to the ground and revolute joint e 32 connects panel (2) to panel (1).Each panel is discretized by 32 triangular elements.The flexible panel elements are represented by edges e 11 and e 12 .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.
The complete deployment of the flexible panels has been achieved in T = 14 seconds where angles ψ 1 and ψ 2 goes from 0 o to 90 o and from 180 o to 0 o , respectively.The imposed drivers at the articulations of the structure have zero velocity and zero acceleration at the beginning where ψ 1 (t)=π/2 and ψ 2 (t)=0 when t ≥ 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 = 0s) through complete deployment (at t = 14s) and rebound effects of the panels (until t = 20s).Since MapleSimII is restricted to the simulation of flexible links only, figure 8 provides a comparison with the software (Adams/flex, 2011) for the deflection of the centers of the two panels considering the first four vibration modes.

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.
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 four-bar 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 that are consistent with the constraint.It is, therefore, necessary to introduce the equation associated with the physical constraint.
Through graph-theoretic methods, the state-of-the-art of general multibody programs has advanced to the point where the flexibility of bodies combined to multi-domain systems can be simulated.Perhaps the most important requirement of a GT general purpose multibody computer program is the quality of the interfaces for the input and output of data.Hence, the self-formulating aspect of all programs exploiting GT methods will become very important for a productive utilization.Other features of GT methods worth mentioning for the future are the portability of the GT algorithms which should be able to adapt easily to all computer environments and the compatibility between different databases of multidisciplinary programs.

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

Fig. 2
Fig. 2. a) Planar four-bar mechanism and b) graph representation of the system

Fig. 3 .
Fig. 3. Ideal kinematic joints a) spherical, b) revolute, c) cylindrical and d) prismatic five terminal constrained equations are [ e x e y e z ] T δr 3 =[000 ] T ; [ e x e y ] T δθ 3 =[00 ] T (23) where δθ 3 represents the relative virtual angular rotation between the bodies.A cylindrical joint between a pair of bodies is shown in figure 3(c).It permits relative translation and relative rotation between bodies about a common axis.The four terminal constrained equations are [ e x e y ] T δr 4 =[00 ] T ; [ e x e y ] T δθ 4 =[00 ] T (24) A prismatic joint, shown in figure 3(d), 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 [ e x e y ] T δr 4 =[00 ] T ; [ e x e y e z ] T δθ 4 =[000 ] T (25)

Fig. 4 .
Fig. 4. Flexible arm element 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 ) a flexible arm element N 1 .

Fig
Fig. 5. a) Flexible beam linear finite element b) Flexible shell triangular finite element Vectors e j 1 , e j 2 and e j 3 , in figure5, form a reference surface base for the j th element which is independent of the nodal numbering.Hence, the transfert matrix between the frame ℜ(e 1 , e 2 , e 3 ) to the body reference frame ℜ( x, ȳ, z) can be written as Q j =(e 1 , e 2 , e 3 ) j .Then, the

Fig. 6 .
Fig. 6.Deflection of center of a) link AB and b) link BC and end of the simulation.The angular drivers e 61 and e 62 are given by,ψ 1 (t)= π 2T t − T 2π sin( 2πt T )(71)

Fig. 7 Fig. 8 .
Fig. 7. a) Deployment of the panels and b) graph representation of the system