The effect of multi-directional nanocomposite materials on the vibrational response of thick shell panels with finite length and rested on two-parameter elastic foundations

The main purpose of this paper is to investigate the effect of bidirectional continuously graded nanocomposite materials on free vibration of thick shell panels rested on elastic foundations. The elastic foundation is considered as a Pasternak model after adding a shear layer to the Winkler model. The panels reinforced by randomly oriented straight single-walled carbon nanotubes are considered. The volume fractions of SWCNTs are assumed to be graded not only in the radial direction, but also in axial direction of the curved panel. This study presents a 2-D sixparameter power-law distribution for CNTs volume fraction of 2-D continuously graded nanocomposite that gives designers a powerful tool for flexible designing of structures under multi-functional requirements. The benefit of using generalized power-law distribution is to illustrate and present useful results arising from symmetric, asymmetric and classic profiles. The material properties are determined in terms of local volume fractions and material properties by Mori–Tanaka scheme. The 2-D differential quadrature method as an efficient numerical tool is used to discretize governing equations and to implement boundary conditions. The fast rate of convergence of the method is shown and results are compared against existing results in literature. Some new results for natural frequencies of the shell are prepared, which include the effects of elastic coefficients of foundation, boundary conditions, material and geometrical parameters. The interesting results indicate that a graded nanocomposite volume fraction in two directions has a higher capability to reduce the natural frequency than conventional 1-D functionally graded nanocomposite materials.


