New Analytical Model Used in Finite Element Analysis of Solids Mechanics

In classical mechanics, determining the governing equations of motion using finite element analysis (FEA) of an elastic multibody system (MBS) leads to a system of second order differential equations. To integrate this, it must be transformed into a system of first-order equations. However, this can also be achieved directly and naturally if Hamilton’s equations are used. The paper presents this useful alternative formalism used in conjunction with the finite element method for MBSs. The motion equations in the very general case of a three-dimensional motion of an elastic solid are obtained. To illustrate the method, two examples are presented. A comparison between the integration times in the two cases presents another possible advantage of applying this method.


Introduction
The analytical methods of mechanics, clearly stated in Lagrange's book (1788) [1], were alternative formulations to Newton-Euler's equations of motion (EoM). Although initially the interest was more theoretical, these methods proved their usefulness when the complexity of the studied systems began to increase. Lagrangian mechanics experienced a continuous development in the 19th century and various formulations appeared in this context, such as Gibbs-Appell's (1879, 1899) methods [2,3], Hamilton's equations (1835) [4], Maggi's equations (1896) [5] and, in the middle of the last century, Kane's equations (1924) [6] (essentially a new form of writing Maggi's equations [7,8]). The interest in these possibilities of expressing the evolution equations of a mechanical system, starting from the fundamental principles was, for a long time, a theoretical interest, without significant practical consequences. A flexible multibody system (FMS) represents a group of rigid and deformable elements, each of which can undergo large translational and rotational movements. The components may come into contact with each other or with the environment. The classic kinematic linkages implied in engineering applications are: cylindrical and spherical joints, prismatic couplings, flat connections, screws, gears or cams. These elements can be connected in a closed loop (closed kinematic chain) as in classical mechanism theory or open (open kinematic chain) as in the case of manipulators. The flexible multibody dynamics (FMD) of an MBS impose computational strategies that are used to determine the dynamic response that involves knowing the time history of the mechanical behavior of the system, deformations, stresses due to external forces, liaisons and contour conditions.
In recent decades, however, the continuous industrialization that led to the emergence of machines and equipment operating at high speeds and accelerations, acting with intense forces, required the continuous improvement of the models and computational methods used. As a result, a reconsideration of the analytical methods previously used for the study of multibody systems followed and results were obtained which demonstrated that, in special cases, the use of one of the above methods may be more advantageous, both in terms of modeling and from the point of view of the calculation volume and, consequently, of the necessary computational time [7][8][9].
It is obvious that the need to solve practical problems has led to the development of the field and harmonized strategies. The excellent and well-documented work [9] makes a classification and ordering of the methods used. The modeling of flexible elements, within an MBS system, is one of the important steps for such a study.
In the analysis of an MBS, an inertial reference system is used as a global reference system to describe the motion of the entire system. Within the system, intermediate reference systems are linkage attached to each flexible element (component). This flexible elastic element has an average general motion in space. In this case, the movement of the flexible element in relation to the local reference system is mainly due to the elastic deformation. Researchers mainly use two types of local reference systems: floating and corotational frames. The first type, the floating frame, considers an average rigid body motion and follows this hypothetical rigid body. The second type, the corotational frame, considers an individual finite element and follows its average rigid body motion. However, it is possible to not use intermediate reference systems and the deformations reported directly to the global coordinate system. Our presentation uses the second type of reference frame.
We make a brief review of the main analytical methods used in FMD.
The use of Lagrange's equations allows for an algorithmization of dynamics problems, the steps necessary to be completed being [10,11]: • choosing generalized coordinates (which do not need to be Cartesian coordinates); • calculating kinetic energy; • establishing generalized forces.
After that, some pre-established differentiation calculations are performed, finally reaching the equations of motion of the studied mechanical system, which are second-order differential equations that allow for the determination of the unknown generalized coordinates. As a conclusion, the basic notions studied using analytical mechanics are the generalized coordinates (or velocities), the kinetic energy, the potential energy and the work.
Known since the early twentieth century, the equations called Gibbs-Appell (GA) equations [10] have rarely been used by researchers to solve practical engineering problems related to the evolution of mechanical systems. They have advantages in terms of the number of the equations needed for modeling and computational effort. The GA method was established as a convenient procedure for the study of non-holonomic systems. To apply this method, a new notion is introduced, with which researchers are less familiar, namely "acceleration energy", which will replace the Lagrangian one. The advantages of applying the method have been revealed by the latest studies on the dynamic analysis of mechanical systems that operate at high speeds and accelerations and operate with high forces. MBSs are very suitable for using the method [12]. The equations finally obtained coincide with those obtained by applying Lagrange's equations, but there is the advantage that the number of operations to be performed is smaller. The use of the method shows the possibility of eliminating Lagrange multipliers and, consequently, reducing the number of computation operations [13].
Due to these advantages, it may become a procedure and researchers are beginning to use it with predilection in the analysis of modern mechanical systems [14].
Maggi's equations are relatively recent, and represent an alternative form to other formalisms in analytical mechanics [10,11]. The feedback linearization technique in the case of control systems with which robots and manipulators are equipped naturally leads to the application of these types of equations. They are a simple way to obtain the equations of motion of the system in the case of non-holonomic liaisons. For large systems, defined by a large number of degrees of freedom (DOF), the application of the method allows for the shortening of the computational time required for the analysis. However, the method is not familiar to researchers and, as a result, is less used. Classical methods of obtaining EoM, such as Newton-Euler (NE) or Lagrange (LE), involve the liaisons forces. If the systems are large, this means a high number of calculus and, as a consequence, a high modeling and computation effort [15][16][17].
The method of Kane's equations, which is relatively new, has begun to be successfully applied in recent decades in the modeling and calculation of MBSs [18]. However, if the system has a large number of bodies, with many constraints, the use of this method becomes problematic in determining the dynamic response of such a system [19,20]. In essence, the Kane equations method eliminates some of the disadvantages of classical LE and NE methods and can be used for both holonomic and non-holonomic constraint systems. Kane's equations are an alternative form of Maggi's equations [21]. Some interesting applications of Kane's equations can be found in [22].
The analytical classic formalism in the study of MBSs was obtained using Lagrange's equations as starting point. The Lagrangian methods are an efficient way to carry out the dynamics analysis of a flexible mechanical system. The algorithms based on these equations introduce, through a set of multipliers, the liaison forces in the dynamic equations-usually in a linear form. However, there are also formalisms that have been developed, starting from Hamilton's equations. These are written in terms of generalized displacements and momenta. Hamiltonian mechanics achieve this precisely, providing a result of the modeling a system of differential equations of the first order with 2n unknowns. The first n unknowns are the generalized coordinates, the others being the generalized momenta [23,24]. The Hamilton method could have the advantage of providing us with a system of first-order equations, a system that can be used directly for numerical solving, without the need for the prior processing of obtained systems [25].

Preliminary
The first papers applying FEA to the study of an elastic MBS first considered a one-dimensional finite element with a two-dimensional motion [26][27][28]. The method used with predilection in all these works was the method of Lagrange's equations. This method involves the calculation of well-known notions of mechanics (kinetic energy, potential and work). The complexity of the studied systems has increased over time, as well as the types of finite elements used. It went from the study of two-dimensional motion to the study of three-dimensional motion and increasingly complex two-dimensional and three-dimensional finite elements were used [29][30][31]. Different aspects of MBSs with elastic elements have been analyzed lately, but the method of obtaining the equations of motion was also the method of Lagrange's equations [32][33][34][35].
Consider one single three-dimensional finite element referred to a local coordinate system, moving together with the element. The velocity of an arbitrary point of the element, M, after deformation, when it becomes M', is [36]: where: r M G is the velocity of the arbitrary chosen point; r O G is the velocity of the origin of the mobile reference frame;

•
[R] is the rotation matrix (transform a vector from the mobile reference frame to the fixed reference frame); • [N] is the matrix of the shape functions depending on the type of finite element used; • with the index L, a size is noted (vector, matrix), with the components expressed in the local reference frame, and with G, a size with the components expressed in the global reference frame.
The kinetic energy of the moving element is: where ρ is denoted the density of the material. The potential energy has the well-known form: The matrix [k] represents the rigidity of the element [18]. Concentrated forces q L give the work: and the vector of volume forces p = p(x, y, z) , the work: Now it is possible to define the Lagrangian equation for a moving finite element: The expressions of all these notions, in the case of the FEA of an elastic MBS, are presented in Appendix A.
The definition of the generalized momentum is: It results in: from where it is possible to obtain the vector of velocities . δ L expressed as function of p L : which is useful in the further considerations. The Lagrange equations, for a single finite element with three motions, are [31]: [m] ..

Hamiltonian Equations for a Finite Element
Hamilton's formalism has been used by many researchers in the analysis of elastic systems and Hamilton's equations have advantages in the case of mechanical systems with symmetries, which allow for obtaining prime integrals and, as a result, the size of the system decreases [37][38][39][40]. In our case, the Hamiltonian equation becomes: where, for the Lagrangian, the rel. (6) is used.
Hamilton's equations are [10]: In Table 1, with the results presented in the Appendix A, Hamiltonian derivatives are shown. Finally, the motion equations in the Hamiltonian formulation, as a set of first-order differential systems, are obtained:

Term of Hamiltonian Equation {
.
Hamilton's equations are attractive due to their simplicity and symmetry. Theoretically, these equations were analyzed from different points of view in different fields, from analytical mechanics to the geometry of vector spaces [9,23,37]. So, the way to solve it is the numerical way. In this case, the equations have an advantage because they are first-order differential equations and can be easily programmed and solved. In the case of second-order systems of differential equations (obtained if the Lagrange equation or other alternate analytical formulation is used), for the numerical solution, the transition to a first-order system of differential equations must be made [41]. Hamilton's equations give us these equations directly.

Examples
An analysis of the advantages that Hamilton's method could offer, compared to other classical methods and especially to Lagrange's method, is difficult to do and would probably require a systematic and distinct effort. However, in the following, we will try to make a division between the two methods, applied in a few simple cases.

a. Rotating beam
Consider a bar of length equal to L, mass m and cross section A, which rotates about an axis (Figure 1), perpendicular to the bar at its end. The beam is divided into n finite elements, the bar nodes are numbered from 1 to n + 1. In our calculus, we study the axial deformation of the beam. Equation (10), obtained using the method of Lagrange's equations, is used to write a program in Matlab for calculating axial deformations. In order to obtain the set of second-order differential equations applicable in this case, the parameters considered in the studied application are introduced. In order to apply Hamilton's equations (first-order differential equations) for the studied bar, Equation (13) is used. In this case, too, a program is written in Matlab. A comparison between the results using these two programs is presented.
the Lagrange equation or other alternate analytical formulation is used), for the numerical solution, the transition to a first-order system of differential equations must be made [41]. Hamilton's equations give us these equations directly.
In conclusion, Lagrange's equations are second-order differential equations which, in order to be solved, are to be transformed into a system of first-order differential equations. Hamilton's equations are first-order differential equations. In the form in which they are obtained, they can be used directly, without transformations, in a specific subroutine.

Examples
An analysis of the advantages that Hamilton's method could offer, compared to other classical methods and especially to Lagrange's method, is difficult to do and would probably require a systematic and distinct effort. However, in the following, we will try to make a division between the two methods, applied in a few simple cases.

a. Rotating beam
Consider a bar of length equal to L, mass m and cross section A, which rotates about an axis (Figure 1), perpendicular to the bar at its end. The beam is divided into n finite elements, the bar nodes are numbered from 1 to n + 1. In our calculus, we study the axial deformation of the beam. Equation (10), obtained using the method of Lagrange's equations, is used to write a program in Matlab for calculating axial deformations. In order to obtain the set of second-order differential equations applicable in this case, the parameters considered in the studied application are introduced. In order to apply Hamilton's equations (first-order differential equations) for the studied bar, Equation (13) is used. In this case, too, a program is written in Matlab. A comparison between the results using these two programs is presented.
. A sinusoidal excitation function will act at the right end of the beam. The results of the integration procedures are presented in Figures 2-4. The two formalisms applied offer the same results and the integration seems to be stable. What is different is the integration time. In Figure 5, the diagram of the computing time considering the two methods is presented. In our application, the length is L = 1 m, the mass is 1 kg, the section A = 0.2 cm 2 and the Young's modulus is 210 GPa. The beam rotates with an angular speed of 15 rad/s. The excitation axial force has a period of T = π/100. A sinusoidal excitation function will act at the right end of the beam. The results of the integration procedures are presented in Figures 2-4. The two formalisms applied offer the same results and the integration seems to be stable. What is different is the integration time. In Figure 5, the diagram of the computing time considering the two methods is presented. In our application, the length is L = 1 m, the mass is 1 kg, the section A = 0.

b. Beam in a general planar motion
An elastic beam of a mechanism with two degrees of freedom was considered. The beam had a planar motion. In order to facilitate the solution of the problem, the successive positions of some markers fixed on the bar were registered. Based on these records, the speeds and accelerations of the points where the markers were fixed were calculated and, based on them, the fields of the speed and acceleration for the beam were determined [22,35]. These results, obtained experimentally, were used as input data in Equations of motion (10) and (13). With the help of these values, integrations of these equations of motion were made.
Based on the mechanism presented in Figure 6, the positions of the markers were registered. Figure 7 shows some of the records. A coordinate system fixed to the bar was used to study the motion of different points on the beam. Figure 8 presents the acceleration field for two points on the bar ( Table 2). In order to determine the computational times for performing an integration over a relatively short time interval (0.1 s), the values indicated in Table 3 were obtained.

b. Beam in a general planar motion
An elastic beam of a mechanism with two degrees of freedom was considered. The beam had a planar motion. In order to facilitate the solution of the problem, the successive positions of some markers fixed on the bar were registered. Based on these records, the speeds and accelerations of the points where the markers were fixed were calculated and, based on them, the fields of the speed and acceleration for the beam were determined [22,35]. These results, obtained experimentally, were used as input data in Equations of motion (10) and (13). With the help of these values, integrations of these equations of motion were made.
Based on the mechanism presented in Figure 6, the positions of the markers were registered. Figure 7 shows some of the records. A coordinate system fixed to the bar was used to study the motion of different points on the beam. Figure 8 presents the acceleration field for two points on the bar ( Table 2). In order to determine the computational times for performing an integration over a relatively short time interval (0.1 s), the values indicated in Table 3 were obtained.

b. Beam in a general planar motion
An elastic beam of a mechanism with two degrees of freedom was considered. The beam had a planar motion. In order to facilitate the solution of the problem, the successive positions of some markers fixed on the bar were registered. Based on these records, the speeds and accelerations of the points where the markers were fixed were calculated and, based on them, the fields of the speed and acceleration for the beam were determined [22,35]. These results, obtained experimentally, were used as input data in Equations of motion (10) and (13). With the help of these values, integrations of these equations of motion were made.
Based on the mechanism presented in Figure 6, the positions of the markers were registered. Figure 7 shows some of the records. A coordinate system fixed to the bar was used to study the motion of different points on the beam. Figure 8 presents the acceleration field for two points on the bar ( Table 2). In order to determine the computational times for performing an integration over a relatively short time interval (0.1 s), the values indicated in Table 3 were obtained.               Following the analysis of the two simple cases presented in this section, it can be concluded that Hamilton's method seems to be more economical in terms of the computation time involved. However, the differences are not spectacular. With the increase in the studied system size, the difference between the two methods is accentuated. For this reason, it can be assumed that, for large systems, the method of Hamilton's equations could become interesting.
A study dedicated to this problem will clarify more precisely the differences in the application of the two methods. The simple comparison, made in this paper, seems to indicate that the differences will increase if the dimensions of the systems increase, which would be consistent with the observation that computer times depend, for this type of application, exponentially on the size of the system.

c. Analysis of the influence of different parameters on the obtained solutions
Equations (14) and (17) were considered in order to determine the displacement of the beam presented in Figure 1.
A series of simulations of the bar rotating around a vertical axis, perpendicular to the bar axis, were made to determine the influence of different parameters on the behavior of the bar and the possibility of obtaining results. Thus, when comparing the computer time using Lagrange's equations and Hamilton's equations, when the accuracy of obtaining the solution is similar, the calculation time for the second method is 6-11% lower. The graphical representations of the two methods are similar. If we vary the section of the bar and consider different rotations of the bar, we obtain the graphs in  Following the analysis of the two simple cases presented in this section, it can be concluded that Hamilton's method seems to be more economical in terms of the computation time involved. However, the differences are not spectacular. With the increase in the studied system size, the difference between the two methods is accentuated. For this reason, it can be assumed that, for large systems, the method of Hamilton's equations could become interesting.
A study dedicated to this problem will clarify more precisely the differences in the application of the two methods. The simple comparison, made in this paper, seems to indicate that the differences will increase if the dimensions of the systems increase, which would be consistent with the observation that computer times depend, for this type of application, exponentially on the size of the system.

c. Analysis of the influence of different parameters on the obtained solutions
Equations (14) and (17) were considered in order to determine the displacement of the beam presented in Figure 1.
A series of simulations of the bar rotating around a vertical axis, perpendicular to the bar axis, were made to determine the influence of different parameters on the behavior of the bar and the possibility of obtaining results. Thus, when comparing the computer time using Lagrange's equations and Hamilton's equations, when the accuracy of obtaining the solution is similar, the calculation time for the second method is 6-11% lower. The graphical representations of the two methods are similar. If we vary the section of the bar and consider different rotations of the bar, we obtain the graphs in Figures 9-12.           Varying the excitation frequency, we obtain the results presented in Figure 13. Varying the excitation frequency, we obtain the results presented in Figure 13.

Conclusions and Discussions
The main step in the FEA of an elastic MBS is to obtain the equations of motion for the finite elements chosen for discretization. If the formalisms of the consecrated analytical mechanics are used, a system of second-order differential equations is obtained. Once these equations are written, the next step involves assembling the equations and solving them. These are properly done according the usual and well-known procedures. The system of equations that is finally obtained is a system of second-order differential equations, with the nodal coordinates being unknown. To solve this system, it must first be transformed into a first-order system by introducing additional unknowns. This transformation is a costly process in terms of computational time.
Hamilton's equations offer, from the beginning, a system of first-order differential equations, with the unknowns consisting of generalized coordinates and generalized momenta. This proves to be an advantage as it is no longer necessary to transform the system of second-order equations into a Varying the excitation frequency, we obtain the results presented in Figure 13.

Conclusions and Discussions
The main step in the FEA of an elastic MBS is to obtain the equations of motion for the finite elements chosen for discretization. If the formalisms of the consecrated analytical mechanics are used, a system of second-order differential equations is obtained. Once these equations are written, the next step involves assembling the equations and solving them. These are properly done according the usual and well-known procedures. The system of equations that is finally obtained is a system of second-order differential equations, with the nodal coordinates being unknown. To solve this system, it must first be transformed into a first-order system by introducing additional unknowns. This transformation is a costly process in terms of computational time.
Hamilton's equations offer, from the beginning, a system of first-order differential equations, with the unknowns consisting of generalized coordinates and generalized momenta. This proves to be an advantage as it is no longer necessary to transform the system of second-order equations into a

Conclusions and Discussions
The main step in the FEA of an elastic MBS is to obtain the equations of motion for the finite elements chosen for discretization. If the formalisms of the consecrated analytical mechanics are used, a system of second-order differential equations is obtained. Once these equations are written, the next step involves assembling the equations and solving them. These are properly done according the usual and well-known procedures. The system of equations that is finally obtained is a system of second-order differential equations, with the nodal coordinates being unknown. To solve this system, it must first be transformed into a first-order system by introducing additional unknowns. This transformation is a costly process in terms of computational time.
Hamilton's equations offer, from the beginning, a system of first-order differential equations, with the unknowns consisting of generalized coordinates and generalized momenta. This proves to be an advantage as it is no longer necessary to transform the system of second-order equations into a system of first-order equations. As a result, the computational effort involved is smaller. The introduction of additional unknowns is done in the first step of modeling when the generalized moment is introduced. This method therefore proves to be better in terms of the simplicity of the solution.
Using Hamilton's equations to obtain the equations of motion for a finite element (thus an approximation of the equations of motion of the continuous medium) has an advantage in terms of the number of formal operations involved and the computational time required to solve. The paper applies this formalism for writing equations of motion in the case of the three-dimensional motion of an elastic linear solid and makes some illustrative applications for applying the method to two relatively simple examples.
In addition to these obvious advantages, which make it easier to obtain systems of differential equations in the models used, the comparison we made in the paper seems to indicate that computer times, applying this formalism, are smaller, which is an additional advantage.
where {ε} contains the distinct component of the strain tensor; {σ} contains the distinct component of the stress tensor.
Using the Hooke law and some elementary transformation, we finally obtain: Concentrated forces acting in knots, q L and volume forces p = p(x, y, z) give a mechanical work: and: The formula for the generalized momentum is: