Evaluation of shear and membrane locking in refined hierarchical shell finite elements for laminated structures

Shear and membrane locking phenomena are fundamental issues of shell finite element models. A family of refined shell elements for laminated structures has been developed in the framework of Carrera Unified Formulation, including hierarchical elements based on higher-order Legendre polynomial expansions. These hierarchical elements were reported to be relatively less prone to locking phenomena, yet an exhaustive evaluation of them regarding the mitigation of shear and membrane locking on laminated shells is still essential. In the present article, numerically efficient integration schemes for hierarchical elements, including also reduced and selective integration procedures, are discussed and evaluated through single-element p-version finite element models. Both shear and membrane locking are assessed quantitatively through the estimation of strain energy components. The numerical results show that the fully integrated hierarchical shell elements can overcome the shear and membrane locking effectively when a sufficiently high polynomial degree is reached. Reduced and selective integration schemes can help with the mitigation of locking on lower-order hierarchical shell elements.


Introduction
Shell structures, especially composite laminated shells, have been widely used in modern engineering due to their high efficiency in holding loads. A variety of shell theories have been proposed, including the Classical Lamination Theory (CLT) based on Kirchhoff-Love hypothesis [1][2][3], the First-Order Deformation Theory (FSDT) based on the Mindlin-Reissner assumption [4,5], a series of Higher-Order Theories (HOT) [6][7][8], and a variety of refined shell theories generated in the framework of Carrera Unified Formulation (CUF) [9]. Refined shell theories have been implemented in the Finite Element (FE) method and can provide solutions with great accuracy [10,11]. Towards implementing numerically efficient simulation methods for multi-layered structures, hierarchical functions have been adopted in the construction of refined shell elements in the framework of CUF [12]. Nevertheless, the effectiveness of hierarchical functions concerning the cal elements using Naghdi shell model on isotropic structures was reported by Chinosi et al. [38]. The selective and reduced integration schemes on hierarchical elements were primarily discussed by Della Croce and Scapolla [39]. Carrera et al. [12] reported a comparison between shell elements with MITC and plain hierarchical shape functions in the analysis of laminated structures.
The numerical efficiency of the hierarchical elements has been reported by many researchers [34,35,40,41]. In the framework of CUF, this type of hierarchical functions has been used on refined beam [42], plate [43], and shell [12] finite element models for multi-layered structures. Via refined hierarchical 2D elements, the FE models can be mathematically enriched on both the kinematic and shape function levels, leading to an adaptable refinement FE approach with optimal numerical efficiency. This article aims to present the evaluation of hierarchical elements concerning the mitigation of shear and membrane locking phenomena in the analysis of multi-layered shells. In the following sections, we first introduce the multi-layered shell models and CUF briefly. Then, an energy decomposition method is presented for refined shell models based on CUF. Thirdly, numerically efficient full, reduced, and selective integration schemes are discussed. Finally, the mitigation of shear locking and membrane locking is demonstrate through two numerical examples, respectively. The numerical evaluations are compared against available analytical solutions in the literature.

Refined shell finite element formulation
Preliminaries of multi-layered shell model Figure 1 shows a typical laminated shell structure with curvatures. This geometry can be described in the orthogonal curvilinear reference system (α, β, z), in which α and β indicate the two in-plane directions and z the through-thickness direction which is usually measured with reference to the middle surface. The infinitesimal in-plane area dS and volume dV can be written as: where d is the infinitesimal in-plane area on the middle surface, and H α , H β and H z are: in which R α and R β are the principal radii of the middle surface, A and B the coefficients of the first fundamental form of . The present work considers only shells with constant curvatures, for which A = B = 1. For more details about shell formulations, the reader is referred to [44,45]. The strains and stresses defined in the curvilinear reference system are: The strains can be obtained through the geometrical relations: = bu (6) wherein u = {u, v, w} T is the displacement vector, and b the differential operators matrix, whose explicit expression reads: The stresses can be attained from the constitutive equations as follows: (8) in which C is the material coefficients matrix which is obtained by transforming the original form C 0 from the material coordinate system (1, 2, 3) to the global system (α, β, z). The orthotropic material coefficients are characterized by nine independent coefficients, namely Young's moduli, shear moduli, and Poisson ratios [7].