Introduction
Layered composite materials, due to their thermal and mechanical merits compared to single-composed materials, have been widely used for a variety of engineering applications. However, owing to the sharp discontinuity in the material properties at interfaces between two different layers, there may exist stress concentrations causing severe material failure (Weissenbek et al. 1997). Functionally graded materials are heterogeneous composite materials, in which the material properties vary continuously from one interface to the other. The advantage of using these materials is that they can survive in high thermal gradient environment, while maintaining their structural integrity. Typically, an FGM is made of a ceramic and a metal for the purpose of thermal protection against large temperature gradients. The ceramic material provides a high-temperature resistance due to its low thermal conductivity, while the ductile metal constituent prevents fracture due to its greater toughness. FGMs are now developed for general use as structural elements in extremely high temperature environments. A listing of different applications can be found in Forum (1991). Most of the studies on FGMs have been restricted to thermal stress analysis, thermal buckling, fracture mechanics and optimization (Cho and Tinsley oden 2000;Chunyu et al. 2001;Lanhe 2004). Some researches (Loy et al. 1999;Pradhan et al. 2000;Han et al. 2001;Ng et al. 2001) are based on the classical shell theory, i.e., neglecting the effect of transverse shear deformation. The application of this theory to moderately thick or thick shell structures can lead to serious errors. Using the first and higher order shear deformation theories, some modifications are done to include the effects of transverse shear deformation. In this study, the problem formulations were based on the higher order shear deformation shell theories. Yang and Shen (2003) proposed a semi-analytical approach based on Reddy's higher order shear deformation shell theory, for free vibration and dynamic instability of FGM cylindrical panels under combined static and periodic axial forces and thermal loads. Free vibration and stability of functionally graded shallow shells according to a 2-D higher order deformation theory were investigated by Matsunaga (2008). Free vibration analysis of functionally graded curved panels was carried out using a higher order formulation by Pradyumna and Bandyopadhyay (2008). They used a C0 finite element formulation to carry out the analysis. Using a 2-D higher order shear deformation theory, vibration and buckling analyses of simply supported circular cylindrical shells made of functionally graded materials (FGMs) were studied by Matsunaga (2009). He used the method of power series expansion of continuous displacement components to solve the problem. In all of the above studies the variation of the radius through the thickness was not considered and the problem formulations were based on the constant mean radius of curvature.
Two-dimensional theories reduce the dimensions of problems from three to two by introducing some assumptions in mathematical modeling leading to simpler expressions and derivation of solutions. However, these simplifications inherently bring errors and therefore may lead to unreliable results for relatively thick panels. As a result, three-dimensional analysis of panels not only provides realistic results, but also allows further physical insights, which cannot otherwise be predicted by the twodimensional analysis. There are some studies on free vibration analysis of isotropic and composite panels and shells based on the three-dimensional elasticity formulation (Chern and Chao 2000).
Structures resting on elastic foundations with different shapes, sizes, and thickness variations and boundary conditions have been the subject of investigations, and those play an important role in aerospace, marine, civil, mechanical, electronic and nuclear engineering problems. For example, plates and shells are used in various kinds of industrial applications such as the analysis of reinforced concrete pavement of roads, airport runways and foundations of buildings. The Pasternak model (also referred to as the two-parameter model) was widely used to describe the mechanical behavior of the foundation, in which the well-known Winkler model is a special case.
The most serious deficiency of the Winkler foundation model is to have 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 two-parameter model, in which the shear stiffness of the foundation is considered. The evident importance in practical applications, investigations on the dynamic characteristics of FGM plates and panels on elastic foundations are still limited in number. Yas and Tahouneh (2012) investigated the free vibration analysis of thick FG annular plates on elastic foundations via differential quadrature method based on the three-dimensional elasticity theory and Tahouneh and Yas (2012) investigated the free vibration analysis of thick FG annular sector plates on Pasternak elastic foundations using DQM. Tahouneh et al. (2013) studied free vibration characteristics of annular continuous grading fiber reinforced (CGFR) plates resting on elastic foundations using DQM. More recently, (Tahouneh and Naei 2014) achieved the natural frequencies of thick multi-directional functionally graded rectangular plates resting on a two-parameter elastic foundation via 2-D differential quadrature method, The proposed rectangular plates had two opposite edges simply supported, while all possible combinations of free, simply supported and clamped boundary conditions were applied to the other two edges. Farid et al. (2010) presented free vibration analysis of initially stressed thick simply supported functionally graded curved panel resting on twoparameter elastic foundation (Pasternak model), subjected in thermal environment was studied using the three-dimensional elasticity formulation. Tahouneh (2014) investigated free vibration analysis of continuous grading fiber reinforced (CGFR) FG annular sector plates on two-parameter elastic foundations under various boundary conditions, based on the three-dimensional theory of elasticity. The plates with simply supported radial edges and arbitrary boundary conditions on their circular edges were considered.
Recently, nanocomposites have significant importance for engineering applications that require high levels of structural performance and multi-functionality. Carbon nanotubes (CNTs) have demonstrated exceptional mechanical, thermal and electrical properties. These materials are considered as one of the most promising reinforcement materials for high performance structural and multi-functional composites with vast application potentials (Esawi and Farag 2007;Thostenson et al. 2001). Most studies on carbon nanotube-reinforced composites (CNTRCs) have focused on their material properties (Esawi and Farag 2007;Thostenson et al. 2001;Dai 2002;Kang et al. 2006;Lau et al. 2006). Gojny et al. (2005) focused on the evaluation of the different types of the CNTs applied, their influence on the mechanical properties of epoxy-based nanocomposites and the relevance of surface fictionalization. Fidelus et al. (2005) investigated thermo-mechanical properties of epoxy-based nanocomposites based on low weight fraction of randomly oriented single-and multi-walled CNTs. Han and Elliott (2007) determined the elastic modulus of composite structures under CNTs reinforcement by molecular dynamic simulation and investigated the effect of volume fraction of SWNTs on mechanical properties of nanocomposites. Manchado et al. (2005) blended small amounts of arc-SWNT into isotactic polypropylene and observed the modulus increase from 0.85 to 1.19 GPa at 0.75 wt%. In addition, the strength increased from 31 to 36 MPa by 0.5 wt%. Both properties were observed to fall off at higher loading levels. These investigations and (Mokashi et al. 2007;Zhu et al. 2007) have shown that the addition of small amount of carbon nanotube in the matrix can considerably improve the mechanical, electrical and thermal properties of polymeric composites. This behavior, combined with their low density makes them suitable for transport industries, especially for aeronautic and aerospace applications where the reduction of weight is crucial in order to reduce the fuel consumption.
The properties of the CNT-reinforced composites (CNTRCs) depend on a variety of parameters including CNT geometry and the inter-phase between the matrix and CNT. Interfacial bonding in the inter-phase region between embedded CNT and its surrounding polymer is a crucial issue for the load transferring and reinforcement phenomena Shokrieh and Rafiee (2010). The traditional approach to fabricating nanocomposites implies that the nanotube is distributed either uniformly or randomly such that the resulting mechanical, thermal, or physical properties do not vary spatially at the macroscopic level. Experimental and numerical studies concerning CNTRCs have shown that distributing CNTs uniformly as the reinforcements in the matrix can achieve moderate improvement of the mechanical properties only (Seidel and Lagoudas 2006). This is mainly due to the weak interface between the CNTs and the matrix where a significant material property mismatch exists. The concept of FGM can be utilized for the management of a material's microstructure, so that the vibrational behavior of a plate/shell structure reinforced by CNTs can be improved. According to a comprehensive survey of literature, the authors found that there are few research studies on the mechanical behavior of functionally graded CNTRC structures. For the first time, Shen (2009) suggested 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.
Due to intrinsic complexity of the formulations based on the three-dimensional elasticity, powerful numerical methods are needed to solve the governing equations. The differential quadrature method (DQM) is a relatively new numerical technique in structural analysis. A review of the early developments in the differential quadrature method can be found in papers by (Bert and Malik 1997). This paper is motivated by the lack of studies in the technical literature concerning to the three-dimensional vibration analysis of thick bidirectional nanocomposite curved panels resting on a two-parameter elastic foundation reinforced by randomly oriented straight single-walled carbon nanotubes CNTs. To the authors' best knowledge, research on the vibration of thick curved panels reinforced by randomly oriented straight singlewalled carbon nanotubes which are graded in both direction including axial and radial directions has not been seen until now. The volume fractions of randomly straight single-walled carbon nanotubes SWCNTs are assumed to be graded in the thickness and also axial directions of the curved panels. The direct application of CNT properties in micromechanics models for predicting material properties of the nanotube/polymer composite is inappropriate without taking into account the effects associated with the significant size difference between a nanotube and a typical carbon fiber (Odegard et al. 2003). In other words, continuum micromechanics equations cannot capture the scale difference between the nano and micro-levels. In order to overcome this limitation, a virtual equivalent fiber consisting of nanotube and its interphase which is perfectly bonded to surrounding resin is applied.
This study presents a novel 2-D six-parameter powerlaw distribution for CNTs volume fraction of 2-D functionally graded nanocomposite materials that gives designers a powerful tool for flexible designing of structures under multi-functional requirements. Various material profiles along the radial and axial directions are illustrated by using the 2-D power-law distribution. The effective material properties at a point are determined in terms of the local volume fractions and the material properties by the Mori-Tanaka scheme. A sensitivity analysis is performed, and the natural frequencies are calculated for different sets of boundary conditions and different combinations of the geometric, material, and foundation parameters. Therefore, very complex combinations of the material properties, boundary conditions, and foundation stiffness are considered in the present semianalytical solution approach.

