Kane’s Method-Based Simulation and Modeling Robots with Elastic Elements, Using Finite Element Method

Abstract

The Lagrange’s equation remains the most used method by researchers to determine the finite element motion equations in the case of elasto-dynamic analysis of a multibody system (MBS). However, applying this method requires the calculation of the kinetic energy of an element and then a series of differentiations that involve a great computational effort. The last decade has shown an increased interest of researchers in the study of multibody systems (MBS) using alternative analytical methods, aiming to simplify the description of the model and the solution of the systems of obtained equations. The method of Kane’s equations is one possibility to do this and, in the paper, we applied this method in the study of a MBS applying finite element analysis (FEA). The number of operations involved is lower than in the case of Lagrange’s equations and Kane’s equations are little used previously in conjunction with the finite element method (FEM). Results are obtained regardless of the type of finite element used. The shape functions will determine the final form of the matrix coefficients in the equations. The results are applied in the case of a planar mechanism with two degrees of freedom.

1. Introduction

The dynamic analysis of a multibody system (MBS) with elastic elements requires, as a first step, time response of a single finite element, obtained via the equations of motion. These equations will depend on the model chosen, the shape function used and the particular type of motion. Once these equations are determined, we move on to the next necessary steps, which consist in refer these equations to a unique coordinate system, assembling the equations of motion and finally solving them. Analytical mechanics offers a large number of formalisms for determining the equations of motion. In general, classical methods, well verified, are used for all these steps. Obtaining the equations of motion is generally a laborious and not a simple step, due to the volume of calculations involved [1]. Some researches presents a comparison between various analytical methods [2,3,4].

For this reason, choosing a method that gives us advantages in terms of the volume and time of calculations is an important objective for the researchers. In general, in the case of using the finite element method (FEM), the method used almost exclusively in all researches is the method of Lagrange’s equations. The advantage of using this method lies in the homogeneity of the writing, the ease of approaching large systems, with large number of degrees of freedom (DOF) and, an important thing, a familiarity of the researchers with this method. However, analytical mechanics offers equivalent methods of writing the equations of motion. For the same purpose, Gibbs–Appell equations, Hamilton equations, Maggi’s equations, Jacobs equations and other equivalent forms can be used [1,5]. Either of these methods ultimately gives the researchers the final form of these equations, which is the same in all cases. For a researcher, however, the problem is posed to determine which of these methods is the most appropriate to the studied problem, from the point of view of the ease of modeling, of the solving algorithm and the computation time needed to study the studied mechanical system. It is, therefore, natural to look for one among them that would allow us to obtain these equations easily.

Based on the studied literature, one can form the opinion that Lagrange’s equations have the advantage of the researchers’ familiarity with this method and with a relatively simple model [6,7]. When applied to holonomic systems using coordinate derivatives, the momentum form of Kane’s equations is an efficient alternative to Lagrange’s equations, providing a momentum formulation without the need to assemble and differentiate the system kinetic co-energy function [8].

The disadvantage is the need to determine Lagrange’s multipliers from the system of differential and algebraic equations (DAE) obtained, which, in the case of large systems, can lead to a large number of calculations and large computation times [9]. The Gibbs–Appell [10,11] method seems to have the advantage of economical computing time dedicated to modeling and computing itself. The disadvantage is the introduction of a notion with which the researchers are a little used, namely the energy of accelerations. Maggi’s equations become very useful when analyzing MBS in which the links between the elements are non-holonomic constraints. In recent years, the method has begun to be applied more often by researchers, especially due to the need to study robots and manipulators widely applied in the manufacturing industry. Hamiltonian formalism can be profitable from the point of view of numerical calculation, since its application ultimately leads to a system of differential equations of the first order, as opposed to the other methods that ultimately give differential equations of the second order. This can be an advantage in the case of numerical calculation because the computing subroutines solve numerically differential equations of the first order. However, the complexity of intermediate calculations is high, which represents a disadvantage of the method. The method of Kane’s equations is equivalent to Maggi’s formalism and has begun to be used more in the last decade, also related to the need for robot and manipulator analysis. In this paper, we try to apply the Kane equations to solve the problem of modeling and obtaining the motion equations in the case of the use of FEM in the dynamical analysis of a MBS with elastic elements. Recent years have shown a high interest in this method.

We assume that in the coming years all these methods mentioned here, but also others, equivalent, from analytical mechanics, will be re-evaluated in the context of the numerical possibilities to solve such problems and will be used in commercial software those methods that will prove those greater advantages. It is assumed that the fields of applicability of the mentioned methods will be expanded.

