Use of independent volume parameters in the development of new large displacement ANCF triangular plate/shell elements

In this paper, a new unified kinematic description, obtained from Bezier geometry using linear mapping and position vector gradients associated with three independent parameters, is used to develop large displacement plate/shell finite elements (FE). Contrary to the conventional FE method, in the approach developed in this paper based on the absolute nodal coordinate formulation (ANCF), no distinction is made between plate and shell structures. The proposed ANCF triangular plate/shell elements have 12 coordinates per node: three position coordinates and nine position gradient coordinates that define vectors tangent to coordinate lines at the nodes. The fundamental differences between the conventional FE and the new ANCF parameterizations are highlighted. In this investigation, two different parameterizations, each of which employs independent coordinates, are used. In the first parameterization, called volume parameterization, coordinate lines along the sides of the triangular element in the straight (un-deformed) configuration are used in order to facilitate the development of closed-form cubic shape functions. In the second parameterization, called Cartesian parameterization, coordinate lines along the global axes of the structure (body) coordinate system are used to facilitate the element assembly. The element transformation between the volume and the Cartesian parameterizations is developed and used to define the structure inertia and elastic forces. Three new fully parameterized ANCF triangular plate/shell elements are developed in this investigation: a four-node mixed-coordinate element (FNMC) and two three-node elements (TN1 and TN2). All the elements developed in this investigation lead to a constant mass matrix and zero Coriolis and centrifugal forces. A non-incremental total Lagrangian procedure is used for the numerical solution of the nonlinear equations of motion. The performance of the proposed ANCF triangular plate/shell elements is analyzed by comparison with the ANCF rectangular plate element and conventional three-node linear (TNL) and six-node quadratic (SNQ) triangular plate elements.


Introduction
In this section, a brief literature review of conventional plate and shell finite elements, the features of ANCF plate/shell elements that distinguish them from conventional elements, and the goal, scope, and contributions of this investigation are discussed.

Literature review
In the conventional FE, structural mechanics, and vibration theory formulations; plates and shells are described using different geometric approaches and different displacement fields . Noor provided a comprehensive survey of books and monographs on plate and shell elements [1]. Reddy presented a unified treatment of plates and shell elements commonly used in aerospace, automotive, and civil engineering applications [2]. Armenakas et al. [3] used plate and shell elements to study the vibration of thin-walled structures subjected to severe dynamic loading conditions. Ashwell and Gallagher edited a series of conference papers on the FE theory of plate and shell elements and discussed the meshing algorithms for modeling curved geometry [4]. Axelrad used a local approach for the stability analysis of plate/shell large deformations [5]. Bieger used analytical methods based on series expansions to investigate circular cylindrical shells subjected to concentrated loads [6]. Billington studied complex three-dimensional shell structures [7], while Brush et al. [8] examined the buckling and stability of bar, plate, and shell elements using energy methods. Calladine investigated the coupling between the bending and stretching of shell structures under static loading [9]. Cox analyzed the buckling of plate and shell elements using continuum mechanics and energy methods [10]. Dym discussed important elasticity problems applicable to shell theory [11], while Flugge analyzed both plate and shell elements using modern analytical and computational methods [12]. Gould studied thin shell structures using a new FE approach for modeling curved geometric shapes [13][14][15]. Heyman investigated the static equilibrium of thin shell structures with complex curved geometry [16]. Hinton and Owen developed a computer program for the FE solution of static and dynamic problems of thin plate and shell structures considering elastoplastic and geometrically nonlinear anisotropic effects [17]. Huang devel-oped an FE software for the static and dynamic analyses of plate and shell structures using an approach designed for alleviating shear and membrane locking that occurs in degenerated Mindlin curved elements [18]. Hughes and Hinton discussed the plate and shell element technology for developing efficient computational algorithms [19,20]. Kelkar and Sewell examined fundamental analysis and design issues related to shell structures [21]. Kratzig and Onate edited a series of conference papers on the nonlinear dynamics and the computational mechanics of plate and shell structures [22]. Kraus studied the load-carrying capacity and analyzed the dynamic behavior of thin shell structures [23]. Kuhn examined the stresses of aircraft shell structures using analytical and computational methods [24]. Timoshenko and Woinowsky-Krieger made significant contributions to the development of plate/shell theory of elasticity [25]. Ugural studied the stresses in the case of small and large deformations of plate and shell structures [26]. Werner investigated the complex vibration behavior of plates, shells, curved membranes, rings, and thin structures by using both classical and computational methods [27].
In the classical FE theory, distinction is always made between plates and shells, and shell geometry cannot be systematically obtained from conventional plate element geometry because of the nature of the nodal coordinates and/or polynomial interpolations used. Use of infinitesimal rotations that impose restriction on the magnitude of rotations, employing vector transformations that differ from the transformations used for the gradients required to properly create accurate geometry, and/or treatment of the slopes as vectors and not strictly as gradients are among the reasons that arbitrary initially curved shell geometry cannot be systematically obtained from conventional plate elements. These geometric restrictions can be alleviated by using control points as in the case of computational geometry methods or by using position vector gradients with the proper gradient tensor transformation in order to correctly capture discontinuities.

Plate/shell elements
It is the objective of this study to introduce a new unified approach, based on the absolute nodal coordinate formulation (ANCF), for modeling three-dimensional triangular plate/shell structures using the same geometric representation and same displacement field. The proposed ANCF triangular plate/shell elements allow for correctly capturing finite rotations, large deformations, and complex geometry that cannot be easily captured using the conventional FE methods; some of which employ infinitesimal rotations as nodal coordinates. Capturing large displacements and complex geometry is necessary for accurate analysis of both structural and multibody system (MBS) applications [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43]. MBS computational challenges, in particular, can be addressed using ANCF elements that impose no restrictions on the amount of rotation or deformation within the element, lead to a constant mass matrix and zero centrifugal and Coriolis forces, correctly account for the dynamic coupling between the rigid-body motion and the elastic deformation, and allow for using a nonincremental procedure for the numerical solution of the equations of motion [44]. In recent years, several ANCF elements have been proposed, validated numerically and experimentally, and successfully used in modeling complex MBS applications [45][46][47][48][49][50][51][52][53][54][55][56][57]].