Problem description
In this section, a virtual equivalent fiber consisting of a nanotube and its inter-phase which is perfectly bonded to surrounding resin is introduced to obtain the mechanical properties of the carbon nanotube/polymer composite by using the results of multi-scale FEM Shokrieh and Rafiee (2010). The equivalent fiber for SWCNT with chiral index of (10, 10) is a solid cylinder with diameter of 1.424 nm. The inverse of the rule of mixture is used to calculate material properties of equivalent fiber (Tsai et al. 2003): where V EF and V M are longitudinal modulus of equivalent fiber, transverse modulus of equivalent fiber, shear modulus of equivalent fiber, Poisson's ratio of equivalent fiber, longitudinal modulus of composites, transverse modulus of composites, shear modulus of composites, Poisson's ratio of composites, modulus of matrix, shear modulus of matrix, Poisson's ratio of matrix, volume fraction of the equivalent fiber and volume fraction of the matrix, respectively. E LC , G C , and E TC are obtained from multiscale FEM or molecular dynamics (MD) simulations. It should be mentioned that the volume fraction of the equivalent fiber is assumed to be 7.5 % (Shokrieh and Rafiee 2010) and Poly {(mphenylenevinylene)-co-[(2,5dioctoxy-p-phenyle) vinylene]}, referred to as (PmPV), is selected as a matrix material: The material properties adopted for equivalent fiber are (Shokrieh and Rafiee 2010):

