On the Stiffness Analysis and Elastodynamics of Parallel Kinematic Machines

robotics important of engineering of as & mathematics and mechanism The in been steadily increasing during the decades. This concern has directly the development of the novel theoretical research and This new book provides about fundamental topics of serial and parallel manipulators such as kinematics & dynamics modeling, optimization, control algorithms and design strategies. I would like to thank all have the book chapters with their valuable novel ideas and current developments.

for large rotations and small deformations, while the ANCF is preferred for large rotations and large deformations. Notwithstanding, both the two formulations usually require great computational efforts to be extended to complex robots or to optimization techniques involving the whole robot's workspace. In this chapter we propose a formulation to study the elastostatics and elastodynamics of PKMs. The method is linear and tries to combine some feature existing in the literature to build a solid framework, the outlines of which are described in Section 2. We start from a MSA approach based on the minimization of the strain energy in which all joints are introduced by means of constraints between nodal displacements. In this way we avoid Lagrange multipliers and the introduction of joints becomes straightforward. Unlike VJM methods based on lumped stiffness, we use 3D Euler beams to simulate links and to distribute stiffness to the flexible structure of a robot, as described in Section 3. The same set of nodal displacement arrays is used to obtain the generalized stiffness and inertia matrix, the latter being lumped or distributed, as discussed in Sections 4 and 5. Besides, the flexibility of the proposed method allows us to provide useful extension to compliant PKMs with joint flexures. The ease in setting up, the direct control of joints and the speed of execution make the procedure adapt for optimization routines, as described in Section 6. Finally, some feasible applications are described in Section 7 studying an articulated PKMs with four dofs.

Outlines of the algorithm
Before introducing the reader to the proposed methodology, we have to clarify the reasons of this work. The first question that might arise is: Why use a new method without recurring to well confirmed formulations existing in the literature? The right answer is essentially tied to the simplicity of the proposed formulation. The method is addressed to designers that want to implement software to study the elastodynamics of PKMs, as well as to students and researchers that wish to create their own customized algorithm. The author experienced that the elastodynamics study of a complex robotic system is not a trivial issue, mainly when dealing with some features concerning the optimization of performances in terms of elasticity or admissible range of eigenfrequencies. Performing these analyzes often needs for a global optimization all over the workspace of a robot too cumbersome to be faced with conventional formulations. Further, considering that the elastodynamics behavior of a robot changes according to its pose, an accurate analysis should be carried out several times to capture the right response inside the robot's workspace, drastically affecting the computational cost. The focused issue becomes even more complicated when constrained optimization routines, based on indices of elastostatics/dynamics performances, are implemented to work in all the workspace. In the latter case the search for local minimum configurations needs for a simple and robust mathematical framework adapt to iterative procedures. Working with ordinary resources, in terms of CPU speed and memory, can make optimization prohibitive unless the complex flexible multibody formulations would be simplified to meet requirements. Therefore, the second question might be: Why simplify a complex formulation and do not create a simpler one instead? The sought formulation is what the author is going to explain in this chapter.
Let us start considering a generic PKM. It is essentially a complex robotic system with parallel kinematics in which one or more limbs connect a base platform to a moving platform, as shown in Fig. 1. The latter contains a tool, often referred to as end-effector, necessary to perform a certain task or, sometimes, as in the case of flight simulators or assembly stations, the end effector is the moving platform itself. The limbs connecting the two platforms are composed of links constrained by joints. A limb can be a serial kinematic chain or an articulated linkage with one or more closed kinematic chains or loops.
In performing our analysis we have to choose what is flexible and what is rigid as well. Generally, each body is flexible and the notion of rigid body is an abstraction that becomes a good approximation if strains of the structure are small when compared to displacements. Thus, considering a link flexible or rigid depends on many aspects as: material, geometry, wrenches involved in the process. Besides, some tasks might need for high precision to be accomplished, then an accurate analysis should take care of deformations to fulfill the requirements without gross errors. Here, we model the MP and BP as rigid bodies because, for the most part of industrial PKMs, these are usually one order of magnitude stiffer than the remaining links. The latter will be either flexible or rigid depending on the assumptions made by the designer. As already pointed out in the Introduction, we use 3D Euler beams to represent links, even if the treatment can be extended to superelements, as reported in (22). The formulation is linear and only small deformations are considered. Given a starting pose of the MP, the IKP allows us to find the robot's configuration; then, the elastodynamics-statics-is performed around this undeformed configuration. It might be useful to lock the actuated joints in order to avoid rigid movements and to isolate only the flexible modes. Below, we summarize all steps necessary to find the elastodynamics equations of a PKM: Individuate the couples of bodies linked by a joint and distinguish among three cases:  rigid-flexible, flexible-flexible and flexible-rigid. 6. Find the equations expressing dependent nodal-array in terms of independent nodal-arrays, then find the generalized stiffness and inertia matrices. The latter change whether lumped or distributed method is used.
7. Introduce the global array q of independent nodal coordinates. This array contains all the independent nodal arrays in the order defined by the reader.
8. Expand all matrices expressing each generalized stiffness and inertia matrix in terms of q by means of the Boolean matrices B 1 and B 2 . Then, sum all contributes to find the generalized stiffness matrix K PKM and inertia matrix M PKM of the PKM. 9.
Introduce the array f of generalized nodal wrenches and, finally, write the elastodynamics equations.

