Three-Dimensional Exact Free Vibration Analysis of Spherical, Cylindrical, and Flat One-Layered Panels

The paper proposes a three-dimensional elastic analysis of the free vibration problem of one-layered spherical, cylindrical, and flat panels. The exact solution is developed for the differential equations of equilibrium written in orthogonal curvilinear coordinates for the free vibrations of simply supported structures.These equations consider an exact geometry for shells without simplifications. The main novelty is the possibility of a general formulation for different geometries. The equations written in general orthogonal curvilinear coordinates allow the analysis of spherical shell panels and they automatically degenerate into cylindrical shell panel, cylindrical closed shell, and plate cases. Results are proposed for isotropic and orthotropic structures. An exhaustive overview is given of the vibrationmodes for a number of thickness ratios, imposed wave numbers, geometries, embeddedmaterials, and angles of orthotropy.These results can also be used as reference solutions to validate two-dimensional models for plates and shells in both analytical and numerical form (e.g., closed solutions, finite element method, differential quadrature method, and global collocation method).


Introduction
Exact solutions for the three-dimensional analysis of plates and shells have been developed by several researchers. In particular, three-dimensional static and dynamic analysis of cylinders and cylindrical/spherical shells has received considerable attention in the past few years. The development of plate and shell elements highlights the importance of research on the design of these geometries because they are fundamental for the analysis of general structures in the engineering field. For this reason, plate and shell elements need accurate validation. Three-dimensional exact solutions allow such validations and checks to be made, and they also give further details about three-dimensional effects and their importance. In the literature, works about exact three-dimensional solutions analyze separately the various geometries and they do not give a general overview of plate and shell elements. The present paper aims to fill this gap by proposing a general formulation for the equations of motion in orthogonal curvilinear coordinates that is valid for square and rectangular plates, cylindrical shell panels, spherical shell panels, and cylinders. A general overview is given for those readers interested in both plate and shell analysis. This paper exactly solves the equations of motion in general curvilinear orthogonal coordinates including an exact geometry for shell structures without simplifications. To the best of the present author's knowledge, this is the first time that this solution is proposed by means of the exponential matrix method for the three-dimensional elastic free vibration analysis of plates and shells.
The next paragraph discusses the most relevant works about three-dimensional shell analysis. Chen et al. [1] exactly studied the coupled free vibrations of a transversely isotropic cylindrical shell embedded in an elastic medium. The behavior of the cylindrical shell was analyzed based on threedimensional elasticity and three displacement functions were chosen to represent three displacement components to decouple the three-dimensional equations of motion of a transversely isotropic body. Fan and Zhang [2] presented the exact solutions for statics, dynamics, and buckling of thick open laminated cylindrical shells by means of the Cayley-Hamilton theorem. The state equation for orthotropy was established in a cylindrical coordinate system without any assumption about displacement models and stress distribution. Gasemzadeh et al. [3] studied free vibrations of cylindrical shells on the basis of three-dimensional exact theory. Simply supported boundary conditions were studied and extensive frequency parameters were obtained by solving frequency equations. The free vibrations of simply supported cross-ply cylindrical and doubly curved laminates were investigated in [4]. The three-dimensional equations of motion, expressed in terms of displacements, were reduced to a system of coupled ordinary differential equations, that were solved using the power series method. Sharma et al. [5] studied the three-dimensional free vibrations of a homogenous isotropic, viscothermoelastic hollow sphere whose surfaces were subjected to stress-free, thermally insulated or isothermal boundary conditions. The governing partial differential equations were solved into a coupled system of ordinary differential equations. Fröbenious matrix method was employed to obtain the solution. The same method was also used for the exact three-dimensional vibration analysis of a transradially isotropic, thermoelastic solid sphere [6]. Soldatos and Ye [7] proposed exact, three-dimensional, free vibration analysis of angle-ply laminated cylinders. Results were presented and discussed for relatively thick cylinders having a regular symmetric or a regular antisymmetric angleply lay-up. Armenakas et al. [8] proposed a self-contained treatment of the problem of plane harmonic waves propagation along a hollow circular cylinder in the framework of the three-dimensional theory of elasticity. The results presented frequencies of free vibrations of such cylinders for a wide range of axial and circumferential wave-lengths and of shell dimensions. An interesting comparison between frequencies obtained from a refined two-dimensional analysis, a shear deformation theory, the Flügge theory, and an exact elasticity analysis were proposed in [9]. Further details about the Flügge classical thin shell theory concerning the free vibrations of cylindrical shells with elastic boundary conditions can be found in [10]. Other comparisons between twodimensional closed form solutions and available exact 3D elastic and analytical solutions for the free vibration analysis of simply supported and clamped homogenous isotropic circular cylindrical shells were also proposed in [11]. Exact elasticity solutions were extended to functionally graded cylindrical shells in [12]. The three-dimensional linear elastodynamics equations, simplified to the case of generalized plane strain deformation in the axial direction, were solved using suitable displacement functions that identically satisfy the boundary conditions. Loy and Lam [13] proposed a layerwise approach to study the vibration of thick circular cylindrical shells on the basis of three-dimensional theory of elasticity. The governing equation was obtained using an energy minimization principle. Another interesting extension was proposed by Wang et al. [14], who investigated the free vibrations of magnetoelectroelastic cylindrical panels based on the three-dimensional theory. Further results about three-dimensional analysis of shells that are not given in closed form were shown by Efraim and Eisenberger [15] for the dynamic stiffness matrix method and by Kang and Leissa [16] and Liew et al. [17] for the three-dimensional Ritz method for vibration of spherical shells.
Three-dimensional analysis of plates is usually performed by using simpler equations that do not allow the analysis of shell geometries. Demasi [18] proposed a Navier-type method for finding the exact three-dimensional solution for isotropic thick and thin rectangular plates. The method used the mixed form of Hooke Law to write the boundary conditions on the top and bottom surfaces of the plate directly in terms of transverse stresses. Aimmanee and Batra [19] proposed an analytical solution for free vibrations of a simply supported rectangular plate made of an incompressible homogeneous linear elastic isotropic material. Some frequencies missing in previous analytical solutions for plates made of compressible materials were identified. Some of the in-plane modes of free vibration of a simply supported rectangular plate that were missed in previous exact solutions were also proposed in [20]. A three-dimensional linear elastic, small deformation theory obtained by the direct method was developed for the free vibration of simply supported, homogeneous, isotropic, thick rectangular plates in [21]. The same method was used in [22] for the flexure of simply supported homogeneous, isotropic, thick rectangular plates under arbitrary loading. This solution, in terms of infinite series, was formally exact and yielded accurate numerical results without undue effort. Useful comparisons between two-dimensional models and an exact three-dimensional solution for the free vibrations of a simply supported rectangular orthotropic thick plate are shown in [23]. On the basis of three-dimensional elasticity, Ye [24] presented a free vibration analysis of cross-ply laminated rectangular plates with clamped boundaries. The analysis was based on a recursive solution suitable for three-dimensional vibration analysis of simply supported plates. Comparisons between 2D-displacement-based models and exact results of the linear three-dimensional elasticity were also proposed in [25] for natural frequencies, displacement, and stress quantities in multilayered plates. Exact three-dimensional elasticity theory was investigated in [26] for isosceles triangular plates by means of a global three-dimensional Ritz formulation. The three-dimensional elasticity free vibration analysis of a circular plate was analyzed in [27] by means of the Ritz method with a set of orthogonal polynomial series to approximate the spatial displacements. Theoretical analysis of high frequency vibrations is indispensable in a variety of engineering designs. Zhao et al. [28] introduced a novel computational approach, the discrete singular convolution (DSC) algorithm, for high frequency vibration analysis of plate structures. A completely independent approach, the Levy method, was employed to provide exact solutions to validate the proposed method. The importance of high frequency analysis of multilayered composite plates was also confirmed in [29]. Further works about the analysis of high frequency vibrations of structures are given in [30], where the discrete singular convolution (DSC) algorithm was used and then compared with a completely independent approach (the Levy method). In [31] the first nine frequency parameters of circular and annular plates with variable thickness and combined boundary conditions were computed for different thickness to radius ratios. Three-dimensional elasticity theory was used in [31] and the Ritz method was applied to derive the eigenvalue equation. Xing and Liu [32] proposed the separation of variables to solve the Hamiltonian dual form of eigenvalue problem for transverse free vibrations of thin plates, and formulation of the natural modes in closed form was performed. The closed form natural mode exactly satisfied the governing equation of the eigenvalue problem of thin plate and was applicable for any types of boundary conditions. The extension of 3D exact analysis to functionally graded plates was performed in [33][34][35]. Exact closed form solutions (based on three-dimensional elasticity theory) were carried out in [33] for both inplane and out-of-plane free vibration of thick homogeneous simply supported rectangular plates coated by a functionally graded (FG) layer. Vel and Batra [34] presented a threedimensional exact solution for free and forced vibrations of simply supported functionally graded rectangular plates and displacement functions that identically satisfy boundary conditions. Xu and Zhou [35] used the three-dimensional elasticity theory to derive the general expressions for the displacements and stresses of the plate under static loads, which exactly satisfy the governing differential equations and the simply supported boundary conditions at four edges of the plate. The extension of three-dimensional solutions to plates embedding piezoelectric layers was given in [36][37][38][39][40][41]. The three-dimensional exact solution by Haojiang et al. [36] allowed the free vibration analysis of a circular plate made of piezoelectric material. Baillargeon and Vel [37] obtained an exact three-dimensional solution for the cylindrical bending vibration of simply supported laminated composite plates with an embedded piezoelectric shear actuator. The bending and free vibration of an imperfectly bonded orthotropic piezoelectric rectangular laminate were investigated in [38] using a three-dimensional state-space approach. Cheng et al. [39] presented an exact three-dimensional method for the distribution of mechanical and electric quantities in interior points of a symmetrically inhomogeneous and laminated piezoelectric plate in the framework of linear theory of piezoelectricity. An exact three-dimensional (3D) piezothermoelasticity solution was presented in [40] for static, free vibration and steady state harmonic response of simply supported cross-ply piezoelectric (hybrid) laminated rectangular plates with interlaminar bonding imperfections. Zhong and Shang [41] presented an exact three-dimensional analysis for a functionally graded piezoelectric rectangular plate that was simply supported and grounded along its four edges. The state equations of the functionally graded piezoelectric plate were developed based on the state space approach.
The revision of the literature about three-dimensional analysis of shells and plates demonstrates that there are a variety of interesting works concerning plate or shell geometry, and they include static and dynamic analysis, functionally graded, composite and piezoelectric materials, displacement, and mixed models. The present work gives a general formulation for all the geometries (square and rectangular plates, cylindrical and spherical shell panels, and cylindrical closed shells). The equations of motion for the dynamic case are written in general orthogonal curvilinear coordinates by using an exact geometry for shells. The system of second order differential equations is reduced to a system of first order differential equations, and afterwards it is exactly solved by using the exponential matrix method. Similar approaches have been used in [25] for the three-dimensional analysis of plates in rectilinear orthogonal coordinates and in [7] for an exact, three-dimensional, free vibration analysis of angle-ply laminated cylinders in cylindrical coordinates. The equations of motion written in orthogonal curvilinear coordinates allow general exact solutions for simply supported plate and shell geometries with constant radii of curvature. These equations are exactly solved, for the first time in this paper, by means of the exponential matrix method.

