Vibration Analysis of MultilinkManipulators Based on Timoshenko Beam Theory

Timoshenko’s theory is adopted in order to accurately describe the freely vibrating dynamics of a multilink flexible manipulator. It is herein presented an analytical modelling strategy that extends previous works through a more refined model which accounts for elastic complicating effects along with lumped inertial loads which are typically mounted on joints of manipulators; in this regard, more accurate results are provided. The eigenproblem is presented from an analytical point of view through a matrix formulation, thus providing an essentially closed formula. Apart from the limitations of the implementing calculator, the formulation can take into account an arbitrary number of links in an arbitrary settled configuration, thus allowing relevant analytical analysis and avoiding the need to recur to nonimmediate numerical schemes. Once the analytical model is introduced, solutions are compared to both those achieved by previous models and those obtained by a finite elements method.


Introduction
Within the frame of robotic applications, high-performance research has recently pushed designers to realize lighter multi-link manipulators; these requirements, indeed, ensure the attainment of high velocities by keeping energy consumptions low and smaller actuators/brakes.The result is a more manageable system which, however, presents nonnegligible problems related to its higher deformability along with large displacements.The high deformability is mainly responsible for induced vibrations during the operations, and, therefore, the need to implement a suitable controller, aimed at significantly reducing the vibrations, becomes a logical consequence.Moreover, when large displacements are contextually involved within the mathematical model, further complications arise, for example, nonlinear straindisplacement relations, noninertial frame effects, timevarying boundary conditions.
In recent decades, many authors have dealt with the above-mentioned problem, and, in this regard, a certain number of models, along with currently open mathematical questions, have been identified.A perusal of these references shows two main problems which the researchers are faced with: (i) the need to implement the dynamic analysis of a local vibrating behaviour and (ii) the capability of extending such a local vibratory analysis to domains larger than those close to an established geometrical configuration.The former problem is essentially related to the description of vibrating linear systems around their stable and static configuration, whilst the latter problem involves the description of large displacements.It is in this mentioned context that a controlling strategy must be implemented.
Several researchers have, in the past, analysed robotic manipulators with deformable links in order to investigate the related dynamics and/or controllers.Significant contributions can be attributed to [1,2], which are related to single link manipulators carrying a payload at one end; in particular, Bellezza et al. [1] deal with an Euler-Bernoulli beam model whilst White and Heppler [2] consider the case of one only flexible slewing beam based on Timoshenko's theory.Recently, Chaolan et al. [3] analyzed only one flexible hub-beam system by accounting for the influence of shear and axial deformation and thus by covering the complexity of dealing with large deformations of the particular system investigated (i.e., one hub-beam).
Ower and Van de Vegte [4] based on an Euler-Bernoulli beam model, along with a Lagrangian approach, presented a flexible manipulator model; their analysis did not involve large displacements, rather the vibration is investigated in the neighbourhood of an established configuration; in this regard the problem is linearized and, therefore, a linear transfer function matrix has a sense in the discussed frequency domain.
Tomei and Tornambe [5], again on the base of an Euler-Bernoulli beam, developed an approximate model through a Lagrangian approach and solved the equations by using the Ritz method.The authors also reported the set of explicit equations for one-and two-link models only.
De Luca and Siciliano [6] used a similar order of approximation to that of Tomei and Tornambe [5] and, in addition, assumed simplifying end mass-boundary conditions in order to settle the set of governing equations.In spite of the promise to investigate large displacements, the equations of the model of De Luca and Siciliano [6] were used to analyze the vibratory behaviour of only two links in the neighbourhood of an established configuration; this was mainly caused by the fact that the mode shapes became time-varying for the mass boundary conditions, indeed, time varying.
Chen [7] established the generalized dynamic model for a planar n-link flexible manipulator basing the relevant analysis on a similar order of approximation of the previous two discussed references.
Lee [8] mentions a certain incompatibility in dealing with the modelling of flexible robots with the bending mechanism through the conventional Lagrangian approach.In particular, Lee [8] raises a question regarding certain free link elongations and the need to include this influence in order to reduce modelling inaccuracy, which was the case of previous works.Therefore, in the hypothesis of an Euler-Bernoulli beam, Lee [8] proposes a new link deflection model to fix the mentioned incompatibilities.Moreover, Zhang et al. [9] pointed out that the assumed-mode method leads to inaccurate models for flexible links because mode shapes change continuously as the solution proceeds in time (i.e., as mentioned beforehand, large displacements involve nonlinear time-varying systems); in this regard Zhang et al. [9] stressed the need to develop an accurate model in order to investigate and design reliable control techniques, and, therefore, obtained a system of partial differential equations for a two-link Euler-Bernoulli manipulator by using the principle of Hamilton.
In summary, the analysis of the relevant literature essentially shows models of manipulators developed on the base of Euler-Bernoulli beams; a part of these references introduce nonsharable approximations and/or solve the relevant equations by using the assumed-mode method within the frame of time-varying systems.Within the frame of the difficulties related to the mathematical formulation of the problem herein being dealt with, Pascal [10] could also deserve further attention.In any case, a full mathematical description/solution of (i) local vibrations (ii) superimposed on large displacements of flexible multi-link manipulators cannot be considered a closed question yet neither in its single aspects (i) or (ii) nor on the whole (i) and (ii).
This work aims at introducing a mathematical model which is able to analytically describe the freely vibrating dynamics of a robotic manipulator assembled through multiple links around static stable configurations established a priori.In this regard, apart from the model approximations related to the classical (Euler-Bernoulli) or uniform shear deformable theory (Timoshenko) in linear elasticity; the governing equations do not contain any approximation.The mathematical formulation in this paper is introduced in such a manner that a general configuration for multilink can be easily settled by recursively assembling blocks in matrix equations.The model herein dealt with must be considered superior to the previously published models at least for being able to deal with complicating effects such as uniform transversal shear effects along with concentrated masses on the joints.To the best of the authors' knowledge, only the works by Milford and Asokanthan [11] and di Castri et al. [12] could be considered of comparable accuracy and, in spite of this conviction, it will be herein shown, by numerical comparisons, as the model proposed in this paper can either contain the previous two works [11,12] as particular cases, and can recover large discrepancies occurring from the model presented in [11].
It is stressed that the modal characterization of the system presented by Milford and Asokanthan [11] derived from partial differential equations governing the free vibrations of a two-link flexible manipulator only and neglected the influence of the beam-by-beam axial dynamics.In this paper the modal analysis of a multi-link flexible manipulator is presented by using the Timoshenko beam theory.Moreover, an analytical treatment of the boundary conditions and the corresponding formulation also includes the local axial dynamics of the structure.The axial dynamics could be neglected at a single level beam (as a minor and decoupled influence), but it should be included in the general configuration of a robotic manipulator due to its global coupling with the flexural behaviour.Therefore, the present model, being able to efficiently and accurately extract eigendata, aims at founding bases for simulating the dynamic behavior (free and forced) of multilink flexible manipulators, thus also allowing to design control strategies for the same structures.Generally speaking, the author is required to choose between Timoshenko and Euler-Bernoulli formulations, depending on the desired accuracy on predictions but keeping in mind that axial contributions should never be neglected because they are coupled with the flexural ones.
The modelling procedure is based on a modular mathematical description; each link corresponds to a block in the matrix formulation of the system, so that the addition or the elimination of that link is equivalent to adding or deleting its corresponding block.This new introduced formulation enables the extraction of exact modal data without the need to use numerical methods as, for example, the finite element method which also requires convergence tests.Flexibilities are assumed to be distributed along the links only, whilst joints are treated as rigid; the effect of this assumption on modal data is discussed through several simulations.
A final numerical example is carried out by referring to the structural parameters of a typical industrial manipulator with three links used for the handling of loads.Simulation results point out how some differences can exist between the present model and the simplified ones.In order to validate the model herein presented, for all the analyzed cases the calculated frequencies are compared with those achieved through a finite element package.Peculiar mode shapes are also accurately depicted in order to offer a fruitful comparison with respect to future results coming from more accurate or approximate models.
The paper is organized as follows.In Section 2 the multi-link manipulator model, the geometrical and mass properties of the links, and the adopted reference systems are presented; moreover, the partial differential equations describing the system, the boundary conditions, and the resolving technique used to calculate modal data in detail are derived and discussed.Section 3 reports tables and comments about the comparison between the proposed model and those in [11,12], while Section 4 is dedicated to the industrial manipulator with three links by providing analytical and numerical results, and graphical representations of the natural modes.Finally, in Section 5 conclusions are drawn.

