Two-Dimensional Finite Element in General Plane Motion Used in the Analysis of Multi-Body Systems

A description of the motion equation of a single two-dimensional finite element used to model a multi-body system with elastic elements is made in the article. To establish them, the Lagrange’s equations are used. Obtaining the dynamic response of a system with deformable components has become important for technical applications in recent decades. These engineering applications are characterized by high applied loads and high acceleration and velocities. A study of such mechanical systems leads to the identification of different mechanical phenomena (due to high deformations, resonance phenomena, and stability). Coriolis effects and relative motions significantly modify the motion equations and, implicitly, the dynamic response. These effects are highlighted in this paper for plane motion.


Introduction
Finite Element Analysis (FEA) of elastic bodies with overall rigid motion began in the 1970s. A number of papers analyzed different cases and types of finite elements used in this type of analysis [1][2][3]. One-dimensional and three-dimensional finite elements have been developed with plane motion or three-dimensional motion [4,5]. The field has remained interesting for researchers so far, with many papers on different aspects of the issue being published. The new field of multi-body systems (MBS), which has developed in recent decades, has reached the study of mechanical systems with (some) linear/nonlinear elastic elements. Why is this type of analysis necessary? If the velocities, accelerations, or loads that a constitutive element of a multi-body system supports are high, the elasticity of the elements may cause instability phenomena. The classical hypothesis of the rigid elements, commonly used in the dynamic analysis and synthesis of multi-body systems, no longer corresponds to the reality. The phenomena of resonance and loss of stability represent classic forms of manifestation of the elasticity properties. The classical approach involves the use of continuous mechanical models. Unfortunately, this type of analysis leads to a nonlinear system of differential equations that cannot be solved analytically in common engineering applications. The most convenient approach is FEA, a time-validated and frequently used method for analyzing mechanical systems with elastic elements. However, classical models of static or quasi-static systems do not correspond to the analysis of multi-body systems, which have rigid motion. The relatively large amount of motion between elements and Coriolis effects lead to the existence of additional terms in motion equations that cannot be neglected. These terms are usually strongly nonlinear. As a consequence, proper modeling of the problem is required [6,7]. The many aspects of a one-dimensional problem are described in [8].
In the following research, we develop a method for a two-dimensional finite element with plane motion. The motion equations for this type of element are established.

Two-Dimensional Model
Let us consider a two-dimensional finite element with a body with the general rigid plane motion of an MBS. To determine the governing equations for this element, the Lagrange's equations are used. The main phase in this type of approach is to obtain the Lagrangian for the one two-dimensional finite element. The terms necessary to build the Lagrangian are the kinetic energy of the two-dimensional finite element, the internal energy, and the work done by the concentrated and distributed loads. The shape function (determined in every case by the type of finite element chosen and the hypothesis made) will finally determine the form of the matrix coefficients and, as a consequence, the differential system of the governing motion equations. A basic hypothesis in this study is that the deformations of the elements of MBS are small enough not to influence the plane rigid motion of the system. Both the problem of the plane multi-body rigid motion of the system and velocities and the acceleration field distribution for each element can be solved using the classical method (see [9][10][11]).
The finite element chosen is related to the local reference frame Ox 1 x 2 (participating with the finite element in the general plane rigid motion) and to the global reference frame O'X 1 X 2 . It is considered that the vector velocity and acceleration of the origin O related to the global reference system are known from a previous dynamical analysis. All the bodies are considered to be rigid as is the angular velocity and angular acceleration of the local reference system. In this way, the fields of velocity and acceleration are known for each point of each constitutive element of the MBS [12][13][14][15][16][17][18].
For one single two-dimensional finite element, the generalized independent coordinates can be the nodal displacements, depending on the shape function and hypothesis used to express the displacement of a point in the finite element. The final number of independent coordinates will depend on this hypothesis and shape function.
Consider the displacement {δ(u, v)} of an arbitrary point M of the domain of the finite element. If we use the shape function N ij and the vector of the nodal displacements δ e,j , the displacements in the local coordinate system are where the nodal displacement vector of the finite element is labeled e, and δ e,j depends on the independent coordinate chosen to define the element. The number of independent coordinates is p. For a triangular finite element we have The two lines of the shape function matrix N correspond to the displacements u and v and are named N (u) = N (1) and N (v) = N (2) : The displacement of the point M(x 1 ,x 2 ), which is arbitrarily chosen, becomes, after deformation, M'(x' 1 ,x' 2 ) and can be expressed through its components: or with respect to the global reference system: or