Geometrical and Constitutive Equations
A shell is a three-dimensional body bounded by two closely spaced curved surfaces where one dimension (the distance between the two surfaces) is small in comparison with the other two dimensions in the plane directions. The middle surface of the shell is the locus of points which lie midway between these surfaces. The distance between the surfaces measured along the normal to the middle surface is the thickness of the shell at that point [42,43]. Geometry and the reference system are indicated in Figure 1 where curvilinear orthogonal coordinates ( , , ) are shown. The threedimensional space (the space in which the three independent variables , , and vary) is an Euclidean space. and are the principal radii of curvature along the coordinates and , respectively. is the thickness coordinate that varies from −ℎ/2 to +ℎ/2 (where ℎ is the thickness of the structure as shown in Figure 1).
In this work, we will focus on shells with constant radii of curvature (e.g., cylindrical and spherical geometries). The geometrical relations written for shells with constant radii of curvature are a particular case of the strain-displacement equations of three-dimensional theory of elasticity in orthogonal curvilinear coordinates (see also [44][45][46]): Shock and Vibration = , where , V, and are the displacement components in , , and direction, respectively. Symbol indicates the partial derivatives; the parametric coefficients for shells, with constant radii of curvature, that use the orthogonal curvilinear coordinate system are General geometrical relations for spherical shells degenerate into geometrical relations for cylindrical shells when or is infinite (with or equals one), and they degenerate into geometrical relations for plates when both and are infinite (with = = 1) [46]. Geometrical relations (see (2)) are substituted in threedimensional linear elastic constitutive equations in orthogonal curvilinear coordinates ( , , ) for orthotropic material in the structural reference system (see [47,48]). Partial derivatives / , / , and / are indicated with subscripts , , , , and , . A closed form solution of differential equations of equilibrium for shells is obtained when coefficients 16

