Authors
^{1} Department of Mechanical Engineering, Islamshahr Branch, Islamic Azad University, Tehran, Iran
^{2} Composite Materials and Technology Center, Tehran, Iran
Abstract
Keywords
A Semianalytical Solution for 3D Dynamic Analysis of Thick Continuously Graded Carbon Nanotubereinforced Annular Plates Resting on a Twoparameter Elastic Foundation
V. Tahouneh^{*a}, J. Eskandari Jam ^{b}
^{a} Department of Mechanical Engineering, Islamshahr Branch, Islamic Azad University, Tehran, Iran
^{b }Composite Materials and Technology Center, Tehran, Iran
paper INFO 

ABSTRACT 
Paper history: Received 30 July 2014 Received in revised form 16 November 2014 Accepted 21 November 2014 
The The main objective of this research paper is to present 3D elasticity solution for free vibration analysis of elastically supported continuously graded carbon nanotubereinforced (CGCNTR) annular plates. The volume fractions of oriented, straight singlewalled carbon nanotubes (SWCNTs) are assumed to be graded in the thickness direction. An equivalent continuum model based on the EshelbyMoriTanaka approach is employed to estimate the effective constitutive law of the elastic isotropic medium (matrix) with oriented, straight carbon nanotubes (CNTs). A semianalytical approach composed of 2D differential quadrature method and series solution is adopted to solve the equations of motion. The novelty of the present work is to exploit EshelbyMoriTanaka approach in order to reveal the impacts of the volume fractions of oriented CNTs and different CNTs distributions on the vibrational characteristics of CGCNTR annular plates. 



Keywords: Threedimensional free vibration Continuously graded carbon nanotubereinforced Annular plates Twoparameter elastic foundations 


