Vibration Analysis of Multilink Manipulators Based on Timoshenko Beam Theory

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.

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 𝑂₁ (we say 𝜃₁) always has null value. For a chosen configuration, the manipulator is fixed at point 𝑂₁ and internal blocked joints at points 𝑂_𝑖 avoid relative rotations between links.

The coordinate system (𝑂₀,𝑋₀,𝑌₀,𝑍₀) 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 𝑋₁ and 𝑋₀ axes, that is, 𝜃₁, is always null for the above-mentioned reasons, so that frame (𝑂₀,𝑋₀,𝑌₀,𝑍₀) and frame (𝑂₁,𝑋₁,𝑌₁,𝑍₁) 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𝜕𝑀_𝑖𝑥_𝑖,𝑡𝜕𝑥_𝑖+𝑇_𝑖𝑥_𝑖,𝑡=𝜌_𝑖𝐼_𝑖𝜕²𝜓_𝑖𝑥_𝑖,𝑡𝜕𝑡²𝜕𝑇_𝑖,𝑥_𝑖,𝑡𝜕𝑥_𝑖−𝜌_𝑖𝑆_𝑖𝜕²𝑣_𝑖𝑥_𝑖,𝑡𝜕𝑡²𝜕𝑁_𝑖=−𝑞_𝑖𝑥_𝑖,𝑡,𝑥_𝑖,𝑡𝜕𝑥_𝑖−𝜌_𝑖𝑆_𝑖𝜕²𝑢_𝑖𝑥_𝑖,𝑡𝜕𝑡²=−𝑝_𝑖𝑥_𝑖,𝑡.(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𝑖+𝑘_𝑖𝐺_𝑖𝑆_𝑖𝜕𝑣^𝑖𝑥_𝑖,𝑡𝜕𝑥_𝑖=𝜌_𝑖−𝜓_𝑖𝑥_𝑖,𝑡𝐼_𝑖𝜕²𝜓_𝑖𝑥_𝑖,𝑡𝜕𝑡²,(4)𝑘_𝑖𝐺_𝑖𝑆_𝑖𝜕𝜓^𝑖𝑥_𝑖,𝑡𝜕𝑥_𝑖−𝜕²𝑣_𝑖𝑥_𝑖,𝑡𝜕𝑥^2𝑖+𝜌^𝑖𝑆_𝑖𝜕²𝑣_𝑖𝑥_𝑖,𝑡𝜕𝑡²=𝑞_𝑖𝑥_𝑖,𝑡,(5)𝐸_𝑖𝑆_𝑖𝜕²𝑢_𝑖𝑥_𝑖,𝑡𝜕𝑥^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𝑖+𝑘_𝑖𝐺_𝑖𝑆_𝑖𝜕𝑣^𝑖𝑥_𝑖,𝑡𝜕𝑥_𝑖=𝜌_𝑖−𝜓_𝑖𝑥_𝑖,𝑡𝐼_𝑖𝜕²𝜓_𝑖𝑥_𝑖,𝑡𝜕𝑡²𝑘,(7)_𝑖𝐺_𝑖𝑆_𝑖𝜕𝜓^𝑖𝑥_𝑖,𝑡𝜕𝑥_𝑖−𝜕²𝑣_𝑖𝑥_𝑖,𝑡𝜕𝑥^2𝑖+𝜌^𝑖𝑆_𝑖𝜕²𝑣_𝑖𝑥_𝑖,𝑡𝜕𝑡²𝐸=0,(8)_𝑖𝑆_𝑖𝜕²𝑢_𝑖𝑥_𝑖,𝑡𝜕𝑥^2𝑖−𝜌_𝑖𝑆_𝑖𝜕²𝑢_𝑖𝑥_𝑖,𝑡𝜕𝑡²=0.(9)

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

Now, assuming separation of variables in the form𝜓_𝑖𝑥_𝑖,𝑡=Ψ𝑢_𝑖𝑖𝑥_𝑖e^j𝜔𝑡,𝑣_𝑖𝑥_𝑖,𝑡=𝑉_𝑖𝑥_𝑖e^j𝜔𝑡,𝑥_𝑖,𝑡=𝑈_𝑖𝑥_𝑖e^j𝜔𝑡(12) with natural frequency 𝜔, (11), (10), and (9) can be, respectively, expressed as ordinary differential equations in the following form: 𝐸_𝑖𝐼_𝑖𝜌_𝑖𝑆_𝑖Ψ^{𝑖}𝑥_𝑖+𝜔²𝐼_𝑖𝑆_𝑖1+𝑘𝐸_𝑖𝐺^𝑖_𝑖Ψ^𝑖𝑥_𝑖−𝜔²+𝜔₄Ψ_𝑖𝑥_𝑖𝜌_𝑖𝐼_𝑖𝑘_𝑖𝐺_𝑖𝑆_𝑖Ψ_𝑖𝑥_𝑖=0,(13)𝐸_𝑖𝐼_𝑖𝜌_𝑖𝑆_𝑖𝑉^{𝑖}𝑥_𝑖+𝜔²𝐼_𝑖𝑆_𝑖1+𝑘𝐸_𝑖𝐺^𝑖_𝑖𝑉^𝑖𝑥_𝑖−𝜔²+𝜔₄𝑉_𝑖𝑥_𝑖𝜌_𝑖𝐼_𝑖𝑘_𝑖𝐺_𝑖𝑆_𝑖𝑉_𝑖𝑥_𝑖=0,(14)𝐸_𝑖𝑆_𝑖𝑈_𝑖𝑥_𝑖+𝜔²𝜌_𝑖𝑆_𝑖𝑈_𝑖𝑥_𝑖=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 𝑂₁, 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)𝑂₁ section: 𝜓𝑣₁𝑥₁,𝑡||_𝑥1=0=0,(16)₁𝑢𝑥₁,𝑡||_𝑥1=0=0,(17)₁𝑥₁,𝑡||_𝑥1=0=0,(18)(ii)𝑂_𝑖 section (𝑖=2,…,𝑛): 𝜓𝑣_𝑖−1𝑥_𝑖−1,𝑡||_{𝑥𝑖−1=𝐿𝑖−1}−𝜓_𝑖𝑥_𝑖,𝑡||_𝑥𝑖=0=0,(19)_𝑖−1𝑥_𝑖−1,𝑡||_{𝑥𝑖−1=𝐿𝑖−1}sin𝜃_𝑖+𝑢_𝑖−1𝑥_𝑖−1,𝑡||_{𝑥𝑖−1=𝐿𝑖−1}cos𝜃_𝑖𝑣−𝑢_𝑖𝑥_𝑖,𝑡||_𝑥𝑖=0=0,_𝑖−1𝑥_𝑖−1,𝑡||_{𝑥𝑖−1=𝐿𝑖−1}cos𝜃_𝑖−𝑢_𝑖−1𝑥_𝑖−1,𝑡||_{𝑥𝑖−1=𝐿𝑖−1}sin𝜃_𝑖𝐸−𝑣_𝑖𝑥_𝑖,𝑡||_𝑥𝑖=0=0,_𝑖−1𝐼_𝑖−1𝜕𝜓_𝑖−1𝑥_𝑖−1,𝑡𝜕𝑥_𝑖−1||||_{𝑥𝑖−1=𝐿𝑖−1}+𝐽_𝑖−1𝜕²𝜓_𝑖−1𝑥_𝑖−1,𝑡𝜕𝑡²||||_{𝑥𝑖−1=𝐿𝑖−1}−𝐸_𝑖𝐼_𝑖𝜕𝜓_𝑖𝑥_𝑖,𝑡𝜕𝑥_𝑖||||_𝑥𝑖=0𝑘=0,_𝑖−1𝐺_𝑖−1𝑆_𝑖−1𝜕𝑣^𝑖−1𝑥_𝑖−1,𝑡𝜕𝑥_𝑖−1−𝜓_𝑖−1𝑥_𝑖−1,𝑡|||||_{𝑥𝑖−1=𝐿𝑖−1}+𝑀_𝑖−1𝜕²𝑣_𝑖−1𝑥_𝑖−1,𝑡𝜕𝑡²||||_{𝑥𝑖−1=𝐿𝑖−1}−𝐸_𝑖𝑆_𝑖𝜕𝑢_𝑖𝑥_𝑖,𝑡𝜕𝑥_𝑖||||_𝑥𝑖=0sin𝜃_𝑖−𝑘_𝑖𝐺_𝑖𝑆_𝑖𝜕𝑣^𝑖𝑥_𝑖,𝑡𝜕𝑥_𝑖𝐸−𝜓_𝑖𝑥_𝑖,𝑡|||||_𝑥𝑖=0cos𝜃_𝑖=0,_𝑖−1𝑆_𝑖−1𝜕𝑢_𝑖−1(𝑥_𝑖−1,𝑡)𝜕𝑥_𝑖−1||||_{𝑥𝑖−1=𝐿𝑖−1}+𝑀_𝑖−1𝜕²𝑢_𝑖−1(𝑥_𝑖−1,𝑡)𝜕𝑡²||||_{𝑥𝑖−1=𝐿𝑖−1}−𝐸_𝑖𝑆_𝑖𝜕𝑢_𝑖(𝑥_𝑖,𝑡)𝜕𝑥_𝑖||||_𝑥𝑖=0cos𝜃_𝑖+𝑘_𝑖𝐺_𝑖𝑆_𝑖𝜕𝑣^𝑖(𝑥_𝑖,𝑡)𝜕𝑥_𝑖−𝜓_𝑖(𝑥_𝑖,𝑡)||||_𝑥𝑖=0sin𝜃_𝑖=0,(20)(iii)EE section (end-effector) 𝐸_𝑛𝐼_𝑛𝜕𝜓_𝑛𝑥_𝑛,𝑡𝜕𝑥_𝑛||||_{𝑥𝑛=𝐿𝑛}+𝐽_𝑛𝜕²𝜓_𝑛𝑥_𝑛,𝑡𝜕𝑡²||||_{𝑥𝑛=𝐿𝑛}𝑘=0,_𝑛𝐺_𝑛𝑆_𝑛𝜕𝑣^𝑛𝑥_𝑛,𝑡𝜕𝑥_𝑛−𝜓_𝑛𝑥_𝑛,𝑡|||||_{𝑥𝑛=𝐿𝑛}+𝑀_𝑛𝜕²𝑣_𝑛𝑥_𝑛,𝑡𝜕𝑡²||||_{𝑥𝑛=𝐿𝑛}𝐸=0,_𝑛𝑆_𝑛𝜕𝑢_𝑛𝑥_𝑛,𝑡𝜕𝑥_𝑛||||_{𝑥𝑛=𝐿𝑛}+𝑀_𝑛𝜕²𝑢_𝑛𝑥_𝑛,𝑡𝜕𝑡²||||_{𝑥𝑛=𝐿𝑛}=0.(21)

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

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