Composites reinforced with aligned, straight CNTs
Following the standard MT derivation, one can develop the expression for effective composite stiffness C. This is obtained by using an equivalent fiber having the effective CNT properties in the MT approach which is given as (Shi et al. 2004): where f r and f m are the fiber and matrix volume fractions, respectively. C m is the stiffness tensor of the matrix material; C r is the stiffness tensor of the equivalent fiber; I is the forth order identity tensor and A r is the dilute strainconcentration tensor of the rth phase for the fiber which is given as: where S is Eshelby's tensor, as given by (Eshelby 1957) and (Mura 1982). The terms enclosed by angle brackets in Eq.
(2) represent the average value of the term over all orientations defined by transformation from the local fiber coordinates (O- 1). Assume axis x 2 as the direction along the aligned nanotube. The elastic properties of the nanocomposite are determined from the average strain obtained in the representative volume element. The matrix is assumed to be elastic and isotropic, with Young's modulus E m and Poisson's ratio t m . Each straight CNT is modeled as a long fiber with transversely isotropic elastic properties and has a stiffness matrix given by Eq. (1). Therefore, the composite is also transversely isotropic, with five independent elastic constants. The substitution of nonvanishing components of the Eshelby tensor S for a straight, long fiber along the x 2 -direction (Shi et al. 2004) in Eq. (3) gives the dilute mechanical strain concentration tensor. Then, the substitution of Eq. (3) into Eq. (2) gives the tensor of effective elastic moduli of the composite reinforced by aligned, straight CNTs. The axial and transverse Young's modulus of the composite can be calculated from the Hill's elastic modulus as (Shi et al. 2004): where k, l, m and n are its plane-strain bulk modulus normal to the fiber direction, cross-modulus, transverse shear modulus, axial modulus and axial shear modulus, respectively, and can be found in the Appendix. As mentioned before, the CNTs are transversely isotropic and have a stiffness matrix given below: 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 5 where E L ; E T ; E Z ; G TZ ; G ZL ; G LT ; t LT ; t LZ ; t TZ are material properties of the equivalent fiber which can be determined from the inverse of the rule of mixture.

Composites reinforced with randomly oriented, straight CNTs
The effective properties of composites with randomly oriented non-clustered CNTs, such as in Fig. 1, are studied in this section. The resulting effective properties for the randomly oriented CNT composite are isotropic, despite the CNTs having transversely isotropic effective properties. The orientation of a straight CNT is characterized by two Euler angles a and b, as shown in Fig. 1. When CNTs are completely randomly oriented in the matrix, the composite is then isotropic, and its bulk modulus k and shear modulus G are derived as: where a r ; b r ; d r and g r can be found in the Appendix. The effective Young's modulus E and Poisson's ratio t of the composite is given by: Functionally graded carbon nanotube-reinforced Consider a bidirectional nanocomposite curved panel rested on two-parameter elastic foundations as shown in Fig. 2. A cylindrical coordinate system (r, h, z) is used to label the material point of the panel. The inner surface is continuously in contact with an elastic medium that acts as an elastic foundation represented by the Winkler/Pasternak model with K w and K g that are Winkler and shear coefficients of Pasternak foundation, respectively. One of the well-known power-law distributions which is widely considered by the researchers is three-or four-parameter power-law distribution. The benefit of using such power-law distributions is to illustrate and present useful results arising from symmetric and asymmetric profiles. Consider V c (volume fraction of the CNTs) in form of f(z) 9 g(r), f(z) and g(r) are both the three-parameter power-law distribution. They can be used to illustrate symmetric, asymmetric and classical profiles along the axial and radial directions of the curved panels, respectively. So by considering V c as f(z) 9 g(r), one can present a 2-D six-parameter power-law distribution which is useful to illustrate different types of volume fraction profiles, including classical-classical, symmetric-symmetric and classical-symmetric in both directions.
In order to investigate 3-D dynamic response of thick bidirectional nanocomposite curved panels resting on a two-parameter elastic foundation, it is assumed that the volume fraction of the CNTs follows a 2-D six-parameter power-law distribution: where the radial volume fraction index c r , and the parameters a r , b r and the axial volume fraction index c z , and the parameters a z , b z govern the material variation profile through the radial and axial directions, respectively. The volume fractions V a and V b , which have values that range from 0 to 1, denote the maximum and minimum volume fraction of CNTs. With assumption V b = 1 and V a = 0.3, some material profiles in the radial  Fig. 3, The CNTs volume fraction decreases in the axial direction from 1 at g z = -0.5 to 0 at g z = 0.5. With another choice of the parameters a z , b z , a r and b r , it is possible to obtain volume fraction profiles along the radial and axial directions of the panel as shown in Fig. 4. This figure shows a classical profile versus g r and a symmetric profile versus g z . As observed, volume fraction on the lower edge (g z = -0.5) is the same as that on the upper edge (g z = 0.5). Figure 5 illustrates symmetric profiles through the radial and axial directions obtained by setting b r ¼ b z ¼ 2; and a r ¼ a z ¼ 1. In the following, we have compared several different volume fraction profiles of conventional 1-D and 2-D continuously graded nanocomposite with appropriate choice of the radial and axial parameters of the 2-D six-parameter power-law distribution, as shown in Table 1. It should be noted that the notation classical-symmetric indicates that the 2-D nanocomposite curved panel has classical and symmetric volume fraction profiles in the radial and axial directions, respectively.