Conventional plate elements
Although a large number of quadrilateral plate and shell elements have been proposed in the conventional FE and ANCF literature [58][59][60], these rectangular elements are not suitable for meshing complex geometry with irregular shapes, and therefore, the development of triangular plate/shell elements becomes necessary. Conventional planar and spatial triangular plate elements are widely used for solving static and dynamic problems [61]. Furthermore, the availability of several automatic meshing algorithms facilitates the development of complex and high fidelity triangular meshes. Another important reason for the popularity of the triangular elements in the analysis of machines and structures is the possibility of performing mesh refinements in regions of stress concentrations that are of particular interest in the FE durability analysis. In structural mechanics, several types of triangular plate elements having transversal gradients as extensible directors are available [62]. While extensible directors are not often treated as position vector gradients throughout the solution algorithms, two conventional triangular plate elements are of particular interest in this investigation and can be distinguished according to their order of interpolation. These triangular plate elements are the three-node linear (TNL) element and the six-node quadratic (SNQ) element. The shape functions of the TNL and SNQ triangular plate elements are presented in the Appendix of the paper. An interesting property of the triangular plate elements that are based on transversal slopes is the fact that their kinematic description allows using a complete set of polynomial functions for the definition of the basis functions expressed in terms of dependent dimensionless volume coordinates with additional linear basis functions expressed in terms of dimensionless thickness coordinate. For a given order of interpolation, this important feature ensures the existence of a linear mapping between the Bezier triangular basis functions and the triangular element shape functions. In this investigation, three new ANCF fully parameterized plate/shell elements are developed using the simple geometric concepts introduced in this section.

ANCF plate elements
In recent years, several ANCF triangular plate elements have been proposed. Using a set of basis functions originally developed by Specht and Morley [63,64], Dmitrochenko and Mikkola proposed two new ANCF triangular plate elements by considering an incomplete set of cubic shape functions [65]. The first triangular plate element developed by Dmitrochenko and Mikkola is based on the basis functions proposed by Specht and employs constant coefficients which depend on the geometry of the particular triangular plate element considered [63]. This ANCF triangular plate element can effectively describe large deformations but it can be inefficient for large meshes because the pre-evaluation of the element shape functions at the Gauss quadrature nodes must be repeated for each element in the triangular mesh, resulting in an increased computational cost. The second triangular plate element developed by Dmitrochenko and Mikkola is based on the basis functions proposed by Morley and employs as nodal coordinates absolute displacement vectors and a particular set of slope vectors normal to the sides of the triangular plate element leading to a set of nonconformal shape functions [64]. Because this triangular plate element produces meshes that are non-conformal along the element sides, it can introduce significant errors in predicting large deformations [66]. A new fully parameterized ANCF triangular plate element was also proposed by Mohamed in order to capture the plate thickness deformation as well as the transversal shear deformation [67]. In the ANCF triangular plate element proposed by Mohamed, three additional shape functions were added to the nine shape functions of the ANCF thin triangular plate element proposed by Dmitrochenko and Mikkola [65]. Olshevskiy et al. [68] proposed two new types of three-and four-node planar plate elements obtained from spatial ANCF plate elements. Dmitrochenko and Pogorelov [69] proposed a way to generate new ANCF rectangular and triangular elements employing displacement fields of conventional elements used in structural mechanics. However, all the ANCF triangular plate elements cited above were developed using a kinematic representation based on the conventional triangular plate elements, and have complex shape function expressions. Furthermore, the shape functions of some of these elements are not reported in a compact closed-form and numerical methods were used for their evaluation [70].
In addition to the triangular plate elements obtained from the conventional elements, Chang et al. [71] proposed three new ANCF triangular shell elements based on Bezier triangular patches. These ANCF triangular elements were referred to as the thin triangular shell element, the low-order thick triangular shell element, and the high-order thick triangular shell element. Because, in this case, the basis functions are independent of the nodal Cartesian coordinates in the parametric domain, concise expressions for the shape functions are obtained compared with the triangular elements obtained from the conventional finite elements. However, in the kinematic description of the triangular plate elements proposed by Chang et al. [71], the geometric interpretation of the gradient vectors associated with the parametric coordinates is not obvious. Furthermore, only geometric continuity associated with the tangential direction of the element position field can be imposed at the element interfaces, and consequently, the structural continuity of the position vector gradients cannot be enforced leading to discontinuous strains and stresses at the mesh interface nodes. Recently, new types of planar ANCF triangular elements and spatial ANCF tetrahedral elements with closed-form shape function expressions were developed starting with Bezier basis functions [72,73]. Using a new approach, the geometric interpretation of the ANCF position vector gradients is utilized in order to facilitate the development of closed-form shape function expressions for these new elements.

Scope of this investigation
The main goal of this paper is to address challenges associated with the geometric description of the threedimensional triangular plate elements, propose new ANCF triangular plate/shell elements that differ from the ones presented in the literature, and evaluate the performance of these new elements by comparing with existing conventional and ANCF elements. As shown in this paper, complete cubic polynomials of the triangular plate elements have thirty-nine coefficients, or equivalently thirteen Bezier control points. Therefore, in addition to the position and gradient coordinates, the use of internal middle surface nodes and/or curvature coordinates is required to formulate an ANCF triangular plate element based on a complete set of polynomial basis functions. On the other hand, a lower-order triangular plate element with position and gradient nodal coordinates can be systematically developed by applying linear constraint equations. These constraint equations allow identifying dependent variables and reducing the number of the element coordinates. Alternatively, an incomplete cubic polynomial representation can be employed from the outset for obtaining a lowerorder triangular plate element that employs only position and gradient coordinates as nodal coordinates. In this study, these two systematic approaches are used for developing two new lower-order ANCF triangular plate elements and a new high-order ANCF triangular plate element based on Bezier geometry. Specifically, three new ANCF fully parameterized triangular plate elements are developed in this investigation; an ANCF four-node mixed-coordinate triangular plate element (FNMC), a three-node triangular plate element (TN1) obtained from the FNMC element by imposing a set of linear algebraic equations, and another three-node triangular plate element (TN2) obtained by using incomplete polynomials from the outset. All the proposed ANCF triangular plate elements are based on a full parameterization; have global nodal positions and gradient coordinates; ensure continuity of displacement, gradient, and rotation fields at the nodes; allow for easily creating arbitrary shell geometry by simply changing the values of the nodal coordinates in the reference stress-free configuration; and have geometry that is related to B-spline and NURBS (non-uniform rational B-splines) by a linear mapping [74]. The performance of the elements developed in this investigation is also evaluated by comparing their results and convergence characteristics with those of the previously developed conventional linear and quadratic triangular plate elements and the ANCF fully parameterized rectangular plate element available in the general purpose MBS software SIGMA/SAMS.
The structure and organization of this paper are as follows. In Sect. 2, the volume and Cartesian parameterizations used in this investigation to obtain compact and closed-form shape function expressions and to systematically assemble the mesh elements are introduced. The coordinate transformation between the volume and Cartesian gradients is developed in Sect. 3. In Sect. 4, the position field of the new ANCF/FNMC triangular plate element is obtained from a complete cubic Bezier interpolation. In Sect. 5, the lower-order ANCF/TN1 element shape functions are derived by imposing linear constraints on the center node of the ANCF/FNMC element. In Sect. 6, the lower-order ANCF/TN2 element shape functions are directly obtained from an incomplete set of cubic basis functions. The isoparametric properties of the proposed ANCF triangular plate elements are discussed in Sect. 7. In Sect. 8, the procedure for developing the equations of motion of the proposed ANCF elements is described. In Sect. 9, numerical results obtained using different types of triangular plate elements are analyzed and compared in order to evaluate the performance and convergence of the proposed elements. In Sect. 10, summary and conclusions of this investigation are provided.

Element parameterization
The continuum mechanics position vector gradients, defined by differentiation of the position vector field with respect to parameters (fibers or coordinate lines), play a fundamental role in the developments presented in this paper. In order to understand the fundamental differences between the triangular plate and shell elements used in the conventional FE literature and the ANCF fully parameterized triangular plate elements developed in this study, the conventional FE coordinate parameterization is analyzed in details in this section. To this end, a set of coordinate lines associated with the material fibers of a continuum defined by the Cartesian that the position field of a three-dimensional continuum r (x, t) = r 1 r 2 r 3 T can be written as follows: The use of the Cartesian parameterization is necessary to assemble triangular plate elements having gradient vectors as nodal coordinates. On the other hand, it is more convenient to use other dimensionless parameters for describing the element geometry. For this purpose, consider an arbitrary triangular plate element whose corner nodes are ordered counterclockwise and labeled 1, 2, and 3, as shown in Fig. 1. The positions of the element vertices, located on the element mid-surface, are defined in the global coordinate system by the vectors v k = x k y k z k T , k = 1, 2, 3, with x k , y k , and z k are the Cartesian coordinates at vertex k of the triangular plate element. The Cartesian coordinates x and the element corner node position vectors v k , k = 1, 2, 3 are defined in the same global coordinate system. The vector x can be written in terms of four dimensionless volume coordinates ξ= ξ η ζ χ T as As shown in Fig. 1, while the first three dimensionless volume coordinates ξ, η, and ζ define the coordinates of an arbitrary point on the plate mid-surface, the fourth dimensionless volume coordinate χ defines the location of an arbitrary plate surface from the mid-surface. It should be noted that the volume coordinates ξ, η, and ζ , associated with the continuum coordinate lines, are not independent parameters. The dimensionless coor- As shown in Fig. 2, each of the dimensionless volume coordinates ξ, η, and ζ ranges from zero to one, while the dimensionless volume coordinate χ ranges from − 1 to + 1. The coordinates ξ, η, and ζ form straight lines normal to the sides of the triangular plate element, whereas the coordinate χ forms straight lines orthogonal to the mid-surface of the element. Furthermore, for a general triangular plate element, the contour lines (isolines) are normal to the coordinate lines ξ, η, and ζ , while a set of contour planes (isoplanes) is associated with the dimensionless coordinate χ . The equation ξ = c, for example, describes a set of straight lines parallel to the triangular plate side identified by the vertices 2 and 3. Similar comments applies to the dimension-less volume coordinates η and ζ , as shown in Fig. 3.
On the other hand, χ = c can be geometrically interpreted as a set of straight planes called contour planes (isoplanes) parallel to the triangular plate mid-surface identified by the vertices 1, 2, and 3, as shown in Fig. 3. Let P be an arbitrary point inside the triangular plate whose vertices are relabeled for convenience as A, B, and C, which correspond to the vertex numbers 1, 2, and 3, respectively, as shown in Fig. 4, where is the volume of the triangular plate element, and 1 , 2 , and 3 represent the volumes of the three triangular plates P BC, PC A, and P AB. The first three volume coordinates ξ, η, and ζ can be written, respectively, as ξ = 1 / , η = 2 / , and ζ = 3 / . Furthermore, considering a triangular plate element of thickness W , one can define a transversal director vector w orthogonal to the triangular plate mid-surface in the straight (un-deformed) configuration as follows: where v 2,1 = v 2 − v 1 , v 3,1 = v 3 − v 1 , and v k , k = 1, 2, 3 are the global position vectors of the triangular plate corner nodes. It is important to note that transversal directors used in the FE literature do not always represent position vector gradients, and the cross product of two position vector gradients does not in general lead to another position vector gradient. The Cartesian coordinates x = x y z T and the volume coordinates ξ = ξ η ζ χ T are related by or equivalently: where x k , y k , and z k represent the Cartesian coordinates of the triangular plate vertex k and w i , i = 1, 2, 3 are the Cartesian components of the transversal director vector w. One can show, using Eq. 5, that the mid-surface triangular plate element corner nodes (x 1 , y 1 , z 1 ) , (x 2 , y 2 , z 2 ), and (x 3 , y 3 , z 3 ) correspond to the volume coordinates (1, 0, 0, 0) , (0, 1, 0, 0), and (0, 0, 1, 0), respectively. Using Eqs. 3 and 5, one can write ⎡ The inverse representation is ξ = V −1x , where V denotes the volume of the triangular plate element, which is equal to the determinant of the vertex matrix V, and the constant coefficients It follows that This equation can be written compactly as In this equation N is In the straight (un-deformed) configuration, it can be shown that the vectors normal to the triangular plate sides opposite to the triangular plate vertices 1, 2, and 3, respectively, while the vector n 4 = L 4,1 L 4,2 L 4,3 T is a vector normal to the triangular plate mid-surface. As mentioned before, the dependent volume coordinate lines ξ, η, and ζ are normal to the element sides, while the independent volume parameter χ is normal to the plate mid-surface, 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 /V = n T 1 dx/V , which represents an infinitesimal change of the dimensionless volume parameter ξ . Since an infinitesimal change of ξ can be expressed as dξ = n T 1 dx/V , the projection of the Cartesian vector dx = dx dy dz T along the vector n 1 = L 1,1 L 1,2 L 1,3 T is given by dξ and is orthogonal to the element side defined by the vertices 2 and 3. Similar comments apply to the dimensionless volume coordinates η and ζ . On the other hand, using Eq. 9 one can write dχ = L 4,1 dx + L 4,2 dy + L 4,3 dz /V = n T 4 dx/V , which represents an infinitesimal change of the volume parameter χ . Since one can write χ in terms of the vectors n 4 and dx as dχ = n T 4 dx/V , the projection of the Cartesian vector dx = dx dy dz T along the vector n 4 = L 4,1 L 4,2 L 4,3 T is given by dχ and is orthogonal to the element middle surface. Therefore, using Eqs. 9 and 12, one can write dξ = Ndx/V , which demonstrates that the infinitesimal change dξ can be defined using the product of N and dx. Because the vectors n 1 , n 2 , and n 3 that form the rows of N are normal to the element faces, n 4 is normal to the element mid-surface, the volume coordinate lines ξ, η, and ζ are normal to the element sides, and χ is along a direction normal to the element mid-surface; the gradients r ξ , r η , r ζ , and r χ , tangent to the coordinate lines ξ, η, ζ , and χ , are normal to the element sides associated with ξ, η, and ζ , and to the element mid-surface associated with χ , respectively. As will be discussed in Sect. 3, special definitions of the corner node coordinate lines will be used to facilitate the element assembly.

Position vector gradient transformations
In order to be able to use a standard FE assembly procedure with the proposed ANCF triangular plate mesh, a new set of volume parameters must be used to define volume gradient vectors parallel to the triangular plate sides. To be able to develop these position vector gradients, a different geometric definition of the volume coordinate lines is obtained by using the algebraic constraint equation ξ + η + ζ = 1 to eliminate a dependent volume parameter at each of the three element nodes.

Node parameterization
Consider the first triangular element node, denoted as node 1 and shown in Fig. 1. The first dimensionless volume coordinate ξ can be written as a linear function of the other two independent volume coordinates η i and ζ i as ξ = 1 − η i − ζ i , where the subscript i refers to independent parameters. The volume parameter χ , on the other hand, is always considered an independent parameter denoted as χ i . Eliminating ξ from Eq. 5, one obtains where Equivalently, the preceding equation can be written as where It is important to observe the change of the geometric meaning of the volume parameters as the result of substituting Eq. 3 into Eq. 5. For instance, as shown in Fig. 5, when ξ is eliminated because of the constraint equations between the element volume coordinates, the geometric interpretation of the independent volume parameters η i and ζ i changes. This important geometric concept is also described by the equation In this case, η i coordinate lines are a set of straight lines parallel to the element side defined by the vertices 2 and 1 (2-1 side), whereas the ζ i coordinate lines are straight lines parallel to the 3 -1 side. On the other hand, the χ i coordinate lines are always a set of straight lines orthogonal to the mid-surface of the plate. Therefore, an infinitesimal change in the Cartesian parameter vector x is the result of infinitesimal changes in the independent η i , ζ i , and χ i volume parameters along, respectively, the vectors v 2,1 , v 3,1 , and w, which can also be made unit vectors. At the first node, a constant transformation between the volume gradient vectors r η i , r ζ i , and r χ i , and the Cartesian gradient vectors r x , r y , and r z can be written as By including the position coordinates, the preceding equation can be written as The nodal coordinates at the first node in terms of the volume and Cartesian gradients given byē 1 = , respectively, are related byē 1 = T 1p1 . In this equation,T 1 is a 12 × 12 square transformation matrix defined as and I is the 3 × 3 identity matrix. Using the Cartesian nodal coordinates allows performing a standard FE assembly. On the other hand, using the independent volume parameterization, the vector of Cartesian coordinates can be written as Therefore, an infinitesimal change in the Cartesian coordinates x is due to infinitesimal changes in the independent volume coordinatesη i ,ζ i , andχ i . These independent infinitesimal changes are, respectively, along the three unit vectorŝ v 2,1 ,v 3,1 , andŵ which are not necessarily orthogonal unit vectors. The Cartesian gradients associated with the element straight (un-deformed) configura-tion are r x = 1 0 0 T , r y = 0 1 0 T , and r z = 0 0 1 T . Therefore, one can write rη i rζ i rχ i = ∂x/∂η i ∂x/∂ζ i ∂x/∂χ i = v 2,1v3,1ŵ since it can be shown that rη i =v 2,1 , rζ i =v 3,1 , and rχ i =ŵ. This equation demonstrates that, in the straight (undeformed) configuration, the volume gradients are tangent to the element sides that intersect at node 1 and to the normal of the plate mid-surface, respectively, as shown in Fig. 5. As previously discussed, the dimensionless coordinates η i and ζ i are independent volume parameters which can be geometrically interpreted as coordinate lines associated with straight lines parallel to the element sides. Therefore, the geometric meaning of the independent volume parameters used in this investigation is different from the geometric interpretation of the redundant parameters used in the conventional FE literature. Applying the same procedure for nodes 2 and 3, by writing η = 1 − ζ i − ξ i and ζ = 1 − ξ i − η i , respectively, one can define, as shown in Fig. 5, the transformation matrices at nodes 2 and 3 The transformations obtained in this section, which will facilitate the derivation of closed-form shape functions for the ANCF triangular plate element, also make the geometric interpretation of the gradient vectors used more obvious.

Element transformation matrix
Using the three transformation matrices developed in this section, and denoting the total vector of volume and Cartesian nodal coordinates asē = ē 1 T ē 2 T ē 3 T T , andp = p 1 T p 2 T p 3 T T , respectively, one can write a linear relationship between the vectors of nodal coordinates asē =Tp. In this equation, the transformation matrixT has dimensions 36 × 36 and is defined in a block diagonal form asT In the following section, a kinematic description based on Bezier geometry and a set of cubic basis functions are used in the development of the shape functions of the four-node mixed-coordinate (FNMC) ANCF triangular plate element. Using the FNMC shape functions, the shape functions of the three-node ANCF element (TN1) can be derived by using a set of linear constraint equations. Another three-node ANCF triangular plate element (TN2) is developed by using incomplete polynomials from the outset. It can be shown that the proposed FNMC and TN triangular elements are isoparametric, lead to exact representation of the rigidbody motion, define a constant mass matrix, and define complex shell geometry by simply changing the values of the nodal coordinates in the reference configuration thereby eliminating the distinction between plate and shell geometries as it is the case when conventional FE formulations are used.

Complete cubic polynomial four-node element
In this section, the shape functions of the ANCF/FNMC triangular plate element shown in Fig. 6 are developed using a special type of a cubic Bezier triangle on which the global position of an arbitrary point can be written as r = m k=1 g k q k , where m = 13, g k are the Bezier triangle basis functions defined in terms of the volume coordinates ξ, η, ζ , and χ , and q k are the Bezier control points. The Bezier polynomial terms g k , k = 1, 2, . . . , 10, are the terms which appear in the expansion of (ξ + η + ζ ) 3 , while the Bezier polynomial terms g k , k = 11, 12, 13, are the product of ξ, η, and ζ with χ . The cubic Bezier basis functions used in this investigation are The cubic Bezier triangle description can be used with a constant coordinate transformation to define the FNMC element assumed displacement field. The relationship between the cubic Bezier control points and the FNMC element nodal coordinates can be expressed as Using these linear algebraic equations, the FNMC position field is defined as r = Se, where S and e are, respectively, the FNMC shape function matrix and vector of nodal coordinates defined as and e = r 1 T r 1 where I is the 3 × 3 identity matrix, and r ξ i = ∂r/∂ξ i , r η i = ∂r/∂η i , r ζ i = ∂r/∂ζ i , and r χ i = ∂r/∂χ i are the nodal gradients defined by differentiation with respect to the independent volume parameters along the element sides and along the normal to the element mid-surface. The FNMC shape functions s k , k = 1, 2, . . . , 13, are The superscripts 1, 2, and 3 refer to corner nodes in the definition of the nodal coordinate vector, while superscript 4 refers to the FNMC center node as shown in Fig. 6. The element corner nodes are ordered counterclockwise in a manner that follows the right-hand rule.
The FNMC nodal coordinates are defined as In this equation, the volume gradient partial derivatives are performed using the volume coordinate constraint equation in the Bezier displacement field. The proposed FNMC element is an isoparametric finite element, describes correctly rigid-body motion, and can be used to develop shell geometry by using an appropriate value for the vector e in the reference configuration.

Elimination of the internal node
As discussed in the preceding section, the FNMC complete cubic polynomials have 39 coefficients, or equivalently 13 control points. For this reason, an internal node was introduced to allow for the use of a complete polynomial representation. Nonetheless, a lower-order three-node element can be systematically obtained by writing the position coordinates of the internal node in terms of the coordinates of a set of material point located on the sides of the element. This procedure leads to linear constraint equations that can be written as r 10 = 1 9 9 i=1 r i , where r i , i = 1, 2, . . . , 9 are global positions of material points distributed on the sides of the FNMC element and are defined as r (1, 0, 0, 0) , 0, 1, 0, 0) , Using the resulting algebraic constraint equations, a linear transformation can be developed and used to define the TN1 shape functions and position field which can be written asr =Sē, whereS andē are, respectively, the shape function matrix and vector of nodal coordinates defined as S = s 1 Is 2 Is 3 Is 4 Is 5 Is 6 Is 7 Is 8 Is 9 Is 10 Is 11 Is 12 I and where I is the 3 × 3 identity matrix, and r ξ i = ∂r/∂ξ i , r η i = ∂r/∂η i , r ζ i = ∂r/∂ζ i , and r χ i = ∂r/∂χ i are the element nodal gradients associated with independent volume parameters along the element sides and along the normal to the element midsurface, as previously discussed. The shape functions s k , k = 1, 2, . . . , 12, of the TN1 element are defined in as The TN1 triangular plate element has 36 nodal coordinates which consist of 9 position coordinates and 27 gradient coordinates. The vector of element coordinatesē which has the volume gradients is related to the vector of structural coordinatesp which has the Cartesian gradients by the constant transformationē =Tp, where the transformationT is defined in Sect. 3. The procedure described in this section, in which algebraic equations are used to eliminate the internal node, allows for systematically obtaining a symmetric incomplete polynomial representation without the need for using a trial and error process to identify the appropriate basis functions. As demonstrated in the following section, an alternate set of incomplete cubic Bezier basis functions can be used to obtain directly the shape functions using a linear mapping.

Use of incomplete cubic polynomial representation
In order to develop the shape functions of the TN2 element shown in Fig. 6, a special type of Bezier basis functions and a set of control points are used to define the global position vector of an arbitrary point as r = 12 k=1 h k q k , where h k are the Bezier basis or blending functions defined in terms of the volume coordinates ξ, η, ζ , and χ , and q k are the control points. The incomplete cubic Bezier basis functions used in this investigation are In this case, the relationship between the incomplete cubic Bezier control points and the TN2 element nodal coordinates can be written as Using these linear algebraic equations, the TN2 position field can be expressed as r =Sē, whereS andē are the TN2 shape function matrix and vector of nodal coordinates given, respectively, bȳS = s 1 Is 2 Is 3 Is 4 Is 5 Is 6 Is 7 Is 8 Is 9 Is 10 Is 11 Is 12 I where the nodal gradients have the same interpretation as in the case of the TN1 element. The TN2 shape functionss k , k = 1, 2, . . . , 12, arē The linear transformations between the incomplete cubic Bezier basis functions h k , k = 1, 2, . . . , 12, and the shape functions of the TN2 triangular plate element s k , k = 1, 2, . . . , 12, are given bȳs

Fig. 7 Continuum configurations
The TN2 triangular plate element has 36 nodal coordinates; 9 position and 27 gradient coordinates. As in the case of the TN1 element, the coordinate vectorē is related to the vector of structural coordinatesp by the linear equationē =Tp, whereT =T.

Isoparametric property and shell geometry
In this section, two important features of the higherorder ANCF elements proposed in this investigation are discussed. The first is the isoparametric property and the second is the shell geometry.

Isoparametric property
Isoparametric elements, which are capable of describing the element geometry and displacement using the same shape function matrix, correctly describe rigidbody motion and can be used in a non-incremental solution procedure. As shown in Fig. 7, three general configuration of a continuum body can be identified, namely a straight (un-deformed) configuration, a reference (curved) stress-free configuration, and a current (deformed) configuration. In the current configuration, one can write the element position field as r = Se = STp, where e and p are the coordinate vectors associated with the volume and Cartesian parame-terizations, respectively. In the reference configuration, the definition of the element position field is given by r 0 = Se 0 = STp 0 , where e 0 and p 0 are appropriate coordinate vectors associated with the volume and Cartesian parameterizations, respectively, which define the stress-free reference configuration and allow for describing a general curved geometry. In the straight (un-deformed) configuration, the element position field is x = Se s = STp s , where e s and p s are nodal coordinate vectors associated with the volume and Cartesian parameterization, respectively, that define the geometry of the element in the straight (un-deformed) configuration. One can then write J = (∂r/∂r 0 ) = (∂r/∂x) (∂x/∂r 0 ) = (∂r/∂x) (∂r 0 /∂x) −1 = J e J −1 0 , where J is the matrix of position vector gradients, and J e = ∂r/∂x and J 0 = ∂r 0 /∂x, as shown in Fig. 7, are matrices of position vector gradients defined using the straight configuration coordinate lines. A simple example can be used to show that all the ANCF triangular plate elements developed in this investigation are isoparametric. In order to demonstrate the isoparameteric property, consider the straight (un-deformed) configuration of the ANCF/FNMC element shown in The coordinates x k , y k , and z k are the Cartesian coordinates of the vertex k of the element. Using the relation ξ + η + ζ = 1, the ANCF/FNMC shape functions, and the expression of the transformation matrix T, one can show that In this equation, w 1 , w 2 , and w 3 are the Cartesian coordinates of the triangular plate director vector associated with the element straight (un-deformed) configuration. Equation 40 shows that the element straight (un-deformed) configuration can be represented by the ANCF/FNMC element shape functions using an appropriate set of nodal coordinates. Similar procedures can be used to demonstrate the isoparametric property of the ANCF/TN1 and ANCF/TN2 elements.

Shell geometry
In the conventional FE literature, there is a clear distinction between the assumed displacement fields of plates and shells. The main reason for making this distinction is attributed to the fact that most conventional plate and shell elements employ infinitesimal rotations as nodal coordinates. Because of the use of infinitesimal rotations, these plate and shell elements are not isoparametric elements, and as a consequence, their displacement fields cannot be used to describe arbitrary initial curved geometry. Furthermore, because of the use of infinitesimal rotations, the element displacement fields cannot be related to the more accurate computational geometry B-spline and NURBS representations using a linear mapping. Special techniques are often used in the conventional FE literature to develop shell elements [1,58,62]. Because the governing partial differential equations of plate and shell bending involve fourth-order derivative with respect to the spatial coordinates, the use of cubic polynomials is important in order to be able to correctly capture the effect of curvature changes within the element. The use of lowerorder conventional isoparametric elements which do not ensure the continuity of the rotation field at the nodal points may require the use of a very large number of elements to capture accurately the bending deformations. Nonetheless, even with the use of a very fine mesh, the rotations, stresses, and curvatures remain discontinuous at the nodal points. Another serious limitation of the conventional plate and shell elements is the restriction imposed on the definition of the elastic forces. For the conventional elements, the elastic forces cannot be easily formulated using a general continuum mechanics approach or general material models, and therefore, specialized theories with simplifying assumptions are often used [19,20,75].
For the ANCF elements proposed in this investigation, there is no distinction between the assumed displacement fields of plates and shells. By simply changing the values of the nodal coordinates in the reference configuration, initially curved and stress-free shell structures with complex shapes can be obtained. The shell reference configuration can simply be defined by the ANCF displacement field as X s = r 0 = S (x) e 0 , where e 0 is the vector of nodal coordinates in the reference configuration. By using the definition J = J e J −1 0 , as discussed in the preceding section, a stress-free initially curved configurations can be systematically achieved. Table 1 shows examples of some of the shell geometry that can be obtained with one ANCF triangular plate element. Other complex shell geometries can be created by simply changing the orientation and magnitude of the position vector gradients. Furthermore, the proposed ANCF triangular plate element assumed displacement fields are related to the B-spline and NURBS representations using a linear mapping, and therefore, they can be effectively used in the successful integration of computer-aided design and analysis (I-CAD-A). Furthermore, the generalized ANCF triangular plate element stress forces can be formulated using a general continuum mechanics approach and general material models.

Element equations of motion
In this section, the D'Alembert-Lagrange principle is used to develop the FE equations of motion which are based on a total Lagrangian description and allow     [76]. While the ANCF/FNMC element is used as an example, the same procedure can be used with the ANCF/TN1 and ANCF/TN2 elements.

Element kinematics
Using the FNMC position field, r = Se = STp, one can write δe = Tδp = STδp. The virtual change in the FNMC Cartesian gradients 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: where S x , S y , and S z are, respectively, the spatial derivatives of the matrix of shape functions with respect to the Cartesian coordinates x, y, and z. The evaluation of the Cartesian derivatives S x , S y , and S z is necessary for the formulation of the element elastic forces.

Element mass matrix
The virtual work of the FNMC element inertia forces is defined as δW i = − V ρr T δr |J 0 | dV , where ρ is the element mass density defined in the stress-free reference configuration,r is the absolute acceleration vector of an arbitrary element material point, |J 0 | is the determinant of the gradient matrix defined in the stress-free reference configuration, and V is the volume in the straight (un-deformed) configuration. By using the FNMC kinematic equations, one can write: where the element mass matrix M = T T V ρS T S |J 0 | dV ) T is a constant, symmetric, and positive definite matrix, and consequently, the centrifugal and Coriolis inertia forces are identically zero. The definition of the element mass matrix associated with the Cartesian coordinate vector allows for performing a standard FE assembly to obtain the mesh mass matrix.

Element external forces
One can write the virtual work of the external forces as δW e = V f T e δr |J 0 | dV , where f e is a vector of body forces such as the gravity and magnetic forces. Using the FNMC kinematic equations, the virtual work of the external forces can be rewritten as where Q e = T T V S T f e |J 0 | dV defines the vector of generalized external forces associated with the element generalized nodal coordinates. The use of the Cartesian gradients in the coordinate vector p allows for a standard FE generalized force assembly.

Element elastic forces
Using a general continuum mechanics approach, one can write the virtual work of the element elastic forces as δW s = − V σ : δε |J 0 | dV , where σ denotes the symmetric second Piola-Kirchhoff stress tensor, ε represents the symmetric Green-Lagrange strain tensor, and the symbol : defines tensor double contraction. All the elements proposed in this investigation allow using general constitutive models in the elastic force formulation. Using the FNMC kinematics, one can write the virtual work of the elastic forces as where Q s = −T T V σ : ∂ε ∂e T |J 0 | dV is the nonlinear vector of the element generalized elastic forces which can be evaluated using Gauss quadrature rules. A standard FE assembly procedure can be adopted to define the body elastic forces because of the use of Cartesian gradients in p.