Carrera Unified Formulation (CUF) for refined 2D models
Through Carrera Unified Formulation (CUF), the displacement field of a shell structure can be assumed as: where u τ (α, β) are the in-plane displacement vectors, and functions F τ (z) are related to the theories of structures (TOS). Since F τ (z) depends only on the thickness coordinates, they are also referred to as thickness functions. By increasing the polynomial order of these thickness functions, the shell kinematic models can be refined. Both Equivalent-Single-Layer (ESL) and Layer-Wise (LW) models can be accounted for in this framework, as elaborated by Carrera et al. [46]. The FSDT [47] can be treated as a particular case of the HOT models adopting Taylor expansions (TE). In the LW model framework, Lagrangian polynomial expansions (LE) can be used to formulate kinematics with only translational degrees of freedom. More discussions about these two types of refined kinematic models can be found in the work of Carrera et al. [46]. By using LE, the interfacial continuity of transverse stresses can be approximately achieved when the thickness functions are adequately refined, as demonstrated by Carrera et al. [26]. When the FE discretization is introduced, the in-plane displacements of a shell structure are approximated through the shape functions N i (α, β) through: (10) in which u iτ are nodal unknown variables. The above expression leads to FE formulation in the framework of CUF: where δ indicates the virtual variation. The above expression is compact through the use of repeated indexes. By applying the Principle of Virtual Displacements (PVD), the general expressions of the stiffness matrix and load vector of the FE model, namely the Fundamental Neuclei (FNs), can be obtained. The explicit expressions of the FNs are given in [12]. As expounded by Carrera et al. [46], the CUF-type FE formulation is independent of the kinematic assumptions adopted and is a general framework for the development of refined FE models. For more details about the derivation of shell FE formulations in the framework of CUF, the reader is referred to the work of Carrera et al. [46].

Decomposition of strain energy in refined shell models
For a general laminated shell structure, the strain energy can be decomposed into four parts as follows: where E pn represents the in-plane normal energy, E ps the in-plane shear energy, E zs the transverse shear energy, and E zz the thickness stretch energy. The transverse shear energy allows us to evaluate the shear locking effects in shell elements. The introduction of the thickness stretch energy makes it convenient to assess the performance of the adopted structural theory. Note that the above decomposition applies to arbitrarily laminated shell structures.
To calculate the strain energy components, their corresponding stiffness matrices are necessary. These matrices can be obtained through standard FE procedure in the framework of CUF. Taking the transverse shear energy E zs as an example, by recalling the PVD, one has: δE zs = V (δε αz σ αz + δε βz σ βz )dV = δu js · k zs ijτ s · u iτ (16) wherein k zs ijτ s is the FNs for the transverse shear stiffness matrix of the element. By substituting the displacement approximations (Eq. 11) into the geometrical relations (Eq. 6), and considering the constitutive equations (Eq. 8), one obtains: thus the FNs for the transverse shear stiffness matrix k zs itτ s can be obtained as: Through the assembly of k zs ijτ s according to the standard procedure of FE formulation in CUF framework, the transverse shear stiffness matrix K zs can be obtained. When the displacement solutions are obtained, the transverse shear strain energy can be calculated through: The in-plane normal stiffness matrix K pn , in-plane shear stiffness matrix K ps , and outof-plane normal stiffness matrix K zz can be achieved accordingly by means of the following FNs: wherein b pn , b ps , and b zz are the sub-matrices of the differential operators matrix b as in Eq. 7, and their explicit expressions are: and their corresponding material coefficients matrices (sub-matrices of the material coefficients matrix C) are as follows: C ps = C 61 C 62 C 63 C 64 C 65 C 66 (28) In the end, the complete stiffness FNs can be obtained as the summation of these terms as: If the multi-layered shell has symmetric lamination properties, the neutral surface of bending will coincide with the geometrical middle surface, and the in-plane normal strain energy, as in Eq. 12, can be further decomposed into membrane energy E memb and bending energy E bend conveniently through: wherein ε 0 αα and ε 0 ββ are the normal strains due to the mid-surface straining, and they can be attained by means of: By following the procedure described before, k memb ijτ s , the FNs for the membrane stiffness matrix K memb , can be derived. The bending energy can be then obtained through: This separation of membrane and bending energy components provides the convenience to evaluate the existence of membrane locking and better understand the structural responses. It should be noted that E pn and E ps are both in-plane strain energy components, however since in laminated plates and shells the calculation E ps does not dependent on a specific neutral surface as the membrane and bending energy components do, it is considered apart in the present article.