The basic formulations
The mechanical constitutive relation that relates the stresses to the strains is as follows: r r r h r z s zh s rz s rh 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 Variations of the symmetric volume fraction profile along the radial and axial directions of the curved panels In the absence of body forces, the governing equations are as follows: Strain-displacement relations are expressed as: where u r , u h and u z are radial, circumferential and axial displacement components, respectively. Upon substitution Eq. (11) into (9) and then into (10), the equations of motion in terms of displacement components with infinitesimal deformations can be written as: The boundary conditions at the concave and convex surfaces, r = r i and r o , respectively, can be described as follows: In this investigation, three different types of classical boundary conditions at edges z = 0 and L z of the panel can be stated as follows: -Simply supported (S): -Clamped (C): -Free (F): For the curved panels with simply supported at one pair of opposite edges, the displacement components can be expanded in terms of trigonometric functions in the direction normal to these edges. In this work, it is assumed that the edges h = 0 and h = U are simply supported.
where m is the circumferential wave number, x is the natural frequency and i (= ffiffiffiffiffiffi ffi À1 p ) is the imaginary number. Substituting for displacement components from Eq. (19) into Eqs. (12,13,14), one gets Equation (12): Equation (13): Equation (14): The geometrical and natural boundary conditions stated in Eq. (15) can also be simplified as follows Equation (15): The boundary conditions stated in Eqs. (16,17,18) can also be simplified; however, for the sale of brevity, they are not shown here.

