Two-Dimensional Equivalent Models in the Analysis of a Multibody Elastic System Using the Finite Element Analysis

Abstract

Analytical mechanics provides methods for analyzing multibody systems with mathematically equivalent elastic elements. The paper analyzes several of these models, highlighting the advantages and disadvantages offered by each of these methods. The main methods used by the researchers are described in a unitary form, presenting the methods of obtaining the evolution equations in each of these cases, mentioning the strengths and weaknesses of each method. The equations of Lagrange, Gibbs–Appell, Kane, Maggi, and Hamilton are analyzed for the particular case of two-dimensional systems, which present certain particularities that facilitate the analysis.

1. Introduction

The problems of the dynamic analysis of multibody systems (MBSs), arising as a result of the increase in the working speeds of these systems, the forces that require them, and the variety of conditions in which they must work, have determine the results of numerous works related to these problems and the definition of a new field of research, namely, the study of MBS systems. The elasticity of the elements proved to significantly influence the behavior of these systems and their existence cannot be neglected in most MBSs that have elastic elements in their compositions. Additionally, since the finite element method (FEM) is at a very advanced stage in its development, the results and experience achieved in the application of this method can prove to be of major importance in the study of MBSs with elastic elements. The application of the method in the particular framework of MBS presents particularities that must be taken into account in the case of FEM applications. The main step in the case of the MBS study using FEMs is the writing of motion equations, a stage where new terms appear in the evolution equations corresponding to the effects produced by the accelerations of elements in a general rigid motion. These new terms represent the specificity in writing these equations, and the problem that arises in modeling concerns the selection of the most suitable approximations to obtain equations that describe reality as faithfully as possible. Analytical mechanics is involved in writing these equations, which offers alternative methods to obtain these equations [1]. Analytical mechanics is used because its methods offer maximum generality in approaching any MBS problems. Additionally, modern systems must be studied considering a multitude of elastic or rigid elements, characterized by numerous parameters. The studied systems are not at all simple and involve consistent design and numerical simulation efforts. In these circumstances, it is obvious that any advantage that a chosen method can offer, for specific applications, regarding the modeling, number of calculation operations involved, and costs generated by the necessary simulations, is welcome. We would like to be able to use the methods of writing equations of motion allowing us to obtain efficient algorithms and easy-to-use software. Analytical mechanics allows the use of procedures with a high degree of generality. In this way, different applications can be treated in a unified way. Obviously, the number of degrees of freedom of such a system plays a very important role in modeling and simulations, and therefore the existence of procedures that allow the unitary and general treatment of any system, provided by Analytical mechanics, is extremely useful in the study of such problems [2,3,4].

Analytical mechanics also provides several equivalent formulations for the fundamental laws of mechanics. The results obtained using any of these forms are identical and, for this reason, any of these forms can be used in the analysis. The researcher uses one of these methods, depending on his/her experience and the concrete problem he/she has to study. In this way, he/she can use the advantages offered by one method or another.

The standard method used by most researchers in the current practice is Lagrange’s equations. This is due to the fact that this method ensures a sufficient degree of generality for most problems encountered in engineering practice and uses very simple mechanics notions familiar to researchers (such as kinetic energy, potential energy, work, or momentum). Therefore, it can very well describe the constraints to which the system is subjected. In the specialized literature, there are many practical applications where this method is applied; we only mention [5,6,7,8] for clarification.

Analytical mechanics provides freedom in choosing the method for obtaining the equations of motion that, for certain engineering applications and in certain circumstances, can present considerable advantages in dynamic analyses. These advantages have been noticed and used by some of the researchers, thus making modeling easier and reducing the analysis time [9,10,11]. Obviously, the use of equivalent methods also led to the research of the possibility of using the most suitable numerical analysis software, so as to ensure an optimal approach to the analysis. Of course, within these methods, the FEM plays a primary role [12,13,14].

The equivalent formulas offered by analytical mechanics must be assisted in their applications by appropriate numerical methods. Such methods that help reduce the costs related to the numerical analysis are presented in [15]. Symbolic formalism represents an important tool in reducing the time dedicated to modeling such systems [16].

An MBS has a general rigid movement, over which the deformations of elastic bodies are superimposed. An attempt at standardization in order to facilitate the writing of algorithms in MBSs is presented in [17]. Applications of MBS procedures appear in the most diverse fields [18]. Models for solving MBS problems in different fields are presented in [19,20]. The studies for such systems and the analysis of some engineering applications are performed in [21,22].

The practical applications for systems currently used in the industry require the study of complex systems, and the modeling and simulation of such systems involve the mobilization of significant resources and the high costs. Reducing the size of such systems is desirable. The presentation of a way to achieve this outcome is shown in [23], where reduced order models are used. Different techniques for the symbolic writing of motion equations and reducing the time required for modeling have been developed in the last decade. Two methods presenting efficient models for a multibody system with elastic elements are presented in [24,25]. Classical methods applied in the MBS are described in [26,27]. A high-precision formulation for a 3D beam element is presented in [28], and the mathematical methods used to study such problems are studied in [29,30].

The main step in the FEM analysis of an MBS system is the writing of motion equations for a single finite element, taking into account the chosen shape functions and general rigid motion in which the system is trained. This problem has been studies for a long time and various results have been obtained [31]. Classical procedures used to solve these problems are presented in [32,33], using different solutions for different applications [34,35]. The use of a composite material for MBS fabrication is analyzed in [36]. Usually, the application studies beam elements [37,38,39]. A known method of analytical mechanics is used in [40] to model the system.

Two-dimensional elements are frequently used in practical applications, many engineering systems being planar, or moving in a plane or being modeled as planar elements. An example where two-dimensional thin plates are used in engineering applications is presented in [41]. In the mentioned work, the natural frequencies of the system were studied in order to optimize it. The method used for modeling was the transfer matrix method for the multibody system. The theoretical bases for the use of this method were also presented; some concrete examples were presented that justified the presented study. An example of a two-dimensional finite element analysis of the wing skin of an aircraft is presented in [42]. The studied plane had six-degrees of freedom regarding its rigid movement over which vibrations, due to the elasticity of the wings, were superimposed. The effect of the elasticity of the plane’s elements on the overall movement was studied. The study allowed for decisions to be made during the design stage of the plane. A shell element with geometrical features of rotating blades was developed in [43]. The proposed model was validated by comparing its predictions with the problems studied in the reference works. A linear dynamic model using an FEM for large structures was proposed in [44]. The advantages of the presented model were highlighted in an application for the synthesis of the dynamic model of a spacecraft. The ship was equipped with two mobile and flexible solar arrays. The paper analyzed the use of an FEM in the case of the study of an MBS using two-dimensional shell finite elements. An innovative finite element with a thin hyper-elastic shell was presented in [45]. The proposal was based on the Kirchhoff–Love theory. With the help of this finite element, the influence of the use of different constitutive models in obtaining equations for static and dynamic analyses was studied. An interesting absolute nodal coordinate formulation was proposed for the study of belt-drive systems using the FEM [46]. The proposed method allowed a significant reduction in the degrees of freedom of the system. The application of the Gibbs–Appell method to two-dimensional finite elements was presented in [47,48]. Other interesting results can be found in [49,50].

