Vibration and Buckling Analysis of Functionally Graded Plates Using New Eight-Unknown Higher Order Shear Deformation Theory

In this paper a new eight-unknown higher order shear deformation theory is proposed to study the buckling and free vibration of functionally graded (FG) material plates. The theory bases on full twelve-unknown higher order shear deformation theory, simultaneously satisfies zero transverse shear stress at the top and bottom surfaces of FG plates. Equations of motion are derived from Hamil-ton’s principle. The critical buckling load and the vibration natural frequency are analyzed. The accuracy of present analytical solution is con ﬁ rmed by comparing the present results with those available in existing literature. The effect of power law index of functionally graded material, side-to-thickness ratio on buckling and free vibration responses of FG plates is investigated.


NOMENCLATURE (continuation)
x y z u v w 0 0 0 , , , , , q q q * * * * * * higher-order terms of displacements in the Taylor series expansion w , ŵ natural frequency and non-dimensional natural frequency cr cr N N ,  critical buckling load and non-dimensional critical buckling load ij Q stiffness coefficients of FG plates x y xy N 0 0 0 , N , N in-plane pre-buckling loads 1 2 , g g in-plane load parameters

INTRODUCTION
Ever since invented by Japanese scientists in the 80s of the last century, Functionally Graded Materials (FGMs) has been more and more widely applied in many fields such as aircraft industry, nuclear industry, civil engineering, automotive, biomechanics, optics… Typical FGMs are composed of ceramic and metal materials.Ceramic provides high temperature resistance while metals have high toughness; thus FGMs are usually used in the manufacture of heat-resistance structural components such as airplane fuselages or walls in nuclear reaction plants.…Theunderstanding of the behavior of FGM, therefore, is very much desired.Studying the static and dynamic behavior of FGM structures has become an interesting topic for researchers around the world.
There have been many computational models and methods of calculation proposed for FGM plates.The classical plate theory (CPT) that bases on Kirchhoff-Love's assumption, is only suitable for thin plates since it ignores the effects of transverse shear deformation.For moderately thick plates, numerical results calculated using CPT yield lower deflection, higher natural frequency and buckling load in comparison with experimental results.In order to correct this inaptitude, the first-order shear deformation theories (FSDT) have been initially proposed by Reissner and further developed by Mindlin.Although FSDT describes more realistic behavior of thin to moderately thick plates, the parabolic distribution of transverse shear stress through the thickness of the plate is not properly reflected, thus the shear correction factor is introduced.The determination of this factor is not simple as it depends on the loading, boundary condition, materials etc… To avoid using shear correction factor, higher order shear deformation theories (HSDTs) are proposed.Based on third order shear deformation theory with five displacement unknowns, Reddy (2000) developed analytical and finite element solutions for static and dynamic analysis of functionally graded rectangular plates.The formulation accounts for the thermo-mechanical coupling, time dependency, and the von Kárman-type geometric non-linearity.Bodaghi et al. (2010) used Reddy's third order shear deformation theory and Levy-type solution for buckling analysis of thick functionally graded rectangular plates.Also with five displacement unknowns, Zenkour (2006) used his generalized shear deformation theory to study static behaviors of simply supported functionally graded rectangular plate subjected to a transverse uniform load.Employing finite element method based on nine Latin American Journal of Solids and Structures 13 (2016) 456-477 unknowns higher order shear deformation theory, Pandya and Kant (1988) investigated deflections, in-plane and inter-laminar stresses of thick laminated composite plates.This displacement model assumed non-linear and constant variation of in-plane and transverse displacement, respectively, through the plate thickness.Using eleven-unknown displacement field and finite element method, Talha and Singh (2010) analyzed static response and natural frequency of functionally graded plates.Higher order terms of the displacement field are determined by vanishing the transverse shear stresses on the top and bottom surfaces of the plate.Kim and Reddy (2013) also used acouple of stress-based third-ordertheory with elevenunknowns to analyze the bending, vibration and buckling behaviors of FG plates by analytical method.Based on the higher-order refined theories, Jha et al. (2012) presented analytical solutions for free vibration analysis of simply supported rectangular functionally graded plates.This HSDT introduces twelve displacement unknowns, and correctly describes the quadratic distribution of transverse normal strain across the thickness although the values at the top and bottom are non-zero.A comprehensive review of the various methods employed to study the static, dynamic and stability behaviors of functionally graded plates can be found in the work of Swaminathan et al. (2015).The review focuses on comparing the stress, vibration and buckling characteristics of FGM plates using different theories.It is observed that most of the above mentioned HSDTs require additional computation efforts due to the additional unknowns introduced to them (usually nine, eleven or thirteen unknowns depending on each particular theory).
In this paper, a new higher order displacement field based on twelve-unknown higher order shear deformation theory is developed to analyze the free vibration and buckling of functionally graded plates.The new form is dictated by the satisfaction of vanishing transverse shear stress at the top and bottom surfaces of the plate.With this proposed higher order displacement field, the number of displacement unknowns reduces from twelve to eight, thus savingcomputational time and optimizing the storage capacity of computers.The accuracy of the present theory is verified by comparisonwith previous studies.

