Xing Yufeng,Xiang Wei

aInstitute of Solid Mechanics,Beihang University,Beijing 100083,China

bSchool of Mechanical Engineering,Southwest Jiaotong University,Chengdu 610031,China

Analytical solution methods for eigenbuckling of symmetric cross-ply composite laminates

Xing Yufenga,*,Xiang Weia,b

aInstitute of Solid Mechanics,Beihang University,Beijing 100083,China

bSchool of Mechanical Engineering,Southwest Jiaotong University,Chengdu 610031,China

Accurate computation;Analytical solution;Composite laminates;Eigenbuckling;Separation of variables;State space

Based on the first-order shear deformation theory,this paper explores the analytical methods for eigenbuckling of symmetric cross-ply rectangular composite laminates with a pair of parallel edges simply supported and the remaining two edges arbitrarily constrained.The main contribution of present paper lies in two aspects:one is to present a simple and effective analytical method,namely,the separation-of-variables method,which can generate the closed-form buckling solutions without any computational difficulty;the other is to incorporate the accurate computation method of exponential matrix into the state space technique to avoid the inevitable numerically illconditioned problems reported in several literatures.The results obtained via both analytical methods are identical,and a good agreement with their counterparts in literature is observed.The separation-of-variables method can generate exact solutions within 1 s,which is impossible if the state space method is employed.Besides,the combination of the accurate computation method of exponential matrix and the state space method greatly improves the computational efficiency and gives correct results compared with the straightforward use of state space method.

1.Introduction

Laminated composite plates are increasingly used in aeronautical,mechanical,civil and marine structures.The buckling analysis of composite plates is essential for reliable and eff icient structural design.

Various theories,generally classified into single-layer and multilayer(layer-wise)plate theory,have been proposed to deal with eigenbuckling of composite laminates.Single-layer theories,in which a composite laminate is treated as an equivalent orthotropic and homogeneous single layer,are often used to predict global response characteristics such as maximum deflections,natural frequencies and critical buckling loads,while multilayer theories are indispensable to model local effects such as delamination.Thus,for the purpose of present study,only single-layer plate theories are discussed.

The classical laminated plate theory based on the Kirchhoff assumptions was employed for buckling analysis of orthotropic plates by Das,1Harik and Ekambaram,2Bao et al.,3Hwang and Lee,4and others.This theory neglects the transverse shear strains,thus overpredicting the buckling loads or even yielding incorrect results especially for composite plates with relatively high ratio of in-plane elastic modulus to transverse shear modulus.

Composite plates are highly flexible in transverse shear;therefore the transverse shear deformation must be taken into account to achieve an accurate representation of laminated plate behavior.Many shear deformation plate theories,which can be extended for analysis of laminated composite plates,are documented in Refs.5–15.

The first-order shear deformation plate theory(FSDT),commonly known as the Mindlin plate theory(MPT),7is representative and one of the most widely used shear deformation theory and provides the best compromise between accuracy and computational efficiency.

Higher-order shear deformation plate theories(HSDT),which involve higher order terms in Taylor’s expansion of the displacements in the thickness coordinate,have been proposed by many researchers including Librescu,9Levinson,10Bhimaraddiand Stevens,11Reddy,12Ren,13Kantand Pandya,14Shimpi and Patel15and so on.The third-order shear deformation theory by Reddy12is representative,in which the shear stress function is parabolic through the thickness and satisfies the stress-free boundary conditions;therefore,shear correction factor is not required.

Considering symmetric cross-ply composite laminates as an equivalent single layer,the standard procedure of determining critical buckling loads is to solve a mathematical eigenvalue problem governed by characteristic differential equations and boundary conditions.However,the exact solutions are limited to a relatively few cases wherein at least two parallel edges are simply supported.For other combinations of edge conditions,approximate procedures such as the Ritz,Galerkin,superposition,finite element and finite difference methods16–23should be used.Numerical methods are effective in analyzing plates of arbitrary boundary conditions,physical properties and loading configurations.