Mathematical background and key concepts
In a mechanical system flexible bodies storage and exchange energy like a tank is able to storage and supply water. The energy associated to deformation is termed strain energy while the aptitude of a body for deformation is tied to the property of stiffness. For a continuum body the stiffness is distributed and the strain energy changes according to the variation of the displacement field of its points. In a discrete flexible element inner points' displacements depend on displacements of some points or nodes.

Nodal and joint-array
Robotic links can be well described by means of 3D-beams, that is, flexible elements with two end-nodes, as shown in Fig. 2. The expression of the strain energy depends on the kind of displacements chosen for rotations: slopes, Euler angles, quaternions, and so on, and on the entity of displacement: large or small. Here, we choose a linear formulation based on small displacements and Euler angles to describe rotations. Based on these assumptions, the strain energy is a positive-definite quadratic form in the nodal displacement coordinates of the end-nodes arrays of the beam, where a nodal displacement array u = u x u y u z u ϕ u θ u ψ T has six scalar displacements, three translational and three rotational. Hereafter, subscripts and superscripts will be, respectively, referred to the beam and to one of the two end-nodes of a beam. Beams belonging to the same link are contiguous, while joints couple beams belonging to two adjoining links. To express the kinematic bond existing between the two nodes, or sections, coupled by the joint, a constraint equation is introduced in which the nodal displacements of the two nodes are tied through the means of an array of joint displacements. Figure 2 describes two beams linked by a prismatic joint of axis parallel to the unit vector e.T h et w o nodal arrays u 2 1 and u 1 2 are tied together by means of the translational displacement s of the prismatic joint P and by the 6-dimensional joint-array h P = e T 0 T T , i.e., We stress that eq.(1) describes a constraint among displacements, then deformations, and not nodal coordinates. In general, H j is a 6 × m(j) joint-matrix where m(j) is the dimension of the joint-array θ j . The joint-matrix H j and the joint-array θ j depend on the type of joint: the former containing unit vectors indicating geometric axes, the latter containing joint displacements, either linear s j , for translations, or angular θ j , for rotations. Below, two more examples of joints are provided: Revolute joint: where e is the unit vector along the axis of the revolute joint R and θ is the angular displacement about the said axis.
Universal joint: where e 1 and e 2 are the unit vectors along the axes of the universal joint U and θ 1 and θ 2 are the angular displacements about the axes of U.
Other joints, i.e. cylindrical, spherical and so on, can be created combining together elementary prismatic and/or revolute joints. The described constraint equations are used to consider joints contribute to elastodynamics in a direct way without recurring to Lagrangian multipliers to introduce joint constraints, (10). In the next section the above equation will be used to obtain joint displacements in terms of nodal displacements.

