Document Type : Research Paper
Authors
Department of Mechanical Engineering, Mashhad Branch, Islamic Azad University, Mashhad, Iran
Abstract
Keywords

Mechanics of Advanced Composite Structures 7 (2020) 297 – 311


Semnan University 
Mechanics of Advanced Composite Structures journal homepage: http://MACS.journals.semnan.ac.ir 
Buckling Analysis of Functionally Graded Cylindrical Shells Under Mechanical and Thermal Loads by Dynamic Relaxation Method
M.E. Golmakani^{*}, M. Moravej, M. Sadeghian
Department of Mechanical Engineering, Mashhad Branch, Islamic Azad University, Mashhad, Iran
Keywords: 

ABSTRACT 
Buckling FG shell Thermal Gradient DR method 
In this study, using the dynamic relaxation method, nonlinear mechanical and thermal buckling behaviors of functionally graded cylindrical shells were studied based on firstorder shear deformation theory (FSDT). It was assumed that material properties of the constituent components of the FG shell vary continuously along the thickness direction based on simple powerlaw and MoriTanaka distribution methods separately. An axial compressive load and thermal gradient were applied to the shell incrementally so that in each load step the incremental form of governing equations were solved by the DR method combined with the finite difference (FD) discretization method to obtain the critical buckling load. After convergence of the code in the first increment, the latter load step was added to the former so that the program could be repeated again. Afterwards, the critical buckling load was achieved from the mechanical/ thermal loaddisplacement curves. In order to validate the present method, the results were compared with other papers and the Abaqus finite element results. Finally, the effects of different boundary conditions, grading index, rule of mixture, radiustothickness ratio and lengthtoradius ratio were investigated on the mechanical and thermal buckling loads. 
Functionally graded materials (FGMs) are a kind of unique composite material typically made from a mixture of ceramic and metal. The changes in these components result in an inhomogeneous microstructure which leads to gradual variations in the macroscopic properties of the material. Through this special property, the ceramic component can boost the thermal resistance while the metallic composition can enhance the fracture toughness. Recently, FG cylindrical shells have been used in various industries such as: aerospace, nuclear and medical applications [1]. Despite this evident importance of circular cylindrical FG shells, there have been few investigations on the buckling behavior of these structures in comparison with plate structures or other types of shells. In this regard, the dynamic thermoelastic response of FG cylinders and plates was presented by Reddy and Chin [2]. Also, Li and Batra analyzed the buckling behavior of axially compressed simply supported thin circular cylindrical shells with an FG middle layer [3]. Moreover, buckling behaviors of FG cylindrical shells subjected to pure bending load were investigated by Huang et al. [4]. In another study, Shariyat studied the dynamic buckling of a prestressed, suddenlyheated imperfect FGM cylindrical shell and dynamic buckling of a mechanicallyloaded imperfect FGM cylindrical shells in a thermal environment, with temperature dependent properties [5]. Additionally, based on the nonlinear large deflection theory of cylindrical shells as well as the Donnell assumptions, Huang and Han presented nonlinear buckling and post buckling analysis of axially compressed functionally graded cylindrical shells using the Ritz energy method [6]. They also considered the buckling behaviors of axially compressed functionally graded cylindrical shells with geometrical imperfections utilizing the Donnell shell theory and the nonlinear straindisplacement relations of large deformations [7]. In further publications by Shahsiah and Eslami, studies were carried out on the buckling analysis of functionally graded cylindrical shells under two types of thermal loads with simply supported boundary conditions based on the firstorder shell theory (FSDT) and the Sanders kinematic equations [8]. A formulation for the free vibration and buckling of FG cylindrical shells subjected to combined static and periodic axial loadings were presented by Ebrahimi and Sepiani based on FSDT and the classical shell theory (CST) [9]. Furthermore, Shen studied the postbuckling analysis for an FG thin cylindrical shell of finite length subjected to external pressure and thermal environments [10]. Shen and Noda analyzed the postbuckling of shear deformable FG cylindrical shells of finite length subjected to combined axial and radial mechanical loads in thermal environments [11]. Based on Galerkin’s method, Mirzavand and Eslami studied the buckling analysis of imperfect FG cylindrical shells under axial compression in thermal environments based on classical shell theory and the Sanders nonlinear kinematic relations [12]. Nonlinear response of imperfect eccentricallystiffened FGM thin circular cylindrical shells surrounded by elastic foundations and subjected to axial compression was presented by Duc and Thang [13]. Khazaeinejad and Najafizadeh considered analytical solutions of the buckling behavior of FG cylindrical shells subjected to three types of mechanical loads using the FSDT [14]. In addition, some studies of nonlinear bending of circular/annular FG plates/disks subjected to mechanical or thermomechanical loadings based on the DR method with FD technique were conducted by Golmakani et al., solving the nonincremental form of governing equations [1519]. Based on a similar method and similar loadings, they also investigated the large deflection behavior of stiffened/ unstiffened FG sector plates [20, 21] and general theta ply laminated plates [22, 23]. Eigenvalue buckling of a multilayered FG cylindrical shells reinforced with graphene sheets was studied by Wang et al. by using finite element method (FEM) [24]. Also, Wang et al. analyzed the torsional buckling of FG cylindrical shells reinforced with GPLs by using FEM. They observed that the enhancement of the number of layers resulted in notable decline of stress gradient between neighboring layers [25]. Yiwen et al. studied thermal buckling of functionally graded orthotropic cylindrical shells analytically using Reissner's shell theory [26]. Trabelsi et al., investigated the thermal buckling of functionally graded plates and cylindrical shells using the modified First Order Shear deformation theory [27]. The nonlinear analytical torsional/ thermal buckling and postbuckling of multilayer FGM cylindrical shells were studied by Nam et al. [28]. Also, thermomechanical buckling and postbuckling of cylindrical shell with FG coatings were analyzed by Thang et al. They considered the classical shell theory based on the von Kármán assumptions, Galerkin’s method and Airy stress function to obtain the closedform solution [29]. Golmakani et al. studied buckling analysis of moderately thick FG cylindrical panels subjected to axial compression in different boundary conditions [30]. RezaieePajand et al. investigated thermomechanical buckling and postbuckling of functionally graded shells based on FSDT and Voigt’s model by finite element method (FEM) [31]. Wang et al., studied the buckling of graphene platelets (GPL) reinforced composite cylindrical shells with cutouts using the finite element method (FEM). In their article, they modified Halpin–Tsai micromechanics model and determined Young’s modulus of the composites, also using the rule of mixture, approximated the mass density and Poisson’s ratio of the composites [32]. Zghal et al. studied linear static analysis of FG carbon nanotubereinforced plate and shell structures [33]. Trabelsi et al., studied thermal postbuckling analysis of functionally graded material structures based on a modified FSDT via the finite element method and a large displacement was described by Green–Lagrange nonlinear strains [34]. Nonlinear thermal buckling of imperfect cylindrical shells using a continuumbased semianalytical finite element formulation was used by Alijani et al. [35]. Zghal et al. studied free vibration of functionally graded composite shell structures reinforced by carbon nanotubes based on the discrete double directors shell finite element formulation [36]. Zghal investigated the mechanical buckling of functionally graded materials and carbon nanotubesreinforced composite plates and curved panels based on a double director’s finite element shell model [37]. Furtermore, the dynamic analysis and forced vibration of functionally graded carbon nanotubesreinforced composite shell structures (FGCNTRC) based on a linear discrete double director’s finite element model were studied by Frikha et al. [38]. They also, investigated nonlinear deflections analysis of thin FGCNTRC shell structures based on a discrete form of Kirchhoff finite element model and the displacement field was approximated by four nodes and three node finite elements [39]. Nonlinear bending of nanocomposites reinforced by graphenenanotubes with finite shell element and membrane enhancement was studied by Zghal et al. based on the Thirdorder shear deformation theory [40].
But despite significant contributions to investigation of buckling behavior of cylindrical shells in the previous years, to the best of the authors’ knowledge, up to now the nonlinear mechanical and thermal buckling of FG cylindrical shells with various boundary conditions have not been studied based on FSDT. Hence, the present paper is concerned with the further development of the mechanical and thermal buckling analysis of FG cylindrical shells for clamped and simply supported boundary conditions based on FSDT and large deflection von Kármán equations. The DR method combined with the finite difference (FD) discretization technique is employed to solve these incremental formulations. However, the buckling behavior of FG cylindrical shells is not considered by the DR method, so far. In this paper, the critical buckling load is predicted based on mechanical/thermal load–displacement curve obtained by solving the incremental form of nonlinear equilibrium equations. In order to accurately predict the elastic properties of actual FGM’s, the Mori–Tanaka scheme is applied. Furthermore, the results of this theory are compared with the powerlaw distribution (simple rule of mixture) which indicates a significant difference. Also, unlike previous investigations, the nonlinear temperature distribution is considered along the thickness direction for the purpose of thermal buckling analysis. The results are compared with some references and those obtained by the Abaqus finite element software. Finally, numerical results for critical buckling load and critical temperature difference are presented for various boundary conditions, two different rules of mixture, grading indices, radius to thickness and lengthtoradius ratios.
Figure 1 shows an FG cylindrical shell with radius R, thickness h, and length L in the cylindrical coordinate system (x, θ, z).
The FG shell is considered a mixture of ceramics and metals with the material properties of the composition varying continuously and smoothly through the thickness of the shell. As mentioned above, there are different models which show the variation in the mechanical and thermal properties of the FGMs. The powerlaw distribution of the volume fraction, is the most common type. Based on the powerlaw model, the material property (the effective values of Young’s modulus , heat conductivity coefficient and the coefficient of thermal expansion through the thickness of the shell can be expressed as [41]:
(1) 
where subscripts and indicate the metallic and ceramic constituents, respectively; and are the ceramic and metal volume fractions, respectively, and follow as [41]:

(2) 
(3) 
where z is the thickness coordinate ( ) and grading index k dictates the material variation profile through the shell thickness. Since the prediction of the macroscopic stressstrain response of FGMs is related to the description of their complex microstructural behavior represented by the interaction between the constituents, using a micromechanicsbased method such as MoriTanaka’s theory; a selfconsistent scheme can represent the realistic prediction for the behavior of the FGMs [33]. The effective values of bulk modulus, , shear modulus, , thermal conductivity coefficient, and thermal Expansion coefficient, , of the functionally gradient material based on the Mori–Tanaka homogenization method are [42 44]:

(4) 

(5) 
(6) 


(7) 
Where:

(8) 
Thus, the effective values of and can be computed as follows:

(9) 

(10) 
The temperature variation is assumed to occur only in the direction of thickness. For the mentioned onedimensional temperature field, the study is designed so that the outer ceramic surface is exposed to higher temperatures compared to the inner metal surface. In this work, two types of linear and nonlinear temperature distributions are considered for thermal buckling analysis of FG cylindrical shell.
Fig. 1. FG cylindrical shell in the cylindrical coordinate system
In some works, in order to obtain the thermal buckling load, the linear temperature variation is assumed to be along the thickness direction. According to the linear temperature distribution, the following linear function is considered for thermal load distribution along the thickness [8, 45]:

(11) 
where at and at .
In this case, the temperature distribution along the thickness can be defined by solving the onedimensional Fourier equation of heat conduction:

(12) 
The nonlinear temperature function obtained from Eq. (12) as following [13]:

(13) 
Since the DR method combined with the finite difference (FD) discretization method is employed to solve the equations, the integrations of Eq. (13) are computed numerically by discretizing the shell along the thickness direction.
In this study based on the FSDT, thickness effects on the buckling load are considered. According to the FSDT, the results are reliable for thin to moderately thick shells and the buckling load can be obtained for a variety of thicknesstolength ratios of shells. The displacement field based on the FSDT in the cylindrical coordinate system (x, θ, z) is as:

(14) 
where , and are the displacements corresponding to the coordinate system and are functions of the spatial coordinates; , , and are the middle surface displacements and , describe the rotations about the and axes, respectively (see Fig. 1). As stated, for obtaining the buckling load by the DR method, the equilibrium equations should be derived in the incremental form. Thus, all of the following governing equations are derived in the incremental form of variables. Based on the incremental nonlinear von Kármán strain–displacement relations, the strain components compatible with the displacement field of Eq. (14) are as follow:

(15) 
Using the Hooke’s law, the incremental constitutive thermoelastic relations can be defined by:

(16) 
The stress and moment resultants ( ) can be achieved by utilizing relevant integration through the thickness:



(17) 


By substituting Eqs. (15) and (16) into Eqs. (17), the incremental form of the constitutive relations in terms of displacement field are as follow:

(18) 
(19) 
Where and are the extensional, coupling, bending, and shear stiffness, respectively, which are obtained by:
, (i,j=1, 2,6)

(20) 
In which is the shear correction factor introduced by Reddy [45] and is equal to 5/6. According to the principle of minimum potential energy, the following force equilibrium equations in incremental form can be computed:

(21) 
It is evident that for mechanical buckling analysis the thermal terms must be raised. Substituting resultant forces and moments derived in Eqs. (17), (18) into Eq. (20) leads to a set of nonlinear displacement equilibrium equations in incremental form. As an example, the first equation of (20) is described in detail:
(22) 
In this paper, it is considered that mechanical buckling of the FG circular cylindrical shell is subjected to uniformly distributed axial compressive load . Whereas for thermal buckling analysis the FG cylinder is only under a thermal gradient along the thickness direction. Also, for both of mechanical and thermal buckling, the clamped and simply supported boundary conditions are applied around the circumference edges. These boundary conditions are set out below in terms of constraints on displacements, stress resultants and stress couples at :
(a) For Mechanical Buckling Analysis
Clamped—inplane movable:

(23) 
Simply supported—inplane movable:

(24) 
(b) For Thermal Buckling Analysis
Clamped—inplane movable:

(25) 
Simply supported—inplane movable:

(26) 
Since, solving the set of nonlinear equilibrium equations are very complicated and are not amenable to a closed form solution, in this study, the DR method with a finite difference discretization scheme was used. The DR method is a strong and reliable method for analyzing the nonlinear bending and buckling problems [1523, 45]. The DR method is an explicit iterative technique which is used to transfer a boundary value problem into timestepping initial value problem. This aim is obtained by adding artificial inertia and damping forces to the right side of Eqs. (21) as:

(27) 
In Eq. (27) LHS = lefthand side and , ( ) are elements of diagonal fictitious mass and damping matrices and , respectively. Here, to assure the numerical stability, the element of diagonal matrix M is obtained by the Gershgörin theorem as [47, 48]:

(28) 
where superscript n indicates the nth iteration and is the increment of fictitious time with its value assumed to be 1. Also, is the element of stiffness matrix which is achieved by:

(29) 
where is the lefthandside of the equilibrium equations (21). Furthermore, by applying the Rayleigh principle to each node, the instant critical damping factor for node i at the n^{th} iteration is as follows [49]:

(30) 
Hence, to make the elements of diagonal fictitious damping matrix , various values for diverse nodes are obtained at each direction as [49]:

(31) 
Finally, the velocity and acceleration terms should be substituted with the following equivalent central finitedifference expressions:

(32) 

(33) 
By replacing Eqs. (32) and (33) into the righthand side of Eq. (27), the equilibrium equations can be rearranged into an initial value format as:

(34) 
By integrating the velocities at the end of each load step, the incremental displacements can be computed:
(35) 
In order to compute the critical buckling load from the loaddisplacement curve, the total displacements of each load must be obtained. For this goal, the computed incremental displacements in each load step should be added to the displacements determined from the previous load steps as follows:
(36) 
It is evident that critical buckling load is a specified load in which a large amount of displacement occurs compared to the previous load steps.
Because of the terms of governing equations in the displacement field are quite long, for example, Eqs. at (34) have been written based on the force equilibrium equation (Eqs. (21)). In this paper, however, the developed numerical code is based on displacement equations. Therefore, the displacement equilibrium equations and Eqs. (34) to (36) with the related boundary conditions in their finite difference forms, constitute the set of equations for the sequential DR approach. For simplicity purposes, the DR algorithm, which is explained in [15, 16], is omitted. In order to clarify the present method for finding the buckling load of FG cylindrical shell with clamped boundary condition, the loaddisplacement curve is plotted and is shown at Fig 2. As can be seen, at a point in the graph when disproportionate increase in displacement occurs, the buckling behavior takes place.
In this analysis, the metal and ceramic phases of the FG shell are considered to be made of aluminum and alumina, respectively, in which the elasticity modules, thermal conductivity and thermal expansion coefficients, respectively, are , , for the former and these values are , , for the latter. Furthermore, the Poisson’s ratios of metal and ceramic are and , respectively. Also, the shell thickness is considered as .
At first, the results are carried out for the mechanical buckling behavior and then the thermal buckling analyses are presented. For both mechanical and thermal buckling analysis, the effects of grading indices, radiustothickness ratios, lengthtoradius ratios and boundary conditions are studied on the buckling load in detail. Obviously, the FSDT cannot present correct responses for thick shell (i.e. R/h=5) but the main goal of considering larger radiustothickness ratios is studying this ratio on the buckling behavior qualitatively which as it was observed in other papers, the trends of results of FSDT are similar to HSDT ones for higher radiustothickness ratios.
6.1 Mechanical Buckling Analysis
For verification of the present approach and relations of the critical buckling load, the obtained results are compared with references related with the mechanical buckling analysis of a circular cylindrical isotropic shell, and can be observed at Table 1.
As seen, there is a very good agreement between the current solution and those obtained by Refs. [14] and [34] for mechanical buckling of isotropic shell. In order to verify the current results for buckling analysis of FG shell, some comparison studies have been conducted between the present solution and the ones obtained by the Abaqus software [50] for different boundary conditions in Tables 2 and 3. As shown, the current results are verified for buckling analysis of FG shell. In order to model the shell structure, element S4R has been taken into the Abaqus software [50]. Also, according to the axisymmetric assumption of the cylindrical shell, 10 elements are considered in the thickness direction and 50 elements are assumed in the longitudinal direction. Definitely, the mesh size is selected using mesh sensitivity diagram to achieve the optimal number of elements that result in the least error and time consumption.
In Tables 46 the values of critical buckling load for different ratios of lengthtoradius, radiustothickness and grading indices are presented for different boundary conditions based on the simple powerlaw and MoriTanaka distributions. As it is expected, with increase of grading index and the tendency of material properties towards metal, the critical buckling load decreases based on both distribution models. In addition, for both of the distribution models the difference of critical buckling load among various grading indices is greater in clamped boundary condition rather than the simply supported one. This fact is seen in radiustothickness ratios of 5 to 30 and in many different ratios of lengthtoradius.
Fig. 2: Loaddisplacement curve for FG cylindrical Shell with clamped boundary condition.
While, for radiustothickness ratios above 30, not much of a difference is observed between the two boundary conditions. In this way, the greater lengthtoradius ratios result in decreasing the effect of grading index on the critical buckling load variations for both of the boundary conditions.
According to Tables 46, it can be concluded that at lengthtoradius ratios of 0.5, 1 and 5 with material grading indices of 1 and 2, the critical buckling load in the MoriTanaka model is smaller than that of a simplepower law. Furthermore, considering the lengthtoradius ratios of 0.5 and 1 and material grading indices of 5 and 10, the critical buckling load in different boundary conditions is greater for the case of the MoriTanaka model in comparison with that of a simple powerlaw. Though, for the case in which the lengthtoradius ratio is 5 and with radiustothickness ratios ranging from 5 to 20 along with the same material grading indices as the previous case, a lower critical buckling load in both boundary conditions exists in the MoriTanaka model rather than the simple powerlaw. This is, however, when the radiustothickness ratio is raised from 30 to 300 then the critical buckling load in the MoriTanaka model is higher than the one in simple powerlaw.
Table 1. Comparison of the critical buckling loads for simply supported (S) isotropic cylindrical shells under axial loads (MN) between the present solution and Refs. [42] and [14].


0.5 

1 

5 

Material 
Present study 
Ref. [34] 
Ref. [14] 

Present study 
Ref. [33] 
Ref. [14] 

Present study 
Ref. [34] 
Ref. [14] 

Aluminium 
5 
0.300 
0.292 
0.294 

0.280 
0.250 
0.271 

0.240 
0.239 
0.247 
10 
0.258 
0.258 
0.258 

0.258 
0.256 
0.258 

0.253 
0.258 
0.256 

20 
0.300 
0.286 
0.293 

0.271 
0.270 
0.270 

0.261 
0.251 
0.261 

30 
0.280 
0.295 
0.289 

0.262 
0.265 
0.264 

0.260 
0.260 
0.263 

100 
0.265 
0.270 
0.266 

0.265 
0.270 
0.266 

0.265 
0.266 
0.265 

300 
0.264 
0.267 
0.266 

0.264 
0.267 
0.266 

0.264 
0.266 
0.266 

Alumina 
5 
1.590 
1.583 
1.598 

1.473 
1.460 
1.472 

1.333 
1.295 
1.341 
10 
1.391 
1.401 
1.403 

1.380 
1.394 
1.403 

1.390 
1.390 
1.392 

20 
1.600 
1.546 
1.594 

1.472 
1.470 
1.468 

1.411 
1.391 
1.417 

30 
1.571 
1.611 
1.566 

1.440 
1.455 
1.435 

1.426 
1.417 
1.426 

100 
1.428 
1.470 
1.443 

1.428 
1.470 
1.443 

1.428 
1.445 
1.439 

300 
1.430 
1.458 
1.443 

1.430 
1.455 
1.443 

1.430 
1.434 
1.443 
Table 2. Comparison of the critical buckling load (MN) of the FG cylindrical shells based on the power law model with simply supported boundary conditions and lengthtoradius ratio of L / R = 0.5 with those of results obtained from the Abaqus software [50].
Volume fraction power 
R/h 

10 
5 
2 
1 

Abaqus [50] 
Present study 
Abaqus [50] 
Present study 
Abaqus [50] 
Present study 
Abaqus [50] 
Present study 

0.452 
0.446 
0534 
0.526 
0.669 
0.664 
0.854 
0.866 
5 
0.395 
0.392 
0.476 
0.456 
0.589 
0.600 
0.753 
0.760 
10 
0.442 
0.434 
0.513 
0.521 
0.701 
0.700 
0.868 
0.909 
20 
0.444 
0.450 
0.526 
0.532 
0.637 
0.667 
0.827 
0.827 
30 
0.408 
0.400 
0.470 
0.473 
0.626 
0.601 
0.785 
0.787 
100 
0.405 
0.400 
0.466 
0.470 
0.630 
0.609 
0.786 
0.788 
300 
Table 3. Comparison of the critical buckling load (MN) of the FG cylindrical shells based on the power law model with clamped boundary conditions and lengthtoradius ratio of L / R = 0.5 with those of results obtained from the Abaqus software [50].
Volume fraction power 
R/h 

10 
5 
2 
1 

Abaqus [50] 
Present study 
Abaqus [50] 
Present study 
Abaqus [50] 
Present study 
Abaqus [50] 
Present study 

0.694 
0.622 
0.860 
0.857 
1.165 
1.173 
1.687 
1.570 
5 
0.532 
0.533 
0.857 
0.854 
1.156 
1.159 
1.533 
1.532 
10 
0.673 
0.673 
0.784 
0.783 
0.970 
0.964 
1.319 
1.317 
20 
0.586 
0.580 
0.665 
0.673 
0.835 
0.837 
0.980 
0.977 
30 
0.411 
0.415 
0.472 
0.476 
0.612 
0.615 
0.789 
0.790 
100 
0.408 
0.410 
0.469 
0.472 
0.616 
0.618 
0.792 
0.793 
300 
According to Tables 46, it can be concluded that at lengthtoradius ratios of 0.5, 1 and 5 with material grading indices of 1 and 2, the critical buckling load in the MoriTanaka model is smaller than that of a simplepower law. Furthermore, considering the lengthtoradius ratios of 0.5 and 1 and material grading indices of 5 and 10, the critical buckling load in different boundary conditions is greater for the case of the MoriTanaka model in comparison with that of a simple powerlaw. Though, for the case in which the lengthtoradius ratio is 5 and with radiustothickness ratios ranging from 5 to 20 along with the same material grading indices as the previous case, a lower critical buckling load in both boundary conditions exists in the MoriTanaka model rather than the simple powerlaw. This is, however, when the radiustothickness ratio is raised from 30 to 300 then the critical buckling load in the MoriTanaka model is higher than the one in simple powerlaw. As indicated in Tables 46, it is also obvious that based on the powerlaw model, in simply supported boundary conditions and lengthtoradius ratios of 0.5 and 1, for various ranges of grading indices with an increase in radiustothickness ratio from 5 to 10 the critical buckling load decreases whereas with a continuation of increasing in the radiustothickness ratio from 10 to 20 the critical buckling load grows opposite to the preceding case.
Table 4. Critical buckling loads (MN) of the FG cylindrical shells for various material distributions, indices, boundary conditions, .
Boundary conditions 
k=1 

k=2 

k=5 

k=10 

P.L 
M.T 

P.L 
M.T 

P.L 
M.T 

P.L 
M.T 

Simply supported

5 
0.866 
0.727 

0.664 
0.627 

0.526 
0.523 

0.446 
0.500 
10 
0.760 
0.705 

0.600 
0.588 

0.456 
0.517 

0.392 
0.512 

20 
0.909 
0.773 

0.700 
0.669 

0.521 
0.571 

0.434 
0.548 

30 
0.827 
0.782 

0.667 
0.714 

0.532 
0.620 

0.450 
0.583 

100 
0.787 
0.697 

0.601 
0.610 

0.473 
0.550 

0.400 
0.540 

300 
0.788 
0.694 

0.609 
0.610 

0.470 
0.550 

0.400 
0.535 

Clamp supported 
5 
1.570 
1.334 

1.173 
1.092 

0.857 
0.815 

0.622 
0.734 
10 
1.532 
1.545 

1.159 
1.142 

0.854 
0.967 

0.533 
0.944 

20 
1.317 
1.043 

0.964 
0.923 

0.783 
0.821 

0.673 
0.800 

30 
0.977 
0.947 

0.837 
0.838 

0.673 
0.714 

0.580 
0.710 

100 
0.790 
0.730 

0.615 
0.640 

0.476 
0.587 

0.415 
0.570 

300 
0.793 
0.710 

0.618 
0.630 

0.472 
0.578 

0.410 
0.571 
Table 5. Critical buckling loads (MN) of the FG cylindrical shells for different material distributions, indices, boundary conditions, .
Boundary conditions 
k=1 

k=2 

k=5 

k=10 


P.L 
M.T 

P.L 
M.T 

P.L 
M.T 

P.L 
M.T 

Simply supported

5 
0.825 
0.649 

0.640 
0.587 

0.500 
0.485 

0.400 
0.476 

10 
0.730 
0.696 

0.588 
0.583 

0.457 
0.519 

0.386 
0.500 

20 
0.780 
0.696 

0.610 
0.595 

0.479 
0.532 

0.421 
0.522 

30 
0.769 
0.699 

0.615 
0.593 

0.462 
0.533 

0.402 
0.524 

100 
0.787 
0.697 

0.601 
0.600 

0.473 
0.540 

0.400 
0.526 

300 
0.788 
0.694 

0.609 
0.610 

0.470 
0.540 

0.400 
0.528 

Clamp supported 
5 
0.950 
0.842 

0.865 
0.714 

0.570 
0.592 

0.493 
0.551 

10 
0.827 
0.769 

0.700 
0.651 

0.569 
0.586 

0.474 
0.545 

20 
0.783 
0.692 

0.625 
0.600 

0.508 
0.538 

0.432 
0.524 

30 
0.790 
0.706 

0.628 
0.600 

0.500 
0.539 

0.429 
0.522 

100 
0.793 
0.708 

0.612 
0.609 

0.476 
0.540 

0.415 
0.528 

300 
0.793 
0.689 

0.616 
0.605 

0.472 
0.545 

0.410 
0.530 

Table 6. Critical buckling loads (MN) of the FG cylindrical shells for different material distributions, indices, boundary conditions and
Boundary conditions 
k=1 

k=2 

k=5 

k=10 


P.L 
M.T 

P.L 
M.T 

P.L 
M.T 

P.L 
M.T 

Simply supported

5 
0.733 
0.476 

0.571 
0.416 

0.431 
0.370 

0.376 
0.357 

10 
0.781 
0.434 

0.590 
0.384 

0.454 
0.332 

0.386 
0.326 

20 
0.777 
0.545 

0.611 
0.461 

0.455 
0.414 

0.390 
0.400 

30 
0.774 
0.615 

0.596 
0.552 

0.462 
0.490 

0.385 
0.470 

100 
0.787 
0.666 

0.601 
0.570 

0.473 
0.507 

0.400 
0.490 

300 
0.788 
0.685 

0.609 
0.588 

0.470 
0.533 

0.400 
0.520 

Clamp supported 
5 
0.750 
0.489 

0.612 
0.424 

0.469 
0.380 

0.382 
0.367 

10 
0.749 
0.435 

0.594 
0.385 

0.467 
0.340 

0.381 
0.335 

20 
0.842 
0.556 

0.639 
0.474 

0.480 
0.416 

0.400 
0.400 

30 
0.783 
0.625 

0.612 
0.555 

0.472 
0.500 

0.415 
0.470 

100 
0.789 
0.666 

0.615 
0.571 

0.476 
0.505 

0.415 
0.490 

300 
0.793 
0.685 

0.618 
0.588 

0.472 
0.529 

0.410 
0.528 

However, for grading indices of k=1, 2 with increase of radiustothickness ratio from 20 to 30 the critical buckling load decreases again. While for k=5, 10, by increasing the ratio from 20 to 30 the buckling load increases. Now, if we leave the lengthtoradius ratio constant at 5 and increase the radiustothickness ratio from 5 to 20, the critical buckling load will encounter a growth. It is obvious that for the higher ratios of radiustothickness there are not any significant differences between the results for various ranges of lengthtoradius ratios and grading indices.
In another case, considering the clamped boundary condition and lengthtoradius ratios of 0.5 and 1, increase of radiustothickness ratio from 5 to 100 will yield a reduction in critical buckling load. Moreover, at lengthtoradius ratio of 5, the increase in radiustothickness ratio from 5 to 10 will create a reduction of critical buckling load, however, the greater increase of the same ratio from 10 to 20 will cause a consequent growth in critical buckling load. Clearly, in higher values of radiustothickness ratios the obtained results do not change considerably for different ratios of lengthtoradius.
As seen in Tables 46, the effects of major parameters on the critical buckling load based on the MoriTanaka distribution, indicate the different behaviors compared to the ones given by the powerlaw model for some cases. As observed in Table 4, according to the MoriTanaka model, in the simply supported case at lengthtoradius ratio of 0.5, the critical buckling load decreases when an increase in the radiustothickness ratio is applied from 5 to 10 while it increases when the mentioned ratio ranges from 10 to 30. However, with increase of this ratio from 30 to 100 the obtained results are decreased again. Besides, considering the implementation of the clamed boundary condition, an increase of critical buckling load is seen by manipulating the radiustothickness ratio from 5 to 10, in contrast to the range of 10 to 30 which results in a decrease of the critical buckling load. A higher increase of the aforementioned ratio from 100 to 300 makes not much of a variation in critical buckling load for both of distribution models and boundary conditions. Also, the presented data from Table 5 based on the MoriTanaka model indicate that at a lengthtoradius ratio of 1, the critical buckling load grows as a consequence of enlarging the radiustothickness ratio from 5 to 20 under a simply supported boundary condition while it goes down in the case of clamped boundary condition. Regarding the previous case, by adding up the radiustothickness ratio from 20 to 300, not an observable change will appear in critical buckling load based on both rules of mixture. In Table 6, at lengthtoradius ratio of 5, increasing the radiustothickness ratio from 5 to 10 causes the critical buckling load to decrease under any boundary conditions, however, adding the mentioned ratio up to 30 makes it go up as well. It needs to be pointed out that similar to the previous cases, based on both distribution models, increase of the radiustothickness ratio from 100 to 300 does not create a great impact on critical buckling load no matter what type of boundary condition is implemented.
Based on the presented results of the Tables 46, it can be concluded that in different boundary conditions and concerning various grading indices, raising lengthtoradius ratio decreases the critical buckling load, as long as radiustothickness ratio ranges from 5 to 30 while for the ratios beyond 30, an increase in the lengthtoradius ratio has a small effect on critical buckling load. Moreover, considering the trend of variations within the critical buckling load as presented in Tables 46, it is inferred that as a result of increasing lengthtoradius ratio values of critical buckling load under both simply supported and clamped boundary conditions they are likely to converge into a common value so whereas at lengthtoradius ratio of 5 or in other words for long shells the values of critical buckling load does not differ greatly in both clamped and simply supported boundary conditions. It must be noted that the mentioned behaviors are more considerable for higher values of grading indices of the FG shell rather than the lower ones. It is concluded from the tables that nondimensional buckling loads predicted by a powerlaw model are larger than that of the Mori–Tanaka homogenization scheme except for some cases of aspect ratios. This can be concluded from the stiffness of the FG shell determined based on the different material distribution models and also the temperature distribution which is affected by the geometry of the shell.
6.2 Thermal Buckling Analysis
The thermal loading is applied in a manner that the temperature of metal surface is assumed to be constant at T_{m} =20 C while the temperature of ceramic surface increases incrementally, so that the nonlinear temperature variation is assumed along the thickness direction. The main objective, here, is to determine the critical temperature difference ∆T_{cr}=(T_{c}T_{m}) which causes the buckling. In order to verify the current solution for the thermal buckling behavior, a comparison study has also been carried out which is shown in Tables 7 and 8 between the DR results and the ones reported by Shahsiah and Eslami [8] for simply supported FG cylindrical shell under linear distribution of thermal gradient. According to [8], in this case the Poisson’s ratio is considered to be constant at ϑ=0.3 and the shell thickness is h=0.01m. Relying on the precision of the results from Tables 7 and 8 and the accuracy within solutions for FG cylindrical shells, in the following there will be a presentation of the results in correspondence with FG cylindrical shells imposed to nonlinear thermal distribution, unless stated otherwise, with different boundary conditions along with various geometrical parameters and grading indices.
In Fig. 3, a diagram has been drawn which represents the variation of the material grading index on the critical temperature difference at radiustothickness ratio of 10 and lengthtoradius ratio of 0.5 with implementation of different boundary conditions. It is observed that as the material grading index goes up, the critical temperature difference decreases. Additionally, the most significant variation occurs when the material grading index is zero. As seen, the critical temperature differences of FG cylindrical shells (R/ℎ=10, L/R=0.5) with a material grading index of k=0.5 and k=1.0 are different for simply supported and clamped boundary conditions which can be originated from rigidity of boundary conditions. In fact, increasing the rigidity of boundary conditions may cause the more significant trend of reduction in clamped boundary condition with respect to the simply supported one.
Table 7. Comparison the critical temperature difference ( ) of the simply supported FG cylindrical shell (k=1, R=0.5m) based on linear temperature distribution obtained by the DR method to the results reported by Ref. [8].
L/R 
h/R 

0.5 
0.3 
0.15 

Present study 
Ref. [8] 
Present study 
Ref. [8] 
Present study 
Ref. [8] 

40 
40 
80 
100 
140 
140 
4.60E03 

60 
80 
120 
120 
210 
200 
6.40E03 

80 
80 
180 
180 
320 
320 
8.22E03 

100 
100 
240 
240 
420 
400 
1.00E02 

Table 8. Comparison the critical temperature difference ( ) of the simply supported FG cylindrical shell (k=1, =0.5) based on linear temperature distribution obtained by the DR method to the results reported by Ref. [8].
h(m) 
R(m) 

0.01 
0.007 
0.005 

Present study 
Ref. [8] 
Present study 
Ref. [8] 
Present study 
Ref. [8] 

1000 
1000 
300 
300 
80 
100 
0.625 

520 
500 
170 
180 
55 
60 
0.90 

330 
340 
120 
120 
80 
60 
1.18 

230 
260 
90 
100 
60 
60 
1.45 

180 
220 
70 
80 
60 
60 
1.73 

160 
160 
60 
60 
60 
60 
2.00 

Fig. 3. Critical temperature difference versus k for FG cylindrical shell (R⁄h=10, L⁄R=0.5) with different boundary conditions based on linear and nonlinear thermal distributions.
Figure 4 depicts the critical temperature difference for different radiustothickness ratios with regards to a material grading index of K=0.5 and lengthtoradius ratio of L⁄R=0.5 for both clamped and simply supported boundary conditions. As noticed, due to an increase of radiustothickness ratio the critical temperature difference decreases for both boundary conditions.
In order to consider the effect of boundary conditions on the critical temperature difference, Fig. 5 is presented for different radiustothickness ratios of simply supported and clamped FG cylindrical shells with a material grading index of k=0.5 and lengthtoradius ratio of L⁄R=1 based on nonlinear thermal distributions. It is discerned that the critical temperature difference in a clamped boundary condition is more noticeable than that of a simply supported one. Also, at high radiustothickness ratios, no major difference is identified between clamped and simply supported boundary conditions.
The difference between the distribution models of the simple powerlaw and the MoriTanaka equation for predicting the critical temperature difference of simply supported and clamped FG cylindrical shell is illustrated in Fig. 6 for different ratios of the material grading index with respect to radiustothickness ratio of R⁄h=5 and lengthtoradius ratio of L⁄R=1. It is identifiable that the critical temperature variation is greater when the FG material is modeled by the Mori Tanaka theory compared to the case of a simple powerlaw.
In Table 9 the critical temperature difference is shown for different values of grading indices and ratios of lengthtoradius and radiustothickness. As indicated, with increase of the lengthtoradius ratio the critical temperature difference decreases. However, in higher values of radiustothickness ratios (R⁄h=100) the increase of lengthtoradius ratio does not have any significant effects on the results.
Fig. 4. Critical temperature difference versus for FG cylindrical shell ( , ) with different boundary conditions based on linear and nonlinear thermal distributions.
Fig. 5. Critical temperature difference versus for the FG cylindrical shell ( , ) with different boundary conditions based on nonlinear thermal distributions.
In the present paper, the mechanical and thermal buckling analysis of FG cylindrical shells have been considered for clamped and simply supported boundary conditions based on FSDT and large deflection von Kármán equations. The mechanical and thermal properties of the constituent components of the FG shells have been assumed to vary continuously along the thickness direction according to the simple powerlaw and the Mori–Tanaka distributions. To predict precisely the elastic properties of FG shell, the variable Poisson’s ratio is considered along the thickness direction. Also, for thermal buckling analysis the nonlinear temperature distribution is assumed. The critical buckling load is predicted based on mechanical/thermal load–displacement curve obtained by solving the incremental form of nonlinear equilibrium equations. The DR method in conjunction with the central finite difference discretization technique is used to solve the incremental formulations. Finally, after verification of the present solutions, a detailed parametric study is carried out to investigate the effects of boundary conditions, rules of mixture, grading indices, radius to thickness and lengthtoradius ratios on the mechanical and thermal buckling loads.
Fig. 6. Critical temperature difference of the (a) simply supported and (b) clamped FG cylindrical shell ( , ) vs for different material model based on nonlinear thermal distribution
Table 9. Critical temperature difference ( ) of the FG cylindrical shells for different grading indices, boundary conditions and geometries based on simple powerlaw method and nonlinear temperature distribution.

k=0 

k=1 

k=5 

Boundary conditions 



Clamp supported

5 
600 
600 

480 
400 

360 
300 

10 
580 
480 

380 
300 

300 
240 

20 
580 
400 

220 
140 

180 
120 

30 
440 
380 

160 
120 

140 
80 

100 
180 
180 

80 
80 

60 
60 

Simply supported 
5 
520 
440 

440 
320 

340 
260 

10 
400 
360 

300 
220 

240 
180 

20 
400 
260 

160 
140 

140 
120 

30 
340 
280 

140 
120 

100 
80 

100 
140 
140 

60 
60 

60 
60 

References
[1] Sofiyev AH. About an approach to the determination of the critical time of viscoelastic functionally graded cylindrical shells. Compos B Eng 2019; 156:156–165.
[2] Reddy JN, Chin, CD. Thermomechanical analysis of functionally graded cylinders and plates. J Therm Stresses 1998; 26:593–626.
[3] Li SR, Batra RC. Buckling of axially compressed thin cylindrical shells with functionally graded middle layer. ThinWalled Struct 2006; 44:1039–1047.
[4] Huang H, Han Q, Wei D. Buckling of FGM cylindrical shells subjected to pure bending load. Comp Struct 2011; 93:2945–2952.
[5] Shariyat M. Dynamic thermal buckling of suddenly heated temperaturedependent FGM cylindrical shells, under combined axial compression and external pressure. Int J Solids Struct 2011; 93:2945–2952.
[6] Huang H, Han Q. Nonlinear elastic buckling and postbuckling of axially compressed functionally graded cylindrical shells. Int J Mec Sci 2009; 51:500–507.
[7] Huang H, Han Q. Buckling of imperfect functionally graded cylindrical shells under axial compression. Euro J Mech A/Solids 2008; 27:1026–1036.
[8] Shahsiah R, Eslami MR. Thermal buckling of functionally graded cylindrical shell. J Therm Stresses 2003; 26:277–294.
[9] Ebrahimi F, Sepiani HA. Vibration and buckling analysis of cylindrical shells made of functionally graded materials under combined static and periodic axial forces. Advanced Composites Letters 2010; 19:67–74.
[10] Shen HS. Postbuckling analysis of pressureloaded functionally graded cylindrical shells in thermal environments. Engng Struct 2003; 25:487–497.
[11] Shen HS, Noda N. Postbuckling of FGM cylindrical shell under combined axial and radial mechanical loads in thermal environments. Int J Solids Struct 2005; 42:4641–4662.
[12] Mirzavand B, Eslami MR. Thermoelastic stability of imperfect functionally graded cylindrical shells. J Mech Mater struct 2008; 3:1561–1572.
[13] Duc ND, Thang PT. Nonlinear response of imperfect eccentrically stiffened ceramicmetalceramic FGM thin circular cylindrical shells surrounded on elastic foundations and subjected to axial compression. Comp Struct 2014; 110:200–206.
[14] Khazaeinejad P, Najafizadeh MM. Mechanical buckling of cylindrical shells with varying material properties. J Mech Eng Sci 2010; 224:1551–1557.
[15] Golmakani ME, Kadkhodayan M. Nonlinear bending analysis of annular FGM plates using higherorder shear deformation plate theories. Compos Struct 2011; 93: 973–982.
[16] Golmakani ME, Kadkhodayan M. Large deflection analysis of circular and annular FGM plates under thermomechanical loadings with temperaturedependent properties. Compos Part B 2011; 42: 614–25.
[17] Golmakani ME. Large deflection thermoelastic analysis of shear deformable functionally graded variable thickness rotating disk. Compos Part B 2013; 45:1143–55.
[18] Golmakani ME. Nonlinear bending analysis of ringstiffened functionally graded circular plates under mechanical and thermal loadings. Int J Mec Sci 2014; 79: 130–142.
[19] Golmakani ME, Kadkhodayan M. An investigation into the thermoelastic analysis of circular and annular functionally graded material plates. Mech Adv Mater Struct 2014; 21: 1–13.
[20] Golmakani ME, Kadkhodayan M. Large deflection thermoelastic analysis of functionally graded stiffened annular sector plates. Int J Mech Sci 2013; 69: 94–106.
[21] Golmakani ME, Alamatian J. Large deflection analysis of shear deformable radially functionally graded sector plates on twoparameter elastic foundations. Euro J Mech A/Solids 2013; 42: 251–265.
[22] Alamatian J, Golmakani ME. Large deflection analysis of the moderately thick general theta ply laminated plates on nonlinear elastic foundation with various boundary conditions. Mech Res Commun 2013; 51: 78–85.
[23] Golmakani ME, Mehrabian M. Nonlinear bending analysis of ringstiffened circular and annular general angleply laminated plates with various boundary conditions. Mech Res Commun 2014; 59: 42–50.
[24] Wang Y, Feng C, Zhao Z, Yang J. Eigenvalue Buckling of Functionally Graded Cylindrical Shells Reinforced with Graphene Platelets (GPL). Compos Struct 2017; 202: 3846
[25] Wang Y, Feng C, Zhao Z, Lu F, Yang J. Torsional buckling of graphene platelets (GPLs) reinforced functionally graded cylindrical shell with cutout. Compo Struct 2018; 197: 72–79.
[26] Yiwen Ni, Zhenzhen Tong, Dalun Rong, Zhenhuan Zhou, Xinsheng Xu. Accurate thermal buckling analysis of functionally graded orthotropic cylindrical shells under the symplectic framework. ThinWalled Struct 2018; 129: 1–9.
[27] Trabelsi S, Frikha A, Zghal S, Dammak F. A modified FSDTbased four nodes finite shell element for thermal buckling analysis of functionally graded plates and cylindrical shells. Eng. Struct 2019; 178: 444–459.
[28] Nam VH, Phuong NT, Van Minh K, Hieu PT. Nonlinear thermomechanical buckling and postbuckling of multilayer FGM cylindrical shell reinforced by spiral stiffeners surrounded by elastic foundation subjected to torsional loads. EUR J MECH ASOLID 2018;72: 393406.
[29] Thang P.T, Dinh Duc N., NguyenThoi T. Thermomechanical buckling and postbuckling of cylindrical shell with functionally graded coatings and reinforced by stringers. Aerosp Sci Technol 2017;66: 392401.
[30] Golmakani ME, Sadraee Far MN, Moravej M. Dynamic relaxation method for nonlinear buckling analysis of moderately thick FG cylindrical panels with various boundary conditions. JMST 2016; 30: 5565–5575.
[31] RezaieePajand M, Pourhekmat D, Arabi E. Thermomechanical stability analysis of functionally graded shells. Eng. Struct 2019; 178:1–11.
[32] Wang, Y., Feng, C., Zhao, Z., & Yang, J. Buckling of graphene platelet reinforced composite cylindrical shell with cutout. INT J STRUCT STAB DY 2018; 18: 1850040.
[33] Zghal S.,Frikha A., Dammak F.,Static analysis of functionally graded carbon nanotubereinforced plate and shell structures. Compos Struct 2017,176: 11071123.
[34] Trabelsi S., Frikha A., Zghal S., Dammak F., Thermal postbuckling analysis of functionally graded material structures using a modified FSDT. IJMS 2018;144: 7489.
[35] Alijani A., Darvizeh M., Darvizeh A., Ansari R., On nonlinear thermal buckling analysis of cylindrical shells. ThinWalled Struct 2015; 95: 170182
[36] Zghal S., Frikha A. Dammak F., Free vibration analysis of carbon nanotubereinforced functionally graded composite shell structures. APPL MATH MODEL 2018; 53:132155.
[37] Zghal S., Frikha A. Dammak F., Dammak, Mechanical buckling analysis of functionally graded powerbased and carbon nanotubesreinforced composite plates and curved panels. COMPOS PART BENG, 2018;150:165183
[38] Frikha A. Zghal S. Dammak F., Dynamic analysis of functionally graded carbon nanotubesreinforced plate and shell structures using a double directors finite shell element. AEROSP SCI TECHNOL 2018; 78: 438451.
[39] Frikha A. Zghal S. Dammak F., Finite rotation three and four nodes shell elements for functionally graded carbon nanotubesreinforced thin composite shells analysis. COMPUT METHOD APPL M 2018; 329: 289311.
[40] Zghal S., Frikha A. Dammak F., Nonlinear bending analysis of nanocomposites reinforced by graphene nanotubes with finite shell element and membrane enhancement. Eng. Struct 2018; 158: 95109.
[41] Reddy JN. Analysis of functionally graded plates, INT J NUMER METH ENG 2000; 47:663684.
[42] Klusemann B, Svendsen B. Homogenization methods for multiphase elastic composites: Comparisons and benchmarks. Technische Mechanik 2010; 30:374–86.
[43] Prakash T, Singha MK, Ganapathi M. Thermal postbuckling analysis of FGM skew plates. Eng Struct 2008; 30:22–32.
[44] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall 1973; 21:571–4.
[45] Reddy JN. Mechanics of Laminated Composite Plates and Shells, Second Edition, CRC Press, New York; 2004.
[46] Kadkhodayan M, Zhang LC, Sowerby R. Analysis of wrinkling and buckling of elastic plates by DXDR method. Comput Struct 1997; 65:561–74.
[47] Underwood P. Dynamic relaxation, in computational methods for transient analysis, Chapter 5. Amsterdam, Elsevier; 1983.
[48] Zhang LC, Yu TX. Modified adaptive dynamic relaxation method and application to elastic–plastic bending and wrinkling of circular plate. Comput Struct 1989; 33:609–14.
[49] Zhang LC, Kadkhodayan M, Mai YW. Development of the maDR method. Comput Struct 1994; 52:1–8.
[50] Abaqus. Ver 6.101, Dassualt Systems Inc.; 2010.