Equilibrium Equations
The three differential equations of equilibrium written for the case of free vibration analysis of spherical shells with constant radii of curvature and are here given (the most general form for variable radii of curvature can be found in [45,49]): where is the mass density and,V , and̈indicate the second temporal derivative of the three displacement components. Equation (3) is introduced in (4), and the closed form is obtained for simply supported shells and plates made of isotropic material or orthotropic material with 0 ∘ or 90 ∘ orthotropic angle (in both cases 16 where , , and are the displacement amplitudes in , , and directions, respectively. is the coefficient of the imaginary unit, = 2 is the circular frequency where is the frequency value, and is the time. In coefficients = / and = / , and are the half-wave numbers and and are the shell dimensions in and directions, respectively. Equilibrium equations are for spherical shell panels; they automatically degenerate into equilibrium equations for cylindrical closed/open shell panels when or is infinite (with or equals one) and into equilibrium equations for plates when and are infinite (with and equal one). In this way a unique and general closed formulation is possible for any geometry:   [34] and Srinivas et al. [22]. is the number of fictitious layers and is the order of expansion for the exponential matrix.   Simply supported aluminum cylinder with Poisson ratio ] = 0.3 and thickness ratio ℎ/ = 0.12. Lowest exact natural circular frequencies in no-dimensional form = (ℎ/ )√ / for imposed axial half-wave number = 1 and several halfwave numbers in circumferential direction. Comparison between present three-dimensional analysis and three-dimensional analysis by Armenakas et al. [8] as adapted in Bhimaraddi [9]. is the number of fictitious layers and is the order of expansion for the exponential matrix.  [8,9] 0.27491 0.27849 0.28447 0.29287

