About this Journal Submit a Manuscript Table of Contents
Journal of Robotics
Volume 2011 (2011), Article ID 890258, 21 pages
Research Article

Vibration Analysis of Multilink Manipulators Based on Timoshenko Beam Theory

Dipartimento di Ingegneria dell'Innovazione, Università del Salento, 73100 Lecce, Italy

Received 11 October 2010; Revised 6 March 2011; Accepted 7 March 2011

Academic Editor: Danica Kragic

Copyright © 2011 Carmelo di Castri and Arcangelo Messina. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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.

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 𝑛-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 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.

2. Theoretical Analysis

2.1. Problem Formulation

Based on the enumeration of 𝑛 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 𝐿𝑖, cross-sectional area 𝑆𝑖, volumetric density 𝜌𝑖, Young’s modulus of elasticity 𝐸𝑖, cross-sectional moment of inertia 𝐼𝑖, shear modulus 𝐺𝑖 and shear correction factor 𝑘𝑖, for 𝑖=1,,𝑛, where 𝑛 is the number of links. Each link presents a rigid body of mass 𝑀𝑖 and moment of inertia 𝐽𝑖 at its tip that for 𝑖=1,,𝑛1 accounts for the presence of possible joint motors; 𝑖=𝑛 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 point 𝑂1 (we say 𝜃1) always has null value. For a chosen configuration, the manipulator is fixed at point 𝑂1 and internal blocked joints at points 𝑂𝑖 avoid relative rotations between links.

The coordinate system (𝑂0,𝑋0,𝑌0,𝑍0) is the global one whilst the coordinate system (𝑂𝑖,𝑋𝑖,𝑌𝑖,𝑍𝑖) is associated to link 𝑖. The origin of frame (𝑂𝑖,𝑋𝑖,𝑌𝑖,𝑍𝑖) is at point 𝑂𝑖 and the 𝑋𝑖-axis coincides with the undeformed axis of the 𝑖th link. The angle between the 𝑋1 and 𝑋0 axes, that is, 𝜃1, is always null for the above-mentioned reasons, so that frame (𝑂0,𝑋0,𝑌0,𝑍0) and frame (𝑂1,𝑋1,𝑌1,𝑍1) will be coincident; the angle between the 𝑋𝑖 and 𝑋𝑖1 axes is denoted with 𝜃𝑖, and the set {𝜃𝑖𝑖=2,,𝑛} 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 𝑃𝑖 of the 𝑖th link, with respect to the correspondent reference system (𝑂𝑖,𝑋𝑖,𝑌𝑖,𝑍𝑖), can be represented by [𝑥𝑖+𝑢𝑖(𝑥𝑖,𝑡)𝑣𝑖(𝑥𝑖,𝑡)0]𝑇, where the two functions 𝑢𝑖(𝑥𝑖,𝑡) and 𝑣𝑖(𝑥𝑖,𝑡), respectively, represent the axial and transverse displacements undergone by the considered point.

2.2. Mathematical Modelling

In the analysis of the 𝑖th link, let 𝜓𝑖(𝑥𝑖,𝑡) be the slope of the deflection curve when the shearing force is neglected and 𝛽𝑖(𝑥𝑖,𝑡) the angle of shear in the same cross-section; then the total slope of the deformed axis is given by [13]𝜕𝑣𝑖𝑥𝑖,𝑡𝜕𝑥𝑖=𝜓𝑖𝑥𝑖,𝑡+𝛽𝑖𝑥𝑖,,𝑡(1) where 𝑣𝑖(𝑥𝑖,𝑡) is the transverse deflection of the 𝑖th link and 𝜓𝑖(𝑥𝑖,𝑡) 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) where 𝑢𝑖(𝑥𝑖,𝑡) is the axial displacement function, 𝑀𝑖(𝑥𝑖,𝑡) is the bending moment of the ith Timoshenko link, and 𝑇𝑖(𝑥𝑖,𝑡) and 𝑁𝑖(𝑥𝑖,𝑡) are the transverse shear force and the axial force of the link, respectively. Considering the link subjected to the transverse distributed load 𝑞𝑖(𝑥𝑖,𝑡) and to the axial distributed load  𝑝𝑖(𝑥𝑖,𝑡), the dynamic equilibrium equations of stress are given by𝜕𝑀𝑖𝑥𝑖,𝑡𝜕𝑥𝑖+𝑇𝑖𝑥𝑖,𝑡=𝜌𝑖𝐼𝑖𝜕2𝜓𝑖𝑥𝑖,𝑡𝜕𝑡2,𝜕𝑇𝑖𝑥𝑖,𝑡𝜕𝑥𝑖𝜌𝑖𝑆𝑖𝜕2𝑣𝑖𝑥𝑖,𝑡𝜕𝑡2=𝑞𝑖𝑥𝑖,,𝑡𝜕𝑁𝑖𝑥𝑖,𝑡𝜕𝑥𝑖𝜌𝑖𝑆𝑖𝜕2𝑢𝑖𝑥𝑖,𝑡𝜕𝑡2=𝑝𝑖𝑥𝑖.,𝑡(3)

