An Equivalent Layer-Wise Approach for the Free Vibration Analysis of Thick and Thin Laminated and Sandwich Shells

The main purpose of the paper is to present an innovative higher-order structural theory to accurately evaluate the natural frequencies of laminated composite shells. A new kinematic model is developed starting from the theoretical framework given by a unified formulation. The kinematic expansion is taken as a free parameter, and the three-dimensional displacement field is described by using alternatively the Legendre or Lagrange polynomials, following the key points of the most typical Layer-wise (LW) approaches. The structure is considered as a unique body and all the geometric and mechanical properties are evaluated on the shell middle surface, following the idea of the well-known Equivalent Single Layer (ESL) models. For this purpose, the name Equivalent Layer-Wise (ELW) is introduced to define the present approach. The governing equations are solved numerically by means of the Generalized Differential Quadrature (GDQ) method and the solutions are compared with the results available in the literature or obtained through a commercial finite element program. Due to the generality of the current method, several boundary conditions and various mechanical and geometric configurations are considered. Finally, it should be underlined that different doubly-curved surfaces may be considered following the mathematical framework given by the differential geometry.


Introduction
The development of refined structural models or higher-order theories currently represents one of the most active areas for many researchers in the structural mechanics field.Their interest in this subject is mainly due to the fact that the use of advanced materials has shown inadequacy of the classical or first-order theories to capture the effective mechanical behavior of this class of materials, such as laminated composites or functionally graded materials [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20].In these particular circumstances, the fundamental hypotheses of a lower-order model could lead to inaccurate results since some effects are neglected or not even considered.A typical situation is the case of structures made of a combination of highly heterogeneous materials in the lamination scheme.As highlighted in the work by Librescu and Reddy [21], many refined approaches arose to predict a more coherent structural behavior.For the sake of completeness, the readers can find some examples of these approaches in [22][23][24][25][26][27][28][29][30][31][32][33].It should be noted that the aforementioned papers were mainly focused on plates and shells.
The insufficiency of classical plate and shell theories can also be highlighted in the dynamic analysis.In order to evaluate the natural frequencies of composite structures correctly, several higher-order shear deformation theories (HSDTs) were developed.In general, two different approaches are followed, namely, equivalent single layer (ESL) and layer-wise (LW) models (see Reddy [1] for a brief review).The former is able to study a composite structure by referring each geometric and mechanical parameter on the middle surface of the structure, which is taken as reference domain for the governing equations.In a similar manner, the degrees of freedom of the problem under consideration are evaluated on the middle surface, independently from the enrichment of the kinematic model [33][34][35][36][37][38][39][40][41][42][43][44][45][46].On the other hand, the latter approach takes into account each layer (or ply) that composes the structure separately from the others.In general, a LW model defines the kinematic expansion, and consequently the degrees of freedom, along the thickness of each layer.As a result, this approach is able to deal with continuous displacements but characterized by discontinuous derivatives at the interfaces [46][47][48][49][50][51][52][53][54][55][56].
The works [57][58][59] represent a turning point in the advancement of HSDTs.In fact, they provided the basis for the development of a unified formulation (UF), which is able to analyze and use several structural models within a single framework.As a matter of fact, the order of kinematic expansion can be arbitrarily chosen in order to obtain the desired theory that is more appropriate to deal with a particular physical phenomenon.In the review paper [57] a complete treatise about both ESL and LW approaches is presented.It is important to cite also the works [58,59], where the so-called zig-zag theories for layered structures were discussed.For the sake of completeness, it should be recalled that a zig-zag theory is actually an ESL model embedded with a specific function that allows to consider continuous through-the-thickness displacements characterized by different slope at the interface between two plies, similar to a LW theory [31][32][33].This aspect is commonly denoted as zig-zag effect and it is generally given by the choice of different materials along the transverse direction of the structure.For instance, this effect is particularly evident when a soft-core is included into the lamination scheme of a generic plate or shell (sandwich structure).In general, the Murakami's function is employed for this purpose.
Recently, a huge number of scientific papers based on the UF is available in the literature.For the sake of completeness, the works [60][61][62][63][64][65][66][67][68][69][70][71][72][73][74] can be considered as examples of higher-order ESL models.On the other hand, the papers [75][76][77][78][79][80][81][82][83] can be cited as the proof of the use of this unified formulation to develop higher-order LW approaches.It should be noticed that the UF was generalized by Demasi [84], who developed the so-called Generalized Unified Formulation (GUF).The research provided by Demasi represents the starting point for the development of more general theories.For this purpose, the work by D'Ottavio [85] must also be mentioned.He proposed an optimized approach able to consider the structure as the sum of packages of layers, each of them governed by a peculiar structural model.This approach was named Sublaminate Generalized Unified Formulation (S-GUF).Analogously, the paper by Fazzolari [86][87][88] and by Fazzolari and Banerjee [89] represent a significant extension of the GUF.Their works aimed to extend the contributions first given by Demasi in [84] for plates to the structural analysis of beams and shells.At this point, it is the authors' intention to apologize for the omissions in the present review of those papers that could be considered as important steps in the development of HSDTs.
The main aim of the present work is to introduce a new higher-order structural model based on the aforementioned unified formulation for the structural analysis of doubly-curved laminated composite plates and shells, which has some elements in common with both the ESL and LW approaches.According to the proposed methodology, the overall mechanical properties of the composites are computed on the shell middle surface, taking into account all the layers that compose the structure.It should be noticed that there is no limitation on the choice of the constituents, as well as their orientation and stacking sequence, due to the general features of this approach.The displacement field is written as a function of an arbitrary order of kinematic enrichment, whereas the Legendre and the Lagrange polynomials are employed to describe the kinematic expansion.The degrees of freedom of the problem are not defined on the shell middle surface as in a typical ESL model, but they are related to specific points placed along the thickness of the structure according to what is typically done for each Appl.Sci.2017, 7, 17 3 of 34 layer by the LW theory.The term Equivalent Layer-Wise (ELW) is introduced to define this approach.Once the governing equations are obtained, they are numerically solved by means of the Generalized Differential Quadrature (GDQ) method [90], whose main features are illustrated in detail by Tornabene et al. in the review paper [91].The accuracy, reliability and stability characteristic of this numerical tool can be checked considering the excellent results shown in the works [39][40][41][61][62][63][64][65][66][67][68][69][70][71][72][73][74][79][80][81][82][83][92][93][94][95][96][97][98][99][100][101][102], which have been obtained in different kinds of structural problems related to plates and shells.