Solution for the System of Differential Equations.
The system of second order differential equations can be reduced to a system of first order differential equations by using the method described in [50,51]. This methodology allows the following system to be obtained: Equation (8) can be written in a compact form as 7 Table 3: Assessment 2 (part II). Simply supported aluminum cylinder with Poisson ratio ] = 0.3 and thickness ratio ℎ/ = 0.18. Lowest exact natural circular frequencies in no-dimensional form = (ℎ/ )√ / for imposed axial half-wave number = 1 and several halfwave numbers in circumferential direction. Comparison between present three-dimensional analysis and three-dimensional analysis by Armenakas et al. [8] as adapted in Bhimaraddi [9]. is the number of fictitious layers and is the order of expansion for the exponential matrix.  [8,9] 0.50338 0.50937 0.51934 0.53325 where U/̃= U and U = [ ]. Equation (9) can be written as In the case of plate geometry coefficients 3 , 4 , 9 , 10 , 13 , 14 , and 18 are zero because the radii of curvature and are infinite. The other coefficients 1 , 2 , 5 , 6 , 7 , 8 , 11 , 12 , 15 , 16 , 17 , and 19 are constant because parametric coefficients = = 1 and they do not depend on the thickness coordinatẽ; therefore matrices D, A, and A * are constant in the plate case. The solution of (12) can be written as [51,52] wherẽis the thickness coordinate referred to the bottom of the structure that has values from 0 at the bottom to ℎ at the top (see Figure 2). The exponential matrix is calculated with = ℎ as where I is the 6 × 6 identity matrix. This expansion has a fast convergence as indicated in [53] and it is not time consuming from the computational point of view. In the case of onelayered plate there is not a necessity to impose the continuity conditions at the interfaces for displacements and transverse shear/normal stresses; for this reason the transfer matrix T is equal to the identity matrix I: Equation (15) allows the calculation of the matrix H as Simply supported orthotropic cylinder with thickness ratio ℎ/ = 0.3 and orthotropic angle equal to 0 ∘ or 90 ∘ . First three exact natural circular frequencies in no-dimensional form = ( / )√ / 12 for imposed half-wave numbers = 0 and = 1.
Comparison between present three-dimensional analysis and three-dimensional analysis by Soldatos and Ye [7]. is the number of fictitious layers and is the order of expansion for the exponential matrix. I  II  III  I  II  III

Mode
In the case of shell geometry matrices D, A, and A * are not constant because of parametric coefficients = (1 + / ) = (1 + (̃− ℎ/2)/ ) and = (1 + / ) = (1 + (̃− ℎ/2)/ ) that depend on or̃coordinate (see Figure 2). A first method could be the use of hypothesis / = / = 0 (it is valid only for very thin shells) that means = = 1. In this case the solution is the same already seen for the plate because matrices D, A, and A * are constant. This method is not used in the benchmarks proposed in this paper because it is an approximation that is valid only for very thin shells, and it does not consider the exact geometry of the structure. However, some comparisons between the method using exact geometry and the method using = = 1 will be given at the end of the section for results. The second method (used in this paper) is the introduction of several fictitious layers where and can exactly be calculated with thẽthickness coordinate of the middle of the fictitious layers; in this way (12), (13), and (14) are valid for each layer because A * is here constant. The exponential matrix A * * is calculated in each fictitious layer with the thickness ℎ of the layer. The results proposed in this paper for the benchmarks are given with a number of fictitious layers = = 100; this value guarantees the exact convergence to the correct geometry of shells for each case investigated as demonstrated in the section for the preliminary results. In the case of layers for shell geometry, it is necessary to calculate − 1 transfer matrices T −1, by using for each interface the following conditions: that means each displacement and transverse stress component at the top ( ) of the −1 layer is equal to displacement and transverse stress component at the bottom ( ) of the layer. The calculated T −1, matrices allow linking U at the bottom of the layer with U at the top of the − 1 layer: U at the top of the layer is linked with U at the bottom of the same layer by means of the exponential matrix A * * :  Equation (19) can recursively be introduced in (20) for the − 1 interfaces to obtain by defining the matrix H for the multi-fictitious-layer shells, it is possible to write (21) in analogy with (17): It links U calculated at the top of the last fictitious layer with U calculated at the bottom of the first layer. Matrices A * * are constant because they are evaluated with , , , and calculated in the midsurface Ω 0 of the shell and with and calculated in the middle of each layer. Matrices T −1, are constant because they are evaluated with , , , and calculated in the midsurface Ω 0 of the shell and with and calculated at each fictitious interface. Equation (17) is valid for one-layered plates and (22) is valid for one-layered shells divided in fictitious layers to correctly consider the geometry. The structures are simply supported and free stresses at the top and at the bottom; this feature means: Equation (23) imposed at the top of the structure has , , , and calculated in the midsurface Ω 0 of the shell, and and are calculated at the top of the shell for̃= ℎ: Equation (23) imposed at the bottom of the structure has , , , and calculated in the midsurface Ω 0 of the shell, and and are calculated at the bottom of the shell for̃= 0: Equations (26) in matrix form are (U (ℎ ) means U calculated at the top of the shell) Equations (28)- (29) in compact form to express the free stress state at the top and bottom of the shell are in the case of plate, there are no fictitious layers: Equation (22) can be substituted in (30): in this way (33) are grouped in the following system: and introducing the (6 × 6) E matrix, (34) becomes: Equation (35) Equation (36) needs to find the roots of a higher order polynomial in = 2 . For each pair of half-wave numbers ( , ) a certain number of circular frequencies are obtained depending on the order chosen for the exponential matrices A * * and the number of fictitious layers.