Strain energy and stiffness matrix
The strain energy , the strain energy assumes the following expression The 12 × 12 stiffness matrix K i of a 3D Euler beam with circular cross-section, and for the case of homogeneous and isotropic material, depends on geometrical and stiffness parameters as: cross section area A,lengthofthebeamL, torsional constant J,massmomentofinertiaI,the Young module E and shear module G. For our purposes, it is convenient to divide K i into blocks, i.e.
in which K 2,1 i =(K 1,2 i ) T and the other blocks are defined as where the two diagonal block matrices, respectively, refer to the nodal displacements of the two end-nodes while the extra-diagonal blocks refer to the coupling among nodal displacements of different nodes: in fact, each entry of the generic 6 × 6 block-matrix K l,m i can be thought as a force-torque-at the lth-node of the ith-beam when a unit displacement-rotation-is applied to the mth-node.

Rigid body displacement
Let us consider a rigid body with center of mass at the point G and a generic point P inside its volume, besides let d P be the vector pointing from G to the point P. If the body can accomplish only small displacements/rotations, let p be the displacement of point G and r,b et h ea x i a l vector of the small rotation matrix R, (9). Upon these assumptions, the displacement array u G of G is defined as u G = p T r T T . the displacement array u P can be expressed in terms of u G by means of the following equation: where 1 and O, respectively, are the 3 × 3 identity-and zero-matrices and D P is the Cross-Product Matrix (C.P.M.) of d P , (23).

Stiffness matrix determination
Let us consider a linkage of flexible beams and rigid bodies connected by means of joints. The strain energy V(u 1 1 , u 2 1 , u 1 2 ,...,u G1 , u G2 ,...,θ 1 , θ 2 ,...) of this system depends on nodal displacement arrays of both rigid and flexible parts and on joint displacement arrays. In the following, we will show how to express V in terms of a reduced set of independent nodal coordinates. In this process, we will start from three elementary blocks composed of rigid and flexible bodies that will be combined to build a generic linkage, for instance, the limb of a PKM. The first and the third case pertain the coupling between a rigid body and a flexible beam by means of a joint; the second case describes the coupling of two beams belonging to two different links. The choice to use two cases to describe the rigid body-flexible body connection is only due to convenience to follow the order of bodies from the base to the moving platform of a PKM. The reader might recur to a unique case to simplify the treatment. The last part of the section is devoted to some insight on joints' stiffness and feasible application to flexures.

Case a) Rigid body-flexible beam
The array u 1 2 can be expressed in terms of u 2 1 and θ recalling eq.(1), while u 2 1 , in turn, is tied to u 1 1 from eq. (7): upon combining both the equations, we obtain where G depends on d. By substituting the previous equation into eq.(8) we find that V a = V(u 1 1 , u 2 2 , θ). If the joint is passive, its displacement θ depends on the elastic properties of the system and, therefore, on the two displacements u 1 1 and u 2 2 : thus, it implies that V a is only function of the two mentioned array, i.e. V a = V(u 1 1 , u 2 2 ). In order to obtain the law for θ we minimize the strain energy V a w. r. t. θ: where 0 6 is the six-dimensional zero array. After rearrangements and simplifications, we find where Then, by substituting eq.(11) into eq.(9), we obtain Let us define the 12-dimensional arrayũ a = u 1 and let us substitute the above expression into V a , thereby obtaining where K a is the 12 × 12 stiffness matrix sought.

Case b) Flexible body-flexible body
For the case b) of Fig.4 two beams are coupled by a joint. The strain energy V b is function of the nodal displacement arrays of the two flexible bodies, i.e.
The four arrays are not all independent as the following equation stands: The strain energy is, thus, dependent on u 1 1 , u 1 2 , u 2 2 and θ: Here, we choose to minimize V b w. r. t. θ, as for the case a), and u 1 2 .T h e reader should notice that our choice is not unique, it is only a particular reduction process necessary for our purposes, it means that the reader might develop a treatment in which u 1 2 is not dependent anymore. Now, by imposing that the derivative of V b w. r. t. θ vanishes, we obtain Following the previous condition, even the dependent nodal-array u 1 ; therefore, by applying the chain-rule, we obtain The array u 1 2 is then written in terms of u 1 1 , u 2 2 ,i.e.
where the 6 × 6 matrices G 1 and G 2 have the following expressions: The joint-arrays θ j becomes, The dependent nodal-array u 2 1 is obtained by substituting eq.(21) and eq.(23) into eq. (17): where K b is the 12 × 12 stiffness matrix.