In Section 2 of the paper was given an overview of the results obtained using various formalisms of analytical mechanics. The classical methods with their advantages and disadvantages were presented. In order to present the original results, which can be found in Section 4, it is necessary to present some preliminary notions. This presentation can be found in Section 3. To emphasize the possibility of using the method for such problems, an application to a mechanism with two degrees of freedom is presented. The conclusion and discussion will be presented in Section 6.

2. A Brief Overview of the Applied Formalism in (FEA) of (MBS)

Lagrange’s method: The Lagrange’s equations (LE) method is, at the moment, the most used method for writing differential equations of motion for mechanical systems with a high number of degrees of freedom (DOF). The construction of Lagrange’s function involves the calculation of kinetic energy, potential energy and mechanical work of a mechanical system. These notions are often used in applications and researchers are familiar with them. It is also the reason why this method is mainly used in the works dedicated to the analysis of elastic multi-body systems with the finite element method (FEM). Most of the articles in the field used the method in the analysis of one, two-and three-dimensional mechanical systems, regardless of the type of finite element chosen for the study.

Using the LE method the motion equations of a single finite element is determined. This represents a first stage in an analysis of a multibody system (MBS) using FEM [6,7,8].

The next step it is represented by the procedure of assembling the equations of motion, which involves the elimination of Lagrange’s multipliers (depending on the liaison forces in the common nodes to several elements). Numerous equivalent methods of analytical mechanics can be used for the assembly procedure, such as: Lagrange’s equations method, Maggi’s equations method, Gibbs–Appell formalism, Kane’s equations method or other formulations. As in the case of determining the equations of motion for a single finite element, the LE method was used too, with predilection [5].

The first papers that studied these systems using FEM dealt with systems with elastic elements that can be modeled with one-dimensional finite elements. In the 1970s there are lot of paper analyzing such kind of problems and developing the domain [11,12,13,14,15,16]. Naturally, the method was then extended for two and three-dimensional elements [17,18,19].

Different aspects of the behavior of such systems have been studied by researchers in recent years [20]. We can mention that almost all researches used the LE method, which proved to be a familiar and convenient method to solve the problems. However, there are equivalent formulations in analytical mechanics and we can put the question, to what extent these methods can replace LE in the analysis of systems with elastic elements [21].

Gibbs–Appell equations. Although known since the beginning of the century XX, equations called Gibbs–Appell (GA) and discovered independently by Gibbs (1879) [22] and Appell (1899) [23] have been rarely used in mechanical systems study. However, they have obvious advantages in terms of the number of equations to be written. This is obviously smaller if compared to LE and it results in advantages in terms of modeling effort and calculation time required. A main disadvantage is the fact that the researchers are less familiar to use GA equations. Usually the GA method represents a more convenient procedure to study non-holonomic systems. Basically the method consists in replacing the well-known Lagrangian with the “energy of accelerations”.

In the last period, multiplying the studies regarding the multibody systems, with high working speed and loaded with great forces, the advantages of applying the Gibbs–Appell equations have been found. There are some papers using this method in the field of MBS [24,25,26]. The equations that are obtained by applying this method coincide with those obtained by applying the method of LE but the number of operations required to be performed is smaller. This method is an application of Gauss’s principle of least constraint.

This method has proven useful for a wide range of problems. It has been successfully used in the dynamic analysis of some rigid systems and considering also the quasi-velocities [27]. The equations of motion can be obtained in this formalism more easily for linear and nonlinear systems with holonomic or non-holonomic liaisons. It becomes possible to eliminate Lagrange’s multipliers and, as a consequence, to reduce the number of unknowns and the calculations to be performed. In [28], such an example is presented for a robot with flexible elements. For the advantages offered the method tends to become a procedure. There are several papers that present the advantage of using this formalism, as a result of the fact that the number of unknowns is reduced, by eliminating Lagrange’s multipliers [29,30,31,32].

Maggi’s formalism: Maggi’s equations, which have been applied little and are relatively recent (1896), are an alternative formalism to the other formalisms in analytical mechanics, presented in several papers [1,33].