Definition of the Geometry
As highlighted by Kraus in his book [103], each shell structure is a three-dimensional body which is bounded by two close surfaces, whose distance measured along the normal direction defines the thickness h of the shell.The middle surface is clearly equidistant from these two external surfaces and can be taken as reference domain for the governing equations.Its importance is evident, since its define the shape of the structure under consideration and rules the mechanical behavior of a shell structure if a two-dimensional structural model is considered.By means of the differential geometry, several kinds of shells can be studied in a general and unified manner.For instance, the same theoretical approach is used to define the geometry of doubly-curved, singly-curved and degenerate shells (or plates).With reference to Figure 1, in which a generic shell element is depicted, the position vector R(α 1 , α 2 , ζ) is introduced to identify each point of the three-dimensional medium as follows where the dimensionless quantity z = 2ζ/h ∈ [−1, 1] denotes the distance of the point P from its projection P on the reference domain.what is typically done for each layer by the LW theory.The term Equivalent Layer-Wise (ELW) is introduced to define this approach.Once the governing equations are obtained, they are numerically solved by means of the Generalized Differential Quadrature (GDQ) method [90], whose main features are illustrated in detail by Tornabene et al. in the review paper [91].The accuracy, reliability and stability characteristic of this numerical tool can be checked considering the excellent results shown in the works [39][40][41][61][62][63][64][65][66][67][68][69][70][71][72][73][74][79][80][81][82][83][92][93][94][95][96][97][98][99][100][101][102], which have been obtained in different kinds of structural problems related to plates and shells.

Definition of the Geometry
As highlighted by Kraus in his book [103], each shell structure is a three-dimensional body which is bounded by two close surfaces, whose distance measured along the normal direction defines the thickness h of the shell.The middle surface is clearly equidistant from these two external surfaces and can be taken as reference domain for the governing equations.Its importance is evident, since its define the shape of the structure under consideration and rules the mechanical behavior of a shell structure if a two-dimensional structural model is considered.By means of the differential geometry, several kinds of shells can be studied in a general and unified manner.For instance, the same theoretical approach is used to define the geometry of doubly-curved, singly-curved and degenerate shells (or plates).With reference to Figure 1, in which a generic shell element is depicted, the position vector  

,,
   R is introduced to identify each point of the three-dimensional medium as follows where the dimensionless quantity denotes the distance of the point P from its projection P on the reference domain. denote the lines of main curvature of the middle surface, whereas  is the normal direction that can be specified by the following outward unit normal vector   On the other hand, r(α 1 , α 2 ) is the position vector that describe the middle surface, which assumes different meanings according to the investigated geometry.A curvilinear orthogonal coordinates system O α 1 α 2 ζ must be defined on the shell middle surface.In particular, the coordinates α 1 , α 2 denote the lines of main curvature of the middle surface, whereas ζ is the normal direction that can be specified by the following outward unit normal vector n(α 1 , α 2 ) where the symbol "×" represents the vector product.For conciseness purposes, the notation r ,i = ∂r/∂α i , for i = 1, 2, is introduced.For the sake of completeness, it should be highlighted that the curvilinear coordinates α 1 , α 2 vary according to the considered structure.For example, they assume the following meaning α 1 = ϕ, α 2 = ϑ for a doubly-curved shell of revolution, whereas the notation α 1 = ϕ, α 2 = y is introduced for a singly-curved shell of translation.In the case of a flat panel, such as a rectangular plate, they can be assumed equal to α 1 = x, α 2 = y.This aspect is to underline the effectiveness of the differential geometry.In general, the position vector of any reference surface r(α 1 , α 2 ) is expressed as where f i (α 1 , α 2 ), for i = 1, 2, 3 are arbitrary functions needed to describe the surface under consideration.On the other hand, the symbols e i , for i = 1, 2, 3, denote the unit vectors of axes of the global reference system O x 1 x 2 x 3 as depicted in Figure 1.In the following sections, the position vector r(α 1 , α 2 ) will be specified for the various shell structures under consideration.Further details concerning these features, as well as a wide range of shell structures, can be found in the recent books by Tornabene et al. [2,3].In general, the curvilinear coordinates in hand must be bounded by specific values to define a finite domain.Without loss of generality, a three-dimensional shell is defined if the following limits are specified It should be noticed that the total thickness of the shell h is defined as the summation of the thickness of each layer (or lamina) h k , if a laminated composite material is taken into account.For a lamination scheme that consists of l layers as the one shown in Figure 1, the overall thickness is given by It is evident that the symbol k is used to define all the geometric and mechanical properties related to the k-th layer.The differential geometry provides also the definition of the well-known Lamè parameters A 1 (α 1 , α 2 ) and A 2 (α 1 , α 2 ), which can be evaluated once the position vector r(α 1 , α 2 ) is specified as follows where the symbol "•" is introduced to specify the scalar product.Finally, the main radii of curvature of the reference surface R 1 (α 1 , α 2 ) and R 2 (α 1 , α 2 ) can be computed using the definition shown below • n (7) in which the notation r ,ii = ∂ 2 r/∂α 2 i , for i = 1, 2, is employed.Expressions (7) are valid only if the curvilinear coordinates are orthogonal and principal.The three-dimensional effect of a shell structure is taken into account by the parameters H 1 (α 1 , α 2 , ζ) and H 2 (α 1 , α 2 , ζ), which are defined as follows

Higher-Order Equivalent Layer-Wise Approach
As already illustrated in the introduction, the well-known UF represents an extremely valid and efficient tool to analyze a huge variety of HSDTs, which are classified in general as Equivalent Single Layer (ESL) theories [61][62][63][64][65][66][67][68][69][70][71][72][73][74] and Layer-Wise (LW) models [79][80][81][82][83].The former approach establishes that all the mechanical and geometric parameters are evaluated on the shell middle surface, whereas each layer is analyzed independently from the others in the latter method.The main novelty introduced by the present work is a new higher-order formulation which keeps the main features of the ESL approaches but uses the kinematic expansion of the displacement field commonly employed in the LW models.For this purpose, the term Equivalent Layer-Wise (ELW) is introduced.Let us consider the three-dimensional displacements They can be defined in each point of the shell structure according to the following expressions (9) in which F τ = F τ (ζ), for τ = 0, 1, . . ., N, N + 1, specify the thickness functions related to the τ-th order of kinematic expansion.On the other hand, u 3 (α 1 , α 2 ), for τ = 0, 1, . . ., N, N + 1, are the degrees of freedom of the current model and denote the displacements of the section placed at the height ζ τ of the shell, as it can be noticed from Figure 2.This aspect represents the main difference in comparison with the ESL approach, where each degree of freedom is always evaluated on the shell middle surface.In other words, the τ-th order of kinematic expansion is strictly related to the number of points along the shell thickness (Figure 2).

