Abstract

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 -continuity at the interface between adjacent elements as well as -continuity 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.

1. 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-to-weight 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 [26], two-dimensional 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 [911].

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 [1214] 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, 1517] 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 [1821] 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 [2327] 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 polynomial-refinement is adopted, highly reliable locking free numerical solutions may be obtained in a purely displacement-based 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.

2. Formulation of Cylindrical Shell Element with Hierarchical Shape Function

2.1. Geometry and Displacement Fields

In cylindrical coordinate shown in Figure 1, based on function theory with continuity at the interface between adjacent layers, dimensional reduction is carried out by incorporating the first-order 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 element such as three sets of nodal translation components (, , 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 asIn (1), indices and refer to nodal and modal contributions, and 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.

2.2. Hierarchical Shape Functions

For definition of the two-dimensional shape functions, firstly one-dimensional hierarchical shape functions (, , and ) on the basis of standard coordinates are presented. The first two of these shape functions in the -direction can be defined asThe 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 two-dimensional 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

2.3. 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 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 formHere 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 bywhere 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. and 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.

2.4. Finite Element Formulation

The displacement fields defined by (1) can be represented in the following general form: 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 workThe internal virtual strain iswith and being the strain and stress tensors used in (8).

If the virtual displacements are defined as the virtual strains can be written aswhere 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 on the basis of the local coordinate system is obtained by the following transformation from material axes to local axes using the transformation matrix :

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

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

3.2. Isotropic Cylindrical Shell Panels with Uniformly Distributed Load

Next, an isotropic ( and ) cylindrical shell panel is considered. Figure 3 shows geometry shape of the panel and loading condition () of which specific values are as follows: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 . 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  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 () and thus the effect of selective reduced integration is worthless.

3.3. Pinched Cylinder with Rigid End Diaphragms

This is another well-known benchmark problem. Figure 5 shows a cylinder with rigid end diaphragms subjected to two point loads . Geometry data in Figure 5 is given byFor material conditions, isotropic material ( and ) 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 from [32] which used the shell theory without transverse shear deformation. Cho and Roh [33] proposed by including the transverse shear deformation. All results of present analysis and the conventional finite elements approach to the reference value 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  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 means the deflection at the point A and refers to the displacement in the direction at the point B. 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. 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 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.

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

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of the paper.

Acknowledgment

This work was supported by National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST) (NRF-2013R1A1A2057756).