A new trigonometric zigzag theory for buckling and free vibration

A new trigonometric zigzag theory for buckling and free vibration analysis of laminated composite and sandwichplatesRosalin Sahoo,B.N.Singh ⇑Department of Aerospace Engineering,Indian Institute of Technology Kharagpur,Kharagpur 721302,Indiaa r t i c l e i n f o Article history:Available online 28June 2014Keywords:Zigzag theory Free vibration BucklingFinite element method Laminated composite Sandwich platesa b s t r a c tIn this work,a new trigonometric zigzag theory is proposed for the analysis of laminated composite and sandwich plates.The theory is based upon shear strain shape function assuming non-linear distribution of transverse shear stresses.It satisfies the necessary conditions of inter-laminar stress continuity at the layer interfaces as well as the condition of zero transverse shear stresses at the top and bottom surfaces of the plates.The theory has the same number of unknown field variables as that of FSDT and the number of unknowns is layer independent which makes the solution computationally more efficient.A C 0continu-ous isoparametric serendipity element is employed to solve the discrete eigenvalue equations arising in both the problems.The accuracy and the efficiency of the theory are thus demonstrated through numer-ical experiments on the free vibration and stability analysis of laminated composite and sandwich plates with different modular ratio,aspect ratio,span to thickness ratio,loading and boundary conditions,ply orientations,lay-up number,eigen modes,etc.Higher modes with corresponding mode shapes of vibra-tion and buckling are also represented for laminated composite and sandwich plates.An elegant perfor-mance of the proposed theory is observed when it is validated with relevant numerical examples.Ó2014Elsevier Ltd.All rights reserved.1.IntroductionThe encompassing growth of laminated composite plates in var-ious weight sensitive industries such as aerospace,marine,naval and civil engineering constructions is due to its high strength/stiff-ness to weight ratio,better energy absorption capability and flex-ibility in design over the conventional materials.However,these structures are weak in shear due to their low shear modulus compared to extensional rigidity.Hence,the effect of shear defor-mation is quite significant and it becomes more complex in case of sandwich construction due to the intricate variation of their through-thickness properties.The exact three dimensional (3D)elasticity solutions have the ability to capture its structural response adequately [1–5].However,such elasticity solutions are restricted to certain geometries,loading and boundary conditions.The computational efforts for 3D finite element solutions are not adequate and often intractable.Furthermore,although the com-posite materials carry many desirable properties,they present challenging problems in the understanding of their structural behavior.Hence,the studies involving the stability behaviors and natural vibrations of those structures are of great importance and essential in predicting their structural responses.Thus,the devel-opment of an appropriate mathematical two dimensional (2D)model to represent the behavior of laminated composite and sand-wich plates has drawn a considerable amount of attention.An extensive review of the literature on the development of modern plate theories have recently been carried out [6–11].The shear deformation plate theories can be broadly classified into displacement-based and stress-based theories.Basset’s theory [12]appears to have been the first to pioneer the displacement based shear deformation theory which inform that the displace-ments can be expanded in power series of the thickness coordi-nate.Following Basset’s work,Reissner [13]and Mindlin [14]developed first order shear deformation theory (FSDT)which is also called Reissner–Mindlin’s plate theory.This theory over-predicts the free vibration and buckling load since it assumes uni-form transverse shear strain over the entire plate thickness [15].The FSDT requires a shear correction factor to compensate for the actual parabolic variation of the shear stress which fluctuates with the geometric parameters,loading and boundary conditions [16].The displacement-based higher order theories (HSDT),where the components of displacements are expanded in power series of thickness coordinate were later formulated [17–19]as further improvement.However,the increase in the order of variation for the in-plane displacements across the thickness is unable to make/10.1016/pstruct.2014.05.0020263-8223/Ó2014Elsevier Ltd.All rights reserved.⇑Corresponding author.Tel.:+913222283026;fax:+913222255303.E-mail address:bnsingh@aero.iitkgp.ernet.in (B.N.Singh).any significant improvement in representing the behavior of mul-tilayered laminates.The non-polynomial displacement based higher order theories[20–25]were suggested using different shear strain functions i.e.,trigonometric,hyperbolic and exponential functions.Touratier[21]proposed a trigonometric higher order shear deformation theory which includes sinus function originally developed by Levy[22]and improved by Stein[23].Soldatos[24] used hyperbolic function of sinus for the development of HSDT.A new exponential function proposed as a shear stress function for laminated composite beam by Karama et al.[25]which was later extended to laminated plates[26].The theory proposed by Karama et al.[26]is modified by Aydogdu[27]with the introduction of a transverse shear stress parameter.A new hyperbolic shear defor-mation theory was proposed by Meiche et al.[28]for buckling and vibration analysis of functionally graded sandwich plates. Mantari et al.[29]also proposed new shear strain functions such as tangent for the analytical solution of laminated and sandwich plates.Grover et al.[30,31]proposed trigonometric higher order theories for the analytical and numerical analysis of laminated composite and sandwich plates considering different shear strain functions such as inverse hyperbolic sinus function,inverse cotan-gent function and secant function.The vibration and buckling anal-ysis of laminated composite and sandwich plates based on HSDT are jotted down[31–37].The all above discussed theories fall under the equivalent single layer(ESL)two-dimensional theories where the deformation of plate is expressed in terms of unknown parameters of the reference plane,i.e.middle plane.The ESL theo-ries predict accurate results for global responses such as funda-mental vibration frequencies,and critical buckling loads. However,the major deficiency of the ESL theories in modeling composite laminated plates is that the transverse strain compo-nents are continuous across layer interfaces due to different values of shear rigidity at the adjacent layers.However,the practical sce-nario is just the reverse,i.e.,the shear strain is discontinuous and the shear stress is continuous at the layer interfaces.To overcome these disparities,layer-wise(LW)laminate theories are developed that provides a framework upon which the unknown displacement components are taken at all the layer interfaces,that give a zigzag pattern of through-thickness variation for the in-plane displace-ment to represent the desired shear strain discontinuity at the layer interfaces.Lu and Liu[38],and Robbins and Reddy[39]have considered different displacementfield in each layer maintaining displacement continuity across the layers.The trigonometric sinus function[21]was employed in LWT proposed by Arya et al.[40], Shimpi and Ghugal[41],Shimpi and Aynapure[42]for the analysis of composite beams.Mantari[43,44]approached trigonometric LWT for the analysis of composite plates.Thefinite element anal-ysis of Vidal and Polit[45],Gaudenzi et al.[46],Pai and Palazotto [47]and the work of Roque et al.[48]using mesh less method is quite significant.Ferreira et al.[49]used radial basis collocation function to analyze static,buckling and vibration responses of the plate.The framework of stress based shear deformation plate theories was developed by Reissner[50].The theory assumes a linear distri-bution of the in-plane normal and shear stresses through the thick-ness of the plate.Toledano and Murakami[51]and Carrera[52] analyzed mixed layer wise models.Yuan and Dawe[53]developed splinefinite strip method for predicting the natural frequencies of vibration and the buckling stresses of rectangular sandwich panels. Dafedar et al.[54]predicted an overall and local buckling loads sandwich plates.Shufrin et al.[55]studied a semi analytical method for the buckling analysis of laminated rectangular plates. Rao and Desai[56]and Rao et al.[57]presented a semi-analytical solution to evaluate the natural frequency based on LW theory. Although LW theories are quite good in predicting the natural frequencies and buckling load of laminated composite plates,it requires a huge computational effort as the numbers of unknowns involved in the formulation are layer dependent.The drawbacks of LW theories are counteracted by the intro-duction of zigzag(ZZ)theory which makes the number of variables independent of the layers by satisfying the condition of transverse shear stress continuity at the layer interfaces.Some improvement on this theory has been done by Di Sciuva[58],Bhaskar and Vara-dan[59],Lee and Liu[60],and Cho and Parmerter[61]by satisfying the condition of zero transverse shear stresses at the top and bot-tom surface of the plates.The ZZ theories provide parabolic through thickness variation of transverse shear strains and zero shear strains at the top and bottom surfaces of the laminate.A comprehensive review on the zigzag theories is presented by Carrera[62]for the analysis of multilayered laminated plates and shells,in which there is a discussion about three basic theories i.e.,Lekhnitskii multilayered theory,Ambartsumian multilayered theory and Reissner multilayered theory.Moreover,Carrera[63]and Demasi[64]considered higher order terms in the displacementfield,using Murakami’s ZZ func-tion[65].Wu et al.[66],Rodrigues et al.[67]and Neves et al.[68]also developed an approximate solution based on Murkami’s ZZ theory for the analysis of the laminated composite plates. Owing to the restrictions of closed form analytical solution in the orientation of geometry,loading and boundaries,thefinite element method(FEM)is approached to be the most significant one among the available computational methods in modeling of structures. Hence,Di Sciuva[69],Cho and Parmerter[70],Chakrabarti and Sheikh[71],Kulkarni and Kapuria[72],Pandit et al.[73–75],Singh et al.[76],Chalak et al.[77,78],Khandelwal et al.[79]and Sahoo and Singh[80,81]analyzed thefinite element representation of laminated composites based on the ZZ theory.It is observed from the discussed literature reviews in the above paragraphs that the developed zigzag theories with power series expansion predicts quite accurate results at the cost of high com-putational effort as more number of unknownfield variables are considered.Hence,there is a genuine need of development of a computationally efficient plate model that satisfies all the required criteria of ZZ theory with minimum number of variables.It can also be pointed out from the past reviews that there is no FEM result regarding the prediction of free vibration and buckling load param-eters of laminated composite and sandwich plates are published yet based on the theory that unites the zigzag concept with trigo-nometric shear strain functions.The above observations have inspired the author to develop a new displacement based zigzag plate theory where the trigonometric function specifically the secant function is assumed as the shear strain shape function which gives non-linear distribution of transverse shear stresses and zero transverse normal strain.Furthermore,the proposed zigzag theory with secant function(ZZTSF)fulfills not only the traction free boundary conditions on the top and bottom surfaces of the laminated structures but also satisfies the inter-laminar stress continuity at the layer interfaces obviating the need of shear correction factor.The ZZTSF requiresfive number offield variables as that of FSDT to complete the whole analysis which makes the theory computationally more sound compared to the available zig-zag theory that consist Taylor’s power series expansion.A compu-tationally efficient C0continuousfinite element is employed and an eight noded isoparametric biquadratic element is used to solve the eigenvalue equations for vibration and stability analysis. Numerical results for free vibration and buckling load parameters are evaluated for laminated composite and sandwich plates fulfill-ing the convergence criteria.The effect of different features such as boundary conditions,aspect and span to thickness ratio and load-ing conditions on the fundamental frequency and on buckling load parameters of laminated composite and sandwich structures are investigated.Higher modes of free vibration,buckling andR.Sahoo,B.N.Singh/Composite Structures117(2014)316–332317corresponding mode shapes are also presented for some numerical problem.The accuracy and precision of proposed model is mani-fested by comparing the evaluated results with the three-dimen-sional elasticity solutions and with the existing results based on other shear deformation theories as well.The effectualness of ZZTSF is ascertained by engrossing the minimum percentage difference among other deformation theories with respect to the 3D elasticity solution with minimum computational efforts.2.Mathematical modelingA rectangular laminated plate of dimensions(aÂbÂh)is com-posed offinite number of linearly elastic,orthotropic layers bonded together as given in Fig.1.A Cartesian co-ordinate system (x–y–z)is associated with the reference plane of the plate.2.1.Displacement modelThe chosen displacementfield is assumed to be the combina-tion of mid-plane displacements,rotations of the normal to the mid-plane varies with thickness,a linear zigzag and shear strain shape functions varying with layers as in Eq.(1).The out of plane transverse displacement is considered to be constant throughout the thickness.The unit step function is included to handle the zigzag concept.Fig.2indicates the through thickness variation of the in-plane displacements which is the combination of a linear zigzag model with different slopes and the parameter X(k)in each layer.However,the newly developed displacement model is repre-sented as follows:Uðx;y;z;tÞ¼u0Àz@w0þX n uÀ1i¼1zÀz uiÀÁH zÀz uiÀÁa ixuþX n lÀ1j¼1zÀz ljHÀzþz lja jxlþgðzÞþXðkÞzh ib xVðx;y;z;tÞ¼v0Àz@w0@yþX n uÀ1i¼1zÀz uiÀÁH zÀz uiÀÁa iyuþX n lÀ1j¼1zÀz ljHÀzþz lja jylþgðzÞþXðkÞzh ibyWðx;y;z;tÞ¼w0ð1ÞWhere u0,v0,w0,are the displacements of a point on the mid plane along x,y and z directions and b x,b y are the shear deformations at the mid plane about the y-and x-axis respectively.n u and n l are number of upper and lower layers respectively.a i xu;a i yu;a j xl;a j yu are the slopes of i th and j th layer corresponding to upper and lowerlayers respectively and H zÀz uiÀÁand HÀzþz ljare the Heavisidestep functions.Zero tangential stress-free boundary conditions are automatically satisfied at top and bottom surfaces of the plates (s x,z,s yz=0at z=±h/2)in the present kinematics.The considered value of the function g(z)is z sec(rz/h).However,for the top and bottom surfaces of the laminated plates the value of parameter X is assumed to beÀsec(r/2)[1+(r/2)tan(r/2)].The optimized value of transverse shear stress parameter is considered as0.1which is ascertained by the inverse method in the post processing step and the detail information is provided in Section4.The variation of parameter XðkÞ;a i xu;a i yu;a j xl;a j yl occurs with the change in layers and is evaluated by including the inter-laminar continuity condi-tions s i xz¼s iþ1xz and s i yz¼s iþ1yzand at the layer interfaces.Theshear strain shape functions of different models are presented in Table1.Fig.1.Schematic diagram of a laminate.Structures117(2014)316–3322.2.Constitutive equationsThe Hooke’s law relating the stress–strain behavior of k th layerof an orthotropic lamina with respect to structural axis system (x –y –z )may be presented asf r ijg 5Â1¼½Q ij ðk Þf e ij g 5Â1ð2Þwheref rij g 5Â1¼r xx r yy r xy r yz r xz ½ T;f e ij g 5Â1¼e xx e yy½e xy e yz e xz T and ½Q ij ðk Þis the reduced transformed rigidity matrix ofk th layer of multilayered plate.The rigidity matrix ½Q ij ðk Þcan be evaluated in terms of material properties and fibre orientation a .The shear strains c ij are considered as two times of corresponding mathematical strain components i.e.,c ij =2e ij .3.Finite element discretization 3.1.Continuity requirementAn efficient finite element model needs to fulfill the continuity requirement and to select a proper element.Presently,the actual displacement of ZZTSF requires C 1continuity at the element inter-faces during its finite element implementation as it carries first order derivatives of transverse displacement as shown in Eq.(1).In order to circumvent the C 1continuity requirement associated with the transverse displacement,the derivatives of transverse dis-placements are defined in terms of separate field variables.There-fore,five unknown variables with C 1continuity is transformed into seven number of field variables with C 0continuity.Hence,the modified form of displacement field for an arbitrary composite plate can be represented as:U ðx ;y ;z ;t Þ¼u 0Àz h x þX n u À1i ¼1z Àz u iÀÁH z Àz u i ÀÁa ixuþX n l À1j ¼1z Àz l j H Àz þz l j a j xl þg ðz ÞþX ðk Þz h i b xV ðx ;y ;z ;t Þ¼v 0Àz h y þX n u À1i ¼1z Àz u i ÀÁH z Àz ui ÀÁa iyuþX n l À1j ¼1z Àz l j H Àz þz l j a j yl þg ðz ÞþX ðk Þz h ib yW ðx ;y ;z ;t Þ¼w 0ð3a Þwhere h x ¼@w 0and h y ¼@w 0.The incorporation of various approaches such as concepts of independent interpolations,penalty functions and hermitian func-tions are used to approximate the transverse displacement in theformulations.However,in the present study,some artificial con-straints have been imposed in the preceding substitution for com-pensating the assumption,which are enforced variationaly through the penalty approach being in line with Ref.[82].@w 0@xÀh x ¼0;@w 0@yÀh y ¼0ð3b ÞThe displacement vector based on the modified displacement field(Eq.(3a))is rewritten asf dg ¼u 0v 0w 0b x b y h x h yÈÉTð3c Þ3.2.Selection of the elementAn eight noded isoparametric biquadratic quadrilateral seren-dipity element (Fig.3)with seven degrees of freedom (u 0,v 0,w 0,b x ,b y ,h x ,h y )per node is employed for the finite element approxi-mation of the ZZTSF.The element has an arbitrary rectangular geometry in x -y coordinate system.The element is mapped to n Àg plane to have a regular rectangular geometry.Consequently,for a discrete finite element domain,the field variables and ele-ment geometry can be expressed in terms of N i ,that is the interpo-lation function [83]associated with the node as follows:d ¼X n i ¼1N i d i ;x ¼Xn i ¼1N i x i ;y ¼X n i ¼1N i y ið4ÞWhere is the generalized field variable,d i is the value of field vari-able corresponding to node i ,x ,y are the generalized Cartesian coor-dinate,x i ,y i are the corresponding coordinate value at i th node and n is the number of nodes per element that is eight in the present study.3.3.Linear strain–displacement relationsThe linear strain–displacement relationships are stated asf eg ¼U ;xV ;y U ;y þV ;x V ;z þW ;y U ;z þW ;x ½Tð5ÞWhere,comma (,)denotes the partial differentiation with respect to the coordinate subscripts that follow.The linear strains given in Eq.(5)can be presented in terms of generalized strains and the generalized strains are the functions of generalized field variables.f eg ¼½H f e gð6ÞWhere the elements of [H ]includes the functions of z ,unit step functions and the parameter X as prescribed in Appendix A .Again,f eg ¼e 01e 02e06k 11k 12k 16k 21k 32k 46k 56e 04e05k 64k 75no Tð7ÞFurthermore,the generalized strains {e }that appeared in Eq.(6)may be expressed in terms of unknowns containing nodal field vari-ables asfollows:Table 1Shear strain shape functions of various shear deformation theories.Theory g (z )XReddy [19]z À4z 33h 0Touratier [21],Void and Polit [45]hp sin p z hÀÁ0Soldatos [24]h sin h z hcos h 12ÀÁ¼1:1276Karama et al.[25]ze À2ðz =h Þ20Aydogdu [27]z a À2ðz =h Þ2=ln a0Mieche et al.[28]0.6626hp sin h p z ÀÁÀ0.6626Mantari et al.[29]tan mzÀm sec 2mh 2Grover et al.[31]sin h À1rz ÀÁÀ2r hffiffiffiffiffiffiffiffir 2þ4p ;r ¼3:0Presentz sec (rz /h )Àsec (r /2)[1+(r /2)tan(r /2)];r =0.1R.Sahoo,B.N.Singh /Composite Structures 117(2014)316–332319f e g¼½B f d gð8ÞWhere[B]is the strain–displacement matrix in the Cartesian coor-dinate system and{d}is the displacement vector consists of the fun-damental variables for structural deformations.3.4.Geometric strain–displacement relationsThe buckling of a plate may be assumed to be a geometrical nonlinear problem for its accurate prediction.The generalized geo-metric strain vector for the present case may be illustrated asf e G g¼12W2;xþ12U2;xþ12V2;x1W2;yþ1U2;yþ1V2;yW;x W;yþU;x U;yþV;x V;y264375¼12½A G f h G gð9ÞWhere½A G ¼W;x0U;x0V;x00W;y0U;y0V;yW;y W;x U;y U;x V;y V;x264375ð10Þf h G g¼W;x W;y U;x U;y V;x V;y½ Tð11ÞAnd a comma(,)indicates the partial differentiation with respect to the coordinate subscripts that follow(as mentioned before).Again,using Eq.(8),the preceding equation may be stated as f h G g¼½H G f e g¼½H G ½B G f d gð12ÞWhere[H G]is also a function of z,unit step functions and the parameter X like that of[H]which are informed in Appendix A. 3.5.Energy calculationsThe total potential energy of the system may be expressed as P e¼TÀU SÀU inÀW extÀU peð13ÞWhere T is the total kinetic energy,U s is the strain energy due to deformation,U in is the potential energy due to applied in-plane compressive loads(r xx,r yy,r xy)is the work done by the external transverse load and U pe is the potential energy due to the artificial constraints as given in Eq.(3b).3.5.1.Kinetic energy of laminateThe displacement components at any point within the plate may be expressed in terms of reference plane parameters d asL¼U V Wf g T¼Q dð14aÞwhere Q is a function of z,unit step functions and the parameter X which is given in Appendix A.With the above equation,the velocity components at any point may be expressed as_L¼Q_dð14bÞWhere the dot(.)indicates derivative with respect to time‘t’.The kinetic energyðTÞof the vibrating laminated plate can be expressed asT¼12X nk¼1Zvf_L g T q k f_L g d v¼12X nk¼1ZZ_d T½I _d dxdyð15ÞWhere,½I ¼P nk¼1RQ T q k Qdz where q k is the mass density of the k thlayer.3.5.2.Strain energy of laminateThe strain energy for n number of layers of composite plate is given byU S¼12X nk¼1Zvf eg T½r d vð16aÞUsing Eqs.(2)and(5),the strain energy can be rewritten asU S¼12X nk¼1ZZZf eg T½ Q ðkÞf e g dxdydz¼1X nk¼1ZZf eg T½D f e g dxdyð16bÞWhere½D ¼X nk¼1Z½H T½ QðkÞ½H dzð16cÞ3.5.3.Potential energy due to applied in-plane loadsThe potential energy due to applied in-plane compressive loads using Eqs.(9)–(12),is expressed asU in¼1X nk¼1ZZZf h Gg T½S k fh G g dxdydz¼1X nk¼1ZZf eg T½G f e g dxdyð17aÞWhere,½G ¼P nk¼1R½H G T½S k ½H G dz.And the stress matrix[S k]may be expressed in terms of in-plane stress components(r xx,r yy,s xy)of the k th layer as½S k ¼r xx s xy0000sxyryy000000r xx s xy0000s xy r yy000000r xx s xy0000s xy r yy266666666437777777753.5.4.Work done due to external transverse loadThe work done by a distributed transverse static load of inten-sity p(x,y)isW ext¼ZZwpdxdyð17bÞWhere½p ¼½00p0000 considering p0as the trans-verse load on the plate.3.5.5.Penalty approachThe potential energy due to the artificial constraints U pe is expressed in terms of penalty parameter,‘c’andfield variables supported with artificial constraints using Eq.(3b),as follows:U pe¼c Z Z@w0Àh xT@w0Àhxþ@w0Àh yT@w0Àhy"#dxdyð18ÞBy the use of Eqs.(5)and(18)can be stated as following@w0@xÀh x&'¼P x d;@w0@yÀh y&'¼P y dð19ÞThen Eq.(18)becomes U pe¼c2R Rd T P TxP x dþd T P TyP y dh idxdy where,P x and P y are the vector containing shape functions and their deriva-tives.Also,is the generalized displacement vector defined in Eq. (3c).To make the proper balance between ill-conditioning of the320R.Sahoo,B.N.Singh/Composite Structures117(2014)316–332augmented stiffness matrix and loss of accuracy,the proper choice of penalty parameter ‘c ’is quite important.Thus,the value of penalty parameter ‘c ’is considered to be 106for all of the subsequent anal-ysis which is observed from different numerical experiences.The elemental potential energy as provided in Eq.(13)may be rewritten with the help of above equations,asP e ¼12f _dg T½M e f _d g À12f dg T ½K e f d g À12f dg T ½K Ge f d g À12f d g T f P e g À12f dg T ½K pe f d g ð20ÞWhere½M e ¼ZZ ½N T½I ½N dxdy ;½K e ¼ZZ½B T ½D ½B dxdy ;f P eg ¼ZZ ½N w T pdxdy½K Ge ¼ZZ½B G T½G ½B G dxdy ;½K pe ¼c ZZ P T x P x þP Ty P ydxdyThe [K e ],[K pe ]are element stiffness matrix and penalty matrixwhich are added for each element whereas [K Ge ],{P e }and [M e ]are the element geometric stiffness matrix,element load vectors and element mass matrix respectively while the terms without subscripts correspond to the global erning equationsThe governing equation of motion of the plate system can beobtained by applying Hamilton’s principle d R t2t 1P dt ¼0which finally yields the following equation½M f €d g þ½K þK G f d g ¼f 0g ð21a ÞNow the governing equations for free vibration and buckling analy-sis can be stated in the form of generalized eigenvalue problem as follows:3.6.1.Free vibration analysisWhere k is the critical load of buckling and [K ],[K G ],{U }are the global linear stiffness matrix,global geometric stiffness matrix and global displacement vector defining the buckling mode shape respectively.3.6.3.Static analysis½K f d g ¼f P g ð21d ÞWhere [K ],{d }and {P }are the global linear stiffness matrix,global displacement and global force vector respectively.4.Numerical examples and discussionNumerical computations followed with the free un-damped transverse vibration and stability analysis of laminated composite and sandwich plates were carried out to ensure the effectiveness and efficiency of newly developed ZZTSF.The selective integration scheme based on Gauss quadrature rules i.e.,reduced integration scheme (2Â2)for shear terms whereas full integration rule (3Â3)for membrane,flexure and inertia terms contributions,are adopted throughout the analysis.The encoded computational results in MATLAB environment are presented and compared with three dimensional elasticity results,analytical,FEM results and some already published results in the form of different tables and figures.The optimized value of parameter r is ascertained in post-processing technique by adopting the reverse method [27].Though the value of shear stress parameter has been chosen according to Ref.[30],ZZTSF also possesses results with minimum error compared to the 3D elasticity solution checked for a wide range of problems for the same value of r that is 0.1.The whole procedure for the estimation of r value is repeated for various span-thickness ratio,boundary conditions and structural analysis.However,the estimation of r parameter through vibration analysis and static analysis is presented in this study.The Fig.4indicates the estimation of r parameter through vibration analysis as well as static analysis of simply supported four layered symmetric cross-ply laminated plate.The non-dimensional fundamental fre- of r parameter through vibration and static analysis of simply supported laminated thick plate R.Sahoo,B.N.Singh /Composite Structures 117(2014)316–332321。