𝑉Ψ_𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}−Ψ_𝑖(𝑥_𝑖)||_𝑥𝑖=0=0_𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}sin𝜃_𝑖+𝑈_𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}cos𝜃_𝑖𝑉−𝑈_𝑖(𝑥_𝑖)||_𝑥𝑖=0=0,_𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}cos𝜃_𝑖−𝑈_𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}sin𝜃_𝑖𝐸−𝑉_𝑖(𝑥_𝑖)||_𝑥𝑖=0=0,_𝑖−1𝐼_𝑖−1Ψ^𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}−𝜔²𝐽_𝑖−1Ψ_𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}𝑘−𝐸_𝑖𝐼_𝑖Ψ_𝑖(𝑥_𝑖)||_𝑥𝑖=0=0,_𝑖−1𝐺_𝑖−1𝑆_𝑖−1𝑉^𝑖−1(𝑥_𝑖−1)−Ψ_𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}−𝜔²𝑀_𝑖−1𝑉_𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}−𝐸_𝑖𝑆_𝑖𝑈^𝑖(𝑥_𝑖)||_𝑥𝑖=0sin𝜃_𝑖𝐸−𝑘_𝑖𝐺_𝑖𝑆_𝑖𝑉^𝑖(𝑥_𝑖)−Ψ_𝑖(𝑥_𝑖)||_𝑥𝑖=0cos𝜃_𝑖=0,_𝑖−1𝑆_𝑖−1𝑈^𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}−+𝑘𝜔²_𝑖𝑀_𝑖−1𝑈_𝑖−1(𝑥_𝑖−1)||_{𝑥𝑖−1=𝐿𝑖−1}−𝐸_𝑖𝑆_𝑖𝑈^𝑖(𝑥_𝑖)||_𝑥𝑖=0cos𝜃_𝑖,𝐺_𝑖𝑆_𝑖𝑉^𝑖(𝑥_𝑖)−Ψ_𝑖(𝑥_𝑖)||_𝑥𝑖=0sin𝜃_𝑖=0,(25)(iii)EE section (end-effector):𝑘𝐸_𝑛𝐼_𝑛Ψ^𝑛(𝑥_𝑛)||_{𝑥𝑛=𝐿𝑛}−𝜔²𝐽_𝑛Ψ_𝑛(𝑥_𝑛)||_{𝑥𝑛=𝐿𝑛}=0,_𝑛𝐺_𝑛𝑆_𝑛𝐸𝑉^𝑛(𝑥_𝑛)−Ψ_𝑛(𝑥_𝑛)||_{𝑥𝑛=𝐿𝑛}−𝜔²𝑀_𝑛𝑉_𝑛(𝑥_𝑛)||_{𝑥𝑛=𝐿𝑛}=0,_𝑛𝑆_𝑛𝑈^𝑛(𝑥_𝑛)||_{𝑥𝑛=𝐿𝑛}−𝜔²𝑀_𝑛𝑈_𝑛(𝑥_𝑛)||_{𝑥𝑛=𝐿𝑛}=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: 𝐸_𝑖𝐼_𝑖𝑉^{𝑖}𝑥_𝑖+𝜔²𝜌_𝑖𝐼_𝑖+𝜔₂1+𝑘𝐸_𝑖𝐺^𝑖_𝑖𝑉^𝑖𝑥_𝑖𝜔²𝜌^2𝑖𝐼_𝑖𝑘_𝑖𝐺_𝑖−𝜌_𝑖𝑆_𝑖𝑉^𝑖𝑥_𝑖=0.(27)

