**Sorin Vlase 1,2, Iuliu Negrean 2,3 , Marin Marin <sup>4</sup> and Maria Luminit,a Scutaru 1,\***


Received: 20 January 2020; Accepted: 13 February 2020; Published: 23 February 2020

**Abstract:** When analyzing the dynamic behavior of multi-body elastic systems, a commonly used method is the finite element method conjunctively with Lagrange's equations. The central problem when approaching such a system is determining the equations of motion for a single finite element. The paper presents an alternative method of calculation theses using the Gibbs–Appell (GA) formulation, which requires a smaller number of calculations and, as a result, is easier to apply in practice. For this purpose, the energy of the accelerations for one single finite element is calculated, which will be used then in the GA equations. This method can have advantages in applying to the study of multi-body systems with elastic elements and in the case of robots and manipulators that have in their composition some elastic elements. The number of differentiation required when using the Gibbs–Appell method is smaller than if the Lagrange method is used which leads to a smaller number of operations to obtain the equations of motion.

**Keywords:** Gibbs–Appell; energy of accelerations; finite element; nonlinear system; elastic elements; analytical dynamics; robotics

### **1. Introduction**

The common method to obtain the equations of motion for a finite element belonging to an elastic element of a multibody system is represented by Lagrange equations. In this way, the most important step in the analysis of the multi-body systems with elastic elements is solved, the rest of the procedures being the classical ones currently used in the finite element method FEA and verified by practice. So, it is possible to determine the dynamic response of a single finite element. Using the known assembly methods and introducing loads and proper boundary conditions, the set of differential equations describing the behavior of the entire elastic system is obtained.

In this process, the type of finite element chosen will determine the shape functions that will be used and which will give the final form of the coefficients of the matrix appearing in the equations of motion. In such an analysis, it is considered that the deformations are small enough to not have influence on the general, rigid motion of the multibody system. Thus, a large class of systems has been studied, having in their structure elastic elements modeled with one-dimensional, two-dimensional or three-dimensional finite elements. One-dimensional finite elements were among the first studied. We can mention here [1–7]. The method was then extended naturally to bi- and three-dimensional elements [8–10].

The basic hypothesis used was that the rigid motion of the multibody system is known and, as a consequence, the velocity and acceleration field of each solid of every element of the system. It is also considered that this rigid motion does not influence the elastic behavior of each element.

Customary Lagrange's equations were used to determine the dynamic response of a single element. For this, it is necessary to determine, as the first step, the Lagrangian composed by kinetic energy, internal energy and the mechanical work of the external loads. The studied elements were considered linear elastic. Thus, the evolution equations of a finite element are obtained. The use of Lagrange's equations proved to be a useful and practical method and was, generally, the only method used to solve this type of problem. Recent works in the field use Lagrange's equations for the study of multi-body systems with elastic elements [11–17]. Theoretically, however, other methods of solving these dynamic systems can be applied, which will ultimately lead to the same equations of motion. In this paper, we will deal with the use of Gibbs–Appell (GA) equations, which have advantages in terms of the number of operations required to obtain the final results and thus prove more economical.

One problem which led to the GA equations was to deduce the differential equations of motion, which are easier to apply in the case of non-holonomic systems but which can also be applied to the holonomic systems. This method was discovered by Gibbs in 1879 [18] and, independently, by Appell in 1899 [19]. The function of Lagrange or Hamilton is replaced, in this formalism, by the energy of accelerations and represents an application of the well-known Gauss principle of least constraint being equivalent to the other formulations currently used in classical mechanics.

The modern mechanical systems are characterized by ultra-fast movements and velocities and include serial structures of robots as well. The fast movements are generally characterized by a time variation law for accelerations and lead to the occurrence of accelerations of higher orders. The last decades are characterized by a sustained development of the use of robots in the industry. This involves studies that allow good control of such systems, especially when working with high speeds and forces. In the case of such systems, which give very high working speeds and especially important speed variations, accelerations become significant. In these cases, the GA method becomes interesting. This type of formalism is useful for a wide class of engineering problems. For example, it can be used to analyze rigid dynamical systems, considering quasi-velocities [20]. The equations of motion can be easily obtained, both for holonomic systems and for non-holonomic systems. It is thus possible to eliminate at this step Lagrange's multipliers and thus reduce the number of unknowns with which they work. The method can be applied both to linear and nonlinear problems. An example is presented in [21] where are obtained the motion equations for an N-link robot, with flexible elements. Experimental procedures validate the result obtained using the GA formalism. In the last period, this formalism tends to become a usual procedure. There are many papers presenting this kind of approach, like [22], where an application consists of an N-flexible-link manipulator used in the cooperative flexible mobile manipulator. The main advantage of recursive GA formulation is used to determine the motion equations without classical use of Lagrange multipliers. Some new methods to study the dynamic behavior of the multibody system have, as the central piece in approach, this formalism. As an example, in [23], such a system is defined using a classic set of coupled differential equations and a set of algebraic equations for expressing the liaisons (constraints). The use of Lagrange formalism asks a large number of derivatives in the governing equations. In these cases, the use of the recursive GA formalism becomes easier to apply. In this mode, it is possible to obtain the equations of motion automatically. In [24], the main advantage of this formulation used in applications is that it is not necessary to eliminate Lagrange's multipliers and, as a consequence, the computational complexity decrease. Identically the formalism is also used in [25]. Dynamical response and motion equations can be obtained via recursive GA formulation and represent the way used to reduce the computational effort, significant for this application. The main advantage is fewer mathematical calculations compared to other formulations. There are a lot of papers using this formalism in the study of the dynamical system [26]. Interesting results are presented in [27–29]. The present paper only aims at obtaining the equations of motion for a three-dimensional finite element with a general motion.

The problems regarding the accuracy of the results and the accuracy of the models used in the case of the mechanisms, which have significant displacements and deformations, were not developed in the paper. These very interesting and current study domains themselves represent a distinct field of research. In this paper, we only dealt with the presentation of an alternative method, which presents advantages of calculation, of determining the equations of motion for a single element. The other operations used then in the method of finite elements were considered the ones currently used in practice.

#### **2. Gibbs–Appell Formalism**

In the case of scleronomic liaisons, the acceleration of a point *i* of a system of material points having *n* degree of freedom has, generally, the expression:

$$\overline{a}\_{l} = \sum\_{k=1}^{n} \sum\_{j=1}^{n} \frac{\partial^2 \overline{r}\_{l}}{\partial q\_k \partial q\_j} \dot{q}\_k \dot{q}\_j + \sum\_{k=1}^{n} \frac{\partial \overline{r}\_{l}}{\partial q\_k} \ddot{q}\_k \tag{1}$$