KINEMATICS
The displacement components u (x,y,z), v(x,y,z) and w(x,y,z) at any point in the plate can be expanded in Taylor's series in terms of the thickness coordinate as (Jha -2012) where denote the displacements of a point along the (x,y, z) coordinates.u v w 0 0 0 , , are corresponding displacements of a point on the mid-plane.x q , y q and z q are rotations of transverse normal to the mid-plane about the y-axis, x-axis and z-axis, respectively.
x y u v w 0 0 0 , , , , q q the higher-order terms in the Taylor series expansion and they represent higher-order transverse crosssectionaldeformation modes.
Latin American Journal of Solids and Structures 13 (2016) 456-477 For plates under bending, the transverse shear stresses xz s , yz s must be vanished at the top and bottom surfaces.These conditions lead to the requirement that the corresponding transverse strains on these surfaces to be zero.From , we obtain: ; ; ; .
The displacement field (1) becomes: with: Using the strain-displacement relations of the theory of elasticity, the linear strains are obtained: where: e e e g q k k k k q q q q e e e g q q q q q q k k k q q q q g g q q k k q q g g q q k k q q In the above formulas, a comma followed by x or y denotes differentiation with respect to the coordinates x or yrespectively.

CONSTITUTIVE EQUATION
Consider a simply supported linearly elastic rectangular FG plate of uniform thickness h as shown in Figure 1.The Poisson's ratio n is assumed to be constant across the plate thickness.The Young's modulus, the mass densityof the FG plate is assumed to follow the power law distribution alongthe thickness, and expressed as (Reddy -2000): In the above formula, subscript c refers to ceramic material and subscript m refers to metal material of the FG plate.It is clear from the expression that the top surface ( z h / 2 = ) of the FG plate is ceramic-rich and the bottom ( z h/ 2 = - ) is metal-rich in constituents.The stress-strain relationship for the FG plate can be written as: in which ( ) .

EQUILIBRIUM EQUATIONS
Hamilton's principle is used to derive the equations of motion.The principle can be stated in analytical form as: ( ) where δU is the variation of strain energy; δV is the variation of external work; and δKis the variation of kinetic energy.The variation of strain energy of the plate can be calculated by: ( ) The variation of work done by in-planeand transverse loads is given by: where and z q + is the transverse load at the top surface of the plate, is the transverse displacement of any point on the top surface of the plate; x y xy N 0 0 0 ,N ,N are in-plane pre- buckling applied loads.The variation of kinetic energy of the plate can be written in the form: where the dot-superscript indicates the differentiation with respect to the time variable t; ( ) is the mass density.

