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

Authors: Jae S. Ahn; Seung H. Yang; Kwang S. Woo
Journal: Mathematical Problems in Engineering
 (2015)

---

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


---

## Body

## 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 [2–6], 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 [9–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–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–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–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–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 polynomial-refinement is adopted, highly reliable locking free numerical solutions may be obtained in a purely displacement-based setting [28]. The first p-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 p-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 Figure1, based on Cr1 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 x and θ variables over the four vertex nodes only, as shown in Figure 1.Figure 1
Geometric configuration and coordinate system of the proposed element.Also, the deformation at any point in the laminated shell is based on three displacement fields for a quadrilateral subparametricCxθ0 element such as three sets of nodal translation components (ui, vi, and wi) and modal translation components (αj, βj, and χj), and two sets of nodal rotation components (ρi and λi) and modal rotation components (φj and ψj). The three displacement fields can be defined as(1)ux,θ,r=Siui+zSiρi+Bjαj+rBjϕj,vx,θ,r=Sivi+zSiλi+Bjβj+rBjψj,wx,θ,r=Siwi+Bjχj,i=1,2,3,4;j=5,6,…,n^+4.In (1), indices i and j 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, Si and Bj 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 (L1, L2, and Lj) on the basis of standard coordinates are presented. The first two of these shape functions in the ξ-direction can be defined as(2)L1ξ=0.51-ξ,L2ξ=0.51+ξ.The shape functions corresponding to modal variables are defined in terms of integrals of Legendre polynomials, as shown below [30]:(3)Li+1ξ=2i-12∫-1ξ12i-1i-1!di-1dwi-1w2-1i-1dwi=2,3,4,….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, Si and Bj, 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: (4)S1ξ,η=L1ξN1η,S2ξ,η=L2ξN1η,S3ξ,η=L2ξN2η,S4ξ,η=L1ξN2η.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 f is dependent on the degree of polynomial (p-level), as shown in the following: (5)fsidep=4p-1,finternal(p)=0inp≤3p-2p-32inp≥4.The shape functions for side modal variables for any p-level are shown in (6). The superscripts in (6) refer to side numbers: (6)Bi1ξ,η=L1ηLi+1ξ,Bi2ξ,η=L2ξLi+1η,Bi3ξ,η=L2ηLi+1ξ,Bi4ξ,η=L1ξLi+1η,with2≤i≤p.The shape functions for internal modal variables are obtained by (7)Binternalξ,η=LiξLjηi+j=p+2i,j=3,4,….


### 2.3. Strain and Stress Relations

The three strain components can be denoted by membrane(m), bending (b), and shear (t) strain, respectively:(8)ε^xε^θγ^xθγ^xrγ^θr=ε^xmε^θmγ^xθmγ^xrtγ^θrt+rε^xbε^θbγ^xθb00,where (9)ε^xmε^θmγ^xθm=∂Si∂ξ00∂Bj∂ξ000κ∂Si∂ηκSi0κ∂Bj∂ηκBjκ∂Si∂η∂Si∂ξ0κ∂Bj∂η∂Bj∂ξ0uiviwiαjβjχj,ε^xbε^θbγ^xθb=0∂Si∂ξ0∂Bj∂ξκ∂Si∂η0κ∂Bj∂η0∂Si∂ξκ∂Si∂η∂Bj∂ξκ∂Bj∂ηλiρiψjϕj,γ^xrtγ^θrt=0∂Si∂ξSi00∂Bj∂ξBi0-κSiκ∂Si∂η0Si-κBjκ∂Bj∂η0Bjviwiλiρiβjχjψjϕj.The constitutive relationship for the laminate with respect to the reference surface can then be expressed as (10)σ^=E8×8ε^,where E8×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(11)σ^=N3×1M3×1Q2×1.Here letters N, M, and Q refer to resultants of membrane stresses, bending stresses, and transverse shear stresses. Based on (10) and (11), the stress resultants can be defined by(12)N3×1=∑l=1n∫rlbottomrltopDij5×5lεxm+rεxbεθm+rεθbγxθm+rγxθbdr;i,j=1,2,3,M3×1=∑l=1n∫rlbottomrltopDij5×5lrεxm+r2εxbrεθm+r2εθbrγxθm+r2γxθbdr;i,j=1,2,3,Q2×1=K∑l=1n∫rlbottomrltopDij5×5lγxrtγθrtdri,j=4,5,where elasticity matrix [D] represents anisotropy with three mutually orthogonal planes of symmetry. The elasticity matrix includes shear correction factors K 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. rlbottom and rltop are distances from the reference surface to bottom and top surfaces of lamina. l and n 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: (13)u~i=H^1×niu^n×1ii=1,2,and3,where the matrix H^ indicates hierarchical shape functions for nodal and modal variables u^ with total number of variables denoted by n. The finite element equation for each model can be expressed by using the principle of virtual work(14)δUε-δW=0.The internal virtual strain is(15)δUε=∫Vδε^Tσ^dVwith ε^ and σ^ being the strain and stress tensors used in (8).If the virtual displacements are defined as(16)δu=H^δu^the virtual strains can be written as(17)δε=B^δu^,where B^ is the strain-displacement matrix. The external virtual work takes the following form:(18)δW=δu^TFa+∫Sδu^TFbdS+∫Aδu^TFcdA.Here the superscripts a, b, c, and d 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 (19)∫VδB^TD^B^dVu^=δu^TFa+∫Sδu^TFbdS+∫Aδu^TFcdA+∫Vδu^TFddV.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 X:(20)D^=XTDX.