When substituting expressions (2) in these latter equations, we get a system of partial differential equations the unknown axial displacement 𝑢𝑖(𝑥𝑖,𝑡), for the unknown deflection 𝑣𝑖(𝑥𝑖,𝑡), and the unknown rotation 𝜓𝑖(𝑥𝑖,𝑡) describing the dynamics of the 𝑖th link: 𝐸𝑖𝐼𝑖𝜕2𝜓𝑖𝑥𝑖,𝑡𝜕𝑥2𝑖+𝑘𝑖𝐺𝑖𝑆𝑖𝜕𝑣𝑖𝑥𝑖,𝑡𝜕𝑥𝑖𝜓𝑖𝑥𝑖,𝑡=𝜌𝑖𝐼𝑖𝜕2𝜓𝑖𝑥𝑖,𝑡𝜕𝑡2,(4)𝑘𝑖𝐺𝑖𝑆𝑖𝜕𝜓𝑖𝑥𝑖,𝑡𝜕𝑥𝑖𝜕2𝑣𝑖𝑥𝑖,𝑡𝜕𝑥2𝑖+𝜌𝑖𝑆𝑖𝜕2𝑣𝑖𝑥𝑖,𝑡𝜕𝑡2=𝑞𝑖𝑥𝑖,,𝑡(5)𝐸𝑖𝑆𝑖𝜕2𝑢𝑖𝑥𝑖,𝑡𝜕𝑥2𝑖𝜌𝑖𝑆𝑖𝜕2𝑢𝑖𝑥𝑖,𝑡𝜕𝑡2=𝑝𝑖𝑥𝑖,𝑡.(6)

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 𝑣𝑖(𝑥𝑖,𝑡) and the other mode is the transverse shearing deformation 𝛽𝑖(𝑥𝑖,𝑡), as indirectly measured by (1) and the bending rotation 𝜓𝑖(𝑥𝑖,𝑡). Finally, (6) is the mode of axial deformation of the link, as measured by 𝑢𝑖(𝑥𝑖,𝑡).

For each link of the manipulator there will be a system of partial differential equations describing the deformations of the link; if we consider 𝑛 links, 3×𝑛 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 forces 𝑞𝑖(𝑥𝑖,𝑡) and 𝑝𝑖(𝑥𝑖,𝑡) are zero, and, in this regard, the dynamics of the 𝑖th link is described by the following partial differential equations:𝐸𝑖𝐼𝑖𝜕2𝜓𝑖𝑥𝑖,𝑡𝜕𝑥2𝑖+𝑘𝑖𝐺𝑖𝑆𝑖𝜕𝑣𝑖𝑥𝑖,𝑡𝜕𝑥𝑖𝜓𝑖𝑥𝑖,𝑡=𝜌𝑖𝐼𝑖𝜕2𝜓𝑖𝑥𝑖,𝑡𝜕𝑡2,𝑘(7)𝑖𝐺𝑖𝑆𝑖𝜕𝜓𝑖𝑥𝑖,𝑡𝜕𝑥𝑖𝜕2𝑣𝑖𝑥𝑖,𝑡𝜕𝑥2𝑖+𝜌𝑖𝑆𝑖𝜕2𝑣𝑖𝑥𝑖,𝑡𝜕𝑡2𝐸=0,(8)𝑖𝑆𝑖𝜕2𝑢𝑖𝑥𝑖,𝑡𝜕𝑥2𝑖𝜌𝑖𝑆𝑖𝜕2𝑢𝑖𝑥𝑖,𝑡𝜕𝑡2=0.(9)

Coupled (7) and (8) can be reduced to a single equation both for 𝑣𝑖(𝑥𝑖,𝑡) and 𝜓𝑖(𝑥𝑖,𝑡), having [14]𝐸𝑖𝐼𝑖𝜌𝑖𝑆𝑖𝜕4𝑣𝑖𝑥𝑖,𝑡𝜕𝑥4𝑖𝐼𝑖𝑆𝑖𝐸1+𝑖𝑘𝑖𝐺𝑖𝜕4𝑣𝑖𝑥𝑖,𝑡𝜕𝑥2𝑖𝜕𝑡2+𝜕2𝑣𝑖𝑥𝑖,𝑡𝜕𝑡2+𝜌𝑖𝐼𝑖𝑘𝑖𝐺𝑖𝑆𝑖𝜕4𝑣𝑖𝑥𝑖,𝑡𝜕𝑡4=0(10) for the transverse deflection 𝑣𝑖(𝑥𝑖,𝑡) and 𝐸𝑖𝐼𝑖𝜌𝑖𝑆𝑖𝜕4𝜓𝑖𝑥𝑖,𝑡𝜕𝑥4𝑖𝐼𝑖𝑆𝑖𝐸1+𝑖𝑘𝑖𝐺𝑖𝜕4𝜓𝑖𝑥𝑖,𝑡𝜕𝑥2𝑖𝜕𝑡2+𝜕2𝜓𝑖𝑥𝑖,𝑡𝜕𝑡2+𝜌𝑖𝐼𝑖𝑘𝑖𝐺𝑖𝑆𝑖𝜕4𝜓𝑖𝑥𝑖,𝑡𝜕𝑡4=0(11) for the rotation 𝜓𝑖(𝑥𝑖,𝑡).