Integration schemes for hierarchical elements
This section addresses the efficient integration schemes of 2D hierarchical elements with full, reduced, and selectively reduced integration. According to Szabó et al. [35,41], the hierarchical shape functions for 2D elements can be classified into nodal modes, edge modes, and surface modes, which can be expressed in a unified manner as: in which the basis polynomials φ m (ξ ) and φ n (η) are decided by their corresponding mode and the polynomial degree, as it is illustrated in Fig. 2.
In the FNs of stiffness matrix as given by Li et al. [12], the contribution of the shape functions to the stiffness matrix accounts for the following integrals: where · · · represents · · · dξ dη. Among these terms, for given i and j combination, the polynomials with the highest order is N i · N j . Consider the product of N i and N j : Thus in the ξ and η directions, the highest polynomial orders are m + r and n + s, respectively. For simplicity, in the present work, the same set of Gauss points are used to calculate the above integrals of given N i and N j .

Full integration scheme
In Gauss-Legendre quadrature, n Gauss points guarantee the exact integration of a polynomial of order 2n − 1. For the exact integration of N i N j , the least number of Gauss points used in the ξ direction, N GX , should be: wherein N is an arbitrary positive integer. The above expression also applies to the calculation of number of Gauss points in the η direction, N GY , for an exact integration. For classical Lagrangian elements, since all the shape functions have the same polynomial order in both ξ and η directions, the scheme of Gauss points can be uniformly decided. Differently, for hierarchical elements, the required number of Gauss points varies according to different combinations of shape functions. Also, in practice, for given N i and N j , the same set of Gauss integration points can be used for the nine integrals in Eq. 36, which is determined by the highest polynomial order given by N i · N j . Figure 3 presents two examples of the Gauss points distribution when full integration is used on hierarchical elements. Figure 3a shows the numerical calculation of N 13 · N 14 needs 3 × 3 Gauss points. N 13 · N 15 requires 5 × 2 Gauss points for its full integration, as illustrated in Fig. 3b.

Reduced integration scheme
The reduced integration technique on the hierarchical elements can be used in the following manners: applying the reduced integration to polynomials with the highest order among all N i · N j combinations, and using full integration for all the lower-order polyno- where the highest order of the polynomials to be integrated is 2p, and p Gauss points are needed for the reduced integration. The product polynomials to be integrated in the other direction are of the second order and are fully integrated by using two Gauss points. Meanwhile, all the lower-order terms should be exactly integrated. The hierarchical element with p = 4 can be taken as an example. The polynomials to be integrated with the highest order in the ξ direction are N 13 · N 13 , N 13 · N 15 , N 15 · N 13 , and N 15 · N 15 . Those with the highest order in the η direction are N 14 · N 14 , N 14 · N 16 , N 16 · N 14 , and N 16 · N 16 . When the reduced integration scheme is adopted, 4 × 2 Gauss points should be used for the first group of polynomials (see Fig. 4a), and a 2 × 4 mesh of Gauss points for the second set (see Fig. 4b). For this fourth-order hierarchical element, the integration schemes that should be used for different blocks have been indicated in Fig. 5. Note that each block represented by a square is a sub-matrix K ij .

