# On Equations of Motion of Elastic Linkages by FEM

Authors: Michał Hać
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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 rigid-body motion. In this paper, in contrast to this assumption, a model of the shape function specifically dedicated to the rigid-body motion is presented, and its influence on elastic motion is included in equations of motion; the inertia matrix related to the rigid-body 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 rigid-body motion on the vibration of links in the case of closed-loop and open-loop mechanisms. The results of numerical simulation show that for transient analysis and for some specific conditions (e.g., starting range, open-loop mechanisms) the influence of assumed shape functions on vibration response can be quite significant.


---

## Body

## 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 high-speed 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 special-purpose finite element methods [1]. In these studies, the response of a mechanism was typically obtained by using so-called “superposition theory”. In this method, rigid-body 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 rigid-body 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 rigid-body acceleration vector on the right-hand side of equations of motions.Relying on the previously described method, however, does not yield accurate results when high-speed 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 rigid-body motion and the elastic motion to influence each other, and both rigid-body 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 high-speed 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 flexible-body systems is the Lagrange’s approach [16]. However, in comparison with Lagrange's equations of motion the use of the Gibbs-Appel 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 rigid-body 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 rigid-body 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 rigid-body 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 rigid-body 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 rigid-body displacements (which are previously calculated) and solved as the unknowns of the system. In this way the rigid-body motion is not influenced by the elastic motion, but the vibration of links appears due to centrifugal forces generated in rigid-body motion of a mechanism.Figure1 presents a general planar beam element in three frames of reference. Frame XY 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 xy frame is located at a node point of a finite element and it remains parallel to the XY system. The frames xy and ξη are updated continually as the element moves.Figure 1
Displacements of a finite element in global and local coordinates.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 rigid-body motion and to the points1′′ and 2′′ due to the elastic deformation. The components of nodal displacement vectors 11′ and 22′ can be expressed in the relative coordinate system xy by the following vector:
(1){s0}T=[u01,w01,u02,w02]
or in the local coordinate system  ξη  by
(2){δ0}T=[p01,v01,p02,v02],
where u01,w02,u01,w02 are displacements of nodes 1 and 2 in  x and y direction, respectively,  p01,p02,v01,v02 are displacements of nodes 1 and 2 in  ξ and  η direction, respectively.The displacement vectorr0 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
(3)r0≡{r0}=[r0x,r0y]T=[Ne0]{δ0}=[Ne0][T0]{s0},
where [Ne0] is the shape function for the rigid body motion which can be expressed as [6]
(4)[Ne0]=[0.500.5001-ξ0ξ],
where 0≤ξ≤1, and [T0] is the following transformation matrix:
(5)[T0]=[cosαsinα00-sinαcosα0000cosαsinα00-sinαcosα],
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 rigid-body motion (4) can be presented in other forms as shown in [22].The position of the general point in the global coordinate system (see Figure1) is denoted by r and can be written as
(6)r=R0+q,
where R0 is the displacement connected with the rigid-body motion and q is an elastic deformation vector in the local coordinate system  ξη. Vector R0 can be defined as a sum of the initial position vector rc and displacement vector r0 of the general point due to rigid-body motion
(7)R0=rc+r0.
By differentiating (6) with respect to time the velocity and acceleration of the general point can be obtained
(8)r˙=R˙0+ω-×q+q˙=r˙0+ω[Ω]q+q˙,r¨=r¨0+2ω[Ω]q˙-ω2q+q¨+ω˙[Ω]q,
where ω-=col[0,0,ω] is the absolute angular velocity of the ξη  frame and [Ω] is Boolean-type operator matrix defined as
(9)[Ω]=[0-110].
In derivation of (8) the constancy of the initial position vector rc (see Figure 1) has been taken into account which implies the following identities: R˙0=r˙0,  R¨0=r¨0.Now let us consider the elastic motion which describes the vibrations of links. The components of nodal displacement vectors1′1′′ and 2′2′′ due to elastic deformation can be expressed in the xy and ξη systems by the vectors {s} and {δ}, respectively,
(10){s}T=[u1,w1,Θ1,u2,w2,Θ2],
where:  u1,w1,u2,w2 are displacements of nodes 1 and 2 in x and y direction, respectively, and  Θ1,Θ2  are angular deformations of nodes 1 and 2,
(11){δ}T=[p1,v1,Θ1,p2,v2,Θ2],
where  p1,p2,v1,v2 are displacements of nodes 1 and 2 in ξ and  η direction, respectively.The vector{δ} is transformed into {s} by
(12){δ}=[T]{s},
where [T] is the transformation matrix
(13)[T]=[cosαsinα0000-sinαcosα0000001000000cosαsinα0000-sinαcosα0000001].
Introducing the matrix notation essential in the finite element formulation, the displacement vector q can be expressed in terms of the displacements of the element's endpoints:
(14)q≡{q}=[Ne]{δ}=[Ne][T]{s}.
The shape function [Ne] for elastic displacement is as follows:(15)Ne=[1-ξ00ξ0001-3ξ2+2ξ3L(ξ-2ξ2+ξ3)03ξ2-2ξ3L(ξ3-ξ2)],where0≤ξ≤1.Equations of motion are obtained by making use of Gibbs-Appel equations [11]:
(16)∂Ge∂{δ¨}+∂Ue∂{δ}={Q},
where {δ} represents the generalized nodal degree of freedom, Ue is the potential energy, {Q} are the generalized forces acting on the element, and Ge 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:(17)Ge=12∫0Lr¨Tr¨dm=∫0L(({δ¨0}T[Ne0]T+2ω{δ˙}T[Ne]T[Ω]T-ω2{δ}T[Ne]T+{δ¨}T[Ne]T+ω˙{δ}T[Ne]T[Ω]T)×([Ne0]{δ¨0}+2ω[Ω][Ne]{δ˙}-ω2[Ne]{δ}+[Ne]{δ¨}+ω˙[Ω][Ne]{δ}{δ¨0}){δ¨0}T)dm,
where [Ne0] and [Ne] are the shape function for the rigid body motion (4) and elastic motion (15), respectively.The strain energyUe of the element is associated with vector {δ} and is given by
(18)Ue=12{δ}T[Ke]{δ},
where the element stiffness matrix [Ke] is in the same form as for beam elements.Substituting expressions on Gibbs function and strain energy of an element into Gibbs-Appel equations (16), the equations of motion for one element are obtained in the following form:
(19)[Me]{δ¨}+[Cge]{δ˙}+([Ke]+[Kce]){δ}={Fe}-[Me0]{δ¨0},
where [Me] is the element inertia matrix,  [Kce] is the centrifugal matrix, [Cge] is the gyroscopic matrix, and  [Me0′] is the coefficient matrix defined by the following notation:
(20){δ˙}T[Me0′]{δ˙0}=12∫0L({δ˙0}T[Ne0]T[Ne]{δ˙}+{δ˙}T[Ne]T[Ne0]{δ˙0})dm,
and {Fe} 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[C], the equations of motion for the system are stated as [20]
(21)[M]{x¨}+[C]{x˙}+[Cg]{x˙}+([K]+[Kc]){x}={F}-[M0]{x¨0},
where coefficient matrices are global matrices obtained from appropriate element matrices, {F} represents generalized forces, and {x¨}, {x˙}, and {x} represent acceleration, velocity, and displacement vectors (in nodal points). Vectors {x} and {x0} are assembled forms of element vectors {s} and {s0} 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
(22)[M]{x¨}+[C]{x˙}+[K]{x}={F}-[M0]{x¨0}.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 rigid-body translation. Thus the shape function for elastic displacement [Ne] (15) can be used to determine rigid-body displacement, that is,
(23)r0=[Ne]{δ0′},
where vector {δ0′} is similar to the one given by (11), but the displacements are due to rigid-body motion.By comparison to (3) it is obviously seen that the rigid-body shape function [Ne0] appearing in the Gibbs function (17) is now changed by the elastic shape function [Ne]. The coefficient matrix denoted previously by [Me0′] (see (20)) is now equal to [Me], that is,
(24)[Me]=∫0L[Ne]T[Ne]dm.
The assumption formulated in (23) leads to equations of motion in which the inertia matrix related to the rigid-body 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]:
(25)[M]{x¨}+[C]{x˙}+[K]{x}={F}-[M]{x¨0}.
By comparing (22) and (25) it can be seen that, assuming the shape function for rigid-body motion (4), the inertia matrices related to the rigid-body acceleration vector at right-hand side of equations of motion differ considerably in the two cases considered.The solution{x} 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 rigid-body motion shape function on vibration of links, two numerical examples are considered: the first one is an example of a closed-chain mechanism, a four-bar linkage, and the second one represents an open-chain mechanism of two-link planar manipulator. It is also interesting to compare the proposed method with the commonly used assumption of applying element-type (in our case beam finite element) shape function to rigid-body 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 rigid-body motion given by (4).Case 2.
By solving commonly used equations of motion of type (25).