Now, assuming separation of variables in the form𝜓𝑖𝑥𝑖,𝑡=Ψ𝑖𝑥𝑖ej𝜔𝑡,𝑣𝑖𝑥𝑖,𝑡=𝑉𝑖𝑥𝑖ej𝜔𝑡,𝑢𝑖𝑥𝑖,𝑡=𝑈𝑖𝑥𝑖ej𝜔𝑡(12) with natural frequency 𝜔, (11), (10), and (9) can be, respectively, expressed as ordinary differential equations in the following form: 𝐸𝑖𝐼𝑖𝜌𝑖𝑆𝑖Ψ𝑖𝑥𝑖+𝜔2𝐼𝑖𝑆𝑖𝐸1+𝑖𝑘𝑖𝐺𝑖Ψ𝑖𝑥𝑖𝜔2Ψ𝑖𝑥𝑖+𝜔4𝜌𝑖𝐼𝑖𝑘𝑖𝐺𝑖𝑆𝑖Ψ𝑖𝑥𝑖=0,(13)𝐸𝑖𝐼𝑖𝜌𝑖𝑆𝑖𝑉𝑖𝑥𝑖+𝜔2𝐼𝑖𝑆𝑖𝐸1+𝑖𝑘𝑖𝐺𝑖𝑉𝑖𝑥𝑖𝜔2𝑉𝑖𝑥𝑖+𝜔4𝜌𝑖𝐼𝑖𝑘𝑖𝐺𝑖𝑆𝑖𝑉𝑖𝑥𝑖=0,(14)𝐸𝑖𝑆𝑖𝑈𝑖𝑥𝑖+𝜔2𝜌𝑖𝑆𝑖𝑈𝑖𝑥𝑖=0,(15) with unknown parameter 𝜔 and unknown functions Ψ𝑖(𝑥𝑖), 𝑉𝑖(𝑥𝑖) and 𝑈𝑖(𝑥𝑖).

Now, given the above 3×𝑛 relations, we need 6×𝑛 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 𝑂𝑖 depend on the joint variables 𝜃𝑖, 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 𝑂1, where three equations will be written, the inner joints 𝑂𝑖(𝑖=1,2,,𝑛) between link 𝑖 and link 𝑖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)𝑂1 section: 𝜓1𝑥1||,𝑡𝑥1=0𝑣=0,(16)1𝑥1||,𝑡𝑥1=0𝑢=0,(17)1𝑥1||,𝑡𝑥1=0=0,(18)(ii)𝑂𝑖 section (𝑖=2,,𝑛): 𝜓𝑖1𝑥𝑖1||,𝑡𝑥𝑖1=𝐿𝑖1𝜓𝑖𝑥𝑖||,𝑡𝑥𝑖=0𝑣=0,(19)𝑖1𝑥𝑖1||,𝑡𝑥𝑖1=𝐿𝑖1sin𝜃𝑖+𝑢𝑖1𝑥𝑖1||,𝑡𝑥𝑖1=𝐿𝑖1cos𝜃𝑖𝑢𝑖𝑥𝑖||,𝑡𝑥𝑖=0𝑣=0,𝑖1𝑥𝑖1||,𝑡𝑥𝑖1=𝐿𝑖1cos𝜃𝑖𝑢𝑖1𝑥𝑖1||,𝑡𝑥𝑖1=𝐿𝑖1sin𝜃𝑖𝑣𝑖𝑥𝑖||,𝑡𝑥𝑖=0𝐸=0,𝑖1𝐼𝑖1𝜕𝜓𝑖1𝑥𝑖1,𝑡𝜕𝑥𝑖1||||𝑥𝑖1=𝐿𝑖1+𝐽𝑖1𝜕2𝜓𝑖1𝑥𝑖1,𝑡𝜕𝑡2||||𝑥𝑖1=𝐿𝑖1𝐸𝑖𝐼𝑖𝜕𝜓𝑖𝑥𝑖,𝑡𝜕𝑥𝑖||||𝑥𝑖=0𝑘=0,𝑖1𝐺𝑖1𝑆𝑖1𝜕𝑣𝑖1𝑥𝑖1,𝑡𝜕𝑥𝑖1𝜓𝑖1𝑥𝑖1|||||,𝑡𝑥𝑖1=𝐿𝑖1+𝑀𝑖1𝜕2𝑣𝑖1𝑥𝑖1,𝑡𝜕𝑡2||||𝑥𝑖1=𝐿𝑖1𝐸𝑖𝑆𝑖𝜕𝑢𝑖𝑥𝑖,𝑡𝜕𝑥𝑖||||𝑥𝑖=0sin𝜃𝑖𝑘𝑖𝐺𝑖𝑆𝑖𝜕𝑣𝑖𝑥𝑖,𝑡𝜕𝑥𝑖𝜓𝑖𝑥𝑖|||||,𝑡𝑥𝑖=0cos𝜃𝑖𝐸=0,𝑖1𝑆𝑖1𝜕𝑢𝑖1(𝑥𝑖1,𝑡)𝜕𝑥𝑖1||||𝑥𝑖1=𝐿𝑖1+𝑀𝑖1𝜕2𝑢𝑖1(𝑥𝑖1,𝑡)𝜕𝑡2||||𝑥𝑖1=𝐿𝑖1𝐸𝑖𝑆𝑖𝜕𝑢𝑖(𝑥𝑖,𝑡)𝜕𝑥𝑖||||𝑥𝑖=0cos𝜃𝑖+𝑘𝑖𝐺𝑖𝑆𝑖𝜕𝑣𝑖(𝑥𝑖,𝑡)𝜕𝑥𝑖𝜓𝑖(𝑥𝑖||||,𝑡)𝑥𝑖=0sin𝜃𝑖=0,(20)(iii)EE section (end-effector) 𝐸𝑛𝐼𝑛𝜕𝜓𝑛𝑥𝑛,𝑡𝜕𝑥𝑛||||𝑥𝑛=𝐿𝑛+𝐽𝑛𝜕2𝜓𝑛𝑥𝑛,𝑡𝜕𝑡2||||𝑥𝑛=𝐿𝑛𝑘=0,𝑛𝐺𝑛𝑆𝑛𝜕𝑣𝑛𝑥𝑛,𝑡𝜕𝑥𝑛𝜓𝑛𝑥𝑛|||||,𝑡𝑥𝑛=𝐿𝑛+𝑀𝑛𝜕2𝑣𝑛𝑥𝑛,𝑡𝜕𝑡2||||𝑥𝑛=𝐿𝑛𝐸=0,𝑛𝑆𝑛𝜕𝑢𝑛𝑥𝑛,𝑡𝜕𝑥𝑛||||𝑥𝑛=𝐿𝑛+𝑀𝑛𝜕2𝑢𝑛𝑥𝑛,𝑡𝜕𝑡2||||𝑥𝑛=𝐿𝑛=0.(21)