2-D DQM solution of governing equations
It is difficult to solve analytically the equations of motion, if it is not impossible. Hence, one should use an approximate method to find a solution. Here, the differential quadrature method (DQM) is employed. One can compare DQM solution procedure with the other two widely used traditional methods for plate analysis, i.e., Rayleigh-Ritz 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 elements of stiffness and mass matrices are evaluated directly. But in Rayleigh-Ritz and FEMs, the weak form of the governing equations should be developed and the boundary conditions are satisfied in the weak form.
Generally by doing so larger number of integrals with increasing amount of differentiation should be done to arrive at the element matrices. In addition, the number of degrees of freedom will be 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, is approximated as a weighted linear sum of the sampling points in the domain of that variable. In order to illustrate the DQ approximation, consider a function f ðn; gÞ defined on a rectangular domain 0 n a and 0 g b. Let in the given domain, the function values be known or desired on a grid of sampling points. According to DQM method, the rth derivative of the function f ðn; gÞ can be approximated as: A nðrÞ im f mj for i ¼ 1; 2; . . .; N n and r ¼ 1; 2; . . .; N n À 1 where N n represents the total number of nodes along the ndirection. From this equation one can deduce that the important components of DQM approximations are the weighting coefficients ðA nðrÞ ij Þ and the choice of sampling points. In order to determine the weighting coefficients, a set of test functions should be used in Eq. (24). The weighting coefficients for the first-order derivatives in ndirection are thus determined as (Bellman and Casti 1971): where The weighting coefficients of the second-order derivative can be obtained in the matrix form (Bellman and Casti 1971): In a similar manner, the weighting coefficients for the gdirection 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 was demonstrated that non-uniform grid points gives a better result with the same number of equally spaced grid points (Bellman and Casti 1971). It is shown (Shu and Wang 1999) that one of the best options for obtaining grid points is Chebyshev-Gauss-Lobatto quadrature points: for i ¼ 1; 2; . . .; N n ; ð28:1; 2Þ j ¼ 1; 2; . . .; N g ;where N n and N g are the total number of nodes along the nand g-directions, respectively. At this stage, the DQ method can be applied to discretize the equations of motion (20, 21, 22). As a result, at each domain grid point (r i ,z j ) with i = 2,…, N r -1 and j = 2,…, N z -1, the discretized equations take the following forms. Equation (20): Equation (21): Equation (22): where A r ij , A z ij and B r ij , B z ij are the first-and second-order DQ weighting coefficients in the r-and z-directions, respectively. The DQ method can also be applied to discretize the boundary conditions at r = r i and r o as follows. Equation (23): where i = 1 at r = r i and i = N r at r = r o , and j = 1,2,…,N z ; also d ij is the Kronecker delta. The boundary conditions at z = 0 and L z stated in Eqs. (16,17,18), become Eq. (16): • Simply supported (S): Equation (17): • Clamped (C): Equation (18): • Free (F): In the above equations i = 2,…,N r -1; also j = 1 at z = 0 and j = N z at z = L z .
In order to carry out the eigenvalue analysis, the domain and boundary nodal displacements should be 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 can be represented in the matrix form as:Equations of motion, Eqs. (29, 30, 31): Boundary conditions, Eq. (32) and Eqs. (33,34,35): Eliminating the boundary degrees of freedom in Eq. (36) using Eq. (37), this equation becomes where The above eigenvalue system of equations can be solved to find the natural frequencies and mode shapes of the curved panel.

Numerical results and discussion
Convergence and comparison studies Due to lack of appropriate results for free vibration of CGCNTR cylindrical panels reinforced by oriented CNTs for direct comparison, validation of the presented formulation is conducted in two ways. Firstly, the results are compared with those of FGM composite cylindrical panels and then, the results of the presented formulations are given in the form of convergence studies with respect to N r and N z , the number of discrete points distributed along the radial and axial directions, respectively. To validate the proposed approach its convergence and accuracy are demonstrated via different examples. The obtained natural frequencies based on the three-dimensional elasticity formulation are compared with those of the power series expansion method for both FGM curved panels with and without elastic foundations (Matsunaga 2008;Pradyumna and Bandyopadhyay 2008;Farid et al. 2010). In these studies the material properties of functionally graded materials are assumed as follows: • Metal (Aluminum, Al): Subscripts M and C refer to the metal and ceramic constituents which denote the material properties of the outer and inner surfaces of the panel, respectively. To validate the analysis, results for FGM cylindrical shells are compared with similar ones in the literature, as shown in Table 2. The comparison shows that the present results agreed well with those in the literatures. Besides the fast rate of convergence of the method being quite evident, it is found that only 13 grid points (N r = N z = 13) along the radial and axial directions can yield accurate results. Further validation of the present results for isotropic FGM cylindrical panel is shown in Table 3. In this table, comparison is made for different L z /R and L z /h ratios, and it is observed there is good agreement between the results. As another example, the convergence and accuracy of the method are investigated by evaluating the first three natural frequency parameters of the FG curved panel resting on Pasternak foundations. The non-dimensional forms of the elastic foundation coefficients are defined as K w = k w R/G c and K g = k g /(G c R) in which G c is the shear modulus of elasticity of the ceramic layer. The results are prepared for different thickness-to-mean radius ratios and different values of the DQ grid points along the radial and axial directions, respectively, are shown in Table 4. Also, one can see that excellent agreement exists between the results.