System equations of motion
The system equations of motion of the ANCF/FNMC element mesh can be developed using the principle of virtual work which can be stated 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. Applying this principle, the element equations of motion can be written as Mp = Q e + Q s . As mentioned previously, the use of the Cartesian gradients r x , r y , and r z in the generalized coordinate vector p allows for standard FE assembly of the element vectors and matrices.

Numerical results and discussion
In this section, several FE models are developed in order to demonstrate the use and evaluate the performance of the proposed FNMC, TN1, and TN2 elements. The numerical results obtained using these elements are compared with the results of the conventional three-node linear (TNL) and six-node quadratic (SNQ) triangular plate elements [77]. The shape functions of the conventional TNL and SNQ elements are presented in the Appendix of the paper. While these conventional elements ensure only the continuity of the position field at the nodes, the proposed ANCF elements guarantee the continuity of both the position and gradient fields, and therefore, the continuity of the rotation, strain, and stress fields is also ensured. The conventional TNL and SNQ elements are implemented using a total Lagrangian approach and the results are obtained using a non-incremental solution procedure [78,79]. Additionally, the numerical results obtained using the proposed ANCF triangular plate elements are compared with the numerical results of the ANCF fully parameterized rectangular plate element developed by Mikkola and Shabana [34] and implemented in the general purpose MBS computer program SIGMA/SAMS (Systematic Integration of Geometric Modeling and Analysis for the Simulation of Articulated Mechanical Systems). It is important, however, to point out that because of the low-order of polynomial interpolation of the conventional elements, the curvature within these elements is either constant or identically zero, and therefore, these conventional elements are not suitable for describing large deformations and extreme bending scenarios.
In this investigation, a MATLAB computer program was developed in order to obtain the numerical results presented in this section for the conventional and proposed ANCF elements. In order to correctly capture geometric nonlinearities, a general continuum mechanics approach is employed for evaluating the generalized elastic forces. The Saint Venant-Kirchhoff hyperelastic constitutive law and a standard Gauss quadrature procedure based on a full integration scheme are used for the evaluation of the generalized elastic forces [58]. For the ANCF elements developed in this investigation, full integration requires combination of 12 triangular Gauss integration points distributed along the plate mid-surface; this corresponds to the standard sixth-order Gauss integration quadrature scheme for triangular elements, with 2 linear Gauss integration points along the thickness, which corresponds to the standard third-order Gauss integration quadrature scheme for linear elements. The choice of the type and number of integration points was made based on numerical experimentation. It was found that using a larger number of quadrature points does not improve the accuracy of the numerical results, whereas using a smaller number of integration points can lead to less accurate results. The locations of the Gauss integration points used for performing the numerical integration of the proposed ANCF fully parameterized triangular plate elements are shown in Fig. 8, while the corresponding quadrature nodes and weights are presented in Table 2. The numerical integration of the equations of motion is performed using the explicit multistep fourth-  order Adams-Bashforth numerical integration method with a constant time step [80].

Rigid-body motion
The stiff rectangular pendulum example shown in Fig. 9 is used to demonstrate that both the conventional and the proposed ANCF elements can correctly describe rigid-body motion. The rectangular pendulum is subjected to a constant gravitational force field and is assumed to be initially in a horizontal configuration. The pendulum has length L = 0.1 m, width H = 0.1 m, thickness W = 0.01 m, mass density ρ = 7860 kg/m 3 , Poisson ratio ν = 0.27, and Young's modulus E = 2.1 × 10 9 N/m 2 . The gravity acceleration considered is g = 9.81 m/s 2 . Based on the assumed inertia and geometry properties, the pendulum  compare the results of this numerical example with the results obtained using a rigid-body pendulum model, a large elastic modulus is used for the flexible pendulum models. Figure 10 shows the FE mesh used for modeling the stiff pendulum system. Only one rectangular macro-element, shown in Fig. 11, which consists of two triangular plate elements, is used. As shown in Fig. 10, two revolute joints located at the points A 1 and A 2 are used to constrain the flexible pendulum motion with respect to the ground. On the other hand, a revolute joint having a joint axis passing through the two points A 1 and A 2 is used to constrain the motion of the rigid pendulum. The constraint equations for the flexible pendulum are linear algebraic equations, and therefore, are systematically eliminated at a preprocessing stage. In addition to the proposed FNMC, TN1, and TN2 element models, two other models based on the conventional TNL and SNQ elements are developed. The duration of the dynamic simulation is assumed T = 1.0 s, while the minimum time step used in the numerical integration algorithm is t = 5.0 × 10 −6 s. The vertical displacement of the plate pendulum tip point B obtained using the stiff flexible and rigid-body models is shown in Fig. 12. The results presented in this figure show a good agreement, demonstrating that the ANCF elements developed in this investigation can correctly describe rigid-body motion.

Patch test
In order to perform the patch test [58,62], the rectangular plate shown in Fig. 13 with a mesh of 10 triangular plate elements, shown in Fig. 14, is used. The model geometric and material properties are the same as previously presented in this section and the nodal coordinates considered are shown in Table 3. In order to eliminate all the rigid-body modes, a minimum set of displacement and gradient boundary conditions are used for constraining the rectangular plate to the ground at the points A 1 and A 2 . A modal analysis was performed in order to obtain the natural frequencies, and no zero spurious eigenvalues were found. Therefore, the states of constant strain can be correctly captured by all the triangular plate elements considered in this investiga-  tion. This simple example demonstrates that the ANCF elements developed in this paper pass the patch test, are compatible elements, and converge when the number of elements is properly increased [61].