Planar problems or problems that can be modeled in this way are numerous in the engineering field and, as such, a detailed study of this type of problems is required. This paper highlights the particularities of using two-dimensional finite elements in the analysis, methodically analyzing the main methods offered by analytical mechanics for the study of the problem. Obviously, the method of Lagrange’s equations remains the most used method; however, equivalent formulas can prove their usefulness in the case of common applications. The most used methods are Gibbs–Appell (GA) formalism, Maggi’s equations, Kane’s equations, and Hamilton’s equations, which is why they are analyzed in this work.

2. Basic Notion and Notations

2.1. Basic Kinematics and Dynamics Notions and Notations

For the ease of explanation, some notions and notations were introduced, most of them widely used and accepted by researchers. These notations were used uniformly to present alternative methods of the analysis and writing of motion equations [51]. In all these methods, a single finite element was analyzed, relative to a local reference system, which participated in the rigid movement together with the considered finite element, having, at the same time, an elastic deformation. The movement of the entire system was related to a global reference system. The speed and acceleration of the local reference system O, ${\overline{v}}_o(v_{O1}i+v_{O2}j)$ and ${\overline{a}}_o(a_{O1}i+a_{O2}j)$, respectively, were considered as known. We examined a planar motion, with the elastic deformations in the same plane. As a result, the mobile reference system had angular velocity and angular acceleration values, which were denoted by $\overline{\omega }=\omega \overline{k}$ and $\overline{\epsilon }=\epsilon \overline{k}$, respectively.

The indexes L (from local) and G (from global) indicate whether a quantity is expressed relative to the local or global reference systems, respectively. In the case of the planar motion, angle θ of the rotation of the local reference system in relation to the global fixed reference system is defined by rotation matrix R, which in our case is:

$$\left[R\right]=\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right].$$

(1)

An arbitrary vector is transformed according to the following relations:

$$a_{1,G}=a_{1,L}\cos \theta -a_{2,L}\hspace{0.17em}\sin \theta \hspace{1em};\hspace{1em}{a}_{2,G}=a_{1,L}\sin \theta +a_{2,L}\hspace{0.17em}\cos \theta .$$

(2)

Then, we used ${\overline{r}}_O$ to denote the position vector of the origin of the mobile reference frame related to the global fixed reference system with the origin in O₁. Its components in the global reference frame were $(X_{O,1},X_{O,2})$ and $(x_{O,1},x_{O,2})$ in the local reference frame. Similarly ${\overline{r}}_M$, the position vector of the point M before deformation, related to the global fixed reference system, had the components $(X_{M,1},X_{M,2})$ in the global reference frame and $(x_{M,1},x_{M,2})$ in the local reference frame, and ${\overline{r}}_{M\prime }$ represented the position vector of point M after deformation, when transformed into M’, with the components $(X_{M\prime ,1},X_{M\prime ,2})$, respectively $(x_{M\prime ,1},x_{M\prime ,2})$. The position vector of point M with respect to origin O was $\overline{r}$, with the components $(X_1,X_2)$, respectively $(x_1,x_2)$, in the two reference frames, global and local, and the displacement vector of point M when became M’, $\overline{u}$, included the components $(u,v)$.

An arbitrary point M of the finite element became, following deformation, point M’. Its coordinates are (in the global reference frame):

$$\begin{array}{l} \\ X_{M\prime ,1}=X_{O,1}+(x_1+u)\cos \theta -(x_2+v)\sin \theta \hspace{1em};\hspace{1em}\\ X_{M\prime ,2}=X_{O,2}+(x_1+u)\sin \theta +(x_2+v)\cos \theta \hspace{0.17em}.\end{array}.$$

(3)

The relation between independent generalized coordinates and the displacement of a current point was approximated in the FEM through the shape functions $N_{1j}$, $N_{2j}\hspace{0.17em},\hspace{1em}j=1,p$, by the following relation:

$$\begin{array}{l}u=N_{1j}{\delta }_j\hspace{1em};\hspace{1em}j=\overline{1,p}\\ v=N_{2j}{\delta }_j\hspace{1em};\hspace{1em}j=\overline{1,p}\end{array}$$

(4)

where:

$${\left\{\delta \right\}}_L=\left\{\begin{array}{c}{\delta }_1\\ {\delta }_2\\ ⋮\\ \\ {\delta }_p\end{array}\right\}$$

(5)

is the vector of the independent generalized coordinates. $\delta$₁, $\delta$₂, …., δ_p are the independent coordinates and p is the number of these coordinates. The displacements are approximated using different shape functions defined by the type of finite element used. Considering Equation (4), Equation (3) becomes:

$$\begin{array}{l} \\ X_{1,M\prime }=X_{1,O}+\left(x_1+N_{1j}{\delta }_j\right)\cos \theta -\left(x_2+N_{2j}{\delta }_j\right)\sin \theta \hspace{1em};\hspace{1em}\\ X_{2,M\prime }=X_{2,O}+\left(x_1+N_{1j}{\delta }_j\right)\sin \theta +\left(x_2+N_{2j}{\delta }_j\right)\cos \theta \hspace{1em};\hspace{1em}j=\overline{1,p}.\end{array}.$$

(6)

By differentiation, the components of the velocity vector of M’ can be obtained:

$$\begin{array}{l} \\ {\dot{X}}_{M\prime ,1}={\dot{X}}_{O,1}-\omega \left(x_1+N_{1j}{\delta }_j\right)\sin \theta +\left(x_1+N_{1j}{\dot{\delta }}_j\right)\cos \theta -\omega \left(x_2+N_{2j}{\delta }_j\right)\cos \theta -\left(x_2+N_{2j}{\dot{\delta }}_j\right)\sin \theta \hspace{1em};\hspace{1em}\\ {\dot{X}}_{M\prime ,2}={\dot{X}}_{O,2}+\omega \left(x_1+N_{1j}{\delta }_j\right)\cos \theta +\left(x_1+N_{1j}{\dot{\delta }}_j\right)\sin \theta -\omega \left(x_2+N_{2j}{\delta }_j\right)\sin \theta +\left(x_2+N_{2j}{\dot{\delta }}_j\right)\cos \theta \hspace{1em};\hspace{1em}j=\overline{1,p}.\end{array}.$$

(7)

and the acceleration:

$$\begin{array}{l} \\ {\ddot{X}}_{M\prime ,1}={\ddot{X}}_{O,1}-\epsilon \left(x_1+N_{1j}{\delta }_j\right)\sin \theta -{\omega }^2\left(x_1+N_{1j}{\delta }_j\right)\cos \theta -2\omega \left(x_1+N_{1j}{\dot{\delta }}_j\right)\sin \theta +\left(x_1+N_{1j}{\ddot{\delta }}_j\right)\cos \theta -\\ -\epsilon \left(x_2+N_{2j}{\delta }_j\right)\cos \theta +{\omega }^2\left(x_2+N_{2j}{\delta }_j\right)\sin \theta -2\omega \left(x_2+N_{2j}{\dot{\delta }}_j\right)\cos \theta -\left(x_2+N_{2j}{\ddot{\dot{\delta }}}_j\right)\sin \theta \hspace{1em};\hspace{1em}\\ {\ddot{X}}_{M\prime ,2}={\ddot{Y}}_{O,2}+\epsilon \left(x_1+N_{1j}{\delta }_j\right)\cos \theta -{\omega }^2\left(x_1+N_{1j}{\delta }_j\right)\sin \theta +2\omega \left(x_1+N_{1j}{\dot{\delta }}_j\right)\cos \theta +\left(x_1+N_{1j}{\ddot{\delta }}_j\right)\sin \theta -\\ -\epsilon \left(x_2+N_{2j}{\delta }_j\right)\sin \theta -{\omega }^2\left(x_2+N_{2j}{\delta }_j\right)\cos \theta -2\omega \left(x_2+N_{2j}{\dot{\delta }}_j\right)\sin \theta +\left(x_2+N_{2j}{\ddot{\delta }}_j\right)\cos \theta \hspace{1em};\hspace{1em}j=\overline{1,p}.\end{array}$$

(8)

In the local coordinate system, we have, successively:

$$\begin{array}{l} \\ x_{M\prime ,1}=x_{O,1}+x_1+N_{1j}{\delta }_j\hspace{1em};\hspace{1em}\\ x_{M\prime ,2}=x_{O,2}+x_2+N_{2j}{\delta }_j\hspace{1em};\hspace{1em}j=\overline{1,p}.\end{array}$$

(9)

$$\begin{array}{l}{\dot{x}}_{M\prime ,1}={\dot{x}}_{O,1}-\omega x_2-\omega N_{2r}{\delta }_r+N_{1r}{\dot{\delta }}_r\hspace{1em};\hspace{1em}\\ {\dot{x}}_{M\prime ,2}={\dot{x}}_{O,2}+\omega x_1+\omega N_{1r}{\delta }_r+N_{2r}{\dot{\delta }}_r\hspace{1em};\hspace{1em}r=\overline{1,p}.\end{array}.$$

(10)

2.2. Kinetic Energy

For a single finite element, the kinetic energy is:

$$\begin{gathered} E_C=\frac{1}{2}{\displaystyle \underset{V}{\int }\rho \left[{\left({\dot{x}}_{M\prime ,1}\right)}^2+{\left({\dot{x}}_{M\prime ,2}\right)}^2\right]}dV= \\ \frac{1}{2}{\displaystyle \underset{V}{\int }\rho \left[{\left({\dot{x}}_{O,1}-\omega x_2-\omega N_{2r}{\delta }_r+N_{1r}{\dot{\delta }}_r\right)}^2+{\left({\dot{x}}_{O,1}+\omega x_1+\omega N_{1r}{\delta }_r+N_{2r}{\dot{\delta }}_r\right)}^2\right]}dV \\ \begin{array}{l}=\frac{1}{2}\left({\dot{x}}_{O,1}^2+{\dot{x}}_{O,2}^2\right){\displaystyle \underset{V}{\int }\rho }dV+\frac{1}{2}{\omega }^2{\displaystyle \underset{V}{\int }\rho \left(x_1^2+x_2^2\right)}dV+\frac{1}{2}{\omega }^2{\delta }_t{\delta }_r{\displaystyle \underset{V}{\int }\rho \left( N_{2r}{N}_{2t}+N_{1r}{N}_{1t}\right)}dV+\frac{1}{2}{\dot{\delta }}_r{\dot{\delta }}_t{\displaystyle \underset{V}{\int }\rho \left(N_{1r}{N}_{1t}+N_{2r}{N}_{2t}\right)}dV+\\ -\omega \left({\dot{x}}_{O,1}{\displaystyle \underset{V}{\int }\rho x_1} dV-\omega {\dot{x}}_{O,2}{\displaystyle \underset{V}{\int }\rho x_2 }dV\right)-\omega {\delta }_r\left({\dot{x}}_{O,1}{\displaystyle \underset{V}{\int }\rho N_{2r}}dV+{\dot{x}}_{O,2}{\displaystyle \underset{V}{\int }\rho N_{1r}}dV\right)+\left({\dot{x}}_{O,1}{\displaystyle \underset{V}{\int }\rho N_{1r}}dV+y{\dot{x}}_{O,2} {\dot{\delta }}_r{\displaystyle \underset{V}{\int }\rho N_{2r}}dV\right){\dot{\delta }}_r+\\ +{\omega }^2{\delta }_r\left({\displaystyle \underset{V}{\int }\rho N_{2r}{x}_2 }dV+{\displaystyle \underset{V}{\int }\rho N_{1r}{x}_1 }dV\right)-\omega {\dot{\delta }}_r\left({\displaystyle \underset{V}{\int }\rho x_2{N}_{1r}}dV-{\displaystyle \underset{V}{\int }\rho x_1{N}_{2r}}dV\right)-\omega {\delta }_r{\dot{\delta }}_t\left({\displaystyle \underset{V}{\int }\rho N_{2r}{N}_{1t}}dV-{\displaystyle \underset{V}{\int }\rho N_{1r}{N}_{2t}}dV\right)\end{array} \end{gathered}$$

(11)

With the notations:

$$m={\displaystyle \underset{V}{\int }\rho }dV\hspace{1em};\hspace{1em}{J}_O={\displaystyle \underset{V}{\int }\rho \left(x_1^2+x_2^2\right)}dV\hspace{1em};\hspace{1em}{m}_{rt}={\displaystyle \underset{V}{\int }\rho }{N}_{kr}{N}_{kt}dV,$$

(12)

