Development of shear locking-free shell elements using an enhanced assumed strain formulation

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng 2002; 53:1721–1750 (DOI: 10.1002/nme.360) Development of shear locking-free shell elements using an enhanced assumed strain formulation Jos4eM.A.C4esar de S4a∗;†, Renato M. Natal Jorge‡, Robertt A. Fontes Valente§and Pedro M. Almeida Areias Department of Mechanical Engineering;Faculty of Engineering;University of Porto; Rua Dr. Roberto Frias;4200-465;Porto;Portugal SUMMARY The degenerated approach for shell elements of Ahmad and co-workers is revisited in this paper. To avoid transverse shear locking e;ects in four-node bilinear elements, an alternative formulation based on the enhanced assumed strain (EAS) method of Simo and Rifai is proposed directed towards the transverse shear terms of the strain <eld. In the <rst part of the work the analysis of the null transverse shear strain subspace for the degenerated element and also for the selective reduced integration (SRI) and assumed natural strain (ANS) formulations is carried out. Locking e;ects are then justi<ed by the inability of the null transverse shear strain subspace, implicitly de<ned by a given <nite element, to properly reproduce the required displacement patterns. Illustrating the proposed approach, a remarkably simple single-element test is described where ANS formulation fails to converge to the correct results, being characterized by the same performance as the degenerated shell element. The adequate enhancement of the null transverse shear strain subspace is provided by the EAS method, enforcing Kirchho; hypothesis for low thickness values and leading to a framework for the development of shear-lockingfree shell elements. Numerical linear elastic tests show improved results obtained with the proposed formulation. Copyright ?2001 John Wiley & Sons, Ltd. KEY WORDS: shell elements; enhanced assumed strain method; transverse shear locking; enhanced transverse shear strain <eld 1. INTRODUCTION In the <nite element analysis of shell structures, a landmark can be established with the work of Ahmad et al. [1] leading to the so-called degenerated shell element concept. This ∗Correspondence to: Jos4e M. A. Cesar de S4a, Department of Mechanical Engineering;Faculty of Engineering; University of Porto;Rue Dr. Roberto Frias, 4200-465, Porto, Portugal †E-mail: [email protected] ‡E-mail: [email protected] §E-mail: [email protected] Contract=grant sponsor: MCT-FCT Portugal; contract=grant numbers: PRAXIS XXI=BD=21662=99 and PRAXIS XXI=BD=18533=98 Contract=grant sponsor: FEDER; contract=grant number: POCTI=33640=EME=2000 Received 2 October 2000 Copyright ?2001 John Wiley & Sons, Ltd. Revised 5 March 2001 1722 J. M. A. C4 ESAR DA S 4 AET AL. formulation is obtained by judiciously imposing some chosen kinematic and mechanical assumptions on the three-dimensional continuum, circumventing the use of classical thin shell theories [2;3]. In Reference [2] Bucalem and Bathe highlighted as some key features of this class of elements the applicability to any shell geometry, the use of C0conforming displacement-based elements and the adoption of ‘engineering’ degrees of freedom such as displacements and rotations. However, it is well known that degenerated formulations possess strong de<ciencies in reproducing the behaviour of thin shell structures, leading to locking phenomena. For lowerorder elements, like the four-node bilinear element based on a displacement formulation, these e;ects are responsible for the complete deterioration of the results for thin shells. The shear locking, a;ecting these elements, can be traced to the fact that no Kirchho;– Love hypothesis is introduced into the displacement interpolations. The transverse shear strains cannot then vanish at all points in the element when it is subjected to a pure bending state [4;5]. Amongst the approaches to overcome this problem, the simplest and earliest of all was the use of uniform or selective reduced integration, as <rstly reported in References [6;7]. In general [2], the uniform reduced integration (URI) leads to spurious zero energy modes although in some cases a correct solution is obtained. For the selective reduced integration (SRI) the same problems are posed, but frequently to a smaller degree. In the speci<c case of plate elements, the work of Belytschko and Tsay [8] pioneered a range of contributions using stabilization procedures for controlling the kinematics modes. Other techniques were proposed to alleviate the excess of sti;ness characterizing locking patterns, most of them being classi<ed as hybrid, mixed or hybrid-mixed methods. For these formulations one or more independent <eld assumptions for strains, stress and=or incompatible displacements can be assumed, in conjunction with the usual displacement approximation. The additional <elds may be discontinuous from element to element so that the corresponding unknowns can be eliminated at the element level [5;9]. These methods include, among others, procedures such as the assumed natural strain (ANS) approach (for bilinear shell elements, see References [10–12]; for higher-order elements see, for example, References [13–15]) and the enhanced assumed strain (EAS) formulation [5;16–19]. The key idea of the ANS method is the replacement, in the minimum potential energy principle, of selected displacement-related strains by independently assumed strain <elds in element natural co-ordinates [20]. The variational basis can be found, for example, in the work of Militello and Felippa [9] and the performance of the elements derived from this formulation relies strongly in choosing correctly the sample strain positions for the assumed interpolation. In addition, some de<ciencies in the case of irregular meshes can be pointed out [5]. Nevertheless, the bilinear MITC4 element of Dvorkin and Bathe [10], falling in that category, is amongst the most employed shell formulations, ‘which sometimes cause diOculties for other assumed strain techniques’ [21]. The EAS method, as <rstly proposed by Simo and Rifai [16], uses a three-<eld mixed functional in terms of displacements, stresses and an enhanced strain <eld, based on the Hu–Washizu variational principle. As discussed in this reference, the total strain <eld is built up as a summation of the (compatible) symmetric gradient of the displacement <eld and the enhanced strain. This latter is not subjected to any interelement continuity requirement, and can be related to an ‘incompatible mode’ <eld. Providing a general formulation, Simo and Rifai presented a bilinear isoparametric plate element with enhanced transverse Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 DEVELOPMENT OF SHEAR-LOCKING-FREE SHELL ELEMENTS 1723 shear strain <eld. For a square con<guration this element proved to be identical to MITC4 element [16]. A deeper application of the enhanced techniques in shell analysis came shortly after with the work of Andel<nger and Ramm [5] in the linear elastic range. The extension to non-linear aspects was carried out in subsequent publications (see, for instance, References [17;18;22;23]). Relating the low performance of MITC4 in the presence of distorted meshes to in-plane de<ciencies, Andel<nger and Ramm developed an enhanced membrane and bending strain <eld composed of seven internal parameters. Despite the improvements obtained, their procedures were not extended to shear strain terms, as they kept the assumed natural strain approach of MITC4 in order to attenuate transverse shear locking e;ects. Within a di;erent philosophy, the present work develops an approach suitable for the treatment of the transverse shear locking phenomenon, relying on the enhanced assumed strain formulation and applied over the degenerated four-node element. Following the methodology of the original work of C4esar de S4a and Owen [24] and, subsequently, C4esar de S4a and Natal Jorge [25] for two-dimensional incompressible problems, shear locking appearance is related to the inability of a basis of the null transverse shear strain subspace, implicitly de<ned by each element formulation, to accurately reproduce the required deformation patterns. A deep analysis of MITC4 and selective integrated (SRI) formulations, leading to possible bases for their null transverse shear strain subspaces, reveals the missing terms on the respective basis for the degenerated element. Performing an enhancement over the covariant shear strain terms of this element, an enlargement of its subspace is obtained, leading to a new class of degenerated shell elements with improved behaviour. The speci<c analysis of the null transverse shear strain subspace carried out in this work also revealed that MITC4 element lacks two missing components in its basis, when compared to SRI elements. This same component is not present in the basis of the degenerated element. In fact, in trying to reproduce the deformation pattern corresponding to this absent component, both elements reveal a strong deterioration of results. In this speci<c, single element, example MITC4 and Ahmad et al. formulations provide the same numerical result, while the SRI approach shows no shear locking. It is worth noting the important insight over locking phenomena brought forward by the work of PitkParanta et al. [26] and Chapelle and Bathe [27]. Based on a detailed analysis of bending and membrane subspaces, those authors developed a class of benchmark tests for the behaviour characterization of shell elements, which will be used in this paper (for details see References [28–30]). An outline of the paper is as follows. In Section 2, the analysis of the null transverse shear strain subspace implicitly de<ned for each <nite element formulation is showed. Three architectures of shell <nite elements are described: the original degenerated approach (applied to the bilinear four-node element), the ANS formulation for MITC4 element and the SRI method. Di;erent subspace dimensions are then obtained for each element, justifying the different performances of each formulation in dealing with transverse shear locking phenomena. Based on these results, in Section 3 suitable enhanced interpolation matrices are proposed, attenuating shear locking e;ects for thin shells. Acting over the four-node degenerated element a class of new elements, with an enhanced transverse shear strain <eld is presented. Following an alternative track a di;erent proposal relying on the enhancement of the MITC4 element Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 1724 J. M. A. C4 ESAR DA S 4 AET AL. and grouping EAS and ANS techniques is proposed. Finally, Section 4 provides numerical test problems stating the performance of the presented elements when compared to previous published results for well-known shell elements available in the literature. 2. BASES FOR THE SUBSPACE OF NULL TRANSVERSE SHEAR STRAINS 2.1. Introduction In the analysis of shear locking phenomena, attention is focused on the transverse components of the strain <eld, de<ned in the natural set of co-ordinates (; ; ). Related to this system, the displacement dependent strain <eld for a given point within an element can be represented in a discrete vector form as E≡Eu=Eu  Eu  Eu  Eu  Eu  Eu  T≡Qdn(1) being Qa matrix grouping di;erential terms, and relating nodal displacements dnwith strains. According to the usual representation of degenerated elements [1;31], the nodal displacement vector includes three translational degrees of freedom (de<ned in the global co-ordinate system) and two independent rotational degrees of freedom (de<ned in a local frame) as dn={uvw 12}|T n(2) for a given node (n). The natural strain components described in Equation (1) can be evaluated, for the linear case [31], from the equation Eu ij=1 2@u @i ·gj+@u @j ·gifor (i; j =1;2;3; 1=; 2=; 3=)(3) based on the displacement vector uand the covariant base vectors gi=@x @i (4) with xbeing the position vector. The importance of the described natural strain <eld is related to the imposition of the Kirchho; constraints for thin shells, requiring that those strains should vanish as thickness values tends to zero. When applying the Mindlin-type hypothesis, this condition cannot be satis<ed at all points over a speci<ed domain. However, in an ‘average’ sense [32] this can be accomplished as follows: SA Eu  dd= 0 (5) SA Eu  dd= 0 (6) where SArepresents the particular area analysed. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 DEVELOPMENT OF SHEAR-LOCKING-FREE SHELL ELEMENTS 1725 The assumption that the shear strain <eld Eu  is a linear function of and, in a similar way, Eu  is a linear function of [32] leads to the equations 0 −1 Eu  d=0;+1 0 Eu  d= 0 (7) 0 −1 Eu  d=0;+1 0 Eu  d= 0 (8) The imposition of these constraints discretely at the Gauss points leads to E E=Qshear de=0(9) with deencompassing the nodal degrees of freedom within a speci<ed element. Following previous works of the authors [24;25;33] it is possible to de<ne Is h={uh∈Uh:Qshear uh=0}(10) where Is his the subspace of possible deformations associated with a null transverse shear strain <eld, and Uhrepresents an approximation of the linear space of all admissible solutions U. To avoid the occurrence of transverse shear-locking e;ects, the discrete <nite element solution uh≡demust lie on the null space of Qshear. In other words, any required displacement solution must be a linear combination from the components of a given basis of Is h[25;33]. If this is not the case, the onset of shear locking is veri<ed. Distinct <nite element formulations, for a given mesh, will lead to distinct null transverse shear strain subspaces. In the next sections an analysis over the bilinear degenerated, mixed interpolated MITC4 and bilinear selective reduced elements is performed, aiming to identify for each one of the elements their respective subspace Is hand, consequently, to clarify their ability (or not) to avoid locking behaviours. 2.2. Analysis of the degenerated element Referring to Equation (9), it is possible to state for a single element (dropping the elemental indices) that E E=Q Q d=0(11) where the former matrix Qshear was decomposed into two sub-matrices Q and Q. In the following analysis, a square element placed in the xy plane is considered. The normal to the surface of the element is taken coincident with the z-axis, which is also responsible for thickness orientation. The element represents then a reference square in the natural co-ordinate system, with representative dimensions Dref =[−1;+1] ×[−1;+1] = {(; )∈ 2;−16; 6+1}(12) Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 1726 J. M. A. C4 ESAR DA S 4 AET AL. With the following con<guration, the in-plane {uv}components of the displacement vector (2) do not interfere with the transverse shear strain <eld. It is obvious that any deformation in which only the in-plane degrees of freedom are active will belong to Is h. Therefore, for the sake of simplicity, these degrees of freedom will be ‘removed’ from the generalized nodal displacement vector (2), leading to the ‘reduced’ vector dr={wn=1 1|n=1 2|n=1 (:::)wn=4 1|n=4 2|n=4}T(13) speci<c to the bilinear element and associated with the condition E E=Q Q dr=0(14) in replacement of (11). The next step refers to the extension of this equality over all the sample points used in the numerical integration. For the four-node bilinear element, with a complete 2 ×2 Gauss integration, 8 restrictions in the form of (14) are introduced and two expressions can be obtained as                E(I;I) E(II;II) E(III;III) E(IV;IV)                =       Q(I;I)T Q(II;II)T Q(III;III)T Q(IV;IV)T          4×12 dr=0(15) and              E(I;I) E(II;II) E(III;III) E(IV;IV)              =       Q(I;I)T Q(II;II)T Q(III;III)T Q(IV;IV)T          4×12 dr=0(16) where each pair (k; k) represents the natural co-ordinates of a speci<c integration point (k=I;::: ;IV). Also, Q(k; k) and Q(k; k) represent vectors (12 components) referring to the same integration point. It is worth noting that Qshear evaluation does not depend on the numerical integration over the thickness and, therefore, the orientation was not taken into account in the following. Therefore a reduced displacement <eld, in the sense of (13) and satisfying Equations (15) and (16) (simultaneously avoiding the trivial solution dr=0) is searched. After these considerations, the analysis of the degenerated four-node element, characterized by a transverse shear strain <eld derivable from (3), provided a matrix Qshear with the following characteristics: rank(Qshear)=8;therefore nullity(Qshear) = 4 (17) Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 DEVELOPMENT OF SHEAR-LOCKING-FREE SHELL ELEMENTS 1727 Figure 1. Basis components for the null transverse shear strain subspace; degenerated (D) elements (one element mesh). As a consequence, the element can reproduce four fundamental displacement patterns. In case a given displacement cannot be obtained as a linear combination of these four terms, transverse shear locking is veri<ed. The displacement con<gurations forming a basis of the null transverse shear strain subspace for the degenerated element Is h|Dare schematically represented in Figure 1. The arrows in the <gure refer to activated degrees of freedom. The displacement vector for the whole element (following the representation in Equation (13)) can be represented as follows: D1 ⇒dr={1;0;0;1;0;0;1;0;0;1;0;0}T D2 ⇒dr={0;1;0;−1;1;0;−1;1;0;0;1;0}T D3 ⇒dr={−1;0;1;−1;0;1;0;0;1;0;0;1}T D4 ⇒dr={1;2;−2;−1;2;2;1;−2;2;−1;−2;−2}T (18) The element is able, then, to reproduce these displacement patterns and also any linear combinations between them. 2.3. Analysis of the mixed interpolated bilinear element The described approach is now applied into the four-node bilinear element of Dvorkin and Bathe [10]. The concept behind the element formulation is well established in the <nite element literature (see, for instance, References [2;10;12;14;31] for details) being only summarized here. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 1728 J. M. A. C4 ESAR DA S 4 AET AL. The assumed natural strain approach (or, as originally stated, mixed interpolation of tensorial components) applied to the four-node element consists in the reduction of the bilinear transverse shear <elds to a constant linear one, based on sampling points for each of the two shear strain components [5]. In a general sense, the interpolated transverse shear strains can be de<ned as EANS  (; ; )= 2  p=1 N p(; )Eu (p; p;0) (19) EANS  (; ; )= 2  q=1 Nq(; )Eu (q; q;0) (20) function of the displacement evaluated strains of Equation (3). In the latter equations, (p) and (q) represent additional sampling points with natural co-ordinates (p; p) and (q; q), along with their respective shape functions (N p) and (Nq). The analysis of Equations (15) and (16) for this element reveals a matrix Qshear with the following characteristics: rank(Qshear)=4;therefore nullity(Qshear) = 8 (21) The null transverse shear strain subspace Is h|Mfor the mixed interpolated MITC4 element should, then, be represented by a basis with 8 terms, being a possible candidate represented in Figure 2, with the respective displacement vectors as follows: M1 ⇒dr={1;0;0;1;0;0;1;0;0;1;0;0}T M2 ⇒dr={0;1;0;−1;1;0;−1;1;0;0;1;0}T M3 ⇒dr={−1;0;1;−1;0;1;0;0;1;0;0;1}T M4 ⇒dr={1;2;−2;−1;2;2;1;−2;2;−1;−2;−2}T M5 ⇒dr={0;1;0;0;−1;0;0;0;0;0;0;0}T M6 ⇒dr={0;0;0;0;0;0;0;−1;0;0;1;0}T M7 ⇒dr={0;0;−1;0;0;0;0;0;0;0;0;1}T M8 ⇒dr={0;0;0;0;0;−1;0;0;1;0;0;0}T (22) It can be noted that the <rst four terms of MITC4 subspace basis are exactly the same of the degenerated element. The remaining four terms do not belong to the null transverse shear strain subspace of the degenerated element. 2.4. The selective reduced element Starting with the degenerated element, an additional analysis is performed over the selective reduced integration procedure [6;7]. For this purpose, a reduced integration (one point quadrature) is adopted for the calculation of the transverse shear strain terms. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 DEVELOPMENT OF SHEAR-LOCKING-FREE SHELL ELEMENTS 1729 Figure 2. Basis components for the null transverse shear strain subspace; MITC4 (M) elements (one element mesh). The analysis of this formulation with Equations (15) and (16) leads to a matrix Qshear characterized by rank(Qshear)=2;therefore nullity(Qshear) = 10 (23) Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 1736 J. M. A. C4 ESAR DA S 4 AET AL. Using the enhancement terms obtained from the subspace analysis, it was veri<ed that an improved performance can be obtained including in (32) the cross derivatives of the bubble function in order to and , leading to a new interpolation matrix of the form M(; )|6=    @N @ 0@N @ 0@N @ @N @ 0 0@N @ 0@N @ 0@N @ @N @    (34) The matrix M(; )|6can then be related, within the context of the present work, to an optimal choice for the enhanced interpolation, with the internal variable <eld now being represented as a vector with six components. Similarly to before, the elements obtained from this formulation can be termed S4E6 and S4E6 P7, as only shear or shear plus in-plane enhancement is present, respectively. 3.3. A di0erent approach:the MITC4 E2 shell element Changing the approach adopted herein, the MITC4 element of Dvorkin and Bathe [10] is now focused. As detailed in the analysis of the null transverse shear strain subspace described in Section 2.3, the chosen basis of the element has two missing components when compared to SRI formulation, namely the (S9) and (S10) terms of Figure 3. As detailed in the beginning of this section, these two components can be included into the deformation <eld with an enhanced process using two internal variables (Equation (30)). In this way, it is possible to improve the performance of the MITC4 shell element via the EAS method over the transverse shear strains. The shell <nite element so obtained is termed MITC4 E2, following the previous notation, and possesses a total transverse strain <eld in the form E E =EANS  EANS  +U E U E2 (35) with the assumed natural strain terms being as in Equations (19) and (20) and the interpolation matrix used in the enhanced strain <eld de<ned in Equation (30). As it would be expected, a new analysis of Example 2:2 using this element provides the same result as the SRI formulation. 4. EXAMPLES FOR ELEMENTS ASSESSMENT In order to assess the performance of the proposed shell elements, several problems selected from the literature are analysed. The obtained results are compared to those previously published using well-known shell element formulations. A list of the analysed shell elements and the abbreviations used to identify them henceforth is contained in Table I. 4.1. Morley’s 30◦Skew plate The skew, simply supported plate originally proposed by Morley [34] is analysed in order to test the sensitivity of the proposed elements to mesh distortions. The geometric, material and load data (L= 100;h=1;E=10 5;=0:3, uniform pressure load Q= 1) are based on the work of Andel<nger and Ramm [5], being represented in Figure 7. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 DEVELOPMENT OF SHEAR-LOCKING-FREE SHELL ELEMENTS 1737 Table I. List of referenced and proposed shell elements. Name Description AHMA4 Bilinear degenerated shell element, based on the approach of Ahmad et al. [1] EAS7 ANS In-plane bilinear enhanced shell element; Andel<nger and Ramm [5] MITC4 Mixed interpolated bilinear shell element; Dvorkin and Bathe [10] SIMO 1990 Enhanced shear formulation, bilinear element; Simo et al. [16] Present formulations S4E4 Bilinear degenerated shell element; enhanced transverse shear formulation with 4 internal variables S4E4 P7 Bilinear degenerated shell element; enhanced transverse shear formulation with 4 internal variables and enhanced in-plane formulation of Andel<nger and Ramm [5] S4E6 Bilinear degenerated shell element; enhanced transverse shear formulation with 6 internal variables S4E6 P7 Bilinear degenerated shell element; enhanced transverse shear formulation with 6 internal variables and enhanced in-plane formulation of Andel<nger and Ramm [5] MITC4 E2 Original mixed interpolated element from Dvorkin and Bathe [10]; enhanced transverse shear formulation with 2 internal variables Figure 7. Morley’s skew plate. Still following this work, the Kirchho; reference solution of 4.455 presented by Morley is replaced by the value 4.640, as even for the length to thickness ratio of 100 shear deformation e;ects cannot be neglected. The present results are compared with the ones obtained by Andel<nger and Ramm (for their EAS7 ANS and Dvorkin and Bathe MITC4 elements). Additionally, the results obtained with the enhanced plate element proposed by Simo and Rifai [16], as implemented by Yunhua and Eriksson in a recent work [20] are showed. The values for the central point deVection of the plate are represented in Table II and graphically in Figure 8. For the sake of completeness, the results obtained with the four-node degenerated shell element are reproduced. The elements with 6 enhanced variables (S4E6 and S4E6 P7) and also the MITC4 E2 element show remarkably good results even for coarse meshes, leading to the same numerical values. The elements with 4 enhanced variables (S4E4 and S4E4 P7) reveal some sensitivity to coarse meshes, converging to the performance of the other elements as the mesh is re<ned. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 1738 J. M. A. C4 ESAR DA S 4 AET AL. Table II. Central point (C) deVection of Morley’s skew plate. Mesh MITC4 EAS7 ANS SIMO 1990 AHMA4 Ref. solution 4×4 3.9182 4.2122 3.9305 0.1151 4.6400 8×8 3.8991 4.2239 3.9841 0.3841 4.6400 16 ×16 4.1875 4.3738 4.2727 0.9882 4.6400 32 ×32 4.4098 4.4827 4.4668 1.9110 4.6400 Mesh S4E4 S4E4 P7 S4E6 S4E6 P7 MITC4 E2 4×4 1.9220 1.9220 4.5090 4.5090 4.5090 8×8 3.4530 3.4530 4.4380 4.4380 4.4380 16 ×16 4.2600 4.2600 4.4820 4.4820 4.4820 32 ×32 4.5310 4.5310 4.5610 4.5610 4.5610 Figure 8. Central point deVection for Morley’s skew plate. 4.2. Scordelis–Lo roof problem The geometry of this cylindrical shell, according to Reference [5], is represented in Figure 9, along with the relevant material properties. This data set is also dealt with in Reference [35], reporting to the original work of Scordelis and Lo [36]. The physical basis of the problem is a deeply arched roof supported only by diaphragms at its curved edges (as an aircraft hanger), deforming under its own weight. The geometry is such that the centre point of the roof moves upwards under the self-weight (directed downward) load. For this problem the vertical displacement at the midpoint of the free edge (point A in the <gure) is taken into account, and a converged numerical solution of 0.3024 is adopted for comparison purposes [35]. Due to symmetry, only one quarter of the structure is considered (dashed zone in the <gure). Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 DEVELOPMENT OF SHEAR-LOCKING-FREE SHELL ELEMENTS 1739 Figure 9. Scordelis–Lo roof. Figure 10. Scordelis–Lo roof. Normalized displacement value at midpoint of free edge (point A). With these considerations, the obtained results are presented graphically in Figure 10, with the normalized values (ratio between <nite elements and analytical solutions) listed in Table III. The results of the presented elements are plotted against the MITC4 element and also in comparison with the EAS7 ANS element. From the <gure and table it is evident the better performance of S4E6 P7 element which, even for the coarser mesh, leads virtually to the correct result. Comparing with EAS7 ANS element, it is patent the improvement obtained with the shear enhancement of the present work in comparison with the assumed natural strain method. Once again it can be veri<ed the good performance of the other enhanced elements proposed, with S4E6 and MITC4 E2 results being almost as good as those obtained with Dvorkin and Bathe element. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 1740 J. M. A. C4 ESAR DA S 4 AET AL. Table III. Scordelis–Lo roof. w(f.e.m.)=w(analytical)—midpoint of free edge Elements 5 nodes per side 7 nodes per side 9 nodes per side MITC4 0.93651 ∗0.97315 EAS7 ANS 1.04068 ∗1.00625 AHMA4 0.06806 0.10152 0.13439 S4E4 0.53472 0.83565 0.92593 S4E4 P7 0.57771 0.87864 0.95602 S4E6 0.90542 0.95866 0.96925 S4E6 P7 1.00099 1.01323 1.00165 MITC4 E2 0.90046 0.95800 0.96892 ∗Not available. Figure 11. Partly clamped hyperbolic paraboloid. One half of the structure meshed with 16 ×8 four-node elements. 4.3. Bending-dominated test I – Partly clamped hyperbolic paraboloid In analysing the asymptotic behaviour of shell elements as the thickness tends to zero, it is possible to de<ne problems as being membrane or bending dominated. When this distinction is not clear the asymptotic deformation state is said to pertain to an intermediate state [26]. Following the analysis of Malinen and PitkParanta [30], bending-dominated problems are known to be hard tests for standard <nite element formulations, which provide a poor approximation for inextensional displacements, with numerical failure likely to occur. A particular interesting bending-dominated test problem was introduced in the work of Chapelle and Bathe [27] and further developed in Reference [28]. The problem consists of a hyperbolic paraboloid shell structure, loaded by self-weight and clamped along one side, as represented in Figure 11. The geometric, material and load data are the adopted in Reference [28], and only one half of the surface needs to be considered in the analysis. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 DEVELOPMENT OF SHEAR-LOCKING-FREE SHELL ELEMENTS 1741 Figure 12. Convergence in strain energy for the hyperbolic paraboloid problem: h=L =1=100. For the <nite element mesh, sequences of N×N=2 elements (N=4;8;16;32;64), where Nis the subdivisions along the X-axis, were adopted. Following the proposal of Bathe et al. [28], thickness–length (h=L) relations of 1=100, 1=1000 and 1=10 000 were considered. For this problem there is no analytical solution, and reference values for the total strain energy, previously obtained by Bathe et al. with a re<ned mesh of high-order shell elements, were used [28]. Graphically, we represent the strain energy error of the <nite element solution in comparison to that reference value, i.e. Er=E−Eh E=1−Eh E(36) with Ehbeing the strain energy of the <nite element approximation and Ebeing its reference value. The performance of the proposed elements can be veri<ed in Figures 12–14, according to the thickness values adopted. The results of MITC4 element, as presented in [28], are also reproduced for comparison. From the <gures, it can be seen that the proposed elements performed quite well when compared with the MITC4 element. However, elements S4E4 and S4E4 P7 revealed a strong sensitivity to the decrease in thickness. For the ratio h=L =1=1000 and with the maximum mesh re<nement allowed (N= 64) the convergence error is about 10 per cent, reaching unacceptable values for the relation h=L = 10 000. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 1742 J. M. A. C4 ESAR DA S 4 AET AL. Figure 13. Convergence in strain energy for the hyperbolic paraboloid problem: h=L =1=1000. In contrast the shell element S4E6 P7 obtained the best results for the ratio h=L = 100, with elements S4E6 and MITC4 E2 showing similar results. For ratios h=L of 1=1000 and 1=10 000 these elements revealed an excess of Vexibility with the lower re<nement meshes (N= 4 and 8), showing however an improved behaviour as the mesh is re<ned. 4.4. Bending-dominated test II – Free cylindrical shell In this example, a cylindrical shell of uniform thickness h, length 2Land radius Ris loaded by an axially constant pressure distribution. The pressure varies with the angular distance from the apex (see Figure 1 in Reference [28]) following the function p(’)=P0cos(2’) (37) where the amplitude factor P0is scaled with the thickness as P0=p0h3(38) being p0a constant. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 DEVELOPMENT OF SHEAR-LOCKING-FREE SHELL ELEMENTS 1743 Figure 14. Convergence in strain energy for the hyperbolic paraboloid problem: h=L =1=10 000. The cylinder possesses free ends with no kinematical constraints, being self-supported by the load pattern. The geometry, material and load data are the ones provided by the study of Bathe et al. in Reference [28]. This problem has been extensively described by the group of Juhani PitkParanta, with the publication of an analytical locking insight for various shell models early in 1995 [26]. The reference solutions for the test case are the ones published in the paper of Malinen and PitkParanta [30]. In this simulation, 1=8 of the structure were modelled by uniform mesh patterns of N×2N elements, being N(N=4;8;16;32) the number of elements along the axial direction. As for the previous example, thickness-to-length ratios of 1=100, 1=1000 and 1=10 000 were employed, and a strain energy analysis was carried out with the results presented in Figures 15–17. In this example all the proposed elements are performed in a very similar way. For the thickness–length ratio of 1=100 the MITC4 element behaviour is clearly superior, as for the most re<ned mesh (32 ×64 elements) the proposed elements presents an error of 1.0 per cent, while Dvorkin and Bathe element leads almost to the correct solution. However, as thickness values decrease, all elements tend to converge to the same performance. It is particularly interesting to note that even the 4 enhanced parameters elements S4E4 and S4E4 P7, in Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 1744 J. M. A. C4 ESAR DA S 4 AET AL. Figure 15. Convergence in strain energy for the free cylinder problem: h=L =1=100. Figure 16. Convergence in strain energy for the free cylinder problem: h=L =1=1000. opposition with the previous example, converged to a correct value for the lower thickness values. 4.5. Membrane-dominated test I – Clamped cylindrical shell The geometry and loading for this example are the same as for the previous one, but now the cylinder presents clamped ends. This modi<cation leads to a membrane-dominated problem, where pure bending is inhibited but still conducting to a well-posed problem (References [26–29]). Although the load pattern remains the same, it is scaled in a way consistent of the membrane characteristic of the asymptotic behaviour. Doing so, Equation (38) is replaced by P0=p0h(39) that is, the load amplitude now being a linear function of thickness values. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750 DEVELOPMENT OF SHEAR-LOCKING-FREE SHELL ELEMENTS 1745 Figure 17. Convergence in strain energy for the free cylinder problem: h=L =1=10 000. Figure 18. Convergence in strain energy for the clamped cylinder problem: h=L =1=100. Using the same mesh con<guration as for the free cylinder problem before, the strain energy is evaluated for every element presented and for the ratios thickness to length described before. Reference solutions come from PitkParanta et al. [26;30] and MITC4 results, as presented by Bathe et al. [28], are reproduced and the <nal results are indicated in Figures 18–20. As thickness values tend to zero, the strain energy error of all elements increases. Nevertheless elements MITC4 E2, S4E6 and S4E6 P7 present similar results as Dvorkin and Bathe element. For the coarser mesh (N= 4) MITC4 element provides the best performance with di;erences between this element and the proposed ones vanishing as re<nement is increased. It is also noticeable the lower convergence rate of S4E4 and S4E4 P7 elements. Copyright ?2001 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2002; 53:1721–1750