© 2014 Published by Semnan University Press. All rights reserved. 
Recently, Recently, Nanocomposites have significant importance for engineering applications that require high levels of structural performance and multifunctionality. Carbon Nanotubes (CNTs) have demonstrated exceptional mechanical, thermal and electrical properties. These materials are considered as one of the most promising reinforcement materials for highperformance structural and multifunctional composites with vast application potentials [1]. A detailed summary of the mechanical properties of CNTs can be found in [2]. The exceptional mechanical properties of CNTs have shown great promise for a wide variety of applications, such as nanotransistors, nanofillers, semiconductors, hydrogen storage devices, structural materials, molecular sensors, fieldemissionbased displays, and fuel cells [3]. The addition of nanosized fibers or nanofillers, such as CNTs, can further increase the merits of polymer composites [4]. These nanocomposites, easily processed due to the small diameter of the CNTs, exhibit unique properties [5,6], such as enhanced modulus and tensile strength, high thermal stability and good environmental resistance. This behavior, combined with their low density makes CNTs suitable for a broad range of technological sectors such as telcommunications, electronics [7] and transport dustries, especially for aeronautic and aerospace applications where the reduction in weight is crucial in order to reduce the fuel consumption. For example, Qian et al. [8] showed that the addition of 1 wt.% (i.e. 1% by weight) multiwall CNT to polystyrene resulted in 3642% and ~25% increase in the elastic modulus and the break stress of the nanocomposite properties, respectively. In addition, Yokozeki et al. [9] reported the retardation of the onset of matrix cracking in the composite laminates containing the cupstacked CNTs compared to those without the cupstacked CNTs. Most studies on CNTReinforced Composites (CNTRCs) have focused on their material properties [10–14]. Shen [15] for the first time suggests that the nonlinear bending behavior can be considerably improved through the use of a functionally graded distribution of CNTs in the matrix. He introduced the CNT efficiency parameter to account load transfer between the nanotube and polymeric phases. Compressive postbuckling and thermal buckling behavior of functionally graded nanocomposite plates reinforced by aligned, straight SingleWalled CNTs (SWCNTs) subjected to inplane temperature variation was reported by Shen and Zhu [16] and Shen and Zhang [17]. They find that in some cases the CNTRC plate with intermediate CNT volume fraction does not have intermediate buckling temperature and initial thermal postbuckling strength. Moreover, Ke et al. [18] investigated the nonlinear free vibration of functionally graded CNTRC Timoshenko beams. They find that both linear and nonlinear frequencies of functionally graded CNTRC beam with symmetrical distribution of CNTs are higher than those of beams with uniform or unsymmetrical distribution of CNTs. To the best of authors’ knowledge the review of open literature showed that the studies on functionally graded CNTRCs were restricted to nanocomposite structures having graded aligned, straight CNTs in the thickness direction, and effective material properties of CNTRCs were estimated through the extended rule of mixture. This fact motivates us to employ the EshelbyMoriTanaka approach for calculating the elastic stiffness properties of nanocomposite materials reinforced by graded oriented, straight CNTs. Plates resting on elastic foundations have found considerable applications in structural engineering problems. Reinforcedconcrete pavements of highways, airport runways, foundation of storage tanks, swimming pools, and deep walls together with foundation slabs of buildings are wellknown direct applications of these kinds of plates. The underlying layers are modeled by a Winklertype elastic foundation. The most serious deficiency of the Winkler foundation model is having no interaction between the springs. In other words, the springs in this model are assumed to be independent and unconnected. The Winkler foundation model is fairly improved by adopting the Pasternak foundation model, a twoparameter model, in which the shear stiffness of the foundation is considered. A closedform solution for the vibration frequencies of simply supported Mindlin plates on Pasternak foundations and subjected to biaxial initial stresses was presented by Xiang et al. [19]. The buckling load of Mindlin plates on Pasternak foundations was obtained in terms of the thin plate solution. Based on firstorder shear deformation plate theory, the buckling and vibration analysis of moderately thick laminates on Pasternak foundations was presented by Xiang et al. [20]. The effects of foundation parameters, transverse shear deformation, and rotary inertia and the number of layers on the buckling and vibration of crossply laminates were examined. Wang et al. [21] presented relationships between the buckling loads of simply supported plates on a Pasternak foundation determined by classical Kirchhoff plate theory, Reissner–Mindlin plate theory, and Reddy plate theory. The vibration of polar orthotropic circular plates on an elastic foundation was investigated by Gupta et al. [22]. The Mindlin shear deformable plate theory was employed and the Chebyshev collocation method was applied to obtain the frequency parameters of the circular plates. Ju et al. [23] developed a finite element model to study the vibration of Mindlin plates with multiple stepped variations in thickness and resting on nonhomogeneous elastic foundations. Gupta et al. [24,25] studied the effect of elastic foundation on axisymmetric vibrations of polar orthotropic circular plates of variable thickness by taking approximating polynomials in Rayleigh–Ritz method. Laura and Gutierrez [26] analyzed the free vibration of a solid circular plate of linearly varying thickness attached to Winkler foundation using the Ritz method. Matsunaga [27] analyzes the natural frequencies and buckling stresses of FG plates using a higher order shear deformation theory which are based on the through the thickness series expansion of the displacement components. Zhou et al. [28] used Ritz method to analyze the freevibration characteristics of rectangular thick plates resting on elastic foundations. Matsunaga [29] investigated a twodimensional, higherorder theory for analyzing the thick simply supported rectangular plates resting on elastic foundations. Tahouneh and Yas [30] investigated the free vibration analysis of thick FG annular sector plates on Pasternak elastic foundations using DQM. Tahouneh and Yas [31] studied free vibration characteristics of thick multidirectional functionally graded annular sector plates under various boundary condition using 2D differential quadrature method. They show that a graded ceramic volume fraction in two directions has a higher capability to reduce the natural frequency than conventional 1D FGM. Tahouneh et al. [32] investigated the effect of continuous grading fiber reinforced on the vibrational response of thick annular plates using DQ method. More recently, Tahouneh [33] studied the threedimensional free vibration analysis of Continuous Grading Fiber Reinforced (CGFR) sector plates with simply supported radial edges and arbitrary boundary conditions on their circular edges using 2D DQ method. Jam et al. [34] studied free vibrations of threeparameter functionally graded plates resting on Pasternak foundations using GDQ method. Nie and Zhong [35] investigated threedimensional vibration of FG circular plates using semianalytical method. Dong [36] developed a threedimensional free vibration analysis of FG annular plates using the ChebyshevRitz method. Cheng and Batra [37] used Reddy’s thirdorder plate theory to study steady state vibrations and buckling of a simply supported functionally gradient isotropic polygonal plate resting on a Pasternak elastic foundation and subjected to uniform inplane hydrostatic loads. Malekzadeh [38] studied free vibration analyses of functionally graded plates on elastic foundations based on the threedimensional elasticity. In structural mechanics, one of the most popular semianalytical methods is Differential Quadrature Method (DQM) [39], remarkable success of which is demonstrated by many researchers in vibration analysis of plates, shells, and beams. Liu and Liew [40], and Liew and Liu [41] presented DQM for free vibration analysis of Mindlin isotropic circular and annular sector plates with various types of boundary conditions. A new version of the DQM was extended by Wang and Wang [42] to analyze the free vibration of thin circular sector plates with six combinations of boundary conditions. Liew et al. [43] employed DQM for free vibration analysis of moderately thick plates on Winkler foundation. Gupta et al. [44] studied the free vibration analysis of nonhomogeneous circular plate of nonlinear thickness variation using the DQM. Nie and Zhong [45] studied the free vibration of FG plates without elastic foundation using DQM. They assume the material properties of the FG plate have an exponentlaw variation along the thickness, radial direction or both directions. The mathematical fundamental and recent developments of differential quadrature method as well as its major applications in engineering are discussed in detail in the book by Shu [46].This paper is motivated by the lack of studies in the technical literature concerning to the threedimensional vibration analysis of a Continuously Graded Carbon NanotubeReinforced (CGCNTR) annular plates resting on a twoparameter foundation. To the authors’ best knowledge, research on the vibration of thickness a continuously graded carbon nanotubereinforced (CGCNTR) annular plates resting on a twoparameter foundation based on the threedimensional theory of elasticity has not been seen yet. In this study, the volume fractions of oriented, straight SingleWalled Carbon Nanotubes (SWCNTs) are assumed to be graded in the thickness direction. An equivalent continuum model based on the EshelbyMoriTanaka approach is employed to estimate the effective constitutive law of the elastic isotropic medium (matrix) with dispersed elastic inhomogeneities (oriented CNTs). A sensitive analysis is performed, and the natural frequencies are calculated for different sets of boundary conditions and different combinations of the geometric, and foundational parameters. Therefore, very complex combinations of the material properties, boundary conditions, and foundation stiffness are considered in the present semianalytical solution approach.
2. Problem Description
The geometric configuration of a CGCNTR annular thick plate is shown in Fig. 1. a, b and h are outer/inner radius and thickness of the plate, respectively. The plate is supported by an elastic foundation with Winkler’s (normal) and Pasternak’s (shear) coefficients. The deformations defined with reference to a cylindrical coordinate system(r, θ, z) are u_{r} ,u_{θ} and u_{z} in the r , θ and z directions, respectively. We assume that the CGCNTR annular plate is made of a mixture of oriented, straight SWCNT, graded distribution in the thickness direction, and polymer matrix which is assumed to be isotropic [47].
2.1. Estimation of Effective Material Properties of CNTRC
In this study, we exploit an equivalent continuum model based on the EshelbyMoriTanaka approach in order to estimate the effective constitutive law of the elastic isotropic medium (matrix) with dispersed elastic inhomogeneities (carbon nanotubes). The major step towards modeling materials with fully dispersed inhomogeneities was taken by Mori and Tanaka [48]. In particular, the presence of multiple inclusions and boundary conditions and their interactions are accounted by them. Giordano et al. [49] used the homogenization procedure, based on the Eshelby theory, under small deformations and small volume fractions of the embedded phases, to determine the bulk and shear moduli and Landau coefficients of the composite material. Previous studies have established the validity of the EshelbyMoriTanaka approach in determining the effective properties of composites reinforced with misaligned, carbon fibres, and with carbon nanotubes [14,50,51]. In this paper, the proposed model is framed with the Eshelby theory for elastic inclusions. The original theory of Eshelby [52,53] is restricted to one single inclusion in a semiinfinite elastic, homogeneous and isotropic medium. The theory, generalized by MoriTanaka, allows extending the original approach to the practical case of multiple inhomogeneities embedded into a finite domain. The EshelbyMoriTanaka approach, based on the equivalent elastic inclusion idea of Eshelby and the concept of average stress in the matrix according to MoriTanaka, is also known as the equivalent inclusionaverage stress method [54].
2.1.1. Nanocomposite Reinforced by Aligned, Straight CNTs
a linear elastic polymer matrix reinforced by a large number of dispersed straight CNTs is considered. First, we consider a polymer composite reinforced with aligned and straight CNTs. According to Benveniste’srevisitation [50], the following equation of the effective elastic tensor is obtained:
(1) 
Figure 1. The sketch of an elastically supported thick continuously graded carbon nanotube reinforced annular plate and setup of the coordinate system. 
Where f_{r} and f_{m} are the fiber and matrix volume fractions, respectively, I is the identity tensor, C_{m} is the stiffness tensor of the matrix material, C_{r} is the stiffness tensor of the equivalent fiber, and A_{r} is the dilute mechanical strain concentration tensor for the fiber which is obtained through the following formula:
(2) 
The tensor S is Eshelby’s tensor, as given by Eshelby [52] and Mura [55]. The terms enclosed with angle brackets in Eq. (1) represent the average value of the term over all orientations defined by transformation from the local fiber coordinates to the global coordinates , as it is shown in Fig. 2. The matrix is assumed to be elastic and isotropic, with Young’s modulus E_{m} and Poisson’s ratio . Each straight CNT is modeled as a long fiber with transversely isotropic elastic properties. Therefore, the composite is also transversely isotropic. The substitution of nonvanishing components of the Eshelby tensor S for a straight, long fiber along the  direction in Eq. (2) gives the dilute mechanical strain concentration tensor. Then inserting A_{r} into Eq. (1) gives the tensor of effective elastic moduli of the composite reinforced by aligned and straight CNTs.
In particular, the Hill’s elastic moduli are found as [47]:
(3) 

