A Numerical Investigation on the Natural Frequencies of Fgm Sandwich Shells with Variable Thickness by the Local Generalized Differential Quadrature Method

The main aim of the present paper is to solve numerically the free vibration problem of sandwich shell structures with variable thickness and made of Functionally Graded Materials (FGMs). Several Higher-order Shear Deformation Theories (HSDTs), defined by a unified formulation, are employed in the study. The FGM structures are characterized by variable mechanical properties due to the through-the-thickness variation of the volume fraction distribution of the two constituents and the arbitrary thickness profile. A four-parameter power law expression is introduced to describe the FGMs, whereas general relations are used to define the thickness variation, which can affect both the principal coordinates of the shell reference domain. A local scheme of the Generalized Differential Quadrature (GDQ) method is employed as numerical tool. The natural frequencies are obtained varying the exponent of the volume fraction distributions using higher-order theories based on a unified formulation. The structural models considered are two-dimensional and require less degrees of freedom when compared to the corresponding three-dimensional finite element (FE) models, which require a huge number of elements to describe the same geometries accurately. A comparison of the present results with the FE solutions is carried out for the isotropic cases only, whereas the numerical results available in the literature are used to prove the validity as well as accuracy of the current approach in dealing with FGM structures characterized by a variable thickness profile.


Introduction
Shells are structural elements commonly employed by engineers, architects and designers to fulfil particular structural requirements.In fact, due to their peculiar shape and the curvature effect they are characterized by high-level stiffness, which allows to carry external loads in an efficient way.Analogously, their dynamic behavior is considerably affected by these curved geometries.Examples of shell structures can be easily seen everywhere, for example, large roofs, boat hulls, car bodies, fuselages, as well as several mechanical components, are all shell structures.Thus, it is evident that their importance is well-known in many fields, such as automotive, mechanical, civil, and aerospace engineering, architecture, aeronautical and naval industries [1][2][3].Nevertheless, assessment of their advantages due to the curvature effect are made difficult by lack of accurate analytical description of the curved middle surface.To the best of authors' knowledge, this issue can be overcome by means of the differential geometry which provide extremely simple and efficient tools to define these surfaces.As highlighted by Kraus in his book [4], shell structures can be obtained by a proper parametric description of their middle surfaces.As a consequence, doubly-curved geometries characterized by a punctual variation of the main radii of curvature can be easily taken into account.
The excellent features of shells are enhanced by the introduction of innovative materials characterized by superior mechanical properties.For example, by using composite materials it is possible to manufacture light-weight structures with higher strength, reducing the amount of material needed to build them.Consequently, the main aim of many designers and engineers is to find the best mechanical configuration to solve particular structural problems, by combining in an optimal way several layers of composite materials.Thus, laminated structures became more popular due to this aspect, as it can be noted from the huge amount of related papers available in the literature .It should be noted that laminates could be affected by some typical issues, such as delamination and stress concentrations at the interfaces, since various materials are used as constituents of the layers.In order to avoid problems associated with material mismatch at the layer interfaces, the class of Functionally Graded Materials (FGMs) was introduced.By assigning a continuous gradual variation of the mechanical properties along a specified direction, these composites do not show discontinuities in the material.As a consequence, the residual stresses and the stress concentrations that commonly affect a laminated structure can be reduced by mixing two or more constituents according to a specific graded distribution of the volume fraction.The wide use of FGMs is proven by the large number of papers published in the literature .For the sake of completeness, it should be noted that the same idea of graded materials is currently employed to define also the volume fractions of nanoparticles such as Carbon Nanotubes (CNTs) in nanocomposites [73][74][75].
The mechanical response of each structural element is obviously affected by the choice of the materials and the distribution of the various constituents.Nevertheless, the shape of these structures plays an important role as well.In fact, many designers aim to find the optimal geometry to achieve desired structural requirements.For this purpose, it is possible to modify the buckling behavior, change the bending response, and influence the vibration characteristics by varying the shape of a particular structure.Among the various approaches of structural optimization processes, the assignment of a variable thickness allows to modify the stiffness of the structure by redistributing the materials within the reference domain without increasing the structural weight.In other words, variable thickness can enhance the mechanical behavior of a structure by increasing the stiffness in those parts of the structure where the stresses are high and reducing the amount of material where it is not needed.
Several works concerning variable thickness structure are available in the literature, most of them are related to isotropic materials.A partial bet relevant literature review of this topic is presented here.Since the present paper is focused on the free vibration analysis, only the studies that were performed the same kind of structural analysis are mentioned.The work by Mizusawa [76] must be cited since it presented many investigations concerning isotropic plates characterized by different thickness profiles.In his paper, a complete analysis was performed for various geometries and different boundary conditions.In a similar manner, Shufrin and Eisenberger [77] computed the natural frequencies of plates with variable thickness by using both first-order and higher-order structural models.On the other hand, a variable kinematic Ritz approach was employed by Dozio and Carrera [78] to evaluate the vibrational behavior of quadrilateral plates with arbitrary thickness.The exact element method was employed to find the axisymmetric vibration frequencies of isotropic circular and annular plates with variable thickness by Eisenberger and Jabareen [79].On the other hand, Wu and Liu [80] solved numerically the same problem by using the generalized differential quadrature method, whereas a finite element method was proposed by Liang et al. [81] for the same purpose.The free vibration analysis of circular sandwich plates made of isotropic layers with a parabolic thickness variation was carried out by Lal and Rani [82].All the aforementioned papers are related to flat panels.The following works are focused on isotropic shell structures.The free vibrations of circular cylindrical shells with variable thickness were investigated by Duan and Koh [83] and by El-Kaabazi and Kennedy [84] by means of an analytical procedure and the Wittrick-Williams algorithm, respectively.The works by Kang and Leissa [85,86] and by Leissa and Kang [87] represented a great contribution to this topic for the results given on thick spherical and paraboloidal shells of revolution with variable thickness obtained by a three-dimensional analysis.The dynamic stiffness method was employed by Efraim and Eisenberger [88] to compute the natural frequencies of thick spherical shell panels with variable thickness for various boundary restraints, taking into account both the effects of transverse shear stresses and rotary inertia.Finally, Jiang and Redekop [89] studied the free vibrations of orthotropic toroidal shells described by the Sanders-Budiansky equations by means of a semi-analytical differential quadrature method.In this circumstance, the variable thickness profile is defined by a sinusoidal variation and the solutions are obtained for different geometric configurations.In this paper an orthotropic medium was considered.At this point, some works concerning FGM structures with variable thickness can be cited.Firstly, the exact element method was used by Efraim and Eisenberger [90] to provide some benchmark solutions to the free vibration problem of FGM annular plates.Several thickness profiles, as well as volume fraction distributions, were considered.The same kind of FGM structures resting on the Pasternak elastic foundation were investigated also by Hosseini-Hashemi et al. [91] and by Tajeddini et al. [92].In the first paper, the differential quadrature method was used to solve numerically the governing equations based on the classical plate theory, whereas the free vibration problem was solved by means of the polynomial-Ritz method in the second paper.
The main aim of the present paper is to investigate the natural frequencies of shell structures with variable thickness made of FGMs.A unified formulation is employed to derive the governing equations for several Higher-order Shear Deformation Theories (HSDTs) .As highlighted by Librescu and Reddy [115], the introduction of innovative materials in the stacking sequence of the structure could generate some effects that classical and first-order theories that are commonly neglected.For the sake of completeness, it should be noted that the same observation was illustrated also in the works [116][117][118][119][120][121].In the present paper, the field equations, as well as the corresponding boundary conditions are solved by means of the Generalized Differential Quadrature (GDQ) method developed by Shu [122].A local approach which does not take into account all the sampling points of the domain is employed to approximate the derivatives of the governing equations, following the adaptive procedure illustrated by Wang [123].Further details concerning the numerical method can be found in [124][125][126][127][128]. It should be underlined that the authors have considered the topic of variable thickness shells in the works [104,109], where laminated composite structures are investigated.In particular, a global version of the GDQ method is employed in the paper by Bacciocchi et al. [109], in which several comparisons with the results available literature allowed to validate both the structural model and the numerical technique.On the other hand, the convergence and accuracy features of the local approach of the GDQ are studied in the work by Tornabene et al. [104], where a general formulation to define arbitrary smooth thickness variations was presented.Nevertheless, in these papers the mechanical properties of the materials were assumed constant in each layer.The novelty introduced in this paper is to employ FGM plies.Thus, all the structures are characterized by a continuous gradual variation of the mechanical properties along the shell thickness.In addition, several lamination schemes and different through-the-thickness distributions of the constituents are considered to define sandwich configurations, without discontinuities at the interfaces between two contiguous layers.In the aforementioned papers, in fact, all the laminae are orthotropic and their mechanical parameters do not depend on the thickness coordinate.As far as the thickness profiles are concerned, the general approach shown in the work by Tornabene et al. [104] represents the mathematical tool employed even in this paper to describe the thickness variation by means of several functions, such as sinusoidal or power-law.
The present work is structured as follows.The mechanics of FGMs is described preliminarily.Then, some details concerning the definition of the shell geometry and the higher-order structural Appl.Sci.2017, 7, 131 4 of 39 models are presented.The main aspects of the local GDQ method are briefly illustrated, followed by some numerical results.Finally, some concluding remarks are presented.

