A 3-Node Co-Rotational Triangular Finite Element for Non-Smooth, Folded and Multi-Shell Laminated Composite Structures

Based on the first-order shear deformation theory, a 3-node co-rotational triangular finite element formulation is developed for large deformation modeling of non-smooth, folded and multi-shell laminated composite structures. The two smaller components of the mid-surface normal vector of shell at a node are defined as nodal rotational variables in the co-rotational local coordinate system. In the global coordinate system, two smaller components of one vector, together with the smallest or second smallest component of another vector, of an orthogonal triad at a node on a non-smooth intersection of plates and/or shells are defined as rotational variables, whereas the two smaller components of the mid-surface normal vector at a node on the smooth part of the plate or shell (away from non-smooth intersections) are defined as rotational variables. All these vectorial rotational variables can be updated in an additive manner during an incremental solution procedure, and thus improve the computational efficiency in the nonlinear solution of these composite shell structures. Due to the commutativity of all nodal variables in calculating of the second derivatives of the local nodal variables with respect to global nodal variables, and the second derivatives of the strain energy functional with respect to local nodal variables, symmetric tangent stiffness matrices in local and global coordinate systems are obtained. To overcome shear locking, the assumed transverse shear strains obtained from the line-integration approach are employed. The reliability and computational accuracy of the present 3-node triangular shell finite element are verified through modeling two patch tests, several smooth and non-smooth laminated composite shells undergoing large displacements and large rotations.

excellent damage tolerance, superior fatigue response characteristics, and good damping behaviors under dynamic loads. By choosing an appropriate combination of reinforcement and matrix material, manufacturers can produce properties that exactly fit the requirements of a particular structure design [1]. The mechanical properties of laminated composite structures are sensitive to the lamination scheme and the ply orientation angle, so they often show unique responses even under simple loading conditions and geometric configurations. Furthermore, the anisotropic constitutive responses and the complexity of shell geometries, such as non-smooth and folded shell structures, make it challenging to perform accurate structural analysis, especially when large deformations are involved. Therefore, the development of reliable and efficient finite element methods for laminated composite shell structures are important [1][2][3][4][5].
Various computational formulations have been proposed for modeling composite shells and plates, and can be broadly classified into three categories: (1) single-director theories with anisotropic constitutive relations and (2) multi-director theories, which include multi-layer formulations, within which each layer has a single director with its own anisotropic constitutive relation, and (3) 3-D continuum theories. Examples of single-director theories include the classical laminated plate/shell theory (CLPT), the first-order shear deformation laminated plate/shell theory (FSDT), and the higher-order shear deformation laminated plate/shell theories (HSDT). Examples of multi-director theories include the layer-wise theory (LWT) and the zig-zag theory (ZZT). Examples of 3-D continuum theories include solid-shell formulations with a single layer or a multilayer structure.
In the single-director category, the CLPT is restricted to thin shell structures, as the effects of transverse shear strains and thickness strains are ignored. Based on the CLPT, Madenci et al. [6] proposed a free-formulation-based 3-node flat triangular shell element for geometrically nonlinear analysis of thin composite shells. Kapania et al. [7] presented a 3-node triangular flat shell element by combining a discrete Kirchhoff plate bending element with a membrane element for linear static, free vibration and thermal analysis of laminated plates and shells. Bisegna et al. [8] proposed a co-rotational triangular facet shell element for geometric nonlinear analysis of thin piezo-actuated structures.
Also in the single-director category, under the FSDT, the transverse shear strains are assumed to remain constant through the thickness, and the shell normal does not need to remain perpendicular to the mid-surface after deformation, while the inextensibility of transverse shell normal is assumed. Based on this theory, Peng et al. [9] proposed a meshfree method for bending analysis of folded laminated plates. Pham et al. [10] presented a combination of the edge-based smoothed finite element method (ES-FEM) and the three-node triangular elements (MITC3) for static responses and free vibration of laminated composite shells. Truong-Thi et al. [11] presented an extension of the cell-based smoothed discrete shear gap method using three-node triangular elements for the static and free vibration analyses of carbon nano-tube reinforced composite plates. Zhang et al. [12,13] developed an eight-node quadrilateral plate element with five mechanical degrees of freedom and one electric degree of freedom for static and dynamic analyses of piezoelectric integrated carbon nanotube reinforced functionally graded composite structures. Kreja et al. [14] presented isoparametric eight-node Serendipity-type shell finite elements to check the relevance of five-and six-parameter variants for large rotation plate and shell problems. The FSDT provides a balance between computational efficiency and accuracy for the global structural behaviors of thin and moderately thick laminated composite shells, while the local effects (e.g., inter-laminar stress distribution between layers, delamination, etc.) are often difficult to capture.
The HSDT provides a more accurate description of the transverse shear stress distributions by introducing more independent displacement parameters. Within the HSDT, Chen et al. [15] presented a refined three-node triangular element satisfying the requirement of C 1 weak-continuity. Tran et al. [16,17] proposed an edge-based smoothing discrete shear gap method using 3-node triangular elements combined with a C 0 -type higher-order shear deformation theory for static, free vibration and buckling analyses of laminated composite plates. Jin et al. [18] proposed a computationally efficient C 0 -type 3-node triangular plate element with linear interpolation functions for the analysis of multi-layered composite plates based on the mixed global-local higher-order theory.
In the multi-director category, the LWT assumes a layer-wise deformation pattern, and it can predict the interlaminar stresses accurately. Liu et al. [19] employed a layer-wise three-node triangular shell element for modeling the opening and shear modes of delamination. Phung-Van et al. [20] presented an extension of the cell-based smoothed discrete shear gap method using three-node triangular elements for dynamic responses of sandwich and laminated composite plates. Marjanović et al. [21] presented a triangular layered finite element with delamination degrees of freedom for composite shells based on the generalized LWT. For large deformation and large overall motion, Vu-Quoc et al. [22,23] developed the dynamic formulations for geometrically-exact multilayered composite beams, plates, and shells with ply drop-offs. However, layer-wise models are computationally expensive since the number of unknowns depends on the number of the layers of the laminates.
The ZZT describes a piecewise continuous displacement field in the plate thickness direction and fulfills interlaminar continuity of transverse stresses at each layer interface. Carrera [24] provides a comprehensive review of the ZZT. Versino et al. [25] developed six-and three-node triangular plate elements for homogeneous, multilayer composite and sandwich plates based on the refined zigzag theory. Nguyen et al. [26] developed a three-node triangular finite element for visco-elastic composite laminates based on a high-order zigzag theory. Wu et al. [27] proposed an efficient three-node triangular element with linear shape functions to model sandwich plates based on the refined higher-order zig-zag model in conjunction with the three-field Hu-Washizu variational principle. Liang et al. [28,29] proposed an efficient zigzag kinematic model for composite laminates with multiple alternating stiff-soft layers, which has been realized within a corotational shell element utilizing additional zigzag DOFs that are not subject to the corotational transformations, making use of a 2D local shell system over the surface of the structure.
The 3D continuum-based theory accounts for fully 3D constitutive behaviors, so the interlaminar stress of composite laminates can be effectively captured. Houmat [30] studied the free vibration of variable stiffness laminated composite plates using the 3D elasticity theory and the p-version finite element method. Ye et al. [31] used the scaled boundary finite element method to analyze the bending behaviors of the angle-ply composite laminated cylindrical shells based on 3D theory of elasticity. Kumari et al. [32] adopted the Reissner variational principle and the extended Kantorovich method to study the bending problem of composite cylindrical shells with arbitrarily support boundary conditions. Vu-Quoc et al. [33] developed multilayered composite solid-shell formulations, for both static and dynamic analyses, that could accommodate 3-D constitutive relations without a need to impose zero transverse normal stresses. Fan et al. [34] developed lowest-order (8-node hexahedral), and higher-order (32-node hexahedral) 3-D continuum solidshell elements, based on the theory of 3D solid mechanics, for static and dynamic analyses of composite laminates.
Despite these developments, numerical formulations based on the theories of HSDT, LWT, ZZT and 3D continuum-based theory often lead to high computational costs, which is a major concern in their practical applications. In the present study, a 3-node co-rotational triangular composite shell finite element is developed based on FSDT, where vectorial rotational variables are employed as rotational variables, two smaller components of one vector, together with the smallest or second smallest component of another vector, of an orthogonal triad initially oriented along the global coordinate system axes at each node on a non-smooth intersection of plates and/or shells are defined as vectorial rotational variables, while two smaller components of the mid-surface normal vector of shell at other nodes are defined as vectorial rotational variables. The resulting element tangent stiffness matrices are symmetric owing to the commutativity of nodal variables in calculating the second derivatives of strain energy with respect to local nodal variables and the second derivatives of the local nodal variables with respect to global nodal variables. Using such vectorial rotational variables, triangular and quadrilateral shell elements have been developed for large displacement and large rotation analyses of smooth shell structures [35][36][37][38], as well as non-smooth shells made of isotropic elastic materials [39,40]. The 3-node finite element formulation proposed in the present study is capable of modeling both smooth and non-smooth laminated composite shell structures experiencing large deformations.
There have been continuous efforts in developing triangular shell finite elements with high computational accuracy and convergence [41,42]. It is well-known that triangular finite elements may suffer from locking phenomena with thin shell/plate thickness, which results in the deterioration of computational accuracy and convergence. Up to now, there is no optimal 3-node triangular shell element, the main reason is that the strain distributions derived from displacement shape functions are often wildly different than expected. Barlow points, at which the strains derived from displacement shape functions are "correct" for all desired deformation modes and all permitted initial element shapes, can be found in a quadrilateral element, and they can be employed in calculating the element tangent stiffness matrix to eliminate locking phenomena. There exist no Barlow points, however, in a triangular shell element, thus there is no set of integration points that provides zero membrane strain for all cases of pure bending about axes parallel to each of the three edges of a triangular element. Most methods for eliminating locking problems of quadrilateral shell elements do not work well in triangular elements, such as reduced integration technique [43] and Hellinger-Reissner mixed formula method [44]. Macneal [45] proposed a line integration method to overcome locking problems, where assumed membrane strains and assumed shear strains are calculated respectively from the edge member membrane strains and the edge member transverse shear strains. Bletzinger [46] presented a discrete shear gap method (DSG) which utilizes only the usual displacement and rotational degrees of freedom at the nodes. Kim et al. [47] developed a 3-node macro triangular element using assumed natural strain method (ANS) for geometrically non-linear analysis of plates and shells. The ANS formula and the macro element strategy can reduce the locking effect and preserve its advantages in preprocessing. Argyris et al. [48,49] proposed a facet triangular shell element (TRIC) for nonlinear dynamic and elasto-plastic shell analyses using the natural mode method. Lee et al. [50] presented a Mixed Interpolation of Tensorial Components approach (MITC) for shell elements, which have been extended to MITC3 and MITC3+ shell finite elements by introducing interpolation cover functions and cubic bubble function for the rotations [51,52]. Cai et al. [53] proposed a lockingfree discrete shear triangular plate element without any numerical expediencies such as the reduce integration and assumed strains method. In the present study, transverse shear strains of the shell finite element are replaced with assumed shear strains obtained using the line integration approach [35,45] to overcome shear locking.
The outline of the paper is as follows. Section 2 describes the co-rotational framework and the kinematics of the 3-node triangular composite shell element. Section 3 presents the composite shell finite element formulation in the co-rotational local co-ordinate system. Section 4 gives the transformation relationship between the local and global responses. Several numerical examples are analyzed in Section 5 to verify the numerical accuracy of the present finite element. Conclusions are presented in Section 6.

Co-Rotational Framework
The co-rotational framework for the 3-node shell finite element is depicted in Fig. 1. The origin of the local co-rotational coordinate system of the element coincides with Node 1 in the initial configuration, and the coordinate system rotates with the element's rigid body rotation, but does not deform with the element. Although the adopted definition is not invariant to nodal ordering [54], the associated variability is negligible for large-displacement small strain problems. To define the local system, we firstly calculate the vectors v 120 and v 130 : where X i0 (i = 1, 2, 3) is the position of Node i in the initial global coordinate system. The base vectors of the local coordinate system in the initial configuration are calculated as follows: In the deformed configuration, the vectors v 12 and v 13 are defined as where d i (i = 1, 2, 3) is the translational displacement vector of Node i in the global coordinate system. Accordingly, the base vectors of the current local coordinate system can be obtained In the global coordinate system, we define the following vector that consists of all the global DOFs (degrees of freedom) for each element: where d T i = U i V i W i are three translational displacement DOFs of Node i, and n T gi denotes the vectorial rotational DOFs of Node i. For any node of smooth shells or located away from non-smooth shells intersections, n T gi = p i,n p i,m is defined to represent the vectorial rotational DOFs, in which p i,n and p i,m denote the two smaller global components of the shell director p i . The remaining component of the shell director p i is defined as where s 4 = ±1 takes the same sign as that of p i,l in the previous loading step, and {n, m, l} denotes the circular permutation of {X , Y , Z}.
On the other hand, if Node i is located on the intersections of non-smooth shells, the vectorial rotational DOFs of Node i is defined as n T gi = e iy,n e iy,m e iz,n . Here, e iy,n and e iy,m are the two smaller components of e iy , and e iz,n is the smallest or second smallest component of e iz , where e iy and e iz are two orientation vectors of an orthogonal triad at Node i. At each incremental loading step, the following relationships can be derived: The remaining components of the vectors e iy and e iz can be calculated by using e iy,n , e iy,m and e iz,n according to Eqs. (7a)-(7c). The sign symbols (s 1 , s 3 ) defined as 1 or −1 are respectively equal to the signs of components e iy,l and e iz,m of the last incremental step, and we have s 2 = −s 1 · s 3 . The vector e ix is calculated as the cross-product of e iy and e iz : Since the norm of a unit vector is identical to 1, defining the vectorial rotational variables as above can avoid ill-conditioning in updating the mid-surface normal vector at a node on the smooth part of the plate or shell (away from non-smooth intersections) or orientation vectors of an orthogonal triad at a node on a non-smooth intersection of plates or shells by properly controlling the size of loading step in a nonlinear incremental solution procedure.
There are 15 degrees of freedom per element in the local coordinate system where t T i = u i v i w i are three translational displacement components, θ T i = r i,x r i,y are vectorial rotational components, r i,x and r i,y are the components of the shell director at Node i along the x-axis and y-axis directions in the local coordinate system.
The relationship between the local and global translational displacements can be expressed as follows: where, At any node of smooth shells or any node away from non-smooth shell intersections, the relationship between the shell directors expressed in the local and global coordinate systems can be described as At any node on intersections of non-smooth shells, the following relationships between the shell directors in the local and global coordinate systems hold: For convenience, e ix0 ,e iy0 and e iz0 in the initial configuration are chosen to be coincident with the base vectors of the global coordinate system.

Description of the Element Geometry and Kinematics
The finite element shape functions are expressed in the natural coordinate system as follows: The displacement at any point of the element is calculated as follows: In the initial configuration, the shell director at Node i is calculated as the cross product of the tangent vectors along two natural coordinate axes. where To ensure the uniqueness of the shell director at any node shared by multiple adjacent elements in smooth shell regions, the following averaging procedure is adopted: The Green-Lagrange strain of the nonlinear shallow shell theory is adopted. For convenience, the strain vector is divided into membrane strain vector ε m , bending strain vector z l χ and transverse shear strain vector γ .
where, z l is the coordinate of the material point along the shell thickness direction, a is the thickness of the element, and

Laminated Composite Shell Finite Element Formulation in the Co-Rotational Local Coordinate System
The potential energy of a 3-node triangular composite shell finite element is defined as where W e is the work done by the external load, C is the elastic constitutive matrix. In the present study, we consider a laminated shell structure composed of long-fiber-reinforced materials. The transversely isotropic elastic constitutive matrix of each individual lamina is expressed in the material coordinate system as where, μ 12 and μ 21 are the Poisson's ratios, E 1 and E 2 are the elastic moduli, G 12 , G 23 and G 31 are the shear moduli [1,14].
By enforcing the stationarity condition to the potential energy, we have which yields the element internal force vector: where u L , f ext and f int are element nodal displacement vector, external load vector and internal force vector in the co-rotational local coordinate system, respectively. B m , B b and B γ contain the first-order derivatives of the membrane strain, bending strain, and shear strain, respectively. Details of these matrices are given in Appendix A.
By taking the first-order derivative of the internal force vector with respect to the local nodal variables, a symmetric element tangent stiffness matrix is obtained For elastic laminated shell elements, the integration along the shell thickness direction is decoupled from the integration on the mid-surface, so Eqs. (27) and (28) involving volume integrals can be rewritten in the following surface integral form: where D eq1 , D eq2 and D eq3 are the equivalent elastic matrices: n denotes the number of laminae; ζ i is the dimensionless coordinate of the interface between Layer i−1 and Layer i relative to the mid-surface: D ci is the elastic matrix of the ith layer lamina: where C i is the transversely isotropic elastic matrix of Layer i in the material coordinate system, which is obtained by replacing the material properties of Eq. (25) with those of the ith layer; T i is the transformation matrix between the material coordinate system of the ith layer lamina and the co-rotational element local coordinate system: where ϕ i is defined as the ply orientation angle measured between the reinforcement fiber of the ith layer lamina and the co-rotational local coordinate system of the composite element. It is noted that, for symmetrical laminates, the third integrals of Eqs. (29) and (30) are equal to zero.
To alleviate shear locking phenomenon, the assumed transverse shear strain vector and its first-order derivatives with respect to local nodal variables are employed. The modified line integration approach [35,45] is adopted to calculate the assumed shear strains, and by replacing the conforming transverse shear strain and its first-order derivatives with the assumed strains and its first-order derivatives with respect to local nodal variables, the following assumed strain finite element formulation is obtained: whereγ is the assumed transverse shear strain vector,B γ is its first-order derivative with respect to local nodal variables u L (refer to Eqs. (30)- (41) in reference [35]).