Higher-Order Equivalent Layer-Wise Approach
As already illustrated in the introduction, the well-known UF represents an extremely valid and efficient tool to analyze a huge variety of HSDTs, which are classified in general as Equivalent Single Layer (ESL) theories [61][62][63][64][65][66][67][68][69][70][71][72][73][74] and Layer-Wise (LW) models [79][80][81][82][83].The former approach establishes that all the mechanical and geometric parameters are evaluated on the shell middle surface, whereas each layer is analyzed independently from the others in the latter method.The main novelty introduced by the present work is a new higher-order formulation which keeps the main features of the ESL approaches but uses the kinematic expansion of the displacement field commonly employed in the LW models.For this purpose, the term Equivalent Layer-Wise (ELW) is introduced.Let us consider the three-dimensional displacements ,, ,, They can be defined in each point of the shell structure according to the following expressions in which

 
FF    , for 0,1,..., , 1 NN   , specify the thickness functions related to the  -th order of kinematic expansion.On the other hand, , and , for 0,1,..., , 1 NN   , are the degrees of freedom of the current model and denote the displacements of the section placed at the height   of the shell, as it can be noticed from Figure 2.This aspect represents the main difference in comparison with the ESL approach, where each degree of freedom is always evaluated on the shell middle surface.In other words, the  -th order of kinematic expansion is strictly related to the number of points along the shell thickness (Figure 2).For the sake of clarity, it should be noticed that the displacements        For the sake of clarity, it should be noticed that the displacements u 3 linked to the order of expansion τ = 0 are the displacements of the points laying on the shell bottom surface (for ζ = −h/2).On the contrary, the displacements u  Similarly to the ESL approach, the (N + 1)-th order of kinematic expansion is always related to the well-known Murakami's function, which can be introduced to capture the zig-zag effect that can be noticed especially if a sandwich structure is analyzed.Having in mind Figure 1, in which a general lamination scheme is depicted, the Murakami's function Z = Z(ζ) can be defined as follows Z = (−1) k z k (10) The non-dimensional parameter z k ∈ [−1, 1] is assumed equal to the following expression in which ζ k represents the boundary coordinate of the k-th layer along the normal direction.Thus, if a soft-core effect must be studied, one gets F N+1 = Z.The reader can find further details concerning the use and the peculiar features of the Murakami's function in the works [57][58][59].Let us consider the thickness functions related to the other orders of kinematic expansions.The first possible choice consists in assuming them equal to the corresponding Legendre polynomials.In other words, they can be defined as follows in which the symbol L τ+1 is introduced to specify the τ-th order Legendre polynomial.They can be computed conveniently in the closed interval z k ∈ [−1, 1] by using the following recursive formula (as shown in [3]) for τ = 3, 4, . . ., N, assuming L 1 (z k ) = 1 and L 2 (z k ) = z k .For clarity purposes, let us consider a fourth order of kinematic expansion (N = 4).The displacement field (9) becomes 1 + F 2 u (2) 3 + F 2 u (2) if the Murakami's function is embedded in the model.Due to expressions (12) the first five Legendre polynomials are required As a consequence, the thickness functions needed for the kinematic expansion at issue take the following aspect It should be pointed out that F 5 is evidently the Murakami's function (F 5 = F N+1 = Z).If the Legendre polynomial are used, the coordinates ζ τ which define the kinematic expansion can be located by evaluating the roots of the Legendre polynomial of highest order.
Alternatively, the Lagrange polynomials can be used as thickness functions simply assuming F τ (ζ) = l τ+1 (ζ), in which l τ+1 (ζ) denotes the Lagrange polynomial of order N defined as follows for τ = 0, 1, . . ., N − 1, N. The coordinates ζ m , for m = 0, 1, . . ., N − 1, N, can be chosen taking into account an arbitrary distribution of points in the closed interval In the present paper, a uniform (or equally spaced) grid distribution is considered for the sake of simplicity.If the Lagrange polynomials are employed as thickness functions, one gets For the sake of completeness, it should be noted that the displacement field (9) without the Murakami's function is equivalent to the displacement approximation of the sampling surfaces (SaS) formulation presented in the works [104,105], if the Lagrange polynomials ( 17)-( 18) are employed as thickness functions.
As in the previous approach, the (N + 1)-th thickness function is still given by the Murakami's function (F N+1 = Z).Using both the Legendre polynomials and the Lagrange ones, the thickness functions have the following properties with τ, m = 0, 1, 2, . . ., N. As a consequence, the relations shown below can be easily deduced This means that the three-dimensional displacement evaluated at the height coordinate linked to the τ-th order of kinematic expansion coincides with the corresponding degree of freedom of the model, except for the Murakami's function.At this point it should be specified also that the kinematic expansion is not affected by either the lamination scheme or the placement of the various layers that compose the laminate, since this aspect is considered only when the elastic coefficients of the composite are evaluated as it will be shown in the following paragraphs.Finally, it should be noticed that the same kinematic model just illustrated is typically used for the LW approach to describe the displacement field of each lamina [79][80][81][82][83].A peculiar notation can be introduced now to specify univocally the various higher-order ELW models that are employed in the following.Each theory is identified by the maximum order of kinematic expansion N. In particular, the structural theories used in the current work are listed below for various orders of expansion where "EL" specifies that an ELW approach is considered, whereas "D" means that the governing equations will be written in terms of generalized displacements.If the Murakami's function must be embedded in the model, the following theories are obtained for several orders of kinematic expansion where "Z" stands clearly for the Murakami's function (zig-zag effect).
Once the kinematic model is defined, both the constitutive relations and the motion equations are formally equal to the ones used for the ESL approach [61][62][63][64][65][66][67][68][69][70][71][72][73][74].Thus, in the following only the main aspects of the model are presented for conciseness purposes.At this point, the vector u (τ) = u (τ) (α 1 , α 2 ) that collects the degrees of freedom for each order τ of kinematic expansion can be conveniently introduced for τ = 0, 1, . . ., N, N + 1.The generalized strain components for each order τ of kinematic expansion can be collected in the corresponding algebraic vector ε (τ) = ε (τ) (α 1 , α 2 ) defined as follows The generalized strains can be directly related to the degrees of freedom of the problem for τ = 0, 1, . . ., N, N + 1 according to the following equation in vector notation The differential operator D Ω assumes the following matrix notation where A 1 , A 2 are the Lamè parameters (6), and R 1 , R 2 the principal radii of curvature (7) of the shell middle surface.
The complete treatise concerning the definition of the generalized strain components, as well as further details about their meaning, can be found in the books by Tornabene et al. [2,3].As highlighted (27) Due to the duality between strains and stresses, the generalized stresses (or stress resultants) can be related directly to the degrees of freedom of the model for τ = 0, 1, . . ., N, N + 1 according to the following equation where the constitutive operator A (τη) is introduced for each order of kinematic expansion (for τ, η = 0, 1, 2, . . ., N, N + 1).It should be noticed that A (τη) represents the stiffness matrix of the considered shell structure.In particular, if a laminated composite structure made of l orthotropic layers is considered, the constitutive operator takes the following aspect A (τη) Each term collected in the stiffness matrix can be evaluated by means of the following definitions for τ, η = 0, 1, 2, . . ., N, N + 1, n, m = 1, 2, 3, 4, 5, 6, and p, q = 0, 1, 2. It should be noticed that the symbols τ, η specify the thickness function order, whereas τ, η denote that the derivatives of the corresponding thickness functions F τ , F η with respect to ζ must be evaluated.It is pointed out that this evaluation can be easily performed since the definition of the thickness functions is shown in extended form in ( 12) and ( 18), respectively.Quantities B (k) nm are defined as follows for n, m = 1, 2, 3, 6 and for n, m = 4, 5.This distinction is performed to introduce the shear correction factor κ = 1/χ that has to be included in the model if the chosen order of kinematic expansion does not provide a parabolic profile of the through-the-thickness shear stresses.Typically, a constant value is used to define the shear correction factor κ. In the present paper, the value χ = 1.2 is used for this purpose.Nevertheless, it should be noticed that a parabolic function can be chosen as shear correction function as illustrated in the work by Tornabene and Reddy [92].On the other hand, quantities nm , for n, m = 1, 2, 3, 4, 5, 6, represent the elastic constant of the material of the k-th layer.They can assume different meanings depending on the case.If the strain along the shell thickness is neglected, the plane-stress-reduced elastic coefficients are needed (E nm ).On the contrary, if the stretching effect is taken into account the non-reduced elastic coefficients must be used (E nm ).The trace over these symbols means that these quantities must be evaluated in the local reference system O α 1 α 2 ζ.In other words, the constitutive equations of each lamina must be transferred into the geometric coordinate system by means of the proper transformation laws which take into account the different orientation that can be given to an orthotropic medium.In fact, since each layer can be arbitrarily oriented with respect to the laminate reference system as shown in Figure 1, in which the symbols α, β, γ, δ are introduced for this purpose, the constitutive relations of each lamina must be converted into the local reference system O α 1 α 2 ζ.Further details concerning these relations can be found in the book by Reddy [1].For the sake of completeness, it should be specified that in the present work the mechanical properties of each layer k are fully described by means of the engineering constants of the materials, which are 23 for an orthotropic medium, and E (k) , G (k) , ν (k) for an isotropic one.In any case, the lamination scheme is identified by the notation (α, β, . . . ,γ, . . . ,δ), in which α, β, γ, δ represent the orientation of each ply.It should be recalled that in the current approach each layer is assumed to be linear elastic.In any circumstance, the use of both the shear correction factor and the plane stress-reduced elastic coefficients is specified by adding a proper notation to the acronym of the structural theory.In particular, the superscript χ = 1.2 is added to specify the use of the constant value for the shear correction factor, whereas the subscript RS stands for "Reduced Stiffness" and means that the reduced elastic coefficients are employed.
The motion equations are obtained by means of the well-known Hamilton's variational principle.For the sake of conciseness, only the equations at issue are shown.For each order of kinematic expansion τ, the free vibration problem is governed by the following matrix equation The equilibrium operator D * Ω takes the following aspect in which A 1 , A 2 represent the Lamè parameters defined in (6), and R 1 , R 2 are the principal radii of curvature of the shell middle surface shown in (7).On the contrary, the matrix M (τη) collects the inertia terms I (τη) for every order of kinematic expansion τ, η = 0, 1, 2, . . ., N, N + 1 as follows If ρ (k) denotes the mass density of the k-th ply, the definition shown below is used to evaluate the inertia masses I (τη) of the laminated composite in hand It is clear that .. u (η) represents the generalized acceleration component vector defined as follows . .
for η = 0, 1, 2, . . ., N, N + 1.By inserting definition (28) into the motion equation shown in (33), the governing system can be conveniently written as a function of the degrees of freedom of the problem.As a consequence, one gets where the so-called fundamental operator L (τη) = D * Ω A (τη) D Ω is introduced.In matrix form, it assumes the following aspect For conciseness purposes, the definition of each term L (τη) f g , for f , g = 1, 2, 3, is omitted in the present paper.Nevertheless, the reader can find their complete meaning in the work by Tornabene et al. [61].It should be noticed that they are formally equal to the ones employed in the ESL approach.Relation (38) represents the fundamental nucleus of the ELW approach and it identify a set of three differential equations for each order of kinematic expansions.Thus, 3(N + 2) equations must be solved to obtain the solution of this structural problem.In order to fully characterize the problem at issue, the proper set of boundary conditions must be enforced along each external edge.For a clamped edge (denoted by letter "C"), these conditions are written in terms of displacements.On the contrary, specific prescription are given to the generalized stress resultants for a free edge (indicated by letter "F").Let us consider Figure 1 to identify each external edge.In particular, if an edge is defined by These edges are the Southern (S) and the Northern (N) ones in Figure 1.On the other hand, if the coordinates of an edge are This is clearly the case of the Western (W) and Eastern (E) edges.In the present paper, the boundary conditions are specified for the panel under consideration following the sequence WSEN.For instance, this means that the acronym FCCF specifies that the Western and the Northern edges are free, whereas the other two are completely clamped.