Case c) Flexible body-rigid body
The case c) describes a beam coupled to a rigid body by means of a joint, as shown in Figure 5.
The expressions involving the case c) can be easily obtained by extension of those of the case a), hence, we write only the final results for brevity: The strain energy V c associated to the beam becomes is the 12-dimensional array of independent nodal displacements and K c is the 12 × 12 generalized stiffness matrix of the case c).

Joint's stiffness and flexure joints
The final part of the present section is devoted to show some feasible extension of the formulation to flexure mechanisms. Similar concepts can be applied even to ordinary joints to take into account joint stiffness. In order to reproduce the counterparts of mechanical joints, in a continuous structure it is a common strategy to recur to flexure joints, in fact, zones where the geometry and shape are designed to increase the compliance along specified degrees of freedom (dofs). An ideal flexure joint should allow for only motions along its dofs, while withstanding to remaining motions along its degrees of constraint (docs). Figure 6 shows the cases of ideal and real flexure revolute joints. While for the ideal case the nodal displacements u 2 1 and u 1 2 of the two flexible beams (1) and (2), coupled by the joint, are tied by the usual constraint equation for the real case H ≡ 1 6 . The joint displacement array is now denoted with δθ. Referring to Figure 7, let us consider three different configurations of the flexure joints: an undeformed and unloaded configuration in which the joint-array is θ o ; a preloaded initial configuration with the joint-array θ i and a final configuration in which θ f = θ i + δθ,b e i n g δθ an array of small displacements around the initial joint-array θ i . For the following explanation, we refer only to the ideal case of Figure 6. Let K θ be the m(j) × m(j) joint stiffness matrix and let w be the generic wrench acting on the flexure joint; then, for the three configurations, we can write where Δθ f = Δθ i + δθ. The strain energy V θ f of the flexure joint in its final configuration simply reduces to The above expressions can be simplified considering θ o = 0 6 for the unloaded configuration.
We do not further discuss on flexure joints, leaving the reader to derive three new cases, similar to those discussed above, taking into account joint stiffness.

Mass matrix determination
In this section two ways to include masses/inertias are discussed: the lumped approach and the distributed approach. The former concentrates masses on nodes of both rigid and flexible parts; the latter considers the real distribution of masses inside beams. As will be explained in the text, the two methods produce good results, particularly when an accurate degree of partitioning is chosen for flexible bodies. Finally, we focus our attention on a way to consider joints with mass and inertia.

Lumped approach
Reducing mass and inertia of a rigid body to a particular point, the center of mass, without changing dynamic properties of the system, is a common procedure. On the contrary, for flexible bodies every reduction process is an approximation generating mistakes. Let us refer to  and inertia of each beam is concentrated on its end-nodes. Here, we choose a symmetric distribution in which the mass m is divided by two, but the reader can use another distribution according to the case to be examined, as instance beams with varying cross section. The first node has a mass of m/2 while any other node, but the fourth, bears a mass m because it receives half a mass from each of the two contiguous beams coupled at its section. The fourth beam of the link is attached to a joint. According to what explained in the previous section, even the mass is concentrated only on the independent joints: it means that the mass of the last beam has to be concentrated on the fourth node, thereby the latter carrying a mass equal to 1/2m + m = 3/2m. Analogous arguments, not reported here for conciseness, may be used to describe inertias. The lumped approach concentrates masses and inertias on all the independent nodes, belonging to both rigid and flexible bodies. Mathematically, a 6 × 6massdyad M i is associated at the ith-independent node, defined as where m i and J i are the mass and the 3 × 3 inertia matrix, respectively, of either a rigid body, if the independent node is at the center of mass of the said body, or of a beam's end-node. The generalized inertia matrix M PKM of a PKM is readily derived upon assembling in diagonal blocks the previous inertia dyads, following the order chosen to enumerate all independent nodes inside the global array of independent nodal coordinates, hence