Parametric studies
After demonstrating the convergence and accuracy of the present method, parametric studies for 3-D vibration analysis of thick curved panels resting on a two-parameter elastic foundation reinforced by randomly oriented straight single-walled carbon nanotubes for various CNTs volume fraction distribution, length-to-mean radius ratio, elastic coefficients of foundation and different combinations of free, simply supported and clamped boundary conditions along the axial direction of the curved panel, are computed. The boundary conditions of the panel are specified by the letter symbols, for example, S-C-S-F denotes a curved panel with edges h = 0 and U simply supported (S), edge z = 0 clamped (C) and edge z = L z free (F). The non-dimensional natural frequency, Winkler and shearing layer elastic coefficients are as follows: where q m , E m and G m represent the mass density, Young's modulus and shear modulus of the matrix, respectively. The effect of the Winkler elastic coefficient on the fundamental frequency parameters for different boundary conditions is shown in Figs. 6, 7 and 8. It is observed that the fundamental frequency parameters converge with increasing Winkler elastic coefficient of the foundation. According to these figures, the lowest frequency parameter is obtained by using classical-classical volume fraction profile. On the contrary, the 1-D FG panel with symmetric volume fraction profile has the maximum value of the frequency parameter. Therefore, a graded CNTs volume fraction in two directions has higher capabilities to reduce the frequency parameter than conventional 1-D nanocomposite. It is also observed from Figs. 6, 7 and 8, for the large values of Winkler elastic coefficient, the shearing layer elastic coefficient has less effect and the results become independent of it, in other words the non-dimensional natural frequencies converge with increasing Winkler foundation stiffness.
The influence of shearing layer elastic coefficient on the non-dimensional natural frequency for S-C-S-C, S-C-S-S and S-F-S-F bidirectional nanocomposite curved panel resting on a two-parameter elastic foundation, is shown in Figs. 9, 10 and 11. It is observed that the variation of Winkler elastic coefficient has little effect on the non-dimensional natural frequency at different values of shearing layer elastic coefficient. It is clear that with increasing the shearing layer elastic coefficient of the foundation, the frequency parameters increase to some limit values and for  the large values of shearing layer elastic coefficient, the frequency parameters become independent of it. The variations of fundamental frequency parameters of bidirectional nanocomposite curved panels resting on an elastic foundation with length-to-mean radius ratio (L z /R) for different types of volume fraction profiles are depicted in Figs. 12, 13 and 14. It can also be inferred from these figures that the frequency is greatly influenced in that fundamental frequency parameter decreases steadily as length-to-mean radius ratio (L z /R) becomes larger and remains almost unaltered for the large values of length-tomean radius ratio. As can be seen from this figure, for the all length-to-mean radius ratio (L z /R), classical-classical volume fraction profile has the lowest frequencies followed by classical-symmetric, classical, symmetric-symmetric and symmetric profiles.
The variations of fundamental frequency parameters of bidirectional nanocomposite curved panels with length-to- Fig. 8 Variations of fundamental frequency parameters of a bidirectional S-F-S-F nanocomposite curved panels resting on a twoparameter elastic foundation with Winkler elastic coefficient for different volume fraction profiles (K g = 100, R/h = L z /R = 3.5, c r = 2, U = 135°) Fig. 9 Variations of fundamental frequency parameters of a bidirectional S-C-S-C nanocomposite curved panels resting on a twoparameter elastic foundation with the shearing layer elastic coefficient (R/h = L z /R = 3.5, c r = c z = 2, a r = a z = 0, U = 135°) Fig. 10 Variations of fundamental frequency parameters of a bidirectional S-C-S-S nanocomposite curved panels resting on a twoparameter elastic foundation with the shearing layer elastic coefficient (R/h = L z /R = 3.5, c r = c z = 2, a r = a z = 0, U = 135°) Fig. 11 Variations of fundamental frequency parameters of a bidirectional S-F-S-F nanocomposite curved panels resting on a twoparameter elastic foundation with the shearing layer elastic coefficient (R/h = L z /R = 3.5, c r = c z = 2, a r = a z = 0, U = 135°) mean radius ratio (L z /R), and the volume fraction index through the radial direction of the panels for S-F-S-F boundary conditions are shown in Fig. 15, by considering ða r ¼ a z ¼ 0; c z ¼ 2; K w ¼ K g ¼ 100Þ for classical-classical 2-D nanocomposite curved panels. Confirming the effect of length-to-mean radius ratio on the natural frequency already shown in the Figs. 12, 13 and 14, it is found that the frequency parameter decreases by increasing the radial volume fraction index ðc r Þ. This behavior is also observed for other boundary conditions, not shown here for brevity.