(4) 

(5) 

(6) 

(7) 
Where k, l, m, n, and p are Hill’s elastic moduli of the composite; k is the planestrain bulk modulus normal to the fiber direction, n is the uniaxial tension modulus in the fiber direction, l is the associated cross modulus, m and p are the shear moduli in planes normal and parallel to the fiber direction, respectively. In addition, k_{r},l_{r} , m_{r} ,n_{r} , and p_{r} are the Hill’s elastic moduli for the reinforcing phase (CNTs). The elastic moduli parallel and normal to CNTs are related to Hill’s elastic moduli through what follows:
(8) 
2.1.2. Nanocomposite Reinforced by Oriented, Straight CNTs
In this section, the influence of oriented, straight CNTs is investigated. The orientation of a straight CNT is characterized using two Euler angles α and β, as shown in Figure 2. The orientation distribution of CNTs in the CNTRC is characterized using a probability density function for oriented nanotubes in which case the composite is isotropic. The base vectors and of the global and the local coordinate systems are related via the transformation matrix g as follows:
(9) 
Figure 2. Representative volume element (RVE) including straight CNTs. 
Where g is given by
(10) 
The orientation distribution of CNTs in a composite is characterized using a probability density function P(α, β) satisfying the normalization condition as the following [47]:
(11) 
If CNTs are completely oriented, the density function is the following:
(12) 
According to the MoriTanaka method, the strain and the stress of the CNT are related to the stress of matrix by the following equations:
(13) 
Where the strain concentration tensor is given by Eq. 2. Then the average strain and stress in all oriented CNTs are written as the following:
(14) 
Using The angle brackets represent the average over special orientations. Using the average theorems and in conjunction with the effective constitutive relation , one can get the effective modulus of the composite according to Eq. (1). When CNTs are completely oriented in the matrix, the composite is then isotropic, and its bulk modulus K and shear modulus G are derived as what follows:
(15) 

(16) 
Where
(17) 

(18) 

(19) 

(20) 
Where K_{m} and G_{m} are the bulk and shear moduli of the matrix, respectively. In addition, k_{r}, m_{r}, n_{r} and l_{r} are the Hill’s elastic moduli for the reinforcing phase. The effective Young’s modulus and Poisson’s ratio of the material are obtained using the following equations:
(21) 

(22) 
In addition, V_{f} and V_{m} are the volume fractions of the CNTs and the matrix, which satisfy the relationship of V_{f }+ V_{m }=1 imilarly, mass density is calculated using the following equation:
(23) 
Where and are the mass density of the CNTs and the matrix, respectively. In order to examine the effect of different CNTs distribution on the free vibration characteristics of CGCNTR annular plates resting on elastic foundations, various types of material profiles are considered through the plate thickness . In this work, we assume only linear distribution of CNTs volume fraction for the different types of the CGCNTR annular plate, as follow:
(24) 

(25) 

(26) 

(27) 
Where is the volume fraction of CNTs [11,15,56] that is calculated by the mass fraction of nanotubes, , assuming two phases and no trapped air, using the following equation [11]:
(28) 
Where is the ratio of nanotube to matrix density. It is worth noting that corresponds to the uniformly distributed CNTR annular plate referred to as CNTRUD. With defined in Eq. (28), both the CGCNTR plate and CNTRUD plate have the same value of CNTs mass fraction. For type V, the top surface of the plate is CNTrich, (Fig. 3).
Figure 3. Variations of CNTs volume fractions through the thickness of the plate for different types of CNT distribution. 
As can be seen in Fig. 3, for type , the distribution of CNTs reinforcements is inversed and the bottom surface of the plate is CNTrich, referred to as CGCNTR . For type X, a midplane symmetric graded distribution of CNTs reinforcements is achieved and both top and bottom surfaces are CNTrich, referred to as CGCNTRX. For type, the distribution of CNTs reinforcements is inversed and both top and bottom surfaces are CNTpoor, whereas the reference surface is CNTrich, referred to as CGCNTR.
3. Teoretical Formulations
The mechanical constitutive relations that relate the stresses to the strains are as follows [54]:
(29) 
Where and are the Lame constants, is the infinitesimal strain tensor and is the Kronecker delta. In the absence of body forces, the equations of motion are as follows:
(30) 
The infinitesimal strain tensor is related to the displacements as follows:
(31) 
Where u_{r}, u_{θ} and u_{z} are displacement components along the r, θ and z directions, respectively. Moreover, ε_{r}, ε_{θ}, ε_{z}, γ_{θz}, γ_{rθ} and γ_{rz} are strain components. Upon substitution Eq. (31) with (29) and then with (30), the following equations of motion are obtained in terms of displacement components:
(32)
(33)