Vibration Modes.
A certain number of circular frequencies are found when half-wave numbers and are imposed in the structures. For each frequency , it is possible to find the vibration mode through the thickness in terms of the three displacement components. If the frequency is substituted in (6 × 6) matrix E, this last matrix has six eigenvalues. We are interested to the null space of matrix E that means finding the (6 × 1) eigenvector related to the minimum of the six eigenvalues proposed. This null space is, for the chosen frequency , the vector U calculated at the bottom: means the transpose of the vector and the subscript means that the null space is calculated for the circular frequency .
It is possible to find U (̃) (with the three displacement components (̃), (̃), and (̃) through the thickness) for each layer by using (19)-(22); the thickness coordinatẽcan assume all the values from the bottom to the top of the structure. For the plate case the procedure is simpler because there are no fictitious layers and the number equals 1.

Results
The three-dimensional exact solution proposed has been validated by means of a comparison with three assessments given in the literature for a one-layered isotropic plate [22,34], one-layered isotropic cylinder [8,9], and one-layered orthotropic cylinder [7]. After this preliminary validation and the appropriate convergence study, the method has been used to investigate the free vibration analysis of onelayered isotropic and orthotropic square and rectangular plates, cylindrical shell panels, cylinders, and spherical shell panels (see Figure 3). The effects of the curvature terms in and parametric coefficients have also been investigated.

Preliminary Assessments.
The assessment proposed in Vel and Batra [34] compares their 3D solution with the 3D solution by Srinivas et al. [22]. The square plate is simply supported embedding one aluminium alloy layer with Young modulus = 70 GPa, Poisson ratio ] = 0.3, and mass density = 2702 kg/m 3 . The thickness value is ℎ = 1. Two different thickness ratios are considered, that is, /ℎ = 10 and /ℎ = √ 10. The first six circular frequencies in no-dimensional form = ( 2 /ℎ)√ / are given in Table 1 for imposed half-wave numbers = = 1. The 3D exact solution proposed in this paper coincides with those proposed by Vel and Batra [34] and by Srinivas et al. [22] for both thickness ratios and for each vibration mode (from the first to the sixth one). A convergence study has been proposed where the number of fictitious layers and the order of expansion for the exponential matrix have conveniently been chosen. For the plate case one fictitious layer is used because there is no effect of parametric coefficients and . For = 1, an order of expansion = 18 for the exponential matrix always gives the correct solution for each vibration mode and thickness ratio /ℎ. The calculation of the exponential matrix is not time consuming. An alternative method could be the use of fictitious layers and a lower order of expansion = 3 for the exponential matrix. In this last case, = 100 fictitious layers always give the correct solution for each   vibration mode and thickness ratio. The use of = 100 layers is not time consuming because the E matrix has always 6 × 6 dimension for any value of layers, even if the proposed method uses a layerwise approach. Solutions with = 1 and = 18 or with = 100 and = 3 are coincident for the plate case and their time solution needs a few seconds. In conclusion, the three-dimensional exact solution presented here is correct for the free vibration analysis of one-layered plates even when the structures are thick and higher order vibration modes are considered.
The assessment proposed by Khalili et al. [11] compared their 2D closed form solutions for isotropic circular cylindrical shells with the three-dimensional analysis by Armenakas et al. [8] adapted from Bhimaraddi [9]. The cylinder is isotropic and one-layered with Poisson ratio ] = 0.3. The structure is simply supported and it has radius of curvature = 1 and infinite radius of curvature in direction. The circular dimension is = 2 ; the thickness ℎ is 0.12 and 0.18; that means thickness ratios are ℎ/ = 0.12 and ℎ/ = 0.18, respectively. Several ratios cylinder length/radius of curvature ( / ) are investigated in the case of = 2, 1, 0.5, 0.25. The fundamental circular frequency in no-dimensional form = (ℎ/ )√ / is given in Tables 2 and 3 (thickness ratio ℎ/ equals 0.12 ad 0.18, resp.) for axial half-wave number = 1 and several circular halfwave numbers = 2, 4, 6, 8. The three-dimensional exact solution proposed here is coincident with the 3D analysis by Armenakas et al. [8] for each half-wave number, thickness ratio, and cylinder length proposed. For the case of shell, the use of fictitious layers is necessary to approximate the parametric coefficients and . The convergence study is shown for a value = 3 for the expansion of the exponential matrix.
= 100 fictitious layers always give the correct solution for each thickness ratio ℎ/ (see Tables 2 and  3), length of the cylinder, and half-wave numbers in circumferential direction. The present method is not time consuming; = 100 layers solution gives the correct results in a few seconds. In conclusion, the solution proposed can be considered as validated for one-layered isotropic shells.
The assessment proposed by Soldatos and Ye [7] gives three-dimensional results for multilayered orthotropic cylinders(− / + ) . A comparison with the present method is made for the cases of orthotropic angle equal to 0 ∘ or 90 ∘ (they allow the closed form solution and it also means one-layered structure). The cylinder is simply supported with a great degree of orthotropy ( 1 / 2 = 1 / 3 = 40, 12 / 2 = 13 / 2 = 0.6, 23 / 2 = 0.5, ] 12 = ] 13 = ] 23 = 0.25, and homogeneous mass density ). The length is 10 and the radius of curvature equals 10. The circular dimension is = 2 . The cylinder is very thick (ℎ = 3); that means thickness ratio is ℎ/ = 0.3. In Table 4 the first three circular frequencies in no-dimensional form = ( / )√ / 12 are given for half-wave numbers = 0 and = 1, and orthotropic angle equals 0 ∘ or 90 ∘ . The 3D solution by Soldatos and Ye [7] is correctly obtained for each mode and angle by using = 100 fictitious layers and = 3 order of expansion. ( = 100, = 3) solution is not heavy from the computational point of view because the E matrix has always 6 × 6 dimension even if the method uses a layerwise approach. In conclusion, the solution is also validated for orthotropic single-layered thick shells.
All the new benchmarks proposed in Section 4.2 will use ( = 1, = 18) solution for plate cases and ( = 100, = 3) solution for shell cases. Both solutions are very convenient from the computational point of view (they need only a few seconds) and they allow the correct analysis for each thick and thin isotropic/orthotropic plate and shell.