$$\begin{array}{l}{m}_{ij,rt}={\displaystyle \underset{V}{\int }\rho N_{ir}{N}_{jt}}dV\hspace{1em};\hspace{1em}{S}_1={\displaystyle \underset{V}{\int }\rho }{x}_{1}dV\hspace{1em};\hspace{1em}{S}_2={\displaystyle \underset{V}{\int }\rho }{x}_{2}dV\hspace{1em};\hspace{1em}\\ m_{O,kr}^I={\displaystyle \underset{V}{\int }\rho N_{kr}}dV\hspace{1em};\hspace{1em}{m}_{1,mr}={\displaystyle \underset{V}{\int }\rho \hspace{0.17em}{x}_1{N}_{mr}}dV\hspace{1em};\hspace{1em}{m}_{2,mr}={\displaystyle \underset{V}{\int }\rho \hspace{0.17em}{x}_2{N}_{mr}}dV\hspace{1em};\hspace{1em}\end{array},$$

(13)

It exists the relation:

$${\alpha }_{ij}{\alpha }_{jk}={\delta }_{ik}$$

(14)

where ${\delta }_{ik}$ is the Kronecker delta.

The kinetic energy has the following expression:

$$\begin{array}{l}{E}_c=\frac{1}{2}m\left({\dot{x}}_{O,1}^2+{\dot{x}}_{O,1}^2\right)+\frac{1}{2}{\omega }^2{J}_O+\frac{1}{2}{\omega }^2{\delta }_t{\delta }_r{m}_{rt}+\frac{1}{2}{\dot{\delta }}_r{\dot{\delta }}_t{m}_{rt}^{}+\\ -\omega \left({\dot{x}}_{O,1}{S}_2- {\dot{x}}_{O,2}{S}_2\right)-\omega {\delta }_r\left({\dot{x}}_{O,1}{m}_{O,2r}^I+{\dot{x}}_{O,2}{m}_{O,1r}^I\right)+\left({\dot{x}}_{O,1}{m}_{O,1r}^I+{\dot{x}}_{O,2} {\dot{\delta }}_r{m}_{O,2r}^I\right){\dot{\delta }}_r+\\ +{\omega }^2{\delta }_r\left(m_{2,2r}+m_{1,1r}\right)-\omega {\dot{\delta }}_r\left(m_{2,1r}-m_{1,2r}\right)-\omega {\delta }_r{\dot{\delta }}_t\left(m_{12,rt}-m_{21,rt}\right)\end{array},$$

(15)

2.3. Potential Energy

The potential energy of the finite element is:

$$E_p=\frac{1}{2}{\displaystyle {\int }_V\left({\sigma }_{11}{\epsilon }_{11}+2{\sigma }_{12}{\epsilon }_{12}+{\sigma }_{22}{\epsilon }_{22}\right)}dV,$$

(16)

In this relation, ${\epsilon }_{ij}$ is the plane strain tensor and ${\sigma }_{ij}$ the plane stress tensor. The generalized Hooke law has the well -known formulation:

$${\sigma }_{ij}=H_{ijkl}{\epsilon }_{kl}$$

(17)

In the case of the plane strain’s state for an isotropic material, Equation (17) becomes:

$${\sigma }_{11}=\frac{E}{1-2\mu }{\epsilon }_{11}\hspace{1em};\hspace{1em}{\sigma }_{22}=\frac{E}{1-2\mu }{\epsilon }_{22}\hspace{1em};\hspace{1em}{\sigma }_{12}=\frac{E}{2\left(1-\mu \right)}{\epsilon }_{12}$$

(18)

The strains in the plane state of deformation are:

$${\epsilon }_{11}=\frac{\partial u}{\partial x_1}\hspace{1em};\hspace{1em}{\epsilon }_{22}=\frac{\partial v}{\partial x_2}\hspace{1em};\hspace{1em}{\epsilon }_{12}=\frac{1}{2}\left(\frac{\partial u}{\partial x_2}+\frac{\partial v}{\partial x_1}\right)$$

(19)

$${\epsilon }_{11}=\frac{\partial u}{\partial x_1}=\frac{\partial N_{1r}}{\partial x_1}{\delta }_r\hspace{1em};\hspace{1em}{\epsilon }_{22}=\frac{\partial v}{\partial x_2}=\frac{\partial N_{2r}}{\partial x_2}{\delta }_r\hspace{1em};\hspace{1em}{\epsilon }_{12}=\frac{1}{2}\left(\frac{\partial N_{1r}}{\partial x_2}+\frac{\partial N_{2r}}{\partial x_1}\right){\delta }_r$$

(20)

$${\sigma }_{11}=\frac{E}{1-2\mu }\frac{\partial N_{1r}}{\partial x_1}{\delta }_r\hspace{1em};\hspace{1em}{\sigma }_{22}=\frac{E}{1-2\mu }\frac{\partial N_{2r}}{\partial x_2}{\delta }_r\hspace{1em};\hspace{1em}{\sigma }_{12}=\frac{E}{4\left(1-\mu \right)}\left(\frac{\partial N_{1r}}{\partial x_2}+\frac{\partial N_{2r}}{\partial x_1}\right){\delta }_r$$

(21)

Using Equations (17) and (18), we obtained:

$$E_p=\frac{1}{2}{\displaystyle {\int }_V{\sigma }_{ij}{\epsilon }_{ij}}dV=\frac{1}{2}{\displaystyle {\int }_V\left[\frac{E}{1-2\mu }\frac{\partial N_{1r}}{\partial x_1}\frac{\partial N_{1t}}{\partial x_1}+\frac{E}{1-2\mu }\frac{\partial N_{2r}}{\partial x_2}\frac{\partial N_{2t}}{\partial x_2}+\frac{E}{8\left(1-\mu \right)}\left(\frac{\partial N_{1r}}{\partial x_2}+\frac{\partial N_{2r}}{\partial x_1}\right)\left(\frac{\partial N_{1t}}{\partial x_2}+\frac{\partial N_{2t}}{\partial x_1}\right)\right]dV}{\delta }_r{\delta }_t=\frac{1}{2}{m}_{rt}{\delta }_r{\delta }_t$$

(22)

where the notation:

$$m_{rt}=\frac{1}{2}{\displaystyle {\int }_V\left[\frac{E}{1-2\mu }\frac{\partial N_{1r}}{\partial x_1}\frac{\partial N_{1t}}{\partial x_1}+\frac{E}{1-2\mu }\frac{\partial N_{2r}}{\partial x_2}\frac{\partial N_{2t}}{\partial x_2}+\frac{E}{8\left(1-\mu \right)}\left(\frac{\partial N_{1r}}{\partial x_2}+\frac{\partial N_{2r}}{\partial x_1}\right)\left(\frac{\partial N_{1t}}{\partial x_2}+\frac{\partial N_{2t}}{\partial x_1}\right)\right]}.$$

(23)

is performed.

2.4. Work

Considering the generalized concentrated forces $q_i$ $i=\overline{1,p}$, they produce the following equation:

$$W^c=q_i{\delta }_i\hspace{1em};\hspace{1em}i=\overline{1,p}.$$

(24)

