Development of ANCF tetrahedral finite elements for the nonlinear dynamics of flexible structures

In this paper, methods for developing isoparametric tetrahedral finite elements (FE) based on the absolute nodal coordinate formulation (ANCF) are presented. The proposed ANCF tetrahedral elements have twelve coordinates per node that include three position and nine gradient coordinates. The fundamental differences between the coordinate parametrizations used for conventional finite elements and the coordinate parametrizations employed for the proposed ANCF tetrahedral elements are discussed. Two different parametric definitions are introduced: a volume parametrization based on coordinate lines along the sides of the tetrahedral element in the straight (un-deformed) configuration and a Cartesian parametrization based on coordinate lines directed along the global axes. The volume parametrization facilitates the development of a concise set of shape functions in a closed form, and the Cartesian parametrization serves as a unique standard for the element assembly. A linear mapping based on the Bezier geometry is used to systematically define the cubic position fields of ANCF tetrahedral elements: the complete polynomial-based eight-node mixed-coordinate and the incomplete polynomial-based four-node ANCF tetrahedral elements. An element transformation matrix that defines the relationship between the volume and Cartesian parametrizations is developed and used to convert the parametric gradients to structure gradients, thereby allowing for the use of a standard FE assembly procedure. A general computational approach is employed to formulate the generalized inertia, external, and elastic forces. The performance of the proposed ANCF tetrahedral elements is evaluated by comparison with the conventional linear and quadratic tetrahedral elements and also with the ANCF brick element. In the case of small deformations, the numerical results obtained show that all the tetrahedral elements considered can correctly produce rigid body motion. In the case of large deformations, on the other hand, the solutions of all the elements considered are in good agreement, provided that appropriate mesh sizes are used.

form, and the Cartesian parametrization serves as a unique standard for the element assembly. A linear mapping based on the Bezier geometry is used to systematically define the cubic position fields of ANCF tetrahedral elements: the complete polynomialbased eight-node mixed-coordinate and the incomplete polynomial-based four-node ANCF tetrahedral elements. An element transformation matrix that defines the relationship between the volume and Cartesian parametrizations is developed and used to convert the parametric gradients to structure gradients, thereby allowing for the use of a standard FE assembly procedure. A general computational approach is employed to formulate the generalized inertia, external, and elastic forces. The performance of the proposed ANCF tetrahedral elements is evaluated by comparison with the conventional linear and quadratic tetrahedral elements and also with the ANCF brick element. In the case of small deformations, the numerical results obtained show that all the tetrahedral elements considered can correctly produce rigid body motion. In the case of large deformations, on the other hand, the solutions of all the elements considered are in good agreement, provided that appropriate mesh sizes are used. Keywords Flexible multibody system dynamics · Absolute nodal coordinate formulation · Bezier tetrahedral patch · Isoparametric ANCF tetrahedral finite elements · Position and gradient continuity conditions