## 2.1. Geometry and Displacement Fields

In cylindrical coordinate shown in Figure1, based on Cr1 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 x and θ variables over the four vertex nodes only, as shown in Figure 1.Figure 1
Geometric configuration and coordinate system of the proposed element.Also, the deformation at any point in the laminated shell is based on three displacement fields for a quadrilateral subparametricCxθ0 element such as three sets of nodal translation components (ui, vi, and wi) and modal translation components (αj, βj, and χj), and two sets of nodal rotation components (ρi and λi) and modal rotation components (φj and ψj). The three displacement fields can be defined as(1)ux,θ,r=Siui+zSiρi+Bjαj+rBjϕj,vx,θ,r=Sivi+zSiλi+Bjβj+rBjψj,wx,θ,r=Siwi+Bjχj,i=1,2,3,4;j=5,6,…,n^+4.In (1), indices i and j 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, Si and Bj 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 (L1, L2, and Lj) on the basis of standard coordinates are presented. The first two of these shape functions in the ξ-direction can be defined as(2)L1ξ=0.51-ξ,L2ξ=0.51+ξ.The shape functions corresponding to modal variables are defined in terms of integrals of Legendre polynomials, as shown below [30]:(3)Li+1ξ=2i-12∫-1ξ12i-1i-1!di-1dwi-1w2-1i-1dwi=2,3,4,….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, Si and Bj, 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: (4)S1ξ,η=L1ξN1η,S2ξ,η=L2ξN1η,S3ξ,η=L2ξN2η,S4ξ,η=L1ξN2η.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 f is dependent on the degree of polynomial (p-level), as shown in the following: (5)fsidep=4p-1,finternal(p)=0inp≤3p-2p-32inp≥4.The shape functions for side modal variables for any p-level are shown in (6). The superscripts in (6) refer to side numbers: (6)Bi1ξ,η=L1ηLi+1ξ,Bi2ξ,η=L2ξLi+1η,Bi3ξ,η=L2ηLi+1ξ,Bi4ξ,η=L1ξLi+1η,with2≤i≤p.The shape functions for internal modal variables are obtained by (7)Binternalξ,η=LiξLjηi+j=p+2i,j=3,4,….


## 2.3. Strain and Stress Relations

