Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 54, 1, pp. 135-147, Warsaw 2016 DOI: 10.15632/jtam-pl.54.1.135 THREE-DIMENSIONAL THERMAL BUCKLING ANALYSIS OF FUNCTIONALLY GRADED CYLINDRICAL PANELS USING DIFFERENTIAL QUADRATURE METHOD (DQM) Seyed A. Ahmadi, Hadi Pourshahsavari Babol University of Technology, Department of Mechanical Engineering, Babol, Iran e-mail: ali ahmadi1366@yahoo.com Thermal buckling analysis of functionally graded cylindrical panels subjected to various conditions is discussed in this paper. Buckling governing equations are solved using the differential quadrature method. It is assumed that the mechanical properties of the panel are graded through thickness according to a power function of the thickness variable. The panel is assumed to be under the action of three types of thermal loading including uniform temperature rise and variable temperature rise in the axial and radial direction. In the present study, the effects of power law index, panel angle, different thermal load conditions and geometric parameters on the buckling behavior of functionally graded curvedpanels are studied.The results obtained through the presentmethod are compared to the finite element solutions and the reported results in the literature. A desirable compatibility is concluded. Keywords: thermal buckling, curvedpanel, functionally gradedmaterial, differential quadra- ture method 1. Introduction Due to special mechanical properties, circular cylindrical panels are widely used in engineering structures such as pressure vessels, nuclear reactors, spacecrafts and jet engine exhausts. Due to the increasing demands for heat-resisting, energy absorbing, light-weight elements and high structural performance requirements in extremely high temperature environments and high- speed industries such as fusion reactors, aircraft and aerospace structures the use of special materials with high thermal and mechanical resistance has gained much popularity by many researchers. The applications of functionally graded materials (FGMs) have attracted much attention in the past two decades since they were first reported by Koizumi (1993). FGMs are composite materials, microscopically inhomogeneous, in which mechanical properties vary smoothly and continuously from one surface to the other. Themain advantage of FGMs is that the ceramic component provides high temperature resistance due to its low thermal conductivity while the metal component prevents fracture induced by thermal stresses due to the high- temperature gradient in a very short period of time. When these are subjected to a thermal loading, the determination of thermal buckling capacity of these structures is important to achieve an optimized design in cost and weight. Buckling analyses of various structureswere carried out bymany researchers. A review of re- search on the buckling response of plates and shells in a temperature environmentwas presented byThornton (1993). He did some research on thermal buckling of plates and shells. In his work, he described elastic thermal buckling of metallic as well as composite plates and shells. Murphy and Ferreira (2001) investigated thermal buckling analysis of imperfect flat plates based on the energy consideration. They showed the ratio of the critical temperature for a perfect rectangular plate to that of an imperfect plate as a function of the initial imperfection amplitude. Mahayni (1966) studied thermal buckling behavior of doubly curved isotropic panels using Galerkin’s 136 S.A. Ahmadi, H. Pourshahsavari method. Chang and Chui (1991) carried out bifurcation buckling analysis of composites under the action of uniform temperature change using higher order transverse shear deformation the- ory and the finite element method. Earlier, the Differential Quadrature Method introduced by Jang et al. (1989), was applied only to rectangular plates and lately it was considered for shells. Mirfakhraei and Redekop (1998) used the Differential Quadrature Method to study buckling behavior of circular cylindrical shells. Alibeigloo and Kani (2010) and Haftchenari et al. (2007) used this method to study cylindrical shells as well. The study of structures of functionally gradedmaterials has received considerable attention in recent years. Buckling of functionally graded plates under thermal loads was studied by Ja- vaheri and Eslami (2002b). They used classical plate theory and obtained nonlinear equilibrium and linear stability equations using variational formulations. Shahsiah and Eslami (2003) con- sidered effects of various temperature distributions on thermal buckling of simply supported FG cylindrical shells, using the first order shear deformation theory, however the temperature dependency of material properties was not included. Thermoelastic stability of FG cylindrical shells subjected to various thermal load conditions was studied by Wu et al. (2005). Thermal buckling analysis of functionally graded plates considering simply supported boundary condi- tions by using the first shear deformation theory was carried out by Wu (2004). He reached the stability equation of functionally graded shells usingDonnell’s shell theory and presented its closed-form solution.Buckling analysis ofFGplates usingahigher order theorywaspresentedby Javaheri and Eslami (2002a). It was shown that higher order shear deformation theory accura- tely predicts the buckling behavior, whereas the classical plate theory overestimates the critical loads. Breivik (1997) discussed the buckling response of composite cylindrical panels under the action of mechanical and thermal loading. Zhao et al. (2007) and Zhao and Liew (2010) used the element-free kp-Ritz method for thermal and mechanical buckling analysis of functionally graded cylindrical shells. They obtained three-dimensional buckling equations of the shell based on theDonnell shell theory andpresented a closed form solution to predict buckling loads caused by thermal loads and critical edge displacement in the longitudinal direction. In this paper, buckling analysis of cylindrical panels made of a functionally gradedmaterial subjected to three types of thermal loading is investigated. To obtain the buckling load of the cylindrical panels, the Differential QuadratureMethod (DQM) is used to discretize differential equations obtained based on the second Piola-Kirchhoff stress tensor using three-dimensional theory of elasticity byAkbari Alashti andAhmadi (2014). Thematerial properties are assumed to be temperature independent and vary continuously along the thickness according to a power law function while Poisson’s ratio of the material is taken to be constant. Effects of various parameters including panel curvature, grading index, various thermal load conditions and geo- metric ratios on the buckling behavior of the curved panels are investigated. Numerical results are validated against finite element calculations and results that are available in the offered literature. 2. Governing equation for buckling Consider a thick cylindrical panel made of ceramic and metallic materials with the inner ra- dius R1, mid-surface radius a, thickness h and length L. The geometric parameters and the cylindrical coordinate system. i.e. r, θ and x-coordinates are shown in Fig. 1. The components of the displacement field in this coordinate system are expressed as w, v andu, respectively. Assume that thematerial is isotropic, inhomogeneouswithYoung’smodulus varying continuously in the thickness direction, i.e. from ceramic in the inner layer to metallic in the outer layer according to the following formula Vm = (2z+h 2h )K Vc+Vm =1 (2.1) Three-dimensional thermal buckling analysis... 137 Fig. 1. Geometry of a cylindrical panel where Vc and Vm represent the volume fractions of the ceramic and metallic constituent and K denotes the volume fraction index that indicates thematerial variation profile through theFG shell thickness. Thus, the Young modulus in the radial direction is assumed to vary according to the power law in the following forms E(z) =Ec+Emc (2z+h 2h )K Emc =Em−Ec (2.2) whereEm andEc denote the elasticmodulusof themetal andceramic, respectively.Thematerial composition varies smoothly from the outer surface (z=h/2) of the shell as metal to the inner surface (z=−h/2) as ceramic.Material properties of the shell are assumed to be independent of the temperature field and Poisson’s ratio is considered to be constant throughout the thickness of the shell. In order to calculate buckling loads of panels, the buckling equations obtained by Akbari Alashti and Ahmadi (2014) are used. In this work, also the finite element linear or bifurcation buckling analysis of the cylindrical panel usingANSYS suite of program is carried out. The eigen buckling analysis predicts theore- tical buckling strength of a shellmade of a linear elasticmaterial. This analysis is used to predict the bifurcation point on an F-U diagram using a linearized model of the elastic structure. It is a technique used to determine buckling pressures at which the structure becomes unstable and their corresponding bucklingmode shapes. The basic form of the eigen buckling analysis is Kφ=λiSφ (2.3) where K, φi, λi and S are the structural stiffness matrix, eigenvector, eigenvalues and stress stiffness matrix, respectively. Eight noded quadrilateral shell elements, namely Shell281, are used to model the thick cy- lindrical shell. The elements can handlemembrane, bending and transverse shear effects and are able to form the curvilinear surface satisfactorily. The elements are suitable for modeling of the layer and have the stress stiffening, large deflection and large strain capabilities. Boundary conditions of shell panels are defined using equilibrium equations. For the initial and perturbed equilibrium positions, we have σrr ( a+ h 2 ,θ ) =σrr ( a− h 2 ,θ ) =0 τrθ ( a+ h 2 ,θ ) = τrθ ( a− h 2 ,θ ) =0 τrx ( a+ h 2 ,θ ) = τrx ( a− h 2 ,θ ) =0 (2.4) 138 S.A. Ahmadi, H. Pourshahsavari Boundary conditions at the panel edges are defined as: — up and down edge, x=0,L Simply supported: w= v=σ′xx =0 Clamped: w= v=u=0 (2.5) — lateral edges, θ=0,β Simply supported: w=σ′θθ =u=0 Clamped: w= v=u=0 (2.6) 3. Calculation of buckling load In this work, two types of panels are considered: Case 1. Thepanel is assumed tobe simply supportedat lateral edges and clamped at two ends. Therefore, thermal variation causes no axial stress on the panel,Nθ =0. Case 2. We assume that the panel has clamped boundary conditions at all edges. For this case thermal loading cause axial and circumferential stresses at the panel walls, Nx 6= 0, Nθ 6=0. Bysubstitutingthecomponents of thedisplacementfield in the stress-strainand linear strain- displacement equations and the resulted expression in the buckling equations, the equilibrium equations are defined in terms of components of the displacement field. In the present work, a polynomial expansion based on the Differential QuadratureMethod applied by Bellman and Casti (1971) is used to discretize and solve the obtained buckling equations. According to this method, the first order derivative of the function f(x) can be approximated as a linear sum of all functional values in the domain df dx ∣ ∣ ∣ ∣ ∣ x=xi = N ∑ j=1 w (1) ij f(xj) for i=1,2, . . . ,N (3.1) wherew (1) ij is theweighting coefficient andN denotes the number of grid pointsxi in thedomain. There are differentmethods for calculation of the weighting coefficients matrix, see Shu (2000). Here, the weighting coefficients of the first order derivatives are defined based on the Lagrange interpolation polynomials as w (1) ij = M(1)(xi) (xi−xj)M(1)(xj) for i 6= j w (1) ii = M(2)(xi) 2M(1)(xj) (3.2) where M(1)(xi)= N ∏ k=1 k 6=i,j (xi−xk) N(xi,xj)=M (1)(xi)δij M(2)(x)=N(2)(x,xk)(x−xk)+2N (1)(x,xk) (3.3) and for higher order derivatives, we have w (r) ij = r ( w (1) ij w (r−1) ii − w (r−1) ij xi−xj ) for i,j =1,2, . . . ,N r=2,3, . . . ,N−1 w(r)ii=− N ∑ j=1,j 6=i w(r)ij (3.4) Three-dimensional thermal buckling analysis... 139 Now, applying the above formulation to the buckling equations, we have G2 N ∑ l=1 a (2) i,l wl,j,k+ G2 r N ∑ l=1 a (1) i,l wl,j,k+G(z) Q ∑ n=1 b (2) j,nwi,n,k+G1 N ∑ l=1 Q ∑ n=1 a (1) i,l b (1) j,nul,n,k − G2 r2 wi,j,k+ G1 r N ∑ l=1 M ∑ m=1 a (1) i,l c (1) k,m vl,j,m+ G(z) r2 M ∑ m=1 c (2) k,m wi,j,m+ G3 r2 M ∑ m=1 c (1) k,m vi,j,m +σ0x Q ∑ n=1 b (2) j,nwi,n,k+σ 0 θθ 2 r2 M ∑ m=1 c (1) k,m vi,j,m+σ 0 θθ 1 r2 Bi,j,k−σ 0 θθ 1 r2 M ∑ m=1 c (2) k,m wi,j,m =0 G(z) N ∑ l=1 a (2) i,l vl,j,k+ G1 r N ∑ l=1 M ∑ m=1 c (1) k,m a (1) i,l wl,j,m+ G3 r2 M ∑ m=1 c (1) k,m wi,j,m+G(z) Q ∑ n=1 b (2) j,nvi,n,k + G(z) r N ∑ l=1 a (1) i,l vl,j,k− G(z) r2 vi,j,k+ G2 r2 M ∑ m=1 c (2) k,m vi,j,m+ G1 r Q ∑ n=1 M ∑ m=1 c (1) k,m b (1) j,nui,n,m +σ0x Q ∑ n=1 b (2) j,nvi,n,k−σ 0 θθ 2 r2 M ∑ m=1 c (1) k,m wi,j,m+σ 0 θθ 1 r2 vi,j,k−σ 0 θθ 1 r2 M ∑ m=1 c (2) k,m vi,j,m =0 G(z) N ∑ l=1 a (2) i,l ul,j,k+G2 Q ∑ n=1 b (2) j,nvi,n,k+G1 N ∑ l=1 Q ∑ n=1 b (1) j,na (1) i,l wl,n,k+ G(z) r N ∑ l=1 a (1) i,l ul,j,k + G(z) r2 M ∑ m=1 c (2) k,m ui,j,m+ G1 r M ∑ m=1 c (1) k,m vi,j,m+ G1 r Q ∑ n=1 b (1) j,nwi,n,k +σ0x Q ∑ n=1 b (2) j,nui,n,k−σ 0 θθ 1 r2 M ∑ m=1 c (2) k,m ui,j,m =0 (3.5) where G1 =G(z)+λ(z) G2 =2G(z)+λ(z) G3 =3G(z)+λ(z) and a (k) ij , b (k) ij and c (k) ij denote the weighting coefficients of the k-th order derivative in the r, θ and x-direction, respectively; N, Q and M are grid point numbers in the r, θ and x-direction, respectively. The critical value of the buckling load is obtained by solving the set of equations presented in the matrix form as [ BB BD DB DD ]      db u v w      =σ [ 0 0 DBG DDG ] (3.6) where the sub-matrices BB, BD and DBG, DD, DB, DDG are found from the boundary con- ditions and governing equations, respectively. Equation (3.6) is transformed into the standard eigenvalue equation, as ( −DBGB −1 B BD +DDG )−1( −DBB −1 B BD+DD ) [u v w]T−σI[u v w]T =0 (3.7) fromwhich, the eigenvalues ofσ can be found.The smallest value ofσ is found tobe thebuckling load. 140 S.A. Ahmadi, H. Pourshahsavari 4. Thermal loading 4.1. Uniform temperature rise The temperature changes uniformly through the thickness and remains constant in the longi- tudinal and circumferential directions of the panel. This thermal variations induces only normal stress, and the parameter Φ is defined as σ= N h N =− Φ 1−ν (4.1) and Φ= h 2 ∫ −h 2 [ Em+Ecm (2z+h 2h )K][ αm+αcm (2z+h 2h )K] ∆T(x,θ,z) dz ⇒ Φ= ( Ecαch+ [Ec(αm−αc)+αc(Em−Ec)]h K+1 + (αm−αc)(Em−Ec)h 2K+1 ) ∆Tcr (4.2) Substituting buckling stress obtained by numerical solution into Eq. (4.1) and (4.2), helps us to obtain the thermal buckling load ∆Tcr. 4.2. Non-uniform temperature rise in the axial direction In this case, the assumed temperature varies in the longitudinal direction according to the following formula T =∆T (x L )n +Tm ∆T =Tc−Tm n> 0 (4.3) whereTm is the temperatureat themetal surfaceof thepanels.According to theabove equations, axial stresses caused by the temperature rise have the same variation in this direction. The critical stresses are obtainedby considering the effects of this loading in thediscretized governing equations and then, the buckling temperatures are achieved using equations (4.1) and (4.2). 4.3. Non-uniform temperature rise in the radial direction The functionally graded materials are designed in order to resist against high temperature rise by ceramic, so the temperature change will be quite different at the two sides of FGM structures. The temperature distribution across the thickness is a function of the z coordinate as follows T =∆T ( − z h + 1 2 )q +Tm − h 2