Introduction
The goal of this investigation is the development of ANCF tetrahedral finite elements (FE) for the nonlinear dynamics and vibration analysis of multibody system (MBS) applications. The development of these new elements is necessary in order to correctly capture large finite rotations and large deformations of flexible system components with complex geometries. In the past three decades, MBS dynamics has emerged as an independent research field which has its own challenges including modeling complex engineering and physics systems that consist of interconnected bodies and joints. The bodies and the joints in a MBS application can be rigid and/or flexible [1][2][3][4][5][6][7][8][9][10][11]. For such complex systems, it is necessary to adopt a general analysis approach which can correctly capture geometric nonlinearities and allows for efficient formulation of the algebraic constraint equations that define mechanical joints. The absolute nodal coordinate formulation (ANCF), which can be used with a non-incremental solution procedure, can correctly describe finite rotations and large deformations, and therefore, it is suited for implementation in general MBS algorithms [12][13][14]. In recent years, ANCF elements have been validated numerically and experimentally and successfully used in modeling complex applications [15][16][17][18][19][20][21][22][23][24][25][26]. ANCF elements have a constant mass matrix, a constant generalized gravity force vector, and zero Coriolis and centrifugal forces regardless of the magnitude of displacement. Because the ANCF mass matrix is constant, a Cholesky transformation can be used at a preprocessing stage to obtain an identity inertia matrix associated with the ANCF Cholesky coordinates, leading to an optimum sparse matrix structure of the MBS equations of motion [27].
In structural mechanics, conventional tetrahedral elements are widely used because of their simplicity, versatility, and ability to describe complex geometries. Furthermore, the availability of several tetrahedral automatic meshing algorithms, including algorithms based on the Delaunay triangulation method, allows for automatic construction of complex tetrahedral meshes that can represent almost any threedimensional geometry with high fidelity. In the FE literature, three types of tetrahedral elements that are of interest in this investigation can be distinguished according to the order of the interpolation used in the definition of the element shape functions: the four-node linear tetrahedral element (FNL), the ten-node quadratic (TNQ) tetrahedral element, and the twentynode cubic (TNC) tetrahedral element. An interesting property of the tetrahedral elements is the fact that their kinematic representation allows for the use of a complete set of monomials for the polynomial expansion of the basis functions. This feature makes the number of the monomial terms used in the polynomial expansion based on the Pascal tetrahedron always the same as the number of the tetrahedral nodes [28]. That is, for a given order of interpolation, a linear mapping between a Bezier tetrahedral patch basis functions and the shape functions of the tetrahedral element always exists [29]. Therefore, for a given order of interpolation, the shape functions of the tetrahedral elements can be directly constructed starting with a Bezier tetrahedral patch using a set of volume coordinates. This simple idea is one of the key geometry concepts employed in this investigation to develop two new ANCF tetrahedral elements.
In the last two decades, while a large amount of research has been devoted to the development of new ANCF beams, plate, and shell elements, very few ANCF hexahedral and tetrahedral elements have been proposed. Because the conventional hexahedral and tetrahedral elements do not use position vector gradients as nodal coordinates, the continuity of the gradient and rotation fields at the nodal points is not ensured. On the other hand, ANCF elements ensure the continuity of the rotation, strain, and stress fields [30]. While ANCF solid (brick) elements were proposed [31,32], only recently Olshevskiy et al. [32] proposed new cubic ANCF hexahedral and tetrahedral elements. While the ANCF brick element proposed by Olshevskiy et al. [32] uses a complete set of position vector gradients as nodal coordinates and ensures the continuity of the gradients at the nodes, Olshevskiy et al. did not address the development of the ANCF tetrahedral element in their paper. Subsequently, two eight-node ANCF brick elements were used in order to simulate the liquid sloshing phenomenon by Wei et al. [33]. The first ANCF brick element is obtained using an incomplete cubic polynomial representation, while the second element is obtained directly from a cubic B-spline volume representation. More recently, a new ANCF tetrahedral element was proposed by Mohamed [34] who used Cartesian gradients from the outset to define a cubic set of polynomials. Nonetheless, efficient quadrature rules based on the Gauss integration method for tetra-hedral elements were not employed for this type of ANCF tetrahedral element, and the volume integration of the system mass matrix, the vector of generalized elastic forces, and the generalized gravity force vector were performed using conventional numerical quadrature procedures based on Newton-Cotes formulas, such as the trapezoidal rule or the Simpson method [35]. For these reasons, the ANCF tetrahedral element proposed by Mohamed [34] requires expensive computations for the evaluation of the generalized elastic forces and exhibits a slow convergence rate.
The goal of this investigation is to address geometry challenges associated with the description of the tetrahedral element kinematics and propose approaches for developing ANCF tetrahedral elements that differ from the ones presented in the literature. As shown in this paper, complete cubic tetrahedral polynomials have sixty coefficients, or equivalently twenty Bezier control points. Therefore, an ANCF tetrahedral element based on a complete polynomial representation requires the use of mixed coordinates, some of which can be located on the element faces, or the use of curvature coordinates in addition to the position and gradient coordinates. Nonetheless, a lowerorder element that employs only position and gradient coordinates at the tetrahedral nodes can be systematically developed by imposing set of algebraic equations that automatically reduce the element dimension [36]. This systematic approach is used in this investigation to develop a new lower-order ANCF tetrahedral element from another new higher-order ANCF tetrahedral element obtained using Bezier geometry. Specifically, the use of the approaches proposed in this study is demonstrated by developing two threedimensional ANCF tetrahedral elements, an eightnode mixed-coordinate (ENMC) and a four-node (FN) ANCF tetrahedral elements. The new ANCF tetrahedral elements employ as nodal coordinates global nodal positions and gradient coordinates and have displacement fields that ensure the gradient and rotation field continuity and are consistent with the Bezier geometry widely used to construct tetrahedral patches in computer-aided design (CAD) programs. Therefore, this investigation also contributes to the integration of computer-aided design and analysis (I-CAD-A) and is an important step in the development of new mechanics-based CAD/analysis system that will allow the use of a unified geometry and analysis mesh from the outset. The performance of the new ANCF tetrahedral elements is evaluated by comparing their solutions with the solutions obtained using the conventional linear and quadratic tetrahedral elements and with the solutions obtained using the ANCF brick element implemented in the general-purpose MBS program SIGMA/SAMS (Systematic Integration of Geometric Modeling and Analysis for the Simulation of Articulated Mechanical Systems). As shown in the paper, if appropriate mesh sizes are used, a good agreement is obtained, in general, between the numerical results of the conventional tetrahedral, the ANCF brick, and the proposed ANCF tetrahedral element models.
The algebraic constraint equation method used in this investigation to obtain a lower-order tetrahedral element can be used in the design of new ANCF elements that require the use of incomplete polynomials. This method can be applied from the outset to convert complete polynomials to incomplete polynomials that have symmetric structure and suited for ANCF elements with certain number and type of nodal coordinates that require the use of smaller number of basis functions. In order to avoid trials and errors in identifying such incomplete polynomials, the method of algebraic constraint equations used in developing the FN element presented in this investigation can be systematically used from the outset by imposing the algebraic constraint equations to reduce the number of each complete cubic polynomial basis functions from 20 to 16 and to systematically define the incomplete polynomial as discussed in the appendix of this paper.
This paper is organized as follows: In Sect. 2, the tetrahedral element parametrization used in this investigation is introduced. In particular, the volume gradient vectors used to obtain compact and closed form expressions of the element shape functions are introduced. The Cartesian gradients used to assemble the new ANCF tetrahedral elements are also defined in this section. In Sect. 3, the transformation between the volume and Cartesian gradients is defined. In Sect. 4, the position field of the new ANCF/ENMC tetrahedral elements is obtained, while in Sect. 5, the lowerorder ANCF/FN element shape functions are derived by imposing constraints on some nodes of the element. The isoparametric properties of the two new ANCF finite elements are discussed in Sect. 6. In Sect. 7, the procedure for developing the element equations of motion is described. In Sect. 8, numerical results obtained using different models based on different tetrahedral elements are analyzed in order to perform a comparative study and evaluate the dynamic behavior and convergence characteristics of the proposed elements. In Sect. 9, summary and conclusions drawn from this investigation are provided.