The three strain components can be denoted by membrane(m), bending (b), and shear (t) strain, respectively:(8)ε^xε^θγ^xθγ^xrγ^θr=ε^xmε^θmγ^xθmγ^xrtγ^θrt+rε^xbε^θbγ^xθb00,where (9)ε^xmε^θmγ^xθm=∂Si∂ξ00∂Bj∂ξ000κ∂Si∂ηκSi0κ∂Bj∂ηκBjκ∂Si∂η∂Si∂ξ0κ∂Bj∂η∂Bj∂ξ0uiviwiαjβjχj,ε^xbε^θbγ^xθb=0∂Si∂ξ0∂Bj∂ξκ∂Si∂η0κ∂Bj∂η0∂Si∂ξκ∂Si∂η∂Bj∂ξκ∂Bj∂ηλiρiψjϕj,γ^xrtγ^θrt=0∂Si∂ξSi00∂Bj∂ξBi0-κSiκ∂Si∂η0Si-κBjκ∂Bj∂η0Bjviwiλiρiβjχjψjϕj.The constitutive relationship for the laminate with respect to the reference surface can then be expressed as (10)σ^=E8×8ε^,where E8×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(11)σ^=N3×1M3×1Q2×1.Here letters N, M, and Q refer to resultants of membrane stresses, bending stresses, and transverse shear stresses. Based on (10) and (11), the stress resultants can be defined by(12)N3×1=∑l=1n∫rlbottomrltopDij5×5lεxm+rεxbεθm+rεθbγxθm+rγxθbdr;i,j=1,2,3,M3×1=∑l=1n∫rlbottomrltopDij5×5lrεxm+r2εxbrεθm+r2εθbrγxθm+r2γxθbdr;i,j=1,2,3,Q2×1=K∑l=1n∫rlbottomrltopDij5×5lγxrtγθrtdri,j=4,5,where elasticity matrix [D] represents anisotropy with three mutually orthogonal planes of symmetry. The elasticity matrix includes shear correction factors K 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. rlbottom and rltop are distances from the reference surface to bottom and top surfaces of lamina. l and n 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: (13)u~i=H^1×niu^n×1ii=1,2,and3,where the matrix H^ indicates hierarchical shape functions for nodal and modal variables u^ with total number of variables denoted by n. The finite element equation for each model can be expressed by using the principle of virtual work(14)δUε-δW=0.The internal virtual strain is(15)δUε=∫Vδε^Tσ^dVwith ε^ and σ^ being the strain and stress tensors used in (8).If the virtual displacements are defined as(16)δu=H^δu^the virtual strains can be written as(17)δε=B^δu^,where B^ is the strain-displacement matrix. The external virtual work takes the following form:(18)δW=δu^TFa+∫Sδu^TFbdS+∫Aδu^TFcdA.Here the superscripts a, b, c, and d 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 (19)∫VδB^TD^B^dVu^=δu^TFa+∫Sδu^TFbdS+∫Aδu^TFcdA+∫Vδu^TFddV.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 X:(20)D^=XTDX.


## 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 p-refinement, occasionally h-refinement, while conventional finite elements choose only h-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 pressurep  (=6410/π) as shown in Figure 2 are estimated. In the problem with two layers, thickness of two layers is identical. Length s of the shell is 20, the radius R is 20, and the total thickness t is 1. Each layer is a unidirectional fiber reinforced composite with the following material constants:(21)E1=7.5×106;E2=2×106;G12=1.25×106;G13=G23=0.625×106;ν12=0.25,where subscript 1 signifies the direction parallel to fibers and 2 and 3 the transverse direction. In these problems, the number 1 means x direction, the number 2 is θ direction, and the number 3 refers to r 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.Figure 2
Orthotropic clamped cylindrical shells with internal pressure.Table1 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 p-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.Table 1
Maximum displacements of the cylindrical shell with one layer.
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
Present analysis
p-level = 3
0.3798
—
—
—
—
p-level = 4
0.3746
—
—
—
—
p-level = 5
0.3748
—
—
—
—
p-level = 6
0.3748
—
—
—
—
p-level = 7
0.3748
—
—
—
—
p-level = 8
0.3748
—
—
—
—
p-level = 9
0.3748
—
—
—
—
p-level = 10
0.3748
—
—
—
—
4N
F
0.1478
0.2816
0.3484
0.3681
0.3731
R
0.1574
0.4072
0.3844
0.3772
0.3754
S
0.1574
0.4072
0.3844
0.3772
0.3754
8N
F
0.3220
0.3729
0.3747
0.3748
0.3748
R
0.4179
0.3770
0.3750
0.3748
0.3748
S
0.4179
0.3770
0.3750
0.3748
0.3748
9N
F
0.3129
0.3728
0.3747
0.3748
0.3748
R
0.4179
0.3770
0.3749
0.3748
0.3748
S
0.4179
0.3770
0.3749
0.3748
0.3748Table 2
Maximum deflections of the cylindrical laminate shell with two layers.
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
Present analysis
p-level = 3
0.1860
—
—
—
—
p-level = 4
0.1782
—
—
—
—
p-level = 5
0.1795
—
—
—
—
p-level = 6
0.1794
—
—
—
—
p-level = 7
0.1794
p-level = 8
0.1794
—
—
—
—
p-level = 9
0.1794
—
—
—
—
p-level = 10
0.1794
4N
F
0.1105
0.1705
0.1811
0.1804
0.1797
R
0.1240
0.2290
0.1882
0.1814
0.1799
S
0.1240
0.2290
0.1882
0.1814
0.1799
8N
F
0.1829
0.1830
0.1797
0.1795
0.1794
R
0.2174
0.1796
0.1794
0.1794
—
S
0.2174
0.1796
0.1794
0.1794
—
9N
F
0.1829
0.1830
0.1797
0.1795
0.1794
R
0.2174
0.1796
0.1794
0.1794
—
S
0.2174
0.1796
0.1794
0.1794
—