4.2.
Benchmarks. Ten different benchmarks are proposed to show a complete overview of the free vibration analysis of one-layered plates and shells. Five different geometries are taken into account in Figure 3. All these structures are simply supported. The square plate has dimensions = = 1 and thickness ratios /ℎ = 100, 50, 10, 5. The rectangular plate has dimensions = 1 and = 3 and the thickness ratios are /ℎ = 100, 50, 10, 5. The cylindrical shell panel has a radius of curvature = 10 and an infinite radius of curvature in direction. The dimensions are = ( /3) and = 20, and the thickness ratios are /ℎ = 1000, 100, 10, 5. The cylinder has the same radii of curvature of the cylindrical shell panel, but it is closed in circumferential direction that means = 2 . The other dimension is = 100 and the thickness ratios are /ℎ = 1000, 100, 10, 5. The last geometry is the spherical shell panel with radii of curvature = = 10, dimensions = = ( /3) , and thickness ratios /ℎ = 1000, 100, 10, 5. Each geometry includes two materials: one isotropic layer in aluminium alloy Al2024 with Young modulus = 73 GPa, Poisson ratio ] = 0.3, and mass density = 2800 kg/m 3 and one orthotropic layer in Gr/Ep with Young modulus components 1 = 132.38 GPa and 2 = 3 = 10.756 GPa, shear modulus components  Table 5), one-layered orthotropic square plate (see Table 6), one-layered isotropic rectangular plate (see Table 7), one-layered orthotropic rectangular plate (see Table 8), one-layered isotropic cylindrical shell panel (see Table 9), one-layered orthotropic cylindrical shell panel (see Table 10), one-layered isotropic cylinder (see Table 11), one-layered orthotropic cylinder (see Table 12), one-layered isotropic spherical shell panel (see Table 13), and one-layered orthotropic spherical shell panel (see Table 14). The first three circular frequencies in no-dimensional form = ( 2 /ℎ)√ / 2 are calculated in Tables 5-14 for various pairs of half-wave numbers and several thickness ratios. The vibration modes plotted in Figures 4-12 are given in terms of no-dimensional values such as * = /| max |, V * = V/|V max |, * = /| max |, and * =̃/ℎ. Table 5 presents a square isotropic plate (benchmark 1). It shows the first three vibration modes for each possible combination of half-wave numbers ( , ) from (0, 1) until to (2,2). Both thick and thin plates are analyzed. The plate is square and isotropic; therefore results for half-wave numbers (0, 1) are equal to those for (1, 0), and the same considerations are made for the half-wave numbers (1, 2) and (2, 1). The first three vibration modes through the thickness of a thick plate ( /ℎ = 10) in terms of the three displacement components * , V * , and * are given in Figure 4 for half-wave numbers (1, 0), (1, 1), (2,1), and (2, 2). The second vibration mode is an in-plane mode for each pair of half-wave numbers because transverse displacement * through the thickness is zero. The first mode has a constant transverse displacement * through the thickness and linear in-plane displacements * and V * . The transverse displacement is not perfectly constant for higher values of half-wave number (effect of higher order modes). The in-plane displacement V * is zero for the first case because the half-wave number = 0. The third vibration mode always has linear transverse displacement * through the thickness for each pair of half-wave numbers. In-plane displacements * and V * are constant (they are not perfectly constant for higher values of half-wave numbers; it is a typical effect of higher order modes). The in-plane displacement V * is zero for the first case because of zero half-wave number . Table 6 shows a square orthotropic plate (benchmark 2) and lists the first three vibration modes for each possible combination of half-wave numbers ( , ) from (0, 1) until to (2, 2) and for two different fibre orientation angles (0 ∘ and 90 ∘ ). For the angle = 0 ∘ results for half-wave numbers (0, 1) are not equal to those for (1, 0), and the same considerations are made for the half-wave numbers (1, 2) and (2, 1). This is due to the in-plane orthotropy. The same considerations are valid for the angle = 90 ∘ . It is important to notice that frequencies for (0, 1) and = 0 ∘ are equal to those for (1, 0) and = 90 ∘ and frequencies for (1, 0) and = 0 ∘ are equal to those for (0, 1) and = 90 ∘ . The same considerations remain valid for (1, 2) and (2, 1). These properties are valid for square plates and they are not true for rectangular plates. Figures 5 and 6 show the first three vibration modes through the thickness in terms of the three displacement components * , V * , and * for half-wave numbers (1, 0), (1, 1), (2, 1), and (2, 2) and for a thick plate ( /ℎ = 10) with fiber orientation = 0 ∘ or = 90 ∘ . The considerations already made for the first benchmark in Figure 4 are also valid in this case. The orthotropy also generates a coupling between in-plane and out-of-plane behavior that gives displacements that are not constant or linear though the thickness even when low halfwave numbers are considered. Moreover, thickness modes are more complicated for higher half-wave numbers. For these reasons a 2D model will have more difficulties to obtain the 3D behavior of the benchmark 2 with respect to the 3D behavior of the isotropic benchmark 1. The change of orthotropic angle has an important effect on the frequency values and a negligible effect on the vibration modes through the thickness. Tables 7 and 8 show the same cases already seen in Tables 5 and 6, but a rectangular plate is analyzed in place of a square plate. This geometry gives a difference in Table 7 (benchmark 3) between modes with (0, 1) and those with (1, 0) and between modes with (1, 2) and those with (2, 1); in fact, the two in-plane dimensions and are not the same ( = 3 ). This difference due to the in-plane dimensions is also important in Table 8 (benchmark 4) where frequencies for (0, 1) and = 0 ∘ are not equal to the frequencies for (1, 0) and = 90 ∘ and and frequencies for (1, 0) and = 0 ∘ are not equal to the frequencies for (0, 1) and = 90 ∘ . The same considerations can be made for higher half-wave numbers (1, 2) and (2, 1). Table 9 gives the first three vibration modes of the isotropic cylindrical shell panel (benchmark 5) for several combinations of half-wave numbers ( , ). The shell is isotropic but the frequencies are always different for each couple of ( , ) because of the geometry with different and dimensions. The same considerations are also valid for benchmark 6 in Table 10 (the orthotropic cylindrical shell  panel). The geometry with different and dimensions does not give any similarity when angle changes in combination with the half-wave numbers. The orthotropic effects give a coupling between in-plane and out-of-plane behavior that will require higher order 2D models to obtain the 3D analysis in terms of frequencies and vibration modes. Tables 11 and 12 show results for isotropic and orthotropic one-layered cylinders, respectively. Results do not add further considerations with respect to those already given for the open cylindrical shell. The main difference is in terms of vibration modes through the thickness as shown in Figure 7 for isotropic cylinder (benchmark 7) and in Figures 8 and 9 for orthotropic cylinder (benchmark 8 with orthotropic angle equals 0 ∘ or 90 ∘ ). In Figure 7 the first mode for each pair of ( , ) always has constant displacements * and * , and a quasi-constant V * displacement, but in the case = 0 and = 1 the only component different from zero is * . The second mode always has constant displacements V * and * and a linear * displacement, but in the case = 0 and = 1 * is zero and * and V * are constant. The third mode has quasi-constant * , V * , and * displacements, but in the case = 0 and = 1 the component * is zero. In Figure 8 the first three modes are given for the orthotropic case with = 0 ∘ . The considerations about the first vibration mode for each pair of ( , ) are the same as those already given for the isotropic benchmark. The second and third modes are different with respect to the isotropic case because of the coupling due to the in-plane orthotropy (see in particular * displacement for the second mode and V * displacement for the third mode). The change of angle of orthotropy gives large differences in terms of frequency values and small differences in terms of vibration modes (see the comparison between Figure 8 for = 0 ∘ and Figure 9 for = 90 ∘ ).
Tables 13 and 14 investigate isotropic spherical shell panel (benchmark 9) and orthotropic spherical shell panel (benchmark 10), respectively. The shell geometry considered here has the same radii of curvature in directions and and the same dimensions and . For this reason the considerations made for the square plate (Tables 5 and 6) are also valid in these two last cases: results for half-wave numbers (0, 1) are equal to those for (1, 0), and the same similarity is valid for the half-wave numbers (1, 2) and (2, 1) (isotropic shell). For the angle = 0 ∘ , results for half-wave numbers (0, 1) are not equal to those for (1, 0), and the same considerations are made for the half-wave numbers (1, 2) and (2, 1). This is due to the in-plane orthotropy. Similar considerations are also valid for the angle = 90 ∘ . It is important to notice that frequencies for (0, 1) and = 0 ∘ are equal to the frequencies for (1, 0) and = 90 ∘ and frequencies for (1, 0) and = 0 ∘ are equal to those for (0, 1) and = 90 ∘ . The same feature is confirmed for (1, 2) and (2, 1). Figure 10 shows the vibration modes for isotropic case, while Figures  11 and 12 show the vibration modes for = 0 ∘ and = 90 ∘ orthotropic cases, respectively. The curvature effect is very important and gives rather complicated modes through the thickness (in particular for the orthotropic spherical shells in Figures 11 and 12). These demonstrate that the Table 14: Benchmark 10. Simply supported orthotropic spherical shell panel. First three vibration modes in terms of no-dimensional circular frequency = ( 2 /ℎ)√ / 2 for several combinations of half-wave numbers and . = 100 fictitious layers and = 3 order of expansion for the exponential matrix.  use of 3D analysis or refined 2D models is fundamental to calculate these frequencies and modes. The change of angle of orthotropy has a large influence on the frequency values, but this influence is smaller for vibration modes in Figures 11  and 12.
The parametric coefficient effects, in particular the effect in cylindrical shell panels and both and effects in spherical shell panels, are investigated in Tables 15 and  16. Table 15 shows the comparison between the present 3D solution using the exact geometry and the same solution using = 1 (which means / ≃ 0) for the benchmark 5 concerning the isotropic cylindrical shell panel. The approximation = 1 is valid for very thin shells ( /ℎ equals 1000 or 100) even if it could have some problems for particular vibration modes and higher values of half-wave numbers ( , ). For thick shells the use of the exact geometry without the approximation = 1 is mandatory. Table 16 shows the comparison between the present 3D solution using the exact geometry and the same solution using the approximation = = 1 (which means / = / ≃ 0) for the benchmark 10 concerning the orthotropic spherical shell panel. The approximation = = 1 is valid for thin shells ( /ℎ equals 1000 or 100), for low values of half-wave numbers, and for the lower vibration modes. The approximation gives significative errors for thick shells where the hypothesis / = / ≃ 0 is not valid.