Functionally Graded Materials
It is well-known that FGMs allows to define a continuous gradual variation of the mechanical properties of two isotropic constituents, in general one ceramic and one metal, from the bottom to the top surfaces of the layer in which they are employed.In the present paper, the rule of mixture is used to evaluate the mechanical properties of these composites.For the sake of simplicity, let us assume that the k-th layer of the structure is made of FGM and the two constituents are ceramic and metallic materials, respectively.Each mechanical quantity related to the first constituent is denoted by the subscript "C", whereas the subscript "M" is used to specify the properties of the metal.Since the two constituents are assumed isotropic and linear elastic, they are fully characterized through the Young's modulus (E C , E M ), Poisson's ratio (ν C , ν M ), and the density (ρ C , ρ M ).By hypothesis, these composites are isotropic and non-homogenous along the thickness direction, which is denoted by ζ.Thus, the rule of mixture provides the following relations where , and ρ (k) (ζ) stand for the overall mechanical properties of the k-th FGM layer.On the other hand, V M (ζ) represent the volume fraction distributions of the ceramic and the metal, respectively.These two quantities are related as follows Thus, Equations ( 1)-( 3) could be rewritten also as a function only of the volume fraction of one constituent.Several expressions are available in the literature to define the continuous and gradual variation of the mechanical properties along the thickness of the FGM layer [67].In the current paper, a four-parameter power law is introduced for this purpose.Let us assume that the volume fraction is strictly related to the materials located to the external surfaces of the k-th FGM layer.With reference to Figure 1, in which a generic stacking sequence made of FGM and isotropic layers is depicted, the k-th ply is identified by the thickness coordinates ζ k and ζ k+1 , which denote respectively the coordinates of its bottom and top surfaces in each point of the reference domain.
On the other hand, h k represents the thickness of the k-th layer in a generic point of the domain.
If the constituent at the top surface is the ceramic material, the volume fraction where a (k) , b (k) , c (k) and p (k) are the four parameters that characterize the distribution at issue.As specified in the paper by Tornabene [60], their value must be chosen accurately so that Equation ( 4) is satisfied for any value of p (k) .On the other hand, if the constituent at the top surface is the metallic material the volume fraction It should be noted that homogeneous and isotropic material can be obtained as a particular case of the FGM, by setting properly the value of the exponent p (k) .In other words, by setting p (k) = 0 or p (k) = ∞ it is possible to model an isotropic layer.In addition, it should be specified that sandwich structures can be obtained by choosing properly the volume fraction distributions of the constituents in the various plies.As a consequence, the composites are not characterized by those discontinuities at the layer interfaces that typically identify conventional laminates.At this point, a peculiar notation must be introduced to specify univocally the volume fraction distribution used to define the gradual variation of the mechanical properties.In the following, the acronyms shown below are employed for this purpose.In particular, if the Equation ( 5) is employed, the k-th FGM layer is denoted by On the other hand, the notation is used to specify the use of Equation ( 6).Finally, it is important to underline that E (k) (ζ), ν (k) (ζ), and ρ (k) (ζ) are all affected by the gradation laws specified in Equation ( 5) or (6).
by the subscript " C ", whereas the subscript " M " is used to specify the properties of the metal.Since the two constituents are assumed isotropic and linear elastic, they are fully characterized through the Young's modulus ( , ν ν ), and the density ( , C M ρ ρ ).By hypothesis, these composites are isotropic and non-homogenous along the thickness direction, which is denoted by ζ .Thus, the rule of mixture provides the following relations where V ζ represent the volume fraction distributions of the ceramic and the metal, respectively.These two quantities are related as follows Thus, Equations ( 1)-( 3) could be rewritten also as a function only of the volume fraction of one constituent.Several expressions are available in the literature to define the continuous and gradual variation of the mechanical properties along the thickness of the FGM layer [67].In the current paper, a four-parameter power law is introduced for this purpose.Let us assume that the volume fraction is strictly related to the materials located to the external surfaces of the k -th FGM layer.With reference to Figure 1, in which a generic stacking sequence made of FGM and isotropic layers is depicted, the k -th ply is identified by the thickness coordinates k