Recent years have shown a great interest of researchers in applying Maggi’s equations to different engineering problems. In the case of modern technical systems, used in industry, such as robots and manipulators, there are situations in which the control system must be developed and designed [34,35]. If we work with nonlinear systems, the technique of linearizing the feedback will naturally lead to the application of these equations. Maggi’s formalism was found to be a simple and stable way to determine the equations that provide the dynamic response of an MBS with constraints. An important step in applying this method is the choice of independent coordinates. The method is considered by some researchers [36,37,38] as the most efficient way to solve Lagrange’s index-1 equations. The existence of the non-holonomic liaisons leads, in a naturally way, to the use of Maggi’s equations. For systems with a large number of DOF, this approach can prove useful from the point of view of computing time. Despite these advantages, the method is a little familiar to the researchers and was, as a consequence, little used. The major advantage is due to the possibility of approaching large systems, with those appearing in current engineering applications. In these situations, if a classical calculation method, such as Newton–Euler or Lagrange, is used, a step would be to determine the liaison forces. For large systems, this can lead to laborious calculations and, as a consequence, a high modeling and computing effort. Maggi’s method also offers a justification of the classical assembly methods, applied empirically in finite element analysis (FEA).

Hamilton’s method: The use of one of the classical methods of determining the equations of motion of an elastic system naturally leads to a system of differential equations of the second order in which the unknowns are represented by the generalized coordinates, q_i. To be applied, a commercial program for solving this second order system with n unknowns implies its writing, by introducing additional unknowns (represented by the generalized velocities), in the form of a system of 2n differential equations of the first order. Hamiltonian mechanics achieves this precisely, offering as a result of modeling a system of differential equations from first order with 2n unknown. The first n unknowns in this case are the generalized coordinates, the others being the generalized impulse momentum.

Introducing this notion, the conjugated canonical moment represents a difference from Lagrange’s method (or an equivalent method from the point of view of analytical mechanics). It follows, therefore, that Hamilton’s formalism offers a system of first order equations. Hamilton method could have the advantage of giving us a system of first order equations, a system which can be used directly for numerical solving, without the need for prior processing of equations [5].

Kane’s equations: In the last decades, Kane’s method and equations have been successfully applied to the modeling and numerical analysis of MBS. The motion equations proposed by Kane have been successfully applied to numerical analysis and simulation of multibody systems [39,40]. However, their use becomes problematic if the MBS is composed of a large number of bodies between which there are many links. For such a situation, the complexity of the equations of motion increases substantially and the determination of the response over time becomes much more difficult. Some examples are presented in [41,42,43,44,45]. The Kane method eliminates the disadvantages of methods of classical mechanics (Newton–Euler and Langrange) and can be used for both holonomic and non-holonomic constraint systems [45,46,47]. This method is also called the Lagrange form of d Alembert’s principle.

Generally, the equations Maggi and Kane are associated, representing essentially the same method. Kane’s equations have been shown to be equivalent to Maggi’s equations. These methods are very good when the constraints are non- holonomic, but may not work in the case of systems with holonomic constraints [33].

One of the first applications of Kane’s equations was made to flexible complex structures [48]. Their application has been found to allow efficient computation of such systems, especially in robots and manipulators [49,50]. The equations that are finally obtained will be equal with the DOF of the system. They do not contain the forces of liaisons. Once integrated the equations of motion, the liaison forces (if their knowledge is necessary) can be obtained, by relatively simple calculations [50].

To our knowledge, Kane’s method was not yet applied to obtain the dynamic response of a finite element in the case of the study of a MBS with elastic elements. This is what we are trying to accomplish in the present work.

Some interesting application of Kane equation can be found in [51,52].

3. Preliminary Kinematics and Kinetics of Finite Elements

In order to be able to apply Kane’s equations to a single finite element, it is necessary to know some kinematic elements. The speeds and accelerations of the different points of the finite element must be determined according to the nodal coordinates, considered as independent coordinates. This has been studied in the specialized literature [6,17,40], which is why we present the relationships that we will use in the paper. The displacement field is described by a set of shape functions.

The velocity of a point M of an elastic solid (Figure 1) is offered by [17,18]:

$${\left\{v_{M^{\prime }}\right\}}_G={\left\{{\dot{r}}_{M^{\prime }}\right\}}_G={\left\{v_O\right\}}_G+\left[\dot{R}\right]{\left\{r\right\}}_L+\left[\dot{R}\right]\left[N\right]{\left\{\delta \right\}}_L+\left[R\right]\left[N\right]{\left\{\dot{\delta }\right\}}_L$$

(1)

where ${\left\{v_O\right\}}_G={\left\{{\dot{r}}_O\right\}}_G$ it represents the velocity vector of the origin of the mobile reference system with component expressed in the global system, ${\left\{r\right\}}_L$ the position vector of the current point with the coordinates expressed in the mobile system, ${\left\{\delta \right\}}_L$ the vector of nodal coordinates, $\left[N\right]$ the shape function matrix ($\left\{u\right\}=\left[N\right]{\left\{\delta \right\}}_L$) and $[R]$ represents the rotation matrix that makes the transition from the mobile reference system to the fixed reference system. $[\dot{R}]$ is the derivative of $\left[R\right]$ with respect to the time.

