A serendipity plate element free of modeling deficiencies for the analysis of laminated composites

Composite Structures 154 (2016) 150–171

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

A serendipity plate element free of modeling deficiencies for the analysis of laminated composites João Elias Abdalla Filho a,⇑, Ivan Moura Belo b, John O. Dow a a b

Department of Civil, Environmental and Architectural Engineering, University of Colorado, Boulder, CO, USA Mechanical Engineering Program, Federal Technological University of Paraná (UTFPR), Curitiba, Paraná, Brazil

a r t i c l e

i n f o

Article history: Received 29 May 2016 Accepted 19 July 2016 Available online 25 July 2016 Keywords: Laminated composite plates Serendipity plate element Strain gradient notation Shear locking Spurious zero-energy modes Parasitic shear

a b s t r a c t An eight-node serendipity element free of shear locking and spurious zero-energy modes is formulated to model laminated composite plate problems. The element is based on a first-order shear deformation theory and on the equivalent lamina assumption. Stresses are calculated through the thickness of the plate. As the model is only capable of representing transverse shear strains and stresses as constants, while their actual variations are parabolic, a shear correction factor is used. The element is formulated using strain gradient notation, which is a physically interpretable notation that allows for a detailed a-priori analysis of the finite element model. The element’s shear strain polynomials are inspected, and the spurious terms which are responsible for shear locking are identified. The element is corrected by simply removing such spurious terms from those shear strain expansions. Further, the compatibility modes are also clearly identified and maintained in the shear strain expansions in order to prevent the introduction of spurious zero-energy modes. Numerical results show the shear locking effects caused by the spurious terms on displacement and transverse stress solutions. They also show that properly refined meshes composed of corrected elements provide solutions which converge rather well for moderately thick to thin plates. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Research on modeling and analysis of laminated composite structures has been very active in at least the last four decades. The concern to accurately represent the actual behavior of this kind of structure has led researchers to develop analytical and numerical models every time more refined. Thorough reviews on theories and finite element models for laminated composites have been presented in the literature [1–4]. Carrera [5] discusses theories for laminated composites in the perspective of the correct representation of the continuity of the transverse displacement and the equilibrium of interlaminar transverse stress components. The author indicates that layerwise models are necessary for an accurate evaluation of the normal transverse stresses. Another work presents a review on the different methods, both analytical and numerical, available for estimating transverse and interlaminar stresses in laminated composite plates and shells [6]. An early

⇑ Corresponding author at: Mechanical Engineering Program, Pontifical Catholic University of Paraná, Brazil, and Civil Engineering Program, Federal Technological University of Paraná, Brazil. E-mail address: [email protected] (J.E. Abdalla Filho). http://dx.doi.org/10.1016/j.compstruct.2016.07.042 0263-8223/Ó 2016 Elsevier Ltd. All rights reserved.

attempt to formulate a finite element for laminated composites made by Mawenya and Davies [7] resulted is an element based on first-order shear deformation theory (FSDT) with rotations varying from layer to layer. Another early example is the work by Panda and Natarajan [8] where it is presented an element based on the eight-node degenerate shell developed by Ahmad et al. [9]. As the first-order shear deformation theory is viewed as a limiting theory in representing accurately the behavior of laminated composites, several other works have proposed high-order deformation theories (HSDT) [10–14]. An overview of the relationships between classical and shear deformation theories has been presented [15]. Computational models ranging from simple to refined have been developed to perform numerical evaluation of all those theories [14,16]. The concern with the development of accurate laminated composite finite element models obviously has included avoiding shear locking and other spurious mechanisms. Isoparametric elements require some form of reduced-order integration of the stiffness matrix. Mawenya and Davies [7], for instance, followed the recommendation of using a 2
2 Gaussian integration provided by Zienkiewicz et al. [17]. It was not known at the time that such a strategy would not be very efficient as it would be proved later by other researchers. The works of Hughes and

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

co-workers with different plate elements shed light over the shear locking problem as well as spurious zero-energy modes (rank deficiency) and made significant progress in developing selective reduced-integration procedures for isoparametric elements [18–20]. However, as selective reduced integration of Lagrange elements introduces spurious zero-energy modes, the Heterosis element was devised to circumvent such problem at the expense of a more complex formulation [21]. Alternative formulation procedures have been devised seeking the development of better plate and shell elements. Among those we cite the field consistency approach [22] and the mixed interpolation of tensorial components [23]. Laminated composite shell elements have been developed using the latter [24]. Also, an approach of adding a shear correction factor to transverse shear strain expressions developed by Verwoerd and Kok [25] has been applied to three-node triangular laminated composite plate elements based on FSDT to suppress shear locking effects [26]. Zhang and Yang [27] introduce triangular and quadrilateral plate and shell elements for linear and nonlinear analysis of thin to moderately thick laminates. Such elements are shear locking-free due to the use of Timoshenko’s beam functions to represent deflections and rotations of the elements sides. A C1 six-node triangle is developed using trigonometric functions to represent transverse shear stresses. Polynomial orders for the transverse displacement and for the rotations are selected such that transverse shear strain compatibility is attained, thus avoiding shear locking naturally [28]. Kim et al. [29] formulate a FSDT shell element within the co-rotational framework for geometrically nonlinear analysis. The ANS (Assumed Natural Strains) method is employed, precluding shear locking effects, and thus allowing good behavior for thin and thick problems. Further advances include layerwise models [30], and elements based on Reddy´s simple higherorder shear deformation theory [31,32]. Moleiro et al. [33] develop a layerwise finite element model based on a mixed least-squares formulation which enables interlaminar continuity of displacements and transverse stresses in the form of Co continuous functions. Mantari and Guedes Soares [34] introduce a generalized layerwise displacement-based HSDT and the corresponding finite element in which the number of DOFs is limited for being independent of the number of layers. Pandey and Pradyumna [35] present a C° layerwise eight-node finite element model for three-layered composite plates where first-order displacement fields are assumed for the top and bottom layers while a higher-order displacement field is assumed for the middle layer. Also, there are contributions made using the Refined Zigzag Theory (RZT), which has been presented by Tessler et al. [36] recently. RZT is based on FSDT and employ piecewise-linear zigzag functions that provide better representations of the deformations states of transverse shear-flexible plates. Those authors [37] present two- and a three-node beam elements that avoid shear locking effects by employing anisoparametric interpolations to approximate the kinematic variables that are necessary to model planar deformations. The same group [38], following the same line of the previous referenced work, formulate RZT-based six- and three-node triangular plate elements, which also employ shape functions based on anisoparametric interpolations. Further, Eijo et al. [39] present the formulation of an isoparametric four-node C° quadrilateral plate element based on the RZT theory. Shear locking is avoided by using an assumed linear strain shear field. Finally, another approach is the one adopted by Natarajan et al. [40] which combines Carrera’s Unified Formulation (CUF) and the cell-based smoothed finite element method. A four-node quadrilateral based on the field consistency requirement is employed to eliminate shear locking. This paper describes the development of a FSDT-based eightnode serendipity plate finite element using strain gradient notation and makes an assessment of its capabilities in the analysis

151

of laminated composite plate problems. Spurious terms which are the cause of shear locking are identified and removed. 2. On the strain gradient notation Strain gradient notation (Dow [41]) is a physically interpretable notation which relates displacements to the kinematics quantities of the continuum. It has been used as an alternative procedure for formulating finite elements and also finite difference templates. The relationships between displacements and the kinematics quantities are derived through a procedure which identifies the physical contents of the polynomial coefficients. An important feature of strain gradient notation is that the modeling characteristics of the finite element become apparent to the developer since the early steps of the formulation. This allows for the sources of modeling errors to be identified and consequently removed from the finite element shear strain polynomial expansions prior to the formation of the element stiffness matrix. In this paper, sources of shear locking (parasitic shear terms) are precisely identified and removed from the shear strain polynomial expansions of the eight-node serendipity plate element. Causes of spurious zero energy modes are explained, and this deficiency is avoided by correctly eliminating shear locking. Parasitic or spurious shear terms which are present in the shear strain polynomial expansions of the serendipity plate element are identified by inspection. It is demonstrated both theoretically that they are flexural terms which cause locking of the model by increasing the element´s shear strain energy unduly when the plate undergoes bending. Such terms are then simply removed from the shear strain polynomials, avoiding shear locking to occur. It is also demonstrated that spurious zero-energy modes are not introduced into the model by recognizing and not removing the compatibility modes. Such modes can be easily confused with shear locking terms and, therefore, be inadvertently removed. This is the limitation of reduced-order integration schemes in attempting to correct elements for locking. Along with eliminating legitimate spurious terms which are responsible for locking, those techniques also eliminate compatibility modes, thus introducing spurious zero-energy modes. Hence, the use of strain gradient notation has the advantages that locking is taken care of correctly and a-priori of the formation of the stiffness matrix and of the computer implementation, and that the formulated element is of correct rank since no spurious zero-energy modes are introduced. Formulation associated to strain gradient notation is described in the next two sections. 3. First-order shear deformation theory The kinematic relations according to the Reissner-Mindlin or first-order plate theory may be written as follows:

uðx; y; zÞ ¼ uo ðx; yÞ þ zhx ðx; yÞ

ð1Þ

v ðx; y; zÞ ¼ v o ðx; yÞ  zhy ðx; yÞ

ð2Þ

wðx; yÞ ¼ w

ð3Þ

hx ðx; yÞ ¼

@uðx; y; zÞ @z

hy ðx; yÞ ¼ 

@ v ðx; y; zÞ @z

ð4Þ ð5Þ

where u and v are in-plane displacements along the x and y directions, respectively, w is the out-of-plane (normal to the middle surface) displacement, hx and hy are rotations in the x and y directions, respectively (or around the y and x axes, respectively), and u0 and v0 are the middle surface´s in-plane displacements along the x and y

152

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

directions, respectively. The displacement w is independent of rotations, which allows transverse shear deformation representation. This feature is particularly important for modeling thick plates and laminated composites. Strains are defined as follows:

8 > > > > > > <

8

9

9

8

Xk

fe

egk dXk

gTk ½Qk f

ð7Þ

where k is a typical lamina, n is the total number of laminae, {e}k is the strain vector of lamina k, [Q]k is the constitutive properties matrix of lamina k, and Xk is the volume of lamina k. At this point of the presentation, the foregoing description is written in terms of strain gradient notation. As mentioned in the previous section, strain gradient notation is mainly characterized by relating displacements and strains to the kinematic quantities of the continuum in a physically interpretable fashion [41]. Such kinematic quantities are rigid body modes, strains, and first- and higher-order derivatives of strains, and they are generally referred to as strain gradients. The relations of displacements and strains to strain gradients are given in symbolic form, respectively, by:

fdg ¼ ½/fesg g

ð8Þ

feg ¼ ½T sg fesg g

ð9Þ

where [/] and [Tsg] are the corresponding transformation matrices, and {esg} is the strain gradients vector. This vector contains the set of independent deformation modes that the model is capable of representing. Matrix [/] is comprised of linearly independent vectors, each associated to a strain gradient component, describing a specific deformation pattern of the model. Specific contents of these vectors and matrices are case-dependent, and they will be specified for the eight-node serendipity plate element in the following section. Eqs. (8) and (9) are combined to eliminate vector {esg}, and the resulting expression is substituted into Eq. (7) to yield:

! n Z X  T 1 T T U ¼ fdg ½/ T sg k ½Qk fT sg gk dXk ½/1 fdg 2 k¼1 Xk

ð10Þ

which is an expression of the strain energy for laminated composites in strain gradient notation. The quantity between parentheses is called strain energy matrix and it is represented by [UM]. The elements of its principal diagonal contain the quantities of strain energy associated with the pure strain modes of the laminate. The other elements of the matrix contain the quantities of energy associated with the coupling between the various strain modes. Matrix [UM] may be written as:

½U M  ¼

Z X n Z

Zk

A k¼1

Z k1

½T sg Tk ½Qk ½T sg k dZ k dA

n X ðQ ij Þk ðZ k  Z k1 Þ

Bij ¼

n   1X ðQ Þ Z 2  Z 2k1 2 k¼1 ij k k

ð13Þ

Dij ¼

n   1X ðQ Þ Z 3  Z 3k1 3 k¼1 ij k k

ð14Þ

Aij ¼ K

n X ðQ ij Þk ðZ k  Z k1 Þ

ð15Þ

k¼1

where A is the membrane stiffness, B is the membrane-bending coupling stiffness, D is the bending stiffness, and A⁄ is the transverse shear stiffness and it is subjected to a shear correction factor K. The determination of the shear correction factor for laminates is still an open issue as it depends on lamination scheme, geometry, and material properties [2]. The value 5/6, which is very accurate for homogeneous, isotropic plates, is most commonly used. This value is adopted in the numerical analyses performed in this paper to allow for comparisons to the analytical solutions employed as reference solutions [2]. Finally, as the strain energy in terms of the stiffness matrix is given by:

U¼

1 T fdg ½Kfdg 2

ð16Þ

the general expression of the stiffness matrix in strain gradient notation turns out to be:

½K ¼ ½/T ½U M ½/1

ð17Þ

4. The strain gradient notation serendipity plate finite element The formulation of an eight-node serendipity plate finite element in strain gradient notation is presented in this section. The finite element is shown in Fig. 1. Five degrees of freedom are associated to each node; namely, in-plane displacements uo and vo, the out-of-plane displacement w, and rotations hx and hy. The displacement polynomial expansions in strain gradient notation, according to Eqs. (1)–(5), are obtained by substituting the unknown coefficients by their physical contents. The procedure to identify physical contents of finite element displacement polynomials is described in detailed by Dow [41]. The procedure has also been detailed in the formulation of a four-node plate element by Abdalla Filho et al. [42]. The strain gradient notation displacement polynomials for the present element are presented below:

ð11Þ

where the volume integral has been broken into an integral over the area of the middle surface of the laminate and an integral over its thickness. This line integral is carried out as the sum of the integrals over the thicknesses of the various laminae. The integration limits zk1 e zk represent the bottom and top coordinates of a typical

ð12Þ

k¼1

ð6Þ

where the first vector contains membrane strains and the second vector contains plate bending strains. The strain energy of the laminate is the sum of the strain energies of its laminae. Hence,

Z

Aij ¼

9

ex > u0;x zhx ; x > > > > > > > > > > > > > > > > > > > > > > > ey > zhy ; y v 0;y > > > < < = > = = > cxy ¼ u0;y þ v 0;x þ zðhx ; y  hy ; xÞ > > > > > > > > > > > > > > > > cyz > w; y  hy > 0 > > > > > > > > > > > > > > > > > > : ; : ; ; : cxz 0 w; x þ hx

n 1X U¼ 2 k¼1

lamina k, respectively. The integration over the thickness of the laminate yields its stiffness quantities:

Fig. 1. Eight-node serendipity rectangular plate finite element.

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171





uðx; y; zÞ ¼ ðuÞo þ ðex Þo x þ cxy =2  r y þ ðex;x Þo x2 =2 o     þ ex;y o xy þ cxy;y =2  ey;x y2 =2 þ ðex;xy Þo x2 y=2 o   þ ex;yy o xy2 =2 þ ðhx þ cxz =2Þo z þ ðex;z Þo xz   þ cxy;z  cyz;x þ cxz;y yz=2 þ ðex;xz Þo x2 z=2 o     2 þ cxy;yz  ey;xz y z=2 þ ex;yz o xyz o     þ ex;xyz o x2 yz=2 þ ex;yyz o xy2 z=2

refinement. Proceeding on with the formulation, the elastic strain expansions for the element are obtained in strain gradient notation through application of Eq. (6):

ð18Þ

v ðx; y; zÞ ¼ ðv Þo þ ðcxy =2 þ rÞo x þ ðey Þo y þ ðcxy;x  ex;y Þo x2 =2

ð19Þ

þ ðcyz;y  ey;z Þo y2 =2 þ ðcxz;xy  ex;yz Þo x2 y=2

ð26Þ

cxz ¼ ðcxz Þo þ ðcxz;x Þo x þ ðcxz;y Þo y þ ðcxz;xy Þo xy þ ðcxy;yz Þo y2 =2 ð20Þ

þ ðcyz;xy Þo y2 =2 þ ðex;xz Þo x2 =2  ðey;xz Þo y2 þ ðex;xyz Þo x2 y=2 þ ðex;yyz Þo xy2 =2