The generalized volume forces $q_i^{\ast }$$i=\overline{1,p}$, similarly present the following equation:

$$W^{\mathrm{d}}=q_i^{\ast }{\delta }_i\hspace{1em};\hspace{1em}i=\overline{1,p},$$

(25)

Therefore, the total work becomes:

$$W=\left(W^c+W^{\mathrm{d}}\right)=\left(q_i+q_i^{\ast }\right){\delta }_i\hspace{1em};\hspace{1em}i=\overline{1,p},$$

(26)

2.5. Lagrangian

The classical form of the Lagrangian formula for one finite element is [51]:

$$L=E_c-E_p+W,$$

(27)

Using Equations (15) and (21)–(23), these expressions are:

$$\begin{array}{l}L=\frac{1}{2}m\left({\dot{x}}_{O,1}^2+{\dot{x}}_{O,2}^2\right)+\frac{1}{2}{\omega }^2{J}_O+\frac{1}{2}{\omega }^2{\delta }_t{\delta }_r{m}_{rt}+\frac{1}{2}{\dot{\delta }}_r{\dot{\delta }}_t{m}_{rt}^{}+\\ -\omega \left({\dot{x}}_{O,1}{S}_2- {\dot{x}}_{O,2}{S}_1\right)-\omega {\delta }_r\left({\dot{x}}_O{m}_{O,2r}^I+{\dot{x}}_{O,2}{m}_{O,1r}^I\right)+\left({\dot{x}}_{O,1}{m}_{O,1r}^I+{\dot{x}}_{O,2} {\dot{\delta }}_r{m}_{O,2r}^I\right){\dot{\delta }}_r+\\ +{\omega }^2{\delta }_r\left(m_{2,2r}+m_{1,1r}\right)-\omega {\dot{\delta }}_r\left(m_{2,1r}-m_{1,2r}\right)-\omega {\delta }_r{\dot{\delta }}_t\left(m_{12,rt}-m_{21,rt}\right)\\ -k_{rt}{\delta }_r{\delta }_t+q_r{\delta }_r+q_r^{\ast }{\delta }_r\hspace{1em};\hspace{1em}r,t=\overline{1,p}\end{array},$$

2.6. Momentum

The momentum vector has following components:

$$p_{r,L}=\frac{\partial L}{\partial {\dot{\delta }}_r}.$$

(28)

It obtains:

$$p_{r,L}=\frac{\partial L}{\partial {\dot{\delta }}_r}=m_{rt}{\dot{\delta }}_t+\left({\dot{x}}_O{m}_{O,1r}^I+{\dot{y}}_O m_{O,2r}^I\right)-\omega \left(m_{y,1r}-m_{x,2r}\right)-\omega \left(m_{xy,rt}-m_{yx,rt}\right){\delta }_r\hspace{1em};\hspace{1em}r,t=\overline{1,p}.$$

(29)

The inverse matrix $m_{ur}^{\ast }$ is chosen so that:

$$m_{ur}^{\ast }{m}_{rt}={\delta }_{ut}\hspace{1em};\hspace{1em}u,r,t=\overline{1,p}$$

(30)

where ${\delta }_{ut}$ is the Kronecker delta, and premultiplying Equation (28) with $m_{ur}^{\ast }$ produces:

$$\begin{array}{l}{\dot{\delta }}_u=m_{ur}^{\ast }{p}_{r,L}-m_{ur}^{\ast }\left({\dot{x}}_O{m}_{O,1r}^I+{\dot{y}}_O m_{O,2r}^I\right)+\omega m_{ur}^{\ast }\left(m_{y,1r}-m_{x,2r}\right)+\omega m_{ur}^{\ast }\left(m_{xy,rt}-m_{yx,rt}\right){\delta }_r\hspace{1em};\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}u,r,t=\overline{1,p}\end{array}$$

(31)

2.7. Hamiltonian

The classical Hamiltonian form is the following:

$$H=\frac{\partial L}{\partial {\dot{\delta }}_r}{\dot{\delta }}_r-L.$$

(32)

In Equation (32), we integrated Equation (27) for the Lagrangian calculation.

2.8. Energy of Accelerations

An important method offered by analytical mechanics as an alternative form to Lagrange’s equations is the Gibbs–Appell equation. This method also implies the energy of acceleration. This is defined as:

$$E_a=\frac{1}{2}{\displaystyle \underset{V}{\int }\rho a^{2}dV}.$$

(33)

Thus, differentiating Equation (7):

$$\begin{array}{l} \\ {\ddot{X}}_{M\prime ,1}={\ddot{X}}_{O,1}-\epsilon \left(x_1+N_{1j}{\delta }_j\right)\sin \theta -{\omega }^2\left(x_1+N_{1j}{\delta }_j\right)\cos \theta -2\omega \left(x_1+N_{1j}{\dot{\delta }}_j\right)\sin \theta +\left(x_1+N_{1j}{\ddot{\delta }}_j\right)\cos \theta -\\ \epsilon \left(x_2+N_{2j}{\delta }_j\right)\cos \theta +{\omega }^2\left(x_2+N_{2j}{\delta }_j\right)\sin \theta -2\omega \left(x_2+N_{2j}{\dot{\delta }}_j\right)\cos \theta -\left(x_2+N_{2j}{\ddot{\dot{\delta }}}_j\right)\sin \theta \hspace{1em};\hspace{1em}\\ {\ddot{X}}_{M\prime ,2}={\ddot{X}}_{O,2}+\epsilon \left(x_1+N_{1j}{\delta }_j\right)\cos \theta -{\omega }^2\left(x_1+N_{1j}{\delta }_j\right)\sin \theta +2\omega \left(x_1+N_{1j}{\dot{\delta }}_j\right)\cos \theta +\left(x_1+N_{1j}{\ddot{\delta }}_j\right)\sin \theta \\ -\epsilon \left(x_2+N_{2j}{\delta }_j\right)\sin \theta -{\omega }^2\left(x_2+N_{2j}{\delta }_j\right)\cos \theta -2\omega \left(x_2+N_{2j}{\dot{\delta }}_j\right)\sin \theta +\left(x_2+N_{2j}{\ddot{\delta }}_j\right)\cos \theta \hspace{1em};\hspace{1em}j=\overline{1,p}.\end{array}.$$

(34)

In the local coordinate system:

$$\begin{gathered} \begin{array}{l}{\ddot{x}}_{M\prime ,1}={\ddot{X}}_{M\prime ,1}\cos \theta +{\ddot{X}}_{M\prime ,2}\sin \theta ={\ddot{X}}_{O,1}\cos \theta +{\ddot{X}}_{O,2}\sin \theta -{\omega }^2\left(x_1+N_{1j}{\delta }_j\right)+\\ \left(x_1+N_{1j}{\ddot{\delta }}_j\right)\cos \theta -\epsilon \left(x_2+N_{2j}{\delta }_j\right)-2\omega \left(x_2+N_{2j}{\dot{\delta }}_j\right)=\\ {\ddot{x}}_{O,1}+\left(x_1+N_{1j}{\ddot{\delta }}_j\right)-{\omega }^2\left(x_1+N_{1j}{\delta }_j\right)-\epsilon \left(x_2+N_{2j}{\delta }_j\right)-2\omega \left(x_2+N_{2j}{\dot{\delta }}_j\right).\end{array} \\ \begin{array}{l}{\ddot{x}}_{M\prime ,2}=-{\ddot{X}}_{M\prime ,1}\sin \theta +{\ddot{X}}_{M\prime ,2}\cos \theta =-{\ddot{X}}_{O,1}\sin \theta +{\ddot{X}}_{O,2}\cos \theta +\epsilon \left(x_1+N_{1j}{\delta }_j\right)\sin \theta +\\ 2\omega \left(x_1+N_{1j}{\dot{\delta }}_j\right)\sin \theta -{\omega }^2\left(x_2+N_{2j}{\delta }_j\right)+\left(x_2+N_{2j}{\ddot{\delta }}_j\right)=\\ {\ddot{x}}_{O,2}+\left(x_2+N_{2j}{\ddot{\delta }}_j\right)-{\omega }^2\left(x_2+N_{2j}{\delta }_j\right)+\epsilon \left(x_1+N_{1j}{\delta }_j\right)+2\omega \left(x_1+N_{1j}{\dot{\delta }}_j\right)\end{array}. \end{gathered}$$

(35)

and the acceleration can be obtained via the following relation:

$$a^2={\ddot{x}}_{M\prime ,1}+{\ddot{x}}_{M\prime ,2}.$$

(36)

3. Evolution Equations for the Finite Element Method

3.1. Lagrange’s Equations

The use of Lagrange’s equations proved to be useful in the case of solving difficult problems, such as those involved in an MBS. An advantage was the use of scalar quantities instead of vector ones. To use Lagrange’s equations, the expressions of kinetic energy and potential energy were used, and the determination of the generalized forces was performed. Although the method is very old, it has not yet lost its importance. The known form of Lagrange’s equation is:

$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{\delta }}_i}\right)-\frac{\partial L}{\partial {\delta }_i}=0\hspace{1em};\hspace{1em}i=\overline{1,p}.$$

(37)

Applying the procedure involved by the Lagrange method, we achieve the following results:

$$\begin{gathered} \frac{\partial L}{\partial {\dot{\delta }}_r}={\dot{\delta }}_t{m}_{rt}^{}+\left({\dot{x}}_{O,1}{m}_{O,1r}^I+{\dot{x}}_{O,2} m_{O,2r}^I\right)-\omega \left(m_{2,1r}-m_{1,2r}\right)-\omega {\delta }_r\left(m_{12,rt}-m_{21,rt}\right)\hspace{1em};\hspace{1em}r,t=\overline{1,p} \\ \frac{d}{dt}\frac{\partial L}{\partial {\dot{\delta }}_r}={\ddot{\delta }}_t{m}_{rt}^{}+\left({\ddot{x}}_{O,1}{m}_{O,1r}^I+{\ddot{x}}_{O,2} m_{O,2r}^I\right)-\epsilon \left(m_{2,1r}-m_{1,2r}\right)-\epsilon {\delta }_r\left(m_{12,rt}-m_{21,rt}\right)-\omega {\dot{\delta }}_r\left(m_{12,rt}-m_{21,rt}\right)\hspace{1em};\hspace{1em}r,t=\overline{1,p} \end{gathered}$$

(38)

$$\begin{gathered} \begin{array}{l}\frac{\partial L}{\partial {\delta }_r}={\omega }^2{\delta }_t{m}_{rt}-\omega \left({\dot{x}}_{O,1}{m}_{O,2r}^I+{\dot{x}}_{O,2}{m}_{O,1r}^I\right)+{\omega }^2\left(m_{2,2r}+m_{1,1r}\right)-\omega {\dot{\delta }}_t\left(m_{12,rt}-m_{21,rt}\right)\\ -k_{rt}{\delta }_t+q_r+q_r^{\ast }\hspace{1em};\hspace{1em}r,t=\overline{1,p}\end{array} \\ \frac{d}{dt}\frac{\partial L}{\partial {\dot{\delta }}_r}-\frac{\partial L}{\partial {\delta }_r}= \\ \begin{array}{l}{\ddot{\delta }}_t{m}_{rt}^{}+2\omega {\dot{\delta }}_r\left(m_{12,rt}-m_{21,rt}\right)+\left[k_{rt}-\epsilon {\delta }_r\left(m_{12,rt}-m_{21,rt}\right)-{\omega }^2{m}_{rt}\right]{\delta }_t+\left({\ddot{x}}_{O,1}{m}_{O,1r}^I+{\ddot{x}}_{O,2} m_{O,2r}^I\right)-\epsilon \left(m_{2,1r}-m_{1,2r}\right)\hspace{1em}\\ +\omega \left({\dot{x}}_{O,1}{m}_{O,2r}^I+{\dot{x}}_{O,2}{m}_{O,1r}^I\right)-{\omega }^2\left(m_{2,2r}+m_{1,1r}\right)-q_r{\delta }_r-q_r^{\ast }{\delta }_r=0\hspace{1em};\hspace{1em}r,t=\overline{1,p}\end{array} \end{gathered}$$

(39)

The evolution equation for the finite element is:

$$\begin{array}{l}{m}_{rt}{\ddot{\delta }}_t+2\omega {\dot{\delta }}_r\left(m_{12,rt}-m_{21,rt}\right)+\left[k_{rt}-\epsilon \left(m_{12,rt}-m_{21,rt}\right)-{\omega }^2{\delta }_t{m}_{rt}\right]{\delta }_t=-\left({\ddot{x}}_{O,1}{m}_{O,1r}^I+{\ddot{x}}_{O,2} m_{O,2r}^I\right)+\epsilon \left(m_{2,1r}-m_{1,2r}\right)\hspace{1em}\\ -\omega \left({\dot{x}}_{O,1}{m}_{O,2r}^I+{\dot{x}}_{O,2}{m}_{O,1r}^I\right)+{\omega }^2\left(m_{2,2r}+m_{1,1r}\right)-q_r-q_r^{\ast }=0\hspace{1em};\hspace{1em}r,t=\overline{1,p}\end{array}.$$

(40)

3.2. Hamilton’s Method