Selectively reduced integration scheme
In the selectively reduced integration method, the low-order integration is only applied to those terms related to the transverse shear stiffness. This technique is aimed to reduce the transverse shear stiffness to alleviate the shear locking phenomenon. These components can be determined by considering the FNs of the transverse shear stiffness matrix in Eq. 19. Also, the sub-matrices of the stiffness matrix K ij that should be selectively integrated follow the same rule as the reduced integration technique as discussed in the "Reduced integration scheme" section.    Figure 6 reports the eigenvalues of the stiffness matrices of hierarchical elements (p = 1, 2, 3) with reduced and selective integration. Note that hierarchical elements with p = 1 are equivalent to standard Q4 (four-node quadrilateral Lagrangian) element. Figure 6a shows that when p = 1, both reduced and selective integration schemes lead to more than six zero eigenvalues numerically, which means the elements are not robust. From Fig. 6b, it can be observed that for polynomial degree p = 2, the reduced integration leads to two spurious modes (which is equal to the number of thickness functions used), and the element with selective integration has exactly six rigid-body modes. When the polynomial degree is further increased to p = 3, the spurious modes are eliminated on the elements with reduced integration, see Fig. 6c. The results demonstrate that, it can be guaranteed that there is zero spurious mode for the reduced integrated hierarchical elements when p ≥ 3, and no spurious mode exists for selective integration when p ≥ 2.

Results and discussion
This section presents two numerical examples on laminated shells considering a wide range of aspect ratio. Single-element FE models are used with p-version refinements. It is obvious that an element with only linear shape functions is not adequate for the modeling, thus the refinement of the shape functions starts from p = 2. The polynomial degree is increased till the chosen convergence threshold is reached. Two kinds of TOS are used in the numerical modeling, namely the FSDT and LE4 (LW model with fourth-order Lagrange polynomials in each layer, see Carrera et al. [46]). These two theories are compared through elements with full integration. Then, with LE4 theory, different integration techniques on the hierarchical elements are compared, including full integration (FULL), reduced integration (REDI), and selective integration (SELI). Besides the displacement and stress evaluations, the strain energy components are also reported. The numerical results are compared against available analytical reference solutions.

Three-layered cylindrical shells under distributed pressure
Three-layered cylindrical shells with symmetric lamination (0 • /90 • /0 • ) are studied. The analytical solutions were presented by Varadan and Bhaskar [48]. The three layers have equal thickness h/3. Radius-to-thickness ratios R β /h ranging from 2 to 500 are investigated. The cylindrical shells are simply supported on the two ends and subjected to transverse distributed pressure on the inner surface. The load distribution follows: Fig. 7 The three-layered simply supported cylindrical shells under inner distributed pressure where L is the cylinder length and R β the middle surface radius, and L = 4R β . Figure 7a directions of the fibers, respectively. By taking advantage of the symmetric features in the cylinder axial direction and the cyclic conditions in the circumferential direction, 1/16 of the structure is taken to build the FE model, as indicated by the shaded zone in Fig. 7. The displacement and stress results are non-dimensionalized through: The 1/16 FE models contain only one element, and the numerical accuracy is improved by increasing the polynomial degree gradually when the relative difference compared to the one-order-lower case is less than 1% regarding the deflection and stresses as well as the energy components. Considering the load distribution, this benchmark is quite challenging for a single-element model. Table 1 summarizes the converged solutions for each radius-to-thickness ratio value. In general, as the shell structure gets thinner, a higher polynomial degree is required to achieve the desired accuracy. The displacement and stress evaluations obtained with LE4 kinematics agree well with the reference solutions given by Varadan and Bhaskar [48]. The accuracy of the FE results of the thick shell (R β /h = 2) can be further improved by refining the thickness functions (TOS), as reported by Li et al. [12]. FSDT leads to good estimation of displacement and in-plane stresses for the thinner shells (R β /h = 50, 100, 500), but fails in other stresses. Also, unlike the LE4 kinematic model,  Fig. 8 Convergence regarding the normalized deflection of FE models for the three-layered cylindrical shells under distributed pressure, for various radius-to-thickness ratio R β /h FSDT ignores the stretch effects in the thickness direction which may play an essential role in thick shells, such as the R β /h = 2 and R β /h = 4 cases in this numerical example.