Definition of the Geometry
Let us consider a generic doubly-curved shell element with variable thickness as the one schematically depicted in Figure 2. It is completely identified by means of the position vector R(α 1 , α 2 , ζ) defined as follows where z = 2ζ/h(α 1 , α 2 ), for z ∈ [−1, 1], is a non-dimensional parameter which represents the distance between the generic point P within the three-dimensional shell body and its projection P , located on the middle surface.It should be noted that O α 1 α 2 ζ denotes the local reference coordinate system of the shell element.In particular, α 1 , α 2 are the principal curvilinear coordinates of the shell middle surface, whereas ζ specify the coordinate along the normal directions, which is identified by the outward unit normal vector n(α 1 , α 2 ).The following boundary values must be specified to define the three-dimensional domain of the shell.
From the definition shown in Equation ( 9), it is clear that the thickness of the shell h(α 1 , α 2 ) is arbitrary and can vary point by point.
From the definition shown in Equation ( 9), it is clear that the thickness of the shell ( ) , h α α is arbitrary and can vary point by point.It should be recalled that its overall value is given by the sum of the thickness of each ply if the structure is made of several layers (Figure 1).Thus, one gets ( ) ( ) where l stands for the number of layers.In the present paper, the thickness variation assumes the following aspect in which 0 h is the thickness reference value, whereas ( ) ( ) are arbitrary functions that describe the thickness profile along the principal curvilinear coordinates 1 2 , α α .It should be pointed out that there is no restriction on the choice of these functions, as illustrated in details in the work by Tornabene et al. [105].Nevertheless, they must represent a smooth variation and their effect has to involve each layer.Thus, the parameters 0,k h , for 1, 2,..., k l = , are introduced to specify the thickness reference value of each lamina.For completeness purposes, the reader can refer to the work by Tornabene et al. [105] for a general treatise on the variable thickness.In connection with Equation ( 9), , α α r is the vector that allows to identify each point of the shell reference surface, which corresponds to the middle surface of the structure.In general, its expression can be written as follows It should be recalled that its overall value is given by the sum of the thickness of each ply if the structure is made of several layers (Figure 1).Thus, one gets where l stands for the number of layers.In the present paper, the thickness variation assumes the following aspect in which h 0 is the thickness reference value, whereas f (α 1 ), g(α 2 ) are arbitrary functions that describe the thickness profile along the principal curvilinear coordinates α 1 , α 2 .It should be pointed out that there is no restriction on the choice of these functions, as illustrated in details in the work by Tornabene et al. [104].Nevertheless, they must represent a smooth variation and their effect has to involve each layer.Thus, the parameters h 0,k , for k = 1, 2, ..., l, are introduced to specify the thickness reference value of each lamina.For completeness purposes, the reader can refer to the work by Tornabene et al. [104] for a general treatise on the variable thickness.In connection with Equation (9), r(α 1 , α 2 ) is the vector that allows to identify each point of the shell reference surface, which corresponds to the middle surface of the structure.In general, its expression can be written as follows (15) in which f i (α 1 , α 2 ), for i = 1, 2, 3, are generic functions that depend on the shell middle surface, whereas e i , for i = 1, 2, 3, are the unit vectors which denote the main axes of the global reference system O x 1 x 2 x 3 shown in Figure 2. Equation ( 15) is fundamental to evaluate each geometric quantity related to the shell under consideration.For this purpose, its derivatives with respect to the principal curvilinear coordinates must be also computed.As far as the first-order derivatives are concerned, one gets whereas the second-order derivatives assume the following aspect At this point, an analytic expression can be found for the unit normal vector n(α 1 , α 2 ) where "∧" denotes the vector product.Analogously, the definition of the Lamè parameters where "•" denotes the scalar product.Finally, the middle surface of a doubly-curved shell is characterized by two principal radii of curvature R 1 (α 1 , α 2 ), R 2 (α 1 , α 2 ) that vary point by point.
The following expressions can be used to define them Equations ( 23) and ( 24) are valid since α 1 , α 2 are principal and orthogonal by hypothesis.Once definition of R 1 (α 1 , α 2 ), R 2 (α 1 , α 2 ) is specified, the following quantities can be also computed It should be specified that the parameters are required to take into account the three-dimensional size of a generic shell structure.Differential geometry is the mathematical tool that allows to compute all the geometric quantities defined in the present section.