Using solutions (12) into (16)–(21), the boundary conditions for the analyzed reference configuration of the manipulator reduce to

(i)𝑂1 section: Ψ1𝑥1||𝑥1=0𝑉=0,(22)1𝑥1||𝑥1=0𝑈=0,(23)1𝑥1||𝑥1=0=0,(24)(ii)𝑂𝑖 section (𝑖=2,,𝑛):

Ψ𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1Ψ𝑖(𝑥𝑖)||𝑥𝑖=0𝑉=0𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1sin𝜃𝑖+𝑈𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1cos𝜃𝑖𝑈𝑖(𝑥𝑖)||𝑥𝑖=0𝑉=0,𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1cos𝜃𝑖𝑈𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1sin𝜃𝑖𝑉𝑖(𝑥𝑖)||𝑥𝑖=0𝐸=0,𝑖1𝐼𝑖1Ψ𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1𝜔2𝐽𝑖1Ψ𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1𝐸𝑖𝐼𝑖Ψ𝑖(𝑥𝑖)||𝑥𝑖=0𝑘=0,𝑖1𝐺𝑖1𝑆𝑖1𝑉𝑖1(𝑥𝑖1)Ψ𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1𝜔2𝑀𝑖1𝑉𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1𝐸𝑖𝑆𝑖𝑈𝑖(𝑥𝑖)||𝑥𝑖=0sin𝜃𝑖𝑘𝑖𝐺𝑖𝑆𝑖𝑉𝑖(𝑥𝑖)Ψ𝑖(𝑥𝑖)||𝑥𝑖=0cos𝜃𝑖𝐸=0,𝑖1𝑆𝑖1𝑈𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1𝜔2𝑀𝑖1𝑈𝑖1(𝑥𝑖1)||𝑥𝑖1=𝐿𝑖1𝐸𝑖𝑆𝑖𝑈𝑖(𝑥𝑖)||𝑥𝑖=0cos𝜃𝑖,+𝑘𝑖𝐺𝑖𝑆𝑖𝑉𝑖(𝑥𝑖)Ψ𝑖(𝑥𝑖)||𝑥𝑖=0sin𝜃𝑖=0,(25)(iii)EE section (end-effector):𝐸𝑛𝐼𝑛Ψ𝑛(𝑥𝑛)||𝑥𝑛=𝐿𝑛𝜔2𝐽𝑛Ψ𝑛(𝑥𝑛)||𝑥𝑛=𝐿𝑛𝑘=0,𝑛𝐺𝑛𝑆𝑛𝑉𝑛(𝑥𝑛)Ψ𝑛(𝑥𝑛)||𝑥𝑛=𝐿𝑛𝜔2𝑀𝑛𝑉𝑛(𝑥𝑛)||𝑥𝑛=𝐿𝑛𝐸=0,𝑛𝑆𝑛𝑈𝑛(𝑥𝑛)||𝑥𝑛=𝐿𝑛𝜔2𝑀𝑛𝑈𝑛(𝑥𝑛)||𝑥𝑛=𝐿𝑛=0.(26)