hx ðx; yÞ ¼ ðhx þ cxz =2Þo þ ðex;z Þo x þ ðcxy;z  cyz;x þ cxz;y Þo y=2 þ ðex;xz Þo x2 =2 þ ðcxy;yz  ey;xz Þo y2 =2 þ ðex;yz Þo xy ð21Þ

hy ðx; yÞ ¼ ðcyz =2  hy Þo  ðcxy;z þ cyz;x  cxz;y Þo x=2  ðey;z Þo y  ðcxy;xz  ex;yz Þo x2 =2  ðey;yz Þo y2 =2  ðey;xz Þo xy ð22Þ

These expansions contain the set of fourty strain states that the element is capable of representing. These strain states are rigid body modes, constant strains, normal strain gradients and shear strain gradients, which are listed below:

Constant strains :

cyz ¼ ðcyz Þo þ ðcyz;x Þo x þ ðcyz;y Þo y þ ðcyz;xy Þo xy þ ðcxy;xz Þo x2 =2 þ ðey;xyz Þo xy2 =2

þ ðcxz;x  ex;z Þo x2 =2 þ ðcxy;z þ cyz;x þ cxz;y Þo xy=2

Rigid body modes :

ð25Þ

þ ðcxz;xy Þo x2 =2 þ ðex;yz Þo x2 þ ðey;yz Þo y2 =2 þ ðey;xxz Þo x2 y=2

wðx; yÞ ¼ ðwÞo þ ðcxz =2  hx Þo x þ ðcyz =2 þ hy Þo y

 ðey;xxz Þo x2 y=2  ðey;xyz Þo xy2 =2

ð24Þ

þ ðey;xyz Þo y2 z=2

þ ðcxy;xz  ex;yz Þo x2 z=2 þ ðey;yz Þo y2 z=2 þ ðey;xz Þo xyz

þ ðex;xyz Þo x2 y=2 þ ðex;yyz Þo xy2 =2

ey ¼ ðey Þo þ ðey;x Þo x þ ðey;y Þo y þ ðey;xx Þo x2 =2 þ ðey;xy Þo xy þ ðey;z Þo z þ ðey;yz Þo yz þ ðey;xz Þo xz þ ðey;xxz Þo x2 z=2 þ ðey;xyz Þo xyz

þ ðey;xy Þo y2 =2 þ ðex;xyz Þo x2 z=2 þ ðex;yyz þ ey;xxz Þo xyz

þ ðcxy;z þ cyz;x  cxz;y Þo xz=2 þ ðey;z Þo yz

þ ðcyz;xy  ey;xz Þo y x=2

ð23Þ

þ ðcxy;yz Þo yz þ ðex;xy Þo x2 =2 þ ðex;yy þ ey;xx Þo xy

þ ðey;xy Þo xy2 =2 þ ðcyz =2  hy Þo z

2

ex ¼ ðex Þo þ ðex;x Þo x þ ðex;y Þo y þ ðex;yy Þo y2 =2 þ ðex;z Þo z þ ðex;xz Þo xz þ ðex;yz Þo yz þ ðex;xyz Þo xyz þ ðex;yyz Þo y2 z=2

cxy ¼ ðcxy Þo þ ðcxy;x Þo x þ ðcxy;y Þo y þ ðcxy;z Þo z þ ðcxy;xz Þo xz

þ ðey;x Þo xy þ ðey;y Þo y2 =2 þ ðey;xx Þo x2 y=2

þ ðey;xxz Þo x2 yz=2 þ ðey;xyz Þo xy2 z=2

153

ðuÞo ðv Þo ðrÞo ðwÞo ðhx Þo ðhy Þo ðex Þo ðey Þo ðcxy Þo ðcxz Þo ðcyz Þo

Normal strain gradients : ðex;x Þo ðex;y Þo ðex;z Þo ðex;xy Þo ðex;yy Þo ðex;xz Þo ðex;yz Þo ðex;xyz Þo ðex;yyz Þo ðey;x Þo ðey;y Þo ðey;z Þo ðey;xy Þo ðey;xz Þo ðey;xx Þo ðey;yz Þo ðey;xxz Þo ðey;xyz Þo Shear strain gradients : ðcxy;x Þo ðcxy;y Þo ðcxy;z Þo ðcxy;xz Þo ðcxy;yz Þo ðcxz;x Þo ðcxz;y Þo ðcxz;xy Þo ðcyz;x Þo ðcyz;y Þo ðcyz;xy Þo It is expected that the reader who is not familiarized with the strain gradient notation has started to realize by now how unique and powerful the notation is as it reveals the physical contents of the finite element polynomials coefficients, which would be otherwise merely represented by symbols were a conventional notation employed. For instance, the reader can see that the element is capable of representing all rigid body modes and all constant strain states according to the plate element kinematics, which are the only quantities expected to be represented in the limit of

ð27Þ

Examination of these strain expressions will reveal the presence of spurious terms a-priori. A strain expansion should be only composed of terms which are physically related to the strain quantity they represent. Examination of the expansions for normal strains, Eqs. (23) and (24), shows that all coefficients are associated to their respective normal strains. Therefore, all those terms are legitimate as they correctly contribute to the representation of normal strains. They belong to the Taylor series expansions corresponding to those strains. However, examination of the expansions for shear strains, Eqs. (25)–(27), reveals that besides legitimate terms, these expansions contain terms which are related to normal strains. In general, normal strain terms in shear strain expansions are spurious terms which will cause the increase of the shear strain energy of the element unduly when activated during deformation. That is, they will cause shear locking. For instance, shear locking will occur when the plate is subjected to a loading condition capable of activating the eighth term of Eq. (26), i.e. ey;yz . This strain gradient represents the variation of ey along directions y and z (through-thethickness) simultaneously, which is legitimate in Eq. (24), but spurious in Eq. (26) as it is adding to shear strain energy when it should only add to the normal strain energy. Spurious terms are introduced due to the use of incompatible displacement polynomials. In the present case, the membrane displacement polynomials selected are incomplete, third order polynomials. They do not contain the terms x3 and y3. This causes the presence of spurious terms in the in-plane shear strain expansion cxy. Further, the order of the rotation polynomial expansions and the order of the out-of-plane displacement polynomial expansion are the same, which is inconsistent with Reissner-Mindlin plate theory. The reason is that the transverse shear strains are defined as the sum of the rotations and the first derivatives of the out-ofplane displacement (see Eq. (6)). These polynomials would be consistent if the out-of-plane displacement polynomial were one order higher than the order of the rotation polynomials. This inconsistency causes the presence of spurious terms in the transverse shear strain expansions cyz and cxz. Now, based on these arguments and after examining the strain expansions above, one might be

154

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

tempted to remove all normal strain gradient terms from the shear strain expansions to produce a well-formulated element. Such action, however, would remove shear locking effects at the expense of introducing spurious zero-energy modes. Careful scrutiny of the shear strain expansions reveals that not all their normal strain gradient terms are spurious. Some are actually terms that are needed by the element as they are recognized as the terms of compatibility equations of the elasticity theory. That is, they are compatible modes. This will be shown in the following for each shear strain expansion. The expansion of the in-plane shear strain cxy, Eq. (25), contains eight normal strain gradient terms. Four of these terms are needed in the formulation because they give rise to the following compatibility equations:

5. Numerical applications The present element is implemented in LAMFEM [43], a FORTRAN finite element code, in two versions; namely, one version containing the spurious terms and the other after elimination of the spurious terms. In order to assess the performance of the element and to validate the procedure of eliminating the sources of shear locking a-priori, solutions of four selected application problems using both versions of the element are compared. Further, the corrected model is validated by comparing numerical solutions with results obtained from analytical solutions and from isoparametric models solutions when available [2].

cxy;xy ¼ ex;yy þ ey;xx

ð28Þ

5.1. Application problem #1

cxy;xyz ¼ ex;yyz þ ey;xxz

ð29Þ