Tetrahedral element parametrization
In general, the position vector gradients are determined by differentiation of the position vector with respect to coordinate lines or parameters. Therefore, the gradient vectors are tangent to these coordinate lines which are discussed in this section when tetrahedral elements are used. As shown in this section, proper understanding of the element parametrization is necessary in order to develop the new ANCF tetrahedral elements which require the use of parametrization that conceptually differs from the parametrization used in the conventional FE literature. In order to understand the fundamental differences between the parametrization used in the development of the new ANCF tetrahedral elements and the parametrization used in the conventional FE literature, the latter is first reviewed in this section. The position field of a three-dimensional continuum body r = r 1 r 2 r 3 T can always be written in terms of an assumed set of coordinate lines defined by the Cartesian parameters x = x y z T as follows: In the case of tetrahedral elements, the use of the Cartesian parametrization is necessary when assembling elements which have gradient vectors as nodal coordinates. However, it is more convenient to use other parameters for describing the tetrahedral element geometry. Consider an arbitrary tetrahedron whose corner nodes are ordered counterclockwise following the right-hand rule and are labelled as 1, 2, 3, and 4, as shown in Fig. 1. The positions of the tetrahedron vertices are defined in the global reference system by the vectors v k = x k y k z k T , k = 1, 2, 3, 4, where x k , y k , and z k represent the Cartesian coordinates of the tetrahedron vertex k. The vector of Cartesian coordinates x is associated with the global reference system in which the position vectors of the tetrahedron corner nodes v k , k = 1, 2, 3, 4 are defined. The Cartesian parameter vector x can be expressed in terms of four dimensionless volume coordinates ξ = ξ η ζ χ T as The set of volume coordinates ξ, η, ζ , and χ , shown in Fig. 1, are not independent parameters and are related by the following algebraic equation: The volume coordinates whose sum is unity range from zero to one in the manner shown in Fig. 2 [29,37]. Because gradient vectors are used in the kinematics of the ANCF tetrahedral elements, it is important to understand the geometric meaning of the volume coordinates which form straight lines normal to the faces of the tetrahedral element, as shown in Fig. 2. On the other hand, for a general tetrahedral element, the contour planes or isoplanes are orthogonal to the directions of the coordinate lines. For example, the equation ξ = c represents a set of straight planes parallel to the tetrahedron face identified by the vertices 2, 3, and 4. Similar comments apply to η, ζ , and χ , as shown in Fig. 3. Relabelling the element vertices, for convenience, as A, B, C, and D, which correspond, respectively, to the tetrahedron vertices 1, 2, 3, and 4, as shown in Fig. 4, an arbitrary point P inside the tetrahedron ABC D divides the tetrahedron into four sub-tetrahedrons P BC D, PC D A, P D AB, and P ABC. Let be the volume of the tetrahedral element, and 1 , 2 , 3 , and 4 are the volumes of the four sub-tetrahedrons P BC D, PC D A, P D AB, and P ABC, respectively. It follows that the volume coordinates can be written as ξ = 1 / , the volume coordinates ξ = ξ η ζ χ T can be written as x = ξ v 1 + ηv 2 + ζ v 3 + χ v 4 , or equivalently: where x k , y k , and z k represent the Cartesian coordinates of the tetrahedron vertex k. These equations define a direct relationship between the three independent Cartesian coordinates x, y, and z and the four dependent volume coordinates ξ , η, ζ , and χ . From these equations, it is clear that the four corner nodes (x 1 , y 1 , z 1 ), (x 2 , y 2 , z 2 ), (x 3 , y 3 , z 3 ), and (x 4 , y 4 , z 4 ) of the tetrahedral element correspond to the volume coordinates (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), and (0, 0, 0, 1), respectively. Combining the preceding two equations, one has wherex = 1 x y z T and the vertex matrix V is defined as: The inverse representation is ξ = V −1x , where where V = denotes the volume of the tetrahedral element which can be simply obtained from the determinant of the vertex matrix V as V = |V|/6 and the constants , and L 4,3 are defined in Appendix A of the paper. The equation ξ = V −1x can be written explicitly as Or equivalently as It follows that This equation can be rewritten in a matrix form as follows: where N is a 4 × 3 matrix defined as: It can be shown that, in the straight (un-deformed) configuration, the vectors n 1 = L 1,1 L 1,2 L 1,3 T , and n 4 = L 4,1 L 4,2 L 4,3 T are vectors normal to the tetrahedral faces opposite to the tetrahedral vertices 1, 2, 3, and 4, respectively. It is important to note that the coordinate lines associated with the dependent volume parameters ξ, η, ζ , and χ are normal to the faces of the tetrahedron as it is clear from Eq. (9). For example, using Eq. (9) one can write dξ = L 1,1 dx + L 1,2 dy + L 1,3 dz /(6V ) = n T 1 dx/(6V ), which shows that dξ is the projection of the Cartesian vector dx = dx dy dz T along the vector T normal to the face of the tetrahedron defined by the vertices 2, 3, and 4. Similar comments apply to η, ζ , and χ since from Eq. (9), dξ = Ndx/(6V ) which shows that dξ can be obtained from the product of the matrix N and dx. It is important to note that the geometric interpretation of the directions of the coordinate lines does not conflict with the equation (9) is obtained from Eq. (5). Therefore, an arbitrary change in the volume coordinates dξ consistent with the algebraic equation ξ + η + ζ + χ = 1 yields an arbitrary change in the Cartesian coordinates dx according to Because the vectors n 1 , n 2 , n 3 , and n 4 that form the matrix N are normal to the faces of the tetrahedral element, the volume coordinate lines are normal to the faces of the tetrahedral element, as previously mentioned. It follows that the position vector gradients r ξ , r η , r ζ , and r χ are tangent to the coordinate lines ξ, η, ζ , and χ , and there-fore, these gradient vectors are normal to the tetrahedral faces if ξ, η, ζ , and χ remain dependent. In order to be able to assemble the elements with gradient vectors as nodal coordinates, special definitions of the coordinate lines are employed for the corner nodes of the ANCF tetrahedral elements, as discussed in the following section.