Table 1 Displacement and stress evaluation on the three-layered cylindrical shells with various
The convergence of the normalized deflections for each aspect ratio is shown in Fig. 8, in whichw =w/w ref andw ref is the reference deflection solution. Figure 9 reports the convergence of the FE model regarding the strain energy error, which is calculated by taking the converged solution employing LE4 kinematics with full integration as the reference. It can be observed that for shells with different R β /h values, convergence is achieved at different polynomial degrees. Generally, the thinner the shell structure is, the higher polynomial order will be required. Compared to the full integration scheme, the . 9 Convergence regarding the strain energy of FE models for the three-layered cylindrical shells under distributed pressure, for various radius-to-thickness ratio R β /h reduced integration technique can help to increase the accuracy in the low-order cases, but the eventual convergence is reached at the same time with full integration. Note that for this single-element model, the detected spurious modes in reduced integrated elements with p = 2 are not observed to be a significant problem. Notably, the selectively reduced integration improves the accuracy of the single-element FE model with p = 2 and leads to results quite close to those obtained with full integration in the higher-order cases (p ≥ 3). When the numerical convergence is reached, all the three kinds of integration schemes lead to results that agree well with the reference solutions.  Fig. 11 Energy components versus radius-to-thickness ratio R β /h on the three-layered cylindrical shells under distributed pressure Figure 10 shows the variation of the ratio of strain energy components with the increase of the polynomial degree of the hierarchical element. Regarding Fig. 10: • It can be found that for the transverse shear energy, the FSDT and LE4 models with full integration have the same trend. • For R β /h = 2, 4 and 10, the in-plane shear energy is less than 1%, which can be neglected (see Fig. 11c); as the radius-to-thickness ratio increases, the E ps becomes more significant and is plotted for comparison in Fig. 10. • The disagreement of FSDT and LE4 in Fig. 10 is due to that the thickness stretch effects are accounted in LE4 model but ignored by FSDT. When the thickness stretch energy is negligible (less than 1% for R β /h = 50, 100, 500), the transverse shear energy obtained with FSDT and LE4 is quite close, which demonstrates that the kinematic assumptions do not affect the shear locking in thin shells. • The fully integrated lower-order hierarchical elements (p = 2, 3, 4) suffer from locking on the thinner shells (R β /h = 50, 100, 500), and this locking can be overcome by increasing the polynomial degree p without using any locking-mitigation techniques. • When the reduced and selective integration schemes are employed, the shear locking phenomenon on the elements with p = 2 can be greatly alleviated. However, these techniques become less influential when the polynomial degree is further increased, as shown in Fig. 10d-f. Since the newly introduced shape functions lead to improved accuracy, the higher the polynomial degree is, the less necessary the reduced integrated polynomials will become. This effect is more evident for selective integration. • It should be pointed out that, models with these three integration schemes will converge to comparable solutions when the polynomial order is sufficiently high. • Figure 10 clearly demonstrates that shear locking is the dominant locking phenomenon for this numerical example.
The variation of the energy components with the radius-to-thickness ratio R β /h is summarized in Fig. 11. This variation provides a comprehensive understanding of the structural responses when the shell thickness decreases. It can be observed that the membrane energy ratio keeps increasing monotonically with the reduction of the shell thickness, and the ratios of transverse shear energy and thickness stretch energy decrease and approach zero when the shell is very thin (R β /h = 100, 500). In general, the energy ratio of the in-plane shear strains increases when the shell gets thinner. The bending energy is significant for moderate-thin shells. To sum up, the transverse strain energy components (transverse shear and thickness stretch) become less dominant with the decrease of the shell thickness.