From the characteristic equation associated to (27) it is possible to calculate the following two solutions for each 𝜔frequency:𝜆^2𝑖1=−𝜌^𝑖𝐼_𝑖𝜔²1+𝐸_𝑖/𝑘_𝑖𝐺_𝑖+√Ω2𝐸_𝑖𝐼_𝑖,(28)𝜆^2𝑖2=−𝜌^𝑖𝐼_𝑖𝜔²1+𝐸_𝑖/𝑘_𝑖𝐺_𝑖−√Ω2𝐸_𝑖𝐼_𝑖,(29) where Ω is a positive or zero quantity defined byΩ=𝜔⁴𝜌^2𝑖𝐼^2𝑖1+𝑘𝐸_𝑖𝐺^𝑖_i²−4𝐸_𝑖𝐼_𝑖𝜔²𝜔²𝜌^2𝑖𝐼_𝑖𝑘_𝑖𝐺_𝑖−𝜌_𝑖𝑆_𝑖=𝜔₄𝜌^2𝑖𝐼^2𝑖1−𝑘𝐸_𝑖𝐺^𝑖_i²+4𝜔²𝜌_𝑖𝑆_𝑖𝐸_𝑖𝐼_𝑖.(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Ψ_𝑖𝑥_𝑖=+𝐴_𝑖1cosh𝜆_𝑖1𝑥_𝑖+𝐴_𝑖2sinh𝜆_𝑖1𝑥_𝑖𝐴_𝑖3cos𝜆_𝑖2𝑥_𝑖+𝐴_𝑖4sin𝜆_𝑖2𝑥_𝑖,(31)𝑉_𝑖𝑥_𝑖=+𝐴𝐴_𝑖1c_𝑖3osh𝜆_𝑖1𝑥_𝑖+𝐴_𝑖2sinh𝜆_𝑖1𝑥_𝑖cos𝜆_𝑖2𝑥_𝑖+𝐴_𝑖4sin𝜆_𝑖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)