Position vector gradient transformations
In order to create an ANCF tetrahedral mesh using a standard FE assembly procedure, a particular set of parametric gradients at the nodes defined by differentiation with respect to the volume coordinates is used. For a given node of the tetrahedral element, these special volume gradients define directions parallel to the tetrahedral sides and not normal to the tetrahedral faces as the ones used in the FE literature. To this end, a different definition of the tetrahedron coordinate lines is used for each node of the tetrahedral element by properly defining the dependent volume coordinate using the constraint equation ξ + η + ζ + χ = 1. Consequently, the geometric interpretation of the coordinate lines associated with the volume coordinates ξ, η, ζ , and χ , previously introduced in the preceding section, changes according to the particular parametrization used for each of the four nodes.

Node parametrization
Consider the first node of the tetrahedral element, denoted as node 1, shown in Fig. 1. The first volume coordinate ξ can be expressed as a linear function of the other three independent volume coordinates η i , ζ i , and where the subscript i is used to refer to independent volume parameters. This equation can be used to eliminate ξ from Eq. (4) and write Or equivalently: and v 2, It is important to note that, when the first homogeneous coordinate ξ is eliminated using the constraint equation between the tetrahedral volume coordinates, the geometric interpretation of the coordinate lines associated with the independent volume parameters η i , ζ i , and χ i changes accordingly, as shown in Fig. 5a, and this is also clear from the equation In this case the coordinate lines associated with the volume coordinate η i are a set of straight lines parallel to the tetrahedron side defined by the vertices 2 and 1 (2-1 side); the coordinate lines associated with the volume coordinate ζ i are a set of straight lines parallel to the tetrahedron 3-1 side, whereas the coordinate lines associated with the volume coordinate χ i are a set of straight lines parallel to the 4-1 side. A change in the Cartesian parameter vector x, therefore, is the result of a change in the η i , ζ i , and χ i parameters along, respectively, the tetrahedron sides defined by the vectors v 2,1 , v 3,1 , and v 4,1 which can also be made unit vectors. Consequently, in this case, the parametric gradients r η i , r ζ i , and r χ i are vectors tangent to the tetrahedron sides 1-2, 1-3, and 1-4, respectively. Using this definition of the coordinate lines for the first node of the ANCF tetrahedron element, a linear transformation between the volume gradient vectors r η i , r ζ i , and r χ i defined with respect to the independent set of volume parameters and the Cartesian nodal gradient vectors r x , r y , and r z can be written as This equation can be rewritten, by adding an identity transformation for the position vector of the first node of the ANCF tetrahedron element, as: For the first node of the ANCF tetrahedral element, the nodal coordinate vector e 1 , with gradients defined by differentiation with respect to the independent volume gradients, can be written in terms of another nodal coordinate vector , with gradients defined by differentiation with respect to the Cartesian parameters, as e 1 = T 1 p 1 , where the transformation matrix T 1 is defined as where I is the 3 × 3 identity matrix. The nodal coordinates expressed using the Cartesian parametrization allow for using a standard FE assembly procedure. As previously mentioned, using the independent volume parametrization, one can write This equation shows that a change in the Cartesian coordinate x is the result of changes in the independent volume coordinatesη i ,ζ i , andχ i along, respectively, the three unit vectorŝ v 2,1 ,v 3,1 , andv 4,1 that are not necessarily orthogonal unit vectors. Furthermore, since in the straight (un-deformed) configuration the Cartesian gradients can be written as r x = 1 0 0 T , r y = 0 1 0 T , and r z = 0 0 1 T , one can show that rη i =v 2,1 , 1 , which shows that, in the straight (un-deformed) configuration, the volume gradients, when independent volume parameters are used, are tangent to the three sides of the tetrahedral element that intersect at node 1, as shown in Fig. 5a. It is important to note that, since the independent volume parameters defined in this section are associated with straight lines tangent to the sides of the tetrahedral element, they have different geometric meaning when compared to the redundant volume parameters used in the FE literature. Furthermore, following the same procedure for nodes 2, 3, and 4, by respectively, one can define a set of transformation matrices at nodes 2, 3, and 4, as shown in Fig. 5b-d, given by: where The transformations obtained in this section allow for expressing the independent volume gradient vectors in terms of the Cartesian gradient vectors. Such transformations will facilitate the derivation of the shape functions for the ANCF tetrahedral element in closed and compact forms. These transformations also make the physical meaning of the gradient vectors used more obvious.

Element transformation matrix
Using the four transformation matrices developed in the previous subsection, the total vector of nodal coordinates for the ANCF tetrahedral element is defined as e = e 1 T e 2 T e 3 T e 4 T T . The vector of nodal coordinates e includes volume gradients and can be expressed in terms of the vector of element coordinates p = p 1 T p 2 T p 3 T p 4 T T which includes the Cartesian gradients. The vector of element coordinates e based on the volume gradients and the vector of element coordinates p based on the Cartesian gradients are related by a linear transformation e = Tp in which the transformation matrix T is a 48×48 block diagonal matrix that is defined as The use of the coordinate vector p allows for conveniently assembling the tetrahedral element mesh using a standard FE assembly procedure. In the following section, the cubic Bezier tetrahedron geometry is used to develop the shape functions of the ANCF eight-node mixed-coordinate (ENMC) tetrahedral element. Using the ENMC tetrahedral element shape functions, the shape functions of the fournode (FN) ANCF tetrahedral element can be derived by employing a set of linear constraint equations on the global position of the tetrahedral center nodes. It can be shown that both the proposed ENMC and FN tetrahedral elements are isoparametric elements and lead to an exact representation of the rigid body motion.