where *r<sup>i</sup>* is the position vector of the point and *q<sup>j</sup>* , *j* = 1, *n* represent the independent coordinates. In the following, the dot placed above a scalar, vector or matrix size means the derivative of the respective size with respect to time. The notion energy of accelerations for a system of material points is defined as [27]:

$$S = \frac{1}{2} \sum\_{i=1}^{N} m\_i a\_i^2 \tag{2}$$

The form of this expression is similar to the expression of kinetic energy but the mechanical significance is different. In the case of a solid body, the definition can be extended to the whole domain of the body as:

$$S = \frac{1}{2} \int\_{V} \rho a^2 dV \tag{3}$$

where the acceleration of a certain point of the solid is similar to Equation (1):

$$\overline{a} = \sum\_{k=1}^{n} \sum\_{j=1}^{n} \frac{\partial^2 \overline{r}}{\partial q\_k \partial q\_j} \dot{q}\_k \dot{q}\_j + \sum\_{k=1}^{n} \frac{\partial \overline{r}}{\partial q\_k} \ddot{q}\_k \tag{3a}$$

Introducing Equation (3a) in Equation (3) we get:

*S* = <sup>1</sup> 2 R *V* ρ*a* <sup>2</sup>*dV* = <sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* . *qk* . *q<sup>j</sup>* + P*n k*=1 ∂*r* ∂*q<sup>k</sup>* .. *qk* 2 *dV* = = <sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* . *qk* . *q<sup>j</sup>* + P*n k*=1 ∂*r* ∂*q<sup>k</sup>* .. *qk* 2 *dV* = = <sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* . *qk* . *qj* 2 + 2 P*n k*=1 P*n j*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* . *qk* . *qj* P*n k*=1 ∂*r* ∂*q<sup>k</sup>* .. *qk* ! + P*n k*=1 ∂*r* ∂*q<sup>k</sup>* .. *qk* !2 2 *dV* = = <sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 P*n l*=1 P*n m*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* ∂ 2 *r* ∂*ql*∂*q<sup>m</sup>* . *qk* . *qj* . *ql* . *q<sup>m</sup>* + P*n k*=1 P*n j*=1 P*n l*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* ∂*r* ∂*q<sup>l</sup>* . *qk* . *qj* .. *q<sup>l</sup>* + P*n k*=1 P*n j*=1 ∂*r* ∂*q<sup>k</sup>* ∂*r* ∂*q<sup>j</sup>* .. *qk* .. *qj dV* = = <sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 P*n l*=1 P*n m*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* ∂ 2 *r* ∂*ql*∂*q<sup>m</sup>* . *qk* . *qj* . *ql* . *qm dV*+ +<sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 P*n l*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* ∂*r* ∂*q<sup>l</sup>* . *qk* . *qj* .. *ql dV* + <sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 ∂*r* ∂*q<sup>k</sup>* ∂*r* ∂*q<sup>j</sup>* .. *qk* .. *qj dV* = = <sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 P*n l*=1 P*n m*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* ∂ 2 *r* ∂*ql*∂*q<sup>m</sup> dV* . *qk* . *qj* . *ql* . *<sup>q</sup>m*+<sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 P*n l*=1 ∂ 2 *r* ∂*qk*∂*q<sup>j</sup>* ∂*r* ∂*q<sup>l</sup> dV* . *qk* . *qj* .. *ql*+ +<sup>1</sup> 2 R *V* ρ P*n k*=1 P*n j*=1 ∂*r* ∂*q<sup>k</sup>* ∂*r* ∂*q<sup>j</sup> dV* .. *qk* .. *q<sup>j</sup>* = *Eao*( . *q*) + *Ea*1( . *q*, .. *q*) + *Ea*2( .. *q*) (4)

where:

$$E\_{ao}(\dot{q}) = \frac{1}{2} \left( \int\_{V} \rho \sum\_{k=1}^{n} \sum\_{j=1}^{n} \sum\_{l=1}^{n} \sum\_{m=1}^{n} \frac{\partial^2 \tilde{r}}{\partial q\_k \partial q\_j} \frac{\partial^2 \tilde{r}}{\partial q\_l \partial q\_m} dV \right) \dot{q}\_k \dot{q}\_j \dot{q}\_l \dot{q}\_m \tag{5}$$

contains terms that have no accelerations,

$$E\_{a1}(\dot{q}, \ddot{q}) = \frac{1}{2} \left( \int\_{V} \rho \sum\_{k=1}^{n} \sum\_{j=1}^{n} \sum\_{l=1}^{n} \frac{\partial^2 \tilde{r}}{\partial q\_k \partial q\_j} \frac{\partial \tilde{r}}{\partial q\_l} dV \right) \dot{q}\_k \dot{q}\_j \ddot{q}\_l \tag{6}$$

represents the terms in which the accelerations intervene linearly and:

$$E\_{d2}(\ddot{\boldsymbol{q}}) = \frac{1}{2} \left| \int\_{V} \rho \sum\_{k=1}^{n} \sum\_{j=1}^{n} \frac{\partial \overline{r}}{\partial q\_{k}} \frac{\partial \overline{r}}{\partial q\_{j}} dV \right| \ddot{q}\_{k} \ddot{q}\_{j} \tag{7}$$

represents the quadratic terms in accelerations.

GA equations are given by relations:

$$\frac{\partial \mathcal{S}}{\partial \dot{\bar{q}}\_{j}} = Q\_{j} \quad j = \overline{1, n} \tag{8}$$

where *Q<sup>j</sup>* it represents the generalized force corresponding to the generalized coordinate *q<sup>j</sup>* .

#### **3. Basic Hypothesis**

Now consider a finite element modeling a domain of an elastic element from the structure of a multibody system. The aim of the paper is to determine the motion equations of this element. The generalized coordinates will be the displacements of the nodes. The displacements of an arbitrary point of the continuous domain modeled by the finite element shall be expressed using a finite number of nodal coordinates. The finite element shall be related to the mobile reference frame, which participates in the general rigid motion (Figure 1). With *vo*( . *Xo*, . *Yo*, . *Zo*), we shall note the velocity, and with *ao*( .. *Xo*, .. *Yo*, .. *Zo*), the acceleration of the origin of the local reference frame relative to the global reference frame OXYZ. We shall note with ω(ω*x*, ω*y*, ω*z*) the angular velocity and with ε(ε*x*, ε*y*, ε*z*) the angular acceleration. Generally, a multibody system consists of several solids. Every solid is characterized by these two vectors, and the number of pairs of such vectors is equal to the number of distinct bodies. The relations between the components of a vector in the local coordinate system {*t*}*<sup>L</sup>* and the components of the same vector in a global coordinate system {*t*}*<sup>G</sup>* are obtained using a matrix of rotation [*R*]:

$$\{t\}\_G = [\mathbb{R}]\{t\}\_L \tag{9}$$