Generalized Differential Quadrature Method
In the present paper, the governing equations shown in matrix form in (38) are solved by means of the GDQ method.According to this numerical technique, it is possible to approximate the n-th derivative of a sufficiently smooth function at a generic point within the reference domain as a weighted linear sum of the values that the function assumes in some chosen grid points.Starting from this definition, it is clear that the strong form of the governing equations is solved.Let us consider a two-dimensional domain, as the one needed to describe any doubly-curved surface, in which a generic function f (α 1 , α 2 ) is defined.Three types of derivative can be computed at a generic point α 1i , α 2j for i = 1, 2, . . ., I N and j = 1, 2, . . ., I M : the n-th derivative with respect to the first coordinate α 1 the m-th derivative with respect to the second coordinate α 2 and the mixed derivative of order (n + m) Appl.Sci.2017, 7, 17 13 of 34 It is pointed out that I N , I M denote the overall number of discrete points within the reference domain.The symbols ς α 2 (jl) are introduce to specify the weighting coefficients for the derivatives at issue.It is clear that the key points of the present numerical approach are the choice of a proper discrete grid point distribution and the evaluation of the weighting coefficients.If the Lagrange polynomials are chosen to approximate the function, there is no restriction on the choice of the grid points.This aspect represents one of the main advantages of the GDQ method.In addition, the weighting coefficients can be easily evaluated if the recursive formulae given by Shu are employed [90].Further details concerning this technique and its major applications can be found in the review paper by Tornabene et al. [91].As far as the choice of the grid points is concerned, several distributions can be used, as highlighted in the work [91].Some examples of grid distributions typically employed in many structural problems are listed in Table 1.It should be recalled that a discrete grid must be applied along the two principal coordinates α 1 , α 2 as shown in the following relations Table 1.Eight different kinds of grid distributions.The indices k and N specify respectively i and I N , when α 1 direction is taken into account.On the other hand, they correspond to j and I M when α 2 direction is considered.Finally, if the grid distribution has to be applied along the coordinate ζ, one gets k = m and N = I T .It should be noted that L N denotes the Legendre polynomials.

Generalized Integral Quadrature Method
A numerical approach has to be introduced to evaluate the integrals needed to compute the elastic constants nm (pq) defined in (30).In the present paper, the Generalized Integral Quadrature (GIQ) method is employed for this purpose.In fact, a numerical solution to the integrals in (30) must be found since it is impossible to solve them analytically due to their definition.According to this approach, the integral of a smooth function Appl.Sci.2017, 7, 17 14 of 34 can be evaluated as a weighted linear sum of the values that the function f (ζ) assumes in each point of the reference domain.In other words, this definition can be written as follows where I T represents the number of discrete points defined within the reference domain, whereas w ek , w lk are the weighting coefficients for the integration.These last quantities can be computed directly from the definitions of the corresponding weighting coefficients of the first order derivatives, as illustrated in depth in the review paper by Tornabene et al. [91].It is clear that the GIQ method shares the same fundamentals with the GDQ method.Thus, the same recursive formulae given by Shu can be used once again to evaluate the weighting coefficients at issue, once an arbitrary grid distribution is chosen to discretize the reference domain.It should be noticed that differently from the well-known Gaussian integration schemes, this method is able to compute the weighting coefficients without any restriction on the position of the points within the reference domain.Therefore, all the different kinds of grid distributions listed in Table 1 can be used with this aim.For the sake of clarity, it should be specified that the one-dimensional domain in which the integrals in (30) must be evaluated is defined along the coordinate ζ of the local reference system.In other words, the integration procedure is performed through the shell thickness.As a consequence, the discrete nodes in which the functions are evaluated are placed following the definition shown below for m = 1, 2, . . ., I T .The values of r m , for m = 1, 2, . . ., I T , can be found in Table 1, as well.For the sake of simplicity, in the present paper the Chebyshev-Gauss-Lobatto grid distribution is employed in any numerical test with I T = 25 due to its accuracy and stability features.