### 3.1. Four-Bar 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 Figure2.Figure 2
Four-bar linkage with nodes.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 {x}) consists of eighteen elements: displacements u and w in the X and Y directions, respectively, of moving nodes 2, 3, 4, 5, and 6 and the nodal deformation angles Θ(26){x}T=[u2,w2,Θ2,u3,w3,Θ31,Θ32,u4,w4,Θ4,u5,w5,Θ51,Θ52,u6,w6,Θ6,Θ7].
The system matrices for elastic vibration [M], [C], and [K] are of dimensions 18×18. The order of the nodal displacement vector for the rigid-body motion {x0} differs depending on the case considered. In (22) it consists of ten elements, that is, the displacements of moving nodes in the X and Y direction, respectively, due to the rigid-body motion (i.e., all displacements are with the index 0):
(27){x0}T=[u02,w02,u03,w03,u04,w04,u05,w05,u06,w06].
In this way the matrix [M0′] in (22) is of dimension 18×10 due to different size of vectors {x} and {x0}. In the case of (25) the rigid-body displacement vector {x0} consists of 18 elements since it stands near global inertia matrix [M] formulated for flexible elements. Thus the vector {x0} contains elements appearing in the beam finite elements: in addition to displacements there are angular deformations which in the case of rigid-body motion are set to be equal to 0. The rigid-body displacement vector {x0} in the case of (25) is as follows:
(28){x0}T=[u02,w02,0,u03,w03,0,0,u04,w04,0,u05,w05,0,0,u06,w06,0,0].
It is seen that the dimension of the rigid-body vector in the proposed method is lower (compare (27) and (28)) and assuming the shape function for rigid-body 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 Table1 [3].Table 1
Four-bar linkage parameters [3].
Parameters
Crank
Coupler
Follower
Length, m
0.10795
0.2794
0.27051
Cross-section, m2
1.07
×
10
-
4
0.406
×
10
-
4
0.406
×
10
-
4
Bending moments of inertia, m4
1.62
×
10
-
10
8.67
×
10
-
12
8.67
×
10
-
12
Distance between ground pivots:   0.254 m.
Lumped mass at node 3:   0.0450 kg
Lumped mass at node 5:   0.0375 kg
Modulus of elasticity:       0.71× 1011 N/m2
Mass density:       2710 kg/m3The input torqueMt applied to the crank of the mechanism is assumed to be as follows:
(29)Mt=Mt0+Mt1sinα,
where α is the crank angle. The data given for the numerical simulation are Mt0=1 Nm, Mt1=0.5 Nm. The initial conditions (crank angle and crank angular velocity) are α=0 rad, ω=0 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 four-bar linkage, the differences in equations of motion do not significantly influence the results. The midspan bending strains for the coupler (Figure3) 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.Figure 3
Midspan bending strains for a coupler with crank angle.Figure 4
Midspan bending strains for a follower with crank angle.Figure 5
Difference between two cases in midspan coupler strains as a function of the crank angle.Figure 6
Difference between two cases in midspan coupler strains as a function of the crank angle.The difference in equations of motion (22) and (25) representing Cases 1 and 2, respectively, is only in the inertia force vector on the right-hand sides of these equations. For (22) this vector is equal to Fc0=[M0]{x¨0}, and for (25) Fc=[M]{x¨0}. 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.,  Fc-Fc0) is presented; its value is approximately 1% of the value of component forces.Figure 7
Inertia force versus the crank angle.Figure 8
Difference in inertia force between two cases versus the crank angle.