### 3.2. Isotropic Cylindrical Shell Panels with Uniformly Distributed Load

Next, an isotropic (E=0.45×106 and ν=0.3) cylindrical shell panel is considered. Figure 3 shows geometry shape of the panel and loading condition (p=0.04) of which specific values are as follows:(22)α=0.2;R=100;s=20;t=0.125.A quarter of panel is considered as a computational region shown in Figure 3 by specifying symmetry conditions with respect to x and θ axes. At x=s, 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 p-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 (S/t=160) and thus the effect of selective reduced integration is worthless.Table 3
Maximum vertical deflection (-w×102) at the center of cylindrical panels.
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
Present analysis
p-level =3
0.4903
1.0048
p-level = 4
0.5410
1.1534
p-level = 5
1.2248
1.1403
p-level = 6
1.2525
1.1358
p-level = 7
1.1175
1.1349
p-level = 8
1.1455
1.1349
p-level = 9
1.1320
1.1349
p-level = 10
1.1359
1.1349
4N
F
0.0082
0.0309
0.1075
0.3378
0.7456
R
0.0111
1.7227
1.2376
1.1577
1.1404
S
0.0110
1.7154
1.2400
1.1584
1.1406
8N
F
0.0461
0.7435
1.1745
1.1428
1.1362
R
1.6496
0.9573
1.1375
1.1350
1.1349
S
1.6399
0.9563
1.1374
1.1350
1.1349
9N
F
0.0463
1.3469
1.1721
1.1427
1.1362
R
1.6668
1.1444
1.1352
1.1349
1.1349
S
1.6532
1.1434
1.1352
1.1349
1.1349Figure 3
Clamped cylindrical shell panels with uniform load.Figure 4
Convergence of maximum displacement with respect to NDF.
(a)
Convergence characteristics by full integration
(b)
Convergence characteristics by reduced and selective reduced integration


### 3.3. Pinched Cylinder with Rigid End Diaphragms