Free Vibration Analysis
The free vibrations of a doubly-curved laminated composite shell can be evaluated starting from the fundamental equation (38) by means of the GDQ method.A solution to the fundamental system can be found in the following form (49) in which the vector T collects the mode shapes for each order of kinematic expansion τ.On the other hand, the symbol ω is introduce to specify the corresponding circular frequencies, which allow to compute the natural frequencies as f = ω/2π.Once the first and second order derivatives with respect to the time variable t of definition (49) are evaluated and substituted into relation (38), one gets for τ = 0, 1, 2, . . ., N, N + 1.It should be recalled that equation ( 50) is the fundamental equation system which rules the free vibrations of any laminated composite shell.The system (50) can be easily solved by means of the GDQ method illustrated in the previous section.Once each derivative is discretized following definitions ( 42)-( 44), the governing system can be conveniently written in its discrete form as follows where K denotes the global stiffness matrix, M the inertia matrix, and δ the mode shape vector.Relation ( 51) is clearly a generalized linear eigenvalue problem, whose eigenvalues are the circular frequencies by definition.It is clear that the corresponding eigenvectors allow to define the mode shapes of the structures under consideration.The well-known kinematic condensation of non-domain degrees of freedom, the quantities related to the internal points of the domain (denoted by the subscript "d") can be separated from the ones linked to the boundaries (specified by the subscript "b").In this manner, the initial problem is noticeably simplified and can be written as follows where K denotes the global stiffness matrix, M the inertia matrix, and δ the mode shape vector.Relation ( 51) is clearly a generalized linear eigenvalue problem, whose eigenvalues are the circular frequencies by definition.It is clear that the corresponding eigenvectors allow to define the mode shapes of the structures under consideration.The well-known kinematic condensation of nondomain degrees of freedom, the quantities related to the internal points of the domain (denoted by the subscript " d ") can be separated from the ones linked to the boundaries (specified by the subscript " b ").In this manner, the initial problem is noticeably simplified and can be written as follows from which the following relation is carried out assuming that I is the identity matrix.The vector δ that collects the degrees of freedom is clearly subdivided into  represents the vector related to the internal degrees of freedom.The same partition is performed for both the stiffness matrix K and the mass matrix M .The solution of the system (53) allows to obtain the circular frequencies and consequently the natural frequencies.

Numerical Results
In the present section, several applications are presented to show the validity of both the present theoretical approach and the numerical methods for the evaluation of the free vibrations of various laminated composite shells.For this purpose, a preliminary convergence analysis is performed for two different mechanical configurations.Then, six shell structures are investigated.The results are compared with the natural frequencies available in the literature and the solutions obtained by means of a Finite Element model (FEM).The GDQ analysis is performed through a MATLAB code [106], whereas the commercial software Strand 7 is employed to obtain the FEM solution.It should be noticed that a specific tool has been developed to export the geometry of the considered doublycurved shells into the commercial finite element code, as specified in the work [98].For the sake of completeness, it should be specified that a different FEM model is used according to the investigated structure.For thick plates and shells, a three-dimensional model made of 20-node hexahedrons ("Hexa20") is employed.On the other hand, a two-dimensional model made of 9-node quadrilateral elements ("Quad9") is used for thin structures.The notations 3D-FEM and 2D-FEM are introduced to denote these solutions, respectively.

Convergence Analysis
A flat panel is considered in order to perform the convergence analysis.The length of each edge of the plate is equal to , whereas the thickness is given by 0.1m h  .Two different lamination schemes are taken into account in order to prove the capability of the current approach to deal with different mechanical configurations.In particular, an isotropic plate made of Aluminum is considered firstly.Then, the same structure is taken with   90 / 0 / 90 / 0 / 90 as lamination scheme.In (52) from which the following relation is carried out assuming that I is the identity matrix.The vector δ that collects the degrees of freedom is clearly subdivided into δ b and δ d , respectively.For the sake of clarity, δ b denotes the vector of the degrees of freedom along the boundary, whereas δ d represents the vector related to the internal degrees of freedom.The same partition is performed for both the stiffness matrix K and the mass matrix M.
The solution of the system (53) allows to obtain the circular frequencies and consequently the natural frequencies.

Numerical Results
In the present section, several applications are presented to show the validity of both the present theoretical approach and the numerical methods for the evaluation of the free vibrations of various laminated composite shells.For this purpose, a preliminary convergence analysis is performed for two different mechanical configurations.Then, six shell structures are investigated.The results are compared with the natural frequencies available in the literature and the solutions obtained by means of a Finite Element model (FEM).The GDQ analysis is performed through a MATLAB code [106], whereas the commercial software Strand 7 is employed to obtain the FEM solution.It should be noticed that a specific tool has been developed to export the geometry of the considered doubly-curved shells into the commercial finite element code, as specified in the work [98].For the sake of completeness, it should be specified that a different FEM model is used according to the investigated structure.For thick plates and shells, a three-dimensional model made of 20-node hexahedrons ("Hexa20") is employed.On the other hand, a two-dimensional model made of 9-node quadrilateral elements ("Quad9") is used for thin structures.The notations 3D-FEM and 2D-FEM are introduced to denote these solutions, respectively.Table 2. Mechanical properties of the materials employed in the numerical analyses.

