Research Article  Open Access
Michał Hać, "On Equations of Motion of Elastic Linkages by FEM", Mathematical Problems in Engineering, vol. 2013, Article ID 648706, 10 pages, 2013. https://doi.org/10.1155/2013/648706
On Equations of Motion of Elastic Linkages by FEM
Abstract
Discussion on equations of motion of planar flexible mechanisms is presented in this paper. The finite element method (FEM) is used for obtaining vibrational analysis of links. In derivation of dynamic equations it is commonly assumed that the shape function of elastic motion can represent rigidbody motion. In this paper, in contrast to this assumption, a model of the shape function specifically dedicated to the rigidbody motion is presented, and its influence on elastic motion is included in equations of motion; the inertia matrix related to the rigidbody acceleration vector depends on both shape functions of the elastic and rigid elements. The numerical calculations are conducted in order to determine the influence of the assumed shape function for rigidbody motion on the vibration of links in the case of closedloop and openloop mechanisms. The results of numerical simulation show that for transient analysis and for some specific conditions (e.g., starting range, openloop mechanisms) the influence of assumed shape functions on vibration response can be quite significant.
1. Introduction
In recent years considerable attention has been given to the analysis of flexible mechanisms. The need for taking flexibility of linkages into account is due to much higher speed of operation and the greater restrictions on weight and power requirements of mechanism members. Consequently, for highspeed operation, both rigid body and elastic effects should be included in the mechanism design process in order to reduce dynamic reactions and allow the linkage to perform its prescribed kinematic functions.
Many papers published in recent years involved applications of specialpurpose finite element methods [1]. In these studies, the response of a mechanism was typically obtained by using socalled “superposition theory”. In this method, rigidbody motion is first determined, and then the unknown elastic displacements are solved. In the derivation of equations of motion an assumption is made that the shape function for rigidbody motion is the same as for elastic motion [2–6]. As a result, the system inertia matrix for elastic displacement vector appears also together with the rigidbody acceleration vector on the righthand side of equations of motions.
Relying on the previously described method, however, does not yield accurate results when highspeed systems are considered, since it does not provide for the mutual dynamic coupling of both the rigid and elastic motions. Analytical procedure developed by Song and Haug [7] introduced the “one pass” method which models both the large motions and small elastic displacements of the links simultaneously. However, the interconnections of the bodies are described by a large set of constraint equations formulated for each type of joints. This procedure increases the dimension of the problem considerably and yields to the set of equations composed of both differential and algebraic equations, which are not easy to solve. Later works by Nagarajan and Turcic [8, 9] and others [10, 11] allowed for the rigidbody motion and the elastic motion to influence each other, and both rigidbody degrees of freedom and elastic degrees of freedom are considered as generalized coordinates in the derivation.
The later works were directed towards developing higher order element model [12] or to include the influence of geometric nonlinearities in highspeed applications [13–15]. The paper of Shaker and Ghosal [14] dealt with nonlinear modeling of planar manipulators; equations of motion were obtained using Lagrangian formulation and were nondimensionalized using two characteristic velocities. These values allowed them to predict whether the nonlinear analysis was needed or the linear case produced satisfactory results. Korayem et al. [15] developed the Lagrange approach to flexible mobile manipulators for determining maximum allowable dynamic load.
The most often method used for obtaining equations of motion for flexiblebody systems is the Lagrange’s approach [16]. However, in comparison with Lagrange's equations of motion the use of the GibbsAppel equations [17, 18] simplifies the derivation of governing dynamic equations. In this second method the second derivative of the generalized vector is firstly calculated and then the appropriate second partial derivative of the Gibbs function is obtained. In Lagrange's equations of motion the sequence is quite opposite, and when the general point displacement is presented in a vector form, Lagrange's equations are not particularly well suited.
In the present paper, a form of equations of motion used in the superposition theory of elastic linkages is discussed. The model of a shape function for the planar rigidbody motion proposed in [19, 20] is used. In this way, in contrast to many other studies [2–6], the shape function for elastic motion is not used to describe an arbitrary large rigidbody translation. Moreover, as proven in the Appendix, the shape function for planar beam elements is only an approximation of displacement distribution along rigid finite elements and for exact description of rigid motion the shape function proposed in this paper should be used. The use of different shape functions for elastic and rigid motions implies that the inertia matrix, standing by rigidbody acceleration vector in the equations of motion of flexible mechanisms, depends on both shape functions of elastic and rigid elements. The influence of assumed shape functions on vibrational response of mechanism links is discussed in the paper.
2. Equations of Motion
The coupling between the nonlinear rigidbody motions and the linear small elastic deformation stands as the main problem in the solution of the dynamics of flexible mechanisms. In the superposition theory the elastic displacements are separated from the rigidbody displacements (which are previously calculated) and solved as the unknowns of the system. In this way the rigidbody motion is not influenced by the elastic motion, but the vibration of links appears due to centrifugal forces generated in rigidbody motion of a mechanism.
Figure 1 presents a general planar beam element in three frames of reference. Frame is fixed to the ground and serves as the global coordinate system. The element oriented frame is a rotating reference frame and its axis is parallel to the undeformed center line of the element. The origin of the frame is located at a node point of a finite element and it remains parallel to the system. The frames and are updated continually as the element moves.
Let us assume that during the motion of the mechanism a finite element has changed its position determined by endpoints 1 and 2 to the points 1′ and 2′ (dotted) due to the rigidbody motion and to the points and due to the elastic deformation. The components of nodal displacement vectors 11′ and 22′ can be expressed in the relative coordinate system by the following vector: or in the local coordinate system by where are displacements of nodes 1 and 2 in and direction, respectively, are displacements of nodes 1 and 2 in and direction, respectively.
The displacement vector of a general point belonging to the finite element can be expressed in terms of the displacements of the element’s endpoints using the formula where is the shape function for the rigid body motion which can be expressed as [6] where , and is the following transformation matrix: where is the angle between and the global coordinate systems. Since the longitudinal displacements of nodes in the case of rigid motion are equal, the shape function for the rigidbody motion (4) can be presented in other forms as shown in [22].
The position of the general point in the global coordinate system (see Figure 1) is denoted by and can be written as where is the displacement connected with the rigidbody motion and is an elastic deformation vector in the local coordinate system . Vector can be defined as a sum of the initial position vector and displacement vector of the general point due to rigidbody motion By differentiating (6) with respect to time the velocity and acceleration of the general point can be obtained where is the absolute angular velocity of the frame and is Booleantype operator matrix defined as In derivation of (8) the constancy of the initial position vector (see Figure 1) has been taken into account which implies the following identities: , .
Now let us consider the elastic motion which describes the vibrations of links. The components of nodal displacement vectors and due to elastic deformation can be expressed in the and systems by the vectors and , respectively, where: are displacements of nodes 1 and 2 in and direction, respectively, and are angular deformations of nodes 1 and 2, where are displacements of nodes 1 and 2 in and direction, respectively.
The vector is transformed into by where is the transformation matrix Introducing the matrix notation essential in the finite element formulation, the displacement vector can be expressed in terms of the displacements of the element's endpoints: The shape function for elastic displacement is as follows:
where .
Equations of motion are obtained by making use of GibbsAppel equations [11]: where represents the generalized nodal degree of freedom, is the potential energy, are the generalized forces acting on the element, and is the Gibbs function or “energy of acceleration”. The energy of dissipation is not involved in formulated equations of motion; in case of structural damping the damping matrices can be introduced in system equations of motion as a linear combination of system inertia and stiffness matrices.
The Gibbs function is expressed as follows: where and are the shape function for the rigid body motion (4) and elastic motion (15), respectively.
The strain energy of the element is associated with vector and is given by where the element stiffness matrix is in the same form as for beam elements.
Substituting expressions on Gibbs function and strain energy of an element into GibbsAppel equations (16), the equations of motion for one element are obtained in the following form: where is the element inertia matrix, is the centrifugal matrix, is the gyroscopic matrix, and is the coefficient matrix defined by the following notation: and are generalized forces acting on finite element.
At the level of element equations of motion there exist forces from adjacent links which at the system size level are internal forces and occur in positive and negative pairs which cancel. Taking into account the boundary conditions and introducing damping matrix , the equations of motion for the system are stated as [20] where coefficient matrices are global matrices obtained from appropriate element matrices, represents generalized forces, and , , and represent acceleration, velocity, and displacement vectors (in nodal points). Vectors and are assembled forms of element vectors and with regard to boundary conditions.
By neglecting the rigid body and elastic motion coupling terms, the equations of motion (21) are obtained in a simplified form as
In most studies on the flexible multibody formulation [2–6] it is assumed that the element shape function adopted for elastic motion can be used to describe an arbitrary large rigidbody translation. Thus the shape function for elastic displacement (15) can be used to determine rigidbody displacement, that is, where vector is similar to the one given by (11), but the displacements are due to rigidbody motion.
By comparison to (3) it is obviously seen that the rigidbody shape function appearing in the Gibbs function (17) is now changed by the elastic shape function . The coefficient matrix denoted previously by (see (20)) is now equal to , that is, The assumption formulated in (23) leads to equations of motion in which the inertia matrix related to the rigidbody acceleration vector is the same as the system inertia matrix. In simplified form the equations of motion in global coordinate system are as follows [3]: By comparing (22) and (25) it can be seen that, assuming the shape function for rigidbody motion (4), the inertia matrices related to the rigidbody acceleration vector at righthand side of equations of motion differ considerably in the two cases considered.
The solution of matrix (22) and (25) is obtained by using Newmark method for integrating the equations. In the numerical examples, the gravity and damping effects are not taken into account.
3. Illustrative Examples
In order to determine the influence of the assumed rigidbody motion shape function on vibration of links, two numerical examples are considered: the first one is an example of a closedchain mechanism, a fourbar linkage, and the second one represents an openchain mechanism of twolink planar manipulator. It is also interesting to compare the proposed method with the commonly used assumption of applying elementtype (in our case beam finite element) shape function to rigidbody motion. Thus, the calculations are conducted for two cases.
Case 1. By solving equations of motion (22) which have been formulated for the shape function for the rigidbody motion given by (4).
Case 2. By solving commonly used equations of motion of type (25).
3.1. FourBar Linkage
All three moving members of the linkage are assumed to be deformable. A total of six elements were employed—each link is represented by two finite elements. The corresponding placement of nodes is shown in Figure 2.
There are nodes (numbers 3 and 5) at the revolute connection between the crank and the coupler, and between the coupler and the follower, which are modelled by assuming two independent rotational d.o.f. (left and right) at each of these nodes. Boundary conditions of zero translational displacement were imposed at nodes 1 and 7.
During numerical simulations the equations of motion given by (22) and (25) were solved. In both cases the vector of the unknown functions (nodal displacement vector ) consists of eighteen elements: displacements and in the and directions, respectively, of moving nodes 2, 3, 4, 5, and 6 and the nodal deformation angles The system matrices for elastic vibration , , and are of dimensions . The order of the nodal displacement vector for the rigidbody motion differs depending on the case considered. In (22) it consists of ten elements, that is, the displacements of moving nodes in the and direction, respectively, due to the rigidbody motion (i.e., all displacements are with the index 0): In this way the matrix in (22) is of dimension due to different size of vectors and . In the case of (25) the rigidbody displacement vector consists of 18 elements since it stands near global inertia matrix formulated for flexible elements. Thus the vector contains elements appearing in the beam finite elements: in addition to displacements there are angular deformations which in the case of rigidbody motion are set to be equal to 0. The rigidbody displacement vector in the case of (25) is as follows: It is seen that the dimension of the rigidbody vector in the proposed method is lower (compare (27) and (28)) and assuming the shape function for rigidbody motion in the form of elastic motion leads to artificial substitution of these angle deformations for rigid motion as equal to 0.
The geometrical data taken for numerical simulation is given in Table 1 [3].