(34) 
Equations (32 and 33) represent the inplane equations of motion along the r and θaxes, respectively; and Eq. (34) is the transverse or outofplane equation of motion. The related boundary conditions at z=h/2 and h/2 are as follows:
at z=h/2:
(35) 
at z=h/2:
(36) 
Where is the component of stress tensor; and are Winkler and shearing layer elastic coefficients of the foundation. Different types of classical boundary conditions at the circular edges of the plate can be stated as
Clamped(r=b)  Clamped(r=a): 
(37) 
at r=a
at r=b
Simply supported(r=b)  Clamped(r=a): 
(38) 
at r=b
at r=a
Simply supported(r=b)  Clamped(r=a): 
(39) 
at r=a
at r=b
4. DQM Solution for Equations of Motion and Boundary Conditions
It is necessary to develop appropriate methods for investigating the mechanical responses of CGCNTR structures. But, due to the complexity of the problem, it is difficult to obtain the exact solution. In this paper, the Differential Quadrature Method (DQM) approach is used to solve the governing equations of CGCNTR annular plates.
One can compare DQM solution procedure with the other two widely used traditional methods for plate analysis, i.e., RayleighRitz method and FEM. The main difference between the DQM and the other methods is how the governing equations are discretized. In DQM the governing equations and boundary conditions are directly discretized, and thus the elements of stiffness and mass matrices are evaluated directly. But in RayleighRitz and FEMs, the weak form of the governing equations is developed and the boundary conditions are satisfied in the weak form. Generally, a larger number of integrals with increasing amount of differentiation should be done to arrive at the element matrices. Also, the number of degrees of freedom is increased for an acceptable accuracy.
The basic idea of the DQM is the derivative of a function, with respect to a space variable at a given sampling point, which is approximated as a weighted linear sum of the sampling points in the domain of that variable. In order to illustrate the DQ approximation, a function defined on a rectangular domain and is considered. The function values are known or desired on a grid of sampling points in the given domain. According to DQM method, the rth derivative of the function is approximated as what follows:
(40) 
For i=1, 2,…, and r =1,2,…,
From this equation one can deduce that the important components of DQM approximations are the weighting coefficients and the choice of sampling points. In order to determine the weighting coefficients a set of test functions are used in Eq. (40). The weighting coefficients for the firstorder derivatives in  direction are thus determined as what follows [58]:
(41) 
where
(42) 
The weighting coefficients of the secondorder derivative can be obtained as the matrix form [58]:
(43) 
In a similar manner, the weighting coefficients for the direction can be obtained.
The natural and simplest choice of the grid points is equally spaced points in the direction of the coordinate axes of computational domain. It is demonstrated that nonuniform grid points give a better result with the same number of equally spaced grid points [58]. It is shown [59] that one of the best options for obtaining grid points is Chebyshev–Gauss–Lobatto quadrature points:
(44) 
For i = 1, 2,…, ; j = 1,2,…,
Using the geometrical periodicity of the plate, the displacement components for the free vibration analysis are represented as the following:
(45) 
Where m (=0,1,…, ) is the circumferential wave number; is the natural frequency and i (= ) is the imaginary number. It is obvious that m=0 means axisymmetric vibration. At this stage the DQ rules are employed to discretize the free vibration equations and the related boundary conditions. Substituting the displacement components for Eq. (45) and then using the DQ rules for the spatial derivatives, the discretized form of the equations of motion at each domain grid point with (j = 2,3,…, N_{r})1) and ( k =2,3,…, N_{z}1) is obtained as Eq. (32):
(46) 
Eq. (33):
(47) 
Eq. (34):
(48) 
Where , and , are the firstorder and secondorder DQ weighting coefficients in the r and z directions, respectively. In addition, U_{rmjk}, U_{θmjk} and U_{zmjk} represent the displacement components of the node defined by r=r_{j} and Z=Z_{k }. Also, N_{r} and N_{z }represent the total number of nodes through the radial and thickness of the plate, respectively. In a similar manner the boundary conditions can be discretized. For this purpose, using Eq. (45) and the DQ discretization rules for spatial derivatives, the boundary conditions at z = h/2 and h/2 , Eq. (35) become,
at z = h/2


(49) 
Eq. (36):
at z =h/2
(50) 
Where k = 1 at z =h/2 and k =N_{z} at z =h/2, and j = 1, 2, . . ., .
The boundary conditions at r = b and a state in equations (3739) become,
Simply supported (S):
(51) 
Clamped (C):
(52) 
Free (F):
(53) 
In the above equations k = 2, . . ., N_{z}1; also j = 1 at r = b and j = N_{r} at r = a. In order to carry out the eigenvalue analysis, the domain and boundary nodal displacements are separated. In vector forms, they are denoted as {d} and {b}, respectively. Based on this definition, the discretized form of the equations of motion and the related boundary conditions is represented in the matrix form as Equations of motion (4648):
(54) 
Boundary conditions (49, 50) and (5153):
(55) 
Eliminating the boundary degrees of freedom in Eq. (54) using Eq. (55), this equation is obtained as follows:
(56) 
Where . The above eigenvalue system of equations can be solved to find the natural frequencies and mode shapes of the plate.
5. Numerical Results and Discussion
5.1. Convergence and Comparison Studies
Due to lack of appropriate results for free vibration of Continuously Graded Carbon NanotubeReinforced (CGCNTR) annular plates resting on a twoparameter foundation for direct comparison, validation of the presented formulation is conducted in two ways. Firstly, the results are compared with those of 1D conventional functionally graded annular plates, and then, the results of the presented formulations are given in the form of convergence studies with respect to N_{x} and N_{y}, the number of discrete points distributed along the thickness and width of the plate, respectively.
As a first example, it is assumed that the material properties have the following exponential distributions in the thickness direction of the plate:
(57) 
Ceramic (Alumina, Al_{2}O_{3}):
Where the superscript C refers to the material properties of the bottom surface and is the material property graded index.
In Table 1, the first nondimensional natural frequency parameters for the simply supportedclamped FG annular plates are compared with those of Nie and Zhong [35] and Dong [36].
As the second example, the first three nondimensional frequencies for FG annular plates with clamped inner and outer edges for different circumferential wave number (m) are compared with those of the threedimensional elasticity solution of Nie and Zhong [45] in Table 2.
As the third example, based on the power law distribution, the Young’s modulus E and the mass density are assumed to be in terms of a power law distribution as follows:
(58) 