Stiffness Matrix of the 3-Node Shell Finite Element in the Global Coordinate System
The internal force f G of the 3-node triangular composite shell element in the global coordinate system can be obtained from the local element internal force f int as follows: where T is a transformation matrix consisting of the first-order derivatives of the local nodal variables with respect to the global nodal variables, which can be readily determined from Eqs. (10), (13a), (13b) and (14a), (14b): The element tangent stiffness K TG in the global coordinate system can be obtained by calculating the first-order derivatives of the global internal force vector with respect to the global nodal variables: where the sub-matrices in (40) involving the second-order derivatives are given in Appendix B. Considering the commutativity of the global nodal variables in the differentiation of (40), the second term in the right-hand side of (39) is symmetric, and since the first term in the right-hand side of (39) is also symmetric, the resulting element tangent stiffness matrix in the global coordinate system is symmetric.

Numerical Examples
To verify the reliability and computational accuracy of present 3-node co-rotational triangular composite shell element, two patch tests, and several smooth and non-smooth composite shell problems are solved, and the solutions are compared to numerical results from literatures [55][56][57][58][59]. In the following examples, the present 3-node triangular shell finite element using assumed transverse shear strains to replace conforming transverse shear strains is designated as "CR3 T" element. To show the element's spatial isotropy (i.e., the element stiffness matrices are independent of nodal ordering [41,50]

Patch Tests
Patch tests for the membrane behavior and the transverse out-of-plane bending behavior of plate and shell elements were suggested by MacNeal et al. [55], whereas rectangular plate with a length L = 0.24, width W = 0.12 and thickness h = 0.001, Young's modulus E = 10 6 , and Poisson's ratio μ = 0.25, is considered.
In the membrane patch test, the displacements u, υ, wat the boundary nodes of the rectangular plate (Fig. 3) are prescribed by u = 10 −3 (X + Y/2), υ = 10 −3 (Y + (X/2)) and w = 0. If this patch test is solved as a linear problem, the displacements at any point of the plate can also be calculated from above equations. The theoretical solution of this plate is a constant in-plane membrane stress field given by σ x = σ y = 1333, τ xy = 400 10 elements are employed in modeling of the rectangular plate (Fig. 3), and the CR3T element passes the membrane patch test exactly. To construct a constant stress state of the plate under out-of-plane bending, the displacements (u, υ, w) and rotations (θ X , θ Y ) at any point of the plate mid-surface are prescribed by For a linear problem, the theoretical solution for the stresses at the top and bottom surfaces of the plate is σ x = σ y = ±0.667, τ xy = ±0.200 In the present triangular shell element formulation, vectorial rotational variables are defined. These can be calculated from the prescribed rotations, The rectangular plate is meshed into triangular elements (Fig. 3), and the linear solution obtained by using 10 CR3T elements agrees exactly with the theoretical solution [55].