**Figure 1.** Three-dimensional finite element. **Figure 1.** Three-dimensional finite element.

If an arbitrary point M of the finite element has a displacement *<sup>L</sup> f* changing into M', we may If an arbitrary point M of the finite element has a displacement *f L* changing into M', we may write:

$$\{r\_{M'}\}\_G = \{r\_O\}\_G + [R]\{ (r)\_L + \{f\}\_L \} \tag{10}$$

 *<sup>M</sup> <sup>G</sup> <sup>O</sup> <sup>G</sup> <sup>L</sup> <sup>L</sup> r r R r f* ' (10) where *<sup>M</sup> <sup>G</sup> r* ' is the position vector of point *M*' with components expressed in the global coordinate where {*rM*0}*<sup>G</sup>* is the position vector of point *M*' with components expressed in the global coordinate system. The continuous displacement field noted *f L* is approximated, using the finite element method, by relation:

$$\{f\}\_{L} = [\mathbf{N}(\mathbf{x}, y, z)] |\delta\rangle\_{L} \tag{11}$$

method, by relation: *<sup>L</sup> <sup>L</sup> f N*(*x*, *y*,*z*) (11) where the matrix elements [*N*(*x*, *y*, *z*)], the shape functions, will depend on the type of the finite element chosen and {δ}*<sup>L</sup>* are the independent nodal coordinates. The velocity of point M' will be:

$$
\langle \boldsymbol{\eta}\_{\rm MF} \rangle\_{\rm G} = \left\langle \dot{\boldsymbol{\tau}}\_{\rm MF} \right\rangle\_{\rm G} = \left\langle \dot{\boldsymbol{\tau}}\_{\rm O} \right\rangle\_{\rm G} + \left[ \dot{\boldsymbol{R}} \right] \langle \boldsymbol{\eta} \rangle\_{\rm L} + \left[ \dot{\boldsymbol{R}} \right] [\boldsymbol{\aleph}] \langle \boldsymbol{\delta} \rangle\_{\rm L} + \left[ \boldsymbol{R} \right] [\boldsymbol{\aleph}] \left\langle \dot{\boldsymbol{\delta}} \right\rangle\_{\rm L} \tag{12}
$$

element chosen and and the acceleration will be:

write:

where

*E*

$$\begin{aligned} \left\{ a\_{M'} \right\}\_{G} &= \left\{ \bar{r}\_{O} \right\}\_{G} + \left[ \bar{R} \right] \left\Vert r \right\Vert\_{L} + \left[ \bar{R} \right] \left\Vert f \right\Vert\_{L} + 2 \left[ \bar{R} \right] \left\Vert \left\{ \dot{f} \right\rbrack\_{L} + \left[ R \right] \left\Vert \ddot{f} \right\Vert\_{L} = \\ \left\{ \bar{r}\_{O} \right\}\_{G} &+ \left[ \bar{R} \right] \left\Vert r \right\Vert\_{L} + \left[ \bar{R} \right] \left\Vert N \right\Vert \left\Vert \delta \right\Vert\_{L} + 2 \left[ \bar{R} \right] \left\Vert N \right\Vert \left\Vert \dot{\delta} \right\Vert\_{L} + \left[ \bar{R} \right] \left\Vert N \right\Vert \left\Vert \ddot{\delta} \right\Vert\_{L} \end{aligned} \tag{13}$$

 *<sup>O</sup><sup>G</sup> <sup>L</sup> RN <sup>L</sup> RN <sup>L</sup> RN <sup>L</sup> r R r* 2 (13) These expressions will be used to determine the energy of accelerations. In the future, the index *G* mark the vector with the components expresses in the global reference system, the index L mark the vector with the components expresses in the local reference frame and the non-indexed vectors are considered to be written in the local coordinate system. The derivatives of the positional matrix *R* occur in the previous relations. These derivatives will define angular velocities and angular These expressions will be used to determine the energy of accelerations. In the future, the index *G* mark the vector with the components expresses in the global reference system, the index L mark the vector with the components expresses in the local reference frame and the non-indexed vectors are considered to be written in the local coordinate system. The derivatives of the positional matrix [*R*] occur in the previous relations. These derivatives will define angular velocities and angular accelerations. From the orthogonality conditions written for the positional matrix we have:

$$\begin{aligned} \left[\mathbb{R}\right] \left[\mathbb{R}\right]^T = \left[\mathbb{R}\right]^T \left[\mathbb{R}\right] = \left[\mathbb{E}\right] = \left[\begin{array}{cc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right] \end{aligned} \tag{14}$$

 0 0 1 where [*E*] is the unit matrix, the following relation can be obtained:

'

$$\begin{aligned} [\boldsymbol{\omega}]\_{\boldsymbol{G}} = \begin{bmatrix} \dot{\boldsymbol{R}} \end{bmatrix} [\boldsymbol{\aleph}]^T = \begin{bmatrix} 0 & -\omega\_{\rm \mathcal{G}G} & \omega\_{\rm \mathcal{G}G} \\ \omega\_{\rm \mathcal{G}G} & 0 & -\omega\_{\rm \mathcal{G}G} \\ -\omega\_{\rm \mathcal{G}G} & \omega\_{\rm \mathcal{G}G} & 0 \end{bmatrix} \end{aligned} \tag{15}$$

*yG xG*

which is a skew symmetric matrix angular velocity, with components in global coordinate system, corresponding to the angular velocity:

$$\{\omega\}\_{\rm G} = \left\{ \begin{array}{c} \omega\_{\rm xG} \\ \omega\_{\rm yG} \\ \omega\_{\rm zG} \end{array} \right\} \tag{16}$$

Obviously, we have, too:

$$[\boldsymbol{\omega}]\_{\boldsymbol{L}} = [\boldsymbol{\mathcal{R}}]^{\mathrm{T}}[\dot{\boldsymbol{\mathcal{R}}}] = \begin{bmatrix} 0 & -\boldsymbol{\omega}\_{\boldsymbol{z}\boldsymbol{L}} & \boldsymbol{\omega}\_{\boldsymbol{y}\boldsymbol{L}} \\ \boldsymbol{\omega}\_{\boldsymbol{z}\boldsymbol{L}} & 0 & -\boldsymbol{\omega}\_{\boldsymbol{x}\boldsymbol{L}} \\ -\boldsymbol{\omega}\_{\boldsymbol{y}\boldsymbol{L}} & \boldsymbol{\omega}\_{\boldsymbol{x}\boldsymbol{L}} & 0 \end{bmatrix} \text{and } \{\boldsymbol{\omega}\}\_{\boldsymbol{L}} = \begin{Bmatrix} \boldsymbol{\omega}\_{\boldsymbol{x}\boldsymbol{L}} \\ \boldsymbol{\omega}\_{\boldsymbol{y}\boldsymbol{L}} \\ \boldsymbol{\omega}\_{\boldsymbol{z}\boldsymbol{L}} \end{Bmatrix} \tag{17}$$

The angular acceleration operator is defined by:

$$\begin{aligned} \begin{bmatrix} \varepsilon \end{bmatrix}\_{\mathbf{G}} = \begin{bmatrix} \dot{\boldsymbol{w}} \end{bmatrix}\_{\mathbf{G}} = \begin{bmatrix} 0 & -\varepsilon\_{\mathbf{z}\mathbf{G}} & \varepsilon\_{\mathbf{y}\mathbf{G}} \\ \varepsilon\_{\mathbf{z}\mathbf{G}} & 0 & -\varepsilon\_{\mathbf{x}\mathbf{G}} \\ -\varepsilon\_{\mathbf{y}\mathbf{G}} & \varepsilon\_{\mathbf{x}\mathbf{G}} & 0 \end{bmatrix} = \begin{bmatrix} 0 & -\dot{\boldsymbol{w}}\_{\mathbf{z}\mathbf{G}} & \dot{\boldsymbol{w}}\_{\mathbf{y}\mathbf{G}} \\ \dot{\boldsymbol{w}}\_{\mathbf{z}\mathbf{G}} & 0 & -\dot{\boldsymbol{w}}\_{\mathbf{x}\mathbf{G}} \\ -\dot{\boldsymbol{w}}\_{\mathbf{y}\mathbf{G}} & \dot{\boldsymbol{w}}\_{\mathbf{x}\mathbf{G}} & 0 \end{bmatrix} = \begin{bmatrix} \ddot{\boldsymbol{R}} \end{bmatrix} \begin{bmatrix} \boldsymbol{R} \end{bmatrix}^{T} + \begin{bmatrix} \dot{\boldsymbol{R}} \end{bmatrix} \begin{bmatrix} \dot{\boldsymbol{R}} \end{bmatrix}^{T} \end{aligned} \tag{18}$$

and:

$$\begin{aligned} \begin{bmatrix} \varepsilon \end{bmatrix}\_{L} = \begin{bmatrix} \dot{\omega} \end{bmatrix}\_{L} = \begin{bmatrix} 0 & -\varepsilon\_{zL} & \varepsilon\_{yL} \\ \varepsilon\_{zL} & 0 & -\varepsilon\_{xL} \\ -\varepsilon\_{yL} & \varepsilon\_{xL} & 0 \end{bmatrix} = \begin{bmatrix} 0 & -\dot{\omega}\_{zL} & \dot{\omega}\_{yL} \\ \dot{\omega}\_{zL} & 0 & -\dot{\omega}\_{xL} \\ -\dot{\omega}\_{yL} & \dot{\omega}\_{xL} & 0 \end{bmatrix} = \begin{bmatrix} \dot{R} \end{bmatrix}^{T} \begin{bmatrix} \dot{R} \end{bmatrix} + \begin{bmatrix} R \end{bmatrix}^{T} \begin{bmatrix} \ddot{R} \end{bmatrix} \end{aligned} \tag{19}$$

It will result after some elementary calculus:

$$\left[\ddot{\boldsymbol{R}}\right]\left[\boldsymbol{R}\right]^{T} = \left[\boldsymbol{\varepsilon}\right]\_{\boldsymbol{G}} - \left[\dot{\boldsymbol{R}}\right]\left[\dot{\boldsymbol{R}}\right]^{T} = \left[\boldsymbol{\varepsilon}\right]\_{\boldsymbol{G}} - \left[\dot{\boldsymbol{R}}\right]\left[\boldsymbol{R}\right]^{T}\left(\left[\dot{\boldsymbol{R}}\right]\left[\boldsymbol{R}\right]^{T}\right)^{T} = \left[\boldsymbol{\varepsilon}\right]\_{\boldsymbol{G}} + \left[\boldsymbol{\omega}\right]\_{\boldsymbol{G}}\left[\boldsymbol{\omega}\right]\_{\boldsymbol{G}}\tag{20}$$

whereas:

$$[\boldsymbol{\omega}]\_{\mathcal{G}} = -[\boldsymbol{\omega}]\_{\mathcal{G}}^T \tag{21}$$

We have, too:

$$\boldsymbol{\dot{R}} \cdot [\boldsymbol{\mathcal{R}}]^T \begin{bmatrix} \ddot{\boldsymbol{R}} \end{bmatrix} = [\boldsymbol{\varepsilon}]\_L - \begin{bmatrix} \dot{\boldsymbol{R}} \end{bmatrix}^T [\boldsymbol{\mathcal{R}}] [\boldsymbol{\mathcal{R}}]^T \begin{bmatrix} \dot{\boldsymbol{R}} \end{bmatrix} = [\boldsymbol{\varepsilon}]\_G - \begin{pmatrix} [\boldsymbol{\mathcal{R}}]^T [\dot{\boldsymbol{R}}] \end{pmatrix}^T [\boldsymbol{\mathcal{R}}]^T [\dot{\boldsymbol{R}}] = [\boldsymbol{\varepsilon}]\_L + [\boldsymbol{\omega}]\_L [\boldsymbol{\omega}]\_L \tag{22}$$

These relations will be useful in the following considerations.

#### **4. Motion Equations**

In the following, the equations of motion for a finite element will be established using the GA formalism. For this, it is necessary to first determine the energy of the accelerations of a finite element. If we use Expression (13) determined previously for the acceleration of a point, the energy of the accelerations for the considered finite element is given by the expression [30]:

*E<sup>a</sup>* = <sup>1</sup> 2 R *V* ρ*a* 2 *<sup>M</sup>*0*dV* <sup>=</sup> <sup>1</sup> 2 R *V* ρ{*aM*0} *T* {*aM*0}*dV* = = <sup>1</sup> 2 R *V* ρ n .. *rO* o*T G* + {*r*} *T L* h.. *R* i*T* + {δ} *T L* [*N*] *T* h.. *R* i*T* + 2 n. δ o*T L* [*N*] *T* h . *R* i*T* + n.. δ o*T L* [*N*] *T* [*R*] *T x* n .. *rO* o *G* + h.. *R* i {*r*}*<sup>L</sup>* + h.. *R* i [*N*]{δ}*<sup>L</sup>* + 2 h . *R* i [*N*] n. δ o *L* + [*R*][*N*] n.. δ o *L dV* = 1 2 R *V* ρ n .. *rO* o*T G* n .. *rO* o *G* + 2{*r*} *T L* h.. *R* i*T*n .. *rO* o *G* + 2{δ} *T L* [*N*] *T* h.. *R* i*T*n .. *rO* o *G* + 4 n. δ o*T L* [*N*] *T* h . *R* i*T*n .. *rO* o *G* + 2 n.. δ o*T L* [*N*] *T* [*R*] *T* n .. *rO* o *G* + 1 2 R *V* ρ {*r*} *T L* h.. *R* i*T*h.. *R* i {*r*}*<sup>L</sup>* + 2{δ} *T L* [*N*] *T* h.. *R* i*T*h.. *R* i {*r*}*<sup>L</sup>* + 4 n. δ o*T L* [*N*] *T* h . *R* i*T*h.. *R* i {*r*}*<sup>L</sup>* + 2 n.. δ o*T L* [*N*] *T* [*R*] *T* h.. *R* i {*r*}*<sup>L</sup> x* 1 2 R *V* ρ {δ} *T L* [*N*] *T* h.. *R* i*T*h.. *R* i [*N*]{δ}*<sup>L</sup>* + 4 n. δ o*T L* [*N*] *T* h . *R* i*T*h.. *R* i [*N*]{δ}*<sup>L</sup>* + 2 n.. δ o*T L* [*N*] *T* [*R*] *T* h.. *R* i [*N*]{δ}*<sup>L</sup>* + 1 2 R *V* ρ 4 n. δ o*T L* [*N*] *T* h . *R* i*T*h . *R* i [*N*] n. δ o *L* + 4 n.. δ o*T L* [*N*] *T* [*R*] *T* h . *R* i [*N*] n. δ o *L* + n.. δ o*T L* [*N*] *T* [*R*] *T* [*R*][*N*] n.. δ o *L* (23)

The quadratic terms can be easily identified in the accelerations of the independent coordinates n.. δ o *L* , noted with *Ea*2:

$$E\_{d2} = \frac{1}{2} \int\_{V} \rho \Big( \left[ \ddot{\boldsymbol{\delta}} \right]\_{L}^{T} [\boldsymbol{N}]^{T} \left[ \boldsymbol{R} \right]^{T} [\boldsymbol{R}] [\boldsymbol{N}] \Big| \boldsymbol{\delta} \Big)\_{L} dV = \frac{1}{2} \Big[ \ddot{\boldsymbol{\delta}} \Big|\_{L}^{T} \Big( \int\_{V} \rho [\boldsymbol{N}]^{T} \left[ \boldsymbol{N} \right] dV \Big) \Big| \boldsymbol{\delta} \Big]\_{L} \tag{24}$$

the terms in which these accelerations appear linear *Ea*1:

*Ea*<sup>1</sup> = R *V* ρ n.. δ o*T L* [*N*] *T* [*R*] *T* n .. *r<sup>O</sup>* o *G* + n.. δ o*T L* [*N*] *T* [*R*] *T* h.. *R* i {*r*}*<sup>L</sup>* + n.. δ o*T L* [*N*] *T* [*R*] *T* h.. *R* i [*N*]{δ}*<sup>L</sup>* + 2 n.. δ o*T L* [*N*] *T* [*R*] *T* h . *R* i [*N*] n. δ o *L dV* = n.. δ o*T L* R *V* ρ [*N*] *T* n .. *r<sup>O</sup>* o *L* + [*N*] *T* ([ε]*<sup>L</sup>* + [ω]*<sup>L</sup>* [ω]*<sup>L</sup>* ){*r*}*<sup>L</sup>* + [*N*] *T* ([ε]*<sup>L</sup>* + [ω]*<sup>L</sup>* [ω]*<sup>L</sup>* )[*N*]{δ}*<sup>L</sup>* + 2[*N*] *T* [ω]*<sup>L</sup>* [*N*] n. δ o *L dV* (25)

and the terms that do not contain at all accelerations and which do not have importance in GA equations, *Ea*0:

*Ea*<sup>0</sup> = <sup>1</sup> 2 R *V* ρ n .. *rO* o*T G* n .. *rO* o *G* + 2{*r*} *T L* h.. *R* i*T*n .. *rO* o *G* + 2{δ} *T L* [*N*] *T* h.. *R* i*T*n .. *rO* o *G* + 4 n. δ o*T L* [*N*] *T* h . *R* i*T*n .. *rO* o *G* + 1 2 R *V* ρ {*r*} *T L* h.. *R* i*T*h.. *R* i {*r*}*<sup>L</sup>* + 2{δ} *T L* [*N*] *T* h.. *R* i*T*h.. *R* i {*r*}*<sup>L</sup>* + 4 n. δ o*T L* [*N*] *T* h . *R* i*T*h.. *R* i {*r*}*<sup>L</sup> dV*+ 1 2 R *V* ρ n. δ o*T L* [*N*] *T* h.. *R* i*T*h.. *R* i [*N*]{δ}*<sup>L</sup>* + 4{δ} *T L* [*N*] *T* h . *R* i*T*h.. *R* i [*N*]{δ}*<sup>L</sup>* + 4 n. δ o*T L* [*N*] *T* h . *R* i*T*h . *R* i [*N*] n. δ o *L dV* (26)

The potential energy (internal work) has a classical form:

$$E\_p = \frac{1}{2} \int\_V \langle \sigma \rangle^T \langle \varepsilon \rangle dV \tag{27}$$

where {σ} is the stress tensor and {ε} the strains tensor:

$$\begin{aligned} \{\sigma\} = \begin{Bmatrix} \sigma\_{\rm xx} \\ \sigma\_{yy} \\ \sigma\_{zz} \\ \tau\_{xy} \\ \tau\_{yz} \\ \tau\_{zx} \end{Bmatrix}; \{\varepsilon\} = \begin{Bmatrix} \varepsilon\_{\rm xx} \\ \varepsilon\_{yy} \\ \varepsilon\_{zz} \\ 2\varepsilon\_{xy} \\ 2\varepsilon\_{yz} \\ 2\varepsilon\_{zx} \end{Bmatrix} \end{aligned} \tag{28}$$

Writing the Hooke law as follows:

$$\{\sigma\} = [D] \{\varepsilon\} \tag{29}$$

and of the differential relations between strains and finite deformations:

$$\begin{aligned} \{\varepsilon\} = [a] \{f\} = [a] [N] \{\delta\} \end{aligned} \tag{30}$$

where [*a*] represents the differentiation operator [31] and, using Equations (13), (15) and (16), the internal work becomes [31]:

$$E\_p = \frac{1}{2} \int\_V \langle \delta \rangle\_L^T [N]^T [a]^T [D]^T [a] [N] \delta \rangle dV = \frac{1}{2} [\delta]\_L^T \left( \int\_V \left[ N \right]^T [a]^T [D]^T [a] [N] dV \right) \delta \rangle\_L \tag{31}$$

where [*k*] is the classical stiffness matrix:

$$\hat{\mathbf{p}}[k] = \bigcup\_{\mathbf{V}} \text{[N]}^T [a]^T [D]^T [a] [N] dV \tag{32}$$

*Symmetry* **2020**, *12*, 321

If noted with *p* = *p*(*x*, *y*, *z*) , the vector of body forces per unit volume, then the external work of these is:

$$\mathcal{W} = \bigcap\_{V} \{ p(\mathbf{x}, y, z) \}\_{L}^{T} \{ f \}\_{L} dV = \left( \int\_{V} \{ p(\mathbf{x}, y, z) \}\_{L}^{T} [N] dV \right) \vert \delta \rangle\_{L} = \{ q^{\*} \}\_{L}^{T} \vert \delta \rangle\_{L} \tag{33}$$

The concentrated forces *q L* in the nodes of the element will give an external work:

$$\mathcal{W}^c = \{q\}\_L^T \{\delta\}\_L \tag{34}$$

The equations of motion are obtained by applying the GA equations previously presented:

$$\left\{\frac{\partial E\_a}{\partial \ddot{\delta}}\right\}\_L - \{Q\}\_L = 0\tag{35}$$

By n ∂*E* ∂δ*e* o , we understand:

$$\begin{Bmatrix} \frac{\partial E}{\partial X} \end{Bmatrix} = \begin{Bmatrix} \frac{\partial E}{\partial \mathbf{x}\_1} \\ \frac{\partial E}{\partial \mathbf{x}\_2} \\ \vdots \\ \frac{\partial E}{\partial \mathbf{x}\_n} \end{Bmatrix} \\ \text{where } \mathbf{: \{X\}} = \begin{Bmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \vdots \\ \mathbf{x}\_n \end{Bmatrix} \tag{36}$$

In our case:

$$E\_a = E\_{a0}(\dot{q}) + E\_{a1}(\dot{q}, \ddot{q}) + E\_{a2}(\ddot{q})\tag{37}$$

and:

$$\langle Q \rangle\_L = [k] \langle \delta \rangle\_L + \langle q \rangle\_L + \langle q^\* \rangle\_L \tag{38}$$

After a series of elementary calculations and rearranging of terms, we get the equations of motion for the finite element considered:

$$\frac{\partial E\_{d2}}{\partial \left\{ \ddot{\delta} \right\}\_{L}} = \left( \int\_{V} \rho \left[ \mathbf{N} \right]^{T} \left[ \mathbf{N} \right] \,dV \right) \left\| \ddot{\delta} \right\|\_{L} = \left[ m \right] \left\{ \ddot{\delta} \right\}\_{L} \tag{39}$$

We shall note with h *N*(1) i , the row *i* of matrix [*N*] and with:

$$\begin{aligned} \left[m\right] &= \int\_{V} \rho \left[\mathbf{N}\right]^{T} \left[\mathbf{N}\right] dV = \int\_{V} \rho \left[\begin{array}{cc} \mathbf{N}\_{(1)}^{T} & \mathbf{N}\_{(2)}^{T} & \mathbf{N}\_{(3)}^{T} \\ & \mathbf{N}\_{(2)}^{T} & \mathbf{N}\_{(3)}^{T} \end{array}\right] \begin{cases} \mathbf{N}\_{(1)} \\ \mathbf{N}\_{(2)} \\ \mathbf{N}\_{(3)} \end{cases} \right] dV = \\ &= \int\_{V} \rho \left[\mathbf{N}\_{(1)}^{T} \mathbf{N}\_{(1)} + \mathbf{N}\_{(2)}^{T} \mathbf{N}\_{(2)} + \mathbf{N}\_{(3)}^{T} \mathbf{N}\_{(3)}\right] dV = \\ &= \int\_{V} \rho \mathbf{N}\_{(1)}^{T} \mathbf{N}\_{(1)} \, dV + \int\_{V} \rho \mathbf{N}\_{(2)}^{T} \mathbf{N}\_{(2)} \, dV + \int\_{V} \rho \mathbf{N}\_{(3)}^{T} \mathbf{N}\_{(3)} \, dV = \left[m\_{11}\right] + \left[m\_{22}\right] + \left[m\_{33}\right] \end{aligned} \tag{40}$$

where:

$$\mathbb{E}\left[m\_{\bar{l}\bar{l}}\right] = \int\_{V} \rho \mathbf{N}\_{(\bar{l})}^{T} \mathbf{N}\_{(\bar{l})} \,dV \tag{41}$$

therefore:

$$\begin{aligned} \left[m\_{11}\right] = \int\_{V} \rho \mathbf{N}\_{(1)}^{T} \mathbf{N}\_{(1)} \, dV, \ \left[m\_{22}\right] = \int\_{V} \rho \mathbf{N}\_{(2)}^{T} \mathbf{N}\_{(2)} \, dV, \ \left[m\_{33}\right] = \int\_{V} \rho \mathbf{N}\_{(3)}^{T} \mathbf{N}\_{(3)} \, dV \end{aligned} \tag{42}$$

$$\begin{split} \frac{\partial \mathbb{E}\_{\mathcal{Q}}}{\partial \left[\dot{\boldsymbol{\delta}}\right]\_{L}} &= \left( \int\_{V} \rho \left[ \mathbf{N} \right]^{T} dV \right) \left[ \ddot{r}\_{O} \right]\_{L} + \left( \int\_{V} \rho \left[ \mathbf{N} \right]^{T} \left( \left[ \boldsymbol{\varepsilon} \right]\_{L} + \left[ \boldsymbol{\omega} \right]\_{L} \left[ \boldsymbol{\omega} \right]\_{L} \right) \left[ r \right]\_{L} dV \right) + \\ & \left( \int\_{V} \rho \left[ \mathbf{N} \right]^{T} \left( \left[ \boldsymbol{\varepsilon} \right]\_{L} + \left[ \boldsymbol{\omega} \right]\_{L} \left[ \boldsymbol{\omega} \right]\_{L} \right) \left[ \mathbf{N} \right] dV \right) \delta \rangle\_{L} + 2 \Big( \int\_{V} \rho \left[ \mathbf{N} \right]^{T} \left[ \boldsymbol{\omega} \right]\_{L} \left[ \mathbf{N} \right] dV \right) \Big| \delta \rangle\_{L} = \\ & \left[ \mathbf{m} \right]^{T} \left[ \left[ \ddot{r}\_{O} \right]\_{L} + \left[ q^{i}(\boldsymbol{\omega}) \right]\_{L} + \left[ q^{i}(\boldsymbol{\varepsilon}) \right]\_{L} + \left( \left[ k(\boldsymbol{\omega}) \right] + \left[ k(\boldsymbol{\varepsilon}) \right] \right) \left[ \delta \right]\_{L} + \left[ c \right] \Big| \delta \right]\_{L} \end{split} \tag{43}$$

where it is denoted:

$$\begin{aligned} \begin{bmatrix} m\_O^i \\ \end{bmatrix} &= \int\_V \rho [\mathbf{N}]^T dV; \left\{ \dot{\boldsymbol{\eta}}^i(\varepsilon) \right\}\_L = \int\_V \rho [\mathbf{N}]^T [\varepsilon]\_L \langle \mathbf{r} \rangle\_L dV; \\ \begin{aligned} \{\dot{\boldsymbol{\eta}}^i(\omega)\}\_L &= \int\_V \rho [\mathbf{N}]^T [\boldsymbol{\omega}]\_L [\boldsymbol{\omega}]\_L \langle \mathbf{r} \rangle\_L dV; [\boldsymbol{k}(\varepsilon)] = \int\_V \rho [\mathbf{N}]^T [\boldsymbol{\varepsilon}] [\mathbf{N}] dV; \\ \left[ \boldsymbol{k}(\varepsilon) \right] &= \int\_V \rho [\mathbf{N}]^T [\boldsymbol{\omega}]\_L [\boldsymbol{\omega}]\_L [\mathbf{N}] dV; [\boldsymbol{c}] = \int\_V \rho [\mathbf{N}]^T [\boldsymbol{\omega}]\_L [\mathbf{N}] dV; \\ \left\{ \boldsymbol{m}\_{\text{ix}} \right\} &= \int\_V \rho [\mathbf{N}\_{(i)}]^T \mathbf{x} dV \quad ; \quad \left\{ \boldsymbol{m}\_{\text{iy}} \right\} = \int\_V \rho \left[\mathbf{N}\_{(i)}\right]^T \mathbf{y} dV \quad ; \quad \left\{ \boldsymbol{m}\_{\text{iz}} \right\} = \int\_V \rho \left[\mathbf{N}\_{(i)}\right]^T \boldsymbol{z} dV \end{aligned} \end{aligned} \tag{44}$$

The term *<sup>E</sup>a*<sup>0</sup> does not contain <sup>n</sup>.. δ o *L* and, consequently:

$$\frac{\partial E\_{a0}}{\partial \left| \ddot{\delta} \right|\_{L}} = 0 \tag{45}$$

Finally, the motion equations can be written in the condensed form:

$$\left[ [m] \begin{matrix} \bar{\delta} \end{matrix} \right]\_L + [\varepsilon] \begin{matrix} \dot{\delta} \end{matrix} \Big|\_L + ([k] + [k(\varepsilon)] + [k(\omega)]) \langle \delta \rangle\_L = [q]\_L + [q^\*]\_L - \left\{ q^i(\varepsilon) \right\}\_L - \left\{ q^i(\omega) \right\}\_L - \left[ m^i\_O \right] \left| \bar{r}\_O \right\rangle\_L \tag{46}$$

These equations are related to the local system of coordinates. Similar formulas can be obtained, and we consider the global reference frame. The matrix coefficients can be calculated after choosing the shape functions and the independent nodal coordinates for expressing the displacement of a point.

#### **5. Conclusions and Discussions**

In the classical Lagrangian formalism, to determine the equations of motion for a three-dimensional finite element, the main step of the calculation is the Lagrangian, based essentially on the kinetic energy of the deformable finite element. The Lagrangian is:

$$L = E\_{\mathbb{C}} - E\_p + \mathcal{W} + \mathcal{W}^{\mathbb{C}} \tag{47}$$

where *Ep*, *W*, *W<sup>c</sup>* are given, respectively, by the Relations (31), (33) and (34) and the kinetic energy has the expression:

*E<sup>c</sup>* = <sup>1</sup> 2 R *V* ρ*v* 2 *<sup>M</sup>*0*dV* <sup>=</sup> <sup>1</sup> 2 R *V* ρ{*vM*0} *T* {*vM*0}*dV* = = <sup>1</sup> 2 R *V* ρ n . *rO* o*T G* + {*r*} *T L* h . *R* i*T* + {δ} *T L* [*N*] *T* h . *R* i*T* + n. δ o*T L* [*N*] *T* [*R*] *T x* n . *rO* o *G* + h . *R* i {*r*}*<sup>L</sup>* + h . *R* i [*N*]{δ}*<sup>L</sup>* + [*R*][*N*] n. δ o *L dV* = 1 2 R *V* ρ n . *rO* o*T G* n . *rO* o *G* + 2 n . *rO* o*T G* h . *R* i {*r*}*<sup>L</sup>* + 2 n . *rO* o*T G* h . *R* i [*N*]{δ}*<sup>L</sup>* + 2 n . *rO* o*T G* [*R*][*N*] n. δ o *L d V*+ 1 2 R *V* ρ {*r*} *T L* h . *R* i*T*h . *R* i {*r*}*<sup>L</sup>* + 2{*r*} *T L* h . *R* i*T*h . *R* i [*N*]{δ}*<sup>L</sup>* + 2{*r*} *T L* h . *R* i*T* [*R*][*N*] n. δ o *L dV* + 1 2 R *V* ρ {δ} *T L* [*N*] *T* h . *R* i*T*h . *R* i [*N*]{δ}*<sup>L</sup>* + 2{δ} *T L* [*N*] *T* h . *R* i*T* [*R*][*N*] n. δ o *L* + n. δ o*T L* [*N*] *T* [*R*] *T* [*R*][*N*] n. δ o *L dV* (48)

Lagrange's equations are:

$$\frac{d}{dt}\left\{\frac{\partial L}{\partial \dot{\delta}\_{\varepsilon}}\right\} - \left\{\frac{\partial L}{\partial \delta\_{\varepsilon}}\right\} = 0\tag{49}$$

Leaving aside the other derivatives, which have a low weight in the computing economy, the term most significant in this formulation is the term kinetic energy. In this expression, there are four terms containing the vector <sup>n</sup>. δ o *L* and four terms containing the vector {δ}*<sup>L</sup>* . If we look at the expression of Lagrange's equations, it follows that we will have four differentiation operations for <sup>n</sup> ∂*L*/∂ . δ*e* o , another four for {∂*L*/∂δ*e*} and then the derivatives with respect to the time of the terms <sup>n</sup> ∂*L*/∂ . δ*e* o . *L* / *e*. Using the Gibbs–Appell equations that have the form of Equation (31), we have to differentiate the energy of the accelerations given by the Relations (24)–(26) and (37). The role of kinetic energy is

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 10 of 13

and then the derivatives with respect to the time of the terms

*L* / *e*

Using the Gibbs–Appell equations that have the form of Equation (31), we have to differentiate the energy of the accelerations given by the Relations (24)–(26) and (37). The role of kinetic energy is play in these equations by the energy of accelerations <sup>n</sup>.. δ o *L* . Note that of the three terms of kinetic energy, only two contain the vector of accelerations. The size *Ea*<sup>2</sup> has a single term that contains this vector and *<sup>E</sup>a*<sup>1</sup> has four terms that contain the vector <sup>n</sup>.. δ o *L* and which are differentiated. The total number of differentiation to be made in this case is five. play in these equations by the energy of accelerations *<sup>L</sup>* . Note that of the three terms of kinetic energy, only two contain the vector of accelerations. The size *Ea*<sup>2</sup> has a single term that contains this vector and *Ea*<sup>1</sup> has four terms that contain the vector *<sup>L</sup>* and which are differentiated. The total number of differentiation to be made in this case is five. The conclusion from these two approaches is that the number of operations to be performed

The conclusion from these two approaches is that the number of operations to be performed when using the Gibbs–Appell equations is lower than when Lagrange's method is applied. The use of this method has the advantage that we are familiar with the calculation of kinetic energy of a solid and to determine the energy of the accelerations it is sufficient to replace the speeds, in the kinetic energy formula, with the accelerations. The use of these equations will result in an economy in the computational effort we make to determine the equations of motion for a single finite element. when using the Gibbs–Appell equations is lower than when Lagrange's method is applied. The use of this method has the advantage that we are familiar with the calculation of kinetic energy of a solid and to determine the energy of the accelerations it is sufficient to replace the speeds, in the kinetic energy formula, with the accelerations. The use of these equations will result in an economy in the computational effort we make to determine the equations of motion for a single finite element. An application was made for a rod within a wind pumping mechanism (Figure 2). The

An application was made for a rod within a wind pumping mechanism (Figure 2). The mechanism proposed for the study is a mechanism with two degrees of freedom of type "kinematic chain closed by inertia" [32]. mechanism proposed for the study is a mechanism with two degrees of freedom of type "kinematic chain closed by inertia" [32]

(**a**) (**b**)

**Figure 2.** (**a**) Experimental test stand and (**b**) sketch of the mechanism with the positioning of the markers (dimensions are in mm). 1—engine; 2—rotating disk; 3—invertor; 4—pumps; 5—valves; 6—rod; 7—pendulum rod; 8—weight. **Figure 2.** (**a**) Experimental test stand and (**b**) sketch of the mechanism with the positioning of the markers (dimensions are in mm). 1—engine; 2—rotating disk; 3—invertor; 4—pumps; 5—valves; 6—rod; 7—pendulum rod; 8—weight.

Table 1 presents a comparison between the number of differentiation required when using Lagrange's method and the number of differentiation required if the Gibbs–Appell method is applied. It is noted that the number of these operations is reduced to less than half. For a large number of finite elements used, in the case of complex structures, the reduction can be significant and can reduce the computational effort. We mention that after obtaining the matrix coefficients of the system of differential equations, the next procedures in the two methods become identical, so Table 1 presents a comparison between the number of differentiation required when using Lagrange's method and the number of differentiation required if the Gibbs–Appell method is applied. It is noted that the number of these operations is reduced to less than half. For a large number of finite elements used, in the case of complex structures, the reduction can be significant and can reduce the computational effort. We mention that after obtaining the matrix coefficients of the system of differential equations, the next procedures in the two methods become identical, so also the necessary computer times.

also the necessary computer times. **Table 1.** Comparison between the Lagrange and Gibbs–Appell methods. Number of differentiations. **Number of Finite Elements Lagrange Gibbs–Appell** 5 288 120 Another problem that arises is the accuracy of the results. In the case of finite element analysis of a mechanism that can have large displacements, specific problems arise, very delicate, related to the accuracy of the results obtained and the correctness of the models used. During the present paper, we did not deal with these aspects, especially interesting and exciting. By themselves, these aspects define a field of research. The purpose of the paper was to provide an alternative way of writing the equations of motion, which would have the advantage of a smaller number of operations to perform.

10 528 220 15 768 320 20 1008 420 25 1248 520 The test mechanism on which the considerations were made, presented in Figure 2, can be considered, at this level as regards the accuracy of the results, accessible through the classical methods. In Figures 3 and 4, we presented the space of the possible positions that can be occupied by the elements of this mechanism, to justify our hypothesis.


**Table 1.** Comparison between the Lagrange and Gibbs–Appell methods. Number of differentiations. Another problem that arises is the accuracy of the results. In the case of finite element analysis Another problem that arises is the accuracy of the results. In the case of finite element analysis

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 11 of 13

*Symmetry* **2020**, *12*, x FOR PEER REVIEW 11 of 13

30 1488 620 40 1968 820

30 1488 620 40 1968 820

**Figure 3.** Sketch of the motion of the two degree of freedom system for the angular speed 16 rad/s. **Figure 3.** Sketch of the motion of the two degree of freedom system for the angular speed 16 rad/s. **Figure 3.** Sketch of the motion of the two degree of freedom system for the angular speed 16 rad/s.

**Figure 4.** Sketch of the motion of the two degree of freedom system for the angular speed 19 rad/s. **Figure 4.** Sketch of the motion of the two degree of freedom system for the angular speed 19 rad/s. **Figure 4.** Sketch of the motion of the two degree of freedom system for the angular speed 19 rad/s.

**Author Contributions:** Conceptualization, S.V. and I.N.; methodology, S.V. and M.L.S.; software, M.M.; validation, S.V., I.N. and M.M.; formal analysis, S.V. and M.L.S.; investigation, M.L.S.; resources, S.V.; writing—original draft preparation, S.V. and I.N.; writing—review and editing, M.L.S.; visualization, M.M.; supervision, I.N. and M.M.; project administration, S.V.; funding acquisition, S.V., I.N., M.M. and M.L.S. All authors have read and agreed to the published version of the manuscript. **Author Contributions:** Conceptualization, S.V. and I.N.; methodology, S.V. and M.L.S.; software, M.M.; validation, S.V., I.N. and M.M.; formal analysis, S.V. and M.L.S.; investigation, M.L.S.; resources, S.V.; writing—original draft preparation, S.V. and I.N.; writing—review and editing, M.L.S.; visualization, M.M.; supervision, I.N. and M.M.; project administration, S.V.; funding acquisition, S.V., I.N., M.M. and M.L.S. All authors have read and agreed to the published version of the manuscript. **Author Contributions:** Conceptualization, S.V. and I.N.; methodology, S.V. and M.L.S.; software, M.M.; validation, S.V., I.N. and M.M.; formal analysis, S.V. and M.L.S.; investigation, M.L.S.; resources, S.V.; writing—original draft preparation, S.V. and I.N.; writing—review and editing, M.L.S.; visualization, M.M.; supervision, I.N. and M.M.;project administration, S.V.; funding acquisition, S.V., I.N., M.M. and M.L.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding. **Funding:** This research received no external funding. **Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflicts of interest **Conflicts of Interest:** The authors declare no conflicts of interest **Conflicts of Interest:** The authors declare no conflict of interest.

#### **References References References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