(59) 
Where h is the thickness of the plate and is the power law index which takes values greater than or equal to zero. Subscripts M and C refer to the metal and ceramic constituents which denote the material property of the top and bottom surface of the plate, respectively. The material properties are as follows:
Ceramic(Alumina, Al_{2}O_{3}):
In Tables 3 and 4, the results for FG annular plates are compared with those of Dong [36] for different values of the power law index and circumferential wave number (m). According to the data presented in the abovementioned tables, excellent solution agreements are observed between the present method and those of the other methods.
Based on the above studies, a numerical value of N_{r} = N_{r} =17 is used for the next studies.
After demonstrating the convergence and accuracy of the method, the results of parametric studies for 3D vibration analysis of Continuously Graded Carbon NanotubeReinforced (CGCNTR) annular plates resting on an elastic foundation for different CNTs distributions and various thickness to outer radius ratio (h/a) and different combinations of free, simply supported and clamped boundary conditions at the circular edges, are computed. The nondimensional natural frequency, Winkler and shearing layer elastic coefficients are as follows:
Table 1. Convergence results of the first nondimensional natural frequency parameters for FG annular plates with simply supported (r=b) and clamped (r=a) edges 

N_{r}=N_{z} 
wave number (m) 
λ 



1 
5 
10 
15 

7 
0 
0.1886 
0.1331 
0.0784 
0.0529 
9 

0.1873 
0.1318 
0.0783 
0.0536 
11 

0.1872 
0.1316 
0.0782 
0.0534 
13 17 

0.1872 0.1870 
0.1314 0.1315 
0.0782 0.0781 
0.0535 0.0534 
Ref. [36] 

0.1871 
0.1315 
0.0780 
0.0536 
Ref. [35] 

0.1936 
 
 
 
7 
1 
0.1801 
0.1313 
0.0733 
0.0475 
9 

0.1972 
0.1394 
0.0809 
0.0576 
11 

0.1990 
0.1401 
0.0821 
0.0579 
13 17 

0.1990 0.1993 
0.1401 0.1402 
0.0852 0.0842 
0.0581 0.0582 
Ref. [36] 

0.1994 
0.1402 
0.0840 
0.0582 
Ref. [35] 

0.2050 
 
 
 
7 
2 
0.2744 
0.1955 
0.1227 
0.0851 
9 

0.2748 
0.1968 
0.1202 
0.0842 
11 

0.2785 
0.1973 
0.1201 
0.0832 
13 17 

0.2783 0.2782 
0.1969 0.1967 
0.1201 0.1187 
0.0831 0.0823 
Ref. [36] 

0.2781 
0.1967 
0.1184 
0.0820 
Ref. [35] 

0.2684 
 
 
 
7 
3 
0.3831 
0.277 
0.1715 
0.1188 
9 

0.3824 
0.2765 
0.1697 
0.1184 
11 

0.3824 
0.2757 
0.1696 
0.1180 
13 17 

0.3819 0.3819 
0.2757 0.2752 
0.1692 0.1692 
0.1181 0.1182 
Ref. [36] 

0.3819 
0.2751 
0.1693 
0.1182 
Ref. [35] 

 
 
 
 
(60) 

(61) 
where and are mechanical properties of matrix. In this work, Poly (methyl methacrylate), referred to as PMMA, is selected for the matrix, and the material properties of which are assumed to be, and [60,61]. The (10,10)
SWCNTs are selected as reinforcements. The material properties of the (10,10) SWCNTs used here from Refs. [17,60,61] are as follows (at room temperature, 300 K) :
Table 2. Convergence results of the first three nondimensional frequencies for FG annular plates with clamped  clamped edges
wave number (m) 
Number of the discrete points along the radial and thickness directions while using DQM 

7 
9 
11 
13 
17 
Ref.[45] 
Ansys[1] 

0 
0.094 
0.0856 
0.0816 
0.0801 
0.0806 
0.0807 
0.0810 
1 
0.1006 
0.0896 
0.0844 
0.0831 
0.0838 
0.0837 
0.0839 
2 
0.1147 
0.1027 
0.0977 
0.0955 
0.0961 
0.0961 
0.0963 
1. Ref. [45]
Table 3. Convergence study of the first five nondimensional natural frequency parameters for free vibration of a clampedclamped FG annular plate. 

N_{r}=N_{z} 
wave number (m) 

7 9 11 13 17 Ref. [36] 7 9 11 13 17 Ref. [36] 7 9 11 13 17 Ref. [36] 7 9 11 13 17 Ref. [36] 
0
1
2
3 
8.177 8.201 8.208 8.210 8.213 8.214 8.303 8.322 8.327 8.329 8.332 8.333 8.849 8.861 8.863 8.865 8.868 8.869 9.901 9.906 9.919 9.921 9.923 9.924 
13.912 13.875 13.867 13.870 13.872 13.872 9.696 9.689 9.688 9.688 9.689 9.689 11.160 11.147 11.146 11.145 11.145 11.145 12.693 12.681 12.670 12.673 12.673 12.672 
15.516 15.511 15.511 15.511 15.515 15.514 13.803 13.769 13.767 13.765 13.766 13.766 13.842 13.814 13.812 13.810 13.811 13.810 14.423 14.399 14.402 14.407 14.407 14.407 
19.446 19.481 19.484 19.485 19.485 19.485 14.885 14.853 14.851 14.850 14.849 14.850 15.638 15.615 15.615 15.614 15.614 15.615 16.390 16.422 16.451 16.453 16.456 16.455 
20.108 20.158 20.162 20.164 20.166 20.167 15.546 15.533 15.533 15.533 15.536 15.535 16.561 16.548 16.549 16.549 16.550 16.550 17.699 17.714 17.718 17.720 17.721 17.721 
Table 4.Convergence study of the first five nondimensional natural frequency parameters for free vibration of a clampedclamped FG annular plate.
N_{r}=N_{z} 
wave number (m) 