where𝑡_𝑖1=1𝜆_𝑖1𝜆^2𝑖1+𝜔²𝜌_𝑖𝑡_𝑖2𝑘_𝑖𝐺_𝑖,=1𝜆_𝑖2𝜆^2𝑖2−𝜔²𝜌_𝑖𝑘_𝑖𝐺_𝑖.(35)

Moreover, inserting (34) into (31), the final form for the solution Ψ_𝑖(𝑥_𝑖) is given byΨ_𝑖𝑥_𝑖=−𝑡_𝑖1𝑡𝐴_𝑖1𝑖2sinh𝜆_𝑖1𝑥_𝑖+𝑡_𝑖1𝐴_𝑖2cosh𝜆_𝑖1𝑥_𝑖𝐴_𝑖3sin𝜆_𝑖2𝑥_𝑖+𝑡_𝑖2𝐴_𝑖4cos𝜆_𝑖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𝐜=𝐴₁₁𝐴₁₂𝐴₁₃𝐴₁₄𝐵₁₁𝐵₁₂⋯𝐴_𝑖1𝐴_𝑖2𝐴_𝑖3𝐴_𝑖4𝐵_𝑖1𝐵_𝑖2⋯𝐴_𝑛1𝐴_𝑛2𝐴_𝑛3𝐴_𝑛4𝐵_𝑛1𝐵_𝑛2^T(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 𝑂₁ and EE, respectively; these blocks have dimensions 3×6 since both represent three boundary conditions and, respectively, concern the modal coefficients (𝐴₁₁,𝐴₁₂,𝐴₁₃,𝐴₁₄,𝐵₁₁,𝐵₁₂) 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 (𝑂₀,𝑋₀,𝑌₀,𝑍₀); 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𝜃_𝑖−sin𝜃_𝑖00co0s𝜃_𝑖10⎞⎟⎟⎟⎠,(42) whilst 𝐝^𝑖𝑖−1=[^𝐿_𝑖−1⁰⁰]^𝑇 denotes the position of the origin 𝑂_𝑖 of coordinate frame i relative to coordinate frame 𝑖−1; by definition, we set 𝐿₀ 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 (𝑂₀,𝑋₀,𝑌₀,𝑍₀) through a succession of consecutive coordinate transformations which give the following global transformation:𝐩⁰=𝐀⁰¹𝐀¹²⋯𝐀^𝑖𝑖−1𝐩^𝑖.(44)

Relation (44) expresses the existing mapping between the coordinate reference frame (𝑂_𝑖,𝑋_𝑖,𝑌_𝑖,𝑍_𝑖) and the global reference frame (𝑂₀,𝑋₀,𝑌₀,𝑍₀), 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 𝐿₁/ℎ and 𝐿_𝑇/ℎ are used as reference parameters to characterize each simulation, 𝐿₁ being the original value indicated in Table 1, 𝐿_𝑇=𝐿₁+𝐿₂ 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 ℎ₁ and ℎ₂ (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 ℎ₂, which is always the minimum between the two actual thickness values.

Tables 2–6 report the first ten natural frequencies calculated through the three models when the slenderness ratios 𝐿₁/ℎ 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; (𝐿₁/ℎ=833.3, 𝐿_𝑇/ℎ=1733).

Table 3 

First ten natural frequencies (Hz); two-link manipulator; (𝐿₁/ℎ=83.33, 𝐿_𝑇/ℎ=173.3).

Table 4 

First ten natural frequencies (Hz); two-link manipulator; (𝐿₁/ℎ=33.33, 𝐿_𝑇/ℎ=69.33).

Table 5 

First ten natural frequencies (Hz); two-link manipulator; (𝐿₁/ℎ=16.67, 𝐿_𝑇/ℎ=34.67).

Table 6 

First ten natural frequencies (Hz); two-link manipulator; (𝐿₁/ℎ=8.33, 𝐿_𝑇/ℎ=17.33).

For each table, three nominal configurations of the two-link manipulator are considered, that is, when 𝜃₂ 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 2–6 reveals the following interesting findings.

Above any consideration, the model based on Timoshenko theory herein introduced always provides lower frequencies than the counterparts predicted by Milford and Asokanthan [11] and di Castri et al. [12]. This can be verified especially when thicker beams are taken into account (100, 50 and 25 times the original ℎ₁ and ℎ₂ 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 𝜃₂); 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 𝜃₂ 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 2–6 have been computed through analytical formulations which assume rigid joints only. In particular, this can be easily verified for the matrix formulation (39) herein proposed by simply noting that geometric boundary conditions (16) and (19) express perfect cantilevered relations for all the joints.

In order to investigate the effect on modal predictions when the above assumptions are not ensured, simulations have been carried out by employing apposite finite element models of the two-link manipulator of Tables 2–6. Both mechanical and electromechanical stiffnesses of each joint have been modeled as torsional springs having stiffness 𝑘_𝑇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 2–6 has been performed; graphical representations of the obtained results are reported by Table 7. Rows and columns represent mode numbers and postures, respectively. Inside each rectangle (for chosen mode 𝑖 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 (10¹¹ 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., 𝑘_𝑇 = 10¹¹ Nm/rad) toward the pinned boundary condition natural frequencies remain almost constant for each mode up to 𝑘_𝑇∼=10⁶ 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 (𝐿₁/ℎ=8.33). Significant frequency changes occur for 𝑘_𝑇 lower than 10⁶ 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 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 (𝑇-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 𝜃₂ and 𝜃₃ 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).

(a)

(b)

(c)

(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^∘.

(a)

(b)

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.

Appendices

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

𝐌_{𝐎𝑖,𝑖−𝟏}(𝜔)=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝^{𝑡i−1,1𝜆sinh𝑖−1,1𝐿𝑖−1𝑡𝑖−1,1𝜆cosh𝑖−1,1𝐿𝑖−1𝜆cosh𝑖−1,1𝐿𝑖−1𝜃sin𝑖𝜆sinh𝑖−1,1𝐿𝑖−1𝜃sin𝑖𝜆cosh𝑖−1,1𝐿𝑖−1𝜃cos𝑖𝜆sinh𝑖−1,1𝐿𝑖−1𝜃cos𝑖𝐸𝑖−1𝐼𝑖−1𝜆𝑖−1,1𝑡𝑖−1,1𝜆cosh𝑖−1,1L𝑖−1𝐸𝑖−1𝐼𝑖−1𝜆𝑖−1,1𝑡𝑖−1,1𝜆sinh𝑖−1,1𝐿𝑖−1−𝜔}_{2𝐽𝑖−1𝑡𝑖−1,1𝜆sinh𝑖−1,1𝐿𝑖−1−𝜔2𝐽𝑖−1𝑡𝑖−1,1𝜆cosh𝑖−1,1𝐿𝑖−1𝑘𝑖−1𝐺𝑖−1𝑆𝑖−1𝜆𝑖−1,1−𝑡𝑖−1,1𝜆sinh𝑖−1,1𝐿𝑖−1𝑘𝑖−1𝐺𝑖−1𝑆𝑖−1𝜆𝑖−1,1−𝑡𝑖−1,1𝜆cosh𝑖−1,1𝐿𝑖−1−𝜔2𝑀𝑖−1𝜆cosh𝑖−1,1𝐿𝑖−1−𝜔2𝑀𝑖−1𝜆sinh𝑖−1,1𝐿𝑖−100−𝑡𝑖−1,2𝜆sin𝑖−1,2𝐿𝑖−1𝑡𝑖−1,2𝜆cos𝑖−1,2𝐿𝑖−1𝜆cos𝑖−1,2𝐿𝑖−1𝜃sin𝑖𝜆sin𝑖−1,2𝐿𝑖−1𝜃sin𝑖𝜆cos𝑖−1,2𝐿𝑖−1𝜃cos𝑖𝜆sin𝑖−1,2𝐿𝑖−1𝜃cos𝑖−𝐸𝑖−1𝐼𝑖−1𝜆𝑖−1,2𝑡𝑖−1,2𝜆cos𝑖−1,2𝐿𝑖−1−𝐸𝑖−1𝐼𝑖−1𝜆𝑖−1,2𝑡𝑖−1,2𝜆sin𝑖−1,2𝐿𝑖−1+𝜔2𝐽𝑖−1𝑡𝑖−1,2𝜆sin𝑖−1,2𝐿𝑖−1−𝜔2𝐽𝑖−1𝑡𝑖−1,2𝜆cos𝑖−1,2𝐿𝑖−1𝑘𝑖−1𝐺𝑖−1𝑆𝑖−1𝑡𝑖−1,2−𝜆𝑖−1,2𝜆sin𝑖−1,2𝐿𝑖−1𝑘𝑖−1𝐺𝑖−1𝑆𝑖−1𝜆𝑖−1,2−𝑡𝑖−1,2𝜆cos𝑖−1,2𝐿𝑖−1−𝜔2𝑀𝑖−1𝜆cos𝑖−1,2𝐿𝑖−1−𝜔2𝑀𝑖−1𝜆sin𝑖−1,2𝐿𝑖−1𝜔0000cos𝑐𝑖−1𝐿𝑖−1𝜃cos𝑖𝜔sin𝑐𝑖−1𝐿𝑖−1𝜃cos𝑖𝜔−cos𝑐𝑖−1𝐿𝑖−1𝜃sin𝑖𝜔−sin𝑐𝑖−1𝐿𝑖−1𝜃sin𝑖0000−𝐸𝑖−1𝑆𝑖−1𝜔𝑐𝑖−1𝜔sin𝑐𝑖−1𝐿𝑖−1𝐸𝑖−1𝑆𝑖−1𝜔𝑐𝑖−1𝜔cos𝑐𝑖−1𝐿𝑖−1−𝜔2𝑀𝑖−1𝜔cos𝑐𝑖−1𝐿𝑖−1−𝜔2𝑀𝑖−1𝜔sin𝑐𝑖−1𝐿𝑖−1}⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠,(A.1)𝐌_{𝐎𝑖,𝑖}(𝜔)=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝0−𝑡^𝑖1−𝐸_𝑖−0100−001−00𝑡_𝑖2−001000𝐼0_𝑖𝜆_𝑖1𝑡_𝑖1−𝑘_𝑖0𝐸_𝑖𝐼_𝑖𝜆_𝑖2𝑡_𝑖2000𝐺_𝑖𝑆_𝑖𝜆_𝑖1−𝑡_𝑖1cos𝜃_𝑖0−𝑘_𝑖𝐺_𝑖𝑆_𝑖𝜆_𝑖2−𝑡_𝑖2cos𝜃_𝑖0−𝐸_𝑖𝑆_𝑖0𝑘_𝑖𝜔𝑐_𝑖sin𝜃^𝑖𝐺_𝑖𝑆_𝑖𝜆_𝑖1−𝑡_𝑖1sin𝜃_𝑖0𝑘_𝑖𝐺_𝑖𝑆_𝑖𝜆_𝑖2−𝑡_𝑖2sin𝜃_𝑖0−𝐸_𝑖𝑆_𝑖(A𝜔𝑐_𝑖.2)cos𝜃^𝑖⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠,𝐌_𝐎𝟏(𝜔)=⎛⎜⎜⎜⎝0𝑡¹¹1000100𝑡00₁₂010000⎞⎟⎟⎟⎠,(A.3)𝐌_𝐎𝐄𝐄(𝜔)=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝𝐸^𝑛−𝜔𝐼_𝑛2𝜆_𝑛1cosh𝜆_𝑛1𝐿_𝑛𝐸_𝑛𝐼_𝑛𝜆_𝑛1sinh𝜆_𝑛1𝐿_𝑛𝑘_𝑛𝐽_𝑛𝑡_𝑛1sinh𝜆_𝑛1𝐿_𝑛−𝜔²𝐽_𝑛𝑡_𝑛1cosh𝜆_𝑛1𝐿_𝑛𝐺_𝑛𝑆−_𝑛𝜔₂𝜆_𝑛1−𝑡_𝑛1sinh𝜆_𝑛1𝐿_𝑛𝑘_𝑛𝐺_𝑛𝑆_𝑛𝜆_𝑛1−𝑡_𝑛1cosh𝜆_𝑛1𝐿_𝑛−𝐸_𝑛𝑀_𝑛co0sh𝜆_𝑛1𝐿_𝑛−𝜔²𝑀_𝑛sin0h𝜆_𝑛1𝐿_𝑛+𝜔₂𝐼_𝑛𝜆_𝑛2cos𝜆_𝑛2𝐿_𝑛−𝐸_𝑛𝐼_𝑛𝜆_𝑛2sin𝜆_𝑛2𝐿_𝑛𝑘_𝑛𝐽_𝑛𝑡_𝑛2sinh𝜆_𝑛2𝐿_𝑛−𝜔²𝐽_𝑛𝑡_𝑛2cos𝜆_𝑛2𝐿_𝑛𝐺_𝑛𝑆−_𝑛𝜔₂𝑡_𝑛2−𝜆_𝑛2sin𝜆_𝑛2𝐿_𝑛𝑘_𝑛𝐺_𝑛𝑆_𝑛𝜆_𝑛2−𝑡_𝑛2cos𝜆_𝑛2𝐿_𝑛−𝐸_𝑛𝑀_𝑛co000s𝜆_𝑛2𝐿_𝑛−𝜔²00𝑀_𝑛si0n𝜆_𝑛2𝐿_𝑛𝑆_𝑛𝑐𝜔_𝑛sin𝑐_𝑛𝜔𝐿_𝑛𝐸^𝑛𝑆_𝑛−𝜔₂𝑐𝜔_𝑛cos𝑐_𝑛𝜔𝐿_𝑛𝑀_𝑛cos𝜔𝑐_𝑛𝐿_𝑛−𝜔²𝑀_𝑛sin𝜔𝑐_𝑛𝐿_𝑛⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠.(A.4)

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

𝐌_{𝐎𝑖,𝑖−1}(𝜔)=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝^{−𝛼𝑖−1𝛼sin𝑖−1𝐿𝑖−1𝛼𝑖−1𝛼cos𝑖−1𝐿𝑖−1𝛼𝑖−1𝛼sinh𝑖−1𝐿𝑖−1𝛼cos𝑖−1𝐿𝑖−1𝜃sin𝑖𝛼sin𝑖−1𝐿𝑖−1𝜃sin𝑖𝛼cosh𝑖−1𝐿𝑖−1𝜃sin𝑖𝛼cos𝑖−1𝐿𝑖−1𝜃cos𝑖𝛼sin𝑖−1𝐿𝑖−1𝜃cos𝑖𝛼cosh𝑖−1𝐿𝑖−1𝜃cos𝑖−𝐸}_𝑖−1^𝐼_𝑖−1^{𝛼2𝑖−1𝛼cos}_𝑖−1^𝐿_{𝑖−1+𝜔2𝐽𝑖−1𝛼𝑖−1𝛼sin𝑖−1𝐿𝑖−1−𝐸𝑖−1}^𝐼_𝑖−1^{𝛼2𝑖−1𝛼sin}_𝑖−1^𝐿_{𝑖−1−𝜔2𝐽𝑖−1𝛼𝑖−1𝛼cos𝑖−1𝐿𝑖−1𝐸𝑖−1}^𝐼_𝑖−1^{𝛼2𝑖−1𝛼cosh}_𝑖−1^𝐿_{𝑖−1−𝜔2𝐽𝑖−1𝛼𝑖−1𝛼sinh𝑖−1𝐿𝑖−1𝐸𝑖−1𝐼𝑖−1𝛼3𝑖−1𝛼sin𝑖−1𝐿𝑖−1−𝜔2𝑀𝑖−1𝛼𝑖−1𝛼cos𝑖−1𝐿𝑖−1𝐸𝑖−1𝐼𝑖−1𝛼3𝑖−1𝛼cos𝑖−1𝐿𝑖−1−𝜔2𝑀𝑖−1𝛼𝑖−1𝛼sin𝑖−1𝐿𝑖−1−𝐸𝑖−1𝐼𝑖−1𝛼3𝑖−1𝛼sinh𝑖−1𝐿𝑖−1−𝜔2𝑀𝑖−1𝛼𝑖−1𝛼cosh𝑖−1𝐿𝑖−1𝛼000𝑖−1𝛼cosh𝑖−1𝐿𝑖−1𝛼00sinh𝑖−1𝐿𝑖−1𝜃sin𝑖𝜔cos𝑐𝑖−1𝐿𝑖−1𝜃cosi𝜔sin𝑐𝑖−1𝐿𝑖−1𝜃cosi𝛼sinh𝑖−1𝐿𝑖−1𝜃sin𝑖𝜔−cos𝑐𝑖−1𝐿𝑖−1𝜃sini𝜔−sin𝑐𝑖−1𝐿𝑖−1𝜃sini𝐸𝑖−1𝐼𝑖−1𝛼2𝑖−1𝛼sinh𝑖−1𝐿𝑖−1−𝜔2𝐽𝑖−1𝛼𝑖−1𝛼cosh𝑖−1𝐿𝑖−100−𝐸𝑖−1𝐼𝑖−1𝛼3𝑖−1𝛼cosh𝑖−1𝐿𝑖−1−𝜔2𝑀𝑖−1𝛼𝑖−1𝛼sinh𝑖−1𝐿𝑖−1000−𝐸𝑖−1𝑆𝑖−1𝜔𝑐𝑖−1𝜔sin𝑐𝑖−1𝐿𝑖−1−𝜔2𝑀𝑖−1𝜔cos𝑐𝑖−1𝐿𝑖−1𝐸𝑖−1𝑆𝑖−1𝜔𝑐𝑖−1𝜔cos𝑐𝑖−1𝐿𝑖−1−𝜔2𝑀𝑖−1𝜔sin𝑐𝑖−1𝐿𝑖−1}⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠,(B.1)𝐌_{𝐎𝑖,𝑖}(𝜔)=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝0−𝛼^𝑖𝐸−0100−001−00𝛼_𝑖−001000_𝑖𝐼0_𝑖𝛼^2𝑖−𝐸_𝑖0−𝐸_𝑖𝐼_𝑖𝛼^2𝑖000𝐼_𝑖𝛼^3𝑖cos𝜃_𝑖0𝐸_𝑖𝐼_𝑖𝛼^3𝑖cos𝜃_𝑖0−𝐸_𝑖𝑆_𝑖0𝐸_𝑖𝜔𝑐_𝑖sin𝜃^𝑖𝐼_𝑖𝛼^3𝑖sin𝜃_𝑖0−𝐸_𝑖𝐼_𝑖𝛼^3𝑖sin𝜃_𝑖0−𝐸_𝑖𝑆_𝑖𝜔𝑐_𝑖cos𝜃^𝑖⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠,(B.2)𝐌_𝐎𝟏(𝜔)=⎛⎜⎜⎜⎝010100010100001000⎞⎟⎟⎟⎠,(B.3)𝐌_𝐎𝐄𝐄(𝜔)=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝−𝐸^𝑛+𝜔₂𝐼_𝑛𝛼^2𝑛cos𝛼_𝑛𝐿_𝑛𝐽_𝑛𝛼_𝑛sin𝛼_𝑛𝐿_𝑛−𝐸^𝑛−𝜔₂𝐼_𝑛𝛼^2𝑛sin𝛼_𝑛𝐿_𝑛𝐽_𝑛𝛼_𝑛cos𝛼_𝑛𝐿_𝑛𝐸^𝑛−𝜔₂𝐼_𝑛𝛼^2𝑛cosh𝛼_𝑛𝐿_𝑛−𝐸_𝑛𝐽_𝑛𝛼_𝑛sinh𝛼_𝑛𝐿_𝑛−𝜔₂𝐼_𝑛𝛼^3𝑛sin𝛼_𝑛𝐿_𝑛𝑀_𝑛cos𝛼_𝑛𝐿_𝑛𝐸^𝑛−𝜔𝐼_𝑛2𝛼^3𝑛cos𝛼_𝑛𝐿_𝑛𝑀_𝑛sin𝛼_𝑛𝐿_𝑛−𝐸^𝑛−𝜔₂𝐼_𝑛𝛼^3𝑛sinh𝛼_𝑛𝐿_𝑛𝐸00𝑀_𝑛co0_𝑛sh𝛼_𝑛𝐿_𝑛−𝜔₂𝐼_𝑛𝛼^2𝑛sinh𝛼_𝑛𝐿_𝑛−𝐸_𝑛𝐽_𝑛𝛼_𝑛cosh𝛼_𝑛𝐿_𝑛00−𝜔₂𝐼_𝑛𝛼^3𝑛cosh𝛼_𝑛𝐿_𝑛𝑀_𝑛sin0hℎ𝛼_𝑛𝐿_𝑛−𝐸_𝑛00𝑆_𝑛−𝜔₂𝑐𝜔_𝑛sin𝑐_𝑛𝜔𝐿_𝑛𝑀_𝑛cos𝜔𝑐_𝑛𝐿_𝑛𝐸^𝑛𝑆_𝑛−𝜔₂𝑐𝜔_𝑛cos𝑐_𝑛𝜔𝐿_𝑛𝑀_𝑛sin𝜔𝑐_𝑛𝐿_𝑛⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠,(B.4)𝛼_i=⁴𝜔²𝜌_iS_iE_iI_i.(B.5)

Funding

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