THE TRAPEZOIDAL FINITE ELEMENT IN ABSOLUTE COORDINATES FOR DYNAMIC MODELING OF AUTOMOTIVE TIRE AND AIR SPRING BELLOWS. PART I: EQUATIONS OF MOTION

Equations of motion of a finite element in absolute coordinates including mass matrix, generalized inertia and internal forces are derived. A trapezoidal element for dynamic models of flexible shells in the shape of surface of revolution is considered. The element can be used for modeling dynamics of automotive tire and air spring bellows and some other flexible elements of transport systems undergoing large elastic deflections.


INTRODUCTION
Including flexible bodies in dynamics models of multibody systems is a very useful approach in simulation of transport systems such as railway, monorail or automotive vehicles. This feature is available in any commercial multibody software. The Craig-Bampton method [1] is the most frequently used technique for modeling flexible bodies, allowing large spatial displacements but small deformations and deflection. This method is characterized by efficient decrease in the degrees of freedom in simulations of flexible body dynamics. As a rule, it requires the generation of a linear finite element model of a body in one of the external FEA programs. An alternative method is offered by the absolute nodal coordinate formulation (ANCF) [2], the advances of which are linked to the correct modeling of flexible bodies undergoing large deflections. Use of the absolute node coordinates allows direct and non-iterative-derived explicit equations for such motions of flexible bodies. Typical examples are a conveyor belt running over a pulley or a drill string rotating in a curved borehole. Some beam, plate and shell elements derived with the ANCF are described in [3][4][5][6].
In this paper, we consider a mathematical model of the dynamics of a flexible body subjected to arbitrary 3D displacements, small strains and large flexible deflections. Dynamic equations are derived in coordinates of node-fixed frames, and the Craig-Bampton approach is used for the generation of equations of motion of each of the finite elements. Large displacements are allowed for the body, but within a single element, both the strains and the flexible deflections are considered to be small. In fact, we apply the Craig-Bampton procedure to finite elements instead of the whole flexible body. On the one hand, this leads to a considerable growth of the degrees of freedom and CPU expenses during the simulation in comparison with the original Craig-Bampton method. On the other hand, this formalism allows the correct description of both the gross 3D body motion and the large deflections of flexible structures under the assumption of small strains.
The presented method differs from the usual ANCF approach in the type of node coordinates as well as in the expressions for equations of motion. The expected advances of the method in comparison with the ANCF consist of decreasing the computational costs and in a transparent structure of the elements of the equations of motion, which is important for computation of Jacobian matrices of internal elastic and damping forces.

D. Pogorelov, A. Rodikov
In Section 2, we derive equations of motion of a free finite element in terms of accelerations of the node-fixed systems of coordinates. In Section 3 of this paper, a model of a four-node finite element in the shape of an isosceles trapezoid is considered. The element is suitable for simulation of dynamics in absolute nodal coordinates of flexible shells, which are surfaces of revolution. In the future, we plan to apply it for simulation of the dynamics of tires, air spring bellows and some other flexible elements of transport vehicles.
Verification tests based on computation of natural frequencies and modes for some flexible bodies, computation of thin-plate buckling and comparison of large deflection results for rectangular and annular plates with theoretical ones and with the results of other authors are considered in Part II of this paper.

EQUATIONS OF MOTION OF A FINITE ELEMENT IN ABSOLUTE COORDINATES
The main idea for derivation of equations of motion of a flexible body in multibody system dynamics consists of splitting its spatial movement on a gross motion of a reference frame and small elastic deformations [7]. It is supposed that the flexible body is divided into a set of finite elements. Consider a finite element, which moves freely in 3D space, Fig. 1. Let us assume that n nodes are connected with the element. Positions of the node-fixed systems of coordinates SCi relative to the inertial frame of reference SCO are determined by the radius-vectors and the direct cosine (rotation) matrices . Hereafter in the paper, the abbreviation SC is used for the system of coordinate designations. It is supposed that the element position and deformed state are uniquely determined by the positions of the node-fixed frames. In a particular case, the node positions can be defined by the radius vectors only, and the rotation matrices are not used. The floating frame SCF is a system of coordinates connected with the element. SCF follows the gross motion. Small deflections of element points relative to the SCF depending on the point coordinates relative to SCF describe the deformed state. To derive the equations of motion of the element in the absolute coordinates, we should first express the radius-vector and rotation matrix of SCF in terms of the SCi coordinates. Several methods are available for selection of the floating frame position. The corotational description is widely used for geometric nonlinear analysis of shells and beams with large displacements and rotations and small strain [8][9][10][11]. Use of an 'average' node position as the floating frame is described in [12,13]. The simplest method consists of connection of SCF with one of the SCi, i.e., the floating frame coincides with one of the node systems of coordinates. This method is applied to simulation of The trapezoidal finite element in absolute coordinates for…