Higher-Order Equivalent Single Layer Approach
A unified formulation is employed in the present paper to define several HSDTs in a compact manner.In particular, an Equivalent Single Layer (ESL) approach is taken into account.The three-dimensional displacements for a generic composite shell made of several layers can be written as a function of a generic order of kinematic expansion τ.
Where the generalized displacements u (τ) represent the degrees of freedom of the problem and are defined on the shell middle surface.They can be conveniently collected in the algebraic vector u (τ) (α 1 , α 2 , t) for each order τ of kinematic expansion.
On the other hand, F τ (ζ) denotes the shear function (or thickness functions) linked to the τ-th order of kinematic expansion and can assumes different meanings.In the current work, the power function ζ τ , for τ = 0, 1, 2, ..., N, is chosen for this purpose.The (N + 1)-st order of kinematic expansion denotes the Murakami's function that is employed to model the zig-zag effect [94][95][96].For the laminated shell element shown in Figure 1, the Murakami's function is defined as follows: in which the index k is introduced to specify the k-th ply.Thus, the (N + 1)-st shear function is assumed equal to Equation ( 29) if the zig-zag effect is included in the structural theory (in other words, one gets F N+1 = Z).It should be noted that each HSDT is identified by the maximum order of kinematic expansion N. If the power function ζ τ is employed to define all the thickness functions of the model, the following HSDT can be introduced.
When the zig-zag effect is contemplated, the same order theories are specified as follows: With reference to the Equations ( 30) and (31), it should be noted that the letter "E" specifies an ESL model, "D" means that the fundamental equations will be written in terms of generalized displacements, and "Z" stands for the Murakami's function.At this point it is important to underline that the displacement field of the well-known first-order shear deformation theory can be obtained from Equation ( 27) by choosing properly the shear functions and the degrees of freedom.In the following, this theory will be indicated by the notation FSDT (First-order Shear Deformation Theory).Analogously, the term FSDTZ is introduced to identify the FSDT embedded with the Murakami's function.
At this point, the generalized strain components can be defined.For the sake of conciseness, they are collected in the corresponding algebraic vector ε (τ) (α 1 , α 2 , t) for each order τ of kinematic expansion.The generalized strain components assume the compact matrix form shown below for each order of kinematic expansion τ = 0, 1, 2, ..., N, N + 1 where D Ω is the differential operator defined as The internal stress resultants for each order of kinematic expansion τ = 0, 1, 2, ..., N, N + 1 are collected in the algebraic vector S (τ) (α 1 , α 2 , t) which assume the following aspect.
Due to the duality between stresses and strains, the generalized stress resultants can be directly related to the corresponding generalized displacements.In compact matrix form, one gets for τ = 0, 1, 2, ..., N, N + 1, in which A (τη) is the matrix of the elastic constants.Equation ( 36) is of extremely importance since it gives also the definition of the stress resultants involved in those boundary conditions that affect the stresses, such as in the free edges.In general, the matrix A (τη)  assumes the following aspect for a shell made of l layers.
11( 20 where for τ, η = 0, 1, 2, ..., N, N + 1, n, m = 1, 2, 3, 4, 5, 6 and p, q = 0, 1, 2. For the sake of conciseness, the symbols τ, η denote the shear function order, whereas τ, η stand for the corresponding derivatives with respect to ζ.On the other hand, B nm specifies the elastic constant of the material.It assumes the following definition for n, m = 1, 2, 3, 6 whereas for n, m = 4, 5 one gets where κ = 1/χ is the shear correction factor.If the structural model requires this correction, the value of χ is set equal to χ = 1.2, otherwise the unitary value is employed (χ = 1).Equations ( 39) and ( 40) are general and allows to introduce also the plane stress hypothesis.In fact, E nm denotes the elastic constant of the material of the k-th ply.If the plane stress hypothesis is required, the reduced elastic coefficients must be used (E nm ).On the other hand, the non-reduced coefficients are employed (E nm ).These aspects will be specified properly by adding specific subscripts and superscripts to the notation used to denote the structural model.For instance, the Reissner-Mindlin theory needs both the shear correction factor and the plane stress-reduced elastic coefficients.As a consequence, it will be indicated as FSDT χ=1.2 RS , where RS stands for "reduced stiffness".Finally, it should be specified that the quantities E (k) nm depend on the third coordinate ζ if the layer is made of FGM.In other words, one gets A set of three motion equations for each order of kinematic expansion τ = 0, 1, 2, ..., N, N + 1 is obtained.In compact matrix form, they assume the following aspect The equilibrium operator D * Ω is defined as follows: On the other hand, the mass matrix M (τη) is defined below: for τ, η = 0, 1, 2, ..., N, N + 1.If ρ (k) stands for the mass density of the k-th layer, the inertia masses I (τη) in Equation ( 43) are computed for each order of kinematic expansion τ, η = 0, 1, 2, ..., N, N + 1 as follows With reference to the motion Equation ( 41), it should be specified that ..
is the algebraic vector that collects the second-order derivatives of the generalized displacements Equation ( 28) with respect to the time variable t.The system of governing Equation ( 41) can be expressed as a function of the degrees of freedom by inserting Equation ( 36) into (41).The fundamental equations are obtained for each order of kinematic expansion τ, η = 0, 1, 2, ..., N, N + 1 where L (τη) is the fundamental operator defined as follows The terms collected into Equation ( 46) can be computed by applying the following definition for τ, η = 0, 1, 2, ..., N, N + 1.It should be pointed out that Equation ( 45) define a set of 3 × (N + 2) motion equations.Since it represents a set of partial differential equations, the proper boundary conditions must be enforced to find the solution.At this point, Figure 2 can be taken as a reference to identify the coordinates of each edge along which the boundary conditions are applied.Let us consider fully clamped edges (specified by "C") or free edges (denoted by "F").The boundary conditions assume the following aspect . This is the case of the South (S) and North (N) edges.The same boundary conditions can be written as follows: . This is the case of the West (W) and East (E) edges.In the next sections related to the numerical applications, the boundary conditions applied along each edge will be specified following the order WSEN.Finally, it should be noted that another set of compatibility conditions must be enforced if the considered structure has a common edge, as in the case of complete shells of revolution or toroidal panels.If the closing edge involves the coordinates α 1 = α 0 1 and α 1 = α 1 1 , the following conditions must be added in terms of stress resultants and displacements On the other hand, if the closing edge is identified by α 2 = α 0 2 and α 2 = α 1 2 , the compatibility conditions in terms of stress resultants become whereas the kinematic conditions assume the aspect written below Finally, it should be specified that the presented theoretical approach is employed to analyze the mechanical behavior of thick and moderately thick shells that are defined by the following limitations where R min and L min stand for the minimum value of the radii of curvature and the minimum size, respectively.

Numerical Aspects
In this section, the numerical techniques employed in the present work are briefly illustrated.In particular, the GDQ method is used to approximate the derivatives in order to solve the strong form of the governing equations and to obtain the geometric parameters required to define the shell geometry.On the other hand, the integrals needed to evaluate the elastic constants of the stiffness matrix are computed by means of the Generalized Integral Quadrature (GIQ).For further details concerning these numerical approaches, the reader can refer to the review paper by Tornabene et al. [124].

Local Generalized Differential Quadrature Method
The Local Generalized Differential Quadrature (LGDQ) method is a numerical tool that allows to approximate the derivatives of a generic function without considering all the points of the reference domain.In other words, local interpolating basis functions are taken into account instead of their corresponding global ones, if compared to the well-known GDQ method.The main features of the LGDQ are explained in this paragraph considering a one-dimensional domain for conciseness purposes.The n-th derivative of a sufficiently smooth function f (x) evaluated at a generic point x i can be computed as a weighted linear sum of the function values at some selected grid points, considering only a limited part of the whole domain, which is globally made up of I N discrete points.From the mathematical point of view, one gets for i = 1, 2, ..., I N , where the symbol ς (n) ij denotes the weighting coefficients for derivative in hand.It is clear that the evaluation of these weighting coefficients and the choice of a proper grid distribution represent the key points of the present technique.The limit values of the sum N 1 , N 2 depend on the position of the point x i where the derivative is computed.The following relations must be introduced to define their value for i = 1, 2, ..., I N .The symbol I PN stands for the number of grid points that must be taken into account on the left and on the right of the point x i itself.It should be noted that the value of I PN has to be set a priori.Due to the definitions shown in Equation ( 56), the values of N 1 , N 2 is variable and assumes different meanings according to the position of the point x i within the domain.Figure 3 is introduced now to explain more clearly the significance of Equation ( 56).Let us assume that x i is placed in the center of the domain so that on both its left and its right the number of points is greater than I PN .As it can be noted from Figure 3a, the derivative shown in Equation ( 55) is computed considering (2I PN + 1) points, since N 1 = i − I PN and N 2 = i + I PN .Let us consider the circumstance in which x i is located next to the limits of the domain (or corresponds also to the boundary points).
In these case, the value of I PN can be greater than the number of points that actually appear on the left (Figure 3b) or on the right (Figure 3c) of As far as the evaluation of the weighting coefficients ς (n) ij is concerned, the formulas shown in the works [104,123] can be used.For the first-order derivatives, one gets ς On the other hand, the following recursive relations must be employed for the n-th order derivatives The formulation just presented can be easily extended to the case of two-dimensional domain following the same considerations shown in the works [108,109].It should be noted that in this case a discrete grid distribution must be applied along both the principal coordinate directions.Let us assume that the domain is bounded by the limitations shown in Equations ( 10)- (12).The discrete points within the shell reference domain can be placed as follows for i = 1, 2, ..., I N and j = 1, 2, ..., I M , where I N , I M denote the total number of points along α 1 , α 2 , respectively.Due to the general features of the present approach, several grid distributions can be considered to define the value of r i , r j .Due to the accuracy and convergence characteristics shown in the work [104], the points are placed within the domain by evaluating the roots of the Chebyshev polynomials of the second kind.In other words, r i , r j assume the following aspect for i = 1, 2, ..., I N , j = 1, 2, ..., I M , and r ∈ [−1, 1].It should be noted that a local interval of points must be defined also along the second principal direction of the two-dimensional domain.For this purpose, the parameter I PM is introduced with the same meaning of I PN .For the sake of simplicity, in the present work the value of I PM corresponds to I PN .In general, the use of the LGDQ should be more advisable to deal with those physical problems in which the properties of the domain change point by point.With this aim, a substantial number of discrete points could be required to define accurately these features, as in the case of shells with variable thickness presented in this work.As a consequence, the cost could be high from the computational point of view if a global approach is employed.Nevertheless, the LGDQ represents an efficient tool to reduce the computational effort since the matrix that collects the coefficients is banded, as explained in the literature [102,104].

Generalized Integral Quadrature Method
The GIQ method is a numerical tool that allows to evaluate integrals.Let us consider a smooth function f (x) defined in a closed domain [a, b].In addition, let us assume that I T points are placed to discretize the domain itself, so that one gets a = x 1 , x 2 , . . ., x I T −1 , x I T = b.As a consequence, I T − 1 subdomains are defined.The continuous function f (x) can be approximated by the Lagrange polynomials of order I T − 1, since I T points are employed to discretize the domain.The integral of the approximated function in the closed interval x i , x j is defined as follows In other words, the numerical integral of f (x) is a weighted linear sum of the values that the function f (x) assumes in each point of the reference domain.It should be noted that this procedure involves the values that the function assumes in the points placed outside the integration domain.As in the GDQ method, w ij k are the weighting coefficients required for the integration, which can be computed as shown below.Let us assume that the function f (x) is approximated by a polynomial expansion as in which a i , for i = 1, 2, ..., I T , denote the arbitrary constants of the polynomial.An auxiliary function F(x) must be introduced so that Due to Equation ( 68), the function F(x) can be defined as in which I f (x, c) is the integrand function that assumes the following aspect Appl.Sci.2017, 7, 131 16 of 39 On the other hand, c denotes an arbitrary constant.Due to the Weierstrass theorem and the properties of a linear vector space, the expression below represents an approximation of the function I f (x, c) where l j (x) specifies the Lagrange polynomials.The first derivative of I f (x, c) with respect to x can be easily evaluated at the point x i dI f (x, c) dx where l (1) j (x i ) stands for the first-order derivative with respect to x of the Lagrange polynomials computed at x i .It should be noted that l (1) j (x i ) stands for the first-order derivative with respect to x of the Lagrange polynomials computed at x i and coincides with the weighting coefficients of the first order derivative ς (1) ij given by the GDQ method [124].By applying the differential quadrature law, Equation ( 73) becomes in which the terms l k x j = δ kj represent the Kronecker delta.The following results is achieved by comparing Equations ( 73) and ( 74) f (x, c) where the weighting coefficients assume the form for i = j, and ii + for i = j.The arbitrary constant c must be set equal to c = x I T + 10 −10 due to accuracy and stability reasons [124].At this point, Equation ( 74) can be conveniently written in matrix form in which ς (1) is the matrix that collects the weighting coefficients, and the following vectors are introduced I (1) f (x 2 , c), . . . ,I Recalling Equation ( 70), vector I f can be rewritten as Thus, the integrals shown in Equation ( 81) can be evaluated as follows where W = ς (1) −1 is the matrix of the weighting coefficients for the integrals that can be computed by inverting the corresponding matrix of the weighting coefficients for the first-order derivatives.One gets for i = 1, 2, ..., I T .Finally, the integral in Equation ( 67) is obtained where the weighting coefficients for the GIQ method are given by It should be clear that they can be evaluated directly from the expressions of the weighting coefficients used to approximate the first order derivatives [124].A proper grid distribution must be introduced to discretize the reference domain.Since in the present work the integrals must be evaluated along the shell thickness to compute the elastic coefficients of the stiffness matrix, the points along the coordinate ζ can be located as for m = 1, 2, ..., I T .Equation ( 86) represents the well-known Chebyshev-Gauss-Lobatto grid distribution.It is important to notice that each node defined by Equation ( 86) depends on the principal coordinates α 1 , α 2 since the thickness of the considered structures is variable.The value of I T is fixed equal to 51 for each numerical analysis presented in the following sections.

Free Vibration Analysis
The governing system of equations can be written in the following form for τ = 0, 1, 2, ..., N, N + 1, if the solution is conveniently assumed equal to where U (τ) is the vector of the mode shapes defined as follows for each order of kinematic expansion τ.It should be noted that ω denotes the corresponding circular frequencies, from which it is possible to evaluate the natural frequencies as f = ω/2π.
At this point, Equation ( 87) can be solved by means of the LGDQ method.In other words, the derivatives that appear in the governing equation can be approximated through Equation ( 55) and written in their discrete form.Thus, the governing system of equations assume the following aspect in which K and M represent the stiffness and the mass matrices, respectively.On the other hand, δ is the vector that collects the vibration spatial amplitude values.It is well-known that Equation ( 90) represents a generalized linear eigenvalue problem.The eigenvalue ω k are evidently the circular frequencies of the considered structure.Once their value is computed, the corresponding eigenvectors δ k that represent the mode shapes can be evaluated, too.It should be specified that the maximum number of eigenvalues, as well as eigenvectors, depends on the number of points I N , I M chosen to discretize the reference domain along the two principal directions α 1 , α 2 , respectively.It is point out that the size of system of discrete Equation ( 90) where N denotes the maximum order of kinematic expansion.It should be noted that the size of the problem can be reduced by applying the well-known kinematic condensation of non-domain degrees of freedom in order to separate the quantities associated with the central points of the domain from the ones related to the boundaries.

Numerical Results
The four structures with variable thickness made of FGMs depicted in Figure 4 are considered in the present section to evaluate their natural frequencies by the LGDQ method.As previously illustrated, the reference domain (or middle surface) is described by the corresponding position vector r(α 1 , α 2 ), whereas a proper expression is introduced to define the smooth thickness variation of the shells.With reference to Figure 4, it should be specified that the shell cross-sections affected by the thickness variation are represented, too.To this end, the local coordinate reference system is also depicted to understand more clearly the thickness variation (the principal axes are specified by the colors magenta, cyan and green, respectively).Different stacking sequences and volume fraction distributions are taken into account to characterize the mechanical features of the composites (Figure 5).The two isotropic constituents and their corresponding properties are listed in Table 1 for conciseness purposes.Firstly, the proposed method is validated by the comparison with the results available in the literature.Then, several parametric investigations are performed to show the effect of the variation of the volume fraction distributions on the natural frequencies.In these circumstances, the results are compared with the ones obtained by a three-dimensional finite element model (3D-FEM) for the limit cases.In other words, the commercial code Abaqus is employed only to evaluate the natural frequencies of composite structures made of isotropic layers.For this purpose, the three-dimensional model is obtained by using 20-node hexahedrons elements (C3D20).On the other hand, all the results shown in the following tables are obtained by the LGDQ method implemented in a MATLAB code [128].In the present solution, the degrees of freedom are equal to 3(N + 2) × (I N × I M ), or 3(N + 1) × (I N × I M ) in case the Murakami's function is neglected in the HSDT.
model is obtained by using 20-node hexahedrons elements (C3D20).On the other hand, all the results shown in the following tables are obtained by the LGDQ method implemented in a MATLAB code [128].In the present solution, the degrees of freedom are equal to      (a) ( )

Comparison with the Literature
The FGM annular plate defined in the work by Efraim and Eisenberger [90] is considered in the current section to validate the proposed approach.The plate middle surface is described by the following position vector where x, ϑ are the principal coordinate of the surface, assuming x ∈ [0, L] and ϑ ∈ [0, 2π].The inner radius is denoted by R i , whereas the outer one can be computed as R out = R i + L. A linear variation is applied along the radial direction to define the thickness profile in which H is a constant parameter.The following relations must be taken into account to describe completely the plate geometry For the sake of completeness, both the structure and its radial cross-section are depicted in Figure 4a.The upper surface of the plate is made of Silicon nitride and a linear variation of the mechanical properties is assumed to describe the functionally graded composite layer.On the other hand, the lower surface is completely made of Stainless steel.The first 20 natural frequencies (apart from free body motions) are shown in Table 2 in their dimensionless form As far as the boundary conditions are concerned, the structure is taken completely free (FF).The LGDQ solutions are obtained by setting I N = I M = 41 and I PN = I PM = 14.Three different structural models are considered: FSDT χ=1.2 RS , ED2 χ=1.2 , and ED3.It should be noted that both the first and second order theories require the shear correction factor χ = 1.2, whereas the hypothesis of plane stress is take into account only in the FSDT.From the results shown in Table 2, it can be noted that the numerical values are in good agreement with the reference solutions for each model.

Elliptic Cylinder
The next application aims to evaluate the natural frequencies of a FC elliptic cylinder made of two FGM layers of equal thickness.The geometry under consideration is depicted in Figure 4b.The stacking sequence, as well as the volume fraction distribution of the constituents, are shown in Figure 5a for different values of the exponent p (1) = p (2) = p.The middle surface of the shell is given by the following position vector where a = 1 m and b = 1.5 m are the semi diameters of the ellipse.The coordinates α 1 , y are bounded respectively by the limitations α 1 ∈ α 0 1 , α 1 1 and y ∈ [0, L], with α 0 1 = −π/2, α 1 1 = 3π/2, and L = 4 m.A sine-wave variation is applied along the first coordinate α 1 to define the thickness profile.Mathematically speaking, the thickness is described by the following expression with h 0 = 0.2 m.The first 10 natural frequencies for the considered structure are shown in Table 3 for different values of the parameter p ∈ [0, ∞], varying the structural model.As in the previous example, the FSDT χ=1.2 RS , ED2 χ=1.2 , and ED3 are employed.The same considerations concerning the shear correction factor and the reduced elastic stiffnesses are still valid.
It should be noted that the boundary values of p denote a structure completely made of Silicon nitride (p = 0) and Stainless steel (p = ∞), respectively.In this cases, the comparison with the 3D-FEM solutions is performed with excellent results for all the structural theories.It should be noted that the finite element model is made of 69680 brick elements (924690 degrees of freedom).On the other hand, the LGDQ solutions is obtained by setting I N = 51, I M = 31, and I PN = I PM = 16.These values are chosen to describe accurately both the elliptical shape and the thickness variation.For completeness purposes, the variation of the first natural frequency is depicted graphically in Figure 6 as a function of the exponent p.It is evident that for lower values of p the frequencies are highly affected by variation of the volume fraction distribution of the constituents.On the other hand, for p > 4 ÷ 6 the variation in terms of natural frequencies is considerably restricted.The same behavior can be observed for each structural model.In the same plots, the black dashed lines represent the extreme cases of the isotropic homogeneous shell (p = 0 and p = ∞) and the corresponding frequencies are evaluated by the finite element software.For the sake of conciseness, Figure 6 shows only the variation of the first natural frequency, but it should be specified that higher modes exhibit the same trend.It is evident that the variation of the natural frequencies is bounded by these two lines.Finally, the first three mode shapes for a FC elliptic cylinder with variable thickness are shown in Figure 7 assuming p = 1, which corresponds to a linear distribution of the volume fraction.From these images, it can be noted that the boundary conditions are well-enforced.In this paper, the colormap used for the mode shapes stands for the intensity of the modal displacements.A zero displacement is specified by a dark blue color, whereas a red tone represents bigger movements.Table 3.First 10 natural frequencies f (Hz) for a FC elliptic cylinder with variable thickness, varying the value of the exponents p (1) = p (2) = p of the volume fraction distributions.The LGDQ solution is obtained for several HSDTs by setting I N = 51, I M = 31 and I PN = I PM = 16.As far as the boundary conditions are concerned, the structure is taken completely free (FF).The LGDQ solutions are obtained by setting 41  2, it can be noted that the numerical values are in good agreement with the reference solutions for each model.

Elliptic Cylinder
The next application aims to evaluate the natural frequencies of a FC elliptic cylinder made of two FGM layers of equal thickness.The geometry under consideration is depicted in Figure 4b.The stacking sequence, as well as the volume fraction distribution of the constituents, are shown in Figure 5a for different values of the exponent (1)   (2) The middle surface of the shell is given by the following position vector ( ) where 1m a = and 1.5 m b = are the semi diameters of the ellipse.The coordinates 1 , y α are bounded respectively by the limitations = , and 4m L = . A sine-wave variation is applied along the first coordinate 1 α to define the thickness profile.Mathematically speaking, the thickness is described by the following expression ( )   lines.Finally, the first three mode shapes for a FC elliptic cylinder with variable thickness are shown in Figure 7 assuming 1 p = , which corresponds to a linear distribution of the volume fraction.From these images, it can be noted that the boundary conditions are well-enforced.In this paper, the colormap used for the mode shapes stands for the intensity of the modal displacements.A zero displacement is specified by a dark blue color, whereas a red tone represents bigger movements.

Doubly-curved Shell of Revolution with Hyperbolic Meridian
The FC hyperbolic meridian shell, depicted in Figure 4c, is a doubly-curved shell of revolution whose reference surface is defined a branch of a hyperbola.The surface at issue is described by the following position vector The coordinates , ϕ ϑ are defined in the intervals

Doubly-Curved Shell of Revolution with Hyperbolic Meridian
The FC hyperbolic meridian shell, depicted in Figure 4c, is a doubly-curved shell of revolution whose reference surface is defined a branch of a hyperbola.The surface at issue is described by the following position vector The coordinates ϕ, ϑ are defined in the intervals ϕ ∈ [ϕ 0 , ϕ 1 ] and ϑ ∈ [ϑ 0 , ϑ 1 ], where ϕ 0 = 1.741190, ϕ 1 = 1.212026, ϑ 0 = 0, ϑ 1 = 2π.On the other hand, the parallel radius R 0 (ϕ) assumes the aspect shown below where k = √ 1 + a 2 /b 2 .It should be specified that a is the parallel radius of the throat section of the structure, and b = aD/ √ d 2 − a 2 is the characteristic dimension of the shell, in which a = 1 m, d = 2 m, D = 4 m.The meaning of these parameters, as well as the ones needed to sketch the following structure, can be found in Figure 8.For further details concerning the present geometry the reader can refer to the book by Tornabene et al. [2].With reference to Figure 8a, the characteristic dimension of the hyperbola can be defined alternatively as b = aC/ √ c 2 − a 2 , where c is the parallel radius at the top section of the structure and C represents its coordinate that in the present case is equal to C = 1 m.The value of c can be computed by means of Equation (99).The reference value of thickness is h 0 = 0.2 m and it is measured along the top edge.A cubic variation is applied along the first coordinate ϕ and it can be defined as follows: sequence depicted in Figure 5b, where the volume fraction distributions of the two constituents are also specified for several values of the   It should be noted that the same thickness variation affects each ply, as it can be observed in the cross-section of Figure 4c.The structure is made of five layers with different thickness reference values.The two external laminae are isotropic and fully ceramic (h 0,1 = h 0,5 = 0.02 m), whereas the inner isotropic core is made of metal and its thickness is equal to h 0,3 = 0.06 m.The remaining two plies have equal thickness h 0,2 = h 0,4 = 0.05 m and are graded as it can be noted from the stacking sequence depicted in Figure 5b, where the volume fraction distributions of the two constituents are also specified for several values of the exponent p (2) = p (4) = p ∈ [0, ∞].In a similar way to the previous example, the boundary values of p define respectively a layer fully made of Silicon nitride (p = 0) and Stainless steel (p = ∞).The results of this parametric study are given in Table 4 for the FSDT χ=1. 2  RS , ED2 χ=1.2 , and ED3 models.
Due to the presence of various layers with different mechanical properties, the influence of the Murakami's function is also investigated.Thus, the aforementioned structural theories are also embedded with this function.The results related to the FSDTZ χ=1.2 RS , EDZ2 χ=1.2 , and EDZ3 models are shown in Table 5.Even in this circumstance, the shear correction factor is used for the first and second order models, whereas the plane stress hypothesis is considered only in the Reissner-Mindlin theory.It is easy to notice that the Murakami's function does not significantly affect the results.Both in Tables 4 and 5, the 3D-FEM solutions are also included for comparison purposes.It should be noted that the finite element model is made of 52,734 brick elements (685,824 degrees of freedom).Excellent agreement with the reference solutions is achieved in each circumstance.The LGDQ results are obtained by setting I N = I M = 31 and I PN = I PM = 16, and the natural frequencies are written for different values of p.It should be noted that Tables 4 and 5 show the first 18 natural frequencies.In fact, some couples of natural frequencies are equal due to the symmetry of the structure, in terms of both geometry and materials.For the sake of clarity, the variation of the first six natural frequencies as a function of the exponent p are depicted in graphical form in Figure 9 for all the considered structural models, with and without the Murakami's function.As in the previous case, a remarkable difference in the structural response can be obtained for lower values of the parameter p. Finally, the first six mode shapes for the doubly-curved shell of revolution with hyperbolic meridian with variable thickness are shown in Figure 10 assuming a linear variation for the volume fraction distributions of the two constituents in the FGM layers (p = 1).
1st-2nd frequencies 3rd-4th frequencies 5th-6th frequencies   Table 6.First 10 natural frequencies f (Hz) for a CFCF doubly-curved panel of translation with variable thickness, varying the value of the exponents (2)  p p = of the volume fraction distribution.
The LGDQ solution is obtained for several HSDTs by setting 31
It is well-known that a parabolic curve is completely defined by the definition of three points with abscissas a, c, d and the relative distance b between the points of coordinates a and d.In particular, a and c represent the abscissas of the extreme points of the parabolic arch as it can be deducted from Figure 8b.The current geometry can be obtained by setting a = 3 m, b = 1 m, c = −3 m, and d = 0 m.For extra details concerning the present shape, the book by Tornabene et al. [2] can be taken into account as a reference.
As it can be noted from Figure 5c, the structure is made of three layers.In particular, the lower ply is isotropic and fully metal (Stainless steel), and the upper one is completely ceramic (Silicon nitride).These two layers have the same reference thickness h 0,1 = h 0,3 = 0.05 m.On the other hand, the central core is graded between these two limit cases and its reference thickness is equal to h 0,2 = 0.2 m.Thus, the overall shell thickness that can be assumed as a reference is h 0 = 0.3 m.Varying the exponent p (2) = p ∈ [0, ∞] of the volume fraction distribution, several mechanical configurations can be obtained.
As far as the thickness profile is concerned, the thickness variation affects both the two principal coordinates α 1 , α 2 , as it can be easily observed in the cross-sections shown in Figure 4d.In particular, a quadratic variation is applied along α 1 , whereas a sine-wave variation is applied on the second direction.Mathematically speaking, the thickness of the shell at issue is defined as follows Even in this circumstance, the three structural models previously employed are taken with and without the Murakami's function.The natural frequencies related to the FSDT χ=1.2 RS , ED2 χ=1.2 , and ED3 are shown in Table 6, whereas the ones obtained by using the FSDTZ χ=1.2 RS , EDZ2 χ=1.2 , and EDZ3are written in Table 7, for several values of the parameter p.For the sake of completeness, the 3D-FEM solutions, which are really close to the LGDQ ones, are also included in these tables.

Conclusions
The GDQ method has been applied according to a local scheme to compute natural frequencies of shell structures with variable thickness made of FGMs.The thickness profiles have been defined through a number of smooth functions (linear, power-law, sinusoidal, and their combinations), whereas a four-parameter power law function has been adopted to describe the through-the-thickness volume fraction distribution of the two constituents.As a consequence, the considered shells have been characterized by a variation of the mechanical properties both along the thickness and within the reference domain.It should be noted that the higher-order based structural models employed in this paper are two-dimensional and have allowed to evaluate accurately the natural frequencies of variable thickness structures.The same geometries have been analyzed, for verification purposes, also by means of a FEM commercial code for specific mechanical configurations.In this circumstance, a three-dimensional model made of many brick elements has been required to model the variable thickness profiles, considerably increasing the computational effort.In addition, the cost of computation has been further reduced by the use of the LGDQ method, which does not take in account all the discrete points within the domain.Consequently, the thickness variation has been described precisely by setting more grid points, without affecting the computational resources.Finally, a parametric investigation has been carried out to analyze the effect of the graded mechanical properties along the shell thickness varying the exponent of the volume fraction distributions, for several structural theories defined by different orders of kinematic expansion.The validity of the current approach has been proven for FGM structures with variable thickness through the comparison with the results available in the literature.To sum up, it has been shown that the LGDQ approach can be considered as an efficient method to solve the free vibration problem of shell structures with variable thickness with a reduced computational cost if compared to the corresponding FEM models, which required a three-dimensional description characterized by a higher number of degrees of freedom to obtain comparable results.

k
ν ζ , and ( ) ( ) k ρ ζ stand for the overall mechanical properties of the k -th FGM layer.On the other hand,

ζ and 1
k ζ + , which denote respectively the coordinates of its bottom and top surfaces in each point of the reference domain.

Figure 1 .
Figure 1.Stacking sequence of a laminated composite structure with variable thickness made of several Functionally Graded Materials (FGM) or isotropic layers.

Figure 1 .
Figure 1.Stacking sequence of a laminated composite structure with variable thickness made of several Functionally Graded Materials (FGM) or isotropic layers.

Figure 2 .
Figure 2. Generic doubly-curved shell element with variable thickness: edge identification, global and local reference coordinate system representation.

Figure 2 .
Figure 2. Generic doubly-curved shell element with variable thickness: edge identification, global and local reference coordinate system representation.

39 Figure 3 .Figure 3 .
Figure 3. One-dimensional scheme for the Local Generalized Differential Quadrature (LGDQ) method considering different positions of the i -th point where the derivative is computed: (a) point placed in the center of the domain; (b) point located next to the limits of the domain; (c) point placed on the boundary.

Quadratic thickness variation along 1 Sine-wave thickness variation along 2 Figure 4 .
Figure 4. Definition of the shell structures and graphical representation of their cross-sections with variable thickness: (a) annular plate [90]; (b) elliptic cylinder; (c) doubly-curved shell of revolution with hyperbolic meridian; (d) doubly-curved panel of translation (a parabola slides over another parabola).

Figure 4 .
Figure 4. Definition of the shell structures and graphical representation of their cross-sections with variable thickness: (a) annular plate [90]; (b) elliptic cylinder; (c) doubly-curved shell of revolution with hyperbolic meridian; (d) doubly-curved panel of translation (a parabola slides over another parabola).

Figure 5 .
Figure 5. Volume fraction distributions and lamination schemes for the shells with variable thickness under consideration: (a) Elliptic cylinder; (b) Doubly-curved shell of revolution with hyperbolic meridian; (c) Doubly-curved panel of translation.

Figure 5 .
Figure 5. Volume fraction distributions and lamination schemes for the shells with variable thickness under consideration: (a) Elliptic cylinder; (b) Doubly-curved shell of revolution with hyperbolic meridian; (c) Doubly-curved panel of translation.

2 ED2
χ = , and ED3 .It should be noted that both the first and second order theories require the shear correction factor 1.2 χ = , whereas the hypothesis of plane stress is take into account only in the FSDT.From the results shown in Table 97) with 0 0.2 m h = .The first 10 natural frequencies for the considered structure are shown in Table 3 for different values of the parameter [ ] 0, p ∈ ∞ , varying the structural model.As in the previous example, and ED3 are employed.The same considerations concerning the shear correction factor and the reduced elastic stiffnesses are still valid.