Convergence Analysis
A flat panel is considered in order to perform the convergence analysis.The length of each edge of the plate is equal to L x = L y = 2m, whereas the thickness is given by h = 0.1m.Two different lamination schemes are taken into account in order to prove the capability of the current approach to deal with different mechanical configurations.In particular, an isotropic plate made of Aluminum is considered firstly.Then, the same structure is taken with (90/0/90/0/90) as lamination scheme.In this circumstance, each ply has the same thickness and is made of the same material (Graphite-Epoxy).For conciseness purposes, the mechanical properties of the materials employed in this section are summarized in Table 2.As far as the boundary conditions are concerned, the structure under consideration is fully clamped along each edge (CCCC).Figure 3 shows the convergence characteristics of the first three natural frequencies as a function of the total number of grid points I N = I M varying the grid distribution.It can be easily noticed that the current approach reaches the maximum level of convergence with a number of point equal to I N = I M = 11, for both the two lamination schemes.In a similar manner, the various grids lead to similar tendency.It should be noticed that the graphs depicted in Figure 3  this circumstance, each ply has the same thickness and is made of the same material (Graphite-Epoxy).
For conciseness purposes, the mechanical properties of the materials employed in this section are summarized in Table 2.As far as the boundary conditions are concerned, the structure under consideration is fully clamped along each edge (CCCC).Figure 3 shows the convergence characteristics of the first three natural frequencies as a function of the total number of grid points NM II  varying the grid distribution.It can be easily noticed that the current approach reaches the maximum level of convergence with a number of point equal to 11 NM II  , for both the two lamination schemes.In a similar manner, the various grids lead to similar tendency.It should be noticed that the graphs depicted in Figure 3 are related to the   The same graphs could be achieved also for higher-order theories, whose results are just summarized in Table 3 (for the isotropic case) and in Table 4 (for the laminated one) for conciseness purposes.It should be noticed that the natural frequencies presented in Tables 3 and 4 are obtained considering the Cheb-Gau-Lob grid distribution with I N = I M = 21.In addition, the 3D-FEM solution (9000 brick elements "Hexa20") is also included in these tables.The GDQ solution is evidently in good agreement with the one used as a reference.For the sake of completeness, the analysis is carried out considering HSDTs with and without the Murakami's function for the laminated composite plate (Table 4).In this circumstance, since the mechanical properties of each layer are the same, the benefits introduced by embedding this function are negligible.It should be noticed also that the reduced elastic coefficients are used only for a first order kinematic expansion (N = 1).On the other hand, a constant value of the shear correction factor equal to χ = 1.2 is used for both the first-order and second order theories (N = 1, 2) due to the reasons illustrated in the previous sections.