Distributed approach
The distributed approach considers the true distribution of mass inside a flexible beam. Particularly, it is possible to find a 12 × 12 matrix M i referred to the twelve nodal coordinates of a beam's end-nodes. This matrix can be divided into blocks, whose entries are reported below, as already done for the stiffness matrix K i ,i.e.
where ρ, L i , A i , I x , I y and I z , respectively, are the density, the length, the cross section area and the mass moments of inertia for unit of length of the ith-beam. The matrix M i is associated to the kinetic energy T i of a beam, the latter being defined as a quadratic forms into the time-derivativesu 1 i ,u 2 i of the nodal displacement arrays u 1 i and u 2 i : In the previous section we have found how dependent nodal displacement arrays are expressed in terms of independent ones. Let us consider the time-derivative of both sides of eq.(13), we writeu with obvious meaning of all terms. The three mass matrices, above defined, are referred to the time-derivativesũ a ,ũ b andũ c of the independent nodal displacement arraysũ a ,ũ b andũ c .T o clarify some doubt that might arise let us refer to Fig. 3. In this case the mass matrix of the rigid body, see eq.(35), is referred to its center of mass with displacement array u 1 1 .T h ed i s t r i b u t e d approach first allows us to find the 12 × 12 mass matrix M 2 of the beam, then the latter is expressed in terms of only independent displacements by means of M a . Obviously, M 2 and M a refer to the same object, the beam; but M a distributes M 2 into the independent nodes with displacement arrays u 1 1 and u 2 2 . This implies that the center of mass of the rigid body carries its own mass/inertia and part of the mass/inertia of the beam. Similar conclusions may be sought for the cases b) and c).

Joint's mass and inertia
In this final subsection we describe a method to consider mass and inertia of joints. For convenience, let us refer to Fig. 4. The joint is split into two parts, one belonging to the first beam, the remaining to the second beam. Let M L and M R , where capital letters stand for left and right, be the mass dyads of the two half-parts of the joint, defined as The kinetic energy T J of the joint can be written in terms of the above matrices M L and M R and of the time-derivatives of the nodal arraysu 2 1 ,u 1 2 ,thus where M J is the 12 × 12 generalized inertia matrix of the joint expressed in terms of independent nodal displacements, respectively, defined as for the three cases in exam.