Figure 6 .
Figure 6.Variation of the first natural frequency for a free-clamped FC elliptic cylinder with variable thickness as a function of the exponent p .Similar graphs can be obtained for the second and third frequencies.

Figure 6 .
Figure 6.Variation of the first natural frequency for a free-clamped FC elliptic cylinder with variable thickness as a function of the exponent p. Similar graphs can be obtained for the second and third frequencies.

Figure 7 .
Figure 7. First three mode shapes for a FC elliptic cylinder with variable thickness assuming 1 p = as exponent.

Figure 7 .
Figure 7. First three mode shapes for a FC elliptic cylinder with variable thickness assuming p = 1 as exponent.

2 ED2
. In a similar way to the previous example, the boundary values of p define respectively a layer fully made of Silicon nitride ( 0 p = ) and Stainless steel ( p = ∞ ).The results of this parametric study are given in Table4for the χ = , and ED3 models.

Figure 8 .
Figure 8. Meaning of the geometric parameters required to sketch two different curves: (a) hyperbola; (b) parabola.

Figure 8 .
Figure 8. Meaning of the geometric parameters required to sketch two different curves: (a) hyperbola; (b) parabola.

Figure 9 .
Figure 9. Variation of the first six natural frequencies for a CF doubly-curved shell of revolution with hyperbolic meridian with variable thickness as a function of the exponent p .