To start, a thin cantilever plate problem (edge length-tothickness relation a/h = 100) is analyzed which demonstrates with great clarity the severe shear locking caused by the spurious terms present in the eight-node serendipity plate element. The plate is a two-layer angle ply laminate with lamination scheme 45°/+45° subjected to two point loads applied at the corners of the free end. Fig. 2 shows the lay-out of this problem. Five uniform meshes are used to analyze the problem (2
1, 4
2, 8
4, 16
8, and 32
16). Figs. 3(a) and 3(b) are representative of the transverse displacement w solutions of the plate modeled with elements containing the spurious terms and elements corrected for the spurious terms, respectively, corresponding to the finer mesh (512 elements). It can be observed that the former is much stiffer than the latter, evidencing strong shear locking effect. According to the numerical solutions associated to these figures, the maximum displacements calculated for both models are 0.08 m and 0.18 m, respectively, being notorious that shear locking causes the solution to be greatly underestimated. Fig. 3(c) shows the behaviors of the two models with mesh refinement. The numerical solutions containing spurious terms are represented by with/PS while those not containing spurious terms are represented by wout/PS, where PS stands for parasitic shear. Even though a reference solution is not available, the plots show that the solutions without spurious terms converge asymptotically.

and, hence, the only spurious terms present in the in-plane shear strain expansion are ðex;xy Þo ; ðey;xy Þo ; ðex;xyz Þo ; ðey;xyz Þo . Proceeding on with the a-priori analysis, it is seen that the expansion of the transverse shear strain cyz, Eq. (26), contains four normal strain gradient terms, one gradient of the in-plane shear strain and one gradient of the transverse shear strain in the plane xz, which appear to be spurious for not belonging to the Taylor series expansion of cyz. However, according to the following compatibility equation:

cyz;xx ¼ cxy;xz þ cxz;xy  2ex;yz

ð30Þ

the only spurious terms are ðey;yz Þo ; ðey;xxz Þo ; and ðey;xyz Þo . Finally, it is seen that in the expansion of the transverse shear strain cxz, Eq. (27), there are four normal strain gradient terms, one gradient of the in-plane shear strain, and one gradient of the transverse shear strain in the plane yz, which appear to be allien to the expansion. However, according to the following compatibility equation:

cxz;yy ¼ cxy;yz þ cyz;xy  2ey;xz

ð31Þ

the only spurious terms are ðex;xz Þo ; ðex;xyz Þo ; and ðex;yyz Þo . The terms on the left-hand sides of the compatibility equations above belong to the Taylor series expansions of the corresponding shear strains. Consequently, the terms on the right-hand sides of those equations are compatible modes. Thus, those terms must be kept, and the only terms to be eliminated are those recognized as spurious terms. Proceeding on with the elimination, results in the corrected shear strain expansions, which are:

cxy ¼ ðcxy Þo þ ðcxy;x Þo x þ ðcxy;y Þo y þ ðcxy;z Þo z þ ðcxy;xz Þo xz þ ðcxy;yz Þo yz þ ðex;yy þ ey;xx Þo xy þ ðex;yyz þ ey;xxz Þo xyz

ð32Þ

cyz ¼ ðcyz Þo þ ðcyz;x Þo x þ ðcyz;y Þo y þ ðcyz;xy Þo xy þ ðcxy;xz Þo x2 =2 þ ðcxz;xy Þo x2 =2  ðex;yz Þo x2

ð33Þ

cxz ¼ ðcxz Þo þ ðcxz;x Þo x þ ðcxz;y Þo y þ ðcxz;xy Þo xy þ ðcxy;yz Þo y2 =2 þ ðcyz;xy Þo y2 =2  ðey;xz Þo y2

ð34Þ

These shear strain expansions do not contain any spurious terms and the compatible modes have been kept. Removal of the spurious terms ensures that shear locking will not occur, and the presence of the compatible modes prevents the introduction of spurious zero-energy modes. Therefore, these two measures must be taken to guarantee that spurious mechanisms will not occur in the eight-node serendipity plate element during analysis. As all spurious terms are present in shear strain polynomials, they may be referred to as parasitic shear terms, as it has been done in the literature [17,41]. The eight-node plate element will be validated in the next section.

5.2. Application problem #2 The second application problem is a square simply supported (SS1 supports) cross-ply laminated composite plate subjected to a uniform load of value qo = 10 N/m2. The laminate is symmetric with lamination scheme 0°/90°/0°, and with edge lengths equal to 1.0 m (a = b = 1.0 m) as shown in Fig. 4. All layers are of graphite-epoxy with the following mechanical properties: Ex = 175 GPa, Ey = 175 GPa, Gxy = 3.5 GPa, Gxz = 3.5 GPa, Gyz = 1.4 GPa, mxy = 0.25. The plate is solved with edge length-to-thickness ratios a/h = 10, 100, and 1000, and using four uniform meshes; namely, 2
2, 4
4, 8
8, and 16
16, both with elements containing the spurious terms (with/PS), and with corrected elements (wout/PS). Numerical solutions for through-the-thickness stresses are compared to analytical solutions [2] referred to as FSDT in the plots. Stress results are averaged nodal values and they are nondimensionalized through the following relations [2]: ! ! !  2 2 2 h h h h r xx ¼ rxx 2 r yy ¼ ryy 2 sxy ¼ sxy 2 sxz ¼ sxz aqo a qo a qo a qo  h syz ¼ syz ð35Þ aqo In-plane normal stresses rxx and ryy are calculated at the center point of the plate, and in-plane shear stresses sxy (only for a/h = 1000) at a corner node while transverse shear stresses sxz

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

155

Fig. 2. Cantilever plate subjected to concentrated loads at the free-end. Anti-symmetric angle-ply laminate 45°/+45°. Problem definition.

Fig. 3(a). Cantilever plate subjected to concentrated loads at the free-end. Antisymmetric angle-ply laminate 45°/+45°. Transverse displacement solutions with spurious terms.

Fig. 3(c). Cantilever plate subjected to concentrated loads at the free-end. Antisymmetric angle-ply laminate 45°/+45°. Convergence of transverse displacement solutions with and without spurious terms.

Fig. 3(b). Cantilever plate subjected to concentrated loads at the free-end. Antisymmetric angle-ply laminate 45°/+45°. Transverse displacement solutions without spurious terms.

and syz are calculated at borders middle points. Table 1 contains the highest percent error with respect to the analytical solution in each solution.

First, solutions for the thick plate (a/h = 10) are presented. Figs. 5 (a)–5(l) show the results of these analyses. In general, as depicted by these figures, the effects of shear locking although evident, are not very significant (as they are in thin plates as will be shown subsequently) and they appear to be greatly attenuated by mesh refinement. Figs. 5(a) and 5(b) show that both models with and without the spurious terms converge to the analytical solution for rxx, respectively. The model with spurious terms is less accurate for the coarser meshes. The convergence behavior of the models can be seen in Fig. 5(c). Figs. 5(d) and 5(e) present the results for ryy. Contrary to what was expected, the model with spurious terms presents better

156

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

Fig. 4. Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load. Symmetric cross-ply laminate 0°/90°/0°. Problem definition.

Fig. 5(a). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Normal stresses rxx computed through the thickness of the laminate with spurious terms.

solutions for the coarser meshes. Nevertheless, both models tend to the analytical solution within very good accuracy with mesh refinement. Fig. 5(f) shows that both solutions converge from an

upper bound, and that, although less accurate for coarser meshes, the solution provided by the model without spurious terms converges monotonically. Finally, Figs. 5(g)–5(i), and Figs. 5(j)–5(l) show the solutions for the transverse shear stresses sxz and syz, respectively. It is observed that shear locking effects are stronger in the coarser meshes solutions. The corresponding percent errors in Table 1 are significant. Errors for the models with spurious terms are 39.63% and 104.57%, respectively. These errors drop to 17.51% and 48.73%, respectively, after the removal of the spurious terms. Further, it is also observed that mesh refinement attenuates shear locking, and that both models converge within acceptable levels of error to the analytical solution. The finer mesh for sxz solutions contains errors of only 1.05% and 0.60%, while the finer mesh for syz solutions contains errors of 5.41% and 2.64%. The convergence plots of Figs. 5(i) and 5(l) show that solutions for transverse shear stresses, with and without spurious terms, converge asymptotically to the analytical solutions from lower bounds. Next, solutions for a thin plate (a/h = 100) are presented in Figs. 6(a)–6(l), and also in Table 1. In general, these results show that the model containing the spurious terms either has convergence delayed or prevented whereas the model without the spurious terms converges well and quickly to the analytical solutions. Fig. 6(a) shows that the spurious terms delay convergence of the normal stresses rxx, and that even the finer mesh (16
16)