Although research regarding the extension of various numerical methods to buckling analysis of composite plates is increasing rapidly,a relative lack of theoretical value and substantial progress exists.In contrast,the development of analytical method is rare,but analytical solution is essential for understanding the physical behavior of plates.

The state space method,which has been widely used in several literatures,24–36is a powerful analytical technique to generateLevy-typesolutionsforeigenbuckling ofcomposite laminates.However,the main drawback of this method is the inevitable computational difficulty caused by the calculation of exponential matrix.As was pointed out by Khdeir33that while shearing deformation theories are used,straightforward application of the state-space concept yields numerically illconditionedproblemsasthelaminatethicknessisreduced.Software with double precision floating point calculation ability may fail to work out the results and the computation time and workload are unsatisfactory.To overcome this computational difficulty,various mathematical methods including the decomposition of matrix and the modified Gram-Schmidt orthonormalization procedure were discussed34–36for free vibration and buckling problems,and these methods are not the straightforward use of state space method although they are effective.

In this context,this paper is to develop analytical methods to overcome or even totally avoid computational difficulty for the straightforward use of state space method.The main contribution of present paper lies in two aspects.First,a simple and effective analytical method,which can totally avoid any computational difficulty,is proposed.This method,referred to as separation-of-variables method in this paper,affords explicit eigenvalue equations and closed-form buckling solutions for Levy-type composite plates.Furthermore,it takes less than 1 s to work out the results and the workload is rather small.Second,a combination of the accurate computation method of exponential matrix and the state space concept is introduced to completely solve the inevitable numerically illconditioned problems aforementioned,greatly reducing computation time and workload.

Closed-form buckling solutions of symmetric cross-ply rectangular laminates based on FSDT are obtained by both analytical methods and a comparison study is conducted to validate both methods.

2.Governing equations and boundary conditions

A laminated rectangular composite plate of lengtha,widthband uniform thicknesshis considered oriented so that its mid-plane surface contains thex-andy-axis of a Cartesian coordinate system(x,y,z),as shown in Fig.1.

The displacement field of the first-order shear deformation theory is

whereu0(x,y),v0(x,y)andw0(x,y)denote the displacements of a point(x,y)in mid-plane;ψxand ψyare the rotations of a normal line with respect toyandxcoordinates,respectively.

The strain field for the assumed displacement field follows immediately as

in which

The governing equations and boundary conditions can be derived from thefollowing potentialenergy variational principle:

where Π,UandVdenote total potential energy,strain energy and potential energy due to applied loads,respectively;δ is the symbol of variation.

whereNxandNyare constant initial in-plane edge loads per unit length in thex-andy-directions,respectively.Integrating Eq.(4)in parts,we have

Collecting the coefficients of δu0,δv0,δw,δψxand δψy,we obtain the following differential equations of motion and boundary conditions:

where the stress resultantsNi,Mi(i=1,2,6)andQi(i=1,2)are defined as

The constitutive equations for each layer can be expressed in the material axes as

whereQij(i,j=1,2,6)are the plane-stress reduced elastic constants in the material axes of the layer

The stiffness equations relating in-plane stress resultants to strains are given as

The coefficients in Eq.(13)are given as

whereQij(i,j=1,2,6)are the elastic constants after transformation into the plate(laminate)coordinate.

It can be found from Eq.(13)that,for symmetrically laminated plates,Bijvanish,leaving bending and stretching uncoupled.For cross-ply plates,we haveA16=A26=D16=D26=C45=0.Thussymmetriccross-plylaminatesare governed by a single-layer orthotropic plate theory,the constitutive equations of which are

It is assumed thatNx= χ1Ncr,Ny= χ2Ncr,where χ1and χ2are defined as scaling parameters andNcris the critical buckling load.And the following parameters are defined for brevity:

Then the governing equations can be rewritten as