143.
flexible beam structures in [14]. In this paper, we apply the Craig-Bampton formalism for localization of SCF.
Consider first the element as a rigid body when the elastic deformations are zeroes. In this case, the SCF must follow the element motion so that coordinates of the element points relative to the SCF are constant. The SCF position relative to the rigid element can be arbitrary; for instance, it can coincide with the principal central axes of inertia of the element. Now, let us add small translations of the node systems of coordinates SCi relative to the floating frame so that are the vectors of node i relative displacements and are the vectors of small rotations of SCi relative to SCF. For small rotations, the components of vector are equal to the angles of rotation of SCi about the axes of SCF.
Following the standard procedure of the finite element method, introduce an approximation of the elastic displacements of the element points relative to the floating frame using a matrix of appropriate shape functions. The matrix consists of 2n blocks of size each corresponding to the translational and rotational degrees of freedom in node i (1) With the matrix of shape functions, the deflections of the element points depend linearly on the node displacements (2) Here, I is the 3´ 3 identity matrix, are the vectors from the SCF origin to the nodes and is the matrix of the node displacements (3) The superscript T denotes the transpose of a matrix.
In general, the nodal displacements include movements of the element as a rigid one. For example, the displacement (4) is the element shift along the x axis on . Excluding these movements by a linear transformation Before deriving equations of motion, a method for computation of the transformation matrix should be described. Following the Craig-Bampton procedure, consider the mass and stiffness matrices of the free finite element derived with the classical linear FEM for the given node coordinates and shape function . For example, the mass matrix is derived from the integral over the element (7) Solution of the natural frequency problem yields the natural frequencies and the free-free modes , j=1…6n. Six frequencies corresponding to the rigid body motions are zeroes. Using the normalization of the modes and skipping the modes related to zero frequencies, we obtain the desired transformation matrix (9) with the following properties: (10) where denotes a diagonal matrix. Consider the principle of virtual work for generation of equations of motion of the finite element [7], taking into account the work of inertia and elastic forces only (11) Here and a are the virtual displacement and acceleration of a point, respectively, whose position relative to SCF in Fig. 1 is determined by the vector r; dm and dV are the mass and volume of an infinitesimal part of the element; and are the elements of the stress and strain tensors.
The position of an arbitrary point of the element relative to SCO according to Fig. 1 is (12) which gives the expressions for the virtual displacement and acceleration (13) Here and below, the ~ sign over the vector indicates a skew-symmetric matrix generated by the vector and used for the matrix notation of the vector product, for example (14) are the vectors of angular velocity, angular acceleration and the virtual vector of rotation of SCF, respectively. Substitution of these formulas into Eq. (11) and equating to zero terms for independent virtual displacements yield the equations of motion of the free element.
Some small terms are omitted for the sake of simplicity: the dependence of the inertia coefficients (elements of the mass matrix) on elastic coordinates as well as the terms in the inertia forces depending on .
(15) Here, are the mass and the matrix of inertia tensor of the element, respectively, and is the radius-vector from the SCF origin to the element center of mass, and the following matrices are introduced:  (15) and (18) in the matrix form (19) we obtain the equations of motion of the element in the absolute nodal accelerations as follows: (20) It is clear that both the mass matrix of the element and the transformation matrix are variable and must be computed at each step of the integration process. To speed up the computation, we note that the matrix in SCF is constant, and the matrix in SCF is constant as well if . If we neglect the dependence of the mass matrix on the elastic coordinates, it will be constant in SCF and can be computed once. Nevertheless, the transformation of the mass matrix to SCO during the simulation is required at each integration step. Generation of the global equations of motion for a meshed flexible body based on Eq. (20) written for each of the finite elements can be done in a standard manner and consists of summation and positioning of blocks of the mass matrices and generalized forces according to the mesh topology and node numbering.
Taking into account geometric nonlinearity within one element leads to the element stiffness matrix of the following form: Here, is the linear stiffness matrix, is the geometric stiffness matrix, which is a function of stresses in the element, and is the stiffness matrix of large displacements, which takes into account the nonlinear terms of the Green-Lagrange strain tensor. Matrices and for the plate element are described in [15]. The above transformation matrix is computed taking into account the linear matrix , whereas the full matrix is used for evaluation of the generalized internal forces instead of in Eq. (15).
As shown in Part II of the paper, large deflections as well as buckling of thin-walled flexible bodies can be successfully modeled without nonlinear terms , of the stiffness matrix. Adding these matrices could improve the accuracy of the element and allows us to decrease the number of finite elements in the models.

