# Vibration Analysis of Multilink Manipulators Based on Timoshenko Beam Theory

Authors: Carmelo di Castri; Arcangelo Messina
Journal: Journal of Robotics
 (2011)

---

## Abstract

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.


---

## Body

## 1. 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 strain-displacement relations, noninertial frame effects, time-varying 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 non-linear 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 establisheda 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 Section2 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.


## 2. Theoretical Analysis

### 2.1. Problem Formulation

Based on the enumeration ofn 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 Li, cross-sectional area Si, volumetric density ρi, Young’s modulus of elasticity Ei, cross-sectional moment of inertia Ii, shear modulus Gi and shear correction factor ki, for i=1,…,n, where n is the number of links. Each link presents a rigid body of mass Mi and moment of inertia Ji 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.Figure 1
The analyzed multibody system.Figure 2
Notation describing of the multi-link flexible manipulators.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 pointO1 (we say θ1) always has null value. For a chosen configuration, the manipulator is fixed at point O1 and internal blocked joints at points Oi avoid relative rotations between links.The coordinate system (O0,X0,Y0,Z0) is the global one whilst the coordinate system (Oi,Xi,Yi,Zi) is associated to link i. The origin of frame (Oi,Xi,Yi,Zi) is at point Oi and the Xi-axis coincides with the undeformed axis of the ith link. The angle between the X1 and X0 axes, that is, θ1, is always null for the above-mentioned reasons, so that frame (O0,X0,Y0,Z0) and frame (O1,X1,Y1,Z1) will be coincident; the angle between the Xi and Xi-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 pointPi of the ith link, with respect to the correspondent reference system (Oi,Xi,Yi,Zi), can be represented by [xi+ui(xi,t)vi(xi,t)0]T, where the two functions ui(xi,t) and vi(xi,t), respectively, represent the axial and transverse displacements undergone by the considered point.


### 2.2. Mathematical Modelling

