ESL Based Cylindrical Shell Elements with Hierarchical Shape Functions for Laminated Composite Shells

We introduce higher-order cylindrical shell element based on ESL (equivalent single-layer) theory for the analysis of laminated composite shells. The proposed elements are formulated by the dimensional reduction technique from three-dimensional solid to two-dimensional cylindrical surface with plane stress assumption. It allows the first-order shear deformation and considers anisotropic materials due to fiber orientation. The element displacement approximation is established by the integrals of Legendre polynomials with hierarchical concept to ensure the C-continuity at the interface between adjacent elements as well as Ccontinuity at the interface between adjacent layers. For geometry mapping, cylindrical coordinate is adopted to implement the exact mapping of curved shell configuration with a constant curvature with respect to any direction in the plane. The verification and characteristics of the proposed element are investigated through the analyses of three cylindrical shell problems with different shapes, loadings, and boundary conditions.


Introduction
Shell structures are three-dimensional structures with any curvature, thin in one direction and long in the other two directions.In engineering design, they are among the most significant and ubiquitous structural components.Applications of them include pressure vessels, the bodies of automobiles and airplanes, bridges, buildings, roofs, the hulls of ships and submarines, and many other structures.Particularly, increasing application of laminated composite shell is evident in a variety of engineering structures and manufactured components, because of the well-recognized characteristics of superior strength-to-weight, stiffness-to-weight, and cost-toweight ratios, compared to conventional materials.While the laminated composite materials provide the design flexibility to achieve desirable stiffness and strength through the choice of lamination scheme, the anisotropic constitution of laminated composite structures often results in stress concentrations near material and geometric discontinuities that can lead to damage in the form of delamination, adhesive bond separation, and matrix cracking.Recently, these problems have been mitigated by replacing conventionally used laminated composites with functionally graded materials where the materials properties are gradually varied at microscopic scale in the thickness direction [1].
Finite element methods are versatile numerical tools to solve differential equations related to physical phenomena.In the finite element applications for shell analysis, some types of shell elements are currently available.They are flat facet element, shell theory-based element, degenerated shell element, and solid-shell element and so forth.For the flat facet element, it does not have any curvature.Thus, the curved surface is approximately explained by the combination of several elements.It is very simple in formulations and has still been used for engineering applications.However, it cannot explain bending-stretching coupling behavior in an element level.On the other hand, the shell theory-based element with curvatures can handle the bending-stretching coupling properly.Also, the degenerated shell element can be used for arbitrary shapes of shell surface.Based on the flat facet element or shell theory-based element [2][3][4][5][6], twodimensional shell elements have been introduced.Acknowledging the need of three-dimensional shell elements, several formulations have been presented based on a degenerated shell concept [7,8].Since 1990s, some solid-shell approaches, which have some benefits as compared to the degenerated shell element because of the simplicity of their kinematics, have been introduced by some researchers [9][10][11].
Meanwhile, a number of innovative approaches have been put forward for the analysis of laminated structural systems, to extend the capabilities of laminated anisotropic composites.As far as two-dimensional modeling is concerned, it is assumed that displacement components are continuously differentiable through the thickness regardless of the layer boundaries.Representatives of the theories are known as classical lamination theory (CLT) and first-order shear deformation theory (FSDT).Both of these models [12][13][14] are known as equivalent single-layer theories (ESLT) based on certain assumptions concerning the kinematics of deformation or stress across the total thickness.Although FSDT provides a sufficiently accurate description of the global laminate response for thin to moderately thick plates, it cannot allow direct calculation of transverse stresses with acceptable accuracy.So a number of higher-order theories [12,[15][16][17] have been put forward using successively third-to higher-degree polynomials and other functions with continuous derivatives to yield more accurate interlaminar stress distributions.The deficiency of the theories has led to layerwise models in which the variation of displacement functions across the thickness is assumed for each layer separately.The layerwise models [18][19][20][21] require displacement continuity at layer interfaces.Such characterization of laminated systems can generally exhibit a rapid change of slopes of displacement fields at layer interfaces, often termed as the zig-zag effect.In order to satisfy the interlaminar continuity of transverse stresses at each layer interface, appropriate functional continuities are required for transverse displacements and stresses [22].The number of modal degrees of freedom in normal layerwise models depends on the number of layers in the laminated system.In conventional finite element analysis based mostly on Lagrangian two-dimensional shape functions, the layerwise models can satisfy displacement continuity but not stress.It is thus true that normal layerwise models would be too expensive when it is intended to comply with transverse normal stress continuity.Thus, multiple model approaches [23][24][25][26][27] have also been attempted to reduce the overall number of modal degrees of freedom by optimizing computation process for maximum solution accuracy within a particular subregion of interest only and in the process reducing the computational effort.
It is well-known that low-order finite element implementation for shells suffers from various forms of locking whenever purely displacement-based formulations are employed.In recent years, the issue of locking has been most prominently addressed through the use of low-order finite technology using mixed variational principles.The assumed strain and enhanced strain formulations are among the successful low-order implementations.High-order finite element implementations have also been advocated in recent years as a means of eliminating the locking phenomena completely.Most notably, whenever a sufficient degree of polynomialrefinement is adopted, highly reliable locking free numerical solutions may be obtained in a purely displacementbased setting [28].The first -version formulation related to shells, one of high-order approaches, was reported by Woo and Basu [29] who presented a cylindrical shell element formulation in the cylindrical coordinates associated with a suitable transfinite mapping function to represent the curved geometry.In this paper, we address the finite element formulation for the laminated cylindrical shell behavior using the -version approach.The approach assumes that a heterogeneous laminated shell stacked with several laminae is treated as a shell element using hierarchic interpolation functions.Thus, characteristics of the proposed approach are presented in detail.Since higher-order Lagrange shape functions cannot be used due to excessive round-off errors, all approximate functions for displacement fields are derived in terms of integrals of Legendre polynomials which are orthogonal in the energy norm.

Formulation of Cylindrical Shell Element with Hierarchical Shape Function
2.1.Geometry and Displacement Fields.In cylindrical coordinate shown in Figure 1, based on  1  function theory with continuity at the interface between adjacent layers, dimensional reduction is carried out by incorporating the firstorder shear deformation for bending behavior and the plane stress condition for membrane action.For geometry and displacement fields, the curvilinear coordinate system is considered in reference to the middle surface of laminated shells keeping a constant curvature  with respect to any direction in the two-dimensional plane.Geometry fields on a surface defined by two axes are expressed by linear interpolation between  and  variables over the four vertex nodes only, as shown in Figure 1.
Also, the deformation at any point in the laminated shell is based on three displacement fields for a quadrilateral subparametric  0  element such as three sets of nodal translation components (  , V  , and   ) and modal translation components (  ,   , and   ), and two sets of nodal rotation components (  and   ) and modal rotation components (  and   ).The three displacement fields can be defined as In ( 1), indices  and  refer to nodal and modal contributions, and n is the number of modal variables which is zero when degree of polynomial approximation is one.In addition,   and   refer to two-dimensional shape functions in terms of standard coordinates (, ) associated with nodal and modal variables which are defined in next section.

Hierarchical Shape Functions.
For definition of the twodimensional shape functions, firstly one-dimensional hierarchical shape functions ( 1 ,  2 , and   ) on the basis of standard coordinates are presented.The first two of these shape functions in the -direction can be defined as The shape functions corresponding to modal variables are defined in terms of integrals of Legendre polynomials, as shown below [30]: For and -directions, analogous expressions are obtained by replacing  by  and , respectively, in (2) and (3).One-dimensional shape functions are used to construct twodimensional shape functions,   and   , in the , -plane.In other words, two-dimensional shape functions at four corner nodes can be obtained by a product of two one-dimensional nodal functions in the following: Next, the modal variables of two-dimensional shape functions of laminated systems are classified into two groups such as "side modes" and "internal modes" as shown in Figure 1.
The number of modal variables  is dependent on the degree of polynomial (-level), as shown in the following: The shape functions for side modal variables for any -level are shown in (6).The superscripts in (6) refer to side numbers: The shape functions for internal modal variables are obtained by

Strain and Stress
Relations.The three strain components can be denoted by membrane (), bending (), and shear () strain, respectively: where The constitutive relationship for the laminate with respect to the reference surface can then be expressed as where [] 8×8 is the constitutive matrix including membrane, bending, and transverse shear force resultants to reference surface strains and curvatures.The stress resultant vector {σ} will be of the form Here letters , , and  refer to resultants of membrane stresses, bending stresses, and transverse shear stresses.Based on (10) and (11), the stress resultants can be defined by where elasticity matrix [] represents anisotropy with three mutually orthogonal planes of symmetry.The elasticity matrix includes shear correction factors  in order to allow for the error resulting from the use of transverse shear strain energy on an average basis which depends on lamina properties and the lamination scheme. bottom  and  top  are distances from the reference surface to bottom and top surfaces of lamina.ℓ and  are the number of laminas making up the thickness of laminated shells.
where the matrix [ Ĥ] indicates hierarchical shape functions for nodal and modal variables û with total number of variables denoted by .The finite element equation for each model can be expressed by using the principle of virtual work The internal virtual strain is with {ε} and {σ} being the strain and stress tensors used in (8).
If the virtual displacements are defined as the virtual strains can be written as where [ B] is the strain-displacement matrix.The external virtual work takes the following form: Here the superscripts , , , and  refer to nodal and modal forces, side forces, surface forces, and body forces, respectively.Based on these definitions, the virtual work equation shown in ( 14) can be rewritten as Here constitutive matrix [ D] on the basis of the local coordinate system is obtained by the following transformation from material axes to local axes using the transformation matrix []:

Numerical Examples
The performance of proposed cylindrical shell element with hierarchical shape function is investigated and compared with results obtained by some conventional finite element methods [31] via the following numerical examples in which all units of parameters are expressed by nondimensional values.Three types of conventional finite elements are considered such as 4-node element (4N) with linear Lagrangian polynomials, 8-node element (8N) often called serendipity element, and 9-node element (9N) with quadratic Lagrangian polynomials.Also all conventional finite elements adopt different integration techniques, separately, like full (F), reduced (R), and selective reduced (S) integration.For present analysis and conventional finite elements, identical linear mapping based on cylindrical coordinate is implemented.For convergence, present analysis is implemented based on dominantly -refinement, occasionally ℎ-refinement, while conventional finite elements choose only ℎ-refinement in which all calculated values have four or five significant digits.

Orthotropic Cylindrical Shells with Internal Pressure.
Firstly, displacements of cylindrical shells with one layer (0 ∘ ) and two layers (0 ∘ /90 ∘ from inner surface) subject to internal pressure  (=6410/) as shown in Figure 2 are estimated.
In the problem with two layers, thickness of two layers is identical.Length  of the shell is 20, the radius  is 20, and the total thickness  is 1.Each layer is a unidirectional fiber reinforced composite with the following material constants: where subscript 1 signifies the direction parallel to fibers and 2 and 3 the transverse direction.In these problems, the number 1 means  direction, the number 2 is  direction, and the number 3 refers to  direction, respectively.Also, clamped boundary conditions are specified at both ends of the shell and only octant of the shell is considered to take advantage of three-way symmetry.Table 1 gives the maximum displacement results in the cylindrical shell with one layer.All converged values are almost identical except those of 4-node element case.Present analysis using only one element has the converged value from -level = 4, while 8-node and 9-node elements require 8 × 8 mesh.Table 2 shows the maximum displacement results with two layers.It is noted that similar characteristics are obtained.Also, Tables 1 and 2 show the converged values shaded according to analysis types.A quarter of panel is considered as a computational region shown in Figure 3 by specifying symmetry conditions with respect to  and  axes.At  = , clamped conditions are applied and free boundary is defined at  = 0. Figures 4(a) and 4(b) show the variation of maximum vertical deflection with respect to the number of degrees of freedom (NDF) obtained by all elements considered.The converged values are all identical except 4NF which requires more fine mesh refinement to obtain a converged value.For the conventional finite elements, it is seen that convergence rates with full integration are slower than those elements with reduced or selective reduced integration because of shear locking problem.Present analysis with full integration gives much faster convergence than the other elements to keep same converged value.Table 3 shows specific values with respect to analysis types.For present analysis, the maximum displacement is converged from -level = 6 or 7 with 2 × 2 mesh.Also, for the conventional finite elements with 8N and 9N, there is little difference according to integration techniques.It is noticed that the shear behavior is not significant as the thickness ratio becomes very thin (/ = 160) and thus the effect of selective reduced integration is worthless.

Pinched Cylinder with Rigid
End Diaphragms.This is another well-known benchmark problem. Figure 5 shows For material conditions, isotropic material ( = 3 × 10 6 and ] = 0.3) is firstly considered.For a numerical analysis, only octant of the entire shell is modelled taking advantage of three-way symmetry as shown in Figure 5.This problem is known as one of bench mark problems with bending-dominant behavior in the thin limit.Thus transverse shear locking and membrane locking occur as the shell is getting thinner.Table 4 shows the maximum radial displacement at the point of load application.The analytical  solution is −1.8248 × 10 −5 from [32] which used the shell theory without transverse shear deformation.Cho and Roh [33] proposed −1.8541 × 10 −5 by including the transverse shear deformation.All results of present analysis and the conventional finite elements approach to the reference value (−1.8541 × 10 −5 ) as the number of elements is increased.
For present analysis, it is difficult to obtain the converged value when only -refinement is implemented.Also, it is seen in Table 4 that the converged solution of present analysis requires ℎand -refinement simultaneously.For -level ≥ 7 in the 8×8 mesh, the relative error of maximum displacement is less than 1.0% as compared with the reference value.Next, the cylindrical laminated shell with two layers is analyzed.The thickness of two layers is identical and total thickness  of the cylinder is fixed as 3.0.Each layer consists of a unidirectional fiber reinforced composite with following material constants: where subscript  signifies the direction parallel to the fibers,  is the transverse direction, and ]  is the Poisson's ratio for strain in the  direction under uniaxial normal stress in the  direction.  is an elastic modulus in the transverse direction and  is an elastic shear modulus.The normalized quantities of interest are defined as follows: where  A means the deflection at the point A and  B refers to the displacement in the  direction at the point B.    convergence rates appear for those values of present analysis.Also, it is seen that the converged values of 8NR, 8NS, and 9NS and present analysis gives good agreement.