Conclusion remarks
In this research work, free vibration of thick bidirectional nanocomposite curved panels resting on a two-parameter elastic is investigated based on three-dimensional theory of Fig. 12 Variations of fundamental frequency parameters of twodimensional continuously graded S-C-S-C nanocomposite curved panels resting on an elastic foundation with L z /R ratio for different volume fraction profiles (K w = K g = 100, R/h = 3.5, c z = 2, U = 135°) Fig. 13 Variations of fundamental frequency parameters of twodimensional continuously graded S-C-S-S nanocomposite curved panels resting on an elastic foundation with L z /R ratio for different volume fraction profiles (K w = K g = 100, R/h = 3.5, c z = 2, U = 135°) Fig. 14 Variations of fundamental frequency parameters of twodimensional continuously graded S-F-S-F nanocomposite curved panels resting on an elastic foundation with L z /R ratio for different volume fraction profiles (K w = K g = 100, R/h = 3.5, c z = 2, U = 135°) Fig. 15 Variations of fundamental frequency parameters of bidirectional nanocomposite curved panels with length-to-mean radius ratio (L z /R), and the volume fraction index through the radial direction of the panels for S-F-S-F boundary condition (K w = K g = 100, R/ h = 3.5, c z = 2, a r = a z = 0, U = 135°) elasticity. The elastic foundation is considered as a Pasternak model with adding a shear layer to the Winkler model. Three complicated equations of motion for the curved panel under consideration are semi-analytically solved by using 2-D differential quadrature method. Using the 2-D differential quadrature method along the radial and axial directions, allows one to deal with curved panel with arbitrary thickness distribution of material properties and also to implement the effects of the elastic foundations as a boundary condition on the lower surface of the curved panel efficiently and in an exact manner. The volume fractions of randomly oriented straight single-walled carbon nanotubes (SWCNTs) are assumed to be graded not only in the radial direction, but also in axial direction of the curved panel. The direct application of CNTs properties in micromechanics models for predicting material properties of the nanotube/polymer composite is inappropriate without taking into account the effects associated with the significant size difference between a nanotube and a typical carbon fiber. In other words, continuum micromechanics equations cannot capture the scale difference between the nano and micro-levels. In order to overcome this limitation, a virtual equivalent fiber consisting of nanotube and its inter-phase which is perfectly bonded to surrounding resin is applied. In this research work, an equivalent continuum model based on the Eshelby-Mori-Tanaka approach is employed to estimate the effective constitutive law of the elastic isotropic medium (matrix) with oriented straight CNTs. The effects of elastic foundation stiffness parameters, various geometrical parameters on the vibration characteristics of CGCNTR curved panel, are investigated, and also, different types of volume fraction profiles along the radial and axial directions of the panels and elastic coefficients of foundation of bidirectional curved panels resting on a two-parameter elastic foundation are studied. Moreover, vibration behavior of 2-D continuously graded nanocomposite panels are compared with conventional one-dimensional nanocomposite panels. From this study, some conclusions can be made: • 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, in other words the non-dimensional natural frequencies converge with increasing Winkler foundation stiffness. • The results show that the variation of Winkler elastic coefficient has little effect on the non-dimensional natural frequency at different values of shearing layer elastic coefficient. It is clear that with increasing the shearing layer elastic coefficient of the foundation, the frequency parameters increase to some limit values and for the large values of shearing layer elastic coefficient; the frequency parameters become independent of it. • The frequency parameter decreases rapidly with the increase of the length-to-mean radius ratio and then remains almost unaltered for the long cylindrical panel (L z /R [ 5). • The interesting results show that the lowest magnitude frequency parameter is obtained by using a classicalclassical volume fraction profile. It can be concluded that a graded CNTs volume fraction in two directions has higher capabilities to reduce the natural frequency than a conventional 1-D nanocomposite. • It is found that the frequency parameter decreases by increasing the radial volume fraction index ðc r Þ. • For the all length-to-mean radius ratio (L z /R), classicalclassical volume fraction profile has the lowest frequencies followed by classical-symmetric, classical, symmetric-symmetric and symmetric profiles.
Based on the achieved results, using 2-D six-parameter power-law distribution leads to a more flexible design so that maximum or minimum value of natural frequency can be obtained in a required manner.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. a r ¼ 3 K m þ G m ð Þþk r þ l r 3 k r þ G m ð Þ b r ¼ 1 5