With the help of Lagrange’s equations, the motion equations of the MBS were determined, obtaining a system of second-order differential equations. In order to be able to solve this system, it must be transformed into a system of differential equations of the first order by introducing new unknowns, but of a double dimension. In the case of using the Hamilton formalism, the equations of motion were directly obtained as a system of differential equations of the first order, of a double dimension. In this case, the unknowns are the generalized coordinates and the canonical conjugate momenta:

$$p_{i,L}=\frac{\partial L}{\partial {\delta }_i}\hspace{1em};\hspace{1em}i=\overline{1,p}.$$

(41)

Therefore, the main difference between Lagrange’s and Hamilton’s methods is the use of canonical conjugated momenta instead of generalized velocities. The main advantage of applying these equations was to directly obtain a system of first-order equations. These could be solved directly; the usual commercial software included moduli dedicated to solving first-order differential equations. To solve a second-order system, it is necessary to introduce new functions in order to reduce the order of the system [52,53,54,55,56,57].

The classical and well-known Hamilton’s equations are:

$${\dot{\delta }}_r=\frac{\partial H}{\partial p_{r,L}}\hspace{1em};\hspace{1em}{\dot{p}}_{r,L}=-\frac{\partial H}{{\delta }_r}.$$

(42)

Considering Equations (41) and (42), we achieved:

$$\begin{gathered} {\dot{\delta }}_r=m_{ru}^{\ast }{p}_{u,L}-m_{ru}^{\ast }{m}_{O,iu}{\dot{x}}_{O,i}-m_{ru}^{\ast }{m}_{k,iu}{\alpha }_{ij}{\dot{\alpha }}_{jk}-m_{ru}^{\ast }{m}_{ik,ut}{\alpha }_{ij}{\dot{\alpha }}_{jk}{\delta }_t\hspace{0.17em}\hspace{1em};\hspace{1em}u,r,t=\overline{1,p}; \\ \begin{array}{l}{\dot{p}}_{r,L}=\frac{\partial L}{\partial {\delta }_r}={\omega }^2{\delta }_t{m}_{rt}-\omega \left({\dot{x}}_{O,2}{m}_{O,2r}^I+{\dot{x}}_{O,2}{m}_{O,1r}^I\right)+{\omega }^2\left(m_{2,2r}+m_{1,1r}\right)-\omega {\dot{\delta }}_t\left(m_{12,rt}-m_{21,rt}\right)\\ -k_{rt}{\delta }_t+q_r+q_r^{\ast }\hspace{1em};\hspace{1em}r,t=\overline{1,p}\end{array} \end{gathered}$$

(43)

Equation (43) represents a Hamilton equation, a system of 2p first-order equations.

3.3. Gibbs–Appell Equation

The formalism Gibbs–Appell equation is reconsidered at present in the context of its necessity to analyze an MBS with elastic elements. The main notion used in Gibbs–Appell equations is the acceleration energy, obtained in Equation (33). The form of the Gibbs–Appell equation is [58,59,60,61,62]:

$$\frac{\partial E_a}{\partial {\ddot{\delta }}_r}=Q_r\hspace{1em}r=\overline{1,p}.$$

(44)

where $Q_r=q_r+q_r^{\ast }$ (see Section 2.4). After performing the calculations in Equation (33), we identified the following three terms:

• $E_{a2}$ is the part of the acceleration energy containing quadratic values:

$$E_{a2}=\frac{1}{2}{m}_{rt}{\ddot{\delta }}_r{\ddot{\delta }}_t\hspace{1em}r,t=\overline{1,p};$$

(45)

• $E_{a1}$ represents the part of the accelerations energy containing linear values:

$$E_{a1}={\ddot{x}}_{O,i}{\ddot{\delta }}_r{m}_{O,ir}^I+\left({\dot{\alpha }}_{ji}{\dot{\alpha }}_{jk}+{\alpha }_{ij}{\ddot{\alpha }}_{jk}\right){\ddot{\delta }}_r{m}_{k,mr}+\left({\dot{\alpha }}_{ji}{\dot{\alpha }}_{jk}+{\alpha }_{ij}{\ddot{\alpha }}_{jk}\right){\delta }_r{\ddot{\delta }}_t{m}_{kr,mt}+2{\alpha }_{ji}{\dot{\alpha }}_{jk}{\dot{\delta }}_r{\ddot{\delta }}_t{m}_{kr,mt}dV$$

(46)

• $E_{a0}$ does not have any terms containing generalized accelerations. This part played no role in obtaining the Gibbs–Appell equation and was not of any interest to us.

With the three parts, the energy of acceleration is the following:

$$E_a=E_{ao}+E_{a1}+E_{a2};$$

(47)

Equation (44) can be written as:

$$\frac{\partial \left(E_{a1}+E_{a2}\right)}{\partial {\ddot{\delta }}_r}=Q_r\hspace{1em}r=\overline{1,p};$$

(48)

The generalized force vector in our description is:

$$Q_{r,L}=k_{rt}{\delta }_t+q_r+q_r^{\ast }\hspace{1em};\hspace{1em}r,t=\overline{1,p};$$

(49)

We obtained:

$$\frac{\partial E_{a2}}{\partial {\ddot{\delta }}_r}=m_{rt}{\ddot{\delta }}_t\hspace{1em}r,t=\overline{1,p};$$

(50)

$$\begin{array}{l}\frac{\partial E_{a1}}{\partial {\ddot{\delta }}_r}=2\omega {\dot{\delta }}_r\left(m_{12,rt}-m_{21,rt}\right)+\left[-\epsilon {\delta }_r\left(m_{12,rt}-m_{21,rt}\right)-{\omega }^2{m}_{rt}\right]{\delta }_t+\\ \left({\ddot{x}}_{O,1}{m}_{O,1r}^I+{\ddot{x}}_{O,2} m_{O,2r}^I\right)-\epsilon \left(m_{2,1r}-m_{1,2r}\right)+\omega \left({\dot{x}}_{O,1}{m}_{O,2r}^I+{\dot{x}}_{O,2}{m}_{O,1r}^I\right)-\\ {\omega }^2\left(m_{2,2r}+m_{1,1r}\right)\hspace{1em};\hspace{1em}r,t=\overline{1,p}\end{array}$$

(51)

It is obvious that:

$$\frac{\partial E_{a0}}{\partial {\delta }_r}=0\hspace{1em};\hspace{1em}r=\overline{1,p}.$$

(52)

Finally, the result is Equation (40).

Since this method required a smaller number of differentials, as in the case of Lagrange’s equations, the number of calculations decreased. Therefore,, the time and cost of the modeling and simulation processes decreased.

3.4. Maggi’s Equation

The role of Maggi’s equation and its application to related problems are presented in [63]. The generalized independent coordinates of the system are q₁, q₂,…., q_N; between them, there are m linear relationships:

$$a_{ij}(q_1,q_2,\dots \hspace{1em},q_N,t)\hspace{0.17em}{\dot{q}}_j+b_i(q_1,q_2,\dots \hspace{1em},q_N,t)=0\hspace{1em}\hspace{1em},\hspace{1em}i=\overline{1,m}\hspace{0.17em};\hspace{1em}j=\overline{1,N}.$$

(53)

The classic form of Maggi’s equation is:

$$a_{kj}\left[\left(\frac{d}{dt}\left(\frac{\partial E_c}{\partial {\dot{q}}_k}\right)-\frac{\partial E_c}{\partial q_k}\right)-Q_k\right] =0\hspace{1em};\hspace{1em}j=\overline{1,N-m}\hspace{0.17em};\hspace{1em}k=\overline{1,N},$$

(54)

The set of p = N-m independent equations represents Maggi’s equations. Using the previous notation, it was possible to apply these equations to a single finite element to finally obtain Equation (40). Maggi’s equation represents another form of the Gibbs–Appell equation.

3.5. Kane’s Equations

Starting with the following equations:

$$\left({\overline{F}}_i-m_i{\overline{a}}_i\right)\delta {\overline{r}}_i=0\hspace{1em},\hspace{1em}i=\overline{1,N},$$

(55)

written for a mechanical system of N material points, described by a number of p generalized coordinates, we obtained:

$${\displaystyle \sum _{i=1}^N\left({\overline{F}}_i-m_i{\overline{a}}_i\right)\frac{\partial {\overline{r}}_i}{\partial q_k}}=0\hspace{1em};\hspace{1em}k=\overline{1,p},$$

(56)

The following relation was used [64,65,66]:

$$\frac{\partial {\overline{r}}_i}{\partial q_k}=\frac{\partial {\overline{v}}_i}{\partial {\dot{q}}_k}\hspace{1em};\hspace{1em}k=\overline{1,p},i=1,N;$$

(57)

Introducing Equations (57) in (56), we achieved:

$${\displaystyle \sum _{i=1}^N\left({\overline{F}}_i-m_i{\overline{a}}_i\right)\frac{\partial {\overline{v}}_i}{\partial {\dot{q}}_k}}=0\hspace{1em};\hspace{1em}k=\overline{1,p}.$$

(58)

If noted:

$$\frac{\partial {\overline{v}}_i}{\partial {\dot{q}}_k}=\frac{\partial {\overline{v}}_i}{\partial u_k}={\overline{v}}_i^{\left(k\right)}\hspace{1em};\hspace{1em}k=\overline{1,p}\hspace{1em};\hspace{1em}i=\overline{1,N},$$

(59)

Equation (58) becomes:

$${\displaystyle \sum _{i=1}^N{\overline{F}}_i^{}\frac{\partial {\overline{v}}_i}{\partial u_k}}={\displaystyle \sum _{i=1}^N{m}_i{\overline{a}}_i\frac{\partial {\overline{v}}_i}{\partial u_k}}\hspace{1em};\hspace{1em}k=\overline{1,p}\hspace{1em};\hspace{1em}i=\overline{1,N}$$

(60)

${\overline{F}}_i^{}$ being the external forces acting in the nodes.

Then, for an elastic finite element considered as a solid, Equation (60) becomes:

$$\hspace{1em}{\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}dm$$

(61)

Thus, it was possible to obtain, after some calculations, Equation (40).

4. Conclusions and Discussions

The use of an FEM in the analysis of an MBS with elastic elements is an increasingly common procedure in the analysis of mechanical systems of high complexity. The FEM has well-established and proven procedures for numerous applications and is extremely useful for such problems. The difficulty of obtaining the equations of motion lies in the complexity and particularities of the system being studied. There are several equivalent formalisms that allow the equations of motion to be written down, equivalent to each other. The analysis of an MBS with plane movement is a common objective in engineering applications, which is why the separate study of these systems is required. The particularities involved in planar motion help to obtain these equations in a simpler and easier form for the researcher to use. In this sense, this paper analyzed the method of writing the equations of motion using alternative methods from analytical mechanics. Among them, we listed the Lagrange, Gibbs–Appell, Hamilton, Kane, and Maggi equations.

The methods offered by analytical mechanics and previously presented for a two-dimensional finite element finally provided the equations of motion for a single finite element. These methods were equivalent from the point of view of the mechanical description, and the problem that a researcher posed was the choice of the most suitable method for the concrete application studied. The choice was determined by the previous experience of the researcher and his familiarity with the concepts presented in the paper.

• Lagrange’s equations is the method currently used by most researchers. We presented the motivation for this choice in the Introduction Section. The main advantages were familiarizing the researchers with the method and the fundamental notions used, the high degree of generality, the possibility of easy introduction into an algorithm, and its simplicity.
• Gibbs–Appell’s equations predicted the advantage of requiring a smaller number of differentiation operations. It was an economic advantage that, in the case of systems with a number of degrees of freedom, resulted in to reduced modeling and simulation costs. The main difficulty was the use of the notion of energy of accelerations, a concept with which the researchers were not very familiar. The GA method is a little-used method that has been reconsidered in recent years, due to the need to provide researchers with methods that ensure efficiency in terms of the time required to design a complex mechanical system.
• Maggi’s equations are essentially a development of the Lagrangian formalism. The study of non-holonomic systems lends itself very well to the application of this method. The multipliers are removed from the motion equations using a projection operator (orthogonal complement matrix). Knowing the kinetic and potential energy values and the links that exist between the nodes of the finite element network allows the equations of motion to be easily obtained. This formalism also provides a justification for the multiplier elimination procedures used empirically in FEM software.
• Kane’s equation method is very similar to Maggi’s equation method, from which it originates. It has recently been used in the automation industry and industrial robot applications. Kane’s equations represent a successful alternative, having the advantage of being economical in the study of systems with many degrees of freedom. The method represents a natural alternative for non-holonomic systems. The need to approach complex mechanical systems, which must operate at high speeds and in difficult conditions, encourages the alternative methods of description, among which Kane’s method is included, to be re-evaluated.
• Hamilton’s equations are starting to be reconsidered by researchers in the technological context in which we find ourselves. They are highly simple. The Hamiltonian is a scalar expressed in terms of generalized coordinates and their conjugate moments, and the degree of generality is very high. Moreover, the second-order equations of motion that are obtained in all other methods are replaced by first-order differential equations. It is true that the number of differential equations doubles, but by introducing new quantities with a physical significance (conjugate momentum).

From this analysis, it can be concluded that alternative methods from analytical mechanics should be reconsidered, since, for certain types of applications, they can present advantages in terms of modeling and numerical procedures. This occurs due to the current technological context, which requires increasingly complex systems and, as a result, even is more difficult to study and model, involving numerous parameters and long numerical analysis times.