Equations (13)–(15), expressed for all the 𝑛 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 𝑈𝑖(𝑥𝑖) satisfying the boundary conditions, it is possible to obtain natural frequencies and mode shapes of the system for an assigned set {𝜃𝑖𝑖=2,,𝑛}. To this end, let us consider (14) in the following form: 𝐸𝑖𝐼𝑖𝑉𝑖𝑥𝑖+𝜔2𝜌𝑖𝐼𝑖𝐸1+𝑖𝑘𝑖𝐺𝑖𝑉𝑖𝑥𝑖+𝜔2𝜔2𝜌2𝑖𝐼𝑖𝑘𝑖𝐺𝑖𝜌𝑖𝑆𝑖𝑉𝑖𝑥𝑖=0.(27)

From the characteristic equation associated to (27) it is possible to calculate the following two solutions for each 𝜔frequency:𝜆2𝑖1=𝜌𝑖𝐼𝑖𝜔2𝐸1+𝑖/𝑘𝑖𝐺𝑖+Ω2𝐸𝑖𝐼𝑖,(28)𝜆2𝑖2=𝜌𝑖𝐼𝑖𝜔2𝐸1+𝑖/𝑘𝑖𝐺𝑖Ω2𝐸𝑖𝐼𝑖,(29) where Ω is a positive or zero quantity defined byΩ=𝜔4𝜌2𝑖𝐼2𝑖𝐸1+𝑖𝑘𝑖𝐺i24𝐸𝑖𝐼𝑖𝜔2𝜔2𝜌2𝑖𝐼𝑖𝑘𝑖𝐺𝑖𝜌𝑖𝑆𝑖=𝜔4𝜌2𝑖𝐼2𝑖𝐸1𝑖𝑘𝑖𝐺i2+4𝜔2𝜌𝑖𝑆𝑖𝐸𝑖𝐼𝑖.(30)

It is observed that for 𝜔=[(𝑘𝑖𝐺𝑖𝑆𝑖)/(𝜌𝑖𝐼𝑖)]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 ±𝜆𝑖1 and two imaginary solutions ±𝑗𝜆𝑖2, whereas, for 𝜔>𝜔lim, four imaginary solutions ±𝑗𝜆𝑖1 and ±𝑗𝜆𝑖2 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Ψ𝑖𝑥𝑖=𝐴𝑖1𝜆cosh𝑖1𝑥𝑖+𝐴𝑖2𝜆sinh𝑖1𝑥𝑖+𝐴𝑖3𝜆cos𝑖2𝑥𝑖+𝐴𝑖4𝜆sin𝑖2𝑥𝑖,(31)𝑉𝑖𝑥𝑖=𝐴𝑖1𝜆cosh𝑖1𝑥𝑖+𝐴𝑖2𝜆sinh𝑖1𝑥𝑖+𝐴𝑖3𝜆cos𝑖2𝑥𝑖+𝐴𝑖4𝜆sin𝑖2𝑥𝑖,(32)𝑈𝑖𝑥𝑖=𝐵𝑖1cos𝜔𝑥𝑖𝑐𝑖+𝐵𝑖2sin𝜔𝑥𝑖𝑐𝑖,(33) where 𝑐𝑖=(𝐸𝑖/𝜌𝑖)1/2 is the well-known longitudinal wave speed.

The coefficients 𝐴𝑖𝑘 and 𝐴𝑖𝑘 (𝑘=1,,4) are not independent, but related through (8); in fact, substituting solutions (12) into it for 𝑣𝑖(𝑥𝑖,𝑡) and 𝜓𝑖(𝑥𝑖,𝑡) and using expressions (31) and (32), we obtain the following relations:𝐴𝑖1=𝑡𝑖1𝐴𝑖2,𝐴𝑖2=𝑡𝑖1𝐴𝑖1,𝐴𝑖3=𝑡𝑖2𝐴𝑖4,𝐴𝑖4=𝑡𝑖2𝐴𝑖3.(34)


Moreover, inserting (34) into (31), the final form for the solution Ψ𝑖(𝑥𝑖) is given byΨ𝑖𝑥𝑖=𝑡𝑖1𝐴𝑖1𝜆sinh𝑖1𝑥𝑖+𝑡𝑖1𝐴𝑖2𝜆cosh𝑖1𝑥𝑖𝑡𝑖2𝐴𝑖3𝜆sin𝑖2𝑥𝑖+𝑡𝑖2𝐴𝑖4𝜆cos𝑖2𝑥𝑖.(36)