Eight-node mixed-coordinate (ENMC) element
In order to develop the shape functions of the ENMC tetrahedral element shown in Fig. 6a, a special type of Bezier tetrahedron which is defined using a cubic interpolation and a set of control points is considered. The global position vector of an arbitrary point located on the cubic Bezier tetrahedral patch can be written as r = m k=1 g k q k , wherem = 20, g k are the basis or blending functions of the cubic Bezier tetrahedron defined in terms of the set of volume coordinates ξ, η, ζ , and χ , and q k are the control points associated with the control net of the Bezier tetrahedron. The polynomial terms used for the interpolation of the cubic Bezier tetrahedral patch can be directly defined using the terms in the expansion of (ξ + η + ζ + χ ) 3 . The Bezier basis functions used in this investigation are The polynomial interpolation of the Bezier tetrahedral patch can be used as the assumed displacement field for the ENMC tetrahedral element, and the corresponding set of shape functions can be obtained using a linear transformation of coordinates. The relationship between the Bezier control points and the ENMC element nodal coordinates can be written as 27 q 1 + q 2 + q 3 + 3q 5 + 3q 6 + 3q 7 + 3q 8 + 3q 9 + 3q 10 + 6q 17 r 6 = 1 27 q 1 + q 2 + q 4 + 3q 5 + 3q 6 + 3q 11 + 3q 12 + 3q 13 + 3q 14 + 6q 18 r 7 = 1 27 q 2 + q 3 + q 4 + 3q 7 + 3q 8 + 3q 13 + 3q 14 + 3q 15 + 3q 16 + 6q 19 r 8 = 1 27 q 1 + q 3 + q 4 + 3q 9 +3q 10 + 3q 11 + 3q 12 + 3q 15 + 3q 16 + 6q 20 Using these linear algebraic relationships, the cubic Bezier tetrahedral patch can be converted into the ANCF/ENMC tetrahedral element position field which can be expressed as r =Sē, whereS is the shape function matrix of the ENMC tetrahedral element andē is the vector of the element nodal coordinates given, respectively, bȳ and where I is the 3 × 3 identity matrix and r ξ = ∂r/∂ξ , r η = ∂r/∂η, r ζ = ∂r/∂ζ , and r χ = ∂r/∂χ are the position vector gradients defined by differentiation with respect to independent volume parameters along the element sides as previously discussed. The ENMC shape functionss k , k = 1, 2, . . . , 20, ares The element shape function matrixS contains a set of interpolating functions which allow for expressing the position of an arbitrary point on the ANCF tetrahedral element as a linear combination of the element nodal coordinatesē. In the definition of the nodal coordinate vectorē, superscripts 1, 2, 3, and 4 refer to the element corner nodes ordered counterclockwise and following the right-hand rule, while superscripts 5, 6, 7, and 8 refer to the center nodes of the ENMC tetrahedral element as shown in Fig. 6a. The ENMC nodal coordinates are defined as In this equation, the computation of the partial derivatives that define the volume gradients is performed by substituting the constraint equations between the volume coordinates in the displacement field of the cubic Bezier tetrahedral patch. Since in the nodal coordinate vectorē, absolute positions and volume gradients at the nodes are used as generalized coordinates, an exact representation of rigid body motion is achieved by the shape functions of the ENMC tetrahedral element. Therefore, the proposed ENMC tetrahedral element is isoparametric and describes correctly rigid body motion. The relationship between the ENMC shape functions and Bezier basis functions is presented in Appendix A of this paper. The ENMC tetrahedral element has 60 nodal coordinates which consist of 24 position coordinates and 36 gradient coordinates. As explained in the preceding section, a linear transformation can be used for the ENMC tetrahedral element to write the coordinate vectorē and the element vector of structural coordinatesp that contain the Cartesian gradients r x , r y , and r z asē =Tp, where the transformation matrixT is a 60 × 60 block diagonal matrix defined as In this equation, the nodal transformation matrices are defined asT 1 = T 1 ,T 2 = T 2 ,T 3 = T 3 ,T 4 = T 4 , T 5 = I,T 6 = I,T 7 = I,T 8 = I, where I is the 3 × 3 identity matrix, as shown in the preceding section. It is important to note that, when the Cartesian gradients are used in a way consistent with the tetrahedral element geometry as explained in the paper, the shape functions of the ENMC tetrahedral finite element can be obtained in a compact closed form.