TRAPEZOIDAL FINITE ELEMENT
Let us consider a flexible shell, whose shape is a surface of revolution. This type of shell can be used for modeling automotive tires and air spring bellows. It is clear that a revolution surface can be easily meshed into plane isosceles trapezoids as shown in Fig. 2.

147.
Consider a mathematical model of a 3D uniform trapezoidal finite element using the flexible thinplate theory. The geometry parameters of the element are shown in Fig. 3. Let us introduce the Cartesian system of coordinates, whose origin coincides with the trapezoid center, and the z axis is perpendicular to the element plane. Four nodes are located in the trapezoid vertices. Each of the nodes has six degrees of freedom, corresponding to the Cartesian coordinates and the angles of rotation about the x, y, z axes Here and below in this section, the subscript i is the node number according to Fig. 3. The node numbers are denoted in this figure as circles.
Let a, b be the median and the height of the trapezoid, and D be a half of the difference in the length of the median and the upper side, We use the known polynomials as the in-plane interpolation functions (25) For the bending degrees of freedom, we consider two variants of the shape functions. The first variant (V1) uses the products of the Hermite polynomials for a beam element [16] (26) Here, are the following polynomials: (27) The second variant (V2) of the bending shape functions [16] corresponds to the following polynomials: (28) The matrices for the isosceles trapezoid element have the same structure (23), and most of the shape functions for the isosceles trapezoid are equal to those for the rectangle  . A square grid with a 0.5m step is drawn in the xy plane on each of the pictures. It is important for the continuity condition that in the case of a combined rotation about x and y axes in node k resulting in the rotation about the trapezoid inclined side, the displacement in the z direction along the rotation axis vanishes, Fig. 5. Now, the mass and stiffness matrices can be computed according to the standard methods [16,17]. For example, the expressions for one of the elements of the mass and stiffness matrix according to the classical theory of a thin uniform plate are   Here, h is the plate thickness, E is the Young's modulus and n is the Poisson's ratio. The elements of the mass and stiffness matrix as well as all other matrices in Eq. (15) can be computed both analytically and numerically.
Free-free modes of the isosceles trapezoid for the bending frequencies are shown in Fig. 6.
Modeling of tire and air spring bellows requires a multilayer finite element, which can be developed with two different approaches [18]. The first method is based on applying different material models to each of the layers; the second method uses a special material model for a multilayer shell. Application of these methods to automotive tires can be found in [18,19]. Equations of motion of flexible bodies are generated in Universal Mechanism software according to the algorithms formulated above. The bodies can have the following shapes (Fig. 2, 7): rectangular plate (Fig 7a), annular plate (Fig 7b,c), torus (Fig 7d,e), circular cylinder and cone (Fig. 7f,g,h) and surface of rotation (Fig 2). All of these flexible bodies can be modeled with the rectangular or isosceles trapezoid finite elements.

CONCLUSIONS
Equations of motion for a finite element in accelerations of a node-fixed system of coordinates derived in Sect. 2 allow simulation of flexible bodies with large gross motion, small strains and large flexible deflections. This approach is a kind of absolute nodal coordinate formulation. It is characterized by a clear structure of equations, which is important for computation of Jacobian matrices of elastic and damping forces required by an implicit solver. The number of floating point operations for the evaluation of equation terms is relatively small. The process of evaluation of equations of motion during the simulation can be easily parallelized.
The method of generation of equations of motion proposed in this paper can be applied to simulation of one or several flexible tires in multibody models of automotive vehicles as well as to analysis of air spring dynamics taking into account flexible bellows and a large change of the spring height. Of course, the finite element model should be modified to take into account the tire and bellows inflation pressure as well as the anisotropic properties of the element, which will be done in future research. Along with these objects, the approach can be applied for modeling many other flexible elements of transport systems subjected to large deflections such as leaf springs, rope systems, conveyor belts, and so on.
The trapezoidal finite element in absolute coordinates for… 151.