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

Zhongxue Li,Jiawei Ji,Loc Vu-Quoc,Bassam A.Izzuddin and Xin Zhuo

1Department of Civil Engineering,Zhejiang University,Hangzhou,310058,China

2Aerospace Engineering,University of Illinois at Urbana-Champaign,Urbana,IL 61801,USA

3Department of Civil and Environmental Engineering,Imperial College London,London,SW7 2BU,UK

ABSTRACT 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.

KEYWORDS Co-rotational approach; 3-node triangular finite element; laminated composite shells; folded and multi-shell structures; vectorial rotational variable; line integration approach; large deformation analysis

1 Introduction

Laminated composite shell structures are extensively used in pressure vessels, aircraft and spacecraft, automotive and other industries due to their high strength- and stiffness-to-mass ratios,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–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 C1weak-continuity.Tran et al.[16,17] proposed an edge-based smoothing discrete shear gap method using 3-node triangular elements combined with a C0-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 C0-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′c 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–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.

2 Co-Rotational Framework and Kinematics of the 3-Node Shell Finite Element

2.1 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 v120and v130:

Figure 1: Illustration of the co-rotational framework

where Xi0(i= 1, 2, 3)is the position of Nodeiin 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 v12and v13are defined as

where di(i= 1, 2, 3)is the translational displacement vector of Nodeiin 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:

wheres4=±1 takes the same sign as that ofpi,lin the previous loading step, and {n,m,l} denotes the circular permutation of {X,Y,Z}.

On the other hand, if Nodeiis located on the intersections of non-smooth shells, the vectorial rotational DOFs of Nodeiis defined asHere,eiy,nandeiy,mare the two smaller components of eiy, andeiz,nis the smallest or second smallest component of eiz, where eiyand eizare two orientation vectors of an orthogonal triad at Nodei.At each incremental loading step, the following relationships can be derived:

The remaining components of the vectors eiyand eizcan be calculated by usingeiy,n,eiy,mandeiz,naccording to Eqs.(7a)–(7c).The sign symbols(s1,s3)defined as 1 or ?1 are respectively equal to the signs of componentseiy,landeiz,mof the last incremental step, and we haves2=?s1·s3.The vector eixis calculated as the cross-product of eiyand eiz:

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

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:

where

For convenience, eix0,eiy0and eiz0in the initial configuration are chosen to be coincident with the base vectors of the global coordinate system.

2.2 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 Nodeiis 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 vectorzlχand transverse shear strain vectorγ.

where,zlis the coordinate of the material point along the shell thickness direction,ais the thickness of the element, and

3 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

whereWeis 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,μ12andμ21are the Poisson’s ratios, E1and E2are the elastic moduli, G12, G23and G31are the shear moduli [1,14].

By enforcing the stationarity condition to the potential energy, we have

which yields the element internal force vector:

where uL, fextand fintare element nodal displacement vector, external load vector and internal force vector in the co-rotational local coordinate system, respectively.Bm, Bband 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 Deq1, Deq2and Deq3are the equivalent elastic matrices:

ndenotes the number of laminae;ζiis the dimensionless coordinate of the interface between Layeri?1 and Layerirelative to the mid-surface:

Dciis the elastic matrix of theith layer lamina:

where Ciis the transversely isotropic elastic matrix of Layeriin the material coordinate system,which is obtained by replacing the material properties of Eq.(25)with those of theith layer; Tiis the transformation matrix between the material coordinate system of theith layer lamina and the co-rotational element local coordinate system:

where?iis defined as the ply orientation angle measured between the reinforcement fiber of theith 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:

4 Stiffness Matrix of the 3-Node Shell Finite Element in the Global Coordinate System

The internal force fGof the 3-node triangular composite shell element in the global coordinate system can be obtained from the local element internal force fintas 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 KTGin 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.

5 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–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]), two alternative nodal ordering schemes are considered in 7 examples of Examples 5.2, where triangular elements used for a quadrilateral patch are employed with different nodal numbering as illustrated in Fig.2.