Free vibrations of Laminated Composite Shells
In this paragraph, the free vibrations are investigated for the six structures depicted in Figure 4, in which the GDQ grid and local coordinate reference systems are shown for each geometry.
The first two laminated shells are taken from the work by Wang et al. [13] and they are employed to compare the results given by the current approach with the ones obtained by using the Fourier-Ritz Method.The reference surface of the toro-circular panel at issue (Figure 4a,b) can be described by the following position vector r(ϕ, ϑ) where ϕ, ϑ are the curvilinear coordinates of the surface.R = 3m is the radius of the circular meridian and R b represents an offset from the revolution axis.The first toro-circular panel is obtained by setting R b = 1.5m, whereas the second one can be achieved by using R b = −1.5m.In both these circumstances, the surfaces are bounded assuming the limitations ϕ ∈ [π/3, 2π/3] and ϑ ∈ [−π/3, π/3], whereas the thickness of the shells is h = 0.2m.In order to perform the validation test, the following mechanical properties are employed for each orthotropic ply:  5 and 6, the excellent agreement with the reference solution can be easily noticed independently from the HSDTs, for each boundary condition and lamination scheme.It should be noted that in the reference work by Wang [13], a modified Fourier series approach in conjunction with a Ritz method is used to obtain a semi-analytical solution based on the First-order Shear Deformation Theory (FSDT).The next structure is a CCFC singly-curved cylindrical panel whose profile is defined by a catenary section (Figure 4c).It can be obtained by sliding a catenary curve on an inclined straight line, whose slope is equal to α = π/12.The reference surface can be described by the following position vector r(ϕ, y) r(ϕ, y) = (d arcsinh(tan ϕ) + y sin α sin ϕ) e 1 − y cos α e 2 + + (d(cos h(arcsinh(tan ϕ)) − 1) − y sin α cos ϕ)e 3 (55) in which ϕ, y are the curvilinear coordinates of the surface, assuming ϕ ∈ [−1.131728,The first 10 natural frequencies of this structure are shown in Table 7 for different higher-order ELW models.As in the previous tests, each order of kinematic expansion is taken with and without the Murakami's function to investigate its contribution.It is easy to notice that in the present case all the theories give the same results, independently from this function.The accuracy of the numerical technique is proven also by the comparison with the 3D-FEM solution (8000 brick elements "Hexa20"), which is nearly equal to the proposed ones.It should be recalled that the same considerations illustrated above concerning the use of both the shear correction factor and the reduced elastic coefficients are still valid.Finally, the first nine mode shapes for the shell structure at issue are depicted in Figure 5.
The following structure is a CCCC doubly-curved revolution shell, whose reference surface is generated by the rotation about the revolution axis of a branch of a catenary (Figure 4d).r(ϕ, ϑ) = a cosh arcsinh 1 tan ϕ cos ϑ e 1 − a cosh arcsinh 1 tan ϕ sin ϑ e 2 + + a arcsinh 1  tan ϕ e 3 (56) where ϕ, ϑ are the curvilinear coordinates of the surface defined in the intervals ϕ ∈ [2.626124, 0.515469] and ϑ ∈ [−π/3, π/3], with a = 0.75 m.The lamination scheme is chosen to define a so-called sandwich structure in which the core exhibits mechanical properties that are considerably different from the ones that characterize the external face-sheets, especially along the transverse direction.In other words, if the stiffness of the constituent materials along the transverse direction is really different the zig-zag effect is generated.As a consequence, the Murakami's function is embedded in the kinematic model to represent the effective behavior of the structure.The lamination scheme here is given by   .With reference to the mechanical properties of Table 2, it is easy to notice that  3 (1) =  3 (3) ≫  3 (2) .Thus, the shell falls into the class of sandwich structures.The first 10 natural frequencies are shown in Table 8 and are evaluated by using the Cheb-Gau-Lob grid distribution with 31 NM II  .In the same table, the 3D-FEM solution is also reported for comparison purposes.It should be noticed that the first two orders of kinematic expansion ( 1, 2 N  ) are taken into account with and without the shear correction factor.From Table 8 it can be noted that the use of the shear correction factor 1.2   leads to a solution that are slightly different from the others.On the other hand, if the correction is not applied (which means 1   ) the natural frequencies are similar to the ones obtained for higher-order theories, as well as the ones given by the 3D-FEM (11560 brick elements "Hexa20").This aspect has been noticed by the authors in their previous work [69], where the same considerations were made for the corresponding theories based on the ESL approach embedded with the Murakami's function.Thus, the shear correction factor should be neglected in those structural theories with the Murakami's function, if a sandwich structure is analyzed.Finally, the related first nine mode shapes can be found in Figure 6.In other words, if the stiffness of the constituent materials along the transverse direction is really different the zig-zag effect is generated.As a consequence, the Murakami's function is embedded in the kinematic model to represent the effective behavior of the structure.The lamination scheme here is given by (45/Core/ − 45), assuming that the external orthotropic layers have equal thickness h 1 = h 3 = 0.02 m and are made of Graphite-Epoxy.On the other hand, the isotropic core is made of Ceramic Foam and it is defined by a thickness equal to h 2 = 0.06 m.With reference to the mechanical properties of Table 2, it is easy to notice that E (1) 3 .Thus, the shell falls into the class of sandwich structures.The first 10 natural frequencies are shown in Table 8 and are evaluated by using the Cheb-Gau-Lob grid distribution with I N = I M = 31.In the same table, the 3D-FEM solution is also reported for comparison purposes.It should be noticed that the first two orders of kinematic expansion (N = 1, 2) are taken into account with and without the shear correction factor.From Table 8 it can be noted that the use of the shear correction factor χ = 1.2 leads to a solution that are slightly different from the others.On the other hand, if the correction is not applied (which means χ = 1) the natural frequencies are similar to the ones obtained for higher-order theories, as well as the ones given by the 3D-FEM (11560 brick elements "Hexa20").This aspect has been noticed by the authors in their previous work [69], where the same considerations were made for the corresponding theories based on the ESL approach embedded with the Murakami's function.Thus, the shear correction factor should be neglected in those structural theories with the Murakami's function, if a sandwich structure is analyzed.Finally, the related first nine mode shapes can be found in Figure 6.The next structure is a CFCF (clamped-free-clamped-free) doubly-curved panel of translation, obtained by sliding an elliptic profile (defined by the semi diameters a, b) over a parabola (described through the parameters a, b, d).It should be noticed that this surface is obtained by moving the elliptic curve along the parabolic curve, keeping the plane containing the ellipse always orthogonal to the parabola (Figure 4e).The position vector r(α 1 , α 2 ), which defines the doubly-curved translational surface in hand, can be written as a function of the curvilinear coordinates α 1 , α where the quantities expressed as a function of the coordinate α 1 are related to the parabola (curve which the other curve is sliding on), whereas the geometric parameters that depend on α 2 are linked to the ellipse (curve that is moving on the other one).As far as the parabolic branch is concerned, one gets x where k α 1 = a 2 − d 2 /b is the characteristic parameter of the curve at issue, assuming the following values a = 2 m, b = 1 m, c = −2 m, d = 0 m.Further details concerning the geometric features of the parabolic arch, as well as the meaning of each parameter, can be found in the work [68].On the other hand, the following terms are related to the elliptic curve in which the characteristic parameter of the ellipse is given by k α 2 = a/b, with a = 1 m and b = 1.5 m.The doubly-curved panel of translation is defined by setting α 1 ∈ [−0.785398, 0.785398] and α 2 ∈ [−5π/12, 5π/12].The total thickness of the structure is equal to h = 0.1m and the lamination scheme is (0/30/60/90).Each ply made of Graphite-Epoxy has the same thickness.The GDQ solution is obtained by using the Cheb-Gau-Lob grid distribution and the number of grid points is equal to I N = 35 and I M = 25.The natural frequencies are shown in Table 9 for several higher-order ELW approaches, with and without the Murakami's function.In the same table, the 3D-FEM solution is also given for comparison purposes.The use of the shear correction factor and the reduced elastic coefficients follows the same considerations illustrated before.It should be noticed that the analyses are carried out considering both the two types of Thickness functions defined in the previous sections.
In the first row, the natural frequencies are obtained by means of the Legendre Polynomials, whereas in the second raw the results are related to the Lagrange ones.It can be noticed that the solutions are not affected by this choice, since the natural frequencies are the same independently from the approach.Even in this circumstance, the solutions are in excellent agreement with the ones provided by the FEM (9600 brick elements "Hexa20").The first nine mode shapes are depicted in Figure 7.It is important to noticed that the deformed shapes follow the boundary conditions applied to the structure.It should be noted that all the structures considered until now are included in the category of thick and moderately thick shells.It can be noticed from the proposed results match perfectly even in this case.The Murakami's function is not taken into account due to the reduced value of thickness that allows to neglect the possible zig-zag effect.The first nine mode shapes are depicted in Figure 8.This last example has clearly proven also the accuracy of the present method for the free vibration analysis of thin shells.It can be noticed from the proposed results match perfectly even in this case.The Murakami's function is not taken into account due to the reduced value of thickness that allows to neglect the possible zig-zag effect.The first nine mode shapes are depicted in Figure 8.This last example has clearly proven also the accuracy of the present method for the free vibration analysis of thin shells.
Table 10.First 10 natural frequencies ([Hz]) for a CCCC thin laminated helicoidal surface with (0/90/0/90) as lamination scheme for several higher-order ELW approaches.The Cheb-Gau-Lob grid distribution is used for the numerical solution assuming

Comparison with the LW and ESL Approaches
In this last paragraph, the results obtained by means of the proposed ELW approach are compared to the ones achievable through the ESL and LW models.For this purpose, the laminated conical shell described in the works [61,83] is considered.The same work can be taken as a reference for both the geometric and mechanical properties of the laminated composite shell at issue.For completeness purposes, the conical surface in hand is given by the following position vector The GDQ method is employed to obtain the numerical solution even in this circumstance, for all the considered approaches.Since various orders of kinematic expansion are taken into account, the comparison should be performed by considering the corresponding higher-order displacement field for each approach (ESL, ELW, and LW).In particular, the ESL solutions are taken from the paper [61], whereas the work [83] is the reference for the LW ones.The same grid distributions and total values of discrete points specified in the works [61,83] are used also in the present paper.The natural

Comparison with the LW and ESL Approaches
In this last paragraph, the results obtained by means of the proposed ELW approach are compared to the ones achievable through the ESL and LW models.For this purpose, the laminated conical shell described in the works [61,83] is considered.The same work can be taken as a reference for both the geometric and mechanical properties of the laminated composite shell at issue.For completeness purposes, the conical surface in hand is given by the following position vector r(x, ϑ) = R 0 (x) cos ϑe 1 − R 0 (x) sin ϑe 2 + x sin ϕe 3 (63) where R 0 (x) = R + x cos ϕ.The constant angle ϕ = π/3 represents the slope of the straight meridian profile, whereas R = 1 m is the radius measured at the top of the shell.The two-dimensional domain is bounded by the limitations x ∈ [0, L x ] and ϑ ∈ [ϑ 0 , ϑ 1 ], with L x = 3 m, ϑ 0 = 0 and ϑ 1 = 2π.The shell is made of two layers of Graphite-Epoxy characterized by the same thickness value h 1 = h 2 = 0.15 m.On the other hand, the lamination scheme is given by (−45/45).As far as the boundary conditions are concerned, the conical shell is completely clamped along the top edge, whereas the lower circumference is free.The GDQ method is employed to obtain the numerical solution even in this circumstance, for all the considered approaches.Since various orders of kinematic expansion are taken into account, the comparison should be performed by considering the corresponding higher-order displacement field for each approach (ESL, ELW, and LW).In particular, the ESL solutions are taken from the paper [61], whereas the work [83] is the reference for the LW ones.The same grid distributions and total values of discrete points specified in the works [61,83] are used also in the present paper.The natural frequencies of the conical shell are shown in Table 11, where the solutions obtained through a three-dimensional finite element model are presented too.It should be noted that both the ESL and ELW models are taken with and without the Murakami's function.In general, the ELW approach provides values of natural frequencies placed between the ones given by the corresponding ESL and LW models.In other words, the results obtained by the ELW approach are closer to the LW solutions, as well as to the 3D-FEM, than the corresponding ESL ones.Nevertheless, the computational cost of the ELW model is reduced if compared to the one that characterizes the LW approach, since the proposed methodology does not require to consider each layer separately.For the sake of completeness, the reader can refer to the papers [61,83] to understand the meaning of the ESL and LW models introduced in Table 11.

Conclusions
A new higher-order approach is presented to analyze the natural frequencies of thin and thick laminated composite shells.The method, named Equivalent Layer-Wise (ELW) method, is used to evaluate the overall mechanical and geometric properties directly on the shell middle surface, as commonly done in the Equivalent Single Layer (ESL) approaches, and the displacement field has been described by using the same kinematic expansions employed in the most typical Layer-Wise (LW) models, but without considering each layer separately.The governing equation are solved numerically by the GDQ method.This new approach has been validated by comparison with the results available in the literature and those obtained using a commercial finite element code.It has been illustrated that the solutions are in excellent agreement with the reference ones for several boundary conditions and various mechanical configurations.When compared to the corresponding three-dimensional finite element model, the present methodology requires a lower number of degrees of freedom to reach the same level of accuracy.In addition, since the differential geometry is employed to describe the reference doubly-curved surface, several geometric shapes have been easily obtained and analyzed.
In fact, the domain has been defined using only one element, whereas the FEM models required many brick or plate elements to describe the same geometry accurately.It should be noted that the present methodology provides a physical meaning to the degrees of freedom of the problem for each order of kinematic expansion, different from the ESL approach, in which the generalized displacements related to higher-order expansions represent a mathematical parametrization of these quantities.Furthermore, the structural analysis of a composite structure has not required the same computational resources needed by the LW methods.Thus, the computational time has been noticeably reduced when compared to the LW models.It should be noted also that the use of the Murakami's function has been also investigated by taking into account a sandwich structure.The linear static analysis of laminated composite shell structures can be seen as a future development of the present approach.Finally, it should be highlighted that the presented methodology could be used to deal with some typical aerospace applications which require an efficient tool to analyze composite structures characterized by complex geometries, with different material and thickness distributions.

Figure 1 . 12 O
Figure 1.Representation of a generic doubly-curved shell element, edge identification and lamination scheme.On the other hand,   12 ,  r is the position vector that describe the middle surface, which assumes different meanings according to the investigated geometry.A curvilinear orthogonal coordinates system 12 O     must be defined on the shell middle surface.In particular, the

Figure 1 .
Figure 1.Representation of a generic doubly-curved shell element, edge identification and lamination scheme.

Figure 2 .
Figure 2. Kinematic expansion and degrees of freedom related to the order for the Equivalent Layer-Wise (ELW) approach.

0 
 are the displacements of the points laying on the shell bottom surface (for

Figure 2 .
Figure 2. Kinematic expansion and degrees of freedom related to the order for the Equivalent Layer-Wise (ELW) approach.
related to the τ = N order of kinematic expansion and represent the displacements of the top surface of the shell.

Figure 3 .
Figure 3. Convergence characteristics of the first three frequencies for a completely clamped (CCCC) square plate ( 2m xy LL  ) of thickness 0.1m h  for two different lamination schemes as a function of the number of grid points NM II  , varying the grid distribution.The ELW theory is the

Figure 3 .
Figure 3. Convergence characteristics of the first three frequencies for a completely clamped (CCCC) square plate (L x = L y = 2 m) of thickness h = 0.1m for two different lamination schemes as a function of the number of grid points I N = I M , varying the grid distribution.The ELW theory is the ELD1 χ=1.2 RS in both the circumstances.(a) and (b) are related to the first frequency, (c) and (d) to the second one and (e) and (f) to the third one.
GPa, G 23 = 5 GPa, ν 12 = ν 13 = ν 23 = 0.23, and ρ = 1500 kg/m 3 .The natural frequencies are written in their non-dimensional form Ω = 2π f R ρ/E 2 in Table5for the toro-circular panel with R b = 1.5 m, whereas the results shown in Table6are related to the second surface (R b = −1.5 m).Several lamination schemes and boundary conditions are taken into account.As far as the GDQ solutions are concerned, the natural frequencies are computed by using the Cheb-Gau-Lob grid distribution assuming I N = 21 and I M = 31 for the first case, and I N = I M = 25 for the second one.As in the previous example, various orders of kinematic expansion are considered for the ELW approach without embedding the Murakami's function.From Tables

Figure 4 .
Figure 4. Definition of the geometry of six shell structures: Generalized Differential Quadrature (GDQ) grid and local coordinate reference system representation.

Figure 4 .
Figure 4. Definition of the geometry of six shell structures: Generalized Differential Quadrature (GDQ) grid and local coordinate reference system representation.

Figure 6 .
Figure 6.First nine mode shapes for a CCCC sandwich panel of revolution with   45 / Core / 45  as

Figure 6 .
Figure 6.First nine mode shapes for a CCCC sandwich panel of revolution with (45/Core/ − 45) as lamination scheme.

Figure 7 .
Figure 7. First nine mode shapes for a CFCF laminated panel of translation with   0 / 30 / 60 / 90 as

Figure 8 .
Figure 8.First nine mode shapes for a CCCC thin laminated helicoidal shell with   0 / 90 / 0 / 90 as measured at the top of the shell.The two-dimensional domain is bounded by the limitations of two layers of Graphite-Epoxy characterized by the same thickness value 12 0.15m hh  .On the other hand, the lamination scheme is given by   45 / 45  .As far as the boundary conditions are concerned, the conical shell is completely clamped along the top edge, whereas the lower circumference is free.

Table 2 .
Mechanical properties of the materials employed in the numerical analyses.

Table 3 .
First 10 natural frequencies ([Hz]) for a CCCC isotropic square plate (L x = L y = 2m) of thickness h = 0.1m for several higher-order ELW approaches.The Cheb-Gau-Lob grid distribution is used for the numerical solution assuming I N = I M = 21.

Table 4 .
First 10 natural frequencies ([Hz]) for a CCCC laminated square plate (L x = L y = 2m) of thickness h = 0.1m for several higher-order ELW approaches, with and without the Murakami's function.The Cheb-Gau-Lob grid distribution is used for the numerical solution assuming I N = I M = 25.

Table 5 .
Comparison of the frequency parameter Ω = 2π f R ρ/E 2 for a laminated composite toro-circular panel (R b = 1.5m) with different boundary conditions and lamination schemes.The ELW solutions are obtained by means of the GDQ method, using the Cheb-Gau-Lob grid distribution and assuming I N = 21 and I M = 31.

Table 6 .
Comparison of the frequency parameter Ω = 2π f R ρ/E 2 for a laminated composite toro-circular panel (R b = −1.5m)with different boundary conditions and lamination schemes.The ELW solutions are obtained by means of the GDQ method, using the Cheb-Gau-Lob grid distribution and assuming I N = I M = 25.

Table 7 .
First 10 natural frequencies ([Hz]) for a CCFC laminated cylindrical surface with (90/30/45/60/0) as lamination scheme for several higher-order ELW approaches, with and without the Murakami's function.The Cheb-Gau-Lob grid distribution is used for the numerical solution assuming I N = 31 and I M = 21.

Table 8 .
First 10 natural frequencies ([Hz]) for a CCCC sandwich panel of revolution with

Table 8 .
First 10 natural frequencies ([Hz]) for a CCCC sandwich panel of revolution with (45/Core/ − 45) as lamination scheme for several higher-order ELW approaches embedded with the Murakami's function.The Cheb-Gau-Lob grid distribution is used for the numerical solution assuming I N = I M = 31.

Table 11 .
[61,83]son of the natural frequencies for a laminated composite conical shell[61,83].The solutions obtained through the present approach are compared to the results given by the corresponding ESL and Layer-Wise (LW) models.