### 3.2. Two-Link Planar Manipulator

The second example concerns the open-loop mechanism of a two-link planar manipulator (Figure9). 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 (m3 = 0.5 kg), the moment of inertia at the hub I3=0.04 kg m2, and the nominal payload  m4=0.3 kg are taken into account.Table 2
Two-link manipulator parameters [21].
Parameters
Link I
Link II
Length, m.
L
1
=
1
L
2
=
1
Cross-section area, m2
F
1
=
4
×
10
-
4
F
2
=
4
×
10
-
4
Moment of inertia, m4
I
1
=
16
/
3
×
10
-
8
I
2
=
16
/
3
×
10
-
8
Modulus of elasticity, N/m2 E=2.1×1011
Mass density, kg/m3
ρ=7800Figure 9
A two-link planar manipulator.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:M1=-M2=10Nm. The vibration analysis is conducted by solving equations of motion (22) and (25) which correspond to the two cases considered. The initial conditions are α=0 rad, β=0.5 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{x} in (22) or (25) consists of 12 elements: displacements u  and v in the X and Y directions, respectively of moving nodes 2, 3, 4, and 5, and nodal deformation angles Θ:
(30){x}T=[u2,w2,Θ2,u3,w3,Θ3,u4,w4,Θ4,u5,w5,Θ5].
From the reason described in Section 3.1 the nodal displacement vector for the rigid-body motion in the case of (22) is equal to
(31){x0}T=[u02,w02,u03,w03,u04,w04,u05,w05]
and in the case of (25) of the form
(32){x0}T=[u02,w02,0,u03,w03,0,u04,w04,0,u05,w05,0].
The system matrices [M], [C], and [K] in (22) and (25) are of dimensions 12×12, and the matrix [M0′] in (22) is of dimensions 12×8.The horizontal and vertical displacements of manipulator’s tip are presented in Figures10 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 Y 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.Figure 10
Displacement of manipulator’s tip inX direction for cases: I (solid line) and II (dashed line) versus time.Figure 11
Displacement of manipulator’s tip inY direction for cases: I (solid line) and II (dashed line) versus time.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 Figure12. 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.Figure 12
Displacement of manipulator’s tip (no payload) inY direction for cases: I (solid line) and II (dashed line) versus time.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 strain-displacement relationship is taken into account [14, 20]; the additional geometric stiffness matrix [kG] appears in dynamic equations, and it should be added to the element stiffness matrix [k]. For a beam element this matrix can be presented in the following form [23]:
(33)[kG]=FL[0065Symm.0L102L21500000-65-L100650L10-L2300-L102L215],
where  F 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 two-link manipulator. The results of midspan bending strains for link I and II are presented in Figures13 and 14, respectively, for linear (solid line) and nonlinear (dashed line) models.Figure 13
Midpoint bending strains for link I for linear (solid line) and nonlinear (dashed line) cases versus time.Figure 14
Midpoint bending strains for link II for linear (solid line) and nonlinear (dashed line) cases versus time.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.