The specific formulations of the boundary conditions for edgesy=0 andy=bare

(1)Simply supported edge(S)

(2)Clamped edge(C)

(3)Free edge(F)

Hereinafter,S,C and F denote simply supported,clamped and free edges,respectively.For the purpose of describing boundary conditions briefly,the notation SCSF is used to denote a plate with edgesx=0 andx=asimply supported,y=0 clamped,andy=bfree.

3.Analytical solutions

In this paper,it is assumed that the two opposite edgesx=0 andx=aare simply supported and the other two edges are arbitrarily constrained such as free,simply supported and clamped.To obtain exact buckling solutions for composite laminates,two methods,namely,separation-of-variables method and state space method,are employed.And the detailed procedure and a further discussion on the efficiency of two methods are presented in this section.

3.1.Separation-of-variables method

For convenience,Eq.(17)is rewritten as

By elimination method,Eq.(21)can be transformed into

The following solution for the transverse deflection is adopted to satisfy the boundary conditions of the two simply supported edgesx=0 andx=a.

The roots of Eq.(25)can be written in the following form:

Then the mode functionsw, ψxand ψyare expressed in terms of the eigenvalues as

The coefficientsgjandhj(j=1,2,3)can be derived by substituting Eq.(28)into Eq.(17).

Taking boundary conditions of the remaining two edgesy=0 andy=binto account,we can easily obtain the closed-form explicit eigenvalue equation and the coefficientsAjandBjof mode functions.Here we take the case SSSC as an example to show the detailed procedure.Substituting Eq.(28)into the boundary conditions ofy=0(S)andy=b(C)as shown in Eqs.(18)and(19),we have

On condition that Eq.(31)has nontrivial solution,the eigenvalue equation for this case SSSC is

And from the set of algebraic equations(31),we can derive the integral coefficientsAi,Bi(i=1,2,3)of mode functions.Tables 1 and 2 present the explicit eigenvalue equations and the integral coefficients of mode functions for different combinations of boundary conditions,respectively.Note that ‘(a)’represents anti-symmetrical mode shapes, ‘(s)’symmetrical mode shapes,and the following symbols are used for brevity:

The critical buckling factors can be obtained by solving the transcendental eigenvalue equation in conjunction with the characteristic equation as shown in Eq.(25).An interactive procedure is needed.Given an initial value of critical buckling loadNcr,and through only a few iterations,we can work out the critical buckling factor without any calculation difficulty.However,the drawback of this method is that it is rather sensitive to the initial value especially when the thickness grows.It should be stressed that the floating point calculation ability(with 16 significant digits)can perfectly meet requirements,and the calculation time is within 1 s.

3.2.State space method

The following representation for the displacement quantities is chosen to automatically satisfy the boundary conditions of the two edges normal to thex-axis.

Table 1 Eigenvalue equations for different boundary conditions(BC).

Table 2 Integral constants Ai,Bi(i=1,2,3)for different boundary conditions.

whereW(y),X(y)andY(y)are unknown functions to be determined.

Substitute Eq.(34)into Eq.(17),and the governing differential equations can be transformed into a system of linear,ordinary differential equations as

where the primes indicate ordinary differentiation with respect toy.Applying the concept of state space,we can convert Eq.(35)into matrix form:

in which ψ =[W W′X X′Y Y′]T,and H is a matrix of dimension 6×6 with the following form:

The elements of matrix H are listed as follows,while the ones that are not given equal zeroes.

A general solution of Eq.(36)is given by

in which c is a constant column vector associated with the boundary conditions of the remaining two edges.

It is noteworthy that the computation of matrix exponential eHywill meet difficulties,and thus in this paper the accurate computation method of exponential matrix37is employed to avoid serious round-off error.The accurate computation of exponential matrix has two key points,namely,the additional theorem of exponential function such as the 2Nalgorithm,and keeping the incremental part of the exponential matrix,rather than the total value.

The incremental part of the exponential matrix Tais given as