Linearized elastodynamics equations
In this section we will derive the generalized inertia and stiffness matrices necessary to write the linearized elastodynamics equations. As described in the two previous sections stiffness and inertia matrices of rigid bodies and flexible beams are referred to the independent nodal displacements of the system. Now, one might ask how to combine different matrices to build the global stiffness and inertia matrices. In order to solve this issue let us introduce two Boolean matrices B 1 and B 2 able to map a 6-dimensional displacement array u i or a 12-dimensional displacement arrayũ in terms of a 6n-dimensional global array q = u 1 u 2 ... u n containing all the independent displacement arrays of the robot. It is important that the reader would define the order in which every single array u i appears inside q.T h e6× 6n matrix B 1 and the 12 × 6n matrix B 2 are defined as: The above expressions allow us to convert each nodal displacement array in term of q.A s instance, let us take into exam the strain energy V b of eq.(27), it simply turns into w h e r ei and j indicate the position indices of u 1 1 and u 2 2 inside q. From V b we can extract a new 6n × 6n matrix K b , that is a stiffness matrix expanded to the final dimension of the problem, defined as Likewise, the generic expanded 6n × 6n mass matrix M i of M i , see eq.(35), can be written as of nodal displacement coordinates, these matrices can be summed and combined to obtain the global generalized matrices K PKM and M PKM of the robot. The way to use B 1 and B 2 will be shown in detail in the case study of the next section. Let us introduce the 6n-dimensional array f of generalized nodal wrenches, i.e. nodal forces and torques, then, the system of linearized elastodynamics equations becomes The previous system may be used to solve statics around a starting posture of the robot. The homogeneous part of eq.(50), i.e. the left side, may be used to find eigenfrequencies and eigenmodes of the robot, or the zero-input response. As already pointed out, the natural frequencies of the robot change with regard to the pose that the robot is attaining. To determine how stiffness or natural frequencies vary inside the workspace local and global performance indices may be introduced to investigate elastodynamic behavior of PKMs. As instance, let us consider the first natural frequency f 1 mapped inside a robot's workspace. The latter can be analyzed to differentiate areas with high range of frequency from areas in which the elastodynamic performances worsen. The multidimensional integral Ω 1 of f 1 , extended to the whole workspace W or part of it, can be a good global index to show how global performances increase or decrease changing some geometrical, structural or inertial parameter, i.e.
where, in the approximated discrete formula of the right side, W is divided into n w hypercubes W i , while the frequency f 1i is calculated at the center point of W i . We conclude this section considering an useful application of the method to find singularities. It is well known that as a robot approaches to a singularity configuration it gains mobility. As a consequence, its generalized stiffness matrix becomes singular and its first eigenfrequency goes to zero. Analyzing the variation of f 1 i n s i d et h ew o r k s p a c em a yb eu s e f u lt oa v o i d zones close to singularity or with low values of frequency. We stress that singularities directly come from kinematics: stiffness and inertial parameters, respectively influencing the elasticity and the dynamics of a mechanical system, can only amplify or reduce, without creating or nullifying, the effects of the kinematics. Different considerations can be made for the boundary surface of the workspace. For the latter case, a robot undergoes poses in which one or more legs are completely extended, thus the robot loses dofs associated to those directions normal to the tangent plane at the boundary surface of the workspace. If the displacement of the MP for the corresponding eigenmode, when evaluated at a pose lying on the boundary surface, is normal to the said tangent plane the robot is virtually locked, it means that its stiffness along that direction is high, thus reflecting into an increment of frequency. On the contrary, if the displacement of the MP lies onto the tangent plane its in-plane stiffness, and its natural frequency as well, in theory goes to zero. All remaining cases are included between zero and a maximum value depending on the projection of displacements on the tangent plane. This important conclusion can be summarized by the following statement: the value of a given natural frequency of a PKM on the boundary surface of its workspace can be interpreted through the projection of the moving platform's displacements, for the corresponding eigenmode, on the tangent plane to the said surface. The previous statement is strict only when a three dimensional workspace, e.g. the constant orientation workspace, can be considered. In other cases, one should recur to concepts of projection based on twist theory, beyond the scopes of this chapter.

Case study
In this section we apply the method to a complex PKM with 4-dofs. The robot, shown in Fig. 9, is composed of four legs connecting the base frame to a moving platform. Due to the particular kind of constraints generated by the limbs, the moving platform can accomplish three translations and a rotation around the vertical axis. This motion, even referred to as Schönflies motion, is typical of SCARA robots and it is largely used in industry for pick and place applications. Each leg is composed of a distal link connected to the base frame by means of an actuated revolute joint and to a parallelogram linkage, i.e. a Π-joint, by means of a revolute joint. The parallelogram, in turn, is linked to a proximal link by another revolute. Finally, the proximal link is coupled to the MP by a vertical axis revolute joint. All revolute joints, apart from the one fixed to the inertial frame, are passive. We do not further go inside the analysis of the kinematics, citing (24) for more explanations and for any detail pertaining the robot's geometry.