Four-node (FN) tetrahedral element
The development presented in the preceding section demonstrates that complete cubic tetrahedral polynomials have sixty coefficients, or equivalently twenty control points. Therefore, an ANCF tetrahedral element having four corner nodes and based on a complete polynomial representation requires the use of mixed coordinates, some of which can be located on the element faces, or the use of curvature coordinates in addition to the position and gradient coordinates. Nonetheless, a lower-order element that employs only position and gradient coordinates at the tetrahedral nodes can be systematically developed by imposing set of algebraic equations that automatically reduce the element dimension. This is the procedure used in this section to develop the new lower-order ANCF/FN element proposed in this investigation.
The shape functions of the FN tetrahedral element shown in Fig. 6b can be obtained from the shape functions of the ENMC tetrahedral element by imposing a set of linear constraint equations between the center nodes on the tetrahedral faces and a set of material points distributed on the sides of the tetrahedral element. These linear constraint equations can be written as r 5 = 1 9 r 1 5 + r 2 5 + r 3 5 + r 4 5 + r 5 5 + r 6 5 + r 7 5 + r 8 5 + r 9 5 r 6 = 1 9 r 1 6 + r 2 6 + r 3 6 + r 4 6 + r 5 6 + r 6 6 + r 7 6 + r 8 6 + r 9 6 r 7 = 1 9 r 1 7 + r 2 7 + r 3 7 + r 4 7 + r 5 7 + r 6 7 + r 7 7 + r 8 7 + r 9 7 r 8 = 1 9 r 1 8 + r 2 8 + r 3 8 + r 4 8 + r 5 8 + r 6 8 + r 7 8 + r 8 8 where r h k , k = 5, 6, 7, 8, h = 1, 2, . . . , 9 are global positions of material points located on the sides of the ENMC tetrahedral element and are defined as Using these linear algebraic constraint equations for the center nodes of the ENMC tetrahedral element, a linear transformation can be obtained and used to define the shape functions of the FN tetrahedral element. The position field of the FN tetrahedral element can then be written as r = Se, where S is the FN element shape function matrix and e is the vector of the element nodal coordinates given, respectively, by: and where I is the 3 × 3 identity matrix, r ξ = ∂r/∂ξ , r η = ∂r/∂η, r ζ = ∂r/∂ζ , and r χ = ∂r/∂χ are the position vector gradients defined using the volume parametrization in which the coordinate lines are parallel to the tetrahedral element sides. The shape functions s k , k = 1, 2, . . . , 16, of the FN tetrahedral element are defined in a closed form as s 1 = ξ ξ 2 + 3ξ (η + ζ + χ ) + 2 (ηζ + χ (η + ζ )) , The FN tetrahedral element has 48 nodal coordinates which consist of 12 position coordinates and 36 gradient coordinates. The vector of element coordinates e and the vector of structural coordinates p which has the Cartesian gradients r x , r y , and r z are related by a linear transformation e = Tp in which the transformation matrix T is previously defined in the paper. The shape functions of the FN and ENMC elements associated with the Cartesian gradients are presented in Appendix A of the paper. The procedure described in this section can be used from the outset to obtain incomplete polynomials that have symmetric structure and designed for new ANCF elements with certain number and type of nodal coordinates that require the use of incomplete polynomials. In order to avoid trials and errors in identifying such incomplete polynomials, the method of algebraic constraint equations used in developing the FN element presented in this section can be systematically employed from the outset by imposing the algebraic constraint equations to reduce the number of each complete cubic polynomial basis functions from 20 to 16 and to systematically define an incomplete polynomial symmetric in the Cartesian coordinates x, y, and z as discussed in the appendix of this paper.

Isoparametric property
It can be shown that the ENMC and FN elements proposed in this investigation are isoparametric elements that correctly describe rigid body motion and can be effectively used in a non-incremental solution framework for the numerical solution of the equations of motion. In general, a finite element can be considered as isoparametric if the same shape functions can be used to describe the element geometry and displacement. Therefore, the shape functions of an isoparametric finite element can be used to represent the geometry of the element in the straight (un-deformed) configuration, in the reference (curved) configuration, and in the current (deformed) configuration, as shown in Fig. 7. The current configuration of a continuum body is a general deformed configuration in which the element position field can be written as r = Se = STp, where e and p are, respectively, the vectors of volume and Cartesian nodal coordinates which define the element current configuration. In the reference configuration, the element position field can be written as r 0 = Se 0 = STp 0 , where e 0 and p 0 are, respectively, appropriate vectors of volume and Cartesian nodal coordinates which define the stress-free reference configuration and allow for describing a general curved geometry. In the straight (un-deformed) configuration, one can write x = Se s = STp s , where e s and p s are, respectively, vectors of volume and Cartesian nodal coordinates that define the geometry of the ele- In these equations, x k , y k , and z k represent the Cartesian coordinates of the tetrahedron vertex k associated with the element straight (un-deformed) configuration.
Using the shape functions of the FN tetrahedral element, the expression of the transformation matrix T, and considering the constraint between the tetrahedral coordinates ξ + η + ζ + χ = 1, one can show that This equation demonstrates that the straight (undeformed) configuration of a tetrahedral element can be captured by the shape functions of the FN tetrahedral element.

Element equations of motion
The equations of motion of the ANCF tetrahedral elements are derived in this investigation using the d'Alembert-Lagrange principle of virtual work based on a total Lagrangian formulation that allows for the use of a non-incremental solution procedure [14,38].

Element kinematics
The displacement field of the two ANCF tetrahedral elements developed in this investigation can be written as r = Se. The nodal coordinate vector e includes volume gradients which define position vector gradients parallel to the element sides. While the use of these volume gradients facilitates obtaining closed form expressions of the tetrahedral element shape functions, these volume gradients are not suitable for using a standard FE assembly procedure. In order to overcome this problem and be able to impose the connectivity conditions between different tetrahedral elements using simple linear constraint equations, the Cartesian or body parametrization of the gradient coordinates can be used using the linear coordinate transformation e = Tp, as previously discussed in this paper. Because the transformation matrix T is constant, one has δe = Tδp and δr = STδp. The virtual change in the Cartesian gradients of the FN tetrahedral element can be written as or equivalently δr x = S x Tδp, δr y = S y Tδp, and δr z = S z Tδp, where In these equations, S x , S y , and S z are the spatial derivatives of the shape function matrix with respect to the Cartesian coordinates x, y, and z, respectively. The evaluation of the Cartesian derivatives of the shape function matrix S x , S y , and S z of the tetrahedral element is necessary in the formulation of the element generalized elastic forces.

Element mass matrix
The virtual work of the inertia forces of the ANCF elements can be written as δW i = − V ρr T δr |J 0 | dV , where ρ is the mass density in the reference configuration,r is the absolute acceleration vector of a material point on the element, |J 0 | is the determinant of the matrix of position vector gradients in the reference configuration obtained by differentiation with respect to the coordinate lines of the straight configuration, and V is the element volume in the straight (un-deformed) configuration. Using the element kinematic equations, the virtual work of the inertia forces can be rewritten as where the element mass matrix M = T T V ρS T S |J 0 | dV ) T is constant, symmetric, and positive definite regardless of the amount of element displacement and rotation. It follows that the centrifugal and Coriolis inertial forces are identically zero. By using the Cartesian gradients in the vector p, a standard FE assembly can be used to define the body mass matrix.