where η=y/2N,Nis an arbitrary integer,and I6is an unit matrix of dimension 6×6.

Execute program statement Ta=2Ta+T2aNtimes,and the exponential matrix eHywith high precision is generated as

Note that,in order to avoid calculating difficulty,thex-axis of coordinate system should be moved to the horizontal axis of symmetry of the plate so that the remaining two edges are denoted asy=±b/2.Substituting the appropriate boundary conditions of the edgesy=±b/2 into Eq.(39)yields a homogeneous system of equations as

where K is a 6×6 matrix.

The buckling loadNcrcan be determined by setting the determinant of K to zero.It should be noted that this solution procedure cannot provide buckling load directly because the undetermined buckling loadNcris included in matrix K.Thus an iterative numerical procedure,which has been stated in several literatures,25,32needs to be used to calculate the buckling load.

3.3.Discussion

The separation-of-variables method and state space method can deal with the elastic buckling of composite laminates effectively and yield the same critical buckling factors,which will be validated in the following section.Additionally,both methods haveidentical physical essence,which isto exactly solve amathematical eigenvalue problem governed by differential equations and boundary conditions.However,in the separation-ofvariables method,this problem is to solve a transcendental equation,while the key of the state space method is to make a homogeneous system of equations yield nontrivial solutions by setting the determinant of the coefficient matrix to zero.

It should be noted that when the separation-of-variables method is employed,the calculation difficulty can be totally avoided.The floating point calculation ability can perfectly meet requirements,and the calculation time is more than satisfactory.

However,computation difficulty caused by calculation of matrix exponential might be encountered when state space method is applied.As Xiang et al.25pointed out,when software with double precision floating point calculation ability,such as some versions of Fortran and Mathematica,is employed,incorrect solutions might be generated for plates with large aspect ratios and small thickness ratios.Although this problem is eventually solved by using computer software with more significant digits(more than 100 significant digits if necessary)in floating point calculation,the computation time and workload are not satisfactory.Thus,in present paper,an improvement to the traditional state space method is made by applying the accurate computation method of exponential matrix aforementioned,and only forty significant digits are needed in MATLAB,thus greatly improving computational efficiency,which will be verified in the next section.

4.Numerical results and comparison

In this section,numerical results are presented for composite square laminates with two parallel edges simply supported and the other two edges having arbitrary boundary conditions.Firstly,a three-layered symmetric cross-ply plate with the stacking sequence 90°/0°/90°is considered.Each layer has the same thickness,and the engineering constants of each layer are

By employing the two analytical methods aforementioned,the exact critical buckling factorskunder biaxial(χ1=-1,χ2=-1)and uniaxial(χ1=-1,χ2=0)loadings are generated,and tabulated in Tables 3 and 4,respectively.A comparison is performed between present results and their counterparts in literatures written by other researchers.Note that in Tables 3 and 4,the separation-of-variables method and the state space method are denoted as separation-ofvariables method(SVM)and state space method(SSM),respectively.

It can be clearly seen from Table 3 that the SVM results and SSM results of present analysis,as well as published results by Xiang et al.25using FSDT,are completely identical,which demonstrates the validity of present approaches.We can also find an excellent agreement between the present analytical results and published results by Khdeir and Libreseu26based on HSDT.Table 3 also reveals that HSDT generally gives critical buckling loads higher than FSDT.

It is pointed out in the last section that straightforward application of state space method yields numerically illconditioned problems,and therefore computer software with 100 significant points is needed to generate correct solutions.However,after the introduction of accurate computationmethod of exponential matrix,only 40 significant digits are needed to meet computational requirements.

Table 3 Comparison of buckling factorssymmetric cross-ply square laminates subjected to biaxial compression(χ1=-1,χ2=-1).

Table 3 Comparison of buckling factorssymmetric cross-ply square laminates subjected to biaxial compression(χ1=-1,χ2=-1).