Table 1 Percent errors in nondimensionalized stress components of symmetric cross-ply square plate (0°/90°/0°). Side-to-thickness relations a/h = 10 and 100. Errors (%)  xx r

(0°/90°/0°)

sxz

 yy r

syz

a/h

Meshes

With/PS

Wout/PS

With/PS

Wout/PS

With/PS

Wout/PS

With/PS

Wout/PS

10

2
2 4
4 8
8 16
16

18.10 4.70 1.09 0.05

6.00 1.40 0.25 0.26

1.89 3.32 1.24 0.72

17.45 4.46 1.24 0.72

39.63 15.54 5.45 1.05

17.51 9.15 3.74 0.60

104.57 54.30 21.18 5.41

48.73 26.71 11.62 2.64

100

2
2 4
4 8
8 16
16

93.06 75.92 43.10 14.90

4.38 1.10 0.17 0.29

90.27 67.07 25.91 0.87

27.30 5.75 1.55 0.82

196.97 173.05 102.11 36.82

15.50 8.52 3.49 0.47

386.42 384.74 291.33 158.89

53.68 29.50 12.86 2.99

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

157

Fig. 5(b). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Normal stresses rxx computed through the thickness of the laminate without spurious terms.

Fig. 5(e). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Normal stresses ryy computed through the thickness of the laminate without spurious terms.

Fig. 5(c). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Convergence of normal stresses rxx solutions with and without spurious terms.

Fig. 5(f). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Convergence of normal stresses ryy solutions with and without spurious terms.

Fig. 5(d). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Normal stresses ryy computed through the thickness of the laminate with spurious terms.

Fig. 5(g). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Transverse shear stresses sxz computed through the thickness of the laminate with spurious terms.

158

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

Fig. 5(h). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Transverse shear stresses sxz computed through the thickness of the laminate without spurious terms.

Fig. 5(k). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Transverse shear stresses syz computed through the thickness of the laminate without spurious terms.

Fig. 5(i). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Convergence of transverse shear stresses sxz solutions with and without spurious terms.

Fig. 5(l). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Convergence of transverse shear stresses syz solutions with and without spurious terms.

Fig. 5(j). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Symmetric crossply laminate 0°/90°/0°. Transverse shear stresses syz computed through the thickness of the laminate with spurious terms.

Fig. 6(a). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Normal stresses rxx computed through the thickness of the laminate with spurious terms.

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

159

Fig. 6(b). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Normal stresses rxx computed through the thickness of the laminate without spurious terms.

Fig. 6(e). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Normal stresses ryy computed through the thickness of the laminate without spurious terms.

Fig. 6(c). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Convergence of normal stresses rxx solutions with and without spurious terms.

Fig. 6(f). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Convergence of normal stresses ryy solutions with and without spurious terms.

Fig. 6(d). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Normal stresses ryy computed through the thickness of the laminate with spurious terms.

Fig. 6(g). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Transverse shear stresses sxz computed through the thickness of the laminate with spurious terms.

160

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

Fig. 6(h). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Transverse shear stresses sxz computed through the thickness of the laminate without spurious terms.

Fig. 6(k). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Transverse shear stresses syz computed through the thickness of the laminate without spurious terms.

Fig. 6(i). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Convergence of transverse shear stresses sxz solutions with and without spurious terms.

Fig. 6(l). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Convergence of transverse shear stresses syz solutions with and without spurious terms.

Fig. 6(j). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Symmetric crossply laminate 0°/90°/0°. Transverse shear stresses syz computed through the thickness of the laminate with spurious terms.

presents results which are still far from the analytical results. Errors in Table 1 range from 93.06% to 14.90%. On the other hand, Fig. 6(b) shows that after the removal of the spurious terms, convergence occurs very quickly and well to the analytical solutions. The highest percent error is only 4.38% and drops to 0.29% with refinement. The convergence plots in Fig. 6(c) help to demonstrate these solutions behaviors. Figs. 6(d) and 6(e) show that both models with and without the spurious terms converge to the analytical solution for ryy. The finer meshes errors registered in Table 1 are only 0.87% and 0.82%, respectively. However, convergence occurs rather quickly using the model corrected for the spurious terms. The coarser mesh error of the model with spurious terms is 90.27% while the error of the model without spurious terms is only 27.30%. These observations are confirmed by the convergence plots shown in Fig. 6(f). Figs. 6(g)–6(l) show the solutions for the transverse shear stresses sxz and syz. The spurious terms cause a strong delay in convergence as demonstrated in Fig. 6(g) and in Fig. 6(j). The coarser meshes errors are 386.42% and 196.97% for both stresses, respectively, dropping to 158.89% and 36.82%, respectively. It is also

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

observed that coarser meshes produce results that are qualitatively erroneous (observe negative values). After the removal of the spurious terms, as shown in Fig. 6(h) and in Fig. 6(k), solutions for both transverse shear stresses converge rather well to the analytical solutions. Errors for syz range from 53.68% to 2.99%, while errors for sxz range from 15.50% to 0.47%. The convergence plots in Fig. 6(i) and in Fig. 6(l) reinforce these observations. Finally, solutions for a thinner plate (a/h = 1000) are presented in Figs. 7(a)–7(o), and percentage errors are listed in Table 2. These set of results also include those for the in-plane shear stresses sxy, which have not been evaluated for the previous plates. The behavior of these solutions are very much alike those of the plate with a/h = 100, so the authors understand that a thorough discussion would be very repetitive. However, the reader must observe, and this can be readily seen by the error values in Table 2, that shear locking effects are stronger in this thinner plate, as it should be expected from theoretical arguments. See that the error values are very high and do not seem to be attenuated by mesh refinement for stresses rxx, ryy, and sxy. Further, results are even worse for the transverse shear stresses sxz and syz as the errors increase with mesh refinement. Thus, the complete removal of the spurious terms without introduction of spurious zero-energy modes that is allowed by strain gradient notation is even more paramount in this situation.

161

Fig. 7(b). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Normal stresses rxx computed through the thickness of the laminate without spurious terms.

5.3. Application problem #3 The third application problem is a square simply supported (SS2 supports) anti-symmetric angle-ply laminated composite plate subjected to a uniform load of value qo = 10 N/m2. The laminate is composed of eight laminae of equal thicknesses and with lamination scheme [45°/+45°]4, and its sides are 1.0 m in length (a = b = 1.0 m) as shown in Fig. 8. The material properties of the previous problem are also employed here for each layer. Also, the same edge length-to-thickness ratios and finite element meshes are employed. Stress results are averaged nodal values and they are nondimensionalized according to Eq. (35). In-plane normal stresses rxx are calculated at the center point of the plate, and in-plane shear stresses sxy are calculated at a corner node, while transverse shear stresses sxz are calculated at a border middle point. Table 3 contains the highest percent error of each solution. First, the thick plate (a/h = 10) is analyzed. Figs. 9(a)–9(i) show the through-the-thickness stresses solutions. Figs. 9(a) and 9(b)

Fig. 7(c). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Convergence of normal stresses rxx solutions with and without spurious terms.

Fig. 7(a). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Normal stresses rxx computed through the thickness of the laminate with spurious terms.

Fig. 7(d). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Normal stresses ryy computed through the thickness of the laminate with spurious terms.

162

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

Fig. 7(e). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Normal stresses ryy computed through the thickness of the laminate without spurious terms.

Fig. 7(h). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Transverse shear stresses sxy computed through the thickness of the laminate without spurious terms.

Fig. 7(f). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Convergence of normal stresses ryy solutions with and without spurious terms.

Fig. 7(i). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Convergence of transverse shear stresses sxy solutions with and without spurious terms.

Fig. 7(g). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Transverse shear stresses sxy computed through the thickness of the laminate with spurious terms.

Fig. 7(j). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Transverse shear stresses sxz computed through the thickness of the laminate with spurious terms.

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

163

Fig. 7(k). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Transverse shear stresses sxz computed through the thickness of the laminate without spurious terms.

Fig. 7(n). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Transverse shear stresses syz computed through the thickness of the laminate without spurious terms.