Laminated Cylindrical Shell
Laminated cylindrical shells subjected to a point load at the central point B are studied. The shell geometry is shown in Fig. 4, where two straight edges are simply supported, length L = 508 mm, radius R = 2540 mm, central angle 2β = 0.2, and thickness t = 6.35 or 12.7 mm.
We consider two groups of lamination schemes consisting of 12 layers and 48 layers, respectively. The first group includes four 12-layer lamination shells arranged as (0 , respectively, where the subscript "s" indicates that the 48 laminas are arranged symmetrically with respect to the shell mid-surface, and the number "3" or "6" before the letter "s" denote that 3 or 6 groups of laminae arranged as those in the parentheses. The material properties of these laminated cylindrical shells are given in Tab. 1.  Due to symmetry, only a quarter of these laminated cylindrical shells are studied by using 8 × 8 × 2 or 16 × 16 × 2 CR3T elements. The load-displacement curves obtained by the present CR3T elements using two different nodal numbering schemes are depicted in Figs. 5-11, which agree very well with the results obtained by Gal and Levy using triangular shell elements built from the linear membrane constant strain triangle and the DKT flat triangular plate element [56] and by Li et al. [36] using 9-node co-rotational quadrilateral composite shell elements (SQCS element). As can be clearly seen from these simulation results, varying the lamination schemes and shell thicknesses can lead to dramatically different structural responses of composite shells.

Stiffened Doubly Curved Cylindrical Panel
A stiffened doubly curved cylindrical panel is subjected to a lateral force at the midpoint of the free curved edge (Fig. 16).

Laminated Sickle-Shell Problem
A laminated cantilever sickle shell is subjected to a lateral force at the free end, as shown in Fig [58] using respectively two quadrature elements (with 10 integration points along the length and 7 integration points along the width) incorporating thickness stretch and drilling rotations and 256 × 10S4R elements of ABAQUS [59] are also reported in these figures. It can be seen that the results of the present triangular shell finite element formulation are in good agreement with them, although slight difference exists because the thickness stretch is neglected in the present shell formulation.

Conclusions
A 3-node co-rotational triangular composite shell element for large deformation analysis of smooth, folded and multi-shell laminated composite structures is proposed. Different from other existing elements using traditional rotational variables, vectorial rotational variables are employed in the present element under a co-rotational framework. All nodal variables are additive in the nonlinear solution procedure, and the global tangent stiffness matrix is symmetric, which enhances the computational efficiency and saves computer storage resource. To overcome shear locking phenomenon, the conforming transverse shear strains are replaced with assumed transverse shear strains by using the line integration method. The computational performance of the developed finite element formulation is demonstrated through solving several smooth and non-smooth laminated composite shell structural problems.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix A. First-and second-order derivatives of strains with respect to local nodal variables
The first-order derivatives of membrane strains with respect to local nodal variables: where is the component at the jth row and the kth column of the inversed Jacobian matrix. N i.ξ and N i.η are the first-order derivatives of the shape function of node i with respect to ξ and η.
The first-order derivatives of shear strains with respect to local nodal variables: The first-order derivatives of bending strain with respect to local nodal variables:

Appendix B. Sub-matrices of the transformation matrix T and its first-order derivatives with respect to global nodal variables
Sub-matrices of the transformation matrix T are give as follows: where, Case 1: If Node i is in smooth shells or away from non-smooth shell intersections: In Eq. (B8), where p k,X , p k,Y , p k,Z are three components of the shell director p k along the directions of global coordinate axes; p k,n , p k,m are two vectorial rotational variables of Node k, which are the two smallest components among p k,X , p k,Y , p k,Z .
where c 0 = 1 − e 2 iy,n − e 2 iz,n , ∂c 0 ∂e iy,n = −e iy,n c 0 , ∂c 0 ∂e iz,n = −e iz,n c 0 (B30-32) The first-order derivatives of the transformation matrix T with respect to the global nodal variables are as follows: In Eqs. (B33)-(B35), Case 1: If Node i is in smooth shells or away from non-smooth shell intersections,