Application of the method to a generic leg
In order to apply the proposed method each leg is decomposed into three modules. The first module includes the BP and the proximal link, the second the parallelogram linkage, while the last includes the distal link and the MP. We recur to some simplifications that do not change the meaning of the treatment. In general, the two bases of a PKM are one order stiffer than the links composing the structure and can be modeled as rigid without loss of accuracy. Even the short input and output links of a Π-joint can be considered rigid when compared to the slender coupler links. As verified by FEM, the said assumptions are good approximations that simplify the final model introducing rigid-flexible and flexible-rigid connections inside and between each module. Let us analyze the first module shown in Fig. 10. The distal link is split into three flexible bodies, the choice of discretization being absolutely arbitrary. In this simple case the end-bodies F 1 and F 3 of the link are coupled to two rigid bodies, i.e., the base BP and the input link I of the Π-joint, by means of two revolute joints of axes parallel to the unit vector . The first actuated revolute joint is considered locked to perform the elastodynamics analysis, thereby, u 1 1 ≡ u 1 BP ≡ u 1 . In this way the generalized stiffness matrix of the robot is not singular as rigid body motions of the robot are deleted and the analysis is performed at a frozen configuration. Even for the other three fixed connections at points 2 and 3, similar equations stand: u 2 1 ≡ u 2 2 ≡ u 2 and u 3 2 ≡ u 3 3 ≡ u 3 . Finally, for the displacement array of the input link R,wecanwrite:u R ≡ u 4 . By substituting into the previous expressions, we deriveũ 1 . The local stiffness matrix K i of the ith-body F i is expressed into the global reference frame. In Figure 10 the first local frame {x 1 , y 1 , z 1 }, with origin at point 1, and the base global reference frame {x, y, z}, with origin at point O, are displayed. The set of independent displacements are defined inside the independent displacement array q M1 of the first module: q M1 = u 1 u 2 u 3 u 4 ,w h e r eu 4 is the nodal displacement array of the center of mass of the rigid body I. While the first and the second term, i.e. V F 1 and V F 2 ,ofthe strain energy V m1 contain only independent displacements, the last term V F 3 is function of the dependent nodal array u 4 3 .T h efl e x i b l eb o d yF 3 , indeed, is coupled by a revolute joint to the rigid body I, thus, the case c), flexible-rigid, can be applied to express its strain energy only in terms of independent displacements. Following the notation above introduced, we write: whereũ c3 ≡ u 3 T u 4 T T and K c3 has been defined in eq.(31). Notice that in order to find K c3 we have to use the 6-dimensional joint array h R (e); besides, the position vector d going Then, by introducing the 12 × n qm1 binary-entry matrix B 2 (•, •) mapping the independent nodal-arrays in terms of the array q M1 , the expression of the generalized stiffness matrix K M1 can be calculated as Here, the matrix K M1 is expressed in terms of q M1 . The reader must notice that to obtain K PKM each matrix has to be expressed in terms of a final n q -dimensional array q including all the independent nodal displacements. In order to define q, as the modules of a leg are in series, the last node of a module coincides with the first one of the next module and it must be denoted by a unique nodal displacement array. Extension to the final dimension of q is readily obtained by substituting into eq.(54) a new 12 × n q binary-entry matrix B 2 (•, •) mapping the independent nodal-arrays in terms of the global array q.
The second module is shown in Fig. 11(a). In this case two rigid-flexible and two flexible-rigid connections must be used to describe the couplings between flexible couplers and input-output rigid bodies I and O. The joint array h R (e) remains the same for all the four revolute joints as a matter of fact of the parallelogram architecture. The third module, shown in Fig. 11(b), includes one rigid-flexible and one flexible-rigid connection, respectively, between the distal link and the output link O of the Π-joint and the distal link and the MP.
As displayed in the same figure, the two revolute joints to be used to find the joint arrays h R have axes of unit vectors e 1  stiffness matrices of the two modules as can be readily obtained following procedures similar to the case of the first module.

Comparison to FEA software and validation
In this subsection the proposed elastodynamics model is compared to FEA results for validation. Two FEA models, with increasing complexity, are implemented in order to establish the nearness of the examined case to the real case. The first model, developed by the commercial software Marc © , considers only 3D-beams and rigid bodies, while joints are modeled by means of relative degrees of freedom between common nodes belonging to two coupled bodies. This model perfectly fits the simplifications of the method in exam. The second model, developed by the commercial software Ansys © , describes a robot with a complex structure closer to reality. Links are solids while joints have finite burdens and provide their function by means of the coupling of surfaces and screws.