## 3.1. Four-Bar 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 Figure2.Figure 2
Four-bar linkage with nodes.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 {x}) consists of eighteen elements: displacements u and w in the X and Y directions, respectively, of moving nodes 2, 3, 4, 5, and 6 and the nodal deformation angles Θ(26){x}T=[u2,w2,Θ2,u3,w3,Θ31,Θ32,u4,w4,Θ4,u5,w5,Θ51,Θ52,u6,w6,Θ6,Θ7].
The system matrices for elastic vibration [M], [C], and [K] are of dimensions 18×18. The order of the nodal displacement vector for the rigid-body motion {x0} differs depending on the case considered. In (22) it consists of ten elements, that is, the displacements of moving nodes in the X and Y direction, respectively, due to the rigid-body motion (i.e., all displacements are with the index 0):
(27){x0}T=[u02,w02,u03,w03,u04,w04,u05,w05,u06,w06].
In this way the matrix [M0′] in (22) is of dimension 18×10 due to different size of vectors {x} and {x0}. In the case of (25) the rigid-body displacement vector {x0} consists of 18 elements since it stands near global inertia matrix [M] formulated for flexible elements. Thus the vector {x0} contains elements appearing in the beam finite elements: in addition to displacements there are angular deformations which in the case of rigid-body motion are set to be equal to 0. The rigid-body displacement vector {x0} in the case of (25) is as follows:
(28){x0}T=[u02,w02,0,u03,w03,0,0,u04,w04,0,u05,w05,0,0,u06,w06,0,0].
It is seen that the dimension of the rigid-body vector in the proposed method is lower (compare (27) and (28)) and assuming the shape function for rigid-body 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 Table1 [3].Table 1
Four-bar linkage parameters [3].
Parameters
Crank
Coupler
Follower
Length, m
0.10795
0.2794
0.27051
Cross-section, m2
1.07
×
10
-
4
0.406
×
10
-
4
0.406
×
10
-
4
Bending moments of inertia, m4
1.62
×
10
-
10
8.67
×
10
-
12
8.67
×
10
-
12
Distance between ground pivots:   0.254 m.
Lumped mass at node 3:   0.0450 kg
Lumped mass at node 5:   0.0375 kg
Modulus of elasticity:       0.71× 1011 N/m2
Mass density:       2710 kg/m3The input torqueMt applied to the crank of the mechanism is assumed to be as follows:
(29)Mt=Mt0+Mt1sinα,
where α is the crank angle. The data given for the numerical simulation are Mt0=1 Nm, Mt1=0.5 Nm. The initial conditions (crank angle and crank angular velocity) are α=0 rad, ω=0 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 four-bar linkage, the differences in equations of motion do not significantly influence the results. The midspan bending strains for the coupler (Figure3) 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.Figure 3
Midspan bending strains for a coupler with crank angle.Figure 4
Midspan bending strains for a follower with crank angle.Figure 5
Difference between two cases in midspan coupler strains as a function of the crank angle.Figure 6
Difference between two cases in midspan coupler strains as a function of the crank angle.The difference in equations of motion (22) and (25) representing Cases 1 and 2, respectively, is only in the inertia force vector on the right-hand sides of these equations. For (22) this vector is equal to Fc0=[M0]{x¨0}, and for (25) Fc=[M]{x¨0}. 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.,  Fc-Fc0) is presented; its value is approximately 1% of the value of component forces.Figure 7
Inertia force versus the crank angle.Figure 8
Difference in inertia force between two cases versus the crank angle.