Now, the expressions for Ψ𝑖(𝑥𝑖), 𝑉𝑖(𝑥𝑖), and 𝑈𝑖(𝑥𝑖) 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×𝑛 equations with 6×𝑛 unknown coefficients 𝐴𝑖𝑘 and 𝐵𝑖𝑧 (𝑖=1,,𝑛, 𝑧=1,2, and 𝑘=1,,4). These equations can be expressed in matrix form as 𝐌(𝜔)𝐜=𝟎,(37) where𝐴𝐜=11𝐴12𝐴13𝐴14𝐵11𝐵12𝐴𝑖1𝐴𝑖2𝐴𝑖3𝐴𝑖4𝐵𝑖1𝐵𝑖2𝐴𝑛1𝐴𝑛2𝐴𝑛3𝐴𝑛4𝐵𝑛1𝐵𝑛2T(38) is the 6×𝑛 vector of modal coefficients characterizing the eigenfunctions of the assigned differential eigenvalue problem, whereas𝐌𝐌(𝜔)=𝐎𝟏(𝜔)𝟑×𝟔[𝟎]𝟑×𝟔[𝟎]𝟑×𝟔𝐌𝐎𝟐,𝟏(𝜔)𝟔×𝟔𝐌𝐎𝟐,𝟐(𝜔)𝟔×𝟔[𝟎]𝟔×𝟔[𝟎]𝟔×𝟔[𝟎]𝟔×𝟔𝐌𝐎𝑖.𝑖1(𝜔)𝟔×𝟔𝐌𝐎𝑖,𝑖(𝜔)𝟔×𝟔[𝟎]𝟔×𝟔[𝟎]𝟔×𝟔[𝟎]𝟔×𝟔𝐌𝐎𝑛,𝑛1(𝜔)𝟔×𝟔𝐌𝐎𝑛,𝑛(𝜔)𝟔×𝟔[𝟎]𝟑×𝟔[𝟎]𝟑×𝟔𝐌𝐎𝐄𝐄(𝜔)𝟑×𝟔(39)

is a square matrix of order 6×𝑛 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, 𝐌𝐎𝟏(𝜔) and 𝐌𝐎𝐄𝐄(𝜔) are blocks corresponding with the two extreme sections of the structure 𝑂1 and EE, respectively; these blocks have dimensions 3×6 since both represent three boundary conditions and, respectively, concern the modal coefficients (𝐴11,𝐴12,𝐴13,𝐴14,𝐵11,𝐵12) and (𝐴𝑛1,𝐴𝑛2,𝐴𝑛3,𝐴𝑛4,𝐵𝑛1,𝐵𝑛2).

𝐌𝐎𝑖.𝑖1(𝜔) and 𝐌𝐎𝑖.𝑖(𝜔) are 6×6 blocks and can be interpreted as a matrix representation of the six boundary conditions at the inner blocked joint 𝑂𝑖 between link 𝑖-1 and link 𝑖, concerning the twelve modal coefficients (𝐴𝑠1,𝐴𝑠2,𝐴𝑠3,𝐴𝑠4,𝐵𝑠1,𝐵𝑠2), where 𝑠=𝑖1, 𝑖. The blocks of 𝐌(𝜔) and all their coefficients are explicitly reported in Appendix A.

The homogeneous system (37) admits nontrivial solutions for the unknown coefficients 𝐴𝑖𝑘 and 𝐵𝑖𝑧, (𝑖=1,,𝑛, 𝑧=1,2, and 𝑘=1,,4) once the determinant of the characteristic matrix 𝐌(𝜔) is settled to be zero:[]det𝐌(𝜔)=0.(40)

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 𝜔𝑗 (𝑗=1,2,), recognized as natural frequencies. Corresponding to each of these roots there are the eigenfunctions 𝑉𝑖(𝑗)(𝑥𝑖) and 𝑈𝑖(𝑗)(𝑥𝑖) that represent the 𝑗th axial-transverse mode shape, and the eigenfunction Ψ𝑖(𝑗)(𝑥𝑖) that represents the 𝑗th cross-section rotation mode shape. Based on relations (32), (33), and (36), the required modal coefficients 𝐴𝑖𝑘(𝑗) and 𝐵𝑖𝑧(𝑗) for the above-mentioned eigenfunctions are obtained by substituting the natural frequency 𝜔𝑗 into the characteristic matrix and solving the homogeneous system of (37).

Referring to the axial and transverse local components of deformation 𝑉𝑖(𝑗)(𝑥𝑖) and 𝑈𝑖(𝑗)(𝑥𝑖), let us suppose to represent the deformed shape of the multi-link manipulator in the global coordinate frame (𝑂0,𝑋0,𝑌0,𝑍0); for this scope, it is primarily observed that the homogeneous transformation matrix (43) settles the coordinate transformation between frames (𝑂𝑖,𝑋𝑖,𝑌𝑖,𝑍𝑖) and (𝑂𝑖1,𝑋𝑖1,𝑌𝑖1,𝑍𝑖1), for 𝑖=1,,𝑛, that is [15],𝐀𝑖𝑖1=𝐑𝑖𝑖1𝐝𝑖𝑖1𝟎𝑇1,(41) where 𝐑𝑖𝑖1 is the rotation matrix describing the orientation of frame 𝑖 relative to frame 𝑖1 and has the following form:𝐑𝑖𝑖1=cos𝜃𝑖sin𝜃𝑖0sin𝜃𝑖cos𝜃𝑖0001,(42) whilst 𝐝𝑖𝑖1=[𝐿𝑖100]𝑇 denotes the position of the origin 𝑂𝑖 of coordinate frame i relative to coordinate frame 𝑖1; by definition, we set 𝐿0 equal to zero. Now, the position of any point of the ith deformed link is expressed in the coordinate frame (𝑂𝑖,𝑋𝑖,𝑌𝑖,𝑍𝑖) through the homogeneous position vector 𝐩𝑖=𝑥𝑖+𝑈𝑖(𝑗)𝑥𝑖𝑉𝑖(𝑗)𝑥𝑖01.(43)