Problem Formulation.
Based on the enumeration of n links of a flexible multibody system (Figure 1), the link-by-link nomenclature is represented in Figure 2 with respect to its undeformed configuration.The manipulator is restricted in the horizontal plane.The links are of homogenous, isotropic and linearly elastic material with constant cross-sections; they have length L i , cross-sectional area S i , volumetric density ρ i , Young's modulus of elasticity E i , cross-sectional moment of inertia I i , shear modulus G i and shear correction factor k i , for i = 1, . . ., n, where n is the number of links.Each link presents a rigid body of mass M i and moment of inertia J i at its tip that for i = 1, . . ., n − 1 accounts for the presence of possible joint motors; i = n can represent a payload.The effects deriving from the geometrical offset between the joint axis and the corresponding link attachment are neglected in the mathematical formulation.The known uniform shear deformable assumptions are taken into account within the small deformation of linear elasticity.
The derivation of the modal data for each configuration of the manipulator is dealt with by blocking the joints and solving the corresponding differential eigenvalue problem.In particular, the obtained results can depend on the relative orientations between links; therefore, without loss of generality, it will be assumed that the joint variable at point O 1 (we say θ 1 ) always has null value.For a chosen configuration, the manipulator is fixed at point O 1 and internal blocked joints at points O i avoid relative rotations between links.
The analyzed multibody system.
The coordinate system (O 0 , X 0 , Y 0 , Z 0 ) is the global one whilst the coordinate system (O i , X i , Y i , Z i ) is associated to link i.The origin of frame (O i , X i , Y i , Z i ) is at point O i and the X i -axis coincides with the undeformed axis of the ith link.The angle between the X 1 and X 0 axes, that is, θ 1 , is always null for the above-mentioned reasons, so that frame (O 0 , X 0 , Y 0 , Z 0 ) and frame (O 1 , X 1 , Y 1 , Z 1 ) will be coincident; the angle between the X i and X i−1 axes is denoted with θ i , and the set {θ i : i = 2, . . ., n} will define the configuration of interest with respect to which natural frequencies and mode shapes will be obtained.
In the assumption of small deformations, the coordinates of a generic point P i of the ith link, with respect to the corre- T , where the two functions u i (x i , t) and v i (x i , t), respectively, represent the axial and transverse displacements undergone by the considered point.