In the paper, the quantities that are expressed in the global coordinate system will be indexed with G, and if they are expressed in the local coordinate system they will be noticed by L. The transformation of a vector is performed by the matrix: $\left[R\right]$:

$${\left\{t\right\}}_G=\left[R\right]{\left\{t\right\}}_L$$

(2)

The matrix $\left[R\right]$ is orthonormal, thus, it results:

$$\left[R\right]{\left[R\right]}^T={\left[R\right]}^T\left[R\right]=\left[E\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$$

(3)

By differentiating (3) it obtains:

$$\left[\dot{R}\right]{\left[R\right]}^T+\left[R\right]{\left[\dot{R}\right]}^T=\left[\dot{R}\right]{\left[R\right]}^T+{\left(\left[\dot{R}\right]{\left[R\right]}^T\right)}^T=0$$

(4)

and:

$${\left[\dot{R}\right]}^T\left[R\right]+{\left[R\right]}^T\left[\dot{R}\right]={\left[\dot{R}\right]}^T\left[R\right]+{\left({\left[\dot{R}\right]}^T\left[R\right]\right)}^T=\left[0\right]$$

(5)

The notations:

$${\left[\omega \right]}_G=\left[\dot{R}\right]{\left[R\right]}^T=-\left[R\right]{\left[\dot{R}\right]}^T=\left[\begin{array}{ccc}0& -{\omega }_Z& {\omega }_Y\\ {\omega }_Z& 0& -{\omega }_X\\ -{\omega }_Y& {\omega }_X& 0\end{array}\right]$$

(6)

and:

$${\left[\omega \right]}_L={\left[R\right]}^T{\left[\omega \right]}_G\left[R\right]={\left[R\right]}^T\left[\dot{R}\right]=-{\left[\dot{R}\right]}^T\left[R\right]=\left[\begin{array}{ccc}0& -{\omega }_z& {\omega }_y\\ {\omega }_z& 0& -{\omega }_x\\ -{\omega }_y& {\omega }_x& 0\end{array}\right]$$

(7)

indicate the skew symmetric operator angular velocity, expresses in global reference frame ${\left[\omega \right]}_G$, and in the local reference frame ${\left[\omega \right]}_L$. The angular velocity vector of the solid is $\overline{\omega }$. There are the relations:

$${\left[\omega \right]}_G=-{\left[\omega \right]}_G^T\mathrm{and}{\left[\omega \right]}_L=-{\left[\omega \right]}_L^T$$

(8)

The operator angular acceleration is:

$${\left[\epsilon \right]}_G={\left[\dot{\omega }\right]}_G=\left[\ddot{R}\right]{\left[R\right]}^T+\left[\dot{R}\right]{\left[\dot{R}\right]}^T=-\left[\dot{R}\right]{\left[\dot{R}\right]}^T-\left[R\right]{\left[\ddot{R}\right]}^T=\left[\begin{array}{ccc}0& -{\epsilon }_Z& {\epsilon }_Y\\ {\epsilon }_Z& 0& -{\epsilon }_X\\ -{\epsilon }_Y& {\epsilon }_X& 0\end{array}\right]$$

(9)

in the global reference frame and:

$${\left[\epsilon \right]}_L={\left[\dot{\omega }\right]}_L={\left[\dot{R}\right]}^T\left[\dot{R}\right]+{\left[R\right]}^T\left[\ddot{R}\right]=-{\left[\dot{R}\right]}^T\left[\dot{R}\right]-{\left[\ddot{R}\right]}^T\left[R\right]=\left[\begin{array}{ccc}0& -{\epsilon }_z& {\epsilon }_y\\ {\epsilon }_z& 0& -{\epsilon }_x\\ -{\epsilon }_y& {\epsilon }_x& 0\end{array}\right]$$

(10)

in a local reference frame. After some elementary calculus it results:

$$\begin{array}{l}\left[\ddot{R}\right]{\left[R\right]}^T={\left[\epsilon \right]}_G-\left[\dot{R}\right]{\left[\dot{R}\right]}^T={\left[\epsilon \right]}_G+{\left[\omega \right]}_G{\left[\omega \right]}_G\\ =\left[\begin{array}{ccc}-\left({\omega }_Y^2+{\omega }_Z^2\right)& -{\epsilon }_Z+{\omega }_X{\omega }_Y& -{\epsilon }_Y+{\omega }_X{\omega }_Z\\ {\epsilon }_Z+{\omega }_X{\omega }_Y& -\left({\omega }_Z^2+{\omega }_X^2\right)& -{\epsilon }_X+{\omega }_Y{\omega }_Z\\ {\epsilon }_Y+{\omega }_X{\omega }_Z& {\epsilon }_X+{\omega }_Y{\omega }_Z& -\left({\omega }_X^2+{\omega }_Y^2\right)\end{array}\right]\end{array}$$

(11)

$$\begin{array}{l}{\left[R\right]}^T\left[\ddot{R}\right]={\left[\epsilon \right]}_L-{\left[\dot{R}\right]}^T\left[\dot{R}\right]={\left[\epsilon \right]}_L+{\left[\omega \right]}_L{\left[\omega \right]}_L\\ =\left[\begin{array}{ccc}-\left({\omega }_y^2+{\omega }_z^2\right)& -{\epsilon }_z+{\omega }_x{\omega }_y& -{\epsilon }_y+{\omega }_x{\omega }_z\\ {\epsilon }_z+{\omega }_x{\omega }_y& -\left({\omega }_z^2+{\omega }_x^2\right)& -{\epsilon }_x+{\omega }_y{\omega }_z\\ {\epsilon }_y+{\omega }_x{\omega }_z& {\epsilon }_x+{\omega }_y{\omega }_z& -\left({\omega }_x^2+{\omega }_y^2\right)\end{array}\right]\end{array}$$

(12)

4. Kane’s Formalism Applied to MBS

The Kane’s equations are presented in the Appendix A. Starting from the basic equations:

$${\displaystyle \sum _{i=1}^N\left({\overline{F}}_i-m_i{\overline{a}}_i\right)\delta }{\overline{r}}_i=0,$$

(13)

we obtain, for an elastic finite element considered as a solid, the system of equations:

$${\displaystyle \sum _{i=1}^N{\overline{F}}_i\frac{\partial {\overline{v}}_i}{\partial {\dot{q}}_k}}={\displaystyle \underset{V}{\int }\overline{a}}\frac{\partial \overline{v}}{\partial {\dot{q}}_k}dmk=\overline{1,n}$$

(14)

The acceleration of an arbitrary point is obtained by differentiating the expression of velocity (1):

$${\left\{a_{M^{\prime }}\right\}}_G={\left\{{\ddot{r}}_{M^{\prime }}\right\}}_G={\left\{a_O\right\}}_G+\left[\ddot{R}\right]{\left\{r\right\}}_L+\left[\ddot{R}\right]\left[N\right]{\left\{\delta \right\}}_L+2\left[\dot{R}\right]\left[N\right]{\left\{\dot{\delta }\right\}}_L+\left[R\right]\left[N\right]{\left\{\ddot{\delta }\right\}}_L$$

(15)

Differentiating (1), we have:

$$\frac{\partial {\left\{{v^{\prime }}_M\right\}}_G}{\partial {\left\{\dot{\delta }\right\}}_L}=\left[R\right]\left[N\right]\mathrm{and}{\left(\frac{\partial {\left\{{v^{\prime }}_M\right\}}_G}{\partial {\left\{\dot{\delta }\right\}}_L}\right)}^T={\left[N\right]}^T{\left[R\right]}^T$$

(16)

The expression $\frac{\partial \left\{T\right\}}{\partial \left\{\delta \right\}}$, if $\left\{T\right\}=\left[A\right]\left\{\delta \right\}$ it is $\left[A\right]$, result used previously. Now we can written:

$$\begin{array}{c}{\left(\frac{\partial {\left\{v_{M^{\prime }}\right\}}_G}{\partial {\left\{\dot{\delta }\right\}}_G}\right)}^T{\left\{a_{M^{\prime }}\right\}}_G={\left[N\right]}^T{\left[R\right]}^T\left({\left\{a_O\right\}}_G+\left[\ddot{R}\right]{\left\{r\right\}}_L+\left[\ddot{R}\right]\left[N\right]{\left\{\delta \right\}}_L+2\left[\dot{R}\right]\left[N\right]{\left\{\dot{\delta }\right\}}_L+\left[R\right]\left[N\right]{\left\{\ddot{\delta }\right\}}_L\right)=\\ =\left({\displaystyle \underset{V}{\int }{\left[N\right]}^T\rho dV}\right){\left\{a_O\right\}}_L+\left({\displaystyle \underset{V}{\int }{\left[N\right]}^T{\left[R\right]}^T\left[\ddot{R}\right]{\left\{r\right\}}_L^T\rho dV}\right)+\left({\displaystyle \underset{V}{\int }{\left[N\right]}^T{\left[R\right]}^T\left[\ddot{R}\right]\left[N\right]\rho dV}\right){\left\{\delta \right\}}_L+\\ +2\left({\displaystyle \underset{V}{\int }{\left[N\right]}^T{\left[R\right]}^T\left[\dot{R}\right]\left[N\right]\rho dV}\right){\left\{\dot{\delta }\right\}}_L+\left({\displaystyle \underset{V}{\int }{\left[N\right]}^T\left[N\right]\rho dV}\right){\left\{\ddot{\delta }\right\}}_L=\\ =\left({\displaystyle \underset{V}{\int }\rho {\left[N\right]}^{T}dV}\right){\left\{a_O\right\}}_L+{\displaystyle \underset{V}{\int }{\left[N\right]}^T\left({\left[\epsilon \right]}_L+{\left[\omega \right]}_L{\left[\omega \right]}_L\right){\left\{r\right\}}_L\rho dV}+\\ +\left({\displaystyle \underset{V}{\int }{\left[N\right]}^T\left({\left[\epsilon \right]}_L+{\left[\omega \right]}_L{\left[\omega \right]}_L\right)\left[N\right]\rho dV}\right){\left\{\delta \right\}}_L+2\left({\displaystyle \underset{V}{\int }\rho {\left[N\right]}^T{\left[\omega \right]}_L\left[N\right]dV}\right){\left\{\dot{\delta }\right\}}_L+\left({\displaystyle \underset{V}{\int }{\left[N\right]}^T\left[N\right]\rho dV}\right){\left\{\ddot{\delta }\right\}}_L\hfill \\ =-\left[m_O^i\right]{\left\{a_O\right\}}_L-\left\{Q^i(\omega )\right\}-\left\{Q^i(\epsilon )\right\}+\left[k\left(\omega \right)\right]{\left\{\delta \right\}}_L+\left[k\left(\epsilon \right)\right]{\left\{\delta \right\}}_L+\left[c\left(\omega \right)\right]{\left\{\dot{\delta }\right\}}_L+\left[m\right]{\left\{\ddot{\delta }\right\}}_L\hfill \end{array}$$

(17)

The following notations were made:

$$\begin{gathered} \left\{Q^i\left(\epsilon \right)\right\}={\displaystyle \underset{V}{\int }{\left[N\right]}^T{\left[\epsilon \right]}_L{\left\{r\right\}}_L\rho dV};\left\{Q^i\left(\omega \right)\right\}={\displaystyle \underset{V}{\int }{\left[N\right]}^T{\left[\omega \right]}_L{\left[\omega \right]}_L{\left\{r\right\}}_L\rho dV}; \\ \left[k\left(\epsilon \right)\right]={\displaystyle \underset{V}{\int }{\left[N\right]}^T{\left[\epsilon \right]}_L\left[N\right]\rho dV};\left[k\left(\omega \right)\right]={\displaystyle \underset{V}{\int }{\left[N\right]}^T{\left[\omega \right]}_L{\left[\omega \right]}_L\left[N\right]\rho dV}; \\ \left[c\right]=2{\displaystyle \underset{V}{\int }\rho {\left[N\right]}^T{\left[\omega \right]}_L\left[N\right]dV};\left[m\right]={\displaystyle \underset{V}{\int }{\left[N\right]}^T\left[N\right]\rho dV};\left[m_O^i\right]={\displaystyle {\int }_V\rho \left[N\right]d}V \end{gathered}$$

(18)

The generalized forces $F_k$ act in nodal points of the finite element. The velocities of these points are among the components of the velocities vector ${\left\{\dot{\delta }\right\}}_G$. It obtains:

$${\left(\frac{\partial {\left\{v_{M^{\prime }}\right\}}_G}{\partial {\left\{\dot{\delta }\right\}}_G}\right)}^T\left\{F_k\right\}\rho dV=\left[k\right]{\left\{\delta \right\}}_L-{\left\{Q^{ext}\right\}}_L-{\left\{Q^{liaison}\right\}}_L$$

(19)

${\left\{Q^{ext}\right\}}_L$ represents the generalized nodal forces vector acting in the node, $\left[k\right]{\left\{\delta \right\}}_L$ represents the elastic force between given by the stiffness matrix $\left[k\right]$, ${\left\{Q^{liaison}\right\}}_L$ represents liaison forces in the nodal points (these forces there are in the common nodes of two neighboring finite elements).

It results the motion equations:

$$\left[m\right]{\left\{\ddot{\delta }\right\}}_L+\left[c\right]{\left\{\dot{\delta }\right\}}_L+\left(\left[k\right]+\left[k\left(\epsilon \right)\right]+\left[k\left(\omega \right)\right]\right){\left\{\delta \right\}}_L=\left[m_O^i\right]{\left\{a_O\right\}}_L+{\left\{Q^i\left(\epsilon \right)\right\}}_L+{\left\{Q^i\left(\omega \right)\right\}}_L+{\left\{Q^{ext}\right\}}_L+{\left\{Q^{nodal}\right\}}_L$$

(20)

5. An Application to a Two Degrees of Freedom Mechanism Used in a Wind Water Pump

We consider the mechanical system of Figure 2, having two degrees of freedom [53]. The B C’ beam is an elastic bar, the other elements of the mechanism being sufficiently massive to be considered rigid (experimental prototype is presented in Figure 3).

In a first approximation, to obtain the liaison forces, appearing in B and C ‘, the mechanism is considered consisting of rigid elements considered and, in this hypothesis, the forces are calculated.

We write the equation of vector contour for the mechanism:

$$\overrightarrow{AB}+\overrightarrow{B{C}^{\prime }}+\overrightarrow{C^{\prime }C}+\overrightarrow{CD}+\overrightarrow{DE}+\overrightarrow{EP}+\overrightarrow{PA}=0$$

(21)

With the notations used in the Figure 2 the relations are obtained:

$$\begin{gathered} l_1\cos \alpha \overrightarrow{i}+l_1\sin \alpha \overrightarrow{j}+l_2\cos \beta \overrightarrow{i}-l_2\sin \beta \overrightarrow{j}+l_3\cos \theta i+l\sin \theta j- \\ -l_4\sin \theta \overrightarrow{i}+l_4\cos \theta \overrightarrow{j}+e\overrightarrow{j}-X_D\overrightarrow{i}-D\overrightarrow{j}=0 \end{gathered}$$

(22)

If the terms are grouped according to the directions of the two verses, the boundary conditions are obtained:

$$\{\begin{array}{l}{l}_1\cos \alpha +l_2\cos \beta +l_3\cos \theta -l_4\sin \theta -X_D=0\\ l_1\sin \alpha -l_2\sin \beta +l_3\sin \theta +l_4\cos \theta +e-D=0\end{array}$$

(23)

If we differentiate the obtained system (23) (see Appendix B) it results:

$$\left\{a\right\}=\left\{\begin{array}{c}{\ddot{x}}_{A_1}\\ {\ddot{y}}_{A_1}\\ {\epsilon }_1\\ {\ddot{x}}_{B_1}\\ {\ddot{y}}_{B_1}\\ {\epsilon }_2\\ {\ddot{x}}_C\\ {\ddot{y}}_C\\ {\epsilon }_3\\ {\ddot{x}}_D\end{array}\right\}=\left[C_1\right]\left\{\begin{array}{c}{\epsilon }_1\\ {\ddot{x}}_D\end{array}\right\}+\left[C_2\right]\left\{\begin{array}{c}{\omega }_1^2\\ {\omega }_1{\dot{x}}_D\\ {\dot{x}}_D^2\end{array}\right\}$$

(24)

where $\left\{a\right\}$ is the vector of the accelerations of the mass centers of the elements of the mechanism and the angular accelerations, the two independent generalized coordinates are the rotation angle of the crank AB, $\alpha$ and the displacement of the slide $x_D$ and it is noted: $\dot{\alpha }={\omega }_1$, ${\dot{\omega }}_1={\epsilon }_1$. It is, therefore, hypothesized that the elasticity of the elements of the mechanism does not influence its rigid movement. The mechanism will be decomposed into four bodies (Figure 4) in order to determine the equations of motion:

Because there are many parameters that can vary, integrating the equations of motion of the mechanism, considering all the rigid elements is difficult and will lead to results far from reality. This is due to the fact that a number of factors such as joints that appear in the joints, the engine couple, and the frictions that occur in the joints which cannot be accurately determined. For this reason, we will use a method of recording the rigid movement of the mechanism through an optical system, determining from these records the law of motion of the two independent coordinates [53,54,55,56]. Knowing these quantities, the forces of liaison in joints can be determined. These forces are necessary for us to be able to write the equations of motion for the BC’ elastic beam, determined previously, formally, using Kane’s equations. After these a computation of the free vibrations of the obtained equation system will be made.

Relationships (13–16) show that in order to determine the equations of motion of a finite element it is necessary to know the velocities and accelerations of a point according to the nodal coordinates. Based on these, Kane’s equations can be written, using simple matrix operations. In the case of applying Lagrange’s equations, it is necessary to write first the kinetic energy, which already involves multiplication operations between numerous terms, followed by several differentiations, which again involves more numerous operations. This is why Kane’s equations can be an alternative to simplify calculations for obtaining equations of motion.

Considering the lever from Figure 5 it is possible to write the Equation (20) for this and to analyze obtained equations. It is possible to integrate the differential equations obtained using a commercial soft. In our problem we did a qualitative analysis of the mechanical system, computing the eigenpairs of the beam integrated into mechanism. The first two eigenvalues and their variation versus time are presented in Figure 6.

Figure 6 represents the first two eigenvalues for the beam BC’, divided in 10 finite elements, for an average angular speed of the engine of 140 rpm.

It is difficult to make a comparison between the different methods of determining the equations of motion of an elastic element of a MBS system. However, we will try to compare the number of differentiations needed to obtain these equations. Thus,

Table 1. Comparison between the Lagrange, Gibbs–Appell [38] and Kane methods. Number of differentiations.

Number of Finite Elements	Lagrange	Gibbs–Appell	Kane
5	288	120	10
10	528	220	20
15	768	320	30
20	1008	420	40
25	1248	520	50
30	1488	620	60
40	1968	820	80

presents a comparison between the number of differentiations required when using the Lagrange method, the number of differentiations in the case of the Gibbs–Appell formalism [38] and if the Kane equations are applied. Starting from the velocity formula, it is found that in the case of this method we need a differentiation of the acceleration as a function of time and a partial derivative, so two differentiation operations. It is observed that the number of these operations is substantially reduced compared to Lagrange’s method. In the case of illustrative application presented in the paper, it does not seem to be a major advantage. But if we refer to a complex application, with a large number of finite elements used, the reduction can be significant and the calculation time and effort can be significantly reduced. We note that after the matrix coefficients will be obtained, regardless of the method by which they were obtained, the following procedures are identical in all cases, so they will not influence the necessary calculation times.

6. Conclusions and Discussions

The study of a MBS with elastic elements using FEM leads to the necessity of determining the equations of motion for a single finite element, taking into account its type and the interpolation functions used. It is the main step in the first stage of this study. This can be done using several methods of analytical mechanics, equivalent between them and with which the same results are obtained. We can use Lagrange’s equations, Newton–Euler equations, Kane’s method or equivalent formulations of these formalisms. The use of any of these methods presents, from the user’s point of view, advantages and disadvantages, which have been highlighted by different authors and which have been discussed in the Introduction.

The formalism Newton–Euler method allows, for example, to obtain formal equations independent of geometric, inertial properties or liaisons. But it also has a disadvantage, which for some systems with large numbers of DOF becomes significant, namely that the liaison forces and moments must be computed. For relatively simple systems, however, the method is useful and intuitive. However, if we study a complex system with a high number of DOF, Lagrange’s method becomes convenient, which also presents the advantage of a familiarity of researchers with it. The method is suitable for generalizations and large systems. The disadvantage is to determine the generalized liaison forces. Due to these advantages the method us preferred by some software developers (ADAMS, DADS, DYMAC). But there are also developers who use equivalent formulations of this formalism (SD-EXACT, NBOD2, SD/FAST) materialized in Kane’s equations, which essentially represent a formalism equivalent to Maggi’s equations. The application of these equations in the modeling of FEM has not yet been used, to our knowledge. As a result, we set out to see to what extent these equations can be applied to the modeling of an elastic linear finite element, resulting in a computational economy.

In conclusion, the paper aims to present how Kane’s equations can be applied to obtain the equations of motion for a beam, as an elastic element of a mechanism, which is in planar motion, if the FEM is used. Kane’s method has also been applied to MBS but using classical methods of analysis. The paper aims to use Kane’s equations in the procedures involved using FEA. An application comes to exemplify the method for the case of a real mechanism.

The paper studied the application of this formalism in the case of a MBS system with elastic elements, to solve the main stage in such an approach, namely the determination of the equations of motion for a finite element. These equations are dependent on the type of finite element chosen and the type of motion (one, two or three-dimensional). This opens the possibility of applying this method to the development of software at the finite element level. Kane’s equations are an economical and simple alternative to the problem of determining the dynamic response of a single finite element. The method is suitable to apply to MBS systems with non-holonomic constraints.