Fig. 7(l). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Convergence of transverse shear stresses sxz solutions with and without spurious terms.

Fig. 7(o). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Convergence of transverse shear stresses syz solutions with and without spurious terms.

Fig. 7(m). Square simply-supported (SS1 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Symmetric cross-ply laminate 0°/90°/0°. Transverse shear stresses syz computed through the thickness of the laminate with spurious terms.

show the normal stresses rxx solutions. It is observed that shear locking effects are strong in the coarser mesh solution containing the spurious terms. The typical discontinuity between laminae is not depicted. This qualitative error is removed in further refinements. Refinement also provides convergence to the analytical solution of both models. Although the percent error is generally larger for the solutions provided by the model without spurious terms as registered in Table 3, it tends to a small value for both models (1.89% and 3.12%). Fig. 9(c) shows that both models converge to the analytical solution with mesh refinement. Figs. 9(d)–9(f) show the transverse shear stresses sxz solutions. It can be seen that solutions provided by the model without the spurious terms are always more accurate, although the solutions from the model with spurious terms also converge to the analytical solution. Percent error in Table 3 ranges from 83.32% to 4.56% for the model with the spurious terms, and from 39.58% to 3.03% for the model without the spurious terms. Therefore, refinement attenuates the shear locking effects. Further, according to Fig. 9 (f), solutions converge asymptotically from lower bounds. Figs. 9(g)–9(i) show the in-plane shear stresses sxy solutions. The finer meshes solutions provided by the model without the

164

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

Table 2 Percent errors in nondimensionalized stress components of symmetric cross-ply square plate (0°/90°/0°). Side-to-thickness relation a/h = 1000. Errors (%)  xx r

(0°/90°/0°)

sxy

 yy r

sxz

syz

a/h

Meshes

With/PS

Wout/PS

With/PS

Wout/PS

With/PS

Wout/PS

With/PS

Wout/PS

With/PS

Wout/PS

1000

2
2 4
4 8
8 16
16

99,93 99,68 98,73 95,10

4,35 1,10 0,17 0,29

99,89 99,56 98,25 93,23

27,43 5,74 1,55 0,82

99,94 99,73 98,92 95,79

36,85 14,76 4,91 1,11

11,56 28,14 36,21 32,41

15,48 8,49 3,48 0,46

209,69 268,63 302,74 310,86

53,72 29,53 12,89 3,01

Fig. 8. Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load. Anti-symmetric angle-ply laminate [45°/+45°]4. Problem definition.

Table 3 Percent errors in nondimensionalized stress components of antisymmetric angle-ply square plate (45°/+45°)4. Under uniform load. Errors (%)

sxz

 xx r

(45°/+45°)4

sxy

a/h

Meshes

With/PS

Wout/PS

With/PS

Wout/PS

With/PS

Wout/PS

10

2
2 4
4 8
8 16
16

36.84 3.01 2.96 1.89

29.99 27.18 8.48 3.12

83.32 38.86 16.65 4.56

39.58 21.89 10.88 3.03

51.73 8.32 10.76 18.11

65.81 20.99 0.41 8.49

100

2
2 4
4 8
8 16
16

97.52 86.93 61.09 26.82

26.89 24.79 6.27 2.46

217.20 237.08 196.98 101.84

38.59 19.59 11.34 3.48

98.18 89.10 62.77 22.46

67.96 26.26 4.49 5.80

1000

2
2 4
4 8
8 16
16

97.63 87.52 62.83 30.09

21.21 19.21 1.52 2.12

215.57 235.18 195.63 101.82

39.44 20.71 12.57 4.82

98.38 90.31 66.89 31.05

71.51 34.42 15.07 5.92

spurious terms are significantly better than the solutions provided by the model with the spurious terms. However, the error of 8.49% registered in Table 3 is still high, and convergence is not monotonic as the previous error is only 0.41%. Fig. 9(i) clearly shows that convergence is not monotonic for either model, and that coarse mesh solutions from the model with spurious terms are superior, which might be misleading. Next, a thin plate (a/h = 100) is analyzed. Figs. 10(a)–10(i) show the through-the-thickness stresses solutions. According to Fig. 10 (a), the spurious terms delay the convergence of the normal stresses rxx, indicating that more effort would be required to overcome the shear locking effects. In turn, Fig. 10(b) shows that the rate of convergence increases after the removal of the spurious terms. Percent error in Table 3 ranges from 97.52% to 26.82% for the model with the spurious terms, and from 26.89% to 2.46% for the model

without the spurious terms. Convergence plots in Fig. 10(c) show that the solutions without the effects of the spurious terms converge faster and in an asymptotic fashion from an upper bound. The spurious terms also prevent convergence of the transverse shear stresses sxz as shown in Fig. 10(d). However, convergence is allowed to occur after the removal of the spurious terms as shown in Fig. 10(e). Percent error in Table 3 ranges from 217.20% to 101.84% for the model with the spurious terms, and from 38.59% to 3.48% for the model without the spurious terms. The difference between the two sets of solutions is shown very clearly by the convergence plots in Fig. 10(f). The solutions without the spurious terms converge asymptotically from a lower bound. Fig. 10(g) shows that the spurious terms prevent convergence of the in-plane transverse shear stresses sxy whereas Fig. 10(h) shows that the solutions improve after the removal of the spurious terms.

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

165

Fig. 9(a). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Anti-symmetric angle-ply laminate [45°/+45°]4. Normal stresses rxx computed through the thickness of the laminate with spurious terms.

Fig. 9(d). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Anti-symmetric angle-ply laminate [45°/+45°]4. Transverse shear stresses sxz computed through the thickness of the laminate with spurious terms.

Fig. 9(b). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Anti-symmetric angle-ply laminate [45°/+45°]4. Normal stresses rxx computed through the thickness of the laminate without spurious terms.

Fig. 9(e). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Anti-symmetric angle-ply laminate [45°/+45°]4. Transverse shear stresses sxz computed through the thickness of the laminate without spurious terms.

Fig. 9(c). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Anti-symmetric angle-ply laminate [45°/+45°]4. Convergence of normal stresses rxx solutions with and without spurious terms.

Fig. 9(f). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Anti-symmetric angle-ply laminate [45°/+45°]4. Convergence of transverse shear stresses sxz solutions with and without spurious terms.

166

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

Fig. 9(g). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Anti-symmetric angle-ply laminate [45°/+45°]4. In-plane shear stresses sxy computed through the thickness of the laminate with spurious terms.

Fig. 10(a). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Anti-symmetric angle-ply laminate [45°/+45°]4. Normal stresses rxx computed through the thickness of the laminate with spurious terms.

Fig. 9(h). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Anti-symmetric angle-ply laminate [45°/+45°]4. In-plane shear stresses sxy computed through the thickness of the laminate without spurious terms.

Fig. 10(b). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Anti-symmetric angle-ply laminate [45°/+45°]4. Normal stresses rxx computed through the thickness of the laminate without spurious terms.

Fig. 9(i). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 10. Anti-symmetric angle-ply laminate [45°/+45°]4. Convergence of In-plane shear stresses sxy solutions with and without spurious terms.

Fig. 10(c). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Anti-symmetric angle-ply laminate [45°/+45°]4. Convergence of normal stresses rxx solutions with and without spurious terms.

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

167

Fig. 10(d). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Anti-symmetric angle-ply laminate [45°/+45°]4. Transverse shear stresses sxz computed through the thickness of the laminate with spurious terms.

Fig. 10(g). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Anti-symmetric angle-ply laminate [45°/+45°]4. In-plane shear stresses sxy computed through the thickness of the laminate with spurious terms.

Fig. 10(e). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Anti-symmetric angle-ply laminate [45°/+45°]4. Transverse shear stresses sxz computed through the thickness of the laminate without spurious terms.

Fig. 10(h). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Anti-symmetric angle-ply laminate [45°/+45°]4. In-plane shear stresses sxy computed through the thickness of the laminate without spurious terms.

Fig. 10(f). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Anti-symmetric angle-ply laminate [45°/+45°]4. Convergence of transverse shear stresses sxz solutions with and without spurious terms.

Fig. 10(i). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 100. Anti-symmetric angle-ply laminate [45°/+45°]4. Convergence of In-plane shear stresses sxy solutions with and without spurious terms.

168

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

Fig. 11(a). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Anti-symmetric angle-ply laminate [45°/+45°]4. Normal stresses rxx computed through the thickness of the laminate with spurious terms.

Fig. 11(b). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Anti-symmetric angle-ply laminate [45°/+45°]4. Normal stresses rxx computed through the thickness of the laminate without spurious terms.

Fig. 11(c). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Anti-symmetric angle-ply laminate [45°/+45°]4. Convergence of normal stresses rxx solutions with and without spurious terms.

Fig. 11(d). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Anti-symmetric angle-ply laminate [45°/+45°]4. Transverse shear stresses sxz computed through the thickness of the laminate with spurious terms.

Fig. 11(e). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Anti-symmetric angle-ply laminate [45°/+45°]4. Transverse shear stresses sxz computed through the thickness of the laminate without spurious terms.

Fig. 11(f). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Anti-symmetric angle-ply laminate [45°/+45°]4. Convergence of transverse shear stresses sxz solutions with and without spurious terms.

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

169

Percent error in Table 3 ranges from 98.18% to 22.46% for the model with the spurious terms, and from 67.96% to 5.80% for the model without the spurious terms. It is observed, however, that the finer mesh results at the lower laminae are larger than the respective analytical results. Convergence plots in Fig. 10(i) show the superiority of the solutions without spurious terms, and also that they do not converge asymptotically. Analysis of a thinner plate (a/h = 1000) is also performed. Stresses results are shown in Figs. 11(a)–11(i). As in Application #2, this is a very thin plate and the effects of shear locking are to be more pronounced. In fact, solutions for rxx and sxy show to converge very slowly due to the presence of the spurious terms. However, more critical are the solutions for the transverse shear stresses sxz where percent errors range from 215% to 101% when spurious terms are present. The corrected model presents solutions which converge monotonically to the analytical solution within reasonably good accuracy. Fig. 11(g). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Anti-symmetric angle-ply laminate [45°/+45°]4. In-plane shear stresses sxy computed through the thickness of the laminate with spurious terms.

5.4. Application problem #4 The final example problem consists of a square, simplysupported, four-laminae symmetric cross-ply laminated plate with lamination scheme [0°/90°/90°/0°] which is subjected to a distributed sinusoidal load, as shown in Fig. 12. Such loading is given by the following expression:

qðx; yÞ ¼ qo sin

Fig. 11(h). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Anti-symmetric angle-ply laminate [45°/+45°]4. In-plane shear stresses sxy computed through the thickness of the laminate without spurious terms.

Fig. 11(i). Square simply-supported (SS2 boundary conditions) plate subjected to a uniformly distributed load with side-to-thickness ratio a/h = 1000. Anti-symmetric angle-ply laminate [45°/+45°]4. Convergence of In-plane shear stresses sxy solutions with and without spurious terms.

px a

sin

py b

ð36Þ

where q0 = 10 N/m2, and a and b are the edged lengths of the plate (a = b = 1.0 m). The material properties for each lamina are the same used in the previous two applications. This example presents a comparison of the corrected strain gradient serendipity element with the isoparametric element in three versions, namely; (i) with complete or full numerical integration (C), (ii) with uniform reduced integration (R), and (iii) with selective reduced integration (S). Results for these isoparametric elements are in Reddy [2], and they are only available for uniform meshes composed of 16 elements (4
4) and 64 elements (8
8). Also, our discussion is limited to the plate with aspect ratio a/h = 100. Although limited, these results show a comparison in performances of the strain gradient and isoparametric elements. Results of transverse displacement and stress components are presented in Table 4. Plots of stress components distribution are not usefull in this case because solutions are very close to each other in regular scale. In Table 4, the models employed are generally represented by nQ8-e where n is either 4 or 8 to indicate a 4
4 or an 8
8 mesh, Q8 indicates that the element is an eight-node quadrilateral, and e can be C, R, or S for isoparametric element with complete, reduced or selective integration, respectively, or yet e can be SG for strain gradient element. Analytical (Navier) solution is included for percent error estimates, which are enclosed in parentheses below the computed quantities. The isoparametric elements provided better results than the strain gradient element when the coarser mesh (4
4) was employed. Errors of the 4Q8-SG are significantly larger than others’. However, such discrepancy is reduced according to the results of the next refinement. Comparison is made only between the isoparametric element with selective reduced integration 8Q8-S and the strain gradient element 8Q8-SG. Errors associated to both models are within 1%, which can be considered low error level. In view of this result, the conclusion is that both models are comparable in accuracy. However, results provided by the strain gradient model are higher than the analytical solution values, which is the reason why errors have been indicated with negative sign. Thus, two further analyses using the strain gradient model have been performed, namely; 16Q8-SG and 32Q8-SG, with the purpose of

170

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

Fig. 12. Square simply-supported (SS2 boundary conditions) plate subjected to a sinusoidally distributed load with side-to-thickness ratio a/h = 100. Symmetric cross-ply laminate 0°/90°/90°/0°. Problem definition.

Table 4 Results in nondimensionalized stress components of symmetric cross-ply square plate (0°/90°/90°/0°) under sinusoidal load. Side-to-thickness relation a/h = 100. Comparison of strain gradient element and isoparametric elements.

[0°/90°/90°/0°] a/h = 100

Model

W
102

 xx r

 yy r

sxz

syz

4Q8-C

0.4143 (4.470) 0.4319 (0.415) 0.4319 (0.415) 0.4080 (5.920)

0.4900 (8.950) 0.5214 (3.122) 0.5214 (3.122) 0.5000 (7.098)

0.2435 (9.940) 0.2621 (3.070) 0.2620 (3.106) 0.2512 (3.122)

0.4300 (3.370) 0.4350 (2.247) 0.4330 (2.697) 0.4070 (8.539)

0.1000 (0.990) 0.1020 (0.990) 0.1020 (0.990) 0.0920 (8.910)

0.4336 (0.023) 0.4372 (0.800) 0.4372 (0.800) 0.4370 (0.761) 0.4337

0.5344 (0.706) 0.5386 (0.074) 0.5385 (0.056) 0.5384 (0.037) 0.5382

0.2685 (0.703) 0.2707 (0.110) 0.2706 (0.074) 0.2705 (0.037) 0.2704

0.4410 (0.899) 0.4450 (0) 0.4452 (0.045) 0.4445 (0.112) 0.4450

0.1000 (0.990) 0.1020 (0.990) 0.1025 (1.49) 0.1017 (0.693) 0.1010

4Q8-R 4Q8-S 4Q8-SG

8Q8-S 8Q8-SG 16Q8-SG 32Q8-SG Navier solution

observing the convergence characteristics of the element in this particular problem. Results, which are also included in Table 4, indicate that there is not monotonic convergence (oscillations occur), but there is very good accuracy.

6. Conclusion The formulation of an eight-node serendipity element using strain gradient notation for the analysis of laminated composite plates has been presented. Inspection of the physically interpretable shear strain polynomial expansions of the element has revealed the presence of various spurious terms which are responsible for shear locking. These spurious terms have been clearly identified as normal strain gradient terms which unduly appear in the shear strain expansions. Once identified, the spurious terms

have been removed from those strain expansions, thus eliminating a-priori the potential for shear locking manifestation. It is important to note that those terms are only spurious in the shear strain expressions. They are legitimate terms in the normal strain expansions, from which they cannot be removed, and they are legitimate part of the displacement expansions, forming the kinematic basis of the element. In addition, compatible mode terms, which might be mistaken as spurious, have been clearly identified by the recognition that they comprise compatibility equations. Those terms have been retained in order to avoid the introduction of spurious zero-energy modes (or rank deficiency). Thick and thin rectangular plates have been analyzed for different lamination schemes and boundary conditions. In each case, stress solutions provided by the model containing the spurious terms have been compared to stress solutions provided by the corrected model. Those numerical analyses have demonstrated that

J.E. Abdalla Filho et al. / Composite Structures 154 (2016) 150–171

the identified spurious terms are the cause of shear locking. Those analyses have also demonstrated the effectiveness of the procedure employed to eliminate the spurious terms as, in general, solutions provided by the corrected model converged monotonically and faster to the analytical solutions. Further, numerical solutions have shown that mostly transverse shear stress solutions provided by the model containing the spurious terms might not converge adequately to the correct solutions as the results oscillate or contain great levels of error such as those in Tables 1–3. Last but not least, the proposed model has provided results within very good accuracy so that it can be considered an efficient element for analyzing laminated composite plate structures. The use of strain gradient notation may be viewed as advantageous as it allows spurious terms responsible for shear locking to be identified and eliminated a-priori. The authors claim to have demonstrated theoretically the sources of shear locking and spurious zero-energy modes in the eight-node serendipity plate element, and to have built an efficient FSDT element to analyze laminated composite plates. As further steps, having considered the literature review performed, these authors have envisioned (i) the development of triangular elements using the same procedure adopted in this work; (ii) the association of the refined zigzag theory (RZT) to strain gradient notation plate elements; (iii) the development of strain gradient notation plate elements for nonlinear analysis of laminated composites; and, (iv) further, an attempt to associate strain gradient notation to a higher-order theory for laminated composites might provide a fruitful result. Acknowledgements The first author acknowledges CAPES (a Foundation affiliated with the Ministry of Education of Brazil) for the scholarship granted during his stay at the University of Colorado Boulder in 2015/16 as a Visiting Research Scholar, where this work was finalized. References [1] Reddy JN, Averill RC. Advances in the modeling of laminated plates. Comput Syst Eng 1991;2(5–6):541–55. [2] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. 2nd ed. Boca Raton, FL: CRC Press; 2004. [3] Ghugal YM, Shimpi RP. A review of refined shear deformation theories of isotropic and anisotropic laminated plates. J Reinf Plast Compos 2002;21 (9):775–805. [4] Reddy JN, Robbins Jr DH. Theories and computational models for composite laminates. Appl Mech Rev 1994;47(6):147–65. [5] Carrera E. C0z requirements – models for the two dimensional analysis of multilayered structures. Compos Struct 1997;37(3–4):373–83. [6] Kant T, Swaminathan K. Estimation of transverse/interlaminar stresses in laminated composites – a selective review and survey of current developments. Compos Struct 2000;49:65–75. [7] Mawenya AS, Davies JD. Finite element bending analysis of multilayer plates. Int J Numer Methods Eng 1974;8:215–25. [8] Panda SC, Natarajan R. Finite element analysis of laminated composite plates. Int J Numer Methods Eng 1979;14:69–79. [9] Ahmad S, Irons BM, Zienkiewicz OC. Analysis of thick and thin shell structures by curved finite elements. Int J Numer Methods Eng 1970;2:419–51. [10] Lo KH, Cristensen RM, Wu EM. A high-order theory of plate deformation, part 1: homogeneous plates. J Appl Mech 1977;44(4):663–8. [11] Lo KH, Cristensen RM, Wu EM. A high-order theory of plate deformation, part 2: laminated plates. J Appl Mech 1977;44(4):669–76. [12] Singh G, Rao GV. A discussion on simple third-order theories and elasticity approaches for flexure of laminated plates. Struct Eng Mech 1995;3 (2):121–33. [13] Bose P, Reddy JN. Analysis of composite plates using various plate theories, part 1: formulation and analytical solutions. Struct Eng Mech 1998;6 (6):583–612.

171

[14] Bose P, Reddy JN. Analysis of composite plates using various plate theories, part 2: finite element model and numerical results. Struct Eng Mech 1998;6 (7):727–46. [15] Reddy JN, Wang CM. An overview of the relationships between solutions of the classical and shear deformation plate theories. Compos Sci Technol 2000;60:2327–35. [16] Reddy JN. On refined computational models of composite laminates. Int J Numer Methods Eng 1989;27:361–82. [17] Zienkiewicz OC, Taylor RL, Too JM. Reduced integration technique in general analysis of plates and shells. Int J Numer Methods Eng 1971;3:275–90. [18] Hughes TJR, Taylor RL, Kanoknukulchai W. A simple and efficient finite element for plate bending. Int J Numer Methods Eng 1977;11:1529–43. [19] Hughes TJR, Cohen M, Haroun M. Reduced and selective integration techniques in finite element analysis of plates. Nucl Eng Des 1978;46:203–22. [20] Hughes TJR, Tezduyar TE. Finite elements based upon Mindlin plate theory with particular reference to the four-node bilinear isoparametric element. J Appl Mech 1981;48:587–96. [21] Hughes TJR, Cohen M. The ‘‘heterosis” finite element for plate bending. Comput Struct 1978;9:445–50. [22] Prathap G. A field-consistency approach to plate elements. Struct Eng Mech 1997;5(6):853–65. [23] Bathe KJ, Dvorkin E. A four-node plate bending element based on Mindlin/ Reissner plate theory and a mixed interpolation. Int J Numer Methods Eng 1985;21:367–83. [24] Brank B, Carrera E. A family of shear-deformable shell finite elements for composite structures. Comput Struct 2000;76:287–97. [25] Verwoerd MH, Kok AWM. A shear locking free six-node Mindlin plate bending element. Comput Struct 1990;36(3):547–51. [26] Kabir HR. A shear-locking free robust isoparametric three-node triangular finite element for moderately-thick and thin arbitrarily laminated plates. Comput Struct 1995;57(4):589–97. [27] Zhang YX, Yang CH. A family of simple and robust finite elements for linear and geometrically nonlinear analysis of laminated composite plates. Compos Struct 2006;75:545–52. [28] Dau F, Polit O, Touratier M. C1 plate and shell finite elements for geometrically nonlinear analysis of multilayered structures. Comput Struct 2006;84:1264–74. [29] Kim KD, Lee CS, Han SC. A 4-node co-rotational ANS shell element for laminated composite structures. Compos Struct 2007;80:234–52. [30] Botello S, Onate E, Canet JM. A layer-wise triangle for analysis of laminated composite plates and shells. Comput Struct 1999;70:635–46. [31] Sheikh AH, Chakrabarti A. A new plate bending element based on higher-order shear deformation theory for the analysis of composite plates. Finite Elem Anal Des 2003;39:883–903. [32] Aagaah MR, Mahinfalah M, Jazar GN. Linear static analysis and finite element modeling for laminated composite plates using third order shear deformation theory. Compos Struct 2003;62:27–39. [33] Moleiro F, Mota Soares CM, Mota Soares CA, Reddy JN. Layerwise mixed leastsquares finite element models for static and free vibration analysis of multilayered composite plates. Compos Struct 2010;92:2328–38. [34] Mantari JL, Guedes SoaresC. Generalized layerwise HSDT and finite element formulation for symmetric laminated and sandwich composite plates. Compos Struct 2013;105:319–31. [35] Pandey S, Prandyumna S. A new C0 higher-order layerwise finite element formulation for the analysis of laminated and sandwich plates. Compos Struct 2015;131:1–16. [36] Tessler A, Di Sciuva M, Gherlone M. A consistent refinement of first-order shear deformation theory for laminated composite and sandwich plates using improved zigzag kinematics. J Mech Mater Struct 2010;5(2):341–67. [37] Gherlone M, Tessler A, Di Sciuva M. C0 beam elements based on the refined zigzag theory for multilayered composite and sandwich laminates. Compos Struct 2011;93(11):2882–94. [38] Versino D, Gherlone M, Matone M, Di Sciuva M, Tessler A. C0 triangular elements based on the refined zigzag theory for multilayer composite and sandwich plates. Compos B Eng 2013;44(1):218–30. [39] Eijo A, Onate E, Oller S. A four-noded quadrilateral element for composite laminated plates/shells using the refined zigzag theory. Int J Numer Methods Eng 2013;95(8):631–60. [40] Natarajan S, Ferreira AJM, Bordas SPA, Carrera E, Cinefra M. Analysis of composite plates by a unified formulation – cell based smoothed finite element method and field consistent elements. Compos Struct 2013;105:75–81. [41] Dow JO. A unified approach to the finite element method and error analysis procedures. San Diego, CA: Academic Press; 1999. [42] Abdalla Filho JE, Belo IM, Pereira MS. A laminated composite plate finite element a-priori corrected for locking. Struct Eng Mech 2008;28(5):603–33. [43] Abdalla Filho JE. Qualitative and discretization error analysis of laminated composite plate models [Doctoral thesis]. University of Colorado Boulder; 1992.