Mathematical Modelling.
In the analysis of the ith link, let ψ i (x i , t) be the slope of the deflection curve when the shearing force is neglected and β i (x i , t) the angle of shear in the same cross-section; then the total slope of the deformed axis is given by [13] where v i (x i , t) is the transverse deflection of the ith link and ψ i (x i , t) is the cross-section rotation produced by the the bending component of deflection.According to the Timoshenko beam theory, the stress-displacement relations are given by the following equations: where u i (x i , t) is the axial displacement function, M i (x i , t) is the bending moment of the ith Timoshenko link, and T i (x i , t) and N i (x i , t) are the transverse shear force and the axial force of the link, respectively.Considering the link subjected to the transverse distributed load q i (x i , t) and to the axial distributed load p i (x i , t), the dynamic equilibrium equations of stress are given by ( When substituting expressions (2) in these latter equations, we get a system of partial differential equations the unknown axial displacement u i (x i , t), for the unknown deflection v i (x i , t), and the unknown rotation ψ i (x i , t) describing the dynamics of the ith link: In particular, (4) and ( 5) are coupled through the dynamics of the Timoshenko beam theory and describes two coupled modes of deformation: one mode of deformation is the transverse deflection of the link as measured by v i (x i , t) and the other mode is the transverse shearing deformation β i (x i , t), as indirectly measured by (1) and the bending rotation ψ i (x i , t).Finally, (6) is the mode of axial deformation of the link, as measured by u i (x i , t).
For each link of the manipulator there will be a system of partial differential equations describing the deformations of the link; if we consider n links, 3 × n of such equations will describe the overall system dynamics.By formulating the exact boundary conditions, in the case of free vibration, it will be possible to achieve the system differential eigenvalue problem associated to a generic configuration, whose eigenvalues and eigenfunctions will represent the sought modal data.

Modal Analysis.
The free vibration problem of the system is now considered; for such a case, the distributed forces q i (x i , t) and p i (x i , t) are zero, and, in this regard, the dynamics of the ith link is described by the following partial differential equations: Coupled ( 7) and ( 8) can be reduced to a single equation both for v i (x i , t) and ψ i (x i , t), having [14] for the transverse deflection v i (x i , t) and for the rotation ψ i (x i , t).

Journal of Robotics 5
Now, assuming separation of variables in the form with natural frequency ω, (11), (10), and ( 9) can be, respectively, expressed as ordinary differential equations in the following form: with unknown parameter ω and unknown functions Now, given the above 3 × n relations, we need 6 × n boundary conditions to completely define the differential eigenvalue problem corresponding to an arbitrary configuration of the multi-link manipulator; in fact, as will be seen later, such a number of boundary conditions is equal to the number of coefficients characterizing the spatial solutions of ( 13), (14), and (15).
We observe that, for a chosen configuration, the boundary conditions are formulated by writing relations between forces and moments in correspondence to each joint section; since the equilibrium conditions at inner blocked joints O i depend on the joint variables θ i , it results that boundary conditions are configuration dependent.
Figure 2 shows all the specified sections of the manipulator for the formulation of the boundary conditions; they correspond to the fixed joint at O 1 , where three equations will be written, the inner joints O i (i = 1, 2, . . ., n) between link i and link i − 1, where we will have six equations for each joint, and the free end of the structure, correspondent to the end effecter, for which three other equations will be derived.Such equations are herein immediately explicated: (i) O 1 section: (iii) EE section (end-effector) Using solutions ( 12) into ( 16)-( 21), the boundary conditions for the analyzed reference configuration of the manipulator reduce to (i) O 1 section: Figure 3: Rectangular hollow section of the adopted industrial manipulator.
(ii) O i section (i = 2, . . ., n): (iii) EE section (end-effector): Equations ( 13)- (15), expressed for all the n links of the structure, and the boundary conditions ( 22)-( 26) define the differential eigenvalue problem for the multi-link manipulator in the chosen nominal configuration; by determining the constant ω such that the problem admits nontrivial solutions , and U i (x i ) satisfying the boundary conditions, it is possible to obtain natural frequencies and mode shapes of the system for an assigned set {θ i : i = 2, . . ., n}.To this end, let us consider (14) in the following form: From the characteristic equation associated to (27) it is possible to calculate the following two solutions for each ω frequency: where Ω is a positive or zero quantity defined by It is observed that for ω = [(k i G i S i )/(ρ i I i )] 1/2 , solution (28) is zero; this limit value of ω, denoted by ω lim , delimits two different vibratory behaviors.In fact, for ω < ω lim , (28) and (29), respectively, give two real solutions ±λ i1 and two imaginary solutions ± jλ i2 , whereas, for ω > ω lim , four imaginary solutions ± jλ i1 and ± jλ i2 are obtained.As can be easily verified, the frequency spectrum represented by ω < ω lim is very broad if it is calculated with typical parameters of numerous flexible manipulators; this allows comprising many cases of real interest.For this reason, the resolving procedure for the derivation of the modal data is explicitly reported in the following, referring to the first of the above "frequency behaviour conditions"; nevertheless, it is remarked that for the case ω > ω lim , natural frequencies and mode shapes should be derived in a very similar manner, by only changing the general solutions of ( 13) and (14).Now, for ω < ω lim , ( 13)-( 15) have the solutions where c i = (E i /ρ i ) 1/2 is the well-known longitudinal wave speed.
The coefficients A ik and A ik (k = 1, . . ., 4) are not independent, but related through (8); in fact, substituting solutions (12) into it for v i (x i , t) and ψ i (x i , t) and using expressions (31) and (32), we obtain the following relations: where Moreover, inserting (34) into (31), the final form for the solution Ψ i (x i ) is given by Now, the expressions for Ψ i (x i ), V i (x i ), and U i (x i ) and the boundary conditions can be used to derive the frequency equation of the manipulator; indeed, substituting (32), (33), and (36) into ( 22)-( 26) we obtain 6 × n equations with 6 × n unknown coefficients A ik and B iz (i = 1, . . ., n, z = 1, 2, and k = 1, . . ., 4).These equations can be expressed in matrix form as is the 6 × n vector of modal coefficients characterizing the eigenfunctions of the assigned differential eigenvalue problem, whereas is a square matrix of order 6 × n and represents the system characteristic matrix, associated with the chosen nominal configuration of the multi-link manipulator.Such a matrix, whose elements are functions of the natural frequency ω, is a block matrix partitioned as in (39).In particular, M O1 (ω)   M Oi.i−1 (ω) and M Oi.i (ω) are 6 × 6 blocks and can be interpreted as a matrix representation of the six boundary conditions at the inner blocked joint O i between link i-1 and link i, concerning the twelve modal coefficients (A s1 , A s2 , A s3 , A s4 , B s1 , B s2 ), where s = i − 1, i.The blocks of M(ω) and all their coefficients are explicitly reported in Appendix A.
The homogeneous system (37) admits nontrivial solutions for the unknown coefficients A ik and B iz , (i = 1, . . ., n, z = 1, 2, and k = 1, . . ., 4) once the determinant of the characteristic matrix M(ω) is settled to be zero: Equation ( 40) is a transcendental equation whose algebraic expression can be obtained through a symbolic code, but it requires a root-finding algorithm to find its denumerable infinite number of roots ω j ( j = 1, 2, . ..), recognized as natural frequencies.Corresponding to each of these roots there are the eigenfunctions V i( j) (x i ) and U i( j) (x i ) that represent the jth axial-transverse mode shape, and the eigenfunction Ψ i( j) (x i ) that represents the jth cross-section rotation mode shape.Based on relations (32), (33), and (36), the required modal coefficients A ik( j) and B iz( j) for the abovementioned eigenfunctions are obtained by substituting the natural frequency ω j into the characteristic matrix and solving the homogeneous system of (37).
Referring to the axial and transverse local components of deformation V i( j) (x i ) and U i( j) (x i ), let us suppose to represent the deformed shape of the multi-link manipulator in the global coordinate frame (O 0 , X 0 , Y 0 , Z 0 ); for this scope, it is primarily observed that the homogeneous transformation matrix (43) settles the coordinate transformation between frames where R i−1 i is the rotation matrix describing the orientation of frame i relative to frame i − 1 and has the following form: T denotes the position of the origin O i of coordinate frame i relative to coordinate frame i − 1; by definition, we set L 0 equal to zero.Now, the position of any point of the ith deformed link is expressed in the coordinate frame (O i , X i , Y i , Z i ) through the homogeneous position vector Finally, such a point is directly expressed in the global coordinate frame (O 0 , X 0 , Y 0 , Z 0 ) through a succession of consecutive coordinate transformations which give the following global transformation: Relation (44) expresses the existing mapping between the coordinate reference frame (O i , X i , Y i , Z i ) and the global reference frame (O 0 , X 0 , Y 0 , Z 0 ), regarding the deformation components of the n links, and it allows representation of the jth axial-transverse mode shape of the multi-link manipulator for the chosen configuration, after V i( j) (x i ) and U i( j) (x i ) have been derived from the corresponding differential eigenvalue problem.

Comparison between Analytical Models for a Two-Link Flexible Manipulator
The model presented in this paper is here compared with two existing analytical models developed by Milford and Asokanthan [11] and di Castri et al. [12] for a two-link flexible manipulator.In particular, Milford and Asokanthan [11] refer to a flexible structure and take into account bending contributions only associated to the Euler-Bernoulli beam theory, whilst di Castri et al. [12] takes into account axial beam-by-beam deformations along with the mentioned bending contributions.Moreover, Milford and Asokanthan [11] evaluated only the first three natural frequencies by using an extremely flexible manipulator.The comparison illustrated in this section points out that the proposed model, according to the Timoshenko beam theory, allows to include the other two models as particular cases; moreover, this section also shows that the model based on the Euler-Bernoulli assumptions along with axial contribution can provide a significant superior accuracy to that concerning previous models.
Only when an extremely reduced thickness of the links is used, the three compared models give the same results; differently, once the thickness increases, errors become remarkable and nonnegligible, also exceeding 100%.In this regard, the new matrix formulation introduced by (39) represents an operative mathematical means for executing reliable vibration analysis of this kind of structure.We remark that all data presented in this section refer to natural frequencies without being concerned with a mode shape analysis.
To realize the comparison, the three analytical models were numerically implemented in MATLAB R2007a.The    model presented by Milford and Asokanthan [11] was implemented by using the equations reported into their work; the model presented in this paper and the one by di Castri et al. [12] were implemented through (39), which for this case of two-link manipulator (n = 2) takes the following form: In particular, for the model herein derived in Section 2, according to the Timoshenko beam theory, the blocks in (45) are those reported in (A.1)-(A.4),whereas for the model by di Castri et al. [12]  In order to carry out the related numerical comparisons, the geometric and mass properties of the two-link flexible manipulator are settled as indicated in Table 1.These parameters refer to the extremely flexible manipulator presented by Milford and Asokanthan [11] and are here used for the verification of the predicted data.
For a more significant comparison between models, various simulations have been carried out by varying the slenderness of the structure; more precisely, the two ratios L 1 /h and L T /h are used as reference parameters to characterize each simulation, L 1 being the original value indicated in Table 1, L T = L 1 + L 2 the total manipulator length, and, finally, h the referential thickness value.In order to derive modal data as a function of the slenderness, the two thicknesses h 1 and h 2 (Table 1) are multiplied in each case for a common thickening factor between 1 and 100; in any case, h, indicated into the following numerical tables/simulations, refers to h 2 , which is always the minimum between the two actual thickness values.
Tables 2-6 report the first ten natural frequencies calculated through the three models when the slenderness ratios L 1 /h and L T /h of the manipulator are, respectively, calculated using a common thickening factor equal to 1 (original structure), 10, 25, 50 and 100.In the tables, T Table 4: First ten natural frequencies (Hz); two-link manipulator; (L 1 /h = 33.33,L T /h = 69.33).denotes the model based on the Timoshenko theory, whilst EB † and EB ‡ , respectively, refer to the models by di Castri et al. [12] and Milford and Asokanthan [11].The Timoshenko shear coefficient k i (i = 1, 2), assuming a Poisson's ratio of 0.3, is taken from Cowper [16].
For each table, three nominal configurations of the twolink manipulator are considered, that is, when θ 2 is equal to 0 • , 45 • , and 90 • ; for each configuration, besides the ten natural frequencies, the occurring natural frequency changes are reported, calculated as where EB refers both to EB † and EB ‡ .
A perusal of Tables 2-6 reveals the following interesting findings.
Above any consideration, the model based on Timoshenko theory herein introduced always provides lower frequencies than the counterparts predicted by Milford and Asokanthan [11] and di Castri et al. [12].This can be verified especially when thicker beams are taken into account (100, 50 and 25 times the original h 1 and h 2 of Table 1).It can be of a certain value to notice that the shear deformation theory provides an angle-correction dependence (i.e., the discrepancy Δ depends on the particular configuration established by θ 2 ); this dependence is clear for the case of Table 6 but decreases and becomes inconsistent for thinner beams (Table 2).
Interestingly, a significant discrepancy must be noted not only between the T-model and EB-models by Milford and Asokanthan [11] and di Castri et al. [12] but also between both these latter EB-models, as made explicit in this work through (B.1)-(B.5).In particular, the numerical comparison between the Euler-Bernoulli models shows discrepancies which are extremely large for thicker beams and for an angle θ 2 different by zero; in this regard, significant discrepancies (major than 2% up to the ∼20%) can even be noticed when the simpler model EB † (if compared to the Timoshenko model) is able to provide frequency values with acceptable accuracy (lower than ∼2%).Such percent discrepancies (between the EB models) are also evaluated for values close to decades and can become negligible for extremely flexible manipulators only (L/h > 80 in Tables 2  and 3).
It is stressed, therefore, how the numerical comparison shows the ability of both the T-model and EB-model by di Castri et al. [12] to contain the model by Milford and Asokanthan [11] as a particular case, whilst this latter is advisable to be used only for extremely flexible manipulators; in particular, in the cases herein investigated errors lower than 15% have been achieved for a two-link manipulator having a global length to thickness ratio higher than 80.

Effect of Joint Flexibilities on Modal Data.
All the natural frequencies used for numerical comparisons of Tables 2-6 have been computed through analytical formulations which assume rigid joints only.In particular, this can be easily verified for the matrix formulation (39) herein proposed by simply noting that geometric boundary conditions ( 16) and (19) express perfect cantilevered relations for all the joints.
In order to investigate the effect on modal predictions when the above assumptions are not ensured, simulations have been carried out by employing apposite finite element models of the two-link manipulator of Tables 2-6.Both mechanical and electromechanical stiffnesses of each joint have been modeled as torsional springs having stiffness k T1 for the shoulder joint and k T2 for the elbow one.As in the previous simulations, the first ten natural frequencies have been computed for different slenderness ratios of the robotic structure and for different postures; the simulations have been carried out through finite element models by using Timoshenko beam elements B22 available in the software package ABAQUS R. 6.6 [17]; manipulator meshes are made in order to achieve convergence of finite element predictions to analytical ones.
A torsional stiffness-dependent analysis for each posture and slenderness ratio of Tables 2-6 has been performed; graphical representations of the obtained results are reported by Table 7. Rows and columns represent mode numbers and postures, respectively.Inside each rectangle (for chosen mode i and posture) the prediction on the ith natural frequency is represented versus torsional stiffness k T1 and k T2 changes; the slenderness ratios of Table 2 ( ), Table 3 (×), Table 4 ( ), Table 5 ( * ) and Table 6 (♦) are taken into account.An initial value of stiffness (10 11 Nm/rad) has been calibrated in such a way that finite element predictions match those of rigid joints.
From trends of Table 7, one can note that moving from rigid condition (i.e., k T = 10 11 Nm/rad) toward the pinned boundary condition natural frequencies remain almost constant for each mode up to k T ∼= 10 6 Nm/rad.Only one case differently behaves, that is, when thicker links of Table 6 constitute the manipulator; however, this last case represents an extremely reduced slenderness ratio (L 1 /h = 8.33).Significant frequency changes occur for k T lower than 10 6 Nm/rad although this is not the case for manipulator extremely slender (i.e., L/h = 833).
In conclusion, the influence of flexibility joint on modal predictions must be taken into account along with the slenderness ratios; Table 7 is an attempt aimed at showing the measure of such adependence.Therefore, based on the fact that the proposed matrix formulation regards rigid boundary conditions ( 16)-(21), Table 7 warns the need to introduce flexibilities at joints when extremely stiff manipulators are taken into account.To this end, although the present model can be extensively used for several cases (e.g., Tables 2-5), the authors do not exclude future generalization accounting for flexible joints.

Modal Analysis of an Industrial Manipulator
In this section another comparison between the presented model (T-model) and the model by di Castri et al. [12] is carried out, referring to a typical industrial manipulator with three links; this kind of structure can be met in industrial plants for the handling of loads and their vibration analysis can be of interest for many and various applications.Since the model by di Castri et al. [12] was developed for a twolink flexible manipulator, it has been easily extended to the three-link case by using the new matrix formulation (39) and the equations reported in Appendix B, which have been derived in this work for a generic multilink Euler-Bernoulli manipulator.
Exact modal data were derived through both analytical models and a comparison between the first ten natural frequencies and mode shapes has been reported.A parallel numerical analysis with a finite element software package was carried out to validate the predicted analytical values.
We assume that the industrial manipulator is made of three identical links in aluminium with rectangular hollow         section (Figure 3) whilst its mass and geometric properties are reported in Table 8; M i and J i (i = 1, 2, 3), respectively, represent mass and inertial loads due to the actuators and the payload acting on the system.The comparison between the T-model and the EBmodel is based on four nominal configurations of the industrial manipulator, chosen among its admissible infinite configurations as illustrated in Figure 4; such configurations are obtained by varying the joint variables θ 2 and θ 3 for the desired values.In particular, we settled θ 2,3 = 45 • , 45 • for Figure 4 For each configuration, natural frequencies and mode shapes are derived for both exact models T and EB by using the following system characteristic matrix: whose blocks are, respectively, reported in Appendix A and B.
It is observed that, while natural frequencies are directly calculated through the frequency equation ( 40), mode shapes are plotted by previously obtaining the system eigenfunctions and then following the procedure indicated by ( 41)-(44).
Figures 5(a) and 5(b) report the first ten exact natural frequencies and the translational part of the mode shapes, the occurring natural frequency changes and the finite element results.More precisely, the left column of both the figures reports EB mode shapes and corresponding natural frequencies, calculated through the model by di Castri et al. [12] and finite elements as in parentheses calculated through Euler-Bernoulli beam elements; differently, the right column of the same figures refers to the Timoshenko beam theory results.In the middle of the figures there are the occurring natural frequency changes, calculated through (46).All modes are ordered by frequency.Finally, Tables 9,  10, and 11 report all analytical and numerical results for the other three chosen configurations of the manipulator, without plotting the mode shapes.
Finite element analysis of the structure is carried out by using ABAQUS R. 6.6, using B23 beam elements for the EBmodel mesh, and B22 beam elements for the T-model mesh.
Both meshes consist of 75 elements for each link of the industrial manipulator, yielding a total of 225 elements in the finite element models, value established after a convergence analysis.
As we can see from all tables, finite element obtained natural frequencies agree in all cases of interest with the exact ones; this confirms the validity of the analytical approach presented in this work and of the matrix formulation given by (39).It can be considered a new mathematical tool for reliable vibration analysis of multi-link industrial manipulator without recurring to other models or software packages for numerical investigations.Finally, only referring to the exact results, it is observed that natural frequency changes between the T-model and EB-model are not negligible, with a peak at about 17% for some configurations; therefore, in some situations, as, for example, when forced response is required or when a modal control has to be actuated, the presented model can amend to the lack of accuracy of traditional Euler-Bernoulli models, giving accurate data also for industrial structures.

Conclusions
The configuration dependency of natural frequencies and mode shapes of a planar multi-link flexible manipulator has been investigated.The study has been based on the Timoshenko beam theory and has regarded the formulation and the resolution of the exact partial differential equations describing the system.Detailed boundary conditions along with the differential eigenvalue problem corresponding to each configuration of the manipulator have been derived.A new block matrix formulation herein introduced allows to derive frequency equation associated to a chosen configuration and to calculate corresponding exact modal data.This new formulation allows an absolutely fast and effective analytical methodology for performing the analysis of interest with an arbitrary number of links; a link-toblock relation is established and it is possible to modify the structure by simply adding or deleting a block in the characteristic matrix of the system.The proposed analytical method has been compared against one of the most relevant models proposed in the literature, developed by taking into account only transverse deformations of the links and neglecting the axial ones, which studies the configuration dependency of modal data for a two-link Euler-Bernoulli flexible manipulator.The present work shows that, as the previous kind of approximation (only transverse deformations) can lead to several tenths of a percentage error in predictions with respect to the exact ones, an axial-transverse dynamic characterization for links must be adopted for better and reliable results, as herein proposed by the present authors.The comparison has been carried out after an Euler-Bernoulli multilink formulation has been derived as new; referring to the twolink structure, several slenderness-dependent simulations, carried out for more postures of the manipulator, have stressed (i) the improvement in predictions when the Timoshenko beam theory is adopted and (ii) the abovementioned discrepancies for Euler-Bernoulli formulations whether axial contributions are neglected or not.Moreover, the basic assumption of rigid joints has been validated through extensive simulations aimed at evaluating frequency changes versus flexible joints.
Finally, a case study for a typical three-link industrial manipulator has been investigated; according to the matrix characterizations developed in this work, modal data have been computed through Timoshenko and Euler-Bernoulli formulations for different postures of the robotic structure.It has been pointed out that non negligible differences can result for low modes too; thus, the present work has to be viewed as an efficient tool for designers to rapidly have useful knowledge on free vibration properties of manipulator for structural optimization or control purposes, without recurring to special-purpose software.Finite element analyses have also validated all the analytical models herein proposed.

Figure 2 :
Figure 2: Notation describing of the multi-link flexible manipulators.
(B.1)-(B.5)are used, referring to the Euler-Bernoulli formulation.It is remarked that (B.1)-(B.5)have herein been derived as new and, therefore, are of a certain value.

Table 1 :
Parameters of the flexible manipulator test.
[11]ugh the model by di Castri et al.[12]‡ Through implementation of the model by Milford and Asokanthan[11].

Table 8 :
Parameters of the adopted industrial manipulator.