Conclusions
The aim of this study is to show the efficiency of the proposed cylindrical shell element with hierarchical shape function for the analysis of laminated shells with isotropic and orthotropic materials.Based on this study, the following main conclusions can be drawn.
(1) Numerical examples demonstrate that the proposed shell element shows a rapid convergence rate and a stable numerical stability regardless of geometry, loading, and boundary condition as compared with the conventional finite element approaches.
(2) The proposed element can endure very thin thickness ratio without any special numerical integration technique.It may be concluded that membrane and shear locking are considerably free over -level = 4 or 5.
(3) From these results, it is necessary to develop the hierarchical shell elements with arbitrary geometry by applying the advanced mapping technique.
Finally, although this study deals with classical shell problems, it is important to note that the proposed shell elements based on high-order approaches can be extended to applications of various shell problems with cutout and functionally graded materials which many researchers have recently been interested in.

Figure 1 :
Figure 1: Geometric configuration and coordinate system of the proposed element.

2. 4 .
Finite Element Formulation.The displacement fields defined by (1) can be represented in the following general form:
Panels with Uniformly Distributed Load.Next, an isotropic ( = 0.45×106 and ] = 0.3) cylindrical shell panel is considered.Figure3shows geometry shape of the panel and loading condition ( = 0.04) of which specific values are as follows: Convergence characteristics by reduced and selective reduced integration

Figure 4 :
Figure 4: Convergence of maximum displacement with respect to NDF.

Figures 6 (
a) and 6(b) show the convergence characteristics of  and  displacements, respectively.It is seen that the converged values of all elements are almost same.

Figure 6 :
Figure 6: Convergence of displacements in the laminated shell.

Table 1 :
Maximum displacements of the cylindrical shell with one layer.

Table 2 :
Maximum deflections of the cylindrical laminate shell with two layers.

Table 5
contains the normalized deflections at the point A. The values of present analysis are close to 1.27 by ℎ-refinement, while the values of 8NR/S and 9NR/S reach 1.24 and 1.31.The difference between present solutions and ℎ-refinement analyses is within 2∼3%.Table 6 shows displacements in  direction at point B. Unlike deflection values at point A, high

Table 5 :
Maximum radial deflection at the point A ().

Table 6 :
Displacement in  direction at the point B ().