7 9 11 13 17 Ref. [36] 7 9 11 13 17 Ref. [36] 7 9 11 13 17 Ref. [36] 7 9 11 13 17 Ref. [36] 
0
1
2
3 
10.063 10.087 10.094 10.096 10.098 10.099 10.237 10.256 10.261 10.263 10.267 10.266 10.917 10.929 10.931 10.933 10.936 10.937 12.178 12.228 12.241 12.243 12.247 12.246 
18.379 18.342 18.333 18.336 18.338 18.338 12.343 12.336 12.335 12.335 12.336 12.336 14.407 14.394 14.393 14.392 14.392 14.392 16.685 16.673 16.662 16.665 16.664 16.664 
19.726 19.720 19.721 19.721 19.723 19.724 18.229 18.195 18.193 18.191 18.191 18.192 18.479 18.451 18.449 18.447 18.448 18.448 19.325 19.301 19.304 19.309 19.310 19.310 
24.456 24.421 24.424 24.425 24.427 24.426 18.615 18.583 18.580 18.578 18.579 18.578 18.514 18.490 19.490 19.489 19.489 19.490 20.630 20.662 20.691 20.693 20.694 20.695 
25.726 25.786 25.790 25.792 25.794 25.794 19.676 19.653 19.649 19.649 19.651 19.651 20.715 20.702 20.703 20.703 20.704 20.704 22.402 22.413 22.418 22.420 22.421 22.421 
5.2. Parametric Studies
After demonstrating the convergence and accuracy of the method, the results of parametric studies for 3D vibration analysis of elastically supported thick CGCNTR annular plates reinforced by oriented CNTs for different CNTs distributions and various innertoouter radius ratio (b/a) and different combinations of free, simply supported, and clamped boundary conditions at the circular edges, are computed. Figures 4, 5 and 6 show the effect of the CNTs volume fraction on CGCNTRV to CNTRUD fundamental frequency ratio, , of the nanocomposite annular plates for different values of a/h ratio and boundary condtions including ClampedClamped, Simply SupportedClamped, and FreeClamped at the circular edges. Three different values of the CNTs volume fraction and 0.28 are taken into account. Correspondingly, the CNTs mass fractions are and 0.321, respectively, by taking the density of CNT and the density of matrix in Eq. (28). It can be seen that the discrepancies between the frequencies for the plates with continuously graded and uniform distribution of CNTs increase with the increase in the CNTs volume fraction . This figure also shows/In this figure,it is also shown that the discrepancies between the frequencies decrease with the increase in the a/h ratio.
Figure 4. Variation of the ratio of the ClampedClamped nanocomposite annular plates resting on a twoparameter elastic foundation for different values of a/h ratio and CNTs volume fraction (Kg=10, Kw=100, b/a=0.4) 
The variation of ratio of the nanocomposite annular plates with b/a and h/a ratios is shown in Figures 7, 8 and 9. As it is observed, the ratio decreases rapidly with the increase in b/a ratio and then remains almost unaltered for the b/a>7.
Figure 5.Variation of the ratio of the ClampedSimply supported nanocomposite annular plates resting on a twoparameter elastic foundation for different values of a/h ratio and CNTs volume fraction (Kg=10, Kw=100, b/a=0.4)

Figure 6. Variation of the ratio of the ClampedFree nanocomposite annular plates resting on a twoparameter elastic foundation for different values of a/h ratio and CNTs volume fraction (Kg=10, Kw=100, b/a=0.4)

Figure 7. Variation of the ratio of the ClampedClamped nanocomposite annular plates resting on a twoparameter elastic foundation for different values of b/a and h/a ratios (Kg=10, Kw=100, Vf*=0.12)

Figure 8. Variation of the ratio of the ClampedSimply supported nanocomposite annular plates resting on a twoparameter elastic foundation for different values of b/a and h/a ratios (Kg=10, Kw=100, Vf*=0.12)

Figure 9. Variation of the ratio of the ClampedFree nanocomposite annular plates resting on a twoparameter elastic foundation for different values of b/a and h/a ratios (Kg=10, Kw=100, Vf*=0.12) 
Figure 10. Variation of the e ratio of the nanocomposite annular plates resting on a twoparameter elastic foundation for different boundary conditions and different values of b/a ratio (Kg=10, Kw=100, Vf*=0.12, h/a=0.2)