Notes: SVM — Seperation-of-variable method; SSM — State space method.

BC Method h/b 0.05 0.10 0.20 SSSS SVM 13.1472 10.2024 5.4920 SSM 13.1472 10.2024 5.4920 Ref.26 13.1850 10.2590 5.5260 SSSC SVM 17.9053 11.6015 5.8763 SSM 17.9053 11.6015 5.8763 Ref.26 17.9170 11.6340 5.8880 SCSC SVM 24.5426 13.2898 6.1096 SSM 24.5426 13.2898 6.1096 Ref.26 24.4520 13.2880 6.1630 SSSF SVM 1.9717 1.8971 1.6587 SSM 1.9717 1.8971 1.6587 Ref.25 1.9717 1.8971 SCSF SVM 4.4695 3.7548 2.8182 SSM 4.4695 3.7548 2.8182 Ref.25 4.4694 3.9586 SFSF SVM 1.9301 1.8268 1.7037 SSM 1.9301 1.8268 1.7037 Ref.25 1.9301 1.8268 Notes:SVM—Seperation-of-variable method;SSM—State space method.

Table 4 Comparison of buckling factorssymmetric cross-ply square laminates subjected to uniaxial compression(χ1=-1,χ2=0).

Table 4 Comparison of buckling factorssymmetric cross-ply square laminates subjected to uniaxial compression(χ1=-1,χ2=0).

BC Method k h/b=0.05 h/b=0.10 h/b=0.20 SSSS Present 31.1959 22.3151 10.9839 SSSC Present 50.6142 28.7820 12.0922 SCSC Present 78.9066 36.6185 11.7379 SSSF Present 8.0309 7.5179 5.7379 Ref.25 8.0310 7.5179 SCSF Present 12.6748 11.1973 7.3561 Ref.25 12.6740 11.1970 SFSF Present 3.5057 3.3055 2.7504 Ref.25 3.5056 3.3055

Table 5 Computation time of buckling factors of bi-axially loaded(χ1=-1,χ2=-1)square laminates(h/a=0.1).

To show that the computation efficiency is improved by the incorporation of accurate computation method of exponential matrix,we set the number of significant digits of MATLAB to 40 and 100 separately,and calculate the critical buckling factors of bi-axially loaded square laminates by using the same computer and running the same program.The runtime is presented in Table 5,from which one can see that runtime is greatly shortened when the number of significant digits is set to 40,thus verifying the improvement of computation eff iciency.Thecomputation timeofseparation-of-variables method is also given in Table 5 for comparison.We can find that separation-of-variables method has incomparable advantages in computation efficiency over state space method.

Note that the runtime in Table 5 with 40 and 100 significant digits is under the condition that the accurate computation method of exponential matrix is involved,otherwise straightforward computation of exponential matrix with even 100 significant digits will take an unbearable long time,and is hard to work out correct solutions.

Buckling mode shapes by two analytical methods for Levytype laminates at thickness ratioh/a=0.05 under uniaxial and biaxial loadings are given in Fig.2.

Secondly,a single orthotropic layer of the following material property is considered:

Table 6 Comparison of buckling factorsfor orthotropic plates subjected to biaxial compression(χ1=-1,χ2=-1).

Table 6 Comparison of buckling factorsfor orthotropic plates subjected to biaxial compression(χ1=-1,χ2=-1).

E1/E2 BC Method k h/b=0.05 h/b=0.10 h/b=0.20 10 SSSS Present 5.3100 4.6367 2.8319 Ref.38 5.2924 4.5650 2.7547 SSSC Present 7.4502 5.7130 3.0705 Ref.38 7.3938 5.6086 3.0139 SCSC Present 10.4475 7.0952 3.3530 Ref.38 10.3767 7.0006 3.3260 40 SSSS Present 9.3049 6.6325 3.2822 Ref.38 9.2368 6.5168 3.2280 SSSC Present 12.3540 7.6616 3.3778 Ref.38 12.2197 7.5105 3.3611 SCSC Present 15.3177 8.3023 3.4503 Ref.38 15.2219 8.2124 3.4665