## 3.2. Two-Link Planar Manipulator

The second example concerns the open-loop mechanism of a two-link planar manipulator (Figure9). 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 (m3 = 0.5 kg), the moment of inertia at the hub I3=0.04 kg m2, and the nominal payload  m4=0.3 kg are taken into account.Table 2
Two-link manipulator parameters [21].
Parameters
Link I
Link II
Length, m.
L
1
=
1
L
2
=
1
Cross-section area, m2
F
1
=
4
×
10
-
4
F
2
=
4
×
10
-
4
Moment of inertia, m4
I
1
=
16
/
3
×
10
-
8
I
2
=
16
/
3
×
10
-
8
Modulus of elasticity, N/m2 E=2.1×1011
Mass density, kg/m3
ρ=7800Figure 9
A two-link planar manipulator.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:M1=-M2=10Nm. The vibration analysis is conducted by solving equations of motion (22) and (25) which correspond to the two cases considered. The initial conditions are α=0 rad, β=0.5 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{x} in (22) or (25) consists of 12 elements: displacements u  and v in the X and Y directions, respectively of moving nodes 2, 3, 4, and 5, and nodal deformation angles Θ:
(30){x}T=[u2,w2,Θ2,u3,w3,Θ3,u4,w4,Θ4,u5,w5,Θ5].
From the reason described in Section 3.1 the nodal displacement vector for the rigid-body motion in the case of (22) is equal to
(31){x0}T=[u02,w02,u03,w03,u04,w04,u05,w05]
and in the case of (25) of the form
(32){x0}T=[u02,w02,0,u03,w03,0,u04,w04,0,u05,w05,0].
The system matrices [M], [C], and [K] in (22) and (25) are of dimensions 12×12, and the matrix [M0′] in (22) is of dimensions 12×8.The horizontal and vertical displacements of manipulator’s tip are presented in Figures10 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 Y 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.Figure 10
Displacement of manipulator’s tip inX direction for cases: I (solid line) and II (dashed line) versus time.Figure 11
Displacement of manipulator’s tip inY direction for cases: I (solid line) and II (dashed line) versus time.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 Figure12. 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.Figure 12
Displacement of manipulator’s tip (no payload) inY direction for cases: I (solid line) and II (dashed line) versus time.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 strain-displacement relationship is taken into account [14, 20]; the additional geometric stiffness matrix [kG] appears in dynamic equations, and it should be added to the element stiffness matrix [k]. For a beam element this matrix can be presented in the following form [23]:
(33)[kG]=FL[0065Symm.0L102L21500000-65-L100650L10-L2300-L102L215],
where  F 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 two-link manipulator. The results of midspan bending strains for link I and II are presented in Figures13 and 14, respectively, for linear (solid line) and nonlinear (dashed line) models.Figure 13
Midpoint bending strains for link I for linear (solid line) and nonlinear (dashed line) cases versus time.Figure 14
Midpoint bending strains for link II for linear (solid line) and nonlinear (dashed line) cases versus time.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 kineto-elastodynamic 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 rigid-body 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 rigid-body acceleration is different for both approaches.The method presented here is illustrated for the cases of planar open- and closed-loop linkages. The results of numerical simulation show that in the case of considered four-bar linkage the differences in equations of motion do not considerably influence the results. An opposite observation is made in the case of an open-loop mechanism—two-link 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, steady-state response is analysed) the two approaches presented in this paper should give quite converged results.



---

*Source: 648706-2013-04-08.xml*