It is also observed that when the h/a ratio becomes bigger, the discrepancies between the frequencies for CGCNTRV and CNTRUD annualr plates become larger. In Fig. 10, the effect of various boundary conditions on the ratio of the nanocomposite annular plates with for different values of b/a ratio is depicted. It can also be inferred from Fig. 10 that the ClampedClamped CGCNTR annular plate has the highest , whereas the FreeClamped one has the lowest ratio values, implying that the discrepancies between the frequencies of CGCNTR and CNTRUD annular plate with greater supporting rigidity are lower.In addition, Fig. 10 reveals that effects of the boundary conditions on the ratio diminish as b/a ratio increases.
In Fig. 11 the effects of variation of wave number (m) on the frequency parameters of ClampedClamped CGCNTRV annular plate with for different values of a/h ratio are demonstrated. According to Fig. 11, the general behavior of the frequency parameters of CGCNTR annular plate for all a/h ratios is that the frequency parameters converge only in the range beyond that of the fundamental frequency parameters. This means that the effects of a/h ratio are more prominent at low wave numbers, particularly those in the range beyond that of the fundamental frequency parameters, rather than at high wave numbers. As it is shown in Fig. 11, when the wave number increases the discrepancies between the frequency parameters for the different values of a/h ratio become larger. This behavior is also observed at other boundary conditions that again are not shown here for the sake of brevity. Fig. 12 shows the effects of variation of the Winkler elastic coefficient on the fundamental frequency parameters of the ClampedClamped nanocomposite annular plate and on different values of shearing layer elastic coefficient. It is clear that with increasing the elastic coefficients of the foundation, the frequency parameters increase to some limited values. It is observed for the large values of Winkler elastic coefficient, the shearing layer elastic coefficient has less effect and the results become independent of it. The influence of shearing layer elastic coefficient on the fundamental frequency parameters is shown in Fig. 13. One can see that the Winkler elastic coefficient has little effect on the fundamental frequency parameters at different values of shearing layer elastic coefficient. This behavior is also observed at other boundary conditions, but, for the sake of brevity, only this type of the boundary condition is considered.
In Fig. 11 the effects of variation of wave number (m) on the frequency parameters of ClampedClamped CGCNTRV annular plate with for different values of a/h ratio are demonstrated. According to Fig. 11, the general behavior of the frequency parameters of CGCNTR annular plate for all a/h ratios is that the frequency parameters converge only in the range beyond that of the fundamental frequency parameters. This means that the effects of a/h ratio are more prominent at low wave numbers, particularly those in the range beyond that of the fundamental frequency parameters, rather than at high wave numbers. As it is shown in Fig. 11, when the wave number increases the discrepancies between the frequency parameters for the different values of a/h ratio become larger. This behavior is also observed at other boundary conditions that again are not shown here for the sake of brevity. Fig. 12 shows the effects of variation of the Winkler elastic coefficient on the fundamental frequency parameters of the ClampedClamped nanocomposite annular plate and on different values of shearing layer elastic coefficient. It is clear that with increasing the elastic coefficients of the foundation, the frequency parameters increase to some limited values. It is observed for the large values of Winkler elastic coefficient, the shearing layer elastic coefficient has less effect and the results become independent of it. The influence of shearing layer elastic coefficient on the fundamental frequency parameters is shown in Fig. 13. One can see that the Winkler elastic coefficient has little effect on the fundamental frequency parameters at different values of shearing layer elastic coefficient. This behavior is also observed at other boundary conditions, but, for the sake of brevity, only this type of the boundary condition is considered.
6. Conclusion Remarks
In the present work, differential quadrature method is employed to obtain a highly accurate semianalytical solution for free vibration of nanocomposite annular plates resting on a twoparameter elastic foundation under various boundary conditions. The study is carried out based on the threedimensional, linear and small strain elasticity theory. The volume fractions of oriented, straight SingleWalled Carbon Nanotubes (SWCNTs) are assumed to be graded in the thickness direction. The EshelbyMoriTanaka approach is used to estimate the effective constitutive law of the elastic isotropic medium (matrix) with oriented, straight CNTs. The impacts of the volume fractions of oriented CNTs, different CNTs distributions, geometrical parameters and elastic coefficients of foundation on the vibrational characteristics of elastically supported thick annular plates are investigated. The following conclusions can be made from this study:
Figure 11. Variation of the circumferential wave number (m) of the ClampedClamped nanocomposite annular plates resting on a twoparameter elastic foundation for different values of a/h ratio (Kg=10, Kw=100, Vf*=0.12, b/a=0.4) 
Figure 12. Variations of fundamental frequency parameters of the ClampedClamped nanocomposite annular plate resting on a twoparameter elastic foundation with Winkler and different shearing layer elastic coefficients (a/h=2, Vf*=0.12, b/a=0.4) 
Figure 13. Variations of fundamental frequency parameters of the ClampedClamped nanocomposite annular plate versus the shearing layer elastic coefficient for different Winkler elastic coefficients (a/h=2, Vf*=0.12, b/a=0.4) 
Refferences
[1] Esawi AMK and Farag MM. Carbon nanotube reinforced composites: potential and current challenges. Mater Des 2007; 28: 23942401.
[2] Salvetat D and Rubio A. Mechanical properties of carbon nanotubes: a fiber digest for beginners. Carbon 2002; 40: 17291734.
[3] Endo M, Hayashi T, Kim YA, Terrones M and Dresselhaus MS. Applications of carbon nanotubes in the twentyfirst century. Phil Trans R SocLondA 2004; 362: 22232238.
[4] Wernik JM and Meguid SA. Multiscale modeling of the nonlinear response of nanoreinforced polymers. ActaMech 2011; 217: 116.
[5] Thostenson ET, Ren ZF and Chou TW. Advances in the Science and Technology of Carbon Nanotubes and their Composites. A Review Compos Sci Technol 2001; 61: 18991912.
[6] Moniruzzaman M and Winey KI. Polymer nanocomposites containing carbon nanotubes. Macromolecules 2006; 39: 51945205.
[7] Valter B, Ram MK and Nicolini C. Synthesis of multiwalled carbon nanotubes and poly (oanisidine) nanocomposite material: fabrication and characterization od its langmuirschaefer films, Langmuir. 2002; 18: 15351541.
[8] Qian D, Dickey EC, Andrews R and Rantell T. Load transfer and deformation mechanisms in carbon nanotubepolystyrene composites. ApplPhys Lett 2000; 76: 28682870.
[9] Yokozeki T, Iwahori Y and Ishiwata S. Matrix cracking behaviors in carbon fiber/epoxy laminates filled with cupstacked carbon nanotubes (CSCNTs). Composites Part A 2007; 38: 917924.
[10] Hu N, Fukunaga H, Lu C, Kameyama M and Yan B. Prediction of elastic properties of carbon nanotube reinforced composites. Proc R Soc A 2005; 461: 16851910.
[11] Fidelus JD, Wiesel E, Gojny FH, Schulte K and Wagner HD. Thermomechanical properties of randomly oriented carbon/epoxy nanocomposites. Composites Part A 2005; 36: 15551561.
[12] Bonnet P, Sireude D, Garnier B and Chauvet O. Thermal properties and percolation in carbon nanotubepolymer composites. Appl Phys 2007; 91: 20192030.
[13] Han Y and Elliott J. Molecular dynamics simulations of the elastic properties of polymer/carbon nanotube composites. Comput Mater Sci 2007; 39: 315323.
[14] Odegard GM, Gates TS, Wise KE, Park C and Siochi EJ. Constitutive modelling of nanotube reinforced polymer composites. Compos SciTechnol 2003; 63: 16711687.
[15] Shen HS. Nonlinear bending of functionally graded carbon nanotubereinforced composite plates in thermal environments. Compos Struct 2009; 91: 919.
[16] Shen HS and Zhu ZH. Buckling and postbucklingbehavior of functionally graded nanotubereinforced composite plates in thermal environments. Comput Mater Continua 2010; 18: 155182.
[17] Shen HS and Zhang CL. Thermal buckling and postbucklingbehavior of functionally graded carbon nanotubereinforced composite plates. Mater Des 2010; 31: 34033411.
[18] Ke LL, Yang J and Kitipornchai S. Nonlinear free vibration of functionally graded carbon nanotube reinforced composite beams. Compos Struct2010; 92: 676683.
[19] Xiang Y, Kitipornchai S and Liew KM. Buckling and vibration of thick laminates on Pasternak foundations. Eng Mech ASCE 1996; 122: 5463.
[20] Xiang Y, Wang CM and Kitipornchai S. Exact vibration solution for initially stressed Mindlin plates on Pasternak foundations. Int J MechSci1994; 36: 311316.
[21] Wang CM, Kitipornchai S and Xiang Y. Relationships between buckling loads of Kirchhoff, Mindlin, and Reddy polygonal plates. Eng Mech ASCE 1997; 123: 11341137.
[22] Gupta US, Lal R and Sagar R. Effect of an elastic foundation on axisymmetric vibrations of polar orthotropic Mindlin circular plates. Indian J Pure Appl Math 1994; 25: 13171326.
[23] Ju F and Lee HPKH. Free vibration of plates with stepped variations in thickness on nonhomogeneous elastic foundations. Sound Vib 1995; 183: 533545.
[24] Gupta US, Lal R and Jain SK. Effect of elastic foundation on axisymmetric vibrations of polar orthotropic circular plates of variable thickness. Sound Vib 1990; 139: 503513.
[25] Gupta US and Ansari AH. Effect of elastic foundation on axisymmetric vibrations of polar orthotropic linearly tapered circular plates. Sound Vib 2002; 254: 411426.
[26] Laura PAA and Gutierrez RH. Free vibrations of a solid circular plate of linearly varying thickness and attached to Winkler foundation. Sound Vib 1991; 144: 149161.
[27] Matsunaga H. Free vibration and stability of functionally graded plates according to a 2D higherorder deformation theory. J Compos Struct 2008; 82: 499512.
[28] Zhou D, Cheung YK, Lo SH and Au FTK. Threedimensional vibration analysis of rectangular thick plates on Pasternak foundation. Int J Numer Methods Eng 2004; 59: 13131334.
[29] Matsunaga H. Vibration and stability of thick plates on elastic foundations. Eng Mech ASCE 2000; 126: 2734.
[30] Tahouneh V and Yas MH. 3D free vibration analysis of thick functionally graded annular sector plates on Pasternak elastic foundation via 2D differential quadrature method. Acta Mech2012; 223: 18791897.
[31] Tahouneh V and Yas MH. Semianalytical solution for threedimensional vibration analysis of thick multidirectional functionally graded annular sector plates under various boundary conditions. J EngMech ASCE 2013; In Press.
[32] Tahouneh V., Yas M.H., Tourang H., Kabirian M. Semianalytical solution for three dimensional vibration of thick continuous grading fiber reinforced (CGFR) annular plates on Pasternak elastic foundations with arbitrary boundary conditions on their circularedges. Meccanica,2013,48:13131336.
[33] Tahouneh V. Free vibration analysis of thick CGFR annular sector plates resting on elastic foundations. Struct Eng Mech 2014, Vol. 50, No. 6, 773796.
[34] Jam J. E., Kamarian S., Pourasghar A., Seidi J. Free Vibrations of ThreeParameter Functionally Graded Plates Resting on Pasternak Foundations, Solid Mechanics 2012, Vol. 4, No. 1, pp. 5974 (IAU, Arak BranchIran)
[35] Nie GJ and Zhong Z. Semianalytical solution for threedimensional vibration of functionally graded circular plates. Comput Methods Appl Mech Eng 2007; 196: 49014910.
[36] Dong CY. Threedimensional free vibration analysis of functionally graded annular plates using the Chebyshev–Ritz method. Mater Des 2008; 29: 15181525.
[37] Cheng ZQ and Batra RC. Exact correspondence between eigenvalues of membranes and functionally graded simply supported polygonal plates. Sound Vib 2000; 229: 879895.
[38] Malekzadeh P. Threedimensional free vibrations analysis of thick functionally graded plates on elastic foundations. Compos Struct 2008; 89: 367373.
[39] Bellman R and Casti J. Differential quadrature and long term integration. Math Anal Appl1971; 34: 235238.
[40] Liu FL and Liew KM. Free vibration analysis of Mindlin sector plates numerical solutions by differential quadrature method. Comput MethodsApplMechEng 1999; 177: 77–92.
[41] Liew KM and Liu FL. Differential quadrature method for vibration analysis of shear deformable annular sector plates. Sound Vib 2000; 230: 335356.
[42] Wang X and Wang Y. Free vibration analyses of thin sector plates by the new version of differential quadrature method. Comput Methods ApplMechEng 2004; 193: 39573971.
[43] Liew KM, Han JB, Xiao ZM and Du H. Differential quadrature method for Mindlin plates on Winkler foundation. Int J MechSci 1996; 38: 405421.
[44] Gupta US, Lal R and Sharma S. Vibration analysis of nonhomogeneous circular plate of nonlinear thickness variation by differential quadrature method. Sound Vib 2006; 298: 892906.
[45] Nie GJ and Zhong Z. Dynamic analysis of multidirectional functionally graded annular plates. Appl Math Model 2010; 34: 608616.
[46] Shu C. Differential quadrature and its application in engineering. Springer, 2000; Berlin.
[47] Shi DL, Feng XQ, Huang YY, Hwang KC and Gao H. The effect of nanotube waviness and agglomeration on the elastic property of carbon nanotubereinforced composites. Eng Mater Technol 2004; 126: 250257.
[48] Mori T and Tanaka K. Average Stress in Matrix and Average Elastic Energy of Materials WithMisfitting Inclusions.ActaMetall 1973; 21: 571574.
[49] Giordano S, Palla PL and Colombo L. Nonlinear elasticity of composite materials, Landau coefficients in dispersions of spherical and cylindrical inclusions. EurPhys J B 2009; 68: 89101.
[50] Benveniste Y. A new approach to the application of MoriTanaka’s theory in composite materials. Mech Mater 1987; 6: 147157.
[51] Chen CH and Cheng CH. Effective elastic moduli of misoriented shortfiber composites. Int J Solids Struct 1996; 33: 25192539.
[52] Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc R Soc A 1957; 241: 376396.
[53] Eshelby JD. The elastic field outside an ellipsoidal inclusion. Proc R Soc A 1959; 252: 561–569.
[54] Formica G, Lacarbonara W and Alessi R. Vibrations of carbon nanotubereinforced composites. Sound Vib 2010; 329: 18751889.
[55] Mura T. Micromechanics of defects in solids. The Hague: MartinusNijhoff; 1982.
[56] Shen HS. Postbuckling of nanotubereinforced composite cylindrical shells in thermal environments. Part I: Axiallyloaded shells, Compos Struct2011; 93: 20962108.
[57] Fung YC and Tong P. Classical and Computational Solid Mechanics. World Scientific Publishing Co, Singapore, 2001.
[58] Bert CW and Malik M. Differential quadrature method in computational mechanics: a review. Appl Mech Rev 1996; 49: 128.
[59] Shu C and Wang CM. Treatment of mixed and nonuniform boundary conditions in GDQ vibration analysis of rectangular plates.EngStruct1999; 21: 125134.
[60] Shen HS. Thermal buckling and postbucklingbehavior of functionally graded carbon nanotube reinforced composite cylindrical shells. Compos Part B: Engineering, Article in Press, doi:10.1016/j.compositesb.2011.10.004; 2011.
[61] Wang ZX and Shen HS. Nonlinear vibration of nanotubereinforced composite plates in thermal environments. Comp Mater Sci 2011; 50: 23192330.