NAVIER'S SOLUTION
Consider a simply supported rectangular FG plate with length a, width b under in-plane loads in two ).The associated simply supported boundary conditions are as follows: At edge x = 0 and x = a: At edge y = 0 and y = b: Latin American Journal of Solids and Structures 13 (2016) 456-477 Following Navier's solution procedure, the displacement variables are chosen to satisfy the above simply supported boundary condition with the form(for the buckling and vibration problems, the transverse load is set to be zero): where: i 1 =is the imaginary unit.
where the elements of matrix [S], [M] are defined in the Appendix, and: The system of Eq. (13a) maybe used to obtain the solutions of the buckling problem soft he FG plates by dropping all the inertia terms ( For convenience, the following non-dimensional forms are used: ( ) In order to emphasize the efficiency of present eight-unknown HSDT, the calculated results are compared with other shear deformation theories.The following models of shear deformation theories are used in this section: HSDT-12: HSDT-5: HSDT-4: where

Buckling Analysis
Example 1. Functionally gradedAl/Al2O3 square (b/a=1) and rectangular (b/a=2) plates subjected to biaxial compression (γ1 = -1, γ2 = -1) are considered.Table 1 gives some numerical results showing the accuracy of the present non-dimensional buckling loads with various values of side-to-thickness ratio.The obtained results based on proposed HSDT are compared with the results of Thai and Choi (2012) whichwere based on an efficient and simple refined theory.The theory, which Thai and Choi used is similar with the classical plate theory in many aspects; it accounts for a quadratic variation of the transverse shear strains across the thickness and satisfies the zero traction boundary conditions on the top and bottom surfaces of the plate.A good agreement between the results is observed.Example 2. In this example, a moderately thick (a/h = 10) rectangular (b/a = 2) FG plate with different values of power law index p is examined.Table 2 contains the non-dimensional buckling loads calculated by present and various shear deformation theories: first-order shear deformation theory with 5 unknowns (FSDT), third-order shear deformation theory with 5 unknowns (HSDT-5), simple higher-order shear deformation theory with 4 unknowns (HSDT-4), and full higher-order shear deformation theory with 12 unknowns (HSDT-12).Fig. 2   Example 3. Thick and thin rectangular FG plates (p = 5)with side-to-thickness ratio varies from 5 to 100 are analyzed using present HSDT and various shear deformation theories.The non-dimensional buckling loads cr N ˆ under uniaxial and bi-axial compression are presented in Table 3. Fig. 3 shows a variation of non-dimensional buckling loads cr N ˆ with respect to side-to-thickness ratio a/hof rectangular FG plates (b/a = 2, p = 5).It can be seen that the non-dimensional buckling load increases by the increase of thickness ratioa/h, and the variation of the non-dimensional buckling loadbecomes significant for thick plate.The difference in the results obtained using proposed HSDT and the rest of HSDT increases with a decreases in the value of the side-to-thickness ratioa/h.It is emphasized that the proposed HSDT model contains a fewer number of unknowns than those associated with the full HSDT theory.However, an excellent agreement between the results predicted by present HSDT and full HSDT with 12 displacement unknowns also can be observed, and FSDT overestimates the buckling loads of FG thick plate as it neglects the thickness stretching effect.

Free Vibration Analysis
Example 4. The next verification is performed for moderately thick and thick FG square plates.Different values of power law index are considered.The non-dimensional fundamental frequencies are given in Table 4. Obtained results are compared with solutions using first-order and higher-order shear deformation theories provided by Hosseini-Hashemi (2011), and sinusoidal shear deformation theory provided by Thai (2013).It can be seen that the difference between the results is very small.Example 5. Non-dimensional frequencies ŵ of moderately thick rectangular plates (b/a = 2, a/h = 10) for different values of power law index p and various modes of vibration are presented in Table 5. Figure 4 illustrates the variation of non-dimensional fundamental frequency (m=n=1) with respect to power law index p.As can be seen from the presented results, the non-dimensional natural frequencydecreases with increasing value of power law index p.It is basically due to the fact that Young's modulus of ceramic is higher than metal.For the same value of p, the non-dimensional natural frequency increases for higher modes.Figure 4 also shows that the non-dimensional naturalfrequency decreases significantly when p is small.Table 6 shows non-dimensional frequencies ŵ of thin to thick rectangular FG plates (b/a = 2, p = 5) for different values of side-to-thickness ratio a/hand various modes of vibration.Figure 5 depict the variation of non-dimensional fundamental frequency (m=n=1) with respect to side-thickness-ratio a/h.
Similarly it is observed that the non-dimensional frequency decreases as the side-to-thickness ratio decreases.The fall in non-dimensional frequency is observed up to around a/h =20, beyond this no changes in non-dimensional frequency are distinguished.All above obtained results are studied using different plate theories.From table 5 and 6, it is apparent that the present proposed HSDT and full HSDT with 12 displacement unknowns give almost identical results for all values of power law index p and side-to-thickness ratio a/h.This emphasizes again, the benefits of the proposed HSDT in comparison with the full HSDT, as the proposed HSDT uses fewer displacement unknowns but requires less computational effort.

CONCLUSIONS
The new eight-unknown HSDT is proposed based on full twelve-unknown HSDT and satisfies vanishing transverse stresses at the top and bottom surface of FG plates.The accuracy of numerical solutions has been validated against existing results in available literatures.The effects of the side-tothickness ratio and the power law index of constituent volume fraction on the buckling loads and on the natural frequencies are also discussed.The results show that the buckling loads increase, and the natural frequencies decrease significantly with increasing power law index.It can be observed by the presented results that the gradation of the constitutive components is an important parameter Latin American Journal of Solids and Structures 13 (2016) 456-477 for buckling and free vibration analysis of FG plates.The present formulation for FG plates involves less computation compared to full twelve-unknown higher-order shear deformation theory,while gives identical results as full twelve-unknown higher-order shear deformation theory.The numerical results of critical buckling loads and natural frequencies should serve as a reference for any other analytical/computational model of FG plates.

Latin
, , , , , , 1, , , , , , .; , and thickness of the plate p volume fraction index u, v, w in-plane and transverse displacement at mid-plane of the plate u of the mid-plane in the x; y; z directions

Figure 1 :
Figure 1: Geometry of FG plate with positive set of reference axes.

0w
= ), and the solutions of the free vibration problems of the plates by removing in-plane loads ( With self-developed Matlab's code, various numerical examples are presented and discussed for verifying the accuracy and efficiency of the present theory in predicting the buckling and vibration responses of simply supported FG plates.The considered FG plate composesof aluminum (as metal) and alumina (as ceramic) with the following material properties: Al:Em = 70 GPa, νm = 0.3, ρm = 2702 kg/m 3 ; Al2O3: Ec = 380 GPa, νc= 0.3, ρc= 3800 kg/m 3 .

Figure 2 :
Figure 2: The variation of non-dimensional critical buckling load cr N ˆ of rectangular plate versus power law index p (b/a = 2, a/h = 10) with various shear deformation theories.

Figure 3 :
Figure 3: The variation of non-dimensional critical buckling load N ˆ of rectangular plate versus thickness ratio a/h (b/a = 2, p = 5) with various shear deformation theories.

Figure 4 :
Figure 4: Comparison of the variation of non-dimensional fundamental frequency ŵ of square plate versus power law index p.

Figure 5 :
Figure 5: Comparison of the variation of non-dimensional fundamental frequency ŵ of square plate versus thickness ratio a/h. :

Table 2 :
variation of the non-dimensional buckling load is considerable when p is small, and the fullyceramic plate gives the largest critical buckling load.An excellent agreement between the results predicted by present HSDT and full HSDT with 12 displacement unknowns for all values of power law indexis shown.Comparison of non-dimensional critical buckling load exhibits a variation of non-dimensional buckling loads cr N ˆ versus power law index p of rectangular FG plates (b/a = 2, a/h = 10)with various types of loading.It is observed that the non-dimensional critical buckling load decreases as p increases, Latin American Journal of Solids and Structures 13 (2016) 456-477 the cr N ˆof plates under different loading types with different values of power-law index p (b/a = 2, a/h = 10).Latin American Journal of Solids and Structures 13 (2016) 456-477

Table 3 :
Comparison of non-dimensional critical buckling load N ˆ of plates under different loading types with various values of side-to-thickness ratio a/h (b/a = 2, p = 5).

Table 4 :
Comparison of non-dimensional fundamental frequency ŵ of square plate.

Table 5 :
Non-dimensional frequency ŵ of plates with different values of power-law index p (b/a = 2, a/h = 10).

Table 6 :
Non-dimensional frequency ŵ of plates with different values of side-to-thickness a/h (b/a = 2, p = 5). D