Statics
We first compared the static deformation of the robot when an external force of 1000 N, along the x-direction, is applied on the end-effector, when the latter is positioned at its home pose with angle of rotation of MP equal to zero. Figure 11 shows the displacements of our model when compared to Ansys © model, while in Tab. 1 the translational displacements of the end-effector center node are reported. It can be observed how the first MARC © model perfectly fit to the results of our method. The relative error grows, but still remaining limited, when displacements are compared to Ansys © model. The reason is well understood as it comes from the use of solid bodies and joints to simulate the robot's structure.

Natural modes and frequencies
As in the case of statics, the elastodynamics model of the 3T1R robot is used to compare natural modes and frequencies to FEA software. We report only the comparison to Ansys © ,asmore indicative of a feasible application of the method to a real case. Figure 12 shows the first four natural modes obtained when the robot is attaining its home pose. In Table 2  with Ansys © model. Some discrepancies still occur due to the same reasons explained for the static case and to the use of flexible bodies to simulate all links and platforms in Ansys © . This choice has been taken in order to avoid asymmetric contacts between rigid and flexible solids leading to convergence mistakes. In turn, we used a stiffer material to approximate rigid behavior of the MP and rigid links.

Distribution of frequencies inside the workspace
In order to provide a further useful extension of the proposed method, the first natural frequency is calculated and plotted inside the constant orientation workspace of the robot (25). We have chosen elementary cubes of side 5 cm to discretize the workspace of the robot while the angle of rotation of the MP is θ z = 0. It can be observed that two privileged diagonals divide the workspace into four symmetric areas. These directions coincide with the two diagonals of the squared BP of Fig. 9, meaning that, at the home pose, geometric symmetry reflects itself into the elastic behavior of the robot. The boundary of the constant orientation workspace shows areas with high range of frequency along with other areas in which the first natural frequency reaches values close to zero (Hz). As already outlined in the text, this behavior is due to the MP's displacement (deformation) in correspondence of the first eigenmode: if at a certain pose, in which the analysis is performed, the displacement is along a doc of the MP, the ensuing frequency will be high; conversely, if it is along a dof of the MP, the frequency will come near zero.
[Hz] Fig. 13. Distribution of the first natural frequency inside the robot's constant orientation workspace.

Conclusions
This chapter has discussed a method, based on the Matrix Structural Analysis, to study the linearized elastodynamics of PKMs. Base and moving platforms are considered rigid, while links can be modeled as rigid or flexible parts, the latter being decomposed into two or more flexible bodies. Here, we used 3D Euler beams but the method can be extended to superelements with two end-nodes. Joints are directly included, without recurring to Lagrange multipliers, by means of kinematic constraints between nodal displacement arrays. Three cases have been taken into account to model the rigid-flexible, flexible-flexible and flexible-rigid coupling of bodies by means joints. Each case yields equations, linking dependent, independent nodal displacement arrays and joint displacements as well, to be used to find generalized stiffness and inertia matrices. The latter are then combined as elementary blocks to find the global matrices of the whole system. Some useful extension to compliant mechanism has been introduced, while two strategy, the lumped and the distributed one, have been explained to include mass/inertia into the model. Feasible applications of the method pertain: the study of natural frequencies inside the robot's workspace by means of local and global indices, the singularity finding, the optimization of elastodynamic performances varying geometric, structural or inertial parameters. Finally, the method has been applied to an articulated four-dofs PKMs with Schönflies motions. A modular approach is used to split each of the four legs into three modules.
Results, compared to commercial software, revealed good accuracy in determining natural frequency range, while drastically reducing the time of computation avoiding the annoying and time-consuming FEM meshing routines. The robotics is an important part of modern engineering and is related to a group of branches such as electric & electronics, computer, mathematics and mechanism design. The interest in robotics has been steadily increasing during the last decades. This concern has directly impacted the development of the novel theoretical research areas and products. This new book provides information about fundamental topics of serial and parallel manipulators such as kinematics & dynamics modeling, optimization, control algorithms and design strategies. I would like to thank all authors who have contributed the book chapters with their valuable novel ideas and current developments.