Figure 9 .
Figure 9. Variation of the first six natural frequencies for a CF doubly-curved shell of revolution with hyperbolic meridian with variable thickness as a function of the exponent p.

Figure 10 .
Figure 10.First six mode shapes for a CF doubly-curved shell of revolution with hyperbolic meridian with variable thickness assuming 1 p = as exponent.

Figure 10 .
Figure 10.First six mode shapes for a CF doubly-curved shell of revolution with hyperbolic meridian with variable thickness assuming p = 1 as exponent.

Figure 11 .
Figure 11.Variation of the first three natural frequencies for a CFCF doubly-curved panel of translation with variable thickness as a function of the exponent p .

Figure 11 .
Figure 11.Variation of the first three natural frequencies for a CFCF doubly-curved panel of translation with variable thickness as a function of the exponent p.

Figure 12 .
Figure 12.First three mode shapes for a CFCF doubly-curved panel of translation with variable thickness assuming 1 p = as exponent.

Figure 12 .
Figure 12.First three mode shapes for a CFCF doubly-curved panel of translation with variable thickness assuming p = 1 as exponent.

Table 2 .
[90]t 20 frequency parameters Ω for a completely free (FF) annular plate with variable thickness[90].The LGDQ solution is obtained for several Higher-order Shear Deformation Theories (HSDTs) by setting I N = I M = 41 and I PN = I PM = 14.