The input torque applied to the crank of the mechanism is assumed to be as follows: where is the crank angle. The data given for the numerical simulation are Nm, Nm. The initial conditions (crank angle and crank angular velocity) are rad, rad/sec. Simulation time is from 0 to 0.25 sec and assures the crank to perform a full rotation.
The results of computer simulation show that, in the case of considered fourbar linkage, the differences in equations of motion do not significantly influence the results. The midspan bending strains for the coupler (Figure 3) and for the follower (Figure 4) are very similar for the two cases considered, and the relative difference in results does not exceed 3%. Thus, it would hardly be visible in the figures, and additional charts are provided in order to show differences in midspan strains: in Figure 5 for the coupler and in Figure 6 for the follower.
The difference in equations of motion (22) and (25) representing Cases 1 and 2, respectively, is only in the inertia force vector on the righthand sides of these equations. For (22) this vector is equal to , and for (25) . It was of interest to investigate how close to each other are these forces. In Figure 7, the inertia force as function of the crank angle is presented; in both cases the results are very close and the differences are not visible in the figure. Thus in the next Figure 8, the difference between these forces (i.e., ) is presented; its value is approximately 1% of the value of component forces.
3.2. TwoLink Planar Manipulator
The second example concerns the openloop mechanism of a twolink planar manipulator (Figure 9). The mechanical basic characteristics of this manipulator are the same as in [21] and are given in Table 2. In addition, the following initial properties are the mass at the joint 3 ( = 0.5 kg), the moment of inertia at the hub kg m^{2}, and the nominal payload kg are taken into account.