This is another well-known benchmark problem. Figure5 shows a cylinder with rigid end diaphragms subjected to two point loads (P=1). Geometry data in Figure 5 is given by(23)S=600;R=300;t=3.0.For material conditions, isotropic material (E=3×106 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.Figure 5
Pinched cylinder with rigid diaphragms.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. Table4 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 p-refinement is implemented. Also, it is seen in Table 4 that the converged solution of present analysis requires h- and p-refinement simultaneously. For p-level ≥ 7 in the 8×8 mesh, the relative error of maximum displacement is less than 1.0% as compared with the reference value.Table 4
Maximum radial deflection (-w×105) at R/t = 100.
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
20 × 20
Present analysis
p-level = 3
0.0674
0.3540
1.2273
1.7090
1.8172
1.8293
p-level = 4
0.1489
0.8821
1.5689
1.7825
1.8350
1.8404
p-level = 5
0.3503
1.3250
1.6937
1.8147
1.8419
1.8452
p-level = 6
0.7655
1.5009
1.7542
1.8289
1.8458
1.8481
p-level = 7
1.1906
1.6180
1.7937
1.8376
1.8485
1.8503
p-level = 8
1.3670
1.7037
1.8162
1.8426
1.8504
1.8520
p-level = 9
1.4789
1.7560
1.8284
1.8457
1.8519
1.8533
p-level = 10
1.5897
1.7883
1.8356
1.8479
1.8532
1.8545
4N
F
0.0098
0.0254
0.0608
0.1282
0.2785
0.3603
R
0.0207
0.3263
1.9281
1.8453
1.8600
1.8632
S
0.0213
0.6723
1.3929
1.6093
1.7762
1.8056
8N
F
0.0369
0.0864
0.3891
1.1772
1.6691
1.7378
R
0.1477
1.0976
1.6444
1.7951
1.8307
1.8365
S
0.1263
0.8565
1.5541
1.7799
1.8282
1.8350
9N
F
0.0414
0.0925
0.4184
1.2238
1.6845
1.7476
R
0.1886
1.7458
1.8451
1.8596
1.8677
1.8702
S
0.1657
1.2428
1.7046
1.8359
1.8633
1.8675Next, the cylindrical laminated shell with two layers is analyzed. The thickness of two layers is identical and total thicknesst of the cylinder is fixed as 3.0. Each layer consists of a unidirectional fiber reinforced composite with following material constants:(24)EL=25ET;GLT=0.5ET;GTT=0.2ET;νLT=νTT=0.25,where subscript L signifies the direction parallel to the fibers, T is the transverse direction, and νLT is the Poisson’s ratio for strain in the T direction under uniaxial normal stress in the L direction. ET is an elastic modulus in the transverse direction and G is an elastic shear modulus. The normalized quantities of interest are defined as follows:(25)w-=-10wAE1t3PR2;u-=-10uBE1t2PR,where wA means the deflection at the point A and uB refers to the displacement in the x direction at the point B. Figures 6(a) and 6(b) show the convergence characteristics of w and u 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 hp-refinement, while the values of 8NR/S and 9NR/S reach 1.24 and 1.31. The difference between present solutions and h-refinement analyses is within 2~3%. Table 6 shows displacements in x 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.Table 5
Maximum radial deflection at the point A (w-).
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
20 × 20
Present analysis
p-level = 3
0.0895
0.2352
0.7834
1.1560
1.2305
1.2393
p-level = 4
0.1366
0.5876
1.0602
1.2086
1.2438
1.2491
p-level = 5
0.3234
0.8951
1.1517
1.2297
1.2518
1.2559
p-level = 6
0.6028
1.0129
1.1894
1.2408
1.2573
1.2610
p-level = 7
0.8076
1.0906
1.2156
1.2486
1.2618
1.2651
p-level = 8
0.8935
1.1485
1.2310
1.2538
1.2654
1.2687
p-level = 9
0.9759
1.1854
1.2404
1.2579
1.2686
1.2718
p-level = 10
1.0577
1.2085
1.2467
1.2611
1.2714
1.2745
4N
F
0.0033
0.0169
0.0539
0.1391
0.3169
0.4038
R
0.0349
0.2760
1.2670
1.2693
1.2915
1.2978
S
0.0351
−0.0530
−0.6145
5.6803
1.4209
1.3725
8N
F
0.0581
0.1218
0.2915
0.7831
1.1546
1.1990
R
0.1446
0.8093
1.1058
1.2021
1.2359
1.2421
S
0.1228
0.7113
1.0710
1.1967
1.2349
1.2414
9N
F
0.0630
0.1314
0.3095
0.8082
1.1669
1.2086
R
0.1705
1.1902
1.2647
1.2878
1.3069
1.3130
S
0.1446
0.7906
1.2049
1.2767
1.3041
1.3112Table 6
Displacement inx direction at the point B (u-).
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16  × 16
20 × 20
Present analysis
p-level = 3
0.1627
0.3403
0.9709
—
—
—
p-level = 4
0.2256
0.9344
1.1661
—
—
—
p-level = 5
0.5338
1.2390
1.1712
—
—
—
p-level = 6
0.9668
1.1987
1.1725
—
—
—
p-level = 7
1.2283
1.1667
1.1714
—
—
—
p-level = 8
1.2321
1.1742
1.1711
—
—
—
p-level = 9
1.1539
1.1694
1.1711
—
—
—
p-level = 10
1.1458
1.1703
1.1711
—
—
—
4N
F
−0.0020
0.0040
0.0387
0.1132
0.2539
0.3253
R
−0.0006
0.3669
1.2780
1.1706
1.1692
1.1697
S
−0.0714
−0.1049
0.0660
0.7326
1.2604
1.2261
8N
F
0.0937
0.1678
0.3530
0.8797
1.1372
1.1558
R
0.2010
1.0642
1.1594
1.1709
1.1711
1.1711
S
0.1522
0.7955
1.1710
1.1704
1.1710
1.1711
9N
F
0.1022
0.1628
0.3519
0.8797
1.1372
1.1558
R
0.2082
0.9531
1.1692
1.1692
1.1705
1.1707
S
0.2039
0.8190
1.1716
1.1707
1.1710
1.1711Figure 6
Convergence of displacements in the laminated shell.
(a)
w
- at point A
(b)
u
- at point B


## 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 pressurep  (=6410/π) as shown in Figure 2 are estimated. In the problem with two layers, thickness of two layers is identical. Length s of the shell is 20, the radius R is 20, and the total thickness t is 1. Each layer is a unidirectional fiber reinforced composite with the following material constants:(21)E1=7.5×106;E2=2×106;G12=1.25×106;G13=G23=0.625×106;ν12=0.25,where subscript 1 signifies the direction parallel to fibers and 2 and 3 the transverse direction. In these problems, the number 1 means x direction, the number 2 is θ direction, and the number 3 refers to r 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.Figure 2
Orthotropic clamped cylindrical shells with internal pressure.Table1 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 p-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.Table 1
Maximum displacements of the cylindrical shell with one layer.
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
Present analysis
p-level = 3
0.3798
—
—
—
—
p-level = 4
0.3746
—
—
—
—
p-level = 5
0.3748
—
—
—
—
p-level = 6
0.3748
—
—
—
—
p-level = 7
0.3748
—
—
—
—
p-level = 8
0.3748
—
—
—
—
p-level = 9
0.3748
—
—
—
—
p-level = 10
0.3748
—
—
—
—
4N
F
0.1478
0.2816
0.3484
0.3681
0.3731
R
0.1574
0.4072
0.3844
0.3772
0.3754
S
0.1574
0.4072
0.3844
0.3772
0.3754
8N
F
0.3220
0.3729
0.3747
0.3748
0.3748
R
0.4179
0.3770
0.3750
0.3748
0.3748
S
0.4179
0.3770
0.3750
0.3748
0.3748
9N
F
0.3129
0.3728
0.3747
0.3748
0.3748
R
0.4179
0.3770
0.3749
0.3748
0.3748
S
0.4179
0.3770
0.3749
0.3748
0.3748Table 2
Maximum deflections of the cylindrical laminate shell with two layers.
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
Present analysis
p-level = 3
0.1860
—
—
—
—
p-level = 4
0.1782
—
—
—
—
p-level = 5
0.1795
—
—
—
—
p-level = 6
0.1794
—
—
—
—
p-level = 7
0.1794
p-level = 8
0.1794
—
—
—
—
p-level = 9
0.1794
—
—
—
—
p-level = 10
0.1794
4N
F
0.1105
0.1705
0.1811
0.1804
0.1797
R
0.1240
0.2290
0.1882
0.1814
0.1799
S
0.1240
0.2290
0.1882
0.1814
0.1799
8N
F
0.1829
0.1830
0.1797
0.1795
0.1794
R
0.2174
0.1796
0.1794
0.1794
—
S
0.2174
0.1796
0.1794
0.1794
—
9N
F
0.1829
0.1830
0.1797
0.1795
0.1794
R
0.2174
0.1796
0.1794
0.1794
—
S
0.2174
0.1796
0.1794
0.1794
—


## 3.2. Isotropic Cylindrical Shell Panels with Uniformly Distributed Load

Next, an isotropic (E=0.45×106 and ν=0.3) cylindrical shell panel is considered. Figure 3 shows geometry shape of the panel and loading condition (p=0.04) of which specific values are as follows:(22)α=0.2;R=100;s=20;t=0.125.A quarter of panel is considered as a computational region shown in Figure 3 by specifying symmetry conditions with respect to x and θ axes. At x=s, 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 p-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 (S/t=160) and thus the effect of selective reduced integration is worthless.Table 3
Maximum vertical deflection (-w×102) at the center of cylindrical panels.
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
Present analysis
p-level =3
0.4903
1.0048
p-level = 4
0.5410
1.1534
p-level = 5
1.2248
1.1403
p-level = 6
1.2525
1.1358
p-level = 7
1.1175
1.1349
p-level = 8
1.1455
1.1349
p-level = 9
1.1320
1.1349
p-level = 10
1.1359
1.1349
4N
F
0.0082
0.0309
0.1075
0.3378
0.7456
R
0.0111
1.7227
1.2376
1.1577
1.1404
S
0.0110
1.7154
1.2400
1.1584
1.1406
8N
F
0.0461
0.7435
1.1745
1.1428
1.1362
R
1.6496
0.9573
1.1375
1.1350
1.1349
S
1.6399
0.9563
1.1374
1.1350
1.1349
9N
F
0.0463
1.3469
1.1721
1.1427
1.1362
R
1.6668
1.1444
1.1352
1.1349
1.1349
S
1.6532
1.1434
1.1352
1.1349
1.1349Figure 3
Clamped cylindrical shell panels with uniform load.Figure 4
Convergence of maximum displacement with respect to NDF.
(a)
Convergence characteristics by full integration
(b)
Convergence characteristics by reduced and selective reduced integration


## 3.3. Pinched Cylinder with Rigid End Diaphragms

This is another well-known benchmark problem. Figure5 shows a cylinder with rigid end diaphragms subjected to two point loads (P=1). Geometry data in Figure 5 is given by(23)S=600;R=300;t=3.0.For material conditions, isotropic material (E=3×106 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.Figure 5
Pinched cylinder with rigid diaphragms.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. Table4 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 p-refinement is implemented. Also, it is seen in Table 4 that the converged solution of present analysis requires h- and p-refinement simultaneously. For p-level ≥ 7 in the 8×8 mesh, the relative error of maximum displacement is less than 1.0% as compared with the reference value.Table 4
Maximum radial deflection (-w×105) at R/t = 100.
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
20 × 20
Present analysis
p-level = 3
0.0674
0.3540
1.2273
1.7090
1.8172
1.8293
p-level = 4
0.1489
0.8821
1.5689
1.7825
1.8350
1.8404
p-level = 5
0.3503
1.3250
1.6937
1.8147
1.8419
1.8452
p-level = 6
0.7655
1.5009
1.7542
1.8289
1.8458
1.8481
p-level = 7
1.1906
1.6180
1.7937
1.8376
1.8485
1.8503
p-level = 8
1.3670
1.7037
1.8162
1.8426
1.8504
1.8520
p-level = 9
1.4789
1.7560
1.8284
1.8457
1.8519
1.8533
p-level = 10
1.5897
1.7883
1.8356
1.8479
1.8532
1.8545
4N
F
0.0098
0.0254
0.0608
0.1282
0.2785
0.3603
R
0.0207
0.3263
1.9281
1.8453
1.8600
1.8632
S
0.0213
0.6723
1.3929
1.6093
1.7762
1.8056
8N
F
0.0369
0.0864
0.3891
1.1772
1.6691
1.7378
R
0.1477
1.0976
1.6444
1.7951
1.8307
1.8365
S
0.1263
0.8565
1.5541
1.7799
1.8282
1.8350
9N
F
0.0414
0.0925
0.4184
1.2238
1.6845
1.7476
R
0.1886
1.7458
1.8451
1.8596
1.8677
1.8702
S
0.1657
1.2428
1.7046
1.8359
1.8633
1.8675Next, the cylindrical laminated shell with two layers is analyzed. The thickness of two layers is identical and total thicknesst of the cylinder is fixed as 3.0. Each layer consists of a unidirectional fiber reinforced composite with following material constants:(24)EL=25ET;GLT=0.5ET;GTT=0.2ET;νLT=νTT=0.25,where subscript L signifies the direction parallel to the fibers, T is the transverse direction, and νLT is the Poisson’s ratio for strain in the T direction under uniaxial normal stress in the L direction. ET is an elastic modulus in the transverse direction and G is an elastic shear modulus. The normalized quantities of interest are defined as follows:(25)w-=-10wAE1t3PR2;u-=-10uBE1t2PR,where wA means the deflection at the point A and uB refers to the displacement in the x direction at the point B. Figures 6(a) and 6(b) show the convergence characteristics of w and u 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 hp-refinement, while the values of 8NR/S and 9NR/S reach 1.24 and 1.31. The difference between present solutions and h-refinement analyses is within 2~3%. Table 6 shows displacements in x 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.Table 5
Maximum radial deflection at the point A (w-).
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16 × 16
20 × 20
Present analysis
p-level = 3
0.0895
0.2352
0.7834
1.1560
1.2305
1.2393
p-level = 4
0.1366
0.5876
1.0602
1.2086
1.2438
1.2491
p-level = 5
0.3234
0.8951
1.1517
1.2297
1.2518
1.2559
p-level = 6
0.6028
1.0129
1.1894
1.2408
1.2573
1.2610
p-level = 7
0.8076
1.0906
1.2156
1.2486
1.2618
1.2651
p-level = 8
0.8935
1.1485
1.2310
1.2538
1.2654
1.2687
p-level = 9
0.9759
1.1854
1.2404
1.2579
1.2686
1.2718
p-level = 10
1.0577
1.2085
1.2467
1.2611
1.2714
1.2745
4N
F
0.0033
0.0169
0.0539
0.1391
0.3169
0.4038
R
0.0349
0.2760
1.2670
1.2693
1.2915
1.2978
S
0.0351
−0.0530
−0.6145
5.6803
1.4209
1.3725
8N
F
0.0581
0.1218
0.2915
0.7831
1.1546
1.1990
R
0.1446
0.8093
1.1058
1.2021
1.2359
1.2421
S
0.1228
0.7113
1.0710
1.1967
1.2349
1.2414
9N
F
0.0630
0.1314
0.3095
0.8082
1.1669
1.2086
R
0.1705
1.1902
1.2647
1.2878
1.3069
1.3130
S
0.1446
0.7906
1.2049
1.2767
1.3041
1.3112Table 6
Displacement inx direction at the point B (u-).
Analysis types
1 × 1
2 × 2
4 × 4
8 × 8
16  × 16
20 × 20
Present analysis
p-level = 3
0.1627
0.3403
0.9709
—
—
—
p-level = 4
0.2256
0.9344
1.1661
—
—
—
p-level = 5
0.5338
1.2390
1.1712
—
—
—
p-level = 6
0.9668
1.1987
1.1725
—
—
—
p-level = 7
1.2283
1.1667
1.1714
—
—
—
p-level = 8
1.2321
1.1742
1.1711
—
—
—
p-level = 9
1.1539
1.1694
1.1711
—
—
—
p-level = 10
1.1458
1.1703
1.1711
—
—
—
4N
F
−0.0020
0.0040
0.0387
0.1132
0.2539
0.3253
R
−0.0006
0.3669
1.2780
1.1706
1.1692
1.1697
S
−0.0714
−0.1049
0.0660
0.7326
1.2604
1.2261
8N
F
0.0937
0.1678
0.3530
0.8797
1.1372
1.1558
R
0.2010
1.0642
1.1594
1.1709
1.1711
1.1711
S
0.1522
0.7955
1.1710
1.1704
1.1710
1.1711
9N
F
0.1022
0.1628
0.3519
0.8797
1.1372
1.1558
R
0.2082
0.9531
1.1692
1.1692
1.1705
1.1707
S
0.2039
0.8190
1.1716
1.1707
1.1710
1.1711Figure 6
Convergence of displacements in the laminated shell.
(a)
w
- at point A
(b)
u
- at point B


## 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 overp-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.



---

*Source: 676181-2015-02-26.xml*