Convergence analysis
In order to check the element convergence in the case of finite rotations and large displacements, the initially straight flexible rectangular pendulum shown in Fig. 9 is considered. The pendulum length, width, thickness, and modulus of elasticity considered in this numerical example are L = 1.    show, respectively, the FE mesh and the rectangular macro-element used. The TNL, SNQ, TN1, TN2, and FNMC element models are developed using a symmetric mesh made of rectangular macro-elements. The  total number of triangular plate elements that form each mesh is N e = 8N 2 m , where N m is the number of macro-elements, each of which is constructed from 8 triangular plate elements. Because the convergence rate is not the same for the five types of triangular plate elements considered, different numbers of elements are used to obtain convergent solutions for different plate models. It was found that the displacement and strain convergence requires a 16 macro-element (N e = 8 × 16 2 = 2048) for the TNL element; a 14 macro-element (N e = 8 × 14 2 = 1568) for the SNQ element; a 14 macro-element mesh (N e = 8 × 14 2 = 1568) for the TN1 element; a 14 macro-element mesh (N e = 8 × 14 2 = 1568) for the TN2 element; and a 12 macro-element mesh (N e = 8 × 12 2 = 1152) for the FNMC element. However, in the case of the conventional TNL and SNQ models, finer meshes cannot achieve geometric convergence because of the gradient discontinuities at the element interface nodes that lead to rotation, strain, and stress discontinuities. In order to perform a comparative convergence rate study of the elements considered in this investigation, the root-mean-square errors of the numerical results are evaluated considering a set of reference solutions and the resulting normalized values are reported in Tables  4, 5, 6, 7, and 8. The normalized root-mean-square errors of the numerical results are computed for the complete simulation interval considering the displacements and the normal strains ε x x , ε yy , and ε zz . The reference solution corresponds to the 16 macro-element mesh for the TNL element; 14 macro-element mesh for the SNQ element; 14 macro-element mesh for the TN1 element; 14 macro-element mesh for the TN2 element; and 12 macro-element mesh for the FNMC element. As expected, for all elements considered, the normalized root-mean-square errors decrease as the number of elements increases. The TNL element has the slowest convergence rate, while the TN1 and TN2 elements have a convergence rate comparable with the convergence rate of the SNQ element. The FNMC element, on the other hand, has the fastest convergence rate. The difference in the convergence rates can be attributed to the fact that the TN1 and TN2 elements are based on an incomplete set of cubic polynomials, whereas the shape functions of the FNMC element are obtained from a complete set of cubic polynomial basis functions.   The numerical results obtained using the proposed ANCF elements are also compared with the numerical results obtained using the ANCF rectangular plate element developed by Mikkola and Shabana [34]. To this end, a 10×10 mesh, a 12×12 mesh, and a 14×14 mesh made of ANCF rectangular plate elements are developed and it was found that the displacement and strain convergence is achieved for this ANCF element using the 14 × 14 element model. Figures 17, 18, 19, and 20 show, respectively, the vertical displacement of point B on the pendulum under the effect of gravity, and the normal strains ε x x , ε yy , and ε zz at the pendulum center point C when convergence is achieved for the TNL and SNQ element models. Because nodal strain continuity is not ensured by the conventional triangular plate elements, averaging techniques were adopted to improve the accuracy of the strain results of these conventional elements. Figures 21, 22, 23, and 24 show, respectively, the point B vertical displacement under the effect of gravity, and the normal strains ε x x , ε yy , and ε zz at the pendulum center point C when convergence is achieved for the ANCF rectangular plate, TN1, TN2, and FNMC element models. Since the ANCF triangular plate elements guarantee the continuity of the rotation, strain, and stress fields at the nodes, there is no need for using strain averaging techniques. The numerical results presented in this section demonstrate that there is, in general, a good agreement between the numerical solutions obtained using different triangular plate element models. However, the strain results obtained using the TNL and SNQ element models show slightly different trends because such elements do not guarantee the continuity of the rotation, strain, and stress fields at the nodal points.

Continuity conditions
In order to demonstrate the gradient and strain continuity that characterize ANCF elements, the flexible straight rectangular pendulum shown in Fig. 25 is modeled using a mesh of one macro-element composed of eight FNMC triangular plate elements. The dynamic simulation is carried out using the geometric and material properties previously used. Figures 26, 27, and 28 show the normal strains ε x x , ε yy , and ε zz at the pendulum center point C, while Figs. 29, 30, and 31 show the normal stresses σ x x , σ yy , and σ zz at the same point C predicted using the coordinates of all the eight FNMC elements of the rectangular macro-element. The numerical results presented in Figs. 26,27,28,29,30, and 31 clearly demonstrate the continuity of the normal strains and stresses; similar continuity results are obtained for the shear strains and stresses. The proposed TN1 and TN2 elements also ensure the gradient, strain, and stress continuity at the nodes.

Quasi-periodic motion
In order to demonstrate that the proposed ANCF elements can capture the quasi-periodic motion that characterizes the dynamic simulation of the flexible pendulum system, the straight rectangular pendulum shown in Fig. 9 is modeled using a mesh of one macro-element composed of eight FNMC triangular plate elements. The dynamic simulation is carried out using the geometric and material properties used previously except for the elasticity modulus and the time interval which are assumed to be E = 2.1 × 10 7 N/m 2 and T = 10.0 s. Figure 32 shows the vertical displacement of

Summary and conclusions
In the conventional FE literature, distinction is often made between plate and shell elements as the results of using rotations as nodal coordinates. The use of this nomial representation, three-node triangular plate element (TN1) obtained from the FNMC element by eliminating the internal node, and three-node element (TN2) obtained by using an incomplete cubic polynomial representation from the outset. A cubic polynomial interpolation and a constant transformation between the ANCF nodal coordinates and the Bezier control points are used to develop the ANCF fully parameterized triangular plate elements developed in this study. All the proposed ANCF elements are isoparametric; describe correctly arbitrary rigid-body motion; lead to constant mass matrix and zero Coriolis and centrifugal forces; allow for consistently representing large finite rotations; ensure continuity of the gradients, rotations, strains, and stresses at the nodal points; and impose no restrictions on the value of the deformation within the finite element. Two coordinate parameterizations are used for describing the kinematics of the proposed elements, namely volume and Cartesian parameterizations. In order to be able to use a standard FE assembly procedure, the volume position vector gradients are first used to develop closed-form and compact element shape function expressions. The numerical results obtained using the proposed elements are compared with the numerical results obtained using the conventional elements (TNL and SNQ) as well as an ANCF fully parameterized rectangular plate element.
The results obtained in this study using the proposed elements showed, in general, good convergence characteristics. TNL and SNQ elements have the nodal position vector r k and the slope vector r k z as nodal coordinates. The linear shape functions s k , k = 1, 2, . . . , 6, of the conventional TNL element can be explicitly written in terms of the set of volume coordinates ξ, η, ζ , and χ as follows: where W is the TNL element thickness. The pairs of shape functions s i and s j , i = 1, 3, 5, j = 2, 4, 6, are, respectively, associated with the TNL nodal positions and slopes of the corner nodes k, k = 1, 2, 3.