Element external forces
The virtual work of the external forces δW e = V f T e δr |J 0 | dV , where f e is vector of distributed external forces such as the gravity and magnetic forces, can be written as where Q e = T T V S T f e |J 0 | dV is the generalized external force vector associated with the element nodal coordinates. By using the Cartesian gradients as nodal coordinates in the coordinate vector p, a standard FE assembly can be used to obtain the body generalized external forces.

Element elastic forces
The virtual work of the element elastic forces is δW s = − V σ : δε |J 0 | dV , where σ is the symmetric second Piola-Kirchhoff stress tensor, ε is the symmetric Green-Lagrange strain tensor, and the symbol : refers to tensor double contraction. In the formulation of the elastic forces, general constitutive models can be used with the ANCF tetrahedral elements developed in this investigation. Using the element kinematic equations, the virtual work of the elastic forces can be rewritten as where Q s = −T T V σ : ∂ε ∂e T |J 0 | dV is vector of the element generalized elastic forces. Because this vector is a highly nonlinear function of the nodal coordinates, the integral must be calculated using a numerical integration procedure based on the Gauss quadrature rules. By using the Cartesian gradients as nodal coordinates in the vector p, a standard FE assembly procedure can be used to obtain the body elastic forces.

Equations of motion
The system equations of motion of the ANCF tetrahedral mesh can be obtained using the virtual work principle, which can be written as δW i + δW e + δW s = 0, where δW i is the virtual work of the element inertia forces, δW e is the virtual work of the element external forces, and δW s is the virtual work of the element elastic forces. Using this approach, the element equations of motion can be obtained as Mp = Q e + Q s . As previously mentioned, because Cartesian gradients r x , r y , and r z are used in the generalized coordinate vector p, the vectors and matrices that appear in this equation can be assembled using a standard FE assembly process in order to obtain the equations of motion of the ANCF tetrahedral element mesh.

Numerical results and discussion
In this section, several tetrahedral element models are developed in order to evaluate the performance and demonstrate the use of the proposed ANCF four-node (FN) and eight-node mixed-coordinate (ENMC) tetrahedral elements. The numerical results obtained using these cubic ANCF tetrahedral elements, which ensure continuity of the positions and gradients at the nodes, are compared with the results obtained using the conventional four-node linear (FNL) and the ten-node quadratic (TNQ) tetrahedral elements, which guarantee only the position continuity at the nodes. In order to perform a consistent comparative study, the conventional FNL and TNQ tetrahedral elements are implemented using a total Lagrangian non-incremental solution procedure [14,38]. Furthermore, for the purpose . This comparative study shows, for the most part, a good agreement for the displacements, while slight differences in the strains between the conventional tetrahedral and the ANCF tetrahedral element models are observed. The linear and quadratic conventional tetrahedral elements do not ensure continuity of the rotations, strains, and stresses at the nodal points. Furthermore, the curvatures within these elements are either constant or identically zero. Therefore, the conventional tetrahedral elements are not well suited for bending problems. The numerical results presented in this section are obtained using a MATLAB computer program. A continuum mechanics approach and nonlinear displacement-strain relationships are used to evaluate the generalized elastic forces, and as a result, geometric nonlinearities are taken into account. The material model assumed is the Saint Venant-Kirchhoff hyperelastic constitutive law. The generalized elastic forces are numerically calculated by using the standard Gauss quadrature procedure with full integration [37].

Rigid body motion
The rectangular pendulum example, shown in Fig. 8, is used to demonstrate that the proposed ANCF tetrahedral elements are capable of correctly describing rigid body motion. The pendulum used in this example is assumed to be initially horizontal and to fall under the effect of its own weight. The pendulum has length L =    inertia of I x x = 0.655 kg m 2 , I yy = 1.31 kg m 2 , and I zz = 0.655 kg m 2 . In this numerical example, a large modulus of elasticity is used in order to ensure small deformations, allowing for comparison with the results obtained using a rigid pendulum model. The symmetric mesh of the flexible pendulum is shown in Fig. 9. To construct a symmetric tetrahedral mesh for the flexible pendulum, two types of asymmetric macro-elements made of 5 tetrahedral elements are used, as shown in Fig. 10. Figure 11 shows the complete prismatic macroelement composed of 40 tetrahedral elements which is employed as a building block for creating the tetrahedral element mesh. Since the elastic modulus used in this test is large, only one prismatic module composed of 40 tetrahedral elements is used to create the FE model of the stiff pendulum. The flexible pendulum is constrained to the ground with three spherical joints located at points A 1 , A 2 , and A 3 , while the rigid pen-dulum is constrained to the ground with only one revolute joint located at point A 2 , as shown in Fig. 8. In this numerical example, the constraint equations are linear algebraic equations, systematically eliminated at a preprocessing stage leading to a set of ordinary differential equations of motion. In addition to the FN and ENMC tetrahedral element models, two other models based on the conventional FNL and TNQ tetrahedral elements are developed. The numerical integration method used to solve the resulting ordinary differential equations of motion is the Gear algorithm, which is a sixth-order implicit linear multistep scheme with a constant step size [35]. The time interval considered for the numerical simulations is T = 2.0 s, whereas the minimum time step used is t = 5 × 10 −6 s. Figure 12 shows the vertical displacement of the pendulum tip point B obtained using the stiff flexible and rigid body models. The results obtained using the stiff and rigid pendu- lum models are in good agreement, which demonstrates that the proposed ANCF tetrahedral elements correctly describe rigid body motion and yield consistent mass distribution.

