Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 56, 4, pp. 1097-1108, Warsaw 2018 DOI: 10.15632/jtam-pl.56.4.1097 FINITE ELEMENT STUDY ON THERMAL BUCKLING OF FUNCTIONALLY GRADED PIEZOELECTRIC BEAMS CONSIDERING INVERSE EFFECTS Reza Nasirzadeh, Bashir Behjat, Mahsa Kharazi Mechanical Engineering Faculty, Sahand University of Technology, Tabriz, Iran e-mail: behjat@sut.ac.ir In this article, the bucklingbehavior andbifurcationpoint ofFunctionallyGradedPiezoelec- tric (FGP) beams are investigatedbased onEuler-Bernoulli beam theory. The finite element method is employed to model the beam in thermal environment. The material properties of the beam are considered to vary gradually in the thickness direction and the beam is subjected to electrical and thermal loading. In this paper, direct and inverse piezoelectric effects are considered andbuckling of the beam in the sensor state is investigated.By solving the eigenvalue problem, the buckling load of the FGP beam is obtained and the effect of various parameters such as power law index, temperature, applied voltage and beam aspect ratio on the buckling load are investigated. The results show that the boundary conditions are the main factor that affects the buckling load of the FGP beam. Keywords: FGPMbeam, buckling load, inverse piezoelectric effect, thermal environment 1. Introduction In the recent years, smart materials have found more applications as sensors and actuators in structures. Piezoelectric materials are the most significant materials that are used in structures for the purpose of monitoring and controlling as sensors and actuators. They are used in elec- tromechanical, medical and aerospace industries (Takagi et al., 2003). Also Functionally Graded Materials (FGMs) have received more attention in the recent decades because of their especial behavior in thermal environments. These materials can endure high thermal stresses and ther- mal shocks because of smooth distributions of thematerials through thickness. New functionally graded piezoelectric materials are introduced by combining the concept of so called piezoelectric materials and FGMs. There are many reports on related topics that have been presented in the last years and someof themarementionedhere.Kapuria andAlam(2004) introducedan efficient coupled one-dimensional nonlinear zigzag theory for buckling analysis of hybrid piezoelectric be- ams under electromechanical loading. They used an analytical solution for buckling of simply supported beams. Li et al. (2006) investigated thermal post-buckling of FGM Timoshenko be- ams subjected to a transversely non-uniform temperature rise. They used the shooting method to solve the problem, and the thermal buckling and post-buckling response of FGMTimoshenko beamswithfixed-fixededgeswere obtained. JeromeandGanesan (2010) usedageneralizedplane strain finite element formulation for the buckling analysis of piezocomposite beams. They used a 2D finite element formulation to improve accuracy in prediction of the buckling load of the piezocomposite. Kiani and Eslami (2010) investigated the buckling load of functionally graded Euler-Bernoulli beams in thermal environments. In thiswork, the effect of three types of thermal loading was investigated on the critical buckling load of FGM beams.Wattanasakulpong et al. (2011) investigated thermal buckling of FGM beams using third-order shear deformation beam theory. The Ritz method was adopted to obtain the eigenvalue problem of thermal buckling in various types of immovable boundary conditions. Kiani et al. (2011) investigated thermo- -electrical buckling of piezoelectric functionally graded beams based onTimoshenko theory. The 1098 R. Nasirzadeh et al. electric field in the piezoelectric layer was assumed constant through the thickness. The results were obtained for three types of thermal loading. Fallah andAghdam (2011) suggested a simple analytical expression for large amplitude free vibration and a post-buckling analysis of functio- nally graded beams resting on elastic foundation. Euler-Bernoulli assumptionswith vonKarman type nonlinear strain-displacement relations were used to derive the governing equation of mo- tion.Komijani et al. (2012) investigated non-linear thermo-electrical stability of anFGPMbeam using Timoshenko beam theory. All of thermo-electro-mechanical properties were assumed to vary in the thickness direction and expressedbyapower lawdistribution.TheRitz finite element method was used to solve the governing equation. Fu et al. (2012) analyzed buckling and free vibration of the FGM beam with two clamped ends and surface-bonded piezoelectric actuators in thermal environments. The governing equation was derived by Hamilton’s principle. Rahimi et al. (2013) studied the free vibration and post-buckling behavior of functionally graded be- ams using Timoshenko beam theory. They investigated the post-buckling behavior of the FGM beams and obtained an exact solution based on the applied axial load. Li and Batra (2013) achieved an analytical relations between the critical buckling load of a FGMTimoshenko beam and the corresponding homogeneous Euler-Bernoulli beam subjected to axial and compressive load for various boundary conditions. Rafiee et al. (2013) investigated nonlinear thermal buc- kling of carbon nanotube reinforced composite beams with surface-bonded piezoelectric layers. The governing equations of the beamwere obtained based on the Euler-Bernoulli beam theory and considering von Karman type geometric nonlinearity. The critical temperature at which buckling occurred was obtained in that article. Esfahani et al. (2013) studied thermal buckling and post-buckling analysis of FGMTimoshenko beams resting on the non-linear elastic founda- tion. Different types of boundary conditions such as clamped, simply supported and roller edges were considered, and the generalized differential quadrature method was employed to solve the problem.Komijani et al. (2014) analyzed buckling, post-buckling and small amplitudevibrations of functionally graded beams resting on a nonlinear elastic foundation and subjected to in-plane thermal loads. Themicrostructural length scale based on themodified couple stress theory was used to derive the governing nonlinear equilibrium equations of the FGM beam. Nasirzadeh et al. (2014) used an exact solution to investigate the stability of FGP beams based on the Euler- -Bernoulli beamtheory.Ghiasian et al. (2015) studieddynamicbucklingof theFGMTimoshenko beam subjected to sudden uniform temperature rise considering imperfection. The analysis was performedwith an assumption of temperature dependency of each thermo-mechanical property of the FGM beam. The obtained non-linear algebraic equations were solved via the Newton- -Raphson iterative scheme. Li and Qiao (2015) studied buckling and post-buckling behavior of shear deformable anisotropic laminated composite beams with initial imperfection subjected to axial compression. The governing equations were based on the higher order shear deformation beam theory with a von Karmann type nonlinearity. In that paper, composite beams with the fixed-fixed, fixed-hinged, and hinged-hinged boundary conditions were considered. The results were obtained by combining Newton’s iterative method and Galerkin’s method. Chen et al. (2015) investigated the elastic buckling and static bending of shear deformable functionally gra- ded (FG) porous beams based on theTimoshenko beam theory. Thepartial differential equation was derived based on Hamilton’s principle, and the Ritz method was employed to obtain the critical buckling load of porous beams. In this paper, the buckling analysis of FGP beam is investigated considering the inverse piezoelectric effect. The finite element method is employed to model the beam. The effect of two different types of thermal loading on the buckling load of the beam is investigated. The main aim of the study is to evaluate the influence of functional grading of the properties on the bucklingbehavior of thebeams in thermal environments.Various numerical results are presented in graphical forms to give an insight into the influence of material composition, loading type and boundary condition on the buckling load and bifurcation point of the FGP beam. Finite element study on thermal buckling of functionally graded... 1099 2. Theoretical formulation 2.1. Geometry and material definition The beam has length l, height and width h and b, respectively, which is shown in Fig. 1. The coordinate system of the FGPM beam oriented in such a way that the x-axis is in the longitudinal direction in the middle of the beam and the z-axis and y-axis are in the thickness of the beam up to the plate and perpendicular to it, respectively (Fig. 1). Fig. 1. Definition of geometry and coordinate system of the FGP beam With regard to gradual changes in thematerial properties through the thickness of the beam, on the basis of the power distribution law, the effective properties of material can be defined as follows P(z) = PL+PUL (1 2 + z h ) (2.1) where PUL = PU −PL and PL, PU are the properties in the bottom and upper surface of the FGP beam, respectively. 2.2. Kinematic and constitutive relations In this study, considering Euler-Bernoulli beam assumptions, the displacement field is pre- sented as follows (Bathe, 1996) û(x,z) = u(x)−zw,x ŵ(x,z) = w(x) (2.2) where ŵ(x), û(x) are beam displacements in the z and x coordinates in a general point, and w(x), u(x) show the displacement of the beam in mid-plane, and w,x is the beam slope. The displacement field based on Euler-Bernoulli beam theory allows only the axial component of the strain component εx and the other components to be zero. The strain-displacement relations considering non-linear von Karman strain can be defined for the normal strain as εx = u,x+ 1 2 (w,x) 2−zw,xx (2.3) The characteristic equation of the FGP beamunder thermal, electrical andmechanical forces in matrix form is expressed as follows σ=Qε−eTE−Qα∆θ D= eε+kE+P∆θ (2.4) in which D, E, σ represent vectors of electric displacement, electric field and stress, respecti- vely. Also the matrices P, k, e, Q and α represents the pyroelectric vector, dielectric matrix, piezoelectric matrix, elastic stiffness matrix and thermal expansion vector. 1100 R. Nasirzadeh et al. 2.3. Thermal and electrical loading In this paper, the effect of two different types of thermal loading (uniform temperature rise and linear temperature rise) on thebuckling loadof thebeamis studied. In thefirst case, the total volume of the beam is subjected to the temperature rise of ∆θ = θ−θ0 where θ0 is the reference temperature. In the second case, the upper layer of the beam is subjected to temperature rise but the lower surface remains in the reference temperature. In this case, by using the narrow beam assumptions, the heat transfer equation can be solved as follows (Kiani et al., 2011) θ = θL+θUL (1 2 + z h ) θUL = θU −θL ∆θ = θ−θ0 (2.5) where θL and θU are the lower and upper surface temperature of the beam, respectively. Considering both the direct and reverse piezoelectric effect, the electrical potential function can be assumed as (Komijani et al., 2013) V (x,z)= cos(βz)ϕ(x)+ V0z h (2.6) here V0 is the applied electrical potential to the beam,which is a constant value, andϕ(x) shows variation of the electric potential in theaxial direction andβ = π/h.Theelectric field isdescribed as follows Ex =−cos(βz)ϕ,x Ey =0 Ez = β sin(βz)ϕ(x)−E0 (2.7) where E0 = V0 h (2.8) 2.4. Governing equations Thegoverning equation of the beamconsidering the reverse effect undermechanical, thermal and electrical loads will be derived by using the principle of minimumpotential. Accordingly, it can be written as δΠ = δU + δWext =0 (2.9) where U is the total potential energy of the beam and Wext is the work of external forces. The total potential energy of the piezoelectric beam can be defined as follows δU = ∫∫∫ V [δεTσ− δETD− δθS] dV (2.10) It is obvious that δθ =0 (2.11) The electric field variation by considering Eq. (2.7) can be obtained as δEx =−cos(βz)δϕ,x δEy =0 δEz = β sin(βz)δϕ (2.12) The linear part of the strain and the nonlinear part of it can be written as ε= εL δε= δεL+ δεNL (2.13) Finite element study on thermal buckling of functionally graded... 1101 where δεL = [(δu,x−zδw,xx),0,0,0,0,0] T δεNL = [w,xδw,x,0,0,0,0,0] T (2.14) By replacing relations (2.11) to (2.14) in relation (2.10), the total potential energy can be divided into two types of elastic potential energy and geometrical potential energy in tension. These equations can be written as δUela = ∫∫∫ V [δεTL(QεL−e T E−Qα∆θ)] dV + ∫∫∫ V [−δET(eεL+kE+P∆θ)] dV δUgeo = ∫∫∫ V [δεTNL(QεL−e T E−Qα∆θ)] dV (2.15) The work done by external forces acting on the beam can be written as follows δWex =−δu T ( fp+ ∫ A fS dA ) − δϕT ∫ A qA dA (2.16) In which qA, fS, fp represents charge density, surface forces and concentrated load, respectively. By replacing relations (2.15) and (2.16) in the principle ofminimumpotential energy (Eq. (2.9)), the total potential energy will be achieved as δ∓= ∫∫∫ V [δεTL(QεL−e T E−Qα∆θ)] dV + ∫∫∫ V [−δET(eεL+kE+P∆θ)] dV + ∫∫∫ V [δεTNL(QεL−e T E−Qα∆θ)] dV − δuT ( fp+ ∫ A fS dA ) − δϕT ∫ A qA dA (2.17) 2.5. Finite element formulation To investigate stability behavior of the FGPbeam, a two-node beam element is used (Fig. 2) tomodel the beam. Each node in the element has 4 degrees of freedom as ϕ, w0,x, w0, u0 which represent the electric potential, slope, deflection and axial displacement of the beam and as nodal variables. Fig. 2. Definition of the 2-node beam element Dimensionless variables instead of the coordinate variable are used as follows ξ = x a η = z b (2.18) By using the Hermite and Lagrange interpolation functions, the axial displacement, vertical deflection and electrical potential in the beam can be written as u0 = 2∑ i=1 ψiu e 0i w0 = 4∑ i=1 ψiΦ e i ϕ0 = 2∑ i=1 ψiϕ e 0i (2.19) 1102 R. Nasirzadeh et al. inwhichψi andψidenote theLagrangeandHermite interpolation functions, respectively (Bathe, 1996). Based on the above relationships, the linear and nonlinear strain vector and the electric field vector are expressed based on nodal variables as follows εL =BLu e ε ′ NL =BNLu e E=Bϕϕ e+E0 (2.20) where E0 =    0 0 −E0    E=    − cos(βz)ϕ,x 0 β sin(βz)ϕ(x)    +    0 0 −E0    BL =    ψ1,x zψ1,xx zψ2,xx ψ2,x zψ3,xx zψ4,xx 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0    BNL =    0 ψ1,x ψ2,x 0 ψ3,x ψ4,x 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0    Bϕ =    −cos(βz)ψ1,x −cos(βz)ψ2,x 0 0 β sin(βz)ψ1 β sin(βz)ψ2    (2.21) By using relation (2.20) and replacing it in equation (2.17), the final relation of the FGP beam inmatrix form can be obtained as (Kuu−σ0KGuu)u e t −Kuϕϕ e t =Fm+Fuθ Kϕuu e t +Kϕϕϕ e t =Fϕ+Fϕθ (2.22) in which σ0 represents the negative pre-stress in the beam andKuϕ,KGuu,Kuu andKϕϕ show the piezoelectric stiffness, geometric stiffness, elastic stiffness matrix and dielectric permittivity matrix.Fuθ,Fϕ,Fm andFϕθ are the thermal expansion, electrical, mechanical andpyro-electric load vectors in the local coordinate system. They can be defined as below Kuu = ∫∫∫ V B T LQBL dV = b ∑ element 1∫ −1 1∫ −1 B T LQBL|J| dξ dη KGuu = ∫∫∫ V B T NLσ0BNL dV = b ∑ element 1∫ −1 1∫ −1 B T NLσ0BNL|J| dξ dη Kϕu = ∫∫∫ V B T ϕeBL];dV = b ∑ element 1∫ −1 1∫ −1 B T ϕeBL|J| dξ dη Finite element study on thermal buckling of functionally graded... 1103 Kϕϕ = ∫∫∫ V B T ϕkBϕ dV = b ∑ element 1∫ −1 1∫ −1 B T ϕkBϕ|J| dξ dη (2.23) Fm = N∑ i=1 N T uH Tfpi+ ∫ A N T uH TfS dA Fuθ = ∫∫∫ V B T LQα∆θ dV Fϕ = ∫ A N T ϕqA dA Fϕθ = ∫∫∫ V −BTϕP∆θ dV Fϕϕ = ∫∫∫ V −BTϕkE0 dV By simplifying relation (2.21)3 based on the electrical potential and combiningwithEq. (2.21)2, the final equation of the FGP beamwill be achieved as (Kuu+KuϕK −1 ϕϕKϕu−KGuu)u e t =Fm+Fuθ+KuϕK −1 ϕϕ(Fϕ+Fϕθ) (2.24) The above equation expresses the final governing relation of the FGP beam using the finite elementmodel based on the Euler-Bernoulli theory under thermo-electro-mechanical loading. It is obvious that based on homogenous and nonhomogeneous equations, the behavior of the beam will be changed and this is affected by the boundary condition of the beam. It means that if the boundary condition is in such a way that the problem change to the eigen-value problem, the beam will buckle, which can be found by solving the eigenvalue problem. Otherwise, if Eq. (2.24) has a non-zero value on the right hand side of the relation, it will be an ordinary equation and represent only the bending problem. 3. Results and discussion 3.1. Comparison studies In this Section, to verify the solution procedure, the results obtained in this paper are compared with data reported in the literature. The buckling load of the FGP beam is listed in Table 1 and the results are comparedwith the values reported byKomijani et al. (2013). The boundary condition of the beam is clamped in both ends and there is no electrical and thermal loading on the beam. The eigenvalue problem of Eq. (2.24) is solved and the results are shown in Table 1. As it is seen from this table, the values of the buckling load of the FGP beam are in good agreement with the data reported by Komijani et al. (2013). Table 1. Buckling load of the FGP beam for various power law indexes, L/h = 60, ∆θ = 0, V0 =0 Power law index 0 0.2 2 10 Komijani et al. (2013) 6.892e+4 7.171e+4 7.919e+4 8.317e+4 Present study 6.603e+4 7.067e+4 8.010e+4 8.319e+4 3.2. Parametric studies In this Section, the stability and buckling load of FGPbeamswill be investigated. The effect of parameters such as temperature and electrical fields, power law index and aspect ratio of the beamwill be discussed. 1104 R. Nasirzadeh et al. Figure 3 depicts the buckling load of the beam considering the direct and revers piezoelectric effects based on different power law indexes in various aspect ratios. It is seen that the buckling load increases with a decrease in the aspect ratio. On the other hand, the critical axial load increases versus the power law index. This is explained as that the beam is stiffer for higher power law indexes considering the elastic modulus. Also the sensor-actuator state has a greater buckling load than the actuator state. This can be explained by considering the sensor stiffness terms in Eq. (2.21)1. Fig. 3. The critical buckling load of the FGP beam in sensor-actuator and sensor states versus the power law index, ∆θ =0, V0 =0(v) Fig. 4. Buckling load of the beam versus the power law index for various constant and linear temperature rise, L/h=25, V0 =500(v) Finite element study on thermal buckling of functionally graded... 1105 Figure 4 compares the effect of uniform and linear temperature fields on the buckling load versus the power law index. The beamhas a clamped-clamped boundary condition and is under electrical loading (L/h =25, V0 =500(v). This figure shows that by increasing the temperature, the buckling load rises. Also, by increasing the power law index, the buckling load of the beam increases. This is explained by the difference in the elastic modulus and thermal expansion coef- ficient of the beammaterials, herePZT-4 andPZT-5H. It isworth to note that the buckling load in thermal environment depends on the aspect ratio, elastic and thermal expansion coefficients. In the next figures, the parameters kUT and kLT mean the power law index for uniform temperature rise, linear temperature rise, and k means the power law index under electrical loading. Figure 5 compares the effect of uniformand linear temperature changes on the buckling load of the clamped-clamped beam. The beam has a constant aspect ratio and is under electrical loading (L/h =25, V0 =500(v)). This figure shows that by increasing the power law index, the buckling load decreases. Also, for a constant power law index, the uniform temperature rise has a greater effect on the buckling load than the linear one. This is explained by considering the temperature distribution through thickness of the beam. Fig. 5. The effect of uniform and linear temperature changes on the buckling load of the beam Figure 6 depicts the buckling load of the beam versus the power law index under electrical loading. The beam has an aspect ratio L/h = 25 and is under uniform and linear thermal loading. It is inferred from this figure that the applied voltage has not a considerable effect on the buckling load. However, these results show that by increasing the power law index, the buckling load decreases because of a lower elastic modulus. Figure 7 shows the deflection of the mid-plane of the beam for various power law indexes and different aspect ratios. The behavior of the beam is very similar to the buckling analysis, and the axial force in the bifurcation point increases with an increase in the power law index and decreases with an increase in the aspect ratio. Figures 8 and 9 show the bifurcation point of the beam for simply-clamped boundary condi- tions. These figures show that the bifurcation point for the simply-clamped boundary conditions takes a larger axial load that the simply-roller one. 1106 R. Nasirzadeh et al. Fig. 6. The buckling load of the beam versus power law index for various electrical loadings Fig. 7. Deflection of the FGP beam versus axial load for aspect ratios and simply-roller boundary conditions, ∆θ =0, V0 =0(v), k =0.5 Fig. 8. Deflection of the FGP beam versus axial load for various power law indexes and simply-clamped boundary conditions, ∆θ =0, V0 =0(v), L/h=25 Finite element study on thermal buckling of functionally graded... 1107 Fig. 9. Deflection of the FGP beam versus axial load for aspect ratios and simply-clamped boundary conditions, ∆θ =0, V0 =0(v), k =0.5 4. Conclusions In this paper, buckling and bifurcation point of the Functionally Graded Piezoelectric (FGP) beams are investigated based on the Euler-Bernoulli beam theory in a thermal environment by using the finite elementmethod. All mechanical and electrical properties of the beam are consi- dered to change gradually in the thickness direction of the beam, and the inverse piezoelectric behavior of the beam is investigated. In addition, the effects of various parameters of the beam such as the power law index, temperature, applied voltage and the aspect ratio on the buckling andbifurcation point of the beamare investigated. The results show that because of the coupled electro-mechanical nature of the FGP beams, some boundary conditions have stable behavior while some others are unstable in thermal environments. Also, considering the inverse effect of the piezoelectricity increases the buckling load of the FGP beam. References 1. Bathe K.J., 1996,Finite Element Procedures, Prentice Hall, NewYork 2. Chen D., Yang J., Kitipornchai S., 2015, Elastic buckling and static bending of shear defor- mable functionally graded porous beam,Composite Structures, 133, 54-61 3. Esfahani S.E., Kiani Y., Eslami M.R., 2013,Non-linear thermal stability analysis of tempera- ture dependent FGM beams supported on non-linear hardening elastic foundations, International Journal of Mechanical Sciences, 69, 10-20 4. Fallah A., Aghdam M.M., 2011, Nonlinear free vibration and post-buckling analysis of functio- nally graded beams on nonlinear elastic foundation,European Journal of Mechanics-A/Solids, 30, 571-583 5. FuY.,Wang J.,MaoY., 2012,Nonlinear analysis of buckling, free vibration and dynamic stabi- lity for the piezoelectric functionally graded beams in thermal environment,Applied Mathematical Modelling, 36, 4324-4340 6. Ghiasian S.E., Kiani Y., Eslami M.R., 2015, Nonlinear thermal dynamic buckling of FGM beams,European Journal of Mechanics A/Solids, 54, 232-242 7. Jerome R., Ganesan N., 2010, New generalized plane strain FE formulation for the buckling analysis of piezocomposite beam, Finite Elements in Analysis and Design, 46, 896-904 1108 R. Nasirzadeh et al. 8. Kapuria S., Alam N., 2004, Zigzag theory for buckling of hybrid piezoelectric beams under electromechanical loads, International Journal of Mechanical Sciences, 46, 1-25 9. Kiani Y., Eslami M., 2010, Thermal buckling analysis of functionally graded material beams, International Journal of Mechanics and Materials in Design, 6, 229-238 10. Kiani Y., Rezaei M., Taheri S., Eslami M., 2011, Thermo-electrical buckling of piezoelectric functionally gradedmaterialTimoshenkobeams, International Journal ofMechanics andMaterials in Design, 7, 185-197 11. Komijani M., Esfahani S.E., Reddy J.N., Liu Y.P., Eslami M.R., 2014, Nonlinear thermal stability and vibration of pre/post-buckled temperature- andmicrostructure-dependent functional- ly graded beams resting on elastic foundation,Composite Structures, 112, 292-307 12. KomijaniM., Kiani Y., EslamiM., 2012,Non-linear thermoelectrical stability analysis of func- tionally graded piezoelectric material beams, Journal of Intelligent Material Systems and Structu- res, 24, 399-410 13. Komijani M., Reddy J.N., Ferreira A.J.M., 2013, Nonlinear stability and vibration of pre/post-buckledmicrostructure-dependent FGPMactuators,Meccanica, 49, 1-17 14. Li S.R., Batra R.C., 2013, Relations between buckling loads of functionally gradedTimoshenko and homogeneous Euler-Bernoulli beams,Composite Structures, 95, 5-9 15. Li S,-R., ZhangJ.-H., ZhaoY.-G., 2006,Thermal post-buckling of functionally gradedmaterial Timoshenko beams,Applied Mathematics and Mechanics (English Edition), 27, 803-810 16. Li Z.-M., Qiao P., 2015, Buckling and postbuckling behavior of shear deformable anisotropic laminated beams with initial geometric imperfections subjected to axial compression,Engineering Structures, 85, 277-292 17. Nasirzadeh R., Behjat B., Kharazi M., 2014, Stability of FGP beams under thermo- electromechanical loading, International Journal of Material Science Innovations, 2, 164-177 18. Rafiee M., Yang J., Kitipornchai S., 2013, Thermal bifurcation buckling of piezoelectric carbon nanotube reinforced composite beams,Computers and Mathematics with Applications, 66, 1147-1160 19. Rahimi G., Gazor M., Hemmatnezhad M., Toorani H., 2013, On the postbuckling and free vibrations of FGTimoshenko beams,Composite Structures, 95, 247-253 20. TakagiK., Li J.F.,YokoyamaS.,WatanabeR., 2003,Fabricationand evaluation ofPZT/PT piezoelectric composites and functionally graded actuators, Journal of the European Ceramic So- ciety, 23, 1577-1583 21. WattanasakulpongN., GangadharaPrustyB., KellyD.W., 2011,Thermal buckling and elastic vibration of third-order shear deformable functionally graded beams, International Journal of Mechanical Sciences, 53, 734-743 Manuscript received June 13, 2017; accepted for print April 7, 2018