The results are presented in Table 6 and a comparison is made between present results based on FSDT and those in Ref.38using layer-wisetheory,and an excellentagreementis observed.

5.Conclusions

The eigenbuckling behavior of symmetric cross-ply laminates subjected to uniaxial and biaxial uniform loadings with one pair of opposite edges simply supported was investigated in this paper.The first-order shear deformation theory was adopted to take transverse shear deformation into consideration.

A simple and effective analyticalmethod,namely,separation-of-variables method,which can totally avoid any computational difficulty,was presented.This method affords explicit eigenvalue equations and closed-form buckling solutions for Levy-type composite plates within 1 s.Moreover,the accurate computation method of exponential matrix was incorporated into the state space concept to overcome the numerically ill-conditioned problems,greatly reducing computation time and workload.

Closed-form buckling solutions of symmetric cross-ply rectangular laminates obtained by both analytical methods are identical,and also in good agreement with their counterparts in literatures;this comparison validated present methods.

Acknowledgements

This study was co-supported by the National Natural Science Foundation of China(Nos.11372021 and 11672019)and Research Fund for the Doctoral Program of Higher Education of China(No.20131102110039).

1.Das Y.Buckling of rectangular orthotropic plates.Appl Sci Res1963;11(1):97–103.

2.Harik I,Ekambaram R.Elastic stability of orthotropic plates.Thin-Wall Struct1988;6(5):405–16.

3.Bao G,Jiang W,Roberts JC.Analytic and finite element solutions for bending and buckling of orthotropic rectangular plates.Int J Solids Struct1997;34(14):1797–822.

4.Hwang I,Lee J.Buckling of orthotropic plates under various inplane loads.KSCE J Civil Eng2006;10(5):349–56.

5.Yang PC,Norris CH,Stavsky Y.Elastic wave propagation in heterogeneous plates.Int J Solids Struct1966;2(4):665–84.

6.Reissner E.The effect of transverse shear deformation on the bending of elastic plate.ASME J Appl Mech1945;12:69–77.

7.Mindlin RD.Influence rotatory inertia and shear in flexural motion ofisotropic,elasticplates.ASMEJApplMech1951;18:31–8.

8.Whitney JM,Pagano NJ.Shear deformation in heterogeneous anisotropic plates.ASME J Appl Mech1970;37:1031–6.

9.Librescu L.On the theory of anisotropic elastic shells and plates.Int J Solids Struct1967;3(1):53–68.

10.Levinson M.An accurate,simple theory of the statics and dynamicsofelasticplates.MechResCommun1980;7(6):343–50.

11.Bhimaraddi A,Stevens L.A higher order theory for free vibration of orthotropic,homogeneous,and laminated rectangular plates.J Appl Mech1984;51(1):195–8.

12.Reddy JN.A simple higher-order theory for laminated composite plates.J Appl Mech1984;51:745–52.

13.Ren J.A new theory of laminated plate.Compos Sci Technol1986;26(3):225–39.

14.Kant T,Pandya B.A simple finite element formulation of a higher-order theory for unsymmetrically laminated composite plates.Compos Struct1988;9(3):215–46.

15.Shimpi RP,Patel HG.A two variable refined plate theory for orthotropic plate analysis.Int J Solids Struct2006;43(22–23):6783–99.

16.Chai GB.Buckling of generally laminated composite plates with various edge support conditions.Compos Struct1994;29(3):299–310.

17.Dickinson SM,Di Blasio A.On the use of orthogonal polynomials in the Reyleigh-Ritz method for the study of the flexural vibration and buckling of isotropic and orthotropic rectangular plates.J Sound Vib1986;108(1):51–62.

18.Gorman DJ,Ding W.Accurate free vibration analysis of laminated symmetric cross-ply rectangular plates by the superposition-Galerkin method.ComposStruct1996;34(4):387–95.