Finally, such a point is directly expressed in the global coordinate frame (𝑂0,𝑋0,𝑌0,𝑍0) through a succession of consecutive coordinate transformations which give the following global transformation:𝐩0=𝐀01𝐀12𝐀𝑖𝑖1𝐩𝑖.(44)

Relation (44) expresses the existing mapping between the coordinate reference frame (𝑂𝑖,𝑋𝑖,𝑌𝑖,𝑍𝑖) and the global reference frame (𝑂0,𝑋0,𝑌0,𝑍0), regarding the deformation components of the 𝑛 links, and it allows representation of the 𝑗th axial-transverse mode shape of the multi-link manipulator for the chosen configuration, after 𝑉𝑖(𝑗)(𝑥𝑖) and 𝑈𝑖(𝑗)(𝑥𝑖) 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 (𝑛=2) takes the following form:𝐌𝐌(𝜔)=𝐎𝟏(𝜔)𝟑×𝟔[𝟎]𝟑×𝟔𝐌𝐎𝟐,𝟏(𝜔)𝟔×𝟔𝐌𝐎𝟐,𝟐(𝜔)𝟔×𝟔[𝟎]𝟑×𝟔𝐌𝐎𝐄𝐄(𝜔)𝟑×𝟔.(45)

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] (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 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.

Table 1: Parameters of the flexible manipulator test.

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 𝐿1/ and 𝐿𝑇/ are used as reference parameters to characterize each simulation, 𝐿1 being the original value indicated in Table 1, 𝐿𝑇=𝐿1+𝐿2 the total manipulator length, and, finally, the referential thickness value. In order to derive modal data as a function of the slenderness, the two thicknesses 1 and 2 (Table 1) are multiplied in each case for a common thickening factor between 1 and 100; in any case, , indicated into the following numerical tables/simulations, refers to 2, which is always the minimum between the two actual thickness values.

Tables 26 report the first ten natural frequencies calculated through the three models when the slenderness ratios 𝐿1/ and 𝐿𝑇/ of the manipulator are, respectively, calculated using a common thickening factor equal to 1 (original structure), 10, 25, 50 and 100. In the tables, 𝑇 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 𝑘𝑖 (𝑖=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 Table 1; (𝐿1/=833.3, 𝐿𝑇/=1733).
Table 3: First ten natural frequencies (Hz); two-link manipulator; (𝐿1/=83.33, 𝐿𝑇/=173.3).
Table 4: First ten natural frequencies (Hz); two-link manipulator; (𝐿1/=33.33, 𝐿𝑇/=69.33).
Table 5: First ten natural frequencies (Hz); two-link manipulator; (𝐿1/=16.67, 𝐿𝑇/=34.67).
Table 6: First ten natural frequencies (Hz); two-link manipulator; (𝐿1/=8.33, 𝐿𝑇/=17.33).

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Δ𝑗%=frequency𝐸𝐵frequency𝑇frequency𝑇100,𝑗=1,,10,(46) where EB refers both to EB and EB.

A perusal of Tables 26 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 1 and 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 𝑇-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 𝐸𝐵 (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 (𝐿/>80 in Tables 2 and 3).

It is stressed, therefore, how the numerical comparison shows the ability of both the 𝑇-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 Tables 26 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 26. Both mechanical and electromechanical stiffnesses of each joint have been modeled as torsional springs having stiffness 𝑘𝑇1 for the shoulder joint and 𝑘𝑇2 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 26 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 𝑖 and posture) the prediction on the 𝑖th natural frequency is represented versus torsional stiffness 𝑘𝑇1 and 𝑘𝑇2 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.

From trends of Table 7, one can note that moving from rigid condition (i.e., 𝑘𝑇 = 1011 Nm/rad) toward the pinned boundary condition natural frequencies remain almost constant for each mode up to 𝑘𝑇=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 (𝐿1/=8.33). Significant frequency changes occur for 𝑘𝑇 lower than 106 Nm/rad although this is not the case for manipulator extremely slender (i.e., 𝐿/=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 25), 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 (𝑇-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 (Figure 3) whilst its mass and geometric properties are reported in Table 8; 𝑀𝑖 and  𝐽𝑖 (𝑖=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.
Figure 3: Rectangular hollow section of the adopted industrial manipulator.

The comparison between the 𝑇-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).

Figure 4: Test scenarios: (a) 𝜃2,3=45, 45°; (b) 𝜃2,3=45, 90°; (c) 𝜃2,3=90, −90°; (d) 𝜃2,3=90, −90°.

For each configuration, natural frequencies and mode shapes are derived for both exact models 𝑇 and EB by using the following system characteristic matrix:𝐌𝐌(𝜔)=𝐎𝟏(𝜔)𝟑×𝟔[𝟎]𝟑×𝟔[𝟎]𝟑×𝟔𝐌𝐎𝟐,𝟏(𝜔)𝟔×𝟔𝐌𝐎𝟐,𝟐(𝜔)𝟔×𝟔[𝟎]𝟔×𝟔[𝟎]𝟔×𝟔𝐌𝐎𝟑,𝟐(𝜔)𝟔×𝟔𝐌𝐎𝟑,𝟑(𝜔)𝟔×𝟔[𝟎]𝟑×𝟔[𝟎]𝟑×𝟔𝐌𝐎𝐄𝐄(𝜔)𝟑×𝟔,(47) 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.

Table 9: First ten natural frequencies (Hz); see Table 7; 𝜃2,3=45,90.
Table 10: First ten natural frequencies (Hz); see Table 7; 𝜃2,3=90,90.
Table 11: First ten natural frequencies (Hz); see Table 7; 𝜃2,3=90,90.
Figure 5: (a) Modes (Hz) of an industrial manipulator (see Table 8). (FEM: ABAQUS 6.6, B23, B22) left column: EB-model; right column: 𝑇-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: 𝑇-model: 𝜃2=𝜃3=45.

Finite element analysis of the structure is carried out by using ABAQUS R. 6.6, using B23 beam elements for the EB-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 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.

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.


A. Blocks of the Characteristic Matrix (Timoshenko Beam Model)


B. Blocks of the Characteristic Matrix (Euler-Bernoulli Beam Model)



This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.


  1. F. Bellezza, L. Lanari, and G. Ulivi, “Exact modeling of the slewing flexible link,” in Proceedings of the IEEE International Conference on Robotics and Automation, pp. 734–739, Cincinnati, Ohio, USA, May 1999.
  2. M. W. D. White and G. R. Heppler, “Timoshenko model of a flexible slewing link,” in Proceedings of the American Control Conference, pp. 2815–2819, Seattle, Wash, USA, June 1995.
  3. Y. Chaolan, H. Jiazhen, and C. Guoping, “Modeling study of a flexible hub-beam system with large motion and with considering the effect of shear deformation,” Journal of Sound and Vibration, vol. 295, no. 1-2, pp. 282–293, 2006. View at Publisher · View at Google Scholar
  4. J. C. Ower and J. Van de Vegte, “Classical control design for a flexible manipulator: modeling and control system design,” IEEE Journal of Robotics and Automation, vol. 3, no. 5, pp. 485–489, 1987.
  5. P. Tomei and A. Tornambe, “Approximate modeling of robots having elastic links,” IEEE Transactions on Systems, Man and Cybernetics, vol. 18, no. 5, pp. 831–840, 1988. View at Publisher · View at Google Scholar
  6. A. De Luca and B. Siciliano, “Closed-form dynamic model of planar multilink lightweight robots,” IEEE Transactions on Systems, Man and Cybernetics, vol. 21, no. 4, pp. 826–839, 1991. View at Publisher · View at Google Scholar · View at MathSciNet
  7. W. Chen, “Dynamic modeling of multi-link flexible robotic manipulators,” Computers and Structures, vol. 79, no. 2, pp. 183–195, 2001. View at Publisher · View at Google Scholar
  8. H. H. Lee, “New dynamic modeling of flexible-link robots,” Journal of Dynamic Systems, Measurement and Control, vol. 127, no. 2, pp. 307–309, 2005. View at Publisher · View at Google Scholar
  9. X. Zhang, W. Xu, S. S. Nair, and V. S. Chellaboina, “PDE modeling and control of a flexible two-link manipulator,” IEEE Transactions on Control Systems Technology, vol. 13, no. 2, pp. 301–312, 2005. View at Publisher · View at Google Scholar
  10. M. Pascal, “Some Open Problems in Dynamic Analysis of Flexible Multibody Systems,” Multibody System Dynamics, vol. 5, no. 4, pp. 315–334, 2001. View at Publisher · View at Google Scholar
  11. R. I. Milford and S. F. Asokanthan, “Configuration dependent eigenfrequencies for a two-link flexible manipulator: experimental verification,” Journal of Sound and Vibration, vol. 222, no. 2, pp. 191–207, 1999.
  12. C. di Castri, A. Messina, and G. Reina, “Modal analysis of a two link flexible manipulator,” in Proceedings of the RAAD 18th International Workshop of Robotics in Alpe-Adria-Danube Region, Brasov, Romania, May 2009.
  13. S. Timoshenko, Vibration Problems in Engineering, D. Van Nostrand Company, Inc., New York, NY, USA, 2nd edition, 1937.
  14. K. F. Graff, Wave Motion in Elastic Solids, Dover Publications, New York, NY, USA, 1991.
  15. B. Siciliano and O. Khatib, Springer Handbook of Robotics, Springer, Berlin, Germany, 2008.
  16. G. R. Cowper, “The shear coefficient in Timoshenko’s beam theory,” Journal of Applied Mechanics, vol. 33, pp. 335–340, 1966.
  17. Simulia, ABAQUS Version 6.6 HTML Documentation, 2006.