Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 648706, 10 pages

http://dx.doi.org/10.1155/2013/648706

Research Article

On Equations of Motion of Elastic Linkages by FEM

Michał Hać

Institute of Machine Design Fundamentals, Warsaw University of Technology, Narbutta 84, 02-524 Warsaw, Poland

Received 26 November 2012; Revised 26 February 2013; Accepted 27 February 2013

Academic Editor: M. Onder Efe

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.

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.

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.

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.

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 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 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 Figure 1) is denoted by and can be written as

where is the displacement connected with the rigid-body 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 rigid-body 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 Boolean-type 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 Gibbs-Appel 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 Gibbs-Appel 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 rigid-body translation. Thus the shape function for elastic displacement (15) can be used to determine rigid-body displacement, that is,

where vector 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 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 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]:

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 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 Figure 2.

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 ) 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 rigid-body 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 rigid-body 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 rigid-body 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 rigid-body motion are set to be equal to 0. The rigid-body displacement vector in the case of (25) is as follows:

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 Table 1 [3].

Table 1: Four-bar linkage parameters [3].

Parameters
Crank
Coupler
Follower

Length, m
0.10795
0.2794
0.27051

Cross-section, m²

Bending moments of inertia, m⁴

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  10¹¹ N/m²

Mass density: 2710 kg/m³

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 four-bar 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.

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

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 (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², and the nominal payload  kg are taken into account.

Table 2: Two-link manipulator parameters [21].

Parameters
Link I
Link II

Length, m.

Cross-section area, m²

Moment of inertia, m⁴

Modulus of elasticity, N/m²    

Mass density, kg/m³    

Figure 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: . 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 rigid-body 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.

Figure 10: Displacement of manipulator’s tip in direction for cases: I (solid line) and II (dashed line) versus time.

Figure 11: Displacement of manipulator’s tip in 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 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.

Figure 12: Displacement of manipulator’s tip (no payload) in 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 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 two-link 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.

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.

Appendix

In this appendix the derivation of displacements of a finite element in the rigid-body 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.

Figure 15: Planar rigid element.

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 rigid-body 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 rigid-body displacements for is correct, the result for transverse displacement of the general point for is incorrect.

Using the shape function for the rigid-body 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 rigid-body 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 rigid-body motion is only an approximation and that the shape function for rigid-body 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.

A. Midha, A. G. Erdman, and D. A. Frohrib, “Finite element approach to mathematical modeling of high-speed elastic linkages,” Mechanism and Machine Theory, vol. 13, no. 6, pp. 603–618, 1978.

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.

Z. Yang and J. P. Sadler, “On issues of elastic-rigid coupling in finite element modeling of high-speed machines,” Mechanism and Machine Theory, vol. 35, no. 1, pp. 71–82, 2000.

P. Fanghella, C. Galletti, and G. Torre, “An explicit independent-coordinate formulation for the equations of motion of flexible multibody systems,” Mechanism and Machine Theory, vol. 38, no. 5, pp. 417–437, 2003.

W. Hou and X. Zhang, “Dynamic analysis of flexible linkage mechanisms under uniform temperature change,” Journal of Sound and Vibration, vol. 319, no. 1-2, pp. 570–592, 2009.

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.

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.

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.

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.

X. Zhang, H. Liu, and Y. Shen, “Finite dynamic element analysis for high-speed flexible linkage mechanisms,” Computers and Structures, vol. 60, no. 5, pp. 787–796, 1996.

J. B. Gunda, R. K. Gupta, and R. Ganguli, “Hybrid stiff-string-polynomial basis functions for vibration analysis of high speed rotating beams,” Computers and Structures, vol. 87, no. 3-4, pp. 254–265, 2009.

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.

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.

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. 5-6, pp. 606–617, 2008.

P. Fanghella, C. Galletti, and G. Torre, “An explicit independent-coordinate formulation for the equations of motion of flexible multibody systems,” Mechanism and Machine Theory, vol. 38, no. 5, pp. 417–437, 2003.

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.

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.

M. Hać, “Dynamic analysis of flexible mechanisms by finite element method,” Machine Dynamics Problems, vol. 14, pp. 7–91, 1996.

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.

M. Hać, “On equations of motion of flexible mechanisms,” Mechanics Research Communication, vol. 20, pp. 399–403, 1993.

J. S. Przemieniecki, Theory of Matrix Structural Analysis, McGraw-Hill, London, UK, 1968.