Flexible body motion
The straight rectangular pendulum example, shown in Fig. 8, is used to assess the performance of the FN and ENMC tetrahedral elements in the case of large finite rotations and large deformations. In this numerical example, the model data are as previously presented except for the pendulum width, thickness, and modulus of elasticity which are, respectively, assumed to be H = 0.1 m, W = 0.1 m, and E = 2.1 × 10 6 N/m 2 . The time interval considered for the numerical simulations is T = 1.0 s, whereas the minimum time step used is t = 6.25 × 10 −6 s. The numerical integration method used to solve the resulting ordinary differential equations of motion is the Adams-Bashforth algorithm, which is a sixth-order explicit linear multistep scheme with a constant step size [35]. First, a convergence analysis and a comparative study with conventional tetrahedral elements are carried out. To this end, FE models based on the FNL, TNQ, ANCF/FN, and ANCF/ENMC tetrahedral elements are developed using a symmetric mesh made of prismatic macroelements as shown in Figs. 9, 10, and 11. If N m is the number of modules constructed using 40 tetrahedral elements and employed to create symmetric tetrahedral meshes, the total number of tetrahedral ele-

Summary and conclusions
In this investigation, two new ANCF tetrahedral elements: eight-node mixed-coordinate (ENMC) and fournode (FN) tetrahedral elements, are developed and their performance is numerically evaluated. Both tetrahedral elements are based on a cubic polynomial interpolation obtained using the linear relationship between the ANCF nodal coordinates and the Bezier control points that define a cubic Bezier tetrahedral patch. As explained in the paper, the ANCF tetrahedral elements developed are isoparametric elements, can consistently describe arbitrary rigid body motion, allow for repre-

Fig. 22
Normal strain ε zz at the flexible pendulum center (circle ANCF brick element, square ANCF/FN tetrahedral element, diamond ANCF/ENMC tetrahedral element senting large finite rotations, and impose no restrictions on the amount of deformation within the finite element. Unlike the conventional tetrahedral elements used in the FE literature, the proposed ANCF tetrahedral elements ensure the continuity of the gradients, rotations, strains, and stresses at the nodal points at the element corners. The position vector gradients based on the volume parametrization are used to obtain compact and closed form expressions of the element shape functions. These volume gradients are systematically converted into Cartesian position vector gradients in order to allow for the use of a standard FE assembly procedure. The proposed ANCF tetrahedral elements lead to a constant mass matrix and zero Coriolis and centrifugal inertia forces. The numerical results obtained using the proposed new FN and ENMC tetrahedral elements are compared with the numerical results obtained using the conventional four-node linear (FN) and tennode quadratic (TNQ) tetrahedral elements, as well as the results obtained using the ANCF brick element.
The results obtained in this study showed, in general, good convergence characteristics of the proposed ANCF tetrahedral elements. When compared to the conventional tetrahedral elements, the proposed ANCF elements do not require any special form of strain-and stress-averaging techniques. In particular, it is demonstrated that complex deformation shapes can be modeled using a smaller number of ANCF/ENMC tetrahedral elements when compared to the conventional linear and quadratic (FNL and TNQ) tetrahedral elements.
It is important also to point out that the method of algebraic constraint equations used to obtain the FN incomplete polynomials from the complete polynomials has several advantages that include eliminating the trial and errors in identifying the incomplete polynomials, ensuring that the obtained incomplete polynomials are symmetric in the coordinates, and making clear the geometric restrictions imposed in the process of reducing the number of the basis functions. s * The above FN shape functions can be written in terms of the ENMC shape functions as ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (A.5) To obtain the above equations, the following constant coefficients need to be defined: where x k , y k , and z k represent the Cartesian coordinates of the tetrahedron vertex k, and V = V 1 +V 2 +V 3 +V 4 . The tetrahedron volume coordinates ξ, η, ζ , and χ can be expressed as a linear combination of the tetrahedral Cartesian coordinates x, y, and z as: ξ = 1 6V 6V 1 + L 1,1 x + L 1,2 y + L 1,3 z η = 1 6V 6V 2 + L 2,1 x + L 2,2 y + L 2,3 z ζ = 1 6V 6V 3 + L 3,1 x + L 3,2 y + L 3,3 z χ = 1 6V 6V 4 + L 4,1 x + L 4,2 y + L 4,3 z It can also be shown that the ENMC and FN shape functions written in terms of the Cartesian gradientss * k and s * k can be expressed in terms of the ENMC and FN shape functions associated with the volume gradients s k and s k asS * =ST and S * = ST, where the transformation matricesT and T are previously defined in the paper.
In developing new ANCF finite elements with certain number and type of nodal coordinates, the use of incomplete polynomials may be necessary. In order to avoid trials and errors in identifying such incomplete polynomials and obtain symmetric structure in x, y, and z, the method of algebraic constraint equations used in developing the FN element presented in this paper can be systematically used. For the FN element, one can use from the outset, the algebraic constraint equations to reduce the number of each polynomial basis function from 20 to 16 and to systematically define the incomplete polynomial which has the following basis expressed in terms of the volume coordinates: h 1 = ξ ξ 2 + 4ηχ + 4ζ (η + χ ) , h 2 = η η 2 + 4ξχ + 4ζ (ξ + χ ) , h 3 = ζ ζ 2 + 4ξχ + 4η (ξ + χ) , h 4 = χ χ 2 + 4ηξ + 4ζ (η + ξ ) , ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (A.14) Using the relationship between the Cartesian and volume coordinates (Eq. A.13), the basis functions can be written in terms of the Cartesian coordinates x, y, and z as h i = 20 j=1 a i, j b j , where b 1 b 2 . . . b 20 = 1 x y z x 2 y 2 z 2 x y yz xz x 3 y 3 z 3 x 2 y y 2 x y 2 z z 2 y x 2 z z 2 x xyz (A. 15) and the coefficients a i, j , i = 1, 2, . . . , 16, j = 1, 2, . . . , 20, can be systematically defined. One can show that the incomplete polynomial defined by the basis functions h i = 20 j=1 a i, j b j has the linear terms 1, x, y, and z that ensure that the rigid body motion can be correctly described, and such a polynomial will have a symmetric structure in x, y, and z.