In the analysis of theith link, let ψi(xi,t) be the slope of the deflection curve when the shearing force is neglected and βi(xi,t) the angle of shear in the same cross-section; then the total slope of the deformed axis is given by [13](1)∂vi(xi,t)∂xi=ψi(xi,t)+βi(xi,t),
where vi(xi,t) is the transverse deflection of the ith link and ψi(xi,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:(2)Mi(xi,t)=EiIi∂ψi(xi,t)∂xi,Ti(xi,t)=kiGiSi[∂vi(xi,t)∂xi-ψi(xi,t)],Ni(xi,t)=EiSi∂ui(xi,t)∂xi,
where ui(xi,t) is the axial displacement function, Mi(xi,t) is the bending moment of the ith Timoshenko link, and Ti(xi,t) and Ni(xi,t) are the transverse shear force and the axial force of the link, respectively. Considering the link subjected to the transverse distributed load qi(xi,t) and to the axial distributed load  pi(xi,t),  the dynamic equilibrium equations of stress are given by(3)∂Mi(xi,t)∂xi+Ti(xi,t)=ρiIi∂2ψi(xi,t)∂t2,∂Ti(xi,t)∂xi-ρiSi∂2vi(xi,t)∂t2=-qi(xi,t),∂Ni(xi,t)∂xi-ρiSi∂2ui(xi,t)∂t2=-pi(xi,t).When substituting expressions (2) in these latter equations, we get a system of partial differential equations the unknown axial displacement ui(xi,t), for the unknown deflection vi(xi,t), and the unknown rotation ψi(xi,t) describing the dynamics of the ith link: (4)EiIi∂2ψi(xi,t)∂xi2+kiGiSi(∂vi(xi,t)∂xi-ψi(xi,t))=ρiIi∂2ψi(xi,t)∂t2,(5)kiGiSi(∂ψi(xi,t)∂xi-∂2vi(xi,t)∂xi2)+ρiSi∂2vi(xi,t)∂t2=qi(xi,t),(6)EiSi∂2ui(xi,t)∂xi2-ρiSi∂2ui(xi,t)∂t2=-pi(xi,t).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 vi(xi,t) and the other mode is the transverse shearing deformation βi(xi,t), as indirectly measured by (1) and the bending rotation ψi(xi,t). Finally, (6) is the mode of axial deformation of the link, as measured by ui(xi,t).For each link of the manipulator there will be a system of partial differential equations describing the deformations of the link; if we considern 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.


### 2.3. Modal Analysis

The free vibration problem of the system is now considered; for such a case, the distributed forcesqi(xi,t)  and pi(xi,t) are zero, and, in this regard, the dynamics of the ith link is described by the following partial differential equations:(7)EiIi∂2ψi(xi,t)∂xi2+kiGiSi(∂vi(xi,t)∂xi-ψi(xi,t))=ρiIi∂2ψi(xi,t)∂t2,(8)kiGiSi(∂ψi(xi,t)∂xi-∂2vi(xi,t)∂xi2)+ρiSi∂2vi(xi,t)∂t2=0,(9)EiSi∂2ui(xi,t)∂xi2-ρiSi∂2ui(xi,t)∂t2=0.Coupled (7) and (8) can be reduced to a single equation both for vi(xi,t) and ψi(xi,t), having [14](10)EiIiρiSi∂4vi(xi,t)∂xi4-IiSi(1+EikiGi)∂4vi(xi,t)∂xi2∂t2+∂2vi(xi,t)∂t2+ρiIikiGiSi∂4vi(xi,t)∂t4=0
for the transverse deflection vi(xi,t) and (11)EiIiρiSi∂4ψi(xi,t)∂xi4-IiSi(1+EikiGi)∂4ψi(xi,t)∂xi2∂t2+∂2ψi(xi,t)∂t2+ρiIikiGiSi∂4ψi(xi,t)∂t4=0
for the rotation ψi(xi,t).Now, assuming separation of variables in the form(12)ψi(xi,t)=Ψi(xi)ejωt,vi(xi,t)=Vi(xi)ejωt,ui(xi,t)=Ui(xi)ejωt
with natural frequency ω, (11), (10), and (9) can be, respectively, expressed as ordinary differential equations in the following form: (13)EiIiρiSiΨi′′′′(xi)+ω2IiSi(1+EikiGi)Ψi′′(xi)-ω2Ψi(xi)+ω4ρiIikiGiSiΨi(xi)=0,(14)EiIiρiSiVi′′′′(xi)+ω2IiSi(1+EikiGi)Vi′′(xi)-ω2Vi(xi)+ω4ρiIikiGiSiVi(xi)=0,(15)EiSiUi′′(xi)+ω2ρiSiUi(xi)=0,
with unknown parameter ω and unknown functions Ψi(xi), Vi(xi) and Ui(xi).Now, given the above3×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 jointsOi depend on the joint variables θi, it results that boundary conditions are configuration dependent.Figure2 shows all the specified sections of the manipulator for the formulation of the boundary conditions; they correspond to the fixed joint at O1, where three equations will be written, the inner joints Oi(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)
O1 section:
(16)ψ1(x1,t)|x1=0=0,(17)v1(x1,t)|x1=0=0,(18)u1(x1,t)|x1=0=0,(ii)
Oi section (i=2,…,n):
(19)ψi-1(xi-1,t)|xi-1=Li-1-ψi(xi,t)|xi=0=0,(20)vi-1(xi-1,t)|xi-1=Li-1sinθi+ui-1(xi-1,t)|xi-1=Li-1cos⁡θi-ui(xi,t)|xi=0=0,vi-1(xi-1,t)|xi-1=Li-1cos⁡θi-ui-1(xi-1,t)|xi-1=Li-1sinθi-vi(xi,t)|xi=0=0,Ei-1Ii-1∂ψi-1(xi-1,t)∂xi-1|xi-1=Li-1+Ji-1∂2ψi-1(xi-1,t)∂t2|xi-1=Li-1-EiIi∂ψi(xi,t)∂xi|xi=0=0,ki-1Gi-1Si-1[∂vi-1(xi-1,t)∂xi-1-ψi-1(xi-1,t)]|xi-1=Li-1+Mi-1∂2vi-1(xi-1,t)∂t2|xi-1=Li-1-EiSi∂ui(xi,t)∂xi|xi=0sinθi-kiGiSi[∂vi(xi,t)∂xi-ψi(xi,t)]|xi=0cos⁡θi=0,Ei-1Si-1∂ui-1(xi-1,t)∂xi-1|xi-1=Li-1+Mi-1∂2ui-1(xi-1,t)∂t2|xi-1=Li-1-EiSi∂ui(xi,t)∂xi|xi=0cos⁡θi+kiGiSi[∂vi(xi,t)∂xi-ψi(xi,t)]|xi=0sinθi=0,(iii)
EE section (end-effector)(21)EnIn∂ψn(xn,t)∂xn|xn=Ln+Jn∂2ψn(xn,t)∂t2|xn=Ln=0,knGnSn[∂vn(xn,t)∂xn-ψn(xn,t)]|xn=Ln+Mn∂2vn(xn,t)∂t2|xn=Ln=0,EnSn∂un(xn,t)∂xn|xn=Ln+Mn∂2un(xn,t)∂t2|xn=Ln=0.Using solutions (12) into (16)–(21), the boundary conditions for the analyzed reference configuration of the manipulator reduce to(i)
O1 section:
(22)Ψ1(x1)|x1=0=0,(23)V1(x1)|x1=0=0,(24)U1(x1)|x1=0=0,(ii)
Oi section (i=2,…,n):(25)Ψi-1(xi-1)|xi-1=Li-1-Ψi(xi)|xi=0=0Vi-1(xi-1)|xi-1=Li-1sinθi+Ui-1(xi-1)|xi-1=Li-1cos⁡θi-Ui(xi)|xi=0=0,Vi-1(xi-1)|xi-1=Li-1cos⁡θi-Ui-1(xi-1)|xi-1=Li-1sinθi-Vi(xi)|xi=0=0,Ei-1Ii-1Ψi-1′(xi-1)|xi-1=Li-1-ω2Ji-1Ψi-1(xi-1)|xi-1=Li-1-EiIiΨi(xi)|xi=0=0,ki-1Gi-1Si-1[Vi-1′(xi-1)-Ψi-1(xi-1)]|xi-1=Li-1-ω2Mi-1Vi-1(xi-1)|xi-1=Li-1-EiSiUi′(xi)|xi=0sinθi-kiGiSi[Vi′(xi)-Ψi(xi)]|xi=0cos⁡θi=0,Ei-1Si-1Ui-1′(xi-1)|xi-1=Li-1-ω2Mi-1Ui-1(xi-1)|xi-1=Li-1-EiSiUi′(xi)|xi=0cos⁡θi,+kiGiSi[Vi′(xi)-Ψi(xi)]|xi=0sinθi=0,(iii)
EE section (end-effector):(26)EnInΨn′(xn)|xn=Ln-ω2JnΨn(xn)|xn=Ln=0,knGnSn[Vn′(xn)-Ψn(xn)]|xn=Ln-ω2MnVn(xn)|xn=Ln=0,EnSnUn′(xn)|xn=Ln-ω2MnUn(xn)|xn=Ln=0.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 Ψi(xi), Vi(xi), and Ui(xi) 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:
(27)EiIiVi′′′′(xi)+ω2ρiIi(1+EikiGi)Vi′′(xi)+ω2(ω2ρi2IikiGi-ρiSi)Vi(xi)=0.From the characteristic equation associated to (27) it is possible to calculate the following two solutions for each ωfrequency:(28)λi12=-ρiIiω2(1+(Ei/kiGi))+Ω2EiIi,(29)λi22=-ρiIiω2(1+(Ei/kiGi))-Ω2EiIi,
where Ω is a positive or zero quantity defined by(30)Ω=ω4ρi2Ii2(1+EikiGi)2-4EiIiω2(ω2ρi2IikiGi-ρiSi)=ω4ρi2Ii2(1-EikiGi)2+4ω2ρiSiEiIi.It is observed that forω=[(kiGiSi)/(ρiIi)]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(31)Ψi(xi)=A̅i1cosh⁡(λi1xi)+A̅i2sinh⁡(λi1xi)+A̅i3cos⁡(λi2xi)+A̅i4sin(λi2xi),(32)Vi(xi)=Ai1cosh⁡(λi1xi)+Ai2sinh⁡(λi1xi)+Ai3cos⁡(λi2xi)+Ai4sin(λi2xi),(33)Ui(xi)=Bi1cos⁡(ωxici)+Bi2sin(ωxici),
where ci=(Ei/ρi)1/2 is the well-known longitudinal wave speed.The coefficientsA̅ik and Aik (k=1,…,4) are not independent, but related through (8); in fact, substituting solutions (12) into it for vi(xi,t) and ψi(xi,t) and using expressions (31) and (32), we obtain the following relations:(34)A̅i1=ti1Ai2,A̅i2=ti1Ai1,A̅i3=ti2Ai4,A̅i4=-ti2Ai3.where(35)ti1=1λi1(λi12+ω2ρikiGi),ti2=1λi2(λi22-ω2ρikiGi).Moreover, inserting (34) into (31), the final form for the solution Ψi(xi) is given by(36)Ψi(xi)=ti1Ai1sinh⁡(λi1xi)+ti1Ai2cosh⁡(λi1xi)-ti2Ai3sin(λi2xi)+ti2Ai4cos⁡(λi2xi).Now, the expressions forΨi(xi), Vi(xi), and Ui(xi) 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 Aik and Biz (i=1,…,n, z=1,2, and k=1,…,4). These equations can be expressed in matrix form as (37)M(ω)c=0,
where(38)c=[A11A12A13A14B11B12⋯Ai1Ai2Ai3Ai4Bi1Bi2⋯An1An2An3An4Bn1Bn2]T
is the 6×n vector of modal coefficients characterizing the eigenfunctions of the assigned differential eigenvalue problem, whereas(39)M(ω)=(⌊MO1(ω)⌋3×6[0]3×6⋯⋯⋯[0]3×6[MO2,1(ω)]6×6[MO2,2(ω)]6×6[0]6×6[0]6×6[0]6×6⋱⋮⋮[MOi.i-1(ω)]6×6[MOi,i(ω)]6×6⋱⋮⋮⋱[0]6×6[0]6×6[0]6×6[MOn,n-1(ω)]6×6[MOn,n(ω)]6×6[0]3×6⋯⋯⋯[0]3×6[MOEE(ω)]3×6)is a square matrix of order6×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, MO1(ω) and MOEE(ω) are blocks corresponding with the two extreme sections of the structure O1 and EE, respectively; these blocks have dimensions 3×6 since both represent three boundary conditions and, respectively, concern the modal coefficients (A11,A12,A13,A14,B11,B12) and (An1,An2,An3,An4,Bn1,Bn2).MOi.i-1(ω) and MOi.i(ω) are 6×6 blocks and can be interpreted as a matrix representation of the six boundary conditions at the inner blocked joint Oi between link i-1 and link i, concerning the twelve modal coefficients (As1,As2,As3,As4,Bs1,Bs2), 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 Aik and Biz, (i=1,…,n, z=1,2, and k=1,…,4) once the determinant of the characteristic matrix M(ω) is settled to be zero:(40)det⁡[M(ω)]=0.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 Vi(j)(xi) and Ui(j)(xi) that represent the jth axial-transverse mode shape, and the eigenfunction Ψi(j)(xi) that represents the jth cross-section rotation mode shape. Based on relations (32), (33), and (36), the required modal coefficients Aik(j) and Biz(j) for the above-mentioned 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 deformationVi(j)(xi) and Ui(j)(xi), let us suppose to represent the deformed shape of the multi-link manipulator in the global coordinate frame (O0,X0,Y0,Z0); for this scope, it is primarily observed that the homogeneous transformation matrix (43) settles the coordinate transformation between frames (Oi,Xi,Yi,Zi) and (Oi-1,Xi-1,Yi-1,Zi-1), for i=1,…,n, that is [15],(41)Aii-1=(Rii-1dii-10T1),
where Rii-1 is the rotation matrix describing the orientation of frame i relative to frame i-1 and has the following form:(42)Rii-1=(cos⁡θi-sinθi0sinθicos⁡θi0001),
whilst dii-1=[Li-100]T denotes the position of the origin Oi of coordinate frame i relative to coordinate frame i-1; by definition, we set L0 equal to zero. Now, the position of any point of the ith deformed link is expressed in the coordinate frame (Oi,Xi,Yi,Zi) through the homogeneous position vector (43)pi=[xi+Ui(j)(xi)Vi(j)(xi)01].Finally, such a point is directly expressed in the global coordinate frame (O0,X0,Y0,Z0) through a succession of consecutive coordinate transformations which give the following global transformation:(44)p0=A10A21⋯Aii-1pi.Relation (44) expresses the existing mapping between the coordinate reference frame (Oi,Xi,Yi,Zi) and the global reference frame (O0,X0,Y0,Z0), 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 Vi(j)(xi) and Ui(j)(xi) have been derived from the corresponding differential eigenvalue problem.


## 2.1. Problem Formulation

Based on the enumeration ofn 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 Li, cross-sectional area Si, volumetric density ρi, Young’s modulus of elasticity Ei, cross-sectional moment of inertia Ii, shear modulus Gi and shear correction factor ki, for i=1,…,n, where n is the number of links. Each link presents a rigid body of mass Mi and moment of inertia Ji 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.Figure 1
The analyzed multibody system.Figure 2
Notation describing of the multi-link flexible manipulators.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 pointO1 (we say θ1) always has null value. For a chosen configuration, the manipulator is fixed at point O1 and internal blocked joints at points Oi avoid relative rotations between links.The coordinate system (O0,X0,Y0,Z0) is the global one whilst the coordinate system (Oi,Xi,Yi,Zi) is associated to link i. The origin of frame (Oi,Xi,Yi,Zi) is at point Oi and the Xi-axis coincides with the undeformed axis of the ith link. The angle between the X1 and X0 axes, that is, θ1, is always null for the above-mentioned reasons, so that frame (O0,X0,Y0,Z0) and frame (O1,X1,Y1,Z1) will be coincident; the angle between the Xi and Xi-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 pointPi of the ith link, with respect to the correspondent reference system (Oi,Xi,Yi,Zi), can be represented by [xi+ui(xi,t)vi(xi,t)0]T, where the two functions ui(xi,t) and vi(xi,t), respectively, represent the axial and transverse displacements undergone by the considered point.


## 2.2. Mathematical Modelling

In the analysis of theith link, let ψi(xi,t) be the slope of the deflection curve when the shearing force is neglected and βi(xi,t) the angle of shear in the same cross-section; then the total slope of the deformed axis is given by [13](1)∂vi(xi,t)∂xi=ψi(xi,t)+βi(xi,t),
where vi(xi,t) is the transverse deflection of the ith link and ψi(xi,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:(2)Mi(xi,t)=EiIi∂ψi(xi,t)∂xi,Ti(xi,t)=kiGiSi[∂vi(xi,t)∂xi-ψi(xi,t)],Ni(xi,t)=EiSi∂ui(xi,t)∂xi,
where ui(xi,t) is the axial displacement function, Mi(xi,t) is the bending moment of the ith Timoshenko link, and Ti(xi,t) and Ni(xi,t) are the transverse shear force and the axial force of the link, respectively. Considering the link subjected to the transverse distributed load qi(xi,t) and to the axial distributed load  pi(xi,t),  the dynamic equilibrium equations of stress are given by(3)∂Mi(xi,t)∂xi+Ti(xi,t)=ρiIi∂2ψi(xi,t)∂t2,∂Ti(xi,t)∂xi-ρiSi∂2vi(xi,t)∂t2=-qi(xi,t),∂Ni(xi,t)∂xi-ρiSi∂2ui(xi,t)∂t2=-pi(xi,t).When substituting expressions (2) in these latter equations, we get a system of partial differential equations the unknown axial displacement ui(xi,t), for the unknown deflection vi(xi,t), and the unknown rotation ψi(xi,t) describing the dynamics of the ith link: (4)EiIi∂2ψi(xi,t)∂xi2+kiGiSi(∂vi(xi,t)∂xi-ψi(xi,t))=ρiIi∂2ψi(xi,t)∂t2,(5)kiGiSi(∂ψi(xi,t)∂xi-∂2vi(xi,t)∂xi2)+ρiSi∂2vi(xi,t)∂t2=qi(xi,t),(6)EiSi∂2ui(xi,t)∂xi2-ρiSi∂2ui(xi,t)∂t2=-pi(xi,t).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 vi(xi,t) and the other mode is the transverse shearing deformation βi(xi,t), as indirectly measured by (1) and the bending rotation ψi(xi,t). Finally, (6) is the mode of axial deformation of the link, as measured by ui(xi,t).For each link of the manipulator there will be a system of partial differential equations describing the deformations of the link; if we considern 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.


## 2.3. Modal Analysis

The free vibration problem of the system is now considered; for such a case, the distributed forcesqi(xi,t)  and pi(xi,t) are zero, and, in this regard, the dynamics of the ith link is described by the following partial differential equations:(7)EiIi∂2ψi(xi,t)∂xi2+kiGiSi(∂vi(xi,t)∂xi-ψi(xi,t))=ρiIi∂2ψi(xi,t)∂t2,(8)kiGiSi(∂ψi(xi,t)∂xi-∂2vi(xi,t)∂xi2)+ρiSi∂2vi(xi,t)∂t2=0,(9)EiSi∂2ui(xi,t)∂xi2-ρiSi∂2ui(xi,t)∂t2=0.Coupled (7) and (8) can be reduced to a single equation both for vi(xi,t) and ψi(xi,t), having [14](10)EiIiρiSi∂4vi(xi,t)∂xi4-IiSi(1+EikiGi)∂4vi(xi,t)∂xi2∂t2+∂2vi(xi,t)∂t2+ρiIikiGiSi∂4vi(xi,t)∂t4=0
for the transverse deflection vi(xi,t) and (11)EiIiρiSi∂4ψi(xi,t)∂xi4-IiSi(1+EikiGi)∂4ψi(xi,t)∂xi2∂t2+∂2ψi(xi,t)∂t2+ρiIikiGiSi∂4ψi(xi,t)∂t4=0
for the rotation ψi(xi,t).Now, assuming separation of variables in the form(12)ψi(xi,t)=Ψi(xi)ejωt,vi(xi,t)=Vi(xi)ejωt,ui(xi,t)=Ui(xi)ejωt
with natural frequency ω, (11), (10), and (9) can be, respectively, expressed as ordinary differential equations in the following form: (13)EiIiρiSiΨi′′′′(xi)+ω2IiSi(1+EikiGi)Ψi′′(xi)-ω2Ψi(xi)+ω4ρiIikiGiSiΨi(xi)=0,(14)EiIiρiSiVi′′′′(xi)+ω2IiSi(1+EikiGi)Vi′′(xi)-ω2Vi(xi)+ω4ρiIikiGiSiVi(xi)=0,(15)EiSiUi′′(xi)+ω2ρiSiUi(xi)=0,
with unknown parameter ω and unknown functions Ψi(xi), Vi(xi) and Ui(xi).Now, given the above3×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 jointsOi depend on the joint variables θi, it results that boundary conditions are configuration dependent.Figure2 shows all the specified sections of the manipulator for the formulation of the boundary conditions; they correspond to the fixed joint at O1, where three equations will be written, the inner joints Oi(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)
O1 section:
(16)ψ1(x1,t)|x1=0=0,(17)v1(x1,t)|x1=0=0,(18)u1(x1,t)|x1=0=0,(ii)
Oi section (i=2,…,n):
(19)ψi-1(xi-1,t)|xi-1=Li-1-ψi(xi,t)|xi=0=0,(20)vi-1(xi-1,t)|xi-1=Li-1sinθi+ui-1(xi-1,t)|xi-1=Li-1cos⁡θi-ui(xi,t)|xi=0=0,vi-1(xi-1,t)|xi-1=Li-1cos⁡θi-ui-1(xi-1,t)|xi-1=Li-1sinθi-vi(xi,t)|xi=0=0,Ei-1Ii-1∂ψi-1(xi-1,t)∂xi-1|xi-1=Li-1+Ji-1∂2ψi-1(xi-1,t)∂t2|xi-1=Li-1-EiIi∂ψi(xi,t)∂xi|xi=0=0,ki-1Gi-1Si-1[∂vi-1(xi-1,t)∂xi-1-ψi-1(xi-1,t)]|xi-1=Li-1+Mi-1∂2vi-1(xi-1,t)∂t2|xi-1=Li-1-EiSi∂ui(xi,t)∂xi|xi=0sinθi-kiGiSi[∂vi(xi,t)∂xi-ψi(xi,t)]|xi=0cos⁡θi=0,Ei-1Si-1∂ui-1(xi-1,t)∂xi-1|xi-1=Li-1+Mi-1∂2ui-1(xi-1,t)∂t2|xi-1=Li-1-EiSi∂ui(xi,t)∂xi|xi=0cos⁡θi+kiGiSi[∂vi(xi,t)∂xi-ψi(xi,t)]|xi=0sinθi=0,(iii)
EE section (end-effector)(21)EnIn∂ψn(xn,t)∂xn|xn=Ln+Jn∂2ψn(xn,t)∂t2|xn=Ln=0,knGnSn[∂vn(xn,t)∂xn-ψn(xn,t)]|xn=Ln+Mn∂2vn(xn,t)∂t2|xn=Ln=0,EnSn∂un(xn,t)∂xn|xn=Ln+Mn∂2un(xn,t)∂t2|xn=Ln=0.Using solutions (12) into (16)–(21), the boundary conditions for the analyzed reference configuration of the manipulator reduce to(i)
O1 section:
(22)Ψ1(x1)|x1=0=0,(23)V1(x1)|x1=0=0,(24)U1(x1)|x1=0=0,(ii)
Oi section (i=2,…,n):(25)Ψi-1(xi-1)|xi-1=Li-1-Ψi(xi)|xi=0=0Vi-1(xi-1)|xi-1=Li-1sinθi+Ui-1(xi-1)|xi-1=Li-1cos⁡θi-Ui(xi)|xi=0=0,Vi-1(xi-1)|xi-1=Li-1cos⁡θi-Ui-1(xi-1)|xi-1=Li-1sinθi-Vi(xi)|xi=0=0,Ei-1Ii-1Ψi-1′(xi-1)|xi-1=Li-1-ω2Ji-1Ψi-1(xi-1)|xi-1=Li-1-EiIiΨi(xi)|xi=0=0,ki-1Gi-1Si-1[Vi-1′(xi-1)-Ψi-1(xi-1)]|xi-1=Li-1-ω2Mi-1Vi-1(xi-1)|xi-1=Li-1-EiSiUi′(xi)|xi=0sinθi-kiGiSi[Vi′(xi)-Ψi(xi)]|xi=0cos⁡θi=0,Ei-1Si-1Ui-1′(xi-1)|xi-1=Li-1-ω2Mi-1Ui-1(xi-1)|xi-1=Li-1-EiSiUi′(xi)|xi=0cos⁡θi,+kiGiSi[Vi′(xi)-Ψi(xi)]|xi=0sinθi=0,(iii)
EE section (end-effector):(26)EnInΨn′(xn)|xn=Ln-ω2JnΨn(xn)|xn=Ln=0,knGnSn[Vn′(xn)-Ψn(xn)]|xn=Ln-ω2MnVn(xn)|xn=Ln=0,EnSnUn′(xn)|xn=Ln-ω2MnUn(xn)|xn=Ln=0.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 Ψi(xi), Vi(xi), and Ui(xi) 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:
(27)EiIiVi′′′′(xi)+ω2ρiIi(1+EikiGi)Vi′′(xi)+ω2(ω2ρi2IikiGi-ρiSi)Vi(xi)=0.From the characteristic equation associated to (27) it is possible to calculate the following two solutions for each ωfrequency:(28)λi12=-ρiIiω2(1+(Ei/kiGi))+Ω2EiIi,(29)λi22=-ρiIiω2(1+(Ei/kiGi))-Ω2EiIi,
where Ω is a positive or zero quantity defined by(30)Ω=ω4ρi2Ii2(1+EikiGi)2-4EiIiω2(ω2ρi2IikiGi-ρiSi)=ω4ρi2Ii2(1-EikiGi)2+4ω2ρiSiEiIi.It is observed that forω=[(kiGiSi)/(ρiIi)]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(31)Ψi(xi)=A̅i1cosh⁡(λi1xi)+A̅i2sinh⁡(λi1xi)+A̅i3cos⁡(λi2xi)+A̅i4sin(λi2xi),(32)Vi(xi)=Ai1cosh⁡(λi1xi)+Ai2sinh⁡(λi1xi)+Ai3cos⁡(λi2xi)+Ai4sin(λi2xi),(33)Ui(xi)=Bi1cos⁡(ωxici)+Bi2sin(ωxici),
where ci=(Ei/ρi)1/2 is the well-known longitudinal wave speed.The coefficientsA̅ik and Aik (k=1,…,4) are not independent, but related through (8); in fact, substituting solutions (12) into it for vi(xi,t) and ψi(xi,t) and using expressions (31) and (32), we obtain the following relations:(34)A̅i1=ti1Ai2,A̅i2=ti1Ai1,A̅i3=ti2Ai4,A̅i4=-ti2Ai3.where(35)ti1=1λi1(λi12+ω2ρikiGi),ti2=1λi2(λi22-ω2ρikiGi).Moreover, inserting (34) into (31), the final form for the solution Ψi(xi) is given by(36)Ψi(xi)=ti1Ai1sinh⁡(λi1xi)+ti1Ai2cosh⁡(λi1xi)-ti2Ai3sin(λi2xi)+ti2Ai4cos⁡(λi2xi).Now, the expressions forΨi(xi), Vi(xi), and Ui(xi) 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 Aik and Biz (i=1,…,n, z=1,2, and k=1,…,4). These equations can be expressed in matrix form as (37)M(ω)c=0,
where(38)c=[A11A12A13A14B11B12⋯Ai1Ai2Ai3Ai4Bi1Bi2⋯An1An2An3An4Bn1Bn2]T
is the 6×n vector of modal coefficients characterizing the eigenfunctions of the assigned differential eigenvalue problem, whereas(39)M(ω)=(⌊MO1(ω)⌋3×6[0]3×6⋯⋯⋯[0]3×6[MO2,1(ω)]6×6[MO2,2(ω)]6×6[0]6×6[0]6×6[0]6×6⋱⋮⋮[MOi.i-1(ω)]6×6[MOi,i(ω)]6×6⋱⋮⋮⋱[0]6×6[0]6×6[0]6×6[MOn,n-1(ω)]6×6[MOn,n(ω)]6×6[0]3×6⋯⋯⋯[0]3×6[MOEE(ω)]3×6)is a square matrix of order6×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, MO1(ω) and MOEE(ω) are blocks corresponding with the two extreme sections of the structure O1 and EE, respectively; these blocks have dimensions 3×6 since both represent three boundary conditions and, respectively, concern the modal coefficients (A11,A12,A13,A14,B11,B12) and (An1,An2,An3,An4,Bn1,Bn2).MOi.i-1(ω) and MOi.i(ω) are 6×6 blocks and can be interpreted as a matrix representation of the six boundary conditions at the inner blocked joint Oi between link i-1 and link i, concerning the twelve modal coefficients (As1,As2,As3,As4,Bs1,Bs2), 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 Aik and Biz, (i=1,…,n, z=1,2, and k=1,…,4) once the determinant of the characteristic matrix M(ω) is settled to be zero:(40)det⁡[M(ω)]=0.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 Vi(j)(xi) and Ui(j)(xi) that represent the jth axial-transverse mode shape, and the eigenfunction Ψi(j)(xi) that represents the jth cross-section rotation mode shape. Based on relations (32), (33), and (36), the required modal coefficients Aik(j) and Biz(j) for the above-mentioned 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 deformationVi(j)(xi) and Ui(j)(xi), let us suppose to represent the deformed shape of the multi-link manipulator in the global coordinate frame (O0,X0,Y0,Z0); for this scope, it is primarily observed that the homogeneous transformation matrix (43) settles the coordinate transformation between frames (Oi,Xi,Yi,Zi) and (Oi-1,Xi-1,Yi-1,Zi-1), for i=1,…,n, that is [15],(41)Aii-1=(Rii-1dii-10T1),
where Rii-1 is the rotation matrix describing the orientation of frame i relative to frame i-1 and has the following form:(42)Rii-1=(cos⁡θi-sinθi0sinθicos⁡θi0001),
whilst dii-1=[Li-100]T denotes the position of the origin Oi of coordinate frame i relative to coordinate frame i-1; by definition, we set L0 equal to zero. Now, the position of any point of the ith deformed link is expressed in the coordinate frame (Oi,Xi,Yi,Zi) through the homogeneous position vector (43)pi=[xi+Ui(j)(xi)Vi(j)(xi)01].Finally, such a point is directly expressed in the global coordinate frame (O0,X0,Y0,Z0) through a succession of consecutive coordinate transformations which give the following global transformation:(44)p0=A10A21⋯Aii-1pi.Relation (44) expresses the existing mapping between the coordinate reference frame (Oi,Xi,Yi,Zi) and the global reference frame (O0,X0,Y0,Z0), 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 Vi(j)(xi) and Ui(j)(xi) have been derived from the corresponding differential eigenvalue problem.


## 3. 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:(45)M(ω)=([MO1(ω)]3×6[0]3×6[MO2,1(ω)]6×6[MO2,2(ω)]6×6[0]3×6[MOEE(ω)]3×6).In particular, for the model herein derived in Section2, 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] (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.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 Table1. 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.Table 1
Parameters of the flexible manipulator test.
ParameterSymbolUnitsLink 1Link 2LengthLim0.50.54Thicknesshimm0.610.6Widthbicm5.15.1Massmikg0.1210.129StiffnessEiIiN m20.2030.193Tip load massMikg0.610.224Tip load inertiaJikg m21.5·10-31.3·10-4For a more significant comparison between models, various simulations have been carried out by varying the slenderness of the structure; more precisely, the two ratiosL1/h and LT/h are used as reference parameters to characterize each simulation, L1 being the original value indicated in Table 1, LT=L1+L2 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 h1 and h2 (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 h2, which is always the minimum between the two actual thickness values.Tables2–6 report the first ten natural frequencies calculated through the three models when the slenderness ratios L1/h and LT/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 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 ki (i=1,2), assuming a Poisson’s ratio of 0.3, is taken from Cowper [16].Table 2
First ten natural frequencies (Hz); two-link manipulator, see Table1; (L1/h=833.3, LT/h=1733).
θ2 = 0°θ2 = 45°θ2 = 90°ModeEB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡10.20020.20020.20020.000.000.21140.21140.21140.000.000.24900.24900.24900.000.0020.95570.95570.95570.000.000.82690.82690.82690.000.000.65310.65310.65310.000.0034.89014.89014.89010.000.004.85584.85584.85580.000.004.82554.82554.82550.000.0049.07449.07449.07440.000.008.99248.99248.99240.000.008.91618.91618.91610.000.00514.50014.50014.4990.010.0114.46014.46014.4600.000.0014.42814.42814.4280.000.00620.65020.65020.6500.000.0020.57920.57920.5780.000.0020.53520.53520.5340.000.00736.44736.44736.4460.000.0036.29836.29836.2970.000.0036.15036.15036.1490.000.00836.74636.74636.7450.000.0036.75636.75636.7550.000.0036.81236.81236.8110.000.00962.87262.87262.8700.000.0062.75962.76062.7570.000.0062.68262.68262.6800.000.001070.91770.91770.9140.000.0070.87070.87070.8670.000.0070.84270.84270.8390.000.00†Through the model by di Castri et al. [12] ‡Through implementation of the model by Milford and Asokanthan [11].Table 3
First ten natural frequencies (Hz); two-link manipulator; (L1/h=83.33, LT/h=173.3).
θ2=  0°θ2 =  45°θ2 =  90°ModeEB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡13.93443.93443.93430.000.004.19074.19074.19060.000.005.09875.09875.09860.000.00221.55021.55021.5460.020.0217.92517.92617.9230.010.0213.55413.55413.5520.010.01370.72070.72070.6950.040.0468.90168.90868.8780.030.0468.12068.13068.0980.030.054130.00130.00129.910.070.07122.38122.41122.290.070.10117.10117.11117.010.080.095218.26218.26218.070.090.09214.44214.52214.260.080.12212.37212.43212.190.080.116317.00317.00316.590.130.13308.98309.23308.560.140.22304.46304.57304.030.140.187425.34425.34424.490.200.20418.99419.31418.170.200.27415.96416.06415.150.200.228545.20545.20544.050.210.21536.96537.95535.810.210.40534.21534.97533.050.220.369745.10745.10742.100.400.40734.87736.49731.950.400.62731.71732.18728.790.400.4710837.96837.96835.260.320.32826.94830.08824.250.330.71825.38827.93822.680.330.64†Through the model by di Castri et al. [12].‡Through implementation of the model by Milford and Asokanthan [11].Table 4
First ten natural frequencies (Hz); two-link manipulator; (L1/h=33.33, LT/h=69.33).
θ2 = 0°θ2 = 45°θ2 = 90°ModeEB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡110.93110.93110.9290.020.0211.67511.67511.6730.020.0214.36314.36314.3600.020.02263.39963.39963.3340.100.1051.70551.71651.6530.100.1238.34838.35338.3090.100.113190.90190.90190.430.250.25186.38186.49185.950.230.29185.36185.54184.940.230.324356.79356.79355.140.460.46322.97323.92321.400.490.78303.66303.97302.120.510.615602.80602.80598.840.660.66587.48589.46583.770.640.97582.51584.21578.870.630.926882.58882.58873.791.01.0831.00839.66822.291.12.1817.71819.88808.781.11.471162.5Axial1162.50.0—1162.31178.91149.21.12.61162.31168.21148.71.21.781205.51205.51190.71.21.21466.91534.01448.41.35.91505.21519.11482.21.62.591571.61571.61548.71.51.51713.91953.41699.70.8151901.21937.21867.81.83.7101998.11998.11953.32.32.32012.62331.01973.02.0182017.82321.71998.50.9716†Through the model by di Castri et al. [12].‡Through implementation of the model by Milford and Asokanthan [11].Table 5
First ten natural frequencies (Hz); two-link manipulator; (L1/h=16.67, LT/h=34.67).
θ2 = 0°θ2 = 45°θ2 = 90°ModeEB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡122.77122.77122.7560.070.0724.34924.35024.3330.070.0730.10730.11030.0850.070.082135.87135.87135.280.440.44109.58109.68109.120.420.5180.51480.55780.1780.420.473395.12395.12391.111.01.0386.31387.26382.590.971.2385.43387.00381.780.961.44750.85750.85736.831.91.9656.26666.58643.751.93.5615.82618.74603.332.12.651206.2axial1206.20.0—1203.81233.31174.82.55.01203.71224.61173.52.64.461261.01261.01226.12.82.81540.51729.71496.82.9161650.51679.61579.74.56.371855.01855.01779.44.24.21965.72505.11918.42.5312117.02484.52083.41.61982570.12570.12440.45.35.32583.33263.22460.95.0332604.73220.22504.84.02993393.13393.13176.06.86.83267.14163.03080.36.1353177.24128.42973.46.939103567.6axial3567.60.0—3782.35096.13701.12.2383959.85066.63788.94.534†Through the model by di Castri et al. [12].‡Through implementation of the model by Milford and Asokanthan [11].Table 6
First ten natural frequencies (Hz); two-link manipulator; (L1/h=8.33, LT/h=17.33).
θ2 = 0°θ2 = 45°θ2 = 90°ModeEB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡EB†EB‡TΔ†Δ‡146.54046.54046.4170.260.2649.79649.80149.6610.270.2861.73661.76261.5500.300.342282.42282.42277.521.81.8225.58226.53221.881.62.1165.04165.42162.321.61.93805.90805.90773.814.14.1782.63791.77753.863.65.0779.33793.67751.293.55.641229.9axial1229.90.0—1223.51356.81154.65.1181221.21250.71132.47.11051551.21551.21443.17.57.51841.42536.31780.62.4422057.52522.51981.33.02762589.52589.52326.611112621.13515.92362.77.3492737.43400.62515.96.53573653.2none3294.011—3444.15198.83109.76.4673271.85161.42887.07.57983828.93828.93653.24.84.84124.36708.33906.93.2724345.46603.44141.63.15995336.25336.24425.321215253.18728.34403.59.7985165.38663.74372.59.298106146.9none5590.410—6316.7107995461.07.9986548.3107165305.612102†Through the model by di Castri et al. [12].‡Through implementation of the model by Milford and Asokanthan [11].For each table, three nominal configurations of the two-link 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(46)Δj%=frequencyEB-frequencyTfrequencyT⋅100,j=1,…,10,
where EB refers both to EB† and EB‡.A perusal of Tables2–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 h1 and h2 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 theT-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 theT-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.


### 3.1. Effect of Joint Flexibilities on Modal Data

All the natural frequencies used for numerical comparisons of Tables2–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 Tables2–6. Both mechanical and electromechanical stiffnesses of each joint have been modeled as torsional springs having stiffness kT1 for the shoulder joint and kT2 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 Tables2–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 kT1 and kT2 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 (1011 Nm/rad) has been calibrated in such a way that finite element predictions match those of rigid joints.Table 7
Natural frequencies as a function of joint flexibilities†. Slenderness ratios: ○ Table 2;  × Table 3; □ Table 4; * Table 5;  ◊ Table 6.
Mode12345678910†ABAQUS R.6.6, B22 element.From trends of Table7, one can note that moving from rigid condition (i.e., kT = 1011 Nm/rad) toward the pinned boundary condition natural frequencies remain almost constant for each mode up to kT~=106 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 (L1/h=8.33). Significant frequency changes occur for kT lower than 106 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; Table7 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.


## 3.1. Effect of Joint Flexibilities on Modal Data

All the natural frequencies used for numerical comparisons of Tables2–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 Tables2–6. Both mechanical and electromechanical stiffnesses of each joint have been modeled as torsional springs having stiffness kT1 for the shoulder joint and kT2 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 Tables2–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 kT1 and kT2 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 (1011 Nm/rad) has been calibrated in such a way that finite element predictions match those of rigid joints.Table 7
Natural frequencies as a function of joint flexibilities†. Slenderness ratios: ○ Table 2;  × Table 3; □ Table 4; * Table 5;  ◊ Table 6.
Mode12345678910†ABAQUS R.6.6, B22 element.From trends of Table7, one can note that moving from rigid condition (i.e., kT = 1011 Nm/rad) toward the pinned boundary condition natural frequencies remain almost constant for each mode up to kT~=106 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 (L1/h=8.33). Significant frequency changes occur for kT lower than 106 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; Table7 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.


## 4. 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 two-link 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 (Figure3) whilst its mass and geometric properties are reported in Table 8; Mi and  Ji (i=1,2,3), respectively, represent mass and inertial loads due to the actuators and the payload acting on the system.Table 8
Parameters of the adopted industrial manipulator.
ParameterSymbolUnitsLinki (i=1,2)Link 3LengthLim1.51.5Width of sectionhicm1515Height of sectionbicm2020Wall thicknessticm11Young’s modulusEiGPa7070Mass densityρikg/m327002700Poisson’s ratioνinone0.340.34Tip load massMikg1035Tip load inertiaJikg m20.30.5Figure 3
Rectangular hollow section of the adopted industrial manipulator.The comparison between theT-model and the EB-model 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(a), θ2,3=-45°,90° for Figure 4(b), θ2,3=90°, −90° for Figure 4(c), and θ2,3=-90°, −90° for Figure 4(d).Test scenarios: (a)θ2,3=45°, 45°; (b) θ2,3=-45°, 90°; (c) θ2,3=90°, −90°; (d) θ2,3=-90°, −90°.
(a)(b)(c)(d)For each configuration, natural frequencies and mode shapes are derived for both exact modelsT and EB by using the following system characteristic matrix:(47)M(ω)=(⌊MO1(ω)⌋3×6[0]3×6[0]3×6[MO2,1(ω)]6×6[MO2,2(ω)]6×6[0]6×6[0]6×6[MO3,2(ω)]6×6[MO3,3(ω)]6×6[0]3×6[0]3×6[MOEE(ω)]3×6),
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).Figures5(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.Table 9
First ten natural frequencies (Hz); see Table7; θ2,3=-45°,90°.
Analytical solutionFinite Element solutionModeEBTΔ )%)ABAQUS R.6.6, B23ABAQUS R.6.6, B2216.18556.17010.256.18556.1701225.59625.2221.4825.59625.222359.14756.5944.5159.14756.5944208.91199.874.52208.91199.875309.93280.4710.5309.93280.476398.82368.748.16398.82368.747534.23496.227.66534.23496.228647.83571.6513.3647.83571.659762.02702.558.46762.02702.5510890.88802.9011.0890.88802.90Table 10
First ten natural frequencies (Hz); see Table7; θ2,3=90°,-90°.
Analytical solutionFinite Element solutionModeEBTΔ )%)ABAQUS R.6.6, B23ABAQUS R.6.6, B2216.66406.64960.226.66406.6496231.47230.4823.2531.47230.482338.76837.6842.8838.76837.6844218.72208.784.76218.72208.785305.61277.5610.1305.61277.566434.66376.1215.6434.66376.127567.22543.014.46567.22543.018655.55592.0010.7655.55592.009760.44710.077.09760.44710.0710888.40797.1911.4888.40797.19Table 11
First ten natural frequencies (Hz); see Table7; θ2,3=-90°,-90°.
Analytical solutionFinite element solutionModeEBTΔ )%)ABAQUS R.6.6, B23ABAQUS R.6.6, B2218.85528.83220.268.85528.8322217.81717.5331.6217.81717.533350.83648.7774.2250.83648.7774218.16207.775.00218.16207.775318.97290.709.72318.97290.706426.96366.3016.6426.96366.307568.09541.764.86568.09541.768641.25588.199.02641.25588.199760.88708.357.42760.88708.3510875.15790.1710.8875.15790.17(a) Modes (Hz) of an industrial manipulator (see Table8). (FEM: ABAQUS 6.6, B23, B22) left column: EB-model; right column: T-model: θ2=θ3=45°. (b) Modes (Hz) of an industrial manipulator (see Table 8). (FEM: ABAQUS 6.6, B23, B22) left column: EB-model; right column: T-model: θ2=θ3=45°.
(a)(b)Finite element analysis of the structure is carried out by using ABAQUS R. 6.6, using B23 beam elements for theEB-model 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 theT-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.


## 5. 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-to-block 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 two-link 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 above-mentioned 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.



---

*Source: 890258-2011-07-05.xml*