Simply supported cylindrical panel under bending
This numerical example consists of three-layered cylindrical panels that undergo bending, as shown in Fig. 12. The cylindrical panels have radius R β = 10, mid-surface arch length b = R β · π/20 in the β direction, and length L = 4.0 along the cylinder axis (α direction). The radius-to-thickness ratios considered include R β /h = 10, 100, 1000, and 10,000. The materials used are the same as in "Three-layered cylindrical shells under distributed pressure" section. The lamination sequence is (0 • /90 • /0 • ), and the thicknesses of the three layers are h 4 , h 2 , and h 4 , separately. As illustrated in Fig. 12b, the cylindrical panels are simply supported on the two ends along the cylinder axis, and free on the other two edges. The simple supports follow: The structure is imposed to constant pressure load p 0 on the top surface. The vertical displacement w is non-dimensionalized through the following parameters: By making use of the symmetric features of the boundary conditions, a 1/4 FE model with one rectangular hierarchical element is employed, as demonstrated in Fig. 12a. The one-element model is refined by increasing its polynomial degree of the hierarchical shape functions until the relative difference of two neighboring orders is less than 0.5% regarding the displacement evaluationw. On the whole, p = 7 is sufficient to guarantee the convergence. The displacement evaluationw and strain components estimation obtained through hierarchical elements with p = 7 are listed in Table 2. Figures 13 and 14  strain energy, respectively. Since no reference is available in the literature, the reference solutions are given by FE models with full integration and p = 8. It can be observed that for the relative thick shells with R β /h = 10 and 100, the reduced integration leads to less accurate displacement estimation and lower convergence compared with full and selective integration schemes. Meanwhile, for the thinner shells with R β /h = 1000 and 10,000, reduced integration models give better accuracy than the other two integration schemes in the lower-order hierarchical shell elements (p = 2 for R β /h = 1000, and p = 2, 3 for R β /h = 10, 000). It seems that reduced integration "unnecessarily" leads to over-soft bending stiffness when no noticeable locking is present in this numerical example. The displacement and energy estimations obtained through selective integration show almost no difference from those of full integration when the element order p ≥ 4. The variation of the relevant energy terms with respect to the element polynomial degree is reported in Fig. 15. For the thick shells with R β = 10, no apparent locking is observed on the low-order hierarchical elements. As the shell becomes thinner, locking becomes significant on low-order elements. For the moderate-thick shell with R β /h = 100, locking occurs on the fully integrated elements for p = 2 and it decreases rapidly as the polynomial degree goes higher. This locking is mainly shear locking, which can be observed through the comparison of bending energy and transverse shear energy in Fig. 15b. On the other hand, for the thin shells with R β /h = 1000 and 10,000, significant membrane locking emerges on the low-order element (p = 2, 3) with full integration, which also drops rapidly with the increase of the polynomial degree. For all the relatively thin shells (R β /h = (a) (b) (c) Fig. 16 Ratio of different energy components versus radius-to-thickness ratio on the three-layered cylindrical panel under bending 100, 1000, and 10,000), reduced integration is able to mitigate the membrane locking effectively which occurs on low-order hierarchical elements, as shown in Fig. 15b-d. The selective integration technique also leads to improved energy estimation, yet less effective compared to reduced integration in this case. To sum up, when the numerical convergence is achieved via p-refinement, all the three integration schemes give equivalent locking-free solutions. From Table 2, it can also be observed that when the hierarchical elements are sufficiently refined, all the integration schemes give results in considerable agreement with each other. The variation of different energy components concerning the increase of the element polynomial order is plotted in Fig. 16. As expected, the strain energy components due to the out-of-plane strains, namely the transverse shear strain energy and thickness stretch energy, decrease rapidly with the increase of the radius-to-thickness ratio and approach zero on thin shells. Also, it can be found that the membrane energy and in-plane shear energy are absent in this numerical example. On the thin shells with R β /h = 1000 and 10,000, the strain energy contains only the bending component.

Conclusions
In this article, the performance of hierarchical elements in the numerical analysis of laminated shell structures is assessed. The numerically efficient integration schemes for hierarchical shell elements are discussed. Proper energy decomposition is proposed and adopted in the evaluation of shear and membrane locking phenomena. Through the numerical investigation with single-element p-version finite element models, the following conclusions can be drawn: • Hierarchical shell elements with full integration are capable of overcoming both shear and membrane locking via plain p-refinement for structures with an aspect ratio in a wide range without using special locking mitigation techniques; • Hierarchical shell elements with polynomial degree p ≥ 3 employing reduced integration scheme and p ≥ 2 adopting selective integration technique are robust; • Reduced and selectively reduced integration schemes are helpful with the lower-order hierarchical elements to alleviate the shear and membrane locking phenomena; • When numerical convergence is reached with a higher-order polynomial degree, the three integration schemes, namely the full, reduced and selective integration, lead to the same accuracy; • With hierarchical shell elements, the physically zero transverse shear strain energy on thin shells can be achieved numerically.
The above evaluation demonstrates the high efficiency of hierarchical shell elements in the analysis of laminated shell structures. Only rectangular elements are tested in the present article. Distorted meshes should also be considered in the future.