Abstract
This paper deals with the hierarchical laminated shell elements with nonsensitivity to adverse conditions for linear static analysis of cylindrical problems. Displacement approximation of the elements is established by high-order shape functions using the integrals of Legendre polynomials to ensure continuity at the interface between adjacent elements. For exact linear mapping of cylindrical shell problems, cylindrical coordinate is adopted. To find global response of laminated composite shells, equivalent single-layer theory is also considered. Thus, the proposed elements are formulated by the dimensional reduction from three-dimensional solid to two-dimensional plane which allows the first-order shear deformation and considers anisotropy due to fiber orientation. The sensitivity tests are implemented to show robustness of the present elements with respect to severe element distortions, very high aspect ratios of elements, and very large radius-to-thickness ratios of shells. In addition, this element has investigated whether material conditions such as isotropic and orthotropic properties may affect the accuracy as the element distortion ratio is increased. The robustness of present element has been compared with that of several shell elements available in ANSYS program.
1. Introduction
Finite element methods generally involve finding approximate solutions in a space of piecewise polynomials of degree, which is often designated by the letter , on a grid of mesh size . In general, conventional finite elements based on mesh refinement have given reliable solutions when discretizing mesh is regular, while poor performance would be obtained when the element geometry is distorted. Also, it is known well that the interpolation precision of quadrilateral finite elements deteriorates if the element geometry considerably differentiated itself with a square. In this sense loss of element performance is commonly associated with a gradual increase of the stiffness, leading to sort of locking the element’s response. It brings about the inability of the element to find a good approximation of the solution. Moreover things would turn for the worse when conventional low-order finite element formulation is used in thin plates or shells since shear or membrane locking phenomena arise. Various robust schemes have been suggested for the problems involving locking. One possibility is to use higher-order elements. It is known that the locking completely eliminates certain types of locking [1]. The high-order finite element implementations have also been advocated in recent years as a means of eliminating the locking phenomena completely. The relevant works have been implemented by some researchers [2–5]. The issue of locking using lower-order elements has been most prominently addressed through the use of low-order finite technology using some mixed variational principles [6–8]. However, these depend upon reformulating the problems in special ways which have not been required in higher-order elements. In this paper, we will address only the finite element formulation for laminated shell behavior using the -version approach. The first successful -version formulation related to shells was reported by Woo and Basu [9] who presented the hierarchical -shell element formulation in the cylindrical coordinates associated with a suitable transfinite mapping function to represent the curved geometry. This element showed not only a constant strain and a rigid body motion from patch and eigenvalue tests but also a strong robustness with regard to very large aspect ratio and severely distorted mesh. The proposed hierarchical -shell element can also be successfully extended to singularity problems including a crack and a cut-out. Surana and Sorem [10] developed a three-dimensional curved shell element with a higher-order hierarchical displacement approximation in the thickness direction of the shell. The extension of this idea to laminated plates and shells was accomplished by the same authors [11]. Laminated plates and shells have usually been analyzed by the use of ESL (equivalent single-layer) theories based on either the classical Kirchhoff-Love hypothesis or first-order shear deformation theories. A drawback of ESL theories in modeling composite laminated plates and shells is that the transversal strain components are continuous across interfaces between dissimilar materials for the perfectly bonded layers. Therefore, the transverse stress components are discontinuous at the layer interfaces. On the other hand, layerwise theories are considerably more accurate than the preceding theories. Thus, Ahn and Woo [12, 13] extended the -version approach to efficient “mixed model analysis,” often called “global-local approach” to obtain interlaminar stresses at free edges in the laminated composite plate under extension and flexure on the basis of layerwise theories. The aim of this study is to present a simple high-order hierarchical -element for laminated composite shells using ESL concepts prior to layerwise theories in the cylindrical coordinates. This approach is computationally less expensive as compared to those obtained by three-dimensional elasticity solutions and layerwise finite element models.
In general, the most important symptoms of accuracy failure in modern finite elements are spurious mechanism, locking, elementary defects like violation of rigid body property, and invariance to node numbering. Also, parameters which affect accuracy are loading, element geometry, problem geometry, material properties, material anisotropy, and so on. Because of these reasons, governmental concern in the USA for the accuracy of finite element analysis is evidenced by NRC (Nuclear Regulatory Commission) requirement for structural analysis computer program validation, and the formation of NAFEM (National Agency for Finite Element Methods and Standards) in the United Kingdom.
In this study, the sensitivity test has been carried out to verify the robustness of proposed element in relation to distortion effect of mesh, very high aspect ratio, and very large radius-to-thickness ratio. In addition to these, the orthotropic cylindrical shell stacking with different fiber orientations is tested to check whether the anisotropy of materials may affect the accuracy as the element distortion ratio is increased. The numerical solutions obtained by present element are compared with several shell elements available in ANSYS program [14].
2. Formulation of Hierarchical Laminated Shell Element
2.1. Hierarchical Shape Functions for Quadrilaterals
A standard quadrilateral element is shown in Figure 1. Two-dimensional shape functions in the element can be classified into three groups such as nodal, edge, and internal shape functions. The nodal shape functions can be defined as and denote the local coordinate of the th node in the Figure 1. Also, the edge shape functions with any order , which vanish in all other edges, are defined separately for each individual edge. ConsiderThe quantity refers to the Kronecker tensor. and denote the local coordinate of both nodes in the th edge shown in the Figure 1. The functions which have one independent variable of the two variables and are defined as follows:where the functions are the Legendre polynomials. These functions with integrals of Legendre polynomials are well suited for computer implementation and have very favorable properties from the point of view of numerical stability [15]. Lastly, there are internal shape functions in . These can be written as
2.2. Geometry Fields
One element with a curved surface is on a cylindrical shell shown in Figure 2(a). Linear mapping on a general Cartesian coordinate is not accurate enough for the element. In the case of degenerated shell concept based on -FEM which is normally adopted in commercial codes, a nodal coordinate system defined at each nodal point with any reference surface is considered [16]. However, the mapping technique explains the curved surface approximately, not exactly, in which the size of the elements is made small enough to follow the curved surface as closely as desired.
In -FEM, this option is not acceptable and, to avoid geometrical sources of error, exact mapping of a curved surface is necessary because one of the sources of efficiency of -FEM is the use of as few elements as possible. Woo and Basu [9] proposed the transfinite mapping with trigonometric polynomial interpolation projects based on blending function method of -FEM. The mapping technique ought to be considered for circular or oval hole to highlight modeling simplicity of -FEM. However, as surfaces on cylindrical shells are only considered in the present study, the simple linear mapping on cylindrical coordinate shown in Figure 2(a) is enough. Figures 2(a) and 2(b) present the mapping concept of the computational domain with the shape functions of four nodes in (1) aforementioned. For the mapping process of the curved surfaces on cylindrical shells with constant curvature radius along an axial direction, the equation can be written as
2.3. Displacement Fields
For modeling of composite structures, ESL elements (Figure 3) derived from function theories [17] in a laminated system are obtained by dimensional reduction from three-dimensional solid to two-dimensional representative surface, satisfying plane stress theory for membrane behavior and first-order shear deformation theory for bending behavior.
Thus, in the case of a quadrilateral subparametric element, there are three translational displacement components and two rotational components corresponding to each vertex, side, and internal mode. For ESL model of an element with any -level, the two in-plane displacement components can be defined asHere, the shape function functions are defined in (1), (2), and (4). Also, and are the nodal variables for translation and rotation, corresponding to the four nodes. Likewise, the values and denote translational and rotational nodal variables corresponding to side and internal modes, respectively. The variable denotes the distance from the reference surface in the thickness direction of a point of interest on standard coordinate. The transverse displacement field of a point in the ESL element with any -level is defined asIn this equation, refers to the nodal variable for the four nodes and and correspond to the nodal variables for side and internal modes, respectively. Equation (7) connotes that the lateral displacement is independent of thickness coordinate, which is contrary to in-plane displacement components. This choice signifies that the transverse normal stress is zero. The constitutive relationship for the ESL element with respect to the local element coordinate system in reference surface can be defined aswhere refers to the full elasticity matrix of a -layered ESL element and is composed of submatrices as shown below. ConsiderHere, is the extensional elasticity matrix, is the bending elasticity matrix, is the coupling between bending and extensional elasticity matrices, is the transverse-shear elasticity matrix, and denotes a null matrix. The summation accounts for the contributions from all the layers. For a typical layer , the submatrices are calculated using its elasticity matrix allowing for anisotropy with three mutually orthogonal planes of symmetry. The elasticity matrix with reference to the elemental coordinate system is obtained from the elasticity matrix referred to the material axes by using the coordinate transformation matrix . The resulting transformation relationship takes the following triple product standard form:In (10), 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. In the case of heterogeneous plates, this factor is based on equilibrium equations and strain energy components used in [18]. The individual strain components with respect to , , and axes shown in (8) can be identified aswhere superscripts , , and refer to membrane, bending, and transverse strain components, respectively. The total strain vector for an arbitrary point will have five components as
For a point at a distance from the reference surface, the flexural components of planar strains can be separated as follows:Based on the strain vector in (8), the general strain vector for an ESL element can be expressed aswhere is the strain matrix composed of derivatives of hierarchical shape functions and is the vector of variables of the ESL element with any number of degrees of freedom related to order. The corresponding stress resultants in any layer can now be expressed aswhere , , and are membrane force, bending moment, and transversal shear resultants, respectively. and are the -axes of bottom and top surfaces of a typical layer with respect to the reference surface.
3. Sensitivity to Geometric Parameters
In this study, the sensitivity test has been carried out to verify the robustness of present element with respect to three critical conditions including severe element distortion, very high aspect ratio, and very large thickness ratio. The cylindrical shell with length, radius, and thickness dimensions specified by , , and , respectively, is subjected to a pinch load as shown in Figure 4. Due to the symmetry condition, one octant of the shell is modeled by 2 × 2 mesh design. The sensitivity to input parameters by the present element with different -levels has been compared with several shell elements available in ANSYS program which are two -convergent elements of SHELL181 (4-node) and SHELL281 (8-node) and one -convergent element of SHELL150. In the case of SHELL150 element, the shape functions are formulated by Lagrangian polynomials and curvilinear mapping technique is adopted in Cartesian coordinate. On the other hand, linear exact mapping in cylindrical coordinate is used in the present -FEM. For three cases of sensitivity tests, two material properties are considered in (16) and (17) as follows: Case 1. Isotropy is as follows: Case 2. Orthotropy is as follows:Also, the numerical results are normalized using the classical solution or the estimated exact solution. Thus, the normalized central deflection is used in this study which is defined inHere, is the maximum deflection in the center of shells that is obtained by numerical analysis, and represents the estimated exact solution.
(a) Pinched cylinder
(b) Finite element modeling in computational domain
3.1. Effect of Element Distortion
In this section, the sensitivity of the present element has been tested with respect to severe element distortion for both isotropic and orthotropic cases. The material conditions are assumed to be isotropic if there is no mention about material property. To check the element accuracy for element distortion effects, we employ the element distortion ratios denoted by to the central node of the 2 × 2 mesh when is fixed as 53. Thus the element distortion ratios are defined by (19). In the cases of SHELL181 and SHELL281, two mesh types (2 × 2 and 4 × 4) are investigated to model one octant of the shell exploiting symmetry:
The convergence characteristics of the present elements with respect to different -levels are plotted in Figure 5 when and . From this figure, it is seen that the convergence begins from by using the hierarchical shell elements as well as the SHELL150 elements based on -FEM. Those results are compared with the analytical solution by Takemoto and Cook [19]. The present solution virtually converges to the normalized exact solution denoted by 1.0. However, the numerical solution by SHELL150 elements is converged to upper bound of 1.2. This tendency is mainly due to the use of different shape functions and mapping technique for curved boundary. Linear mapping technique on rectangular coordinate is considered as interpolation functions of geometry fields for SHELL150 elements. Therefore, SHELL150 elements require mesh refinement for curved shapes to improve accuracy although they are -FEM. In Table 1, the numerical results obtained by the present elements are tabulated for with different element distortion ratios and compared with the references. It is noted that the hierarchical shell elements based on 2 × 2 mesh with show an excellent behavior even for extremely distorted mesh of . All results are very close to the normalized displacement 1.0 which represents the approximate solutions equal to exact solution. However, other solutions by SHELL181 and SHELL281 show poor accuracy and behave relatively stiff and converge very slowly even though the meshes are refined up to 16 elements (4 × 4 mesh). Thus, it is concluded that numerical solutions in references by ANSYS program degrade significantly with increasing element distortion. These remarks can be reconfirmed by the graphical illustration as shown in Figure 6. From this figure, it is advised that the element distortion ratio should not exceed roughly 10% for 4-node and 20% for 8-node -convergent shell elements to get good displacement results. On the contrary, the effect of element distortion on the numerical results by the present element is plotted in Figure 7 as the element distortion ratios vary from 0.0% to 90% on 2 × 2 mesh design. It is observed that the present element shows strong robustness to severe element distortion even for extreme case of , especially for .
Next example is another pinched cylindrical problem composed of orthotropic materials defined in (17). The orthotropic pinched cylindrical shell stacking with different fiber orientations is tested in order to check whether material conditions may affect the accuracy as the element distortion ratio is increased. For this purpose, the laminated cylindrical shell with a pinch load is modeled by two layers with cross-ply (0°/90°). Similar to the isotropic problem previously, the present elements may endure severe element distortion up to from . The details of numerical solutions are presented in Table 2 as the -level increases from 1 to 10. It is noticed that the present shell elements based on -FEM exhibit strong robustness with respect to severe element distortion regardless of the anisotropy of materials.
3.2. Effect of Aspect Ratios of Elements
In general, Robinson [20] and Macneal and Harder [21] suggested a single element test in which response is examined as the element aspect ratio is changed. It is noticed that element aspect ratios should not exceed roughly seven for good displacement results and roughly three for good stress results. Of course, this conclusion can be changed according to types of problems. However, since large aspect ratio is one of the sources causing numerical errors, it is desirable to make the aspect ratio one, especially in the vicinity of a singular point. The aspect ratio is defined by (20), and and are denoted in Figure 8. ConsiderThe finite element results obtained by present elements with are presented in Table 3 with respect to aspect ratio and compared with reference values by SHELL181, SHELL281, and SHELL150 . All results are based on 2 × 2 mesh design. When the aspect ratio is 1.0, the normalized transverse deflections at the center by SHELL281, SHELL 150, and present elements are close to an exact solution. However, the results by SHELL181 elements fail to converge. Precisely speaking, the relative errors of present elements to an exact solution are below 3% until . On the other hand, the relative errors of SHELL281 and SHELL150 show approximately more than 10% from the aspect ratio . The graphical solution is also presented in Figure 9 to easily understand the deviation of numerical errors in reference to the exact solution. The tolerance of present elements is represented in Figure 10 as the aspect ratios are increased from 1 to 1000. It is seen that the hierarchical shell element tolerates the large aspect ratios up to 1000 under 3% accuracy of relative error.
3.3. Effect of Thickness Ratios
In order to investigate the membrane and shear locking phenomena of hierarchical shell elements, the central deflections have been calculated under very large radius-to-thickness ratios . Thus, the tolerance of numerical solutions has been exhibited in Table 4 in regard to ratios ranging from 10 (thick) to 1000 (extremely thin). Present solutions show good agreement with those by 8-node shell element denoted by SHELL281 available in ANSYS program. It is noted that both elements avoid the shear locking regardless of extremely thin case.
4. Conclusions
In the mesh design of the finite element analysis of shells, the sources of numerical errors such as severe element distortions, very high aspect ratios, and very large radius-to-thickness ratios cause the numerical instability of a stiffness matrix. The proposed hierarchical shell elements based on -FEM have been tested from the point of view of the element robustness. The obtained results and subsequent future works can be summarized as follows.(1)It is observed that the present elements show strong robustness to severe element distortion even for extreme case of , especially for .(2)It is noticed that the hierarchical shell elements based on -FEM exhibit a strong robustness with respect to severe element distortion regardless of cross-ply lamination scheme of materials.(3)It is seen that the proposed elements tolerate the large aspect ratio up to 1000 under 3% accuracy of relative error.(4)It is noted that the present elements avoid the shear and membrane locking regardless of the thin cases extremely denoted by .(5)In this paper, performance of the present elements is only limited to cylindrical shells based on linear analysis. From these results, in future it is necessary to investigate geometrical nonlinearity, and critical shell shapes such as spherical, hyperbolic shells, which would differ with cylindrical shells.
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 the 2014 Yeungnam University Research Grant.