19.Bhaskar K,Sivaram A.Untruncated infinite series superposition method for accurate flexural analysis of isotropic/orthotropic rectangular plates with arbitrary edge conditions.Compos Struct2008;83(1):83–92.

20.Kshirsagar S,Bhaskar K.Accurate and elegant free vibration and buckling studies of orthotropic rectangular plates using untruncated infinite series.J Sound Vib2008;314(3–5):837–50.

21.Tsay CS,Reddy JN.Bending,stability and free vibrations of thin orthotropic plates by simplified mixed finite elements.J Sound Vib1978;59(2):307–11.

22.Putcha NS,Reddy JN.Stability and natural vibration analysis of laminated plates by using a mixed element based on a refined plate theory.J Sound Vib1986;104(2):285–300.

23.Yuan S,Jin Y.Computation of elastic buckling loads of rectangular thin plates using the extended Kantorovich method.Comput Struct1998;66(6):861–7.

24.Reddy JN,Khdeir AA,Librescu L.Levy type solutions for symmetrically laminated rectangular plates using first-order shear deformation theory.J Appl Mech1987;54(3):740–2.

25.Xiang Y,Liew KM,Kitipornchai S.Exact buckling solutions for composite laminates:Proper free edge conditions under in-plane loadings.Acta Mech1996;117(1):115–28.

26.Khdeir AA,Libreseu L.Analysis of symmetric cross-ply elastic plates using a higher-order theory.Part II:Buckling and free vibration.Compos Struct1988;9(4):259–77.

27.Khdeir AA.Free vibration and buckling of symmetric cross-ply laminated plates by an exact method.J Sound Vib1988;126(3):447–61.

28.Khdeir AA,Reddy JN,Librescu L.Analytical solutions of a refined shear-deformation theory for rectangular composite plates.Int J Solids Struct1987;23(10):1447–63.

29.Liew KM,Xiang Y,Kitipornchai S.Analytical buckling solutions for Mindlin plates involving free edges.Int J Mech Sci1996;38(10):1127–38.

30.Reddy JN,Khdeir AA.Buckling and vibration of laminated composite plates using various plate theories.AIAA J1989;27(12):1808–17.

31.Chen WC,Liu WH.Deflection and free vibrations of laminated plates-Levy-type solutions.Int J Mech Sci1990;32(9):779–93.

32.Thai HT,Kim SE.Levy-type solution for buckling analysis of orthotropic plates based on two variable refined plate theory.Compos Struct2011;93(7):1738–46.

33.Khdeir AA.A remark on the state-space concept applied to bending,buckling and free vibration of composite laminate.Comput Struct1996;59(5):813–7.

34.Khdeir AA,Reddy JN.Free vibration of cross-ply laminated beams with arbitrary boundary conditions.J Eng Sci1994;32(12):1971–80.

35.Nosier A,Reddy JN.Vibration and stability analyses of cross-ply laminated circular cylindrical shells.J Sound Vib1992;157(1):139–59.

36.Hadian J,Nayfeh AH.Free vibration and buckling of sheardeformable cross-ply laminated plates using the state-space concept.Comput Struct1993;48(4):677–93.

37.Zhong WX.On precise integration method.J Comput Appl Math2004;163(1):59–78.

38.Vescovini R,Dozio L.Exact refined buckling solutions for laminated plates under uniaxial and biaxial loads.Compos Struct2015;127:356–68.

4 July 2016;revised 20 July 2016;accepted 5 August 2016

Available online 22 December 2016

?2016 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is anopenaccessarticleundertheCCBY-NC-NDlicense(http://creativecommons.org/licenses/by-nc-nd/4.0/).

*Corresponding author.

E-mail addresses:xingyf@buaa.edu.cn(Y.Xing),xiangweivee@163.com(W.Xiang).

Peer review under responsibility of Editorial Committee of CJA.