Figure 2: Two nodal numbering schemes, (a)Scheme 1; (b)Scheme 2

5.1 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 lengthL=0.24, widthW=0.12 and thicknessh=0.001, Young’s modulusE=106, and Poisson’s ratioμ=0.25, is considered.

In the membrane patch test, the displacementsu,υ,wat the boundary nodes of the rectangular plate (Fig.3)are prescribed byu= 10?3(X + Y/2),υ= 10?3(Y +(X/2))andw= 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.

Figure 3: Patch tests for in-plane membrane/out-of-plane bending plates

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].

5.2 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 ()and(, respectively.Here, 0?denotes the circumferential direction of the cylinder, and 90?denotes the cylinder’s axial direction.For instance, (/...)indicates that there are 4 laminas at 90 degrees followed by 4 laminas at 0 degrees, etc.).The second group consists of three 48-layer laminations arranged asand (, 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.

Figure 4: Laminated cylindrical shell

Table 1: Material properties of the laminated cylindrical shells

Figure 5: Load-displacement curves of the laminated cylindrical shell ()with t =6.35 mm

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.

Figure 6: Load-displacement curves of the laminated cylindrical shell ()with t = 6.35 mm

Figure 7: Load-displacement curves of the laminated cylindrical shell ( with t =12.7 mm

5.3 A Laminated Channel Section Cantilever

A laminated channel section cantilever is subjected to a concentrated load at the upper corner point, as shown in Fig.12, whereL= 36 in,a= 2 in,b= 6 in, h = 0.06 in, andP(f)=fPrefwith Pref= 100 lb.Two cross-ply lamination schemes are considered, including (0?/90?/0?)and(90?/0?/90?), where the 0?and 90?directions are respectively along the t1-axis and the t2-axis of the local material coordinate system depicted in Fig.12.The composite material parameters areE1= 107lb/in2,E2= 4×105lb/in2,G12=G13= 2×105lb/in2,G23= 8×104lb/in2,μ12=μ21=0.333.

Figure 8: Load-displacement curves of the laminated cylindrical shell ( with t = 12.7 mm

Figure 9: Load-displacement curves of the laminated cylindrical shell S with t =6.35 mm

The laminated channel section cantilever is modeled using respectively (4 + 4 + 6)× 72 ×2 and (6 + 6 + 9)× 108 × 2 CR3T element meshes, where the web is discretized respectively by 6×72×2 and 9×108×2 elements, and the upper and the lower flanges are discretized by 4×72×2 and 6×108×2 elements, respectively.The load-deflection curves at the upper corner point are presented in Fig.13.For comparison, the results from Chró′scielewski et al.[57] using(4 + 4 + 6)× 72 CAMe16 elements (Lagrange family of 16-node displacement-rotation-based shell elements)are also plotted in this figure.

Figure 10: Load-displacement curves of the laminated cylindrical shell ()3S with t= 6.35 mm

Figure 11: Load-displacement curves of the laminated cylindrical shell ()6S with t =6.35 mm

Figure 12: Geometry and loading of the three-layer channel section cantilever

Figure 13: Load-displacement curves of the three-layer channel section cantilever

The deformed shapes of the channel section cantilever with different lamination schemes are presented in Figs.14 and 15, where the vertical displacements at the upper corner of two cantilevers are up to 5.0 in.

Figure 14: Deformed shape of the laminated channel section cantilever (0?/90?/0?)

Figure 15: Deformed shape of the laminated channel section cantilever (90?/0?/90?)

5.4 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).The geometry is described byL= 2 m,a= 45?,R= 1 m,H= 0.4 m,h0= 0.01 m,Pref= 1 MN,P(λ)=λPref.Two lamination schemes are considered as (0?/90?/90?/0?)and (90?/0?/0?/90?), where the 0?and 90?directions are respectively along the t1-axis and t2-axis of local material coordinate system depicted in Fig.16.The material properties areE1= 105MPa,E2= 7×103MPa,G12=G13= 4×103MPa,G23= 3.2×103MPa,μ12=μ21= 0.25.

Figure 16: Geometry and loading of the stiffened doubly curved cylindrical panel

The stiffened doubly curved cylindrical panel with lamination scheme (0?/90?/90?/0?)is modeled using (30 + 30 + 14)×40×2 CR3T element meshes, where the flat panel is discretized by 14×40×2 elements, and two pieces of cylindrical panels are discretized by 30×40×2 elements,respectively.The stiffened doubly curved cylindrical panel with lamination scheme (90?/0?/0?/90?)is modeled using (25 + 25 + 12)×40×2 CR3T element meshes, where the flat panel is discretized by 12×40×2 elements, and two pieces of cylindrical panels are discretized by 25×40×2 elements,respectively.The load-deflection curves at Point A of the free curved edge are plotted in Fig.17.For comparison, the results from Chró′scielewski et al.[57] using (10 + 10 + 4)×12 CAMe16 elements are also reported in this figure.

Figure 17: Load-displacement curves of the stiffened doubly curved cylindrical panel

The deformed shapes of the stiffened doubly curved cylindrical panel are presented in Figs.18 and 19, which correspond to the states when the displacement of Point A is up to WA= 0.2167 m for lamination scheme (0?/90?/90?/0?)and WA= 0.1692 m for lamination scheme (90?/0?/0?/90?)in Fig.17, respectively.

Figure 18: Deformed shape of the stiffened doubly curved cylindrical panel (0?/90?/90?/0?)

Figure 19: Deformed shape of the stiffened doubly curved cylindrical panel (90?/0?/0?/90?)

5.5 Laminated Sickle-Shell Problem

A laminated cantilever sickle shell is subjected to a lateral force at the free end, as shown in Fig.20, where lengthL= 10, widthB= 1, thicknessa= 0.01, and the radius of the semicylinderR= 5.Two lamination schemes are considered as (0?/90?/0?)and (90?/0?/90?), where the 0?direction is along the flat panel’s longitudinal direction and the semi-cylinder panel’s circumferential direction, and the 90?direction is along the transverse direction.The dimensionless material properties areE1= 6×107,E2= 3×107,G12=G13= 1.49×107,G23= 1.15×107,μ12=μ21= 0.3.

The cantilever sickle shell is modeled using (40 + 40)×4×2 CR3T elements, where the flat panel and the semi-cylinder are discretized by 40×4×2 elements, respectively.The load-deflection curves at the midpoint of the free end are plotted in Figs.21 and 22 for the two lamination schemes.For comparison, the results from Zhang et al.[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.

Figure 20: Cantilever sickle shell subject to a lateral tip load

Figure 21: Load-displacement curves of the cantilever sickle shell (0?/90?/0?)

The deformed shapes of the sickle shell at different levels of the lateral tip load are presented in Figs.23 and 24, showing large deflections and large rotations of the composite shell structure.

Figure 22: Load-displacement curves of the cantilever sickle shell (90?/0?/90?)

Figure 23: Deformed shapes of the cantilever sickle shell (0?/90?/0?)

Figure 24: Deformed shapes of the cantilever sickle shell (90?/0?/90?)

6 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.

Funding Statement:This work was supported by National Natural Science Foundation of China under Grant 11672266.

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

(j,k= 1, 2)is the component at thejth row and thekth column of the inversed Jacobian matrix.Ni.ξandNi.ηare the first-order derivatives of the shape function of nodeiwith 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 Nodeiis in smooth shells or away from non-smooth shell intersections:

In Eq.(B8),

wherepk,X,pk,Y,pk,Zare three components of the shell director pkalong the directions of global coordinate axes;pk,n,pk,mare two vectorial rotational variables of Nodek, which are the two smallest components amongpk,X,pk,Y,pk,Z.

Case 2: If Nodeiis located at an intersection of non-smooth shells,

In Eq.(B15),

In Eqs.(B16)–(B18),

where

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 Nodeiis in smooth shells or away from non-smooth shell intersections,

In Eqs.(B36)–(B39),

Case 2: If Nodeiis located at an intersection of non-smooth shells,

where