Lagrangian of an Element
Using (6), the velocity of point M' can be obtained after differentiation with respect to the time: .
The kinetic energy for one single element due to this velocity is where t is the thickness of the element, and ρ is the mass density.
In plane motion, the rotation of the local system of coordinates is expressed by the plane rotation tensor r ij . The columns of this tensor define the positions of the unit vectors of the local reference frame Oxyz. For plane rotation, these coefficients can be expressed as The ortho-normality condition of these unit vectors leads to r ij r k j = r jk r ji = δ ij , i jk = 1, 2 where δ ij is the Kronecker delta (derived from the general three-dimensional transformation). By differentiation, the following equation will result: . r ij r k j + r ij . r k j = 0, i, k = 1, 2. Denote r ij r k j , .
the skew-symmetric tensor angular velocity. The relation (12) becomes To this corresponds the angular velocity vector components defined by The angular velocity vector and the angular acceleration vector have, in plane motion, the same components in both the local reference system and the global reference system.
The angular acceleration skew symmetric operator it is defined as .
The angular acceleration vector is defined by The other components are null. We shall have . . .
from where ..
The internal energy stored is obtained via the relation From (20), it is possible to obtain the stiffness matrix depending on the hypothesis and type of finite element (and consequently, the shape functions) chosen. Here, we present a case as an example; any other case can be treated in the same manner. The relationship between stress and strain is If we choose, for example, a rectangular finite element with sides a and b, a displacement field of the form can be used to obtain the stiffness tensor. The boundary conditions (x = 0, ) lead to the following fields for the displacements: Applying the classical laws of linear elasticity [2], the relation between the strain and nodal displacement are, in this case, Using (20), (21), and (24), after some elementary calculus, it is possible to obtain the internal energy of the form (see [2]) The external work of the concentrated load is W = q e,j δ e,j , i = 1, 2, j = 1, p.
Here, q ej is the vector of the concentrated loads acting in nodes.
The external work of the distributed loads is where the notation is used. We have all the data to build the Lagrangian [8,19] or, taking into account the relations (9) δ e,m dA − 1 2 δ e,i k e,ij δ e,j + q * eL,j δ e,j + q e,i δ e,i . (30)

Motion Equations
Theorem: The motion equations written in the local coordinate system for two-dimensional finite element take the form . ..
We obtain the motion equations (an extended presentation can be found in Appendix A): ..
x jo . (35) Using the notation mentioned above, it is possible to obtain (30).

Conclusions
In engineering practice, there are frequent cases in which the elements of a device or a machine are required to operate with high forces or have high accelerations and speeds. In these cases, the elasticity of the elements can be manifested in such a way as to negatively influence the working process of the assembly. An analysis of the MBS motion in such conditions is required. In this work, we determined the motion equations for a two-dimensional finite element with plane motion. This case often occurs in practice. Most of the technical systems have elements operating in planar motion. In this case, a system of differential equations with matrix coefficients that can be strongly nonlinear is obtained. The mass tensor and the classical rigidity tensor are symmetrical matrices, while the Coriolis effects tensor is skew-symmetric. There are also terms that change the classical stiffness of the element and additional inertial terms. If we use FEA to obtain the motion equations for one single element in plane motion, we obtain a system of second-order differential equations, as presented in (30). Generally, this system is strongly nonlinear, with the matrix coefficients of the system depending on time and on the position of the system. For a usual engineering application, some aspects of these equations can cause difficulties. We highlight some properties of this system:

•
The classical inertia tensor m e,ij is a symmetrical tensor; • The damping tensor c ω e,ij is a skew symmetric tensor; this represents the effect of the accelerations due to the Coriolis effects (relative motions with respect to the mobile reference co-ordinate system; • The stiffness tensor k e,ij is a symmetric tensor; this tensor is modify by additional terms depending on the general plane rotation of the element, becoming k e,ij + k ε e,ij + k ω 2 e,ij ; • The vector of the generalized loads contains some supplementary terms due to inertia of finite elements being in rigid motion; these are −q ε e,i − q ω 2 e,i − m o e,ij ..
x jo .
The common method used to solve this system is to linearize this system, considering the tensor coefficients as being constant for very short time intervals (rigid motion freezing). In this way, it is possible to obtain a system of differential equations with constant coefficients. To solve this, usual and well-known methods can be used. There can be singularities due to the inertia term affecting the stiffness matrix [20].