Conclusions
The differential equations of equilibrium in orthogonal curvilinear coordinates for the free vibrations of simply supported structures have been exactly solved in three-dimensional form by using the exponential matrix method. The proposed general 3D formulation uses an exact geometry for shells, and it allows results for spherical, open cylindrical, closed cylindrical and flat panels to be obtained. The first three vibration modes have been investigated for several geometries, both isotropic and orthotropic layers, various thickness ratios, and imposed half-wave numbers. The modes plotted through the thickness make it possible to recognize the most complicated cases and these results will be useful benchmarks to validate future refined 2D models. The method is simple and intuitive and it will be extended to multilayered structures (also embedding functionally graded layers). This extension will give a global three-dimensional overview of the free vibration problem of one-layered and multilayered plates and shells.     Third mode for m = 2 and n = 2 u * * w * u * * w * u * * w * First mode for m = 2 and n = 1 First mode for m = 2 and n = 2 First mode for m = 2 and n = 3 Second mode for m = 0 and n = 1 Second mode for m = 2 and n = 1 Second mode for m = 2 and n = 2 Second mode for m = 2 and n = 3 Third mode for m = 0 and n = 1 Third mode for m = 2 and n = 1 Third mode for m = 2 and n = 2 Third mode for m = 2 and n = 3 u * * w * u * * w * u * * w * Figure 9: Benchmark 8, simply supported orthotropic cylinder ( = 90 ∘ ). First three vibration modes in terms of displacement components through the thickness for thickness ratio /ℎ = 10 and several combinations of half-wave numbers. Third mode for m = 2 and n = 1 First mode for m = 2 and n = 2 Second mode for m = 2 and n = 2 Third mode for m = 2 and n = 2 u * * w * u * * w * u * * w * Figure 10: Benchmark 9, simply supported isotropic spherical shell panel. First three vibration modes in terms of displacement components through the thickness for thickness ratio /ℎ = 5 and several combinations of half-wave numbers.

Shock and Vibration
The present method has a fast convergence and it is not heavy from the computational point of view for each geometry (plates and shells). This feature is an important advantage with respect to other methods proposed in the literature that are valid only for a chosen geometry.