The vibration of links arises due to the motion of the manipulator excited by external torques applied to the manipulator links at the joints. It was assumed that the values of torques are as follows: . The vibration analysis is conducted by solving equations of motion (22) and (25) which correspond to the two cases considered. The initial conditions are rad, rad, zero angular velocities; the time of simulation is equal to 0.25 sec.
Taking into account the elastic vibration of links the nodal displacement vector in (22) or (25) consists of 12 elements: displacements and in the and directions, respectively of moving nodes 2, 3, 4, and 5, and nodal deformation angles : From the reason described in Section 3.1 the nodal displacement vector for the rigidbody motion in the case of (22) is equal to and in the case of (25) of the form The system matrices , , and in (22) and (25) are of dimensions , and the matrix in (22) is of dimensions .
The horizontal and vertical displacements of manipulator’s tip are presented in Figures 10 and 11, respectively. Now the difference in results for the two cases considered is considerable and can be clearly seen. In particular the displacements in direction for starting range differ significantly. It should be noted, however, that the calculations have been performed for transient motion with sudden torques applied to both links of the manipulator.
In order to analyze the influence of payload on displacements of manipulator’s tip, the numerical simulation with no payload was conducted. The results in vertical direction are presented in Figure 12. It can be seen that for the given data, due to lower inertial forces in the no payload case, the increase of vibration response is observed.
The previous results have been obtained for linear model of planar beam elements which assumes that the transverse displacements are independent of the axial displacements or forces. If, in derivation of equations of motion, the nonlinear straindisplacement relationship is taken into account [14, 20]; the additional geometric stiffness matrix appears in dynamic equations, and it should be added to the element stiffness matrix . For a beam element this matrix can be presented in the following form [23]: where is equal to the axial force in the element.
The additional numerical simulations have been performed in order to evaluate the influence of nonlinear effects on dynamic response of the considered twolink manipulator. The results of midspan bending strains for link I and II are presented in Figures 13 and 14, respectively, for linear (solid line) and nonlinear (dashed line) models.
The results show that in the case considered the geometric stiffness matrix has very small influence on the solution of equations of motion. However, in the case of manipulators with the low actuating torque frequency [14] the difference in results can be significant.
4. Conclusions
A finite element model for analysing the kinetoelastodynamic transient characteristics of a flexible mechanism is developed in the present paper. In multibody dynamics, there is a direct relationship between the selection of mode shapes in assumed coordinate system and the vibration response of constrained deformable bodies. By assuming the proposed shape function for rigidbody motion, the vibration equations of motion are derived, and the comparative analysis with existing commonly used approach is presented. The difference in the equations of motion between the former and proposed approaches is in the coupling between elastic and rigid body modes; namely, the inertia force vector corresponding to the rigidbody acceleration is different for both approaches.
The method presented here is illustrated for the cases of planar open and closedloop linkages. The results of numerical simulation show that in the case of considered fourbar linkage the differences in equations of motion do not considerably influence the results. An opposite observation is made in the case of an openloop mechanism—twolink planar manipulator—the differences in results between two discussed approaches are significant. However, in this case the severe excitation conditions were assumed; the transient response was examined after stepped torques applied to the links. Apart from the starting range, the solutions are very close to each other, which leads to the conclusion that for stable conditions of work (e.g., mechanisms circulating at constant speed, steadystate response is analysed) the two approaches presented in this paper should give quite converged results.
Appendix
In this appendix the derivation of displacements of a finite element in the rigidbody motion for two different shape functions, beam finite element (15), and rigid element (4) is presented.
Let us consider the rigid finite element of length in the local coordinate frame shown in Figure 15.
The components of the nodal displacement vector are expressed in the local coordinate system by (2). The displacement vector of a general point belonging to the finite element and determined by coordinate (see Figure 15) can be expressed in terms of the displacements of the element's endpoints using the formula (3) For rigid elements the displacements in longitudinal direction of all points belonging to the element are equal; that is, the relationship holds for any . In the transverse direction the displacement of any point is a linear combination of transverse displacements of element’s endpoints. Thus, for example, displacements of a general point determined by coordinates and are as follows (note that nondimensional coordinate appearing in the shape function formulation is equal to , so the calculations are provided for and ): If we assume that the shape function for elastic displacement (15) can be used to determine rigidbody displacements, the formula of (23) is used and the general point displacements are formulated as
It is easy to show that, for the given values of , the following displacements are obtained: It is seen that while the formulation for rigidbody displacements for is correct, the result for transverse displacement of the general point for is incorrect.
Using the shape function for the rigidbody motion (4) the general point displacements are formulated as (A.1) and the longitudinal and transverse displacements for the given values of are It is seen that for both values of the obtained results are now correct and are in fact the same as expected in the rigidbody motion and shown in (A.2). The only difference is in the formulation of but in view of the longitudinal deformation in the rigid motion being constant (i.e., ), the obtained formulas give the same results.
From this short derivation one can conclude that the assumption that the shape function for beam finite element can convey also rigidbody motion is only an approximation and that the shape function for rigidbody motion (4) should be used for exact description.
References
 S. K. Dwivedy and P. Eberhard, “Dynamic analysis of flexible manipulators, a literature review,” Mechanism and Machine Theory, vol. 41, no. 7, pp. 749–777, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 A. Midha, A. G. Erdman, and D. A. Frohrib, “Finite element approach to mathematical modeling of highspeed elastic linkages,” Mechanism and Machine Theory, vol. 13, no. 6, pp. 603–618, 1978. View at: Publisher Site  Google Scholar
 D. A. Turcic and A. Midha, “Generalized equations of motion for the dynamic analysis of elastic mechanism systems,” AMSE Journal of Dynamic Systems, Measurement and Control, vol. 106, no. 4, pp. 243–248, 1984. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 Z. Yang and J. P. Sadler, “On issues of elasticrigid coupling in finite element modeling of highspeed machines,” Mechanism and Machine Theory, vol. 35, no. 1, pp. 71–82, 2000. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 P. Fanghella, C. Galletti, and G. Torre, “An explicit independentcoordinate formulation for the equations of motion of flexible multibody systems,” Mechanism and Machine Theory, vol. 38, no. 5, pp. 417–437, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 W. Hou and X. Zhang, “Dynamic analysis of flexible linkage mechanisms under uniform temperature change,” Journal of Sound and Vibration, vol. 319, no. 12, pp. 570–592, 2009. View at: Publisher Site  Google Scholar
 J. O. Song and E. J. Haug, “Dynamic analysis of planar flexible mechanisms,” Computer Methods in Applied Mechanics and Engineering, vol. 24, no. 3, pp. 359–381, 1980. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 S. Nagarajan and D. A. Turcic, “Lagrangian formulation of the equations of motion for elastic mechanisms with mutual dependence between rigid body and elastic motions—part I: element level equations,” ASME Journal of Dynamic Systems, Measurement and Control, vol. 112, no. 2, pp. 203–214, 1990. View at: Publisher Site  Google Scholar
 S. Nagarajan and D. A. Turcic, “Lagrangian formulation of the equations of motion for elastic mechanisms with mutual dependence between rigid body and elastic motions—part II: system equations,” ASME Journal of Dynamic Systems, Measurement and Control, vol. 112, no. 2, pp. 215–224, 1990. View at: Publisher Site  Google Scholar
 M. Hać, “Dynamics of flexible mechanisms with mutual dependence between rigid body motion and longitudinal deformation of links,” Mechanism and Machine Theory, vol. 30, no. 6, pp. 837–847, 1995. View at: Publisher Site  Google Scholar
 X. Zhang, H. Liu, and Y. Shen, “Finite dynamic element analysis for highspeed flexible linkage mechanisms,” Computers and Structures, vol. 60, no. 5, pp. 787–796, 1996. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. B. Gunda, R. K. Gupta, and R. Ganguli, “Hybrid stiffstringpolynomial basis functions for vibration analysis of high speed rotating beams,” Computers and Structures, vol. 87, no. 34, pp. 254–265, 2009. View at: Publisher Site  Google Scholar
 R. Avilés, A. Hernández, E. Amezua, and O. Altuzarra, “Kinematic analysis of linkages based in finite elements and the geometric stiffness matrix,” Mechanism and Machine Theory, vol. 43, no. 8, pp. 964–983, 2008. View at: Publisher Site  Google Scholar
 M. C. Shaker and A. Ghosal, “Nonlinear modeling of flexible manipulators using nondimensional variables,” Journal of Computational and Nonlinear Dynamics, vol. 1, no. 2, pp. 123–134, 2006. View at: Publisher Site  Google Scholar
 M. H. Korayem, A. Heidari, and A. Nikoobin, “Maximum allowable dynamic load of flexible mobile manipulators using finite element approach,” International Journal of Advanced Manufacturing Technology, vol. 36, no. 56, pp. 606–617, 2008. View at: Publisher Site  Google Scholar
 P. Fanghella, C. Galletti, and G. Torre, “An explicit independentcoordinate formulation for the equations of motion of flexible multibody systems,” Mechanism and Machine Theory, vol. 38, no. 5, pp. 417–437, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 L. A. Pars, A Treatise on Analytical Dynamics, Heinemann, London, UK, 1965.
 M. Hać and K. From, “Design of retraction mechanism of aircraft landing gear,” Mechanics and Mechanical Engineering, vol. 12, no. 4, pp. 357–373, 2008. View at: Google Scholar
 M. Hać and J. Osiński, “Finite element formulation of rigid body motion in dynamic analysis of mechanisms,” Computers and Structures, vol. 57, no. 2, pp. 213–217, 1995. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Hać, “Dynamic analysis of flexible mechanisms by finite element method,” Machine Dynamics Problems, vol. 14, pp. 7–91, 1996. View at: Google Scholar
 J. N. Bricout, J. C. Debus, and P. Micheau, “A finite element model for the dynamics of flexible manipulators,” Mechanism and Machine Theory, vol. 25, no. 1, pp. 119–128, 1990. View at: Publisher Site  Google Scholar
 M. Hać, “On equations of motion of flexible mechanisms,” Mechanics Research Communication, vol. 20, pp. 399–403, 1993. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. S. Przemieniecki, Theory of Matrix Structural Analysis, McGrawHill, London, UK, 1968.
Copyright
Copyright © 2013 Michał Hać. 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.